CFPP$ NOCONCUR R SUBROUTINE MONINP(KM,ntrac,A,B,TAU,RTG, & U1,V1,T1,Q1, & PSTAR,RBSOIL,CD,CH,FM,FH,TSEA,QSS,DPHI,SPD1,KPBL, & SI,DEL,SL,SLK,RCL,DELTIM,THOUR, & DUSFC,DVSFC,DTSFC,DQSFC,HPBL,HGAMT,HGAMQ &, HEAT, EVAP,alpha) cc USE MACHINE , ONLY : kind_phys implicit none cc include 'constant.h' cc integer iprt,is,iun,k,kk,km,kmpbl,kpbl,lond,ntrac,ilat,ilon real(kind=kind_phys) aphi16,aphi5,bet1,betaq,betat,betaw,bvf2 real(kind=kind_phys) cd,cfac,ch,conq,cont,conw,conwrc,deltim real(kind=kind_phys) dk,dkmax,dkmin,dphi,dq1,dqsfc,dsdz2,dsdzq real(kind=kind_phys) dsdzt,dsig,dt,dthe1,dtodsd,dtodsu,dtsfc real(kind=kind_phys) dusfc,dvsfc,dw2,dw2min,evap,fh,fm,g,gamcrq real(kind=kind_phys) gamcrt,gocp,gor,heat,hgamq,hgamt,hol,hpbl real(kind=kind_phys) pfac,phih,phim,prmax,prmin,prnum,pstar real(kind=kind_phys) qmin,qss,qtend,rbcr,rbdn,rbint,rbsoil,rbup real(kind=kind_phys) rcl,rdt,rdz,rdzt1,ri,rimin,rl2,rlam,rocp real(kind=kind_phys) rone,rzero,sfcfrac,sflux,shr2,spd1,spdk2 real(kind=kind_phys) sri,the1,the1v,thekv,thermal,thesv,thour real(kind=kind_phys) ti,tsea,ttend,tvd,tvu,ustar,utend,vk,vk2 real(kind=kind_phys) vpert,vtend,wscale,xkzo,zfac,zfmin,zk,zl1 real(kind=kind_phys) tem,alpha cc csela %INCLUDE DBMONINP; csela PARAMETER(CP=#CP,G=#G,RD=#RD,RV=#RV,HVAP=#HVAP) PARAMETER(g=grav) PARAMETER(GOR=G/RD,GOCP=G/CP,ROCP=RD/CP) PARAMETER(CONT=1000.*CP/G,CONQ=1000.*HVAP/G,CONW=1000./G) PARAMETER(RLAM=150.,VK=0.4,VK2=VK*VK,PRMIN=1.0,PRMAX=4.) PARAMETER(DW2MIN=0.0001,DKMIN=1.0,DKMAX=1000.,RIMIN=-100.) ccc PARAMETER(RBCR=0.5,CFAC=7.8,PFAC=2.0,SFCFRAC=0.1) PARAMETER(RBCR=0.25,CFAC=7.8,PFAC=2.0,SFCFRAC=0.1) PARAMETER(QMIN=1.E-8,XKZO=1.0,ZFMIN=1.E-8,APHI5=5.,APHI16=16.) c PARAMETER(GAMCRT=3.,GAMCRQ=2.E-3) PARAMETER(GAMCRT=3.,GAMCRQ=0.) PARAMETER(IUN=84) real(kind=kind_phys) A(KM),B(KM),TAU(KM),RTG(KM*ntrac), & U1(KM),V1(KM),T1(KM),Q1(KM*ntrac), & SI(KM+1),DEL(KM),CL(KM),SL(KM),SLK(KM) real(kind=kind_phys) DZOT(KM-1),RDZT(KM-1), 1 ZI(KM+1),ZL(KM), 2 DKU(KM-1),DKT(KM-1), 3 AL(KM-1),AD(KM),AU(KM-1), 4 A1(KM),A2(KM*ntrac) logical pblflg, sfcflg, stable C C----------------------------------------------------------------------- C 601 FORMAT(1X,' MONINP LAT LON STEP HOUR ',3I6,F6.1) 602 FORMAT(1X,' K',' Z',' T',' TH', 1 ' TVH',' Q',' U',' V', 2 ' SP') 603 FORMAT(1X,I5,8F9.1) 604 FORMAT(1X,' SFC',9X,F9.1,18X,F9.1) 605 FORMAT(1X,' K ZL SPD2 THEKV THE1V' 1 ,' THERMAL RBUP') 606 FORMAT(1X,I5,6F8.2) 607 FORMAT(1X,' KPBL HPBL FM FH HGAMT', 1 ' HGAMQ WS USTAR CD CH') 608 FORMAT(1X,I5,9F8.2) 609 FORMAT(1X,' K PR DKT DKU ',I5,3F8.2) 610 FORMAT(1X,' K PR DKT DKU ',I5,3F8.2,' L2 RI T2', 1 ' SR2 ',2F8.2,2E10.2) C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C COMPUTE PRELIMINARY VARIABLES C IPRT = 0 IF(IPRT.EQ.1) THEN CCC LATD = 0 LOND = 0 ELSE CCC LATD = 0 LOND = 0 ENDIF C DT = 2. * DELTIM RDT = 1. / DT RDZT1 = GOR * SL(1) / DEL(1) KMPBL = KM / 2 C DO K = 1,KM-1 RDZT(K) = GOR * SI(K+1) / (SL(K) - SL(K+1)) DZOT(K) = LOG(SI(K+1) / SI(K)) / GOR ENDDO C ZI(1) = 0. DO K = 1, KM - 1 ZI(K+1) = ZI(K) - T1(K) * DZOT(K) ENDDO C CCC IF(LAT.EQ.LATD) THEN CCC I = LOND CCC WRITE(IUN,602) CCC DO K = KM-1,1,-1 CCC ZLKP1 = (ZI(K) + ZI(K+1))/2. CCC THETA = T1(K)/SLK(K) CCC THETAV = T1(K)/SLK(K)*(1.+FV*MAX(Q1(K),QMIN)) CCC SPEED = SQRT(MAX(RCL*(U1(k)**2+V1(k)**2),1.)) CCC WRITE(IUN,603) K,ZLKP1,T1(k)-273.15,THETA,THETAV, CCC 1 Q1(k)*1000.,U1(k),V1(k),SPEED CCC ENDDO CCC WRITE(IUN,604) TSEA,QSS*1000. CCC ENDIF C DUSFC = 0. DVSFC = 0. DTSFC = 0. DQSFC = 0. HGAMT = 0. HGAMQ = 0. WSCALE = 0. KPBL = 1 HPBL = ZI(2) PBLFLG = .TRUE. SFCFLG = .TRUE. IF(RBSOIL.GT.0.0) SFCFLG = .FALSE. C BET1 = DT*RDZT1*SPD1/T1(1) BETAW = BET1*CD BETAT = BET1*CH BETAQ = DPHI*BETAT C ZL1 = 0.-(T1(1)+TSEA)/2.*LOG(SL(1))/GOR USTAR = SQRT(CD*SPD1**2) C THESV = TSEA*(1.+FV*MAX(QSS,QMIN)) THE1 = T1(1)/SLK(1) THE1V = THE1*(1.+FV*MAX(Q1(1),QMIN)) THERMAL = THE1V DTHE1 = (THE1-TSEA) DQ1 = (MAX(Q1(1),QMIN) - MAX(QSS,QMIN)) HEAT = -CH*SPD1*DTHE1 EVAP = -CH*SPD1*DQ1 C C C COMPUTE THE FIRST GUESS OF PBL HEIGHT C STABLE = .FALSE. ZL(1) = ZL1 RBUP = RBSOIL DO K = 2, KMPBL IF(.NOT.STABLE) THEN RBDN = RBUP ZL(k) = ZL(K-1) - (T1(k)+T1(K-1))/2 * & LOG(SL(K)/SL(K-1)) / GOR THEKV = T1(k)/SLK(K)*(1.+FV*MAX(Q1(k),QMIN)) SPDK2 = MAX(RCL*(U1(k)**2+V1(k)**2),1.) RBUP = (THEKV-THE1V)*(G*ZL(k)/THE1V)/SPDK2 KPBL = K STABLE = RBUP.GT.RBCR ENDIF ENDDO C K = KPBL IF(RBDN.GE.RBCR) THEN RBINT = 0. ELSEIF(RBUP.LE.RBCR) THEN RBINT = 1. ELSE RBINT = (RBCR-RBDN)/(RBUP-RBDN) ENDIF HPBL = ZL(K-1) + RBINT*(ZL(k)-ZL(K-1)) IF(HPBL.LT.ZI(KPBL)) KPBL = KPBL - 1 C HOL = MAX(RBSOIL*FM*FM/FH,RIMIN) IF(SFCFLG) THEN HOL = MIN(HOL,-ZFMIN) ELSE HOL = MAX(HOL,ZFMIN) ENDIF C HOL = HOL*HPBL/ZL1*SFCFRAC IF(SFCFLG) THEN ! PHIM = (1.-APHI16*HOL)**(-1./4.) ! PHIH = (1.-APHI16*HOL)**(-1./2.) TEM = 1.0 / (1. - APHI16*HOL) PHIH = SQRT(TEM) PHIM = SQRT(PHIH) ELSE PHIM = (1.+APHI5*HOL) PHIH = PHIM ENDIF WSCALE = USTAR/PHIM WSCALE = MIN(WSCALE,USTAR*APHI16) WSCALE = MAX(WSCALE,USTAR/APHI5) C C COMPUTE THE SURFACE VARIABLES FOR PBL HEIGHT ESTIMATION C UNDER UNSTABLE CONDITIONS C SFLUX = HEAT + EVAP*FV*THE1 IF(SFCFLG.AND.SFLUX.GT.0.0) THEN HGAMT = MIN(CFAC*HEAT/WSCALE,GAMCRT) HGAMQ = MIN(CFAC*EVAP/WSCALE,GAMCRQ) VPERT = HGAMT + FV*THE1*HGAMQ VPERT = MIN(VPERT,GAMCRT) THERMAL = THERMAL + MAX(VPERT,0.) HGAMT = MAX(HGAMT,0.0) HGAMQ = MAX(HGAMQ,0.0) ELSE PBLFLG = .FALSE. ENDIF C IF(PBLFLG) THEN KPBL = 1 HPBL = ZI(2) ENDIF C C ENHANCE THE PBL HEIGHT BY CONSIDERING THE THERMAL C IF(PBLFLG) THEN STABLE = .FALSE. RBUP = RBSOIL ENDIF DO K = 2, KMPBL IF(.NOT.STABLE.AND.PBLFLG) THEN RBDN = RBUP ZL(k) = ZL(K-1) - (T1(k)+T1(K-1))/2 * & LOG(SL(K)/SL(K-1)) / GOR THEKV = T1(k)/SLK(K)*(1.+FV*MAX(Q1(k),QMIN)) SPDK2 = MAX(RCL*(U1(k)**2+V1(k)**2),1.) RBUP = (THEKV-THERMAL)*(G*ZL(k)/THE1V)/SPDK2 KPBL = K STABLE = RBUP.GT.RBCR ENDIF ENDDO C IF(PBLFLG) THEN K = KPBL IF(RBDN.GE.RBCR) THEN RBINT = 0. ELSEIF(RBUP.LE.RBCR) THEN RBINT = 1. ELSE RBINT = (RBCR-RBDN)/(RBUP-RBDN) ENDIF HPBL = ZL(K-1) + RBINT*(ZL(k)-ZL(K-1)) IF(HPBL.LT.ZI(KPBL)) KPBL = KPBL - 1 IF(KPBL.LE.1) PBLFLG = .FALSE. ENDIF C C COMPUTE DIFFUSION COEFFICIENTS BELOW PBL C ! if (hpbl .eq. zl1) then ! print *,' hpbl=',hpbl,' kpbl=',kpbl ! *,' zl1=',zl1,' ilat=',ilat,' ilon=',ilon ! print *,' u1=',u1 ! print *,' v1=',v1 ! print *,' t1=',t1 ! print *,' q1=',q1 ! print *,' tsea=',tsea,' sl1=',sl(1),' gor=',gor ! endif DO K = 1, KMPBL IF(KPBL.GT.K) THEN PRNUM = (PHIH/PHIM+CFAC*VK*.1) PRNUM = MIN(PRNUM,PRMAX) PRNUM = MAX(PRNUM,PRMIN) ZFAC = MAX((1.-(ZI(K+1)-ZL1)/ 1 (HPBL-ZL1)), ZFMIN) DKU(k) = XKZO + WSCALE*VK*ZI(K+1) 1 *ALPHA*ZFAC**PFAC DKT(k) = DKU(k)/PRNUM DKU(k) = MIN(DKU(k),DKMAX) DKU(k) = MAX(DKU(k),DKMIN) DKT(k) = MIN(DKT(k),DKMAX) DKT(k) = MAX(DKT(k),DKMIN) ENDIF ENDDO C C COMPUTE DIFFUSION COEFFICIENTS OVER PBL (FREE ATMOSPHERE) C DO K = 1, KM-1 IF(K.GE.KPBL) THEN TI =0.5*(T1(k)+T1(K+1)) RDZ =RDZT(K)/TI DW2 =RCL*((U1(k)-U1(K+1))**2+(V1(k)-V1(K+1))**2) SHR2 =MAX(DW2,DW2MIN)*RDZ**2 TVD =T1(k)*(1.+FV*MAX(Q1(k),QMIN)) TVU =T1(K+1)*(1.+FV*MAX(Q1(K+1),QMIN)) BVF2 =G*(GOCP+RDZ*(TVU-TVD))/TI RI =MAX(BVF2/SHR2,RIMIN) ZK =VK*ZI(K+1) RL2 =(ZK*RLAM/(RLAM+ZK))**2 DK =RL2*SQRT(SHR2) IF(RI.LT.0.) THEN ! UNSTABLE REGIME SRI = SQRT(-RI) DKU(k) = XKZO + DK*(1+8.*(-RI)/(1+1.746*SRI)) DKT(k) = XKZO + DK*(1+8.*(-RI)/(1+1.286*SRI)) ELSE ! STABLE REGIME DKT(k) = XKZO + DK/(1+5.*RI)**2 PRNUM = 1.0 + 2.1*RI PRNUM = MIN(PRNUM,PRMAX) DKU(k) = (DKT(k)-XKZO)*PRNUM + XKZO ENDIF C DKU(k) = MIN(DKU(k),DKMAX) DKU(k) = MAX(DKU(k),DKMIN) DKT(k) = MIN(DKT(k),DKMAX) DKT(k) = MAX(DKT(k),DKMIN) C CCC IF(I.EQ.LOND.AND.LAT.EQ.LATD) THEN CCC PRNUM = DKU(k)/DKT(k) CCC WRITE(IUN,610) K,PRNUM,DKT(k),DKU(k),RL2,RI, CCC 1 BVF2,SHR2 CCC ENDIF C ENDIF ENDDO C C COMPUTE TRIDIAGONAL MATRIX ELEMENTS FOR HEAT AND MOISTURE C AD(1) = 1. A1(1) = T1(1)-BETAT*(T1(1)/SLK(1)-TSEA) A2(1) = Q1(1)-BETAQ*(MAX(Q1(1),QMIN)-MAX(QSS,QMIN)) if(ntrac.ge.2) then do k = 2, ntrac is = (k-1) * km A2(1+is) = q1(1+is) enddo endif C DO K = 1,KM-1 DTODSD = DT/DEL(K) DTODSU = DT/DEL(K+1) DSIG = SL(K)-SL(K+1) RDZ = RDZT(K)*2./(T1(k)+T1(K+1)) IF(PBLFLG.AND.K.LT.KPBL) THEN DSDZT = DSIG*DKT(k)*RDZ*(GOCP-HGAMT/HPBL) DSDZQ = DSIG*DKT(k)*RDZ*(-HGAMQ/HPBL) A2(k) = A2(k)+DTODSD*DSDZQ A2(k+1) = Q1(k+1)-DTODSU*DSDZQ ELSE DSDZT = DSIG*DKT(k)*RDZ*(GOCP) A2(k+1) = Q1(k+1) ENDIF DSDZ2 = DSIG*DKT(k)*RDZ*RDZ AU(k) = -DTODSD*DSDZ2 AL(k) = -DTODSU*DSDZ2 AD(k) = AD(k)-AU(k) AD(k+1) = 1.-AL(k) A1(k) = A1(k)+DTODSD*DSDZT A1(k+1) = T1(k+1)-DTODSU*DSDZT ENDDO if(ntrac.ge.2) then do kk = 2, ntrac is = (kk-1) * km do k = 1, km - 1 A2(k+1+is) = q1(k+1+is) enddo enddo endif C C SOLVE TRIDIAGONAL PROBLEM FOR HEAT AND MOISTURE C CALL TRIDIN(KM,ntrac,AL,AD,AU,A1,A2,AU,A1,A2) C C RECOVER TENDENCIES OF HEAT AND MOISTURE C DO K = 1,KM TTEND = (A1(k)-T1(k))*RDT QTEND = (A2(k)-Q1(k))*RDT TAU(k) = TAU(k)+TTEND RTG(k) = RTG(k)+QTEND DTSFC = DTSFC+CONT*DEL(K)*PSTAR*TTEND DQSFC = DQSFC+CONQ*DEL(K)*PSTAR*QTEND ENDDO if(ntrac.ge.2) then do kk = 2, ntrac is = (kk-1) * km do k = 1, km QTEND = (A2(K+is)-Q1(K+is))*RDT RTG(K+is) = RTG(K+is)+QTEND enddo enddo endif C C COMPUTE TRIDIAGONAL MATRIX ELEMENTS FOR MOMENTUM C AD(1) = 1.+BETAW A1(1) = U1(1) A2(1) = V1(1) C DO K = 1,KM-1 DTODSD = DT/DEL(K) DTODSU = DT/DEL(K+1) DSIG = SL(K)-SL(K+1) RDZ = RDZT(K)*2./(T1(k)+T1(k+1)) DSDZ2 = DSIG*DKU(k)*RDZ*RDZ AU(k) = -DTODSD*DSDZ2 AL(k) = -DTODSU*DSDZ2 AD(k) = AD(k)-AU(k) AD(k+1)= 1.-AL(k) A1(k+1)= U1(k+1) A2(k+1)= V1(k+1) ENDDO C C SOLVE TRIDIAGONAL PROBLEM FOR MOMENTUM C CALL TRIDI2(KM,AL,AD,AU,A1,A2,AU,A1,A2) C C RECOVER TENDENCIES OF MOMENTUM C CONWRC = CONW*SQRT(RCL) DO K = 1,KM UTEND = (A1(k)-U1(k))*RDT VTEND = (A2(k)-V1(k))*RDT B(k)= B(k)+UTEND A(k)= A(k)+VTEND DUSFC = DUSFC+CONWRC*DEL(K)*PSTAR*UTEND DVSFC = DVSFC+CONWRC*DEL(K)*PSTAR*VTEND ENDDO C RETURN END CFPP$ NOCONCUR R C----------------------------------------------------------------------- SUBROUTINE TRIDI2(N,CL,CM,CU,R1,R2,AU,A1,A2) csela %INCLUDE DBTRIDI2; cc USE MACHINE , ONLY : kind_phys implicit none integer k,n real(kind=kind_phys) fk cc real(kind=kind_phys) CL(2:N),CM(N),CU(N-1),R1(N),R2(N), & AU(N-1),A1(N),A2(N) C----------------------------------------------------------------------- FK=1./CM(1) AU(1)=FK*CU(1) A1(1)=FK*R1(1) A2(1)=FK*R2(1) DO K=2,N-1 FK=1./(CM(k)-CL(k)*AU(k-1)) AU(k)=FK*CU(k) A1(k)=FK*(R1(k)-CL(k)*A1(k-1)) A2(k)=FK*(R2(k)-CL(k)*A2(k-1)) ENDDO FK=1./(CM(N)-CL(N)*AU(N-1)) A1(N)=FK*(R1(N)-CL(N)*A1(N-1)) A2(N)=FK*(R2(N)-CL(N)*A2(N-1)) DO K=N-1,1,-1 A1(k)=A1(k)-AU(k)*A1(k+1) A2(k)=A2(k)-AU(k)*A2(k+1) ENDDO C----------------------------------------------------------------------- RETURN END CFPP$ NOCONCUR R C----------------------------------------------------------------------- SUBROUTINE TRIDIN(N,nt,CL,CM,CU,R1,R2,AU,A1,A2) csela %INCLUDE DBTRIDI2; cc USE MACHINE , ONLY : kind_phys implicit none integer is,k,kk,n,nt real(kind=kind_phys) fk cc real(kind=kind_phys) CL(2:N),CM(N),CU(N-1),R1(N),R2(N*nt), & AU(N-1),A1(N),A2(N*nt) C----------------------------------------------------------------------- FK=1./CM(1) AU(1)=FK*CU(1) A1(1)=FK*R1(1) do k = 1, nt is = (k-1) * n FK=1./CM(1) a2(1+is) = fk * r2(1+is) enddo DO K=2,N-1 FK=1./(CM(k)-CL(k)*AU(k-1)) AU(k)=FK*CU(k) A1(k)=FK*(R1(k)-CL(k)*A1(k-1)) ENDDO do kk = 1, nt is = (kk-1) * n DO K=2,N-1 FK=1./(CM(k)-CL(k)*AU(k-1)) A2(k+is)=FK*(R2(k+is)-CL(k)*A2(k+is-1)) ENDDO ENDDO FK=1./(CM(N)-CL(N)*AU(N-1)) A1(N)=FK*(R1(N)-CL(N)*A1(N-1)) do k = 1, nt is = (k-1) * n FK=1./(CM(n)-CL(n)*AU(n-1)) A2(n+is)=FK*(R2(n+is)-CL(n)*A2(n+is-1)) enddo DO K=N-1,1,-1 A1(k)=A1(k)-AU(k)*A1(k+1) ENDDO do kk = 1, nt is = (kk-1) * n DO K=n-1,1,-1 A2(k+is)=A2(k+is)-AU(k)*A2(k+is+1) ENDDO ENDDO C----------------------------------------------------------------------- RETURN END