source: LMDZ6/branches/Ocean_skin/libf/phylmd/lscp_tools_mod.F90 @ 4051

Last change on this file since 4051 was 4013, checked in by lguez, 3 years ago

Sync latest trunk changes to Ocean_skin
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
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.