source: LMDZ5/trunk/libf/phylmd/wake.F @ 1907

Last change on this file since 1907 was 1907, checked in by lguez, 10 years ago

Added a copyright property to every file of the distribution, except
for the fcm files (which have their own copyright). Use svn propget on
a file to see the copyright. For instance:

$ svn propget copyright libf/phylmd/physiq.F90
Name of program: LMDZ
Creation date: 1984
Version: LMDZ5
License: CeCILL version 2
Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
See the license file in the root directory

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