source: LMDZ6/trunk/libf/phylmd/lscp_tools_mod.F90 @ 4062

Last change on this file since 4062 was 4059, checked in by oboucher, 3 years ago

Audran Borella's parametrisation for ice supersaturation
activated with flag_ice_sursat (FALSE by default)
File size: 10.5 KB
Line 
1MODULE LSCP_TOOLS_MOD
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   
18   
19    IMPLICIT NONE
20
21    INCLUDE "nuage.h"
22    INCLUDE "fisrtilp.h"
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
31    REAL, INTENT(OUT), DIMENSION(klon) :: velo    ! fallspeed velocity of crystals [m/s]
32
33
34    INTEGER i
35    REAL logvm,iwcg,tempc,phpa,cvel,dvel,fallv_tun
36    REAL m2ice, m2snow, vmice, vmsnow
37    REAL aice, bice, asnow, bsnow
38   
39
40    DO i=1,klon
41
42        IF (ptconv(i)) THEN
43            fallv_tun=ffallv_con
44        ELSE
45            fallv_tun=ffallv_lsc
46        ENDIF
47
48        tempc=temp(i)-273.15 ! celcius temp
49        iwcg=iwc(i)*1000.    ! iwc in g/m3
50        phpa=pres(i)/100.    ! pressure in hPa
51
52    IF (iflag_vice .EQ. 1) THEN
53        ! so-called 'empirical parameterization' in Stubenrauch et al. 2019
54        if (tempc .GE. -60.0) then
55
56            logvm= -0.0000414122*tempc*tempc*log(iwcg)-0.00538922*tempc*log(iwcg) &
57                    -0.0516344*log(iwcg)+0.00216078*tempc + 1.9714   
58            velo(i)=exp(logvm)
59        else
60            velo(i)=65.0*(iwcg**0.2)*(150./phpa)**0.15
61        endif
62       
63        velo(i)=fallv_tun*velo(i)/100.0 ! from cm/s to m/s
64        dvel=0.2
65        cvel=fallv_tun*65.0*(rho(i)**0.2)*(150./phpa)**0.15
66
67    ELSE IF (iflag_vice .EQ. 2) THEN
68        ! so called  PSDM empirical coherent bulk ice scheme in Stubenrauch et al. 2019
69        aice=0.587
70        bice=2.45
71        asnow=0.0444
72        bsnow=2.1
73       
74        m2ice=((iwcg*0.001/aice)/(exp(13.6-bice*7.76+0.479*bice**2)* &
75                exp((-0.0361+bice*0.0151+0.00149*bice**2)*tempc)))   &
76                **(1./(0.807+bice*0.00581+0.0457*bice**2))
77
78        vmice=100.*1042.4*exp(13.6-(bice+1)*7.76+0.479*(bice+1.)**2)*exp((-0.0361+&
79                 (bice+1.)*0.0151+0.00149*(bice+1.)**2)*tempc)&
80                *(m2ice**(0.807+(bice+1.)*0.00581+0.0457*(bice+1.)**2))/(iwcg*0.001/aice)
81
82       
83        vmice=vmice*((1000./phpa)**0.2)
84     
85        m2snow=((iwcg*0.001/asnow)/(exp(13.6-bsnow*7.76+0.479*bsnow**2)* &
86               exp((-0.0361+bsnow*0.0151+0.00149*bsnow**2)*tempc)))         &
87               **(1./(0.807+bsnow*0.00581+0.0457*bsnow**2))
88
89
90        vmsnow=100.*14.3*exp(13.6-(bsnow+.416)*7.76+0.479*(bsnow+.416)**2)&
91                  *exp((-0.0361+(bsnow+.416)*0.0151+0.00149*(bsnow+.416)**2)*tempc)&
92                  *(m2snow**(0.807+(bsnow+.416)*0.00581+0.0457*(bsnow+.416)**2))/(iwcg*0.001/asnow)
93
94        vmsnow=vmsnow*((1000./phpa)**0.35)
95
96        velo(i)=fallv_tun*min(vmsnow,vmice)/100. ! to m/s
97        dvel=0.2
98        cvel=velo(i)/((iwc(i)*rho(i))**dvel)
99
100    ELSE
101        ! By default, fallspeed velocity of ice crystals according to Heymsfield & Donner 1990
102        velo(i) = fallv_tun*3.29/2.0 * ((iwc(i))**0.16)
103        dvel=0.16
104        cvel=fallv_tun*3.29/2.0*(rho(i)**0.16)
105    ENDIF
106
107    ENDDO
108
109END SUBROUTINE FALLICE_VELOCITY
110!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
111
112!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
113SUBROUTINE ICEFRAC_LSCP(klon, temp, sig, 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
125    IMPLICIT NONE
126
127
128    INCLUDE "YOMCST.h"
129    INCLUDE "nuage.h"
130    INCLUDE "clesphys.h"
131
132
133  ! nuage.h contains:
134  ! t_glace_min: if T < Tmin, the cloud is only made of water ice
135  ! t_glace_max: if T > Tmax, the cloud is only made of liquid water
136  ! exposant_glace: controls the sharpness of the transition
137
138    INTEGER, INTENT(IN)                 :: klon       ! number of horizontal grid points
139    REAL, INTENT(IN), DIMENSION(klon)   :: temp       ! temperature
140    REAL, INTENT(IN), DIMENSION(klon)   :: sig
141    REAL, INTENT(OUT), DIMENSION(klon)  :: icefrac
142    REAL, INTENT(OUT), DIMENSION(klon)  :: dicefracdT
143
144
145    INTEGER i
146    REAL    sig0,www,tmin_tmp,liqfrac_tmp
147    REAL    Dv, denomdep,beta,qsi,dqsidt
148    INTEGER exposant_glace_old
149    REAL t_glace_min_old
150    LOGICAL ice_thermo
151
152    sig0=0.8
153    t_glace_min_old = RTT - 15.0
154    ice_thermo = (iflag_ice_thermo .EQ. 1).OR.(iflag_ice_thermo .GE. 3)
155    IF (ice_thermo) THEN
156      exposant_glace_old = 2
157    ELSE
158      exposant_glace_old = 6
159    ENDIF
160
161
162! calculation of icefrac and dicefrac/dT
163
164    DO i=1,klon
165
166    IF (iflag_t_glace.EQ.1) THEN
167            ! Transition to ice close to surface for T<Tmax
168            ! w=1 at the surface and 0 for sig < sig0
169            www=(max(sig(i)-sig0,0.))/(1.-sig0)
170    ELSEIF (iflag_t_glace.GE.2) THEN
171            ! No convertion to ice close to surface
172            www = 0.
173    ENDIF
174     
175    tmin_tmp=www*t_glace_max+(1.-www)*t_glace_min
176    liqfrac_tmp=  (temp(i)-tmin_tmp) / (t_glace_max-tmin_tmp)
177    liqfrac_tmp = MIN(MAX(liqfrac_tmp,0.0),1.0)
178       
179    IF (iflag_t_glace.GE.3) THEN
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
186                 dicefracdT(i)=0.
187            ENDIF
188
189    ELSE
190            icefrac(i) = (1.0-liqfrac_tmp)**exposant_glace
191            IF (icefrac(i) .GT.0.) THEN
192                 dicefracdT(i)= exposant_glace * (icefrac(i)**(exposant_glace-1.)) &
193                           / (t_glace_min - t_glace_max)
194            ENDIF
195
196            IF ((icefrac(i).EQ.0).OR.(icefrac(i).EQ.1)) THEN
197                 dicefracdT(i)=0.
198            ENDIF
199
200    ENDIF
201   
202    ENDDO
203
204    RETURN
205
206END SUBROUTINE ICEFRAC_LSCP
207!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
208
209
210
211SUBROUTINE CALC_QSAT_ECMWF(temp,qtot,pressure,tref,phase,flagth,qs,dqs)
212!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
213    ! Calculate qsat following ECMWF method
214!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
215
216    IMPLICIT NONE
217
218    include "YOMCST.h"
219    include "YOETHF.h"
220    include "FCTTRE.h"
221
222    REAL, INTENT(IN) :: temp     ! temperature in K
223    REAL, INTENT(IN) :: qtot     ! total specific water in kg/kg
224    REAL, INTENT(IN) :: pressure ! pressure in Pa
225    REAL, INTENT(IN) :: tref     ! reference temperature in K
226    LOGICAL, INTENT(IN) :: flagth     ! flag for qsat calculation for thermals
227
228    INTEGER, INTENT(IN) :: phase
229    ! phase: 0=depend on temperature sign (temp>tref -> liquid, temp<tref, solid)
230    !        1=liquid
231    !        2=solid
232
233    REAL, INTENT(OUT) :: qs      ! saturation specific humidity [kg/kg]
234    REAL, INTENT(OUT) :: dqs     ! derivation of saturation specific humidity wrt T
235
236    REAL delta, cor, cvm5
237   
238    IF (phase .EQ. 1) THEN
239        delta=0.
240    ELSEIF (phase .EQ. 2) THEN
241        delta=1.
242    ELSE
243        delta=MAX(0.,SIGN(1.,tref-temp))
244    ENDIF
245
246    IF (flagth) THEN
247    cvm5=R5LES*(1.-delta) + R5IES*delta
248    ELSE
249    cvm5 = R5LES*RLVTT*(1.-delta) + R5IES*RLSTT*delta
250    cvm5 = cvm5 /RCPD/(1.0+RVTMP2*(qtot))
251    ENDIF
252
253    qs= R2ES*FOEEW(temp,delta)/pressure
254    qs=MIN(0.5,qs)
255    cor=1./(1.-RETV*qs)
256    qs=qs*cor
257    dqs= FOEDE(temp,delta,cvm5,qs,cor)
258
259END SUBROUTINE CALC_QSAT_ECMWF
260!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
261
262
263!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
264SUBROUTINE CALC_GAMMASAT(temp,qtot,pressure,gammasat,dgammasatdt)
265
266!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
267    ! programme that calculates the gammasat parameter that determines the
268    ! homogeneous condensation thresholds for cold (<0oC) clouds
269    ! condensation at q>gammasat*qsat
270    ! Etienne Vignon, March 2021
271!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
272
273
274    IMPLICIT NONE
275
276    include "YOMCST.h"
277    include "YOETHF.h"
278    include "FCTTRE.h"
279    include "nuage.h"
280
281
282    REAL, INTENT(IN) :: temp     ! temperature in K
283    REAL, INTENT(IN) :: qtot     ! total specific water in kg/kg
284
285    REAL, INTENT(IN) :: pressure ! pressure in Pa
286
287    REAL, INTENT(OUT) :: gammasat           ! coefficient to multiply qsat with to calculate saturation
288    REAL, INTENT(OUT) :: dgammasatdt       ! derivative of gammasat wrt temperature
289
290    REAL qsi,qsl,fac,dqsl,dqsi,fcirrus
291    REAL, PARAMETER :: acirrus=2.349
292    REAL, PARAMETER :: bcirrus=259.0
293
294
295        CALL CALC_QSAT_ECMWF(temp,qtot,pressure,RTT,1,.false.,qsl,dqsl)
296        CALL CALC_QSAT_ECMWF(temp,qtot,pressure,RTT,2,.false.,qsi,dqsi)
297
298        IF (temp .GE. RTT) THEN
299            ! warm clouds: condensation at saturation wrt liquid
300            gammasat=1.
301            dgammasatdt=0.
302
303        ELSEIF ((temp .LT. RTT) .AND. (temp .GT. t_glace_min)) THEN
304           
305            IF (iflag_gammasat .GE. 2) THEN         
306               gammasat=qsl/qsi
307               dgammasatdt=(dqsl*qsi-dqsi*qsl)/qsi/qsi
308            ELSE
309               gammasat=1.
310               dgammasatdt=0.
311            ENDIF
312
313        ELSE
314
315            IF (iflag_gammasat .GE.1) THEN
316            ! homogeneous freezing of aerosols, according to
317            ! Koop, 2000 and Karcher 2008, QJRMS
318            ! 'Cirrus regime'
319               fcirrus=acirrus-temp/bcirrus
320               IF (fcirrus .LT. qsl/qsi) THEN
321                  gammasat=qsl/qsi
322                  dgammasatdt=(dqsl*qsi-dqsi*qsl)/qsi/qsi
323               ELSE
324                  gammasat=fcirrus
325                  dgammasatdt=-1.0/bcirrus
326               ENDIF
327           
328            ELSE
329
330               gammasat=1.
331               dgammasatdt=0.
332
333            ENDIF
334
335        ENDIF
336   
337END SUBROUTINE CALC_GAMMASAT
338
339
340!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
341
342END MODULE LSCP_TOOLS_MOD
343
344
Note: See TracBrowser for help on using the repository browser.