Changeset 2635 for LMDZ5/trunk


Ignore:
Timestamp:
Sep 27, 2016, 6:02:46 PM (8 years ago)
Author:
jyg
Message:

Some cleaning in the wake routines: (i) the wake
number per unit area (wdens) is now a state
variable (held constant for the time being);
(ii) wake state variable changes are computed in
subroutine 'physiq' if iflag_wake_tend=1 (it is
computed within wake routines if
iflag_wake_tend=0, consistent with earlier
versions); (iii) the new routine 'add_wake_tend'
adds tendencies to wake state variables; (iv)
tendencies due to various processes (pbl, wakes,
thermals) are named and added separately.

Location:
LMDZ5/trunk/libf/phylmd
Files:
1 added
9 edited

Legend:

Unmodified
Added
Removed

  • LMDZ5/trunk/libf/phylmd/YOETHF.h

    r2043 r2635  
    2424     &               RALFDCP,RTWAT,RTBER,RTBERCU,                       &
    2525     &               RTICE,RTICECU,RTWAT_RTICE_R,RTWAT_RTICECU_R,RKOOP1,&
    26      &               RKOOP2
     26     &               RKOOP2
    2727
    2828!$OMP THREADPRIVATE(/YOETHF/)

  • LMDZ5/trunk/libf/phylmd/calwake.F90

    r2346 r2635  
    22! $Id$
    33
    4 SUBROUTINE calwake(paprs, pplay, dtime, t, q, omgb, dt_dwn, dq_dwn, m_dwn, &
    5     m_up, dt_a, dq_a, sigd, wdt_pbl, wdq_pbl, udt_pbl, udq_pbl, wake_deltat, &
    6     wake_deltaq, wake_dth, wake_h, wake_s, wake_dens, wake_pe, wake_fip, &
    7     wake_gfl, dt_wake, dq_wake, wake_k, undi_t, undi_q, wake_omgbdth, &
    8     wake_dp_omgb, wake_dtke, wake_dqke, wake_dtpbl, wake_dqpbl, wake_omg, &
    9     wake_dp_deltomg, wake_spread, wake_cstar, wake_d_deltat_gw, wake_ddeltat, &
    10     wake_ddeltaq)
     4SUBROUTINE calwake(iflag_wake_tend, paprs, pplay, dtime, &
     5    t, q, omgb, &
     6    dt_dwn, dq_dwn, m_dwn, m_up, dt_a, dq_a, &
     7    sigd, &
     8    wake_deltat, wake_deltaq, wake_s, wake_dens, &
     9    wake_dth, wake_h, &
     10    wake_pe, wake_fip, wake_gfl, &
     11    dt_wake, dq_wake, wake_k, t_x, q_x, wake_omgbdth, &
     12    wake_dp_omgb, &
     13    wake_dtke, wake_dqke, &
     14    wake_omg, wake_dp_deltomg, &
     15    wake_spread, wake_cstar, wake_d_deltat_gw, &
     16    wake_ddeltat, wake_ddeltaq, wake_ds, wake_ddens)
    1117  ! **************************************************************
    1218  ! *
    … …  
    2632  ! Arguments
    2733  ! ----------
    28 
    29   INTEGER i, l, ktopw(klon)
    30   REAL dtime
    31 
    32   REAL paprs(klon, klev+1), pplay(klon, klev)
    33   REAL t(klon, klev), q(klon, klev), omgb(klon, klev)
    34   REAL dt_dwn(klon, klev), dq_dwn(klon, klev), m_dwn(klon, klev)
    35   REAL m_up(klon, klev)
    36   REAL dt_a(klon, klev), dq_a(klon, klev)
    37   REAL wdt_pbl(klon, klev), wdq_pbl(klon, klev)
    38   REAL udt_pbl(klon, klev), udq_pbl(klon, klev)
    39   REAL wake_deltat(klon, klev), wake_deltaq(klon, klev)
    40   REAL dt_wake(klon, klev), dq_wake(klon, klev)
    41   REAL wake_d_deltat_gw(klon, klev)
    42   REAL wake_h(klon), wake_s(klon)
    43   REAL wake_dth(klon, klev)
    44   REAL wake_pe(klon), wake_fip(klon), wake_gfl(klon)
    45   REAL undi_t(klon, klev), undi_q(klon, klev)
    46   REAL wake_omgbdth(klon, klev), wake_dp_omgb(klon, klev)
    47   REAL wake_dtke(klon, klev), wake_dqke(klon, klev)
    48   REAL wake_dtpbl(klon, klev), wake_dqpbl(klon, klev)
    49   REAL wake_omg(klon, klev), wake_dp_deltomg(klon, klev)
    50   REAL wake_spread(klon, klev), wake_cstar(klon)
    51   REAL wake_ddeltat(klon, klev), wake_ddeltaq(klon, klev)
    52   REAL d_deltatw(klon, klev), d_deltaqw(klon, klev)
    53   INTEGER wake_k(klon)
    54   REAL sigd(klon)
    55   REAL wake_dens(klon)
     34  ! Input
     35  ! ----
     36  INTEGER,                       INTENT (IN)         :: iflag_wake_tend
     37  REAL,                          INTENT (IN)         :: dtime
     38  REAL, DIMENSION(klon, klev),   INTENT (IN)         :: pplay
     39  REAL, DIMENSION(klon, klev+1), INTENT (IN)         :: paprs
     40  REAL, DIMENSION(klon, klev),   INTENT (IN)         :: t, q, omgb
     41  REAL, DIMENSION(klon, klev),   INTENT (IN)         :: dt_dwn, dq_dwn
     42  REAL, DIMENSION(klon, klev),   INTENT (IN)         :: m_up, m_dwn
     43  REAL, DIMENSION(klon, klev),   INTENT (IN)         :: dt_a, dq_a
     44  REAL, DIMENSION(klon),         INTENT (IN)         :: sigd
     45  ! Input/Output
     46  ! ------------
     47  REAL, DIMENSION(klon, klev),   INTENT (INOUT)      :: wake_deltat, wake_deltaq
     48  REAL, DIMENSION(klon),         INTENT (INOUT)      :: wake_s
     49  REAL, DIMENSION(klon),         INTENT (INOUT)      :: wake_dens
     50  ! Output
     51  ! ------
     52  REAL, DIMENSION(klon, klev),   INTENT (OUT)        :: dt_wake, dq_wake
     53  INTEGER, DIMENSION(klon),      INTENT (OUT)        :: wake_k
     54  REAL, DIMENSION(klon, klev),   INTENT (OUT)        :: wake_d_deltat_gw
     55  REAL, DIMENSION(klon),         INTENT (OUT)        :: wake_h
     56  REAL, DIMENSION(klon, klev),   INTENT (OUT)        :: wake_dth
     57  REAL, DIMENSION(klon),         INTENT (OUT)        :: wake_pe, wake_fip, wake_gfl
     58  REAL, DIMENSION(klon, klev),   INTENT (OUT)        :: t_x, q_x
     59  REAL, DIMENSION(klon, klev),   INTENT (OUT)        :: wake_omgbdth, wake_dp_omgb
     60  REAL, DIMENSION(klon, klev),   INTENT (OUT)        :: wake_dtke, wake_dqke
     61  REAL, DIMENSION(klon, klev),   INTENT (OUT)        :: wake_omg, wake_dp_deltomg
     62  REAL, DIMENSION(klon, klev),   INTENT (OUT)        :: wake_spread
     63  REAL, DIMENSION(klon),         INTENT (OUT)        :: wake_cstar
     64  REAL, DIMENSION(klon, klev),   INTENT (OUT)        :: wake_ddeltat, wake_ddeltaq
     65  REAL, DIMENSION(klon),         INTENT (OUT)        :: wake_ds, wake_ddens
     66
    5667
    5768  ! Variable internes
    5869  ! -----------------
    59 
    60   REAL aire
    61   REAL p(klon, klev), ph(klon, klev+1), pi(klon, klev)
    62   REAL te(klon, klev), qe(klon, klev), omgbe(klon, klev+1)
    63   REAL dtdwn(klon, klev), dqdwn(klon, klev)
    64   REAL dta(klon, klev), dqa(klon, klev)
    65   REAL wdtpbl(klon, klev), wdqpbl(klon, klev)
    66   REAL udtpbl(klon, klev), udqpbl(klon, klev)
    67   REAL amdwn(klon, klev), amup(klon, klev)
    68   REAL dtw(klon, klev), dqw(klon, klev), dth(klon, klev)
    69   REAL d_deltat_gw(klon, klev)
    70   REAL dtls(klon, klev), dqls(klon, klev)
    71   REAL tu(klon, klev), qu(klon, klev)
    72   REAL hw(klon), sigmaw(klon), wape(klon), fip(klon), gfl(klon)
    73   REAL omgbdth(klon, klev+1), dp_omgb(klon, klev)
    74   REAL dtke(klon, klev), dqke(klon, klev)
    75   REAL dtpbl(klon, klev), dqpbl(klon, klev)
    76   REAL omg(klon, klev+1), dp_deltomg(klon, klev), spread(klon, klev)
    77   REAL cstar(klon)
    78   REAL sigd0(klon), wdens(klon)
    79 
    80   REAL rdcp
     70  INTEGER                                            :: i, l
     71  REAL                                               :: aire
     72  REAL, DIMENSION(klon, klev)                        :: p,  pi
     73  REAL, DIMENSION(klon, klev+1)                      ::  ph, omgbe
     74  REAL, DIMENSION(klon, klev)                        :: te, qe
     75  REAL, DIMENSION(klon, klev)                        :: dtdwn, dqdwn
     76  REAL, DIMENSION(klon, klev)                        :: dta, dqa
     77  REAL, DIMENSION(klon, klev)                        :: amdwn, amup
     78  REAL, DIMENSION(klon, klev)                        :: dtw, dqw, dth
     79  REAL, DIMENSION(klon, klev)                        :: dtls, dqls
     80  REAL, DIMENSION(klon, klev)                        :: tx, qx
     81  REAL, DIMENSION(klon)                              :: hw, wape, fip, gfl
     82  REAL, DIMENSION(klon)                              :: sigmaw, wdens
     83  REAL, DIMENSION(klon, klev+1)                      :: omgbdth
     84  REAL, DIMENSION(klon, klev)                        :: dp_omgb
     85  REAL, DIMENSION(klon, klev)                        :: dtke, dqke
     86  REAL, DIMENSION(klon, klev+1)                      :: omg
     87  REAL, DIMENSION(klon, klev)                        :: dp_deltomg, spread
     88  REAL, DIMENSION(klon)                              :: cstar
     89  REAL, DIMENSION(klon)                              :: sigd0
     90  INTEGER, DIMENSION(klon)                           :: ktopw
     91  REAL, DIMENSION(klon, klev)                        :: d_deltat_gw
     92  REAL, DIMENSION(klon, klev)                        :: d_deltatw, d_deltaqw
     93  REAL, DIMENSION(klon)                              :: d_sigmaw, d_wdens
     94
     95  REAL                                               :: rdcp
     96
    8197
    8298  ! print *, '-> calwake, wake_s ', wake_s(1)
    … …  
    104120      dta(i, l) = dt_a(i, l)
    105121      dqa(i, l) = dq_a(i, l)
    106       wdtpbl(i, l) = wdt_pbl(i, l)
    107       wdqpbl(i, l) = wdq_pbl(i, l)
    108       udtpbl(i, l) = udt_pbl(i, l)
    109       udqpbl(i, l) = udq_pbl(i, l)
    110122    END DO
    111123  END DO
    … …  
    125137  END DO
    126138
     139  DO i = 1, klon
     140    hw(i) = wake_h(i)
     141  END DO
     142!
     143!    Make a copy of state variables
    127144  DO l = 1, klev
    128145    DO i = 1, klon
    … …  
    132149  END DO
    133150
    134   DO l = 1, klev
    135     DO i = 1, klon
    136       dtls(i, l) = dt_wake(i, l)
    137       dqls(i, l) = dq_wake(i, l)
    138     END DO
    139   END DO
    140 
    141   DO i = 1, klon
    142     hw(i) = wake_h(i)
     151  DO i = 1, klon
    143152    sigmaw(i) = wake_s(i)
     153  END DO
     154
     155  DO i = 1, klon
     156    wdens(i) = wake_dens(i)
    144157  END DO
    145158
    … …  
    166179  END DO
    167180
    168   CALL wake(p, ph, pi, dtime, sigd0, te, qe, omgbe, dtdwn, dqdwn, amdwn, &
    169     amup, dta, dqa, wdtpbl, wdqpbl, udtpbl, udqpbl, dtw, dqw, dth, hw, &
    170     sigmaw, wape, fip, gfl, dtls, dqls, ktopw, omgbdth, dp_omgb, wdens, tu, &
    171     qu, dtke, dqke, dtpbl, dqpbl, omg, dp_deltomg, spread, cstar, &
    172     d_deltat_gw, d_deltatw, d_deltaqw)
    173 
     181  CALL wake(p, ph, pi, dtime, &
     182    te, qe, omgbe, &
     183    dtdwn, dqdwn, amdwn, amup, dta, dqa, &
     184    sigd0, &
     185    dtw, dqw, sigmaw, wdens, &                                   ! state variables
     186    dth, hw, wape, fip, gfl, &
     187    dtls, dqls, ktopw, omgbdth, dp_omgb, tx, qx, &
     188    dtke, dqke, omg, dp_deltomg, spread, cstar, &
     189    d_deltat_gw, &
     190    d_deltatw, d_deltaqw, d_sigmaw, d_wdens)                     ! tendencies
     191
     192!
    174193  DO l = 1, klev
    175194    DO i = 1, klon
    176195      IF (ktopw(i)>0) THEN
    177         wake_deltat(i, l) = dtw(i, l)
    178         wake_deltaq(i, l) = dqw(i, l)
    179196        wake_d_deltat_gw(i, l) = d_deltat_gw(i, l)
    180197        wake_omgbdth(i, l) = omgbdth(i, l)
    … …  
    182199        wake_dtke(i, l) = dtke(i, l)
    183200        wake_dqke(i, l) = dqke(i, l)
    184         wake_dtpbl(i, l) = dtpbl(i, l)
    185         wake_dqpbl(i, l) = dqpbl(i, l)
    186201        wake_omg(i, l) = omg(i, l)
    187202        wake_dp_deltomg(i, l) = dp_deltomg(i, l)
    188203        wake_spread(i, l) = spread(i, l)
    189204        wake_dth(i, l) = dth(i, l)
    190         dt_wake(i, l) = dtls(i, l)
    191         dq_wake(i, l) = dqls(i, l)
    192         undi_t(i, l) = tu(i, l)
    193         undi_q(i, l) = qu(i, l)
    194         wake_ddeltat(i, l) = d_deltatw(i, l)
    195         wake_ddeltaq(i, l) = d_deltaqw(i, l)
     205        dt_wake(i, l) = dtls(i, l)*dtime         ! derivative -> tendency
     206        dq_wake(i, l) = dqls(i, l)*dtime         ! derivative -> tendency
     207        t_x(i, l) = tx(i, l)
     208        q_x(i, l) = qx(i, l)
    196209      ELSE
    197         wake_deltat(i, l) = 0.
    198         wake_deltaq(i, l) = 0.
    199210        wake_d_deltat_gw(i, l) = 0.
    200211        wake_omgbdth(i, l) = 0.
    … …  
    202213        wake_dtke(i, l) = 0.
    203214        wake_dqke(i, l) = 0.
    204         wake_dtpbl(i, l) = 0.
    205         wake_dqpbl(i, l) = 0.
    206215        wake_omg(i, l) = 0.
    207216        wake_dp_deltomg(i, l) = 0.
    … …  
    210219        dt_wake(i, l) = 0.
    211220        dq_wake(i, l) = 0.
    212         undi_t(i, l) = te(i, l)
    213         undi_q(i, l) = qe(i, l)
    214         wake_ddeltat(i, l) = 0.
    215         wake_ddeltaq(i, l) = 0.
     221        t_x(i, l) = te(i, l)
     222        q_x(i, l) = qe(i, l)
    216223      END IF
    217224    END DO
    … …  
    220227  DO i = 1, klon
    221228    wake_h(i) = hw(i)
    222     wake_s(i) = sigmaw(i)
    223229    wake_pe(i) = wape(i)
    224230    wake_fip(i) = fip(i)
    … …  
    226232    wake_k(i) = ktopw(i)
    227233    wake_cstar(i) = cstar(i)
    228     wake_dens(i) = wdens(i)
    229   END DO
     234  END DO
     235
     236!  Tendencies of state variables
     237  DO l = 1, klev
     238    DO i = 1, klon
     239      IF (ktopw(i)>0) THEN
     240        wake_ddeltat(i, l) = d_deltatw(i, l)*dtime
     241        wake_ddeltaq(i, l) = d_deltaqw(i, l)*dtime
     242      ELSE
     243        wake_ddeltat(i, l) = -wake_deltat(i, l)
     244        wake_ddeltaq(i, l) = -wake_deltaq(i, l)
     245      END IF
     246    END DO
     247  END DO
     248  DO i = 1, klon
     249    IF (ktopw(i)>0) THEN
     250      wake_ds(i) = d_sigmaw(i)*dtime
     251      wake_ddens(i) = d_wdens(i)*dtime
     252    ELSE
     253      wake_ds(i)   = -wake_s(i)
     254      wake_ddens(i)= -wake_dens(i)
     255    END IF
     256  END DO
     257
     258!jyg< 
     259  IF (iflag_wake_tend .EQ. 0) THEN
     260!  Update State variables
     261    DO l = 1, klev
     262      DO i = 1, klon
     263        IF (ktopw(i)>0) THEN
     264          wake_deltat(i, l) = dtw(i, l)
     265          wake_deltaq(i, l) = dqw(i, l)
     266        ELSE
     267          wake_deltat(i, l) = 0.
     268          wake_deltaq(i, l) = 0.
     269        END IF
     270      END DO
     271    END DO
     272    DO i = 1, klon
     273      wake_s(i) = sigmaw(i)
     274      wake_dens(i) = wdens(i)
     275    END DO
     276  ENDIF
     277!>jyg
    230278
    231279  RETURN
    232280END SUBROUTINE calwake
     281

  • LMDZ5/trunk/libf/phylmd/dyn1d/lmdz1d.F90

    r2611 r2635  
    1616       ftsol, pbl_tke, pctsrf, radsol, rain_fall, snow_fall, ratqs, &
    1717       rnebcon, rugoro, sig1, w01, solaire_etat0, sollw, sollwdown, &
    18        solsw, t_ancien, q_ancien, u_ancien, v_ancien, wake_cstar, wake_deltaq, &
    19        wake_deltat, wake_delta_pbl_TKE, delta_tsurf, wake_fip, wake_pe, &
    20        wake_s, zgam, &
    21        zmax0, zmea, zpic, zsig, &
     18       solsw, t_ancien, q_ancien, u_ancien, v_ancien, wake_cstar, &
     19       wake_delta_pbl_TKE, delta_tsurf, wake_fip, wake_pe, &
     20       wake_deltaq, wake_deltat, wake_s, wake_dens, &
     21       zgam, zmax0, zmea, zpic, zsig, &
    2222       zstd, zthe, zval, ale_bl, ale_bl_trig, alp_bl
    23    use dimphy
    24    use surface_data, only : type_ocean,ok_veget
    25    use pbl_surface_mod, only : ftsoil, pbl_surface_init,                    &
    26      &                            pbl_surface_final
    27       use fonte_neige_mod, only : fonte_neige_init, fonte_neige_final
    28 
    29    use infotrac ! new
    30    use control_mod
     23   USE dimphy
     24   USE surface_data, only : type_ocean,ok_veget
     25   USE pbl_surface_mod, only : ftsoil, pbl_surface_init, &
     26                                 pbl_surface_final
     27   USE fonte_neige_mod, only : fonte_neige_init, fonte_neige_final
     28
     29   USE infotrac ! new
     30   USE control_mod
    3131   USE indice_sol_mod
    3232   USE phyaqua_mod
    … …  
    787787        wake_pe = 0.
    788788        wake_s = 0.
     789        wake_dens = 0.
    789790        ale_bl = 0.
    790791        ale_bl_trig = 0.
    … …  
    810811! t_ancien,q_ancien,,frugs(:,is_oce),clwcon(:,1),rnebcon(:,1),ratqs(:,1)
    811812! run_off_lic_0,pbl_tke(:,1:klev,nsrf), zmax0,f0,sig1,w01
    812 ! wake_deltat,wake_deltaq,wake_s,wake_cstar,wake_fip,wake_delta_pbl_tke(:,1:klev,nsrf)
     813! wake_deltat,wake_deltaq,wake_s,wake_dens,wake_cstar,
     814! wake_fip,wake_delta_pbl_tke(:,1:klev,nsrf)
    813815!
    814816! NB2: The content of the startphy.nc file depends on some flags defined in

  • LMDZ5/trunk/libf/phylmd/phyetat0.F90

    r2569 r2635  
    1717       solsw, t_ancien, u_ancien, v_ancien, w01, wake_cstar, wake_deltaq, &
    1818       wake_deltat, wake_delta_pbl_TKE, delta_tsurf, wake_fip, wake_pe, &
    19        wake_s, zgam, zmax0, zmea, zpic, zsig, &
     19       wake_s, wake_dens, zgam, zmax0, zmea, zpic, zsig, &
    2020       zstd, zthe, zval, ale_bl, ale_bl_trig, alp_bl, u10m, v10m
    2121  USE geometry_mod, ONLY : longitude_deg, latitude_deg
    … …  
    377377  found=phyetat0_get(klev,wake_deltat,"WAKE_DELTAT","Delta T wake/env",0.)
    378378  found=phyetat0_get(klev,wake_deltaq,"WAKE_DELTAQ","Delta hum. wake/env",0.)
    379   found=phyetat0_get(1,wake_s,"WAKE_S","WAKE_S",0.)
     379  found=phyetat0_get(1,wake_s,"WAKE_S","Wake frac. area",0.)
     380  found=phyetat0_get(1,wake_dens,"WAKE_DENS","Wake num. /unit area",0.)
    380381  found=phyetat0_get(1,wake_cstar,"WAKE_CSTAR","WAKE_CSTAR",0.)
    381382  found=phyetat0_get(1,wake_pe,"WAKE_PE","WAKE_PE",0.)

  • LMDZ5/trunk/libf/phylmd/phyredem.F90

    r2569 r2635  
    1818                                v_ancien, clwcon, rnebcon, ratqs, pbl_tke,   &
    1919                                wake_delta_pbl_tke, zmax0, f0, sig1, w01,    &
    20                                 wake_deltat, wake_deltaq, wake_s, wake_cstar,&
     20                                wake_deltat, wake_deltaq, wake_s, wake_dens, &
     21                                wake_cstar,                                  &
    2122                                wake_pe, wake_fip, fm_therm, entr_therm,     &
    2223                                detr_therm, Ale_bl, Ale_bl_trig, Alp_bl,     &
    … …  
    258259  CALL put_field("WAKE_DELTAQ", "WAKE_DELTAQ", wake_deltaq)
    259260
    260   CALL put_field("WAKE_S", "WAKE_S", wake_s)
     261  CALL put_field("WAKE_S", "Wake frac. area", wake_s)
     262
     263  CALL put_field("WAKE_DENS", "Wake num. /unit area", wake_dens)
    261264
    262265  CALL put_field("WAKE_CSTAR", "WAKE_CSTAR", wake_cstar)

  • LMDZ5/trunk/libf/phylmd/phys_local_var_mod.F90

    r2607 r2635  
    253253      REAL,ALLOCATABLE,SAVE,DIMENSION(:) :: sens, flwp, fiwp
    254254!$OMP THREADPRIVATE(sens, flwp, fiwp)
    255       REAL,ALLOCATABLE,SAVE,DIMENSION(:) :: ale_wake, alp_wake, bils
    256 !$OMP THREADPRIVATE(ale_wake, alp_wake, bils)
     255!!
     256!!         Wake variables
     257      REAL,ALLOCATABLE,SAVE,DIMENSION(:)            :: ale_wake, alp_wake
     258!$OMP THREADPRIVATE(ale_wake, alp_wake)           
     259      REAL,ALLOCATABLE,SAVE,DIMENSION(:)            :: wake_h
     260!$OMP THREADPRIVATE(wake_h)                       
     261    REAL,ALLOCATABLE,SAVE,DIMENSION(:,:)            :: wake_omg
     262!$OMP THREADPRIVATE(wake_omg)                     
     263    REAL, SAVE, ALLOCATABLE,DIMENSION(:,:)          :: d_deltat_wk, d_deltaq_wk
     264!$OMP THREADPRIVATE(d_deltat_wk, d_deltaq_wk)
     265      REAL,ALLOCATABLE,SAVE,DIMENSION(:)            :: d_s_wk, d_dens_wk
     266!$OMP THREADPRIVATE(d_s_wk, d_dens_wk)
     267    REAL, SAVE, ALLOCATABLE,DIMENSION(:,:)          :: d_deltat_wk_gw, d_deltaq_wk_gw
     268!$OMP THREADPRIVATE(d_deltat_wk_gw, d_deltaq_wk_gw)
     269    REAL, SAVE, ALLOCATABLE,DIMENSION(:,:)          :: d_deltat_vdf, d_deltaq_vdf
     270!$OMP THREADPRIVATE(d_deltat_vdf, d_deltaq_vdf)
     271!!!      REAL,ALLOCATABLE,SAVE,DIMENSION(:)          :: d_s_vdf, d_dens_vdf
     272!!!OMP THREADPRIVATE(d_s_vdf, d_dens_vdf)
     273    REAL, SAVE, ALLOCATABLE,DIMENSION(:,:)          :: d_deltat_the, d_deltaq_the
     274!$OMP THREADPRIVATE(d_deltat_the, d_deltaq_the)
     275!!!      REAL,ALLOCATABLE,SAVE,DIMENSION(:)          :: d_s_the, d_dens_the
     276!!!OMP THREADPRIVATE(d_s_the, d_dens_the)
     277      REAL,ALLOCATABLE,SAVE,DIMENSION(:,:)           :: d_deltat_ajs_cv, d_deltaq_ajs_cv
     278!$OMP THREADPRIVATE(d_deltat_ajs_cv, d_deltaq_ajs_cv)                       
     279!!         End of Wake variables
     280!!
     281      REAL,ALLOCATABLE,SAVE,DIMENSION(:) :: bils
     282!$OMP THREADPRIVATE(bils)
    257283      REAL,ALLOCATABLE,SAVE,DIMENSION(:) :: cdragm, cdragh
    258284!$OMP THREADPRIVATE(cdragm, cdragh)
    … …  
    312338      REAL,ALLOCATABLE,SAVE,DIMENSION(:,:) :: dqvdf_x, dqvdf_w
    313339!$OMP THREADPRIVATE(dqvdf_x, dqvdf_w)
    314       REAL,ALLOCATABLE,SAVE,DIMENSION(:,:) :: undi_tke, wake_tke
    315 !$OMP THREADPRIVATE(undi_tke, wake_tke)
    316340! Variables supplémentaires dans physiq.F relative au splitting de la surface
    317341      REAL,ALLOCATABLE,SAVE,DIMENSION(:,:,:) :: pbl_tke_input
    … …  
    330354!>jyg+nrlmd
    331355  !
    332       REAL,ALLOCATABLE,SAVE,DIMENSION(:) :: wake_h, wbeff, zmax_th, zq2m, zt2m
    333 !$OMP THREADPRIVATE(wake_h, wbeff, zmax_th, zq2m, zt2m)
     356      REAL,ALLOCATABLE,SAVE,DIMENSION(:) :: wbeff, zmax_th, zq2m, zt2m
     357!$OMP THREADPRIVATE(wbeff, zmax_th, zq2m, zt2m)
    334358      REAL,ALLOCATABLE,SAVE,DIMENSION(:) :: zt2m_min_mon, zt2m_max_mon
    335359!$OMP THREADPRIVATE(zt2m_min_mon, zt2m_max_mon)
    … …  
    365389      REAL,ALLOCATABLE,SAVE,DIMENSION(:,:) :: ref_liq_pi, ref_ice_pi
    366390!$OMP THREADPRIVATE(ref_liq_pi, ref_ice_pi)
    367       REAL,ALLOCATABLE,SAVE,DIMENSION(:,:) :: wake_omg, zx_rh
    368 !$OMP THREADPRIVATE(wake_omg, zx_rh)
     391      REAL,ALLOCATABLE,SAVE,DIMENSION(:,:) :: zx_rh
     392!$OMP THREADPRIVATE(zx_rh)
    369393      REAL,ALLOCATABLE,SAVE,DIMENSION(:,:) :: pmflxr, pmflxs, prfl, psfl, fraca
    370394!$OMP THREADPRIVATE(pmflxr, pmflxs, prfl, psfl, fraca)
    … …  
    530554      ALLOCATE(tal1(klon), pal1(klon), pab1(klon), pab2(klon))
    531555      ALLOCATE(ptstar(klon),pt0(klon),slp(klon))
    532       ALLOCATE(ale_wake(klon), alp_wake(klon), bils(klon))
     556!!
     557!!          Wake variables
     558      ALLOCATE(ale_wake(klon), alp_wake(klon))
     559      ALLOCATE(wake_h(klon))
     560      ALLOCATE(wake_omg(klon, klev))
     561      ALLOCATE(d_deltat_wk(klon, klev), d_deltaq_wk(klon, klev))
     562      ALLOCATE(d_s_wk(klon), d_dens_wk(klon))
     563      ALLOCATE(d_deltat_wk_gw(klon, klev), d_deltaq_wk_gw(klon, klev))
     564      ALLOCATE(d_deltat_vdf(klon, klev), d_deltaq_vdf(klon, klev))
     565!!      ALLOCATE( d_s_vdf(klon), d_dens_vdf(klon))
     566      ALLOCATE(d_deltat_the(klon, klev), d_deltaq_the(klon, klev))
     567!!      ALLOCATE( d_s_the(klon), d_dens_the(klon))
     568      ALLOCATE(d_deltat_ajs_cv(klon, klev), d_deltaq_ajs_cv(klon, klev))
     569!!         End of wake variables
     570!!
     571      ALLOCATE(bils(klon))
    533572      ALLOCATE(cdragm(klon), cdragh(klon), cldh(klon), cldl(klon))
    534573      ALLOCATE(cldm(klon), cldq(klon), cldt(klon), qsat2m(klon))
    … …  
    561600      ALLOCATE(dtvdf_x(klon,klev), dtvdf_w(klon,klev))
    562601      ALLOCATE(dqvdf_x(klon,klev), dqvdf_w(klon,klev))
    563       ALLOCATE(undi_tke(klon,klev), wake_tke(klon,klev))
    564602      ALLOCATE(pbl_tke_input(klon,klev+1,nbsrf))
    565603      ALLOCATE(t_therm(klon,klev), q_therm(klon,klev),u_therm(klon,klev), v_therm(klon,klev))
    … …  
    568606      ALLOCATE(kh(klon), kh_x(klon), kh_w(klon))
    569607!
    570       ALLOCATE(wake_h(klon), wbeff(klon), zmax_th(klon))
     608      ALLOCATE(wbeff(klon), zmax_th(klon))
    571609      ALLOCATE(zq2m(klon), zt2m(klon), weak_inversion(klon))
    572610      ALLOCATE(zt2m_min_mon(klon), zt2m_max_mon(klon))
    … …  
    589627      ALLOCATE(ref_liq(klon, klev), ref_ice(klon, klev), theta(klon, klev))
    590628      ALLOCATE(ref_liq_pi(klon, klev), ref_ice_pi(klon, klev))
    591       ALLOCATE(zphi(klon, klev), wake_omg(klon, klev), zx_rh(klon, klev))
     629      ALLOCATE(zphi(klon, klev), zx_rh(klon, klev))
    592630      ALLOCATE(pmfd(klon, klev), pmfu(klon, klev))
    593631
    … …  
    741779      DEALLOCATE(tal1, pal1, pab1, pab2)
    742780      DEALLOCATE(ptstar, pt0, slp)
    743       DEALLOCATE(ale_wake, alp_wake, bils)
     781!
     782      DEALLOCATE(ale_wake, alp_wake)
     783      DEALLOCATE(wake_h)
     784      DEALLOCATE(wake_omg)
     785      DEALLOCATE(d_deltat_wk, d_deltaq_wk)
     786      DEALLOCATE(d_s_wk, d_dens_wk)
     787      DEALLOCATE(d_deltat_wk_gw, d_deltaq_wk_gw)
     788      DEALLOCATE(d_deltat_vdf, d_deltaq_vdf)
     789!!      DEALLOCATE( d_s_vdf, d_dens_vdf)
     790      DEALLOCATE(d_deltat_the, d_deltaq_the)
     791!!      DEALLOCATE( d_s_the, d_dens_the)
     792      DEALLOCATE(d_deltat_ajs_cv, d_deltaq_ajs_cv)
     793!
     794      DEALLOCATE(bils)
    744795      DEALLOCATE(cdragm, cdragh, cldh, cldl)
    745796      DEALLOCATE(cldm, cldq, cldt, qsat2m)
    … …  
    770821      DEALLOCATE(dtvdf_x, dtvdf_w)
    771822      DEALLOCATE(dqvdf_x, dqvdf_w)
    772       DEALLOCATE(undi_tke, wake_tke)
    773823      DEALLOCATE(pbl_tke_input)
    774824      DEALLOCATE(t_therm, q_therm, u_therm, v_therm)
    … …  
    777827      DEALLOCATE(kh, kh_x, kh_w)
    778828!
    779       DEALLOCATE(wake_h, wbeff, zmax_th)
     829      DEALLOCATE(wbeff, zmax_th)
    780830      DEALLOCATE(zq2m, zt2m, weak_inversion)
    781831      DEALLOCATE(zt2m_min_mon, zt2m_max_mon)
    … …  
    798848      DEALLOCATE(ref_liq, ref_ice, theta)
    799849      DEALLOCATE(ref_liq_pi, ref_ice_pi)
    800       DEALLOCATE(zphi, wake_omg, zx_rh)
     850      DEALLOCATE(zphi, zx_rh)
    801851      DEALLOCATE(pmfd, pmfu)
    802852

  • LMDZ5/trunk/libf/phylmd/phys_state_var_mod.F90

    r2499 r2635  
    203203      REAL,ALLOCATABLE,SAVE :: cin(:)
    204204!$OMP THREADPRIVATE(cin)
    205 ! ftd : differential heating between wake and environment
     205! ftd : convective heating due to unsaturated downdraughts
    206206      REAL,ALLOCATABLE,SAVE :: ftd(:,:)
    207207!$OMP THREADPRIVATE(ftd)
    208 ! fqd : differential moistening between wake and environment
     208! fqd : convective moistening due to unsaturated downdraughts
    209209      REAL,ALLOCATABLE,SAVE :: fqd(:,:)     
    210210!$OMP THREADPRIVATE(fqd)
    … …  
    232232! wake_deltat : ecart de temperature avec la zone non perturbee
    233233! wake_deltaq : ecart d'humidite avec la zone non perturbee
     234! wake_s      : fraction surfacique occupee par la poche froide
     235! wake_dens   : number of wakes per unit area
     236! wake_occ    : occurence of wakes (= 1 if wakes occur, =0 otherwise)
    234237! wake_Cstar  : vitesse d'etalement de la poche
    235 ! wake_s      : fraction surfacique occupee par la poche froide
    236238! wake_pe     : wake potential energy - WAPE
    237239! wake_fip    : Gust Front Impinging power - ALP
    238 ! dt_wake, dq_wake: LS tendencies due to wake
    239240      REAL,ALLOCATABLE,SAVE :: wake_deltat(:,:)
    240241!$OMP THREADPRIVATE(wake_deltat)
    241242      REAL,ALLOCATABLE,SAVE :: wake_deltaq(:,:)
    242243!$OMP THREADPRIVATE(wake_deltaq)
     244      REAL,ALLOCATABLE,SAVE :: wake_s(:)
     245!$OMP THREADPRIVATE(wake_s)
     246      REAL,ALLOCATABLE,SAVE :: wake_dens(:)
     247!$OMP THREADPRIVATE(wake_dens)
    243248      REAL,ALLOCATABLE,SAVE :: wake_Cstar(:)
    244249!$OMP THREADPRIVATE(wake_Cstar)
    245       REAL,ALLOCATABLE,SAVE :: wake_s(:)
    246 !$OMP THREADPRIVATE(wake_s)
    247250      REAL,ALLOCATABLE,SAVE :: wake_pe(:)
    248251!$OMP THREADPRIVATE(wake_pe)
    249252      REAL,ALLOCATABLE,SAVE :: wake_fip(:)
    250253!$OMP THREADPRIVATE(wake_fip)
    251       REAL,ALLOCATABLE,SAVE :: dt_wake(:,:)
    252 !$OMP THREADPRIVATE(dt_wake)
    253       REAL,ALLOCATABLE,SAVE :: dq_wake(:,:)
    254 !$OMP THREADPRIVATE(dq_wake)
    255254!
    256255!jyg<
    … …  
    524523      ALLOCATE(wght_th(klon,klev))
    525524      ALLOCATE(wake_deltat(klon,klev), wake_deltaq(klon,klev))
    526       ALLOCATE(wake_Cstar(klon), wake_s(klon))
     525      ALLOCATE(wake_s(klon), wake_dens(klon))
     526      ALLOCATE(wake_Cstar(klon))
    527527      ALLOCATE(wake_pe(klon), wake_fip(klon))
    528       ALLOCATE(dt_wake(klon,klev), dq_wake(klon,klev))
    529528!jyg<
    530529      ALLOCATE(wake_delta_pbl_TKE(klon,klev+1,nbsrf+1))
    … …  
    664663      deallocate(lalim_conv, wght_th)
    665664      deallocate(wake_deltat, wake_deltaq)
    666       deallocate(wake_Cstar, wake_s, wake_pe, wake_fip)
    667       deallocate(dt_wake, dq_wake)
     665      deallocate(wake_s, wake_dens)
     666      deallocate(wake_Cstar, wake_pe, wake_fip)
    668667!jyg<
    669668      deallocate(wake_delta_pbl_TKE)

  • LMDZ5/trunk/libf/phylmd/physiq_mod.F90

    r2630 r2635  
    133133       dtvdf_x, dtvdf_w, &
    134134       dqvdf_x, dqvdf_w, &
    135        undi_tke, wake_tke, &
    136135       pbl_tke_input, &
    137136       t_therm, q_therm, u_therm, v_therm, &
    … …  
    141140       !
    142141       ale_wake, alp_wake, &
    143        wake_h, wbeff, zmax_th, &
     142       wake_h, wake_omg, &
     143                       ! tendencies of delta T and delta q:
     144       d_deltat_wk, d_deltaq_wk, &         ! due to wakes
     145       d_deltat_wk_gw, d_deltaq_wk_gw, &   ! due to wake induced gravity waves
     146       d_deltat_vdf, d_deltaq_vdf, &       ! due to vertical diffusion
     147       d_deltat_the, d_deltaq_the, &       ! due to thermals
     148       d_deltat_ajs_cv, d_deltaq_ajs_cv, & ! due to dry adjustment of (w) before convection
     149                       ! tendencies of wake fractional area and wake number per unit area:
     150       d_s_wk,  d_dens_wk, &             ! due to wakes
     151!!!       d_s_vdf, d_dens_vdf, &            ! due to vertical diffusion
     152!!!       d_s_the, d_dens_the, &            ! due to thermals
     153       !                                 
     154       wbeff, zmax_th, &
    144155       sens, flwp, fiwp,  &
    145156       ale_bl_stat,alp_bl_conv,alp_bl_det,  &
    … …  
    157168       ref_liq, ref_ice, theta,  &
    158169       ref_liq_pi, ref_ice_pi,  &
    159        zphi, wake_omg, zx_rh,  &
     170       zphi, zx_rh,  &
    160171       pmfd, pmfu,  &
    161172       !
    … …  
    529540    !RC
    530541    ! Variables li\'ees \`a la poche froide (jyg et rr)
    531     ! Version diagnostique pour l'instant : pas de r\'etroaction sur
    532     ! la convection
    533 
    534     REAL t_wake(klon,klev),q_wake(klon,klev) ! wake pour la convection
     542
     543    INTEGER,  SAVE               :: iflag_wake_tend  ! wake: if =0, then wake state variables are
     544                                                     ! updated within calwake
     545    !$OMP THREADPRIVATE(iflag_wake_tend)
     546    REAL t_w(klon,klev),q_w(klon,klev) ! temperature and moisture profiles in the wake region
     547    REAL t_x(klon,klev),q_x(klon,klev) ! temperature and moisture profiles in the off-wake region
    535548
    536549    REAL wake_dth(klon,klev)        ! wake : temp pot difference
    537550
    538     REAL wake_d_deltat_gw(klon,klev)! wake : delta T tendency due to
    539     ! Gravity Wave (/s)
    540551    REAL wake_omgbdth(klon,klev)    ! Wake : flux of Delta_Theta
    541552    ! transported by LS omega
    … …  
    546557    REAL wake_dqKE(klon,klev)       ! Wake : differential moistening
    547558    ! (wake - unpertubed) CONV
    548     REAL wake_dtPBL(klon,klev)      ! Wake : differential heating
    549     ! (wake - unpertubed) PBL
    550     REAL wake_dqPBL(klon,klev)      ! Wake : differential moistening
    551     ! (wake - unpertubed) PBL
    552     REAL wake_ddeltat(klon,klev),wake_ddeltaq(klon,klev)
    553559    REAL wake_dp_deltomg(klon,klev) ! Wake : gradient vertical de wake_omg
    554560    REAL wake_spread(klon,klev)     ! spreading term in wake_delt
    … …  
    557563    !
    558564    INTEGER wake_k(klon)            ! Wake sommet
    559     !
    560     REAL t_undi(klon,klev)          ! temperature moyenne dans la zone
    561     ! non perturbee
    562     REAL q_undi(klon,klev)          ! humidite moyenne dans la zone
    563     ! non perturbee
    564565    !
    565566    !jyg<
    … …  
    568569
    569570    REAL wake_gfl(klon)             ! Gust Front Length
    570     REAL wake_dens(klon)
     571!!!    REAL wake_dens(klon)         ! moved to phys_state_var_mod
    571572    !
    572573    !
    573574    REAL dt_dwn(klon,klev)
    574575    REAL dq_dwn(klon,klev)
    575     REAL wdt_PBL(klon,klev)
    576     REAL udt_PBL(klon,klev)
    577     REAL wdq_PBL(klon,klev)
    578     REAL udq_PBL(klon,klev)
    579576    REAL M_dwn(klon,klev)
    580577    REAL M_up(klon,klev)
    … …  
    589586    REAL, SAVE :: alp_offset
    590587    !$OMP THREADPRIVATE(alp_offset)
    591 
    592     ! !!
    593     !=================================================================
    594     !         PROVISOIRE : DECOUPLAGE PBL/WAKE
    595     !         --------------------------------
    596     REAL wake_deltat_sav(klon,klev)
    597     REAL wake_deltaq_sav(klon,klev)
    598     !=================================================================
    599588
    600589    !
    … …  
    815804    ! tendance nulles
    816805    REAL, dimension(klon,klev):: du0, dv0, dt0, dq0, dql0, dqi0
     806    REAL, dimension(klon)     :: dsig0, ddens0
     807    INTEGER, dimension(klon)  :: wkoccur1
    817808    !
    818809    ! Flag pour pouvoir ne pas ajouter les tendances.
    … …  
    12091200       CALL getin_p('ratqsp0',ratqsp0)
    12101201       CALL getin_p('ratqsdp',ratqsdp)
     1202       iflag_wake_tend = 0
     1203       CALL getin_p('iflag_wake_tend',iflag_wake_tend)
    12111204    ENDIF
    12121205
    … …  
    16981691    dql0(:,:)=0.
    16991692    dqi0(:,:)=0.
     1693    dsig0(:) = 0.
     1694    ddens0(:) = 0.
     1695    wkoccur1(:)=1
    17001696    !
    17011697    ! Mettre a zero des variables de sortie (pour securite)
    … …  
    21352131       ENDIF
    21362132       ! !!
    2137        !=================================================================
    2138        !         PROVISOIRE : DECOUPLAGE PBL/WAKE
    2139        !         --------------------------------
    2140        !
    2141        !!      wake_deltat_sav(:,:)=wake_deltat(:,:)
    2142        !!      wake_deltaq_sav(:,:)=wake_deltaq(:,:)
    2143        !!      wake_deltat(:,:)=0.
    2144        !!      wake_deltaq(:,:)=0.
    2145        !=================================================================
    21462133       !>jyg+nrlmd
    21472134       !
    … …  
    22182205            )
    22192206       !
    2220        !=================================================================
    2221        !         PROVISOIRE : DECOUPLAGE PBL/WAKE
    2222        !         --------------------------------
    2223        !
    2224        !!      wake_deltat(:,:)=wake_deltat_sav(:,:)
    2225        !!      wake_deltaq(:,:)=wake_deltaq_sav(:,:)
    2226        !=================================================================
    2227        !
    22282207       !  Add turbulent diffusion tendency to the wake difference variables
    22292208       IF (mod(iflag_pbl_split,2) .NE. 0) THEN
    2230           wake_deltat(:,:) = wake_deltat(:,:) + (d_t_vdf_w(:,:)-d_t_vdf_x(:,:))
    2231           wake_deltaq(:,:) = wake_deltaq(:,:) + (d_q_vdf_w(:,:)-d_q_vdf_x(:,:))
     2209!jyg<
     2210          d_deltat_vdf(:,:) = d_t_vdf_w(:,:)-d_t_vdf_x(:,:)
     2211          d_deltaq_vdf(:,:) = d_q_vdf_w(:,:)-d_q_vdf_x(:,:)
     2212          CALL add_wake_tend &
     2213             (d_deltat_vdf, d_deltaq_vdf, dsig0, ddens0, wkoccur1, 'vdf', abortphy)
     2214       ELSE
     2215          d_deltat_vdf(:,:) = 0.
     2216          d_deltaq_vdf(:,:) = 0.
     2217!>jyg
    22322218       ENDIF
    22332219
    … …  
    24042390       !=======================================================================
    24052391       !ajout pour la parametrisation des poches froides: calcul de
    2406        !t_wake et t_undi: si pas de poches froides, t_wake=t_undi=t_seri
    2407        do k=1,klev
    2408           do i=1,klon
    2409              if (iflag_wake>=1) then
    2410                 t_wake(i,k) = t_seri(i,k) &
     2392       !t_w et t_x: si pas de poches froides, t_w=t_x=t_seri
     2393       if (iflag_wake>=1) then
     2394         do k=1,klev
     2395            do i=1,klon
     2396                t_w(i,k) = t_seri(i,k) &
    24112397                     +(1-wake_s(i))*wake_deltat(i,k)
    2412                 q_wake(i,k) = q_seri(i,k) &
     2398                q_w(i,k) = q_seri(i,k) &
    24132399                     +(1-wake_s(i))*wake_deltaq(i,k)
    2414                 t_undi(i,k) = t_seri(i,k) &
     2400                t_x(i,k) = t_seri(i,k) &
    24152401                     -wake_s(i)*wake_deltat(i,k)
    2416                 q_undi(i,k) = q_seri(i,k) &
     2402                q_x(i,k) = q_seri(i,k) &
    24172403                     -wake_s(i)*wake_deltaq(i,k)
    2418              else
    2419                 t_wake(i,k) = t_seri(i,k)
    2420                 q_wake(i,k) = q_seri(i,k)
    2421                 t_undi(i,k) = t_seri(i,k)
    2422                 q_undi(i,k) = q_seri(i,k)
    2423              endif
    2424           enddo
    2425        enddo
     2404            enddo
     2405         enddo
     2406       else
     2407                t_w(:,:) = t_seri(:,:)
     2408                q_w(:,:) = q_seri(:,:)
     2409                t_x(:,:) = t_seri(:,:)
     2410                q_x(:,:) = q_seri(:,:)
     2411       endif
    24262412       !
    24272413       !jyg<
    … …  
    24322418          IF (ok_adjwk) THEN
    24332419             limbas(:) = 1
    2434              CALL ajsec(paprs, pplay, t_wake, q_wake, limbas, &
     2420             CALL ajsec(paprs, pplay, t_w, q_w, limbas, &
    24352421                  d_t_adjwk, d_q_adjwk)
    24362422          ENDIF
    … …  
    24392425             DO i=1,klon
    24402426                IF (wake_s(i) .GT. 1.e-3) THEN
    2441                    t_wake(i,k) = t_wake(i,k) + d_t_adjwk(i,k)
    2442                    q_wake(i,k) = q_wake(i,k) + d_q_adjwk(i,k)
    2443                    wake_deltat(i,k) = wake_deltat(i,k) + d_t_adjwk(i,k)
    2444                    wake_deltaq(i,k) = wake_deltaq(i,k) + d_q_adjwk(i,k)
     2427                   t_w(i,k) = t_w(i,k) + d_t_adjwk(i,k)
     2428                   q_w(i,k) = q_w(i,k) + d_q_adjwk(i,k)
     2429                   d_deltat_ajs_cv(i,k) = d_t_adjwk(i,k)
     2430                   d_deltaq_ajs_cv(i,k) = d_q_adjwk(i,k)
     2431                ELSE
     2432                   d_deltat_ajs_cv(i,k) = 0.
     2433                   d_deltaq_ajs_cv(i,k) = 0.
    24452434                ENDIF
    24462435             ENDDO
    24472436          ENDDO
     2437          CALL add_wake_tend &
     2438              (d_deltat_ajs_cv, d_deltaq_ajs_cv, dsig0, ddens0, wkoccur1, 'ajs_cv', abortphy)
    24482439       ENDIF ! (iflag_wake>=1)
    24492440       !>jyg
    … …  
    24862477          !jyg   iflag_con est dans clesphys
    24872478          !c          CALL concvl (iflag_con,iflag_clos,
    2488           clw=0.
    24892479          CALL concvl (iflag_clos, &
    2490                dtime, paprs, pplay, k_upper_cv, t_undi,q_undi, &
    2491                t_wake,q_wake,wake_s, &
     2480               dtime, paprs, pplay, k_upper_cv, t_x,q_x, &
     2481               t_w,q_w,wake_s, &
    24922482               u_seri,v_seri,tr_seri,nbtr_tmp, &
    24932483               ALE,ALP, &
    … …  
    27052695          ENDDO
    27062696       ENDDO
    2707        !nrlmd+jyg<
    2708        DO k=1,klev
    2709           DO i=1,klon
    2710              wdt_PBL(i,k) =  0.
    2711              wdq_PBL(i,k) =  0.
    2712              udt_PBL(i,k) =  0.
    2713              udq_PBL(i,k) =  0.
    2714           ENDDO
    2715        ENDDO
    2716        !
    2717        IF (mod(iflag_pbl_split,2) .EQ. 1) THEN
    2718           DO k=1,klev
    2719              DO i=1,klon
    2720                 wdt_PBL(i,k) = wdt_PBL(i,k) + d_t_vdf_w(i,k)/dtime
    2721                 wdq_PBL(i,k) = wdq_PBL(i,k) + d_q_vdf_w(i,k)/dtime
    2722                 udt_PBL(i,k) = udt_PBL(i,k) + d_t_vdf_x(i,k)/dtime
    2723                 udq_PBL(i,k) = udq_PBL(i,k) + d_q_vdf_x(i,k)/dtime
    2724                 !!        dt_dwn(i,k)  = dt_dwn(i,k) + d_t_vdf_w(i,k)/dtime
    2725                 !!        dq_dwn(i,k)  = dq_dwn(i,k) + d_q_vdf_w(i,k)/dtime
    2726                 !!        dt_a  (i,k)    = dt_a(i,k) + d_t_vdf_x(i,k)/dtime
    2727                 !!        dq_a  (i,k)    = dq_a(i,k) + d_q_vdf_x(i,k)/dtime
    2728              ENDDO
    2729           ENDDO
    2730        ENDIF
    2731        IF (mod(iflag_pbl_split/2,2) .EQ. 1) THEN
    2732           DO k=1,klev
    2733              DO i=1,klon
    2734                 !!        dt_dwn(i,k)  = dt_dwn(i,k) + 0.
    2735                 !!        dq_dwn(i,k)  = dq_dwn(i,k) + 0.
    2736                 !!        dt_a(i,k)   = dt_a(i,k)   + d_t_ajs(i,k)/dtime
    2737                 !!        dq_a(i,k)   = dq_a(i,k)   + d_q_ajs(i,k)/dtime
    2738                 udt_PBL(i,k)   = udt_PBL(i,k)   + d_t_ajs(i,k)/dtime
    2739                 udq_PBL(i,k)   = udq_PBL(i,k)   + d_q_ajs(i,k)/dtime
    2740              ENDDO
    2741           ENDDO
    2742        ENDIF
    2743        !>nrlmd+jyg
    27442697
    27452698       IF (iflag_wake==2) THEN
    … …  
    27702723       !
    27712724       !calcul caracteristiques de la poche froide
    2772        call calWAKE (paprs,pplay,dtime &
    2773             ,t_seri,q_seri,omega &
    2774             ,dt_dwn,dq_dwn,M_dwn,M_up &
    2775             ,dt_a,dq_a,sigd &
    2776             ,wdt_PBL,wdq_PBL &
    2777             ,udt_PBL,udq_PBL &
    2778             ,wake_deltat,wake_deltaq,wake_dth &
    2779             ,wake_h,wake_s,wake_dens &
    2780             ,wake_pe,wake_fip,wake_gfl &
    2781             ,dt_wake,dq_wake &
    2782             ,wake_k, t_undi,q_undi &
    2783             ,wake_omgbdth,wake_dp_omgb &
    2784             ,wake_dtKE,wake_dqKE &
    2785             ,wake_dtPBL,wake_dqPBL &
    2786             ,wake_omg,wake_dp_deltomg &
    2787             ,wake_spread,wake_Cstar,wake_d_deltat_gw &
    2788             ,wake_ddeltat,wake_ddeltaq)
     2725       call calWAKE (iflag_wake_tend, paprs, pplay, dtime, &
     2726            t_seri, q_seri, omega,  &
     2727            dt_dwn, dq_dwn, M_dwn, M_up,  &
     2728            dt_a, dq_a,  &
     2729            sigd,  &
     2730            wake_deltat, wake_deltaq, wake_s, wake_dens,  &
     2731            wake_dth, wake_h,  &
     2732            wake_pe, wake_fip, wake_gfl,  &
     2733            d_t_wake, d_q_wake,  &
     2734            wake_k, t_x, q_x,  &
     2735            wake_omgbdth, wake_dp_omgb,  &
     2736            wake_dtKE, wake_dqKE,  &
     2737            wake_omg, wake_dp_deltomg,  &
     2738            wake_spread, wake_Cstar, d_deltat_wk_gw,  &
     2739            d_deltat_wk, d_deltaq_wk, d_s_wk, d_dens_wk)
    27892740       !
    27902741       !-----------------------------------------------------------------------
    27912742       ! ajout des tendances des poches froides
    2792        ! Faire rapidement disparaitre l'ancien dt_wake pour garder un d_t_wake
    2793        ! coherent avec les autres d_t_...
    2794        d_t_wake(:,:)=dt_wake(:,:)*dtime
    2795        d_q_wake(:,:)=dq_wake(:,:)*dtime
    27962743       CALL add_phys_tend(du0,dv0,d_t_wake,d_q_wake,dql0,dqi0,paprs,'wake', &
    27972744            abortphy,flag_inhib_tend)
    27982745       !------------------------------------------------------------------------
     2746
     2747!       Increment Wake state variables
     2748       IF (iflag_wake_tend .GT. 0.) THEN
     2749
     2750         CALL add_wake_tend &
     2751            (d_deltat_wk, d_deltaq_wk, d_s_wk, d_dens_wk, wake_k, &
     2752             'wake', abortphy)
     2753
     2754       ENDIF   ! (iflag_wake_tend .GT. 0.)
    27992755
    28002756    endif  ! (iflag_wake>=1)
    … …  
    29192875                DO i=1,klon
    29202876                   !
    2921                    wake_deltat(i,k) = wake_deltat(i,k) - d_t_ajs(i,k)
    2922                    wake_deltaq(i,k) = wake_deltaq(i,k) - d_q_ajs(i,k)
    2923                    !
    2924                    !!!t_seri(i,k) = t_therm(i,k) + wake_s(i)*wake_deltat(i,k)
    2925                    !!!q_seri(i,k) = q_therm(i,k) + wake_s(i)*wake_deltaq(i,k)
     2877                   d_deltat_the(i,k) = - d_t_ajs(i,k)
     2878                   d_deltaq_the(i,k) = - d_q_ajs(i,k)
    29262879                   !
    29272880                   d_u_ajs(i,k) = d_u_ajs(i,k)*(1.-wake_s(i))
    … …  
    29322885                ENDDO
    29332886             ENDDO
    2934 !!!          ELSE
    2935 !!!             DO k=1,klev
    2936 !!!                DO i=1,klon
    2937 !!!                   t_seri(i,k) = t_therm(i,k)
    2938 !!!                   q_seri(i,k) = q_therm(i,k)
    2939 !!!                ENDDO
    2940 !!!             ENDDO
    29412887          ENDIF
     2888!
     2889          CALL add_wake_tend &
     2890              (d_deltat_the, d_deltaq_the, dsig0, ddens0, wkoccur1, 'the', abortphy)
     2891!
    29422892          !
    29432893          CALL add_phys_tend(d_u_ajs,d_v_ajs,d_t_ajs,d_q_ajs,  &
    29442894                             dql0,dqi0,paprs,'thermals', abortphy,flag_inhib_tend)
    29452895          !
    2946           !>jyg
    2947 !jyg<
    29482896!
    29492897          CALL alpale_th( dtime, lmax_th, t_seri, cell_area,  &
    … …  
    29522900                          alp_bl, alp_bl_stat, &
    29532901                          proba_notrig, random_notrig)
     2902          !>jyg
    29542903
    29552904          ! ------------------------------------------------------------------

  • LMDZ5/trunk/libf/phylmd/wake.F90

    r2495 r2635  
    22! $Id$
    33
    4 SUBROUTINE wake(p, ph, pi, dtime, sigd_con, &
     4SUBROUTINE wake(p, ph, pi, dtime, &
    55                te0, qe0, omgb, &
    66                dtdwn, dqdwn, amdwn, amup, dta, dqa, &
    7                 wdtpbl, wdqpbl, udtpbl, udqpbl, &
    8                 deltatw, deltaqw, dth, hw, sigmaw, wape, fip, gfl, &
    9                 dtls, dqls, ktopw, omgbdth, dp_omgb, wdens, tu, qu, &
    10                 dtke, dqke, dtpbl, dqpbl, omg, dp_deltomg, spread, cstar, &
    11                 d_deltat_gw, d_deltatw2, d_deltaqw2)
     7                sigd_con, &
     8                deltatw, deltaqw, sigmaw, wdens, &                          ! state variables
     9                dth, hw, wape, fip, gfl, &
     10                dtls, dqls, ktopw, omgbdth, dp_omgb, tu, qu, &
     11                dtke, dqke, omg, dp_deltomg, spread, cstar, &
     12                d_deltat_gw, &
     13                d_deltatw2, d_deltaqw2, d_sigmaw2, d_wdens2)                ! tendencies
    1214
    1315
    … …  
    3335  ! le declenchement de nouvelles colonnes convectives.
    3436
    35   ! Variables d'etat : deltatw    : ecart de temperature wake-undisturbed
    36   ! area
    37   ! deltaqw    : ecart d'humidite wake-undisturbed area
    38   ! sigmaw     : fraction d'aire occupee par la poche.
     37  ! State variables :
     38  ! deltatw    : temperature difference between wake and off-wake regions
     39  ! deltaqw    : specific humidity difference between wake and off-wake regions
     40  ! sigmaw     : fractional area covered by wakes.
     41  ! wdens      : number of wakes per unit area
    3942
    4043  ! Variable de sortie :
    … …  
    5356  ! omg    : Delta_omg =vertical velocity diff. wake-undist. (Pa/s)
    5457  ! dp_deltomg  : vertical gradient of omg (s-1)
    55   ! spread  : spreading term in dt_wake and dq_wake
     58  ! spread  : spreading term in d_t_wake and d_q_wake
    5659  ! deltatw     : updated temperature difference (T_w-T_u).
    5760  ! deltaqw     : updated humidity difference (q_w-q_u).
    … …  
    126129  REAL, DIMENSION (klon, klev),     INTENT(IN)          :: te0, qe0
    127130  REAL, DIMENSION (klon, klev),     INTENT(IN)          :: dtdwn, dqdwn
    128   REAL, DIMENSION (klon, klev),     INTENT(IN)          :: wdtpbl, wdqpbl, udtpbl, udqpbl ! UNUSED
    129131  REAL, DIMENSION (klon, klev),     INTENT(IN)          :: amdwn, amup
    130132  REAL, DIMENSION (klon, klev),     INTENT(IN)          :: dta, dqa
    … …  
    133135  !
    134136  ! Input/Output
     137  ! State variables
    135138  REAL, DIMENSION (klon, klev),     INTENT(INOUT)       :: deltatw, deltaqw
    136139  REAL, DIMENSION (klon),           INTENT(INOUT)       :: sigmaw
     140  REAL, DIMENSION (klon),           INTENT(INOUT)       :: wdens
    137141
    138142  ! Sorties
    … …  
    143147  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: dtls, dqls
    144148  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: dtke, dqke
    145   REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: dtpbl, dqpbl
    146149  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: spread
    147   REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: d_deltatw2, d_deltaqw2
    148150  REAL, DIMENSION (klon, klev+1),   INTENT(OUT)         :: omgbdth, omg
    149151  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: dp_omgb, dp_deltomg
    150152  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: d_deltat_gw
    151153  REAL, DIMENSION (klon),           INTENT(OUT)         :: hw, wape, fip, gfl, cstar
    152   REAL, DIMENSION (klon),           INTENT(OUT)         :: wdens
    153154  INTEGER, DIMENSION (klon),        INTENT(OUT)         :: ktopw
     155  ! Tendencies of state variables
     156  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: d_deltatw2, d_deltaqw2
     157  REAL, DIMENSION (klon),           INTENT(OUT)         :: d_sigmaw2, d_wdens2
    154158
    155159  ! Variables internes
    … …  
    157161
    158162  ! Variables à fixer
    159   REAL alon
    160   LOGICAL, SAVE :: first = .TRUE.
     163  REAL                                                  :: alon
     164  LOGICAL, SAVE                                         :: first = .TRUE.
    161165  !$OMP THREADPRIVATE(first)
    162   REAL, SAVE ::  stark, wdens_ref, coefgw, alpk, crep_upper, crep_sol 
     166  REAL, SAVE                                            :: stark, wdens_ref, coefgw, alpk
     167  REAL, SAVE                                            :: crep_upper, crep_sol 
    163168  !$OMP THREADPRIVATE(stark, wdens_ref, coefgw, alpk, crep_upper, crep_sol)
    164169
    165   REAL delta_t_min
    166   INTEGER nsub
    167   REAL dtimesub
    168   REAL sigmad, hwmin, wapecut
    169   REAL :: sigmaw_max
    170   REAL :: dens_rate
    171   REAL wdens0
     170  REAL                                                  :: delta_t_min
     171  INTEGER                                               :: nsub
     172  REAL                                                  :: dtimesub
     173  REAL                                                  :: sigmad, hwmin, wapecut
     174  REAL                                                  :: sigmaw_max
     175  REAL                                                  :: dens_rate
     176  REAL                                                  :: wdens0
    172177  ! IM 080208
    173   LOGICAL, DIMENSION (klon) :: gwake
     178  LOGICAL, DIMENSION (klon)                             :: gwake
    174179
    175180  ! Variables de sauvegarde
    176   REAL, DIMENSION (klon, klev) :: deltatw0
    177   REAL, DIMENSION (klon, klev) :: deltaqw0
    178   REAL, DIMENSION (klon, klev) :: te, qe
    179   REAL, DIMENSION (klon) :: sigmaw0, sigmaw1
     181  REAL, DIMENSION (klon, klev)                          :: deltatw0
     182  REAL, DIMENSION (klon, klev)                          :: deltaqw0
     183  REAL, DIMENSION (klon, klev)                          :: te, qe
     184  REAL, DIMENSION (klon)                                :: sigmaw0, sigmaw1
    180185
    181186  ! Variables pour les GW
    182   REAL, DIMENSION (klon) :: ll
    183   REAL, DIMENSION (klon, klev) :: n2
    184   REAL, DIMENSION (klon, klev) :: cgw
    185   REAL, DIMENSION (klon, klev) :: tgw
     187  REAL, DIMENSION (klon)                                :: ll
     188  REAL, DIMENSION (klon, klev)                          :: n2
     189  REAL, DIMENSION (klon, klev)                          :: cgw
     190  REAL, DIMENSION (klon, klev)                          :: tgw
    186191
    187192  ! Variables liées au calcul de hw
    188   REAL, DIMENSION (klon) :: ptop_provis, ptop, ptop_new
    189   REAL, DIMENSION (klon) :: sum_dth
    190   REAL, DIMENSION (klon) :: dthmin
    191   REAL, DIMENSION (klon) :: z, dz, hw0
    192   INTEGER, DIMENSION (klon) :: ktop, kupper
     193  REAL, DIMENSION (klon)                                :: ptop_provis, ptop, ptop_new
     194  REAL, DIMENSION (klon)                                :: sum_dth
     195  REAL, DIMENSION (klon)                                :: dthmin
     196  REAL, DIMENSION (klon)                                :: z, dz, hw0
     197  INTEGER, DIMENSION (klon)                             :: ktop, kupper
    193198
    194199  ! Sub-timestep tendencies and related variables
    195   REAL d_deltatw(klon, klev), d_deltaqw(klon, klev)
    196   REAL d_te(klon, klev), d_qe(klon, klev)
    197   REAL d_sigmaw(klon), alpha(klon)
    198   REAL q0_min(klon), q1_min(klon)
    199   LOGICAL wk_adv(klon), ok_qx_qw(klon)
    200   REAL epsilon
     200  REAL, DIMENSION (klon, klev)                          :: d_deltatw, d_deltaqw
     201  REAL, DIMENSION (klon, klev)                          :: d_te, d_qe
     202  REAL, DIMENSION (klon)                                :: d_sigmaw, alpha
     203  REAL, DIMENSION (klon)                                :: q0_min, q1_min
     204  LOGICAL, DIMENSION (klon)                             :: wk_adv, ok_qx_qw
     205  REAL, SAVE                                            :: epsilon
     206  !$OMP THREADPRIVATE(epsilon)
    201207  DATA epsilon/1.E-15/
    202208
    203209  ! Autres variables internes
    204   INTEGER isubstep, k, i
    205 
    206   REAL, DIMENSION (klon) :: sum_thu, sum_tu, sum_qu, sum_thvu
    207   REAL, DIMENSION (klon) :: sum_dq, sum_rho
    208   REAL, DIMENSION (klon) :: sum_dtdwn, sum_dqdwn
    209   REAL, DIMENSION (klon) :: av_thu, av_tu, av_qu, av_thvu
    210   REAL, DIMENSION (klon) :: av_dth, av_dq, av_rho
    211   REAL, DIMENSION (klon) :: av_dtdwn, av_dqdwn
    212 
    213   REAL, DIMENSION (klon, klev) :: rho, rhow
    214   REAL, DIMENSION (klon, klev+1) :: rhoh
    215   REAL, DIMENSION (klon, klev) :: rhow_moyen
    216   REAL, DIMENSION (klon, klev) :: zh
    217   REAL, DIMENSION (klon, klev+1) :: zhh
    218   REAL, DIMENSION (klon, klev) :: epaisseur1, epaisseur2
    219 
    220   REAL, DIMENSION (klon, klev) :: the, thu
    221 
    222   ! REAL, DIMENSION(klon,klev) :: d_deltatw, d_deltaqw
    223 
    224   REAL, DIMENSION (klon, klev+1) :: omgbw
    225   REAL, DIMENSION (klon) :: pupper
    226   REAL, DIMENSION (klon) :: omgtop
    227   REAL, DIMENSION (klon, klev) :: dp_omgbw
    228   REAL, DIMENSION (klon) :: ztop, dztop
    229   REAL, DIMENSION (klon, klev) :: alpha_up
    230 
    231   REAL, DIMENSION (klon) :: rre1, rre2
    232   REAL :: rrd1, rrd2
    233   REAL, DIMENSION (klon, klev) :: th1, th2, q1, q2
    234   REAL, DIMENSION (klon, klev) :: d_th1, d_th2, d_dth
    235   REAL, DIMENSION (klon, klev) :: d_q1, d_q2, d_dq
    236   REAL, DIMENSION (klon, klev) :: omgbdq
    237 
    238   REAL, DIMENSION (klon) :: ff, gg
    239   REAL, DIMENSION (klon) :: wape2, cstar2, heff
    240 
    241   REAL, DIMENSION (klon, klev) :: crep
    242 
    243   REAL, DIMENSION (klon, klev) :: ppi
     210  INTEGER                                               ::isubstep, k, i
     211
     212  REAL                                                  :: sigmaw_targ
     213
     214  REAL, DIMENSION (klon)                                :: sum_thu, sum_tu, sum_qu, sum_thvu
     215  REAL, DIMENSION (klon)                                :: sum_dq, sum_rho
     216  REAL, DIMENSION (klon)                                :: sum_dtdwn, sum_dqdwn
     217  REAL, DIMENSION (klon)                                :: av_thu, av_tu, av_qu, av_thvu
     218  REAL, DIMENSION (klon)                                :: av_dth, av_dq, av_rho
     219  REAL, DIMENSION (klon)                                :: av_dtdwn, av_dqdwn
     220
     221  REAL, DIMENSION (klon, klev)                          :: rho, rhow
     222  REAL, DIMENSION (klon, klev+1)                        :: rhoh
     223  REAL, DIMENSION (klon, klev)                          :: rhow_moyen
     224  REAL, DIMENSION (klon, klev)                          :: zh
     225  REAL, DIMENSION (klon, klev+1)                        :: zhh
     226  REAL, DIMENSION (klon, klev)                          :: epaisseur1, epaisseur2
     227
     228  REAL, DIMENSION (klon, klev)                          :: the, thu
     229
     230  REAL, DIMENSION (klon, klev+1)                        :: omgbw
     231  REAL, DIMENSION (klon)                                :: pupper
     232  REAL, DIMENSION (klon)                                :: omgtop
     233  REAL, DIMENSION (klon, klev)                          :: dp_omgbw
     234  REAL, DIMENSION (klon)                                :: ztop, dztop
     235  REAL, DIMENSION (klon, klev)                          :: alpha_up
     236
     237  REAL, DIMENSION (klon)                                :: rre1, rre2
     238  REAL                                                  :: rrd1, rrd2
     239  REAL, DIMENSION (klon, klev)                          :: th1, th2, q1, q2
     240  REAL, DIMENSION (klon, klev)                          :: d_th1, d_th2, d_dth
     241  REAL, DIMENSION (klon, klev)                          :: d_q1, d_q2, d_dq
     242  REAL, DIMENSION (klon, klev)                          :: omgbdq
     243
     244  REAL, DIMENSION (klon)                                :: ff, gg
     245  REAL, DIMENSION (klon)                                :: wape2, cstar2, heff
     246                                                       
     247  REAL, DIMENSION (klon, klev)                          :: crep
     248                                                       
     249  REAL, DIMENSION (klon, klev)                          :: ppi
    244250
    245251  ! cc nrlmd
    246   REAL, DIMENSION (klon) :: death_rate, nat_rate
    247   REAL, DIMENSION (klon, klev) :: entr
    248   REAL, DIMENSION (klon, klev) :: detr
     252  REAL, DIMENSION (klon)                                :: death_rate, nat_rate
     253  REAL, DIMENSION (klon, klev)                          :: entr
     254  REAL, DIMENSION (klon, klev)                          :: detr
     255
     256  REAL, DIMENSION(klon)                                 :: sigmaw_in   ! pour les prints
    249257
    250258  ! -------------------------------------------------------------------------
    … …  
    312320  ! Les densites peuvent evoluer si les poches debordent
    313321  ! (voir au tout debut de la boucle sur les substeps)
    314   wdens = wdens_ref
     322  wdens(:) = wdens_ref
    315323
    316324  ! print*,'stark',stark
    … …  
    347355    END DO
    348356  END DO
     357  DO i = 1, klon
     358   sigmaw_in(i) = sigmaw(i)
     359  END DO
    349360  ! sigmaw1=sigmaw
    350361  ! IF (sigd_con.GT.sigmaw1) THEN
    … …  
    353364  DO i = 1, klon
    354365    ! c      sigmaw(i) = amax1(sigmaw(i),sigd_con(i))
    355     sigmaw(i) = amax1(sigmaw(i), sigmad)
    356     sigmaw(i) = amin1(sigmaw(i), 0.99)
     366!jyg<
     367!!    sigmaw(i) = amax1(sigmaw(i), sigmad)
     368!!    sigmaw(i) = amin1(sigmaw(i), 0.99)
     369    sigmaw_targ = min(max(sigmaw(i), sigmad),0.99)
     370    d_sigmaw2(i) = sigmaw_targ - sigmaw(i)
     371    sigmaw(i) = sigmaw_targ
     372!>jyg
    357373    sigmaw0(i) = sigmaw(i)
    358374    wape(i) = 0.
    359375    wape2(i) = 0.
    360376    d_sigmaw(i) = 0.
     377    d_wdens2(i) = 0.
    361378    ktopw(i) = 0
    362379  END DO
    … …  
    433450        n2(i, k) = 0
    434451      ELSE
    435         n2(i, k) = amax1(0., -rg**2/the(i,k)*rho(i,k)*(the(i,k+1)-the(i, &
    436           k-1))/(p(i,k+1)-p(i,k-1)))
     452        n2(i, k) = amax1(0., -rg**2/the(i,k)*rho(i,k)*(the(i,k+1)-the(i,k-1))/ &
     453                             (p(i,k+1)-p(i,k-1)))
    437454      END IF
    438455      zh(i, k) = (zhh(i,k)+zhh(i,k+1))/2
    … …  
    504521      IF (dth(i,k)>-delta_t_min .AND. dth(i,k-1)<-delta_t_min .AND. &
    505522          ptop_provis(i)==ph(i,1)) THEN
    506         ptop_provis(i) = ((dth(i,k)+delta_t_min)*p(i,k-1)-(dth(i, &
    507           k-1)+delta_t_min)*p(i,k))/(dth(i,k)-dth(i,k-1))
     523        ptop_provis(i) = ((dth(i,k)+delta_t_min)*p(i,k-1)- &
     524                          (dth(i,k-1)+delta_t_min)*p(i,k))/(dth(i,k)-dth(i,k-1))
    508525      END IF
    509526    END DO
    … …  
    602619        deltatw(i, k) = 0.
    603620        deltaqw(i, k) = 0.
     621        d_deltatw2(i,k) = -deltatw0(i,k)
     622        d_deltaqw2(i,k) = -deltaqw0(i,k)
    604623      END IF
    605624    END DO
    … …  
    671690    av_dqdwn(i) = sum_dqdwn(i)/hw0(i)
    672691
    673     wape(i) = -rg*hw0(i)*(av_dth(i)+epsim1*(av_thu(i)*av_dq(i)+av_dth(i)*av_qu(i &
    674       )+av_dth(i)*av_dq(i)))/av_thvu(i)
     692    wape(i) = -rg*hw0(i)*(av_dth(i)+ &
     693        epsim1*(av_thu(i)*av_dq(i)+av_dth(i)*av_qu(i)+av_dth(i)*av_dq(i)))/av_thvu(i)
     694
    675695  END DO
    676696
    … …  
    686706        deltaqw(i, k) = 0.
    687707        dth(i, k) = 0.
     708        d_deltatw2(i,k) = -deltatw0(i,k)
     709        d_deltaqw2(i,k) = -deltaqw0(i,k)
    688710      END IF
    689711    END DO
    … …  
    695717      cstar(i) = 0.
    696718      hw(i) = hwmin
    697       sigmaw(i) = amax1(sigmad, sigd_con(i))
     719!jyg<
     720!!      sigmaw(i) = amax1(sigmad, sigd_con(i))
     721      sigmaw_targ = max(sigmad, sigd_con(i))
     722      d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i)
     723      sigmaw(i) = sigmaw_targ
     724!>jyg
    698725      fip(i) = 0.
    699726      gwake(i) = .FALSE.
    … …  
    708735  ! --------------------------
    709736  DO i = 1, klon
    710     q0_min(i) = min((qe(i,1)-sigmaw(i)*deltaqw(i,1)), (qe(i, &
    711       1)+(1.-sigmaw(i))*deltaqw(i,1)))
     737    q0_min(i) = min((qe(i,1)-sigmaw(i)*deltaqw(i,1)), &
     738                    (qe(i,1)+(1.-sigmaw(i))*deltaqw(i,1)))
    712739  END DO
    713740  DO k = 2, klev
    714741    DO i = 1, klon
    715       q1_min(i) = min((qe(i,k)-sigmaw(i)*deltaqw(i,k)), (qe(i, &
    716         k)+(1.-sigmaw(i))*deltaqw(i,k)))
     742      q1_min(i) = min((qe(i,k)-sigmaw(i)*deltaqw(i,k)), &
     743                      (qe(i,k)+(1.-sigmaw(i))*deltaqw(i,k)))
    717744      IF (q1_min(i)<=q0_min(i)) THEN
    718745        q0_min(i) = q1_min(i)
    … …  
    752779      IF (wk_adv(i) .AND. cstar(i)>0.01) THEN
    753780        omg(i, kupper(i)+1) = -rg*amdwn(i, kupper(i)+1)/sigmaw(i) + &
    754           rg*amup(i, kupper(i)+1)/(1.-sigmaw(i))
    755         wdens0 = (sigmaw(i)/(4.*3.14))*((1.-sigmaw(i))*omg(i,kupper(i)+1)/(( &
    756           ph(i,1)-pupper(i))*cstar(i)))**(2)
     781                               rg*amup(i, kupper(i)+1)/(1.-sigmaw(i))
     782        wdens0 = (sigmaw(i)/(4.*3.14))* &
     783          ((1.-sigmaw(i))*omg(i,kupper(i)+1)/((ph(i,1)-pupper(i))*cstar(i)))**(2)
    757784        IF (wdens(i)<=wdens0*1.1) THEN
    758785          wdens(i) = wdens0
    … …  
    770797      IF (wk_adv(i)) THEN
    771798        gfl(i) = 2.*sqrt(3.14*wdens(i)*sigmaw(i))
    772         sigmaw(i) = amin1(sigmaw(i), sigmaw_max)
    773       END IF
    774     END DO
     799!jyg<
     800!!        sigmaw(i) = amin1(sigmaw(i), sigmaw_max)
     801        sigmaw_targ = min(sigmaw(i), sigmaw_max)
     802        d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i)
     803        sigmaw(i) = sigmaw_targ
     804!>jyg
     805      END IF
     806    END DO
     807
    775808    DO i = 1, klon
    776809      IF (wk_adv(i)) THEN
    … …  
    783816          death_rate(i) = 0.
    784817        END IF
     818
    785819        d_sigmaw(i) = gfl(i)*cstar(i)*dtimesub - death_rate(i)*sigmaw(i)* &
    786820          dtimesub
    … …  
    798832      END IF
    799833    END DO
    800 
    801834
    802835    ! calcul de la difference de vitesse verticale poche - zone non perturbee
    … …  
    10061039
    10071040          d_deltatw(i, k) = dtimesub/(ph(i,k)-ph(i,k+1))* &
    1008             (rrd1*omg(i,k)*sigmaw(i)*d_th1(i,k)-rrd2*omg(i,k+1)*(1.-sigmaw( &
    1009             i))*d_th2(i,k+1)-(1.-alpha_up(i,k))*omgbdth(i,k)-alpha_up(i,k+1)* &
    1010             omgbdth(i,k+1))*ppi(i, k)
     1041            (rrd1*omg(i,k)*sigmaw(i)*d_th1(i,k) - &
     1042             rrd2*omg(i,k+1)*(1.-sigmaw(i))*d_th2(i,k+1)- &
     1043             (1.-alpha_up(i,k))*omgbdth(i,k)- &
     1044             alpha_up(i,k+1)*omgbdth(i,k+1))*ppi(i, k)
    10111045          ! print*,'d_deltatw=',d_deltatw(i,k)
    10121046
    10131047          d_deltaqw(i, k) = dtimesub/(ph(i,k)-ph(i,k+1))* &
    1014             (rrd1*omg(i,k)*sigmaw(i)*d_q1(i,k)-rrd2*omg(i,k+1)*(1.-sigmaw( &
    1015             i))*d_q2(i,k+1)-(1.-alpha_up(i,k))*omgbdq(i,k)-alpha_up(i,k+1)* &
    1016             omgbdq(i,k+1))
     1048            (rrd1*omg(i,k)*sigmaw(i)*d_q1(i,k)- &
     1049             rrd2*omg(i,k+1)*(1.-sigmaw(i))*d_q2(i,k+1)- &
     1050             (1.-alpha_up(i,k))*omgbdq(i,k)- &
     1051             alpha_up(i,k+1)*omgbdq(i,k+1))
    10171052          ! print*,'d_deltaqw=',d_deltaqw(i,k)
    10181053
    … …  
    10241059          ! C
    10251060          ! -----------------------------------------------------------------
    1026           d_te(i, k) = dtimesub*((rre1(i)*omg(i,k)*sigmaw(i)*d_th1(i, &
    1027             k)-rre2(i)*omg(i,k+1)*(1.-sigmaw(i))*d_th2(i,k+1))/(ph(i,k)-ph(i, &
    1028             k+1)) &                ! cc nrlmd     $
    1029                                    ! -sigmaw(i)*(1.-sigmaw(i))*dth(i,k)*dp_deltomg(i,k)
    1030             -sigmaw(i)*(1.-sigmaw(i))*dth(i,k)*(omg(i,k)-omg(i,k+1))/(ph(i, &
    1031             k)-ph(i,k+1)) &        ! cc
    1032             )*ppi(i, k)
    1033 
    1034           d_qe(i, k) = dtimesub*((rre1(i)*omg(i,k)*sigmaw(i)*d_q1(i, &
    1035             k)-rre2(i)*omg(i,k+1)*(1.-sigmaw(i))*d_q2(i,k+1))/(ph(i,k)-ph(i, &
    1036             k+1)) &                ! cc nrlmd     $
    1037                                    ! -sigmaw(i)*(1.-sigmaw(i))*deltaqw(i,k)*dp_deltomg(i,k)
    1038             -sigmaw(i)*(1.-sigmaw(i))*deltaqw(i,k)*(omg(i,k)-omg(i, &
    1039             k+1))/(ph(i,k)-ph(i,k+1)) & ! cc
    1040             )
    1041           ! cc nrlmd
     1061          d_te(i, k) = dtimesub*((rre1(i)*omg(i,k)*sigmaw(i)*d_th1(i,k)- &
     1062                                  rre2(i)*omg(i,k+1)*(1.-sigmaw(i))*d_th2(i,k+1))/ &
     1063                                 (ph(i,k)-ph(i,k+1)) &
     1064                                 -sigmaw(i)*(1.-sigmaw(i))*dth(i,k)*(omg(i,k)-omg(i,k+1))/ &
     1065                                 (ph(i,k)-ph(i,k+1)) )*ppi(i, k)
     1066
     1067          d_qe(i, k) = dtimesub*((rre1(i)*omg(i,k)*sigmaw(i)*d_q1(i,k)- &
     1068                                  rre2(i)*omg(i,k+1)*(1.-sigmaw(i))*d_q2(i,k+1))/ &
     1069                                 (ph(i,k)-ph(i,k+1)) &
     1070                                 -sigmaw(i)*(1.-sigmaw(i))*deltaqw(i,k)*(omg(i,k)-omg(i,k+1))/ &
     1071                                 (ph(i,k)-ph(i,k+1)) )
    10421072        ELSE IF (wk_adv(i) .AND. k==kupper(i)) THEN
    1043           d_te(i, k) = dtimesub*((rre1(i)*omg(i,k)*sigmaw(i)*d_th1(i, &
    1044             k)/(ph(i,k)-ph(i,k+1))))*ppi(i, k)
    1045 
    1046           d_qe(i, k) = dtimesub*((rre1(i)*omg(i,k)*sigmaw(i)*d_q1(i, &
    1047             k)/(ph(i,k)-ph(i,k+1))))
     1073          d_te(i, k) = dtimesub*(rre1(i)*omg(i,k)*sigmaw(i)*d_th1(i,k)/(ph(i,k)-ph(i,k+1)))*ppi(i, k)
     1074
     1075          d_qe(i, k) = dtimesub*(rre1(i)*omg(i,k)*sigmaw(i)*d_q1(i,k)/(ph(i,k)-ph(i,k+1)))
    10481076
    10491077        END IF
    … …  
    10671095          crep(i, k) = crep_sol*(ph(i,kupper(i))-ph(i,k))/ &
    10681096            (ph(i,kupper(i))-ph(i,1))
    1069           crep(i, k) = crep(i, k) + crep_upper*(ph(i,1)-ph(i,k))/(p(i,1)-ph(i &
    1070             ,kupper(i)))
     1097          crep(i, k) = crep(i, k) + crep_upper*(ph(i,1)-ph(i,k))/ &
     1098            (p(i,1)-ph(i,kupper(i)))
    10711099
    10721100
    … …  
    10911119          ! print*,'dtKE= ',dtKE(i,k),' dqKE= ',dqKE(i,k)
    10921120
    1093 !jyg<
    1094 !!
    1095 !!---------------------------------------------------------------
    1096 !! The change of delta_T due to PBL (vertical diffusion plus thermal plumes)
    1097 !! is accounted for by the PBL and the Thermals schemes. It is now set to zero
    1098 !! within the Wake scheme.
    1099 !!---------------------------------------------------------------
    1100           dtPBL(i,k) = 0.
    1101           dqPBL(i,k) = 0.
    1102 !
    1103 !!          dtPBL(i,k)=wdtPBL(i,k) - udtPBL(i,k)
    1104 !!          dqPBL(i,k)=wdqPBL(i,k) - udqPBL(i,k)
    1105 !
    1106 !!          dtpbl(i, k) = (wdtpbl(i,k)/sigmaw(i)-udtpbl(i,k)/(1.-sigmaw(i)))
    1107 !!          dqpbl(i, k) = (wdqpbl(i,k)/sigmaw(i)-udqpbl(i,k)/(1.-sigmaw(i)))
    1108           ! print*,'dtPBL= ',dtPBL(i,k),' dqPBL= ',dqPBL(i,k)
    1109 !>jyg
    11101121!
    11111122
    … …  
    11451156
    11461157          IF (dtimesub*tgw(i,k)<1.E-10) THEN
    1147             d_deltatw(i, k) = dtimesub*(ff(i)+dtke(i,k)+dtpbl(i,k) & ! cc
    1148                                                                      ! $
    1149                                                                      ! -spread(i,k)*deltatw(i,k)
    1150               -entr(i,k)*deltatw(i,k)/sigmaw(i)-(death_rate(i)*sigmaw( &
    1151               i)+detr(i,k))*deltatw(i,k)/(1.-sigmaw(i)) & ! cc
    1152               -tgw(i,k)*deltatw(i,k))
     1158            d_deltatw(i, k) = dtimesub*(ff(i)+dtke(i,k) - &
     1159               entr(i,k)*deltatw(i,k)/sigmaw(i) - &
     1160               (death_rate(i)*sigmaw(i)+detr(i,k))*deltatw(i,k)/(1.-sigmaw(i)) - & ! cc
     1161               tgw(i,k)*deltatw(i,k) )
    11531162          ELSE
    1154             d_deltatw(i, k) = 1/tgw(i, k)*(1-exp(-dtimesub*tgw(i, &
    1155               k)))*(ff(i)+dtke(i,k)+dtpbl(i,k) & ! cc     $
    1156                                                  ! -spread(i,k)*deltatw(i,k)
    1157               -entr(i,k)*deltatw(i,k)/sigmaw(i)-(death_rate(i)*sigmaw( &
    1158               i)+detr(i,k))*deltatw(i,k)/(1.-sigmaw(i)) & ! cc
    1159               -tgw(i,k)*deltatw(i,k))
     1163            d_deltatw(i, k) = 1/tgw(i, k)*(1-exp(-dtimesub*tgw(i,k)))* &
     1164               (ff(i)+dtke(i,k) - &
     1165                entr(i,k)*deltatw(i,k)/sigmaw(i) - &
     1166                (death_rate(i)*sigmaw(i)+detr(i,k))*deltatw(i,k)/(1.-sigmaw(i)) - &
     1167                tgw(i,k)*deltatw(i,k) )
    11601168          END IF
    11611169
    … …  
    11641172          gg(i) = d_deltaqw(i, k)/dtimesub
    11651173
    1166           d_deltaqw(i, k) = dtimesub*(gg(i)+dqke(i,k)+dqpbl(i,k) & ! cc     $
    1167                                                                    ! -spread(i,k)*deltaqw(i,k))
    1168             -entr(i,k)*deltaqw(i,k)/sigmaw(i)-(death_rate(i)*sigmaw(i)+detr( &
    1169             i,k))*deltaqw(i,k)/(1.-sigmaw(i)))
     1174          d_deltaqw(i, k) = dtimesub*(gg(i)+dqke(i,k) - &
     1175            entr(i,k)*deltaqw(i,k)/sigmaw(i) - &
     1176            (death_rate(i)*sigmaw(i)+detr(i,k))*deltaqw(i,k)/(1.-sigmaw(i)))
    11701177          ! cc
    11711178
    … …  
    12391246      IF (wk_adv(i)) THEN
    12401247        sigmaw(i) = sigmaw(i) + d_sigmaw(i)
     1248!jyg<
     1249        d_sigmaw2(i) = d_sigmaw2(i) + d_sigmaw(i)
     1250!>jyg
    12411251      END IF
    12421252    END DO
    … …  
    12581268        IF (wk_adv(i) .AND. ptop_provis(i)==ph(i,1) .AND. &
    12591269            dth(i,k)>-delta_t_min .AND. dth(i,k-1)<-delta_t_min) THEN
    1260           ptop_provis(i) = ((dth(i,k)+delta_t_min)*p(i,k-1)-(dth(i, &
    1261             k-1)+delta_t_min)*p(i,k))/(dth(i,k)-dth(i,k-1))
     1270          ptop_provis(i) = ((dth(i,k)+delta_t_min)*p(i,k-1) - &
     1271                            (dth(i,k-1)+delta_t_min)*p(i,k))/(dth(i,k)-dth(i,k-1))
    12621272        END IF
    12631273      END DO
    … …  
    13311341        IF (wk_adv(i) .AND. k<=ktop(i) .AND. ptop_new(i)==ptop(i) .AND. &
    13321342            dth(i,k)>-delta_t_min .AND. dth(i,k-1)<-delta_t_min) THEN
    1333           ptop_new(i) = ((dth(i,k)+delta_t_min)*p(i,k-1)-(dth(i, &
    1334             k-1)+delta_t_min)*p(i,k))/(dth(i,k)-dth(i,k-1))
     1343          ptop_new(i) = ((dth(i,k)+delta_t_min)*p(i,k-1) - &
     1344                         (dth(i,k-1)+delta_t_min)*p(i,k))/(dth(i,k)-dth(i,k-1))
    13351345        END IF
    13361346      END DO
    … …  
    13591369          deltatw(i, k) = 0.
    13601370          deltaqw(i, k) = 0.
     1371          d_deltatw2(i,k) = -deltatw0(i,k)
     1372          d_deltaqw2(i,k) = -deltaqw0(i,k)
    13611373        END IF
    13621374      END DO
    … …  
    14491461        av_dqdwn(i) = sum_dqdwn(i)/hw0(i)
    14501462
    1451         wape(i) = -rg*hw0(i)*(av_dth(i)+epsim1*(av_thu(i)*av_dq(i)+av_dth(i)* &
    1452           av_qu(i)+av_dth(i)*av_dq(i)))/av_thvu(i)
     1463        wape(i) = -rg*hw0(i)*(av_dth(i)+epsim1*(av_thu(i)*av_dq(i) + &
     1464                              av_dth(i)*av_qu(i)+av_dth(i)*av_dq(i)))/av_thvu(i)
    14531465      END IF
    14541466    END DO
    … …  
    14631475            deltaqw(i, k) = 0.
    14641476            dth(i, k) = 0.
     1477            d_deltatw2(i,k) = -deltatw0(i,k)
     1478            d_deltaqw2(i,k) = -deltaqw0(i,k)
    14651479          END IF
    14661480        END IF
    … …  
    14741488          cstar(i) = 0.
    14751489          hw(i) = hwmin
    1476           sigmaw(i) = max(sigmad, sigd_con(i))
     1490!jyg<
     1491!!          sigmaw(i) = max(sigmad, sigd_con(i))
     1492          sigmaw_targ = max(sigmad, sigd_con(i))
     1493          d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i)
     1494          sigmaw(i) = sigmaw_targ
     1495!>jyg
    14771496          fip(i) = 0.
    14781497          gwake(i) = .FALSE.
    … …  
    14861505  END DO ! end sub-timestep loop
    14871506
    1488   ! -----------------------------------------------------------------
    1489   ! Get back to tendencies per second
    1490 
    1491   DO k = 1, klev
    1492     DO i = 1, klon
    1493 
    1494       ! cc nrlmd        IF ( wk_adv(i) .AND. k .LE. kupper(i)) THEN
    1495       IF (ok_qx_qw(i) .AND. k<=kupper(i)) THEN
    1496         ! cc
    1497         dtls(i, k) = dtls(i, k)/dtime
    1498         dqls(i, k) = dqls(i, k)/dtime
    1499         d_deltatw2(i, k) = d_deltatw2(i, k)/dtime
    1500         d_deltaqw2(i, k) = d_deltaqw2(i, k)/dtime
    1501         d_deltat_gw(i, k) = d_deltat_gw(i, k)/dtime
    1502         ! c      print*,'k,dqls,omg,entr,detr',k,dqls(i,k),omg(i,k),entr(i,k)
    1503         ! c     $         ,death_rate(i)*sigmaw(i)
    1504       END IF
    1505     END DO
    1506   END DO
    15071507
    15081508
    … …  
    16321632      av_dqdwn(i) = sum_dqdwn(i)/hw0(i)
    16331633
    1634       wape2(i) = -rg*hw0(i)*(av_dth(i)+epsim1*(av_thu(i)*av_dq(i)+av_dth(i)* &
    1635         av_qu(i)+av_dth(i)*av_dq(i)))/av_thvu(i)
     1634      wape2(i) = -rg*hw0(i)*(av_dth(i)+epsim1*(av_thu(i)*av_dq(i) + &
     1635                             av_dth(i)*av_qu(i)+av_dth(i)*av_dq(i)))/av_thvu(i)
    16361636    END IF
    16371637  END DO
     1638
     1639
    16381640
    16391641  ! Prognostic variable update
    … …  
    16501652        deltaqw(i, k) = 0.
    16511653        dth(i, k) = 0.
     1654        d_deltatw2(i,k) = -deltatw0(i,k)
     1655        d_deltaqw2(i,k) = -deltaqw0(i,k)
    16521656      END IF
    16531657    END DO
    … …  
    16631667        cstar2(i) = 0.
    16641668        hw(i) = hwmin
    1665         sigmaw(i) = amax1(sigmad, sigd_con(i))
     1669!jyg<
     1670!!      sigmaw(i) = amax1(sigmad, sigd_con(i))
     1671      sigmaw_targ = max(sigmad, sigd_con(i))
     1672      d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i)
     1673      sigmaw(i) = sigmaw_targ
     1674!>jyg
    16661675        fip(i) = 0.
    16671676        gwake(i) = .FALSE.
    … …  
    17261735        deltatw(i, k) = 0.
    17271736        deltaqw(i, k) = 0.
     1737        d_deltatw2(i,k) = -deltatw0(i,k)
     1738        d_deltaqw2(i,k) = -deltaqw0(i,k)
    17281739      END IF
    17291740    END DO
    … …  
    17341745    IF (((wape(i)>=wape2(i)) .AND. (wape2(i)<=1.0)) .OR. (ktopw(i)<=2) .OR. &
    17351746        .NOT. ok_qx_qw(i)) THEN
     1747      ktopw(i) = 0
    17361748      wape(i) = 0.
    17371749      cstar(i) = 0.
    … …  
    17491761    ! c     $          wape(i),wape2(i),ktopw(i),OK_qx_qw(i)
    17501762  END DO
     1763
     1764  ! -----------------------------------------------------------------
     1765  ! Get back to tendencies per second
     1766
     1767  DO k = 1, klev
     1768    DO i = 1, klon
     1769
     1770      ! cc nrlmd        IF ( wk_adv(i) .AND. k .LE. kupper(i)) THEN
     1771      IF (ok_qx_qw(i) .AND. k<=kupper(i)) THEN
     1772        ! cc
     1773        dtls(i, k) = dtls(i, k)/dtime
     1774        dqls(i, k) = dqls(i, k)/dtime
     1775        d_deltatw2(i, k) = d_deltatw2(i, k)/dtime
     1776        d_deltaqw2(i, k) = d_deltaqw2(i, k)/dtime
     1777        d_deltat_gw(i, k) = d_deltat_gw(i, k)/dtime
     1778        ! c      print*,'k,dqls,omg,entr,detr',k,dqls(i,k),omg(i,k),entr(i,k)
     1779        ! c     $         ,death_rate(i)*sigmaw(i)
     1780      END IF
     1781    END DO
     1782  END DO
     1783!jyg<
     1784  DO i = 1, klon
     1785    d_sigmaw2(i) = d_sigmaw2(i)/dtime
     1786    d_wdens2(i) = d_wdens2(i)/dtime
     1787  ENDDO
     1788!>jyg
     1789
    17511790
    17521791
    … …  
    17921831      IF (wk_adv(i)) THEN
    17931832        x = qe(i, k) + (zeta(i,k)-sigmaw(i))*deltaqw(i, k) + d_qe(i, k) + &
    1794           (zeta(i,k)-sigmaw(i))*d_deltaqw(i, k) - d_sigmaw(i)*(deltaqw(i,k)+ &
    1795           d_deltaqw(i,k))
     1833          (zeta(i,k)-sigmaw(i))*d_deltaqw(i, k) - d_sigmaw(i) * &
     1834          (deltaqw(i,k)+d_deltaqw(i,k))
    17961835        a = -d_sigmaw(i)*d_deltaqw(i, k)
    17971836        b = d_qe(i, k) + (zeta(i,k)-sigmaw(i))*d_deltaqw(i, k) - &
    … …  
    18071846          ELSE
    18081847            IF (a>0.) THEN
    1809               alpha1(i) = 0.9*min((2.*c)/(-b+sqrt(discrim)), (-b+sqrt(discrim &
    1810                 ))/(2.*a))
     1848              alpha1(i) = 0.9*min( (2.*c)/(-b+sqrt(discrim)), &
     1849                                   (-b+sqrt(discrim))/(2.*a) )
    18111850            ELSE IF (a==0.) THEN
    18121851              alpha1(i) = 0.9*(-c/b)
    18131852            ELSE
    18141853              ! print*,'a,b,c discrim',a,b,c discrim
    1815               alpha1(i) = 0.9*max((2.*c)/(-b+sqrt(discrim)), (-b+sqrt(discrim &
    1816                 ))/(2.*a))
     1854              alpha1(i) = 0.9*max( (2.*c)/(-b+sqrt(discrim)), &
     1855                                   (-b+sqrt(discrim))/(2.*a))
    18171856            END IF
    18181857          END IF
    … …  
    18271866
    18281867
     1868
     1869
Note: See TracChangeset for help on using the changeset viewer.