CFPP$ NOCONCUR R SUBROUTINE MONINP(IX,IM,KM,ntrac,DV,DU,TAU,RTG, & U1,V1,T1,Q1, & PSK,RBSOIL,FM,FH,TSEA,QSS,HEAT,EVAP,STRESS,SPD1,KPBL, ! & PSK,RBSOIL,CD,CH,FM,FH,TSEA,QSS,DPHI,SPD1,KPBL, & PRSI,DEL,PRSL,PRSLK,PHII,PHIL,RCL,DELTIM, & DUSFC,DVSFC,DTSFC,DQSFC,HPBL,HGAMT,HGAMQ) ! USE MACHINE , ONLY : kind_phys USE PHYSCONS, grav => con_g, RD => con_RD, CP => con_CP &, HVAP => con_HVAP, ROG => con_ROG, FV => con_FVirt implicit none ! ! include 'constant.h' ! ! ! Arguments ! integer IX, IM, KM, ntrac, KPBL(IM) ! real(kind=kind_phys) DELTIM real(kind=kind_phys) DV(IM,KM), DU(IM,KM), & TAU(IM,KM), RTG(IM,KM,ntrac), & U1(IX,KM), V1(IX,KM), & T1(IX,KM), Q1(IX,KM,ntrac), & PSK(IM), RBSOIL(IM), ! & CD(IM), CH(IM), & FM(IM), FH(IM), & TSEA(IM), QSS(IM), & SPD1(IM), ! & DPHI(IM), SPD1(IM), & PRSI(IX,KM+1), DEL(IX,KM), & PRSL(IX,KM), PRSLK(IX,KM), & PHII(IX,KM+1), PHIL(IX,KM), & RCL(IM), DUSFC(IM), & dvsfc(IM), dtsfc(IM), & DQSFC(IM), HPBL(IM), & HGAMT(IM), hgamq(IM) ! ! Locals ! integer i,iprt,is,iun,k,kk,kmpbl,lond ! real(kind=kind_phys) betaq(IM), betat(IM), betaw(IM), real(kind=kind_phys) evap(IM), heat(IM), phih(IM), & phim(IM), rbdn(IM), rbup(IM), & the1(IM), stress(im), beta(im), & the1v(IM), thekv(IM), thermal(IM), & thesv(IM), ustar(IM), wscale(IM) ! & thesv(IM), ustar(IM), wscale(IM), zl1(IM) ! real(kind=kind_phys) RDZT(IM,KM-1), & ZI(IM,KM+1), ZL(IM,KM), & DKU(IM,KM-1), DKT(IM,KM-1), & AL(IM,KM-1), AD(IM,KM), & AU(IM,KM-1), A1(IM,KM), & A2(IM,KM*ntrac), THETA(IM,KM) logical pblflg(IM), sfcflg(IM), stable(IM) ! real(kind=kind_phys) aphi16, aphi5, bet1, bvf2, & cfac, conq, cont, conw, & conwrc, dk, dkmax, dkmin, & dq1, dsdz2, dsdzq, dsdzt, & dsig, dt, dthe1, dtodsd, & dtodsu, dw2, dw2min, g, & gamcrq, gamcrt, gocp, gor, gravi, & hol, pfac, prmax, prmin, prinv, & prnum, qmin, qtend, rbcr, & rbint, rdt, rdz, ! & rbint, rdt, rdz, rdzt1, & ri, rimin, rl2, rlam, rlamun, & rone, rzero, sfcfrac, & sflux, shr2, spdk2, sri, & tem, ti, ttend, tvd, & tvu, utend, vk, vk2, & vpert, vtend, xkzo(im,km), zfac, & zfmin, zk, tem1, xkzm cc parameter (gravi=1.0/grav) PARAMETER(g=grav) PARAMETER(GOR=G/RD,GOCP=G/CP) PARAMETER(CONT=1000.*CP/G,CONQ=1000.*HVAP/G,CONW=1000./G) ! PARAMETER(RLAM=150.0,VK=0.4,VK2=VK*VK,PRMIN=1.0,PRMAX=4.) PARAMETER(RLAM=30.0,VK=0.4,VK2=VK*VK,PRMIN=1.0,PRMAX=4.) ! PARAMETER(RLAM=50.0,VK=0.4,VK2=VK*VK,PRMIN=1.0,PRMAX=4.) PARAMETER(DW2MIN=0.0001,DKMIN=0.0,DKMAX=1000.,RIMIN=-100.) PARAMETER(RBCR=0.25,CFAC=7.8,PFAC=2.0,SFCFRAC=0.1) ! PARAMETER(RBCR=0.5,CFAC=7.8,PFAC=2.0,SFCFRAC=0.1) ! PARAMETER(QMIN=1.E-8,XKZM=3.0,ZFMIN=1.E-8,APHI5=5.,APHI16=16.) ! PARAMETER(QMIN=1.E-8,XKZM=2.0,ZFMIN=1.E-8,APHI5=5.,APHI16=16.) PARAMETER(QMIN=1.E-8,XKZM=1.0,ZFMIN=1.E-8,APHI5=5.,APHI16=16.) ! PARAMETER(QMIN=1.E-8,XKZM=0.5,ZFMIN=1.E-8,APHI5=5.,APHI16=16.) ! PARAMETER(QMIN=1.E-8,XKZM=0.25,ZFMIN=1.E-8,APHI5=5.,APHI16=16.) ! PARAMETER(QMIN=1.E-8,XKZM=0.10,ZFMIN=1.E-8,APHI5=5.,APHI16=16.) ! PARAMETER(QMIN=1.E-8,XKZM=0.0,ZFMIN=1.E-8,APHI5=5.,APHI16=16.) ! PARAMETER(GAMCRT=3.,GAMCRQ=2.E-3) PARAMETER(GAMCRT=3.,GAMCRQ=0., RLAMUN=150.0) ! PARAMETER(GAMCRT=3.,GAMCRQ=0., RLAMUN=30.0) PARAMETER(IUN=84) ! 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 if (IX .lt. im) stop ! ! 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 KMPBL = KM / 2 ! do k=1,km do i=1,im zi(i,k) = phii(i,k) * gravi zl(i,k) = phil(i,k) * gravi enddo enddo ! do k=1,kmpbl do i=1,im theta(i,k) = t1(i,k) * psk(i) / prslk(i,k) enddo enddo C DO K = 1,KM-1 DO I=1,IM RDZT(I,K) = 1.0 / (ZL(I,K+1) - ZL(I,K)) ! RDZT(I,K) = GOR * PRSI(I,K+1) / (PRSL(I,K) - PRSL(I,K+1)) ! if (prsi(i,1) .gt. 60.0) then ! tem1 = max((prsi(i,k+1)-30.0)/(prsi(i,1)-30.0), 0.0) ! else ! tem1 = 0.0 ! endif ! xkzo(i,k) = xkzm * tem1 * tem1 !! tem1 = (zi(i,k+1) - zi(i,1)) * 0.002 !! xkzo(i,k) = xkzm * exp(-tem1) ENDDO ENDDO C DO I = 1,IM DUSFC(I) = 0. DVSFC(I) = 0. DTSFC(I) = 0. DQSFC(I) = 0. HGAMT(I) = 0. HGAMQ(I) = 0. WSCALE(I) = 0. KPBL(I) = 1 HPBL(I) = ZI(I,2) PBLFLG(I) = .TRUE. SFCFLG(I) = .TRUE. IF(RBSOIL(I).GT.0.0) SFCFLG(I) = .FALSE. ENDDO !! DO I=1,IM ! RDZT1 = GOR * prSL(i,1) / DEL(i,1) ! BET1 = DT*RDZT1*SPD1(I)/T1(I,1) ! BETA(I) = DT*RDZT1/T1(I,1) BETA(I) = DT / (zi(i,2)-zi(i,1)) ! BETAW(I) = BET1*CD(I) ! BETAT(I) = BET1*CH(I) ! BETAQ(I) = DPHI(I)*BETAT(I) ENDDO C DO I=1,IM ! ZL1(i) = 0.-(T1(I,1)+TSEA(I))/2.*LOG(PRSL(I,1)/PRSI(I,1))*ROG ! USTAR(I) = SQRT(CD(I)*SPD1(I)**2) USTAR(I) = SQRT(STRESS(I)) ENDDO C DO I=1,IM THESV(I) = TSEA(I)*(1.+FV*MAX(QSS(I),QMIN)) THE1(I) = THETA(I,1) THE1V(I) = THE1(I)*(1.+FV*MAX(Q1(I,1,1),QMIN)) THERMAL(I) = THE1V(I) ! DTHE1 = (THE1(I)-TSEA(I)) ! DQ1 = (MAX(Q1(I,1,1),QMIN) - MAX(QSS(I),QMIN)) ! HEAT(I) = -CH(I)*SPD1(I)*DTHE1 ! EVAP(I) = -CH(I)*SPD1(I)*DQ1 ENDDO C C C COMPUTE THE FIRST GUESS OF PBL HEIGHT C DO I=1,IM STABLE(I) = .FALSE. ! ZL(i,1) = ZL1(i) RBUP(I) = RBSOIL(I) ENDDO DO K = 2, KMPBL DO I = 1, IM IF(.NOT.STABLE(I)) THEN RBDN(I) = RBUP(I) ! ZL(I,k) = ZL(I,K-1) - (T1(i,k)+T1(i,K-1))/2 * ! & LOG(PRSL(I,K)/PRSL(I,K-1)) * ROG THEKV(I) = THETA(i,k)*(1.+FV*MAX(Q1(i,k,1),QMIN)) SPDK2 = MAX(RCL(i)*(U1(i,k)**2+V1(i,k)**2),1.) RBUP(I) = (THEKV(I)-THE1V(I))*(G*ZL(I,k)/THE1V(I))/SPDK2 KPBL(I) = K STABLE(I) = RBUP(I).GT.RBCR ENDIF ENDDO ENDDO C DO I = 1,IM K = KPBL(I) IF(RBDN(I).GE.RBCR) THEN RBINT = 0. ELSEIF(RBUP(I).LE.RBCR) THEN RBINT = 1. ELSE RBINT = (RBCR-RBDN(I))/(RBUP(I)-RBDN(I)) ENDIF HPBL(I) = ZL(I,K-1) + RBINT*(ZL(I,K)-ZL(I,K-1)) IF(HPBL(I).LT.ZI(I,KPBL(I))) KPBL(I) = KPBL(I) - 1 ENDDO !! DO I=1,IM HOL = MAX(RBSOIL(I)*FM(I)*FM(I)/FH(I),RIMIN) IF(SFCFLG(I)) THEN HOL = MIN(HOL,-ZFMIN) ELSE HOL = MAX(HOL,ZFMIN) ENDIF C ! HOL = HOL*HPBL(I)/ZL1(I)*SFCFRAC HOL = HOL*HPBL(I)/ZL(I,1)*SFCFRAC IF(SFCFLG(I)) THEN ! PHIM = (1.-APHI16*HOL)**(-1./4.) ! PHIH = (1.-APHI16*HOL)**(-1./2.) TEM = 1.0 / (1. - APHI16*HOL) PHIH(I) = SQRT(TEM) PHIM(I) = SQRT(PHIH(I)) ELSE PHIM(I) = (1.+APHI5*HOL) PHIH(I) = PHIM(I) ENDIF WSCALE(I) = USTAR(I)/PHIM(I) WSCALE(I) = MIN(WSCALE(I),USTAR(I)*APHI16) WSCALE(I) = MAX(WSCALE(I),USTAR(I)/APHI5) ENDDO C C COMPUTE THE SURFACE VARIABLES FOR PBL HEIGHT ESTIMATION C UNDER UNSTABLE CONDITIONS C DO I = 1,IM SFLUX = HEAT(I) + EVAP(I)*FV*THE1(I) IF(SFCFLG(I).AND.SFLUX.GT.0.0) THEN HGAMT(I) = MIN(CFAC*HEAT(I)/WSCALE(I),GAMCRT) HGAMQ(I) = MIN(CFAC*EVAP(I)/WSCALE(I),GAMCRQ) VPERT = HGAMT(I) + FV*THE1(I)*HGAMQ(I) VPERT = MIN(VPERT,GAMCRT) THERMAL(I) = THERMAL(I) + MAX(VPERT,0.) HGAMT(I) = MAX(HGAMT(I),0.0) HGAMQ(I) = MAX(HGAMQ(I),0.0) ELSE PBLFLG(I) = .FALSE. ENDIF ENDDO C DO I = 1,IM IF(PBLFLG(I)) THEN KPBL(I) = 1 HPBL(I) = ZI(I,2) ENDIF ENDDO C C ENHANCE THE PBL HEIGHT BY CONSIDERING THE THERMAL C DO I = 1, IM IF(PBLFLG(I)) THEN STABLE(I) = .FALSE. RBUP(I) = RBSOIL(I) ENDIF ENDDO DO K = 2, KMPBL DO I = 1, IM IF(.NOT.STABLE(I).AND.PBLFLG(I)) THEN RBDN(I) = RBUP(I) ! ZL(I,k) = ZL(I,K-1) - (T1(i,k)+T1(i,K-1))/2 * ! & LOG(PRSL(I,K)/PRSL(I,K-1)) * ROG THEKV(I) = THETA(i,k)*(1.+FV*MAX(Q1(i,k,1),QMIN)) SPDK2 = MAX(RCL(i)*(U1(i,k)**2+V1(i,k)**2),1.) RBUP(I) = (THEKV(I)-THERMAL(I))*(G*ZL(I,k)/THE1V(I))/SPDK2 KPBL(I) = K STABLE(I) = RBUP(I).GT.RBCR ENDIF ENDDO ENDDO C DO I = 1,IM IF(PBLFLG(I)) THEN K = KPBL(I) IF(RBDN(I).GE.RBCR) THEN RBINT = 0. ELSEIF(RBUP(I).LE.RBCR) THEN RBINT = 1. ELSE RBINT = (RBCR-RBDN(I))/(RBUP(I)-RBDN(I)) ENDIF HPBL(I) = ZL(I,K-1) + RBINT*(ZL(I,k)-ZL(I,K-1)) IF(HPBL(I).LT.ZI(I,KPBL(I))) KPBL(I) = KPBL(I) - 1 IF(KPBL(I).LE.1) PBLFLG(I) = .FALSE. ENDIF ENDDO !! ! DO K = 1,KM-1 ! DO I=1,IM ! if (rbsoil(i) .gt. 0.0 .or. (.not. pblflg(i))) then ! tem1 = max((prsi(i,k+1)-30.0)/(prsi(i,1)-30.0), 0.0) ! xkzo(i,k) = xkzm * tem1 * tem1 !! tem1 = (zi(i,k+1) - zi(i,1)) * 0.002 !! xkzo(i,k) = xkzm * exp(-tem1) ! else ! xkzo(i,k) = 0.0 ! ! if (pblflg(i)) then ! if (sfcflg(i)) then ! xkzo(i,k) = 0.0 ! else ! tem1 = (zi(i,k+1) - zi(i,1)) * 0.0005 ! tem1 = (100.0 - prsi(i,k+1)) * 0.075 ! tem1 = 100.0 - prsi(i,k+1) ! tem1 = max(0.0, 100.0 - prsi(i,k+1)) ! tem1 = tem1 * tem1 * 0.00075 ! tem1 = tem1 * tem1 * 0.001 ! tem1 = tem1 * tem1 * 0.0011 ! tem1 = tem1 * tem1 * 0.0012 ! ! tem1 = 1.0 - prsi(i,k+1) / prsi(i,1) ! tem1 = tem1 * tem1 * 5.0 ! tem1 = tem1 * tem1 * 7.5 ! tem1 = tem1 * tem1 * 10.0 ! tem1 = tem1 * tem1 * 12.0 ! ! xkzo(i,k) = xkzm * min(1.0, exp(-tem1)) ! if (xkzo(i,k) .lt. 0.01) xkzo(i,k) = 0.0 ! endif ! ENDDO ! ENDDO ! DO K = 1,KM-1 DO I=1,IM tem1 = 1.0 - prsi(i,k+1) / prsi(i,1) tem1 = tem1 * tem1 * 10.0 xkzo(i,k) = xkzm * min(1.0, exp(-tem1)) ENDDO ENDDO !! C C COMPUTE DIFFUSION COEFFICIENTS BELOW PBL C DO K = 1, KMPBL DO I=1,IM IF(KPBL(I).GT.K) THEN PRINV = 1.0 / (PHIH(I)/PHIM(I)+CFAC*VK*.1) PRINV = MIN(PRINV,PRMAX) PRINV = MAX(PRINV,PRMIN) ! ZFAC = MAX((1.-(ZI(I,K+1)-ZL1(I))/ ! 1 (HPBL(I)-ZL1(I))), ZFMIN) ZFAC = MAX((1.-(ZI(I,K+1)-ZL(I,1))/ 1 (HPBL(I)-ZL(I,1))), ZFMIN) DKU(i,k) = XKZO(i,k) + WSCALE(I)*VK*ZI(I,K+1) 1 * ZFAC**PFAC DKT(i,k) = DKU(i,k)*PRINV DKU(i,k) = MIN(DKU(i,k),DKMAX) DKU(i,k) = MAX(DKU(i,k),DKMIN) DKT(i,k) = MIN(DKT(i,k),DKMAX) DKT(i,k) = MAX(DKT(i,k),DKMIN) ENDIF ENDDO ENDDO C C COMPUTE DIFFUSION COEFFICIENTS OVER PBL (FREE ATMOSPHERE) C DO K = 1, KM-1 DO I=1,IM IF(K.GE.KPBL(I)) THEN ! TI = 0.5*(T1(i,k)+T1(i,K+1)) TI = 2.0 / (T1(i,k)+T1(i,K+1)) ! RDZ = RDZT(I,K)/TI ! RDZ = RDZT(I,K) * TI RDZ = RDZT(I,K) DW2 = RCL(i)*((U1(i,k)-U1(i,K+1))**2 & + (V1(i,k)-V1(i,K+1))**2) SHR2 = MAX(DW2,DW2MIN)*RDZ*RDZ TVD = T1(i,k)*(1.+FV*MAX(Q1(i,k,1),QMIN)) TVU = T1(i,K+1)*(1.+FV*MAX(Q1(i,K+1,1),QMIN)) ! BVF2 = G*(GOCP+RDZ*(TVU-TVD))/TI BVF2 = G*(GOCP+RDZ*(TVU-TVD)) * TI RI = MAX(BVF2/SHR2,RIMIN) ZK = VK*ZI(I,K+1) ! RL2 = (ZK*RLAM/(RLAM+ZK))**2 ! DK = RL2*SQRT(SHR2) ! RL2 = ZK*RLAM/(RLAM+ZK) ! DK = RL2*RL2*SQRT(SHR2) IF(RI.LT.0.) THEN ! UNSTABLE REGIME RL2 = ZK*RLAMUN/(RLAMUN+ZK) DK = RL2*RL2*SQRT(SHR2) SRI = SQRT(-RI) DKU(i,k) = XKZO(i,k) + DK*(1+8.*(-RI)/(1+1.746*SRI)) DKT(i,k) = XKZO(i,k) + DK*(1+8.*(-RI)/(1+1.286*SRI)) ELSE ! STABLE REGIME RL2 = ZK*RLAM/(RLAM+ZK) ! tem = rlam * sqrt(0.01*prsi(i,k)) ! RL2 = ZK*tem/(tem+ZK) DK = RL2*RL2*SQRT(SHR2) DKT(i,k) = XKZO(i,k) + DK/(1+5.*RI)**2 PRNUM = 1.0 + 2.1*RI PRNUM = MIN(PRNUM,PRMAX) DKU(i,k) = (DKT(i,k)-XKZO(i,k))*PRNUM + XKZO(i,k) ENDIF C DKU(i,k) = MIN(DKU(i,k),DKMAX) DKU(i,k) = MAX(DKU(i,k),DKMIN) DKT(i,k) = MIN(DKT(i,k),DKMAX) DKT(i,k) = MAX(DKT(i,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 ENDDO C C COMPUTE TRIDIAGONAL MATRIX ELEMENTS FOR HEAT AND MOISTURE C DO I=1,IM AD(I,1) = 1. A1(I,1) = T1(i,1) + BETA(i) * HEAT(I) A2(I,1) = Q1(i,1,1) + BETA(i) * EVAP(I) ! A1(I,1) = T1(i,1)-BETAT(I)*(THETA(i,1)-TSEA(I)) ! A2(I,1) = Q1(i,1,1)-BETAQ(I)* ! & (MAX(Q1(i,1,1),QMIN)-MAX(QSS(I),QMIN)) ENDDO if(ntrac.ge.2) then do k = 2, ntrac is = (k-1) * km do i = 1, im A2(I,1+is) = Q1(i,1,k) enddo enddo endif C DO K = 1,KM-1 DO I = 1,IM DTODSD = DT/DEL(I,K) DTODSU = DT/DEL(I,K+1) DSIG = PRSL(I,K)-PRSL(I,K+1) ! RDZ = RDZT(I,K)*2./(T1(i,k)+T1(i,K+1)) RDZ = RDZT(I,K) tem1 = DSIG * DKT(i,k) * RDZ IF(PBLFLG(I).AND.K.LT.KPBL(I)) THEN ! DSDZT = DSIG*DKT(i,k)*RDZ*(GOCP-HGAMT(I)/HPBL(I)) ! DSDZQ = DSIG*DKT(i,k)*RDZ*(-HGAMQ(I)/HPBL(I)) tem = 1.0 / HPBL(I) DSDZT = tem1 * (GOCP-HGAMT(I)*tem) DSDZQ = tem1 * (-HGAMQ(I)*tem) A2(I,k) = A2(I,k)+DTODSD*DSDZQ A2(I,k+1) = Q1(i,k+1,1)-DTODSU*DSDZQ ELSE ! DSDZT = DSIG*DKT(i,k)*RDZ*(GOCP) DSDZT = tem1 * GOCP A2(I,k+1) = Q1(i,k+1,1) ENDIF ! DSDZ2 = DSIG*DKT(i,k)*RDZ*RDZ DSDZ2 = tem1 * RDZ AU(I,k) = -DTODSD*DSDZ2 AL(I,k) = -DTODSU*DSDZ2 AD(I,k) = AD(I,k)-AU(I,k) AD(I,k+1) = 1.-AL(I,k) A1(I,k) = A1(I,k)+DTODSD*DSDZT A1(I,k+1) = T1(i,k+1)-DTODSU*DSDZT ENDDO ENDDO if(ntrac.ge.2) then do kk = 2, ntrac is = (kk-1) * km do k = 1, km - 1 do i = 1, im A2(I,k+1+is) = Q1(i,k+1,kk) enddo enddo enddo endif C C SOLVE TRIDIAGONAL PROBLEM FOR HEAT AND MOISTURE C CALL TRIDIN(IM,KM,ntrac,AL,AD,AU,A1,A2,AU,A1,A2) C C RECOVER TENDENCIES OF HEAT AND MOISTURE C DO K = 1,KM DO I = 1,IM TTEND = (A1(I,k)-T1(i,k))*RDT QTEND = (A2(I,k)-Q1(i,k,1))*RDT TAU(i,k) = TAU(i,k)+TTEND RTG(I,k,1) = RTG(i,k,1)+QTEND DTSFC(I) = DTSFC(I)+CONT*DEL(I,K)*TTEND DQSFC(I) = DQSFC(I)+CONQ*DEL(I,K)*QTEND ENDDO ENDDO if(ntrac.ge.2) then do kk = 2, ntrac is = (kk-1) * km do k = 1, km do i = 1, im QTEND = (A2(I,K+is)-Q1(i,K,kk))*RDT RTG(i,K,kk) = RTG(i,K,kk)+QTEND enddo enddo enddo endif C C COMPUTE TRIDIAGONAL MATRIX ELEMENTS FOR MOMENTUM C DO I=1,IM ! AD(I,1) = 1.+BETAW(I) AD(I,1) = 1.0 + BETA(i) * STRESS(I) / SPD1(I) A1(I,1) = U1(i,1) A2(I,1) = V1(i,1) ! AD(I,1) = 1.0 ! tem = 1.0 + BETA(I) * STRESS(I) / SPD1(I) ! A1(I,1) = U1(i,1) * tem ! A2(I,1) = V1(i,1) * tem ENDDO C DO K = 1,KM-1 DO I=1,IM DTODSD = DT/DEL(I,K) DTODSU = DT/DEL(I,K+1) DSIG = PRSL(I,K)-PRSL(I,K+1) ! RDZ = RDZT(I,K)*2./(T1(i,k)+T1(i,k+1)) RDZ = RDZT(I,K) DSDZ2 = DSIG*DKU(i,k)*RDZ*RDZ AU(I,k) = -DTODSD*DSDZ2 AL(I,k) = -DTODSU*DSDZ2 AD(I,k) = AD(I,k)-AU(I,k) AD(I,k+1) = 1.-AL(I,k) A1(I,k+1) = U1(i,k+1) A2(I,k+1) = V1(i,k+1) ENDDO ENDDO C C SOLVE TRIDIAGONAL PROBLEM FOR MOMENTUM C CALL TRIDI2(IM,KM,AL,AD,AU,A1,A2,AU,A1,A2) C C RECOVER TENDENCIES OF MOMENTUM C DO K = 1,KM DO I = 1,IM CONWRC = CONW*SQRT(RCL(i)) UTEND = (A1(I,k)-U1(i,k))*RDT VTEND = (A2(I,k)-V1(i,k))*RDT DU(i,k) = DU(i,k)+UTEND DV(i,k) = DV(i,k)+VTEND DUSFC(I) = DUSFC(I)+CONWRC*DEL(I,K)*UTEND DVSFC(I) = DVSFC(I)+CONWRC*DEL(I,K)*VTEND ENDDO ENDDO !! RETURN END CFPP$ NOCONCUR R C----------------------------------------------------------------------- SUBROUTINE TRIDI2(L,N,CL,CM,CU,R1,R2,AU,A1,A2) csela %INCLUDE DBTRIDI2; cc USE MACHINE , ONLY : kind_phys implicit none integer k,n,l,i real(kind=kind_phys) fk cc real(kind=kind_phys) CL(L,2:N),CM(L,N),CU(L,N-1),R1(L,N),R2(L,N), & AU(L,N-1),A1(L,N),A2(L,N) C----------------------------------------------------------------------- DO I=1,L FK = 1./CM(I,1) AU(I,1) = FK*CU(I,1) A1(I,1) = FK*R1(I,1) A2(I,1) = FK*R2(I,1) ENDDO DO K=2,N-1 DO I=1,L FK = 1./(CM(I,K)-CL(I,K)*AU(I,K-1)) AU(I,K) = FK*CU(I,K) A1(I,K) = FK*(R1(I,K)-CL(I,K)*A1(I,K-1)) A2(I,K) = FK*(R2(I,K)-CL(I,K)*A2(I,K-1)) ENDDO ENDDO DO I=1,L FK = 1./(CM(I,N)-CL(I,N)*AU(I,N-1)) A1(I,N) = FK*(R1(I,N)-CL(I,N)*A1(I,N-1)) A2(I,N) = FK*(R2(I,N)-CL(I,N)*A2(I,N-1)) ENDDO DO K=N-1,1,-1 DO I=1,L A1(I,K) = A1(I,K)-AU(I,K)*A1(I,K+1) A2(I,K) = A2(I,K)-AU(I,K)*A2(I,K+1) ENDDO ENDDO C----------------------------------------------------------------------- RETURN END CFPP$ NOCONCUR R C----------------------------------------------------------------------- SUBROUTINE TRIDIN(L,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,l,i real(kind=kind_phys) fk(L) cc real(kind=kind_phys) CL(L,2:N), CM(L,N), CU(L,N-1), & R1(L,N), R2(L,N*nt), & AU(L,N-1), A1(L,N), A2(L,N*nt), & FKK(L,2:N-1) C----------------------------------------------------------------------- DO I=1,L FK(I) = 1./CM(I,1) AU(I,1) = FK(I)*CU(I,1) A1(I,1) = FK(I)*R1(I,1) ENDDO do k = 1, nt is = (k-1) * n do i = 1, l a2(i,1+is) = fk(I) * r2(i,1+is) enddo enddo DO K=2,N-1 DO I=1,L FKK(I,K) = 1./(CM(I,K)-CL(I,K)*AU(I,K-1)) AU(I,K) = FKK(I,K)*CU(I,K) A1(I,K) = FKK(I,K)*(R1(I,K)-CL(I,K)*A1(I,K-1)) ENDDO ENDDO do kk = 1, nt is = (kk-1) * n DO K=2,N-1 DO I=1,L A2(I,K+is) = FKK(I,K)*(R2(I,K+is)-CL(I,K)*A2(I,K+is-1)) ENDDO ENDDO ENDDO DO I=1,L FK(I) = 1./(CM(I,N)-CL(I,N)*AU(I,N-1)) A1(I,N) = FK(I)*(R1(I,N)-CL(I,N)*A1(I,N-1)) ENDDO do k = 1, nt is = (k-1) * n do i = 1, l A2(I,N+is) = FK(I)*(R2(I,N+is)-CL(I,N)*A2(I,N+is-1)) enddo enddo DO K=N-1,1,-1 DO I=1,L A1(I,K) = A1(I,K) - AU(I,K)*A1(I,K+1) ENDDO ENDDO do kk = 1, nt is = (kk-1) * n DO K=n-1,1,-1 DO I=1,L A2(I,K+is) = A2(I,K+is) - AU(I,K)*A2(I,K+is+1) ENDDO ENDDO ENDDO C----------------------------------------------------------------------- RETURN END