Ignore:
Timestamp:
Jul 23, 2024, 8:38:56 AM (5 months ago)
Author:
abarral
Message:

Handle CPP_INLANDSIS in lmdz_cppkeys_wrapper.F90
Remove obsolete key wrgrads_thermcell

Location:
LMDZ6/branches/Amaury_dev/libf
Files:
2 edited
2 moved

Legend:

Unmodified
Added
Removed

  • LMDZ6/branches/Amaury_dev/libf/misc/lmdz_cppkeys_wrapper.F90

    r5101 r5102  
    55! and imported through USE ..., ONLY: ... elsewhere
    66! CPP keys supported (key -> fortran variables associated):
    7 !      NC_DOUBLE -> nf90_format
    8 !      CPP_PHYS  -> CPPKEY_PHYS
    9 !      INCA      -> CPPKEY_INCA
    10 !      CPP_StratAer-> CPPKEY_STRATAER
    11 !      CPP_DUST  -> CPPKEY_DUST
     7!      NC_DOUBLE      -> nf90_format
     8!      CPP_PHYS       -> CPPKEY_PHYS
     9!      INCA           -> CPPKEY_INCA
     10!      CPP_StratAer   -> CPPKEY_STRATAER
     11!      CPP_DUST       -> CPPKEY_DUST
     12!      CPP_INLANDSIS  -> CPPKEY_INLANDSIS
    1213! ---------------------------------------------
    1314
    … …  
    1718  IMPLICIT NONE; PRIVATE
    1819  PUBLIC nf90_format, CPPKEY_PHYS, CPPKEY_INCA, CPPKEY_STRATAER, CPPKEY_DUST, &
    19       CPPKEY_DEBUGIO
     20      CPPKEY_DEBUGIO, CPPKEY_INLANDSIS
    2021
    2122#ifdef NC_DOUBLE
    … …  
    5556#endif
    5657
     58#ifdef CPP_INLANDSIS
     59  LOGICAL, PARAMETER :: CPPKEY_INLANDSIS = .TRUE.
     60#else
     61  LOGICAL, PARAMETER :: CPPKEY_INLANDSIS = .FALSE.
     62#endif
     63
    5764END MODULE lmdz_cppkeys_wrapper

  • LMDZ6/branches/Amaury_dev/libf/phylmd/lmdz_thermcell_plume.f90

    r5101 r5102  
    11MODULE lmdz_thermcell_plume
    22
    3 ! $Id: thermcell_plume.F90 3074 2017-11-15 13:31:44Z fhourdin $
     3  ! $Id: thermcell_plume.F90 3074 2017-11-15 13:31:44Z fhourdin $
    44
    55CONTAINS
    66
    7       SUBROUTINE thermcell_plume(itap,ngrid,nlay,ptimestep,ztv,zthl,po,zl,rhobarz,  &
    8              zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot,  &
    9              lalim,f0,detr_star,entr_star,f_star,csc,ztva,  &
    10              ztla,zqla,zqta,zha,zw2,w_est,ztva_est,zqsatth,lmix,lmix_bis,linter &
    11              ,lev_out,lunout1,igout)
    12 !     &           ,lev_out,lunout1,igout,zbuoy,zbuoyjam)
    13 !--------------------------------------------------------------------------
    14 ! Auhtors : Catherine Rio, Frédéric Hourdin, Arnaud Jam
    15 
    16 !thermcell_plume: calcule les valeurs de qt, thetal et w dans l ascendance
    17 !   This versions starts from a cleaning of thermcell_plume_6A (2019/01/20)
    18 !   thermcell_plume_6A is activate for flag_thermas_ed < 10
    19 !   thermcell_plume_5B for flag_thermas_ed < 20
    20 !   thermcell_plume for flag_thermals_ed>= 20
    21 !   Various options are controled by the flag_thermals_ed parameter
    22 !   = 20 : equivalent to thermcell_plume_6A with flag_thermals_ed=8
    23 !   = 21 : the Jam strato-cumulus modif is not activated in detrainment
    24 !   = 29 : an other way to compute the modified buoyancy (to be tested)
    25 !--------------------------------------------------------------------------
    26 
    27        USE lmdz_thermcell_ini, ONLY: prt_level,fact_thermals_ed_dz,iflag_thermals_ed,RLvCP,RETV,RG
    28        USE lmdz_thermcell_ini, ONLY: fact_epsilon, betalpha, afact, fact_shell
    29        USE lmdz_thermcell_ini, ONLY: detr_min, entr_min, detr_q_coef, detr_q_power
    30        USE lmdz_thermcell_ini, ONLY: mix0, thermals_flag_alim
    31        USE lmdz_thermcell_alim, ONLY: thermcell_alim
    32        USE lmdz_thermcell_qsat, ONLY: thermcell_qsat
    33 
    34 
    35        IMPLICIT NONE
    36 
    37       integer,intent(in) :: itap,lev_out,lunout1,igout,ngrid,nlay
    38       real,intent(in) :: ptimestep
    39       real,intent(in),dimension(ngrid,nlay) :: ztv
    40       real,intent(in),dimension(ngrid,nlay) :: zthl
    41       real,intent(in),dimension(ngrid,nlay) :: po
    42       real,intent(in),dimension(ngrid,nlay) :: zl
    43       real,intent(in),dimension(ngrid,nlay) :: rhobarz
    44       real,intent(in),dimension(ngrid,nlay+1) :: zlev
    45       real,intent(in),dimension(ngrid,nlay+1) :: pplev
    46       real,intent(in),dimension(ngrid,nlay) :: pphi
    47       real,intent(in),dimension(ngrid,nlay) :: zpspsk
    48       real,intent(in),dimension(ngrid) :: f0
    49 
    50       integer,intent(out) :: lalim(ngrid)
    51       real,intent(out),dimension(ngrid,nlay) :: alim_star
    52       real,intent(out),dimension(ngrid) :: alim_star_tot
    53       real,intent(out),dimension(ngrid,nlay) :: detr_star
    54       real,intent(out),dimension(ngrid,nlay) :: entr_star
    55       real,intent(out),dimension(ngrid,nlay+1) :: f_star
    56       real,intent(out),dimension(ngrid,nlay) :: csc
    57       real,intent(out),dimension(ngrid,nlay) :: ztva
    58       real,intent(out),dimension(ngrid,nlay) :: ztla
    59       real,intent(out),dimension(ngrid,nlay) :: zqla
    60       real,intent(out),dimension(ngrid,nlay) :: zqta
    61       real,intent(out),dimension(ngrid,nlay) :: zha
    62       real,intent(out),dimension(ngrid,nlay+1) :: zw2
    63       real,intent(out),dimension(ngrid,nlay+1) :: w_est
    64       real,intent(out),dimension(ngrid,nlay) :: ztva_est
    65       real,intent(out),dimension(ngrid,nlay) :: zqsatth
    66       integer,intent(out),dimension(ngrid) :: lmix(ngrid)
    67       integer,intent(out),dimension(ngrid) :: lmix_bis(ngrid)
    68       real,intent(out),dimension(ngrid) :: linter(ngrid)
    69 
    70 
    71       REAL,dimension(ngrid,nlay+1) :: wa_moy
    72       REAL,dimension(ngrid,nlay) :: entr,detr
    73       REAL,dimension(ngrid,nlay) :: ztv_est
    74       REAL,dimension(ngrid,nlay) :: zqla_est
    75       REAL,dimension(ngrid,nlay) :: zta_est
    76       REAL,dimension(ngrid) :: ztemp,zqsat
    77       REAL zdw2,zdw2bis
    78       REAL zw2modif
    79       REAL zw2fact,zw2factbis
    80       REAL,dimension(ngrid,nlay) :: zeps
    81 
    82       REAL,dimension(ngrid) :: wmaxa
    83 
    84       INTEGER ig,l,k,lt,it,lm,nbpb
    85 
    86       real,dimension(ngrid,nlay) :: zbuoy,gamma,zdqt
    87       real zdz,zalpha,zw2m
    88       real,dimension(ngrid,nlay) :: zbuoyjam,zdqtjam
    89       real zdz2,zdz3,lmel,entrbis,zdzbis
    90       real,dimension(ngrid) :: d_temp
    91       real ztv1,ztv2,factinv,zinv,zlmel
    92       real zlmelup,zlmeldwn,zlt,zltdwn,zltup
    93       real atv1,atv2,btv1,btv2
    94       real ztv_est1,ztv_est2
    95       real zcor,zdelta,zcvm5,qlbef
    96       real zbetalpha, coefzlmel
    97       real eps
    98       logical Zsat
    99       LOGICAL,dimension(ngrid) :: active,activetmp
    100       REAL fact_gamma,fact_gamma2,fact_epsilon2
    101 
    102 
    103       REAL,dimension(ngrid,nlay) :: c2
    104 
    105       if (ngrid==1) print*,'THERMCELL PLUME MODIFIE 2014/07/11'
    106       Zsat=.false.
    107 ! Initialisation
    108 
    109 
    110       zbetalpha=betalpha/(1.+betalpha)
    111 
    112 
    113 ! Initialisations des variables r?elles
    114 if (1==1) then
    115       ztva(:,:)=ztv(:,:)
    116       ztva_est(:,:)=ztva(:,:)
    117       ztv_est(:,:)=ztv(:,:)
    118       ztla(:,:)=zthl(:,:)
    119       zqta(:,:)=po(:,:)
    120       zqla(:,:)=0.
    121       zha(:,:) = ztva(:,:)
    122 else
    123       ztva(:,:)=0.
    124       ztv_est(:,:)=0.
    125       ztva_est(:,:)=0.
    126       ztla(:,:)=0.
    127       zqta(:,:)=0.
    128       zha(:,:) =0.
    129 endif
    130 
    131       zqla_est(:,:)=0.
    132       zqsatth(:,:)=0.
    133       zqla(:,:)=0.
    134       detr_star(:,:)=0.
    135       entr_star(:,:)=0.
    136       alim_star(:,:)=0.
    137       alim_star_tot(:)=0.
    138       csc(:,:)=0.
    139       detr(:,:)=0.
    140       entr(:,:)=0.
    141       zw2(:,:)=0.
    142       zbuoy(:,:)=0.
    143       zbuoyjam(:,:)=0.
    144       gamma(:,:)=0.
    145       zeps(:,:)=0.
    146       w_est(:,:)=0.
    147       f_star(:,:)=0.
    148       wa_moy(:,:)=0.
    149       linter(:)=1.
    150 !     linter(:)=1.
    151 ! Initialisation des variables entieres
    152       lmix(:)=1
    153       lmix_bis(:)=2
    154       wmaxa(:)=0.
    155 
    156 
    157 !-------------------------------------------------------------------------
    158 ! On ne considere comme actif que les colonnes dont les deux premieres
    159 ! couches sont instables.
    160 !-------------------------------------------------------------------------
    161 
    162       active(:)=ztv(:,1)>ztv(:,2)
    163       d_temp(:)=0. ! Pour activer un contraste de temperature a la base
    164                    ! du panache
    165 !  Cet appel pourrait être fait avant thermcell_plume dans thermcell_main
    166       CALL thermcell_alim(thermals_flag_alim,ngrid,nlay,ztv,d_temp,zlev,alim_star,lalim)
    167 
    168 !------------------------------------------------------------------------------
    169 ! Calcul dans la premiere couche
    170 ! On decide dans cette version que le thermique n'est actif que si la premiere
    171 ! couche est instable.
    172 ! Pourrait etre change si on veut que le thermiques puisse se d??clencher
    173 ! dans une couche l>1
    174 !------------------------------------------------------------------------------
    175 do ig=1,ngrid
    176 ! Le panache va prendre au debut les caracteristiques de l'air contenu
    177 ! dans cette couche.
    178     if (active(ig)) then
    179     ztla(ig,1)=zthl(ig,1)
    180     zqta(ig,1)=po(ig,1)
    181     zqla(ig,1)=zl(ig,1)
    182 !cr: attention, prise en compte de f*(1)=1
    183     f_star(ig,2)=alim_star(ig,1)
    184     zw2(ig,2)=2.*RG*(ztv(ig,1)-ztv(ig,2))/ztv(ig,2)  &
    185 &                     *(zlev(ig,2)-zlev(ig,1))  &
    186 &                     *0.4*pphi(ig,1)/(pphi(ig,2)-pphi(ig,1))
    187     w_est(ig,2)=zw2(ig,2)
     7  SUBROUTINE thermcell_plume(itap, ngrid, nlay, ptimestep, ztv, zthl, po, zl, rhobarz, &
     8          zlev, pplev, pphi, zpspsk, alim_star, alim_star_tot, &
     9          lalim, f0, detr_star, entr_star, f_star, csc, ztva, &
     10          ztla, zqla, zqta, zha, zw2, w_est, ztva_est, zqsatth, lmix, lmix_bis, linter &
     11          , lev_out, lunout1, igout)
     12    !     &           ,lev_out,lunout1,igout,zbuoy,zbuoyjam)
     13    !--------------------------------------------------------------------------
     14    ! Auhtors : Catherine Rio, Frédéric Hourdin, Arnaud Jam
     15
     16    !thermcell_plume: calcule les valeurs de qt, thetal et w dans l ascendance
     17    !   This versions starts from a cleaning of thermcell_plume_6A (2019/01/20)
     18    !   thermcell_plume_6A is activate for flag_thermas_ed < 10
     19    !   thermcell_plume_5B for flag_thermas_ed < 20
     20    !   thermcell_plume for flag_thermals_ed>= 20
     21    !   Various options are controled by the flag_thermals_ed parameter
     22    !   = 20 : equivalent to thermcell_plume_6A with flag_thermals_ed=8
     23    !   = 21 : the Jam strato-cumulus modif is not activated in detrainment
     24    !   = 29 : an other way to compute the modified buoyancy (to be tested)
     25    !--------------------------------------------------------------------------
     26
     27    USE lmdz_thermcell_ini, ONLY: prt_level, fact_thermals_ed_dz, iflag_thermals_ed, RLvCP, RETV, RG
     28    USE lmdz_thermcell_ini, ONLY: fact_epsilon, betalpha, afact, fact_shell
     29    USE lmdz_thermcell_ini, ONLY: detr_min, entr_min, detr_q_coef, detr_q_power
     30    USE lmdz_thermcell_ini, ONLY: mix0, thermals_flag_alim
     31    USE lmdz_thermcell_alim, ONLY: thermcell_alim
     32    USE lmdz_thermcell_qsat, ONLY: thermcell_qsat
     33
     34    IMPLICIT NONE
     35
     36    integer, intent(in) :: itap, lev_out, lunout1, igout, ngrid, nlay
     37    real, intent(in) :: ptimestep
     38    real, intent(in), dimension(ngrid, nlay) :: ztv
     39    real, intent(in), dimension(ngrid, nlay) :: zthl
     40    real, intent(in), dimension(ngrid, nlay) :: po
     41    real, intent(in), dimension(ngrid, nlay) :: zl
     42    real, intent(in), dimension(ngrid, nlay) :: rhobarz
     43    real, intent(in), dimension(ngrid, nlay + 1) :: zlev
     44    real, intent(in), dimension(ngrid, nlay + 1) :: pplev
     45    real, intent(in), dimension(ngrid, nlay) :: pphi
     46    real, intent(in), dimension(ngrid, nlay) :: zpspsk
     47    real, intent(in), dimension(ngrid) :: f0
     48
     49    integer, intent(out) :: lalim(ngrid)
     50    real, intent(out), dimension(ngrid, nlay) :: alim_star
     51    real, intent(out), dimension(ngrid) :: alim_star_tot
     52    real, intent(out), dimension(ngrid, nlay) :: detr_star
     53    real, intent(out), dimension(ngrid, nlay) :: entr_star
     54    real, intent(out), dimension(ngrid, nlay + 1) :: f_star
     55    real, intent(out), dimension(ngrid, nlay) :: csc
     56    real, intent(out), dimension(ngrid, nlay) :: ztva
     57    real, intent(out), dimension(ngrid, nlay) :: ztla
     58    real, intent(out), dimension(ngrid, nlay) :: zqla
     59    real, intent(out), dimension(ngrid, nlay) :: zqta
     60    real, intent(out), dimension(ngrid, nlay) :: zha
     61    real, intent(out), dimension(ngrid, nlay + 1) :: zw2
     62    real, intent(out), dimension(ngrid, nlay + 1) :: w_est
     63    real, intent(out), dimension(ngrid, nlay) :: ztva_est
     64    real, intent(out), dimension(ngrid, nlay) :: zqsatth
     65    integer, intent(out), dimension(ngrid) :: lmix(ngrid)
     66    integer, intent(out), dimension(ngrid) :: lmix_bis(ngrid)
     67    real, intent(out), dimension(ngrid) :: linter(ngrid)
     68
     69    REAL, dimension(ngrid, nlay + 1) :: wa_moy
     70    REAL, dimension(ngrid, nlay) :: entr, detr
     71    REAL, dimension(ngrid, nlay) :: ztv_est
     72    REAL, dimension(ngrid, nlay) :: zqla_est
     73    REAL, dimension(ngrid, nlay) :: zta_est
     74    REAL, dimension(ngrid) :: ztemp, zqsat
     75    REAL zdw2, zdw2bis
     76    REAL zw2modif
     77    REAL zw2fact, zw2factbis
     78    REAL, dimension(ngrid, nlay) :: zeps
     79
     80    REAL, dimension(ngrid) :: wmaxa
     81
     82    INTEGER ig, l, k, lt, it, lm, nbpb
     83
     84    real, dimension(ngrid, nlay) :: zbuoy, gamma, zdqt
     85    real zdz, zalpha, zw2m
     86    real, dimension(ngrid, nlay) :: zbuoyjam, zdqtjam
     87    real zdz2, zdz3, lmel, entrbis, zdzbis
     88    real, dimension(ngrid) :: d_temp
     89    real ztv1, ztv2, factinv, zinv, zlmel
     90    real zlmelup, zlmeldwn, zlt, zltdwn, zltup
     91    real atv1, atv2, btv1, btv2
     92    real ztv_est1, ztv_est2
     93    real zcor, zdelta, zcvm5, qlbef
     94    real zbetalpha, coefzlmel
     95    real eps
     96    logical Zsat
     97    LOGICAL, dimension(ngrid) :: active, activetmp
     98    REAL fact_gamma, fact_gamma2, fact_epsilon2
     99
     100    REAL, dimension(ngrid, nlay) :: c2
     101
     102    if (ngrid==1) print*, 'THERMCELL PLUME MODIFIE 2014/07/11'
     103    Zsat = .false.
     104    ! Initialisation
     105
     106    zbetalpha = betalpha / (1. + betalpha)
     107
     108
     109    ! Initialisations des variables r?elles
     110    if (1==1) then
     111      ztva(:, :) = ztv(:, :)
     112      ztva_est(:, :) = ztva(:, :)
     113      ztv_est(:, :) = ztv(:, :)
     114      ztla(:, :) = zthl(:, :)
     115      zqta(:, :) = po(:, :)
     116      zqla(:, :) = 0.
     117      zha(:, :) = ztva(:, :)
     118    else
     119      ztva(:, :) = 0.
     120      ztv_est(:, :) = 0.
     121      ztva_est(:, :) = 0.
     122      ztla(:, :) = 0.
     123      zqta(:, :) = 0.
     124      zha(:, :) = 0.
    188125    endif
    189 enddo
    190 
    191 !==============================================================================
    192 !boucle de calcul de la vitesse verticale dans le thermique
    193 !==============================================================================
    194 do l=2,nlay-1
    195 !==============================================================================
    196 
    197 
    198 ! On decide si le thermique est encore actif ou non
    199 ! AFaire : Il faut sans doute ajouter entr_star a alim_star dans ce test
    200     do ig=1,ngrid
    201        active(ig)=active(ig) &
    202 &                 .and. zw2(ig,l)>1.e-10 &
    203 &                 .and. f_star(ig,l)+alim_star(ig,l)>1.e-10
     126
     127    zqla_est(:, :) = 0.
     128    zqsatth(:, :) = 0.
     129    zqla(:, :) = 0.
     130    detr_star(:, :) = 0.
     131    entr_star(:, :) = 0.
     132    alim_star(:, :) = 0.
     133    alim_star_tot(:) = 0.
     134    csc(:, :) = 0.
     135    detr(:, :) = 0.
     136    entr(:, :) = 0.
     137    zw2(:, :) = 0.
     138    zbuoy(:, :) = 0.
     139    zbuoyjam(:, :) = 0.
     140    gamma(:, :) = 0.
     141    zeps(:, :) = 0.
     142    w_est(:, :) = 0.
     143    f_star(:, :) = 0.
     144    wa_moy(:, :) = 0.
     145    linter(:) = 1.
     146    !     linter(:)=1.
     147    ! Initialisation des variables entieres
     148    lmix(:) = 1
     149    lmix_bis(:) = 2
     150    wmaxa(:) = 0.
     151
     152
     153    !-------------------------------------------------------------------------
     154    ! On ne considere comme actif que les colonnes dont les deux premieres
     155    ! couches sont instables.
     156    !-------------------------------------------------------------------------
     157
     158    active(:) = ztv(:, 1)>ztv(:, 2)
     159    d_temp(:) = 0. ! Pour activer un contraste de temperature a la base
     160    ! du panache
     161    !  Cet appel pourrait être fait avant thermcell_plume dans thermcell_main
     162    CALL thermcell_alim(thermals_flag_alim, ngrid, nlay, ztv, d_temp, zlev, alim_star, lalim)
     163
     164    !------------------------------------------------------------------------------
     165    ! Calcul dans la premiere couche
     166    ! On decide dans cette version que le thermique n'est actif que si la premiere
     167    ! couche est instable.
     168    ! Pourrait etre change si on veut que le thermiques puisse se d??clencher
     169    ! dans une couche l>1
     170    !------------------------------------------------------------------------------
     171    do ig = 1, ngrid
     172      ! Le panache va prendre au debut les caracteristiques de l'air contenu
     173      ! dans cette couche.
     174      if (active(ig)) then
     175        ztla(ig, 1) = zthl(ig, 1)
     176        zqta(ig, 1) = po(ig, 1)
     177        zqla(ig, 1) = zl(ig, 1)
     178        !cr: attention, prise en compte de f*(1)=1
     179        f_star(ig, 2) = alim_star(ig, 1)
     180        zw2(ig, 2) = 2. * RG * (ztv(ig, 1) - ztv(ig, 2)) / ztv(ig, 2)  &
     181                & * (zlev(ig, 2) - zlev(ig, 1))  &
     182                & * 0.4 * pphi(ig, 1) / (pphi(ig, 2) - pphi(ig, 1))
     183        w_est(ig, 2) = zw2(ig, 2)
     184      endif
    204185    enddo
    205186
    206 
    207 
    208 !---------------------------------------------------------------------------
    209 ! calcul des proprietes thermodynamiques et de la vitesse de la couche l
    210 ! sans tenir compte du detrainement et de l'entrainement dans cette
    211 ! couche
    212 ! C'est a dire qu'on suppose
    213 ! ztla(l)=ztla(l-1) et zqta(l)=zqta(l-1)
    214 ! Ici encore, on doit pouvoir ajouter entr_star (qui peut etre calculer
    215 ! avant) a l'alimentation pour avoir un calcul plus propre
    216 !---------------------------------------------------------------------------
    217 
    218    ztemp(:)=zpspsk(:,l)*ztla(:,l-1)
    219    call thermcell_qsat(ngrid,active,pplev(:,l),ztemp,zqta(:,l-1),zqsat(:))
    220     do ig=1,ngrid
    221 !       print*,'active',active(ig),ig,l
    222         if(active(ig)) then
    223         zqla_est(ig,l)=max(0.,zqta(ig,l-1)-zqsat(ig))
    224         ztva_est(ig,l) = ztla(ig,l-1)*zpspsk(ig,l)+RLvCp*zqla_est(ig,l)
    225         zta_est(ig,l)=ztva_est(ig,l)
    226         ztva_est(ig,l) = ztva_est(ig,l)/zpspsk(ig,l)
    227         ztva_est(ig,l) = ztva_est(ig,l)*(1.+RETV*(zqta(ig,l-1)  &
    228         -zqla_est(ig,l))-zqla_est(ig,l))
    229  
    230 
    231 !Modif AJAM
    232 
    233         zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
    234         zdz=zlev(ig,l+1)-zlev(ig,l)         
    235         lmel=fact_thermals_ed_dz*zlev(ig,l)
    236 !        lmel=0.09*zlev(ig,l)
    237         zlmel=zlev(ig,l)+lmel
    238         zlmelup=zlmel+(zdz/2)
    239         zlmeldwn=zlmel-(zdz/2)
    240 
    241         lt=l+1
    242         zlt=zlev(ig,lt)
    243         zdz3=zlev(ig,lt+1)-zlt
    244         zltdwn=zlt-zdz3/2
    245         zltup=zlt+zdz3/2
    246          
    247 !=========================================================================
    248 ! 3. Calcul de la flotabilite modifie par melange avec l'air au dessus
    249 !=========================================================================
    250 
    251 !--------------------------------------------------
    252         lt=l+1
    253         zlt=zlev(ig,lt)
    254         zdz2=zlev(ig,lt)-zlev(ig,l)
    255 
    256         do while (lmel>zdz2)
    257            lt=lt+1
    258            zlt=zlev(ig,lt)
    259            zdz2=zlt-zlev(ig,l)
    260         enddo
    261         zdz3=zlev(ig,lt+1)-zlt
    262         zltdwn=zlev(ig,lt)-zdz3/2
    263         zlmelup=zlmel+(zdz/2)
    264         coefzlmel=Min(1.,(zlmelup-zltdwn)/zdz)
    265         zbuoyjam(ig,l)=1.*RG*(coefzlmel*(ztva_est(ig,l)- &
    266      ztv(ig,lt))/ztv(ig,lt)+(1.-coefzlmel)*(ztva_est(ig,l)- &
    267      ztv(ig,lt-1))/ztv(ig,lt-1))+0.*zbuoy(ig,l)
    268 
    269 !------------------------------------------------
    270 !AJAM:nouveau calcul de w? 
    271 !------------------------------------------------
    272         zdz=zlev(ig,l+1)-zlev(ig,l)
    273         zdzbis=zlev(ig,l)-zlev(ig,l-1)
    274         zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
    275         zw2fact=fact_epsilon*2.*zdz/(1.+betalpha)
    276         zw2factbis=fact_epsilon*2.*zdzbis/(1.+betalpha)
    277         zdw2=afact*zbuoy(ig,l)/fact_epsilon
    278         zdw2bis=afact*zbuoy(ig,l-1)/fact_epsilon
    279         lm=Max(1,l-2)
    280         w_est(ig,l+1)=Max(0.0001,exp(-zw2fact)*(w_est(ig,l)-zdw2)+zdw2)
    281        endif
     187    !==============================================================================
     188    !boucle de calcul de la vitesse verticale dans le thermique
     189    !==============================================================================
     190    do l = 2, nlay - 1
     191      !==============================================================================
     192
     193
     194      ! On decide si le thermique est encore actif ou non
     195      ! AFaire : Il faut sans doute ajouter entr_star a alim_star dans ce test
     196      do ig = 1, ngrid
     197        active(ig) = active(ig) &
     198                &                 .and. zw2(ig, l)>1.e-10 &
     199                &                 .and. f_star(ig, l) + alim_star(ig, l)>1.e-10
     200      enddo
     201
     202
     203
     204      !---------------------------------------------------------------------------
     205      ! calcul des proprietes thermodynamiques et de la vitesse de la couche l
     206      ! sans tenir compte du detrainement et de l'entrainement dans cette
     207      ! couche
     208      ! C'est a dire qu'on suppose
     209      ! ztla(l)=ztla(l-1) et zqta(l)=zqta(l-1)
     210      ! Ici encore, on doit pouvoir ajouter entr_star (qui peut etre calculer
     211      ! avant) a l'alimentation pour avoir un calcul plus propre
     212      !---------------------------------------------------------------------------
     213
     214      ztemp(:) = zpspsk(:, l) * ztla(:, l - 1)
     215      call thermcell_qsat(ngrid, active, pplev(:, l), ztemp, zqta(:, l - 1), zqsat(:))
     216      do ig = 1, ngrid
     217        !       print*,'active',active(ig),ig,l
     218        if(active(ig)) then
     219          zqla_est(ig, l) = max(0., zqta(ig, l - 1) - zqsat(ig))
     220          ztva_est(ig, l) = ztla(ig, l - 1) * zpspsk(ig, l) + RLvCp * zqla_est(ig, l)
     221          zta_est(ig, l) = ztva_est(ig, l)
     222          ztva_est(ig, l) = ztva_est(ig, l) / zpspsk(ig, l)
     223          ztva_est(ig, l) = ztva_est(ig, l) * (1. + RETV * (zqta(ig, l - 1)  &
     224                  - zqla_est(ig, l)) - zqla_est(ig, l))
     225
     226
     227          !Modif AJAM
     228
     229          zbuoy(ig, l) = RG * (ztva_est(ig, l) - ztv(ig, l)) / ztv(ig, l)
     230          zdz = zlev(ig, l + 1) - zlev(ig, l)
     231          lmel = fact_thermals_ed_dz * zlev(ig, l)
     232          !        lmel=0.09*zlev(ig,l)
     233          zlmel = zlev(ig, l) + lmel
     234          zlmelup = zlmel + (zdz / 2)
     235          zlmeldwn = zlmel - (zdz / 2)
     236
     237          lt = l + 1
     238          zlt = zlev(ig, lt)
     239          zdz3 = zlev(ig, lt + 1) - zlt
     240          zltdwn = zlt - zdz3 / 2
     241          zltup = zlt + zdz3 / 2
     242
     243          !=========================================================================
     244          ! 3. Calcul de la flotabilite modifie par melange avec l'air au dessus
     245          !=========================================================================
     246
     247          !--------------------------------------------------
     248          lt = l + 1
     249          zlt = zlev(ig, lt)
     250          zdz2 = zlev(ig, lt) - zlev(ig, l)
     251
     252          do while (lmel>zdz2)
     253            lt = lt + 1
     254            zlt = zlev(ig, lt)
     255            zdz2 = zlt - zlev(ig, l)
     256          enddo
     257          zdz3 = zlev(ig, lt + 1) - zlt
     258          zltdwn = zlev(ig, lt) - zdz3 / 2
     259          zlmelup = zlmel + (zdz / 2)
     260          coefzlmel = Min(1., (zlmelup - zltdwn) / zdz)
     261          zbuoyjam(ig, l) = 1. * RG * (coefzlmel * (ztva_est(ig, l) - &
     262                  ztv(ig, lt)) / ztv(ig, lt) + (1. - coefzlmel) * (ztva_est(ig, l) - &
     263                  ztv(ig, lt - 1)) / ztv(ig, lt - 1)) + 0. * zbuoy(ig, l)
     264
     265          !------------------------------------------------
     266          !AJAM:nouveau calcul de w?
     267          !------------------------------------------------
     268          zdz = zlev(ig, l + 1) - zlev(ig, l)
     269          zdzbis = zlev(ig, l) - zlev(ig, l - 1)
     270          zbuoy(ig, l) = RG * (ztva_est(ig, l) - ztv(ig, l)) / ztv(ig, l)
     271          zw2fact = fact_epsilon * 2. * zdz / (1. + betalpha)
     272          zw2factbis = fact_epsilon * 2. * zdzbis / (1. + betalpha)
     273          zdw2 = afact * zbuoy(ig, l) / fact_epsilon
     274          zdw2bis = afact * zbuoy(ig, l - 1) / fact_epsilon
     275          lm = Max(1, l - 2)
     276          w_est(ig, l + 1) = Max(0.0001, exp(-zw2fact) * (w_est(ig, l) - zdw2) + zdw2)
     277        endif
     278      enddo
     279
     280
     281      !-------------------------------------------------
     282      !calcul des taux d'entrainement et de detrainement
     283      !-------------------------------------------------
     284
     285      do ig = 1, ngrid
     286        if (active(ig)) then
     287
     288          !          zw2m=max(0.5*(w_est(ig,l)+w_est(ig,l+1)),0.1)
     289          zw2m = w_est(ig, l + 1)
     290          zdz = zlev(ig, l + 1) - zlev(ig, l)
     291          zbuoy(ig, l) = RG * (ztva_est(ig, l) - ztv(ig, l)) / ztv(ig, l)
     292          zalpha = f0(ig) * f_star(ig, l) / sqrt(w_est(ig, l + 1)) / rhobarz(ig, l)
     293          zdqt(ig, l) = max(zqta(ig, l - 1) - po(ig, l), 0.) / po(ig, l)
     294
     295          !=========================================================================
     296          ! 4. Calcul de l'entrainement et du detrainement
     297          !=========================================================================
     298
     299          detr_star(ig, l) = f_star(ig, l) * zdz             &
     300                  * (mix0 * 0.1 / (zalpha + 0.001)               &
     301                          + MAX(detr_min, -afact * zbetalpha * zbuoyjam(ig, l) / zw2m   &
     302                                  + detr_q_coef * (zdqt(ig, l) / zw2m)**detr_q_power))
     303
     304          if (iflag_thermals_ed == 20) then
     305            entr_star(ig, l) = f_star(ig, l) * zdz * (&
     306                    mix0 * 0.1 / (zalpha + 0.001)               &
     307                            + zbetalpha * MAX(entr_min, &
     308                            afact * zbuoyjam(ig, l) / zw2m - fact_epsilon))
     309          else
     310            entr_star(ig, l) = f_star(ig, l) * zdz * (&
     311                    mix0 * 0.1 / (zalpha + 0.001)               &
     312                            + zbetalpha * MAX(entr_min, &
     313                            afact * zbuoy(ig, l) / zw2m - fact_epsilon))
     314          endif
     315
     316          ! En dessous de lalim, on prend le max de alim_star et entr_star pour
     317          ! alim_star et 0 sinon
     318          if (l<lalim(ig)) then
     319            alim_star(ig, l) = max(alim_star(ig, l), entr_star(ig, l))
     320            entr_star(ig, l) = 0.
     321          endif
     322          f_star(ig, l + 1) = f_star(ig, l) + alim_star(ig, l) + entr_star(ig, l)  &
     323                  - detr_star(ig, l)
     324
     325        endif
     326      enddo
     327
     328
     329      !============================================================================
     330      ! 5. calcul de la vitesse verticale en melangeant Tl et qt du thermique
     331      !===========================================================================
     332
     333      activetmp(:) = active(:) .and. f_star(:, l + 1)>1.e-10
     334      do ig = 1, ngrid
     335        if (activetmp(ig)) then
     336          Zsat = .false.
     337          ztla(ig, l) = (f_star(ig, l) * ztla(ig, l - 1) + &
     338                  (alim_star(ig, l) + entr_star(ig, l)) * zthl(ig, l))  &
     339                  / (f_star(ig, l + 1) + detr_star(ig, l))
     340          zqta(ig, l) = (f_star(ig, l) * zqta(ig, l - 1) + &
     341                  (alim_star(ig, l) + entr_star(ig, l)) * po(ig, l))  &
     342                  / (f_star(ig, l + 1) + detr_star(ig, l))
     343
     344        endif
     345      enddo
     346
     347      ztemp(:) = zpspsk(:, l) * ztla(:, l)
     348      call thermcell_qsat(ngrid, activetmp, pplev(:, l), ztemp, zqta(:, l), zqsatth(:, l))
     349      do ig = 1, ngrid
     350        if (activetmp(ig)) then
     351          ! on ecrit de maniere conservative (sat ou non)
     352          !          T = Tl +Lv/Cp ql
     353          zqla(ig, l) = max(0., zqta(ig, l) - zqsatth(ig, l))
     354          ztva(ig, l) = ztla(ig, l) * zpspsk(ig, l) + RLvCp * zqla(ig, l)
     355          ztva(ig, l) = ztva(ig, l) / zpspsk(ig, l)
     356          !on rajoute le calcul de zha pour diagnostiques (temp potentielle)
     357          zha(ig, l) = ztva(ig, l)
     358          ztva(ig, l) = ztva(ig, l) * (1. + RETV * (zqta(ig, l)  &
     359                  - zqla(ig, l)) - zqla(ig, l))
     360          zbuoy(ig, l) = RG * (ztva(ig, l) - ztv(ig, l)) / ztv(ig, l)
     361          zdz = zlev(ig, l + 1) - zlev(ig, l)
     362          zdzbis = zlev(ig, l) - zlev(ig, l - 1)
     363          zeps(ig, l) = (entr_star(ig, l) + alim_star(ig, l)) / (f_star(ig, l) * zdz)
     364          zw2fact = fact_epsilon * 2. * zdz / (1. + betalpha)
     365          zw2factbis = fact_epsilon * 2. * zdzbis / (1. + betalpha)
     366          zdw2 = afact * zbuoy(ig, l) / (fact_epsilon)
     367          zdw2bis = afact * zbuoy(ig, l - 1) / (fact_epsilon)
     368          zw2(ig, l + 1) = Max(0.0001, exp(-zw2fact) * (zw2(ig, l) - zdw2) + zdw2)
     369        endif
     370      enddo
     371
     372      if (prt_level>=20) print*, 'coucou calcul detr 460: ig, l', ig, l
     373
     374      !===========================================================================
     375      ! 6. initialisations pour le calcul de la hauteur du thermique, de l'inversion et de la vitesse verticale max
     376      !===========================================================================
     377
     378      nbpb = 0
     379      do ig = 1, ngrid
     380        if (zw2(ig, l + 1)>0. .and. zw2(ig, l + 1)<1.e-10) then
     381          !               stop'On tombe sur le cas particulier de thermcell_dry'
     382          !               print*,'On tombe sur le cas particulier de thermcell_plume'
     383          nbpb = nbpb + 1
     384          zw2(ig, l + 1) = 0.
     385          linter(ig) = l + 1
     386        endif
     387
     388        if (zw2(ig, l + 1)<0.) then
     389          linter(ig) = (l * (zw2(ig, l + 1) - zw2(ig, l))  &
     390                  - zw2(ig, l)) / (zw2(ig, l + 1) - zw2(ig, l))
     391          zw2(ig, l + 1) = 0.
     392          !+CR:04/05/12:correction calcul linter pour calcul de zmax continu
     393        elseif (f_star(ig, l + 1)<0.) then
     394          linter(ig) = (l * (f_star(ig, l + 1) - f_star(ig, l))  &
     395                  - f_star(ig, l)) / (f_star(ig, l + 1) - f_star(ig, l))
     396          zw2(ig, l + 1) = 0.
     397          !fin CR:04/05/12
     398        endif
     399
     400        wa_moy(ig, l + 1) = sqrt(zw2(ig, l + 1))
     401
     402        if (wa_moy(ig, l + 1)>wmaxa(ig)) then
     403          !   lmix est le niveau de la couche ou w (wa_moy) est maximum
     404          !on rajoute le calcul de lmix_bis
     405          if (zqla(ig, l)<1.e-10) then
     406            lmix_bis(ig) = l + 1
     407          endif
     408          lmix(ig) = l + 1
     409          wmaxa(ig) = wa_moy(ig, l + 1)
     410        endif
     411      enddo
     412
     413      if (nbpb>0) then
     414        print*, 'WARNING on tombe ', nbpb, ' x sur un pb pour l=', l, ' dans thermcell_plume'
     415      endif
     416
     417      !=========================================================================
     418      ! FIN DE LA BOUCLE VERTICALE
    282419    enddo
    283 
    284 
    285 !-------------------------------------------------
    286 !calcul des taux d'entrainement et de detrainement
    287 !-------------------------------------------------
    288 
    289      do ig=1,ngrid
    290         if (active(ig)) then
    291 
    292 !          zw2m=max(0.5*(w_est(ig,l)+w_est(ig,l+1)),0.1)
    293           zw2m=w_est(ig,l+1)
    294           zdz=zlev(ig,l+1)-zlev(ig,l)
    295           zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
    296           zalpha=f0(ig)*f_star(ig,l)/sqrt(w_est(ig,l+1))/rhobarz(ig,l)
    297           zdqt(ig,l)=max(zqta(ig,l-1)-po(ig,l),0.)/po(ig,l)
    298 
    299 !=========================================================================
    300 ! 4. Calcul de l'entrainement et du detrainement
    301 !=========================================================================
    302 
    303           detr_star(ig,l)=f_star(ig,l)*zdz             &
    304        *( mix0 * 0.1 / (zalpha+0.001)               &
    305        + MAX(detr_min, -afact*zbetalpha*zbuoyjam(ig,l)/zw2m   &
    306        + detr_q_coef*(zdqt(ig,l)/zw2m)**detr_q_power))
    307 
    308           if ( iflag_thermals_ed == 20 ) then
    309              entr_star(ig,l)=f_star(ig,l)*zdz* (         &
    310             mix0 * 0.1 / (zalpha+0.001)               &
    311           + zbetalpha*MAX(entr_min,                   &
    312           afact*zbuoyjam(ig,l)/zw2m - fact_epsilon))
    313           else
    314              entr_star(ig,l)=f_star(ig,l)*zdz* (         &
    315             mix0 * 0.1 / (zalpha+0.001)               &
    316           + zbetalpha*MAX(entr_min,                   &
    317           afact*zbuoy(ig,l)/zw2m - fact_epsilon))
    318           endif
    319          
    320 ! En dessous de lalim, on prend le max de alim_star et entr_star pour
    321 ! alim_star et 0 sinon
    322         if (l<lalim(ig)) then
    323           alim_star(ig,l)=max(alim_star(ig,l),entr_star(ig,l))
    324           entr_star(ig,l)=0.
    325         endif
    326         f_star(ig,l+1)=f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l)  &
    327                 -detr_star(ig,l)
    328 
    329       endif
    330    enddo
    331 
    332 
    333 !============================================================================
    334 ! 5. calcul de la vitesse verticale en melangeant Tl et qt du thermique
    335 !===========================================================================
    336 
    337    activetmp(:)=active(:) .and. f_star(:,l+1)>1.e-10
    338    do ig=1,ngrid
    339        if (activetmp(ig)) then
    340            Zsat=.false.
    341            ztla(ig,l)=(f_star(ig,l)*ztla(ig,l-1)+  &
    342               (alim_star(ig,l)+entr_star(ig,l))*zthl(ig,l))  &
    343               /(f_star(ig,l+1)+detr_star(ig,l))
    344            zqta(ig,l)=(f_star(ig,l)*zqta(ig,l-1)+  &
    345               (alim_star(ig,l)+entr_star(ig,l))*po(ig,l))  &
    346               /(f_star(ig,l+1)+detr_star(ig,l))
    347 
    348         endif
     420    !=========================================================================
     421
     422    !on recalcule alim_star_tot
     423    do ig = 1, ngrid
     424      alim_star_tot(ig) = 0.
    349425    enddo
    350 
    351    ztemp(:)=zpspsk(:,l)*ztla(:,l)
    352    call thermcell_qsat(ngrid,activetmp,pplev(:,l),ztemp,zqta(:,l),zqsatth(:,l))
    353    do ig=1,ngrid
    354       if (activetmp(ig)) then
    355 ! on ecrit de maniere conservative (sat ou non)
    356 !          T = Tl +Lv/Cp ql
    357            zqla(ig,l)=max(0.,zqta(ig,l)-zqsatth(ig,l))
    358            ztva(ig,l) = ztla(ig,l)*zpspsk(ig,l)+RLvCp*zqla(ig,l)
    359            ztva(ig,l) = ztva(ig,l)/zpspsk(ig,l)
    360 !on rajoute le calcul de zha pour diagnostiques (temp potentielle)
    361            zha(ig,l) = ztva(ig,l)
    362            ztva(ig,l) = ztva(ig,l)*(1.+RETV*(zqta(ig,l)  &
    363                 -zqla(ig,l))-zqla(ig,l))
    364            zbuoy(ig,l)=RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)
    365            zdz=zlev(ig,l+1)-zlev(ig,l)
    366            zdzbis=zlev(ig,l)-zlev(ig,l-1)
    367            zeps(ig,l)=(entr_star(ig,l)+alim_star(ig,l))/(f_star(ig,l)*zdz)
    368            zw2fact=fact_epsilon*2.*zdz/(1.+betalpha)
    369            zw2factbis=fact_epsilon*2.*zdzbis/(1.+betalpha)
    370            zdw2= afact*zbuoy(ig,l)/(fact_epsilon)
    371            zdw2bis= afact*zbuoy(ig,l-1)/(fact_epsilon)
    372            zw2(ig,l+1)=Max(0.0001,exp(-zw2fact)*(zw2(ig,l)-zdw2)+zdw2)
    373       endif
    374    enddo
    375 
    376    if (prt_level>=20) print*,'coucou calcul detr 460: ig, l',ig, l
    377 
    378 !===========================================================================
    379 ! 6. initialisations pour le calcul de la hauteur du thermique, de l'inversion et de la vitesse verticale max
    380 !===========================================================================
    381 
    382    nbpb=0
    383    do ig=1,ngrid
    384             if (zw2(ig,l+1)>0. .and. zw2(ig,l+1)<1.e-10) then
    385 !               stop'On tombe sur le cas particulier de thermcell_dry'
    386 !               print*,'On tombe sur le cas particulier de thermcell_plume'
    387                 nbpb=nbpb+1
    388                 zw2(ig,l+1)=0.
    389                 linter(ig)=l+1
    390             endif
    391 
    392         if (zw2(ig,l+1)<0.) then
    393            linter(ig)=(l*(zw2(ig,l+1)-zw2(ig,l))  &
    394                  -zw2(ig,l))/(zw2(ig,l+1)-zw2(ig,l))
    395            zw2(ig,l+1)=0.
    396 !+CR:04/05/12:correction calcul linter pour calcul de zmax continu
    397         elseif (f_star(ig,l+1)<0.) then
    398            linter(ig)=(l*(f_star(ig,l+1)-f_star(ig,l))  &
    399                  -f_star(ig,l))/(f_star(ig,l+1)-f_star(ig,l))
    400            zw2(ig,l+1)=0.
    401 !fin CR:04/05/12
    402         endif
    403 
    404            wa_moy(ig,l+1)=sqrt(zw2(ig,l+1))
    405 
    406         if (wa_moy(ig,l+1)>wmaxa(ig)) then
    407 !   lmix est le niveau de la couche ou w (wa_moy) est maximum
    408 !on rajoute le calcul de lmix_bis
    409             if (zqla(ig,l)<1.e-10) then
    410                lmix_bis(ig)=l+1
    411             endif
    412             lmix(ig)=l+1
    413             wmaxa(ig)=wa_moy(ig,l+1)
    414         endif
    415    enddo
    416 
    417    if (nbpb>0) then
    418    print*,'WARNING on tombe ',nbpb,' x sur un pb pour l=',l,' dans thermcell_plume'
    419    endif
    420 
    421 !=========================================================================
    422 ! FIN DE LA BOUCLE VERTICALE
    423       enddo
    424 !=========================================================================
    425 
    426 !on recalcule alim_star_tot
    427        do ig=1,ngrid
    428           alim_star_tot(ig)=0.
    429        enddo
    430        do ig=1,ngrid
    431           do l=1,lalim(ig)-1
    432           alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
    433           enddo
    434        enddo
    435        
    436 
    437         if (prt_level>=20) print*,'coucou calcul detr 470: ig, l', ig, l
    438 
    439 #undef wrgrads_thermcell
    440 #ifdef wrgrads_thermcell
    441          call wrgradsfi(1,nlay,entr_star(igout,1:nlay),'esta      ','esta      ')
    442          call wrgradsfi(1,nlay,detr_star(igout,1:nlay),'dsta      ','dsta      ')
    443          call wrgradsfi(1,nlay,zbuoy(igout,1:nlay),'buoy      ','buoy      ')
    444          call wrgradsfi(1,nlay,zdqt(igout,1:nlay),'dqt      ','dqt      ')
    445          call wrgradsfi(1,nlay,w_est(igout,1:nlay),'w_est     ','w_est     ')
    446          call wrgradsfi(1,nlay,w_est(igout,2:nlay+1),'w_es2     ','w_es2     ')
    447          call wrgradsfi(1,nlay,zw2(igout,1:nlay),'zw2A      ','zw2A      ')
    448 #endif
    449 
    450 
    451  RETURN
    452      end
     426    do ig = 1, ngrid
     427      do l = 1, lalim(ig) - 1
     428        alim_star_tot(ig) = alim_star_tot(ig) + alim_star(ig, l)
     429      enddo
     430    enddo
     431
     432    if (prt_level>=20) print*, 'coucou calcul detr 470: ig, l', ig, l
     433    RETURN
     434  end
    453435END MODULE lmdz_thermcell_plume

  • LMDZ6/branches/Amaury_dev/libf/phylmd/lmdz_thermcell_plume_6A.f90

    r5101 r5102  
     1! $Id$
     2
     3
    14MODULE lmdz_thermcell_plume_6A
    25
    3 ! $Id$
    4 
    56CONTAINS
    67
    7       SUBROUTINE thermcell_plume_6A(itap,ngrid,nlay,ptimestep,ztv,zthl,po,zl,rhobarz,  &
    8              zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot,  &
    9              lalim,f0,detr_star,entr_star,f_star,csc,ztva,  &
    10              ztla,zqla,zqta,zha,zw2,w_est,ztva_est,zqsatth,lmix,lmix_bis,linter &
    11              ,lev_out,lunout1,igout)
    12 !     &           ,lev_out,lunout1,igout,zbuoy,zbuoyjam)
    13 !--------------------------------------------------------------------------
    14 !thermcell_plume: calcule les valeurs de qt, thetal et w dans l ascendance
    15 !--------------------------------------------------------------------------
    16 
    17        USE lmdz_thermcell_ini, ONLY: prt_level,fact_thermals_ed_dz,iflag_thermals_ed,RLvCP,RETV,RG
    18        USE lmdz_thermcell_ini, ONLY: fact_epsilon, betalpha, afact, fact_shell
    19        USE lmdz_thermcell_ini, ONLY: detr_min, entr_min, detr_q_coef, detr_q_power
    20        USE lmdz_thermcell_ini, ONLY: mix0, thermals_flag_alim
    21        USE lmdz_thermcell_alim, ONLY: thermcell_alim
    22        USE lmdz_thermcell_qsat, ONLY: thermcell_qsat
    23 
    24 
    25        IMPLICIT NONE
    26 
    27       integer,intent(in) :: itap,lev_out,lunout1,igout,ngrid,nlay
    28       real,intent(in) :: ptimestep
    29       real,intent(in),dimension(ngrid,nlay) :: ztv
    30       real,intent(in),dimension(ngrid,nlay) :: zthl
    31       real,intent(in),dimension(ngrid,nlay) :: po
    32       real,intent(in),dimension(ngrid,nlay) :: zl
    33       real,intent(in),dimension(ngrid,nlay) :: rhobarz
    34       real,intent(in),dimension(ngrid,nlay+1) :: zlev
    35       real,intent(in),dimension(ngrid,nlay+1) :: pplev
    36       real,intent(in),dimension(ngrid,nlay) :: pphi
    37       real,intent(in),dimension(ngrid,nlay) :: zpspsk
    38       real,intent(in),dimension(ngrid) :: f0
    39 
    40       integer,intent(out) :: lalim(ngrid)
    41       real,intent(out),dimension(ngrid,nlay) :: alim_star
    42       real,intent(out),dimension(ngrid) :: alim_star_tot
    43       real,intent(out),dimension(ngrid,nlay) :: detr_star
    44       real,intent(out),dimension(ngrid,nlay) :: entr_star
    45       real,intent(out),dimension(ngrid,nlay+1) :: f_star
    46       real,intent(out),dimension(ngrid,nlay) :: csc
    47       real,intent(out),dimension(ngrid,nlay) :: ztva
    48       real,intent(out),dimension(ngrid,nlay) :: ztla
    49       real,intent(out),dimension(ngrid,nlay) :: zqla
    50       real,intent(out),dimension(ngrid,nlay) :: zqta
    51       real,intent(out),dimension(ngrid,nlay) :: zha
    52       real,intent(out),dimension(ngrid,nlay+1) :: zw2
    53       real,intent(out),dimension(ngrid,nlay+1) :: w_est
    54       real,intent(out),dimension(ngrid,nlay) :: ztva_est
    55       real,intent(out),dimension(ngrid,nlay) :: zqsatth
    56       integer,intent(out),dimension(ngrid) :: lmix
    57       integer,intent(out),dimension(ngrid) :: lmix_bis
    58       real,intent(out),dimension(ngrid) :: linter
    59 
    60       REAL zdw2,zdw2bis
    61       REAL zw2modif
    62       REAL zw2fact,zw2factbis
    63       REAL,dimension(ngrid,nlay) :: zeps
    64 
    65       REAL, dimension(ngrid) ::    wmaxa
    66 
    67       INTEGER ig,l,k,lt,it,lm
    68       integer nbpb
    69 
    70       real,dimension(ngrid,nlay) :: detr
    71       real,dimension(ngrid,nlay) :: entr
    72       real,dimension(ngrid,nlay+1) :: wa_moy
    73       real,dimension(ngrid,nlay) :: ztv_est
    74       real,dimension(ngrid) :: ztemp,zqsat
    75       real,dimension(ngrid,nlay) :: zqla_est
    76       real,dimension(ngrid,nlay) :: zta_est
    77 
    78       real,dimension(ngrid,nlay) :: zbuoy,gamma,zdqt
    79       real zdz,zalpha,zw2m
    80       real,dimension(ngrid,nlay) :: zbuoyjam,zdqtjam
    81       real zbuoybis,zdz2,zdz3,lmel,entrbis,zdzbis
    82       real, dimension(ngrid) :: d_temp
    83       real ztv1,ztv2,factinv,zinv,zlmel
    84       real zlmelup,zlmeldwn,zlt,zltdwn,zltup
    85       real atv1,atv2,btv1,btv2
    86       real ztv_est1,ztv_est2
    87       real zcor,zdelta,zcvm5,qlbef
    88       real zbetalpha, coefzlmel
    89       real eps
    90       logical Zsat
    91       LOGICAL,dimension(ngrid) :: active,activetmp
    92       REAL fact_gamma,fact_gamma2,fact_epsilon2
    93       REAL coefc
    94       REAL,dimension(ngrid,nlay) :: c2
    95 
    96       if (ngrid==1) print*,'THERMCELL PLUME MODIFIE 2014/07/11'
    97       Zsat=.false.
    98 ! Initialisation
    99 
    100 
    101       zbetalpha=betalpha/(1.+betalpha)
    102 
    103 
    104 ! Initialisations des variables r?elles
    105 if (1==1) then
    106       ztva(:,:)=ztv(:,:)
    107       ztva_est(:,:)=ztva(:,:)
    108       ztv_est(:,:)=ztv(:,:)
    109       ztla(:,:)=zthl(:,:)
    110       zqta(:,:)=po(:,:)
    111       zqla(:,:)=0.
    112       zha(:,:) = ztva(:,:)
    113 else
    114       ztva(:,:)=0.
    115       ztv_est(:,:)=0.
    116       ztva_est(:,:)=0.
    117       ztla(:,:)=0.
    118       zqta(:,:)=0.
    119       zha(:,:) =0.
    120 endif
    121 
    122       zqla_est(:,:)=0.
    123       zqsatth(:,:)=0.
    124       zqla(:,:)=0.
    125       detr_star(:,:)=0.
    126       entr_star(:,:)=0.
    127       alim_star(:,:)=0.
    128       alim_star_tot(:)=0.
    129       csc(:,:)=0.
    130       detr(:,:)=0.
    131       entr(:,:)=0.
    132       zw2(:,:)=0.
    133       zbuoy(:,:)=0.
    134       zbuoyjam(:,:)=0.
    135       gamma(:,:)=0.
    136       zeps(:,:)=0.
    137       w_est(:,:)=0.
    138       f_star(:,:)=0.
    139       wa_moy(:,:)=0.
    140       linter(:)=1.
    141 !     linter(:)=1.
    142 ! Initialisation des variables entieres
    143       lmix(:)=1
    144       lmix_bis(:)=2
    145       wmaxa(:)=0.
    146 
    147 ! Initialisation a 0  en cas de sortie dans replay
    148       zqsat(:)=0.
    149       zta_est(:,:)=0.
    150       zdqt(:,:)=0.
    151       zdqtjam(:,:)=0.
    152       c2(:,:)=0.
    153 
    154 
    155 !-------------------------------------------------------------------------
    156 ! On ne considere comme actif que les colonnes dont les deux premieres
    157 ! couches sont instables.
    158 !-------------------------------------------------------------------------
    159 
    160       active(:)=ztv(:,1)>ztv(:,2)
    161       d_temp(:)=0. ! Pour activer un contraste de temperature a la base
    162                    ! du panache
    163 !  Cet appel pourrait être fait avant thermcell_plume dans thermcell_main
    164       CALL thermcell_alim(thermals_flag_alim,ngrid,nlay,ztv,d_temp,zlev,alim_star,lalim)
    165 
    166 !------------------------------------------------------------------------------
    167 ! Calcul dans la premiere couche
    168 ! On decide dans cette version que le thermique n'est actif que si la premiere
    169 ! couche est instable.
    170 ! Pourrait etre change si on veut que le thermiques puisse se d??clencher
    171 ! dans une couche l>1
    172 !------------------------------------------------------------------------------
    173 do ig=1,ngrid
    174 ! Le panache va prendre au debut les caracteristiques de l'air contenu
    175 ! dans cette couche.
    176     if (active(ig)) then
    177     ztla(ig,1)=zthl(ig,1)
    178     zqta(ig,1)=po(ig,1)
    179     zqla(ig,1)=zl(ig,1)
    180 !cr: attention, prise en compte de f*(1)=1
    181     f_star(ig,2)=alim_star(ig,1)
    182     zw2(ig,2)=2.*RG*(ztv(ig,1)-ztv(ig,2))/ztv(ig,2)  &
    183 &                     *(zlev(ig,2)-zlev(ig,1))  &
    184 &                     *0.4*pphi(ig,1)/(pphi(ig,2)-pphi(ig,1))
    185     w_est(ig,2)=zw2(ig,2)
     8  SUBROUTINE thermcell_plume_6A(itap, ngrid, nlay, ptimestep, ztv, zthl, po, zl, rhobarz, &
     9          zlev, pplev, pphi, zpspsk, alim_star, alim_star_tot, &
     10          lalim, f0, detr_star, entr_star, f_star, csc, ztva, &
     11          ztla, zqla, zqta, zha, zw2, w_est, ztva_est, zqsatth, lmix, lmix_bis, linter &
     12          , lev_out, lunout1, igout)
     13    !     &           ,lev_out,lunout1,igout,zbuoy,zbuoyjam)
     14    !--------------------------------------------------------------------------
     15    !thermcell_plume: calcule les valeurs de qt, thetal et w dans l ascendance
     16    !--------------------------------------------------------------------------
     17
     18    USE lmdz_thermcell_ini, ONLY: prt_level, fact_thermals_ed_dz, iflag_thermals_ed, RLvCP, RETV, RG
     19    USE lmdz_thermcell_ini, ONLY: fact_epsilon, betalpha, afact, fact_shell
     20    USE lmdz_thermcell_ini, ONLY: detr_min, entr_min, detr_q_coef, detr_q_power
     21    USE lmdz_thermcell_ini, ONLY: mix0, thermals_flag_alim
     22    USE lmdz_thermcell_alim, ONLY: thermcell_alim
     23    USE lmdz_thermcell_qsat, ONLY: thermcell_qsat
     24
     25    IMPLICIT NONE
     26
     27    integer, intent(in) :: itap, lev_out, lunout1, igout, ngrid, nlay
     28    real, intent(in) :: ptimestep
     29    real, intent(in), dimension(ngrid, nlay) :: ztv
     30    real, intent(in), dimension(ngrid, nlay) :: zthl
     31    real, intent(in), dimension(ngrid, nlay) :: po
     32    real, intent(in), dimension(ngrid, nlay) :: zl
     33    real, intent(in), dimension(ngrid, nlay) :: rhobarz
     34    real, intent(in), dimension(ngrid, nlay + 1) :: zlev
     35    real, intent(in), dimension(ngrid, nlay + 1) :: pplev
     36    real, intent(in), dimension(ngrid, nlay) :: pphi
     37    real, intent(in), dimension(ngrid, nlay) :: zpspsk
     38    real, intent(in), dimension(ngrid) :: f0
     39
     40    integer, intent(out) :: lalim(ngrid)
     41    real, intent(out), dimension(ngrid, nlay) :: alim_star
     42    real, intent(out), dimension(ngrid) :: alim_star_tot
     43    real, intent(out), dimension(ngrid, nlay) :: detr_star
     44    real, intent(out), dimension(ngrid, nlay) :: entr_star
     45    real, intent(out), dimension(ngrid, nlay + 1) :: f_star
     46    real, intent(out), dimension(ngrid, nlay) :: csc
     47    real, intent(out), dimension(ngrid, nlay) :: ztva
     48    real, intent(out), dimension(ngrid, nlay) :: ztla
     49    real, intent(out), dimension(ngrid, nlay) :: zqla
     50    real, intent(out), dimension(ngrid, nlay) :: zqta
     51    real, intent(out), dimension(ngrid, nlay) :: zha
     52    real, intent(out), dimension(ngrid, nlay + 1) :: zw2
     53    real, intent(out), dimension(ngrid, nlay + 1) :: w_est
     54    real, intent(out), dimension(ngrid, nlay) :: ztva_est
     55    real, intent(out), dimension(ngrid, nlay) :: zqsatth
     56    integer, intent(out), dimension(ngrid) :: lmix
     57    integer, intent(out), dimension(ngrid) :: lmix_bis
     58    real, intent(out), dimension(ngrid) :: linter
     59
     60    REAL zdw2, zdw2bis
     61    REAL zw2modif
     62    REAL zw2fact, zw2factbis
     63    REAL, dimension(ngrid, nlay) :: zeps
     64
     65    REAL, dimension(ngrid) :: wmaxa
     66
     67    INTEGER ig, l, k, lt, it, lm
     68    integer nbpb
     69
     70    real, dimension(ngrid, nlay) :: detr
     71    real, dimension(ngrid, nlay) :: entr
     72    real, dimension(ngrid, nlay + 1) :: wa_moy
     73    real, dimension(ngrid, nlay) :: ztv_est
     74    real, dimension(ngrid) :: ztemp, zqsat
     75    real, dimension(ngrid, nlay) :: zqla_est
     76    real, dimension(ngrid, nlay) :: zta_est
     77
     78    real, dimension(ngrid, nlay) :: zbuoy, gamma, zdqt
     79    real zdz, zalpha, zw2m
     80    real, dimension(ngrid, nlay) :: zbuoyjam, zdqtjam
     81    real zbuoybis, zdz2, zdz3, lmel, entrbis, zdzbis
     82    real, dimension(ngrid) :: d_temp
     83    real ztv1, ztv2, factinv, zinv, zlmel
     84    real zlmelup, zlmeldwn, zlt, zltdwn, zltup
     85    real atv1, atv2, btv1, btv2
     86    real ztv_est1, ztv_est2
     87    real zcor, zdelta, zcvm5, qlbef
     88    real zbetalpha, coefzlmel
     89    real eps
     90    logical Zsat
     91    LOGICAL, dimension(ngrid) :: active, activetmp
     92    REAL fact_gamma, fact_gamma2, fact_epsilon2
     93    REAL coefc
     94    REAL, dimension(ngrid, nlay) :: c2
     95
     96    if (ngrid==1) print*, 'THERMCELL PLUME MODIFIE 2014/07/11'
     97    Zsat = .false.
     98    ! Initialisation
     99
     100    zbetalpha = betalpha / (1. + betalpha)
     101
     102
     103    ! Initialisations des variables r?elles
     104    if (1==1) then
     105      ztva(:, :) = ztv(:, :)
     106      ztva_est(:, :) = ztva(:, :)
     107      ztv_est(:, :) = ztv(:, :)
     108      ztla(:, :) = zthl(:, :)
     109      zqta(:, :) = po(:, :)
     110      zqla(:, :) = 0.
     111      zha(:, :) = ztva(:, :)
     112    else
     113      ztva(:, :) = 0.
     114      ztv_est(:, :) = 0.
     115      ztva_est(:, :) = 0.
     116      ztla(:, :) = 0.
     117      zqta(:, :) = 0.
     118      zha(:, :) = 0.
    186119    endif
    187 enddo
    188 
    189 !==============================================================================
    190 !boucle de calcul de la vitesse verticale dans le thermique
    191 !==============================================================================
    192 do l=2,nlay-1
    193 !==============================================================================
    194 
    195 
    196 ! On decide si le thermique est encore actif ou non
    197 ! AFaire : Il faut sans doute ajouter entr_star a alim_star dans ce test
    198     do ig=1,ngrid
    199        active(ig)=active(ig) &
    200 &                 .and. zw2(ig,l)>1.e-10 &
    201 &                 .and. f_star(ig,l)+alim_star(ig,l)>1.e-10
     120
     121    zqla_est(:, :) = 0.
     122    zqsatth(:, :) = 0.
     123    zqla(:, :) = 0.
     124    detr_star(:, :) = 0.
     125    entr_star(:, :) = 0.
     126    alim_star(:, :) = 0.
     127    alim_star_tot(:) = 0.
     128    csc(:, :) = 0.
     129    detr(:, :) = 0.
     130    entr(:, :) = 0.
     131    zw2(:, :) = 0.
     132    zbuoy(:, :) = 0.
     133    zbuoyjam(:, :) = 0.
     134    gamma(:, :) = 0.
     135    zeps(:, :) = 0.
     136    w_est(:, :) = 0.
     137    f_star(:, :) = 0.
     138    wa_moy(:, :) = 0.
     139    linter(:) = 1.
     140    !     linter(:)=1.
     141    ! Initialisation des variables entieres
     142    lmix(:) = 1
     143    lmix_bis(:) = 2
     144    wmaxa(:) = 0.
     145
     146    ! Initialisation a 0  en cas de sortie dans replay
     147    zqsat(:) = 0.
     148    zta_est(:, :) = 0.
     149    zdqt(:, :) = 0.
     150    zdqtjam(:, :) = 0.
     151    c2(:, :) = 0.
     152
     153
     154    !-------------------------------------------------------------------------
     155    ! On ne considere comme actif que les colonnes dont les deux premieres
     156    ! couches sont instables.
     157    !-------------------------------------------------------------------------
     158
     159    active(:) = ztv(:, 1)>ztv(:, 2)
     160    d_temp(:) = 0. ! Pour activer un contraste de temperature a la base
     161    ! du panache
     162    !  Cet appel pourrait être fait avant thermcell_plume dans thermcell_main
     163    CALL thermcell_alim(thermals_flag_alim, ngrid, nlay, ztv, d_temp, zlev, alim_star, lalim)
     164
     165    !------------------------------------------------------------------------------
     166    ! Calcul dans la premiere couche
     167    ! On decide dans cette version que le thermique n'est actif que si la premiere
     168    ! couche est instable.
     169    ! Pourrait etre change si on veut que le thermiques puisse se d??clencher
     170    ! dans une couche l>1
     171    !------------------------------------------------------------------------------
     172    do ig = 1, ngrid
     173      ! Le panache va prendre au debut les caracteristiques de l'air contenu
     174      ! dans cette couche.
     175      if (active(ig)) then
     176        ztla(ig, 1) = zthl(ig, 1)
     177        zqta(ig, 1) = po(ig, 1)
     178        zqla(ig, 1) = zl(ig, 1)
     179        !cr: attention, prise en compte de f*(1)=1
     180        f_star(ig, 2) = alim_star(ig, 1)
     181        zw2(ig, 2) = 2. * RG * (ztv(ig, 1) - ztv(ig, 2)) / ztv(ig, 2)  &
     182                & * (zlev(ig, 2) - zlev(ig, 1))  &
     183                & * 0.4 * pphi(ig, 1) / (pphi(ig, 2) - pphi(ig, 1))
     184        w_est(ig, 2) = zw2(ig, 2)
     185      endif
    202186    enddo
    203187
    204 
    205 
    206 !---------------------------------------------------------------------------
    207 ! calcul des proprietes thermodynamiques et de la vitesse de la couche l
    208 ! sans tenir compte du detrainement et de l'entrainement dans cette
    209 ! couche
    210 ! C'est a dire qu'on suppose
    211 ! ztla(l)=ztla(l-1) et zqta(l)=zqta(l-1)
    212 ! Ici encore, on doit pouvoir ajouter entr_star (qui peut etre calculer
    213 ! avant) a l'alimentation pour avoir un calcul plus propre
    214 !---------------------------------------------------------------------------
    215 
    216    ztemp(:)=zpspsk(:,l)*ztla(:,l-1)
    217    call thermcell_qsat(ngrid,active,pplev(:,l),ztemp,zqta(:,l-1),zqsat(:))
    218     do ig=1,ngrid
    219 !       print*,'active',active(ig),ig,l
    220         if(active(ig)) then
    221         zqla_est(ig,l)=max(0.,zqta(ig,l-1)-zqsat(ig))
    222         ztva_est(ig,l) = ztla(ig,l-1)*zpspsk(ig,l)+RLvCp*zqla_est(ig,l)
    223         zta_est(ig,l)=ztva_est(ig,l)
    224         ztva_est(ig,l) = ztva_est(ig,l)/zpspsk(ig,l)
    225         ztva_est(ig,l) = ztva_est(ig,l)*(1.+RETV*(zqta(ig,l-1)  &
    226         -zqla_est(ig,l))-zqla_est(ig,l))
    227  
    228 
    229 !Modif AJAM
    230 
    231         zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
    232         zdz=zlev(ig,l+1)-zlev(ig,l)         
    233         lmel=fact_thermals_ed_dz*zlev(ig,l)
    234 !        lmel=0.09*zlev(ig,l)
    235         zlmel=zlev(ig,l)+lmel
    236         zlmelup=zlmel+(zdz/2)
    237         zlmeldwn=zlmel-(zdz/2)
    238 
    239         lt=l+1
    240         zlt=zlev(ig,lt)
    241         zdz3=zlev(ig,lt+1)-zlt
    242         zltdwn=zlt-zdz3/2
    243         zltup=zlt+zdz3/2
    244          
    245 !=========================================================================
    246 ! 3. Calcul de la flotabilite modifie par melange avec l'air au dessus
    247 !=========================================================================
    248 
    249 !--------------------------------------------------
    250         if (iflag_thermals_ed<8) then
    251 !--------------------------------------------------
    252 !AJ052014: J'ai remplac?? la boucle do par un do while
    253 ! afin de faire moins de calcul dans la boucle
    254 !--------------------------------------------------
     188    !==============================================================================
     189    !boucle de calcul de la vitesse verticale dans le thermique
     190    !==============================================================================
     191    do l = 2, nlay - 1
     192      !==============================================================================
     193
     194
     195      ! On decide si le thermique est encore actif ou non
     196      ! AFaire : Il faut sans doute ajouter entr_star a alim_star dans ce test
     197      do ig = 1, ngrid
     198        active(ig) = active(ig) &
     199                &                 .and. zw2(ig, l)>1.e-10 &
     200                &                 .and. f_star(ig, l) + alim_star(ig, l)>1.e-10
     201      enddo
     202
     203
     204
     205      !---------------------------------------------------------------------------
     206      ! calcul des proprietes thermodynamiques et de la vitesse de la couche l
     207      ! sans tenir compte du detrainement et de l'entrainement dans cette
     208      ! couche
     209      ! C'est a dire qu'on suppose
     210      ! ztla(l)=ztla(l-1) et zqta(l)=zqta(l-1)
     211      ! Ici encore, on doit pouvoir ajouter entr_star (qui peut etre calculer
     212      ! avant) a l'alimentation pour avoir un calcul plus propre
     213      !---------------------------------------------------------------------------
     214
     215      ztemp(:) = zpspsk(:, l) * ztla(:, l - 1)
     216      call thermcell_qsat(ngrid, active, pplev(:, l), ztemp, zqta(:, l - 1), zqsat(:))
     217      do ig = 1, ngrid
     218        !       print*,'active',active(ig),ig,l
     219        if(active(ig)) then
     220          zqla_est(ig, l) = max(0., zqta(ig, l - 1) - zqsat(ig))
     221          ztva_est(ig, l) = ztla(ig, l - 1) * zpspsk(ig, l) + RLvCp * zqla_est(ig, l)
     222          zta_est(ig, l) = ztva_est(ig, l)
     223          ztva_est(ig, l) = ztva_est(ig, l) / zpspsk(ig, l)
     224          ztva_est(ig, l) = ztva_est(ig, l) * (1. + RETV * (zqta(ig, l - 1)  &
     225                  - zqla_est(ig, l)) - zqla_est(ig, l))
     226
     227
     228          !Modif AJAM
     229
     230          zbuoy(ig, l) = RG * (ztva_est(ig, l) - ztv(ig, l)) / ztv(ig, l)
     231          zdz = zlev(ig, l + 1) - zlev(ig, l)
     232          lmel = fact_thermals_ed_dz * zlev(ig, l)
     233          !        lmel=0.09*zlev(ig,l)
     234          zlmel = zlev(ig, l) + lmel
     235          zlmelup = zlmel + (zdz / 2)
     236          zlmeldwn = zlmel - (zdz / 2)
     237
     238          lt = l + 1
     239          zlt = zlev(ig, lt)
     240          zdz3 = zlev(ig, lt + 1) - zlt
     241          zltdwn = zlt - zdz3 / 2
     242          zltup = zlt + zdz3 / 2
     243
     244          !=========================================================================
     245          ! 3. Calcul de la flotabilite modifie par melange avec l'air au dessus
     246          !=========================================================================
     247
     248          !--------------------------------------------------
     249          if (iflag_thermals_ed<8) then
     250            !--------------------------------------------------
     251            !AJ052014: J'ai remplac?? la boucle do par un do while
     252            ! afin de faire moins de calcul dans la boucle
     253            !--------------------------------------------------
    255254            do while (zlmelup>zltup)
    256                lt=lt+1
    257                zlt=zlev(ig,lt)
    258                zdz3=zlev(ig,lt+1)-zlt
    259                zltdwn=zlt-zdz3/2
    260                zltup=zlt+zdz3/2       
     255              lt = lt + 1
     256              zlt = zlev(ig, lt)
     257              zdz3 = zlev(ig, lt + 1) - zlt
     258              zltdwn = zlt - zdz3 / 2
     259              zltup = zlt + zdz3 / 2
    261260            enddo
    262 !--------------------------------------------------
    263 !AJ052014: Si iflag_thermals_ed<8 (par ex 6), alors
    264 ! on cherche o?? se trouve l'altitude d'inversion
    265 ! en calculant ztv1 (interpolation de la valeur de
    266 ! theta au niveau lt en utilisant les niveaux lt-1 et
    267 ! lt-2) et ztv2 (interpolation avec les niveaux lt+1
    268 ! et lt+2). Si theta r??ellement calcul??e au niveau lt
    269 ! comprise entre ztv1 et ztv2, alors il y a inversion
    270 ! et on calcule son altitude zinv en supposant que ztv(lt)
    271 ! est une combinaison lineaire de ztv1 et ztv2.
    272 ! Ensuite, on calcule la flottabilite en comparant
    273 ! la temperature de la couche l a celle de l'air situe
    274 ! l+lmel plus haut, ce qui necessite de savoir quel fraction
    275 ! de cet air est au-dessus ou en-dessous de l'inversion   
    276 !--------------------------------------------------
    277             atv1=(ztv(ig,lt-1)-ztv(ig,lt-2))/(zlev(ig,lt-1)-zlev(ig,lt-2))
    278             btv1=(ztv(ig,lt-2)*zlev(ig,lt-1)-ztv(ig,lt-1)*zlev(ig,lt-2)) &
    279             /(zlev(ig,lt-1)-zlev(ig,lt-2))
    280             atv2=(ztv(ig,lt+2)-ztv(ig,lt+1))/(zlev(ig,lt+2)-zlev(ig,lt+1))
    281             btv2=(ztv(ig,lt+1)*zlev(ig,lt+2)-ztv(ig,lt+2)*zlev(ig,lt+1)) &
    282             /(zlev(ig,lt+2)-zlev(ig,lt+1))
    283 
    284              ztv1=atv1*zlt+btv1
    285              ztv2=atv2*zlt+btv2
    286 
    287              if (ztv(ig,lt)>ztv1.and.ztv(ig,lt)<ztv2) then
    288 
    289 !--------------------------------------------------
    290 !AJ052014: D??calage de zinv qui est entre le haut
    291 !          et le bas de la couche lt
    292 !--------------------------------------------------
    293                 factinv=(ztv2-ztv(ig,lt))/(ztv2-ztv1)
    294                 zinv=zltdwn+zdz3*factinv
    295 
    296          
    297                 if (zlmeldwn>=zinv) then
    298                    ztv_est(ig,l)=atv2*zlmel+btv2
    299                    zbuoyjam(ig,l)=fact_shell*RG*(ztva_est(ig,l)-ztv_est(ig,l))/ztv_est(ig,l) &
    300                       +(1.-fact_shell)*zbuoy(ig,l)
    301                 elseif (zlmelup>=zinv) then
    302                  ztv_est2=atv2*0.5*(zlmelup+zinv)+btv2
    303                    ztv_est1=atv1*0.5*(zinv+zlmeldwn)+btv1
    304                    ztv_est(ig,l)=((zlmelup-zinv)/zdz)*ztv_est2+((zinv-zlmeldwn)/zdz)*ztv_est1
    305 
    306                    zbuoyjam(ig,l)=fact_shell*RG*(((zlmelup-zinv)/zdz)*(ztva_est(ig,l)- &
    307               ztv_est2)/ztv_est2+((zinv-zlmeldwn)/zdz)*(ztva_est(ig,l)- &
    308               ztv_est1)/ztv_est1)+(1.-fact_shell)*zbuoy(ig,l)
    309 
    310                 else
    311                    ztv_est(ig,l)=atv1*zlmel+btv1
    312                    zbuoyjam(ig,l)=fact_shell*RG*(ztva_est(ig,l)-ztv_est(ig,l))/ztv_est(ig,l) &
    313                              +(1.-fact_shell)*zbuoy(ig,l)
    314                 endif
    315 
    316              else ! if (ztv(ig,lt).gt.ztv1.and.ztv(ig,lt).lt.ztv2) then
    317 
    318                 if (zlmeldwn>zltdwn) then
    319                    zbuoyjam(ig,l)=fact_shell*RG*((ztva_est(ig,l)- &
    320                   ztv(ig,lt))/ztv(ig,lt))+(1.-fact_shell)*zbuoy(ig,l)
    321                 else
    322                    zbuoyjam(ig,l)=fact_shell*RG*(((zlmelup-zltdwn)/zdz)*(ztva_est(ig,l)- &
    323                   ztv(ig,lt))/ztv(ig,lt)+((zltdwn-zlmeldwn)/zdz)*(ztva_est(ig,l)- &
    324                   ztv(ig,lt-1))/ztv(ig,lt-1))+(1.-fact_shell)*zbuoy(ig,l)
    325 
    326                 endif
    327 
    328 !          zbuoyjam(ig,l)=fact_shell*RG*(((zlmelup-zltdwn)/zdz)*(ztva_est(ig,l)- &
    329 !    &          ztv1)/ztv1+((zltdwn-zlmeldwn)/zdz)*(ztva_est(ig,l)- &
    330 !    &          ztv(ig,lt-1))/ztv(ig,lt-1))+(1.-fact_shell)*zbuoy(ig,l)
    331 !         zdqt(ig,l)=Max(0.,((lmel+zdz3-zdz2)/zdz3)*(zqta(ig,l-1)- &
    332 !    &          po(ig,lt))/po(ig,lt)+((zdz2-lmel)/zdz3)*(zqta(ig,l-1)- &
    333 !     &          po(ig,lt-1))/po(ig,lt-1))
    334           endif ! if (ztv(ig,lt).gt.ztv1.and.ztv(ig,lt).lt.ztv2) then
    335 
    336         else  !   if (iflag_thermals_ed.lt.8) then
    337            lt=l+1
    338            zlt=zlev(ig,lt)
    339            zdz2=zlev(ig,lt)-zlev(ig,l)
    340 
    341            do while (lmel>zdz2)
    342              lt=lt+1
    343              zlt=zlev(ig,lt)
    344              zdz2=zlt-zlev(ig,l)
    345            enddo
    346            zdz3=zlev(ig,lt+1)-zlt
    347            zltdwn=zlev(ig,lt)-zdz3/2
    348            zlmelup=zlmel+(zdz/2)
    349            coefzlmel=Min(1.,(zlmelup-zltdwn)/zdz)
    350            zbuoyjam(ig,l)=1.*RG*(coefzlmel*(ztva_est(ig,l)- &
    351             ztv(ig,lt))/ztv(ig,lt)+(1.-coefzlmel)*(ztva_est(ig,l)- &
    352             ztv(ig,lt-1))/ztv(ig,lt-1))+0.*zbuoy(ig,l)
    353         endif !   if (iflag_thermals_ed.lt.8) then
    354 
    355 !------------------------------------------------
    356 !AJAM:nouveau calcul de w? 
    357 !------------------------------------------------
    358               zdz=zlev(ig,l+1)-zlev(ig,l)
    359               zdzbis=zlev(ig,l)-zlev(ig,l-1)
    360               zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
    361 
    362               zw2fact=fact_epsilon*2.*zdz/(1.+betalpha)
    363               zw2factbis=fact_epsilon*2.*zdzbis/(1.+betalpha)
    364               zdw2=afact*zbuoy(ig,l)/fact_epsilon
    365               zdw2bis=afact*zbuoy(ig,l-1)/fact_epsilon
    366 !              zdw2bis=0.5*(zdw2+zdw2bis)
    367               lm=Max(1,l-2)
    368 !              zdw2=(afact/fact_epsilon)*((zdz/zdzbis)*zbuoy(ig,l) &
    369 !    &              +((zdzbis-zdz)/zdzbis)*zbuoy(ig,l-1))
    370 !              zdw2bis=(afact/fact_epsilon)*((zdz/zdzbis)*zbuoy(ig,l-1) &
    371 !    &              +((zdzbis-zdz)/zdzbis)*zbuoy(ig,l-1))
    372 !             w_est(ig,l+1)=Max(0.0001,exp(-zw2fact)*(w_est(ig,l)-zdw2)+zdw2)
    373 !             w_est(ig,l+1)=(zdz/zdzbis)*Max(0.0001,exp(-zw2fact)* &
    374 !    &                     (w_est(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
    375 !    &                     Max(0.0001,exp(-zw2factbis)*(w_est(ig,l-1)-zdw2bis)+zdw2)
    376 !              w_est(ig,l+1)=Max(0.0001,(1-exp(-zw2fact))*zdw2+w_est(ig,l)*exp(-zw2fact))
    377 
    378 !--------------------------------------------------
    379 !AJ052014: J'ai remplac? w_est(ig,l) par zw2(ig,l)
    380 !--------------------------------------------------
    381          if (iflag_thermals_ed==8) then
    382 ! Ancienne version
    383 !             w_est(ig,l+1)=Max(0.0001,(zdz/zdzbis)*(exp(-zw2fact)* &
    384 !    &                     (w_est(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
    385 !    &                     (exp(-zw2factbis)*(w_est(ig,l-1)-zdw2bis)+zdw2))
    386 
    387             w_est(ig,l+1)=Max(0.0001,exp(-zw2fact)*(w_est(ig,l)-zdw2)+zdw2)
    388 
    389 ! Nouvelle version Arnaud
    390          else
    391 !             w_est(ig,l+1)=Max(0.0001,(zdz/zdzbis)*(exp(-zw2fact)* &
    392 !    &                     (w_est(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
    393 !    &                     (exp(-zw2factbis)*(w_est(ig,l-1)-zdw2bis)+zdw2))
    394 
    395             w_est(ig,l+1)=Max(0.0001,exp(-zw2fact)*(w_est(ig,l)-zdw2bis)+zdw2)
    396 
    397 !             w_est(ig,l+1)=Max(0.0001,(zdz/(zdzbis+zdz))*(exp(-zw2fact)* &
    398 !    &                     (w_est(ig,l)-zdw2bis)+zdw2)+(zdzbis/(zdzbis+zdz))* &
    399 !    &                     (exp(-zw2factbis)*(w_est(ig,l-1)-zdw2bis)+zdw2bis))
    400 
    401 
    402 
    403 !            w_est(ig,l+1)=Max(0.0001,(w_est(ig,l)+zdw2bis*zw2fact)*exp(-zw2fact))
    404 
    405 !             w_est(ig,l+1)=Max(0.0001,(zdz/zdzbis)*(zw2(ig,l)+zdw2*zw2fact)*exp(-zw2fact)+ &
    406 !    &                      (zdzbis-zdz)/zdzbis*(zw2(ig,l-1)+zdw2bis*zw2factbis)*exp(-zw2factbis))
    407 
    408 !             w_est(ig,l+1)=Max(0.0001,exp(-zw2factbis)*(w_est(ig,l-1)-zdw2bis)+zdw2)
    409 
    410          endif
    411 
    412 
    413          if (iflag_thermals_ed<6) then
    414              zalpha=f0(ig)*f_star(ig,l)/sqrt(w_est(ig,l+1))/rhobarz(ig,l)
    415 !              fact_epsilon=0.0005/(zalpha+0.025)**0.5
    416 !              fact_epsilon=Min(0.003,0.0004/(zalpha)**0.5)
    417               fact_epsilon=0.0002/(zalpha+0.1)
    418               zw2fact=fact_epsilon*2.*zdz/(1.+betalpha)
    419               zw2factbis=fact_epsilon*2.*zdzbis/(1.+betalpha)
    420               zdw2=afact*zbuoy(ig,l)/fact_epsilon
    421               zdw2bis=afact*zbuoy(ig,l-1)/fact_epsilon
    422 !              w_est(ig,l+1)=Max(0.0001,(zw2(ig,l)+zdw2*zw2fact)*exp(-zw2fact))
    423 
    424 !              w_est(ig,l+1)=Max(0.0001,(zdz/zdzbis)*(exp(-zw2fact)* &
    425 !    &                     (zw2(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
    426 !    &                     (exp(-zw2factbis)*(zw2(ig,l-1)-zdw2bis)+zdw2))
    427 
    428             w_est(ig,l+1)=Max(0.0001,exp(-zw2fact)*(w_est(ig,l)-zdw2bis)+zdw2)
    429 
    430 
    431          endif
    432 !--------------------------------------------------
    433 !AJ052014: J'ai comment? ce if plus n?cessaire puisqu'
    434 !on fait max(0.0001,.....)
    435 !--------------------------------------------------         
    436 
    437 !             if (w_est(ig,l+1).lt.0.) then
    438 !               w_est(ig,l+1)=zw2(ig,l)
    439 !                w_est(ig,l+1)=0.0001
    440 !             endif
    441 
    442        endif
     261            !--------------------------------------------------
     262            !AJ052014: Si iflag_thermals_ed<8 (par ex 6), alors
     263            ! on cherche o?? se trouve l'altitude d'inversion
     264            ! en calculant ztv1 (interpolation de la valeur de
     265            ! theta au niveau lt en utilisant les niveaux lt-1 et
     266            ! lt-2) et ztv2 (interpolation avec les niveaux lt+1
     267            ! et lt+2). Si theta r??ellement calcul??e au niveau lt
     268            ! comprise entre ztv1 et ztv2, alors il y a inversion
     269            ! et on calcule son altitude zinv en supposant que ztv(lt)
     270            ! est une combinaison lineaire de ztv1 et ztv2.
     271            ! Ensuite, on calcule la flottabilite en comparant
     272            ! la temperature de la couche l a celle de l'air situe
     273            ! l+lmel plus haut, ce qui necessite de savoir quel fraction
     274            ! de cet air est au-dessus ou en-dessous de l'inversion
     275            !--------------------------------------------------
     276            atv1 = (ztv(ig, lt - 1) - ztv(ig, lt - 2)) / (zlev(ig, lt - 1) - zlev(ig, lt - 2))
     277            btv1 = (ztv(ig, lt - 2) * zlev(ig, lt - 1) - ztv(ig, lt - 1) * zlev(ig, lt - 2)) &
     278                    / (zlev(ig, lt - 1) - zlev(ig, lt - 2))
     279            atv2 = (ztv(ig, lt + 2) - ztv(ig, lt + 1)) / (zlev(ig, lt + 2) - zlev(ig, lt + 1))
     280            btv2 = (ztv(ig, lt + 1) * zlev(ig, lt + 2) - ztv(ig, lt + 2) * zlev(ig, lt + 1)) &
     281                    / (zlev(ig, lt + 2) - zlev(ig, lt + 1))
     282
     283            ztv1 = atv1 * zlt + btv1
     284            ztv2 = atv2 * zlt + btv2
     285
     286            if (ztv(ig, lt)>ztv1.and.ztv(ig, lt)<ztv2) then
     287
     288              !--------------------------------------------------
     289              !AJ052014: D??calage de zinv qui est entre le haut
     290              !          et le bas de la couche lt
     291              !--------------------------------------------------
     292              factinv = (ztv2 - ztv(ig, lt)) / (ztv2 - ztv1)
     293              zinv = zltdwn + zdz3 * factinv
     294
     295              if (zlmeldwn>=zinv) then
     296                ztv_est(ig, l) = atv2 * zlmel + btv2
     297                zbuoyjam(ig, l) = fact_shell * RG * (ztva_est(ig, l) - ztv_est(ig, l)) / ztv_est(ig, l) &
     298                        + (1. - fact_shell) * zbuoy(ig, l)
     299              elseif (zlmelup>=zinv) then
     300                ztv_est2 = atv2 * 0.5 * (zlmelup + zinv) + btv2
     301                ztv_est1 = atv1 * 0.5 * (zinv + zlmeldwn) + btv1
     302                ztv_est(ig, l) = ((zlmelup - zinv) / zdz) * ztv_est2 + ((zinv - zlmeldwn) / zdz) * ztv_est1
     303
     304                zbuoyjam(ig, l) = fact_shell * RG * (((zlmelup - zinv) / zdz) * (ztva_est(ig, l) - &
     305                        ztv_est2) / ztv_est2 + ((zinv - zlmeldwn) / zdz) * (ztva_est(ig, l) - &
     306                        ztv_est1) / ztv_est1) + (1. - fact_shell) * zbuoy(ig, l)
     307
     308              else
     309                ztv_est(ig, l) = atv1 * zlmel + btv1
     310                zbuoyjam(ig, l) = fact_shell * RG * (ztva_est(ig, l) - ztv_est(ig, l)) / ztv_est(ig, l) &
     311                        + (1. - fact_shell) * zbuoy(ig, l)
     312              endif
     313
     314            else ! if (ztv(ig,lt).gt.ztv1.and.ztv(ig,lt).lt.ztv2) then
     315
     316              if (zlmeldwn>zltdwn) then
     317                zbuoyjam(ig, l) = fact_shell * RG * ((ztva_est(ig, l) - &
     318                        ztv(ig, lt)) / ztv(ig, lt)) + (1. - fact_shell) * zbuoy(ig, l)
     319              else
     320                zbuoyjam(ig, l) = fact_shell * RG * (((zlmelup - zltdwn) / zdz) * (ztva_est(ig, l) - &
     321                        ztv(ig, lt)) / ztv(ig, lt) + ((zltdwn - zlmeldwn) / zdz) * (ztva_est(ig, l) - &
     322                        ztv(ig, lt - 1)) / ztv(ig, lt - 1)) + (1. - fact_shell) * zbuoy(ig, l)
     323
     324              endif
     325
     326              !          zbuoyjam(ig,l)=fact_shell*RG*(((zlmelup-zltdwn)/zdz)*(ztva_est(ig,l)- &
     327              !    &          ztv1)/ztv1+((zltdwn-zlmeldwn)/zdz)*(ztva_est(ig,l)- &
     328              !    &          ztv(ig,lt-1))/ztv(ig,lt-1))+(1.-fact_shell)*zbuoy(ig,l)
     329              !         zdqt(ig,l)=Max(0.,((lmel+zdz3-zdz2)/zdz3)*(zqta(ig,l-1)- &
     330              !    &          po(ig,lt))/po(ig,lt)+((zdz2-lmel)/zdz3)*(zqta(ig,l-1)- &
     331              !     &          po(ig,lt-1))/po(ig,lt-1))
     332            endif ! if (ztv(ig,lt).gt.ztv1.and.ztv(ig,lt).lt.ztv2) then
     333
     334          else  !   if (iflag_thermals_ed.lt.8) then
     335            lt = l + 1
     336            zlt = zlev(ig, lt)
     337            zdz2 = zlev(ig, lt) - zlev(ig, l)
     338
     339            do while (lmel>zdz2)
     340              lt = lt + 1
     341              zlt = zlev(ig, lt)
     342              zdz2 = zlt - zlev(ig, l)
     343            enddo
     344            zdz3 = zlev(ig, lt + 1) - zlt
     345            zltdwn = zlev(ig, lt) - zdz3 / 2
     346            zlmelup = zlmel + (zdz / 2)
     347            coefzlmel = Min(1., (zlmelup - zltdwn) / zdz)
     348            zbuoyjam(ig, l) = 1. * RG * (coefzlmel * (ztva_est(ig, l) - &
     349                    ztv(ig, lt)) / ztv(ig, lt) + (1. - coefzlmel) * (ztva_est(ig, l) - &
     350                    ztv(ig, lt - 1)) / ztv(ig, lt - 1)) + 0. * zbuoy(ig, l)
     351          endif !   if (iflag_thermals_ed.lt.8) then
     352
     353          !------------------------------------------------
     354          !AJAM:nouveau calcul de w?
     355          !------------------------------------------------
     356          zdz = zlev(ig, l + 1) - zlev(ig, l)
     357          zdzbis = zlev(ig, l) - zlev(ig, l - 1)
     358          zbuoy(ig, l) = RG * (ztva_est(ig, l) - ztv(ig, l)) / ztv(ig, l)
     359
     360          zw2fact = fact_epsilon * 2. * zdz / (1. + betalpha)
     361          zw2factbis = fact_epsilon * 2. * zdzbis / (1. + betalpha)
     362          zdw2 = afact * zbuoy(ig, l) / fact_epsilon
     363          zdw2bis = afact * zbuoy(ig, l - 1) / fact_epsilon
     364          !              zdw2bis=0.5*(zdw2+zdw2bis)
     365          lm = Max(1, l - 2)
     366          !              zdw2=(afact/fact_epsilon)*((zdz/zdzbis)*zbuoy(ig,l) &
     367          !    &              +((zdzbis-zdz)/zdzbis)*zbuoy(ig,l-1))
     368          !              zdw2bis=(afact/fact_epsilon)*((zdz/zdzbis)*zbuoy(ig,l-1) &
     369          !    &              +((zdzbis-zdz)/zdzbis)*zbuoy(ig,l-1))
     370          !             w_est(ig,l+1)=Max(0.0001,exp(-zw2fact)*(w_est(ig,l)-zdw2)+zdw2)
     371          !             w_est(ig,l+1)=(zdz/zdzbis)*Max(0.0001,exp(-zw2fact)* &
     372          !    &                     (w_est(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
     373          !    &                     Max(0.0001,exp(-zw2factbis)*(w_est(ig,l-1)-zdw2bis)+zdw2)
     374          !              w_est(ig,l+1)=Max(0.0001,(1-exp(-zw2fact))*zdw2+w_est(ig,l)*exp(-zw2fact))
     375
     376          !--------------------------------------------------
     377          !AJ052014: J'ai remplac? w_est(ig,l) par zw2(ig,l)
     378          !--------------------------------------------------
     379          if (iflag_thermals_ed==8) then
     380            ! Ancienne version
     381            !             w_est(ig,l+1)=Max(0.0001,(zdz/zdzbis)*(exp(-zw2fact)* &
     382            !    &                     (w_est(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
     383            !    &                     (exp(-zw2factbis)*(w_est(ig,l-1)-zdw2bis)+zdw2))
     384
     385            w_est(ig, l + 1) = Max(0.0001, exp(-zw2fact) * (w_est(ig, l) - zdw2) + zdw2)
     386
     387            ! Nouvelle version Arnaud
     388          else
     389            !             w_est(ig,l+1)=Max(0.0001,(zdz/zdzbis)*(exp(-zw2fact)* &
     390            !    &                     (w_est(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
     391            !    &                     (exp(-zw2factbis)*(w_est(ig,l-1)-zdw2bis)+zdw2))
     392
     393            w_est(ig, l + 1) = Max(0.0001, exp(-zw2fact) * (w_est(ig, l) - zdw2bis) + zdw2)
     394
     395            !             w_est(ig,l+1)=Max(0.0001,(zdz/(zdzbis+zdz))*(exp(-zw2fact)* &
     396            !    &                     (w_est(ig,l)-zdw2bis)+zdw2)+(zdzbis/(zdzbis+zdz))* &
     397            !    &                     (exp(-zw2factbis)*(w_est(ig,l-1)-zdw2bis)+zdw2bis))
     398
     399
     400
     401            !            w_est(ig,l+1)=Max(0.0001,(w_est(ig,l)+zdw2bis*zw2fact)*exp(-zw2fact))
     402
     403            !             w_est(ig,l+1)=Max(0.0001,(zdz/zdzbis)*(zw2(ig,l)+zdw2*zw2fact)*exp(-zw2fact)+ &
     404            !    &                      (zdzbis-zdz)/zdzbis*(zw2(ig,l-1)+zdw2bis*zw2factbis)*exp(-zw2factbis))
     405
     406            !             w_est(ig,l+1)=Max(0.0001,exp(-zw2factbis)*(w_est(ig,l-1)-zdw2bis)+zdw2)
     407
     408          endif
     409
     410          if (iflag_thermals_ed<6) then
     411            zalpha = f0(ig) * f_star(ig, l) / sqrt(w_est(ig, l + 1)) / rhobarz(ig, l)
     412            !              fact_epsilon=0.0005/(zalpha+0.025)**0.5
     413            !              fact_epsilon=Min(0.003,0.0004/(zalpha)**0.5)
     414            fact_epsilon = 0.0002 / (zalpha + 0.1)
     415            zw2fact = fact_epsilon * 2. * zdz / (1. + betalpha)
     416            zw2factbis = fact_epsilon * 2. * zdzbis / (1. + betalpha)
     417            zdw2 = afact * zbuoy(ig, l) / fact_epsilon
     418            zdw2bis = afact * zbuoy(ig, l - 1) / fact_epsilon
     419            !              w_est(ig,l+1)=Max(0.0001,(zw2(ig,l)+zdw2*zw2fact)*exp(-zw2fact))
     420
     421            !              w_est(ig,l+1)=Max(0.0001,(zdz/zdzbis)*(exp(-zw2fact)* &
     422            !    &                     (zw2(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
     423            !    &                     (exp(-zw2factbis)*(zw2(ig,l-1)-zdw2bis)+zdw2))
     424
     425            w_est(ig, l + 1) = Max(0.0001, exp(-zw2fact) * (w_est(ig, l) - zdw2bis) + zdw2)
     426
     427          endif
     428          !--------------------------------------------------
     429          !AJ052014: J'ai comment? ce if plus n?cessaire puisqu'
     430          !on fait max(0.0001,.....)
     431          !--------------------------------------------------
     432
     433          !             if (w_est(ig,l+1).lt.0.) then
     434          !               w_est(ig,l+1)=zw2(ig,l)
     435          !                w_est(ig,l+1)=0.0001
     436          !             endif
     437
     438        endif
     439      enddo
     440
     441
     442      !-------------------------------------------------
     443      !calcul des taux d'entrainement et de detrainement
     444      !-------------------------------------------------
     445
     446      do ig = 1, ngrid
     447        if (active(ig)) then
     448
     449          !          zw2m=max(0.5*(w_est(ig,l)+w_est(ig,l+1)),0.1)
     450          zw2m = w_est(ig, l + 1)
     451          !          zw2m=zw2(ig,l)
     452          zdz = zlev(ig, l + 1) - zlev(ig, l)
     453          zbuoy(ig, l) = RG * (ztva_est(ig, l) - ztv(ig, l)) / ztv(ig, l)
     454          !          zbuoybis=zbuoy(ig,l)+RG*0.1/300.
     455          zbuoybis = zbuoy(ig, l)
     456          zalpha = f0(ig) * f_star(ig, l) / sqrt(w_est(ig, l + 1)) / rhobarz(ig, l)
     457          zdqt(ig, l) = max(zqta(ig, l - 1) - po(ig, l), 0.) / po(ig, l)
     458
     459
     460          !          entr_star(ig,l)=f_star(ig,l)*zdz*zbetalpha*MAX(0.,  &
     461          !    &     afact*zbuoybis/zw2m - fact_epsilon )
     462
     463          !          entr_star(ig,l)=MAX(0.,f_star(ig,l)*zdz*zbetalpha*  &
     464          !    &     afact*zbuoybis/zw2m - fact_epsilon )
     465
     466
     467
     468          !          zbuoyjam(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
     469
     470          !=========================================================================
     471          ! 4. Calcul de l'entrainement et du detrainement
     472          !=========================================================================
     473
     474          !          entr_star(ig,l)=f_star(ig,l)*zdz*zbetalpha*MAX(0.,  &
     475          !    &     afact*zbuoyjam(ig,l)/zw2m - fact_epsilon )
     476          !          entrbis=entr_star(ig,l)
     477
     478          if (iflag_thermals_ed<6) then
     479            fact_epsilon = 0.0002 / (zalpha + 0.1)
     480          endif
     481
     482          detr_star(ig, l) = f_star(ig, l) * zdz             &
     483                  * (mix0 * 0.1 / (zalpha + 0.001)               &
     484                          + MAX(detr_min, -afact * zbetalpha * zbuoyjam(ig, l) / zw2m   &
     485                                  + detr_q_coef * (zdqt(ig, l) / zw2m)**detr_q_power))
     486
     487          !          detr_star(ig,l)=(zdz/zdzbis)*detr_star(ig,l)+ &
     488          !    &                          ((zdzbis-zdz)/zdzbis)*detr_star(ig,l-1)
     489
     490          zbuoy(ig, l) = RG * (ztva_est(ig, l) - ztv(ig, l)) / ztv(ig, l)
     491
     492          entr_star(ig, l) = f_star(ig, l) * zdz * (&
     493                  mix0 * 0.1 / (zalpha + 0.001)               &
     494                          + zbetalpha * MAX(entr_min, &
     495                          afact * zbuoyjam(ig, l) / zw2m - fact_epsilon))
     496
     497
     498          !          entr_star(ig,l)=f_star(ig,l)*zdz* (         &
     499          !    &       mix0 * 0.1 / (zalpha+0.001)               &
     500          !    &     + MAX(entr_min,                   &
     501          !    &     zbetalpha*afact*zbuoyjam(ig,l)/zw2m - fact_epsilon +  &
     502          !    &     detr_q_coef*(zdqt(ig,l)/zw2m)**detr_q_power))
     503
     504
     505          !          entr_star(ig,l)=(zdz/zdzbis)*entr_star(ig,l)+ &
     506          !    &                          ((zdzbis-zdz)/zdzbis)*entr_star(ig,l-1)
     507
     508          !          entr_star(ig,l)=Max(0.,f_star(ig,l)*zdz*zbetalpha*  &
     509          !    &     afact*zbuoy(ig,l)/zw2m &
     510          !    &     - 1.*fact_epsilon)
     511
     512
     513          ! En dessous de lalim, on prend le max de alim_star et entr_star pour
     514          ! alim_star et 0 sinon
     515          if (l<lalim(ig)) then
     516            alim_star(ig, l) = max(alim_star(ig, l), entr_star(ig, l))
     517            entr_star(ig, l) = 0.
     518          endif
     519          !        if (l.lt.lalim(ig).and.alim_star(ig,l)>alim_star(ig,l-1)) then
     520          !          alim_star(ig,l)=entrbis
     521          !        endif
     522
     523          !        print*,'alim0',zlev(ig,l),entr_star(ig,l),detr_star(ig,l),zw2m,zbuoy(ig,l),f_star(ig,l)
     524          ! Calcul du flux montant normalise
     525          f_star(ig, l + 1) = f_star(ig, l) + alim_star(ig, l) + entr_star(ig, l)  &
     526                  - detr_star(ig, l)
     527
     528        endif
     529      enddo
     530
     531
     532      !============================================================================
     533      ! 5. calcul de la vitesse verticale en melangeant Tl et qt du thermique
     534      !===========================================================================
     535
     536      activetmp(:) = active(:) .and. f_star(:, l + 1)>1.e-10
     537      do ig = 1, ngrid
     538        if (activetmp(ig)) then
     539          Zsat = .false.
     540          ztla(ig, l) = (f_star(ig, l) * ztla(ig, l - 1) + &
     541                  (alim_star(ig, l) + entr_star(ig, l)) * zthl(ig, l))  &
     542                  / (f_star(ig, l + 1) + detr_star(ig, l))
     543          zqta(ig, l) = (f_star(ig, l) * zqta(ig, l - 1) + &
     544                  (alim_star(ig, l) + entr_star(ig, l)) * po(ig, l))  &
     545                  / (f_star(ig, l + 1) + detr_star(ig, l))
     546
     547        endif
     548      enddo
     549
     550      ztemp(:) = zpspsk(:, l) * ztla(:, l)
     551      call thermcell_qsat(ngrid, activetmp, pplev(:, l), ztemp, zqta(:, l), zqsatth(:, l))
     552      do ig = 1, ngrid
     553        if (activetmp(ig)) then
     554          ! on ecrit de maniere conservative (sat ou non)
     555          !          T = Tl +Lv/Cp ql
     556          zqla(ig, l) = max(0., zqta(ig, l) - zqsatth(ig, l))
     557          ztva(ig, l) = ztla(ig, l) * zpspsk(ig, l) + RLvCp * zqla(ig, l)
     558          ztva(ig, l) = ztva(ig, l) / zpspsk(ig, l)
     559          !on rajoute le calcul de zha pour diagnostiques (temp potentielle)
     560          zha(ig, l) = ztva(ig, l)
     561          ztva(ig, l) = ztva(ig, l) * (1. + RETV * (zqta(ig, l)  &
     562                  - zqla(ig, l)) - zqla(ig, l))
     563          zbuoy(ig, l) = RG * (ztva(ig, l) - ztv(ig, l)) / ztv(ig, l)
     564          zdz = zlev(ig, l + 1) - zlev(ig, l)
     565          zdzbis = zlev(ig, l) - zlev(ig, l - 1)
     566          zeps(ig, l) = (entr_star(ig, l) + alim_star(ig, l)) / (f_star(ig, l) * zdz)
     567          !!!!!!!          fact_epsilon=0.002
     568          zw2fact = fact_epsilon * 2. * zdz / (1. + betalpha)
     569          zw2factbis = fact_epsilon * 2. * zdzbis / (1. + betalpha)
     570          zdw2 = afact * zbuoy(ig, l) / (fact_epsilon)
     571          zdw2bis = afact * zbuoy(ig, l - 1) / (fact_epsilon)
     572          !              zdw2=(afact/fact_epsilon)*((zdz/zdzbis)*zbuoy(ig,l) &
     573          !    &              +((zdzbis-zdz)/zdzbis)*zbuoy(ig,l-1))
     574          !              lm=Max(1,l-2)
     575          !              zdw2bis=(afact/fact_epsilon)*((zdz/zdzbis)*zbuoy(ig,l-1) &
     576          !    &              +((zdzbis-zdz)/zdzbis)*zbuoy(ig,l-1))
     577          !            zw2(ig,l+1)=Max(0.0001,exp(-zw2fact)*(zw2(ig,l)-zdw2bis)+zdw2)
     578          !            zw2(ig,l+1)=Max(0.0001,(zdz/zdzbis)*(zw2(ig,l)+zdw2*zw2fact)*exp(-zw2fact)+ &
     579          !     &                   (zdzbis-zdz)/zdzbis*(zw2(ig,l-1)+zdw2bis*zw2factbis)*exp(-zw2factbis))
     580          !            zw2(ig,l+1)=Max(0.0001,(zw2(ig,l)+zdw2*zw2fact)*exp(-zw2fact))
     581          !            zw2(ig,l+1)=Max(0.0001,(zdz/zdzbis)*(exp(-zw2fact)* &
     582          !    &                     (zw2(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
     583          !    &                     (exp(-zw2factbis)*(zw2(ig,l-1)-zdw2bis)+zdw2))
     584          if (iflag_thermals_ed==8) then
     585            zw2(ig, l + 1) = Max(0.0001, exp(-zw2fact) * (zw2(ig, l) - zdw2) + zdw2)
     586          else
     587            zw2(ig, l + 1) = Max(0.0001, exp(-zw2fact) * (zw2(ig, l) - zdw2bis) + zdw2)
     588          endif
     589          !            zw2(ig,l+1)=Max(0.0001,(zdz/(zdz+zdzbis))*(exp(-zw2fact)* &
     590          !    &                     (zw2(ig,l)-zdw2)+zdw2bis)+(zdzbis/(zdz+zdzbis))* &
     591          !    &                     (exp(-zw2factbis)*(zw2(ig,l-1)-zdw2bis)+zdw2bis))
     592
     593          if (iflag_thermals_ed<6) then
     594            zalpha = f0(ig) * f_star(ig, l) / sqrt(zw2(ig, l + 1)) / rhobarz(ig, l)
     595            !           fact_epsilon=0.0005/(zalpha+0.025)**0.5
     596            !           fact_epsilon=Min(0.003,0.0004/(zalpha)**0.5)
     597            fact_epsilon = 0.0002 / (zalpha + 0.1)**1
     598            zw2fact = fact_epsilon * 2. * zdz / (1. + betalpha)
     599            zw2factbis = fact_epsilon * 2. * zdzbis / (1. + betalpha)
     600            zdw2 = afact * zbuoy(ig, l) / (fact_epsilon)
     601            zdw2bis = afact * zbuoy(ig, l - 1) / (fact_epsilon)
     602
     603            !            zw2(ig,l+1)=Max(0.0001,(zdz/zdzbis)*(exp(-zw2fact)* &
     604            !    &                     (zw2(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
     605            !    &                     (exp(-zw2factbis)*(zw2(ig,l-1)-zdw2bis)+zdw2))
     606            !            zw2(ig,l+1)=Max(0.0001,(zw2(ig,l)+zdw2*zw2fact)*exp(-zw2fact))
     607            zw2(ig, l + 1) = Max(0.0001, exp(-zw2fact) * (zw2(ig, l) - zdw2bis) + zdw2)
     608
     609          endif
     610
     611        endif
     612      enddo
     613
     614      if (prt_level>=20) print*, 'coucou calcul detr 460: ig, l', ig, l
     615
     616      !===========================================================================
     617      ! 6. initialisations pour le calcul de la hauteur du thermique, de l'inversion et de la vitesse verticale max
     618      !===========================================================================
     619
     620      nbpb = 0
     621      do ig = 1, ngrid
     622        if (zw2(ig, l + 1)>0. .and. zw2(ig, l + 1)<1.e-10) then
     623          !               stop'On tombe sur le cas particulier de thermcell_dry'
     624          !               print*,'On tombe sur le cas particulier de thermcell_plume'
     625          nbpb = nbpb + 1
     626          zw2(ig, l + 1) = 0.
     627          linter(ig) = l + 1
     628        endif
     629
     630        if (zw2(ig, l + 1)<0.) then
     631          linter(ig) = (l * (zw2(ig, l + 1) - zw2(ig, l))  &
     632                  - zw2(ig, l)) / (zw2(ig, l + 1) - zw2(ig, l))
     633          zw2(ig, l + 1) = 0.
     634          !+CR:04/05/12:correction calcul linter pour calcul de zmax continu
     635        elseif (f_star(ig, l + 1)<0.) then
     636          linter(ig) = (l * (f_star(ig, l + 1) - f_star(ig, l))  &
     637                  - f_star(ig, l)) / (f_star(ig, l + 1) - f_star(ig, l))
     638          zw2(ig, l + 1) = 0.
     639          !fin CR:04/05/12
     640        endif
     641
     642        wa_moy(ig, l + 1) = sqrt(zw2(ig, l + 1))
     643
     644        if (wa_moy(ig, l + 1)>wmaxa(ig)) then
     645          !   lmix est le niveau de la couche ou w (wa_moy) est maximum
     646          !on rajoute le calcul de lmix_bis
     647          if (zqla(ig, l)<1.e-10) then
     648            lmix_bis(ig) = l + 1
     649          endif
     650          lmix(ig) = l + 1
     651          wmaxa(ig) = wa_moy(ig, l + 1)
     652        endif
     653      enddo
     654
     655      if (nbpb>0) then
     656        print*, 'WARNING on tombe ', nbpb, ' x sur un pb pour l=', l, ' dans thermcell_plume'
     657      endif
     658
     659      !=========================================================================
     660      ! FIN DE LA BOUCLE VERTICALE
    443661    enddo
    444 
    445 
    446 !-------------------------------------------------
    447 !calcul des taux d'entrainement et de detrainement
    448 !-------------------------------------------------
    449 
    450      do ig=1,ngrid
     662    !=========================================================================
     663
     664    !on recalcule alim_star_tot
     665    do ig = 1, ngrid
     666      alim_star_tot(ig) = 0.
     667    enddo
     668    do ig = 1, ngrid
     669      do l = 1, lalim(ig) - 1
     670        alim_star_tot(ig) = alim_star_tot(ig) + alim_star(ig, l)
     671      enddo
     672    enddo
     673
     674    if (prt_level>=20) print*, 'coucou calcul detr 470: ig, l', ig, l
     675    RETURN
     676  end
     677
     678
     679  SUBROUTINE thermcell_plume_5B(itap, ngrid, nlay, ptimestep, ztv, zthl, po, zl, rhobarz, &
     680          zlev, pplev, pphi, zpspsk, alim_star, alim_star_tot, &
     681          lalim, f0, detr_star, entr_star, f_star, csc, ztva, &
     682          ztla, zqla, zqta, zha, zw2, w_est, ztva_est, zqsatth, lmix, lmix_bis, linter &
     683          , lev_out, lunout1, igout)
     684    !&           ,lev_out,lunout1,igout,zbuoy,zbuoyjam)
     685
     686    !--------------------------------------------------------------------------
     687    !thermcell_plume: calcule les valeurs de qt, thetal et w dans l ascendance
     688    ! Version conforme a l'article de Rio et al. 2010.
     689    ! Code ecrit par Catherine Rio, Arnaud Jam et Frederic Hourdin
     690    !--------------------------------------------------------------------------
     691
     692    USE lmdz_thermcell_ini, ONLY: prt_level, fact_thermals_ed_dz, iflag_thermals_ed, RLvCP, RETV, RG
     693    USE lmdz_thermcell_qsat, ONLY: thermcell_qsat
     694    IMPLICIT NONE
     695
     696    INTEGER itap
     697    INTEGER lunout1, igout
     698    INTEGER ngrid, nlay
     699    REAL ptimestep
     700    REAL ztv(ngrid, nlay)
     701    REAL zthl(ngrid, nlay)
     702    REAL, INTENT(IN) :: po(ngrid, nlay)
     703    REAL zl(ngrid, nlay)
     704    REAL rhobarz(ngrid, nlay)
     705    REAL zlev(ngrid, nlay + 1)
     706    REAL pplev(ngrid, nlay + 1)
     707    REAL pphi(ngrid, nlay)
     708    REAL zpspsk(ngrid, nlay)
     709    REAL alim_star(ngrid, nlay)
     710    REAL f0(ngrid)
     711    INTEGER lalim(ngrid)
     712    integer lev_out                           ! niveau pour les print
     713    integer nbpb
     714
     715    real alim_star_tot(ngrid)
     716
     717    REAL ztva(ngrid, nlay)
     718    REAL ztla(ngrid, nlay)
     719    REAL zqla(ngrid, nlay)
     720    REAL zqta(ngrid, nlay)
     721    REAL zha(ngrid, nlay)
     722
     723    REAL detr_star(ngrid, nlay)
     724    REAL coefc
     725    REAL entr_star(ngrid, nlay)
     726    REAL detr(ngrid, nlay)
     727    REAL entr(ngrid, nlay)
     728
     729    REAL csc(ngrid, nlay)
     730
     731    REAL zw2(ngrid, nlay + 1)
     732    REAL w_est(ngrid, nlay + 1)
     733    REAL f_star(ngrid, nlay + 1)
     734    REAL wa_moy(ngrid, nlay + 1)
     735
     736    REAL ztva_est(ngrid, nlay)
     737    REAL zqla_est(ngrid, nlay)
     738    REAL zqsatth(ngrid, nlay)
     739    REAL zta_est(ngrid, nlay)
     740    REAL zbuoyjam(ngrid, nlay)
     741    REAL ztemp(ngrid), zqsat(ngrid)
     742    REAL zdw2
     743    REAL zw2modif
     744    REAL zw2fact
     745    REAL zeps(ngrid, nlay)
     746
     747    REAL linter(ngrid)
     748    INTEGER lmix(ngrid)
     749    INTEGER lmix_bis(ngrid)
     750    REAL    wmaxa(ngrid)
     751
     752    INTEGER ig, l, k
     753
     754    real zdz, zbuoy(ngrid, nlay), zalpha, gamma(ngrid, nlay), zdqt(ngrid, nlay), zw2m
     755    real zbuoybis
     756    real zcor, zdelta, zcvm5, qlbef, zdz2
     757    real betalpha, zbetalpha
     758    real eps, afact
     759    logical Zsat
     760    LOGICAL active(ngrid), activetmp(ngrid)
     761    REAL fact_gamma, fact_epsilon, fact_gamma2, fact_epsilon2
     762    REAL c2(ngrid, nlay)
     763    Zsat = .false.
     764    ! Initialisation
     765
     766    fact_epsilon = 0.002
     767    betalpha = 0.9
     768    afact = 2. / 3.
     769
     770    zbetalpha = betalpha / (1. + betalpha)
     771
     772
     773    ! Initialisations des variables reeles
     774    if (1==1) then
     775      ztva(:, :) = ztv(:, :)
     776      ztva_est(:, :) = ztva(:, :)
     777      ztla(:, :) = zthl(:, :)
     778      zqta(:, :) = po(:, :)
     779      zha(:, :) = ztva(:, :)
     780    else
     781      ztva(:, :) = 0.
     782      ztva_est(:, :) = 0.
     783      ztla(:, :) = 0.
     784      zqta(:, :) = 0.
     785      zha(:, :) = 0.
     786    endif
     787
     788    zqla_est(:, :) = 0.
     789    zqsatth(:, :) = 0.
     790    zqla(:, :) = 0.
     791    detr_star(:, :) = 0.
     792    entr_star(:, :) = 0.
     793    alim_star(:, :) = 0.
     794    alim_star_tot(:) = 0.
     795    csc(:, :) = 0.
     796    detr(:, :) = 0.
     797    entr(:, :) = 0.
     798    zw2(:, :) = 0.
     799    zbuoy(:, :) = 0.
     800    zbuoyjam(:, :) = 0.
     801    gamma(:, :) = 0.
     802    zeps(:, :) = 0.
     803    w_est(:, :) = 0.
     804    f_star(:, :) = 0.
     805    wa_moy(:, :) = 0.
     806    linter(:) = 1.
     807    !     linter(:)=1.
     808    ! Initialisation des variables entieres
     809    lmix(:) = 1
     810    lmix_bis(:) = 2
     811    wmaxa(:) = 0.
     812    lalim(:) = 1
     813
     814
     815    !-------------------------------------------------------------------------
     816    ! On ne considere comme actif que les colonnes dont les deux premieres
     817    ! couches sont instables.
     818    !-------------------------------------------------------------------------
     819    active(:) = ztv(:, 1)>ztv(:, 2)
     820
     821    !-------------------------------------------------------------------------
     822    ! Definition de l'alimentation
     823    !-------------------------------------------------------------------------
     824    do l = 1, nlay - 1
     825      do ig = 1, ngrid
     826        if (ztv(ig, l)> ztv(ig, l + 1) .and. ztv(ig, 1)>=ztv(ig, l)) then
     827          alim_star(ig, l) = MAX((ztv(ig, l) - ztv(ig, l + 1)), 0.)  &
     828                  * sqrt(zlev(ig, l + 1))
     829          lalim(ig) = l + 1
     830          alim_star_tot(ig) = alim_star_tot(ig) + alim_star(ig, l)
     831        endif
     832      enddo
     833    enddo
     834    do l = 1, nlay
     835      do ig = 1, ngrid
     836        if (alim_star_tot(ig) > 1.e-10) then
     837          alim_star(ig, l) = alim_star(ig, l) / alim_star_tot(ig)
     838        endif
     839      enddo
     840    enddo
     841    alim_star_tot(:) = 1.
     842
     843
     844
     845    !------------------------------------------------------------------------------
     846    ! Calcul dans la premiere couche
     847    ! On decide dans cette version que le thermique n'est actif que si la premiere
     848    ! couche est instable.
     849    ! Pourrait etre change si on veut que le thermiques puisse se d??clencher
     850    ! dans une couche l>1
     851    !------------------------------------------------------------------------------
     852    do ig = 1, ngrid
     853      ! Le panache va prendre au debut les caracteristiques de l'air contenu
     854      ! dans cette couche.
     855      if (active(ig)) then
     856        ztla(ig, 1) = zthl(ig, 1)
     857        zqta(ig, 1) = po(ig, 1)
     858        zqla(ig, 1) = zl(ig, 1)
     859        !cr: attention, prise en compte de f*(1)=1
     860        f_star(ig, 2) = alim_star(ig, 1)
     861        zw2(ig, 2) = 2. * RG * (ztv(ig, 1) - ztv(ig, 2)) / ztv(ig, 2)  &
     862                & * (zlev(ig, 2) - zlev(ig, 1))  &
     863                & * 0.4 * pphi(ig, 1) / (pphi(ig, 2) - pphi(ig, 1))
     864        w_est(ig, 2) = zw2(ig, 2)
     865      endif
     866    enddo
     867
     868    !==============================================================================
     869    !boucle de calcul de la vitesse verticale dans le thermique
     870    !==============================================================================
     871    do l = 2, nlay - 1
     872      !==============================================================================
     873
     874
     875      ! On decide si le thermique est encore actif ou non
     876      ! AFaire : Il faut sans doute ajouter entr_star a alim_star dans ce test
     877      do ig = 1, ngrid
     878        active(ig) = active(ig) &
     879                &                 .and. zw2(ig, l)>1.e-10 &
     880                &                 .and. f_star(ig, l) + alim_star(ig, l)>1.e-10
     881      enddo
     882
     883
     884
     885      !---------------------------------------------------------------------------
     886      ! calcul des proprietes thermodynamiques et de la vitesse de la couche l
     887      ! sans tenir compte du detrainement et de l'entrainement dans cette
     888      ! couche
     889      ! C'est a dire qu'on suppose
     890      ! ztla(l)=ztla(l-1) et zqta(l)=zqta(l-1)
     891      ! Ici encore, on doit pouvoir ajouter entr_star (qui peut etre calculer
     892      ! avant) a l'alimentation pour avoir un calcul plus propre
     893      !---------------------------------------------------------------------------
     894
     895      ztemp(:) = zpspsk(:, l) * ztla(:, l - 1)
     896      call thermcell_qsat(ngrid, active, pplev(:, l), ztemp, zqta(:, l - 1), zqsat(:))
     897
     898      do ig = 1, ngrid
     899        !       print*,'active',active(ig),ig,l
     900        if(active(ig)) then
     901          zqla_est(ig, l) = max(0., zqta(ig, l - 1) - zqsat(ig))
     902          ztva_est(ig, l) = ztla(ig, l - 1) * zpspsk(ig, l) + RLvCp * zqla_est(ig, l)
     903          zta_est(ig, l) = ztva_est(ig, l)
     904          ztva_est(ig, l) = ztva_est(ig, l) / zpspsk(ig, l)
     905          ztva_est(ig, l) = ztva_est(ig, l) * (1. + RETV * (zqta(ig, l - 1)  &
     906                  - zqla_est(ig, l)) - zqla_est(ig, l))
     907
     908          !------------------------------------------------
     909          !AJAM:nouveau calcul de w?
     910          !------------------------------------------------
     911          zdz = zlev(ig, l + 1) - zlev(ig, l)
     912          zbuoy(ig, l) = RG * (ztva_est(ig, l) - ztv(ig, l)) / ztv(ig, l)
     913
     914          zw2fact = fact_epsilon * 2. * zdz / (1. + betalpha)
     915          zdw2 = (afact) * zbuoy(ig, l) / (fact_epsilon)
     916          w_est(ig, l + 1) = Max(0.0001, exp(-zw2fact) * (w_est(ig, l) - zdw2) + zdw2)
     917
     918          if (w_est(ig, l + 1)<0.) then
     919            w_est(ig, l + 1) = zw2(ig, l)
     920          endif
     921        endif
     922      enddo
     923
     924
     925      !-------------------------------------------------
     926      !calcul des taux d'entrainement et de detrainement
     927      !-------------------------------------------------
     928
     929      do ig = 1, ngrid
    451930        if (active(ig)) then
    452931
    453 !          zw2m=max(0.5*(w_est(ig,l)+w_est(ig,l+1)),0.1)
    454           zw2m=w_est(ig,l+1)
    455 !          zw2m=zw2(ig,l)
    456           zdz=zlev(ig,l+1)-zlev(ig,l)
    457           zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
    458 !          zbuoybis=zbuoy(ig,l)+RG*0.1/300.
    459           zbuoybis=zbuoy(ig,l)
    460           zalpha=f0(ig)*f_star(ig,l)/sqrt(w_est(ig,l+1))/rhobarz(ig,l)
    461           zdqt(ig,l)=max(zqta(ig,l-1)-po(ig,l),0.)/po(ig,l)
    462 
    463          
    464 !          entr_star(ig,l)=f_star(ig,l)*zdz*zbetalpha*MAX(0.,  &
    465 !    &     afact*zbuoybis/zw2m - fact_epsilon )
    466 
    467 !          entr_star(ig,l)=MAX(0.,f_star(ig,l)*zdz*zbetalpha*  &
    468 !    &     afact*zbuoybis/zw2m - fact_epsilon )
    469 
    470 
    471 
    472 !          zbuoyjam(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
    473 
    474 !=========================================================================
    475 ! 4. Calcul de l'entrainement et du detrainement
    476 !=========================================================================
    477 
    478 !          entr_star(ig,l)=f_star(ig,l)*zdz*zbetalpha*MAX(0.,  &
    479 !    &     afact*zbuoyjam(ig,l)/zw2m - fact_epsilon )
    480 !          entrbis=entr_star(ig,l)
    481 
    482           if (iflag_thermals_ed<6) then
    483           fact_epsilon=0.0002/(zalpha+0.1)
     932          zw2m = max(0.5 * (w_est(ig, l) + w_est(ig, l + 1)), 0.1)
     933          zw2m = w_est(ig, l + 1)
     934          zdz = zlev(ig, l + 1) - zlev(ig, l)
     935          zbuoy(ig, l) = RG * (ztva_est(ig, l) - ztv(ig, l)) / ztv(ig, l)
     936          !          zbuoybis=zbuoy(ig,l)+RG*0.1/300.
     937          zbuoybis = zbuoy(ig, l)
     938          zalpha = f0(ig) * f_star(ig, l) / sqrt(w_est(ig, l + 1)) / rhobarz(ig, l)
     939          zdqt(ig, l) = max(zqta(ig, l - 1) - po(ig, l), 0.) / po(ig, l)
     940
     941          entr_star(ig, l) = f_star(ig, l) * zdz * zbetalpha * MAX(0., &
     942                  afact * zbuoybis / zw2m - fact_epsilon)
     943
     944          detr_star(ig, l) = f_star(ig, l) * zdz                        &
     945                  * MAX(1.e-3, -afact * zbetalpha * zbuoy(ig, l) / zw2m          &
     946                          + 0.012 * (zdqt(ig, l) / zw2m)**0.5)
     947
     948          ! En dessous de lalim, on prend le max de alim_star et entr_star pour
     949          ! alim_star et 0 sinon
     950          if (l<lalim(ig)) then
     951            alim_star(ig, l) = max(alim_star(ig, l), entr_star(ig, l))
     952            entr_star(ig, l) = 0.
    484953          endif
    485          
    486 
    487 
    488           detr_star(ig,l)=f_star(ig,l)*zdz             &
    489        *( mix0 * 0.1 / (zalpha+0.001)               &
    490        + MAX(detr_min, -afact*zbetalpha*zbuoyjam(ig,l)/zw2m   &
    491        + detr_q_coef*(zdqt(ig,l)/zw2m)**detr_q_power))
    492 
    493 !          detr_star(ig,l)=(zdz/zdzbis)*detr_star(ig,l)+ &
    494 !    &                          ((zdzbis-zdz)/zdzbis)*detr_star(ig,l-1)
    495 
    496           zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
    497 
    498           entr_star(ig,l)=f_star(ig,l)*zdz* (         &
    499          mix0 * 0.1 / (zalpha+0.001)               &
    500        + zbetalpha*MAX(entr_min,                   &
    501        afact*zbuoyjam(ig,l)/zw2m - fact_epsilon))
    502 
    503 
    504 !          entr_star(ig,l)=f_star(ig,l)*zdz* (         &
    505 !    &       mix0 * 0.1 / (zalpha+0.001)               &
    506 !    &     + MAX(entr_min,                   &
    507 !    &     zbetalpha*afact*zbuoyjam(ig,l)/zw2m - fact_epsilon +  &
    508 !    &     detr_q_coef*(zdqt(ig,l)/zw2m)**detr_q_power))
    509 
    510 
    511 !          entr_star(ig,l)=(zdz/zdzbis)*entr_star(ig,l)+ &
    512 !    &                          ((zdzbis-zdz)/zdzbis)*entr_star(ig,l-1)
    513 
    514 !          entr_star(ig,l)=Max(0.,f_star(ig,l)*zdz*zbetalpha*  &     
    515 !    &     afact*zbuoy(ig,l)/zw2m &
    516 !    &     - 1.*fact_epsilon)
    517 
    518          
    519 ! En dessous de lalim, on prend le max de alim_star et entr_star pour
    520 ! alim_star et 0 sinon
    521         if (l<lalim(ig)) then
    522           alim_star(ig,l)=max(alim_star(ig,l),entr_star(ig,l))
    523           entr_star(ig,l)=0.
    524         endif
    525 !        if (l.lt.lalim(ig).and.alim_star(ig,l)>alim_star(ig,l-1)) then
    526 !          alim_star(ig,l)=entrbis
    527 !        endif
    528 
    529 !        print*,'alim0',zlev(ig,l),entr_star(ig,l),detr_star(ig,l),zw2m,zbuoy(ig,l),f_star(ig,l)
    530 ! Calcul du flux montant normalise
    531       f_star(ig,l+1)=f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l)  &
    532                 -detr_star(ig,l)
    533 
     954
     955          ! Calcul du flux montant normalise
     956          f_star(ig, l + 1) = f_star(ig, l) + alim_star(ig, l) + entr_star(ig, l)  &
     957                  - detr_star(ig, l)
     958
     959        endif
     960      enddo
     961
     962
     963      !----------------------------------------------------------------------------
     964      !calcul de la vitesse verticale en melangeant Tl et qt du thermique
     965      !---------------------------------------------------------------------------
     966      activetmp(:) = active(:) .and. f_star(:, l + 1)>1.e-10
     967      do ig = 1, ngrid
     968        if (activetmp(ig)) then
     969          Zsat = .false.
     970          ztla(ig, l) = (f_star(ig, l) * ztla(ig, l - 1) + &
     971                  (alim_star(ig, l) + entr_star(ig, l)) * zthl(ig, l))  &
     972                  / (f_star(ig, l + 1) + detr_star(ig, l))
     973          zqta(ig, l) = (f_star(ig, l) * zqta(ig, l - 1) + &
     974                  (alim_star(ig, l) + entr_star(ig, l)) * po(ig, l))  &
     975                  / (f_star(ig, l + 1) + detr_star(ig, l))
     976
     977        endif
     978      enddo
     979
     980      ztemp(:) = zpspsk(:, l) * ztla(:, l)
     981      call thermcell_qsat(ngrid, activetmp, pplev(:, l), ztemp, zqta(:, l), zqsatth(:, l))
     982
     983      do ig = 1, ngrid
     984        if (activetmp(ig)) then
     985          ! on ecrit de maniere conservative (sat ou non)
     986          !          T = Tl +Lv/Cp ql
     987          zqla(ig, l) = max(0., zqta(ig, l) - zqsatth(ig, l))
     988          ztva(ig, l) = ztla(ig, l) * zpspsk(ig, l) + RLvCp * zqla(ig, l)
     989          ztva(ig, l) = ztva(ig, l) / zpspsk(ig, l)
     990          !on rajoute le calcul de zha pour diagnostiques (temp potentielle)
     991          zha(ig, l) = ztva(ig, l)
     992          ztva(ig, l) = ztva(ig, l) * (1. + RETV * (zqta(ig, l)  &
     993                  - zqla(ig, l)) - zqla(ig, l))
     994          zbuoy(ig, l) = RG * (ztva(ig, l) - ztv(ig, l)) / ztv(ig, l)
     995          zdz = zlev(ig, l + 1) - zlev(ig, l)
     996          zeps(ig, l) = (entr_star(ig, l) + alim_star(ig, l)) / (f_star(ig, l) * zdz)
     997
     998          zw2fact = fact_epsilon * 2. * zdz / (1. + betalpha)
     999          zdw2 = afact * zbuoy(ig, l) / (fact_epsilon)
     1000          zw2(ig, l + 1) = Max(0.0001, exp(-zw2fact) * (zw2(ig, l) - zdw2) + zdw2)
     1001        endif
     1002      enddo
     1003
     1004      if (prt_level>=20) print*, 'coucou calcul detr 460: ig, l', ig, l
     1005
     1006      !---------------------------------------------------------------------------
     1007      !initialisations pour le calcul de la hauteur du thermique, de l'inversion et de la vitesse verticale max
     1008      !---------------------------------------------------------------------------
     1009
     1010      nbpb = 0
     1011      do ig = 1, ngrid
     1012        if (zw2(ig, l + 1)>0. .and. zw2(ig, l + 1)<1.e-10) then
     1013          !               stop'On tombe sur le cas particulier de thermcell_dry'
     1014          !               print*,'On tombe sur le cas particulier de thermcell_plume'
     1015          nbpb = nbpb + 1
     1016          zw2(ig, l + 1) = 0.
     1017          linter(ig) = l + 1
     1018        endif
     1019
     1020        if (zw2(ig, l + 1)<0.) then
     1021          linter(ig) = (l * (zw2(ig, l + 1) - zw2(ig, l))  &
     1022                  - zw2(ig, l)) / (zw2(ig, l + 1) - zw2(ig, l))
     1023          zw2(ig, l + 1) = 0.
     1024        elseif (f_star(ig, l + 1)<0.) then
     1025          linter(ig) = (l * (f_star(ig, l + 1) - f_star(ig, l))  &
     1026                  - f_star(ig, l)) / (f_star(ig, l + 1) - f_star(ig, l))
     1027          !           print*,"linter plume", linter(ig)
     1028          zw2(ig, l + 1) = 0.
     1029        endif
     1030
     1031        wa_moy(ig, l + 1) = sqrt(zw2(ig, l + 1))
     1032
     1033        if (wa_moy(ig, l + 1)>wmaxa(ig)) then
     1034          !   lmix est le niveau de la couche ou w (wa_moy) est maximum
     1035          !on rajoute le calcul de lmix_bis
     1036          if (zqla(ig, l)<1.e-10) then
     1037            lmix_bis(ig) = l + 1
     1038          endif
     1039          lmix(ig) = l + 1
     1040          wmaxa(ig) = wa_moy(ig, l + 1)
     1041        endif
     1042      enddo
     1043
     1044      if (nbpb>0) then
     1045        print*, 'WARNING on tombe ', nbpb, ' x sur un pb pour l=', l, ' dans thermcell_plume'
    5341046      endif
    535    enddo
    536 
    537 
    538 !============================================================================
    539 ! 5. calcul de la vitesse verticale en melangeant Tl et qt du thermique
    540 !===========================================================================
    541 
    542    activetmp(:)=active(:) .and. f_star(:,l+1)>1.e-10
    543    do ig=1,ngrid
    544        if (activetmp(ig)) then
    545            Zsat=.false.
    546            ztla(ig,l)=(f_star(ig,l)*ztla(ig,l-1)+  &
    547               (alim_star(ig,l)+entr_star(ig,l))*zthl(ig,l))  &
    548               /(f_star(ig,l+1)+detr_star(ig,l))
    549            zqta(ig,l)=(f_star(ig,l)*zqta(ig,l-1)+  &
    550               (alim_star(ig,l)+entr_star(ig,l))*po(ig,l))  &
    551               /(f_star(ig,l+1)+detr_star(ig,l))
    552 
    553         endif
     1047
     1048      !=========================================================================
     1049      ! FIN DE LA BOUCLE VERTICALE
    5541050    enddo
    555 
    556    ztemp(:)=zpspsk(:,l)*ztla(:,l)
    557    call thermcell_qsat(ngrid,activetmp,pplev(:,l),ztemp,zqta(:,l),zqsatth(:,l))
    558    do ig=1,ngrid
    559       if (activetmp(ig)) then
    560 ! on ecrit de maniere conservative (sat ou non)
    561 !          T = Tl +Lv/Cp ql
    562            zqla(ig,l)=max(0.,zqta(ig,l)-zqsatth(ig,l))
    563            ztva(ig,l) = ztla(ig,l)*zpspsk(ig,l)+RLvCp*zqla(ig,l)
    564            ztva(ig,l) = ztva(ig,l)/zpspsk(ig,l)
    565 !on rajoute le calcul de zha pour diagnostiques (temp potentielle)
    566            zha(ig,l) = ztva(ig,l)
    567            ztva(ig,l) = ztva(ig,l)*(1.+RETV*(zqta(ig,l)  &
    568                 -zqla(ig,l))-zqla(ig,l))
    569            zbuoy(ig,l)=RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)
    570            zdz=zlev(ig,l+1)-zlev(ig,l)
    571            zdzbis=zlev(ig,l)-zlev(ig,l-1)
    572            zeps(ig,l)=(entr_star(ig,l)+alim_star(ig,l))/(f_star(ig,l)*zdz)
    573 !!!!!!!          fact_epsilon=0.002
    574             zw2fact=fact_epsilon*2.*zdz/(1.+betalpha)
    575             zw2factbis=fact_epsilon*2.*zdzbis/(1.+betalpha)
    576             zdw2= afact*zbuoy(ig,l)/(fact_epsilon)
    577             zdw2bis= afact*zbuoy(ig,l-1)/(fact_epsilon)
    578 !              zdw2=(afact/fact_epsilon)*((zdz/zdzbis)*zbuoy(ig,l) &
    579 !    &              +((zdzbis-zdz)/zdzbis)*zbuoy(ig,l-1))
    580 !              lm=Max(1,l-2)
    581 !              zdw2bis=(afact/fact_epsilon)*((zdz/zdzbis)*zbuoy(ig,l-1) &
    582 !    &              +((zdzbis-zdz)/zdzbis)*zbuoy(ig,l-1))
    583 !            zw2(ig,l+1)=Max(0.0001,exp(-zw2fact)*(zw2(ig,l)-zdw2bis)+zdw2)
    584 !            zw2(ig,l+1)=Max(0.0001,(zdz/zdzbis)*(zw2(ig,l)+zdw2*zw2fact)*exp(-zw2fact)+ &
    585 !     &                   (zdzbis-zdz)/zdzbis*(zw2(ig,l-1)+zdw2bis*zw2factbis)*exp(-zw2factbis))
    586 !            zw2(ig,l+1)=Max(0.0001,(zw2(ig,l)+zdw2*zw2fact)*exp(-zw2fact))
    587 !            zw2(ig,l+1)=Max(0.0001,(zdz/zdzbis)*(exp(-zw2fact)* &
    588 !    &                     (zw2(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
    589 !    &                     (exp(-zw2factbis)*(zw2(ig,l-1)-zdw2bis)+zdw2))
    590             if (iflag_thermals_ed==8) then
    591             zw2(ig,l+1)=Max(0.0001,exp(-zw2fact)*(zw2(ig,l)-zdw2)+zdw2)
    592             else
    593             zw2(ig,l+1)=Max(0.0001,exp(-zw2fact)*(zw2(ig,l)-zdw2bis)+zdw2)
    594             endif
    595 !            zw2(ig,l+1)=Max(0.0001,(zdz/(zdz+zdzbis))*(exp(-zw2fact)* &
    596 !    &                     (zw2(ig,l)-zdw2)+zdw2bis)+(zdzbis/(zdz+zdzbis))* &
    597 !    &                     (exp(-zw2factbis)*(zw2(ig,l-1)-zdw2bis)+zdw2bis))
    598 
    599 
    600            if (iflag_thermals_ed<6) then
    601            zalpha=f0(ig)*f_star(ig,l)/sqrt(zw2(ig,l+1))/rhobarz(ig,l)
    602 !           fact_epsilon=0.0005/(zalpha+0.025)**0.5
    603 !           fact_epsilon=Min(0.003,0.0004/(zalpha)**0.5)
    604            fact_epsilon=0.0002/(zalpha+0.1)**1
    605             zw2fact=fact_epsilon*2.*zdz/(1.+betalpha)
    606             zw2factbis=fact_epsilon*2.*zdzbis/(1.+betalpha)
    607             zdw2= afact*zbuoy(ig,l)/(fact_epsilon)
    608             zdw2bis= afact*zbuoy(ig,l-1)/(fact_epsilon)
    609 
    610 !            zw2(ig,l+1)=Max(0.0001,(zdz/zdzbis)*(exp(-zw2fact)* &
    611 !    &                     (zw2(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
    612 !    &                     (exp(-zw2factbis)*(zw2(ig,l-1)-zdw2bis)+zdw2))
    613 !            zw2(ig,l+1)=Max(0.0001,(zw2(ig,l)+zdw2*zw2fact)*exp(-zw2fact))
    614             zw2(ig,l+1)=Max(0.0001,exp(-zw2fact)*(zw2(ig,l)-zdw2bis)+zdw2)
    615 
    616            endif
    617 
    618 
    619       endif
    620    enddo
    621 
    622         if (prt_level>=20) print*,'coucou calcul detr 460: ig, l',ig, l
    623 
    624 !===========================================================================
    625 ! 6. initialisations pour le calcul de la hauteur du thermique, de l'inversion et de la vitesse verticale max
    626 !===========================================================================
    627 
    628    nbpb=0
    629    do ig=1,ngrid
    630             if (zw2(ig,l+1)>0. .and. zw2(ig,l+1)<1.e-10) then
    631 !               stop'On tombe sur le cas particulier de thermcell_dry'
    632 !               print*,'On tombe sur le cas particulier de thermcell_plume'
    633                 nbpb=nbpb+1
    634                 zw2(ig,l+1)=0.
    635                 linter(ig)=l+1
    636             endif
    637 
    638         if (zw2(ig,l+1)<0.) then
    639            linter(ig)=(l*(zw2(ig,l+1)-zw2(ig,l))  &
    640                  -zw2(ig,l))/(zw2(ig,l+1)-zw2(ig,l))
    641            zw2(ig,l+1)=0.
    642 !+CR:04/05/12:correction calcul linter pour calcul de zmax continu
    643         elseif (f_star(ig,l+1)<0.) then
    644            linter(ig)=(l*(f_star(ig,l+1)-f_star(ig,l))  &
    645                  -f_star(ig,l))/(f_star(ig,l+1)-f_star(ig,l))
    646            zw2(ig,l+1)=0.
    647 !fin CR:04/05/12
    648         endif
    649 
    650            wa_moy(ig,l+1)=sqrt(zw2(ig,l+1))
    651 
    652         if (wa_moy(ig,l+1)>wmaxa(ig)) then
    653 !   lmix est le niveau de la couche ou w (wa_moy) est maximum
    654 !on rajoute le calcul de lmix_bis
    655             if (zqla(ig,l)<1.e-10) then
    656                lmix_bis(ig)=l+1
    657             endif
    658             lmix(ig)=l+1
    659             wmaxa(ig)=wa_moy(ig,l+1)
    660         endif
    661    enddo
    662 
    663    if (nbpb>0) then
    664    print*,'WARNING on tombe ',nbpb,' x sur un pb pour l=',l,' dans thermcell_plume'
    665    endif
    666 
    667 !=========================================================================
    668 ! FIN DE LA BOUCLE VERTICALE
    669       enddo
    670 !=========================================================================
    671 
    672 !on recalcule alim_star_tot
    673        do ig=1,ngrid
    674           alim_star_tot(ig)=0.
    675        enddo
    676        do ig=1,ngrid
    677           do l=1,lalim(ig)-1
    678           alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
    679           enddo
    680        enddo
    681        
    682 
    683         if (prt_level>=20) print*,'coucou calcul detr 470: ig, l', ig, l
    684 
    685 #undef wrgrads_thermcell
    686 #ifdef wrgrads_thermcell
    687          call wrgradsfi(1,nlay,entr_star(igout,1:nlay),'esta      ','esta      ')
    688          call wrgradsfi(1,nlay,detr_star(igout,1:nlay),'dsta      ','dsta      ')
    689          call wrgradsfi(1,nlay,zbuoy(igout,1:nlay),'buoy      ','buoy      ')
    690          call wrgradsfi(1,nlay,zdqt(igout,1:nlay),'dqt      ','dqt      ')
    691          call wrgradsfi(1,nlay,w_est(igout,1:nlay),'w_est     ','w_est     ')
    692          call wrgradsfi(1,nlay,w_est(igout,2:nlay+1),'w_es2     ','w_es2     ')
    693          call wrgradsfi(1,nlay,zw2(igout,1:nlay),'zw2A      ','zw2A      ')
    694 #endif
    695 
    696 
    697  RETURN
    698      end
    699 
    700 
    701 
    702 
    703 
    704 
    705 
    706 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    707 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    708 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    709 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    710 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    711 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    712  SUBROUTINE thermcell_plume_5B(itap,ngrid,nlay,ptimestep,ztv,zthl,po,zl,rhobarz,  &
    713 &           zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot,  &
    714 &           lalim,f0,detr_star,entr_star,f_star,csc,ztva,  &
    715 &           ztla,zqla,zqta,zha,zw2,w_est,ztva_est,zqsatth,lmix,lmix_bis,linter &
    716 &           ,lev_out,lunout1,igout)
    717 !&           ,lev_out,lunout1,igout,zbuoy,zbuoyjam)
    718 
    719 !--------------------------------------------------------------------------
    720 !thermcell_plume: calcule les valeurs de qt, thetal et w dans l ascendance
    721 ! Version conforme a l'article de Rio et al. 2010.
    722 ! Code ecrit par Catherine Rio, Arnaud Jam et Frederic Hourdin
    723 !--------------------------------------------------------------------------
    724 
    725       USE lmdz_thermcell_ini, ONLY: prt_level,fact_thermals_ed_dz,iflag_thermals_ed,RLvCP,RETV,RG
    726        USE lmdz_thermcell_qsat, ONLY: thermcell_qsat
    727       IMPLICIT NONE
    728 
    729       INTEGER itap
    730       INTEGER lunout1,igout
    731       INTEGER ngrid,nlay
    732       REAL ptimestep
    733       REAL ztv(ngrid,nlay)
    734       REAL zthl(ngrid,nlay)
    735       REAL, INTENT(IN) :: po(ngrid,nlay)
    736       REAL zl(ngrid,nlay)
    737       REAL rhobarz(ngrid,nlay)
    738       REAL zlev(ngrid,nlay+1)
    739       REAL pplev(ngrid,nlay+1)
    740       REAL pphi(ngrid,nlay)
    741       REAL zpspsk(ngrid,nlay)
    742       REAL alim_star(ngrid,nlay)
    743       REAL f0(ngrid)
    744       INTEGER lalim(ngrid)
    745       integer lev_out                           ! niveau pour les print
    746       integer nbpb
    747    
    748       real alim_star_tot(ngrid)
    749 
    750       REAL ztva(ngrid,nlay)
    751       REAL ztla(ngrid,nlay)
    752       REAL zqla(ngrid,nlay)
    753       REAL zqta(ngrid,nlay)
    754       REAL zha(ngrid,nlay)
    755 
    756       REAL detr_star(ngrid,nlay)
    757       REAL coefc
    758       REAL entr_star(ngrid,nlay)
    759       REAL detr(ngrid,nlay)
    760       REAL entr(ngrid,nlay)
    761 
    762       REAL csc(ngrid,nlay)
    763 
    764       REAL zw2(ngrid,nlay+1)
    765       REAL w_est(ngrid,nlay+1)
    766       REAL f_star(ngrid,nlay+1)
    767       REAL wa_moy(ngrid,nlay+1)
    768 
    769       REAL ztva_est(ngrid,nlay)
    770       REAL zqla_est(ngrid,nlay)
    771       REAL zqsatth(ngrid,nlay)
    772       REAL zta_est(ngrid,nlay)
    773       REAL zbuoyjam(ngrid,nlay)
    774       REAL ztemp(ngrid),zqsat(ngrid)
    775       REAL zdw2
    776       REAL zw2modif
    777       REAL zw2fact
    778       REAL zeps(ngrid,nlay)
    779 
    780       REAL linter(ngrid)
    781       INTEGER lmix(ngrid)
    782       INTEGER lmix_bis(ngrid)
    783       REAL    wmaxa(ngrid)
    784 
    785       INTEGER ig,l,k
    786 
    787       real zdz,zbuoy(ngrid,nlay),zalpha,gamma(ngrid,nlay),zdqt(ngrid,nlay),zw2m
    788       real zbuoybis
    789       real zcor,zdelta,zcvm5,qlbef,zdz2
    790       real betalpha,zbetalpha
    791       real eps, afact
    792       logical Zsat
    793       LOGICAL active(ngrid),activetmp(ngrid)
    794       REAL fact_gamma,fact_epsilon,fact_gamma2,fact_epsilon2
    795       REAL c2(ngrid,nlay)
    796       Zsat=.false.
    797 ! Initialisation
    798 
    799       fact_epsilon=0.002
    800       betalpha=0.9
    801       afact=2./3.           
    802 
    803       zbetalpha=betalpha/(1.+betalpha)
    804 
    805 
    806 ! Initialisations des variables reeles
    807 if (1==1) then
    808       ztva(:,:)=ztv(:,:)
    809       ztva_est(:,:)=ztva(:,:)
    810       ztla(:,:)=zthl(:,:)
    811       zqta(:,:)=po(:,:)
    812       zha(:,:) = ztva(:,:)
    813 else
    814       ztva(:,:)=0.
    815       ztva_est(:,:)=0.
    816       ztla(:,:)=0.
    817       zqta(:,:)=0.
    818       zha(:,:) =0.
    819 endif
    820 
    821       zqla_est(:,:)=0.
    822       zqsatth(:,:)=0.
    823       zqla(:,:)=0.
    824       detr_star(:,:)=0.
    825       entr_star(:,:)=0.
    826       alim_star(:,:)=0.
    827       alim_star_tot(:)=0.
    828       csc(:,:)=0.
    829       detr(:,:)=0.
    830       entr(:,:)=0.
    831       zw2(:,:)=0.
    832       zbuoy(:,:)=0.
    833       zbuoyjam(:,:)=0.
    834       gamma(:,:)=0.
    835       zeps(:,:)=0.
    836       w_est(:,:)=0.
    837       f_star(:,:)=0.
    838       wa_moy(:,:)=0.
    839       linter(:)=1.
    840 !     linter(:)=1.
    841 ! Initialisation des variables entieres
    842       lmix(:)=1
    843       lmix_bis(:)=2
    844       wmaxa(:)=0.
    845       lalim(:)=1
    846 
    847 
    848 !-------------------------------------------------------------------------
    849 ! On ne considere comme actif que les colonnes dont les deux premieres
    850 ! couches sont instables.
    851 !-------------------------------------------------------------------------
    852       active(:)=ztv(:,1)>ztv(:,2)
    853 
    854 !-------------------------------------------------------------------------
    855 ! Definition de l'alimentation
    856 !-------------------------------------------------------------------------
    857       do l=1,nlay-1
    858          do ig=1,ngrid
    859             if (ztv(ig,l)> ztv(ig,l+1) .and. ztv(ig,1)>=ztv(ig,l) ) then
    860                alim_star(ig,l)=MAX((ztv(ig,l)-ztv(ig,l+1)),0.)  &
    861                          *sqrt(zlev(ig,l+1))
    862                lalim(ig)=l+1
    863                alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
    864             endif
    865          enddo
    866       enddo
    867       do l=1,nlay
    868          do ig=1,ngrid
    869             if (alim_star_tot(ig) > 1.e-10 ) then
    870                alim_star(ig,l)=alim_star(ig,l)/alim_star_tot(ig)
    871             endif
    872          enddo
    873       enddo
    874       alim_star_tot(:)=1.
    875 
    876 
    877 
    878 !------------------------------------------------------------------------------
    879 ! Calcul dans la premiere couche
    880 ! On decide dans cette version que le thermique n'est actif que si la premiere
    881 ! couche est instable.
    882 ! Pourrait etre change si on veut que le thermiques puisse se d??clencher
    883 ! dans une couche l>1
    884 !------------------------------------------------------------------------------
    885 do ig=1,ngrid
    886 ! Le panache va prendre au debut les caracteristiques de l'air contenu
    887 ! dans cette couche.
    888     if (active(ig)) then
    889     ztla(ig,1)=zthl(ig,1)
    890     zqta(ig,1)=po(ig,1)
    891     zqla(ig,1)=zl(ig,1)
    892 !cr: attention, prise en compte de f*(1)=1
    893     f_star(ig,2)=alim_star(ig,1)
    894     zw2(ig,2)=2.*RG*(ztv(ig,1)-ztv(ig,2))/ztv(ig,2)  &
    895 &                     *(zlev(ig,2)-zlev(ig,1))  &
    896 &                     *0.4*pphi(ig,1)/(pphi(ig,2)-pphi(ig,1))
    897     w_est(ig,2)=zw2(ig,2)
    898     endif
    899 enddo
    900 
    901 !==============================================================================
    902 !boucle de calcul de la vitesse verticale dans le thermique
    903 !==============================================================================
    904 do l=2,nlay-1
    905 !==============================================================================
    906 
    907 
    908 ! On decide si le thermique est encore actif ou non
    909 ! AFaire : Il faut sans doute ajouter entr_star a alim_star dans ce test
    910     do ig=1,ngrid
    911        active(ig)=active(ig) &
    912 &                 .and. zw2(ig,l)>1.e-10 &
    913 &                 .and. f_star(ig,l)+alim_star(ig,l)>1.e-10
     1051    !=========================================================================
     1052
     1053    !on recalcule alim_star_tot
     1054    do ig = 1, ngrid
     1055      alim_star_tot(ig) = 0.
    9141056    enddo
    915 
    916 
    917 
    918 !---------------------------------------------------------------------------
    919 ! calcul des proprietes thermodynamiques et de la vitesse de la couche l
    920 ! sans tenir compte du detrainement et de l'entrainement dans cette
    921 ! couche
    922 ! C'est a dire qu'on suppose
    923 ! ztla(l)=ztla(l-1) et zqta(l)=zqta(l-1)
    924 ! Ici encore, on doit pouvoir ajouter entr_star (qui peut etre calculer
    925 ! avant) a l'alimentation pour avoir un calcul plus propre
    926 !---------------------------------------------------------------------------
    927 
    928    ztemp(:)=zpspsk(:,l)*ztla(:,l-1)
    929    call thermcell_qsat(ngrid,active,pplev(:,l),ztemp,zqta(:,l-1),zqsat(:))
    930 
    931     do ig=1,ngrid
    932 !       print*,'active',active(ig),ig,l
    933         if(active(ig)) then
    934         zqla_est(ig,l)=max(0.,zqta(ig,l-1)-zqsat(ig))
    935         ztva_est(ig,l) = ztla(ig,l-1)*zpspsk(ig,l)+RLvCp*zqla_est(ig,l)
    936         zta_est(ig,l)=ztva_est(ig,l)
    937         ztva_est(ig,l) = ztva_est(ig,l)/zpspsk(ig,l)
    938         ztva_est(ig,l) = ztva_est(ig,l)*(1.+RETV*(zqta(ig,l-1)  &
    939         -zqla_est(ig,l))-zqla_est(ig,l))
    940 
    941 !------------------------------------------------
    942 !AJAM:nouveau calcul de w? 
    943 !------------------------------------------------
    944               zdz=zlev(ig,l+1)-zlev(ig,l)
    945               zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
    946 
    947               zw2fact=fact_epsilon*2.*zdz/(1.+betalpha)
    948               zdw2=(afact)*zbuoy(ig,l)/(fact_epsilon)
    949               w_est(ig,l+1)=Max(0.0001,exp(-zw2fact)*(w_est(ig,l)-zdw2)+zdw2)
    950  
    951 
    952              if (w_est(ig,l+1)<0.) then
    953                 w_est(ig,l+1)=zw2(ig,l)
    954              endif
    955        endif
     1057    do ig = 1, ngrid
     1058      do l = 1, lalim(ig) - 1
     1059        alim_star_tot(ig) = alim_star_tot(ig) + alim_star(ig, l)
     1060      enddo
    9561061    enddo
    9571062
    958 
    959 !-------------------------------------------------
    960 !calcul des taux d'entrainement et de detrainement
    961 !-------------------------------------------------
    962 
    963      do ig=1,ngrid
    964         if (active(ig)) then
    965 
    966           zw2m=max(0.5*(w_est(ig,l)+w_est(ig,l+1)),0.1)
    967           zw2m=w_est(ig,l+1)
    968           zdz=zlev(ig,l+1)-zlev(ig,l)
    969           zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
    970 !          zbuoybis=zbuoy(ig,l)+RG*0.1/300.
    971           zbuoybis=zbuoy(ig,l)
    972           zalpha=f0(ig)*f_star(ig,l)/sqrt(w_est(ig,l+1))/rhobarz(ig,l)
    973           zdqt(ig,l)=max(zqta(ig,l-1)-po(ig,l),0.)/po(ig,l)
    974 
    975          
    976           entr_star(ig,l)=f_star(ig,l)*zdz*  zbetalpha*MAX(0.,  &
    977        afact*zbuoybis/zw2m - fact_epsilon )
    978 
    979 
    980           detr_star(ig,l)=f_star(ig,l)*zdz                        &
    981        *MAX(1.e-3, -afact*zbetalpha*zbuoy(ig,l)/zw2m          &
    982        + 0.012*(zdqt(ig,l)/zw2m)**0.5 )
    983          
    984 ! En dessous de lalim, on prend le max de alim_star et entr_star pour
    985 ! alim_star et 0 sinon
    986         if (l<lalim(ig)) then
    987           alim_star(ig,l)=max(alim_star(ig,l),entr_star(ig,l))
    988           entr_star(ig,l)=0.
    989         endif
    990 
    991 ! Calcul du flux montant normalise
    992       f_star(ig,l+1)=f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l)  &
    993                 -detr_star(ig,l)
    994 
    995       endif
    996    enddo
    997 
    998 
    999 !----------------------------------------------------------------------------
    1000 !calcul de la vitesse verticale en melangeant Tl et qt du thermique
    1001 !---------------------------------------------------------------------------
    1002    activetmp(:)=active(:) .and. f_star(:,l+1)>1.e-10
    1003    do ig=1,ngrid
    1004        if (activetmp(ig)) then
    1005            Zsat=.false.
    1006            ztla(ig,l)=(f_star(ig,l)*ztla(ig,l-1)+  &
    1007               (alim_star(ig,l)+entr_star(ig,l))*zthl(ig,l))  &
    1008               /(f_star(ig,l+1)+detr_star(ig,l))
    1009            zqta(ig,l)=(f_star(ig,l)*zqta(ig,l-1)+  &
    1010               (alim_star(ig,l)+entr_star(ig,l))*po(ig,l))  &
    1011               /(f_star(ig,l+1)+detr_star(ig,l))
    1012 
    1013         endif
    1014     enddo
    1015 
    1016    ztemp(:)=zpspsk(:,l)*ztla(:,l)
    1017    call thermcell_qsat(ngrid,activetmp,pplev(:,l),ztemp,zqta(:,l),zqsatth(:,l))
    1018 
    1019    do ig=1,ngrid
    1020       if (activetmp(ig)) then
    1021 ! on ecrit de maniere conservative (sat ou non)
    1022 !          T = Tl +Lv/Cp ql
    1023            zqla(ig,l)=max(0.,zqta(ig,l)-zqsatth(ig,l))
    1024            ztva(ig,l) = ztla(ig,l)*zpspsk(ig,l)+RLvCp*zqla(ig,l)
    1025            ztva(ig,l) = ztva(ig,l)/zpspsk(ig,l)
    1026 !on rajoute le calcul de zha pour diagnostiques (temp potentielle)
    1027            zha(ig,l) = ztva(ig,l)
    1028            ztva(ig,l) = ztva(ig,l)*(1.+RETV*(zqta(ig,l)  &
    1029                 -zqla(ig,l))-zqla(ig,l))
    1030            zbuoy(ig,l)=RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)
    1031            zdz=zlev(ig,l+1)-zlev(ig,l)
    1032            zeps(ig,l)=(entr_star(ig,l)+alim_star(ig,l))/(f_star(ig,l)*zdz)
    1033 
    1034             zw2fact=fact_epsilon*2.*zdz/(1.+betalpha)
    1035             zdw2=afact*zbuoy(ig,l)/(fact_epsilon)
    1036             zw2(ig,l+1)=Max(0.0001,exp(-zw2fact)*(zw2(ig,l)-zdw2)+zdw2)
    1037       endif
    1038    enddo
    1039 
    1040         if (prt_level>=20) print*,'coucou calcul detr 460: ig, l',ig, l
    1041 
    1042 !---------------------------------------------------------------------------
    1043 !initialisations pour le calcul de la hauteur du thermique, de l'inversion et de la vitesse verticale max
    1044 !---------------------------------------------------------------------------
    1045 
    1046    nbpb=0
    1047    do ig=1,ngrid
    1048             if (zw2(ig,l+1)>0. .and. zw2(ig,l+1)<1.e-10) then
    1049 !               stop'On tombe sur le cas particulier de thermcell_dry'
    1050 !               print*,'On tombe sur le cas particulier de thermcell_plume'
    1051                 nbpb=nbpb+1
    1052                 zw2(ig,l+1)=0.
    1053                 linter(ig)=l+1
    1054             endif
    1055 
    1056         if (zw2(ig,l+1)<0.) then
    1057            linter(ig)=(l*(zw2(ig,l+1)-zw2(ig,l))  &
    1058                  -zw2(ig,l))/(zw2(ig,l+1)-zw2(ig,l))
    1059            zw2(ig,l+1)=0.
    1060         elseif (f_star(ig,l+1)<0.) then
    1061            linter(ig)=(l*(f_star(ig,l+1)-f_star(ig,l))  &
    1062                  -f_star(ig,l))/(f_star(ig,l+1)-f_star(ig,l))
    1063 !           print*,"linter plume", linter(ig)
    1064            zw2(ig,l+1)=0.
    1065         endif
    1066 
    1067            wa_moy(ig,l+1)=sqrt(zw2(ig,l+1))
    1068 
    1069         if (wa_moy(ig,l+1)>wmaxa(ig)) then
    1070 !   lmix est le niveau de la couche ou w (wa_moy) est maximum
    1071 !on rajoute le calcul de lmix_bis
    1072             if (zqla(ig,l)<1.e-10) then
    1073                lmix_bis(ig)=l+1
    1074             endif
    1075             lmix(ig)=l+1
    1076             wmaxa(ig)=wa_moy(ig,l+1)
    1077         endif
    1078    enddo
    1079 
    1080    if (nbpb>0) then
    1081    print*,'WARNING on tombe ',nbpb,' x sur un pb pour l=',l,' dans thermcell_plume'
    1082    endif
    1083 
    1084 !=========================================================================
    1085 ! FIN DE LA BOUCLE VERTICALE
    1086       enddo
    1087 !=========================================================================
    1088 
    1089 !on recalcule alim_star_tot
    1090        do ig=1,ngrid
    1091           alim_star_tot(ig)=0.
    1092        enddo
    1093        do ig=1,ngrid
    1094           do l=1,lalim(ig)-1
    1095           alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
    1096           enddo
    1097        enddo
    1098        
    1099 
    1100         if (prt_level>=20) print*,'coucou calcul detr 470: ig, l', ig, l
    1101 
    1102 #undef wrgrads_thermcell
    1103 #ifdef wrgrads_thermcell
    1104          call wrgradsfi(1,nlay,entr_star(igout,1:nlay),'esta      ','esta      ')
    1105          call wrgradsfi(1,nlay,detr_star(igout,1:nlay),'dsta      ','dsta      ')
    1106          call wrgradsfi(1,nlay,zbuoy(igout,1:nlay),'buoy      ','buoy      ')
    1107          call wrgradsfi(1,nlay,zdqt(igout,1:nlay),'dqt      ','dqt      ')
    1108          call wrgradsfi(1,nlay,w_est(igout,1:nlay),'w_est     ','w_est     ')
    1109          call wrgradsfi(1,nlay,w_est(igout,2:nlay+1),'w_es2     ','w_es2     ')
    1110          call wrgradsfi(1,nlay,zw2(igout,1:nlay),'zw2A      ','zw2A      ')
    1111 #endif
    1112 
    1113 
    1114      return
    1115      end
     1063    if (prt_level>=20) print*, 'coucou calcul detr 470: ig, l', ig, l
     1064  end
    11161065END MODULE lmdz_thermcell_plume_6A

  • LMDZ6/branches/Amaury_dev/libf/phylmd/surf_landice_mod.F90

    r5101 r5102  
    1 
    21MODULE surf_landice_mod
    3  
     2
    43  IMPLICIT NONE
    54
    65CONTAINS
    76
    8 !****************************************************************************************
     7  !****************************************************************************************
    98
    109  SUBROUTINE surf_landice(itime, dtime, knon, knindex, &
    11        rlon, rlat, debut, lafin, &
    12        rmu0, lwdownm, albedo, pphi1, &
    13        swnet, lwnet, tsurf, p1lay, &
    14        cdragh, cdragm, precip_rain, precip_snow, precip_bs, temp_air, spechum, &
    15        AcoefH, AcoefQ, BcoefH, BcoefQ, &
    16        AcoefU, AcoefV, BcoefU, BcoefV, &
    17        AcoefQBS, BcoefQBS, &
    18        ps, u1, v1, gustiness, rugoro, pctsrf, &
    19        snow, qsurf, qsol, qbs1, agesno, &
    20        tsoil, z0m, z0h, SFRWL, alb_dir, alb_dif, evap, fluxsens, fluxlat, fluxbs, &
    21        tsurf_new, dflux_s, dflux_l, &
    22        alt, slope, cloudf, &
    23        snowhgt, qsnow, to_ice, sissnow, &
    24        alb3, runoff, &
    25        flux_u1, flux_v1 &
     10            rlon, rlat, debut, lafin, &
     11            rmu0, lwdownm, albedo, pphi1, &
     12            swnet, lwnet, tsurf, p1lay, &
     13            cdragh, cdragm, precip_rain, precip_snow, precip_bs, temp_air, spechum, &
     14            AcoefH, AcoefQ, BcoefH, BcoefQ, &
     15            AcoefU, AcoefV, BcoefU, BcoefV, &
     16            AcoefQBS, BcoefQBS, &
     17            ps, u1, v1, gustiness, rugoro, pctsrf, &
     18            snow, qsurf, qsol, qbs1, agesno, &
     19            tsoil, z0m, z0h, SFRWL, alb_dir, alb_dif, evap, fluxsens, fluxlat, fluxbs, &
     20            tsurf_new, dflux_s, dflux_l, &
     21            alt, slope, cloudf, &
     22            snowhgt, qsnow, to_ice, sissnow, &
     23            alb3, runoff, &
     24            flux_u1, flux_v1 &
    2625#ifdef ISO
    2726        ,xtprecip_rain, xtprecip_snow,xtspechum,Rland_ice &
    2827        ,xtsnow,xtsol,xtevap &
    2928#endif               
    30       )
     29            )
    3130
    3231    USE dimphy
    33     USE geometry_mod,     ONLY: longitude,latitude
    34     USE surface_data,     ONLY: type_ocean, calice, calsno, landice_opt, iflag_albcalc
    35     USE fonte_neige_mod,  ONLY: fonte_neige,run_off_lic,fqcalving_global,ffonte_global,fqfonte_global,runofflic_global
    36     USE cpl_mod,          ONLY: cpl_send_landice_fields
     32    USE geometry_mod, ONLY: longitude, latitude
     33    USE surface_data, ONLY: type_ocean, calice, calsno, landice_opt, iflag_albcalc
     34    USE fonte_neige_mod, ONLY: fonte_neige, run_off_lic, fqcalving_global, ffonte_global, fqfonte_global, runofflic_global
     35    USE cpl_mod, ONLY: cpl_send_landice_fields
    3736    USE calcul_fluxs_mod
    3837    USE phys_local_var_mod, ONLY: zxrhoslic, zxustartlic, zxqsaltlic
    39     USE phys_output_var_mod, ONLY: snow_o,zfra_o
    40 #ifdef ISO   
     38    USE phys_output_var_mod, ONLY: snow_o, zfra_o
     39#ifdef ISO
    4140    USE fonte_neige_mod,  ONLY: xtrun_off_lic
    4241    USE infotrac_phy,     ONLY: ntiso,niso
    … …  
    4746#endif
    4847#endif
    49  
    50 !FC
     48
     49    !FC
    5150    USE ioipsl_getin_p_mod, ONLY: getin_p
    5251    USE lmdz_blowing_snow_ini, ONLY: c_esalt_bs, zeta_bs, pbst_bs, prt_bs, rhoice_bs, rhohard_bs
    5352    USE lmdz_blowing_snow_ini, ONLY: rhofresh_bs, tau_eqsalt_bs, tau_dens0_bs, tau_densmin_bs
    54 #ifdef CPP_INLANDSIS
    55     USE surf_inlandsis_mod,  ONLY: surf_inlandsis
    56 #endif
     53    USE surf_inlandsis_mod, ONLY: surf_inlandsis
     54    USE lmdz_cppkeys_wrapper, ONLY: CPPKEY_INLANDSIS
    5755
    5856    USE indice_sol_mod
    5957
    60 !    INCLUDE "indicesol.h"
     58    !    INCLUDE "indicesol.h"
    6159    INCLUDE "dimsoil.h"
    6260    INCLUDE "YOMCST.h"
    6361    INCLUDE "clesphys.h"
    6462
    65 ! Input variables
    66 !****************************************************************************************
    67     INTEGER, INTENT(IN)                           :: itime, knon
    68     INTEGER, DIMENSION(klon), INTENT(in)          :: knindex
    69     REAL, INTENT(in)                              :: dtime
    70     REAL, DIMENSION(klon), INTENT(IN)             :: swnet ! net shortwave radiance
    71     REAL, DIMENSION(klon), INTENT(IN)             :: lwnet ! net longwave radiance
    72     REAL, DIMENSION(klon), INTENT(IN)             :: tsurf
    73     REAL, DIMENSION(klon), INTENT(IN)             :: p1lay
    74     REAL, DIMENSION(klon), INTENT(IN)             :: cdragh, cdragm
    75     REAL, DIMENSION(klon), INTENT(IN)             :: precip_rain, precip_snow, precip_bs
    76     REAL, DIMENSION(klon), INTENT(IN)             :: temp_air, spechum
    77     REAL, DIMENSION(klon), INTENT(IN)             :: AcoefH, AcoefQ
    78     REAL, DIMENSION(klon), INTENT(IN)             :: BcoefH, BcoefQ
    79     REAL, DIMENSION(klon), INTENT(IN)             :: AcoefU, AcoefV, BcoefU, BcoefV
    80     REAL, DIMENSION(klon), INTENT(IN)             :: AcoefQBS, BcoefQBS
    81     REAL, DIMENSION(klon), INTENT(IN)             :: ps
    82     REAL, DIMENSION(klon), INTENT(IN)             :: u1, v1, gustiness, qbs1
    83     REAL, DIMENSION(klon), INTENT(IN)             :: rugoro
    84     REAL, DIMENSION(klon,nbsrf), INTENT(IN)      :: pctsrf
     63    ! Input variables
     64    !****************************************************************************************
     65    INTEGER, INTENT(IN) :: itime, knon
     66    INTEGER, DIMENSION(klon), INTENT(in) :: knindex
     67    REAL, INTENT(in) :: dtime
     68    REAL, DIMENSION(klon), INTENT(IN) :: swnet ! net shortwave radiance
     69    REAL, DIMENSION(klon), INTENT(IN) :: lwnet ! net longwave radiance
     70    REAL, DIMENSION(klon), INTENT(IN) :: tsurf
     71    REAL, DIMENSION(klon), INTENT(IN) :: p1lay
     72    REAL, DIMENSION(klon), INTENT(IN) :: cdragh, cdragm
     73    REAL, DIMENSION(klon), INTENT(IN) :: precip_rain, precip_snow, precip_bs
     74    REAL, DIMENSION(klon), INTENT(IN) :: temp_air, spechum
     75    REAL, DIMENSION(klon), INTENT(IN) :: AcoefH, AcoefQ
     76    REAL, DIMENSION(klon), INTENT(IN) :: BcoefH, BcoefQ
     77    REAL, DIMENSION(klon), INTENT(IN) :: AcoefU, AcoefV, BcoefU, BcoefV
     78    REAL, DIMENSION(klon), INTENT(IN) :: AcoefQBS, BcoefQBS
     79    REAL, DIMENSION(klon), INTENT(IN) :: ps
     80    REAL, DIMENSION(klon), INTENT(IN) :: u1, v1, gustiness, qbs1
     81    REAL, DIMENSION(klon), INTENT(IN) :: rugoro
     82    REAL, DIMENSION(klon, nbsrf), INTENT(IN) :: pctsrf
    8583#ifdef ISO
    8684    REAL, DIMENSION(ntiso,klon), INTENT(IN)       :: xtprecip_rain, xtprecip_snow
    … …  
    8987
    9088
    91     LOGICAL,  INTENT(IN)                          :: debut   !true if first step
    92     LOGICAL,  INTENT(IN)                          :: lafin   !true if last step
    93     REAL, DIMENSION(klon), INTENT(IN)             :: rlon, rlat
    94     REAL, DIMENSION(klon), INTENT(IN)             :: rmu0
    95     REAL, DIMENSION(klon), INTENT(IN)             :: lwdownm !ylwdown
    96     REAL, DIMENSION(klon), INTENT(IN)             :: albedo  !mean albedo
    97     REAL, DIMENSION(klon), INTENT(IN)             :: pphi1   
    98     REAL, DIMENSION(klon), INTENT(IN)             :: alt   !mean altitude of the grid box 
    99     REAL, DIMENSION(klon), INTENT(IN)             :: slope   !mean slope in grid box 
    100     REAL, DIMENSION(klon), INTENT(IN)             :: cloudf  !total cloud fraction
    101 
    102 ! In/Output variables
    103 !****************************************************************************************
    104     REAL, DIMENSION(klon), INTENT(INOUT)          :: snow, qsol
    105     REAL, DIMENSION(klon), INTENT(INOUT)          :: agesno
     89    LOGICAL, INTENT(IN) :: debut   !true if first step
     90    LOGICAL, INTENT(IN) :: lafin   !true if last step
     91    REAL, DIMENSION(klon), INTENT(IN) :: rlon, rlat
     92    REAL, DIMENSION(klon), INTENT(IN) :: rmu0
     93    REAL, DIMENSION(klon), INTENT(IN) :: lwdownm !ylwdown
     94    REAL, DIMENSION(klon), INTENT(IN) :: albedo  !mean albedo
     95    REAL, DIMENSION(klon), INTENT(IN) :: pphi1
     96    REAL, DIMENSION(klon), INTENT(IN) :: alt   !mean altitude of the grid box
     97    REAL, DIMENSION(klon), INTENT(IN) :: slope   !mean slope in grid box
     98    REAL, DIMENSION(klon), INTENT(IN) :: cloudf  !total cloud fraction
     99
     100    ! In/Output variables
     101    !****************************************************************************************
     102    REAL, DIMENSION(klon), INTENT(INOUT) :: snow, qsol
     103    REAL, DIMENSION(klon), INTENT(INOUT) :: agesno
    106104    REAL, DIMENSION(klon, nsoilmx), INTENT(INOUT) :: tsoil
    107105#ifdef ISO
    … …  
    111109
    112110
    113 ! Output variables
    114 !****************************************************************************************
    115     REAL, DIMENSION(klon), INTENT(OUT)            :: qsurf
    116     REAL, DIMENSION(klon), INTENT(OUT)            :: z0m, z0h
    117 !albedo SB >>>
    118 !    REAL, DIMENSION(klon), INTENT(OUT)            :: alb1  ! new albedo in visible SW interval
    119 !    REAL, DIMENSION(klon), INTENT(OUT)            :: alb2  ! new albedo in near IR interval
    120     REAL, DIMENSION(6), INTENT(IN)                :: SFRWL
    121     REAL, DIMENSION(klon,nsw), INTENT(OUT)        :: alb_dir,alb_dif
    122 !albedo SB <<<
    123     REAL, DIMENSION(klon), INTENT(OUT)            :: evap, fluxsens, fluxlat
    124     REAL, DIMENSION(klon), INTENT(OUT)            :: fluxbs
    125     REAL, DIMENSION(klon), INTENT(OUT)            :: tsurf_new
    126     REAL, DIMENSION(klon), INTENT(OUT)            :: dflux_s, dflux_l     
    127     REAL, DIMENSION(klon), INTENT(OUT)            :: flux_u1, flux_v1
    128 
    129     REAL, DIMENSION(klon), INTENT(OUT)           :: alb3
    130     REAL, DIMENSION(klon), INTENT(OUT)           :: qsnow   !column water in snow [kg/m2]
    131     REAL, DIMENSION(klon), INTENT(OUT)           :: snowhgt !Snow height (m)
    132     REAL, DIMENSION(klon), INTENT(OUT)           :: to_ice
    133     REAL, DIMENSION(klon), INTENT(OUT)           :: sissnow
    134     REAL, DIMENSION(klon), INTENT(OUT)           :: runoff  !Land ice runoff
     111    ! Output variables
     112    !****************************************************************************************
     113    REAL, DIMENSION(klon), INTENT(OUT) :: qsurf
     114    REAL, DIMENSION(klon), INTENT(OUT) :: z0m, z0h
     115    !albedo SB >>>
     116    !    REAL, DIMENSION(klon), INTENT(OUT)            :: alb1  ! new albedo in visible SW interval
     117    !    REAL, DIMENSION(klon), INTENT(OUT)            :: alb2  ! new albedo in near IR interval
     118    REAL, DIMENSION(6), INTENT(IN) :: SFRWL
     119    REAL, DIMENSION(klon, nsw), INTENT(OUT) :: alb_dir, alb_dif
     120    !albedo SB <<<
     121    REAL, DIMENSION(klon), INTENT(OUT) :: evap, fluxsens, fluxlat
     122    REAL, DIMENSION(klon), INTENT(OUT) :: fluxbs
     123    REAL, DIMENSION(klon), INTENT(OUT) :: tsurf_new
     124    REAL, DIMENSION(klon), INTENT(OUT) :: dflux_s, dflux_l
     125    REAL, DIMENSION(klon), INTENT(OUT) :: flux_u1, flux_v1
     126
     127    REAL, DIMENSION(klon), INTENT(OUT) :: alb3
     128    REAL, DIMENSION(klon), INTENT(OUT) :: qsnow   !column water in snow [kg/m2]
     129    REAL, DIMENSION(klon), INTENT(OUT) :: snowhgt !Snow height (m)
     130    REAL, DIMENSION(klon), INTENT(OUT) :: to_ice
     131    REAL, DIMENSION(klon), INTENT(OUT) :: sissnow
     132    REAL, DIMENSION(klon), INTENT(OUT) :: runoff  !Land ice runoff
    135133#ifdef ISO
    136134    REAL, DIMENSION(ntiso,klon), INTENT(OUT)     :: xtevap     
    … …  
    138136!    fonte_neige
    139137#endif
    140  
    141 
    142 ! Local variables
    143 !****************************************************************************************
    144     REAL, DIMENSION(klon)    :: soilcap, soilflux
    145     REAL, DIMENSION(klon)    :: cal, beta, dif_grnd
    146     REAL, DIMENSION(klon)    :: zfra, alb_neig
    147     REAL, DIMENSION(klon)    :: radsol
    148     REAL, DIMENSION(klon)    :: u0, v0, u1_lay, v1_lay, ustar
    149     INTEGER                  :: i,j,nt
    150     REAL, DIMENSION(klon)    :: fqfonte,ffonte
    151     REAL, DIMENSION(klon)    :: run_off_lic_frac
     138
     139
     140    ! Local variables
     141    !****************************************************************************************
     142    REAL, DIMENSION(klon) :: soilcap, soilflux
     143    REAL, DIMENSION(klon) :: cal, beta, dif_grnd
     144    REAL, DIMENSION(klon) :: zfra, alb_neig
     145    REAL, DIMENSION(klon) :: radsol
     146    REAL, DIMENSION(klon) :: u0, v0, u1_lay, v1_lay, ustar
     147    INTEGER :: i, j, nt
     148    REAL, DIMENSION(klon) :: fqfonte, ffonte
     149    REAL, DIMENSION(klon) :: run_off_lic_frac
    152150#ifdef ISO       
    153151    REAL, PARAMETER          :: t_coup = 273.15
    … …  
    167165
    168166
    169     REAL, DIMENSION(klon)    :: emis_new                  !Emissivity
    170     REAL, DIMENSION(klon)    :: swdown,lwdown
    171     REAL, DIMENSION(klon)    :: precip_snow_adv, snow_adv !Snow Drift precip./advection (not used in inlandsis)
    172     REAL, DIMENSION(klon)    :: erod                      !erosion of surface snow (flux, kg/m2/s like evap)
    173     REAL, DIMENSION(klon)    :: zsl_height, wind_velo     !surface layer height, wind spd
    174     REAL, DIMENSION(klon)    :: dens_air,  snow_cont_air  !air density; snow content air
    175     REAL, DIMENSION(klon)    :: alb_soil                  !albedo of underlying ice
    176     REAL, DIMENSION(klon)    :: pexner                    !Exner potential
    177     REAL                     :: pref
    178     REAL, DIMENSION(klon,nsoilmx) :: tsoil0               !modif
    179     REAL                          :: dtis                ! subtimestep
    180     LOGICAL                       :: debut_is, lafin_is  ! debut and lafin for inlandsis
    181 
    182     CHARACTER (len = 20)                      :: modname = 'surf_landice'
    183     CHARACTER (len = 80)                      :: abort_message
    184 
    185 
    186     REAL,DIMENSION(klon) :: alb1,alb2
    187     REAL,DIMENSION(klon) :: precip_totsnow, evap_totsnow
    188     REAL, DIMENSION (klon,6) :: alb6
    189     REAL                   :: esalt
    190     REAL                   :: lambdasalt,fluxsalt, csalt, nunu, aa, bb, cc
    191     REAL                   :: tau_dens, maxerosion
    192     REAL, DIMENSION(klon)  :: ws1, rhod, rhos, ustart0, ustart, qsalt, hsalt
    193     REAL, DIMENSION(klon)  :: fluxbs_1, fluxbs_2, bsweight_fresh
     167    REAL, DIMENSION(klon) :: emis_new                  !Emissivity
     168    REAL, DIMENSION(klon) :: swdown, lwdown
     169    REAL, DIMENSION(klon) :: precip_snow_adv, snow_adv !Snow Drift precip./advection (not used in inlandsis)
     170    REAL, DIMENSION(klon) :: erod                      !erosion of surface snow (flux, kg/m2/s like evap)
     171    REAL, DIMENSION(klon) :: zsl_height, wind_velo     !surface layer height, wind spd
     172    REAL, DIMENSION(klon) :: dens_air, snow_cont_air  !air density; snow content air
     173    REAL, DIMENSION(klon) :: alb_soil                  !albedo of underlying ice
     174    REAL, DIMENSION(klon) :: pexner                    !Exner potential
     175    REAL :: pref
     176    REAL, DIMENSION(klon, nsoilmx) :: tsoil0               !modif
     177    REAL :: dtis                ! subtimestep
     178    LOGICAL :: debut_is, lafin_is  ! debut and lafin for inlandsis
     179
     180    CHARACTER (len = 20) :: modname = 'surf_landice'
     181    CHARACTER (len = 80) :: abort_message
     182
     183    REAL, DIMENSION(klon) :: alb1, alb2
     184    REAL, DIMENSION(klon) :: precip_totsnow, evap_totsnow
     185    REAL, DIMENSION (klon, 6) :: alb6
     186    REAL :: esalt
     187    REAL :: lambdasalt, fluxsalt, csalt, nunu, aa, bb, cc
     188    REAL :: tau_dens, maxerosion
     189    REAL, DIMENSION(klon) :: ws1, rhod, rhos, ustart0, ustart, qsalt, hsalt
     190    REAL, DIMENSION(klon) :: fluxbs_1, fluxbs_2, bsweight_fresh
    194191    LOGICAL, DIMENSION(klon) :: ok_remaining_freshsnow
    195     REAL  :: ta1, ta2, ta3, z01, z02, z03, coefa, coefb, coefc, coefd
    196 
    197 
    198 ! End definition
    199 !****************************************************************************************
    200 !FC
    201 !FC
    202    REAL,SAVE :: alb_vis_sno_lic
    203   !$OMP THREADPRIVATE(alb_vis_sno_lic)
    204    REAL,SAVE :: alb_nir_sno_lic
    205   !$OMP THREADPRIVATE(alb_nir_sno_lic)
    206   LOGICAL, SAVE :: firstcall = .TRUE.
    207   !$OMP THREADPRIVATE(firstcall)
    208 
    209 
    210 !FC firtscall initializations
    211 !******************************************************************************************
     192    REAL :: ta1, ta2, ta3, z01, z02, z03, coefa, coefb, coefc, coefd
     193
     194
     195    ! End definition
     196    !****************************************************************************************
     197    !FC
     198    !FC
     199    REAL, SAVE :: alb_vis_sno_lic
     200    !$OMP THREADPRIVATE(alb_vis_sno_lic)
     201    REAL, SAVE :: alb_nir_sno_lic
     202    !$OMP THREADPRIVATE(alb_nir_sno_lic)
     203    LOGICAL, SAVE :: firstcall = .TRUE.
     204    !$OMP THREADPRIVATE(firstcall)
     205
     206
     207    !FC firtscall initializations
     208    !******************************************************************************************
    212209#ifdef ISO
    213210#ifdef ISOVERIF
    … …  
    222219#endif
    223220
    224   IF (firstcall) THEN
    225   alb_vis_sno_lic=0.77
    226   CALL getin_p('alb_vis_sno_lic',alb_vis_sno_lic)
    227            PRINT*, 'alb_vis_sno_lic',alb_vis_sno_lic
    228   alb_nir_sno_lic=0.77
    229   CALL getin_p('alb_nir_sno_lic',alb_nir_sno_lic)
    230            PRINT*, 'alb_nir_sno_lic',alb_nir_sno_lic
    231  
    232   firstcall=.false.
    233   ENDIF
    234 !******************************************************************************************
    235 
    236 ! Initialize output variables
     221    IF (firstcall) THEN
     222      alb_vis_sno_lic = 0.77
     223      CALL getin_p('alb_vis_sno_lic', alb_vis_sno_lic)
     224      PRINT*, 'alb_vis_sno_lic', alb_vis_sno_lic
     225      alb_nir_sno_lic = 0.77
     226      CALL getin_p('alb_nir_sno_lic', alb_nir_sno_lic)
     227      PRINT*, 'alb_nir_sno_lic', alb_nir_sno_lic
     228
     229      firstcall = .false.
     230    ENDIF
     231    !******************************************************************************************
     232
     233    ! Initialize output variables
    237234    alb3(:) = 999999.
    238235    alb2(:) = 999999.
    239236    alb1(:) = 999999.
    240     fluxbs(:)=0. 
     237    fluxbs(:) = 0.
    241238    runoff(:) = 0.
    242 !****************************************************************************************
    243 ! Calculate total absorbed radiance at surface
    244 
    245 !****************************************************************************************
     239    !****************************************************************************************
     240    ! Calculate total absorbed radiance at surface
     241
     242    !****************************************************************************************
    246243    radsol(:) = 0.0
    247244    radsol(1:knon) = swnet(1:knon) + lwnet(1:knon)
    248245
    249 !****************************************************************************************
    250 
    251 !****************************************************************************************
    252 !  landice_opt = 0 : soil_model, calcul_flux, fonte_neige, ... 
    253 !  landice_opt = 1  : prepare and call INterace Lmdz SISvat (INLANDSIS)
    254 !****************************************************************************************
    255 
     246    !****************************************************************************************
     247
     248    !****************************************************************************************
     249    !  landice_opt = 0 : soil_model, calcul_flux, fonte_neige, ...
     250    !  landice_opt = 1  : prepare and call INterace Lmdz SISvat (INLANDSIS)
     251    !****************************************************************************************
    256252
    257253    IF (landice_opt == 1) THEN
    258254
    259 !****************************************************************************************   
    260 ! CALL to INLANDSIS interface
    261 !****************************************************************************************
    262 #ifdef CPP_INLANDSIS
     255      !****************************************************************************************
     256      ! CALL to INLANDSIS interface
     257      !****************************************************************************************
     258      IF (CPPKEY_INLANDSIS) THEN
    263259
    264260#ifdef ISO
    … …  
    266262#endif
    267263
    268         debut_is=debut
    269         lafin_is=.false.
     264        debut_is = debut
     265        lafin_is = .false.
    270266        ! Suppose zero surface speed
    271         u0(:)            = 0.0
    272         v0(:)            = 0.0
    273 
     267        u0(:) = 0.0
     268        v0(:) = 0.0
    274269
    275270        CALL calcul_flux_wind(knon, dtime, &
    276          u0, v0, u1, v1, gustiness, cdragm, &
    277          AcoefU, AcoefV, BcoefU, BcoefV, &
    278          p1lay, temp_air, &
    279          flux_u1, flux_v1)
    280 
    281        
    282        ! Set constants and compute some input for SISVAT
    283        ! = 1000 hPa
    284        ! and calculate incoming flux for SW and LW interval: swdown, lwdown
    285        swdown(:)        = 0.0
    286        lwdown(:)        = 0.0
    287        snow_cont_air(:) = 0.  ! the snow content in air is not a prognostic variable of the model     
    288        alb_soil(:)      = 0.4 ! before albedo(:) but here it is the ice albedo that we have to set
    289        ustar(:)         = 0.
    290        pref             = 100000.       
    291        DO i = 1, knon
    292           swdown(i)        = swnet(i)/(1-albedo(i))
    293           lwdown(i)        = lwdownm(i)
    294           wind_velo(i)     = u1(i)**2 + v1(i)**2
    295           wind_velo(i)     = wind_velo(i)**0.5
    296           pexner(i)        = (p1lay(i)/pref)**(RD/RCPD)
    297           dens_air(i)      = p1lay(i)/RD/temp_air(i)  ! dry air density
    298           zsl_height(i)    = pphi1(i)/RG     
    299           tsoil0(i,:)      = tsoil(i,:) 
    300           ustar(i)= (cdragm(i)*(wind_velo(i)**2))**0.5   
    301        END DO
    302        
    303 
    304 
    305         dtis=dtime
    306 
    307           IF (lafin) THEN
    308             lafin_is=.true.
    309           END IF
    310 
    311           CALL surf_inlandsis(knon, rlon, rlat, knindex, itime, dtis, debut_is, lafin_is,&
    312             rmu0, swdown, lwdown, albedo, pexner, ps, p1lay, precip_rain, precip_snow,   &
    313             zsl_height, wind_velo, ustar, temp_air, dens_air, spechum, tsurf,&
    314             rugoro, snow_cont_air, alb_soil, alt, slope, cloudf, &
    315             radsol, qsol, tsoil0, snow, zfra, snowhgt, qsnow, to_ice, sissnow,agesno,   &
    316             AcoefH, AcoefQ, BcoefH, BcoefQ, cdragm, cdragh, &
    317             run_off_lic, fqfonte, ffonte, evap, erod, fluxsens, fluxlat,dflux_s, dflux_l, &
    318             tsurf_new, alb1, alb2, alb3, alb6, &
    319             emis_new, z0m, z0h, qsurf)
    320 
    321           debut_is=.false.
     271                u0, v0, u1, v1, gustiness, cdragm, &
     272                AcoefU, AcoefV, BcoefU, BcoefV, &
     273                p1lay, temp_air, &
     274                flux_u1, flux_v1)
     275
     276
     277        ! Set constants and compute some input for SISVAT
     278        ! = 1000 hPa
     279        ! and calculate incoming flux for SW and LW interval: swdown, lwdown
     280        swdown(:) = 0.0
     281        lwdown(:) = 0.0
     282        snow_cont_air(:) = 0.  ! the snow content in air is not a prognostic variable of the model
     283        alb_soil(:) = 0.4 ! before albedo(:) but here it is the ice albedo that we have to set
     284        ustar(:) = 0.
     285        pref = 100000.
     286        DO i = 1, knon
     287          swdown(i) = swnet(i) / (1 - albedo(i))
     288          lwdown(i) = lwdownm(i)
     289          wind_velo(i) = u1(i)**2 + v1(i)**2
     290          wind_velo(i) = wind_velo(i)**0.5
     291          pexner(i) = (p1lay(i) / pref)**(RD / RCPD)
     292          dens_air(i) = p1lay(i) / RD / temp_air(i)  ! dry air density
     293          zsl_height(i) = pphi1(i) / RG
     294          tsoil0(i, :) = tsoil(i, :)
     295          ustar(i) = (cdragm(i) * (wind_velo(i)**2))**0.5
     296        END DO
     297
     298        dtis = dtime
     299
     300        IF (lafin) THEN
     301          lafin_is = .true.
     302        END IF
     303
     304        CALL surf_inlandsis(knon, rlon, rlat, knindex, itime, dtis, debut_is, lafin_is, &
     305                rmu0, swdown, lwdown, albedo, pexner, ps, p1lay, precip_rain, precip_snow, &
     306                zsl_height, wind_velo, ustar, temp_air, dens_air, spechum, tsurf, &
     307                rugoro, snow_cont_air, alb_soil, alt, slope, cloudf, &
     308                radsol, qsol, tsoil0, snow, zfra, snowhgt, qsnow, to_ice, sissnow, agesno, &
     309                AcoefH, AcoefQ, BcoefH, BcoefQ, cdragm, cdragh, &
     310                run_off_lic, fqfonte, ffonte, evap, erod, fluxsens, fluxlat, dflux_s, dflux_l, &
     311                tsurf_new, alb1, alb2, alb3, alb6, &
     312                emis_new, z0m, z0h, qsurf)
     313
     314        debut_is = .false.
    322315
    323316
    … …  
    325318
    326319        ! for consistency with standard LMDZ, add calving to run_off_lic
    327         run_off_lic(:)=run_off_lic(:) + to_ice(:)
     320        run_off_lic(:) = run_off_lic(:) + to_ice(:)
    328321
    329322        DO i = 1, knon
    330            ffonte_global(knindex(i),is_lic)    = ffonte(i)
    331            fqfonte_global(knindex(i),is_lic)  = fqfonte(i)! net melting= melting - refreezing
    332            fqcalving_global(knindex(i),is_lic) = to_ice(i) ! flux
    333            runofflic_global(knindex(i)) = run_off_lic(i)
     323          ffonte_global(knindex(i), is_lic) = ffonte(i)
     324          fqfonte_global(knindex(i), is_lic) = fqfonte(i)! net melting= melting - refreezing
     325          fqcalving_global(knindex(i), is_lic) = to_ice(i) ! flux
     326          runofflic_global(knindex(i)) = run_off_lic(i)
    334327        ENDDO
    335328        ! Here, we assume that the calving term is equal to the to_ice term
    336329        ! (no ice accumulation)
    337330
    338 
    339 #else
    340        abort_message='Pb de coherence: landice_opt = 1  mais CPP_INLANDSIS = .false.'
    341        CALL abort_physic(modname,abort_message,1)
    342 #endif
    343 
    344 
    345     ELSE
    346 
    347 !****************************************************************************************
    348 ! Soil calculations
    349 
    350 !****************************************************************************************
    351 
    352     ! EV: use calbeta
    353     CALL calbeta(dtime, is_lic, knon, snow, qsol, beta, cal, dif_grnd)
    354 
    355 
    356     ! use soil model and recalculate properly cal
    357     IF (soil_model) THEN
    358        CALL soil(dtime, is_lic, knon, snow, tsurf, qsol, &
    359    longitude(knindex(1:knon)), latitude(knindex(1:knon)), tsoil, soilcap, soilflux)
    360        cal(1:knon) = RCPD / soilcap(1:knon)
    361        radsol(1:knon)  = radsol(1:knon) + soilflux(1:knon)
    362     ELSE
    363        cal = RCPD * calice
    364        WHERE (snow > 0.0) cal = RCPD * calsno
    365     ENDIF
    366 
    367 
    368 !****************************************************************************************
    369 ! Calulate fluxes
    370 
    371 !****************************************************************************************
    372 !    beta(:) = 1.0
    373 !    dif_grnd(:) = 0.0
    374 
    375 ! Suppose zero surface speed
    376     u0(:)=0.0
    377     v0(:)=0.0
    378     u1_lay(:) = u1(:) - u0(:)
    379     v1_lay(:) = v1(:) - v0(:)
    380 
    381     CALL calcul_fluxs(knon, is_lic, dtime, &
    382          tsurf, p1lay, cal, beta, cdragh, cdragh, ps, &
    383          precip_rain, precip_snow, snow, qsurf,  &
    384          radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, gustiness, &
    385          1.,AcoefH, AcoefQ, BcoefH, BcoefQ, &
    386          tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)
     331      ELSE
     332        abort_message = 'Pb de coherence: landice_opt = 1  mais CPP_INLANDSIS = .false.'
     333        CALL abort_physic(modname, abort_message, 1)
     334      END IF
     335
     336    ELSE
     337
     338      !****************************************************************************************
     339      ! Soil calculations
     340
     341      !****************************************************************************************
     342
     343      ! EV: use calbeta
     344      CALL calbeta(dtime, is_lic, knon, snow, qsol, beta, cal, dif_grnd)
     345
     346
     347      ! use soil model and recalculate properly cal
     348      IF (soil_model) THEN
     349        CALL soil(dtime, is_lic, knon, snow, tsurf, qsol, &
     350                longitude(knindex(1:knon)), latitude(knindex(1:knon)), tsoil, soilcap, soilflux)
     351        cal(1:knon) = RCPD / soilcap(1:knon)
     352        radsol(1:knon) = radsol(1:knon) + soilflux(1:knon)
     353      ELSE
     354        cal = RCPD * calice
     355        WHERE (snow > 0.0) cal = RCPD * calsno
     356      ENDIF
     357
     358
     359      !****************************************************************************************
     360      ! Calulate fluxes
     361
     362      !****************************************************************************************
     363      !    beta(:) = 1.0
     364      !    dif_grnd(:) = 0.0
     365
     366      ! Suppose zero surface speed
     367      u0(:) = 0.0
     368      v0(:) = 0.0
     369      u1_lay(:) = u1(:) - u0(:)
     370      v1_lay(:) = v1(:) - v0(:)
     371
     372      CALL calcul_fluxs(knon, is_lic, dtime, &
     373              tsurf, p1lay, cal, beta, cdragh, cdragh, ps, &
     374              precip_rain, precip_snow, snow, qsurf, &
     375              radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, gustiness, &
     376              1., AcoefH, AcoefQ, BcoefH, BcoefQ, &
     377              tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)
    387378
    388379#ifdef ISO
    … …  
    411402#endif         
    412403
    413     CALL calcul_flux_wind(knon, dtime, &
    414          u0, v0, u1, v1, gustiness, cdragm, &
    415          AcoefU, AcoefV, BcoefU, BcoefV, &
    416          p1lay, temp_air, &
    417          flux_u1, flux_v1)
    418 
    419 
    420 !****************************************************************************************
    421 ! Calculate albedo
    422 
    423 !****************************************************************************************
    424 
    425 !IM: plusieurs choix/tests sur l'albedo des "glaciers continentaux"
    426 !       alb1(1 : knon)  = 0.6 !IM cf FH/GK
    427 !       alb1(1 : knon)  = 0.82
    428 !       alb1(1 : knon)  = 0.77 !211003 Ksta0.77
    429 !       alb1(1 : knon)  = 0.8 !KstaTER0.8 & LMD_ARMIP5
    430 !IM: KstaTER0.77 & LMD_ARMIP6   
    431 
    432 ! Attantion: alb1 and alb2 are not the same!
    433     alb1(1:knon)  = alb_vis_sno_lic
    434     alb2(1:knon)  = alb_nir_sno_lic
    435 
    436 
    437 !****************************************************************************************
    438 ! Rugosity
    439 
    440 !****************************************************************************************
    441 
    442 if (z0m_landice > 0.) then
    443     z0m(1:knon) = z0m_landice
    444     z0h(1:knon) = z0h_landice
    445 else
    446     ! parameterization of z0=f(T) following measurements in Adelie Land by Amory et al 2018
    447     coefa = 0.1658 !0.1862 !Ant
    448     coefb = -50.3869 !-55.7718 !Ant
    449     ta1 = 253.15 !255. Ant
    450     ta2 = 273.15
    451     ta3 = 273.15+3
    452     z01 = exp(coefa*ta1 + coefb) !~0.2 ! ~0.25 mm
    453     z02 = exp(coefa*ta2 + coefb) !~6  !~7 mm
    454     z03 = z01
    455     coefc = log(z03/z02)/(ta3-ta2)
    456     coefd = log(z03)-coefc*ta3
    457     do j=1,knon
    458       if (temp_air(j) < ta1) then
    459         z0m(j) = z01
    460       else if (temp_air(j)>=ta1 .and. temp_air(j)<ta2) then
    461         z0m(j) = exp(coefa*temp_air(j) + coefb)
    462       else if (temp_air(j)>=ta2 .and. temp_air(j)<ta3) then
    463         ! if st > 0, melting induce smooth surface
    464         z0m(j) = exp(coefc*temp_air(j) + coefd)
     404      CALL calcul_flux_wind(knon, dtime, &
     405              u0, v0, u1, v1, gustiness, cdragm, &
     406              AcoefU, AcoefV, BcoefU, BcoefV, &
     407              p1lay, temp_air, &
     408              flux_u1, flux_v1)
     409
     410
     411      !****************************************************************************************
     412      ! Calculate albedo
     413
     414      !****************************************************************************************
     415
     416      !IM: plusieurs choix/tests sur l'albedo des "glaciers continentaux"
     417      !       alb1(1 : knon)  = 0.6 !IM cf FH/GK
     418      !       alb1(1 : knon)  = 0.82
     419      !       alb1(1 : knon)  = 0.77 !211003 Ksta0.77
     420      !       alb1(1 : knon)  = 0.8 !KstaTER0.8 & LMD_ARMIP5
     421      !IM: KstaTER0.77 & LMD_ARMIP6
     422
     423      ! Attantion: alb1 and alb2 are not the same!
     424      alb1(1:knon) = alb_vis_sno_lic
     425      alb2(1:knon) = alb_nir_sno_lic
     426
     427
     428      !****************************************************************************************
     429      ! Rugosity
     430
     431      !****************************************************************************************
     432
     433      if (z0m_landice > 0.) then
     434        z0m(1:knon) = z0m_landice
     435        z0h(1:knon) = z0h_landice
    465436      else
    466         z0m(j) = z03
     437        ! parameterization of z0=f(T) following measurements in Adelie Land by Amory et al 2018
     438        coefa = 0.1658 !0.1862 !Ant
     439        coefb = -50.3869 !-55.7718 !Ant
     440        ta1 = 253.15 !255. Ant
     441        ta2 = 273.15
     442        ta3 = 273.15 + 3
     443        z01 = exp(coefa * ta1 + coefb) !~0.2 ! ~0.25 mm
     444        z02 = exp(coefa * ta2 + coefb) !~6  !~7 mm
     445        z03 = z01
     446        coefc = log(z03 / z02) / (ta3 - ta2)
     447        coefd = log(z03) - coefc * ta3
     448        do j = 1, knon
     449          if (temp_air(j) < ta1) then
     450            z0m(j) = z01
     451          else if (temp_air(j)>=ta1 .and. temp_air(j)<ta2) then
     452            z0m(j) = exp(coefa * temp_air(j) + coefb)
     453          else if (temp_air(j)>=ta2 .and. temp_air(j)<ta3) then
     454            ! if st > 0, melting induce smooth surface
     455            z0m(j) = exp(coefc * temp_air(j) + coefd)
     456          else
     457            z0m(j) = z03
     458          endif
     459          z0h(j) = z0m(j)
     460        enddo
     461
    467462      endif
    468       z0h(j)=z0m(j)
    469     enddo
    470 
    471 endif   
    472  
    473 
    474 !****************************************************************************************
    475 ! Simple blowing snow param
    476 !****************************************************************************************
    477 ! we proceed in 2 steps:
    478 ! first we erode - if possible -the accumulated snow during the time step
    479 ! then we update the density of the underlying layer and see if we can also erode
    480 ! this layer
    481 
    482 
    483    if (ok_bs) then
    484        fluxbs(:)=0.
    485        do j=1,knon
    486           ws1(j)=(u1(j)**2+v1(j)**2)**0.5
    487           ustar(j)=(cdragm(j)*(u1(j)**2+v1(j)**2))**0.5
    488           rhod(j)=p1lay(j)/RD/temp_air(j)
    489           ustart0(j) =(log(2.868)-log(1.625))/0.085*sqrt(cdragm(j))
    490        enddo
    491 
    492        ! 1st step: erosion of fresh snow accumulated during the time step
    493        do j=1, knon
    494        if (precip_snow(j) > 0.) then
    495            rhos(j)=rhofresh_bs
    496            ! blowing snow flux formula used in MAR
    497            ustart(j)=ustart0(j)*exp(max(rhoice_bs/rhofresh_bs-rhoice_bs/rhos(j),0.))*exp(max(0.,rhos(j)-rhohard_bs))
    498            ! we have multiplied by exp to prevent erosion when rhos>rhohard_bs
    499            ! computation of qbs at the top of the saltation layer
    500            ! default formulation from MAR model (Amory et al. 2021, Gallee et al. 2001)
    501            esalt=1./(c_esalt_bs*max(1.e-6,ustar(j)))
    502            hsalt(j)=0.08436*(max(1.e-6,ustar(j))**1.27)
    503            qsalt(j)=(max(ustar(j)**2-ustart(j)**2,0.))/(RG*hsalt(j))*esalt
    504            ! calculation of erosion (flux positive towards the surface here)
    505            ! consistent with implicit resolution of turbulent mixing equation
    506            ! Nemoto and Nishimura 2004 show that steady-state saltation is achieved within a time tau_eqsalt_bs of about 10s
    507            ! we thus prevent snowerosion (snow particle transfer from the saltation layer to the first model level)
    508            ! integrated over tau_eqsalt_bs to exceed the total mass of snow particle in the saltation layer
    509            ! (rho*qsalt*hsalt)
    510            ! during this first step we also lower bound the erosion to the amount of fresh snow accumulated during the time step
    511            maxerosion=min(precip_snow(j),hsalt(j)*qsalt(j)*rhod(j)/tau_eqsalt_bs)
    512 
    513            fluxbs_1(j)=rhod(j)*ws1(j)*cdragh(j)*zeta_bs*(AcoefQBS(j)-qsalt(j)) &
    514                    / (1.-rhod(j)*ws1(j)*cdragh(j)*zeta_bs*BcoefQBS(j)*dtime)
    515            fluxbs_1(j)=max(-maxerosion,fluxbs_1(j))
    516 
    517            if (precip_snow(j) > abs(fluxbs_1(j))) then
    518                ok_remaining_freshsnow(j)=.true.
    519                bsweight_fresh(j)=1.
    520            else
    521                ok_remaining_freshsnow(j)=.false.
    522                bsweight_fresh(j)=exp(-(abs(fluxbs_1(j))-precip_snow(j))/precip_snow(j))
    523            endif
    524        else
    525            ok_remaining_freshsnow(j)=.false.
    526            fluxbs_1(j)=0.
    527            bsweight_fresh(j)=0.
    528        endif
    529        enddo
    530 
    531 
    532        ! we now compute the snow age of the overlying layer (snow surface after erosion of the fresh snow accumulated during the time step)
    533        ! this is done through the routine albsno
    534        CALL albsno(klon,knon,dtime,agesno(:),alb_neig(:), precip_snow(:)+fluxbs_1(:))
    535 
    536        ! 2nd step:
    537        ! computation of threshold friction velocity
    538        ! which depends on surface snow density
    539        do j = 1, knon
    540         if (ok_remaining_freshsnow(j)) then
    541            fluxbs_2(j)=0.
    542         else
    543            ! we start eroding the underlying layer
    544            ! estimation of snow density
    545            ! snow density increases with snow age and
    546            ! increases even faster in case of sedimentation of blowing snow or rain
    547            tau_dens=max(tau_densmin_bs, tau_dens0_bs*exp(-abs(precip_bs(j))/pbst_bs - &
    548                     abs(precip_rain(j))/prt_bs)*exp(-max(tsurf(j)-RTT,0.)))
    549            rhos(j)=rhofresh_bs+(rhohard_bs-rhofresh_bs)*(1.-exp(-agesno(j)*86400.0/tau_dens))
    550            ! blowing snow flux formula used in MAR
    551            ustart(j)=ustart0(j)*exp(max(rhoice_bs/rhofresh_bs-rhoice_bs/rhos(j),0.))*exp(max(0.,rhos(j)-rhohard_bs))
    552            ! we have multiplied by exp to prevent erosion when rhos>rhohard_bs
    553            ! computation of qbs at the top of the saltation layer
    554            ! default formulation from MAR model (Amory et al. 2021, Gallee et al. 2001)
    555            esalt=1./(c_esalt_bs*max(1.e-6,ustar(j)))
    556            hsalt(j)=0.08436*(max(1.e-6,ustar(j))**1.27)
    557            qsalt(j)=(max(ustar(j)**2-ustart(j)**2,0.))/(RG*hsalt(j))*esalt
    558            ! calculation of erosion (flux positive towards the surface here)
    559            ! consistent with implicit resolution of turbulent mixing equation
    560            ! Nemoto and Nishimura 2004 show that steady-state saltation is achieved within a time tau_eqsalt_bs of about 10s
    561            ! we thus prevent snowerosion (snow particle transfer from the saltation layer to the first model level)
    562            ! integrated over tau_eqsalt_bs to exceed the total mass of snow particle in the saltation layer
    563            ! (rho*qsalt*hsalt)
    564            maxerosion=hsalt(j)*qsalt(j)*rhod(j)/tau_eqsalt_bs
    565            fluxbs_2(j)=rhod(j)*ws1(j)*cdragh(j)*zeta_bs*(AcoefQBS(j)-qsalt(j)) &
    566                    / (1.-rhod(j)*ws1(j)*cdragh(j)*zeta_bs*BcoefQBS(j)*dtime)
    567            fluxbs_2(j)=max(-maxerosion,fluxbs_2(j))
    568          endif
    569        enddo
    570 
    571 
    572 
    573 
    574        ! final flux and outputs       
    575         do j=1, knon
    576               ! total flux is the erosion of fresh snow +
    577               ! a fraction of the underlying snow (if all the fresh snow has been eroded)
    578               ! the calculation of the fraction is quite delicate since we do not know
    579               ! how much time was needed to erode the fresh snow. We assume that this time
    580               ! is dt*exp(-(abs(fluxbs1)-precipsnow)/precipsnow)=dt*bsweight_fresh
    581 
    582               fluxbs(j)=fluxbs_1(j)+fluxbs_2(j)*(1.-bsweight_fresh(j))
    583               i = knindex(j)
    584               zxustartlic(i) = ustart(j)
    585               zxrhoslic(i) = rhos(j)
    586               zxqsaltlic(i)=qsalt(j)
     463
     464
     465      !****************************************************************************************
     466      ! Simple blowing snow param
     467      !****************************************************************************************
     468      ! we proceed in 2 steps:
     469      ! first we erode - if possible -the accumulated snow during the time step
     470      ! then we update the density of the underlying layer and see if we can also erode
     471      ! this layer
     472
     473      if (ok_bs) then
     474        fluxbs(:) = 0.
     475        do j = 1, knon
     476          ws1(j) = (u1(j)**2 + v1(j)**2)**0.5
     477          ustar(j) = (cdragm(j) * (u1(j)**2 + v1(j)**2))**0.5
     478          rhod(j) = p1lay(j) / RD / temp_air(j)
     479          ustart0(j) = (log(2.868) - log(1.625)) / 0.085 * sqrt(cdragm(j))
    587480        enddo
    588481
    589 
    590   else ! not ok_bs
    591   ! those lines are useful to calculate the snow age
    592        CALL albsno(klon,knon,dtime,agesno(:),alb_neig(:), precip_snow(:))
    593 
    594   endif ! if ok_bs
    595 
    596 
    597 
    598 !****************************************************************************************
    599 ! Calculate snow amount
    600 
    601 !****************************************************************************************
    602     IF (ok_bs) THEN
    603       precip_totsnow(:)=precip_snow(:)+precip_bs(:)
    604       evap_totsnow(:)=evap(:)-fluxbs(:) ! flux bs is positive towards the surface (snow erosion)
    605     ELSE
    606       precip_totsnow(:)=precip_snow(:)
    607       evap_totsnow(:)=evap(:)
    608     ENDIF
    609    
    610  
    611     CALL fonte_neige(knon, is_lic, knindex, dtime, &
    612          tsurf, precip_rain, precip_totsnow, &
    613          snow, qsol, tsurf_new, evap_totsnow &
     482        ! 1st step: erosion of fresh snow accumulated during the time step
     483        do j = 1, knon
     484          if (precip_snow(j) > 0.) then
     485            rhos(j) = rhofresh_bs
     486            ! blowing snow flux formula used in MAR
     487            ustart(j) = ustart0(j) * exp(max(rhoice_bs / rhofresh_bs - rhoice_bs / rhos(j), 0.)) * exp(max(0., rhos(j) - rhohard_bs))
     488            ! we have multiplied by exp to prevent erosion when rhos>rhohard_bs
     489            ! computation of qbs at the top of the saltation layer
     490            ! default formulation from MAR model (Amory et al. 2021, Gallee et al. 2001)
     491            esalt = 1. / (c_esalt_bs * max(1.e-6, ustar(j)))
     492            hsalt(j) = 0.08436 * (max(1.e-6, ustar(j))**1.27)
     493            qsalt(j) = (max(ustar(j)**2 - ustart(j)**2, 0.)) / (RG * hsalt(j)) * esalt
     494            ! calculation of erosion (flux positive towards the surface here)
     495            ! consistent with implicit resolution of turbulent mixing equation
     496            ! Nemoto and Nishimura 2004 show that steady-state saltation is achieved within a time tau_eqsalt_bs of about 10s
     497            ! we thus prevent snowerosion (snow particle transfer from the saltation layer to the first model level)
     498            ! integrated over tau_eqsalt_bs to exceed the total mass of snow particle in the saltation layer
     499            ! (rho*qsalt*hsalt)
     500            ! during this first step we also lower bound the erosion to the amount of fresh snow accumulated during the time step
     501            maxerosion = min(precip_snow(j), hsalt(j) * qsalt(j) * rhod(j) / tau_eqsalt_bs)
     502
     503            fluxbs_1(j) = rhod(j) * ws1(j) * cdragh(j) * zeta_bs * (AcoefQBS(j) - qsalt(j)) &
     504                    / (1. - rhod(j) * ws1(j) * cdragh(j) * zeta_bs * BcoefQBS(j) * dtime)
     505            fluxbs_1(j) = max(-maxerosion, fluxbs_1(j))
     506
     507            if (precip_snow(j) > abs(fluxbs_1(j))) then
     508              ok_remaining_freshsnow(j) = .true.
     509              bsweight_fresh(j) = 1.
     510            else
     511              ok_remaining_freshsnow(j) = .false.
     512              bsweight_fresh(j) = exp(-(abs(fluxbs_1(j)) - precip_snow(j)) / precip_snow(j))
     513            endif
     514          else
     515            ok_remaining_freshsnow(j) = .false.
     516            fluxbs_1(j) = 0.
     517            bsweight_fresh(j) = 0.
     518          endif
     519        enddo
     520
     521
     522        ! we now compute the snow age of the overlying layer (snow surface after erosion of the fresh snow accumulated during the time step)
     523        ! this is done through the routine albsno
     524        CALL albsno(klon, knon, dtime, agesno(:), alb_neig(:), precip_snow(:) + fluxbs_1(:))
     525
     526        ! 2nd step:
     527        ! computation of threshold friction velocity
     528        ! which depends on surface snow density
     529        do j = 1, knon
     530          if (ok_remaining_freshsnow(j)) then
     531            fluxbs_2(j) = 0.
     532          else
     533            ! we start eroding the underlying layer
     534            ! estimation of snow density
     535            ! snow density increases with snow age and
     536            ! increases even faster in case of sedimentation of blowing snow or rain
     537            tau_dens = max(tau_densmin_bs, tau_dens0_bs * exp(-abs(precip_bs(j)) / pbst_bs - &
     538                    abs(precip_rain(j)) / prt_bs) * exp(-max(tsurf(j) - RTT, 0.)))
     539            rhos(j) = rhofresh_bs + (rhohard_bs - rhofresh_bs) * (1. - exp(-agesno(j) * 86400.0 / tau_dens))
     540            ! blowing snow flux formula used in MAR
     541            ustart(j) = ustart0(j) * exp(max(rhoice_bs / rhofresh_bs - rhoice_bs / rhos(j), 0.)) * exp(max(0., rhos(j) - rhohard_bs))
     542            ! we have multiplied by exp to prevent erosion when rhos>rhohard_bs
     543            ! computation of qbs at the top of the saltation layer
     544            ! default formulation from MAR model (Amory et al. 2021, Gallee et al. 2001)
     545            esalt = 1. / (c_esalt_bs * max(1.e-6, ustar(j)))
     546            hsalt(j) = 0.08436 * (max(1.e-6, ustar(j))**1.27)
     547            qsalt(j) = (max(ustar(j)**2 - ustart(j)**2, 0.)) / (RG * hsalt(j)) * esalt
     548            ! calculation of erosion (flux positive towards the surface here)
     549            ! consistent with implicit resolution of turbulent mixing equation
     550            ! Nemoto and Nishimura 2004 show that steady-state saltation is achieved within a time tau_eqsalt_bs of about 10s
     551            ! we thus prevent snowerosion (snow particle transfer from the saltation layer to the first model level)
     552            ! integrated over tau_eqsalt_bs to exceed the total mass of snow particle in the saltation layer
     553            ! (rho*qsalt*hsalt)
     554            maxerosion = hsalt(j) * qsalt(j) * rhod(j) / tau_eqsalt_bs
     555            fluxbs_2(j) = rhod(j) * ws1(j) * cdragh(j) * zeta_bs * (AcoefQBS(j) - qsalt(j)) &
     556                    / (1. - rhod(j) * ws1(j) * cdragh(j) * zeta_bs * BcoefQBS(j) * dtime)
     557            fluxbs_2(j) = max(-maxerosion, fluxbs_2(j))
     558          endif
     559        enddo
     560
     561
     562
     563
     564        ! final flux and outputs
     565        do j = 1, knon
     566          ! total flux is the erosion of fresh snow +
     567          ! a fraction of the underlying snow (if all the fresh snow has been eroded)
     568          ! the calculation of the fraction is quite delicate since we do not know
     569          ! how much time was needed to erode the fresh snow. We assume that this time
     570          ! is dt*exp(-(abs(fluxbs1)-precipsnow)/precipsnow)=dt*bsweight_fresh
     571
     572          fluxbs(j) = fluxbs_1(j) + fluxbs_2(j) * (1. - bsweight_fresh(j))
     573          i = knindex(j)
     574          zxustartlic(i) = ustart(j)
     575          zxrhoslic(i) = rhos(j)
     576          zxqsaltlic(i) = qsalt(j)
     577        enddo
     578
     579      else ! not ok_bs
     580        ! those lines are useful to calculate the snow age
     581        CALL albsno(klon, knon, dtime, agesno(:), alb_neig(:), precip_snow(:))
     582
     583      endif ! if ok_bs
     584
     585
     586
     587      !****************************************************************************************
     588      ! Calculate snow amount
     589
     590      !****************************************************************************************
     591      IF (ok_bs) THEN
     592        precip_totsnow(:) = precip_snow(:) + precip_bs(:)
     593        evap_totsnow(:) = evap(:) - fluxbs(:) ! flux bs is positive towards the surface (snow erosion)
     594      ELSE
     595        precip_totsnow(:) = precip_snow(:)
     596        evap_totsnow(:) = evap(:)
     597      ENDIF
     598
     599      CALL fonte_neige(knon, is_lic, knindex, dtime, &
     600      tsurf, precip_rain, precip_totsnow, &
     601      snow, qsol, tsurf_new, evap_totsnow &
    614602#ifdef ISO   
    615603    ,fq_fonte_diag,fqfonte_diag,snow_evap_diag,fqcalving_diag     &
    616604    ,max_eau_sol_diag,runoff_diag,run_off_lic_diag,coeff_rel_diag &
    617605#endif
    618      )
     606      )
    619607
    620608
    … …  
    641629
    642630#endif
    643    
    644     WHERE (snow(1 : knon) < 0.0001) agesno(1 : knon) = 0.
    645     zfra(1:knon) = MAX(0.0,MIN(1.0,snow(1:knon)/(snow(1:knon)+10.0))) 
    646 
     631
     632      WHERE (snow(1:knon) < 0.0001) agesno(1:knon) = 0.
     633      zfra(1:knon) = MAX(0.0, MIN(1.0, snow(1:knon) / (snow(1:knon) + 10.0)))
    647634
    648635    END IF ! landice_opt
    649636
    650637
    651 !****************************************************************************************
    652 ! Send run-off on land-ice to coupler if coupled ocean.
    653 ! run_off_lic has been calculated in fonte_neige or surf_inlandsis
    654 ! If landice_opt>=2, corresponding call is done from surf_land_orchidee
    655 !****************************************************************************************
     638    !****************************************************************************************
     639    ! Send run-off on land-ice to coupler if coupled ocean.
     640    ! run_off_lic has been calculated in fonte_neige or surf_inlandsis
     641    ! If landice_opt>=2, corresponding call is done from surf_land_orchidee
     642    !****************************************************************************************
    656643    IF (type_ocean=='couple' .AND. landice_opt < 2) THEN
    657        ! Compress fraction where run_off_lic is active (here all pctsrf(is_lic))
    658        run_off_lic_frac(:)=0.0
    659        DO j = 1, knon
    660           i = knindex(j)
    661           run_off_lic_frac(j) = pctsrf(i,is_lic)
    662        ENDDO
    663 
    664        CALL cpl_send_landice_fields(itime, knon, knindex, run_off_lic, run_off_lic_frac)
     644      ! Compress fraction where run_off_lic is active (here all pctsrf(is_lic))
     645      run_off_lic_frac(:) = 0.0
     646      DO j = 1, knon
     647        i = knindex(j)
     648        run_off_lic_frac(j) = pctsrf(i, is_lic)
     649      ENDDO
     650
     651      CALL cpl_send_landice_fields(itime, knon, knindex, run_off_lic, run_off_lic_frac)
    665652    ENDIF
    666653
    667  ! transfer runoff rate [kg/m2/s](!) to physiq for output
    668     runoff(1:knon)=run_off_lic(1:knon)/dtime
    669 
    670        snow_o=0.
    671        zfra_o = 0.
    672        DO j = 1, knon
    673            i = knindex(j)
    674            snow_o(i) = snow(j)
    675            zfra_o(i) = zfra(j)
    676        ENDDO
    677 
    678 
    679 !albedo SB >>>
    680      select case(NSW)
    681      case(2)
    682        alb_dir(1:knon,1)=alb1(1:knon)
    683        alb_dir(1:knon,2)=alb2(1:knon)
    684      case(4)
    685        alb_dir(1:knon,1)=alb1(1:knon)
    686        alb_dir(1:knon,2)=alb2(1:knon)
    687        alb_dir(1:knon,3)=alb2(1:knon)
    688        alb_dir(1:knon,4)=alb2(1:knon)
    689      case(6)
    690        alb_dir(1:knon,1)=alb1(1:knon)
    691        alb_dir(1:knon,2)=alb1(1:knon)
    692        alb_dir(1:knon,3)=alb1(1:knon)
    693        alb_dir(1:knon,4)=alb2(1:knon)
    694        alb_dir(1:knon,5)=alb2(1:knon)
    695        alb_dir(1:knon,6)=alb2(1:knon)
    696 
    697        IF ((landice_opt == 1) .AND. (iflag_albcalc == 2)) THEN
    698        alb_dir(1:knon,1)=alb6(1:knon,1)
    699        alb_dir(1:knon,2)=alb6(1:knon,2)
    700        alb_dir(1:knon,3)=alb6(1:knon,3)
    701        alb_dir(1:knon,4)=alb6(1:knon,4)
    702        alb_dir(1:knon,5)=alb6(1:knon,5)
    703        alb_dir(1:knon,6)=alb6(1:knon,6)
    704        ENDIF
    705 
    706      end select
    707 alb_dif=alb_dir
    708 !albedo SB <<<
    709 
     654    ! transfer runoff rate [kg/m2/s](!) to physiq for output
     655    runoff(1:knon) = run_off_lic(1:knon) / dtime
     656
     657    snow_o = 0.
     658    zfra_o = 0.
     659    DO j = 1, knon
     660      i = knindex(j)
     661      snow_o(i) = snow(j)
     662      zfra_o(i) = zfra(j)
     663    ENDDO
     664
     665
     666    !albedo SB >>>
     667    select case(NSW)
     668    case(2)
     669      alb_dir(1:knon, 1) = alb1(1:knon)
     670      alb_dir(1:knon, 2) = alb2(1:knon)
     671    case(4)
     672      alb_dir(1:knon, 1) = alb1(1:knon)
     673      alb_dir(1:knon, 2) = alb2(1:knon)
     674      alb_dir(1:knon, 3) = alb2(1:knon)
     675      alb_dir(1:knon, 4) = alb2(1:knon)
     676    case(6)
     677      alb_dir(1:knon, 1) = alb1(1:knon)
     678      alb_dir(1:knon, 2) = alb1(1:knon)
     679      alb_dir(1:knon, 3) = alb1(1:knon)
     680      alb_dir(1:knon, 4) = alb2(1:knon)
     681      alb_dir(1:knon, 5) = alb2(1:knon)
     682      alb_dir(1:knon, 6) = alb2(1:knon)
     683
     684      IF ((landice_opt == 1) .AND. (iflag_albcalc == 2)) THEN
     685        alb_dir(1:knon, 1) = alb6(1:knon, 1)
     686        alb_dir(1:knon, 2) = alb6(1:knon, 2)
     687        alb_dir(1:knon, 3) = alb6(1:knon, 3)
     688        alb_dir(1:knon, 4) = alb6(1:knon, 4)
     689        alb_dir(1:knon, 5) = alb6(1:knon, 5)
     690        alb_dir(1:knon, 6) = alb6(1:knon, 6)
     691      ENDIF
     692
     693    end select
     694    alb_dif = alb_dir
     695    !albedo SB <<<
    710696
    711697  END SUBROUTINE surf_landice
    712698
    713 !****************************************************************************************
     699  !****************************************************************************************
    714700
    715701END MODULE surf_landice_mod
Note: See TracChangeset for help on using the changeset viewer.