MODULE module_sf_mynn !------------------------------------------------------------------- !Modifications implemented by Joseph Olson NOAA/GSD/AMB - CU/CIRES !for WRFv3.4 and WRFv3.4.1: ! ! BOTH LAND AND WATER: !1) Calculation of stability parameter (z/L) taken from Li et al. (2010 BLM) ! for first iteration of first time step; afterwards, exact calculation. !2) Fixed isflux=0 option to turn off scalar fluxes, but keep momentum ! fluxes for idealized studies (credit: Anna Fitch). !3) Kinematic viscosity now varies with temperature !4) Uses Monin-Obukhov flux-profile relationships more consistent with ! those used in the MYNN PBL code. !5) Allows negative QFX, similar to MYJ scheme ! ! LAND only: !1) iz0tlnd option is now available with the following options: ! (default) =0: Zilitinkevich (1995) with Czil=0.1, ! =1: Czil_new (modified according to Chen & Zhang 2008) ! =2: Modified Yang et al (2002, 2008) - generalized for all landuse ! =3: constant zt = z0/7.4 (original form; Garratt 1992) !2) Relaxed u* minimum from 0.1 to 0.01 ! ! WATER only: !1) isftcflx option is now available with the following options: ! (default) =0: z0, zt, and zq from COARE3.0 (Fairall et al 2003) ! =1: z0 from Davis et al (2008), zt & zq from COARE3.0 ! =2: z0 from Davis et al (2008), zt & zq from Garratt (1992) ! =3: z0 from Taylor and Yelland (2004), zt and zq from COARE3.0 ! =4: z0 from Zilitinkevich (2001), zt & zq from COARE3.0 ! ! SNOW/ICE only: !1) Added Andreas (2002) snow/ice parameterization for thermal and ! moisture roughness to help reduce the cool/moist bias in the arctic ! region. ! !NOTE: This code was primarily tested in combination with the RUC LSM. ! Performance with the Noah (or other) LSM is relatively unknown. !------------------------------------------------------------------- USE module_model_constants, only: & &p1000mb, cp, xlv, ep_2 USE module_sf_sfclay, ONLY: sfclayinit USE module_bl_mynn, only: tv0, mym_condensation !------------------------------------------------------------------- IMPLICIT NONE !------------------------------------------------------------------- REAL, PARAMETER :: xlvcp=xlv/cp, ep_3=1.-ep_2 REAL, PARAMETER :: wmin=0.1 ! Minimum wind speed REAL, PARAMETER :: zm2h=7.4 ! = z_0m/z_0h REAL, PARAMETER :: charnock=0.016, bvisc=1.5e-5, z0hsea=5.e-5 REAL, PARAMETER :: VCONVC=1.0 REAL, PARAMETER :: SNOWZ0=0.012 REAL, DIMENSION(0:1000 ),SAVE :: PSIMTB,PSIHTB CONTAINS !------------------------------------------------------------------- SUBROUTINE mynn_sf_init_driver(allowed_to_read) LOGICAL, INTENT(in) :: allowed_to_read !Fill the PSIM and PSIH tables. The subroutine "sfclayinit" !can be found in module_sf_sfclay.F. This subroutine returns !the forms from Dyer and Hicks (1974). CALL sfclayinit(allowed_to_read) END SUBROUTINE mynn_sf_init_driver !------------------------------------------------------------------- SUBROUTINE SFCLAY_mynn(U3D,V3D,T3D,QV3D,P3D,dz8w, & CP,G,ROVCP,R,XLV,PSFC,CHS,CHS2,CQS2,CPM, & ZNT,UST,PBLH,MAVAIL,ZOL,MOL,REGIME,PSIM,PSIH, & XLAND,HFX,QFX,LH,TSK,FLHC,FLQC,QGH,QSFC,RMOL, & U10,V10,TH2,T2,Q2, & GZ1OZ0,WSPD,BR,ISFFLX,DX, & SVP1,SVP2,SVP3,SVPT0,EP1,EP2, & KARMAN,EOMEG,STBOLT, & itimestep,ch,th3d,pi3d,qc3d, & tsq,qsq,cov,qcg, & !JOE-add output ! z0zt_ratio,BulkRi,wstar,qstar,resist,logres, & ! Rreynolds,niters,psixrat,psitrat, & !JOE-end ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte, & ustm,ck,cka,cd,cda,isftcflx,iz0tlnd ) !------------------------------------------------------------------- IMPLICIT NONE !------------------------------------------------------------------- !-- U3D 3D u-velocity interpolated to theta points (m/s) !-- V3D 3D v-velocity interpolated to theta points (m/s) !-- T3D temperature (K) !-- QV3D 3D water vapor mixing ratio (Kg/Kg) !-- P3D 3D pressure (Pa) !-- dz8w dz between full levels (m) !-- CP heat capacity at constant pressure for dry air (J/kg/K) !-- G acceleration due to gravity (m/s^2) !-- ROVCP R/CP !-- R gas constant for dry air (J/kg/K) !-- XLV latent heat of vaporization for water (J/kg) !-- PSFC surface pressure (Pa) !-- ZNT roughness length (m) !-- UST u* in similarity theory (m/s) !-- USTM u* in similarity theory (m/s) w* added to WSPD. This is ! used to couple with TKE scheme but not in MYNN. ! (as of now, USTM = UST in this version) !-- PBLH PBL height from previous time (m) !-- MAVAIL surface moisture availability (between 0 and 1) !-- ZOL z/L height over Monin-Obukhov length !-- MOL T* (similarity theory) (K) !-- REGIME flag indicating PBL regime (stable, unstable, etc.) !-- PSIM similarity stability function for momentum !-- PSIH similarity stability function for heat !-- XLAND land mask (1 for land, 2 for water) !-- HFX upward heat flux at the surface (W/m^2) !-- QFX upward moisture flux at the surface (kg/m^2/s) !-- LH net upward latent heat flux at surface (W/m^2) !-- TSK surface temperature (K) !-- FLHC exchange coefficient for heat (W/m^2/K) !-- FLQC exchange coefficient for moisture (kg/m^2/s) !-- CHS heat/moisture exchange coefficient for LSM (m/s) !-- QGH lowest-level saturated mixing ratio !-- U10 diagnostic 10m u wind !-- V10 diagnostic 10m v wind !-- TH2 diagnostic 2m theta (K) !-- T2 diagnostic 2m temperature (K) !-- Q2 diagnostic 2m mixing ratio (kg/kg) !-- GZ1OZ0 log(z/z0) where z0 is roughness length !-- WSPD wind speed at lowest model level (m/s) !-- BR bulk Richardson number in surface layer !-- ISFFLX isfflx=1 for surface heat and moisture fluxes !-- DX horizontal grid size (m) !-- SVP1 constant for saturation vapor pressure (=0.6112 kPa) !-- SVP2 constant for saturation vapor pressure (=17.67 dimensionless) !-- SVP3 constant for saturation vapor pressure (=29.65 K) !-- SVPT0 constant for saturation vapor pressure (=273.15 K) !-- EP1 constant for virtual temperature (Rv/Rd - 1) (dimensionless) !-- EP2 constant for spec. hum. calc (Rd/Rv = 0.622) (dimensionless) !-- EP3 constant for spec. hum. calc (1 - Rd/Rv = 0.378 ) (dimensionless) !-- KARMAN Von Karman constant !-- EOMEG angular velocity of earth's rotation (rad/s) !-- STBOLT Stefan-Boltzmann constant (W/m^2/K^4) !-- ck enthalpy exchange coeff at 10 meters !-- cd momentum exchange coeff at 10 meters !-- cka enthalpy exchange coeff at the lowest model level !-- cda momentum exchange coeff at the lowest model level !-- isftcflx =0: z0, zt, and zq from COARE3.0 (Fairall et al 2003) ! (water =1: z0 from Davis et al (2008), zt & zq from COARE3.0 ! only) =2: z0 from Davis et al (2008), zt & zq from Garratt (1992) ! =3: z0 from Taylor and Yelland (2004), zt and zq from COARE3.0 ! =4: z0 from Zilitinkevich (2001), zt & zq from COARE3.0 !-- iz0tlnd =0: Zilitinkevich (1995) with Czil=0.1, ! (land =1: Czil_new (modified according to Chen & Zhang 2008) ! only) =2: Modified Yang et al (2002, 2008) - generalized for all landuse ! =3: constant zt = z0/7.4 (Garratt 1992) !-- 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 !-- its start index for i in tile !-- ite end index for i in tile !-- 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 !------------------------------------------------------------------- INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte ! INTEGER, INTENT(IN ) :: ISFFLX REAL, INTENT(IN ) :: SVP1,SVP2,SVP3,SVPT0 REAL, INTENT(IN ) :: EP1,EP2,KARMAN,EOMEG,STBOLT ! REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) , & INTENT(IN ) :: dz8w REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) , & INTENT(IN ) :: QV3D, & P3D, & T3D, & QC3D,& th3d,pi3d,tsq,qsq,cov INTEGER, INTENT(in) :: itimestep REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) ::& & qcg REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) ::& & ch REAL, DIMENSION( ims:ime, jms:jme ) , & INTENT(IN ) :: MAVAIL, & PBLH, & XLAND, & TSK REAL, DIMENSION( ims:ime, jms:jme ) , & INTENT(OUT ) :: U10, & V10, & TH2, & T2, & !JOE-use value from LSM Q2, & Q2 !JOE-moved down below QSFC ! REAL, DIMENSION( ims:ime, jms:jme ) , & INTENT(INOUT) :: REGIME, & HFX, & QFX, & LH, & MOL,RMOL,QSFC !m the following 5 are change to memory size ! REAL, DIMENSION( ims:ime, jms:jme ) , & INTENT(INOUT) :: GZ1OZ0,WSPD,BR, & PSIM,PSIH REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) , & INTENT(IN ) :: U3D, & V3D REAL, DIMENSION( ims:ime, jms:jme ) , & INTENT(IN ) :: PSFC REAL, DIMENSION( ims:ime, jms:jme ) , & INTENT(INOUT) :: ZNT, & ZOL, & UST, & CPM, & CHS2, & CQS2, & CHS REAL, DIMENSION( ims:ime, jms:jme ) , & INTENT(INOUT) :: FLHC,FLQC REAL, DIMENSION( ims:ime, jms:jme ) , & INTENT(INOUT) :: QGH !JOE-begin ! REAL, DIMENSION( ims:ime, jms:jme ) , & ! INTENT(OUT) :: z0zt_ratio, & ! BulkRi,wstar,qstar,resist,logres, & ! Rreynolds,niters,psixrat,psitrat !JOE-end REAL, INTENT(IN ) :: CP,G,ROVCP,R,XLV,DX REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ) , & INTENT(OUT) :: ck,cka,cd,cda,ustm INTEGER, OPTIONAL, INTENT(IN ) :: ISFTCFLX, IZ0TLND !----------- LOCAL VARS ----------------------------------- REAL, DIMENSION( its:ite ) :: U1D, & V1D, & QV1D, & P1D, & T1D,qc1d REAL, DIMENSION( its:ite ) :: dz8w1d REAL, DIMENSION( its:ite ) :: vt1,vq1 REAL, DIMENSION(kts:kts+1) :: thl, qw, vt, vq REAL :: ql INTEGER :: I,J,K !----------------------------------------------------------- DO J=jts,jte DO i=its,ite dz8w1d(I) = dz8w(i,kts,j) ENDDO DO i=its,ite U1D(i) =U3D(i,kts,j) V1D(i) =V3D(i,kts,j) QV1D(i)=QV3D(i,kts,j) QC1D(i)=QC3D(i,kts,j) P1D(i) =P3D(i,kts,j) T1D(i) =T3D(i,kts,j) ENDDO IF (itimestep==1) THEN DO i=its,ite vt1(i)=0. vq1(i)=0. UST(i,j)=MAX(0.025*SQRT(U1D(i)*U1D(i) + V1D(i)*V1D(i)),0.001) MOL(i,j)=0. ! Tstar !qstar(i,j)=0.0 ENDDO ELSE DO i=its,ite do k = kts,kts+1 ql = qc3d(i,k,j)/(1.+qc3d(i,k,j)) qw(k) = qv3d(i,k,j)/(1.+qv3d(i,k,j)) + ql thl(k) = th3d(i,k,j)-xlvcp*ql/pi3d(i,k,j) end do ! NOTE: The last grid number is kts+1 instead of kte. CALL mym_condensation (kts,kts+1, & & dz8w(i,kts:kts+1,j), & & thl(kts:kts+1), qw(kts:kts+1), & & p3d(i,kts:kts+1,j),& & pi3d(i,kts:kts+1,j), & & tsq(i,kts:kts+1,j), & & qsq(i,kts:kts+1,j), & & cov(i,kts:kts+1,j), & & vt(kts:kts+1), vq(kts:kts+1)) vt1(i) = vt(kts) vq1(i) = vq(kts) ENDDO ENDIF CALL SFCLAY1D_mynn(J,U1D,V1D,T1D,QV1D,P1D,dz8w1d, & CP,G,ROVCP,R,XLV,PSFC(ims,j),CHS(ims,j),CHS2(ims,j),& CQS2(ims,j),CPM(ims,j),PBLH(ims,j), RMOL(ims,j), & ZNT(ims,j),UST(ims,j),MAVAIL(ims,j),ZOL(ims,j), & MOL(ims,j),REGIME(ims,j),PSIM(ims,j),PSIH(ims,j), & XLAND(ims,j),HFX(ims,j),QFX(ims,j),TSK(ims,j), & U10(ims,j),V10(ims,j),TH2(ims,j),T2(ims,j), & Q2(ims,j),FLHC(ims,j),FLQC(ims,j),QGH(ims,j), & QSFC(ims,j),LH(ims,j), & GZ1OZ0(ims,j),WSPD(ims,j),BR(ims,j),ISFFLX,DX, & SVP1,SVP2,SVP3,SVPT0,EP1,EP2,KARMAN,EOMEG,STBOLT, & ch(ims,j),vt1,vq1,qc1d,qcg(ims,j),& itimestep,& !JOE-begin ! z0zt_ratio(ims,j),BulkRi(ims,j),wstar(ims,j),qstar(ims,j), & ! resist(ims,j),logres(ims,j),Rreynolds(ims,j),niters(ims,j), & ! psixrat(ims,j),psitrat(ims,j), & !JOE-end ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte & ,isftcflx,iz0tlnd, & USTM(ims,j),CK(ims,j),CKA(ims,j), & CD(ims,j),CDA(ims,j) & ) ENDDO END SUBROUTINE SFCLAY_MYNN !------------------------------------------------------------------- SUBROUTINE SFCLAY1D_mynn(J,UX,VX,T1D,QV1D,P1D,dz8w1d, & CP,G,ROVCP,R,XLV,PSFCPA,CHS,CHS2,CQS2,CPM,PBLH,RMOL, & ZNT,UST,MAVAIL,ZOL,MOL,REGIME,PSIM,PSIH, & XLAND,HFX,QFX,TSK, & U10,V10,TH2,T2,Q2,FLHC,FLQC,QGH, & QSFC,LH,GZ1OZ0,WSPD,BR,ISFFLX,DX, & SVP1,SVP2,SVP3,SVPT0,EP1,EP2, & KARMAN,EOMEG,STBOLT, & ch,vt1,vq1,qc1d,qcg, & itimestep, & !JOE-add ! zratio,BRi,wstar,qstar,resist,logres, & ! Rreynolds,niters,psixrat,psitrat, & !JOE-end ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte, & isftcflx, iz0tlnd, & ustm,ck,cka,cd,cda ) !------------------------------------------------------------------- IMPLICIT NONE !------------------------------------------------------------------- REAL, PARAMETER :: XKA=2.4E-5 !molecular diffusivity REAL, PARAMETER :: PRT=1. !prandlt number INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte, & J ! INTEGER, INTENT(in) :: itimestep INTEGER, INTENT(IN ) :: ISFFLX REAL, INTENT(IN ) :: SVP1,SVP2,SVP3,SVPT0 REAL, INTENT(IN ) :: EP1,EP2,KARMAN,EOMEG,STBOLT ! REAL, DIMENSION( ims:ime ) , & INTENT(IN ) :: MAVAIL, & PBLH, & XLAND, & TSK ! REAL, DIMENSION( ims:ime ) , & INTENT(IN ) :: PSFCPA REAL, DIMENSION( ims:ime ) , & INTENT(INOUT) :: REGIME, & HFX, & QFX, & MOL,RMOL !m the following 5 are changed to memory size--- ! REAL, DIMENSION( ims:ime ) , & INTENT(INOUT) :: GZ1OZ0,WSPD,BR, & PSIM,PSIH REAL, DIMENSION( ims:ime ) , & INTENT(INOUT) :: ZNT, & ZOL, & UST, & CPM, & CHS2, & CQS2, & CHS !JOE-add REAL, DIMENSION( its:ite ) :: zratio,BRi,wstar,qstar,& resist,logres,Rreynolds, & niters,psixrat,psitrat ! REAL, DIMENSION( ims:ime ) , & ! INTENT(OUT) :: zratio,BRi,wstar,qstar, & ! resist,logres,Rreynolds, & ! niters,psixrat,psitrat !JOE-end REAL, DIMENSION( ims:ime ) , & INTENT(INOUT) :: FLHC,FLQC REAL, DIMENSION( ims:ime ) , & INTENT(INOUT) :: QGH,QSFC REAL, DIMENSION( ims:ime ) , & INTENT(OUT) :: U10,V10, & !JOE-make qsfc inout (moved up) TH2,T2,Q2,QSFC,LH TH2,T2,Q2,LH REAL, INTENT(IN) :: CP,G,ROVCP,R,XLV,DX ! MODULE-LOCAL VARIABLES, DEFINED IN SUBROUTINE SFCLAY REAL, DIMENSION( its:ite ), INTENT(IN ) :: dz8w1d REAL, DIMENSION( its:ite ), INTENT(IN ) :: UX, & VX, & QV1D, & P1D, & T1D,qc1d REAL, DIMENSION( ims:ime ), INTENT(IN) :: qcg REAL, DIMENSION( ims:ime ), INTENT(INOUT) :: ch REAL, DIMENSION( its:ite ), INTENT(IN) :: vt1,vq1 REAL, OPTIONAL, DIMENSION( ims:ime ) , & INTENT(OUT) :: ck,cka,cd,cda,ustm INTEGER, OPTIONAL, INTENT(IN ) :: ISFTCFLX, IZ0TLND ! LOCAL VARS REAL, DIMENSION( its:ite ) :: z_t,z_q REAL :: thl1,sqv1,sqc1,exner1,sqvg,sqcg,vv,ww REAL, DIMENSION( its:ite ) :: ZA, & THVX,ZQKL, & THX,QX, & PSIH2, & PSIM2, & PSIH10, & PSIM10, & GZ2OZ0, & GZ10OZ0, & WSPDI ! REAL, DIMENSION( its:ite ) :: RHOX,GOVRTH ! REAL, DIMENSION( its:ite) :: SCR4 REAL, DIMENSION( its:ite ) :: THGB, PSFC, QSFCMR REAL, DIMENSION( its:ite ) :: GZ2OZt,GZ10OZt,GZ1OZt ! INTEGER :: N,I,K,KK,L,NZOL,NK,NZOL2,NZOL10, ITER INTEGER, PARAMETER :: ITMAX=5 REAL :: PL,THCON,TVCON,E1 REAL :: ZL,TSKV,DTHVDZ,DTHVM,VCONV,RZOL,RZOL2,RZOL10,ZOL2,ZOL10 REAL :: DTG,PSIX,DTTHX,DTHDZ,PSIX10,PSIT,PSIT2,PSIT10, & PSIQ,PSIQ2,PSIQ10 REAL :: FLUXC,VSGD real :: restar,VISC,psilim,DQG,OLDUST,OLDTST !------------------------------------------------------------------- !----CONVERT GROUND TEMPERATURE TO POTENTIAL TEMPERATURE: DO I=its,ite ! PSFC cmb (or kPa) PSFC(I)=PSFCPA(I)/1000. THGB(I)=TSK(I)*(100./PSFC(I))**ROVCP ENDDO ! ! SCR4(I,K) STORES EITHER TEMPERATURE OR VIRTUAL TEMPERATURE, ! DEPENDING ON IFDRY (CURRENTLY NOT USED, SO SCR4 == TVX). DO 30 I=its,ite ! PL cmb PL=P1D(I)/1000. THCON=(100./PL)**ROVCP THX(I)=T1D(I)*THCON SCR4(I)=T1D(I) THVX(I)=THX(I) QX(I)=0. 30 CONTINUE ! INITIALIZE SOME VARIABLES HERE: DO I=its,ite niters(I)=0. QGH(I)=0. CPM(I)=CP IF (itimestep .LE. 1) THEN qstar(I)=0.0 ENDIF ENDDO ! IF(IDRY.EQ.1)GOTO 80 DO 50 I=its,ite QX(I)=QV1D(I)/(1.+QV1D(I)) !CONVERT TO SPEC HUM TVCON=(1.+EP1*QX(I)) THVX(I)=THX(I)*TVCON SCR4(I)=T1D(I)*TVCON 50 CONTINUE ! DO 60 I=its,ite IF (TSK(I) .LT. 273.15) THEN !SATURATION VAPOR PRESSURE WRT ICE (SVP1=.6112; 10*mb) E1=SVP1*EXP(4648*(1./273.15 - 1./TSK(I)) - & 11.64*LOG(273.15/TSK(I)) + 0.02265*(273.15 - TSK(I))) ELSE !SATURATION VAPOR PRESSURE WRT WATER (Bolton 1980) E1=SVP1*EXP(SVP2*(TSK(I)-SVPT0)/(TSK(I)-SVP3)) ENDIF QSFC(I)=EP2*E1/(PSFC(I)-ep_3*E1) !specific humidity QSFCMR(I)=EP2*E1/(PSFC(I)-E1) !mixing ratio !FOR LAND POINTS, QSFC can come from previous time step (in LSM) !if(xland(i).gt.1.5 .or. QSFC(i).le.0.0) QSFC(I)=EP2*E1/(PSFC(I)-ep_3*E1) ! QGH CHANGED TO USE LOWEST-LEVEL AIR TEMP CONSISTENT WITH MYJSFC CHANGE ! Q2SAT = QGH IN LSM IF (TSK(I) .LT. 273.15) THEN !SATURATION VAPOR PRESSURE WRT ICE E1=SVP1*EXP(4648*(1./273.15 - 1./T1D(I)) - & 11.64*LOG(273.15/T1D(I)) + 0.02265*(273.15 - T1D(I))) ELSE !SATURATION VAPOR PRESSURE WRT WATER (Bolton 1980) E1=SVP1*EXP(SVP2*(T1D(I)-SVPT0)/(T1D(I)-SVP3)) ENDIF PL=P1D(I)/1000. QGH(I)=EP2*E1/(PL-ep_3*E1) !specific humidity !QGH(I)=EP2*E1/(PL-E1) !mixing ratio CPM(I)=CP*(1.+0.84*QX(I)/(1.-qx(i))) 60 CONTINUE 80 CONTINUE !-----COMPUTE THE HEIGHT OF FULL- AND HALF-SIGMA LEVELS ABOVE GROUND ! LEVEL, AND THE LAYER THICKNESSES. DO I=its,ite RHOX(I)=PSFC(I)*1000./(R*SCR4(I)) ZQKL(I)=dz8w1d(I) !first full-sigma level ZA(I)=0.5*ZQKL(I) !first half-sigma level GOVRTH(I)=G/THX(I) ENDDO DO I=its,ite WSPD(I)=SQRT(UX(I)*UX(I)+VX(I)*VX(I)) !account for partial condensation exner1=(p1d(i)/p1000mb)**ROVCP sqc1=qc1d(i)/(1.+qc1d(i)) !convert to spec hum. sqv1=qx(i) thl1=THX(I)-xlvcp/exner1*sqc1 sqvg=qsfc(i) sqcg=qcg(i)/(1.+qcg(i)) !convert to spec hum. vv = thl1-THGB(I) ww = mavail(i)*(sqv1-sqvg) + (sqc1-sqcg) TSKV=THGB(I)*(1.+EP1*QSFC(I)*MAVAIL(I)) DTHDZ=(THX(I)-THGB(I)) !DTHVDZ=(THVX(I)-TSKV) DTHVDZ= (vt1(i) + 1.0)*vv + (vq1(i) + tv0)*ww !-------------------------------------------------------- ! Calculate the convective velocity scale (WSTAR) and ! subgrid-scale velocity (VSGD) following Beljaars (1995, QJRMS) ! and Mahrt and Sun (1995, MWR), respectively !------------------------------------------------------- ! VCONV = 0.25*sqrt(g/tskv*pblh(i)*dthvm) ! Use Beljaars over land, old MM5 (Wyngaard) formula over water IF (xland(i).lt.1.5) then !LAND (xland == 1) fluxc = max(hfx(i)/rhox(i)/cp & + ep1*tskv*qfx(i)/rhox(i),0.) WSTAR(I) = vconvc*(g/TSK(i)*pblh(i)*fluxc)**.33 ELSE !WATER (xland == 2) IF(-DTHVDZ.GE.0)THEN DTHVM=-DTHVDZ ELSE DTHVM=0. ENDIF !JOE-the Wyngaard formula is ~3 times larger than the Beljaars !formula, so switch to Beljaars for water, but use VCONVC = 1.25, !as in the COARE3.0 bulk parameterizations. !WSTAR(I) = 2.*SQRT(DTHVM) fluxc = max(hfx(i)/rhox(i)/cp & + ep1*tskv*qfx(i)/rhox(i),0.) WSTAR(I) = 1.25*(g/TSK(i)*pblh(i)*fluxc)**.33 ENDIF !-------------------------------------------------------- ! Mahrt and Sun low-res correction ! (for 13 km ~ 0.37 m/s; for 3 km == 0 m/s) !-------------------------------------------------------- VSGD = 0.32 * (max(dx/5000.-1.,0.))**.33 WSPD(I)=SQRT(WSPD(I)*WSPD(I)+WSTAR(I)*WSTAR(I)+vsgd*vsgd) WSPD(I)=MAX(WSPD(I),wmin) !-------------------------------------------------------- ! CALCULATE THE BULK RICHARDSON NUMBER OF SURFACE LAYER, ! ACCORDING TO AKB(1976), EQ(12). !-------------------------------------------------------- BR(I)=GOVRTH(I)*ZA(I)*DTHVDZ/(WSPD(I)*WSPD(I)) !SET LIMITS ACCORDING TO Li et al. (2010) Boundary-Layer Meteorol (p.158) BR(I)=MAX(BR(I),-2.0) BR(I)=MIN(BR(I),1.0) BRi(I)=BR(I) !new variable for output - BR is not a "state" variable. ! IF PREVIOUSLY UNSTABLE, DO NOT LET INTO REGIMES 1 AND 2 (STABLE) !if (itimestep .GT. 1) THEN ! IF(MOL(I).LT.0.)BR(I)=MIN(BR(I),0.0) !ENDIF !IF(I .eq. 2)THEN ! write(*,1006)"BR:",BR(I)," fluxc:",fluxc," vt1:",vt1(i)," vq1:",vq1(i) ! write(*,1007)"XLAND:",XLAND(I)," WSPD:",WSPD(I)," DTHVDZ:",DTHVDZ," WSTAR:",WSTAR(I) !ENDIF ENDDO 1006 format(A,F7.3,A,f9.4,A,f9.5,A,f9.4) 1007 format(A,F2.0,A,f6.2,A,f7.3,A,f7.2) !-------------------------------------------------------------------- !-------------------------------------------------------------------- !--- BEGIN ITERATION LOOP (ITMAX=5); USUALLY CONVERGES IN TWO PASSES !-------------------------------------------------------------------- !-------------------------------------------------------------------- DO I=its,ite ITER = 1 DO WHILE (ITER .LE. ITMAX) niters(I)=ITER !COMPUTE KINEMATIC VISCOSITY VISC=(1.32+0.009*(T1D(I)-273.15))*1.E-5 IF((XLAND(I)-1.5).GE.0)THEN !-------------------------------------- ! WATER !-------------------------------------- !COMPUTE KINEMATIC VISCOSITY VISC=(1.32+0.009*(T1D(I)-273.15))*1.E-5 !-------------------------------------- !CALCULATE z0 (znt) !-------------------------------------- IF ( PRESENT(ISFTCFLX) ) THEN IF ( ISFTCFLX .EQ. 0 ) THEN !NAME OF SUBROUTINE IS MISLEADING - ACTUALLY VARIABLE CHARNOCK !PARAMETER FROM COARE3.0: CALL charnock_1955(ZNT(i),UST(i),WSPD(i),visc) ELSEIF ( ISFTCFLX .EQ. 1 .OR. ISFTCFLX .EQ. 2 ) THEN CALL davis_etal_2008(ZNT(i),UST(i)) ELSEIF ( ISFTCFLX .EQ. 3 ) THEN CALL Taylor_Yelland_2001(ZNT(i),UST(i),WSPD(i)) ELSEIF ( ISFTCFLX .EQ. 4 ) THEN CALL charnock_1955(ZNT(i),UST(i),WSPD(i),visc) ENDIF ELSE !DEFAULT TO COARE 3.0 CALL charnock_1955(ZNT(i),UST(i),WSPD(i),visc) ENDIF !COMPUTE ROUGHNESS REYNOLDS NUMBER (restar) USING NEW ZNT ! AHW: Garrattt formula: Calculate roughness Reynolds number ! Kinematic viscosity of air (linear approx to ! temp dependence at sea level) restar=MAX(ust(i)*ZNT(i)/visc, 0.1) !-------------------------------------- !CALCULATE z_t and z_q !-------------------------------------- IF ( PRESENT(ISFTCFLX) ) THEN IF ( ISFTCFLX .EQ. 0 ) THEN CALL fairall_2001(z_t(i),z_q(i),restar,UST(i),visc) ELSEIF ( ISFTCFLX .EQ. 1 ) THEN CALL fairall_2001(z_t(i),z_q(i),restar,UST(i),visc) ELSEIF ( ISFTCFLX .EQ. 2 ) THEN CALL garratt_1992(z_t(i),z_q(i),ZNT(i),restar,XLAND(I)) ELSEIF ( ISFTCFLX .EQ. 3 ) THEN CALL fairall_2001(z_t(i),z_q(i),restar,UST(i),visc) ELSEIF ( ISFTCFLX .EQ. 4 ) THEN CALL zilitinkevich_1995(ZNT(i),z_t(i),z_q(i),restar,& UST(I),KARMAN,XLAND(I),IZ0TLND) ENDIF ELSE !DEFAULT TO COARE 3.0 CALL fairall_2001(z_t(i),z_q(i),restar,UST(i),visc) ENDIF ELSE !-------------------------------------- ! LAND !-------------------------------------- !COMPUTE ROUGHNESS REYNOLDS NUMBER (restar) USING DEFAULT ZNT VISC=(1.32+0.009*(T1D(I)-273.15))*1.E-5 restar=MAX(ust(i)*ZNT(i)/visc, 0.1) !-------------------------------------- !GET z_t and z_q !-------------------------------------- !CHECK FOR SNOW/ICE POINTS OVER LAND IF ( ZNT(i) .LE. SNOWZ0 .AND. TSK(I) .LE. 273.15 ) THEN CALL Andreas_2002(ZNT(i),restar,z_t(i),z_q(i)) ELSE IF ( PRESENT(IZ0TLND) ) THEN IF ( IZ0TLND .LE. 1 ) THEN CALL zilitinkevich_1995(ZNT(i),z_t(i),z_q(i),restar,& UST(I),KARMAN,XLAND(I),IZ0TLND) ELSEIF ( IZ0TLND .EQ. 2 ) THEN CALL Yang_2008(ZNT(i),z_t(i),z_q(i),UST(i),MOL(I),& qstar(I),restar,visc,XLAND(I)) ELSEIF ( IZ0TLND .EQ. 3 ) THEN !Original MYNN in WRF-ARW used this form: CALL garratt_1992(z_t(i),z_q(i),ZNT(i),restar,XLAND(I)) ENDIF ELSE !DEFAULT TO ZILITINKEVICH WITH CZIL = 0.1 CALL zilitinkevich_1995(ZNT(i),z_t(i),z_q(i),restar,& UST(I),KARMAN,XLAND(I),0) ENDIF ENDIF ENDIF zratio(i)=znt(i)/z_t(i) Rreynolds(I)=restar !ADD RESISTANCE (SOMEWHAT FOLLOWING JIMENEZ ET AL. (2012)) TO PROTECT AGAINST !EXCESSIVE FLUXES WHEN USING A LOW FIRST MODEL LEVEL (ZA < 10 m). !Formerly: GZ1OZ0(I)= LOG(ZA(I)/ZNT(I)) GZ1OZ0(I)= LOG((ZA(I)+ZNT(I))/ZNT(I)) GZ1OZt(I)= LOG((ZA(I)+z_t(i))/z_t(i)) GZ2OZ0(I)= LOG((2.0+ZNT(I))/ZNT(I)) GZ2OZt(I)= LOG((2.0+z_t(i))/z_t(i)) GZ10OZ0(I)=LOG((10.+ZNT(I))/ZNT(I)) GZ10OZt(I)=LOG((10.+z_t(i))/z_t(i)) !-------------------------------------------------------------------- !--- DIAGNOSE BASIC PARAMETERS FOR THE APPROPRIATE STABILITY CLASS: ! ! THE STABILITY CLASSES ARE DETERMINED BY BR (BULK RICHARDSON NO.). ! ! CRITERIA FOR THE CLASSES ARE AS FOLLOWS: ! ! 1. BR .GE. 0.2; ! REPRESENTS NIGHTTIME STABLE CONDITIONS (REGIME=1), ! ! 2. BR .LT. 0.2 .AND. BR .GT. 0.0; ! REPRESENTS DAMPED MECHANICAL TURBULENT CONDITIONS ! (REGIME=2), ! ! 3. BR .EQ. 0.0 ! REPRESENTS FORCED CONVECTION CONDITIONS (REGIME=3), ! ! 4. BR .LT. 0.0 ! REPRESENTS FREE CONVECTION CONDITIONS (REGIME=4). ! !-------------------------------------------------------------------- IF (BR(I) .GT. 0.2) THEN !=================================================== !---CLASS 1; STABLE (NIGHTTIME) CONDITIONS: !=================================================== REGIME(I)=1. !COMPUTE z/L !CALL Li_etal_2010(ZOL(I),BR(I),ZA(I)/ZNT(I),zratio(I)) IF (ITER .EQ. 1 .AND. itimestep .LE. 1) THEN CALL Li_etal_2010(ZOL(I),BR(I),ZA(I)/ZNT(I),zratio(I)) ELSE ZOL(I)=ZA(I)*KARMAN*9.81*MOL(I)/(THX(I)*MAX(UST(I),0.001)**2) ZOL(I)=MAX(ZOL(I),0.0) ZOL(I)=MIN(ZOL(I),20.0) ENDIF !COMPUTE PSIM and PSIH IF((XLAND(I)-1.5).GE.0)THEN ! WATER !CALL PSI_Suselj_Sood_2010(PSIM(I),PSIH(I),ZOL(I)) !produces neg TKE !CALL PSI_Beljaars_Holtslag_1991(PSIM(I),PSIH(I),ZOL(I)) !CALL PSI_Businger_1971(PSIM(I),PSIH(I),ZOL(I)) CALL PSI_DyerHicks(PSIM(I),PSIH(I),ZOL(I), z_t(I), ZNT(I), ZA(I)) ELSE ! LAND !CALL PSI_Beljaars_Holtslag_1991(PSIM(I),PSIH(I),ZOL(I)) !CALL PSI_Businger_1971(PSIM(I),PSIH(I),ZOL(I)) !CALL PSI_Zilitinkevich_Esau_2007(PSIM(I),PSIH(I),ZOL(I)) CALL PSI_DyerHicks(PSIM(I),PSIH(I),ZOL(I), z_t(I), ZNT(I), ZA(I)) ENDIF !LOWER LIMIT ON PSI IN STABLE CONDITIONS psilim = -10. !JOE: this limit will be hit for z/L > 2, but ! appears to be necessary to control "runaway cooling" ! in the polar regions.. PSIM(I)=MAX(PSIM(I),psilim) PSIH(I)=MAX(PSIH(I),psilim) PSIM10(I)=10./ZA(I)*PSIM(I) PSIM10(I)=MAX(PSIM10(I),psilim) PSIH10(I)=PSIM10(I) PSIM2(I)=2./ZA(I)*PSIM(I) PSIM2(I)=MAX(PSIM2(I),psilim) PSIH2(I)=PSIM2(I) RMOL(I) = ZOL(I)/ZA(I) !1.0/L ELSEIF(BR(I) .GT. 0. .AND. BR(I) .LE. 0.2) THEN !======================================================== !---CLASS 2; DAMPED MECHANICAL TURBULENCE: !======================================================== REGIME(I)=2. !COMPUTE z/L !CALL Li_etal_2010(ZOL(I),BR(I),ZA(I)/ZNT(I),zratio(I)) IF (ITER .EQ. 1 .AND. itimestep .LE. 1) THEN CALL Li_etal_2010(ZOL(I),BR(I),ZA(I)/ZNT(I),zratio(I)) ELSE ZOL(I)=ZA(I)*KARMAN*9.81*MOL(I)/(THX(I)*MAX(UST(I),0.001)**2) ZOL(I)=MAX(ZOL(I),0.0) ZOL(I)=MIN(ZOL(I),5.0) ENDIF !COMPUTE PSIM and PSIH IF((XLAND(I)-1.5).GE.0)THEN ! WATER !CALL PSI_Suselj_Sood_2010(PSIM(I),PSIH(I),ZOL(I)) !CALL PSI_Beljaars_Holtslag_1991(PSIM(I),PSIH(I),ZOL(I)) !CALL PSI_Businger_1971(PSIM(I),PSIH(I),ZOL(I)) CALL PSI_DyerHicks(PSIM(I),PSIH(I),ZOL(I), z_t(I), ZNT(I), ZA(I)) ELSE ! LAND !CALL PSI_Beljaars_Holtslag_1991(PSIM(I),PSIH(I),ZOL(I)) !CALL PSI_Businger_1971(PSIM(I),PSIH(I),ZOL(I)) !CALL PSI_Zilitinkevich_Esau_2007(PSIM(I),PSIH(I),ZOL(I)) CALL PSI_DyerHicks(PSIM(I),PSIH(I),ZOL(I), z_t(I), ZNT(I), ZA(I)) ENDIF !LOWER LIMIT ON PSI IN WEAKLY STABLE CONDITIONS psilim = -10. !JOE: this limit is never hit in this regime. ! LOWER LIMIT ON PSI IN STABLE CONDITIONS PSIM(I)=MAX(PSIM(I),psilim) PSIH(I)=MAX(PSIH(I),psilim) PSIM10(I)=10./ZA(I)*PSIM(I) PSIM10(I)=MAX(PSIM10(I),psilim) PSIH10(I)=PSIM10(I) PSIM2(I)=2./ZA(I)*PSIM(I) PSIM2(I)=MAX(PSIM2(I),psilim) PSIH2(I)=PSIM2(I) ! 1.0 over Monin-Obukhov length RMOL(I)= ZOL(I)/ZA(I) ELSEIF(BR(I) .EQ. 0.) THEN !========================================================= !-----CLASS 3; FORCED CONVECTION/NEUTRAL: !========================================================= REGIME(I)=3. PSIM(I)=0.0 PSIH(I)=PSIM(I) PSIM10(I)=0. PSIH10(I)=PSIM10(I) PSIM2(I)=0. PSIH2(I)=PSIM2(I) !ZOL(I)=0. IF(UST(I) .LT. 0.01)THEN ZOL(I)=BR(I)*GZ1OZ0(I) ELSE ZOL(I)=KARMAN*GOVRTH(I)*ZA(I)*MOL(I)/(UST(I)*UST(I)) ENDIF RMOL(I) = ZOL(I)/ZA(I) ELSEIF(BR(I) .LT. 0.)THEN !========================================================== !-----CLASS 4; FREE CONVECTION: !========================================================== REGIME(I)=4. !COMPUTE z/L !CALL Li_etal_2010(ZOL(I),BR(I),ZA(I)/ZNT(I),zratio(I)) IF (ITER .EQ. 1 .AND. itimestep .LE. 1) THEN CALL Li_etal_2010(ZOL(I),BR(I),ZA(I)/ZNT(I),zratio(I)) ELSE ZOL(I)=ZA(I)*KARMAN*9.81*MOL(I)/(THX(I)*MAX(UST(I),0.001)**2) ZOL(I)=MAX(ZOL(I),-10.0) ZOL(I)=MIN(ZOL(I),0.0) ENDIF ZOL10=10./ZA(I)*ZOL(I) ZOL2=2./ZA(I)*ZOL(I) ZOL(I)=MIN(ZOL(I),0.) ZOL(I)=MAX(ZOL(I),-9.9999) ZOL10=MIN(ZOL10,0.) ZOL10=MAX(ZOL10,-9.9999) ZOL2=MIN(ZOL2,0.) ZOL2=MAX(ZOL2,-9.9999) NZOL=INT(-ZOL(I)*100.) RZOL=-ZOL(I)*100.-NZOL NZOL10=INT(-ZOL10*100.) RZOL10=-ZOL10*100.-NZOL10 NZOL2=INT(-ZOL2*100.) RZOL2=-ZOL2*100.-NZOL2 !COMPUTE PSIM and PSIH IF((XLAND(I)-1.5).GE.0)THEN ! WATER !CALL PSI_Suselj_Sood_2010(PSIM(I),PSIH(I),ZOL(I)) !CALL PSI_Hogstrom_1996(PSIM(I),PSIH(I),ZOL(I), z_t(I), ZNT(I), ZA(I)) !CALL PSI_Businger_1971(PSIM(I),PSIH(I),ZOL(I)) CALL PSI_DyerHicks(PSIM(I),PSIH(I),ZOL(I), z_t(I), ZNT(I), ZA(I)) ELSE ! LAND !CALL PSI_Hogstrom_1996(PSIM(I),PSIH(I),ZOL(I), z_t(I), ZNT(I), ZA(I)) !CALL PSI_Businger_1971(PSIM(I),PSIH(I),ZOL(I)) CALL PSI_DyerHicks(PSIM(I),PSIH(I),ZOL(I), z_t(I), ZNT(I), ZA(I)) ENDIF PSIM10(I)=PSIMTB(NZOL10)+RZOL10*(PSIMTB(NZOL10+1)-PSIMTB(NZOL10)) PSIH10(I)=PSIHTB(NZOL10)+RZOL10*(PSIHTB(NZOL10+1)-PSIHTB(NZOL10)) PSIM2(I)=PSIMTB(NZOL2)+RZOL2*(PSIMTB(NZOL2+1)-PSIMTB(NZOL2)) PSIH2(I)=PSIHTB(NZOL2)+RZOL2*(PSIHTB(NZOL2+1)-PSIHTB(NZOL2)) !---LIMIT PSIH AND PSIM IN THE CASE OF THIN LAYERS AND !---HIGH ROUGHNESS. THIS PREVENTS DENOMINATOR IN FLUXES !---FROM GETTING TOO SMALL PSIH(I)=MIN(PSIH(I),0.9*GZ1OZ0(I)) PSIM(I)=MIN(PSIM(I),0.9*GZ1OZ0(I)) PSIH2(I)=MIN(PSIH2(I),0.9*GZ2OZ0(I)) PSIM10(I)=MIN(PSIM10(I),0.9*GZ10OZ0(I)) RMOL(I) = ZOL(I)/ZA(I) ENDIF !------------------------------------------------------------ !-----COMPUTE THE FRICTIONAL VELOCITY: !------------------------------------------------------------ ! ZA(1982) EQS(2.60),(2.61). GZ1OZ0(I) =LOG((ZA(I)+ZNT(I))/ZNT(I)) GZ10OZ0(I)=LOG((10.+ZNT(I))/ZNT(I)) PSIX=GZ1OZ0(I)-PSIM(I) PSIX10=GZ10OZ0(I)-PSIM10(I) ! TO PREVENT OSCILLATIONS AVERAGE WITH OLD VALUE OLDUST = UST(I) UST(I)=0.5*UST(I)+0.5*KARMAN*WSPD(I)/PSIX !NON-AVERAGED: UST(I)=KARMAN*WSPD(I)/PSIX ! Compute u* without vconv for use in HFX calc when isftcflx > 0 WSPDI(I)=SQRT(UX(I)*UX(I)+VX(I)*VX(I)) IF ( PRESENT(USTM) ) THEN USTM(I)=0.5*USTM(I)+0.5*KARMAN*WSPDI(I)/PSIX ENDIF IF ((XLAND(I)-1.5).LT.0.) THEN !LAND !JOE: UST(I)=MAX(UST(I),0.1) UST(I)=MAX(UST(I),0.01) !Relaxing this limit !Keep ustm = ust over land. USTM(I)=UST(I) ENDIF !------------------------------------------------------------ !-----COMPUTE THE THERMAL AND MOISTURE RESISTANCE (PSIQ AND PSIT): !------------------------------------------------------------ ! LOWER LIMIT ADDED TO PREVENT LARGE FLHC IN SOIL MODEL ! ACTIVATES IN UNSTABLE CONDITIONS WITH THIN LAYERS OR HIGH Z0 GZ1OZt(I)= LOG((ZA(I)+z_t(i))/z_t(i)) GZ2OZt(I)= LOG((2.0+z_t(i))/z_t(i)) !PSIT=MAX(GZ1OZ0(I)-PSIH(I),2.) PSIT=MAX(LOG((ZA(I)+z_t(i))/z_t(i))-PSIH(I) ,2.0) PSIT2=MAX(LOG((2.0+z_t(i))/z_t(i))-PSIH2(I) ,2.0) resist(I)=PSIT logres(I)=GZ1OZt(I) PSIQ=MAX(LOG((za(i)+z_q(i))/z_q(I))-PSIH(I) ,2.0) PSIQ2=MAX(LOG((2.0+z_q(i))/z_q(I))-PSIH2(I) ,2.0) !CARLSON AND BOLAND (1978): IF((XLAND(I)-1.5).GE.0)THEN ZL=ZNT(I) ELSE ZL=0.01 !PSIQ =MAX(LOG(KARMAN*UST(I)*ZA(I)/XKA + ZA(I)/ZL)-PSIH(I),2.0) !PSIQ2=MAX(LOG(KARMAN*UST(I)*2./XKA + 2./ZL)-PSIH2(I) ,2.0) ENDIF !---------------------------------------------------- !COMPUTE THE TEMPERATURE SCALE (or FRICTION TEMPERATURE, T*) !---------------------------------------------------- DTG=THX(I)-THGB(I) OLDTST=MOL(I) MOL(I)=KARMAN*DTG/PSIT/PRT !t_star(I) = -HFX(I)/(UST(I)*CPM(I)*RHOX(I)) !t_star(I) = MOL(I) !---------------------------------------------------- !COMPUTE THE MOISTURE SCALE (or q*) DQG=(QX(i)-qsfc(i))*1000. !(kg/kg -> g/kg) qstar(I)=KARMAN*DQG/PSIQ/PRT !----------------------------------------------------- !COMPUTE DIAGNOSTICS !----------------------------------------------------- !COMPUTE 10 M WNDS !----------------------------------------------------- ! If the lowest model level is close to 10-m, use it ! instead of the flux-based diagnostic formula. if (ZA(i) .gt. 7.0 .and. ZA(i) .lt. 13.0) then U10(I)=UX(I) V10(I)=VX(I) else U10(I)=UX(I)*PSIX10/PSIX V10(I)=VX(I)*PSIX10/PSIX endif psixrat(I)=PSIX10/PSIX psitrat(I)=PSIT2/PSIT !----------------------------------------------------- !COMPUTE 2m T, TH, AND Q !THESE WILL BE OVERWRITTEN FOR LAND POINTS IN THE LSM !----------------------------------------------------- TH2(I)=THGB(I)+DTG*PSIT2/PSIT !*** BE CERTAIN THAT THE 2-M THETA IS BRACKETED BY !*** THE VALUES AT THE SURFACE AND LOWEST MODEL LEVEL. ! !IF (THX(I)>THGB(I) .AND. (TH2(I)THX(I)) .OR. & ! THX(I)THGB(I) .OR. TH2(I)0.2:",zL zL = MIN(zL,20.) !LIMITS ACCORDING TO Li et al (2010), THIER !FIGUE 1C. zL = MAX(zL,1.) ENDIF return END SUBROUTINE Li_etal_2010 !-------------------------------------------------------------------- END MODULE module_sf_mynn