source: LMDZ6/branches/Ocean_skin/libf/phylmd/lscp_mod.F90 @ 5185

Last change on this file since 5185 was 4368, checked in by lguez, 2 years ago

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