diff --git a/CMakeLists.txt b/CMakeLists.txt index b8d3c3e18..8e6785c71 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -191,15 +191,17 @@ elseif (${CMAKE_Fortran_COMPILER_ID} STREQUAL "Intel") ${CMAKE_CURRENT_SOURCE_DIR}/physics/cu_gf_sh.F90 ${CMAKE_CURRENT_SOURCE_DIR}/physics/module_bl_mynn.F90 ${CMAKE_CURRENT_SOURCE_DIR}/physics/module_MYNNPBL_wrapper.F90 - ${CMAKE_CURRENT_SOURCE_DIR}/physics/module_sf_mynn.F90 ${CMAKE_CURRENT_SOURCE_DIR}/physics/module_MYNNSFC_wrapper.F90 - ${CMAKE_CURRENT_SOURCE_DIR}/physics/module_MYNNrad_pre.F90 - ${CMAKE_CURRENT_SOURCE_DIR}/physics/module_MYNNrad_post.F90 ${CMAKE_CURRENT_SOURCE_DIR}/physics/module_mp_thompson_make_number_concentrations.F90 ${CMAKE_CURRENT_SOURCE_DIR}/physics/module_SF_JSFC.F90 ${CMAKE_CURRENT_SOURCE_DIR}/physics/module_BL_MYJPBL.F90 PROPERTIES COMPILE_FLAGS "-r8 -ftz") + # Reduce optimization for module_sf_mynn.F90 (to avoid an apparent compiler bug with Intel 18 on Hera) + SET_SOURCE_FILES_PROPERTIES(${CMAKE_CURRENT_SOURCE_DIR}/physics/module_sf_mynn.F90 + PROPERTIES COMPILE_FLAGS "${CMAKE_Fortran_FLAGS_OPT} -O1") + list(APPEND SCHEMES_SFX_OPT ${CMAKE_CURRENT_SOURCE_DIR}/physics/module_sf_mynn.F90) + # Replace -xHost or -xCORE-AVX2 with -xCORE-AVX-I for certain files set(CMAKE_Fortran_FLAGS_LOPT1 ${CMAKE_Fortran_FLAGS_OPT}) string(REPLACE "-xHOST" "-xCORE-AVX-I" diff --git a/physics/GFS_debug.F90 b/physics/GFS_debug.F90 index df56cc069..3bb50d9ef 100644 --- a/physics/GFS_debug.F90 +++ b/physics/GFS_debug.F90 @@ -225,6 +225,7 @@ subroutine GFS_diagtoscreen_run (Model, Statein, Stateout, Sfcprop, Coupling, call print_var(mpirank,omprank, blkno, 'Tbd%acv' , Tbd%acv) call print_var(mpirank,omprank, blkno, 'Tbd%acvb' , Tbd%acvb) call print_var(mpirank,omprank, blkno, 'Tbd%acvt' , Tbd%acvt) + call print_var(mpirank,omprank, blkno, 'Tbd%hpbl' , Tbd%hpbl) if (Model%do_sppt) then call print_var(mpirank,omprank, blkno, 'Tbd%dtdtr' , Tbd%dtdtr) call print_var(mpirank,omprank, blkno, 'Tbd%dtotprcp' , Tbd%dtotprcp) @@ -294,7 +295,6 @@ subroutine GFS_diagtoscreen_run (Model, Statein, Stateout, Sfcprop, Coupling, call print_var(mpirank,omprank, blkno, 'Diag%dpt2m ', Diag%dpt2m) call print_var(mpirank,omprank, blkno, 'Diag%zlvl ', Diag%zlvl) call print_var(mpirank,omprank, blkno, 'Diag%psurf ', Diag%psurf) - call print_var(mpirank,omprank, blkno, 'Diag%hpbl ', Diag%hpbl) call print_var(mpirank,omprank, blkno, 'Diag%pwat ', Diag%pwat) call print_var(mpirank,omprank, blkno, 'Diag%t1 ', Diag%t1) call print_var(mpirank,omprank, blkno, 'Diag%q1 ', Diag%q1) @@ -310,6 +310,7 @@ subroutine GFS_diagtoscreen_run (Model, Statein, Stateout, Sfcprop, Coupling, call print_var(mpirank,omprank, blkno, 'Diag%tdomzr ', Diag%tdomzr) call print_var(mpirank,omprank, blkno, 'Diag%tdomip ', Diag%tdomip) call print_var(mpirank,omprank, blkno, 'Diag%tdoms ', Diag%tdoms) + ! CCPP/RUC only if (Model%lsm == Model%lsm_ruc) then call print_var(mpirank,omprank, blkno, 'Diag%wet1 ', Sfcprop%wetness) else @@ -345,6 +346,7 @@ subroutine GFS_diagtoscreen_run (Model, Statein, Stateout, Sfcprop, Coupling, if(Model%lradar) then call print_var(mpirank,omprank, blkno, 'Diag%refl_10cm ', Diag%refl_10cm) end if + ! CCPP/MYNNPBL only if (Model%do_mynnedmf) then call print_var(mpirank,omprank, blkno, 'Diag%edmf_a ', Diag%edmf_a) call print_var(mpirank,omprank, blkno, 'Diag%edmf_w ', Diag%edmf_w) diff --git a/physics/GFS_surface_composites.F90 b/physics/GFS_surface_composites.F90 index 9636eb384..0060e1a7b 100644 --- a/physics/GFS_surface_composites.F90 +++ b/physics/GFS_surface_composites.F90 @@ -28,10 +28,11 @@ end subroutine GFS_surface_composites_pre_finalize subroutine GFS_surface_composites_pre_run (im, frac_grid, flag_cice, cplflx, landfrac, lakefrac, oceanfrac, & frland, dry, icy, lake, ocean, wet, cice, cimin, zorl, zorlo, zorll, zorl_ocn, & zorl_lnd, zorl_ice, snowd, snowd_ocn, snowd_lnd, snowd_ice, tprcp, tprcp_ocn, & - tprcp_lnd, tprcp_ice, uustar, uustar_lnd, uustar_ice, weasd, weasd_ocn, & - weasd_lnd, weasd_ice, ep1d_ice, tsfc, tsfco, tsfcl, tsfc_ocn, tsfc_lnd, & - tsfc_ice, tisfc, tice, tsurf, tsurf_ocn, tsurf_lnd, tsurf_ice, gflx_ice, & - tgice, islmsk, semis_rad, semis_ocn, semis_lnd, semis_ice, & + tprcp_lnd, tprcp_ice, uustar, uustar_ocn, uustar_lnd, uustar_ice, & + weasd, weasd_ocn, weasd_lnd, weasd_ice, ep1d_ice, tsfc, tsfco, tsfcl, tsfc_ocn,& + tsfc_lnd, tsfc_ice, tisfc, tice, tsurf, tsurf_ocn, tsurf_lnd, tsurf_ice, & + gflx_ice, tgice, islmsk, semis_rad, semis_ocn, semis_lnd, semis_ice, & + qss, qss_ocn, qss_lnd, qss_ice, hflx, hflx_ocn, hflx_lnd, hflx_ice, & min_lakeice, min_seaice, errmsg, errflg) implicit none @@ -45,12 +46,13 @@ subroutine GFS_surface_composites_pre_run (im, frac_grid, flag_cice, cplflx, lan real(kind=kind_phys), dimension(im), intent(in ) :: landfrac, lakefrac, oceanfrac real(kind=kind_phys), dimension(im), intent(inout) :: cice real(kind=kind_phys), dimension(im), intent( out) :: frland - real(kind=kind_phys), dimension(im), intent(in ) :: zorl, snowd, tprcp, uustar, weasd + real(kind=kind_phys), dimension(im), intent(in ) :: zorl, snowd, tprcp, uustar, weasd, qss, hflx real(kind=kind_phys), dimension(im), intent(inout) :: zorlo, zorll, tsfc, tsfco, tsfcl, tisfc, tsurf real(kind=kind_phys), dimension(im), intent(inout) :: snowd_ocn, snowd_lnd, snowd_ice, tprcp_ocn, & tprcp_lnd, tprcp_ice, zorl_ocn, zorl_lnd, zorl_ice, tsfc_ocn, tsfc_lnd, tsfc_ice, tsurf_ocn, & - tsurf_lnd, tsurf_ice, uustar_lnd, uustar_ice, weasd_ocn, weasd_lnd, weasd_ice, ep1d_ice, gflx_ice + tsurf_lnd, tsurf_ice, uustar_ocn, uustar_lnd, uustar_ice, weasd_ocn, weasd_lnd, weasd_ice, & + qss_ocn, qss_lnd, qss_ice, hflx_ocn, hflx_lnd, hflx_ice, ep1d_ice, gflx_ice real(kind=kind_phys), dimension(im), intent( out) :: tice real(kind=kind_phys), intent(in ) :: tgice integer, dimension(im), intent(in ) :: islmsk @@ -145,6 +147,7 @@ subroutine GFS_surface_composites_pre_run (im, frac_grid, flag_cice, cplflx, lan tprcp_lnd(i) = tprcp(i) tprcp_ice(i) = tprcp(i) if (wet(i)) then ! Water + uustar_ocn(i) = uustar(i) zorl_ocn(i) = zorlo(i) tsfc_ocn(i) = tsfco(i) tsurf_ocn(i) = tsfco(i) @@ -153,6 +156,8 @@ subroutine GFS_surface_composites_pre_run (im, frac_grid, flag_cice, cplflx, lan weasd_ocn(i) = zero snowd_ocn(i) = zero semis_ocn(i) = 0.984d0 + qss_ocn(i) = qss(i) + hflx_ocn(i) = hflx(i) endif if (dry(i)) then ! Land uustar_lnd(i) = uustar(i) @@ -162,6 +167,8 @@ subroutine GFS_surface_composites_pre_run (im, frac_grid, flag_cice, cplflx, lan tsurf_lnd(i) = tsfcl(i) snowd_lnd(i) = snowd(i) semis_lnd(i) = semis_rad(i) + qss_lnd(i) = qss(i) + hflx_lnd(i) = hflx(i) end if if (icy(i)) then ! Ice uustar_ice(i) = uustar(i) @@ -173,6 +180,8 @@ subroutine GFS_surface_composites_pre_run (im, frac_grid, flag_cice, cplflx, lan ep1d_ice(i) = zero gflx_ice(i) = zero semis_ice(i) = 0.95d0 + qss_ice(i) = qss(i) + hflx_ice(i) = hflx(i) end if enddo diff --git a/physics/GFS_surface_composites.meta b/physics/GFS_surface_composites.meta index 74c6b9575..bf613e160 100644 --- a/physics/GFS_surface_composites.meta +++ b/physics/GFS_surface_composites.meta @@ -262,6 +262,15 @@ kind = kind_phys intent = in optional = F +[uustar_ocn] + standard_name = surface_friction_velocity_over_ocean + long_name = surface friction velocity over ocean + units = m s-1 + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F [uustar_lnd] standard_name = surface_friction_velocity_over_land long_name = surface friction velocity over land @@ -495,6 +504,78 @@ kind = kind_phys intent = inout optional = F +[qss] + standard_name = surface_specific_humidity + long_name = surface air saturation specific humidity + units = kg kg-1 + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = in + optional = F +[qss_ocn] + standard_name = surface_specific_humidity_over_ocean + long_name = surface air saturation specific humidity over ocean + units = kg kg-1 + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[qss_lnd] + standard_name = surface_specific_humidity_over_land + long_name = surface air saturation specific humidity over land + units = kg kg-1 + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[qss_ice] + standard_name = surface_specific_humidity_over_ice + long_name = surface air saturation specific humidity over ice + units = kg kg-1 + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[hflx] + standard_name = kinematic_surface_upward_sensible_heat_flux + long_name = kinematic surface upward sensible heat flux + units = K m s-1 + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = in + optional = F +[hflx_ocn] + standard_name = kinematic_surface_upward_sensible_heat_flux_over_ocean + long_name = kinematic surface upward sensible heat flux over ocean + units = K m s-1 + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[hflx_lnd] + standard_name = kinematic_surface_upward_sensible_heat_flux_over_land + long_name = kinematic surface upward sensible heat flux over land + units = K m s-1 + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[hflx_ice] + standard_name = kinematic_surface_upward_sensible_heat_flux_over_ice + long_name = kinematic surface upward sensible heat flux over ice + units = K m s-1 + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F [min_lakeice] standard_name = lake_ice_minimum long_name = minimum lake ice value diff --git a/physics/module_MYNNSFC_wrapper.F90 b/physics/module_MYNNSFC_wrapper.F90 index 5471c4825..82cdbca76 100644 --- a/physics/module_MYNNSFC_wrapper.F90 +++ b/physics/module_MYNNSFC_wrapper.F90 @@ -3,9 +3,15 @@ MODULE mynnsfc_wrapper + USE module_sf_mynn + contains subroutine mynnsfc_wrapper_init () + + ! initialize tables for psih and psim (stable and unstable) + CALL PSI_INIT + end subroutine mynnsfc_wrapper_init subroutine mynnsfc_wrapper_finalize () @@ -19,46 +25,56 @@ end subroutine mynnsfc_wrapper_finalize !! #endif !###=================================================================== -SUBROUTINE mynnsfc_wrapper_run( & - & ix,im,levs, & - & iter,flag_init,flag_restart, & - & delt,dx, & - & u, v, t3d, qvsh, qc, prsl, phii,& - & exner, tsq, qsq, cov, sh3d, & - & el_pbl, qc_bl, cldfra_bl, & - & ps, PBLH, slmsk, TSK, & - & QSFC, snowd, & - & zorl,UST,USTM, ZOL,MOL,RMOL, & - & fm, fh, fm10, fh2, WSPD, br, ch,& - & HFLX, QFX, LH, FLHC, FLQC, & - & U10, V10, TH2, T2, Q2, & - & wstar, CHS2, CQS2, & - & cda, cka, stress, & +SUBROUTINE mynnsfc_wrapper_run( & + & ix,im,levs, & + & itimestep,iter, & + & flag_init,flag_restart,lsm, & + & delt,dx, & + & u, v, t3d, qvsh, qc, prsl, phii, & + & exner, ps, PBLH, slmsk, & + & wet, dry, icy, & !intent(in) + & tskin_ocn, tskin_lnd, tskin_ice, & !intent(in) + & tsurf_ocn, tsurf_lnd, tsurf_ice, & !intent(in) + & qsfc_ocn, qsfc_lnd, qsfc_ice, & !intent(in) + & snowh_ocn, snowh_lnd, snowh_ice, & !intent(in) + & znt_ocn, znt_lnd, znt_ice, & !intent(inout) + & ust_ocn, ust_lnd, ust_ice, & !intent(inout) + & cm_ocn, cm_lnd, cm_ice, & !intent(inout) + & ch_ocn, ch_lnd, ch_ice, & !intent(inout) + & rb_ocn, rb_lnd, rb_ice, & !intent(inout) + & stress_ocn,stress_lnd,stress_ice, & !intent(inout) + & fm_ocn, fm_lnd, fm_ice, & !intent(inout) + & fh_ocn, fh_lnd, fh_ice, & !intent(inout) + & fm10_ocn, fm10_lnd, fm10_ice, & !intent(inout) + & fh2_ocn, fh2_lnd, fh2_ice, & !intent(inout) + & QSFC, qsfc_ruc, USTM, ZOL, MOL, & + & RMOL, WSPD, ch, HFLX, QFLX, LH, & + & FLHC, FLQC, & + & U10, V10, TH2, T2, Q2, & + & wstar, CHS2, CQS2, & ! & CP, G, ROVCP, R, XLV, & ! & SVP1, SVP2, SVP3, SVPT0, & ! & EP1,EP2,KARMAN, & - & icloud_bl, bl_mynn_cloudpdf, & & lprnt, errmsg, errflg ) ! should be moved to inside the mynn: use machine , only : kind_phys -! use funcphys, only : fpvs - - use physcons, only : cp => con_cp, & - & g => con_g, & - & r_d => con_rd, & - & r_v => con_rv, & - & cpv => con_cvap, & - & cliq => con_cliq, & - & Cice => con_csol, & - & rcp => con_rocp, & - & XLV => con_hvap, & - & XLF => con_hfus, & - & EP_1 => con_fvirt, & - & EP_2 => con_eps - - USE module_sf_mynn, only : SFCLAY_mynn + +! use physcons, only : cp => con_cp, & +! & g => con_g, & +! & r_d => con_rd, & +! & r_v => con_rv, & +! & cpv => con_cvap, & +! & cliq => con_cliq, & +! & Cice => con_csol, & +! & rcp => con_rocp, & +! & XLV => con_hvap, & +! & XLF => con_hfus, & +! & EP_1 => con_fvirt, & +! & EP_2 => con_eps + +! USE module_sf_mynn, only : SFCLAY_mynn !------------------------------------------------------------------- implicit none @@ -73,50 +89,13 @@ SUBROUTINE mynnsfc_wrapper_run( & real(kind=kind_phys), parameter :: SVP3 = 29.65 real(kind=kind_phys), parameter :: SVPT0 = 273.15 -!------------------------------------------------------------------- -!For WRF: -!------------------------------------------------------------------- -! USE module_model_constants, only: & -! &karman, g, p1000mb, & -! &cp, r_d, r_v, rcp, xlv, xlf, xls, & -! &svp1, svp2, svp3, svpt0, ep_1, ep_2, rvovrd, & -! &cpv, cliq, cice - -!------------------------------------------------------------------- -!For reference -! REAL , PARAMETER :: karman = 0.4 -! REAL , PARAMETER :: g = 9.81 -! REAL , PARAMETER :: r_d = 287. -! REAL , PARAMETER :: cp = 7.*r_d/2. -! REAL , PARAMETER :: r_v = 461.6 -! REAL , PARAMETER :: cpv = 4.*r_v -! REAL , PARAMETER :: cliq = 4190. -! REAL , PARAMETER :: Cice = 2106. -! REAL , PARAMETER :: rcp = r_d/cp -! REAL , PARAMETER :: XLS = 2.85E6 -! REAL , PARAMETER :: XLV = 2.5E6 -! REAL , PARAMETER :: XLF = 3.50E5 -! REAL , PARAMETER :: p1000mb = 100000. -! REAL , PARAMETER :: rvovrd = r_v/r_d -! REAL , PARAMETER :: SVP1 = 0.6112 -! REAL , PARAMETER :: SVP2 = 17.67 -! REAL , PARAMETER :: SVP3 = 29.65 -! REAL , PARAMETER :: SVPT0 = 273.15 -! REAL , PARAMETER :: EP_1 = R_v/R_d-1. -! REAL , PARAMETER :: EP_2 = R_d/R_v - REAL, PARAMETER :: xlvcp=xlv/cp, xlscp=(xlv+xlf)/cp, ev=xlv, rd=r_d, & - &rk=cp/rd, svp11=svp1*1.e3, p608=ep_1, ep_3=1.-ep_2, g_inv=1/g + &rk=cp/rd, svp11=svp1*1.e3, p608=ep_1, ep_3=1.-ep_2, g_inv=1./g character(len=*), intent(out) :: errmsg integer, intent(out) :: errflg -! NAMELIST OPTIONS (INPUT): - INTEGER, INTENT(IN) :: & - & bl_mynn_cloudpdf, & - & icloud_bl - !MISC CONFIGURATION OPTIONS INTEGER, PARAMETER :: & & spp_pbl = 0, & @@ -127,168 +106,215 @@ SUBROUTINE mynnsfc_wrapper_run( & !MYNN-1D REAL :: delt INTEGER :: im, ix, levs - INTEGER :: iter, k, i, itimestep + INTEGER :: iter, k, i, itimestep, lsm LOGICAL :: flag_init,flag_restart,lprnt INTEGER :: IDS,IDE,JDS,JDE,KDS,KDE, & & IMS,IME,JMS,JME,KMS,KME, & & ITS,ITE,JTS,JTE,KTS,KTE -!MYNN-3D - real(kind=kind_phys), dimension(im,levs+1) :: phii - real(kind=kind_phys), dimension(im,levs) :: & - & exner, PRSL, & - & u, v, t3d, qvsh, qc, & - & Sh3D, EL_PBL, EXCH_H, & - & qc_bl, cldfra_bl, & - & Tsq, Qsq, Cov - !LOCAL + real(kind=kind_phys), dimension(im,levs+1), & + & intent(in) :: phii + real(kind=kind_phys), dimension(im,levs), & + & intent(in) :: exner, PRSL, & + & u, v, t3d, qvsh, qc + real(kind=kind_phys), dimension(im,levs) :: & - & dz, rho, th, qv, & - & pattern_spp_pbl + & pattern_spp_pbl, dz, th, qv + + logical, dimension(im), intent(in) :: wet, dry, icy + + real(kind=kind_phys), dimension(im), intent(in) :: & + & tskin_ocn, tskin_lnd, tskin_ice, & + & tsurf_ocn, tsurf_lnd, tsurf_ice, & + & snowh_ocn, snowh_lnd, snowh_ice + + real(kind=kind_phys), dimension(im), intent(inout) :: & + & znt_ocn, znt_lnd, znt_ice, & + & ust_ocn, ust_lnd, ust_ice, & + & cm_ocn, cm_lnd, cm_ice, & + & ch_ocn, ch_lnd, ch_ice, & + & rb_ocn, rb_lnd, rb_ice, & + & stress_ocn,stress_lnd,stress_ice, & + & fm_ocn, fm_lnd, fm_ice, & + & fh_ocn, fh_lnd, fh_ice, & + & fm10_ocn, fm10_lnd, fm10_ice, & + & fh2_ocn, fh2_lnd, fh2_ice, & + & qsfc_ocn, qsfc_lnd, qsfc_ice !MYNN-2D - real(kind=kind_phys), dimension(im) :: & - & dx, pblh, slmsk, tsk, qsfc, ps, & - & zorl, ust, ustm, hflx, qfx, br, wspd, snowd, & + real(kind=kind_phys), dimension(im), intent(in) :: & + & dx, pblh, slmsk, ps + + real(kind=kind_phys), dimension(im), intent(inout) :: & + & ustm, hflx, qflx, wspd, qsfc, qsfc_ruc, & & FLHC, FLQC, U10, V10, TH2, T2, Q2, & & CHS2, CQS2, rmol, zol, mol, ch, & - & fm, fh, fm10, fh2, & - & lh, cda, cka, stress, wstar + & lh, wstar !LOCAL real, dimension(im) :: & - & qcg, hfx, znt, ts, snowh, psim, psih, & - & chs, ck, cd, mavail, regime, xland, GZ1OZ0 + & hfx, znt, ts, psim, psih, & + & chs, ck, cd, mavail, xland, GZ1OZ0, & + & cpm, qgh, qfx ! Initialize CCPP error handling variables errmsg = '' errflg = 0 - if (lprnt) then - write(0,*)"==============================================" - write(0,*)"in mynn surface layer wrapper..." - write(0,*)"flag_init=",flag_init - write(0,*)"flag_restart=",flag_restart - write(0,*)"iter=",iter - endif - - ! If initialization is needed and mynnsfc_wrapper is called - ! in a subcycling loop, then test for (flag_init==.T. .and. iter==1); - ! initialization in sfclay_mynn is triggered by itimestep == 1 - ! DH* TODO: Use flag_restart to distinguish which fields need - ! to be initialized and which are read from restart files - if (flag_init.and.iter==1) then - itimestep = 1 - else - itimestep = 2 - endif - - !prep MYNN-only variables - do k=1,levs - do i=1,im - dz(i,k)=(phii(i,k+1) - phii(i,k))*g_inv - th(i,k)=t3d(i,k)/exner(i,k) - !qc(i,k)=MAX(qgrs(i,k,ntcw),0.0) - qv(i,k)=qvsh(i,k)/(1.0 - qvsh(i,k)) - rho(i,k)=prsl(i,k)/(r_d*t3d(i,k)) !gt0(i,k)) - pattern_spp_pbl(i,k)=0.0 - enddo - enddo - do i=1,im - if (slmsk(i)==1. .or. slmsk(i)==2.)then !sea/land/ice mask (=0/1/2) in FV3 - xland(i)=1.0 !but land/water = (1/2) in SFCLAY_mynn - else - xland(i)=2.0 - endif -! ust(i) = sqrt(stress(i)) - !ch(i)=0.0 - HFX(i)=hflx(i)*rho(i,1)*cp - !QFX(i)=evap(i) - !wstar(i)=0.0 - qcg(i)=0.0 - snowh(i)=snowd(i)*800. !mm -> m - znt(i)=zorl(i)*0.01 !cm -> m? - ts(i)=tsk(i)/exner(i,1) !theta -! qsfc(i)=qss(i) -! ps(i)=pgr(i) -! wspd(i)=wind(i) - mavail(i)=1.0 !???? - enddo - - if (lprnt) then - write(0,*)"CALLING SFCLAY_mynn; input:" - print*,"T:",t3d(1,1),t3d(1,2),t3d(1,3) - print*,"TH:",th(1,1),th(1,2),th(1,3) - print*,"rho:",rho(1,1),rho(1,2),rho(1,3) - print*,"u:",u(1,1:3) - !print*,"qv:",qv(1,1:3,1) - print*,"p:",prsl(1,1)," snowh=",snowh(1) - print*,"dz:",dz(1,1)," qsfc=",qsfc(1) - print*,"rmol:",rmol(1)," ust:",ust(1) - print*,"Tsk:",tsk(1)," Thetasurf:",ts(1) - print*,"HFX:",hfx(1)," qfx",qfx(1) - print*,"qsfc:",qsfc(1)," ps:",ps(1) - print*,"wspd:",wspd(1),"br=",br(1) - print*,"znt:",znt(1)," delt=",delt - print*,"im=",im," levs=",levs - print*,"flag_init=",flag_init !," ntcw=",ntcw!," ntk=",ntk - print*,"flag_restart=",flag_restart !," ntcw=",ntcw!," ntk=",ntk - print*,"iter=",iter - !print*,"ncld=",ncld," ntrac(gq0)=",ntrac - print*,"zlvl(1)=",dz(1,1)*0.5 - print*,"PBLH=",pblh(1)," xland=",xland(1) - endif - - - CALL SFCLAY_mynn( & - u3d=u,v3d=v,t3d=t3d,qv3d=qv,p3d=prsl,dz8w=dz, & - CP=cp,G=g,ROVCP=rcp,R=r_d,XLV=xlv, & - PSFCPA=ps,CHS=chs,CHS2=chs2,CQS2=cqs2, & - ZNT=znt,UST=ust,PBLH=pblh,MAVAIL=mavail, & - ZOL=zol,MOL=mol,REGIME=regime,psim=psim,psih=psih, & - psix=fm,psit=fh,psix10=fm10,psit2=fh2, & -! fm=psix,fh=psit,fm10=psix10,fh2=psit2, & - XLAND=xland,HFX=hfx,QFX=qfx,LH=lh,TSK=tsk, & - FLHC=flhc,FLQC=flqc,QSFC=qsfc,RMOL=rmol, & - U10=u10,V10=v10,TH2=th2,T2=t2,Q2=q2,SNOWH=snowh, & - GZ1OZ0=GZ1OZ0,WSPD=wspd,BR=br,ISFFLX=isfflx,DX=dx, & - SVP1=svp1,SVP2=svp2,SVP3=svp3,SVPT0=svpt0, & - EP1=ep_1,EP2=ep_2,KARMAN=karman, & - itimestep=itimestep,ch=ch, & - th3d=th,pi3d=exner,qc3d=qc,rho3d=rho, & - tsq=tsq,qsq=qsq,cov=cov,sh3d=sh3d,el_pbl=el_pbl, & - qcg=qcg,wstar=wstar, & - icloud_bl=icloud_bl,qc_bl=qc_bl,cldfra_bl=cldfra_bl, & - spp_pbl=spp_pbl,pattern_spp_pbl=pattern_spp_pbl, & - ids=1,ide=im, jds=1,jde=1, kds=1,kde=levs, & - ims=1,ime=im, jms=1,jme=1, kms=1,kme=levs, & - its=1,ite=im, jts=1,jte=1, kts=1,kte=levs, & - ustm=ustm, ck=ck, cka=cka, cd=cd, cda=cda, & - isftcflx=isftcflx, iz0tlnd=iz0tlnd, & - bl_mynn_cloudpdf=bl_mynn_cloudpdf ) - - - ! POST MYNN SURFACE LAYER (INTERSTITIAL) WORK: - do i = 1, im - hflx(i)=hfx(i)/(rho(i,1)*cp) - !QFX(i)=evap(i) - zorl(i)=znt(i)*100. !m -> cm - stress(i) = ust(i)**2 +! if (lprnt) then +! write(0,*)"==============================================" +! write(0,*)"in mynn surface layer wrapper..." +! write(0,*)"flag_init=",flag_init +! write(0,*)"flag_restart=",flag_restart +! write(0,*)"iter=",iter +! endif + + ! prep MYNN-only variables + do k=1,2 !levs + do i=1,im + dz(i,k)=(phii(i,k+1) - phii(i,k))*g_inv + th(i,k)=t3d(i,k)/exner(i,k) + !qc(i,k)=MAX(qgrs(i,k,ntcw),0.0) + qv(i,k)=qvsh(i,k)/(1.0 - qvsh(i,k)) + pattern_spp_pbl(i,k)=0.0 enddo - - - if (lprnt) then - print* - print*,"finished with mynn_surface layer; output:" - print*,"xland=",xland(1)," cda=",cda(1) - print*,"rmol:",rmol(1)," ust:",ust(1) - print*,"Tsk:",tsk(1)," Thetasurf:",ts(1) - print*,"HFX:",hfx(1)," qfx",qfx(1) - print*,"qsfc:",qsfc(1)," ps:",ps(1) - print*,"wspd:",wspd(1)," br=",br(1) - print*,"znt:",znt(1),"pblh:",pblh(1) - print*,"FLHC=",FLHC(1)," CHS=",CHS(1) - print* - endif + enddo + do i=1,im + if (slmsk(i)==1. .or. slmsk(i)==2.)then !sea/land/ice mask (=0/1/2) in FV3 + xland(i)=1.0 !but land/water = (1/2) in SFCLAY_mynn + else + xland(i)=2.0 + endif + qgh(i)=0.0 + !snowh(i)=snowd(i)*800. !mm -> m + !znt_lnd(i)=znt_lnd(i)*0.01 !cm -> m + !znt_ocn(i)=znt_ocn(i)*0.01 !cm -> m + !znt_ice(i)=znt_ice(i)*0.01 !cm -> m + ! DH* do the following line only if wet(i)? + ts(i)=tskin_ocn(i)/exner(i,1) !theta + ! *DH + mavail(i)=1.0 !???? + cpm(i)=cp + enddo + + ! cm -> m + where (dry) znt_lnd=znt_lnd*0.01 + where (wet) znt_ocn=znt_ocn*0.01 + where (icy) znt_ice=znt_ice*0.01 + +! if (lprnt) then +! write(0,*)"CALLING SFCLAY_mynn; input:" +! write(0,*)"T:",t3d(1,1),t3d(1,2),t3d(1,3) +! write(0,*)"TH:",th(1,1),th(1,2),th(1,3) +! write(0,*)"u:",u(1,1:3) +! write(0,*)"v:",v(1,1:3) +! !write(0,*)"qv:",qv(1,1:3,1) +! write(0,*)"p:",prsl(1,1) +! write(0,*)"dz:",dz(1,1)," qsfc=",qsfc(1)," rmol:",rmol(1) +! write(0,*)" land water ice" +! write(0,*)dry(1),wet(1),icy(1) +! write(0,*)"ust:",ust_lnd(1),ust_ocn(1),ust_ice(1) +! write(0,*)"Tsk:",tskin_lnd(1),tskin_ocn(1),tskin_ice(1) +! write(0,*)"Tsurf:",tsurf_lnd(1),tsurf_ocn(1),tsurf_ice(1) +! write(0,*)"Qsfc:",qsfc_lnd(1),qsfc_ocn(1),qsfc_ice(1) +! write(0,*)"sno:",snowh_lnd(1),snowh_ocn(1),snowh_ice(1) +! write(0,*)"znt:",znt_lnd(1),znt_ocn(1),znt_ice(1) +! !write(0,*)"HFX:",hfx(1)," qfx",qfx(1) +! write(0,*)"qsfc:",qsfc(1)," ps:",ps(1) +! write(0,*)"wspd:",wspd(1),"rb=",rb_ocn(1) +! write(0,*)"delt=",delt," im=",im," levs=",levs +! write(0,*)"flag_init=",flag_init +! write(0,*)"flag_restart=",flag_restart +! write(0,*)"iter=",iter +! write(0,*)"zlvl(1)=",dz(1,1)*0.5 +! write(0,*)"PBLH=",pblh(1)," xland=",xland(1) +! endif + + + CALL SFCLAY_mynn( & + u3d=u,v3d=v,t3d=t3d,qv3d=qv,p3d=prsl,dz8w=dz, & + th3d=th,pi3d=exner,qc3d=qc, & + PSFCPA=ps,PBLH=pblh,MAVAIL=mavail,XLAND=xland,DX=dx, & + CP=cp,G=g,ROVCP=rcp,R=r_d,XLV=xlv, & + SVP1=svp1,SVP2=svp2,SVP3=svp3,SVPT0=svpt0, & + EP1=ep_1,EP2=ep_2,KARMAN=karman, & + ISFFLX=isfflx,isftcflx=isftcflx,LSM=lsm, & + iz0tlnd=iz0tlnd,itimestep=itimestep,iter=iter, & + wet=wet, dry=dry, icy=icy, & !intent(in) + tskin_ocn=tskin_ocn, tskin_lnd=tskin_lnd, tskin_ice=tskin_ice, & !intent(in) + tsurf_ocn=tsurf_ocn, tsurf_lnd=tsurf_lnd, tsurf_ice=tsurf_ice, & !intent(in) + qsfc_ocn=qsfc_ocn, qsfc_lnd=qsfc_lnd, qsfc_ice=qsfc_ice, & !intent(in) + snowh_ocn=snowh_ocn, snowh_lnd=snowh_lnd, snowh_ice=snowh_ice, & !intent(in) + znt_ocn=znt_ocn, znt_lnd=znt_lnd, znt_ice=znt_ice, & !intent(inout) + ust_ocn=ust_ocn, ust_lnd=ust_lnd, ust_ice=ust_ice, & !intent(inout) + cm_ocn=cm_ocn, cm_lnd=cm_lnd, cm_ice=cm_ice, & !intent(inout) + ch_ocn=ch_ocn, ch_lnd=ch_lnd, ch_ice=ch_ice, & !intent(inout) + rb_ocn=rb_ocn, rb_lnd=rb_lnd, rb_ice=rb_ice, & !intent(inout) + stress_ocn=stress_ocn,stress_lnd=stress_lnd,stress_ice=stress_ice, & !intent(inout) + fm_ocn=fm_ocn, fm_lnd=fm_lnd, fm_ice=fm_ice, & !intent(inout) + fh_ocn=fh_ocn, fh_lnd=fh_lnd, fh_ice=fh_ice, & !intent(inout) + fm10_ocn=fm10_ocn, fm10_lnd=fm10_lnd, fm10_ice=fm10_ice, & !intent(inout) + fh2_ocn=fh2_ocn, fh2_lnd=fh2_lnd, fh2_ice=fh2_ice, & !intent(inout) + ch=ch,CHS=chs,CHS2=chs2,CQS2=cqs2,CPM=cpm, & + ZNT=znt,USTM=ustm,ZOL=zol,MOL=mol,RMOL=rmol, & + psim=psim,psih=psih, & + HFLX=hflx,HFX=hfx,QFLX=qflx,QFX=qfx,LH=lh,FLHC=flhc,FLQC=flqc, & + QGH=qgh,QSFC=qsfc,QSFC_RUC=qsfc_ruc, & + U10=u10,V10=v10,TH2=th2,T2=t2,Q2=q2, & + GZ1OZ0=GZ1OZ0,WSPD=wspd,wstar=wstar, & + spp_pbl=spp_pbl,pattern_spp_pbl=pattern_spp_pbl, & + ids=1,ide=im, jds=1,jde=1, kds=1,kde=levs, & + ims=1,ime=im, jms=1,jme=1, kms=1,kme=levs, & + its=1,ite=im, jts=1,jte=1, kts=1,kte=levs ) + + + !! POST MYNN SURFACE LAYER (INTERSTITIAL) WORK: + !do i = 1, im + ! !* Taken from sfc_nst.f + ! !* ch = surface exchange coeff heat & moisture(m/s) im + ! !* rch(i) = rho_a(i) * cp * ch(i) * wind(i) + ! !* hflx(i) = rch(i) * (tsurf(i) - theta1(i)) !K m s-1 + ! !* hflx(i)=hfx(i)/(rho(i,1)*cp) - now calculated inside module_sf_mynn.F90 + ! !* Taken from sfc_nst.f + ! !* evap(i) = elocp * rch(i) * (qss(i) - q0(i)) !kg kg-1 m s-1 + ! !NOTE: evap & qflx will be solved for later + ! !qflx(i)=QFX(i)/ + ! !evap(i)=QFX(i) !or /rho ?? + ! ! DH* note - this could be automated (CCPP knows how to convert m to cm) + ! znt_lnd(i)=znt_lnd(i)*100. !m -> cm + ! znt_ocn(i)=znt_ocn(i)*100. + ! znt_ice(i)=znt_ice(i)*100. + !enddo + + ! m -> cm + where (dry) znt_lnd=znt_lnd*100. + where (wet) znt_ocn=znt_ocn*100. + where (icy) znt_ice=znt_ice*100. + +! if (lprnt) then +! write(0,*) +! write(0,*)"finished with mynn_surface layer; output:" +! write(0,*)" land water ice" +! write(0,*)dry(1),wet(1),icy(1) +! write(0,*)"ust:",ust_lnd(1),ust_ocn(1),ust_ice(1) +! write(0,*)"Tsk:",tskin_lnd(1),tskin_ocn(1),tskin_ice(1) +! write(0,*)"Tsurf:",tsurf_lnd(1),tsurf_ocn(1),tsurf_ice(1) +! write(0,*)"Qsfc:",qsfc_lnd(1),qsfc_ocn(1),qsfc_ice(1) +! write(0,*)"sno:",snowh_lnd(1),snowh_ocn(1),snowh_ice(1) +! write(0,*)"znt (cm):",znt_lnd(1),znt_ocn(1),znt_ice(1) +! write(0,*)"cm:",cm_lnd(1),cm_ocn(1),cm_ice(1) +! write(0,*)"ch:",ch_lnd(1),ch_ocn(1),ch_ice(1) +! write(0,*)"fm:",fm_lnd(1),fm_ocn(1),fm_ice(1) +! write(0,*)"fh:",fh_lnd(1),fh_ocn(1),fh_ice(1) +! write(0,*)"rb:",rb_lnd(1),rb_ocn(1),rb_ice(1) +! write(0,*)"xland=",xland(1)," wstar:",wstar(1) +! write(0,*)"HFX:",hfx(1)," qfx:",qfx(1) +! write(0,*)"HFLX:",hflx(1)," evap:",evap(1) +! write(0,*)"qsfc:",qsfc(1)," ps:",ps(1)," wspd:",wspd(1) +! write(0,*)"ZOL:",ZOL(1)," rmol=",rmol(1) +! write(0,*)"psim:",psim(1)," psih=",psih(1)," pblh:",pblh(1) +! write(0,*)"FLHC=",FLHC(1)," CHS=",CHS(1) +! write(0,*) +! endif END SUBROUTINE mynnsfc_wrapper_run diff --git a/physics/module_MYNNSFC_wrapper.meta b/physics/module_MYNNSFC_wrapper.meta index 2f877075c..b12837233 100644 --- a/physics/module_MYNNSFC_wrapper.meta +++ b/physics/module_MYNNSFC_wrapper.meta @@ -25,6 +25,14 @@ type = integer intent = in optional = F +[itimestep] + standard_name = index_of_time_step + long_name = current number of time steps + units = index + dimensions = () + type = integer + intent = in + optional = F [iter] standard_name = ccpp_loop_counter long_name = loop counter for subcycling loops in CCPP @@ -49,6 +57,14 @@ type = logical intent = in optional = F +[lsm] + standard_name = flag_for_land_surface_scheme + long_name = flag for land surface model + units = flag + dimensions = () + type = integer + intent = in + optional = F [delt] standard_name = time_step_for_physics long_name = time step for physics @@ -139,116 +155,149 @@ kind = kind_phys intent = in optional = F -[tsq] - standard_name = t_prime_squared - long_name = temperature fluctuation squared - units = K2 - dimensions = (horizontal_dimension,vertical_dimension) +[ps] + standard_name = surface_air_pressure + long_name = surface pressure + units = Pa + dimensions = (horizontal_dimension) type = real kind = kind_phys intent = in optional = F -[qsq] - standard_name = q_prime_squared - long_name = water vapor fluctuation squared - units = kg2 kg-2 - dimensions = (horizontal_dimension,vertical_dimension) +[PBLH] + standard_name = atmosphere_boundary_layer_thickness + long_name = PBL thickness + units = m + dimensions = (horizontal_dimension) type = real kind = kind_phys intent = in optional = F -[cov] - standard_name = t_prime_q_prime - long_name = covariance of temperature and moisture - units = K kg kg-1 - dimensions = (horizontal_dimension,vertical_dimension) +[slmsk] + standard_name = sea_land_ice_mask_real + long_name = landmask: sea/land/ice=0/1/2 + units = flag + dimensions = (horizontal_dimension) type = real kind = kind_phys intent = in optional = F -[el_pbl] - standard_name = mixing_length - long_name = mixing length in meters - units = m - dimensions = (horizontal_dimension,vertical_dimension) +[wet] + standard_name = flag_nonzero_wet_surface_fraction + long_name = flag indicating presence of some ocean or lake surface area fraction + units = flag + dimensions = (horizontal_dimension) + type = logical + intent = in + optional = F +[dry] + standard_name = flag_nonzero_land_surface_fraction + long_name = flag indicating presence of some land surface area fraction + units = flag + dimensions = (horizontal_dimension) + type = logical + intent = in + optional = F +[icy] + standard_name = flag_nonzero_sea_ice_surface_fraction + long_name = flag indicating presence of some sea ice surface area fraction + units = flag + dimensions = (horizontal_dimension) + type = logical + intent = in + optional = F +[tskin_ocn] + standard_name = surface_skin_temperature_over_ocean_interstitial + long_name = surface skin temperature over ocean (temporary use as interstitial) + units = K + dimensions = (horizontal_dimension) type = real kind = kind_phys intent = in optional = F -[Sh3D] - standard_name = stability_function_for_heat - long_name = stability function for heat - units = none - dimensions = (horizontal_dimension,vertical_dimension) +[tskin_lnd] + standard_name = surface_skin_temperature_over_land_interstitial + long_name = surface skin temperature over land (temporary use as interstitial) + units = K + dimensions = (horizontal_dimension) type = real kind = kind_phys intent = in optional = F -[QC_BL] - standard_name = subgrid_cloud_mixing_ratio_pbl - long_name = subgrid cloud cloud mixing ratio from PBL scheme - units = kg kg-1 - dimensions = (horizontal_dimension,vertical_dimension) +[tskin_ice] + standard_name = surface_skin_temperature_over_ice_interstitial + long_name = surface skin temperature over ice (temporary use as interstitial) + units = K + dimensions = (horizontal_dimension) type = real kind = kind_phys intent = in optional = F -[CLDFRA_BL] - standard_name = subgrid_cloud_fraction_pbl - long_name = subgrid cloud fraction from PBL scheme - units = frac - dimensions = (horizontal_dimension,vertical_dimension) +[tsurf_ocn] + standard_name = surface_skin_temperature_after_iteration_over_ocean + long_name = surface skin temperature after iteration over ocean + units = K + dimensions = (horizontal_dimension) type = real kind = kind_phys intent = in optional = F -[ps] - standard_name = surface_air_pressure - long_name = surface pressure - units = Pa +[tsurf_lnd] + standard_name = surface_skin_temperature_after_iteration_over_land + long_name = surface skin temperature after iteration over land + units = K dimensions = (horizontal_dimension) type = real kind = kind_phys intent = in optional = F -[PBLH] - standard_name = atmosphere_boundary_layer_thickness - long_name = PBL thickness - units = m +[tsurf_ice] + standard_name = surface_skin_temperature_after_iteration_over_ice + long_name = surface skin temperature after iteration over ice + units = K dimensions = (horizontal_dimension) type = real kind = kind_phys intent = in optional = F -[slmsk] - standard_name = sea_land_ice_mask_real - long_name = landmask: sea/land/ice=0/1/2 - units = flag +[qsfc_ocn] + standard_name = surface_specific_humidity_over_ocean + long_name = surface air saturation specific humidity over ocean + units = kg kg-1 dimensions = (horizontal_dimension) type = real kind = kind_phys - intent = in + intent = inout optional = F -[tsk] - standard_name = surface_skin_temperature - long_name = surface temperature - units = K +[qsfc_lnd] + standard_name = surface_specific_humidity_over_land + long_name = surface air saturation specific humidity over land + units = kg kg-1 dimensions = (horizontal_dimension) type = real kind = kind_phys - intent = in + intent = inout optional = F -[qsfc] - standard_name = surface_specific_humidity - long_name = surface air saturation specific humidity +[qsfc_ice] + standard_name = surface_specific_humidity_over_ice + long_name = surface air saturation specific humidity over ice units = kg kg-1 dimensions = (horizontal_dimension) type = real kind = kind_phys + intent = inout + optional = F +[snowh_ocn] + standard_name = surface_snow_thickness_water_equivalent_over_ocean + long_name = water equivalent snow depth over ocean + units = mm + dimensions = (horizontal_dimension) + type = real + kind = kind_phys intent = in optional = F -[snowd] - standard_name = surface_snow_thickness_water_equivalent +[snowh_lnd] + standard_name = surface_snow_thickness_water_equivalent_over_land long_name = water equivalent snow depth over land units = mm dimensions = (horizontal_dimension) @@ -256,114 +305,348 @@ kind = kind_phys intent = in optional = F -[zorl] - standard_name = surface_roughness_length - long_name = surface roughness length in cm +[snowh_ice] + standard_name = surface_snow_thickness_water_equivalent_over_ice + long_name = water equivalent snow depth over ice + units = mm + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = in + optional = F +[znt_ocn] + standard_name = surface_roughness_length_over_ocean_interstitial + long_name = surface roughness length over ocean (temporary use as interstitial) units = cm dimensions = (horizontal_dimension) type = real kind = kind_phys intent = inout optional = F -[ust] - standard_name = surface_friction_velocity - long_name = boundary layer parameter +[znt_lnd] + standard_name = surface_roughness_length_over_land_interstitial + long_name = surface roughness length over land (temporary use as interstitial) + units = cm + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[znt_ice] + standard_name = surface_roughness_length_over_ice_interstitial + long_name = surface roughness length over ice (temporary use as interstitial) + units = cm + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[ust_ocn] + standard_name = surface_friction_velocity_over_ocean + long_name = surface friction velocity over ocean units = m s-1 dimensions = (horizontal_dimension) type = real kind = kind_phys intent = inout optional = F -[ustm] - standard_name = surface_friction_velocity_drag - long_name = friction velocity isolated for momentum only +[ust_lnd] + standard_name = surface_friction_velocity_over_land + long_name = surface friction velocity over land units = m s-1 dimensions = (horizontal_dimension) type = real kind = kind_phys intent = inout optional = F -[zol] - standard_name = surface_stability_parameter - long_name = monin obukhov surface stability parameter +[ust_ice] + standard_name = surface_friction_velocity_over_ice + long_name = surface friction velocity over ice + units = m s-1 + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[cm_ocn] + standard_name = surface_drag_coefficient_for_momentum_in_air_over_ocean + long_name = surface exchange coeff for momentum over ocean units = none dimensions = (horizontal_dimension) type = real kind = kind_phys intent = inout optional = F -[mol] - standard_name = theta_star - long_name = temperature flux divided by ustar (temperature scale) - units = K +[cm_lnd] + standard_name = surface_drag_coefficient_for_momentum_in_air_over_land + long_name = surface exchange coeff for momentum over land + units = none dimensions = (horizontal_dimension) type = real kind = kind_phys intent = inout optional = F -[rmol] - standard_name = reciprocal_of_obukhov_length - long_name = one over obukhov length - units = m-1 +[cm_ice] + standard_name = surface_drag_coefficient_for_momentum_in_air_over_ice + long_name = surface exchange coeff for momentum over ice + units = none dimensions = (horizontal_dimension) type = real kind = kind_phys intent = inout optional = F -[fm] - standard_name = Monin_Obukhov_similarity_function_for_momentum - long_name = Monin-Obukhov similarity parameter for momentum +[ch_ocn] + standard_name = surface_drag_coefficient_for_heat_and_moisture_in_air_over_ocean + long_name = surface exchange coeff heat & moisture over ocean units = none dimensions = (horizontal_dimension) type = real kind = kind_phys intent = inout optional = F -[fh] - standard_name = Monin_Obukhov_similarity_function_for_heat - long_name = Monin-Obukhov similarity parameter for heat +[ch_lnd] + standard_name = surface_drag_coefficient_for_heat_and_moisture_in_air_over_land + long_name = surface exchange coeff heat & moisture over land units = none dimensions = (horizontal_dimension) type = real kind = kind_phys intent = inout optional = F -[fm10] - standard_name = Monin_Obukhov_similarity_function_for_momentum_at_10m - long_name = Monin-Obukhov similarity parameter for momentum +[ch_ice] + standard_name = surface_drag_coefficient_for_heat_and_moisture_in_air_over_ice + long_name = surface exchange coeff heat & moisture over ice units = none dimensions = (horizontal_dimension) type = real kind = kind_phys intent = inout optional = F -[fh2] - standard_name = Monin_Obukhov_similarity_function_for_heat_at_2m - long_name = Monin-Obukhov similarity parameter for heat +[rb_ocn] + standard_name = bulk_richardson_number_at_lowest_model_level_over_ocean + long_name = bulk Richardson number at the surface over ocean units = none dimensions = (horizontal_dimension) type = real kind = kind_phys intent = inout optional = F -[wspd] - standard_name = wind_speed_at_lowest_model_layer - long_name = wind speed at lowest model level +[rb_lnd] + standard_name = bulk_richardson_number_at_lowest_model_level_over_land + long_name = bulk Richardson number at the surface over land + units = none + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[rb_ice] + standard_name = bulk_richardson_number_at_lowest_model_level_over_ice + long_name = bulk Richardson number at the surface over ice + units = none + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[stress_ocn] + standard_name = surface_wind_stress_over_ocean + long_name = surface wind stress over ocean + units = m2 s-2 + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[stress_lnd] + standard_name = surface_wind_stress_over_land + long_name = surface wind stress over land + units = m2 s-2 + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[stress_ice] + standard_name = surface_wind_stress_over_ice + long_name = surface wind stress over ice + units = m2 s-2 + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[fm_ocn] + standard_name = Monin_Obukhov_similarity_function_for_momentum_over_ocean + long_name = Monin-Obukhov similarity function for momentum over ocean + units = none + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[fm_lnd] + standard_name = Monin_Obukhov_similarity_function_for_momentum_over_land + long_name = Monin-Obukhov similarity function for momentum over land + units = none + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[fm_ice] + standard_name = Monin_Obukhov_similarity_function_for_momentum_over_ice + long_name = Monin-Obukhov similarity function for momentum over ice + units = none + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[fh_ocn] + standard_name = Monin_Obukhov_similarity_function_for_heat_over_ocean + long_name = Monin-Obukhov similarity function for heat over ocean + units = none + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[fh_lnd] + standard_name = Monin_Obukhov_similarity_function_for_heat_over_land + long_name = Monin-Obukhov similarity function for heat over land + units = none + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[fh_ice] + standard_name = Monin_Obukhov_similarity_function_for_heat_over_ice + long_name = Monin-Obukhov similarity function for heat over ice + units = none + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[fm10_ocn] + standard_name = Monin_Obukhov_similarity_function_for_momentum_at_10m_over_ocean + long_name = Monin-Obukhov similarity parameter for momentum at 10m over ocean + units = none + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[fm10_lnd] + standard_name = Monin_Obukhov_similarity_function_for_momentum_at_10m_over_land + long_name = Monin-Obukhov similarity parameter for momentum at 10m over land + units = none + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[fm10_ice] + standard_name = Monin_Obukhov_similarity_function_for_momentum_at_10m_over_ice + long_name = Monin-Obukhov similarity parameter for momentum at 10m over ice + units = none + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[fh2_ocn] + standard_name = Monin_Obukhov_similarity_function_for_heat_at_2m_over_ocean + long_name = Monin-Obukhov similarity parameter for heat at 2m over ocean + units = none + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[fh2_lnd] + standard_name = Monin_Obukhov_similarity_function_for_heat_at_2m_over_land + long_name = Monin-Obukhov similarity parameter for heat at 2m over land + units = none + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[fh2_ice] + standard_name = Monin_Obukhov_similarity_function_for_heat_at_2m_over_ice + long_name = Monin-Obukhov similarity parameter for heat at 2m over ice + units = none + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[qsfc] + standard_name = surface_specific_humidity + long_name = surface air saturation specific humidity + units = kg kg-1 + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[qsfc_ruc] + standard_name = water_vapor_mixing_ratio_at_surface + long_name = water vapor mixing ratio at surface + units = kg kg-1 + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[ustm] + standard_name = surface_friction_velocity_drag + long_name = friction velocity isolated for momentum only units = m s-1 dimensions = (horizontal_dimension) type = real kind = kind_phys intent = inout optional = F -[br] - standard_name = bulk_richardson_number_at_lowest_model_level - long_name = bulk Richardson number at the surface +[zol] + standard_name = surface_stability_parameter + long_name = monin obukhov surface stability parameter units = none dimensions = (horizontal_dimension) type = real kind = kind_phys intent = inout optional = F +[mol] + standard_name = theta_star + long_name = temperature flux divided by ustar (temperature scale) + units = K + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[rmol] + standard_name = reciprocal_of_obukhov_length + long_name = one over obukhov length + units = m-1 + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[wspd] + standard_name = wind_speed_at_lowest_model_layer + long_name = wind speed at lowest model level + units = m s-1 + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F [ch] standard_name = surface_drag_wind_speed_for_momentum_in_air long_name = momentum exchange coefficient @@ -382,7 +665,7 @@ kind = kind_phys intent = inout optional = F -[QFX] +[qflx] standard_name = kinematic_surface_upward_latent_heat_flux long_name = kinematic surface upward latent heat flux units = kg kg-1 m s-1 @@ -490,49 +773,6 @@ kind = kind_phys intent = inout optional = F -[cda] - standard_name = surface_drag_coefficient_for_momentum_in_air - long_name = surface exchange coeff for momentum - units = none - dimensions = (horizontal_dimension) - type = real - kind = kind_phys - intent = inout - optional = F -[cka] - standard_name = surface_drag_coefficient_for_heat_and_moisture_in_air - long_name = surface exchange coeff heat & moisture - units = none - dimensions = (horizontal_dimension) - type = real - kind = kind_phys - intent = inout - optional = F -[stress] - standard_name = surface_wind_stress - long_name = surface wind stress - units = m2 s-2 - dimensions = (horizontal_dimension) - type = real - kind = kind_phys - intent = inout - optional = F -[bl_mynn_cloudpdf] - standard_name = cloudpdf - long_name = flag to determine which cloud PDF to use - units = flag - dimensions = () - type = integer - intent = in - optional = F -[icloud_bl] - standard_name = couple_sgs_clouds_to_radiation_flag - long_name = flag for coupling sgs clouds to radiation - units = flag - dimensions = () - type = integer - intent = in - optional = F [lprnt] standard_name = flag_print long_name = control flag for diagnostic print out diff --git a/physics/module_sf_mynn.F90 b/physics/module_sf_mynn.F90 index 70b98363d..73ef5e1fb 100644 --- a/physics/module_sf_mynn.F90 +++ b/physics/module_sf_mynn.F90 @@ -8,59 +8,63 @@ MODULE module_sf_mynn !------------------------------------------------------------------- !Modifications implemented by Joseph Olson NOAA/GSD/AMB - CU/CIRES -!for WRFv3.4, v3.4.1, v3.5.1, v3.6, v3.7.1, and v3.9: +!The following overviews the current state of this scheme:: ! ! 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 isfflx=0 option to turn off scalar fluxes, but keep momentum +! for first iteration of first time step; afterwards, exact calculation +! using basically the same iterative technique in the module_sf_sfclayrev.F, +! which leverages Pedro Jimenez's code, and is adapted for MYNN. +!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 +!3) Kinematic viscosity varies with temperature according to Andreas (1989). +!4) Uses the blended Monin-Obukhov flux-profile relationships COARE (Fairall +! et al 2003) for the unstable regime (a blended mix of Dyer-Hicks 1974 and +! Grachev et al (2000). Uses Cheng and Brutsaert (2005) for stable conditions. +!5) The following overviews the namelist variables that control the +! aerodynamic roughness lengths (over water) and the thermal and moisture +! roughness lengths (defaults are recommended): ! ! LAND only: -!1) iz0tlnd option is now available with the following options: -! (default) =0: Zilitinkevich (1995) +! "iz0tlnd" namelist option is used to select the following options: +! (default) =0: Zilitinkevich (1995); Czil now set to 0.085 ! =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) -! =4: Pan et al. (1994) with RUC mods for z_q, zili for z_t -!2) Relaxed u* minimum from 0.1 to 0.01 ! ! WATER only: -!1) isftcflx option is now available with the following options: +! "isftcflx" namelist option is used to select the following options: ! (default) =0: z0, zt, and zq from the COARE algorithm. Set COARE_OPT (below) to ! 3.0 (Fairall et al. 2003, default) ! 3.5 (Edson et al 2013) ! =1: z0 from Davis et al (2008), zt & zq from COARE 3.0/3.5 ! =2: z0 from Davis et al (2008), zt & zq from Garratt (1992) ! =3: z0 from Taylor and Yelland (2004), zt and zq from COARE 3.0/3.5 -! =4: z0 from Zilitinkevich (2001), zt & zq from COARE 3.0/3.5 ! ! 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. Also added a z0 mod for snow (Andreas et al. 2005, BLM), which +! Andreas (2002) snow/ice parameterization for thermal and +! moisture roughness is used over all gridpoints with snow deeper than +! 0.1 m. This algorithm calculates a z0 for snow (Andreas et al. 2005, BLM), +! which is only used as part of the thermal and moisture roughness +! length calculation, not to directly impact the surface winds. ! ! Misc: -! 2) added a more elaborate diagnostic for u10 & V10 for high vertical resolution -! model configurations. +!1) Added a more elaborate diagnostic for u10 & V10 for high vertical resolution +! model configurations but for most model configurations with depth of +! the lowest half-model level near 10 m, a neutral-log diagnostic is used. ! -! New for v3.9: -! - option for stochastic parameter perturbations (SPP) +!2) Option to activate stochastic parameter perturbations (SPP), which +! perturb z0, zt, and zq, along with many other parameters in the MYNN- +! EDMF scheme. ! !NOTE: This code was primarily tested in combination with the RUC LSM. ! Performance with the Noah (or other) LSM is relatively unknown. !------------------------------------------------------------------- !For WRF ! USE module_model_constants, only: & -! &g, p1000mb, cp, xlv, ep_2, r_d, r_v, rcp, cpv +! & p1000mb, ep_2 ! - USE module_bl_mynn, only: tv0, b1, b2, p608, ev, rd, & !, mym_condensation - &esat_blend, xl_blend, qsat_blend - +!For non-WRF use physcons, only : cp => con_cp, & & g => con_g, & & r_d => con_rd, & @@ -89,52 +93,81 @@ MODULE module_sf_mynn REAL , PARAMETER :: p1000mb = 100000. ! REAL , PARAMETER :: EP_2 = r_d/r_v - - REAL, PARAMETER :: xlvcp=xlv/cp, ep_3=1.-ep_2 REAL, PARAMETER :: wmin=0.1 ! Minimum wind speed REAL, PARAMETER :: VCONVC=1.25 + REAL, PARAMETER :: onethird = 1./3. + REAL, PARAMETER :: sqrt3 = 1.7320508075688773 + REAL, PARAMETER :: atan1 = 0.785398163397 !in radians REAL, PARAMETER :: SNOWZ0=0.011 REAL, PARAMETER :: COARE_OPT=3.0 ! 3.0 or 3.5 !For debugging purposes: - LOGICAL, PARAMETER :: debug_code = .false. + INTEGER, PARAMETER :: debug_code = 0 !0: no extra ouput + !1: some step-by-step output + !2: everything - heavy I/O + LOGICAL, PARAMETER :: compute_diag = .false. + LOGICAL, PARAMETER :: compute_flux = .false. !shouldn't need compute + ! these in FV3. They will be written over anyway. + ! Computing the fluxes here is leftover from the WRF world. + + REAL, DIMENSION(0:1000 ),SAVE :: psim_stab,psim_unstab, & + psih_stab,psih_unstab CONTAINS !------------------------------------------------------------------- !>\ingroup module_sf_mynn_mod -!> 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). +!> Fill the PSIM and PSIH tables. The subroutine "psi_init" was leveraged from +!! module_sf_sfclayrev.F, leveraging the work from Pedro Jimenez. +!! This subroutine returns a blended form from Dyer and Hicks (1974) +!! and Grachev et al (2000) for unstable conditions and the form +!! from Cheng and Brutsaert (2005) for stable conditions. + SUBROUTINE mynn_sf_init_driver(allowed_to_read) LOGICAL, INTENT(in) :: allowed_to_read -! CALL sfclayinit(allowed_to_read) + CALL psi_init END SUBROUTINE mynn_sf_init_driver !------------------------------------------------------------------- !>\ingroup module_sf_mynn_mod !! This subroutine - SUBROUTINE SFCLAY_mynn( & - U3D,V3D,T3D,QV3D,P3D,dz8w, & - CP,G,ROVCP,R,XLV,PSFCPA,CHS,CHS2,CQS2, & - ZNT,UST,PBLH,MAVAIL,ZOL,MOL,REGIME, & - PSIM,PSIH,PSIX,PSIX10,PSIT,PSIT2, & - XLAND,HFX,QFX,LH,TSK,FLHC,FLQC,QSFC,RMOL, & - U10,V10,TH2,T2,Q2,SNOWH, & - GZ1OZ0,WSPD,BR,ISFFLX,DX, & - SVP1,SVP2,SVP3,SVPT0,EP1,EP2, & - KARMAN,itimestep,ch,th3d,pi3d,qc3d,rho3d, & - tsq,qsq,cov,sh3d,el_pbl,qcg,wstar, & - icloud_bl,qc_bl,cldfra_bl, & - spp_pbl,pattern_spp_pbl, & - 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, & - bl_mynn_cloudpdf) + SUBROUTINE SFCLAY_mynn( & + U3D,V3D,T3D,QV3D,P3D,dz8w, & !in + th3d,pi3d,qc3d, & !in + PSFCPA,PBLH,MAVAIL,XLAND,DX, & !in + CP,G,ROVCP,R,XLV, & !in + SVP1,SVP2,SVP3,SVPT0,EP1,EP2,KARMAN, & !in + ISFFLX,isftcflx,lsm,iz0tlnd, & !in + itimestep,iter, & !in + wet, dry, icy, & !intent(in) + tskin_ocn, tskin_lnd, tskin_ice, & !intent(in) + tsurf_ocn, tsurf_lnd, tsurf_ice, & !intent(in) + qsfc_ocn, qsfc_lnd, qsfc_ice, & !intent(in) + snowh_ocn, snowh_lnd, snowh_ice, & !intent(in) + ZNT_ocn, ZNT_lnd, ZNT_ice, & !intent(inout) + UST_ocn, UST_lnd, UST_ice, & !intent(inout) + cm_ocn, cm_lnd, cm_ice, & !intent(inout) + ch_ocn, ch_lnd, ch_ice, & !intent(inout) + rb_ocn, rb_lnd, rb_ice, & !intent(inout) + stress_ocn,stress_lnd,stress_ice, & !intent(inout) + fm_ocn, fm_lnd, fm_ice, & !intent(inout) + fh_ocn, fh_lnd, fh_ice, & !intent(inout) + fm10_ocn, fm10_lnd, fm10_ice, & !intent(inout) + fh2_ocn, fh2_lnd, fh2_ice, & !intent(inout) + CH,CHS,CHS2,CQS2,CPM, & + ZNT,USTM,ZOL,MOL,RMOL, & + PSIM,PSIH, & + HFLX,HFX,QFLX,QFX,LH,FLHC,FLQC, & + QGH,QSFC,QSFC_RUC, & + U10,V10,TH2,T2,Q2, & + GZ1OZ0,WSPD,WSTAR, & + spp_pbl,pattern_spp_pbl, & + ids,ide, jds,jde, kds,kde, & + ims,ime, jms,jme, kms,kme, & + its,ite, jts,jte, kts,kte ) !------------------------------------------------------------------- IMPLICIT NONE !------------------------------------------------------------------- @@ -143,7 +176,6 @@ SUBROUTINE SFCLAY_mynn( & !-- T3D 3D temperature (K) !-- QV3D 3D water vapor mixing ratio (Kg/Kg) !-- P3D 3D pressure (Pa) -!-- RHO3D 3D density (kg/m3) !-- dz8w 3D 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) @@ -166,7 +198,11 @@ SUBROUTINE SFCLAY_mynn( & !-- 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) +! HFX = HFLX * rho * cp +!-- HFLX upward temperature flux at the surface (K m s^-1) !-- QFX upward moisture flux at the surface (kg/m^2/s) +! QFX = QFLX * rho +!-- QFLX upward moisture flux at the surface (kg kg-1 m s-1) !-- LH net upward latent heat flux at surface (W/m^2) !-- TSK surface temperature (K) !-- FLHC exchange coefficient for heat (W/m^2/K) @@ -202,22 +238,10 @@ SUBROUTINE SFCLAY_mynn( & ! (water =1: z0 from Davis et al (2008), zt & zq from COARE3.0/3.5 ! only) =2: z0 from Davis et al (2008), zt & zq from Garratt (1992) ! =3: z0 from Taylor and Yelland (2004), zt and zq from COARE 3.0/3.5 -! =4: z0 from Zilitinkevich (2001), zt & zq from COARE 3.0/3.5 -!-- iz0tlnd =0: Zilitinkevich (1995) with Czil=0.10, +!-- iz0tlnd =0: Zilitinkevich (1995) with Czil=0.085, ! (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) -! =4: Pan et al (1994) for zq; ZIlitintevich for zt -!-- bl_mynn_cloudpdf =0: Mellor & Yamada -! =1: Kuwano et al. -!-- el_pbl = mixing length from PBL scheme (meters) -!-- Sh3d = Stability finction for heat (unitless) -!-- cov = T'q' from PBL scheme -!-- tsq = T'T' from PBL scheme -!-- qsq = q'q' from PBL scheme -!-- icloud_bl = namelist option for subgrid scale cloud/radiation feedback -!-- qc_bl = subgrid scale (bloundary layer) clouds -!-- cldfra_bl = subgridscale cloud fraction ! !-- ids start index for i in domain !-- ide end index for i in domain @@ -243,16 +267,14 @@ SUBROUTINE SFCLAY_mynn( & INTEGER, INTENT(IN) :: ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte - INTEGER, INTENT(IN) :: itimestep + INTEGER, INTENT(IN) :: itimestep,iter REAL, INTENT(IN) :: SVP1,SVP2,SVP3,SVPT0 REAL, INTENT(IN) :: EP1,EP2,KARMAN REAL, INTENT(IN) :: CP,G,ROVCP,R,XLV !,DX !NAMELIST OPTIONS: - INTEGER, INTENT(IN) :: ISFFLX - INTEGER, OPTIONAL, INTENT(IN) :: ISFTCFLX, IZ0TLND,& - bl_mynn_cloudpdf,& - icloud_bl - INTEGER, INTENT(IN),OPTIONAL :: spp_pbl + INTEGER, INTENT(IN) :: ISFFLX, LSM + INTEGER, OPTIONAL, INTENT(IN) :: ISFTCFLX, IZ0TLND + INTEGER, OPTIONAL, INTENT(IN) :: spp_pbl !=================================== ! 3D VARIABLES @@ -264,11 +286,10 @@ SUBROUTINE SFCLAY_mynn( & T3D, & QC3D, & U3D,V3D, & - RHO3D,th3d,pi3d,tsq,qsq,cov,sh3d,el_pbl + th3d,pi3d - REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) :: qc_bl, & - cldfra_bl - REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN),OPTIONAL ::pattern_spp_pbl + REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), OPTIONAL, & + INTENT(IN) :: pattern_spp_pbl !=================================== ! 2D VARIABLES !=================================== @@ -276,85 +297,84 @@ SUBROUTINE SFCLAY_mynn( & INTENT(IN ) :: MAVAIL, & PBLH, & XLAND, & - TSK, & - QCG, & PSFCPA, & - SNOWH, & DX REAL, DIMENSION( ims:ime, jms:jme ) , & INTENT(OUT ) :: U10,V10, & TH2,T2,Q2 - REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ) , & - INTENT(OUT) :: ck,cka,cd,cda,ustm -! + REAL, DIMENSION( ims:ime, jms:jme ) , & - INTENT(INOUT) :: REGIME, & - HFX, & - QFX, & + INTENT(INOUT) :: HFLX,HFX, & + QFLX,QFX, & LH, & MOL,RMOL, & QSFC, & + QGH, & ZNT, & ZOL, & - UST, & + USTM, & + CPM, & CHS2, & CQS2, & CHS, & CH, & FLHC,FLQC, & - GZ1OZ0,WSPD,BR, & + GZ1OZ0,WSPD, & PSIM,PSIH, & - WSTAR, & - PSIX,PSIX10,PSIT,PSIT2 + WSTAR + + LOGICAL, DIMENSION( ims:ime ), INTENT(IN) :: & + & wet, dry, icy + + REAL, DIMENSION( ims:ime ), INTENT(IN) :: & + & tskin_ocn, tskin_lnd, tskin_ice, & + & tsurf_ocn, tsurf_lnd, tsurf_ice, & + & snowh_ocn, snowh_lnd, snowh_ice + + REAL, DIMENSION( ims:ime), INTENT(INOUT) :: & + & ZNT_ocn, ZNT_lnd, ZNT_ice, & + & UST_ocn, UST_lnd, UST_ice, & + & cm_ocn, cm_lnd, cm_ice, & + & ch_ocn, ch_lnd, ch_ice, & + & rb_ocn, rb_lnd, rb_ice, & + & stress_ocn,stress_lnd,stress_ice, & + & fm_ocn, fm_lnd, fm_ice, & + & fh_ocn, fh_lnd, fh_ice, & + & fm10_ocn, fm10_lnd, fm10_ice, & + & fh2_ocn, fh2_lnd, fh2_ice, & + & qsfc_ocn, qsfc_lnd, qsfc_ice, & + & qsfc_ruc !ADDITIONAL OUTPUT !JOE-begin - REAL, DIMENSION( ims:ime, jms:jme ) :: z0zt_ratio, & - BulkRi,qstar,resist,logres -!JOE-end + REAL, DIMENSION( ims:ime, jms:jme ) :: qstar +!JOE-end !=================================== ! 1D LOCAL ARRAYS !=================================== - REAL, DIMENSION( its:ite ) :: U1D, & - V1D, & + REAL, DIMENSION( its:ite ) :: U1D,V1D, & !level1 winds U1D2,V1D2, & !level2 winds QV1D, & P1D, & T1D,QC1D, & - RHO1D, & dz8w1d, & !level 1 height dz2w1d !level 2 height REAL, DIMENSION( its:ite ) :: rstoch1D - ! VARIABLE FOR PASSING TO MYM_CONDENSATION - REAL, DIMENSION(kts:kts+1 ) :: dummy1,dummy2,dummy3,dummy4, & - dummy5,dummy6,dummy7,dummy8, & - dummy9,dummy10,dummy11, & - dummy12,dummy13,dummy14 - - REAL, DIMENSION( its:ite ) :: vt1,vq1 - REAL, DIMENSION(kts:kts+1) :: thl, qw, vt, vq - REAL :: ql - INTEGER :: I,J,K,itf,jtf,ktf !----------------------------------------------------------- -!joe -test printing of constants: -! print*,"cp=", cp -! print*,"g=", g -! print*,"Rd=", r_d -! print*,"Rv=", r_v -! print*,"cpc=", cpv -! print*,"cliq=", cliq -! print*,"cice=", Cice -! print*,"rcp=", rcp -! print*,"xlv=", XLV -! print*,"xlf=", XLF -! print*,"ep1=", EP_1 -! print*,"ep2=", EP_2 + IF (debug_code >= 1) THEN + write(*,*)"======= printing of constants:" + write(*,*)"cp=", cp," g=", g + write(*,*)"Rd=", r_d," Rv=", r_v, " cpc=", cpv + write(*,*)"cliq=", cliq," cice=", Cice," rcp=", rcp + write(*,*)"xlv=", XLV," xlf=", XLF + write(*,*)"ep1=", EP_1, " ep2=", EP_2 + ENDIF itf=ite !MIN0(ite,ide-1) jtf=jte !MIN0(jte,jde-1) @@ -373,7 +393,6 @@ SUBROUTINE SFCLAY_mynn( & QC1D(i)=QC3D(i,kts,j) P1D(i) =P3D(i,kts,j) T1D(i) =T3D(i,kts,j) - RHO1D(i)=RHO3D(i,kts,j) if (spp_pbl==1) then rstoch1D(i)=pattern_spp_pbl(i,kts,j) else @@ -381,104 +400,69 @@ SUBROUTINE SFCLAY_mynn( & endif ENDDO - IF (itimestep==1) THEN + IF (itimestep==1 .AND. iter==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) + !Everything here is used before calculated + UST_OCN(i)=MAX(0.04*SQRT(U1D(i)*U1D(i) + V1D(i)*V1D(i)),0.001) + UST_LND(i)=MAX(0.04*SQRT(U1D(i)*U1D(i) + V1D(i)*V1D(i)),0.001) + UST_ICE(i)=MAX(0.04*SQRT(U1D(i)*U1D(i) + V1D(i)*V1D(i)),0.001) MOL(i,j)=0. ! Tstar QSFC(i,j)=QV3D(i,kts,j)/(1.+QV3D(i,kts,j)) + QSFC_OCN(i)=QSFC(i,j) + QSFC_LND(i)=QSFC(i,j) + QSFC_ICE(i)=QSFC(i,j) qstar(i,j)=0.0 + QFX(i,j)=0. + HFX(i,j)=0. + QFLX(i,j)=0. + HFLX(i,j)=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) - dummy1(k)=dz8w(i,k,j) - dummy2(k)=thl(k) - dummy3(k)=qw(k) - dummy4(k)=p3d(i,k,j) - dummy5(k)=pi3d(i,k,j) - dummy6(k)=tsq(i,k,j) - dummy7(k)=qsq(i,k,j) - dummy8(k)=cov(i,k,j) - dummy9(k)=Sh3d(i,k,j) - dummy10(k)=el_pbl(i,k,j) - dummy14(k)=th3d(i,k,j) - if(icloud_bl > 0) then - dummy11(k)=qc_bl(i,k,j) - dummy12(k)=cldfra_bl(i,k,j) - else - dummy11(k)=0.0 - dummy12(k)=0.0 - endif - dummy13(k)=0.0 !sgm + IF (LSM == 3) THEN + DO i=its,ite + QSFC_LND(i)=QSFC_RUC(i) ENDDO - - ! NOTE: The last grid number is kts+1 instead of kte. - CALL mym_condensation (kts,kts+1, dx(i,j),& - & dummy1,dummy2,dummy3, & - & dummy4,dummy5,dummy6, & - & dummy7,dummy8,dummy9, & - & dummy10,bl_mynn_cloudpdf,& - & dummy11,dummy12, & - & PBLH(i,j),HFX(i,j), & - & vt(kts:kts+1), vq(kts:kts+1), & - & dummy14,dummy13) - -! ! NOTE: The last grid number is kts+1 instead of kte. -! CALL mym_condensation (kts,kts+1, dx, & -! & 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), & -! & Sh3d(i,kts:kts+1,j), & !JOE - cloud PDF testing -! & el_pbl(i,kts:kts+1,j), & !JOE - cloud PDF testing -! & bl_mynn_cloudpdf, & !JOE - cloud PDF testing -! & qc_bl2D(i,kts:kts+1), & !JOE-subgrid BL clouds -! & cldfra_bl2D(i,kts:kts+1),& !JOE-subgrid BL clouds -! & PBLH(i,j),HFX(i,j), & !JOE-subgrid BL clouds -! & vt(kts:kts+1), vq(kts:kts+1), & - ! & th,sgm) - vt1(i) = vt(kts) - vq1(i) = vq(kts) - ENDDO + ENDIF ENDIF - CALL SFCLAY1D_mynn( & - J,U1D,V1D,T1D,QV1D,P1D,dz8w1d,rho1d, & - U1D2,V1D2,dz2w1d, & - CP,G,ROVCP,R,XLV,PSFCPA(ims,j),CHS(ims,j),CHS2(ims,j),& - CQS2(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), & - PSIX(ims,j),PSIX10(ims,j),PSIT(ims,j),PSIT2(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),SNOWH(ims,j), & - QSFC(ims,j),LH(ims,j), & - GZ1OZ0(ims,j),WSPD(ims,j),BR(ims,j),ISFFLX,DX(ims,j),& - SVP1,SVP2,SVP3,SVPT0,EP1,EP2,KARMAN, & - ch(ims,j),vt1,vq1,qc1d,qcg(ims,j), & - itimestep, & -!JOE-begin additional output - z0zt_ratio(ims,j),wstar(ims,j), & - qstar(ims,j),resist(ims,j),logres(ims,j), & -!JOE-end - spp_pbl,rstoch1D, & - 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) & - ) + CALL SFCLAY1D_mynn( & + J,U1D,V1D,T1D,QV1D,P1D,dz8w1d, & + U1D2,V1D2,dz2w1d, & + PSFCPA(ims,j),PBLH(ims,j),MAVAIL(ims,j), & + XLAND(ims,j),DX(ims,j), & + CP,G,ROVCP,R,XLV,SVP1,SVP2,SVP3,SVPT0, & + EP1,EP2,KARMAN, & + ISFFLX,isftcflx,iz0tlnd,itimestep,iter, & + wet, dry, icy, & !intent(in) + tskin_ocn, tskin_lnd, tskin_ice, & !intent(in) + tsurf_ocn, tsurf_lnd, tsurf_ice, & !intent(in) + qsfc_ocn, qsfc_lnd, qsfc_ice, & !intent(in) + snowh_ocn, snowh_lnd, snowh_ice, & !intent(in) + ZNT_ocn, ZNT_lnd, ZNT_ice, & !intent(inout) + UST_ocn, UST_lnd, UST_ice, & !intent(inout) + cm_ocn, cm_lnd, cm_ice, & !intent(inout) + ch_ocn, ch_lnd, ch_ice, & !intent(inout) + rb_ocn, rb_lnd, rb_ice, & !intent(inout) + stress_ocn, stress_lnd, stress_ice, & !intent(inout) + fm_ocn, fm_lnd, fm_ice, & !intent(inout) + fh_ocn, fh_lnd, fh_ice, & !intent(inout) + fm10_ocn, fm10_lnd, fm10_ice, & !intent(inout) + fh2_ocn, fh2_lnd, fh2_ice, & + ch(ims,j),CHS(ims,j),CHS2(ims,j),CQS2(ims,j), & + CPM(ims,j), & + ZNT(ims,j),USTM(ims,j),ZOL(ims,j), & + MOL(ims,j),RMOL(ims,j), & + PSIM(ims,j),PSIH(ims,j), & + HFLX(ims,j),HFX(ims,j),QFLX(ims,j),QFX(ims,j), & + LH(ims,j),FLHC(ims,j),FLQC(ims,j), & + QGH(ims,j),QSFC(ims,j), & + U10(ims,j),V10(ims,j),TH2(ims,j),T2(ims,j),Q2(ims,j),& + GZ1OZ0(ims,j),WSPD(ims,j),wstar(ims,j), & + spp_pbl,rstoch1D, & + ids,ide, jds,jde, kds,kde, & + ims,ime, jms,jme, kms,kme, & + its,ite, jts,jte, kts,kte & + ) ENDDO @@ -486,29 +470,43 @@ END SUBROUTINE SFCLAY_MYNN !------------------------------------------------------------------- !>\ingroup module_sf_mynn_mod -!! This subroutine calculates - SUBROUTINE SFCLAY1D_mynn( & - J,U1D,V1D,T1D,QV1D,P1D,dz8w1d,rho1d, & - U1D2,V1D2,dz2w1d, & - CP,G,ROVCP,R,XLV,PSFCPA,CHS,CHS2,CQS2, & - PBLH,RMOL,ZNT,UST,MAVAIL,ZOL,MOL,REGIME, & - PSIM,PSIH,PSIX,PSIX10,PSIT,PSIT2, & - XLAND,HFX,QFX,TSK, & - U10,V10,TH2,T2,Q2,FLHC,FLQC,SNOWH, & - QSFC,LH,GZ1OZ0,WSPD,BR,ISFFLX,DX, & - SVP1,SVP2,SVP3,SVPT0,EP1,EP2, & - KARMAN,ch,vt1,vq1,qc1d,qcg, & - itimestep, & -!JOE-additional output - zratio,wstar,qstar,resist,logres, & -!JOE-end - spp_pbl,rstoch1D, & - 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 & - ) +!! This subroutine calculates u*, z/L, and the exchange coefficients +!! which are passed to subsequent scheme to calculate the fluxes. +!! This scheme has options to calculate the fluxes and near-surface +!! diagnostics, as was needed in WRF, but these are skipped for FV3. + SUBROUTINE SFCLAY1D_mynn( & + J,U1D,V1D,T1D,QV1D,P1D,dz8w1d,U1D2,V1D2,dz2w1d, & + PSFCPA,PBLH,MAVAIL,XLAND,DX, & + CP,G,ROVCP,R,XLV,SVP1,SVP2,SVP3,SVPT0, & + EP1,EP2,KARMAN, & + ISFFLX,isftcflx,iz0tlnd,itimestep,iter, & + wet, dry, icy, & !intent(in) + tskin_ocn, tskin_lnd, tskin_ice, & !intent(in) + tsurf_ocn, tsurf_lnd, tsurf_ice, & !intent(in) + qsfc_ocn, qsfc_lnd, qsfc_ice, & !intent(in) + snowh_ocn, snowh_lnd, snowh_ice, & !intent(in) + ZNT_ocn, ZNT_lnd, ZNT_ice, & !intent(inout) + UST_ocn, UST_lnd, UST_ice, & !intent(inout) + cm_ocn, cm_lnd, cm_ice, & !intent(inout) + ch_ocn, ch_lnd, ch_ice, & !intent(inout) + rb_ocn, rb_lnd, rb_ice, & !intent(inout) + stress_ocn, stress_lnd, stress_ice, & !intent(inout) + psix_ocn, psix_lnd, psix_ice, & !=fm, intent(inout) + psit_ocn, psit_lnd, psit_ice, & !=fh, intent(inout) + psix10_ocn, psix10_lnd, psix10_ice, & !=fm10, intent(inout) + psit2_ocn, psit2_lnd, psit2_ice, & !=fh2, intent(inout) + ch,CHS,CHS2,CQS2,CPM, & + ZNT,USTM,ZOL,MOL,RMOL, & + PSIM,PSIH, & + HFLX,HFX,QFLX,QFX,LH,FLHC,FLQC, & + QGH,QSFC, & + U10,V10,TH2,T2,Q2, & + GZ1OZ0,WSPD,wstar, & + spp_pbl,rstoch1D, & + ids,ide, jds,jde, kds,kde, & + ims,ime, jms,jme, kms,kme, & + its,ite, jts,jte, kts,kte & + ) !------------------------------------------------------------------- IMPLICIT NONE @@ -518,7 +516,7 @@ SUBROUTINE SFCLAY1D_mynn( & INTEGER, INTENT(IN) :: ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte, & - J, itimestep + J, itimestep, iter REAL, PARAMETER :: XKA=2.4E-5 !molecular diffusivity REAL, PARAMETER :: PRT=1. !prandlt number @@ -538,34 +536,52 @@ SUBROUTINE SFCLAY1D_mynn( & REAL, DIMENSION( ims:ime ), INTENT(IN) :: MAVAIL, & PBLH, & XLAND, & - TSK, & PSFCPA, & - QCG, & - SNOWH, DX + DX REAL, DIMENSION( its:ite ), INTENT(IN) :: U1D,V1D, & U1D2,V1D2, & QV1D,P1D, & - T1D,QC1d, & - dz8w1d,dz2w1d, & - RHO1D, & - vt1,vq1 + T1D, & + dz8w1d, & + dz2w1d - REAL, DIMENSION( ims:ime ), INTENT(INOUT) :: REGIME, & - HFX,QFX,LH, & + REAL, DIMENSION( ims:ime ), INTENT(INOUT) :: HFLX,HFX, & + QFLX,QFX,LH, & MOL,RMOL, & - QSFC, & + QGH,QSFC, & ZNT, & ZOL, & - UST, & + CPM, & CHS2,CQS2, & CHS,CH, & FLHC,FLQC, & GZ1OZ0, & WSPD, & - BR, & - PSIM,PSIH, & - PSIX,PSIX10,PSIT,PSIT2 + PSIM, & + PSIH, & + USTM + + LOGICAL, DIMENSION( ims:ime ), INTENT(IN) :: & + & wet, dry, icy + + REAL, DIMENSION( ims:ime ), INTENT(in) :: & + & tskin_ocn, tskin_lnd, tskin_ice, & + & tsurf_ocn, tsurf_lnd, tsurf_ice, & + & snowh_ocn, snowh_lnd, snowh_ice + + REAL, DIMENSION( ims:ime ), INTENT(inout) :: & + & ZNT_ocn, ZNT_lnd, ZNT_ice, & + & UST_ocn, UST_lnd, UST_ice, & + & cm_ocn, cm_lnd, cm_ice, & + & ch_ocn, ch_lnd, ch_ice, & + & rb_ocn, rb_lnd, rb_ice, & + & stress_ocn,stress_lnd,stress_ice, & + & psix_ocn, psix_lnd, psix_ice, & + & psit_ocn, psit_lnd, psit_ice, & + & psix10_ocn,psix10_lnd,psix10_ice, & + & psit2_ocn, psit2_lnd, psit2_ice, & + & qsfc_ocn, qsfc_lnd, qsfc_ice REAL, DIMENSION( its:ite ), INTENT(IN) :: rstoch1D @@ -573,18 +589,13 @@ SUBROUTINE SFCLAY1D_mynn( & REAL, DIMENSION( ims:ime ), INTENT(OUT) :: U10,V10, & TH2,T2,Q2 - REAL, OPTIONAL, DIMENSION( ims:ime ) , & - INTENT(OUT) :: ck,cka,cd,cda,ustm !-------------------------------------------- !JOE-additinal output - REAL, DIMENSION( ims:ime ) :: zratio,wstar,qstar, & - resist,logres + REAL, DIMENSION( ims:ime ) :: wstar,qstar !JOE-end !---------------------------------------------------------------- ! LOCAL VARS !---------------------------------------------------------------- - REAL :: thl1,sqv1,sqc1,exner1,sqvg,sqcg,vv,ww - REAL, DIMENSION(its:ite) :: & ZA, & !Height of lowest 1/2 sigma level(m) ZA2, & !Height of 2nd lowest 1/2 sigma level(m) @@ -592,146 +603,349 @@ SUBROUTINE SFCLAY1D_mynn( & TH1D, & !Theta at lowest 1/2 sigma (K) TC1D, & !T at lowest 1/2 sigma (Celsius) TV1D, & !Tv at lowest 1/2 sigma (K) + RHO1D, & !density at lowest 1/2 sigma level QVSH, & !qv at lowest 1/2 sigma (spec humidity) - PSIH2,PSIM2, & !M-O stability functions at z=2 m - PSIH10,PSIM10, & !M-O stability functions at z=10 m - WSPDI, & - CPM, & - z_t,z_q, & !thermal & moisture roughness lengths - ZNTstoch, & + PSIH2, & !M-O stability functions at z=2 m + PSIM10, & !M-O stability functions at z=10 m + PSIH10, & !M-O stability functions at z=10 m + WSPDI, & GOVRTH, & !g/theta - THGB, & !theta at ground - THVGB, & !theta-v at ground PSFC, & !press at surface (Pa/1000) QSFCMR, & !qv at surface (mixing ratio, kg/kg) - GZ2OZ0, & !LOG((2.0+ZNT(I))/ZNT(I)) - GZ10OZ0, & !LOG((10.+ZNT(I))/ZNT(I)) - GZ2OZt, & !LOG((2.0+z_t(i))/z_t(i)) - GZ10OZt, & !LOG((10.+z_t(i))/z_t(i)) - GZ1OZt !LOG((ZA(I)+z_t(i))/z_t(i)) - - INTEGER :: N,I,K,L,NZOL,NK,NZOL2,NZOL10, ITER, yesno - INTEGER, PARAMETER :: ITMAX=1 - - REAL :: PL,THCON,TVCON,E1 - REAL :: DTHVDZ,DTHVM,VCONV,RZOL,RZOL2,RZOL10,ZOL2,ZOL10 - REAL :: DTG,DTTHX,DTHDZ,PSIT10,PSIQ,PSIQ2,PSIQ10 + THCON, & !conversion from temp to theta + zratio_lnd, zratio_ice, zratio_ocn, & !z0/zt + TSK_lnd, TSK_ice, TSK_ocn, & !absolute temperature + THSK_lnd, THSK_ice, THSK_ocn, & !theta + THVSK_lnd, THVSK_ice, THVSK_ocn, & !theta-v + GZ1OZ0_lnd, GZ1OZ0_ice, GZ1OZ0_ocn, & !LOG((ZA(I)+ZNT(i))/ZNT(i)) + GZ1OZt_lnd, GZ1OZt_ice, GZ1OZt_ocn, & !LOG((ZA(I)+ZT(i))/ZT(i)) + GZ2OZ0_lnd, GZ2OZ0_ice, GZ2OZ0_ocn, & !LOG((2.0+ZNT(I))/ZNT(I)) + GZ2OZt_lnd, GZ2OZt_ice, GZ2OZt_ocn, & !LOG((2.0+ZT(I))/ZT(I)) + GZ10OZ0_lnd, GZ10OZ0_ice, GZ10OZ0_ocn, & !LOG((10.+ZNT(I))/ZNT(I)) + GZ10OZt_lnd, GZ10OZt_ice, GZ10OZt_ocn, & !LOG((10.+ZT(I))/ZT(I)) + ZNTstoch_lnd, ZNTstoch_ice, ZNTstoch_ocn, & + ZT_lnd, ZT_ice, ZT_ocn, & + ZQ_lnd, ZQ_ice, ZQ_ocn, & + PSIQ_lnd, PSIQ_ice, PSIQ_ocn, & + PSIQ2_lnd, PSIQ2_ice, PSIQ2_ocn, & + QSFCMR_lnd, QSFCMR_ice, QSFCMR_ocn + + INTEGER :: N,I,K,L,yesno + + REAL :: PL,E1,TABS + REAL :: WSPD_lnd, WSPD_ice, WSPD_ocn + REAL :: DTHVDZ,DTHVM,VCONV,ZOL2,ZOL10,ZOLZA,ZOLZ0 + REAL :: DTG,DTTHX,PSIQ,PSIQ2,PSIQ10,PSIT10 REAL :: FLUXC,VSGD REAL :: restar,VISC,DQG,OLDUST,OLDTST - REAL, PARAMETER :: psilim = -10. ! ONLY AFFECTS z/L > 2.0 + !------------------------------------------------------------------- + IF (debug_code >= 1) THEN + write(0,*)"ITIMESTEP=",ITIMESTEP," iter=",iter + DO I=its,ite + write(0,*)"=== imortant input to mynnsfclayer, i:", i + IF (dry(i)) THEN + write(0,*)"dry=",dry(i)," pblh=",pblh(i)," tsk=", tskin_lnd(i),& + " tsurf=", tsurf_lnd(i)," qsfc=", qsfc_lnd(i)," znt=", znt_lnd(i),& + " ust=", ust_lnd(i)," snowh=", snowh_lnd(i),"psfcpa=",PSFCPA(i), & + " dz=",dz8w1d(i)," qflx=",qflx(i)," hflx=",hflx(i)," hpbl=",pblh(i) + ENDIF + IF (icy(i)) THEN + write(0,*)"icy=",icy(i)," pblh=",pblh(i)," tsk=", tskin_ice(i),& + " tsurf=", tsurf_ice(i)," qsfc=", qsfc_ice(i)," znt=", znt_ice(i),& + " ust=", ust_ice(i)," snowh=", snowh_ice(i),"psfcpa=",PSFCPA(i), & + " dz=",dz8w1d(i)," qflx=",qflx(i)," hflx=",hflx(i)," hpbl=",pblh(i) + ENDIF + IF (wet(i)) THEN + write(0,*)"wet=",wet(i)," pblh=",pblh(i)," tsk=", tskin_ocn(i),& + " tsurf=", tsurf_ocn(i)," qsfc=", qsfc_ocn(i)," znt=", znt_ocn(i),& + " ust=", ust_ocn(i)," snowh=", snowh_ocn(i),"psfcpa=",PSFCPA(i), & + " dz=",dz8w1d(i)," qflx=",qflx(i)," hflx=",hflx(i)," hpbl=",pblh(i) + ENDIF + ENDDO + ENDIF DO I=its,ite - ! CONVERT GROUND & LOWEST LAYER TEMPERATURE TO POTENTIAL TEMPERATURE: - ! PSFC cmb + ! PSFC ( in cmb) is used later in saturation checks PSFC(I)=PSFCPA(I)/1000. - THGB(I)=TSK(I)*(100./PSFC(I))**ROVCP !(K) - ! PL cmb - PL=P1D(I)/1000. - THCON=(100./PL)**ROVCP - TH1D(I)=T1D(I)*THCON !(Theta, K) + ! DEFINE SKIN TEMPERATURES FOR LAND/WATER/ICE + TSK_lnd(I) = 0.5 * (tsurf_lnd(i)+tskin_lnd(i)) + TSK_ice(I) = 0.5 * (tsurf_ice(i)+tskin_ice(i)) + TSK_ocn(I) = 0.5 * (tsurf_ocn(i)+tskin_ocn(i)) + QVSH(I)=QV1D(I)/(1.+QV1D(I)) !CONVERT TO SPEC HUM (kg/kg) + THCON(I)=(100000./PSFCPA(I))**ROVCP + ENDDO + + DO I=its,ite + ! CONVERT SKIN TEMPERATURES TO POTENTIAL TEMPERATURE: + THSK_lnd(I) = TSK_lnd(I)*THCON(I) !(K) + THSK_ice(I) = TSK_ice(I)*THCON(I) !(K) + THSK_ocn(I) = TSK_ocn(I)*THCON(I) !(K) + ENDDO + + DO I=its,ite + ! CONVERT SKIN POTENTIAL TEMPERATURES TO VIRTUAL POTENTIAL TEMPERATURE: + THVSK_lnd(I) = THSK_lnd(I)*(1.+EP1*QVSH(I)) !(K) + THVSK_ice(I) = THSK_ice(I)*(1.+EP1*QVSH(I)) !(K) + THVSK_ocn(I) = THSK_ocn(I)*(1.+EP1*QVSH(I)) !(K) + ENDDO + + DO I=its,ite + ! CONVERT LOWEST LAYER TEMPERATURE TO POTENTIAL TEMPERATURE: + TH1D(I)=T1D(I)*THCON(I) !(Theta, K) TC1D(I)=T1D(I)-273.15 !(T, Celsius) + ENDDO + DO I=its,ite ! CONVERT TO VIRTUAL TEMPERATURE - QVSH(I)=QV1D(I)/(1.+QV1D(I)) !CONVERT TO SPEC HUM (kg/kg) - TVCON=(1.+EP1*QVSH(I)) - THV1D(I)=TH1D(I)*TVCON !(K) - TV1D(I)=T1D(I)*TVCON !(K) + THV1D(I)=TH1D(I)*(1.+EP1*QVSH(I)) !(K) + TV1D(I)=T1D(I)*(1.+EP1*QVSH(I)) !(K) + ENDDO - !RHO1D(I)=PSFCPA(I)/(R*TV1D(I)) !now using value calculated in sfc driver + DO I=its,ite + RHO1D(I)=PSFCPA(I)/(R*TV1D(I)) !now using value calculated in sfc driver ZA(I)=0.5*dz8w1d(I) !height of first half-sigma level ZA2(I)=dz8w1d(I) + 0.5*dz2w1d(I) !height of 2nd half-sigma level GOVRTH(I)=G/TH1D(I) ENDDO DO 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))) + QFX(i)=QFLX(i)*RHO1D(I) + HFX(i)=HFLX(i)*RHO1D(I)*cp + ENDDO + + IF (debug_code ==2) THEN + !write(*,*)"ITIMESTEP=",ITIMESTEP + DO I=its,ite + write(*,*)"=== derived quantities in mynn sfc layer, i:", i + write(*,*)" land, ice, water" + write(*,*)"dry=",dry(i)," icy=",icy(i)," wet=",wet(i) + write(*,*)"tsk=", tsk_lnd(i),tsk_ice(i),tsk_ocn(i) + write(*,*)"thvsk=", thvsk_lnd(i),thvsk_ice(i),thvsk_ocn(i) + write(*,*)"THV1D=", THV1D(i)," TV1D=",TV1D(i) + write(*,*)"RHO1D=", RHO1D(i)," GOVRTH=",GOVRTH(i) + ENDDO + ENDIF + + DO I=its,ite + + IF (ITIMESTEP == 1) THEN + IF (wet(i)) THEN + IF (TSK_ocn(I) .LT. 273.15) THEN + !SATURATION VAPOR PRESSURE WRT ICE (SVP1=.6112; 10*mb) + E1=SVP1*EXP(4648*(1./273.15 - 1./TSK_ocn(I)) - & + & 11.64*LOG(273.15/TSK_ocn(I)) + 0.02265*(273.15 - TSK_ocn(I))) + ELSE + !SATURATION VAPOR PRESSURE WRT WATER (Bolton 1980) + E1=SVP1*EXP(SVP2*(TSK_ocn(I)-SVPT0)/(TSK_ocn(i)-SVP3)) + ENDIF + QSFC_ocn(I)=EP2*E1/(PSFC(I)-ep_3*E1) !specific humidity + QSFCMR_ocn(I)=EP2*E1/(PSFC(I)-E1) !mixing ratio + ENDIF + IF (dry(i)) THEN + TABS = 0.5*(TSK_lnd(I) + T1D(I)) + IF (TABS .LT. 273.15) THEN + !SATURATION VAPOR PRESSURE WRT ICE (SVP1=.6112; 10*mb) + E1=SVP1*EXP(4648*(1./273.15 - 1./TABS) - & + & 11.64*LOG(273.15/TABS) + 0.02265*(273.15 - TABS)) + ELSE + !SATURATION VAPOR PRESSURE WRT WATER (Bolton 1980) + E1=SVP1*EXP(SVP2*(TABS-SVPT0)/(TABS-SVP3)) + ENDIF + QSFC_lnd(I)=EP2*E1/(PSFC(I)-ep_3*E1) !specific humidity + QSFC_lnd(I)=0.5*(QSFC_lnd(I) + QSFC(I)) + QSFCMR_lnd(I)=QSFC_lnd(I)/(1.-QSFC_lnd(I)) !mixing ratio + ENDIF + IF (icy(i)) THEN + IF (TSK_ice(I) .LT. 273.15) THEN + !SATURATION VAPOR PRESSURE WRT ICE (SVP1=.6112; 10*mb) + E1=SVP1*EXP(4648*(1./273.15 - 1./TSK_ice(I)) - & + & 11.64*LOG(273.15/TSK_ice(I)) + 0.02265*(273.15 - TSK_ice(I))) + ELSE + !SATURATION VAPOR PRESSURE WRT WATER (Bolton 1980) + E1=SVP1*EXP(SVP2*(TSK_ice(I)-SVPT0)/(TSK_ice(i)-SVP3)) + ENDIF + QSFC_ice(I)=EP2*E1/(PSFC(I)-ep_3*E1) !specific humidity + QSFCMR_ice(I)=EP2*E1/(PSFC(I)-E1) !mixing ratio + ENDIF + ELSE - !SATURATION VAPOR PRESSURE WRT WATER (Bolton 1980) - E1=SVP1*EXP(SVP2*(TSK(I)-SVPT0)/(TSK(I)-SVP3)) - ENDIF - !FOR LAND POINTS, QSFC can come from LSM, ONLY RECOMPUTE OVER WATER - IF (xland(i).gt.1.5 .or. QSFC(i).le.0.0) THEN !WATER - QSFC(I)=EP2*E1/(PSFC(I)-ep_3*E1) !specific humidity - QSFCMR(I)=EP2*E1/(PSFC(I)-E1) !mixing ratio - ELSE !LAND - QSFCMR(I)=QSFC(I)/(1.-QSFC(I)) + + ! Use what comes out of the LSM, NST, and CICE + IF (wet(i)) QSFCMR_ocn(I)=QSFC_ocn(I)/(1.-QSFC_ocn(I)) + IF (dry(i)) QSFCMR_lnd(I)=QSFC_lnd(I)/(1.-QSFC_lnd(I)) + IF (icy(i)) QSFCMR_ice(I)=QSFC_ice(I)/(1.-QSFC_ice(I)) + ENDIF - IF (TSK(I) .LT. 273.15) THEN + ! QGH CHANGED TO USE LOWEST-LEVEL AIR TEMP + ! Q2SAT = QGH IN LSM + IF (T1D(I) .LT. 273.15) THEN !SATURATION VAPOR PRESSURE WRT ICE - E1=SVP1*EXP(4648*(1./273.15 - 1./T1D(I)) - & + 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*QV1D(I)) ENDDO + IF (debug_code == 2) THEN + write(*,*)"ITIMESTEP=",ITIMESTEP + DO I=its,ite + if (wet(i)) then + write(*,*)"==== q-bombs, i:",i," wet" + write(*,*)"QSFC_ocn=", QSFC_ocn(I)," QSFCMR_ocn=", QSFCMR_ocn(I)," QGH=",QGH(I) + endif + if(dry(i)) then + write(*,*)"==== q-bombs, i:",i," dry" + write(*,*)"QSFC_lnd=", QSFC_lnd(I)," QSFCMR_lnd=", QSFCMR_lnd(I)," QGH=",QGH(I) + endif + if(icy(i)) then + write(*,*)"==== q-bombs, i:",i," ice" + write(*,*)"QSFC_ice=", QSFC_ice(I)," QSFCMR_ice=", QSFCMR_ice(I)," QGH=",QGH(I) + endif + ENDDO + ENDIF + DO I=its,ite - WSPD(I)=SQRT(U1D(I)*U1D(I)+V1D(I)*V1D(I)) - - !account for partial condensation - exner1=(p1d(I)/p1000mb)**ROVCP - sqc1=qc1d(I)/(1.+qc1d(I)) !lowest mod level cloud water spec hum - sqv1=QVSH(I) !lowest mod level water vapor spec hum - thl1=TH1D(I)-xlvcp/exner1*sqc1 - sqvg=qsfc(I) !sfc water vapor spec hum - sqcg=qcg(I)/(1.+qcg(I)) !sfc cloud water spec hum - - vv = thl1-THGB(I) - !TGS:ww = mavail(I)*(sqv1-sqvg) + (sqc1-sqcg) - ww = (sqv1-sqvg) + (sqc1-sqcg) - - !TGS:THVGB(I)=THGB(I)*(1.+EP1*QSFC(I)*MAVAIL(I)) - THVGB(I)=THGB(I)*(1.+EP1*QSFC(I)) - - DTHDZ=(TH1D(I)-THGB(I)) - DTHVDZ=(THV1D(I)-THVGB(I)) - !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 - !------------------------------------------------------- - ! Use Beljaars over land and water - fluxc = max(hfx(i)/RHO1D(i)/cp & - & + ep1*THVGB(I)*qfx(i)/RHO1D(i),0.) - WSTAR(I) = vconvc*(g/TSK(i)*pblh(i)*fluxc)**.33 - - !-------------------------------------------------------- - ! Mahrt and Sun low-res correction - ! (for 13 km ~ 0.37 m/s; for 3 km == 0 m/s) - !-------------------------------------------------------- - VSGD = 0.32 * (max(dx(i)/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),-20.0) - BR(I)=MIN(BR(I),2.0) - + ! DH* 20200401 - note. A weird bug in Intel 18 on hera prevents using the + ! normal -O2 optimization in REPRO and PROD mode for this file. Not reproducible + ! by every user, the bug manifests itself in the resulting wind speed WSPD(I) + ! being -99.0 despite the assignments in lines 932 and 933. *DH + WSPD(I)=SQRT(U1D(I)*U1D(I)+V1D(I)*V1D(I)) + WSPD_ocn = -99. + WSPD_ice = -99. + WSPD_lnd = -99. + + IF (wet(i)) THEN + DTHVDZ=(THV1D(I)-THVSK_ocn(I)) + !-------------------------------------------------------- + ! Calculate the convective velocity scale (WSTAR) and + ! subgrid-scale velocity (VSGD) following Beljaars (1995, QJRMS) + ! and Mahrt and Sun (1995, MWR), respectively + !------------------------------------------------------- + fluxc = max(hfx(i)/RHO1D(i)/cp & + & + ep1*THVSK_ocn(I)*qfx(i)/RHO1D(i),0.) + !WSTAR(I) = vconvc*(g/TSK(i)*pblh(i)*fluxc)**onethird + WSTAR(I) = vconvc*(g/TSK_ocn(i)*pblh(i)*fluxc)**onethird + !-------------------------------------------------------- + ! Mahrt and Sun low-res correction - modified for water points (halved) + ! (for 13 km ~ 0.18 m/s; for 3 km == 0 m/s) + !-------------------------------------------------------- + VSGD = MIN( 0.16 * (max(dx(i)/5000.-1.,0.))**onethird , 0.25) + WSPD_ocn=SQRT(WSPD(I)*WSPD(I)+WSTAR(I)*WSTAR(I)+vsgd*vsgd) + WSPD_ocn=MAX(WSPD_ocn,wmin) + !-------------------------------------------------------- + ! CALCULATE THE BULK RICHARDSON NUMBER OF SURFACE LAYER, + ! ACCORDING TO AKB(1976), EQ(12). + !-------------------------------------------------------- + rb_ocn(I)=GOVRTH(I)*ZA(I)*DTHVDZ/(WSPD_ocn*WSPD_ocn) + IF (ITIMESTEP == 1) THEN + rb_ocn(I)=MAX(rb_ocn(I),-2.0) + rb_ocn(I)=MIN(rb_ocn(I), 2.0) + ELSE + rb_ocn(I)=MAX(rb_ocn(I),-50.0) + rb_ocn(I)=MIN(rb_ocn(I), 50.0) + ENDIF + ENDIF ! end water point + + IF (dry(i)) THEN + DTHVDZ=(THV1D(I)-THVSK_lnd(I)) + !-------------------------------------------------------- + ! Calculate the convective velocity scale (WSTAR) and + ! subgrid-scale velocity (VSGD) following Beljaars (1995, QJRMS) + ! and Mahrt and Sun (1995, MWR), respectively + !------------------------------------------------------- + fluxc = max(hfx(i)/RHO1D(i)/cp & + & + ep1*THVSK_lnd(I)*qfx(i)/RHO1D(i),0.) + !WSTAR(I) = vconvc*(g/TSK(i)*pblh(i)*fluxc)**onethird + !increase height scale, assuming that the non-local transoport + !from the mass-flux (plume) mixing exceedsd the PBLH. + WSTAR(I) = vconvc*(g/TSK_lnd(i)*MIN(1.5*pblh(i),4000.)*fluxc)**onethird + !-------------------------------------------------------- + ! Mahrt and Sun low-res correction + ! (for 13 km ~ 0.37 m/s; for 3 km == 0 m/s) + !-------------------------------------------------------- + VSGD = MIN( 0.32 * (max(dx(i)/5000.-1.,0.))**onethird , 0.5) + WSPD_lnd=SQRT(WSPD(I)*WSPD(I)+WSTAR(I)*WSTAR(I)+vsgd*vsgd) + WSPD_lnd=MAX(WSPD_lnd,wmin) + !-------------------------------------------------------- + ! CALCULATE THE BULK RICHARDSON NUMBER OF SURFACE LAYER, + ! ACCORDING TO AKB(1976), EQ(12). + !-------------------------------------------------------- + rb_lnd(I)=GOVRTH(I)*ZA(I)*DTHVDZ/(WSPD_lnd*WSPD_lnd) + !From Tilden Meyers: + !IF (rb_lnd(I) .GE 0.0) THEN + ! ust_lnd(i)=WSPD_lnd*0.1/(1.0 + 10.0*rb_lnd(I)) + !ELSE + ! ust_lnd(i)=WSPD_lnd*0.1*(1.0 - 10.0*rb_lnd(I))**onethird + !ENDIF + IF (ITIMESTEP == 1) THEN + rb_lnd(I)=MAX(rb_lnd(I),-2.0) + rb_lnd(I)=MIN(rb_lnd(I), 2.0) + ELSE + rb_lnd(I)=MAX(rb_lnd(I),-50.0) + rb_lnd(I)=MIN(rb_lnd(I), 50.0) + ENDIF + ENDIF ! end land point + + IF (icy(i)) THEN + DTHVDZ=(THV1D(I)-THVSK_ice(I)) + !-------------------------------------------------------- + ! Calculate the convective velocity scale (WSTAR) and + ! subgrid-scale velocity (VSGD) following Beljaars (1995, QJRMS) + ! and Mahrt and Sun (1995, MWR), respectively + !------------------------------------------------------- + fluxc = max(hfx(i)/RHO1D(i)/cp & + & + ep1*THVSK_ice(I)*qfx(i)/RHO1D(i),0.) + !WSTAR(I) = vconvc*(g/TSK(i)*pblh(i)*fluxc)**onethird + !increase height scale, assuming that the non-local transport + !from the mass-flux (plume) mixing exceedsd the PBLH. + WSTAR(I) = vconvc*(g/TSK_ice(i)*MIN(1.5*pblh(i),4000.)*fluxc)**onethird + !-------------------------------------------------------- + ! Mahrt and Sun low-res correction + ! (for 13 km ~ 0.37 m/s; for 3 km == 0 m/s) + !-------------------------------------------------------- + VSGD = MIN( 0.32 * (max(dx(i)/5000.-1.,0.))**onethird , 0.5) + WSPD_ice=SQRT(WSPD(I)*WSPD(I)+WSTAR(I)*WSTAR(I)+vsgd*vsgd) + WSPD_ice=MAX(WSPD_ice,wmin) + !-------------------------------------------------------- + ! CALCULATE THE BULK RICHARDSON NUMBER OF SURFACE LAYER, + ! ACCORDING TO AKB(1976), EQ(12). + !-------------------------------------------------------- + rb_ice(I)=GOVRTH(I)*ZA(I)*DTHVDZ/(WSPD_ice*WSPD_ice) + IF (ITIMESTEP == 1) THEN + rb_ice(I)=MAX(rb_ice(I),-2.0) + rb_ice(I)=MIN(rb_ice(I), 2.0) + ELSE + rb_ice(I)=MAX(rb_ice(I),-50.0) + rb_ice(I)=MIN(rb_ice(I), 50.0) + ENDIF + ENDIF ! end ice point + + !NOW CONDENSE THE POSSIBLE WSPD VALUES BY TAKING THE MAXIMUM + WSPD(I) = MAX(WSPD_ice,WSPD_ocn) + WSPD(I) = MAX(WSPD_lnd,WSPD(I)) + + IF (debug_code >= 1) THEN + write(*,*)"===== After rb calc in mynn sfc layer:" + write(*,*)"ITIMESTEP=",ITIMESTEP + write(*,*)"WSPD=", WSPD(I)," WSTAR=", WSTAR(I)," vsgd=",vsgd + IF (icy(i))write(*,*)"rb_ice=", rb_ice(I)," DTHVDZ=",DTHVDZ + IF (wet(i))write(*,*)"rb_ocn=", rb_ocn(I)," DTHVDZ=",DTHVDZ + IF (dry(i))write(*,*)"rb_lnd=", rb_lnd(I)," DTHVDZ=",DTHVDZ + ENDIF + ! 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) @@ -739,626 +953,1125 @@ SUBROUTINE SFCLAY1D_mynn( & !-------------------------------------------------------------------- !-------------------------------------------------------------------- -!--- BEGIN ITERATION LOOP (ITMAX=5); USUALLY CONVERGES IN TWO PASSES +!--- BEGIN I-LOOP !-------------------------------------------------------------------- !-------------------------------------------------------------------- DO I=its,ite - ITER = 1 - DO WHILE (ITER .LE. ITMAX) - - !COMPUTE KINEMATIC VISCOSITY (m2/s) Andreas (1989) CRREL Rep. 89-11 - !valid between -173 and 277 degrees C. - VISC=1.326e-5*(1. + 6.542e-3*TC1D(I) + 8.301e-6*TC1D(I)*TC1D(I) & - - 4.84e-9*TC1D(I)*TC1D(I)*TC1D(I)) - - IF((XLAND(I)-1.5).GE.0)THEN - !-------------------------------------- - ! WATER - !-------------------------------------- - ! CALCULATE z0 (znt) - !-------------------------------------- - IF ( PRESENT(ISFTCFLX) ) THEN - IF ( ISFTCFLX .EQ. 0 ) THEN - IF (COARE_OPT .EQ. 3.0) THEN - !COARE 3.0 (MISLEADING SUBROUTINE NAME) - CALL charnock_1955(ZNT(i),UST(i),WSPD(i),visc,ZA(I)) - ELSE - !COARE 3.5 - CALL edson_etal_2013(ZNT(i),UST(i),WSPD(i),visc,ZA(I)) - ENDIF - 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 - IF (COARE_OPT .EQ. 3.0) THEN - !COARE 3.0 (MISLEADING SUBROUTINE NAME) - CALL charnock_1955(ZNT(i),UST(i),WSPD(i),visc,ZA(I)) - ELSE - !COARE 3.5 - CALL edson_etal_2013(ZNT(i),UST(i),WSPD(i),visc,ZA(I)) - ENDIF + !COMPUTE KINEMATIC VISCOSITY (m2/s) Andreas (1989) CRREL Rep. 89-11 + !valid between -173 and 277 degrees C. + VISC=1.326e-5*(1. + 6.542e-3*TC1D(I) + 8.301e-6*TC1D(I)*TC1D(I) & + - 4.84e-9*TC1D(I)*TC1D(I)*TC1D(I)) + + IF (wet(i)) THEN + !-------------------------------------- + ! WATER + !-------------------------------------- + ! CALCULATE z0 (znt) + !-------------------------------------- + IF (debug_code >= 1) THEN + write(*,*)"=============Input to ZNT over water:" + write(*,*)"u*:",UST_ocn(i)," wspd=",WSPD(i)," visc=",visc," za=",ZA(I) + ENDIF + IF ( PRESENT(ISFTCFLX) ) THEN + IF ( ISFTCFLX .EQ. 0 ) THEN + IF (COARE_OPT .EQ. 3.0) THEN + !COARE 3.0 (MISLEADING SUBROUTINE NAME) + CALL charnock_1955(ZNT_ocn(i),UST_ocn(i),WSPD(i),visc,ZA(I)) + ELSE + !COARE 3.5 + CALL edson_etal_2013(ZNT_ocn(i),UST_ocn(i),WSPD(i),visc,ZA(I)) ENDIF - ELSE - !DEFAULT TO COARE 3.0/3.5 + ELSEIF ( ISFTCFLX .EQ. 1 .OR. ISFTCFLX .EQ. 2 ) THEN + CALL davis_etal_2008(ZNT_ocn(i),UST_ocn(i)) + ELSEIF ( ISFTCFLX .EQ. 3 ) THEN + CALL Taylor_Yelland_2001(ZNT_ocn(i),UST_ocn(i),WSPD(i)) + ELSEIF ( ISFTCFLX .EQ. 4 ) THEN IF (COARE_OPT .EQ. 3.0) THEN - !COARE 3.0 - CALL charnock_1955(ZNT(i),UST(i),WSPD(i),visc,ZA(I)) + !COARE 3.0 (MISLEADING SUBROUTINE NAME) + CALL charnock_1955(ZNT_ocn(i),UST_ocn(i),WSPD(i),visc,ZA(I)) ELSE !COARE 3.5 - CALL edson_etal_2013(ZNT(i),UST(i),WSPD(i),visc,ZA(I)) + CALL edson_etal_2013(ZNT_ocn(i),UST_ocn(i),WSPD(i),visc,ZA(I)) ENDIF ENDIF + ELSE + !DEFAULT TO COARE 3.0/3.5 + IF (COARE_OPT .EQ. 3.0) THEN + !COARE 3.0 + CALL charnock_1955(ZNT_ocn(i),UST_ocn(i),WSPD(i),visc,ZA(I)) + ELSE + !COARE 3.5 + CALL edson_etal_2013(ZNT_ocn(i),UST_ocn(i),WSPD(i),visc,ZA(I)) + ENDIF + ENDIF - ! add stochastic perturbaction of ZNT - if (spp_pbl==1) then - ZNTstoch(I) = MAX(ZNT(I) + 1.5 * ZNT(I) * rstoch1D(i), 1e-6) - else - ZNTstoch(I) = ZNT(I) - endif + ! add stochastic perturbation of ZNT + if (spp_pbl==1) then + ZNTstoch_ocn(I) = MAX(ZNT_ocn(I) + ZNT_ocn(I)*1.0*rstoch1D(i), 1e-6) + else + ZNTstoch_ocn(I) = ZNT_ocn(I) + 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)*ZNTstoch(i)/visc, 0.1) - - !-------------------------------------- - !CALCULATE z_t and z_q - !-------------------------------------- - IF ( PRESENT(ISFTCFLX) ) THEN - IF ( ISFTCFLX .EQ. 0 ) THEN - IF (COARE_OPT .EQ. 3.0) THEN - CALL fairall_etal_2003(z_t(i),z_q(i),restar,UST(i),visc) - ELSE - !presumably, this will be published soon, but hasn't yet - CALL fairall_etal_2014(z_t(i),z_q(i),restar,UST(i),visc,rstoch1D(i),spp_pbl) - ENDIF - ELSEIF ( ISFTCFLX .EQ. 1 ) THEN - IF (COARE_OPT .EQ. 3.0) THEN - CALL fairall_etal_2003(z_t(i),z_q(i),restar,UST(i),visc) - ELSE - CALL fairall_etal_2014(z_t(i),z_q(i),restar,UST(i),visc,rstoch1D(i),spp_pbl) - ENDIF - ELSEIF ( ISFTCFLX .EQ. 2 ) THEN - CALL garratt_1992(z_t(i),z_q(i),ZNTstoch(i),restar,XLAND(I)) - ELSEIF ( ISFTCFLX .EQ. 3 ) THEN - IF (COARE_OPT .EQ. 3.0) THEN - CALL fairall_etal_2003(z_t(i),z_q(i),restar,UST(i),visc) - ELSE - CALL fairall_etal_2014(z_t(i),z_q(i),restar,UST(i),visc,rstoch1D(i),spp_pbl) - ENDIF - ELSEIF ( ISFTCFLX .EQ. 4 ) THEN - CALL zilitinkevich_1995(ZNTstoch(i),z_t(i),z_q(i),restar,& - UST(I),KARMAN,XLAND(I),IZ0TLND,spp_pbl,rstoch1D(i)) + IF (debug_code >= 1) THEN + write(*,*)"==========Output ZNT over water:" + write(*,*)"ZNT:",ZNTstoch_ocn(i) + 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_ocn(i)*ZNTstoch_ocn(i)/visc, 0.1) + + !-------------------------------------- + !CALCULATE z_t and z_q + !-------------------------------------- + IF (debug_code >= 1) THEN + write(*,*)"=============Input to ZT over water:" + write(*,*)"u*:",UST_ocn(i)," restar=",restar," visc=",visc + ENDIF + + IF ( PRESENT(ISFTCFLX) ) THEN + IF ( ISFTCFLX .EQ. 0 ) THEN + IF (COARE_OPT .EQ. 3.0) THEN + CALL fairall_etal_2003(ZT_ocn(i),ZQ_ocn(i),restar,UST_ocn(i),visc,& + rstoch1D(i),spp_pbl) + ELSE + !presumably, this will be published soon, but hasn't yet + CALL fairall_etal_2014(ZT_ocn(i),ZQ_ocn(i),restar,UST_ocn(i),visc,& + rstoch1D(i),spp_pbl) ENDIF - ELSE - !DEFAULT TO COARE 3.0/3.5 + ELSEIF ( ISFTCFLX .EQ. 1 ) THEN IF (COARE_OPT .EQ. 3.0) THEN - CALL fairall_etal_2003(z_t(i),z_q(i),restar,UST(i),visc) + CALL fairall_etal_2003(ZT_ocn(i),ZQ_ocn(i),restar,UST_ocn(i),visc,& + rstoch1D(i),spp_pbl) ELSE - CALL fairall_etal_2014(z_t(i),z_q(i),restar,UST(i),visc,rstoch1D(i),spp_pbl) + CALL fairall_etal_2014(ZT_ocn(i),ZQ_ocn(i),restar,UST_ocn(i),visc,& + rstoch1D(i),spp_pbl) + ENDIF + ELSEIF ( ISFTCFLX .EQ. 2 ) THEN + CALL garratt_1992(ZT_ocn(i),ZQ_ocn(i),ZNTstoch_ocn(i),restar,2.0) + ELSEIF ( ISFTCFLX .EQ. 3 ) THEN + IF (COARE_OPT .EQ. 3.0) THEN + CALL fairall_etal_2003(ZT_ocn(i),ZQ_ocn(i),restar,UST_ocn(i),visc,& + rstoch1D(i),spp_pbl) + ELSE + CALL fairall_etal_2014(ZT_ocn(i),ZQ_ocn(i),restar,UST_ocn(i),visc,& + rstoch1D(i),spp_pbl) ENDIF ENDIF - ELSE + !DEFAULT TO COARE 3.0/3.5 + IF (COARE_OPT .EQ. 3.0) THEN + CALL fairall_etal_2003(ZT_ocn(i),ZQ_ocn(i),restar,UST_ocn(i),visc,& + rstoch1D(i),spp_pbl) + ELSE + CALL fairall_etal_2014(ZT_ocn(i),ZQ_ocn(i),restar,UST_ocn(i),visc,& + rstoch1D(i),spp_pbl) + ENDIF + ENDIF + IF (debug_code >= 1) THEN + write(*,*)"=============Output ZT & ZQ over water:" + write(*,*)"ZT:",ZT_ocn(i)," ZQ:",ZQ_ocn(i) + ENDIF - ! add stochastic perturbaction of ZNT - if (spp_pbl==1) then - ZNTstoch(I) = MAX(ZNT(I) + 1.5 * ZNT(I) * rstoch1D(i), 1e-6) - else - ZNTstoch(I) = ZNT(I) - endif + GZ1OZ0_ocn(I)= LOG((ZA(I)+ZNTstoch_ocn(I))/ZNTstoch_ocn(I)) + GZ1OZt_ocn(I)= LOG((ZA(I)+ZT_ocn(i))/ZT_ocn(i)) + GZ2OZ0_ocn(I)= LOG((2.0+ZNTstoch_ocn(I))/ZNTstoch_ocn(I)) + GZ2OZt_ocn(I)= LOG((2.0+ZT_ocn(i))/ZT_ocn(i)) + GZ10OZ0_ocn(I)=LOG((10.+ZNTstoch_ocn(I))/ZNTstoch_ocn(I)) + GZ10OZt_ocn(I)=LOG((10.+ZT_ocn(i))/ZT_ocn(i)) + zratio_ocn(i)=ZNTstoch_ocn(I)/ZT_ocn(I) !need estimate for Li et al. + + ENDIF !end water point + + IF (dry(I)) THEN + + ! add stochastic perturbaction of ZNT + if (spp_pbl==1) then + ZNTstoch_lnd(I) = MAX(ZNT_lnd(I) + ZNT_lnd(I)*1.0*rstoch1D(i), 1e-6) + else + ZNTstoch_lnd(I) = ZNT_lnd(I) + endif + + !-------------------------------------- + ! LAND + !-------------------------------------- + !COMPUTE ROUGHNESS REYNOLDS NUMBER (restar) USING DEFAULT ZNT + restar=MAX(ust_lnd(i)*ZNTstoch_lnd(i)/visc, 0.1) + + !-------------------------------------- + !GET z_t and z_q + !-------------------------------------- + IF (snowh_lnd(i) > 50.) THEN ! (mm) Treat as snow cover - use Andreas + CALL Andreas_2002(ZNTstoch_lnd(i),visc,ust_lnd(i),ZT_lnd(i),ZQ_lnd(i)) + ELSE + IF ( PRESENT(IZ0TLND) ) THEN + IF ( IZ0TLND .LE. 1 ) THEN + CALL zilitinkevich_1995(ZNTstoch_lnd(i),ZT_lnd(i),ZQ_lnd(i),restar,& + UST_lnd(I),KARMAN,1.0,IZ0TLND,spp_pbl,rstoch1D(i)) + ELSEIF ( IZ0TLND .EQ. 2 ) THEN + CALL Yang_2008(ZNTSTOCH_lnd(i),ZT_lnd(i),ZQ_lnd(i),UST_lnd(i),MOL(I),& + qstar(I),restar,visc) + ELSEIF ( IZ0TLND .EQ. 3 ) THEN + !Original MYNN in WRF-ARW used this form: + CALL garratt_1992(ZT_lnd(i),ZQ_lnd(i),ZNTSTOCH_lnd(i),restar,1.0) + ENDIF + ELSE + !DEFAULT TO ZILITINKEVICH + CALL zilitinkevich_1995(ZNTSTOCH_lnd(i),ZT_lnd(i),ZQ_lnd(i),restar,& + UST_lnd(I),KARMAN,1.0,0,spp_pbl,rstoch1D(i)) + ENDIF + ENDIF - !-------------------------------------- - ! LAND - !-------------------------------------- - !COMPUTE ROUGHNESS REYNOLDS NUMBER (restar) USING DEFAULT ZNT - restar=MAX(ust(i)*ZNTstoch(i)/visc, 0.1) - - !-------------------------------------- - !GET z_t and z_q - !-------------------------------------- - !CHECK FOR SNOW/ICE POINTS OVER LAND - !IF ( ZNTSTOCH(i) .LE. SNOWZ0 .AND. TSK(I) .LE. 273.15 ) THEN - IF ( SNOWH(i) .GE. 0.1) THEN - CALL Andreas_2002(ZNTSTOCH(i),visc,ust(i),z_t(i),z_q(i)) + GZ1OZ0_lnd(I)= LOG((ZA(I)+ZNTstoch_lnd(I))/ZNTstoch_lnd(I)) + GZ1OZt_lnd(I)= LOG((ZA(I)+ZT_lnd(i))/ZT_lnd(i)) + GZ2OZ0_lnd(I)= LOG((2.0+ZNTstoch_lnd(I))/ZNTstoch_lnd(I)) + GZ2OZt_lnd(I)= LOG((2.0+ZT_lnd(i))/ZT_lnd(i)) + GZ10OZ0_lnd(I)=LOG((10.+ZNTstoch_lnd(I))/ZNTstoch_lnd(I)) + GZ10OZt_lnd(I)=LOG((10.+ZT_lnd(i))/ZT_lnd(i)) + zratio_lnd(i)=ZNTstoch_lnd(I)/ZT_lnd(I) !need estimate for Li et al. + + ENDIF !end land point + + IF (icy(I)) THEN + + ! add stochastic perturbaction of ZNT + if (spp_pbl==1) then + ZNTstoch_ice(I) = MAX(ZNT_ice(I) + ZNT_ice(I)*1.0*rstoch1D(i), 1e-6) + else + ZNTstoch_ice(I) = ZNT_ice(I) + endif + + !-------------------------------------- + ! ICE + !-------------------------------------- + !COMPUTE ROUGHNESS REYNOLDS NUMBER (restar) USING DEFAULT ZNT + restar=MAX(ust_ice(i)*ZNTstoch_ice(i)/visc, 0.1) + !-------------------------------------- + !GET z_t and z_q + !-------------------------------------- + CALL Andreas_2002(ZNTstoch_ice(i),visc,ust_ice(i),ZT_ice(i),ZQ_ice(i)) + + GZ1OZ0_ice(I)= LOG((ZA(I)+ZNTstoch_ice(I))/ZNTstoch_ice(I)) + GZ1OZt_ice(I)= LOG((ZA(I)+ZT_ice(i))/ZT_ice(i)) + GZ2OZ0_ice(I)= LOG((2.0+ZNTstoch_ice(I))/ZNTstoch_ice(I)) + GZ2OZt_ice(I)= LOG((2.0+ZT_ice(i))/ZT_ice(i)) + GZ10OZ0_ice(I)=LOG((10.+ZNTstoch_ice(I))/ZNTstoch_ice(I)) + GZ10OZt_ice(I)=LOG((10.+ZT_ice(i))/ZT_ice(i)) + zratio_ice(i)=ZNTstoch_ice(I)/ZT_ice(I) !need estimate for Li et al. + + ENDIF !end ice point + + !Capture a representative ZNT + IF (dry(i)) THEN + ZNT(i)=ZNTstoch_lnd(I) + ELSEIF (wet(i)) THEN + ZNT(i)=ZNTstoch_ocn(I) + ELSEIF (icy(i)) THEN + ZNT(i)=ZNTstoch_ice(I) + ENDIF + + !-------------------------------------------------------------------- + !--- DIAGNOSE STABILITY FUNCTIONS FOR THE APPROPRIATE STABILITY CLASS: + ! THE STABILITY CLASSES ARE DETERMINED BY THE BULK RICHARDSON NUMBER. + !-------------------------------------------------------------------- + + IF (wet(i)) THEN + IF (rb_ocn(I) .GT. 0.0) THEN + + !COMPUTE z/L first guess: + IF (itimestep .LE. 1) THEN + CALL Li_etal_2010(ZOL(I),rb_ocn(I),ZA(I)/ZNTstoch_ocn(I),zratio_ocn(I)) ELSE - IF ( PRESENT(IZ0TLND) ) THEN - IF ( IZ0TLND .LE. 1 .OR. IZ0TLND .EQ. 4) THEN - !IF IZ0TLND==4, THEN PSIQ WILL BE RECALCULATED USING - !PAN ET AL (1994), but PSIT FROM ZILI WILL BE USED. - CALL zilitinkevich_1995(ZNTSTOCH(i),z_t(i),z_q(i),restar,& - UST(I),KARMAN,XLAND(I),IZ0TLND,spp_pbl,rstoch1D(i)) - ELSEIF ( IZ0TLND .EQ. 2 ) THEN - CALL Yang_2008(ZNTSTOCH(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),ZNTSTOCH(i),restar,XLAND(I)) - ENDIF - ELSE - !DEFAULT TO ZILITINKEVICH - CALL zilitinkevich_1995(ZNTSTOCH(i),z_t(i),z_q(i),restar,& - UST(I),KARMAN,XLAND(I),0,spp_pbl,rstoch1D(i)) - ENDIF + ZOL(I)=ZA(I)*KARMAN*G*MOL(I)/(TH1D(I)*MAX(UST_ocn(I)*UST_ocn(I),0.0001)) + ZOL(I)=MAX(ZOL(I),0.0) + ZOL(I)=MIN(ZOL(I),50.) + ENDIF + + IF (debug_code >= 1) THEN + IF (ZNTstoch_ocn(i) < 1E-8 .OR. Zt_ocn(i) < 1E-10) THEN + write(0,*)"===(wet) capture bad input in mynn sfc layer, i=:",i + write(0,*)"rb=", rb_ocn(I)," ZNT=", ZNTstoch_ocn(i)," ZT=",Zt_ocn(i) + write(0,*)" tsk=", tskin_ocn(i)," prev z/L=",ZOL(I),& + " tsurf=", tsurf_ocn(i)," qsfc=", qsfc_ocn(i)," znt=", znt_ocn(i),& + " ust=", ust_ocn(i)," snowh=", snowh_ocn(i),"psfcpa=",PSFCPA(i), & + " dz=",dz8w1d(i)," qflx=",qflx(i)," hflx=",hflx(i)," hpbl=",pblh(i) + ENDIF + ENDIF + !Use Pedros iterative function to find z/L + zol(I)=zolri(rb_ocn(I),ZA(I),ZNTstoch_ocn(I),ZT_ocn(I),ZOL(I)) + ZOL(I)=MAX(ZOL(I),0.0) + ZOL(I)=MIN(ZOL(I),50.) + + zolz0 = zol(I)*ZNTstoch_ocn(I)/ZA(I) ! z0/L + zolza = zol(I)*(za(I)+ZNTstoch_ocn(I))/za(I) ! (z+z0/L + zol10 = zol(I)*(10.+ZNTstoch_ocn(I))/za(I) ! (10+z0)/L + zol2 = zol(I)*(2.+ZNTstoch_ocn(I))/za(I) ! (2+z0)/L + + !COMPUTE PSIM and PSIH + !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),ZT_ocn(I),ZNTstoch_ocn(I),ZA(I)) + !CALL PSI_CB2005(PSIM(I),PSIH(I),zolza,zolz0) + ! or use tables + psim(I)=psim_stable(zolza)-psim_stable(zolz0) + psih(I)=psih_stable(zolza)-psih_stable(zolz0) + psim10(I)=psim_stable(zol10)-psim_stable(zolz0) + psih10(I)=psih_stable(zol10)-psih_stable(zolz0) + psih2(I)=psih_stable(zol2)-psih_stable(zolz0) + + ! 1.0 over Monin-Obukhov length + RMOL(I)= ZOL(I)/ZA(I) + + ELSEIF(rb_ocn(I) .EQ. 0.) THEN + !========================================================= + !-----CLASS 3; FORCED CONVECTION/NEUTRAL: + !========================================================= + + PSIM(I)=0.0 + PSIH(I)=PSIM(I) + PSIM10(I)=0. + PSIH10(I)=0. + PSIH2(I)=0. + + ZOL(I) =0. + RMOL(I) =0. + + ELSEIF(rb_ocn(I) .LT. 0.)THEN + !========================================================== + !-----CLASS 4; FREE CONVECTION: + !========================================================== + + !COMPUTE z/L first guess: + IF (itimestep .LE. 1) THEN + CALL Li_etal_2010(ZOL(I),rb_ocn(I),ZA(I)/ZNTstoch_ocn(I),zratio_ocn(I)) + ELSE + ZOL(I)=ZA(I)*KARMAN*G*MOL(I)/(TH1D(I)*MAX(UST_ocn(I)*UST_ocn(I),0.001)) + ZOL(I)=MAX(ZOL(I),-50.0) + ZOL(I)=MIN(ZOL(I),0.0) + ENDIF + + IF (debug_code >= 1) THEN + IF (ZNTstoch_ocn(i) < 1E-8 .OR. Zt_ocn(i) < 1E-10) THEN + write(0,*)"===(wet) capture bad input in mynn sfc layer, i=:",i + write(0,*)"rb=", rb_ocn(I)," ZNT=", ZNTstoch_ocn(i)," ZT=",Zt_ocn(i) + write(0,*)" tsk=", tskin_ocn(i)," wstar=",wstar(i)," prev z/L=",ZOL(I),& + " tsurf=", tsurf_ocn(i)," qsfc=", qsfc_ocn(i)," znt=", znt_ocn(i),& + " ust=", ust_ocn(i)," snowh=", snowh_ocn(i),"psfcpa=",PSFCPA(i), & + " dz=",dz8w1d(i)," qflx=",qflx(i)," hflx=",hflx(i)," hpbl=",pblh(i) + ENDIF ENDIF + !Use Pedros iterative function to find z/L + zol(I)=zolri(rb_ocn(I),ZA(I),ZNTstoch_ocn(I),ZT_ocn(I),ZOL(I)) + ZOL(I)=MAX(ZOL(I),-50.0) + ZOL(I)=MIN(ZOL(I),0.0) + + zolz0 = zol(I)*ZNTstoch_ocn(I)/ZA(I) ! z0/L + zolza = zol(I)*(za(I)+ZNTstoch_ocn(I))/za(I) ! (z+z0/L + zol10 = zol(I)*(10.+ZNTstoch_ocn(I))/za(I) ! (10+z0)/L + zol2 = zol(I)*(2.+ZNTstoch_ocn(I))/za(I) ! (2+z0)/L + + !COMPUTE PSIM and PSIH + !CALL PSI_Suselj_Sood_2010(PSIM(I),PSIH(I),ZOL(I)) + !CALL PSI_Hogstrom_1996(PSIM(I),PSIH(I),ZOL(I), ZT_ocn(I), ZNTstoch_ocn(I), ZA(I)) + !CALL PSI_Businger_1971(PSIM(I),PSIH(I),ZOL(I)) + !CALL PSI_DyerHicks(PSIM(I),PSIH(I),ZOL(I),ZT_ocn(I),ZNTstoch_ocn(I),ZA(I)) + ! use tables + psim(I)=psim_unstable(zolza)-psim_unstable(zolz0) + psih(I)=psih_unstable(zolza)-psih_unstable(zolz0) + psim10(I)=psim_unstable(zol10)-psim_unstable(zolz0) + psih10(I)=psih_unstable(zol10)-psih_unstable(zolz0) + psih2(I)=psih_unstable(zol2)-psih_unstable(zolz0) + + !---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*GZ1OZt_ocn(I)) + PSIM(I)=MIN(PSIM(I),0.9*GZ1OZ0_ocn(I)) + PSIH2(I)=MIN(PSIH2(I),0.9*GZ2OZt_ocn(I)) + PSIM10(I)=MIN(PSIM10(I),0.9*GZ10OZ0_ocn(I)) + PSIH10(I)=MIN(PSIH10(I),0.9*GZ10OZt_ocn(I)) + + RMOL(I) = ZOL(I)/ZA(I) ENDIF - zratio(i)=zntstoch(i)/z_t(i) - - !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)/ZNTstoch(I)) - GZ1OZ0(I)= LOG((ZA(I)+ZNTstoch(I))/ZNTstoch(I)) - GZ1OZt(I)= LOG((ZA(I)+z_t(i))/z_t(i)) - GZ2OZ0(I)= LOG((2.0+ZNTstoch(I))/ZNTstoch(I)) - GZ2OZt(I)= LOG((2.0+z_t(i))/z_t(i)) - GZ10OZ0(I)=LOG((10.+ZNTstoch(I))/ZNTstoch(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.0) THEN - IF (BR(I) .GT. 0.2) THEN - !---CLASS 1; STABLE (NIGHTTIME) CONDITIONS: - REGIME(I)=1. - ELSE - !---CLASS 2; DAMPED MECHANICAL TURBULENCE: - REGIME(I)=2. - ENDIF - !COMPUTE z/L - !CALL Li_etal_2010(ZOL(I),BR(I),ZA(I)/ZNTstoch(I),zratio(I)) -! IF (ITER .EQ. 1 .AND. itimestep .LE. 1) THEN - CALL Li_etal_2010(ZOL(I),BR(I),ZA(I)/ZNTstoch(I),zratio(I)) -! ELSE -! ZOL(I)=ZA(I)*KARMAN*G*MOL(I)/(TH1D(I)*MAX(UST(I)*UST(I),0.0001)) -! ZOL(I)=MAX(ZOL(I),0.0) -! ZOL(I)=MIN(ZOL(I),2.) -! 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),ZNTstoch(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),ZNTstoch(I),ZA(I)) - ENDIF - - ! LOWER LIMIT ON PSI IN STABLE CONDITIONS - PSIM(I)=MAX(PSIM(I),psilim) - PSIH(I)=MAX(PSIH(I),psilim) - PSIM10(I)=MAX(10./ZA(I)*PSIM(I), psilim) - PSIH10(I)=MAX(10./ZA(I)*PSIH(I), psilim) - PSIM2(I)=MAX(2./ZA(I)*PSIM(I), psilim) - PSIH2(I)=MAX(2./ZA(I)*PSIH(I), psilim) - ! 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)/(MAX(UST(I)*UST(I),0.001)) - 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)/ZNTstoch(I),zratio(I)) - !IF (ITER .EQ. 1 .AND. itimestep .LE. 1) THEN - CALL Li_etal_2010(ZOL(I),BR(I),ZA(I)/ZNTstoch(I),zratio(I)) - !ELSE - ! ZOL(I)=ZA(I)*KARMAN*G*MOL(I)/(TH1D(I)*MAX(UST(I)*UST(I),0.001)) - ! ZOL(I)=MAX(ZOL(I),-19.999) - ! 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),-19.9999) - ZOL10=MIN(ZOL10,0.) - ZOL10=MAX(ZOL10,-19.9999) - ZOL2=MIN(ZOL2,0.) - ZOL2=MAX(ZOL2,-19.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), ZNTstoch(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),ZNTstoch(I),ZA(I)) - ELSE - ! LAND - !CALL PSI_Hogstrom_1996(PSIM(I),PSIH(I),ZOL(I), z_t(I), ZNTstoch(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),ZNTstoch(I),ZA(I)) - ENDIF - - PSIM10(I)=10./ZA(I)*PSIM(I) - PSIH10(I)=10./ZA(I)*PSIH(I) - PSIM2(I)=2./ZA(I)*PSIM(I) - PSIH2(I)=2./ZA(I)*PSIH(I) - - !---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*GZ1OZt(I)) !JOE: less restricitive over forest/urban. - 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*GZ2OZt(I)) !JOE: less restricitive over forest/urban. - PSIH2(I)=MIN(PSIH2(I),0.9*GZ2OZ0(I)) - PSIM2(I)=MIN(PSIM2(I),0.9*GZ2OZ0(I)) - PSIM10(I)=MIN(PSIM10(I),0.9*GZ10OZ0(I)) - PSIH10(I)=MIN(PSIH10(I),0.9*GZ10OZ0(I)) - - RMOL(I) = ZOL(I)/ZA(I) - - ENDIF - - !------------------------------------------------------------ - !-----COMPUTE THE FRICTIONAL VELOCITY: - !------------------------------------------------------------ - ! ZA(1982) EQS(2.60),(2.61). - PSIX(I)=GZ1OZ0(I)-PSIM(I) - PSIX10(I)=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(I) - !NON-AVERAGED: UST(I)=KARMAN*WSPD(I)/PSIX(I) - - ! Compute u* without vconv for use in HFX calc when isftcflx > 0 - WSPDI(I)=MAX(SQRT(U1D(I)*U1D(I)+V1D(I)*V1D(I)), wmin) - IF ( PRESENT(USTM) ) THEN - USTM(I)=0.5*USTM(I)+0.5*KARMAN*WSPDI(I)/PSIX(I) - ENDIF + ! CALCULATE THE RESISTANCE: + PSIX_ocn(I) =MAX(GZ1OZ0_ocn(I)-PSIM(I) , 1.0) ! = fm + PSIX10_ocn(I)=MAX(GZ10OZ0_ocn(I)-PSIM10(I), 1.0) ! = fm10 + PSIT_ocn(I) =MAX(GZ1OZt_ocn(I)-PSIH(I) , 1.0) ! = fh + PSIT2_ocn(I) =MAX(GZ2OZt_ocn(I)-PSIH2(I) , 1.0) ! = fh2 + PSIQ_ocn(I) =MAX(LOG((ZA(I)+ZQ_ocn(i))/ZQ_ocn(I))-PSIH(I) ,1.0) + PSIQ2_ocn(I) =MAX(LOG((2.0+ZQ_ocn(i))/ZQ_ocn(I))-PSIH2(I) ,1.0) - IF ((XLAND(I)-1.5).LT.0.) THEN !LAND - UST(I)=MAX(UST(I),0.005) !Further relaxing this limit - no need to go lower - !Keep ustm = ust over land. - IF ( PRESENT(USTM) ) USTM(I)=UST(I) - ENDIF + ENDIF ! end water points - !------------------------------------------------------------ - !-----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(I) =MAX(GZ1OZt(I)-PSIH(I) ,1.) - PSIT2(I)=MAX(GZ2OZt(I)-PSIH2(I),1.) - resist(I)=PSIT(I) - logres(I)=GZ1OZt(I) - - PSIQ=MAX(LOG((ZA(I)+z_q(i))/z_q(I))-PSIH(I) ,1.0) - PSIQ2=MAX(LOG((2.0+z_q(i))/z_q(I))-PSIH2(I) ,1.0) - - IF((XLAND(I)-1.5).LT.0)THEN !Land only - IF ( IZ0TLND .EQ. 4 ) THEN - CALL Pan_etal_1994(PSIQ,PSIQ2,UST(I),PSIH(I),PSIH2(I),& - & KARMAN,ZA(I)) - ENDIF - ENDIF + IF (dry(i)) THEN + IF (rb_lnd(I) .GT. 0.0) THEN - !---------------------------------------------------- - !COMPUTE THE TEMPERATURE SCALE (or FRICTION TEMPERATURE, T*) - !---------------------------------------------------- - !DTG=TH1D(I)-THGB(I) !SWITCH TO THETA-V - DTG=THV1D(I)-THVGB(I) - OLDTST=MOL(I) - MOL(I)=KARMAN*DTG/PSIT(I)/PRT - !t_star(I) = -HFX(I)/(UST(I)*CPM(I)*RHO1D(I)) - !t_star(I) = MOL(I) - !---------------------------------------------------- - !COMPUTE THE MOISTURE SCALE (or q*) - DQG=(QVSH(i)-qsfc(i))*1000. !(kg/kg -> g/kg) - qstar(I)=KARMAN*DQG/PSIQ/PRT - - !CHECK FOR CONVERGENCE - IF (ITER .GE. 2) THEN - !IF (ABS(OLDUST-UST(I)) .lt. 0.01) THEN - IF (ABS(OLDTST-MOL(I)) .lt. 0.01) THEN - ITER = ITER+ITMAX - ENDIF + !COMPUTE z/L first guess: + IF (itimestep .LE. 1) THEN + CALL Li_etal_2010(ZOL(I),rb_lnd(I),ZA(I)/ZNTstoch_lnd(I),zratio_lnd(I)) + ELSE + ZOL(I)=ZA(I)*KARMAN*G*MOL(I)/(TH1D(I)*MAX(UST_lnd(I)*UST_lnd(I),0.0001)) + ZOL(I)=MAX(ZOL(I),0.0) + ZOL(I)=MIN(ZOL(I),50.) + ENDIF - !IF () THEN - ! print*,"ITER:",ITER - ! write(*,1001)"REGIME:",REGIME(I)," z/L:",ZOL(I)," U*:",UST(I)," Tstar:",MOL(I) - ! write(*,1002)"PSIM:",PSIM(I)," PSIH:",PSIH(I)," W*:",WSTAR(I)," DTHV:",THV1D(I)-THVGB(I) - ! write(*,1003)"CPM:",CPM(I)," RHO1D:",RHO1D(I)," L:",ZOL(I)/ZA(I)," DTH:",TH1D(I)-THGB(I) - ! write(*,1004)"Z0/Zt:",zratio(I)," Z0:",ZNTstoch(I)," Zt:",z_t(I)," za:",za(I) - ! write(*,1005)"Re:",restar," MAVAIL:",MAVAIL(I)," QSFC(I):",QSFC(I)," QVSH(I):",QVSH(I) - ! print*,"VISC=",VISC," Z0:",ZNTstoch(I)," T1D(i):",T1D(i) - ! write(*,*)"=============================================" - !ENDIF - ENDIF + IF (debug_code >= 1) THEN + IF (ZNTstoch_lnd(i) < 1E-8 .OR. Zt_lnd(i) < 1E-10) THEN + write(0,*)"===(land) capture bad input in mynn sfc layer, i=:",i + write(0,*)"rb=", rb_lnd(I)," ZNT=", ZNTstoch_lnd(i)," ZT=",Zt_lnd(i) + write(0,*)" tsk=", tskin_lnd(i)," prev z/L=",ZOL(I),& + " tsurf=", tsurf_lnd(i)," qsfc=", qsfc_lnd(i)," znt=", znt_lnd(i),& + " ust=", ust_lnd(i)," snowh=", snowh_lnd(i),"psfcpa=",PSFCPA(i), & + " dz=",dz8w1d(i)," qflx=",qflx(i)," hflx=",hflx(i)," hpbl=",pblh(i) + ENDIF + ENDIF + !Use Pedros iterative function to find z/L + zol(I)=zolri(rb_lnd(I),ZA(I),ZNTstoch_lnd(I),ZT_lnd(I),ZOL(I)) + ZOL(I)=MAX(ZOL(I),0.0) + ZOL(I)=MIN(ZOL(I),50.) + + zolz0 = zol(I)*ZNTstoch_lnd(I)/ZA(I) ! z0/L + zolza = zol(I)*(za(I)+ZNTstoch_lnd(I))/za(I) ! (z+z0/L + zol10 = zol(I)*(10.+ZNTstoch_lnd(I))/za(I) ! (10+z0)/L + zol2 = zol(I)*(2.+ZNTstoch_lnd(I))/za(I) ! (2+z0)/L + + !COMPUTE PSIM and PSIH + !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),ZT_lnd(I),ZNTstoch_lnd(I),ZA(I)) + !CALL PSI_CB2005(PSIM(I),PSIH(I),zolza,zolz0) + psim(I)=psim_stable(zolza)-psim_stable(zolz0) + psih(I)=psih_stable(zolza)-psih_stable(zolz0) + psim10(I)=psim_stable(zol10)-psim_stable(zolz0) + psih10(I)=psih_stable(zol10)-psih_stable(zolz0) + psih2(I)=psih_stable(zol2)-psih_stable(zolz0) + + ! 1.0 over Monin-Obukhov length + RMOL(I)= ZOL(I)/ZA(I) + + ELSEIF(rb_lnd(I) .EQ. 0.) THEN + !========================================================= + !-----CLASS 3; FORCED CONVECTION/NEUTRAL: + !========================================================= + + PSIM(I)=0.0 + PSIH(I)=PSIM(I) + PSIM10(I)=0. + PSIH10(I)=0. + PSIH2(I)=0. + + ZOL(I) =0. + RMOL(I) =0. + + ELSEIF(rb_lnd(I) .LT. 0.)THEN + !========================================================== + !-----CLASS 4; FREE CONVECTION: + !========================================================== + + !COMPUTE z/L first guess: + IF (itimestep .LE. 1) THEN + CALL Li_etal_2010(ZOL(I),rb_lnd(I),ZA(I)/ZNTstoch_lnd(I),zratio_lnd(I)) + ELSE + ZOL(I)=ZA(I)*KARMAN*G*MOL(I)/(TH1D(I)*MAX(UST_lnd(I)*UST_lnd(I),0.001)) + ZOL(I)=MAX(ZOL(I),-50.0) + ZOL(I)=MIN(ZOL(I),0.0) + ENDIF + + IF (debug_code >= 1) THEN + IF (ZNTstoch_lnd(i) < 1E-8 .OR. Zt_lnd(i) < 1E-10) THEN + write(0,*)"===(land) capture bad input in mynn sfc layer, i=:",i + write(0,*)"rb=", rb_lnd(I)," ZNT=", ZNTstoch_lnd(i)," ZT=",Zt_lnd(i) + write(0,*)" tsk=", tskin_lnd(i)," wstar=",wstar(i)," prev z/L=",ZOL(I),& + " tsurf=", tsurf_lnd(i)," qsfc=", qsfc_lnd(i)," znt=", znt_lnd(i),& + " ust=", ust_lnd(i)," snowh=", snowh_lnd(i),"psfcpa=",PSFCPA(i), & + " dz=",dz8w1d(i)," qflx=",qflx(i)," hflx=",hflx(i)," hpbl=",pblh(i) + ENDIF + ENDIF + !Use Pedros iterative function to find z/L + zol(I)=zolri(rb_lnd(I),ZA(I),ZNTstoch_lnd(I),ZT_lnd(I),ZOL(I)) + ZOL(I)=MAX(ZOL(I),-50.0) + ZOL(I)=MIN(ZOL(I),0.0) + + zolz0 = zol(I)*ZNTstoch_lnd(I)/ZA(I) ! z0/L + zolza = zol(I)*(za(I)+ZNTstoch_lnd(I))/za(I) ! (z+z0/L + zol10 = zol(I)*(10.+ZNTstoch_lnd(I))/za(I) ! (10+z0)/L + zol2 = zol(I)*(2.+ZNTstoch_lnd(I))/za(I) ! (2+z0)/L + + !COMPUTE PSIM and PSIH + !CALL PSI_Hogstrom_1996(PSIM(I),PSIH(I),ZOL(I), ZT_lnd(I), ZNTstoch_lnd(I), ZA(I)) + !CALL PSI_Businger_1971(PSIM(I),PSIH(I),ZOL(I)) + !CALL PSI_DyerHicks(PSIM(I),PSIH(I),ZOL(I),ZT_lnd(I),ZNTstoch_lnd(I),ZA(I)) + ! use tables + psim(I)=psim_unstable(zolza)-psim_unstable(zolz0) + psih(I)=psih_unstable(zolza)-psih_unstable(zolz0) + psim10(I)=psim_unstable(zol10)-psim_unstable(zolz0) + psih10(I)=psih_unstable(zol10)-psih_unstable(zolz0) + psih2(I)=psih_unstable(zol2)-psih_unstable(zolz0) + + !---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*GZ1OZt_lnd(I)) + PSIM(I)=MIN(PSIM(I),0.9*GZ1OZ0_lnd(I)) + PSIH2(I)=MIN(PSIH2(I),0.9*GZ2OZt_lnd(I)) + PSIM10(I)=MIN(PSIM10(I),0.9*GZ10OZ0_lnd(I)) + PSIH10(I)=MIN(PSIH10(I),0.9*GZ10OZt_lnd(I)) + + RMOL(I) = ZOL(I)/ZA(I) - ITER = ITER + 1 + ENDIF - ENDDO ! end ITERATION-loop + ! CALCULATE THE RESISTANCE: + PSIX_lnd(I) =MAX(GZ1OZ0_lnd(I)-PSIM(I), 1.0) + PSIX10_lnd(I)=MAX(GZ10OZ0_lnd(I)-PSIM10(I), 1.0) + PSIT_lnd(I) =MAX(GZ1OZt_lnd(I)-PSIH(I) , 1.0) + PSIT2_lnd(I) =MAX(GZ2OZt_lnd(I)-PSIH2(I), 1.0) + PSIQ_lnd(I) =MAX(LOG((ZA(I)+ZQ_lnd(i))/ZQ_lnd(I))-PSIH(I) ,1.0) + PSIQ2_lnd(I) =MAX(LOG((2.0+ZQ_lnd(i))/ZQ_lnd(I))-PSIH2(I) ,1.0) - ENDDO ! end i-loop + ENDIF ! end land points - 1000 format(A,F6.1, A,f6.1, A,f5.1, A,f7.1) - 1001 format(A,F2.0, A,f10.4,A,f5.3, A,f11.5) - 1002 format(A,f7.2, A,f7.2, A,f7.2, A,f10.3) - 1003 format(A,f7.2, A,f7.2, A,f10.3,A,f10.3) - 1004 format(A,f11.3,A,f9.7, A,f9.7, A,f6.2, A,f10.3) - 1005 format(A,f9.2,A,f6.4,A,f7.4,A,f7.4) + IF (icy(i)) THEN + IF (rb_ice(I) .GT. 0.0) THEN - !---------------------------------------------------------- - ! COMPUTE SURFACE HEAT AND MOISTURE FLUXES - !---------------------------------------------------------- - DO I=its,ite + !COMPUTE z/L first guess: + IF (itimestep .LE. 1) THEN + CALL Li_etal_2010(ZOL(I),rb_ice(I),ZA(I)/ZNTstoch_ice(I),zratio_ice(I)) + ELSE + ZOL(I)=ZA(I)*KARMAN*G*MOL(I)/(TH1D(I)*MAX(UST_ice(I)*UST_ice(I),0.0001)) + ZOL(I)=MAX(ZOL(I),0.0) + ZOL(I)=MIN(ZOL(I),50.) + ENDIF + + IF (debug_code >= 1) THEN + IF (ZNTstoch_ice(i) < 1E-8 .OR. Zt_ice(i) < 1E-10) THEN + write(0,*)"===(ice) capture bad input in mynn sfc layer, i=:",i + write(0,*)"rb=", rb_ice(I)," ZNT=", ZNTstoch_ice(i)," ZT=",Zt_ice(i) + write(0,*)" tsk=", tskin_ice(i)," prev z/L=",ZOL(I),& + " tsurf=", tsurf_ice(i)," qsfc=", qsfc_ice(i)," znt=", znt_ice(i),& + " ust=", ust_ice(i)," snowh=", snowh_ice(i),"psfcpa=",PSFCPA(i), & + " dz=",dz8w1d(i)," qflx=",qflx(i)," hflx=",hflx(i)," hpbl=",pblh(i) + ENDIF + ENDIF + !Use Pedros iterative function to find z/L + zol(I)=zolri(rb_ice(I),ZA(I),ZNTstoch_ice(I),ZT_ice(I),ZOL(I)) + ZOL(I)=MAX(ZOL(I),0.0) + ZOL(I)=MIN(ZOL(I),50.) + + zolz0 = zol(I)*ZNTstoch_ice(I)/ZA(I) ! z0/L + zolza = zol(I)*(za(I)+ZNTstoch_ice(I))/za(I) ! (z+z0/L + zol10 = zol(I)*(10.+ZNTstoch_ice(I))/za(I) ! (10+z0)/L + zol2 = zol(I)*(2.+ZNTstoch_ice(I))/za(I) ! (2+z0)/L + + !COMPUTE PSIM and PSIH + !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),ZT_ice(I),ZNTstoch_ice(I),ZA(I)) + !CALL PSI_CB2005(PSIM(I),PSIH(I),zolza,zolz0) + psim(I)=psim_stable(zolza)-psim_stable(zolz0) + psih(I)=psih_stable(zolza)-psih_stable(zolz0) + psim10(I)=psim_stable(zol10)-psim_stable(zolz0) + psih10(I)=psih_stable(zol10)-psih_stable(zolz0) + psih2(I)=psih_stable(zol2)-psih_stable(zolz0) + + ! 1.0 over Monin-Obukhov length + RMOL(I)= ZOL(I)/ZA(I) + + ELSEIF(rb_ice(I) .EQ. 0.) THEN + !========================================================= + !-----CLASS 3; FORCED CONVECTION/NEUTRAL: + !========================================================= + + PSIM(I)=0.0 + PSIH(I)=PSIM(I) + PSIM10(I)=0. + PSIH10(I)=0. + PSIH2(I)=0. + + ZOL(I) =0. + RMOL(I) =0. + + ELSEIF(rb_ice(I) .LT. 0.)THEN + !========================================================== + !-----CLASS 4; FREE CONVECTION: + !========================================================== + + !COMPUTE z/L first guess: + IF (itimestep .LE. 1) THEN + CALL Li_etal_2010(ZOL(I),rb_ice(I),ZA(I)/ZNTstoch_ice(I),zratio_ice(I)) + ELSE + ZOL(I)=ZA(I)*KARMAN*G*MOL(I)/(TH1D(I)*MAX(UST_ice(I)*UST_ice(I),0.001)) + ZOL(I)=MAX(ZOL(I),-50.0) + ZOL(I)=MIN(ZOL(I),0.0) + ENDIF + + IF (debug_code >= 1) THEN + IF (ZNTstoch_ice(i) < 1E-8 .OR. Zt_ice(i) < 1E-10) THEN + write(0,*)"===(ice) capture bad input in mynn sfc layer, i=:",i + write(0,*)"rb=", rb_ice(I)," ZNT=", ZNTstoch_ice(i)," ZT=",Zt_ice(i) + write(0,*)" tsk=", tskin_ice(i)," wstar=",wstar(i)," prev z/L=",ZOL(I),& + " tsurf=", tsurf_ice(i)," qsfc=", qsfc_ice(i)," znt=", znt_ice(i),& + " ust=", ust_ice(i)," snowh=", snowh_ice(i),"psfcpa=",PSFCPA(i), & + " dz=",dz8w1d(i)," qflx=",qflx(i)," hflx=",hflx(i)," hpbl=",pblh(i) + ENDIF + ENDIF + !Use Pedros iterative function to find z/L + zol(I)=zolri(rb_ice(I),ZA(I),ZNTstoch_ice(I),ZT_ice(I),ZOL(I)) + ZOL(I)=MAX(ZOL(I),-50.0) + ZOL(I)=MIN(ZOL(I),0.0) + + zolz0 = zol(I)*ZNTstoch_ice(I)/ZA(I) ! z0/L + zolza = zol(I)*(za(I)+ZNTstoch_ice(I))/za(I) ! (z+z0/L + zol10 = zol(I)*(10.+ZNTstoch_ice(I))/za(I) ! (10+z0)/L + zol2 = zol(I)*(2.+ZNTstoch_ice(I))/za(I) ! (2+z0)/L + + !COMPUTE PSIM and PSIH + !CALL PSI_Hogstrom_1996(PSIM(I),PSIH(I),ZOL(I), ZT_ice(I), ZNTstoch_ice(I), ZA(I)) + !CALL PSI_Businger_1971(PSIM(I),PSIH(I),ZOL(I)) + !CALL PSI_DyerHicks(PSIM(I),PSIH(I),ZOL(I),ZT_ice(I),ZNTstoch_ice(I),ZA(I)) + ! use tables + psim(I)=psim_unstable(zolza)-psim_unstable(zolz0) + psih(I)=psih_unstable(zolza)-psih_unstable(zolz0) + psim10(I)=psim_unstable(zol10)-psim_unstable(zolz0) + psih10(I)=psih_unstable(zol10)-psih_unstable(zolz0) + psih2(I)=psih_unstable(zol2)-psih_unstable(zolz0) + + !---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*GZ1OZt_ice(I)) + PSIM(I)=MIN(PSIM(I),0.9*GZ1OZ0_ice(I)) + PSIH2(I)=MIN(PSIH2(I),0.9*GZ2OZt_ice(I)) + PSIM10(I)=MIN(PSIM10(I),0.9*GZ10OZ0_ice(I)) + PSIH10(I)=MIN(PSIH10(I),0.9*GZ10OZt_ice(I)) + + RMOL(I) = ZOL(I)/ZA(I) - !For computing the diagnostics and fluxes (below), whether the fluxes - !are turned off or on, we need the following: - PSIX(I)=GZ1OZ0(I)-PSIM(I) - PSIX10(I)=GZ10OZ0(I)-PSIM10(I) + ENDIF - PSIT(I) =MAX(GZ1OZt(I)-PSIH(I), 1.0) - PSIT2(I)=MAX(GZ2OZt(I)-PSIH2(I),1.0) - PSIT10=MAX(GZ10OZ0(I)-PSIH10(I), 1.0) - - PSIQ=MAX(LOG((ZA(I)+z_q(i))/z_q(I))-PSIH(I) ,1.0) - PSIQ2=MAX(LOG((2.0+z_q(i))/z_q(I))-PSIH2(I) ,1.0) - PSIQ10=MAX(GZ10OZ0(I)-PSIH10(I),1.0) + ! CALCULATE THE RESISTANCE: + PSIX_ice(I) =MAX(GZ1OZ0_ice(I)-PSIM(I) , 1.0) + PSIX10_ice(I)=MAX(GZ10OZ0_ice(I)-PSIM10(I), 1.0) + PSIT_ice(I) =MAX(GZ1OZt_ice(I)-PSIH(I) , 1.0) + PSIT2_ice(I) =MAX(GZ2OZt_ice(I)-PSIH2(I) , 1.0) + PSIQ_ice(I) =MAX(LOG((ZA(I)+ZQ_ice(i))/ZQ_ice(I))-PSIH(I) ,1.0) + PSIQ2_ice(I) =MAX(LOG((2.0+ZQ_ice(i))/ZQ_ice(I))-PSIH2(I) ,1.0) + + ENDIF ! end ice points + + !------------------------------------------------------------ + !-----COMPUTE THE FRICTIONAL VELOCITY: + !------------------------------------------------------------ + + IF (wet(I)) THEN + ! TO PREVENT OSCILLATIONS AVERAGE WITH OLD VALUE + OLDUST = UST_ocn(I) + UST_ocn(I)=0.5*UST_ocn(I)+0.5*KARMAN*WSPD(I)/PSIX_ocn(I) + !NON-AVERAGED: + !UST_ocn(I)=KARMAN*WSPD(I)/PSIX_ocn(I) + stress_ocn(i)=ust_ocn(i)**2 + + ! Compute u* without vconv for use in HFX calc when isftcflx > 0 + WSPDI(I)=MAX(SQRT(U1D(I)*U1D(I)+V1D(I)*V1D(I)), wmin) + USTM(I)=0.5*USTM(I)+0.5*KARMAN*WSPDI(I)/PSIX_ocn(I) + + ENDIF ! end water points + + IF (dry(I)) THEN + ! TO PREVENT OSCILLATIONS AVERAGE WITH OLD VALUE + OLDUST = UST_lnd(I) + UST_lnd(I)=0.5*UST_lnd(I)+0.5*KARMAN*WSPD(I)/PSIX_lnd(I) + !NON-AVERAGED: + !UST_lnd(I)=KARMAN*WSPD(I)/PSIX_lnd(I) + !From Tilden Meyers: + !IF (rb_lnd(I) .GE 0.0) THEN + ! ust_lnd(i)=WSPD_lnd*0.1/(1.0 + 10.0*rb_lnd(I)) + !ELSE + ! ust_lnd(i)=WSPD_lnd*0.1*(1.0 - 10.0*rb_lnd(I))**onethird + !ENDIF + UST_lnd(I)=MAX(UST_lnd(I),0.005) + stress_lnd(i)=ust_lnd(i)**2 + + !set ustm = ust over land. + USTM(I)=UST_lnd(I) + ENDIF ! end water points + + IF (icy(I)) THEN + ! TO PREVENT OSCILLATIONS AVERAGE WITH OLD VALUE + OLDUST = UST_ice(I) + UST_ice(I)=0.5*UST_ice(I)+0.5*KARMAN*WSPD(I)/PSIX_ice(I) + !NON-AVERAGED: + !UST_ice(I)=KARMAN*WSPD(I)/PSIX_ice(I) + UST_ice(I)=MAX(UST_ice(I),0.005) + stress_ice(i)=ust_ice(i)**2 + + !Set ustm = ust over ice. + USTM(I)=UST_ice(I) + ENDIF ! end ice points + + !---------------------------------------------------- + !----COMPUTE THE TEMPERATURE SCALE (a.k.a. FRICTION TEMPERATURE, T*, or MOL) + !----AND COMPUTE THE MOISTURE SCALE (or q*) + !---------------------------------------------------- + + IF (wet(I)) THEN + DTG=THV1D(I)-THVSK_ocn(I) + OLDTST=MOL(I) + MOL(I)=KARMAN*DTG/PSIT_ocn(I)/PRT + !t_star(I) = -HFX(I)/(UST(I)*CPM(I)*RHO1D(I)) + !t_star(I) = MOL(I) + !---------------------------------------------------- + DQG=(QVSH(i)-qsfc_ocn(i))*1000. !(kg/kg -> g/kg) + qstar(I)=KARMAN*DQG/PSIQ_ocn(I)/PRT + ENDIF + + IF (dry(I)) THEN + DTG=THV1D(I)-THVSK_lnd(I) + OLDTST=MOL(I) + MOL(I)=KARMAN*DTG/PSIT_lnd(I)/PRT + !t_star(I) = -HFX(I)/(UST(I)*CPM(I)*RHO1D(I)) + !t_star(I) = MOL(I) + !---------------------------------------------------- + DQG=(QVSH(i)-qsfc_lnd(i))*1000. !(kg/kg -> g/kg) + qstar(I)=KARMAN*DQG/PSIQ_lnd(I)/PRT + ENDIF + + IF (icy(I)) THEN + DTG=THV1D(I)-THVSK_ice(I) + OLDTST=MOL(I) + MOL(I)=KARMAN*DTG/PSIT_ice(I)/PRT + !t_star(I) = -HFX(I)/(UST(I)*CPM(I)*RHO1D(I)) + !t_star(I) = MOL(I) + !---------------------------------------------------- + DQG=(QVSH(i)-qsfc_ice(i))*1000. !(kg/kg -> g/kg) + qstar(I)=KARMAN*DQG/PSIQ_ice(I)/PRT + ENDIF + + ENDDO ! end i-loop + + IF (debug_code == 2) THEN + DO I=its,ite + IF(wet(i))write(*,*)"==== AT END OF MAIN LOOP, i=",i, "(wet)" + IF(dry(i))write(*,*)"==== AT END OF MAIN LOOP, i=",i, "(land)" + IF(icy(i))write(*,*)"==== AT END OF MAIN LOOP, i=",i, "(ice)" + write(*,*)"z/L:",ZOL(I)," wspd:",wspd(I)," Tstar:",MOL(I) + IF(wet(i))write(*,*)"PSIM:",PSIM(I)," PSIH:",PSIH(I)," W*:",WSTAR(I),& + " DTHV:",THV1D(I)-THVSK_ocn(I) + IF(dry(i))write(*,*)"PSIM:",PSIM(I)," PSIH:",PSIH(I)," W*:",WSTAR(I),& + " DTHV:",THV1D(I)-THVSK_lnd(I) + IF(icy(i))write(*,*)"PSIM:",PSIM(I)," PSIH:",PSIH(I)," W*:",WSTAR(I),& + " DTHV:",THV1D(I)-THVSK_ice(i) + write(*,*)"CPM:",CPM(I)," RHO1D:",RHO1D(I)," q*:",qstar(I)," T*:",MOL(I) + IF(wet(i))write(*,*)"U*:",UST_ocn(I)," Z0:",ZNTstoch_ocn(I)," Zt:",zt_ocn(I) + IF(dry(i))write(*,*)"U*:",UST_lnd(I)," Z0:",ZNTstoch_lnd(I)," Zt:",zt_lnd(I) + IF(icy(i))write(*,*)"U*:",UST_ice(I)," Z0:",ZNTstoch_ice(I)," Zt:",zt_ice(I) + write(*,*)"hfx:",HFX(I)," MAVAIL:",MAVAIL(I)," QVSH(I):",QVSH(I) + write(*,*)"=============================================" + ENDDO ! end i-loop + ENDIF + + !---------------------------------------------------------- + ! COMPUTE SURFACE HEAT AND MOISTURE FLUXES + !---------------------------------------------------------- + DO I=its,ite - IF (ISFFLX .LT. 1) THEN + IF (ISFFLX .LT. 1) THEN QFX(i) = 0. - HFX(i) = 0. + HFX(i) = 0. + HFLX(i) = 0. FLHC(I) = 0. FLQC(I) = 0. LH(I) = 0. CHS(I) = 0. CH(I) = 0. CHS2(i) = 0. - CQS2(i) = 0. - IF(PRESENT(ck) .and. PRESENT(cd) .and. & - &PRESENT(cka) .and. PRESENT(cda)) THEN - Ck(I) = 0. - Cd(I) = 0. - Cka(I)= 0. - Cda(I)= 0. - ENDIF + CQS2(i) = 0. + ch_ocn(I)= 0. + cm_ocn(I)= 0. + ch_lnd(I)= 0. + cm_lnd(I)= 0. + ch_ice(I)= 0. + cm_ice(I)= 0. + ELSE - IF((XLAND(I)-1.5).LT.0)THEN !LAND Only - IF ( IZ0TLND .EQ. 4 ) THEN - CALL Pan_etal_1994(PSIQ,PSIQ2,UST(I),PSIH(I),PSIH2(I),& - & KARMAN,ZA(I)) + IF (dry(i)) THEN + + !------------------------------------------ + ! CALCULATE THE EXCHANGE COEFFICIENTS FOR HEAT (FLHC) + ! AND MOISTURE (FLQC) + !------------------------------------------ + FLQC(I)=RHO1D(I)*MAVAIL(I)*UST_lnd(I)*KARMAN/PSIQ_lnd(i) + FLHC(I)=RHO1D(I)*CPM(I)*UST_lnd(I)*KARMAN/PSIT_lnd(I) + + IF (compute_flux) THEN + !---------------------------------- + ! COMPUTE SURFACE MOISTURE FLUX: + !---------------------------------- + !QFX(I)=FLQC(I)*(QSFCMR_lnd(I)-QV1D(I)) + QFX(I)=FLQC(I)*(QSFC_lnd(I)-QV1D(I)) + QFX(I)=MAX(QFX(I),-0.02) !allows small neg QFX + LH(i)=XLV*QFX(i) + QFLX(i)=QFX(i)/RHO1D(i) + + !---------------------------------- + ! COMPUTE SURFACE HEAT FLUX: + !---------------------------------- + HFX(I)=FLHC(I)*(THSK_lnd(I)-TH1D(I)) + HFX(I)=MAX(HFX(I),-250.) + HFLX(I)=HFX(I)/(RHO1D(I)*cpm(I)) ENDIF - ENDIF - !------------------------------------------ - ! CALCULATE THE EXCHANGE COEFFICIENTS FOR HEAT (FLHC) - ! AND MOISTURE (FLQC) - !------------------------------------------ - FLQC(I)=RHO1D(I)*MAVAIL(I)*UST(I)*KARMAN/PSIQ - FLHC(I)=RHO1D(I)*CPM(I)*UST(I)*KARMAN/PSIT(I) - !OLD WAY: - !DTTHX=ABS(TH1D(I)-THGB(I)) - !IF(DTTHX.GT.1.E-5)THEN - ! FLHC(I)=CPM(I)*RHO1D(I)*UST(I)*MOL(I)/(TH1D(I)-THGB(I)) - !ELSE - ! FLHC(I)=0. - !ENDIF - - !---------------------------------- - ! COMPUTE SURFACE MOISTURE FLUX: - !---------------------------------- - QFX(I)=FLQC(I)*(QSFCMR(I)-QV1D(I)) - !JOE: QFX(I)=MAX(QFX(I),0.) !originally did not allow neg QFX - QFX(I)=MAX(QFX(I),-0.02) !allows small neg QFX, like MYJ - LH(I)=XLV*QFX(I) - - !---------------------------------- - ! COMPUTE SURFACE HEAT FLUX: - !---------------------------------- - IF(XLAND(I)-1.5.GT.0.)THEN !WATER - HFX(I)=FLHC(I)*(THGB(I)-TH1D(I)) - IF ( PRESENT(ISFTCFLX) ) THEN - IF ( ISFTCFLX.NE.0 ) THEN - ! AHW: add dissipative heating term - HFX(I)=HFX(I)+RHO1D(I)*USTM(I)*USTM(I)*WSPDI(I) + !TRANSFER COEFF FOR SOME LSMs: + !CHS(I)=UST(I)*KARMAN/(ALOG(KARMAN*UST(I)*ZA(I) & + ! /XKA+ZA(I)/ZL)-PSIH(I)) + CHS(I)=UST_lnd(I)*KARMAN/PSIT_lnd(I) + + !THESE ARE USED FOR 2-M DIAGNOSTICS ONLY + CQS2(I)=UST_lnd(I)*KARMAN/PSIQ2_lnd(i) + CHS2(I)=UST_lnd(I)*KARMAN/PSIT2_lnd(I) + + QSFC(I)=QSFC_lnd(I) + + ELSEIF (wet(i)) THEN + + !------------------------------------------ + ! CALCULATE THE EXCHANGE COEFFICIENTS FOR HEAT (FLHC) + ! AND MOISTURE (FLQC) + !------------------------------------------ + FLQC(I)=RHO1D(I)*MAVAIL(I)*UST_ocn(I)*KARMAN/PSIQ_ocn(i) + FLHC(I)=RHO1D(I)*CPM(I)*UST_ocn(I)*KARMAN/PSIT_ocn(I) + + IF (compute_flux) THEN + !---------------------------------- + ! COMPUTE SURFACE MOISTURE FLUX: + !---------------------------------- + !QFX(I)=FLQC(I)*(QSFCMR_ocn(I)-QV1D(I)) + QFX(I)=FLQC(I)*(QSFC_ocn(I)-QV1D(I)) + QFX(I)=MAX(QFX(I),-0.02) !allows small neg QFX + LH(I)=XLV*QFX(I) + QFLX(i)=QFX(i)/RHO1D(i) + + !---------------------------------- + ! COMPUTE SURFACE HEAT FLUX: + !---------------------------------- + HFX(I)=FLHC(I)*(THSK_ocn(I)-TH1D(I)) + IF ( PRESENT(ISFTCFLX) ) THEN + IF ( ISFTCFLX.NE.0 ) THEN + ! AHW: add dissipative heating term + HFX(I)=HFX(I)+RHO1D(I)*USTM(I)*USTM(I)*WSPDI(I) + ENDIF ENDIF + HFLX(I)=HFX(I)/(RHO1D(I)*cpm(I)) ENDIF - ELSEIF(XLAND(I)-1.5.LT.0.)THEN !LAND - HFX(I)=FLHC(I)*(THGB(I)-TH1D(I)) - HFX(I)=MAX(HFX(I),-250.) - ENDIF - !CHS(I)=UST(I)*KARMAN/(ALOG(KARMAN*UST(I)*ZA(I) & - ! /XKA+ZA(I)/ZL)-PSIH(I)) + !TRANSFER COEFF FOR SOME LSMs: + !CHS(I)=UST(I)*KARMAN/(ALOG(KARMAN*UST(I)*ZA(I) & + ! /XKA+ZA(I)/ZL)-PSIH(I)) + CHS(I)=UST_ocn(I)*KARMAN/PSIT_ocn(I) + + !THESE ARE USED FOR 2-M DIAGNOSTICS ONLY + CQS2(I)=UST_ocn(I)*KARMAN/PSIQ2_ocn(i) + CHS2(I)=UST_ocn(I)*KARMAN/PSIT2_ocn(I) + + QSFC(I)=QSFC_ocn(I) + + ELSEIF (icy(i)) THEN + + !------------------------------------------ + ! CALCULATE THE EXCHANGE COEFFICIENTS FOR HEAT (FLHC) + ! AND MOISTURE (FLQC) + !------------------------------------------ + FLQC(I)=RHO1D(I)*MAVAIL(I)*UST_ice(I)*KARMAN/PSIQ_ice(i) + FLHC(I)=RHO1D(I)*CPM(I)*UST_ice(I)*KARMAN/PSIT_ice(I) + + IF (compute_flux) THEN + !---------------------------------- + ! COMPUTE SURFACE MOISTURE FLUX: + !---------------------------------- + !QFX(I)=FLQC(I)*(QSFCMR_ice(I)-QV1D(I)) + QFX(I)=FLQC(I)*(QSFC_ice(I)-QV1D(I)) + QFX(I)=MAX(QFX(I),-0.02) !allows small neg QFX + LH(I)=XLF*QFX(I) + QFLX(i)=QFX(i)/RHO1D(i) + + !---------------------------------- + ! COMPUTE SURFACE HEAT FLUX: + !---------------------------------- + HFX(I)=FLHC(I)*(THSK_ice(I)-TH1D(I)) + HFX(I)=MAX(HFX(I),-250.) + HFLX(I)=HFX(I)/(RHO1D(I)*cpm(I)) + ENDIF - CHS(I)=UST(I)*KARMAN/PSIT(I) + !TRANSFER COEFF FOR SOME LSMs: + !CHS(I)=UST(I)*KARMAN/(ALOG(KARMAN*UST(I)*ZA(I) & + ! /XKA+ZA(I)/ZL)-PSIH(I)) + CHS(I)=UST_ice(I)*KARMAN/PSIT_ice(I) - ! The exchange coefficient for cloud water is assumed to be the - ! same as that for heat. CH is multiplied by WSPD. + !THESE ARE USED FOR 2-M DIAGNOSTICS ONLY + CQS2(I)=UST_ice(I)*KARMAN/PSIQ2_ice(i) + CHS2(I)=UST_ice(I)*KARMAN/PSIT2_ice(I) - !ch(i)=chs(i) - ch(i)=flhc(i)/( cpm(i)*RHO1D(i) ) + QSFC(I)=QSFC_ice(I) - !THESE ARE USED FOR 2-M DIAGNOSTICS ONLY - CQS2(I)=UST(I)*KARMAN/PSIQ2 - CHS2(I)=UST(I)*KARMAN/PSIT2(I) + ENDIF - IF(PRESENT(ck) .and. PRESENT(cd) .and. & - &PRESENT(cka) .and. PRESENT(cda)) THEN - Ck(I)=(karman/psix10(I))*(karman/psiq10) - Cd(I)=(karman/psix10(I))*(karman/psix10(I)) - Cka(I)=(karman/psix(I))*(karman/psiq) - Cda(I)=(karman/psix(I))*(karman/psix(I)) + IF (debug_code >= 1) THEN + write(*,*)"QFX=",QFX(I),"FLQC=",FLQC(I) + if(icy(i))write(*,*)"ice, MAVAIL:",MAVAIL(I)," u*=",UST_ice(I)," psiq=",PSIQ_ice(i) + if(dry(i))write(*,*)"lnd, MAVAIL:",MAVAIL(I)," u*=",UST_lnd(I)," psiq=",PSIQ_lnd(i) + if(wet(i))write(*,*)"ocn, MAVAIL:",MAVAIL(I)," u*=",UST_ocn(I)," psiq=",PSIQ_ocn(i) ENDIF - ENDIF !end ISFFLX option + ! The exchange coefficient for cloud water is assumed to be the + ! same as that for heat. CH is multiplied by WSPD. + ch(i)=flhc(i)/( cpm(i)*RHO1D(i) ) + + !----------------------------------------- + !--- COMPUTE EXCHANGE COEFFICIENTS FOR FV3 + !----------------------------------------- + IF (wet(i)) THEN + ch_ocn(I)=(karman/psix_ocn(I))*(karman/psit_ocn(i)) + cm_ocn(I)=(karman/psix_ocn(I))*(karman/psix_ocn(I)) + ENDIF + IF (dry(i)) THEN + ch_lnd(I)=(karman/psix_lnd(I))*(karman/psit_lnd(i)) + cm_lnd(I)=(karman/psix_lnd(I))*(karman/psix_lnd(I)) + ENDIF + IF (icy(i)) THEN + ch_ice(I)=(karman/psix_ice(I))*(karman/psit_ice(i)) + cm_ice(I)=(karman/psix_ice(I))*(karman/psix_ice(I)) + ENDIF - !----------------------------------------------------- - !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) .le. 7.0) then - ! high vertical resolution - if(ZA2(i) .gt. 7.0 .and. ZA2(i) .lt. 13.0) then - !use 2nd model level - U10(I)=U1D2(I) - V10(I)=V1D2(I) + ENDIF !end ISFFLX option +ENDDO ! end i-loop + +IF (compute_diag) then + DO I=its,ite + !----------------------------------------------------- + !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) .le. 7.0) then + ! high vertical resolution + if(ZA2(i) .gt. 7.0 .and. ZA2(i) .lt. 13.0) then + !use 2nd model level + U10(I)=U1D2(I) + V10(I)=V1D2(I) + else + IF (dry(i)) THEN + !U10(I)=U1D(I)*PSIX10_lnd(I)/PSIX_lnd(I) + !V10(I)=V1D(I)*PSIX10_lnd(I)/PSIX_lnd(I) + !use neutral-log: + U10(I)=U1D(I)*log(10./ZNTstoch_lnd(I))/log(ZA(I)/ZNTstoch_lnd(I)) + V10(I)=V1D(I)*log(10./ZNTstoch_lnd(I))/log(ZA(I)/ZNTstoch_lnd(I)) + ELSEIF (wet(i)) THEN + U10(I)=U1D(I)*log(10./ZNTstoch_ocn(I))/log(ZA(I)/ZNTstoch_ocn(I)) + V10(I)=V1D(I)*log(10./ZNTstoch_ocn(I))/log(ZA(I)/ZNTstoch_ocn(I)) + ELSEIF (icy(i)) THEN + U10(I)=U1D(I)*log(10./ZNTstoch_ice(I))/log(ZA(I)/ZNTstoch_ice(I)) + V10(I)=V1D(I)*log(10./ZNTstoch_ice(I))/log(ZA(I)/ZNTstoch_ice(I)) + ENDIF + endif + elseif (ZA(i) .gt. 7.0 .and. ZA(i) .lt. 13.0) then + !moderate vertical resolution + IF (dry(i)) THEN + !U10(I)=U1D(I)*PSIX10_lnd(I)/PSIX_lnd(I) + !V10(I)=V1D(I)*PSIX10_lnd(I)/PSIX_lnd(I) + !use neutral-log: + U10(I)=U1D(I)*log(10./ZNTstoch_lnd(I))/log(ZA(I)/ZNTstoch_lnd(I)) + V10(I)=V1D(I)*log(10./ZNTstoch_lnd(I))/log(ZA(I)/ZNTstoch_lnd(I)) + ELSEIF (wet(i)) THEN + U10(I)=U1D(I)*log(10./ZNTstoch_ocn(I))/log(ZA(I)/ZNTstoch_ocn(I)) + V10(I)=V1D(I)*log(10./ZNTstoch_ocn(I))/log(ZA(I)/ZNTstoch_ocn(I)) + ELSEIF (icy(i)) THEN + U10(I)=U1D(I)*log(10./ZNTstoch_ice(I))/log(ZA(I)/ZNTstoch_ice(I)) + V10(I)=V1D(I)*log(10./ZNTstoch_ice(I))/log(ZA(I)/ZNTstoch_ice(I)) + ENDIF else - U10(I)=U1D(I)*log(10./ZNTstoch(I))/log(ZA(I)/ZNTstoch(I)) - V10(I)=V1D(I)*log(10./ZNTstoch(I))/log(ZA(I)/ZNTstoch(I)) + ! very coarse vertical resolution + IF (dry(i)) THEN + U10(I)=U1D(I)*PSIX10_lnd(I)/PSIX_lnd(I) + V10(I)=V1D(I)*PSIX10_lnd(I)/PSIX_lnd(I) + ELSEIF (wet(i)) THEN + U10(I)=U1D(I)*PSIX10_ocn(I)/PSIX_ocn(I) + V10(I)=V1D(I)*PSIX10_ocn(I)/PSIX_ocn(I) + ELSEIF (icy(i)) THEN + U10(I)=U1D(I)*PSIX10_ice(I)/PSIX_ice(I) + V10(I)=V1D(I)*PSIX10_ice(I)/PSIX_ice(I) + ENDIF endif - elseif(ZA(i) .gt. 7.0 .and. ZA(i) .lt. 13.0) then - !moderate vertical resolution - !U10(I)=U1D(I)*PSIX10(I)/PSIX(I) - !V10(I)=V1D(I)*PSIX10(I)/PSIX(I) - !use neutral-log: - U10(I)=U1D(I)*log(10./ZNTstoch(I))/log(ZA(I)/ZNTstoch(I)) - V10(I)=V1D(I)*log(10./ZNTstoch(I))/log(ZA(I)/ZNTstoch(I)) - else - ! very coarse vertical resolution - U10(I)=U1D(I)*PSIX10(I)/PSIX(I) - V10(I)=V1D(I)*PSIX10(I)/PSIX(I) - endif - - !----------------------------------------------------- - !COMPUTE 2m T, TH, AND Q - !THESE WILL BE OVERWRITTEN FOR LAND POINTS IN THE LSM - !----------------------------------------------------- - DTG=TH1D(I)-THGB(I) - TH2(I)=THGB(I)+DTG*PSIT2(I)/PSIT(I) - !*** BE CERTAIN THAT THE 2-M THETA IS BRACKETED BY - !*** THE VALUES AT THE SURFACE AND LOWEST MODEL LEVEL. - IF ((TH1D(I)>THGB(I) .AND. (TH2(I)TH1D(I))) .OR. & - (TH1D(I)THGB(I) .OR. TH2(I)QSFCMR(I) .AND. (Q2(I)QV1D(I))) .OR. & - (QV1D(I)QSFCMR(I) .OR. Q2(I)THSK_lnd(I) .AND. (TH2(I)TH1D(I))) .OR. & + (TH1D(I)THSK_lnd(I) .OR. TH2(I)THSK_ocn(I) .AND. (TH2(I)TH1D(I))) .OR. & + (TH1D(I)THSK_ocn(I) .OR. TH2(I)THSK_ice(I) .AND. (TH2(I)TH1D(I))) .OR. & + (TH1D(I)THSK_ice(I) .OR. TH2(I) 1200. .OR. HFX(I) < -700.)THEN + IF (compute_flux) THEN + IF (HFX(I) > 1200. .OR. HFX(I) < -700.)THEN print*,"SUSPICIOUS VALUES IN MYNN SFCLAYER",& - ITER-ITMAX," ITERATIONS",I,J, "HFX: ",HFX(I) + I,J, "HFX: ",HFX(I) yesno = 1 + ENDIF + IF (LH(I) > 1200. .OR. LH(I) < -700.)THEN + print*,"SUSPICIOUS VALUES IN MYNN SFCLAYER",& + I,J, "LH: ",LH(I) + yesno = 1 + ENDIF ENDIF - IF (LH(I) > 1200. .OR. LH(I) < -700.)THEN + IF (wet(i)) THEN + IF (UST_ocn(I) < 0.0 .OR. UST_ocn(I) > 4.0 )THEN print*,"SUSPICIOUS VALUES IN MYNN SFCLAYER",& - ITER-ITMAX," ITERATIONS",I,J, "LH: ",LH(I) + I,J, "UST_ocn: ",UST_ocn(I) yesno = 1 + ENDIF ENDIF - IF (UST(I) < 0.0 .OR. UST(I) > 4.0 )THEN - print*,"SUSPICIOUS VALUES IN MYNN SFCLAYER",& - ITER-ITMAX," ITERATIONS",I,J, "UST: ",UST(I) + IF (dry(i)) THEN + IF (UST_lnd(I) < 0.0 .OR. UST_lnd(I) > 4.0 )THEN + print*,"SUSPICIOUS VALUES IN MYNN SFCLAYER",& + I,J, "UST_lnd: ",UST_lnd(I) yesno = 1 + ENDIF + ENDIF + IF (icy(i)) THEN + IF (UST_ice(I) < 0.0 .OR. UST_ice(I) > 4.0 )THEN + print*,"SUSPICIOUS VALUES IN MYNN SFCLAYER",& + I,J, "UST_ice: ",UST_ice(I) + yesno = 1 + ENDIF ENDIF IF (WSTAR(I)<0.0 .OR. WSTAR(I) > 6.0)THEN print*,"SUSPICIOUS VALUES IN MYNN SFCLAYER",& - ITER-ITMAX," ITERATIONS",I,J, "WSTAR: ",WSTAR(I) + I,J, "WSTAR: ",WSTAR(I) yesno = 1 ENDIF IF (RHO1D(I)<0.0 .OR. RHO1D(I) > 1.6 )THEN print*,"SUSPICIOUS VALUES IN MYNN SFCLAYER",& - ITER-ITMAX," ITERATIONS",I,J, "rho: ",RHO1D(I) + I,J, "rho: ",RHO1D(I) yesno = 1 ENDIF - IF (QSFC(I)*1000. <0.0 .OR. QSFC(I)*1000. >40.)THEN + IF (dry(i)) THEN + IF (QSFC_lnd(I)*1000. <0.0 .OR. QSFC_lnd(I)*1000. >40.)THEN print*,"SUSPICIOUS VALUES IN MYNN SFCLAYER",& - ITER-ITMAX," ITERATIONS",I,J, "QSFC: ",QSFC(I) + I,J, "QSFC_lnd: ",QSFC_lnd(I) yesno = 1 + ENDIF ENDIF IF (PBLH(I)<0. .OR. PBLH(I)>6000.)THEN - print*,"SUSPICIOUS VALUES IN MYNN SFCLAYER",& - ITER-ITMAX," ITERATIONS",I,J, "PBLH: ",PBLH(I) + print*,"SUSPICIOUS VALUES IN MYNN SFCLAYER",& + I,J, "PBLH: ",PBLH(I) yesno = 1 ENDIF IF (yesno == 1) THEN - print*," OTHER INFO:" - write(*,1001)"REGIME:",REGIME(I)," z/L:",ZOL(I)," U*:",UST(I),& + IF (wet(i)) THEN + print*," OTHER INFO over water:" + print*,"z/L:",ZOL(I)," U*:",UST_ocn(I)," Tstar:",MOL(I) + print*,"PSIM:",PSIM(I)," PSIH:",PSIH(I)," W*:",WSTAR(I),& + " DTHV:",THV1D(I)-THVSK_ocn(I) + print*,"CPM:",CPM(I)," RHO1D:",RHO1D(I)," L:",& + ZOL(I)/ZA(I)," DTH:",TH1D(I)-THSK_ocn(I) + print*," Z0:",ZNTstoch_ocn(I)," Zt:",ZT_ocn(I)," za:",za(I) + print*,"MAVAIL:",MAVAIL(I)," QSFC_ocn(I):",& + QSFC_ocn(I)," QVSH(I):",QVSH(I) + print*,"PSIX=",PSIX_ocn(I)," T1D(i):",T1D(i) + write(*,*)"=============================================" + ENDIF + IF (dry(i)) THEN + print*," OTHER INFO over land:" + print*,"z/L:",ZOL(I)," U*:",UST_lnd(I),& + " Tstar:",MOL(I) + print*,"PSIM:",PSIM(I)," PSIH:",PSIH(I)," W*:",WSTAR(I),& + " DTHV:",THV1D(I)-THVSK_lnd(I) + print*,"CPM:",CPM(I)," RHO1D:",RHO1D(I)," L:",& + ZOL(I)/ZA(I)," DTH:",TH1D(I)-THSK_lnd(I) + print*," Z0:",ZNTstoch_lnd(I)," Zt:",ZT_lnd(I)," za:",za(I) + print*," MAVAIL:",MAVAIL(I)," QSFC_lnd(I):",& + QSFC_lnd(I)," QVSH(I):",QVSH(I) + print*,"PSIX=",PSIX_lnd(I)," T1D(i):",T1D(i) + write(*,*)"=============================================" + ENDIF + IF (icy(i)) THEN + print*," OTHER INFO:" + print*,"z/L:",ZOL(I)," U*:",UST_ice(I),& " Tstar:",MOL(I) - write(*,1002)"PSIM:",PSIM(I)," PSIH:",PSIH(I)," W*:",WSTAR(I),& - " DTHV:",THV1D(I)-THVGB(I) - write(*,1003)"CPM:",CPM(I)," RHO1D:",RHO1D(I)," L:",& - ZOL(I)/ZA(I)," DTH:",TH1D(I)-THGB(I) - write(*,1004)"Z0/Zt:",zratio(I)," Z0:",ZNTstoch(I)," Zt:",z_t(I),& - " za:",za(I) - write(*,1005)"Re:",restar," MAVAIL:",MAVAIL(I)," QSFC(I):",& - QSFC(I)," QVSH(I):",QVSH(I) - print*,"PSIX=",PSIX(I)," Z0:",ZNTstoch(I)," T1D(i):",T1D(i) - write(*,*)"=============================================" + print*,"PSIM:",PSIM(I)," PSIH:",PSIH(I)," W*:",WSTAR(I),& + " DTHV:",THV1D(I)-THVSK_ice(I) + print*,"CPM:",CPM(I)," RHO1D:",RHO1D(I)," L:",& + ZOL(I)/ZA(I)," DTH:",TH1D(I)-THSK_ice(I) + print*," Z0:",ZNTstoch_ice(I)," Zt:",ZT_ice(I)," za:",za(I) + print*," MAVAIL:",MAVAIL(I)," QSFC_ice(I):",& + QSFC_ice(I)," QVSH(I):",QVSH(I) + print*,"PSIX=",PSIX_ice(I)," T1D(i):",T1D(i) + write(*,*)"=============================================" + ENDIF ENDIF - ENDIF - - ENDDO !end i-loop + ENDDO ! end i-loop + ENDIF ! end debug option END SUBROUTINE SFCLAY1D_mynn !------------------------------------------------------------------- @@ -1411,23 +2124,20 @@ SUBROUTINE zilitinkevich_1995(Z_0,Zt,Zq,restar,ustar,KARMAN,& IF ( IZ0TLND2 .EQ. 1 ) THEN CZIL = 10.0 ** ( -0.40 * ( Z_0 / 0.07 ) ) ELSE - CZIL = 0.075 !0.10 + CZIL = 0.085 !0.075 !0.10 END IF Zt = Z_0*EXP(-KARMAN*CZIL*SQRT(restar)) - Zt = MIN( Zt, Z_0/2.) + Zt = MIN( Zt, 0.75*Z_0) Zq = Z_0*EXP(-KARMAN*CZIL*SQRT(restar)) - Zq = MIN( Zq, Z_0/2.) + Zq = MIN( Zq, 0.75*Z_0) -! perturb thermal and moisture roughness lenth by +/-50% -! uses same perturbation pattern for perturbing cloud fraction -! and turbulent mixing length (module_sf_mynn.F), but -! twice the amplitude; -! multiplication with -1.0 anticorrelates patterns +! stochastically perturb thermal and moisture roughness length. +! currently set to half the amplitude: if (spp_pbl==1) then - Zt = Zt + Zt * 2.0 * rstoch - Zt = MAX(Zt, 0.001) + Zt = Zt + Zt * 0.5 * rstoch + Zt = MAX(Zt, 0.0001) Zq = Zt endif @@ -1437,60 +2147,26 @@ SUBROUTINE zilitinkevich_1995(Z_0,Zt,Zq,restar,ustar,KARMAN,& END SUBROUTINE zilitinkevich_1995 !-------------------------------------------------------------------- -!>\ingroup module_sf_mynn_mod -!! This subroutine returns the resistance (PSIQ) for moisture -!! exchange. This is a modified form originating from Pan et al.. -!! (1994) but modified according to tests in both the RUC model. -!! and WRF-ARW. Note that it is very similar to Carlson and -!! Boland (1978) model (include below in comments) but has an -!! extra molecular layer (a third layer) instead of two layers. - SUBROUTINE Pan_etal_1994(PSIQ,PSIQ2,ustar,psih,psih2,KARMAN,Z1) - - IMPLICIT NONE - REAL, INTENT(IN) :: Z1,ustar,KARMAN,psih,psih2 - REAL, INTENT(OUT) :: psiq,psiq2 - REAL, PARAMETER :: Cpan=1.0 !was 20.8 in Pan et al 1994 - REAL, PARAMETER :: ZL=0.01 - REAL, PARAMETER :: ZMUs=0.2E-3 - REAL, PARAMETER :: XKA = 2.4E-5 - - !PAN et al. (1994): 3-layer model, as in paper: - !ZMU = Cpan*XKA/(KARMAN*UST(I)) - !PSIQ =MAX(KARMAN*ustar*ZMU/XKA + LOG((KARMAN*ustar*ZL + XKA)/XKA + & - ! & Z1/ZL) - PSIH,2.0) - !PSIQ2=MAX(KARMAN*ustar*ZMU/XKA + LOG((KARMAN*ustar*ZL + XKA)/XKA + & - ! & 2./ZL) - PSIH2,2.0) - !MODIFIED FORM: - PSIQ =MAX(KARMAN*ustar*ZMUs/XKA + LOG((KARMAN*ustar*Z1)/XKA + & - & Z1/ZL) - PSIH,2.0) - PSIQ2=MAX(KARMAN*ustar*ZMUs/XKA + LOG((KARMAN*ustar*2.0)/XKA + & - & 2./ZL) - PSIH2,2.0) - - !CARLSON AND BOLAND (1978): 2-layer model - !PSIQ =MAX(LOG(KARMAN*ustar*Z1/XKA + Z1/ZL)-PSIH ,2.0) - !PSIQ2=MAX(LOG(KARMAN*ustar*2./XKA + 2./ZL)-PSIH2 ,2.0) - - END SUBROUTINE Pan_etal_1994 -!-------------------------------------------------------------- -!>\ingroup module_sf_mynn_mod -!! This formulation for roughness length was designed to match. -!! the labratory experiments of Donelan et al. (2004). -!! This is an update version from Davis et al. 2008, which -!! corrects a small-bias in Z_0 (AHW real-time 2012). SUBROUTINE davis_etal_2008(Z_0,ustar) + !a.k.a. : Donelan et al. (2004) + !This formulation for roughness length was designed to match + !the labratory experiments of Donelan et al. (2004). + !This is an update version from Davis et al. 2008, which + !corrects a small-bias in Z_0 (AHW real-time 2012). + IMPLICIT NONE REAL, INTENT(IN) :: ustar REAL, INTENT(OUT) :: Z_0 REAL :: ZW, ZN1, ZN2 REAL, PARAMETER :: G=9.81, OZO=1.59E-5 - !OLD FORM: Z_0 = 10.*EXP(-10./(ustar**(1./3.))) + !OLD FORM: Z_0 = 10.*EXP(-10./(ustar**onethird)) !NEW FORM: ZW = MIN((ustar/1.06)**(0.3),1.0) ZN1 = 0.011*ustar*ustar/G + OZO - ZN2 = 10.*exp(-9.5*ustar**(-.3333)) + & + ZN2 = 10.*exp(-9.5*ustar**(-onethird)) + & 0.11*1.5E-5/AMAX1(ustar,0.01) Z_0 = (1.0-ZW) * ZN1 + ZW * ZN2 @@ -1623,11 +2299,12 @@ END SUBROUTINE garratt_1992 !!(1992, p. 102), is available for flows with Ren < 2. !! !!This is for use over water only. - SUBROUTINE fairall_etal_2003(Zt,Zq,Ren,ustar,visc) + SUBROUTINE fairall_etal_2003(Zt,Zq,Ren,ustar,visc,rstoch,spp_pbl) IMPLICIT NONE - REAL, INTENT(IN) :: Ren,ustar,visc - REAL, INTENT(OUT) :: Zt,Zq + REAL, INTENT(IN) :: Ren,ustar,visc,rstoch + INTEGER, INTENT(IN):: spp_pbl + REAL, INTENT(OUT) :: Zt,Zq IF (Ren .le. 2.) then @@ -1645,6 +2322,11 @@ SUBROUTINE fairall_etal_2003(Zt,Zq,Ren,ustar,visc) ENDIF + if (spp_pbl==1) then + Zt = Zt + Zt * 0.5 * rstoch + Zq = Zt + endif + Zt = MIN(Zt,1.0e-4) Zt = MAX(Zt,2.0e-9) @@ -1673,8 +2355,8 @@ SUBROUTINE fairall_etal_2014(Zt,Zq,Ren,ustar,visc,rstoch,spp_pbl) Zq = Zt IF (spp_pbl ==1) THEN - Zt = MAX(Zt + Zt*2.0*rstoch,2.0e-9) - Zq = MAX(Zt + Zt*2.0*rstoch,2.0e-9) + Zt = MAX(Zt + Zt*0.5*rstoch,2.0e-9) + Zq = MAX(Zt + Zt*0.5*rstoch,2.0e-9) ELSE Zt = MAX(Zt,2.0e-9) Zq = MAX(Zt,2.0e-9) @@ -1708,10 +2390,10 @@ END SUBROUTINE fairall_etal_2014 !!Zt was reduced too much for low-moderate positive heat fluxes. !! !!This should only be used over land! - SUBROUTINE Yang_2008(Z_0,Zt,Zq,ustar,tstar,qst,Ren,visc,landsea) + SUBROUTINE Yang_2008(Z_0,Zt,Zq,ustar,tstar,qst,Ren,visc) IMPLICIT NONE - REAL, INTENT(IN) :: Z_0, Ren, ustar, tstar, qst, visc, landsea + REAL, INTENT(IN) :: Z_0, Ren, ustar, tstar, qst, visc REAL :: ht, &! roughness height at critical Reynolds number tstar2, &! bounded T*, forced to be non-positive qstar2, &! bounded q*, forced to be non-positive @@ -1994,12 +2676,31 @@ SUBROUTINE PSI_Suselj_Sood_2010(psi_m, psi_h, zL) END SUBROUTINE PSI_Suselj_Sood_2010 !-------------------------------------------------------------------- -!>\ingroup module_sf_mynn_mod -!>This subroutine returns a more robust z/L that best matches -!! the z/L from Hogstrom (1996) for unstable conditions and Beljaars -!! and Holtslag (1991) for stable conditions. + SUBROUTINE PSI_CB2005(psim1,psih1,zL,z0L) + + ! This subroutine returns the stability functions based off + ! of Cheng and Brutseart (2005, BLM), for use in stable conditions only. + ! The returned values are the combination of psi((za+zo)/L) - psi(z0/L) + + IMPLICIT NONE + REAL, INTENT(IN) :: zL,z0L + REAL, INTENT(OUT) :: psim1,psih1 + + psim1 = -6.1*LOG(zL + (1.+ zL**2.5)**0.4) & + -6.1*LOG(z0L + (1.+ z0L**2.5)**0.4) + psih1 = -5.5*log(zL + (1.+ zL**1.1)**0.90909090909) & + -5.5*log(z0L + (1.+ z0L**1.1)**0.90909090909) + + return + + END SUBROUTINE PSI_CB2005 +!-------------------------------------------------------------------- SUBROUTINE Li_etal_2010(zL, Rib, zaz0, z0zt) + !This subroutine returns a more robust z/L that best matches + !the z/L from Hogstrom (1996) for unstable conditions and Beljaars + !and Holtslag (1991) for stable conditions. + IMPLICIT NONE REAL, INTENT(OUT) :: zL REAL, INTENT(IN) :: Rib, zaz0, z0zt @@ -2054,393 +2755,235 @@ SUBROUTINE Li_etal_2010(zL, Rib, zaz0, z0zt) return END SUBROUTINE Li_etal_2010 - !------------------------------------------------------------------- -!>\ingroup module_sf_mynn_mod -!! This subroutine adds pbl modules so they can be optimized in pbl code - SUBROUTINE mym_condensation (kts,kte, & - & dx, dz, & - & thl, qw, & - & p,exner, & - & tsq, qsq, cov, & - & Sh, el, bl_mynn_cloudpdf,& - & qc_bl1D, cldfra_bl1D, & - & PBLH1,HFX1, & - & Vt, Vq, th, sgm) + REAL function zolri(ri,za,z0,zt,zol1) -!------------------------------------------------------------------- + ! This iterative algorithm was taken from the revised surface layer + ! scheme in WRF-ARW, written by Pedro Jimenez and Jimy Dudhia and + ! summarized in Jimenez et al. (2012, MWR). This function was adapted + ! to input the thermal roughness length, zt, (as well as z0) because + ! zt is necessary input for the Dyer-Hicks functions used in MYNN. - INTEGER, INTENT(IN) :: kts,kte, bl_mynn_cloudpdf - REAL, INTENT(IN) :: dx,PBLH1,HFX1 - REAL, DIMENSION(kts:kte), INTENT(IN) :: dz - REAL, DIMENSION(kts:kte), INTENT(IN) :: p,exner, thl, qw, & - &tsq, qsq, cov, th - - REAL, DIMENSION(kts:kte), INTENT(INOUT) :: vt,vq,sgm - - REAL, DIMENSION(kts:kte) :: qmq,alp,a,bet,b,ql,q1,cld,RH - REAL, DIMENSION(kts:kte), INTENT(OUT) :: qc_bl1D,cldfra_bl1D - DOUBLE PRECISION :: t3sq, r3sq, c3sq - - REAL :: qsl,esat,qsat,tlk,qsat_tl,dqsl,cld0,q1k,eq1,qll,& - &q2p,pt,rac,qt,t,xl,rsl,cpm,cdhdz,Fng,qww,alpha,beta,bb,ls_min,ls,wt - INTEGER :: i,j,k - - REAL :: erf - - !JOE: NEW VARIABLES FOR ALTERNATE SIGMA - REAL::dth,dtl,dqw,dzk - REAL, DIMENSION(kts:kte), INTENT(IN) :: Sh,el - - !JOE: variables for BL clouds - REAL::zagl,cld9,damp,edown,RHcrit,RHmean,RHsum,RHnum,Hshcu,PBLH2,ql_limit - REAL, PARAMETER :: Hfac = 3.0 !cloud depth factor for HFX (m^3/W) - REAL, PARAMETER :: HFXmin = 50.0 !min W/m^2 for BL clouds - REAL :: RH_00L, RH_00O, phi_dz, lfac - REAL, PARAMETER :: cdz = 2.0 - REAL, PARAMETER :: mdz = 1.5 - - !JAYMES: variables for tropopause-height estimation - REAL :: theta1, theta2, ht1, ht2 - INTEGER :: k_tropo - - REAL, PARAMETER :: rr2=0.7071068, rrp=0.3989423 - - k_tropo=5 - - zagl = 0. - - SELECT CASE(bl_mynn_cloudpdf) - - CASE (0) ! ORIGINAL MYNN PARTIAL-CONDENSATION SCHEME - - DO k = kts,kte-1 - t = th(k)*exner(k) - - !SATURATED VAPOR PRESSURE - esat = esat_blend(t) - !SATURATED SPECIFIC HUMIDITY - qsl=ep_2*esat/(p(k)-ep_3*esat) - !dqw/dT: Clausius-Clapeyron - dqsl = qsl*ep_2*ev/( rd*t**2 ) - !RH (0 to 1.0) - RH(k)=MAX(MIN(1.0,qw(k)/MAX(1.E-8,qsl)),0.001) - - alp(k) = 1.0/( 1.0+dqsl*xlvcp ) - bet(k) = dqsl*exner(k) - - !NOTE: negative bl_mynn_cloudpdf will zero-out the stratus subgrid clouds - ! at the end of this subroutine. - !Sommeria and Deardorff (1977) scheme, as implemented - !in Nakanishi and Niino (2009), Appendix B - t3sq = MAX( tsq(k), 0.0 ) - r3sq = MAX( qsq(k), 0.0 ) - c3sq = cov(k) - c3sq = SIGN( MIN( ABS(c3sq), SQRT(t3sq*r3sq) ), c3sq ) - r3sq = r3sq +bet(k)**2*t3sq -2.0*bet(k)*c3sq - !DEFICIT/EXCESS WATER CONTENT - qmq(k) = qw(k) -qsl - !ORIGINAL STANDARD DEVIATION: limit e-6 produces ~10% more BL clouds - !than e-10 - sgm(k) = SQRT( MAX( r3sq, 1.0d-10 )) - !NORMALIZED DEPARTURE FROM SATURATION - q1(k) = qmq(k) / sgm(k) - !CLOUD FRACTION. rr2 = 1/SQRT(2) = 0.707 - cld(k) = 0.5*( 1.0+erf( q1(k)*rr2 ) ) - - END DO - - CASE (1, -1) !ALTERNATIVE FORM (Nakanishi & Niino 2004 BLM, eq. B6, and - !Kuwano-Yoshida et al. 2010 QJRMS, eq. 7): - DO k = kts,kte-1 - t = th(k)*exner(k) - !SATURATED VAPOR PRESSURE - esat = esat_blend(t) - !SATURATED SPECIFIC HUMIDITY - qsl=ep_2*esat/(p(k)-ep_3*esat) - !dqw/dT: Clausius-Clapeyron - dqsl = qsl*ep_2*ev/( rd*t**2 ) - !RH (0 to 1.0) - RH(k)=MAX(MIN(1.0,qw(k)/MAX(1.E-8,qsl)),0.001) - - alp(k) = 1.0/( 1.0+dqsl*xlvcp ) - bet(k) = dqsl*exner(k) - - if (k .eq. kts) then - dzk = 0.5*dz(k) - else - dzk = 0.5*( dz(k) + dz(k-1) ) - end if - dth = 0.5*(thl(k+1)+thl(k)) - 0.5*(thl(k)+thl(MAX(k-1,kts))) - dqw = 0.5*(qw(k+1) + qw(k)) - 0.5*(qw(k) + qw(MAX(k-1,kts))) - sgm(k) = SQRT( MAX( (alp(k)**2 * MAX(el(k)**2,0.1) * & - b2 * MAX(Sh(k),0.03))/4. * & - (dqw/dzk - bet(k)*(dth/dzk ))**2 , 1.0e-10) ) - qmq(k) = qw(k) -qsl - q1(k) = qmq(k) / sgm(k) - cld(k) = 0.5*( 1.0+erf( q1(k)*rr2 ) ) - END DO - - CASE (2, -2) - !Diagnostic statistical scheme of Chaboureau and Bechtold (2002), JAS - !JAYMES- this added 27 Apr 2015 - DO k = kts,kte-1 - t = th(k)*exner(k) - !SATURATED VAPOR PRESSURE - esat = esat_blend(t) - !SATURATED SPECIFIC HUMIDITY - qsl=ep_2*esat/(p(k)-ep_3*esat) - !dqw/dT: Clausius-Clapeyron - dqsl = qsl*ep_2*ev/( rd*t**2 ) - !RH (0 to 1.0) - RH(k)=MAX(MIN(1.0,qw(k)/MAX(1.E-8,qsl)),0.001) - - alp(k) = 1.0/( 1.0+dqsl*xlvcp ) - bet(k) = dqsl*exner(k) - - xl = xl_blend(t) ! obtain latent heat - - tlk = thl(k)*(p(k)/p1000mb)**rcp ! recover liquid temp (tl) from thl - qsat_tl = qsat_blend(tlk,p(k)) ! get saturation water vapor mixing ratio - ! at tl and p - - rsl = xl*qsat_tl / (r_v*tlk**2) ! slope of C-C curve at t = tl - ! CB02, Eqn. 4 - - cpm = cp + qw(k)*cpv ! CB02, sec. 2, para. 1 - - a(k) = 1./(1. + xl*rsl/cpm) ! CB02 variable "a" - - qmq(k) = a(k) * (qw(k) - qsat_tl) ! saturation deficit/excess; - ! the numerator of Q1 - - b(k) = a(k)*rsl ! CB02 variable "b" - - dtl = 0.5*(thl(k+1)*(p(k+1)/p1000mb)**rcp + tlk) & - & - 0.5*(tlk + thl(MAX(k-1,kts))*(p(MAX(k-1,kts))/p1000mb)**rcp) - - dqw = 0.5*(qw(k+1) + qw(k)) - 0.5*(qw(k) + qw(MAX(k-1,kts))) - - if (k .eq. kts) then - dzk = 0.5*dz(k) - else - dzk = 0.5*( dz(k) + dz(k-1) ) - end if - - cdhdz = dtl/dzk + (g/cpm)*(1.+qw(k)) ! expression below Eq. 9 - ! in CB02 - zagl = zagl + dz(k) - ls_min = MIN(MAX(zagl,25.),300.) ! Let this be the minimum possible length scale: - ! 25 m < ls_min(=zagl) < 300 m - lfac=MIN(4.25+dx/4000.,6.) ! A dx-dependent multiplier for the master length scale: - ! lfac(750 m) = 4.4 - ! lfac(3 km) = 5.0 - ! lfac(13 km) = 6.0 - - ls = MAX(MIN(lfac*el(k),900.),ls_min) ! Bounded: ls_min < ls < 900 m - ! Note: CB02 use 900 m as a constant free-atmosphere length scale. - ! Above 300 m AGL, ls_min remains 300 m. For dx = 3 km, the - ! MYNN master length scale (el) must exceed 60 m before ls - ! becomes responsive to el, otherwise ls = ls_min = 300 m. - - sgm(k) = MAX(1.e-10, 0.225*ls*SQRT(MAX(0., & ! Eq. 9 in CB02: - & (a(k)*dqw/dzk)**2 & ! < 1st term in brackets, - & -2*a(k)*b(k)*cdhdz*dqw/dzk & ! < 2nd term, - & +b(k)**2 * cdhdz**2))) ! < 3rd term - ! CB02 use a multiplier of 0.2, but 0.225 is chosen - ! based on tests - - q1(k) = qmq(k) / sgm(k) ! Q1, the normalized saturation - - cld(k) = MAX(0., MIN(1., 0.5+0.36*ATAN(1.55*q1(k)))) ! Eq. 7 in CB02 - - END DO - END SELECT - - zagl = 0. - RHsum=0. - RHnum=0. - RHmean=0.1 !initialize with small value for small PBLH cases - damp =0 - PBLH2=MAX(10.,PBLH1) - - SELECT CASE(bl_mynn_cloudpdf) - - CASE (-1 : 1) ! ORIGINAL MYNN PARTIAL-CONDENSATION SCHEME - ! OR KUWANO ET AL. - DO k = kts,kte-1 - t = th(k)*exner(k) - q1k = q1(k) - zagl = zagl + dz(k) - !q1=0. - !cld(k)=0. - - !COMPUTE MEAN RH IN PBL (NOT PRESSURE WEIGHTED). - IF (zagl < PBLH2 .AND. PBLH2 > 400.) THEN - RHsum=RHsum+RH(k) - RHnum=RHnum+1.0 - RHmean=RHsum/RHnum - ENDIF - RHcrit = 1. - 0.35*(1.0 - (MAX(250.- MAX(HFX1,HFXmin),0.0)/200.)**2) - if (HFX1 > HFXmin) then - cld9=MIN(MAX(0., (rh(k)-RHcrit)/(1.1-RHcrit)), 1.)**2 - else - cld9=0.0 - endif + IMPLICIT NONE + REAL, INTENT(IN) :: ri,za,z0,zt,zol1 + REAL :: x1,x2,fx1,fx2 + INTEGER :: n - edown=PBLH2*.1 - !Vary BL cloud depth (Hshcu) by mean RH in PBL and HFX - !(somewhat following results from Zhang and Klein (2013, JAS)) - Hshcu=200. + (RHmean+0.5)**1.5*MAX(HFX1,0.)*Hfac - if (zagl < PBLH2-edown) then - damp=MIN(1.0,exp(-ABS(((PBLH2-edown)-zagl)/edown))) - elseif(zagl >= PBLH2-edown .AND. zagl < PBLH2+Hshcu)then - damp=1. - elseif (zagl >= PBLH2+Hshcu)then - damp=MIN(1.0,exp(-ABS((zagl-(PBLH2+Hshcu))/500.))) - endif - cldfra_bl1D(k)=cld9*damp - !cldfra_bl1D(k)=cld(k) ! JAYMES: use this form to retain the Sommeria-Deardorff value - - !use alternate cloud fraction to estimate qc for use in BL clouds-radiation - eq1 = rrp*EXP( -0.5*q1k*q1k ) - qll = MAX( cldfra_bl1D(k)*q1k + eq1, 0.0 ) - !ESTIMATED LIQUID WATER CONTENT (UNNORMALIZED) - ql (k) = alp(k)*sgm(k)*qll - if(cldfra_bl1D(k)>0.01 .and. ql(k)<1.E-6)ql(k)=1.E-6 - qc_bl1D(k)=ql(k)*damp - !now recompute estimated lwc for PBL scheme's use - !qll IS THE NORMALIZED LIQUID WATER CONTENT (Sommeria and - !Deardorff (1977, eq 29a). rrp = 1/(sqrt(2*pi)) = 0.3989 - eq1 = rrp*EXP( -0.5*q1k*q1k ) - qll = MAX( cld(k)*q1k + eq1, 0.0 ) - !ESTIMATED LIQUID WATER CONTENT (UNNORMALIZED) - ql (k) = alp(k)*sgm(k)*qll - - q2p = xlvcp/exner(k) - pt = thl(k) +q2p*ql(k) ! potential temp - - !qt is a THETA-V CONVERSION FOR TOTAL WATER (i.e., THETA-V = qt*THETA) - qt = 1.0 +p608*qw(k) -(1.+p608)*ql(k) - rac = alp(k)*( cld(k)-qll*eq1 )*( q2p*qt-(1.+p608)*pt ) - - !BUOYANCY FACTORS: wherever vt and vq are used, there is a - !"+1" and "+tv0", respectively, so these are subtracted out here. - !vt is unitless and vq has units of K. - vt(k) = qt-1.0 -rac*bet(k) - vq(k) = p608*pt-tv0 +rac - - !To avoid FPE in radiation driver, when these two quantities are multiplied by eachother, - ! add limit to qc_bl and cldfra_bl: - IF (QC_BL1D(k) < 1E-6 .AND. ABS(CLDFRA_BL1D(k)) > 0.01) QC_BL1D(k)= 1E-6 - IF (CLDFRA_BL1D(k) < 1E-2)THEN - CLDFRA_BL1D(k)=0. - QC_BL1D(k)=0. - ENDIF + if (ri.lt.0.)then + x1=zol1 - 0.02 !-5. + x2=0. + else + x1=0. + x2=zol1 + 0.02 !5. + endif - END DO - CASE ( 2, -2) - ! JAYMES- this option added 8 May 2015 - ! The cloud water formulations are taken from CB02, Eq. 8. - ! "fng" represents the non-Gaussian contribution to the liquid - ! water flux; these formulations are from Cuijpers and Bechtold - ! (1995), Eq. 7. CB95 also draws from Bechtold et al. 1995, - ! hereafter BCMT95 - DO k = kts,kte-1 - t = th(k)*exner(k) - q1k = q1(k) - zagl = zagl + dz(k) - IF (q1k < 0.) THEN - ql (k) = sgm(k)*EXP(1.2*q1k-1) - ELSE IF (q1k > 2.) THEN - ql (k) = sgm(k)*q1k - ELSE - ql (k) = sgm(k)*(EXP(-1.) + 0.66*q1k + 0.086*q1k**2) - ENDIF + n=0 + fx1=zolri2(x1,ri,za,z0,zt) + fx2=zolri2(x2,ri,za,z0,zt) + Do While (abs(x1 - x2) > 0.01 .and. n < 5) + if(abs(fx2).lt.abs(fx1))then + x1=x1-fx1/(fx2-fx1)*(x2-x1) + fx1=zolri2(x1,ri,za,z0,zt) + zolri=x1 + else + x2=x2-fx2/(fx2-fx1)*(x2-x1) + fx2=zolri2(x2,ri,za,z0,zt) + zolri=x2 + endif + n=n+1 + !print*," n=",n," x1=",x1," x2=",x2 + enddo + + if (n==5 .and. abs(x1 - x2) >= 0.01) then + !print*,"iter FAIL, n=",n," Ri=",ri," z/L=",zolri + !Tests results: fails convergence ~ 0.07 % of the time + !set approximate values: + if (ri.lt.0.)then + zolri=ri*5. + else + zolri=ri*8. + endif + !else + ! print*,"iter OK, n=",n," Ri=",ri," z/L=",zolri + endif - !Next, adjust our initial estimates of cldfra and ql based - !on tropopause-height and PBLH considerations - !JAYMES: added 4 Nov 2016 - if ((cld(k) .gt. 0.) .or. (ql(k) .gt. 0.)) then - if (k .le. k_tropo) then - !At and below tropopause: impose an upper limit on ql; assume that - !a maximum of 0.5 percent supersaturation in water vapor can be - !available for cloud production - ql_limit = 0.005 * qsat_blend( th(k)*exner(k), p(k) ) - ql(k) = MIN( ql(k), ql_limit ) - else - !Above tropopause: eliminate subgrid clouds from CB scheme - cld(k) = 0. - ql(k) = 0. - endif - endif + return + end function +!------------------------------------------------------------------- + REAL function zolri2(zol2,ri2,za,z0,zt) - !Buoyancy-flux-related calculations follow... - ! "Fng" represents the non-Gaussian transport factor - ! (non-dimensional) from from Bechtold et al. 1995 - ! (hereafter BCMT95), section 3(c). Their suggested - ! forms for Fng (from their Eq. 20) are: - ! For purposes of the buoyancy flux in stratus, we will use Fng = 1 - Fng = 1. - - xl = xl_blend(t) - bb = b(k)*t/th(k) ! bb is "b" in BCMT95. Their "b" differs from - ! "b" in CB02 (i.e., b(k) above) by a factor - ! of T/theta. Strictly, b(k) above is formulated in - ! terms of sat. mixing ratio, but bb in BCMT95 is - ! cast in terms of sat. specific humidity. The - ! conversion is neglected here. - qww = 1.+0.61*qw(k) - alpha = 0.61*th(k) - beta = (th(k)/t)*(xl/cp) - 1.61*th(k) - - vt(k) = qww - cld(k)*beta*bb*Fng - 1. - vq(k) = alpha + cld(k)*beta*a(k)*Fng - tv0 - ! vt and vq correspond to beta-theta and beta-q, respectively, - ! in NN09, Eq. B8. They also correspond to the bracketed - ! expressions in BCMT95, Eq. 15, since (s*ql/sigma^2) = cldfra*Fng - ! The "-1" and "-tv0" terms are included for consistency with - ! the legacy vt and vq formulations (above). - - ! increase the cloud fraction estimate below PBLH+1km - if (zagl .lt. PBLH2+1000.) cld(k) = MIN( 1., 1.8*cld(k) ) - ! return a cloud condensate and cloud fraction for icloud_bl option: - cldfra_bl1D(k) = cld(k) - qc_bl1D(k) = ql(k) - - !To avoid FPE in radiation driver, when these two quantities are multiplied by eachother, - ! add limit to qc_bl and cldfra_bl: - IF (QC_BL1D(k) < 1E-6 .AND. ABS(CLDFRA_BL1D(k)) > 0.01) QC_BL1D(k)= 1E-6 - IF (CLDFRA_BL1D(k) < 1E-2)THEN - CLDFRA_BL1D(k)=0. - QC_BL1D(k)=0. - ENDIF + ! INPUT: ================================= + ! zol2 - estimated z/L + ! ri2 - calculated bulk Richardson number + ! za - 1/2 depth of first model layer + ! z0 - aerodynamic roughness length + ! zt - thermal roughness length + ! OUTPUT: ================================ + ! zolri2 - updated estimate of z/L - END DO + IMPLICIT NONE + REAL, INTENT(IN) :: ri2,za,z0,zt + REAL, INTENT(INOUT) :: zol2 + REAL :: zol20,zol3,psim1,psih1,psix2,psit2 - END SELECT !end cloudPDF option + if(zol2*ri2 .lt. 0.)zol2=0. ! limit zol2 - must be same sign as ri2 - !FOR TESTING PURPOSES ONLY, ISOLATE ON THE MASS-CLOUDS. - IF (bl_mynn_cloudpdf .LT. 0) THEN - DO k = kts,kte-1 - cldfra_bl1D(k) = 0.0 - qc_bl1D(k) = 0.0 - END DO - ENDIF + zol20=zol2*z0/za ! z0/L + zol3=zol2+zol20 ! (z+z0)/L + + if (ri2.lt.0) then + !CALL PSI_DyerHicks(psim1,psih1,zol3,zt,z0,za) + psix2=log((za+z0)/z0)-(psim_unstable(zol3)-psim_unstable(zol20)) + psit2=log((za+zt)/zt)-(psih_unstable(zol3)-psih_unstable(zol20)) + !psix2=log((za+z0)/z0)-psim1 + !psit2=log((za+zt)/zt)-psih1 + else + !CALL PSI_DyerHicks(psim1,psih1,zol2,zt,z0,za) + !CALL PSI_CB2005(psim1,psih1,zol3,zol20) + psix2=log((za+z0)/z0)-(psim_stable(zol3)-psim_stable(zol20)) + psit2=log((za+zt)/zt)-(psih_stable(zol3)-psih_stable(zol20)) + !psix2=log((za+z0)/z0)-psim1 + !psit2=log((za+zt)/zt)-psih1 + endif + + zolri2=zol2*psit2/psix2**2 - ri2 - cld(kte) = cld(kte-1) - ql(kte) = ql(kte-1) - vt(kte) = vt(kte-1) - vq(kte) = vq(kte-1) - qc_bl1D(kte)=0. - cldfra_bl1D(kte)=0. + return + end function +!==================================================================== + SUBROUTINE psi_init - RETURN + INTEGER :: N + REAL :: zolf - END SUBROUTINE mym_condensation + 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) + ! unstable function tables + zolf = -float(n)*0.01 + psim_unstab(n)=psim_unstable_full(zolf) + psih_unstab(n)=psih_unstable_full(zolf) + ENDDO + + END SUBROUTINE psi_init ! ================================================================== +! ... integrated similarity functions ... +! + REAL function psim_stable_full(zolf) + REAL :: zolf + + !psim_stable_full=-6.1*log(zolf+(1+zolf**2.5)**(1./2.5)) + psim_stable_full=-6.1*log(zolf+(1+zolf**2.5)**0.4) + 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)) + psih_stable_full=-5.3*log(zolf+(1+zolf**1.1)**0.9090909090909090909) + + return + end function + + REAL function psim_unstable_full(zolf) + REAL :: zolf,x,ym,psimc,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.) + psimk=2.*ALOG(0.5*(1+X))+ALOG(0.5*(1+X*X))-2.*ATAN(X)+2.*atan1 + + ym=(1.-10.*zolf)**onethird + !psimc=(3./2.)*log((ym**2.+ym+1.)/3.)-sqrt(3.)*ATAN((2.*ym+1)/sqrt(3.))+4.*ATAN(1.)/sqrt(3.) + psimc=1.5*log((ym**2 + ym+1.)*onethird)-sqrt3*ATAN((2.*ym+1)/sqrt3)+4.*atan1/sqrt3 + + psim_unstable_full=(psimk+zolf**2*(psimc))/(1+zolf**2.) + + return + end function + + REAL function psih_unstable_full(zolf) + REAL :: zolf,y,yh,psihc,psihk + + y=(1.-16.*zolf)**.5 + !psihk=2.*log((1+y)/2.) + psihk=2.*log((1+y)*0.5) + + yh=(1.-34.*zolf)**onethird + !psihc=(3./2.)*log((yh**2.+yh+1.)/3.)-sqrt(3.)*ATAN((2.*yh+1)/sqrt(3.))+4.*ATAN(1.)/sqrt(3.) + psihc=1.5*log((yh**2.+yh+1.)*onethird)-sqrt3*ATAN((2.*yh+1)/sqrt3)+4.*atan1/sqrt3 + + psih_unstable_full=(psihk+zolf**2*(psihc))/(1+zolf**2) + + return + end function +!================================================================= +! look-up table functions +!================================================================= + REAL function psim_stable(zolf) + integer :: nzol + real :: rzol,zolf + + nzol = int(zolf*100.) + rzol = zolf*100. - nzol + if(nzol+1 .le. 1000)then + psim_stable = psim_stab(nzol) + rzol*(psim_stab(nzol+1)-psim_stab(nzol)) + else + psim_stable = psim_stable_full(zolf) + endif + + return + end function + + REAL function psih_stable(zolf) + integer :: nzol + real :: rzol,zolf + + nzol = int(zolf*100.) + rzol = zolf*100. - nzol + if(nzol+1 .le. 1000)then + psih_stable = psih_stab(nzol) + rzol*(psih_stab(nzol+1)-psih_stab(nzol)) + else + psih_stable = psih_stable_full(zolf) + endif + + return + end function + + REAL function psim_unstable(zolf) + integer :: nzol + real :: rzol,zolf + + nzol = int(-zolf*100.) + rzol = -zolf*100. - nzol + if(nzol+1 .le. 1000)then + psim_unstable = psim_unstab(nzol) + rzol*(psim_unstab(nzol+1)-psim_unstab(nzol)) + else + psim_unstable = psim_unstable_full(zolf) + endif + + return + end function + + REAL function psih_unstable(zolf) + integer :: nzol + real :: rzol,zolf + + nzol = int(-zolf*100.) + rzol = -zolf*100. - nzol + if(nzol+1 .le. 1000)then + psih_unstable = psih_unstab(nzol) + rzol*(psih_unstab(nzol+1)-psih_unstab(nzol)) + else + psih_unstable = psih_unstable_full(zolf) + endif + + return + end function +!======================================================================== END MODULE module_sf_mynn