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

Last change on this file since 4380 was 4380, checked in by evignon, 21 months ago

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