Changeset 5255 for LMDZ6/trunk


Ignore:
Timestamp:
Oct 22, 2024, 3:35:01 PM (2 months ago)
Author:
crisi
Message:

a few updates in isotopes

Location:
LMDZ6/trunk/libf/phylmdiso
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • LMDZ6/trunk/libf/phylmdiso/calwake.F90

    r4783 r5255  
    66    dt_dwn, dq_dwn, m_dwn, m_up, dt_a, dq_a, wgen, &
    77    sigd, Cin, &
    8     wake_deltat, wake_deltaq, wake_s, awake_dens, wake_dens, &
     8    wake_deltat, wake_deltaq, wake_s, awake_s, wake_dens, awake_dens, &
    99    wake_dth, wake_h, &
    1010    wake_pe, wake_fip, wake_gfl, &
    … …  
    1414    wake_omg, wake_dp_deltomg, &
    1515    wake_spread, wake_cstar, wake_d_deltat_gw, &
    16     wake_ddeltat, wake_ddeltaq, wake_ds, awake_ddens, wake_ddens &
     16    wake_ddeltat, wake_ddeltaq, wake_ds, awake_ds, wake_ddens, awake_ddens &
    1717#ifdef ISO
    1818     &               ,xt,dxt_dwn,dxt_a &
    … …  
    6767  ! ------------
    6868  REAL, DIMENSION(klon, klev),   INTENT (INOUT)      :: wake_deltat, wake_deltaq
    69   REAL, DIMENSION(klon),         INTENT (INOUT)      :: wake_s
    70   REAL, DIMENSION(klon),         INTENT (INOUT)      :: awake_dens, wake_dens
     69  REAL, DIMENSION(klon),         INTENT (INOUT)      :: wake_s, awake_s
     70  REAL, DIMENSION(klon),         INTENT (INOUT)      :: wake_dens, awake_dens
    7171#ifdef ISO
    7272  REAL, DIMENSION(ntraciso,klon, klev),   INTENT (INOUT)      :: wake_deltaxt
    … …  
    8888  REAL, DIMENSION(klon),         INTENT (OUT)        :: wake_cstar
    8989  REAL, DIMENSION(klon, klev),   INTENT (OUT)        :: wake_ddeltat, wake_ddeltaq
    90   REAL, DIMENSION(klon),         INTENT (OUT)        :: wake_ds, awake_ddens, wake_ddens
     90  REAL, DIMENSION(klon),         INTENT (OUT)        :: wake_ds, awake_ds, wake_ddens, awake_ddens
    9191#ifdef ISO
    9292  REAL, DIMENSION(ntraciso,klon, klev),   INTENT (OUT)        :: dxt_wake
    … …  
    115115  REAL, DIMENSION(klon, klev)                        :: tx, qx
    116116  REAL, DIMENSION(klon)                              :: hw, wape, fip, gfl
    117   REAL, DIMENSION(klon)                              :: sigmaw, awdens, wdens
     117  REAL, DIMENSION(klon)                              :: sigmaw, asigmaw, wdens, awdens
    118118  REAL, DIMENSION(klon, klev)                        :: omgbdth
    119119  REAL, DIMENSION(klon, klev)                        :: dp_omgb
    … …  
    126126  REAL, DIMENSION(klon, klev)                        :: d_deltat_gw
    127127  REAL, DIMENSION(klon, klev)                        :: d_deltatw, d_deltaqw
    128   REAL, DIMENSION(klon)                              :: d_sigmaw, d_awdens, d_wdens
     128  REAL, DIMENSION(klon)                              :: d_sigmaw, d_asigmaw, d_wdens, d_awdens
    129129#ifdef ISO
    130130  REAL, DIMENSION(ntraciso,klon, klev)                        :: xte
    … …  
    142142
    143143  IF (prt_level >= 10) THEN
    144     print *, '-> calwake, wake_s, wgen input ', wake_s(1), wgen(1)
     144    print *, '-> calwake, wake_s, awake_s, wgen input ', wake_s(1), awake_s(1), wgen(1)
    145145  ENDIF
    146146
    … …  
    191191d_deltaqw(:,:) = 0.
    192192d_sigmaw(:) = 0.
     193d_asigmaw(:) = 0.
     194d_wdens(:) = 0.
    193195d_awdens(:) = 0.
    194 d_wdens(:) = 0.
    195196#ifdef ISO
    196197dxtls(:,:,:) = 0.
    … …  
    227228
    228229  DO i = 1, klon
    229     sigmaw(i) = wake_s(i)
    230   END DO
    231 
    232   DO i = 1, klon
     230    sigmaw(i)  = wake_s(i)
     231    asigmaw(i) = awake_s(i)
     232  END DO
     233
     234  DO i = 1, klon
     235    wdens(i) = max(0., wake_dens(i))
    233236    awdens(i) = max(0., awake_dens(i))
    234     wdens(i) = max(0., wake_dens(i))
    235237  END DO
    236238
    … …  
    272274#endif
    273275
    274   CALL wake(znatsurf, p, ph, pi, dtime, &
     276
     277  CALL wake(klon,klev,znatsurf, p, ph, pi, dtime, &
    275278    te, qe, omgbe, &
    276279    dtdwn, dqdwn, amdwn, amup, dta, dqa, wgen, &
    277280    sigd0, Cin, &
    278     dtw, dqw, sigmaw, awdens, wdens, &                                   ! state variables
     281    dtw, dqw, sigmaw, asigmaw, wdens, awdens, &                      ! state variables
    279282    dth, hw, wape, fip, gfl, &
    280283    dtls, dqls, ktopw, omgbdth, dp_omgb, tx, qx, &
    281284    dtke, dqke, omg, dp_deltomg, spread, cstar, &
    282285    d_deltat_gw, &
    283     d_deltatw, d_deltaqw, d_sigmaw, d_awdens, d_wdens &                     ! tendencies
     286    d_deltatw, d_deltaqw, d_sigmaw, d_asigmaw, d_wdens, d_awdens &      ! tendencies
    284287#ifdef ISO
    285288     , xte,dxtdwn,dxta,dxtw &
    … …  
    386389    IF (ktopw(i)>0) THEN
    387390      wake_ds(i) = d_sigmaw(i)*dtime
     391      awake_ds(i) = d_asigmaw(i)*dtime
    388392      awake_ddens(i) = d_awdens(i)*dtime
    389393      wake_ddens(i) = d_wdens(i)*dtime
    390394    ELSE
    391       wake_ds(i)   = -wake_s(i)
    392       wake_ddens(i)= -wake_dens(i)
     395      wake_ds(i)    = -wake_s(i)
     396      awake_ds(i)   = -awake_s(i)
     397      wake_ddens(i) = -wake_dens(i)
     398      awake_ddens(i)= -awake_dens(i)
    393399    END IF
    394400  END DO
    … …  
    421427    DO i = 1, klon
    422428      wake_s(i) = sigmaw(i)
     429      awake_s(i) = asigmaw(i)
    423430      awake_dens(i) = awdens(i)
    424431      wake_dens(i) = wdens(i)
    … …  
    435442  ENDIF  ! (first)
    436443!>jyg
     444  IF (prt_level >= 10) THEN
     445    print *, 'calwake ->, wake_s, awake_s ', wake_s(1), awake_s(1)
     446  ENDIF
    437447
    438448  RETURN
    439449END SUBROUTINE calwake
    440450
     451

  • LMDZ6/trunk/libf/phylmdiso/isotopes_routines_mod.F90

    r5199 r5255  
    1882218822        if ((iso_HDO.gt.0).and.(ixt.eq.iso_HDO)) then
    1882318823            if (q.gt.ridicule) then
    18824                     write(*,*) 'xt,q=',xt,q
    18825                     write(*,*) 'alpha=',alpha
    18826                     write(*,*) 'toce,kcin,h0=',toce,kcin,h0
    18827                     write(*,*) 'RMerlivat=',RMerlivat
     18824                    !write(*,*) 'xt,q=',xt,q
     18825                    !write(*,*) 'alpha=',alpha
     18826                    !write(*,*) 'toce,kcin,h0=',toce,kcin,h0
     18827                    !write(*,*) 'RMerlivat=',RMerlivat
    1882818828                call iso_verif_aberrant_encadre( xt/q, 'isotopes_routines_mod 18930b: iso_init_ideal')
    1882918829            endif

  • LMDZ6/trunk/libf/phylmdiso/lmdz_wake.F90

    r4594 r5255  
    11MODULE lmdz_wake
    2 ! $Id: wake.F90 3648 2020-03-16 15:36:59Z jghattas $
     2
     3! $Id: lmdz_wake.F90 4908 2024-04-15 17:30:55Z jyg $
     4
    35CONTAINS
    4 SUBROUTINE wake(znatsurf, p, ph, pi, dtime, &
    5                 te0, qe0, omgb, &
     6
     7SUBROUTINE wake(klon,klev,znatsurf, p, ph, pi, dtime, &
     8                tb0, qb0, omgb, &
    69                dtdwn, dqdwn, amdwn, amup, dta, dqa, wgen, &
    710                sigd_con, Cin, &
    8                 deltatw, deltaqw, sigmaw, awdens, wdens, &                  ! state variables
     11                deltatw, deltaqw, sigmaw, asigmaw, wdens, awdens, &                  ! state variables
    912                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_awdens2, d_wdens2 &     ! tendencies
     13                dtls, dqls, ktopw, omgbdth, dp_omgb, tx, qx, &
     14                dtke, dqke, omg, dp_deltomg, wkspread, cstar, &
     15                d_deltat_gw, &                                                      ! tendencies
     16                d_deltatw2, d_deltaqw2, d_sigmaw2, d_asigmaw2, d_wdens2, d_awdens2 &     ! tendencies
    1417
    1518#ifdef ISO
    16                      ,xte0,dxtdwn,dxta,deltaxtw &
    17                      ,dxtls,xtu,dxtke,d_deltaxtw2 &
     19                     ,xtb0,dxtdwn,dxta,deltaxtw &
     20                     ,dxtls,xtx,dxtke,d_deltaxtw2 &
    1821#endif                 
    1922                ) 
    … …  
    2932  ! **************************************************************
    3033
    31   USE ioipsl_getin_p_mod, ONLY : getin_p
    32   USE dimphy
    33   use mod_phys_lmdz_para
    34   USE print_control_mod, ONLY: prt_level
     34
     35  USE lmdz_wake_ini , ONLY : wake_ini
     36  USE lmdz_wake_ini , ONLY : prt_level,epsim1,RG,RD
     37  USE lmdz_wake_ini , ONLY : stark, wdens_ref, coefgw, alpk, wk_pupper
     38  USE lmdz_wake_ini , ONLY : crep_upper, crep_sol, tau_cv, rzero, aa0, flag_wk_check_trgl
     39  USE lmdz_wake_ini , ONLY : ok_bug_gfl
     40  USE lmdz_wake_ini , ONLY : iflag_wk_act, iflag_wk_check_trgl, iflag_wk_pop_dyn, wdensinit, wdensthreshold
     41  USE lmdz_wake_ini , ONLY : sigmad, hwmin, wapecut, cstart, sigmaw_max, dens_rate, epsilon_loc
     42  USE lmdz_wake_ini , ONLY : iflag_wk_profile
     43  USE lmdz_wake_ini , ONLY : smallestreal,wk_nsub
     44
     45
    3546#ifdef ISO
    3647  USE infotrac_phy, ONLY : ntraciso=>ntiso
    … …  
    5364  ! deltaqw    : specific humidity difference between wake and off-wake regions
    5465  ! sigmaw     : fractional area covered by wakes.
     66  ! asigmaw    : fractional area covered by active wakes.
    5567  ! wdens      : number of wakes per unit area
     68  ! awdens     : number of active wakes per unit area
    5669
    5770  ! Variable de sortie :
    … …  
    7184  ! omg    : Delta_omg =vertical velocity diff. wake-undist. (Pa/s)
    7285  ! dp_deltomg  : vertical gradient of omg (s-1)
    73   ! spread  : spreading term in d_t_wake and d_q_wake
     86  ! wkspread  : spreading term in d_t_wake and d_q_wake
    7487  ! deltatw     : updated temperature difference (T_w-T_u).
    7588  ! deltaqw     : updated humidity difference (q_w-q_u).
    7689  ! sigmaw      : updated wake fractional area.
     90  ! asigmaw     : updated active wake fractional area.
    7791  ! d_deltat_gw : delta T tendency due to GW
    7892
    … …  
    8094
    8195  ! aire : aire de la maille
    82   ! te0  : temperature dans l'environnement  (K)
    83   ! qe0  : humidite dans l'environnement     (kg/kg)
     96  ! tb0  : horizontal average of temperature  (K)
     97  ! qb0  : horizontal average of humidity   (kg/kg)
    8498  ! omgb : vitesse verticale moyenne sur la maille (Pa/s)
    8599  ! dtdwn: source de chaleur due aux descentes (K/s)
    … …  
    101115  ! Variables internes :
    102116
    103   ! rhow : masse volumique de la poche froide
    104   ! rho  : environment density at P levels
    105   ! rhoh : environment density at Ph levels
    106   ! te   : environment temperature | may change within
    107   ! qe   : environment humidity    | sub-time-stepping
    108   ! the  : environment potential temperature
    109   ! thu  : potential temperature in undisturbed area
    110   ! tu   :  temperature  in undisturbed area
    111   ! qu   : humidity in undisturbed area
     117  ! rho  : mean density at P levels
     118  ! rhoh : mean density at Ph levels
     119  ! tb   : mean temperature | may change within
     120  ! qb   : mean humidity    | sub-time-stepping
     121  ! thb  : mean potential temperature
     122  ! thx  : potential temperature in (x) area
     123  ! tx   : temperature  in (x) area
     124  ! qx   : humidity in (x) area
    112125  ! dp_omgb: vertical gradient og LS omega
    113126  ! omgbw  : wake average vertical omega
    … …  
    115128  ! omgbdq : flux of Delta_q transported by LS omega
    116129  ! dth  : potential temperature diff. wake-undist.
    117   ! th1  : first pot. temp. for vertical advection (=thu)
     130  ! th1  : first pot. temp. for vertical advection (=thx)
    118131  ! th2  : second pot. temp. for vertical advection (=thw)
    119132  ! q1   : first humidity for vertical advection
    120133  ! q2   : second humidity for vertical advection
    121   ! d_deltatw   : terme de redistribution pour deltatw
    122   ! d_deltaqw   : terme de redistribution pour deltaqw
    123   ! deltatw0   : deltatw initial
    124   ! deltaqw0   : deltaqw initial
     134  ! d_deltatw   : redistribution term for deltatw
     135  ! d_deltaqw   : redistribution term for deltaqw
     136  ! deltatw0   : initial deltatw
     137  ! deltaqw0   : initial deltaqw
    125138  ! hw0    : wake top hight (defined as the altitude at which deltatw=0)
    126139  ! amflux : horizontal mass flux through wake boundary
    127140  ! wdens_ref: initial number of wakes per unit area (3D) or per
    128   ! unit length (2D), at the beginning of each time step
    129   ! Tgw    : 1 sur la période de onde de gravité
    130   ! Cgw    : vitesse de propagation de onde de gravité
    131   ! LL     : distance entre 2 poches
     141  !            unit length (2D), at the beginning of each time step
     142  ! Tgw    : 1 sur la periode de onde de gravite
     143  ! Cgw    : vitesse de propagation de onde de gravite
     144  ! LL     : distance between 2 wakes
    132145
    133146  ! -------------------------------------------------------------------------
    134   ! Déclaration de variables
     147  ! Declaration de variables
    135148  ! -------------------------------------------------------------------------
    136149
    137   include "YOMCST.h"
    138   include "cvthermo.h"
    139150
    140151  ! Arguments en entree
    141152  ! --------------------
    142153
     154  INTEGER, INTENT(IN) :: klon,klev
    143155  INTEGER, DIMENSION (klon),        INTENT(IN)          :: znatsurf
    144156  REAL, DIMENSION (klon, klev),     INTENT(IN)          :: p, pi
    … …  
    146158  REAL, DIMENSION (klon, klev),     INTENT(IN)          :: omgb
    147159  REAL,                             INTENT(IN)          :: dtime
    148   REAL, DIMENSION (klon, klev),     INTENT(IN)          :: te0, qe0
     160  REAL, DIMENSION (klon, klev),     INTENT(IN)          :: tb0, qb0
    149161  REAL, DIMENSION (klon, klev),     INTENT(IN)          :: dtdwn, dqdwn
    150162  REAL, DIMENSION (klon, klev),     INTENT(IN)          :: amdwn, amup
    … …  
    154166  REAL, DIMENSION (klon),           INTENT(IN)          :: Cin
    155167#ifdef ISO
    156   REAL, DIMENSION (ntraciso,klon, klev),     INTENT(IN)          :: xte0
     168  REAL, DIMENSION (ntraciso,klon, klev),     INTENT(IN)          :: xtb0
    157169  REAL, DIMENSION (ntraciso,klon, klev),     INTENT(IN)          :: dxtdwn
    158170  REAL, DIMENSION (ntraciso,klon, klev),     INTENT(IN)          :: dxta
    … …  
    164176  REAL, DIMENSION (klon, klev),     INTENT(INOUT)       :: deltatw, deltaqw
    165177  REAL, DIMENSION (klon),           INTENT(INOUT)       :: sigmaw
     178  REAL, DIMENSION (klon),           INTENT(INOUT)       :: asigmaw
     179  REAL, DIMENSION (klon),           INTENT(INOUT)       :: wdens
    166180  REAL, DIMENSION (klon),           INTENT(INOUT)       :: awdens
    167   REAL, DIMENSION (klon),           INTENT(INOUT)       :: wdens
    168181#ifdef ISO
    169182  REAL, DIMENSION (ntraciso,klon, klev),     INTENT(INOUT)       :: deltaxtw
    … …  
    174187
    175188  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: dth
    176   REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: tu, qu
     189  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: tx, qx
    177190  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: dtls, dqls
    178191  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: dtke, dqke
    179   REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: spread    !  unused (jyg)
     192  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: wkspread    !  unused (jyg)
    180193  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: omgbdth, omg
    181194  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: dp_omgb, dp_deltomg
    182   REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: d_deltat_gw
    183195  REAL, DIMENSION (klon),           INTENT(OUT)         :: hw, wape, fip, gfl, cstar
    184196  INTEGER, DIMENSION (klon),        INTENT(OUT)         :: ktopw
    185   ! Tendencies of state variables
     197  ! Tendencies of state variables (2 is appended to the names of fields which are the cumul of fields
     198  !                                 computed at each sub-timestep; e.g. d_wdens2 is the cumul of d_wdens)
     199  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: d_deltat_gw
    186200  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: d_deltatw2, d_deltaqw2
    187   REAL, DIMENSION (klon),           INTENT(OUT)         :: d_sigmaw2, d_awdens2, d_wdens2
     201  REAL, DIMENSION (klon),           INTENT(OUT)         :: d_sigmaw2, d_asigmaw2, d_wdens2, d_awdens2
     202
    188203#ifdef ISO
    189   REAL, DIMENSION (ntraciso,klon, klev),     INTENT(OUT)         :: xtu
     204  REAL, DIMENSION (ntraciso,klon, klev),     INTENT(OUT)         :: xtx
    190205  REAL, DIMENSION (ntraciso,klon, klev),     INTENT(OUT)         :: dxtls
    191206  REAL, DIMENSION (ntraciso,klon, klev),     INTENT(OUT)         :: dxtke
    … …  
    196211  ! -------------------
    197212
    198   ! Variables à fixer
    199   INTEGER, SAVE                                         :: igout
    200   !$OMP THREADPRIVATE(igout)
    201   LOGICAL, SAVE                                         :: first = .TRUE.
    202   !$OMP THREADPRIVATE(first)
    203 !jyg<
    204 !!  REAL, SAVE                                            :: stark, wdens_ref, coefgw, alpk
    205   REAL, SAVE, DIMENSION(2)                              :: wdens_ref
    206   REAL, SAVE                                            :: stark, coefgw, alpk
    207 !>jyg
    208   REAL, SAVE                                            :: crep_upper, crep_sol 
    209   !$OMP THREADPRIVATE(stark, wdens_ref, coefgw, alpk, crep_upper, crep_sol)
    210 
    211   REAL, SAVE                                            :: tau_cv
    212   !$OMP THREADPRIVATE(tau_cv)
    213   REAL, SAVE                                            :: rzero, aa0 ! minimal wake radius and area
    214   !$OMP THREADPRIVATE(rzero, aa0)
    215 
    216   LOGICAL, SAVE                                         :: flag_wk_check_trgl
    217   !$OMP THREADPRIVATE(flag_wk_check_trgl)
    218   INTEGER, SAVE                                         :: iflag_wk_check_trgl
    219   !$OMP THREADPRIVATE(iflag_wk_check_trgl)
    220   INTEGER, SAVE                                         :: iflag_wk_pop_dyn
    221   !$OMP THREADPRIVATE(iflag_wk_pop_dyn)
     213  ! Variables a fixer
    222214
    223215  REAL                                                  :: delta_t_min
    224   INTEGER                                               :: nsub
    225216  REAL                                                  :: dtimesub
    226   REAL, SAVE                                            :: wdensmin
    227   !$OMP THREADPRIVATE(wdensmin)
    228   REAL, SAVE                                            :: sigmad, hwmin, wapecut, cstart
    229   !$OMP THREADPRIVATE(sigmad, hwmin, wapecut, cstart)
    230   REAL, SAVE                                            :: sigmaw_max
    231   !$OMP THREADPRIVATE(sigmaw_max) 
    232   REAL, SAVE                                            :: dens_rate
    233   !$OMP THREADPRIVATE(dens_rate)
    234217  REAL                                                  :: wdens0
    235218  ! IM 080208
    … …  
    239222  REAL, DIMENSION (klon, klev)                          :: deltatw0
    240223  REAL, DIMENSION (klon, klev)                          :: deltaqw0
    241   REAL, DIMENSION (klon, klev)                          :: te, qe
    242 !!  REAL, DIMENSION (klon)                                :: sigmaw1
     224  REAL, DIMENSION (klon, klev)                          :: tb, qb
    243225#ifdef ISO
    244226  REAL, DIMENSION (ntraciso,klon, klev)                          :: deltaxtw0
    245   REAL, DIMENSION (ntraciso,klon, klev)                          :: xte
    246 #endif
    247 
    248   ! Variables liees a la dynamique de population
     227  REAL, DIMENSION (ntraciso,klon, klev)                          :: xtb
     228#endif
     229
     230  ! Variables liees a la dynamique de population 1
    249231  REAL, DIMENSION(klon)                                 :: act
    250232  REAL, DIMENSION(klon)                                 :: rad_wk, tau_wk_inv
    251233  REAL, DIMENSION(klon)                                 :: f_shear
    252234  REAL, DIMENSION(klon)                                 :: drdt
    253   REAL, DIMENSION(klon)                                 :: d_sig_gen, d_sig_death, d_sig_col
     235 
     236  ! Variables liees a la dynamique de population 2
     237  REAL, DIMENSION(klon)                                 :: cont_fact 
     238 
     239  ! Variables liees a la dynamique de population 3
     240  REAL, DIMENSION(klon)                                 :: arad_wk, irad_wk
     241 
     242!!  REAL, DIMENSION(klon)                                 :: d_sig_gen, d_sig_death, d_sig_col
    254243  REAL, DIMENSION(klon)                                 :: wape1_act, wape2_act
    255244  LOGICAL, DIMENSION (klon)                             :: kill_wake
    256   INTEGER, SAVE                                         :: iflag_wk_act
    257   !$OMP THREADPRIVATE(iflag_wk_act)
    258245  REAL                                                  :: drdt_pos
    259246  REAL                                                  :: tau_wk_inv_min
     247  ! Some components of the tendencies of state variables 
     248  REAL, DIMENSION (klon)                                :: d_sig_gen2, d_sig_death2, d_sig_col2, d_sig_spread2, d_sig_bnd2
     249  REAL, DIMENSION (klon)                                :: d_asig_death2, d_asig_aicol2, d_asig_iicol2, d_asig_spread2, d_asig_bnd2
     250  REAL, DIMENSION (klon)                                :: d_dens_gen2, d_dens_death2, d_dens_col2, d_dens_bnd2
     251  REAL, DIMENSION (klon)                                :: d_adens_death2, d_adens_icol2, d_adens_acol2, d_adens_bnd2
    260252
    261253  ! Variables pour les GW
    … …  
    266258
    267259  ! Variables liees au calcul de hw
    268   REAL, DIMENSION (klon)                                :: ptop_provis, ptop, ptop_new
     260  REAL, DIMENSION (klon)                                :: ptop
    269261  REAL, DIMENSION (klon)                                :: sum_dth
    270262  REAL, DIMENSION (klon)                                :: dthmin
    … …  
    278270  ! Sub-timestep tendencies and related variables
    279271  REAL, DIMENSION (klon, klev)                          :: d_deltatw, d_deltaqw
    280   REAL, DIMENSION (klon, klev)                          :: d_te, d_qe
    281   REAL, DIMENSION (klon)                                :: d_awdens, d_wdens
    282   REAL, DIMENSION (klon)                                :: d_sigmaw, alpha
     272  REAL, DIMENSION (klon, klev)                          :: d_tb, d_qb
     273  REAL, DIMENSION (klon)                                :: d_wdens, d_awdens, d_sigmaw, d_asigmaw
     274  REAL, DIMENSION (klon)                                :: d_sig_gen, d_sig_death, d_sig_col, d_sig_spread, d_sig_bnd
     275  REAL, DIMENSION (klon)                                :: d_asig_death, d_asig_aicol, d_asig_iicol, d_asig_spread, d_asig_bnd
     276  REAL, DIMENSION (klon)                                :: d_dens_gen, d_dens_death, d_dens_col, d_dens_bnd
     277  REAL, DIMENSION (klon)                                :: d_adens_death, d_adens_icol, d_adens_acol, d_adens_bnd
     278  REAL, DIMENSION (klon)                                :: agfl              !! gust front length of active wakes
     279                                                                             !!  per unit area
     280  REAL, DIMENSION (klon)                                :: alpha, alpha_tot
    283281  REAL, DIMENSION (klon)                                :: q0_min, q1_min
    284282  LOGICAL, DIMENSION (klon)                             :: wk_adv, ok_qx_qw
    285283#ifdef ISO
    286284  REAL, DIMENSION (ntraciso,klon, klev)                          :: d_deltaxtw
    287   REAL, DIMENSION (ntraciso,klon, klev)                          :: d_xte
    288 #endif
    289   REAL, SAVE                                            :: epsilon
    290   !$OMP THREADPRIVATE(epsilon)
    291   DATA epsilon/1.E-15/
     285  REAL, DIMENSION (ntraciso,klon, klev)                          :: d_xtb
     286#endif
    292287
    293288  ! Autres variables internes
    294   INTEGER                                               ::isubstep, k, i
    295 
     289  INTEGER                                               ::isubstep, k, i, igout
     290
     291  REAL                                                  :: wdensmin
     292
     293  REAL                                                  :: sigmaw_targ
    296294  REAL                                                  :: wdens_targ
    297   REAL                                                  :: sigmaw_targ
    298 
    299   REAL, DIMENSION (klon)                                :: sum_thu, sum_tu, sum_qu, sum_thvu
    300   REAL, DIMENSION (klon)                                :: sum_dq, sum_rho
     295  REAL                                                  :: d_sigmaw_targ
     296  REAL                                                  :: d_wdens_targ
     297
     298  REAL, DIMENSION (klon)                                :: sum_thx, sum_tx, sum_qx, sum_thvx
     299  REAL, DIMENSION (klon)                                :: sum_dq
    301300  REAL, DIMENSION (klon)                                :: sum_dtdwn, sum_dqdwn
    302   REAL, DIMENSION (klon)                                :: av_thu, av_tu, av_qu, av_thvu
    303   REAL, DIMENSION (klon)                                :: av_dth, av_dq, av_rho
     301  REAL, DIMENSION (klon)                                :: av_thx, av_tx, av_qx, av_thvx
     302  REAL, DIMENSION (klon)                                :: av_dth, av_dq
    304303  REAL, DIMENSION (klon)                                :: av_dtdwn, av_dqdwn
    305 ! pas besoin d'isos ici
    306 
    307   REAL, DIMENSION (klon, klev)                          :: rho, rhow
     304
     305  REAL, DIMENSION (klon, klev)                          :: rho
    308306  REAL, DIMENSION (klon, klev+1)                        :: rhoh
    309   REAL, DIMENSION (klon, klev)                          :: rhow_moyen
    310307  REAL, DIMENSION (klon, klev)                          :: zh
    311308  REAL, DIMENSION (klon, klev+1)                        :: zhh
    312   REAL, DIMENSION (klon, klev)                          :: epaisseur1, epaisseur2
    313 
    314   REAL, DIMENSION (klon, klev)                          :: the, thu
     309
     310  REAL, DIMENSION (klon, klev)                          :: thb, thx
    315311
    316312  REAL, DIMENSION (klon, klev)                          :: omgbw
    … …  
    347343  REAL, DIMENSION (klon, klev)                          :: detr
    348344
    349   REAL, DIMENSION(klon)                                 :: sigmaw_in             ! pour les prints
    350   REAL, DIMENSION(klon)                                 :: awdens_in, wdens_in   ! pour les prints
     345  REAL, DIMENSION(klon)                                 :: sigmaw_in, asigmaw_in ! pour les prints
     346  REAL, DIMENSION(klon)                                 :: wdens_in, awdens_in   ! pour les prints
     347
     348!!!  LOGICAL                                               :: phys_sub=.true.
     349  LOGICAL                                               :: phys_sub=.false.
     350
     351  LOGICAL                                               :: first_call=.true.
     352
     353
     354  !!-- variables liees au nouveau calcul de ptop et hw
     355  REAL, DIMENSION (klon, klev)                          :: int_dth
     356  REAL, DIMENSION (klon, klev)                          :: zzz, dzzz
     357  REAL                                                  :: epsil
     358  REAL, DIMENSION (klon)                                :: ptop1
     359  INTEGER, DIMENSION (klon)                             :: ktop1
     360  REAL, DIMENSION (klon)                                :: omega
     361  REAL, DIMENSION (klon)                                :: h_zzz
     362
     363!print*,'WAKE LJYFz'
    351364
    352365  ! -------------------------------------------------------------------------
    353366  ! Initialisations
    354367  ! -------------------------------------------------------------------------
    355 
    356   ! print*, 'wake initialisations'
    357 !#ifdef ISOVERIF
    358 !        write(*,*) 'wake 358: entree'
    359 !#endif
    360 
    361   ! Essais d'initialisation avec sigmaw = 0.02 et hw = 10.
    362   ! -------------------------------------------------------------------------
    363 
    364 !!  DATA wapecut, sigmad, hwmin/5., .02, 10./
    365 !!  DATA wapecut, sigmad, hwmin/1., .02, 10./
    366   DATA sigmad, hwmin/.02, 10./
    367 !!  DATA wdensmin/1.e-12/
    368   DATA wdensmin/1.e-14/
    369   ! cc nrlmd
    370   DATA sigmaw_max/0.4/
    371   DATA dens_rate/0.1/
    372   ! cc
    373   ! Longueur de maille (en m)
    374   ! -------------------------------------------------------------------------
    375 
    376368  ! ALON = 3.e5
    377369  ! alon = 1.E6
    … …  
    380372  f_shear(:) = 1.       ! 0. for strong shear, 1. for weak shear
    381373
    382 
    383374  ! Configuration de coefgw,stark,wdens (22/02/06 by YU Jingmei)
    384375
    385   ! coefgw : Coefficient pour les ondes de gravité
     376  ! coefgw : Coefficient pour les ondes de gravite
    386377  ! stark : Coefficient k dans Cstar=k*sqrt(2*WAPE)
    387   ! wdens : Densité surfacique de poche froide
     378  ! wdens : Densite surfacique de poche froide
    388379  ! -------------------------------------------------------------------------
    389380
    … …  
    398389  ! alpk = 0.5
    399390  ! alpk = 0.05
    400 
    401  if (first) then
    402 
    403   igout = klon/2+1/klon
    404 
    405   crep_upper = 0.9
    406   crep_sol = 1.0
    407 
    408   ! Get wapecut from parameter file
    409   wapecut = 1.
    410   CALL getin_p('wapecut', wapecut)
    411 
    412   ! cc nrlmd Lecture du fichier wake_param.data
    413   stark=0.33
    414   CALL getin_p('stark',stark)
    415   cstart = stark*sqrt(2.*wapecut)
    416 
    417   alpk=0.25
    418   CALL getin_p('alpk',alpk)
    419 !jyg<
    420 !!  wdens_ref=8.E-12
    421 !!  CALL getin_p('wdens_ref',wdens_ref)
    422   wdens_ref(1)=8.E-12
    423   wdens_ref(2)=8.E-12
    424   CALL getin_p('wdens_ref_o',wdens_ref(1))    !wake number per unit area ; ocean
    425   CALL getin_p('wdens_ref_l',wdens_ref(2))    !wake number per unit area ; land
    426 !>jyg
    427391!
    428 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    429 !!!!!!!!!  Population dynamics parameters    !!!!!!!!!!!!!!!!!!!!!!!!!!!!
    430 !------------------------------------------------------------------------
    431 
    432   iflag_wk_pop_dyn = 0
    433   CALL getin_p('iflag_wk_pop_dyn',iflag_wk_pop_dyn) ! switch between wdens prescribed
    434                                                     ! and wdens prognostic
    435   iflag_wk_act = 0
    436   CALL getin_p('iflag_wk_act',iflag_wk_act) ! 0: act(:)=0.
    437                                             ! 1: act(:)=1.
    438                                             ! 2: act(:)=f(Wape)
    439 
    440   rzero = 5000.
    441   CALL getin_p('rzero_wk', rzero)
    442   aa0 = 3.14*rzero*rzero
     392 igout = klon/2+1/klon
    443393!
    444   tau_cv = 4000.
    445   CALL getin_p('tau_cv', tau_cv)
    446 
    447 !------------------------------------------------------------------------
    448 
    449   coefgw=4.
    450   CALL getin_p('coefgw',coefgw)
    451 
    452   WRITE(*,*) 'stark=', stark
    453   WRITE(*,*) 'alpk=', alpk
    454 !jyg<
    455 !!  WRITE(*,*) 'wdens_ref=', wdens_ref
    456   WRITE(*,*) 'wdens_ref_o=', wdens_ref(1)
    457   WRITE(*,*) 'wdens_ref_l=', wdens_ref(2)
    458 !>jyg
    459   WRITE(*,*) 'iflag_wk_pop_dyn=',iflag_wk_pop_dyn
    460   WRITE(*,*) 'iflag_wk_act',iflag_wk_act
    461   WRITE(*,*) 'coefgw=', coefgw
    462 
    463   flag_wk_check_trgl=.false.
    464   CALL getin_p('flag_wk_check_trgl ', flag_wk_check_trgl)
    465   WRITE(*,*) 'flag_wk_check_trgl=', flag_wk_check_trgl
    466   WRITE(*,*) 'flag_wk_check_trgl OBSOLETE. Utilisr iflag_wk_check_trgl plutot'
    467   iflag_wk_check_trgl=0 ; IF (flag_wk_check_trgl) iflag_wk_check_trgl=1
    468   CALL getin_p('iflag_wk_check_trgl ', iflag_wk_check_trgl)
    469   WRITE(*,*) 'iflag_wk_check_trgl=', iflag_wk_check_trgl
    470 
    471   first=.false.
    472  endif
     394!   sub-time-stepping parameters
     395  dtimesub = dtime/wk_nsub
     396!
     397IF (first_call) THEN
     398!!#define IOPHYS_WK
     399#undef IOPHYS_WK
     400#ifdef IOPHYS_WK
     401  IF (phys_sub) THEN
     402    call iophys_ini(dtimesub)
     403  ELSE
     404    call iophys_ini(dtime)
     405  ENDIF
     406#endif
     407  first_call = .false.
     408ENDIF   !(first_call)
    473409
    474410 IF (iflag_wk_pop_dyn == 0) THEN
    … …  
    483419  !>jyg
    484420 ENDIF  ! (iflag_wk_pop_dyn == 0)
     421!
     422 IF (iflag_wk_pop_dyn >=1) THEN
     423   IF (iflag_wk_pop_dyn == 3) THEN
     424     wdensmin = wdensthreshold
     425   ELSE
     426     wdensmin = wdensinit
     427   ENDIF
     428 ENDIF
    485429
    486430  ! print*,'stark',stark
    … …  
    492436  ! -------------------------------------------------------------------------
    493437
    494   delta_t_min = 0.2
     438   delta_t_min = 0.2
    495439
    496440  ! 1. - Save initial values, initialize tendencies, initialize output fields
    … …  
    503447!!      deltatw0(i, k) = deltatw(i, k)
    504448!!      deltaqw0(i, k) = deltaqw(i, k)
    505 !!      te(i, k) = te0(i, k)
    506 !!      qe(i, k) = qe0(i, k)
     449!!      tb(i, k) = tb0(i, k)
     450!!      qb(i, k) = qb0(i, k)
    507451!!      dtls(i, k) = 0.
    508452!!      dqls(i, k) = 0.
    509453!!      d_deltat_gw(i, k) = 0.
    510 !!      d_te(i, k) = 0.
    511 !!      d_qe(i, k) = 0.
     454!!      d_tb(i, k) = 0.
     455!!      d_qb(i, k) = 0.
    512456!!      d_deltatw(i, k) = 0.
    513457!!      d_deltaqw(i, k) = 0.
    … …  
    521465      deltatw0(:,:) = deltatw(:,:)
    522466      deltaqw0(:,:) = deltaqw(:,:)
    523       te(:,:) = te0(:,:)
    524       qe(:,:) = qe0(:,:)
     467      tb(:,:) = tb0(:,:)
     468      qb(:,:) = qb0(:,:)
    525469      dtls(:,:) = 0.
    526470      dqls(:,:) = 0.
    527471      d_deltat_gw(:,:) = 0.
    528       d_te(:,:) = 0.
    529       d_qe(:,:) = 0.
     472      d_tb(:,:) = 0.
     473      d_qb(:,:) = 0.
    530474      d_deltatw(:,:) = 0.
    531475      d_deltaqw(:,:) = 0.
    532476      d_deltatw2(:,:) = 0.
    533       d_deltaqw2(:,:) = 0. 
     477      d_deltaqw2(:,:) = 0.
    534478#ifdef ISO
    535479      deltaxtw0(:,:,:)= deltaxtw(:,:,:)
    536       xte(:,:,:) = xte0(:,:,:)
     480      xtb(:,:,:) = xtb0(:,:,:)
    537481      dxtls(:,:,:) = 0.
    538       d_xte(:,:,:) = 0.
     482      d_xtb(:,:,:) = 0.
    539483      d_deltaxtw(:,:,:) = 0.
    540484      d_deltaxtw2(:,:,:)=0.
    541       xt1(:,:,:) = 0.
    542       xt2(:,:,:)=0.   
    543       ! init non indispensable mais facilite verif iso
    544       q1(:,:)=0.0
    545       q2(:,:)=0.0
    546 #endif
     485#endif
     486
     487      d_sig_gen2(:)   = 0.
     488      d_sig_death2(:) = 0.
     489      d_sig_col2(:)   = 0.
     490      d_sig_spread2(:)= 0.
     491      d_asig_death2(:) = 0.
     492      d_asig_iicol2(:) = 0.
     493      d_asig_aicol2(:) = 0.
     494      d_asig_spread2(:)= 0.
     495      d_asig_bnd2(:) = 0.
     496      d_asigmaw2(:) = 0.
     497!
     498      d_dens_gen2(:)   = 0.
     499      d_dens_death2(:) = 0.
     500      d_dens_col2(:)   = 0.
     501      d_dens_bnd2(:) = 0.
     502      d_wdens2(:) = 0.
     503      d_adens_bnd2(:) = 0.
     504      d_awdens2(:) = 0.
     505      d_adens_death2(:) = 0.
     506      d_adens_icol2(:) = 0.
     507      d_adens_acol2(:) = 0.
    547508
    548509      IF (iflag_wk_act == 0) THEN
    … …  
    555516!!   sigmaw_in(i) = sigmaw(i)
    556517!!  END DO
    557    sigmaw_in(:) = sigmaw(:)
     518   sigmaw_in(:)  = sigmaw(:)
     519   asigmaw_in(:) = asigmaw(:)
    558520!>jyg
    559 
    560   ! sigmaw1=sigmaw
    561   ! IF (sigd_con.GT.sigmaw1) THEN
    562   ! print*, 'sigmaw,sigd_con', sigmaw, sigd_con
    563   ! ENDIF
    564   IF (iflag_wk_pop_dyn >=1) THEN
    565     DO i = 1, klon
    566       wdens_targ = max(wdens(i),wdensmin)
    567       d_wdens2(i) = wdens_targ - wdens(i)
    568       wdens(i) = wdens_targ
    569     END DO
    570   ELSE
    571     DO i = 1, klon
    572       d_awdens2(i) = 0.
    573       d_wdens2(i) = 0.
    574     END DO
    575   ENDIF  ! (iflag_wk_pop_dyn >=1)
    576 !
    577   DO i = 1, klon
    578     ! c      sigmaw(i) = amax1(sigmaw(i),sigd_con(i))
    579 !jyg<
    580 !!    sigmaw(i) = amax1(sigmaw(i), sigmad)
    581 !!    sigmaw(i) = amin1(sigmaw(i), 0.99)
    582     sigmaw_targ = min(max(sigmaw(i), sigmad),0.99)
    583     d_sigmaw2(i) = sigmaw_targ - sigmaw(i)
    584     sigmaw(i) = sigmaw_targ
    585 !>jyg
    586   END DO
    587 
    588521!
    589522  IF (iflag_wk_pop_dyn >= 1) THEN
    … …  
    595528  ENDIF  ! (iflag_wk_pop_dyn >= 1)
    596529
     530
     531  ! sigmaw1=sigmaw
     532  ! IF (sigd_con.GT.sigmaw1) THEN
     533  ! print*, 'sigmaw,sigd_con', sigmaw, sigd_con
     534  ! ENDIF
     535  IF (iflag_wk_pop_dyn >= 1) THEN
     536    DO i = 1, klon
     537      d_dens_gen2(i)   = 0.
     538      d_dens_death2(i) = 0.
     539      d_dens_col2(i)   = 0.
     540      d_awdens2(i) = 0.
     541      IF (wdens(i) < wdensthreshold) THEN
     542  !!      wdens_targ = max(wdens(i),wdensmin)
     543        wdens_targ = max(wdens(i),wdensinit)
     544        d_dens_bnd2(i) = wdens_targ - wdens(i)
     545        d_wdens2(i) = wdens_targ - wdens(i)
     546        wdens(i) = wdens_targ
     547      ELSE
     548        d_dens_bnd2(i) = 0.
     549        d_wdens2(i) = 0.
     550      ENDIF  !! (wdens(i) < wdensthreshold)
     551    END DO
     552    IF (iflag_wk_pop_dyn >= 2) THEN
     553      DO i = 1, klon 
     554        IF (awdens(i) < wdensthreshold) THEN
     555!!          wdens_targ = min(max(awdens(i),wdensmin),wdens(i))
     556            wdens_targ = min(max(awdens(i),wdensinit),wdens(i))
     557            d_adens_bnd2(i) = wdens_targ - awdens(i)
     558            d_awdens2(i) = wdens_targ - awdens(i)
     559            awdens(i) = wdens_targ
     560        ELSE
     561            wdens_targ = min(awdens(i), wdens(i))
     562            d_adens_bnd2(i) = wdens_targ - awdens(i)
     563            d_awdens2(i) = wdens_targ - awdens(i)
     564            awdens(i) = wdens_targ
     565        ENDIF
     566      END DO
     567    ENDIF ! (iflag_wk_pop_dyn >= 2)
     568  ELSE 
     569    DO i = 1, klon
     570      d_awdens2(i) = 0.
     571      d_wdens2(i) = 0.
     572    END DO
     573  ENDIF  ! (iflag_wk_pop_dyn >= 1)
     574!
     575  DO i = 1, klon
     576    sigmaw_targ = min(max(sigmaw(i), sigmad),0.99)
     577    d_sig_bnd2(i) = sigmaw_targ - sigmaw(i)
     578    d_sigmaw2(i) = sigmaw_targ - sigmaw(i)
     579    sigmaw(i) = sigmaw_targ
     580  END DO
     581!
     582  IF (iflag_wk_pop_dyn == 3) THEN
     583     DO i = 1, klon 
     584        IF ((wdens(i)-awdens(i)) <= smallestreal) THEN
     585          sigmaw_targ = sigmaw(i)
     586        ELSE
     587          sigmaw_targ = min(max(asigmaw(i),sigmad),sigmaw(i))
     588        ENDIF
     589        d_asig_bnd2(i) = sigmaw_targ - asigmaw(i)
     590        d_asigmaw2(i) = sigmaw_targ - asigmaw(i)
     591        asigmaw(i) = sigmaw_targ
     592     END DO
     593  ENDIF ! (iflag_wk_pop_dyn == 3)
     594
    597595  wape(:) = 0.
    598596  wape2(:) = 0.
    599597  d_sigmaw(:) = 0.
     598  d_asigmaw(:) = 0.
    600599  ktopw(:) = 0
    601600!
    602601!<jyg
    603602dth(:,:) = 0.
    604 tu(:,:) = 0.
    605 qu(:,:) = 0.
     603tx(:,:) = 0.
     604qx(:,:) = 0.
    606605dtke(:,:) = 0.
    607606dqke(:,:) = 0.
    608 spread(:,:) = 0.
     607wkspread(:,:) = 0.
    609608omgbdth(:,:) = 0.
    610609omg(:,:) = 0.
    … …  
    624623omgbdq(:,:) = 0.
    625624#ifdef ISO
    626 xtu(:,:,:) = 0.
     625xtx(:,:,:) = 0.
    627626dxtke(:,:,:) = 0.
    628627omgbdxt(:,:,:) = 0.
    … …  
    653652    ktop(i) = 0
    654653    kupper(i) = 0
    655     sum_thu(i) = 0.
    656     sum_tu(i) = 0.
    657     sum_qu(i) = 0.
    658     sum_thvu(i) = 0.
     654    sum_thx(i) = 0.
     655    sum_tx(i) = 0.
     656    sum_qx(i) = 0.
     657    sum_thvx(i) = 0.
    659658    sum_dth(i) = 0.
    660659    sum_dq(i) = 0.
    661     sum_rho(i) = 0.
    662660    sum_dtdwn(i) = 0.
    663661    sum_dqdwn(i) = 0.
    664662
    665     av_thu(i) = 0.
    666     av_tu(i) = 0.
    667     av_qu(i) = 0.
    668     av_thvu(i) = 0.
     663    av_thx(i) = 0.
     664    av_tx(i) = 0.
     665    av_qx(i) = 0.
     666    av_thvx(i) = 0.
    669667    av_dth(i) = 0.
    670668    av_dq(i) = 0.
    671     av_rho(i) = 0.
    672669    av_dtdwn(i) = 0.
    673670    av_dqdwn(i) = 0.
    674671    ! pas besoin d'isos ici
    675672  END DO
    676 
    677 
    678 #ifdef ISOVERIF
    679       ! TODO     
    680 #endif
    681673
    682674  ! Distance between wakes
    … …  
    688680  DO k = 1, klev
    689681    DO i = 1, klon
    690       ! write(*,*)'wake 1',i,k,rd,te(i,k)
    691       rho(i, k) = p(i, k)/(rd*te(i,k))
     682      ! write(*,*)'wake 1',i,k,RD,tb(i,k)
     683      rho(i, k) = p(i, k)/(RD*tb(i,k))
    692684      ! write(*,*)'wake 2',rho(i,k)
    693685      IF (k==1) THEN
    694         ! write(*,*)'wake 3',i,k,rd,te(i,k)
    695         rhoh(i, k) = ph(i, k)/(rd*te(i,k))
    696         ! write(*,*)'wake 4',i,k,rd,te(i,k)
     686        ! write(*,*)'wake 3',i,k,rd,tb(i,k)
     687        rhoh(i, k) = ph(i, k)/(RD*tb(i,k))
     688        ! write(*,*)'wake 4',i,k,rd,tb(i,k)
    697689        zhh(i, k) = 0
    698690      ELSE
    699         ! write(*,*)'wake 5',rd,(te(i,k)+te(i,k-1))
    700         rhoh(i, k) = ph(i, k)*2./(rd*(te(i,k)+te(i,k-1)))
     691        ! write(*,*)'wake 5',rd,(tb(i,k)+tb(i,k-1))
     692        rhoh(i, k) = ph(i, k)*2./(RD*(tb(i,k)+tb(i,k-1)))
    701693        ! write(*,*)'wake 6',(-rhoh(i,k)*RG)+zhh(i,k-1)
    702         zhh(i, k) = (ph(i,k)-ph(i,k-1))/(-rhoh(i,k)*rg) + zhh(i, k-1)
     694        zhh(i, k) = (ph(i,k)-ph(i,k-1))/(-rhoh(i,k)*RG) + zhh(i, k-1)
    703695      END IF
    704696      ! write(*,*)'wake 7',ppi(i,k)
    705       the(i, k) = te(i, k)/ppi(i, k)
    706       thu(i, k) = (te(i,k)-deltatw(i,k)*sigmaw(i))/ppi(i, k)
    707       tu(i, k) = te(i, k) - deltatw(i, k)*sigmaw(i)
    708       qu(i, k) = qe(i, k) - deltaqw(i, k)*sigmaw(i)
    709       ! write(*,*)'wake 8',(rd*(te(i,k)+deltatw(i,k)))
    710       rhow(i, k) = p(i, k)/(rd*(te(i,k)+deltatw(i,k)))
     697      thb(i, k) = tb(i, k)/ppi(i, k)
     698      thx(i, k) = (tb(i,k)-deltatw(i,k)*sigmaw(i))/ppi(i, k)
     699      tx(i, k) = tb(i, k) - deltatw(i, k)*sigmaw(i)
     700      qx(i, k) = qb(i, k) - deltaqw(i, k)*sigmaw(i)
     701      ! write(*,*)'wake 8',(RD*(tb(i,k)+deltatw(i,k)))
    711702      dth(i, k) = deltatw(i, k)/ppi(i, k)
    712703#ifdef ISO
    713704        do ixt=1,ntraciso
    714           xtu(ixt,i,k)  =  xte(ixt,i,k) - deltaxtw(ixt,i,k)*sigmaw(i)
     705          xtx(ixt,i,k)  =  xtb(ixt,i,k) - deltaxtw(ixt,i,k)*sigmaw(i)
    715706        enddo !do ixt=1,ntraciso
    716707#ifdef ISOVERIF
    717708        if (iso_eau.gt.0) then
    718709              call iso_verif_egalite(deltaqw(i,k),deltaxtw(iso_eau,i,k),'wake 723a')
    719               call iso_verif_egalite(qu(i,k),xtu(iso_eau,i,k),'wake 723b')
     710              call iso_verif_egalite(qx(i,k),xtx(iso_eau,i,k),'wake 723b')
    720711        endif
    721712        if (iso_HDO.gt.0) then
    722                if (iso_verif_aberrant_enc_choix_nostop(xtu(iso_hdo,i,k),qu(i,k),ridicule,deltalim, &
    723                         'wake 723c xtu').eq.1) then
     713               if (iso_verif_aberrant_enc_choix_nostop(xtx(iso_hdo,i,k),qx(i,k),ridicule,deltalim, &
     714                        'wake 723c xtx').eq.1) then
    724715                     stop
    725716               endif
    … …  
    730721  END DO
    731722
    732 
    733 
    734723  DO k = 1, klev - 1
    735724    DO i = 1, klon
    … …  
    737726        n2(i, k) = 0
    738727      ELSE
    739         n2(i, k) = amax1(0., -rg**2/the(i,k)*rho(i,k)*(the(i,k+1)-the(i,k-1))/ &
     728        n2(i, k) = amax1(0., -RG**2/thb(i,k)*rho(i,k)*(thb(i,k+1)-thb(i,k-1))/ &
    740729                             (p(i,k+1)-p(i,k-1)))
    741730      END IF
    … …  
    754743  END DO
    755744
    756   ! Calcul de la masse volumique moyenne de la colonne   (bdlmd)
    757   ! -----------------------------------------------------------------
    758 
    759   DO k = 1, klev
    760     DO i = 1, klon
    761       epaisseur1(i, k) = 0.
    762       epaisseur2(i, k) = 0.
    763     END DO
    764   END DO
    765 
    766   DO i = 1, klon
    767     epaisseur1(i, 1) = -(ph(i,2)-ph(i,1))/(rho(i,1)*rg) + 1.
    768     epaisseur2(i, 1) = -(ph(i,2)-ph(i,1))/(rho(i,1)*rg) + 1.
    769     rhow_moyen(i, 1) = rhow(i, 1)
    770   END DO
    771 
    772   DO k = 2, klev
    773     DO i = 1, klon
    774       epaisseur1(i, k) = -(ph(i,k+1)-ph(i,k))/(rho(i,k)*rg) + 1.
    775       epaisseur2(i, k) = epaisseur2(i, k-1) + epaisseur1(i, k)
    776       rhow_moyen(i, k) = (rhow_moyen(i,k-1)*epaisseur2(i,k-1)+rhow(i,k)* &
    777         epaisseur1(i,k))/epaisseur2(i, k)
    778     END DO
    779   END DO
    780 
    781 
     745 
    782746  ! Choose an integration bound well above wake top
    783747  ! -----------------------------------------------------------------
    784748
    785   ! Pupper = 50000.  ! melting level
    786   ! Pupper = 60000.
    787   ! Pupper = 80000.  ! essais pour case_e
    788   DO i = 1, klon
    789     pupper(i) = 0.6*ph(i, 1)
    790     pupper(i) = max(pupper(i), 45000.)
    791     ! cc        Pupper(i) = 60000.
    792   END DO
    793 
    794 
    795749  ! Determine Wake top pressure (Ptop) from buoyancy integral
    796750  ! --------------------------------------------------------
    797751
    798   ! -1/ Pressure of the level where dth becomes less than delta_t_min.
    799 
    800   DO i = 1, klon
    801     ptop_provis(i) = ph(i, 1)
    802   END DO
    803   DO k = 2, klev
    804     DO i = 1, klon
    805 
    806       ! IM v3JYG; ptop_provis(i).LT. ph(i,1)
    807 
    808       IF (dth(i,k)>-delta_t_min .AND. dth(i,k-1)<-delta_t_min .AND. &
    809           ptop_provis(i)==ph(i,1)) THEN
    810         ptop_provis(i) = ((dth(i,k)+delta_t_min)*p(i,k-1)- &
    811                           (dth(i,k-1)+delta_t_min)*p(i,k))/(dth(i,k)-dth(i,k-1))
    812       END IF
    813     END DO
    814   END DO
    815 
    816   ! -2/ dth integral
    817 
    818   DO i = 1, klon
    819     sum_dth(i) = 0.
    820     dthmin(i) = -delta_t_min
    821     z(i) = 0.
    822   END DO
    823 
    824   DO k = 1, klev
    825     DO i = 1, klon
    826       dz(i) = -(amax1(ph(i,k+1),ptop_provis(i))-ph(i,k))/(rho(i,k)*rg)
    827       IF (dz(i)>0) THEN
    828         z(i) = z(i) + dz(i)
    829         sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i)
    830         dthmin(i) = amin1(dthmin(i), dth(i,k))
    831       END IF
    832     END DO
    833   END DO
    834 
    835   ! -3/ height of triangle with area= sum_dth and base = dthmin
    836 
    837   DO i = 1, klon
    838     hw0(i) = 2.*sum_dth(i)/amin1(dthmin(i), -0.5)
    839     hw0(i) = amax1(hwmin, hw0(i))
    840   END DO
    841 
    842   ! -4/ now, get Ptop
    843 
    844   DO i = 1, klon
    845     z(i) = 0.
    846     ptop(i) = ph(i, 1)
    847   END DO
    848 
    849   DO k = 1, klev
    850     DO i = 1, klon
    851       dz(i) = amin1(-(ph(i,k+1)-ph(i,k))/(rho(i,k)*rg), hw0(i)-z(i))
    852       IF (dz(i)>0) THEN
    853         z(i) = z(i) + dz(i)
    854         ptop(i) = ph(i, k) - rho(i, k)*rg*dz(i)
    855       END IF
    856     END DO
    857   END DO
    858 
    859   IF (prt_level>=10) THEN
    860     PRINT *, 'wake-2, ptop_provis(igout), ptop(igout) ', ptop_provis(igout), ptop(igout)
    861   ENDIF
    862 
    863 
    864   ! -5/ Determination de ktop et kupper
    865 
    866   DO k = klev, 1, -1
    867     DO i = 1, klon
    868       IF (ph(i,k+1)<ptop(i)) ktop(i) = k
    869       IF (ph(i,k+1)<pupper(i)) kupper(i) = k
    870     END DO
    871   END DO
    872 
    873   ! On evite kupper = 1 et kupper = klev
    874   DO i = 1, klon
    875     kupper(i) = max(kupper(i), 2)
    876     kupper(i) = min(kupper(i), klev-1)
    877   END DO
    878 
    879 
    880   ! -6/ Correct ktop and ptop
    881 
    882   DO i = 1, klon
    883     ptop_new(i) = ptop(i)
    884   END DO
    885   DO k = klev, 2, -1
    886     DO i = 1, klon
    887       IF (k<=ktop(i) .AND. ptop_new(i)==ptop(i) .AND. &
    888           dth(i,k)>-delta_t_min .AND. dth(i,k-1)<-delta_t_min) THEN
    889         ptop_new(i) = ((dth(i,k)+delta_t_min)*p(i,k-1)-(dth(i, &
    890           k-1)+delta_t_min)*p(i,k))/(dth(i,k)-dth(i,k-1))
    891       END IF
    892     END DO
    893   END DO
    894 
    895   DO i = 1, klon
    896     ptop(i) = ptop_new(i)
    897   END DO
    898 
    899   DO k = klev, 1, -1
    900     DO i = 1, klon
    901       IF (ph(i,k+1)<ptop(i)) ktop(i) = k
    902     END DO
    903   END DO
    904 
    905   IF (prt_level>=10) THEN
    906     PRINT *, 'wake-3, ktop(igout), kupper(igout) ', ktop(igout), kupper(igout)
     752   Do i=1, klon
     753       wk_adv(i) = .True.
     754   Enddo
     755   Call pkupper (klon, klev, ptop, ph, p, pupper, kupper, &
     756                    dth, hw0, rho, delta_t_min, &
     757                    ktop, wk_adv, h_zzz, ptop1, ktop1)
     758 
     759   !!print'("pkupper APPEL ",7i6)',0,int(ptop/100.),int(ptop1/100.),int(pupper/100.),ktop,ktop1,kupper
     760   
     761   IF (prt_level>=10) THEN
     762     PRINT *, 'wake-3, ktop(igout), kupper(igout) ', ktop(igout), kupper(igout)
    907763  ENDIF
    908764
    … …  
    915771        deltaqw(i, k) = 0.
    916772        d_deltatw2(i,k) = -deltatw0(i,k)
     773        d_deltaqw2(i,k) = -deltaqw0(i,k)
    917774#ifdef ISO
    918775        do ixt=1,ntraciso
    … …  
    939796  ! --------------------------------------
    940797
    941   ! Initialize sum_thvu to 1st level virt. pot. temp.
     798  ! Initialize sum_thvx to 1st level virt. pot. temp.
    942799
    943800  DO i = 1, klon
    944801    z(i) = 1.
    945802    dz(i) = 1.
    946     sum_thvu(i) = thu(i, 1)*(1.+epsim1*qu(i,1))*dz(i)
     803    sum_thvx(i) = thx(i, 1)*(1.+epsim1*qx(i,1))*dz(i)
    947804    sum_dth(i) = 0.
    948805  END DO
    … …  
    950807  DO k = 1, klev
    951808    DO i = 1, klon
    952       dz(i) = -(amax1(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*rg)
     809      dz(i) = -(amax1(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*RG)
    953810      IF (dz(i)>0) THEN
     811              ! LJYF : ecriture pas sympa avec un tableau z(i) qui n'est pas utilise come tableau
    954812        z(i) = z(i) + dz(i)
    955         sum_thu(i) = sum_thu(i) + thu(i, k)*dz(i)
    956         sum_tu(i) = sum_tu(i) + tu(i, k)*dz(i)
    957         sum_qu(i) = sum_qu(i) + qu(i, k)*dz(i)
    958         sum_thvu(i) = sum_thvu(i) + thu(i, k)*(1.+epsim1*qu(i,k))*dz(i)
     813        sum_thx(i) = sum_thx(i) + thx(i, k)*dz(i)
     814        sum_tx(i) = sum_tx(i) + tx(i, k)*dz(i)
     815        sum_qx(i) = sum_qx(i) + qx(i, k)*dz(i)
     816        sum_thvx(i) = sum_thvx(i) + thx(i, k)*(1.+epsim1*qx(i,k))*dz(i)
    959817        sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i)
    960818        sum_dq(i) = sum_dq(i) + deltaqw(i, k)*dz(i)
    961         sum_rho(i) = sum_rho(i) + rhow(i, k)*dz(i)
    962819        sum_dtdwn(i) = sum_dtdwn(i) + dtdwn(i, k)*dz(i)
    963820        sum_dqdwn(i) = sum_dqdwn(i) + dqdwn(i, k)*dz(i)
    … …  
    979836
    980837  DO i = 1, klon
    981     av_thu(i) = sum_thu(i)/hw0(i)
    982     av_tu(i) = sum_tu(i)/hw0(i)
    983     av_qu(i) = sum_qu(i)/hw0(i)
    984     av_thvu(i) = sum_thvu(i)/hw0(i)
     838    av_thx(i) = sum_thx(i)/hw0(i)
     839    av_tx(i) = sum_tx(i)/hw0(i)
     840    av_qx(i) = sum_qx(i)/hw0(i)
     841    av_thvx(i) = sum_thvx(i)/hw0(i)
    985842    ! av_thve = sum_thve/hw0
    986843    av_dth(i) = sum_dth(i)/hw0(i)
    987844    av_dq(i) = sum_dq(i)/hw0(i)
    988     av_rho(i) = sum_rho(i)/hw0(i)
    989845    av_dtdwn(i) = sum_dtdwn(i)/hw0(i)
    990846    av_dqdwn(i) = sum_dqdwn(i)/hw0(i)
    991847
    992     wape(i) = -rg*hw0(i)*(av_dth(i)+ &
    993         epsim1*(av_thu(i)*av_dq(i)+av_dth(i)*av_qu(i)+av_dth(i)*av_dq(i)))/av_thvu(i)
     848    wape(i) = -RG*hw0(i)*(av_dth(i)+ &
     849        epsim1*(av_thx(i)*av_dq(i)+av_dth(i)*av_qx(i)+av_dth(i)*av_dq(i)))/av_thvx(i)
    994850
    995851  END DO
     852#ifdef IOPHYS_WK
     853  IF (.not.phys_sub) CALL iophys_ecrit('wape_a',1,'wape_a','J/kg',wape)
     854#endif
    996855
    997856  ! 2.2 Prognostic variable update
    … …  
    1020879  DO i = 1, klon
    1021880    IF (wape(i)<0.) THEN
     881!!      sigmaw(i) = amax1(sigmad, sigd_con(i))
     882      sigmaw_targ = max(sigmad, sigd_con(i))
     883      d_sig_bnd2(i) = d_sig_bnd2(i) + sigmaw_targ - sigmaw(i)
     884      d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i)
     885      sigmaw(i) = sigmaw_targ
     886    ENDIF  !!  (wape(i)<0.)
     887  ENDDO
     888  !
     889  IF (iflag_wk_pop_dyn == 3) THEN
     890    DO i = 1, klon
     891      IF (wape(i)<0.) THEN
     892        sigmaw_targ = max(sigmad, sigd_con(i))
     893        d_asig_bnd2(i) = d_asig_bnd2(i) + sigmaw_targ - asigmaw(i)
     894        d_asigmaw2(i) = d_asigmaw2(i) + sigmaw_targ - asigmaw(i)
     895        asigmaw(i) = sigmaw_targ
     896      ENDIF  !!  (wape(i)<0.)
     897    ENDDO
     898  ENDIF  !!  (iflag_wk_pop_dyn == 3)
     899
     900  DO i = 1, klon
     901    IF (wape(i)<0.) THEN
    1022902      wape(i) = 0.
    1023903      cstar(i) = 0.
    1024904      hw(i) = hwmin
    1025 !jyg<
    1026 !!      sigmaw(i) = amax1(sigmad, sigd_con(i))
    1027       sigmaw_targ = max(sigmad, sigd_con(i))
    1028       d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i)
    1029       sigmaw(i) = sigmaw_targ
    1030 !>jyg
    1031905      fip(i) = 0.
    1032906      gwake(i) = .FALSE.
    … …  
    1037911    END IF
    1038912  END DO
    1039 
     913  !
    1040914
    1041915  ! Check qx and qw positivity
    1042916  ! --------------------------
    1043917  DO i = 1, klon
    1044     q0_min(i) = min((qe(i,1)-sigmaw(i)*deltaqw(i,1)),  &
    1045                     (qe(i,1)+(1.-sigmaw(i))*deltaqw(i,1)))
     918    q0_min(i) = min((qb(i,1)-sigmaw(i)*deltaqw(i,1)),  &
     919                    (qb(i,1)+(1.-sigmaw(i))*deltaqw(i,1)))
    1046920  END DO
    1047921  DO k = 2, klev
    1048922    DO i = 1, klon
    1049       q1_min(i) = min((qe(i,k)-sigmaw(i)*deltaqw(i,k)), &
    1050                       (qe(i,k)+(1.-sigmaw(i))*deltaqw(i,k)))
     923      q1_min(i) = min((qb(i,k)-sigmaw(i)*deltaqw(i,k)), &
     924                      (qb(i,k)+(1.-sigmaw(i))*deltaqw(i,k)))
    1051925      IF (q1_min(i)<=q0_min(i)) THEN
    1052926        q0_min(i) = q1_min(i)
    … …  
    1058932    ok_qx_qw(i) = q0_min(i) >= 0.
    1059933    alpha(i) = 1.
     934    alpha_tot(i) = 1.
    1060935  END DO
    1061936
    … …  
    1070945  ! -----------------
    1071946
    1072   nsub = 10
    1073   dtimesub = dtime/nsub
    1074 
    1075   ! ------------------------------------------------------------
    1076   DO isubstep = 1, nsub
    1077     ! ------------------------------------------------------------
    1078 
     947!    wk_nsub and dtimesub definitions moved to begining of routine.
     948!!  wk_nsub = 10
     949!!  dtimesub = dtime/wk_nsub
     950
     951 
     952  ! ------------------------------------------------------------------------
     953  ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
     954  ! ------------------------------------------------------------------------
     955  !
     956  DO isubstep = 1, wk_nsub
     957  !
     958  ! ------------------------------------------------------------------------
    1079959    ! wk_adv is the logical flag enabling wake evolution in the time advance
    1080960    ! loop
    … …  
    1085965      PRINT *, 'wake-4.1, isubstep,wk_adv(igout),cstar(igout),wape(igout), ptop(igout) ', &
    1086966                          isubstep,wk_adv(igout),cstar(igout),wape(igout), ptop(igout)
     967     
    1087968    ENDIF
    1088969
    1089970    ! cc nrlmd   Ajout d'un recalcul de wdens dans le cas d'un entrainement
    1090     ! négatif de ktop à kupper --------
    1091     ! cc           On calcule pour cela une densité wdens0 pour laquelle on
     971    ! negatif de ktop a kupper --------
     972    ! cc           On calcule pour cela une densite wdens0 pour laquelle on
    1092973    ! aurait un entrainement nul ---
    1093974    !jyg<
    … …  
    1096977    ! des descentes unsaturees. Nous faisons alors l'hypothese que la
    1097978    ! convection profonde cree directement de nouvelles poches, sans passer
    1098     ! par les thermiques. La nouvelle valeur de wdens est alors imposée.
     979    ! par les thermiques. La nouvelle valeur de wdens est alors imposee.
    1099980
    1100981    DO i = 1, klon
    … …  
    1102983      ! c     $           isubstep,wk_adv(i),cstar(i),wape(i)
    1103984      IF (wk_adv(i) .AND. cstar(i)>0.01) THEN
    1104         omg(i, kupper(i)+1) = -rg*amdwn(i, kupper(i)+1)/sigmaw(i) + &
    1105                                rg*amup(i, kupper(i)+1)/(1.-sigmaw(i))
     985        IF ( iflag_wk_profile == 0 ) THEN
     986           omg(i, kupper(i)+1)=-RG*amdwn(i, kupper(i)+1)/sigmaw(i) + &
     987                               RG*amup(i, kupper(i)+1)/(1.-sigmaw(i))
     988        ELSE
     989           omg(i, kupper(i)+1)=0.
     990        ENDIF
    1106991        wdens0 = (sigmaw(i)/(4.*3.14))* &
    1107992          ((1.-sigmaw(i))*omg(i,kupper(i)+1)/((ph(i,1)-pupper(i))*cstar(i)))**(2)
    … …  
    1112997        IF (wdens(i)<=wdens0*1.1) THEN
    1113998          IF (iflag_wk_pop_dyn >= 1) THEN
     999             d_dens_bnd2(i) = d_dens_bnd2(i) + wdens0 - wdens(i)
    11141000             d_wdens2(i) = d_wdens2(i) + wdens0 - wdens(i)
    11151001          ENDIF
    … …  
    11191005    END DO
    11201006
     1007    IF (iflag_wk_pop_dyn == 0 .AND. ok_bug_gfl) THEN
     1008!!--------------------------------------------------------
     1009!!Bug : computing gfl and rad_wk before changing sigmaw
     1010!!      This bug exists only for iflag_wk_pop_dyn=0. Otherwise, gfl and rad_wk
     1011!!      are computed within  wake_popdyn
     1012!!--------------------------------------------------------
     1013      DO i = 1, klon
     1014        IF (wk_adv(i)) THEN
     1015          gfl(i) = 2.*sqrt(3.14*wdens(i)*sigmaw(i))
     1016          rad_wk(i) = sqrt(sigmaw(i)/(3.14*wdens(i)))
     1017        END IF
     1018      END DO
     1019    ENDIF   ! (iflag_wk_pop_dyn == 0 .AND. ok_bug_gfl)
     1020!!--------------------------------------------------------
     1021
    11211022    DO i = 1, klon
    11221023      IF (wk_adv(i)) THEN
    1123         gfl(i) = 2.*sqrt(3.14*wdens(i)*sigmaw(i))
    1124         rad_wk(i) = sqrt(sigmaw(i)/(3.14*wdens(i)))
    1125 !jyg<
    1126 !!        sigmaw(i) = amin1(sigmaw(i), sigmaw_max)
    11271024        sigmaw_targ = min(sigmaw(i), sigmaw_max)
    1128         d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i)
     1025        d_sig_bnd2(i) = d_sig_bnd2(i) + sigmaw_targ - sigmaw(i)
     1026        d_sigmaw2(i)  = d_sigmaw2(i)  + sigmaw_targ - sigmaw(i)
    11291027        sigmaw(i) = sigmaw_targ
    1130 !>jyg
    1131       END IF
    1132     END DO
     1028      END IF
     1029    END DO
     1030
     1031    IF (iflag_wk_pop_dyn == 0 .AND. .NOT.ok_bug_gfl) THEN
     1032!!--------------------------------------------------------
     1033!!Fix : computing gfl and rad_wk after changing sigmaw
     1034!!--------------------------------------------------------
     1035      DO i = 1, klon
     1036        IF (wk_adv(i)) THEN
     1037          gfl(i) = 2.*sqrt(3.14*wdens(i)*sigmaw(i))
     1038          rad_wk(i) = sqrt(sigmaw(i)/(3.14*wdens(i)))
     1039        END IF
     1040      END DO
     1041    ENDIF   ! (iflag_wk_pop_dyn == 0 .AND. .NOT.ok_bug_gfl)
     1042!!--------------------------------------------------------
    11331043
    11341044    IF (iflag_wk_pop_dyn >= 1) THEN
    1135 !    The variable "death_rate" is significant only when iflag_wk_pop_dyn = 0.
    1136 !    Here, it has to be set to zero.
     1045  !  The variable "death_rate" is significant only when iflag_wk_pop_dyn = 0.
     1046  !  Here, it has to be set to zero.
    11371047      death_rate(:) = 0.
    1138 
    1139       IF (iflag_wk_act ==2) THEN
     1048    ENDIF
     1049 
     1050    IF (iflag_wk_pop_dyn >= 3) THEN
    11401051      DO i = 1, klon
    11411052        IF (wk_adv(i)) THEN
    1142           wape1_act(i) = abs(cin(i))
    1143           wape2_act(i) = 2.*wape1_act(i) + 1.
    1144           act(i) = min(1., max(0., (wape(i)-wape1_act(i)) / (wape2_act(i)-wape1_act(i)) ))
    1145         ENDIF  ! (wk_adv(i))
    1146       ENDDO
    1147       ENDIF  ! (iflag_wk_act ==2)
    1148 
    1149 
    1150       DO i = 1, klon
    1151         IF (wk_adv(i)) THEN
    1152 !!          tau_wk(i) = max(rad_wk(i)/(3.*cstar(i))*((cstar(i)/cstart)**1.5 - 1), 100.)
    1153           tau_wk_inv(i) = max( (3.*cstar(i))/(rad_wk(i)*((cstar(i)/cstart)**1.5 - 1)), 0.)
    1154           tau_wk_inv_min = min(tau_wk_inv(i), 1./dtimesub)
    1155           drdt(i) = (cstar(i) - wgen(i)*(sigmaw(i)/wdens(i)-aa0)/gfl(i)) / &
    1156                     (1 + 2*f_shear(i)*(2.*sigmaw(i)-aa0*wdens(i)) - 2.*sigmaw(i))
    1157 !!                    (1 - 2*sigmaw(i)*(1.-f_shear(i)))
    1158           drdt_pos=max(drdt(i),0.)
    1159 
    1160 !!          d_wdens(i) = ( wgen(i)*(1.+2.*(sigmaw(i)-sigmad)) &
    1161 !!                     - wdens(i)*tau_wk_inv_min &
    1162 !!                     - 2.*gfl(i)*wdens(i)*Cstar(i) )*dtimesub
    1163           d_awdens(i) = ( wgen(i) - (1./tau_cv)*(awdens(i) - act(i)*wdens(i)) )*dtimesub
    1164           d_wdens(i) = ( wgen(i) - (wdens(i)-awdens(i))*tau_wk_inv_min -  &
    1165                          2.*wdens(i)*gfl(i)*drdt_pos )*dtimesub
    1166           d_wdens(i) = max(d_wdens(i), wdensmin-wdens(i))
    1167 
    1168 !!          d_sigmaw(i) = ( (1.-2*f_shear(i)*sigmaw(i))*(gfl(i)*Cstar(i)+wgen(i)*sigmad/wdens(i)) &
    1169 !!                      + 2.*f_shear(i)*wgen(i)*sigmaw(i)**2/wdens(i) &
    1170 !!                      - sigmaw(i)*tau_wk_inv_min )*dtimesub
    1171           d_sig_gen(i) = wgen(i)*aa0
    1172           d_sig_death(i) = - sigmaw(i)*(1.-awdens(i)/wdens(i))*tau_wk_inv_min
    1173 !!          d_sig_col(i) = - 2*f_shear(i)*sigmaw(i)*gfl(i)*drdt_pos
    1174           d_sig_col(i) = - 2*f_shear(i)*(2.*sigmaw(i)-wdens(i)*aa0)*gfl(i)*drdt_pos
    1175           d_sigmaw(i) = ( d_sig_gen(i) + d_sig_death(i) + d_sig_col(i) + gfl(i)*cstar(i) )*dtimesub
    1176           d_sigmaw(i) = max(d_sigmaw(i), sigmad-sigmaw(i))
     1053         sigmaw_targ = min(asigmaw(i), sigmaw_max)
     1054         d_asig_bnd2(i) = d_asig_bnd2(i) + sigmaw_targ - asigmaw(i)
     1055         d_asigmaw2(i)  = d_asigmaw2(i)  + sigmaw_targ - asigmaw(i)
     1056         asigmaw(i) = sigmaw_targ
    11771057        ENDIF
    11781058      ENDDO
    1179 
    1180       IF (prt_level >= 10) THEN
    1181         print *,'wake, cstar(1), cstar(1)/cstart, rad_wk(1), tau_wk_inv(1), drdt(1) ', &
    1182                        cstar(1), cstar(1)/cstart, rad_wk(1), tau_wk_inv(1), drdt(1)
    1183         print *,'wake, wdens(1), awdens(1), act(1), d_awdens(1) ', &
    1184                        wdens(1), awdens(1), act(1), d_awdens(1)
    1185         print *,'wake, wgen, -(wdens-awdens)*tau_wk_inv, -2.*wdens*gfl*drdt_pos, d_wdens ', &
    1186                        wgen(1), -(wdens(1)-awdens(1))*tau_wk_inv(1), -2.*wdens(1)*gfl(1)*drdt_pos, d_wdens(1)
    1187         print *,'wake, d_sig_gen(1), d_sig_death(1), d_sig_col(1), d_sigmaw(1) ', &
    1188                        d_sig_gen(1), d_sig_death(1), d_sig_col(1), d_sigmaw(1)
    1189       ENDIF
     1059    ENDIF
     1060 
     1061!!--------------------------------------------------------
     1062!!--------------------------------------------------------
     1063    IF (iflag_wk_pop_dyn == 1) THEN
     1064  !
     1065     CALL wake_popdyn_1 (klon, klev, dtime, cstar, tau_wk_inv, wgen, wdens, awdens, sigmaw, &
     1066                  wdensmin, &
     1067                  dtimesub, gfl, rad_wk, f_shear, drdt_pos, &
     1068                  d_awdens, d_wdens, d_sigmaw, &
     1069                  iflag_wk_act, wk_adv, cin, wape, &
     1070                  drdt, &
     1071                  d_dens_gen, d_dens_death, d_dens_col, d_dens_bnd, &
     1072                  d_sig_gen, d_sig_death, d_sig_col, d_sig_spread, d_sig_bnd, &
     1073                  d_wdens_targ, d_sigmaw_targ)
     1074                     
    11901075   
    1191     ELSE  !  (iflag_wk_pop_dyn >= 1)
    1192 
     1076!!--------------------------------------------------------
     1077    ELSEIF (iflag_wk_pop_dyn == 2) THEN
     1078  !
     1079     CALL wake_popdyn_2 ( klon, klev, wk_adv, dtimesub, wgen, &
     1080                             wdensmin, &
     1081                             sigmaw, wdens, awdens, &   !! state variables
     1082                             gfl, cstar, cin, wape, rad_wk, &
     1083                             d_sigmaw, d_wdens, d_awdens, &  !! tendencies                             
     1084                             cont_fact, &
     1085                             d_sig_gen, d_sig_death, d_sig_col, d_sig_spread, d_sig_bnd, &
     1086                             d_dens_gen, d_dens_death, d_dens_col, d_dens_bnd, &
     1087                             d_adens_death, d_adens_icol, d_adens_acol, d_adens_bnd )
     1088sigmaw=sigmaw-d_sigmaw
     1089wdens=wdens-d_wdens
     1090awdens=awdens-d_awdens
     1091
     1092!!--------------------------------------------------------
     1093    ELSEIF (iflag_wk_pop_dyn == 3) THEN
     1094#ifdef IOPHYS_WK
     1095    IF (phys_sub) THEN
     1096     CALL iophys_ecrit('ptop',1,'ptop','Pa',ptop)
     1097     CALL iophys_ecrit('sigmaw',1,'sigmaw','',sigmaw)
     1098     CALL iophys_ecrit('asigmaw',1,'asigmaw','',asigmaw)
     1099     CALL iophys_ecrit('wdens',1,'wdens','1/m2',wdens)
     1100     CALL iophys_ecrit('awdens',1,'awdens','1/m2',awdens)
     1101     CALL iophys_ecrit('rad_wk',1,'rad_wk','m',rad_wk)
     1102     CALL iophys_ecrit('arad_wk',1,'arad_wk','m',arad_wk)
     1103     CALL iophys_ecrit('irad_wk',1,'irad_wk','m',irad_wk)
     1104    ENDIF
     1105#endif
     1106  !
     1107     CALL wake_popdyn_3 ( klon, klev, phys_sub, wk_adv, dtimesub, wgen, &
     1108                             wdensmin, &
     1109                             sigmaw, asigmaw, wdens, awdens, &   !! state variables
     1110                             gfl, agfl, cstar, cin, wape, &
     1111                             rad_wk, arad_wk, irad_wk, &
     1112                             d_sigmaw, d_asigmaw, d_wdens, d_awdens, &  !! tendencies                             
     1113                             d_sig_gen, d_sig_death, d_sig_col, d_sig_spread, d_sig_bnd, &
     1114                             d_asig_death, d_asig_aicol, d_asig_iicol, d_asig_spread, d_asig_bnd, &
     1115                             d_dens_gen, d_dens_death, d_dens_col, d_dens_bnd, &
     1116                             d_adens_death, d_adens_icol, d_adens_acol, d_adens_bnd )
     1117sigmaw=sigmaw-d_sigmaw
     1118asigmaw=asigmaw-d_asigmaw
     1119wdens=wdens-d_wdens
     1120awdens=awdens-d_awdens
     1121   
     1122!!--------------------------------------------------------
     1123    ELSEIF (iflag_wk_pop_dyn == 0) THEN
     1124   
    11931125    ! cc nrlmd
    11941126
    11951127      DO i = 1, klon
    11961128        IF (wk_adv(i)) THEN
    1197           ! cc nrlmd          Introduction du taux de mortalité des poches et
     1129
     1130          ! cc nrlmd          Introduction du taux de mortalite des poches et
    11981131          ! test sur sigmaw_max=0.4
    11991132          ! cc         d_sigmaw(i) = gfl(i)*Cstar(i)*dtimesub
    … …  
    12201153      END DO
    12211154
    1222     ENDIF   !  (iflag_wk_pop_dyn >= 1)
    1223 
    1224 
     1155    ENDIF   !  (iflag_wk_pop_dyn == 1) ... ELSEIF (iflag_wk_pop_dyn == 0)
     1156!!--------------------------------------------------------
     1157!!--------------------------------------------------------
     1158
     1159#ifdef IOPHYS_WK
     1160    IF (phys_sub) THEN
     1161     CALL iophys_ecrit('wdensa',1,'wdensa','m',wdens)
     1162     CALL iophys_ecrit('awdensa',1,'awdensa','m',awdens)
     1163     CALL iophys_ecrit('sigmawa',1,'sigmawa','m',sigmaw)
     1164     CALL iophys_ecrit('asigmawa',1,'asigmawa','m',asigmaw)
     1165    ENDIF
     1166#endif
    12251167    ! calcul de la difference de vitesse verticale poche - zone non perturbee
    12261168    ! IM 060208 differences par rapport au code initial; init. a 0 dp_deltomg
    12271169    ! IM 060208 et omg sur les niveaux de 1 a klev+1, alors que avant l'on definit
    1228     ! IM 060208 au niveau k=1..?
     1170    ! IM 060208 au niveau k=1...
    12291171    !JYG 161013 Correction : maintenant omg est dimensionne a klev.
    12301172    DO k = 1, klev
    … …  
    12541196      DO i = 1, klon
    12551197        IF (wk_adv(i) .AND. k<=ktop(i)) THEN
    1256           dz(i) = -(ph(i,k)-ph(i,k-1))/(rho(i,k-1)*rg)
     1198          dz(i) = -(ph(i,k)-ph(i,k-1))/(rho(i,k-1)*RG)
    12571199          z(i) = z(i) + dz(i)
    12581200          dp_deltomg(i, k) = dp_deltomg(i, 1)
    … …  
    12641206    DO i = 1, klon
    12651207      IF (wk_adv(i)) THEN
    1266         dztop(i) = -(ptop(i)-ph(i,ktop(i)))/(rho(i,ktop(i))*rg)
     1208        dztop(i) = -(ptop(i)-ph(i,ktop(i)))/(rho(i,ktop(i))*RG)
    12671209        ztop(i) = z(i) + dztop(i)
    12681210        omgtop(i) = dp_deltomg(i, 1)*ztop(i)
    … …  
    12821224    DO i = 1, klon
    12831225      IF (wk_adv(i)) THEN
    1284         omgtop(i) = -rho(i, ktop(i))*rg*omgtop(i)
    1285         dp_deltomg(i, 1) = omgtop(i)/(ptop(i)-ph(i,1))
     1226        omgtop(i) = -rho(i, ktop(i))*RG*omgtop(i)
     1227!! LJYF        dp_deltomg(i, 1) = omgtop(i)/(ptop(i)-ph(i,1))
     1228        dp_deltomg(i, 1) = omgtop(i)/min(ptop(i)-ph(i,1),-smallestreal)
    12861229      END IF
    12871230    END DO
    … …  
    12901233      DO i = 1, klon
    12911234        IF (wk_adv(i) .AND. k<=ktop(i)) THEN
    1292           omg(i, k) = -rho(i, k)*rg*omg(i, k)
     1235          omg(i, k) = -rho(i, k)*RG*omg(i, k)
    12931236          dp_deltomg(i, k) = dp_deltomg(i, 1)
    12941237        END IF
    … …  
    13001243    DO i = 1, klon
    13011244      IF (wk_adv(i) .AND. kupper(i)>ktop(i)) THEN
    1302         omg(i, kupper(i)+1) = -rg*amdwn(i, kupper(i)+1)/sigmaw(i) + &
    1303           rg*amup(i, kupper(i)+1)/(1.-sigmaw(i))
     1245        IF ( iflag_wk_profile == 0 ) THEN
     1246           omg(i, kupper(i)+1) =-RG*amdwn(i, kupper(i)+1)/sigmaw(i) + &
     1247          RG*amup(i, kupper(i)+1)/(1.-sigmaw(i))
     1248        ELSE
     1249           omg(i, kupper(i)+1) = 0.
     1250        ENDIF
    13041251        dp_deltomg(i, kupper(i)) = (omgtop(i)-omg(i,kupper(i)+1))/ &
    13051252          (ptop(i)-pupper(i))
    … …  
    13081255
    13091256    ! c      DO i=1,klon
    1310     ! c        print*,'Pente entre 0 et kupper (référence)'
     1257    ! c        print*,'Pente entre 0 et kupper (reference)'
    13111258    ! c     $           ,omg(i,kupper(i)+1)/(pupper(i)-ph(i,1))
    13121259    ! c        print*,'Pente entre ktop et kupper'
    … …  
    13841331        IF (wk_adv(i) .AND. k<=kupper(i)+1) THEN
    13851332          dth(i, k) = deltatw(i, k)/ppi(i, k)
    1386           th1(i, k) = the(i, k) - sigmaw(i)*dth(i, k) ! undisturbed area
    1387           th2(i, k) = the(i, k) + (1.-sigmaw(i))*dth(i, k) ! wake
    1388           q1(i, k) = qe(i, k) - sigmaw(i)*deltaqw(i, k) ! undisturbed area
    1389           q2(i, k) = qe(i, k) + (1.-sigmaw(i))*deltaqw(i, k) ! wake
     1333          th1(i, k) = thb(i, k) - sigmaw(i)*dth(i, k) ! undisturbed area
     1334          th2(i, k) = thb(i, k) + (1.-sigmaw(i))*dth(i, k) ! wake
     1335          q1(i, k) = qb(i, k) - sigmaw(i)*deltaqw(i, k) ! undisturbed area
     1336          q2(i, k) = qb(i, k) + (1.-sigmaw(i))*deltaqw(i, k) ! wake
    13901337#ifdef ISO
    13911338          do ixt=1,ntraciso
    1392             xt1(ixt,i,k) = xte(ixt,i,k) -sigmaw(i)     *deltaxtw(ixt,i,k) ! undisturbed area
    1393             xt2(ixt,i,k) = xte(ixt,i,k) +(1.-sigmaw(i))*deltaxtw(ixt,i,k) ! wake
     1339            xt1(ixt,i,k) = xtb(ixt,i,k) -sigmaw(i)     *deltaxtw(ixt,i,k) ! undisturbed area
     1340            xt2(ixt,i,k) = xtb(ixt,i,k) +(1.-sigmaw(i))*deltaxtw(ixt,i,k) ! wake
    13941341          enddo
    13951342#endif
    … …  
    14691416    END DO
    14701417
    1471     IF (prt_level>=10) THEN
    1472       PRINT *, 'wake-4.3, th1(igout,k) ', (k,th1(igout,k), k=1,klev)
    1473       PRINT *, 'wake-4.3, th2(igout,k) ', (k,th2(igout,k), k=1,klev)
    1474       PRINT *, 'wake-4.3, dth(igout,k) ', (k,dth(igout,k), k=1,klev)
    1475       PRINT *, 'wake-4.3, omgbdth(igout,k) ', (k,omgbdth(igout,k), k=1,klev)
     1418!!    IF (prt_level>=10) THEN
     1419    IF (prt_level>=10 .and. wk_adv(igout)) THEN
     1420      PRINT *, 'wake-4.3, th1(igout,k) ', (k,th1(igout,k), k=1,kupper(igout))
     1421      PRINT *, 'wake-4.3, th2(igout,k) ', (k,th2(igout,k), k=1,kupper(igout))
     1422      PRINT *, 'wake-4.3, dth(igout,k) ', (k,dth(igout,k), k=1,kupper(igout))
     1423      PRINT *, 'wake-4.3, omgbdth(igout,k) ', (k,omgbdth(igout,k), k=1,kupper(igout))
    14761424    ENDIF
    14771425
    … …  
    14891437             (1.-alpha_up(i,k))*omgbdth(i,k)- &
    14901438             alpha_up(i,k+1)*omgbdth(i,k+1))*ppi(i, k)
    1491 !           print*,'d_deltatw=', k, d_deltatw(i,k)
     1439!           print*,'d_d,k_ptop_provis(i)eltatw=', k, d_deltatw(i,k)
    14921440
    14931441          d_deltaqw(i, k) = dtimesub/(ph(i,k)-ph(i,k+1))* &
    … …  
    15231471          ! C
    15241472          ! -----------------------------------------------------------------
    1525           d_te(i, k) = dtimesub*((rre1(i)*omg(i,k)*sigmaw(i)*d_th1(i,k)- &
     1473          d_tb(i, k) = dtimesub*((rre1(i)*omg(i,k)*sigmaw(i)*d_th1(i,k)- &
    15261474                                  rre2(i)*omg(i,k+1)*(1.-sigmaw(i))*d_th2(i,k+1))/ &
    15271475                                 (ph(i,k)-ph(i,k+1)) &
    … …  
    15291477                                 (ph(i,k)-ph(i,k+1)) )*ppi(i, k)
    15301478
    1531           d_qe(i, k) = dtimesub*((rre1(i)*omg(i,k)*sigmaw(i)*d_q1(i,k)- &
     1479          d_qb(i, k) = dtimesub*((rre1(i)*omg(i,k)*sigmaw(i)*d_q1(i,k)- &
    15321480                                  rre2(i)*omg(i,k+1)*(1.-sigmaw(i))*d_q2(i,k+1))/ &
    15331481                                 (ph(i,k)-ph(i,k+1)) &
    … …  
    15361484#ifdef ISO
    15371485        do ixt=1,ntraciso
    1538          d_xte(ixt,i,k) =  dtimesub*( &
     1486         d_xtb(ixt,i,k) =  dtimesub*( &
    15391487             ( rre1(i)*omg(i,k  )*sigmaw(i)     *d_xt1(ixt,i,k) &
    15401488              -rre2(i)*omg(i,k+1)*(1.-sigmaw(i))*d_xt2(ixt,i,k+1) ) &
    … …  
    15461494#endif
    15471495        ELSE IF (wk_adv(i) .AND. k==kupper(i)) THEN
    1548           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)
    1549 
    1550           d_qe(i, k) = dtimesub*(rre1(i)*omg(i,k)*sigmaw(i)*d_q1(i,k)/(ph(i,k)-ph(i,k+1)))
     1496          d_tb(i, k) = dtimesub*(rre1(i)*omg(i,k)*sigmaw(i)*d_th1(i,k)/(ph(i,k)-ph(i,k+1)))*ppi(i, k)
     1497
     1498          d_qb(i, k) = dtimesub*(rre1(i)*omg(i,k)*sigmaw(i)*d_q1(i,k)/(ph(i,k)-ph(i,k+1)))
    15511499
    15521500#ifdef ISO
    15531501        do ixt=1,ntraciso
    1554          d_xte(ixt,i,k) =  dtimesub*( &
     1502         d_xtb(ixt,i,k) =  dtimesub*( &
    15551503            ( rre1(i)*omg(i,k  )*sigmaw(i)  *d_xt1(ixt,i,k)     &   
    15561504             /(Ph(i,k)-Ph(i,k+1))) &
    … …  
    16021550
    16031551
    1604           ! Coefficient de répartition
     1552          ! Coefficient de repartition
    16051553
    16061554          crep(i, k) = crep_sol*(ph(i,kupper(i))-ph(i,k))/ &
    16071555            (ph(i,kupper(i))-ph(i,1))
    16081556          crep(i, k) = crep(i, k) + crep_upper*(ph(i,1)-ph(i,k))/ &
    1609             (p(i,1)-ph(i,kupper(i)))
     1557            (ph(i,1)-ph(i,kupper(i)))
    16101558
    16111559
    … …  
    16461594!
    16471595
    1648           ! cc nrlmd          Prise en compte du taux de mortalité
    1649           ! cc               Définitions de entr, detr
     1596          ! cc nrlmd          Prise en compte du taux de mortalite
     1597          ! cc               Definitions de entr, detr
    16501598!jyg<
    16511599!!            detr(i, k) = 0.
    … …  
    16571605                         sigmaw(i)*(1.-sigmaw(i))*dp_deltomg(i, k)   
    16581606!>jyg
    1659             spread(i, k) = (entr(i,k)-detr(i,k))/sigmaw(i)
    1660 
    1661           ! cc        spread(i,k) =
     1607            wkspread(i, k) = (entr(i,k)-detr(i,k))/sigmaw(i)
     1608
     1609          ! cc        wkspread(i,k) =
    16621610          ! (1.-sigmaw(i))*dp_deltomg(i,k)+gfl(i)*Cstar(i)/
    16631611          ! cc     $  sigmaw(i)
    16641612
    16651613
    1666           ! ajout d'un effet onde de gravité -Tgw(k)*deltatw(k) 03/02/06 YU
     1614          ! ajout d'un effet onde de gravite -Tgw(k)*deltatw(k) 03/02/06 YU
    16671615          ! Jingmei
    16681616
    … …  
    16761624          ! Sans GW
    16771625
    1678           ! deltatw(k)=deltatw(k)+dtimesub*(ff+dtKE(k)-spread(k)*deltatw(k))
     1626          ! deltatw(k)=deltatw(k)+dtimesub*(ff+dtKE(k)-wkspread(k)*deltatw(k))
    16791627
    16801628          ! GW formule 1
    16811629
    16821630          ! deltatw(k) = deltatw(k)+dtimesub*
    1683           ! $         (ff+dtKE(k) - spread(k)*deltatw(k)-Tgw(k)*deltatw(k))
     1631          ! $         (ff+dtKE(k) - wkspread(k)*deltatw(k)-Tgw(k)*deltatw(k))
    16841632
    16851633          ! GW formule 2
    … …  
    17411689    ! Scale tendencies so that water vapour remains positive in w and x.
    17421690
    1743     CALL wake_vec_modulation(klon, klev, wk_adv, epsilon, qe, d_qe, deltaqw, &
     1691    CALL wake_vec_modulation(klon, klev, wk_adv, epsilon_loc, qb, d_qb, deltaqw, &
    17441692      d_deltaqw, sigmaw, d_sigmaw, alpha)
     1693    !
     1694    ! Alpha_tot = Product of all the alpha's
     1695    DO i = 1, klon
     1696      IF (wk_adv(i)) THEN
     1697        alpha_tot(i) = alpha_tot(i)*alpha(i)   
     1698      END IF
     1699    END DO
    17451700
    17461701    ! cc nrlmd
    … …  
    17531708      DO i = 1, klon
    17541709        IF (wk_adv(i) .AND. k<=kupper(i)) THEN
    1755           d_te(i, k) = alpha(i)*d_te(i, k)
    1756           d_qe(i, k) = alpha(i)*d_qe(i, k)
     1710          d_tb(i, k) = alpha(i)*d_tb(i, k)
     1711          d_qb(i, k) = alpha(i)*d_qb(i, k)
    17571712          d_deltatw(i, k) = alpha(i)*d_deltatw(i, k)
    17581713          d_deltaqw(i, k) = alpha(i)*d_deltaqw(i, k)
    … …  
    17601715#ifdef ISO
    17611716        do ixt=1,ntraciso
    1762           d_xte(ixt,i,k)=alpha(i)*d_xte(ixt,i,k)
     1717          d_xtb(ixt,i,k)=alpha(i)*d_xtb(ixt,i,k)
    17631718          d_deltaxtw(ixt,i,k)=alpha(i)*d_deltaxtw(ixt,i,k)
    17641719        enddo !do ixt=1,ntraciso
    … …  
    17791734      DO i = 1, klon
    17801735        IF (wk_adv(i) .AND. k<=kupper(i)) THEN
    1781           dtls(i, k) = dtls(i, k) + d_te(i, k)
    1782           dqls(i, k) = dqls(i, k) + d_qe(i, k)
     1736          dtls(i, k) = dtls(i, k) + d_tb(i, k)
     1737          dqls(i, k) = dqls(i, k) + d_qb(i, k)
    17831738#ifdef ISO
    17841739        do ixt=1,ntraciso
    1785           dxtls(ixt,i,k)=dxtls(ixt,i,k)+d_xte(ixt,i,k)
     1740          dxtls(ixt,i,k)=dxtls(ixt,i,k)+d_xtb(ixt,i,k)
    17861741        enddo !do ixt=1,ntraciso
    17871742#endif
    … …  
    18171772      DO i = 1, klon
    18181773        IF (wk_adv(i) .AND. k<=kupper(i)) THEN
    1819           te(i, k) = te0(i, k) + dtls(i, k)
    1820           qe(i, k) = qe0(i, k) + dqls(i, k)
    1821           the(i, k) = te(i, k)/ppi(i, k)
     1774          tb(i, k) = tb0(i, k) + dtls(i, k)
     1775          qb(i, k) = qb0(i, k) + dqls(i, k)
     1776          thb(i, k) = tb(i, k)/ppi(i, k)
    18221777          deltatw(i, k) = deltatw(i, k) + d_deltatw(i, k)
    18231778          deltaqw(i, k) = deltaqw(i, k) + d_deltaqw(i, k)
    18241779          dth(i, k) = deltatw(i, k)/ppi(i, k)
    1825           ! c      print*,'k,qx,qw',k,qe(i,k)-sigmaw(i)*deltaqw(i,k)
    1826           ! c     $        ,qe(i,k)+(1-sigmaw(i))*deltaqw(i,k)
     1780          ! c      print*,'k,qx,qw',k,qb(i,k)-sigmaw(i)*deltaqw(i,k)
     1781          ! c     $        ,qb(i,k)+(1-sigmaw(i))*deltaqw(i,k)
    18271782#ifdef ISO
    18281783        do ixt=1,ntraciso
    1829           xte(ixt,i,k) = xte0(ixt,i,k) + dxtls(ixt,i,k)
     1784          xtb(ixt,i,k) = xtb0(ixt,i,k) + dxtls(ixt,i,k)
    18301785          deltaxtw(ixt,i,k) = deltaxtw(ixt,i,k)+d_deltaxtw(ixt,i,k)
    18311786        enddo !do ixt=1,ntraciso
    … …  
    18481803        DO k= 1,klev
    18491804          DO i = 1,klon
    1850               call iso_verif_egalite_choix(xte(iso_eau,i,k), &
    1851                 qe(i,k),'wake 1379',errmax,errmaxrel)
     1805              call iso_verif_egalite_choix(xtb(iso_eau,i,k), &
     1806                qb(i,k),'wake 1379',errmax,errmaxrel)
    18521807          enddo ! DO i = 1,klon     
    18531808        enddo ! DO k= 1,klev
    … …  
    18551810      if (iso_hdo.gt.0) then
    18561811      call iso_verif_aberrant_enc_vect2D( &
    1857             xte,qe, &
    1858             'wake 1456, xte apres modifs',ntraciso,klon,klev)
     1812            xtb,qb, &
     1813            'wake 1456, xtb apres modifs',ntraciso,klon,klev)
    18591814!      call iso_verif_aberrant_enc_vect2D_ns(
    18601815!     :       deltaxtw,deltaqw,
    … …  
    18681823      DO i = 1, klon
    18691824        IF (wk_adv(i)) THEN
    1870           awdens(i) = awdens(i) + d_awdens(i)
     1825          d_sig_gen2(i)   = d_sig_gen2(i)   + d_sig_gen(i)
     1826          d_sig_death2(i) = d_sig_death2(i) + d_sig_death(i)
     1827          d_sig_col2(i)   = d_sig_col2(i)   + d_sig_col(i)
     1828          d_sig_spread2(i)= d_sig_spread2(i)+ d_sig_spread(i)
     1829          d_sig_bnd2(i)   = d_sig_bnd2(i)   + d_sig_bnd(i)
     1830        END IF
     1831      END DO
     1832!  Bounds
     1833      DO i = 1, klon
     1834        IF (wk_adv(i)) THEN
     1835          sigmaw_targ = max(sigmaw(i),sigmad)
     1836          d_sig_bnd2(i) = d_sig_bnd2(i) + sigmaw_targ - sigmaw(i)
     1837          d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i)
     1838          sigmaw(i) = sigmaw_targ
     1839        END IF
     1840      END DO
     1841!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! wdens  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     1842!  Cumulatives
     1843      DO i = 1, klon
     1844        IF (wk_adv(i)) THEN
    18711845          wdens(i) = wdens(i) + d_wdens(i)
    1872           d_awdens2(i) = d_awdens2(i) + d_awdens(i)
    18731846          d_wdens2(i) = d_wdens2(i) + d_wdens(i)
     1847          d_dens_gen2(i)   = d_dens_gen2(i)   + d_dens_gen(i)
     1848          d_dens_death2(i) = d_dens_death2(i) + d_dens_death(i)
     1849          d_dens_col2(i)   = d_dens_col2(i)   + d_dens_col(i)
     1850          d_dens_bnd2(i)   = d_dens_bnd2(i)   + d_dens_bnd(i)
    18741851        END IF
    18751852      END DO
     1853!  Bounds
    18761854      DO i = 1, klon
    18771855        IF (wk_adv(i)) THEN
    18781856          wdens_targ = max(wdens(i),wdensmin)
     1857          d_dens_bnd2(i) = d_dens_bnd2(i) + wdens_targ - wdens(i)
    18791858          d_wdens2(i) = d_wdens2(i) + wdens_targ - wdens(i)
    18801859          wdens(i) = wdens_targ
    1881 !
     1860        END IF
     1861      END DO
     1862!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! awdens !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     1863!  Cumulatives
     1864      DO i = 1, klon
     1865        IF (wk_adv(i)) THEN
     1866          awdens(i) = awdens(i) + d_awdens(i)
     1867          d_awdens2(i) = d_awdens2(i) + d_awdens(i)
     1868        END IF
     1869      END DO
     1870!  Bounds
     1871      DO i = 1, klon
     1872        IF (wk_adv(i)) THEN
    18821873          wdens_targ = min( max(awdens(i),0.), wdens(i) )
     1874          d_adens_bnd2(i) = d_adens_bnd2(i) + wdens_targ - awdens(i)
    18831875          d_awdens2(i) = d_awdens2(i) + wdens_targ - awdens(i)
    18841876          awdens(i) = wdens_targ
    18851877        END IF
    18861878      END DO
    1887       DO i = 1, klon
    1888         IF (wk_adv(i)) THEN
    1889           sigmaw_targ = max(sigmaw(i),sigmad)
    1890           d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i)
    1891           sigmaw(i) = sigmaw_targ
    1892         END IF
    1893       END DO
     1879!
     1880      IF (iflag_wk_pop_dyn >= 2) THEN
     1881!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! awdens again for iflag_wk_pop_dyn >= 2!!!!!!
     1882!  Cumulatives
     1883        DO i = 1, klon
     1884           IF (wk_adv(i)) THEN
     1885               d_adens_death2(i)   = d_adens_death2(i)   + d_adens_death(i)
     1886               d_adens_icol2(i)   = d_adens_icol2(i)   + d_adens_icol(i)
     1887               d_adens_acol2(i)   = d_adens_acol2(i)   + d_adens_acol(i)
     1888               d_adens_bnd2(i)   = d_adens_bnd2(i)   + d_adens_bnd(i)         
     1889           END IF
     1890        END DO
     1891!  Bounds
     1892        DO i = 1, klon
     1893           IF (wk_adv(i)) THEN
     1894               wdens_targ = min( max(awdens(i),0.), wdens(i) )
     1895               d_adens_bnd2(i) = d_adens_bnd2(i) + wdens_targ - awdens(i)
     1896               awdens(i) = wdens_targ
     1897           END IF
     1898        END DO
     1899!
     1900        IF (iflag_wk_pop_dyn == 3) THEN
     1901!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! asigmaw for iflag_wk_pop_dyn = 3!!!!!!
     1902!  Cumulatives
     1903          DO i = 1, klon
     1904             IF (wk_adv(i)) THEN
     1905                 asigmaw(i) = asigmaw(i) + d_asigmaw(i)
     1906                 d_asigmaw2(i) = d_asigmaw2(i) + d_asigmaw(i)
     1907                 d_asig_death2(i)   = d_asig_death2(i)   + d_asig_death(i)
     1908                 d_asig_spread2(i)  = d_asig_spread2(i)  + d_asig_spread(i)
     1909                 d_asig_iicol2(i)   = d_asig_iicol2(i)   + d_asig_iicol(i)
     1910                 d_asig_aicol2(i)   = d_asig_aicol2(i)   + d_asig_aicol(i)
     1911                 d_asig_bnd2(i)     = d_asig_bnd2(i)     + d_asig_bnd(i)         
     1912             END IF
     1913          END DO
     1914!  Bounds
     1915          DO i = 1, klon
     1916             IF (wk_adv(i)) THEN
     1917   !   asigmaw lower bound set to sigmad/2 in order to allow asigmaw values lower than sigmad.
     1918   !!             sigmaw_targ = min(max(asigmaw(i),sigmad),sigmaw(i))
     1919                sigmaw_targ = min(max(asigmaw(i),sigmad/2.),sigmaw(i))
     1920                d_asig_bnd2(i) = d_asig_bnd2(i) + sigmaw_targ - asigmaw(i)
     1921                d_asigmaw2(i) = d_asigmaw2(i) + sigmaw_targ - asigmaw(i)
     1922                asigmaw(i) = sigmaw_targ
     1923             END IF
     1924          END DO
     1925
     1926#ifdef IOPHYS_WK
     1927    IF (phys_sub) THEN
     1928     CALL iophys_ecrit('wdensb',1,'wdensb','m',wdens)
     1929     CALL iophys_ecrit('awdensb',1,'awdensb','m',awdens)
     1930     CALL iophys_ecrit('sigmawb',1,'sigmawb','m',sigmaw)
     1931     CALL iophys_ecrit('asigmawb',1,'asigmawb','m',asigmaw)
     1932!
     1933     call iophys_ecrit('d_wdens2',1,'d_wdens2','',d_wdens2)
     1934     call iophys_ecrit('d_dens_gen2',1,'d_dens_gen2','',d_dens_gen2)
     1935     call iophys_ecrit('d_dens_death2',1,'d_dens_death2','',d_dens_death2)
     1936     call iophys_ecrit('d_dens_col2',1,'d_dens_col2','',d_dens_col2)
     1937     call iophys_ecrit('d_dens_bnd2',1,'d_dens_bnd2','',d_dens_bnd2)
     1938!
     1939     call iophys_ecrit('d_awdens2',1,'d_awdens2','',d_awdens2)
     1940     call iophys_ecrit('d_adens_death2',1,'d_adens_death2','',d_adens_death2)
     1941     call iophys_ecrit('d_adens_icol2',1,'d_adens_icol2','',d_adens_icol2)
     1942     call iophys_ecrit('d_adens_acol2',1,'d_adens_acol2','',d_adens_acol2)
     1943     call iophys_ecrit('d_adens_bnd2',1,'d_adens_bnd2','',d_adens_bnd2)
     1944!
     1945     CALL iophys_ecrit('d_sigmaw2',1,'d_sigmaw2','',d_sigmaw2)
     1946     CALL iophys_ecrit('d_sig_gen2',1,'d_sig_gen2','m',d_sig_gen2)
     1947     CALL iophys_ecrit('d_sig_spread2',1,'d_sig_spread2','',d_sig_spread2)
     1948     CALL iophys_ecrit('d_sig_col2',1,'d_sig_col2','',d_sig_col2)
     1949     CALL iophys_ecrit('d_sig_death2',1,'d_sig_death2','',d_sig_death2)
     1950     CALL iophys_ecrit('d_sig_bnd2',1,'d_sig_bnd2','',d_sig_bnd2)
     1951!
     1952     CALL iophys_ecrit('d_asigmaw2',1,'d_asigmaw2','',d_asigmaw2)
     1953     CALL iophys_ecrit('d_asig_spread2',1,'d_asig_spread2','m',d_asig_spread2)
     1954     CALL iophys_ecrit('d_asig_aicol2',1,'d_asig_aicol2','m',d_asig_aicol2)
     1955     CALL iophys_ecrit('d_asig_iicol2',1,'d_asig_iicol2','m',d_asig_iicol2)
     1956     CALL iophys_ecrit('d_asig_death2',1,'d_asig_death2','m',d_asig_death2)
     1957     CALL iophys_ecrit('d_asig_bnd2',1,'d_asig_bnd2','m',d_asig_bnd2)
     1958    ENDIF
     1959#endif
     1960        ENDIF ! (iflag_wk_pop_dyn == 3)
     1961      ENDIF ! (iflag_wk_pop_dyn >= 2)
    18941962    ENDIF  ! (iflag_wk_pop_dyn >= 1)
    1895 !>jyg
    1896 
    1897 
    1898     ! Determine Ptop from buoyancy integral
    1899     ! ---------------------------------------
    1900 
    1901     ! -     1/ Pressure of the level where dth changes sign.
    1902 
    1903     DO i = 1, klon
    1904       IF (wk_adv(i)) THEN
    1905         ptop_provis(i) = ph(i, 1)
    1906       END IF
    1907     END DO
    1908 
    1909     DO k = 2, klev
    1910       DO i = 1, klon
    1911         IF (wk_adv(i) .AND. ptop_provis(i)==ph(i,1) .AND. &
    1912             dth(i,k)>-delta_t_min .AND. dth(i,k-1)<-delta_t_min) THEN
    1913           ptop_provis(i) = ((dth(i,k)+delta_t_min)*p(i,k-1) - &
    1914                             (dth(i,k-1)+delta_t_min)*p(i,k))/(dth(i,k)-dth(i,k-1))
    1915         END IF
    1916       END DO
    1917     END DO
    1918 
    1919     ! -     2/ dth integral
    1920 
    1921     DO i = 1, klon
    1922       IF (wk_adv(i)) THEN !!! nrlmd
    1923         sum_dth(i) = 0.
    1924         dthmin(i) = -delta_t_min
    1925         z(i) = 0.
    1926       END IF
    1927     END DO
    1928 
    1929     DO k = 1, klev
    1930       DO i = 1, klon
    1931         IF (wk_adv(i)) THEN
    1932           dz(i) = -(amax1(ph(i,k+1),ptop_provis(i))-ph(i,k))/(rho(i,k)*rg)
    1933           IF (dz(i)>0) THEN
    1934             z(i) = z(i) + dz(i)
    1935             sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i)
    1936             dthmin(i) = amin1(dthmin(i), dth(i,k))
    1937           END IF
    1938         END IF
    1939       END DO
    1940     END DO
    1941 
    1942     ! -     3/ height of triangle with area= sum_dth and base = dthmin
    1943 
    1944     DO i = 1, klon
    1945       IF (wk_adv(i)) THEN
    1946         hw(i) = 2.*sum_dth(i)/amin1(dthmin(i), -0.5)
    1947         hw(i) = amax1(hwmin, hw(i))
    1948       END IF
    1949     END DO
    1950 
    1951     ! -     4/ now, get Ptop
    1952 
    1953     DO i = 1, klon
    1954       IF (wk_adv(i)) THEN !!! nrlmd
    1955         ktop(i) = 0
    1956         z(i) = 0.
    1957       END IF
    1958     END DO
    1959 
    1960     DO k = 1, klev
    1961       DO i = 1, klon
    1962         IF (wk_adv(i)) THEN
    1963           dz(i) = amin1(-(ph(i,k+1)-ph(i,k))/(rho(i,k)*rg), hw(i)-z(i))
    1964           IF (dz(i)>0) THEN
    1965             z(i) = z(i) + dz(i)
    1966             ptop(i) = ph(i, k) - rho(i, k)*rg*dz(i)
    1967             ktop(i) = k
    1968           END IF
    1969         END IF
    1970       END DO
    1971     END DO
    1972 
    1973     ! 4.5/Correct ktop and ptop
    1974 
    1975     DO i = 1, klon
    1976       IF (wk_adv(i)) THEN
    1977         ptop_new(i) = ptop(i)
    1978       END IF
    1979     END DO
    1980 
    1981     DO k = klev, 2, -1
    1982       DO i = 1, klon
    1983         ! IM v3JYG; IF (k .GE. ktop(i)
    1984         IF (wk_adv(i) .AND. k<=ktop(i) .AND. ptop_new(i)==ptop(i) .AND. &
    1985             dth(i,k)>-delta_t_min .AND. dth(i,k-1)<-delta_t_min) THEN
    1986           ptop_new(i) = ((dth(i,k)+delta_t_min)*p(i,k-1) - &
    1987                          (dth(i,k-1)+delta_t_min)*p(i,k))/(dth(i,k)-dth(i,k-1))
    1988         END IF
    1989       END DO
    1990     END DO
    1991 
    1992 
    1993     DO i = 1, klon
    1994       IF (wk_adv(i)) THEN
    1995         ptop(i) = ptop_new(i)
    1996       END IF
    1997     END DO
    1998 
    1999     DO k = klev, 1, -1
    2000       DO i = 1, klon
    2001         IF (wk_adv(i)) THEN !!! nrlmd
    2002           IF (ph(i,k+1)<ptop(i)) ktop(i) = k
    2003         END IF
    2004       END DO
    2005     END DO
     1963
     1964
     1965
     1966   Call pkupper (klon, klev, ptop, ph, p, pupper, kupper, &
     1967                    dth, hw, rho, delta_t_min, &
     1968                    ktop, wk_adv, h_zzz, ptop1, ktop1)
     1969   !! print'("pkupper APPEL ",7i6)',isubstep,int(ptop/100.),int(ptop1/100.),int(pupper/100.),ktop,ktop1,kupper
    20061970
    20071971    ! 5/ Set deltatw & deltaqw to 0 above kupper
    … …  
    20281992    DO i = 1, klon
    20291993      IF (wk_adv(i)) THEN !!! nrlmd
    2030         sum_thu(i) = 0.
    2031         sum_tu(i) = 0.
    2032         sum_qu(i) = 0.
    2033         sum_thvu(i) = 0.
     1994        sum_thx(i) = 0.
     1995        sum_tx(i) = 0.
     1996        sum_qx(i) = 0.
     1997        sum_thvx(i) = 0.
    20341998        sum_dth(i) = 0.
    20351999        sum_dq(i) = 0.
    2036         sum_rho(i) = 0.
    20372000        sum_dtdwn(i) = 0.
    20382001        sum_dqdwn(i) = 0.
    20392002
    2040         av_thu(i) = 0.
    2041         av_tu(i) = 0.
    2042         av_qu(i) = 0.
    2043         av_thvu(i) = 0.
     2003        av_thx(i) = 0.
     2004        av_tx(i) = 0.
     2005        av_qx(i) = 0.
     2006        av_thvx(i) = 0.
    20442007        av_dth(i) = 0.
    20452008        av_dq(i) = 0.
    2046         av_rho(i) = 0.
    20472009        av_dtdwn(i) = 0.
    20482010        av_dqdwn(i) = 0.
    … …  
    20532015    ! --------------------------------------
    20542016
    2055     ! Initialize sum_thvu to 1st level virt. pot. temp.
     2017    ! Initialize sum_thvx to 1st level virt. pot. temp.
    20562018
    20572019    DO i = 1, klon
    … …  
    20592021        z(i) = 1.
    20602022        dz(i) = 1.
    2061         sum_thvu(i) = thu(i, 1)*(1.+epsim1*qu(i,1))*dz(i)
     2023        sum_thvx(i) = thx(i, 1)*(1.+epsim1*qx(i,1))*dz(i)
    20622024        sum_dth(i) = 0.
    20632025      END IF
    … …  
    20672029      DO i = 1, klon
    20682030        IF (wk_adv(i)) THEN !!! nrlmd
    2069           dz(i) = -(max(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*rg)
     2031          dz(i) = -(max(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*RG)
    20702032          IF (dz(i)>0) THEN
    20712033            z(i) = z(i) + dz(i)
    2072             sum_thu(i) = sum_thu(i) + thu(i, k)*dz(i)
    2073             sum_tu(i) = sum_tu(i) + tu(i, k)*dz(i)
    2074             sum_qu(i) = sum_qu(i) + qu(i, k)*dz(i)
    2075             sum_thvu(i) = sum_thvu(i) + thu(i, k)*(1.+epsim1*qu(i,k))*dz(i)
     2034            sum_thx(i) = sum_thx(i) + thx(i, k)*dz(i)
     2035            sum_tx(i) = sum_tx(i) + tx(i, k)*dz(i)
     2036            sum_qx(i) = sum_qx(i) + qx(i, k)*dz(i)
     2037            sum_thvx(i) = sum_thvx(i) + thx(i, k)*(1.+epsim1*qx(i,k))*dz(i)
    20762038            sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i)
    20772039            sum_dq(i) = sum_dq(i) + deltaqw(i, k)*dz(i)
    2078             sum_rho(i) = sum_rho(i) + rhow(i, k)*dz(i)
    20792040            sum_dtdwn(i) = sum_dtdwn(i) + dtdwn(i, k)*dz(i)
    20802041            sum_dqdwn(i) = sum_dqdwn(i) + dqdwn(i, k)*dz(i)
    … …  
    21002061    DO i = 1, klon
    21012062      IF (wk_adv(i)) THEN !!! nrlmd
    2102         av_thu(i) = sum_thu(i)/hw0(i)
    2103         av_tu(i) = sum_tu(i)/hw0(i)
    2104         av_qu(i) = sum_qu(i)/hw0(i)
    2105         av_thvu(i) = sum_thvu(i)/hw0(i)
     2063        av_thx(i) = sum_thx(i)/hw0(i)
     2064        av_tx(i) = sum_tx(i)/hw0(i)
     2065        av_qx(i) = sum_qx(i)/hw0(i)
     2066        av_thvx(i) = sum_thvx(i)/hw0(i)
    21062067        av_dth(i) = sum_dth(i)/hw0(i)
    21072068        av_dq(i) = sum_dq(i)/hw0(i)
    2108         av_rho(i) = sum_rho(i)/hw0(i)
    21092069        av_dtdwn(i) = sum_dtdwn(i)/hw0(i)
    21102070        av_dqdwn(i) = sum_dqdwn(i)/hw0(i)
    21112071
    2112         wape(i) = -rg*hw0(i)*(av_dth(i)+epsim1*(av_thu(i)*av_dq(i) + &
    2113                               av_dth(i)*av_qu(i)+av_dth(i)*av_dq(i)))/av_thvu(i)
    2114       END IF
    2115     END DO
     2072        wape(i) = -RG*hw0(i)*(av_dth(i)+epsim1*(av_thx(i)*av_dq(i) + &
     2073                              av_dth(i)*av_qx(i)+av_dth(i)*av_dq(i)))/av_thvx(i)
     2074      END IF
     2075    END DO
     2076
    21162077
    21172078    ! Filter out bad wakes
    … …  
    21462107!!          sigmaw(i) = max(sigmad, sigd_con(i))
    21472108          sigmaw_targ = max(sigmad, sigd_con(i))
     2109          d_sig_bnd2(i) = d_sig_bnd2(i) + sigmaw_targ - sigmaw(i)
    21482110          d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i)
    21492111          sigmaw(i) = sigmaw_targ
     2112!
     2113          d_asig_bnd2(i) = d_asig_bnd2(i) + sigmaw_targ - asigmaw(i)
     2114          d_asigmaw2(i) = d_asigmaw2(i) + sigmaw_targ - asigmaw(i)
     2115          asigmaw(i) = sigmaw_targ
    21502116!>jyg
    21512117          fip(i) = 0.
    … …  
    21572123      END IF
    21582124    END DO
    2159 
    2160   END DO ! end sub-timestep loop
    2161 
     2125  !
     2126  ! ------------------------------------------------------------------------
     2127  !
     2128  END DO   ! isubstep end sub-timestep loop
     2129  !
     2130  ! ------------------------------------------------------------------------
     2131  ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
     2132  ! ------------------------------------------------------------------------
     2133  !
     2134
     2135#ifdef IOPHYS_WK
     2136    IF (.not.phys_sub) CALL iophys_ecrit('wape_b',1,'wape_b','J/kg',wape)
     2137#endif
    21622138  IF (prt_level>=10) THEN
    21632139    PRINT *, 'wake-5, sigmaw(igout), cstar(igout), wape(igout), ptop(igout) ', &
    … …  
    21792155      ! cc
    21802156      z(i) = 0.
    2181       sum_thu(i) = 0.
    2182       sum_tu(i) = 0.
    2183       sum_qu(i) = 0.
    2184       sum_thvu(i) = 0.
     2157      sum_thx(i) = 0.
     2158      sum_tx(i) = 0.
     2159      sum_qx(i) = 0.
     2160      sum_thvx(i) = 0.
    21852161      sum_dth(i) = 0.
    21862162      sum_half_dth(i) = 0.
    21872163      sum_dq(i) = 0.
    2188       sum_rho(i) = 0.
    21892164      sum_dtdwn(i) = 0.
    21902165      sum_dqdwn(i) = 0.
    21912166
    2192       av_thu(i) = 0.
    2193       av_tu(i) = 0.
    2194       av_qu(i) = 0.
    2195       av_thvu(i) = 0.
     2167      av_thx(i) = 0.
     2168      av_tx(i) = 0.
     2169      av_qx(i) = 0.
     2170      av_thvx(i) = 0.
    21962171      av_dth(i) = 0.
    21972172      av_dq(i) = 0.
    2198       av_rho(i) = 0.
    21992173      av_dtdwn(i) = 0.
    22002174      av_dqdwn(i) = 0.
    … …  
    22112185      IF (ok_qx_qw(i)) THEN
    22122186        ! cc
    2213         rho(i, k) = p(i, k)/(rd*te(i,k))
     2187        rho(i, k) = p(i, k)/(RD*tb(i,k))
    22142188        IF (k==1) THEN
    2215           rhoh(i, k) = ph(i, k)/(rd*te(i,k))
     2189          rhoh(i, k) = ph(i, k)/(RD*tb(i,k))
    22162190          zhh(i, k) = 0
    22172191        ELSE
    2218           rhoh(i, k) = ph(i, k)*2./(rd*(te(i,k)+te(i,k-1)))
    2219           zhh(i, k) = (ph(i,k)-ph(i,k-1))/(-rhoh(i,k)*rg) + zhh(i, k-1)
    2220         END IF
    2221         the(i, k) = te(i, k)/ppi(i, k)
    2222         thu(i, k) = (te(i,k)-deltatw(i,k)*sigmaw(i))/ppi(i, k)
    2223         tu(i, k) = te(i, k) - deltatw(i, k)*sigmaw(i)
    2224         qu(i, k) = qe(i, k) - deltaqw(i, k)*sigmaw(i)
    2225         rhow(i, k) = p(i, k)/(rd*(te(i,k)+deltatw(i,k)))
     2192          rhoh(i, k) = ph(i, k)*2./(RD*(tb(i,k)+tb(i,k-1)))
     2193          zhh(i, k) = (ph(i,k)-ph(i,k-1))/(-rhoh(i,k)*RG) + zhh(i, k-1)
     2194        END IF
     2195        thb(i, k) = tb(i, k)/ppi(i, k)
     2196        thx(i, k) = (tb(i,k)-deltatw(i,k)*sigmaw(i))/ppi(i, k)
     2197        tx(i, k) = tb(i, k) - deltatw(i, k)*sigmaw(i)
     2198        qx(i, k) = qb(i, k) - deltaqw(i, k)*sigmaw(i)
    22262199        dth(i, k) = deltatw(i, k)/ppi(i, k)
    22272200#ifdef ISO
    22282201        do ixt=1,ntraciso
    2229           xtu(ixt,i,k)  =  xte(ixt,i,k) - deltaxtw(ixt,i,k)*sigmaw(i)
     2202          xtx(ixt,i,k)  =  xtb(ixt,i,k) - deltaxtw(ixt,i,k)*sigmaw(i)
    22302203        enddo !do ixt=1,ntraciso
    22312204#endif
    … …  
    22382211      if (iso_hdo.gt.0) then
    22392212      call iso_verif_aberrant_enc_vect2D( &
    2240               xtu,qu, &
     2213              xtx,qx, &
    22412214              'wake 1834, apres modifs',ntraciso,klon,klev)
    22422215      endif   
    … …  
    22472220  ! -----------------------------------------------------------
    22482221
    2249   ! Initialize sum_thvu to 1st level virt. pot. temp.
     2222  ! Initialize sum_thvx to 1st level virt. pot. temp.
    22502223
    22512224  DO i = 1, klon
    … …  
    22562229      dz(i) = 1.
    22572230      dz_half(i) = 1.
    2258       sum_thvu(i) = thu(i, 1)*(1.+epsim1*qu(i,1))*dz(i)
     2231      sum_thvx(i) = thx(i, 1)*(1.+epsim1*qx(i,1))*dz(i)
    22592232      sum_dth(i) = 0.
    22602233    END IF
    … …  
    22662239      IF (ok_qx_qw(i)) THEN
    22672240        ! cc
    2268         dz(i) = -(amax1(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*rg)
    2269         dz_half(i) = -(amax1(ph(i,k+1),0.5*(ptop(i)+ph(i,1)))-ph(i,k))/(rho(i,k)*rg)
     2241        dz(i) = -(amax1(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*RG)
     2242        dz_half(i) = -(amax1(ph(i,k+1),0.5*(ptop(i)+ph(i,1)))-ph(i,k))/(rho(i,k)*RG)
    22702243        IF (dz(i)>0) THEN
    22712244          z(i) = z(i) + dz(i)
    2272           sum_thu(i) = sum_thu(i) + thu(i, k)*dz(i)
    2273           sum_tu(i) = sum_tu(i) + tu(i, k)*dz(i)
    2274           sum_qu(i) = sum_qu(i) + qu(i, k)*dz(i)
    2275           sum_thvu(i) = sum_thvu(i) + thu(i, k)*(1.+epsim1*qu(i,k))*dz(i)
     2245          sum_thx(i) = sum_thx(i) + thx(i, k)*dz(i)
     2246          sum_tx(i) = sum_tx(i) + tx(i, k)*dz(i)
     2247          sum_qx(i) = sum_qx(i) + qx(i, k)*dz(i)
     2248          sum_thvx(i) = sum_thvx(i) + thx(i, k)*(1.+epsim1*qx(i,k))*dz(i)
    22762249          sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i)
    22772250          sum_dq(i) = sum_dq(i) + deltaqw(i, k)*dz(i)
    2278           sum_rho(i) = sum_rho(i) + rhow(i, k)*dz(i)
    22792251          sum_dtdwn(i) = sum_dtdwn(i) + dtdwn(i, k)*dz(i)
    22802252          sum_dqdwn(i) = sum_dqdwn(i) + dqdwn(i, k)*dz(i)
    … …  
    23062278    IF (ok_qx_qw(i)) THEN
    23072279      ! cc
    2308       av_thu(i) = sum_thu(i)/hw0(i)
    2309       av_tu(i) = sum_tu(i)/hw0(i)
    2310       av_qu(i) = sum_qu(i)/hw0(i)
    2311       av_thvu(i) = sum_thvu(i)/hw0(i)
     2280      av_thx(i) = sum_thx(i)/hw0(i)
     2281      av_tx(i) = sum_tx(i)/hw0(i)
     2282      av_qx(i) = sum_qx(i)/hw0(i)
     2283      av_thvx(i) = sum_thvx(i)/hw0(i)
    23122284      av_dth(i) = sum_dth(i)/hw0(i)
    23132285      av_dq(i) = sum_dq(i)/hw0(i)
    2314       av_rho(i) = sum_rho(i)/hw0(i)
    23152286      av_dtdwn(i) = sum_dtdwn(i)/hw0(i)
    23162287      av_dqdwn(i) = sum_dqdwn(i)/hw0(i)
    23172288
    2318       wape2(i) = -rg*hw0(i)*(av_dth(i)+epsim1*(av_thu(i)*av_dq(i) + &
    2319                              av_dth(i)*av_qu(i)+av_dth(i)*av_dq(i)))/av_thvu(i)
     2289      wape2(i) = -RG*hw0(i)*(av_dth(i)+epsim1*(av_thx(i)*av_dq(i) + &
     2290                             av_dth(i)*av_qx(i)+av_dth(i)*av_dq(i)))/av_thvx(i)
    23202291    END IF
    23212292  END DO
    2322 
     2293#ifdef IOPHYS_WK
     2294  IF (.not.phys_sub) CALL iophys_ecrit('wape2_a',1,'wape2_a','J/kg',wape2)
     2295#endif
    23232296
    23242297
    … …  
    23502323    END DO
    23512324  END IF
     2325#ifdef IOPHYS_WK
     2326  IF (.not.phys_sub) CALL iophys_ecrit('wape2_b',1,'wape2_b','J/kg',wape2)
     2327#endif
    23522328
    23532329
    … …  
    23842360!!      sigmaw(i) = amax1(sigmad, sigd_con(i))
    23852361      sigmaw_targ = max(sigmad, sigd_con(i))
     2362      d_sig_bnd2(i) = d_sig_bnd2(i) + sigmaw_targ - sigmaw(i)
    23862363      d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i)
    23872364      sigmaw(i) = sigmaw_targ
     2365!
     2366      d_asig_bnd2(i) = d_asig_bnd2(i) + sigmaw_targ - asigmaw(i)
     2367      d_asigmaw2(i) = d_asigmaw2(i) + sigmaw_targ - asigmaw(i)
     2368      asigmaw(i) = sigmaw_targ
    23882369!>jyg
    23892370        fip(i) = 0.
    … …  
    23942375        gwake(i) = .TRUE.
    23952376      END IF
    2396     END IF
     2377#ifdef IOPHYS_WK
     2378  IF (.not.phys_sub) CALL iophys_ecrit('cstar2',1,'cstar2','J/kg',cstar2)
     2379#endif
     2380    END IF  ! (ok_qx_qw(i))
    23972381  END DO
    23982382
    … …  
    24262410    END IF
    24272411  END DO
    2428 
     2412    IF (iflag_wk_pop_dyn >= 3) THEN
     2413#ifdef IOPHYS_WK
     2414      IF (.not.phys_sub) THEN
     2415       CALL iophys_ecrit('fip',1,'fip','J/kg',fip)
     2416       CALL iophys_ecrit('hw',1,'hw','J/kg',hw)
     2417       CALL iophys_ecrit('ptop',1,'ptop','J/kg',ptop)
     2418       CALL iophys_ecrit('wdens',1,'wdens','J/kg',wdens)
     2419       CALL iophys_ecrit('awdens',1,'awdens','m',awdens)
     2420       CALL iophys_ecrit('sigmaw',1,'sigmaw','m',sigmaw)
     2421       CALL iophys_ecrit('asigmaw',1,'asigmaw','m',asigmaw)
     2422!   
     2423       CALL iophys_ecrit('rad_wk',1,'rad_wk','J/kg',rad_wk)
     2424       CALL iophys_ecrit('arad_wk',1,'arad_wk','J/kg',arad_wk)
     2425       CALL iophys_ecrit('irad_wk',1,'irad_wk','J/kg',irad_wk)
     2426!   
     2427       call iophys_ecrit('d_wdens2',1,'d_wdens2','',d_wdens2)
     2428       call iophys_ecrit('d_dens_gen2',1,'d_dens_gen2','',d_dens_gen2)
     2429       call iophys_ecrit('d_dens_death2',1,'d_dens_death2','',d_dens_death2)
     2430       call iophys_ecrit('d_dens_col2',1,'d_dens_col2','',d_dens_col2)
     2431       call iophys_ecrit('d_dens_bnd2',1,'d_dens_bnd2','',d_dens_bnd2)
     2432!   
     2433       call iophys_ecrit('d_awdens2',1,'d_awdens2','',d_awdens2)
     2434       call iophys_ecrit('d_adens_death2',1,'d_adens_death2','',d_adens_death2)
     2435       call iophys_ecrit('d_adens_icol2',1,'d_adens_icol2','',d_adens_icol2)
     2436       call iophys_ecrit('d_adens_acol2',1,'d_adens_acol2','',d_adens_acol2)
     2437       call iophys_ecrit('d_adens_bnd2',1,'d_adens_bnd2','',d_adens_bnd2)
     2438!   
     2439       CALL iophys_ecrit('d_sigmaw2',1,'d_sigmaw2','',d_sigmaw2)
     2440       CALL iophys_ecrit('d_sig_gen2',1,'d_sig_gen2','m',d_sig_gen2)
     2441       CALL iophys_ecrit('d_sig_spread2',1,'d_sig_spread2','',d_sig_spread2)
     2442       CALL iophys_ecrit('d_sig_col2',1,'d_sig_col2','',d_sig_col2)
     2443       CALL iophys_ecrit('d_sig_death2',1,'d_sig_death2','',d_sig_death2)
     2444       CALL iophys_ecrit('d_sig_bnd2',1,'d_sig_bnd2','',d_sig_bnd2)
     2445!   
     2446       CALL iophys_ecrit('d_asigmaw2',1,'d_asigmaw2','',d_asigmaw2)
     2447       CALL iophys_ecrit('d_asig_spread2',1,'d_asig_spread2','m',d_asig_spread2)
     2448       CALL iophys_ecrit('d_asig_aicol2',1,'d_asig_aicol2','m',d_asig_aicol2)
     2449       CALL iophys_ecrit('d_asig_iicol2',1,'d_asig_iicol2','m',d_asig_iicol2)
     2450       CALL iophys_ecrit('d_asig_death2',1,'d_asig_death2','m',d_asig_death2)
     2451       CALL iophys_ecrit('d_asig_bnd2',1,'d_asig_bnd2','m',d_asig_bnd2)
     2452      ENDIF  ! (.not.phys_sub)
     2453#endif
     2454    ENDIF  ! (iflag_wk_pop_dyn >= 3)
    24292455  ! Limitation de sigmaw
    24302456
    … …  
    24412467    DO i = 1, klon
    24422468      kill_wake(i) = ((wape(i)>=wape2(i)) .AND. (wape2(i)<=wapecut)) .OR. (ktopw(i)<=2) .OR. &
    2443           .NOT. ok_qx_qw(i) .OR. (wdens(i) < 2.*wdensmin)
     2469          .NOT. ok_qx_qw(i) .OR. (wdens(i) < wdensthreshold)
     2470!!          .NOT. ok_qx_qw(i) .OR. (wdens(i) < 2.*wdensmin)
    24442471    ENDDO
    24452472  ELSE  ! (iflag_wk_pop_dyn >= 1)
    … …  
    24812508      wape(i) = 0.
    24822509      cstar(i) = 0.
    2483 !!jyg   Outside subroutine "Wake" hw, wdens and sigmaw are zero when there are no wakes
     2510!!jyg   Outside subroutine "Wake" hw, wdens sigmaw and asigmaw are zero when there are no wakes
    24842511!!      hw(i) = hwmin                       !jyg
    24852512!!      sigmaw(i) = sigmad                  !jyg
    24862513      hw(i) = 0.                            !jyg
    24872514      fip(i) = 0.
     2515!
    24882516!!      sigmaw(i) = 0.                        !jyg
    24892517      sigmaw_targ = 0.
    2490       d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i)
     2518      d_sig_bnd2(i) = d_sig_bnd2(i) + sigmaw_targ - sigmaw(i)
     2519!!      d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i)
     2520      d_sigmaw2(i) = sigmaw_targ - sigmaw_in(i)      ! _in = correction jyg 20220124
    24912521      sigmaw(i) = sigmaw_targ
     2522!
     2523      IF (iflag_wk_pop_dyn >= 3) THEN
     2524        sigmaw_targ = 0.
     2525        d_asig_bnd2(i) = d_asig_bnd2(i) + sigmaw_targ - asigmaw(i)
     2526!!        d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i)
     2527        d_asigmaw2(i) = sigmaw_targ - asigmaw_in(i)      ! _in = correction jyg 20220124
     2528        asigmaw(i) = sigmaw_targ
     2529      ELSE
     2530        asigmaw(i) = 0.
     2531      ENDIF ! (iflag_wk_pop_dyn >= 3)
     2532!
    24922533      IF (iflag_wk_pop_dyn >= 1) THEN
    24932534!!        awdens(i) = 0.
    24942535!!        wdens(i) = 0.
    24952536        wdens_targ = 0.
    2496         d_wdens2(i) = wdens_targ - wdens(i)
     2537        d_dens_bnd2(i) = d_dens_bnd2(i) + wdens_targ - wdens(i)
     2538!!        d_wdens2(i) = wdens_targ - wdens(i)
     2539        d_wdens2(i) = wdens_targ - wdens_in(i)      ! jyg 20220916
    24972540        wdens(i) = wdens_targ
    24982541        wdens_targ = 0.
    2499         d_awdens2(i) = wdens_targ - awdens(i)
     2542!!jyg: bug fix : the d_adens_bnd2 computation must be before the update of awdens.
     2543        IF (iflag_wk_pop_dyn >= 2) THEN
     2544            d_adens_bnd2(i) = d_adens_bnd2(i) + wdens_targ - awdens(i)
     2545        ENDIF ! (iflag_wk_pop_dyn >= 2)
     2546!!        d_awdens2(i) = wdens_targ - awdens(i)
     2547        d_awdens2(i) = wdens_targ - awdens_in(i)    ! jyg 20220916
    25002548        awdens(i) = wdens_targ
     2549!!        IF (iflag_wk_pop_dyn == 2) THEN
     2550!!            d_adens_bnd2(i) = d_adens_bnd2(i) + wdens_targ - awdens(i)
     2551!!        ENDIF ! (iflag_wk_pop_dyn == 2)
    25012552      ENDIF  ! (iflag_wk_pop_dyn >= 1)
    25022553    ELSE  ! (kill_wake(i))
    … …  
    25122563                      wape(igout),wape2(igout),ktopw(igout),OK_qx_qw(igout)
    25132564  ENDIF
     2565#ifdef IOPHYS_WK
     2566  IF (.not.phys_sub) CALL iophys_ecrit('wape_c',1,'wape_c','J/kg',wape)
     2567#endif
    25142568
    25152569
    … …  
    25432597  END DO
    25442598!jyg<
    2545   DO i = 1, klon
    2546     d_sigmaw2(i) = d_sigmaw2(i)/dtime
    2547     d_awdens2(i) = d_awdens2(i)/dtime
    2548     d_wdens2(i) = d_wdens2(i)/dtime
    2549   ENDDO
     2599  IF (iflag_wk_pop_dyn >= 1) THEN
     2600    DO i = 1, klon
     2601        IF (ok_qx_qw(i)) THEN
     2602      d_sig_gen2(i) = d_sig_gen2(i)/dtime
     2603      d_sig_death2(i) = d_sig_death2(i)/dtime
     2604      d_sig_col2(i) = d_sig_col2(i)/dtime
     2605      d_sig_spread2(i) = d_sig_spread2(i)/dtime
     2606      d_sig_bnd2(i) = d_sig_bnd2(i)/dtime
     2607      d_sigmaw2(i) = d_sigmaw2(i)/dtime
     2608!
     2609      d_dens_gen2(i) = d_dens_gen2(i)/dtime
     2610      d_dens_death2(i) = d_dens_death2(i)/dtime
     2611      d_dens_col2(i) = d_dens_col2(i)/dtime
     2612      d_dens_bnd2(i) = d_dens_bnd2(i)/dtime
     2613      d_awdens2(i) = d_awdens2(i)/dtime
     2614      d_wdens2(i) = d_wdens2(i)/dtime
     2615        ENDIF
     2616    ENDDO
     2617    IF (iflag_wk_pop_dyn >= 2) THEN
     2618      DO i = 1, klon
     2619        IF (ok_qx_qw(i)) THEN
     2620        d_adens_death2(i) = d_adens_death2(i)/dtime
     2621        d_adens_icol2(i) = d_adens_icol2(i)/dtime
     2622        d_adens_acol2(i) = d_adens_acol2(i)/dtime
     2623        d_adens_bnd2(i) = d_adens_bnd2(i)/dtime
     2624        ENDIF
     2625      ENDDO
     2626      IF (iflag_wk_pop_dyn == 3) THEN
     2627       DO i = 1, klon
     2628          IF (ok_qx_qw(i)) THEN
     2629        d_asig_death2(i)  = d_asig_death2(i)/dtime
     2630        d_asig_iicol2(i)  = d_asig_iicol2(i)/dtime
     2631        d_asig_aicol2(i)  = d_asig_aicol2(i)/dtime
     2632        d_asig_spread2(i) = d_asig_spread2(i)/dtime
     2633        d_asig_bnd2(i) = d_asig_bnd2(i)/dtime
     2634          ENDIF
     2635       ENDDO
     2636      ENDIF ! (iflag_wk_pop_dyn == 3) 
     2637    ENDIF ! (iflag_wk_pop_dyn >= 2) 
     2638  ENDIF  ! (iflag_wk_pop_dyn >= 1)
     2639 
    25502640!>jyg
    25512641
    2552 
    2553 
    2554   RETURN
     2642 RETURN
    25552643END SUBROUTINE wake
    25562644
    2557 SUBROUTINE wake_vec_modulation(nlon, nl, wk_adv, epsilon, qe, d_qe, deltaqw, &
     2645SUBROUTINE wake_vec_modulation(nlon, nl, wk_adv, epsilon_loc, qb, d_qb, deltaqw, &
    25582646    d_deltaqw, sigmaw, d_sigmaw, alpha)
    25592647  ! ------------------------------------------------------
    2560   ! D\'etermination du coefficient alpha tel que les tendances
     2648  ! Dtermination du coefficient alpha tel que les tendances
    25612649  ! corriges alpha*d_G, pour toutes les grandeurs G, correspondent
    25622650  ! a une humidite positive dans la zone (x) et dans la zone (w).
    … …  
    25652653
    25662654  ! Input
    2567   REAL qe(nlon, nl), d_qe(nlon, nl)
     2655  REAL qb(nlon, nl), d_qb(nlon, nl)
    25682656  REAL deltaqw(nlon, nl), d_deltaqw(nlon, nl)
    25692657  REAL sigmaw(nlon), d_sigmaw(nlon)
    … …  
    25762664  REAL alpha1(nlon)
    25772665  REAL x, a, b, c, discrim
    2578   REAL epsilon
    2579   ! DATA epsilon/1.e-15/
     2666  REAL epsilon_loc
    25802667  INTEGER i,k
    25812668
    … …  
    25922679    DO i = 1, nlon
    25932680      IF (wk_adv(i)) THEN
    2594         x = qe(i, k) + (zeta(i,k)-sigmaw(i))*deltaqw(i, k) + d_qe(i, k) + &
     2681        x = qb(i, k) + (zeta(i,k)-sigmaw(i))*deltaqw(i, k) + d_qb(i, k) + &
    25952682          (zeta(i,k)-sigmaw(i))*d_deltaqw(i, k) - d_sigmaw(i) * &
    25962683          (deltaqw(i,k)+d_deltaqw(i,k))
    25972684        a = -d_sigmaw(i)*d_deltaqw(i, k)
    2598         b = d_qe(i, k) + (zeta(i,k)-sigmaw(i))*d_deltaqw(i, k) - &
     2685        b = d_qb(i, k) + (zeta(i,k)-sigmaw(i))*d_deltaqw(i, k) - &
    25992686          deltaqw(i, k)*d_sigmaw(i)
    2600         c = qe(i, k) + (zeta(i,k)-sigmaw(i))*deltaqw(i, k) + epsilon
     2687        c = qb(i, k) + (zeta(i,k)-sigmaw(i))*deltaqw(i, k) + epsilon_loc
    26012688        discrim = b*b - 4.*a*c
    26022689        ! print*, 'x, a, b, c, discrim', x, a, b, c, discrim
    2603         IF (a+b>=0.) THEN !! Condition suffisante pour la positivité de ovap
     2690        IF (a+b>=0.) THEN !! Condition suffisante pour la positivite de ovap
    26042691          alpha1(i) = 1.
    26052692        ELSE
    … …  
    26272714END SUBROUTINE wake_vec_modulation
    26282715
     2716
     2717
     2718SUBROUTINE pkupper (klon, klev, ptop, ph, p, pupper, kupper, &
     2719                    dth, hw_, rho, delta_t_min_in, &
     2720                    ktop, wk_adv, h_zzz, ptop1, ktop1)
     2721
     2722USE lmdz_wake_ini , ONLY : wk_pupper
     2723USE lmdz_wake_ini , ONLY : RG
     2724USE lmdz_wake_ini , ONLY : hwmin
     2725USE lmdz_wake_ini , ONLY : iflag_wk_new_ptop, wk_delta_t_min, wk_frac_int_delta_t
     2726USE lmdz_wake_ini , ONLY : wk_int_delta_t_min
     2727
     2728IMPLICIT NONE
     2729
     2730INTEGER,                              INTENT(IN) :: klon,klev
     2731REAL,       DIMENSION (klon,klev+1) , INTENT(IN) :: ph, p
     2732REAL,       DIMENSION (klon,klev+1) , INTENT(IN) :: rho
     2733LOGICAL,    DIMENSION (klon)        , INTENT(IN) :: wk_adv
     2734REAL,       DIMENSION (klon,klev+1) , INTENT(IN) :: dth
     2735REAL,                                 INTENT(IN) :: delta_t_min_in
     2736
     2737
     2738REAL,       DIMENSION (klon)  , INTENT(OUT)        :: hw_
     2739REAL,       DIMENSION (klon)  , INTENT(OUT)        :: ptop
     2740INTEGER,    DIMENSION (klon)  , INTENT(OUT)        :: Ktop
     2741REAL,       DIMENSION (klon)  , INTENT(OUT)        :: pupper
     2742INTEGER,    DIMENSION (klon)  , INTENT(OUT)        :: kupper
     2743REAL,       DIMENSION (klon)  , INTENT(OUT)        :: h_zzz       !!
     2744REAL,       DIMENSION (klon)  , INTENT(OUT)        :: Ptop1      !!
     2745INTEGER,    DIMENSION (klon)  , INTENT(OUT)        :: ktop1      !!
     2746
     2747INTEGER :: i,k
     2748
     2749LOGICAL,    DIMENSION (klon)       :: wk_active
     2750REAL                               :: delta_t_min
     2751REAL,     DIMENSION (klon)         :: dthmin
     2752REAL,     DIMENSION (klon)         :: ptop_provis,ptop_new
     2753REAL,     DIMENSION (klon)         :: z, dz
     2754REAL,     DIMENSION (klon)         :: sum_dth
     2755
     2756INTEGER,     DIMENSION (klon)                     :: k_ptop_provis
     2757REAL,     DIMENSION (klon)                     :: zk_ptop_provis
     2758REAL,     DIMENSION (klon)                        :: omega        !!
     2759REAL,     DIMENSION (klon,klev+1)                 :: int_dth      !!
     2760REAL,     DIMENSION (klon,klev+1)                 :: dzz          !!
     2761REAL,     DIMENSION (klon,klev+1)                 :: zzz          !!
     2762REAL,     DIMENSION (klon)                 :: frac_int_dth          !!
     2763REAL                                              :: ddd!!
     2764
     2765
     2766INTEGER, SAVE :: ipas=0
     2767
     2768
     2769
     2770!INTEGER, SAVE :: compte=0
     2771
     2772! LJYF : a priori z, dz sum_dth sont aussi des variables internes
     2773! Les eliminer apres verification convergence numerique
     2774
     2775!compte=compte+1
     2776!print*,'compte=',compte
     2777
     2778    ! Determine Ptop from buoyancy integral
     2779    ! ---------------------------------------
     2780
     2781    ! -     1/ Pressure of the level where dth changes sign.
     2782    !print*,'WAKE LJYF'
     2783
     2784
     2785if (iflag_wk_new_ptop==0) then
     2786    delta_t_min=delta_t_min_in
     2787else
     2788    delta_t_min=wk_delta_t_min
     2789endif
     2790
     2791    DO i = 1, klon
     2792        ptop_provis(i) = ph(i, 1)
     2793        k_ptop_provis(i) = 1
     2794    END DO
     2795
     2796    DO k = 2, klev
     2797      DO i = 1, klon
     2798        IF (wk_adv(i) .AND. ptop_provis(i)==ph(i,1) .AND. &
     2799! LJYF changer :           dth(i,k)>=-delta_t_min .AND. dth(i,k-1)<-delta_t_min) THEN
     2800            dth(i,k)>-delta_t_min .AND. dth(i,k-1)<-delta_t_min) THEN
     2801          ptop_provis(i) = ((dth(i,k)+delta_t_min)*p(i,k-1) - &
     2802                            (dth(i,k-1)+delta_t_min)*p(i,k))/(dth(i,k)-dth(i,k-1))
     2803          k_ptop_provis(i) = k
     2804        END IF
     2805      END DO
     2806    END DO
     2807
     2808
     2809
     2810    ! -     2/ dth integral
     2811
     2812    DO i = 1, klon
     2813      IF (wk_adv(i)) THEN !!! nrlmd
     2814        sum_dth(i) = 0.
     2815        dthmin(i) = -delta_t_min
     2816        z(i) = 0.
     2817      END IF
     2818    END DO
     2819
     2820    DO k = 1, klev
     2821      DO i = 1, klon
     2822        IF (wk_adv(i)) THEN
     2823          dz(i) = -(amax1(ph(i,k+1),ptop_provis(i))-ph(i,k))/(rho(i,k)*RG)
     2824          IF (dz(i)>0) THEN
     2825            z(i) = z(i) + dz(i)
     2826            sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i)
     2827            dthmin(i) = amin1(dthmin(i), dth(i,k))
     2828          END IF
     2829        END IF
     2830      END DO
     2831    END DO
     2832
     2833    ! -     3/ height of triangle with area= sum_dth and base = dthmin
     2834
     2835    DO i = 1, klon
     2836      IF (wk_adv(i)) THEN
     2837        hw_(i) = 2.*sum_dth(i)/amin1(dthmin(i), -0.5)
     2838        hw_(i) = amax1(hwmin, hw_(i))
     2839      END IF
     2840    END DO
     2841
     2842    ! -     4/ now, get Ptop
     2843
     2844    DO i = 1, klon
     2845      IF (wk_adv(i)) THEN !!! nrlmd
     2846        ktop(i) = 0
     2847        z(i) = 0.
     2848      END IF
     2849    END DO
     2850
     2851    DO k = 1, klev
     2852      DO i = 1, klon
     2853        IF (wk_adv(i)) THEN
     2854          dz(i) = amin1(-(ph(i,k+1)-ph(i,k))/(rho(i,k)*RG), hw_(i)-z(i))
     2855          IF (dz(i)>0) THEN
     2856            z(i) = z(i) + dz(i)
     2857            ptop(i) = ph(i, k) - rho(i, k)*RG*dz(i)
     2858            ktop(i) = k
     2859          END IF
     2860        END IF
     2861      END DO
     2862    END DO
     2863
     2864    ! 4.5/Correct ktop and ptop
     2865
     2866    DO i = 1, klon
     2867        ptop_new(i) = ptop(i)
     2868    END DO
     2869
     2870    DO k = klev, 2, -1
     2871      DO i = 1, klon
     2872        ! IM v3JYG; IF (k .GE. ktop(i)
     2873        IF (wk_adv(i) .AND. k<=ktop(i) .AND. ptop_new(i)==ptop(i) .AND. &
     2874! LJYF changer :           dth(i,k)>=-delta_t_min .AND. dth(i,k-1)<-delta_t_min) THEN
     2875            dth(i,k)>=-delta_t_min .AND. dth(i,k-1)<-delta_t_min) THEN
     2876          ptop_new(i) = ((dth(i,k)+delta_t_min)*p(i,k-1) - &
     2877                         (dth(i,k-1)+delta_t_min)*p(i,k))/(dth(i,k)-dth(i,k-1))
     2878        END IF
     2879      END DO
     2880    END DO
     2881
     2882
     2883    DO i = 1, klon
     2884        ptop(i) = ptop_new(i)
     2885    END DO
     2886
     2887    DO k = klev, 1, -1
     2888      DO i = 1, klon
     2889        IF (wk_adv(i)) THEN !!! nrlmd
     2890          IF (ph(i,k+1)<ptop(i)) ktop(i) = k
     2891        END IF
     2892      END DO
     2893    END DO
     2894 
     2895!  IF (prt_level>=10) THEN
     2896!    PRINT *, 'wake-3, ktop(igout), kupper(igout) ', ktop(igout), kupper(igout)
     2897!  ENDIF
     2898
     2899    ! -----------------------------------------------------------------------
     2900    ! nouveau calcul de hw et ptop
     2901    ! -----------------------------------------------------------------------
     2902!if (iflag_wk_new_ptop>0) then
     2903do i=1,klon
     2904   ptop1(i)=ph(i,1)
     2905   ktop1(i)=1
     2906   h_zzz(i)=0.
     2907enddo
     2908   
     2909IF (iflag_wk_new_ptop/=0) THEN
     2910   
     2911    int_dth(1:klon,1:klev+1)=0.
     2912    DO i = 1, klon
     2913       IF (wk_adv(i)) THEN
     2914          int_dth(i,1) = 0.
     2915      END IF
     2916    END DO
     2917   
     2918    if (abs(iflag_wk_new_ptop) == 1 ) then
     2919        DO k = 2, klev+1
     2920           Do i = 1, klon
     2921              IF (wk_adv(i)) THEN
     2922                 if (k<=k_ptop_provis(i)) then
     2923                      ddd=dth(i,k-1)*(ph(i,k-1) - max(ptop_provis(i),ph(i,k)))
     2924                      !ddd=dth(i,k-1)*(ph(i,k-1) - ph(i,k))
     2925                 else
     2926                      ddd=0.
     2927                 endif             
     2928                 int_dth(i,k) = int_dth(i,k-1) + ddd
     2929              !ELSE
     2930              !   int_dth(i,k) = 0.
     2931              END IF
     2932           END DO
     2933        END DO
     2934    else
     2935        k_ptop_provis(:)=klev+1
     2936        dthmin(:)=dth(:,1)
     2937        ! calcul de l'int??grale de dT * dP jusqu'au dernier
     2938        ! niveau avec dT<0. (en s'assurant qu'on a bien un
     2939        ! dT negatif plus bas)
     2940        DO k = 1, klev
     2941           DO i = 1, klon
     2942              dthmin(i)=min(dthmin(i),dth(i,k))
     2943              ddd=dth(i,k)*(ph(i,k)-ph(i,k+1))
     2944              if (dthmin(i)<0.) then
     2945                  if (k>=k_ptop_provis(i)) then
     2946                      ddd=0.
     2947                  else if (dth(i,k)>=0.) then
     2948                      ddd=0.
     2949                      k_ptop_provis(i)=k+1
     2950                  endif
     2951              endif
     2952              int_dth(i,k+1) = int_dth(i,k)+ ddd
     2953           ENDDO
     2954        ENDDO
     2955
     2956        DO i = 1, klon
     2957           if ( k_ptop_provis(i)==klev+1 .or. .not. wk_adv(i)) then
     2958                k_ptop_provis(i)=1
     2959           endif
     2960        ENDDO
     2961    endif ! (abs(iflag_wk_new_ptop) == 1 )
     2962   ! print*, 'xxx, int_dth', (k,int_dth(1,k),k=1,klev)
     2963   ! print*, 'xxx, k_ptop_provis', k_ptop_provis(1)
     2964   
     2965
     2966 
     2967    ! On se limite ?? des poches avec integrale dT * dp < -wk_int_delta_t_min
     2968    do i=1,klon
     2969          if (int_dth(i,k_ptop_provis(i)) > -wk_int_delta_t_min .or. k_ptop_provis(i)==1) then
     2970          !if (1==0) then
     2971             wk_active(i)=.false.
     2972             ptop(i)=ph(i,1)
     2973             ktop(i)=1
     2974             hw_(i)=0.
     2975          else
     2976             wk_active(i)=wk_adv(i)
     2977          endif
     2978    enddo
     2979
     2980    DO i=1,klon
     2981       IF (wk_active(i)) THEN
     2982        frac_int_dth(i)=wk_frac_int_delta_t*int_dth(i,k_ptop_provis(i))
     2983       ENDIF
     2984    ENDDO
     2985    DO k = 1,klev
     2986       DO i =1, klon
     2987!         print*,ipas,'yyy ',k,int_dth(i,k),frac_int_dth(i)
     2988          IF (wk_active(i)) THEN
     2989            IF (int_dth(i,k)>=frac_int_dth(i)) THEN
     2990              ktop1(i) = min(k, k_ptop_provis(i))
     2991              !ktop1(i) = k
     2992              !print*,ipas,'yyy ktop1= ',ktop1
     2993            ENDIF
     2994          ENDIF
     2995       END DO
     2996    END DO
     2997    !print*, 'LAMINE'
     2998   
     2999    DO i = 1, klon
     3000       IF (wk_active(i)) THEN
     3001           !print*, ipas,'xxx1, int_dth(i,ktop1(i)), frac_int_dth(i), int_dth(i,ktop1(i)+1) ',ktop1
     3002           ddd=int_dth(i,ktop1(i)+1)-int_dth(i,ktop1(i))
     3003           if (ddd==0.) then
     3004              omega(i)=0.
     3005           else
     3006              omega(i) = (frac_int_dth(i) - int_dth(i,ktop1(i)))/ddd
     3007           endif
     3008           !! print*,'OMEGA ',omega(i)
     3009       END IF
     3010    END DO
     3011   
     3012    !! print*, 'xxx'
     3013    DO i = 1, klon
     3014       IF (wk_active(i)) THEN
     3015      ! print*, 'xxx, int_dth(i,ktop1(i)), frac_int_dth(i), int_dth(i,ktop1(i)+1) ', &
     3016      !               int_dth(i,ktop1(i)), frac_int_dth(i), int_dth(i,ktop1(i)+1)
     3017      ! print*, 'xxx, omega(i), ph(i,ktop1(i)), ph(i,ktop1(i)+1) ',  &
     3018      !e               omega(i), ph(i,ktop1(i)), ph(i,ktop1(i)+1)
     3019          ptop1(i) = min((1 - omega(i))*ph(i,ktop1(i)) + omega(i)*ph(i,ktop1(i)+1), ph(i,1))
     3020      END IF
     3021    END DO
     3022   
     3023    DO i=1, klon
     3024       IF (wk_active(i)) THEN
     3025           zzz(i, 1) = 0
     3026       END IF
     3027     END DO
     3028    DO k = 1, klev
     3029       DO i = 1, klon
     3030           IF (wk_active(i)) THEN         
     3031              dzz(i,k) = (ph(i,k) - ph(i,k+1))/(rho(i,k)*RG)
     3032              zzz(i,k+1) = zzz(i,k) + dzz(i,k)
     3033           END IF
     3034       END DO
     3035    END DO
     3036   
     3037    DO i =1, klon
     3038       IF (wk_active(i)) THEN
     3039           h_zzz(i) = max((1- omega(i))*zzz(i,ktop1(i)) + omega(i)*zzz(i,ktop1(i)+1), hwmin)
     3040       END IF
     3041    END DO
     3042
     3043
     3044ENDIF ! (iflag_wk_new_ptop/=0)
     3045
     3046!if (iflag_wk_new_ptop==2) then
     3047IF (iflag_wk_new_ptop>0) THEN
     3048   do i=1,klon
     3049      ptop(i)=ptop1(i)
     3050      ktop(i)=ktop1(i)
     3051      hw_(i)=h_zzz(i)
     3052   enddo
     3053
     3054!endif
     3055ENDIF
     3056
     3057 kupper = 0
     3058 
     3059IF (wk_pupper<1.) THEN
     3060 ! Choose an integration bound well above wake top
     3061  ! -----------------------------------------------------------------
     3062
     3063  ! Pupper = 50000.  ! melting level
     3064  ! Pupper = 60000.
     3065  ! Pupper = 80000.  ! essais pour case_e
     3066  DO i = 1, klon
     3067  !  pupper(i) = 0.6*ph(i, 1)
     3068    pupper(i) = wk_pupper*ph(i, 1)
     3069    pupper(i) = max(pupper(i), 45000.)
     3070    ! cc        Pupper(i) = 60000.
     3071  END DO
     3072
     3073ELSE
     3074  DO i=1, klon
     3075     ! pupper(i) = wk_pupper*ptop(i)+(1.-wk_pupper)*ph(i, 1)
     3076     !  pupper(i) = min( wk_pupper*ptop(i)+(1.-wk_pupper)*ph(i, 1) , ptop(i)-50.)
     3077      pupper(i) = min( wk_pupper*ptop(i)+(1.-wk_pupper)*ph(i, 1) , ptop(i)-5000.)
     3078  END DO
     3079END IF
     3080 
     3081  ! -5/ Determination de kupper
     3082
     3083  DO k = klev, 1, -1
     3084    DO i = 1, klon
     3085      IF (ph(i,k+1)<pupper(i)) kupper(i) = k
     3086    END DO
     3087  END DO
     3088
     3089  ! On evite kupper = 1 et kupper = klev
     3090  DO i = 1, klon
     3091    kupper(i) = max(kupper(i), 2)
     3092    kupper(i) = min(kupper(i), klev-1)
     3093  END DO
     3094  !---------- FIN nouveau calcul hw et ptop -------------------------------------
     3095
     3096IF (iflag_wk_new_ptop==999) THEN
     3097    DO i = 1, klon
     3098    hw_(i)=0.
     3099    ptop(i)=ph(i,1)
     3100    Ktop(i)=1
     3101    pupper(i)=ph(i,2)
     3102    kupper(i)=2
     3103    h_zzz(i)=0.
     3104    Ptop1(i)=ph(i,1)
     3105    ENDDO
     3106ENDIF
     3107
     3108zk_ptop_provis=k_ptop_provis
     3109
     3110    RETURN
     3111END SUBROUTINE pkupper
     3112
     3113
     3114SUBROUTINE wake_popdyn_1(klon, klev, dtime, cstar, tau_wk_inv, wgen, wdens, awdens, sigmaw, &
     3115                  wdensmin, &
     3116                  dtimesub, gfl, rad_wk, f_shear, drdt_pos, &
     3117                  d_awdens, d_wdens, d_sigmaw, &
     3118                  iflag_wk_act, wk_adv, cin, wape, &
     3119                  drdt, &
     3120                  d_dens_gen, d_dens_death, d_dens_col, d_dens_bnd, &
     3121                  d_sig_gen, d_sig_death, d_sig_col, d_sig_spread, d_sig_bnd, &
     3122                  d_wdens_targ, d_sigmaw_targ)
     3123               
     3124
     3125  USE lmdz_wake_ini , ONLY : wake_ini
     3126  USE lmdz_wake_ini , ONLY : prt_level,RG
     3127  USE lmdz_wake_ini , ONLY : stark, wdens_ref
     3128  USE lmdz_wake_ini , ONLY : tau_cv, rzero, aa0
     3129!!  USE lmdz_wake_ini , ONLY : iflag_wk_pop_dyn, wdensmin
     3130  USE lmdz_wake_ini , ONLY : iflag_wk_pop_dyn
     3131  USE lmdz_wake_ini , ONLY : sigmad, cstart, sigmaw_max
     3132 
     3133IMPLICIT NONE
     3134
     3135  INTEGER, INTENT(IN)                                   :: klon,klev
     3136  LOGICAL, DIMENSION (klon),        INTENT(IN)          :: wk_adv
     3137  REAL,                             INTENT(IN)          :: dtime
     3138  REAL,                             INTENT(IN)          :: dtimesub
     3139  REAL,                             INTENT(IN)          :: wdensmin
     3140  REAL, DIMENSION (klon),           INTENT(IN)          :: wgen
     3141  REAL, DIMENSION (klon),           INTENT(IN)          :: wdens
     3142  REAL, DIMENSION (klon),           INTENT(IN)          :: awdens
     3143  REAL, DIMENSION (klon),           INTENT(IN)          :: sigmaw
     3144  REAL, DIMENSION (klon),           INTENT(IN)          :: cstar
     3145  REAL, DIMENSION (klon),           INTENT(IN)          :: cin, wape
     3146  REAL, DIMENSION (klon),           INTENT(IN)          :: f_shear
     3147  INTEGER,                          INTENT(IN)          :: iflag_wk_act
     3148
     3149 
     3150  !
     3151 
     3152  ! Tendencies of state variables (2 is appended to the names of fields which are the cumul of fields
     3153  !                                 computed at each sub-timestep; e.g. d_wdens2 is the cumul of d_wdens)
     3154  REAL, DIMENSION (klon),           INTENT(OUT)         :: rad_wk
     3155  REAL, DIMENSION (klon),           INTENT(OUT)         :: gfl
     3156  REAL, DIMENSION (klon),           INTENT(OUT)         :: d_sigmaw, d_awdens, d_wdens
     3157  REAL, DIMENSION (klon),           INTENT(OUT)         :: drdt
     3158  ! Some components of the tendencies of state variables 
     3159  REAL, DIMENSION (klon),           INTENT(OUT)         :: d_sig_gen, d_sig_death, d_sig_col, d_sig_bnd
     3160  REAL, DIMENSION (klon),           INTENT(OUT)         :: d_sig_spread
     3161  REAL, DIMENSION (klon),           INTENT(OUT)         :: d_dens_gen, d_dens_death, d_dens_col, d_dens_bnd
     3162  REAL,                             INTENT(OUT)         :: d_wdens_targ, d_sigmaw_targ
     3163 
     3164 
     3165  REAL                                                  :: delta_t_min
     3166  INTEGER                                               :: i, k
     3167  REAL                                                  :: wdens0
     3168  ! IM 080208
     3169  LOGICAL, DIMENSION (klon)                             :: gwake
     3170 
     3171   ! Variables liees a la dynamique de population
     3172  REAL, DIMENSION(klon)                                 :: act
     3173  REAL, DIMENSION(klon)                                 :: tau_wk_inv
     3174  REAL, DIMENSION(klon)                                 :: wape1_act, wape2_act
     3175  LOGICAL, DIMENSION (klon)                             :: kill_wake
     3176  REAL                                                  :: drdt_pos
     3177  REAL                                                  :: tau_wk_inv_min
     3178 
     3179     
     3180
     3181      IF (iflag_wk_act == 0) THEN
     3182        act(:) = 0.
     3183      ELSEIF (iflag_wk_act == 1) THEN
     3184        act(:) = 1.
     3185      ELSEIF (iflag_wk_act ==2) THEN
     3186      DO i = 1, klon
     3187        IF (wk_adv(i)) THEN
     3188          wape1_act(i) = abs(cin(i))
     3189          wape2_act(i) = 2.*wape1_act(i) + 1.
     3190          act(i) = min(1., max(0., (wape(i)-wape1_act(i)) / (wape2_act(i)-wape1_act(i)) ))
     3191        ENDIF  ! (wk_adv(i))
     3192      ENDDO
     3193      ENDIF  ! (iflag_wk_act ==2)
     3194
     3195      DO i = 1, klon
     3196        IF (wk_adv(i)) THEN
     3197          rad_wk(i) = max( sqrt(sigmaw(i)/(3.14*wdens(i))) , rzero)
     3198          gfl(i)  = 2.*sqrt(3.14*wdens(i)*sigmaw(i))
     3199        END IF
     3200      END DO
     3201
     3202      DO i = 1, klon
     3203        IF (wk_adv(i)) THEN
     3204!!          tau_wk(i) = max(rad_wk(i)/(3.*cstar(i))*((cstar(i)/cstart)**1.5 - 1), 100.)
     3205          tau_wk_inv(i) = max( (3.*cstar(i))/(rad_wk(i)*((cstar(i)/cstart)**1.5 - 1)), 0.)
     3206          tau_wk_inv_min = min(tau_wk_inv(i), 1./dtimesub)
     3207          drdt(i) = (cstar(i) - wgen(i)*(sigmaw(i)/wdens(i)-aa0)/gfl(i)) / &
     3208                    (1 + 2*f_shear(i)*(2.*sigmaw(i)-aa0*wdens(i)) - 2.*sigmaw(i))
     3209!!                    (1 - 2*sigmaw(i)*(1.-f_shear(i)))
     3210          drdt_pos=max(drdt(i),0.)
     3211
     3212!!          d_wdens(i) = ( wgen(i)*(1.+2.*(sigmaw(i)-sigmad)) &
     3213!!                     - wdens(i)*tau_wk_inv_min &
     3214!!                     - 2.*gfl(i)*wdens(i)*Cstar(i) )*dtimesub
     3215!jyg+mlt<
     3216          d_awdens(i) = ( wgen(i) - (1./tau_cv)*(awdens(i) - act(i)*wdens(i)) )*dtimesub
     3217          d_dens_gen(i) = wgen(i)
     3218          d_dens_death(i) = - (wdens(i)-awdens(i))*tau_wk_inv_min
     3219          d_dens_col(i) =  -2.*wdens(i)*gfl(i)*drdt_pos
     3220          d_dens_gen(i) =  d_dens_gen(i)*dtimesub
     3221          d_dens_death(i) = d_dens_death(i)*dtimesub
     3222          d_dens_col(i) =  d_dens_col(i)*dtimesub
     3223
     3224          d_wdens(i) = d_dens_gen(i)+d_dens_death(i)+d_dens_col(i)
     3225!!          d_wdens(i) = ( wgen(i) - (wdens(i)-awdens(i))*tau_wk_inv_min -  &
     3226!!                         2.*wdens(i)*gfl(i)*drdt_pos )*dtimesub
     3227!>jyg+mlt
     3228!
     3229!jyg<
     3230          d_wdens_targ = max(d_wdens(i), wdensmin-wdens(i))
     3231!!          d_dens_bnd(i) = d_dens_bnd(i) + d_wdens_targ - d_wdens(i)
     3232          d_dens_bnd(i) = d_wdens_targ - d_wdens(i)
     3233          d_wdens(i) = d_wdens_targ
     3234!!          d_wdens(i) = max(d_wdens(i), wdensmin-wdens(i))
     3235!>jyg
     3236
     3237!jyg+mlt<
     3238!!          d_sigmaw(i) = ( (1.-2*f_shear(i)*sigmaw(i))*(gfl(i)*Cstar(i)+wgen(i)*sigmad/wdens(i)) &
     3239!!                      + 2.*f_shear(i)*wgen(i)*sigmaw(i)**2/wdens(i) &
     3240!!                      - sigmaw(i)*tau_wk_inv_min )*dtimesub
     3241          d_sig_gen(i) = wgen(i)*aa0
     3242          d_sig_death(i) = - sigmaw(i)*(1.-awdens(i)/wdens(i))*tau_wk_inv_min
     3243!!       
     3244         
     3245          d_sig_col(i) = - 2*f_shear(i)*sigmaw(i)*gfl(i)*drdt_pos
     3246          d_sig_col(i) = - 2*f_shear(i)*(2.*sigmaw(i)-wdens(i)*aa0)*gfl(i)*drdt_pos
     3247          d_sig_spread(i) = gfl(i)*cstar(i)
     3248          d_sig_gen(i) =  d_sig_gen(i)*dtimesub
     3249          d_sig_death(i) = d_sig_death(i)*dtimesub
     3250          d_sig_col(i) =  d_sig_col(i)*dtimesub
     3251          d_sig_spread(i) =  d_sig_spread(i)*dtimesub
     3252          d_sigmaw(i) =  d_sig_gen(i) + d_sig_death(i) + d_sig_col(i) + d_sig_spread(i)
     3253!>jyg+mlt
     3254!
     3255!jyg<
     3256          d_sigmaw_targ = max(d_sigmaw(i), sigmad-sigmaw(i))
     3257!!          d_sig_bnd(i) = d_sig_bnd(i) + d_sigmaw_targ - d_sigmaw(i)
     3258!!          d_sig_bnd_provis(i) = d_sigmaw_targ - d_sigmaw(i)
     3259          d_sig_bnd(i) = d_sigmaw_targ - d_sigmaw(i)
     3260          d_sigmaw(i) = d_sigmaw_targ
     3261!!          d_sigmaw(i) = max(d_sigmaw(i), sigmad-sigmaw(i))
     3262!>jyg
     3263        ENDIF
     3264      ENDDO
     3265
     3266      IF (prt_level >= 10) THEN
     3267        print *,'wake, cstar(1), cstar(1)/cstart, rad_wk(1), tau_wk_inv(1), drdt(1) ', &
     3268                       cstar(1), cstar(1)/cstart, rad_wk(1), tau_wk_inv(1), drdt(1)
     3269        print *,'wake, wdens(1), awdens(1), act(1), d_awdens(1) ', &
     3270                       wdens(1), awdens(1), act(1), d_awdens(1)
     3271        print *,'wake, wgen, -(wdens-awdens)*tau_wk_inv, -2.*wdens*gfl*drdt_pos, d_wdens ', &
     3272                       wgen(1), -(wdens(1)-awdens(1))*tau_wk_inv(1), -2.*wdens(1)*gfl(1)*drdt_pos, d_wdens(1)
     3273        print *,'wake, d_sig_gen(1), d_sig_death(1), d_sig_col(1), d_sigmaw(1) ', &
     3274                       d_sig_gen(1), d_sig_death(1), d_sig_col(1), d_sigmaw(1)
     3275      ENDIF
     3276   
     3277    RETURN
     3278    END SUBROUTINE wake_popdyn_1
     3279   
     3280    SUBROUTINE wake_popdyn_2 ( klon, klev, wk_adv, dtimesub, wgen, &
     3281                             wdensmin, &
     3282                             sigmaw, wdens, awdens, &   !! states variables
     3283                             gfl, cstar, cin, wape, rad_wk, &
     3284                             d_sigmaw, d_wdens, d_awdens, &  !! tendences
     3285                             cont_fact, &
     3286                             d_sig_gen, d_sig_death, d_sig_col, d_sig_spread, d_sig_bnd, &
     3287                             d_dens_gen, d_dens_death, d_dens_col, d_dens_bnd, &
     3288                             d_adens_death, d_adens_icol, d_adens_acol, d_adens_bnd )
     3289                             
     3290                                             
     3291
     3292  USE lmdz_wake_ini , ONLY : wake_ini
     3293  USE lmdz_wake_ini , ONLY : prt_level,RG
     3294  USE lmdz_wake_ini , ONLY : stark, wdens_ref
     3295  USE lmdz_wake_ini , ONLY : tau_cv, rzero, aa0
     3296!!  USE lmdz_wake_ini , ONLY : iflag_wk_pop_dyn, wdensmin
     3297  USE lmdz_wake_ini , ONLY : iflag_wk_pop_dyn
     3298  USE lmdz_wake_ini , ONLY : sigmad, cstart, sigmaw_max
     3299 
     3300IMPLICIT NONE
     3301
     3302  INTEGER, INTENT(IN)                                   :: klon,klev
     3303  LOGICAL, DIMENSION (klon),        INTENT(IN)          :: wk_adv
     3304  REAL,                             INTENT(IN)          :: dtimesub
     3305  REAL,                             INTENT(IN)          :: wdensmin
     3306  REAL, DIMENSION (klon),           INTENT(IN)          :: wgen      !! B = birth rate of wakes
     3307  REAL, DIMENSION (klon),           INTENT(INOUT)       :: sigmaw    !! sigma = fractional area of wakes
     3308  REAL, DIMENSION (klon),           INTENT(INOUT)       :: wdens     !! D = number of wakes per unit area
     3309  REAL, DIMENSION (klon),           INTENT(INOUT)       :: awdens    !! A = number of active wakes per unit area
     3310  REAL, DIMENSION (klon),           INTENT(IN)          :: cstar     !! C* = spreading velocity of wakes
     3311  REAL, DIMENSION (klon),           INTENT(IN)          :: cin, wape  ! RM : A Faire disparaitre
     3312
     3313! 
     3314  REAL, DIMENSION (klon),           INTENT(OUT)         :: rad_wk    !! r = wake radius
     3315  REAL, DIMENSION (klon),           INTENT(OUT)         :: gfl       !! Lg = gust front lenght per unit area
     3316  REAL, DIMENSION (klon),           INTENT(OUT)         :: d_sigmaw, d_wdens, d_awdens
     3317  REAL, DIMENSION (klon),           INTENT(OUT)         :: cont_fact  !! RM facteur de contact = 2 pi * rad * C*
     3318  ! Some components of the tendencies of state variables 
     3319  REAL, DIMENSION (klon),           INTENT(OUT)         :: d_sig_gen, d_sig_death, d_sig_col, d_sig_spread, d_sig_bnd
     3320  REAL, DIMENSION (klon),           INTENT(OUT)         :: d_dens_gen, d_dens_death, d_dens_col, d_dens_bnd
     3321  REAL, DIMENSION (klon),           INTENT(OUT)         :: d_adens_death, d_adens_icol, d_adens_acol, d_adens_bnd
     3322
     3323
     3324!! internal variables
     3325 
     3326  INTEGER                                               :: i, k
     3327  REAL, DIMENSION (klon)                                :: tau_wk_inv      !! tau = life time of wakes
     3328  REAL                                                  :: tau_wk_inv_min
     3329  REAL, DIMENSION (klon)                                :: tau_prime       !! tau_prime = life time of actives wakes
     3330  REAL                                                  :: d_wdens_targ, d_sigmaw_targ
     3331 
     3332
     3333!! Equations
     3334!! dD/dt = B - (D-A)/tau - f D^2
     3335!! dA/dt = B - A/tau_prime + f (D-A)^2 - f A^2
     3336!! dsigma/dt = B a0 - sigma/D (D-A)/tau + Lg C* - f (D-A)^2 (sigma/D-a0)
     3337!!
     3338!! f = 2 (B (a0-sigma/D) + Lg C*) / (2 (D-A)^2 (2 sigma/D-a0) + D (1-2 sigma))
     3339
     3340
     3341      DO i = 1, klon
     3342        IF (wk_adv(i)) THEN
     3343          rad_wk(i) = max( sqrt(sigmaw(i)/(3.14*wdens(i))) , rzero)
     3344          gfl(i)  = 2.*sqrt(3.14*wdens(i)*sigmaw(i))
     3345        END IF
     3346      END DO
     3347
     3348
     3349      DO i = 1, klon
     3350        IF (wk_adv(i)) THEN
     3351!!          tau_wk(i) = max(rad_wk(i)/(3.*cstar(i))*((cstar(i)/cstart)**1.5 - 1), 100.)
     3352          tau_wk_inv(i) = max( (3.*cstar(i))/(rad_wk(i)*((cstar(i)/cstart)**1.5 - 1)), 0.)
     3353          tau_wk_inv_min = min(tau_wk_inv(i), 1./dtimesub)
     3354          tau_prime(i) = tau_cv
     3355!!          cont_fact(i) = 2.*(wgen(i)*(aa0-sigmaw(i)/wdens(i)) + gfl(i)*cstar(i)) / &
     3356!!                             (2.*(wdens(i)-awdens(i))**2*(2.*sigmaw(i)/wdens(i) - aa0) + wdens(i)*(1.-2.*sigmaw(i)))
     3357!!          cont_fact(i) = 2.*3.14*rad_wk(i)*cstar(i)     ! bug
     3358!!          cont_fact(i) = 4.*3.14*rad_wk(i)*cstar(i)
     3359          cont_fact(i) = 2.*gfl(i)*cstar(i)/wdens(i)
     3360
     3361          d_sig_gen(i) = wgen(i)*aa0
     3362          d_sig_death(i) = - sigmaw(i)*(1.-awdens(i)/wdens(i))*tau_wk_inv_min
     3363          d_sig_col(i) = - cont_fact(i)*(wdens(i)-awdens(i))**2*(2.*sigmaw(i)/wdens(i)-aa0)
     3364          d_sig_spread(i) = gfl(i)*cstar(i)
     3365!
     3366          d_sig_gen(i) =  d_sig_gen(i)*dtimesub
     3367          d_sig_death(i) = d_sig_death(i)*dtimesub
     3368          d_sig_col(i) =  d_sig_col(i)*dtimesub
     3369          d_sig_spread(i) =  d_sig_spread(i)*dtimesub
     3370          d_sigmaw(i) =  d_sig_gen(i) + d_sig_death(i) + d_sig_col(i) + d_sig_spread(i)
     3371
     3372         
     3373          d_sigmaw_targ = max(d_sigmaw(i), sigmad-sigmaw(i))
     3374!!          d_sig_bnd(i) = d_sig_bnd(i) + d_sigmaw_targ - d_sigmaw(i)
     3375!!          d_sig_bnd_provis(i) = d_sigmaw_targ - d_sigmaw(i)
     3376          d_sig_bnd(i) = d_sigmaw_targ - d_sigmaw(i)
     3377          d_sigmaw(i) = d_sigmaw_targ
     3378!!          d_sigmaw(i) = max(d_sigmaw(i), sigmad-sigmaw(i))
     3379         
     3380         
     3381          d_dens_gen(i) = wgen(i)
     3382          d_dens_death(i) = - (wdens(i)-awdens(i))*tau_wk_inv_min
     3383          d_dens_col(i) =  - cont_fact(i)*wdens(i)**2
     3384!
     3385          d_dens_gen(i) =  d_dens_gen(i)*dtimesub
     3386          d_dens_death(i) = d_dens_death(i)*dtimesub
     3387          d_dens_col(i) =  d_dens_col(i)*dtimesub
     3388          d_wdens(i) = d_dens_gen(i) + d_dens_death(i) + d_dens_col(i)
     3389
     3390
     3391          d_adens_death(i) = -awdens(i)/tau_prime(i)
     3392          d_adens_icol(i) = cont_fact(i)*(wdens(i)-awdens(i))**2
     3393          d_adens_acol(i) = - cont_fact(i)*awdens(i)**2
     3394!
     3395          d_adens_death(i) =  d_adens_death(i)*dtimesub
     3396          d_adens_icol(i) =   d_adens_icol(i)*dtimesub
     3397          d_adens_acol(i) =   d_adens_acol(i)*dtimesub
     3398          d_awdens(i) =   d_dens_gen(i) + d_adens_death(i) + d_adens_icol(i) + d_adens_acol(i)     
     3399           
     3400!!
     3401          d_wdens_targ = max(d_wdens(i), wdensmin-wdens(i))
     3402!!          d_dens_bnd(i) = d_dens_bnd(i) + d_wdens_targ - d_wdens(i)
     3403          d_dens_bnd(i) = d_wdens_targ - d_wdens(i)
     3404          d_wdens(i) = d_wdens_targ
     3405         
     3406          d_wdens_targ = min(max(d_awdens(i),-awdens(i)), wdens(i)-awdens(i))
     3407!!          d_dens_bnd(i) = d_dens_bnd(i) + d_wdens_targ - d_wdens(i)
     3408          d_adens_bnd(i) = d_wdens_targ - d_awdens(i)
     3409          d_awdens(i) = d_wdens_targ
     3410
     3411
     3412
     3413        ENDIF
     3414      ENDDO
     3415
     3416      IF (prt_level >= 10) THEN
     3417        print *,'wake, cstar(1), cstar(1)/cstart, rad_wk(1), tau_wk_inv(1), cont_fact(1) ', &
     3418                       cstar(1), cstar(1)/cstart, rad_wk(1), tau_wk_inv(1), cont_fact(1)
     3419        print *,'wake, wdens(1), awdens(1), d_awdens(1) ', &
     3420                       wdens(1), awdens(1), d_awdens(1)
     3421        print *,'wake, d_sig_gen(1), d_sig_death(1), d_sig_col(1), d_sigmaw(1) ', &
     3422                       d_sig_gen(1), d_sig_death(1), d_sig_col(1), d_sigmaw(1)
     3423      ENDIF
     3424sigmaw=sigmaw+d_sigmaw
     3425wdens=wdens+d_wdens
     3426awdens=awdens+d_awdens
     3427
     3428    RETURN
     3429    END SUBROUTINE wake_popdyn_2 
     3430 
     3431    SUBROUTINE wake_popdyn_3 ( klon, klev, phys_sub, wk_adv, dtimesub, wgen, &
     3432                             wdensmin, &
     3433                             sigmaw, asigmaw, wdens, awdens, &                       !! state variables
     3434                             gfl, agfl, cstar, cin, wape, &
     3435                             rad_wk, arad_wk, irad_wk, &
     3436                             d_sigmaw, d_asigmaw, d_wdens, d_awdens, &               !! tendencies
     3437                             d_sig_gen, d_sig_death, d_sig_col, d_sig_spread, d_sig_bnd, &
     3438                             d_asig_death, d_asig_aicol, d_asig_iicol, d_asig_spread, d_asig_bnd, &
     3439                             d_dens_gen, d_dens_death, d_dens_col, d_dens_bnd, &
     3440                             d_adens_death, d_adens_icol, d_adens_acol, d_adens_bnd )
     3441                             
     3442                                             
     3443
     3444  USE lmdz_wake_ini , ONLY : wake_ini
     3445  USE lmdz_wake_ini , ONLY : prt_level,RG
     3446  USE lmdz_wake_ini , ONLY : stark, wdens_ref
     3447  USE lmdz_wake_ini , ONLY : tau_cv, rzero, aa0
     3448!!  USE lmdz_wake_ini , ONLY : iflag_wk_pop_dyn, wdensmin
     3449  USE lmdz_wake_ini , ONLY : iflag_wk_pop_dyn
     3450  USE lmdz_wake_ini , ONLY : sigmad, cstart, sigmaw_max
     3451  USE lmdz_wake_ini , ONLY : smallestreal
     3452 
     3453IMPLICIT NONE
     3454
     3455  INTEGER, INTENT(IN)                                   :: klon,klev
     3456  LOGICAL,                          INTENT(IN)          :: phys_sub
     3457  LOGICAL, DIMENSION (klon),        INTENT(IN)          :: wk_adv
     3458  REAL,                             INTENT(IN)          :: dtimesub
     3459  REAL,                             INTENT(IN)          :: wdensmin
     3460  REAL, DIMENSION (klon),           INTENT(IN)          :: wgen      !! B = birth rate of wakes
     3461  REAL, DIMENSION (klon),           INTENT(INOUT)       :: sigmaw    !! sigma = fractional area of wakes
     3462  REAL, DIMENSION (klon),           INTENT(INOUT)       :: asigmaw   !! sigma = fractional area of active wakes
     3463  REAL, DIMENSION (klon),           INTENT(INOUT)       :: wdens     !! D = number of wakes per unit area
     3464  REAL, DIMENSION (klon),           INTENT(INOUT)       :: awdens    !! A = number of active wakes per unit area
     3465  REAL, DIMENSION (klon),           INTENT(IN)          :: cstar     !! C* = spreading velocity of wakes
     3466  REAL, DIMENSION (klon),           INTENT(IN)          :: cin, wape  ! RM : A Faire disparaitre
     3467
     3468! 
     3469  REAL, DIMENSION (klon),           INTENT(OUT)         :: rad_wk    !! r = mean wake radius
     3470  REAL, DIMENSION (klon),           INTENT(OUT)         :: arad_wk    !! r_A = wake radius of active wakes
     3471  REAL, DIMENSION (klon),           INTENT(OUT)         :: irad_wk    !! r_I = wake radius of inactive wakes
     3472  REAL, DIMENSION (klon),           INTENT(OUT)         :: gfl       !! Lg = gust front length per unit area
     3473  REAL, DIMENSION (klon),           INTENT(OUT)         :: agfl      !! LgA = gust front length of active wakes
     3474                                                                     !!  per unit area
     3475  REAL, DIMENSION (klon),           INTENT(OUT)         :: d_sigmaw, d_asigmaw, d_wdens, d_awdens
     3476  ! Some components of the tendencies of state variables 
     3477  REAL, DIMENSION (klon),           INTENT(OUT)         :: d_sig_gen, d_sig_death, d_sig_col, d_sig_spread, d_sig_bnd
     3478  REAL, DIMENSION (klon),           INTENT(OUT)         :: d_asig_death, d_asig_aicol, d_asig_iicol, d_asig_spread, d_asig_bnd
     3479  REAL, DIMENSION (klon),           INTENT(OUT)         :: d_dens_gen, d_dens_death, d_dens_col, d_dens_bnd
     3480  REAL, DIMENSION (klon),           INTENT(OUT)         :: d_adens_death, d_adens_acol, d_adens_icol, d_adens_bnd
     3481
     3482
     3483!! internal variables
     3484 
     3485  INTEGER                                               :: i, k
     3486  REAL, DIMENSION (klon)                                :: iwdens, isigmaw !! inactive wake density and fractional area
     3487!!  REAL, DIMENSION (klon)                                :: d_arad, d_irad
     3488  REAL, DIMENSION (klon)                                :: igfl            !! LgI = gust front length of inactive wakes
     3489                                                                           !!  per unit area
     3490  REAL, DIMENSION (klon)                                :: s_wk            !! mean area of individual wakes
     3491  REAL, DIMENSION (klon)                                :: as_wk           !! mean area of individual active wakes
     3492  REAL, DIMENSION (klon)                                :: is_wk           !! mean area of individual inactive wakes
     3493  REAL, DIMENSION (klon)                                :: tau_wk_inv      !! tau = life time of wakes
     3494  REAL                                                  :: tau_wk_inv_min
     3495  REAL, DIMENSION (klon)                                :: tau_prime       !! tau_prime = life time of actives wakes
     3496  REAL                                                  :: d_wdens_targ, d_sigmaw_targ
     3497
     3498
     3499!! Equations
     3500!! ---------
     3501!! Gust fronts:
     3502!! Lg_A = 2 pi r_A A
     3503!! Lg_I = 2 pi r_I I
     3504!! Lg   = 2 pi r   D
     3505!!
     3506!! Areas:
     3507!! s = pi r^2
     3508!! s_A = pi r_A^2
     3509!! s_I = pi r_I^2
     3510!!
     3511!! Life expectancy:
     3512!! tau_I = 3 C* ((C*/C*t)^3/2 - 1) / r_I
     3513!!
     3514!! Time deratives:
     3515!! dD/dt = B - (D-A)/tau_I - 2 Lg C* D
     3516!! dA/dt = B - A/tau_A     + 2 Lg_I C* (D-A) - 2 Lg_A C* A
     3517!! dsigma/dt = B a0 - sigma_I/tau_I + Lg C* - 2 Lg_I C* (D-A) (2 s_I - a0)
     3518!! dsigma_A/dt = B a0 - sigma_A/tau_A + Lg_A C* + (Lg_A I + Lg_I A) C* s_I + 2 Lg_I C* I a0
     3519!!
     3520
     3521      DO i = 1, klon
     3522        IF (wk_adv(i)) THEN
     3523         iwdens(i) = wdens(i) - awdens(i)
     3524         isigmaw(i) = sigmaw(i) - asigmaw(i)
     3525!
     3526         arad_wk(i) = max( sqrt(asigmaw(i)/(3.14*awdens(i))) , rzero)
     3527         irad_wk(i) = max( sqrt((sigmaw(i)-asigmaw(i))/  &
     3528                           (3.14*max(smallestreal,(wdens(i)-awdens(i))))), rzero)
     3529         rad_wk(i) = (awdens(i)*arad_wk(i)+(wdens(i)-awdens(i))*irad_wk(i))/wdens(i)
     3530!
     3531         s_wk(i) = 3.14*rad_wk(i)**2
     3532         as_wk(i) = 3.14*arad_wk(i)**2
     3533         is_wk(i) = 3.14*irad_wk(i)**2
     3534!
     3535         gfl(i)  = 2.*sqrt(3.14*wdens(i)*sigmaw(i))
     3536         agfl(i) = 2.*sqrt(3.14*awdens(i)*asigmaw(i))
     3537         igfl(i) = gfl(i) - agfl(i)
     3538        ENDIF
     3539      ENDDO
     3540
     3541
     3542      DO i = 1, klon
     3543        IF (wk_adv(i)) THEN
     3544          tau_wk_inv(i) = max( (3.*cstar(i))/(irad_wk(i)*((cstar(i)/cstart)**1.5 - 1)), 0.)
     3545          tau_wk_inv_min = min(tau_wk_inv(i), 1./dtimesub)
     3546          tau_prime(i) = tau_cv
     3547
     3548          d_sig_gen(i) = wgen(i)*aa0
     3549          d_sig_death(i) = - isigmaw(i)*tau_wk_inv_min
     3550          d_sig_col(i) = - 2.*igfl(i)*cstar(i)*iwdens(i)*(2.*is_wk(i)-aa0)
     3551          d_sig_spread(i) = gfl(i)*cstar(i)
     3552!
     3553          d_sig_gen(i) =  d_sig_gen(i)*dtimesub
     3554          d_sig_death(i) = d_sig_death(i)*dtimesub
     3555          d_sig_col(i) =  d_sig_col(i)*dtimesub
     3556          d_sig_spread(i) =  d_sig_spread(i)*dtimesub
     3557          d_sigmaw(i) =  d_sig_gen(i) + d_sig_death(i) + d_sig_col(i) + d_sig_spread(i)
     3558#ifdef IOPHYS_WK
     3559          IF (phys_sub) call iophys_ecrit('d_sigmaw0',1,'d_sigmaw0','',d_sigmaw)
     3560#endif
     3561
     3562         
     3563          d_sigmaw_targ = max(d_sigmaw(i), sigmad-sigmaw(i))
     3564!!          d_sig_bnd(i) = d_sig_bnd(i) + d_sigmaw_targ - d_sigmaw(i)
     3565!!          d_sig_bnd_provis(i) = d_sigmaw_targ - d_sigmaw(i)
     3566          d_sig_bnd(i) = d_sigmaw_targ - d_sigmaw(i)
     3567          d_sigmaw(i) = d_sigmaw_targ
     3568!!          d_sigmaw(i) = max(d_sigmaw(i), sigmad-sigmaw(i))
     3569#ifdef IOPHYS_WK
     3570          IF (phys_sub) THEN
     3571             call iophys_ecrit('tauwk_inv',1,'tau_wk_inv_min','',tau_wk_inv_min)
     3572             call iophys_ecrit('d_sigmaw',1,'d_sigmaw','',d_sigmaw)
     3573             call iophys_ecrit('d_sig_gen',1,'d_sig_gen','',d_sig_gen)
     3574             call iophys_ecrit('d_sig_death',1,'d_sig_death','',d_sig_death)
     3575             call iophys_ecrit('d_sig_col',1,'d_sig_col','',d_sig_col)
     3576             call iophys_ecrit('d_sig_spread',1,'d_sig_spread','',d_sig_spread)
     3577             call iophys_ecrit('d_sig_bnd',1,'d_sig_bnd','',d_sig_bnd)
     3578          ENDIF
     3579#endif
     3580          d_asig_death(i) = - asigmaw(i)/tau_prime(i)
     3581          d_asig_aicol(i) = (agfl(i)*iwdens(i) + igfl(i)*awdens(i))*cstar(i)*is_wk(i)
     3582          d_asig_iicol(i) = 2.*igfl(i)*cstar(i)*iwdens(i)*aa0
     3583          d_asig_spread(i) = agfl(i)*cstar(i)
     3584!
     3585          d_asig_death(i) = d_asig_death(i)*dtimesub
     3586          d_asig_aicol(i) =  d_asig_aicol(i)*dtimesub
     3587          d_asig_iicol(i) =  d_asig_iicol(i)*dtimesub
     3588          d_asig_spread(i) =  d_asig_spread(i)*dtimesub
     3589          d_asigmaw(i) =  d_sig_gen(i) + d_asig_death(i) + d_asig_aicol(i) + d_asig_iicol(i) + d_asig_spread(i)
     3590#ifdef IOPHYS_WK
     3591          IF (phys_sub) call iophys_ecrit('d_asigmaw0',1,'d_asigmaw0','',d_asigmaw)
     3592#endif
     3593
     3594          d_sigmaw_targ = min(max(d_asigmaw(i),-asigmaw(i)), sigmaw(i)-asigmaw(i))
     3595!!          d_dens_bnd(i) = d_dens_bnd(i) + d_sigmaw_targ - d_sigmaw(i)
     3596          d_asig_bnd(i) = d_sigmaw_targ - d_asigmaw(i)
     3597          d_asigmaw(i) = d_sigmaw_targ
     3598#ifdef IOPHYS_WK
     3599          IF (phys_sub) THEN
     3600             call iophys_ecrit('d_asigmaw',1,'d_asigmaw','',d_asigmaw)
     3601             call iophys_ecrit('d_asig_death',1,'d_asig_death','',d_asig_death)
     3602             call iophys_ecrit('d_asig_aicol',1,'d_asig_aicol','',d_asig_aicol)
     3603             call iophys_ecrit('d_asig_iicol',1,'d_asig_iicol','',d_asig_iicol)
     3604             call iophys_ecrit('d_asig_spread',1,'d_asig_spread','',d_asig_spread)
     3605             call iophys_ecrit('d_asig_bnd',1,'d_asig_bnd','',d_asig_bnd)
     3606          ENDIF
     3607#endif
     3608          d_dens_gen(i) = wgen(i)
     3609          d_dens_death(i) = - iwdens(i)*tau_wk_inv_min
     3610          d_dens_col(i) =  - 2.*gfl(i)*cstar(i)*wdens(i)
     3611!
     3612          d_dens_gen(i) =  d_dens_gen(i)*dtimesub
     3613          d_dens_death(i) = d_dens_death(i)*dtimesub
     3614          d_dens_col(i) =  d_dens_col(i)*dtimesub
     3615          d_wdens(i) = d_dens_gen(i) + d_dens_death(i) + d_dens_col(i)
     3616!!
     3617          d_wdens_targ = max(d_wdens(i), wdensmin-wdens(i))
     3618!!          d_dens_bnd(i) = d_dens_bnd(i) + d_wdens_targ - d_wdens(i)
     3619          d_dens_bnd(i) = d_wdens_targ - d_wdens(i)
     3620          d_wdens(i) = d_wdens_targ
     3621#ifdef IOPHYS_WK
     3622    IF (phys_sub) THEN
     3623        call iophys_ecrit('d_wdens',1,'d_wdens','',d_wdens)
     3624        call iophys_ecrit('d_dens_gen',1,'d_dens_gen','',d_dens_gen)
     3625        call iophys_ecrit('d_dens_death',1,'d_dens_death','',d_dens_death)
     3626        call iophys_ecrit('d_dens_col',1,'d_dens_col','',d_dens_col)
     3627    ENDIF
     3628#endif
     3629
     3630          d_adens_death(i) = -awdens(i)/tau_prime(i)
     3631          d_adens_icol(i) =   2.*igfl(i)*cstar(i)*iwdens(i)
     3632          d_adens_acol(i)  = - 2.*agfl(i)*cstar(i)*awdens(i)
     3633!
     3634          d_adens_death(i) =  d_adens_death(i)*dtimesub
     3635          d_adens_icol(i) =   d_adens_icol(i)*dtimesub
     3636          d_adens_acol(i)  =   d_adens_acol(i)*dtimesub
     3637          d_awdens(i) =   d_dens_gen(i) + d_adens_death(i) + d_adens_icol(i) + d_adens_acol(i)     
     3638#ifdef IOPHYS_WK
     3639    IF (phys_sub) THEN
     3640        call iophys_ecrit('d_awdens',1,'d_awdens','',d_awdens)
     3641        call iophys_ecrit('d_adens_death',1,'d_adens_death','',d_adens_death)
     3642        call iophys_ecrit('d_adens_icol',1,'d_adens_icol','',d_adens_icol)
     3643        call iophys_ecrit('d_adens_acol',1,'d_adens_acol','',d_adens_acol)
     3644    ENDIF
     3645#endif
     3646          d_wdens_targ = min(max(d_awdens(i),-awdens(i)), wdens(i)-awdens(i))
     3647!!          d_dens_bnd(i) = d_dens_bnd(i) + d_wdens_targ - d_wdens(i)
     3648          d_adens_bnd(i) = d_wdens_targ - d_awdens(i)
     3649          d_awdens(i) = d_wdens_targ
     3650
     3651!!          d_irad(i) = (d_sigmaw(i)-d_asigmaw(i)-isigmaw(i)*(d_wdens(i)-awdens(i))/iwdens(i)) / &
     3652!!                      max(smallestreal,(2.*3.14*iwdens(i)*irad_wk(i)))
     3653!!          d_arad(i) = (d_asigmaw(i)-asigmaw(i)*d_awdens(i)/awdens(i)) / &
     3654!!                      max(smallestreal,(2.*3.14*awdens(i)*arad_wk(i)))
     3655!!          d_irad(i) = d_irad(i)*dtimesub
     3656!!          d_arad(i) = d_arad(i)*dtimesub
     3657!!          call iophys_ecrit('d_irad',1,'d_irad','m',d_irad)
     3658!!          call iophys_ecrit('d_airad',1,'d_arad','m',d_arad)
     3659!!
     3660        ENDIF
     3661      ENDDO
     3662
     3663      IF (prt_level >= 10) THEN
     3664        print *,'wake, cstar(1), cstar(1)/cstart, rad_wk(1), tau_wk_inv(1), gfl(1) ', &
     3665                       cstar(1), cstar(1)/cstart, rad_wk(1), tau_wk_inv(1), gfl(1)
     3666        print *,'wake, wdens(1), awdens(1), d_awdens(1) ', &
     3667                       wdens(1), awdens(1), d_awdens(1)
     3668        print *,'wake, d_sig_gen(1), d_sig_death(1), d_sig_col(1), d_sigmaw(1) ', &
     3669                       d_sig_gen(1), d_sig_death(1), d_sig_col(1), d_sigmaw(1)
     3670      ENDIF
     3671sigmaw=sigmaw+d_sigmaw
     3672asigmaw=asigmaw+d_asigmaw
     3673wdens=wdens+d_wdens
     3674awdens=awdens+d_awdens
     3675
     3676    RETURN
     3677    END SUBROUTINE wake_popdyn_3 
     3678
    26293679END MODULE lmdz_wake

  • LMDZ6/trunk/libf/phylmdiso/physiq_mod.F90

    r5253 r5255  
    329329       d_deltat_ajs_cv, d_deltaq_ajs_cv, & ! due to dry adjustment of (w) before convection
    330330                       ! tendencies of wake fractional area and wake number per unit area:
    331        d_s_wk, d_dens_a_wk,  d_dens_wk, &  ! due to wakes
     331       d_s_wk, d_s_a_wk, d_dens_a_wk,  d_dens_wk, &  ! due to wakes
    332332!!!       d_s_vdf, d_dens_a_vdf, d_dens_vdf, & ! due to vertical diffusion
    333333!!!       d_s_the, d_dens_a_the, d_dens_the, & ! due to thermals
    … …  
    14141414    endif                                                                  !
    14151415    !======================================================================!
     1416#ifdef ISO
     1417#ifdef ISOVERIF
     1418        write(*,*) 'physiq 1421b: qx(1,1,:)=',qx(1,1,:)
     1419#endif
     1420#endif
    14161421
    14171422
    … …  
    44454450               dt_a, dq_a, cv_gen,  &
    44464451               sigd, cin,  &
    4447                wake_deltat, wake_deltaq, wake_s, awake_dens, wake_dens,  &
     4452               wake_deltat, wake_deltaq, wake_s, awake_s, wake_dens, awake_dens,  &
    44484453               wake_dth, wake_h,  &
    44494454!!               wake_pe, wake_fip, wake_gfl,  &
    … …  
    44554460               wake_omg, wake_dp_deltomg,  &
    44564461               wake_spread, wake_Cstar, d_deltat_wk_gw,  &
    4457                d_deltat_wk, d_deltaq_wk, d_s_wk, d_dens_a_wk, d_dens_wk &
     4462               d_deltat_wk, d_deltaq_wk, d_s_wk, d_s_a_wk, d_dens_wk, d_dens_a_wk &
    44584463#ifdef ISO
    44594464     &               ,xt_seri,dxt_dwn,dxt_a &
    … …  
    44954500
    44964501         CALL add_wake_tend &
    4497             (d_deltat_wk, d_deltaq_wk, d_s_wk, d_dens_a_wk, d_dens_wk, wake_k, &
     4502            (d_deltat_wk, d_deltaq_wk, d_s_wk, d_s_a_wk, d_dens_wk, d_dens_a_wk, wake_k, &
    44984503             'wake', abortphy &
    44994504#ifdef ISO
    … …  
    47704775               ,zqla,ztva &
    47714776#ifdef ISO         
    4772      &      ,xt_env,d_xt_ajs &
     4777     &      ,xt_therm,d_xt_ajs &
    47734778#ifdef DIAGISO         
    47744779     &      ,q_the,xt_the &
    … …  
    65576562                IF(flag_verbose_strataer) print *,'IN physiq_mod: min max d_q_emiss=',&
    65586563                     minval(d_q_emiss),maxval(d_q_emiss)
     6564#ifdef ISO
     6565       abort_message='isotopes pas encore dans emission volc H2O'
     6566       CALL abort_physic(modname,abort_message,1)
     6567#endif
    65596568
    65606569                CALL add_phys_tend(du0, dv0, dt0, d_q_emiss, dql0, dqi0, dqbs0, paprs, &
    … …  
    69486957
    69496958          print*,'avt add phys rep',abortphy
    6950 
     6959#ifdef ISO
     6960       abort_message='isotopes pas encore dans rep'
     6961       CALL abort_physic(modname,abort_message,1)
     6962#endif
    69516963     CALL add_phys_tend &
    69526964            (du0,dv0,dt0,d_q_rep,d_ql_rep,d_qi_rep,dqbs0,paprs,&
Note: See TracChangeset for help on using the changeset viewer.