source: trunk/LMDZ.TITAN.old/libf/phytitan/orodrag.F @ 3094

Last change on this file since 3094 was 1530, checked in by emillour, 9 years ago

Venus and Titan GCMs:
Updates in the physics to keep up with updates in LMDZ5 (up to
LMDZ5 trunk, rev 2350) concerning dynamics/physics separation:

  • Adapted makelmdz and makelmdz_fcm script to stop if trying to compile 1d model or newstart or start2archive in parallel.
  • got rid of references to "dimensions.h" in physics. Within physics packages, use nbp_lon (=iim), nbp_lat (=jjmp1) and nbp_lev (=llm) from module mod_grid_phy_lmdz (in phy_common) instead. Only partially done for Titan, because of many hard-coded commons; a necessary first step will be to clean these up (using modules).

EM
File size: 9.7 KB
Line 
1      SUBROUTINE orodrag( nlon,nlev
2     i                 , kgwd,  kdx, ktest
3     r                 , ptsphy
4     r                 , paphm1,papm1,pgeom1,pn2m1,ptm1,pum1,pvm1
5     r                 , pmea, pstd, psig, pgam, pthe, ppic, pval
6c outputs
7     r                 , pulow,pvlow
8     r                 , pvom,pvol,pte )
9     
10      use dimphy
11      IMPLICIT NONE
12
13c
14c
15c**** *orodrag* - does the SSO drag  parametrization.
16c
17c     purpose.
18c     --------
19c
20c     this routine computes the physical tendencies of the
21c     prognostic variables u,v  and t due to  vertical transports by
22c     subgridscale orographically excited gravity waves, and to
23c     low level blocked flow drag.
24c
25c**   interface.
26c     ----------
27c          called from *drag_noro*.
28c
29c          the routine takes its input from the long-term storage:
30c          u,v,t and p at t-1.
31c
32c        explicit arguments :
33c        --------------------
34c     ==== inputs ===
35c nlon----input-I-Total number of horizontal points that get into physics
36c nlev----input-I-Number of vertical levels
37c
38c kgwd- -input-I: Total nb of points where the orography schemes are active
39c ktest--input-I: Flags to indicate active points
40c kdx----input-I: Locate the physical location of an active point.
41c ptsphy--input-R-Time-step (s)
42c paphm1--input-R: pressure at model 1/2 layer
43c papm1---input-R: pressure at model layer
44c pgeom1--input-R: Altitude of layer above ground
45c pn2m1---input-R-Brunt-Vaisala freq.^2 at 1/2 layers
46c ptm1, pum1, pvm1--R-: t, u and v
47c pmea----input-R-Mean Orography (m)
48C pstd----input-R-SSO standard deviation (m)
49c psig----input-R-SSO slope
50c pgam----input-R-SSO Anisotropy
51c pthe----input-R-SSO Angle
52c ppic----input-R-SSO Peacks elevation (m)
53c pval----input-R-SSO Valleys elevation (m)
54
55      integer nlon,nlev,kgwd
56      real ptsphy
57
58c     ==== outputs ===
59c pulow, pvlow -output-R: Low-level wind
60c
61c pte -----output-R: T tendency
62c pvom-----output-R: U tendency
63c pvol-----output-R: V tendency
64c
65c
66c Implicit Arguments:
67c ===================
68c
69c klon-common-I: Number of points seen by the physics
70c klev-common-I: Number of vertical layers
71c
72c     method.
73c     -------
74c
75c     externals.
76c     ----------
77Coff  integer ismin, ismax
78Coff  external ismin, ismax
79c
80c     reference.
81c     ----------
82c
83c     author.
84c     -------
85c     m.miller + b.ritter   e.c.m.w.f.     15/06/86.
86c
87c     f.lott + m. miller    e.c.m.w.f.     22/11/94
88c-----------------------------------------------------------------------
89c
90c
91#include "YOMCST.h"
92#include "YOEGWD.h"
93
94c-----------------------------------------------------------------------
95c
96c*       0.1   arguments
97c              ---------
98c
99c
100      real  pte(nlon,nlev),
101     *      pvol(nlon,nlev),
102     *      pvom(nlon,nlev),
103     *      pulow(nlon),
104     *      pvlow(nlon)
105      real  pum1(nlon,nlev),
106     *      pvm1(nlon,nlev),
107     *      ptm1(nlon,nlev),
108     *      pmea(nlon),pstd(nlon),psig(nlon),
109     *      pgam(nlon),pthe(nlon),ppic(nlon),pval(nlon),
110     *      pgeom1(nlon,nlev),pn2m1(nlon,nlev),
111     *      papm1(nlon,nlev),
112     *      paphm1(nlon,nlev+1)
113c
114      integer  kdx(nlon),ktest(nlon)
115c-----------------------------------------------------------------------
116c
117c*       0.2   local arrays
118c              ------------
119      integer  isect(klon),
120     *         icrit(klon),
121     *         ikcrith(klon),
122     *         ikenvh(klon),
123     *         iknu(klon),
124     *         iknu2(klon),
125     *         ikcrit(klon),
126     *         ikhlim(klon)
127c
128      real   ztau(klon,klev+1),
129     *       zstab(klon,klev+1),
130     *       zvph(klon,klev+1),
131     *       zrho(klon,klev+1),
132     *       zri(klon,klev+1),
133     *       zpsi(klon,klev+1),
134     *       zzdep(klon,klev)
135      real   zdudt(klon),
136     *       zdvdt(klon),
137     *       zdtdt(klon),
138     *       zdedt(klon),
139     *       zvidis(klon),
140     *       ztfr(klon),
141     *       znu(klon),
142     *       zd1(klon),
143     *       zd2(klon),
144     *       zdmod(klon)
145
146
147c local quantities:
148
149      integer jl,jk,ji
150      real ztmst,zdelp,ztemp,zforc,ztend,rover               
151      real zb,zc,zconb,zabsv,zzd1,ratio,zbet,zust,zvst,zdis
152   
153c
154c------------------------------------------------------------------
155c
156c*         1.    initialization
157c                --------------
158c
159c        print *,' in orodrag'
160 100  continue
161c
162c     ------------------------------------------------------------------
163c
164c*         1.1   computational constants
165c                -----------------------
166c
167 110  continue
168c
169c     ztmst=twodt
170c     if(nstep.eq.nstart) ztmst=0.5*twodt
171      ztmst=ptsphy
172c     ------------------------------------------------------------------
173c
174 120  continue
175c
176c     ------------------------------------------------------------------
177c
178c*         1.3   check whether row contains point for printing
179c                ---------------------------------------------
180c
181 130  continue
182c
183c     ------------------------------------------------------------------
184c
185c*         2.     precompute basic state variables.
186c*                ---------- ----- ----- ----------
187c*                define low level wind, project winds in plane of
188c*                low level wind, determine sector in which to take
189c*                the variance and set indicator for critical levels.
190c
191
192  200 continue
193c
194      do jk=1,klev
195       zstab(:,jk) = pn2m1(:,jk)
196      enddo
197c
198      call orosetup
199     *     ( nlon, nlev , ktest
200     *     , ikcrit, ikcrith, icrit, isect, ikhlim, ikenvh,iknu,iknu2
201     *     , paphm1, papm1 , pum1   , pvm1 , ptm1 , pgeom1, zstab, pstd
202     *     , zrho  , zri   , ztau , zvph , zpsi, zzdep
203     *     , pulow, pvlow
204     *     , pthe,pgam,pmea,ppic,pval,znu  ,zd1,  zd2,  zdmod )
205
206c
207c
208c
209c***********************************************************
210c
211c
212c*         3.      compute low level stresses using subcritical and
213c*                 supercritical forms.computes anisotropy coefficient
214c*                 as measure of orographic twodimensionality.
215c
216  300 continue
217c
218      call gwstress
219     *    ( nlon  , nlev
220     *    , ikcrit, isect, ikhlim, ktest, ikcrith, icrit, ikenvh, iknu
221     *    , zrho  , zstab, zvph  , pstd,  psig, pmea, ppic, pval
222     *    , ztfr   , ztau
223     *    , pgeom1,pgam,zd1,zd2,zdmod,znu)
224
225c
226c
227c*         4.      compute stress profile including
228c                  trapped waves, wave breaking,
229c                  linear decay in stratosphere.
230c
231  400 continue
232c
233c
234
235      call gwprofil
236     *       (  nlon , nlev
237     *       , kgwd   , kdx  , ktest
238     *       , ikcrit, ikcrith, icrit  , ikenvh, iknu
239     *       ,iknu2 , paphm1, zrho   , zstab , ztfr   , zvph
240     *       , zri   , ztau
241 
242     *       , zdmod , znu    , psig  , pgam , pstd , ppic , pval)
243
244c
245c*         5.      Compute tendencies from waves stress profile.
246c                  Compute low level blocked flow drag.
247c*                 --------------------------------------------
248c
249  500 continue
250
251     
252c
253c  explicit solution at all levels for the gravity wave
254c  implicit solution for the blocked levels
255
256      do 510 jl=kidia,kfdia
257      zvidis(jl)=0.0
258      zdudt(jl)=0.0
259      zdvdt(jl)=0.0
260      zdtdt(jl)=0.0
261  510 continue
262c
263
264      do 524 jk=1,klev
265c
266
267C  WAVE STRESS
268C-------------
269c
270c
271      do 523 ji=kidia,kfdia
272
273      if(ktest(ji).eq.1) then
274
275      zdelp=paphm1(ji,jk+1)-paphm1(ji,jk)
276      ztemp=-rg*(ztau(ji,jk+1)-ztau(ji,jk))/(zvph(ji,klev+1)*zdelp)
277
278      zdudt(ji)=(pulow(ji)*zd1(ji)-pvlow(ji)*zd2(ji))*ztemp/zdmod(ji)
279      zdvdt(ji)=(pvlow(ji)*zd1(ji)+pulow(ji)*zd2(ji))*ztemp/zdmod(ji)
280c
281c Control Overshoots
282c
283
284      if(jk.ge.ntop)then
285        rover=0.10
286        if(abs(zdudt(ji)).gt.rover*abs(pum1(ji,jk))/ztmst)
287     C    zdudt(ji)=rover*abs(pum1(ji,jk))/ztmst*
288     C              zdudt(ji)/(abs(zdudt(ji))+1.E-10)
289        if(abs(zdvdt(ji)).gt.rover*abs(pvm1(ji,jk))/ztmst)
290     C    zdvdt(ji)=rover*abs(pvm1(ji,jk))/ztmst*
291     C              zdvdt(ji)/(abs(zdvdt(ji))+1.E-10)
292      endif
293
294      rover=0.25
295      zforc=sqrt(zdudt(ji)**2+zdvdt(ji)**2)       
296      ztend=sqrt(pum1(ji,jk)**2+pvm1(ji,jk)**2)/ztmst                     
297
298      if(zforc.ge.rover*ztend)then
299        zdudt(ji)=rover*ztend/zforc*zdudt(ji)
300        zdvdt(ji)=rover*ztend/zforc*zdvdt(ji)
301      endif
302c
303c BLOCKED FLOW DRAG:
304C -----------------
305c
306      if(jk.gt.ikenvh(ji)) then
307         zb=1.0-0.18*pgam(ji)-0.04*pgam(ji)**2
308         zc=0.48*pgam(ji)+0.3*pgam(ji)**2
309         zconb=2.*ztmst*gkwake*psig(ji)/(4.*pstd(ji))
310         zabsv=sqrt(pum1(ji,jk)**2+pvm1(ji,jk)**2)/2.
311         zzd1=zb*cos(zpsi(ji,jk))**2+zc*sin(zpsi(ji,jk))**2
312         ratio=(cos(zpsi(ji,jk))**2+pgam(ji)*sin(zpsi(ji,jk))**2)/
313     *   (pgam(ji)*cos(zpsi(ji,jk))**2+sin(zpsi(ji,jk))**2)
314         zbet=max(0.,2.-1./ratio)*zconb*zzdep(ji,jk)*zzd1*zabsv
315c
316c OPPOSED TO THE WIND
317c
318         zdudt(ji)=-pum1(ji,jk)/ztmst
319         zdvdt(ji)=-pvm1(ji,jk)/ztmst
320c
321c PERPENDICULAR TO THE SSO MAIN AXIS:
322C                           
323cmod     zdudt(ji)=-(pum1(ji,jk)*cos(pthe(ji)*rpi/180.)
324cmod *              +pvm1(ji,jk)*sin(pthe(ji)*rpi/180.))
325cmod *              *cos(pthe(ji)*rpi/180.)/ztmst
326cmod     zdvdt(ji)=-(pum1(ji,jk)*cos(pthe(ji)*rpi/180.)
327cmod *              +pvm1(ji,jk)*sin(pthe(ji)*rpi/180.))
328cmod *              *sin(pthe(ji)*rpi/180.)/ztmst
329C
330         zdudt(ji)=zdudt(ji)*(zbet/(1.+zbet))
331         zdvdt(ji)=zdvdt(ji)*(zbet/(1.+zbet))
332      end if
333      pvom(ji,jk)=zdudt(ji)
334      pvol(ji,jk)=zdvdt(ji)
335      zust=pum1(ji,jk)+ztmst*zdudt(ji)
336      zvst=pvm1(ji,jk)+ztmst*zdvdt(ji)
337      zdis=0.5*(pum1(ji,jk)**2+pvm1(ji,jk)**2-zust**2-zvst**2)
338      zdedt(ji)=zdis/ztmst
339      zvidis(ji)=zvidis(ji)+zdis*zdelp
340c VENUS ATTENTION: CP VARIABLE
341      zdtdt(ji)=zdedt(ji)/rcpd
342c
343c  NO TENDENCIES ON TEMPERATURE .....
344c
345c  Instead of, pte(ji,jk)=zdtdt(ji), due to mechanical dissipation
346c
347      pte(ji,jk)=0.0
348
349      endif
350
351  523 continue
352  524 continue
353c
354c
355  501 continue
356
357      return
358      end
359
Note: See TracBrowser for help on using the repository browser.