!WRF:MODEL_LAYER:PHYSICS ! MODULE module_bl_acm ! USE module_model_constants REAL, PARAMETER :: RIC = 0.25 ! critical Richardson number REAL, PARAMETER :: CRANKP = 0.5 ! CRANK-NIC PARAMETER CONTAINS !------------------------------------------------------------------- SUBROUTINE ACMPBL(XTIME, DTPBL, ZNW, SIGMAH, & U3D, V3D, PP3D, DZ8W, TH3D, T3D, & QV3D, QC3D, QI3D, RR3D, & !! Moisture UST, HFX, QFX, TSK, & PSFC, EP1, G, & ROVCP, RD, CPD, & PBLH, KPBL2D, REGIME, & GZ1OZ0, WSPD, PSIM, MUT, & RUBLTEN, RVBLTEN, RTHBLTEN, & RQVBLTEN, RQCBLTEN, RQIBLTEN, & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte) !********************************************************************** ! THIS MODULE COMPUTES VERTICAL MIXING IN AND ABOVE THE PBL ACCORDING TO ! THE ASYMMETRICAL CONVECTIVE MODEL, VERSION 2 (ACM2), WHICH IS A COMBINED ! LOCAL NON-LOCAL CLOSURE SCHEME BASED ON THE ORIGINAL ACM (PLEIM AND CHANG 1992) ! ! REFERENCES: ! Pleim (2007) A combined local and non-local closure model for the atmospheric ! boundary layer. Part1: Model description and testing. ! JAMC, 46, 1383-1395 ! Pleim (2007) A combined local and non-local closure model for the atmospheric ! boundary layer. Part2: Application and evaluation in a mesoscale ! meteorology model. JAMC, 46, 1396-1409 ! ! REVISION HISTORY: ! AX 3/2005 - developed WRF version based on the MM5 PX LSM ! RG and JP 7/2006 - Finished WRF adaptation ! !********************************************************************** ! ARGUMENT LIST: ! !... Inputs: !-- XTIME Time since simulation start (min) !-- DTPBL PBL time step !-- ZNW Sigma at full layer !-- SIGMAH Sigma at half layer !-- U3D 3D u-velocity interpolated to theta points (m/s) !-- V3D 3D v-velocity interpolated to theta points (m/s) !-- PP3D Pressure at half levels (Pa) !-- DZ8W dz between full levels (m) !-- TH3D Potential Temperature (K) !-- T3D Temperature (K) !-- QV3D 3D water vapor mixing ratio (Kg/Kg) !-- QC3D 3D cloud mixing ratio (Kg/Kg) !-- QI3D 3D ice mixing ratio (Kg/Kg) !-- RR3D 3D dry air density (kg/m^3) !-- UST Friction Velocity (m/s) !-- HFX Upward heat flux at the surface (w/m^2) !-- QFX Upward moisture flux at the surface (Kg/m^2/s) !-- TSK Surface temperature (K) !-- PSFC Pressure at the surface (Pa) !-- EP1 Constant for virtual temperature (r_v/r_d-1) (dimensionless) !-- G Gravity (m/s^2) !-- ROVCP r/cp !-- RD gas constant for dry air (j/kg/k) !-- CPD heat capacity at constant pressure for dry air (j/kg/k) !-- GZ1OZ0 log(z/z0) where z0 is roughness length !-- WSPD wind speed at lowest model level (m/s) !-- PSIM similarity stability function for momentum !-- MUT Total Mu : Psfc - Ptop !-- ids start index for i in domain !-- ide end index for i in domain !-- jds start index for j in domain !-- jde end index for j in domain !-- kds start index for k in domain !-- kde end index for k in domain !-- ims start index for i in memory !-- ime end index for i in memory !-- jms start index for j in memory !-- jme end index for j in memory !-- kms start index for k in memory !-- kme end index for k in memory !-- jts start index for j in tile !-- jte end index for j in tile !-- kts start index for k in tile !-- kte end index for k in tile ! !... Outputs: !-- PBLH PBL height (m) !-- KPBL2D K index for PBL layer !-- REGIME Flag indicating PBL regime (stable, unstable, etc.) !-- RUBLTEN U tendency due to PBL parameterization (m/s^2) !-- RVBLTEN V tendency due to PBL parameterization (m/s^2) !-- RTHBLTEN Theta tendency due to PBL parameterization (K/s) !-- RQVBLTEN Qv tendency due to PBL parameterization (kg/kg/s) !-- RQCBLTEN Qc tendency due to PBL parameterization (kg/kg/s) !-- RQIBLTEN Qi tendency due to PBL parameterization (kg/kg/s) ! !------------------------------------------------------------------- IMPLICIT NONE !.......Arguments ! DECLARATIONS - INTEGER INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte, XTIME ! DECLARATIONS - REAL REAL, INTENT(IN) :: DTPBL, EP1, & G, ROVCP, RD, CPD REAL, DIMENSION( kms:kme ), INTENT(IN) :: ZNW, SIGMAH REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), & INTENT(IN) :: U3D, V3D, & PP3D, DZ8W, T3D, & QV3D, QC3D, QI3D, & RR3D, TH3D REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: UST REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: & HFX, QFX, TSK, & PSFC INTEGER, DIMENSION( ims:ime, jms:jme ), & INTENT(OUT ) :: KPBL2D REAL, DIMENSION( ims:ime, jms:jme ), & INTENT(INOUT) :: PBLH REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: REGIME REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: & psim, & gz1oz0, & wspd , mut REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), & INTENT(INOUT) :: RUBLTEN, RVBLTEN, & RTHBLTEN, RQVBLTEN, & RQCBLTEN, RQIBLTEN !... Local Variables !... Integer INTEGER :: J, K !... Real REAL, DIMENSION( kts:kte ) :: DSIGH, DSIGHI, DSIGFI REAL, DIMENSION( 0:kte ) :: SIGMAF REAL RDT REAL, PARAMETER :: KARMAN = 0.4 !... RDT = 1.0 / DTPBL DO K = 1, kte SIGMAF(K-1) = ZNW(K) ENDDO SIGMAF(kte) = 0.0 DO K = 1, kte DSIGH(K) = SIGMAF(K) - SIGMAF(K-1) DSIGHI(K) = 1.0 / DSIGH(K) ENDDO DO K = kts,kte-1 DSIGFI(K) = 1.0 / (SIGMAH(K+1) - SIGMAH(K)) ENDDO DSIGFI(kte) = DSIGFI(kte-1) do j = jts,jte CALL ACM2D(j=J,xtime=XTIME, dtpbl=DTPBL, sigmaf=SIGMAF, sigmah=SIGMAH & ,dsigfi=DSIGFI,dsighi=DSIGHI,dsigh=DSIGH & ,us=u3d(ims,kms,j),vs=v3d(ims,kms,j) & ,theta=th3d(ims,kms,j),tt=t3d(ims,kms,j) & ,qvs=qv3d(ims,kms,j),qcs=qc3d(ims,kms,j) & ,qis=qi3d(ims,kms,j),dzf=DZ8W(ims,kms,j) & ,densx=RR3D(ims,kms,j) & ,utnp=rublten(ims,kms,j),vtnp=rvblten(ims,kms,j) & ,ttnp=rthblten(ims,kms,j),qvtnp=rqvblten(ims,kms,j) & ,qctnp=rqcblten(ims,kms,j),qitnp=rqiblten(ims,kms,j) & ,cpd=cpd,g=g,rovcp=rovcp,rd=rd,rdt=rdt & ,psfcpa=psfc(ims,j),ust=ust(ims,j) & ,pbl=pblh(ims,j) & ,regime=regime(ims,j),psim=psim(ims,j) & ,hfx=hfx(ims,j),qfx=qfx(ims,j) & ,tg=tsk(ims,j),gz1oz0=gz1oz0(ims,j) & ,wspd=wspd(ims,j) ,klpbl=kpbl2d(ims,j) & ,mut=mut(ims,j) & ,ep1=ep1,karman=karman & ,ids=ids,ide=ide, jds=jds,jde=jde, kds=kds,kde=kde & ,ims=ims,ime=ime, jms=jms,jme=jme, kms=kms,kme=kme & ,its=its,ite=ite, jts=jts,jte=jte, kts=kts,kte=kte ) enddo END SUBROUTINE ACMPBL !-------------------------------------------------------------------- SUBROUTINE ACM2D(j,XTIME, DTPBL, sigmaf, sigmah & ,dsigfi,dsighi,dsigh & ,us,vs,theta,tt,qvs,qcs,qis & ,dzf,densx,utnp,vtnp,ttnp,qvtnp,qctnp,qitnp & ,cpd,g,rovcp,rd,rdt,psfcpa,ust & ,pbl,regime,psim & ,hfx,qfx,tg,gz1oz0,wspd ,klpbl & ,mut & ,ep1,karman & ,ids,ide, jds,jde, kds,kde & ,ims,ime, jms,jme, kms,kme & ,its,ite, jts,jte, kts,kte ) IMPLICIT NONE !.......Arguments !... Real REAL, DIMENSION( 0:kte ), INTENT(IN) :: SIGMAF REAL, DIMENSION( kms:kme ), INTENT(IN) :: SIGMAH REAL, DIMENSION( kts:kte ), INTENT(IN) :: DSIGH, DSIGHI, DSIGFI REAL , INTENT(IN) :: DTPBL, G, RD,ep1,karman,CPD,ROVCP,RDT REAL , DIMENSION( ims:ime ), INTENT(INOUT) :: PBL, UST REAL , DIMENSION( ims:ime, kms:kme ), INTENT(IN) :: US,VS, THETA, TT, & QVS, QCS, QIS, DENSX REAL, DIMENSION( ims:ime, kms:kme ), intent(in) :: DZF REAL, DIMENSION( ims:ime, kms:kme ), intent(inout) :: utnp, & vtnp, & ttnp, & qvtnp, & qctnp, & qitnp real, dimension( ims:ime ), intent(in ) :: psfcpa real, dimension( ims:ime ), intent(in ) :: tg real, dimension( ims:ime ), intent(inout) :: regime real, dimension( ims:ime ), intent(in) :: wspd, psim, gz1oz0 real, dimension( ims:ime ), intent(in) :: hfx, qfx real, dimension( ims:ime ), intent(in) :: mut !... Integer INTEGER, DIMENSION( ims:ime ), INTENT(OUT):: KLPBL INTEGER, INTENT(IN) :: XTIME integer, intent(in ) :: ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte, j !-------------------------------------------------------------------- !--Local INTEGER I, K INTEGER :: KPBLHT INTEGER, DIMENSION( its:ite ) :: KPBLH, NOCONV !... Real REAL :: TVCON, WSS, TCONV, TH1, TOG, DTMP, WSSQ REAL :: psix, THV1 REAL, DIMENSION( its:ite ) :: FINT, PSTAR, CPAIR REAL, DIMENSION( its:ite, kts:kte ) :: THETAV, RIB, & EDDYZ, UX, VX, THETAX, & QVX, QCX, QIX, ZA REAL, DIMENSION( its:ite, 0:kte ) :: ZF REAL, DIMENSION( its:ite) :: WST, TST, QST, USTM, TSTV REAL, DIMENSION( its:ite ) :: PBLSIG, MOL REAL :: FINTT, ZMIX, UMIX, VMIX REAL :: TMPFX, TMPVTCON, TMPP, TMPTHS, TMPTH1, TMPVCONV, WS1, DTH REAL :: A,TST12,RL,ZFUNC ! REAL, PARAMETER :: KARMAN = 0.4 !... Integer INTEGER :: KL, jtf, ktf, itf, KMIX !... !-----initialize vertical tendencies and ! do i = its,ite do k = kts,kte utnp(i,k) = 0. vtnp(i,k) = 0. ttnp(i,k) = 0. enddo enddo ! do k = kts,kte do i = its,ite qvtnp(i,k) = 0. enddo enddo ! do k = kts,kte do i = its,ite qctnp(i,k) = 0. qitnp(i,k) = 0. enddo enddo ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!! Compute Micromet Scaling variables, not availiable in WRF for ACM DO I = its,ite cpair(i) = CPD * (1.0 + 0.84 * QVS(I,1)) ! J/(K KG) TMPFX = HFX(I) / (cpair(i) * DENSX(I,1)) TMPVTCON = 1.0 + EP1 * QVS(I,1) ! COnversion factor for virtual temperature WS1 = SQRT(US(I,1)**2 + VS(I,1)**2) ! Level 1 wind speed TST(I) = -TMPFX / UST(I) QST(I) = -QFX(I) / (UST(I)*DENSX(I,1)) USTM(I) = UST(I) * WS1 / wspd(i) THV1 = TMPVTCON * THETA(I,1) TSTV(I) = TST(I)*TMPVTCON + THV1*EP1*QST(I) MOL(I) = THV1 * UST(i)**2/(KARMAN*G*TSTV(I)) WST(I) = UST(I) * (PBL(I)/(KARMAN*ABS(MOL(I)))) ** 0.333333 PSTAR(I) = MUT(I)/1000. ! P* in cb ENDDO !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !... Compute PBL height !... compute the height of full- and half-sigma level above ground level DO I = its,ite ZF(I,0) = 0.0 KLPBL(I) = 1 ENDDO DO K = kts, kte DO I = its,ite ZF(I,K) = DZF(I,K) + ZF(I,K-1) ZA(I,K) = 0.5 * (ZF(I,K) + ZF(I,K-1)) ENDDO ENDDO DO K = kts, kte DO I = its,ite TVCON = 1.0 + EP1 * QVS(I,K) THETAV(I,K) = THETA(I,K) * TVCON ENDDO ENDDO !... COMPUTE PBL WHERE RICHARDSON NUMBER = RIC (0.25) HOLTSLAG ET AL 1990 DO 100 I = its,ite if(MOL(I).LT.0.0) then WSS = (UST(I) ** 3 + 0.6 * WST(I) ** 3) ** 0.33333 TMPFX = -UST(I) * TSTV(I) TCONV = 8.5 * TMPFX / WSS TH1 = THETAV(I,1) + TCONV else TH1 = THETAV(I,1) endif 99 KMIX = 1 DO K = 1,kte DTMP = THETAV(I,K) - TH1 IF (DTMP.LT.0.0) KMIX = K ENDDO IF(KMIX.GT.1) THEN FINTT = (TH1 - THETAV(I,KMIX)) / (THETAV(I,KMIX+1) & - THETAV(I,KMIX)) ZMIX = FINTT * (ZA(I,KMIX+1)-ZA(I,KMIX)) + ZA(I,KMIX) UMIX = FINTT * (US(I,KMIX+1)-US(I,KMIX)) + US(I,KMIX) VMIX = FINTT * (VS(I,KMIX+1)-VS(I,KMIX)) + VS(I,KMIX) ELSE ZMIX = ZA(I,1) UMIX = US(I,1) VMIX = VS(I,1) ENDIF DO K = KMIX,kte DTMP = THETAV(I,K) - TH1 TOG = 0.5 * (THETAV(I,K) + TH1) / G WSSQ = (US(I,K)-UMIX)**2 & + (VS(I,K)-VMIX)**2 IF (KMIX == 1) WSSQ = WSSQ + 100.*UST(I)*UST(I) WSSQ = MAX( WSSQ, 0.1 ) RIB(I,K) = ABS(ZA(I,K)-ZMIX) * DTMP / (TOG * WSSQ) IF (RIB(I,K) .GE. RIC) GO TO 201 ENDDO print *,' RIB never exceeds RIC, RIB(i,kte) = ',rib(i,5), & ' THETAV(i,1) = ',thetav(i,1),' MOL=',mol(i), & ' TCONV = ',TCONV,' WST = ',WST(I), & ' KMIX = ',kmix,' UST = ',UST(I), & ' TST = ',TST(I),' U,V = ',US(I,1),VS(I,1), & ' I,J=',I,J STOP 201 CONTINUE KPBLH(I) = K 100 CONTINUE DO I = its,ite IF (KPBLH(I) .NE. 1) THEN !---------INTERPOLATE BETWEEN LEVELS -- jp 7/93 FINT(I) = (RIC - RIB(I,KPBLH(I)-1)) / (RIB(I,KPBLH(I)) - & RIB(I,KPBLH(I)-1)) IF (FINT(I) .GT. 0.5) THEN KPBLHT = KPBLH(I) FINT(I) = FINT(I) - 0.5 ELSE KPBLHT = KPBLH(I) - 1 FINT(I) = FINT(I) + 0.5 ENDIF PBL(I) = FINT(I) * (ZF(I,KPBLHT) - ZF(I,KPBLHT-1)) + & ZF(I,KPBLHT-1) KLPBL(I) = KPBLHT PBLSIG(I) = FINT(I) * DSIGH(KPBLHT) + SIGMAF(KPBLHT-1) ! sigma at PBL height ELSE KLPBL(I) = 1 PBL(I) = ZF(I,1) PBLSIG(I) = SIGMAF(1) ENDIF ENDDO DO I = its,ite NOCONV(I) = 0 ! Check for CBL and identify conv. vs. non-conv cells IF (PBL(I) / MOL(I) .LT. -0.02 .AND. KLPBL(I) .GT. 3 & .AND. THETAV(I,1) .GT. THETAV(I,2) .AND. XTIME .GT. 1) THEN NOCONV(I) = 1 REGIME(I) = 4.0 ! FREE CONVECTIVE - ACM ENDIF ENDDO !... Calculate Kz CALL EDDYX(DTPBL, ZF, ZA, MOL, PBL, UST, & US, VS, TT, THETAV, DENSX, PSTAR, & QVS, QCS, QIS, DSIGFI, G, RD, CPAIR, & EDDYZ, its,ite, kts,kte,ims,ime, kms,kme) CALL ACM(DTPBL, PSTAR, NOCONV, SIGMAF, DSIGH, DSIGHI, J, & KLPBL, PBL, PBLSIG, MOL, UST, & TST, QST, USTM, EDDYZ, DENSX, & US, VS, THETA, QVS, QCS, QIS, & UX, VX, THETAX, QVX, QCX, QIX, & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte) !... Calculate tendency due to PBL parameterization DO K = kts, kte DO I = its, ite UTNP(I,K) = UTNP(I,K) + (UX(I,K) - US(I,K)) * RDT VTNP(I,K) = VTNP(I,K) + (VX(I,K) - VS(I,K)) * RDT TTNP(I,K) = TTNP(I,K) + (THETAX(I,K) - THETA(I,K)) * RDT QVTNP(I,K) = QVTNP(I,K) + (QVX(I,K) - QVS(I,K)) * RDT QCTNP(I,K) = QCTNP(I,K) + (QCX(I,K) - QCS(I,K)) * RDT QITNP(I,K) = QITNP(I,K) + (QIX(I,K) - QIS(I,K)) * RDT ENDDO ENDDO ! IF(J.EQ.36)PRINT *,' PBL,THETA,THETAX,QVS,QVX=',PBL(20),THETA(20,1),THETAX(20,1),QVS(20,1),QVX(20,1) ! IF(J.EQ.36)PRINT *,' UST,ustm,TG,TST=',UST(20),ustm(20),Tg(20),TST(20) ! IF(J.EQ.36)PRINT *,' HFX,MOL,WST,WSPD=',HFX(20),MOL(20),WST(20),wspd(20) ! IF(J.EQ.36)then ! i=20 ! do k=1,kte ! PRINT *,' qvten,uten,vten=',QVTNP(I,K),UTNP(I,K),VTNP(I,K) ! PRINT *,' k,thten,th,thx,edy=',k,TTNP(I,K),THETA(20,k),THETAX(20,k),EDDYZ(20,K) ! enddo ! ENDIF END SUBROUTINE ACM2D !------------------------------------------------------------------- SUBROUTINE ACMINIT(RUBLTEN,RVBLTEN,RTHBLTEN,RQVBLTEN, & RQCBLTEN,RQIBLTEN,P_QI,P_FIRST_SCALAR, & restart, allowed_to_read , & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) !********************************************************************** ! ! This subroutine is for preparing ACM PBL variables. ! Called from module_physics_init.F ! ! REVISION HISTORY: ! AX 3/2005 - Originally developed !********************************************************************** ! ARGUMENT LIST: ! !------------------------------------------------------------------- IMPLICIT NONE ! LOGICAL , INTENT(IN) :: restart , allowed_to_read INTEGER , INTENT(IN) :: ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte INTEGER , INTENT(IN) :: P_QI,P_FIRST_SCALAR ! REAL , DIMENSION( kms:kme ), INTENT(IN) :: SHALF REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT) :: & RUBLTEN, & RVBLTEN, & RTHBLTEN, & RQVBLTEN, & RQCBLTEN, & RQIBLTEN !... Local Variables INTEGER :: i, j, k, itf, jtf, ktf ! jtf=min0(jte,jde-1) ktf=min0(kte,kde-1) itf=min0(ite,ide-1) IF(.not.restart)THEN DO j=jts,jtf DO k=kts,ktf DO i=its,itf RUBLTEN(i,k,j)=0. RVBLTEN(i,k,j)=0. RTHBLTEN(i,k,j)=0. RQVBLTEN(i,k,j)=0. RQCBLTEN(i,k,j)=0. ENDDO ENDDO ENDDO ENDIF IF (P_QI .ge. P_FIRST_SCALAR .and. .not.restart) THEN DO j=jts,jtf DO k=kts,ktf DO i=its,itf RQIBLTEN(i,k,j)=0. ENDDO ENDDO ENDDO ENDIF END SUBROUTINE acminit !------------------------------------------------------------------- SUBROUTINE EDDYX(DTPBL, ZF, ZA, MOL, PBL, UST, & US, VS, TT, THETAV, DENSX, PSTAR, & QVS, QCS, QIS, DSIGFI, G, RD, CPAIR, & EDDYZ, its,ite, kts,kte,ims,ime,kms,kme ) !********************************************************************** ! Two methods for computing Kz: ! 1. Boundary scaling similar to Holtslag and Boville (1993) ! 2. Local Kz computed as function of local Richardson # and vertical ! wind shear, similar to LIU & CARROLL (1996) ! !********************************************************************** ! !-- DTPBL time step of the minor loop for the land-surface/pbl model !-- ZF height of full sigma level !-- ZA height of half sigma level !-- MOL Monin-Obukhov length in 1D form !-- PBL PBL height in 1D form !-- UST friction velocity U* in 1D form (m/s) !-- US U wind !-- VS V wind !-- TT temperature !-- THETAV potential virtual temperature !-- DENSX dry air density (kg/m^3) !-- PSTAR P*=Psfc-Ptop !-- QVS water vapor mixing ratio (Kg/Kg) !-- QCS cloud mixing ratio (Kg/Kg) !-- QIS ice mixing ratio (Kg/Kg) !-- DSIGFI inverse of sigma layer delta !-- G gravity !-- RD gas constant for dry air (j/kg/k) !-- CPAIR specific heat of moist air (M^2 S^-2 K^-1) !-- EDDYZ eddy diffusivity KZ !----------------------------------------------------------------------- IMPLICIT NONE !.......Arguments !... Integer INTEGER, INTENT(IN) :: its,ite, kts,kte,ims,ime,kms,kme !... Real REAL , DIMENSION( ims:ime ), INTENT(IN) :: PBL, UST REAL , INTENT(IN) :: DTPBL, G, RD REAL , DIMENSION( kts:kte ), INTENT(IN) :: DSIGFI REAL , DIMENSION( its:ite ), INTENT(IN) :: MOL, PSTAR, CPAIR REAL , DIMENSION( ims:ime, kms:kme ), INTENT(IN) :: US,VS, TT, & QVS, QCS, QIS, DENSX REAL, DIMENSION( its:ite, kts:kte ), INTENT(IN) :: ZA, THETAV REAL, DIMENSION( its:ite, 0:kte ) , INTENT(IN) :: ZF REAL , DIMENSION( its:ite, kts:kte ), INTENT(OUT) :: EDDYZ !.......Local variables !... Integer INTEGER :: ILX, KL, KLM, K, I !... Real REAL :: ZOVL, PHIH, WT, ZSOL, ZFUNC, DZF, SS, GOTH, EDYZ REAL :: RI, QMEAN, TMEAN, XLV, ALPH, CHI, ZK, SQL, DENSF,kzo !... Parameters REAL, PARAMETER :: RV = 461.5 REAL, PARAMETER :: RC = 0.25 REAL, PARAMETER :: RLAM = 80.0 REAL, PARAMETER :: GAMH = 16.0 !15.0 ! Holtslag and Boville (1993) REAL, PARAMETER :: BETAH = 5.0 ! Holtslag and Boville (1993) REAL, PARAMETER :: KARMAN = 0.4 REAL, PARAMETER :: EDYZ0 = 0.01 ! New Min Kz ! REAL, PARAMETER :: EDYZ0 = 0.1 !-- IMVDIF imvdif=1 for moist adiabat vertical diffusion INTEGER, PARAMETER :: imvdif = 1 ! ILX = ite KL = kte KLM = kte - 1 DO K = kts,KLM DO I = its,ILX EDYZ = 0.0 ZOVL = 0.0 DZF = ZA(I,K+1) - ZA(I,K) ! kzo = min(0.001 * dzf,EDYZ0) kzo = 0.001 * dzf ! kzo = EDYZ0 !------------------------------------------------- IF (ZF(I,K) .LT. PBL(I)) THEN ZOVL = ZF(I,K) / MOL(I) IF (ZOVL .LT. 0.0) THEN IF (ZF(I,K) .LT. 0.1 * PBL(I)) THEN PHIH = 1.0 / SQRT(1.0 - GAMH * ZOVL) WT = UST(I) / PHIH ELSE ZSOL = 0.1 * PBL(I) / MOL(I) PHIH = 1.0 / SQRT(1.0 - GAMH * ZSOL) WT = UST(I) / PHIH ENDIF ELSE IF (ZOVL .LT. 1.0) THEN PHIH = 1.0 + BETAH * ZOVL WT = UST(I) / PHIH ELSE PHIH = BETAH + ZOVL WT = UST(I) / PHIH ENDIF ZFUNC = ZF(I,K) * (1.0 - ZF(I,K) / PBL(I)) ** 2 EDYZ = KARMAN * WT * ZFUNC ! EDYZ = AMAX1(EDYZ,EDYZ0) ENDIF ! PRINT *,I,K,TT(I,K) !-------------------------------------------------------------------------- ! ! RC = 0.257 * DZF ** 0.175 SS = ((US(I,K+1) - US(I,K)) ** 2 + (VS(I,K+1) - VS(I,K)) ** 2) & / (DZF * DZF) + 1.0E-9 GOTH = 2.0 * G / (THETAV(I,K+1) + THETAV(I,K)) RI = GOTH * (THETAV(I,K+1) - THETAV(I,K)) / (DZF * SS) ! PRINT *,I,K,TT(I,K), RI, IMVDIF ! !-- Adjustment to vert diff in Moist air if(imvdif.eq.1)then IF ((QCS(I,K)+QIS(I,K)) .GT. 0.01E-3 .OR. (QCS(I,K+1)+ & QIS(I,K+1)) .GT. 0.01E-3) THEN QMEAN = 0.5 * (QVS(I,K) + QVS(I,K+1)) TMEAN = 0.5 * (TT(I,K) + TT(I,K+1)) XLV = (2.501 - 0.00237 * (TMEAN - 273.15)) * 1.E6 ALPH = XLV * QMEAN / RD / TMEAN CHI = XLV * XLV * QMEAN / CPAIR(I) / RV / TMEAN / TMEAN RI = (1.0 + ALPH) * (RI -G * G / SS / TMEAN / CPAIR(I) * & ((CHI - ALPH) / (1.0 + CHI))) ENDIF endif ZK = 0.4 * ZF(I,K) SQL = (ZK * RLAM / (RLAM + ZK)) ** 2 IF (RI .GE. RC) THEN EDDYZ(I,K) = kzo !EDYZ0 ELSE IF (RI .GE. 0.0) THEN EDDYZ(I,K) = kzo + SQRT(SS) * (1.- RI / RC) ** 2 * SQL ELSE EDDYZ(I,K) = kzo + SQRT(SS * (1.0 - 25.0 * RI)) * SQL ENDIF IF(EDYZ.GT.EDDYZ(I,K)) THEN EDDYZ(I,K) = EDYZ ENDIF EDDYZ(I,K) = MIN(1000.0,EDDYZ(I,K)) EDDYZ(I,K) = MAX(kzo,EDDYZ(I,K)) DENSF = 0.5 * (DENSX(I,K+1) + DENSX(I,K)) EDDYZ(I,K) = EDDYZ(I,K) * (DENSF * G / PSTAR(I)) ** 2 * & DTPBL * DSIGFI(K)*1E-6 ENDDO ! for I loop ENDDO ! for k loop ! DO I = its,ILX EDDYZ(I,KL) = 0.0 ! EDDYZ(I,KLM) -- changed jp 3/08 ENDDO END SUBROUTINE EDDYX !------------------------------------------------------------------- SUBROUTINE ACM (DTPBL, PSTAR, NOCONV, SIGMAF, DSIGH, DSIGHI, JX, & KLPBL, PBL, PBLSIG, MOL, UST, & TST, QST, USTM, EDDYZ, DENSX, & US, VS, THETA, QVS, QCS, QIS, & UX, VX, THETAX, QVX, QCX, QIX, & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte) !********************************************************************** ! PBL model called the Asymmetric Convective Model, Version 2 (ACM2) ! -- See top of module for summary and references ! !---- REVISION HISTORY: ! AX 3/2005 - developed WRF version based on ACM2 in the MM5 PX LSM ! JP and RG 8/2006 - updates ! !********************************************************************** ! ARGUMENTS: !-- DTPBL PBL time step !-- PSTAR Psurf - Ptop in cb !-- NOCONV If free convection =0, no; =1, yes !-- SIGMAF Sigma for full layer !-- DSIGH Sigma thickness !-- DSIGHI Inverse of sigma thickness !-- JX N-S index !-- KLPBL PBL level at K index !-- PBL PBL height in m !-- PBLSIG Sigma level for PBL !-- MOL Monin-Obukhov length in 1D form !-- UST U* in 1D form !-- TST Theta* in 1D form !-- QST Q* in 1D form !-- USTM U* for computation of momemtum flux !-- EDDYZ eddy diffusivity KZ !-- DENSX dry air density (kg/m^3) !-- US U wind !-- VS V wind !-- THETA potential temperature !-- QVS water vapor mixing ratio (Kg/Kg) !-- QCS cloud mixing ratio (Kg/Kg) !-- QIS ice mixing ratio (Kg/Kg) !-- UX new U wind !-- VX new V wind !-- THETAX new potential temperature !-- QVX new water vapor mixing ratio (Kg/Kg) !-- QCX new cloud mixing ratio (Kg/Kg) !-- QIX new ice mixing ratio (Kg/Kg) !----------------------------------------------------------------------- IMPLICIT NONE !.......Arguments !... Integer INTEGER, INTENT(IN) :: ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte, JX INTEGER, DIMENSION( its:ite ), INTENT(IN) :: NOCONV INTEGER, DIMENSION( ims:ime ), INTENT(IN) :: KLPBL !... Real REAL , DIMENSION( ims:ime ), INTENT(IN) :: PBL, UST REAL , INTENT(IN) :: DTPBL REAL , DIMENSION( its:ite ), INTENT(IN) :: PSTAR, PBLSIG, & MOL, TST, & QST, USTM REAL , DIMENSION( kts:kte ), INTENT(IN) :: DSIGHI, DSIGH REAL , DIMENSION( 0:kte ), INTENT(IN) :: SIGMAF REAL , DIMENSION( its:ite, kts:kte ), INTENT(INOUT) :: EDDYZ REAL , DIMENSION( ims:ime, kms:kme ), INTENT(IN) :: US,VS, THETA, & QVS, QCS, QIS, DENSX REAL , DIMENSION( its:ite, kts:kte ), INTENT(OUT) :: UX, VX, THETAX, & QVX, QCX, QIX !.......Local variables !... Parameters INTEGER, PARAMETER :: NSP = 6 ! !......ACM2 Parameters ! INTEGER, PARAMETER :: IFACM = 0 ! REAL, PARAMETER :: G1000 = 9.8 * 1.0E-3 REAL, PARAMETER :: XX = 0.5 ! FACTOR APPLIED TO CONV MIXING TIME STEP REAL, PARAMETER :: KARMAN = 0.4 !... Integer INTEGER :: ILX, KL, KLM, I, K, NSPX, NLP, NL, JJ, L INTEGER :: KCBLMX INTEGER, DIMENSION( its:ite ) :: KCBL !... Real REAL :: G1000I, MBMAX, HOVL, MEDDY, MBAR REAL :: EKZ, RZ, FM, WSPD, DTS, DTRAT, F1 REAL, DIMENSION( its:ite ) :: PSTARI, FSACM, RAH, DTLIM REAL, DIMENSION( kts:kte, its:ite) :: MBARKS, MDWN REAL, DIMENSION( 1:NSP, its:ite ) :: FS, BCBOTN REAL, DIMENSION( kts:kte ) :: XPLUS, XMINUS REAL DELC REAL, DIMENSION( 1:NSP,its:ite,kts:kte ) :: VCI REAL, DIMENSION( kts:kte ) :: AI, BI, CI, EI !, Y REAL, DIMENSION( 1:NSP,kts:kte ) :: DI, UI ! !--Start Exicutable ---- ILX = ite KL = kte KLM = kte - 1 G1000I = 1.0 / G1000 KCBLMX = 0 MBMAX = 0.0 !---COMPUTE ACM MIXING RATE DO I = its, ILX DTLIM(I) = DTPBL PSTARI(I) = 1.0 / PSTAR(I) KCBL(I) = 1 FSACM(I) = 0.0 IF (NOCONV(I) .EQ. 1) THEN KCBL(I) = KLPBL(I) !-------MBARKS IS UPWARD MIXING RATE; MDWN IS DOWNWARD MIXING RATE !--New couple ACM & EDDY------------------------------------------------------------- HOVL = -PBL(I) / MOL(I) FSACM(I) = 1./(1.+((KARMAN/(HOVL))**0.3333)/(0.72*KARMAN)) MEDDY = EDDYZ(I,1) / (DTPBL * (PBLSIG(I) - SIGMAF(1))) MBAR = MEDDY * FSACM(I) DO K = kts,KCBL(I)-1 EDDYZ(I,K) = EDDYZ(I,K) * (1.0 - FSACM(I)) ENDDO ! if(i.eq.100.and.jx.eq.43) PRINT *,' Edy,meddy,mbar=', EDDYZ(I,1),MEDDY,MBAR MBMAX = AMAX1(MBMAX,MBAR) DO K = kts+1,KCBL(I) MBARKS(K,I) = MBAR MDWN(K,I) = MBAR * (PBLSIG(I) - SIGMAF(K-1)) * DSIGHI(K) ENDDO MBARKS(1,I) = MBAR MBARKS(KCBL(I),I) = MDWN(KCBL(I),I) MDWN(KCBL(I)+1,I) = 0.0 ENDIF ENDDO ! end of I loop ! if(NOCONV(20).EQ.1) print *,' KCBL,PBLSIG=',KCBL(20),PBLSIG(20),SIGMAF(KCBL(20)-1),DSIGHI(KCBL(20)) ! do k = kts,klm ! if(NOCONV(20).EQ.1.and.MBAR.lt.1.0e-3)print *,' k,MBARKS,MDWN=',k,MBARKS(k,20),MDWN(k,20) ! enddo DO K = kts,KLM DO I = its,ILX EKZ = EDDYZ(I,K) / DTPBL * DSIGHI(K) DTLIM(I) = AMIN1(0.75 / EKZ,DTLIM(I)) ENDDO ENDDO DO I = its,ILX IF (NOCONV(I) .EQ. 1) THEN KCBLMX = AMAX0(KLPBL(I),KCBLMX) RZ = (SIGMAF(KCBL(I)) - SIGMAF(1)) * DSIGHI(1) DTLIM(I) = AMIN1(XX / (MBARKS(1,I) * RZ),DTLIM(I)) ENDIF ENDDO DO I = its,ILX ENDDO ! DO K = kts,KL DO I = its,ILX VCI(1,I,K) = THETA(I,K) VCI(2,I,K) = QVS(I,K) VCI(3,I,K) = US(I,K) VCI(4,I,K) = VS(I,K) ! -- Also mix cloud water and ice if necessary ! IF (IMOISTX.NE.1.AND.IMOISTX.NE.3) THEN !!! Check other PBL models VCI(5,I,K) = QCS(I,K) VCI(6,I,K) = QIS(I,K) ENDDO ENDDO NSPX=6 DO I = its,ILX ! RAH(I) = RA(I,J) + 3.976 / UST(I) FS(1,I) = -UST(I) * TST(I) * DENSX(I,1) * PSTARI(I) FS(2,I) = -UST(I) * QST(I) * DENSX(I,1) * PSTARI(I) FM = -USTM(I) * USTM(I) * DENSX(I,1) * PSTARI(I) WSPD = SQRT(US(I,1) * US(I,1) + VS(I,1) * VS(I,1)) + 1.E-9 FS(3,I) = FM * US(I,1) / WSPD FS(4,I) = FM * VS(I,1) / WSPD FS(5,I) = 0.0 FS(6,I) = 0.0 ! SURFACE FLUXES OF CLOUD WATER AND ICE = 0 ENDDO ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! DO I = its,ILX !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! NLP = INT(DTPBL / DTLIM(I) + 1.0) DTS = (DTPBL / NLP) DTRAT = DTS / DTPBL DO NL = 1,NLP ! LOOP OVER SUB TIME LOOP ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !-- COMPUTE ARRAY ELEMENTS THAT ARE INDEPENDANT OF SPECIES DO K = kts,KL AI(K) = 0.0 BI(K) = 0.0 CI(K) = 0.0 EI(K) = 0.0 ENDDO DO K = 2, KCBL(I) EI(K-1) = -CRANKP * MDWN(K,I) * DTS * DSIGH(K) * DSIGHI(K-1) BI(K) = 1.0 + CRANKP * MDWN(K,I) * DTS AI(K) = -CRANKP * MBARKS(K,I) * DTS ENDDO EI(1) = EI(1) -EDDYZ(I,1) * CRANKP * DSIGHI(1 )* DTRAT AI(2) = AI(2) -EDDYZ(I,1) * CRANKP * DSIGHI(2) * DTRAT DO K = KCBL(I)+1, KL BI(K) = 1.0 ENDDO DO K = 2,KL XPLUS(K) = EDDYZ(I,K) * DSIGHI(K) * DTRAT XMINUS(K) = EDDYZ(I,K-1) * DSIGHI(K) * DTRAT CI(K) = - XMINUS(K) * CRANKP EI(K) = EI(K) - XPLUS(K) * CRANKP BI(K) = BI(K) + XPLUS(K) * CRANKP + XMINUS(K) * CRANKP ENDDO IF (NOCONV(I) .EQ. 1) THEN BI(1) = 1.0 + CRANKP * MBARKS(1,I) * (PBLSIG(I) - SIGMAF(1)) * & DTS * DSIGHI(1) + EDDYZ(I,1) * DSIGHI(1) * CRANKP * DTRAT ELSE BI(1) = 1.0 + EDDYZ(I,1) * DSIGHI(1) * CRANKP * DTRAT ENDIF DO K = 1,KL DO L = 1,NSPX DI(L,K) = 0.0 ENDDO ENDDO ! !** COMPUTE TENDENCY OF CBL CONCENTRATIONS - SEMI-IMPLICIT SOLUTION DO K = 2,KCBL(I) DO L = 1,NSPX DELC = DTS * (MBARKS(K,I) * VCI(L,I,1) - MDWN(K,I) * & VCI(L,I,K) + DSIGH(K+1) * DSIGHI(K) * & MDWN(K+1,I) * VCI(L,I,K+1)) DI(L,K) = VCI(L,I,K) + (1.0 - CRANKP) * DELC ENDDO ENDDO DO K = KCBL(I)+1, KL DO L = 1,NSPX DI(L,K) = VCI(L,I,K) ENDDO ENDDO DO K = 2,KL IF (K .EQ. KL) THEN DO L = 1,NSPX DI(L,K) = DI(L,K) - (1.0 - CRANKP) * XMINUS(K) * & (VCI(L,I,K) - VCI(L,I,K-1)) ENDDO ELSE DO L = 1,NSPX DI(L,K) = DI(L,K) + (1.0 - CRANKP) * XPLUS(K) * & (VCI(L,I,K+1) - VCI(L,I,K)) - & (1.0 - CRANKP) * XMINUS(K) * & (VCI(L,I,K) - VCI(L,I,K-1)) ENDDO ENDIF ENDDO IF (NOCONV(I) .EQ. 1) THEN DO L = 1,NSPX F1 = -G1000I * (MBARKS(1,I) * & (PBLSIG(I) - SIGMAF(1)) * VCI(L,I,1) - & MDWN(2,I) * VCI(L,I,2) * DSIGH(2)) DI(L,1) = VCI(L,I,1) - G1000 * (FS(L,I) - (1.0 - CRANKP) & * F1) * DSIGHI(1) * DTS ENDDO ELSE DO L = 1,NSPX DI(L,1) = VCI(L,I,1) - G1000 * FS(L,I) * DSIGHI(1) * DTS ENDDO ENDIF DO L = 1,NSPX DI(L,1) = DI(L,1) + (1.0 - CRANKP) * EDDYZ(I,1) * DSIGHI(1) & * DTRAT * (VCI(L,I,2) - VCI(L,I,1)) ENDDO IF ( NOCONV(I) .EQ. 1 ) THEN CALL MATRIX (AI, BI, CI, DI, EI, UI, KL, NSPX) ELSE CALL TRI (CI, BI, EI, DI, UI, KL, NSPX) END IF ! !-- COMPUTE NEW THETAV AND Q DO K = 1,KL DO L = 1,NSPX VCI(L,I,K) = UI(L,K) ENDDO ENDDO ! ENDDO ! END I LOOP ENDDO ! END SUB TIME LOOP !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! DO K = kts,KL DO I = its,ILX THETAX(I,K) = VCI(1,I,K) QVX(I,K) = VCI(2,I,K) UX(I,K) = VCI(3,I,K) VX(I,K) = VCI(4,I,K) ! IF( (I.EQ.61 .OR. I.EQ.15) .AND. JX.EQ.22) THEN ! PRINT *,'TEST -->',I,K,THETAX(I,K) ! ENDIF ENDDO ENDDO DO K = kts,KL DO I = its,ILX QCX(I,K) = VCI(5,I,K) QIX(I,K) = VCI(6,I,K) ENDDO ENDDO ! ENDIF END SUBROUTINE ACM !-------------------------------------------------------- SUBROUTINE MATRIX(A,B,C,D,E,X,KL,NSP) IMPLICIT NONE ! !-- Bordered band diagonal matrix solver for ACM2 !-- ACM2 Matrix is in this form: ! B1 E1 ! A2 B2 E2 ! A3 C3 B3 E3 ! A4 C4 B4 E4 ! A5 C5 B5 E5 ! A6 C6 B6 !--Upper Matrix is ! U11 U12 ! U22 U23 ! U33 U34 ! U44 U45 ! U55 U56 ! U66 !--Lower Matrix is: ! 1 ! L21 1 ! L31 L32 1 ! L41 L42 L43 1 ! L51 L52 L53 L54 1 ! L61 L62 L63 L64 L65 1 !--------------------------------------------------------- !...Arguments INTEGER, INTENT(IN) :: KL INTEGER, INTENT(IN) :: NSP REAL A(KL),B(KL),E(KL) REAL C(KL),D(NSP,KL),X(NSP,KL) !...Locals REAL Y(NSP,KL),AIJ,SUM REAL L(KL,KL),UII(KL),UIIP1(KL),RUII(KL) INTEGER I,J,V !-- Define Upper and Lower matrices L(1,1) = 1. ! U(1,1) = B(1) UII(1) = B(1) RUII(1) = 1./UII(1) DO I = 2, KL L(I,I) = 1. L(I,1) = A(I)/B(1) ! U(I-1,I) =E(I-1) UIIP1(I-1)=E(I-1) IF(I.GE.3) THEN DO J = 2,I-1 IF(I.EQ.J+1) THEN AIJ = C(I) ELSE AIJ = 0. ENDIF L(I,J) = (AIJ-L(I,J-1)*E(J-1))/ & (B(J)-L(J,J-1)*E(J-1)) ENDDO ENDIF ENDDO DO I = 2,KL ! U(I,I) = B(I)-L(I,I-1)*E(I-1) UII(I) = B(I)-L(I,I-1)*E(I-1) RUII(I) = 1./UII(I) ENDDO !-- Forward sub for Ly=d DO V= 1, NSP Y(V,1) = D(V,1) DO I=2,KL SUM = D(V,I) DO J=1,I-1 SUM = SUM - L(I,J)*Y(V,J) ENDDO Y(V,I) = SUM ENDDO ENDDO !-- Back sub for Ux=y DO V= 1, NSP ! X(V,KL) = Y(V,KL)/U(KL,KL) X(V,KL) = Y(V,KL)*RUII(KL) ENDDO DO I = KL-1,1,-1 DO V= 1, NSP ! X(V,I) = (Y(V,I)-U(I,I+1)*X(V,I+1))/U(I,I) X(V,I) = (Y(V,I)-UIIP1(I)*X(V,I+1))*RUII(I) ENDDO ENDDO END SUBROUTINE MATRIX !::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: SUBROUTINE TRI ( L, D, U, B, X,KL,NSP) !----------------------------------------------------------------------- ! FUNCTION: ! Solves tridiagonal system by Thomas algorithm. ! The associated tri-diagonal system is stored in 3 arrays ! D : diagonal ! L : sub-diagonal ! U : super-diagonal ! B : right hand side function ! X : return solution from tridiagonal solver ! [ D(1) U(1) 0 0 0 ... 0 ] ! [ L(2) D(2) U(2) 0 0 ... . ] ! [ 0 L(3) D(3) U(3) 0 ... . ] ! [ . . . . . ] X(i) = B(i) ! [ . . . . 0 ] ! [ . . . . ] ! [ 0 L(n) D(n) ] !----------------------------------------------------------------------- IMPLICIT NONE ! Arguments: INTEGER, INTENT(IN) :: KL INTEGER, INTENT(IN) :: NSP REAL L( KL ) ! subdiagonal REAL D(KL) ! diagonal REAL U( KL ) ! superdiagonal REAL B(NSP,KL ) ! R.H. side REAL X( NSP,KL ) ! solution ! Local Variables: REAL GAM( KL ) REAL BET INTEGER V, K ! Decomposition and forward substitution: BET = 1.0 / D( 1 ) DO V = 1, NSP X( V,1 ) = BET * B(V,1 ) ENDDO DO K = 2, KL GAM(K ) = BET * U( K-1 ) BET = 1.0 / ( D( K ) - L( K ) * GAM( K ) ) DO V = 1, NSP X( V, K ) = BET * ( B( V,K ) - L( K ) * X( V,K-1 ) ) ENDDO END DO ! Back-substitution: DO K = KL - 1, 1, -1 DO V = 1, NSP X( V,K ) = X( V,K ) - GAM( K+1 ) * X( V,K+1 ) END DO ENDDO END SUBROUTINE TRI END MODULE module_bl_acm