SUBROUTINE DELDIFS(RTE,WE,QME,XE,YE,TEME, X RTO,WO,QMO,XO,YO,TEMO,DELTIM,SL, X LS_NODE,COEF00,K_LEVEL,hybrid,gen_coord_hybrid) ! use resol_def use layout1 use coordinate_def ! hmhj use deldifs_def use physcons, rerth => con_rerth, rd => con_rd, cp => con_cp IMPLICIT NONE ! logical hybrid, gen_coord_hybrid REAL(KIND=KIND_EVOD) RTE(LEN_TRIE_LS,2,LEVS,NTRAC) &, WE(LEN_TRIE_LS,2) &, QME(LEN_TRIE_LS,2) &, XE(LEN_TRIE_LS,2) &, YE(LEN_TRIE_LS,2) &, TEME(LEN_TRIE_LS,2) &, PE(LEN_TRIE_LS,2) ! REAL(KIND=KIND_EVOD) RTO(LEN_TRIO_LS,2,LEVS,NTRAC) &, WO(LEN_TRIO_LS,2) &, QMO(LEN_TRIO_LS,2) &, XO(LEN_TRIO_LS,2) &, YO(LEN_TRIO_LS,2) &, TEMO(LEN_TRIO_LS,2) &, PO(LEN_TRIO_LS,2) ! REAL(KIND=KIND_EVOD) DELTIM, SL(LEVS) ! INTEGER LS_NODE(LS_DIM,3) ! !CMR LS_NODE(1,1) ... LS_NODE(LS_MAX_NODE,1) : VALUES OF L !CMR LS_NODE(1,2) ... LS_NODE(LS_MAX_NODE,2) : VALUES OF JBASEV !CMR LS_NODE(1,3) ... LS_NODE(LS_MAX_NODE,3) : VALUES OF JBASOD ! REAL(KIND=KIND_EVOD) COEF00(LEVS,NTRAC) ! INTEGER K_LEVEL ! INTEGER I,IS,IT,JDEL,JDELH,K,KD,KU INTEGER L,LOCL,N,N0,ND,NP,NPD ! INTEGER INDEV INTEGER INDOD integer indev1,indev2 integer indod1,indod2 real(kind=kind_evod), parameter :: rkappa = cp / rd REAL(KIND=KIND_EVOD) DN1,REALVAL,RTNP,DF_DK,FACT &, SLRD0,FTRD1,RFACT,RFACTRD,RTRD1,FSHK, tem ! REAL(KIND=KIND_EVOD), parameter :: CONS0=0.0, CONS1=1.0, CONS2=2.0 ! INTEGER INDLSEV,JBASEV INTEGER INDLSOD,JBASOD ! INCLUDE 'function2' ! ! print *,' enter deldifs_fd ' ! hmhj !...................................................................... ! IF(K_LEVEL.EQ.0) THEN ! CALL COUNTPERF(0,15,0.) !! allocate(RTRD(LEVS),RTHK(LEVS),sf(levs)) ALLOCATE ( DNE(LEN_TRIE_LS) ) ALLOCATE ( DNO(LEN_TRIO_LS) ) ALLOCATE ( BKLY(levs) ) ! hmhj ALLOCATE ( CKLY(levs) ) ! hmhj BKLY(:) = 1.0 CKLY(:) = 0.0 if (gen_coord_hybrid) then ! hmhj DO k=1,LEVS ! hmhj ! hmhj ak5, bk5, ck5 in gen_coord_hybrid is the same order as model index BKLY(k)=0.5*(bk5(k)+bk5(k+1)) ! hmhj CKLY(k)=0.5*(ck5(k)+ck5(k+1))*rkappa/thref(k) ! hmhj if( me.eq.0 ) ! hmhj & print*,'sl bkly ckly in deldif=',k,sl(k),bkly(k),ckly(k)! hmhj enddo ! hmhj else if (hybrid) then ! hmhj DO k=1,LEVS ! hmhj sl(k) go bottom to top but bk(k) go top to bottom BKLY(k)=0.5*(bk5(levs-k+1)+bk5(levs-k+2))/SL(k) if( me.eq.0 ) print*,'sl bkly in deldif=',k,sl(k),bkly(k) enddo endif ! N0=0 ! MAXIMUM WAVENUMBER FOR ZERO DIFFUSION JDEL=8 ! ORDER OF DIFFUSION (EVEN POWER TO RAISE DEL) NP=JCAP FSHK=1.0 ! EXTRA HEIGHT-DEPENDENT DIFFUSION FACTOR PER SCALE HEIGHT IF(JCAP.GT.170) THEN ! RECIPROCAL OF TIME SCALE OF DIFFUSION AT REFERENCE WAVENUMBER NP RTNP=(JCAP/170.)**4*1.1/3600 FSHK=2.2 ! EXTRA HEIGHT-DEPENDENT DIFFUSION FACTOR PER SCALE HEIGHT ELSEIF(JCAP.EQ.170) THEN ! RECIPROCAL OF TIME SCALE OF DIFFUSION AT REFERENCE WAVENUMBER NP RTNP=4*3.E15/(RERTH**4)*FLOAT(80*81)**2 ELSEIF(JCAP.EQ.126) THEN ! hmhj ! BELOW HAS BEEN TESTED IN SIGMA-THETA FOR 2 YEAR CFS RUN ! hmhj ! FSHK=1.5 ! EXTRA HEIGHT-DEPENDENT DIFFUSION FACTOR PER SCALE HEIGHT FSHK=2.2 ! EXTRA HEIGHT-DEPENDENT DIFFUSION FACTOR PER SCALE HEIGHT RTNP=4*3.E15/(RERTH**4)*FLOAT(80*81)**2 ! hmhj ELSE ! RECIPROCAL OF TIME SCALE OF DIFFUSION AT REFERENCE WAVENUMBER NP RTNP=1*3.E15/(RERTH**4)*FLOAT(80*81)**2 ENDIF ! IF (ME.EQ.0) THEN PRINT 6,RTNP,NP,N0,JDEL 6 FORMAT(' HORIZONTAL DIFFUSION PARAMETERS'/ & ' EFFECTIVE ',6PF10.3,' MICROHERTZ AT WAVENUMBER ',I4/ & ' MAXIMUM WAVENUMBER FOR ZERO DIFFUSION ',I4/ & ' ORDER OF DIFFUSION ',I2) ENDIF ! SLRD0=0.002 ! SIGMA LEVEL AT WHICH TO BEGIN RAYLEIGH DAMPING RTRD1=1./(5*86400) ! RECIPROCAL OF TIME SCALE PER SCALE HEIGHT ! ABOVE BEGINNING SIGMA LEVEL FOR RAYLEIGH DAMPING ! DO K=1,LEVS IF(SL(K).LT.SLRD0) THEN if (k .gt. levr) then RTRD(K)=RTRD1*LOG(SLRD0/SL(K)) ** 2 else RTRD(K)=RTRD1*LOG(SLRD0/SL(K)) endif ELSE RTRD(K)=0 ENDIF RTHK(K)=(SL(K))**LOG(1/FSHK) ENDDO ! JDELH=JDEL/2 NPD=MAX(NP-N0,0) REALVAL=NPD*(NPD+1) DN1=CONS2*RTNP/REALVAL**JDELH ! !...................................................................... ! DO LOCL=1,LS_MAX_NODE L=LS_NODE(LOCL,1) JBASEV=LS_NODE(LOCL,2) INDEV=INDLSEV(L,L) DO N=L,JCAP,2 ND=MAX(N-N0,0) REALVAL=ND*(ND+1) DNE(INDEV)=DN1*REALVAL**JDELH INDEV=INDEV+1 ENDDO ENDDO ! !...................................................................... ! DO LOCL=1,LS_MAX_NODE L=LS_NODE(LOCL,1) JBASEV=LS_NODE(LOCL,2) if (mod(L,2).eq.mod(jcap+1,2)) then DNE(INDLSEV(JCAP+1,L))=CONS0 ! SET THE EVEN (N-L) TERMS OF THE TOP ROW TO ZERO ENDIF ENDDO ! !...................................................................... ! DO LOCL=1,LS_MAX_NODE L=LS_NODE(LOCL,1) JBASOD=LS_NODE(LOCL,3) INDOD=INDLSOD(L+1,L) DO N=L+1,JCAP,2 ND=MAX(N-N0,0) REALVAL=ND*(ND+1) DNO(INDOD)=DN1*REALVAL**JDELH INDOD=INDOD+1 ENDDO ENDDO ! !...................................................................... ! DO LOCL=1,LS_MAX_NODE L=LS_NODE(LOCL,1) JBASOD=LS_NODE(LOCL,3) if (mod(L,2).ne.mod(jcap+1,2)) then DNO(INDLSOD(JCAP+1,L))=CONS0 ! SET THE ODD (N-L) TERMS OF THE TOP ROW TO ZERO ENDIF ENDDO ! !...................................................................... ! DO K=1,LEVS KD=MAX(K-1,1) KU=MIN(K+1,LEVS) SF(K)=SL(K)/(SL(KU)-SL(KD))/SQRT(CONS2) !CONSTANT ENDDO ! CALL COUNTPERF(1,15,0.) !! RETURN ENDIF ! !...................................................................... ! CALL COUNTPERF(0,13,0.) !! K=K_LEVEL !! ! ! TEM = COEF00(K,1) * BKLY(K) ! DO LOCL=1,LS_MAX_NODE L=LS_NODE(LOCL,1) JBASEV=LS_NODE(LOCL,2) IF (L.EQ.0) THEN N0=2 ELSE N0=L ENDIF indev1 = indlsev(N0,L) if (mod(L,2).eq.mod(jcap+1,2)) then indev2 = indlsev(jcap+1,L) else indev2 = indlsev(jcap ,L) endif !! DO N = N0, JCAP+1, 2 DO INDEV = indev1 , indev2 FACT = DELTIM*DNE(INDEV)*RTHK(K) RFACT = CONS1/(CONS1+FACT) RFACTRD = CONS1/(CONS1+FACT+DELTIM*RTRD(K)) WE(INDEV,1) = WE(INDEV,1)*RFACTRD WE(INDEV,2) = WE(INDEV,2)*RFACTRD XE(INDEV,1) = XE(INDEV,1)*RFACTRD XE(INDEV,2) = XE(INDEV,2)*RFACTRD RTE(INDEV,1,1,1) = RTE(INDEV,1,1,1)*RFACT RTE(INDEV,2,1,1) = RTE(INDEV,2,1,1)*RFACT PE(INDEV,1)=BKLY(K)*QME(INDEV,1)+CKLY(K)*TEME(INDEV,1) ! hmhj PE(INDEV,2)=BKLY(K)*QME(INDEV,2)+CKLY(K)*TEME(INDEV,2) ! hmhj YE(INDEV,1) = ( YE(INDEV,1)+ X FACT*COEF00(K,1)* PE(INDEV,1) )*RFACT ! hmhj YE(INDEV,2) = ( YE(INDEV,2) + X FACT*COEF00(K,1)* PE(INDEV,2) )*RFACT ! hmhj ENDDO ENDDO ! !...................................................................... ! ! DO L = 0, JCAP DO LOCL=1,LS_MAX_NODE L=LS_NODE(LOCL,1) JBASOD=LS_NODE(LOCL,3) indod1 = indlsod(L+1,L) if (mod(L,2).eq.mod(jcap+1,2)) then indod2 = indlsod(jcap ,L) else indod2 = indlsod(jcap+1,L) endif ! DO N = L+1, JCAP+1, 2 DO INDOD = indod1 , indod2 FACT = DELTIM*DNO(INDOD)*RTHK(K) RFACT = CONS1/(CONS1+FACT) RFACTRD = CONS1/(CONS1+FACT+DELTIM*RTRD(K)) WO(INDOD,1) = WO(INDOD,1)*RFACTRD WO(INDOD,2) = WO(INDOD,2)*RFACTRD XO(INDOD,1) = XO(INDOD,1)*RFACTRD XO(INDOD,2) = XO(INDOD,2)*RFACTRD RTO(INDOD,1,1,1) = RTO(INDOD,1,1,1)*RFACT RTO(INDOD,2,1,1) = RTO(INDOD,2,1,1)*RFACT PO(INDOD,1)=BKLY(K)*QMO(INDOD,1)+CKLY(K)*TEMO(INDOD,1) ! hmhj PO(INDOD,2)=BKLY(K)*QMO(INDOD,2)+CKLY(K)*TEMO(INDOD,2) ! hmhj YO(INDOD,1) = ( YO(INDOD,1)+ X FACT*COEF00(K,1)* PO(INDOD,1) )*RFACT ! hmhj YO(INDOD,2) = ( YO(INDOD,2)+ X FACT*COEF00(K,1)* PO(INDOD,2) )*RFACT ! hmhj ENDDO ENDDO ! !...................................................................... ! if (ntrac .gt. 1) then do it=2,ntrac ! ! DO L = 0, JCAP DO LOCL=1,LS_MAX_NODE L=LS_NODE(LOCL,1) JBASEV=LS_NODE(LOCL,2) IF (L.EQ.0) THEN N0=2 ELSE N0=L ENDIF indev1 = indlsev(N0,L) if (mod(L,2).eq.mod(jcap+1,2)) then indev2 = indlsev(jcap+1,L) else indev2 = indlsev(jcap ,L) endif ! DO N = N0, JCAP+1, 2 DO INDEV = indev1 , indev2 FACT = DELTIM*DNE(INDEV)*RTHK(K) RFACT = CONS1/(CONS1+FACT) ! RTE(INDEV,1,1,IT) = RTE(INDEV,1,1,IT) * RFACT ! RTE(INDEV,2,1,IT) = RTE(INDEV,2,1,IT) * RFACT PE(INDEV,1)=BKLY(K)*QME(INDEV,1)+CKLY(K)*TEME(INDEV,1) ! hmhj PE(INDEV,2)=BKLY(K)*QME(INDEV,2)+CKLY(K)*TEME(INDEV,2) ! hmhj RTE(INDEV,1,1,IT) = ( RTE(INDEV,1,1,IT)+ X FACT*COEF00(K,it)* PE(INDEV,1) )*RFACT ! hmhj RTE(INDEV,2,1,IT) = ( RTE(INDEV,2,1,IT)+ X FACT*COEF00(K,it)* PE(INDEV,2) )*RFACT ! hmhj ENDDO ENDDO ! !...................................................................... ! ! DO L = 0, JCAP DO LOCL=1,LS_MAX_NODE L=LS_NODE(LOCL,1) JBASOD=LS_NODE(LOCL,3) indod1 = indlsod(L+1,L) if (mod(L,2).eq.mod(jcap+1,2)) then indod2 = indlsod(jcap ,L) else indod2 = indlsod(jcap+1,L) endif ! DO N = L+1, JCAP+1, 2 DO INDOD = indod1 , indod2 FACT = DELTIM*DNO(INDOD)*RTHK(K) RFACT = CONS1/(CONS1+FACT) ! RTO(INDOD,1,1,IT) = RTO(INDOD,1,1,IT) * RFACT ! RTO(INDOD,2,1,IT) = RTO(INDOD,2,1,IT) * RFACT ! PO(INDOD,1)=BKLY(K)*QMO(INDOD,1)+CKLY(K)*TEMO(INDOD,1) ! hmhj PO(INDOD,2)=BKLY(K)*QMO(INDOD,2)+CKLY(K)*TEMO(INDOD,2) ! hmhj RTO(INDOD,1,1,IT) = ( RTO(INDOD,1,1,IT)+ X FACT*COEF00(K,it)* PO(INDOD,1) )*RFACT ! hmhj RTO(INDOD,2,1,IT) = ( RTO(INDOD,2,1,IT)+ X FACT*COEF00(K,it)* PO(INDOD,2) )*RFACT ! hmhj ENDDO ENDDO enddo ! Tracer do loop end endif ! CALL COUNTPERF(1,13,0.) ! print *,' leave deldifs_fd ' ! hmhj !! RETURN END