MODULE orodrag_mod IMPLICIT NONE CONTAINS SUBROUTINE ORODRAG( ngrid,nlayer, kgwd, kgwdim, kdx, ktest, ptimestep, & pplev,pplay,zgeom,pt,pu, pv, pvar, psig, pgam, pthe, & !OUTPUTS PULOW,PVLOW, pdudt,pdvdt,pdtdt) !-------------------------------------------------------------------------------- ! The main routine ORODRAG that does the orographic gravity wave parameterization. ! This routine computes the physical tendencies of the prognostic variables U,V, and ! T due to vertical transport by subgridscale orographically excited gravity waves ! M.MILLER + B.RITTER E.C.M.W.F. 15/06/86. ! F.LOTT + M. MILLER E.C.M.W.F. 22/11/94 !-- ! The purpose of this routine is: ! 1) call OROSETUP to get the parameters such as the top altitude(level) of the blocked ! flow and other critical levels. ! 2) call GWSTRESS and GWPROFILE to get the gw stress. ! 3) calculate the zonal/meridional wind tendency based on the GW stress and mountain drag ! on the flow. ! 4) calculate the temperature tendency based on the differtial of square of the wind ! between t and t+dt. ! Rewrited by J.LIU LMDZ.COMMON/phymars/libf 29/01/2022 !-- ! REFERENCE. ! 1. SEE ECMWF RESEARCH DEPARTMENT DOCUMENTATION OF THE "I.F.S." ! 2. Lott, F., & Miller, M. J. (1997). A new subgrid‐scale orographic drag parametrization: ! Its formulation and testing.Quarterly Journal of the Royal cMeteorological Society, ! 123(537), 101-127. !----------------------------------------------------------------------- use dimradmars_mod, only: ndomainsz USE gwstress_mod, ONLY: gwstress USE gwprofil_mod, ONLY: gwprofil USE comcstfi_h, ONLY: g, cpp USE yoegwd_h, ONLY: gkwake implicit none ! 0. DECLARATIONS: ! 0.1 INPUTS INTEGER, intent(in):: ngrid ! Number of atmospheric columns INTEGER, intent(in):: nlayer ! Number of atmospheric layers INTEGER, intent(in):: kgwd ! Number of points at which the scheme is called INTEGER, intent(in):: kgwdim ! kgwdim=max(1,kgwd) INTEGER, intent(in):: kdx(ndomainsz) ! Points at which to call the scheme. INTEGER, intent(in):: ktest(ndomainsz) ! Map of calling points REAL, intent(in):: pu(ndomainsz,nlayer) ! Zonal wind at full levels(m/s) (has been inversed by DRAG_NORO, =inv_pu) REAL, intent(in):: pv(ndomainsz,nlayer) ! Meridional winds at full levels(m/s)(has been inversed by DRAG_NORO, =inv_pv) REAL, intent(in):: pt(ndomainsz,nlayer) ! Temperature at full levels(m/s) (has been inversed by DRAG_NORO=inv_pt) REAL, intent(in):: pvar(ndomainsz) ! Sub-grid scale standard deviation REAL, intent(in):: psig(ndomainsz) ! Sub-grid scale slope REAL, intent(in):: pgam(ndomainsz) ! Sub-grid scale anisotropy REAL, intent(in):: pthe(ndomainsz) ! Sub-grid scale principal axes angle REAL, intent(in):: zgeom(ndomainsz,nlayer) ! Geopotential height of full levels REAL, intent(in):: pplay(ndomainsz,nlayer) ! Pressure at full levels(Pa) (has been inversed by DRAG_NORO=inv_pplay) REAL, intent(in):: pplev(ndomainsz,nlayer+1) ! Pressure at 1/2 levels(Pa) (has been inversed by DRAG_NORO=inv_pplev) REAL, intent(in):: ptimestep ! Time step of the Physics(s) ! 0.2 OUTPUTS REAL, intent(out):: pdtdt(ndomainsz,nlayer) ! Tendency on temperature (K/s) due to orographic gravity waves REAL, intent(out):: pdvdt(ndomainsz,nlayer) ! Tendency on meridional wind (m/s/s) due to orographic gravity waves REAL, intent(out):: pdudt(ndomainsz,nlayer) ! Tendency on zonal wind (m/s/s) due to orographic gravity waves REAL, intent(out):: PULOW(ndomainsz) ! Low level zonal wind REAL, intent(out):: PVLOW(ndomainsz) ! Low level meridional wind !0.3 LOCAL ARRAYS ! INTEGER ISECT(ndomainsz) ! not used ? INTEGER ICRIT(ndomainsz) ! Critical layer where orographic GW breaks INTEGER IKCRITH(ndomainsz) ! Dynamical mixing height for the breaking of gravity waves INTEGER IKenvh(ndomainsz) ! Top of the blocked layer INTEGER IKNU(ndomainsz) ! 4*pvar layer INTEGER IKNU2(ndomainsz) ! 3*pvar layer INTEGER IKCRIT(ndomainsz) ! Top of low level flow height ! INTEGER IKHLIM(ndomainsz) ! not used? integer ji,jk,jl integer ilevp1 !=nlayer+1 just used once. Maybe directly use nlayer+1 to replace in the future ! real ztauf(ndomainsz,nlayer+1) ! not used real zdelp ! =pplev(nlayer+1)-pplev(nlayer) dp of two 1/2 levels real zb ! B(pgam), is a seceond order polynomial fit of SS ANISOTRPY pgam(Reference 2) real zc ! C(pgam), is a seceond order polynomial fit of SS ANISOTRPY pgam(Reference 2) real zbet ! Equation (16) blocked flow drag real zconb ! Middle part of the equation (16) real zabsv ! Tail of equation (16),U*|U|/2+V*|V|/2 real zzd1 ! Second tail part of the equation (16) real ratio ! Aspect ratio of the obstacle (mountain),equation (14) real zust ! U(t+dt) real zvst ! V(t+dt) real zdis ! =1/2*[U(t)*|U(t)|+V(t)*|V(t)|-U(t+dt)*|U(t+dt)|-V(t+dt)*|V(t+dt)|] real ztemp ! T=-g*dtao/U_fistlevel*dp temperature due to GW stress? REAL ZTAU(ndomainsz,nlayer+1) ! Output of subrountine GWPROFILE: Gravtiy wave stress REAL BV(ndomainsz,nlayer+1) ! Brunt–Väisälä frequency at 1/2 level (output from OROSETUP input for GWSTRESS and GWPROFIL) REAL ZVPH(ndomainsz,nlayer+1) ! Low level wind speed U_H in Ref.2 REAL ZRHO(ndomainsz,nlayer+1) ! density calculated (output) by the OROSETUP rountines(but not used by ORODRAG) ! used as input for GWSTRESS and GWPROFIL REAL PRI(ndomainsz,nlayer+1) ! Mean flow richardson number at 1/2 level REAL ZpsI(ndomainsz,nlayer+1) ! The angle between the incident flow direction and the normal ridge direction of the mountain REAL Zzdep(ndomainsz,nlayer) ! dp by full level REAL ZDUDT(ndomainsz) ! du/dt due to oro-gw REAL ZDVDT(ndomainsz) ! dv/dt due to oro-gw REAL ZDTDT(ndomainsz) ! dT/dt due to oro-gw REAL ZDEDT(ndomainsz) ! zdis/ptimestemp, zdedt/Cp=dT/dt REAL ZVIDIS(ndomainsz) ! in an invalid line, not used REAL Znu(ndomainsz) ! A critical value see equation 9(Since it only compute inside OROSETUP,Maybe needs delete in future) REAL Zd1(ndomainsz) ! Bcos^2(psi)-Csin^2(psi) see equation 17 or 18 in Ref.2 REAL Zd2(ndomainsz) ! (B-C)sin(psi)cos(psi) see equation 17 or 18 in Ref.2 REAL Zdmod(ndomainsz) ! sqrt(zd1^2+zd2^2) !------------------------------------------------------------------------------------ ! 1.INITIALIZATION (NOT USED ) !------------------------------------------------------------------------------------ ! 100 CONTINUE ! continue tag without source, maybe need delete in future !* 1.1 Computional constants ! 110 CONTINUE ! continue tag without source, maybe need delete in future ! kfdia=ndomainsz ! ZTMST=TWODT ! IF(NSTEP.EQ.NSTART) ZTMST=0.5*TWODT ! nlayerM1=nlayer-1 ! ZTMST=ptimestep ! ZRTMST=1./ptimestep ! 120 CONTINUE ! continue tag without source, maybe need delete in future !* 1.3 Check whether row contains point for printing ! 130 CONTINUE ! continue tag without source, maybe need delete in future !------------------------------------------------------------------------------------ ! 2. Precompute basic state variables . !------------------------------------------------------------------------------------ ! 200 CONTINUE ! continue tag without source, maybe need delete in future ! Define low level wind,project winds in plane of low level wind, determine sector ! in which to take the variance and set indicator for critical levels CALL OROSETUP( ngrid, nlayer, ktest, pplev, pplay, pu, pv, pt, zgeom, & pvar,pthe, pgam, & !output in capital IKCRIT, IKCRITH, ICRIT, IKENVH,IKNU,IKNU2, & !ISECT, IKHLIM, not used ZRHO,PRI,BV,ZTAU,ZVPH,ZPSI,ZZDEP,ZNU,ZD1,ZD2,ZDMOD, & PULOW, PVLOW) !------------------------------------------------------------------------------------ ! 3. Computes low level stresses using subcritical and super critical forms. ! Computes anisotropy coefficient as measure of orographic two-dimensionality !------------------------------------------------------------------------------------ ! 300 CONTINUE ! continue tag without source, maybe need delete in future ! Computes low level stresses CALL GWSTRESS(ngrid,nlayer,ktest,zrho, BV, pvar,psig,zgeom,zdmod,& ! Notice that this 3 variables are actually not used due to illness lines ICRIT,IKNU,ZVPH, & ! not used variables ! IKCRIT,ISECT,IKHLIM, IKCRITH,IKENVH,PVAR1,pgam,zd1,zd2,znu, & ! not defined not used variables ! ZTFR !in(as 0.0)-output: ZTAU ) !------------------------------------------------------------------------------------ ! 4. Compute stress profile. !------------------------------------------------------------------------------------ ! 400 CONTINUE ! continue tag without source, maybe need delete in future ! Compute stress profile. CALL GWPROFIL( ngrid, nlayer, kgwd ,kdx, ktest, & IKCRITH, ICRIT, IKENVH, IKNU,IKNU2,pplev, & ZRHO, BV, ZVPH, PRI,zdmod, psig , pvar, & !not used variables !IKCRIT,ZTAUf,ZTFR,znu,pgam, ! in-output ZTAU ) !------------------------------------------------------------------------------------ ! 5. Compute tendencies. !------------------------------------------------------------------------------------ ! 500 CONTINUE ! continue tag without source, maybe need delete in future ! EXplicit solution at all levels for gravity wave ! Implicit solution for the blocked levels ! DO JL=KIDIA,KFDIA ! Initialization the tendencies DO JL=1,ndomainsz ZVIDIS(JL)=0.0 !An invalid lines contain this variable ZDUDT(JL)=0.0 ZDVDT(JL)=0.0 ZDTDT(JL)=0.0 ENDDO ! all the calculation of equation (16-18) are in the flowing loop ! 1) if the wind flow over the mountain [Either the mountain is too low or the ! flow higher than the blocked level], then calculate the wind tendencies by oro-GW stress directly ! 2) if the flow level lower than the blocked level, then calculate the wind tendencies by obstacle drag ! ILEVP1=nlayer+1 DO JK=1,nlayer DO JL=1,kgwd ! Pick up the position of each calling points JI=kdx(JL) ! dp ZDELP=pplev(Ji,JK+1)-pplev(Ji,JK) ! T=-g*dtao/U_fistlevel*dp temperature due to GW stress? ZTEMP=-g*(ZTAU(Ji,JK+1)-ZTAU(Ji,JK))/(ZVPH(Ji,ILEVP1)*ZDELP) ! du/dt=(U_low*tao1-V_low*tao2)*T/(tao1^2+tao2^2)*0.5 ZDUDT(JI)=(PULOW(JI)*Zd1(ji)-pvlow(ji)*zd2(ji))*ztemp/zdmod(ji) ! dv/dt=(V_low*tao1+U_low*tao2)*T/(tao1^2+tao2^2)*0.5 ZDVDT(JI)=(PVLOW(JI)*Zd1(ji)+pulow(ji)*zd2(ji))*ztemp/zdmod(ji) if(jk.ge.ikenvh(ji)) then ! IF the level lower than(notice here the level is inversed thus it is lower) ! the top of the blocked layer (commonlly is the highest top of a mountain.) ! Then use equation (16) to calculate the mountain drag on the incident flow ! Coefficient B and C Ref.2 zb=1.0-0.18*pgam(ji)-0.04*pgam(ji)**2 zc=0.48*pgam(ji)+0.3*pgam(ji)**2 ! Middle part of the equation (16) zconb=2.*ptimestep*GKWAKE*psig(ji)/(4.*pvar(ji)) ! Tail part of the equation (16),U*|U|/2+V*|V|/2 zabsv=sqrt(pu(JI,JK)**2+pv(JI,JK)**2)/2. ! Second tail part of the equation (16) zzd1=zb*cos(zpsi(ji,jk))**2+zc*sin(zpsi(ji,jk))**2 ! Aspect ratio of the obstacle(topography), equation(14) ratio=(cos(zpsi(ji,jk))**2+pgam(ji)*sin(zpsi(ji,jk))**2)/(pgam(ji)*cos(zpsi(ji,jk))**2+sin(zpsi(ji,jk))**2) ! Equation (16) zbet=max(0.,2.-1./ratio)*zconb*zzdep(ji,jk)*zzd1*zabsv ! du/dt and dv/dt due to the mountain drag zdudt(ji)=-pu(ji,jk)/ptimestep zdvdt(ji)=-pv(ji,jk)/ptimestep zdudt(ji)=zdudt(ji)*(zbet/(1.+zbet)) zdvdt(ji)=zdvdt(ji)*(zbet/(1.+zbet)) end if ! wind tendencies due to the orographic gravity wave (due to stress(flow over) and drag (blocked)) pdudt(JI,JK)=ZDUDT(JI) pdvdt(JI,JK)=ZDVDT(JI) ! winds at t+dt (oro-gw induced increaments are added) ZUST=pu(JI,JK)+ptimestep*ZDUDT(JI) ZVST=pv(JI,JK)+ptimestep*ZDVDT(JI) ! =1/2*[U(t)*|U(t)|+V(t)*|V(t)|-U(t+dt)*|U(t+dt)|-V(t+dt)*|V(t+dt)|] ZDIS=0.5*(pu(JI,JK)**2+pv(JI,JK)**2-ZUST**2-ZVST**2) ZDEDT(JI)=ZDIS/ptimestep ! This line no use maybe delete in the future ZVIDIS(JI)=ZVIDIS(JI)+ZDIS*ZDELP ! Temperature tendencies due to the orographic gravity wave (due to stress(flow over) and drag (blocked)) ZDTDT(JI)=ZDEDT(JI)/cpp pdtdt(JI,JK)=ZDTDT(JI) ENDDO !JL=1,kgwd end DO !JK=1,nlayer END SUBROUTINE ORODRAG END MODULE orodrag_mod