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

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