source: LMDZ6/branches/contrails/libf/phylmd/lmdz_lscp.f90 @ 5428

Last change on this file since 5428 was 5413, checked in by evignon, 9 days ago

fin de l'externalisation des precips dans lscp.

  1. Borella, E Vignon
File size: 49.0 KB
Line 
1MODULE lmdz_lscp
2
3IMPLICIT NONE
4
5CONTAINS
6
7!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
8SUBROUTINE lscp(klon,klev,dtime,missing_val,            &
9     paprs, pplay, omega, temp, qt, ql_seri, qi_seri,   &
10     ptconv, ratqs, sigma_qtherm,                       &
11     d_t, d_q, d_ql, d_qi, rneb, rneblsvol,             &
12     pfraclr, pfracld,                                  &
13     cldfraliq, sigma2_icefracturb,mean_icefracturb,    &
14     radocond, radicefrac, rain, snow,                  &
15     frac_impa, frac_nucl, beta,                        &
16     prfl, psfl, rhcl, qta, fraca,                      &
17     tv, pspsk, tla, thl, iflag_cld_th,                 &
18     iflag_ice_thermo, distcltop, temp_cltop,           &
19     tke, tke_dissip,                                   &
20     cell_area,                                         &
21     cf_seri, rvc_seri, u_seri, v_seri,                 &
22     qsub, qissr, qcld, subfra, issrfra, gamma_cond,    &
23     dcf_sub, dcf_con, dcf_mix,          &
24     dqi_adj, dqi_sub, dqi_con, dqi_mix, dqvc_adj,      &
25     dqvc_sub, dqvc_con, dqvc_mix, qsatl, qsati,        &
26     Tcontr, qcontr, qcontr2, fcontrN, fcontrP, dcf_avi,&
27     dqi_avi, dqvc_avi, flight_dist, flight_h2o,        &
28     cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv, &
29     qraindiag, qsnowdiag, dqreva, dqssub, dqrauto,     &
30     dqrcol, dqrmelt, dqrfreez, dqsauto, dqsagg, dqsrim,&
31     dqsmelt, dqsfreez)
32
33!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
34! Authors: Z.X. Li (LMD), J-L Dufresne (LMD), C. Rio (LMD), J-Y Grandpeix (LMD)
35!          A. JAM (LMD), J-B Madeleine (LMD), E. Vignon (LMD), L. Touzze-Peiffert (LMD)
36!--------------------------------------------------------------------------------
37! Date: 01/2021
38!--------------------------------------------------------------------------------
39! Aim: Large Scale Clouds and Precipitation (LSCP)
40!
41! This code is a new version of the fisrtilp.F90 routine, which is itself a
42! merge of 'first' (superrsaturation physics, P. LeVan K. Laval)
43! and 'ilp' (il pleut, L. Li)
44!
45! Compared to the original fisrtilp code, lscp
46! -> assumes thermcep = .TRUE. all the time (fisrtilp inconsistent when .FALSE.)
47! -> consider always precipitation thermalisation (fl_cor_ebil>0)
48! -> option iflag_fisrtilp_qsat<0 no longer possible (qsat does not evolve with T)
49! -> option "oldbug" by JYG has been removed
50! -> iflag_t_glace >0 always
51! -> the 'all or nothing' cloud approach is no longer available (cpartiel=T always)
52! -> rectangular distribution from L. Li no longer available
53! -> We always account for the Wegener-Findeisen-Bergeron process (iflag_bergeron = 2 in fisrt)
54!--------------------------------------------------------------------------------
55! References:
56!
57! - Bony, S., & Emanuel, K. A. 2001, JAS, doi: 10.1175/1520-0469(2001)058<3158:APOTCA>2.0.CO;2
58! - Hourdin et al. 2013, Clim Dyn, doi:10.1007/s00382-012-1343-y
59! - Jam et al. 2013, Boundary-Layer Meteorol, doi:10.1007/s10546-012-9789-3
60! - Jouhaud, et al. 2018. JAMES, doi:10.1029/2018MS001379
61! - Madeleine et al. 2020, JAMES, doi:10.1029/2020MS002046
62! - Touzze-Peifert Ludo, PhD thesis, p117-124
63! -------------------------------------------------------------------------------
64! Code structure:
65!
66! P0>     Thermalization of the precipitation coming from the overlying layer
67! P1>     Evaporation of the precipitation (falling from the k+1 level)
68! P2>     Cloud formation (at the k level)
69! P2.A.1> With the PDFs, calculation of cloud properties using the inital
70!         values of T and Q
71! P2.A.2> Coupling between condensed water and temperature
72! P2.A.3> Calculation of final quantities associated with cloud formation
73! P2.B>   Release of Latent heat after cloud formation
74! P3>     Autoconversion to precipitation (k-level)
75! P4>     Wet scavenging
76!------------------------------------------------------------------------------
77! Some preliminary comments (JBM) :
78!
79! The cloud water that the radiation scheme sees is not the same that the cloud
80! water used in the physics and the dynamics
81!
82! During the autoconversion to precipitation (P3 step), radocond (cloud water used
83! by the radiation scheme) is calculated as an average of the water that remains
84! in the cloud during the precipitation and not the water remaining at the end
85! of the time step. The latter is used in the rest of the physics and advected
86! by the dynamics.
87!
88! In summary:
89!
90! Radiation:
91! xflwc(newmicro)+xfiwc(newmicro) =
92!   radocond=lwcon(nc)+iwcon(nc)
93!
94! Notetheless, be aware of:
95!
96! radocond .NE. ocond(nc)
97! i.e.:
98! lwcon(nc)+iwcon(nc) .NE. ocond(nc)
99! but oliq+(ocond-oliq) .EQ. ocond
100! (which is not trivial)
101!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
102!
103
104! USE de modules contenant des fonctions.
105USE lmdz_cloudth, ONLY : cloudth, cloudth_v3, cloudth_v6, cloudth_mpc
106USE lmdz_lscp_tools, ONLY : calc_qsat_ecmwf, calc_gammasat
107USE lmdz_lscp_tools, ONLY : icefrac_lscp, icefrac_lscp_turb, distance_to_cloud_top
108USE lmdz_lscp_condensation, ONLY : condensation_lognormal, condensation_ice_supersat
109USE lmdz_lscp_precip, ONLY : histprecip_precld, histprecip_postcld
110USE lmdz_lscp_precip, ONLY : poprecip_precld, poprecip_postcld
111
112! Use du module lmdz_lscp_ini contenant les constantes
113USE lmdz_lscp_ini, ONLY : prt_level, lunout, eps
114USE lmdz_lscp_ini, ONLY : seuil_neb, iflag_evap_prec, DDT0
115USE lmdz_lscp_ini, ONLY : ok_radocond_snow, a_tr_sca
116USE lmdz_lscp_ini, ONLY : iflag_cloudth_vert, iflag_t_glace, iflag_fisrtilp_qsat
117USE lmdz_lscp_ini, ONLY : t_glace_min, min_frac_th_cld
118USE lmdz_lscp_ini, ONLY : RCPD, RLSTT, RLVTT, RVTMP2, RTT, RD, RG
119USE lmdz_lscp_ini, ONLY : ok_poprecip, ok_bug_phase_lscp
120USE lmdz_lscp_ini, ONLY : ok_ice_supersat, ok_unadjusted_clouds, iflag_icefrac
121
122IMPLICIT NONE
123
124!===============================================================================
125! VARIABLES DECLARATION
126!===============================================================================
127
128! INPUT VARIABLES:
129!-----------------
130
131  INTEGER,                         INTENT(IN)   :: klon,klev       ! number of horizontal grid points and vertical levels
132  REAL,                            INTENT(IN)   :: dtime           ! time step [s]
133  REAL,                            INTENT(IN)   :: missing_val     ! missing value for output
134
135  REAL, DIMENSION(klon,klev+1),    INTENT(IN)   :: paprs           ! inter-layer pressure [Pa]
136  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: pplay           ! mid-layer pressure [Pa]
137  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: temp            ! temperature (K)
138  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: omega           ! vertical velocity [Pa/s]
139  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: qt              ! total specific humidity (in vapor phase in input) [kg/kg]
140  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: ql_seri         ! liquid specific content [kg/kg]
141  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: qi_seri         ! ice specific content [kg/kg]
142  INTEGER,                         INTENT(IN)   :: iflag_cld_th    ! flag that determines the distribution of convective clouds
143  INTEGER,                         INTENT(IN)   :: iflag_ice_thermo! flag to activate the ice thermodynamics
144                                                                   ! CR: if iflag_ice_thermo=2, only convection is active
145  LOGICAL, DIMENSION(klon,klev),   INTENT(IN)   :: ptconv          ! grid points where deep convection scheme is active
146
147  !Inputs associated with thermal plumes
148
149  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: tv                  ! virtual potential temperature [K]
150  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: qta                 ! specific humidity within thermals [kg/kg]
151  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: fraca               ! fraction of thermals within the mesh [-]
152  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: pspsk               ! exner potential (p/100000)**(R/cp)
153  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: tla                 ! liquid temperature within thermals [K]
154  REAL, DIMENSION(klon,klev+1),    INTENT(IN)   :: tke                 !--turbulent kinetic energy [m2/s2]
155  REAL, DIMENSION(klon,klev+1),    INTENT(IN)   :: tke_dissip          !--TKE dissipation [m2/s3]
156
157  ! INPUT/OUTPUT variables
158  !------------------------
159 
160  REAL, DIMENSION(klon,klev),      INTENT(INOUT)   :: thl          ! liquid potential temperature [K]
161  REAL, DIMENSION(klon,klev),      INTENT(INOUT)   :: ratqs,sigma_qtherm        ! function of pressure that sets the large-scale
162
163
164  ! INPUT/OUTPUT condensation and ice supersaturation
165  !--------------------------------------------------
166  REAL, DIMENSION(klon,klev),      INTENT(INOUT):: cf_seri          ! cloud fraction [-]
167  REAL, DIMENSION(klon,klev),      INTENT(INOUT):: rvc_seri         ! cloudy water vapor to total water vapor ratio [-]
168  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: u_seri           ! eastward wind [m/s]
169  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: v_seri           ! northward wind [m/s]
170  REAL, DIMENSION(klon),           INTENT(IN)   :: cell_area        ! area of each cell [m2]
171
172  ! INPUT/OUTPUT aviation
173  !--------------------------------------------------
174  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: flight_dist      ! Aviation distance flown within the mesh [m/s/mesh]
175  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: flight_h2o       ! Aviation H2O emitted within the mesh [kg H2O/s/mesh]
176 
177  ! OUTPUT variables
178  !-----------------
179
180  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: d_t              ! temperature increment [K]
181  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: d_q              ! specific humidity increment [kg/kg]
182  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: d_ql             ! liquid water increment [kg/kg]
183  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: d_qi             ! cloud ice mass increment [kg/kg]
184  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: rneb             ! cloud fraction [-]
185  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: rneblsvol        ! cloud fraction per unit volume [-] 
186  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: pfraclr          ! precip fraction clear-sky part [-]
187  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: pfracld          ! precip fraction cloudy part [-]
188  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: cldfraliq           ! liquid fraction of cloud [-]
189  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: sigma2_icefracturb  ! Variance of the diagnostic supersaturation distribution (icefrac_turb) [-]
190  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: mean_icefracturb    ! Mean of the diagnostic supersaturation distribution (icefrac_turb) [-]
191  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: radocond         ! condensed water used in the radiation scheme [kg/kg]
192  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: radicefrac       ! ice fraction of condensed water for radiation scheme
193  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: rhcl             ! clear-sky relative humidity [-]
194  REAL, DIMENSION(klon),           INTENT(OUT)  :: rain             ! surface large-scale rainfall [kg/s/m2]
195  REAL, DIMENSION(klon),           INTENT(OUT)  :: snow             ! surface large-scale snowfall [kg/s/m2]
196  REAL, DIMENSION(klon,klev+1),    INTENT(OUT)  :: prfl             ! large-scale rainfall flux in the column [kg/s/m2]
197  REAL, DIMENSION(klon,klev+1),    INTENT(OUT)  :: psfl             ! large-scale snowfall flux in the column [kg/s/m2]
198  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: distcltop        ! distance to cloud top [m]
199  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: temp_cltop       ! temperature of cloud top [K]
200  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: beta             ! conversion rate of condensed water
201
202  ! fraction of aerosol scavenging through impaction and nucleation    (for on-line)
203 
204  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: frac_impa           ! scavenging fraction due tu impaction [-]
205  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: frac_nucl           ! scavenging fraction due tu nucleation [-]
206 
207  ! for condensation and ice supersaturation
208
209  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qsub           !--specific total water content in sub-saturated clear sky region [kg/kg]
210  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qissr          !--specific total water content in supersat region [kg/kg]
211  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qcld           !--specific total water content in cloudy region [kg/kg]
212  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: subfra         !--mesh fraction of subsaturated clear sky [-]   
213  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: issrfra        !--mesh fraction of ISSR [-] 
214  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: gamma_cond     !--coefficient governing the ice nucleation RHi threshold [-]     
215  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dcf_sub        !--cloud fraction tendency because of sublimation [s-1]
216  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dcf_con        !--cloud fraction tendency because of condensation [s-1]
217  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dcf_mix        !--cloud fraction tendency because of cloud mixing [s-1]
218  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqi_adj        !--specific ice content tendency because of temperature adjustment [kg/kg/s]
219  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqi_sub        !--specific ice content tendency because of sublimation [kg/kg/s]
220  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqi_con        !--specific ice content tendency because of condensation [kg/kg/s]
221  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqi_mix        !--specific ice content tendency because of cloud mixing [kg/kg/s]
222  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqvc_adj       !--specific cloud water vapor tendency because of temperature adjustment [kg/kg/s]
223  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqvc_sub       !--specific cloud water vapor tendency because of sublimation [kg/kg/s]
224  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqvc_con       !--specific cloud water vapor tendency because of condensation [kg/kg/s]
225  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqvc_mix       !--specific cloud water vapor tendency because of cloud mixing [kg/kg/s]
226  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qsatl          !--saturation specific humidity wrt liquid [kg/kg]
227  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qsati          !--saturation specific humidity wrt ice [kg/kg] 
228
229  ! for contrails and aviation
230
231  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: Tcontr         !--threshold temperature for contrail formation [K]
232  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qcontr         !--threshold humidity for contrail formation [kg/kg]
233  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qcontr2        !--// (2nd expression more consistent with LMDZ expression of q)
234  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: fcontrN        !--fraction of grid favourable to non-persistent contrails
235  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: fcontrP        !--fraction of grid favourable to persistent contrails
236  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dcf_avi        !--cloud fraction tendency because of aviation [s-1]
237  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqi_avi        !--specific ice content tendency because of aviation [kg/kg/s]
238  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqvc_avi       !--specific cloud water vapor tendency because of aviation [kg/kg/s]
239
240
241  ! for POPRECIP
242
243  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qraindiag      !--DIAGNOSTIC specific rain content [kg/kg]
244  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qsnowdiag      !--DIAGNOSTIC specific snow content [kg/kg]
245  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqreva         !--rain tendendy due to evaporation [kg/kg/s]
246  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqssub         !--snow tendency due to sublimation [kg/kg/s]
247  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqrcol         !--rain tendendy due to collection by rain of liquid cloud droplets [kg/kg/s]
248  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqsagg         !--snow tendency due to collection of lcoud ice by aggregation [kg/kg/s]
249  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqrauto        !--rain tendency due to autoconversion of cloud liquid [kg/kg/s]
250  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqsauto        !--snow tendency due to autoconversion of cloud ice [kg/kg/s]
251  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqsrim         !--snow tendency due to riming [kg/kg/s]
252  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqsmelt        !--snow tendency due to melting [kg/kg/s]
253  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqrmelt        !--rain tendency due to melting [kg/kg/s]
254  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqsfreez       !--snow tendency due to freezing [kg/kg/s]
255  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqrfreez       !--rain tendency due to freezing [kg/kg/s]
256
257  ! for thermals
258
259  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: cloudth_sth      !--mean saturation deficit in thermals
260  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: cloudth_senv     !--mean saturation deficit in environment
261  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: cloudth_sigmath  !--std of saturation deficit in thermals
262  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: cloudth_sigmaenv !--std of saturation deficit in environment
263
264
265  ! LOCAL VARIABLES:
266  !----------------
267  REAL,DIMENSION(klon) :: qice_ini
268  REAL, DIMENSION(klon,klev) :: ctot
269  REAL, DIMENSION(klon,klev) :: ctot_vol
270  REAL, DIMENSION(klon) :: zqs, zdqs
271  REAL :: zdelta
272  REAL, DIMENSION(klon) :: zdqsdT_raw
273  REAL, DIMENSION(klon) :: gammasat,dgammasatdt                   ! coefficient to make cold condensation at the correct RH and derivative wrt T
274  REAL, DIMENSION(klon) :: Tbef,qlbef,DT                          ! temperature, humidity and temp. variation during lognormal iteration
275  REAL :: num,denom
276  REAL :: cste
277  REAL, DIMENSION(klon) :: zfice_th
278  REAL, DIMENSION(klon) :: qcloud, qincloud_mpc
279  REAL, DIMENSION(klon) :: zrfl, zifl
280  REAL, DIMENSION(klon) :: zoliq, zcond, zq, zqn
281  REAL, DIMENSION(klon) :: zoliql, zoliqi
282  REAL, DIMENSION(klon) :: zt
283  REAL, DIMENSION(klon) :: zfice,zneb
284  REAL, DIMENSION(klon) :: dzfice
285  REAL, DIMENSION(klon) :: zfice_turb, dzfice_turb
286  REAL, DIMENSION(klon) :: qtot, qzero
287  ! Variables precipitation energy conservation
288  REAL, DIMENSION(klon) :: zmqc
289  REAL :: zalpha_tr
290  REAL :: zfrac_lessi
291  REAL, DIMENSION(klon) :: zprec_cond
292  REAL, DIMENSION(klon) :: zlh_solid
293  REAL, DIMENSION(klon) :: ztupnew
294  REAL, DIMENSION(klon) :: zqvapclr, zqupnew ! for poprecip evap / subl
295  REAL, DIMENSION(klon) :: zrflclr, zrflcld
296  REAL, DIMENSION(klon) :: ziflclr, ziflcld
297  REAL, DIMENSION(klon) :: znebprecip, znebprecipclr, znebprecipcld
298  REAL, DIMENSION(klon) :: tot_zneb
299  REAL :: qlmpc, qimpc, rnebmpc
300  REAL, DIMENSION(klon) :: zdistcltop, ztemp_cltop
301  REAL, DIMENSION(klon) :: zqliq, zqice, zqvapcl        ! for icefrac_lscp_turb
302
303  ! for quantity of condensates seen by radiation
304  REAL, DIMENSION(klon) :: zradocond, zradoice
305  REAL, DIMENSION(klon) :: zrho_up, zvelo_up
306 
307  ! for condensation and ice supersaturation
308  REAL, DIMENSION(klon) :: qvc, shear
309  REAL :: delta_z
310  !--Added for ice supersaturation (ok_ice_supersat) and contrails (ok_plane_contrails)
311  ! Constants used for calculating ratios that are advected (using a parent-child
312  ! formalism). This is not done in the dynamical core because at this moment,
313  ! only isotopes can use this parent-child formalism. Note that the two constants
314  ! are the same as the one use in the dynamical core, being also defined in
315  ! dyn3d_common/infotrac.F90
316  REAL :: min_qParent, min_ratio
317
318  INTEGER i, k, kk, iter
319  INTEGER, DIMENSION(klon) :: n_i
320  INTEGER ncoreczq
321  INTEGER, DIMENSION(klon,klev) :: mpc_bl_points
322  LOGICAL iftop
323
324  LOGICAL, DIMENSION(klon) :: lognormale
325  LOGICAL, DIMENSION(klon) :: keepgoing
326
327  CHARACTER (len = 20) :: modname = 'lscp'
328  CHARACTER (len = 80) :: abort_message
329 
330
331!===============================================================================
332! INITIALISATION
333!===============================================================================
334
335! Few initial checks
336
337
338IF (iflag_fisrtilp_qsat .LT. 0) THEN
339    abort_message = 'lscp cannot be used with iflag_fisrtilp<0'
340    CALL abort_physic(modname,abort_message,1)
341ENDIF
342
343! Few initialisations
344
345ctot_vol(1:klon,1:klev)=0.0
346rneblsvol(1:klon,1:klev)=0.0
347znebprecip(:)=0.0
348znebprecipclr(:)=0.0
349znebprecipcld(:)=0.0
350mpc_bl_points(:,:)=0
351
352IF (prt_level>9) WRITE(lunout,*) 'NUAGES4 A. JAM'
353
354! AA for 'safety' reasons
355zalpha_tr   = 0.
356zfrac_lessi = 0.
357beta(:,:)= 0.
358
359! Initialisation of variables:
360
361prfl(:,:) = 0.0
362psfl(:,:) = 0.0
363d_t(:,:) = 0.0
364d_q(:,:) = 0.0
365d_ql(:,:) = 0.0
366d_qi(:,:) = 0.0
367rneb(:,:) = 0.0
368pfraclr(:,:)=0.0
369pfracld(:,:)=0.0
370cldfraliq(:,:)=0.
371sigma2_icefracturb(:,:)=0.
372mean_icefracturb(:,:)=0.
373radocond(:,:) = 0.0
374radicefrac(:,:) = 0.0
375frac_nucl(:,:) = 1.0
376frac_impa(:,:) = 1.0
377rain(:) = 0.0
378snow(:) = 0.0
379zfice(:)=0.0
380dzfice(:)=0.0
381zfice_turb(:)=0.0
382dzfice_turb(:)=0.0
383zrfl(:) = 0.0
384zifl(:) = 0.0
385zneb(:) = seuil_neb
386zrflclr(:) = 0.0
387ziflclr(:) = 0.0
388zrflcld(:) = 0.0
389ziflcld(:) = 0.0
390tot_zneb(:) = 0.0
391qzero(:) = 0.0
392zdistcltop(:)=0.0
393ztemp_cltop(:) = 0.0
394ztupnew(:)=0.0
395
396distcltop(:,:)=0.
397temp_cltop(:,:)=0.
398
399!--Ice supersaturation
400gamma_cond(:,:) = 1.
401qissr(:,:)      = 0.
402issrfra(:,:)    = 0.
403dcf_sub(:,:)    = 0.
404dcf_con(:,:)    = 0.
405dcf_mix(:,:)    = 0.
406dqi_adj(:,:)    = 0.
407dqi_sub(:,:)    = 0.
408dqi_con(:,:)    = 0.
409dqi_mix(:,:)    = 0.
410dqvc_adj(:,:)   = 0.
411dqvc_sub(:,:)   = 0.
412dqvc_con(:,:)   = 0.
413dqvc_mix(:,:)   = 0.
414fcontrN(:,:)    = 0.
415fcontrP(:,:)    = 0.
416Tcontr(:,:)     = missing_val
417qcontr(:,:)     = missing_val
418qcontr2(:,:)    = missing_val
419dcf_avi(:,:)    = 0.
420dqi_avi(:,:)    = 0.
421dqvc_avi(:,:)   = 0.
422qvc(:)          = 0.
423shear(:)        = 0.
424min_qParent     = 1.e-30
425min_ratio       = 1.e-16
426
427!-- poprecip
428qraindiag(:,:)= 0.
429qsnowdiag(:,:)= 0.
430dqreva(:,:)   = 0.
431dqrauto(:,:)  = 0.
432dqrmelt(:,:)  = 0.
433dqrfreez(:,:) = 0.
434dqrcol(:,:)   = 0.
435dqssub(:,:)   = 0.
436dqsauto(:,:)  = 0.
437dqsrim(:,:)   = 0.
438dqsagg(:,:)   = 0.
439dqsfreez(:,:) = 0.
440dqsmelt(:,:)  = 0.
441zqupnew(:)    = 0.
442zqvapclr(:)   = 0.
443
444
445
446!c_iso: variable initialisation for iso
447
448
449!===============================================================================
450!              BEGINNING OF VERTICAL LOOP FROM TOP TO BOTTOM
451!===============================================================================
452
453ncoreczq=0
454
455DO k = klev, 1, -1
456
457    qice_ini = qi_seri(:,k)
458
459    IF (k.LE.klev-1) THEN
460       iftop=.false.
461    ELSE
462       iftop=.true.
463    ENDIF
464
465    ! Initialisation temperature and specific humidity
466    ! temp(klon,klev) is not modified by the routine, instead all changes in temperature are made on zt
467    ! at the end of the klon loop, a temperature incremtent d_t due to all processes
468    ! (thermalization, evap/sub incoming precip, cloud formation, precipitation processes) is calculated
469    ! d_t = temperature tendency due to lscp
470    ! The temperature of the overlying layer is updated here because needed for thermalization
471    DO i = 1, klon
472        zt(i)=temp(i,k)
473        zq(i)=qt(i,k)
474        IF (.not. iftop) THEN
475           ztupnew(i)  = temp(i,k+1) + d_t(i,k+1)
476           zqupnew(i)  = qt(i,k+1) + d_q(i,k+1) + d_ql(i,k+1) + d_qi(i,k+1)
477           !--zqs(i) is the saturation specific humidity in the layer above
478           zqvapclr(i) = MAX(0., qt(i,k+1) + d_q(i,k+1) - rneb(i,k+1) * zqs(i))
479        ENDIF
480        !c_iso init of iso
481    ENDDO
482
483    ! --------------------------------------------------------------------
484    ! P1> Precipitation processes, before cloud formation:
485    !     Thermalization of precipitation falling from the overlying layer AND
486    !     Precipitation evaporation/sublimation/melting
487    !---------------------------------------------------------------------
488   
489    !================================================================
490    ! Flag for the new and more microphysical treatment of precipitation from Atelier Nuage (R)
491    IF ( ok_poprecip .OR. ( iflag_evap_prec .EQ. 6 ) ) THEN
492
493      CALL poprecip_precld(klon, dtime, iftop, paprs(:,k), paprs(:,k+1), pplay(:,k), &
494                        zt, ztupnew, zq, zmqc, znebprecipclr, znebprecipcld, &
495                        zqvapclr, zqupnew, &
496                        zrfl, zrflclr, zrflcld, &
497                        zifl, ziflclr, ziflcld, &
498                        dqreva(:,k), dqssub(:,k) &
499                        )
500
501    !================================================================
502    ELSE
503
504      CALL histprecip_precld(klon, dtime, iftop, paprs(:,k), paprs(:,k+1), pplay(:,k), &
505                        zt, ztupnew, zq, zmqc, zneb, znebprecip, znebprecipclr, &
506                        zrfl, zrflclr, zrflcld, &
507                        zifl, ziflclr, ziflcld, &
508                        dqreva(:,k), dqssub(:,k) &
509                        )
510
511    ENDIF ! (ok_poprecip)
512   
513    ! Calculation of qsat, L/Cp*dqsat/dT and ncoreczq counter
514    !-------------------------------------------------------
515
516     qtot(:)=zq(:)+zmqc(:)
517     CALL calc_qsat_ecmwf(klon,zt,qtot,pplay(:,k),RTT,0,.false.,zqs,zdqs)
518     DO i = 1, klon
519            zdelta = MAX(0.,SIGN(1.,RTT-zt(i)))
520            zdqsdT_raw(i) = zdqs(i)*RCPD*(1.0+RVTMP2*zq(i)) / (RLVTT*(1.-zdelta) + RLSTT*zdelta)
521            IF (zq(i) .LT. 1.e-15) THEN
522                ncoreczq=ncoreczq+1
523                zq(i)=1.e-15
524            ENDIF
525            ! c_iso: do something similar for isotopes
526
527     ENDDO
528   
529    ! --------------------------------------------------------------------
530    ! P2> Cloud formation
531    !---------------------------------------------------------------------
532    !
533    ! Unlike fisrtilp, we always assume a 'fractional cloud' approach
534    ! i.e. clouds occupy only a fraction of the mesh (the subgrid distribution
535    ! is prescribed and depends on large scale variables and boundary layer
536    ! properties)
537    ! The decrease in condensed part due tu latent heating is taken into
538    ! account
539    ! -------------------------------------------------------------------
540
541        ! P2.1> With the PDFs (log-normal, bigaussian)
542        ! cloud properties calculation with the initial values of t and q
543        ! ----------------------------------------------------------------
544
545        ! initialise gammasat and qincloud_mpc
546        gammasat(:)=1.
547        qincloud_mpc(:)=0.
548
549        IF (iflag_cld_th.GE.5) THEN
550               ! Cloud cover and content in meshes affected by shallow convection,
551               ! are retrieved from a bi-gaussian distribution of the saturation deficit
552               ! following Jam et al. 2013
553
554               IF (iflag_cloudth_vert.LE.2) THEN
555                  ! Old version of Arnaud Jam
556
557                    CALL cloudth(klon,klev,k,tv,                  &
558                         zq,qta,fraca,                            &
559                         qcloud,ctot,pspsk,paprs,pplay,tla,thl, &
560                         ratqs,zqs,temp,                              &
561                         cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
562
563
564                ELSEIF (iflag_cloudth_vert.GE.3 .AND. iflag_cloudth_vert.LE.5) THEN
565                   ! Default version of Arnaud Jam
566                       
567                    CALL cloudth_v3(klon,klev,k,tv,                        &
568                         zq,qta,fraca,                                     &
569                         qcloud,ctot,ctot_vol,pspsk,paprs,pplay,tla,thl, &
570                         ratqs,sigma_qtherm,zqs,temp, &
571                         cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
572
573
574                ELSEIF (iflag_cloudth_vert.EQ.6) THEN
575                   ! Jean Jouhaud's version, with specific separation between surface and volume
576                   ! cloud fraction Decembre 2018
577
578                    CALL cloudth_v6(klon,klev,k,tv,                        &
579                         zq,qta,fraca,                                     &
580                         qcloud,ctot,ctot_vol,pspsk,paprs,pplay,tla,thl, &
581                         ratqs,zqs,temp, &
582                         cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
583
584                ELSEIF (iflag_cloudth_vert .EQ. 7) THEN
585                   ! Updated version of Arnaud Jam (correction by E. Vignon) + adapted treatment
586                   ! for boundary-layer mixed phase clouds
587                    CALL cloudth_mpc(klon,klev,k,mpc_bl_points,zt,zq,qta(:,k),fraca(:,k), &
588                                     pspsk(:,k),paprs(:,k+1),paprs(:,k),pplay(:,k), tla(:,k), &
589                                     ratqs(:,k),qcloud,qincloud_mpc,zfice_th,ctot(:,k),ctot_vol(:,k), &
590                                     cloudth_sth(:,k),cloudth_senv(:,k),cloudth_sigmath(:,k),cloudth_sigmaenv(:,k))
591
592                ENDIF
593
594
595                DO i=1,klon
596                    rneb(i,k)=ctot(i,k)
597                    rneblsvol(i,k)=ctot_vol(i,k)
598                    zqn(i)=qcloud(i)
599                    !--AB grid-mean vapor in the cloud - we assume saturation adjustment
600                    qvc(i) = rneb(i,k) * zqs(i)
601                ENDDO
602
603        ENDIF
604
605        IF (iflag_cld_th .LE. 4) THEN
606           
607                ! lognormal
608            lognormale(:) = .TRUE.
609
610        ELSEIF (iflag_cld_th .GE. 6) THEN
611           
612                ! lognormal distribution when no thermals
613            lognormale(:) = fraca(:,k) < min_frac_th_cld
614
615        ELSE
616                ! When iflag_cld_th=5, we always assume
617                ! bi-gaussian distribution
618            lognormale(:) = .FALSE.
619       
620        ENDIF
621
622        DT(:) = 0.
623        n_i(:)=0
624        Tbef(:)=zt(:)
625        qlbef(:)=0.
626
627        ! Treatment of non-boundary layer clouds (lognormale)
628        ! We iterate here to take into account the change in qsat(T) and ice fraction
629        ! during the condensation process
630        ! the increment in temperature is calculated using the first principle of
631        ! thermodynamics (enthalpy conservation equation in a mesh composed of a cloud fraction
632        ! and a clear sky fraction)
633        ! note that we assume that the vapor in the cloud is at saturation for this calculation     
634
635        DO iter=1,iflag_fisrtilp_qsat+1
636
637                keepgoing(:) = .FALSE.
638
639                DO i=1,klon
640
641                    ! keepgoing = .true. while convergence is not satisfied
642
643                    IF (((ABS(DT(i)).GT.DDT0) .OR. (n_i(i) .EQ. 0)) .AND. lognormale(i)) THEN
644
645                        ! if not convergence:
646                        ! we calculate a new iteration
647                        keepgoing(i) = .TRUE.
648
649                        ! P2.2.1> cloud fraction and condensed water mass calculation
650                        ! Calculated variables:
651                        ! rneb : cloud fraction
652                        ! zqn : total water within the cloud
653                        ! zcond: mean condensed water within the mesh
654                        ! rhcl: clear-sky relative humidity
655                        !---------------------------------------------------------------
656
657                        ! new temperature that only serves in the iteration process:
658                        Tbef(i)=Tbef(i)+DT(i)
659
660                        ! Rneb, qzn and zcond for lognormal PDFs
661                        qtot(i)=zq(i)+zmqc(i)
662
663                      ENDIF
664
665                  ENDDO
666
667                  ! Calculation of saturation specific humidity and ice fraction
668                  CALL calc_qsat_ecmwf(klon,Tbef,qtot,pplay(:,k),RTT,0,.false.,zqs,zdqs)
669                  CALL calc_gammasat(klon,Tbef,qtot,pplay(:,k),gammasat,dgammasatdt)
670                  ! saturation may occur at a humidity different from qsat (gamma qsat), so gamma correction for dqs
671                  zdqs(:) = gammasat(:)*zdqs(:)+zqs(:)*dgammasatdt(:)
672                  ! cloud phase determination
673                  IF (iflag_t_glace.GE.4) THEN
674                  ! For iflag_t_glace GE 4 the phase partition function dependends on temperature AND distance to cloud top
675                       CALL distance_to_cloud_top(klon,klev,k,temp,pplay,paprs,rneb,zdistcltop,ztemp_cltop)
676                  ENDIF
677
678                  CALL icefrac_lscp(klon, zt, iflag_ice_thermo, zdistcltop,ztemp_cltop,zfice,dzfice)
679
680                  !--AB Activates a condensation scheme that allows for
681                  !--ice supersaturation and contrails evolution from aviation
682                  IF (ok_ice_supersat) THEN
683
684                    !--Calculate the shear value (input for condensation and ice supersat)
685                    DO i = 1, klon
686                      !--Cell thickness [m]
687                      delta_z = ( paprs(i,k) - paprs(i,k+1) ) / RG / pplay(i,k) * Tbef(i) * RD
688                      IF ( iftop ) THEN
689                        ! top
690                        shear(i) = SQRT( ( (u_seri(i,k) - u_seri(i,k-1)) / delta_z )**2. &
691                                       + ( (v_seri(i,k) - v_seri(i,k-1)) / delta_z )**2. )
692                      ELSEIF ( k .EQ. 1 ) THEN
693                        ! surface
694                        shear(i) = SQRT( ( (u_seri(i,k+1) - u_seri(i,k)) / delta_z )**2. &
695                                       + ( (v_seri(i,k+1) - v_seri(i,k)) / delta_z )**2. )
696                      ELSE
697                        ! other layers
698                        shear(i) = SQRT( ( ( (u_seri(i,k+1) + u_seri(i,k)) / 2. &
699                                           - (u_seri(i,k) + u_seri(i,k-1)) / 2. ) / delta_z )**2. &
700                                       + ( ( (v_seri(i,k+1) + v_seri(i,k)) / 2. &
701                                           - (v_seri(i,k) + v_seri(i,k-1)) / 2. ) / delta_z )**2. )
702                      ENDIF
703                    ENDDO
704
705                    !---------------------------------------------
706                    !--   CONDENSATION AND ICE SUPERSATURATION  --
707                    !---------------------------------------------
708
709                    CALL condensation_ice_supersat( &
710                        klon, dtime, missing_val, &
711                        pplay(:,k), paprs(:,k), paprs(:,k+1), &
712                        cf_seri(:,k), rvc_seri(:,k), ql_seri(:,k), qi_seri(:,k), &
713                        shear, tke_dissip(:,k), cell_area, &
714                        Tbef, zq, zqs, gammasat, ratqs(:,k), keepgoing, &
715                        rneb(:,k), zqn, qvc, issrfra(:,k), qissr(:,k), &
716                        dcf_sub(:,k), dcf_con(:,k), dcf_mix(:,k), &
717                        dqi_adj(:,k), dqi_sub(:,k), dqi_con(:,k), dqi_mix(:,k), &
718                        dqvc_adj(:,k), dqvc_sub(:,k), dqvc_con(:,k), dqvc_mix(:,k), &
719                        Tcontr(:,k), qcontr(:,k), qcontr2(:,k), fcontrN(:,k), fcontrP(:,k), &
720                        flight_dist(:,k), flight_h2o(:,k), &
721                        dcf_avi(:,k), dqi_avi(:,k), dqvc_avi(:,k))
722
723
724                  ELSE
725                  !--generalised lognormal condensation scheme (Bony and Emanuel 2001)
726
727                   CALL condensation_lognormal( &
728                       klon, Tbef, zq, zqs, gammasat, ratqs(:,k), &
729                       keepgoing, rneb(:,k), zqn, qvc)
730
731
732                  ENDIF ! .NOT. ok_ice_supersat
733
734                  DO i=1,klon
735                      IF (keepgoing(i)) THEN
736
737                        ! If vertical heterogeneity, change fraction by volume as well
738                        IF (iflag_cloudth_vert.GE.3) THEN
739                            ctot_vol(i,k)=rneb(i,k)
740                            rneblsvol(i,k)=ctot_vol(i,k)
741                        ENDIF
742
743
744                       ! P2.2.2> Approximative calculation of temperature variation DT
745                       !           due to condensation.
746                       ! Calculated variables:
747                       ! dT : temperature change due to condensation
748                       !---------------------------------------------------------------
749
750               
751                        IF (zfice(i).LT.1) THEN
752                            cste=RLVTT
753                        ELSE
754                            cste=RLSTT
755                        ENDIF
756                       
757                        ! LEA_R : check formule
758                        IF ( ok_unadjusted_clouds ) THEN
759                          !--AB We relax the saturation adjustment assumption
760                          !-- qvc (grid-mean vapor in cloud) is calculated by the condensation scheme
761                          IF ( rneb(i,k) .GT. eps ) THEN
762                            qlbef(i) = MAX(0., zqn(i) - qvc(i) / rneb(i,k))
763                          ELSE
764                            qlbef(i) = 0.
765                          ENDIF
766                        ELSE
767                          qlbef(i)=max(0.,zqn(i)-zqs(i))
768                        ENDIF
769
770                        num = -Tbef(i)+zt(i)+rneb(i,k)*((1-zfice(i))*RLVTT &
771                              +zfice(i)*RLSTT)/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))*qlbef(i)
772                        denom = 1.+rneb(i,k)*((1-zfice(i))*RLVTT+zfice(i)*RLSTT)/cste*zdqs(i) &
773                              -(RLSTT-RLVTT)/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))*rneb(i,k)      &
774                              *qlbef(i)*dzfice(i)
775                        ! here we update a provisory temperature variable that only serves in the iteration
776                        ! process
777                        DT(i)=num/denom
778                        n_i(i)=n_i(i)+1
779
780                    ENDIF ! end keepgoing
781
782                ENDDO     ! end loop on i
783
784        ENDDO         ! iter=1,iflag_fisrtilp_qsat+1
785
786        ! P2.2> Final quantities calculation
787        ! Calculated variables:
788        !   rneb : cloud fraction
789        !   zcond: mean condensed water in the mesh
790        !   zqn  : mean water vapor in the mesh
791        !   zfice: ice fraction in clouds
792        !   zt   : temperature
793        !   rhcl : clear-sky relative humidity
794        ! ----------------------------------------------------------------
795
796
797        ! For iflag_t_glace GE 4 the phase partition function dependends on temperature AND distance to cloud top
798        IF (iflag_t_glace.GE.4) THEN
799           CALL distance_to_cloud_top(klon,klev,k,temp,pplay,paprs,rneb,zdistcltop,ztemp_cltop)
800           distcltop(:,k)=zdistcltop(:)
801           temp_cltop(:,k)=ztemp_cltop(:)
802        ENDIF
803
804        ! Partition function depending on temperature
805        CALL icefrac_lscp(klon, zt, iflag_ice_thermo, zdistcltop, ztemp_cltop, zfice, dzfice)
806
807        ! Partition function depending on tke for non shallow-convective clouds
808        IF (iflag_icefrac .GE. 1) THEN
809           CALL icefrac_lscp_turb(klon, dtime, Tbef, pplay(:,k), paprs(:,k), paprs(:,k+1), omega(:,k), qice_ini, ziflcld, zqn, &
810           rneb(:,k), tke(:,k), tke_dissip(:,k), zqliq, zqvapcl, zqice, zfice_turb, dzfice_turb, cldfraliq(:,k),sigma2_icefracturb(:,k), mean_icefracturb(:,k))
811        ENDIF
812
813        ! Water vapor update, Phase determination and subsequent latent heat exchange
814        DO i=1, klon
815            ! Overwrite phase partitioning in boundary layer mixed phase clouds when the
816            ! iflag_cloudth_vert=7 and specific param is activated
817            IF (mpc_bl_points(i,k) .GT. 0) THEN
818                zcond(i) = MAX(0.0,qincloud_mpc(i))*rneb(i,k)
819                ! following line is very strange and probably wrong
820                rhcl(i,k)= (zqs(i)+zq(i))/2./zqs(i)
821                ! water vapor update and partition function if thermals
822                zq(i) = zq(i) - zcond(i)       
823                zfice(i)=zfice_th(i)
824            ELSE
825                ! Checks on rneb, rhcl and zqn
826                IF (rneb(i,k) .LE. 0.0) THEN
827                    zqn(i) = 0.0
828                    rneb(i,k) = 0.0
829                    zcond(i) = 0.0
830                    rhcl(i,k)=zq(i)/zqs(i)
831                ELSE IF (rneb(i,k) .GE. 1.0) THEN
832                    zqn(i) = zq(i)
833                    rneb(i,k) = 1.0
834                    IF ( ok_unadjusted_clouds ) THEN
835                      !--AB We relax the saturation adjustment assumption
836                      !-- qvc (grid-mean vapor in cloud) is calculated by the condensation scheme
837                      zcond(i) = MAX(0., zqn(i) - qvc(i))
838                    ELSE
839                      zcond(i) = MAX(0.0,zqn(i)-zqs(i))
840                    ENDIF
841                    rhcl(i,k)=1.0
842                ELSE
843                    IF ( ok_unadjusted_clouds ) THEN
844                      !--AB We relax the saturation adjustment assumption
845                      !-- qvc (grid-mean vapor in cloud) is calculated by the condensation scheme
846                      zcond(i) = MAX(0., zqn(i) * rneb(i,k) - qvc(i))
847                    ELSE
848                      zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)
849                    ENDIF
850                    ! following line is very strange and probably wrong:
851                    rhcl(i,k)=(zqs(i)+zq(i))/2./zqs(i)
852                    ! Overwrite partitioning for non shallow-convective clouds if iflag_icefrac>1 (icefrac turb param)
853                    IF (iflag_icefrac .GE. 1) THEN
854                        IF (lognormale(i)) THEN 
855                           zcond(i)  = zqliq(i) + zqice(i)
856                           zfice(i)=zfice_turb(i)
857                           rhcl(i,k) = zqvapcl(i) * rneb(i,k) + (zq(i) - zqn(i)) * (1.-rneb(i,k))
858                        ENDIF
859                    ENDIF
860                ENDIF
861
862                ! water vapor update
863                zq(i) = zq(i) - zcond(i)
864
865            ENDIF
866
867            ! c_iso : routine that computes in-cloud supersaturation
868            ! c_iso condensation of isotopes (zcond, zsursat, zfice, zq in input)
869
870            ! temperature update due to phase change
871            zt(i) = zt(i) + (1.-zfice(i))*zcond(i) &
872            &     * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i))) &
873                  +zfice(i)*zcond(i) * RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i)))
874        ENDDO
875
876        ! If vertical heterogeneity, change volume fraction
877        IF (iflag_cloudth_vert .GE. 3) THEN
878          ctot_vol(1:klon,k)=min(max(ctot_vol(1:klon,k),0.),1.)
879          rneblsvol(1:klon,k)=ctot_vol(1:klon,k)
880        ENDIF
881
882
883    ! ----------------------------------------------------------------
884    ! P3> Precipitation processes, after cloud formation
885    !     - precipitation formation, melting/freezing
886    ! ----------------------------------------------------------------
887
888    ! Initiate the quantity of liquid and solid condensates
889    ! Note that in the following, zcond is the total amount of condensates
890    ! including newly formed precipitations (i.e., condensates formed by the
891    ! condensation process), while zoliq is the total amount of condensates in
892    ! the cloud (i.e., on which precipitation processes have applied)
893    DO i = 1, klon
894      zoliq(i)  = zcond(i)
895      zoliql(i) = zcond(i) * ( 1. - zfice(i) )
896      zoliqi(i) = zcond(i) * zfice(i)
897      ! c_iso : initialisation of zoliq* also for isotopes
898    ENDDO
899
900    !================================================================
901    ! Flag for the new and more microphysical treatment of precipitation from Atelier Nuage (R)
902    IF (ok_poprecip) THEN
903
904      CALL poprecip_postcld(klon, dtime, paprs(:,k), paprs(:,k+1), pplay(:,k), &
905                            ctot_vol(:,k), ptconv(:,k), &
906                            zt, zq, zoliql, zoliqi, zfice, &
907                            rneb(:,k), znebprecipclr, znebprecipcld, &
908                            zrfl, zrflclr, zrflcld, &
909                            zifl, ziflclr, ziflcld, &
910                            qraindiag(:,k), qsnowdiag(:,k), dqrauto(:,k), &
911                            dqrcol(:,k), dqrmelt(:,k), dqrfreez(:,k), &
912                            dqsauto(:,k), dqsagg(:,k), dqsrim(:,k), &
913                            dqsmelt(:,k), dqsfreez(:,k) &
914                            )
915
916    !================================================================
917    ELSE
918
919      CALL histprecip_postcld(klon, dtime, iftop, paprs(:,k), paprs(:,k+1), pplay(:,k), &
920                            ctot_vol(:,k), ptconv(:,k), zdqsdT_raw, &
921                            zt, zq, zoliq, zoliql, zoliqi, zcond, zfice, zmqc, &
922                            rneb(:,k), znebprecipclr, znebprecipcld, &
923                            zneb, tot_zneb, zrho_up, zvelo_up, &
924                            zrfl, zrflclr, zrflcld, zifl, ziflclr, ziflcld, &
925                            zradocond, zradoice, dqrauto(:,k), dqsauto(:,k) &
926                            )
927
928    ENDIF ! ok_poprecip
929
930    ! End of precipitation processes after cloud formation
931    ! ----------------------------------------------------
932
933    !----------------------------------------------------------------------
934    ! P4> Calculation of cloud condensates amount seen by radiative scheme
935    !----------------------------------------------------------------------
936
937    DO i=1,klon
938
939      IF (ok_poprecip) THEN
940        IF (ok_radocond_snow) THEN
941          radocond(i,k) = zoliq(i)
942          zradoice(i) = zoliqi(i) + qsnowdiag(i,k)
943        ELSE
944          radocond(i,k) = zoliq(i)
945          zradoice(i) = zoliqi(i)
946        ENDIF
947      ELSE
948        radocond(i,k) = zradocond(i)
949      ENDIF
950
951      ! calculate the percentage of ice in "radocond" so cloud+precip seen
952      ! by the radiation scheme
953      IF (radocond(i,k) .GT. 0.) THEN
954        radicefrac(i,k)=MIN(MAX(zradoice(i)/radocond(i,k),0.),1.)
955      ENDIF
956    ENDDO
957
958    ! ----------------------------------------------------------------
959    ! P5> Wet scavenging
960    ! ----------------------------------------------------------------
961
962    !Scavenging through nucleation in the layer
963
964    DO i = 1,klon
965       
966        IF(zcond(i).GT.zoliq(i)+1.e-10) THEN
967            beta(i,k) = (zcond(i)-zoliq(i))/zcond(i)/dtime
968        ELSE
969            beta(i,k) = 0.
970        ENDIF
971
972        zprec_cond(i) = MAX(zcond(i)-zoliq(i),0.0)*(paprs(i,k)-paprs(i,k+1))/RG
973
974        IF (rneb(i,k).GT.0.0.AND.zprec_cond(i).GT.0.) THEN
975
976            IF (temp(i,k) .GE. t_glace_min) THEN
977                zalpha_tr = a_tr_sca(3)
978            ELSE
979                zalpha_tr = a_tr_sca(4)
980            ENDIF
981
982            zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i))
983            frac_nucl(i,k)= 1.-zneb(i)*zfrac_lessi
984
985            ! Nucleation with a factor of -1 instead of -0.5
986            zfrac_lessi = 1. - EXP(-zprec_cond(i)/zneb(i))
987
988        ENDIF
989
990    ENDDO
991
992    ! Scavenging through impaction in the underlying layer
993
994    DO kk = k-1, 1, -1
995
996        DO i = 1, klon
997
998            IF (rneb(i,k).GT.0.0.AND.zprec_cond(i).GT.0.) THEN
999
1000                IF (temp(i,kk) .GE. t_glace_min) THEN
1001                    zalpha_tr = a_tr_sca(1)
1002                ELSE
1003                    zalpha_tr = a_tr_sca(2)
1004                ENDIF
1005
1006              zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i))
1007              frac_impa(i,kk)= 1.-zneb(i)*zfrac_lessi
1008
1009            ENDIF
1010
1011        ENDDO
1012
1013    ENDDO
1014   
1015    !------------------------------------------------------------
1016    ! P6 > write diagnostics and outputs
1017    !------------------------------------------------------------
1018   
1019    !--AB Write diagnostics and tracers for ice supersaturation
1020    IF ( ok_ice_supersat ) THEN
1021      CALL calc_qsat_ecmwf(klon,zt,qzero,pplay(:,k),RTT,1,.false.,qsatl(:,k),zdqs)
1022      CALL calc_qsat_ecmwf(klon,zt,qzero,pplay(:,k),RTT,2,.false.,qsati(:,k),zdqs)
1023
1024      DO i = 1, klon
1025
1026        IF ( zoliq(i) .LE. 0. ) THEN
1027          !--If everything was precipitated, the remaining empty cloud is dissipated
1028          !--and everything is transfered to the subsaturated clear sky region
1029          rneb(i,k) = 0.
1030          qvc(i) = 0.
1031        ENDIF
1032
1033        cf_seri(i,k) = rneb(i,k)
1034
1035        IF ( .NOT. ok_unadjusted_clouds ) THEN
1036          qvc(i) = zqs(i) * rneb(i,k)
1037        ENDIF
1038        IF ( zq(i) .GT. min_qParent ) THEN
1039          rvc_seri(i,k) = qvc(i) / zq(i)
1040        ELSE
1041          rvc_seri(i,k) = min_ratio
1042        ENDIF
1043        !--The MIN barrier is NEEDED because of:
1044        !-- 1) very rare pathological cases of the lsc scheme (rvc = 1. + 1e-16 sometimes)
1045        !-- 2) the thermal scheme does NOT guarantee that qvc <= qvap (or even qincld <= qtot)
1046        !--The MAX barrier is a safeguard that should not be activated
1047        rvc_seri(i,k) = MIN(MAX(rvc_seri(i,k), 0.), 1.)
1048
1049        !--Diagnostics
1050        gamma_cond(i,k) = gammasat(i)
1051        subfra(i,k) = 1. - cf_seri(i,k) - issrfra(i,k)
1052        qsub(i,k) = zq(i) - qvc(i) - qissr(i,k)
1053        qcld(i,k) = qvc(i) + zoliq(i)
1054      ENDDO
1055    ENDIF
1056
1057    ! Outputs:
1058    !-------------------------------
1059    ! Precipitation fluxes at layer interfaces
1060    ! + precipitation fractions +
1061    ! temperature and water species tendencies
1062    DO i = 1, klon
1063        psfl(i,k)=zifl(i)
1064        prfl(i,k)=zrfl(i)
1065        pfraclr(i,k)=znebprecipclr(i)
1066        pfracld(i,k)=znebprecipcld(i)
1067        d_q(i,k) = zq(i) - qt(i,k)
1068        d_t(i,k) = zt(i) - temp(i,k)
1069
1070        IF (ok_bug_phase_lscp) THEN
1071           d_ql(i,k) = (1-zfice(i))*zoliq(i)
1072           d_qi(i,k) = zfice(i)*zoliq(i)
1073        ELSE
1074           d_ql(i,k) = zoliql(i)
1075           d_qi(i,k) = zoliqi(i)   
1076        ENDIF
1077
1078    ENDDO
1079
1080
1081ENDDO ! loop on k from top to bottom
1082
1083
1084  ! Rain or snow at the surface (depending on the first layer temperature)
1085  DO i = 1, klon
1086      snow(i) = zifl(i)
1087      rain(i) = zrfl(i)
1088      ! c_iso final output
1089  ENDDO
1090
1091  IF (ncoreczq>0) THEN
1092      WRITE(lunout,*)'WARNING : ZQ in LSCP ',ncoreczq,' val < 1.e-15.'
1093  ENDIF
1094
1095
1096RETURN
1097
1098END SUBROUTINE lscp
1099!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1100
1101END MODULE lmdz_lscp
Note: See TracBrowser for help on using the repository browser.