diff --git a/Registry/Registry.EM_COMMON b/Registry/Registry.EM_COMMON index d0065b6063..b960b676ed 100644 --- a/Registry/Registry.EM_COMMON +++ b/Registry/Registry.EM_COMMON @@ -1050,6 +1050,10 @@ i1 real PSIM ij misc 1 - - "PSIM i1 real PSIH ij misc 1 - - "PSIH" "SIMILARITY FUNCTION FOR HEAT" "" state real FM ij misc 1 - - "FM" "INTEGRATED FUNCTION FOR MOMENTUM" "" state real FH ij misc 1 - - "FH" "INTEGRATED FUNCTION FOR HEAT" "" +#$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$-Inserted by PN, MS, PS, and SKM-$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ +state real FM_hat ij misc 1 - - "FM_HAT" " " "" +state real FH_hat ij misc 1 - - "FH_HAT" "" "" +#$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ state real WSPD ij misc 1 - - "WSPD" "Wind Speed At Lowest Model Level (may contain vconv)" "m s-1" i1 real GZ1OZ0 ij misc 1 - - "GZ1OZ0" "LOG OF Z1 over Z0" "" state real BR ij misc 1 - - "BR" "Bulk Richardson" "" @@ -2656,6 +2660,10 @@ rconfig integer isftcflx namelist,physics 1 0 rconfig integer iz0tlnd namelist,physics 1 0 h "iz0tlnd" "switch to control land thermal roughness length" "" rconfig logical kim_tofd namelist,physics max_domains .false. rh "kim_tofd" "turbulent form drag (tofd) option for kim gravity wave drag" "" rconfig real tofd_factor namelist,physics max_domains .003 rh "tofd_factor" "factor in kim tofd scheme" "" +#!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$-Inserted by PN, MS, PS, and SKM-$$$$$$$$$$$$$$$$$$$$$$$$$$$ +rconfig integer psimh_stab_opt namelist,physics 1 0 - "psimh_stab_op t" "0: CB05, 1: GR07, 2: SR20, 3: GRY20, 4: H90, 5: BH91, 6: D91" "" +rconfig integer psimh_unstab_opt namelist,physics 1 0 - "psimh_unstab_o pt" "0: F96, 1: BD71, 2: CL73, 3: KY90 functions" "" +#$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ rconfig real shadlen namelist,physics 1 25000. - "shadow_length" "maximum length of orographic shadow" "m" rconfig integer slope_rad namelist,physics max_domains 0 - "slope_rad" "1: use slope-dependent radiation, 0:not" "" rconfig integer topo_shading namelist,physics max_domains 0 - "topo_shading" "1: apply topographic shading to radiation, 0:not" "" @@ -3195,6 +3203,9 @@ package gfssfcscheme sf_sfclay_physics==3 - - package qnsesfcscheme sf_sfclay_physics==4 - state:u10e,v10e package mynnsfcscheme sf_sfclay_physics==5 - - package pxsfcscheme sf_sfclay_physics==7 - - +#!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$-Inserted by PN, MS, PS, and SKM-$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ +package sfclayrev_IITD_scheme sf_sfclay_physics==8 - - +#$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ package temfsfcscheme sf_sfclay_physics==10 - state:wm_temf package idealscmsfcscheme sf_sfclay_physics==89 - - package sfclayscheme sf_sfclay_physics==91 - - diff --git a/dyn_em/module_first_rk_step_part1.F b/dyn_em/module_first_rk_step_part1.F index 2f1809c815..01fc4863fe 100644 --- a/dyn_em/module_first_rk_step_part1.F +++ b/dyn_em/module_first_rk_step_part1.F @@ -687,6 +687,10 @@ SUBROUTINE first_rk_step_part1 ( grid , config_flags & & ,USTM=grid%ustm ,CK=grid%ck ,CKA=grid%cka & & ,CD=grid%cd ,CDA=grid%cda & & ,ISFTCFLX=config_flags%isftcflx, IZ0TLND=config_flags%iz0tlnd & +!$$$$$$$$$$$$$$$$$$$$-Inserted by PN, MS, PS, and SKM-$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ + & ,PSIMH_STAB_OPT=config_flags%psimh_stab_opt & + & ,PSIMH_UNSTAB_OPT=config_flags%psimh_unstab_opt & +!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ & ,SF_OCEAN_PHYSICS=config_flags%sf_ocean_physics & & ,OML_HML0=config_flags%oml_hml0 ,OML_GAMMA=config_flags%oml_gamma & & ,TML=grid%tml, T0ML=grid%t0ml, HML=grid%hml, H0ML=grid%h0ml & @@ -1128,6 +1132,9 @@ SUBROUTINE first_rk_step_part1 ( grid , config_flags & & ,spp_lsm=config_flags%spp_lsm,pattern_spp_lsm=grid%pattern_spp_lsm & !SPP & ,field_sf=grid%field_sf & !SPP & ,spp_pbl=config_flags%spp_pbl,pattern_spp_pbl=grid%pattern_spp_pbl & !SPP +!$$$$$$$$$$$$$$$$$$$$$$$$$$-Inserted by PN, MS, PS, and SKM-$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ + & ,FM_hat=grid%fm_hat, FH_hat=grid%fh_hat & +!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ ! WRF-Solar EPS: Noah LSM & ,multi_perturb=config_flags% multi_perturb & & ,pert_noah=config_flags%pert_noah & diff --git a/dyn_em/start_em.F b/dyn_em/start_em.F index 524f338dd3..364c266b3f 100644 --- a/dyn_em/start_em.F +++ b/dyn_em/start_em.F @@ -1296,6 +1296,10 @@ SUBROUTINE start_domain_em ( grid, allowed_to_read & ,grid%NC_CU, grid%NI_CU, grid%NR_CU, grid%NS_CU,grid%CCN_CU & ,grid%alevsiz_cu,grid%num_months,grid%no_src_types_cu,grid%aeromcu & ,grid%aeropcu,grid%EFCG,grid%EFCS,grid%EFIG,grid%EFIS,grid%EFSG,grid%EFSS & +!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$-Inserted by PN, MS, PS, and SKM-$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ + ,grid%psimh_stab_opt & + ,grid%psimh_unstab_opt & +!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ #ifdef CDXWRF ,grid%fracsoil10cm & ,model_config_rec%londbg, model_config_rec%latdbg & diff --git a/main/depend.common b/main/depend.common index 3b3dc21149..cb8fa30fff 100644 --- a/main/depend.common +++ b/main/depend.common @@ -2407,6 +2407,7 @@ module_physics_init.o: \ module_ra_flg.o \ module_sf_sfclay.o \ physics_mmm/sf_sfclayrev.o \ + module_sf_sfclayrev_IITD.o \ module_sf_slab.o \ module_sf_myjsfc.o \ module_sf_mynnsfc_driver.o \ @@ -2696,6 +2697,7 @@ module_surface_driver.o: \ module_sf_noahlsm.o \ module_sf_sfclay.o \ module_sf_sfclayrev.o \ + module_sf_sfclayrev_IITD.o \ module_sf_slab.o \ module_sf_myjsfc.o \ module_sf_qnsesfc.o \ diff --git a/phys/Makefile b/phys/Makefile index 93cc5fcc37..626f5512d7 100644 --- a/phys/Makefile +++ b/phys/Makefile @@ -148,6 +148,7 @@ MODULES = \ module_ra_farms.o \ module_sf_sfclay.o \ module_sf_sfclayrev.o \ + module_sf_sfclayrev_IITD.o \ module_sf_gfs.o \ module_sf_exchcoef.o \ module_sf_gfdl.o \ diff --git a/phys/module_fddagd_driver.F b/phys/module_fddagd_driver.F index ba5fcdfda4..a7d6e73359 100644 --- a/phys/module_fddagd_driver.F +++ b/phys/module_fddagd_driver.F @@ -508,7 +508,10 @@ SUBROUTINE fddagd_driver(itimestep,dt,xtime, & IF( config_flags%sf_sfclay_physics /= sfclayscheme & .AND. config_flags%sf_sfclay_physics /= mynnsfcscheme & .AND. config_flags%sf_sfclay_physics /= pxsfcscheme & - .AND. config_flags%sf_sfclay_physics /= sfclayrevscheme ) THEN + .AND. config_flags%sf_sfclay_physics /= sfclayrevscheme & +!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$-Inserted by PN, MS, PS, and SKM-$$$$$$$$$$$$$$$$$$$$$$ + .AND. config_flags%sf_sfclay_physics /= sfclayrev_IITD_scheme ) THEN +!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ DO j=MAX(j_start(ij)-1,jds),j_end(ij) DO i=MAX(i_start(ij)-1,ids),i_end(ij) IF( pblh(i,j) > z_at_w(i,2,j)-ht(i,j) ) THEN diff --git a/phys/module_fddaobs_driver.F b/phys/module_fddaobs_driver.F index 8a91c8dc25..9e0b0b5414 100644 --- a/phys/module_fddaobs_driver.F +++ b/phys/module_fddaobs_driver.F @@ -235,7 +235,10 @@ SUBROUTINE fddaobs_driver( inest, domid, parid, restart, & IF( config_flags%sf_sfclay_physics /= sfclayscheme & .AND. config_flags%sf_sfclay_physics /= mynnsfcscheme & .AND. config_flags%sf_sfclay_physics /= pxsfcscheme & - .AND. config_flags%sf_sfclay_physics /= sfclayrevscheme ) THEN + .AND. config_flags%sf_sfclay_physics /= sfclayrevscheme & +!$$$$$$$$$$$$$$$$$$$$$-Inserted by PN, MS, PS, and SKM-$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ + .AND. config_flags%sf_sfclay_physics /= sfclayrev_IITD_scheme ) THEN +!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ DO j = jts, jte DO i = its, ite IF( pblh(i,j) > z_at_w(i,2,j)-ht(i,j) ) THEN diff --git a/phys/module_physics_init.F b/phys/module_physics_init.F index 9e604a4f37..23a790950b 100644 --- a/phys/module_physics_init.F +++ b/phys/module_physics_init.F @@ -280,6 +280,9 @@ SUBROUTINE phy_init ( id, config_flags, DT, restart, zfull, zhalf, & ,alevsiz_cu,num_months,no_src_types_cu,aeromcu,aeropcu & ! PSH/TWG 06/10/16 ,EFCG,EFCS,EFIG,EFIS,EFSG,EFSS & ! TWG #endif +!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$-Inserted by PN, MS, PS, and SKM-$$$$$ + , psimh_stab_opt,psimh_unstab_opt & +!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ #ifdef CDXWRF ,fracsoil10cm & ,londbg, latdbg & @@ -442,6 +445,9 @@ SUBROUTINE phy_init ( id, config_flags, DT, restart, zfull, zhalf, & INTEGER , INTENT(OUT) :: STEPRA, STEPBL, STEPCU INTEGER , INTENT(IN) :: JULYR, JULDAY +!$$$$$$$$$$$$$$$$$$$$$$$$$$-Inserted by PN, MS, PS, and SKM-$$$$$$$$$$$$$$$ + INTEGER , INTENT(IN) :: psimh_stab_opt, psimh_unstab_opt +!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ ! cps @@ -1664,6 +1670,9 @@ SUBROUTINE phy_init ( id, config_flags, DT, restart, zfull, zhalf, & ALBEDOsubgrid,LHsubgrid,HFXsubgrid,LWUPsubgrid, & Q2subgrid,SABVsubgrid,SABGsubgrid,NRAsubgrid, & SWUPsubgrid,lhsoi,lhveg,lhtran & +!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$-Inserted by PN, MS, PS, and SKM-$$$$$$$$$$$$$$$$$$$$ + ,psimh_stab_opt,psimh_unstab_opt & +!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ #if (WRF_USE_CLM == 1) && ( WRF_CHEM == 1) ,megan_specifier, megan_factors_file & ,megan_mapped_emisfctrs & @@ -2733,6 +2742,9 @@ SUBROUTINE bl_init(STEPBL,BLDT,DT,RUBLTEN,RVBLTEN,RTHBLTEN, & ALBEDOsubgrid,LHsubgrid,HFXsubgrid,LWUPsubgrid, & Q2subgrid,SABVsubgrid,SABGsubgrid,NRAsubgrid, & SWUPsubgrid,lhsoi,lhveg,lhtran & +!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$-Inserted by PN, MS, PS, and SKM-$$$$$$$$$$$$$$$$$$ $$$$$ + ,psimh_stab_opt, psimh_unstab_opt & +!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ $$$$$ #if (WRF_USE_CLM == 1) && ( WRF_CHEM == 1 ) ,megan_specifier, megan_factors_file & ,megan_mapped_emisfctrs & @@ -2761,6 +2773,9 @@ SUBROUTINE bl_init(STEPBL,BLDT,DT,RUBLTEN,RVBLTEN,RTHBLTEN, & !-------------------------------------------------------------------- USE module_sf_sfclay USE sf_sfclayrev +!$$$$$$$$$$$$$$$$$$$$$$$$-Inserted by PN, MS, PS, and SKM-$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ + USE module_sf_sfclayrev_IITD +!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ USE module_sf_slab USE module_sf_pxsfclay USE module_bl_ysu @@ -3185,7 +3200,9 @@ SUBROUTINE bl_init(STEPBL,BLDT,DT,RUBLTEN,RVBLTEN,RTHBLTEN, & TSK_SAVE LOGICAL, INTENT(IN) :: start_of_simulation INTEGER :: i,j,patch - +!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$-Inserted by PN, MS, PS, and SKM-$$$$$$$$$ + INTEGER, INTENT(IN) :: psimh_stab_opt, psimh_unstab_opt +!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ ! CLM varaibles INTEGER , INTENT(IN) :: maxpatch REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: HT @@ -3328,6 +3345,11 @@ SUBROUTINE bl_init(STEPBL,BLDT,DT,RUBLTEN,RVBLTEN,RTHBLTEN, & xland,LakeModel,lake_depth,lakemask ) END IF isfc = 1 +!$$$$$$$$$$$$$$$$$$$$$$$-Inserted PN, MS, PS, and SKM-$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ + CASE (SFCLAYREV_IITD_SCHEME) + CALL sfclayrevinit_IITD(psimh_stab_opt,psimh_unstab_opt) + isfc = 1 +!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ CASE (PXSFCSCHEME) CALL pxsfclayinit( allowed_to_read ) isfc = 7 diff --git a/phys/module_sf_sfclayrev_IITD.F b/phys/module_sf_sfclayrev_IITD.F new file mode 100755 index 0000000000..eca6d1baa2 --- /dev/null +++ b/phys/module_sf_sfclayrev_IITD.F @@ -0,0 +1,2350 @@ +!WRF:MODEL_LAYER:PHYSICS +! +MODULE module_sf_sfclayrev_IITD + + REAL , PARAMETER :: VCONVC=1. + REAL , PARAMETER :: CZO=0.0185 + REAL , PARAMETER :: OZO=1.59E-5 + + REAL, DIMENSION(0:1000 ),SAVE :: psim_stab,psim_unstab,psih_stab,psih_unstab + +CONTAINS + +!------------------------------------------------------------------- +!$$$$$$$$$$$$$$$$$$$-Inserted by PN, MS, PS, and SKM-$$$$$$$$$$$$$$$ + SUBROUTINE SFCLAYREV_IITD(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, & +!$$$$$$$$$$$$$$$$$$$-Inserted by PN, MS, PS, and SKM-$$$$$$$$$$$$$$$ + FM,FH,FM_hat,FH_hat, & +!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ + 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, & + P1000mb, & + ids,ide, jds,jde, kds,kde, & + ims,ime, jms,jme, kms,kme, & + its,ite, jts,jte, kts,kte, & +!$$$$$$$$$$$$$$$$$$-Inserted by PN, MS, PS, and SKM-$$$$$$$$$$$$$$$$ + ustm,ck,cka,cd,cda,isftcflx,iz0tlnd, & + psimh_stab_opt,psimh_unstab_opt, scm_force_flux) +!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ +!------------------------------------------------------------------- + IMPLICIT NONE +!------------------------------------------------------------------- +! Changes in V3.7 over water surfaces: +! 1. for ZNT/Cd, replacing constant OZO with 0.11*1.5E-5/UST(I) +! the COARE 3.5 (Edson et al. 2013) formulation is also available +! 2. for VCONV, reducing magnitude by half +! 3. for Ck, replacing Carlson-Boland with COARE 3 +!------------------------------------------------------------------- +!-- 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) without vconv correction +! used to couple with TKE scheme +!-- 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 +!-- FM integrated stability function for momentum +!-- FH integrated 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 +!-- QSFC ground 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 (kPa) +!-- SVP2 constant for saturation vapor pressure (dimensionless) +!-- SVP3 constant for saturation vapor pressure (K) +!-- SVPT0 constant for saturation vapor pressure (K) +!-- EP1 constant for virtual temperature (R_v/R_d - 1) (dimensionless) +!-- EP2 constant for specific humidity calculation +! (R_d/R_v) (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, (Charnock and Carlson-Boland); =1, AHW Ck, Cd, =2 Garratt +!-- iz0tlnd =0 Carlson-Boland, =1 Czil_new +!-- 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, INTENT(IN ) :: P1000mb +! + REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) , & + INTENT(IN ) :: dz8w + + REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) , & + INTENT(IN ) :: QV3D, & + P3D, & + T3D + + REAL, DIMENSION( ims:ime, jms:jme ) , & + INTENT(IN ) :: MAVAIL, & + PBLH, & + XLAND, & + TSK + REAL, DIMENSION( ims:ime, jms:jme ) , & + INTENT(OUT ) :: U10, & + V10, & + TH2, & + T2, & + Q2, & + QSFC + +! + REAL, DIMENSION( ims:ime, jms:jme ) , & + INTENT(INOUT) :: REGIME, & + HFX, & + QFX, & + LH, & + MOL,RMOL +!m the following 5 are change to memory size +! + REAL, DIMENSION( ims:ime, jms:jme ) , & + INTENT(INOUT) :: GZ1OZ0,WSPD,BR, & + PSIM,PSIH,FM,FH, & +!$$$$$$$$$$$$$$$$$$-Inserted by PN, MS, PS, and SKM-$$$$$$$$$$$$$$$$ + FM_hat, FH_hat +!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ + 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 + + REAL, INTENT(IN ) :: CP,G,ROVCP,R,XLV,DX + + REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ) , & + INTENT(OUT) :: ck,cka,cd,cda + + REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ) , & + INTENT(INOUT) :: USTM + + INTEGER, OPTIONAL, INTENT(IN ) :: ISFTCFLX, IZ0TLND + INTEGER, OPTIONAL, INTENT(IN ) :: SCM_FORCE_FLUX +!$$$$$$$$$$$$$$$$$-Inserted by PN, MS, PS, and SKM-$$$$$$$$$$$$$$$$$ + INTEGER, INTENT(IN) :: psimh_unstab_opt + INTEGER, INTENT(IN) :: psimh_stab_opt +!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ +! LOCAL VARS + + REAL, DIMENSION( its:ite ) :: U1D, & + V1D, & + QV1D, & + P1D, & + T1D + + REAL, DIMENSION( its:ite ) :: dz8w1d + + INTEGER :: I,J + + DO J=jts,jte + DO i=its,ite + dz8w1d(I) = dz8w(i,1,j) + ENDDO + + DO i=its,ite + U1D(i) =U3D(i,1,j) + V1D(i) =V3D(i,1,j) + QV1D(i)=QV3D(i,1,j) + P1D(i) =P3D(i,1,j) + T1D(i) =T3D(i,1,j) + ENDDO + + ! Sending array starting locations of optional variables may cause + ! troubles, so we explicitly change the call. + + CALL SFCLAYREV1D(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), & +!$$$$$$$$$$$$$$$$$-Inserted by PN, MS, PS, and SKM-$$$$$$$$$$$$$$$$$ + FM(ims,j),FH(ims,j),FM_hat(ims,j),FH_hat(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, & +!$$$$$$$$$$$$$$$$$-Inserted by PN, MS, PS, and SKM-$$$$$$$$$$$$$$$$$ + P1000mb, psimh_stab_opt, psimh_unstab_opt, & +!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ + ids,ide, jds,jde, kds,kde, & + ims,ime, jms,jme, kms,kme, & + its,ite, jts,jte, kts,kte & +#if ( EM_CORE == 1 ) + ,isftcflx,iz0tlnd,scm_force_flux, & + USTM(ims,j),CK(ims,j),CKA(ims,j), & + CD(ims,j),CDA(ims,j) & +#endif + ) + ENDDO + +!$$$$$$$$$$$$$$$$$$-Inserted by PN, MS, PS, and SKM-$$$$$$$$$$$$$$$$$$$ + END SUBROUTINE SFCLAYREV_IITD +!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ + +!------------------------------------------------------------------- + SUBROUTINE SFCLAYREV1D(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,FM,FH,& +!$$$$$$$$$$$$$$$$$$-Inserted by PN, MS, PS, and SKM-$$$$$$$$$$$$$$$$$$$$ + FM_hat,FH_hat,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, & +!$$$$$$$$$$$$$$$$$$-Inserted by PN, MS, PS, and SKM-$$$$$$$$$$$$$$$$$$$$ + P1000mb, psimh_stab_opt, psimh_unstab_opt, & +!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ + ids,ide, jds,jde, kds,kde, & + ims,ime, jms,jme, kms,kme, & + its,ite, jts,jte, kts,kte, & + isftcflx, iz0tlnd,scm_force_flux, & + ustm,ck,cka,cd,cda ) +!------------------------------------------------------------------- + IMPLICIT NONE +!------------------------------------------------------------------- + REAL, PARAMETER :: XKA=2.4E-5 + REAL, PARAMETER :: PRT=1. + + 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 ) :: ISFFLX + REAL, INTENT(IN ) :: SVP1,SVP2,SVP3,SVPT0 + REAL, INTENT(IN ) :: EP1,EP2,KARMAN,EOMEG,STBOLT + REAL, INTENT(IN ) :: P1000mb + +! + 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,FM,FH, & +!$$$$$$$$$$$$$$$$$$$$$$-Inserted by PN, MS, PS, and SKM-$$$$$$$$$$$$ + FM_hat,FH_hat +!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ + + REAL, DIMENSION( ims:ime ) , & + INTENT(INOUT) :: ZNT, & + ZOL, & + UST, & + CPM, & + CHS2, & + CQS2, & + CHS + + REAL, DIMENSION( ims:ime ) , & + INTENT(INOUT) :: FLHC,FLQC + + REAL, DIMENSION( ims:ime ) , & + INTENT(INOUT) :: & + QGH + + REAL, DIMENSION( ims:ime ) , & + INTENT(OUT) :: U10,V10, & + TH2,T2,Q2,QSFC,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 + + REAL, OPTIONAL, DIMENSION( ims:ime ) , & + INTENT(OUT) :: ck,cka,cd,cda + REAL, OPTIONAL, DIMENSION( ims:ime ) , & + INTENT(INOUT) :: USTM + + INTEGER, OPTIONAL, INTENT(IN ) :: ISFTCFLX, IZ0TLND + INTEGER, OPTIONAL, INTENT(IN ) :: SCM_FORCE_FLUX +!$$$$$$$$$$$$$$$$$$$$$-Inserted by PN, MS, PS, and SKM-$$$$$$$$$$$$$ + INTEGER, INTENT(IN) :: psimh_unstab_opt + INTEGER, INTENT(IN) :: psimh_stab_opt +!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ +! LOCAL VARS + + REAL, DIMENSION( its:ite ) :: ZA, & + THVX,ZQKL, & + ZQKLP1, & + THX,QX, & + PSIH2, & + PSIM2, & + PSIH10, & + PSIM10, & + DENOMQ, & + DENOMQ2, & + DENOMT2, & + WSPDI, & + GZ2OZ0, & + GZ10OZ0 +! + REAL, DIMENSION( its:ite ) :: & + RHOX,GOVRTH, & + TGDSA +! + REAL, DIMENSION( its:ite) :: SCR3,SCR4 + REAL, DIMENSION( its:ite ) :: THGB, PSFC +! + INTEGER :: KL + + INTEGER :: N,I,K,KK,L,NZOL,NK,NZOL2,NZOL10 + + REAL :: PL,THCON,TVCON,E1 + REAL :: ZL,TSKV,DTHVDZ,DTHVM,VCONV,RZOL,RZOL2,RZOL10,ZOL2,ZOL10 +!$$$$$$$$$$$$$$$$$$$$$$-Inserted by PN, MS, PS, and SKM-$$$$$$$$$$$$$$$$$$ + REAL :: PSIM_hat, PSIH_hat +!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ + REAL :: DTG,PSIX,DTTHX,PSIX10,PSIT,PSIT2,PSIQ,PSIQ2,PSIQ10 + REAL :: FLUXC,VSGD,Z0Q,VISC,RESTAR,CZIL,GZ0OZQ,GZ0OZT + REAL :: ZW, ZN1, ZN2 +! +! .... paj ... +! + REAL :: zolzz,zol0 +! REAL :: zolri,zolri2 +! REAL :: psih_stable,psim_stable,psih_unstable,psim_unstable +! REAL :: psih_stable_full,psim_stable_full,psih_unstable_full,psim_unstable_full + REAL :: zl2,zl10,z0t + REAL, DIMENSION( its:ite ) :: pq,pq2,pq10 + + +!------------------------------------------------------------------- + KL=kte + + DO i=its,ite +! PSFC cb + PSFC(I)=PSFCPA(I)/1000. + ENDDO +! +!----CONVERT GROUND TEMPERATURE TO POTENTIAL TEMPERATURE: +! + DO 5 I=its,ite + TGDSA(I)=TSK(I) +! PSFC cb +! THGB(I)=TSK(I)*(100./PSFC(I))**ROVCP + THGB(I)=TSK(I)*(P1000mb/PSFCPA(I))**ROVCP + 5 CONTINUE +! +!-----DECOUPLE FLUX-FORM VARIABLES TO GIVE U,V,T,THETA,THETA-VIR., +! T-VIR., QV, AND QC AT CROSS POINTS AND AT KTAU-1. +! +! *** NOTE *** +! THE BOUNDARY WINDS MAY NOT BE ADEQUATELY AFFECTED BY FRICTION, +! SO USE ONLY INTERIOR VALUES OF UX AND VX TO CALCULATE +! TENDENCIES. +! + 10 CONTINUE + +! DO 24 I=its,ite +! UX(I)=U1D(I) +! VX(I)=V1D(I) +! 24 CONTINUE + + 26 CONTINUE + +!.....SCR3(I,K) STORE TEMPERATURE, +! SCR4(I,K) STORE VIRTUAL TEMPERATURE. + + DO 30 I=its,ite +! PL cb + PL=P1D(I)/1000. + SCR3(I)=T1D(I) +! THCON=(100./PL)**ROVCP + THCON=(P1000mb*0.001/PL)**ROVCP + THX(I)=SCR3(I)*THCON + SCR4(I)=SCR3(I) + THVX(I)=THX(I) + QX(I)=0. + 30 CONTINUE +! + DO I=its,ite + QGH(I)=0. + FLHC(I)=0. + FLQC(I)=0. + CPM(I)=CP + ENDDO +! +! IF(IDRY.EQ.1)GOTO 80 + DO 50 I=its,ite + QX(I)=QV1D(I) + TVCON=(1.+EP1*QX(I)) + THVX(I)=THX(I)*TVCON + SCR4(I)=SCR3(I)*TVCON + 50 CONTINUE +! + DO 60 I=its,ite + E1=SVP1*EXP(SVP2*(TGDSA(I)-SVPT0)/(TGDSA(I)-SVP3)) +! for land points QSFC can come from previous time step + if(xland(i).gt.1.5.or.qsfc(i).le.0.0)QSFC(I)=EP2*E1/(PSFC(I)-E1) +! QGH CHANGED TO USE LOWEST-LEVEL AIR TEMP CONSISTENT WITH MYJSFC CHANGE +! Q2SAT = QGH IN LSM + E1=SVP1*EXP(SVP2*(T1D(I)-SVPT0)/(T1D(I)-SVP3)) + PL=P1D(I)/1000. + QGH(I)=EP2*E1/(PL-E1) + CPM(I)=CP*(1.+0.8*QX(I)) + 60 CONTINUE + 80 CONTINUE + +!-----COMPUTE THE HEIGHT OF FULL- AND HALF-SIGMA LEVELS ABOVE GROUND +! LEVEL, AND THE LAYER THICKNESSES. + + DO 90 I=its,ite + ZQKLP1(I)=0. + RHOX(I)=PSFC(I)*1000./(R*SCR4(I)) + 90 CONTINUE +! + DO 110 I=its,ite + ZQKL(I)=dz8w1d(I)+ZQKLP1(I) + 110 CONTINUE +! + DO 120 I=its,ite + ZA(I)=0.5*(ZQKL(I)+ZQKLP1(I)) + 120 CONTINUE +! + DO 160 I=its,ite + GOVRTH(I)=G/THX(I) + 160 CONTINUE + +!-----CALCULATE BULK RICHARDSON NO. OF SURFACE LAYER, ACCORDING TO +! AKB(1976), EQ(12). + + DO 260 I=its,ite + GZ1OZ0(I)=ALOG((ZA(I)+ZNT(I))/ZNT(I)) ! log((z+z0)/z0) + GZ2OZ0(I)=ALOG((2.+ZNT(I))/ZNT(I)) ! log((2+z0)/z0) + GZ10OZ0(I)=ALOG((10.+ZNT(I))/ZNT(I)) ! log((10+z0)z0) + IF((XLAND(I)-1.5).GE.0)THEN + ZL=ZNT(I) + ELSE + ZL=0.01 + ENDIF + WSPD(I)=SQRT(UX(I)*UX(I)+VX(I)*VX(I)) + + TSKV=THGB(I)*(1.+EP1*QSFC(I)) + DTHVDZ=(THVX(I)-TSKV) +! write(35,*) ZNT +! Convective velocity scale Vc and subgrid-scale velocity Vsg +! following Beljaars (1994, QJRMS) and Mahrt and Sun (1995, MWR) +! ... HONG Aug. 2001 +! +! 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 + fluxc = max(hfx(i)/rhox(i)/cp & + + ep1*tskv*qfx(i)/rhox(i),0.) + VCONV = vconvc*(g/tgdsa(i)*pblh(i)*fluxc)**.33 + else + IF(-DTHVDZ.GE.0)THEN + DTHVM=-DTHVDZ + ELSE + DTHVM=0. + ENDIF +! VCONV = 2.*SQRT(DTHVM) +! V3.7: reducing contribution in calm conditions + VCONV = SQRT(DTHVM) + endif +! Mahrt and Sun low-res correction + VSGD = 0.32 * (max(dx/5000.-1.,0.))**.33 + WSPD(I)=SQRT(WSPD(I)*WSPD(I)+VCONV*VCONV+vsgd*vsgd) + WSPD(I)=AMAX1(WSPD(I),0.1) + BR(I)=GOVRTH(I)*ZA(I)*DTHVDZ/(WSPD(I)*WSPD(I)) +! BR(I)=max(BR(I),-2.) +!$$$$$$$$$$$$$$$$$$$-Inserted by PN, MS, PS, and SKM-$$$$$$$$$$$$$$$$$$$$$$$$$$$$ + if(psimh_stab_opt == 1.or.psimh_stab_opt == 2.or.psimh_stab_opt ==5) then + BR(I)=min(BR(I),1.5) + ! write(31,*) BR(I) + + elseif(psimh_stab_opt == 6) then + BR(I)=min(BR(I),1.4) + ! write(311,*) BR(I) + + elseif(psimh_stab_opt == 3) then + BR(I) = min(BR(I),1.2) + ! write(3111,*) BR(I) + + elseif(psimh_stab_opt == 4) then + BR(I) = min(BR(I),0.9) + endif +!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ + +! IF PREVIOUSLY UNSTABLE, DO NOT LET INTO REGIMES 1 AND 2 + IF(MOL(I).LT.0.)BR(I)=AMIN1(BR(I),0.0) +!jdf + RMOL(I)=-GOVRTH(I)*DTHVDZ*ZA(I)*KARMAN +!jdf + + 260 CONTINUE + +! +!-----DIAGNOSE BASIC PARAMETERS FOR THE APPROPRIATED STABILITY CLASS: +! +! +! THE STABILITY CLASSES ARE DETERMINED BY BR (BULK RICHARDSON NO.) +! AND HOL (HEIGHT OF PBL/MONIN-OBUKHOV LENGTH). +! +! CRITERIA FOR THE CLASSES ARE AS FOLLOWS: +! +! 1. BR .GE. 0.0; +! REPRESENTS NIGHTTIME STABLE CONDITIONS (REGIME=1), +! +! 3. BR .EQ. 0.0 +! REPRESENTS FORCED CONVECTION CONDITIONS (REGIME=3), +! +! 4. BR .LT. 0.0 +! REPRESENTS FREE CONVECTION CONDITIONS (REGIME=4). +! +!CCCCC + + DO 320 I=its,ite +! + if (br(I).gt.0) then + if (br(I).gt.250.0) then +!$$$$$$$$$$$$$$$$$$$$-Inserted by PN, MS, PS, and SKM-$$$$$$$$$$$$$$$$$$$$ + zol(I)=zolri(250.0,ZA(I),ZNT(I),psimh_stab_opt,psimh_unstab_opt) + else + zol(I)=zolri(br(I),ZA(I),ZNT(I),psimh_stab_opt,psimh_unstab_opt) + endif + endif +! + if (br(I).lt.0) then + IF(UST(I).LT.0.001)THEN + ZOL(I)=BR(I)*GZ1OZ0(I) + ELSE + if (br(I).lt.-250.0) then + zol(I)=zolri(-250.0,ZA(I),ZNT(I),psimh_stab_opt,psimh_unstab_opt) + else + zol(I)=zolri(br(I),ZA(I),ZNT(I),psimh_stab_opt,psimh_unstab_opt) +!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ + endif + ENDIF +! zol(I) = max(zol(I),-100.) + endif +! +! ... paj: compute integrated similarity functions. +! + zolzz=zol(I)*(za(I)+znt(I))/za(I) ! (z+z0/L + zol10=zol(I)*(10.+znt(I))/za(I) ! (10+z0)/L + zol2=zol(I)*(2.+znt(I))/za(I) ! (2+z0)/L + zol0=zol(I)*znt(I)/za(I) ! z0/L + ZL2=(2.)/ZA(I)*ZOL(I) ! 2/L + ZL10=(10.)/ZA(I)*ZOL(I) ! 10/L + + IF((XLAND(I)-1.5).LT.0.)THEN + ZL=(0.01)/ZA(I)*ZOL(I) ! (0.01)/L + ELSE + ZL=ZOL0 ! z0/L + ENDIF + + IF(BR(I).LT.0.)GOTO 310 ! go to unstable regime (class 4) + IF(BR(I).EQ.0.)GOTO 280 ! go to neutral regime (class 3) +! +!-----CLASS 1; STABLE (NIGHTTIME) CONDITIONS: +! + REGIME(I)=1. +! ... contains various expressions for psim and psih for stable and convective conditions + +!$$$$$$$$$$$$$$$$$$-Inserted by PN, MS, PS, and SKM-$$$$$$$$$$$$$$$$$$$$$$$$$$ + + if(psimh_stab_opt == 4.) then + if(zolzz.le.1.) then + psim(I)=psim_stable(zolzz,psimh_stab_opt) & + -psim_stable(zol0,psimh_stab_opt) + psih(I)=psih_stable(zolzz,psimh_stab_opt) & + -psih_stable(zol0,psimh_stab_opt) + else + psim(I) = -5.-psim_stable(zol0,psimh_stab_opt)+& + (-4.*log(zolzz)-zolzz+1.) + psih(I) = -5.-psih_stable(zol0,psimh_stab_opt)+& + (-4.*log(zolzz)-zolzz+1.) + endif + else + psim(I)=psim_stable(zolzz,psimh_stab_opt)-& + psim_stable(zol0,psimh_stab_opt) + psih(I)=psih_stable(zolzz,psimh_stab_opt)-& + psih_stable(zol0,psimh_stab_opt) + endif +! + if(psimh_stab_opt == 4.) then + if(zol10.le.1.) then + psim10(I)=psim_stable(zol10,psimh_stab_opt)-& + psim_stable(zol0,psimh_stab_opt) + psih10(I)=psih_stable(zol10,psimh_stab_opt)-& + psih_stable(zol0,psimh_stab_opt) + else + psim10(I)=-5.-psim_stable(zol0,psimh_stab_opt)+& + (-4.*log(zol10)-zol10+1.) + psih10(I)=-5.-psih_stable(zol0,psimh_stab_opt)+& + (-4.*log(zol10)-zol10+1.) + endif + else + psim10(I)=psim_stable(zol10,psimh_stab_opt) & + -psim_stable(zol0,psimh_stab_opt) + psih10(I)=psih_stable(zol10,psimh_stab_opt) & + -psih_stable(zol0,psimh_stab_opt) + endif +! + if(psimh_stab_opt == 4.) then + if(zol2.le.1.) then + psim2(I)=psim_stable(zol2,psimh_stab_opt)- & + psim_stable(zol0,psimh_stab_opt) + psih2(I)=psih_stable(zol2,psimh_stab_opt)- & + psih_stable(zol0,psimh_stab_opt) + else + psim2(I)=-5.-psim_stable(zol0,psimh_stab_opt)+& + (-4.*log(zol2)-zol2+1.) + psih2(I)=-5.-psih_stable(zol0,psimh_stab_opt)+& + (-4.*log(zol2)-zol2+1.) + endif + else + psim2(I)=psim_stable(zol2,psimh_stab_opt)- & + psim_stable(zol0,psimh_stab_opt) + psih2(I)=psih_stable(zol2,psimh_stab_opt)- & + psih_stable(zol0,psimh_stab_opt) + endif +!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ + +! ... paj: preparations to compute PSIQ. Follows CB05+Carlson Boland JAM 1978. + +!$$$$$$$$$$$$-Inserted by PN, MS, PS, and SKM-$$$$$$$$$$$$$$$$$$$$$$$$ + if(psimh_stab_opt == 4.) then + if(zol(I).le.1.) then + pq(I)=psih_stable(zol(I),psimh_stab_opt)-& + psih_stable(zl,psimh_stab_opt) + else + pq(I)=-5.-psih_stable(zl,psimh_stab_opt)+& + (-4.*log(zol(I))-zol(I)+1.) + endif + else + pq(I)=psih_stable(zol(I),psimh_stab_opt)-& + psih_stable(zl,psimh_stab_opt) + endif +! + if(psimh_stab_opt == 4.) then + if(zl2.le.1.) then + pq2(I)=psih_stable(zl2,psimh_stab_opt)-& + psih_stable(zl,psimh_stab_opt) + else + pq2(I) = -5.-psih_stable(zl,psimh_stab_opt)+& + (-4.*log(zl2)-zl2+1.) + endif + else + pq2(I)=psih_stable(zl2,psimh_stab_opt)-& + psih_stable(zl,psimh_stab_opt) + endif +! + if(psimh_stab_opt == 4.) then + if(zl10.le.1.) then + pq10(I)=psih_stable(zl10,psimh_stab_opt)-& + psih_stable(zl,psimh_stab_opt) + else + pq10(I)=-5.-psih_stable(zl,psimh_stab_opt)+& + (-4.*log(zl10)-zl10+1.) + endif + else + pq10(I)=psih_stable(zl10,psimh_stab_opt)-& + psih_stable(zl,psimh_stab_opt) + endif +!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ + +! 1.0 over Monin-Obukhov length + RMOL(I)=ZOL(I)/ZA(I) +! + + GOTO 320 +! +!-----CLASS 3; FORCED CONVECTION: +! + 280 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) +! +! paj: preparations to compute PSIQ. +! + pq(I)=PSIH(I) + pq2(I)=PSIH2(I) + pq10(I)=0. +! + ZOL(I)=0. + RMOL(I) = ZOL(I)/ZA(I) + + GOTO 320 +! +!-----CLASS 4; FREE CONVECTION: +! + 310 CONTINUE + REGIME(I)=4. +! +! ... paj: PSIM and PSIH ... +!$$$$$$$$$$$$$$$$$$$$$$$$-Inserted by PN, MS, PS, and SKM-$$$$$$$$$$$$$$$$ + if(psimh_unstab_opt == 3.) then + if(zolzz.gt.-1.574) then + psim(I)=psim_unstable(zolzz,psimh_unstab_opt)& + -psim_unstable(zol0,psimh_unstab_opt) + else + psim(I)=psim_unstable(-1.574,psimh_unstab_opt)-& + psim_unstable(zol0,psimh_unstab_opt)+ & + log(zolzz/(-1.574))-1.14*((- zolzz)**& + (1./3.)-(1.574)**(1./3.)) + endif + else + psim(I)=psim_unstable(zolzz,psimh_unstab_opt)- & + psim_unstable(zol0,psimh_unstab_opt) + endif +! + if(psimh_unstab_opt == 3.) then + if(zolzz.gt.-0.465) then + psih(I)=psih_unstable(zolzz,psimh_unstab_opt)-& + psih_unstable(zol0,psimh_unstab_opt) + else + psih(I)=psih_unstable(-0.465,psimh_unstab_opt)-& + psih_unstable(zol0,psimh_unstab_opt)+ & + log(zolzz/(-0.465))-0.8*((0.465)**(-1./3.) & + -(- zolzz)**(-1./3.)) + endif + else + psih(I)=psih_unstable(zolzz,psimh_unstab_opt)-& + psih_unstable(zol0,psimh_unstab_opt) + endif +! + if(psimh_unstab_opt == 3.) then + if(zol10.gt.-1.574) then + psim10(I)=psim_unstable(zol10,psimh_unstab_opt)-& + psim_unstable(zol0,psimh_unstab_opt) + else + psim10(I)=psim_unstable(-1.574,psimh_unstab_opt)-& + psim_unstable(zol0,psimh_unstab_opt)+ & + log(zol10/(-1.574))-1.14*((- zol10)**(1./3.) & + -(1.574)**(1./3.)) + endif + else + psim10(I)=psim_unstable(zol10,psimh_unstab_opt)-& + psim_unstable(zol0,psimh_unstab_opt) + endif + + if(psimh_unstab_opt == 3.) then + if(zol10.gt.-0.465) then + psih10(I)=psih_unstable(zol10,psimh_unstab_opt)-& + psih_unstable(zol0,psimh_unstab_opt) + else + psih10(I)=psih_unstable(-0.465,psimh_unstab_opt)-& + psih_unstable(zol0,psimh_unstab_opt)+ & + log(zol10/(-0.465))-0.8*((0.465)**(-1./3.) & + -(- zol10)**(-1./3.)) + endif + else + psih10(I)=psih_unstable(zol10,psimh_unstab_opt)-& + psih_unstable(zol0,psimh_unstab_opt) + endif +! + if(psimh_unstab_opt == 3.) then + if(zol2.gt.-1.574) then + psim2(I)=psim_unstable(zol2,psimh_unstab_opt)-& + psim_unstable(zol0,psimh_unstab_opt) + else + psim2(I)=psim_unstable(-1.574,psimh_unstab_opt)-& + psim_unstable(zol0,psimh_unstab_opt)+ & + log(zol2/(-1.574))-1.14*((- zol2)**(1./3.) & + -(1.574)**(1./3.)) + endif + else + psim2(I)=psim_unstable(zol2,psimh_unstab_opt)-& + psim_unstable(zol0,psimh_unstab_opt) + endif +! + if(psimh_unstab_opt == 3.) then + if(zol2.gt.-0.465) then + psih2(I)=psih_unstable(zol2,psimh_unstab_opt)-& + psih_unstable(zol0,psimh_unstab_opt) + else + psih2(I)=psih_unstable(-0.465,psimh_unstab_opt)-& + psih_unstable(zol0,psimh_unstab_opt)+ & + log(zol2/(-0.465))-0.8*((0.465)**(-1./3.) & + -(- zol2)**(-1./3.)) + endif + else + psih2(I)=psih_unstable(zol2,psimh_unstab_opt)-& + psih_unstable(zol0,psimh_unstab_opt) + endif +!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ +! ... paj: preparations to compute PSIQ +!$$$$$$$$$$$$$$$$$-Inserted by PN, MS, PS, and SKM-$$$$$$$$$$$$$$$$$$$$ + if(psimh_unstab_opt == 3.) then + if(zol(I).gt.-0.465) then + pq(I)=psih_unstable(zol(I),psimh_unstab_opt)-& + psih_unstable(zl,psimh_unstab_opt) + else + pq(I)=psih_unstable(-0.465,psimh_unstab_opt)-& + psih_unstable(zl,psimh_unstab_opt)+& + log(zol(I)/(-0.465))-0.8*((0.465)**(-1./3.) & + -(-zol(I))**(-1./3.)) + endif + else + pq(I)=psih_unstable(zol(I),psimh_unstab_opt)-& + psih_unstable(zl,psimh_unstab_opt) + endif +! + if(psimh_unstab_opt == 3.) then + if(zl2.gt.-0.465) then + pq2(I)=psih_unstable(zl2,psimh_unstab_opt)-& + psih_unstable(zl,psimh_unstab_opt) + else + pq2(I)=psih_unstable(-0.465,psimh_unstab_opt)-& + psih_unstable(zl,psimh_unstab_opt)+ & + log(zol(I)/(-0.465))-0.8*((0.465)**(-1./3.) & + -(- zl2)**(-1./3.)) + endif + else + pq2(I)=psih_unstable(zl2,psimh_unstab_opt)-& + psih_unstable(zl,psimh_unstab_opt) + endif +! + if(psimh_unstab_opt == 3.) then + if(zl10.gt.-0.465) then + pq10(I)=psih_unstable(zl10,psimh_unstab_opt)-& + psih_unstable(zl,psimh_unstab_opt) + else + pq10(I)=psih_unstable(-0.465,psimh_unstab_opt)- & + psih_unstable(zl,psimh_unstab_opt)+ & + log(zl10/(-0.465))-0.8*((0.465)**(-1./3.) & + -(- zl10)**(-1./3.)) + endif + else + pq10(I)=psih_unstable(zl10,psimh_unstab_opt)-& + psih_unstable(zl,psimh_unstab_opt) + endif +!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ +!========================================================================= +!---LIMIOT PSIH AND PSIM IN THE CASE OF THIN LAYERS AND HIGH ROUGHNESS +!--- THIS PREVENTS DENOMINATOR IN FLUXES FROM GETTING TOO SMALL + PSIH(I)=AMIN1(PSIH(I),0.9*GZ1OZ0(I)) + PSIM(I)=AMIN1(PSIM(I),0.9*GZ1OZ0(I)) + PSIH2(I)=AMIN1(PSIH2(I),0.9*GZ2OZ0(I)) + PSIM10(I)=AMIN1(PSIM10(I),0.9*GZ10OZ0(I)) +! +! AHW: mods to compute ck, cd + PSIH10(I)=AMIN1(PSIH10(I),0.9*GZ10OZ0(I)) + + RMOL(I) = ZOL(I)/ZA(I) + + 320 CONTINUE +! +!-----COMPUTE THE FRICTIONAL VELOCITY: +! ZA(1982) EQS(2.60),(2.61). +! + DO 330 I=its,ite + DTG=THX(I)-THGB(I) + PSIX=GZ1OZ0(I)-PSIM(I) + PSIX10=GZ10OZ0(I)-PSIM10(I) +!$$$$$$$$$$$$$$$$$$$$$$-Inserted by PN, MS, PS, and SKM-$$$$$$$$$ + PSIM_hat = PSIM(I) + PSIH_hat = PSIH(I) +!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ + +! LOWER LIMIT ADDED TO PREVENT LARGE FLHC IN SOIL MODEL +! ACTIVATES IN UNSTABLE CONDITIONS WITH THIN LAYERS OR HIGH Z0 +! PSIT=AMAX1(GZ1OZ0(I)-PSIH(I),2.) + PSIT=GZ1OZ0(I)-PSIH(I) + PSIT2=GZ2OZ0(I)-PSIH2(I) +! + IF((XLAND(I)-1.5).GE.0)THEN + ZL=ZNT(I) + ELSE + ZL=0.01 + ENDIF +! + PSIQ=ALOG(KARMAN*UST(I)*ZA(I)/XKA+ZA(I)/ZL)-pq(I) + PSIQ2=ALOG(KARMAN*UST(I)*2./XKA+2./ZL)-pq2(I) + +! AHW: mods to compute ck, cd + PSIQ10=ALOG(KARMAN*UST(I)*10./XKA+10./ZL)-pq10(I) + +! V3.7: using Fairall 2003 to compute z0q and z0t over water: +! adapted from module_sf_mynn.F + IF ( (XLAND(I)-1.5).GE.0. ) THEN + VISC=(1.32+0.009*(SCR3(I)-273.15))*1.E-5 + RESTAR=UST(I)*ZNT(I)/VISC + Z0T = (5.5e-5)*(RESTAR**(-0.60)) + Z0T = MIN(Z0T,1.0e-4) + Z0T = MAX(Z0T,2.0e-9) + Z0Q = Z0T + +! following paj: + zolzz=zol(I)*(za(I)+z0t)/za(I) ! (z+z0t)/L + zol10=zol(I)*(10.+z0t)/za(I) ! (10+z0t)/L + zol2=zol(I)*(2.+z0t)/za(I) ! (2+z0t)/L + zol0=zol(I)*z0t/za(I) ! z0t/L +! + if(zol(I).gt.0.) then +!$$$$$$$$$$$$$$$$$$$$$$$$$-Inserted by PN, MS, PS, and SKM-$$$$$$$$$$$$$ + if(psimh_stab_opt == 4.) then + if(zol(I).le.1.) then + psih(I)=psih_stable(zolzz,psimh_stab_opt)& + -psih_stable(zol0,psimh_stab_opt) + psih10(I)=psih_stable(zol10,psimh_stab_opt)& + -psih_stable(zol0,psimh_stab_opt) + psih2(I)=psih_stable(zol2,psimh_stab_opt)& + -psih_stable(zol0,psimh_stab_opt) + else + psih(I)=-5.-psih_stable(zol0,psimh_stab_opt)+& + (-4.*log(zolzz)-zolzz+1.) + psih10(I) = -5.-psih_stable(zol0,psimh_stab_opt)+& + (-4.*log(zol10)-zol10+1.) + psih2(I) = -5.-psih_stable(zol0,psimh_stab_opt)+& + (-4.*log(zol2)-zol2+1.) + endif + else + psih(I)=psih_stable(zolzz,psimh_stab_opt)& + -psih_stable(zol0,psimh_stab_opt) + psih10(I)=psih_stable(zol10,psimh_stab_opt)& + -psih_stable(zol0,psimh_stab_opt) + psih2(I)=psih_stable(zol2,psimh_stab_opt)& + -psih_stable(zol0,psimh_stab_opt) + endif +!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ + else + if(zol(I).eq.0) then + psih(I)=0. + psih10(I)=0. + psih2(I)=0. + else +!$$$$$$$$$$$$$$$$$$$$$$$$-Inserted by PN, MS, PS, and SKM-$$$$$$$$$$$$$$$$$$$ + if(psimh_unstab_opt == 3.) then + if(zolzz.gt.-0.465) then + psih(I)=psih_unstable(zolzz,psimh_unstab_opt)& + -psih_unstable(zol0,psimh_unstab_opt) + else + psih(I)=psih_unstable(-0.465,psimh_unstab_opt)& + -psih_unstable(zol0,psimh_unstab_opt)+& + log(zolzz/(-0.465))-0.8*((0.465)**(-1./3.)-& + (- zolzz)**(-1./3.)) + endif + else + psih(I)=psih_unstable(zolzz,psimh_unstab_opt)& + -psih_unstable(zol0,psimh_unstab_opt) + endif +! + if(psimh_unstab_opt == 3.) then + if(zol10.gt.-0.465) then + psih10(I)=psih_unstable(zol10,psimh_unstab_opt)& + -psih_unstable(zol0,psimh_unstab_opt) + else + psih10(I)=psih_unstable(-0.465,psimh_unstab_opt)& + -psih_unstable(zol0,psimh_unstab_opt)& + +log(zol10/(-0.465))-0.8*((0.465)**(-1./3.)-& + (- zol10)**(-1./3.)) + endif + else + psih10(I)=psih_unstable(zol10,psimh_unstab_opt)& + -psih_unstable(zol0,psimh_unstab_opt) + endif + if(psimh_unstab_opt == 3.) then + if(zol2.gt.-0.465) then + psih2(I)=psih_unstable(zol2,psimh_unstab_opt)& + -psih_unstable(zol0,psimh_unstab_opt) + else + psih2(I)=psih_unstable(-0.465,psimh_unstab_opt)& + -psih_unstable(zol0,psimh_unstab_opt)& + +log(zol2/(-0.465))-0.8*((0.465)**(-1./3.)-& + (- zol2)**(-1./3.)) + endif + else + psih2(I)=psih_unstable(zol2,psimh_unstab_opt)& + -psih_unstable(zol0,psimh_unstab_opt) + endif +!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ + endif + endif + PSIT=ALOG((ZA(I)+z0t)/Z0t)-PSIH(I) + PSIT2=ALOG((2.+z0t)/Z0t)-PSIH2(I) + + zolzz=zol(I)*(za(I)+z0q)/za(I) ! (z+z0q)/L + zol10=zol(I)*(10.+z0q)/za(I) ! (10+z0q)/L + zol2=zol(I)*(2.+z0q)/za(I) ! (2+z0q)/L + zol0=zol(I)*z0q/za(I) ! z0q/L +! + if(zol(I).gt.0.) then +!$$$$$$$$$$$$$$$$$$$$$$$$-Inserted by PN, MS, PS, and SKM-$$$$$$$$$$$$$$$$$$$ + if(psimh_stab_opt == 4.) then + if(zol(I).le.1.) then + psih(I)=psih_stable(zolzz,psimh_stab_opt)& + -psih_stable(zol0,psimh_stab_opt) + psih10(I)=psih_stable(zol10,psimh_stab_opt)& + -psih_stable(zol0,psimh_stab_opt) + psih2(I)=psih_stable(zol2,psimh_stab_opt)& + -psih_stable(zol0,psimh_stab_opt) + else + psih(I)=-5.-psih_stable(zol0,psimh_stab_opt)+& + (-4.*log(zolzz)-zolzz+1.) + psih10(I)=-5.-psih_stable(zol0,psimh_stab_opt)+& + (-4.*log(zol10)-zol10+1.) + psih2(I) = -5.-psih_stable(zol0,psimh_stab_opt)+& + (-4.*log(zol2)-zol2+1.) + endif + else + psih(I)=psih_stable(zolzz,psimh_stab_opt)& + -psih_stable(zol0,psimh_stab_opt) + psih10(I)=psih_stable(zol10,psimh_stab_opt)& + -psih_stable(zol0,psimh_stab_opt) + psih2(I)=psih_stable(zol2,psimh_stab_opt)& + -psih_stable(zol0,psimh_stab_opt) + endif +!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ + else + if (zol(I).eq.0) then + psih(I)=0. + psih10(I)=0. + psih2(I)=0. + else +!$$$$$$$$$$$$$$$$$$$$$$$$-Inserted by PN, MS, PS, and SKM-$$$$$$$$$$$$$$$$$$ + if(psimh_unstab_opt == 3.) then + if(zolzz.gt.-0.465) then + psih(I)=psih_unstable(zolzz,psimh_unstab_opt)& + -psih_unstable(zol0,psimh_unstab_opt) + else + psih(I)=psih_unstable(-0.465,psimh_unstab_opt)& + -psih_unstable(zol0,psimh_unstab_opt)& + +log(zolzz/(-0.465))-0.8*((0.465)**(-1./3.)-& + (- zolzz)**(-1./3.)) + endif + else + psih(I)=psih_unstable(zolzz,psimh_unstab_opt)& + -psih_unstable(zol0,psimh_unstab_opt) + endif +! + if(psimh_unstab_opt == 3.) then + if(zol10.gt.-0.465) then + psih10(I)=psih_unstable(zol10,psimh_unstab_opt)& + -psih_unstable(zol0,psimh_unstab_opt) + else + psih10(I)=psih_unstable(-0.465,psimh_unstab_opt)& + -psih_unstable(zol0,psimh_unstab_opt)& + +log(zol10/(-0.465))-0.8*((0.465)**(-1./3.)-& + (- zol10)**(-1./3.)) + endif + else + psih10(I)=psih_unstable(zol10,psimh_unstab_opt)& + -psih_unstable(zol0,psimh_unstab_opt) + endif +! + if(psimh_unstab_opt == 3.) then + if(zol2.gt.-0.465) then + psih2(I)=psih_unstable(zol2,psimh_unstab_opt)& + -psih_unstable(zol0,psimh_unstab_opt) + else + psih2(I)=psih_unstable(-0.465,psimh_unstab_opt)& + -psih_unstable(zol0,psimh_unstab_opt)+& + log(zol2/(-0.465))-0.8*((0.465)**(-1./3.)-& + (-zol2)**(-1./3.)) + endif + else + psih2(I)=psih_unstable(zol2,psimh_unstab_opt)& + -psih_unstable(zol0,psimh_unstab_opt) + endif +!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ + endif + endif +! + PSIQ=ALOG((ZA(I)+z0q)/Z0q)-PSIH(I) + PSIQ2=ALOG((2.+z0q)/Z0q)-PSIH2(I) + PSIQ10=ALOG((10.+z0q)/Z0q)-PSIH10(I) + ENDIF + + IF ( PRESENT(ISFTCFLX) ) THEN + IF ( ISFTCFLX.EQ.1 .AND. (XLAND(I)-1.5).GE.0. ) THEN +! v3.1 +! Z0Q = 1.e-4 + 1.e-3*(MAX(0.,UST(I)-1.))**2 +! hfip1 +! Z0Q = 0.62*2.0E-5/UST(I) + 1.E-3*(MAX(0.,UST(I)-1.5))**2 +! v3.2 + Z0Q = 1.e-4 +! +! ... paj: recompute psih for z0q +! + zolzz=zol(I)*(za(I)+z0q)/za(I) ! (z+z0q)/L + zol10=zol(I)*(10.+z0q)/za(I) ! (10+z0q)/L + zol2=zol(I)*(2.+z0q)/za(I) ! (2+z0q)/L + zol0=zol(I)*z0q/za(I) ! z0q/L +! + if(zol(I).gt.0.) then +!$$$$$$$$$$$$$$$$$$$$$$$$-Inserted by PN, MS, PS, and SKM-$$$$$$$$$$$$$$$$$$$ + if(psimh_stab_opt == 4.) then + if(zol(I).le.1.) then + psih(I)=psih_stable(zolzz,psimh_stab_opt)& + -psih_stable(zol0,psimh_stab_opt) + psih10(I)=psih_stable(zol10,psimh_stab_opt) & + -psih_stable(zol0,psimh_stab_opt) + psih2(I)=psih_stable(zol2,psimh_stab_opt) & + -psih_stable(zol0,psimh_stab_opt) + else + psih(I)=-5.-psih_stable(zol0,psimh_stab_opt)+& + (-4.*log(zolzz)-zolzz+1.) + psih10(I)=-5.-psih_stable(zol0,psimh_stab_opt)+& + (-4.*log(zol10)-zol10+1.) + psih2(I)=-5.-psih_stable(zol0,psimh_stab_opt)+& + (-4.*log(zol2)-zol2+1.) + endif + else + psih(I)=psih_stable(zolzz,psimh_stab_opt)& + -psih_stable(zol0,psimh_stab_opt) + psih10(I)=psih_stable(zol10,psimh_stab_opt) & + -psih_stable(zol0,psimh_stab_opt) + psih2(I)=psih_stable(zol2,psimh_stab_opt) & + -psih_stable(zol0,psimh_stab_opt) + endif +!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ + else + if (zol(I).eq.0) then + psih(I)=0. + psih10(I)=0. + psih2(I)=0. + else +!$$$$$$$$$$$$$$$$$$$$$$$$-Inserted by PN, MS, PS, and SKM-$$$$$$$$$$$$$$$$ + if(psimh_unstab_opt == 3.) then + if(zolzz.gt.-0.465) then + psih(I)=psih_unstable(zolzz,psimh_unstab_opt) & + -psih_unstable(zol0,psimh_unstab_opt) + else + psih(I)=psih_unstable(-0.465,psimh_unstab_opt) & + -psih_unstable(zol0,psimh_unstab_opt)+& + log(zolzz/(-0.465))-0.8*((0.465)**(-1./3.)-& + (- zolzz)**(-1./3.)) + endif + else + psih(I)=psih_unstable(zolzz,psimh_unstab_opt) & + -psih_unstable(zol0,psimh_unstab_opt) + endif +! + if(psimh_unstab_opt == 3.) then + if(zol10.gt.-0.465) then + psih10(I)=psih_unstable(zol10,psimh_unstab_opt) & + -psih_unstable(zol0,psimh_unstab_opt) + else + psih10(I)=psih_unstable(-0.465,psimh_unstab_opt) & + -psih_unstable(zol0,psimh_unstab_opt)& + +log(zol10/(-0.465))-0.8*((0.465)**(-1./3.)-& + (- zol10)**(-1./3.)) + endif + else + psih10(I)=psih_unstable(zol10,psimh_unstab_opt) & + -psih_unstable(zol0,psimh_unstab_opt) + endif +! + if(psimh_unstab_opt == 3.) then + if(zol2.gt.-0.465) then + psih2(I)=psih_unstable(zol2,psimh_unstab_opt) & + -psih_unstable(zol0,psimh_unstab_opt) + else + psih2(I)=psih_unstable(-0.465,psimh_unstab_opt) & + -psih_unstable(zol0,psimh_unstab_opt)& + +log(zol2/(-0.465))-0.8*((0.465)**(-1./3.)-& + (- zol2)**(-1./3.)) + endif + else + psih2(I)=psih_unstable(zol2,psimh_unstab_opt) & + -psih_unstable(zol0,psimh_unstab_opt) + endif +!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ + endif + endif +! + PSIQ=ALOG((ZA(I)+z0q)/Z0Q)-PSIH(I) + PSIT=PSIQ + PSIQ2=ALOG((2.+z0q)/Z0Q)-PSIH2(I) + PSIQ10=ALOG((10.+z0q)/Z0Q)-PSIH10(I) + PSIT2=PSIQ2 + ENDIF + IF ( ISFTCFLX.EQ.2 .AND. (XLAND(I)-1.5).GE.0. ) THEN +! AHW: Garratt formula: Calculate roughness Reynolds number +! Kinematic viscosity of air (linear approc to +! temp dependence at sea level) +! GZ0OZT and GZ0OZQ are based off formulas from Brutsaert (1975), which +! Garratt (1992) used with values of k = 0.40, Pr = 0.71, and Sc = 0.60 + VISC=(1.32+0.009*(SCR3(I)-273.15))*1.E-5 +!! VISC=1.5E-5 + RESTAR=UST(I)*ZNT(I)/VISC + GZ0OZT=0.40*(7.3*SQRT(SQRT(RESTAR))*SQRT(0.71)-5.) +! +! ... paj: compute psih for z0t for temperature ... +! + z0t=znt(I)/exp(GZ0OZT) +! + zolzz=zol(I)*(za(I)+z0t)/za(I) ! (z+z0t)/L + zol10=zol(I)*(10.+z0t)/za(I) ! (10+z0t)/L + zol2=zol(I)*(2.+z0t)/za(I) ! (2+z0t)/L + zol0=zol(I)*z0t/za(I) ! z0t/L +! + if(zol(I).gt.0.) then +!$$$$$$$$$$$$$$$$$$$$$$$$-Inserted by PN, MS, PS, and SKM-$$$$$$$$$$$$$$$$$$$ + if(psimh_stab_opt == 4.) then + if(zol(I).le.1.) then + psih(I)=psih_stable(zolzz,psimh_stab_opt) & + -psih_stable(zol0,psimh_stab_opt) + psih10(I)=psih_stable(zol10,psimh_stab_opt) & + -psih_stable(zol0,psimh_stab_opt) + psih2(I)=psih_stable(zol2,psimh_stab_opt) & + -psih_stable(zol0,psimh_stab_opt) + else + psih(I)=-5.-psih_stable(zol0,psimh_stab_opt)+& + (-4.*log(zolzz)-zolzz+1.) + psih10(I)=-5.-psih_stable(zol0,psimh_stab_opt)+& + (-4.*log(zol10)-zol10+1.) + psih2(I)=-5.-psih_stable(zol0,psimh_stab_opt)+& + (-4.*log(zol2)-zol2+1.) + endif + else + psih(I)=psih_stable(zolzz,psimh_stab_opt) & + -psih_stable(zol0,psimh_stab_opt) + psih10(I)=psih_stable(zol10,psimh_stab_opt) & + -psih_stable(zol0,psimh_stab_opt) + psih2(I)=psih_stable(zol2,psimh_stab_opt) & + -psih_stable(zol0,psimh_stab_opt) + endif +!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ + else + if (zol(I).eq.0) then + psih(I)=0. + psih10(I)=0. + psih2(I)=0. + else +!$$$$$$$$$$$$$$$$$$$$$$$$-Inserted by PN, MS, PS, and SKM-$$$$$$$$$$$$$$$$$ + if(psimh_unstab_opt == 3.) then + if(zolzz.gt.-0.465) then + psih(I)=psih_unstable(zolzz,psimh_unstab_opt) & + -psih_unstable(zol0,psimh_unstab_opt) + else + psih(I)=psih_unstable(-0.465,psimh_unstab_opt) & + -psih_unstable(zol0,psimh_unstab_opt) & + +log(zolzz/(-0.465))-0.8*((0.465)**(-1./3.)-& + (- zolzz)**(-1./3.)) + endif + else + psih(I)=psih_unstable(zolzz,psimh_unstab_opt) & + -psih_unstable(zol0,psimh_unstab_opt) + endif +! + if(psimh_unstab_opt == 3.) then + if(zol10.gt.-0.465) then + psih10(I)=psih_unstable(zol10,psimh_unstab_opt) & + -psih_unstable(zol0,psimh_unstab_opt) + else + psih10(I)=psih_unstable(-0.465,psimh_unstab_opt) & + -psih_unstable(zol0,psimh_unstab_opt)& + +log(zol10/(-0.465))-0.8*((0.465)**(-1./3.)-& + (- zol10)**(-1./3.)) + endif + else + psih10(I)=psih_unstable(zol10,psimh_unstab_opt) & + -psih_unstable(zol0,psimh_unstab_opt) + endif +! + if(psimh_unstab_opt == 3.) then + if(zol2.gt.-0.465) then + psih2(I)=psih_unstable(zol2,psimh_unstab_opt) & + -psih_unstable(zol0,psimh_unstab_opt) + else + psih2(I)=psih_unstable(-0.465,psimh_unstab_opt) & + -psih_unstable(zol0,psimh_unstab_opt) & + +log(zol2/(-0.465))-0.8*((0.465)**(-1./3.)-& + (- zol2)**(-1./3.)) + endif + else + psih2(I)=psih_unstable(zol2,psimh_unstab_opt) & + -psih_unstable(zol0,psimh_unstab_opt) + endif +!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ + endif + endif +! +! PSIT=GZ1OZ0(I)-PSIH(I)+RESTAR2 +! PSIT2=GZ2OZ0(I)-PSIH2(I)+RESTAR2 + PSIT=ALOG((ZA(I)+z0t)/Z0t)-PSIH(I) + PSIT2=ALOG((2.+z0t)/Z0t)-PSIH2(I) +! + GZ0OZQ=0.40*(7.3*SQRT(SQRT(RESTAR))*SQRT(0.60)-5.) + z0q=znt(I)/exp(GZ0OZQ) +! + zolzz=zol(I)*(za(I)+z0q)/za(I) ! (z+z0q)/L + zol10=zol(I)*(10.+z0q)/za(I) ! (10+z0q)/L + zol2=zol(I)*(2.+z0q)/za(I) ! (2+z0q)/L + zol0=zol(I)*z0q/za(I) ! z0q/L +! + if(zol(I).gt.0.) then +!$$$$$$$$$$$$$$$$$$$$$$$$-Inserted by PN, MS, PS, and SKM-$$$$$$$$$$$$$$$$$$$ + if(psimh_stab_opt == 4.) then + if(zol(I).le.1.) then + psih(I)=psih_stable(zolzz,psimh_stab_opt) & + -psih_stable(zol0,psimh_stab_opt) + psih10(I)=psih_stable(zol10,psimh_stab_opt) & + -psih_stable(zol0,psimh_stab_opt) + psih2(I)=psih_stable(zol2,psimh_stab_opt) & + -psih_stable(zol0,psimh_stab_opt) + else + psih(I)=-5.-psih_stable(zol0,psimh_stab_opt)+& + (-4.*log(zolzz)-zolzz+1.) + psih10(I)=-5.-psih_stable(zol0,psimh_stab_opt)+& + (-4.*log(zol10)-zol10+1.) + psih2(I)=-5.-psih_stable(zol0,psimh_stab_opt)+& + (-4.*log(zol2)-zol2+1.) + endif + else + psih(I)=psih_stable(zolzz,psimh_stab_opt) & + -psih_stable(zol0,psimh_stab_opt) + psih10(I)=psih_stable(zol10,psimh_stab_opt) & + -psih_stable(zol0,psimh_stab_opt) + psih2(I)=psih_stable(zol2,psimh_stab_opt) & + -psih_stable(zol0,psimh_stab_opt) + endif +!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ + else + if (zol(I).eq.0) then + psih(I)=0. + psih10(I)=0. + psih2(I)=0. + else +!$$$$$$$$$$$$$$$$$$$$$$$$-Inserted by PN, MS, PS, and SKM-$$$$$$$$$$$$$$$ + if(psimh_unstab_opt == 3.) then + if(zolzz.gt.-0.465) then + psih(I)=psih_unstable(zolzz,psimh_unstab_opt) & + -psih_unstable(zol0,psimh_unstab_opt) + else + psih(I)=psih_unstable(-0.465,psimh_unstab_opt) & + -psih_unstable(zol0,psimh_unstab_opt) & + +log(zolzz/(-0.465))-0.8*((0.465)**(-1./3.)-& + (- zolzz)**(-1./3.)) + endif + else + psih(I)=psih_unstable(zolzz,psimh_unstab_opt) & + -psih_unstable(zol0,psimh_unstab_opt) + endif +! + if(psimh_unstab_opt == 3.) then + if(zol10.gt.-0.465) then + psih10(I)=psih_unstable(zol10,psimh_unstab_opt) & + -psih_unstable(zol0,psimh_unstab_opt) + else + psih10(I)=psih_unstable(-0.465,psimh_unstab_opt)-& + psih_unstable(zol0,psimh_unstab_opt)+ & + log(zol10/(-0.465))-& + 0.8*((0.465)**(-1./3.)-(- zol10)**(-1./3.)) + endif + else + psih10(I)=psih_unstable(zol10,psimh_unstab_opt) & + -psih_unstable(zol0,psimh_unstab_opt) + endif +! + if(psimh_unstab_opt == 3.) then + if(zol2.gt.-0.465) then + psih2(I)=psih_unstable(zol2,psimh_unstab_opt) & + -psih_unstable(zol0,psimh_unstab_opt) + else + psih2(I)=psih_unstable(-0.465,psimh_unstab_opt)-& + psih_unstable(zol0,psimh_unstab_opt) & + +log(zol2/(-0.465))-& + 0.8*((0.465)**(-1./3.)-(- zol2)**(-1./3.)) + endif + else + psih2(I)=psih_unstable(zol2,psimh_unstab_opt) & + -psih_unstable(zol0,psimh_unstab_opt) + endif +!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ + endif + endif +! + PSIQ=ALOG((ZA(I)+z0q)/Z0q)-PSIH(I) + PSIQ2=ALOG((2.+z0q)/Z0q)-PSIH2(I) + PSIQ10=ALOG((10.+z0q)/Z0q)-PSIH10(I) +! PSIQ=GZ1OZ0(I)-PSIH(I)+2.28*SQRT(SQRT(RESTAR))-2. +! PSIQ2=GZ2OZ0(I)-PSIH2(I)+2.28*SQRT(SQRT(RESTAR))-2. +! PSIQ10=GZ10OZ0(I)-PSIH(I)+2.28*SQRT(SQRT(RESTAR))-2. + ENDIF + ENDIF + IF(PRESENT(ck) .and. PRESENT(cd) .and. PRESENT(cka) .and. PRESENT(cda)) THEN + Ck(I)=(karman/psix10)*(karman/psiq10) + Cd(I)=(karman/psix10)*(karman/psix10) + Cka(I)=(karman/psix)*(karman/psiq) + Cda(I)=(karman/psix)*(karman/psix) + ENDIF + IF ( PRESENT(IZ0TLND) ) THEN + IF ( IZ0TLND.GE.1 .AND. (XLAND(I)-1.5).LE.0. ) THEN + ZL=ZNT(I) +! CZIL RELATED CHANGES FOR LAND + VISC=(1.32+0.009*(SCR3(I)-273.15))*1.E-5 + RESTAR=UST(I)*ZL/VISC +! Modify CZIL according to Chen & Zhang, 2009 if iz0tlnd = 1 +! If iz0tlnd = 2, use traditional value + + IF ( IZ0TLND.EQ.1 ) THEN + CZIL = 10.0 ** ( -0.40 * ( ZL / 0.07 ) ) + ELSE IF ( IZ0TLND.EQ.2 ) THEN + CZIL = 0.1 + END IF +! +! ... paj: compute phish for z0t over land +! + z0t=znt(I)/exp(CZIL*KARMAN*SQRT(RESTAR)) +! + zolzz=zol(I)*(za(I)+z0t)/za(I) ! (z+z0t)/L + zol10=zol(I)*(10.+z0t)/za(I) ! (10+z0t)/L + zol2=zol(I)*(2.+z0t)/za(I) ! (2+z0t)/L + zol0=zol(I)*z0t/za(I) ! z0t/L +! + if(zol(I).gt.0.) then +!$$$$$$$$$$$$$$$$$$$$$$$$-Inserted by PN, MS, PS, and SKM-$$$$$$$$$$$$$$$$$$$ + if(psimh_stab_opt == 4.) then + if(zol(I).le.1.) then + psih(I)=psih_stable(zolzz,psimh_stab_opt) & + -psih_stable(zol0,psimh_stab_opt) + psih10(I)=psih_stable(zol10,psimh_stab_opt) & + -psih_stable(zol0,psimh_stab_opt) + psih2(I)=psih_stable(zol2,psimh_stab_opt) & + -psih_stable(zol0,psimh_stab_opt) + else + psih(I)=-5.-psih_stable(zol0,psimh_stab_opt)+ & + (-4.*log(zolzz)-zolzz+1.) + psih10(I)=-5.-psih_stable(zol0,psimh_stab_opt)+& + (-4.*log(zol10)-zol10+1.) + psih2(I)=-5.-psih_stable(zol0,psimh_stab_opt)+& + (-4.*log(zol2)-zol2+1.) + endif + else + psih(I)=psih_stable(zolzz,psimh_stab_opt) & + -psih_stable(zol0,psimh_stab_opt) + psih10(I)=psih_stable(zol10,psimh_stab_opt) & + -psih_stable(zol0,psimh_stab_opt) + psih2(I)=psih_stable(zol2,psimh_stab_opt) & + -psih_stable(zol0,psimh_stab_opt) + endif +!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ + else + if (zol(I).eq.0) then + psih(I)=0. + psih10(I)=0. + psih2(I)=0. + else +!$$$$$$$$$$$$$$$$$$$$$$$$-Inserted by PN, MS, PS, and SKM-$$$$$$$$$$$$$ + if(psimh_unstab_opt == 3.) then + if(zolzz.gt.-0.465) then + psih(I)=psih_unstable(zolzz,psimh_unstab_opt) & + -psih_unstable(zol0,psimh_unstab_opt) + else + psih(I)=psih_unstable(-0.465,psimh_unstab_opt)-& + psih_unstable(zol0,psimh_unstab_opt) & + +log(zolzz/(-0.465))-& + 0.8*((0.465)**(-1./3.)-(- zolzz)**(-1./3.)) + endif + else + psih(I)=psih_unstable(zolzz,psimh_unstab_opt) & + -psih_unstable(zol0,psimh_unstab_opt) + endif +! + if(psimh_unstab_opt == 3.) then + if(zol10.gt.-0.465) then + psih10(I)=psih_unstable(zol10,psimh_unstab_opt) & + -psih_unstable(zol0,psimh_unstab_opt) + else + psih10(I)=psih_unstable(-0.465,psimh_unstab_opt)-& + psih_unstable(zol0,psimh_unstab_opt) & + +log(zol10/(-0.465))-& + 0.8*((0.465)**(-1./3.)-(- zol10)**(-1./3.)) + endif + else + psih10(I)=psih_unstable(zol10,psimh_unstab_opt) & + -psih_unstable(zol0,psimh_unstab_opt) + endif +! + if(psimh_unstab_opt == 3.) then + if(zol2.gt.-0.465) then + psih2(I)=psih_unstable(zol2,psimh_unstab_opt) & + -psih_unstable(zol0,psimh_unstab_opt) + else + psih2(I)=psih_unstable(-0.465,psimh_unstab_opt)-& + psih_unstable(zol0,psimh_unstab_opt) & + +log(zol2/(-0.465))-& + 0.8*((0.465)**(-1./3.)-(- zol2)**(-1./3.)) + endif + else + psih2(I)=psih_unstable(zol2,psimh_unstab_opt) & + -psih_unstable(zol0,psimh_unstab_opt) + endif +!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ + endif + endif +! + PSIQ=ALOG((ZA(I)+z0t)/Z0t)-PSIH(I) + PSIQ2=ALOG((2.+z0t)/Z0t)-PSIH2(I) + PSIT=PSIQ + PSIT2=PSIQ2 +! +! PSIT=GZ1OZ0(I)-PSIH(I)+CZIL*KARMAN*SQRT(RESTAR) +! PSIQ=GZ1OZ0(I)-PSIH(I)+CZIL*KARMAN*SQRT(RESTAR) +! PSIT2=GZ2OZ0(I)-PSIH2(I)+CZIL*KARMAN*SQRT(RESTAR) +! PSIQ2=GZ2OZ0(I)-PSIH2(I)+CZIL*KARMAN*SQRT(RESTAR) + + ENDIF + ENDIF +! TO PREVENT OSCILLATIONS AVERAGE WITH OLD VALUE + UST(I)=0.5*UST(I)+0.5*KARMAN*WSPD(I)/PSIX +! TKE coupling: compute ust without vconv for use in tke scheme + 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 + + U10(I)=UX(I)*PSIX10/PSIX + V10(I)=VX(I)*PSIX10/PSIX + TH2(I)=THGB(I)+DTG*PSIT2/PSIT + Q2(I)=QSFC(I)+(QX(I)-QSFC(I))*PSIQ2/PSIQ + T2(I) = TH2(I)*(PSFCPA(I)/P1000mb)**ROVCP +! + IF((XLAND(I)-1.5).LT.0.)THEN + UST(I)=AMAX1(UST(I),0.001) + ENDIF + MOL(I)=KARMAN*DTG/PSIT/PRT + DENOMQ(I)=PSIQ + DENOMQ2(I)=PSIQ2 + DENOMT2(I)=PSIT2 + FM(I)=PSIX + FH(I)=PSIT +!$$$$$$$$$$-Inserted by PN, MS, PS, and SKM-$$$$$$$$$$$$$ + FM_hat(I) = PSIM_hat + FH_hat(I) = PSIH_hat +!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ + 330 CONTINUE +! + 335 CONTINUE + +!-----COMPUTE THE SURFACE SENSIBLE AND LATENT HEAT FLUXES: + IF ( PRESENT(SCM_FORCE_FLUX) ) THEN + IF (SCM_FORCE_FLUX.EQ.1) GOTO 350 + ENDIF + DO i=its,ite + QFX(i)=0. + HFX(i)=0. + ENDDO + 350 CONTINUE + + IF (ISFFLX.EQ.0) GOTO 410 + +!-----OVER WATER, ALTER ROUGHNESS LENGTH (ZNT) ACCORDING TO WIND (UST). + + DO 360 I=its,ite + IF((XLAND(I)-1.5).GE.0)THEN +! ZNT(I)=CZO*UST(I)*UST(I)/G+OZO +! Since V3.7 (ref: EC Physics document for Cy36r1) + ZNT(I)=CZO*UST(I)*UST(I)/G+0.11*1.5E-5/UST(I) +! V3.9: Add limit as in isftcflx = 1,2 + ZNT(I)=MIN(ZNT(I),2.85e-3) +! COARE 3.5 (Edson et al. 2013) +! CZC = 0.0017*WSPD(I)-0.005 +! CZC = min(CZC,0.028) +! ZNT(I)=CZC*UST(I)*UST(I)/G+0.11*1.5E-5/UST(I) +! AHW: change roughness length, and hence the drag coefficients Ck and Cd + IF ( PRESENT(ISFTCFLX) ) THEN + IF ( ISFTCFLX.NE.0 ) THEN +! ZNT(I)=10.*exp(-9.*UST(I)**(-.3333)) +! ZNT(I)=10.*exp(-9.5*UST(I)**(-.3333)) +! ZNT(I)=ZNT(I) + 0.11*1.5E-5/AMAX1(UST(I),0.01) +! ZNT(I)=0.011*UST(I)*UST(I)/G+OZO +! ZNT(I)=MAX(ZNT(I),3.50e-5) +! AHW 2012: + ZW = MIN((UST(I)/1.06)**(0.3),1.0) + ZN1 = 0.011*UST(I)*UST(I)/G + OZO + ZN2 = 10.*exp(-9.5*UST(I)**(-.3333)) + & + 0.11*1.5E-5/AMAX1(UST(I),0.01) + ZNT(I)=(1.0-ZW) * ZN1 + ZW * ZN2 + ZNT(I)=MIN(ZNT(I),2.85e-3) + ZNT(I)=MAX(ZNT(I),1.27e-7) + ENDIF + ENDIF + ZL = ZNT(I) + ELSE + ZL = 0.01 + ENDIF + FLQC(I)=RHOX(I)*MAVAIL(I)*UST(I)*KARMAN/DENOMQ(I) +! FLQC(I)=RHOX(I)*MAVAIL(I)*UST(I)*KARMAN/( & +! ALOG(KARMAN*UST(I)*ZA(I)/XKA+ZA(I)/ZL)-PSIH(I)) + DTTHX=ABS(THX(I)-THGB(I)) + IF(DTTHX.GT.1.E-5)THEN + FLHC(I)=CPM(I)*RHOX(I)*UST(I)*MOL(I)/(THX(I)-THGB(I)) +! write(*,1001)FLHC(I),CPM(I),RHOX(I),UST(I),MOL(I),THX(I),THGB(I),I + 1001 format(f8.5,2x,f12.7,2x,f12.10,2x,f12.10,2x,f13.10,2x,f12.8,f12.8,2x,i3) + ELSE + FLHC(I)=0. + ENDIF + 360 CONTINUE + +! +!-----COMPUTE SURFACE MOIST FLUX: +! +! IF(IDRY.EQ.1)GOTO 390 +! + IF ( PRESENT(SCM_FORCE_FLUX) ) THEN + IF (SCM_FORCE_FLUX.EQ.1) GOTO 405 + ENDIF + + DO 370 I=its,ite + QFX(I)=FLQC(I)*(QSFC(I)-QX(I)) + QFX(I)=AMAX1(QFX(I),0.) + LH(I)=XLV*QFX(I) + 370 CONTINUE + +!-----COMPUTE SURFACE HEAT FLUX: +! + 390 CONTINUE + DO 400 I=its,ite + IF(XLAND(I)-1.5.GT.0.)THEN + HFX(I)=FLHC(I)*(THGB(I)-THX(I)) +! IF ( PRESENT(ISFTCFLX) ) THEN +! IF ( ISFTCFLX.NE.0 ) THEN +! AHW: add dissipative heating term (commented out in 3.6.1) +! HFX(I)=HFX(I)+RHOX(I)*USTM(I)*USTM(I)*WSPDI(I) +! ENDIF +! ENDIF + ELSEIF(XLAND(I)-1.5.LT.0.)THEN + HFX(I)=FLHC(I)*(THGB(I)-THX(I)) + HFX(I)=AMAX1(HFX(I),-250.) + ENDIF + 400 CONTINUE + + 405 CONTINUE + + DO I=its,ite + IF((XLAND(I)-1.5).GE.0)THEN + ZL=ZNT(I) + ELSE + ZL=0.01 + ENDIF +!v3.1.1 +! CHS(I)=UST(I)*KARMAN/(ALOG(KARMAN*UST(I)*ZA(I) & +! /XKA+ZA(I)/ZL)-PSIH(I)) + CHS(I)=UST(I)*KARMAN/DENOMQ(I) +! GZ2OZ0(I)=ALOG(2./ZNT(I)) +! PSIM2(I)=-10.*GZ2OZ0(I) +! PSIM2(I)=AMAX1(PSIM2(I),-10.) +! PSIH2(I)=PSIM2(I) +! v3.1.1 +! CQS2(I)=UST(I)*KARMAN/(ALOG(KARMAN*UST(I)*2.0 & +! /XKA+2.0/ZL)-PSIH2(I)) +! CHS2(I)=UST(I)*KARMAN/(GZ2OZ0(I)-PSIH2(I)) + CQS2(I)=UST(I)*KARMAN/DENOMQ2(I) + CHS2(I)=UST(I)*KARMAN/DENOMT2(I) + ENDDO + + 410 CONTINUE +!jdf +! DO I=its,ite +! IF(UST(I).GE.0.1) THEN +! RMOL(I)=RMOL(I)*(-FLHC(I))/(UST(I)*UST(I)*UST(I)) +! ELSE +! RMOL(I)=RMOL(I)*(-FLHC(I))/(0.1*0.1*0.1) +! ENDIF +! ENDDO +!jdf + +! + END SUBROUTINE SFCLAYREV1D + +!$$$$$$$$$$$$$$$$$$$$$$$$-Inserted by PN, MS, PS, and SKM-$$$$$$$$$$$$ + SUBROUTINE sfclayrevinit_IITD(psimh_stab_opt,psimh_unstab_opt) + + INTEGER, INTENT(IN) :: psimh_stab_opt,psimh_unstab_opt + INTEGER :: N + REAL :: zolf + + if(psimh_stab_opt == 0) then + DO N=0,1000 +! stable function tables + zolf = float(n)*0.01 + psim_stab(n)=psim_stable_full(zolf) + psih_stab(n)=psih_stable_full(zolf) + ! write(61,*) psimh_stab_opt + ENDDO + + elseif(psimh_stab_opt == 1) then + DO N=0,1000 + zolf = float(n)*0.01 + psim_stab(n)=psim_stable_full_GR07(zolf) + psih_stab(n)=psih_stable_full_GR07(zolf) + ! write(62,*) psimh_stab_opt + ENDDO + + elseif(psimh_stab_opt == 2) then + DO N=0,1000 + zolf = float(n)*0.01 + psim_stab(n)=psim_stable_full_SR20(zolf) + psih_stab(n)=psih_stable_full_SR20(zolf) + ENDDO + + elseif(psimh_stab_opt == 3) then + DO N=0,1000 + zolf = float(n)*0.01 + psim_stab(n)=psim_stable_full_GRY20(zolf) + psih_stab(n)=psih_stable_full_GRY20(zolf) + ENDDO + + elseif(psimh_stab_opt == 4) then + DO N=0,1000 + zolf = float(n)*0.01 + psim_stab(n)=psim_stable_full_H90(zolf) + psih_stab(n)=psih_stable_full_H90(zolf) + ENDDO + + elseif(psimh_stab_opt == 5) then + DO N=0,1000 + zolf = float(n)*0.01 + psim_stab(n)=psim_stable_full_BH91(zolf) + psih_stab(n)=psih_stable_full_BH91(zolf) + ENDDO + + elseif(psimh_stab_opt == 6) then + DO N=0,1000 + zolf = float(n)*0.01 + psim_stab(n)=psim_stable_full_D91(zolf) + psih_stab(n)=psih_stable_full_D91(zolf) + ! write(63,*) psimh_stab_opt + ENDDO + + elseif(psimh_stab_opt .gt. 6) then + CALL wrf_error_fatal("ERROR: psimh_stab_opt more & + than 6 is not allowed") + endif + +! unstable function tables + if(psimh_unstab_opt == 0) then + DO N=0,1000 + zolf = -float(n)*0.01 + psim_unstab(n)=psim_unstable_full(zolf) + psih_unstab(n)=psih_unstable_full(zolf) + ! write(64) psimh_unstab_opt + ENDDO + + elseif(psimh_unstab_opt == 1) then + DO N=0,1000 + zolf = -float(n)*0.01 + psim_unstab(n)=psim_unstable_full_BD71(zolf) + psih_unstab(n)=psih_unstable_full_BD71(zolf) + ! write(65,*) psimh_unstab_opt + ENDDO + + elseif(psimh_unstab_opt == 2) then + DO N=0,1000 + zolf = -float(n)*0.01 + psim_unstab(n)=psim_unstable_full_CL73(zolf) + psih_unstab(n)=psih_unstable_full_CL73(zolf) + ENDDO + + elseif(psimh_unstab_opt == 3) then + DO N=0,1000 + zolf = -float(n)*0.01 + psim_unstab(n)=psim_unstable_full_KY90(zolf) + psih_unstab(n)=psih_unstable_full_KY90(zolf) + ENDDO + + elseif(psimh_unstab_opt .gt. 3) then + CALL wrf_error_fatal("ERROR: psimh_unstab_opt more & + than 3 is not allowed") + endif + + END SUBROUTINE sfclayrevinit_IITD +!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ +! +!$$$$$$$$$$$$$$$$$$$$$$$$-Inserted by PN, MS, PS, and SKM-$$$$$$$$$$$$$$ + REAL function zolri(ri,z,z0,psimh_stab_opt,psimh_unstab_opt) +!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ + IMPLICIT NONE + REAL, INTENT(IN) :: ri, z, z0 + INTEGER, INTENT(IN) :: psimh_stab_opt,psimh_unstab_opt + REAL :: x1, x2, fx1, fx2 + + if (ri.lt.0.)then + x1=-5. + x2=0. + else + x1=0. + x2=5. + endif + + fx1=zolri2(x1,ri,z,z0,psimh_stab_opt,psimh_unstab_opt) + fx2=zolri2(x2,ri,z,z0,psimh_stab_opt,psimh_unstab_opt) + Do While (abs(x1 - x2) > 0.01) +!! check added for potential divide by zero (2019/11) + if(fx1.eq.fx2)return + if(abs(fx2).lt.abs(fx1))then + x1=x1-fx1/(fx2-fx1)*(x2-x1) + fx1=zolri2(x1,ri,z,z0,psimh_stab_opt,psimh_unstab_opt) + zolri=x1 + else + x2=x2-fx2/(fx2-fx1)*(x2-x1) + fx2=zolri2(x2,ri,z,z0,psimh_stab_opt,psimh_unstab_opt) + zolri=x2 + endif +!! + enddo + + return + end function + +!$$$$$$$$$$$$$$$$$$$$$$$$-Inserted by PN, MS, PS, and SKM-$$$$$$$$$$$$$$$$$$$ + REAL function zolri2(zol2,ri2,z,z0,psimh_stab_opt,psimh_unstab_opt) +!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ + IMPLICIT NONE + INTEGER, INTENT(IN) :: psimh_stab_opt, psimh_unstab_opt + REAL, INTENT(IN) :: ri2, z, z0 + REAL, INTENT(INOUT) :: zol2 + REAL :: zol20, zol3, psix2, psih2 + + if(zol2*ri2 .lt. 0.)zol2=0. ! limit zol2 - must be same sign as ri2 +! + zol20=zol2*z0/z ! z0/L + zol3=zol2+zol20 ! (z+z0)/L +!$$$$$$$$$$$$$$$$$$$$$$$-Inserted by PN, MS, PS, and SKM-$$$$$$$$$$$$ + if(ri2.lt.0) then + if(psimh_unstab_opt == 3.) then + if(zol3.ge.-0.465) then + psih2=log((z+z0)/z0)-(psih_unstable(zol3,& + psimh_unstab_opt)-& + psih_unstable(zol20,psimh_unstab_opt)) + psix2=log((z+z0)/z0)-(psim_unstable(zol3, & + psimh_unstab_opt)-& + psim_unstable(zol20,psimh_unstab_opt)) +! write(71,*) zol3 + elseif (zol3.ge.-1.574.and.zol3.lt.-0.465) then +! zol2 = max(zol2,-4.) + zol3 = max(zol3,-100.) + psix2=log((z+z0)/z0)-(psim_unstable(zol3,& + psimh_unstab_opt)-& + psim_unstable(zol20,psimh_unstab_opt)) + psih2=log((z+z0)/z0)-(psih_unstable(-0.465, & + psimh_unstab_opt)& + -psih_unstable(zol20,psimh_unstab_opt)& + +log(zol3/(-0.465))-& + 0.8*((0.465)**(-1./3.)-(- zol3)**(-1./3.))) + elseif (zol3.lt.-1.574) then + zol3 = max(zol3,-100.) + psix2=log((z+z0)/z0)-(psim_unstable(-1.574,& + psimh_unstab_opt)-& + psim_unstable(zol20,psimh_unstab_opt) & + +log(zol3/(-1.574))- & + 1.14*((- zol3)**(1./3.)-(1.574)**(1./3.))) + psih2=log((z+z0)/z0)-(psih_unstable(-0.465, & + psimh_unstab_opt)& + -psih_unstable(zol20,psimh_unstab_opt) & + +log(zol3/(-0.465))-& + 0.8*((0.465)**(-1./3.)-(- zol3)**(-1./3.))) + endif + else + psih2=log((z+z0)/z0)-(psih_unstable(zol3, & + psimh_unstab_opt)-& + psih_unstable(zol20,psimh_unstab_opt)) + psix2=log((z+z0)/z0)-(psim_unstable(zol3, & + psimh_unstab_opt)-& + psim_unstable(zol20,psimh_unstab_opt)) + endif + else + if(psimh_stab_opt == 4.) then + if(zol3.le.1.) then + psix2=log((z+z0)/z0)-(psim_stable(zol3, & + psimh_stab_opt)-psim_stable(zol20,& + psimh_stab_opt)) + psih2=log((z+z0)/z0)-(psih_stable(zol3, & + psimh_stab_opt)-psih_stable(zol20,& + psimh_stab_opt)) + else + psix2=log((z+z0)/z0)+5.+psim_stable(zol20,& + psimh_stab_opt)+(4.*log(zol3)+zol3-1.) + psih2=log((z+z0)/z0)+5.+psih_stable(zol20, & + psimh_stab_opt)+(4.*log(zol3)+zol3-1.) + endif + else + psix2=log((z+z0)/z0)-(psim_stable(zol3, & + psimh_stab_opt)-psim_stable(zol20,& + psimh_stab_opt)) + psih2=log((z+z0)/z0)-(psih_stable(zol3, & + psimh_stab_opt)-psih_stable(zol20, & + psimh_stab_opt)) + endif + endif +!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ + zolri2=zol2*psih2/psix2**2-ri2 +!! + return + end function + +!--------------------------------------------------------------------------------------- +! ... integrated similarity functions ... +! +!$$$$$$$$$$$$$$$$$$$$$$$$-Inserted by PN, MS, PS, and SKM-$$$$$$$$$$$$$$$$$$$ +!==============================-psimh_stab_opt == 0-========================= + REAL function psim_stable_full(zolf) + REAL :: zolf + psim_stable_full=-6.1*log(zolf+(1.+zolf**2.5)**(1./2.5)) + return + end function + + REAL function psih_stable_full(zolf) + REAL :: zolf + psih_stable_full=-5.3*log(zolf+(1.+zolf**1.1)**(1./1.1)) +! return + ! write(51,*) psih_stable_full + return + end function + +!==============================-psimh_stab_opt == 1-======================== + REAL function psim_stable_full_GR07(zolf) + REAL :: zolf, BB_GR, yG, Y_GR + REAL, PARAMETER :: a_GR = 5. + REAL, PARAMETER :: b_GR = 5./6.5 + BB_GR = ((1-b_GR)/b_GR)**(1./3.) + yG = (1.+zolf)**(1./3.) + + Y_GR=2.*ALOG((yG+BB_GR)/(1.+BB_GR))-ALOG((yG**2- & + yG*BB_GR+BB_GR**2)/(1.-BB_GR+BB_GR**2)) + Y_GR=Y_GR+2*sqrt(3.)*(atan((2.*yG-BB_GR)/ & + (sqrt(3.)*BB_GR))-atan((2.-BB_GR)/ & + (sqrt(3.)*BB_GR))) + Y_GR = -3.*(yG - 1.) + 0.5*BB_GR*Y_GR + psim_stable_full_GR07 = (a_GR/b_GR)*(Y_GR) + return + end function + + REAL function psih_stable_full_GR07(zolf) + REAL :: zolf, BB_GR, yG, Y_GRh + REAL, PARAMETER :: a_GR = 5. + REAL, PARAMETER :: b_GR = 5. + REAL, PARAMETER :: c_GR = 3. + BB_GR = sqrt(c_GR**2 - 4.) + yG = zolf + Y_GRh=log((2*yG+c_GR-BB_GR)/(2*yG+c_GR+BB_GR))-& + log((c_GR-BB_GR)/(c_GR+BB_GR)) + Y_GRh=-0.5*b_GR*log(1.+c_GR*yG+yG**2)+(-a_GR/BB_GR+& + 0.5*b_GR*c_GR/BB_GR)*Y_GRh + psih_stable_full_GR07 = Y_GRh + if(psimh_stab_opt == 1.) then + ! write(52,*) psih_stable_full_GR07 + endif + return + end function + +!==========================-psimh_stab_opt == 2-========================== + + REAL function psim_stable_full_SR20(zolf) + REAL :: zolf, BB_GR, yG, Y_GR + REAL, PARAMETER :: a_GR = 5.4 + REAL, PARAMETER :: b_GR = 0.89 + BB_GR = ((1-b_GR)/b_GR)**(1./3.) + yG = (1.+zolf)**(1./3.) +! yG = zolf + Y_GR=2.*ALOG((yG+BB_GR)/(1.+BB_GR))-& + ALOG((yG**2-yG*BB_GR+& + BB_GR**2)/(1.-BB_GR+BB_GR**2)) + Y_GR=Y_GR+2*sqrt(3.)*(atan((2.*yG-BB_GR) & + /(sqrt(3.)*BB_GR))-atan((2.-BB_GR)/ & + (sqrt(3.)*BB_GR))) + Y_GR = -3.*(yG - 1.) + 0.5*BB_GR*Y_GR + psim_stable_full_SR20 = (a_GR/b_GR)*(Y_GR) + return + end function + + REAL function psih_stable_full_SR20(zolf) + REAL :: zolf, BB_GR, yG, Y_GR + REAL, PARAMETER :: a_GR = 5. + REAL, PARAMETER :: b_GR = 5. + REAL, PARAMETER :: c_GR = 3. + BB_GR = sqrt(c_GR**2 - 4.) + yG = zolf + Y_GR=log((2*yG+c_GR-BB_GR)/(2*yG+c_GR+BB_GR))-& + log((c_GR-BB_GR)/(c_GR+BB_GR)) + Y_GR=-0.5*b_GR*log(1.+c_GR*yG+yG**2)+(-a_GR/BB_GR+& + 0.5*b_GR*c_GR/BB_GR)*Y_GR + psih_stable_full_SR20 = Y_GR + if(psimh_stab_opt == 2.) then + ! write(53,*) psih_full_stable_SR20 + endif + return + end function + +!==========================-psimh_stab_opt == 3-========================= + + REAL function psim_stable_full_GRY20(zolf) + REAL :: zolf +! REAL, PARAMETER :: a_grm = 5. +! REAL, PARAMETER :: b_grm = 0.3 + x_gry = (1.+0.3*zolf)**(1./3.) + psim_stable_full_GRY20 = -3.*(5./0.3)*(x_gry - 1.) + return + end function + + REAL function psih_stable_full_GRY20(zolf) + REAL :: zolf +! REAL, PARAMETER :: a_grh = 5. +! REAL, PARAMETER :: b_grh = 0.4 +! REAL, PARAMETER :: pr_grh = 0.98 + psih_stable_full_GRY20 = -0.98*(5./0.4)*log(1+0.4*zolf) +! if(psimh_opt == 3.) then +! write(74,*) GRY20 +! endif + return + end function + +!========================-psimh_stab_opt == 4-============================ + + REAL function psim_stable_full_H90(zolf) + REAL :: zolf + psim_stable_full_H90 = -5.*zolf + return + end function + + REAL function psih_stable_full_H90(zolf) + REAL :: zolf + psih_stable_full_H90 = -5.*zolf + ! write(54,*) psih_stable_full_H90 + return + end function + +!========================-psimh_stab_opt == 5-============================= + + REAL function psim_stable_full_BH91(zolf) + REAL :: zolf + REAL, PARAMETER :: a=1. , b=0.666, c=5. , d=0.35 + psim_stable_full_BH91 = -(a*zolf + b*(zolf-(c/d))*& + exp(-d*zolf) +(b*c/d)) + return + end function + + REAL function psih_stable_full_BH91(zolf) + REAL :: zolf + REAL, PARAMETER :: a=1. , b=0.666, c=5. , d=0.35 + psih_stable_full_BH91 = -((1.+0.666*a*zolf)**(1.5)+& + b*(zolf - (c/d))*exp(-d*zolf) + (b*c/d)-1.) + ! write(756,*) psih_stable_full_BH91 + return + end function + +!=======================-psimh_stab_opt == 6-================================ + + REAL function psim_stable_full_D91(zolf) + REAL :: zolf + psim_stable_full_D91= 1.-(1.+6.25*zolf)**(4./5.) +!-6.1*log(zolf+(1+zolf**2.5)**(1./2.5)) + return + end function +! + REAL function psih_stable_full_D91(zolf) + REAL :: zolf + psih_stable_full_D91= 1.-(1.+6.25*zolf)**(4./5.) + !-5.3*log(zolf+(1+zolf**1.1)**(1./1.1)) + if(psimh_stab_opt == 6.) then + ! write(55,*) psih_stable_full_D91 + endif + return + end function + +!Integrated similarity functions for unstable conditions +!==============================psimh_unstab_opt == 0-====================== + REAL function psim_unstable_full(zolf) + REAL :: zolf, x, psimk, ym, psimc + x=(1.-16.*zolf)**.25 + psimk=2*ALOG(0.5*(1+X))+ALOG(0.5*(1+X*X))-2.*ATAN(X)+ & + 2.*ATAN(1.) + ym=(1.-10.*zolf)**0.33 + psimc=(3./2.)*log((ym**2.+ym+1.)/3.)-sqrt(3.)* & + ATAN((2.*ym+1)/sqrt(3.))+4.*ATAN(1.)/sqrt(3.) + psim_unstable_full=(psimk+zolf**2*(psimc))/(1+zolf**2.) +! return + ! write(79,*) psim_unstable_full + return + end function + + REAL function psih_unstable_full(zolf) + REAL :: y, zolf, psihk, yh, psihc + y=(1.-16.*zolf)**.5 + psihk=2.*log((1+y)/2.) + + yh=(1.-34.*zolf)**0.33 + psihc=(3./2.)*log((yh**2.+yh+1.)/3.)-sqrt(3.)* & + ATAN((2.*yh+1)/sqrt(3.))+4.*ATAN(1.)/sqrt(3.) +! + psih_unstable_full=(psihk+zolf**2*(psihc))/(1+zolf**2.) + return + end function + +!=============================psimh_unstab_opt == 1-======================== + + REAL function psim_unstable_full_BD71(zolf) + REAL :: zolf, x, psimk + x=(1.-16.*zolf)**.25 + psimk=2*ALOG(0.5*(1+X))+ALOG(0.5*(1+X*X))-2.*ATAN(X)+ & + 2.*ATAN(1.) + psim_unstable_full_BD71=psimk + ! write(72,*) psim_unstable_full_BD71 + return + end function +! + REAL function psih_unstable_full_BD71(zolf) + REAL :: y, zolf, psihk + y=(1.-16.*zolf)**.5 + psihk=2.*log((1+y)/2.) +! + psih_unstable_full_BD71=psihk + return + end function + +!===============================-psimh_unstab_opt ==2-====================== + + REAL function psim_unstable_full_CL73(zolf) + REAL :: zolf, ym, psimc + ym=(1.-10.*zolf)**0.33 + psimc=(3./2.)*log((ym**2.+ym+1.)/3.)-sqrt(3.)* & + ATAN((2.*ym+1)/sqrt(3.))+4.*ATAN(1.)/sqrt(3.) + psim_unstable_full_CL73=psimc + ! write(73,*) psim_unstable_full_CL73 + return + end function +! + REAL function psih_unstable_full_CL73(zolf) + REAL :: zolf, yh, psihc + yh=(1.-34.*zolf)**0.33 + psihc=(3./2.)*log((yh**2.+yh+1.)/3.)-sqrt(3.)* & + ATAN((2.*yh+1)/sqrt(3.))+4.*ATAN(1.)/sqrt(3.) +! + psih_unstable_full_CL73=psihc + return + end function + +!================================-psimh_unstab_opt ==3-===================== + REAL function psim_unstable_full_KY90(zolf) + REAL :: x11, term1, term2, term3, zolf + x11 = (1.-16.*zolf)**(0.25) + term1 = 2.*log((1.+x11)*0.5) + term2 = log((1.+x11*x11)*0.5) + term3 = 2.*ATAN(x11) + psim_unstable_full_KY90 = term1+term2-term3+ATAN(1.)*0.5 + ! write(74,*) psim_unstable_full_KY90 + return + end function + + REAL function psih_unstable_full_KY90(zolf) + REAL :: y11, zolf + y11 = sqrt(1.-16.*zolf) + psih_unstable_full_KY90 = 2.*log((1.+y11)*0.5) + return + end function +!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ + +! look-up table functions +!$$$$$$$$$$$$$$$$$$$$$$$$-Inserted by PN, MS, PS, and SKM-$$$$$$$$$$$$$$$$$$$ + REAL function psim_stable(zolf,psimh_stab_opt) + INTEGER,INTENT(IN) :: psimh_stab_opt + integer :: nzol + real :: rzol + nzol = int(zolf*100.) + rzol = zolf*100. - nzol + if(nzol+1 .lt. 1000)then + psim_stable = psim_stab(nzol) + & + rzol*(psim_stab(nzol+1)-psim_stab(nzol)) + else + if(psimh_stab_opt == 0) then + psim_stable = psim_stable_full(zolf) + elseif(psimh_stab_opt == 1) then + psim_stable = psim_stable_full_GR07(zolf) + elseif(psimh_stab_opt == 2) then + psim_stable = psim_stable_full_SR20(zolf) + elseif(psimh_stab_opt == 3) then + psim_stable = psim_stable_full_GRY20(zolf) + elseif(psimh_stab_opt == 4) then + psim_stable = psim_stable_full_H90(zolf) + elseif(psimh_stab_opt == 5) then + psim_stable = psim_stable_full_BH91(zolf) + elseif(psimh_stab_opt == 6) then + psim_stable = psim_stable_full_D91(zolf) + endif + endif + return + end function + + REAL function psih_stable(zolf,psimh_stab_opt) + INTEGER,INTENT(IN) :: psimh_stab_opt + integer :: nzol + real :: rzol + nzol = int(zolf*100.) + rzol = zolf*100. - nzol + if(nzol+1 .lt. 1000)then + psih_stable = psih_stab(nzol) + & + rzol*(psih_stab(nzol+1)-psih_stab(nzol)) + else + if(psimh_stab_opt == 0) then + psih_stable = psih_stable_full(zolf) + elseif(psimh_stab_opt == 1) then + psih_stable = psih_stable_full_GR07(zolf) + elseif(psimh_stab_opt == 2) then + psih_stable = psih_stable_full_SR20(zolf) + elseif(psimh_stab_opt == 3) then + psih_stable = psih_stable_full_GRY20(zolf) + elseif(psimh_stab_opt == 4) then + psih_stable = psih_stable_full_H90(zolf) + elseif(psimh_stab_opt == 5) then + psih_stable = psih_stable_full_BH91(zolf) + elseif(psimh_stab_opt == 6) then + psih_stable = psih_stable_full_D91(zolf) + endif + endif + return + end function + + REAL function psim_unstable(zolf,psimh_unstab_opt) + INTEGER,INTENT(IN) :: psimh_unstab_opt + integer :: nzol + real :: rzol + nzol = int(-zolf*100.) + rzol = -zolf*100. - nzol + if(nzol+1 .lt. 1000)then + psim_unstable = psim_unstab(nzol) + & + rzol*(psim_unstab(nzol+1)-psim_unstab(nzol)) + else + if(psimh_unstab_opt == 0) then + psim_unstable = psim_unstable_full(zolf) + ! write(66,*) psimh_unstab_opt + elseif(psimh_unstab_opt == 1) then + psim_unstable = psim_unstable_full_BD71(zolf) + elseif(psimh_unstab_opt == 2) then + psim_unstable = psim_unstable_full_CL73(zolf) + elseif(psimh_unstab_opt == 3) then + psim_unstable = psim_unstable_full_KY90(zolf) + endif + endif + return + end function + + REAL function psih_unstable(zolf,psimh_unstab_opt) + INTEGER,INTENT(IN) :: psimh_unstab_opt + integer :: nzol + real :: rzol + nzol = int(-zolf*100.) + rzol = -zolf*100. - nzol + if(nzol+1 .lt. 1000)then + psih_unstable = psih_unstab(nzol) + & + rzol*(psih_unstab(nzol+1)-psih_unstab(nzol)) + else + if(psimh_unstab_opt == 0) then + psih_unstable = psih_unstable_full(zolf) + ! write(67,*) psimh_unstab_opt + elseif(psimh_unstab_opt == 1) then + psih_unstable = psih_unstable_full_BD71(zolf) + elseif(psimh_unstab_opt == 2) then + psih_unstable = psih_unstable_full_CL73(zolf) + elseif(psimh_unstab_opt == 3) then + psih_unstable = psih_unstable_full_KY90(zolf) + endif + endif + return + end function + +END MODULE module_sf_sfclayrev_IITD +!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ diff --git a/phys/module_surface_driver.F b/phys/module_surface_driver.F index 22f96580fb..a2afcc967d 100644 --- a/phys/module_surface_driver.F +++ b/phys/module_surface_driver.F @@ -268,6 +268,9 @@ SUBROUTINE surface_driver( & & ,tml,t0ml,hml,h0ml,huml,hvml,f,tmoml & & ,oml_relaxation_time & & ,ustm,ck,cka,cd,cda,isftcflx,iz0tlnd & +!$$$$$$$$$$$$$$$$$$$$$-Inserted by PN, MS, PS, and SKM-$$$$$$$$$$$$$$$$ + & ,psimh_stab_opt, psimh_unstab_opt & +!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ & ,isurban, mminlu & & ,snotime & & ,rdlai2d & @@ -346,10 +349,15 @@ SUBROUTINE surface_driver( & & ,irrigation,sf_surf_irr_scheme, irr_daily_amount & !IRRIG & ,irr_start_hour,irr_num_hours,irr_start_julianday & & ,irr_end_julianday,irr_freq,irr_ph,irr_rand_field & - & ) +!$$$$$$$$$$$$$$$$$$$$$$-Inserted by PN, MS, PS, and SKM-$$$$$$$$$$$$$$$ + & ,FM_hat, FH_hat ) +!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ USE module_state_description, ONLY : SFCLAYSCHEME & ,SFCLAYREVSCHEME & +!$$$$$$$$$$$$$$$$$$$$$$$$$$-Inserted by PN, MS, PS, and SKM-$$$$$$$$$$$$$ + ,SFCLAYREV_IITD_SCHEME & +!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ ,MYJSFCSCHEME & ,QNSESFCSCHEME & ,MYJPBLSCHEME & @@ -400,6 +408,9 @@ SUBROUTINE surface_driver( & USE module_sf_pxlsm USE module_sf_temfsfclay USE module_sf_sfclayrev +!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$-Inserted by PN, MS, PS, and SKM-$$$$$$$$$ + USE module_sf_sfclayrev_IITD +!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ USE module_sf_noah_seaice_drv #if ( EM_CORE==1) USE module_sf_mynnsfc_driver @@ -848,6 +859,10 @@ SUBROUTINE surface_driver( & REAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT):: PSHLTR REAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT):: FHH REAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT):: FM +!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$-Inserted by PN, MS, PS, and SKM-$$$$$$$$$ + REAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT):: FM_hat + REAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT):: FH_hat +!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ REAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT):: PSIH REAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT):: PSIM REAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT):: Q10 @@ -1099,6 +1114,9 @@ SUBROUTINE surface_driver( & REAL, DIMENSION( ims:ime, jms:jme ), OPTIONAL, INTENT(IN):: diffuse_frac INTEGER, OPTIONAL, INTENT(IN ):: ISFTCFLX,IZ0TLND +!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$-Inserted by PN, MS, PS, and SKM-$$$$$$$$$ + INTEGER, INTENT(IN) :: psimh_stab_opt, psimh_unstab_opt +!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ INTEGER, OPTIONAL, INTENT(IN ):: SF_OCEAN_PHYSICS REAL , OPTIONAL, INTENT(IN ):: OML_HML0 REAL , OPTIONAL, INTENT(IN ):: OML_GAMMA @@ -1967,6 +1985,9 @@ SUBROUTINE surface_driver( & isisfc = ( FRACTIONAL_SEAICE .EQ. 1 .AND. ( & (sf_sfclay_physics .EQ. SFCLAYSCHEME ) .OR. & (sf_sfclay_physics .EQ. SFCLAYREVSCHEME ) .OR. & +!$$$$$$$$$$$$$$$$$$$$$$$$$$$$-Inserted by PN, MS, PS, and SKM-$$$$$$$$$$ + (sf_sfclay_physics .EQ. SFCLAYREV_IITD_SCHEME) .OR. & +!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ (sf_sfclay_physics .EQ. PXSFCSCHEME ) .OR. & (sf_sfclay_physics .EQ. MYJSFCSCHEME ) .OR. & (sf_sfclay_physics .EQ. QNSESFCSCHEME ) .OR. & !emt @@ -2204,6 +2225,64 @@ SUBROUTINE surface_driver( & CALL wrf_error_fatal('Lacking arguments for SFCLAY in surface driver') ENDIF +!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$-Inserted by PN, MS, PS, and SKM-$$$$$$$ + CASE (SFCLAYREV_IITD_SCHEME) +! DX varies spatially in NMM, therefore, SFCLAY cannot be called +! because it takes a scalar DX. NMM passes in a dummy value for this +! scalar. NEEDS FURTHER ATTENTION. JM 20050215 + IF (PRESENT(qv_curr) .AND. & + PRESENT(mol) .AND. PRESENT(regime) .AND. & + .TRUE. ) THEN + CALL wrf_debug( 100, 'in SFCLAY' ) + IF ( FRACTIONAL_SEAICE == 1 ) THEN + CALL SFCLAYREV_SEAICE_WRAPPER(u_phytmp,v_phytmp,t_phy,qv_curr,& + p_phy,dz8w,cp,g,rcp,r_d,xlv,psfc,chs,chs2,cqs2,cpm, & + znt,ust,pblh,mavail,zol,mol,regime,psim,psih,fm,fhh, & + xland,hfx,qfx,lh,tsk,flhc,flqc,qgh,qsfc,rmol, & + u10,v10,th2,t2,q2, & + gz1oz0,wspd,br,isfflx,dx2d, & + svp1,svp2,svp3,svpt0,ep_1,ep_2,karman, & + P1000mb,lakemask, & + kim_tofd, tofd_factor, var2d, & + XICE,SST,TSK_SEA, & + CHS2_SEA,CHS_SEA,CPM_SEA,CQS2_SEA,FLHC_SEA,FLQC_SEA, & + HFX_SEA,LH_SEA,QFX_SEA,QGH_SEA,QSFC_SEA,ZNT_SEA, & + ITIMESTEP,TICE2TSK_IF2COLD,XICE_THRESHOLD, & + ids,ide, jds,jde, kds,kde, & + ims,ime, jms,jme, kms,kme, & + i_start(ij),i_end(ij), j_start(ij),j_end(ij), kts,kte, & + ustm,ck,cka,cd,cda,isftcflx,iz0tlnd, & + shalwater_z0,water_depth, & + scm_force_flux,sf_surface_physics,errmsg,errflg ) + + ELSE + CALL SFCLAYREV_IITD(u_phytmp,v_phytmp,t_phy,qv_curr,& + p_phy,dz8w,cp,g,rcp,r_d,xlv,psfc,chs,chs2,cqs2,cpm, & + znt,ust,pblh,mavail,zol,mol,regime,psim,psih,fm,fhh, & + fm_hat,fh_hat,xland,hfx,qfx,lh,tsk,flhc,flqc,qgh,qsfc,rmol, & + u10,v10,th2,t2,q2, & + gz1oz0,wspd,br,isfflx,dx, & + svp1,svp2,svp3,svpt0,ep_1,ep_2,karman,eomeg,stbolt, & + P1000mb, & + ids,ide, jds,jde, kds,kde, & + ims,ime, jms,jme, kms,kme, & + i_start(ij),i_end(ij), j_start(ij),j_end(ij), kts,kte, & + ustm,ck,cka,cd,cda,isftcflx,iz0tlnd, psimh_stab_opt, & + psimh_unstab_opt, scm_force_flux ) +#if ( EM_CORE==1) + DO j = j_start(ij),j_end(ij) + DO i = i_start(ij),i_end(ij) + ch(i,j) = chs (i,j) +!! ch(i,j) = flhc(i,j)/( cpm(i,j)*rho(i,kts,j) ) + end do + end do +#endif + ENDIF + ELSE + CALL wrf_error_fatal('Lacking arguments for SFCLAY in surface driver') + ENDIF +!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ + CASE (PXSFCSCHEME) IF (PRESENT(qv_curr) .AND. & PRESENT(mol) .AND. PRESENT(regime) .AND. & diff --git a/share/module_check_a_mundo.F b/share/module_check_a_mundo.F index c35d74d885..8a2f95614b 100644 --- a/share/module_check_a_mundo.F +++ b/share/module_check_a_mundo.F @@ -393,6 +393,10 @@ END FUNCTION bep_bem_ngr_u (model_config_rec % sf_sfclay_physics(i) .NE. nosfcscheme .AND. & model_config_rec % sf_sfclay_physics(i) .NE. sfclayscheme .AND. & model_config_rec % sf_sfclay_physics(i) .NE. sfclayrevscheme .AND. & +!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$-Inserted by PN, MS, PS, and SKM-$$$$$$$$$$$$$$ + model_config_rec % sf_sfclay_physics(i) .NE. & + sfclayrev_IITD_scheme .AND. & +!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ model_config_rec % sf_sfclay_physics(i) .NE. mynnsfcscheme ) ) THEN wrf_err_message = '--- ERROR: SMS-3DTKE scheme works with sf_sfclay_physics = 0,1,5,91 ' CALL wrf_message ( wrf_err_message )