source: LMDZ5/branches/testing/libf/phylmd/wake.F90 @ 4674

Last change on this file since 4674 was 2787, checked in by Laurent Fairhead, 8 years ago

Merged trunk changes r2727:2785 into testing branch

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