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

Last change on this file since 5453 was 5453, checked in by aborella, 21 hours ago

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