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

Last change on this file since 5451 was 5443, checked in by evignon, 4 days ago

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