Changeset 5803 for LMDZ6


Ignore:
Timestamp:
Sep 5, 2025, 12:52:13 AM (5 months ago)
Author:
fhourdin
Message:

Work in progress on wake (JYG&FH)

Location:
LMDZ6/trunk/libf/phylmd
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • LMDZ6/trunk/libf/phylmd/calwake.f90

    r5770 r5803  
    11! $Id$
    2 !$gpum horizontal klon nlon
    32MODULE calwake_mod
    43  PRIVATE
    … …  
    5857  USE indice_sol_mod, ONLY: is_oce
    5958  USE print_control_mod, ONLY: lunout, prt_level
    60   USE lmdz_wake, ONLY : wake, wake2
     59  USE lmdz_wake, ONLY : wake, wake2, wake3
    6160  USE alpale_mod, ONLY: iflag_wake
    6261  USE yomcst_mod_h
    … …  
    139138    print *, '-> calwake, wake_s, awake_s, wgen input ', wake_s(1), awake_s(1), wgen(1)
    140139  ENDIF
     140  print*,'NOUVEAU WAKE ================================='
    141141
    142142  rdcp = 1./3.5
    … …  
    256256  ELSE IF (iflag_wake/10 == 2) THEN
    257257    CALL wake2(klon,klev,znatsurf, p, ph, pi, dtime, &
     258      te, qe, omgbe, &
     259      dtdwn, dqdwn, amdwn, amup, dta, dqa, wgen, &
     260      sigd0, Cin, &
     261      dtw, dqw, sigmaw, asigmaw, wdens, awdens, &                      ! state variables
     262      dth, hw, wape, fip, gfl, &
     263      dtls, dqls, ktopw, omgbdth, dp_omgb, tx, qx, &
     264      dtke, dqke, omg, dp_deltomg, spread, cstar, &
     265      d_deltat_gw, &
     266      d_deltatw, d_deltaqw, d_sigmaw, d_asigmaw, d_wdens, d_awdens)      ! tendencies
     267
     268  ELSE IF (iflag_wake/10 == 3) THEN
     269    CALL wake3(klon,klev,znatsurf, p, ph, pi, dtime, &
    258270      te, qe, omgbe, &
    259271      dtdwn, dqdwn, amdwn, amup, dta, dqa, wgen, &

  • LMDZ6/trunk/libf/phylmd/lmdz_wake.f90

    r5775 r5803  
    99  !$OMP THREADPRIVATE(first_call)
    1010
    11   PUBLIC wake, wake2, wake_first
     11  PUBLIC wake, wake2, wake3, wake_first, wake_popdyn_3
    1212
    1313CONTAINS
    … …  
    20342034    DO i = 1, klon
    20352035      IF (ok_qx_qw(i)) THEN
    2036         !! print *,'wake, hw0(i), dthmin(i) ', hw0(i), dthmin(i)
    2037         !! print *,'wake, 2.*sum_dth(i)/(hw0(i)*dthmin(i)) ', &
    2038         !!                2.*sum_dth(i)/(hw0(i)*dthmin(i))
    2039         !! print *,'wake, sum_half_dth(i), sum_dth(i) ', &
    2040         !!                sum_half_dth(i), sum_dth(i)
     2036         !! print *,'wake, hw0(i), dthmin(i) ', hw0(i), dthmin(i)
     2037         !! print *,'wake, 2.*sum_dth(i)/(hw0(i)*dthmin(i)) ', &
     2038         !!               2.*sum_dth(i)/(hw0(i)*dthmin(i))
     2039         !! print *,'wake, sum_half_dth(i), sum_dth(i) ', &
     2040         !!               sum_half_dth(i), sum_dth(i)
    20412041        IF ((hw0(i) < 1.) .or. (dthmin(i) >= -delta_t_min) ) THEN
    20422042          wape2(i) = -1.
    … …  
    37123712    d_deltaq_dadv(:,:) = d_deltaq_dadv(:,:)/dtimesub
    37133713!
    3714     !   For the mean fields: tb and qb the computation of the tendencies due to wakes is
     3714    !   For the mean fields tb and qb the computation of the tendencies due to wakes is
    37153715    !   already complete.
    37163716    d_tb(:,:) = d_tb_dadv(:,:)
    … …  
    44564456    DO i = 1, klon
    44574457      IF (ok_qx_qw(i)) THEN
    4458         !! print *,'wake, hw0(i), dthmin(i) ', hw0(i), dthmin(i)
    4459         !! print *,'wake, 2.*sum_dth(i)/(hw0(i)*dthmin(i)) ', &
    4460         !!                2.*sum_dth(i)/(hw0(i)*dthmin(i))
    4461         !! print *,'wake, sum_half_dth(i), sum_dth(i) ', &
    4462         !!                sum_half_dth(i), sum_dth(i)
    4463         IF ((hw0(i) < 1.) .or. (dthmin(i) >= -delta_t_min) ) THEN
     4458         !!print *,'XXXwake, hw0(i), dthmin(i) ', hw0(i), dthmin(i)
     4459         !!print *,'XXXwake, 2.*sum_dth(i)/(hw0(i)*dthmin(i)) ', &
     4460         !!                  2.*sum_dth(i)/(hw0(i)*dthmin(i))
     4461         !!print *,'XXXwake, sum_half_dth(i), sum_dth(i) ', &
     4462         !!                  sum_half_dth(i), sum_dth(i)
     4463        IF (iflag_wk_check_trgl<=2 .and. ((hw0(i) < 1.) .or. (dthmin(i) >= -delta_t_min)) ) THEN
    44644464          wape2(i) = -1.
    4465           !! print *,'wake, rej 1'
     4465          !!  print *,'XXXwake, rej 1'
     4466        ELSE IF (iflag_wk_check_trgl==3 .and. ((hw0(i) < 1.) .or. (dthmin(i) >= dth(i,ktop(i)))) ) THEN
     4467          wape2(i) = -1.
     4468           !! print *,'XXXwake, rej 1'
    44664469        ELSE IF (iflag_wk_check_trgl==1.AND.abs(2.*sum_dth(i)/(hw0(i)*dthmin(i)) - 1.) > 0.5) THEN
    44674470          wape2(i) = -1.
    4468           !! print *,'wake, rej 2'
     4471           !! print *,'XXXwake, rej 2'
    44694472        ELSE IF (abs(sum_half_dth(i)) < 0.5*abs(sum_dth(i)) ) THEN
    44704473          wape2(i) = -1.
    4471           !! print *,'wake, rej 3'
     4474           !! print *,'XXXwake, rej 3'
    44724475        END IF
    44734476      END IF
    … …  
    47774780 RETURN
    47784781END SUBROUTINE wake2
     4782
     4783
     4784
     4785SUBROUTINE wake3(klon,klev,znatsurf, p, ph, pi, dtime, &
     4786                tb0, qb0, omgb, &
     4787                dtdwn, dqdwn, amdwn, amup, dta, dqa, wgen, &
     4788                sigd_con, Cin, &
     4789                deltatw, deltaqw, sigmaw, asigmaw, wdens, awdens, &                  ! state variables           
     4790                dth, hw, wape, fip, gfl, &
     4791                dtls, dqls, ktopw, omgbdth, dp_omgb, tx, qx, &
     4792!!!                dtke, dqke, omg, dp_deltomg, wkspread, cstar, &   ! changes in notation
     4793                d_deltat_dcv, d_deltaq_dcv, omg, dp_deltomg, wkspread, cstar, &
     4794                d_deltat_gw, &                                                      ! tendencies
     4795                d_deltatw2, d_deltaqw2, d_sigmaw2, d_asigmaw2, d_wdens2, d_awdens2)             ! tendencies
     4796
     4797
     4798  ! **************************************************************
     4799  ! *
     4800  ! WAKE                                                        *
     4801  ! retour a un Pupper fixe                                *
     4802  ! *
     4803  ! written by   :  GRANDPEIX Jean-Yves   09/03/2000            *
     4804  ! modified by :   ROEHRIG Romain        01/29/2007            *
     4805  ! **************************************************************
     4806
     4807
     4808  USE lmdz_wake_ini , ONLY : wake_ini
     4809  USE lmdz_wake_ini , ONLY : prt_level,epsim1,RG,RD
     4810  USE lmdz_wake_ini , ONLY : stark, wdens_ref, coefgw, alpk, wk_pupper
     4811  USE lmdz_wake_ini , ONLY : crep_upper, crep_sol, tau_cv, rzero, aa0, flag_wk_check_trgl
     4812  USE lmdz_wake_ini , ONLY : ok_bug_gfl
     4813  USE lmdz_wake_ini , ONLY : iflag_wk_act, iflag_wk_check_trgl, iflag_wk_pop_dyn, wdensinit, wdensthreshold
     4814  USE lmdz_wake_ini , ONLY : sigmad, hwmin, wapecut, cstart, sigmaw_max, dens_rate, epsilon_loc
     4815  USE lmdz_wake_ini , ONLY : iflag_wk_profile
     4816  USE lmdz_wake_ini , ONLY : smallestreal,wk_nsub
     4817
     4818
     4819  IMPLICIT NONE
     4820  ! ============================================================================
     4821
     4822
     4823  ! But : Decrire le comportement des poches froides apparaissant dans les
     4824  ! grands systemes convectifs, et fournir l'energie disponible pour
     4825  ! le declenchement de nouvelles colonnes convectives.
     4826
     4827  ! State variables :
     4828  ! deltatw    : temperature difference between wake and off-wake regions
     4829  ! deltaqw    : specific humidity difference between wake and off-wake regions
     4830  ! sigmaw     : fractional area covered by wakes.
     4831  ! asigmaw    : fractional area covered by active wakes.
     4832  ! wdens      : number of wakes per unit area
     4833  ! awdens     : number of active wakes per unit area
     4834
     4835  ! Variable de sortie :
     4836
     4837  ! wape : WAke Potential Energy
     4838  ! fip  : Front Incident Power (W/m2) - ALP
     4839  ! gfl  : Gust Front Length per unit area (m-1)
     4840  ! dtls : large scale temperature tendency due to wake
     4841  ! dqls : large scale humidity tendency due to wake
     4842  ! hw   : wake top hight (given by hw*deltatw(1)/2=wape)
     4843  ! dp_omgb : vertical gradient of large scale omega
     4844  ! awdens  : densite de poches actives
     4845  ! wdens   : densite de poches
     4846  ! omgbdth: flux of Delta_Theta transported by LS omega
     4847  ! d_deltat_dcv   : differential heating (wake - unpertubed)
     4848  ! d_deltat_dcv   : differential moistening (wake - unpertubed)
     4849  ! omg    : Delta_omg =vertical velocity diff. wake-undist. (Pa/s)
     4850  ! dp_deltomg  : vertical gradient of omg (s-1)
     4851  ! wkspread  : spreading term in d_t_wake and d_q_wake
     4852  ! deltatw     : updated temperature difference (T_w-T_u).
     4853  ! deltaqw     : updated humidity difference (q_w-q_u).
     4854  ! sigmaw      : updated wake fractional area.
     4855  ! asigmaw     : updated active wake fractional area.
     4856  ! d_deltat_gw : delta T tendency due to GW
     4857
     4858  ! Variables d'entree :
     4859
     4860  ! aire : aire de la maille
     4861  ! tb0  : horizontal average of temperature  (K)
     4862  ! qb0  : horizontal average of humidity   (kg/kg)
     4863  ! omgb : vitesse verticale moyenne sur la maille (Pa/s)
     4864  ! dtdwn: source de chaleur due aux descentes (K/s)
     4865  ! dqdwn: source d'humidite due aux descentes (kg/kg/s)
     4866  ! dta  : source de chaleur due courants satures et detrain  (K/s)
     4867  ! dqa  : source d'humidite due aux courants satures et detra (kg/kg/s)
     4868  ! wgen : number of wakes generated per unit area and per sec (/m^2/s)
     4869  ! amdwn: flux de masse total des descentes, par unite de
     4870  !        surface de la maille (kg/m2/s)
     4871  ! amup : flux de masse total des ascendances, par unite de
     4872  !        surface de la maille (kg/m2/s)
     4873  ! sigd_con:
     4874  ! Cin  : convective inhibition
     4875  ! p    : pressions aux milieux des couches (Pa)
     4876  ! ph   : pressions aux interfaces (Pa)
     4877  ! pi  : (p/p_0)**kapa (adim)
     4878  ! dtime: increment temporel (s)
     4879
     4880  ! Variables internes :
     4881
     4882  ! rho  : mean density at P levels
     4883  ! rhoh : mean density at Ph levels
     4884  ! tb   : mean temperature | may change within
     4885  ! qb   : mean humidity    | sub-time-stepping
     4886  ! thb  : mean potential temperature
     4887  ! thx  : potential temperature in (x) area
     4888  ! tx   : temperature  in (x) area
     4889  ! qx   : humidity in (x) area
     4890  ! dp_omgb: vertical gradient og LS omega
     4891  ! omgbw  : wake average vertical omega
     4892  ! dp_omgbw: vertical gradient of omgbw
     4893  ! omgbdq : flux of Delta_q transported by LS omega
     4894  ! dth  : potential temperature diff. wake-undist.
     4895  ! th1  : first pot. temp. for vertical advection (=thx)
     4896  ! th2  : second pot. temp. for vertical advection (=thw)
     4897  ! q1   : first humidity for vertical advection
     4898  ! q2   : second humidity for vertical advection
     4899  ! d_deltatw   : redistribution term for deltatw
     4900  ! d_deltaqw   : redistribution term for deltaqw
     4901  ! deltatw0   : initial deltatw
     4902  ! deltaqw0   : initial deltaqw
     4903  ! hw0    : wake top hight (defined as the altitude at which deltatw=0)
     4904  ! amflux : horizontal mass flux through wake boundary
     4905  ! wdens_ref: initial number of wakes per unit area (3D) or per
     4906  !            unit length (2D), at the beginning of each time step
     4907  ! Tgw    : 1 sur la periode de onde de gravite
     4908  ! Cgw    : vitesse de propagation de onde de gravite
     4909  ! LL     : distance between 2 wakes
     4910  ! Tgen   : 1 sur le temps caracteristique d'amortissement par les naissances de poches
     4911
     4912  ! -------------------------------------------------------------------------
     4913  ! Declaration de variables
     4914  ! -------------------------------------------------------------------------
     4915
     4916
     4917  ! Arguments en entree
     4918  ! --------------------
     4919
     4920  INTEGER,                          INTENT(IN)          :: klon,klev
     4921  INTEGER, DIMENSION (klon),        INTENT(IN)          :: znatsurf
     4922  REAL, DIMENSION (klon, klev),     INTENT(IN)          :: p, pi
     4923  REAL, DIMENSION (klon, klev+1),   INTENT(IN)          :: ph
     4924  REAL, DIMENSION (klon, klev),     INTENT(IN)          :: omgb
     4925  REAL,                             INTENT(IN)          :: dtime
     4926  REAL, DIMENSION (klon, klev),     INTENT(IN)          :: tb0, qb0
     4927  REAL, DIMENSION (klon, klev),     INTENT(IN)          :: dtdwn, dqdwn
     4928  REAL, DIMENSION (klon, klev),     INTENT(IN)          :: amdwn, amup
     4929  REAL, DIMENSION (klon, klev),     INTENT(IN)          :: dta, dqa
     4930  REAL, DIMENSION (klon),           INTENT(IN)          :: wgen
     4931  REAL, DIMENSION (klon),           INTENT(IN)          :: sigd_con
     4932  REAL, DIMENSION (klon),           INTENT(IN)          :: Cin
     4933
     4934  !
     4935  ! Input/Output
     4936  ! State variables
     4937  REAL, DIMENSION (klon, klev),     INTENT(INOUT)       :: deltatw, deltaqw
     4938  REAL, DIMENSION (klon),           INTENT(INOUT)       :: sigmaw
     4939  REAL, DIMENSION (klon),           INTENT(INOUT)       :: asigmaw
     4940  REAL, DIMENSION (klon),           INTENT(INOUT)       :: wdens
     4941  REAL, DIMENSION (klon),           INTENT(INOUT)       :: awdens
     4942
     4943  ! Sorties
     4944  ! --------
     4945
     4946  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: dth
     4947  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: tx, qx
     4948  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: dtls, dqls
     4949!!  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: dtke, dqke
     4950  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: d_deltat_dcv, d_deltaq_dcv
     4951  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: wkspread    !  unused (jyg)
     4952  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: omgbdth, omg
     4953  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: dp_omgb, dp_deltomg
     4954  REAL, DIMENSION (klon),           INTENT(OUT)         :: hw, wape, fip, gfl, cstar
     4955  INTEGER, DIMENSION (klon),        INTENT(OUT)         :: ktopw
     4956  ! Tendencies of state variables (2 is appended to the names of fields which are the cumul of fields
     4957  !                                 computed at each sub-timestep; e.g. d_wdens2 is the cumul of d_wdens)
     4958  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: d_deltat_gw
     4959  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: d_deltatw2, d_deltaqw2
     4960  REAL, DIMENSION (klon),           INTENT(OUT)         :: d_sigmaw2, d_asigmaw2, d_wdens2, d_awdens2
     4961
     4962  ! Variables internes
     4963  ! -------------------
     4964
     4965  ! Variables a fixer
     4966
     4967  REAL                                                  :: delta_t_min
     4968  REAL                                                  :: dtimesub
     4969  REAL                                                  :: wdens0
     4970  ! IM 080208
     4971  LOGICAL, DIMENSION (klon)                             :: gwake
     4972
     4973  ! Variables de sauvegarde
     4974  REAL, DIMENSION (klon, klev)                          :: deltatw0
     4975  REAL, DIMENSION (klon, klev)                          :: deltaqw0
     4976  REAL, DIMENSION (klon, klev)                          :: tb, qb
     4977
     4978  ! Variables liees a la dynamique de population 1
     4979  REAL, DIMENSION(klon)                                 :: act
     4980  REAL, DIMENSION(klon)                                 :: rad_wk, tau_wk_inv
     4981  REAL, DIMENSION(klon)                                 :: f_shear
     4982  REAL, DIMENSION(klon)                                 :: drdt
     4983 
     4984  ! Variables liees a la dynamique de population 2
     4985  REAL, DIMENSION(klon)                                 :: cont_fact 
     4986 
     4987  ! Variables liees a la dynamique de population 3
     4988  REAL, DIMENSION(klon)                                 :: arad_wk, irad_wk
     4989 
     4990!!  REAL, DIMENSION(klon)                                 :: d_sig_gen, d_sig_death, d_sig_col
     4991  REAL, DIMENSION(klon)                                 :: wape1_act, wape2_act
     4992  LOGICAL, DIMENSION (klon)                             :: kill_wake
     4993  REAL                                                  :: drdt_pos
     4994  REAL                                                  :: tau_wk_inv_min
     4995  ! Some components of the tendencies of state variables 
     4996  REAL, DIMENSION (klon)                                :: d_sig_gen2, d_sig_death2, d_sig_col2, d_sig_spread2, d_sig_bnd2
     4997  REAL, DIMENSION (klon)                                :: d_asig_death2, d_asig_aicol2, d_asig_iicol2, d_asig_spread2, d_asig_bnd2
     4998  REAL, DIMENSION (klon)                                :: d_dens_gen2, d_dens_death2, d_dens_col2, d_dens_bnd2
     4999  REAL, DIMENSION (klon)                                :: d_adens_death2, d_adens_icol2, d_adens_acol2, d_adens_bnd2
     5000
     5001  ! Variables pour les GW
     5002  REAL, DIMENSION (klon)                                :: ll
     5003  REAL, DIMENSION (klon, klev)                          :: n2
     5004  REAL, DIMENSION (klon, klev)                          :: cgw
     5005  REAL, DIMENSION (klon, klev)                          :: tgw
     5006  REAL, DIMENSION (klon, klev)                          :: tgen
     5007
     5008  ! Variables liees au calcul de hw
     5009  REAL, DIMENSION (klon)                                :: ptop
     5010  REAL, DIMENSION (klon)                                :: sum_dth
     5011  REAL, DIMENSION (klon)                                :: dthmin
     5012  REAL, DIMENSION (klon)                                :: z, dz, hw0
     5013  INTEGER, DIMENSION (klon)                             :: ktop, kupper
     5014
     5015  ! Variables liees au test de la forme triangulaire du profil de Delta_theta
     5016  REAL, DIMENSION (klon)                                :: sum_half_dth
     5017  REAL, DIMENSION (klon)                                :: dz_half
     5018
     5019  ! Sub-timestep tendencies and related variables
     5020  REAL, DIMENSION (klon, klev)                          :: d_deltatw, d_deltaqw
     5021  REAL, DIMENSION (klon, klev)                          :: d_deltat_dadv, d_deltaq_dadv
     5022  REAL, DIMENSION (klon, klev)                          :: d_deltat_lsadv, d_deltaq_lsadv
     5023  REAL, DIMENSION (klon, klev)                          :: d_deltat_entrp
     5024  REAL, DIMENSION (klon, klev)                          :: d_deltat_spread, d_deltaq_spread
     5025
     5026  REAL, DIMENSION (klon, klev)                          :: d_tb, d_qb
     5027  REAL, DIMENSION (klon, klev)                          :: d_tb_dadv, d_qb_dadv
     5028  REAL, DIMENSION (klon, klev)                          :: d_tb_spread, d_qb_spread
     5029
     5030  REAL, DIMENSION (klon)                                :: d_wdens, d_awdens, d_sigmaw, d_asigmaw
     5031  REAL, DIMENSION (klon)                                :: d_sig_gen, d_sig_death, d_sig_col, d_sig_spread, d_sig_bnd
     5032  REAL, DIMENSION (klon)                                :: d_asig_death, d_asig_aicol, d_asig_iicol, d_asig_spread, d_asig_bnd
     5033  REAL, DIMENSION (klon)                                :: d_dens_gen, d_dens_death, d_dens_col, d_dens_bnd
     5034  REAL, DIMENSION (klon)                                :: d_adens_death, d_adens_icol, d_adens_acol, d_adens_bnd
     5035  REAL, DIMENSION (klon)                                :: agfl              !! gust front length of active wakes
     5036                                                                             !!  per unit area
     5037  REAL, DIMENSION (klon)                                :: alpha, alpha_tot
     5038  REAL, DIMENSION (klon)                                :: q0_min, q1_min
     5039  LOGICAL, DIMENSION (klon)                             :: wk_adv, ok_qx_qw
     5040
     5041
     5042  ! Autres variables internes
     5043  INTEGER                                               ::isubstep, k, i, igout
     5044
     5045  REAL                                                  :: wdensmin
     5046
     5047  REAL                                                  :: sigmaw_targ
     5048  REAL                                                  :: wdens_targ
     5049  REAL                                                  :: d_sigmaw_targ
     5050  REAL                                                  :: d_wdens_targ
     5051
     5052  REAL, DIMENSION (klon)                                :: dsigspread  !rate of change of sigmaw due to spreading
     5053
     5054  REAL, DIMENSION (klon)                                :: sum_thx, sum_tx, sum_qx, sum_thvx
     5055  REAL, DIMENSION (klon)                                :: sum_dq
     5056  REAL, DIMENSION (klon)                                :: sum_dtdwn, sum_dqdwn
     5057  REAL, DIMENSION (klon)                                :: av_thx, av_tx, av_qx, av_thvx
     5058  REAL, DIMENSION (klon)                                :: av_dth, av_dq
     5059  REAL, DIMENSION (klon)                                :: av_dtdwn, av_dqdwn
     5060
     5061  REAL, DIMENSION (klon, klev)                          :: rho
     5062  REAL, DIMENSION (klon, klev+1)                        :: rhoh
     5063  REAL, DIMENSION (klon, klev)                          :: zh
     5064  REAL, DIMENSION (klon, klev+1)                        :: zhh
     5065
     5066  REAL, DIMENSION (klon, klev)                          :: thb, thx
     5067
     5068  REAL, DIMENSION (klon, klev)                          :: omgbw
     5069  REAL, DIMENSION (klon)                                :: pupper
     5070  REAL, DIMENSION (klon)                                :: omgtop
     5071  REAL, DIMENSION (klon, klev)                          :: dp_omgbw
     5072  REAL, DIMENSION (klon)                                :: ztop, dztop
     5073  REAL, DIMENSION (klon, klev)                          :: alpha_up
     5074
     5075  REAL, DIMENSION (klon)                                :: rre1, rre2
     5076  REAL                                                  :: rrd1, rrd2
     5077  REAL, DIMENSION (klon, klev)                          :: th1, th2, q1, q2
     5078  REAL, DIMENSION (klon, klev)                          :: d_th1, d_th2, d_dth
     5079  REAL, DIMENSION (klon, klev)                          :: d_q1, d_q2, d_dq
     5080  REAL, DIMENSION (klon, klev)                          :: omgbdq
     5081
     5082  REAL, DIMENSION (klon)                                :: wape2, cstar2, heff
     5083                                                       
     5084  REAL, DIMENSION (klon, klev)                          :: crep
     5085                                                       
     5086  REAL, DIMENSION (klon, klev)                          :: ppi
     5087
     5088  ! cc nrlmd
     5089  REAL, DIMENSION (klon)                                :: death_rate
     5090!!  REAL, DIMENSION (klon)                                :: nat_rate
     5091  REAL, DIMENSION (klon, klev)                          :: entr   ! total entrainment into wakes (spread + birth)
     5092  REAL, DIMENSION (klon, klev)                          :: entr_p ! entrainment into wakes (due to births)
     5093  REAL, DIMENSION (klon, klev)                          :: detr   ! detrainment from wakes (due to deaths)
     5094
     5095  REAL, DIMENSION(klon)                                 :: sigmaw_in, asigmaw_in ! pour les prints
     5096  REAL, DIMENSION(klon)                                 :: wdens_in, awdens_in   ! pour les prints
     5097
     5098  !!-- variables liees au nouveau calcul de ptop et hw
     5099  REAL, DIMENSION (klon, klev)                          :: int_dth
     5100  REAL, DIMENSION (klon, klev)                          :: zzz, dzzz
     5101  REAL                                                  :: epsil
     5102  REAL, DIMENSION (klon)                                :: ptop1
     5103  INTEGER, DIMENSION (klon)                             :: ktop1
     5104  REAL, DIMENSION (klon)                                :: omega
     5105  REAL, DIMENSION (klon)                                :: h_zzz
     5106
     5107  !! Bidouilles
     5108  REAL                                                  :: iwkadv
     5109  REAL                                                  :: iokqxqw
     5110  REAL                                                  :: zzzz,zzzzm1
     5111
     5112
     5113!print*,'WAKE LJYFz'
     5114
     5115  ! -------------------------------------------------------------------------
     5116  ! Initialisations
     5117  ! -------------------------------------------------------------------------
     5118  ! ALON = 3.e5
     5119  ! alon = 1.E6
     5120
     5121  !  Provisionnal; to be suppressed when f_shear is parameterized
     5122  f_shear(:) = 1.       ! 0. for strong shear, 1. for weak shear
     5123
     5124  ! Configuration de coefgw,stark,wdens (22/02/06 by YU Jingmei)
     5125
     5126  ! coefgw : Coefficient pour les ondes de gravite
     5127  ! stark : Coefficient k dans Cstar=k*sqrt(2*WAPE)
     5128  ! wdens : Densite surfacique de poche froide
     5129  ! -------------------------------------------------------------------------
     5130
     5131  ! cc nrlmd      coefgw=10
     5132  ! coefgw=1
     5133  ! wdens0 = 1.0/(alon**2)
     5134  ! cc nrlmd      wdens = 1.0/(alon**2)
     5135  ! cc nrlmd      stark = 0.50
     5136  ! CRtest
     5137  ! cc nrlmd      alpk=0.1
     5138  ! alpk = 1.0
     5139  ! alpk = 0.5
     5140  ! alpk = 0.05
     5141!
     5142 igout = klon/2+1/klon 
     5143!
     5144!   sub-time-stepping parameters
     5145  dtimesub = dtime/wk_nsub
     5146!
     5147 IF (iflag_wk_pop_dyn == 0) THEN
     5148  ! Initialisation de toutes des densites a wdens_ref.
     5149  ! Les densites peuvent evoluer si les poches debordent
     5150  ! (voir au tout debut de la boucle sur les substeps)
     5151  !jyg<
     5152  !!  wdens(:) = wdens_ref
     5153   DO i = 1,klon
     5154     wdens(i) = wdens_ref(znatsurf(i)+1)
     5155   ENDDO
     5156  !>jyg
     5157 ENDIF  ! (iflag_wk_pop_dyn == 0)
     5158!
     5159 IF (iflag_wk_pop_dyn >=1) THEN
     5160   IF (iflag_wk_pop_dyn == 3) THEN
     5161     wdensmin = wdensthreshold
     5162   ELSE
     5163     wdensmin = wdensinit
     5164   ENDIF
     5165 ENDIF
     5166
     5167  ! print*,'stark',stark
     5168  ! print*,'alpk',alpk
     5169  ! print*,'wdens',wdens
     5170  ! print*,'coefgw',coefgw
     5171  ! cc
     5172  ! Minimum value for |T_wake - T_undist|. Used for wake top definition
     5173  ! -------------------------------------------------------------------------
     5174
     5175   delta_t_min = 0.2
     5176
     5177  ! 1. - Save initial values, initialize tendencies, initialize output fields
     5178  ! ------------------------------------------------------------------------
     5179
     5180!jyg<
     5181!!  DO k = 1, klev
     5182!!    DO i = 1, klon
     5183!!      ppi(i, k) = pi(i, k)
     5184!!      deltatw0(i, k) = deltatw(i, k)
     5185!!      deltaqw0(i, k) = deltaqw(i, k)
     5186!!      tb(i, k) = tb0(i, k)
     5187!!      qb(i, k) = qb0(i, k)
     5188!!      dtls(i, k) = 0.
     5189!!      dqls(i, k) = 0.
     5190!!      d_deltat_gw(i, k) = 0.
     5191!!      d_tb(i, k) = 0.
     5192!!      d_qb(i, k) = 0.
     5193!!      d_deltatw(i, k) = 0.
     5194!!      d_deltaqw(i, k) = 0.
     5195!!      ! IM 060508 beg
     5196!!      d_deltatw2(i, k) = 0.
     5197!!      d_deltaqw2(i, k) = 0.
     5198!!      ! IM 060508 end
     5199!!    END DO
     5200!!  END DO
     5201      ppi(:,:) = pi(:,:)
     5202      deltatw0(:,:) = deltatw(:,:)
     5203      deltaqw0(:,:) = deltaqw(:,:)
     5204      tb(:,:) = tb0(:,:)
     5205      qb(:,:) = qb0(:,:)
     5206      dtls(:,:) = 0.
     5207      dqls(:,:) = 0.
     5208      d_deltat_gw(:,:) = 0.
     5209
     5210      detr(:,:) = 0.
     5211      entr(:,:) = 0.
     5212      entr_p(:,:) = 0.
     5213
     5214      th1(:,:) = 0.
     5215      th2(:,:) = 0.
     5216      q1(:,:) = 0.
     5217      q2(:,:) = 0.
     5218
     5219      d_tb(:,:) = 0.
     5220      d_tb_dadv(:,:) = 0.
     5221      d_tb_spread(:,:) = 0.
     5222
     5223      d_qb(:,:) = 0.
     5224      d_qb_dadv(:,:) = 0.
     5225      d_qb_spread(:,:) = 0.
     5226
     5227      d_deltatw(:,:) = 0.
     5228      d_deltat_dadv  (:,:) = 0.
     5229      d_deltat_lsadv (:,:) = 0.
     5230      d_deltat_dcv   (:,:) = 0.
     5231      d_deltat_spread(:,:) = 0.
     5232
     5233      d_deltaqw(:,:) = 0.
     5234      d_deltaq_dadv(:,:) = 0.
     5235      d_deltaq_lsadv(:,:) = 0.
     5236      d_deltaq_dcv(:,:) = 0.
     5237      d_deltaq_spread(:,:) = 0.
     5238
     5239      d_deltatw2(:,:) = 0.
     5240      d_deltaqw2(:,:) = 0.
     5241
     5242      d_sig_gen2(:)   = 0.
     5243      d_sig_death2(:) = 0.
     5244      d_sig_col2(:)   = 0.
     5245      d_sig_spread2(:)= 0.
     5246      d_asig_death2(:) = 0.
     5247      d_asig_iicol2(:) = 0.
     5248      d_asig_aicol2(:) = 0.
     5249      d_asig_spread2(:)= 0.
     5250      d_asig_bnd2(:) = 0.
     5251      d_asigmaw2(:) = 0.
     5252!
     5253      d_dens_gen2(:)   = 0.
     5254      d_dens_death2(:) = 0.
     5255      d_dens_col2(:)   = 0.
     5256      d_dens_bnd2(:) = 0.
     5257      d_wdens2(:) = 0.
     5258      d_adens_bnd2(:) = 0.
     5259      d_awdens2(:) = 0.
     5260      d_adens_death2(:) = 0.
     5261      d_adens_icol2(:) = 0.
     5262      d_adens_acol2(:) = 0.
     5263
     5264      IF (iflag_wk_act == 0) THEN
     5265        act(:) = 0.
     5266      ELSEIF (iflag_wk_act == 1) THEN
     5267        act(:) = 1.
     5268      ENDIF
     5269
     5270!!  DO i = 1, klon
     5271!!   sigmaw_in(i) = sigmaw(i)
     5272!!  END DO
     5273   sigmaw_in(:)  = sigmaw(:)
     5274   asigmaw_in(:) = asigmaw(:)
     5275!>jyg
     5276!
     5277  IF (iflag_wk_pop_dyn >= 1) THEN
     5278    awdens_in(:) = awdens(:)
     5279    wdens_in(:) = wdens(:)
     5280!!    wdens(:) = wdens(:) + wgen(:)*dtime
     5281!!    d_wdens2(:) = wgen(:)*dtime
     5282!!  ELSE
     5283  ENDIF  ! (iflag_wk_pop_dyn >= 1)
     5284
     5285
     5286  ! sigmaw1=sigmaw
     5287  ! IF (sigd_con.GT.sigmaw1) THEN
     5288  ! print*, 'sigmaw,sigd_con', sigmaw, sigd_con
     5289  ! ENDIF
     5290  IF (iflag_wk_pop_dyn >= 1) THEN
     5291    DO i = 1, klon
     5292      d_dens_gen2(i)   = 0.
     5293      d_dens_death2(i) = 0.
     5294      d_dens_col2(i)   = 0.
     5295      d_awdens2(i) = 0.
     5296      IF (wdens(i) < wdensthreshold) THEN
     5297  !!      wdens_targ = max(wdens(i),wdensmin)
     5298        wdens_targ = max(wdens(i),wdensinit)
     5299        d_dens_bnd2(i) = wdens_targ - wdens(i)
     5300        d_wdens2(i) = wdens_targ - wdens(i)
     5301        wdens(i) = wdens_targ
     5302      ELSE
     5303        d_dens_bnd2(i) = 0.
     5304        d_wdens2(i) = 0.
     5305      ENDIF  !! (wdens(i) < wdensthreshold)
     5306    END DO
     5307    IF (iflag_wk_pop_dyn >= 2) THEN
     5308      DO i = 1, klon 
     5309        IF (awdens(i) < wdensthreshold) THEN
     5310!!          wdens_targ = min(max(awdens(i),wdensmin),wdens(i))
     5311            wdens_targ = min(max(awdens(i),wdensinit),wdens(i))
     5312            d_adens_bnd2(i) = wdens_targ - awdens(i)
     5313            d_awdens2(i) = wdens_targ - awdens(i)
     5314            awdens(i) = wdens_targ
     5315        ELSE
     5316            wdens_targ = min(awdens(i), wdens(i))
     5317            d_adens_bnd2(i) = wdens_targ - awdens(i)
     5318            d_awdens2(i) = wdens_targ - awdens(i)
     5319            awdens(i) = wdens_targ
     5320        ENDIF
     5321      END DO
     5322    ENDIF ! (iflag_wk_pop_dyn >= 2)
     5323  ELSE 
     5324    DO i = 1, klon
     5325      d_awdens2(i) = 0.
     5326      d_wdens2(i) = 0.
     5327    END DO
     5328  ENDIF  ! (iflag_wk_pop_dyn >= 1)
     5329!
     5330  DO i = 1, klon
     5331    sigmaw_targ = min(max(sigmaw(i), sigmad),0.99)
     5332    d_sig_bnd2(i) = sigmaw_targ - sigmaw(i)
     5333    d_sigmaw2(i) = sigmaw_targ - sigmaw(i)
     5334    sigmaw(i) = sigmaw_targ
     5335  END DO
     5336!
     5337  IF (iflag_wk_pop_dyn == 3) THEN
     5338     DO i = 1, klon 
     5339        IF ((wdens(i)-awdens(i)) <= smallestreal) THEN
     5340          sigmaw_targ = sigmaw(i)
     5341        ELSE
     5342          sigmaw_targ = min(max(asigmaw(i),sigmad),sigmaw(i))
     5343        ENDIF
     5344        d_asig_bnd2(i) = sigmaw_targ - asigmaw(i)
     5345        d_asigmaw2(i) = sigmaw_targ - asigmaw(i)
     5346        asigmaw(i) = sigmaw_targ
     5347     END DO
     5348  ENDIF ! (iflag_wk_pop_dyn == 3)
     5349
     5350  wape(:) = 0.
     5351  wape2(:) = 0.
     5352  d_sigmaw(:) = 0.
     5353  d_asigmaw(:) = 0.
     5354  ktopw(:) = 0
     5355!
     5356!<jyg
     5357dth(:,:) = 0.
     5358tx(:,:) = 0.
     5359qx(:,:) = 0.
     5360d_deltat_dcv(:,:) = 0.
     5361d_deltaq_dcv(:,:) = 0.
     5362wkspread(:,:) = 0.
     5363omgbdth(:,:) = 0.
     5364omg(:,:) = 0.
     5365dp_omgb(:,:) = 0.
     5366dp_deltomg(:,:) = 0.
     5367tgen(:,:) = 0.
     5368hw(:) = 0.
     5369wape(:) = 0.
     5370fip(:) = 0.
     5371gfl(:) = 0.
     5372cstar(:) = 0.
     5373ktopw(:) = 0
     5374!
     5375!  Vertical advection local variables
     5376omgbw(:,:) = 0.
     5377omgtop(:) = 0
     5378dp_omgbw(:,:) = 0.
     5379omgbdq(:,:) = 0.
     5380
     5381!>jyg
     5382!
     5383  IF (prt_level>=10) THEN
     5384    PRINT *, 'wake-1, sigmaw(igout) ', sigmaw(igout)
     5385    PRINT *, 'wake-1, deltatw(igout,k) ', (k,deltatw(igout,k), k=1,klev)
     5386    PRINT *, 'wake-1, deltaqw(igout,k) ', (k,deltaqw(igout,k), k=1,klev)
     5387    PRINT *, 'wake-1, dowwdraughts, amdwn(igout,k) ', (k,amdwn(igout,k), k=1,klev)
     5388    PRINT *, 'wake-1, dowwdraughts, dtdwn(igout,k) ', (k,dtdwn(igout,k), k=1,klev)
     5389    PRINT *, 'wake-1, dowwdraughts, dqdwn(igout,k) ', (k,dqdwn(igout,k), k=1,klev)
     5390    PRINT *, 'wake-1, updraughts, amup(igout,k) ', (k,amup(igout,k), k=1,klev)
     5391    PRINT *, 'wake-1, updraughts, dta(igout,k) ', (k,dta(igout,k), k=1,klev)
     5392    PRINT *, 'wake-1, updraughts, dqa(igout,k) ', (k,dqa(igout,k), k=1,klev)
     5393  ENDIF
     5394
     5395  ! 2. - Prognostic part
     5396  ! --------------------
     5397
     5398
     5399  ! 2.1 - Undisturbed area and Wake integrals
     5400  ! ---------------------------------------------------------
     5401
     5402  DO i = 1, klon
     5403    z(i) = 0.
     5404    ktop(i) = 0
     5405    kupper(i) = 0
     5406    sum_thx(i) = 0.
     5407    sum_tx(i) = 0.
     5408    sum_qx(i) = 0.
     5409    sum_thvx(i) = 0.
     5410    sum_dth(i) = 0.
     5411    sum_dq(i) = 0.
     5412    sum_dtdwn(i) = 0.
     5413    sum_dqdwn(i) = 0.
     5414
     5415    av_thx(i) = 0.
     5416    av_tx(i) = 0.
     5417    av_qx(i) = 0.
     5418    av_thvx(i) = 0.
     5419    av_dth(i) = 0.
     5420    av_dq(i) = 0.
     5421    av_dtdwn(i) = 0.
     5422    av_dqdwn(i) = 0.
     5423  END DO
     5424
     5425  ! Distance between wakes
     5426  DO i = 1, klon
     5427    ll(i) = (1-sqrt(sigmaw(i)))/sqrt(wdens(i))
     5428  END DO
     5429  ! Potential temperatures and humidity
     5430  ! ----------------------------------------------------------
     5431  DO k = 1, klev
     5432    DO i = 1, klon
     5433      ! write(*,*)'wake 1',i,k,RD,tb(i,k)
     5434      rho(i, k) = p(i, k)/(RD*tb(i,k))
     5435      ! write(*,*)'wake 2',rho(i,k)
     5436      IF (k==1) THEN
     5437        ! write(*,*)'wake 3',i,k,rd,tb(i,k)
     5438        rhoh(i, k) = ph(i, k)/(RD*tb(i,k))
     5439        ! write(*,*)'wake 4',i,k,rd,tb(i,k)
     5440        zhh(i, k) = 0
     5441      ELSE
     5442        ! write(*,*)'wake 5',rd,(tb(i,k)+tb(i,k-1))
     5443        rhoh(i, k) = ph(i, k)*2./(RD*(tb(i,k)+tb(i,k-1)))
     5444        ! write(*,*)'wake 6',(-rhoh(i,k)*RG)+zhh(i,k-1)
     5445        zhh(i, k) = (ph(i,k)-ph(i,k-1))/(-rhoh(i,k)*RG) + zhh(i, k-1)
     5446      END IF
     5447      ! write(*,*)'wake 7',ppi(i,k)
     5448      thb(i, k) = tb(i, k)/ppi(i, k)
     5449      thx(i, k) = (tb(i,k)-deltatw(i,k)*sigmaw(i))/ppi(i, k)
     5450      tx(i, k) = tb(i, k) - deltatw(i, k)*sigmaw(i)
     5451      qx(i, k) = qb(i, k) - deltaqw(i, k)*sigmaw(i)
     5452      ! write(*,*)'wake 8',(RD*(tb(i,k)+deltatw(i,k)))
     5453      dth(i, k) = deltatw(i, k)/ppi(i, k)
     5454    END DO
     5455  END DO
     5456
     5457  DO k = 1, klev - 1
     5458    DO i = 1, klon
     5459      IF (k==1) THEN
     5460        n2(i, k) = 0
     5461      ELSE
     5462        n2(i, k) = amax1(0., -RG**2/thb(i,k)*rho(i,k)*(thb(i,k+1)-thb(i,k-1))/ &
     5463                             (p(i,k+1)-p(i,k-1)))
     5464      END IF
     5465      zh(i, k) = (zhh(i,k)+zhh(i,k+1))/2
     5466
     5467      cgw(i, k) = sqrt(n2(i,k))*zh(i, k)
     5468      tgw(i, k) = coefgw*cgw(i, k)/ll(i)
     5469    END DO
     5470  END DO
     5471
     5472  DO i = 1, klon
     5473    n2(i, klev) = 0
     5474    zh(i, klev) = 0
     5475    cgw(i, klev) = 0
     5476    tgw(i, klev) = 0
     5477  END DO
     5478
     5479 
     5480  ! Choose an integration bound well above wake top
     5481  ! -----------------------------------------------------------------
     5482
     5483  ! Determine Wake top pressure (Ptop) from buoyancy integral
     5484  ! --------------------------------------------------------
     5485
     5486   Do i=1, klon
     5487       wk_adv(i) = .True.
     5488   Enddo
     5489   Call pkupper (klon, klev, ptop, ph, p, pupper, kupper, &
     5490                    dth, hw0, rho, delta_t_min, &
     5491                    ktop, wk_adv, h_zzz, ptop1, ktop1)
     5492 
     5493   !!print'("pkupper APPEL ",7i6)',0,int(ptop/100.),int(ptop1/100.),int(pupper/100.),ktop,ktop1,kupper
     5494   
     5495   IF (prt_level>=10) THEN
     5496     PRINT *, 'wake-3, ktop(igout), kupper(igout) ', ktop(igout), kupper(igout)
     5497  ENDIF
     5498
     5499  ! -5/ Set deltatw & deltaqw to 0 above kupper
     5500
     5501  DO k = 1, klev
     5502    DO i = 1, klon
     5503      IF (k>=kupper(i)) THEN
     5504        deltatw(i, k) = 0.
     5505        deltaqw(i, k) = 0.
     5506        d_deltatw2(i,k) = -deltatw0(i,k)
     5507        d_deltaqw2(i,k) = -deltaqw0(i,k)
     5508      END IF
     5509    END DO
     5510  END DO
     5511
     5512
     5513  ! Vertical gradient of LS omega
     5514
     5515  DO k = 1, klev
     5516    DO i = 1, klon
     5517      IF (k<=kupper(i)) THEN
     5518        dp_omgb(i, k) = (omgb(i,k+1)-omgb(i,k))/(ph(i,k+1)-ph(i,k))
     5519      END IF
     5520    END DO
     5521  END DO
     5522
     5523  ! Integrals (and wake top level number)
     5524  ! --------------------------------------
     5525
     5526  ! Initialize sum_thvx to 1st level virt. pot. temp.
     5527
     5528  DO i = 1, klon
     5529    z(i) = 1.
     5530    dz(i) = 1.
     5531    sum_thvx(i) = thx(i, 1)*(1.+epsim1*qx(i,1))*dz(i)
     5532    sum_dth(i) = 0.
     5533  END DO
     5534
     5535  DO k = 1, klev
     5536    DO i = 1, klon
     5537      dz(i) = -(amax1(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*RG)
     5538      IF (dz(i)>0) THEN
     5539              ! LJYF : ecriture pas sympa avec un tableau z(i) qui n'est pas utilise come tableau
     5540        z(i) = z(i) + dz(i)
     5541        sum_thx(i) = sum_thx(i) + thx(i, k)*dz(i)
     5542        sum_tx(i) = sum_tx(i) + tx(i, k)*dz(i)
     5543        sum_qx(i) = sum_qx(i) + qx(i, k)*dz(i)
     5544        sum_thvx(i) = sum_thvx(i) + thx(i, k)*(1.+epsim1*qx(i,k))*dz(i)
     5545        sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i)
     5546        sum_dq(i) = sum_dq(i) + deltaqw(i, k)*dz(i)
     5547        sum_dtdwn(i) = sum_dtdwn(i) + dtdwn(i, k)*dz(i)
     5548        sum_dqdwn(i) = sum_dqdwn(i) + dqdwn(i, k)*dz(i)
     5549      END IF
     5550    END DO
     5551  END DO
     5552
     5553  DO i = 1, klon
     5554    hw0(i) = z(i)
     5555  END DO
     5556
     5557
     5558  ! 2.1 - WAPE and mean forcing computation
     5559  ! ---------------------------------------
     5560
     5561  ! ---------------------------------------
     5562
     5563  ! Means
     5564
     5565  DO i = 1, klon
     5566    av_thx(i) = sum_thx(i)/hw0(i)
     5567    av_tx(i) = sum_tx(i)/hw0(i)
     5568    av_qx(i) = sum_qx(i)/hw0(i)
     5569    av_thvx(i) = sum_thvx(i)/hw0(i)
     5570    ! av_thve = sum_thve/hw0
     5571    av_dth(i) = sum_dth(i)/hw0(i)
     5572    av_dq(i) = sum_dq(i)/hw0(i)
     5573    av_dtdwn(i) = sum_dtdwn(i)/hw0(i)
     5574    av_dqdwn(i) = sum_dqdwn(i)/hw0(i)
     5575
     5576    wape(i) = -RG*hw0(i)*(av_dth(i)+ &
     5577        epsim1*(av_thx(i)*av_dq(i)+av_dth(i)*av_qx(i)+av_dth(i)*av_dq(i)))/av_thvx(i)
     5578
     5579  END DO
     5580IF (CPPKEY_IOPHYS_WK) THEN
     5581  IF (.not.phys_sub) CALL iophys_ecrit('wape_a',1,'wape_a','J/kg',wape)
     5582END IF
     5583
     5584  ! 2.2 Prognostic variable update
     5585  ! ------------------------------
     5586
     5587  ! Filter out bad wakes
     5588
     5589  DO k = 1, klev
     5590    DO i = 1, klon
     5591      IF (wape(i)<0.) THEN
     5592        deltatw(i, k) = 0.
     5593        deltaqw(i, k) = 0.
     5594        dth(i, k) = 0.
     5595        d_deltatw2(i,k) = -deltatw0(i,k)
     5596        d_deltaqw2(i,k) = -deltaqw0(i,k)
     5597      END IF
     5598    END DO
     5599  END DO
     5600
     5601  DO i = 1, klon
     5602    IF (wape(i)<0.) THEN
     5603!!      sigmaw(i) = amax1(sigmad, sigd_con(i))
     5604      sigmaw_targ = max(sigmad, sigd_con(i))
     5605      d_sig_bnd2(i) = d_sig_bnd2(i) + sigmaw_targ - sigmaw(i)
     5606      d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i)
     5607      sigmaw(i) = sigmaw_targ
     5608    ENDIF  !!  (wape(i)<0.)
     5609  ENDDO
     5610  !
     5611  IF (iflag_wk_pop_dyn == 3) THEN
     5612    DO i = 1, klon
     5613      IF (wape(i)<0.) THEN
     5614        sigmaw_targ = max(sigmad, sigd_con(i))
     5615        d_asig_bnd2(i) = d_asig_bnd2(i) + sigmaw_targ - asigmaw(i)
     5616        d_asigmaw2(i) = d_asigmaw2(i) + sigmaw_targ - asigmaw(i)
     5617        asigmaw(i) = sigmaw_targ
     5618      ENDIF  !!  (wape(i)<0.)
     5619    ENDDO
     5620  ENDIF  !!  (iflag_wk_pop_dyn == 3)
     5621
     5622  DO i = 1, klon
     5623    IF (wape(i)<0.) THEN
     5624      wape(i) = 0.
     5625      cstar(i) = 0.
     5626      hw(i) = hwmin
     5627      fip(i) = 0.
     5628      gwake(i) = .FALSE.
     5629    ELSE
     5630      hw(i) = hw0(i)
     5631      cstar(i) = stark*sqrt(2.*wape(i))
     5632      gwake(i) = .TRUE.
     5633    END IF
     5634  END DO
     5635  !
     5636
     5637  ! Check qx and qw positivity
     5638  ! --------------------------
     5639  DO i = 1, klon
     5640    q0_min(i) = min((qb(i,1)-sigmaw(i)*deltaqw(i,1)),  &
     5641                    (qb(i,1)+(1.-sigmaw(i))*deltaqw(i,1)))
     5642  END DO
     5643  DO k = 2, klev
     5644    DO i = 1, klon
     5645      q1_min(i) = min((qb(i,k)-sigmaw(i)*deltaqw(i,k)), &
     5646                      (qb(i,k)+(1.-sigmaw(i))*deltaqw(i,k)))
     5647      IF (q1_min(i)<=q0_min(i)) THEN
     5648        q0_min(i) = q1_min(i)
     5649      END IF
     5650    END DO
     5651  END DO
     5652
     5653  DO i = 1, klon
     5654    ok_qx_qw(i) = q0_min(i) >= 0.
     5655    alpha(i) = 1.
     5656    alpha_tot(i) = 1.
     5657  END DO
     5658
     5659  IF (prt_level>=10) THEN
     5660    PRINT *, 'wake-4, sigmaw(igout), cstar(igout), wape(igout), ktop(igout) ', &
     5661                      sigmaw(igout), cstar(igout), wape(igout), ktop(igout)
     5662  ENDIF
     5663
     5664
     5665  ! C -----------------------------------------------------------------
     5666  ! Sub-time-stepping
     5667  ! -----------------
     5668
     5669!    wk_nsub and dtimesub definitions moved to begining of routine.
     5670!!  wk_nsub = 10
     5671!!  dtimesub = dtime/wk_nsub
     5672
     5673 
     5674  ! ------------------------------------------------------------------------
     5675  ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
     5676  ! ------------------------------------------------------------------------
     5677  !
     5678  DO isubstep = 1, wk_nsub
     5679  !
     5680  ! ------------------------------------------------------------------------
     5681  !
     5682    ! wk_adv is the logical flag enabling wake evolution in the time advance
     5683    ! loop
     5684    DO i = 1, klon
     5685      wk_adv(i) = ok_qx_qw(i) .AND. alpha(i) >= 1.
     5686    END DO
     5687IF (CPPKEY_IOPHYS_WK) THEN
     5688    IF (phys_sub) THEN
     5689     iwkadv=0.
     5690     IF (wk_adv(1)) iwkadv=1.
     5691     iokqxqw=0.
     5692     IF (ok_qx_qw(1)) iokqxqw=1.
     5693     CALL iophys_ecrit('iwkadv',1,'iwkadv','',iwkadv)
     5694     CALL iophys_ecrit('iokqxqw',1,'iokqxqw','',iokqxqw)
     5695     CALL iophys_ecrit('alpha',1,'alpha','',alpha(1))
     5696    ENDIF
     5697END IF
     5698    IF (prt_level>=10) THEN
     5699      PRINT *, 'wake-4.1, isubstep,wk_adv(igout),cstar(igout),wape(igout), ptop(igout) ', &
     5700                          isubstep,wk_adv(igout),cstar(igout),wape(igout), ptop(igout)
     5701     
     5702    ENDIF
     5703
     5704    ! cc nrlmd   Ajout d'un recalcul de wdens dans le cas d'un entrainement
     5705    ! negatif de ktop a kupper --------
     5706    ! cc           On calcule pour cela une densite wdens0 pour laquelle on
     5707    ! aurait un entrainement nul ---
     5708    !jyg<
     5709    ! Dans la configuration avec wdens prognostique, il s'agit d'un cas ou
     5710    ! les poches sont insuffisantes pour accueillir tout le flux de masse
     5711    ! des descentes unsaturees. Nous faisons alors l'hypothese que la
     5712    ! convection profonde cree directement de nouvelles poches, sans passer
     5713    ! par les thermiques. La nouvelle valeur de wdens est alors imposee.
     5714
     5715    DO i = 1, klon
     5716      ! c       print *,' isubstep,wk_adv(i),cstar(i),wape(i) ',
     5717      ! c     $           isubstep,wk_adv(i),cstar(i),wape(i)
     5718      IF (wk_adv(i) .AND. cstar(i)>0.01) THEN
     5719        IF ( iflag_wk_profile == 0 ) THEN
     5720           omg(i, kupper(i)+1)=-RG*amdwn(i, kupper(i)+1)/sigmaw(i) + &
     5721                               RG*amup(i, kupper(i)+1)/(1.-sigmaw(i))
     5722        ELSE
     5723           omg(i, kupper(i)+1)=0.
     5724        ENDIF
     5725        wdens0 = (sigmaw(i)/(4.*3.14))* &
     5726          ((1.-sigmaw(i))*omg(i,kupper(i)+1)/((ph(i,1)-pupper(i))*cstar(i)))**(2)
     5727        IF (prt_level >= 10) THEN
     5728             print*,'omg(i,kupper(i)+1),wdens0,wdens(i),cstar(i), ph(i,1)-pupper(i)', &
     5729                     omg(i,kupper(i)+1),wdens0,wdens(i),cstar(i), ph(i,1)-pupper(i)
     5730        ENDIF
     5731        IF (wdens(i)<=wdens0*1.1) THEN
     5732          IF (iflag_wk_pop_dyn >= 1) THEN
     5733             d_dens_bnd2(i) = d_dens_bnd2(i) + wdens0 - wdens(i)
     5734             d_wdens2(i) = d_wdens2(i) + wdens0 - wdens(i)
     5735          ENDIF
     5736          wdens(i) = wdens0
     5737        END IF
     5738      END IF
     5739    END DO
     5740
     5741    IF (iflag_wk_pop_dyn == 0 .AND. ok_bug_gfl) THEN
     5742!!--------------------------------------------------------
     5743!!Bug : computing gfl and rad_wk before changing sigmaw
     5744!!      This bug exists only for iflag_wk_pop_dyn=0. Otherwise, gfl and rad_wk
     5745!!      are computed within  wake_popdyn
     5746!!--------------------------------------------------------
     5747      DO i = 1, klon
     5748        IF (wk_adv(i)) THEN
     5749          gfl(i) = 2.*sqrt(3.14*wdens(i)*sigmaw(i))
     5750          rad_wk(i) = sqrt(sigmaw(i)/(3.14*wdens(i)))
     5751        END IF
     5752      END DO
     5753    ENDIF   ! (iflag_wk_pop_dyn == 0 .AND. ok_bug_gfl)
     5754!!--------------------------------------------------------
     5755
     5756    DO i = 1, klon
     5757      IF (wk_adv(i)) THEN
     5758        sigmaw_targ = min(sigmaw(i), sigmaw_max)
     5759        d_sig_bnd2(i) = d_sig_bnd2(i) + sigmaw_targ - sigmaw(i)
     5760        d_sigmaw2(i)  = d_sigmaw2(i)  + sigmaw_targ - sigmaw(i)
     5761        sigmaw(i) = sigmaw_targ
     5762      END IF
     5763    END DO
     5764
     5765    IF (iflag_wk_pop_dyn == 0 .AND. .NOT.ok_bug_gfl) THEN
     5766!!--------------------------------------------------------
     5767!!Fix : computing gfl and rad_wk after changing sigmaw
     5768!!--------------------------------------------------------
     5769      DO i = 1, klon
     5770        IF (wk_adv(i)) THEN
     5771          gfl(i) = 2.*sqrt(3.14*wdens(i)*sigmaw(i))
     5772          rad_wk(i) = sqrt(sigmaw(i)/(3.14*wdens(i)))
     5773        END IF
     5774      END DO
     5775    ENDIF   ! (iflag_wk_pop_dyn == 0 .AND. .NOT.ok_bug_gfl)
     5776!!--------------------------------------------------------
     5777
     5778    IF (iflag_wk_pop_dyn >= 1) THEN
     5779  !  The variable "death_rate" is significant only when iflag_wk_pop_dyn = 0.
     5780  !  Here, it has to be set to zero.
     5781      death_rate(:) = 0.
     5782    ENDIF
     5783 
     5784    IF (iflag_wk_pop_dyn >= 3) THEN
     5785      DO i = 1, klon
     5786        IF (wk_adv(i)) THEN
     5787         sigmaw_targ = min(asigmaw(i), sigmaw_max)
     5788         d_asig_bnd2(i) = d_asig_bnd2(i) + sigmaw_targ - asigmaw(i)
     5789         d_asigmaw2(i)  = d_asigmaw2(i)  + sigmaw_targ - asigmaw(i)
     5790         asigmaw(i) = sigmaw_targ
     5791        ENDIF
     5792      ENDDO
     5793    ENDIF
     5794 
     5795!!--------------------------------------------------------
     5796!!--------------------------------------------------------
     5797    IF (iflag_wk_pop_dyn == 1) THEN
     5798  !
     5799     CALL wake_popdyn_1 (klon, klev, dtime, cstar, tau_wk_inv, wgen, wdens, awdens, sigmaw, &
     5800                  wdensmin, &
     5801                  dtimesub, gfl, rad_wk, f_shear, drdt_pos, &
     5802                  d_awdens, d_wdens, d_sigmaw, &
     5803                  iflag_wk_act, wk_adv, cin, wape, &
     5804                  drdt, &
     5805                  d_dens_gen, d_dens_death, d_dens_col, d_dens_bnd, &
     5806                  d_sig_gen, d_sig_death, d_sig_col, d_sig_spread, d_sig_bnd, &
     5807                  d_wdens_targ, d_sigmaw_targ)
     5808                     
     5809   
     5810!!--------------------------------------------------------
     5811    ELSEIF (iflag_wk_pop_dyn == 2) THEN
     5812  !
     5813     CALL wake_popdyn_2 ( klon, klev, wk_adv, dtimesub, wgen, &
     5814                             wdensmin, &
     5815                             sigmaw, wdens, awdens, &   !! state variables
     5816                             gfl, cstar, cin, wape, rad_wk, &
     5817                             d_sigmaw, d_wdens, d_awdens, &  !! tendencies                             
     5818                             cont_fact, &
     5819                             d_sig_gen, d_sig_death, d_sig_col, d_sig_spread, d_sig_bnd, &
     5820                             d_dens_gen, d_dens_death, d_dens_col, d_dens_bnd, &
     5821                             d_adens_death, d_adens_icol, d_adens_acol, d_adens_bnd )
     5822sigmaw=sigmaw-d_sigmaw
     5823wdens=wdens-d_wdens
     5824awdens=awdens-d_awdens
     5825
     5826!!--------------------------------------------------------
     5827    ELSEIF (iflag_wk_pop_dyn == 3) THEN
     5828IF (CPPKEY_IOPHYS_WK) THEN
     5829    IF (phys_sub) THEN
     5830     CALL iophys_ecrit('ptop',1,'ptop','Pa',ptop)
     5831     CALL iophys_ecrit('wape',1,'wape','J/kg',wape)
     5832     CALL iophys_ecrit('wgen',1,'wgen','1/(s.m2)',wgen)
     5833     CALL iophys_ecrit('sigmaw',1,'sigmaw','',sigmaw)
     5834     CALL iophys_ecrit('asigmaw',1,'asigmaw','',asigmaw)
     5835     CALL iophys_ecrit('wdens',1,'wdens','1/m2',wdens)
     5836     CALL iophys_ecrit('awdens',1,'awdens','1/m2',awdens)
     5837    ENDIF
     5838END IF
     5839  !
     5840     CALL wake_popdyn_3 ( klon, klev, phys_sub, wk_adv, dtimesub, wgen, &
     5841                             wdensmin, &
     5842                             sigmaw, asigmaw, wdens, awdens, &   !! state variables
     5843                             gfl, agfl, cstar, cin, wape, &
     5844                             rad_wk, arad_wk, irad_wk, &
     5845                             d_sigmaw, d_asigmaw, d_wdens, d_awdens, &  !! tendencies                             
     5846                             d_sig_gen, d_sig_death, d_sig_col, d_sig_spread, d_sig_bnd, &
     5847                             d_asig_death, d_asig_aicol, d_asig_iicol, d_asig_spread, d_asig_bnd, &
     5848                             d_dens_gen, d_dens_death, d_dens_col, d_dens_bnd, &
     5849                             d_adens_death, d_adens_icol, d_adens_acol, d_adens_bnd )
     5850IF (CPPKEY_IOPHYS_WK) THEN
     5851    IF (phys_sub) THEN
     5852     CALL iophys_ecrit('rad_wk',1,'rad_wk','m',rad_wk)
     5853     CALL iophys_ecrit('arad_wk',1,'arad_wk','m',arad_wk)
     5854     CALL iophys_ecrit('irad_wk',1,'irad_wk','m',irad_wk)
     5855    ENDIF
     5856END IF
     5857sigmaw=sigmaw-d_sigmaw
     5858asigmaw=asigmaw-d_asigmaw
     5859wdens=wdens-d_wdens
     5860awdens=awdens-d_awdens
     5861   
     5862!!--------------------------------------------------------
     5863    ELSEIF (iflag_wk_pop_dyn == 0) THEN
     5864   
     5865    ! cc nrlmd
     5866
     5867      DO i = 1, klon
     5868        IF (wk_adv(i)) THEN
     5869
     5870          ! cc nrlmd          Introduction du taux de mortalite des poches et
     5871          ! test sur sigmaw_max=0.4
     5872          ! cc         d_sigmaw(i) = gfl(i)*Cstar(i)*dtimesub
     5873          IF (sigmaw(i)>=sigmaw_max) THEN
     5874            death_rate(i) = gfl(i)*cstar(i)/sigmaw(i)
     5875          ELSE
     5876            death_rate(i) = 0.
     5877          END IF
     5878   
     5879          d_sigmaw(i) = gfl(i)*cstar(i)*dtimesub - death_rate(i)*sigmaw(i)* &
     5880            dtimesub
     5881          ! $              - nat_rate(i)*sigmaw(i)*dtimesub
     5882          ! c        print*, 'd_sigmaw(i),sigmaw(i),gfl(i),Cstar(i),wape(i),
     5883          ! c     $  death_rate(i),ktop(i),kupper(i)',
     5884          ! c     $              d_sigmaw(i),sigmaw(i),gfl(i),Cstar(i),wape(i),
     5885          ! c     $  death_rate(i),ktop(i),kupper(i)
     5886   
     5887          ! sigmaw(i) =sigmaw(i) + gfl(i)*Cstar(i)*dtimesub
     5888          ! sigmaw(i) =min(sigmaw(i),0.99)     !!!!!!!!
     5889          ! wdens = wdens0/(10.*sigmaw)
     5890          ! sigmaw =max(sigmaw,sigd_con)
     5891          ! sigmaw =max(sigmaw,sigmad)
     5892        END IF
     5893      END DO
     5894
     5895    ENDIF   !  (iflag_wk_pop_dyn == 1) ... ELSEIF (iflag_wk_pop_dyn == 0)
     5896!!--------------------------------------------------------
     5897!!--------------------------------------------------------
     5898
     5899IF (CPPKEY_IOPHYS_WK) THEN
     5900    IF (phys_sub) THEN
     5901     CALL iophys_ecrit('wdensa',1,'wdensa','m',wdens)
     5902     CALL iophys_ecrit('awdensa',1,'awdensa','m',awdens)
     5903     CALL iophys_ecrit('sigmawa',1,'sigmawa','m',sigmaw)
     5904     CALL iophys_ecrit('asigmawa',1,'asigmawa','m',asigmaw)
     5905    ENDIF
     5906END IF
     5907    ! calcul de la difference de vitesse verticale poche - zone non perturbee
     5908    ! IM 060208 differences par rapport au code initial; init. a 0 dp_deltomg
     5909    ! IM 060208 et omg sur les niveaux de 1 a klev+1, alors que avant l'on definit
     5910    ! IM 060208 au niveau k=1...
     5911    !JYG 161013 Correction : maintenant omg est dimensionne a klev.
     5912    DO k = 1, klev
     5913      DO i = 1, klon
     5914        IF (wk_adv(i)) THEN !!! nrlmd
     5915          dp_deltomg(i, k) = 0.
     5916        END IF
     5917      END DO
     5918    END DO
     5919    DO k = 1, klev
     5920      DO i = 1, klon
     5921        IF (wk_adv(i)) THEN !!! nrlmd
     5922          omg(i, k) = 0.
     5923        END IF
     5924      END DO
     5925    END DO
     5926
     5927    DO i = 1, klon
     5928      IF (wk_adv(i)) THEN
     5929        z(i) = 0.
     5930        omg(i, 1) = 0.
     5931        dp_deltomg(i, 1) = -(gfl(i)*cstar(i))/(sigmaw(i)*(1-sigmaw(i)))
     5932      END IF
     5933    END DO
     5934
     5935    DO k = 2, klev
     5936      DO i = 1, klon
     5937        IF (wk_adv(i) .AND. k<=ktop(i)) THEN
     5938          dz(i) = -(ph(i,k)-ph(i,k-1))/(rho(i,k-1)*RG)
     5939          z(i) = z(i) + dz(i)
     5940          dp_deltomg(i, k) = dp_deltomg(i, 1)
     5941          omg(i, k) = dp_deltomg(i, 1)*z(i)
     5942        END IF
     5943      END DO
     5944    END DO
     5945
     5946    DO i = 1, klon
     5947      IF (wk_adv(i)) THEN
     5948        dztop(i) = -(ptop(i)-ph(i,ktop(i)))/(rho(i,ktop(i))*RG)
     5949        ztop(i) = z(i) + dztop(i)
     5950        omgtop(i) = dp_deltomg(i, 1)*ztop(i)
     5951      END IF
     5952    END DO
     5953
     5954    IF (prt_level>=10) THEN
     5955      PRINT *, 'wake-4.2, omg(igout,k) ', (k,omg(igout,k), k=1,klev)
     5956      PRINT *, 'wake-4.2, omgtop(igout), ptop(igout), ktop(igout) ', &
     5957                          omgtop(igout), ptop(igout), ktop(igout)
     5958    ENDIF
     5959
     5960    ! -----------------
     5961    ! From m/s to Pa/s
     5962    ! -----------------
     5963
     5964    DO i = 1, klon
     5965      IF (wk_adv(i)) THEN
     5966        omgtop(i) = -rho(i, ktop(i))*RG*omgtop(i)
     5967!! LJYF        dp_deltomg(i, 1) = omgtop(i)/(ptop(i)-ph(i,1))
     5968        dp_deltomg(i, 1) = omgtop(i)/min(ptop(i)-ph(i,1),-smallestreal)
     5969      END IF
     5970    END DO
     5971
     5972    DO k = 1, klev
     5973      DO i = 1, klon
     5974        IF (wk_adv(i) .AND. k<=ktop(i)) THEN
     5975          omg(i, k) = -rho(i, k)*RG*omg(i, k)
     5976          dp_deltomg(i, k) = dp_deltomg(i, 1)
     5977        END IF
     5978      END DO
     5979    END DO
     5980
     5981    ! raccordement lineaire de omg de ptop a pupper
     5982
     5983    DO i = 1, klon
     5984      IF (wk_adv(i) .AND. kupper(i)>ktop(i)) THEN
     5985        IF ( iflag_wk_profile == 0 ) THEN
     5986           omg(i, kupper(i)+1) =-RG*amdwn(i, kupper(i)+1)/sigmaw(i) + &
     5987          RG*amup(i, kupper(i)+1)/(1.-sigmaw(i))
     5988        ELSE
     5989           omg(i, kupper(i)+1) = 0.
     5990        ENDIF
     5991        dp_deltomg(i, kupper(i)) = (omgtop(i)-omg(i,kupper(i)+1))/ &
     5992          (ptop(i)-pupper(i))
     5993      END IF
     5994    END DO
     5995
     5996    ! c      DO i=1,klon
     5997    ! c        print*,'Pente entre 0 et kupper (reference)'
     5998    ! c     $           ,omg(i,kupper(i)+1)/(pupper(i)-ph(i,1))
     5999    ! c        print*,'Pente entre ktop et kupper'
     6000    ! c     $   ,(omg(i,kupper(i)+1)-omgtop(i))/(pupper(i)-ptop(i))
     6001    ! c      ENDDO
     6002    ! c
     6003    DO k = 1, klev
     6004      DO i = 1, klon
     6005        IF (wk_adv(i) .AND. k>ktop(i) .AND. k<=kupper(i)) THEN
     6006          dp_deltomg(i, k) = dp_deltomg(i, kupper(i))
     6007          omg(i, k) = omgtop(i) + (ph(i,k)-ptop(i))*dp_deltomg(i, kupper(i))
     6008        END IF
     6009      END DO
     6010    END DO
     6011!!    print *,'omg(igout,k) ', (k,omg(igout,k),k=1,klev)
     6012    ! cc nrlmd
     6013    ! c      DO i=1,klon
     6014    ! c      print*,'deltaw_ktop,deltaw_conv',omgtop(i),omg(i,kupper(i)+1)
     6015    ! c      END DO
     6016    ! cc
     6017
     6018
     6019    ! --    Compute wake average vertical velocity omgbw
     6020
     6021
     6022    DO k = 1, klev
     6023      DO i = 1, klon
     6024        IF (wk_adv(i)) THEN
     6025          omgbw(i, k) = omgb(i, k) + (1.-sigmaw(i))*omg(i, k)
     6026        END IF
     6027      END DO
     6028    END DO
     6029    ! --    and its vertical gradient dp_omgbw
     6030
     6031    DO k = 1, klev-1
     6032      DO i = 1, klon
     6033        IF (wk_adv(i)) THEN
     6034          dp_omgbw(i, k) = (omgbw(i,k+1)-omgbw(i,k))/(ph(i,k+1)-ph(i,k))
     6035        END IF
     6036      END DO
     6037    END DO
     6038    DO i = 1, klon
     6039      IF (wk_adv(i)) THEN
     6040          dp_omgbw(i, klev) = 0.
     6041      END IF
     6042    END DO
     6043
     6044    ! --    Upstream coefficients for omgb velocity
     6045    ! --    (alpha_up(k) is the coefficient of the value at level k)
     6046    ! --    (1-alpha_up(k) is the coefficient of the value at level k-1)
     6047    DO k = 1, klev
     6048      DO i = 1, klon
     6049        IF (wk_adv(i)) THEN
     6050          alpha_up(i, k) = 0.
     6051          IF (omgb(i,k)>0.) alpha_up(i, k) = 1.
     6052        END IF
     6053      END DO
     6054    END DO
     6055
     6056    ! Matrix expressing [The,deltatw] from  [Th1,Th2]
     6057
     6058    DO i = 1, klon
     6059      IF (wk_adv(i)) THEN
     6060        rre1(i) = 1. - sigmaw(i)
     6061        rre2(i) = sigmaw(i)
     6062      END IF
     6063    END DO
     6064    rrd1 = -1.
     6065    rrd2 = 1.
     6066
     6067    ! --    Get [Th1,Th2], dth and [q1,q2]
     6068
     6069    DO k = 1, klev
     6070      DO i = 1, klon
     6071        IF (wk_adv(i) .AND. k<=kupper(i)+1) THEN
     6072          dth(i, k) = deltatw(i, k)/ppi(i, k)
     6073          th1(i, k) = thb(i, k) - sigmaw(i)*dth(i, k) ! undisturbed area
     6074          th2(i, k) = thb(i, k) + (1.-sigmaw(i))*dth(i, k) ! wake
     6075          q1(i, k) = qb(i, k) - sigmaw(i)*deltaqw(i, k) ! undisturbed area
     6076          q2(i, k) = qb(i, k) + (1.-sigmaw(i))*deltaqw(i, k) ! wake
     6077        END IF
     6078      END DO
     6079    END DO
     6080
     6081    DO i = 1, klon
     6082      IF (wk_adv(i)) THEN !!! nrlmd
     6083        d_th1(i, 1) = 0.
     6084        d_th2(i, 1) = 0.
     6085        d_dth(i, 1) = 0.
     6086        d_q1(i, 1) = 0.
     6087        d_q2(i, 1) = 0.
     6088        d_dq(i, 1) = 0.
     6089      END IF
     6090    END DO
     6091
     6092    DO k = 2, klev
     6093      DO i = 1, klon
     6094        IF (wk_adv(i) .AND. k<=kupper(i)+1) THEN
     6095          d_th1(i, k) = th1(i, k-1) - th1(i, k)
     6096          d_th2(i, k) = th2(i, k-1) - th2(i, k)
     6097          d_dth(i, k) = dth(i, k-1) - dth(i, k)
     6098          d_q1(i, k) = q1(i, k-1) - q1(i, k)
     6099          d_q2(i, k) = q2(i, k-1) - q2(i, k)
     6100          d_dq(i, k) = deltaqw(i, k-1) - deltaqw(i, k)
     6101        END IF
     6102      END DO
     6103    END DO
     6104
     6105    DO i = 1, klon
     6106      IF (wk_adv(i)) THEN
     6107        omgbdth(i, 1) = 0.
     6108        omgbdq(i, 1) = 0.
     6109      END IF
     6110    END DO
     6111
     6112    DO k = 2, klev
     6113      DO i = 1, klon
     6114        IF (wk_adv(i) .AND. k<=kupper(i)+1) THEN !   loop on interfaces
     6115          omgbdth(i, k) = omgb(i, k)*(dth(i,k-1)-dth(i,k))
     6116          omgbdq(i, k) = omgb(i, k)*(deltaqw(i,k-1)-deltaqw(i,k))
     6117        END IF
     6118      END DO
     6119    END DO
     6120
     6121!!    IF (prt_level>=10) THEN
     6122    IF (prt_level>=10 .and. wk_adv(igout)) THEN
     6123      PRINT *, 'wake-4.3, th1(igout,k) ', (k,th1(igout,k), k=1,kupper(igout))
     6124      PRINT *, 'wake-4.3, th2(igout,k) ', (k,th2(igout,k), k=1,kupper(igout))
     6125      PRINT *, 'wake-4.3, dth(igout,k) ', (k,dth(igout,k), k=1,kupper(igout))
     6126      PRINT *, 'wake-4.3, omgbdth(igout,k) ', (k,omgbdth(igout,k), k=1,kupper(igout))
     6127    ENDIF
     6128
     6129
     6130    ! -----------------------------------------------------------------
     6131          ! Compute redistribution (advective) term
     6132
     6133!     rate of change of sigmaw due to spreading
     6134          dsigspread(:) = gfl(:)*cstar(:)
     6135
     6136          CALL wake_dadv(klon, klev, dtimesub, ph, ppi, wk_adv, kupper, &
     6137                    omg, dp_deltomg, sigmaw, dsigspread, &
     6138                    th2, th1, q2, q1, &
     6139                    d_deltat_dadv, d_deltaq_dadv, d_tb_dadv, d_qb_dadv)
     6140
     6141    ! For the difference fields: convert to change per second in order to combine with the
     6142    ! other terms (d_deltat_ls, d_deltat_cv, d_deltat_gw)
     6143    d_deltat_dadv(:,:) = d_deltat_dadv(:,:)/dtimesub
     6144    d_deltaq_dadv(:,:) = d_deltaq_dadv(:,:)/dtimesub
     6145!
     6146    !   For the mean fields: tb and qb the computation of the tendencies due to wakes is
     6147    !   already complete.
     6148    d_tb(:,:) = d_tb_dadv(:,:)
     6149    d_qb(:,:) = d_qb_dadv(:,:)
     6150
     6151    ! -----------------------------------------------------------------
     6152    DO k = 1, klev-1
     6153      DO i = 1, klon
     6154        IF (wk_adv(i) .AND. k<=kupper(i)-1) THEN
     6155          ! -----------------------------------------------------------------
     6156          d_deltat_lsadv(i, k) = 1./(ph(i,k)-ph(i,k+1))* &
     6157            (-(1.-alpha_up(i,k))*omgbdth(i,k)- &
     6158             alpha_up(i,k+1)*omgbdth(i,k+1))*ppi(i, k)
     6159
     6160          d_deltaq_lsadv(i, k) = 1./(ph(i,k)-ph(i,k+1))* &
     6161            (-(1.-alpha_up(i,k))*omgbdq(i,k)- &
     6162             alpha_up(i,k+1)*omgbdq(i,k+1))
     6163
     6164        END IF
     6165        ! cc
     6166      END DO
     6167    END DO
     6168    ! ------------------------------------------------------------------
     6169
     6170!!    IF (prt_level>=10) THEN
     6171    IF (prt_level>=10 .and. wk_adv(igout)) THEN
     6172      PRINT *, 'wake-4.3, d_deltat_dadv(igout,k) ', (k,d_deltat_dadv(igout,k), k=1,klev)
     6173      PRINT *, 'wake-4.3, d_deltat_lsadv(igout,k) ', (k,d_deltat_lsadv(igout,k), k=1,klev)
     6174      PRINT *, 'wake-4.3, d_deltaq_dadv(igout,k) ', (k,d_deltaq_dadv(igout,k), k=1,klev)
     6175      PRINT *, 'wake-4.3, d_deltaq_lsadv(igout,k) ', (k,d_deltaq_lsadv(igout,k), k=1,klev)
     6176    ENDIF
     6177
     6178    ! Increment state variables
     6179!jyg<
     6180    IF (iflag_wk_pop_dyn >= 1) THEN
     6181      DO k = 1, klev
     6182        DO i = 1, klon
     6183          IF (wk_adv(i) .AND. k<=kupper(i)) THEN
     6184            detr(i,k) = - d_sig_death(i) - d_sig_col(i)     
     6185            entr_p(i,k) = d_sig_gen(i)
     6186          ENDIF
     6187        ENDDO
     6188      ENDDO
     6189      ELSE  ! (iflag_wk_pop_dyn >= 1)
     6190      DO k = 1, klev
     6191        DO i = 1, klon
     6192          IF (wk_adv(i) .AND. k<=kupper(i)) THEN
     6193            detr(i, k) = 0.
     6194   
     6195            entr_p(i, k) = 0.
     6196          ENDIF
     6197        ENDDO
     6198      ENDDO
     6199    ENDIF  ! (iflag_wk_pop_dyn >= 1)
     6200
     6201   
     6202
     6203    DO k = 1, klev
     6204      DO i = 1, klon
     6205        ! cc nrlmd       IF( wk_adv(i) .AND. k .LE. kupper(i)-1) THEN
     6206        IF (wk_adv(i) .AND. k<=kupper(i)) THEN
     6207          ! cc
     6208
     6209
     6210
     6211          ! Coefficient de repartition
     6212
     6213          crep(i, k) = crep_sol*(ph(i,kupper(i))-ph(i,k))/ &
     6214            (ph(i,kupper(i))-ph(i,1))
     6215          crep(i, k) = crep(i, k) + crep_upper*(ph(i,1)-ph(i,k))/ &
     6216            (ph(i,1)-ph(i,kupper(i)))
     6217
     6218
     6219          ! Reintroduce compensating subsidence term.
     6220
     6221          ! dtKE(k)=(dtdwn(k)*Crep(k))/sigmaw
     6222          ! dtKE(k)=dtKE(k)-(dtdwn(k)*(1-Crep(k))+dta(k))
     6223          ! .                   /(1-sigmaw)
     6224          ! dqKE(k)=(dqdwn(k)*Crep(k))/sigmaw
     6225          ! dqKE(k)=dqKE(k)-(dqdwn(k)*(1-Crep(k))+dqa(k))
     6226          ! .                   /(1-sigmaw)
     6227
     6228          ! dtKE(k)=(dtdwn(k)*Crep(k)+(1-Crep(k))*dta(k))/sigmaw
     6229          ! dtKE(k)=dtKE(k)-(dtdwn(k)*(1-Crep(k))+dta(k)*Crep(k))
     6230          ! .                   /(1-sigmaw)
     6231          ! dqKE(k)=(dqdwn(k)*Crep(k)+(1-Crep(k))*dqa(k))/sigmaw
     6232          ! dqKE(k)=dqKE(k)-(dqdwn(k)*(1-Crep(k))+dqa(k)*Crep(k))
     6233          ! .                   /(1-sigmaw)
     6234
     6235          !! Differential heating (d_deltat_dcv) and moistening (d_deltaq_dcv) by deep convection
     6236          d_deltat_dcv(i, k) = (dtdwn(i,k)/sigmaw(i)-dta(i,k)/(1.-sigmaw(i)))
     6237!!          dtke(i, k) = (dtdwn(i,k)/sigmaw(i)-dta(i,k)/(1.-sigmaw(i)))        ! supprime
     6238          d_deltaq_dcv(i, k) = (dqdwn(i,k)/sigmaw(i)-dqa(i,k)/(1.-sigmaw(i)))
     6239!!          dqke(i, k) = (dqdwn(i,k)/sigmaw(i)-dqa(i,k)/(1.-sigmaw(i)))        ! supprime
     6240          ! print*,'dtKE= ',dtKE(i,k),' dqKE= ',dqKE(i,k)
     6241
     6242!
     6243
     6244          ! cc nrlmd          Prise en compte du taux de mortalite
     6245          ! cc               Definitions de entr, detr
     6246!jyg<
     6247!!            detr(i, k) = 0.
     6248!!   
     6249!!            entr(i, k) = detr(i, k) + gfl(i)*cstar(i) + &
     6250!!              sigmaw(i)*(1.-sigmaw(i))*dp_deltomg(i, k)
     6251!!
     6252            entr(i, k) = entr_p(i,k) + gfl(i)*cstar(i) + &
     6253                         sigmaw(i)*(1.-sigmaw(i))*dp_deltomg(i, k)   
     6254            tgen(i,k) = entr_p(i,k)/sigmaw(i)
     6255!>jyg
     6256            wkspread(i, k) = (entr(i,k)-detr(i,k))/sigmaw(i)
     6257
     6258          ! cc        wkspread(i,k) =
     6259          ! (1.-sigmaw(i))*dp_deltomg(i,k)+gfl(i)*Cstar(i)/
     6260          ! cc     $  sigmaw(i)
     6261
     6262
     6263          ! ajout d'un effet onde de gravite -Tgw(k)*deltatw(k) 03/02/06 YU
     6264          ! Jingmei
     6265
     6266          ! write(lunout,*)'wake.F ',i,k, dtimesub,d_deltat_gw(i,k),
     6267          ! &  Tgw(i,k),deltatw(i,k)
     6268          d_deltat_gw(i, k) = d_deltat_gw(i, k) - tgw(i, k)*deltatw(i, k)* &
     6269            dtimesub
     6270          ! write(lunout,*)'wake.F ',i,k, dtimesub,d_deltatw(i,k)
     6271
     6272          ! Sans GW
     6273
     6274          ! deltatw(k)=deltatw(k)+dtimesub*(ff+dtKE(k)-wkspread(k)*deltatw(k))
     6275
     6276          ! GW formule 1
     6277
     6278          ! deltatw(k) = deltatw(k)+dtimesub*
     6279          ! $         (ff+dtKE(k) - wkspread(k)*deltatw(k)-Tgw(k)*deltatw(k))
     6280
     6281          ! GW formule 2
     6282
     6283          !! Entrainment due to spread is supposed to be included in the differential advection
     6284          !! term (d_deltat_dadv); hence only the entrainment due to population dynamics (entr_p)
     6285          !! appears in the expression of d_deltatw.
     6286          IF (dtimesub*(tgw(i,k)+tgen(i,k))<1.E-10) THEN
     6287!!!            d_deltatw(i, k) = dtimesub*(ff(i)+dtke(i,k) - &        ! nouvelle notation
     6288            d_deltatw(i, k) = dtimesub*(d_deltat_dadv(i,k)+d_deltat_lsadv(i,k)+d_deltat_dcv(i,k) - &
     6289               (death_rate(i)*sigmaw(i)+detr(i,k))*deltatw(i,k)/(1.-sigmaw(i)) - & ! cc
     6290               (tgw(i,k)+tgen(i,k))*deltatw(i,k) )
     6291          ELSE
     6292            d_deltatw(i, k) = 1/(tgw(i,k)+tgen(i,k))*(1-exp(-dtimesub*(tgw(i,k)+tgen(i,k))))* &
     6293!!!               (ff(i)+dtke(i,k) - &                                ! nouvelle notation
     6294               (d_deltat_dadv(i,k)+d_deltat_lsadv(i,k)+d_deltat_dcv(i,k) - &
     6295                (death_rate(i)*sigmaw(i)+detr(i,k))*deltatw(i,k)/(1.-sigmaw(i)) - &
     6296                (tgw(i,k)+tgen(i,k))*deltatw(i,k) )
     6297          END IF
     6298
     6299          dth(i, k) = deltatw(i, k)/ppi(i, k)
     6300
     6301          !! Entrainment due to spread is supposed to be included in the differential advection
     6302          !! term (d_deltaq_dadv); hence only the entrainment due to population dynamics (entr_p)
     6303          !! appears in the expression of d_deltaqw.
     6304          IF (dtimesub*tgen(i,k)<1.E-10) THEN
     6305            d_deltaqw(i, k) = dtimesub*(d_deltaq_dadv(i,k)+d_deltaq_lsadv(i,k)+d_deltaq_dcv(i,k) - &
     6306               (death_rate(i)*sigmaw(i)+detr(i,k))*deltaqw(i,k)/(1.-sigmaw(i)) - &
     6307               tgen(i,k)*deltaqw(i,k))
     6308          ELSE
     6309            d_deltaqw(i, k) = 1/tgen(i,k)*(1-exp(-dtimesub*tgen(i,k))) * &
     6310               (d_deltaq_dadv(i,k)+d_deltaq_lsadv(i,k)+d_deltaq_dcv(i,k) - &
     6311               (death_rate(i)*sigmaw(i)+detr(i,k))*deltaqw(i,k)/(1.-sigmaw(i)) - &
     6312               tgen(i,k)*deltaqw(i,k))
     6313          END IF
     6314          ! cc
     6315
     6316          ! cc nrlmd
     6317          ! cc       d_deltatw2(i,k)=d_deltatw2(i,k)+d_deltatw(i,k)
     6318          ! cc       d_deltaqw2(i,k)=d_deltaqw2(i,k)+d_deltaqw(i,k)
     6319          ! cc
     6320        END IF
     6321      END DO
     6322    END DO
     6323
     6324!!    IF (prt_level>=10) THEN
     6325    IF (prt_level>=10 .and. wk_adv(igout)) THEN
     6326      PRINT *, 'wake-4.4, isubstep= ', isubstep,' deltatw(igout,k) ', (k,deltatw(igout,k), k=1,klev)
     6327      PRINT *, 'wake-4.4, isubstep= ', isubstep,' d_deltat_dcv(igout,k) ', (k,d_deltat_dcv(igout,k), k=1,klev)
     6328      PRINT *, 'wake-4.4, isubstep= ', isubstep,' d_deltat_dadv(igout,k) ', (k,d_deltat_dadv(igout,k), k=1,klev)
     6329      PRINT *, 'wake-4.4, isubstep= ', isubstep,' d_deltat_lsadv(igout,k) ', (k,d_deltat_lsadv(igout,k), k=1,klev)
     6330      PRINT *, 'wake-4.4, isubstep= ', isubstep,' tgen(igout,k)*deltatw(igout,k) ', (k,tgen(igout,k)*deltatw(igout,k), k=1,klev)
     6331      PRINT *, 'wake-4.4, isubstep= ', isubstep,' tgw(igout,k)*deltatw(igout,k) ', (k,tgw(igout,k)*deltatw(igout,k), k=1,klev)
     6332      PRINT *, 'wake-4.4, isubstep= ', isubstep,' death_rate(igout) ', death_rate(igout)
     6333      PRINT *, 'wake-4.4, isubstep= ', isubstep,' detr(igout,k) ',  (k,detr(igout,k), k=1,klev)
     6334      PRINT *, 'wake-4.4, isubstep= ', isubstep,' d_deltatw(igout,k) ', (k,d_deltatw(igout,k), k=1,klev)
     6335      PRINT *, 'wake-4.4, isubstep= ', isubstep,' d_deltaqw(igout,k) ', (k,d_deltaqw(igout,k), k=1,klev)
     6336      PRINT *, 'wake-4.4, isubstep= ', isubstep,' d_tb(igout,k) ', (k,d_tb(igout,k), k=1,klev)
     6337      PRINT *, 'wake-4.4, isubstep= ', isubstep,' d_qb(igout,k) ', (k,d_qb(igout,k), k=1,klev)
     6338    ENDIF
     6339
     6340!
     6341IF (CPPKEY_IOPHYS_WK) THEN
     6342    IF (phys_sub) THEN
     6343     DO k = 1,klev
     6344      d_deltat_entrp(:,k) = - entr_p(:,k)*deltatw(:,k)/sigmaw(:)
     6345     ENDDO
     6346     CALL iophys_ecrit('d_deltatw',klev,'d_deltatw','K/s',d_deltatw(:,1:klev))
     6347     CALL iophys_ecrit('d_deltat_dadv',klev,'d_deltat_dadv','K/s',d_deltat_dadv(:,1:klev))
     6348     CALL iophys_ecrit('d_deltat_lsadv',klev,'d_deltat_lsadv','K/s',d_deltat_lsadv(:,1:klev))
     6349     CALL iophys_ecrit('d_deltat_dcv',klev,'d_deltat_dcv','K/s',d_deltat_dcv(:,1:klev))
     6350     CALL iophys_ecrit('d_deltat_entrp',klev,'d_deltat_entrp','K/s',d_deltat_entrp(:,1:klev))
     6351
     6352    ENDIF
     6353END IF
     6354
     6355    ! Scale tendencies so that water vapour remains positive in w and x.
     6356
     6357    CALL wake_vec_modulation(klon, klev, wk_adv, epsilon_loc, qb, d_qb, deltaqw, &
     6358      d_deltaqw, sigmaw, d_sigmaw, alpha)
     6359    !
     6360    ! Alpha_tot = Product of all the alpha's
     6361    DO i = 1, klon
     6362      IF (wk_adv(i)) THEN
     6363        alpha_tot(i) = alpha_tot(i)*alpha(i)   
     6364      END IF
     6365    END DO
     6366
     6367    ! cc nrlmd
     6368    ! c      print*,'alpha'
     6369    ! c      do i=1,klon
     6370    ! c         print*,alpha(i)
     6371    ! c      end do
     6372    ! cc
     6373    DO k = 1, klev
     6374      DO i = 1, klon
     6375        IF (wk_adv(i) .AND. k<=kupper(i)) THEN
     6376          d_tb(i, k) = alpha(i)*d_tb(i, k)
     6377          d_qb(i, k) = alpha(i)*d_qb(i, k)
     6378          d_deltatw(i, k) = alpha(i)*d_deltatw(i, k)
     6379          d_deltaqw(i, k) = alpha(i)*d_deltaqw(i, k)
     6380          d_deltat_gw(i, k) = alpha(i)*d_deltat_gw(i, k)
     6381        END IF
     6382      END DO
     6383    END DO
     6384    DO i = 1, klon
     6385      IF (wk_adv(i)) THEN
     6386        d_sigmaw(i) = alpha(i)*d_sigmaw(i)
     6387      END IF
     6388    END DO
     6389
     6390    ! Update large scale variables and wake variables
     6391    ! IM 060208 manque DO i + remplace DO k=1,kupper(i)
     6392    ! IM 060208     DO k = 1,kupper(i)
     6393    DO k = 1, klev
     6394      DO i = 1, klon
     6395        IF (wk_adv(i) .AND. k<=kupper(i)) THEN
     6396          dtls(i, k) = dtls(i, k) + d_tb(i, k)
     6397          dqls(i, k) = dqls(i, k) + d_qb(i, k)
     6398          ! cc nrlmd
     6399          d_deltatw2(i, k) = d_deltatw2(i, k) + d_deltatw(i, k)
     6400          d_deltaqw2(i, k) = d_deltaqw2(i, k) + d_deltaqw(i, k)
     6401          ! cc
     6402        END IF
     6403      END DO
     6404    END DO
     6405    DO k = 1, klev
     6406      DO i = 1, klon
     6407        IF (wk_adv(i) .AND. k<=kupper(i)) THEN
     6408          tb(i, k) = tb0(i, k) + dtls(i, k)
     6409          qb(i, k) = qb0(i, k) + dqls(i, k)
     6410          thb(i, k) = tb(i, k)/ppi(i, k)
     6411          deltatw(i, k) = deltatw(i, k) + d_deltatw(i, k)
     6412          deltaqw(i, k) = deltaqw(i, k) + d_deltaqw(i, k)
     6413          dth(i, k) = deltatw(i, k)/ppi(i, k)
     6414          ! c      print*,'k,qx,qw',k,qb(i,k)-sigmaw(i)*deltaqw(i,k)
     6415          ! c     $        ,qb(i,k)+(1-sigmaw(i))*deltaqw(i,k)
     6416        END IF
     6417      END DO
     6418    END DO
     6419!
     6420    DO i = 1, klon
     6421      IF (wk_adv(i)) THEN
     6422        sigmaw(i) = sigmaw(i) + d_sigmaw(i)
     6423        d_sigmaw2(i) = d_sigmaw2(i) + d_sigmaw(i)
     6424      END IF
     6425    END DO
     6426!jyg<
     6427    IF (iflag_wk_pop_dyn >= 1) THEN
     6428!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! sigmaw !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     6429!  Cumulatives
     6430      DO i = 1, klon
     6431        IF (wk_adv(i)) THEN
     6432          d_sig_gen2(i)   = d_sig_gen2(i)   + d_sig_gen(i)
     6433          d_sig_death2(i) = d_sig_death2(i) + d_sig_death(i)
     6434          d_sig_col2(i)   = d_sig_col2(i)   + d_sig_col(i)
     6435          d_sig_spread2(i)= d_sig_spread2(i)+ d_sig_spread(i)
     6436          d_sig_bnd2(i)   = d_sig_bnd2(i)   + d_sig_bnd(i)
     6437        END IF
     6438      END DO
     6439!  Bounds
     6440      DO i = 1, klon
     6441        IF (wk_adv(i)) THEN
     6442          sigmaw_targ = max(sigmaw(i),sigmad)
     6443          d_sig_bnd2(i) = d_sig_bnd2(i) + sigmaw_targ - sigmaw(i)
     6444          d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i)
     6445          sigmaw(i) = sigmaw_targ
     6446        END IF
     6447      END DO
     6448!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! wdens  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     6449!  Cumulatives
     6450      DO i = 1, klon
     6451        IF (wk_adv(i)) THEN
     6452          wdens(i) = wdens(i) + d_wdens(i)
     6453          d_wdens2(i) = d_wdens2(i) + d_wdens(i)
     6454          d_dens_gen2(i)   = d_dens_gen2(i)   + d_dens_gen(i)
     6455          d_dens_death2(i) = d_dens_death2(i) + d_dens_death(i)
     6456          d_dens_col2(i)   = d_dens_col2(i)   + d_dens_col(i)
     6457          d_dens_bnd2(i)   = d_dens_bnd2(i)   + d_dens_bnd(i)
     6458        END IF
     6459      END DO
     6460!  Bounds
     6461      DO i = 1, klon
     6462        IF (wk_adv(i)) THEN
     6463          wdens_targ = max(wdens(i),wdensmin)
     6464          d_dens_bnd2(i) = d_dens_bnd2(i) + wdens_targ - wdens(i)
     6465          d_wdens2(i) = d_wdens2(i) + wdens_targ - wdens(i)
     6466          wdens(i) = wdens_targ
     6467        END IF
     6468      END DO
     6469!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! awdens !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     6470!  Cumulatives
     6471      DO i = 1, klon
     6472        IF (wk_adv(i)) THEN
     6473          awdens(i) = awdens(i) + d_awdens(i)
     6474          d_awdens2(i) = d_awdens2(i) + d_awdens(i)
     6475        END IF
     6476      END DO
     6477!  Bounds
     6478      DO i = 1, klon
     6479        IF (wk_adv(i)) THEN
     6480          wdens_targ = min( max(awdens(i),0.), wdens(i) )
     6481          d_adens_bnd2(i) = d_adens_bnd2(i) + wdens_targ - awdens(i)
     6482          d_awdens2(i) = d_awdens2(i) + wdens_targ - awdens(i)
     6483          awdens(i) = wdens_targ
     6484        END IF
     6485      END DO
     6486!
     6487      IF (iflag_wk_pop_dyn >= 2) THEN
     6488!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! awdens again for iflag_wk_pop_dyn >= 2!!!!!!
     6489!  Cumulatives
     6490        DO i = 1, klon
     6491           IF (wk_adv(i)) THEN
     6492               d_adens_death2(i)   = d_adens_death2(i)   + d_adens_death(i)
     6493               d_adens_icol2(i)   = d_adens_icol2(i)   + d_adens_icol(i)
     6494               d_adens_acol2(i)   = d_adens_acol2(i)   + d_adens_acol(i)
     6495               d_adens_bnd2(i)   = d_adens_bnd2(i)   + d_adens_bnd(i)         
     6496           END IF
     6497        END DO
     6498!  Bounds
     6499        DO i = 1, klon
     6500           IF (wk_adv(i)) THEN
     6501               wdens_targ = min( max(awdens(i),0.), wdens(i) )
     6502               d_adens_bnd2(i) = d_adens_bnd2(i) + wdens_targ - awdens(i)
     6503               awdens(i) = wdens_targ
     6504           END IF
     6505        END DO
     6506!
     6507        IF (iflag_wk_pop_dyn == 3) THEN
     6508!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! asigmaw for iflag_wk_pop_dyn = 3!!!!!!
     6509!  Cumulatives
     6510          DO i = 1, klon
     6511             IF (wk_adv(i)) THEN
     6512                 asigmaw(i) = asigmaw(i) + d_asigmaw(i)
     6513                 d_asigmaw2(i) = d_asigmaw2(i) + d_asigmaw(i)
     6514                 d_asig_death2(i)   = d_asig_death2(i)   + d_asig_death(i)
     6515                 d_asig_spread2(i)  = d_asig_spread2(i)  + d_asig_spread(i)
     6516                 d_asig_iicol2(i)   = d_asig_iicol2(i)   + d_asig_iicol(i)
     6517                 d_asig_aicol2(i)   = d_asig_aicol2(i)   + d_asig_aicol(i)
     6518                 d_asig_bnd2(i)     = d_asig_bnd2(i)     + d_asig_bnd(i)         
     6519             END IF
     6520          END DO
     6521!  Bounds
     6522          DO i = 1, klon
     6523             IF (wk_adv(i)) THEN
     6524   !   asigmaw lower bound set to sigmad/2 in order to allow asigmaw values lower than sigmad.
     6525   !!             sigmaw_targ = min(max(asigmaw(i),sigmad),sigmaw(i))
     6526                sigmaw_targ = min(max(asigmaw(i),sigmad/2.),sigmaw(i))
     6527                d_asig_bnd2(i) = d_asig_bnd2(i) + sigmaw_targ - asigmaw(i)
     6528                d_asigmaw2(i) = d_asigmaw2(i) + sigmaw_targ - asigmaw(i)
     6529                asigmaw(i) = sigmaw_targ
     6530             END IF
     6531          END DO
     6532
     6533IF (CPPKEY_IOPHYS_WK) THEN
     6534    IF (phys_sub) THEN
     6535     CALL iophys_ecrit('wdensb',1,'wdensb','m',wdens)
     6536     CALL iophys_ecrit('awdensb',1,'awdensb','m',awdens)
     6537     CALL iophys_ecrit('sigmawb',1,'sigmawb','m',sigmaw)
     6538     CALL iophys_ecrit('asigmawb',1,'asigmawb','m',asigmaw)
     6539!
     6540     call iophys_ecrit('d_wdens2',1,'d_wdens2','',d_wdens2)
     6541     call iophys_ecrit('d_dens_gen2',1,'d_dens_gen2','',d_dens_gen2)
     6542     call iophys_ecrit('d_dens_death2',1,'d_dens_death2','',d_dens_death2)
     6543     call iophys_ecrit('d_dens_col2',1,'d_dens_col2','',d_dens_col2)
     6544     call iophys_ecrit('d_dens_bnd2',1,'d_dens_bnd2','',d_dens_bnd2)
     6545!
     6546     call iophys_ecrit('d_awdens2',1,'d_awdens2','',d_awdens2)
     6547     call iophys_ecrit('d_adens_death2',1,'d_adens_death2','',d_adens_death2)
     6548     call iophys_ecrit('d_adens_icol2',1,'d_adens_icol2','',d_adens_icol2)
     6549     call iophys_ecrit('d_adens_acol2',1,'d_adens_acol2','',d_adens_acol2)
     6550     call iophys_ecrit('d_adens_bnd2',1,'d_adens_bnd2','',d_adens_bnd2)
     6551!
     6552     CALL iophys_ecrit('d_sigmaw2',1,'d_sigmaw2','',d_sigmaw2)
     6553     CALL iophys_ecrit('d_sig_gen2',1,'d_sig_gen2','m',d_sig_gen2)
     6554     CALL iophys_ecrit('d_sig_spread2',1,'d_sig_spread2','',d_sig_spread2)
     6555     CALL iophys_ecrit('d_sig_col2',1,'d_sig_col2','',d_sig_col2)
     6556     CALL iophys_ecrit('d_sig_death2',1,'d_sig_death2','',d_sig_death2)
     6557     CALL iophys_ecrit('d_sig_bnd2',1,'d_sig_bnd2','',d_sig_bnd2)
     6558!
     6559     CALL iophys_ecrit('d_asigmaw2',1,'d_asigmaw2','',d_asigmaw2)
     6560     CALL iophys_ecrit('d_asig_spread2',1,'d_asig_spread2','m',d_asig_spread2)
     6561     CALL iophys_ecrit('d_asig_aicol2',1,'d_asig_aicol2','m',d_asig_aicol2)
     6562     CALL iophys_ecrit('d_asig_iicol2',1,'d_asig_iicol2','m',d_asig_iicol2)
     6563     CALL iophys_ecrit('d_asig_death2',1,'d_asig_death2','m',d_asig_death2)
     6564     CALL iophys_ecrit('d_asig_bnd2',1,'d_asig_bnd2','m',d_asig_bnd2)
     6565    ENDIF
     6566END IF
     6567        ENDIF ! (iflag_wk_pop_dyn == 3)
     6568      ENDIF ! (iflag_wk_pop_dyn >= 2)
     6569    ENDIF  ! (iflag_wk_pop_dyn >= 1)
     6570
     6571
     6572
     6573   Call pkupper (klon, klev, ptop, ph, p, pupper, kupper, &
     6574                    dth, hw, rho, delta_t_min, &
     6575                    ktop, wk_adv, h_zzz, ptop1, ktop1)
     6576   !! print'("pkupper APPEL ",7i6)',isubstep,int(ptop/100.),int(ptop1/100.),int(pupper/100.),ktop,ktop1,kupper
     6577
     6578    ! 5/ Set deltatw & deltaqw to 0 above kupper
     6579
     6580    DO k = 1, klev
     6581      DO i = 1, klon
     6582        IF (wk_adv(i) .AND. k>=kupper(i)) THEN
     6583          deltatw(i, k) = 0.
     6584          deltaqw(i, k) = 0.
     6585          d_deltatw2(i,k) = -deltatw0(i,k)
     6586          d_deltaqw2(i,k) = -deltaqw0(i,k)
     6587        END IF
     6588      END DO
     6589    END DO
     6590
     6591
     6592    ! -------------Cstar computation---------------------------------
     6593    DO i = 1, klon
     6594      IF (wk_adv(i)) THEN !!! nrlmd
     6595        sum_thx(i) = 0.
     6596        sum_tx(i) = 0.
     6597        sum_qx(i) = 0.
     6598        sum_thvx(i) = 0.
     6599        sum_dth(i) = 0.
     6600        sum_dq(i) = 0.
     6601        sum_dtdwn(i) = 0.
     6602        sum_dqdwn(i) = 0.
     6603
     6604        av_thx(i) = 0.
     6605        av_tx(i) = 0.
     6606        av_qx(i) = 0.
     6607        av_thvx(i) = 0.
     6608        av_dth(i) = 0.
     6609        av_dq(i) = 0.
     6610        av_dtdwn(i) = 0.
     6611        av_dqdwn(i) = 0.
     6612      END IF
     6613    END DO
     6614
     6615    ! Integrals (and wake top level number)
     6616    ! --------------------------------------
     6617
     6618    ! Initialize sum_thvx to 1st level virt. pot. temp.
     6619
     6620    DO i = 1, klon
     6621      IF (wk_adv(i)) THEN !!! nrlmd
     6622        z(i) = 1.
     6623        dz(i) = 1.
     6624        sum_thvx(i) = thx(i, 1)*(1.+epsim1*qx(i,1))*dz(i)
     6625        sum_dth(i) = 0.
     6626      END IF
     6627    END DO
     6628
     6629    DO k = 1, klev
     6630      DO i = 1, klon
     6631        IF (wk_adv(i)) THEN !!! nrlmd
     6632          dz(i) = -(max(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*RG)
     6633          IF (dz(i)>0) THEN
     6634            z(i) = z(i) + dz(i)
     6635            sum_thx(i) = sum_thx(i) + thx(i, k)*dz(i)
     6636            sum_tx(i) = sum_tx(i) + tx(i, k)*dz(i)
     6637            sum_qx(i) = sum_qx(i) + qx(i, k)*dz(i)
     6638            sum_thvx(i) = sum_thvx(i) + thx(i, k)*(1.+epsim1*qx(i,k))*dz(i)
     6639            sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i)
     6640            sum_dq(i) = sum_dq(i) + deltaqw(i, k)*dz(i)
     6641            sum_dtdwn(i) = sum_dtdwn(i) + dtdwn(i, k)*dz(i)
     6642            sum_dqdwn(i) = sum_dqdwn(i) + dqdwn(i, k)*dz(i)
     6643          END IF
     6644        END IF
     6645      END DO
     6646    END DO
     6647
     6648    DO i = 1, klon
     6649      IF (wk_adv(i)) THEN !!! nrlmd
     6650        hw0(i) = z(i)
     6651      END IF
     6652    END DO
     6653
     6654
     6655    ! - WAPE and mean forcing computation
     6656    ! ---------------------------------------
     6657
     6658    ! ---------------------------------------
     6659
     6660    ! Means
     6661
     6662    DO i = 1, klon
     6663      IF (wk_adv(i)) THEN !!! nrlmd
     6664        av_thx(i) = sum_thx(i)/hw0(i)
     6665        av_tx(i) = sum_tx(i)/hw0(i)
     6666        av_qx(i) = sum_qx(i)/hw0(i)
     6667        av_thvx(i) = sum_thvx(i)/hw0(i)
     6668        av_dth(i) = sum_dth(i)/hw0(i)
     6669        av_dq(i) = sum_dq(i)/hw0(i)
     6670        av_dtdwn(i) = sum_dtdwn(i)/hw0(i)
     6671        av_dqdwn(i) = sum_dqdwn(i)/hw0(i)
     6672
     6673        wape(i) = -RG*hw0(i)*(av_dth(i)+epsim1*(av_thx(i)*av_dq(i) + &
     6674                              av_dth(i)*av_qx(i)+av_dth(i)*av_dq(i)))/av_thvx(i)
     6675!!print *,'XXXXwake wape(i), hw0(i), av_dth(i), av_thx(i), av_dq(i), av_qx(i), av_thvx(i) ', &
     6676!!                  wape(i), hw0(i), av_dth(i), av_thx(i), av_dq(i), av_qx(i), av_thvx(i)
     6677      END IF
     6678    END DO
     6679
     6680
     6681    ! Filter out bad wakes
     6682
     6683    DO k = 1, klev
     6684      DO i = 1, klon
     6685        IF (wk_adv(i)) THEN !!! nrlmd
     6686          IF (wape(i)<=0.) THEN
     6687            deltatw(i, k) = 0.
     6688            deltaqw(i, k) = 0.
     6689            dth(i, k) = 0.
     6690            d_deltatw2(i,k) = -deltatw0(i,k)
     6691            d_deltaqw2(i,k) = -deltaqw0(i,k)
     6692          END IF
     6693        END IF
     6694      END DO
     6695    END DO
     6696
     6697! FH Tentative pour faire décroitre la fraction et la densité des poches
     6698! (de facon équivalente, comme si le rayon était inchangé) si wape<wapecut
     6699
     6700    DO k = 1, klev
     6701      DO i = 1, klon
     6702        IF (wk_adv(i)) THEN !!! nrlmd
     6703          IF (wape(i)>0. .and. wape(i)<wapecut) THEN
     6704            zzzz=wapecut/wape(i)
     6705            zzzzm1=zzzz-1.
     6706            deltatw(i, k) = zzzz*deltatw(i,k)
     6707            deltaqw(i, k) = zzzz*deltaqw(i,k)
     6708            dth(i, k) = zzzz*dth(i,k)
     6709            d_deltatw2(i,k) =  d_deltatw2(i, k) + zzzzm1*d_deltatw(i, k)
     6710            d_deltaqw2(i,k) =  d_deltaqw2(i, k) + zzzzm1*d_deltaqw(i, k)
     6711            if ( k == 1 ) print*,'zzzz',zzzz,wape,wapecut,isubstep
     6712          END IF
     6713        END IF
     6714      END DO
     6715    END DO
     6716
     6717    DO i = 1, klon
     6718      IF (wk_adv(i)) THEN !!! nrlmd
     6719        IF (wape(i)<=0.) THEN
     6720          wape(i) = 0.
     6721          cstar(i) = 0.
     6722          hw(i) = hwmin
     6723!jyg<
     6724!!          sigmaw(i) = max(sigmad, sigd_con(i))
     6725          sigmaw_targ = max(sigmad, sigd_con(i))
     6726          d_sig_bnd2(i) = d_sig_bnd2(i) + sigmaw_targ - sigmaw(i)
     6727          d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i)
     6728          sigmaw(i) = sigmaw_targ
     6729!
     6730          d_asig_bnd2(i) = d_asig_bnd2(i) + sigmaw_targ - asigmaw(i)
     6731          d_asigmaw2(i) = d_asigmaw2(i) + sigmaw_targ - asigmaw(i)
     6732          asigmaw(i) = sigmaw_targ
     6733!>jyg
     6734          fip(i) = 0.
     6735          gwake(i) = .FALSE.
     6736        ELSE
     6737          IF (wape(i)<wapecut) THEN
     6738
     6739              zzzz=wape(i)/wapecut
     6740              zzzzm1=wape(i)/wapecut-1.
     6741              wape(i)=wapecut
     6742
     6743              d_wdens2(i) = d_wdens2(i) + zzzzm1*wdens(i)
     6744              d_dens_bnd2(i) = d_dens_bnd2(i) + zzzzm1*wdens(i)
     6745              wdens(i) = zzzz*wdens(i)
     6746
     6747              d_sigmaw2(i) = d_sigmaw2(i) + zzzzm1*sigmaw(i)
     6748              d_sig_bnd2(i) = d_sig_bnd2(i) + zzzzm1*sigmaw(i)
     6749              sigmaw(i)=zzzz*sigmaw(i)
     6750
     6751              d_awdens2(i)   = d_awdens2(i)   + zzzzm1*awdens(i)
     6752              IF (iflag_wk_pop_dyn >= 2) THEN
     6753                  d_adens_bnd2(i)   = d_adens_bnd2(i)   + zzzzm1*awdens(i)
     6754              END IF
     6755              awdens(i) = zzzz*awdens(i)
     6756
     6757              IF (iflag_wk_pop_dyn == 3) THEN
     6758                 d_asigmaw2(i) = d_asigmaw2(i) + zzzzm1*asigmaw(i)
     6759                 d_asig_bnd2(i)   = d_asig_bnd2(i)   + zzzzm1*asigmaw(i)
     6760              END IF
     6761              asigmaw(i)=zzzz*asigmaw(i)
     6762           END IF
     6763           gwake(i) = .TRUE.
     6764          END IF
     6765          cstar(i) = stark*sqrt(2.*wape(i))
     6766      END IF
     6767    END DO
     6768  !
     6769  ! ------------------------------------------------------------------------
     6770  !
     6771  END DO   ! isubstep end sub-timestep loop
     6772  !
     6773  ! ------------------------------------------------------------------------
     6774  ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
     6775  ! ------------------------------------------------------------------------
     6776  !
     6777
     6778IF (CPPKEY_IOPHYS_WK) THEN
     6779    IF (.not.phys_sub) CALL iophys_ecrit('wape_b',1,'wape_b','J/kg',wape)
     6780    IF (.not.phys_sub) CALL iophys_ecrit('alpha_mod',1,'alpha_modulation','-',alpha_tot)
     6781END IF
     6782  IF (prt_level>=10) THEN
     6783    PRINT *, 'wake-5, sigmaw(igout), cstar(igout), wape(igout), ptop(igout) ', &
     6784                      sigmaw(igout), cstar(igout), wape(igout), ptop(igout)
     6785  ENDIF
     6786
     6787
     6788  ! ----------------------------------------------------------
     6789  ! Determine wake final state; recompute wape, cstar, ktop;
     6790  ! filter out bad wakes.
     6791  ! ----------------------------------------------------------
     6792
     6793  ! 2.1 - Undisturbed area and Wake integrals
     6794  ! ---------------------------------------------------------
     6795
     6796  DO i = 1, klon
     6797    ! cc nrlmd       if (wk_adv(i)) then !!! nrlmd
     6798    IF (ok_qx_qw(i)) THEN
     6799      ! cc
     6800      z(i) = 0.
     6801      sum_thx(i) = 0.
     6802      sum_tx(i) = 0.
     6803      sum_qx(i) = 0.
     6804      sum_thvx(i) = 0.
     6805      sum_dth(i) = 0.
     6806      sum_half_dth(i) = 0.
     6807      sum_dq(i) = 0.
     6808      sum_dtdwn(i) = 0.
     6809      sum_dqdwn(i) = 0.
     6810
     6811      av_thx(i) = 0.
     6812      av_tx(i) = 0.
     6813      av_qx(i) = 0.
     6814      av_thvx(i) = 0.
     6815      av_dth(i) = 0.
     6816      av_dq(i) = 0.
     6817      av_dtdwn(i) = 0.
     6818      av_dqdwn(i) = 0.
     6819
     6820      dthmin(i) = -delta_t_min
     6821    END IF
     6822  END DO
     6823  ! Potential temperatures and humidity
     6824  ! ----------------------------------------------------------
     6825
     6826  DO k = 1, klev
     6827    DO i = 1, klon
     6828      ! cc nrlmd       IF ( wk_adv(i)) THEN
     6829      IF (ok_qx_qw(i)) THEN
     6830        ! cc
     6831        rho(i, k) = p(i, k)/(RD*tb(i,k))
     6832        IF (k==1) THEN
     6833          rhoh(i, k) = ph(i, k)/(RD*tb(i,k))
     6834          zhh(i, k) = 0
     6835        ELSE
     6836          rhoh(i, k) = ph(i, k)*2./(RD*(tb(i,k)+tb(i,k-1)))
     6837          zhh(i, k) = (ph(i,k)-ph(i,k-1))/(-rhoh(i,k)*RG) + zhh(i, k-1)
     6838        END IF
     6839        thb(i, k) = tb(i, k)/ppi(i, k)
     6840        thx(i, k) = (tb(i,k)-deltatw(i,k)*sigmaw(i))/ppi(i, k)
     6841        tx(i, k) = tb(i, k) - deltatw(i, k)*sigmaw(i)
     6842        qx(i, k) = qb(i, k) - deltaqw(i, k)*sigmaw(i)
     6843        dth(i, k) = deltatw(i, k)/ppi(i, k)
     6844      END IF
     6845    END DO
     6846  END DO
     6847
     6848  ! Integrals (and wake top level number)
     6849  ! -----------------------------------------------------------
     6850
     6851  ! Initialize sum_thvx to 1st level virt. pot. temp.
     6852
     6853  DO i = 1, klon
     6854    ! cc nrlmd       IF ( wk_adv(i)) THEN
     6855    IF (ok_qx_qw(i)) THEN
     6856      ! cc
     6857      z(i) = 1.
     6858      dz(i) = 1.
     6859      dz_half(i) = 1.
     6860      sum_thvx(i) = thx(i, 1)*(1.+epsim1*qx(i,1))*dz(i)
     6861      sum_dth(i) = 0.
     6862    END IF
     6863  END DO
     6864
     6865  DO k = 1, klev
     6866    DO i = 1, klon
     6867      ! cc nrlmd       IF ( wk_adv(i)) THEN
     6868      IF (ok_qx_qw(i)) THEN
     6869        ! cc
     6870        dz(i) = -(amax1(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*RG)
     6871        dz_half(i) = -(amax1(ph(i,k+1),0.5*(ptop(i)+ph(i,1)))-ph(i,k))/(rho(i,k)*RG)
     6872        IF (dz(i)>0) THEN
     6873          z(i) = z(i) + dz(i)
     6874          sum_thx(i) = sum_thx(i) + thx(i, k)*dz(i)
     6875          sum_tx(i) = sum_tx(i) + tx(i, k)*dz(i)
     6876          sum_qx(i) = sum_qx(i) + qx(i, k)*dz(i)
     6877          sum_thvx(i) = sum_thvx(i) + thx(i, k)*(1.+epsim1*qx(i,k))*dz(i)
     6878          sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i)
     6879          sum_dq(i) = sum_dq(i) + deltaqw(i, k)*dz(i)
     6880          sum_dtdwn(i) = sum_dtdwn(i) + dtdwn(i, k)*dz(i)
     6881          sum_dqdwn(i) = sum_dqdwn(i) + dqdwn(i, k)*dz(i)
     6882!
     6883          dthmin(i) = min(dthmin(i), dth(i,k))
     6884        END IF
     6885        IF (dz_half(i)>0) THEN
     6886          sum_half_dth(i) = sum_half_dth(i) + dth(i, k)*dz_half(i)
     6887        END IF
     6888      END IF
     6889    END DO
     6890  END DO
     6891
     6892  DO i = 1, klon
     6893    ! cc nrlmd       IF ( wk_adv(i)) THEN
     6894    IF (ok_qx_qw(i)) THEN
     6895      ! cc
     6896      hw0(i) = z(i)
     6897    END IF
     6898  END DO
     6899
     6900  ! - WAPE and mean forcing computation
     6901  ! -------------------------------------------------------------
     6902
     6903  ! Means
     6904
     6905  DO i = 1, klon
     6906    ! cc nrlmd       IF ( wk_adv(i)) THEN
     6907    IF (ok_qx_qw(i)) THEN
     6908      ! cc
     6909      av_thx(i) = sum_thx(i)/hw0(i)
     6910      av_tx(i) = sum_tx(i)/hw0(i)
     6911      av_qx(i) = sum_qx(i)/hw0(i)
     6912      av_thvx(i) = sum_thvx(i)/hw0(i)
     6913      av_dth(i) = sum_dth(i)/hw0(i)
     6914      av_dq(i) = sum_dq(i)/hw0(i)
     6915      av_dtdwn(i) = sum_dtdwn(i)/hw0(i)
     6916      av_dqdwn(i) = sum_dqdwn(i)/hw0(i)
     6917
     6918      wape2(i) = -RG*hw0(i)*(av_dth(i)+epsim1*(av_thx(i)*av_dq(i) + &
     6919                             av_dth(i)*av_qx(i)+av_dth(i)*av_dq(i)))/av_thvx(i)
     6920    END IF
     6921  END DO
     6922IF (CPPKEY_IOPHYS_WK) THEN
     6923  IF (.not.phys_sub) CALL iophys_ecrit('wape2_a',1,'wape2_a','J/kg',wape2)
     6924END IF
     6925
     6926
     6927  ! Prognostic variable update
     6928  ! ------------------------------------------------------------
     6929
     6930  ! Filter out bad wakes
     6931
     6932  IF (iflag_wk_check_trgl>=1) THEN
     6933    ! Check triangular shape of dth profile
     6934    DO i = 1, klon
     6935      IF (ok_qx_qw(i)) THEN
     6936         !!print *,'XXXwake, hw0(i), dthmin(i) ', hw0(i), dthmin(i)
     6937         !!print *,'XXXwake, 2.*sum_dth(i)/(hw0(i)*dthmin(i)) ', &
     6938         !!                  2.*sum_dth(i)/(hw0(i)*dthmin(i))
     6939         !!print *,'XXXwake, sum_half_dth(i), sum_dth(i) ', &
     6940         !!                  sum_half_dth(i), sum_dth(i)
     6941        IF (iflag_wk_check_trgl<=2 .and. ((hw0(i) < 1.) .or. (dthmin(i) >= -delta_t_min)) ) THEN
     6942          wape2(i) = -1.
     6943          !!  print *,'XXXwake, rej 1'
     6944        ELSE IF (iflag_wk_check_trgl==3 .and. ((hw0(i) < 1.) .or. (dthmin(i) >= dth(i,ktop(i)))) ) THEN
     6945          wape2(i) = -1.
     6946           !! print *,'XXXwake, rej 1'
     6947        ELSE IF (iflag_wk_check_trgl==1.AND.abs(2.*sum_dth(i)/(hw0(i)*dthmin(i)) - 1.) > 0.5) THEN
     6948          wape2(i) = -1.
     6949           !! print *,'XXXwake, rej 2'
     6950        ELSE IF (abs(sum_half_dth(i)) < 0.5*abs(sum_dth(i)) ) THEN
     6951          wape2(i) = -1.
     6952           !! print *,'XXXwake, rej 3'
     6953        END IF
     6954      END IF
     6955    END DO
     6956  END IF
     6957IF (CPPKEY_IOPHYS_WK) THEN
     6958  IF (.not.phys_sub) CALL iophys_ecrit('wape2_b',1,'wape2_b','J/kg',wape2)
     6959END IF
     6960
     6961
     6962  DO k = 1, klev
     6963    DO i = 1, klon
     6964      ! cc nrlmd        IF ( wk_adv(i) .AND. wape2(i) .LT. 0.) THEN
     6965      IF (ok_qx_qw(i) .AND. wape2(i)<0.) THEN
     6966        ! cc
     6967        deltatw(i, k) = 0.
     6968        deltaqw(i, k) = 0.
     6969        dth(i, k) = 0.
     6970        d_deltatw2(i,k) = -deltatw0(i,k)
     6971        d_deltaqw2(i,k) = -deltaqw0(i,k)
     6972      END IF
     6973    END DO
     6974  END DO
     6975
     6976
     6977  DO i = 1, klon
     6978    ! cc nrlmd       IF ( wk_adv(i)) THEN
     6979    IF (ok_qx_qw(i)) THEN
     6980      ! cc
     6981      IF (wape2(i)<0.) THEN
     6982        wape2(i) = 0.
     6983        cstar2(i) = 0.
     6984        hw(i) = hwmin
     6985!jyg<
     6986!!      sigmaw(i) = amax1(sigmad, sigd_con(i))
     6987      sigmaw_targ = max(sigmad, sigd_con(i))
     6988      d_sig_bnd2(i) = d_sig_bnd2(i) + sigmaw_targ - sigmaw(i)
     6989      d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i)
     6990      sigmaw(i) = sigmaw_targ
     6991!
     6992      d_asig_bnd2(i) = d_asig_bnd2(i) + sigmaw_targ - asigmaw(i)
     6993      d_asigmaw2(i) = d_asigmaw2(i) + sigmaw_targ - asigmaw(i)
     6994      asigmaw(i) = sigmaw_targ
     6995!>jyg
     6996        fip(i) = 0.
     6997        gwake(i) = .FALSE.
     6998      ELSE
     6999        IF (prt_level>=10) PRINT *, 'wape2>0'
     7000        cstar2(i) = stark*sqrt(2.*wape2(i))
     7001        gwake(i) = .TRUE.
     7002      END IF
     7003IF (CPPKEY_IOPHYS_WK) THEN
     7004  IF (.not.phys_sub) CALL iophys_ecrit('cstar2',1,'cstar2','J/kg',cstar2)
     7005END IF
     7006    END IF  ! (ok_qx_qw(i))
     7007  END DO
     7008
     7009  DO i = 1, klon
     7010    ! cc nrlmd       IF ( wk_adv(i)) THEN
     7011    IF (ok_qx_qw(i)) THEN
     7012      ! cc
     7013      ktopw(i) = ktop(i)
     7014    END IF
     7015  END DO
     7016
     7017  DO i = 1, klon
     7018    ! cc nrlmd       IF ( wk_adv(i)) THEN
     7019    IF (ok_qx_qw(i)) THEN
     7020      ! cc
     7021      IF (ktopw(i)>0 .AND. gwake(i)) THEN
     7022
     7023        ! jyg1     Utilisation d'un h_efficace constant ( ~ feeding layer)
     7024        ! cc       heff = 600.
     7025        ! Utilisation de la hauteur hw
     7026        ! c       heff = 0.7*hw
     7027        heff(i) = hw(i)
     7028
     7029        fip(i) = 0.5*rho(i, ktopw(i))*cstar2(i)**3*heff(i)*2* &
     7030          sqrt(sigmaw(i)*wdens(i)*3.14)
     7031        fip(i) = alpk*fip(i)
     7032        ! jyg2
     7033      ELSE
     7034        fip(i) = 0.
     7035      END IF
     7036    END IF
     7037  END DO
     7038    IF (iflag_wk_pop_dyn >= 3) THEN
     7039IF (CPPKEY_IOPHYS_WK) THEN
     7040      IF (.not.phys_sub) THEN
     7041       call iophys_ecrit('d_wdens2',1,'d_wdens2','',d_wdens2)
     7042       call iophys_ecrit('d_dens_gen2',1,'d_dens_gen2','',d_dens_gen2)
     7043       call iophys_ecrit('d_dens_death2',1,'d_dens_death2','',d_dens_death2)
     7044       call iophys_ecrit('d_dens_col2',1,'d_dens_col2','',d_dens_col2)
     7045       call iophys_ecrit('d_dens_bnd2',1,'d_dens_bnd2','',d_dens_bnd2)
     7046!   
     7047       call iophys_ecrit('d_awdens2',1,'d_awdens2','',d_awdens2)
     7048       call iophys_ecrit('d_adens_death2',1,'d_adens_death2','',d_adens_death2)
     7049       call iophys_ecrit('d_adens_icol2',1,'d_adens_icol2','',d_adens_icol2)
     7050       call iophys_ecrit('d_adens_acol2',1,'d_adens_acol2','',d_adens_acol2)
     7051       call iophys_ecrit('d_adens_bnd2',1,'d_adens_bnd2','',d_adens_bnd2)
     7052!   
     7053       CALL iophys_ecrit('d_sigmaw2',1,'d_sigmaw2','',d_sigmaw2)
     7054       CALL iophys_ecrit('d_sig_gen2',1,'d_sig_gen2','m',d_sig_gen2)
     7055       CALL iophys_ecrit('d_sig_spread2',1,'d_sig_spread2','',d_sig_spread2)
     7056       CALL iophys_ecrit('d_sig_col2',1,'d_sig_col2','',d_sig_col2)
     7057       CALL iophys_ecrit('d_sig_death2',1,'d_sig_death2','',d_sig_death2)
     7058       CALL iophys_ecrit('d_sig_bnd2',1,'d_sig_bnd2','',d_sig_bnd2)
     7059!   
     7060       CALL iophys_ecrit('d_asigmaw2',1,'d_asigmaw2','',d_asigmaw2)
     7061       CALL iophys_ecrit('d_asig_spread2',1,'d_asig_spread2','m',d_asig_spread2)
     7062       CALL iophys_ecrit('d_asig_aicol2',1,'d_asig_aicol2','m',d_asig_aicol2)
     7063       CALL iophys_ecrit('d_asig_iicol2',1,'d_asig_iicol2','m',d_asig_iicol2)
     7064       CALL iophys_ecrit('d_asig_death2',1,'d_asig_death2','m',d_asig_death2)
     7065       CALL iophys_ecrit('d_asig_bnd2',1,'d_asig_bnd2','m',d_asig_bnd2)
     7066      ENDIF  ! (.not.phys_sub)
     7067END IF
     7068    ENDIF  ! (iflag_wk_pop_dyn >= 3)
     7069  ! Limitation de sigmaw
     7070
     7071  ! cc nrlmd
     7072  ! DO i=1,klon
     7073  ! IF (OK_qx_qw(i)) THEN
     7074  ! IF (sigmaw(i).GE.sigmaw_max) sigmaw(i)=sigmaw_max
     7075  ! ENDIF
     7076  ! ENDDO
     7077  ! cc
     7078
     7079  !jyg<
     7080  IF (iflag_wk_pop_dyn >= 1) THEN
     7081    DO i = 1, klon
     7082      kill_wake(i) = ((wape(i)>=wape2(i)) .AND. (wape2(i)<=wapecut)) .OR. (ktopw(i)<=2) .OR. &
     7083          .NOT. ok_qx_qw(i) .OR. (wdens(i) < wdensthreshold)
     7084!!          .NOT. ok_qx_qw(i) .OR. (wdens(i) < 2.*wdensmin)
     7085    ENDDO
     7086  ELSE  ! (iflag_wk_pop_dyn >= 1)
     7087    DO i = 1, klon
     7088      kill_wake(i) = ((wape(i)>=wape2(i)) .AND. (wape2(i)<=wapecut)) .OR. (ktopw(i)<=2) .OR. &
     7089          .NOT. ok_qx_qw(i)
     7090    ENDDO
     7091  ENDIF  ! (iflag_wk_pop_dyn >= 1)
     7092  !>jyg
     7093
     7094  DO k = 1, klev
     7095    DO i = 1, klon
     7096!!jyg      IF (((wape(i)>=wape2(i)) .AND. (wape2(i)<=wapecut)) .OR. (ktopw(i)<=2) .OR. &
     7097!!jyg          .NOT. ok_qx_qw(i)) THEN
     7098      IF (kill_wake(i)) THEN
     7099        ! cc
     7100        dtls(i, k) = 0.
     7101        dqls(i, k) = 0.
     7102        deltatw(i, k) = 0.
     7103        deltaqw(i, k) = 0.
     7104        d_deltatw2(i,k) = -deltatw0(i,k)
     7105        d_deltaqw2(i,k) = -deltaqw0(i,k)
     7106      END IF  ! (kill_wake(i))
     7107    END DO
     7108  END DO
     7109
     7110  DO i = 1, klon
     7111!!jyg    IF (((wape(i)>=wape2(i)) .AND. (wape2(i)<=wapecut)) .OR. (ktopw(i)<=2) .OR. &
     7112!!jyg        .NOT. ok_qx_qw(i)) THEN
     7113      IF (kill_wake(i)) THEN
     7114      ktopw(i) = 0
     7115      wape(i) = 0.
     7116      cstar(i) = 0.
     7117!!jyg   Outside subroutine "Wake" hw, wdens sigmaw and asigmaw are zero when there are no wakes
     7118!!      hw(i) = hwmin                       !jyg
     7119!!      sigmaw(i) = sigmad                  !jyg
     7120      hw(i) = 0.                            !jyg
     7121      fip(i) = 0.
     7122!
     7123!!      sigmaw(i) = 0.                        !jyg
     7124      sigmaw_targ = 0.
     7125      d_sig_bnd2(i) = d_sig_bnd2(i) + sigmaw_targ - sigmaw(i)
     7126!!      d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i)
     7127      d_sigmaw2(i) = sigmaw_targ - sigmaw_in(i)      ! _in = correction jyg 20220124
     7128      sigmaw(i) = sigmaw_targ
     7129!
     7130      IF (iflag_wk_pop_dyn >= 3) THEN
     7131        sigmaw_targ = 0.
     7132        d_asig_bnd2(i) = d_asig_bnd2(i) + sigmaw_targ - asigmaw(i)
     7133!!        d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i)
     7134        d_asigmaw2(i) = sigmaw_targ - asigmaw_in(i)      ! _in = correction jyg 20220124
     7135        asigmaw(i) = sigmaw_targ
     7136      ELSE
     7137        asigmaw(i) = 0.
     7138      ENDIF ! (iflag_wk_pop_dyn >= 3)
     7139!
     7140      IF (iflag_wk_pop_dyn >= 1) THEN
     7141!!        awdens(i) = 0.
     7142!!        wdens(i) = 0.
     7143        wdens_targ = 0.
     7144        d_dens_bnd2(i) = d_dens_bnd2(i) + wdens_targ - wdens(i)
     7145!!        d_wdens2(i) = wdens_targ - wdens(i)
     7146        d_wdens2(i) = wdens_targ - wdens_in(i)      ! jyg 20220916
     7147        wdens(i) = wdens_targ
     7148        wdens_targ = 0.
     7149!!jyg: bug fix : the d_adens_bnd2 computation must be before the update of awdens.
     7150        IF (iflag_wk_pop_dyn >= 2) THEN
     7151            d_adens_bnd2(i) = d_adens_bnd2(i) + wdens_targ - awdens(i)
     7152        ENDIF ! (iflag_wk_pop_dyn >= 2)
     7153!!        d_awdens2(i) = wdens_targ - awdens(i)
     7154        d_awdens2(i) = wdens_targ - awdens_in(i)    ! jyg 20220916
     7155        awdens(i) = wdens_targ
     7156!!        IF (iflag_wk_pop_dyn == 2) THEN
     7157!!            d_adens_bnd2(i) = d_adens_bnd2(i) + wdens_targ - awdens(i)
     7158!!        ENDIF ! (iflag_wk_pop_dyn == 2)
     7159      ENDIF  ! (iflag_wk_pop_dyn >= 1)
     7160    ELSE  ! (kill_wake(i))
     7161      wape(i) = wape2(i)
     7162      cstar(i) = cstar2(i)
     7163    END IF  ! (kill_wake(i))
     7164    ! c        print*,'wape wape2 ktopw OK_qx_qw =',
     7165    ! c     $          wape(i),wape2(i),ktopw(i),OK_qx_qw(i)
     7166  END DO
     7167
     7168  IF (prt_level>=10) THEN
     7169    PRINT *, 'wake-6, wape wape2 ktopw OK_qx_qw =', &
     7170                      wape(igout),wape2(igout),ktopw(igout),OK_qx_qw(igout)
     7171  ENDIF
     7172IF (CPPKEY_IOPHYS_WK) THEN
     7173  IF (.not.phys_sub) CALL iophys_ecrit('wape_c',1,'wape_c','J/kg',wape)
     7174  IF (iflag_wk_pop_dyn >= 3) THEN
     7175    IF (.not.phys_sub) THEN
     7176     CALL iophys_ecrit('fip',1,'fip','J/kg',fip)
     7177     CALL iophys_ecrit('hw',1,'hw','J/kg',hw)
     7178     CALL iophys_ecrit('ptop',1,'ptop','J/kg',ptop)
     7179     CALL iophys_ecrit('wdens',1,'wdens','J/kg',wdens)
     7180     CALL iophys_ecrit('awdens',1,'awdens','m',awdens)
     7181     CALL iophys_ecrit('sigmaw',1,'sigmaw','m',sigmaw)
     7182     CALL iophys_ecrit('asigmaw',1,'asigmaw','m',asigmaw)
     7183!
     7184     CALL iophys_ecrit('rad_wk',1,'rad_wk','J/kg',rad_wk)
     7185     CALL iophys_ecrit('arad_wk',1,'arad_wk','J/kg',arad_wk)
     7186     CALL iophys_ecrit('irad_wk',1,'irad_wk','J/kg',irad_wk)
     7187    ENDIF  ! (.not.phys_sub)
     7188  ENDIF  ! (iflag_wk_pop_dyn >= 3)
     7189END IF  !(CPPKEY_IOPHYS_WK)
     7190
     7191
     7192  ! -----------------------------------------------------------------
     7193  ! Get back to tendencies per second
     7194
     7195  DO k = 1, klev
     7196    DO i = 1, klon
     7197
     7198      ! cc nrlmd        IF ( wk_adv(i) .AND. k .LE. kupper(i)) THEN
     7199!jyg<
     7200!!      IF (ok_qx_qw(i) .AND. k<=kupper(i)) THEN
     7201      IF (ok_qx_qw(i)) THEN
     7202!>jyg
     7203        ! cc
     7204        dtls(i, k) = dtls(i, k)/dtime
     7205        dqls(i, k) = dqls(i, k)/dtime
     7206        d_deltatw2(i, k) = d_deltatw2(i, k)/dtime
     7207        d_deltaqw2(i, k) = d_deltaqw2(i, k)/dtime
     7208        d_deltat_gw(i, k) = d_deltat_gw(i, k)/dtime
     7209        ! c      print*,'k,dqls,omg,entr,detr',k,dqls(i,k),omg(i,k),entr(i,k)
     7210        ! c     $         ,death_rate(i)*sigmaw(i)
     7211      END IF
     7212    END DO
     7213  END DO
     7214!jyg<
     7215  IF (iflag_wk_pop_dyn >= 1) THEN
     7216    DO i = 1, klon
     7217        IF (ok_qx_qw(i)) THEN
     7218      d_sig_gen2(i) = d_sig_gen2(i)/dtime
     7219      d_sig_death2(i) = d_sig_death2(i)/dtime
     7220      d_sig_col2(i) = d_sig_col2(i)/dtime
     7221      d_sig_spread2(i) = d_sig_spread2(i)/dtime
     7222      d_sig_bnd2(i) = d_sig_bnd2(i)/dtime
     7223      d_sigmaw2(i) = d_sigmaw2(i)/dtime
     7224!
     7225      d_dens_gen2(i) = d_dens_gen2(i)/dtime
     7226      d_dens_death2(i) = d_dens_death2(i)/dtime
     7227      d_dens_col2(i) = d_dens_col2(i)/dtime
     7228      d_dens_bnd2(i) = d_dens_bnd2(i)/dtime
     7229      d_awdens2(i) = d_awdens2(i)/dtime
     7230      d_wdens2(i) = d_wdens2(i)/dtime
     7231        ENDIF
     7232    ENDDO
     7233    IF (iflag_wk_pop_dyn >= 2) THEN
     7234      DO i = 1, klon
     7235        IF (ok_qx_qw(i)) THEN
     7236        d_adens_death2(i) = d_adens_death2(i)/dtime
     7237        d_adens_icol2(i) = d_adens_icol2(i)/dtime
     7238        d_adens_acol2(i) = d_adens_acol2(i)/dtime
     7239        d_adens_bnd2(i) = d_adens_bnd2(i)/dtime
     7240        ENDIF
     7241      ENDDO
     7242      IF (iflag_wk_pop_dyn == 3) THEN
     7243       DO i = 1, klon
     7244          IF (ok_qx_qw(i)) THEN
     7245        d_asig_death2(i)  = d_asig_death2(i)/dtime
     7246        d_asig_iicol2(i)  = d_asig_iicol2(i)/dtime
     7247        d_asig_aicol2(i)  = d_asig_aicol2(i)/dtime
     7248        d_asig_spread2(i) = d_asig_spread2(i)/dtime
     7249        d_asig_bnd2(i) = d_asig_bnd2(i)/dtime
     7250          ENDIF
     7251       ENDDO
     7252      ENDIF ! (iflag_wk_pop_dyn == 3) 
     7253    ENDIF ! (iflag_wk_pop_dyn >= 2) 
     7254  ENDIF  ! (iflag_wk_pop_dyn >= 1)
     7255 
     7256!>jyg
     7257
     7258 RETURN
     7259END SUBROUTINE wake3
     7260
    47797261
    47807262SUBROUTINE wake_vec_modulation(nlon, nl, wk_adv, epsilon_loc, qb, d_qb, deltaqw, &
    … …  
    57368218!!          d_sigmaw(i) = max(d_sigmaw(i), sigmad-sigmaw(i))
    57378219          d_asig_death(i) = - asigmaw(i)/tau_prime(i)
    5738           d_asig_aicol(i) = (agfl(i)*iwdens(i) + igfl(i)*awdens(i))*cstar(i)*is_wk(i)
     8220!!  Bug : factor 2 omitted by mistake (bug found by Lamine Thiam)
     8221!!          d_asig_aicol(i) = (agfl(i)*iwdens(i) + igfl(i)*awdens(i))*cstar(i)*is_wk(i)
     8222          d_asig_aicol(i) = 2.*(agfl(i)*iwdens(i) + igfl(i)*awdens(i))*cstar(i)*is_wk(i)
    57398223          d_asig_iicol(i) = 2.*igfl(i)*cstar(i)*iwdens(i)*aa0
    57408224          d_asig_spread(i) = agfl(i)*cstar(i)
    … …  
    59038387
    59048388!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     8389
     8390print *,'ZZZwake_dadv_IN wk_adv(1) ', wk_adv(1)
     8391print *,'ZZZwake_dadv_IN kupper(1) ', kupper(1)
     8392print *,'ZZZwake_dadv_IN k, thw(1,k), thx(1,k) ', (k, thw(1,k), thx(1,k), k = 1,3)
     8393print *,'ZZZwake_dadv_IN k, deltomg(1,k) ', (k, deltomg(1,k), k = 1,3)
     8394print *,'ZZZwake_dadv_IN k, dp_deltomg(1,k) ', (k, dp_deltomg(1,k), k = 1,3)
     8395print *,'ZZZwake_dadv_IN sigmaw(1) ', sigmaw(1)
     8396print *,'ZZZwake_dadv_IN dsigspread(1) ', dsigspread(1)
    59058397
    59068398    entr_s(:,:) = 0.
    … …  
    62918783  ENDIF! (flag_dadv_implicit)
    62928784
     8785print *,'ZZZwake_dadv k, d_deltat_dadv(1,k) ', (k, d_deltat_dadv(1,k), k = 1,3)
     8786
    62938787    END SUBROUTINE wake_dadv
    62948788
Note: See TracChangeset for help on using the changeset viewer.