source: LMDZ6/branches/Amaury_dev/libf/phylmd/lmdz_lscp.F90 @ 5105

Last change on this file since 5105 was 5105, checked in by abarral, 8 weeks ago

Replace 1DUTILS.h by module lmdz_1dutils.f90
Replace 1DConv.h by module lmdz_old_1dconv.f90 (it's only used by old_* files)
Convert *.F to *.f90
Fix gradsdef.h formatting
Remove unnecessary "RETURN" at the end of functions/subroutines
File size: 64.5 KB
Line 
1MODULE lmdz_lscp
2
3IMPLICIT NONE
4
5CONTAINS
6
7!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
8SUBROUTINE lscp(klon,klev,dtime,missing_val,            &
9     paprs,pplay,temp,qt,qice_save,ptconv,ratqs,        &
10     d_t, d_q, d_ql, d_qi, rneb, rneblsvol, rneb_seri,  &
11     pfraclr, pfracld,                                  &
12     cldfraliq, sigma2_icefracturb,mean_icefracturb,    &
13     radocond, radicefrac, rain, snow,                  &
14     frac_impa, frac_nucl, beta,                        &
15     prfl, psfl, rhcl, qta, fraca,                      &
16     tv, pspsk, tla, thl, iflag_cld_th,                 &
17     iflag_ice_thermo, ok_ice_sursat, qsatl, qsats,     &
18     distcltop, temp_cltop, tke, tke_dissip,            &
19     qclr, qcld, qss, qvc, rnebclr, rnebss, gamma_ss,   &
20     Tcontr, qcontr, qcontr2, fcontrN, fcontrP,         &
21     cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv, &
22     qraindiag, qsnowdiag, dqreva, dqssub, dqrauto,     &
23     dqrcol, dqrmelt, dqrfreez, dqsauto, dqsagg, dqsrim,&
24     dqsmelt, dqsfreez)
25
26!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
27! Authors: Z.X. Li (LMD), J-L Dufresne (LMD), C. Rio (LMD), J-Y Grandpeix (LMD)
28!          A. JAM (LMD), J-B Madeleine (LMD), E. Vignon (LMD), L. Touzze-Peiffert (LMD)
29!--------------------------------------------------------------------------------
30! Date: 01/2021
31!--------------------------------------------------------------------------------
32! Aim: Large Scale Clouds and Precipitation (LSCP)
33
34! This code is a new version of the fisrtilp.F90 routine, which is itself a
35! merge of 'first' (superrsaturation physics, P. LeVan K. Laval)
36! and 'ilp' (il pleut, L. Li)
37
38! Compared to the original fisrtilp code, lscp
39! -> assumes thermcep = .TRUE. all the time (fisrtilp inconsistent when .FALSE.)
40! -> consider always precipitation thermalisation (fl_cor_ebil>0)
41! -> option iflag_fisrtilp_qsat<0 no longer possible (qsat does not evolve with T)
42! -> option "oldbug" by JYG has been removed
43! -> iflag_t_glace >0 always
44! -> the 'all or nothing' cloud approach is no longer available (cpartiel=T always)
45! -> rectangular distribution from L. Li no longer available
46! -> We always account for the Wegener-Findeisen-Bergeron process (iflag_bergeron = 2 in fisrt)
47!--------------------------------------------------------------------------------
48! References:
49
50! - Bony, S., & Emanuel, K. A. 2001, JAS, doi: 10.1175/1520-0469(2001)058<3158:APOTCA>2.0.CO;2
51! - Hourdin et al. 2013, Clim Dyn, doi:10.1007/s00382-012-1343-y
52! - Jam et al. 2013, Boundary-Layer Meteorol, doi:10.1007/s10546-012-9789-3
53! - Jouhaud, et al. 2018. JAMES, doi:10.1029/2018MS001379
54! - Madeleine et al. 2020, JAMES, doi:10.1029/2020MS002046
55! - Touzze-Peifert Ludo, PhD thesis, p117-124
56! -------------------------------------------------------------------------------
57! Code structure:
58
59! P0>     Thermalization of the precipitation coming from the overlying layer
60! P1>     Evaporation of the precipitation (falling from the k+1 level)
61! P2>     Cloud formation (at the k level)
62! P2.A.1> With the PDFs, calculation of cloud properties using the inital
63!         values of T and Q
64! P2.A.2> Coupling between condensed water and temperature
65! P2.A.3> Calculation of final quantities associated with cloud formation
66! P2.B>   Release of Latent heat after cloud formation
67! P3>     Autoconversion to precipitation (k-level)
68! P4>     Wet scavenging
69!------------------------------------------------------------------------------
70! Some preliminary comments (JBM) :
71
72! The cloud water that the radiation scheme sees is not the same that the cloud
73! water used in the physics and the dynamics
74
75! During the autoconversion to precipitation (P3 step), radocond (cloud water used
76! by the radiation scheme) is calculated as an average of the water that remains
77! in the cloud during the precipitation and not the water remaining at the end
78! of the time step. The latter is used in the rest of the physics and advected
79! by the dynamics.
80
81! In summary:
82
83! Radiation:
84! xflwc(newmicro)+xfiwc(newmicro) =
85!   radocond=lwcon(nc)+iwcon(nc)
86
87! Notetheless, be aware of:
88
89! radocond .NE. ocond(nc)
90! i.e.:
91! lwcon(nc)+iwcon(nc) .NE. ocond(nc)
92! but oliq+(ocond-oliq) .EQ. ocond
93! (which is not trivial)
94!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
95
96! USE de modules contenant des fonctions.
97USE lmdz_cloudth, ONLY: cloudth, cloudth_v3, cloudth_v6, cloudth_mpc
98USE lmdz_lscp_tools, ONLY: calc_qsat_ecmwf, calc_gammasat
99USE lmdz_lscp_tools, ONLY: icefrac_lscp, icefrac_lscp_turb
100USE lmdz_lscp_tools, ONLY: fallice_velocity, distance_to_cloud_top
101USE ice_sursat_mod, ONLY: ice_sursat
102USE lmdz_lscp_poprecip, ONLY: poprecip_precld, poprecip_postcld
103
104! Use du module lmdz_lscp_ini contenant les constantes
105USE lmdz_lscp_ini, ONLY: prt_level, lunout
106USE lmdz_lscp_ini, ONLY: seuil_neb, niter_lscp, iflag_evap_prec, t_coup, DDT0, ztfondue, rain_int_min
107USE lmdz_lscp_ini, ONLY: ok_radocond_snow, a_tr_sca, cld_expo_con, cld_expo_lsc
108USE lmdz_lscp_ini, ONLY: iflag_cloudth_vert, iflag_rain_incloud_vol, iflag_t_glace, t_glace_min
109USE lmdz_lscp_ini, ONLY: coef_eva, coef_sub,cld_tau_lsc, cld_tau_con, cld_lc_lsc, cld_lc_con
110USE lmdz_lscp_ini, ONLY: iflag_bergeron, iflag_fisrtilp_qsat, iflag_vice, cice_velo, dice_velo
111USE lmdz_lscp_ini, ONLY: iflag_autoconversion, ffallv_con, ffallv_lsc, min_frac_th_cld
112USE lmdz_lscp_ini, ONLY: RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT, RD, RG
113USE lmdz_lscp_ini, ONLY: ok_poprecip
114USE lmdz_lscp_ini, ONLY: iflag_icefrac
115IMPLICIT NONE
116
117!===============================================================================
118! VARIABLES DECLARATION
119!===============================================================================
120
121! INPUT VARIABLES:
122!-----------------
123
124  INTEGER,                         INTENT(IN)   :: klon,klev       ! number of horizontal grid points and vertical levels
125  REAL,                            INTENT(IN)   :: dtime           ! time step [s]
126  REAL,                            INTENT(IN)   :: missing_val     ! missing value for output
127
128  REAL, DIMENSION(klon,klev+1),    INTENT(IN)   :: paprs           ! inter-layer pressure [Pa]
129  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: pplay           ! mid-layer pressure [Pa]
130  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: temp            ! temperature (K)
131  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: qt              ! total specific humidity (in vapor phase in input) [kg/kg]
132  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: qice_save       ! ice specific from previous time step [kg/kg]
133  INTEGER,                         INTENT(IN)   :: iflag_cld_th    ! flag that determines the distribution of convective clouds
134  INTEGER,                         INTENT(IN)   :: iflag_ice_thermo! flag to activate the ice thermodynamics
135                                                                   ! CR: if iflag_ice_thermo=2, only convection is active
136  LOGICAL,                         INTENT(IN)   :: ok_ice_sursat   ! flag to determine if ice sursaturation is activated
137
138  LOGICAL, DIMENSION(klon,klev),   INTENT(IN)   :: ptconv          ! grid points where deep convection scheme is active
139
140  !Inputs associated with thermal plumes
141
142  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: tv                  ! virtual potential temperature [K]
143  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: qta                 ! specific humidity within thermals [kg/kg]
144  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: fraca               ! fraction of thermals within the mesh [-]
145  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: pspsk               ! exner potential (p/100000)**(R/cp)
146  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: tla                 ! liquid temperature within thermals [K]
147  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: tke                 !--turbulent kinetic energy [m2/s2]
148  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: tke_dissip          !--TKE dissipation [m2/s3]
149
150  ! INPUT/OUTPUT variables
151  !------------------------
152 
153  REAL, DIMENSION(klon,klev),      INTENT(INOUT)   :: thl              ! liquid potential temperature [K]
154  REAL, DIMENSION(klon,klev),      INTENT(INOUT)   :: ratqs            ! function of pressure that sets the large-scale
155
156
157  ! Input sursaturation en glace
158  REAL, DIMENSION(klon,klev),      INTENT(INOUT):: rneb_seri           ! fraction nuageuse en memoire
159 
160  ! OUTPUT variables
161  !-----------------
162
163  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: d_t                 ! temperature increment [K]
164  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: d_q                 ! specific humidity increment [kg/kg]
165  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: d_ql                ! liquid water increment [kg/kg]
166  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: d_qi                ! cloud ice mass increment [kg/kg]
167  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: rneb                ! cloud fraction [-]
168  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: rneblsvol           ! cloud fraction per unit volume [-] 
169  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: pfraclr             ! precip fraction clear-sky part [-]
170  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: pfracld             ! precip fraction cloudy part [-]
171  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: cldfraliq           ! liquid fraction of cloud [-]
172  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: sigma2_icefracturb  ! Variance of the diagnostic supersaturation distribution (icefrac_turb) [-]
173  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: mean_icefracturb    ! Mean of the diagnostic supersaturation distribution (icefrac_turb) [-]
174  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: radocond            ! condensed water used in the radiation scheme [kg/kg]
175  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: radicefrac          ! ice fraction of condensed water for radiation scheme
176  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: rhcl                ! clear-sky relative humidity [-]
177  REAL, DIMENSION(klon),           INTENT(OUT)  :: rain                ! surface large-scale rainfall [kg/s/m2]
178  REAL, DIMENSION(klon),           INTENT(OUT)  :: snow                ! surface large-scale snowfall [kg/s/m2]
179  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qsatl               ! saturation specific humidity wrt liquid [kg/kg]
180  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qsats               ! saturation specific humidity wrt ice [kg/kg] 
181  REAL, DIMENSION(klon,klev+1),    INTENT(OUT)  :: prfl                ! large-scale rainfall flux in the column [kg/s/m2]
182  REAL, DIMENSION(klon,klev+1),    INTENT(OUT)  :: psfl                ! large-scale snowfall flux in the column [kg/s/m2]
183  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: distcltop           ! distance to cloud top [m]
184  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: temp_cltop          ! temperature of cloud top [K]
185  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: beta                ! conversion rate of condensed water
186
187  ! fraction of aerosol scavenging through impaction and nucleation    (for on-line)
188 
189  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: frac_impa           ! scavenging fraction due tu impaction [-]
190  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: frac_nucl           ! scavenging fraction due tu nucleation [-]
191 
192  ! for supersaturation and contrails parameterisation
193 
194  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qclr                ! specific total water content in clear sky region [kg/kg]
195  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qcld                ! specific total water content in cloudy region [kg/kg]
196  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qss                 ! specific total water content in supersat region [kg/kg]
197  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qvc                 ! specific vapor content in clouds [kg/kg]
198  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: rnebclr             ! mesh fraction of clear sky [-]   
199  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: rnebss              ! mesh fraction of ISSR [-] 
200  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: gamma_ss            ! coefficient governing the ice nucleation RHi threshold [-]     
201  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: Tcontr              ! threshold temperature for contrail formation [K]
202  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qcontr              ! threshold humidity for contrail formation [kg/kg]
203  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qcontr2             ! // (2nd expression more consistent with LMDZ expression of q)         
204  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: fcontrN             ! fraction of grid favourable to non-persistent contrails
205  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: fcontrP             ! fraction of grid favourable to persistent contrails
206  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: cloudth_sth         ! mean saturation deficit in thermals
207  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: cloudth_senv        ! mean saturation deficit in environment
208  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: cloudth_sigmath     ! std of saturation deficit in thermals
209  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: cloudth_sigmaenv    ! std of saturation deficit in environment
210
211  ! for POPRECIP
212
213  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qraindiag      !--DIAGNOSTIC specific rain content [kg/kg]
214  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qsnowdiag      !--DIAGNOSTIC specific snow content [kg/kg]
215  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqreva         !--rain tendendy due to evaporation [kg/kg/s]
216  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqssub         !--snow tendency due to sublimation [kg/kg/s]
217  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqrcol         !--rain tendendy due to collection by rain of liquid cloud droplets [kg/kg/s]
218  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqsagg         !--snow tendency due to collection of lcoud ice by aggregation [kg/kg/s]
219  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqrauto        !--rain tendency due to autoconversion of cloud liquid [kg/kg/s]
220  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqsauto        !--snow tendency due to autoconversion of cloud ice [kg/kg/s]
221  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqsrim         !--snow tendency due to riming [kg/kg/s]
222  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqsmelt        !--snow tendency due to melting [kg/kg/s]
223  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqrmelt        !--rain tendency due to melting [kg/kg/s]
224  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqsfreez       !--snow tendency due to freezing [kg/kg/s]
225  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqrfreez       !--rain tendency due to freezing [kg/kg/s]
226
227
228  ! LOCAL VARIABLES:
229  !----------------
230  REAL,DIMENSION(klon) :: qsl, qsi                                ! saturation threshold at current vertical level
231  REAL :: zct, zcl,zexpo
232  REAL, DIMENSION(klon,klev) :: ctot
233  REAL, DIMENSION(klon,klev) :: ctot_vol
234  REAL, DIMENSION(klon) :: zqs, zdqs
235  REAL :: zdelta, zcor, zcvm5
236  REAL, DIMENSION(klon) :: zdqsdT_raw
237  REAL, DIMENSION(klon) :: gammasat,dgammasatdt                   ! coefficient to make cold condensation at the correct RH and derivative wrt T
238  REAL, DIMENSION(klon) :: Tbef,qlbef,DT                          ! temperature, humidity and temp. variation during lognormal iteration
239  REAL :: num,denom
240  REAL :: cste
241  REAL, DIMENSION(klon) :: zpdf_sig,zpdf_k,zpdf_delta             ! lognormal parameters
242  REAL, DIMENSION(klon) :: Zpdf_a,zpdf_b,zpdf_e1,zpdf_e2          ! lognormal intermediate variables
243  REAL :: erf
244  REAL, DIMENSION(klon) :: zfice_th
245  REAL, DIMENSION(klon) :: qcloud, qincloud_mpc
246  REAL, DIMENSION(klon) :: zrfl, zrfln
247  REAL :: zqev, zqevt
248  REAL, DIMENSION(klon) :: zifl, zifln, ziflprev
249  REAL :: zqev0,zqevi, zqevti
250  REAL, DIMENSION(klon) :: zoliq, zcond, zq, zqn
251  REAL, DIMENSION(klon) :: zoliql, zoliqi
252  REAL, DIMENSION(klon) :: zt
253  REAL, DIMENSION(klon,klev) :: zrho
254  REAL, DIMENSION(klon) :: zdz,iwc
255  REAL :: zchau,zfroi
256  REAL, DIMENSION(klon) :: zfice,zneb,znebprecip
257  REAL :: zmelt,zrain,zsnow,zprecip
258  REAL, DIMENSION(klon) :: dzfice
259  REAL, DIMENSION(klon) :: zfice_turb, dzfice_turb
260  REAL :: zsolid
261  REAL, DIMENSION(klon) :: qtot, qzero
262  REAL, DIMENSION(klon) :: dqsl,dqsi
263  REAL :: smallestreal
264  !  Variables for Bergeron process
265  REAL :: zcp, coef1, DeltaT, Deltaq, Deltaqprecl
266  REAL, DIMENSION(klon) :: zqpreci, zqprecl
267  ! Variables precipitation energy conservation
268  REAL, DIMENSION(klon) :: zmqc
269  REAL :: zalpha_tr
270  REAL :: zfrac_lessi
271  REAL, DIMENSION(klon) :: zprec_cond
272  REAL :: zmair
273  REAL :: zcpair, zcpeau
274  REAL, DIMENSION(klon) :: zlh_solid
275  REAL, DIMENSION(klon) :: ztupnew
276  REAL, DIMENSION(klon) :: zqvapclr, zqupnew ! for poprecip evap / subl
277  REAL :: zm_solid         ! for liquid -> solid conversion
278  REAL, DIMENSION(klon) :: zrflclr, zrflcld
279  REAL, DIMENSION(klon) :: d_zrfl_clr_cld, d_zifl_clr_cld
280  REAL, DIMENSION(klon) :: d_zrfl_cld_clr, d_zifl_cld_clr
281  REAL, DIMENSION(klon) :: ziflclr, ziflcld
282  REAL, DIMENSION(klon) :: znebprecipclr, znebprecipcld
283  REAL, DIMENSION(klon) :: tot_zneb, tot_znebn, d_tot_zneb
284  REAL, DIMENSION(klon) :: d_znebprecip_clr_cld, d_znebprecip_cld_clr
285  REAL, DIMENSION(klon,klev) :: velo
286  REAL :: vr, ffallv
287  REAL :: qlmpc, qimpc, rnebmpc
288  REAL, DIMENSION(klon,klev) :: radocondi, radocondl
289  REAL :: effective_zneb
290  REAL, DIMENSION(klon) :: zdistcltop, ztemp_cltop
291  REAL, DIMENSION(klon) :: zqliq, zqice, zqvapcl        ! for icefrac_lscp_turb
292
293  INTEGER i, k, n, kk, iter
294  INTEGER, DIMENSION(klon) :: n_i
295  INTEGER ncoreczq
296  INTEGER, DIMENSION(klon,klev) :: mpc_bl_points
297  LOGICAL iftop
298
299  LOGICAL, DIMENSION(klon) :: lognormale
300  LOGICAL, DIMENSION(klon) :: keepgoing
301
302  CHARACTER (len = 20) :: modname = 'lscp'
303  CHARACTER (len = 80) :: abort_message
304 
305
306!===============================================================================
307! INITIALISATION
308!===============================================================================
309
310! Few initial checks
311
312
313IF (iflag_fisrtilp_qsat < 0) THEN
314    abort_message = 'lscp cannot be used with iflag_fisrtilp<0'
315    CALL abort_physic(modname,abort_message,1)
316ENDIF
317
318! Few initialisations
319
320znebprecip(:)=0.0
321ctot_vol(1:klon,1:klev)=0.0
322rneblsvol(1:klon,1:klev)=0.0
323smallestreal=1.e-9
324znebprecipclr(:)=0.0
325znebprecipcld(:)=0.0
326mpc_bl_points(:,:)=0
327
328IF (prt_level>9) WRITE(lunout,*) 'NUAGES4 A. JAM'
329
330! AA for 'safety' reasons
331zalpha_tr   = 0.
332zfrac_lessi = 0.
333beta(:,:)= 0.
334
335! Initialisation of variables:
336
337prfl(:,:) = 0.0
338psfl(:,:) = 0.0
339d_t(:,:) = 0.0
340d_q(:,:) = 0.0
341d_ql(:,:) = 0.0
342d_qi(:,:) = 0.0
343rneb(:,:) = 0.0
344pfraclr(:,:)=0.0
345pfracld(:,:)=0.0
346cldfraliq(:,:)=0.
347sigma2_icefracturb(:,:)=0.
348mean_icefracturb(:,:)=0.
349radocond(:,:) = 0.0
350radicefrac(:,:) = 0.0
351frac_nucl(:,:) = 1.0
352frac_impa(:,:) = 1.0
353rain(:) = 0.0
354snow(:) = 0.0
355zoliq(:)=0.0
356zfice(:)=0.0
357dzfice(:)=0.0
358zfice_turb(:)=0.0
359dzfice_turb(:)=0.0
360zqprecl(:)=0.0
361zqpreci(:)=0.0
362zrfl(:) = 0.0
363zifl(:) = 0.0
364ziflprev(:)=0.0
365zneb(:) = seuil_neb
366zrflclr(:) = 0.0
367ziflclr(:) = 0.0
368zrflcld(:) = 0.0
369ziflcld(:) = 0.0
370tot_zneb(:) = 0.0
371tot_znebn(:) = 0.0
372d_tot_zneb(:) = 0.0
373qzero(:) = 0.0
374zdistcltop(:)=0.0
375ztemp_cltop(:) = 0.0
376ztupnew(:)=0.0
377!--ice supersaturation
378gamma_ss(:,:) = 1.0
379qss(:,:) = 0.0
380rnebss(:,:) = 0.0
381Tcontr(:,:) = missing_val
382qcontr(:,:) = missing_val
383qcontr2(:,:) = missing_val
384fcontrN(:,:) = 0.0
385fcontrP(:,:) = 0.0
386distcltop(:,:)=0.
387temp_cltop(:,:)=0.
388
389!-- poprecip
390qraindiag(:,:)= 0.
391qsnowdiag(:,:)= 0.
392dqreva(:,:)   = 0.
393dqrauto(:,:)  = 0.
394dqrmelt(:,:)  = 0.
395dqrfreez(:,:) = 0.
396dqrcol(:,:)   = 0.
397dqssub(:,:)   = 0.
398dqsauto(:,:)  = 0.
399dqsrim(:,:)   = 0.
400dqsagg(:,:)   = 0.
401dqsfreez(:,:) = 0.
402dqsmelt(:,:)  = 0.
403zqupnew(:)    = 0.
404zqvapclr(:)   = 0.
405
406
407
408!c_iso: variable initialisation for iso
409
410
411!===============================================================================
412!              BEGINNING OF VERTICAL LOOP FROM TOP TO BOTTOM
413!===============================================================================
414
415ncoreczq=0
416
417DO k = klev, 1, -1
418
419    IF (k<=klev-1) THEN
420       iftop=.FALSE.
421    ELSE
422       iftop=.TRUE.
423    ENDIF
424
425    ! Initialisation temperature and specific humidity
426    ! temp(klon,klev) is not modified by the routine, instead all changes in temperature are made on zt
427    ! at the end of the klon loop, a temperature incremtent d_t due to all processes
428    ! (thermalization, evap/sub incoming precip, cloud formation, precipitation processes) is calculated
429    ! d_t = temperature tendency due to lscp
430    ! The temperature of the overlying layer is updated here because needed for thermalization
431    DO i = 1, klon
432        zt(i)=temp(i,k)
433        zq(i)=qt(i,k)
434        IF (.not. iftop) THEN
435           ztupnew(i)  = temp(i,k+1) + d_t(i,k+1)
436           zqupnew(i)  = qt(i,k+1) + d_q(i,k+1) + d_ql(i,k+1) + d_qi(i,k+1)
437           !--zqs(i) is the saturation specific humidity in the layer above
438           zqvapclr(i) = MAX(0., qt(i,k+1) + d_q(i,k+1) - rneb(i,k+1) * zqs(i))
439        ENDIF
440        !c_iso init of iso
441    ENDDO
442   
443    !================================================================
444    ! Flag for the new and more microphysical treatment of precipitation from Atelier Nuage (R)
445    IF (ok_poprecip) THEN
446
447            CALL poprecip_precld(klon, dtime, iftop, paprs(:,k), paprs(:,k+1), pplay(:,k), &
448                              zt, ztupnew, zq, zmqc, znebprecipclr, znebprecipcld, &
449                              zqvapclr, zqupnew, &
450                              zrfl, zrflclr, zrflcld, &
451                              zifl, ziflclr, ziflcld, &
452                              dqreva(:,k),dqssub(:,k) &
453                             )
454
455    !================================================================
456    ELSE
457
458    ! --------------------------------------------------------------------
459    ! P1> Thermalization of precipitation falling from the overlying layer
460    ! --------------------------------------------------------------------
461    ! Computes air temperature variation due to enthalpy transported by
462    ! precipitation. Precipitation is then thermalized with the air in the
463    ! layer.
464    ! The precipitation should remain thermalized throughout the different
465    ! thermodynamical transformations.
466    ! The corresponding water mass should
467    ! be added when calculating the layer's enthalpy change with
468    ! temperature
469    ! See lmdzpedia page todoan
470    ! todoan: check consistency with ice phase
471    ! todoan: understand why several steps
472    ! ---------------------------------------------------------------------
473   
474    IF (iftop) THEN
475
476        DO i = 1, klon
477            zmqc(i) = 0.
478        ENDDO
479
480    ELSE
481
482        DO i = 1, klon
483 
484            zmair=(paprs(i,k)-paprs(i,k+1))/RG
485            ! no condensed water so cp=cp(vapor+dry air)
486            ! RVTMP2=rcpv/rcpd-1
487            zcpair=RCPD*(1.0+RVTMP2*zq(i))
488            zcpeau=RCPD*RVTMP2
489
490            ! zmqc: precipitation mass that has to be thermalized with
491            ! layer's air so that precipitation at the ground has the
492            ! same temperature as the lowermost layer
493            zmqc(i) = (zrfl(i)+zifl(i))*dtime/zmair
494            ! t(i,k+1)+d_t(i,k+1): new temperature of the overlying layer
495            zt(i) = ( ztupnew(i)*zmqc(i)*zcpeau + zcpair*zt(i) ) &
496             / (zcpair + zmqc(i)*zcpeau)
497
498        ENDDO
499 
500    ENDIF
501
502    ! --------------------------------------------------------------------
503    ! P2> Precipitation evaporation/sublimation/melting
504    ! --------------------------------------------------------------------
505    ! A part of the precipitation coming from above is evaporated/sublimated/melted.
506    ! --------------------------------------------------------------------
507   
508    IF (iflag_evap_prec>=1) THEN
509
510        ! Calculation of saturation specific humidity
511        ! depending on temperature:
512        CALL calc_qsat_ecmwf(klon,zt(:),qzero(:),pplay(:,k),RTT,0,.FALSE.,zqs(:),zdqs(:))
513        ! wrt liquid water
514        CALL calc_qsat_ecmwf(klon,zt(:),qzero(:),pplay(:,k),RTT,1,.FALSE.,qsl(:),dqsl(:))
515        ! wrt ice
516        CALL calc_qsat_ecmwf(klon,zt(:),qzero(:),pplay(:,k),RTT,2,.FALSE.,qsi(:),dqsi(:))
517
518        DO i = 1, klon
519
520            ! if precipitation
521            IF (zrfl(i)+zifl(i)>0.) THEN
522
523            ! LudoTP: we only account for precipitation evaporation in the clear-sky (iflag_evap_prec>=4).
524            ! c_iso: likely important to distinguish cs from neb isotope precipitation
525
526                IF (iflag_evap_prec>=4) THEN
527                    zrfl(i) = zrflclr(i)
528                    zifl(i) = ziflclr(i)
529                ENDIF
530
531                IF (iflag_evap_prec==1) THEN
532                    znebprecip(i)=zneb(i)
533                ELSE
534                    znebprecip(i)=MAX(zneb(i),znebprecip(i))
535                ENDIF
536
537                IF (iflag_evap_prec>4) THEN
538                ! Max evaporation not to saturate the clear sky precip fraction
539                ! i.e. the fraction where evaporation occurs
540                    zqev0 = MAX(0.0, (zqs(i)-zq(i))*znebprecipclr(i))
541                ELSEIF (iflag_evap_prec == 4) THEN
542                ! Max evaporation not to saturate the whole mesh
543                ! Pay attention -> lead to unrealistic and excessive evaporation
544                    zqev0 = MAX(0.0, zqs(i)-zq(i))
545                ELSE
546                ! Max evap not to saturate the fraction below the cloud
547                    zqev0 = MAX(0.0, (zqs(i)-zq(i))*znebprecip(i))
548                ENDIF
549
550                ! Evaporation of liquid precipitation coming from above
551                ! dP/dz=beta*(1-q/qsat)*sqrt(P)
552                ! formula from Sundquist 1988, Klemp & Wilhemson 1978
553                ! LTP: evaporation only in the clear sky part (iflag_evap_prec>=4)
554
555                IF (iflag_evap_prec==3) THEN
556                    zqevt = znebprecip(i)*coef_eva*(1.0-zq(i)/qsl(i)) &
557                    *SQRT(zrfl(i)/max(1.e-4,znebprecip(i)))        &
558                    *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
559                ELSE IF (iflag_evap_prec>=4) THEN
560                     zqevt = znebprecipclr(i)*coef_eva*(1.0-zq(i)/qsl(i)) &
561                    *SQRT(zrfl(i)/max(1.e-8,znebprecipclr(i))) &
562                    *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
563                ELSE
564                    zqevt = 1.*coef_eva*(1.0-zq(i)/qsl(i))*SQRT(zrfl(i)) &
565                    *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
566                ENDIF
567
568                zqevt = MAX(0.0,MIN(zqevt,zrfl(i))) &
569                    *RG*dtime/(paprs(i,k)-paprs(i,k+1))
570
571                ! sublimation of the solid precipitation coming from above
572                IF (iflag_evap_prec==3) THEN
573                    zqevti = znebprecip(i)*coef_sub*(1.0-zq(i)/qsi(i)) &
574                    *SQRT(zifl(i)/max(1.e-4,znebprecip(i))) &
575                    *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
576                ELSE IF (iflag_evap_prec>=4) THEN
577                     zqevti = znebprecipclr(i)*coef_sub*(1.0-zq(i)/qsi(i)) &
578                    *SQRT(zifl(i)/max(1.e-8,znebprecipclr(i))) &
579                    *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
580                ELSE
581                    zqevti = 1.*coef_sub*(1.0-zq(i)/qsi(i))*SQRT(zifl(i)) &
582                    *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
583                ENDIF
584
585                zqevti = MAX(0.0,MIN(zqevti,zifl(i))) &
586                *RG*dtime/(paprs(i,k)-paprs(i,k+1))
587
588                ! A. JAM
589                ! Evaporation limit: we ensure that the layer's fraction below
590                ! the cloud or the whole mesh (depending on iflag_evap_prec)
591                ! does not reach saturation. In this case, we
592                ! redistribute zqev0 conserving the ratio liquid/ice
593
594                IF (zqevt+zqevti>zqev0) THEN
595                    zqev=zqev0*zqevt/(zqevt+zqevti)
596                    zqevi=zqev0*zqevti/(zqevt+zqevti)
597                ELSE
598                    zqev=zqevt
599                    zqevi=zqevti
600                ENDIF
601
602
603                ! New solid and liquid precipitation fluxes after evap and sublimation
604                zrfln(i) = Max(0.,zrfl(i) - zqev*(paprs(i,k)-paprs(i,k+1))  &
605                         /RG/dtime)
606                zifln(i) = Max(0.,zifl(i) - zqevi*(paprs(i,k)-paprs(i,k+1)) &
607                         /RG/dtime)
608
609
610                ! vapor, temperature, precip fluxes update
611                ! vapor is updated after evaporation/sublimation (it is increased)
612                zq(i) = zq(i) - (zrfln(i)+zifln(i)-zrfl(i)-zifl(i)) &
613                * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
614                ! zmqc is the total condensed water in the precip flux (it is decreased)
615                zmqc(i) = zmqc(i) + (zrfln(i)+zifln(i)-zrfl(i)-zifl(i)) &
616                * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
617                ! air and precip temperature (i.e., gridbox temperature)
618                ! is updated due to latent heat cooling
619                zt(i) = zt(i) + (zrfln(i)-zrfl(i))        &
620                * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime    &
621                * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i))) &
622                + (zifln(i)-zifl(i))                      &
623                * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime    &
624                * RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))
625
626                ! New values of liquid and solid precipitation
627                zrfl(i) = zrfln(i)
628                zifl(i) = zifln(i)
629
630                ! c_iso here call_reevap that updates isotopic zrfl, zifl (in inout)
631                !                         due to evap + sublim                                     
632                                           
633
634                IF (iflag_evap_prec>=4) THEN
635                    zrflclr(i) = zrfl(i)
636                    ziflclr(i) = zifl(i)
637                    IF(zrflclr(i) + ziflclr(i)<=0) THEN
638                        znebprecipclr(i) = 0.0
639                    ENDIF
640                    zrfl(i) = zrflclr(i) + zrflcld(i)
641                    zifl(i) = ziflclr(i) + ziflcld(i)
642                ENDIF
643
644                ! c_iso duplicate for isotopes or loop on isotopes
645
646                ! Melting:
647                zmelt = ((zt(i)-RTT)/(ztfondue-RTT)) ! JYG
648                ! precip fraction that is melted
649                zmelt = MIN(MAX(zmelt,0.),1.)
650
651                ! update of rainfall and snowfall due to melting
652                IF (iflag_evap_prec>=4) THEN
653                    zrflclr(i)=zrflclr(i)+zmelt*ziflclr(i)
654                    zrflcld(i)=zrflcld(i)+zmelt*ziflcld(i)
655                    zrfl(i)=zrflclr(i)+zrflcld(i)
656                ELSE
657                    zrfl(i)=zrfl(i)+zmelt*zifl(i)
658                ENDIF
659
660
661                ! c_iso: melting of isotopic precipi with zmelt (no fractionation)
662
663                ! Latent heat of melting because of precipitation melting
664                ! NB: the air + precip temperature is simultaneously updated
665                zt(i)=zt(i)-zifl(i)*zmelt*(RG*dtime)/(paprs(i,k)-paprs(i,k+1)) &
666                *RLMLT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))
667               
668                IF (iflag_evap_prec>=4) THEN
669                    ziflclr(i)=ziflclr(i)*(1.-zmelt)
670                    ziflcld(i)=ziflcld(i)*(1.-zmelt)
671                    zifl(i)=ziflclr(i)+ziflcld(i)
672                ELSE
673                    zifl(i)=zifl(i)*(1.-zmelt)
674                ENDIF
675
676            ELSE
677                ! if no precip, we reinitialize the cloud fraction used for the precip to 0
678                znebprecip(i)=0.
679
680            ENDIF ! (zrfl(i)+zifl(i).GT.0.)
681
682        ENDDO ! loop on klon
683
684    ENDIF ! (iflag_evap_prec>=1)
685
686    ENDIF ! (ok_poprecip)
687
688    ! --------------------------------------------------------------------
689    ! End precip evaporation
690    ! --------------------------------------------------------------------
691   
692    ! Calculation of qsat, L/Cp*dqsat/dT and ncoreczq counter
693    !-------------------------------------------------------
694
695     qtot(:)=zq(:)+zmqc(:)
696     CALL calc_qsat_ecmwf(klon,zt(:),qtot(:),pplay(:,k),RTT,0,.FALSE.,zqs(:),zdqs(:))
697     DO i = 1, klon
698            zdelta = MAX(0.,SIGN(1.,RTT-zt(i)))
699            zdqsdT_raw(i) = zdqs(i)*RCPD*(1.0+RVTMP2*zq(i)) / (RLVTT*(1.-zdelta) + RLSTT*zdelta)
700            IF (zq(i) < 1.e-15) THEN
701                ncoreczq=ncoreczq+1
702                zq(i)=1.e-15
703            ENDIF
704            ! c_iso: do something similar for isotopes
705
706     ENDDO
707   
708    ! --------------------------------------------------------------------
709    ! P2> Cloud formation
710    !---------------------------------------------------------------------
711
712    ! Unlike fisrtilp, we always assume a 'fractional cloud' approach
713    ! i.e. clouds occupy only a fraction of the mesh (the subgrid distribution
714    ! is prescribed and depends on large scale variables and boundary layer
715    ! properties)
716    ! The decrease in condensed part due tu latent heating is taken into
717    ! account
718    ! -------------------------------------------------------------------
719
720        ! P2.1> With the PDFs (log-normal, bigaussian)
721        ! cloud properties calculation with the initial values of t and q
722        ! ----------------------------------------------------------------
723
724        ! initialise gammasat and qincloud_mpc
725        gammasat(:)=1.
726        qincloud_mpc(:)=0.
727
728        IF (iflag_cld_th>=5) THEN
729               ! Cloud cover and content in meshes affected by shallow convection,
730               ! are retrieved from a bi-gaussian distribution of the saturation deficit
731               ! following Jam et al. 2013
732
733               IF (iflag_cloudth_vert<=2) THEN
734                  ! Old version of Arnaud Jam
735
736                    CALL cloudth(klon,klev,k,tv,                  &
737                         zq,qta,fraca,                            &
738                         qcloud,ctot,pspsk,paprs,pplay,tla,thl, &
739                         ratqs,zqs,temp,                              &
740                         cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
741
742
743                ELSEIF (iflag_cloudth_vert>=3 .AND. iflag_cloudth_vert<=5) THEN
744                   ! Default version of Arnaud Jam
745                       
746                    CALL cloudth_v3(klon,klev,k,tv,                        &
747                         zq,qta,fraca,                                     &
748                         qcloud,ctot,ctot_vol,pspsk,paprs,pplay,tla,thl, &
749                         ratqs,zqs,temp, &
750                         cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
751
752
753                ELSEIF (iflag_cloudth_vert==6) THEN
754                   ! Jean Jouhaud's version, with specific separation between surface and volume
755                   ! cloud fraction Decembre 2018
756
757                    CALL cloudth_v6(klon,klev,k,tv,                        &
758                         zq,qta,fraca,                                     &
759                         qcloud,ctot,ctot_vol,pspsk,paprs,pplay,tla,thl, &
760                         ratqs,zqs,temp, &
761                         cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
762
763                ELSEIF (iflag_cloudth_vert == 7) THEN
764                   ! Updated version of Arnaud Jam (correction by E. Vignon) + adapted treatment
765                   ! for boundary-layer mixed phase clouds
766                    CALL cloudth_mpc(klon,klev,k,mpc_bl_points,zt,zq,qta(:,k),fraca(:,k), &
767                                     pspsk(:,k),paprs(:,k+1),paprs(:,k),pplay(:,k), tla(:,k), &
768                                     ratqs(:,k),qcloud,qincloud_mpc,zfice_th,ctot(:,k),ctot_vol(:,k), &
769                                     cloudth_sth(:,k),cloudth_senv(:,k),cloudth_sigmath(:,k),cloudth_sigmaenv(:,k))
770
771                ENDIF
772
773
774                DO i=1,klon
775                    rneb(i,k)=ctot(i,k)
776                    rneblsvol(i,k)=ctot_vol(i,k)
777                    zqn(i)=qcloud(i)
778                ENDDO
779
780        ENDIF
781
782        IF (iflag_cld_th <= 4) THEN
783           
784                ! lognormal
785            lognormale(:) = .TRUE.
786
787        ELSEIF (iflag_cld_th >= 6) THEN
788           
789                ! lognormal distribution when no thermals
790            lognormale(:) = fraca(:,k) < min_frac_th_cld
791
792        ELSE
793                ! When iflag_cld_th=5, we always assume
794                ! bi-gaussian distribution
795            lognormale(:) = .FALSE.
796       
797        ENDIF
798
799        DT(:) = 0.
800        n_i(:)=0
801        Tbef(:)=zt(:)
802        qlbef(:)=0.
803
804        ! Treatment of non-boundary layer clouds (lognormale)
805        ! condensation with qsat(T) variation (adaptation)
806        ! Iterative resolution to converge towards qsat
807        ! with update of temperature, ice fraction and qsat at
808        ! each iteration
809
810        ! todoan -> sensitivity to iflag_fisrtilp_qsat
811        DO iter=1,iflag_fisrtilp_qsat+1
812
813                DO i=1,klon
814
815                    ! keepgoing = .TRUE. while convergence is not satisfied
816                    keepgoing(i)=ABS(DT(i))>DDT0
817
818                    IF ((keepgoing(i) .OR. (n_i(i) == 0)) .AND. lognormale(i)) THEN
819
820                        ! if not convergence:
821
822                        ! P2.2.1> cloud fraction and condensed water mass calculation
823                        ! Calculated variables:
824                        ! rneb : cloud fraction
825                        ! zqn : total water within the cloud
826                        ! zcond: mean condensed water within the mesh
827                        ! rhcl: clear-sky relative humidity
828                        !---------------------------------------------------------------
829
830                        ! new temperature that only serves in the iteration process:
831                        Tbef(i)=Tbef(i)+DT(i)
832
833                        ! Rneb, qzn and zcond for lognormal PDFs
834                        qtot(i)=zq(i)+zmqc(i)
835
836                      ENDIF
837
838                  ENDDO
839
840                  ! Calculation of saturation specific humidity and ice fraction
841                  CALL calc_qsat_ecmwf(klon,Tbef(:),qtot(:),pplay(:,k),RTT,0,.FALSE.,zqs(:),zdqs(:))
842                  CALL calc_gammasat(klon,Tbef(:),qtot(:),pplay(:,k),gammasat(:),dgammasatdt(:))
843                  ! saturation may occur at a humidity different from qsat (gamma qsat), so gamma correction for dqs
844                  zdqs(:) = gammasat(:)*zdqs(:)+zqs(:)*dgammasatdt(:)
845                  ! cloud phase determination
846                  IF (iflag_t_glace>=4) THEN
847                  ! For iflag_t_glace GE 4 the phase partition function dependends on temperature AND distance to cloud top
848                       CALL distance_to_cloud_top(klon,klev,k,temp,pplay,paprs,rneb,zdistcltop,ztemp_cltop)
849                  ENDIF
850
851                  CALL icefrac_lscp(klon, zt(:), iflag_ice_thermo, zdistcltop(:),ztemp_cltop(:),zfice(:),dzfice(:))
852
853                  DO i=1,klon !todoan : check if loop in i is needed
854
855                      IF ((keepgoing(i) .OR. (n_i(i) == 0)) .AND. lognormale(i)) THEN
856
857                        zpdf_sig(i)=ratqs(i,k)*zq(i)
858                        zpdf_k(i)=-sqrt(log(1.+(zpdf_sig(i)/zq(i))**2))
859                        zpdf_delta(i)=log(zq(i)/(gammasat(i)*zqs(i)))
860                        zpdf_a(i)=zpdf_delta(i)/(zpdf_k(i)*sqrt(2.))
861                        zpdf_b(i)=zpdf_k(i)/(2.*sqrt(2.))
862                        zpdf_e1(i)=zpdf_a(i)-zpdf_b(i)
863                        zpdf_e1(i)=sign(min(ABS(zpdf_e1(i)),5.),zpdf_e1(i))
864                        zpdf_e1(i)=1.-erf(zpdf_e1(i))
865                        zpdf_e2(i)=zpdf_a(i)+zpdf_b(i)
866                        zpdf_e2(i)=sign(min(ABS(zpdf_e2(i)),5.),zpdf_e2(i))
867                        zpdf_e2(i)=1.-erf(zpdf_e2(i))
868
869                        IF ((.NOT.ok_ice_sursat).OR.(Tbef(i)>t_glace_min)) THEN
870
871                          IF (zpdf_e1(i)<1.e-10) THEN
872                              rneb(i,k)=0.
873                              zqn(i)=gammasat(i)*zqs(i)
874                          ELSE
875                              rneb(i,k)=0.5*zpdf_e1(i)
876                              zqn(i)=zq(i)*zpdf_e2(i)/zpdf_e1(i)
877                          ENDIF
878
879                          rnebss(i,k)=0.0   !--added by OB (needed because of the convergence loop on time)
880                          fcontrN(i,k)=0.0  !--idem
881                          fcontrP(i,k)=0.0  !--idem
882                          qss(i,k)=0.0      !--idem
883
884                       ELSE ! in case of ice supersaturation by Audran
885
886                        !------------------------------------
887                        ! ICE SUPERSATURATION
888                        !------------------------------------
889
890                        CALL ice_sursat(pplay(i,k), paprs(i,k)-paprs(i,k+1), dtime, i, k, temp(i,k), zq(i), &
891                             gamma_ss(i,k), zqs(i), Tbef(i), rneb_seri(i,k), ratqs(i,k),                 &
892                             rneb(i,k), zqn(i), rnebss(i,k), qss(i,k),                                   &
893                             Tcontr(i,k), qcontr(i,k), qcontr2(i,k), fcontrN(i,k), fcontrP(i,k) )
894
895                        ENDIF ! ((flag_ice_sursat.eq.0).or.(Tbef(i).gt.t_glace_min))
896
897                        ! If vertical heterogeneity, change fraction by volume as well
898                        IF (iflag_cloudth_vert>=3) THEN
899                            ctot_vol(i,k)=rneb(i,k)
900                            rneblsvol(i,k)=ctot_vol(i,k)
901                        ENDIF
902
903
904                       ! P2.2.2> Approximative calculation of temperature variation DT
905                       !           due to condensation.
906                       ! Calculated variables:
907                       ! dT : temperature change due to condensation
908                       !---------------------------------------------------------------
909
910               
911                        IF (zfice(i)<1) THEN
912                            cste=RLVTT
913                        ELSE
914                            cste=RLSTT
915                        ENDIF
916                       
917                        ! LEA_R : check formule
918                        qlbef(i)=max(0.,zqn(i)-zqs(i))
919                        num = -Tbef(i)+zt(i)+rneb(i,k)*((1-zfice(i))*RLVTT &
920                              +zfice(i)*RLSTT)/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))*qlbef(i)
921                        denom = 1.+rneb(i,k)*((1-zfice(i))*RLVTT+zfice(i)*RLSTT)/cste*zdqs(i) &
922                              -(RLSTT-RLVTT)/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))*rneb(i,k)      &
923                              *qlbef(i)*dzfice(i)
924                        ! here we update a provisory temperature variable that only serves in the iteration
925                        ! process
926                        DT(i)=num/denom
927                        n_i(i)=n_i(i)+1
928
929                    ENDIF ! end keepgoing
930
931                ENDDO     ! end loop on i
932
933        ENDDO         ! iter=1,iflag_fisrtilp_qsat+1
934
935        ! P2.3> Final quantities calculation
936        ! Calculated variables:
937        !   rneb : cloud fraction
938        !   zcond: mean condensed water in the mesh
939        !   zqn  : mean water vapor in the mesh
940        !   zfice: ice fraction in clouds
941        !   zt   : temperature
942        !   rhcl : clear-sky relative humidity
943        ! ----------------------------------------------------------------
944
945
946        ! For iflag_t_glace GE 4 the phase partition function dependends on temperature AND distance to cloud top
947        IF (iflag_t_glace>=4) THEN
948           CALL distance_to_cloud_top(klon,klev,k,temp,pplay,paprs,rneb,zdistcltop,ztemp_cltop)
949           distcltop(:,k)=zdistcltop(:)
950           temp_cltop(:,k)=ztemp_cltop(:)
951        ENDIF
952
953        ! Partition function depending on temperature
954        CALL icefrac_lscp(klon, Tbef, iflag_ice_thermo, zdistcltop, ztemp_cltop, zfice, dzfice)
955
956        ! Partition function depending on tke for non shallow-convective clouds
957        IF (iflag_icefrac >= 1) THEN
958
959           CALL icefrac_lscp_turb(klon, dtime, Tbef, pplay(:,k), paprs(:,k), paprs(:,k+1), qice_save(:,k), ziflcld, zqn, &
960           rneb(:,k), tke(:,k), tke_dissip(:,k), zqliq, zqvapcl, zqice, zfice_turb, dzfice_turb, cldfraliq(:,k),sigma2_icefracturb(:,k), mean_icefracturb(:,k))
961
962        ENDIF
963
964        ! Water vapor update, Phase determination and subsequent latent heat exchange
965        DO i=1, klon
966            ! Overwrite phase partitioning in boundary layer mixed phase clouds when the
967            ! iflag_cloudth_vert=7 and specific param is activated
968            IF (mpc_bl_points(i,k) > 0) THEN
969                zcond(i) = MAX(0.0,qincloud_mpc(i))*rneb(i,k)
970                ! following line is very strange and probably wrong
971                rhcl(i,k)= (zqs(i)+zq(i))/2./zqs(i)
972                ! water vapor update and partition function if thermals
973                zq(i) = zq(i) - zcond(i)       
974                zfice(i)=zfice_th(i)
975            ELSE
976                ! Checks on rneb, rhcl and zqn
977                IF (rneb(i,k) <= 0.0) THEN
978                    zqn(i) = 0.0
979                    rneb(i,k) = 0.0
980                    zcond(i) = 0.0
981                    rhcl(i,k)=zq(i)/zqs(i)
982                ELSE IF (rneb(i,k) >= 1.0) THEN
983                    zqn(i) = zq(i)
984                    rneb(i,k) = 1.0
985                    zcond(i) = MAX(0.0,zqn(i)-zqs(i))
986                    rhcl(i,k)=1.0
987                ELSE
988                    zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)
989                    ! following line is very strange and probably wrong:
990                    rhcl(i,k)=(zqs(i)+zq(i))/2./zqs(i)
991                    ! Overwrite partitioning for non shallow-convective clouds if iflag_icefrac>1 (icefrac turb param)
992                    IF (iflag_icefrac >= 1) THEN
993                        IF (lognormale(i)) THEN 
994                           zcond(i)  = zqliq(i) + zqice(i)
995                           zfice(i)=zfice_turb(i)
996                           rhcl(i,k) = zqvapcl(i) * rneb(i,k) + (zq(i) - zqn(i)) * (1.-rneb(i,k))
997                        ENDIF
998                    ENDIF
999                ENDIF
1000
1001                ! water vapor update
1002                zq(i) = zq(i) - zcond(i)
1003
1004            ENDIF
1005
1006            ! c_iso : routine that computes in-cloud supersaturation
1007            ! c_iso condensation of isotopes (zcond, zsursat, zfice, zq in input)
1008
1009            ! temperature update due to phase change
1010            zt(i) = zt(i) + (1.-zfice(i))*zcond(i) &
1011       * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i))) &
1012                  +zfice(i)*zcond(i) * RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i)))
1013        ENDDO
1014
1015        ! If vertical heterogeneity, change volume fraction
1016        IF (iflag_cloudth_vert >= 3) THEN
1017          ctot_vol(1:klon,k)=min(max(ctot_vol(1:klon,k),0.),1.)
1018          rneblsvol(1:klon,k)=ctot_vol(1:klon,k)
1019        ENDIF
1020
1021    ! ----------------------------------------------------------------
1022    ! P3> Precipitation formation
1023    ! ----------------------------------------------------------------
1024
1025    !================================================================
1026    IF (ok_poprecip) THEN
1027
1028      DO i = 1, klon
1029        zoliql(i) = zcond(i) * ( 1. - zfice(i) )
1030        zoliqi(i) = zcond(i) * zfice(i)
1031      ENDDO
1032
1033      CALL poprecip_postcld(klon, dtime, paprs(:,k), paprs(:,k+1), pplay(:,k), &
1034                            ctot_vol(:,k), ptconv(:,k), &
1035                            zt, zq, zoliql, zoliqi, zfice, &
1036                            rneb(:,k), znebprecipclr, znebprecipcld, &
1037                            zrfl, zrflclr, zrflcld, &
1038                            zifl, ziflclr, ziflcld, &
1039                            qraindiag(:,k), qsnowdiag(:,k), dqrauto(:,k), &
1040                            dqrcol(:,k), dqrmelt(:,k), dqrfreez(:,k), &
1041                            dqsauto(:,k), dqsagg(:,k), dqsrim(:,k), &
1042                            dqsmelt(:,k), dqsfreez(:,k) &
1043                            )
1044
1045      DO i = 1, klon
1046        zoliq(i) = zoliql(i) + zoliqi(i)
1047        IF ( zoliq(i) > 0. ) THEN
1048                zfice(i) = zoliqi(i)/zoliq(i)
1049        ELSE 
1050                zfice(i) = 0.0
1051        ENDIF
1052
1053        ! calculation of specific content of condensates seen by radiative scheme
1054        IF (ok_radocond_snow) THEN
1055           radocond(i,k) = zoliq(i)
1056           radocondl(i,k)= radocond(i,k)*(1.-zfice(i))
1057           radocondi(i,k)= radocond(i,k)*zfice(i)+qsnowdiag(i,k)
1058        ELSE
1059           radocond(i,k) = zoliq(i)
1060           radocondl(i,k)= radocond(i,k)*(1.-zfice(i))
1061           radocondi(i,k)= radocond(i,k)*zfice(i)
1062        ENDIF
1063      ENDDO
1064
1065    !================================================================
1066    ELSE
1067
1068    ! LTP:
1069    IF (iflag_evap_prec >= 4) THEN
1070
1071        !Partitionning between precipitation coming from clouds and that coming from CS
1072
1073        !0) Calculate tot_zneb, total cloud fraction above the cloud
1074        !assuming a maximum-random overlap (voir Jakob and Klein, 2000)
1075       
1076        DO i=1, klon
1077                tot_znebn(i) = 1. - (1.-tot_zneb(i))*(1 - max(rneb(i,k),zneb(i))) &
1078                        /(1.-min(zneb(i),1.-smallestreal))
1079                d_tot_zneb(i) = tot_znebn(i) - tot_zneb(i)
1080                tot_zneb(i) = tot_znebn(i)
1081
1082
1083                !1) Cloudy to clear air
1084                d_znebprecip_cld_clr(i) = znebprecipcld(i) - min(rneb(i,k),znebprecipcld(i))
1085                IF (znebprecipcld(i) > 0.) THEN
1086                        d_zrfl_cld_clr(i) = d_znebprecip_cld_clr(i)/znebprecipcld(i)*zrflcld(i)
1087                        d_zifl_cld_clr(i) = d_znebprecip_cld_clr(i)/znebprecipcld(i)*ziflcld(i)
1088                ELSE
1089                        d_zrfl_cld_clr(i) = 0.
1090                        d_zifl_cld_clr(i) = 0.
1091                ENDIF
1092
1093                !2) Clear to cloudy air
1094                d_znebprecip_clr_cld(i) = max(0., min(znebprecipclr(i), rneb(i,k) - d_tot_zneb(i) - zneb(i)))
1095                IF (znebprecipclr(i) > 0) THEN
1096                        d_zrfl_clr_cld(i) = d_znebprecip_clr_cld(i)/znebprecipclr(i)*zrflclr(i)
1097                        d_zifl_clr_cld(i) = d_znebprecip_clr_cld(i)/znebprecipclr(i)*ziflclr(i)
1098                ELSE
1099                        d_zrfl_clr_cld(i) = 0.
1100                        d_zifl_clr_cld(i) = 0.
1101                ENDIF
1102
1103                !Update variables
1104                znebprecipcld(i) = znebprecipcld(i) + d_znebprecip_clr_cld(i) - d_znebprecip_cld_clr(i)
1105                znebprecipclr(i) = znebprecipclr(i) + d_znebprecip_cld_clr(i) - d_znebprecip_clr_cld(i)
1106                zrflcld(i) = zrflcld(i) + d_zrfl_clr_cld(i) - d_zrfl_cld_clr(i)
1107                ziflcld(i) = ziflcld(i) + d_zifl_clr_cld(i) - d_zifl_cld_clr(i)
1108                zrflclr(i) = zrflclr(i) + d_zrfl_cld_clr(i) - d_zrfl_clr_cld(i)
1109                ziflclr(i) = ziflclr(i) + d_zifl_cld_clr(i) - d_zifl_clr_cld(i)
1110
1111                ! c_iso  do the same thing for isotopes precip
1112        ENDDO
1113    ENDIF
1114
1115
1116    ! Autoconversion
1117    ! -------------------------------------------------------------------------------
1118
1119
1120    ! Initialisation of zoliq and radocond variables
1121
1122    DO i = 1, klon
1123            zoliq(i) = zcond(i)
1124            zoliqi(i)= zoliq(i)*zfice(i)
1125            zoliql(i)= zoliq(i)*(1.-zfice(i))
1126            ! c_iso : initialisation of zoliq* also for isotopes
1127            zrho(i,k)  = pplay(i,k) / zt(i) / RD
1128            zdz(i)   = (paprs(i,k)-paprs(i,k+1)) / (zrho(i,k)*RG)
1129            iwc(i)   = 0.
1130            zneb(i)  = MAX(rneb(i,k),seuil_neb)
1131            radocond(i,k) = zoliq(i)/REAL(niter_lscp+1)
1132            radocondi(i,k)=zoliq(i)*zfice(i)/REAL(niter_lscp+1)
1133            radocondl(i,k)=zoliq(i)*(1.-zfice(i))/REAL(niter_lscp+1)
1134    ENDDO
1135
1136       
1137    DO n = 1, niter_lscp
1138
1139        DO i=1,klon
1140            IF (rneb(i,k)>0.0) THEN
1141                iwc(i) = zrho(i,k) * zoliqi(i) / zneb(i) ! in-cloud ice condensate content
1142            ENDIF
1143        ENDDO
1144
1145        CALL fallice_velocity(klon,iwc(:),zt(:),zrho(:,k),pplay(:,k),ptconv(:,k),velo(:,k))
1146
1147        DO i = 1, klon
1148
1149            IF (rneb(i,k)>0.0) THEN
1150
1151                ! Initialization of zrain, zsnow and zprecip:
1152                zrain=0.
1153                zsnow=0.
1154                zprecip=0.
1155                ! c_iso same init for isotopes. Externalisation?
1156
1157                IF (zneb(i)==seuil_neb) THEN
1158                    zprecip = 0.0
1159                    zsnow = 0.0
1160                    zrain= 0.0
1161                ELSE
1162
1163                    IF (ptconv(i,k)) THEN ! if convective point
1164                        zcl=cld_lc_con
1165                        zct=1./cld_tau_con
1166                        zexpo=cld_expo_con
1167                        ffallv=ffallv_con
1168                    ELSE
1169                        zcl=cld_lc_lsc
1170                        zct=1./cld_tau_lsc
1171                        zexpo=cld_expo_lsc
1172                        ffallv=ffallv_lsc
1173                    ENDIF
1174
1175
1176                    ! if vertical heterogeneity is taken into account, we use
1177                    ! the "true" volume fraction instead of a modified
1178                    ! surface fraction (which is larger and artificially
1179                    ! reduces the in-cloud water).
1180
1181                    ! Liquid water quantity to remove: zchau (Sundqvist, 1978)
1182                    ! dqliq/dt=-qliq/tau*(1-exp(-qcin/clw)**2)
1183                    !.........................................................
1184                    IF ((iflag_cloudth_vert>=3).AND.(iflag_rain_incloud_vol==1)) THEN
1185
1186                    ! if vertical heterogeneity is taken into account, we use
1187                    ! the "true" volume fraction instead of a modified
1188                    ! surface fraction (which is larger and artificially
1189                    ! reduces the in-cloud water).
1190                        effective_zneb=ctot_vol(i,k)
1191                    ELSE
1192                        effective_zneb=zneb(i)
1193                    ENDIF
1194
1195
1196                    IF (iflag_autoconversion == 2) THEN
1197                    ! two-steps resolution with niter_lscp=1 sufficient
1198                    ! we first treat the second term (with exponential) in an explicit way
1199                    ! and then treat the first term (-q/tau) in an exact way
1200                        zchau=zoliql(i)*(1.-exp(-dtime/REAL(niter_lscp)*zct &
1201                        *(1.-exp(-(zoliql(i)/effective_zneb/zcl)**zexpo))))
1202                    ELSE
1203                    ! old explicit resolution with subtimesteps
1204                        zchau = zct*dtime/REAL(niter_lscp)*zoliql(i) &
1205                        *(1.0-EXP(-(zoliql(i)/effective_zneb/zcl)**zexpo))
1206                    ENDIF
1207
1208
1209                    ! Ice water quantity to remove (Zender & Kiehl, 1997)
1210                    ! dqice/dt=1/rho*d(rho*wice*qice)/dz
1211                    !....................................
1212                    IF (iflag_autoconversion == 2) THEN
1213                    ! exact resolution, niter_lscp=1 is sufficient but works only
1214                    ! with iflag_vice=0
1215                       IF (zoliqi(i) > 0.) THEN
1216                          zfroi=(zoliqi(i)-((zoliqi(i)**(-dice_velo)) &
1217                          +dice_velo*dtime/REAL(niter_lscp)*cice_velo/zdz(i)*ffallv)**(-1./dice_velo))
1218                       ELSE
1219                          zfroi=0.
1220                       ENDIF
1221                    ELSE
1222                    ! old explicit resolution with subtimesteps
1223                       zfroi = dtime/REAL(niter_lscp)/zdz(i)*zoliqi(i)*velo(i,k)
1224                    ENDIF
1225
1226                    zrain   = MIN(MAX(zchau,0.0),zoliql(i))
1227                    zsnow   = MIN(MAX(zfroi,0.0),zoliqi(i))
1228                    zprecip = MAX(zrain + zsnow,0.0)
1229
1230                ENDIF
1231
1232                IF (iflag_autoconversion >= 1) THEN
1233                   ! debugged version with phase conservation through the autoconversion process
1234                   zoliql(i) = MAX(zoliql(i)-1.*zrain  , 0.0)
1235                   zoliqi(i) = MAX(zoliqi(i)-1.*zsnow  , 0.0)
1236                   zoliq(i)  = MAX(zoliq(i)-zprecip , 0.0)
1237                ELSE
1238                   ! bugged version with phase resetting
1239                   zoliql(i) = MAX(zoliq(i)*(1.-zfice(i))-1.*zrain  , 0.0)
1240                   zoliqi(i) = MAX(zoliq(i)*zfice(i)-1.*zsnow  , 0.0)
1241                   zoliq(i)  = MAX(zoliq(i)-zprecip , 0.0)
1242                ENDIF
1243
1244                ! c_iso: CALL isotope_conversion (for readibility) or duplicate
1245
1246                radocond(i,k) = radocond(i,k) + zoliq(i)/REAL(niter_lscp+1)
1247                radocondl(i,k) = radocondl(i,k) + zoliql(i)/REAL(niter_lscp+1)
1248                radocondi(i,k) = radocondi(i,k) + zoliqi(i)/REAL(niter_lscp+1)
1249
1250            ENDIF ! rneb >0
1251
1252        ENDDO  ! i = 1,klon
1253
1254    ENDDO      ! n = 1,niter
1255
1256    ! Precipitation flux calculation
1257
1258    DO i = 1, klon
1259           
1260            IF (iflag_evap_prec>=4) THEN
1261                ziflprev(i)=ziflcld(i)
1262            ELSE
1263                ziflprev(i)=zifl(i)*zneb(i)
1264            ENDIF
1265
1266            IF (rneb(i,k) > 0.0) THEN
1267
1268               ! CR&JYG: We account for the Wegener-Findeisen-Bergeron process in the precipitation flux:
1269               ! If T<0C, liquid precip are converted into ice, which leads to an increase in
1270               ! temperature DeltaT. The effect of DeltaT on condensates and precipitation is roughly
1271               ! taken into account through a linearization of the equations and by approximating
1272               ! the liquid precipitation process with a "threshold" process. We assume that
1273               ! condensates are not modified during this operation. Liquid precipitation is
1274               ! removed (in the limit DeltaT<273.15-T). Solid precipitation increases.
1275               ! Water vapor increases as well
1276               ! Note that compared to fisrtilp, we always assume iflag_bergeron=2
1277
1278                    zqpreci(i)=(zcond(i)-zoliq(i))*zfice(i)
1279                    zqprecl(i)=(zcond(i)-zoliq(i))*(1.-zfice(i))
1280                    zcp=RCPD*(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i)))
1281                    coef1 = rneb(i,k)*RLSTT/zcp*zdqsdT_raw(i)
1282                    ! Computation of DT if all the liquid precip freezes
1283                    DeltaT = RLMLT*zqprecl(i) / (zcp*(1.+coef1))
1284                    ! T should not exceed the freezing point
1285                    ! that is Delta > RTT-zt(i)
1286                    DeltaT = max( min( RTT-zt(i), DeltaT) , 0. )
1287                    zt(i)  = zt(i) + DeltaT
1288                    ! water vaporization due to temp. increase
1289                    Deltaq = rneb(i,k)*zdqsdT_raw(i)*DeltaT
1290                    ! we add this vaporized water to the total vapor and we remove it from the precip
1291                    zq(i) = zq(i) +  Deltaq
1292                    ! The three "max" lines herebelow protect from rounding errors
1293                    zcond(i) = max( zcond(i)- Deltaq, 0. )
1294                    ! liquid precipitation converted to ice precip
1295                    Deltaqprecl = -zcp/RLMLT*(1.+coef1)*DeltaT
1296                    zqprecl(i) = max( zqprecl(i) + Deltaqprecl, 0. )
1297                    ! iced water budget
1298                    zqpreci(i) = max (zqpreci(i) - Deltaqprecl - Deltaq, 0.)
1299
1300                    ! c_iso : mv here condensation of isotopes + redispatchage en precip
1301
1302                IF (iflag_evap_prec>=4) THEN
1303                    zrflcld(i) = zrflcld(i)+zqprecl(i) &
1304                    *(paprs(i,k)-paprs(i,k+1))/(RG*dtime)
1305                    ziflcld(i) = ziflcld(i)+ zqpreci(i) &
1306                        *(paprs(i,k)-paprs(i,k+1))/(RG*dtime)
1307                    znebprecipcld(i) = rneb(i,k)
1308                    zrfl(i) = zrflcld(i) + zrflclr(i)
1309                    zifl(i) = ziflcld(i) + ziflclr(i)
1310                ELSE
1311                    zrfl(i) = zrfl(i)+ zqprecl(i)        &
1312                    *(paprs(i,k)-paprs(i,k+1))/(RG*dtime)
1313                    zifl(i) = zifl(i)+ zqpreci(i)        &
1314                   *(paprs(i,k)-paprs(i,k+1))/(RG*dtime)
1315                ENDIF
1316                ! c_iso : same for isotopes
1317 
1318           ENDIF ! rneb>0
1319
1320   ENDDO
1321
1322    ! LTP: limit of surface cloud fraction covered by precipitation when the local intensity of the flux is below rain_int_min
1323    ! if iflag_evap_prec>=4
1324    IF (iflag_evap_prec>=4) THEN
1325
1326        DO i=1,klon
1327
1328            IF ((zrflclr(i) + ziflclr(i)) > 0. ) THEN
1329                znebprecipclr(i) = min(znebprecipclr(i),max(zrflclr(i)/ &
1330                (MAX(znebprecipclr(i),seuil_neb)*rain_int_min), ziflclr(i)/(MAX(znebprecipclr(i),seuil_neb)*rain_int_min)))
1331            ELSE
1332                znebprecipclr(i)=0.0
1333            ENDIF
1334
1335            IF ((zrflcld(i) + ziflcld(i)) > 0.) THEN
1336                znebprecipcld(i) = min(znebprecipcld(i), max(zrflcld(i)/ &
1337                (MAX(znebprecipcld(i),seuil_neb)*rain_int_min), ziflcld(i)/(MAX(znebprecipcld(i),seuil_neb)*rain_int_min)))
1338            ELSE
1339                znebprecipcld(i)=0.0
1340            ENDIF
1341        !IF ( ((1-zfice(i))*zoliq(i) .GT. 0.) .AND. (zt(i) .LE. 233.15) ) THEN
1342        !PRINT*,'WARNING LEA OLIQ A <-40°C '
1343        !PRINT*,'zt,Tbef,oliq,oice,cldfraliq,icefrac,rneb',zt(i),Tbef(i),(1-zfice(i))*zoliq(i),zfice(i)*zoliq(i),cldfraliq(i,k),zfice(i),rneb(i,k)
1344        !ENDIF
1345        ENDDO
1346
1347    ENDIF
1348
1349
1350    ENDIF ! ok_poprecip
1351
1352    ! End of precipitation formation
1353    ! --------------------------------
1354
1355
1356    ! Calculation of cloud condensates amount seen by radiative scheme
1357    !-----------------------------------------------------------------
1358
1359    ! Calculation of the concentration of condensates seen by the radiation scheme
1360    ! for liquid phase, we take radocondl
1361    ! for ice phase, we take radocondi if we neglect snowfall, otherwise (ok_radocond_snow=true)
1362    ! we recalculate radocondi to account for contributions from the precipitation flux
1363    ! TODO: for the moment, we deactivate this option when ok_poprecip=.TRUE.
1364
1365    IF ((ok_radocond_snow) .AND. (k < klev) .AND. (.NOT. ok_poprecip)) THEN
1366        ! for the solid phase (crystals + snowflakes)
1367        ! we recalculate radocondi by summing
1368        ! the ice content calculated in the mesh
1369        ! + the contribution of the non-evaporated snowfall
1370        ! from the overlying layer
1371        DO i=1,klon
1372           IF (ziflprev(i) /= 0.0) THEN
1373              radocondi(i,k)=zoliq(i)*zfice(i)+zqpreci(i)+ziflprev(i)/zrho(i,k+1)/velo(i,k+1)
1374           ELSE
1375              radocondi(i,k)=zoliq(i)*zfice(i)+zqpreci(i)
1376           ENDIF
1377           radocond(i,k)=radocondl(i,k)+radocondi(i,k)
1378        ENDDO
1379    ENDIF
1380
1381    ! caculate the percentage of ice in "radocond" so cloud+precip seen by the radiation scheme
1382    DO i=1,klon
1383        IF (radocond(i,k) > 0.) THEN
1384            radicefrac(i,k)=MIN(MAX(radocondi(i,k)/radocond(i,k),0.),1.)
1385        ENDIF
1386    ENDDO
1387
1388    ! ----------------------------------------------------------------
1389    ! P4> Wet scavenging
1390    ! ----------------------------------------------------------------
1391
1392    !Scavenging through nucleation in the layer
1393
1394    DO i = 1,klon
1395       
1396        IF(zcond(i)>zoliq(i)+1.e-10) THEN
1397            beta(i,k) = (zcond(i)-zoliq(i))/zcond(i)/dtime
1398        ELSE
1399            beta(i,k) = 0.
1400        ENDIF
1401
1402        zprec_cond(i) = MAX(zcond(i)-zoliq(i),0.0)*(paprs(i,k)-paprs(i,k+1))/RG
1403
1404        IF (rneb(i,k)>0.0.AND.zprec_cond(i)>0.) THEN
1405
1406            IF (temp(i,k) >= t_glace_min) THEN
1407                zalpha_tr = a_tr_sca(3)
1408            ELSE
1409                zalpha_tr = a_tr_sca(4)
1410            ENDIF
1411
1412            zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i))
1413            frac_nucl(i,k)= 1.-zneb(i)*zfrac_lessi
1414
1415            ! Nucleation with a factor of -1 instead of -0.5
1416            zfrac_lessi = 1. - EXP(-zprec_cond(i)/zneb(i))
1417
1418        ENDIF
1419
1420    ENDDO
1421
1422    ! Scavenging through impaction in the underlying layer
1423
1424    DO kk = k-1, 1, -1
1425
1426        DO i = 1, klon
1427
1428            IF (rneb(i,k)>0.0.AND.zprec_cond(i)>0.) THEN
1429
1430                IF (temp(i,kk) >= t_glace_min) THEN
1431                    zalpha_tr = a_tr_sca(1)
1432                ELSE
1433                    zalpha_tr = a_tr_sca(2)
1434                ENDIF
1435
1436              zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i))
1437              frac_impa(i,kk)= 1.-zneb(i)*zfrac_lessi
1438
1439            ENDIF
1440
1441        ENDDO
1442
1443    ENDDO
1444
1445    !--save some variables for ice supersaturation
1446
1447    DO i = 1, klon
1448        ! for memory
1449        rneb_seri(i,k) = rneb(i,k)
1450
1451        ! for diagnostics
1452        rnebclr(i,k) = 1.0 - rneb(i,k) - rnebss(i,k)
1453
1454        qvc(i,k) = zqs(i) * rneb(i,k)
1455        qclr(i,k) = MAX(1.e-10,zq(i) - qvc(i,k) - qss(i,k))  !--added by OB because of pathologic cases with the lognormal
1456        qcld(i,k) = qvc(i,k) + zcond(i)
1457   ENDDO
1458   !q_sat
1459   CALL calc_qsat_ecmwf(klon,Tbef(:),qzero(:),pplay(:,k),RTT,1,.FALSE.,qsatl(:,k),zdqs(:))
1460   CALL calc_qsat_ecmwf(klon,Tbef(:),qzero(:),pplay(:,k),RTT,2,.FALSE.,qsats(:,k),zdqs(:))
1461
1462
1463
1464    ! Outputs:
1465    !-------------------------------
1466    ! Precipitation fluxes at layer interfaces
1467    ! + precipitation fractions +
1468    ! temperature and water species tendencies
1469    DO i = 1, klon
1470        psfl(i,k)=zifl(i)
1471        prfl(i,k)=zrfl(i)
1472        pfraclr(i,k)=znebprecipclr(i)
1473        pfracld(i,k)=znebprecipcld(i)
1474        d_ql(i,k) = (1-zfice(i))*zoliq(i)
1475        d_qi(i,k) = zfice(i)*zoliq(i)
1476        d_q(i,k) = zq(i) - qt(i,k)
1477        ! c_iso: same for isotopes
1478        d_t(i,k) = zt(i) - temp(i,k)
1479    ENDDO
1480
1481
1482ENDDO
1483
1484
1485  ! Rain or snow at the surface (depending on the first layer temperature)
1486  DO i = 1, klon
1487      snow(i) = zifl(i)
1488      rain(i) = zrfl(i)
1489      ! c_iso final output
1490  ENDDO
1491
1492  IF (ncoreczq>0) THEN
1493      WRITE(lunout,*)'WARNING : ZQ in LSCP ',ncoreczq,' val < 1.e-15.'
1494  ENDIF
1495
1496
1497
1498
1499END SUBROUTINE lscp
1500!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1501
1502END MODULE lmdz_lscp
Note: See TracBrowser for help on using the repository browser.