source: trunk/LMDZ.MARS/libf/phymars/orodrag_mod.F90 @ 3118

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