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

Last change on this file since 1098 was 1059, checked in by lmdzadmin, 16 years ago

Ajout condition supplementaire disparition poche froide (si elle n'atteint pas au - 2eme niveau)
NR/JYG/IM
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 67.0 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)).or.
1247     $      (ktopw(i).le.2)) THEN
1248cIM cf NR/JYG 251108  $     ((wape(i).ge.wape2(i)).and.(wape2(i).le.1.0))) THEN
1249ccc      IF (sigmaw(i).GT.0.9) THEN
1250          dtls(i,k) = 0.
1251          dqls(i,k) = 0.
1252          deltatw(i,k) = 0.
1253          deltaqw(i,k) = 0.
1254        ENDIF
1255       ENDDO
1256      ENDDO
1257c
1258      DO i=1, klon
1259         IF ((sigmaw(i).GT.0.9).or.
1260     $     ((wape(i).ge.wape2(i)).and.(wape2(i).le.1.0))) THEN
1261ccc      IF (sigmaw(i).GT.0.9) THEN
1262         wape(i) = 0.
1263         hw(i) = hwmin
1264         sigmaw(i) = sigmad
1265         fip(i) = 0.
1266        ELSE
1267         wape(i) = wape2(i)
1268        ENDIF
1269      ENDDO
1270c
1271c
1272      RETURN
1273      END
1274      Subroutine WAKE_scal (p,ph,ppi,dtime,sigd_con
1275     :                ,te0,qe0,omgb
1276     :                ,dtdwn,dqdwn,amdwn,amup,dta,dqa
1277     :                ,wdtPBL,wdqPBL,udtPBL,udqPBL
1278     o                ,deltatw,deltaqw,dth,hw,sigmaw,wape,fip,gfl
1279     o                ,dtls,dqls
1280     o                ,ktopw,omgbdth,dp_omgb,wdens
1281     o                ,tu,qu
1282     o                ,dtKE,dqKE
1283     o                ,dtPBL,dqPBL
1284     o                ,omg,dp_deltomg,spread
1285     o                ,Cstar,d_deltat_gw
1286     o                ,d_deltatw2,d_deltaqw2)
1287
1288***************************************************************
1289*                                                             *
1290* WAKE                                                        *
1291*      retour a un Pupper fixe                                *
1292*                                                             *
1293* written by   :  GRANDPEIX Jean-Yves   09/03/2000            *
1294* modified by :   ROEHRIG Romain        01/29/2007            *
1295***************************************************************
1296c
1297      USE dimphy
1298      IMPLICIT none
1299c============================================================================
1300C
1301C
1302C   But : Decrire le comportement des poches froides apparaissant dans les
1303C        grands systemes convectifs, et fournir l'energie disponible pour
1304C        le declenchement de nouvelles colonnes convectives.
1305C
1306C   Variables d'etat : deltatw    : ecart de temperature wake-undisturbed area
1307C                      deltaqw    : ecart d'humidite wake-undisturbed area
1308C                      sigmaw     : fraction d'aire occupee par la poche.
1309C
1310C   Variable de sortie :
1311c
1312c                        wape : WAke Potential Energy
1313c                        fip  : Front Incident Power (W/m2) - ALP
1314c                        gfl  : Gust Front Length per unit area (m-1)
1315C                        dtls : large scale temperature tendency due to wake
1316C                        dqls : large scale humidity tendency due to wake
1317C                        hw   : hauteur de la poche
1318C                     dp_omgb : vertical gradient of large scale omega
1319C                      omgbdth: flux of Delta_Theta transported by LS omega
1320C                      dtKE   : differential heating (wake - unpertubed)
1321C                      dqKE   : differential moistening (wake - unpertubed)
1322C                      omg    : Delta_omg =vertical velocity diff. wake-undist. (Pa/s)
1323C                 dp_deltomg  : vertical gradient of omg (s-1)
1324C                     spread  : spreading term in dt_wake and dq_wake
1325C                 deltatw     : updated temperature difference (T_w-T_u).
1326C                 deltaqw     : updated humidity difference (q_w-q_u).
1327C                 sigmaw      : updated wake fractional area.
1328C                 d_deltat_gw : delta T tendency due to GW
1329c
1330C   Variables d'entree :
1331c
1332c                        aire : aire de la maille
1333c                        te0  : temperature dans l'environnement  (K)
1334C                        qe0  : humidite dans l'environnement     (kg/kg)
1335C                        omgb : vitesse verticale moyenne sur la maille (Pa/s)
1336C                        dtdwn: source de chaleur due aux descentes (K/s)
1337C                        dqdwn: source d'humidite due aux descentes (kg/kg/s)
1338C                        dta  : source de chaleur due courants satures et detrain  (K/s)
1339C                        dqa  : source d'humidite due aux courants satures et detra (kg/kg/s)
1340C                        amdwn: flux de masse total des descentes, par unite de
1341C                                surface de la maille (kg/m2/s)
1342C                        amup : flux de masse total des ascendances, par unite de
1343C                                surface de la maille (kg/m2/s)
1344C                        p    : pressions aux milieux des couches (Pa)
1345C                        ph   : pressions aux interfaces (Pa)
1346C                        ppi  : (p/p_0)**kapa (adim)
1347C                        dtime: increment temporel (s)
1348c
1349C   Variables internes :
1350c
1351c                        rhow : masse volumique de la poche froide
1352C                        rho  : environment density at P levels
1353C                        rhoh : environment density at Ph levels
1354C                        te   : environment temperature | may change within
1355C                        qe   : environment humidity    | sub-time-stepping
1356C                        the  : environment potential temperature
1357C                        thu  : potential temperature in undisturbed area
1358C                        tu   :  temperature  in undisturbed area
1359C                        qu   : humidity in undisturbed area
1360C                      dp_omgb: vertical gradient og LS omega
1361C                      omgbw  : wake average vertical omega
1362C                     dp_omgbw: vertical gradient of omgbw
1363C                      omgbdq : flux of Delta_q transported by LS omega
1364C                        dth  : potential temperature diff. wake-undist.
1365C                        th1  : first pot. temp. for vertical advection (=thu)
1366C                        th2  : second pot. temp. for vertical advection (=thw)
1367C                        q1   : first humidity for vertical advection
1368C                        q2   : second humidity for vertical advection
1369C                     d_deltatw   : terme de redistribution pour deltatw
1370C                     d_deltaqw   : terme de redistribution pour deltaqw
1371C                      deltatw0   : deltatw initial
1372C                      deltaqw0   : deltaqw initial
1373C                      hw0    : hw initial
1374C                      sigmaw0: sigmaw initial
1375C                      amflux : horizontal mass flux through wake boundary
1376C                      wdens  : number of wakes per unit area (3D) or per
1377C                               unit length (2D)
1378C                      Tgw    : 1 sur la période de onde de gravité
1379c                      Cgw    : vitesse de propagation de onde de gravité
1380c                      LL     : distance entre 2 poches
1381
1382c-------------------------------------------------------------------------
1383c          Déclaration de variables
1384c-------------------------------------------------------------------------
1385
1386#include "dimensions.h"
1387cccc#include "dimphy.h"
1388#include "YOMCST.h"
1389#include "cvthermo.h"
1390#include "iniprint.h"
1391
1392c Arguments en entree
1393c--------------------
1394
1395      REAL p(klev),ph(klev+1),ppi(klev)
1396      REAL dtime
1397      REAL te0(klev),qe0(klev)
1398      REAL omgb(klev+1)
1399      REAL dtdwn(klev), dqdwn(klev)
1400      REAL wdtPBL(klev),wdqPBL(klev)
1401      REAL udtPBL(klev),udqPBL(klev)
1402      REAL amdwn(klev), amup(klev)
1403      REAL dta(klev), dqa(klev)
1404      REAL sigd_con
1405
1406c Sorties
1407c--------
1408
1409      REAL deltatw(klev), deltaqw(klev), dth(klev)
1410      REAL tu(klev), qu(klev)
1411      REAL dtls(klev), dqls(klev)
1412      REAL dtKE(klev), dqKE(klev)
1413      REAL dtPBL(klev), dqPBL(klev)
1414      REAL spread(klev)
1415      REAL d_deltatgw(klev)
1416      REAL d_deltatw2(klev), d_deltaqw2(klev)
1417      REAL omgbdth(klev+1), omg(klev+1)
1418      REAL dp_omgb(klev), dp_deltomg(klev)
1419      REAL d_deltat_gw(klev)
1420      REAL hw, sigmaw, wape, fip, gfl, Cstar
1421      INTEGER ktopw
1422
1423c Variables internes
1424c-------------------
1425
1426c Variables à fixer
1427      REAL ALON
1428      REAL coefgw
1429      REAL wdens0, wdens
1430      REAL stark
1431      REAL alpk
1432      REAL delta_t_min
1433      REAL Pupper
1434      INTEGER nsub
1435      REAL dtimesub
1436      REAL sigmad, hwmin
1437
1438c Variables de sauvegarde
1439      REAL deltatw0(klev)
1440      REAL deltaqw0(klev)
1441      REAL te(klev), qe(klev)
1442      REAL sigmaw0, sigmaw1
1443
1444c Variables pour les GW
1445      REAL LL
1446      REAL N2(klev)
1447      REAL Cgw(klev)
1448      REAL Tgw(klev)
1449
1450c Variables liées au calcul de hw
1451      REAL ptop_provis, ptop, ptop_new
1452      REAL sum_dth
1453      REAL dthmin
1454      REAL z, dz, hw0
1455      INTEGER ktop, kupper
1456
1457c Autres variables internes
1458      INTEGER isubstep, k
1459
1460      REAL sum_thu, sum_tu, sum_qu,sum_thvu
1461      REAL sum_dq, sum_rho
1462      REAL sum_dtdwn, sum_dqdwn
1463      REAL av_thu, av_tu, av_qu, av_thvu
1464      REAL av_dth, av_dq, av_rho
1465      REAL av_dtdwn, av_dqdwn
1466
1467      REAL rho(klev), rhoh(klev+1), rhow(klev)
1468      REAL rhow_moyen(klev)
1469      REAL zh(klev), zhh(klev+1)
1470      REAL epaisseur1(klev), epaisseur2(klev)
1471
1472      REAL the(klev), thu(klev)
1473
1474      REAL d_deltatw(klev), d_deltaqw(klev)
1475
1476      REAL omgbw(klev+1), omgtop
1477      REAL dp_omgbw(klev)
1478      REAL ztop, dztop
1479      REAL alpha_up(klev)
1480     
1481      REAL RRe1, RRe2, RRd1, RRd2
1482      REAL Th1(klev), Th2(klev), q1(klev), q2(klev)
1483      REAL D_Th1(klev), D_Th2(klev), D_dth(klev)
1484      REAL D_q1(klev), D_q2(klev), D_dq(klev)
1485      REAL omgbdq(klev)
1486
1487      REAL ff, gg
1488      REAL wape2, Cstar2, heff
1489
1490      REAL Crep(klev)
1491      REAL Crep_upper, Crep_sol
1492
1493C-------------------------------------------------------------------------
1494c         Initialisations
1495c-------------------------------------------------------------------------
1496
1497c      print*, 'wake initialisations'
1498
1499c   Essais d'initialisation avec sigmaw = 0.02 et hw = 10.
1500c-------------------------------------------------------------------------
1501
1502      DATA sigmad, hwmin /.02,10./
1503
1504C Longueur de maille (en m)
1505c-------------------------------------------------------------------------
1506
1507c      ALON = 3.e5
1508      ALON = 1.e6
1509
1510
1511C Configuration de coefgw,stark,wdens (22/02/06 by YU Jingmei)
1512c
1513c      coefgw : Coefficient pour les ondes de gravité
1514c       stark : Coefficient k dans Cstar=k*sqrt(2*WAPE)
1515c       wdens : Densité de poche froide par maille
1516c-------------------------------------------------------------------------
1517
1518      coefgw=10
1519c      coefgw=1
1520c      wdens0 = 1.0/(alon**2)   
1521      wdens = 1.0/(alon**2)       
1522      stark = 0.50
1523cCRtest
1524      alpk=0.1
1525c      alpk = 1.0
1526c      alpk = 0.5
1527c      alpk = 0.05
1528      Crep_upper=0.9
1529      Crep_sol=1.0
1530
1531
1532C Minimum value for |T_wake - T_undist|. Used for wake top definition
1533c-------------------------------------------------------------------------
1534
1535      delta_t_min = 0.2
1536
1537
1538C 1. - Save initial values and initialize tendencies
1539C --------------------------------------------------
1540
1541      DO k=1,klev
1542        deltatw0(k) = deltatw(k)
1543        deltaqw0(k)= deltaqw(k)
1544        te(k) = te0(k)
1545        qe(k) = qe0(k)
1546        dtls(k) = 0.
1547        dqls(k) = 0.
1548        d_deltat_gw(k)=0.
1549        d_deltatw2(k)=0.
1550        d_deltaqw2(k)=0.
1551      ENDDO
1552c      sigmaw1=sigmaw
1553c      IF (sigd_con.GT.sigmaw1) THEN
1554c      print*, 'sigmaw,sigd_con', sigmaw, sigd_con
1555c      ENDIF
1556      sigmaw = max(sigmaw,sigd_con)
1557      sigmaw = max(sigmaw,sigmad)
1558      sigmaw = min(sigmaw,0.99)
1559      sigmaw0 = sigmaw
1560c      wdens=wdens0/(10.*sigmaw)
1561c      IF (sigd_con.GT.sigmaw1) THEN
1562c      print*, 'sigmaw1,sigd1', sigmaw, sigd_con
1563c      ENDIF
1564
1565C 2. - Prognostic part
1566C =========================================================
1567
1568c      print *, 'prognostic wake computation'
1569
1570
1571C 2.1 - Undisturbed area and Wake integrals
1572C ---------------------------------------------------------
1573
1574      z = 0.
1575      ktop=0
1576      kupper = 0
1577      sum_thu = 0.
1578      sum_tu = 0.
1579      sum_qu = 0.
1580      sum_thvu = 0.
1581      sum_dth = 0.
1582      sum_dq = 0.
1583      sum_rho = 0.
1584      sum_dtdwn = 0.
1585      sum_dqdwn = 0.
1586
1587      av_thu = 0.
1588      av_tu =0.
1589      av_qu =0.
1590      av_thvu = 0.
1591      av_dth = 0.
1592      av_dq = 0.
1593      av_rho =0.
1594      av_dtdwn =0.
1595      av_dqdwn = 0.
1596
1597C Potential temperatures and humidity
1598c----------------------------------------------------------
1599
1600      DO k =1,klev
1601        rho(k) = p(k)/(rd*te(k))
1602        IF(k .eq. 1) THEN
1603          rhoh(k) = ph(k)/(rd*te(k))
1604          zhh(k)=0
1605        ELSE
1606          rhoh(k) = ph(k)*2./(rd*(te(k)+te(k-1)))
1607          zhh(k)=(ph(k)-ph(k-1))/(-rhoh(k)*RG)+zhh(k-1)
1608        ENDIF
1609        the(k) = te(k)/ppi(k)
1610        thu(k) = (te(k) - deltatw(k)*sigmaw)/ppi(k)
1611        tu(k) = te(k) - deltatw(k)*sigmaw
1612        qu(k)  =  qe(k) - deltaqw(k)*sigmaw
1613        rhow(k) = p(k)/(rd*(te(k)+deltatw(k)))
1614        dth(k) = deltatw(k)/ppi(k)
1615        LL = (1-sqrt(sigmaw))/sqrt(wdens)       
1616      ENDDO
1617       
1618      DO k = 1, klev-1
1619        IF(k.eq.1) THEN
1620          N2(k)=0
1621        ELSE
1622          N2(k)=max(0.,-RG**2/the(k)*rho(k)*(the(k+1)-the(k-1))
1623     $           /(p(k+1)-p(k-1)))
1624        ENDIF
1625        ZH(k)=(zhh(k)+zhh(k+1))/2
1626
1627        Cgw(k)=sqrt(N2(k))*ZH(k)
1628        Tgw(k)=coefgw*Cgw(k)/LL
1629      ENDDO
1630         
1631      N2(klev)=0
1632      ZH(klev)=0
1633      Cgw(klev)=0
1634      Tgw(klev)=0
1635
1636c  Calcul de la masse volumique moyenne de la colonne
1637c-----------------------------------------------------------------
1638
1639      DO k=1,klev
1640        epaisseur1(k)=0.
1641        epaisseur2(k)=0.
1642      ENDDO
1643
1644      epaisseur1(1)= -(Ph(2)-Ph(1))/(rho(1)*rg)+1.
1645      epaisseur2(1)= -(Ph(2)-Ph(1))/(rho(1)*rg)+1.
1646      rhow_moyen(1) = rhow(1)
1647
1648      DO k = 2, klev
1649        epaisseur1(k)= -(Ph(k+1)-Ph(k))/(rho(k)*rg) +1.
1650        epaisseur2(k)=epaisseur2(k-1)+epaisseur1(k)
1651        rhow_moyen(k) = (rhow_moyen(k-1)*epaisseur2(k-1)+
1652     $                 rhow(k)*epaisseur1(k))/epaisseur2(k)
1653      ENDDO
1654
1655
1656C Choose an integration bound well above wake top
1657c-----------------------------------------------------------------
1658
1659c       Pupper = 50000.  ! melting level
1660       Pupper = 60000.
1661c       Pupper = 70000.
1662
1663
1664C    Determine Wake top pressure (Ptop) from buoyancy integral
1665C-----------------------------------------------------------------
1666
1667c-1/ Pressure of the level where dth becomes less than delta_t_min.
1668
1669      Ptop_provis=ph(1)
1670      DO k= 2,klev
1671        IF (dth(k) .GT. -delta_t_min .and.
1672     $      dth(k-1).LT. -delta_t_min) THEN
1673          Ptop_provis = ((dth(k)+delta_t_min)*p(k-1)
1674     $          - (dth(k-1)+delta_t_min)*p(k)) /(dth(k) - dth(k-1))
1675          GO TO 25
1676        ENDIF
1677      ENDDO
167825    CONTINUE
1679
1680c-2/ dth integral
1681
1682      sum_dth = 0.
1683      dthmin = -delta_t_min
1684      z = 0.
1685
1686      DO k = 1,klev
1687        dz = -(max(ph(k+1),Ptop_provis)-Ph(k))/(rho(k)*rg)
1688        IF (dz .le. 0) GO TO 40
1689        z = z+dz
1690        sum_dth = sum_dth + dth(k)*dz
1691        dthmin = min(dthmin,dth(k))
1692      ENDDO
169340    CONTINUE
1694
1695c-3/ height of triangle with area= sum_dth and base = dthmin
1696
1697      hw0 = 2.*sum_dth/min(dthmin,-0.5)
1698      hw0 = max(hwmin,hw0)
1699
1700c-4/ now, get Ptop
1701
1702      z = 0.
1703      ptop = ph(1)
1704
1705      DO k = 1,klev
1706        dz = min(-(ph(k+1)-Ph(k))/(rho(k)*rg),hw0-z)
1707        IF (dz .le. 0) GO TO 45
1708        z = z+dz
1709        Ptop = Ph(k)-rho(k)*rg*dz
1710      ENDDO
171145    CONTINUE
1712
1713
1714C-5/ Determination de ktop et kupper
1715
1716      DO k=klev,1,-1
1717        IF (ph(k+1) .lt. ptop) ktop=k
1718        IF (ph(k+1) .lt. pupper) kupper=k
1719      ENDDO
1720
1721c-6/ Correct ktop and ptop
1722
1723      Ptop_new=ptop
1724      DO k= ktop,2,-1
1725        IF (dth(k) .GT. -delta_t_min .and.
1726     $      dth(k-1).LT. -delta_t_min) THEN
1727          Ptop_new = ((dth(k)+delta_t_min)*p(k-1)
1728     $          - (dth(k-1)+delta_t_min)*p(k)) /(dth(k) - dth(k-1))
1729          GO TO 225
1730        ENDIF
1731      ENDDO
1732225   CONTINUE
1733
1734      ptop = ptop_new
1735
1736      DO k=klev,1,-1
1737        IF (ph(k+1) .lt. ptop) ktop=k
1738      ENDDO
1739
1740c Set deltatw & deltaqw to 0 above kupper
1741c-----------------------------------------------------------
1742
1743      DO k = kupper,klev
1744        deltatw(k) = 0.
1745        deltaqw(k) = 0.
1746      ENDDO
1747
1748
1749C Vertical gradient of LS omega
1750C------------------------------------------------------------
1751
1752      DO k = 1,kupper
1753        dp_omgb(k) = (omgb(k+1) - omgb(k))/(ph(k+1)-ph(k))
1754      ENDDO
1755
1756
1757C Integrals (and wake top level number)
1758C -----------------------------------------------------------
1759
1760C Initialize sum_thvu to 1st level virt. pot. temp.
1761
1762      z = 1.
1763      dz = 1.
1764      sum_thvu =  thu(1)*(1.+eps*qu(1))*dz
1765      sum_dth = 0.
1766
1767      DO k = 1,klev
1768        dz = -(max(ph(k+1),Ptop)-Ph(k))/(rho(k)*rg)
1769        IF (dz .LE. 0) GO TO 50
1770        z = z+dz
1771        sum_thu = sum_thu + thu(k)*dz
1772        sum_tu = sum_tu + tu(k)*dz
1773        sum_qu = sum_qu + qu(k)*dz
1774        sum_thvu = sum_thvu + thu(k)*(1.+eps*qu(k))*dz
1775        sum_dth = sum_dth + dth(k)*dz
1776        sum_dq = sum_dq + deltaqw(k)*dz
1777        sum_rho = sum_rho + rhow(k)*dz
1778        sum_dtdwn = sum_dtdwn + dtdwn(k)*dz
1779        sum_dqdwn = sum_dqdwn + dqdwn(k)*dz
1780      ENDDO
178150    CONTINUE
1782
1783      hw0 = z
1784
1785C 2.1 - WAPE and mean forcing computation
1786C-------------------------------------------------------------
1787
1788C Means
1789
1790      av_thu = sum_thu/hw0
1791      av_tu = sum_tu/hw0
1792      av_qu = sum_qu/hw0
1793      av_thvu = sum_thvu/hw0
1794c      av_thve = sum_thve/hw0
1795      av_dth = sum_dth/hw0
1796      av_dq = sum_dq/hw0
1797      av_rho = sum_rho/hw0
1798      av_dtdwn = sum_dtdwn/hw0
1799      av_dqdwn = sum_dqdwn/hw0
1800
1801      wape = - rg*hw0*(av_dth
1802     $     + eps*(av_thu*av_dq+av_dth*av_qu+av_dth*av_dq ))/av_thvu
1803
1804C 2.2 Prognostic variable update
1805C ------------------------------------------------------------
1806
1807C Filter out bad wakes
1808
1809      IF ( wape .LT. 0.) THEN
1810        if(prt_level.ge.10) print*,'wape<0'
1811        wape = 0.
1812        hw = hwmin
1813        sigmaw = max(sigmad,sigd_con)
1814        fip = 0.
1815        DO k = 1,klev
1816          deltatw(k) = 0.
1817          deltaqw(k) = 0.
1818          dth(k) = 0.
1819        ENDDO
1820      ELSE
1821        if(prt_level.ge.10) print*,'wape>0'
1822        Cstar = stark*sqrt(2.*wape)
1823      ENDIF
1824
1825C------------------------------------------------------------------
1826C    Sub-time-stepping
1827C------------------------------------------------------------------
1828
1829c      nsub=36
1830      nsub=10
1831      dtimesub=dtime/nsub
1832
1833c------------------------------------------------------------
1834      DO isubstep = 1,nsub
1835c------------------------------------------------------------
1836
1837c        print*,'---------------','substep=',isubstep,'-------------'
1838
1839c  Evolution of sigmaw
1840
1841
1842        gfl = 2.*sqrt(3.14*wdens*sigmaw)           
1843
1844        sigmaw =sigmaw + gfl*Cstar*dtimesub
1845        sigmaw =min(sigmaw,0.99)     !!!!!!!!
1846c        wdens = wdens0/(10.*sigmaw)
1847c        sigmaw =max(sigmaw,sigd_con)
1848c        sigmaw =max(sigmaw,sigmad)
1849
1850c calcul de la difference de vitesse verticale poche - zone non perturbee
1851
1852        z= 0.
1853        dp_deltomg(1:klev)=0.
1854        omg(1:klev+1)=0.
1855
1856        omg(1) = 0.
1857        dp_deltomg(1) = -(gfl*Cstar)/(sigmaw * (1-sigmaw))
1858
1859        DO k=2,ktop
1860          dz = -(Ph(k)-Ph(k-1))/(rho(k-1)*rg)
1861          z = z+dz
1862          dp_deltomg(k)= dp_deltomg(1)
1863          omg(k)= dp_deltomg(1)*z
1864        ENDDO
1865
1866        dztop=-(Ptop-Ph(ktop))/(rho(ktop)*rg)
1867        ztop = z+dztop
1868        omgtop=dp_deltomg(1)*ztop
1869
1870
1871c Conversion de la vitesse verticale de m/s a Pa/s
1872
1873        omgtop = -rho(ktop)*rg*omgtop
1874        dp_deltomg(1) = omgtop/(ptop-ph(1))
1875
1876        DO k = 1,ktop
1877          omg(k) = - rho(k)*rg*omg(k)
1878          dp_deltomg(k) = dp_deltomg(1)
1879        ENDDO
1880
1881c   raccordement lineaire de omg de ptop a pupper
1882
1883      IF (kupper .GT. ktop) THEN
1884        omg(kupper+1) = - Rg*amdwn(kupper+1)/sigmaw
1885     $                + Rg*amup(kupper+1)/(1.-sigmaw)
1886        dp_deltomg(kupper) = (omgtop-omg(kupper+1))/(Ptop-Pupper)
1887        DO k=ktop+1,kupper
1888          dp_deltomg(k) = dp_deltomg(kupper)
1889          omg(k) = omgtop+(ph(k)-Ptop)*dp_deltomg(kupper)
1890        ENDDO
1891      ENDIF
1892
1893c   Compute wake average vertical velocity omgbw
1894
1895      DO k = 1,klev+1
1896        omgbw(k) = omgb(k)+(1.-sigmaw)*omg(k)
1897      ENDDO
1898
1899c  and its vertical gradient dp_omgbw
1900
1901      DO k = 1,klev
1902        dp_omgbw(k) = (omgbw(k+1)-omgbw(k))/(ph(k+1)-ph(k))
1903      ENDDO
1904
1905
1906c  Upstream coefficients for omgb velocity
1907c--    (alpha_up(k) is the coefficient of the value at level k)
1908c--    (1-alpha_up(k) is the coefficient of the value at level k-1)
1909
1910      DO k = 1,klev
1911       alpha_up(k) = 0.
1912       IF (omgb(k) .GT. 0.) alpha_up(k) = 1.
1913      ENDDO
1914
1915c  Matrix expressing [The,deltatw] from  [Th1,Th2]
1916
1917      RRe1 = 1.-sigmaw
1918      RRe2 = sigmaw
1919      RRd1 = -1.
1920      RRd2 = 1.
1921
1922c Get [Th1,Th2], dth and [q1,q2]
1923
1924      DO k = 1,kupper+1
1925        dth(k) = deltatw(k)/ppi(k)
1926        Th1(k) = the(k) - sigmaw     *dth(k)   ! undisturbed area
1927        Th2(k) = the(k) + (1.-sigmaw)*dth(k)   ! wake
1928        q1(k) = qe(k) - sigmaw     *deltaqw(k) ! undisturbed area
1929        q2(k) = qe(k) + (1.-sigmaw)*deltaqw(k) ! wake
1930      ENDDO
1931
1932      D_Th1(1) = 0.
1933      D_Th2(1) = 0.
1934      D_dth(1) = 0.
1935      D_q1(1) = 0.
1936      D_q2(1) = 0.
1937      D_dq(1) = 0.
1938
1939      DO k = 2,kupper+1 !   loop on interfaces
1940        D_Th1(k) = Th1(k-1)-Th1(k)
1941        D_Th2(k) = Th2(k-1)-Th2(k)
1942        D_dth(k) = dth(k-1)-dth(k)
1943        D_q1(k) = q1(k-1)-q1(k)
1944        D_q2(k) = q2(k-1)-q2(k)
1945        D_dq(k) = deltaqw(k-1)-deltaqw(k)
1946      ENDDO
1947
1948      omgbdth(1) = 0.
1949      omgbdq(1) = 0.
1950
1951      DO k = 2,kupper+1  !   loop on interfaces
1952        omgbdth(k) = omgb(k)*(    dth(k-1) -     dth(k))
1953        omgbdq(k)  = omgb(k)*(deltaqw(k-1) - deltaqw(k))
1954      ENDDO
1955
1956
1957c-----------------------------------------------------------------
1958      DO k=1,kupper-1
1959c-----------------------------------------------------------------
1960c
1961c   Compute redistribution (advective) term
1962c
1963         d_deltatw(k) =
1964     $             dtimesub/(Ph(k)-Ph(k+1))*(
1965     $       RRd1*omg(k  )*sigmaw     *D_Th1(k)
1966     $      -RRd2*omg(k+1)*(1.-sigmaw)*D_Th2(k+1)
1967     $      -(1.-alpha_up(k))*omgbdth(k) - alpha_up(k+1)*omgbdth(k+1)
1968     $                      )*ppi(k)
1969c         print*,'d_deltatw=',d_deltatw(k)
1970c
1971         d_deltaqw(k) =
1972     $             dtimesub/(Ph(k)-Ph(k+1))*(
1973     $       RRd1*omg(k  )*sigmaw     *D_q1(k)
1974     $      -RRd2*omg(k+1)*(1.-sigmaw)*D_q2(k+1)
1975     $      -(1.-alpha_up(k))*omgbdq(k) - alpha_up(k+1)*omgbdq(k+1)
1976     $                      )
1977c         print*,'d_deltaqw=',d_deltaqw(k)
1978c
1979c   and increment large scale tendencies
1980c
1981         dtls(k) = dtls(k) +
1982     $               dtimesub*(
1983     $        ( RRe1*omg(k  )*sigmaw     *D_Th1(k)
1984     $         -RRe2*omg(k+1)*(1.-sigmaw)*D_Th2(k+1) )
1985     $               /(Ph(k)-Ph(k+1))
1986     $         -sigmaw*(1.-sigmaw)*dth(k)*dp_deltomg(k)
1987     $                      )*ppi(k)
1988c         print*,'dtls=',dtls(k)
1989c
1990         dqls(k) = dqls(k) +
1991     $               dtimesub*(
1992     $        ( RRe1*omg(k  )*sigmaw     *D_q1(k)
1993     $         -RRe2*omg(k+1)*(1.-sigmaw)*D_q2(k+1) )
1994     $               /(Ph(k)-Ph(k+1))
1995     $         -sigmaw*(1.-sigmaw)*deltaqw(k)*dp_deltomg(k)
1996     $                      )
1997c         print*,'dqls=',dqls(k)
1998
1999c-------------------------------------------------------------------
2000      ENDDO
2001c------------------------------------------------------------------
2002
2003C   Increment state variables
2004
2005      DO k = 1,kupper-1
2006
2007c Coefficient de répartition
2008
2009        Crep(k)=Crep_sol*(ph(kupper)-ph(k))/(ph(kupper)-ph(1))
2010        Crep(k)=Crep(k)+Crep_upper*(ph(1)-ph(k))/(p(1)-ph(kupper))
2011       
2012
2013c Reintroduce compensating subsidence term.
2014
2015c        dtKE(k)=(dtdwn(k)*Crep(k))/sigmaw
2016c        dtKE(k)=dtKE(k)-(dtdwn(k)*(1-Crep(k))+dta(k))
2017c     .                   /(1-sigmaw)
2018c        dqKE(k)=(dqdwn(k)*Crep(k))/sigmaw
2019c        dqKE(k)=dqKE(k)-(dqdwn(k)*(1-Crep(k))+dqa(k))
2020c     .                   /(1-sigmaw)
2021c
2022c        dtKE(k)=(dtdwn(k)*Crep(k)+(1-Crep(k))*dta(k))/sigmaw
2023c        dtKE(k)=dtKE(k)-(dtdwn(k)*(1-Crep(k))+dta(k)*Crep(k))
2024c     .                   /(1-sigmaw)
2025c        dqKE(k)=(dqdwn(k)*Crep(k)+(1-Crep(k))*dqa(k))/sigmaw
2026c        dqKE(k)=dqKE(k)-(dqdwn(k)*(1-Crep(k))+dqa(k)*Crep(k))
2027c     .                   /(1-sigmaw)
2028
2029        dtKE(k)=(dtdwn(k)/sigmaw - dta(k)/(1.-sigmaw))
2030        dqKE(k)=(dqdwn(k)/sigmaw - dqa(k)/(1.-sigmaw))
2031c        print*,'dtKE=',dtKE(k)
2032c        print*,'dqKE=',dqKE(k)
2033c
2034        dtPBL(k)=(wdtPBL(k)/sigmaw - udtPBL(k)/(1.-sigmaw))
2035        dqPBL(k)=(wdqPBL(k)/sigmaw - udqPBL(k)/(1.-sigmaw))
2036c
2037        spread(k) = (1.-sigmaw)*dp_deltomg(k)+gfl*Cstar/sigmaw
2038c        print*,'spread=',spread(k)
2039
2040
2041c ajout d'un effet onde de gravité -Tgw(k)*deltatw(k) 03/02/06 YU Jingmei
2042
2043        d_deltat_gw(k)=d_deltat_gw(k)-Tgw(k)*deltatw(k)* dtimesub
2044c        print*,'d_delta_gw=',d_deltat_gw(k)
2045        ff=d_deltatw(k)/dtimesub
2046
2047c Sans GW
2048c
2049c        deltatw(k)=deltatw(k)+dtimesub*(ff+dtKE(k)-spread(k)*deltatw(k))
2050c
2051c GW formule 1
2052c
2053c        deltatw(k) = deltatw(k)+dtimesub*
2054c     $         (ff+dtKE(k) - spread(k)*deltatw(k)-Tgw(k)*deltatw(k))
2055c
2056c GW formule 2
2057
2058        IF (dtimesub*Tgw(k).lt.1.e-10) THEN
2059          deltatw(k) = deltatw(k)+dtimesub*
2060     $          (ff+dtKE(k)+dtPBL(k)
2061     $          - spread(k)*deltatw(k)-Tgw(k)*deltatw(k))
2062        ELSE
2063           deltatw(k) = deltatw(k)+1/Tgw(k)*(1-exp(-dtimesub*Tgw(k)))*
2064     $          (ff+dtKE(k)+dtPBL(k)
2065     $          - spread(k)*deltatw(k)-Tgw(k)*deltatw(k))
2066        ENDIF
2067   
2068        dth(k) = deltatw(k)/ppi(k)
2069
2070        gg=d_deltaqw(k)/dtimesub
2071
2072       deltaqw(k) = deltaqw(k) +
2073     $         dtimesub*(gg+ dqKE(k)+dqPBL(k) - spread(k)*deltaqw(k))
2074
2075       d_deltatw2(k)=d_deltatw2(k)+d_deltatw(k)
2076       d_deltaqw2(k)=d_deltaqw2(k)+d_deltaqw(k)
2077      ENDDO
2078
2079C   And update large scale variables
2080
2081      DO k = 1,kupper
2082        te(k) = te0(k) + dtls(k)
2083        qe(k) = qe0(k) + dqls(k)
2084        the(k) = te(k)/ppi(k)
2085      ENDDO
2086
2087c     Determine Ptop from buoyancy integral
2088c----------------------------------------------------------------------
2089
2090c-1/ Pressure of the level where dth changes sign.
2091
2092      Ptop_provis=ph(1)
2093
2094      DO k= 2,klev
2095        IF (dth(k) .GT. -delta_t_min .and.
2096     $      dth(k-1).LT. -delta_t_min) THEN
2097          Ptop_provis = ((dth(k)+delta_t_min)*p(k-1)
2098     $          - (dth(k-1)+delta_t_min)*p(k)) /(dth(k) - dth(k-1))
2099        GO TO 65
2100        ENDIF
2101      ENDDO
210265    CONTINUE
2103
2104c-2/ dth integral
2105
2106      sum_dth = 0.
2107      dthmin = -delta_t_min
2108      z = 0.
2109
2110      DO k = 1,klev
2111        dz = -(max(ph(k+1),Ptop_provis)-Ph(k))/(rho(k)*rg)
2112        IF (dz .le. 0) GO TO 70
2113        z = z+dz
2114        sum_dth = sum_dth + dth(k)*dz
2115        dthmin = min(dthmin,dth(k))
2116      ENDDO
211770    CONTINUE
2118
2119c-3/ height of triangle with area= sum_dth and base = dthmin
2120
2121      hw = 2.*sum_dth/min(dthmin,-0.5)
2122      hw = max(hwmin,hw)
2123
2124c-4/ now, get Ptop
2125
2126      ktop = 0
2127      z=0.
2128
2129      DO k = 1,klev
2130        dz = min(-(ph(k+1)-Ph(k))/(rho(k)*rg),hw-z)
2131        IF (dz .le. 0) GO TO 75
2132        z = z+dz
2133        Ptop = Ph(k)-rho(k)*rg*dz
2134        ktop = k
2135      ENDDO
213675    CONTINUE
2137
2138c-5/Correct ktop and ptop
2139
2140      Ptop_new=ptop
2141
2142      DO k= ktop,2,-1
2143        IF (dth(k) .GT. -delta_t_min .and.
2144     $      dth(k-1).LT. -delta_t_min) THEN
2145          Ptop_new = ((dth(k)+delta_t_min)*p(k-1)
2146     $          - (dth(k-1)+delta_t_min)*p(k)) /(dth(k) - dth(k-1))
2147          GO TO 275
2148        ENDIF
2149      ENDDO
2150275   CONTINUE
2151
2152      ptop = ptop_new
2153
2154      DO k=klev,1,-1
2155        IF (ph(k+1) .LT. ptop) ktop=k
2156      ENDDO
2157
2158c-6/ Set deltatw & deltaqw to 0 above kupper
2159
2160      DO k = kupper,klev
2161        deltatw(k) = 0.
2162        deltaqw(k) = 0.
2163      ENDDO
2164
2165c------------------------------------------------------------------
2166      ENDDO      ! end sub-timestep loop
2167C -----------------------------------------------------------------
2168
2169c   Get back to tendencies per second
2170
2171      DO k = 1,kupper-1
2172        dtls(k) = dtls(k)/dtime
2173        dqls(k) = dqls(k)/dtime
2174        d_deltatw2(k)=d_deltatw2(k)/dtime
2175        d_deltaqw2(k)=d_deltaqw2(k)/dtime
2176        d_deltat_gw(k) = d_deltat_gw(k)/dtime
2177      ENDDO
2178
2179C 2.1 - Undisturbed area and Wake integrals
2180C ---------------------------------------------------------
2181
2182      z = 0.
2183      sum_thu = 0.
2184      sum_tu = 0.
2185      sum_qu = 0.
2186      sum_thvu = 0.
2187      sum_dth = 0.
2188      sum_dq = 0.
2189      sum_rho = 0.
2190      sum_dtdwn = 0.
2191      sum_dqdwn = 0.
2192
2193      av_thu = 0.
2194      av_tu =0.
2195      av_qu =0.
2196      av_thvu = 0.
2197      av_dth = 0.
2198      av_dq = 0.
2199      av_rho =0.
2200      av_dtdwn =0.
2201      av_dqdwn = 0.
2202
2203C Potential temperatures and humidity
2204c----------------------------------------------------------
2205
2206      DO k =1,klev
2207        rho(k) = p(k)/(rd*te(k))
2208        IF(k .eq. 1) THEN
2209          rhoh(k) = ph(k)/(rd*te(k))
2210          zhh(k)=0
2211        ELSE
2212          rhoh(k) = ph(k)*2./(rd*(te(k)+te(k-1)))
2213          zhh(k)=(ph(k)-ph(k-1))/(-rhoh(k)*RG)+zhh(k-1)
2214        ENDIF
2215        the(k) = te(k)/ppi(k)
2216        thu(k) = (te(k) - deltatw(k)*sigmaw)/ppi(k)
2217        tu(k) = te(k) - deltatw(k)*sigmaw
2218        qu(k)  =  qe(k) - deltaqw(k)*sigmaw
2219        rhow(k) = p(k)/(rd*(te(k)+deltatw(k)))
2220        dth(k) = deltatw(k)/ppi(k)
2221       
2222      ENDDO
2223
2224C Integrals (and wake top level number)
2225C -----------------------------------------------------------
2226
2227C Initialize sum_thvu to 1st level virt. pot. temp.
2228
2229      z = 1.
2230      dz = 1.
2231      sum_thvu =  thu(1)*(1.+eps*qu(1))*dz
2232      sum_dth = 0.
2233
2234      DO k = 1,klev
2235        dz = -(max(ph(k+1),Ptop)-Ph(k))/(rho(k)*rg)
2236
2237        IF (dz .LE. 0) GO TO 51
2238        z = z+dz
2239        sum_thu = sum_thu + thu(k)*dz
2240        sum_tu = sum_tu + tu(k)*dz
2241        sum_qu = sum_qu + qu(k)*dz
2242        sum_thvu = sum_thvu + thu(k)*(1.+eps*qu(k))*dz
2243        sum_dth = sum_dth + dth(k)*dz
2244        sum_dq = sum_dq + deltaqw(k)*dz
2245        sum_rho = sum_rho + rhow(k)*dz
2246        sum_dtdwn = sum_dtdwn + dtdwn(k)*dz
2247        sum_dqdwn = sum_dqdwn + dqdwn(k)*dz
2248      ENDDO
2249 51   CONTINUE
2250
2251      hw0 = z
2252
2253C 2.1 - WAPE and mean forcing computation
2254C-------------------------------------------------------------
2255
2256C Means
2257
2258      av_thu = sum_thu/hw0
2259      av_tu = sum_tu/hw0
2260      av_qu = sum_qu/hw0
2261      av_thvu = sum_thvu/hw0
2262      av_dth = sum_dth/hw0
2263      av_dq = sum_dq/hw0
2264      av_rho = sum_rho/hw0
2265      av_dtdwn = sum_dtdwn/hw0
2266      av_dqdwn = sum_dqdwn/hw0
2267
2268      wape2 = - rg*hw0*(av_dth
2269     $     + eps*(av_thu*av_dq+av_dth*av_qu+av_dth*av_dq ))/av_thvu
2270
2271
2272C 2.2 Prognostic variable update
2273C ------------------------------------------------------------
2274
2275C Filter out bad wakes
2276
2277      IF ( wape2 .LT. 0.) THEN
2278        if(prt_level.ge.10) print*,'wape2<0'
2279        wape2 = 0.
2280        hw = hwmin
2281        sigmaw = max(sigmad,sigd_con)
2282        fip = 0.
2283        DO k = 1,klev
2284          deltatw(k) = 0.
2285          deltaqw(k) = 0.
2286          dth(k) = 0.
2287        ENDDO
2288      ELSE
2289        if(prt_level.ge.10) print*,'wape2>0'
2290        Cstar2 = stark*sqrt(2.*wape2)
2291
2292      ENDIF
2293
2294      ktopw = ktop
2295
2296      IF (ktopw .gt. 0) then
2297
2298Cjyg1     Utilisation d'un h_efficace constant ( ~ feeding layer)
2299ccc       heff = 600.
2300C      Utilisation de la hauteur hw
2301cc       heff = 0.7*hw
2302       heff = hw
2303
2304       FIP = 0.5*rho(ktopw)*Cstar2**3*heff*2*sqrt(sigmaw*wdens*3.14)
2305       FIP = alpk * FIP
2306Cjyg2
2307       ELSE
2308         FIP = 0.
2309       ENDIF
2310
2311
2312C   Limitation de sigmaw
2313c
2314C   sécurité : si le wake occuppe plus de 90 % de la surface de la maille,
2315C              alors il disparait en se mélangeant à la partie undisturbed
2316
2317      IF ((sigmaw.GT.0.9).or.
2318     .     ((wape.ge.wape2).and.(wape2.le.1.0)).or.(ktopw.le.2)) THEN
2319cIM cf NR/JYG 251108    .     ((wape.ge.wape2).and.(wape2.le.1.0))) THEN
2320c      IF (sigmaw.GT.0.9) THEN
2321        DO k = 1,klev
2322          dtls(k) = 0.
2323          dqls(k) = 0.
2324          deltatw(k) = 0.
2325          deltaqw(k) = 0.
2326        ENDDO
2327        wape = 0.
2328        hw = hwmin
2329        sigmaw = sigmad
2330        fip = 0.
2331      ENDIF
2332
2333      RETURN
2334      END
2335
2336
2337
Note: See TracBrowser for help on using the repository browser.