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

Last change on this file since 5523 was 5520, checked in by evignon, 2 days ago

revert du commit sur rhcl dans lscp

File size: 52.2 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                ! following line is very strange and probably wrong
844                rhcl(i,k)= (zqs(i)+zq(i))/2./zqs(i)
845                ! Correct calculation of clear-sky relative humidity. To activate once
846                ! issues related to zqn>zq in bi-gaussian clouds are corrected
847                !--Relative humidity (no unit) in clear sky, calculated as rh = q / qsat
848                !--This is slighly approximated, the actual formula is
849                !-- rh = q / qsat * (eps + (1-eps)*qsat) / (eps + (1-eps)*q)
850                !--Here, rh = (qtot - qincld * cldfra) / clrfra / qsat
851                !--where (qtot - qincld * cldfra) is the grid-mean clear sky water content
852                ! rhcl(i,k) = ( zq(i) - qincloud_mpc(i) * rneb(i,k) ) / ( 1. - rneb(i,k) ) / zqs(i)
853                ! water vapor update and partition function if thermals
854                zq(i) = zq(i) - zcond(i)       
855                zfice(i)=zfice_th(i)
856            ELSE
857                ! Checks on rneb, rhcl and zqn
858                IF (rneb(i,k) .LE. 0.0) THEN
859                    zqn(i) = 0.0
860                    rneb(i,k) = 0.0
861                    zcond(i) = 0.0
862                    rhcl(i,k)=zq(i)/zqs(i)
863                ELSE IF (rneb(i,k) .GE. 1.0) THEN
864                    zqn(i) = zq(i)
865                    rneb(i,k) = 1.0
866                    IF ( ok_unadjusted_clouds ) THEN
867                      !--AB We relax the saturation adjustment assumption
868                      !-- qvc (grid-mean vapor in cloud) is calculated by the condensation scheme
869                      zcond(i) = MAX(0., zqn(i) - qvc(i))
870                    ELSE
871                      zcond(i) = MAX(0.0,zqn(i)-zqs(i))
872                    ENDIF
873                    rhcl(i,k)=1.0
874                ELSE
875                    IF ( ok_unadjusted_clouds ) THEN
876                      !--AB We relax the saturation adjustment assumption
877                      !-- qvc (grid-mean vapor in cloud) is calculated by the condensation scheme
878                      zcond(i) = MAX(0., zqn(i) * rneb(i,k) - qvc(i))
879                    ELSE
880                      zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)
881                    ENDIF
882                    ! Overwrite partitioning for non shallow-convective clouds if iflag_icefrac>1 (icefrac turb param)
883                    IF (iflag_icefrac .GE. 1) THEN
884                        IF (lognormale(i)) THEN 
885                           zcond(i)  = zqliq(i) + zqice(i)
886                           zfice(i)  = zfice_turb(i)
887                        ENDIF
888                    ENDIF
889
890                    ! following line is very strange and probably wrong
891                    rhcl(i,k)= (zqs(i)+zq(i))/2./zqs(i)
892                    ! Correct calculation of clear-sky relative humidity. To activate once
893                    ! issues related to zqn>zq in bi-gaussian clouds are corrected
894                    !--Relative humidity (no unit) in clear sky, calculated as rh = q / qsat
895                    !--This is slighly approximated, the actual formula is
896                    !-- rh = q / qsat * (eps + (1-eps)*qsat) / (eps + (1-eps)*q)
897                    !--Here, rh = (qtot - qincld * cldfra) / clrfra / qsat
898                    !--where (qtot - qincld * cldfra) is the grid-mean clear sky water content
899                    ! rhcl(i,k) = ( zq(i) - zqn(i) * rneb(i,k) ) / ( 1. - rneb(i,k) ) / zqs(i)
900                    ! Overwrite partitioning for non shallow-convective clouds if iflag_icefrac>1 (icefrac turb param)
901
902                ENDIF
903
904                ! water vapor update
905                zq(i) = zq(i) - zcond(i)
906
907            ENDIF
908           
909           
910            ! temperature update due to phase change
911            zt(i) = zt(i) + (1.-zfice(i))*zcond(i) &
912            &     * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i))) &
913                  +zfice(i)*zcond(i) * RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i)))
914        ENDDO
915
916        ! If vertical heterogeneity, change volume fraction
917        IF (iflag_cloudth_vert .GE. 3) THEN
918          ctot_vol(1:klon,k)=min(max(ctot_vol(1:klon,k),0.),1.)
919          rneblsvol(1:klon,k)=ctot_vol(1:klon,k)
920        ENDIF
921
922
923    ! ----------------------------------------------------------------
924    ! P3> Precipitation processes, after cloud formation
925    !     - precipitation formation, melting/freezing
926    ! ----------------------------------------------------------------
927
928    ! Initiate the quantity of liquid and solid condensates
929    ! Note that in the following, zcond is the total amount of condensates
930    ! including newly formed precipitations (i.e., condensates formed by the
931    ! condensation process), while zoliq is the total amount of condensates in
932    ! the cloud (i.e., on which precipitation processes have applied)
933    DO i = 1, klon
934      zoliq(i)  = zcond(i)
935      zoliql(i) = zcond(i) * ( 1. - zfice(i) )
936      zoliqi(i) = zcond(i) * zfice(i)
937    ENDDO
938
939    !================================================================
940    ! Flag for the new and more microphysical treatment of precipitation from Atelier Nuage (R)
941    IF (ok_poprecip) THEN
942
943      CALL poprecip_postcld(klon, dtime, paprs(:,k), paprs(:,k+1), pplay(:,k), &
944                            ctot_vol(:,k), ptconv(:,k), &
945                            zt, zq, zoliql, zoliqi, zfice, &
946                            rneb(:,k), znebprecipclr, znebprecipcld, &
947                            zrfl, zrflclr, zrflcld, &
948                            zifl, ziflclr, ziflcld, &
949                            qraindiag(:,k), qsnowdiag(:,k), dqrauto(:,k), &
950                            dqrcol(:,k), dqrmelt(:,k), dqrfreez(:,k), &
951                            dqsauto(:,k), dqsagg(:,k), dqsrim(:,k), &
952                            dqsmelt(:,k), dqsfreez(:,k) &
953                            )
954      DO i = 1, klon
955          zoliq(i) = zoliql(i) + zoliqi(i)
956      ENDDO
957
958    !================================================================
959    ELSE
960
961      CALL histprecip_postcld(klon, dtime, iftop, paprs(:,k), paprs(:,k+1), pplay(:,k), &
962                            ctot_vol(:,k), ptconv(:,k), zdqsdT_raw, &
963                            zt, zq, zoliq, zoliql, zoliqi, zcond, zfice, zmqc, &
964                            rneb(:,k), znebprecipclr, znebprecipcld, &
965                            zneb, tot_zneb, zrho_up, zvelo_up, &
966                            zrfl, zrflclr, zrflcld, zifl, ziflclr, ziflcld, &
967                            zradocond, zradoice, dqrauto(:,k), dqsauto(:,k) &
968                            )
969
970    ENDIF ! ok_poprecip
971
972    ! End of precipitation processes after cloud formation
973    ! ----------------------------------------------------
974
975    !----------------------------------------------------------------------
976    ! P4> Calculation of cloud condensates amount seen by radiative scheme
977    !----------------------------------------------------------------------
978
979    DO i=1,klon
980
981      IF (ok_poprecip) THEN
982        IF (ok_radocond_snow) THEN
983          radocond(i,k) = zoliq(i)
984          zradoice(i) = zoliqi(i) + qsnowdiag(i,k)
985        ELSE
986          radocond(i,k) = zoliq(i)
987          zradoice(i) = zoliqi(i)
988        ENDIF
989      ELSE
990        radocond(i,k) = zradocond(i)
991      ENDIF
992
993      ! calculate the percentage of ice in "radocond" so cloud+precip seen
994      ! by the radiation scheme
995      IF (radocond(i,k) .GT. 0.) THEN
996        radicefrac(i,k)=MIN(MAX(zradoice(i)/radocond(i,k),0.),1.)
997      ENDIF
998    ENDDO
999
1000    ! ----------------------------------------------------------------
1001    ! P5> Wet scavenging
1002    ! ----------------------------------------------------------------
1003
1004    !Scavenging through nucleation in the layer
1005
1006    DO i = 1,klon
1007       
1008        IF(zcond(i).GT.zoliq(i)+1.e-10) THEN
1009            beta(i,k) = (zcond(i)-zoliq(i))/zcond(i)/dtime
1010        ELSE
1011            beta(i,k) = 0.
1012        ENDIF
1013
1014        zprec_cond(i) = MAX(zcond(i)-zoliq(i),0.0)*(paprs(i,k)-paprs(i,k+1))/RG
1015
1016        IF (rneb(i,k).GT.0.0.AND.zprec_cond(i).GT.0.) THEN
1017
1018            IF (temp(i,k) .GE. t_glace_min) THEN
1019                zalpha_tr = a_tr_sca(3)
1020            ELSE
1021                zalpha_tr = a_tr_sca(4)
1022            ENDIF
1023
1024            zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/MAX(rneb(i,k),seuil_neb))
1025            frac_nucl(i,k)= 1.-MAX(rneb(i,k),seuil_neb)*zfrac_lessi
1026
1027            ! Nucleation with a factor of -1 instead of -0.5
1028            zfrac_lessi = 1. - EXP(-zprec_cond(i)/MAX(rneb(i,k),seuil_neb))
1029
1030        ENDIF
1031
1032    ENDDO
1033
1034    ! Scavenging through impaction in the underlying layer
1035
1036    DO kk = k-1, 1, -1
1037
1038        DO i = 1, klon
1039
1040            IF (rneb(i,k).GT.0.0.AND.zprec_cond(i).GT.0.) THEN
1041
1042                IF (temp(i,kk) .GE. t_glace_min) THEN
1043                    zalpha_tr = a_tr_sca(1)
1044                ELSE
1045                    zalpha_tr = a_tr_sca(2)
1046                ENDIF
1047
1048              zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/MAX(rneb(i,k),seuil_neb))
1049              frac_impa(i,kk)= 1.-MAX(rneb(i,k),seuil_neb)*zfrac_lessi
1050
1051            ENDIF
1052
1053        ENDDO
1054
1055    ENDDO
1056   
1057    !------------------------------------------------------------
1058    ! P6 > write diagnostics and outputs
1059    !------------------------------------------------------------
1060   
1061    !--AB Write diagnostics and tracers for ice supersaturation
1062    IF ( ok_ice_supersat ) THEN
1063      CALL calc_qsat_ecmwf(klon,zt,qzero,pplay(:,k),RTT,1,.false.,qsatl(:,k),zdqs)
1064      CALL calc_qsat_ecmwf(klon,zt,qzero,pplay(:,k),RTT,2,.false.,qsati(:,k),zdqs)
1065
1066      DO i = 1, klon
1067
1068        IF ( zoliq(i) .LE. 0. ) THEN
1069          !--If everything was precipitated, the remaining empty cloud is dissipated
1070          !--and everything is transfered to the subsaturated clear sky region
1071          rneb(i,k) = 0.
1072          qvc(i) = 0.
1073        ENDIF
1074
1075        cf_seri(i,k) = rneb(i,k)
1076
1077        IF ( .NOT. ok_unadjusted_clouds ) THEN
1078          qvc(i) = zqs(i) * rneb(i,k)
1079        ENDIF
1080        IF ( zq(i) .GT. min_qParent ) THEN
1081          rvc_seri(i,k) = qvc(i) / zq(i)
1082        ELSE
1083          rvc_seri(i,k) = min_ratio
1084        ENDIF
1085        !--The MIN barrier is NEEDED because of:
1086        !-- 1) very rare pathological cases of the lsc scheme (rvc = 1. + 1e-16 sometimes)
1087        !-- 2) the thermal scheme does NOT guarantee that qvc <= qvap (or even qincld <= qtot)
1088        !--The MAX barrier is a safeguard that should not be activated
1089        rvc_seri(i,k) = MIN(MAX(rvc_seri(i,k), 0.), 1.)
1090
1091        !--Diagnostics
1092        gamma_cond(i,k) = gammasat(i)
1093        subfra(i,k) = 1. - cf_seri(i,k) - issrfra(i,k)
1094        qsub(i,k) = zq(i) - qvc(i) - qissr(i,k)
1095        qcld(i,k) = qvc(i) + zoliq(i)
1096      ENDDO
1097    ENDIF
1098
1099    ! Outputs:
1100    !-------------------------------
1101    ! Precipitation fluxes at layer interfaces
1102    ! + precipitation fractions +
1103    ! temperature and water species tendencies
1104    DO i = 1, klon
1105        psfl(i,k)=zifl(i)
1106        prfl(i,k)=zrfl(i)
1107        pfraclr(i,k)=znebprecipclr(i)
1108        pfracld(i,k)=znebprecipcld(i)
1109        d_q(i,k) = zq(i) - qt(i,k)
1110        d_t(i,k) = zt(i) - temp(i,k)
1111
1112        IF (ok_bug_phase_lscp) THEN
1113           d_ql(i,k) = (1-zfice(i))*zoliq(i)
1114           d_qi(i,k) = zfice(i)*zoliq(i)
1115        ELSE
1116           d_ql(i,k) = zoliql(i)
1117           d_qi(i,k) = zoliqi(i)   
1118        ENDIF
1119
1120    ENDDO
1121
1122
1123ENDDO ! loop on k from top to bottom
1124
1125
1126  ! Rain or snow at the surface (depending on the first layer temperature)
1127  DO i = 1, klon
1128      snow(i) = zifl(i)
1129      rain(i) = zrfl(i)
1130      ! c_iso final output
1131  ENDDO
1132
1133  IF (ncoreczq>0) THEN
1134      WRITE(lunout,*)'WARNING : ZQ in LSCP ',ncoreczq,' val < 1.e-15.'
1135  ENDIF
1136
1137
1138RETURN
1139
1140END SUBROUTINE lscp
1141!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1142
1143END MODULE lmdz_lscp
Note: See TracBrowser for help on using the repository browser.