source: LMDZ5/branches/testing/libf/phylmd/wake.F90 @ 2056

Last change on this file since 2056 was 1999, checked in by Laurent Fairhead, 11 years ago

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