source: LMDZ6/branches/Amaury_dev/libf/phylmd/lmdz_lscp_tools.F90 @ 5411

Last change on this file since 5411 was 5224, checked in by abarral, 3 months ago

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