source: LMDZ4/branches/LMDZ4-dev-20091210/libf/phylmd/wake.F @ 3655

Last change on this file since 3655 was 1276, checked in by musat, 15 years ago

Corrections conservation de l'eau dans la nouvelle physique
JYG/IM

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