diff --git a/Registry/Registry.EM_COMMON b/Registry/Registry.EM_COMMON index 64dfddc9c1..447ea7d7b2 100644 --- a/Registry/Registry.EM_COMMON +++ b/Registry/Registry.EM_COMMON @@ -2484,6 +2484,9 @@ rconfig integer bl_pbl_physics namelist,physics max_domains -1 rconfig integer tke_budget namelist,physics max_domains 0 rh "tke_budget" "" "" rconfig integer ysu_topdown_pblmix namelist,physics 1 1 rh "ysu_topdown_pblmix" "" "" rconfig integer shinhong_tke_diag namelist,physics max_domains 0 rh "shinhong_tke_diag" "" "" +rconfig logical shinhong_scu_mixing namelist,physics max_domains .false. rh "shinhong_scu_mixing" "strato-cumulus(scu ) downdraft mixing for shinhong pbl" "" +rconfig logical shinhong_nonlocal_flux namelist,physics max_domains .false. rh "shinhong_nonlocal_flux" ".false.:counter-gradient (ysu), .true.:shinhong-LES type nonlocal flux" "" +rconfig logical shinhong_ke_dissipation namelist,physics max_domains .false. rh "shinhong_ke_dissipation" "additional heating due to kinetic energy dissipation" "" rconfig logical bl_mynn_tkeadvect namelist,physics max_domains .false. rh "bl_mynn_tkeadvect" "" "" rconfig integer bl_mynn_cloudpdf namelist,physics 1 2 irh "bl_mynn_cloudpdf" "0:original, 1:Kuwano, 2:Chaboreau-Bechtold" "" rconfig integer bl_mynn_mixlength namelist,physics 1 1 irh "bl_mynn_mixlength" "0:original,1:operational,2:new blending&cloud mix length" "" diff --git a/dyn_em/module_first_rk_step_part1.F b/dyn_em/module_first_rk_step_part1.F index 860fac46a8..b88c578a42 100644 --- a/dyn_em/module_first_rk_step_part1.F +++ b/dyn_em/module_first_rk_step_part1.F @@ -1195,6 +1195,9 @@ SUBROUTINE first_rk_step_part1 ( grid , config_flags & & ,ZNT=grid%znt & & ,ysu_topdown_pblmix=config_flags%ysu_topdown_pblmix & & ,shinhong_tke_diag=config_flags%shinhong_tke_diag & + & ,shinhong_scu_mixing=config_flags%shinhong_scu_mixing & + & ,shinhong_nonlocal_flux=config_flags%shinhong_nonlocal_flux & + & ,shinhong_ke_dissipation=config_flags%shinhong_ke_dissipation & ! paj: topo_wind & ,CTOPO=grid%ctopo,CTOPO2=grid%ctopo2 & ! variables added for BEP diff --git a/main/depend.common b/main/depend.common index 1df2c2b66a..49e166049e 100644 --- a/main/depend.common +++ b/main/depend.common @@ -627,6 +627,10 @@ module_madwrf.o: \ module_bl_ysu.o: \ ccpp_kind_types.o \ physics_mmm/bl_ysu.o + +module_bl_shinhong.o: \ + ccpp_kind_types.o \ + physics_mmm/bl_shinhong.o module_bl_myjpbl.o: \ diff --git a/phys/Makefile b/phys/Makefile index 5c2a345298..6988b5cddf 100644 --- a/phys/Makefile +++ b/phys/Makefile @@ -255,8 +255,9 @@ PHYSMMM_MODULES = \ physics_mmm/mp_wsm6.o \ physics_mmm/mp_wsm6_effectRad.o \ physics_mmm/mp_radar.o \ - physics_mmm/bl_gwdo.o \ - physics_mmm/bl_ysu.o + physics_mmm/bl_gwdo.o \ + physics_mmm/bl_shinhong.o \ + physics_mmm/bl_ysu.o NOAHMP_MODULES = \ noahmp/drivers/wrf/NoahmpWRFinitMod.o \ diff --git a/phys/module_bl_shinhong.F b/phys/module_bl_shinhong.F index ac95882eb4..e5fef311c8 100644 --- a/phys/module_bl_shinhong.F +++ b/phys/module_bl_shinhong.F @@ -1,33 +1,50 @@ -!WRF:model_layer:physics -! -module module_bl_shinhong -! -contains -! -!------------------------------------------------------------------------------- -! - subroutine shinhong(u3d,v3d,th3d,t3d,qv3d,qc3d,qi3d,p3d,p3di,pi3d, & +#define NEED_B4B_DURING_CCPP_TESTING 1 +!================================================================================================================= + module module_bl_shinhong + use ccpp_kind_types,only: kind_phys + use bl_shinhong + + + implicit none + private + public:: shinhong + + + contains + + +!================================================================================================================= + subroutine shinhong(u3d,v3d,t3d,qv3d,qc3d,qi3d,p3d,p3di,pi3d, & rublten,rvblten,rthblten, & - rqvblten,rqcblten,rqiblten,flag_qi, & + rqvblten,rqcblten,rqiblten,flag_qc,flag_qi, & cp,g,rovcp,rd,rovg,ep1,ep2,karman,xlv,rv, & dz8w,psfc, & - znu,znw,p_top, & znt,ust,hpbl,psim,psih, & xland,hfx,qfx,wspd,br, & dt,kpbl2d, & - exch_h, & + exch_h,exch_m, & + wstar,delta, & + shinhong_nonlocal_flux, & + tke,el,corf, & u10,v10, & - shinhong_tke_diag,tke_pbl,el_pbl,corf, & - dx,dy, & + uoce,voce, & + rthraten,shinhong_scu_mixing, & + shinhong_dissi_heating, & + ctopo,ctopo2,dx, & + idiff,flag_bep,frc_urb2d, & + a_u_bep,a_v_bep,a_t_bep, & + a_q_bep, & + a_e_bep,b_u_bep,b_v_bep, & + b_t_bep,b_q_bep, & + b_e_bep,dlg_bep, & + dl_u_bep,sf_bep,vl_bep, & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte, & - !optional - ctopo,ctopo2, & - wstar,delta, & - regime ) + errmsg,errflg & + ) !------------------------------------------------------------------------------- - implicit none + implicit none !------------------------------------------------------------------------------- !-- u3d 3d u-velocity interpolated to theta points (m/s) !-- v3d 3d v-velocity interpolated to theta points (m/s) @@ -40,6 +57,7 @@ subroutine shinhong(u3d,v3d,th3d,t3d,qv3d,qc3d,qi3d,p3d,p3di,pi3d, & !-- p3d 3d pressure (pa) !-- p3di 3d pressure (pa) at interface level !-- pi3d 3d exner function (dimensionless) +!-- rr3d 3d dry air density (kg/m^3) !-- rublten u tendency due to ! pbl parameterization (m/s/s) !-- rvblten v tendency due to @@ -57,29 +75,47 @@ subroutine shinhong(u3d,v3d,th3d,t3d,qv3d,qc3d,qi3d,p3d,p3di,pi3d, & !-- rovcp r/cp !-- rd gas constant for dry air (j/kg/k) !-- rovg r/g -!-- ep1 constant for virtual temperature (r_v/r_d - 1) -!-- ep2 constant for specific humidity calculation -!-- karman von karman constant +!-- dz8w dz between full levels (m) !-- xlv latent heat of vaporization (j/kg) !-- rv gas constant for water vapor (j/kg/k) -!-- dz8w dz between full levels (m) !-- psfc pressure at the surface (pa) -!-- znu eta values on half (mass) levels -!-- znw eta values on full (w) levels -!-- p_top pressure top of the model (pa) !-- znt roughness length (m) -!-- ust u* in similarity theory (m/s) +!-- ust u* in similarity theory (m/s) !-- hpbl pbl height (m) !-- psim similarity stability function for momentum !-- psih similarity stability function for heat !-- xland land mask (1 for land, 2 for water) -!-- hfx upward heat flux at the surface (w/m^2) -!-- qfx upward moisture flux at the surface (kg/m^2/s) +!-- hfx upward heat flux at the surface (w/m^2) +!-- qfx upward moisture flux at the surface (kg/m^2/s) !-- wspd wind speed at lowest model level (m/s) -!-- br bulk richardson number in surface layer !-- u10 u-wind speed at 10 m (m/s) !-- v10 v-wind speed at 10 m (m/s) +!-- uoce sea surface zonal currents (m s-1) +!-- voce sea surface meridional currents (m s-1) +!-- br bulk richardson number in surface layer +!-- dx model grid interval (m) !-- dt time step (s) +!-- rvovrd r_v divided by r_d (dimensionless) +!-- ep1 constant for virtual temperature (r_v/r_d - 1) +!-- ep2 constant for specific humidity calculation +!-- karman von karman constant +!-- idiff diff3d BEP/BEM+BEM diffusion flag +!-- flag_bep flag to use BEP/BEP+BEM +!-- frc_urb2d urban fraction +!-- a_u_bep BEP/BEP+BEM implicit component u-mom +!-- a_v_bep BEP/BEP+BEM implicit component v-mom +!-- a_t_bep BEP/BEP+BEM implicit component pot. temp. +!-- a_q_bep BEP/BEP+BEM implicit component vapor mixing ratio +!-- a_e_bep BEP/BEP+BEM implicit component TKE +!-- b_u_bep BEP/BEP+BEM explicit component u-mom +!-- b_v_bep BEP/BEP+BEM explicit component v-mom +!-- b_t_bep BEP/BEP+BEM explicit component pot.temp. +!-- b_q_bep BEP/BEP+BEM explicit component vapor mixing ratio +!-- b_e_bep BEP/BEP+BEM explicit component TKE +!-- dlg_bep Height above ground Martilli et al. (2002) Eq. 24 +!-- dl_u_bep modified length scale Martilli et al. (2002) Eq. 22 +!-- sf_bep fraction of vertical surface not occupied by buildings +!-- vl_bep volume fraction of grid cell not occupied by buildings !-- ids start index for i in domain !-- ide end index for i in domain !-- jds start index for j in domain @@ -100,2332 +136,364 @@ subroutine shinhong(u3d,v3d,th3d,t3d,qv3d,qc3d,qi3d,p3d,p3di,pi3d, & !-- kte end index for k in tile !------------------------------------------------------------------------------- ! - integer,parameter :: ndiff = 3 - real,parameter :: rcl = 1.0 ! integer, intent(in ) :: ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte - integer, intent(in ) :: shinhong_tke_diag + ! - real, intent(in ) :: dt,cp,g,rovcp,rovg,rd,xlv,rv - real, intent(in ) :: ep1,ep2,karman - real, intent(in ) :: dx,dy + real(kind=kind_phys), intent(in ) :: dt,cp,g,rovcp,rovg,rd,xlv,rv ! - integer, dimension( ims:ime, jms:jme ) , & - intent(out ) :: kpbl2d + real(kind=kind_phys), intent(in ) :: ep1,ep2,karman ! - real, dimension( ims:ime, kms:kme, jms:jme ) , & - intent(in ) :: u3d, & - v3d - real, dimension( ims:ime, kms:kme, jms:jme ) , & + real(kind=kind_phys), dimension( ims:ime, kms:kme, jms:jme ) , & intent(in ) :: qv3d, & qc3d, & qi3d, & p3d, & pi3d, & - th3d, & t3d, & - dz8w - real, dimension( ims:ime, kms:kme, jms:jme ) , & + dz8w, & + rthraten + real(kind=kind_phys), dimension( ims:ime, kms:kme, jms:jme ) , & intent(in ) :: p3di ! - real, dimension( ims:ime, kms:kme, jms:jme ) , & - intent(inout) :: rublten, & + real(kind=kind_phys), dimension( ims:ime, kms:kme, jms:jme ) , & + intent(out ) :: rublten, & rvblten, & rthblten, & rqvblten, & - rqcblten - real, dimension( ims:ime, kms:kme, jms:jme ) , & - intent(inout) :: exch_h - real, dimension( ims:ime, kms:kme, jms:jme ) , & - intent(inout) :: tke_pbl, & - el_pbl + rqcblten, & + rqiblten +! + real(kind=kind_phys), dimension( ims:ime, kms:kme, jms:jme ) , & + intent(out ) :: exch_h, & + exch_m + real(kind=kind_phys), dimension( ims:ime, kms:kme, jms:jme ) , & + intent(out ) :: tke, & + el + real(kind=kind_phys), dimension( ims:ime, jms:jme ) , & + intent(out ) :: wstar + real(kind=kind_phys), dimension( ims:ime, jms:jme ) , & + intent(out ) :: delta + real(kind=kind_phys), dimension( ims:ime, jms:jme ) , & + intent(inout) :: u10, & + v10 + real(kind=kind_phys), dimension( ims:ime, jms:jme ) , & + intent(in ) :: uoce, & + voce ! - real, dimension( ims:ime, jms:jme ) , & + real(kind=kind_phys), dimension( ims:ime, jms:jme ) , & intent(in ) :: xland, & hfx, & qfx, & - corf, & br, & + dx, & + corf, & psfc - real, dimension( ims:ime, jms:jme ) , & + real(kind=kind_phys), dimension( ims:ime, jms:jme ) , & intent(in ) :: & psim, & psih -! - real, dimension( ims:ime, jms:jme ) , & - intent(inout) :: u10, & - v10 - real, dimension( ims:ime, jms:jme ) , & - intent(inout) :: znt, & + real(kind=kind_phys), dimension( ims:ime, jms:jme ) , & + intent(in ) :: znt, & ust, & - hpbl, & - wspd -! - logical, intent(in) :: flag_qi -! -! optional -! - real, dimension( ims:ime, kms:kme, jms:jme ) , & - intent(inout), optional :: rqiblten -! - real, dimension( ims:ime, jms:jme ) , & - intent(inout), optional :: wstar, & - delta - real, dimension( ims:ime, jms:jme ) , & - intent(inout), optional :: regime -! - real, dimension( ims:ime, jms:jme ) , & - intent(in ), optional :: ctopo, & - ctopo2 -! - real, dimension( kms:kme ) , & - intent(in ), optional :: znu, & - znw -! - real, optional, intent(in ) :: p_top -! -! local -! - integer :: i,j,k - real, dimension( its:ite, kts:kte*ndiff ) :: rqvbl2dt, & - qv2d - real, dimension( its:ite, kts:kte ) :: pdh - real, dimension( its:ite, kts:kte+1 ) :: pdhi - real, dimension( its:ite ) :: & - dusfc, & - dvsfc, & - dtsfc, & - dqsfc -! - qv2d(its:ite,:) = 0.0 -! - do j = jts,jte - do k = kts,kte+1 - do i = its,ite - if(k.le.kte)pdh(i,k) = p3d(i,k,j) - pdhi(i,k) = p3di(i,k,j) - enddo - enddo - do k = kts,kte - do i = its,ite - qv2d(i,k) = qv3d(i,k,j) - qv2d(i,k+kte) = qc3d(i,k,j) - if(present(rqiblten)) qv2d(i,k+kte+kte) = qi3d(i,k,j) - enddo - enddo -! - call shinhong2d(J=j,ux=u3d(ims,kms,j),vx=v3d(ims,kms,j) & - ,tx=t3d(ims,kms,j) & - ,qx=qv2d(its,kts) & - ,p2d=pdh(its,kts),p2di=pdhi(its,kts) & - ,pi2d=pi3d(ims,kms,j) & - ,utnp=rublten(ims,kms,j),vtnp=rvblten(ims,kms,j) & - ,ttnp=rthblten(ims,kms,j),qtnp=rqvbl2dt(its,kts),ndiff=ndiff & - ,cp=cp,g=g,rovcp=rovcp,rd=rd,rovg=rovg & - ,xlv=xlv,rv=rv & - ,ep1=ep1,ep2=ep2,karman=karman & - ,dz8w2d=dz8w(ims,kms,j) & - ,psfcpa=psfc(ims,j),znt=znt(ims,j),ust=ust(ims,j) & - ,hpbl=hpbl(ims,j) & - ,regime=regime(ims,j),psim=psim(ims,j) & - ,psih=psih(ims,j),xland=xland(ims,j) & - ,hfx=hfx(ims,j),qfx=qfx(ims,j) & - ,wspd=wspd(ims,j),br=br(ims,j) & - ,dusfc=dusfc,dvsfc=dvsfc,dtsfc=dtsfc,dqsfc=dqsfc & - ,dt=dt,rcl=1.0,kpbl1d=kpbl2d(ims,j) & - ,exch_hx=exch_h(ims,kms,j) & - ,wstar=wstar(ims,j) & - ,delta=delta(ims,j) & - ,u10=u10(ims,j),v10=v10(ims,j) & - ,ctopo=ctopo(ims,j),ctopo2=ctopo2(ims,j) & - ,shinhong_tke_diag=shinhong_tke_diag & - ,tke=tke_pbl(ims,kms,j),el_pbl=el_pbl(ims,kms,j) & - ,corf=corf(ims,j) & - ,dx=dx,dy=dy & - ,ids=ids,ide=ide, jds=jds,jde=jde, kds=kds,kde=kde & - ,ims=ims,ime=ime, jms=jms,jme=jme, kms=kms,kme=kme & - ,its=its,ite=ite, jts=jts,jte=jte, kts=kts,kte=kte ) -! - do k = kts,kte - do i = its,ite - rthblten(i,k,j) = rthblten(i,k,j)/pi3d(i,k,j) - rqvblten(i,k,j) = rqvbl2dt(i,k) - rqcblten(i,k,j) = rqvbl2dt(i,k+kte) - if(present(rqiblten)) rqiblten(i,k,j) = rqvbl2dt(i,k+kte+kte) - enddo - enddo -! - enddo -! - end subroutine shinhong -!------------------------------------------------------------------------------- -! -!------------------------------------------------------------------------------- - subroutine shinhong2d(j,ux,vx,tx,qx,p2d,p2di,pi2d, & - utnp,vtnp,ttnp,qtnp,ndiff, & - cp,g,rovcp,rd,rovg,ep1,ep2,karman,xlv,rv, & - dz8w2d,psfcpa, & - znt,ust,hpbl,psim,psih, & - xland,hfx,qfx,wspd,br, & - dusfc,dvsfc,dtsfc,dqsfc, & - dt,rcl,kpbl1d, & - exch_hx, & - wstar,delta, & - shinhong_tke_diag, & - tke,el_pbl,corf, & - u10,v10, & - ctopo,ctopo2, & - dx,dy, & - ids,ide, jds,jde, kds,kde, & - ims,ime, jms,jme, kms,kme, & - its,ite, jts,jte, kts,kte, & - !optional - regime ) -!------------------------------------------------------------------------------- - implicit none -!------------------------------------------------------------------------------- -! -! the shinhongpbl (shin and hong 2015) is based on the les study of shin -! and hong (2013). the major ingredients of the shinhongpbl are -! 1) the prescribed nonlocal heat transport profile fit to the les and -! 2) inclusion of explicit scale dependency functions for vertical -! transport in convective pbl. -! so, the shinhongpbl works at the gray zone resolution of convective pbl. -! note that honnert et al. (2011) first suggested explicit scale dependency -! function, and shin and hong (2013) further classified the function by -! stability (u*/w*) in convective pbl and calculated the function for -! nonlocal and local transport separately. -! vertical mixing in the stable boundary layer and free atmosphere follows -! hong (2010) and hong et al. (2006), same as the ysupbl scheme. + wspd + real(kind=kind_phys), dimension( ims:ime, jms:jme ) , & + intent(out ) :: hpbl ! -! shinhongpbl: -! coded and implemented by hyeyum hailey shin (ncar) -! summer 2014 -! -! ysupbl: -! coded by song-you hong (yonsei university) and implemented by -! song-you hong (yonsei university) and jimy dudhia (ncar) -! summer 2002 -! -! references: -! shin and hong (2015) mon. wea. rev. -! shin and hong (2013) j. atmos. sci. -! honnert, masson, and couvreux (2011) j. atmos. sci. -! hong (2010) quart. j. roy. met. soc -! hong, noh, and dudhia (2006), mon. wea. rev. -! -!------------------------------------------------------------------------------- -! - real,parameter :: xkzminm = 0.1,xkzminh = 0.01 - real,parameter :: xkzmax = 1000.,rimin = -100. - real,parameter :: rlam = 30.,prmin = 0.25,prmax = 4. - real,parameter :: brcr_ub = 0.0,brcr_sb = 0.25,cori = 1.e-4 - real,parameter :: afac = 6.8,bfac = 6.8,pfac = 2.0,pfac_q = 2.0 - real,parameter :: phifac = 8.,sfcfrac = 0.1 - real,parameter :: d1 = 0.02, d2 = 0.05, d3 = 0.001 - real,parameter :: h1 = 0.33333333, h2 = 0.6666667 - real,parameter :: ckz = 0.001,zfmin = 1.e-8,aphi5 = 5.,aphi16 = 16. - real,parameter :: tmin=1.e-2 - real,parameter :: gamcrt = 3.,gamcrq = 2.e-3 - real,parameter :: xka = 2.4e-5 - integer,parameter :: imvdif = 1 -! -! tunable parameters for tke -! - real,parameter :: epsq2l = 0.01,c_1 = 1.0,gamcre = 0.224 -! -! tunable parameters for prescribed nonlocal transport profile -! - real,parameter :: mltop = 1.0,sfcfracn1 = 0.075 - real,parameter :: nlfrac = 0.7,enlfrac = -0.4 - real,parameter :: a11 = 1.0,a12 = -1.15 - real,parameter :: ezfac = 1.5 - real,parameter :: cpent = -0.4,rigsmax = 100. - real,parameter :: entfmin = 1.0, entfmax = 5.0 -! - integer, intent(in ) :: ids,ide, jds,jde, kds,kde, & - ims,ime, jms,jme, kms,kme, & - its,ite, jts,jte, kts,kte, & - j,ndiff - integer, intent(in ) :: shinhong_tke_diag -! - real, intent(in ) :: dt,rcl,cp,g,rovcp,rovg,rd,xlv,rv - real, intent(in ) :: ep1,ep2,karman - real, intent(in ) :: dx,dy -! - integer, dimension( ims:ime ) , & - intent(out ) :: kpbl1d -! - real, dimension( ims:ime, kms:kme ) , & - intent(in ) :: dz8w2d, & - pi2d - real, dimension( ims:ime, kms:kme ) , & - intent(in ) :: ux, & - vx - real, dimension( ims:ime, kms:kme ) , & - intent(in ) :: tx -! - real, dimension( its:ite, kts:kte*ndiff ) , & - intent(in ) :: qx - real, dimension( its:ite, kts:kte+1 ) , & - intent(in ) :: p2di - real, dimension( its:ite, kts:kte ) , & - intent(in ) :: p2d + real(kind=kind_phys), dimension( ims:ime, kms:kme, jms:jme ) , & + intent(in ) :: u3d, & + v3d ! - real, dimension( ims:ime, kms:kme ) , & - intent(inout) :: utnp, & - vtnp, & - ttnp - real, dimension( ims:ime, kms:kme ) , & - intent(inout) :: exch_hx - real, dimension( ims:ime, kms:kme ) , & - intent(inout) :: tke, & - el_pbl - real, dimension( its:ite, kts:kte*ndiff ) , & - intent(inout) :: qtnp + integer, dimension( ims:ime, jms:jme ) , & + intent(out ) :: kpbl2d ! - real, dimension( ims:ime ) , & - intent(in ) :: xland, & - hfx, & - qfx - real, dimension( ims:ime ) , & - intent(in ) :: br, & - psim, & - psih, & - psfcpa - real, dimension( ims:ime ) , & - intent(in ) :: corf + logical, intent(in) :: flag_qc, & + flag_qi + logical, intent(in) :: shinhong_nonlocal_flux, & + shinhong_scu_mixing, & + shinhong_dissi_heating ! - real, dimension( ims:ime ) , & - intent(inout) :: ust, & - hpbl, & - znt - real, dimension( ims:ime ) , & - intent(inout) :: wspd - real, dimension( ims:ime ) , & - intent(inout) :: u10, & - v10 + integer, intent(in) :: idiff + logical, intent(in) :: flag_bep + real(kind=kind_phys), dimension( ims:ime, kms:kme, jms:jme ) , & + optional , & + intent(in) :: a_u_bep, & + a_v_bep,a_t_bep, & + a_e_bep,b_u_bep, & + a_q_bep,b_q_bep, & + b_v_bep,b_t_bep, & + b_e_bep,dlg_bep, & + dl_u_bep, & + vl_bep,sf_bep + real(kind=kind_phys), dimension(ims:ime,jms:jme) , & + optional , & + intent(in) :: frc_urb2d ! - real, dimension( ims:ime ) , & + real(kind=kind_phys), dimension( ims:ime, jms:jme ) , & optional , & intent(in ) :: ctopo, & ctopo2 - real, dimension( ims:ime ) , & - optional , & - intent(inout) :: regime - real, dimension( its:ite ) , & - intent(out ) :: wstar, & - delta -! -! local vars -! - integer :: n,i,k,l,ic,is,nwmass - integer :: klpbl, kqc, kqi - integer :: lmh,lmxl -! - real :: dt2,rdt,spdk2,fm,fh,hol1,gamfac,vpert,prnum,prnum0 - real :: ss,ri,qmean,tmean,alpha,chi,zk,rl2,dk,sri - real :: brint,dtodsd,dtodsu,rdz,dsdzt,dsdzq,dsdz2,rlamdz - real :: utend,vtend,ttend,qtend - real :: dtstep,govrthv - real :: cont, conq, conw, conwrc - real :: delxy,pu1,pth1,pq1 - real :: dex,hgame_c - real :: zfacdx - real :: amf1,amf2,bmf2,amf3,bmf3,amf4,bmf4,sflux0,snlflux0 - real :: mlfrac,ezfrac,sfcfracn - real :: uwst,uwstx,csfac - real :: prnumfac,bfx0,hfx0,qfx0,delb,dux,dvx, & - dsdzu,dsdzv,wm3,dthx,dqx,wspd10,ross,tem1,dsig,tvcon,conpr, & - prfac,prfac2,phim8z - real :: cenlfrac -! - integer, dimension( its:ite ) :: kpbl - real, dimension( its:ite ) :: hol - real, dimension( its:ite ) :: deltaoh - real, dimension( its:ite ) :: rigs, & - enlfrac2, & - cslen - real, dimension( its:ite ) :: & - rhox, & - govrth, & - zl1,thermal, & - wscale, & - hgamt,hgamq, & - brdn,brup, & - phim,phih, & - dusfc,dvsfc, & - dtsfc,dqsfc, & - prpbl, & - wspd1 - real, dimension( its:ite ) :: & - ust3, & - wstar3, & - hgamu,hgamv, & - wm2, we, & - bfxpbl, & - hfxpbl,qfxpbl, & - ufxpbl,vfxpbl, & - dthvx - real, dimension( its:ite ) :: & - brcr, & - sflux, & - zol1, & - brcr_sbro - real, dimension( its:ite ) :: & - efxpbl, & - hpbl_cbl, & - epshol, & - ct -! - real, dimension( its:ite, kts:kte ) :: xkzm,xkzh, & - f1,f2, & - r1,r2, & - ad,au, & - cu, & - al, & - xkzq, & - zfac - real, dimension( its:ite, kts:kte ) :: & - thx,thvx, & - del, & - dza, & - dzq, & - xkzom, & - xkzoh, & - za - real, dimension( its:ite, kts:kte ) :: & - wscalek - real, dimension( its:ite, kts:kte ) :: & - xkzml,xkzhl, & - zfacent,entfac - real, dimension( its:ite, kts:kte ) :: & - mf, & - zfacmf, & - entfacmf - real, dimension( its:ite, kts:kte ) :: & - q2x, & - hgame2d, & - tflux_e, & - qflux_e, & - tvflux_e - real, dimension( its:ite, kts:kte+1 ) :: zq - real, dimension( its:ite, kts:kte, ndiff ) :: r3,f3 -! - real, dimension( kts:kte ) :: & - uxk,vxk, & - txk,thxk,thvxk, & - q2xk, & - hgame - real, dimension( kts:kte ) :: & - ps1d,pb1d,eps1d,pt1d, & - xkze1d,eflx_l1d,eflx_nl1d, & - ptke1 - real, dimension( kts+1:kte ) :: & - s2,gh,rig,el, & - akmk,akhk, & - mfk,ufxpblk,vfxpblk,qfxpblk - real, dimension( kts:kte+1 ) :: zqk - real, dimension( kts:kte*ndiff ) :: qxk -! - logical, dimension( its:ite ) :: pblflg, & - sfcflg, & - stable - logical, dimension( ndiff ) :: ifvmix -! -!------------------------------------------------------------------------------- -! - klpbl = kte - lmh = 1 - lmxl = 1 -! - cont=cp/g - conq=xlv/g - conw=1./g - conwrc = conw*sqrt(rcl) - conpr = bfac*karman*sfcfrac -! -! k-start index for cloud and rain -! - kqc = 1 + kte - kqi = 1 + kte*2 - nwmass = 3 - ifvmix(:) = .true. -! - do k = kts,kte - do i = its,ite - thx(i,k) = tx(i,k)/pi2d(i,k) - enddo - enddo -! - do k = kts,kte - do i = its,ite - tvcon = (1.+ep1*qx(i,k)) - thvx(i,k) = thx(i,k)*tvcon - enddo - enddo -! - do i = its,ite - tvcon = (1.+ep1*qx(i,1)) - rhox(i) = psfcpa(i)/(rd*tx(i,1)*tvcon) - govrth(i) = g/thx(i,1) - enddo -! -!-----compute the height of full- and half-sigma levels above ground -! level, and the layer thicknesses. -! - do i = its,ite - zq(i,1) = 0. - enddo -! - do k = kts,kte - do i = its,ite - zq(i,k+1) = dz8w2d(i,k)+zq(i,k) - enddo - enddo -! - do k = kts,kte - do i = its,ite - za(i,k) = 0.5*(zq(i,k)+zq(i,k+1)) - dzq(i,k) = zq(i,k+1)-zq(i,k) - del(i,k) = p2di(i,k)-p2di(i,k+1) - enddo - enddo -! - do i = its,ite - dza(i,1) = za(i,1) - enddo -! - do k = kts+1,kte - do i = its,ite - dza(i,k) = za(i,k)-za(i,k-1) - enddo - enddo -! -! -!-----initialize vertical tendencies and -! - utnp(its:ite,:) = 0. - vtnp(its:ite,:) = 0. - ttnp(its:ite,:) = 0. - qtnp(its:ite,:) = 0. -! - do i = its,ite - wspd1(i) = sqrt(ux(i,1)*ux(i,1)+vx(i,1)*vx(i,1))+1.e-9 - enddo -! -!---- compute vertical diffusion -! -! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -! compute preliminary variables -! - dtstep = dt - dt2 = 2.*dtstep - rdt = 1./dt2 -! - do i = its,ite - bfxpbl(i) = 0.0 - hfxpbl(i) = 0.0 - qfxpbl(i) = 0.0 - ufxpbl(i) = 0.0 - vfxpbl(i) = 0.0 - hgamu(i) = 0.0 - hgamv(i) = 0.0 - delta(i) = 0.0 - enddo ! - do i = its,ite - efxpbl(i) = 0.0 - hpbl_cbl(i) = 0.0 - epshol(i) = 0.0 - ct(i) = 0.0 - enddo -! - do i = its,ite - deltaoh(i) = 0.0 - rigs(i) = 0.0 - enlfrac2(i) = 0.0 - cslen(i) = 0.0 - enddo -! - do k = kts,klpbl - do i = its,ite - wscalek(i,k) = 0.0 - enddo - enddo -! - do k = kts,klpbl - do i = its,ite - zfac(i,k) = 0.0 - enddo - enddo -! - do k = kts,kte - do i = its,ite - q2x(i,k) = 2.*tke(i,k) - enddo - enddo -! - do k = kts,kte - do i = its,ite - el_pbl(i,k) = 0.0 - hgame2d(i,k) = 0.0 - tflux_e(i,k) = 0.0 - qflux_e(i,k) = 0.0 - tvflux_e(i,k) = 0.0 - enddo - enddo -! - do k = kts,kte - do i = its,ite - mf(i,k) = 0.0 - zfacmf(i,k) = 0.0 - enddo - enddo -! - do k = kts,klpbl-1 - do i = its,ite - xkzom(i,k) = xkzminm - xkzoh(i,k) = xkzminh - enddo - enddo -! - do i = its,ite - dusfc(i) = 0. - dvsfc(i) = 0. - dtsfc(i) = 0. - dqsfc(i) = 0. - enddo -! - do i = its,ite - hgamt(i) = 0. - hgamq(i) = 0. - wscale(i) = 0. - kpbl(i) = 1 - hpbl(i) = zq(i,1) - hpbl_cbl(i) = zq(i,1) - zl1(i) = za(i,1) - thermal(i)= thvx(i,1) - pblflg(i) = .true. - sfcflg(i) = .true. - sflux(i) = hfx(i)/rhox(i)/cp + qfx(i)/rhox(i)*ep1*thx(i,1) - if(br(i).gt.0.0) sfcflg(i) = .false. - enddo -! -! compute the first guess of pbl height -! - do i = its,ite - stable(i) = .false. - brup(i) = br(i) - brcr(i) = brcr_ub - enddo -! - do k = 2,klpbl - do i = its,ite - if(.not.stable(i))then - brdn(i) = brup(i) - spdk2 = max(ux(i,k)**2+vx(i,k)**2,1.) - brup(i) = (thvx(i,k)-thermal(i))*(g*za(i,k)/thvx(i,1))/spdk2 - kpbl(i) = k - stable(i) = brup(i).gt.brcr(i) - endif - enddo - enddo -! - do i = its,ite - k = kpbl(i) - if(brdn(i).ge.brcr(i))then - brint = 0. - elseif(brup(i).le.brcr(i))then - brint = 1. - else - brint = (brcr(i)-brdn(i))/(brup(i)-brdn(i)) - endif - hpbl(i) = za(i,k-1)+brint*(za(i,k)-za(i,k-1)) - if(hpbl(i).lt.zq(i,2)) kpbl(i) = 1 - if(kpbl(i).le.1) pblflg(i) = .false. - enddo -! - do i = its,ite - fm = psim(i) - fh = psih(i) - zol1(i) = max(br(i)*fm*fm/fh,rimin) - if(sfcflg(i))then - zol1(i) = min(zol1(i),-zfmin) - else - zol1(i) = max(zol1(i),zfmin) - endif - hol1 = zol1(i)*hpbl(i)/zl1(i)*sfcfrac - epshol(i) = hol1 - if(sfcflg(i))then - phim(i) = (1.-aphi16*hol1)**(-1./4.) - phih(i) = (1.-aphi16*hol1)**(-1./2.) - bfx0 = max(sflux(i),0.) - hfx0 = max(hfx(i)/rhox(i)/cp,0.) - qfx0 = max(ep1*thx(i,1)*qfx(i)/rhox(i),0.) - wstar3(i) = (govrth(i)*bfx0*hpbl(i)) - wstar(i) = (wstar3(i))**h1 - else - phim(i) = (1.+aphi5*hol1) - phih(i) = phim(i) - wstar(i) = 0. - wstar3(i) = 0. - endif - ust3(i) = ust(i)**3. - wscale(i) = (ust3(i)+phifac*karman*wstar3(i)*0.5)**h1 - wscale(i) = min(wscale(i),ust(i)*aphi16) - wscale(i) = max(wscale(i),ust(i)/aphi5) - enddo -! -! compute the surface variables for pbl height estimation -! under unstable conditions -! - do i = its,ite - if(sfcflg(i).and.sflux(i).gt.0.0)then - gamfac = bfac/rhox(i)/wscale(i) - hgamt(i) = min(gamfac*hfx(i)/cp,gamcrt) - hgamq(i) = min(gamfac*qfx(i),gamcrq) - vpert = (hgamt(i)+ep1*thx(i,1)*hgamq(i))/bfac*afac - thermal(i) = thermal(i)+max(vpert,0.)*min(za(i,1)/(sfcfrac*hpbl(i)),1.0) - hgamt(i) = max(hgamt(i),0.0) - hgamq(i) = max(hgamq(i),0.0) - brint = -15.9*ust(i)*ust(i)/wspd(i)*wstar3(i)/(wscale(i)**4.) - hgamu(i) = brint*ux(i,1) - hgamv(i) = brint*vx(i,1) - else - pblflg(i) = .false. - endif - enddo -! -! enhance the pbl height by considering the thermal -! - do i = its,ite - if(pblflg(i))then - kpbl(i) = 1 - hpbl(i) = zq(i,1) - endif - enddo -! - do i = its,ite - if(pblflg(i))then - stable(i) = .false. - brup(i) = br(i) - brcr(i) = brcr_ub - endif - enddo -! - do k = 2,klpbl - do i = its,ite - if(.not.stable(i).and.pblflg(i))then - brdn(i) = brup(i) - spdk2 = max(ux(i,k)**2+vx(i,k)**2,1.) - brup(i) = (thvx(i,k)-thermal(i))*(g*za(i,k)/thvx(i,1))/spdk2 - kpbl(i) = k - stable(i) = brup(i).gt.brcr(i) - endif - enddo - enddo -! - do i = its,ite - if(pblflg(i)) then - k = kpbl(i) - if(brdn(i).ge.brcr(i))then - brint = 0. - elseif(brup(i).le.brcr(i))then - brint = 1. - else - brint = (brcr(i)-brdn(i))/(brup(i)-brdn(i)) - endif - hpbl(i) = za(i,k-1)+brint*(za(i,k)-za(i,k-1)) - if(hpbl(i).lt.zq(i,2)) kpbl(i) = 1 - if(kpbl(i).le.1) pblflg(i) = .false. - if (wstar(i) .ne. 0) then - uwst = abs(ust(i)/wstar(i)-0.5) - uwstx = -80.*uwst+14. - csfac = 0.5*(tanh(uwstx)+3.) - else - csfac = 1 - endif - cslen(i) = csfac*hpbl(i) - endif - enddo -! -! stable boundary layer -! - do i = its,ite - hpbl_cbl(i) = hpbl(i) - if((.not.sfcflg(i)).and.hpbl(i).lt.zq(i,2)) then - brup(i) = br(i) - stable(i) = .false. - else - stable(i) = .true. - endif - enddo -! - do i = its,ite - if((.not.stable(i)).and.((xland(i)-1.5).ge.0))then - wspd10 = u10(i)*u10(i) + v10(i)*v10(i) - wspd10 = sqrt(wspd10) - ross = wspd10 / (cori*znt(i)) - brcr_sbro(i) = min(0.16*(1.e-7*ross)**(-0.18),.3) - endif - enddo -! - do i = its,ite - if(.not.stable(i))then - if((xland(i)-1.5).ge.0)then - brcr(i) = brcr_sbro(i) - else - brcr(i) = brcr_sb - endif - endif - enddo -! - do k = 2,klpbl - do i = its,ite - if(.not.stable(i))then - brdn(i) = brup(i) - spdk2 = max(ux(i,k)**2+vx(i,k)**2,1.) - brup(i) = (thvx(i,k)-thermal(i))*(g*za(i,k)/thvx(i,1))/spdk2 - kpbl(i) = k - stable(i) = brup(i).gt.brcr(i) - endif - enddo - enddo -! - do i = its,ite - if((.not.sfcflg(i)).and.hpbl(i).lt.zq(i,2)) then - k = kpbl(i) - if(brdn(i).ge.brcr(i))then - brint = 0. - elseif(brup(i).le.brcr(i))then - brint = 1. - else - brint = (brcr(i)-brdn(i))/(brup(i)-brdn(i)) - endif - hpbl(i) = za(i,k-1)+brint*(za(i,k)-za(i,k-1)) - if(hpbl(i).lt.zq(i,2)) kpbl(i) = 1 - if(kpbl(i).le.1) pblflg(i) = .false. - endif - enddo -! -! scale dependency for nonlocal momentum and moisture transport -! - delxy=sqrt(dx*dy) -! - do i = its,ite - pu1=pu(delxy,cslen(i)) - pq1=pq(delxy,cslen(i)) - if(pblflg(i)) then - hgamu(i) = hgamu(i)*pu1 - hgamv(i) = hgamv(i)*pu1 - hgamq(i) = hgamq(i)*pq1 - endif - enddo -! -! estimate the entrainment parameters -! - delxy=sqrt(dx*dy) -! - do i = its,ite - if(pblflg(i)) then - k = kpbl(i) - 1 - prpbl(i) = 1.0 - wm3 = wstar3(i) + 5. * ust3(i) - wm2(i) = wm3**h2 - bfxpbl(i) = -0.15*thvx(i,1)/g*wm3/hpbl(i) - dthvx(i) = max(thvx(i,k+1)-thvx(i,k),tmin) - dthx = max(thx(i,k+1)-thx(i,k),tmin) - dqx = min(qx(i,k+1)-qx(i,k),0.0) - we(i) = max(bfxpbl(i)/dthvx(i),-sqrt(wm2(i))) - hfxpbl(i) = we(i)*dthx - pq1=pq(delxy,cslen(i)) - qfxpbl(i) = we(i)*dqx*pq1 -! - pu1=pu(delxy,cslen(i)) - dux = ux(i,k+1)-ux(i,k) - dvx = vx(i,k+1)-vx(i,k) - if(dux.gt.tmin) then - ufxpbl(i) = max(prpbl(i)*we(i)*dux*pu1,-ust(i)*ust(i)) - elseif(dux.lt.-tmin) then - ufxpbl(i) = min(prpbl(i)*we(i)*dux*pu1,ust(i)*ust(i)) - else - ufxpbl(i) = 0.0 - endif - if(dvx.gt.tmin) then - vfxpbl(i) = max(prpbl(i)*we(i)*dvx*pu1,-ust(i)*ust(i)) - elseif(dvx.lt.-tmin) then - vfxpbl(i) = min(prpbl(i)*we(i)*dvx*pu1,ust(i)*ust(i)) - else - vfxpbl(i) = 0.0 - endif - delb = govrth(i)*d3*hpbl(i) - delta(i) = min(d1*hpbl(i) + d2*wm2(i)/delb,100.) - delb = govrth(i)*dthvx(i) - deltaoh(i) = d1*hpbl(i) + d2*wm2(i)/delb - deltaoh(i) = max(ezfac*deltaoh(i),hpbl(i)-za(i,kpbl(i)-1)-1.) - deltaoh(i) = min(deltaoh(i) ,hpbl(i)) - if ((dux .ne. 0) .or. (dvx .ne. 0)) then - rigs(i) = govrth(i)*dthvx(i)*deltaoh(i)/(dux**2.+dvx**2.) - else - rigs(i) = rigsmax - endif - rigs(i) = max(min(rigs(i), rigsmax),rimin) - if((rigs(i).gt.0) .and. (abs(rigs(i)+cpent) .le. 1.e-6))then - cenlfrac = entfmax - else - cenlfrac = rigs(i)/(rigs(i)+cpent) - endif - cenlfrac = min(cenlfrac,entfmax) - enlfrac2(i) = max(wm3/wstar3(i)*cenlfrac, entfmin) - enlfrac2(i) = enlfrac2(i)*enlfrac - endif - enddo -! - do k = kts,klpbl - do i = its,ite - if(pblflg(i))then - entfacmf(i,k) = sqrt(((zq(i,k+1)-hpbl(i))/deltaoh(i))**2.) - endif - if(pblflg(i).and.k.ge.kpbl(i))then - entfac(i,k) = ((zq(i,k+1)-hpbl(i))/deltaoh(i))**2. - else - entfac(i,k) = 1.e30 - endif - enddo - enddo -! -! compute diffusion coefficients below pbl -! - do k = kts,klpbl - do i = its,ite - if(k.lt.kpbl(i)) then - zfac(i,k) = min(max((1.-(zq(i,k+1)-zl1(i))/(hpbl(i)-zl1(i))),zfmin),1.) - zfacent(i,k) = (1.-zfac(i,k))**3. - wscalek(i,k) = (ust3(i)+phifac*karman*wstar3(i)*(1.-zfac(i,k)))**h1 - if(sfcflg(i)) then - prfac = conpr - prfac2 = 15.9*wstar3(i)/ust3(i)/(1.+4.*karman*wstar3(i)/ust3(i)) - prnumfac = -3.*(max(zq(i,k+1)-sfcfrac*hpbl(i),0.))**2./hpbl(i)**2. - else - prfac = 0. - prfac2 = 0. - prnumfac = 0. - phim8z = 1.+aphi5*zol1(i)*zq(i,k+1)/zl1(i) - wscalek(i,k) = ust(i)/phim8z - wscalek(i,k) = max(wscalek(i,k),0.001) - endif - prnum0 = (phih(i)/phim(i)+prfac) - prnum0 = max(min(prnum0,prmax),prmin) - xkzm(i,k) = wscalek(i,k)*karman*zq(i,k+1)*zfac(i,k)**pfac - prnum = 1. + (prnum0-1.)*exp(prnumfac) - xkzq(i,k) = xkzm(i,k)/prnum*zfac(i,k)**(pfac_q-pfac) - prnum0 = prnum0/(1.+prfac2*karman*sfcfrac) - prnum = 1. + (prnum0-1.)*exp(prnumfac) - xkzh(i,k) = xkzm(i,k)/prnum - xkzm(i,k) = xkzm(i,k)+xkzom(i,k) - xkzh(i,k) = xkzh(i,k)+xkzoh(i,k) - xkzq(i,k) = xkzq(i,k)+xkzoh(i,k) - xkzm(i,k) = min(xkzm(i,k),xkzmax) - xkzh(i,k) = min(xkzh(i,k),xkzmax) - xkzq(i,k) = min(xkzq(i,k),xkzmax) - endif - enddo - enddo -! -! compute diffusion coefficients over pbl (free atmosphere) -! - do k = kts,kte-1 - do i = its,ite - if(k.ge.kpbl(i)) then - ss = ((ux(i,k+1)-ux(i,k))*(ux(i,k+1)-ux(i,k)) & - +(vx(i,k+1)-vx(i,k))*(vx(i,k+1)-vx(i,k))) & - /(dza(i,k+1)*dza(i,k+1))+1.e-9 - govrthv = g/(0.5*(thvx(i,k+1)+thvx(i,k))) - ri = govrthv*(thvx(i,k+1)-thvx(i,k))/(ss*dza(i,k+1)) -! in cloud - if(imvdif.eq.1.and.nwmass.ge.3)then - if((qx(i,kqc+k-1)+qx(i,kqi+k-1)).gt.0.01e-3 & - .and.(qx(i,kqc+k)+qx(i,kqi+k)).gt.0.01e-3) then - qmean = 0.5*(qx(i,k)+qx(i,k+1)) - tmean = 0.5*(tx(i,k)+tx(i,k+1)) - alpha = xlv*qmean/rd/tmean - chi = xlv*xlv*qmean/cp/rv/tmean/tmean - ri = (1.+alpha)*(ri-g*g/ss/tmean/cp*((chi-alpha)/(1.+chi))) - endif - endif - zk = karman*zq(i,k+1) - rlamdz = min(max(0.1*dza(i,k+1),rlam),300.) - rlamdz = min(dza(i,k+1),rlamdz) - rl2 = (zk*rlamdz/(rlamdz+zk))**2 - dk = rl2*sqrt(ss) - if(ri.lt.0.)then -! unstable regime - ri = max(ri, rimin) - sri = sqrt(-ri) - xkzm(i,k) = dk*(1+8.*(-ri)/(1+1.746*sri)) - xkzh(i,k) = dk*(1+8.*(-ri)/(1+1.286*sri)) - else -! stable regime - xkzh(i,k) = dk/(1+5.*ri)**2 - prnum = 1.0+2.1*ri - prnum = min(prnum,prmax) - xkzm(i,k) = xkzh(i,k)*prnum - endif -! - xkzm(i,k) = xkzm(i,k)+xkzom(i,k) - xkzh(i,k) = xkzh(i,k)+xkzoh(i,k) - xkzm(i,k) = min(xkzm(i,k),xkzmax) - xkzh(i,k) = min(xkzh(i,k),xkzmax) - xkzml(i,k) = xkzm(i,k) - xkzhl(i,k) = xkzh(i,k) - endif - enddo - enddo -! -! prescribe nonlocal heat transport below pbl -! - do i = its,ite - deltaoh(i) = deltaoh(i)/hpbl(i) - enddo -! - delxy=sqrt(dx*dy) - do i = its,ite - mlfrac = mltop-deltaoh(i) - ezfrac = mltop+deltaoh(i) - zfacmf(i,1) = min(max((zq(i,2)/hpbl(i)),zfmin),1.) - sfcfracn = max(sfcfracn1,zfacmf(i,1)) -! - sflux0 = (a11+a12*sfcfracn)*sflux(i) - snlflux0 = nlfrac*sflux0 - amf1 = snlflux0/sfcfracn - if (pblflg(i)) then - amf2 = -snlflux0/(mlfrac-sfcfracn) - bmf2 = -mlfrac*amf2 - endif - if ((deltaoh(i) .eq. 0) .and. (enlfrac2(i) .eq. 0)) then - amf3 = 0. - else - amf3 = snlflux0*enlfrac2(i)/deltaoh(i) - endif - bmf3 = -amf3*mlfrac - hfxpbl(i) = amf3+bmf3 - pth1=pthnl(delxy,cslen(i)) - hfxpbl(i) = hfxpbl(i)*pth1 -! - do k = kts,klpbl - zfacmf(i,k) = max((zq(i,k+1)/hpbl(i)),zfmin) - if(pblflg(i).and.k.lt.kpbl(i)) then - if(zfacmf(i,k).le.sfcfracn) then - mf(i,k) = amf1*zfacmf(i,k) - else if (zfacmf(i,k).le.mlfrac) then - mf(i,k) = amf2*zfacmf(i,k)+bmf2 - endif - mf(i,k) = mf(i,k)+hfxpbl(i)*exp(-entfacmf(i,k)) - mf(i,k) = mf(i,k)*pth1 - endif - enddo - enddo -! -! compute tridiagonal matrix elements for heat -! - do k = kts,kte - do i = its,ite - au(i,k) = 0. - al(i,k) = 0. - ad(i,k) = 0. - f1(i,k) = 0. - enddo - enddo -! - do i = its,ite - ad(i,1) = 1. - f1(i,1) = thx(i,1)-300.+hfx(i)/cont/del(i,1)*dt2 - enddo -! - delxy=sqrt(dx*dy) - do k = kts,kte-1 - do i = its,ite - dtodsd = dt2/del(i,k) - dtodsu = dt2/del(i,k+1) - dsig = p2d(i,k)-p2d(i,k+1) - rdz = 1./dza(i,k+1) - tem1 = dsig*xkzh(i,k)*rdz - if(pblflg(i).and.k.lt.kpbl(i)) then - dsdzt = tem1*(-mf(i,k)/xkzh(i,k)) - f1(i,k) = f1(i,k)+dtodsd*dsdzt - f1(i,k+1) = thx(i,k+1)-300.-dtodsu*dsdzt - elseif(pblflg(i).and.k.ge.kpbl(i).and.entfac(i,k).lt.4.6) then - xkzh(i,k) = -we(i)*dza(i,kpbl(i))*exp(-entfac(i,k)) - xkzh(i,k) = sqrt(xkzh(i,k)*xkzhl(i,k)) - xkzh(i,k) = max(xkzh(i,k),xkzoh(i,k)) - xkzh(i,k) = min(xkzh(i,k),xkzmax) - f1(i,k+1) = thx(i,k+1)-300. - else - f1(i,k+1) = thx(i,k+1)-300. - endif - tem1 = dsig*xkzh(i,k)*rdz - dsdz2 = tem1*rdz - au(i,k) = -dtodsd*dsdz2 - al(i,k) = -dtodsu*dsdz2 -! -! scale dependency for local heat transport -! - zfacdx=0.2*hpbl(i)/zq(i,k+1) - delxy=sqrt(dx*dy)*max(zfacdx,1.0) - pth1=pthl(delxy,hpbl(i)) - if(pblflg(i).and.k.lt.kpbl(i)) then - au(i,k) = au(i,k)*pth1 - al(i,k) = al(i,k)*pth1 - endif - ad(i,k) = ad(i,k)-au(i,k) - ad(i,k+1) = 1.-al(i,k) - exch_hx(i,k+1) = xkzh(i,k) - enddo - enddo -! -! copies here to avoid duplicate input args for tridin -! - do k = kts,kte - do i = its,ite - cu(i,k) = au(i,k) - r1(i,k) = f1(i,k) - enddo - enddo -! - call tridin_ysu(al,ad,cu,r1,au,f1,its,ite,kts,kte,1) -! -! recover tendencies of heat -! - do k = kte,kts,-1 - do i = its,ite - ttend = (f1(i,k)-thx(i,k)+300.)*rdt*pi2d(i,k) - ttnp(i,k) = ttnp(i,k)+ttend - dtsfc(i) = dtsfc(i)+ttend*cont*del(i,k)/pi2d(i,k) - if(k.eq.kte) then - tflux_e(i,k) = ttend*dz8w2d(i,k) - else - tflux_e(i,k) = tflux_e(i,k+1) + ttend*dz8w2d(i,k) - endif - enddo - enddo -! -! compute tridiagonal matrix elements for moisture, clouds, and gases -! - do k = kts,kte - do i = its,ite - au(i,k) = 0. - al(i,k) = 0. - ad(i,k) = 0. - enddo - enddo -! - do ic = 1,ndiff - do i = its,ite - do k = kts,kte - f3(i,k,ic) = 0. - enddo - enddo - enddo -! - do i = its,ite - ad(i,1) = 1. - f3(i,1,1) = qx(i,1)+qfx(i)*g/del(i,1)*dt2 - enddo -! - if(ndiff.ge.2) then - do ic = 2,ndiff - is = (ic-1) * kte - do i = its,ite - f3(i,1,ic) = qx(i,1+is) - enddo - enddo - endif -! - do k = kts,kte-1 - do i = its,ite - if(k.ge.kpbl(i)) then - xkzq(i,k) = xkzh(i,k) - endif - enddo - enddo -! - do k = kts,kte-1 - do i = its,ite - dtodsd = dt2/del(i,k) - dtodsu = dt2/del(i,k+1) - dsig = p2d(i,k)-p2d(i,k+1) - rdz = 1./dza(i,k+1) - tem1 = dsig*xkzq(i,k)*rdz - if(pblflg(i).and.k.lt.kpbl(i)) then - dsdzq = tem1*(-qfxpbl(i)*zfacent(i,k)/xkzq(i,k)) - f3(i,k,1) = f3(i,k,1)+dtodsd*dsdzq - f3(i,k+1,1) = qx(i,k+1)-dtodsu*dsdzq - elseif(pblflg(i).and.k.ge.kpbl(i).and.entfac(i,k).lt.4.6) then - xkzq(i,k) = -we(i)*dza(i,kpbl(i))*exp(-entfac(i,k)) - xkzq(i,k) = sqrt(xkzq(i,k)*xkzhl(i,k)) - xkzq(i,k) = max(xkzq(i,k),xkzoh(i,k)) - xkzq(i,k) = min(xkzq(i,k),xkzmax) - f3(i,k+1,1) = qx(i,k+1) - else - f3(i,k+1,1) = qx(i,k+1) - endif - tem1 = dsig*xkzq(i,k)*rdz - dsdz2 = tem1*rdz - au(i,k) = -dtodsd*dsdz2 - al(i,k) = -dtodsu*dsdz2 -! -! scale dependency for local moisture transport -! - zfacdx=0.2*hpbl(i)/zq(i,k+1) - delxy=sqrt(dx*dy)*max(zfacdx,1.0) - pq1=pq(delxy,hpbl(i)) - if(pblflg(i).and.k.lt.kpbl(i)) then - au(i,k) = au(i,k)*pq1 - al(i,k) = al(i,k)*pq1 - endif - ad(i,k) = ad(i,k)-au(i,k) - ad(i,k+1) = 1.-al(i,k) -! exch_hx(i,k+1) = xkzh(i,k) - enddo - enddo + character(len=*), intent(out) :: errmsg + integer, intent(out) :: errflg +!local + integer :: i,j,k + +!temporary allocation of local chemical species and/or passive tracers that are vertically- +!mixed in subroutine bl_shinhong_run: + + integer, parameter :: nmix = 0 + integer :: n + + real(kind=kind_phys), dimension(ims:ime,kms:kme,jms:jme,nmix):: qmix + real(kind=kind_phys), dimension(ims:ime,kms:kme,jms:jme,nmix):: rqmixblten + + ! Local tile-sized arrays for contiguous data for bl_shinhong_run call. + + real(kind=kind_phys), dimension(its:ite,kts:kte,nmix) :: & + qmix_hv , & + rqmixblten_hv + + real(kind=kind_phys), dimension(its:ite,kts:kte) :: & + u3d_hv , & + v3d_hv , & + t3d_hv , & + qv3d_hv , & + qc3d_hv , & + qi3d_hv , & + p3d_hv , & + pi3d_hv , & + rublten_hv , & + rvblten_hv , & + rthblten_hv , & + rqvblten_hv , & + rqcblten_hv , & + rqiblten_hv , & + dz8w_hv , & + exch_h_hv , & + exch_m_hv , & + rthraten_hv + real(kind=kind_phys), dimension(its:ite,kts:kte) :: & + tke_hv , & + el_hv + + real(kind=kind_phys), dimension(its:ite,kts:kte) :: & + a_u_hv , & + a_v_hv , & + a_t_hv , & + a_e_hv , & + b_u_hv , & + a_q_hv , & + b_q_hv , & + b_v_hv , & + b_t_hv , & + b_e_hv , & + dlg_hv , & + dl_u_hv , & + vlk_hv , & + sfk_hv + real(kind=kind_phys), dimension(its:ite,kts:kte+1) :: & + p3di_hv + + real(kind=kind_phys), dimension(its:ite) :: & + psfc_hv , & + znt_hv , & + ust_hv , & + hpbl_hv , & + psim_hv , & + psih_hv , & + xland_hv , & + dx_hv , & + corf_hv , & + hfx_hv , & + qfx_hv , & + wspd_hv , & + br_hv , & + wstar_hv , & + delta_hv , & + u10_hv , & + v10_hv , & + uoce_hv , & + voce_hv , & + ctopo_hv , & + ctopo2_hv + + integer, dimension(its:ite) :: & + kpbl2d_hv + real(kind=kind_phys), dimension(its:ite) :: & + frcurb_hv + +!----------------------------------------------------------------------------------------------------------------- + + + do j = jts,jte ! - if(ndiff.ge.2) then - do ic = 2,ndiff - is = (ic-1) * kte - do k = kts,kte-1 - do i = its,ite - f3(i,k+1,ic) = qx(i,k+1+is) + ! Assign input data to local tile-sized arrays. + + do n = 1, nmix + do k = kts, kte + do i = its, ite + qmix_hv(i,k,n) = qmix(i,k,j,n) + end do + end do + end do + + do k = kts, kte+1 + do i = its, ite + p3di_hv(i,k) = p3di(i,k,j) + end do + end do + + do k = kts, kte + do i = its, ite + u3d_hv(i,k) = u3d(i,k,j) + v3d_hv(i,k) = v3d(i,k,j) + t3d_hv(i,k) = t3d(i,k,j) + qv3d_hv(i,k) = qv3d(i,k,j) + qc3d_hv(i,k) = qc3d(i,k,j) + qi3d_hv(i,k) = qi3d(i,k,j) + p3d_hv(i,k) = p3d(i,k,j) + pi3d_hv(i,k) = pi3d(i,k,j) + dz8w_hv(i,k) = dz8w(i,k,j) + rthraten_hv(i,k) = rthraten(i,k,j) + end do + end do + + if(present(a_u_bep) .and. present(a_v_bep) .and. present(a_t_bep) .and. & + present(a_q_bep) .and. present(a_e_bep) .and. present(b_u_bep) .and. & + present(b_v_bep) .and. present(b_t_bep) .and. present(b_q_bep) .and. & + present(b_e_bep) .and. present(dlg_bep) .and. present(dl_u_bep) .and. & + present(sf_bep) .and. present(vl_bep) .and. present(frc_urb2d)) then + do k = kts, kte + do i = its,ite + a_u_hv(i,k) = a_u_bep(i,k,j) + a_v_hv(i,k) = a_v_bep(i,k,j) + a_t_hv(i,k) = a_t_bep(i,k,j) + a_q_hv(i,k) = a_q_bep(i,k,j) + a_e_hv(i,k) = a_e_bep(i,k,j) + b_u_hv(i,k) = b_u_bep(i,k,j) + b_v_hv(i,k) = b_v_bep(i,k,j) + b_t_hv(i,k) = b_t_bep(i,k,j) + b_q_hv(i,k) = b_q_bep(i,k,j) + b_e_hv(i,k) = b_e_bep(i,k,j) + dlg_hv(i,k) = dlg_bep(i,k,j) + dl_u_hv(i,k) = dl_u_bep(i,k,j) + vlk_hv(i,k) = vl_bep(i,k,j) + sfk_hv(i,k) = sf_bep(i,k,j) + enddo enddo - enddo - enddo - endif -! -! copies here to avoid duplicate input args for tridin -! - do k = kts,kte - do i = its,ite - cu(i,k) = au(i,k) - enddo - enddo -! - do ic = 1,ndiff - do k = kts,kte - do i = its,ite - r3(i,k,ic) = f3(i,k,ic) - enddo - enddo - enddo -! -! solve tridiagonal problem for moisture, clouds, and gases -! - call tridin_ysu(al,ad,cu,r3,au,f3,its,ite,kts,kte,ndiff) -! -! recover tendencies of heat and moisture -! - do k = kte,kts,-1 - do i = its,ite - qtend = (f3(i,k,1)-qx(i,k))*rdt - qtnp(i,k) = qtnp(i,k)+qtend - dqsfc(i) = dqsfc(i)+qtend*conq*del(i,k) - if(k.eq.kte) then - qflux_e(i,k) = qtend*dz8w2d(i,k) - else - qflux_e(i,k) = qflux_e(i,k+1) + qtend*dz8w2d(i,k) - endif - tvflux_e(i,k) = tflux_e(i,k) + qflux_e(i,k)*ep1*thx(i,k) - enddo - enddo -! - do k = kts,kte - do i = its,ite - if(pblflg(i).and.k.lt.kpbl(i)) then - hgame_c=c_1*0.2*2.5*(g/thvx(i,k))*wstar(i)/(0.25*(q2x(i,k+1)+q2x(i,k))) - hgame_c=min(hgame_c,gamcre) - if(k.eq.kte)then - hgame2d(i,k)=hgame_c*0.5*tvflux_e(i,k)*hpbl(i) - hgame2d(i,k)=max(hgame2d(i,k),0.0) - else - hgame2d(i,k)=hgame_c*0.5*(tvflux_e(i,k)+tvflux_e(i,k+1))*hpbl(i) - hgame2d(i,k)=max(hgame2d(i,k),0.0) - endif - endif - enddo - enddo -! - if(ndiff.ge.2) then - do ic = 2,ndiff - if(ifvmix(ic)) then - is = (ic-1) * kte - do k = kte,kts,-1 - do i = its,ite - qtend = (f3(i,k,ic)-qx(i,k+is))*rdt - qtnp(i,k+is) = qtnp(i,k+is)+qtend - enddo + do i = its, ite + frcurb_hv(i) = frc_urb2d(i,j) enddo - endif - enddo - endif -! -! compute tridiagonal matrix elements for momentum -! - do i = its,ite - do k = kts,kte - au(i,k) = 0. - al(i,k) = 0. - ad(i,k) = 0. - f1(i,k) = 0. - f2(i,k) = 0. - enddo - enddo -! - do i = its,ite -! paj: ctopo=1 if topo_wind=0 (default) -! mchen add this line to make sure NMM can still work with YSU PBL - ad(i,1) = 1.+ctopo(i)*ust(i)**2/wspd1(i)*rhox(i)*g/del(i,1)*dt2 & - *(wspd1(i)/wspd(i))**2 - f1(i,1) = ux(i,1) - f2(i,1) = vx(i,1) - enddo -! - delxy=sqrt(dx*dy) - do k = kts,kte-1 - do i = its,ite - dtodsd = dt2/del(i,k) - dtodsu = dt2/del(i,k+1) - dsig = p2d(i,k)-p2d(i,k+1) - rdz = 1./dza(i,k+1) - tem1 = dsig*xkzm(i,k)*rdz - if(pblflg(i).and.k.lt.kpbl(i))then - dsdzu = tem1*(-hgamu(i)/hpbl(i)-ufxpbl(i)*zfacent(i,k)/xkzm(i,k)) - dsdzv = tem1*(-hgamv(i)/hpbl(i)-vfxpbl(i)*zfacent(i,k)/xkzm(i,k)) - f1(i,k) = f1(i,k)+dtodsd*dsdzu - f1(i,k+1) = ux(i,k+1)-dtodsu*dsdzu - f2(i,k) = f2(i,k)+dtodsd*dsdzv - f2(i,k+1) = vx(i,k+1)-dtodsu*dsdzv - elseif(pblflg(i).and.k.ge.kpbl(i).and.entfac(i,k).lt.4.6) then - xkzm(i,k) = prpbl(i)*xkzh(i,k) - xkzm(i,k) = sqrt(xkzm(i,k)*xkzml(i,k)) - xkzm(i,k) = max(xkzm(i,k),xkzom(i,k)) - xkzm(i,k) = min(xkzm(i,k),xkzmax) - f1(i,k+1) = ux(i,k+1) - f2(i,k+1) = vx(i,k+1) - else - f1(i,k+1) = ux(i,k+1) - f2(i,k+1) = vx(i,k+1) - endif - tem1 = dsig*xkzm(i,k)*rdz - dsdz2 = tem1*rdz - au(i,k) = -dtodsd*dsdz2 - al(i,k) = -dtodsu*dsdz2 -! -! scale dependency for local momentum transport -! - zfacdx=0.2*hpbl(i)/zq(i,k+1) - delxy=sqrt(dx*dy)*max(zfacdx,1.0) - pu1=pu(delxy,hpbl(i)) - if(pblflg(i).and.k.lt.kpbl(i)) then - au(i,k) = au(i,k)*pu1 - al(i,k) = al(i,k)*pu1 - endif - ad(i,k) = ad(i,k)-au(i,k) - ad(i,k+1) = 1.-al(i,k) - enddo - enddo -! -! copies here to avoid duplicate input args for tridin -! - do k = kts,kte - do i = its,ite - cu(i,k) = au(i,k) - r1(i,k) = f1(i,k) - r2(i,k) = f2(i,k) - enddo - enddo -! -! solve tridiagonal problem for momentum -! - call tridi1n(al,ad,cu,r1,r2,au,f1,f2,its,ite,kts,kte,1) -! -! recover tendencies of momentum -! - do k = kte,kts,-1 - do i = its,ite - utend = (f1(i,k)-ux(i,k))*rdt - vtend = (f2(i,k)-vx(i,k))*rdt - utnp(i,k) = utnp(i,k)+utend - vtnp(i,k) = vtnp(i,k)+vtend - dusfc(i) = dusfc(i) + utend*conwrc*del(i,k) - dvsfc(i) = dvsfc(i) + vtend*conwrc*del(i,k) - enddo - enddo -! - do i = its,ite - kpbl1d(i) = kpbl(i) - enddo -! -! paj: ctopo2=1 if topo_wind=0 (default) -! - do i = its,ite - u10(i) = ctopo2(i)*u10(i)+(1-ctopo2(i))*ux(i,1) - v10(i) = ctopo2(i)*v10(i)+(1-ctopo2(i))*vx(i,1) - enddo -! -!---- calculate sgs tke which is consistent with shinhongpbl algorithm -! - if (shinhong_tke_diag.eq.1) then -! - tke_calculation: do i = its,ite - do k = kts+1,kte - s2(k) = 0.0 - gh(k) = 0.0 - rig(k) = 0.0 - el(k) = 0.0 - akmk(k) = 0.0 - akhk(k) = 0.0 - mfk(k) = 0.0 - ufxpblk(k) = 0.0 - vfxpblk(k) = 0.0 - qfxpblk(k) = 0.0 - enddo -! - do k = kts,kte - uxk(k) = 0.0 - vxk(k) = 0.0 - txk(k) = 0.0 - thxk(k) = 0.0 - thvxk(k) = 0.0 - q2xk(k) = 0.0 - hgame(k) = 0.0 - ps1d(k) = 0.0 - pb1d(k) = 0.0 - eps1d(k) = 0.0 - pt1d(k) = 0.0 - xkze1d(k) = 0.0 - eflx_l1d(k) = 0.0 - eflx_nl1d(k) = 0.0 - ptke1(k) = 1.0 - enddo -! - do k = kts,kte+1 - zqk(k) = 0.0 - enddo -! - do k = kts,kte*ndiff - qxk(k) = 0.0 - enddo -! - do k = kts,kte - uxk(k) = ux(i,k) - vxk(k) = vx(i,k) - txk(k) = tx(i,k) - thxk(k) = thx(i,k) - thvxk(k) = thvx(i,k) - q2xk(k) = q2x(i,k) - hgame(k) = hgame2d(i,k) - enddo -! - do k = kts,kte-1 - if(pblflg(i).and.k.le.kpbl(i)) then - zfacdx = 0.2*hpbl(i)/za(i,k) - delxy = sqrt(dx*dy)*max(zfacdx,1.0) - ptke1(k+1) = ptke(delxy,hpbl(i)) - endif - enddo -! - do k = kts,kte+1 - zqk(k) = zq(i,k) - enddo -! - do k = kts,kte*ndiff - qxk(k) = qx(i,k) - enddo -! - do k = kts+1,kte - akmk(k) = xkzm(i,k-1) - akhk(k) = xkzh(i,k-1) - mfk(k) = mf(i,k-1)/xkzh(i,k-1) - ufxpblk(k) = ufxpbl(i)*zfacent(i,k-1)/xkzm(i,k-1) - vfxpblk(k) = vfxpbl(i)*zfacent(i,k-1)/xkzm(i,k-1) - qfxpblk(k) = qfxpbl(i)*zfacent(i,k-1)/xkzq(i,k-1) - enddo -! - if(pblflg(i)) then - k = kpbl(i) - 1 - dex = 0.25*(q2xk(k+2)-q2xk(k)) - efxpbl(i) = we(i)*dex - endif -! - delxy=sqrt(dx*dy) -! -!---- find the mixing length -! - call mixlen(lmh,uxk,vxk,txk,thxk,qxk(kts),qxk(kte+1) & - ,q2xk,zqk,ust(i),corf(i),epshol(i) & - ,s2,gh,rig,el & - ,hpbl(i),kpbl(i),lmxl,ct(i) & - ,hgamu(i),hgamv(i),hgamq(i),pblflg(i) & - ,mfk,ufxpblk,vfxpblk,qfxpblk & - ,ep1,karman,cp & - ,ids,ide,jds,jde,kds,kde & - ,ims,ime,jms,jme,kms,kme & - ,its,ite,jts,jte,kts,kte ) -! -!---- solve for the production/dissipation of the turbulent kinetic energy -! - call prodq2(lmh,dt,ust(i),s2,rig,q2xk,el,zqk,akmk,akhk & - ,uxk,vxk,thxk,thvxk & - ,hgamu(i),hgamv(i),hgamq(i),delxy & - ,hpbl(i),pblflg(i),kpbl(i) & - ,mfk,ufxpblk,vfxpblk,qfxpblk & - ,ep1 & - ,ids,ide,jds,jde,kds,kde & - ,ims,ime,jms,jme,kms,kme & - ,its,ite,jts,jte,kts,kte ) -! -! -!---- carry out the vertical diffusion of turbulent kinetic energy -! - call vdifq(lmh,dt,q2xk,el,zqk & - ,akhk,ptke1 & - ,hgame,hpbl(i),pblflg(i),kpbl(i) & - ,efxpbl(i) & - ,ids,ide,jds,jde,kds,kde & - ,ims,ime,jms,jme,kms,kme & - ,its,ite,jts,jte,kts,kte ) -! -!---- save the new tke and mixing length. -! - do k = kts,kte - q2x(i,k) = amax1(q2xk(k),epsq2l) - tke(i,k) = 0.5*q2x(i,k) - if(k/=kts) el_pbl(i,k) = el(k) ! el is not defined at kte - enddo -! - enddo tke_calculation - endif -! -!---- end of tke calculation -! -! -!---- end of vertical diffusion -! - end subroutine shinhong2d -!------------------------------------------------------------------------------- -! -!------------------------------------------------------------------------------- - subroutine tridi1n(cl,cm,cu,r1,r2,au,f1,f2,its,ite,kts,kte,nt) -!------------------------------------------------------------------------------- - implicit none -!------------------------------------------------------------------------------- -! - integer, intent(in ) :: its,ite, kts,kte, nt -! - real, dimension( its:ite, kts+1:kte+1 ) , & - intent(in ) :: cl -! - real, dimension( its:ite, kts:kte ) , & - intent(in ) :: cm, & - r1 - real, dimension( its:ite, kts:kte,nt ) , & - intent(in ) :: r2 -! - real, dimension( its:ite, kts:kte ) , & - intent(inout) :: au, & - cu, & - f1 - real, dimension( its:ite, kts:kte,nt ) , & - intent(inout) :: f2 -! - real :: fk - integer :: i,k,l,n,it -! -!------------------------------------------------------------------------------- -! - l = ite - n = kte -! - do i = its,l - fk = 1./cm(i,1) - au(i,1) = fk*cu(i,1) - f1(i,1) = fk*r1(i,1) - enddo -! - do it = 1,nt - do i = its,l - fk = 1./cm(i,1) - f2(i,1,it) = fk*r2(i,1,it) - enddo - enddo -! - do k = kts+1,n-1 - do i = its,l - fk = 1./(cm(i,k)-cl(i,k)*au(i,k-1)) - au(i,k) = fk*cu(i,k) - f1(i,k) = fk*(r1(i,k)-cl(i,k)*f1(i,k-1)) - enddo - enddo -! - do it = 1,nt - do k = kts+1,n-1 - do i = its,l - fk = 1./(cm(i,k)-cl(i,k)*au(i,k-1)) - f2(i,k,it) = fk*(r2(i,k,it)-cl(i,k)*f2(i,k-1,it)) - enddo - enddo - enddo -! - do i = its,l - fk = 1./(cm(i,n)-cl(i,n)*au(i,n-1)) - f1(i,n) = fk*(r1(i,n)-cl(i,n)*f1(i,n-1)) - enddo -! - do it = 1,nt - do i = its,l - fk = 1./(cm(i,n)-cl(i,n)*au(i,n-1)) - f2(i,n,it) = fk*(r2(i,n,it)-cl(i,n)*f2(i,n-1,it)) - enddo - enddo -! - do k = n-1,kts,-1 - do i = its,l - f1(i,k) = f1(i,k)-au(i,k)*f1(i,k+1) - enddo - enddo -! - do it = 1,nt - do k = n-1,kts,-1 - do i = its,l - f2(i,k,it) = f2(i,k,it)-au(i,k)*f2(i,k+1,it) - enddo - enddo - enddo -! - end subroutine tridi1n -!------------------------------------------------------------------------------- -! -!------------------------------------------------------------------------------- - subroutine tridin_ysu(cl,cm,cu,r2,au,f2,its,ite,kts,kte,nt) -!------------------------------------------------------------------------------- - implicit none -!------------------------------------------------------------------------------- -! - integer, intent(in ) :: its,ite, kts,kte, nt -! - real, dimension( its:ite, kts+1:kte+1 ) , & - intent(in ) :: cl -! - real, dimension( its:ite, kts:kte ) , & - intent(in ) :: cm - real, dimension( its:ite, kts:kte,nt ) , & - intent(in ) :: r2 -! - real, dimension( its:ite, kts:kte ) , & - intent(inout) :: au, & - cu - real, dimension( its:ite, kts:kte,nt ) , & - intent(inout) :: f2 -! - real :: fk - integer :: i,k,l,n,it -! -!------------------------------------------------------------------------------- -! - l = ite - n = kte -! - do it = 1,nt - do i = its,l - fk = 1./cm(i,1) - au(i,1) = fk*cu(i,1) - f2(i,1,it) = fk*r2(i,1,it) - enddo - enddo -! - do it = 1,nt - do k = kts+1,n-1 - do i = its,l - fk = 1./(cm(i,k)-cl(i,k)*au(i,k-1)) - au(i,k) = fk*cu(i,k) - f2(i,k,it) = fk*(r2(i,k,it)-cl(i,k)*f2(i,k-1,it)) - enddo - enddo - enddo -! - do it = 1,nt - do i = its,l - fk = 1./(cm(i,n)-cl(i,n)*au(i,n-1)) - f2(i,n,it) = fk*(r2(i,n,it)-cl(i,n)*f2(i,n-1,it)) - enddo - enddo -! - do it = 1,nt - do k = n-1,kts,-1 - do i = its,l - f2(i,k,it) = f2(i,k,it)-au(i,k)*f2(i,k+1,it) - enddo - enddo - enddo -! - end subroutine tridin_ysu -!------------------------------------------------------------------------------- -! -!------------------------------------------------------------------------------- - subroutine mixlen(lmh,u,v,t,the,q,cwm,q2,z,ustar,corf,epshol, & - s2,gh,ri,el,hpbl,lpbl,lmxl,ct, & - hgamu,hgamv,hgamq,pblflg, & - mf,ufxpbl,vfxpbl,qfxpbl, & - p608,vkarman,cp, & - ids,ide,jds,jde,kds,kde, & - ims,ime,jms,jme,kms,kme, & - its,ite,jts,jte,kts,kte) -!------------------------------------------------------------------------------- - implicit none -!------------------------------------------------------------------------------- -! qnse model constants -!------------------------------------------------------------------------------- - real,parameter :: blckdr=0.0063,cn=0.75 - real,parameter :: eps1=1.e-12,epsl=0.32,epsru=1.e-7,epsrs=1.e-7 - real,parameter :: el0max=1000.,el0min=1.,elfc=0.23*0.5 - real,parameter :: alph=0.30,beta=1./273.,g=9.81,btg=beta*g - real,parameter :: a1=0.659888514560862645,a2x=0.6574209922667784586 - real,parameter :: b1=11.87799326209552761,b2=7.226971804046074028 - real,parameter :: c1=0.000830955950095854396 - real,parameter :: adnh= 9.*a1*a2x*a2x*(12.*a1+3.*b2)*btg*btg - real,parameter :: adnm=18.*a1*a1*a2x*(b2-3.*a2x)*btg - real,parameter :: bdnh= 3.*a2x*(7.*a1+b2)*btg,bdnm= 6.*a1*a1 -!------------------------------------------------------------------------------- -! free term in the equilibrium equation for (l/q)**2 -!------------------------------------------------------------------------------- - real,parameter :: aeqh=9.*a1*a2x*a2x*b1*btg*btg & - +9.*a1*a2x*a2x*(12.*a1+3.*b2)*btg*btg - real,parameter :: aeqm=3.*a1*a2x*b1*(3.*a2x+3.*b2*c1+18.*a1*c1-b2) & - *btg+18.*a1*a1*a2x*(b2-3.*a2x)*btg -!------------------------------------------------------------------------------- -! forbidden turbulence area -!------------------------------------------------------------------------------- - real,parameter :: requ=-aeqh/aeqm - real,parameter :: epsgh=1.e-9,epsgm=requ*epsgh -!------------------------------------------------------------------------------- -! near isotropy for shear turbulence, ww/q2 lower limit -!------------------------------------------------------------------------------- - real,parameter :: ubryl=(18.*requ*a1*a1*a2x*b2*c1*btg & - +9.*a1*a2x*a2x*b2*btg*btg) & - /(requ*adnm+adnh) - real,parameter :: ubry=(1.+epsrs)*ubryl,ubry3=3.*ubry - real,parameter :: aubh=27.*a1*a2x*a2x*b2*btg*btg-adnh*ubry3 - real,parameter :: aubm=54.*a1*a1*a2x*b2*c1*btg -adnm*ubry3 - real,parameter :: bubh=(9.*a1*a2x+3.*a2x*b2)*btg-bdnh*ubry3 - real,parameter :: bubm=18.*a1*a1*c1 -bdnm*ubry3 - real,parameter :: cubr=1.-ubry3,rcubr=1./cubr -!------------------------------------------------------------------------------- -! k profile constants -!------------------------------------------------------------------------------- - real,parameter :: elcbl=0.77 -!------------------------------------------------------------------------------- -! - integer, intent(in ) :: ids,ide, jds,jde, kds,kde, & - ims,ime, jms,jme, kms,kme, & - its,ite, jts,jte, kts,kte - integer, intent(in ) :: lmh,lmxl,lpbl -! - real, intent(in ) :: p608,vkarman,cp - real, intent(in ) :: hpbl,corf,ustar,hgamu,hgamv,hgamq - real, intent(inout) :: ct,epshol -! - real, dimension( kts:kte ) , & - intent(in ) :: cwm, & - q, & - q2, & - t, & - the, & - u, & - v -! - real, dimension( kts+1:kte ) , & - intent(in ) :: mf, & - ufxpbl, & - vfxpbl, & - qfxpbl -! - real, dimension( kts:kte+1 ) , & - intent(in ) :: z -! - real, dimension( kts+1:kte ) , & - intent(out ) :: el, & - ri, & - gh, & - s2 -! - logical,intent(in) :: pblflg -! -! local vars -! - integer :: k,lpblm - real :: suk,svk,elocp - real :: a,aden,b,bden,aubr,bubr,blmx,el0,eloq2x,ghl,s2l, & - qol2st,qol2un,qdzl,rdz,sq,srel,szq,tem,thm,vkrmz,rlambda, & - rlb,rln,f - real :: ckp - real, dimension( kts:kte ) :: q1, & - en2 - real, dimension( kts+1:kte ) :: dth, & - elm, & - rel -! -!------------------------------------------------------------------------------- -! - elocp=2.72e6/cp - ct=0. -! - do k = kts,kte - q1(k) = 0. - enddo -! - do k = kts+1,kte - dth(k) = the(k)-the(k-1) - enddo -! - do k = kts+2,kte - if(dth(k)>0..and.dth(k-1)<=0.)then - dth(k)=dth(k)+ct - exit - endif - enddo -! -! compute local gradient richardson number -! - do k = kte,kts+1,-1 - rdz=2./(z(k+1)-z(k-1)) - s2l=((u(k)-u(k-1))**2+(v(k)-v(k-1))**2)*rdz*rdz ! s**2 - if(pblflg.and.k.le.lpbl)then - suk=(u(k)-u(k-1))*rdz - svk=(v(k)-v(k-1))*rdz - s2l=(suk-hgamu/hpbl-ufxpbl(k))*suk+(svk-hgamv/hpbl-vfxpbl(k))*svk - endif - s2l=max(s2l,epsgm) - s2(k)=s2l -! - tem=(t(k)+t(k-1))*0.5 - thm=(the(k)+the(k-1))*0.5 - a=thm*p608 - b=(elocp/tem-1.-p608)*thm - ghl=(dth(k)*((q(k)+q(k-1)+cwm(k)+cwm(k-1))*(0.5*p608)+1.) & - +(q(k)-q(k-1)+cwm(k)-cwm(k-1))*a & - +(cwm(k)-cwm(k-1))*b)*rdz ! dtheta/dz - if(pblflg.and.k.le.lpbl)then - ghl=ghl-mf(k)-(hgamq/hpbl+qfxpbl(k))*a - endif - if(abs(ghl)<=epsgh)ghl=epsgh -! - en2(k)=ghl*g/thm ! n**2 - gh(k)=ghl - ri(k)=en2(k)/s2l - enddo -! -! find maximum mixing lengths and the level of the pbl top -! - do k = kte,kts+1,-1 - s2l=s2(k) - ghl=gh(k) - if(ghl>=epsgh)then - if(s2l/ghl<=requ)then - elm(k)=epsl - else - aubr=(aubm*s2l+aubh*ghl)*ghl - bubr= bubm*s2l+bubh*ghl - qol2st=(-0.5*bubr+sqrt(bubr*bubr*0.25-aubr*cubr))*rcubr - eloq2x=1./qol2st - elm(k)=max(sqrt(eloq2x*q2(k)),epsl) - endif - else - aden=(adnm*s2l+adnh*ghl)*ghl - bden= bdnm*s2l+bdnh*ghl - qol2un=-0.5*bden+sqrt(bden*bden*0.25-aden) - eloq2x=1./(qol2un+epsru) ! repsr1/qol2un - elm(k)=max(sqrt(eloq2x*q2(k)),epsl) - endif - enddo -! - do k = lpbl,lmh,-1 - q1(k)=sqrt(q2(k)) - enddo -! - szq=0. - sq =0. - do k = kte,kts+1,-1 - qdzl=(q1(k)+q1(k-1))*(z(k)-z(k-1)) - szq=(z(k)+z(k-1)-z(lmh)-z(lmh))*qdzl+szq - sq=qdzl+sq - enddo -! -! computation of asymptotic l in blackadar formula -! - el0=min(alph*szq*0.5/sq,el0max) - el0=max(el0 ,el0min) -! -! above the pbl top -! - lpblm=min(lpbl+1,kte) - do k = kte,lpblm,-1 - el(k)=(z(k+1)-z(k-1))*elfc - rel(k)=el(k)/elm(k) - enddo -! -! inside the pbl -! - epshol=min(epshol,0.0) - ckp=elcbl*((1.0-8.0*epshol)**(1./3.)) - if(lpbl>lmh)then - do k = lpbl,lmh+1,-1 - vkrmz=(z(k)-z(lmh))*vkarman - if(pblflg) then - vkrmz=ckp*(z(k)-z(lmh))*vkarman - el(k)=vkrmz/(vkrmz/el0+1.) - else - el(k)=vkrmz/(vkrmz/el0+1.) - endif - rel(k)=el(k)/elm(k) - enddo - endif -! - do k = lpbl-1,lmh+2,-1 - srel=min(((rel(k-1)+rel(k+1))*0.5+rel(k))*0.5,rel(k)) - el(k)=max(srel*elm(k),epsl) - enddo -! -! mixing length for the qnse model in stable case -! - f=max(corf,eps1) - rlambda=f/(blckdr*ustar) - do k = kte,kts+1,-1 - if(en2(k)>=0.0)then ! stable case - vkrmz=(z(k)-z(lmh))*vkarman - rlb=rlambda+1./vkrmz - rln=sqrt(2.*en2(k)/q2(k))/cn - el(k)=1./(rlb+rln) - endif - enddo -! - end subroutine mixlen -!------------------------------------------------------------------------------- -! -!------------------------------------------------------------------------------- - subroutine prodq2(lmh,dtturbl,ustar,s2,ri,q2,el,z,akm,akh, & - uxk,vxk,thxk,thvxk, & - hgamu,hgamv,hgamq,delxy, & - hpbl,pblflg,kpbl, & - mf,ufxpbl,vfxpbl,qfxpbl, & - p608, & - ids,ide,jds,jde,kds,kde, & - ims,ime,jms,jme,kms,kme, & - its,ite,jts,jte,kts,kte) -!------------------------------------------------------------------------------- - implicit none -!------------------------------------------------------------------------------- -! - real,parameter :: epsq2l = 0.01,c0 = 0.55,ceps = 16.6,g = 9.81 -! - integer, intent(in ) :: ids,ide, jds,jde, kds,kde, & - ims,ime, jms,jme, kms,kme, & - its,ite, jts,jte, kts,kte - integer, intent(in ) :: lmh,kpbl -! - real, intent(in ) :: p608,dtturbl,ustar - real, intent(in ) :: hgamu,hgamv,hgamq,delxy,hpbl -! - logical, intent(in ) :: pblflg -! - real, dimension( kts:kte ) , & - intent(in ) :: uxk, & - vxk, & - thxk, & - thvxk - real, dimension( kts+1:kte ) , & - intent(in ) :: s2, & - ri, & - akm, & - akh, & - el, & - mf, & - ufxpbl, & - vfxpbl, & - qfxpbl -! - real, dimension( kts:kte+1 ) , & - intent(in ) :: z -! - real, dimension( kts:kte ) , & - intent(inout) :: q2 -! -! local vars -! - integer :: k -! - real :: s2l,q2l,deltaz,akml,akhl,en2,pr,bpr,dis,rc02 - real :: suk,svk,gthvk,govrthvk,pru,prv - real :: thm,disel -! -!------------------------------------------------------------------------------- -! - rc02=2.0/(c0*c0) -! -! start of production/dissipation loop -! - main_integration: do k = kts+1,kte - deltaz=0.5*(z(k+1)-z(k-1)) - s2l=s2(k) - q2l=q2(k) - suk=(uxk(k)-uxk(k-1))/deltaz - svk=(vxk(k)-vxk(k-1))/deltaz - gthvk=(thvxk(k)-thvxk(k-1))/deltaz - govrthvk=g/(0.5*(thvxk(k)+thvxk(k-1))) - akml=akm(k) - akhl=akh(k) - en2=ri(k)*s2l !n**2 - thm=(thxk(k)+thxk(k-1))*0.5 -! -! turbulence production term -! - if(pblflg.and.k.le.kpbl)then - pru=(akml*(suk-hgamu/hpbl-ufxpbl(k)))*suk - prv=(akml*(svk-hgamv/hpbl-vfxpbl(k)))*svk - else - pru=akml*suk*suk - prv=akml*svk*svk - endif - pr=pru+prv -! -! buoyancy production -! - if(pblflg.and.k.le.kpbl)then - bpr=(akhl*(gthvk-mf(k)-(hgamq/hpbl+qfxpbl(k))*p608*thm))*govrthvk - else - bpr=akhl*gthvk*govrthvk - endif -! -! dissipation -! - disel=min(delxy,ceps*el(k)) - dis=(q2l)**1.5/disel -! - q2l=q2l+2.0*(pr-bpr-dis)*dtturbl - q2(k)=amax1(q2l,epsq2l) -! -! end of production/dissipation loop -! - enddo main_integration -! -! lower boundary condition for q2 -! - q2(kts)=amax1(rc02*ustar*ustar,epsq2l) -! - end subroutine prodq2 -!------------------------------------------------------------------------------- -! -!------------------------------------------------------------------------------- - subroutine vdifq(lmh,dtdif,q2,el,z, & - akhk,ptke1, & - hgame,hpbl,pblflg,kpbl, & - efxpbl, & - ids,ide,jds,jde,kds,kde, & - ims,ime,jms,jme,kms,kme, & - its,ite,jts,jte,kts,kte) -!------------------------------------------------------------------------------- - implicit none -!------------------------------------------------------------------------------- -! - real,parameter :: c_k=1.0,esq=5.0 -! - integer, intent(in ) :: ids,ide, jds,jde, kds,kde, & - ims,ime, jms,jme, kms,kme, & - its,ite, jts,jte, kts,kte - integer, intent(in ) :: lmh,kpbl -! - real, intent(in ) :: dtdif,hpbl,efxpbl -! - logical, intent(in ) :: pblflg -! - real, dimension( kts:kte ) , & - intent(in ) :: hgame, & - ptke1 - real, dimension( kts+1:kte ) , & - intent(in ) :: el, & - akhk - real, dimension( kts:kte+1 ) , & - intent(in ) :: z -! - real, dimension( kts:kte ) , & - intent(inout) :: q2 -! -! local vars -! - integer :: k -! - real :: aden,akqs,bden,besh,besm,cden,cf,dtozs,ell,eloq2,eloq4 - real :: elqdz,esh,esm,esqhf,ghl,gml,q1l,rden,rdz - real :: zak -! - real, dimension( kts+1:kte ) :: zfacentk - real, dimension( kts+2:kte ) :: akq, & - cm, & - cr, & - dtoz, & - rsq2 -! -!------------------------------------------------------------------------------- -! -! vertical turbulent diffusion -! - esqhf=0.5*esq - do k = kts+1,kte - zak=0.5*(z(k)+z(k-1)) !zak of vdifq = za(k-1) of shinhong2d - zfacentk(k)=(zak/hpbl)**3.0 - enddo -! - do k = kte,kts+2,-1 - dtoz(k)=(dtdif+dtdif)/(z(k+1)-z(k-1)) - akq(k)=c_k*(akhk(k)/(z(k+1)-z(k-1))+akhk(k-1)/(z(k)-z(k-2))) - akq(k)=akq(k)*ptke1(k) - cr(k)=-dtoz(k)*akq(k) - enddo -! - akqs=c_k*akhk(kts+1)/(z(kts+2)-z(kts)) - akqs=akqs*ptke1(kts+1) - cm(kte)=dtoz(kte)*akq(kte)+1. - rsq2(kte)=q2(kte) -! - do k = kte-1,kts+2,-1 - cf=-dtoz(k)*akq(k+1)/cm(k+1) - cm(k)=-cr(k+1)*cf+(akq(k+1)+akq(k))*dtoz(k)+1. - rsq2(k)=-rsq2(k+1)*cf+q2(k) - if(pblflg.and.k.lt.kpbl) then - rsq2(k)=rsq2(k)-dtoz(k)*(2.0*hgame(k)/hpbl)*akq(k+1)*(z(k+1)-z(k)) & - +dtoz(k)*(2.0*hgame(k-1)/hpbl)*akq(k)*(z(k)-z(k-1)) - rsq2(k)=rsq2(k)-dtoz(k)*2.0*efxpbl*zfacentk(k+1) & - +dtoz(k)*2.0*efxpbl*zfacentk(k) - endif - enddo -! - dtozs=(dtdif+dtdif)/(z(kts+2)-z(kts)) - cf=-dtozs*akq(lmh+2)/cm(lmh+2) -! - if(pblflg.and.((lmh+1).lt.kpbl)) then - q2(lmh+1)=(dtozs*akqs*q2(lmh)-rsq2(lmh+2)*cf+q2(lmh+1) & - -dtozs*(2.0*hgame(lmh+1)/hpbl)*akq(lmh+2)*(z(lmh+2)-z(lmh+1)) & - +dtozs*(2.0*hgame(lmh)/hpbl)*akqs*(z(lmh+1)-z(lmh))) - q2(lmh+1)=q2(lmh+1)-dtozs*2.0*efxpbl*zfacentk(lmh+2) & - +dtozs*2.0*efxpbl*zfacentk(lmh+1) - q2(lmh+1)=q2(lmh+1)/((akq(lmh+2)+akqs)*dtozs-cr(lmh+2)*cf+1.) - else - q2(lmh+1)=(dtozs*akqs*q2(lmh)-rsq2(lmh+2)*cf+q2(lmh+1)) & - /((akq(lmh+2)+akqs)*dtozs-cr(lmh+2)*cf+1.) - endif -! - do k = lmh+2,kte - q2(k)=(-cr(k)*q2(k-1)+rsq2(k))/cm(k) - enddo -! - end subroutine vdifq -!------------------------------------------------------------------------------- -! -!------------------------------------------------------------------------------- - subroutine shinhonginit(rublten,rvblten,rthblten,rqvblten, & - rqcblten,rqiblten, & - tke_pbl, & - p_qi,p_first_scalar, & - restart, allowed_to_read, & - ids, ide, jds, jde, kds, kde, & - ims, ime, jms, jme, kms, kme, & - its, ite, jts, jte, kts, kte ) -!------------------------------------------------------------------------------- - implicit none -!------------------------------------------------------------------------------- -! - real,parameter :: epsq2l = 0.01 - logical , intent(in) :: restart, allowed_to_read - integer , intent(in) :: ids, ide, jds, jde, kds, kde, & - ims, ime, jms, jme, kms, kme, & - its, ite, jts, jte, kts, kte - integer , intent(in) :: p_qi,p_first_scalar - real , dimension( ims:ime , kms:kme , jms:jme ), intent(out) :: & - rublten, & - rvblten, & - rthblten, & - rqvblten, & - rqcblten, & - rqiblten - real , dimension( ims:ime , kms:kme , jms:jme ), intent(out) :: & - tke_pbl - integer :: i, j, k, itf, jtf, ktf -! - jtf = min0(jte,jde-1) - ktf = min0(kte,kde-1) - itf = min0(ite,ide-1) -! - if(.not.restart)then - do j = jts,jtf - do k = kts,ktf - do i = its,itf - rublten(i,k,j) = 0. - rvblten(i,k,j) = 0. - rthblten(i,k,j) = 0. - rqvblten(i,k,j) = 0. - rqcblten(i,k,j) = 0. - tke_pbl(i,k,j) = epsq2l/2. - enddo - enddo - enddo - endif -! - if (p_qi .ge. p_first_scalar .and. .not.restart) then - do j = jts,jtf - do k = kts,ktf - do i = its,itf - rqiblten(i,k,j) = 0. - enddo - enddo - enddo - endif -! - end subroutine shinhonginit -!------------------------------------------------------------------------------- -! -!------------------------------------------------------------------------------- - function pu(d,h) -!------------------------------------------------------------------------------- - implicit none -!------------------------------------------------------------------------------- - real :: pu - real,parameter :: pmin = 0.0,pmax = 1.0 - real,parameter :: a1 = 1.0, a2 = 0.070, a3 = 1.0, a4 = 0.142, a5 = 0.071 - real,parameter :: b1 = 2.0, b2 = 0.6666667 - real :: d,h,doh,num,den -! - if (h .ne. 0) then - doh=d/h - num=a1*(doh)**b1+a2*(doh)**b2 - den=a3*(doh)**b1+a4*(doh)**b2+a5 - pu=num/den - else - pu=1. - endif - pu=max(pu,pmin) - pu=min(pu,pmax) -! - return - end function -!------------------------------------------------------------------------------- -! -!------------------------------------------------------------------------------- - function pq(d,h) -!------------------------------------------------------------------------------- - implicit none -!------------------------------------------------------------------------------- - real :: pq - real,parameter :: pmin = 0.0,pmax = 1.0 - real,parameter :: a1 = 1.0, a2 = -0.098, a3 = 1.0, a4 = 0.106, a5 = 0.5 - real,parameter :: b1 = 2.0 - real :: d,h,doh,num,den -! - if (h .ne. 0) then - doh=d/h - num=a1*(doh)**b1+a2 - den=a3*(doh)**b1+a4 - pq=a5*num/den+(1.-a5) - else - pq=1. - endif - pq=max(pq,pmin) - pq=min(pq,pmax) -! - return - end function -!------------------------------------------------------------------------------- -! -!------------------------------------------------------------------------------- - function pthnl(d,h) -!------------------------------------------------------------------------------- - implicit none -!------------------------------------------------------------------------------- - real :: pthnl - real,parameter :: pmin = 0.0,pmax = 1.0 - real,parameter :: a1 = 1.000, a2 = 0.936, a3 = -1.110, & - a4 = 1.000, a5 = 0.312, a6 = 0.329, a7 = 0.243 - real,parameter :: b1 = 2.0, b2 = 0.875 - real :: d,h,doh,num,den -! - if (h .ne. 0) then - doh=d/h - num=a1*(doh)**b1+a2*(doh)**b2+a3 - den=a4*(doh)**b1+a5*(doh)**b2+a6 - pthnl=a7*num/den+(1.-a7) - else - pthnl=1. - endif - pthnl=max(pthnl,pmin) - pthnl=min(pthnl,pmax) -! - return - end function -!------------------------------------------------------------------------------- -! -!------------------------------------------------------------------------------- - function pthl(d,h) -!------------------------------------------------------------------------------- - implicit none -!------------------------------------------------------------------------------- - real :: pthl - real,parameter :: pmin = 0.0,pmax = 1.0 - real,parameter :: a1 = 1.000, a2 = 0.870, a3 = -0.913, & - a4 = 1.000, a5 = 0.153, a6 = 0.278, a7 = 0.280 - real,parameter :: b1 = 2.0, b2 = 0.5 - real :: d,h,doh,num,den -! - if (h .ne. 0) then - doh=d/h - num=a1*(doh)**b1+a2*(doh)**b2+a3 - den=a4*(doh)**b1+a5*(doh)**b2+a6 - pthl=a7*num/den+(1.-a7) - else - pthl=1. - endif - pthl=max(pthl,pmin) - pthl=min(pthl,pmax) -! - return - end function -!------------------------------------------------------------------------------- -! -!------------------------------------------------------------------------------- - function ptke(d,h) -!------------------------------------------------------------------------------- - implicit none -!------------------------------------------------------------------------------- - real :: ptke - real,parameter :: pmin = 0.0,pmax = 1.0 - real,parameter :: a1 = 1.000, a2 = 0.070, & - a3 = 1.000, a4 = 0.142, a5 = 0.071 - real,parameter :: b1 = 2.0, b2 = 0.6666667 - real :: d,h,doh,num,den -! - if (h .ne. 0) then - doh=d/h - num=a1*(doh)**b1+a2*(doh)**b2 - den=a3*(doh)**b1+a4*(doh)**b2+a5 - ptke=num/den - else - ptke=1. - endif - ptke=max(ptke,pmin) - ptke=min(ptke,pmax) -! - return - end function -!------------------------------------------------------------------------------- -! -!------------------------------------------------------------------------------- -end module module_bl_shinhong -!------------------------------------------------------------------------------- + endif + + do i = its, ite + psfc_hv(i) = psfc(i,j) + znt_hv(i) = znt(i,j) + ust_hv(i) = ust(i,j) + wspd_hv(i) = wspd(i,j) + psim_hv(i) = psim(i,j) + psih_hv(i) = psih(i,j) + xland_hv(i) = xland(i,j) + dx_hv(i) = dx(i,j) + corf_hv(i) = corf(i,j) + hfx_hv(i) = hfx(i,j) + qfx_hv(i) = qfx(i,j) + br_hv(i) = br(i,j) + u10_hv(i) = u10(i,j) + v10_hv(i) = v10(i,j) + uoce_hv(i) = uoce(i,j) + voce_hv(i) = voce(i,j) + ctopo_hv(i) = ctopo(i,j) + ctopo2_hv(i) = ctopo2(i,j) + end do +! + call bl_shinhong_run(ux=u3d_hv,vx=v3d_hv & + ,tx=t3d_hv & + ,qvx=qv3d_hv,qcx=qc3d_hv,qix=qi3d_hv & + ,f_qc=flag_qc,f_qi=flag_qi & + ,nmix=nmix,qmix=qmix_hv & + ,p2d=p3d_hv,p2di=p3di_hv & + ,pi2d=pi3d_hv & + ,utnp=rublten_hv,vtnp=rvblten_hv & + ,ttnp=rthblten_hv,qvtnp=rqvblten_hv & + ,qctnp=rqcblten_hv,qitnp=rqiblten_hv & + ,qmixtnp=rqmixblten_hv & + ,cp=cp,g=g,rovcp=rovcp,rd=rd,rovg=rovg & + ,xlv=xlv,rv=rv & + ,ep1=ep1,ep2=ep2,karman=karman & + ,dz8w2d=dz8w_hv & + ,psfcpa=psfc_hv,znt=znt_hv,ust=ust_hv & + ,hpbl=hpbl_hv & + ,psim=psim_hv & + ,psih=psih_hv,xland=xland_hv & + ,hfx=hfx_hv,qfx=qfx_hv & + ,wspd=wspd_hv,br=br_hv & + ,dt=dt,kpbl1d=kpbl2d_hv & + ,exch_hx=exch_h_hv & + ,exch_mx=exch_m_hv & + ,wstar=wstar_hv & + ,delta=delta_hv & + ,if_shinhong_nonlocal_flux=shinhong_nonlocal_flux & + ,tke=tke_hv,el_pbl=el_hv,corf=corf_hv & + ,u10=u10_hv,v10=v10_hv & + ,uox=uoce_hv,vox=voce_hv & + ,rthraten=rthraten_hv & + ,if_scu_mixing=shinhong_scu_mixing & + ,if_dissipative_heating=shinhong_dissi_heating & + ,ctopo=ctopo_hv,ctopo2=ctopo2_hv & + ,dxmeter=dx_hv & + ,a_u=a_u_hv,a_v=a_v_hv,a_t=a_t_hv,a_q=a_q_hv,a_e=a_e_hv & + ,b_u=b_u_hv,b_v=b_v_hv,b_t=b_t_hv,b_q=b_q_hv,b_e=b_e_hv & + ,sfk=sfk_hv,vlk=vlk_hv,dlu=dl_u_hv,dlg=dlg_hv,frcurb=frcurb_hv & + ,flag_bep=flag_bep & + ,its=its,ite=ite,kte=kte,kme=kme & + ,errmsg=errmsg,errflg=errflg ) +! + ! Assign local data back to full-sized arrays. + ! Only required for the INTENT(OUT) or INTENT(INOUT) arrays. + + do n = 1, nmix + do k = kts, kte + do i = its, ite + rqmixblten(i,k,j,n) = rqmixblten_hv(i,k,n) + end do + end do + end do + + do k = kts, kte + do i = its, ite + rublten(i,k,j) = rublten_hv(i,k) + rvblten(i,k,j) = rvblten_hv(i,k) +#if (NEED_B4B_DURING_CCPP_TESTING == 1) + rthblten(i,k,j) = rthblten_hv(i,k)/pi3d_hv(i,k) +#elif (NEED_B4B_DURING_CCPP_TESTING != 1) + rthblten(i,k,j) = rthblten_hv(i,k) +#endif + rqvblten(i,k,j) = rqvblten_hv(i,k) + rqcblten(i,k,j) = rqcblten_hv(i,k) + rqiblten(i,k,j) = rqiblten_hv(i,k) + exch_h(i,k,j) = exch_h_hv(i,k) + exch_m(i,k,j) = exch_m_hv(i,k) + tke(i,k,j) = tke_hv(i,k) + el(i,k,j) = el_hv(i,k) + end do + end do + + do i = its, ite + u10(i,j) = u10_hv(i) + v10(i,j) = v10_hv(i) + hpbl(i,j) = hpbl_hv(i) + kpbl2d(i,j) = kpbl2d_hv(i) + wstar(i,j) = wstar_hv(i) + delta(i,j) = delta_hv(i) + end do + enddo + + end subroutine shinhong + +!================================================================================================================= + end module module_bl_shinhong +!================================================================================================================= diff --git a/phys/module_pbl_driver.F b/phys/module_pbl_driver.F index f787fbb19d..4aa476c86f 100644 --- a/phys/module_pbl_driver.F +++ b/phys/module_pbl_driver.F @@ -34,6 +34,9 @@ SUBROUTINE pbl_driver( & ,windfarm_wake_model, windfarm_overlap_method & ,ysu_topdown_pblmix & ,shinhong_tke_diag & + ,shinhong_scu_mixing & + ,shinhong_nonlocal_flux & + ,shinhong_ke_dissipation & ! OPTIONAL for TEMF scheme ,te_temf,km_temf,kh_temf & ,shf_temf,qf_temf,uw_temf,vw_temf & @@ -428,6 +431,9 @@ SUBROUTINE pbl_driver( & INTENT(IN ) :: LOWLYR ! LOGICAL, INTENT(IN ) :: warm_rain + LOGICAL, INTENT(IN ) :: shinhong_scu_mixing + LOGICAL, INTENT(IN ) :: shinhong_nonlocal_flux + LOGICAL, INTENT(IN ) :: shinhong_ke_dissipation LOGICAL, INTENT(IN ) :: is_CAMMGMP_used !BSINGH:01/31/2013: Added for CAMUWPBL LOGICAL, OPTIONAL, INTENT(IN ) :: restart,cycling !used by mynn @@ -1311,33 +1317,44 @@ SUBROUTINE pbl_driver( & PRESENT( hol ) ) THEN ! CALL shinhong( & - U3D=u_phytmp,V3D=v_phytmp,TH3D=th_phy,T3D=t_phy & + U3D=u_phytmp,V3D=v_phytmp,T3D=t_phy & ,QV3D=qv_curr,QC3D=qc_curr,QI3D=qi_curr & ,P3D=p_phy,P3DI=p8w,PI3D=pi_phy & ,RUBLTEN=rublten,RVBLTEN=rvblten & ,RTHBLTEN=rthblten,RQVBLTEN=rqvblten & ,RQCBLTEN=rqcblten,RQIBLTEN=rqiblten & - ,FLAG_QI=flag_qi & + ,FLAG_QI=flag_qi,FLAG_QC=flag_qc & ,CP=cp,G=g,ROVCP=rcp,RD=r_D,ROVG=rovg & ,DZ8W=dz8w,XLV=XLV,RV=r_v,PSFC=PSFC & - ,ZNU=znu,ZNW=znw,P_TOP=p_top & ,ZNT=znt,UST=ust,HPBL=pblh & ,PSIM=fm,PSIH=fhh,XLAND=xland & ,HFX=hfx,QFX=qfx & ,U10=u10,V10=v10 & + ,UOCE=uoce,VOCE=voce & ! paj - ,CTOPO=ctopo,CTOPO2=ctopo2 & - ,SHINHONG_TKE_DIAG=shinhong_tke_diag & + ,shinhong_nonlocal_flux=shinhong_nonlocal_flux & + ,shinhong_scu_mixing=shinhong_scu_mixing & + ,shinhong_dissi_heating=shinhong_ke_dissipation & + ,tke=tke_pbl,el=el_pbl,corf=f & + ,CTOPO=ctopo,CTOPO2=ctopo2,dx=dx2d & ,WSPD=wspd,BR=br,DT=dtbl,KPBL2D=kpbl & ,EP1=ep_1,EP2=ep_2,KARMAN=karman & - ,EXCH_H=exch_h,REGIME=regime & -! for grims shallow convection with shinhongpbl + ,EXCH_H=exch_h,EXCH_M=exch_m & + ,RTHRATEN=RTHRATEN & +! for multilayer UCM + ,IDIFF=idiff,FLAG_BEP=flag_bep,FRC_URB2D=frc_urb2d & + ,A_U_BEP=a_u_bep,A_V_BEP=a_v_bep,A_T_BEP=a_t_bep & + ,A_Q_BEP=a_q_bep & + ,A_E_BEP=a_e_bep,B_U_BEP=b_u_bep,B_V_BEP=b_v_bep & + ,B_T_BEP=b_t_bep,B_Q_BEP=b_q_bep & + ,B_E_BEP=b_e_bep,DLG_BEP=dlg_bep & + ,DL_U_BEP=dl_u_bep,SF_BEP=sf_bep,VL_BEP=vl_bep & +! for grims shallow convection ,WSTAR=wstar,DELTA=delta & - ,TKE_PBL=tke_pbl,EL_PBL=el_pbl,CORF=f & + ,errmsg=errmsg,errflg=errflg & ,IDS=ids,IDE=ide,JDS=jds,JDE=jde,KDS=kds,KDE=kde & ,IMS=ims,IME=ime,JMS=jms,JME=jme,KMS=kms,KME=kme & ,ITS=its,ITE=ite,JTS=jts,JTE=jte,KTS=kts,KTE=kte & - ,DX=dx,DY=dy & ) ELSE WRITE ( message , FMT = '(A,7(L1,1X))' ) & diff --git a/phys/module_physics_init.F b/phys/module_physics_init.F index 0c25ac1ec6..9e604a4f37 100644 --- a/phys/module_physics_init.F +++ b/phys/module_physics_init.F @@ -3892,16 +3892,18 @@ SUBROUTINE bl_init(STEPBL,BLDT,DT,RUBLTEN,RVBLTEN,RTHBLTEN, & CASE (SHINHONGSCHEME) if(isfc .ne. 1)CALL wrf_error_fatal & ( 'module_physics_init: Use sf_sfclay_physics= 1 or 91 for this pbl option' ) - IF ((SF_URBAN_PHYSICS.eq.2).OR.(SF_URBAN_PHYSICS.EQ.3)) CALL wrf_error_fatal & - ( 'module_physics_init: use ysu (option1), myj (option 2), or boulac (option 8) with BEP/BEM urban scheme' ) - CALL shinhonginit(RUBLTEN,RVBLTEN,RTHBLTEN,RQVBLTEN,& - RQCBLTEN,RQIBLTEN,TKE_PBL,P_QI, & - PARAM_FIRST_SCALAR, & - restart, & - allowed_to_read , & - ids, ide, jds, jde, kds, kde, & - ims, ime, jms, jme, kms, kme, & - its, ite, jts, jte, kts, kte ) +! if(isfc .ne. 1)CALL wrf_error_fatal & +! ( 'module_physics_init: Use sf_sfclay_physics= 1 or 91 for this pbl option' ) +! IF ((SF_URBAN_PHYSICS.eq.2).OR.(SF_URBAN_PHYSICS.EQ.3)) CALL wrf_error_fatal & +! ( 'module_physics_init: use ysu (option1), myj (option 2), or boulac (option 8) with BEP/BEM urban scheme' ) +! CALL shinhonginit(RUBLTEN,RVBLTEN,RTHBLTEN,RQVBLTEN,& +! RQCBLTEN,RQIBLTEN,TKE_PBL,P_QI, & +! PARAM_FIRST_SCALAR, & +! restart, & +! allowed_to_read , & +! ids, ide, jds, jde, kds, kde, & +! ims, ime, jms, jme, kms, kme, & +! its, ite, jts, jte, kts, kte ) CASE (MRFSCHEME) if(isfc .ne. 1)CALL wrf_error_fatal & ( 'module_physics_init: Use sf_sfclay_physics= 1 or 91 for this pbl option' ) diff --git a/phys/physics_mmm b/phys/physics_mmm index 0ea59b1cd6..2f8c1fa0d6 160000 --- a/phys/physics_mmm +++ b/phys/physics_mmm @@ -1 +1 @@ -Subproject commit 0ea59b1cd673006ee7a9a9958c533a6a0e354243 +Subproject commit 2f8c1fa0d66a54c5de836c49904f524f7a89a1d4