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

Last change on this file since 4799 was 2922, checked in by fhourdin, 7 years ago

Nouvelle option iflag_wk_check_trgl pour tester la bonne sante des poches.
L'ancien flag logique flag_wk_check_trgl est toujours utilisable.
Il vaut y/n si iflag_wk_check_trgl=1/0.
L'option 1 avait tendance à faire disparaitre les poches dans des cas 1D.
L'option 2 est donc précaunisée pour les réglages à venir.
Fredho pour Jean-Yves

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