source: LMDZ5/branches/IPSLCM5A2.1_ISO/libf/phyiso/wake.F90

Last change on this file was 3331, checked in by acozic, 7 years ago

Add modification for isotopes

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