source: LMDZ5/branches/LMDZ5_SPLA/libf/phylmd/wake.F90

Last change on this file was 1992, checked in by lguez, 11 years ago

Converted to free source form files in libf/phylmd which were still in
fixed source form. The conversion was done using the polish mode of
the NAG Fortran Compiler.

In addition to converting to free source form, the processing of the
files also:

-- indented the code (including comments);

-- set Fortran keywords to uppercase, and set all other identifiers
to lower case;

-- added qualifiers to end statements (for example "end subroutine
conflx", instead of "end");

-- changed the terminating statements of all DO loops so that each
loop ends with an ENDDO statement (instead of a labeled continue).

-- replaced #include by include.
  • 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 1992 2014-03-05 13:19:12Z evignon $
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.