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

Last change on this file since 4040 was 3999, checked in by evignon, 3 years ago

commission de la nouvelle routine de condensation
grande echelle simplifiee (lscp, version epuree de fistilp)
et du schema de nuages de phase mixte (en developpement)

La routine lscp n'est active que sous flag
ok_new_lscp=y
File size: 10.5 KB
RevLine 
[3999]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
205    RETURN
206
207END SUBROUTINE ICEFRAC_LSCP
208!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
209
210
211
212SUBROUTINE CALC_QSAT_ECMWF(temp,qtot,pressure,tref,phase,flagth,qs,dqs)
213!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
214    ! Calculate qsat following ECMWF method
215!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
216
217
218    IMPLICIT none
219
220    include "YOMCST.h"
221    include "YOETHF.h"
222    include "FCTTRE.h"
223
224    REAL, INTENT(IN) :: temp     ! temperature in K
225    REAL, INTENT(IN) :: qtot     ! total specific water in kg/kg
226    REAL, INTENT(IN) :: pressure ! pressure in Pa
227    REAL, INTENT(IN) :: tref     ! reference temperature in K
228    LOGICAL, INTENT(IN) :: flagth     ! flag for qsat calculation for thermals
229
230    INTEGER, INTENT(IN) :: phase
231    ! phase: 0=depend on temperature sign (temp>tref -> liquid, temp<tref, solid)
232    !        1=liquid
233    !        2=solid
234
235    REAL, INTENT(OUT) :: qs      ! saturation specific humidity [kg/kg]
236    REAL, INTENT(OUT) :: dqs     ! derivation of saturation specific humidity wrt T
237
238
239    REAL delta, cor, cvm5
240
241   
242    IF (phase .EQ. 1) THEN
243        delta=0.
244    ELSEIF (phase .EQ. 2) THEN
245        delta=1.
246    ELSE
247        delta=MAX(0.,SIGN(1.,tref-temp))
248    ENDIF
249
250    IF (flagth) THEN
251    cvm5=R5LES*(1.-delta) + R5IES*delta
252    ELSE
253    cvm5 = R5LES*RLVTT*(1.-delta) + R5IES*RLSTT*delta
254    cvm5 = cvm5 /RCPD/(1.0+RVTMP2*(qtot))
255    ENDIF
256
257    qs= R2ES*FOEEW(temp,delta)/pressure
258    qs=MIN(0.5,qs)
259    cor=1./(1.-RETV*qs)
260    qs=qs*cor
261    dqs= FOEDE(temp,delta,cvm5,qs,cor)
262
263
264
265END SUBROUTINE CALC_QSAT_ECMWF
266!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
267
268
269!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
270SUBROUTINE CALC_GAMMASAT(temp,qtot,pressure,gammasat,dgammasatdt)
271
272!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
273    ! programme that calculates the gammasat parameter that determines the
274    ! homogeneous condensation thresholds for cold (<0oC) clouds
275    ! condensation at q>gammasat*qsat
276    ! Etienne Vignon, March 2021
277!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
278
279
280    IMPLICIT none
281
282    include "YOMCST.h"
283    include "YOETHF.h"
284    include "FCTTRE.h"
285    include "nuage.h"
286
287
288    REAL, INTENT(IN) :: temp     ! temperature in K
289    REAL, INTENT(IN) :: qtot     ! total specific water in kg/kg
290
291    REAL, INTENT(IN) :: pressure ! pressure in Pa
292
293    REAL, INTENT(OUT) :: gammasat           ! coefficient to multiply qsat with to calculate saturation
294    REAL, INTENT(OUT) :: dgammasatdt       ! derivative of gammasat wrt temperature
295
296    REAL qsi,qsl,fac,dqsl,dqsi,fcirrus
297    REAL, PARAMETER :: acirrus=2.349
298    REAL, PARAMETER :: bcirrus=259.0
299
300
301        CALL CALC_QSAT_ECMWF(temp,qtot,pressure,RTT,1,.false.,qsl,dqsl)
302        CALL CALC_QSAT_ECMWF(temp,qtot,pressure,RTT,2,.false.,qsi,dqsi)
303
304        IF (temp .GE. RTT) THEN
305            ! warm clouds: condensation at saturation wrt liquid
306            gammasat=1.
307            dgammasatdt=0.
308
309        ELSEIF ((temp .LT. RTT) .AND. (temp .GT. t_glace_min)) THEN
310           
311            IF (iflag_gammasat .GE. 2) THEN         
312               gammasat=qsl/qsi
313               dgammasatdt=(dqsl*qsi-dqsi*qsl)/qsi/qsi
314            ELSE
315               gammasat=1.
316               dgammasatdt=0.
317            ENDIF
318
319        ELSE
320
321            IF (iflag_gammasat .GE.1) THEN
322            ! homogeneous freezing of aerosols, according to
323            ! Koop, 2000 and Karcher 2008, QJRMS
324            ! 'Cirrus regime'
325               fcirrus=acirrus-temp/bcirrus
326               IF (fcirrus .LT. qsl/qsi) THEN
327                  gammasat=qsl/qsi
328                  dgammasatdt=(dqsl*qsi-dqsi*qsl)/qsi/qsi
329               ELSE
330                  gammasat=fcirrus
331                  dgammasatdt=-1.0/bcirrus
332               ENDIF
333           
334            ELSE
335
336               gammasat=1.
337               dgammasatdt=0.
338
339            ENDIF
340
341        ENDIF
342   
343
344
345
346END SUBROUTINE CALC_GAMMASAT
347
348
349!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
350
351END MODULE LSCP_TOOLS_MOD
352
353
Note: See TracBrowser for help on using the repository browser.