source: LMDZ6/trunk/libf/phylmd/lscp_mod.F90 @ 4660

Last change on this file since 4660 was 4654, checked in by fhourdin, 15 months ago

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