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

Last change on this file since 5456 was 5456, checked in by aborella, 37 hours ago

Added diagnostics for contrails fraction

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