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

Last change on this file since 4397 was 4397, checked in by evignon, 16 months ago

Correction d'un beau bug dans l'autoconversion des precip dans lscp (herite de firstilp et
donc encore present dans cette routine). On perd la convergence avec les anciennes versions quand on active lscp
File size: 49.5 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!c_iso: variable initialisation for iso
332
333
334!===============================================================================
335!              BEGINNING OF VERTICAL LOOP FROM TOP TO BOTTOM
336!===============================================================================
337
338ncoreczq=0
339
340DO k = klev, 1, -1
341
342    ! Initialisation temperature and specific humidity
343    DO i = 1, klon
344        zt(i)=t(i,k)
345        zq(i)=q(i,k)
346        !c_iso init of iso
347    ENDDO
348   
349    ! --------------------------------------------------------------------
350    ! P0> Thermalization of precipitation falling from the overlying layer
351    ! --------------------------------------------------------------------
352    ! Computes air temperature variation due to latent heat transported by
353    ! precipitation. Precipitation is then thermalized with the air in the
354    ! layer.
355    ! The precipitation should remain thermalized throughout the different
356    ! thermodynamical transformations. The corresponding water mass should
357    ! be added when calculating the layer's enthalpy change with
358    ! temperature
359    ! ---------------------------------------------------------------------
360   
361    IF (k.LE.klev-1) THEN
362
363        DO i = 1, klon
364 
365            zmair(i)=(paprs(i,k)-paprs(i,k+1))/RG
366            ! no condensed water so cp=cp(vapor+dry air)
367            ! RVTMP2=rcpv/rcpd-1
368            zcpair=RCPD*(1.0+RVTMP2*zq(i))
369            zcpeau=RCPD*RVTMP2
370
371            ! zmqc: precipitation mass that has to be thermalized with
372            ! layer's air so that precipitation at the ground has the
373            ! same temperature as the lowermost layer
374            zmqc(i) = (zrfl(i)+zifl(i))*dtime/zmair(i)
375            ! t(i,k+1)+d_t(i,k+1): new temperature of the overlying layer
376            zt(i) = ( (t(i,k+1)+d_t(i,k+1))*zmqc(i)*zcpeau + zcpair*zt(i) ) &
377             / (zcpair + zmqc(i)*zcpeau)
378
379        ENDDO
380 
381    ELSE
382 
383        DO i = 1, klon
384            zmair(i)=(paprs(i,k)-paprs(i,k+1))/RG
385            zmqc(i) = 0.
386        ENDDO
387 
388    ENDIF
389
390    ! --------------------------------------------------------------------
391    ! P1> Precipitation evaporation/sublimation
392    ! --------------------------------------------------------------------
393    ! A part of the precipitation coming from above is evaporated/sublimated.
394    ! --------------------------------------------------------------------
395
396    IF (iflag_evap_prec.GE.1) THEN
397
398        ! Calculation of saturation specific humidity
399        ! depending on temperature:
400        CALL calc_qsat_ecmwf(klon,zt(:),qzero(:),pplay(:,k),RTT,0,.false.,zqs(:),zdqs(:))
401        ! wrt liquid water
402        CALL calc_qsat_ecmwf(klon,zt(:),qzero(:),pplay(:,k),RTT,1,.false.,qsl(:),dqsl(:))
403        ! wrt ice
404        CALL calc_qsat_ecmwf(klon,zt(:),qzero(:),pplay(:,k),RTT,2,.false.,qsi(:),dqsi(:))
405
406        DO i = 1, klon
407
408            ! if precipitation
409            IF (zrfl(i)+zifl(i).GT.0.) THEN
410
411            ! LTP: we only account for precipitation evaporation in the clear-sky (iflag_evap_prec=4).
412            ! c_iso: likely important to distinguish cs from neb isotope precipitation
413
414                IF (iflag_evap_prec.EQ.4) THEN
415                    zrfl(i) = zrflclr(i)
416                    zifl(i) = ziflclr(i)
417                ENDIF
418
419                IF (iflag_evap_prec.EQ.1) THEN
420                    znebprecip(i)=zneb(i)
421                ELSE
422                    znebprecip(i)=MAX(zneb(i),znebprecip(i))
423                ENDIF
424
425                IF (iflag_evap_prec.EQ.4) THEN
426                ! Max evaporation not to saturate the whole mesh
427                    zqev0 = MAX(0.0, zqs(i)-zq(i))
428                ELSE
429                ! Max evap not to saturate the fraction below the cloud
430                    zqev0 = MAX(0.0, (zqs(i)-zq(i))*znebprecip(i))
431                ENDIF
432
433                ! Evaporation of liquid precipitation coming from above
434                ! dP/dz=beta*(1-q/qsat)*sqrt(P)
435                ! formula from Sundquist 1988, Klemp & Wilhemson 1978
436                ! LTP: evaporation only in the clear sky part (iflag_evap_prec=4)
437
438                IF (iflag_evap_prec.EQ.3) THEN
439                    zqevt = znebprecip(i)*coef_eva*(1.0-zq(i)/qsl(i)) &
440                    *SQRT(zrfl(i)/max(1.e-4,znebprecip(i)))        &
441                    *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
442                ELSE IF (iflag_evap_prec.EQ.4) THEN
443                     zqevt = znebprecipclr(i)*coef_eva*(1.0-zq(i)/qsl(i)) &
444                    *SQRT(zrfl(i)/max(1.e-8,znebprecipclr(i))) &
445                    *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
446                ELSE
447                    zqevt = 1.*coef_eva*(1.0-zq(i)/qsl(i))*SQRT(zrfl(i)) &
448                    *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
449                ENDIF
450
451                zqevt = MAX(0.0,MIN(zqevt,zrfl(i))) &
452                    *RG*dtime/(paprs(i,k)-paprs(i,k+1))
453
454                ! sublimation of the solid precipitation coming from above
455                IF (iflag_evap_prec.EQ.3) THEN
456                    zqevti = znebprecip(i)*coef_eva_i*(1.0-zq(i)/qsi(i)) &
457                    *SQRT(zifl(i)/max(1.e-4,znebprecip(i))) &
458                    *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
459                ELSE IF (iflag_evap_prec.EQ.4) THEN
460                     zqevti = znebprecipclr(i)*coef_eva_i*(1.0-zq(i)/qsi(i)) &
461                    *SQRT(zifl(i)/max(1.e-8,znebprecipclr(i))) &
462                    *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
463                ELSE
464                    zqevti = 1.*coef_eva_i*(1.0-zq(i)/qsi(i))*SQRT(zifl(i)) &
465                    *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
466                ENDIF
467
468                zqevti = MAX(0.0,MIN(zqevti,zifl(i))) &
469                *RG*dtime/(paprs(i,k)-paprs(i,k+1))
470
471                ! A. JAM
472                ! Evaporation limit: we ensure that the layer's fraction below
473                ! the cloud or the whole mesh (depending on iflag_evap_prec)
474                ! does not reach saturation. In this case, we
475                ! redistribute zqev0 conserving the ratio liquid/ice
476
477                IF (zqevt+zqevti.GT.zqev0) THEN
478                    zqev=zqev0*zqevt/(zqevt+zqevti)
479                    zqevi=zqev0*zqevti/(zqevt+zqevti)
480                ELSE
481                    zqev=zqevt
482                    zqevi=zqevti
483                ENDIF
484
485
486                ! New solid and liquid precipitation fluxes after evap and sublimation
487                zrfln(i) = Max(0.,zrfl(i) - zqev*(paprs(i,k)-paprs(i,k+1))  &
488                         /RG/dtime)
489                zifln(i) = Max(0.,zifl(i) - zqevi*(paprs(i,k)-paprs(i,k+1)) &
490                         /RG/dtime)
491
492
493                ! vapor, temperature, precip fluxes update
494                zq(i) = zq(i) - (zrfln(i)+zifln(i)-zrfl(i)-zifl(i)) &
495                * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
496                ! precip thermalization
497                zmqc(i) = zmqc(i) + (zrfln(i)+zifln(i)-zrfl(i)-zifl(i)) &
498                * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
499                zt(i) = zt(i) + (zrfln(i)-zrfl(i))        &
500                * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime    &
501                * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i))) &
502                + (zifln(i)-zifl(i))                      &
503                * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime    &
504                * RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))
505
506                ! New values of liquid and solid precipitation
507                zrfl(i) = zrfln(i)
508                zifl(i) = zifln(i)
509
510                ! c_iso here call_reevap that updates isotopic zrfl, zifl (in inout)
511                !                         due to evap + sublim                                     
512                                           
513
514                IF (iflag_evap_prec.EQ.4) THEN
515                    zrflclr(i) = zrfl(i)
516                    ziflclr(i) = zifl(i)
517                    IF(zrflclr(i) + ziflclr(i).LE.0) THEN
518                        znebprecipclr(i) = 0.0
519                    ENDIF
520                    zrfl(i) = zrflclr(i) + zrflcld(i)
521                    zifl(i) = ziflclr(i) + ziflcld(i)
522                ENDIF
523
524                ! c_iso duplicate for isotopes or loop on isotopes
525
526                ! Melting:
527                zmelt = ((zt(i)-273.15)/(ztfondue-273.15)) ! JYG
528                ! precip fraction that is melted
529                zmelt = MIN(MAX(zmelt,0.),1.)
530
531                ! Ice melting
532                IF (iflag_evap_prec.EQ.4) THEN
533                    zrflclr(i)=zrflclr(i)+zmelt*ziflclr(i)
534                    zrflcld(i)=zrflcld(i)+zmelt*ziflcld(i)
535                    zrfl(i)=zrflclr(i)+zrflcld(i)
536                ELSE
537                    zrfl(i)=zrfl(i)+zmelt*zifl(i)
538                ENDIF
539                ! c_iso: melting of isotopic precipi with zmelt (no fractionation)
540
541                ! Latent heat of melting with precipitation thermalization
542                zt(i)=zt(i)-zifl(i)*zmelt*(RG*dtime)/(paprs(i,k)-paprs(i,k+1)) &
543                *RLMLT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))
544               
545                IF (iflag_evap_prec.EQ.4) THEN
546                    ziflclr(i)=ziflclr(i)*(1.-zmelt)
547                    ziflcld(i)=ziflcld(i)*(1.-zmelt)
548                    zifl(i)=ziflclr(i)+ziflcld(i)
549                ELSE
550                    zifl(i)=zifl(i)*(1.-zmelt)
551                ENDIF
552                !c_iso: same for isotopic precip.
553
554            ELSE
555                ! if no precip, we reinitialize the cloud fraction used for the precip to 0
556                znebprecip(i)=0.
557
558            ENDIF ! (zrfl(i)+zifl(i).GT.0.)
559
560        ENDDO
561
562    ENDIF ! (iflag_evap_prec>=1)
563
564    ! --------------------------------------------------------------------
565    ! End precip evaporation
566    ! --------------------------------------------------------------------
567   
568    ! Calculation of qsat, L/Cp*dqsat/dT and ncoreczq counter
569    !-------------------------------------------------------
570
571     qtot(:)=zq(:)+zmqc(:)
572     CALL calc_qsat_ecmwf(klon,zt(:),qtot(:),pplay(:,k),RTT,0,.false.,zqs(:),zdqs(:))
573     DO i = 1, klon
574            zdelta = MAX(0.,SIGN(1.,RTT-zt(i)))
575            zdqsdT_raw(i) = zdqs(i)*RCPD*(1.0+RVTMP2*zq(i)) / (RLVTT*(1.-zdelta) + RLSTT*zdelta)
576            IF (zq(i) .LT. 1.e-15) THEN
577                ncoreczq=ncoreczq+1
578                zq(i)=1.e-15
579            ENDIF
580            ! c_iso: do something similar for isotopes
581
582     ENDDO
583   
584    ! --------------------------------------------------------------------
585    ! P2> Cloud formation
586    !---------------------------------------------------------------------
587    !
588    ! Unlike fisrtilp, we always assume a 'fractional cloud' approach
589    ! i.e. clouds occupy only a fraction of the mesh (the subgrid distribution
590    ! is prescribed and depends on large scale variables and boundary layer
591    ! properties)
592    ! The decrease in condensed part due tu latent heating is taken into
593    ! account
594    ! -------------------------------------------------------------------
595
596        ! P2.1> With the PDFs (log-normal, bigaussian)
597        ! cloud propertues calculation with the initial values of t and q
598        ! ----------------------------------------------------------------
599
600        ! initialise gammasat and qincloud_mpc
601        gammasat(:)=1.
602        qincloud_mpc(:)=0.
603
604        IF (iflag_cld_th.GE.5) THEN
605
606             IF (iflag_mpc_bl .LT. 1) THEN
607
608                IF (iflag_cloudth_vert.LE.2) THEN
609
610                    CALL cloudth(klon,klev,k,ztv,                  &
611                         zq,zqta,fraca,                            &
612                         qcloud,ctot,zpspsk,paprs,pplay,ztla,zthl, &
613                         ratqs,zqs,t)
614
615                ELSEIF (iflag_cloudth_vert.GE.3 .AND. iflag_cloudth_vert.LE.5) THEN
616
617                    CALL cloudth_v3(klon,klev,k,ztv,                        &
618                         zq,zqta,fraca,                                     &
619                         qcloud,ctot,ctot_vol,zpspsk,paprs,pplay,ztla,zthl, &
620                         ratqs,zqs,t)
621
622                !Jean Jouhaud's version, Decembre 2018
623                !-------------------------------------
624
625                ELSEIF (iflag_cloudth_vert.EQ.6) THEN
626
627                    CALL cloudth_v6(klon,klev,k,ztv,                        &
628                         zq,zqta,fraca,                                     &
629                         qcloud,ctot,ctot_vol,zpspsk,paprs,pplay,ztla,zthl, &
630                         ratqs,zqs,t)
631
632                ENDIF
633
634         ELSE
635                ! cloudth_v3 + cold microphysical considerations
636                ! + stationary mixed-phase cloud model
637
638                    CALL cloudth_mpc(klon,klev,k,iflag_mpc_bl,mpc_bl_points, &
639                         t,ztv,zq,zqta,fraca, zpspsk,                      &
640                         paprs,pplay,ztla,zthl,ratqs,zqs,psfl,             &
641                         qcloud,qincloud_mpc,icefrac_mpc,ctot,ctot_vol)
642
643         ENDIF ! iflag_mpc_bl
644
645                DO i=1,klon
646                    rneb(i,k)=ctot(i,k)
647                    rneblsvol(i,k)=ctot_vol(i,k)
648                    zqn(i)=qcloud(i)
649                ENDDO
650
651        ENDIF
652
653        IF (iflag_cld_th .LE. 4) THEN
654           
655                ! lognormal
656            lognormale = .TRUE.
657
658        ELSEIF (iflag_cld_th .GE. 6) THEN
659           
660                ! lognormal distribution when no thermals
661            lognormale = fraca(:,k) < 1e-10
662
663        ELSE
664                ! When iflag_cld_th=5, we always assume
665                ! bi-gaussian distribution
666            lognormale = .FALSE.
667       
668        ENDIF
669
670        DT(:) = 0.
671        n_i(:)=0
672        Tbef(:)=zt(:)
673        qlbef(:)=0.
674
675        ! Treatment of non-boundary layer clouds (lognormale)
676        ! condensation with qsat(T) variation (adaptation)
677        ! Iterative Loop to converge towards qsat
678
679        DO iter=1,iflag_fisrtilp_qsat+1
680
681                DO i=1,klon
682
683                    ! convergence = .true. until when convergence is not satisfied
684                    convergence(i)=ABS(DT(i)).GT.DDT0
685
686                    IF ((convergence(i) .OR. (n_i(i) .EQ. 0)) .AND. lognormale(i)) THEN
687
688                        ! if not convergence:
689
690                        ! P2.2.1> cloud fraction and condensed water mass calculation
691                        ! Calculated variables:
692                        ! rneb : cloud fraction
693                        ! zqn : total water within the cloud
694                        ! zcond: mean condensed water within the mesh
695                        ! rhcl: clear-sky relative humidity
696                        !---------------------------------------------------------------
697
698                        ! new temperature:
699                        Tbef(i)=Tbef(i)+DT(i)
700
701                        ! Rneb, qzn and zcond for lognormal PDFs
702                        qtot(i)=zq(i)+zmqc(i)
703
704                      ENDIF
705
706                  ENDDO
707
708                  ! Calculation of saturation specific humidity and ice fraction
709                  CALL calc_qsat_ecmwf(klon,Tbef(:),qtot(:),pplay(:,k),RTT,0,.false.,zqs(:),zdqs(:))
710                  CALL calc_gammasat(klon,Tbef(:),qtot(:),pplay(:,k),gammasat(:),dgammasatdt(:))
711                  ! saturation may occur at a humidity different from qsat (gamma qsat), so gamma correction for dqs
712                  zdqs(:) = gammasat(:)*zdqs(:)+zqs(:)*dgammasatdt(:)
713                  CALL icefrac_lscp(klon, zt(:),pplay(:,k)/paprs(:,1),zfice(:),dzfice(:))
714
715                  DO i=1,klon
716
717                      IF ((convergence(i) .OR. (n_i(i) .EQ. 0)) .AND. lognormale(i)) THEN
718
719                        zpdf_sig(i)=ratqs(i,k)*zq(i)
720                        zpdf_k(i)=-sqrt(log(1.+(zpdf_sig(i)/zq(i))**2))
721                        zpdf_delta(i)=log(zq(i)/(gammasat(i)*zqs(i)))
722                        zpdf_a(i)=zpdf_delta(i)/(zpdf_k(i)*sqrt(2.))
723                        zpdf_b(i)=zpdf_k(i)/(2.*sqrt(2.))
724                        zpdf_e1(i)=zpdf_a(i)-zpdf_b(i)
725                        zpdf_e1(i)=sign(min(ABS(zpdf_e1(i)),5.),zpdf_e1(i))
726                        zpdf_e1(i)=1.-erf(zpdf_e1(i))
727                        zpdf_e2(i)=zpdf_a(i)+zpdf_b(i)
728                        zpdf_e2(i)=sign(min(ABS(zpdf_e2(i)),5.),zpdf_e2(i))
729                        zpdf_e2(i)=1.-erf(zpdf_e2(i))
730
731                        !--ice sursaturation by Audran
732                        IF ((.NOT.ok_ice_sursat).OR.(Tbef(i).GT.t_glace_min)) THEN
733
734                          IF (zpdf_e1(i).LT.1.e-10) THEN
735                              rneb(i,k)=0.
736                              zqn(i)=gammasat(i)*zqs(i)
737                          ELSE
738                              rneb(i,k)=0.5*zpdf_e1(i)
739                              zqn(i)=zq(i)*zpdf_e2(i)/zpdf_e1(i)
740                          ENDIF
741
742                          rnebss(i,k)=0.0   !--added by OB (needed because of the convergence loop on time)
743                          fcontrN(i,k)=0.0  !--idem
744                          fcontrP(i,k)=0.0  !--idem
745                          qss(i,k)=0.0      !--idem
746
747                        ELSE
748
749                        !------------------------------------
750                        ! ICE SUPERSATURATION
751                        !------------------------------------
752
753                        CALL ice_sursat(pplay(i,k), paprs(i,k)-paprs(i,k+1), dtime, i, k, t(i,k), zq(i), &
754                             gamma_ss(i,k), zqs(i), Tbef(i), rneb_seri(i,k), ratqs(i,k),                 &
755                             rneb(i,k), zqn(i), rnebss(i,k), qss(i,k),                                   &
756                             Tcontr(i,k), qcontr(i,k), qcontr2(i,k), fcontrN(i,k), fcontrP(i,k) )
757
758                        ENDIF ! ((flag_ice_sursat.eq.0).or.(Tbef(i).gt.t_glace_min))
759
760                        ! If vertical heterogeneity, change fraction by volume as well
761                        IF (iflag_cloudth_vert.GE.3) THEN
762                            ctot_vol(i,k)=rneb(i,k)
763                            rneblsvol(i,k)=ctot_vol(i,k)
764                        ENDIF
765
766
767                       ! P2.2.2> Approximative calculation of temperature variation DT
768                       !           due to condensation.
769                       ! Calculated variables:
770                       ! dT : temperature change due to condensation
771                       !---------------------------------------------------------------
772
773               
774                        IF (zfice(i).LT.1) THEN
775                            cste=RLVTT
776                        ELSE
777                            cste=RLSTT
778                        ENDIF
779
780                        qlbef(i)=max(0.,zqn(i)-zqs(i))
781                        num = -Tbef(i)+zt(i)+rneb(i,k)*((1-zfice(i))*RLVTT &
782                              +zfice(i)*RLSTT)/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))*qlbef(i)
783                        denom = 1.+rneb(i,k)*((1-zfice(i))*RLVTT+zfice(i)*RLSTT)/cste*zdqs(i) &
784                              -(RLSTT-RLVTT)/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))*rneb(i,k)      &
785                              *qlbef(i)*dzfice(i)
786                        DT(i)=num/denom
787                        n_i(i)=n_i(i)+1
788
789                    ENDIF ! end convergence
790
791                ENDDO     ! end loop on i
792
793        ENDDO         ! iter=1,iflag_fisrtilp_qsat+1
794
795        ! P2.3> Final quantities calculation
796        ! Calculated variables:
797        !   rneb : cloud fraction
798        !   zcond: mean condensed water in the mesh
799        !   zqn  : mean water vapor in the mesh
800        !   zt   : temperature
801        !   rhcl : clear-sky relative humidity
802        ! ----------------------------------------------------------------
803
804        ! Water vapor update, Phase determination and subsequent latent heat exchange
805
806        ! Partition function in stratiform clouds (will be overwritten in boundary-layer MPCs)
807        CALL icefrac_lscp(klon,zt(:),pplay(:,k)/paprs(:,1),zfice(:), dzfice(:))
808
809        DO i=1, klon
810
811
812            IF (mpc_bl_points(i,k) .GT. 0) THEN
813                zcond(i) = MAX(0.0,qincloud_mpc(i))*rneb(i,k)
814                ! following line is very strange and probably wrong
815                rhcl(i,k)= (zqs(i)+zq(i))/2./zqs(i)
816                ! water vapor update and partition function if thermals
817                zq(i) = zq(i) - zcond(i)       
818                zfice(i)=icefrac_mpc(i,k)
819
820            ELSE
821
822                ! Checks on rneb, rhcl and zqn
823                IF (rneb(i,k) .LE. 0.0) THEN
824                    zqn(i) = 0.0
825                    rneb(i,k) = 0.0
826                    zcond(i) = 0.0
827                    rhcl(i,k)=zq(i)/zqs(i)
828                ELSE IF (rneb(i,k) .GE. 1.0) THEN
829                    zqn(i) = zq(i)
830                    rneb(i,k) = 1.0
831                    zcond(i) = MAX(0.0,zqn(i)-gammasat(i)*zqs(i))
832                    rhcl(i,k)=1.0
833                ELSE
834                    zcond(i) = MAX(0.0,zqn(i)-gammasat(i)*zqs(i))*rneb(i,k)
835                    ! following line is very strange and probably wrong:
836                    rhcl(i,k)=(zqs(i)+zq(i))/2./zqs(i)
837                ENDIF
838
839                ! water vapor update
840                zq(i) = zq(i) - zcond(i)
841
842            ENDIF
843
844            ! c_iso : routine that computes in-cloud supersaturation
845            ! c_iso condensation of isotopes (zcond, zsursat, zfice, zq in input)
846
847            ! temperature update due to phase change
848            zt(i) = zt(i) + (1.-zfice(i))*zcond(i) &
849            &     * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i))) &
850                  +zfice(i)*zcond(i) * RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i)))
851        ENDDO
852
853        ! If vertical heterogeneity, change volume fraction
854        IF (iflag_cloudth_vert .GE. 3) THEN
855          ctot_vol(1:klon,k)=min(max(ctot_vol(1:klon,k),0.),1.)
856          rneblsvol(1:klon,k)=ctot_vol(1:klon,k)
857        ENDIF
858
859    ! ----------------------------------------------------------------
860    ! P3> Precipitation formation
861    ! ----------------------------------------------------------------
862
863    ! LTP:
864    IF (iflag_evap_prec==4) THEN
865
866        !Partitionning between precipitation coming from clouds and that coming from CS
867
868        !0) Calculate tot_zneb, total cloud fraction above the cloud
869        !assuming a maximum-random overlap (voir Jakob and Klein, 2000)
870       
871        DO i=1, klon
872                tot_znebn(i) = 1 - (1-tot_zneb(i))*(1 - max(rneb(i,k),zneb(i))) &
873                        /(1-min(zneb(i),1-smallestreal))
874                d_tot_zneb(i) = tot_znebn(i) - tot_zneb(i)
875                tot_zneb(i) = tot_znebn(i)
876
877
878                !1) Cloudy to clear air
879                d_znebprecip_cld_clr(i) = znebprecipcld(i) - min(rneb(i,k),znebprecipcld(i))
880                IF (znebprecipcld(i) .GT. 0) THEN
881                        d_zrfl_cld_clr(i) = d_znebprecip_cld_clr(i)/znebprecipcld(i)*zrflcld(i)
882                        d_zifl_cld_clr(i) = d_znebprecip_cld_clr(i)/znebprecipcld(i)*ziflcld(i)
883                ELSE
884                        d_zrfl_cld_clr(i) = 0.
885                        d_zifl_cld_clr(i) = 0.
886                ENDIF
887
888                !2) Clear to cloudy air
889                d_znebprecip_clr_cld(i) = max(0., min(znebprecipclr(i), rneb(i,k) - d_tot_zneb(i) - zneb(i)))
890                IF (znebprecipclr(i) .GT. 0) THEN
891                        d_zrfl_clr_cld(i) = d_znebprecip_clr_cld(i)/znebprecipclr(i)*zrflclr(i)
892                        d_zifl_clr_cld(i) = d_znebprecip_clr_cld(i)/znebprecipclr(i)*ziflclr(i)
893                ELSE
894                        d_zrfl_clr_cld(i) = 0.
895                        d_zifl_clr_cld(i) = 0.
896                ENDIF
897
898                !Update variables
899                znebprecipcld(i) = znebprecipcld(i) + d_znebprecip_clr_cld(i) - d_znebprecip_cld_clr(i)
900                znebprecipclr(i) = znebprecipclr(i) + d_znebprecip_cld_clr(i) - d_znebprecip_clr_cld(i)
901                zrflcld(i) = zrflcld(i) + d_zrfl_clr_cld(i) - d_zrfl_cld_clr(i)
902                ziflcld(i) = ziflcld(i) + d_zifl_clr_cld(i) - d_zifl_cld_clr(i)
903                zrflclr(i) = zrflclr(i) + d_zrfl_cld_clr(i) - d_zrfl_clr_cld(i)
904                ziflclr(i) = ziflclr(i) + d_zifl_cld_clr(i) - d_zifl_clr_cld(i)
905
906                ! c_iso  do the same thing for isotopes precip
907        ENDDO
908    ENDIF
909
910    ! Initialisation of zoliq and radliq variables
911
912    DO i = 1, klon
913            zoliq(i) = zcond(i)
914            zoliqi(i)= zoliq(i)*zfice(i)
915            zoliql(i)= zoliq(i)*(1.-zfice(i))
916            ! c_iso : initialisation of zoliq* also for isotopes
917            zrho(i,k)  = pplay(i,k) / zt(i) / RD
918            zdz(i)   = (paprs(i,k)-paprs(i,k+1)) / (zrho(i,k)*RG)
919            iwc(i)   = 0.
920            zneb(i)  = MAX(rneb(i,k),seuil_neb)
921            radliq(i,k) = zoliq(i)/REAL(ninter+1)
922            radicefrac(i,k) = zfice(i)
923            radliqi(i,k)=zoliq(i)*zfice(i)/REAL(ninter+1)
924            radliql(i,k)=zoliq(i)*(1.-zfice(i))/REAL(ninter+1)
925    ENDDO
926
927    ! Autoconversion
928    ! -------------------------------------------------------------------------------
929     
930    DO n = 1, ninter
931
932        DO i=1,klon
933            IF (rneb(i,k).GT.0.0) THEN
934                iwc(i) = zrho(i,k) * zoliqi(i) / zneb(i) ! in-cloud ice condensate content
935            ENDIF
936        ENDDO
937
938        CALL fallice_velocity(klon,iwc(:),zt(:),zrho(:,k),pplay(:,k),ptconv(:,k),velo(:,k))
939
940        DO i = 1, klon
941
942            IF (rneb(i,k).GT.0.0) THEN
943
944                ! Initialization of zrain, zsnow and zprecip:
945                zrain=0.
946                zsnow=0.
947                zprecip=0.
948                ! c_iso same init for isotopes. Externalisation?
949
950                IF (zneb(i).EQ.seuil_neb) THEN
951                    zprecip = 0.0
952                    zsnow = 0.0
953                    zrain= 0.0
954                ELSE
955                    ! water quantity to remove: zchau (Sundqvist, 1978)
956                    ! same thing for the ice: zfroi (Zender & Kiehl, 1997)
957                    IF (ptconv(i,k)) THEN ! if convective point
958                        zcl=cld_lc_con
959                        zct=1./cld_tau_con
960                        zfroi = dtime/REAL(ninter)/zdz(i)*zoliqi(i)*velo(i,k)
961                    ELSE
962                        zcl=cld_lc_lsc
963                        zct=1./cld_tau_lsc
964                        zfroi = dtime/REAL(ninter)/zdz(i)*zoliqi(i) &   ! dqice/dt=1/rho*d(rho*wice*qice)/dz
965                            *velo(i,k)
966                    ENDIF
967
968                    ! if vertical heterogeneity is taken into account, we use
969                    ! the "true" volume fraction instead of a modified
970                    ! surface fraction (which is larger and artificially
971                    ! reduces the in-cloud water).
972
973
974                    IF ((iflag_cloudth_vert.GE.3).AND.(iflag_rain_incloud_vol.EQ.1)) THEN
975                        zchau = zct   *dtime/REAL(ninter) * zoliql(i) &
976                        *(1.0-EXP(-(zoliql(i)/ctot_vol(i,k)/zcl)**2))
977                    ELSE
978                        zchau = zct   *dtime/REAL(ninter) * zoliql(i) &
979                        *(1.0-EXP(-(zoliql(i)/zneb(i)/zcl)**2))        ! dqliq/dt=-qliq/tau*(1-exp(-qcin/clw)**2)
980                    ENDIF
981
982                    zrain   = MIN(MAX(zchau,0.0),zoliql(i))
983                    zsnow   = MIN(MAX(zfroi,0.0),zoliqi(i))
984                    zprecip = MAX(zrain + zsnow,0.0)
985
986                ENDIF
987
988                zoliql(i) = MAX(zoliq(i)*(1.-zfice(i))-1.*zrain  , 0.0)
989                zoliqi(i) = MAX(zoliq(i)*zfice(i)-1.*zsnow  , 0.0)
990                zoliq(i)  = MAX(zoliq(i)-zprecip , 0.0)
991
992                ! c_iso: call isotope_conversion (for readibility) or duplicate
993
994                radliq(i,k) = radliq(i,k) + zoliq(i)/REAL(ninter+1)
995                radliql(i,k) = radliql(i,k) + zoliql(i)/REAL(ninter+1)
996                radliqi(i,k) = radliqi(i,k) + zoliqi(i)/REAL(ninter+1)
997
998            ENDIF ! rneb >0
999
1000        ENDDO  ! i = 1,klon
1001
1002    ENDDO      ! n = 1,niter
1003
1004    ! Precipitation flux calculation
1005
1006    DO i = 1, klon
1007           
1008            IF (iflag_evap_prec.EQ.4) THEN
1009                ziflprev(i)=ziflcld(i)
1010            ELSE
1011                ziflprev(i)=zifl(i)*zneb(i)
1012            ENDIF
1013
1014            IF (rneb(i,k) .GT. 0.0) THEN
1015
1016               ! CR&JYG: We account for the Wegener-Findeisen-Bergeron process in the precipitation flux:
1017               ! If T<0C, liquid precip are converted into ice, which leads to an increase in
1018               ! temperature DeltaT. The effect of DeltaT on condensates and precipitation is roughly
1019               ! taken into account through a linearization of the equations and by approximating
1020               ! the liquid precipitation process with a "threshold" process. We assume that
1021               ! condensates are not modified during this operation. Liquid precipitation is
1022               ! removed (in the limit DeltaT<273.15-T). Solid precipitation increases.
1023               ! Water vapor increases as well
1024               ! Note that compared to fisrtilp, we always assume iflag_bergeron=2
1025
1026                    zqpreci(i)=(zcond(i)-zoliq(i))*zfice(i)
1027                    zqprecl(i)=(zcond(i)-zoliq(i))*(1.-zfice(i))
1028                    zcp=RCPD*(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i)))
1029                    coef1 = rneb(i,k)*RLSTT/zcp*zdqsdT_raw(i)
1030                    ! Computation of DT if all the liquid precip freezes
1031                    DeltaT = RLMLT*zqprecl(i) / (zcp*(1.+coef1))
1032                    ! T should not exceed the freezing point
1033                    ! that is Delta > RTT-zt(i)
1034                    DeltaT = max( min( RTT-zt(i), DeltaT) , 0. )
1035                    zt(i)  = zt(i) + DeltaT
1036                    ! water vaporization due to temp. increase
1037                    Deltaq = rneb(i,k)*zdqsdT_raw(i)*DeltaT
1038                    ! we add this vaporized water to the total vapor and we remove it from the precip
1039                    zq(i) = zq(i) +  Deltaq
1040                    ! The three "max" lines herebelow protect from rounding errors
1041                    zcond(i) = max( zcond(i)- Deltaq, 0. )
1042                    ! liquid precipitation converted to ice precip
1043                    Deltaqprecl = -zcp/RLMLT*(1.+coef1)*DeltaT
1044                    zqprecl(i) = max( zqprecl(i) + Deltaqprecl, 0. )
1045                    ! iced water budget
1046                    zqpreci(i) = max (zqpreci(i) - Deltaqprecl - Deltaq, 0.)
1047
1048                    ! c_iso : mv here condensation of isotopes + redispatchage en precip
1049
1050                IF (iflag_evap_prec.EQ.4) THEN
1051                    zrflcld(i) = zrflcld(i)+zqprecl(i) &
1052                    *(paprs(i,k)-paprs(i,k+1))/(RG*dtime)
1053                    ziflcld(i) = ziflcld(i)+ zqpreci(i) &
1054                        *(paprs(i,k)-paprs(i,k+1))/(RG*dtime)
1055                    znebprecipcld(i) = rneb(i,k)
1056                    zrfl(i) = zrflcld(i) + zrflclr(i)
1057                    zifl(i) = ziflcld(i) + ziflclr(i)
1058                ELSE
1059                    zrfl(i) = zrfl(i)+ zqprecl(i)        &
1060                    *(paprs(i,k)-paprs(i,k+1))/(RG*dtime)
1061                    zifl(i) = zifl(i)+ zqpreci(i)        &
1062                   *(paprs(i,k)-paprs(i,k+1))/(RG*dtime)
1063                ENDIF
1064                ! c_iso : same for isotopes
1065 
1066           ENDIF ! rneb>0
1067
1068   ENDDO
1069
1070    ! LTP: limit of surface cloud fraction covered by precipitation when the local intensity of the flux is below rain_int_min
1071    ! if iflag_evap_pre=4
1072    IF (iflag_evap_prec.EQ.4) THEN
1073
1074        DO i=1,klon
1075
1076            IF ((zrflclr(i) + ziflclr(i)) .GT. 0. ) THEN
1077                znebprecipclr(i) = min(znebprecipclr(i),max(zrflclr(i)/ &
1078                (MAX(znebprecipclr(i),seuil_neb)*rain_int_min), ziflclr(i)/(MAX(znebprecipclr(i),seuil_neb)*rain_int_min)))
1079            ELSE
1080                znebprecipclr(i)=0.0
1081            ENDIF
1082
1083            IF ((zrflcld(i) + ziflcld(i)) .GT. 0.) THEN
1084                znebprecipcld(i) = min(znebprecipcld(i), max(zrflcld(i)/ &
1085                (MAX(znebprecipcld(i),seuil_neb)*rain_int_min), ziflcld(i)/(MAX(znebprecipcld(i),seuil_neb)*rain_int_min)))
1086            ELSE
1087                znebprecipcld(i)=0.0
1088            ENDIF
1089
1090        ENDDO
1091
1092    ENDIF
1093
1094    ! End of precipitation formation
1095    ! --------------------------------
1096
1097    ! Outputs:
1098    ! Precipitation fluxes at layer interfaces
1099    ! and temperature and water species tendencies
1100    DO i = 1, klon
1101        psfl(i,k)=zifl(i)
1102        prfl(i,k)=zrfl(i)
1103        d_ql(i,k) = (1-zfice(i))*zoliq(i)
1104        d_qi(i,k) = zfice(i)*zoliq(i)
1105        d_q(i,k) = zq(i) - q(i,k)
1106        ! c_iso: same for isotopes
1107        d_t(i,k) = zt(i) - t(i,k)
1108    ENDDO
1109
1110    ! Calculation of the concentration of condensates seen by the radiation scheme
1111    ! for liquid phase, we take radliql
1112    ! for ice phase, we take radliqi if we neglect snowfall, otherwise (ok_radliq_snow=true)
1113    ! we recaulate radliqi to account for contributions from the precipitation flux
1114
1115    IF ((ok_radliq_snow) .AND. (k .LT. klev)) THEN
1116        ! for the solid phase (crystals + snowflakes)
1117        ! we recalculate radliqi by summing
1118        ! the ice content calculated in the mesh
1119        ! + the contribution of the non-evaporated snowfall
1120        ! from the overlying layer
1121        DO i=1,klon
1122           IF (ziflprev(i) .NE. 0.0) THEN
1123              radliqi(i,k)=zoliq(i)*zfice(i)+zqpreci(i)+ziflprev(i)/zrho(i,k+1)/velo(i,k+1)
1124           ELSE
1125              radliqi(i,k)=zoliq(i)*zfice(i)+zqpreci(i)
1126           ENDIF
1127           radliq(i,k)=radliql(i,k)+radliqi(i,k)
1128        ENDDO
1129    ENDIF
1130
1131    ! caculate the percentage of ice in "radliq" so cloud+precip seen by the radiation scheme
1132    DO i=1,klon
1133        IF (radliq(i,k) .GT. 0.) THEN
1134            radicefrac(i,k)=MIN(MAX(radliqi(i,k)/radliq(i,k),0.),1.)
1135        ENDIF
1136    ENDDO
1137
1138    ! ----------------------------------------------------------------
1139    ! P4> Wet scavenging
1140    ! ----------------------------------------------------------------
1141
1142    !Scavenging through nucleation in the layer
1143
1144    DO i = 1,klon
1145       
1146        IF(zcond(i).GT.zoliq(i)+1.e-10) THEN
1147            beta(i,k) = (zcond(i)-zoliq(i))/zcond(i)/dtime
1148        ELSE
1149            beta(i,k) = 0.
1150        ENDIF
1151
1152        zprec_cond(i) = MAX(zcond(i)-zoliq(i),0.0)*(paprs(i,k)-paprs(i,k+1))/RG
1153
1154        IF (rneb(i,k).GT.0.0.AND.zprec_cond(i).GT.0.) THEN
1155
1156            IF (t(i,k) .GE. t_glace_min) THEN
1157                zalpha_tr = a_tr_sca(3)
1158            ELSE
1159                zalpha_tr = a_tr_sca(4)
1160            ENDIF
1161
1162            zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i))
1163            frac_nucl(i,k)= 1.-zneb(i)*zfrac_lessi
1164
1165            ! Nucleation with a factor of -1 instead of -0.5
1166            zfrac_lessi = 1. - EXP(-zprec_cond(i)/zneb(i))
1167
1168        ENDIF
1169
1170    ENDDO
1171
1172    ! Scavenging through impaction in the underlying layer
1173
1174    DO kk = k-1, 1, -1
1175
1176        DO i = 1, klon
1177
1178            IF (rneb(i,k).GT.0.0.AND.zprec_cond(i).GT.0.) THEN
1179
1180                IF (t(i,kk) .GE. t_glace_min) THEN
1181                    zalpha_tr = a_tr_sca(1)
1182                ELSE
1183                    zalpha_tr = a_tr_sca(2)
1184                ENDIF
1185
1186              zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i))
1187              frac_impa(i,kk)= 1.-zneb(i)*zfrac_lessi
1188
1189            ENDIF
1190
1191        ENDDO
1192
1193    ENDDO
1194
1195    !--save some variables for ice supersaturation
1196    !
1197    DO i = 1, klon
1198        ! for memory
1199        rneb_seri(i,k) = rneb(i,k)
1200
1201        ! for diagnostics
1202        rnebclr(i,k) = 1.0 - rneb(i,k) - rnebss(i,k)
1203
1204        qvc(i,k) = zqs(i) * rneb(i,k)
1205        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
1206        qcld(i,k) = qvc(i,k) + zcond(i)
1207   ENDDO
1208   !q_sat
1209   CALL calc_qsat_ecmwf(klon,Tbef(:),qzero(:),pplay(:,k),RTT,1,.false.,qsatl(:,k),zdqs(:))
1210   CALL calc_qsat_ecmwf(klon,Tbef(:),qzero(:),pplay(:,k),RTT,2,.false.,qsats(:,k),zdqs(:))
1211
1212ENDDO
1213
1214!======================================================================
1215!                      END OF VERTICAL LOOP
1216!======================================================================
1217
1218  ! Rain or snow at the surface (depending on the first layer temperature)
1219  DO i = 1, klon
1220      snow(i) = zifl(i)
1221      rain(i) = zrfl(i)
1222      ! c_iso final output
1223  ENDDO
1224
1225  IF (ncoreczq>0) THEN
1226      WRITE(lunout,*)'WARNING : ZQ in LSCP ',ncoreczq,' val < 1.e-15.'
1227  ENDIF
1228
1229END SUBROUTINE lscp
1230!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1231
1232END MODULE lscp_mod
Note: See TracBrowser for help on using the repository browser.