source: LMDZ5/branches/AI-cosp/libf/phylmd/wake.F90 @ 4049

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

Physics/dynamics separation:

  • remove all references to dimensions.h from physics. nbp_lon (==iim) , nbp_lat (==jjm+1) and nbp_lev (==llm) from mod_grid_phy_lmdz should be used instead.
  • added module regular_lonlat_mod in phy_common to store information about the global (lon-lat) grid cell boundaries and centers.

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