Changeset 4744 for LMDZ6


Ignore:
Timestamp:
Nov 2, 2023, 10:09:59 AM (14 months ago)
Author:
jyg
Message:

Implementation of a two radii model for wake population dynamics.
It is activated with : iflag_wk_pop_dyn=3

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

Legend:

Unmodified
Added
Removed
  • LMDZ6/trunk/libf/phylmd/add_wake_tend.F90

    r3208 r4744  
    1 SUBROUTINE add_wake_tend(zddeltat, zddeltaq, zds, zddensaw, zddensw, zoccur, text, abortphy)
     1SUBROUTINE add_wake_tend(zddeltat, zddeltaq, zds, zdas, zddensw, zddensaw, zoccur, text, abortphy)
    22!===================================================================
    33! Ajoute les tendances liées aux diverses parametrisations physiques aux
     
    99
    1010USE dimphy, ONLY: klon, klev
    11 USE phys_state_var_mod, ONLY: wake_deltat, wake_deltaq, wake_s, &
    12                               awake_dens, wake_dens
     11USE phys_state_var_mod, ONLY: wake_deltat, wake_deltaq, wake_s, awake_s, &
     12                              wake_dens, awake_dens
    1313
    1414USE print_control_mod, ONLY: prt_level
     
    1818!------------
    1919  REAL, DIMENSION(klon, klev),   INTENT (IN)         :: zddeltat, zddeltaq
    20   REAL, DIMENSION(klon),         INTENT (IN)         :: zds, zddensaw, zddensw
     20  REAL, DIMENSION(klon),         INTENT (IN)         :: zds, zdas, zddensw, zddensaw
    2121  INTEGER, DIMENSION(klon),      INTENT (IN)         :: zoccur
    2222  CHARACTER*(*),                 INTENT (IN)         :: text
     
    5555           IF (zoccur(i) .GE. 1) THEN
    5656             wake_s(i)     = wake_s(i)    + zds(i)
     57             awake_s(i)    = awake_s(i)    + zdas(i)
     58             wake_dens(i)  = wake_dens(i) + zddensw(i)
    5759             awake_dens(i) = awake_dens(i) + zddensaw(i)
    58              wake_dens(i)  = wake_dens(i) + zddensw(i)
    5960           ELSE
    6061             wake_s(i)     = 0.
     62             awake_s(i)    = 0.
     63             wake_dens(i)  = 0.
    6164             awake_dens(i) = 0.
    62              wake_dens(i)  = 0.
    6365           ENDIF   ! (zoccur(i) .GE. 1)
    6466         END DO
  • LMDZ6/trunk/libf/phylmd/calwake.F90

    r4588 r4744  
    66    dt_dwn, dq_dwn, m_dwn, m_up, dt_a, dq_a, wgen, &
    77    sigd, Cin, &
    8     wake_deltat, wake_deltaq, wake_s, awake_dens, wake_dens, &
     8    wake_deltat, wake_deltaq, wake_s, awake_s, wake_dens, awake_dens, &
    99    wake_dth, wake_h, &
    1010    wake_pe, wake_fip, wake_gfl, &
     
    1414    wake_omg, wake_dp_deltomg, &
    1515    wake_spread, wake_cstar, wake_d_deltat_gw, &
    16     wake_ddeltat, wake_ddeltaq, wake_ds, awake_ddens, wake_ddens)
     16    wake_ddeltat, wake_ddeltaq, wake_ds, awake_ds, wake_ddens, awake_ddens)
    1717  ! **************************************************************
    1818  ! *
     
    5252  ! ------------
    5353  REAL, DIMENSION(klon, klev),   INTENT (INOUT)      :: wake_deltat, wake_deltaq
    54   REAL, DIMENSION(klon),         INTENT (INOUT)      :: wake_s
    55   REAL, DIMENSION(klon),         INTENT (INOUT)      :: awake_dens, wake_dens
     54  REAL, DIMENSION(klon),         INTENT (INOUT)      :: wake_s, awake_s
     55  REAL, DIMENSION(klon),         INTENT (INOUT)      :: wake_dens, awake_dens
    5656  ! Output
    5757  ! ------
     
    7070  REAL, DIMENSION(klon),         INTENT (OUT)        :: wake_cstar
    7171  REAL, DIMENSION(klon, klev),   INTENT (OUT)        :: wake_ddeltat, wake_ddeltaq
    72   REAL, DIMENSION(klon),         INTENT (OUT)        :: wake_ds, awake_ddens, wake_ddens
     72  REAL, DIMENSION(klon),         INTENT (OUT)        :: wake_ds, awake_ds, wake_ddens, awake_ddens
    7373
    7474
     
    9191  REAL, DIMENSION(klon, klev)                        :: tx, qx
    9292  REAL, DIMENSION(klon)                              :: hw, wape, fip, gfl
    93   REAL, DIMENSION(klon)                              :: sigmaw, awdens, wdens
     93  REAL, DIMENSION(klon)                              :: sigmaw, asigmaw, wdens, awdens
    9494  REAL, DIMENSION(klon, klev)                        :: omgbdth
    9595  REAL, DIMENSION(klon, klev)                        :: dp_omgb
     
    102102  REAL, DIMENSION(klon, klev)                        :: d_deltat_gw
    103103  REAL, DIMENSION(klon, klev)                        :: d_deltatw, d_deltaqw
    104   REAL, DIMENSION(klon)                              :: d_sigmaw, d_awdens, d_wdens
     104  REAL, DIMENSION(klon)                              :: d_sigmaw, d_asigmaw, d_wdens, d_awdens
    105105
    106106  REAL                                               :: rdcp
     
    108108
    109109  IF (prt_level >= 10) THEN
    110     print *, '-> calwake, wake_s, wgen input ', wake_s(1), wgen(1)
     110    print *, '-> calwake, wake_s, awake_s, wgen input ', wake_s(1), awake_s(1), wgen(1)
    111111  ENDIF
    112112
     
    150150d_deltaqw(:,:) = 0.
    151151d_sigmaw(:) = 0.
     152d_asigmaw(:) = 0.
     153d_wdens(:) = 0.
    152154d_awdens(:) = 0.
    153 d_wdens(:) = 0.
    154155!
    155156
     
    179180
    180181  DO i = 1, klon
    181     sigmaw(i) = wake_s(i)
    182   END DO
    183 
    184   DO i = 1, klon
     182    sigmaw(i)  = wake_s(i)
     183    asigmaw(i) = awake_s(i)
     184  END DO
     185
     186  DO i = 1, klon
     187    wdens(i) = max(0., wake_dens(i))
    185188    awdens(i) = max(0., awake_dens(i))
    186     wdens(i) = max(0., wake_dens(i))
    187189  END DO
    188190
     
    215217    dtdwn, dqdwn, amdwn, amup, dta, dqa, wgen, &
    216218    sigd0, Cin, &
    217     dtw, dqw, sigmaw, awdens, wdens, &                      ! state variables
     219    dtw, dqw, sigmaw, asigmaw, wdens, awdens, &                      ! state variables
    218220    dth, hw, wape, fip, gfl, &
    219221    dtls, dqls, ktopw, omgbdth, dp_omgb, tx, qx, &
    220222    dtke, dqke, omg, dp_deltomg, spread, cstar, &
    221223    d_deltat_gw, &
    222     d_deltatw, d_deltaqw, d_sigmaw, d_awdens, d_wdens)      ! tendencies
     224    d_deltatw, d_deltaqw, d_sigmaw, d_asigmaw, d_wdens, d_awdens)      ! tendencies
    223225
    224226!
     
    281283    IF (ktopw(i)>0) THEN
    282284      wake_ds(i) = d_sigmaw(i)*dtime
     285      awake_ds(i) = d_asigmaw(i)*dtime
    283286      awake_ddens(i) = d_awdens(i)*dtime
    284287      wake_ddens(i) = d_wdens(i)*dtime
    285288    ELSE
    286       wake_ds(i)   = -wake_s(i)
    287       wake_ddens(i)= -wake_dens(i)
     289      wake_ds(i)    = -wake_s(i)
     290      awake_ds(i)   = -awake_s(i)
     291      wake_ddens(i) = -wake_dens(i)
     292      awake_ddens(i)= -awake_dens(i)
    288293    END IF
    289294  END DO
     
    306311    DO i = 1, klon
    307312      wake_s(i) = sigmaw(i)
     313      awake_s(i) = asigmaw(i)
    308314      awake_dens(i) = awdens(i)
    309315      wake_dens(i) = wdens(i)
     
    320326  ENDIF  ! (first)
    321327!>jyg
     328  IF (prt_level >= 10) THEN
     329    print *, 'calwake ->, wake_s, awake_s ', wake_s(1), awake_s(1)
     330  ENDIF
    322331
    323332  RETURN
    324333END SUBROUTINE calwake
    325334
     335
  • LMDZ6/trunk/libf/phylmd/dyn1d/old_lmdz1d.F90

    r4593 r4744  
    1515       solsw, solswfdiff, t_ancien, q_ancien, u_ancien, v_ancien, rneb_ancien, &
    1616       wake_delta_pbl_TKE, delta_tsurf, wake_fip, wake_pe, &
    17        wake_deltaq, wake_deltat, wake_s, wake_dens, &
     17       wake_deltaq, wake_deltat, wake_s, awake_s, wake_dens, &
    1818       awake_dens, cv_gen, wake_cstar, &
    1919       zgam, zmax0, zmea, zpic, zsig, &
     
    905905        wake_pe = 0.
    906906        wake_s = 0.
     907        awake_s = 0.
    907908        wake_dens = 0.
    908909        awake_dens = 0.
     
    938939! t_ancien,q_ancien,,frugs(:,is_oce),clwcon(:,1),rnebcon(:,1),ratqs(:,1)
    939940! run_off_lic_0,pbl_tke(:,1:klev,nsrf), zmax0,f0,sig1,w01
    940 ! wake_deltat,wake_deltaq,wake_s,wake_dens,awake_dens,cv_gen,wake_cstar,
     941! wake_deltat,wake_deltaq,wake_s,awake_s,wake_dens,awake_dens,cv_gen,wake_cstar,
    941942! wake_fip,wake_delta_pbl_tke(:,1:klev,nsrf)
    942943!
  • LMDZ6/trunk/libf/phylmd/dyn1d/scm.F90

    r4593 r4744  
    1111       solsw, solswfdiff, t_ancien, q_ancien, u_ancien, v_ancien, rneb_ancien, &
    1212       wake_delta_pbl_TKE, delta_tsurf, wake_fip, wake_pe, &
    13        wake_deltaq, wake_deltat, wake_s, wake_dens, &
     13       wake_deltaq, wake_deltat, wake_s, awake_s, wake_dens, &
    1414       awake_dens, cv_gen, wake_cstar, &
    1515       zgam, zmax0, zmea, zpic, zsig, &
     
    663663        wake_pe = 0.
    664664        wake_s = 0.
     665        awake_s = 0.
    665666        wake_dens = 0.
    666667        awake_dens = 0.
     
    696697! t_ancien,q_ancien,,frugs(:,is_oce),clwcon(:,1),rnebcon(:,1),ratqs(:,1)
    697698! run_off_lic_0,pbl_tke(:,1:klev,nsrf), zmax0,f0,sig1,w01
    698 ! wake_deltat,wake_deltaq,wake_s,wake_dens,awake_dens,cv_gen,wake_cstar,
     699! wake_deltat,wake_deltaq,wake_s,awake_s,wake_dens,awake_dens,cv_gen,wake_cstar,
    699700! wake_fip,wake_delta_pbl_tke(:,1:klev,nsrf)
    700701!
  • LMDZ6/trunk/libf/phylmd/lmdz_wake.F90

    r4695 r4744  
    99                dtdwn, dqdwn, amdwn, amup, dta, dqa, wgen, &
    1010                sigd_con, Cin, &
    11                 deltatw, deltaqw, sigmaw, awdens, wdens, &                  ! state variables           
     11                deltatw, deltaqw, sigmaw, asigmaw, wdens, awdens, &                  ! state variables           
    1212                dth, hw, wape, fip, gfl, &
    1313                dtls, dqls, ktopw, omgbdth, dp_omgb, tu, qu, &
    1414                dtke, dqke, omg, dp_deltomg, wkspread, cstar, &
    1515                d_deltat_gw, &                                                      ! tendencies
    16                 d_deltatw2, d_deltaqw2, d_sigmaw2, d_awdens2, d_wdens2)             ! tendencies
     16                d_deltatw2, d_deltaqw2, d_sigmaw2, d_asigmaw2, d_wdens2, d_awdens2)             ! tendencies
    1717
    1818
     
    3232  USE lmdz_wake_ini , ONLY : crep_upper, crep_sol, tau_cv, rzero, aa0, flag_wk_check_trgl
    3333  USE lmdz_wake_ini , ONLY : ok_bug_gfl
    34   USE lmdz_wake_ini , ONLY : iflag_wk_act, iflag_wk_check_trgl, iflag_wk_pop_dyn, wdensmin
     34  USE lmdz_wake_ini , ONLY : iflag_wk_act, iflag_wk_check_trgl, iflag_wk_pop_dyn, wdensinit, wdensthreshold
    3535  USE lmdz_wake_ini , ONLY : sigmad, hwmin, wapecut, cstart, sigmaw_max, dens_rate, epsilon_loc
    3636  USE lmdz_wake_ini , ONLY : iflag_wk_profile
     37  USE lmdz_wake_ini , ONLY : smallestreal
    3738
    3839
     
    4950  ! deltaqw    : specific humidity difference between wake and off-wake regions
    5051  ! sigmaw     : fractional area covered by wakes.
     52  ! asigmaw    : fractional area covered by active wakes.
    5153  ! wdens      : number of wakes per unit area
     54  ! awdens     : number of active wakes per unit area
    5255
    5356  ! Variable de sortie :
     
    7174  ! deltaqw     : updated humidity difference (q_w-q_u).
    7275  ! sigmaw      : updated wake fractional area.
     76  ! asigmaw     : updated active wake fractional area.
    7377  ! d_deltat_gw : delta T tendency due to GW
    7478
     
    154158  REAL, DIMENSION (klon, klev),     INTENT(INOUT)       :: deltatw, deltaqw
    155159  REAL, DIMENSION (klon),           INTENT(INOUT)       :: sigmaw
     160  REAL, DIMENSION (klon),           INTENT(INOUT)       :: asigmaw
     161  REAL, DIMENSION (klon),           INTENT(INOUT)       :: wdens
    156162  REAL, DIMENSION (klon),           INTENT(INOUT)       :: awdens
    157   REAL, DIMENSION (klon),           INTENT(INOUT)       :: wdens
    158163
    159164  ! Sorties
     
    173178  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: d_deltat_gw
    174179  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: d_deltatw2, d_deltaqw2
    175   REAL, DIMENSION (klon),           INTENT(OUT)         :: d_sigmaw2, d_awdens2, d_wdens2
     180  REAL, DIMENSION (klon),           INTENT(OUT)         :: d_sigmaw2, d_asigmaw2, d_wdens2, d_awdens2
    176181
    177182  ! Variables internes
     
    191196  REAL, DIMENSION (klon, klev)                          :: deltaqw0
    192197  REAL, DIMENSION (klon, klev)                          :: tenv, qe
    193 !!  REAL, DIMENSION (klon)                                :: sigmaw1
    194198
    195199  ! Variables liees a la dynamique de population 1
     
    202206  REAL, DIMENSION(klon)                                 :: cont_fact 
    203207 
     208  ! Variables liees a la dynamique de population 3
     209  REAL, DIMENSION(klon)                                 :: arad_wk, irad_wk
     210 
    204211!!  REAL, DIMENSION(klon)                                 :: d_sig_gen, d_sig_death, d_sig_col
    205212  REAL, DIMENSION(klon)                                 :: wape1_act, wape2_act
     
    209216  ! Some components of the tendencies of state variables 
    210217  REAL, DIMENSION (klon)                                :: d_sig_gen2, d_sig_death2, d_sig_col2, d_sig_spread2, d_sig_bnd2
     218  REAL, DIMENSION (klon)                                :: d_asig_death2, d_asig_aicol2, d_asig_iicol2, d_asig_spread2, d_asig_bnd2
    211219  REAL, DIMENSION (klon)                                :: d_dens_gen2, d_dens_death2, d_dens_col2, d_dens_bnd2
    212220  REAL, DIMENSION (klon)                                :: d_adens_death2, d_adens_icol2, d_adens_acol2, d_adens_bnd2
     
    232240  REAL, DIMENSION (klon, klev)                          :: d_deltatw, d_deltaqw
    233241  REAL, DIMENSION (klon, klev)                          :: d_tenv, d_qe
    234   REAL, DIMENSION (klon)                                :: d_awdens, d_wdens, d_sigmaw
     242  REAL, DIMENSION (klon)                                :: d_wdens, d_awdens, d_sigmaw, d_asigmaw
    235243  REAL, DIMENSION (klon)                                :: d_sig_gen, d_sig_death, d_sig_col, d_sig_spread, d_sig_bnd
     244  REAL, DIMENSION (klon)                                :: d_asig_death, d_asig_aicol, d_asig_iicol, d_asig_spread, d_asig_bnd
    236245  REAL, DIMENSION (klon)                                :: d_dens_gen, d_dens_death, d_dens_col, d_dens_bnd
    237246  REAL, DIMENSION (klon)                                :: d_adens_death, d_adens_icol, d_adens_acol, d_adens_bnd
     247  REAL, DIMENSION (klon)                                :: agfl              !! gust front length of active wakes
     248                                                                             !!  per unit area
    238249  REAL, DIMENSION (klon)                                :: alpha, alpha_tot
    239250  REAL, DIMENSION (klon)                                :: q0_min, q1_min
     
    242253  ! Autres variables internes
    243254  INTEGER                                               ::isubstep, k, i, igout
     255
     256  REAL                                                  :: wdensmin
    244257
    245258  REAL                                                  :: sigmaw_targ
     
    291304  REAL, DIMENSION (klon, klev)                          :: detr
    292305
    293   REAL, DIMENSION(klon)                                 :: sigmaw_in             ! pour les prints
    294   REAL, DIMENSION(klon)                                 :: awdens_in, wdens_in   ! pour les prints
     306  REAL, DIMENSION(klon)                                 :: sigmaw_in, asigmaw_in ! pour les prints
     307  REAL, DIMENSION(klon)                                 :: wdens_in, awdens_in   ! pour les prints
     308
     309!!!  LOGICAL                                               :: phys_sub=.true.
     310  LOGICAL                                               :: phys_sub=.false.
     311
     312  LOGICAL                                               :: first_call=.true.
    295313
    296314  ! -------------------------------------------------------------------------
     
    302320  !  Provisionnal; to be suppressed when f_shear is parameterized
    303321  f_shear(:) = 1.       ! 0. for strong shear, 1. for weak shear
    304 
    305322
    306323  ! Configuration de coefgw,stark,wdens (22/02/06 by YU Jingmei)
     
    321338  ! alpk = 0.5
    322339  ! alpk = 0.05
    323 !print *,'XXXX dtime input ', dtime
     340!
    324341 igout = klon/2+1/klon
     342!
     343!   sub-time-stepping parameters
     344  nsub = 10
     345  dtimesub = dtime/nsub
     346!
     347IF (first_call) THEN
     348!!#define IOPHYS_WK
     349#undef IOPHYS_WK
     350#ifdef IOPHYS_WK
     351  IF (phys_sub) THEN
     352    call iophys_ini(dtimesub)
     353  ELSE
     354    call iophys_ini(dtime)
     355  ENDIF
     356#endif
     357  first_call = .false.
     358ENDIF   !(first_call)
    325359
    326360 IF (iflag_wk_pop_dyn == 0) THEN
     
    335369  !>jyg
    336370 ENDIF  ! (iflag_wk_pop_dyn == 0)
     371!
     372 IF (iflag_wk_pop_dyn >=1) THEN
     373   IF (iflag_wk_pop_dyn == 3) THEN
     374     wdensmin = wdensthreshold
     375   ELSE
     376     wdensmin = wdensinit
     377   ENDIF
     378 ENDIF
    337379
    338380  ! print*,'stark',stark
     
    385427      d_deltaqw2(:,:) = 0.
    386428
     429      d_sig_gen2(:)   = 0.
     430      d_sig_death2(:) = 0.
     431      d_sig_col2(:)   = 0.
     432      d_sig_spread2(:)= 0.
     433      d_asig_death2(:) = 0.
     434      d_asig_iicol2(:) = 0.
     435      d_asig_aicol2(:) = 0.
     436      d_asig_spread2(:)= 0.
     437      d_asig_bnd2(:) = 0.
     438      d_asigmaw2(:) = 0.
     439!
     440      d_dens_gen2(:)   = 0.
     441      d_dens_death2(:) = 0.
     442      d_dens_col2(:)   = 0.
     443      d_dens_bnd2(:) = 0.
     444      d_wdens2(:) = 0.
     445      d_adens_bnd2(:) = 0.
     446      d_awdens2(:) = 0.
     447      d_adens_death2(:) = 0.
     448      d_adens_icol2(:) = 0.
     449      d_adens_acol2(:) = 0.
     450
    387451      IF (iflag_wk_act == 0) THEN
    388452        act(:) = 0.
     
    394458!!   sigmaw_in(i) = sigmaw(i)
    395459!!  END DO
    396    sigmaw_in(:) = sigmaw(:)
     460   sigmaw_in(:)  = sigmaw(:)
     461   asigmaw_in(:) = asigmaw(:)
    397462!>jyg
    398463!
     
    416481      d_dens_col2(i)   = 0.
    417482      d_awdens2(i) = 0.
    418 !
    419       wdens_targ = max(wdens(i),wdensmin)
    420       d_dens_bnd2(i) = wdens_targ - wdens(i)
    421       d_wdens2(i) = wdens_targ - wdens(i)
    422       wdens(i) = wdens_targ
    423     END DO
    424     IF (iflag_wk_pop_dyn == 2) THEN
    425        DO i = 1, klon 
    426           d_adens_death2(i) = 0.
    427           d_adens_icol2(i) = 0.
    428           d_adens_acol2(i) = 0.
    429 !     
    430           wdens_targ = min(max(awdens(i),0.),wdens(i))
    431           d_adens_bnd2(i) = wdens_targ - awdens(i)
    432           d_awdens2(i) = wdens_targ - awdens(i)
    433           awdens(i) = wdens_targ
    434        END DO
    435     ENDIF ! (iflag_wk_pop_dyn == 2)
     483      IF (wdens(i) < wdensthreshold) THEN
     484  !!      wdens_targ = max(wdens(i),wdensmin)
     485        wdens_targ = max(wdens(i),wdensinit)
     486        d_dens_bnd2(i) = wdens_targ - wdens(i)
     487        d_wdens2(i) = wdens_targ - wdens(i)
     488        wdens(i) = wdens_targ
     489      ELSE
     490        d_dens_bnd2(i) = 0.
     491        d_wdens2(i) = 0.
     492      ENDIF  !! (wdens(i) < wdensthreshold)
     493    END DO
     494    IF (iflag_wk_pop_dyn >= 2) THEN
     495      DO i = 1, klon 
     496        IF (awdens(i) < wdensthreshold) THEN
     497!!          wdens_targ = min(max(awdens(i),wdensmin),wdens(i))
     498            wdens_targ = min(max(awdens(i),wdensinit),wdens(i))
     499            d_adens_bnd2(i) = wdens_targ - awdens(i)
     500            d_awdens2(i) = wdens_targ - awdens(i)
     501            awdens(i) = wdens_targ
     502        ELSE
     503            wdens_targ = min(awdens(i), wdens(i))
     504            d_adens_bnd2(i) = wdens_targ - awdens(i)
     505            d_awdens2(i) = wdens_targ - awdens(i)
     506            awdens(i) = wdens_targ
     507        ENDIF
     508      END DO
     509    ENDIF ! (iflag_wk_pop_dyn >= 2)
    436510  ELSE 
    437511    DO i = 1, klon
     
    442516!
    443517  DO i = 1, klon
    444     ! c      sigmaw(i) = amax1(sigmaw(i),sigd_con(i))
    445 !jyg<
    446 !!    sigmaw(i) = amax1(sigmaw(i), sigmad)
    447 !!    sigmaw(i) = amin1(sigmaw(i), 0.99)
    448     d_sig_gen2(i)   = 0.
    449     d_sig_death2(i) = 0.
    450     d_sig_col2(i)   = 0.
    451     d_sig_spread2(i)= 0.
    452518    sigmaw_targ = min(max(sigmaw(i), sigmad),0.99)
    453519    d_sig_bnd2(i) = sigmaw_targ - sigmaw(i)
    454520    d_sigmaw2(i) = sigmaw_targ - sigmaw(i)
    455 !  print *,'XXXX1 d_sigmaw2(i), sigmaw(i) ', d_sigmaw2(i), sigmaw(i)
    456521    sigmaw(i) = sigmaw_targ
    457 !>jyg
    458   END DO
     522  END DO
     523!
     524  IF (iflag_wk_pop_dyn == 3) THEN
     525     DO i = 1, klon 
     526        IF ((wdens(i)-awdens(i)) <= smallestreal) THEN
     527          sigmaw_targ = sigmaw(i)
     528        ELSE
     529          sigmaw_targ = min(max(asigmaw(i),sigmad),sigmaw(i))
     530        ENDIF
     531        d_asig_bnd2(i) = sigmaw_targ - asigmaw(i)
     532        d_asigmaw2(i) = sigmaw_targ - asigmaw(i)
     533        asigmaw(i) = sigmaw_targ
     534     END DO
     535  ENDIF ! (iflag_wk_pop_dyn == 3)
    459536
    460537  wape(:) = 0.
    461538  wape2(:) = 0.
    462539  d_sigmaw(:) = 0.
     540  d_asigmaw(:) = 0.
    463541  ktopw(:) = 0
    464542!
     
    810888
    811889  END DO
     890#ifdef IOPHYS_WK
     891  IF (.not.phys_sub) CALL iophys_ecrit('wape_a',1,'wape_a','J/kg',wape)
     892#endif
    812893
    813894  ! 2.2 Prognostic variable update
     
    830911  DO i = 1, klon
    831912    IF (wape(i)<0.) THEN
    832       wape(i) = 0.
    833       cstar(i) = 0.
    834       hw(i) = hwmin
    835 !jyg<
    836913!!      sigmaw(i) = amax1(sigmad, sigd_con(i))
    837914      sigmaw_targ = max(sigmad, sigd_con(i))
    838915      d_sig_bnd2(i) = d_sig_bnd2(i) + sigmaw_targ - sigmaw(i)
    839916      d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i)
    840 !  print *,'XXXX2 d_sigmaw2(i), sigmaw(i) ', d_sigmaw2(i), sigmaw(i)
    841917      sigmaw(i) = sigmaw_targ
    842 !>jyg
     918    ENDIF  !!  (wape(i)<0.)
     919  ENDDO
     920  !
     921  IF (iflag_wk_pop_dyn == 3) THEN
     922    DO i = 1, klon
     923      IF (wape(i)<0.) THEN
     924        sigmaw_targ = max(sigmad, sigd_con(i))
     925        d_asig_bnd2(i) = d_asig_bnd2(i) + sigmaw_targ - asigmaw(i)
     926        d_asigmaw2(i) = d_asigmaw2(i) + sigmaw_targ - asigmaw(i)
     927        asigmaw(i) = sigmaw_targ
     928      ENDIF  !!  (wape(i)<0.)
     929    ENDDO
     930  ENDIF  !!  (iflag_wk_pop_dyn == 3)
     931
     932  DO i = 1, klon
     933    IF (wape(i)<0.) THEN
     934      wape(i) = 0.
     935      cstar(i) = 0.
     936      hw(i) = hwmin
    843937      fip(i) = 0.
    844938      gwake(i) = .FALSE.
     
    849943    END IF
    850944  END DO
    851 
     945  !
    852946
    853947  ! Check qx and qw positivity
     
    883977  ! -----------------
    884978
    885   nsub = 10
    886   dtimesub = dtime/nsub
    887 
     979!    nsub and dtimesub definitions moved to begining of routine.
     980!!  nsub = 10
     981!!  dtimesub = dtime/nsub
    888982
    889983 
    890   ! ------------------------------------------------------------
     984  ! ------------------------------------------------------------------------
     985  ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
     986  ! ------------------------------------------------------------------------
     987  !
    891988  DO isubstep = 1, nsub
    892     ! ------------------------------------------------------------
    893   CALL pkupper (klon, klev, ptop, ph, pupper, kupper)
     989  !
     990  ! ------------------------------------------------------------------------
     991  !
     992    CALL pkupper (klon, klev, ptop, ph, pupper, kupper)
    894993 
    895994    !print*, 'ptop, pupper, ktop, kupper', ptop, pupper, ktop, kupper
     
    9431042    END DO
    9441043
    945     IF (ok_bug_gfl) THEN
     1044    IF (iflag_wk_pop_dyn == 0 .AND. ok_bug_gfl) THEN
    9461045!!--------------------------------------------------------
    9471046!!Bug : computing gfl and rad_wk before changing sigmaw
     1047!!      This bug exists only for iflag_wk_pop_dyn=0. Otherwise, gfl and rad_wk
     1048!!      are computed within  wake_popdyn
    9481049!!--------------------------------------------------------
    9491050      DO i = 1, klon
     
    9531054        END IF
    9541055      END DO
    955     ENDIF   ! (ok_bug_gfl)
     1056    ENDIF   ! (iflag_wk_pop_dyn == 0 .AND. ok_bug_gfl)
     1057!!--------------------------------------------------------
    9561058
    9571059    DO i = 1, klon
     
    9591061        sigmaw_targ = min(sigmaw(i), sigmaw_max)
    9601062        d_sig_bnd2(i) = d_sig_bnd2(i) + sigmaw_targ - sigmaw(i)
    961         d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i)
     1063        d_sigmaw2(i)  = d_sigmaw2(i) + sigmaw_targ - sigmaw(i)
    9621064        sigmaw(i) = sigmaw_targ
    9631065      END IF
    9641066    END DO
    9651067
    966     IF (.NOT.ok_bug_gfl) THEN
     1068    IF (iflag_wk_pop_dyn == 0 .AND. .NOT.ok_bug_gfl) THEN
    9671069!!--------------------------------------------------------
    9681070!!Fix : computing gfl and rad_wk after changing sigmaw
     
    9741076        END IF
    9751077      END DO
    976     ENDIF   ! (.NOT.ok_bug_gfl)
    977 
     1078    ENDIF   ! (iflag_wk_pop_dyn == 0 .AND. .NOT.ok_bug_gfl)
     1079!!--------------------------------------------------------
     1080
     1081    IF (iflag_wk_pop_dyn >= 1) THEN
     1082  !  The variable "death_rate" is significant only when iflag_wk_pop_dyn = 0.
     1083  !  Here, it has to be set to zero.
     1084      death_rate(:) = 0.
     1085    ENDIF
     1086 
     1087    IF (iflag_wk_pop_dyn >= 3) THEN
     1088      DO i = 1, klon
     1089        IF (wk_adv(i)) THEN
     1090         sigmaw_targ = min(asigmaw(i), sigmaw_max)
     1091         d_asig_bnd2(i) = d_asig_bnd2(i) + sigmaw_targ - asigmaw(i)
     1092         d_asigmaw2(i)  = d_asigmaw2(i)  + sigmaw_targ - asigmaw(i)
     1093         asigmaw(i) = sigmaw_targ
     1094        ENDIF
     1095      ENDDO
     1096    ENDIF
     1097 
     1098!!--------------------------------------------------------
     1099!!--------------------------------------------------------
    9781100    IF (iflag_wk_pop_dyn == 1) THEN
    979  
     1101  !
    9801102     CALL wake_popdyn_1 (klon, klev, dtime, cstar, tau_wk_inv, wgen, wdens, awdens, sigmaw, &
     1103                  wdensmin, &
    9811104                  dtimesub, gfl, rad_wk, f_shear, drdt_pos, &
    9821105                  d_awdens, d_wdens, d_sigmaw, &
     
    9871110                  d_wdens_targ, d_sigmaw_targ)
    9881111                     
    989 !    The variable "death_rate" is significant only when iflag_wk_pop_dyn = 0.
    990 !    Here, it has to be set to zero.
    991       death_rate(:) = 0.
    9921112   
     1113!!--------------------------------------------------------
    9931114    ELSEIF (iflag_wk_pop_dyn == 2) THEN
     1115  !
    9941116     CALL wake_popdyn_2 ( klon, klev, wk_adv, dtimesub, wgen, &
    995                              sigmaw, wdens, awdens, &   !! states variables
     1117                             wdensmin, &
     1118                             sigmaw, wdens, awdens, &   !! state variables
    9961119                             gfl, cstar, cin, wape, rad_wk, &
    997                              d_sigmaw, d_wdens, d_awdens, &  !! tendences                             
     1120                             d_sigmaw, d_wdens, d_awdens, &  !! tendencies                             
    9981121                             cont_fact, &
    9991122                             d_sig_gen, d_sig_death, d_sig_col, d_sig_spread, d_sig_bnd, &
    10001123                             d_dens_gen, d_dens_death, d_dens_col, d_dens_bnd, &
    10011124                             d_adens_death, d_adens_icol, d_adens_acol, d_adens_bnd )
    1002      death_rate(:) = 0.
    10031125sigmaw=sigmaw-d_sigmaw
    10041126wdens=wdens-d_wdens
    10051127awdens=awdens-d_awdens
     1128
     1129!!--------------------------------------------------------
     1130    ELSEIF (iflag_wk_pop_dyn == 3) THEN
     1131#ifdef IOPHYS_WK
     1132    IF (phys_sub) THEN
     1133     CALL iophys_ecrit('ptop',1,'ptop','Pa',ptop)
     1134     CALL iophys_ecrit('sigmaw',1,'sigmaw','',sigmaw)
     1135     CALL iophys_ecrit('asigmaw',1,'asigmaw','',asigmaw)
     1136     CALL iophys_ecrit('wdens',1,'wdens','1/m2',wdens)
     1137     CALL iophys_ecrit('awdens',1,'awdens','1/m2',awdens)
     1138     CALL iophys_ecrit('rad_wk',1,'rad_wk','m',rad_wk)
     1139     CALL iophys_ecrit('arad_wk',1,'arad_wk','m',arad_wk)
     1140     CALL iophys_ecrit('irad_wk',1,'irad_wk','m',irad_wk)
     1141    ENDIF
     1142#endif
     1143  !
     1144     CALL wake_popdyn_3 ( klon, klev, phys_sub, wk_adv, dtimesub, wgen, &
     1145                             wdensmin, &
     1146                             sigmaw, asigmaw, wdens, awdens, &   !! state variables
     1147                             gfl, agfl, cstar, cin, wape, &
     1148                             rad_wk, arad_wk, irad_wk, &
     1149                             d_sigmaw, d_asigmaw, d_wdens, d_awdens, &  !! tendencies                             
     1150                             d_sig_gen, d_sig_death, d_sig_col, d_sig_spread, d_sig_bnd, &
     1151                             d_asig_death, d_asig_aicol, d_asig_iicol, d_asig_spread, d_asig_bnd, &
     1152                             d_dens_gen, d_dens_death, d_dens_col, d_dens_bnd, &
     1153                             d_adens_death, d_adens_icol, d_adens_acol, d_adens_bnd )
     1154sigmaw=sigmaw-d_sigmaw
     1155asigmaw=asigmaw-d_asigmaw
     1156wdens=wdens-d_wdens
     1157awdens=awdens-d_awdens
    10061158   
     1159!!--------------------------------------------------------
    10071160    ELSEIF (iflag_wk_pop_dyn == 0) THEN
    10081161   
     
    10111164      DO i = 1, klon
    10121165        IF (wk_adv(i)) THEN
     1166
    10131167          ! cc nrlmd          Introduction du taux de mortalite des poches et
    10141168          ! test sur sigmaw_max=0.4
     
    10361190      END DO
    10371191
    1038     ENDIF   !  (iflag_wk_pop_dyn >= 1)
    1039 
    1040 
     1192    ENDIF   !  (iflag_wk_pop_dyn == 1) ... ELSEIF (iflag_wk_pop_dyn == 0)
     1193!!--------------------------------------------------------
     1194!!--------------------------------------------------------
     1195
     1196#ifdef IOPHYS_WK
     1197    IF (phys_sub) THEN
     1198     CALL iophys_ecrit('wdensa',1,'wdensa','m',wdens)
     1199     CALL iophys_ecrit('awdensa',1,'awdensa','m',awdens)
     1200     CALL iophys_ecrit('sigmawa',1,'sigmawa','m',sigmaw)
     1201     CALL iophys_ecrit('asigmawa',1,'asigmawa','m',asigmaw)
     1202    ENDIF
     1203#endif
    10411204    ! calcul de la difference de vitesse verticale poche - zone non perturbee
    10421205    ! IM 060208 differences par rapport au code initial; init. a 0 dp_deltomg
     
    12041367        IF (wk_adv(i) .AND. k<=kupper(i)+1) THEN
    12051368          dth(i, k) = deltatw(i, k)/ppi(i, k)
    1206 ! print *, 'VVVVwake k, the(i,k), dth(i,k), sigmaw(i) ', k, the(i,k), dth(i,k), sigmaw(i)
    12071369          th1(i, k) = the(i, k) - sigmaw(i)*dth(i, k) ! undisturbed area
    12081370          th2(i, k) = the(i, k) + (1.-sigmaw(i))*dth(i, k) ! wake
     
    13551517            (ph(i,kupper(i))-ph(i,1))
    13561518          crep(i, k) = crep(i, k) + crep_upper*(ph(i,1)-ph(i,k))/ &
    1357             (p(i,1)-ph(i,kupper(i)))
     1519            (ph(i,1)-ph(i,kupper(i)))
    13581520
    13591521
     
    15191681        sigmaw(i) = sigmaw(i) + d_sigmaw(i)
    15201682        d_sigmaw2(i) = d_sigmaw2(i) + d_sigmaw(i)
    1521 !  print *,'XXXX4 d_sigmaw2(i), sigmaw(i) ', d_sigmaw2(i), sigmaw(i)
    15221683      END IF
    15231684    END DO
     
    15411702          d_sig_bnd2(i) = d_sig_bnd2(i) + sigmaw_targ - sigmaw(i)
    15421703          d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i)
    1543 !  print *,'XXXX5 d_sigmaw2(i), sigmaw(i) ', d_sigmaw2(i), sigmaw(i)
    15441704          sigmaw(i) = sigmaw_targ
    15451705        END IF
     
    15781738        IF (wk_adv(i)) THEN
    15791739          wdens_targ = min( max(awdens(i),0.), wdens(i) )
     1740          d_adens_bnd2(i) = d_adens_bnd2(i) + wdens_targ - awdens(i)
    15801741          d_awdens2(i) = d_awdens2(i) + wdens_targ - awdens(i)
    15811742          awdens(i) = wdens_targ
     
    15831744      END DO
    15841745!
    1585       IF (iflag_wk_pop_dyn == 2) THEN
    1586 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! awdens again for iflag_wk_pop_dyn = 2!!!!!!
     1746      IF (iflag_wk_pop_dyn >= 2) THEN
     1747!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! awdens again for iflag_wk_pop_dyn >= 2!!!!!!
     1748!  Cumulatives
     1749        DO i = 1, klon
     1750           IF (wk_adv(i)) THEN
     1751               d_adens_death2(i)   = d_adens_death2(i)   + d_adens_death(i)
     1752               d_adens_icol2(i)   = d_adens_icol2(i)   + d_adens_icol(i)
     1753               d_adens_acol2(i)   = d_adens_acol2(i)   + d_adens_acol(i)
     1754               d_adens_bnd2(i)   = d_adens_bnd2(i)   + d_adens_bnd(i)         
     1755           END IF
     1756        END DO
     1757!  Bounds
     1758        DO i = 1, klon
     1759           IF (wk_adv(i)) THEN
     1760               wdens_targ = min( max(awdens(i),0.), wdens(i) )
     1761               d_adens_bnd2(i) = d_adens_bnd2(i) + wdens_targ - awdens(i)
     1762               awdens(i) = wdens_targ
     1763           END IF
     1764        END DO
     1765!
     1766        IF (iflag_wk_pop_dyn == 3) THEN
     1767!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! asigmaw for iflag_wk_pop_dyn = 3!!!!!!
    15871768!  Cumulatives
    15881769          DO i = 1, klon
    15891770             IF (wk_adv(i)) THEN
    1590                  d_adens_death2(i)   = d_adens_death2(i)   + d_adens_death(i)
    1591                  d_adens_icol2(i)   = d_adens_icol2(i)   + d_adens_icol(i)
    1592                  d_adens_acol2(i)   = d_adens_acol2(i)   + d_adens_acol(i)
    1593                  d_adens_bnd2(i)   = d_adens_bnd2(i)   + d_adens_bnd(i)         
     1771                 asigmaw(i) = asigmaw(i) + d_asigmaw(i)
     1772                 d_asigmaw2(i) = d_asigmaw2(i) + d_asigmaw(i)
     1773                 d_asig_death2(i)   = d_asig_death2(i)   + d_asig_death(i)
     1774                 d_asig_spread2(i)  = d_asig_spread2(i)  + d_asig_spread(i)
     1775                 d_asig_iicol2(i)   = d_asig_iicol2(i)   + d_asig_iicol(i)
     1776                 d_asig_aicol2(i)   = d_asig_aicol2(i)   + d_asig_aicol(i)
     1777                 d_asig_bnd2(i)     = d_asig_bnd2(i)     + d_asig_bnd(i)         
    15941778             END IF
    15951779          END DO
     
    15971781          DO i = 1, klon
    15981782             IF (wk_adv(i)) THEN
    1599                  wdens_targ = min( max(awdens(i),0.), wdens(i) )
    1600                  d_adens_bnd2(i) = d_adens_bnd2(i) + wdens_targ - awdens(i)
     1783   !   asigmaw lower bound set to sigmad/2 in order to allow asigmaw values lower than sigmad.
     1784   !!             sigmaw_targ = min(max(asigmaw(i),sigmad),sigmaw(i))
     1785                sigmaw_targ = min(max(asigmaw(i),sigmad/2.),sigmaw(i))
     1786                d_asig_bnd2(i) = d_asig_bnd2(i) + sigmaw_targ - asigmaw(i)
     1787                d_asigmaw2(i) = d_asigmaw2(i) + sigmaw_targ - asigmaw(i)
     1788                asigmaw(i) = sigmaw_targ
    16011789             END IF
    16021790          END DO
    1603       ENDIF ! (iflag_wk_pop_dyn == 2)
     1791
     1792#ifdef IOPHYS_WK
     1793    IF (phys_sub) THEN
     1794     CALL iophys_ecrit('wdensb',1,'wdensb','m',wdens)
     1795     CALL iophys_ecrit('awdensb',1,'awdensb','m',awdens)
     1796     CALL iophys_ecrit('sigmawb',1,'sigmawb','m',sigmaw)
     1797     CALL iophys_ecrit('asigmawb',1,'asigmawb','m',asigmaw)
     1798!
     1799     call iophys_ecrit('d_wdens2',1,'d_wdens2','',d_wdens2)
     1800     call iophys_ecrit('d_dens_gen2',1,'d_dens_gen2','',d_dens_gen2)
     1801     call iophys_ecrit('d_dens_death2',1,'d_dens_death2','',d_dens_death2)
     1802     call iophys_ecrit('d_dens_col2',1,'d_dens_col2','',d_dens_col2)
     1803     call iophys_ecrit('d_dens_bnd2',1,'d_dens_bnd2','',d_dens_bnd2)
     1804!
     1805     call iophys_ecrit('d_awdens2',1,'d_awdens2','',d_awdens2)
     1806     call iophys_ecrit('d_adens_death2',1,'d_adens_death2','',d_adens_death2)
     1807     call iophys_ecrit('d_adens_icol2',1,'d_adens_icol2','',d_adens_icol2)
     1808     call iophys_ecrit('d_adens_acol2',1,'d_adens_acol2','',d_adens_acol2)
     1809     call iophys_ecrit('d_adens_bnd2',1,'d_adens_bnd2','',d_adens_bnd2)
     1810!
     1811     CALL iophys_ecrit('d_sigmaw2',1,'d_sigmaw2','',d_sigmaw2)
     1812     CALL iophys_ecrit('d_sig_gen2',1,'d_sig_gen2','m',d_sig_gen2)
     1813     CALL iophys_ecrit('d_sig_spread2',1,'d_sig_spread2','',d_sig_spread2)
     1814     CALL iophys_ecrit('d_sig_col2',1,'d_sig_col2','',d_sig_col2)
     1815     CALL iophys_ecrit('d_sig_death2',1,'d_sig_death2','',d_sig_death2)
     1816     CALL iophys_ecrit('d_sig_bnd2',1,'d_sig_bnd2','',d_sig_bnd2)
     1817!
     1818     CALL iophys_ecrit('d_asigmaw2',1,'d_asigmaw2','',d_asigmaw2)
     1819     CALL iophys_ecrit('d_asig_spread2',1,'d_asig_spread2','m',d_asig_spread2)
     1820     CALL iophys_ecrit('d_asig_aicol2',1,'d_asig_aicol2','m',d_asig_aicol2)
     1821     CALL iophys_ecrit('d_asig_iicol2',1,'d_asig_iicol2','m',d_asig_iicol2)
     1822     CALL iophys_ecrit('d_asig_death2',1,'d_asig_death2','m',d_asig_death2)
     1823     CALL iophys_ecrit('d_asig_bnd2',1,'d_asig_bnd2','m',d_asig_bnd2)
     1824    ENDIF
     1825#endif
     1826        ENDIF ! (iflag_wk_pop_dyn == 3)
     1827      ENDIF ! (iflag_wk_pop_dyn >= 2)
    16041828    ENDIF  ! (iflag_wk_pop_dyn >= 1)
    16051829
     
    18182042    END DO
    18192043
     2044
    18202045    ! Filter out bad wakes
    18212046
     
    18452070          d_sig_bnd2(i) = d_sig_bnd2(i) + sigmaw_targ - sigmaw(i)
    18462071          d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i)
    1847 !  print *,'XXXX6 d_sigmaw2(i), sigmaw(i) ', d_sigmaw2(i), sigmaw(i)
    18482072          sigmaw(i) = sigmaw_targ
     2073!
     2074          d_asig_bnd2(i) = d_asig_bnd2(i) + sigmaw_targ - asigmaw(i)
     2075          d_asigmaw2(i) = d_asigmaw2(i) + sigmaw_targ - asigmaw(i)
     2076          asigmaw(i) = sigmaw_targ
    18492077!>jyg
    18502078          fip(i) = 0.
     
    18562084      END IF
    18572085    END DO
    1858 
    1859   END DO ! end sub-timestep loop
    1860 
     2086  !
     2087  ! ------------------------------------------------------------------------
     2088  !
     2089  END DO                  ! end sub-timestep loop
     2090  !
     2091  ! ------------------------------------------------------------------------
     2092  ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
     2093  ! ------------------------------------------------------------------------
     2094  !
     2095
     2096#ifdef IOPHYS_WK
     2097    IF (.not.phys_sub) CALL iophys_ecrit('wape_b',1,'wape_b','J/kg',wape)
     2098#endif
    18612099  IF (prt_level>=10) THEN
    18622100    PRINT *, 'wake-5, sigmaw(igout), cstar(igout), wape(igout), ptop(igout) ', &
     
    20042242    END IF
    20052243  END DO
    2006 
     2244#ifdef IOPHYS_WK
     2245  IF (.not.phys_sub) CALL iophys_ecrit('wape2_a',1,'wape2_a','J/kg',wape2)
     2246#endif
    20072247
    20082248
     
    20342274    END DO
    20352275  END IF
     2276#ifdef IOPHYS_WK
     2277  IF (.not.phys_sub) CALL iophys_ecrit('wape2_b',1,'wape2_b','J/kg',wape2)
     2278#endif
    20362279
    20372280
     
    20642307      d_sig_bnd2(i) = d_sig_bnd2(i) + sigmaw_targ - sigmaw(i)
    20652308      d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i)
    2066 !  print *,'XXXX7 d_sigmaw2(i), sigmaw(i) ', d_sigmaw2(i), sigmaw(i)
    20672309      sigmaw(i) = sigmaw_targ
     2310!
     2311      d_asig_bnd2(i) = d_asig_bnd2(i) + sigmaw_targ - asigmaw(i)
     2312      d_asigmaw2(i) = d_asigmaw2(i) + sigmaw_targ - asigmaw(i)
     2313      asigmaw(i) = sigmaw_targ
    20682314!>jyg
    20692315        fip(i) = 0.
     
    20742320        gwake(i) = .TRUE.
    20752321      END IF
    2076     END IF
     2322#ifdef IOPHYS_WK
     2323  IF (.not.phys_sub) CALL iophys_ecrit('cstar2',1,'cstar2','J/kg',cstar2)
     2324#endif
     2325    END IF  ! (ok_qx_qw(i))
    20772326  END DO
    20782327
     
    21062355    END IF
    21072356  END DO
    2108 
     2357    IF (iflag_wk_pop_dyn >= 3) THEN
     2358#ifdef IOPHYS_WK
     2359      IF (.not.phys_sub) THEN
     2360       CALL iophys_ecrit('fip',1,'fip','J/kg',fip)
     2361       CALL iophys_ecrit('hw',1,'hw','J/kg',hw)
     2362       CALL iophys_ecrit('ptop',1,'ptop','J/kg',ptop)
     2363       CALL iophys_ecrit('wdens',1,'wdens','J/kg',wdens)
     2364       CALL iophys_ecrit('awdens',1,'awdens','m',awdens)
     2365       CALL iophys_ecrit('sigmaw',1,'sigmaw','m',sigmaw)
     2366       CALL iophys_ecrit('asigmaw',1,'asigmaw','m',asigmaw)
     2367!   
     2368       CALL iophys_ecrit('rad_wk',1,'rad_wk','J/kg',rad_wk)
     2369       CALL iophys_ecrit('arad_wk',1,'arad_wk','J/kg',arad_wk)
     2370       CALL iophys_ecrit('irad_wk',1,'irad_wk','J/kg',irad_wk)
     2371!   
     2372       call iophys_ecrit('d_wdens2',1,'d_wdens2','',d_wdens2)
     2373       call iophys_ecrit('d_dens_gen2',1,'d_dens_gen2','',d_dens_gen2)
     2374       call iophys_ecrit('d_dens_death2',1,'d_dens_death2','',d_dens_death2)
     2375       call iophys_ecrit('d_dens_col2',1,'d_dens_col2','',d_dens_col2)
     2376       call iophys_ecrit('d_dens_bnd2',1,'d_dens_bnd2','',d_dens_bnd2)
     2377!   
     2378       call iophys_ecrit('d_awdens2',1,'d_awdens2','',d_awdens2)
     2379       call iophys_ecrit('d_adens_death2',1,'d_adens_death2','',d_adens_death2)
     2380       call iophys_ecrit('d_adens_icol2',1,'d_adens_icol2','',d_adens_icol2)
     2381       call iophys_ecrit('d_adens_acol2',1,'d_adens_acol2','',d_adens_acol2)
     2382       call iophys_ecrit('d_adens_bnd2',1,'d_adens_bnd2','',d_adens_bnd2)
     2383!   
     2384       CALL iophys_ecrit('d_sigmaw2',1,'d_sigmaw2','',d_sigmaw2)
     2385       CALL iophys_ecrit('d_sig_gen2',1,'d_sig_gen2','m',d_sig_gen2)
     2386       CALL iophys_ecrit('d_sig_spread2',1,'d_sig_spread2','',d_sig_spread2)
     2387       CALL iophys_ecrit('d_sig_col2',1,'d_sig_col2','',d_sig_col2)
     2388       CALL iophys_ecrit('d_sig_death2',1,'d_sig_death2','',d_sig_death2)
     2389       CALL iophys_ecrit('d_sig_bnd2',1,'d_sig_bnd2','',d_sig_bnd2)
     2390!   
     2391       CALL iophys_ecrit('d_asigmaw2',1,'d_asigmaw2','',d_asigmaw2)
     2392       CALL iophys_ecrit('d_asig_spread2',1,'d_asig_spread2','m',d_asig_spread2)
     2393       CALL iophys_ecrit('d_asig_aicol2',1,'d_asig_aicol2','m',d_asig_aicol2)
     2394       CALL iophys_ecrit('d_asig_iicol2',1,'d_asig_iicol2','m',d_asig_iicol2)
     2395       CALL iophys_ecrit('d_asig_death2',1,'d_asig_death2','m',d_asig_death2)
     2396       CALL iophys_ecrit('d_asig_bnd2',1,'d_asig_bnd2','m',d_asig_bnd2)
     2397      ENDIF  ! (.not.phys_sub)
     2398#endif
     2399    ENDIF  ! (iflag_wk_pop_dyn >= 3)
    21092400  ! Limitation de sigmaw
    21102401
     
    21212412    DO i = 1, klon
    21222413      kill_wake(i) = ((wape(i)>=wape2(i)) .AND. (wape2(i)<=wapecut)) .OR. (ktopw(i)<=2) .OR. &
    2123           .NOT. ok_qx_qw(i) .OR. (wdens(i) < 2.*wdensmin)
     2414          .NOT. ok_qx_qw(i) .OR. (wdens(i) < wdensthreshold)
     2415!!          .NOT. ok_qx_qw(i) .OR. (wdens(i) < 2.*wdensmin)
    21242416    ENDDO
    21252417  ELSE  ! (iflag_wk_pop_dyn >= 1)
     
    21542446      wape(i) = 0.
    21552447      cstar(i) = 0.
    2156 !!jyg   Outside subroutine "Wake" hw, wdens and sigmaw are zero when there are no wakes
     2448!!jyg   Outside subroutine "Wake" hw, wdens sigmaw and asigmaw are zero when there are no wakes
    21572449!!      hw(i) = hwmin                       !jyg
    21582450!!      sigmaw(i) = sigmad                  !jyg
    21592451      hw(i) = 0.                            !jyg
    21602452      fip(i) = 0.
     2453!
    21612454!!      sigmaw(i) = 0.                        !jyg
    21622455      sigmaw_targ = 0.
     
    21642457!!      d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i)
    21652458      d_sigmaw2(i) = sigmaw_targ - sigmaw_in(i)      ! _in = correction jyg 20220124
    2166 !  print *,'XXXX8 d_sigmaw2(i), sigmaw(i) ', d_sigmaw2(i), sigmaw(i)
    21672459      sigmaw(i) = sigmaw_targ
     2460!
     2461      IF (iflag_wk_pop_dyn >= 3) THEN
     2462        sigmaw_targ = 0.
     2463        d_asig_bnd2(i) = d_asig_bnd2(i) + sigmaw_targ - asigmaw(i)
     2464!!        d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i)
     2465        d_asigmaw2(i) = sigmaw_targ - asigmaw_in(i)      ! _in = correction jyg 20220124
     2466        asigmaw(i) = sigmaw_targ
     2467      ELSE
     2468        asigmaw(i) = 0.
     2469      ENDIF ! (iflag_wk_pop_dyn >= 3)
     2470!
    21682471      IF (iflag_wk_pop_dyn >= 1) THEN
    21692472!!        awdens(i) = 0.
     
    21762479        wdens_targ = 0.
    21772480!!jyg: bug fix : the d_adens_bnd2 computation must be before the update of awdens.
    2178         IF (iflag_wk_pop_dyn == 2) THEN
     2481        IF (iflag_wk_pop_dyn >= 2) THEN
    21792482            d_adens_bnd2(i) = d_adens_bnd2(i) + wdens_targ - awdens(i)
    2180         ENDIF ! (iflag_wk_pop_dyn == 2)
     2483        ENDIF ! (iflag_wk_pop_dyn >= 2)
    21812484!!        d_awdens2(i) = wdens_targ - awdens(i)
    21822485        d_awdens2(i) = wdens_targ - awdens_in(i)    ! jyg 20220916
     
    21982501                      wape(igout),wape2(igout),ktopw(igout),OK_qx_qw(igout)
    21992502  ENDIF
     2503#ifdef IOPHYS_WK
     2504  IF (.not.phys_sub) CALL iophys_ecrit('wape_c',1,'wape_c','J/kg',wape)
     2505#endif
    22002506
    22012507
     
    22242530!jyg<
    22252531  IF (iflag_wk_pop_dyn >= 1) THEN
    2226   DO i = 1, klon
    2227       IF (ok_qx_qw(i)) THEN
    2228     d_sig_gen2(i) = d_sig_gen2(i)/dtime
    2229     d_sig_death2(i) = d_sig_death2(i)/dtime
    2230     d_sig_col2(i) = d_sig_col2(i)/dtime
    2231     d_sig_spread2(i) = d_sig_spread2(i)/dtime
    2232     d_sig_bnd2(i) = d_sig_bnd2(i)/dtime
    2233     d_sigmaw2(i) = d_sigmaw2(i)/dtime
    2234 !  print *,'XXXX9 d_sigmaw2(i), sigmaw(i), dtime ', d_sigmaw2(i), sigmaw(i), dtime
    2235 !
    2236     d_dens_gen2(i) = d_dens_gen2(i)/dtime
    2237     d_dens_death2(i) = d_dens_death2(i)/dtime
    2238     d_dens_col2(i) = d_dens_col2(i)/dtime
    2239     d_dens_bnd2(i) = d_dens_bnd2(i)/dtime
    2240     d_awdens2(i) = d_awdens2(i)/dtime
    2241     d_wdens2(i) = d_wdens2(i)/dtime
    2242       ENDIF
    2243   ENDDO
    2244   IF (iflag_wk_pop_dyn == 2) THEN
    2245     DO i = 1, klon
    2246       IF (ok_qx_qw(i)) THEN
    2247     d_adens_death2(i) = d_adens_death2(i)/dtime
    2248     d_adens_icol2(i) = d_adens_icol2(i)/dtime
    2249     d_adens_acol2(i) = d_adens_acol2(i)/dtime
    2250     d_adens_bnd2(i) = d_adens_bnd2(i)/dtime
    2251       ENDIF
     2532    DO i = 1, klon
     2533        IF (ok_qx_qw(i)) THEN
     2534      d_sig_gen2(i) = d_sig_gen2(i)/dtime
     2535      d_sig_death2(i) = d_sig_death2(i)/dtime
     2536      d_sig_col2(i) = d_sig_col2(i)/dtime
     2537      d_sig_spread2(i) = d_sig_spread2(i)/dtime
     2538      d_sig_bnd2(i) = d_sig_bnd2(i)/dtime
     2539      d_sigmaw2(i) = d_sigmaw2(i)/dtime
     2540!
     2541      d_dens_gen2(i) = d_dens_gen2(i)/dtime
     2542      d_dens_death2(i) = d_dens_death2(i)/dtime
     2543      d_dens_col2(i) = d_dens_col2(i)/dtime
     2544      d_dens_bnd2(i) = d_dens_bnd2(i)/dtime
     2545      d_awdens2(i) = d_awdens2(i)/dtime
     2546      d_wdens2(i) = d_wdens2(i)/dtime
     2547        ENDIF
    22522548    ENDDO
    2253    ENDIF ! (iflag_wk_pop_dyn == 2) 
     2549    IF (iflag_wk_pop_dyn >= 2) THEN
     2550      DO i = 1, klon
     2551        IF (ok_qx_qw(i)) THEN
     2552        d_adens_death2(i) = d_adens_death2(i)/dtime
     2553        d_adens_icol2(i) = d_adens_icol2(i)/dtime
     2554        d_adens_acol2(i) = d_adens_acol2(i)/dtime
     2555        d_adens_bnd2(i) = d_adens_bnd2(i)/dtime
     2556        ENDIF
     2557      ENDDO
     2558      IF (iflag_wk_pop_dyn == 3) THEN
     2559       DO i = 1, klon
     2560          IF (ok_qx_qw(i)) THEN
     2561        d_asig_death2(i)  = d_asig_death2(i)/dtime
     2562        d_asig_iicol2(i)  = d_asig_iicol2(i)/dtime
     2563        d_asig_aicol2(i)  = d_asig_aicol2(i)/dtime
     2564        d_asig_spread2(i) = d_asig_spread2(i)/dtime
     2565        d_asig_bnd2(i) = d_asig_bnd2(i)/dtime
     2566          ENDIF
     2567       ENDDO
     2568      ENDIF ! (iflag_wk_pop_dyn == 3) 
     2569    ENDIF ! (iflag_wk_pop_dyn >= 2) 
    22542570  ENDIF  ! (iflag_wk_pop_dyn >= 1)
    22552571 
     
    23872703
    23882704SUBROUTINE wake_popdyn_1(klon, klev, dtime, cstar, tau_wk_inv, wgen, wdens, awdens, sigmaw, &
     2705                  wdensmin, &
    23892706                  dtimesub, gfl, rad_wk, f_shear, drdt_pos, &
    23902707                  d_awdens, d_wdens, d_sigmaw, &
     
    24002717  USE lmdz_wake_ini , ONLY : stark, wdens_ref
    24012718  USE lmdz_wake_ini , ONLY : tau_cv, rzero, aa0
    2402   USE lmdz_wake_ini , ONLY : iflag_wk_pop_dyn, wdensmin
     2719!!  USE lmdz_wake_ini , ONLY : iflag_wk_pop_dyn, wdensmin
     2720  USE lmdz_wake_ini , ONLY : iflag_wk_pop_dyn
    24032721  USE lmdz_wake_ini , ONLY : sigmad, cstart, sigmaw_max
    24042722 
     
    24092727  REAL,                             INTENT(IN)          :: dtime
    24102728  REAL,                             INTENT(IN)          :: dtimesub
     2729  REAL,                             INTENT(IN)          :: wdensmin
    24112730  REAL, DIMENSION (klon),           INTENT(IN)          :: wgen
    24122731  REAL, DIMENSION (klon),           INTENT(IN)          :: wdens
    24132732  REAL, DIMENSION (klon),           INTENT(IN)          :: awdens
    24142733  REAL, DIMENSION (klon),           INTENT(IN)          :: sigmaw
    2415   REAL, DIMENSION (klon),           INTENT(IN)          :: gfl, cstar
     2734  REAL, DIMENSION (klon),           INTENT(IN)          :: cstar
    24162735  REAL, DIMENSION (klon),           INTENT(IN)          :: cin, wape
    2417   REAL, DIMENSION (klon),           INTENT(IN)          :: rad_wk
    24182736  REAL, DIMENSION (klon),           INTENT(IN)          :: f_shear
    24192737  INTEGER,                          INTENT(IN)          :: iflag_wk_act
     
    24242742  ! Tendencies of state variables (2 is appended to the names of fields which are the cumul of fields
    24252743  !                                 computed at each sub-timestep; e.g. d_wdens2 is the cumul of d_wdens)
     2744  REAL, DIMENSION (klon),           INTENT(OUT)         :: rad_wk
     2745  REAL, DIMENSION (klon),           INTENT(OUT)         :: gfl
    24262746  REAL, DIMENSION (klon),           INTENT(OUT)         :: d_sigmaw, d_awdens, d_wdens
    24272747  REAL, DIMENSION (klon),           INTENT(OUT)         :: drdt
     
    24642784      ENDIF  ! (iflag_wk_act ==2)
    24652785
    2466 
    2467       DO i = 1, klon
    2468        ! print*, 'XXX wk_adv(i)', wk_adv(i)
     2786      DO i = 1, klon
     2787        IF (wk_adv(i)) THEN
     2788          rad_wk(i) = max( sqrt(sigmaw(i)/(3.14*wdens(i))) , rzero)
     2789          gfl(i)  = 2.*sqrt(3.14*wdens(i)*sigmaw(i))
     2790        END IF
     2791      END DO
     2792
     2793      DO i = 1, klon
    24692794        IF (wk_adv(i)) THEN
    24702795!!          tau_wk(i) = max(rad_wk(i)/(3.*cstar(i))*((cstar(i)/cstart)**1.5 - 1), 100.)
     
    25062831!!                      - sigmaw(i)*tau_wk_inv_min )*dtimesub
    25072832          d_sig_gen(i) = wgen(i)*aa0
    2508 !!          print*, 'XXX sigmaw(i), awdens(i), wdens(i), tau_wk_inv_min', &
    2509 !!                  sigmaw(i), awdens(i), wdens(i), tau_wk_inv_min
    25102833          d_sig_death(i) = - sigmaw(i)*(1.-awdens(i)/wdens(i))*tau_wk_inv_min
    25112834!!       
     
    25472870   
    25482871    SUBROUTINE wake_popdyn_2 ( klon, klev, wk_adv, dtimesub, wgen, &
     2872                             wdensmin, &
    25492873                             sigmaw, wdens, awdens, &   !! states variables
    25502874                             gfl, cstar, cin, wape, rad_wk, &
     
    25612885  USE lmdz_wake_ini , ONLY : stark, wdens_ref
    25622886  USE lmdz_wake_ini , ONLY : tau_cv, rzero, aa0
    2563   USE lmdz_wake_ini , ONLY : iflag_wk_pop_dyn, wdensmin
     2887!!  USE lmdz_wake_ini , ONLY : iflag_wk_pop_dyn, wdensmin
     2888  USE lmdz_wake_ini , ONLY : iflag_wk_pop_dyn
    25642889  USE lmdz_wake_ini , ONLY : sigmad, cstart, sigmaw_max
    25652890 
     
    25692894  LOGICAL, DIMENSION (klon),        INTENT(IN)          :: wk_adv
    25702895  REAL,                             INTENT(IN)          :: dtimesub
     2896  REAL,                             INTENT(IN)          :: wdensmin
    25712897  REAL, DIMENSION (klon),           INTENT(IN)          :: wgen      !! B = birth rate of wakes
    25722898  REAL, DIMENSION (klon),           INTENT(INOUT)       :: sigmaw    !! sigma = fractional area of wakes
    25732899  REAL, DIMENSION (klon),           INTENT(INOUT)       :: wdens     !! D = number of wakes per unit area
    25742900  REAL, DIMENSION (klon),           INTENT(INOUT)       :: awdens    !! A = number of active wakes per unit area
    2575   REAL, DIMENSION (klon),           INTENT(IN)          :: gfl       !! Lg = gust front lenght per unit area
    25762901  REAL, DIMENSION (klon),           INTENT(IN)          :: cstar     !! C* = spreading velocity of wakes
    25772902  REAL, DIMENSION (klon),           INTENT(IN)          :: cin, wape  ! RM : A Faire disparaitre
    2578   REAL, DIMENSION (klon),           INTENT(IN)          :: rad_wk    !! r = wake radius
    25792903
    25802904
     2905  REAL, DIMENSION (klon),           INTENT(OUT)         :: rad_wk    !! r = wake radius
     2906  REAL, DIMENSION (klon),           INTENT(OUT)         :: gfl       !! Lg = gust front lenght per unit area
    25812907  REAL, DIMENSION (klon),           INTENT(OUT)         :: d_sigmaw, d_wdens, d_awdens
    25822908  REAL, DIMENSION (klon),           INTENT(OUT)         :: cont_fact  !! RM facteur de contact = 2 pi * rad * C*
     
    26052931
    26062932      DO i = 1, klon
    2607        ! print*, 'XXX wk_adv(i)', wk_adv(i)
     2933        IF (wk_adv(i)) THEN
     2934          rad_wk(i) = max( sqrt(sigmaw(i)/(3.14*wdens(i))) , rzero)
     2935          gfl(i)  = 2.*sqrt(3.14*wdens(i)*sigmaw(i))
     2936        END IF
     2937      END DO
     2938
     2939
     2940      DO i = 1, klon
    26082941        IF (wk_adv(i)) THEN
    26092942!!          tau_wk(i) = max(rad_wk(i)/(3.*cstar(i))*((cstar(i)/cstart)**1.5 - 1), 100.)
     
    26873020    END SUBROUTINE wake_popdyn_2 
    26883021 
     3022    SUBROUTINE wake_popdyn_3 ( klon, klev, phys_sub, wk_adv, dtimesub, wgen, &
     3023                             wdensmin, &
     3024                             sigmaw, asigmaw, wdens, awdens, &                       !! state variables
     3025                             gfl, agfl, cstar, cin, wape, &
     3026                             rad_wk, arad_wk, irad_wk, &
     3027                             d_sigmaw, d_asigmaw, d_wdens, d_awdens, &               !! tendencies
     3028                             d_sig_gen, d_sig_death, d_sig_col, d_sig_spread, d_sig_bnd, &
     3029                             d_asig_death, d_asig_aicol, d_asig_iicol, d_asig_spread, d_asig_bnd, &
     3030                             d_dens_gen, d_dens_death, d_dens_col, d_dens_bnd, &
     3031                             d_adens_death, d_adens_icol, d_adens_acol, d_adens_bnd )
     3032                             
     3033                                             
     3034
     3035  USE lmdz_wake_ini , ONLY : wake_ini
     3036  USE lmdz_wake_ini , ONLY : prt_level,RG
     3037  USE lmdz_wake_ini , ONLY : stark, wdens_ref
     3038  USE lmdz_wake_ini , ONLY : tau_cv, rzero, aa0
     3039!!  USE lmdz_wake_ini , ONLY : iflag_wk_pop_dyn, wdensmin
     3040  USE lmdz_wake_ini , ONLY : iflag_wk_pop_dyn
     3041  USE lmdz_wake_ini , ONLY : sigmad, cstart, sigmaw_max
     3042  USE lmdz_wake_ini , ONLY : smallestreal
     3043 
     3044IMPLICIT NONE
     3045
     3046  INTEGER, INTENT(IN)                                   :: klon,klev
     3047  LOGICAL,                          INTENT(IN)          :: phys_sub
     3048  LOGICAL, DIMENSION (klon),        INTENT(IN)          :: wk_adv
     3049  REAL,                             INTENT(IN)          :: dtimesub
     3050  REAL,                             INTENT(IN)          :: wdensmin
     3051  REAL, DIMENSION (klon),           INTENT(IN)          :: wgen      !! B = birth rate of wakes
     3052  REAL, DIMENSION (klon),           INTENT(INOUT)       :: sigmaw    !! sigma = fractional area of wakes
     3053  REAL, DIMENSION (klon),           INTENT(INOUT)       :: asigmaw   !! sigma = fractional area of active wakes
     3054  REAL, DIMENSION (klon),           INTENT(INOUT)       :: wdens     !! D = number of wakes per unit area
     3055  REAL, DIMENSION (klon),           INTENT(INOUT)       :: awdens    !! A = number of active wakes per unit area
     3056  REAL, DIMENSION (klon),           INTENT(IN)          :: cstar     !! C* = spreading velocity of wakes
     3057  REAL, DIMENSION (klon),           INTENT(IN)          :: cin, wape  ! RM : A Faire disparaitre
     3058
     3059
     3060  REAL, DIMENSION (klon),           INTENT(OUT)         :: rad_wk    !! r = mean wake radius
     3061  REAL, DIMENSION (klon),           INTENT(OUT)         :: arad_wk    !! r_A = wake radius of active wakes
     3062  REAL, DIMENSION (klon),           INTENT(OUT)         :: irad_wk    !! r_I = wake radius of inactive wakes
     3063  REAL, DIMENSION (klon),           INTENT(OUT)         :: gfl       !! Lg = gust front length per unit area
     3064  REAL, DIMENSION (klon),           INTENT(OUT)         :: agfl      !! LgA = gust front length of active wakes
     3065                                                                     !!  per unit area
     3066  REAL, DIMENSION (klon),           INTENT(OUT)         :: d_sigmaw, d_asigmaw, d_wdens, d_awdens
     3067  ! Some components of the tendencies of state variables 
     3068  REAL, DIMENSION (klon),           INTENT(OUT)         :: d_sig_gen, d_sig_death, d_sig_col, d_sig_spread, d_sig_bnd
     3069  REAL, DIMENSION (klon),           INTENT(OUT)         :: d_asig_death, d_asig_aicol, d_asig_iicol, d_asig_spread, d_asig_bnd
     3070  REAL, DIMENSION (klon),           INTENT(OUT)         :: d_dens_gen, d_dens_death, d_dens_col, d_dens_bnd
     3071  REAL, DIMENSION (klon),           INTENT(OUT)         :: d_adens_death, d_adens_acol, d_adens_icol, d_adens_bnd
     3072
     3073
     3074!! internal variables
     3075 
     3076  INTEGER                                               :: i, k
     3077  REAL, DIMENSION (klon)                                :: iwdens, isigmaw !! inactive wake density and fractional area
     3078!!  REAL, DIMENSION (klon)                                :: d_arad, d_irad
     3079  REAL, DIMENSION (klon)                                :: igfl            !! LgI = gust front length of inactive wakes
     3080                                                                           !!  per unit area
     3081  REAL, DIMENSION (klon)                                :: s_wk            !! mean area of individual wakes
     3082  REAL, DIMENSION (klon)                                :: as_wk           !! mean area of individual active wakes
     3083  REAL, DIMENSION (klon)                                :: is_wk           !! mean area of individual inactive wakes
     3084  REAL, DIMENSION (klon)                                :: tau_wk_inv      !! tau = life time of wakes
     3085  REAL                                                  :: tau_wk_inv_min
     3086  REAL, DIMENSION (klon)                                :: tau_prime       !! tau_prime = life time of actives wakes
     3087  REAL                                                  :: d_wdens_targ, d_sigmaw_targ
     3088
     3089
     3090!! Equations
     3091!! ---------
     3092!! Gust fronts:
     3093!! Lg_A = 2 pi r_A A
     3094!! Lg_I = 2 pi r_I I
     3095!! Lg   = 2 pi r   D
     3096!!
     3097!! Areas:
     3098!! s = pi r^2
     3099!! s_A = pi r_A^2
     3100!! s_I = pi r_I^2
     3101!!
     3102!! Life expectancy:
     3103!! tau_I = 3 C* ((C*/C*t)^3/2 - 1) / r_I
     3104!!
     3105!! Time deratives:
     3106!! dD/dt = B - (D-A)/tau_I - 2 Lg C* D
     3107!! dA/dt = B - A/tau_A     + 2 Lg_I C* (D-A) - 2 Lg_A C* A
     3108!! dsigma/dt = B a0 - sigma_I/tau_I + Lg C* - 2 Lg_I C* (D-A) (2 s_I - a0)
     3109!! dsigma_A/dt = B a0 - sigma_A/tau_A + Lg_A C* + (Lg_A I + Lg_I A) C* s_I + 2 Lg_I C* I a0
     3110!!
     3111
     3112      DO i = 1, klon
     3113        IF (wk_adv(i)) THEN
     3114         iwdens(i) = wdens(i) - awdens(i)
     3115         isigmaw(i) = sigmaw(i) - asigmaw(i)
     3116!
     3117         arad_wk(i) = max( sqrt(asigmaw(i)/(3.14*awdens(i))) , rzero)
     3118         irad_wk(i) = max( sqrt((sigmaw(i)-asigmaw(i))/  &
     3119                           (3.14*max(smallestreal,(wdens(i)-awdens(i))))), rzero)
     3120         rad_wk(i) = (awdens(i)*arad_wk(i)+(wdens(i)-awdens(i))*irad_wk(i))/wdens(i)
     3121!
     3122         s_wk(i) = 3.14*rad_wk(i)**2
     3123         as_wk(i) = 3.14*arad_wk(i)**2
     3124         is_wk(i) = 3.14*irad_wk(i)**2
     3125!
     3126         gfl(i)  = 2.*sqrt(3.14*wdens(i)*sigmaw(i))
     3127         agfl(i) = 2.*sqrt(3.14*awdens(i)*asigmaw(i))
     3128         igfl(i) = gfl(i) - agfl(i)
     3129        ENDIF
     3130      ENDDO
     3131
     3132
     3133      DO i = 1, klon
     3134        IF (wk_adv(i)) THEN
     3135          tau_wk_inv(i) = max( (3.*cstar(i))/(irad_wk(i)*((cstar(i)/cstart)**1.5 - 1)), 0.)
     3136          tau_wk_inv_min = min(tau_wk_inv(i), 1./dtimesub)
     3137          tau_prime(i) = tau_cv
     3138
     3139          d_sig_gen(i) = wgen(i)*aa0
     3140          d_sig_death(i) = - isigmaw(i)*tau_wk_inv_min
     3141          d_sig_col(i) = - 2.*igfl(i)*cstar(i)*iwdens(i)*(2.*is_wk(i)-aa0)
     3142          d_sig_spread(i) = gfl(i)*cstar(i)
     3143!
     3144          d_sig_gen(i) =  d_sig_gen(i)*dtimesub
     3145          d_sig_death(i) = d_sig_death(i)*dtimesub
     3146          d_sig_col(i) =  d_sig_col(i)*dtimesub
     3147          d_sig_spread(i) =  d_sig_spread(i)*dtimesub
     3148          d_sigmaw(i) =  d_sig_gen(i) + d_sig_death(i) + d_sig_col(i) + d_sig_spread(i)
     3149#ifdef IOPHYS_WK
     3150          IF (phys_sub) call iophys_ecrit('d_sigmaw0',1,'d_sigmaw0','',d_sigmaw)
     3151#endif
     3152
     3153         
     3154          d_sigmaw_targ = max(d_sigmaw(i), sigmad-sigmaw(i))
     3155!!          d_sig_bnd(i) = d_sig_bnd(i) + d_sigmaw_targ - d_sigmaw(i)
     3156!!          d_sig_bnd_provis(i) = d_sigmaw_targ - d_sigmaw(i)
     3157          d_sig_bnd(i) = d_sigmaw_targ - d_sigmaw(i)
     3158          d_sigmaw(i) = d_sigmaw_targ
     3159!!          d_sigmaw(i) = max(d_sigmaw(i), sigmad-sigmaw(i))
     3160#ifdef IOPHYS_WK
     3161          IF (phys_sub) THEN
     3162             call iophys_ecrit('tauwk_inv',1,'tau_wk_inv_min','',tau_wk_inv_min)
     3163             call iophys_ecrit('d_sigmaw',1,'d_sigmaw','',d_sigmaw)
     3164             call iophys_ecrit('d_sig_gen',1,'d_sig_gen','',d_sig_gen)
     3165             call iophys_ecrit('d_sig_death',1,'d_sig_death','',d_sig_death)
     3166             call iophys_ecrit('d_sig_col',1,'d_sig_col','',d_sig_col)
     3167             call iophys_ecrit('d_sig_spread',1,'d_sig_spread','',d_sig_spread)
     3168             call iophys_ecrit('d_sig_bnd',1,'d_sig_bnd','',d_sig_bnd)
     3169          ENDIF
     3170#endif
     3171          d_asig_death(i) = - asigmaw(i)/tau_prime(i)
     3172          d_asig_aicol(i) = (agfl(i)*iwdens(i) + igfl(i)*awdens(i))*cstar(i)*is_wk(i)
     3173          d_asig_iicol(i) = 2.*igfl(i)*cstar(i)*iwdens(i)*aa0
     3174          d_asig_spread(i) = agfl(i)*cstar(i)
     3175!
     3176          d_asig_death(i) = d_asig_death(i)*dtimesub
     3177          d_asig_aicol(i) =  d_asig_aicol(i)*dtimesub
     3178          d_asig_iicol(i) =  d_asig_iicol(i)*dtimesub
     3179          d_asig_spread(i) =  d_asig_spread(i)*dtimesub
     3180          d_asigmaw(i) =  d_sig_gen(i) + d_asig_death(i) + d_asig_aicol(i) + d_asig_iicol(i) + d_asig_spread(i)
     3181#ifdef IOPHYS_WK
     3182          IF (phys_sub) call iophys_ecrit('d_asigmaw0',1,'d_asigmaw0','',d_asigmaw)
     3183#endif
     3184
     3185          d_sigmaw_targ = min(max(d_asigmaw(i),-asigmaw(i)), sigmaw(i)-asigmaw(i))
     3186!!          d_dens_bnd(i) = d_dens_bnd(i) + d_sigmaw_targ - d_sigmaw(i)
     3187          d_asig_bnd(i) = d_sigmaw_targ - d_asigmaw(i)
     3188          d_asigmaw(i) = d_sigmaw_targ
     3189#ifdef IOPHYS_WK
     3190          IF (phys_sub) THEN
     3191             call iophys_ecrit('d_asigmaw',1,'d_asigmaw','',d_asigmaw)
     3192             call iophys_ecrit('d_asig_death',1,'d_asig_death','',d_asig_death)
     3193             call iophys_ecrit('d_asig_aicol',1,'d_asig_aicol','',d_asig_aicol)
     3194             call iophys_ecrit('d_asig_iicol',1,'d_asig_iicol','',d_asig_iicol)
     3195             call iophys_ecrit('d_asig_spread',1,'d_asig_spread','',d_asig_spread)
     3196             call iophys_ecrit('d_asig_bnd',1,'d_asig_bnd','',d_asig_bnd)
     3197          ENDIF
     3198#endif
     3199          d_dens_gen(i) = wgen(i)
     3200          d_dens_death(i) = - iwdens(i)*tau_wk_inv_min
     3201          d_dens_col(i) =  - 2.*gfl(i)*cstar(i)*wdens(i)
     3202!
     3203          d_dens_gen(i) =  d_dens_gen(i)*dtimesub
     3204          d_dens_death(i) = d_dens_death(i)*dtimesub
     3205          d_dens_col(i) =  d_dens_col(i)*dtimesub
     3206          d_wdens(i) = d_dens_gen(i) + d_dens_death(i) + d_dens_col(i)
     3207!!
     3208          d_wdens_targ = max(d_wdens(i), wdensmin-wdens(i))
     3209!!          d_dens_bnd(i) = d_dens_bnd(i) + d_wdens_targ - d_wdens(i)
     3210          d_dens_bnd(i) = d_wdens_targ - d_wdens(i)
     3211          d_wdens(i) = d_wdens_targ
     3212#ifdef IOPHYS_WK
     3213    IF (phys_sub) THEN
     3214        call iophys_ecrit('d_wdens',1,'d_wdens','',d_wdens)
     3215        call iophys_ecrit('d_dens_gen',1,'d_dens_gen','',d_dens_gen)
     3216        call iophys_ecrit('d_dens_death',1,'d_dens_death','',d_dens_death)
     3217        call iophys_ecrit('d_dens_col',1,'d_dens_col','',d_dens_col)
     3218    ENDIF
     3219#endif
     3220
     3221          d_adens_death(i) = -awdens(i)/tau_prime(i)
     3222          d_adens_icol(i) =   2.*igfl(i)*cstar(i)*iwdens(i)
     3223          d_adens_acol(i)  = - 2.*agfl(i)*cstar(i)*awdens(i)
     3224!
     3225          d_adens_death(i) =  d_adens_death(i)*dtimesub
     3226          d_adens_icol(i) =   d_adens_icol(i)*dtimesub
     3227          d_adens_acol(i)  =   d_adens_acol(i)*dtimesub
     3228          d_awdens(i) =   d_dens_gen(i) + d_adens_death(i) + d_adens_icol(i) + d_adens_acol(i)     
     3229#ifdef IOPHYS_WK
     3230    IF (phys_sub) THEN
     3231        call iophys_ecrit('d_awdens',1,'d_awdens','',d_awdens)
     3232        call iophys_ecrit('d_adens_death',1,'d_adens_death','',d_adens_death)
     3233        call iophys_ecrit('d_adens_icol',1,'d_adens_icol','',d_adens_icol)
     3234        call iophys_ecrit('d_adens_acol',1,'d_adens_acol','',d_adens_acol)
     3235    ENDIF
     3236#endif
     3237          d_wdens_targ = min(max(d_awdens(i),-awdens(i)), wdens(i)-awdens(i))
     3238!!          d_dens_bnd(i) = d_dens_bnd(i) + d_wdens_targ - d_wdens(i)
     3239          d_adens_bnd(i) = d_wdens_targ - d_awdens(i)
     3240          d_awdens(i) = d_wdens_targ
     3241
     3242!!          d_irad(i) = (d_sigmaw(i)-d_asigmaw(i)-isigmaw(i)*(d_wdens(i)-awdens(i))/iwdens(i)) / &
     3243!!                      max(smallestreal,(2.*3.14*iwdens(i)*irad_wk(i)))
     3244!!          d_arad(i) = (d_asigmaw(i)-asigmaw(i)*d_awdens(i)/awdens(i)) / &
     3245!!                      max(smallestreal,(2.*3.14*awdens(i)*arad_wk(i)))
     3246!!          d_irad(i) = d_irad(i)*dtimesub
     3247!!          d_arad(i) = d_arad(i)*dtimesub
     3248!!          call iophys_ecrit('d_irad',1,'d_irad','m',d_irad)
     3249!!          call iophys_ecrit('d_airad',1,'d_arad','m',d_arad)
     3250!!
     3251        ENDIF
     3252      ENDDO
     3253
     3254      IF (prt_level >= 10) THEN
     3255        print *,'wake, cstar(1), cstar(1)/cstart, rad_wk(1), tau_wk_inv(1), gfl(1) ', &
     3256                       cstar(1), cstar(1)/cstart, rad_wk(1), tau_wk_inv(1), gfl(1)
     3257        print *,'wake, wdens(1), awdens(1), d_awdens(1) ', &
     3258                       wdens(1), awdens(1), d_awdens(1)
     3259        print *,'wake, d_sig_gen(1), d_sig_death(1), d_sig_col(1), d_sigmaw(1) ', &
     3260                       d_sig_gen(1), d_sig_death(1), d_sig_col(1), d_sigmaw(1)
     3261      ENDIF
     3262sigmaw=sigmaw+d_sigmaw
     3263asigmaw=asigmaw+d_asigmaw
     3264wdens=wdens+d_wdens
     3265awdens=awdens+d_awdens
     3266
     3267    RETURN
     3268    END SUBROUTINE wake_popdyn_3 
     3269
    26893270END MODULE lmdz_wake
  • LMDZ6/trunk/libf/phylmd/lmdz_wake_ini.F90

    r4695 r4744  
    11MODULE lmdz_wake_ini
    22IMPLICIT NONE
    3 
    43
    54  ! ============================================================================
     
    2322!jyg<
    2423!!  REAL, SAVE                                            :: stark, wdens_ref, coefgw, alpk
    25   INTEGER, SAVE, PROTECTED                                         :: prt_level
    26   REAL, SAVE, PROTECTED, DIMENSION(2)                              :: wdens_ref
    27   REAL, SAVE, PROTECTED                                            :: stark, coefgw, alpk, wk_pupper
     24  INTEGER, SAVE, PROTECTED                                    :: prt_level
     25  REAL, SAVE, PROTECTED, DIMENSION(2)                         :: wdens_ref
     26  REAL, SAVE, PROTECTED                                       :: stark, coefgw, alpk, wk_pupper
    2827!>jyg
    29   REAL, SAVE, PROTECTED                                            :: crep_upper, crep_sol 
     28  REAL, SAVE, PROTECTED                                       :: crep_upper, crep_sol 
    3029  !$OMP THREADPRIVATE(stark, wdens_ref, coefgw, alpk, wk_pupper, crep_upper, crep_sol)
    3130
    32   REAL, SAVE, PROTECTED                                            :: tau_cv
     31  REAL, SAVE, PROTECTED                                       :: tau_cv
    3332  !$OMP THREADPRIVATE(tau_cv)
    34   REAL, SAVE, PROTECTED                                            :: rzero, aa0 ! minimal wake radius and area
     33  REAL, SAVE, PROTECTED                                       :: rzero, aa0 ! minimal wake radius and area
    3534  !$OMP THREADPRIVATE(rzero, aa0)
    3635
    37   LOGICAL, SAVE, PROTECTED                                         :: ok_bug_gfl
     36  LOGICAL, SAVE, PROTECTED                                    :: ok_bug_gfl
    3837  !$OMP THREADPRIVATE(ok_bug_gfl)
    39   LOGICAL, SAVE, PROTECTED                                         :: flag_wk_check_trgl
     38  LOGICAL, SAVE, PROTECTED                                    :: flag_wk_check_trgl
    4039  !$OMP THREADPRIVATE(flag_wk_check_trgl)
    41   INTEGER, SAVE, PROTECTED                                         :: iflag_wk_act
     40  INTEGER, SAVE, PROTECTED                                    :: iflag_wk_act
    4241  !$OMP THREADPRIVATE(iflag_wk_act)
    4342
    44   INTEGER, SAVE, PROTECTED                                         :: iflag_wk_check_trgl
     43  INTEGER, SAVE, PROTECTED                                    :: iflag_wk_check_trgl
    4544  !$OMP THREADPRIVATE(iflag_wk_check_trgl)
    46   INTEGER, SAVE, PROTECTED                                         :: iflag_wk_pop_dyn
     45  INTEGER, SAVE, PROTECTED                                    :: iflag_wk_pop_dyn
    4746  !$OMP THREADPRIVATE(iflag_wk_pop_dyn)
    4847
    49   INTEGER, SAVE, PROTECTED                                         :: iflag_wk_profile
     48  INTEGER, SAVE, PROTECTED                                    :: iflag_wk_profile
    5049  !$OMP THREADPRIVATE(iflag_wk_profile)
    5150
    52   REAL, SAVE, PROTECTED                                            :: wdensmin
    53   !$OMP THREADPRIVATE(wdensmin)
    54   REAL, SAVE, PROTECTED                                            :: sigmad, hwmin, wapecut, cstart
     51  REAL, SAVE, PROTECTED                                       :: wdensinit ! Minimum wake density used to restart wakes from a wake-free state
     52  !$OMP THREADPRIVATE(wdensinit)
     53  REAL, SAVE, PROTECTED                                       :: wdensthreshold ! Threshold wake density below which wakes are killed
     54  !$OMP THREADPRIVATE(wdensthreshold)
     55  REAL, SAVE, PROTECTED                                       :: sigmad, hwmin, wapecut, cstart
    5556  !$OMP THREADPRIVATE(sigmad, hwmin, wapecut, cstart)
    56   REAL, SAVE, PROTECTED                                            :: sigmaw_max
     57  REAL, SAVE, PROTECTED                                       :: sigmaw_max
    5758  !$OMP THREADPRIVATE(sigmaw_max) 
    58   REAL, SAVE, PROTECTED                                            :: dens_rate
     59  REAL, SAVE, PROTECTED                                       :: dens_rate
    5960  !$OMP THREADPRIVATE(dens_rate)
    60   REAL, SAVE, PROTECTED                                            :: epsilon_loc
     61  REAL, SAVE, PROTECTED                                       :: epsilon_loc
    6162  !$OMP THREADPRIVATE(epsilon_loc)
    62   REAL, SAVE, PROTECTED                                            :: epsim1,RG,RD
     63  REAL, SAVE, PROTECTED                                       :: epsim1,RG,RD
    6364  !$OMP THREADPRIVATE(epsim1,RG,RD)
     65  REAL, SAVE, PROTECTED                                        ::smallestreal
     66  !$OMP THREADPRIVATE(smallestreal)
    6467
    6568
     
    8992  real, intent(in) :: rg_in,rd_in,rv_in
    9093
     94  smallestreal=tiny(smallestreal)
     95!
    9196  prt_level=prt_lev
    9297  epsilon_loc=1.E-15
    9398  wapecut=1. ! previously 5.
     99!
     100  rzero = 5000.
     101  CALL getin_p('rzero_wk', rzero)
     102  aa0 = 3.14*rzero*rzero
     103!
    94104  ! Essais d'initialisation avec sigmaw = 0.02 et hw = 10.
     105!!  sigmad=0.005
    95106  sigmad=0.02
     107  CALL getin_p('sigmad', sigmad)
    96108  hwmin=10.
    97 !!  DATA wdensmin/1.e-12/
    98   wdensmin=1.e-14
     109!
     110!!wdensthreshold=1.e-12
     111  wdensthreshold=1.e-14
     112  wdensthreshold=2.e-14
     113  CALL getin_p('wdensthreshold', wdensthreshold)
     114!
     115  IF (sigmad < 0.) THEN
     116    sigmad = abs(sigmad)
     117!!    wdensmin=sigmad/(3.14*rzero**2)
     118    wdensinit=sigmad/(3.14*rzero**2)
     119  ELSE
     120    wdensinit = wdensthreshold/2.
     121  ENDIF
     122!
     123!
    99124  ! cc nrlmd
    100125  sigmaw_max=0.4
     
    187212  CALL getin_p('iflag_wk_profile',iflag_wk_profile) ! switch between wdens prescribed
    188213                                                    ! and wdens prognostic
    189   rzero = 5000.
    190   CALL getin_p('rzero_wk', rzero)
    191   aa0 = 3.14*rzero*rzero
    192 !
    193214  tau_cv = 4000.
    194215  CALL getin_p('tau_cv', tau_cv)
     
    223244END SUBROUTINE wake_ini
    224245
     246
     247
    225248END MODULE lmdz_wake_ini
  • LMDZ6/trunk/libf/phylmd/phyetat0_mod.F90

    r4619 r4744  
    2525       solsw, solswfdiff, t_ancien, u_ancien, v_ancien, w01, wake_cstar, wake_deltaq, &
    2626       wake_deltat, wake_delta_pbl_TKE, delta_tsurf, beta_aridity, wake_fip, wake_pe, &
    27        wake_s, wake_dens, awake_dens, cv_gen, zgam, zmax0, zmea, zpic, zsig, &
     27       wake_s, awake_s, wake_dens, awake_dens, cv_gen, zgam, zmax0, zmea, zpic, zsig, &
    2828       zstd, zthe, zval, ale_bl, ale_bl_trig, alp_bl, u10m, v10m, treedrg, &
    2929       ale_wake, ale_bl_stat, ds_ns, dt_ns, delta_sst, delta_sal, dter, dser, &
     
    468468  found=phyetat0_get(wake_deltaq,"WAKE_DELTAQ","Delta hum. wake/env",0.)
    469469  found=phyetat0_get(wake_s,"WAKE_S","Wake frac. area",0.)
     470  found=phyetat0_get(awake_s,"AWAKE_S","Active Wake frac. area",0.)
    470471!jyg<
    471472!  Set wake_dens to -1000. when there is no restart so that the actual
  • LMDZ6/trunk/libf/phylmd/phyredem.F90

    r4613 r4744  
    2222                                v_ancien, clwcon, rnebcon, ratqs, pbl_tke,   &
    2323                                wake_delta_pbl_tke, zmax0, f0, sig1, w01,    &
    24                                 wake_deltat, wake_deltaq, wake_s, wake_dens, &
    25                                 awake_dens, cv_gen,                          &
     24                                wake_deltat, wake_deltaq, wake_s, awake_s,  &
     25                                wake_dens, awake_dens, cv_gen,               &
    2626                                wake_cstar,                                  &
    2727                                wake_pe, wake_fip, fm_therm, entr_therm,     &
     
    301301
    302302    CALL put_field(pass,"WAKE_S", "Wake frac. area", wake_s)
     303
     304    CALL put_field(pass,"AWAKE_S", "Active Wake frac. area", awake_s)
    303305
    304306    CALL put_field(pass,"WAKE_DENS", "Wake num. /unit area", wake_dens)
  • LMDZ6/trunk/libf/phylmd/phys_local_var_mod.F90

    r4742 r4744  
    292292    REAL, SAVE, ALLOCATABLE,DIMENSION(:,:)          :: d_deltat_wk, d_deltaq_wk
    293293!$OMP THREADPRIVATE(d_deltat_wk, d_deltaq_wk)
    294       REAL,ALLOCATABLE,SAVE,DIMENSION(:)            :: d_s_wk, d_dens_a_wk, d_dens_wk
    295 !$OMP THREADPRIVATE(d_s_wk, d_dens_a_wk, d_dens_wk)
     294      REAL,ALLOCATABLE,SAVE,DIMENSION(:)            :: d_s_wk, d_s_a_wk, d_dens_wk, d_dens_a_wk
     295!$OMP THREADPRIVATE(d_s_wk, d_s_a_wk, d_dens_wk, d_dens_a_wk)
    296296    REAL, SAVE, ALLOCATABLE,DIMENSION(:,:)          :: d_deltat_wk_gw, d_deltaq_wk_gw
    297297!$OMP THREADPRIVATE(d_deltat_wk_gw, d_deltaq_wk_gw)
     
    777777      ALLOCATE(wake_omg(klon, klev))
    778778      ALLOCATE(d_deltat_wk(klon, klev), d_deltaq_wk(klon, klev))
    779       ALLOCATE(d_s_wk(klon), d_dens_a_wk(klon), d_dens_wk(klon))
     779      ALLOCATE(d_s_wk(klon), d_s_a_wk(klon), d_dens_wk(klon), d_dens_a_wk(klon))
    780780      ALLOCATE(d_deltat_wk_gw(klon, klev), d_deltaq_wk_gw(klon, klev))
    781781      ALLOCATE(d_deltat_vdf(klon, klev), d_deltaq_vdf(klon, klev))
     
    11061106      DEALLOCATE(wake_omg)
    11071107      DEALLOCATE(d_deltat_wk, d_deltaq_wk)
    1108       DEALLOCATE(d_s_wk, d_dens_a_wk, d_dens_wk)
     1108      DEALLOCATE(d_s_wk, d_s_a_wk, d_dens_wk, d_dens_a_wk)
    11091109      DEALLOCATE(d_deltat_wk_gw, d_deltaq_wk_gw)
    11101110      DEALLOCATE(d_deltat_vdf, d_deltaq_vdf)
  • LMDZ6/trunk/libf/phylmd/phys_state_var_mod.F90

    r4677 r4744  
    263263! wake_deltaq : ecart d'humidite avec la zone non perturbee
    264264! wake_s      : fraction surfacique occupee par la poche froide
     265! awake_s     : surface fraction covered by active wakes
     266! wake_dens   : number of wakes per unit area
    265267! awake_dens  : number of active wakes per unit area
    266 ! wake_dens   : number of wakes per unit area
    267268! cv_gen      : birth rate of cumulonimbus per unit area.
    268269! wake_occ    : occurence of wakes (= 1 if wakes occur, =0 otherwise)
     
    278279!$OMP THREADPRIVATE(wake_deltaxt)
    279280#endif
    280       REAL,ALLOCATABLE,SAVE :: wake_s(:)
    281 !$OMP THREADPRIVATE(wake_s)
    282       REAL,ALLOCATABLE,SAVE :: awake_dens(:), wake_dens(:)
    283 !$OMP THREADPRIVATE(awake_dens, wake_dens)
     281      REAL,ALLOCATABLE,SAVE :: wake_s(:), awake_s(:)
     282!$OMP THREADPRIVATE(wake_s, awake_s)
     283      REAL,ALLOCATABLE,SAVE :: wake_dens(:), awake_dens(:)
     284!$OMP THREADPRIVATE(wake_dens, awake_dens)
    284285      REAL,ALLOCATABLE,SAVE :: cv_gen(:)
    285286!$OMP THREADPRIVATE(cv_gen)
     
    670671      ALLOCATE(wght_th(klon,klev))
    671672      ALLOCATE(wake_deltat(klon,klev), wake_deltaq(klon,klev))
    672       ALLOCATE(wake_s(klon), awake_dens(klon), wake_dens(klon))
     673      ALLOCATE(wake_s(klon), awake_s(klon), wake_dens(klon), awake_dens(klon))
    673674!!      awake_dens = 0.  ! initialized in phyetat0
    674675      ALLOCATE(cv_gen(klon))
     
    864865      DEALLOCATE(lalim_conv, wght_th)
    865866      DEALLOCATE(wake_deltat, wake_deltaq)
    866       DEALLOCATE(wake_s, awake_dens, wake_dens)
     867      DEALLOCATE(wake_s, awake_s, wake_dens, awake_dens)
    867868      DEALLOCATE(cv_gen)
    868869      DEALLOCATE(wake_Cstar, wake_pe, wake_fip)
  • LMDZ6/trunk/libf/phylmd/physiq_mod.F90

    r4742 r4744  
    292292       d_deltat_ajs_cv, d_deltaq_ajs_cv, & ! due to dry adjustment of (w) before convection
    293293                       ! tendencies of wake fractional area and wake number per unit area:
    294        d_s_wk,  d_dens_a_wk,  d_dens_wk, &  ! due to wakes
     294       d_s_wk, d_s_a_wk, d_dens_wk,  d_dens_a_wk, &  ! due to wakes
    295295!!!       d_s_vdf, d_dens_a_vdf, d_dens_vdf, & ! due to vertical diffusion
    296296!!!       d_s_the, d_dens_a_the, d_dens_the, & ! due to thermals
     
    20202020             d_deltaq_ajs_cv(:,:) = 0.
    20212021             d_s_wk(:) = 0.
     2022             d_s_a_wk(:) = 0.
    20222023             d_dens_wk(:) = 0.
     2024             d_dens_a_wk(:) = 0.
    20232025          ENDIF  !  (iflag_wake>=1)
    20242026
     
    28942896          d_deltaq_vdf(:,:) = d_q_vdf_w(:,:)-d_q_vdf_x(:,:)
    28952897          CALL add_wake_tend &
    2896              (d_deltat_vdf, d_deltaq_vdf, dsig0, ddens0, ddens0, wkoccur1, 'vdf', abortphy)
     2898             (d_deltat_vdf, d_deltaq_vdf, dsig0, dsig0, ddens0, ddens0, wkoccur1, 'vdf', abortphy)
    28972899       ELSE
    28982900          d_deltat_vdf(:,:) = 0.
     
    31563158             IF (iflag_adjwk == 2 .AND. OK_bug_ajs_cv) THEN
    31573159               CALL add_wake_tend &
    3158                  (d_deltat_ajs_cv, d_deltaq_ajs_cv, dsig0, ddens0, ddens0, wkoccur1, 'ajs_cv', abortphy)
     3160                 (d_deltat_ajs_cv, d_deltaq_ajs_cv, dsig0, dsig0, ddens0, ddens0, wkoccur1, 'ajs_cv', abortphy)
    31593161             ENDIF  ! (iflag_adjwk == 2 .AND. OK_bug_ajs_cv)
    31603162          ENDIF  ! (iflag_adjwk >= 1)
     
    35313533               dt_a, dq_a, cv_gen,  &
    35323534               sigd, cin,  &
    3533                wake_deltat, wake_deltaq, wake_s, awake_dens, wake_dens,  &
     3535               wake_deltat, wake_deltaq, wake_s, awake_s, wake_dens, awake_dens,  &
    35343536               wake_dth, wake_h,  &
    35353537!!               wake_pe, wake_fip, wake_gfl,  &
     
    35413543               wake_omg, wake_dp_deltomg,  &
    35423544               wake_spread, wake_Cstar, d_deltat_wk_gw,  &
    3543                d_deltat_wk, d_deltaq_wk, d_s_wk, d_dens_a_wk, d_dens_wk)
     3545               d_deltat_wk, d_deltaq_wk, d_s_wk, d_s_a_wk, d_dens_wk, d_dens_a_wk)
    35443546          !
    35453547          !jyg    Reinitialize itapwk when wakes have been called
     
    35603562
    35613563         CALL add_wake_tend &
    3562             (d_deltat_wk, d_deltaq_wk, d_s_wk, d_dens_a_wk, d_dens_wk, wake_k, &
     3564            (d_deltat_wk, d_deltaq_wk, d_s_wk, d_s_a_wk, d_dens_wk, d_dens_a_wk, wake_k, &
    35633565             'wake', abortphy)
    35643566          CALL prt_enerbil('wake',itap)
     
    37373739             IF (ok_bug_split_th) THEN
    37383740               CALL add_wake_tend &
    3739                    (d_deltat_the, d_deltaq_the, dsig0, ddens0, ddens0, wkoccur1, 'the', abortphy)
     3741                   (d_deltat_the, d_deltaq_the, dsig0, dsig0, ddens0, ddens0, wkoccur1, 'the', abortphy)
    37403742             ELSE
    37413743               CALL add_wake_tend &
    3742                    (d_deltat_the, d_deltaq_the, dsig0, ddens0, ddens0, wake_k, 'the', abortphy)
     3744                   (d_deltat_the, d_deltaq_the, dsig0, dsig0, ddens0, ddens0, wake_k, 'the', abortphy)
    37433745             ENDIF
    37443746             CALL prt_enerbil('the',itap)
Note: See TracChangeset for help on using the changeset viewer.