source: LMDZ4/trunk/libf/phylmd/wake.F @ 1042

Last change on this file since 1042 was 974, checked in by lmdzadmin, 16 years ago

Nouvelles versions vectorisees ; on garde versions scalaires; nom _scal
IM
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 66.8 KB
Line 
1      Subroutine WAKE (p,ph,ppi,dtime,sigd_con
2     :                ,te0,qe0,omgb
3     :                ,dtdwn,dqdwn,amdwn,amup,dta,dqa
4     :                ,wdtPBL,wdqPBL,udtPBL,udqPBL
5     o                ,deltatw,deltaqw,dth,hw,sigmaw,wape,fip,gfl
6     o                ,dtls,dqls
7     o                ,ktopw,omgbdth,dp_omgb,wdens
8     o                ,tu,qu
9     o                ,dtKE,dqKE
10     o                ,dtPBL,dqPBL
11     o                ,omg,dp_deltomg,spread
12     o                ,Cstar,d_deltat_gw
13     o                ,d_deltatw2,d_deltaqw2)
14
15***************************************************************
16*                                                             *
17* WAKE                                                        *
18*      retour a un Pupper fixe                                *
19*                                                             *
20* written by   :  GRANDPEIX Jean-Yves   09/03/2000            *
21* modified by :   ROEHRIG Romain        01/29/2007            *
22***************************************************************
23c
24      use dimphy
25      IMPLICIT none
26c============================================================================
27C
28C
29C   But : Decrire le comportement des poches froides apparaissant dans les
30C        grands systemes convectifs, et fournir l'energie disponible pour
31C        le declenchement de nouvelles colonnes convectives.
32C
33C   Variables d'etat : deltatw    : ecart de temperature wake-undisturbed area
34C                      deltaqw    : ecart d'humidite wake-undisturbed area
35C                      sigmaw     : fraction d'aire occupee par la poche.
36C
37C   Variable de sortie :
38c
39c                        wape : WAke Potential Energy
40c                        fip  : Front Incident Power (W/m2) - ALP
41c                        gfl  : Gust Front Length per unit area (m-1)
42C                        dtls : large scale temperature tendency due to wake
43C                        dqls : large scale humidity tendency due to wake
44C                        hw   : hauteur de la poche
45C                     dp_omgb : vertical gradient of large scale omega
46C                      omgbdth: flux of Delta_Theta transported by LS omega
47C                      dtKE   : differential heating (wake - unpertubed)
48C                      dqKE   : differential moistening (wake - unpertubed)
49C                      omg    : Delta_omg =vertical velocity diff. wake-undist. (Pa/s)
50C                 dp_deltomg  : vertical gradient of omg (s-1)
51C                     spread  : spreading term in dt_wake and dq_wake
52C                 deltatw     : updated temperature difference (T_w-T_u).
53C                 deltaqw     : updated humidity difference (q_w-q_u).
54C                 sigmaw      : updated wake fractional area.
55C                 d_deltat_gw : delta T tendency due to GW
56c
57C   Variables d'entree :
58c
59c                        aire : aire de la maille
60c                        te0  : temperature dans l'environnement  (K)
61C                        qe0  : humidite dans l'environnement     (kg/kg)
62C                        omgb : vitesse verticale moyenne sur la maille (Pa/s)
63C                        dtdwn: source de chaleur due aux descentes (K/s)
64C                        dqdwn: source d'humidite due aux descentes (kg/kg/s)
65C                        dta  : source de chaleur due courants satures et detrain  (K/s)
66C                        dqa  : source d'humidite due aux courants satures et detra (kg/kg/s)
67C                        amdwn: flux de masse total des descentes, par unite de
68C                                surface de la maille (kg/m2/s)
69C                        amup : flux de masse total des ascendances, par unite de
70C                                surface de la maille (kg/m2/s)
71C                        p    : pressions aux milieux des couches (Pa)
72C                        ph   : pressions aux interfaces (Pa)
73C                        ppi  : (p/p_0)**kapa (adim)
74C                        dtime: increment temporel (s)
75c
76C   Variables internes :
77c
78c                        rhow : masse volumique de la poche froide
79C                        rho  : environment density at P levels
80C                        rhoh : environment density at Ph levels
81C                        te   : environment temperature | may change within
82C                        qe   : environment humidity    | sub-time-stepping
83C                        the  : environment potential temperature
84C                        thu  : potential temperature in undisturbed area
85C                        tu   :  temperature  in undisturbed area
86C                        qu   : humidity in undisturbed area
87C                      dp_omgb: vertical gradient og LS omega
88C                      omgbw  : wake average vertical omega
89C                     dp_omgbw: vertical gradient of omgbw
90C                      omgbdq : flux of Delta_q transported by LS omega
91C                        dth  : potential temperature diff. wake-undist.
92C                        th1  : first pot. temp. for vertical advection (=thu)
93C                        th2  : second pot. temp. for vertical advection (=thw)
94C                        q1   : first humidity for vertical advection
95C                        q2   : second humidity for vertical advection
96C                     d_deltatw   : terme de redistribution pour deltatw
97C                     d_deltaqw   : terme de redistribution pour deltaqw
98C                      deltatw0   : deltatw initial
99C                      deltaqw0   : deltaqw initial
100C                      hw0    : hw initial
101C                      sigmaw0: sigmaw initial
102C                      amflux : horizontal mass flux through wake boundary
103C                      wdens  : number of wakes per unit area (3D) or per
104C                               unit length (2D)
105C                      Tgw    : 1 sur la période de onde de gravité
106c                      Cgw    : vitesse de propagation de onde de gravité
107c                      LL     : distance entre 2 poches
108
109c-------------------------------------------------------------------------
110c          Déclaration de variables
111c-------------------------------------------------------------------------
112
113#include "dimensions.h"
114#include "YOMCST.h"
115#include "cvthermo.h"
116#include "iniprint.h"
117
118c Arguments en entree
119c--------------------
120
121      REAL, dimension(klon,klev) :: p, ppi
122      REAL, dimension(klon,klev+1) :: ph, omgb
123      REAL dtime
124      REAL, dimension(klon,klev) :: te0,qe0
125      REAL, dimension(klon,klev) :: dtdwn, dqdwn
126      REAL, dimension(klon,klev) :: wdtPBL,wdqPBL
127      REAL, dimension(klon,klev) :: udtPBL,udqPBL
128      REAL, dimension(klon,klev) :: amdwn, amup
129      REAL, dimension(klon,klev) :: dta, dqa
130      REAL, dimension(klon) :: sigd_con
131
132c Sorties
133c--------
134
135      REAL, dimension(klon,klev) :: deltatw, deltaqw, dth
136      REAL, dimension(klon,klev) :: tu, qu
137      REAL, dimension(klon,klev) :: dtls, dqls
138      REAL, dimension(klon,klev) :: dtKE, dqKE
139      REAL, dimension(klon,klev) :: dtPBL, dqPBL
140      REAL, dimension(klon,klev) :: spread
141      REAL, dimension(klon,klev) :: d_deltatgw
142      REAL, dimension(klon,klev) :: d_deltatw2, d_deltaqw2
143      REAL, dimension(klon,klev+1) :: omgbdth, omg
144      REAL, dimension(klon,klev) :: dp_omgb, dp_deltomg
145      REAL, dimension(klon,klev) :: d_deltat_gw
146      REAL, dimension(klon) :: hw, sigmaw, wape, fip, gfl, Cstar
147      INTEGER, dimension(klon) :: ktopw
148
149c Variables internes
150c-------------------
151
152c Variables à fixer
153      REAL ALON
154      REAL coefgw
155      REAL :: wdens0, wdens
156      REAL stark
157      REAL alpk
158      REAL delta_t_min
159      REAL Pupper
160      INTEGER nsub
161      REAL dtimesub
162      REAL sigmad, hwmin
163cIM 080208
164      LOGICAL, dimension(klon) :: gwake
165
166c Variables de sauvegarde
167      REAL, dimension(klon,klev) :: deltatw0
168      REAL, dimension(klon,klev) :: deltaqw0
169      REAL, dimension(klon,klev) :: te, qe
170      REAL, dimension(klon) :: sigmaw0, sigmaw1
171
172c Variables pour les GW
173      REAL, DIMENSION(klon) :: LL
174      REAL, dimension(klon,klev) :: N2
175      REAL, dimension(klon,klev) :: Cgw
176      REAL, dimension(klon,klev) :: Tgw
177
178c Variables liées au calcul de hw
179      REAL, DIMENSION(klon) :: ptop_provis, ptop, ptop_new
180      REAL, DIMENSION(klon) :: sum_dth
181      REAL, DIMENSION(klon) :: dthmin
182      REAL, DIMENSION(klon) :: z, dz, hw0
183      INTEGER, DIMENSION(klon) :: ktop, kupper
184
185c Autres variables internes
186      INTEGER isubstep, k, i
187
188      REAL, DIMENSION(klon) :: sum_thu, sum_tu, sum_qu,sum_thvu
189      REAL, DIMENSION(klon) :: sum_dq, sum_rho
190      REAL, DIMENSION(klon) :: sum_dtdwn, sum_dqdwn
191      REAL, DIMENSION(klon) :: av_thu, av_tu, av_qu, av_thvu
192      REAL, DIMENSION(klon) :: av_dth, av_dq, av_rho
193      REAL, DIMENSION(klon) :: av_dtdwn, av_dqdwn
194
195      REAL, DIMENSION(klon,klev) :: rho, rhow
196      REAL, DIMENSION(klon,klev+1) :: rhoh
197      REAL, DIMENSION(klon,klev) :: rhow_moyen
198      REAL, DIMENSION(klon,klev) :: zh
199      REAL, DIMENSION(klon,klev+1) :: zhh
200      REAL, DIMENSION(klon,klev) :: epaisseur1, epaisseur2
201
202      REAL, DIMENSION(klon,klev) :: the, thu
203
204      REAL, DIMENSION(klon,klev) :: d_deltatw, d_deltaqw
205
206      REAL, DIMENSION(klon,klev+1) :: omgbw
207      REAL, DIMENSION(klon) :: omgtop
208      REAL, DIMENSION(klon,klev) :: dp_omgbw
209      REAL, DIMENSION(klon) :: ztop, dztop
210      REAL, DIMENSION(klon,klev) :: alpha_up
211     
212      REAL, dimension(klon) :: RRe1, RRe2
213      REAL :: RRd1, RRd2
214      REAL, DIMENSION(klon,klev) :: Th1, Th2, q1, q2
215      REAL, DIMENSION(klon,klev) :: D_Th1, D_Th2, D_dth
216      REAL, DIMENSION(klon,klev) :: D_q1, D_q2, D_dq
217      REAL, DIMENSION(klon,klev) :: omgbdq
218
219      REAL, dimension(klon) :: ff, gg
220      REAL, dimension(klon) :: wape2, Cstar2, heff
221
222      REAL, DIMENSION(klon,klev) :: Crep
223      REAL Crep_upper, Crep_sol
224
225C-------------------------------------------------------------------------
226c         Initialisations
227c-------------------------------------------------------------------------
228
229c      print*, 'wake initialisations'
230
231c   Essais d'initialisation avec sigmaw = 0.02 et hw = 10.
232c-------------------------------------------------------------------------
233
234      DATA sigmad, hwmin /.02,10./
235
236C Longueur de maille (en m)
237c-------------------------------------------------------------------------
238
239c      ALON = 3.e5
240      ALON = 1.e6
241
242
243C Configuration de coefgw,stark,wdens (22/02/06 by YU Jingmei)
244c
245c      coefgw : Coefficient pour les ondes de gravité
246c       stark : Coefficient k dans Cstar=k*sqrt(2*WAPE)
247c       wdens : Densité de poche froide par maille
248c-------------------------------------------------------------------------
249
250      coefgw=10
251c      coefgw=1
252c      wdens0 = 1.0/(alon**2)   
253      wdens = 1.0/(alon**2)       
254      stark = 0.50
255cCRtest
256      alpk=0.1
257c      alpk = 1.0
258c      alpk = 0.5
259c      alpk = 0.05
260      Crep_upper=0.9
261      Crep_sol=1.0
262
263
264C Minimum value for |T_wake - T_undist|. Used for wake top definition
265c-------------------------------------------------------------------------
266
267      delta_t_min = 0.2
268
269C 1. - Save initial values and initialize tendencies
270C --------------------------------------------------
271
272      DO k=1,klev
273      DO i=1, klon
274        deltatw0(i,k) = deltatw(i,k)
275        deltaqw0(i,k)= deltaqw(i,k)
276        te(i,k) = te0(i,k)
277        qe(i,k) = qe0(i,k)
278        dtls(i,k) = 0.
279        dqls(i,k) = 0.
280        d_deltat_gw(i,k)=0.
281!IM 060508 beg
282        d_deltatw2(i,k)=0.
283        d_deltaqw2(i,k)=0.
284!IM 060508 end
285      ENDDO
286      ENDDO
287c      sigmaw1=sigmaw
288c      IF (sigd_con.GT.sigmaw1) THEN
289c      print*, 'sigmaw,sigd_con', sigmaw, sigd_con
290c      ENDIF
291      DO i=1, klon
292cc      sigmaw(i) = amax1(sigmaw(i),sigd_con(i))
293      sigmaw(i) = amax1(sigmaw(i),sigmad)
294      sigmaw(i) = amin1(sigmaw(i),0.99)
295      sigmaw0(i) = sigmaw(i)
296      ENDDO
297C
298C
299C 2. - Prognostic part
300C --------------------
301C
302C
303C 2.1 - Undisturbed area and Wake integrals
304C ---------------------------------------------------------
305
306      DO i=1, klon
307      z(i) = 0.
308      ktop(i)=0
309      kupper(i) = 0
310      sum_thu(i) = 0.
311      sum_tu(i) = 0.
312      sum_qu(i) = 0.
313      sum_thvu(i) = 0.
314      sum_dth(i) = 0.
315      sum_dq(i) = 0.
316      sum_rho(i) = 0.
317      sum_dtdwn(i) = 0.
318      sum_dqdwn(i) = 0.
319
320      av_thu(i) = 0.
321      av_tu(i) =0.
322      av_qu(i) =0.
323      av_thvu(i) = 0.
324      av_dth(i) = 0.
325      av_dq(i) = 0.
326      av_rho(i) =0.
327      av_dtdwn(i) =0.
328      av_dqdwn(i) = 0.
329      ENDDO
330c
331c Distance between wakes
332       DO i = 1,klon
333        LL(i) = (1-sqrt(sigmaw(i)))/sqrt(wdens)
334       ENDDO
335C Potential temperatures and humidity
336c----------------------------------------------------------
337      DO k =1,klev
338       DO i=1, klon
339        rho(i,k) = p(i,k)/(rd*te(i,k))
340        IF(k .eq. 1) THEN
341          rhoh(i,k) = ph(i,k)/(rd*te(i,k))
342          zhh(i,k)=0
343        ELSE
344          rhoh(i,k) = ph(i,k)*2./(rd*(te(i,k)+te(i,k-1)))
345          zhh(i,k)=(ph(i,k)-ph(i,k-1))/(-rhoh(i,k)*RG)+zhh(i,k-1)
346        ENDIF
347        the(i,k) = te(i,k)/ppi(i,k)
348        thu(i,k) = (te(i,k) - deltatw(i,k)*sigmaw(i))/ppi(i,k)
349        tu(i,k) = te(i,k) - deltatw(i,k)*sigmaw(i)
350        qu(i,k)  =  qe(i,k) - deltaqw(i,k)*sigmaw(i)
351        rhow(i,k) = p(i,k)/(rd*(te(i,k)+deltatw(i,k)))
352        dth(i,k) = deltatw(i,k)/ppi(i,k)
353       ENDDO
354      ENDDO
355       
356      DO k = 1, klev-1
357      DO i=1, klon
358        IF(k.eq.1) THEN
359          N2(i,k)=0
360        ELSE
361          N2(i,k)=amax1(0.,-RG**2/the(i,k)*rho(i,k)*(the(i,k+1)-
362     $            the(i,k-1))/(p(i,k+1)-p(i,k-1)))
363        ENDIF
364        ZH(i,k)=(zhh(i,k)+zhh(i,k+1))/2
365
366        Cgw(i,k)=sqrt(N2(i,k))*ZH(i,k)
367        Tgw(i,k)=coefgw*Cgw(i,k)/LL(i)
368      ENDDO
369      ENDDO
370
371      DO i=1, klon
372      N2(i,klev)=0
373      ZH(i,klev)=0
374      Cgw(i,klev)=0
375      Tgw(i,klev)=0
376      ENDDO
377
378c  Calcul de la masse volumique moyenne de la colonne   (bdlmd)
379c-----------------------------------------------------------------
380
381      DO k=1,klev
382       DO i=1, klon
383        epaisseur1(i,k)=0.
384        epaisseur2(i,k)=0.
385       ENDDO
386      ENDDO
387
388      DO i=1, klon
389      epaisseur1(i,1)= -(ph(i,2)-ph(i,1))/(rho(i,1)*rg)+1.
390      epaisseur2(i,1)= -(ph(i,2)-ph(i,1))/(rho(i,1)*rg)+1.
391      rhow_moyen(i,1) = rhow(i,1)
392      ENDDO
393
394      DO k = 2, klev
395      DO i=1, klon
396        epaisseur1(i,k)= -(ph(i,k+1)-ph(i,k))/(rho(i,k)*rg) +1.
397        epaisseur2(i,k)=epaisseur2(i,k-1)+epaisseur1(i,k)
398        rhow_moyen(i,k) = (rhow_moyen(i,k-1)*epaisseur2(i,k-1)+
399     $                 rhow(i,k)*epaisseur1(i,k))/epaisseur2(i,k)
400      ENDDO
401      ENDDO
402
403C
404C Choose an integration bound well above wake top
405c-----------------------------------------------------------------
406c
407C       Pupper = 50000.  ! melting level
408       Pupper = 60000.
409c       Pupper = 80000.  ! essais pour case_e
410C
411C    Determine Wake top pressure (Ptop) from buoyancy integral
412C    --------------------------------------------------------
413c
414c-1/ Pressure of the level where dth becomes less than delta_t_min.
415
416      DO i=1,klon
417      ptop_provis(i)=ph(i,1)
418      ENDDO
419      DO k= 2,klev
420      DO i=1,klon
421c
422cIM v3JYG; ptop_provis(i).LT. ph(i,1)
423c
424        IF (dth(i,k) .GT. -delta_t_min .and.
425     $      dth(i,k-1).LT. -delta_t_min .and.
426     $      ptop_provis(i).EQ. ph(i,1)) THEN
427          ptop_provis(i) = ((dth(i,k)+delta_t_min)*p(i,k-1)
428     $          - (dth(i,k-1)+delta_t_min)*p(i,k)) /
429     $          (dth(i,k) - dth(i,k-1))
430        ENDIF
431      ENDDO
432      ENDDO
433
434c-2/ dth integral
435
436      DO i=1,klon
437      sum_dth(i) = 0.
438      dthmin(i) = -delta_t_min
439      z(i) = 0.
440      ENDDO
441
442      DO k = 1,klev
443      DO i=1,klon
444        dz(i) = -(amax1(ph(i,k+1),ptop_provis(i))-Ph(i,k))/(rho(i,k)*rg)
445        IF (dz(i) .gt. 0) THEN
446          z(i) = z(i)+dz(i)
447          sum_dth(i) = sum_dth(i) + dth(i,k)*dz(i)
448          dthmin(i) = amin1(dthmin(i),dth(i,k))
449        ENDIF
450      ENDDO
451      ENDDO
452
453c-3/ height of triangle with area= sum_dth and base = dthmin
454
455      DO i=1,klon
456      hw0(i) = 2.*sum_dth(i)/amin1(dthmin(i),-0.5)
457      hw0(i) = amax1(hwmin,hw0(i))
458      ENDDO
459
460c-4/ now, get Ptop
461
462      DO i=1,klon
463      z(i) = 0.
464      ptop(i) = ph(i,1)
465      ENDDO
466
467      DO k = 1,klev
468      DO i=1,klon
469        dz(i) = amin1(-(ph(i,k+1)-ph(i,k))/(rho(i,k)*rg),hw0(i)-z(i))
470        IF (dz(i) .gt. 0) THEN
471         z(i) = z(i)+dz(i)
472         ptop(i) = ph(i,k)-rho(i,k)*rg*dz(i)
473        ENDIF
474      ENDDO
475      ENDDO
476
477
478C-5/ Determination de ktop et kupper
479
480      DO k=klev,1,-1
481      DO i=1,klon
482        IF (ph(i,k+1) .lt. ptop(i)) ktop(i)=k
483        IF (ph(i,k+1) .lt. pupper) kupper(i)=k
484      ENDDO
485      ENDDO
486
487c-6/ Correct ktop and ptop
488
489      DO i = 1,klon
490        ptop_new(i)=ptop(i)
491      ENDDO
492      DO k= klev,2,-1
493      DO i=1,klon
494        IF (k .LE. ktop(i) .and.
495     $      ptop_new(i) .EQ. ptop(i) .and.
496     $      dth(i,k) .GT. -delta_t_min .and.
497     $      dth(i,k-1).LT. -delta_t_min) THEN
498          ptop_new(i) = ((dth(i,k)+delta_t_min)*p(i,k-1)
499     $          - (dth(i,k-1)+delta_t_min)*p(i,k)) /
500     $          (dth(i,k) - dth(i,k-1))
501        ENDIF
502      ENDDO
503      ENDDO
504
505      DO i=1,klon
506        ptop(i) = ptop_new(i)
507      ENDDO
508
509      DO k=klev,1,-1
510      DO i=1,klon
511        IF (ph(i,k+1) .lt. ptop(i)) ktop(i)=k
512      ENDDO
513      ENDDO
514c
515c-5/ Set deltatw & deltaqw to 0 above kupper
516c
517      DO k = 1,klev
518      DO i=1,klon
519       IF (k.GE. kupper(i)) THEN
520        deltatw(i,k) = 0.
521        deltaqw(i,k) = 0.
522       ENDIF
523      ENDDO
524      ENDDO
525c
526C
527C Vertical gradient of LS omega
528C
529      DO k = 1,klev
530      DO i=1,klon
531       IF (k.LE. kupper(i)) THEN
532        dp_omgb(i,k) = (omgb(i,k+1) - omgb(i,k))/(ph(i,k+1)-ph(i,k))
533       ENDIF
534      ENDDO
535      ENDDO
536C
537C Integrals (and wake top level number)
538C --------------------------------------
539C
540C Initialize sum_thvu to 1st level virt. pot. temp.
541
542      DO i=1,klon
543      z(i) = 1.
544      dz(i) = 1.
545      sum_thvu(i) =  thu(i,1)*(1.+eps*qu(i,1))*dz(i)
546      sum_dth(i) = 0.
547      ENDDO
548
549      DO k = 1,klev
550      DO i=1,klon
551        dz(i) = -(amax1(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*rg)
552        IF (dz(i) .GT. 0) THEN
553         z(i) = z(i)+dz(i)
554         sum_thu(i) = sum_thu(i) + thu(i,k)*dz(i)
555         sum_tu(i) = sum_tu(i) + tu(i,k)*dz(i)
556         sum_qu(i) = sum_qu(i) + qu(i,k)*dz(i)
557         sum_thvu(i) = sum_thvu(i) + thu(i,k)*(1.+eps*qu(i,k))*dz(i)
558         sum_dth(i) = sum_dth(i) + dth(i,k)*dz(i)
559         sum_dq(i) = sum_dq(i) + deltaqw(i,k)*dz(i)
560         sum_rho(i) = sum_rho(i) + rhow(i,k)*dz(i)
561         sum_dtdwn(i) = sum_dtdwn(i) + dtdwn(i,k)*dz(i)
562         sum_dqdwn(i) = sum_dqdwn(i) + dqdwn(i,k)*dz(i)
563        ENDIF
564      ENDDO
565      ENDDO
566c
567      DO i=1,klon
568        hw0(i) = z(i)
569      ENDDO
570c
571C
572C 2.1 - WAPE and mean forcing computation
573C ---------------------------------------
574C
575C ---------------------------------------
576C
577C Means
578
579      DO i=1,klon
580      av_thu(i) = sum_thu(i)/hw0(i)
581      av_tu(i) = sum_tu(i)/hw0(i)
582      av_qu(i) = sum_qu(i)/hw0(i)
583      av_thvu(i) = sum_thvu(i)/hw0(i)
584c      av_thve = sum_thve/hw0
585      av_dth(i) = sum_dth(i)/hw0(i)
586      av_dq(i) = sum_dq(i)/hw0(i)
587      av_rho(i) = sum_rho(i)/hw0(i)
588      av_dtdwn(i) = sum_dtdwn(i)/hw0(i)
589      av_dqdwn(i) = sum_dqdwn(i)/hw0(i)
590
591      wape(i) = - rg*hw0(i)*(av_dth(i)
592     $     + eps*(av_thu(i)*av_dq(i)+av_dth(i)*av_qu(i)+av_dth(i)*
593     $     av_dq(i) ))/av_thvu(i)
594      ENDDO
595C
596C 2.2 Prognostic variable update
597C ------------------------------
598C
599C Filter out bad wakes
600
601      DO k = 1,klev
602       DO i=1,klon
603        IF ( wape(i) .LT. 0.) THEN
604          deltatw(i,k) = 0.
605          deltaqw(i,k) = 0.
606          dth(i,k) = 0.
607        ENDIF
608       ENDDO
609      ENDDO
610c
611      DO i=1,klon
612      IF ( wape(i) .LT. 0.) THEN
613        wape(i) = 0.
614        Cstar(i) = 0.
615        hw(i) = hwmin
616        sigmaw(i) = amax1(sigmad,sigd_con(i))
617        fip(i) = 0.
618        gwake(i) = .FALSE.
619      ELSE
620        Cstar(i) = stark*sqrt(2.*wape(i))
621        gwake(i) = .TRUE.
622      ENDIF
623      ENDDO
624c
625C
626CC -----------------------------------------------------------------
627C    Sub-time-stepping
628C    -----------------
629C
630      nsub=10
631      dtimesub=dtime/nsub
632c
633c------------------------------------------------------------
634      DO isubstep = 1,nsub
635c------------------------------------------------------------
636      DO i=1,klon
637        gfl(i) = 2.*sqrt(3.14*wdens*sigmaw(i))
638      ENDDO
639      DO i=1,klon
640        sigmaw(i) =sigmaw(i) + gfl(i)*Cstar(i)*dtimesub
641        sigmaw(i) =amin1(sigmaw(i),0.99)     !!!!!!!!
642c        wdens = wdens0/(10.*sigmaw)
643c        sigmaw =max(sigmaw,sigd_con)
644c        sigmaw =max(sigmaw,sigmad)
645      ENDDO
646C
647C
648c calcul de la difference de vitesse verticale poche - zone non perturbee
649cIM 060208 differences par rapport au code initial; init. a 0 dp_deltomg
650cIM 060208 et omg sur les niveaux de 1 a klev+1, alors que avant l'on definit
651cIM 060208 au niveau k=1..?
652      dp_deltomg(1:klon,1:klev)=0.
653      DO k= 1,klev+1
654      DO i = 1,klon
655        omg(i,k)=0.
656      ENDDO
657      ENDDO
658c
659      DO i=1,klon
660        z(i)= 0.
661        omg(i,1) = 0.
662        dp_deltomg(i,1) = -(gfl(i)*Cstar(i))/(sigmaw(i) * (1-sigmaw(i)))
663      ENDDO
664c
665      DO k= 2,klev
666      DO i = 1,klon
667       IF( k .LE. ktop(i)) THEN
668          dz(i) = -(ph(i,k)-ph(i,k-1))/(rho(i,k-1)*rg)
669          z(i) = z(i)+dz(i)
670          dp_deltomg(i,k)= dp_deltomg(i,1)
671          omg(i,k)= dp_deltomg(i,1)*z(i)
672       ENDIF
673      ENDDO
674      ENDDO
675c
676      DO i = 1,klon
677        dztop(i)=-(ptop(i)-ph(i,ktop(i)))/(rho(i,ktop(i))*rg)
678        ztop(i) = z(i)+dztop(i)
679        omgtop(i)=dp_deltomg(i,1)*ztop(i)
680      ENDDO
681c
682c        -----------------
683c        From m/s to Pa/s
684c        -----------------
685c
686       DO i=1,klon
687        omgtop(i) = -rho(i,ktop(i))*rg*omgtop(i)
688        dp_deltomg(i,1) = omgtop(i)/(ptop(i)-ph(i,1))
689       ENDDO
690c
691      DO k= 1,klev
692      DO i = 1,klon
693       IF( k .LE. ktop(i)) THEN
694          omg(i,k) = - rho(i,k)*rg*omg(i,k)
695          dp_deltomg(i,k) = dp_deltomg(i,1)
696       ENDIF
697      ENDDO
698      ENDDO
699c
700c   raccordement lineaire de omg de ptop a pupper
701
702      DO i=1,klon
703      IF (kupper(i) .GT. ktop(i)) THEN
704        omg(i,kupper(i)+1) = - Rg*amdwn(i,kupper(i)+1)/sigmaw(i)
705     $                + Rg*amup(i,kupper(i)+1)/(1.-sigmaw(i))
706        dp_deltomg(i,kupper(i)) = (omgtop(i)-omg(i,kupper(i)+1))/
707     $                     (ptop(i)-pupper)
708      ENDIF
709      ENDDO
710c
711      DO k= 1,klev
712      DO i = 1,klon
713       IF( k .GT. ktop(i) .AND. k .LE. kupper(i)) THEN
714          dp_deltomg(i,k) = dp_deltomg(i,kupper(i))
715          omg(i,k) = omgtop(i)+(ph(i,k)-ptop(i))*dp_deltomg(i,kupper(i))
716       ENDIF
717      ENDDO
718      ENDDO
719c
720c--    Compute wake average vertical velocity omgbw
721c
722c
723      DO k = 1,klev+1
724      DO i=1,klon
725        omgbw(i,k) = omgb(i,k)+(1.-sigmaw(i))*omg(i,k)
726      ENDDO
727      ENDDO
728c--    and its vertical gradient dp_omgbw
729c
730      DO k = 1,klev
731      DO i=1,klon
732        dp_omgbw(i,k) = (omgbw(i,k+1)-omgbw(i,k))/(ph(i,k+1)-ph(i,k))
733      ENDDO
734      ENDDO
735C
736c--    Upstream coefficients for omgb velocity
737c--    (alpha_up(k) is the coefficient of the value at level k)
738c--    (1-alpha_up(k) is the coefficient of the value at level k-1)
739      DO k = 1,klev
740      DO i=1,klon
741       alpha_up(i,k) = 0.
742       IF (omgb(i,k) .GT. 0.) alpha_up(i,k) = 1.
743      ENDDO
744      ENDDO
745
746c  Matrix expressing [The,deltatw] from  [Th1,Th2]
747
748      DO i=1,klon
749        RRe1(i) = 1.-sigmaw(i)
750        RRe2(i) = sigmaw(i)
751      ENDDO
752      RRd1 = -1.
753      RRd2 = 1.
754c
755c--    Get [Th1,Th2], dth and [q1,q2]
756c
757      DO k= 1,klev
758      DO i = 1,klon
759       IF(k .LE. kupper(i)+1) THEN
760        dth(i,k) = deltatw(i,k)/ppi(i,k)
761        Th1(i,k) = the(i,k) - sigmaw(i)     *dth(i,k)   ! undisturbed area
762        Th2(i,k) = the(i,k) + (1.-sigmaw(i))*dth(i,k)   ! wake
763        q1(i,k) = qe(i,k) - sigmaw(i)     *deltaqw(i,k) ! undisturbed area
764        q2(i,k) = qe(i,k) + (1.-sigmaw(i))*deltaqw(i,k) ! wake
765       ENDIF
766      ENDDO
767      ENDDO
768
769      DO i=1,klon
770       D_Th1(i,1) = 0.
771       D_Th2(i,1) = 0.
772       D_dth(i,1) = 0.
773       D_q1(i,1) = 0.
774       D_q2(i,1) = 0.
775       D_dq(i,1) = 0.
776      ENDDO
777
778      DO k= 2,klev
779      DO i = 1,klon
780       IF(k .LE. kupper(i)+1) THEN
781        D_Th1(i,k) = Th1(i,k-1)-Th1(i,k)
782        D_Th2(i,k) = Th2(i,k-1)-Th2(i,k)
783        D_dth(i,k) = dth(i,k-1)-dth(i,k)
784        D_q1(i,k) = q1(i,k-1)-q1(i,k)
785        D_q2(i,k) = q2(i,k-1)-q2(i,k)
786        D_dq(i,k) = deltaqw(i,k-1)-deltaqw(i,k)
787       ENDIF
788      ENDDO
789      ENDDO
790
791      DO i=1,klon
792        omgbdth(i,1) = 0.
793        omgbdq(i,1) = 0.
794      ENDDO
795
796      DO k= 2,klev
797      DO i = 1,klon
798       IF(k .LE. kupper(i)+1) THEN  !   loop on interfaces
799        omgbdth(i,k) = omgb(i,k)*(    dth(i,k-1) -     dth(i,k))
800        omgbdq(i,k)  = omgb(i,k)*(deltaqw(i,k-1) - deltaqw(i,k))
801       ENDIF
802      ENDDO
803      ENDDO
804c
805c-----------------------------------------------------------------
806      DO k= 1,klev
807      DO i = 1,klon
808       IF(k .LE. kupper(i)-1) THEN
809c-----------------------------------------------------------------
810c
811c   Compute redistribution (advective) term
812c
813         d_deltatw(i,k) =
814     $             dtimesub/(Ph(i,k)-Ph(i,k+1))*(
815     $       RRd1*omg(i,k  )*sigmaw(i)     *D_Th1(i,k)
816     $      -RRd2*omg(i,k+1)*(1.-sigmaw(i))*D_Th2(i,k+1)
817     $      -(1.-alpha_up(i,k))*omgbdth(i,k) - alpha_up(i,k+1)*
818     $      omgbdth(i,k+1))*ppi(i,k)
819c         print*,'d_deltatw=',d_deltatw(i,k)
820c
821         d_deltaqw(i,k) =
822     $             dtimesub/(Ph(i,k)-Ph(i,k+1))*(
823     $       RRd1*omg(i,k  )*sigmaw(i)     *D_q1(i,k)
824     $      -RRd2*omg(i,k+1)*(1.-sigmaw(i))*D_q2(i,k+1)
825     $      -(1.-alpha_up(i,k))*omgbdq(i,k) - alpha_up(i,k+1)*
826     $      omgbdq(i,k+1))
827c         print*,'d_deltaqw=',d_deltaqw(i,k)
828c
829c   and increment large scale tendencies
830c
831         dtls(i,k) = dtls(i,k) +
832     $               dtimesub*(
833     $        ( RRe1(i)*omg(i,k  )*sigmaw(i)     *D_Th1(i,k)
834     $         -RRe2(i)*omg(i,k+1)*(1.-sigmaw(i))*D_Th2(i,k+1) )
835     $               /(Ph(i,k)-Ph(i,k+1))
836     $         -sigmaw(i)*(1.-sigmaw(i))*dth(i,k)*dp_deltomg(i,k)
837     $                      )*ppi(i,k)
838c         print*,'dtls=',dtls(i,k)
839c
840         dqls(i,k) = dqls(i,k) +
841     $               dtimesub*(
842     $        ( RRe1(i)*omg(i,k  )*sigmaw(i)     *D_q1(i,k)
843     $         -RRe2(i)*omg(i,k+1)*(1.-sigmaw(i))*D_q2(i,k+1) )
844     $               /(Ph(i,k)-Ph(i,k+1))
845     $         -sigmaw(i)*(1.-sigmaw(i))*deltaqw(i,k)*dp_deltomg(i,k)
846     $                      )
847c         print*,'dqls=',dqls(k)
848       ENDIF
849c-------------------------------------------------------------------
850      ENDDO
851      ENDDO
852c------------------------------------------------------------------
853C
854C   Increment state variables
855
856      DO k= 1,klev
857      DO i = 1,klon
858       IF(k .LE. kupper(i)-1) THEN
859c
860c Coefficient de répartition
861
862        Crep(i,k)=Crep_sol*(ph(i,kupper(i))-ph(i,k))/(ph(i,kupper(i))
863     $          -ph(i,1))
864        Crep(i,k)=Crep(i,k)+Crep_upper*(ph(i,1)-ph(i,k))/(p(i,1)-
865     $          ph(i,kupper(i)))
866       
867
868c Reintroduce compensating subsidence term.
869
870c        dtKE(k)=(dtdwn(k)*Crep(k))/sigmaw
871c        dtKE(k)=dtKE(k)-(dtdwn(k)*(1-Crep(k))+dta(k))
872c     .                   /(1-sigmaw)
873c        dqKE(k)=(dqdwn(k)*Crep(k))/sigmaw
874c        dqKE(k)=dqKE(k)-(dqdwn(k)*(1-Crep(k))+dqa(k))
875c     .                   /(1-sigmaw)
876c
877c        dtKE(k)=(dtdwn(k)*Crep(k)+(1-Crep(k))*dta(k))/sigmaw
878c        dtKE(k)=dtKE(k)-(dtdwn(k)*(1-Crep(k))+dta(k)*Crep(k))
879c     .                   /(1-sigmaw)
880c        dqKE(k)=(dqdwn(k)*Crep(k)+(1-Crep(k))*dqa(k))/sigmaw
881c        dqKE(k)=dqKE(k)-(dqdwn(k)*(1-Crep(k))+dqa(k)*Crep(k))
882c     .                   /(1-sigmaw)
883
884        dtKE(i,k)=(dtdwn(i,k)/sigmaw(i) - dta(i,k)/(1.-sigmaw(i)))
885        dqKE(i,k)=(dqdwn(i,k)/sigmaw(i) - dqa(i,k)/(1.-sigmaw(i)))
886c        print*,'dtKE=',dtKE(k)
887c        print*,'dqKE=',dqKE(k)
888c
889        dtPBL(i,k)=(wdtPBL(i,k)/sigmaw(i) - udtPBL(i,k)/(1.-sigmaw(i)))
890        dqPBL(i,k)=(wdqPBL(i,k)/sigmaw(i) - udqPBL(i,k)/(1.-sigmaw(i)))
891c
892        spread(i,k) = (1.-sigmaw(i))*dp_deltomg(i,k)+gfl(i)*Cstar(i)/
893     $  sigmaw(i)
894
895
896c ajout d'un effet onde de gravité -Tgw(k)*deltatw(k) 03/02/06 YU Jingmei
897
898        d_deltat_gw(i,k)=d_deltat_gw(i,k)-Tgw(i,k)*deltatw(i,k)*
899     $  dtimesub
900        ff(i)=d_deltatw(i,k)/dtimesub
901
902c Sans GW
903c
904c        deltatw(k)=deltatw(k)+dtimesub*(ff+dtKE(k)-spread(k)*deltatw(k))
905c
906c GW formule 1
907c
908c        deltatw(k) = deltatw(k)+dtimesub*
909c     $         (ff+dtKE(k) - spread(k)*deltatw(k)-Tgw(k)*deltatw(k))
910c
911c GW formule 2
912
913        IF (dtimesub*Tgw(i,k).lt.1.e-10) THEN
914          deltatw(i,k) = deltatw(i,k)+dtimesub*
915     $          (ff(i)+dtKE(i,k)+dtPBL(i,k)
916     $          - spread(i,k)*deltatw(i,k)-Tgw(i,k)*deltatw(i,k))
917        ELSE
918           deltatw(i,k) = deltatw(i,k)+1/Tgw(i,k)*(1-exp(-dtimesub*
919     $          Tgw(i,k)))*
920     $          (ff(i)+dtKE(i,k)+dtPBL(i,k)
921     $          - spread(i,k)*deltatw(i,k)-Tgw(i,k)*deltatw(i,k))
922        ENDIF
923   
924        dth(i,k) = deltatw(i,k)/ppi(i,k)
925
926        gg(i)=d_deltaqw(i,k)/dtimesub
927
928       deltaqw(i,k) = deltaqw(i,k) +
929     $         dtimesub*(gg(i)+ dqKE(i,k)+dqPBL(i,k) - spread(i,k)*
930     $         deltaqw(i,k))
931
932       d_deltatw2(i,k)=d_deltatw2(i,k)+d_deltatw(i,k)
933       d_deltaqw2(i,k)=d_deltaqw2(i,k)+d_deltaqw(i,k)
934       ENDIF
935      ENDDO
936      ENDDO
937
938C   And update large scale variables
939cIM 060208 manque DO i + remplace DO k=1,kupper(i)
940cIM 060208     DO k = 1,kupper(i)
941      DO k= 1,klev
942      DO i = 1,klon
943       IF(k .LE. kupper(i)) THEN
944        te(i,k) = te0(i,k) + dtls(i,k)
945        qe(i,k) = qe0(i,k) + dqls(i,k)
946        the(i,k) = te(i,k)/ppi(i,k)
947       ENDIF
948      ENDDO
949      ENDDO
950c
951C
952c     Determine Ptop from buoyancy integral
953c     ---------------------------------------
954c
955c-     1/ Pressure of the level where dth changes sign.
956c
957      DO i=1,klon
958      Ptop_provis(i)=ph(i,1)
959      ENDDO
960c
961      DO k= 2,klev
962      DO i=1,klon
963        IF (Ptop_provis(i) .EQ. ph(i,1) .AND.
964     $      dth(i,k) .GT. -delta_t_min .and.
965     $      dth(i,k-1).LT. -delta_t_min) THEN
966          Ptop_provis(i) = ((dth(i,k)+delta_t_min)*p(i,k-1)
967     $          - (dth(i,k-1)+delta_t_min)*p(i,k)) /(dth(i,k)
968     $          - dth(i,k-1))
969        ENDIF
970      ENDDO
971      ENDDO
972c
973c-     2/ dth integral
974c
975      DO i=1,klon
976       sum_dth(i) = 0.
977       dthmin(i) = -delta_t_min
978       z(i) = 0.
979      ENDDO
980
981      DO k = 1,klev
982      DO i=1,klon
983        dz(i) = -(amax1(ph(i,k+1),Ptop_provis(i))-Ph(i,k))/(rho(i,k)*rg)
984        IF (dz(i) .gt. 0) THEN
985         z(i) = z(i)+dz(i)
986         sum_dth(i) = sum_dth(i) + dth(i,k)*dz(i)
987         dthmin(i) = amin1(dthmin(i),dth(i,k))
988        ENDIF
989      ENDDO
990      ENDDO
991c
992c-     3/ height of triangle with area= sum_dth and base = dthmin
993
994      DO i=1,klon
995       hw(i) = 2.*sum_dth(i)/amin1(dthmin(i),-0.5)
996       hw(i) = amax1(hwmin,hw(i))
997      ENDDO
998c
999c-     4/ now, get Ptop
1000c
1001      DO i=1,klon
1002       ktop(i) = 0
1003       z(i)=0.
1004      ENDDO
1005c
1006      DO k = 1,klev
1007      DO i=1,klon
1008        dz(i) = amin1(-(ph(i,k+1)-Ph(i,k))/(rho(i,k)*rg),hw(i)-z(i))
1009        IF (dz(i) .gt. 0) THEN
1010         z(i) = z(i)+dz(i)
1011         Ptop(i) = Ph(i,k)-rho(i,k)*rg*dz(i)
1012         ktop(i) = k
1013        ENDIF
1014      ENDDO
1015      ENDDO
1016c
1017c      4.5/Correct ktop and ptop
1018c
1019      DO i=1,klon
1020        Ptop_new(i)=ptop(i)
1021      ENDDO
1022c
1023      DO k= klev,2,-1
1024      DO i=1,klon
1025cIM v3JYG; IF (k .GE. ktop(i)
1026       IF (k .LE. ktop(i) .AND.
1027     $      ptop_new(i) .EQ. ptop(i) .AND.
1028     $      dth(i,k) .GT. -delta_t_min .and.
1029     $      dth(i,k-1).LT. -delta_t_min) THEN
1030          Ptop_new(i) = ((dth(i,k)+delta_t_min)*p(i,k-1)
1031     $          - (dth(i,k-1)+delta_t_min)*p(i,k)) /(dth(i,k)
1032     $          - dth(i,k-1))
1033        ENDIF
1034      ENDDO
1035      ENDDO
1036c
1037c
1038      DO i=1,klon
1039      ptop(i) = ptop_new(i)
1040      ENDDO
1041
1042      DO k=klev,1,-1
1043      DO i=1,klon
1044        IF (ph(i,k+1) .LT. ptop(i)) ktop(i)=k
1045      ENDDO
1046      ENDDO
1047c
1048c      5/ Set deltatw & deltaqw to 0 above kupper
1049c
1050      DO k = 1,klev
1051      DO i=1,klon
1052        IF (k .GE. kupper(i)) THEN
1053         deltatw(i,k) = 0.
1054         deltaqw(i,k) = 0.
1055        ENDIF
1056      ENDDO
1057      ENDDO
1058c
1059C
1060       ENDDO      ! end sub-timestep loop
1061C
1062C -----------------------------------------------------------------
1063c   Get back to tendencies per second
1064c
1065      DO k = 1,klev
1066      DO i=1,klon
1067        IF (k .LE. kupper(i)-1) THEN
1068         dtls(i,k) = dtls(i,k)/dtime
1069         dqls(i,k) = dqls(i,k)/dtime
1070         d_deltatw2(i,k)=d_deltatw2(i,k)/dtime
1071         d_deltaqw2(i,k)=d_deltaqw2(i,k)/dtime
1072         d_deltat_gw(i,k) = d_deltat_gw(i,k)/dtime
1073        ENDIF
1074      ENDDO
1075      ENDDO
1076c
1077c
1078c----------------------------------------------------------
1079c   Determine wake final state; recompute wape, cstar, ktop;
1080c   filter out bad wakes.
1081c----------------------------------------------------------
1082c
1083C 2.1 - Undisturbed area and Wake integrals
1084C ---------------------------------------------------------
1085
1086      DO i=1,klon
1087        z(i) = 0.
1088        sum_thu(i) = 0.
1089        sum_tu(i) = 0.
1090        sum_qu(i) = 0.
1091        sum_thvu(i) = 0.
1092        sum_dth(i) = 0.
1093        sum_dq(i) = 0.
1094        sum_rho(i) = 0.
1095        sum_dtdwn(i) = 0.
1096        sum_dqdwn(i) = 0.
1097
1098        av_thu(i) = 0.
1099        av_tu(i) =0.
1100        av_qu(i) =0.
1101        av_thvu(i) = 0.
1102        av_dth(i) = 0.
1103        av_dq(i) = 0.
1104        av_rho(i) =0.
1105        av_dtdwn(i) =0.
1106        av_dqdwn(i) = 0.
1107      ENDDO
1108C Potential temperatures and humidity
1109c----------------------------------------------------------
1110
1111      DO k =1,klev
1112      DO i=1,klon
1113        rho(i,k) = p(i,k)/(rd*te(i,k))
1114        IF(k .eq. 1) THEN
1115          rhoh(i,k) = ph(i,k)/(rd*te(i,k))
1116          zhh(i,k)=0
1117        ELSE
1118          rhoh(i,k) = ph(i,k)*2./(rd*(te(i,k)+te(i,k-1)))
1119          zhh(i,k)=(ph(i,k)-ph(i,k-1))/(-rhoh(i,k)*RG)+zhh(i,k-1)
1120        ENDIF
1121        the(i,k) = te(i,k)/ppi(i,k)
1122        thu(i,k) = (te(i,k) - deltatw(i,k)*sigmaw(i))/ppi(i,k)
1123        tu(i,k) = te(i,k) - deltatw(i,k)*sigmaw(i)
1124        qu(i,k)  =  qe(i,k) - deltaqw(i,k)*sigmaw(i)
1125        rhow(i,k) = p(i,k)/(rd*(te(i,k)+deltatw(i,k)))
1126        dth(i,k) = deltatw(i,k)/ppi(i,k)
1127      ENDDO
1128      ENDDO
1129
1130C Integrals (and wake top level number)
1131C -----------------------------------------------------------
1132
1133C Initialize sum_thvu to 1st level virt. pot. temp.
1134
1135      DO i=1,klon
1136        z(i) = 1.
1137        dz(i) = 1.
1138        sum_thvu(i) =  thu(i,1)*(1.+eps*qu(i,1))*dz(i)
1139        sum_dth(i) = 0.
1140      ENDDO
1141
1142      DO k = 1,klev
1143      DO i=1,klon
1144        dz(i) = -(amax1(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*rg)
1145        IF (dz(i) .GT. 0) THEN
1146         z(i) = z(i)+dz(i)
1147         sum_thu(i) = sum_thu(i) + thu(i,k)*dz(i)
1148         sum_tu(i) = sum_tu(i) + tu(i,k)*dz(i)
1149         sum_qu(i) = sum_qu(i) + qu(i,k)*dz(i)
1150         sum_thvu(i) = sum_thvu(i) + thu(i,k)*(1.+eps*qu(i,k))*dz(i)
1151         sum_dth(i) = sum_dth(i) + dth(i,k)*dz(i)
1152         sum_dq(i) = sum_dq(i) + deltaqw(i,k)*dz(i)
1153         sum_rho(i) = sum_rho(i) + rhow(i,k)*dz(i)
1154         sum_dtdwn(i) = sum_dtdwn(i) + dtdwn(i,k)*dz(i)
1155         sum_dqdwn(i) = sum_dqdwn(i) + dqdwn(i,k)*dz(i)
1156        ENDIF
1157      ENDDO
1158      ENDDO
1159c
1160      DO i=1,klon
1161        hw0(i) = z(i)
1162      ENDDO
1163c
1164C 2.1 - WAPE and mean forcing computation
1165C-------------------------------------------------------------
1166
1167C Means
1168
1169      DO i=1, klon
1170        av_thu(i) = sum_thu(i)/hw0(i)
1171        av_tu(i) = sum_tu(i)/hw0(i)
1172        av_qu(i) = sum_qu(i)/hw0(i)
1173        av_thvu(i) = sum_thvu(i)/hw0(i)
1174        av_dth(i) = sum_dth(i)/hw0(i)
1175        av_dq(i) = sum_dq(i)/hw0(i)
1176        av_rho(i) = sum_rho(i)/hw0(i)
1177        av_dtdwn(i) = sum_dtdwn(i)/hw0(i)
1178        av_dqdwn(i) = sum_dqdwn(i)/hw0(i)
1179
1180        wape2(i) = - rg*hw0(i)*(av_dth(i)
1181     $     + eps*(av_thu(i)*av_dq(i)+av_dth(i)*av_qu(i)+
1182     $     av_dth(i)*av_dq(i) ))/av_thvu(i)
1183      ENDDO
1184
1185C 2.2 Prognostic variable update
1186C ------------------------------------------------------------
1187
1188C Filter out bad wakes
1189c
1190      DO k = 1,klev
1191      DO i=1,klon
1192        IF ( wape2(i) .LT. 0.) THEN
1193          deltatw(i,k) = 0.
1194          deltaqw(i,k) = 0.
1195          dth(i,k) = 0.
1196        ENDIF
1197      ENDDO
1198      ENDDO
1199c
1200
1201      DO i=1, klon
1202      IF ( wape2(i) .LT. 0.) THEN
1203        wape2(i) = 0.
1204        Cstar2(i) = 0.
1205        hw(i) = hwmin
1206        sigmaw(i) = amax1(sigmad,sigd_con(i))
1207        fip(i) = 0.
1208        gwake(i) = .FALSE.
1209      ELSE
1210        if(prt_level.ge.10) print*,'wape2>0'
1211        Cstar2(i) = stark*sqrt(2.*wape2(i))
1212        gwake(i) = .TRUE.
1213      ENDIF
1214      ENDDO
1215c
1216      DO i=1, klon
1217        ktopw(i) = ktop(i)
1218      ENDDO
1219c
1220      DO i=1, klon
1221      IF (ktopw(i) .gt. 0 .and. gwake(i)) then
1222
1223Cjyg1     Utilisation d'un h_efficace constant ( ~ feeding layer)
1224ccc       heff = 600.
1225C      Utilisation de la hauteur hw
1226cc       heff = 0.7*hw
1227       heff(i) = hw(i)
1228
1229       FIP(i) = 0.5*rho(i,ktopw(i))*Cstar2(i)**3*heff(i)*2*
1230     $      sqrt(sigmaw(i)*wdens*3.14)
1231       FIP(i) = alpk * FIP(i)
1232Cjyg2
1233       ELSE
1234         FIP(i) = 0.
1235       ENDIF
1236      ENDDO
1237c
1238C   Limitation de sigmaw
1239c
1240C   sécurité : si le wake occuppe plus de 90 % de la surface de la maille,
1241C              alors il disparait en se mélangeant à la partie undisturbed
1242c
1243      DO k = 1,klev
1244       DO i=1, klon
1245         IF ((sigmaw(i).GT.0.9).or.
1246     $     ((wape(i).ge.wape2(i)).and.(wape2(i).le.1.0))) THEN
1247ccc      IF (sigmaw(i).GT.0.9) THEN
1248          dtls(i,k) = 0.
1249          dqls(i,k) = 0.
1250          deltatw(i,k) = 0.
1251          deltaqw(i,k) = 0.
1252        ENDIF
1253       ENDDO
1254      ENDDO
1255c
1256      DO i=1, klon
1257         IF ((sigmaw(i).GT.0.9).or.
1258     $     ((wape(i).ge.wape2(i)).and.(wape2(i).le.1.0))) THEN
1259ccc      IF (sigmaw(i).GT.0.9) THEN
1260         wape(i) = 0.
1261         hw(i) = hwmin
1262         sigmaw(i) = sigmad
1263         fip(i) = 0.
1264        ELSE
1265         wape(i) = wape2(i)
1266        ENDIF
1267      ENDDO
1268c
1269c
1270      RETURN
1271      END
1272      Subroutine WAKE_scal (p,ph,ppi,dtime,sigd_con
1273     :                ,te0,qe0,omgb
1274     :                ,dtdwn,dqdwn,amdwn,amup,dta,dqa
1275     :                ,wdtPBL,wdqPBL,udtPBL,udqPBL
1276     o                ,deltatw,deltaqw,dth,hw,sigmaw,wape,fip,gfl
1277     o                ,dtls,dqls
1278     o                ,ktopw,omgbdth,dp_omgb,wdens
1279     o                ,tu,qu
1280     o                ,dtKE,dqKE
1281     o                ,dtPBL,dqPBL
1282     o                ,omg,dp_deltomg,spread
1283     o                ,Cstar,d_deltat_gw
1284     o                ,d_deltatw2,d_deltaqw2)
1285
1286***************************************************************
1287*                                                             *
1288* WAKE                                                        *
1289*      retour a un Pupper fixe                                *
1290*                                                             *
1291* written by   :  GRANDPEIX Jean-Yves   09/03/2000            *
1292* modified by :   ROEHRIG Romain        01/29/2007            *
1293***************************************************************
1294c
1295      USE dimphy
1296      IMPLICIT none
1297c============================================================================
1298C
1299C
1300C   But : Decrire le comportement des poches froides apparaissant dans les
1301C        grands systemes convectifs, et fournir l'energie disponible pour
1302C        le declenchement de nouvelles colonnes convectives.
1303C
1304C   Variables d'etat : deltatw    : ecart de temperature wake-undisturbed area
1305C                      deltaqw    : ecart d'humidite wake-undisturbed area
1306C                      sigmaw     : fraction d'aire occupee par la poche.
1307C
1308C   Variable de sortie :
1309c
1310c                        wape : WAke Potential Energy
1311c                        fip  : Front Incident Power (W/m2) - ALP
1312c                        gfl  : Gust Front Length per unit area (m-1)
1313C                        dtls : large scale temperature tendency due to wake
1314C                        dqls : large scale humidity tendency due to wake
1315C                        hw   : hauteur de la poche
1316C                     dp_omgb : vertical gradient of large scale omega
1317C                      omgbdth: flux of Delta_Theta transported by LS omega
1318C                      dtKE   : differential heating (wake - unpertubed)
1319C                      dqKE   : differential moistening (wake - unpertubed)
1320C                      omg    : Delta_omg =vertical velocity diff. wake-undist. (Pa/s)
1321C                 dp_deltomg  : vertical gradient of omg (s-1)
1322C                     spread  : spreading term in dt_wake and dq_wake
1323C                 deltatw     : updated temperature difference (T_w-T_u).
1324C                 deltaqw     : updated humidity difference (q_w-q_u).
1325C                 sigmaw      : updated wake fractional area.
1326C                 d_deltat_gw : delta T tendency due to GW
1327c
1328C   Variables d'entree :
1329c
1330c                        aire : aire de la maille
1331c                        te0  : temperature dans l'environnement  (K)
1332C                        qe0  : humidite dans l'environnement     (kg/kg)
1333C                        omgb : vitesse verticale moyenne sur la maille (Pa/s)
1334C                        dtdwn: source de chaleur due aux descentes (K/s)
1335C                        dqdwn: source d'humidite due aux descentes (kg/kg/s)
1336C                        dta  : source de chaleur due courants satures et detrain  (K/s)
1337C                        dqa  : source d'humidite due aux courants satures et detra (kg/kg/s)
1338C                        amdwn: flux de masse total des descentes, par unite de
1339C                                surface de la maille (kg/m2/s)
1340C                        amup : flux de masse total des ascendances, par unite de
1341C                                surface de la maille (kg/m2/s)
1342C                        p    : pressions aux milieux des couches (Pa)
1343C                        ph   : pressions aux interfaces (Pa)
1344C                        ppi  : (p/p_0)**kapa (adim)
1345C                        dtime: increment temporel (s)
1346c
1347C   Variables internes :
1348c
1349c                        rhow : masse volumique de la poche froide
1350C                        rho  : environment density at P levels
1351C                        rhoh : environment density at Ph levels
1352C                        te   : environment temperature | may change within
1353C                        qe   : environment humidity    | sub-time-stepping
1354C                        the  : environment potential temperature
1355C                        thu  : potential temperature in undisturbed area
1356C                        tu   :  temperature  in undisturbed area
1357C                        qu   : humidity in undisturbed area
1358C                      dp_omgb: vertical gradient og LS omega
1359C                      omgbw  : wake average vertical omega
1360C                     dp_omgbw: vertical gradient of omgbw
1361C                      omgbdq : flux of Delta_q transported by LS omega
1362C                        dth  : potential temperature diff. wake-undist.
1363C                        th1  : first pot. temp. for vertical advection (=thu)
1364C                        th2  : second pot. temp. for vertical advection (=thw)
1365C                        q1   : first humidity for vertical advection
1366C                        q2   : second humidity for vertical advection
1367C                     d_deltatw   : terme de redistribution pour deltatw
1368C                     d_deltaqw   : terme de redistribution pour deltaqw
1369C                      deltatw0   : deltatw initial
1370C                      deltaqw0   : deltaqw initial
1371C                      hw0    : hw initial
1372C                      sigmaw0: sigmaw initial
1373C                      amflux : horizontal mass flux through wake boundary
1374C                      wdens  : number of wakes per unit area (3D) or per
1375C                               unit length (2D)
1376C                      Tgw    : 1 sur la période de onde de gravité
1377c                      Cgw    : vitesse de propagation de onde de gravité
1378c                      LL     : distance entre 2 poches
1379
1380c-------------------------------------------------------------------------
1381c          Déclaration de variables
1382c-------------------------------------------------------------------------
1383
1384#include "dimensions.h"
1385cccc#include "dimphy.h"
1386#include "YOMCST.h"
1387#include "cvthermo.h"
1388#include "iniprint.h"
1389
1390c Arguments en entree
1391c--------------------
1392
1393      REAL p(klev),ph(klev+1),ppi(klev)
1394      REAL dtime
1395      REAL te0(klev),qe0(klev)
1396      REAL omgb(klev+1)
1397      REAL dtdwn(klev), dqdwn(klev)
1398      REAL wdtPBL(klev),wdqPBL(klev)
1399      REAL udtPBL(klev),udqPBL(klev)
1400      REAL amdwn(klev), amup(klev)
1401      REAL dta(klev), dqa(klev)
1402      REAL sigd_con
1403
1404c Sorties
1405c--------
1406
1407      REAL deltatw(klev), deltaqw(klev), dth(klev)
1408      REAL tu(klev), qu(klev)
1409      REAL dtls(klev), dqls(klev)
1410      REAL dtKE(klev), dqKE(klev)
1411      REAL dtPBL(klev), dqPBL(klev)
1412      REAL spread(klev)
1413      REAL d_deltatgw(klev)
1414      REAL d_deltatw2(klev), d_deltaqw2(klev)
1415      REAL omgbdth(klev+1), omg(klev+1)
1416      REAL dp_omgb(klev), dp_deltomg(klev)
1417      REAL d_deltat_gw(klev)
1418      REAL hw, sigmaw, wape, fip, gfl, Cstar
1419      INTEGER ktopw
1420
1421c Variables internes
1422c-------------------
1423
1424c Variables à fixer
1425      REAL ALON
1426      REAL coefgw
1427      REAL wdens0, wdens
1428      REAL stark
1429      REAL alpk
1430      REAL delta_t_min
1431      REAL Pupper
1432      INTEGER nsub
1433      REAL dtimesub
1434      REAL sigmad, hwmin
1435
1436c Variables de sauvegarde
1437      REAL deltatw0(klev)
1438      REAL deltaqw0(klev)
1439      REAL te(klev), qe(klev)
1440      REAL sigmaw0, sigmaw1
1441
1442c Variables pour les GW
1443      REAL LL
1444      REAL N2(klev)
1445      REAL Cgw(klev)
1446      REAL Tgw(klev)
1447
1448c Variables liées au calcul de hw
1449      REAL ptop_provis, ptop, ptop_new
1450      REAL sum_dth
1451      REAL dthmin
1452      REAL z, dz, hw0
1453      INTEGER ktop, kupper
1454
1455c Autres variables internes
1456      INTEGER isubstep, k
1457
1458      REAL sum_thu, sum_tu, sum_qu,sum_thvu
1459      REAL sum_dq, sum_rho
1460      REAL sum_dtdwn, sum_dqdwn
1461      REAL av_thu, av_tu, av_qu, av_thvu
1462      REAL av_dth, av_dq, av_rho
1463      REAL av_dtdwn, av_dqdwn
1464
1465      REAL rho(klev), rhoh(klev+1), rhow(klev)
1466      REAL rhow_moyen(klev)
1467      REAL zh(klev), zhh(klev+1)
1468      REAL epaisseur1(klev), epaisseur2(klev)
1469
1470      REAL the(klev), thu(klev)
1471
1472      REAL d_deltatw(klev), d_deltaqw(klev)
1473
1474      REAL omgbw(klev+1), omgtop
1475      REAL dp_omgbw(klev)
1476      REAL ztop, dztop
1477      REAL alpha_up(klev)
1478     
1479      REAL RRe1, RRe2, RRd1, RRd2
1480      REAL Th1(klev), Th2(klev), q1(klev), q2(klev)
1481      REAL D_Th1(klev), D_Th2(klev), D_dth(klev)
1482      REAL D_q1(klev), D_q2(klev), D_dq(klev)
1483      REAL omgbdq(klev)
1484
1485      REAL ff, gg
1486      REAL wape2, Cstar2, heff
1487
1488      REAL Crep(klev)
1489      REAL Crep_upper, Crep_sol
1490
1491C-------------------------------------------------------------------------
1492c         Initialisations
1493c-------------------------------------------------------------------------
1494
1495c      print*, 'wake initialisations'
1496
1497c   Essais d'initialisation avec sigmaw = 0.02 et hw = 10.
1498c-------------------------------------------------------------------------
1499
1500      DATA sigmad, hwmin /.02,10./
1501
1502C Longueur de maille (en m)
1503c-------------------------------------------------------------------------
1504
1505c      ALON = 3.e5
1506      ALON = 1.e6
1507
1508
1509C Configuration de coefgw,stark,wdens (22/02/06 by YU Jingmei)
1510c
1511c      coefgw : Coefficient pour les ondes de gravité
1512c       stark : Coefficient k dans Cstar=k*sqrt(2*WAPE)
1513c       wdens : Densité de poche froide par maille
1514c-------------------------------------------------------------------------
1515
1516      coefgw=10
1517c      coefgw=1
1518c      wdens0 = 1.0/(alon**2)   
1519      wdens = 1.0/(alon**2)       
1520      stark = 0.50
1521cCRtest
1522      alpk=0.1
1523c      alpk = 1.0
1524c      alpk = 0.5
1525c      alpk = 0.05
1526      Crep_upper=0.9
1527      Crep_sol=1.0
1528
1529
1530C Minimum value for |T_wake - T_undist|. Used for wake top definition
1531c-------------------------------------------------------------------------
1532
1533      delta_t_min = 0.2
1534
1535
1536C 1. - Save initial values and initialize tendencies
1537C --------------------------------------------------
1538
1539      DO k=1,klev
1540        deltatw0(k) = deltatw(k)
1541        deltaqw0(k)= deltaqw(k)
1542        te(k) = te0(k)
1543        qe(k) = qe0(k)
1544        dtls(k) = 0.
1545        dqls(k) = 0.
1546        d_deltat_gw(k)=0.
1547        d_deltatw2(k)=0.
1548        d_deltaqw2(k)=0.
1549      ENDDO
1550c      sigmaw1=sigmaw
1551c      IF (sigd_con.GT.sigmaw1) THEN
1552c      print*, 'sigmaw,sigd_con', sigmaw, sigd_con
1553c      ENDIF
1554      sigmaw = max(sigmaw,sigd_con)
1555      sigmaw = max(sigmaw,sigmad)
1556      sigmaw = min(sigmaw,0.99)
1557      sigmaw0 = sigmaw
1558c      wdens=wdens0/(10.*sigmaw)
1559c      IF (sigd_con.GT.sigmaw1) THEN
1560c      print*, 'sigmaw1,sigd1', sigmaw, sigd_con
1561c      ENDIF
1562
1563C 2. - Prognostic part
1564C =========================================================
1565
1566c      print *, 'prognostic wake computation'
1567
1568
1569C 2.1 - Undisturbed area and Wake integrals
1570C ---------------------------------------------------------
1571
1572      z = 0.
1573      ktop=0
1574      kupper = 0
1575      sum_thu = 0.
1576      sum_tu = 0.
1577      sum_qu = 0.
1578      sum_thvu = 0.
1579      sum_dth = 0.
1580      sum_dq = 0.
1581      sum_rho = 0.
1582      sum_dtdwn = 0.
1583      sum_dqdwn = 0.
1584
1585      av_thu = 0.
1586      av_tu =0.
1587      av_qu =0.
1588      av_thvu = 0.
1589      av_dth = 0.
1590      av_dq = 0.
1591      av_rho =0.
1592      av_dtdwn =0.
1593      av_dqdwn = 0.
1594
1595C Potential temperatures and humidity
1596c----------------------------------------------------------
1597
1598      DO k =1,klev
1599        rho(k) = p(k)/(rd*te(k))
1600        IF(k .eq. 1) THEN
1601          rhoh(k) = ph(k)/(rd*te(k))
1602          zhh(k)=0
1603        ELSE
1604          rhoh(k) = ph(k)*2./(rd*(te(k)+te(k-1)))
1605          zhh(k)=(ph(k)-ph(k-1))/(-rhoh(k)*RG)+zhh(k-1)
1606        ENDIF
1607        the(k) = te(k)/ppi(k)
1608        thu(k) = (te(k) - deltatw(k)*sigmaw)/ppi(k)
1609        tu(k) = te(k) - deltatw(k)*sigmaw
1610        qu(k)  =  qe(k) - deltaqw(k)*sigmaw
1611        rhow(k) = p(k)/(rd*(te(k)+deltatw(k)))
1612        dth(k) = deltatw(k)/ppi(k)
1613        LL = (1-sqrt(sigmaw))/sqrt(wdens)       
1614      ENDDO
1615       
1616      DO k = 1, klev-1
1617        IF(k.eq.1) THEN
1618          N2(k)=0
1619        ELSE
1620          N2(k)=max(0.,-RG**2/the(k)*rho(k)*(the(k+1)-the(k-1))
1621     $           /(p(k+1)-p(k-1)))
1622        ENDIF
1623        ZH(k)=(zhh(k)+zhh(k+1))/2
1624
1625        Cgw(k)=sqrt(N2(k))*ZH(k)
1626        Tgw(k)=coefgw*Cgw(k)/LL
1627      ENDDO
1628         
1629      N2(klev)=0
1630      ZH(klev)=0
1631      Cgw(klev)=0
1632      Tgw(klev)=0
1633
1634c  Calcul de la masse volumique moyenne de la colonne
1635c-----------------------------------------------------------------
1636
1637      DO k=1,klev
1638        epaisseur1(k)=0.
1639        epaisseur2(k)=0.
1640      ENDDO
1641
1642      epaisseur1(1)= -(Ph(2)-Ph(1))/(rho(1)*rg)+1.
1643      epaisseur2(1)= -(Ph(2)-Ph(1))/(rho(1)*rg)+1.
1644      rhow_moyen(1) = rhow(1)
1645
1646      DO k = 2, klev
1647        epaisseur1(k)= -(Ph(k+1)-Ph(k))/(rho(k)*rg) +1.
1648        epaisseur2(k)=epaisseur2(k-1)+epaisseur1(k)
1649        rhow_moyen(k) = (rhow_moyen(k-1)*epaisseur2(k-1)+
1650     $                 rhow(k)*epaisseur1(k))/epaisseur2(k)
1651      ENDDO
1652
1653
1654C Choose an integration bound well above wake top
1655c-----------------------------------------------------------------
1656
1657c       Pupper = 50000.  ! melting level
1658       Pupper = 60000.
1659c       Pupper = 70000.
1660
1661
1662C    Determine Wake top pressure (Ptop) from buoyancy integral
1663C-----------------------------------------------------------------
1664
1665c-1/ Pressure of the level where dth becomes less than delta_t_min.
1666
1667      Ptop_provis=ph(1)
1668      DO k= 2,klev
1669        IF (dth(k) .GT. -delta_t_min .and.
1670     $      dth(k-1).LT. -delta_t_min) THEN
1671          Ptop_provis = ((dth(k)+delta_t_min)*p(k-1)
1672     $          - (dth(k-1)+delta_t_min)*p(k)) /(dth(k) - dth(k-1))
1673          GO TO 25
1674        ENDIF
1675      ENDDO
167625    CONTINUE
1677
1678c-2/ dth integral
1679
1680      sum_dth = 0.
1681      dthmin = -delta_t_min
1682      z = 0.
1683
1684      DO k = 1,klev
1685        dz = -(max(ph(k+1),Ptop_provis)-Ph(k))/(rho(k)*rg)
1686        IF (dz .le. 0) GO TO 40
1687        z = z+dz
1688        sum_dth = sum_dth + dth(k)*dz
1689        dthmin = min(dthmin,dth(k))
1690      ENDDO
169140    CONTINUE
1692
1693c-3/ height of triangle with area= sum_dth and base = dthmin
1694
1695      hw0 = 2.*sum_dth/min(dthmin,-0.5)
1696      hw0 = max(hwmin,hw0)
1697
1698c-4/ now, get Ptop
1699
1700      z = 0.
1701      ptop = ph(1)
1702
1703      DO k = 1,klev
1704        dz = min(-(ph(k+1)-Ph(k))/(rho(k)*rg),hw0-z)
1705        IF (dz .le. 0) GO TO 45
1706        z = z+dz
1707        Ptop = Ph(k)-rho(k)*rg*dz
1708      ENDDO
170945    CONTINUE
1710
1711
1712C-5/ Determination de ktop et kupper
1713
1714      DO k=klev,1,-1
1715        IF (ph(k+1) .lt. ptop) ktop=k
1716        IF (ph(k+1) .lt. pupper) kupper=k
1717      ENDDO
1718
1719c-6/ Correct ktop and ptop
1720
1721      Ptop_new=ptop
1722      DO k= ktop,2,-1
1723        IF (dth(k) .GT. -delta_t_min .and.
1724     $      dth(k-1).LT. -delta_t_min) THEN
1725          Ptop_new = ((dth(k)+delta_t_min)*p(k-1)
1726     $          - (dth(k-1)+delta_t_min)*p(k)) /(dth(k) - dth(k-1))
1727          GO TO 225
1728        ENDIF
1729      ENDDO
1730225   CONTINUE
1731
1732      ptop = ptop_new
1733
1734      DO k=klev,1,-1
1735        IF (ph(k+1) .lt. ptop) ktop=k
1736      ENDDO
1737
1738c Set deltatw & deltaqw to 0 above kupper
1739c-----------------------------------------------------------
1740
1741      DO k = kupper,klev
1742        deltatw(k) = 0.
1743        deltaqw(k) = 0.
1744      ENDDO
1745
1746
1747C Vertical gradient of LS omega
1748C------------------------------------------------------------
1749
1750      DO k = 1,kupper
1751        dp_omgb(k) = (omgb(k+1) - omgb(k))/(ph(k+1)-ph(k))
1752      ENDDO
1753
1754
1755C Integrals (and wake top level number)
1756C -----------------------------------------------------------
1757
1758C Initialize sum_thvu to 1st level virt. pot. temp.
1759
1760      z = 1.
1761      dz = 1.
1762      sum_thvu =  thu(1)*(1.+eps*qu(1))*dz
1763      sum_dth = 0.
1764
1765      DO k = 1,klev
1766        dz = -(max(ph(k+1),Ptop)-Ph(k))/(rho(k)*rg)
1767        IF (dz .LE. 0) GO TO 50
1768        z = z+dz
1769        sum_thu = sum_thu + thu(k)*dz
1770        sum_tu = sum_tu + tu(k)*dz
1771        sum_qu = sum_qu + qu(k)*dz
1772        sum_thvu = sum_thvu + thu(k)*(1.+eps*qu(k))*dz
1773        sum_dth = sum_dth + dth(k)*dz
1774        sum_dq = sum_dq + deltaqw(k)*dz
1775        sum_rho = sum_rho + rhow(k)*dz
1776        sum_dtdwn = sum_dtdwn + dtdwn(k)*dz
1777        sum_dqdwn = sum_dqdwn + dqdwn(k)*dz
1778      ENDDO
177950    CONTINUE
1780
1781      hw0 = z
1782
1783C 2.1 - WAPE and mean forcing computation
1784C-------------------------------------------------------------
1785
1786C Means
1787
1788      av_thu = sum_thu/hw0
1789      av_tu = sum_tu/hw0
1790      av_qu = sum_qu/hw0
1791      av_thvu = sum_thvu/hw0
1792c      av_thve = sum_thve/hw0
1793      av_dth = sum_dth/hw0
1794      av_dq = sum_dq/hw0
1795      av_rho = sum_rho/hw0
1796      av_dtdwn = sum_dtdwn/hw0
1797      av_dqdwn = sum_dqdwn/hw0
1798
1799      wape = - rg*hw0*(av_dth
1800     $     + eps*(av_thu*av_dq+av_dth*av_qu+av_dth*av_dq ))/av_thvu
1801
1802C 2.2 Prognostic variable update
1803C ------------------------------------------------------------
1804
1805C Filter out bad wakes
1806
1807      IF ( wape .LT. 0.) THEN
1808        if(prt_level.ge.10) print*,'wape<0'
1809        wape = 0.
1810        hw = hwmin
1811        sigmaw = max(sigmad,sigd_con)
1812        fip = 0.
1813        DO k = 1,klev
1814          deltatw(k) = 0.
1815          deltaqw(k) = 0.
1816          dth(k) = 0.
1817        ENDDO
1818      ELSE
1819        if(prt_level.ge.10) print*,'wape>0'
1820        Cstar = stark*sqrt(2.*wape)
1821      ENDIF
1822
1823C------------------------------------------------------------------
1824C    Sub-time-stepping
1825C------------------------------------------------------------------
1826
1827c      nsub=36
1828      nsub=10
1829      dtimesub=dtime/nsub
1830
1831c------------------------------------------------------------
1832      DO isubstep = 1,nsub
1833c------------------------------------------------------------
1834
1835c        print*,'---------------','substep=',isubstep,'-------------'
1836
1837c  Evolution of sigmaw
1838
1839
1840        gfl = 2.*sqrt(3.14*wdens*sigmaw)           
1841
1842        sigmaw =sigmaw + gfl*Cstar*dtimesub
1843        sigmaw =min(sigmaw,0.99)     !!!!!!!!
1844c        wdens = wdens0/(10.*sigmaw)
1845c        sigmaw =max(sigmaw,sigd_con)
1846c        sigmaw =max(sigmaw,sigmad)
1847
1848c calcul de la difference de vitesse verticale poche - zone non perturbee
1849
1850        z= 0.
1851        dp_deltomg(1:klev)=0.
1852        omg(1:klev+1)=0.
1853
1854        omg(1) = 0.
1855        dp_deltomg(1) = -(gfl*Cstar)/(sigmaw * (1-sigmaw))
1856
1857        DO k=2,ktop
1858          dz = -(Ph(k)-Ph(k-1))/(rho(k-1)*rg)
1859          z = z+dz
1860          dp_deltomg(k)= dp_deltomg(1)
1861          omg(k)= dp_deltomg(1)*z
1862        ENDDO
1863
1864        dztop=-(Ptop-Ph(ktop))/(rho(ktop)*rg)
1865        ztop = z+dztop
1866        omgtop=dp_deltomg(1)*ztop
1867
1868
1869c Conversion de la vitesse verticale de m/s a Pa/s
1870
1871        omgtop = -rho(ktop)*rg*omgtop
1872        dp_deltomg(1) = omgtop/(ptop-ph(1))
1873
1874        DO k = 1,ktop
1875          omg(k) = - rho(k)*rg*omg(k)
1876          dp_deltomg(k) = dp_deltomg(1)
1877        ENDDO
1878
1879c   raccordement lineaire de omg de ptop a pupper
1880
1881      IF (kupper .GT. ktop) THEN
1882        omg(kupper+1) = - Rg*amdwn(kupper+1)/sigmaw
1883     $                + Rg*amup(kupper+1)/(1.-sigmaw)
1884        dp_deltomg(kupper) = (omgtop-omg(kupper+1))/(Ptop-Pupper)
1885        DO k=ktop+1,kupper
1886          dp_deltomg(k) = dp_deltomg(kupper)
1887          omg(k) = omgtop+(ph(k)-Ptop)*dp_deltomg(kupper)
1888        ENDDO
1889      ENDIF
1890
1891c   Compute wake average vertical velocity omgbw
1892
1893      DO k = 1,klev+1
1894        omgbw(k) = omgb(k)+(1.-sigmaw)*omg(k)
1895      ENDDO
1896
1897c  and its vertical gradient dp_omgbw
1898
1899      DO k = 1,klev
1900        dp_omgbw(k) = (omgbw(k+1)-omgbw(k))/(ph(k+1)-ph(k))
1901      ENDDO
1902
1903
1904c  Upstream coefficients for omgb velocity
1905c--    (alpha_up(k) is the coefficient of the value at level k)
1906c--    (1-alpha_up(k) is the coefficient of the value at level k-1)
1907
1908      DO k = 1,klev
1909       alpha_up(k) = 0.
1910       IF (omgb(k) .GT. 0.) alpha_up(k) = 1.
1911      ENDDO
1912
1913c  Matrix expressing [The,deltatw] from  [Th1,Th2]
1914
1915      RRe1 = 1.-sigmaw
1916      RRe2 = sigmaw
1917      RRd1 = -1.
1918      RRd2 = 1.
1919
1920c Get [Th1,Th2], dth and [q1,q2]
1921
1922      DO k = 1,kupper+1
1923        dth(k) = deltatw(k)/ppi(k)
1924        Th1(k) = the(k) - sigmaw     *dth(k)   ! undisturbed area
1925        Th2(k) = the(k) + (1.-sigmaw)*dth(k)   ! wake
1926        q1(k) = qe(k) - sigmaw     *deltaqw(k) ! undisturbed area
1927        q2(k) = qe(k) + (1.-sigmaw)*deltaqw(k) ! wake
1928      ENDDO
1929
1930      D_Th1(1) = 0.
1931      D_Th2(1) = 0.
1932      D_dth(1) = 0.
1933      D_q1(1) = 0.
1934      D_q2(1) = 0.
1935      D_dq(1) = 0.
1936
1937      DO k = 2,kupper+1 !   loop on interfaces
1938        D_Th1(k) = Th1(k-1)-Th1(k)
1939        D_Th2(k) = Th2(k-1)-Th2(k)
1940        D_dth(k) = dth(k-1)-dth(k)
1941        D_q1(k) = q1(k-1)-q1(k)
1942        D_q2(k) = q2(k-1)-q2(k)
1943        D_dq(k) = deltaqw(k-1)-deltaqw(k)
1944      ENDDO
1945
1946      omgbdth(1) = 0.
1947      omgbdq(1) = 0.
1948
1949      DO k = 2,kupper+1  !   loop on interfaces
1950        omgbdth(k) = omgb(k)*(    dth(k-1) -     dth(k))
1951        omgbdq(k)  = omgb(k)*(deltaqw(k-1) - deltaqw(k))
1952      ENDDO
1953
1954
1955c-----------------------------------------------------------------
1956      DO k=1,kupper-1
1957c-----------------------------------------------------------------
1958c
1959c   Compute redistribution (advective) term
1960c
1961         d_deltatw(k) =
1962     $             dtimesub/(Ph(k)-Ph(k+1))*(
1963     $       RRd1*omg(k  )*sigmaw     *D_Th1(k)
1964     $      -RRd2*omg(k+1)*(1.-sigmaw)*D_Th2(k+1)
1965     $      -(1.-alpha_up(k))*omgbdth(k) - alpha_up(k+1)*omgbdth(k+1)
1966     $                      )*ppi(k)
1967c         print*,'d_deltatw=',d_deltatw(k)
1968c
1969         d_deltaqw(k) =
1970     $             dtimesub/(Ph(k)-Ph(k+1))*(
1971     $       RRd1*omg(k  )*sigmaw     *D_q1(k)
1972     $      -RRd2*omg(k+1)*(1.-sigmaw)*D_q2(k+1)
1973     $      -(1.-alpha_up(k))*omgbdq(k) - alpha_up(k+1)*omgbdq(k+1)
1974     $                      )
1975c         print*,'d_deltaqw=',d_deltaqw(k)
1976c
1977c   and increment large scale tendencies
1978c
1979         dtls(k) = dtls(k) +
1980     $               dtimesub*(
1981     $        ( RRe1*omg(k  )*sigmaw     *D_Th1(k)
1982     $         -RRe2*omg(k+1)*(1.-sigmaw)*D_Th2(k+1) )
1983     $               /(Ph(k)-Ph(k+1))
1984     $         -sigmaw*(1.-sigmaw)*dth(k)*dp_deltomg(k)
1985     $                      )*ppi(k)
1986c         print*,'dtls=',dtls(k)
1987c
1988         dqls(k) = dqls(k) +
1989     $               dtimesub*(
1990     $        ( RRe1*omg(k  )*sigmaw     *D_q1(k)
1991     $         -RRe2*omg(k+1)*(1.-sigmaw)*D_q2(k+1) )
1992     $               /(Ph(k)-Ph(k+1))
1993     $         -sigmaw*(1.-sigmaw)*deltaqw(k)*dp_deltomg(k)
1994     $                      )
1995c         print*,'dqls=',dqls(k)
1996
1997c-------------------------------------------------------------------
1998      ENDDO
1999c------------------------------------------------------------------
2000
2001C   Increment state variables
2002
2003      DO k = 1,kupper-1
2004
2005c Coefficient de répartition
2006
2007        Crep(k)=Crep_sol*(ph(kupper)-ph(k))/(ph(kupper)-ph(1))
2008        Crep(k)=Crep(k)+Crep_upper*(ph(1)-ph(k))/(p(1)-ph(kupper))
2009       
2010
2011c Reintroduce compensating subsidence term.
2012
2013c        dtKE(k)=(dtdwn(k)*Crep(k))/sigmaw
2014c        dtKE(k)=dtKE(k)-(dtdwn(k)*(1-Crep(k))+dta(k))
2015c     .                   /(1-sigmaw)
2016c        dqKE(k)=(dqdwn(k)*Crep(k))/sigmaw
2017c        dqKE(k)=dqKE(k)-(dqdwn(k)*(1-Crep(k))+dqa(k))
2018c     .                   /(1-sigmaw)
2019c
2020c        dtKE(k)=(dtdwn(k)*Crep(k)+(1-Crep(k))*dta(k))/sigmaw
2021c        dtKE(k)=dtKE(k)-(dtdwn(k)*(1-Crep(k))+dta(k)*Crep(k))
2022c     .                   /(1-sigmaw)
2023c        dqKE(k)=(dqdwn(k)*Crep(k)+(1-Crep(k))*dqa(k))/sigmaw
2024c        dqKE(k)=dqKE(k)-(dqdwn(k)*(1-Crep(k))+dqa(k)*Crep(k))
2025c     .                   /(1-sigmaw)
2026
2027        dtKE(k)=(dtdwn(k)/sigmaw - dta(k)/(1.-sigmaw))
2028        dqKE(k)=(dqdwn(k)/sigmaw - dqa(k)/(1.-sigmaw))
2029c        print*,'dtKE=',dtKE(k)
2030c        print*,'dqKE=',dqKE(k)
2031c
2032        dtPBL(k)=(wdtPBL(k)/sigmaw - udtPBL(k)/(1.-sigmaw))
2033        dqPBL(k)=(wdqPBL(k)/sigmaw - udqPBL(k)/(1.-sigmaw))
2034c
2035        spread(k) = (1.-sigmaw)*dp_deltomg(k)+gfl*Cstar/sigmaw
2036c        print*,'spread=',spread(k)
2037
2038
2039c ajout d'un effet onde de gravité -Tgw(k)*deltatw(k) 03/02/06 YU Jingmei
2040
2041        d_deltat_gw(k)=d_deltat_gw(k)-Tgw(k)*deltatw(k)* dtimesub
2042c        print*,'d_delta_gw=',d_deltat_gw(k)
2043        ff=d_deltatw(k)/dtimesub
2044
2045c Sans GW
2046c
2047c        deltatw(k)=deltatw(k)+dtimesub*(ff+dtKE(k)-spread(k)*deltatw(k))
2048c
2049c GW formule 1
2050c
2051c        deltatw(k) = deltatw(k)+dtimesub*
2052c     $         (ff+dtKE(k) - spread(k)*deltatw(k)-Tgw(k)*deltatw(k))
2053c
2054c GW formule 2
2055
2056        IF (dtimesub*Tgw(k).lt.1.e-10) THEN
2057          deltatw(k) = deltatw(k)+dtimesub*
2058     $          (ff+dtKE(k)+dtPBL(k)
2059     $          - spread(k)*deltatw(k)-Tgw(k)*deltatw(k))
2060        ELSE
2061           deltatw(k) = deltatw(k)+1/Tgw(k)*(1-exp(-dtimesub*Tgw(k)))*
2062     $          (ff+dtKE(k)+dtPBL(k)
2063     $          - spread(k)*deltatw(k)-Tgw(k)*deltatw(k))
2064        ENDIF
2065   
2066        dth(k) = deltatw(k)/ppi(k)
2067
2068        gg=d_deltaqw(k)/dtimesub
2069
2070       deltaqw(k) = deltaqw(k) +
2071     $         dtimesub*(gg+ dqKE(k)+dqPBL(k) - spread(k)*deltaqw(k))
2072
2073       d_deltatw2(k)=d_deltatw2(k)+d_deltatw(k)
2074       d_deltaqw2(k)=d_deltaqw2(k)+d_deltaqw(k)
2075      ENDDO
2076
2077C   And update large scale variables
2078
2079      DO k = 1,kupper
2080        te(k) = te0(k) + dtls(k)
2081        qe(k) = qe0(k) + dqls(k)
2082        the(k) = te(k)/ppi(k)
2083      ENDDO
2084
2085c     Determine Ptop from buoyancy integral
2086c----------------------------------------------------------------------
2087
2088c-1/ Pressure of the level where dth changes sign.
2089
2090      Ptop_provis=ph(1)
2091
2092      DO k= 2,klev
2093        IF (dth(k) .GT. -delta_t_min .and.
2094     $      dth(k-1).LT. -delta_t_min) THEN
2095          Ptop_provis = ((dth(k)+delta_t_min)*p(k-1)
2096     $          - (dth(k-1)+delta_t_min)*p(k)) /(dth(k) - dth(k-1))
2097        GO TO 65
2098        ENDIF
2099      ENDDO
210065    CONTINUE
2101
2102c-2/ dth integral
2103
2104      sum_dth = 0.
2105      dthmin = -delta_t_min
2106      z = 0.
2107
2108      DO k = 1,klev
2109        dz = -(max(ph(k+1),Ptop_provis)-Ph(k))/(rho(k)*rg)
2110        IF (dz .le. 0) GO TO 70
2111        z = z+dz
2112        sum_dth = sum_dth + dth(k)*dz
2113        dthmin = min(dthmin,dth(k))
2114      ENDDO
211570    CONTINUE
2116
2117c-3/ height of triangle with area= sum_dth and base = dthmin
2118
2119      hw = 2.*sum_dth/min(dthmin,-0.5)
2120      hw = max(hwmin,hw)
2121
2122c-4/ now, get Ptop
2123
2124      ktop = 0
2125      z=0.
2126
2127      DO k = 1,klev
2128        dz = min(-(ph(k+1)-Ph(k))/(rho(k)*rg),hw-z)
2129        IF (dz .le. 0) GO TO 75
2130        z = z+dz
2131        Ptop = Ph(k)-rho(k)*rg*dz
2132        ktop = k
2133      ENDDO
213475    CONTINUE
2135
2136c-5/Correct ktop and ptop
2137
2138      Ptop_new=ptop
2139
2140      DO k= ktop,2,-1
2141        IF (dth(k) .GT. -delta_t_min .and.
2142     $      dth(k-1).LT. -delta_t_min) THEN
2143          Ptop_new = ((dth(k)+delta_t_min)*p(k-1)
2144     $          - (dth(k-1)+delta_t_min)*p(k)) /(dth(k) - dth(k-1))
2145          GO TO 275
2146        ENDIF
2147      ENDDO
2148275   CONTINUE
2149
2150      ptop = ptop_new
2151
2152      DO k=klev,1,-1
2153        IF (ph(k+1) .LT. ptop) ktop=k
2154      ENDDO
2155
2156c-6/ Set deltatw & deltaqw to 0 above kupper
2157
2158      DO k = kupper,klev
2159        deltatw(k) = 0.
2160        deltaqw(k) = 0.
2161      ENDDO
2162
2163c------------------------------------------------------------------
2164      ENDDO      ! end sub-timestep loop
2165C -----------------------------------------------------------------
2166
2167c   Get back to tendencies per second
2168
2169      DO k = 1,kupper-1
2170        dtls(k) = dtls(k)/dtime
2171        dqls(k) = dqls(k)/dtime
2172        d_deltatw2(k)=d_deltatw2(k)/dtime
2173        d_deltaqw2(k)=d_deltaqw2(k)/dtime
2174        d_deltat_gw(k) = d_deltat_gw(k)/dtime
2175      ENDDO
2176
2177C 2.1 - Undisturbed area and Wake integrals
2178C ---------------------------------------------------------
2179
2180      z = 0.
2181      sum_thu = 0.
2182      sum_tu = 0.
2183      sum_qu = 0.
2184      sum_thvu = 0.
2185      sum_dth = 0.
2186      sum_dq = 0.
2187      sum_rho = 0.
2188      sum_dtdwn = 0.
2189      sum_dqdwn = 0.
2190
2191      av_thu = 0.
2192      av_tu =0.
2193      av_qu =0.
2194      av_thvu = 0.
2195      av_dth = 0.
2196      av_dq = 0.
2197      av_rho =0.
2198      av_dtdwn =0.
2199      av_dqdwn = 0.
2200
2201C Potential temperatures and humidity
2202c----------------------------------------------------------
2203
2204      DO k =1,klev
2205        rho(k) = p(k)/(rd*te(k))
2206        IF(k .eq. 1) THEN
2207          rhoh(k) = ph(k)/(rd*te(k))
2208          zhh(k)=0
2209        ELSE
2210          rhoh(k) = ph(k)*2./(rd*(te(k)+te(k-1)))
2211          zhh(k)=(ph(k)-ph(k-1))/(-rhoh(k)*RG)+zhh(k-1)
2212        ENDIF
2213        the(k) = te(k)/ppi(k)
2214        thu(k) = (te(k) - deltatw(k)*sigmaw)/ppi(k)
2215        tu(k) = te(k) - deltatw(k)*sigmaw
2216        qu(k)  =  qe(k) - deltaqw(k)*sigmaw
2217        rhow(k) = p(k)/(rd*(te(k)+deltatw(k)))
2218        dth(k) = deltatw(k)/ppi(k)
2219       
2220      ENDDO
2221
2222C Integrals (and wake top level number)
2223C -----------------------------------------------------------
2224
2225C Initialize sum_thvu to 1st level virt. pot. temp.
2226
2227      z = 1.
2228      dz = 1.
2229      sum_thvu =  thu(1)*(1.+eps*qu(1))*dz
2230      sum_dth = 0.
2231
2232      DO k = 1,klev
2233        dz = -(max(ph(k+1),Ptop)-Ph(k))/(rho(k)*rg)
2234
2235        IF (dz .LE. 0) GO TO 51
2236        z = z+dz
2237        sum_thu = sum_thu + thu(k)*dz
2238        sum_tu = sum_tu + tu(k)*dz
2239        sum_qu = sum_qu + qu(k)*dz
2240        sum_thvu = sum_thvu + thu(k)*(1.+eps*qu(k))*dz
2241        sum_dth = sum_dth + dth(k)*dz
2242        sum_dq = sum_dq + deltaqw(k)*dz
2243        sum_rho = sum_rho + rhow(k)*dz
2244        sum_dtdwn = sum_dtdwn + dtdwn(k)*dz
2245        sum_dqdwn = sum_dqdwn + dqdwn(k)*dz
2246      ENDDO
2247 51   CONTINUE
2248
2249      hw0 = z
2250
2251C 2.1 - WAPE and mean forcing computation
2252C-------------------------------------------------------------
2253
2254C Means
2255
2256      av_thu = sum_thu/hw0
2257      av_tu = sum_tu/hw0
2258      av_qu = sum_qu/hw0
2259      av_thvu = sum_thvu/hw0
2260      av_dth = sum_dth/hw0
2261      av_dq = sum_dq/hw0
2262      av_rho = sum_rho/hw0
2263      av_dtdwn = sum_dtdwn/hw0
2264      av_dqdwn = sum_dqdwn/hw0
2265
2266      wape2 = - rg*hw0*(av_dth
2267     $     + eps*(av_thu*av_dq+av_dth*av_qu+av_dth*av_dq ))/av_thvu
2268
2269
2270C 2.2 Prognostic variable update
2271C ------------------------------------------------------------
2272
2273C Filter out bad wakes
2274
2275      IF ( wape2 .LT. 0.) THEN
2276        if(prt_level.ge.10) print*,'wape2<0'
2277        wape2 = 0.
2278        hw = hwmin
2279        sigmaw = max(sigmad,sigd_con)
2280        fip = 0.
2281        DO k = 1,klev
2282          deltatw(k) = 0.
2283          deltaqw(k) = 0.
2284          dth(k) = 0.
2285        ENDDO
2286      ELSE
2287        if(prt_level.ge.10) print*,'wape2>0'
2288        Cstar2 = stark*sqrt(2.*wape2)
2289
2290      ENDIF
2291
2292      ktopw = ktop
2293
2294      IF (ktopw .gt. 0) then
2295
2296Cjyg1     Utilisation d'un h_efficace constant ( ~ feeding layer)
2297ccc       heff = 600.
2298C      Utilisation de la hauteur hw
2299cc       heff = 0.7*hw
2300       heff = hw
2301
2302       FIP = 0.5*rho(ktopw)*Cstar2**3*heff*2*sqrt(sigmaw*wdens*3.14)
2303       FIP = alpk * FIP
2304Cjyg2
2305       ELSE
2306         FIP = 0.
2307       ENDIF
2308
2309
2310C   Limitation de sigmaw
2311c
2312C   sécurité : si le wake occuppe plus de 90 % de la surface de la maille,
2313C              alors il disparait en se mélangeant à la partie undisturbed
2314
2315      IF ((sigmaw.GT.0.9).or.
2316     .     ((wape.ge.wape2).and.(wape2.le.1.0))) THEN
2317c      IF (sigmaw.GT.0.9) THEN
2318        DO k = 1,klev
2319          dtls(k) = 0.
2320          dqls(k) = 0.
2321          deltatw(k) = 0.
2322          deltaqw(k) = 0.
2323        ENDDO
2324        wape = 0.
2325        hw = hwmin
2326        sigmaw = sigmad
2327        fip = 0.
2328      ENDIF
2329
2330      RETURN
2331      END
2332
2333
2334
Note: See TracBrowser for help on using the repository browser.