source: LMDZ6/branches/contrails/libf/phylmd/lmdz_lscp_tools.f90 @ 5644

Last change on this file since 5644 was 5609, checked in by aborella, 3 months ago
  • changed treatment of prognostic variables for prognostic clouds
  • adapted sedimentation and autoconversion for prognostic cirrus clouds
  • cloud mixing, ice sedimentation and ISSR diagnosis are now consistent with the water vapor PDF
  • simplified assumptions for ice crystals deposition / sublimation
  • first version of the coupling between prognostic cirrus clouds and deep convection
  • added persistent contrail cirrus clouds in radiative diagnostics
File size: 37.0 KB
Line 
1MODULE lmdz_lscp_tools
2
3    IMPLICIT NONE
4
5CONTAINS
6
7!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
8SUBROUTINE FALLICE_VELOCITY(klon,iwc,temp,rho,pres,ptconv,ptpronclds,velo)
9
10    ! Ref:
11    ! Stubenrauch, C. J., Bonazzola, M.,
12    ! Protopapadaki, S. E., & Musat, I. (2019).
13    ! New cloud system metrics to assess bulk
14    ! ice cloud schemes in a GCM. Journal of
15    ! Advances in Modeling Earth Systems, 11,
16    ! 3212–3234. https://doi.org/10.1029/2019MS001642
17   
18    use lmdz_lscp_ini, only: iflag_vice, ffallv_con, ffallv_lsc, ffallv_issr
19    use lmdz_lscp_ini, only: cice_velo, dice_velo
20    use lmdz_lscp_ini, only: ok_bug_ice_fallspeed, eps
21
22    IMPLICIT NONE
23
24    INTEGER, INTENT(IN) :: klon
25    REAL, INTENT(IN), DIMENSION(klon) :: iwc       ! specific ice water content [kg/m3]
26    REAL, INTENT(IN), DIMENSION(klon) :: temp      ! temperature [K]
27    REAL, INTENT(IN), DIMENSION(klon) :: rho       ! dry air density [kg/m3]
28    REAL, INTENT(IN), DIMENSION(klon) :: pres      ! air pressure [Pa]
29    LOGICAL, INTENT(IN), DIMENSION(klon) :: ptconv    ! convective point  [-]
30    LOGICAL, INTENT(IN), DIMENSION(klon) :: ptpronclds! prognostic clouds point  [-]
31
32    REAL, INTENT(OUT), DIMENSION(klon) :: velo    ! fallspeed velocity of crystals [m/s]
33
34
35    INTEGER i
36    REAL logvm,iwcg,tempc,phpa,fallv_tun
37    REAL m2ice, m2snow, vmice, vmsnow
38    REAL aice, bice, asnow, bsnow
39   
40
41    DO i=1,klon
42
43        IF (ptconv(i)) THEN
44            fallv_tun=ffallv_con
45        ELSEIF (ptpronclds(i)) THEN
46            fallv_tun=ffallv_issr
47        ELSE
48            fallv_tun=ffallv_lsc
49        ENDIF
50
51        tempc=temp(i)-273.15 ! celcius temp
52        IF ( ok_bug_ice_fallspeed ) THEN
53            iwcg=MAX(iwc(i)*1000.,1E-3) ! iwc in g/m3. We set a minimum value to prevent from division by 0
54        ELSE
55            ! AB the threshold is way too high, we reduce it
56            iwcg=MAX(iwc(i)*1000.,eps) ! iwc in g/m3. We set a minimum value to prevent from division by 0
57        ENDIF
58        phpa=pres(i)/100.    ! pressure in hPa
59
60    IF (iflag_vice .EQ. 1) THEN
61        ! so-called 'empirical parameterization' in Stubenrauch et al. 2019
62        if (tempc .GE. -60.0) then
63
64            logvm= -0.0000414122*tempc*tempc*log(iwcg)-0.00538922*tempc*log(iwcg) &
65                    -0.0516344*log(iwcg)+0.00216078*tempc + 1.9714   
66            velo(i)=exp(logvm)
67        else
68            velo(i)=65.0*(iwcg**0.2)*(150./phpa)**0.15
69        endif
70       
71        velo(i)=fallv_tun*velo(i)/100.0 ! from cm/s to m/s
72
73    ELSE IF (iflag_vice .EQ. 2) THEN
74        ! so called  PSDM empirical coherent bulk ice scheme in Stubenrauch et al. 2019
75        aice=0.587
76        bice=2.45
77        asnow=0.0444
78        bsnow=2.1
79       
80        m2ice=((iwcg*0.001/aice)/(exp(13.6-bice*7.76+0.479*bice**2)* &
81                exp((-0.0361+bice*0.0151+0.00149*bice**2)*tempc)))   &
82                **(1./(0.807+bice*0.00581+0.0457*bice**2))
83
84        vmice=100.*1042.4*exp(13.6-(bice+1)*7.76+0.479*(bice+1.)**2)*exp((-0.0361+ &
85                 (bice+1.)*0.0151+0.00149*(bice+1.)**2)*tempc) &
86                 *(m2ice**(0.807+(bice+1.)*0.00581+0.0457*(bice+1.)**2))/(iwcg*0.001/aice)
87
88       
89        vmice=vmice*((1000./phpa)**0.2)
90     
91        m2snow=((iwcg*0.001/asnow)/(exp(13.6-bsnow*7.76+0.479*bsnow**2)* &
92               exp((-0.0361+bsnow*0.0151+0.00149*bsnow**2)*tempc)))         &
93               **(1./(0.807+bsnow*0.00581+0.0457*bsnow**2))
94
95
96        vmsnow=100.*14.3*exp(13.6-(bsnow+.416)*7.76+0.479*(bsnow+.416)**2)&
97                  *exp((-0.0361+(bsnow+.416)*0.0151+0.00149*(bsnow+.416)**2)*tempc)&
98                  *(m2snow**(0.807+(bsnow+.416)*0.00581+0.0457*(bsnow+.416)**2))/(iwcg*0.001/asnow)
99
100        vmsnow=vmsnow*((1000./phpa)**0.35)
101        velo(i)=fallv_tun*min(vmsnow,vmice)/100. ! to m/s
102
103    ELSE
104        ! By default, fallspeed velocity of ice crystals according to Heymsfield & Donner 1990
105        velo(i) = fallv_tun*cice_velo*((iwcg/1000.)**dice_velo)
106    ENDIF
107    ENDDO
108
109END SUBROUTINE FALLICE_VELOCITY
110!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
111
112!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
113SUBROUTINE ICEFRAC_LSCP(klon, temp, iflag_ice_thermo, distcltop, temp_cltop, icefrac, dicefracdT)
114!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
115 
116  ! Compute the ice fraction 1-xliq (see e.g.
117  ! Doutriaux-Boucher & Quaas 2004, section 2.2.)
118  ! as a function of temperature
119  ! see also Fig 3 of Madeleine et al. 2020, JAMES
120!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
121
122
123    USE print_control_mod, ONLY: lunout, prt_level
124    USE lmdz_lscp_ini, ONLY: t_glace_min, t_glace_max, exposant_glace, iflag_t_glace
125    USE lmdz_lscp_ini, ONLY : RTT, dist_liq, temp_nowater
126
127    IMPLICIT NONE
128
129
130    INTEGER, INTENT(IN)                 :: klon              ! number of horizontal grid points
131    REAL, INTENT(IN), DIMENSION(klon)   :: temp              ! temperature
132    REAL, INTENT(IN), DIMENSION(klon)   :: distcltop         ! distance to cloud top
133    REAL, INTENT(IN), DIMENSION(klon)   :: temp_cltop        ! temperature of cloud top
134    INTEGER, INTENT(IN)                 :: iflag_ice_thermo
135    REAL, INTENT(OUT), DIMENSION(klon)  :: icefrac
136    REAL, INTENT(OUT), DIMENSION(klon)  :: dicefracdT
137
138
139    INTEGER i
140    REAL    liqfrac_tmp, dicefrac_tmp
141    REAL    Dv, denomdep,beta,qsi,dqsidt
142    LOGICAL ice_thermo
143
144    CHARACTER (len = 20) :: modname = 'lscp_tools'
145    CHARACTER (len = 80) :: abort_message
146
147    IF ((iflag_t_glace.LT.2)) THEN !.OR. (iflag_t_glace.GT.6)) THEN
148       abort_message = 'lscp cannot be used if iflag_t_glace<2 or >6'
149       CALL abort_physic(modname,abort_message,1)
150    ENDIF
151
152    IF (.NOT.((iflag_ice_thermo .EQ. 1).OR.(iflag_ice_thermo .GE. 3))) THEN
153       abort_message = 'lscp cannot be used without ice thermodynamics'
154       CALL abort_physic(modname,abort_message,1)
155    ENDIF
156
157
158    DO i=1,klon
159 
160        ! old function with sole dependence upon temperature
161        IF (iflag_t_glace .EQ. 2) THEN
162            liqfrac_tmp = (temp(i)-t_glace_min) / (t_glace_max-t_glace_min)
163            liqfrac_tmp = MIN(MAX(liqfrac_tmp,0.0),1.0)
164            icefrac(i) = (1.0-liqfrac_tmp)**exposant_glace
165            IF (icefrac(i) .GT.0.) THEN
166                 dicefracdT(i)= exposant_glace * (icefrac(i)**(exposant_glace-1.)) &
167                           / (t_glace_min - t_glace_max)
168            ENDIF
169
170            IF ((icefrac(i).EQ.0).OR.(icefrac(i).EQ.1)) THEN
171                 dicefracdT(i)=0.
172            ENDIF
173
174        ENDIF
175
176        ! function of temperature used in CMIP6 physics
177        IF (iflag_t_glace .EQ. 3) THEN
178            liqfrac_tmp = (temp(i)-t_glace_min) / (t_glace_max-t_glace_min)
179            liqfrac_tmp = MIN(MAX(liqfrac_tmp,0.0),1.0)
180            icefrac(i) = 1.0-liqfrac_tmp**exposant_glace
181            IF ((icefrac(i) .GT.0.) .AND. (liqfrac_tmp .GT. 0.)) THEN
182                dicefracdT(i)= exposant_glace * ((liqfrac_tmp)**(exposant_glace-1.)) &
183                          / (t_glace_min - t_glace_max)
184            ELSE
185                dicefracdT(i)=0.
186            ENDIF
187        ENDIF
188
189        ! for iflag_t_glace .GE. 4, the liquid fraction depends upon temperature at cloud top
190        ! and then decreases with decreasing height
191
192        !with linear function of temperature at cloud top
193        IF (iflag_t_glace .EQ. 4) THEN 
194                liqfrac_tmp = (temp(i)-t_glace_min) / (t_glace_max-t_glace_min)
195                liqfrac_tmp = MIN(MAX(liqfrac_tmp,0.0),1.0)
196                icefrac(i)    =  MAX(MIN(1.,1.0 - liqfrac_tmp*exp(-distcltop(i)/dist_liq)),0.)
197                dicefrac_tmp = - temp(i)/(t_glace_max-t_glace_min)
198                dicefracdT(i) = dicefrac_tmp*exp(-distcltop(i)/dist_liq)
199                IF ((liqfrac_tmp .LE.0) .OR. (liqfrac_tmp .GE. 1)) THEN
200                        dicefracdT(i) = 0.
201                ENDIF
202        ENDIF
203
204        ! with CMIP6 function of temperature at cloud top
205        IF ((iflag_t_glace .EQ. 5) .OR. (iflag_t_glace .EQ. 7)) THEN
206                liqfrac_tmp = (temp(i)-t_glace_min) / (t_glace_max-t_glace_min)
207                liqfrac_tmp =  MIN(MAX(liqfrac_tmp,0.0),1.0)
208                liqfrac_tmp = liqfrac_tmp**exposant_glace
209                icefrac(i)  =  MAX(MIN(1.,1.0 - liqfrac_tmp*exp(-distcltop(i)/dist_liq)),0.)
210                IF ((liqfrac_tmp .LE.0) .OR. (liqfrac_tmp .GE. 1)) THEN
211                        dicefracdT(i) = 0.
212                ELSE
213                        dicefracdT(i) = exposant_glace*((liqfrac_tmp)**(exposant_glace-1.))/(t_glace_min- t_glace_max) &
214                                        *exp(-distcltop(i)/dist_liq)
215                ENDIF
216        ENDIF
217
218        ! with modified function of temperature at cloud top
219        ! to get largere values around 260 K, works well with t_glace_min = 241K
220        IF (iflag_t_glace .EQ. 6) THEN
221                IF (temp(i) .GT. t_glace_max) THEN
222                        liqfrac_tmp = 1.
223                ELSE
224                        liqfrac_tmp = -((temp(i)-t_glace_max) / (t_glace_max-t_glace_min))**2+1.
225                ENDIF
226                liqfrac_tmp  = MIN(MAX(liqfrac_tmp,0.0),1.0)
227                icefrac(i)   = MAX(MIN(1.,1.0 - liqfrac_tmp*exp(-distcltop(i)/dist_liq)),0.)       
228                IF ((liqfrac_tmp .LE.0) .OR. (liqfrac_tmp .GE. 1)) THEN
229                        dicefracdT(i) = 0.
230                ELSE
231                        dicefracdT(i) = 2*((temp(i)-t_glace_max) / (t_glace_max-t_glace_min))/(t_glace_max-t_glace_min) &
232                                  *exp(-distcltop(i)/dist_liq)
233                ENDIF
234        ENDIF
235
236        ! if temperature of cloud top <-40°C,
237        IF (iflag_t_glace .GE. 4) THEN
238                IF ((temp_cltop(i) .LE. temp_nowater) .AND. (temp(i) .LE. t_glace_max)) THEN
239                        icefrac(i) = 1.
240                        dicefracdT(i) = 0.
241                ENDIF
242        ENDIF
243     
244
245     ENDDO ! klon
246     RETURN
247 
248END SUBROUTINE ICEFRAC_LSCP
249!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
250
251
252SUBROUTINE ICEFRAC_LSCP_TURB(klon, dtime, temp, pplay, paprsdn, paprsup, omega, qice_ini, snowcld, qtot_incl, cldfra, tke,   &
253                             tke_dissip, qliq, qvap_cld, qice, icefrac, dicefracdT, cldfraliq, sigma2_icefracturb, mean_icefracturb)
254!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
255  ! Compute the liquid, ice and vapour content (+ice fraction) based
256  ! on turbulence (see Fields 2014, Furtado 2016, Raillard 2025)
257  ! L.Raillard (23/09/24)
258!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
259
260
261   USE lmdz_lscp_ini, ONLY : prt_level, lunout
262   USE lmdz_lscp_ini, ONLY : RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT, RD, RG, RV, RPI
263   USE lmdz_lscp_ini, ONLY : seuil_neb, temp_nowater
264   USE lmdz_lscp_ini, ONLY : naero5, gamma_snwretro, gamma_taud, capa_crystal
265   USE lmdz_lscp_ini, ONLY : eps
266
267   IMPLICIT NONE
268
269   INTEGER,   INTENT(IN)                           :: klon                !--number of horizontal grid points
270   REAL,      INTENT(IN)                           :: dtime               !--time step [s]
271
272   REAL,      INTENT(IN),       DIMENSION(klon)    :: temp                !--temperature
273   REAL,      INTENT(IN),       DIMENSION(klon)    :: pplay               !--pressure in the middle of the layer           [Pa]
274   REAL,      INTENT(IN),       DIMENSION(klon)    :: paprsdn             !--pressure at the bottom interface of the layer [Pa]
275   REAL,      INTENT(IN),       DIMENSION(klon)    :: paprsup             !--pressure at the top interface of the layer    [Pa]
276   REAL,      INTENT(IN),       DIMENSION(klon)    :: omega               !--resolved vertical velocity                    [Pa/s]
277   REAL,      INTENT(IN),       DIMENSION(klon)    :: qtot_incl           !--specific total cloud water in-cloud content   [kg/kg]
278   REAL,      INTENT(IN),       DIMENSION(klon)    :: cldfra              !--cloud fraction in gridbox                     [-]
279   REAL,      INTENT(IN),       DIMENSION(klon)    :: tke                 !--turbulent kinetic energy                      [m2/s2]
280   REAL,      INTENT(IN),       DIMENSION(klon)    :: tke_dissip          !--TKE dissipation                               [m2/s3]
281
282   REAL,      INTENT(IN),       DIMENSION(klon)    :: qice_ini            !--initial specific ice content gridbox-mean     [kg/kg]
283   REAL,      INTENT(IN),       DIMENSION(klon)    :: snowcld
284   REAL,      INTENT(OUT),      DIMENSION(klon)    :: qliq                !--specific liquid content gridbox-mean          [kg/kg]
285   REAL,      INTENT(OUT),      DIMENSION(klon)    :: qvap_cld            !--specific cloud vapor content, gridbox-mean    [kg/kg]
286   REAL,      INTENT(OUT),      DIMENSION(klon)    :: qice                !--specific ice content gridbox-mean             [kg/kg]
287   REAL,      INTENT(OUT),      DIMENSION(klon)    :: icefrac             !--fraction of ice in condensed water            [-]
288   REAL,      INTENT(OUT),      DIMENSION(klon)    :: dicefracdT
289
290   REAL,      INTENT(OUT),      DIMENSION(klon)    :: cldfraliq           !--fraction of cldfra where liquid               [-]
291   REAL,      INTENT(OUT),      DIMENSION(klon)    :: sigma2_icefracturb  !--Sigma2 of the ice supersaturation PDF         [-]
292   REAL,      INTENT(OUT),      DIMENSION(klon)    :: mean_icefracturb    !--Mean of the ice supersaturation PDF           [-]
293
294   REAL, DIMENSION(klon) :: qzero, qsatl, dqsatl, qsati, dqsati           !--specific humidity saturation values
295   INTEGER :: i
296
297   REAL :: qvap_incl, qice_incl, qliq_incl, qiceini_incl                  !--In-cloud specific quantities                  [kg/kg]
298   REAL :: water_vapor_diff                                               !--Water-vapour diffusion coeff in air (f(T,P))  [m2/s]
299   REAL :: air_thermal_conduct                                            !--Thermal conductivity of air (f(T))            [J/m/K/s]
300   REAL :: C0                                                             !--Lagrangian structure function                 [-]
301   REAL :: tau_dissipturb
302   REAL :: tau_phaserelax
303   REAL :: sigma2_pdf, mean_pdf
304   REAL :: ai, bi, B0
305   REAL :: sursat_iceliq
306   REAL :: sursat_equ
307   REAL :: liqfra_max
308   REAL :: sursat_iceext
309   REAL :: nb_crystals                                                    !--number concentration of ice crystals          [#/m3]
310   REAL :: moment1_PSD                                                    !--1st moment of ice PSD
311   REAL :: N0_PSD, lambda_PSD                                             !--parameters of the exponential PSD
312
313   REAL :: rho_ice                                                        !--ice density                                   [kg/m3]
314   REAL :: cldfra1D
315   REAL :: rho_air
316   REAL :: psati                                                          !--saturation vapor pressure wrt ice             [Pa]
317   
318   REAL :: vitw                                                           !--vertical velocity                             [m/s]
319                                                                       
320   C0            = 10.                                                    !--value assumed in Field2014           
321   rho_ice       = 950.
322   sursat_iceext = -0.1
323   qzero(:)      = 0.
324   cldfraliq(:)  = 0.
325   icefrac(:)    = 0.
326   dicefracdT(:) = 0.
327
328   sigma2_icefracturb(:) = 0.
329   mean_icefracturb(:)   = 0.
330
331   !--wrt liquid
332   CALL calc_qsat_ecmwf(klon,temp(:),qzero(:),pplay(:),RTT,1,.false.,qsatl(:),dqsatl(:))
333   !--wrt ice
334   CALL calc_qsat_ecmwf(klon,temp(:),qzero(:),pplay(:),RTT,2,.false.,qsati(:),dqsati(:))
335
336
337    DO i=1,klon
338
339     rho_air  = pplay(i) / temp(i) / RD
340
341     ! because cldfra is intent in, but can be locally modified due to test
342     cldfra1D = cldfra(i)
343     IF (cldfra(i) .LE. 0.) THEN
344        qvap_cld(i)   = 0.
345        qliq(i)       = 0.
346        qice(i)       = 0.
347        cldfraliq(i)  = 0.
348        icefrac(i)    = 0.
349        dicefracdT(i) = 0.
350
351     ! If there is a cloud
352     ELSE
353        IF (cldfra(i) .GE. 1.0) THEN
354           cldfra1D = 1.0
355        END IF
356       
357        ! T>0°C, no ice allowed
358        IF ( temp(i) .GE. RTT ) THEN
359           qvap_cld(i)   = qsatl(i) * cldfra1D
360           qliq(i)       = MAX(0.0,qtot_incl(i)-qsatl(i))  * cldfra1D
361           qice(i)       = 0.
362           cldfraliq(i)  = 1.
363           icefrac(i)    = 0.
364           dicefracdT(i) = 0.
365       
366        ! T<-38°C, no liquid allowed
367        ELSE IF ( temp(i) .LE. temp_nowater) THEN
368           qvap_cld(i)   = qsati(i) * cldfra1D
369           qliq(i)       = 0.
370           qice(i)       = MAX(0.0,qtot_incl(i)-qsati(i)) * cldfra1D
371           cldfraliq(i)  = 0.
372           icefrac(i)    = 1.
373           dicefracdT(i) = 0.
374
375        !---------------------------------------------------------
376        !--             MIXED PHASE TEMPERATURE REGIME         
377        !---------------------------------------------------------
378        !--In the mixed phase regime (-38°C< T <0°C) we distinguish
379        !--3 possible subcases.
380        !--1.  No pre-existing ice 
381        !--2A. Pre-existing ice and no turbulence
382        !--2B. Pre-existing ice and turbulence
383        ELSE
384
385           vitw = -omega(i) / RG / rho_air
386           qiceini_incl  = qice_ini(i) / cldfra1D + snowcld(i) * RG * dtime / ( paprsdn(i) - paprsup(i) ) / cldfra1D
387
388           !--1. No preexisting ice : if vertical motion, fully liquid
389           !--cloud else fully iced cloud
390           IF ( qiceini_incl .LT. eps ) THEN
391              IF ( (vitw .GT. eps) .OR. (tke(i) .GT. eps) ) THEN
392                 qvap_cld(i)   = qsatl(i) * cldfra1D
393                 qliq(i)       = MAX(0.,qtot_incl(i)-qsatl(i)) * cldfra1D
394                 qice(i)       = 0.
395                 cldfraliq(i)  = 1.
396                 icefrac(i)    = 0.
397                 dicefracdT(i) = 0.
398              ELSE
399                 qvap_cld(i)   = qsati(i) * cldfra1D
400                 qliq(i)       = 0.
401                 qice(i)       = MAX(0.,qtot_incl(i)-qsati(i)) * cldfra1D
402                 cldfraliq(i)  = 0.
403                 icefrac(i)    = 1.
404                 dicefracdT(i) = 0.
405              ENDIF
406           
407
408           !--2. Pre-existing ice :computation of ice properties for
409           !--feedback
410           ELSE
411              ai = RG / RD / temp(i) * ( RD * RLSTT / RCPD / RV / temp(i) - 1. )
412
413              sursat_equ    = ai * vitw * tau_phaserelax
414
415              sursat_iceliq = qsatl(i)/qsati(i) - 1.
416              psati         = qsati(i) * pplay(i) / (RD/RV)
417             
418              !--We assume an exponential ice PSD whose parameters
419              !--are computed following Morrison&Gettelman 2008
420              !--Ice number density is assumed equals to INP density
421              !--which is a function of temperature (DeMott 2010) 
422              !--bi and B0 are microphysical function characterizing
423              !--vapor/ice interactions
424              !--tau_phase_relax is the typical time of vapor deposition
425              !--onto ice crystals
426
427              nb_crystals = 1.e3 * 5.94e-5 * ( RTT - temp(i) )**3.33 * naero5**(0.0264*(RTT-temp(i))+0.0033)
428              lambda_PSD  = ( (RPI*rho_ice*nb_crystals) / (rho_air * qiceini_incl ) ) ** (1./3.)
429              N0_PSD      = nb_crystals * lambda_PSD
430              moment1_PSD = N0_PSD/lambda_PSD**2
431
432              !--Formulae for air thermal conductivity and water vapor diffusivity
433              !--comes respectively from Beard and Pruppacher (1971)
434              !--and  Hall and Pruppacher (1976)
435
436              air_thermal_conduct = ( 5.69 + 0.017 * ( temp(i) - RTT ) ) * 1.e-3 * 4.184
437              water_vapor_diff    = 2.11*1e-5 * ( temp(i) / RTT )**1.94 * ( 101325 / pplay(i) )
438             
439              bi = 1./((qsati(i)+qsatl(i))/2.) + RLSTT**2 / RCPD / RV / temp(i)**2
440              B0 = 4. * RPI * capa_crystal * 1. / (  RLSTT**2 / air_thermal_conduct / RV / temp(i)**2  &
441                                                  +  RV * temp(i) / psati / water_vapor_diff  )
442              tau_phaserelax = 1. / (bi * B0 * moment1_PSD )
443             
444              ai = RG / RD / temp(i) * ( RD * RLSTT / RCPD / RV / temp(i) - 1. )
445
446              !--2A. No TKE : stationnary binary solution depending on omega
447              ! If Sequ > Siw liquid cloud, else ice cloud
448              IF ( tke_dissip(i) .LE. eps )  THEN
449                 IF (sursat_equ .GT. sursat_iceliq) THEN
450                    qvap_cld(i)   = qsatl(i) * cldfra1D
451                    qliq(i)       = MAX(0.,qtot_incl(i)-qsatl(i)) * cldfra1D
452                    qice(i)       = 0.
453                    cldfraliq(i)  = 1.
454                    icefrac(i)    = 0.
455                    dicefracdT(i) = 0.
456                 ELSE
457                    qvap_cld(i)   = qsati(i) * cldfra1D
458                    qliq(i)       = 0.
459                    qice(i)       = MAX(0.,qtot_incl(i)-qsati(i)) * cldfra1D
460                    cldfraliq(i)  = 0.
461                    icefrac(i)    = 1.
462                    dicefracdT(i) = 0.
463                 ENDIF
464                 
465              !--2B. TKE and ice : ice supersaturation PDF
466              !--we compute the cloud liquid properties with a Gaussian PDF
467              !--of ice supersaturation F(Si) (Field2014, Furtado2016).
468              !--Parameters of the PDF are function of turbulence and
469              !--microphysics/existing ice.
470              ELSE 
471                     
472                 !--Tau_dissipturb is the time needed for turbulence to decay
473                 !--due to viscosity
474                 tau_dissipturb = gamma_taud * 2. * 2./3. * tke(i) / tke_dissip(i) / C0
475
476                 !--------------------- PDF COMPUTATIONS ---------------------
477                 !--Formulae for sigma2_pdf (variance), mean of PDF in Raillard2025
478                 !--cloud liquid fraction and in-cloud liquid content are given
479                 !--by integrating resp. F(Si) and Si*F(Si)
480                 !--Liquid is limited by the available water vapor trough a
481                 !--maximal liquid fraction
482                 !--qice_ini(i) / cldfra1D = qiceincld without precip
483
484                 liqfra_max = MAX(0., (MIN (1.,( qtot_incl(i) - (qice_ini(i) / cldfra1D) - qsati(i) * (1 + sursat_iceext ) ) / ( qsatl(i) - qsati(i) ) ) ) )
485                 sigma2_pdf = 1./2. * ( ai**2 ) *  2./3. * tke(i) * tau_dissipturb * tau_phaserelax
486                 
487                 mean_pdf = ai * vitw * tau_phaserelax
488                 
489                 cldfraliq(i) = 0.5 * (1. - erf( ( sursat_iceliq - mean_pdf) / (SQRT(2.* sigma2_pdf) ) ) )
490                 IF (cldfraliq(i) .GT. liqfra_max) THEN
491                     cldfraliq(i) = liqfra_max
492                 ENDIF
493                 
494                 qliq_incl = qsati(i) * SQRT(sigma2_pdf) / SQRT(2.*RPI) * EXP( -1.*(sursat_iceliq - mean_pdf)**2. / (2.*sigma2_pdf) )  &
495                           - qsati(i) * cldfraliq(i) * (sursat_iceliq - mean_pdf )
496                 
497                 sigma2_icefracturb(i)= sigma2_pdf
498                 mean_icefracturb(i)  = mean_pdf
499     
500                 !------------ SPECIFIC VAPOR CONTENT AND WATER CONSERVATION  ------------
501
502                 IF ( (qliq_incl .LE. eps) .OR. (cldfraliq(i) .LE. eps) ) THEN
503                     qliq_incl    = 0.
504                     cldfraliq(i) = 0.
505                 END IF
506                 
507                 !--Specific humidity is the max between qsati and the weighted mean between
508                 !--qv in MPC patches and qv in ice-only parts. We assume that MPC parts are
509                 !--always at qsatl and ice-only parts slightly subsaturated (qsati*sursat_iceext+1)
510                 !--The whole cloud can therefore be supersaturated but never subsaturated.
511
512                 qvap_incl = MAX(qsati(i), ( 1. - cldfraliq(i) ) * (sursat_iceext + 1.) * qsati(i) + cldfraliq(i) * qsatl(i) )
513
514                 IF ( qvap_incl  .GE. qtot_incl(i) ) THEN
515                    qvap_incl = qsati(i)
516                    qliq_incl = qtot_incl(i) - qvap_incl
517                    qice_incl = 0.
518
519                 ELSEIF ( (qvap_incl + qliq_incl) .GE. qtot_incl(i) ) THEN
520                    qliq_incl = MAX(0.0,qtot_incl(i) - qvap_incl)
521                    qice_incl = 0.
522                 ELSE
523                    qice_incl = qtot_incl(i) - qvap_incl - qliq_incl
524                 END IF
525
526                 qvap_cld(i)   = qvap_incl * cldfra1D
527                 qliq(i)       = qliq_incl * cldfra1D
528                 qice(i)       = qice_incl * cldfra1D
529                 icefrac(i)    = qice(i) / ( qice(i) + qliq(i) )
530                 dicefracdT(i) = 0.
531
532              END IF ! Enough TKE
533           
534           END IF ! End qini
535
536        END IF ! ! MPC temperature
537
538     END IF ! cldfra
539   
540   ENDDO ! klon
541END SUBROUTINE ICEFRAC_LSCP_TURB
542!
543!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
544
545
546SUBROUTINE CALC_QSAT_ECMWF(klon,temp,qtot,pressure,tref,phase,flagth,qs,dqs)
547!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
548    ! Calculate qsat following ECMWF method
549!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
550
551
552USE yoethf_mod_h
553        USE yomcst_mod_h
554IMPLICIT NONE
555
556
557    include "FCTTRE.h"
558
559    INTEGER, INTENT(IN) :: klon  ! number of horizontal grid points
560    REAL, INTENT(IN), DIMENSION(klon) :: temp     ! temperature in K
561    REAL, INTENT(IN), DIMENSION(klon) :: qtot     ! total specific water in kg/kg
562    REAL, INTENT(IN), DIMENSION(klon) :: pressure ! pressure in Pa
563    REAL, INTENT(IN)                  :: tref     ! reference temperature in K
564    LOGICAL, INTENT(IN) :: flagth     ! flag for qsat calculation for thermals
565    INTEGER, INTENT(IN) :: phase
566    ! phase: 0=depend on temperature sign (temp>tref -> liquid, temp<tref, solid)
567    !        1=liquid
568    !        2=solid
569
570    REAL, INTENT(OUT), DIMENSION(klon) :: qs      ! saturation specific humidity [kg/kg]
571    REAL, INTENT(OUT), DIMENSION(klon) :: dqs     ! derivation of saturation specific humidity wrt T
572
573    REAL delta, cor, cvm5
574    INTEGER i
575
576    DO i=1,klon
577
578    IF (phase .EQ. 1) THEN
579        delta=0.
580    ELSEIF (phase .EQ. 2) THEN
581        delta=1.
582    ELSE
583        delta=MAX(0.,SIGN(1.,tref-temp(i)))
584    ENDIF
585
586    IF (flagth) THEN
587    cvm5=R5LES*(1.-delta) + R5IES*delta
588    ELSE
589    cvm5 = R5LES*RLVTT*(1.-delta) + R5IES*RLSTT*delta
590    cvm5 = cvm5 /RCPD/(1.0+RVTMP2*(qtot(i)))
591    ENDIF
592
593    qs(i)= R2ES*FOEEW(temp(i),delta)/pressure(i)
594    qs(i)=MIN(0.5,qs(i))
595    cor=1./(1.-RETV*qs(i))
596    qs(i)=qs(i)*cor
597    dqs(i)= FOEDE(temp(i),delta,cvm5,qs(i),cor)
598
599    END DO
600
601END SUBROUTINE CALC_QSAT_ECMWF
602!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
603
604
605!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
606SUBROUTINE CALC_GAMMASAT(klon,temp,qtot,pressure,gammasat,dgammasatdt)
607
608!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
609    ! programme that calculates the gammasat parameter that determines the
610    ! homogeneous condensation thresholds for cold (<0oC) clouds
611    ! condensation at q>gammasat*qsat
612    ! Etienne Vignon, March 2021
613!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
614
615    use lmdz_lscp_ini, only: iflag_gammasat, temp_nowater, RTT
616    use lmdz_lscp_ini, only: a_homofreez, b_homofreez, delta_hetfreez
617
618    IMPLICIT NONE
619
620
621    INTEGER, INTENT(IN) :: klon                       ! number of horizontal grid points
622    REAL, INTENT(IN), DIMENSION(klon) :: temp         ! temperature in K
623    REAL, INTENT(IN), DIMENSION(klon) :: qtot         ! total specific water in kg/kg
624
625    REAL, INTENT(IN), DIMENSION(klon) :: pressure     ! pressure in Pa
626
627    REAL, INTENT(OUT), DIMENSION(klon) :: gammasat    ! coefficient to multiply qsat with to calculate saturation
628    REAL, INTENT(OUT), DIMENSION(klon) :: dgammasatdt ! derivative of gammasat wrt temperature
629
630    REAL, DIMENSION(klon) ::  qsi,qsl,dqsl,dqsi
631    REAL  f_homofreez, fac
632
633    INTEGER i
634   
635        CALL CALC_QSAT_ECMWF(klon,temp,qtot,pressure,RTT,1,.false.,qsl,dqsl)
636        CALL CALC_QSAT_ECMWF(klon,temp,qtot,pressure,RTT,2,.false.,qsi,dqsi)
637
638    DO i = 1, klon
639
640        IF ( temp(i) .GE. RTT ) THEN
641            ! warm clouds: condensation at saturation wrt liquid
642            gammasat(i) = 1.
643            dgammasatdt(i) = 0.
644
645        ELSE
646            ! cold clouds: qsi > qsl
647           
648            ! homogeneous freezing of aerosols, according to
649            ! Koop, 2000 and Ren and MacKenzie, 2005 (QJRMS)
650            ! 'Cirrus regime'
651            ! if f_homofreez > qsl / qsi, liquid nucleation
652            ! if f_homofreez < qsl / qsi, homogeneous freezing of aerosols
653            ! Note: f_homofreez = qsl / qsi for temp ~= -38degC
654            f_homofreez = a_homofreez - temp(i) / b_homofreez
655           
656            IF ( iflag_gammasat .GE. 3 ) THEN
657              ! condensation at homogeneous freezing threshold for temp < -38 degC
658              ! condensation at liquid saturation for temp > -38 degC
659              IF ( f_homofreez .LE. qsl(i) / qsi(i) ) THEN
660                gammasat(i) = f_homofreez
661                dgammasatdt(i) = - 1. / b_homofreez
662              ELSE
663                gammasat(i) = qsl(i) / qsi(i)
664                dgammasatdt(i) = ( dqsl(i) * qsi(i) - dqsi(i) * qsl(i) ) / qsi(i) / qsi(i)
665              ENDIF
666
667            ELSEIF ( iflag_gammasat .EQ. 2 ) THEN
668              ! condensation at homogeneous freezing threshold for temp < -38 degC
669              ! condensation at a threshold linearly decreasing between homogeneous
670              ! freezing and ice saturation for -38 degC < temp < temp_nowater
671              ! condensation at ice saturation for temp > temp_nowater
672              ! If temp_nowater = 235.15 K, this is equivalent to iflag_gammasat = 1
673              IF ( f_homofreez .LE. qsl(i) / qsi(i) ) THEN
674                gammasat(i) = f_homofreez
675                dgammasatdt(i) = - 1. / b_homofreez
676              ELSEIF ( temp(i) .LE. temp_nowater ) THEN
677                ! Here, we assume that f_homofreez = qsl / qsi for temp = -38 degC = 235.15 K
678                dgammasatdt(i) = ( a_homofreez - 235.15 / b_homofreez - 1. ) &
679                               / ( 235.15 - temp_nowater )
680                gammasat(i) = dgammasatdt(i) * ( temp(i) - temp_nowater ) + 1.
681              ELSE
682                gammasat(i) = 1.
683                dgammasatdt(i) = 0.
684              ENDIF
685
686            ELSEIF ( iflag_gammasat .EQ. 1 ) THEN
687              ! condensation at homogeneous freezing threshold for temp < -38 degC
688              ! condensation at ice saturation for temp > -38 degC
689              IF ( f_homofreez .LE. qsl(i) / qsi(i) ) THEN
690                gammasat(i) = f_homofreez
691                dgammasatdt(i) = - 1. / b_homofreez
692              ELSE
693                gammasat(i) = 1.
694                dgammasatdt(i) = 0.
695              ENDIF
696
697            ELSE
698              ! condensation at ice saturation for temp < -38 degC
699              ! condensation at ice saturation for temp > -38 degC
700              gammasat(i) = 1.
701              dgammasatdt(i) = 0.
702
703            ENDIF
704
705            ! Note that the delta_hetfreez parameter allows to linearly decrease the
706            ! condensation threshold between the calculated threshold and the ice saturation
707            ! for delta_hetfreez = 1, the threshold is the calculated condensation threshold
708            ! for delta_hetfreez = 0, the threshold is the ice saturation
709            gammasat(i) = ( 1. - delta_hetfreez ) + delta_hetfreez * gammasat(i)
710            dgammasatdt(i) = delta_hetfreez * dgammasatdt(i)
711
712        ENDIF
713   
714    END DO
715
716
717END SUBROUTINE CALC_GAMMASAT
718!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
719
720
721!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
722SUBROUTINE DISTANCE_TO_CLOUD_TOP(klon,klev,k,temp,pplay,paprs,rneb,distcltop1D,temp_cltop)
723!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
724   
725   USE lmdz_lscp_ini, ONLY : rd,rg,tresh_cl
726
727   IMPLICIT NONE
728   
729   INTEGER, INTENT(IN) :: klon,klev                !number of horizontal and vertical grid points
730   INTEGER, INTENT(IN) :: k                        ! vertical index
731   REAL, INTENT(IN), DIMENSION(klon,klev) :: temp  ! temperature in K
732   REAL, INTENT(IN), DIMENSION(klon,klev) :: pplay ! pressure middle layer in Pa
733   REAL, INTENT(IN), DIMENSION(klon,klev+1) :: paprs ! pressure interfaces in Pa
734   REAL, INTENT(IN), DIMENSION(klon,klev) :: rneb  ! cloud fraction
735
736   REAL, INTENT(OUT), DIMENSION(klon) :: distcltop1D  ! distance from cloud top
737   REAL, INTENT(OUT), DIMENSION(klon) :: temp_cltop     ! temperature of cloud top
738   
739   REAL dzlay(klon,klev)
740   REAL zlay(klon,klev)
741   REAL dzinterf
742   INTEGER i,k_top, kvert
743   LOGICAL bool_cl
744
745
746   DO i=1,klon
747         ! Initialization height middle of first layer
748          dzlay(i,1) = Rd * temp(i,1) / rg * log(paprs(i,1)/paprs(i,2))
749          zlay(i,1) = dzlay(i,1)/2
750
751          DO kvert=2,klev
752                 IF (kvert.EQ.klev) THEN
753                       dzlay(i,kvert) = 2*(rd * temp(i,kvert) / rg * log(paprs(i,kvert)/pplay(i,kvert)))
754                 ELSE
755                       dzlay(i,kvert) = rd * temp(i,kvert) / rg * log(paprs(i,kvert)/paprs(i,kvert+1))
756                 ENDIF
757                       dzinterf       = rd * temp(i,kvert) / rg * log(pplay(i,kvert-1)/pplay(i,kvert))
758                       zlay(i,kvert)  = zlay(i,kvert-1) + dzinterf
759           ENDDO
760   ENDDO
761   
762   DO i=1,klon
763          k_top = k
764          IF (rneb(i,k) .LE. tresh_cl) THEN
765                 bool_cl = .FALSE.
766          ELSE
767                 bool_cl = .TRUE.
768          ENDIF
769
770          DO WHILE ((bool_cl) .AND. (k_top .LE. klev))
771          ! find cloud top
772                IF (rneb(i,k_top) .GT. tresh_cl) THEN
773                      k_top = k_top + 1
774                ELSE
775                      bool_cl = .FALSE.
776                      k_top   = k_top - 1
777                ENDIF
778          ENDDO
779          k_top=min(k_top,klev)
780
781          !dist to top is dist between current layer and layer of cloud top (from middle to middle) + dist middle to
782          !interf for layer of cloud top
783          distcltop1D(i) = zlay(i,k_top) - zlay(i,k) + dzlay(i,k_top)/2
784          temp_cltop(i)  = temp(i,k_top)
785   ENDDO ! klon
786
787END SUBROUTINE DISTANCE_TO_CLOUD_TOP
788!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
789
790!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
791FUNCTION GAMMAINC ( p, x )
792
793!*****************************************************************************80
794!
795!! GAMMAINC computes the regularized lower incomplete Gamma Integral
796!
797!  Modified:
798!
799!    20 January 2008
800!
801!  Author:
802!
803!    Original FORTRAN77 version by B Shea.
804!    FORTRAN90 version by John Burkardt.
805!
806!  Reference:
807!
808!    B Shea,
809!    Algorithm AS 239:
810!    Chi-squared and Incomplete Gamma Integral,
811!    Applied Statistics,
812!    Volume 37, Number 3, 1988, pages 466-473.
813!
814!  Parameters:
815!
816!    Input, real X, P, the parameters of the incomplete
817!    gamma ratio.  0 <= X, and 0 < P.
818!
819!    Output, real GAMMAINC, the value of the incomplete
820!    Gamma integral.
821!
822  IMPLICIT NONE
823
824  REAL A
825  REAL AN
826  REAL ARG
827  REAL B
828  REAL C
829  REAL, PARAMETER :: ELIMIT = - 88.0E+00
830  REAL GAMMAINC
831  REAL, PARAMETER :: OFLO = 1.0E+37
832  REAL P
833  REAL, PARAMETER :: PLIMIT = 1000.0E+00
834  REAL PN1
835  REAL PN2
836  REAL PN3
837  REAL PN4
838  REAL PN5
839  REAL PN6
840  REAL RN
841  REAL, PARAMETER :: TOL = 1.0E-14
842  REAL X
843  REAL, PARAMETER :: XBIG = 1.0E+08
844
845  GAMMAINC = 0.0E+00
846
847  IF ( X == 0.0E+00 ) THEN
848    GAMMAINC = 0.0E+00
849    RETURN
850  END IF
851!
852!  IF P IS LARGE, USE A NORMAL APPROXIMATION.
853!
854  IF ( PLIMIT < P ) THEN
855
856    PN1 = 3.0E+00 * SQRT ( P ) * ( ( X / P )**( 1.0E+00 / 3.0E+00 ) &
857    + 1.0E+00 / ( 9.0E+00 * P ) - 1.0E+00 )
858
859    GAMMAINC = 0.5E+00 * ( 1. + ERF ( PN1 ) )
860    RETURN
861
862  END IF
863!
864!  IF X IS LARGE SET GAMMAD = 1.
865!
866  IF ( XBIG < X ) THEN
867    GAMMAINC = 1.0E+00
868    RETURN
869  END IF
870!
871!  USE PEARSON'S SERIES EXPANSION.
872!  (NOTE THAT P IS NOT LARGE ENOUGH TO FORCE OVERFLOW IN ALOGAM).
873!
874  IF ( X <= 1.0E+00 .OR. X < P ) THEN
875
876    ARG = P * LOG ( X ) - X - LOG_GAMMA ( P + 1.0E+00 )
877    C = 1.0E+00
878    GAMMAINC = 1.0E+00
879    A = P
880
881    DO
882
883      A = A + 1.0E+00
884      C = C * X / A
885      GAMMAINC = GAMMAINC + C
886
887      IF ( C <= TOL ) THEN
888        EXIT
889      END IF
890
891    END DO
892
893    ARG = ARG + LOG ( GAMMAINC )
894
895    IF ( ELIMIT <= ARG ) THEN
896      GAMMAINC = EXP ( ARG )
897    ELSE
898      GAMMAINC = 0.0E+00
899    END IF
900!
901!  USE A CONTINUED FRACTION EXPANSION.
902!
903  ELSE
904
905    ARG = P * LOG ( X ) - X - LOG_GAMMA ( P )
906    A = 1.0E+00 - P
907    B = A + X + 1.0E+00
908    C = 0.0E+00
909    PN1 = 1.0E+00
910    PN2 = X
911    PN3 = X + 1.0E+00
912    PN4 = X * B
913    GAMMAINC = PN3 / PN4
914
915    DO
916
917      A = A + 1.0E+00
918      B = B + 2.0E+00
919      C = C + 1.0E+00
920      AN = A * C
921      PN5 = B * PN3 - AN * PN1
922      PN6 = B * PN4 - AN * PN2
923
924      IF ( PN6 /= 0.0E+00 ) THEN
925
926        RN = PN5 / PN6
927
928        IF ( ABS ( GAMMAINC - RN ) <= MIN ( TOL, TOL * RN ) ) THEN
929          EXIT
930        END IF
931
932        GAMMAINC = RN
933
934      END IF
935
936      PN1 = PN3
937      PN2 = PN4
938      PN3 = PN5
939      PN4 = PN6
940!
941!  RE-SCALE TERMS IN CONTINUED FRACTION IF TERMS ARE LARGE.
942!
943      IF ( OFLO <= ABS ( PN5 ) ) THEN
944        PN1 = PN1 / OFLO
945        PN2 = PN2 / OFLO
946        PN3 = PN3 / OFLO
947        PN4 = PN4 / OFLO
948      END IF
949
950    END DO
951
952    ARG = ARG + LOG ( GAMMAINC )
953
954    IF ( ELIMIT <= ARG ) THEN
955      GAMMAINC = 1.0E+00 - EXP ( ARG )
956    ELSE
957      GAMMAINC = 1.0E+00
958    END IF
959
960  END IF
961
962  RETURN
963END FUNCTION GAMMAINC
964!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
965
966END MODULE lmdz_lscp_tools
967
968
Note: See TracBrowser for help on using the repository browser.