source: LMDZ6/branches/LMDZ_cdrag_LSCE/libf/phylmd/lscp_tools_mod.F90 @ 5444

Last change on this file since 5444 was 4669, checked in by Laurent Fairhead, 16 months ago

Merged with trunk revision 4586 corresponding to june 2023 testing
File size: 16.0 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    use lscp_ini_mod, only: iflag_vice, ffallv_con, ffallv_lsc
19    use lscp_ini_mod, 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
33    INTEGER i
34    REAL logvm,iwcg,tempc,phpa,fallv_tun
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
48        iwcg=MAX(iwc(i)*1000.,1E-3) ! iwc in g/m3. We set a minimum value to prevent from division by 0
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
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)
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
96        velo(i) = fallv_tun*cice_velo*((iwcg/1000.)**dice_velo)
97    ENDIF
98    ENDDO
99
100END SUBROUTINE FALLICE_VELOCITY
101!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
102
103!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
104SUBROUTINE ICEFRAC_LSCP(klon, temp, iflag_ice_thermo, distcltop, icefrac, dicefracdT)
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
115    USE lscp_ini_mod, ONLY: t_glace_min, t_glace_max, exposant_glace, iflag_t_glace
116    USE lscp_ini_mod, ONLY : RTT, dist_liq
117
118    IMPLICIT NONE
119
120
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    INTEGER, INTENT(IN)                 :: iflag_ice_thermo
125    REAL, INTENT(OUT), DIMENSION(klon)  :: icefrac
126    REAL, INTENT(OUT), DIMENSION(klon)  :: dicefracdT
127
128
129    INTEGER i
130    REAL    liqfrac_tmp, dicefrac_tmp
131    REAL    Dv, denomdep,beta,qsi,dqsidt
132    LOGICAL ice_thermo
133
134    CHARACTER (len = 20) :: modname = 'lscp_tools'
135    CHARACTER (len = 80) :: abort_message
136
137    IF ((iflag_t_glace.LT.2) .OR. (iflag_t_glace.GT.6)) THEN
138       abort_message = 'lscp cannot be used if iflag_t_glace<2 or >6'
139       CALL abort_physic(modname,abort_message,1)
140    ENDIF
141
142    IF (.NOT.((iflag_ice_thermo .EQ. 1).OR.(iflag_ice_thermo .GE. 3))) THEN
143       abort_message = 'lscp cannot be used without ice thermodynamics'
144       CALL abort_physic(modname,abort_message,1)
145    ENDIF
146
147
148    DO i=1,klon
149 
150        ! old function with sole dependence upon temperature
151        IF (iflag_t_glace .EQ. 2) THEN
152            liqfrac_tmp = (temp(i)-t_glace_min) / (t_glace_max-t_glace_min)
153            liqfrac_tmp = MIN(MAX(liqfrac_tmp,0.0),1.0)
154            icefrac(i) = (1.0-liqfrac_tmp)**exposant_glace
155            IF (icefrac(i) .GT.0.) THEN
156                 dicefracdT(i)= exposant_glace * (icefrac(i)**(exposant_glace-1.)) &
157                           / (t_glace_min - t_glace_max)
158            ENDIF
159
160            IF ((icefrac(i).EQ.0).OR.(icefrac(i).EQ.1)) THEN
161                 dicefracdT(i)=0.
162            ENDIF
163
164        ENDIF
165
166        ! function of temperature used in CMIP6 physics
167        IF (iflag_t_glace .EQ. 3) THEN
168            liqfrac_tmp = (temp(i)-t_glace_min) / (t_glace_max-t_glace_min)
169            liqfrac_tmp = MIN(MAX(liqfrac_tmp,0.0),1.0)
170            icefrac(i) = 1.0-liqfrac_tmp**exposant_glace
171            IF ((icefrac(i) .GT.0.) .AND. (liqfrac_tmp .GT. 0.)) THEN
172                dicefracdT(i)= exposant_glace * ((liqfrac_tmp)**(exposant_glace-1.)) &
173                          / (t_glace_min - t_glace_max)
174            ELSE
175                dicefracdT(i)=0.
176            ENDIF
177        ENDIF
178
179        ! for iflag_t_glace .GE. 4, the liquid fraction depends upon temperature at cloud top
180        ! and then decreases with decreasing height
181
182        !with linear function of temperature at cloud top
183        IF (iflag_t_glace .EQ. 4) THEN 
184                liqfrac_tmp = (temp(i)-t_glace_min) / (t_glace_max-t_glace_min)
185                liqfrac_tmp = MIN(MAX(liqfrac_tmp,0.0),1.0)
186                icefrac(i)    =  MAX(MIN(1.,1.0 - liqfrac_tmp*exp(-distcltop(i)/dist_liq)),0.)
187                dicefrac_tmp = - temp(i)/(t_glace_max-t_glace_min)
188                dicefracdT(i) = dicefrac_tmp*exp(-distcltop(i)/dist_liq)
189                IF ((liqfrac_tmp .LE.0) .OR. (liqfrac_tmp .GE. 1)) THEN
190                        dicefracdT(i) = 0.
191                ENDIF
192        ENDIF
193
194        ! with CMIP6 function of temperature at cloud top
195        IF (iflag_t_glace .EQ. 5) THEN
196                liqfrac_tmp = (temp(i)-t_glace_min) / (t_glace_max-t_glace_min)
197                liqfrac_tmp =  MIN(MAX(liqfrac_tmp,0.0),1.0)
198                liqfrac_tmp = liqfrac_tmp**exposant_glace
199                icefrac(i)  =  MAX(MIN(1.,1.0 - liqfrac_tmp*exp(-distcltop(i)/dist_liq)),0.)
200                IF ((liqfrac_tmp .LE.0) .OR. (liqfrac_tmp .GE. 1)) THEN
201                        dicefracdT(i) = 0.
202                ELSE
203                        dicefracdT(i) = exposant_glace*((liqfrac_tmp)**(exposant_glace-1.))/(t_glace_min- t_glace_max) &
204                                        *exp(-distcltop(i)/dist_liq)
205                ENDIF
206        ENDIF
207
208        ! with modified function of temperature at cloud top
209        ! to get largere values around 260 K, works well with t_glace_min = 241K
210        IF (iflag_t_glace .EQ. 6) THEN
211                IF (temp(i) .GT. t_glace_max) THEN
212                        liqfrac_tmp = 1.
213                ELSE
214                        liqfrac_tmp = -((temp(i)-t_glace_max) / (t_glace_max-t_glace_min))**2+1.
215                ENDIF
216                liqfrac_tmp  = MIN(MAX(liqfrac_tmp,0.0),1.0)
217                icefrac(i)   = MAX(MIN(1.,1.0 - liqfrac_tmp*exp(-distcltop(i)/dist_liq)),0.)       
218                IF ((liqfrac_tmp .LE.0) .OR. (liqfrac_tmp .GE. 1)) THEN
219                        dicefracdT(i) = 0.
220                ELSE
221                        dicefracdT(i) = 2*((temp(i)-t_glace_max) / (t_glace_max-t_glace_min))/(t_glace_max-t_glace_min) &
222                                  *exp(-distcltop(i)/dist_liq)
223                ENDIF
224        ENDIF
225
226
227
228     ENDDO ! klon
229 
230     RETURN
231 
232END SUBROUTINE ICEFRAC_LSCP
233!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
234
235
236
237SUBROUTINE CALC_QSAT_ECMWF(klon,temp,qtot,pressure,tref,phase,flagth,qs,dqs)
238!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
239    ! Calculate qsat following ECMWF method
240!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
241
242
243    IMPLICIT NONE
244
245    include "YOMCST.h"
246    include "YOETHF.h"
247    include "FCTTRE.h"
248
249    INTEGER, INTENT(IN) :: klon  ! number of horizontal grid points
250    REAL, INTENT(IN), DIMENSION(klon) :: temp     ! temperature in K
251    REAL, INTENT(IN), DIMENSION(klon) :: qtot     ! total specific water in kg/kg
252    REAL, INTENT(IN), DIMENSION(klon) :: pressure ! pressure in Pa
253    REAL, INTENT(IN)                  :: tref     ! reference temperature in K
254    LOGICAL, INTENT(IN) :: flagth     ! flag for qsat calculation for thermals
255    INTEGER, INTENT(IN) :: phase
256    ! phase: 0=depend on temperature sign (temp>tref -> liquid, temp<tref, solid)
257    !        1=liquid
258    !        2=solid
259
260    REAL, INTENT(OUT), DIMENSION(klon) :: qs      ! saturation specific humidity [kg/kg]
261    REAL, INTENT(OUT), DIMENSION(klon) :: dqs     ! derivation of saturation specific humidity wrt T
262
263    REAL delta, cor, cvm5
264    INTEGER i
265
266    DO i=1,klon
267
268    IF (phase .EQ. 1) THEN
269        delta=0.
270    ELSEIF (phase .EQ. 2) THEN
271        delta=1.
272    ELSE
273        delta=MAX(0.,SIGN(1.,tref-temp(i)))
274    ENDIF
275
276    IF (flagth) THEN
277    cvm5=R5LES*(1.-delta) + R5IES*delta
278    ELSE
279    cvm5 = R5LES*RLVTT*(1.-delta) + R5IES*RLSTT*delta
280    cvm5 = cvm5 /RCPD/(1.0+RVTMP2*(qtot(i)))
281    ENDIF
282
283    qs(i)= R2ES*FOEEW(temp(i),delta)/pressure(i)
284    qs(i)=MIN(0.5,qs(i))
285    cor=1./(1.-RETV*qs(i))
286    qs(i)=qs(i)*cor
287    dqs(i)= FOEDE(temp(i),delta,cvm5,qs(i),cor)
288
289    END DO
290
291END SUBROUTINE CALC_QSAT_ECMWF
292!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
293
294
295!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
296SUBROUTINE CALC_GAMMASAT(klon,temp,qtot,pressure,gammasat,dgammasatdt)
297
298!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
299    ! programme that calculates the gammasat parameter that determines the
300    ! homogeneous condensation thresholds for cold (<0oC) clouds
301    ! condensation at q>gammasat*qsat
302    ! Etienne Vignon, March 2021
303!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
304
305    use lscp_ini_mod, only: iflag_gammasat, t_glace_min, RTT
306
307    IMPLICIT NONE
308
309
310    INTEGER, INTENT(IN) :: klon                       ! number of horizontal grid points
311    REAL, INTENT(IN), DIMENSION(klon) :: temp         ! temperature in K
312    REAL, INTENT(IN), DIMENSION(klon) :: qtot         ! total specific water in kg/kg
313
314    REAL, INTENT(IN), DIMENSION(klon) :: pressure     ! pressure in Pa
315
316    REAL, INTENT(OUT), DIMENSION(klon) :: gammasat    ! coefficient to multiply qsat with to calculate saturation
317    REAL, INTENT(OUT), DIMENSION(klon) :: dgammasatdt ! derivative of gammasat wrt temperature
318
319    REAL, DIMENSION(klon) ::  qsi,qsl,dqsl,dqsi
320    REAL  fcirrus, fac
321    REAL, PARAMETER :: acirrus=2.349
322    REAL, PARAMETER :: bcirrus=259.0
323
324    INTEGER i
325   
326        CALL CALC_QSAT_ECMWF(klon,temp,qtot,pressure,RTT,1,.false.,qsl,dqsl)
327        CALL CALC_QSAT_ECMWF(klon,temp,qtot,pressure,RTT,2,.false.,qsi,dqsi)
328
329    DO i=1,klon
330
331        IF (temp(i) .GE. RTT) THEN
332            ! warm clouds: condensation at saturation wrt liquid
333            gammasat(i)=1.
334            dgammasatdt(i)=0.
335
336        ELSEIF ((temp(i) .LT. RTT) .AND. (temp(i) .GT. t_glace_min)) THEN
337           
338            IF (iflag_gammasat .GE. 2) THEN         
339               gammasat(i)=qsl(i)/qsi(i)
340               dgammasatdt(i)=(dqsl(i)*qsi(i)-dqsi(i)*qsl(i))/qsi(i)/qsi(i)
341            ELSE
342               gammasat(i)=1.
343               dgammasatdt(i)=0.
344            ENDIF
345
346        ELSE
347
348            IF (iflag_gammasat .GE.1) THEN
349            ! homogeneous freezing of aerosols, according to
350            ! Koop, 2000 and Karcher 2008, QJRMS
351            ! 'Cirrus regime'
352               fcirrus=acirrus-temp(i)/bcirrus
353               IF (fcirrus .LT. qsl(i)/qsi(i)) THEN
354                  gammasat(i)=qsl(i)/qsi(i)
355                  dgammasatdt(i)=(dqsl(i)*qsi(i)-dqsi(i)*qsl(i))/qsi(i)/qsi(i)
356               ELSE
357                  gammasat(i)=fcirrus
358                  dgammasatdt(i)=-1.0/bcirrus
359               ENDIF
360           
361            ELSE
362
363               gammasat(i)=1.
364               dgammasatdt(i)=0.
365
366            ENDIF
367
368        ENDIF
369   
370    END DO
371
372
373END SUBROUTINE CALC_GAMMASAT
374!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
375
376
377!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
378SUBROUTINE DISTANCE_TO_CLOUD_TOP(klon,klev,k,temp,pplay,paprs,rneb,distcltop1D)
379!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
380   
381   USE lscp_ini_mod, ONLY : rd,rg,tresh_cl
382
383   IMPLICIT NONE
384   
385   INTEGER, INTENT(IN) :: klon,klev                !number of horizontal and vertical grid points
386   INTEGER, INTENT(IN) :: k                        ! vertical index
387   REAL, INTENT(IN), DIMENSION(klon,klev) :: temp  ! temperature in K
388   REAL, INTENT(IN), DIMENSION(klon,klev) :: pplay ! pressure middle layer in Pa
389   REAL, INTENT(IN), DIMENSION(klon,klev+1) :: paprs ! pressure interfaces in Pa
390   REAL, INTENT(IN), DIMENSION(klon,klev) :: rneb  ! cloud fraction
391
392   REAL, INTENT(OUT), DIMENSION(klon) :: distcltop1D  ! distance from cloud top
393
394   REAL dzlay(klon,klev)
395   REAL zlay(klon,klev)
396   REAL dzinterf
397   INTEGER i,k_top, kvert
398   LOGICAL bool_cl
399
400
401   DO i=1,klon
402         ! Initialization height middle of first layer
403          dzlay(i,1) = Rd * temp(i,1) / rg * log(paprs(i,1)/paprs(i,2))
404          zlay(i,1) = dzlay(i,1)/2
405
406          DO kvert=2,klev
407                 IF (kvert.EQ.klev) THEN
408                       dzlay(i,kvert) = 2*(rd * temp(i,kvert) / rg * log(paprs(i,kvert)/pplay(i,kvert)))
409                 ELSE
410                       dzlay(i,kvert) = rd * temp(i,kvert) / rg * log(paprs(i,kvert)/paprs(i,kvert+1))
411                 ENDIF
412                       dzinterf       = rd * temp(i,kvert) / rg * log(pplay(i,kvert-1)/pplay(i,kvert))
413                       zlay(i,kvert)  = zlay(i,kvert-1) + dzinterf
414           ENDDO
415   ENDDO
416   
417   k_top = k
418   DO i=1,klon
419          IF (rneb(i,k) .LE. tresh_cl) THEN
420                 bool_cl = .FALSE.
421          ELSE
422                 bool_cl = .TRUE.
423          ENDIF
424
425          DO WHILE ((bool_cl) .AND. (k_top .LE. klev))
426          ! find cloud top
427                IF (rneb(i,k_top) .GT. tresh_cl) THEN
428                      k_top = k_top + 1
429                ELSE
430                      bool_cl = .FALSE.
431                      k_top   = k_top - 1
432                ENDIF
433          ENDDO
434          k_top=min(k_top,klev)
435
436          !dist to top is dist between current layer and layer of cloud top (from middle to middle) + dist middle to
437          !interf for layer of cloud top
438          distcltop1D(i) = zlay(i,k_top) - zlay(i,k) + dzlay(i,k_top)/2
439   ENDDO ! klon
440
441END SUBROUTINE DISTANCE_TO_CLOUD_TOP
442!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
443
444END MODULE LSCP_TOOLS_MOD
445
446
Note: See TracBrowser for help on using the repository browser.