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

Last change on this file since 1243 was 1056, checked in by slebonnois, 11 years ago

SL: Titan runs ! see DOC/chantiers/commit_importants.log

File size: 9.7 KB
RevLine 
[1056]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
13#include "dimensions.h"
14#include "paramet.h"
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)
116c
117      integer  kdx(nlon),ktest(nlon)
118c-----------------------------------------------------------------------
119c
120c*       0.2   local arrays
121c              ------------
122      integer  isect(klon),
123     *         icrit(klon),
124     *         ikcrith(klon),
125     *         ikenvh(klon),
126     *         iknu(klon),
127     *         iknu2(klon),
128     *         ikcrit(klon),
129     *         ikhlim(klon)
130c
131      real   ztau(klon,klev+1),
132     *       zstab(klon,klev+1),
133     *       zvph(klon,klev+1),
134     *       zrho(klon,klev+1),
135     *       zri(klon,klev+1),
136     *       zpsi(klon,klev+1),
137     *       zzdep(klon,klev)
138      real   zdudt(klon),
139     *       zdvdt(klon),
140     *       zdtdt(klon),
141     *       zdedt(klon),
142     *       zvidis(klon),
143     *       ztfr(klon),
144     *       znu(klon),
145     *       zd1(klon),
146     *       zd2(klon),
147     *       zdmod(klon)
148
149
150c local quantities:
151
152      integer jl,jk,ji
153      real ztmst,zdelp,ztemp,zforc,ztend,rover               
154      real zb,zc,zconb,zabsv,zzd1,ratio,zbet,zust,zvst,zdis
155   
156c
157c------------------------------------------------------------------
158c
159c*         1.    initialization
160c                --------------
161c
162c        print *,' in orodrag'
163 100  continue
164c
165c     ------------------------------------------------------------------
166c
167c*         1.1   computational constants
168c                -----------------------
169c
170 110  continue
171c
172c     ztmst=twodt
173c     if(nstep.eq.nstart) ztmst=0.5*twodt
174      ztmst=ptsphy
175c     ------------------------------------------------------------------
176c
177 120  continue
178c
179c     ------------------------------------------------------------------
180c
181c*         1.3   check whether row contains point for printing
182c                ---------------------------------------------
183c
184 130  continue
185c
186c     ------------------------------------------------------------------
187c
188c*         2.     precompute basic state variables.
189c*                ---------- ----- ----- ----------
190c*                define low level wind, project winds in plane of
191c*                low level wind, determine sector in which to take
192c*                the variance and set indicator for critical levels.
193c
194
195  200 continue
196c
197      do jk=1,klev
198       zstab(:,jk) = pn2m1(:,jk)
199      enddo
200c
201      call orosetup
202     *     ( nlon, nlev , ktest
203     *     , ikcrit, ikcrith, icrit, isect, ikhlim, ikenvh,iknu,iknu2
204     *     , paphm1, papm1 , pum1   , pvm1 , ptm1 , pgeom1, zstab, pstd
205     *     , zrho  , zri   , ztau , zvph , zpsi, zzdep
206     *     , pulow, pvlow
207     *     , pthe,pgam,pmea,ppic,pval,znu  ,zd1,  zd2,  zdmod )
208
209c
210c
211c
212c***********************************************************
213c
214c
215c*         3.      compute low level stresses using subcritical and
216c*                 supercritical forms.computes anisotropy coefficient
217c*                 as measure of orographic twodimensionality.
218c
219  300 continue
220c
221      call gwstress
222     *    ( nlon  , nlev
223     *    , ikcrit, isect, ikhlim, ktest, ikcrith, icrit, ikenvh, iknu
224     *    , zrho  , zstab, zvph  , pstd,  psig, pmea, ppic, pval
225     *    , ztfr   , ztau
226     *    , pgeom1,pgam,zd1,zd2,zdmod,znu)
227
228c
229c
230c*         4.      compute stress profile including
231c                  trapped waves, wave breaking,
232c                  linear decay in stratosphere.
233c
234  400 continue
235c
236c
237
238      call gwprofil
239     *       (  nlon , nlev
240     *       , kgwd   , kdx  , ktest
241     *       , ikcrit, ikcrith, icrit  , ikenvh, iknu
242     *       ,iknu2 , paphm1, zrho   , zstab , ztfr   , zvph
243     *       , zri   , ztau
244 
245     *       , zdmod , znu    , psig  , pgam , pstd , ppic , pval)
246
247c
248c*         5.      Compute tendencies from waves stress profile.
249c                  Compute low level blocked flow drag.
250c*                 --------------------------------------------
251c
252  500 continue
253
254     
255c
256c  explicit solution at all levels for the gravity wave
257c  implicit solution for the blocked levels
258
259      do 510 jl=kidia,kfdia
260      zvidis(jl)=0.0
261      zdudt(jl)=0.0
262      zdvdt(jl)=0.0
263      zdtdt(jl)=0.0
264  510 continue
265c
266
267      do 524 jk=1,klev
268c
269
270C  WAVE STRESS
271C-------------
272c
273c
274      do 523 ji=kidia,kfdia
275
276      if(ktest(ji).eq.1) then
277
278      zdelp=paphm1(ji,jk+1)-paphm1(ji,jk)
279      ztemp=-rg*(ztau(ji,jk+1)-ztau(ji,jk))/(zvph(ji,klev+1)*zdelp)
280
281      zdudt(ji)=(pulow(ji)*zd1(ji)-pvlow(ji)*zd2(ji))*ztemp/zdmod(ji)
282      zdvdt(ji)=(pvlow(ji)*zd1(ji)+pulow(ji)*zd2(ji))*ztemp/zdmod(ji)
283c
284c Control Overshoots
285c
286
287      if(jk.ge.ntop)then
288        rover=0.10
289        if(abs(zdudt(ji)).gt.rover*abs(pum1(ji,jk))/ztmst)
290     C    zdudt(ji)=rover*abs(pum1(ji,jk))/ztmst*
291     C              zdudt(ji)/(abs(zdudt(ji))+1.E-10)
292        if(abs(zdvdt(ji)).gt.rover*abs(pvm1(ji,jk))/ztmst)
293     C    zdvdt(ji)=rover*abs(pvm1(ji,jk))/ztmst*
294     C              zdvdt(ji)/(abs(zdvdt(ji))+1.E-10)
295      endif
296
297      rover=0.25
298      zforc=sqrt(zdudt(ji)**2+zdvdt(ji)**2)       
299      ztend=sqrt(pum1(ji,jk)**2+pvm1(ji,jk)**2)/ztmst                     
300
301      if(zforc.ge.rover*ztend)then
302        zdudt(ji)=rover*ztend/zforc*zdudt(ji)
303        zdvdt(ji)=rover*ztend/zforc*zdvdt(ji)
304      endif
305c
306c BLOCKED FLOW DRAG:
307C -----------------
308c
309      if(jk.gt.ikenvh(ji)) then
310         zb=1.0-0.18*pgam(ji)-0.04*pgam(ji)**2
311         zc=0.48*pgam(ji)+0.3*pgam(ji)**2
312         zconb=2.*ztmst*gkwake*psig(ji)/(4.*pstd(ji))
313         zabsv=sqrt(pum1(ji,jk)**2+pvm1(ji,jk)**2)/2.
314         zzd1=zb*cos(zpsi(ji,jk))**2+zc*sin(zpsi(ji,jk))**2
315         ratio=(cos(zpsi(ji,jk))**2+pgam(ji)*sin(zpsi(ji,jk))**2)/
316     *   (pgam(ji)*cos(zpsi(ji,jk))**2+sin(zpsi(ji,jk))**2)
317         zbet=max(0.,2.-1./ratio)*zconb*zzdep(ji,jk)*zzd1*zabsv
318c
319c OPPOSED TO THE WIND
320c
321         zdudt(ji)=-pum1(ji,jk)/ztmst
322         zdvdt(ji)=-pvm1(ji,jk)/ztmst
323c
324c PERPENDICULAR TO THE SSO MAIN AXIS:
325C                           
326cmod     zdudt(ji)=-(pum1(ji,jk)*cos(pthe(ji)*rpi/180.)
327cmod *              +pvm1(ji,jk)*sin(pthe(ji)*rpi/180.))
328cmod *              *cos(pthe(ji)*rpi/180.)/ztmst
329cmod     zdvdt(ji)=-(pum1(ji,jk)*cos(pthe(ji)*rpi/180.)
330cmod *              +pvm1(ji,jk)*sin(pthe(ji)*rpi/180.))
331cmod *              *sin(pthe(ji)*rpi/180.)/ztmst
332C
333         zdudt(ji)=zdudt(ji)*(zbet/(1.+zbet))
334         zdvdt(ji)=zdvdt(ji)*(zbet/(1.+zbet))
335      end if
336      pvom(ji,jk)=zdudt(ji)
337      pvol(ji,jk)=zdvdt(ji)
338      zust=pum1(ji,jk)+ztmst*zdudt(ji)
339      zvst=pvm1(ji,jk)+ztmst*zdvdt(ji)
340      zdis=0.5*(pum1(ji,jk)**2+pvm1(ji,jk)**2-zust**2-zvst**2)
341      zdedt(ji)=zdis/ztmst
342      zvidis(ji)=zvidis(ji)+zdis*zdelp
343c VENUS ATTENTION: CP VARIABLE
344      zdtdt(ji)=zdedt(ji)/rcpd
345c
346c  NO TENDENCIES ON TEMPERATURE .....
347c
348c  Instead of, pte(ji,jk)=zdtdt(ji), due to mechanical dissipation
349c
350      pte(ji,jk)=0.0
351
352      endif
353
354  523 continue
355  524 continue
356c
357c
358  501 continue
359
360      return
361      end
362
Note: See TracBrowser for help on using the repository browser.