source: LMDZ5/branches/LMDZ6_rc0/libf/phylmd/wake.F90 @ 2687

Last change on this file since 2687 was 2160, checked in by Laurent Fairhead, 10 years ago

Merged trunk changes -r2070:2158 into testing branch. Compilation problems introduced by revision r2155 have been corrected by hand

  • 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: 52.0 KB
Line 
1
2! $Id: wake.F90 2160 2014-11-28 15:36:29Z musat $
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!jyg<
1092!!
1093!!---------------------------------------------------------------
1094!! The change of delta_T due to PBL (vertical diffusion plus thermal plumes)
1095!! is accounted for by the PBL and the Thermals schemes. It is now set to zero
1096!! within the Wake scheme.
1097!!---------------------------------------------------------------
1098          dtPBL(i,k) = 0.
1099          dqPBL(i,k) = 0.
1100!
1101!!          dtPBL(i,k)=wdtPBL(i,k) - udtPBL(i,k)
1102!!          dqPBL(i,k)=wdqPBL(i,k) - udqPBL(i,k)
1103!
1104!!          dtpbl(i, k) = (wdtpbl(i,k)/sigmaw(i)-udtpbl(i,k)/(1.-sigmaw(i)))
1105!!          dqpbl(i, k) = (wdqpbl(i,k)/sigmaw(i)-udqpbl(i,k)/(1.-sigmaw(i)))
1106          ! print*,'dtPBL= ',dtPBL(i,k),' dqPBL= ',dqPBL(i,k)
1107!>jyg
1108!
1109
1110          ! cc nrlmd          Prise en compte du taux de mortalité
1111          ! cc               Définitions de entr, detr
1112          detr(i, k) = 0.
1113
1114          entr(i, k) = detr(i, k) + gfl(i)*cstar(i) + &
1115            sigmaw(i)*(1.-sigmaw(i))*dp_deltomg(i, k)
1116
1117          spread(i, k) = (entr(i,k)-detr(i,k))/sigmaw(i)
1118          ! cc        spread(i,k) =
1119          ! (1.-sigmaw(i))*dp_deltomg(i,k)+gfl(i)*Cstar(i)/
1120          ! cc     $  sigmaw(i)
1121
1122
1123          ! ajout d'un effet onde de gravité -Tgw(k)*deltatw(k) 03/02/06 YU
1124          ! Jingmei
1125
1126          ! write(lunout,*)'wake.F ',i,k, dtimesub,d_deltat_gw(i,k),
1127          ! &  Tgw(i,k),deltatw(i,k)
1128          d_deltat_gw(i, k) = d_deltat_gw(i, k) - tgw(i, k)*deltatw(i, k)* &
1129            dtimesub
1130          ! write(lunout,*)'wake.F ',i,k, dtimesub,d_deltatw(i,k)
1131          ff(i) = d_deltatw(i, k)/dtimesub
1132
1133          ! Sans GW
1134
1135          ! deltatw(k)=deltatw(k)+dtimesub*(ff+dtKE(k)-spread(k)*deltatw(k))
1136
1137          ! GW formule 1
1138
1139          ! deltatw(k) = deltatw(k)+dtimesub*
1140          ! $         (ff+dtKE(k) - spread(k)*deltatw(k)-Tgw(k)*deltatw(k))
1141
1142          ! GW formule 2
1143
1144          IF (dtimesub*tgw(i,k)<1.E-10) THEN
1145            d_deltatw(i, k) = dtimesub*(ff(i)+dtke(i,k)+dtpbl(i,k) & ! cc
1146                                                                     ! $
1147                                                                     ! -spread(i,k)*deltatw(i,k)
1148              -entr(i,k)*deltatw(i,k)/sigmaw(i)-(death_rate(i)*sigmaw( &
1149              i)+detr(i,k))*deltatw(i,k)/(1.-sigmaw(i)) & ! cc
1150              -tgw(i,k)*deltatw(i,k))
1151          ELSE
1152            d_deltatw(i, k) = 1/tgw(i, k)*(1-exp(-dtimesub*tgw(i, &
1153              k)))*(ff(i)+dtke(i,k)+dtpbl(i,k) & ! cc     $
1154                                                 ! -spread(i,k)*deltatw(i,k)
1155              -entr(i,k)*deltatw(i,k)/sigmaw(i)-(death_rate(i)*sigmaw( &
1156              i)+detr(i,k))*deltatw(i,k)/(1.-sigmaw(i)) & ! cc
1157              -tgw(i,k)*deltatw(i,k))
1158          END IF
1159
1160          dth(i, k) = deltatw(i, k)/ppi(i, k)
1161
1162          gg(i) = d_deltaqw(i, k)/dtimesub
1163
1164          d_deltaqw(i, k) = dtimesub*(gg(i)+dqke(i,k)+dqpbl(i,k) & ! cc     $
1165                                                                   ! -spread(i,k)*deltaqw(i,k))
1166            -entr(i,k)*deltaqw(i,k)/sigmaw(i)-(death_rate(i)*sigmaw(i)+detr( &
1167            i,k))*deltaqw(i,k)/(1.-sigmaw(i)))
1168          ! cc
1169
1170          ! cc nrlmd
1171          ! cc       d_deltatw2(i,k)=d_deltatw2(i,k)+d_deltatw(i,k)
1172          ! cc       d_deltaqw2(i,k)=d_deltaqw2(i,k)+d_deltaqw(i,k)
1173          ! cc
1174        END IF
1175      END DO
1176    END DO
1177
1178
1179    ! Scale tendencies so that water vapour remains positive in w and x.
1180
1181    CALL wake_vec_modulation(klon, klev, wk_adv, epsilon, qe, d_qe, deltaqw, &
1182      d_deltaqw, sigmaw, d_sigmaw, alpha)
1183
1184    ! cc nrlmd
1185    ! c      print*,'alpha'
1186    ! c      do i=1,klon
1187    ! c         print*,alpha(i)
1188    ! c      end do
1189    ! cc
1190    DO k = 1, klev
1191      DO i = 1, klon
1192        IF (wk_adv(i) .AND. k<=kupper(i)) THEN
1193          d_te(i, k) = alpha(i)*d_te(i, k)
1194          d_qe(i, k) = alpha(i)*d_qe(i, k)
1195          d_deltatw(i, k) = alpha(i)*d_deltatw(i, k)
1196          d_deltaqw(i, k) = alpha(i)*d_deltaqw(i, k)
1197          d_deltat_gw(i, k) = alpha(i)*d_deltat_gw(i, k)
1198        END IF
1199      END DO
1200    END DO
1201    DO i = 1, klon
1202      IF (wk_adv(i)) THEN
1203        d_sigmaw(i) = alpha(i)*d_sigmaw(i)
1204      END IF
1205    END DO
1206
1207    ! Update large scale variables and wake variables
1208    ! IM 060208 manque DO i + remplace DO k=1,kupper(i)
1209    ! IM 060208     DO k = 1,kupper(i)
1210    DO k = 1, klev
1211      DO i = 1, klon
1212        IF (wk_adv(i) .AND. k<=kupper(i)) THEN
1213          dtls(i, k) = dtls(i, k) + d_te(i, k)
1214          dqls(i, k) = dqls(i, k) + d_qe(i, k)
1215          ! cc nrlmd
1216          d_deltatw2(i, k) = d_deltatw2(i, k) + d_deltatw(i, k)
1217          d_deltaqw2(i, k) = d_deltaqw2(i, k) + d_deltaqw(i, k)
1218          ! cc
1219        END IF
1220      END DO
1221    END DO
1222    DO k = 1, klev
1223      DO i = 1, klon
1224        IF (wk_adv(i) .AND. k<=kupper(i)) THEN
1225          te(i, k) = te0(i, k) + dtls(i, k)
1226          qe(i, k) = qe0(i, k) + dqls(i, k)
1227          the(i, k) = te(i, k)/ppi(i, k)
1228          deltatw(i, k) = deltatw(i, k) + d_deltatw(i, k)
1229          deltaqw(i, k) = deltaqw(i, k) + d_deltaqw(i, k)
1230          dth(i, k) = deltatw(i, k)/ppi(i, k)
1231          ! c      print*,'k,qx,qw',k,qe(i,k)-sigmaw(i)*deltaqw(i,k)
1232          ! c     $        ,qe(i,k)+(1-sigmaw(i))*deltaqw(i,k)
1233        END IF
1234      END DO
1235    END DO
1236    DO i = 1, klon
1237      IF (wk_adv(i)) THEN
1238        sigmaw(i) = sigmaw(i) + d_sigmaw(i)
1239      END IF
1240    END DO
1241
1242
1243    ! Determine Ptop from buoyancy integral
1244    ! ---------------------------------------
1245
1246    ! -     1/ Pressure of the level where dth changes sign.
1247
1248    DO i = 1, klon
1249      IF (wk_adv(i)) THEN
1250        ptop_provis(i) = ph(i, 1)
1251      END IF
1252    END DO
1253
1254    DO k = 2, klev
1255      DO i = 1, klon
1256        IF (wk_adv(i) .AND. ptop_provis(i)==ph(i,1) .AND. &
1257            dth(i,k)>-delta_t_min .AND. dth(i,k-1)<-delta_t_min) THEN
1258          ptop_provis(i) = ((dth(i,k)+delta_t_min)*p(i,k-1)-(dth(i, &
1259            k-1)+delta_t_min)*p(i,k))/(dth(i,k)-dth(i,k-1))
1260        END IF
1261      END DO
1262    END DO
1263
1264    ! -     2/ dth integral
1265
1266    DO i = 1, klon
1267      IF (wk_adv(i)) THEN !!! nrlmd
1268        sum_dth(i) = 0.
1269        dthmin(i) = -delta_t_min
1270        z(i) = 0.
1271      END IF
1272    END DO
1273
1274    DO k = 1, klev
1275      DO i = 1, klon
1276        IF (wk_adv(i)) THEN
1277          dz(i) = -(amax1(ph(i,k+1),ptop_provis(i))-ph(i,k))/(rho(i,k)*rg)
1278          IF (dz(i)>0) THEN
1279            z(i) = z(i) + dz(i)
1280            sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i)
1281            dthmin(i) = amin1(dthmin(i), dth(i,k))
1282          END IF
1283        END IF
1284      END DO
1285    END DO
1286
1287    ! -     3/ height of triangle with area= sum_dth and base = dthmin
1288
1289    DO i = 1, klon
1290      IF (wk_adv(i)) THEN
1291        hw(i) = 2.*sum_dth(i)/amin1(dthmin(i), -0.5)
1292        hw(i) = amax1(hwmin, hw(i))
1293      END IF
1294    END DO
1295
1296    ! -     4/ now, get Ptop
1297
1298    DO i = 1, klon
1299      IF (wk_adv(i)) THEN !!! nrlmd
1300        ktop(i) = 0
1301        z(i) = 0.
1302      END IF
1303    END DO
1304
1305    DO k = 1, klev
1306      DO i = 1, klon
1307        IF (wk_adv(i)) THEN
1308          dz(i) = amin1(-(ph(i,k+1)-ph(i,k))/(rho(i,k)*rg), hw(i)-z(i))
1309          IF (dz(i)>0) THEN
1310            z(i) = z(i) + dz(i)
1311            ptop(i) = ph(i, k) - rho(i, k)*rg*dz(i)
1312            ktop(i) = k
1313          END IF
1314        END IF
1315      END DO
1316    END DO
1317
1318    ! 4.5/Correct ktop and ptop
1319
1320    DO i = 1, klon
1321      IF (wk_adv(i)) THEN
1322        ptop_new(i) = ptop(i)
1323      END IF
1324    END DO
1325
1326    DO k = klev, 2, -1
1327      DO i = 1, klon
1328        ! IM v3JYG; IF (k .GE. ktop(i)
1329        IF (wk_adv(i) .AND. k<=ktop(i) .AND. ptop_new(i)==ptop(i) .AND. &
1330            dth(i,k)>-delta_t_min .AND. dth(i,k-1)<-delta_t_min) THEN
1331          ptop_new(i) = ((dth(i,k)+delta_t_min)*p(i,k-1)-(dth(i, &
1332            k-1)+delta_t_min)*p(i,k))/(dth(i,k)-dth(i,k-1))
1333        END IF
1334      END DO
1335    END DO
1336
1337
1338    DO i = 1, klon
1339      IF (wk_adv(i)) THEN
1340        ptop(i) = ptop_new(i)
1341      END IF
1342    END DO
1343
1344    DO k = klev, 1, -1
1345      DO i = 1, klon
1346        IF (wk_adv(i)) THEN !!! nrlmd
1347          IF (ph(i,k+1)<ptop(i)) ktop(i) = k
1348        END IF
1349      END DO
1350    END DO
1351
1352    ! 5/ Set deltatw & deltaqw to 0 above kupper
1353
1354    DO k = 1, klev
1355      DO i = 1, klon
1356        IF (wk_adv(i) .AND. k>=kupper(i)) THEN
1357          deltatw(i, k) = 0.
1358          deltaqw(i, k) = 0.
1359        END IF
1360      END DO
1361    END DO
1362
1363
1364    ! -------------Cstar computation---------------------------------
1365    DO i = 1, klon
1366      IF (wk_adv(i)) THEN !!! nrlmd
1367        sum_thu(i) = 0.
1368        sum_tu(i) = 0.
1369        sum_qu(i) = 0.
1370        sum_thvu(i) = 0.
1371        sum_dth(i) = 0.
1372        sum_dq(i) = 0.
1373        sum_rho(i) = 0.
1374        sum_dtdwn(i) = 0.
1375        sum_dqdwn(i) = 0.
1376
1377        av_thu(i) = 0.
1378        av_tu(i) = 0.
1379        av_qu(i) = 0.
1380        av_thvu(i) = 0.
1381        av_dth(i) = 0.
1382        av_dq(i) = 0.
1383        av_rho(i) = 0.
1384        av_dtdwn(i) = 0.
1385        av_dqdwn(i) = 0.
1386      END IF
1387    END DO
1388
1389    ! Integrals (and wake top level number)
1390    ! --------------------------------------
1391
1392    ! Initialize sum_thvu to 1st level virt. pot. temp.
1393
1394    DO i = 1, klon
1395      IF (wk_adv(i)) THEN !!! nrlmd
1396        z(i) = 1.
1397        dz(i) = 1.
1398        sum_thvu(i) = thu(i, 1)*(1.+eps*qu(i,1))*dz(i)
1399        sum_dth(i) = 0.
1400      END IF
1401    END DO
1402
1403    DO k = 1, klev
1404      DO i = 1, klon
1405        IF (wk_adv(i)) THEN !!! nrlmd
1406          dz(i) = -(max(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*rg)
1407          IF (dz(i)>0) THEN
1408            z(i) = z(i) + dz(i)
1409            sum_thu(i) = sum_thu(i) + thu(i, k)*dz(i)
1410            sum_tu(i) = sum_tu(i) + tu(i, k)*dz(i)
1411            sum_qu(i) = sum_qu(i) + qu(i, k)*dz(i)
1412            sum_thvu(i) = sum_thvu(i) + thu(i, k)*(1.+eps*qu(i,k))*dz(i)
1413            sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i)
1414            sum_dq(i) = sum_dq(i) + deltaqw(i, k)*dz(i)
1415            sum_rho(i) = sum_rho(i) + rhow(i, k)*dz(i)
1416            sum_dtdwn(i) = sum_dtdwn(i) + dtdwn(i, k)*dz(i)
1417            sum_dqdwn(i) = sum_dqdwn(i) + dqdwn(i, k)*dz(i)
1418          END IF
1419        END IF
1420      END DO
1421    END DO
1422
1423    DO i = 1, klon
1424      IF (wk_adv(i)) THEN !!! nrlmd
1425        hw0(i) = z(i)
1426      END IF
1427    END DO
1428
1429
1430    ! - WAPE and mean forcing computation
1431    ! ---------------------------------------
1432
1433    ! ---------------------------------------
1434
1435    ! Means
1436
1437    DO i = 1, klon
1438      IF (wk_adv(i)) THEN !!! nrlmd
1439        av_thu(i) = sum_thu(i)/hw0(i)
1440        av_tu(i) = sum_tu(i)/hw0(i)
1441        av_qu(i) = sum_qu(i)/hw0(i)
1442        av_thvu(i) = sum_thvu(i)/hw0(i)
1443        av_dth(i) = sum_dth(i)/hw0(i)
1444        av_dq(i) = sum_dq(i)/hw0(i)
1445        av_rho(i) = sum_rho(i)/hw0(i)
1446        av_dtdwn(i) = sum_dtdwn(i)/hw0(i)
1447        av_dqdwn(i) = sum_dqdwn(i)/hw0(i)
1448
1449        wape(i) = -rg*hw0(i)*(av_dth(i)+eps*(av_thu(i)*av_dq(i)+av_dth(i)* &
1450          av_qu(i)+av_dth(i)*av_dq(i)))/av_thvu(i)
1451      END IF
1452    END DO
1453
1454    ! Filter out bad wakes
1455
1456    DO k = 1, klev
1457      DO i = 1, klon
1458        IF (wk_adv(i)) THEN !!! nrlmd
1459          IF (wape(i)<0.) THEN
1460            deltatw(i, k) = 0.
1461            deltaqw(i, k) = 0.
1462            dth(i, k) = 0.
1463          END IF
1464        END IF
1465      END DO
1466    END DO
1467
1468    DO i = 1, klon
1469      IF (wk_adv(i)) THEN !!! nrlmd
1470        IF (wape(i)<0.) THEN
1471          wape(i) = 0.
1472          cstar(i) = 0.
1473          hw(i) = hwmin
1474          sigmaw(i) = max(sigmad, sigd_con(i))
1475          fip(i) = 0.
1476          gwake(i) = .FALSE.
1477        ELSE
1478          cstar(i) = stark*sqrt(2.*wape(i))
1479          gwake(i) = .TRUE.
1480        END IF
1481      END IF
1482    END DO
1483
1484  END DO ! end sub-timestep loop
1485
1486  ! -----------------------------------------------------------------
1487  ! Get back to tendencies per second
1488
1489  DO k = 1, klev
1490    DO i = 1, klon
1491
1492      ! cc nrlmd        IF ( wk_adv(i) .AND. k .LE. kupper(i)) THEN
1493      IF (ok_qx_qw(i) .AND. k<=kupper(i)) THEN
1494        ! cc
1495        dtls(i, k) = dtls(i, k)/dtime
1496        dqls(i, k) = dqls(i, k)/dtime
1497        d_deltatw2(i, k) = d_deltatw2(i, k)/dtime
1498        d_deltaqw2(i, k) = d_deltaqw2(i, k)/dtime
1499        d_deltat_gw(i, k) = d_deltat_gw(i, k)/dtime
1500        ! c      print*,'k,dqls,omg,entr,detr',k,dqls(i,k),omg(i,k),entr(i,k)
1501        ! c     $         ,death_rate(i)*sigmaw(i)
1502      END IF
1503    END DO
1504  END DO
1505
1506
1507  ! ----------------------------------------------------------
1508  ! Determine wake final state; recompute wape, cstar, ktop;
1509  ! filter out bad wakes.
1510  ! ----------------------------------------------------------
1511
1512  ! 2.1 - Undisturbed area and Wake integrals
1513  ! ---------------------------------------------------------
1514
1515  DO i = 1, klon
1516    ! cc nrlmd       if (wk_adv(i)) then !!! nrlmd
1517    IF (ok_qx_qw(i)) THEN
1518      ! cc
1519      z(i) = 0.
1520      sum_thu(i) = 0.
1521      sum_tu(i) = 0.
1522      sum_qu(i) = 0.
1523      sum_thvu(i) = 0.
1524      sum_dth(i) = 0.
1525      sum_dq(i) = 0.
1526      sum_rho(i) = 0.
1527      sum_dtdwn(i) = 0.
1528      sum_dqdwn(i) = 0.
1529
1530      av_thu(i) = 0.
1531      av_tu(i) = 0.
1532      av_qu(i) = 0.
1533      av_thvu(i) = 0.
1534      av_dth(i) = 0.
1535      av_dq(i) = 0.
1536      av_rho(i) = 0.
1537      av_dtdwn(i) = 0.
1538      av_dqdwn(i) = 0.
1539    END IF
1540  END DO
1541  ! Potential temperatures and humidity
1542  ! ----------------------------------------------------------
1543
1544  DO k = 1, klev
1545    DO i = 1, klon
1546      ! cc nrlmd       IF ( wk_adv(i)) THEN
1547      IF (ok_qx_qw(i)) THEN
1548        ! cc
1549        rho(i, k) = p(i, k)/(rd*te(i,k))
1550        IF (k==1) THEN
1551          rhoh(i, k) = ph(i, k)/(rd*te(i,k))
1552          zhh(i, k) = 0
1553        ELSE
1554          rhoh(i, k) = ph(i, k)*2./(rd*(te(i,k)+te(i,k-1)))
1555          zhh(i, k) = (ph(i,k)-ph(i,k-1))/(-rhoh(i,k)*rg) + zhh(i, k-1)
1556        END IF
1557        the(i, k) = te(i, k)/ppi(i, k)
1558        thu(i, k) = (te(i,k)-deltatw(i,k)*sigmaw(i))/ppi(i, k)
1559        tu(i, k) = te(i, k) - deltatw(i, k)*sigmaw(i)
1560        qu(i, k) = qe(i, k) - deltaqw(i, k)*sigmaw(i)
1561        rhow(i, k) = p(i, k)/(rd*(te(i,k)+deltatw(i,k)))
1562        dth(i, k) = deltatw(i, k)/ppi(i, k)
1563      END IF
1564    END DO
1565  END DO
1566
1567  ! Integrals (and wake top level number)
1568  ! -----------------------------------------------------------
1569
1570  ! Initialize sum_thvu to 1st level virt. pot. temp.
1571
1572  DO i = 1, klon
1573    ! cc nrlmd       IF ( wk_adv(i)) THEN
1574    IF (ok_qx_qw(i)) THEN
1575      ! cc
1576      z(i) = 1.
1577      dz(i) = 1.
1578      sum_thvu(i) = thu(i, 1)*(1.+eps*qu(i,1))*dz(i)
1579      sum_dth(i) = 0.
1580    END IF
1581  END DO
1582
1583  DO k = 1, klev
1584    DO i = 1, klon
1585      ! cc nrlmd       IF ( wk_adv(i)) THEN
1586      IF (ok_qx_qw(i)) THEN
1587        ! cc
1588        dz(i) = -(amax1(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*rg)
1589        IF (dz(i)>0) THEN
1590          z(i) = z(i) + dz(i)
1591          sum_thu(i) = sum_thu(i) + thu(i, k)*dz(i)
1592          sum_tu(i) = sum_tu(i) + tu(i, k)*dz(i)
1593          sum_qu(i) = sum_qu(i) + qu(i, k)*dz(i)
1594          sum_thvu(i) = sum_thvu(i) + thu(i, k)*(1.+eps*qu(i,k))*dz(i)
1595          sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i)
1596          sum_dq(i) = sum_dq(i) + deltaqw(i, k)*dz(i)
1597          sum_rho(i) = sum_rho(i) + rhow(i, k)*dz(i)
1598          sum_dtdwn(i) = sum_dtdwn(i) + dtdwn(i, k)*dz(i)
1599          sum_dqdwn(i) = sum_dqdwn(i) + dqdwn(i, k)*dz(i)
1600        END IF
1601      END IF
1602    END DO
1603  END DO
1604
1605  DO i = 1, klon
1606    ! cc nrlmd       IF ( wk_adv(i)) THEN
1607    IF (ok_qx_qw(i)) THEN
1608      ! cc
1609      hw0(i) = z(i)
1610    END IF
1611  END DO
1612
1613  ! - WAPE and mean forcing computation
1614  ! -------------------------------------------------------------
1615
1616  ! Means
1617
1618  DO i = 1, klon
1619    ! cc nrlmd       IF ( wk_adv(i)) THEN
1620    IF (ok_qx_qw(i)) THEN
1621      ! cc
1622      av_thu(i) = sum_thu(i)/hw0(i)
1623      av_tu(i) = sum_tu(i)/hw0(i)
1624      av_qu(i) = sum_qu(i)/hw0(i)
1625      av_thvu(i) = sum_thvu(i)/hw0(i)
1626      av_dth(i) = sum_dth(i)/hw0(i)
1627      av_dq(i) = sum_dq(i)/hw0(i)
1628      av_rho(i) = sum_rho(i)/hw0(i)
1629      av_dtdwn(i) = sum_dtdwn(i)/hw0(i)
1630      av_dqdwn(i) = sum_dqdwn(i)/hw0(i)
1631
1632      wape2(i) = -rg*hw0(i)*(av_dth(i)+eps*(av_thu(i)*av_dq(i)+av_dth(i)* &
1633        av_qu(i)+av_dth(i)*av_dq(i)))/av_thvu(i)
1634    END IF
1635  END DO
1636
1637  ! Prognostic variable update
1638  ! ------------------------------------------------------------
1639
1640  ! Filter out bad wakes
1641
1642  DO k = 1, klev
1643    DO i = 1, klon
1644      ! cc nrlmd        IF ( wk_adv(i) .AND. wape2(i) .LT. 0.) THEN
1645      IF (ok_qx_qw(i) .AND. wape2(i)<0.) THEN
1646        ! cc
1647        deltatw(i, k) = 0.
1648        deltaqw(i, k) = 0.
1649        dth(i, k) = 0.
1650      END IF
1651    END DO
1652  END DO
1653
1654
1655  DO i = 1, klon
1656    ! cc nrlmd       IF ( wk_adv(i)) THEN
1657    IF (ok_qx_qw(i)) THEN
1658      ! cc
1659      IF (wape2(i)<0.) THEN
1660        wape2(i) = 0.
1661        cstar2(i) = 0.
1662        hw(i) = hwmin
1663        sigmaw(i) = amax1(sigmad, sigd_con(i))
1664        fip(i) = 0.
1665        gwake(i) = .FALSE.
1666      ELSE
1667        IF (prt_level>=10) PRINT *, 'wape2>0'
1668        cstar2(i) = stark*sqrt(2.*wape2(i))
1669        gwake(i) = .TRUE.
1670      END IF
1671    END IF
1672  END DO
1673
1674  DO i = 1, klon
1675    ! cc nrlmd       IF ( wk_adv(i)) THEN
1676    IF (ok_qx_qw(i)) THEN
1677      ! cc
1678      ktopw(i) = ktop(i)
1679    END IF
1680  END DO
1681
1682  DO i = 1, klon
1683    ! cc nrlmd       IF ( wk_adv(i)) THEN
1684    IF (ok_qx_qw(i)) THEN
1685      ! cc
1686      IF (ktopw(i)>0 .AND. gwake(i)) THEN
1687
1688        ! jyg1     Utilisation d'un h_efficace constant ( ~ feeding layer)
1689        ! cc       heff = 600.
1690        ! Utilisation de la hauteur hw
1691        ! c       heff = 0.7*hw
1692        heff(i) = hw(i)
1693
1694        fip(i) = 0.5*rho(i, ktopw(i))*cstar2(i)**3*heff(i)*2* &
1695          sqrt(sigmaw(i)*wdens(i)*3.14)
1696        fip(i) = alpk*fip(i)
1697        ! jyg2
1698      ELSE
1699        fip(i) = 0.
1700      END IF
1701    END IF
1702  END DO
1703
1704  ! Limitation de sigmaw
1705
1706  ! cc nrlmd
1707  ! DO i=1,klon
1708  ! IF (OK_qx_qw(i)) THEN
1709  ! IF (sigmaw(i).GE.sigmaw_max) sigmaw(i)=sigmaw_max
1710  ! ENDIF
1711  ! ENDDO
1712  ! cc
1713  DO k = 1, klev
1714    DO i = 1, klon
1715
1716      ! cc nrlmd      On maintient désormais constant sigmaw en régime
1717      ! permanent
1718      ! cc      IF ((sigmaw(i).GT.sigmaw_max).or.
1719      IF (((wape(i)>=wape2(i)) .AND. (wape2(i)<=1.0)) .OR. (ktopw(i)<=2) .OR. &
1720          .NOT. ok_qx_qw(i)) THEN
1721        ! cc
1722        dtls(i, k) = 0.
1723        dqls(i, k) = 0.
1724        deltatw(i, k) = 0.
1725        deltaqw(i, k) = 0.
1726      END IF
1727    END DO
1728  END DO
1729
1730  ! cc nrlmd      On maintient désormais constant sigmaw en régime permanent
1731  DO i = 1, klon
1732    IF (((wape(i)>=wape2(i)) .AND. (wape2(i)<=1.0)) .OR. (ktopw(i)<=2) .OR. &
1733        .NOT. ok_qx_qw(i)) THEN
1734      wape(i) = 0.
1735      cstar(i) = 0.
1736      hw(i) = hwmin
1737      sigmaw(i) = sigmad
1738      fip(i) = 0.
1739    ELSE
1740      wape(i) = wape2(i)
1741      cstar(i) = cstar2(i)
1742    END IF
1743    ! c        print*,'wape wape2 ktopw OK_qx_qw =',
1744    ! c     $          wape(i),wape2(i),ktopw(i),OK_qx_qw(i)
1745  END DO
1746
1747
1748  RETURN
1749END SUBROUTINE wake
1750
1751SUBROUTINE wake_vec_modulation(nlon, nl, wk_adv, epsilon, qe, d_qe, deltaqw, &
1752    d_deltaqw, sigmaw, d_sigmaw, alpha)
1753  ! ------------------------------------------------------
1754  ! Dtermination du coefficient alpha tel que les tendances
1755  ! corriges alpha*d_G, pour toutes les grandeurs G, correspondent
1756  ! a une humidite positive dans la zone (x) et dans la zone (w).
1757  ! ------------------------------------------------------
1758
1759
1760  ! Input
1761  REAL qe(nlon, nl), d_qe(nlon, nl)
1762  REAL deltaqw(nlon, nl), d_deltaqw(nlon, nl)
1763  REAL sigmaw(nlon), d_sigmaw(nlon)
1764  LOGICAL wk_adv(nlon)
1765  INTEGER nl, nlon
1766  ! Output
1767  REAL alpha(nlon)
1768  ! Internal variables
1769  REAL zeta(nlon, nl)
1770  REAL alpha1(nlon)
1771  REAL x, a, b, c, discrim
1772  REAL epsilon
1773  ! DATA epsilon/1.e-15/
1774
1775  DO k = 1, nl
1776    DO i = 1, nlon
1777      IF (wk_adv(i)) THEN
1778        IF ((deltaqw(i,k)+d_deltaqw(i,k))>=0.) THEN
1779          zeta(i, k) = 0.
1780        ELSE
1781          zeta(i, k) = 1.
1782        END IF
1783      END IF
1784    END DO
1785    DO i = 1, nlon
1786      IF (wk_adv(i)) THEN
1787        x = qe(i, k) + (zeta(i,k)-sigmaw(i))*deltaqw(i, k) + d_qe(i, k) + &
1788          (zeta(i,k)-sigmaw(i))*d_deltaqw(i, k) - d_sigmaw(i)*(deltaqw(i,k)+ &
1789          d_deltaqw(i,k))
1790        a = -d_sigmaw(i)*d_deltaqw(i, k)
1791        b = d_qe(i, k) + (zeta(i,k)-sigmaw(i))*d_deltaqw(i, k) - &
1792          deltaqw(i, k)*d_sigmaw(i)
1793        c = qe(i, k) + (zeta(i,k)-sigmaw(i))*deltaqw(i, k) + epsilon
1794        discrim = b*b - 4.*a*c
1795        ! print*, 'x, a, b, c, discrim', x, a, b, c, discrim
1796        IF (a+b>=0.) THEN !! Condition suffisante pour la positivité de ovap
1797          alpha1(i) = 1.
1798        ELSE
1799          IF (x>=0.) THEN
1800            alpha1(i) = 1.
1801          ELSE
1802            IF (a>0.) THEN
1803              alpha1(i) = 0.9*min((2.*c)/(-b+sqrt(discrim)), (-b+sqrt(discrim &
1804                ))/(2.*a))
1805            ELSE IF (a==0.) THEN
1806              alpha1(i) = 0.9*(-c/b)
1807            ELSE
1808              ! print*,'a,b,c discrim',a,b,c discrim
1809              alpha1(i) = 0.9*max((2.*c)/(-b+sqrt(discrim)), (-b+sqrt(discrim &
1810                ))/(2.*a))
1811            END IF
1812          END IF
1813        END IF
1814        alpha(i) = min(alpha(i), alpha1(i))
1815      END IF
1816    END DO
1817  END DO
1818
1819  RETURN
1820END SUBROUTINE wake_vec_modulation
1821
1822
Note: See TracBrowser for help on using the repository browser.