source: LMDZ6/trunk/libf/phylmd/lmdz_lscp_tools.F90 @ 5229

Last change on this file since 5229 was 5204, checked in by Laurent Fairhead, 8 weeks ago

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