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

Last change on this file since 4420 was 4420, checked in by evignon, 17 months ago

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