source: trunk/LMDZ.VENUS/libf/phyvenus/orosetup.F @ 3567

Last change on this file since 3567 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: 16.2 KB
RevLine 
[3]1      SUBROUTINE orosetup
2     *         ( nlon   , nlev  , ktest
3     *         , kkcrit, kkcrith, kcrit, ksect , kkhlim
[2047]4     *         , kkenvh, kknu  , kknu2, kkbreak
[3]5     *         , paphm1, papm1 , pum1 , pvm1, ptm1, pgeom1, pstab, pstd
6     *         , prho  , pri   , ptau, pvph, ppsi, pzdep
7     *         , pulow , pvlow 
8     *         , ptheta, pgam, pmea, ppic, pval
9     *         , pnu  ,  pd1  ,  pd2  ,pdmod  )
10C
11c**** *gwsetup*
12c
13c     purpose.
14c     --------
15c     SET-UP THE ESSENTIAL PARAMETERS OF THE SSO DRAG SCHEME:
16C     DEPTH OF LOW WBLOCKED LAYER, LOW-LEVEL FLOW, BACKGROUND
17C     STRATIFICATION.....
18c
19c**   interface.
20c     ----------
21c          from *orodrag*
22c
23c        explicit arguments :
24c        --------------------
25c     ==== inputs ===
26c
27c nlon----input-I-Total number of horizontal points that get into physics
28c nlev----input-I-Number of vertical levels
29c ktest--input-I: Flags to indicate active points
30c
31c ptsphy--input-R-Time-step (s)
32c paphm1--input-R: pressure at model 1/2 layer
33c papm1---input-R: pressure at model layer
34c pgeom1--input-R: Altitude of layer above ground
35c VENUS ATTENTION: CP VARIABLE PSTAB CALCULE EN AMONT DES PARAMETRISATIONS
36c pstab-----R-: Brunt-Vaisala freq.^2 at 1/2 layers (input except klev+1)
37c ptm1, pum1, pvm1--R-: t, u and v
38c pmea----input-R-Mean Orography (m)
39C pstd----input-R-SSO standard deviation (m)
40c psig----input-R-SSO slope
41c pgam----input-R-SSO Anisotropy
42c pthe----input-R-SSO Angle
43c ppic----input-R-SSO Peacks elevation (m)
44c pval----input-R-SSO Valleys elevation (m)
45
46c     ==== outputs ===
47c pulow, pvlow -output-R: Low-level wind
48c kkcrit----I-: Security value for top of low level flow
49c kcrit-----I-: Critical level
50c ksect-----I-: Not used
51c kkhlim----I-: Not used
52c kkenvh----I-: Top of blocked flow layer
53c kknu------I-: Layer that sees mountain peacks
54c kknu2-----I-: Layer that sees mountain peacks above mountain mean
55c kknub-----I-: Layer that sees mountain mean above valleys
56c prho------R-: Density at 1/2 layers
57c pri-------R-: Background Richardson Number, Wind shear measured along GW stress
58c pvph------R-: Wind in  plan of GW stress, Half levels.
59c ppsi------R-: Angle between low level wind and SS0 main axis.
60c pd1-------R-| Compared the ratio of the stress
61c pd2-------R-| that is along the wind to that Normal to it.
62c               pdi define the plane of low level stress
63c               compared to the low level wind.
64c see p. 108 Lott & Miller (1997).                     
65c pdmod-----R-: Norme of pdi
66
67c     === local arrays ===
68c
69c zvpf------R-: Wind projected in the plan of the low-level stress.
70
71c     ==== outputs ===
72c
73c        implicit arguments :   none
74c        --------------------
75c
76c     method.
77c     -------
78c
79c
80c     externals.
81c     ----------
82c
83c
84c     reference.
85c     ----------
86c
87c        see ecmwf research department documentation of the "i.f.s."
88c
89c     author.
90c     -------
91c
92c     modifications.
93c     --------------
94c     f.lott  for the new-gwdrag scheme november 1993
95c
96c-----------------------------------------------------------------------
[101]97      use dimphy
[3]98      implicit none
99
100#include "YOMCST.h"
101#include "YOEGWD.h"
102
103c-----------------------------------------------------------------------
104c
105c*       0.1   arguments
106c              ---------
107c
108      integer nlon,nlev
109      integer kkcrit(nlon),kkcrith(nlon),kcrit(nlon),ksect(nlon),
110     *        kkhlim(nlon),ktest(nlon),kkenvh(nlon)
[2047]111      integer kkbreak(nlon)
[3]112c
113      real paphm1(nlon,klev+1),papm1(nlon,klev),pum1(nlon,klev),
114     *     pvm1(nlon,klev),ptm1(nlon,klev),pgeom1(nlon,klev),
115     *     prho(nlon,klev+1),pri(nlon,klev+1),pstab(nlon,klev+1),
116     *     ptau(nlon,klev+1),pvph(nlon,klev+1),ppsi(nlon,klev+1),
117     *     pzdep(nlon,klev)
118       real pulow(nlon),pvlow(nlon),ptheta(nlon),pgam(nlon),pnu(nlon),
119     *     pd1(nlon),pd2(nlon),pdmod(nlon)
120      real pstd(nlon),pmea(nlon),ppic(nlon),pval(nlon)
121c
122c-----------------------------------------------------------------------
123c
124c*       0.2   local arrays
125c              ------------
126c
127c
128      integer ilevh ,jl,jk,iii
129      real zcons1,zhgeo,zu,zphi
130      real zvt1,zvt2,zdwind,zwind,zdelp
131      real zstabm,zstabp,zrhom,zrhop
132      logical lo
133      logical ll1(klon,klev+1)
134      integer kknu(klon),kknu2(klon),kknub(klon),kknul(klon),
135     *        kentp(klon),ncount(klon) 
136c
137      real zhcrit(klon,klev),zvpf(klon,klev),
138     *     zdp(klon,klev)
139      real znorm(klon),zb(klon),zc(klon),
140     *      zulow(klon),zvlow(klon),znup(klon),znum(klon)
141
142c     ------------------------------------------------------------------
143c
144c*         1.    initialization
145c                --------------
146c
147c       PRINT *,' in orosetup'
148 100  continue
149c
150c     ------------------------------------------------------------------
151c
152c*         1.1   computational constants
153c                -----------------------
154c
155 110  continue
156c
157      ilevh =klev/3
158c
159      zcons1=1./rd
160c
161c     ------------------------------------------------------------------
162c
163c*         2.
164c                --------------
165c
166 200  continue
167c
168c     ------------------------------------------------------------------
169c
170c*         2.1     define low level wind, project winds in plane of
171c*                 low level wind, determine sector in which to take
172c*                 the variance and set indicator for critical levels.
173c
174c
175c
176      do 2001 jl=kidia,kfdia
177      kknu(jl)    =klev
178      kknu2(jl)   =klev
179      kknub(jl)   =klev
180      kknul(jl)   =klev
[2047]181      kkbreak(jl) =klev + 1
[3]182      pgam(jl) =max(pgam(jl),gtsec)
183      ll1(jl,klev+1)=.false.
184 2001 continue
185c
186c Ajouter une initialisation (L. Li, le 23fev99):
187c
188      do jk=klev,ilevh,-1
189      do jl=kidia,kfdia
190      ll1(jl,jk)= .false.
191      ENDDO
192      ENDDO
[2047]193
194c      VENUS: define break for subcritical stress
195c      ----------------------------
196      do jk=klev,ilevh,-1
197      do jl=kidia,kfdia
198      if(ktest(jl).eq.1) then
199      !zhgeo=pgeom1(jl,jk)/rg
200      !!zhcrit(jl,jk)=ppic(jl)
201      !zhcrit(jl,jk)=amin1(ppic(jl)-pmea(jl),pmea(jl)-pval(jl))
202      !ll1(jl,jk)=(zhgeo.gt.zhcrit(jl,jk))
203      !if(ll1(jl,jk) .neqv. ll1(jl,jk+1)) then
204      !  kkbreak(jl)=jk
205      !endif
206
207      !if (paphm1(jl,jk) .ge. 7.e6) kkbreak(jl)=jk
208      !kkbreak(jl) = klev - 2 ! gwd1103
209      !kkbreak(jl) = klev - 4 ! gwd1104
210      !kkbreak(jl) = klev - 3 ! gwd1105
211
212      endif     
213      enddo
214      enddo
215
[3]216c
217c*      define top of low level flow
218c       ----------------------------
219      do 2002 jk=klev,ilevh,-1
220      do 2003 jl=kidia,kfdia
221      if(ktest(jl).eq.1) then
222      lo=(paphm1(jl,jk)/paphm1(jl,klev+1)).ge.gsigcr
223      if(lo) then
224        kkcrit(jl)=jk
225      endif
226      zhcrit(jl,jk)=ppic(jl)-pval(jl)           
227      zhgeo=pgeom1(jl,jk)/rg
228      ll1(jl,jk)=(zhgeo.gt.zhcrit(jl,jk))
229C     if(ll1(jl,jk).xor.ll1(jl,jk+1)) then
230      if(ll1(jl,jk).neqv.ll1(jl,jk+1)) then
231        kknu(jl)=jk
232      endif
233      if(.not.ll1(jl,ilevh))kknu(jl)=ilevh
234      endif
235 2003 continue
236 2002 continue
237      do 2004 jk=klev,ilevh,-1
238      do 2005 jl=kidia,kfdia
239      if(ktest(jl).eq.1) then
[2047]240!      zhcrit(jl,jk)=ppic(jl)-pmea(jl)
[3]241      zhgeo=pgeom1(jl,jk)/rg
[2047]242!      ll1(jl,jk)=(zhgeo.gt.zhcrit(jl,jk))
243      ll1(jl,jk)=(zhgeo.gt.0.5*pstd(jl)) ! TN : do not consider outlier peaks
244!      ll1(jl,jk)=(zhgeo.gt.pstd(jl)) ! TN : do not consider outlier peaks
245!      ll1(jl,jk)=(zhgeo.gt.2*pstd(jl)) ! TN : do not consider outlier peaks
[3]246      if(ll1(jl,jk) .neqv. ll1(jl,jk+1)) then
247        kknu2(jl)=jk
248      endif
249      if(.not.ll1(jl,ilevh))kknu2(jl)=ilevh
250      endif
251 2005 continue
252 2004 continue
[2047]253
254!      do 2104 jk=klev,ilevh,-1
255!      do 2105 jl=kidia,kfdia
256!      if(ktest(jl).eq.1) then
257!      zhgeo=pgeom1(jl,jk)/rg
258!      ll1(jl,jk)=(zhgeo.gt.pstd(jl)) ! TN : do not consider outlier peaks
259!      if(ll1(jl,jk) .neqv. ll1(jl,jk+1)) then
260!        kknul(jl)=jk
261!      endif
262!      if(.not.ll1(jl,ilevh))kknul(jl)=ilevh
263!      endif
264! 2105 continue
265! 2104 continue
266
267
[3]268      do 2006 jk=klev,ilevh,-1
269      do 2007 jl=kidia,kfdia
270      if(ktest(jl).eq.1) then
271      zhcrit(jl,jk)=amin1(ppic(jl)-pmea(jl),pmea(jl)-pval(jl))
272      zhgeo=pgeom1(jl,jk)/rg
273      ll1(jl,jk)=(zhgeo.gt.zhcrit(jl,jk))
274c     if(ll1(jl,jk).xor.ll1(jl,jk+1)) then
275      if(ll1(jl,jk).neqv.ll1(jl,jk+1)) then
276        kknub(jl)=jk
277      endif
278      if(.not.ll1(jl,ilevh))kknub(jl)=ilevh
279      endif
280 2007 continue
281 2006 continue
282c
283      do 2010 jl=kidia,kfdia 
284      if(ktest(jl).eq.1) then
285      kknu(jl)=min(kknu(jl),nktopg)
286      kknu2(jl)=min(kknu2(jl),nktopg)
287      kknub(jl)=min(kknub(jl),nktopg)
[2047]288!      kknul(jl)=klev
[3]289      endif
290 2010 continue     
291c
292 210  continue
293c
294cc*     initialize various arrays
295c
296      do 2107 jl=kidia,kfdia
297      prho(jl,klev+1)  =0.0
298cym correction en attendant mieux
299      prho(jl,1)  =0.0     
300      pstab(jl,klev+1) =0.0
301      pstab(jl,1)      =0.0
302      pri(jl,klev+1)   =9999.0
303      ppsi(jl,klev+1)  =0.0
304      pri(jl,1)        =0.0
305      pvph(jl,1)       =0.0
306      pvph(jl,klev+1)  =0.0
307cym correction en attendant mieux
308cym      pvph(jl,klev)    =0.0
309      pulow(jl)        =0.0
310      pvlow(jl)        =0.0
311      zulow(jl)        =0.0
312      zvlow(jl)        =0.0
313      kkcrith(jl)      =klev
314      kkenvh(jl)       =klev
315      kentp(jl)        =klev
316      kcrit(jl)        =1
317      ncount(jl)       =0
318      ll1(jl,klev+1)   =.false.
319 2107 continue
320c
321c*     define flow density and stratification (rho and N2)
322c      at semi layers.
323c      -------------------------------------------------------
324c
325      do 223 jk=klev,2,-1
326      do 222 jl=kidia,kfdia
327      if(ktest(jl).eq.1) then
328        zdp(jl,jk)=papm1(jl,jk)-papm1(jl,jk-1)
329        prho(jl,jk)=2.*paphm1(jl,jk)*zcons1/(ptm1(jl,jk)+ptm1(jl,jk-1))
330      endif
331  222 continue
332  223 continue
333c      print*,"altitude(m)=",pgeom1(kfdia/2,:)
334c      print*,"pression(Pa)=",papm1(kfdia/2,:)
335c
336c********************************************************************
337c
338c*     define Low level flow (between ground and peacks-valleys)
339c      ---------------------------------------------------------
340      do 2115 jk=klev,ilevh,-1
341      do 2116 jl=kidia,kfdia
342      if(ktest(jl).eq.1)  then
343      if(jk.ge.kknu2(jl).and.jk.le.kknul(jl)) then
344        pulow(jl)=pulow(jl)+pum1(jl,jk)*(paphm1(jl,jk+1)-paphm1(jl,jk))
345        pvlow(jl)=pvlow(jl)+pvm1(jl,jk)*(paphm1(jl,jk+1)-paphm1(jl,jk))
346        pstab(jl,klev+1)=pstab(jl,klev+1)
347     c                   +pstab(jl,jk)*(paphm1(jl,jk+1)-paphm1(jl,jk))
348        prho(jl,klev+1)=prho(jl,klev+1)
349     c                   +prho(jl,jk)*(paphm1(jl,jk+1)-paphm1(jl,jk))
350      end if
351      endif
352 2116 continue
353 2115 continue
354      do 2110 jl=kidia,kfdia
355      if(ktest(jl).eq.1)  then
356      pulow(jl)=pulow(jl)/(paphm1(jl,kknul(jl)+1)-paphm1(jl,kknu2(jl)))
357      pvlow(jl)=pvlow(jl)/(paphm1(jl,kknul(jl)+1)-paphm1(jl,kknu2(jl)))
358      znorm(jl)=max(sqrt(pulow(jl)**2+pvlow(jl)**2),gvsec)
359      pvph(jl,klev+1)=znorm(jl)
360      pstab(jl,klev+1)=pstab(jl,klev+1)
361     c                /(paphm1(jl,kknul(jl)+1)-paphm1(jl,kknu2(jl)))
362      prho(jl,klev+1)=prho(jl,klev+1)
363     c                /(paphm1(jl,kknul(jl)+1)-paphm1(jl,kknu2(jl)))
364      endif
365 2110 continue
366
367c
368c*******  setup orography orientation relative to the low level
369C       wind and define parameters of the Anisotropic wave stress.
370c
371      do 2112 jl=kidia,kfdia
372      if(ktest(jl).eq.1)  then
373        lo=(pulow(jl).lt.gvsec).and.(pulow(jl).ge.-gvsec)
374        if(lo) then
375          zu=pulow(jl)+2.*gvsec
376        else
377          zu=pulow(jl)
378        endif
379        zphi=atan(pvlow(jl)/zu)
380        ppsi(jl,klev+1)=ptheta(jl)*rpi/180.-zphi
381        zb(jl)=1.-0.18*pgam(jl)-0.04*pgam(jl)**2
382        zc(jl)=0.48*pgam(jl)+0.3*pgam(jl)**2
383        pd1(jl)=zb(jl)-(zb(jl)-zc(jl))*(sin(ppsi(jl,klev+1))**2)
384        pd2(jl)=(zb(jl)-zc(jl))*sin(ppsi(jl,klev+1))
385     *                         *cos(ppsi(jl,klev+1))
386        pdmod(jl)=sqrt(pd1(jl)**2+pd2(jl)**2)
387      endif
388 2112 continue
389c
390c  ************ projet flow in plane of lowlevel stress *************
391C  ************ Find critical levels...                 *************
392c
393      do 213 jk=1,klev
394      do 212 jl=kidia,kfdia
395      if(ktest(jl).eq.1)  then
396        zvt1       =pulow(jl)*pum1(jl,jk)+pvlow(jl)*pvm1(jl,jk)
397        zvt2       =-pvlow(jl)*pum1(jl,jk)+pulow(jl)*pvm1(jl,jk)
398        zvpf(jl,jk)=(zvt1*pd1(jl)+zvt2*pd2(jl))/(znorm(jl)*pdmod(jl))
399      endif
400      ptau(jl,jk)  =0.0
401      pzdep(jl,jk) =0.0
402      ppsi(jl,jk)  =0.0
403      ll1(jl,jk)   =.false.
404  212 continue
405  213 continue
406      do 215 jk=2,klev
407      do 214 jl=kidia,kfdia
408      if(ktest(jl).eq.1) then
409        zdp(jl,jk)=papm1(jl,jk)-papm1(jl,jk-1)
410        pvph(jl,jk)=((paphm1(jl,jk)-papm1(jl,jk-1))*zvpf(jl,jk)+
411     *            (papm1(jl,jk)-paphm1(jl,jk))*zvpf(jl,jk-1))
412     *            /zdp(jl,jk)
413        if(pvph(jl,jk).lt.gvsec) then
414          pvph(jl,jk)=gvsec
415          kcrit(jl)=jk
416        endif
417      endif
418  214 continue
419  215 continue
420c
421c*         2.3     mean flow richardson number.
422c
423  230 continue
424c
425      do 232 jk=2,klev
426      do 231 jl=kidia,kfdia
427      if(ktest(jl).eq.1) then
428        zdwind=max(abs(zvpf(jl,jk)-zvpf(jl,jk-1)),gvsec)
429        pri(jl,jk)=pstab(jl,jk)*(zdp(jl,jk)
430     *          /(rg*prho(jl,jk)*zdwind))**2
431        pri(jl,jk)=max(pri(jl,jk),grcrit)
432      endif
433  231 continue
434  232 continue
435 
436c
437c
438c*      define top of 'envelope' layer
439c       ----------------------------
440
441      do 233 jl=kidia,kfdia
442      pnu (jl)=0.0
443      znum(jl)=0.0
444 233  continue
445     
446      do 234 jk=2,klev-1
447      do 234 jl=kidia,kfdia
448     
449      if(ktest(jl).eq.1) then
450       
451      if (jk.ge.kknu2(jl)) then
452         
453            znum(jl)=pnu(jl)
454            zwind=(pulow(jl)*pum1(jl,jk)+pvlow(jl)*pvm1(jl,jk))/
455     *            max(sqrt(pulow(jl)**2+pvlow(jl)**2),gvsec)
456            zwind=max(sqrt(zwind**2),gvsec)
457            zdelp=paphm1(jl,jk+1)-paphm1(jl,jk)
458            zstabm=sqrt(max(pstab(jl,jk  ),gssec))
459            zstabp=sqrt(max(pstab(jl,jk+1),gssec))
460            zrhom=prho(jl,jk  )
461            zrhop=prho(jl,jk+1)
462            pnu(jl) = pnu(jl) + (zdelp/rg)*
463     *            ((zstabp/zrhop+zstabm/zrhom)/2.)/zwind     
464            if((znum(jl).le.gfrcrit).and.(pnu(jl).gt.gfrcrit)
465     *                          .and.(kkenvh(jl).eq.klev))
466     *      kkenvh(jl)=jk
467     
468      endif   
469
470      endif
471     
472 234  continue
473     
474c  calculation of a dynamical mixing height for when the waves
475C  BREAK AT LOW LEVEL: The drag will be repartited over
476C  a depths that depends on waves vertical wavelength,
477C  not just between two adjacent model layers.
478c  of gravity waves:
479
480      do 235 jl=kidia,kfdia
481      znup(jl)=0.0
482      znum(jl)=0.0
483 235  continue
484
485      do 236 jk=klev-1,2,-1
486      do 236 jl=kidia,kfdia
487     
488      if(ktest(jl).eq.1) then
489
490            znum(jl)=znup(jl)
491            zwind=(pulow(jl)*pum1(jl,jk)+pvlow(jl)*pvm1(jl,jk))/
492     *            max(sqrt(pulow(jl)**2+pvlow(jl)**2),gvsec)
493            zwind=max(sqrt(zwind**2),gvsec)
494            zdelp=paphm1(jl,jk+1)-paphm1(jl,jk)
495            zstabm=sqrt(max(pstab(jl,jk  ),gssec))
496            zstabp=sqrt(max(pstab(jl,jk+1),gssec))
497            zrhom=prho(jl,jk  )
498            zrhop=prho(jl,jk+1)
499            znup(jl) = znup(jl) + (zdelp/rg)*
500     *            ((zstabp/zrhop+zstabm/zrhom)/2.)/zwind     
501            if((znum(jl).le.rpi/4.).and.(znup(jl).gt.rpi/4.)
502     *                          .and.(kkcrith(jl).eq.klev))
503     *      kkcrith(jl)=jk
504     
505      endif
506     
507 236  continue
508 
509      do 237 jl=kidia,kfdia
510      if(ktest(jl).eq.1) then
511      kkcrith(jl)=max0(kkcrith(jl),ilevh*2)
512      kkcrith(jl)=max0(kkcrith(jl),kknu(jl))
513      if(kcrit(jl).ge.kkcrith(jl))kcrit(jl)=1
514      endif
515 237  continue         
516c
517c     directional info for flow blocking *************************
518c
519      do 251 jk=1,klev   
520      do 252 jl=kidia,kfdia
521      if(ktest(jl).eq.1) then
522      lo=(pum1(jl,jk).lt.gvsec).and.(pum1(jl,jk).ge.-gvsec)
523      if(lo) then
524        zu=pum1(jl,jk)+2.*gvsec
525      else
526        zu=pum1(jl,jk)
527      endif
528       zphi=atan(pvm1(jl,jk)/zu)
529       ppsi(jl,jk)=ptheta(jl)*rpi/180.-zphi
530      endif
531 252  continue
532 251  continue
533
534c      forms the vertical 'leakiness' **************************
535
536      do 254  jk=ilevh,klev
537      do 253  jl=kidia,kfdia
538      if(ktest(jl).eq.1) then
539      pzdep(jl,jk)=0
540      if(jk.ge.kkenvh(jl).and.kkenvh(jl).ne.klev) then
541        pzdep(jl,jk)=(pgeom1(jl,kkenvh(jl)  )-pgeom1(jl,  jk))/
542     *               (pgeom1(jl,kkenvh(jl)  )-pgeom1(jl,klev))
543      end if
544      endif
545 253  continue
546 254  continue
547
548      return
549      end
550
551
Note: See TracBrowser for help on using the repository browser.