diff --git a/arch/postamble b/arch/postamble index 936f0405c8..aa55662073 100644 --- a/arch/postamble +++ b/arch/postamble @@ -203,6 +203,13 @@ wrfio_esmf : fi $(FC) -o $@ -c $(FCFLAGS) $(OMP) $(MODULE_DIRS) $(PROMOTION) $(FCSUFFIX) $*.f90 +.F90.o: + $(RM) $@ + sed -e "s/^\!.*'.*//" -e "s/^ *\!.*'.*//" $*.F90 > $*.G + $(CPP) -I$(WRF_SRC_ROOT_DIR)/inc $(CPPFLAGS) $(OMPCPP) $*.G > $*.bb + $(SED_FTN) $*.bb | $(CPP) $(TRADFLAG) > $*.f90 + $(RM) $*.G $*.bb + $(FC) -o $@ -c $(FCFLAGS) $(OMP) $(MODULE_DIRS) $(PROMOTION) $(FCSUFFIX) $*.f90 .F.f90: $(RM) $@ @@ -211,6 +218,13 @@ wrfio_esmf : $(CPP) -I$(WRF_SRC_ROOT_DIR)/inc $(CPPFLAGS) $*.H > $@ $(RM) $*.G $*.H +.F90.f90: + $(RM) $@ + sed -e "s/^\!.*'.*//" -e "s/^ *\!.*'.*//" $*.F90 > $*.G + $(SED_FTN) $*.G > $*.H + $(CPP) -I$(WRF_SRC_ROOT_DIR)/inc $(CPPFLAGS) $*.H > $@ + $(RM) $*.G $*.H + .f90.o: $(RM) $@ $(FC) -o $@ -c $(FCFLAGS) $(PROMOTION) $(FCSUFFIX) $*.f90 diff --git a/arch/preamble b/arch/preamble index 4543411e6f..4ae897d496 100644 --- a/arch/preamble +++ b/arch/preamble @@ -17,7 +17,7 @@ SHELL = /bin/sh DEVTOP = `pwd` LIBINCLUDE = . -.SUFFIXES: .F .i .o .f90 .c +.SUFFIXES: .F .i .o .f90 .c .F90 #### Get core settings from environment (set in compile script) #### Note to add a core, this has to be added to. diff --git a/clean b/clean index 6ce8142cd6..2fd6453194 100755 --- a/clean +++ b/clean @@ -3,12 +3,12 @@ set nonomatch -foreach dir ( frame chem share dyn_em phys cmaq main tools wrftladj ) +foreach dir ( frame chem share dyn_em phys phys/physics_mmm cmaq main tools wrftladj ) if ( -d $dir ) then if ( $dir == cmaq ) then ( cd $dir ; echo $dir ; /bin/rm -f *.o *.mod ) >& /dev/null else - ( cd $dir ; echo $dir ; /bin/rm -f core wrf *.f90 *.exe *.kmo *.mod *.o *.obj *.inc *.F90 *.a \ + ( cd $dir ; echo $dir ; /bin/rm -f core wrf *.f90 *.exe *.kmo *.mod *.o *.obj *.inc *.a \ db_* Warnings module_state_description.F module_dm.F gmeta \ wrfdata whatiread rsl.* show_domain* ) >& /dev/null endif diff --git a/main/depend.common b/main/depend.common index 5bcdc1720a..f40adae791 100644 --- a/main/depend.common +++ b/main/depend.common @@ -299,7 +299,7 @@ module_bl_qnsepbl.o: \ module_progtm.o: \ - module_gfs_machine.o + module_gfs_machine.o module_bl_gfs.o: \ @@ -321,6 +321,10 @@ module_bl_mynn_wrapper.o: \ module_bl_mynn_common.o +module_bl_gwdo.o: \ + physics_mmm/bl_gwdo.o + + module_cam_upper_bc.o: \ module_cam_shr_kind_mod.o \ module_cam_support.o @@ -1477,6 +1481,11 @@ module_mp_gsfcgce_4ice_nuwrf.o: \ module_checkerror.o \ module_mp_radar.o \ + +physics_mmm/bl_gwdo.o : \ + ccpp_kind_types.o + + # End of DEPENDENCIES for phys diff --git a/phys/Makefile b/phys/Makefile index dd88587d3c..ffe85ea1c5 100644 --- a/phys/Makefile +++ b/phys/Makefile @@ -233,19 +233,23 @@ DIAGNOSTIC_MODULES_EM = \ module_diag_trad_fields.o \ module_diag_solar.o +PHYSMMM_MODULES = \ + ./physics_mmm/bl_gwdo.o + OBJS = LIBTARGET = physics TARGETDIR = ./ $(LIBTARGET) : - $(MAKE) $(J) non_nmm ; \ + $(MAKE) $(J) non_nmm shared_phys ; \ $(AR) $(ARFLAGS) ../main/$(LIBWRFLIB) $(MODULES) $(OBJS) \ - $(FIRE_MODULES) $(DIAGNOSTIC_MODULES_EM) + $(FIRE_MODULES) $(DIAGNOSTIC_MODULES_EM) $(PHYSMMM_MODULES) include ../configure.wrf -non_nmm : $(MODULES) $(FIRE_MODULES) $(WIND_MODULES) $(OBJS) $(DIAGNOSTIC_MODULES_EM) +non_nmm : $(MODULES) $(FIRE_MODULES) $(OBJS) $(DIAGNOSTIC_MODULES_EM) +shared_phys : $(PHYSMMM_MODULES) submodules : @if [ \( ! -f module_sf_noahmpdrv.F \) -o \( ! -f module_sf_noahmp_glacier.F \) -o \ diff --git a/phys/module_bl_gwdo.F b/phys/module_bl_gwdo.F index c81e67c33e..81026c6404 100644 --- a/phys/module_bl_gwdo.F +++ b/phys/module_bl_gwdo.F @@ -1,21 +1,32 @@ -!WRF:model_layer:physics -! -module module_bl_gwdo -contains -!------------------------------------------------------------------------------- - subroutine gwdo(u3d,v3d,t3d,qv3d,p3d,p3di,pi3d,z, & - rublten,rvblten, & - dtaux3d,dtauy3d,dusfcg,dvsfcg, & - var2d,oc12d,oa2d1,oa2d2,oa2d3,oa2d4,ol2d1,ol2d2,ol2d3,ol2d4, & - sina,cosa,znu,znw,p_top, & - cp,g,rd,rv,ep1,pi, & - dt,dx,kpbl2d,itimestep, & - ids,ide, jds,jde, kds,kde, & - ims,ime, jms,jme, kms,kme, & - its,ite, jts,jte, kts,kte) -!------------------------------------------------------------------------------- - implicit none -!------------------------------------------------------------------------------- +!================================================================================================================= + module module_bl_gwdo + use ccpp_kind_types,only: kind_phys + + use bl_gwdo,only: bl_gwdo_run + + + implicit none + private + public:: gwdo + + + contains + + +!================================================================================================================= + subroutine gwdo(u3d,v3d,t3d,qv3d,p3d,p3di,pi3d,z, & + rublten,rvblten, & + dtaux3d,dtauy3d,dusfcg,dvsfcg, & + var2d,oc12d,oa2d1,oa2d2,oa2d3,oa2d4,ol2d1,ol2d2,ol2d3,ol2d4, & + sina,cosa,znu,znw,p_top, & + cp,g,rd,rv,ep1,pi, & + dt,dx,kpbl2d,itimestep, & + ids,ide, jds,jde, kds,kde, & + ims,ime, jms,jme, kms,kme, & + its,ite, jts,jte, kts,kte, & + errmsg,errflg & + ) +!================================================================================================================= ! !-- u3d 3d u-velocity interpolated to theta points (m/s) !-- v3d 3d v-velocity interpolated to theta points (m/s) @@ -56,672 +67,177 @@ subroutine gwdo(u3d,v3d,t3d,qv3d,p3d,p3di,pi3d,z, & !-- kts start index for k in tile !-- kte end index for k in tile ! -!------------------------------------------------------------------------------- - integer, intent(in ) :: ids,ide, jds,jde, kds,kde, & - ims,ime, jms,jme, kms,kme, & - its,ite, jts,jte, kts,kte - integer, intent(in ) :: itimestep -! - real, intent(in ) :: dt,dx,cp,g,rd,rv,ep1,pi -! - real, dimension( ims:ime, kms:kme, jms:jme ) , & - intent(in ) :: qv3d, & - p3d, & - pi3d, & - t3d, & - z - real, dimension( ims:ime, kms:kme, jms:jme ) , & - intent(in ) :: p3di -! - real, dimension( ims:ime, kms:kme, jms:jme ) , & - intent(inout) :: rublten, & - rvblten - real, dimension( ims:ime, kms:kme, jms:jme ) , & - intent(inout) :: dtaux3d, & - dtauy3d -! - real, dimension( ims:ime, kms:kme, jms:jme ) , & - intent(in ) :: u3d, & - v3d -! - integer, dimension( ims:ime, jms:jme ) , & - intent(in ) :: kpbl2d - real, dimension( ims:ime, jms:jme ) , & - intent(inout ) :: dusfcg, & - dvsfcg -! - real, dimension( ims:ime, jms:jme ) , & - intent(in ) :: var2d, & - oc12d, & - oa2d1,oa2d2,oa2d3,oa2d4, & - ol2d1,ol2d2,ol2d3,ol2d4, & - sina,cosa -! - real, dimension( kms:kme ) , & - optional , & - intent(in ) :: znu, & - znw -! - real, optional, intent(in ) :: p_top -! -!local -! - real, dimension( its:ite, kts:kte ) :: delprsi, & - pdh - real, dimension( its:ite, kts:kte ) :: ugeo, vgeo, dudt, dvdt, dtaux, dtauy - real, dimension( its:ite ) :: dusfc, dvsfc - real, dimension( its:ite, kts:kte+1 ) :: pdhi - real, dimension( its:ite, 4 ) :: oa4, & - ol4 - integer :: i,j,k,kpblmax -! - do k = kts,kte - if (znu(k).gt.0.6) kpblmax = k + 1 - enddo -! - 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 - delprsi(i,k) = pdhi(i,k)-pdhi(i,k+1) -! rotate winds to zonal/meridional - ugeo(i,k) = u3d(i,k,j)*cosa(i,j) - v3d(i,k,j)*sina(i,j) - vgeo(i,k) = u3d(i,k,j)*sina(i,j) + v3d(i,k,j)*cosa(i,j) - dudt(i,k) = 0.0 - dvdt(i,k) = 0.0 - enddo - enddo - do i = its,ite - oa4(i,1) = oa2d1(i,j) - oa4(i,2) = oa2d2(i,j) - oa4(i,3) = oa2d3(i,j) - oa4(i,4) = oa2d4(i,j) - ol4(i,1) = ol2d1(i,j) - ol4(i,2) = ol2d2(i,j) - ol4(i,3) = ol2d3(i,j) - ol4(i,4) = ol2d4(i,j) - enddo - call gwdo2d(dudt=dudt(its,kts),dvdt=dvdt(its,kts) & - ,dtaux2d=dtaux(its,kts),dtauy2d=dtauy(its,kts) & - ,u1=ugeo(its,kts),v1=vgeo(its,kts) & - ,t1=t3d(ims,kms,j),q1=qv3d(ims,kms,j) & - ,del=delprsi(its,kts) & - ,prsi=pdhi(its,kts) & - ,prsl=pdh(its,kts),prslk=pi3d(ims,kms,j) & - ,zl=z(ims,kms,j) & - ,kpblmax=kpblmax & - ,var=var2d(ims,j),oc1=oc12d(ims,j) & - ,oa4=oa4,ol4=ol4 & - ,dusfc=dusfc(its),dvsfc=dvsfc(its) & - ,g_=g,cp_=cp,rd_=rd,rv_=rv,fv_=ep1,pi_=pi & - ,dxmeter=dx,deltim=dt & - ,kpbl=kpbl2d(ims,j),lat=j & - ,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 -! rotate tendencies from zonal/meridional to model grid - rublten(i,k,j) = rublten(i,k,j)+dudt(i,k)*cosa(i,j) + dvdt(i,k)*sina(i,j) - rvblten(i,k,j) = rvblten(i,k,j)-dudt(i,k)*sina(i,j) + dvdt(i,k)*cosa(i,j) - dtaux3d(i,k,j) = dtaux(i,k)*cosa(i,j) + dtauy(i,k)*sina(i,j) - dtauy3d(i,k,j) =-dtaux(i,k)*sina(i,j) + dtauy(i,k)*cosa(i,j) - if(k.eq.kts)then - dusfcg(i,j) = dusfc(i)*cosa(i,j) + dvsfc(i)*sina(i,j) - dvsfcg(i,j) =-dusfc(i)*sina(i,j) + dvsfc(i)*cosa(i,j) - endif - enddo - enddo - enddo -! - end subroutine gwdo -!------------------------------------------------------------------------------- -!------------------------------------------------------------------------------- - subroutine gwdo2d(dudt, dvdt, dtaux2d, dtauy2d, & - u1, v1, t1, q1, & - del, & - prsi, prsl, prslk, zl, & - kpblmax, & - var, oc1, oa4, ol4, dusfc, dvsfc, & - g_, cp_, rd_, rv_, fv_, pi_, & - dxmeter, deltim, kpbl, lat, & - ids, ide, jds, jde, kds, kde, & - ims, ime, jms, jme, kms, kme, & - its, ite, jts, jte, kts, kte) -!------------------------------------------------------------------------------- -! -! abstract : -! this code handles the time tendencies of u v due to the effect of -! mountain induced gravity wave drag from sub-grid scale orography. -! this routine not only treats the traditional upper-level wave breaking due -! to mountain variance (alpert 1988), but also the enhanced -! lower-tropospheric wave breaking due to mountain convexity and asymmetry -! (kim and arakawa 1995). thus, in addition to the terrain height data -! in a model grid gox, additional 10-2d topographic statistics files are -! needed, including orographic standard deviation (var), convexity (oc1), -! asymmetry (oa4) and ol (ol4). these data sets are prepared based on the -! 30 sec usgs orography (hong 1999). the current scheme was implmented as in -! choi and hong (2015), which names kim gwdo since it was developed by -! kiaps staffs for kiaps integrated model system (kim). the scheme -! additionally includes the effects of orographic anisotropy and -! flow-blocking drag. -! coded by song-you hong and young-joon kim and implemented by song-you hong -! -! history log : -! 2015-07-01 hyun-joo choi add flow-blocking drag and orographic anisotropy -! -! references : -! choi and hong (2015), j. geophys. res. -! hong et al. (2008), wea. forecasting -! kim and doyle (2005), q. j. r. meteor. soc. -! kim and arakawa (1995), j. atmos. sci. -! alpet et al. (1988), NWP conference -! hong (1999), NCEP office note 424 -! -! input : -! dudt, dvdt - non-lin tendency for u and v wind component -! u1, v1 - zonal and meridional wind m/sec at t0-dt -! t1 - temperature deg k at t0-dt -! q1 - mixing ratio at t0-dt -! deltim - time step (s) -! del - positive increment of pressure across layer (pa) -! kpblmax, kpbl - vertical index of pbl height -! prslk, zl, prsl, prsi - pressure and height variables -! oa4, ol4, omax, var, oc1 - orographic statistics -! -! output : -! dudt, dvdt - wind tendency due to gwdo -! dtaux2d, dtauy2d - diagnoised orographic gwd -! dusfc, dvsfc - gw stress -! -!------------------------------------------------------------------------------- - implicit none -! - integer , intent(in ) :: lat, kpblmax, & - ids, ide, jds, jde, & - kds, kde, ims, ime, & - jms, jme, kms, kme, & - its, ite, jts, jte, & - kts, kte - integer, dimension(ims:ime) , intent(in ) :: kpbl - real , intent(in ) :: g_, pi_, rd_, rv_, fv_,& - cp_, deltim - real , intent(in ) :: dxmeter - real, dimension(its:ite,kts:kte) , intent(inout) :: dudt, dvdt - real, dimension(its:ite,kts:kte) , intent( out) :: dtaux2d, dtauy2d - real, dimension(its:ite,kts:kte) , intent(in ) :: u1, v1 - real, dimension(ims:ime,kms:kme) , intent(in ) :: t1, q1, prslk, zl -! - real, dimension(its:ite,kts:kte) , intent(in ) :: prsl, del - real, dimension(its:ite,kts:kte+1), intent(in ) :: prsi - real, dimension(its:ite,4) , intent(in ) :: oa4, ol4 -! - real, dimension(ims:ime) , intent(in ) :: var, oc1 - real, dimension(its:ite) , intent( out) :: dusfc, dvsfc -! - real, parameter :: ric = 0.25 ! critical richardson number - real, parameter :: dw2min = 1. - real, parameter :: rimin = -100. - real, parameter :: bnv2min = 1.0e-5 - real, parameter :: efmin = 0.0 - real, parameter :: efmax = 10.0 - real, parameter :: xl = 4.0e4 - real, parameter :: critac = 1.0e-5 - real, parameter :: gmax = 1. - real, parameter :: veleps = 1.0 - real, parameter :: frc = 1.0 - real, parameter :: ce = 0.8 - real, parameter :: cg = 0.5 - integer,parameter :: kpblmin = 2 -! -! local variables -! - integer :: latd,lond - integer :: i,k,lcap,lcapp1,nwd,idir, & - klcap,kp1,ikount,kk -! - real :: fdir,cs,rcsks, & - wdir,ti,rdz,temp,tem2,dw2,shr2,bvf2,rdelks, & - wtkbj,tem,gfobnv,hd,fro,rim,temc,tem1,efact, & - temv,dtaux,dtauy -! - logical, dimension(its:ite) :: ldrag, icrilv, flag,kloop1 - real, dimension(its:ite) :: coefm -! - real, dimension(its:ite) :: taub, xn, yn, ubar, vbar, fr, & - ulow, rulow, bnv, oa, ol, rhobar, & - dtfac, brvf, xlinv, delks,delks1, & - zlowtop,cleff - real, dimension(its:ite,kts:kte+1) :: taup - real, dimension(its:ite,kts:kte-1) :: velco - real, dimension(its:ite,kts:kte) :: bnv2, usqj, taud, rho, vtk, vtj -! - integer, dimension(its:ite) :: kbl, klowtop - integer, parameter :: mdir=8 - integer, dimension(mdir) :: nwdir - data nwdir/6,7,5,8,2,3,1,4/ -! -! variables for flow-blocking drag -! - real, parameter :: frmax = 10. - real, parameter :: olmin = 1.0e-5 - real, parameter :: odmin = 0.1 - real, parameter :: odmax = 10. -! - real :: fbdcd - real :: zblk, tautem - real :: fbdpe, fbdke - real, dimension(its:ite) :: delx, dely - real, dimension(its:ite,4) :: dxy4, dxy4p - real, dimension(4) :: ol4p - real, dimension(its:ite) :: dxy, dxyp, olp, od - real, dimension(its:ite,kts:kte+1) :: taufb -! - integer, dimension(its:ite) :: komax - integer :: kblk -!------------------------------------------------------------------------------- -! -! constants -! - lcap = kte - lcapp1 = lcap + 1 - fdir = mdir / (2.0*pi_) -! -! calculate length of grid for flow-blocking drag -! - delx(its:ite) = dxmeter - dely(its:ite) = dxmeter - dxy4(its:ite,1) = delx(its:ite) - dxy4(its:ite,2) = dely(its:ite) - dxy4(its:ite,3) = sqrt(delx(its:ite)**2. + dely(its:ite)**2.) - dxy4(its:ite,4) = dxy4(its:ite,3) - dxy4p(its:ite,1) = dxy4(its:ite,2) - dxy4p(its:ite,2) = dxy4(its:ite,1) - dxy4p(its:ite,3) = dxy4(its:ite,4) - dxy4p(its:ite,4) = dxy4(its:ite,3) -! - cleff(its:ite) = dxmeter -! -! initialize arrays -! - ldrag = .false. ; icrilv = .false. ; flag = .true. -! - klowtop = 0 ; kbl = 0 -! - dtaux = 0. ; dtauy = 0. ; xn = 0. ; yn = 0. - ubar = 0. ; vbar = 0. ; rhobar = 0. ; ulow = 0. - oa = 0. ; ol = 0. ; taub = 0. -! - usqj = 0. ; bnv2 = 0. ; vtj = 0. ; vtk = 0. - taup = 0. ; taud = 0. ; dtaux2d = 0. ; dtauy2d = 0. -! - dtfac = 1.0 ; xlinv = 1.0/xl -! -! initialize arrays for flow-blocking drag -! - komax = 0 - taufb = 0.0 -! - do k = kts,kte - do i = its,ite - vtj(i,k) = t1(i,k) * (1.+fv_*q1(i,k)) - vtk(i,k) = vtj(i,k) / prslk(i,k) - rho(i,k) = 1./rd_ * prsl(i,k) / vtj(i,k) ! density kg/m**3 - enddo - enddo -! - do i = its,ite - zlowtop(i) = 2. * var(i) - enddo -! - do i = its,ite - kloop1(i) = .true. - enddo -! - do k = kts+1,kte - do i = its,ite - if (kloop1(i).and.zl(i,k)-zl(i,1).ge.zlowtop(i)) then - klowtop(i) = k+1 - kloop1(i) = .false. - endif - enddo - enddo -! - do i = its,ite -! -! determine reference level: 2*var -! - kbl(i) = klowtop(i) - kbl(i) = max(min(kbl(i),kpblmax),kpblmin) - enddo -! -! determine the level of maximum orographic height -! - komax(:) = kbl(:) -! - do i = its,ite - delks(i) = 1.0 / (prsi(i,1) - prsi(i,kbl(i))) - delks1(i) = 1.0 / (prsl(i,1) - prsl(i,kbl(i))) - enddo -! -! compute low level averages within pbl -! - do k = kts,kpblmax - do i = its,ite - if (k.lt.kbl(i)) then - rcsks = del(i,k) * delks(i) - rdelks = del(i,k) * delks(i) - ubar(i) = ubar(i) + rcsks * u1(i,k) ! pbl u mean - vbar(i) = vbar(i) + rcsks * v1(i,k) ! pbl v mean - rhobar(i) = rhobar(i) + rdelks * rho(i,k) ! pbl rho mean - endif - enddo - enddo -! -! figure out low-level horizontal wind direction -! -! nwd 1 2 3 4 5 6 7 8 -! wd w s sw nw e n ne se -! - do i = its,ite - wdir = atan2(ubar(i),vbar(i)) + pi_ - idir = mod(nint(fdir*wdir),mdir) + 1 - nwd = nwdir(idir) - oa(i) = (1-2*int( (nwd-1)/4 )) * oa4(i,mod(nwd-1,4)+1) - ol(i) = ol4(i,mod(nwd-1,4)+1) -! -! compute orographic width along (ol) and perpendicular (olp) the wind direction -! - ol4p(1) = ol4(i,2) - ol4p(2) = ol4(i,1) - ol4p(3) = ol4(i,4) - ol4p(4) = ol4(i,3) - olp(i) = ol4p(mod(nwd-1,4)+1) -! -! compute orographic direction (horizontal orographic aspect ratio) -! - od(i) = olp(i)/max(ol(i),olmin) - od(i) = min(od(i),odmax) - od(i) = max(od(i),odmin) -! -! compute length of grid in the along(dxy) and cross(dxyp) wind directions -! - dxy(i) = dxy4(i,MOD(nwd-1,4)+1) - dxyp(i) = dxy4p(i,MOD(nwd-1,4)+1) - enddo -! -! saving richardson number in usqj for migwdi -! - do k = kts,kte-1 - do i = its,ite - ti = 2.0 / (t1(i,k)+t1(i,k+1)) - rdz = 1./(zl(i,k+1) - zl(i,k)) - tem1 = u1(i,k) - u1(i,k+1) - tem2 = v1(i,k) - v1(i,k+1) - dw2 = tem1*tem1 + tem2*tem2 - shr2 = max(dw2,dw2min) * rdz * rdz - bvf2 = g_*(g_/cp_+rdz*(vtj(i,k+1)-vtj(i,k))) * ti - usqj(i,k) = max(bvf2/shr2,rimin) - bnv2(i,k) = 2.0*g_*rdz*(vtk(i,k+1)-vtk(i,k))/(vtk(i,k+1)+vtk(i,k)) - enddo - enddo -! -! compute the "low level" or 1/3 wind magnitude (m/s) -! - do i = its,ite - ulow(i) = max(sqrt(ubar(i)*ubar(i) + vbar(i)*vbar(i)), 1.0) - rulow(i) = 1./ulow(i) - enddo -! - do k = kts,kte-1 - do i = its,ite - velco(i,k) = 0.5 * ((u1(i,k)+u1(i,k+1)) * ubar(i) & - + (v1(i,k)+v1(i,k+1)) * vbar(i)) - velco(i,k) = velco(i,k) * rulow(i) - if ((velco(i,k).lt.veleps) .and. (velco(i,k).gt.0.)) then - velco(i,k) = veleps - endif - enddo - enddo -! -! no drag when critical level in the base layer -! - do i = its,ite - ldrag(i) = velco(i,1).le.0. - enddo -! -! no drag when velco.lt.0 -! - do k = kpblmin,kpblmax - do i = its,ite - if (k .lt. kbl(i)) ldrag(i) = ldrag(i).or. velco(i,k).le.0. - enddo - enddo -! -! the low level weighted average ri is stored in usqj(1,1; im) -! the low level weighted average n**2 is stored in bnv2(1,1; im) -! this is called bnvl2 in phy_gwd_alpert_sub not bnv2 -! rdelks (del(k)/delks) vert ave factor so we can * instead of / -! - do i = its,ite - wtkbj = (prsl(i,1)-prsl(i,2)) * delks1(i) - bnv2(i,1) = wtkbj * bnv2(i,1) - usqj(i,1) = wtkbj * usqj(i,1) - enddo -! - do k = kpblmin,kpblmax - do i = its,ite - if (k .lt. kbl(i)) then - rdelks = (prsl(i,k)-prsl(i,k+1)) * delks1(i) - bnv2(i,1) = bnv2(i,1) + bnv2(i,k) * rdelks - usqj(i,1) = usqj(i,1) + usqj(i,k) * rdelks - endif - enddo - enddo -! - do i = its,ite - ldrag(i) = ldrag(i) .or. bnv2(i,1).le.0.0 - ldrag(i) = ldrag(i) .or. ulow(i).eq.1.0 - ldrag(i) = ldrag(i) .or. var(i) .le. 0.0 - enddo -! -! set all ri low level values to the low level value -! - do k = kpblmin,kpblmax - do i = its,ite - if (k .lt. kbl(i)) usqj(i,k) = usqj(i,1) - enddo - enddo -! - do i = its,ite - if (.not.ldrag(i)) then - bnv(i) = sqrt( bnv2(i,1) ) - fr(i) = bnv(i) * rulow(i) * var(i) * od(i) - fr(i) = min(fr(i),frmax) - xn(i) = ubar(i) * rulow(i) - yn(i) = vbar(i) * rulow(i) - endif - enddo -! -! compute the base level stress and store it in taub -! calculate enhancement factor, number of mountains & aspect -! ratio const. use simplified relationship between standard -! deviation & critical hgt -! - do i = its,ite - if (.not. ldrag(i)) then - efact = (oa(i) + 2.) ** (ce*fr(i)/frc) - efact = min( max(efact,efmin), efmax ) - coefm(i) = (1. + ol(i)) ** (oa(i)+1.) - xlinv(i) = coefm(i) / cleff(i) - tem = fr(i) * fr(i) * oc1(i) - gfobnv = gmax * tem / ((tem + cg)*bnv(i)) - taub(i) = xlinv(i) * rhobar(i) * ulow(i) * ulow(i) & - * ulow(i) * gfobnv * efact - else - taub(i) = 0.0 - xn(i) = 0.0 - yn(i) = 0.0 - endif - enddo -! -! now compute vertical structure of the stress. -! - do k = kts,kpblmax - do i = its,ite - if (k .le. kbl(i)) taup(i,k) = taub(i) - enddo - enddo -! - do k = kpblmin, kte-1 ! vertical level k loop! - kp1 = k + 1 - do i = its,ite -! -! unstablelayer if ri < ric -! unstable layer if upper air vel comp along surf vel <=0 (crit lay) -! at (u-c)=0. crit layer exists and bit vector should be set (.le.) -! - if (k .ge. kbl(i)) then - icrilv(i) = icrilv(i) .or. ( usqj(i,k) .lt. ric) & - .or. (velco(i,k) .le. 0.0) - brvf(i) = max(bnv2(i,k),bnv2min) ! brunt-vaisala frequency squared - brvf(i) = sqrt(brvf(i)) ! brunt-vaisala frequency - endif - enddo -! - do i = its,ite - if (k .ge. kbl(i) .and. (.not. ldrag(i))) then - if (.not.icrilv(i) .and. taup(i,k) .gt. 0.0 ) then - temv = 1.0 / velco(i,k) - tem1 = coefm(i)/dxy(i)*(rho(i,kp1)+rho(i,k))*brvf(i)*velco(i,k)*0.5 - hd = sqrt(taup(i,k) / tem1) - fro = brvf(i) * hd * temv -! -! rim is the minimum-richardson number by shutts (1985) -! - tem2 = sqrt(usqj(i,k)) - tem = 1. + tem2 * fro - rim = usqj(i,k) * (1.-fro) / (tem * tem) -! -! check stability to employ the 'saturation hypothesis' -! of lindzen (1981) except at tropospheric downstream regions -! - if (rim .le. ric) then ! saturation hypothesis! - if ((oa(i) .le. 0.).or.(kp1 .ge. kpblmin )) then - temc = 2.0 + 1.0 / tem2 - hd = velco(i,k) * (2.*sqrt(temc)-temc) / brvf(i) - taup(i,kp1) = tem1 * hd * hd - endif - else ! no wavebreaking! - taup(i,kp1) = taup(i,k) - endif - endif - endif - enddo - enddo -! - if (lcap.lt.kte) then - do klcap = lcapp1,kte - do i = its,ite - taup(i,klcap) = prsi(i,klcap) / prsi(i,lcap) * taup(i,lcap) - enddo - enddo - endif - do i = its,ite - if (.not.ldrag(i)) then -! -! determine the height of flow-blocking layer -! - kblk = 0 - fbdpe = 0.0 - fbdke = 0.0 - do k = kte, kpblmin, -1 - if (kblk.eq.0 .and. k.le.kbl(i)) then - fbdpe = fbdpe + bnv2(i,k)*(zl(i,kbl(i))-zl(i,k)) & - *del(i,k)/g_/rho(i,k) - fbdke = 0.5*(u1(i,k)**2.+v1(i,k)**2.) -! -! apply flow-blocking drag when fbdpe >= fbdke -! - if (fbdpe.ge.fbdke) then - kblk = k - kblk = min(kblk,kbl(i)) - zblk = zl(i,kblk)-zl(i,kts) - endif - endif +!================================================================================================================= + +!--- input arguments: + 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),dimension(ims:ime,jms:jme):: kpbl2d + + real(kind=kind_phys),intent(in):: dt,cp,g,rd,rv,ep1,pi + real(kind=kind_phys),intent(in),optional:: p_top + + real(kind=kind_phys),intent(in),dimension(kms:kme),optional:: & + znu, & + znw + + real(kind=kind_phys),intent(in),dimension(ims:ime,jms:jme):: & + dx, & + var2d, & + oc12d, & + oa2d1,oa2d2,oa2d3,oa2d4, & + ol2d1,ol2d2,ol2d3,ol2d4, & + sina,cosa + + real(kind=kind_phys),intent(in),dimension(ims:ime,kms:kme,jms:jme):: & + qv3d, & + p3d, & + pi3d, & + t3d, & + u3d, & + v3d, & + z + + real(kind=kind_phys),intent(in),dimension(ims:ime,kms:kme,jms:jme):: & + p3di + +!--- output arguments: + character(len=*),intent(out):: errmsg + + integer,intent(out):: errflg + + real(kind=kind_phys),intent(out),dimension(ims:ime,jms:jme):: & + dusfcg, & + dvsfcg + + real(kind=kind_phys),intent(out),dimension(ims:ime,kms:kme,jms:jme ):: & + dtaux3d, & + dtauy3d + +!--- inout arguments: + real(kind=kind_phys),intent(inout),dimension(ims:ime,kms:kme,jms:jme):: & + rublten, & + rvblten + +!--- local variables and arrays: + integer:: i,j,k + + real(kind=kind_phys),dimension(its:ite):: & + var2d_hv,oc12d_hv,dx_hv,sina_hv,cosa_hv + real(kind=kind_phys),dimension(its:ite):: & + oa2d1_hv,oa2d2_hv,oa2d3_hv,oa2d4_hv,ol2d1_hv,ol2d2_hv,ol2d3_hv,ol2d4_hv + real(kind=kind_phys),dimension(its:ite):: & + dusfcg_hv,dvsfcg_hv + + real(kind=kind_phys),dimension(its:ite,kts:kte):: & + u3d_hv,v3d_hv,t3d_hv,qv3d_hv,pi3d_hv,p3d_hv,z_hv + + real(kind=kind_phys),dimension(its:ite,kts:kte):: & + rublten_hv,rvblten_hv,dtaux3d_hv,dtauy3d_hv + + real(kind=kind_phys),dimension(its:ite,kms:kme):: & + p3di_hv + +!----------------------------------------------------------------------------------------------------------------- + +! Outer j-loop. Allows consistency between WRF and MPAS in the driver. + + do j = jts,jte + + ! All variables for gwdo2d are tile-sized and have only a single + ! horizontal dimension. The _hv suffix refers to "horizontal vertical", + ! a reminder that there is a single horizontal index. Yes, we know that + ! variables that have only a horizontal index are not *really* _hv. + + ! All of the following 3d and 2d variables are declared intent(in) in the + ! gwdo2d subroutine, so there is no need to put the updated values back + ! from the temporary arrays back into the original arrays. + + ! Variables that are INTENT(IN) or INTENT(INOUT) + + ! 3d, interface levels: + do k = kts,kte+1 + do i = its,ite + p3di_hv(i,k) = p3di(i,k,j) enddo - if (kblk.ne.0) then -! -! compute flow-blocking stress -! - fbdcd = max(2.0-1.0/od(i),0.0) - taufb(i,kts) = 0.5*rhobar(i)*coefm(i)/dxmeter**2*fbdcd*dxyp(i) & - *olp(i)*zblk*ulow(i)**2 - tautem = taufb(i,kts)/real(kblk-kts) - do k = kts+1, kblk - taufb(i,k) = taufb(i,k-1) - tautem - enddo -! -! sum orographic GW stress and flow-blocking stress -! - taup(i,:) = taup(i,:) + taufb(i,:) - endif - endif - enddo -! -! calculate - (g)*d(tau)/d(pressure) and deceleration terms dtaux, dtauy -! - do k = kts,kte - do i = its,ite - taud(i,k) = 1. * (taup(i,k+1) - taup(i,k)) * g_ / del(i,k) - enddo - enddo -! -! if the gravity wave drag would force a critical line -! in the lower ksmm1 layers during the next deltim timestep, -! then only apply drag until that critical line is reached. -! - do k = kts,kpblmax-1 - do i = its,ite - if (k .le. kbl(i)) then - if (taud(i,k).ne.0.) & - dtfac(i) = min(dtfac(i),abs(velco(i,k)/(deltim*taud(i,k)))) - endif - enddo - enddo -! - do i = its,ite - dusfc(i) = 0. - dvsfc(i) = 0. - enddo -! - do k = kts,kte - do i = its,ite - taud(i,k) = taud(i,k) * dtfac(i) - dtaux = taud(i,k) * xn(i) - dtauy = taud(i,k) * yn(i) - dtaux2d(i,k) = dtaux - dtauy2d(i,k) = dtauy - dudt(i,k) = dtaux + dudt(i,k) - dvdt(i,k) = dtauy + dvdt(i,k) - dusfc(i) = dusfc(i) + dtaux * del(i,k) - dvsfc(i) = dvsfc(i) + dtauy * del(i,k) - enddo - enddo -! - do i = its,ite - dusfc(i) = (-1./g_) * dusfc(i) - dvsfc(i) = (-1./g_) * dvsfc(i) - enddo -! - return - end subroutine gwdo2d -!------------------------------------------------------------------------------- -!------------------------------------------------------------------------------- + enddo + + ! 3d, layers: + do k = kts,kte + do i = its,ite + rublten_hv(i,k) = rublten(i,k,j) + rvblten_hv(i,k) = rvblten(i,k,j) + 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) + p3d_hv(i,k) = p3d(i,k,j) + pi3d_hv(i,k) = pi3d(i,k,j) + z_hv(i,k) = z(i,k,j) + enddo + enddo + + ! 2d: + do i = its,ite + dx_hv(i) = dx(i,j) + var2d_hv(i) = var2d(i,j) + oc12d_hv(i) = oc12d(i,j) + sina_hv(i) = sina(i,j) + cosa_hv(i) = cosa(i,j) + oa2d1_hv(i) = oa2d1(i,j) + oa2d2_hv(i) = oa2d2(i,j) + oa2d3_hv(i) = oa2d3(i,j) + oa2d4_hv(i) = oa2d4(i,j) + ol2d1_hv(i) = ol2d1(i,j) + ol2d2_hv(i) = ol2d2(i,j) + ol2d3_hv(i) = ol2d3(i,j) + ol2d4_hv(i) = ol2d4(i,j) + enddo + + call bl_gwdo_run(sina=sina_hv,cosa=cosa_hv & + ,rublten=rublten_hv,rvblten=rvblten_hv & + ,dtaux3d=dtaux3d_hv,dtauy3d=dtauy3d_hv & + ,dusfcg=dusfcg_hv,dvsfcg=dvsfcg_hv & + ,uproj=u3d_hv,vproj=v3d_hv & + ,t1=t3d_hv,q1=qv3d_hv & + ,prsi=p3di_hv & + ,prsl=p3d_hv,prslk=pi3d_hv & + ,zl=z_hv & + ,var=var2d_hv,oc1=oc12d_hv & + ,oa2d1=oa2d1_hv, oa2d2=oa2d2_hv & + ,oa2d3=oa2d3_hv, oa2d4=oa2d4_hv & + ,ol2d1=ol2d1_hv, ol2d2=ol2d2_hv & + ,ol2d3=ol2d3_hv, ol2d4=ol2d4_hv & + ,g_=g,cp_=cp,rd_=rd,rv_=rv,fv_=ep1,pi_=pi & + ,dxmeter=dx_hv,deltim=dt & + ,its=its,ite=ite,kte=kte,kme=kte+1 & + ,errmsg=errmsg,errflg=errflg) + + ! Variables that are INTENT(OUT) or INTENT(INOUT): + + ! 3d, layers: + do k = kts,kte + do i = its,ite + rublten(i,k,j) = rublten_hv(i,k) + rvblten(i,k,j) = rvblten_hv(i,k) + dtaux3d(i,k,j) = dtaux3d_hv(i,k) + dtauy3d(i,k,j) = dtauy3d_hv(i,k) + enddo + enddo + + ! 2d: + do i = its,ite + dusfcg(i,j) = dusfcg_hv(i) + dvsfcg(i,j) = dvsfcg_hv(i) + enddo + + enddo ! Outer J-loop + + end subroutine gwdo + +!================================================================================================================= end module module_bl_gwdo +!================================================================================================================= diff --git a/phys/module_pbl_driver.F b/phys/module_pbl_driver.F index 3045adfe02..ea077aaf45 100644 --- a/phys/module_pbl_driver.F +++ b/phys/module_pbl_driver.F @@ -838,6 +838,10 @@ SUBROUTINE pbl_driver( & ! ! END FASDAS ! +! To accommodate shared physics + character*256 :: errmsg + integer :: errflg +! !------------------------------------------------------------------ ! !!!!!!!if using BEP set flag_bep to true @@ -2133,7 +2137,8 @@ SUBROUTINE pbl_driver( & ,ZNU=znu,ZNW=znw,P_TOP=p_top & ,CP=cp,G=g,RD=r_d & ,RV=r_v,EP1=ep_1,PI=3.141592653 & - ,DT=dtbl,DX=dx,KPBL2D=kpbl,ITIMESTEP=itimestep & + ,DT=dtbl,DX=dx2d,KPBL2D=kpbl,ITIMESTEP=itimestep & + ,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 ) diff --git a/phys/physics_mmm/bl_gwdo.F90 b/phys/physics_mmm/bl_gwdo.F90 new file mode 100644 index 0000000000..b314634539 --- /dev/null +++ b/phys/physics_mmm/bl_gwdo.F90 @@ -0,0 +1,649 @@ +!================================================================================================================= + module bl_gwdo + use ccpp_kind_types,only: kind_phys + + implicit none + private + public:: bl_gwdo_run, & + bl_gwdo_init, & + bl_gwdo_finalize + + + contains + + +!================================================================================================================= +!>\section arg_table_bl_gwdo_init +!!\html\include bl_gwdo_init.html +!! + subroutine bl_gwdo_init(errmsg,errflg) +!================================================================================================================= + +!--- output arguments: + character(len=*),intent(out):: errmsg + integer,intent(out):: errflg + +!----------------------------------------------------------------------------------------------------------------- + + errmsg = 'bl_gwdo_init OK' + errflg = 0 + + end subroutine bl_gwdo_init + +!================================================================================================================= +!>\section arg_table_bl_gwdo_finalize +!!\html\include bl_gwdo_finalize.html +!! + subroutine bl_gwdo_finalize(errmsg,errflg) +!================================================================================================================= + +!--- output arguments: + character(len=*),intent(out):: errmsg + integer,intent(out):: errflg + +!----------------------------------------------------------------------------------------------------------------- + + errmsg = 'bl_gwdo_finalize OK' + errflg = 0 + + end subroutine bl_gwdo_finalize + +!================================================================================================================= +!>\section arg_table_bl_gwdo_run +!!\html\include bl_gwdo_run.html +!! + subroutine bl_gwdo_run(sina, cosa, & + rublten,rvblten, & + dtaux3d,dtauy3d, & + dusfcg,dvsfcg, & + uproj, vproj, & + t1, q1, & + prsi, prsl, prslk, zl, & + var, oc1, & + oa2d1, oa2d2, & + oa2d3, oa2d4, & + ol2d1, ol2d2, & + ol2d3, ol2d4, & + g_, cp_, rd_, rv_, fv_, pi_, & + dxmeter, deltim, & + its, ite, kte, kme, & + errmsg, errflg ) +!------------------------------------------------------------------------------- +! +! abstract : +! this code handles the time tendencies of u v due to the effect of +! mountain induced gravity wave drag from sub-grid scale orography. +! this routine not only treats the traditional upper-level wave breaking due +! to mountain variance (alpert 1988), but also the enhanced +! lower-tropospheric wave breaking due to mountain convexity and asymmetry +! (kim and arakawa 1995). thus, in addition to the terrain height data +! in a model grid gox, additional 10-2d topographic statistics files are +! needed, including orographic standard deviation (var), convexity (oc1), +! asymmetry (oa4) and ol (ol4). these data sets are prepared based on the +! 30 sec usgs orography (hong 1999). the current scheme was implmented as in +! choi and hong (2015), which names kim gwdo since it was developed by +! kiaps staffs for kiaps integrated model system (kim). the scheme +! additionally includes the effects of orographic anisotropy and +! flow-blocking drag. +! coded by song-you hong and young-joon kim and implemented by song-you hong +! +! history log : +! 2015-07-01 hyun-joo choi add flow-blocking drag and orographic anisotropy +! +! references : +! choi and hong (2015), j. geophys. res. +! hong et al. (2008), wea. forecasting +! kim and doyle (2005), q. j. r. meteor. soc. +! kim and arakawa (1995), j. atmos. sci. +! alpet et al. (1988), NWP conference +! hong (1999), NCEP office note 424 +! +! input : +! dudt, dvdt - non-lin tendency for u and v wind component +! uproj, vproj - projection-relative U and V m/sec +! u1, v1 - zonal and meridional wind m/sec at t0-dt +! t1 - temperature deg k at t0-dt +! q1 - mixing ratio at t0-dt +! deltim - time step (s) +! del - positive increment of pressure across layer (pa) +! prslk, zl, prsl, prsi - pressure and height variables +! oa4, ol4, omax, var, oc1 - orographic statistics +! +! output : +! dudt, dvdt - wind tendency due to gwdo +! dtaux2d, dtauy2d - diagnoised orographic gwd +! dusfc, dvsfc - gw stress +! +!------------------------------------------------------------------------------- + implicit none +! + integer, parameter :: kts = 1 + integer , intent(in ) :: its, ite, kte, kme + real(kind=kind_phys) , intent(in ) :: g_, pi_, rd_, rv_, fv_,& + cp_, deltim + real(kind=kind_phys), dimension(its:) , intent(in ) :: dxmeter + real(kind=kind_phys), dimension(its:,:) , intent(inout) :: rublten, rvblten + real(kind=kind_phys), dimension(its:,:) , intent( out) :: dtaux3d, dtauy3d + real(kind=kind_phys), dimension(its:) , intent( out) :: dusfcg, dvsfcg + real(kind=kind_phys), dimension(its:) , intent(in ) :: sina, cosa + real(kind=kind_phys), dimension(its:,:) , intent(in ) :: uproj, vproj + real(kind=kind_phys), dimension(its:,:) , intent(in ) :: t1, q1, prslk, zl +! + real(kind=kind_phys), dimension(its:,:) , intent(in ) :: prsl + real(kind=kind_phys), dimension(its:,:) , intent(in ) :: prsi +! + real(kind=kind_phys), dimension(its:) , intent(in ) :: var, oc1, & + oa2d1, oa2d2, oa2d3, oa2d4, & + ol2d1, ol2d2, ol2d3, ol2d4 + character(len=*) , intent( out) :: errmsg + integer , intent( out) :: errflg +! + real(kind=kind_phys), parameter :: ric = 0.25 ! critical richardson number + real(kind=kind_phys), parameter :: dw2min = 1. + real(kind=kind_phys), parameter :: rimin = -100. + real(kind=kind_phys), parameter :: bnv2min = 1.0e-5 + real(kind=kind_phys), parameter :: efmin = 0.0 + real(kind=kind_phys), parameter :: efmax = 10.0 + real(kind=kind_phys), parameter :: xl = 4.0e4 + real(kind=kind_phys), parameter :: critac = 1.0e-5 + real(kind=kind_phys), parameter :: gmax = 1. + real(kind=kind_phys), parameter :: veleps = 1.0 + real(kind=kind_phys), parameter :: frc = 1.0 + real(kind=kind_phys), parameter :: ce = 0.8 + real(kind=kind_phys), parameter :: cg = 0.5 + integer,parameter :: kpblmin = 2 +! +! local variables +! + integer :: kpblmax + integer :: latd,lond + integer :: i,k,lcap,lcapp1,nwd,idir, & + klcap,kp1,ikount,kk +! + real(kind=kind_phys) :: fdir,cs,rcsks, & + wdir,ti,rdz,temp,tem2,dw2,shr2,bvf2,rdelks, & + wtkbj,tem,gfobnv,hd,fro,rim,temc,tem1,efact, & + temv,dtaux,dtauy +! + real(kind=kind_phys), dimension(its:ite,kts:kte) :: dudt, dvdt + real(kind=kind_phys), dimension(its:ite,kts:kte) :: dtaux2d, dtauy2d + real(kind=kind_phys), dimension(its:ite) :: dusfc, dvsfc + logical, dimension(its:ite) :: ldrag, icrilv, flag,kloop1 + real(kind=kind_phys), dimension(its:ite) :: coefm +! + real(kind=kind_phys), dimension(its:ite) :: taub, xn, yn, ubar, vbar, fr, & + ulow, rulow, bnv, oa, ol, rhobar, & + dtfac, brvf, xlinv, delks,delks1, & + zlowtop,cleff + real(kind=kind_phys), dimension(its:ite,kts:kte+1) :: taup + real(kind=kind_phys), dimension(its:ite,kts:kte-1) :: velco + real(kind=kind_phys), dimension(its:ite,kts:kte) :: bnv2, usqj, taud, rho, vtk, vtj + real(kind=kind_phys), dimension(its:ite,kts:kte) :: del + real(kind=kind_phys), dimension(its:ite,kts:kte) :: u1, v1 + real(kind=kind_phys), dimension(its:ite,4) :: oa4, ol4 +! + integer, dimension(its:ite) :: kbl, klowtop + integer, parameter :: mdir=8 + integer, dimension(mdir) :: nwdir + data nwdir/6,7,5,8,2,3,1,4/ +! +! variables for flow-blocking drag +! + real(kind=kind_phys), parameter :: frmax = 10. + real(kind=kind_phys), parameter :: olmin = 1.0e-5 + real(kind=kind_phys), parameter :: odmin = 0.1 + real(kind=kind_phys), parameter :: odmax = 10. +! + real(kind=kind_phys) :: fbdcd + real(kind=kind_phys) :: zblk, tautem + real(kind=kind_phys) :: fbdpe, fbdke + real(kind=kind_phys), dimension(its:ite) :: delx, dely + real(kind=kind_phys), dimension(its:ite,4) :: dxy4, dxy4p + real(kind=kind_phys), dimension(4) :: ol4p + real(kind=kind_phys), dimension(its:ite) :: dxy, dxyp, olp, od + real(kind=kind_phys), dimension(its:ite,kts:kte+1) :: taufb +! + integer, dimension(its:ite) :: komax + integer :: kblk +!------------------------------------------------------------------------------- +! +! constants +! + lcap = kte + lcapp1 = lcap + 1 + fdir = mdir / (2.0*pi_) +! +! initialize CCPP error flag and message +! + errmsg = '' + errflg = 0 +! +! calculate length of grid for flow-blocking drag +! + delx(its:ite) = dxmeter(its:ite) + dely(its:ite) = dxmeter(its:ite) + dxy4(its:ite,1) = delx(its:ite) + dxy4(its:ite,2) = dely(its:ite) + dxy4(its:ite,3) = sqrt(delx(its:ite)**2. + dely(its:ite)**2.) + dxy4(its:ite,4) = dxy4(its:ite,3) + dxy4p(its:ite,1) = dxy4(its:ite,2) + dxy4p(its:ite,2) = dxy4(its:ite,1) + dxy4p(its:ite,3) = dxy4(its:ite,4) + dxy4p(its:ite,4) = dxy4(its:ite,3) +! + cleff(its:ite) = dxmeter(its:ite) +! +! initialize arrays, array syntax is OK for OpenMP since these are local +! + ldrag = .false. ; icrilv = .false. ; flag = .true. +! + klowtop = 0 ; kbl = 0 +! + dtaux = 0. ; dtauy = 0. ; xn = 0. ; yn = 0. + ubar = 0. ; vbar = 0. ; rhobar = 0. ; ulow = 0. + oa = 0. ; ol = 0. ; taub = 0. +! + usqj = 0. ; bnv2 = 0. ; vtj = 0. ; vtk = 0. + taup = 0. ; taud = 0. ; dtaux2d = 0. ; dtauy2d = 0. +! + dtfac = 1.0 ; xlinv = 1.0/xl +! + komax = 0 + taufb = 0.0 +! + do k = kts,kte + do i = its,ite + vtj(i,k) = t1(i,k) * (1.+fv_*q1(i,k)) + vtk(i,k) = vtj(i,k) / prslk(i,k) + + ! Density (kg/m^3) + + rho(i,k) = 1./rd_ * prsl(i,k) / vtj(i,k) + + ! Delta p (positive) between interfaces levels (Pa) + + del(i,k) = prsi(i,k) - prsi(i,k+1) + + ! Earth-relative zonal and meridional winds (m/s) + + u1(i,k) = uproj(i,k)*cosa(i) - vproj(i,k)*sina(i) + v1(i,k) = uproj(i,k)*sina(i) + vproj(i,k)*cosa(i) + + enddo + enddo + +! + do i = its,ite + zlowtop(i) = 2. * var(i) + enddo +! + do i = its,ite + kloop1(i) = .true. + enddo +! + do k = kts+1,kte + do i = its,ite + if(zlowtop(i) .gt. 0.) then + if (kloop1(i).and.zl(i,k)-zl(i,1).ge.zlowtop(i)) then + klowtop(i) = k+1 + kloop1(i) = .false. + endif + endif + enddo + enddo +! + kpblmax = kte + do i = its,ite + kbl(i) = klowtop(i) + kbl(i) = max(min(kbl(i),kpblmax),kpblmin) + enddo +! +! determine the level of maximum orographic height +! + komax(:) = kbl(:) +! + do i = its,ite + delks(i) = 1.0 / (prsi(i,1) - prsi(i,kbl(i))) + delks1(i) = 1.0 / (prsl(i,1) - prsl(i,kbl(i))) + enddo +! +! compute low level averages within pbl +! + do k = kts,kpblmax + do i = its,ite + if (k.lt.kbl(i)) then + rcsks = del(i,k) * delks(i) + rdelks = del(i,k) * delks(i) + ubar(i) = ubar(i) + rcsks * u1(i,k) ! pbl u mean + vbar(i) = vbar(i) + rcsks * v1(i,k) ! pbl v mean + rhobar(i) = rhobar(i) + rdelks * rho(i,k) ! pbl rho mean + endif + enddo + enddo +! +! figure out low-level horizontal wind direction +! +! nwd 1 2 3 4 5 6 7 8 +! wd w s sw nw e n ne se +! + do i = its,ite + oa4(i,1) = oa2d1(i) + oa4(i,2) = oa2d2(i) + oa4(i,3) = oa2d3(i) + oa4(i,4) = oa2d4(i) + ol4(i,1) = ol2d1(i) + ol4(i,2) = ol2d2(i) + ol4(i,3) = ol2d3(i) + ol4(i,4) = ol2d4(i) + wdir = atan2(ubar(i),vbar(i)) + pi_ + idir = mod(nint(fdir*wdir),mdir) + 1 + nwd = nwdir(idir) + oa(i) = (1-2*int( (nwd-1)/4 )) * oa4(i,mod(nwd-1,4)+1) + ol(i) = ol4(i,mod(nwd-1,4)+1) +! +! compute orographic width along (ol) and perpendicular (olp) the wind direction +! + ol4p(1) = ol4(i,2) + ol4p(2) = ol4(i,1) + ol4p(3) = ol4(i,4) + ol4p(4) = ol4(i,3) + olp(i) = ol4p(mod(nwd-1,4)+1) +! +! compute orographic direction (horizontal orographic aspect ratio) +! + od(i) = olp(i)/max(ol(i),olmin) + od(i) = min(od(i),odmax) + od(i) = max(od(i),odmin) +! +! compute length of grid in the along(dxy) and cross(dxyp) wind directions +! + dxy(i) = dxy4(i,MOD(nwd-1,4)+1) + dxyp(i) = dxy4p(i,MOD(nwd-1,4)+1) + enddo +! +! saving richardson number in usqj for migwdi +! + do k = kts,kte-1 + do i = its,ite + ti = 2.0 / (t1(i,k)+t1(i,k+1)) + rdz = 1./(zl(i,k+1) - zl(i,k)) + tem1 = u1(i,k) - u1(i,k+1) + tem2 = v1(i,k) - v1(i,k+1) + dw2 = tem1*tem1 + tem2*tem2 + shr2 = max(dw2,dw2min) * rdz * rdz + bvf2 = g_*(g_/cp_+rdz*(vtj(i,k+1)-vtj(i,k))) * ti + usqj(i,k) = max(bvf2/shr2,rimin) + bnv2(i,k) = 2.0*g_*rdz*(vtk(i,k+1)-vtk(i,k))/(vtk(i,k+1)+vtk(i,k)) + enddo + enddo +! +! compute the "low level" or 1/3 wind magnitude (m/s) +! + do i = its,ite + ulow(i) = max(sqrt(ubar(i)*ubar(i) + vbar(i)*vbar(i)), 1.0) + rulow(i) = 1./ulow(i) + enddo +! + do k = kts,kte-1 + do i = its,ite + velco(i,k) = 0.5 * ((u1(i,k)+u1(i,k+1)) * ubar(i) & + + (v1(i,k)+v1(i,k+1)) * vbar(i)) + velco(i,k) = velco(i,k) * rulow(i) + if ((velco(i,k).lt.veleps) .and. (velco(i,k).gt.0.)) then + velco(i,k) = veleps + endif + enddo + enddo +! +! no drag when critical level in the base layer +! + do i = its,ite + ldrag(i) = velco(i,1).le.0. + enddo +! +! no drag when velco.lt.0 +! + do k = kpblmin,kpblmax + do i = its,ite + if (k .lt. kbl(i)) ldrag(i) = ldrag(i).or. velco(i,k).le.0. + enddo + enddo +! +! the low level weighted average ri is stored in usqj(1,1; im) +! the low level weighted average n**2 is stored in bnv2(1,1; im) +! this is called bnvl2 in phy_gwd_alpert_sub not bnv2 +! rdelks (del(k)/delks) vert ave factor so we can * instead of / +! + do i = its,ite + wtkbj = (prsl(i,1)-prsl(i,2)) * delks1(i) + bnv2(i,1) = wtkbj * bnv2(i,1) + usqj(i,1) = wtkbj * usqj(i,1) + enddo +! + do k = kpblmin,kpblmax + do i = its,ite + if (k .lt. kbl(i)) then + rdelks = (prsl(i,k)-prsl(i,k+1)) * delks1(i) + bnv2(i,1) = bnv2(i,1) + bnv2(i,k) * rdelks + usqj(i,1) = usqj(i,1) + usqj(i,k) * rdelks + endif + enddo + enddo +! + do i = its,ite + ldrag(i) = ldrag(i) .or. bnv2(i,1).le.0.0 + ldrag(i) = ldrag(i) .or. ulow(i).eq.1.0 + ldrag(i) = ldrag(i) .or. var(i) .le. 0.0 + enddo +! +! set all ri low level values to the low level value +! + do k = kpblmin,kpblmax + do i = its,ite + if (k .lt. kbl(i)) usqj(i,k) = usqj(i,1) + enddo + enddo +! + do i = its,ite + if (.not.ldrag(i)) then + bnv(i) = sqrt( bnv2(i,1) ) + fr(i) = bnv(i) * rulow(i) * var(i) * od(i) + fr(i) = min(fr(i),frmax) + xn(i) = ubar(i) * rulow(i) + yn(i) = vbar(i) * rulow(i) + endif + enddo +! +! compute the base level stress and store it in taub +! calculate enhancement factor, number of mountains & aspect +! ratio const. use simplified relationship between standard +! deviation & critical hgt +! + do i = its,ite + if (.not. ldrag(i)) then + efact = (oa(i) + 2.) ** (ce*fr(i)/frc) + efact = min( max(efact,efmin), efmax ) + coefm(i) = (1. + ol(i)) ** (oa(i)+1.) + xlinv(i) = coefm(i) / cleff(i) + tem = fr(i) * fr(i) * oc1(i) + gfobnv = gmax * tem / ((tem + cg)*bnv(i)) + taub(i) = xlinv(i) * rhobar(i) * ulow(i) * ulow(i) & + * ulow(i) * gfobnv * efact + else + taub(i) = 0.0 + xn(i) = 0.0 + yn(i) = 0.0 + endif + enddo +! +! now compute vertical structure of the stress. +! + do k = kts,kpblmax + do i = its,ite + if (k .le. kbl(i)) taup(i,k) = taub(i) + enddo + enddo +! + do k = kpblmin, kte-1 ! vertical level k loop! + kp1 = k + 1 + do i = its,ite +! +! unstablelayer if ri < ric +! unstable layer if upper air vel comp along surf vel <=0 (crit lay) +! at (u-c)=0. crit layer exists and bit vector should be set (.le.) +! + if (k .ge. kbl(i)) then + icrilv(i) = icrilv(i) .or. ( usqj(i,k) .lt. ric) & + .or. (velco(i,k) .le. 0.0) + brvf(i) = max(bnv2(i,k),bnv2min) ! brunt-vaisala frequency squared + brvf(i) = sqrt(brvf(i)) ! brunt-vaisala frequency + endif + enddo +! + do i = its,ite + if (k .ge. kbl(i) .and. (.not. ldrag(i))) then + if (.not.icrilv(i) .and. taup(i,k) .gt. 0.0 ) then + temv = 1.0 / velco(i,k) + tem1 = coefm(i)/dxy(i)*(rho(i,kp1)+rho(i,k))*brvf(i)*velco(i,k)*0.5 + hd = sqrt(taup(i,k) / tem1) + fro = brvf(i) * hd * temv +! +! rim is the minimum-richardson number by shutts (1985) +! + tem2 = sqrt(usqj(i,k)) + tem = 1. + tem2 * fro + rim = usqj(i,k) * (1.-fro) / (tem * tem) +! +! check stability to employ the 'saturation hypothesis' +! of lindzen (1981) except at tropospheric downstream regions +! + if (rim .le. ric) then ! saturation hypothesis! + if ((oa(i) .le. 0.).or.(kp1 .ge. kpblmin )) then + temc = 2.0 + 1.0 / tem2 + hd = velco(i,k) * (2.*sqrt(temc)-temc) / brvf(i) + taup(i,kp1) = tem1 * hd * hd + endif + else ! no wavebreaking! + taup(i,kp1) = taup(i,k) + endif + endif + endif + enddo + enddo +! + if (lcap.lt.kte) then + do klcap = lcapp1,kte + do i = its,ite + taup(i,klcap) = prsi(i,klcap) / prsi(i,lcap) * taup(i,lcap) + enddo + enddo + endif + do i = its,ite + if (.not.ldrag(i)) then +! +! determine the height of flow-blocking layer +! + kblk = 0 + fbdpe = 0.0 + fbdke = 0.0 + do k = kte, kpblmin, -1 + if (kblk.eq.0 .and. k.le.kbl(i)) then + fbdpe = fbdpe + bnv2(i,k)*(zl(i,kbl(i))-zl(i,k)) & + *del(i,k)/g_/rho(i,k) + fbdke = 0.5*(u1(i,k)**2.+v1(i,k)**2.) +! +! apply flow-blocking drag when fbdpe >= fbdke +! + if (fbdpe.ge.fbdke) then + kblk = k + kblk = min(kblk,kbl(i)) + zblk = zl(i,kblk)-zl(i,kts) + endif + endif + enddo + if (kblk.ne.0) then +! +! compute flow-blocking stress +! + fbdcd = max(2.0-1.0/od(i),0.0) + taufb(i,kts) = 0.5*rhobar(i)*coefm(i)/dxmeter(i)**2*fbdcd*dxyp(i) & + *olp(i)*zblk*ulow(i)**2 + tautem = taufb(i,kts)/real(kblk-kts) + do k = kts+1, kblk + taufb(i,k) = taufb(i,k-1) - tautem + enddo +! +! sum orographic GW stress and flow-blocking stress +! + taup(i,:) = taup(i,:) + taufb(i,:) + endif + endif + enddo +! +! calculate - (g)*d(tau)/d(pressure) and deceleration terms dtaux, dtauy +! + do k = kts,kte + do i = its,ite + taud(i,k) = 1. * (taup(i,k+1) - taup(i,k)) * g_ / del(i,k) + enddo + enddo +! +! if the gravity wave drag would force a critical line +! in the lower ksmm1 layers during the next deltim timestep, +! then only apply drag until that critical line is reached. +! + do k = kts,kpblmax-1 + do i = its,ite + if (k .le. kbl(i)) then + if (taud(i,k).ne.0.) & + dtfac(i) = min(dtfac(i),abs(velco(i,k)/(deltim*taud(i,k)))) + endif + enddo + enddo +! + do i = its,ite + dusfc(i) = 0. + dvsfc(i) = 0. + enddo +! + do k = kts,kte + do i = its,ite + taud(i,k) = taud(i,k) * dtfac(i) + dtaux = taud(i,k) * xn(i) + dtauy = taud(i,k) * yn(i) + dtaux2d(i,k) = dtaux + dtauy2d(i,k) = dtauy + dudt(i,k) = dtaux + dvdt(i,k) = dtauy + dusfc(i) = dusfc(i) + dtaux * del(i,k) + dvsfc(i) = dvsfc(i) + dtauy * del(i,k) + enddo + enddo +! + do i = its,ite + dusfc(i) = (-1./g_) * dusfc(i) + dvsfc(i) = (-1./g_) * dvsfc(i) + enddo +! +! rotate tendencies from zonal/meridional back to model grid +! + do k = kts,kte + do i = its,ite + rublten(i,k) = rublten(i,k)+dudt(i,k)*cosa(i) + dvdt(i,k)*sina(i) + rvblten(i,k) = rvblten(i,k)-dudt(i,k)*sina(i) + dvdt(i,k)*cosa(i) + dtaux3d(i,k) = dtaux2d(i,k)*cosa(i) + dtauy2d(i,k)*sina(i) + dtauy3d(i,k) =-dtaux2d(i,k)*sina(i) + dtauy2d(i,k)*cosa(i) + enddo + enddo + do i = its,ite + dusfcg(i) = dusfc(i)*cosa(i) + dvsfc(i)*sina(i) + dvsfcg(i) =-dusfc(i)*sina(i) + dvsfc(i)*cosa(i) + enddo + return + end subroutine bl_gwdo_run + + +!================================================================================================================= + end module bl_gwdo +!================================================================================================================= + diff --git a/wrftladj/module_pbl_driver_ad.F b/wrftladj/module_pbl_driver_ad.F index 27fc22efbe..3001a38490 100644 --- a/wrftladj/module_pbl_driver_ad.F +++ b/wrftladj/module_pbl_driver_ad.F @@ -502,6 +502,10 @@ SUBROUTINE A_PBL_DRIVER(itimestep, dt, u_frame, v_frame, bldt, curr_secs& REAL :: seamask, thsk, zzz, unew, vnew, tnew, qnew, umom, vmom REAL :: z0, z1, z2, w1, w2 !------------------------------------------------------------------ +! For shared physics + REAL, DIMENSION(ims:ime, jms:jme) :: dx2dtmp + character*256 :: errmsg + integer :: errflg ! !!!!!!!if using BEP set flag_bep to true INTEGER :: branch @@ -635,6 +639,7 @@ SUBROUTINE A_PBL_DRIVER(itimestep, dt, u_frame, v_frame, bldt, curr_secs& ELSE CALL PUSHCONTROL1B(1) END IF + dx2dtmp(i,j)=dx END DO CALL PUSHINTEGER4(i - 1) CALL PUSHINTEGER4(ad_from) @@ -677,8 +682,9 @@ SUBROUTINE A_PBL_DRIVER(itimestep, dt, u_frame, v_frame, bldt, curr_secs& & oc12d=oc12d, oa2d1=oa1, oa2d2=oa2, oa2d3=oa3, oa2d4=oa4& & , ol2d1=ol1, ol2d2=ol2, ol2d3=ol3, ol2d4=ol4, & & SINA=sina,COSA=cosa, znu=znu, & +& errmsg= errmsg, errflg=errflg, & & znw=znw, p_top=p_top, cp=cp, g=g, rd=r_d, rv=& -& r_v, ep1=ep_1, pi=3.141592653, dt=dtbl, dx=dx, kpbl2d=& +& r_v, ep1=ep_1, pi=3.141592653, dt=dtbl, dx=dx2dtmp, kpbl2d=& & kpbl, itimestep=itimestep, 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, &