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

Last change on this file since 2311 was 2311, checked in by Ehouarn Millour, 9 years ago

Further modifications to enforce physics/dynamics separation:

  • moved iniprint.h and misc_mod back to dyn3d_common, as these should only be used by dynamics.
  • created print_control_mod in the physics to store flags prt_level, lunout, debug to be local to physics (should be used rather than iniprint.h)
  • created abort_physic.F90 , which does the same job as abort_gcm() did, but should be used instead when in physics.
  • reactivated inifis (turned it into a module, inifis_mod.F90) to initialize physical constants and print_control_mod flags.

EM

  • 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: 53.1 KB
Line 
1
2! $Id: wake.F90 2311 2015-06-25 07:45:24Z emillour $
3
4SUBROUTINE wake(p, ph, pi, dtime, sigd_con, &
5                te0, qe0, omgb, &
6                dtdwn, dqdwn, amdwn, amup, dta, dqa, &
7                wdtpbl, wdqpbl, udtpbl, udqpbl, &
8                deltatw, deltaqw, dth, hw, sigmaw, wape, fip, gfl, &
9                dtls, dqls, ktopw, omgbdth, dp_omgb, wdens, tu, qu, &
10                dtke, dqke, dtpbl, dqpbl, omg, dp_deltomg, spread, cstar, &
11                d_deltat_gw, d_deltatw2, d_deltaqw2)
12
13
14  ! **************************************************************
15  ! *
16  ! WAKE                                                        *
17  ! retour a un Pupper fixe                                *
18  ! *
19  ! written by   :  GRANDPEIX Jean-Yves   09/03/2000            *
20  ! modified by :   ROEHRIG Romain        01/29/2007            *
21  ! **************************************************************
22
23  USE dimphy
24  use mod_phys_lmdz_para
25  USE print_control_mod, ONLY: prt_level
26  IMPLICIT NONE
27  ! ============================================================================
28
29
30  ! But : Decrire le comportement des poches froides apparaissant dans les
31  ! grands systemes convectifs, et fournir l'energie disponible pour
32  ! le declenchement de nouvelles colonnes convectives.
33
34  ! Variables d'etat : deltatw    : ecart de temperature wake-undisturbed
35  ! area
36  ! deltaqw    : ecart d'humidite wake-undisturbed area
37  ! sigmaw     : fraction d'aire occupee par la poche.
38
39  ! Variable de sortie :
40
41  ! wape : WAke Potential Energy
42  ! fip  : Front Incident Power (W/m2) - ALP
43  ! gfl  : Gust Front Length per unit area (m-1)
44  ! dtls : large scale temperature tendency due to wake
45  ! dqls : large scale humidity tendency due to wake
46  ! hw   : hauteur de la poche
47  ! dp_omgb : vertical gradient of large scale omega
48  ! wdens   : densite de poches
49  ! omgbdth: flux of Delta_Theta transported by LS omega
50  ! dtKE   : differential heating (wake - unpertubed)
51  ! dqKE   : differential moistening (wake - unpertubed)
52  ! omg    : Delta_omg =vertical velocity diff. wake-undist. (Pa/s)
53  ! dp_deltomg  : vertical gradient of omg (s-1)
54  ! spread  : spreading term in dt_wake and dq_wake
55  ! deltatw     : updated temperature difference (T_w-T_u).
56  ! deltaqw     : updated humidity difference (q_w-q_u).
57  ! sigmaw      : updated wake fractional area.
58  ! d_deltat_gw : delta T tendency due to GW
59
60  ! Variables d'entree :
61
62  ! aire : aire de la maille
63  ! te0  : temperature dans l'environnement  (K)
64  ! qe0  : humidite dans l'environnement     (kg/kg)
65  ! omgb : vitesse verticale moyenne sur la maille (Pa/s)
66  ! dtdwn: source de chaleur due aux descentes (K/s)
67  ! dqdwn: source d'humidite due aux descentes (kg/kg/s)
68  ! dta  : source de chaleur due courants satures et detrain  (K/s)
69  ! dqa  : source d'humidite due aux courants satures et detra (kg/kg/s)
70  ! amdwn: flux de masse total des descentes, par unite de
71  ! surface de la maille (kg/m2/s)
72  ! amup : flux de masse total des ascendances, par unite de
73  ! surface de la maille (kg/m2/s)
74  ! p    : pressions aux milieux des couches (Pa)
75  ! ph   : pressions aux interfaces (Pa)
76  ! pi  : (p/p_0)**kapa (adim)
77  ! dtime: increment temporel (s)
78
79  ! Variables internes :
80
81  ! rhow : masse volumique de la poche froide
82  ! rho  : environment density at P levels
83  ! rhoh : environment density at Ph levels
84  ! te   : environment temperature | may change within
85  ! qe   : environment humidity    | sub-time-stepping
86  ! the  : environment potential temperature
87  ! thu  : potential temperature in undisturbed area
88  ! tu   :  temperature  in undisturbed area
89  ! qu   : humidity in undisturbed area
90  ! dp_omgb: vertical gradient og LS omega
91  ! omgbw  : wake average vertical omega
92  ! dp_omgbw: vertical gradient of omgbw
93  ! omgbdq : flux of Delta_q transported by LS omega
94  ! dth  : potential temperature diff. wake-undist.
95  ! th1  : first pot. temp. for vertical advection (=thu)
96  ! th2  : second pot. temp. for vertical advection (=thw)
97  ! q1   : first humidity for vertical advection
98  ! q2   : second humidity for vertical advection
99  ! d_deltatw   : terme de redistribution pour deltatw
100  ! d_deltaqw   : terme de redistribution pour deltaqw
101  ! deltatw0   : deltatw initial
102  ! deltaqw0   : deltaqw initial
103  ! hw0    : hw initial
104  ! sigmaw0: sigmaw initial
105  ! amflux : horizontal mass flux through wake boundary
106  ! wdens_ref: initial number of wakes per unit area (3D) or per
107  ! unit length (2D), at the beginning of each time step
108  ! Tgw    : 1 sur la période de onde de gravité
109  ! Cgw    : vitesse de propagation de onde de gravité
110  ! LL     : distance entre 2 poches
111
112  ! -------------------------------------------------------------------------
113  ! Déclaration de variables
114  ! -------------------------------------------------------------------------
115
116  include "dimensions.h"
117  include "YOMCST.h"
118  include "cvthermo.h"
119
120  ! Arguments en entree
121  ! --------------------
122
123  REAL, DIMENSION (klon, klev),     INTENT(IN)          :: p, pi
124  REAL, DIMENSION (klon, klev+1),   INTENT(IN)          :: ph, omgb
125  REAL,                             INTENT(IN)          :: dtime
126  REAL, DIMENSION (klon, klev),     INTENT(IN)          :: te0, qe0
127  REAL, DIMENSION (klon, klev),     INTENT(IN)          :: dtdwn, dqdwn
128  REAL, DIMENSION (klon, klev),     INTENT(IN)          :: wdtpbl, wdqpbl, udtpbl, udqpbl ! UNUSED
129  REAL, DIMENSION (klon, klev),     INTENT(IN)          :: amdwn, amup
130  REAL, DIMENSION (klon, klev),     INTENT(IN)          :: dta, dqa
131  REAL, DIMENSION (klon),           INTENT(IN)          :: sigd_con
132
133  !
134  ! Input/Output
135  REAL, DIMENSION (klon, klev),     INTENT(INOUT)       :: deltatw, deltaqw
136  REAL, DIMENSION (klon),           INTENT(INOUT)       :: sigmaw
137
138  ! Sorties
139  ! --------
140
141  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: dth
142  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: tu, qu
143  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: dtls, dqls
144  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: dtke, dqke
145  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: dtpbl, dqpbl
146  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: spread
147  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: d_deltatw2, d_deltaqw2
148  REAL, DIMENSION (klon, klev+1),   INTENT(OUT)         :: omgbdth, omg
149  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: dp_omgb, dp_deltomg
150  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: d_deltat_gw
151  REAL, DIMENSION (klon),           INTENT(OUT)         :: hw, wape, fip, gfl, cstar
152  REAL, DIMENSION (klon),           INTENT(OUT)         :: wdens
153  INTEGER, DIMENSION (klon),        INTENT(OUT)         :: ktopw
154
155  ! Variables internes
156  ! -------------------
157
158  ! Variables à fixer
159  REAL alon
160  LOGICAL, SAVE :: first = .TRUE.
161  !$OMP THREADPRIVATE(first)
162  REAL, SAVE ::  stark, wdens_ref, coefgw, alpk, crep_upper, crep_sol 
163  !$OMP THREADPRIVATE(stark, wdens_ref, coefgw, alpk, crep_upper, crep_sol)
164  REAL delta_t_min
165  INTEGER nsub
166  REAL dtimesub
167  REAL sigmad, hwmin, wapecut
168  REAL :: sigmaw_max
169  REAL :: dens_rate
170  REAL wdens0
171  ! IM 080208
172  LOGICAL, DIMENSION (klon) :: gwake
173
174  ! Variables de sauvegarde
175  REAL, DIMENSION (klon, klev) :: deltatw0
176  REAL, DIMENSION (klon, klev) :: deltaqw0
177  REAL, DIMENSION (klon, klev) :: te, qe
178  REAL, DIMENSION (klon) :: sigmaw0, sigmaw1
179
180  ! Variables pour les GW
181  REAL, DIMENSION (klon) :: ll
182  REAL, DIMENSION (klon, klev) :: n2
183  REAL, DIMENSION (klon, klev) :: cgw
184  REAL, DIMENSION (klon, klev) :: tgw
185
186  ! Variables liées au calcul de hw
187  REAL, DIMENSION (klon) :: ptop_provis, ptop, ptop_new
188  REAL, DIMENSION (klon) :: sum_dth
189  REAL, DIMENSION (klon) :: dthmin
190  REAL, DIMENSION (klon) :: z, dz, hw0
191  INTEGER, DIMENSION (klon) :: ktop, kupper
192
193  ! Sub-timestep tendencies and related variables
194  REAL d_deltatw(klon, klev), d_deltaqw(klon, klev)
195  REAL d_te(klon, klev), d_qe(klon, klev)
196  REAL d_sigmaw(klon), alpha(klon)
197  REAL q0_min(klon), q1_min(klon)
198  LOGICAL wk_adv(klon), ok_qx_qw(klon)
199  REAL epsilon
200  DATA epsilon/1.E-15/
201
202  ! Autres variables internes
203  INTEGER isubstep, k, i
204
205  REAL, DIMENSION (klon) :: sum_thu, sum_tu, sum_qu, sum_thvu
206  REAL, DIMENSION (klon) :: sum_dq, sum_rho
207  REAL, DIMENSION (klon) :: sum_dtdwn, sum_dqdwn
208  REAL, DIMENSION (klon) :: av_thu, av_tu, av_qu, av_thvu
209  REAL, DIMENSION (klon) :: av_dth, av_dq, av_rho
210  REAL, DIMENSION (klon) :: av_dtdwn, av_dqdwn
211
212  REAL, DIMENSION (klon, klev) :: rho, rhow
213  REAL, DIMENSION (klon, klev+1) :: rhoh
214  REAL, DIMENSION (klon, klev) :: rhow_moyen
215  REAL, DIMENSION (klon, klev) :: zh
216  REAL, DIMENSION (klon, klev+1) :: zhh
217  REAL, DIMENSION (klon, klev) :: epaisseur1, epaisseur2
218
219  REAL, DIMENSION (klon, klev) :: the, thu
220
221  ! REAL, DIMENSION(klon,klev) :: d_deltatw, d_deltaqw
222
223  REAL, DIMENSION (klon, klev+1) :: omgbw
224  REAL, DIMENSION (klon) :: pupper
225  REAL, DIMENSION (klon) :: omgtop
226  REAL, DIMENSION (klon, klev) :: dp_omgbw
227  REAL, DIMENSION (klon) :: ztop, dztop
228  REAL, DIMENSION (klon, klev) :: alpha_up
229
230  REAL, DIMENSION (klon) :: rre1, rre2
231  REAL :: rrd1, rrd2
232  REAL, DIMENSION (klon, klev) :: th1, th2, q1, q2
233  REAL, DIMENSION (klon, klev) :: d_th1, d_th2, d_dth
234  REAL, DIMENSION (klon, klev) :: d_q1, d_q2, d_dq
235  REAL, DIMENSION (klon, klev) :: omgbdq
236
237  REAL, DIMENSION (klon) :: ff, gg
238  REAL, DIMENSION (klon) :: wape2, cstar2, heff
239
240  REAL, DIMENSION (klon, klev) :: crep
241
242  REAL, DIMENSION (klon, klev) :: ppi
243
244  ! cc nrlmd
245  REAL, DIMENSION (klon) :: death_rate, nat_rate
246  REAL, DIMENSION (klon, klev) :: entr
247  REAL, DIMENSION (klon, klev) :: detr
248
249  ! -------------------------------------------------------------------------
250  ! Initialisations
251  ! -------------------------------------------------------------------------
252
253  ! print*, 'wake initialisations'
254
255  ! Essais d'initialisation avec sigmaw = 0.02 et hw = 10.
256  ! -------------------------------------------------------------------------
257
258  DATA wapecut, sigmad, hwmin/5., .02, 10./
259  ! cc nrlmd
260  DATA sigmaw_max/0.4/
261  DATA dens_rate/0.1/
262  ! cc
263  ! Longueur de maille (en m)
264  ! -------------------------------------------------------------------------
265
266  ! ALON = 3.e5
267  alon = 1.E6
268
269
270  ! Configuration de coefgw,stark,wdens (22/02/06 by YU Jingmei)
271
272  ! coefgw : Coefficient pour les ondes de gravité
273  ! stark : Coefficient k dans Cstar=k*sqrt(2*WAPE)
274  ! wdens : Densité de poche froide par maille
275  ! -------------------------------------------------------------------------
276
277  ! cc nrlmd      coefgw=10
278  ! coefgw=1
279  ! wdens0 = 1.0/(alon**2)
280  ! cc nrlmd      wdens = 1.0/(alon**2)
281  ! cc nrlmd      stark = 0.50
282  ! CRtest
283  ! cc nrlmd      alpk=0.1
284  ! alpk = 1.0
285  ! alpk = 0.5
286  ! alpk = 0.05
287
288 if (first) then
289  stark = 0.33
290  alpk = 0.25
291  wdens_ref = 8.E-12
292  coefgw = 4.
293  crep_upper = 0.9
294  crep_sol = 1.0
295
296  ! cc nrlmd Lecture du fichier wake_param.data
297 !$OMP MASTER
298  OPEN (99, FILE='wake_param.data', STATUS='old', FORM='formatted', ERR=9999)
299  READ (99, *, END=9998) stark
300  READ (99, *, END=9998) alpk
301  READ (99, *, END=9998) wdens_ref
302  READ (99, *, END=9998) coefgw
3039998 CONTINUE
304  CLOSE (99)
3059999 CONTINUE
306 !$OMP END MASTER
307  CALL bcast(stark)
308  CALL bcast(alpk)
309  CALL bcast(wdens_ref)
310  CALL bcast(coefgw)
311
312  first=.false.
313 endif
314
315  ! Initialisation de toutes des densites a wdens_ref.
316  ! Les densites peuvent evoluer si les poches debordent
317  ! (voir au tout debut de la boucle sur les substeps)
318  wdens = wdens_ref
319
320  ! print*,'stark',stark
321  ! print*,'alpk',alpk
322  ! print*,'wdens',wdens
323  ! print*,'coefgw',coefgw
324  ! cc
325  ! Minimum value for |T_wake - T_undist|. Used for wake top definition
326  ! -------------------------------------------------------------------------
327
328  delta_t_min = 0.2
329
330  ! 1. - Save initial values and initialize tendencies
331  ! --------------------------------------------------
332
333  DO k = 1, klev
334    DO i = 1, klon
335      ppi(i, k) = pi(i, k)
336      deltatw0(i, k) = deltatw(i, k)
337      deltaqw0(i, k) = deltaqw(i, k)
338      te(i, k) = te0(i, k)
339      qe(i, k) = qe0(i, k)
340      dtls(i, k) = 0.
341      dqls(i, k) = 0.
342      d_deltat_gw(i, k) = 0.
343      d_te(i, k) = 0.
344      d_qe(i, k) = 0.
345      d_deltatw(i, k) = 0.
346      d_deltaqw(i, k) = 0.
347      ! IM 060508 beg
348      d_deltatw2(i, k) = 0.
349      d_deltaqw2(i, k) = 0.
350      ! IM 060508 end
351    END DO
352  END DO
353  ! sigmaw1=sigmaw
354  ! IF (sigd_con.GT.sigmaw1) THEN
355  ! print*, 'sigmaw,sigd_con', sigmaw, sigd_con
356  ! ENDIF
357  DO i = 1, klon
358    ! c      sigmaw(i) = amax1(sigmaw(i),sigd_con(i))
359    sigmaw(i) = amax1(sigmaw(i), sigmad)
360    sigmaw(i) = amin1(sigmaw(i), 0.99)
361    sigmaw0(i) = sigmaw(i)
362    wape(i) = 0.
363    wape2(i) = 0.
364    d_sigmaw(i) = 0.
365    ktopw(i) = 0
366  END DO
367
368
369  ! 2. - Prognostic part
370  ! --------------------
371
372
373  ! 2.1 - Undisturbed area and Wake integrals
374  ! ---------------------------------------------------------
375
376  DO i = 1, klon
377    z(i) = 0.
378    ktop(i) = 0
379    kupper(i) = 0
380    sum_thu(i) = 0.
381    sum_tu(i) = 0.
382    sum_qu(i) = 0.
383    sum_thvu(i) = 0.
384    sum_dth(i) = 0.
385    sum_dq(i) = 0.
386    sum_rho(i) = 0.
387    sum_dtdwn(i) = 0.
388    sum_dqdwn(i) = 0.
389
390    av_thu(i) = 0.
391    av_tu(i) = 0.
392    av_qu(i) = 0.
393    av_thvu(i) = 0.
394    av_dth(i) = 0.
395    av_dq(i) = 0.
396    av_rho(i) = 0.
397    av_dtdwn(i) = 0.
398    av_dqdwn(i) = 0.
399  END DO
400
401  ! Distance between wakes
402  DO i = 1, klon
403    ll(i) = (1-sqrt(sigmaw(i)))/sqrt(wdens(i))
404  END DO
405  ! Potential temperatures and humidity
406  ! ----------------------------------------------------------
407  DO k = 1, klev
408    DO i = 1, klon
409      ! write(*,*)'wake 1',i,k,rd,te(i,k)
410      rho(i, k) = p(i, k)/(rd*te(i,k))
411      ! write(*,*)'wake 2',rho(i,k)
412      IF (k==1) THEN
413        ! write(*,*)'wake 3',i,k,rd,te(i,k)
414        rhoh(i, k) = ph(i, k)/(rd*te(i,k))
415        ! write(*,*)'wake 4',i,k,rd,te(i,k)
416        zhh(i, k) = 0
417      ELSE
418        ! write(*,*)'wake 5',rd,(te(i,k)+te(i,k-1))
419        rhoh(i, k) = ph(i, k)*2./(rd*(te(i,k)+te(i,k-1)))
420        ! write(*,*)'wake 6',(-rhoh(i,k)*RG)+zhh(i,k-1)
421        zhh(i, k) = (ph(i,k)-ph(i,k-1))/(-rhoh(i,k)*rg) + zhh(i, k-1)
422      END IF
423      ! write(*,*)'wake 7',ppi(i,k)
424      the(i, k) = te(i, k)/ppi(i, k)
425      thu(i, k) = (te(i,k)-deltatw(i,k)*sigmaw(i))/ppi(i, k)
426      tu(i, k) = te(i, k) - deltatw(i, k)*sigmaw(i)
427      qu(i, k) = qe(i, k) - deltaqw(i, k)*sigmaw(i)
428      ! write(*,*)'wake 8',(rd*(te(i,k)+deltatw(i,k)))
429      rhow(i, k) = p(i, k)/(rd*(te(i,k)+deltatw(i,k)))
430      dth(i, k) = deltatw(i, k)/ppi(i, k)
431    END DO
432  END DO
433
434  DO k = 1, klev - 1
435    DO i = 1, klon
436      IF (k==1) THEN
437        n2(i, k) = 0
438      ELSE
439        n2(i, k) = amax1(0., -rg**2/the(i,k)*rho(i,k)*(the(i,k+1)-the(i, &
440          k-1))/(p(i,k+1)-p(i,k-1)))
441      END IF
442      zh(i, k) = (zhh(i,k)+zhh(i,k+1))/2
443
444      cgw(i, k) = sqrt(n2(i,k))*zh(i, k)
445      tgw(i, k) = coefgw*cgw(i, k)/ll(i)
446    END DO
447  END DO
448
449  DO i = 1, klon
450    n2(i, klev) = 0
451    zh(i, klev) = 0
452    cgw(i, klev) = 0
453    tgw(i, klev) = 0
454  END DO
455
456  ! Calcul de la masse volumique moyenne de la colonne   (bdlmd)
457  ! -----------------------------------------------------------------
458
459  DO k = 1, klev
460    DO i = 1, klon
461      epaisseur1(i, k) = 0.
462      epaisseur2(i, k) = 0.
463    END DO
464  END DO
465
466  DO i = 1, klon
467    epaisseur1(i, 1) = -(ph(i,2)-ph(i,1))/(rho(i,1)*rg) + 1.
468    epaisseur2(i, 1) = -(ph(i,2)-ph(i,1))/(rho(i,1)*rg) + 1.
469    rhow_moyen(i, 1) = rhow(i, 1)
470  END DO
471
472  DO k = 2, klev
473    DO i = 1, klon
474      epaisseur1(i, k) = -(ph(i,k+1)-ph(i,k))/(rho(i,k)*rg) + 1.
475      epaisseur2(i, k) = epaisseur2(i, k-1) + epaisseur1(i, k)
476      rhow_moyen(i, k) = (rhow_moyen(i,k-1)*epaisseur2(i,k-1)+rhow(i,k)* &
477        epaisseur1(i,k))/epaisseur2(i, k)
478    END DO
479  END DO
480
481
482  ! Choose an integration bound well above wake top
483  ! -----------------------------------------------------------------
484
485  ! Pupper = 50000.  ! melting level
486  ! Pupper = 60000.
487  ! Pupper = 80000.  ! essais pour case_e
488  DO i = 1, klon
489    pupper(i) = 0.6*ph(i, 1)
490    pupper(i) = max(pupper(i), 45000.)
491    ! cc        Pupper(i) = 60000.
492  END DO
493
494
495  ! Determine Wake top pressure (Ptop) from buoyancy integral
496  ! --------------------------------------------------------
497
498  ! -1/ Pressure of the level where dth becomes less than delta_t_min.
499
500  DO i = 1, klon
501    ptop_provis(i) = ph(i, 1)
502  END DO
503  DO k = 2, klev
504    DO i = 1, klon
505
506      ! IM v3JYG; ptop_provis(i).LT. ph(i,1)
507
508      IF (dth(i,k)>-delta_t_min .AND. dth(i,k-1)<-delta_t_min .AND. &
509          ptop_provis(i)==ph(i,1)) THEN
510        ptop_provis(i) = ((dth(i,k)+delta_t_min)*p(i,k-1)-(dth(i, &
511          k-1)+delta_t_min)*p(i,k))/(dth(i,k)-dth(i,k-1))
512      END IF
513    END DO
514  END DO
515
516  ! -2/ dth integral
517
518  DO i = 1, klon
519    sum_dth(i) = 0.
520    dthmin(i) = -delta_t_min
521    z(i) = 0.
522  END DO
523
524  DO k = 1, klev
525    DO i = 1, klon
526      dz(i) = -(amax1(ph(i,k+1),ptop_provis(i))-ph(i,k))/(rho(i,k)*rg)
527      IF (dz(i)>0) THEN
528        z(i) = z(i) + dz(i)
529        sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i)
530        dthmin(i) = amin1(dthmin(i), dth(i,k))
531      END IF
532    END DO
533  END DO
534
535  ! -3/ height of triangle with area= sum_dth and base = dthmin
536
537  DO i = 1, klon
538    hw0(i) = 2.*sum_dth(i)/amin1(dthmin(i), -0.5)
539    hw0(i) = amax1(hwmin, hw0(i))
540  END DO
541
542  ! -4/ now, get Ptop
543
544  DO i = 1, klon
545    z(i) = 0.
546    ptop(i) = ph(i, 1)
547  END DO
548
549  DO k = 1, klev
550    DO i = 1, klon
551      dz(i) = amin1(-(ph(i,k+1)-ph(i,k))/(rho(i,k)*rg), hw0(i)-z(i))
552      IF (dz(i)>0) THEN
553        z(i) = z(i) + dz(i)
554        ptop(i) = ph(i, k) - rho(i, k)*rg*dz(i)
555      END IF
556    END DO
557  END DO
558
559
560  ! -5/ Determination de ktop et kupper
561
562  DO k = klev, 1, -1
563    DO i = 1, klon
564      IF (ph(i,k+1)<ptop(i)) ktop(i) = k
565      IF (ph(i,k+1)<pupper(i)) kupper(i) = k
566    END DO
567  END DO
568
569  ! On evite kupper = 1 et kupper = klev
570  DO i = 1, klon
571    kupper(i) = max(kupper(i), 2)
572    kupper(i) = min(kupper(i), klev-1)
573  END DO
574
575
576  ! -6/ Correct ktop and ptop
577
578  DO i = 1, klon
579    ptop_new(i) = ptop(i)
580  END DO
581  DO k = klev, 2, -1
582    DO i = 1, klon
583      IF (k<=ktop(i) .AND. ptop_new(i)==ptop(i) .AND. &
584          dth(i,k)>-delta_t_min .AND. dth(i,k-1)<-delta_t_min) THEN
585        ptop_new(i) = ((dth(i,k)+delta_t_min)*p(i,k-1)-(dth(i, &
586          k-1)+delta_t_min)*p(i,k))/(dth(i,k)-dth(i,k-1))
587      END IF
588    END DO
589  END DO
590
591  DO i = 1, klon
592    ptop(i) = ptop_new(i)
593  END DO
594
595  DO k = klev, 1, -1
596    DO i = 1, klon
597      IF (ph(i,k+1)<ptop(i)) ktop(i) = k
598    END DO
599  END DO
600
601  ! -5/ Set deltatw & deltaqw to 0 above kupper
602
603  DO k = 1, klev
604    DO i = 1, klon
605      IF (k>=kupper(i)) THEN
606        deltatw(i, k) = 0.
607        deltaqw(i, k) = 0.
608      END IF
609    END DO
610  END DO
611
612
613  ! Vertical gradient of LS omega
614
615  DO k = 1, klev
616    DO i = 1, klon
617      IF (k<=kupper(i)) THEN
618        dp_omgb(i, k) = (omgb(i,k+1)-omgb(i,k))/(ph(i,k+1)-ph(i,k))
619      END IF
620    END DO
621  END DO
622
623  ! Integrals (and wake top level number)
624  ! --------------------------------------
625
626  ! Initialize sum_thvu to 1st level virt. pot. temp.
627
628  DO i = 1, klon
629    z(i) = 1.
630    dz(i) = 1.
631    sum_thvu(i) = thu(i, 1)*(1.+eps*qu(i,1))*dz(i)
632    sum_dth(i) = 0.
633  END DO
634
635  DO k = 1, klev
636    DO i = 1, klon
637      dz(i) = -(amax1(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*rg)
638      IF (dz(i)>0) THEN
639        z(i) = z(i) + dz(i)
640        sum_thu(i) = sum_thu(i) + thu(i, k)*dz(i)
641        sum_tu(i) = sum_tu(i) + tu(i, k)*dz(i)
642        sum_qu(i) = sum_qu(i) + qu(i, k)*dz(i)
643        sum_thvu(i) = sum_thvu(i) + thu(i, k)*(1.+eps*qu(i,k))*dz(i)
644        sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i)
645        sum_dq(i) = sum_dq(i) + deltaqw(i, k)*dz(i)
646        sum_rho(i) = sum_rho(i) + rhow(i, k)*dz(i)
647        sum_dtdwn(i) = sum_dtdwn(i) + dtdwn(i, k)*dz(i)
648        sum_dqdwn(i) = sum_dqdwn(i) + dqdwn(i, k)*dz(i)
649      END IF
650    END DO
651  END DO
652
653  DO i = 1, klon
654    hw0(i) = z(i)
655  END DO
656
657
658  ! 2.1 - WAPE and mean forcing computation
659  ! ---------------------------------------
660
661  ! ---------------------------------------
662
663  ! Means
664
665  DO i = 1, klon
666    av_thu(i) = sum_thu(i)/hw0(i)
667    av_tu(i) = sum_tu(i)/hw0(i)
668    av_qu(i) = sum_qu(i)/hw0(i)
669    av_thvu(i) = sum_thvu(i)/hw0(i)
670    ! av_thve = sum_thve/hw0
671    av_dth(i) = sum_dth(i)/hw0(i)
672    av_dq(i) = sum_dq(i)/hw0(i)
673    av_rho(i) = sum_rho(i)/hw0(i)
674    av_dtdwn(i) = sum_dtdwn(i)/hw0(i)
675    av_dqdwn(i) = sum_dqdwn(i)/hw0(i)
676
677    wape(i) = -rg*hw0(i)*(av_dth(i)+eps*(av_thu(i)*av_dq(i)+av_dth(i)*av_qu(i &
678      )+av_dth(i)*av_dq(i)))/av_thvu(i)
679  END DO
680
681  ! 2.2 Prognostic variable update
682  ! ------------------------------
683
684  ! Filter out bad wakes
685
686  DO k = 1, klev
687    DO i = 1, klon
688      IF (wape(i)<0.) THEN
689        deltatw(i, k) = 0.
690        deltaqw(i, k) = 0.
691        dth(i, k) = 0.
692      END IF
693    END DO
694  END DO
695
696  DO i = 1, klon
697    IF (wape(i)<0.) THEN
698      wape(i) = 0.
699      cstar(i) = 0.
700      hw(i) = hwmin
701      sigmaw(i) = amax1(sigmad, sigd_con(i))
702      fip(i) = 0.
703      gwake(i) = .FALSE.
704    ELSE
705      cstar(i) = stark*sqrt(2.*wape(i))
706      gwake(i) = .TRUE.
707    END IF
708  END DO
709
710
711  ! Check qx and qw positivity
712  ! --------------------------
713  DO i = 1, klon
714    q0_min(i) = min((qe(i,1)-sigmaw(i)*deltaqw(i,1)), (qe(i, &
715      1)+(1.-sigmaw(i))*deltaqw(i,1)))
716  END DO
717  DO k = 2, klev
718    DO i = 1, klon
719      q1_min(i) = min((qe(i,k)-sigmaw(i)*deltaqw(i,k)), (qe(i, &
720        k)+(1.-sigmaw(i))*deltaqw(i,k)))
721      IF (q1_min(i)<=q0_min(i)) THEN
722        q0_min(i) = q1_min(i)
723      END IF
724    END DO
725  END DO
726
727  DO i = 1, klon
728    ok_qx_qw(i) = q0_min(i) >= 0.
729    alpha(i) = 1.
730  END DO
731
732  ! C -----------------------------------------------------------------
733  ! Sub-time-stepping
734  ! -----------------
735
736  nsub = 10
737  dtimesub = dtime/nsub
738
739  ! ------------------------------------------------------------
740  DO isubstep = 1, nsub
741    ! ------------------------------------------------------------
742
743    ! wk_adv is the logical flag enabling wake evolution in the time advance
744    ! loop
745    DO i = 1, klon
746      wk_adv(i) = ok_qx_qw(i) .AND. alpha(i) >= 1.
747    END DO
748
749    ! cc nrlmd   Ajout d'un recalcul de wdens dans le cas d'un entrainement
750    ! négatif de ktop à kupper --------
751    ! cc           On calcule pour cela une densité wdens0 pour laquelle on
752    ! aurait un entrainement nul ---
753    DO i = 1, klon
754      ! c       print *,' isubstep,wk_adv(i),cstar(i),wape(i) ',
755      ! c     $           isubstep,wk_adv(i),cstar(i),wape(i)
756      IF (wk_adv(i) .AND. cstar(i)>0.01) THEN
757        omg(i, kupper(i)+1) = -rg*amdwn(i, kupper(i)+1)/sigmaw(i) + &
758          rg*amup(i, kupper(i)+1)/(1.-sigmaw(i))
759        wdens0 = (sigmaw(i)/(4.*3.14))*((1.-sigmaw(i))*omg(i,kupper(i)+1)/(( &
760          ph(i,1)-pupper(i))*cstar(i)))**(2)
761        IF (wdens(i)<=wdens0*1.1) THEN
762          wdens(i) = wdens0
763        END IF
764        ! c        print*,'omg(i,kupper(i)+1),wdens0,wdens(i),cstar(i)
765        ! c     $     ,ph(i,1)-pupper(i)',
766        ! c     $             omg(i,kupper(i)+1),wdens0,wdens(i),cstar(i)
767        ! c     $     ,ph(i,1)-pupper(i)
768      END IF
769    END DO
770
771    ! cc nrlmd
772
773    DO i = 1, klon
774      IF (wk_adv(i)) THEN
775        gfl(i) = 2.*sqrt(3.14*wdens(i)*sigmaw(i))
776        sigmaw(i) = amin1(sigmaw(i), sigmaw_max)
777      END IF
778    END DO
779    DO i = 1, klon
780      IF (wk_adv(i)) THEN
781        ! cc nrlmd          Introduction du taux de mortalité des poches et
782        ! test sur sigmaw_max=0.4
783        ! cc         d_sigmaw(i) = gfl(i)*Cstar(i)*dtimesub
784        IF (sigmaw(i)>=sigmaw_max) THEN
785          death_rate(i) = gfl(i)*cstar(i)/sigmaw(i)
786        ELSE
787          death_rate(i) = 0.
788        END IF
789        d_sigmaw(i) = gfl(i)*cstar(i)*dtimesub - death_rate(i)*sigmaw(i)* &
790          dtimesub
791        ! $              - nat_rate(i)*sigmaw(i)*dtimesub
792        ! c        print*, 'd_sigmaw(i),sigmaw(i),gfl(i),Cstar(i),wape(i),
793        ! c     $  death_rate(i),ktop(i),kupper(i)',
794        ! c     $                d_sigmaw(i),sigmaw(i),gfl(i),Cstar(i),wape(i),
795        ! c     $  death_rate(i),ktop(i),kupper(i)
796
797        ! sigmaw(i) =sigmaw(i) + gfl(i)*Cstar(i)*dtimesub
798        ! sigmaw(i) =min(sigmaw(i),0.99)     !!!!!!!!
799        ! wdens = wdens0/(10.*sigmaw)
800        ! sigmaw =max(sigmaw,sigd_con)
801        ! sigmaw =max(sigmaw,sigmad)
802      END IF
803    END DO
804
805
806    ! calcul de la difference de vitesse verticale poche - zone non perturbee
807    ! IM 060208 differences par rapport au code initial; init. a 0 dp_deltomg
808    ! IM 060208 et omg sur les niveaux de 1 a klev+1, alors que avant l'on
809    ! definit
810    ! IM 060208 au niveau k=1..?
811    DO k = 1, klev
812      DO i = 1, klon
813        IF (wk_adv(i)) THEN !!! nrlmd
814          dp_deltomg(i, k) = 0.
815        END IF
816      END DO
817    END DO
818    DO k = 1, klev + 1
819      DO i = 1, klon
820        IF (wk_adv(i)) THEN !!! nrlmd
821          omg(i, k) = 0.
822        END IF
823      END DO
824    END DO
825
826    DO i = 1, klon
827      IF (wk_adv(i)) THEN
828        z(i) = 0.
829        omg(i, 1) = 0.
830        dp_deltomg(i, 1) = -(gfl(i)*cstar(i))/(sigmaw(i)*(1-sigmaw(i)))
831      END IF
832    END DO
833
834    DO k = 2, klev
835      DO i = 1, klon
836        IF (wk_adv(i) .AND. k<=ktop(i)) THEN
837          dz(i) = -(ph(i,k)-ph(i,k-1))/(rho(i,k-1)*rg)
838          z(i) = z(i) + dz(i)
839          dp_deltomg(i, k) = dp_deltomg(i, 1)
840          omg(i, k) = dp_deltomg(i, 1)*z(i)
841        END IF
842      END DO
843    END DO
844
845    DO i = 1, klon
846      IF (wk_adv(i)) THEN
847        dztop(i) = -(ptop(i)-ph(i,ktop(i)))/(rho(i,ktop(i))*rg)
848        ztop(i) = z(i) + dztop(i)
849        omgtop(i) = dp_deltomg(i, 1)*ztop(i)
850      END IF
851    END DO
852
853    ! -----------------
854    ! From m/s to Pa/s
855    ! -----------------
856
857    DO i = 1, klon
858      IF (wk_adv(i)) THEN
859        omgtop(i) = -rho(i, ktop(i))*rg*omgtop(i)
860        dp_deltomg(i, 1) = omgtop(i)/(ptop(i)-ph(i,1))
861      END IF
862    END DO
863
864    DO k = 1, klev
865      DO i = 1, klon
866        IF (wk_adv(i) .AND. k<=ktop(i)) THEN
867          omg(i, k) = -rho(i, k)*rg*omg(i, k)
868          dp_deltomg(i, k) = dp_deltomg(i, 1)
869        END IF
870      END DO
871    END DO
872
873    ! raccordement lineaire de omg de ptop a pupper
874
875    DO i = 1, klon
876      IF (wk_adv(i) .AND. kupper(i)>ktop(i)) THEN
877        omg(i, kupper(i)+1) = -rg*amdwn(i, kupper(i)+1)/sigmaw(i) + &
878          rg*amup(i, kupper(i)+1)/(1.-sigmaw(i))
879        dp_deltomg(i, kupper(i)) = (omgtop(i)-omg(i,kupper(i)+1))/ &
880          (ptop(i)-pupper(i))
881      END IF
882    END DO
883
884    ! c      DO i=1,klon
885    ! c        print*,'Pente entre 0 et kupper (référence)'
886    ! c     $           ,omg(i,kupper(i)+1)/(pupper(i)-ph(i,1))
887    ! c        print*,'Pente entre ktop et kupper'
888    ! c     $   ,(omg(i,kupper(i)+1)-omgtop(i))/(pupper(i)-ptop(i))
889    ! c      ENDDO
890    ! c
891    DO k = 1, klev
892      DO i = 1, klon
893        IF (wk_adv(i) .AND. k>ktop(i) .AND. k<=kupper(i)) THEN
894          dp_deltomg(i, k) = dp_deltomg(i, kupper(i))
895          omg(i, k) = omgtop(i) + (ph(i,k)-ptop(i))*dp_deltomg(i, kupper(i))
896        END IF
897      END DO
898    END DO
899    ! cc nrlmd
900    ! c      DO i=1,klon
901    ! c      print*,'deltaw_ktop,deltaw_conv',omgtop(i),omg(i,kupper(i)+1)
902    ! c      END DO
903    ! cc
904
905
906    ! --    Compute wake average vertical velocity omgbw
907
908
909    DO k = 1, klev + 1
910      DO i = 1, klon
911        IF (wk_adv(i)) THEN
912          omgbw(i, k) = omgb(i, k) + (1.-sigmaw(i))*omg(i, k)
913        END IF
914      END DO
915    END DO
916    ! --    and its vertical gradient dp_omgbw
917
918    DO k = 1, klev
919      DO i = 1, klon
920        IF (wk_adv(i)) THEN
921          dp_omgbw(i, k) = (omgbw(i,k+1)-omgbw(i,k))/(ph(i,k+1)-ph(i,k))
922        END IF
923      END DO
924    END DO
925
926    ! --    Upstream coefficients for omgb velocity
927    ! --    (alpha_up(k) is the coefficient of the value at level k)
928    ! --    (1-alpha_up(k) is the coefficient of the value at level k-1)
929    DO k = 1, klev
930      DO i = 1, klon
931        IF (wk_adv(i)) THEN
932          alpha_up(i, k) = 0.
933          IF (omgb(i,k)>0.) alpha_up(i, k) = 1.
934        END IF
935      END DO
936    END DO
937
938    ! Matrix expressing [The,deltatw] from  [Th1,Th2]
939
940    DO i = 1, klon
941      IF (wk_adv(i)) THEN
942        rre1(i) = 1. - sigmaw(i)
943        rre2(i) = sigmaw(i)
944      END IF
945    END DO
946    rrd1 = -1.
947    rrd2 = 1.
948
949    ! --    Get [Th1,Th2], dth and [q1,q2]
950
951    DO k = 1, klev
952      DO i = 1, klon
953        IF (wk_adv(i) .AND. k<=kupper(i)+1) THEN
954          dth(i, k) = deltatw(i, k)/ppi(i, k)
955          th1(i, k) = the(i, k) - sigmaw(i)*dth(i, k) ! undisturbed area
956          th2(i, k) = the(i, k) + (1.-sigmaw(i))*dth(i, k) ! wake
957          q1(i, k) = qe(i, k) - sigmaw(i)*deltaqw(i, k) ! undisturbed area
958          q2(i, k) = qe(i, k) + (1.-sigmaw(i))*deltaqw(i, k) ! wake
959        END IF
960      END DO
961    END DO
962
963    DO i = 1, klon
964      IF (wk_adv(i)) THEN !!! nrlmd
965        d_th1(i, 1) = 0.
966        d_th2(i, 1) = 0.
967        d_dth(i, 1) = 0.
968        d_q1(i, 1) = 0.
969        d_q2(i, 1) = 0.
970        d_dq(i, 1) = 0.
971      END IF
972    END DO
973
974    DO k = 2, klev
975      DO i = 1, klon
976        IF (wk_adv(i) .AND. k<=kupper(i)+1) THEN
977          d_th1(i, k) = th1(i, k-1) - th1(i, k)
978          d_th2(i, k) = th2(i, k-1) - th2(i, k)
979          d_dth(i, k) = dth(i, k-1) - dth(i, k)
980          d_q1(i, k) = q1(i, k-1) - q1(i, k)
981          d_q2(i, k) = q2(i, k-1) - q2(i, k)
982          d_dq(i, k) = deltaqw(i, k-1) - deltaqw(i, k)
983        END IF
984      END DO
985    END DO
986
987    DO i = 1, klon
988      IF (wk_adv(i)) THEN
989        omgbdth(i, 1) = 0.
990        omgbdq(i, 1) = 0.
991      END IF
992    END DO
993
994    DO k = 2, klev
995      DO i = 1, klon
996        IF (wk_adv(i) .AND. k<=kupper(i)+1) THEN !   loop on interfaces
997          omgbdth(i, k) = omgb(i, k)*(dth(i,k-1)-dth(i,k))
998          omgbdq(i, k) = omgb(i, k)*(deltaqw(i,k-1)-deltaqw(i,k))
999        END IF
1000      END DO
1001    END DO
1002
1003    ! -----------------------------------------------------------------
1004    DO k = 1, klev
1005      DO i = 1, klon
1006        IF (wk_adv(i) .AND. k<=kupper(i)-1) THEN
1007          ! -----------------------------------------------------------------
1008
1009          ! Compute redistribution (advective) term
1010
1011          d_deltatw(i, k) = dtimesub/(ph(i,k)-ph(i,k+1))* &
1012            (rrd1*omg(i,k)*sigmaw(i)*d_th1(i,k)-rrd2*omg(i,k+1)*(1.-sigmaw( &
1013            i))*d_th2(i,k+1)-(1.-alpha_up(i,k))*omgbdth(i,k)-alpha_up(i,k+1)* &
1014            omgbdth(i,k+1))*ppi(i, k)
1015          ! print*,'d_deltatw=',d_deltatw(i,k)
1016
1017          d_deltaqw(i, k) = dtimesub/(ph(i,k)-ph(i,k+1))* &
1018            (rrd1*omg(i,k)*sigmaw(i)*d_q1(i,k)-rrd2*omg(i,k+1)*(1.-sigmaw( &
1019            i))*d_q2(i,k+1)-(1.-alpha_up(i,k))*omgbdq(i,k)-alpha_up(i,k+1)* &
1020            omgbdq(i,k+1))
1021          ! print*,'d_deltaqw=',d_deltaqw(i,k)
1022
1023          ! and increment large scale tendencies
1024
1025
1026
1027
1028          ! C
1029          ! -----------------------------------------------------------------
1030          d_te(i, k) = dtimesub*((rre1(i)*omg(i,k)*sigmaw(i)*d_th1(i, &
1031            k)-rre2(i)*omg(i,k+1)*(1.-sigmaw(i))*d_th2(i,k+1))/(ph(i,k)-ph(i, &
1032            k+1)) &                ! cc nrlmd     $
1033                                   ! -sigmaw(i)*(1.-sigmaw(i))*dth(i,k)*dp_deltomg(i,k)
1034            -sigmaw(i)*(1.-sigmaw(i))*dth(i,k)*(omg(i,k)-omg(i,k+1))/(ph(i, &
1035            k)-ph(i,k+1)) &        ! cc
1036            )*ppi(i, k)
1037
1038          d_qe(i, k) = dtimesub*((rre1(i)*omg(i,k)*sigmaw(i)*d_q1(i, &
1039            k)-rre2(i)*omg(i,k+1)*(1.-sigmaw(i))*d_q2(i,k+1))/(ph(i,k)-ph(i, &
1040            k+1)) &                ! cc nrlmd     $
1041                                   ! -sigmaw(i)*(1.-sigmaw(i))*deltaqw(i,k)*dp_deltomg(i,k)
1042            -sigmaw(i)*(1.-sigmaw(i))*deltaqw(i,k)*(omg(i,k)-omg(i, &
1043            k+1))/(ph(i,k)-ph(i,k+1)) & ! cc
1044            )
1045          ! cc nrlmd
1046        ELSE IF (wk_adv(i) .AND. k==kupper(i)) THEN
1047          d_te(i, k) = dtimesub*((rre1(i)*omg(i,k)*sigmaw(i)*d_th1(i, &
1048            k)/(ph(i,k)-ph(i,k+1))))*ppi(i, k)
1049
1050          d_qe(i, k) = dtimesub*((rre1(i)*omg(i,k)*sigmaw(i)*d_q1(i, &
1051            k)/(ph(i,k)-ph(i,k+1))))
1052
1053        END IF
1054        ! cc
1055      END DO
1056    END DO
1057    ! ------------------------------------------------------------------
1058
1059    ! Increment state variables
1060
1061    DO k = 1, klev
1062      DO i = 1, klon
1063        ! cc nrlmd       IF( wk_adv(i) .AND. k .LE. kupper(i)-1) THEN
1064        IF (wk_adv(i) .AND. k<=kupper(i)) THEN
1065          ! cc
1066
1067
1068
1069          ! Coefficient de répartition
1070
1071          crep(i, k) = crep_sol*(ph(i,kupper(i))-ph(i,k))/ &
1072            (ph(i,kupper(i))-ph(i,1))
1073          crep(i, k) = crep(i, k) + crep_upper*(ph(i,1)-ph(i,k))/(p(i,1)-ph(i &
1074            ,kupper(i)))
1075
1076
1077          ! Reintroduce compensating subsidence term.
1078
1079          ! dtKE(k)=(dtdwn(k)*Crep(k))/sigmaw
1080          ! dtKE(k)=dtKE(k)-(dtdwn(k)*(1-Crep(k))+dta(k))
1081          ! .                   /(1-sigmaw)
1082          ! dqKE(k)=(dqdwn(k)*Crep(k))/sigmaw
1083          ! dqKE(k)=dqKE(k)-(dqdwn(k)*(1-Crep(k))+dqa(k))
1084          ! .                   /(1-sigmaw)
1085
1086          ! dtKE(k)=(dtdwn(k)*Crep(k)+(1-Crep(k))*dta(k))/sigmaw
1087          ! dtKE(k)=dtKE(k)-(dtdwn(k)*(1-Crep(k))+dta(k)*Crep(k))
1088          ! .                   /(1-sigmaw)
1089          ! dqKE(k)=(dqdwn(k)*Crep(k)+(1-Crep(k))*dqa(k))/sigmaw
1090          ! dqKE(k)=dqKE(k)-(dqdwn(k)*(1-Crep(k))+dqa(k)*Crep(k))
1091          ! .                   /(1-sigmaw)
1092
1093          dtke(i, k) = (dtdwn(i,k)/sigmaw(i)-dta(i,k)/(1.-sigmaw(i)))
1094          dqke(i, k) = (dqdwn(i,k)/sigmaw(i)-dqa(i,k)/(1.-sigmaw(i)))
1095          ! print*,'dtKE= ',dtKE(i,k),' dqKE= ',dqKE(i,k)
1096
1097!jyg<
1098!!
1099!!---------------------------------------------------------------
1100!! The change of delta_T due to PBL (vertical diffusion plus thermal plumes)
1101!! is accounted for by the PBL and the Thermals schemes. It is now set to zero
1102!! within the Wake scheme.
1103!!---------------------------------------------------------------
1104          dtPBL(i,k) = 0.
1105          dqPBL(i,k) = 0.
1106!
1107!!          dtPBL(i,k)=wdtPBL(i,k) - udtPBL(i,k)
1108!!          dqPBL(i,k)=wdqPBL(i,k) - udqPBL(i,k)
1109!
1110!!          dtpbl(i, k) = (wdtpbl(i,k)/sigmaw(i)-udtpbl(i,k)/(1.-sigmaw(i)))
1111!!          dqpbl(i, k) = (wdqpbl(i,k)/sigmaw(i)-udqpbl(i,k)/(1.-sigmaw(i)))
1112          ! print*,'dtPBL= ',dtPBL(i,k),' dqPBL= ',dqPBL(i,k)
1113!>jyg
1114!
1115
1116          ! cc nrlmd          Prise en compte du taux de mortalité
1117          ! cc               Définitions de entr, detr
1118          detr(i, k) = 0.
1119
1120          entr(i, k) = detr(i, k) + gfl(i)*cstar(i) + &
1121            sigmaw(i)*(1.-sigmaw(i))*dp_deltomg(i, k)
1122
1123          spread(i, k) = (entr(i,k)-detr(i,k))/sigmaw(i)
1124          ! cc        spread(i,k) =
1125          ! (1.-sigmaw(i))*dp_deltomg(i,k)+gfl(i)*Cstar(i)/
1126          ! cc     $  sigmaw(i)
1127
1128
1129          ! ajout d'un effet onde de gravité -Tgw(k)*deltatw(k) 03/02/06 YU
1130          ! Jingmei
1131
1132          ! write(lunout,*)'wake.F ',i,k, dtimesub,d_deltat_gw(i,k),
1133          ! &  Tgw(i,k),deltatw(i,k)
1134          d_deltat_gw(i, k) = d_deltat_gw(i, k) - tgw(i, k)*deltatw(i, k)* &
1135            dtimesub
1136          ! write(lunout,*)'wake.F ',i,k, dtimesub,d_deltatw(i,k)
1137          ff(i) = d_deltatw(i, k)/dtimesub
1138
1139          ! Sans GW
1140
1141          ! deltatw(k)=deltatw(k)+dtimesub*(ff+dtKE(k)-spread(k)*deltatw(k))
1142
1143          ! GW formule 1
1144
1145          ! deltatw(k) = deltatw(k)+dtimesub*
1146          ! $         (ff+dtKE(k) - spread(k)*deltatw(k)-Tgw(k)*deltatw(k))
1147
1148          ! GW formule 2
1149
1150          IF (dtimesub*tgw(i,k)<1.E-10) THEN
1151            d_deltatw(i, k) = dtimesub*(ff(i)+dtke(i,k)+dtpbl(i,k) & ! cc
1152                                                                     ! $
1153                                                                     ! -spread(i,k)*deltatw(i,k)
1154              -entr(i,k)*deltatw(i,k)/sigmaw(i)-(death_rate(i)*sigmaw( &
1155              i)+detr(i,k))*deltatw(i,k)/(1.-sigmaw(i)) & ! cc
1156              -tgw(i,k)*deltatw(i,k))
1157          ELSE
1158            d_deltatw(i, k) = 1/tgw(i, k)*(1-exp(-dtimesub*tgw(i, &
1159              k)))*(ff(i)+dtke(i,k)+dtpbl(i,k) & ! cc     $
1160                                                 ! -spread(i,k)*deltatw(i,k)
1161              -entr(i,k)*deltatw(i,k)/sigmaw(i)-(death_rate(i)*sigmaw( &
1162              i)+detr(i,k))*deltatw(i,k)/(1.-sigmaw(i)) & ! cc
1163              -tgw(i,k)*deltatw(i,k))
1164          END IF
1165
1166          dth(i, k) = deltatw(i, k)/ppi(i, k)
1167
1168          gg(i) = d_deltaqw(i, k)/dtimesub
1169
1170          d_deltaqw(i, k) = dtimesub*(gg(i)+dqke(i,k)+dqpbl(i,k) & ! cc     $
1171                                                                   ! -spread(i,k)*deltaqw(i,k))
1172            -entr(i,k)*deltaqw(i,k)/sigmaw(i)-(death_rate(i)*sigmaw(i)+detr( &
1173            i,k))*deltaqw(i,k)/(1.-sigmaw(i)))
1174          ! cc
1175
1176          ! cc nrlmd
1177          ! cc       d_deltatw2(i,k)=d_deltatw2(i,k)+d_deltatw(i,k)
1178          ! cc       d_deltaqw2(i,k)=d_deltaqw2(i,k)+d_deltaqw(i,k)
1179          ! cc
1180        END IF
1181      END DO
1182    END DO
1183
1184
1185    ! Scale tendencies so that water vapour remains positive in w and x.
1186
1187    CALL wake_vec_modulation(klon, klev, wk_adv, epsilon, qe, d_qe, deltaqw, &
1188      d_deltaqw, sigmaw, d_sigmaw, alpha)
1189
1190    ! cc nrlmd
1191    ! c      print*,'alpha'
1192    ! c      do i=1,klon
1193    ! c         print*,alpha(i)
1194    ! c      end do
1195    ! cc
1196    DO k = 1, klev
1197      DO i = 1, klon
1198        IF (wk_adv(i) .AND. k<=kupper(i)) THEN
1199          d_te(i, k) = alpha(i)*d_te(i, k)
1200          d_qe(i, k) = alpha(i)*d_qe(i, k)
1201          d_deltatw(i, k) = alpha(i)*d_deltatw(i, k)
1202          d_deltaqw(i, k) = alpha(i)*d_deltaqw(i, k)
1203          d_deltat_gw(i, k) = alpha(i)*d_deltat_gw(i, k)
1204        END IF
1205      END DO
1206    END DO
1207    DO i = 1, klon
1208      IF (wk_adv(i)) THEN
1209        d_sigmaw(i) = alpha(i)*d_sigmaw(i)
1210      END IF
1211    END DO
1212
1213    ! Update large scale variables and wake variables
1214    ! IM 060208 manque DO i + remplace DO k=1,kupper(i)
1215    ! IM 060208     DO k = 1,kupper(i)
1216    DO k = 1, klev
1217      DO i = 1, klon
1218        IF (wk_adv(i) .AND. k<=kupper(i)) THEN
1219          dtls(i, k) = dtls(i, k) + d_te(i, k)
1220          dqls(i, k) = dqls(i, k) + d_qe(i, k)
1221          ! cc nrlmd
1222          d_deltatw2(i, k) = d_deltatw2(i, k) + d_deltatw(i, k)
1223          d_deltaqw2(i, k) = d_deltaqw2(i, k) + d_deltaqw(i, k)
1224          ! cc
1225        END IF
1226      END DO
1227    END DO
1228    DO k = 1, klev
1229      DO i = 1, klon
1230        IF (wk_adv(i) .AND. k<=kupper(i)) THEN
1231          te(i, k) = te0(i, k) + dtls(i, k)
1232          qe(i, k) = qe0(i, k) + dqls(i, k)
1233          the(i, k) = te(i, k)/ppi(i, k)
1234          deltatw(i, k) = deltatw(i, k) + d_deltatw(i, k)
1235          deltaqw(i, k) = deltaqw(i, k) + d_deltaqw(i, k)
1236          dth(i, k) = deltatw(i, k)/ppi(i, k)
1237          ! c      print*,'k,qx,qw',k,qe(i,k)-sigmaw(i)*deltaqw(i,k)
1238          ! c     $        ,qe(i,k)+(1-sigmaw(i))*deltaqw(i,k)
1239        END IF
1240      END DO
1241    END DO
1242    DO i = 1, klon
1243      IF (wk_adv(i)) THEN
1244        sigmaw(i) = sigmaw(i) + d_sigmaw(i)
1245      END IF
1246    END DO
1247
1248
1249    ! Determine Ptop from buoyancy integral
1250    ! ---------------------------------------
1251
1252    ! -     1/ Pressure of the level where dth changes sign.
1253
1254    DO i = 1, klon
1255      IF (wk_adv(i)) THEN
1256        ptop_provis(i) = ph(i, 1)
1257      END IF
1258    END DO
1259
1260    DO k = 2, klev
1261      DO i = 1, klon
1262        IF (wk_adv(i) .AND. ptop_provis(i)==ph(i,1) .AND. &
1263            dth(i,k)>-delta_t_min .AND. dth(i,k-1)<-delta_t_min) THEN
1264          ptop_provis(i) = ((dth(i,k)+delta_t_min)*p(i,k-1)-(dth(i, &
1265            k-1)+delta_t_min)*p(i,k))/(dth(i,k)-dth(i,k-1))
1266        END IF
1267      END DO
1268    END DO
1269
1270    ! -     2/ dth integral
1271
1272    DO i = 1, klon
1273      IF (wk_adv(i)) THEN !!! nrlmd
1274        sum_dth(i) = 0.
1275        dthmin(i) = -delta_t_min
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) = -(amax1(ph(i,k+1),ptop_provis(i))-ph(i,k))/(rho(i,k)*rg)
1284          IF (dz(i)>0) THEN
1285            z(i) = z(i) + dz(i)
1286            sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i)
1287            dthmin(i) = amin1(dthmin(i), dth(i,k))
1288          END IF
1289        END IF
1290      END DO
1291    END DO
1292
1293    ! -     3/ height of triangle with area= sum_dth and base = dthmin
1294
1295    DO i = 1, klon
1296      IF (wk_adv(i)) THEN
1297        hw(i) = 2.*sum_dth(i)/amin1(dthmin(i), -0.5)
1298        hw(i) = amax1(hwmin, hw(i))
1299      END IF
1300    END DO
1301
1302    ! -     4/ now, get Ptop
1303
1304    DO i = 1, klon
1305      IF (wk_adv(i)) THEN !!! nrlmd
1306        ktop(i) = 0
1307        z(i) = 0.
1308      END IF
1309    END DO
1310
1311    DO k = 1, klev
1312      DO i = 1, klon
1313        IF (wk_adv(i)) THEN
1314          dz(i) = amin1(-(ph(i,k+1)-ph(i,k))/(rho(i,k)*rg), hw(i)-z(i))
1315          IF (dz(i)>0) THEN
1316            z(i) = z(i) + dz(i)
1317            ptop(i) = ph(i, k) - rho(i, k)*rg*dz(i)
1318            ktop(i) = k
1319          END IF
1320        END IF
1321      END DO
1322    END DO
1323
1324    ! 4.5/Correct ktop and ptop
1325
1326    DO i = 1, klon
1327      IF (wk_adv(i)) THEN
1328        ptop_new(i) = ptop(i)
1329      END IF
1330    END DO
1331
1332    DO k = klev, 2, -1
1333      DO i = 1, klon
1334        ! IM v3JYG; IF (k .GE. ktop(i)
1335        IF (wk_adv(i) .AND. k<=ktop(i) .AND. ptop_new(i)==ptop(i) .AND. &
1336            dth(i,k)>-delta_t_min .AND. dth(i,k-1)<-delta_t_min) THEN
1337          ptop_new(i) = ((dth(i,k)+delta_t_min)*p(i,k-1)-(dth(i, &
1338            k-1)+delta_t_min)*p(i,k))/(dth(i,k)-dth(i,k-1))
1339        END IF
1340      END DO
1341    END DO
1342
1343
1344    DO i = 1, klon
1345      IF (wk_adv(i)) THEN
1346        ptop(i) = ptop_new(i)
1347      END IF
1348    END DO
1349
1350    DO k = klev, 1, -1
1351      DO i = 1, klon
1352        IF (wk_adv(i)) THEN !!! nrlmd
1353          IF (ph(i,k+1)<ptop(i)) ktop(i) = k
1354        END IF
1355      END DO
1356    END DO
1357
1358    ! 5/ Set deltatw & deltaqw to 0 above kupper
1359
1360    DO k = 1, klev
1361      DO i = 1, klon
1362        IF (wk_adv(i) .AND. k>=kupper(i)) THEN
1363          deltatw(i, k) = 0.
1364          deltaqw(i, k) = 0.
1365        END IF
1366      END DO
1367    END DO
1368
1369
1370    ! -------------Cstar computation---------------------------------
1371    DO i = 1, klon
1372      IF (wk_adv(i)) THEN !!! nrlmd
1373        sum_thu(i) = 0.
1374        sum_tu(i) = 0.
1375        sum_qu(i) = 0.
1376        sum_thvu(i) = 0.
1377        sum_dth(i) = 0.
1378        sum_dq(i) = 0.
1379        sum_rho(i) = 0.
1380        sum_dtdwn(i) = 0.
1381        sum_dqdwn(i) = 0.
1382
1383        av_thu(i) = 0.
1384        av_tu(i) = 0.
1385        av_qu(i) = 0.
1386        av_thvu(i) = 0.
1387        av_dth(i) = 0.
1388        av_dq(i) = 0.
1389        av_rho(i) = 0.
1390        av_dtdwn(i) = 0.
1391        av_dqdwn(i) = 0.
1392      END IF
1393    END DO
1394
1395    ! Integrals (and wake top level number)
1396    ! --------------------------------------
1397
1398    ! Initialize sum_thvu to 1st level virt. pot. temp.
1399
1400    DO i = 1, klon
1401      IF (wk_adv(i)) THEN !!! nrlmd
1402        z(i) = 1.
1403        dz(i) = 1.
1404        sum_thvu(i) = thu(i, 1)*(1.+eps*qu(i,1))*dz(i)
1405        sum_dth(i) = 0.
1406      END IF
1407    END DO
1408
1409    DO k = 1, klev
1410      DO i = 1, klon
1411        IF (wk_adv(i)) THEN !!! nrlmd
1412          dz(i) = -(max(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*rg)
1413          IF (dz(i)>0) THEN
1414            z(i) = z(i) + dz(i)
1415            sum_thu(i) = sum_thu(i) + thu(i, k)*dz(i)
1416            sum_tu(i) = sum_tu(i) + tu(i, k)*dz(i)
1417            sum_qu(i) = sum_qu(i) + qu(i, k)*dz(i)
1418            sum_thvu(i) = sum_thvu(i) + thu(i, k)*(1.+eps*qu(i,k))*dz(i)
1419            sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i)
1420            sum_dq(i) = sum_dq(i) + deltaqw(i, k)*dz(i)
1421            sum_rho(i) = sum_rho(i) + rhow(i, k)*dz(i)
1422            sum_dtdwn(i) = sum_dtdwn(i) + dtdwn(i, k)*dz(i)
1423            sum_dqdwn(i) = sum_dqdwn(i) + dqdwn(i, k)*dz(i)
1424          END IF
1425        END IF
1426      END DO
1427    END DO
1428
1429    DO i = 1, klon
1430      IF (wk_adv(i)) THEN !!! nrlmd
1431        hw0(i) = z(i)
1432      END IF
1433    END DO
1434
1435
1436    ! - WAPE and mean forcing computation
1437    ! ---------------------------------------
1438
1439    ! ---------------------------------------
1440
1441    ! Means
1442
1443    DO i = 1, klon
1444      IF (wk_adv(i)) THEN !!! nrlmd
1445        av_thu(i) = sum_thu(i)/hw0(i)
1446        av_tu(i) = sum_tu(i)/hw0(i)
1447        av_qu(i) = sum_qu(i)/hw0(i)
1448        av_thvu(i) = sum_thvu(i)/hw0(i)
1449        av_dth(i) = sum_dth(i)/hw0(i)
1450        av_dq(i) = sum_dq(i)/hw0(i)
1451        av_rho(i) = sum_rho(i)/hw0(i)
1452        av_dtdwn(i) = sum_dtdwn(i)/hw0(i)
1453        av_dqdwn(i) = sum_dqdwn(i)/hw0(i)
1454
1455        wape(i) = -rg*hw0(i)*(av_dth(i)+eps*(av_thu(i)*av_dq(i)+av_dth(i)* &
1456          av_qu(i)+av_dth(i)*av_dq(i)))/av_thvu(i)
1457      END IF
1458    END DO
1459
1460    ! Filter out bad wakes
1461
1462    DO k = 1, klev
1463      DO i = 1, klon
1464        IF (wk_adv(i)) THEN !!! nrlmd
1465          IF (wape(i)<0.) THEN
1466            deltatw(i, k) = 0.
1467            deltaqw(i, k) = 0.
1468            dth(i, k) = 0.
1469          END IF
1470        END IF
1471      END DO
1472    END DO
1473
1474    DO i = 1, klon
1475      IF (wk_adv(i)) THEN !!! nrlmd
1476        IF (wape(i)<0.) THEN
1477          wape(i) = 0.
1478          cstar(i) = 0.
1479          hw(i) = hwmin
1480          sigmaw(i) = max(sigmad, sigd_con(i))
1481          fip(i) = 0.
1482          gwake(i) = .FALSE.
1483        ELSE
1484          cstar(i) = stark*sqrt(2.*wape(i))
1485          gwake(i) = .TRUE.
1486        END IF
1487      END IF
1488    END DO
1489
1490  END DO ! end sub-timestep loop
1491
1492  ! -----------------------------------------------------------------
1493  ! Get back to tendencies per second
1494
1495  DO k = 1, klev
1496    DO i = 1, klon
1497
1498      ! cc nrlmd        IF ( wk_adv(i) .AND. k .LE. kupper(i)) THEN
1499      IF (ok_qx_qw(i) .AND. k<=kupper(i)) THEN
1500        ! cc
1501        dtls(i, k) = dtls(i, k)/dtime
1502        dqls(i, k) = dqls(i, k)/dtime
1503        d_deltatw2(i, k) = d_deltatw2(i, k)/dtime
1504        d_deltaqw2(i, k) = d_deltaqw2(i, k)/dtime
1505        d_deltat_gw(i, k) = d_deltat_gw(i, k)/dtime
1506        ! c      print*,'k,dqls,omg,entr,detr',k,dqls(i,k),omg(i,k),entr(i,k)
1507        ! c     $         ,death_rate(i)*sigmaw(i)
1508      END IF
1509    END DO
1510  END DO
1511
1512
1513  ! ----------------------------------------------------------
1514  ! Determine wake final state; recompute wape, cstar, ktop;
1515  ! filter out bad wakes.
1516  ! ----------------------------------------------------------
1517
1518  ! 2.1 - Undisturbed area and Wake integrals
1519  ! ---------------------------------------------------------
1520
1521  DO i = 1, klon
1522    ! cc nrlmd       if (wk_adv(i)) then !!! nrlmd
1523    IF (ok_qx_qw(i)) THEN
1524      ! cc
1525      z(i) = 0.
1526      sum_thu(i) = 0.
1527      sum_tu(i) = 0.
1528      sum_qu(i) = 0.
1529      sum_thvu(i) = 0.
1530      sum_dth(i) = 0.
1531      sum_dq(i) = 0.
1532      sum_rho(i) = 0.
1533      sum_dtdwn(i) = 0.
1534      sum_dqdwn(i) = 0.
1535
1536      av_thu(i) = 0.
1537      av_tu(i) = 0.
1538      av_qu(i) = 0.
1539      av_thvu(i) = 0.
1540      av_dth(i) = 0.
1541      av_dq(i) = 0.
1542      av_rho(i) = 0.
1543      av_dtdwn(i) = 0.
1544      av_dqdwn(i) = 0.
1545    END IF
1546  END DO
1547  ! Potential temperatures and humidity
1548  ! ----------------------------------------------------------
1549
1550  DO k = 1, klev
1551    DO i = 1, klon
1552      ! cc nrlmd       IF ( wk_adv(i)) THEN
1553      IF (ok_qx_qw(i)) THEN
1554        ! cc
1555        rho(i, k) = p(i, k)/(rd*te(i,k))
1556        IF (k==1) THEN
1557          rhoh(i, k) = ph(i, k)/(rd*te(i,k))
1558          zhh(i, k) = 0
1559        ELSE
1560          rhoh(i, k) = ph(i, k)*2./(rd*(te(i,k)+te(i,k-1)))
1561          zhh(i, k) = (ph(i,k)-ph(i,k-1))/(-rhoh(i,k)*rg) + zhh(i, k-1)
1562        END IF
1563        the(i, k) = te(i, k)/ppi(i, k)
1564        thu(i, k) = (te(i,k)-deltatw(i,k)*sigmaw(i))/ppi(i, k)
1565        tu(i, k) = te(i, k) - deltatw(i, k)*sigmaw(i)
1566        qu(i, k) = qe(i, k) - deltaqw(i, k)*sigmaw(i)
1567        rhow(i, k) = p(i, k)/(rd*(te(i,k)+deltatw(i,k)))
1568        dth(i, k) = deltatw(i, k)/ppi(i, k)
1569      END IF
1570    END DO
1571  END DO
1572
1573  ! Integrals (and wake top level number)
1574  ! -----------------------------------------------------------
1575
1576  ! Initialize sum_thvu to 1st level virt. pot. temp.
1577
1578  DO i = 1, klon
1579    ! cc nrlmd       IF ( wk_adv(i)) THEN
1580    IF (ok_qx_qw(i)) THEN
1581      ! cc
1582      z(i) = 1.
1583      dz(i) = 1.
1584      sum_thvu(i) = thu(i, 1)*(1.+eps*qu(i,1))*dz(i)
1585      sum_dth(i) = 0.
1586    END IF
1587  END DO
1588
1589  DO k = 1, klev
1590    DO i = 1, klon
1591      ! cc nrlmd       IF ( wk_adv(i)) THEN
1592      IF (ok_qx_qw(i)) THEN
1593        ! cc
1594        dz(i) = -(amax1(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*rg)
1595        IF (dz(i)>0) THEN
1596          z(i) = z(i) + dz(i)
1597          sum_thu(i) = sum_thu(i) + thu(i, k)*dz(i)
1598          sum_tu(i) = sum_tu(i) + tu(i, k)*dz(i)
1599          sum_qu(i) = sum_qu(i) + qu(i, k)*dz(i)
1600          sum_thvu(i) = sum_thvu(i) + thu(i, k)*(1.+eps*qu(i,k))*dz(i)
1601          sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i)
1602          sum_dq(i) = sum_dq(i) + deltaqw(i, k)*dz(i)
1603          sum_rho(i) = sum_rho(i) + rhow(i, k)*dz(i)
1604          sum_dtdwn(i) = sum_dtdwn(i) + dtdwn(i, k)*dz(i)
1605          sum_dqdwn(i) = sum_dqdwn(i) + dqdwn(i, k)*dz(i)
1606        END IF
1607      END IF
1608    END DO
1609  END DO
1610
1611  DO i = 1, klon
1612    ! cc nrlmd       IF ( wk_adv(i)) THEN
1613    IF (ok_qx_qw(i)) THEN
1614      ! cc
1615      hw0(i) = z(i)
1616    END IF
1617  END DO
1618
1619  ! - WAPE and mean forcing computation
1620  ! -------------------------------------------------------------
1621
1622  ! Means
1623
1624  DO i = 1, klon
1625    ! cc nrlmd       IF ( wk_adv(i)) THEN
1626    IF (ok_qx_qw(i)) THEN
1627      ! cc
1628      av_thu(i) = sum_thu(i)/hw0(i)
1629      av_tu(i) = sum_tu(i)/hw0(i)
1630      av_qu(i) = sum_qu(i)/hw0(i)
1631      av_thvu(i) = sum_thvu(i)/hw0(i)
1632      av_dth(i) = sum_dth(i)/hw0(i)
1633      av_dq(i) = sum_dq(i)/hw0(i)
1634      av_rho(i) = sum_rho(i)/hw0(i)
1635      av_dtdwn(i) = sum_dtdwn(i)/hw0(i)
1636      av_dqdwn(i) = sum_dqdwn(i)/hw0(i)
1637
1638      wape2(i) = -rg*hw0(i)*(av_dth(i)+eps*(av_thu(i)*av_dq(i)+av_dth(i)* &
1639        av_qu(i)+av_dth(i)*av_dq(i)))/av_thvu(i)
1640    END IF
1641  END DO
1642
1643  ! Prognostic variable update
1644  ! ------------------------------------------------------------
1645
1646  ! Filter out bad wakes
1647
1648  DO k = 1, klev
1649    DO i = 1, klon
1650      ! cc nrlmd        IF ( wk_adv(i) .AND. wape2(i) .LT. 0.) THEN
1651      IF (ok_qx_qw(i) .AND. wape2(i)<0.) THEN
1652        ! cc
1653        deltatw(i, k) = 0.
1654        deltaqw(i, k) = 0.
1655        dth(i, k) = 0.
1656      END IF
1657    END DO
1658  END DO
1659
1660
1661  DO i = 1, klon
1662    ! cc nrlmd       IF ( wk_adv(i)) THEN
1663    IF (ok_qx_qw(i)) THEN
1664      ! cc
1665      IF (wape2(i)<0.) THEN
1666        wape2(i) = 0.
1667        cstar2(i) = 0.
1668        hw(i) = hwmin
1669        sigmaw(i) = amax1(sigmad, sigd_con(i))
1670        fip(i) = 0.
1671        gwake(i) = .FALSE.
1672      ELSE
1673        IF (prt_level>=10) PRINT *, 'wape2>0'
1674        cstar2(i) = stark*sqrt(2.*wape2(i))
1675        gwake(i) = .TRUE.
1676      END IF
1677    END IF
1678  END DO
1679
1680  DO i = 1, klon
1681    ! cc nrlmd       IF ( wk_adv(i)) THEN
1682    IF (ok_qx_qw(i)) THEN
1683      ! cc
1684      ktopw(i) = ktop(i)
1685    END IF
1686  END DO
1687
1688  DO i = 1, klon
1689    ! cc nrlmd       IF ( wk_adv(i)) THEN
1690    IF (ok_qx_qw(i)) THEN
1691      ! cc
1692      IF (ktopw(i)>0 .AND. gwake(i)) THEN
1693
1694        ! jyg1     Utilisation d'un h_efficace constant ( ~ feeding layer)
1695        ! cc       heff = 600.
1696        ! Utilisation de la hauteur hw
1697        ! c       heff = 0.7*hw
1698        heff(i) = hw(i)
1699
1700        fip(i) = 0.5*rho(i, ktopw(i))*cstar2(i)**3*heff(i)*2* &
1701          sqrt(sigmaw(i)*wdens(i)*3.14)
1702        fip(i) = alpk*fip(i)
1703        ! jyg2
1704      ELSE
1705        fip(i) = 0.
1706      END IF
1707    END IF
1708  END DO
1709
1710  ! Limitation de sigmaw
1711
1712  ! cc nrlmd
1713  ! DO i=1,klon
1714  ! IF (OK_qx_qw(i)) THEN
1715  ! IF (sigmaw(i).GE.sigmaw_max) sigmaw(i)=sigmaw_max
1716  ! ENDIF
1717  ! ENDDO
1718  ! cc
1719  DO k = 1, klev
1720    DO i = 1, klon
1721
1722      ! cc nrlmd      On maintient désormais constant sigmaw en régime
1723      ! permanent
1724      ! cc      IF ((sigmaw(i).GT.sigmaw_max).or.
1725      IF (((wape(i)>=wape2(i)) .AND. (wape2(i)<=1.0)) .OR. (ktopw(i)<=2) .OR. &
1726          .NOT. ok_qx_qw(i)) THEN
1727        ! cc
1728        dtls(i, k) = 0.
1729        dqls(i, k) = 0.
1730        deltatw(i, k) = 0.
1731        deltaqw(i, k) = 0.
1732      END IF
1733    END DO
1734  END DO
1735
1736  ! cc nrlmd      On maintient désormais constant sigmaw en régime permanent
1737  DO i = 1, klon
1738    IF (((wape(i)>=wape2(i)) .AND. (wape2(i)<=1.0)) .OR. (ktopw(i)<=2) .OR. &
1739        .NOT. ok_qx_qw(i)) THEN
1740      wape(i) = 0.
1741      cstar(i) = 0.
1742!!jyg   Outside subroutine "Wake" hw and sigmaw are zero when there are no wakes
1743!!      hw(i) = hwmin                       !jyg
1744!!      sigmaw(i) = sigmad                  !jyg
1745      hw(i) = 0.                            !jyg
1746      sigmaw(i) = 0.                        !jyg
1747      fip(i) = 0.
1748    ELSE
1749      wape(i) = wape2(i)
1750      cstar(i) = cstar2(i)
1751    END IF
1752    ! c        print*,'wape wape2 ktopw OK_qx_qw =',
1753    ! c     $          wape(i),wape2(i),ktopw(i),OK_qx_qw(i)
1754  END DO
1755
1756
1757  RETURN
1758END SUBROUTINE wake
1759
1760SUBROUTINE wake_vec_modulation(nlon, nl, wk_adv, epsilon, qe, d_qe, deltaqw, &
1761    d_deltaqw, sigmaw, d_sigmaw, alpha)
1762  ! ------------------------------------------------------
1763  ! Dtermination du coefficient alpha tel que les tendances
1764  ! corriges alpha*d_G, pour toutes les grandeurs G, correspondent
1765  ! a une humidite positive dans la zone (x) et dans la zone (w).
1766  ! ------------------------------------------------------
1767  IMPLICIT NONE
1768
1769  ! Input
1770  REAL qe(nlon, nl), d_qe(nlon, nl)
1771  REAL deltaqw(nlon, nl), d_deltaqw(nlon, nl)
1772  REAL sigmaw(nlon), d_sigmaw(nlon)
1773  LOGICAL wk_adv(nlon)
1774  INTEGER nl, nlon
1775  ! Output
1776  REAL alpha(nlon)
1777  ! Internal variables
1778  REAL zeta(nlon, nl)
1779  REAL alpha1(nlon)
1780  REAL x, a, b, c, discrim
1781  REAL epsilon
1782  ! DATA epsilon/1.e-15/
1783  INTEGER i,k
1784
1785  DO k = 1, nl
1786    DO i = 1, nlon
1787      IF (wk_adv(i)) THEN
1788        IF ((deltaqw(i,k)+d_deltaqw(i,k))>=0.) THEN
1789          zeta(i, k) = 0.
1790        ELSE
1791          zeta(i, k) = 1.
1792        END IF
1793      END IF
1794    END DO
1795    DO i = 1, nlon
1796      IF (wk_adv(i)) THEN
1797        x = qe(i, k) + (zeta(i,k)-sigmaw(i))*deltaqw(i, k) + d_qe(i, k) + &
1798          (zeta(i,k)-sigmaw(i))*d_deltaqw(i, k) - d_sigmaw(i)*(deltaqw(i,k)+ &
1799          d_deltaqw(i,k))
1800        a = -d_sigmaw(i)*d_deltaqw(i, k)
1801        b = d_qe(i, k) + (zeta(i,k)-sigmaw(i))*d_deltaqw(i, k) - &
1802          deltaqw(i, k)*d_sigmaw(i)
1803        c = qe(i, k) + (zeta(i,k)-sigmaw(i))*deltaqw(i, k) + epsilon
1804        discrim = b*b - 4.*a*c
1805        ! print*, 'x, a, b, c, discrim', x, a, b, c, discrim
1806        IF (a+b>=0.) THEN !! Condition suffisante pour la positivité de ovap
1807          alpha1(i) = 1.
1808        ELSE
1809          IF (x>=0.) THEN
1810            alpha1(i) = 1.
1811          ELSE
1812            IF (a>0.) THEN
1813              alpha1(i) = 0.9*min((2.*c)/(-b+sqrt(discrim)), (-b+sqrt(discrim &
1814                ))/(2.*a))
1815            ELSE IF (a==0.) THEN
1816              alpha1(i) = 0.9*(-c/b)
1817            ELSE
1818              ! print*,'a,b,c discrim',a,b,c discrim
1819              alpha1(i) = 0.9*max((2.*c)/(-b+sqrt(discrim)), (-b+sqrt(discrim &
1820                ))/(2.*a))
1821            END IF
1822          END IF
1823        END IF
1824        alpha(i) = min(alpha(i), alpha1(i))
1825      END IF
1826    END DO
1827  END DO
1828
1829  RETURN
1830END SUBROUTINE wake_vec_modulation
1831
1832
Note: See TracBrowser for help on using the repository browser.