source: LMDZ5/trunk/libf/phylmd/wake.F90 @ 2093

Last change on this file since 2093 was 2078, checked in by idelkadi, 10 years ago

Rajout des directives OPENMP pour la lecture des cles et parametres de convection et wakes dans les fichier conv_param.data, ep_param.data et wake_param.data

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