SUBROUTINE DELDIFS_tracers(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 gfs_dyn_MACHINE, only : kind_evod use gfs_dyn_resol_def, only : jcap,levr,levs,ntrac use namelist_dynamics_def, only : hdif_fac, hdif_fac2, slrd0 use gfs_dyn_layout1, only : len_trie_ls,len_trio_ls, & ls_dim,ls_max_node,me use gfs_dyn_coordinate_def, only : bk5,ck5,thref ! hmhj use deldifs_def , only : bkly,ckly,dne,dno,rthk,rtrd,sf use gfs_dyn_physcons, only : 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 &, 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' ! !...................................................................... ! 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 == 0 ) ! hmhj & print*,'gen_cor 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 == 0 ) ! hmhj & print*,'hybrid sl bkly ckly in deldif=',k,sl(k),bkly(k),ckly(k)! hmhj print*,'k_level=',k_level 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 FSHK = 1.0*hdif_fac ! EXTRA HEIGHT-DEPENDENT DIFFUSION FACTOR PER SCALE HEIGHT ! IF(JCAP > 170) THEN ! RECIPROCAL OF TIME SCALE OF DIFFUSION AT REFERENCE WAVENUMBER NP RTNP = hdif_fac2*(JCAP/170.)**4*1.1/3600 ! RTNP =(JCAP/170.)**4*1.1/3600 FSHK = 2.2*hdif_fac ! EXTRA HEIGHT-DEPENDENT DIFFUSION FACTOR PER SCALE HEIGHT ELSEIF(JCAP == 170) THEN ! RECIPROCAL OF TIME SCALE OF DIFFUSION AT REFERENCE WAVENUMBER NP RTNP = hdif_fac2*4*3.E15/(RERTH**4)*FLOAT(80*81)**2 ! RTNP = 4*3.E15/(RERTH**4)*FLOAT(80*81)**2 ELSEIF(JCAP == 126) THEN ! hmhj ! BELOW HAS BEEN TESTED IN SIGMA-THETA FOR 2 YEAR CFS RUN ! hmhj RTNP = hdif_fac2*4*3.E15/(RERTH**4)*FLOAT(80*81)**2 ! hmhj ! RTNP = 4*3.E15/(RERTH**4)*FLOAT(80*81)**2 ! hmhj FSHK = 1.5*hdif_fac ! EXTRA HEIGHT-DEPENDENT DIFFUSION FACTOR PER SCALE HEIGHT ELSE ! RECIPROCAL OF TIME SCALE OF DIFFUSION AT REFERENCE WAVENUMBER NP RTNP = hdif_fac2*1*3.E15/(RERTH**4)*FLOAT(80*81)**2 ! RTNP = 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/ & ' in deldifs ORDER OF DIFFUSION ',I2) ENDIF ! ! SLRD0=0.002 ! SIGMA LEVEL AT WHICH TO BEGIN RAYLEIGH DAMPING ! SLRD0=0.0005 ! 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) < SLRD0) THEN if (k > 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 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 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 ! ! print *,' leave deldifs_tracers ' ! hmhj !! RETURN END