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

Last change on this file since 3026 was 2651, checked in by emillour, 3 years ago

Mars GCM:

  • turn sugwd.F to sugwd.F90 with extra comments
  • turn yoegwd.h into a module

JL+EM
File size: 14.2 KB
RevLine 
[2642]1MODULE orodrag_mod
[1912]2     
[2642]3IMPLICIT NONE
[1912]4     
[2642]5CONTAINS
[1912]6     
[2642]7      SUBROUTINE ORODRAG( ngrid,nlayer, kgwd, kgwdim, kdx, ktest, ptimestep, &
8                  pplev,pplay,zgeom,pt,pu, pv, pvar, psig, pgam, pthe,       &
9                  !OUTPUTS
10                  PULOW,PVLOW, pdudt,pdvdt,pdtdt)
11                 
12     !--------------------------------------------------------------------------------
13     ! The main routine ORODRAG that does the orographic gravity wave parameterization.
14     ! This routine computes the physical tendencies of the prognostic variables U,V, and
15     ! T due to vertical transport by subgridscale orographically excited gravity waves
16     ! M.MILLER + B.RITTER   E.C.M.W.F.     15/06/86.
17     ! F.LOTT + M. MILLER    E.C.M.W.F.     22/11/94
18     !--
19     ! The purpose of this routine is:
20     ! 1) call OROSETUP to get the parameters such as the top altitude(level) of the blocked
21     !    flow and other critical levels.
22     ! 2) call GWSTRESS and GWPROFILE to get the gw stress.
23     ! 3) calculate the zonal/meridional wind tendency based on the GW stress and mountain drag
24     !    on the flow.
25     ! 4) calculate the temperature tendency based on the differtial of square of the wind
26     !    between t and t+dt.
27     ! Rewrited by J.LIU     LMDZ.COMMON/phymars/libf  29/01/2022
28     !--
29     ! REFERENCE.
30     ! 1. SEE ECMWF RESEARCH DEPARTMENT DOCUMENTATION OF THE "I.F.S."
31     ! 2. Lott, F., & Miller, M. J. (1997). A new subgrid‐scale orographic drag parametrization:
32     !    Its formulation and testing.Quarterly Journal of the Royal cMeteorological Society,
33     !    123(537), 101-127.
34     !-----------------------------------------------------------------------
35     
36      use dimradmars_mod, only: ndomainsz
[1912]37      USE gwstress_mod, ONLY: gwstress
[1913]38      USE gwprofil_mod, ONLY: gwprofil
[1912]39      USE comcstfi_h, ONLY: g, cpp
[2651]40      USE yoegwd_h, ONLY: gkwake
[2642]41     
[38]42      implicit none
[2616]43     
[2642]44      ! 0. DECLARATIONS:
45     
46      ! 0.1   INPUTS
47      INTEGER, intent(in):: ngrid          ! Number of atmospheric columns
48      INTEGER, intent(in):: nlayer         ! Number of atmospheric layers
49      INTEGER, intent(in):: kgwd           ! Number of points at which the scheme is called
50      INTEGER, intent(in):: kgwdim         ! kgwdim=max(1,kgwd)
51      INTEGER, intent(in):: kdx(ndomainsz)     ! Points at which to call the scheme.
52      INTEGER, intent(in):: ktest(ndomainsz)   ! Map of calling points   
[38]53
[2642]54      REAL, intent(in):: pu(ndomainsz,nlayer)  ! Zonal wind at full levels(m/s) (has been inversed by DRAG_NORO, =inv_pu)
55      REAL, intent(in):: pv(ndomainsz,nlayer)  ! Meridional winds at full levels(m/s)(has been inversed by DRAG_NORO, =inv_pv)
56      REAL, intent(in):: pt(ndomainsz,nlayer)  ! Temperature at full levels(m/s) (has been inversed by DRAG_NORO=inv_pt)
57      REAL, intent(in):: pvar(ndomainsz)       ! Sub-grid scale standard deviation
58      REAL, intent(in):: psig(ndomainsz)       ! Sub-grid scale slope
59      REAL, intent(in):: pgam(ndomainsz)       ! Sub-grid scale anisotropy
60      REAL, intent(in):: pthe(ndomainsz)       ! Sub-grid scale principal axes angle
61      REAL, intent(in):: zgeom(ndomainsz,nlayer)     ! Geopotential height of full levels
62      REAL, intent(in):: pplay(ndomainsz,nlayer)     ! Pressure at full levels(Pa) (has been inversed by DRAG_NORO=inv_pplay)
63      REAL, intent(in):: pplev(ndomainsz,nlayer+1)   ! Pressure at 1/2 levels(Pa) (has been inversed by DRAG_NORO=inv_pplev) 
64      REAL, intent(in):: ptimestep                   ! Time step of the Physics(s)
[1047]65     
[2642]66      ! 0.2   OUTPUTS
67      REAL, intent(out):: pdtdt(ndomainsz,nlayer)   ! Tendency on temperature (K/s) due to orographic gravity waves
68      REAL, intent(out):: pdvdt(ndomainsz,nlayer)   ! Tendency on meridional wind (m/s/s) due to orographic gravity waves
69      REAL, intent(out):: pdudt(ndomainsz,nlayer)   ! Tendency on zonal wind (m/s/s) due to orographic gravity waves
70      REAL, intent(out):: PULOW(ndomainsz)          ! Low level zonal wind
71      REAL, intent(out):: PVLOW(ndomainsz)          ! Low level meridional wind
[38]72
[2642]73      !0.3   LOCAL ARRAYS             
74!      INTEGER ISECT(ndomainsz)      ! not used ?
75      INTEGER ICRIT(ndomainsz)       ! Critical layer where orographic GW breaks
76      INTEGER IKCRITH(ndomainsz)     ! Dynamical mixing height for the breaking of gravity waves
77      INTEGER IKenvh(ndomainsz)      ! Top of the  blocked layer
78      INTEGER IKNU(ndomainsz)        ! 4*pvar layer
79      INTEGER IKNU2(ndomainsz)       ! 3*pvar layer
80      INTEGER IKCRIT(ndomainsz)      ! Top of low level flow height
81!      INTEGER IKHLIM(ndomainsz)     ! not used?
82     
83      integer ji,jk,jl
84      integer ilevp1                    !=nlayer+1 just used once. Maybe directly use nlayer+1 to replace in the future
85     
86!      real ztauf(ndomainsz,nlayer+1)   ! not used
87      real zdelp                        ! =pplev(nlayer+1)-pplev(nlayer) dp of two 1/2 levels
88      real zb                           ! B(pgam), is a seceond order polynomial fit of SS ANISOTRPY pgam(Reference 2)
89      real zc                           ! C(pgam), is a seceond order polynomial fit of SS ANISOTRPY pgam(Reference 2)
90      real zbet                         ! Equation (16) blocked flow drag
91      real zconb                        ! Middle part of the equation (16)
92      real zabsv                        ! Tail of equation (16),U*|U|/2+V*|V|/2
93      real zzd1                         ! Second tail part of the equation (16)
94      real ratio                        ! Aspect ratio of the obstacle (mountain),equation (14)
95      real zust                         ! U(t+dt)
96      real zvst                         ! V(t+dt)
97      real zdis                         ! =1/2*[U(t)*|U(t)|+V(t)*|V(t)|-U(t+dt)*|U(t+dt)|-V(t+dt)*|V(t+dt)|]
98      real ztemp                        ! T=-g*dtao/U_fistlevel*dp temperature due to GW stress?
99      REAL ZTAU(ndomainsz,nlayer+1)   ! Output of subrountine GWPROFILE: Gravtiy wave stress
100      REAL BV(ndomainsz,nlayer+1)     ! Brunt–Väisälä frequency at 1/2 level (output from OROSETUP input for GWSTRESS and GWPROFIL)
101      REAL ZVPH(ndomainsz,nlayer+1)   ! Low level wind speed U_H in Ref.2
102      REAL ZRHO(ndomainsz,nlayer+1)   ! density calculated (output) by the OROSETUP rountines(but not used by ORODRAG)
103                                      ! used as input for GWSTRESS and GWPROFIL
104      REAL PRI(ndomainsz,nlayer+1)    ! Mean flow richardson number at 1/2 level
105      REAL ZpsI(ndomainsz,nlayer+1)   ! The angle between the incident flow direction and the normal ridge direction of the mountain
106      REAL Zzdep(ndomainsz,nlayer)    ! dp by full level
107      REAL ZDUDT(ndomainsz)           ! du/dt due to oro-gw
108      REAL ZDVDT(ndomainsz)           ! dv/dt due to oro-gw
109      REAL ZDTDT(ndomainsz)           ! dT/dt due to oro-gw
110      REAL ZDEDT(ndomainsz)           ! zdis/ptimestemp, zdedt/Cp=dT/dt
111      REAL ZVIDIS(ndomainsz)          ! in an invalid line, not used
112      REAL Znu(ndomainsz)             ! A critical value see equation 9(Since it only compute inside OROSETUP,Maybe needs delete in future)
113      REAL Zd1(ndomainsz)             ! Bcos^2(psi)-Csin^2(psi) see equation 17 or 18 in Ref.2
114      REAL Zd2(ndomainsz)             ! (B-C)sin(psi)cos(psi)   see equation 17 or 18 in Ref.2
115      REAL Zdmod(ndomainsz)           ! sqrt(zd1^2+zd2^2)
[38]116
[2642]117!------------------------------------------------------------------------------------
118! 1.INITIALIZATION (NOT USED )
119!------------------------------------------------------------------------------------
120! 100  CONTINUE   ! continue tag without source, maybe need delete in future
121!*         1.1   Computional constants
122! 110  CONTINUE ! continue tag without source, maybe need delete in future
123!      kfdia=ndomainsz     
124!     ZTMST=TWODT
125!     IF(NSTEP.EQ.NSTART) ZTMST=0.5*TWODT
126!      nlayerM1=nlayer-1
127!      ZTMST=ptimestep
128!      ZRTMST=1./ptimestep
129! 120  CONTINUE ! continue tag without source, maybe need delete in future
130!*         1.3   Check whether row contains point for printing
131! 130  CONTINUE ! continue tag without source, maybe need delete in future
132 
133!------------------------------------------------------------------------------------
134! 2. Precompute basic state variables .
135!------------------------------------------------------------------------------------     
136! 200  CONTINUE ! continue tag without source, maybe need delete in future
137      ! Define low level wind,project winds in plane of low level wind, determine sector
138      ! in which to take the variance and set indicator for critical levels       
139      CALL OROSETUP( ngrid, nlayer, ktest, pplev, pplay, pu, pv, pt, zgeom,       &
140                           pvar,pthe, pgam,                                       &
141                           !output in capital
142                           IKCRIT, IKCRITH, ICRIT, IKENVH,IKNU,IKNU2,             &
143                           !ISECT, IKHLIM, not used
144                           ZRHO,PRI,BV,ZTAU,ZVPH,ZPSI,ZZDEP,ZNU,ZD1,ZD2,ZDMOD,    &
145                           PULOW, PVLOW)
[38]146
[2642]147!------------------------------------------------------------------------------------
148! 3. Computes low level stresses using subcritical and super critical forms.
149!    Computes anisotropy coefficient as measure of orographic two-dimensionality
150!------------------------------------------------------------------------------------             
151!  300 CONTINUE ! continue tag without source, maybe need delete in future
152      ! Computes low level stresses
153      CALL GWSTRESS(ngrid,nlayer,ktest,zrho, BV, pvar,psig,zgeom,zdmod,&
154                ! Notice that this 3 variables are actually not used due to illness lines
155                 ICRIT,IKNU,ZVPH,                                      &
156                ! not used variables
157 !                IKCRIT,ISECT,IKHLIM, IKCRITH,IKENVH,PVAR1,pgam,zd1,zd2,znu, &
158                ! not defined not used variables
159 !                 ZTFR
160                !in(as 0.0)-output:
161                ZTAU )
162               
163!------------------------------------------------------------------------------------
164! 4. Compute stress profile.
165!------------------------------------------------------------------------------------
166!  400 CONTINUE ! continue tag without source, maybe need delete in future
167      ! Compute stress profile.
168      CALL GWPROFIL( ngrid, nlayer, kgwd ,kdx, ktest, &
169                 IKCRITH, ICRIT, IKENVH, IKNU,IKNU2,pplev,  &
170                 ZRHO, BV, ZVPH, PRI,zdmod, psig , pvar,    &
171                 !not used variables
172                 !IKCRIT,ZTAUf,ZTFR,znu,pgam,
173                 ! in-output
174                 ZTAU      )
[1912]175
[2642]176!------------------------------------------------------------------------------------
177! 5. Compute tendencies.
178!------------------------------------------------------------------------------------
179!  500 CONTINUE ! continue tag without source, maybe need delete in future
180!  EXplicit solution at all levels for gravity wave
181!  Implicit solution for the blocked levels
182!      DO JL=KIDIA,KFDIA
183      ! Initialization the tendencies
184      DO JL=1,ndomainsz
185            ZVIDIS(JL)=0.0  !An invalid lines contain this variable
186            ZDUDT(JL)=0.0
187            ZDVDT(JL)=0.0
188            ZDTDT(JL)=0.0
189      ENDDO
190     
191      ! all the calculation of equation (16-18) are in the flowing loop
192      ! 1) if the wind flow over the mountain [Either the mountain is too low or the
193      !    flow higher than the blocked level], then calculate the wind tendencies by oro-GW stress directly
194      ! 2) if the flow level lower than the blocked level, then calculate the wind tendencies by obstacle drag
195      !   
196      ILEVP1=nlayer+1
197      DO  JK=1,nlayer
198        DO  JL=1,kgwd
199             ! Pick up the position of each calling points
200             JI=kdx(JL)
201             ! dp
202             ZDELP=pplev(Ji,JK+1)-pplev(Ji,JK)
203             ! T=-g*dtao/U_fistlevel*dp temperature due to GW stress?
204             ZTEMP=-g*(ZTAU(Ji,JK+1)-ZTAU(Ji,JK))/(ZVPH(Ji,ILEVP1)*ZDELP)
205             ! du/dt=(U_low*tao1-V_low*tao2)*T/(tao1^2+tao2^2)*0.5
206             ZDUDT(JI)=(PULOW(JI)*Zd1(ji)-pvlow(ji)*zd2(ji))*ztemp/zdmod(ji)
207             ! dv/dt=(V_low*tao1+U_low*tao2)*T/(tao1^2+tao2^2)*0.5
208             ZDVDT(JI)=(PVLOW(JI)*Zd1(ji)+pulow(ji)*zd2(ji))*ztemp/zdmod(ji)
209             if(jk.ge.ikenvh(ji)) then  ! IF the level lower than(notice here the level is inversed thus it is lower)
210                                        ! the top of the blocked layer (commonlly is the highest top of a mountain.)
211                                        ! Then use equation (16) to calculate the mountain drag on the incident flow
212                ! Coefficient B and C Ref.2
213                zb=1.0-0.18*pgam(ji)-0.04*pgam(ji)**2
214                zc=0.48*pgam(ji)+0.3*pgam(ji)**2
215                ! Middle part of the equation (16)
216                zconb=2.*ptimestep*GKWAKE*psig(ji)/(4.*pvar(ji))
217                ! Tail part of the equation (16),U*|U|/2+V*|V|/2
218                zabsv=sqrt(pu(JI,JK)**2+pv(JI,JK)**2)/2.
219                ! Second tail part of the equation (16)
220                zzd1=zb*cos(zpsi(ji,jk))**2+zc*sin(zpsi(ji,jk))**2
221                ! Aspect ratio of the obstacle(topography), equation(14)
222                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)
223                ! Equation (16)
224                zbet=max(0.,2.-1./ratio)*zconb*zzdep(ji,jk)*zzd1*zabsv
225                ! du/dt and dv/dt due to the mountain drag
226                zdudt(ji)=-pu(ji,jk)/ptimestep
227                zdvdt(ji)=-pv(ji,jk)/ptimestep
228                zdudt(ji)=zdudt(ji)*(zbet/(1.+zbet))
229                zdvdt(ji)=zdvdt(ji)*(zbet/(1.+zbet))
230             end if
231             ! wind tendencies due to the orographic gravity wave (due to stress(flow over) and drag (blocked))
232             pdudt(JI,JK)=ZDUDT(JI)
233             pdvdt(JI,JK)=ZDVDT(JI)
234             ! winds at t+dt (oro-gw induced increaments are added)
235             ZUST=pu(JI,JK)+ptimestep*ZDUDT(JI)
236             ZVST=pv(JI,JK)+ptimestep*ZDVDT(JI)
237             ! =1/2*[U(t)*|U(t)|+V(t)*|V(t)|-U(t+dt)*|U(t+dt)|-V(t+dt)*|V(t+dt)|]
238             ZDIS=0.5*(pu(JI,JK)**2+pv(JI,JK)**2-ZUST**2-ZVST**2)
239             ZDEDT(JI)=ZDIS/ptimestep
240             ! This line no use maybe delete in the future
241             ZVIDIS(JI)=ZVIDIS(JI)+ZDIS*ZDELP
242             ! Temperature tendencies due to the orographic gravity wave (due to stress(flow over) and drag (blocked))
243             ZDTDT(JI)=ZDEDT(JI)/cpp
244             pdtdt(JI,JK)=ZDTDT(JI)
245        ENDDO !JL=1,kgwd
246      end DO !JK=1,nlayer
247
[1912]248      END SUBROUTINE ORODRAG
249     
[2642]250END MODULE orodrag_mod
Note: See TracBrowser for help on using the repository browser.