source: LMDZ6/trunk/libf/phylmd/lmdz_lscp.F90 @ 4913

Last change on this file since 4913 was 4913, checked in by evignon, 4 weeks ago

amelioration de poprecip avec:

  • formule du rayon des cristaux
  • qvap du ciel clair pour l'évap des precip

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