source: LMDZ6/trunk/libf/phylmd/lmdz_lscp.f90 @ 5503

Last change on this file since 5503 was 5495, checked in by evignon, 3 days ago

correction du diagnostique de l'humidite relative en ciel clair.
n'affecte que le coating des aerosols.
Audran et Etienne
File size: 51.6 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, zqsl, zdqsl, zqsi, zdqsi
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, zqnl
281  REAL, DIMENSION(klon) :: zoliql, zoliqi
282  REAL, DIMENSION(klon) :: zt
283  REAL, DIMENSION(klon) :: zfice,zneb, znebl
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, qvcl, 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(:)=1.0 ! initialized at 1 as by default we assume mpc to be at ice saturation
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 ) 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                        cf_seri(:,k), rvc_seri(:,k), ql_seri(:,k), qi_seri(:,k), &
497                        zrfl, zrflclr, zrflcld, &
498                        zifl, ziflclr, ziflcld, &
499                        dqreva(:,k), dqssub(:,k) &
500                        )
501
502    !================================================================
503    ELSE
504
505      CALL histprecip_precld(klon, dtime, iftop, paprs(:,k), paprs(:,k+1), pplay(:,k), &
506                        zt, ztupnew, zq, zmqc, zneb, znebprecip, znebprecipclr, &
507                        zrfl, zrflclr, zrflcld, &
508                        zifl, ziflclr, ziflcld, &
509                        dqreva(:,k), dqssub(:,k) &
510                        )
511
512    ENDIF ! (ok_poprecip)
513   
514    ! Calculation of qsat, L/Cp*dqsat/dT and ncoreczq counter
515    !-------------------------------------------------------
516
517     qtot(:)=zq(:)+zmqc(:)
518     CALL calc_qsat_ecmwf(klon,zt,qtot,pplay(:,k),RTT,0,.false.,zqs,zdqs)
519     DO i = 1, klon
520            zdelta = MAX(0.,SIGN(1.,RTT-zt(i)))
521            zdqsdT_raw(i) = zdqs(i)*RCPD*(1.0+RVTMP2*zq(i)) / (RLVTT*(1.-zdelta) + RLSTT*zdelta)
522            IF (zq(i) .LT. 1.e-15) THEN
523                ncoreczq=ncoreczq+1
524                zq(i)=1.e-15
525            ENDIF
526            ! c_iso: do something similar for isotopes
527
528     ENDDO
529   
530    ! --------------------------------------------------------------------
531    ! P2> Cloud formation
532    !---------------------------------------------------------------------
533    !
534    ! Unlike fisrtilp, we always assume a 'fractional cloud' approach
535    ! i.e. clouds occupy only a fraction of the mesh (the subgrid distribution
536    ! is prescribed and depends on large scale variables and boundary layer
537    ! properties)
538    ! The decrease in condensed part due tu latent heating is taken into
539    ! account
540    ! -------------------------------------------------------------------
541
542        ! P2.1> With the PDFs (log-normal, bigaussian)
543        ! cloud properties calculation with the initial values of t and q
544        ! ----------------------------------------------------------------
545
546        ! initialise gammasat and qincloud_mpc
547        gammasat(:)=1.
548        qincloud_mpc(:)=0.
549
550        IF (iflag_cld_th.GE.5) THEN
551               ! Cloud cover and content in meshes affected by shallow convection,
552               ! are retrieved from a bi-gaussian distribution of the saturation deficit
553               ! following Jam et al. 2013
554
555               IF (iflag_cloudth_vert.LE.2) THEN
556                  ! Old version of Arnaud Jam
557
558                    CALL cloudth(klon,klev,k,tv,                  &
559                         zq,qta,fraca,                            &
560                         qcloud,ctot,pspsk,paprs,pplay,tla,thl, &
561                         ratqs,zqs,temp,                              &
562                         cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
563
564
565                ELSEIF (iflag_cloudth_vert.GE.3 .AND. iflag_cloudth_vert.LE.5) THEN
566                   ! Default version of Arnaud Jam
567                       
568                    CALL cloudth_v3(klon,klev,k,tv,                        &
569                         zq,qta,fraca,                                     &
570                         qcloud,ctot,ctot_vol,pspsk,paprs,pplay,tla,thl, &
571                         ratqs,sigma_qtherm,zqs,temp, &
572                         cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
573
574
575                ELSEIF (iflag_cloudth_vert.EQ.6) THEN
576                   ! Jean Jouhaud's version, with specific separation between surface and volume
577                   ! cloud fraction Decembre 2018
578
579                    CALL cloudth_v6(klon,klev,k,tv,                        &
580                         zq,qta,fraca,                                     &
581                         qcloud,ctot,ctot_vol,pspsk,paprs,pplay,tla,thl, &
582                         ratqs,zqs,temp, &
583                         cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
584
585                ELSEIF (iflag_cloudth_vert .EQ. 7) THEN
586                   ! Updated version of Arnaud Jam (correction by E. Vignon) + adapted treatment
587                   ! for boundary-layer mixed phase clouds
588                    CALL cloudth_mpc(klon,klev,k,mpc_bl_points,zt,zq,qta(:,k),fraca(:,k), &
589                                     pspsk(:,k),paprs(:,k+1),paprs(:,k),pplay(:,k), tla(:,k), &
590                                     ratqs(:,k),qcloud,qincloud_mpc,zfice_th,ctot(:,k),ctot_vol(:,k), &
591                                     cloudth_sth(:,k),cloudth_senv(:,k),cloudth_sigmath(:,k),cloudth_sigmaenv(:,k))
592
593                ENDIF
594
595
596                DO i=1,klon
597                    rneb(i,k)=ctot(i,k)
598                    rneblsvol(i,k)=ctot_vol(i,k)
599                    zqn(i)=qcloud(i)
600                    !--AB grid-mean vapor in the cloud - we assume saturation adjustment
601                    qvc(i) = rneb(i,k) * zqs(i)
602                ENDDO
603
604        ENDIF
605
606        IF (iflag_cld_th .LE. 4) THEN
607           
608                ! lognormal
609            lognormale(:) = .TRUE.
610
611        ELSEIF (iflag_cld_th .GE. 6) THEN
612           
613                ! lognormal distribution when no thermals
614            lognormale(:) = fraca(:,k) < min_frac_th_cld
615
616        ELSE
617                ! When iflag_cld_th=5, we always assume
618                ! bi-gaussian distribution
619            lognormale(:) = .FALSE.
620       
621        ENDIF
622
623        DT(:) = 0.
624        n_i(:)=0
625        Tbef(:)=zt(:)
626        qlbef(:)=0.
627
628        ! Treatment of non-boundary layer clouds (lognormale)
629        ! We iterate here to take into account the change in qsat(T) and ice fraction
630        ! during the condensation process
631        ! the increment in temperature is calculated using the first principle of
632        ! thermodynamics (enthalpy conservation equation in a mesh composed of a cloud fraction
633        ! and a clear sky fraction)
634        ! note that we assume that the vapor in the cloud is at saturation for this calculation     
635
636        DO iter=1,iflag_fisrtilp_qsat+1
637
638                keepgoing(:) = .FALSE.
639
640                DO i=1,klon
641
642                    ! keepgoing = .true. while convergence is not satisfied
643
644                    IF (((ABS(DT(i)).GT.DDT0) .OR. (n_i(i) .EQ. 0)) .AND. lognormale(i)) THEN
645
646                        ! if not convergence:
647                        ! we calculate a new iteration
648                        keepgoing(i) = .TRUE.
649
650                        ! P2.2.1> cloud fraction and condensed water mass calculation
651                        ! Calculated variables:
652                        ! rneb : cloud fraction
653                        ! zqn : total water within the cloud
654                        ! zcond: mean condensed water within the mesh
655                        ! rhcl: clear-sky relative humidity
656                        !---------------------------------------------------------------
657
658                        ! new temperature that only serves in the iteration process:
659                        Tbef(i)=Tbef(i)+DT(i)
660
661                        ! Rneb, qzn and zcond for lognormal PDFs
662                        qtot(i)=zq(i)+zmqc(i)
663
664                      ENDIF
665
666                  ENDDO
667
668                  ! Calculation of saturation specific humidity and ice fraction
669                  CALL calc_qsat_ecmwf(klon,Tbef,qtot,pplay(:,k),RTT,0,.false.,zqs,zdqs)
670                 
671                  IF (iflag_icefrac .GE. 3) THEN
672                  ! consider a ice weighted qs to ensure that liquid clouds at T<0 have a consistent cloud fraction
673                  ! and cloud condensed water content. idea following Dietlitcher et al. 2018, GMD
674                  ! we make this option works only for the physically-based and tke-depdenent param (iflag_icefrac>=1)
675                      DO i=1,klon
676                           CALL calc_qsat_ecmwf(klon,Tbef,qtot,pplay(:,k),RTT,1,.false.,zqsl,zdqsl)
677                           CALL calc_qsat_ecmwf(klon,Tbef,qtot,pplay(:,k),RTT,2,.false.,zqsi,zdqsi)
678                           zqs(i)=zfice(i)*zqsi(i)+(1.-zfice(i))*zqsl(i)
679                           zdqs(i)=zfice(i)*zdqsi(i)+zqsi(i)*dzfice(i)+(1.-zfice(i))*zdqsl(i)-zqsl(i)*dzfice(i)
680                      ENDDO
681                  ENDIF
682
683                  CALL calc_gammasat(klon,Tbef,qtot,pplay(:,k),gammasat,dgammasatdt)
684                  ! saturation may occur at a humidity different from qsat (gamma qsat), so gamma correction for dqs
685                  zdqs(:) = gammasat(:)*zdqs(:)+zqs(:)*dgammasatdt(:)
686
687                  ! Cloud condensation based on subgrid pdf
688                  !--AB Activates a condensation scheme that allows for
689                  !--ice supersaturation and contrails evolution from aviation
690                  IF (ok_ice_supersat) THEN
691
692                    !--Calculate the shear value (input for condensation and ice supersat)
693                    DO i = 1, klon
694                      !--Cell thickness [m]
695                      delta_z = ( paprs(i,k) - paprs(i,k+1) ) / RG / pplay(i,k) * Tbef(i) * RD
696                      IF ( iftop ) THEN
697                        ! top
698                        shear(i) = SQRT( ( (u_seri(i,k) - u_seri(i,k-1)) / delta_z )**2. &
699                                       + ( (v_seri(i,k) - v_seri(i,k-1)) / delta_z )**2. )
700                      ELSEIF ( k .EQ. 1 ) THEN
701                        ! surface
702                        shear(i) = SQRT( ( (u_seri(i,k+1) - u_seri(i,k)) / delta_z )**2. &
703                                       + ( (v_seri(i,k+1) - v_seri(i,k)) / delta_z )**2. )
704                      ELSE
705                        ! other layers
706                        shear(i) = SQRT( ( ( (u_seri(i,k+1) + u_seri(i,k)) / 2. &
707                                           - (u_seri(i,k) + u_seri(i,k-1)) / 2. ) / delta_z )**2. &
708                                       + ( ( (v_seri(i,k+1) + v_seri(i,k)) / 2. &
709                                           - (v_seri(i,k) + v_seri(i,k-1)) / 2. ) / delta_z )**2. )
710                      ENDIF
711                    ENDDO
712
713                    !---------------------------------------------
714                    !--   CONDENSATION AND ICE SUPERSATURATION  --
715                    !---------------------------------------------
716
717                    CALL condensation_ice_supersat( &
718                        klon, dtime, missing_val, &
719                        pplay(:,k), paprs(:,k), paprs(:,k+1), &
720                        cf_seri(:,k), rvc_seri(:,k), ql_seri(:,k), qi_seri(:,k), &
721                        shear, tke_dissip(:,k), cell_area, &
722                        Tbef, zq, zqs, gammasat, ratqs(:,k), keepgoing, &
723                        rneb(:,k), zqn, qvc, issrfra(:,k), qissr(:,k), &
724                        dcf_sub(:,k), dcf_con(:,k), dcf_mix(:,k), &
725                        dqi_adj(:,k), dqi_sub(:,k), dqi_con(:,k), dqi_mix(:,k), &
726                        dqvc_adj(:,k), dqvc_sub(:,k), dqvc_con(:,k), dqvc_mix(:,k), &
727                        Tcontr(:,k), qcontr(:,k), qcontr2(:,k), fcontrN(:,k), fcontrP(:,k), &
728                        flight_dist(:,k), flight_h2o(:,k), &
729                        dcf_avi(:,k), dqi_avi(:,k), dqvc_avi(:,k))
730
731
732                  ELSE
733                  !--generalised lognormal condensation scheme (Bony and Emanuel 2001)
734
735                   CALL condensation_lognormal( &
736                       klon, Tbef, zq, zqs, gammasat, ratqs(:,k), &
737                       keepgoing, rneb(:,k), zqn, qvc)
738
739
740                  ENDIF ! .NOT. ok_ice_supersat
741
742                  ! cloud phase determination
743                  IF (iflag_icefrac .GE. 2) THEN
744                     ! phase partitioning depending on temperature. activates here in the iteration process if iflag_icefrac > 2
745                     CALL icefrac_lscp_turb(klon, dtime, Tbef, pplay(:,k), paprs(:,k), paprs(:,k+1), omega(:,k), qice_ini, ziflcld, zqn, &
746                     rneb(:,k), tke(:,k), tke_dissip(:,k), zqliq, zqvapcl, zqice, zfice, dzfice, cldfraliq(:,k),sigma2_icefracturb(:,k), mean_icefracturb(:,k))
747                  ELSE
748                     ! phase partitioning depending on temperature and eventually distance to cloud top
749                     IF (iflag_t_glace.GE.4) THEN
750                       ! For iflag_t_glace GE 4 the phase partition function dependends on temperature AND distance to cloud top
751                       CALL distance_to_cloud_top(klon,klev,k,temp,pplay,paprs,rneb,zdistcltop,ztemp_cltop)
752                     ENDIF
753                     CALL icefrac_lscp(klon, zt, iflag_ice_thermo, zdistcltop,ztemp_cltop,zfice,dzfice)
754                  ENDIF
755
756
757                  DO i=1,klon
758                      IF (keepgoing(i)) THEN
759
760                        ! If vertical heterogeneity, change fraction by volume as well
761                        IF (iflag_cloudth_vert.GE.3) THEN
762                            ctot_vol(i,k)=rneb(i,k)
763                            rneblsvol(i,k)=ctot_vol(i,k)
764                        ENDIF
765
766
767                       ! P2.2.2> Approximative calculation of temperature variation DT
768                       !           due to condensation.
769                       ! Calculated variables:
770                       ! dT : temperature change due to condensation
771                       !---------------------------------------------------------------
772
773               
774                        IF (zfice(i).LT.1) THEN
775                            cste=RLVTT
776                        ELSE
777                            cste=RLSTT
778                        ENDIF
779                       
780                        IF ( ok_unadjusted_clouds ) THEN
781                          !--AB We relax the saturation adjustment assumption
782                          !-- qvc (grid-mean vapor in cloud) is calculated by the condensation scheme
783                          IF ( rneb(i,k) .GT. eps ) THEN
784                            qlbef(i) = MAX(0., zqn(i) - qvc(i) / rneb(i,k))
785                          ELSE
786                            qlbef(i) = 0.
787                          ENDIF
788                        ELSE
789                          qlbef(i)=max(0.,zqn(i)-zqs(i))
790                        ENDIF
791
792                        num = -Tbef(i)+zt(i)+rneb(i,k)*((1-zfice(i))*RLVTT &
793                              +zfice(i)*RLSTT)/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))*qlbef(i)
794                        denom = 1.+rneb(i,k)*((1-zfice(i))*RLVTT+zfice(i)*RLSTT)/cste*zdqs(i) &
795                              -(RLSTT-RLVTT)/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))*rneb(i,k)      &
796                              *qlbef(i)*dzfice(i)
797                        ! here we update a provisory temperature variable that only serves in the iteration
798                        ! process
799                        DT(i)=num/denom
800                        n_i(i)=n_i(i)+1
801
802                    ENDIF ! end keepgoing
803
804                ENDDO     ! end loop on i
805
806        ENDDO         ! iter=1,iflag_fisrtilp_qsat+1
807
808        ! P2.2> Final quantities calculation
809        ! Calculated variables:
810        !   rneb : cloud fraction
811        !   zcond: mean condensed water in the mesh
812        !   zqn  : mean water vapor in the mesh
813        !   zfice: ice fraction in clouds
814        !   zt   : temperature
815        !   rhcl : clear-sky relative humidity
816        ! ----------------------------------------------------------------
817
818
819        ! Cloud phase final determination
820        !--------------------------------
821        ! For iflag_t_glace GE 4 the phase partition function dependends on temperature AND distance to cloud top
822        IF (iflag_t_glace.GE.4) THEN
823           CALL distance_to_cloud_top(klon,klev,k,temp,pplay,paprs,rneb,zdistcltop,ztemp_cltop)
824           distcltop(:,k)=zdistcltop(:)
825           temp_cltop(:,k)=ztemp_cltop(:)
826        ENDIF
827        ! Partition function depending on temperature for all clouds (shallow convective and stratiform)
828        CALL icefrac_lscp(klon, zt, iflag_ice_thermo, zdistcltop, ztemp_cltop, zfice, dzfice)
829
830        ! Partition function depending on tke for non shallow-convective clouds, erase previous estimation
831        IF (iflag_icefrac .GE. 1) THEN
832           CALL icefrac_lscp_turb(klon, dtime, Tbef, pplay(:,k), paprs(:,k), paprs(:,k+1), omega(:,k), qice_ini, ziflcld, zqn, &
833           rneb(:,k), tke(:,k), tke_dissip(:,k), zqliq, zqvapcl, zqice, zfice_turb, dzfice_turb, cldfraliq(:,k),sigma2_icefracturb(:,k), mean_icefracturb(:,k))
834        ENDIF
835
836        ! Water vapor update, subsequent latent heat exchange for each cloud type
837        !------------------------------------------------------------------------
838        DO i=1, klon
839            ! Overwrite phase partitioning in boundary layer mixed phase clouds when the
840            ! iflag_cloudth_vert=7 and specific param is activated
841            IF (mpc_bl_points(i,k) .GT. 0) THEN
842                zcond(i) = MAX(0.0,qincloud_mpc(i))*rneb(i,k)
843                !--Relative humidity (no unit) in clear sky, calculated as rh = q / qsat
844                !--This is slighly approximated, the actual formula is
845                !-- rh = q / qsat * (eps + (1-eps)*qsat) / (eps + (1-eps)*q)
846                !--Here, rh = (qtot - qincld * cldfra) / clrfra / qsat
847                !--where (qtot - qincld * cldfra) is the grid-mean clear sky water content
848                rhcl(i,k) = ( zq(i) - qincloud_mpc(i) * rneb(i,k) ) / ( 1. - rneb(i,k) ) / zqs(i)
849                ! water vapor update and partition function if thermals
850                zq(i) = zq(i) - zcond(i)       
851                zfice(i)=zfice_th(i)
852            ELSE
853                ! Checks on rneb, rhcl and zqn
854                IF (rneb(i,k) .LE. 0.0) THEN
855                    zqn(i) = 0.0
856                    rneb(i,k) = 0.0
857                    zcond(i) = 0.0
858                    rhcl(i,k)=zq(i)/zqs(i)
859                ELSE IF (rneb(i,k) .GE. 1.0) THEN
860                    zqn(i) = zq(i)
861                    rneb(i,k) = 1.0
862                    IF ( ok_unadjusted_clouds ) THEN
863                      !--AB We relax the saturation adjustment assumption
864                      !-- qvc (grid-mean vapor in cloud) is calculated by the condensation scheme
865                      zcond(i) = MAX(0., zqn(i) - qvc(i))
866                    ELSE
867                      zcond(i) = MAX(0.0,zqn(i)-zqs(i))
868                    ENDIF
869                    rhcl(i,k)=1.0
870                ELSE
871                    IF ( ok_unadjusted_clouds ) THEN
872                      !--AB We relax the saturation adjustment assumption
873                      !-- qvc (grid-mean vapor in cloud) is calculated by the condensation scheme
874                      zcond(i) = MAX(0., zqn(i) * rneb(i,k) - qvc(i))
875                    ELSE
876                      zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)
877                    ENDIF
878                    !--Relative humidity (no unit) in clear sky, calculated as rh = q / qsat
879                    !--This is slighly approximated, the actual formula is
880                    !-- rh = q / qsat * (eps + (1-eps)*qsat) / (eps + (1-eps)*q)
881                    !--Here, rh = (qtot - qincld * cldfra) / clrfra / qsat
882                    !--where (qtot - qincld * cldfra) is the grid-mean clear sky water content
883                    rhcl(i,k) = ( zq(i) - zqn(i) * rneb(i,k) ) / ( 1. - rneb(i,k) ) / zqs(i)
884                    ! Overwrite partitioning for non shallow-convective clouds if iflag_icefrac>1 (icefrac turb param)
885                    IF (iflag_icefrac .GE. 1) THEN
886                        IF (lognormale(i)) THEN 
887                           zcond(i)  = zqliq(i) + zqice(i)
888                           zfice(i)  = zfice_turb(i)
889                           rhcl(i,k) = zqvapcl(i) * rneb(i,k) + (zq(i) - zqn(i)) * (1.-rneb(i,k))
890                        ENDIF
891                    ENDIF
892                ENDIF
893
894                ! water vapor update
895                zq(i) = zq(i) - zcond(i)
896
897            ENDIF
898           
899           
900            ! temperature update due to phase change
901            zt(i) = zt(i) + (1.-zfice(i))*zcond(i) &
902            &     * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i))) &
903                  +zfice(i)*zcond(i) * RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i)))
904        ENDDO
905
906        ! If vertical heterogeneity, change volume fraction
907        IF (iflag_cloudth_vert .GE. 3) THEN
908          ctot_vol(1:klon,k)=min(max(ctot_vol(1:klon,k),0.),1.)
909          rneblsvol(1:klon,k)=ctot_vol(1:klon,k)
910        ENDIF
911
912
913    ! ----------------------------------------------------------------
914    ! P3> Precipitation processes, after cloud formation
915    !     - precipitation formation, melting/freezing
916    ! ----------------------------------------------------------------
917
918    ! Initiate the quantity of liquid and solid condensates
919    ! Note that in the following, zcond is the total amount of condensates
920    ! including newly formed precipitations (i.e., condensates formed by the
921    ! condensation process), while zoliq is the total amount of condensates in
922    ! the cloud (i.e., on which precipitation processes have applied)
923    DO i = 1, klon
924      zoliq(i)  = zcond(i)
925      zoliql(i) = zcond(i) * ( 1. - zfice(i) )
926      zoliqi(i) = zcond(i) * zfice(i)
927    ENDDO
928
929    !================================================================
930    ! Flag for the new and more microphysical treatment of precipitation from Atelier Nuage (R)
931    IF (ok_poprecip) THEN
932
933      CALL poprecip_postcld(klon, dtime, paprs(:,k), paprs(:,k+1), pplay(:,k), &
934                            ctot_vol(:,k), ptconv(:,k), &
935                            zt, zq, zoliql, zoliqi, zfice, &
936                            rneb(:,k), znebprecipclr, znebprecipcld, &
937                            zrfl, zrflclr, zrflcld, &
938                            zifl, ziflclr, ziflcld, &
939                            qraindiag(:,k), qsnowdiag(:,k), dqrauto(:,k), &
940                            dqrcol(:,k), dqrmelt(:,k), dqrfreez(:,k), &
941                            dqsauto(:,k), dqsagg(:,k), dqsrim(:,k), &
942                            dqsmelt(:,k), dqsfreez(:,k) &
943                            )
944      DO i = 1, klon
945          zoliq(i) = zoliql(i) + zoliqi(i)
946      ENDDO
947
948    !================================================================
949    ELSE
950
951      CALL histprecip_postcld(klon, dtime, iftop, paprs(:,k), paprs(:,k+1), pplay(:,k), &
952                            ctot_vol(:,k), ptconv(:,k), zdqsdT_raw, &
953                            zt, zq, zoliq, zoliql, zoliqi, zcond, zfice, zmqc, &
954                            rneb(:,k), znebprecipclr, znebprecipcld, &
955                            zneb, tot_zneb, zrho_up, zvelo_up, &
956                            zrfl, zrflclr, zrflcld, zifl, ziflclr, ziflcld, &
957                            zradocond, zradoice, dqrauto(:,k), dqsauto(:,k) &
958                            )
959
960    ENDIF ! ok_poprecip
961
962    ! End of precipitation processes after cloud formation
963    ! ----------------------------------------------------
964
965    !----------------------------------------------------------------------
966    ! P4> Calculation of cloud condensates amount seen by radiative scheme
967    !----------------------------------------------------------------------
968
969    DO i=1,klon
970
971      IF (ok_poprecip) THEN
972        IF (ok_radocond_snow) THEN
973          radocond(i,k) = zoliq(i)
974          zradoice(i) = zoliqi(i) + qsnowdiag(i,k)
975        ELSE
976          radocond(i,k) = zoliq(i)
977          zradoice(i) = zoliqi(i)
978        ENDIF
979      ELSE
980        radocond(i,k) = zradocond(i)
981      ENDIF
982
983      ! calculate the percentage of ice in "radocond" so cloud+precip seen
984      ! by the radiation scheme
985      IF (radocond(i,k) .GT. 0.) THEN
986        radicefrac(i,k)=MIN(MAX(zradoice(i)/radocond(i,k),0.),1.)
987      ENDIF
988    ENDDO
989
990    ! ----------------------------------------------------------------
991    ! P5> Wet scavenging
992    ! ----------------------------------------------------------------
993
994    !Scavenging through nucleation in the layer
995
996    DO i = 1,klon
997       
998        IF(zcond(i).GT.zoliq(i)+1.e-10) THEN
999            beta(i,k) = (zcond(i)-zoliq(i))/zcond(i)/dtime
1000        ELSE
1001            beta(i,k) = 0.
1002        ENDIF
1003
1004        zprec_cond(i) = MAX(zcond(i)-zoliq(i),0.0)*(paprs(i,k)-paprs(i,k+1))/RG
1005
1006        IF (rneb(i,k).GT.0.0.AND.zprec_cond(i).GT.0.) THEN
1007
1008            IF (temp(i,k) .GE. t_glace_min) THEN
1009                zalpha_tr = a_tr_sca(3)
1010            ELSE
1011                zalpha_tr = a_tr_sca(4)
1012            ENDIF
1013
1014            zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/MAX(rneb(i,k),seuil_neb))
1015            frac_nucl(i,k)= 1.-MAX(rneb(i,k),seuil_neb)*zfrac_lessi
1016
1017            ! Nucleation with a factor of -1 instead of -0.5
1018            zfrac_lessi = 1. - EXP(-zprec_cond(i)/MAX(rneb(i,k),seuil_neb))
1019
1020        ENDIF
1021
1022    ENDDO
1023
1024    ! Scavenging through impaction in the underlying layer
1025
1026    DO kk = k-1, 1, -1
1027
1028        DO i = 1, klon
1029
1030            IF (rneb(i,k).GT.0.0.AND.zprec_cond(i).GT.0.) THEN
1031
1032                IF (temp(i,kk) .GE. t_glace_min) THEN
1033                    zalpha_tr = a_tr_sca(1)
1034                ELSE
1035                    zalpha_tr = a_tr_sca(2)
1036                ENDIF
1037
1038              zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/MAX(rneb(i,k),seuil_neb))
1039              frac_impa(i,kk)= 1.-MAX(rneb(i,k),seuil_neb)*zfrac_lessi
1040
1041            ENDIF
1042
1043        ENDDO
1044
1045    ENDDO
1046   
1047    !------------------------------------------------------------
1048    ! P6 > write diagnostics and outputs
1049    !------------------------------------------------------------
1050   
1051    !--AB Write diagnostics and tracers for ice supersaturation
1052    IF ( ok_ice_supersat ) THEN
1053      CALL calc_qsat_ecmwf(klon,zt,qzero,pplay(:,k),RTT,1,.false.,qsatl(:,k),zdqs)
1054      CALL calc_qsat_ecmwf(klon,zt,qzero,pplay(:,k),RTT,2,.false.,qsati(:,k),zdqs)
1055
1056      DO i = 1, klon
1057
1058        IF ( zoliq(i) .LE. 0. ) THEN
1059          !--If everything was precipitated, the remaining empty cloud is dissipated
1060          !--and everything is transfered to the subsaturated clear sky region
1061          rneb(i,k) = 0.
1062          qvc(i) = 0.
1063        ENDIF
1064
1065        cf_seri(i,k) = rneb(i,k)
1066
1067        IF ( .NOT. ok_unadjusted_clouds ) THEN
1068          qvc(i) = zqs(i) * rneb(i,k)
1069        ENDIF
1070        IF ( zq(i) .GT. min_qParent ) THEN
1071          rvc_seri(i,k) = qvc(i) / zq(i)
1072        ELSE
1073          rvc_seri(i,k) = min_ratio
1074        ENDIF
1075        !--The MIN barrier is NEEDED because of:
1076        !-- 1) very rare pathological cases of the lsc scheme (rvc = 1. + 1e-16 sometimes)
1077        !-- 2) the thermal scheme does NOT guarantee that qvc <= qvap (or even qincld <= qtot)
1078        !--The MAX barrier is a safeguard that should not be activated
1079        rvc_seri(i,k) = MIN(MAX(rvc_seri(i,k), 0.), 1.)
1080
1081        !--Diagnostics
1082        gamma_cond(i,k) = gammasat(i)
1083        subfra(i,k) = 1. - cf_seri(i,k) - issrfra(i,k)
1084        qsub(i,k) = zq(i) - qvc(i) - qissr(i,k)
1085        qcld(i,k) = qvc(i) + zoliq(i)
1086      ENDDO
1087    ENDIF
1088
1089    ! Outputs:
1090    !-------------------------------
1091    ! Precipitation fluxes at layer interfaces
1092    ! + precipitation fractions +
1093    ! temperature and water species tendencies
1094    DO i = 1, klon
1095        psfl(i,k)=zifl(i)
1096        prfl(i,k)=zrfl(i)
1097        pfraclr(i,k)=znebprecipclr(i)
1098        pfracld(i,k)=znebprecipcld(i)
1099        d_q(i,k) = zq(i) - qt(i,k)
1100        d_t(i,k) = zt(i) - temp(i,k)
1101
1102        IF (ok_bug_phase_lscp) THEN
1103           d_ql(i,k) = (1-zfice(i))*zoliq(i)
1104           d_qi(i,k) = zfice(i)*zoliq(i)
1105        ELSE
1106           d_ql(i,k) = zoliql(i)
1107           d_qi(i,k) = zoliqi(i)   
1108        ENDIF
1109
1110    ENDDO
1111
1112
1113ENDDO ! loop on k from top to bottom
1114
1115
1116  ! Rain or snow at the surface (depending on the first layer temperature)
1117  DO i = 1, klon
1118      snow(i) = zifl(i)
1119      rain(i) = zrfl(i)
1120      ! c_iso final output
1121  ENDDO
1122
1123  IF (ncoreczq>0) THEN
1124      WRITE(lunout,*)'WARNING : ZQ in LSCP ',ncoreczq,' val < 1.e-15.'
1125  ENDIF
1126
1127
1128RETURN
1129
1130END SUBROUTINE lscp
1131!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1132
1133END MODULE lmdz_lscp
Note: See TracBrowser for help on using the repository browser.