source: trunk/LMDZ.MARS/libf/phymars/orodrag_mod.F @ 1934

Last change on this file since 1934 was 1913, checked in by emillour, 7 years ago

Mars GCM:

  • Forgotten in previous commit: gwprofil.F -> gwprofil_mod.F (here also the

size of an argument, rho, was incorrect in caller orodrag).

  • Turned newsedim.F into a module newsedim_mod.F
  • Adapted co2cloud.F and improvedCO2clouds.F to not use "newunit" to open file

(it is perfectly legitimate F2008 Fortran, but older compiler such as gfortran
on local LMD machines are not there yet).
EM
File size: 8.4 KB
Line 
1      MODULE orodrag_mod
2     
3      IMPLICIT NONE
4     
5      CONTAINS
6     
7      SUBROUTINE ORODRAG( klon,klev
8     I                 , KGWD, KGWDIM, KDX, KTEST
9     R                 , PTSPHY
10     R                 , PAPHM1,PAPM1,PGEOM1,PTM1,PUM1
11     R                 , PVM1, PVAROR, PSIG, PGAMMA, PTHETA
12C OUTPUTS
13     R                 , PULOW,PVLOW
14     R                 , PVOM,PVOL,PTE )
15C
16C
17C**** *ORODRAG* - DOES THE GRAVITY WAVE PARAMETRIZATION.
18C
19C     PURPOSE.
20C     --------
21C
22C          THIS ROUTINE COMPUTES THE PHYSICAL TENDENCIES OF THE
23C          PROGNOSTIC VARIABLES U,V  AND T DUE TO  VERTICAL TRANSPORTS BY
24C          SUBGRIDSCALE OROGRAPHICALLY EXCITED GRAVITY WAVES
25C
26C     EXPLICIT ARGUMENTS :
27C     --------------------
28C
29C     INPUT :
30C
31C     NLON               : NUMBER OF HORIZONTAL GRID POINTS
32C     NLEV               : NUMBER OF LEVELS
33C     KGWD               : NUMBER OF POINTS AT WHICH THE SCHEME IS CALLED
34C     KGWDIM             : NUMBER OF POINTS AT WHICH THE SCHEME IS CALLED
35C     KDX(NLON)          : POINTS AT WHICH TO CALL THE SCHEME
36C     KTEST(NLON)        : MAP OF CALLING POINTS
37C     PTSPHY             : LENGTH OF TIME STEP
38C     PAPHM1(NLON,NLEV+1): PRESSURE AT MIDDLE LEVELS
39C     PAPM1(NLON,NLEV)   : PRESSURE ON MODEL LEVELS
40C     PGEOM1(NLON,NLEV)  : GEOPOTENTIAL HIEGHT OF MODEL LEVELS
41C     PTM1(NLON,NLEV)    : TEMPERATURE
42C     PUM1(NLON,NLEV)    : ZONAL WIND
43C     PVM1(NLON,NLEV)    : MERIDIONAL WIND
44C     PVAROR(NLON)       : SUB-GRID SCALE STANDARD DEVIATION
45C     PSIG(NLON)         : SUB-GRID SCALE SLOPE
46C     PGAMMA(NLON)       : SUB-GRID SCALE ANISOTROPY
47C     PTHETA(NLON)       : SUB-GRID SCALE PRINCIPAL AXES ANGLE
48C
49C     OUTPUT :
50C
51C     PULOW(NLON)        : LOW LEVEL ZONAL WIND
52C     PVLOW(NLON)        : LOW LEVEL MERIDIONAL WIND
53C     PVOM(NLON,NLEV)    : ZONAL WIND TENDENCY
54C     PVOL(NLON,NLEV)    : MERIDIONAL WIND TENDENCY
55C     PTE(NLON,NLEV)     : TEMPERATURE TENDENCY
56C
57C     IMPLICIT ARGUMENTS :
58C     --------------------
59C
60C     comcstfi.h
61C     yoegwd.h
62C
63C     METHOD.
64C     -------
65C
66C     EXTERNALS.
67C     ----------
68C
69C     REFERENCE.
70C     ----------
71C
72C     AUTHOR.
73C     -------
74C     M.MILLER + B.RITTER   E.C.M.W.F.     15/06/86.
75C
76C     F.LOTT + M. MILLER    E.C.M.W.F.     22/11/94
77C-----------------------------------------------------------------------
78      use dimradmars_mod, only: ndlo2
79      USE gwstress_mod, ONLY: gwstress
80      USE gwprofil_mod, ONLY: gwprofil
81      USE comcstfi_h, ONLY: g, cpp
82      implicit none
83C
84C
85      integer klon,klev,kidia
86      parameter(kidia=1)
87      integer, save :: kfdia ! =NDLO2
88
89      include "yoegwd.h"
90C-----------------------------------------------------------------------
91C
92C*       0.1   ARGUMENTS
93C              ---------
94C
95C
96      REAL  PTE(NDLO2,klev),
97     *      PVOL(NDLO2,klev),
98     *      PVOM(NDLO2,klev),
99     *      PULOW(NDLO2),
100     *      PVLOW(NDLO2)
101      REAL  PUM1(NDLO2,klev),
102     *      PVM1(NDLO2,klev),
103     *      PTM1(NDLO2,klev),
104     *      PVAROR(NDLO2),PSIG(NDLO2),PGAMMA(NDLO2),PTHETA(NDLO2),
105     *      PGEOM1(NDLO2,klev),
106     *      PAPM1(NDLO2,klev),
107     *      PAPHM1(NDLO2,klev+1)
108C
109      integer kgwd,kgwdim
110      real ptsphy
111      INTEGER  KDX(NDLO2),KTEST(NDLO2)
112C-----------------------------------------------------------------------
113C
114C*       0.2   LOCAL ARRAYS
115C              ------------
116      INTEGER  ISECT(NDLO2),
117     *         ICRIT(NDLO2),
118     *         IKCRITH(NDLO2),
119     *         IKenvh(NDLO2),
120     *         IKNU(NDLO2),
121     *         IKNU2(NDLO2),
122     *         IKCRIT(NDLO2),
123     *         IKHLIM(NDLO2)
124      integer ji,jk,jl,klevm1,ilevp1
125C      real gkwake
126      real ztmst,pvar(NDLO2,4),ztauf(NDLO2,klev+1)
127      real zrtmst,zdelp,zb,zc,zbet
128      real zconb,zabsv,zzd1,ratio,zust,zvst,zdis,ztemp
129C
130      REAL   ZTAU(NDLO2,klev+1),
131     *       ZSTAB(NDLO2,klev+1),
132     *       ZVPH(NDLO2,klev+1),
133     *       ZRHO(NDLO2,klev+1),
134     *       ZRI(NDLO2,klev+1),
135     *       ZpsI(NDLO2,klev+1),
136     *       Zzdep(NDLO2,klev)
137      REAL   ZDUDT(NDLO2),
138     *       ZDVDT(NDLO2),
139     *       ZDTDT(NDLO2),
140     *       ZDEDT(NDLO2),
141     *       ZVIDIS(NDLO2),
142     *       ZTFR(NDLO2),
143     *       Znu(NDLO2),
144     *       Zd1(NDLO2),
145     *       Zd2(NDLO2),
146     *       Zdmod(NDLO2)
147C
148C------------------------------------------------------------------
149C
150C*         1.    INITIALIZATION
151C                --------------
152C
153 100  CONTINUE
154C
155C     ------------------------------------------------------------------
156C
157C*         1.1   COMPUTATIONAL CONSTANTS
158C                -----------------------
159C
160 110  CONTINUE
161C
162      kfdia=NDLO2
163     
164c     ZTMST=TWODT
165c     IF(NSTEP.EQ.NSTART) ZTMST=0.5*TWODT
166      KLEVM1=KLEV-1
167      ZTMST=PTSPHY
168      ZRTMST=1./ZTMST
169C     ------------------------------------------------------------------
170C
171 120  CONTINUE
172C
173C     ------------------------------------------------------------------
174C
175C*         1.3   CHECK WHETHER ROW CONTAINS POINT FOR PRINTING
176C                ---------------------------------------------
177C
178 130  CONTINUE
179C
180C     ------------------------------------------------------------------
181C
182C*         2.     PRECOMPUTE BASIC STATE VARIABLES.
183C*                ---------- ----- ----- ----------
184C*                DEFINE LOW LEVEL WIND, PROJECT WINDS IN PLANE OF
185C*                LOW LEVEL WIND, DETERMINE SECTOR IN WHICH TO TAKE
186C*                THE VARIANCE AND SET INDICATOR FOR CRITICAL LEVELS.
187C
188  200 CONTINUE
189C
190C
191C
192      CALL OROSETUP
193     *     ( klon, klev , KTEST
194     *     , IKCRIT, IKCRITH, ICRIT, ISECT, IKHLIM, ikenvh,IKNU,iknu2
195     *     , PAPHM1, PAPM1 , PUM1   , PVM1 , PTM1 , PGEOM1, pvaror
196     *     , ZRHO  , ZRI   , ZSTAB  , ZTAU , ZVPH , zpsi, zzdep
197     *     , PULOW, PVLOW
198     *     , ptheta,pgamma,znu  ,zd1,  zd2,  zdmod                )
199C
200C
201C
202C***********************************************************
203C
204C
205C*         3.      COMPUTE LOW LEVEL STRESSES USING SUBCRITICAL AND
206C*                 SUPERCRITICAL FORMS.COMPUTES ANISOTROPY COEFFICIENT
207C*                 AS MEASURE OF OROGRAPHIC TWODIMENSIONALITY.
208C
209  300 CONTINUE
210C
211      CALL GWSTRESS
212     *    ( klon  , klev
213     *    , IKCRIT, ISECT, IKHLIM, KTEST, IKCRITH, ICRIT, ikenvh, IKNU
214     *    , ZRHO  , ZSTAB, ZVPH  , PVAR , pvaror,  psig
215     *    , ZTFR   , ZTAU
216     *    , pgeom1,pgamma,zd1,zd2,zdmod,znu)
217C
218C*         4.      COMPUTE STRESS PROFILE.
219C*                 ------- ------ --------
220C
221  400 CONTINUE
222C
223C
224      CALL GWPROFIL
225     *       (  klon , klev
226     *       , kgwd   , kdx  , KTEST
227     *       , IKCRIT, IKCRITH, ICRIT  , ikenvh, IKNU
228     *       ,iknu2 , pAPHM1, ZRHO   , ZSTAB , ZTFR   , ZVPH
229     *       , ZRI   , ZTAU   , ztauf
230     *       , zdmod , znu    , psig  , pgamma , pvaror   )
231C
232C
233C*         5.      COMPUTE TENDENCIES.
234C*                 -------------------
235C
236  500 CONTINUE
237C
238C  EXPLICIT SOLUTION AT ALL LEVELS FOR THE GRAVITY WAVE
239C  IMPLICIT SOLUTION FOR THE BLOCKED LEVELS
240
241      DO 510 JL=KIDIA,KFDIA
242      ZVIDIS(JL)=0.0
243      ZDUDT(JL)=0.0
244      ZDVDT(JL)=0.0
245      ZDTDT(JL)=0.0
246  510 CONTINUE
247C
248      ILEVP1=KLEV+1
249C
250C
251      DO 524 JK=1,klev
252C
253CDIR$ IVDEP
254C
255C      GKWAKE=0.5
256C
257C     NOW SET IN SUGWD.F
258C
259      DO 523 JL=1,KGWD
260      JI=KDX(JL)
261      ZDELP=pAPHM1(Ji,JK+1)-pAPHM1(Ji,JK)
262      ZTEMP=-g*(ZTAU(Ji,JK+1)-ZTAU(Ji,JK))/(ZVPH(Ji,ILEVP1)*ZDELP)
263      ZDUDT(JI)=(PULOW(JI)*Zd1(ji)-pvlow(ji)*zd2(ji))*ztemp/zdmod(ji)
264      ZDVDT(JI)=(pvLOW(JI)*Zd1(ji)+pulow(ji)*zd2(ji))*ztemp/zdmod(ji)
265      if(jk.ge.ikenvh(ji)) then
266         zb=1.0-0.18*pgamma(ji)-0.04*pgamma(ji)**2
267         zc=0.48*pgamma(ji)+0.3*pgamma(ji)**2
268         zconb=2.*ztmst*GKWAKE*psig(ji)/(4.*pvaror(ji))
269         zabsv=sqrt(PUM1(JI,JK)**2+PVM1(JI,JK)**2)/2.
270         zzd1=zb*cos(zpsi(ji,jk))**2+zc*sin(zpsi(ji,jk))**2
271         ratio=(cos(zpsi(ji,jk))**2+pgamma(ji)*sin(zpsi(ji,jk))**2)/
272     *   (pgamma(ji)*cos(zpsi(ji,jk))**2+sin(zpsi(ji,jk))**2)
273         zbet=max(0.,2.-1./ratio)*zconb*zzdep(ji,jk)*zzd1*zabsv
274         zdudt(ji)=-pum1(ji,jk)/ztmst
275         zdvdt(ji)=-pvm1(ji,jk)/ztmst
276         zdudt(ji)=zdudt(ji)*(zbet/(1.+zbet))
277         zdvdt(ji)=zdvdt(ji)*(zbet/(1.+zbet))
278      end if
279      PVOM(JI,JK)=ZDUDT(JI)
280      PVOL(JI,JK)=ZDVDT(JI)
281      ZUST=PUM1(JI,JK)+ZTMST*ZDUDT(JI)
282      ZVST=PVM1(JI,JK)+ZTMST*ZDVDT(JI)
283      ZDIS=0.5*(PUM1(JI,JK)**2+PVM1(JI,JK)**2-ZUST**2-ZVST**2)
284      ZDEDT(JI)=ZDIS/ZTMST
285      ZVIDIS(JI)=ZVIDIS(JI)+ZDIS*ZDELP
286      ZDTDT(JI)=ZDEDT(JI)/cpp
287      PTE(JI,JK)=ZDTDT(JI)
288
289  523 CONTINUE
290
291  524 CONTINUE
292C
293C
294
295      END SUBROUTINE ORODRAG
296     
297      END MODULE orodrag_mod
Note: See TracBrowser for help on using the repository browser.