source: LMDZ6/trunk/libf/phylmd/lmdz_lscp.F90 @ 4830

Last change on this file since 4830 was 4830, checked in by evignon, 3 months ago

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