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

Last change on this file since 3807 was 3739, checked in by emillour, 7 weeks ago

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