source: trunk/LMDZ.VENUS/libf/phyvenus/orodrag.F @ 3094

Last change on this file since 3094 was 2047, checked in by slebonnois, 6 years ago

SL: VENUS, ajout des modifs apportees par Thomas Navarro pour la parametrisation des ondes de gravite orographiques

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