source: LMDZ6/branches/LMDZ_ECRad/libf/phylmd/lmdz_lscp_tools.F90 @ 5134

Last change on this file since 5134 was 4727, checked in by idelkadi, 14 months ago

Merged trunk changes -r4488:4726 LMDZ_ECRad branch
File size: 16.6 KB
Line 
1MODULE lmdz_lscp_tools
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 lmdz_lscp_ini, only: iflag_vice, ffallv_con, ffallv_lsc
19    use lmdz_lscp_ini, 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, temp_cltop, 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 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
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    REAL, INTENT(IN), DIMENSION(klon)   :: temp_cltop        ! temperature of cloud top
125    INTEGER, INTENT(IN)                 :: iflag_ice_thermo
126    REAL, INTENT(OUT), DIMENSION(klon)  :: icefrac
127    REAL, INTENT(OUT), DIMENSION(klon)  :: dicefracdT
128
129
130    INTEGER i
131    REAL    liqfrac_tmp, dicefrac_tmp
132    REAL    Dv, denomdep,beta,qsi,dqsidt
133    LOGICAL ice_thermo
134
135    CHARACTER (len = 20) :: modname = 'lscp_tools'
136    CHARACTER (len = 80) :: abort_message
137
138    IF ((iflag_t_glace.LT.2) .OR. (iflag_t_glace.GT.6)) THEN
139       abort_message = 'lscp cannot be used if iflag_t_glace<2 or >6'
140       CALL abort_physic(modname,abort_message,1)
141    ENDIF
142
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
147
148
149    DO i=1,klon
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)
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
165        ENDIF
166
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
179
180        ! for iflag_t_glace .GE. 4, the liquid fraction depends upon temperature at cloud top
181        ! and then decreases with decreasing height
182
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
196        IF (iflag_t_glace .EQ. 5) THEN
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
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
234
235
236     ENDDO ! klon
237 
238     RETURN
239 
240END SUBROUTINE ICEFRAC_LSCP
241!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
242
243
244
245SUBROUTINE CALC_QSAT_ECMWF(klon,temp,qtot,pressure,tref,phase,flagth,qs,dqs)
246!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
247    ! Calculate qsat following ECMWF method
248!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
249
250
251    IMPLICIT NONE
252
253    include "YOMCST.h"
254    include "YOETHF.h"
255    include "FCTTRE.h"
256
257    INTEGER, INTENT(IN) :: klon  ! number of horizontal grid points
258    REAL, INTENT(IN), DIMENSION(klon) :: temp     ! temperature in K
259    REAL, INTENT(IN), DIMENSION(klon) :: qtot     ! total specific water in kg/kg
260    REAL, INTENT(IN), DIMENSION(klon) :: pressure ! pressure in Pa
261    REAL, INTENT(IN)                  :: tref     ! reference temperature in K
262    LOGICAL, INTENT(IN) :: flagth     ! flag for qsat calculation for thermals
263    INTEGER, INTENT(IN) :: phase
264    ! phase: 0=depend on temperature sign (temp>tref -> liquid, temp<tref, solid)
265    !        1=liquid
266    !        2=solid
267
268    REAL, INTENT(OUT), DIMENSION(klon) :: qs      ! saturation specific humidity [kg/kg]
269    REAL, INTENT(OUT), DIMENSION(klon) :: dqs     ! derivation of saturation specific humidity wrt T
270
271    REAL delta, cor, cvm5
272    INTEGER i
273
274    DO i=1,klon
275
276    IF (phase .EQ. 1) THEN
277        delta=0.
278    ELSEIF (phase .EQ. 2) THEN
279        delta=1.
280    ELSE
281        delta=MAX(0.,SIGN(1.,tref-temp(i)))
282    ENDIF
283
284    IF (flagth) THEN
285    cvm5=R5LES*(1.-delta) + R5IES*delta
286    ELSE
287    cvm5 = R5LES*RLVTT*(1.-delta) + R5IES*RLSTT*delta
288    cvm5 = cvm5 /RCPD/(1.0+RVTMP2*(qtot(i)))
289    ENDIF
290
291    qs(i)= R2ES*FOEEW(temp(i),delta)/pressure(i)
292    qs(i)=MIN(0.5,qs(i))
293    cor=1./(1.-RETV*qs(i))
294    qs(i)=qs(i)*cor
295    dqs(i)= FOEDE(temp(i),delta,cvm5,qs(i),cor)
296
297    END DO
298
299END SUBROUTINE CALC_QSAT_ECMWF
300!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
301
302
303!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
304SUBROUTINE CALC_GAMMASAT(klon,temp,qtot,pressure,gammasat,dgammasatdt)
305
306!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
307    ! programme that calculates the gammasat parameter that determines the
308    ! homogeneous condensation thresholds for cold (<0oC) clouds
309    ! condensation at q>gammasat*qsat
310    ! Etienne Vignon, March 2021
311!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
312
313    use lmdz_lscp_ini, only: iflag_gammasat, t_glace_min, RTT
314
315    IMPLICIT NONE
316
317
318    INTEGER, INTENT(IN) :: klon                       ! number of horizontal grid points
319    REAL, INTENT(IN), DIMENSION(klon) :: temp         ! temperature in K
320    REAL, INTENT(IN), DIMENSION(klon) :: qtot         ! total specific water in kg/kg
321
322    REAL, INTENT(IN), DIMENSION(klon) :: pressure     ! pressure in Pa
323
324    REAL, INTENT(OUT), DIMENSION(klon) :: gammasat    ! coefficient to multiply qsat with to calculate saturation
325    REAL, INTENT(OUT), DIMENSION(klon) :: dgammasatdt ! derivative of gammasat wrt temperature
326
327    REAL, DIMENSION(klon) ::  qsi,qsl,dqsl,dqsi
328    REAL  fcirrus, fac
329    REAL, PARAMETER :: acirrus=2.349
330    REAL, PARAMETER :: bcirrus=259.0
331
332    INTEGER i
333   
334        CALL CALC_QSAT_ECMWF(klon,temp,qtot,pressure,RTT,1,.false.,qsl,dqsl)
335        CALL CALC_QSAT_ECMWF(klon,temp,qtot,pressure,RTT,2,.false.,qsi,dqsi)
336
337    DO i=1,klon
338
339        IF (temp(i) .GE. RTT) THEN
340            ! warm clouds: condensation at saturation wrt liquid
341            gammasat(i)=1.
342            dgammasatdt(i)=0.
343
344        ELSEIF ((temp(i) .LT. RTT) .AND. (temp(i) .GT. t_glace_min)) THEN
345           
346            IF (iflag_gammasat .GE. 2) THEN         
347               gammasat(i)=qsl(i)/qsi(i)
348               dgammasatdt(i)=(dqsl(i)*qsi(i)-dqsi(i)*qsl(i))/qsi(i)/qsi(i)
349            ELSE
350               gammasat(i)=1.
351               dgammasatdt(i)=0.
352            ENDIF
353
354        ELSE
355
356            IF (iflag_gammasat .GE.1) THEN
357            ! homogeneous freezing of aerosols, according to
358            ! Koop, 2000 and Karcher 2008, QJRMS
359            ! 'Cirrus regime'
360               fcirrus=acirrus-temp(i)/bcirrus
361               IF (fcirrus .LT. qsl(i)/qsi(i)) THEN
362                  gammasat(i)=qsl(i)/qsi(i)
363                  dgammasatdt(i)=(dqsl(i)*qsi(i)-dqsi(i)*qsl(i))/qsi(i)/qsi(i)
364               ELSE
365                  gammasat(i)=fcirrus
366                  dgammasatdt(i)=-1.0/bcirrus
367               ENDIF
368           
369            ELSE
370
371               gammasat(i)=1.
372               dgammasatdt(i)=0.
373
374            ENDIF
375
376        ENDIF
377   
378    END DO
379
380
381END SUBROUTINE CALC_GAMMASAT
382!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
383
384
385!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
386SUBROUTINE DISTANCE_TO_CLOUD_TOP(klon,klev,k,temp,pplay,paprs,rneb,distcltop1D,temp_cltop)
387!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
388   
389   USE lmdz_lscp_ini, ONLY : rd,rg,tresh_cl
390
391   IMPLICIT NONE
392   
393   INTEGER, INTENT(IN) :: klon,klev                !number of horizontal and vertical grid points
394   INTEGER, INTENT(IN) :: k                        ! vertical index
395   REAL, INTENT(IN), DIMENSION(klon,klev) :: temp  ! temperature in K
396   REAL, INTENT(IN), DIMENSION(klon,klev) :: pplay ! pressure middle layer in Pa
397   REAL, INTENT(IN), DIMENSION(klon,klev+1) :: paprs ! pressure interfaces in Pa
398   REAL, INTENT(IN), DIMENSION(klon,klev) :: rneb  ! cloud fraction
399
400   REAL, INTENT(OUT), DIMENSION(klon) :: distcltop1D  ! distance from cloud top
401   REAL, INTENT(OUT), DIMENSION(klon) :: temp_cltop     ! temperature of cloud top
402   
403   REAL dzlay(klon,klev)
404   REAL zlay(klon,klev)
405   REAL dzinterf
406   INTEGER i,k_top, kvert
407   LOGICAL bool_cl
408
409
410   DO i=1,klon
411         ! Initialization height middle of first layer
412          dzlay(i,1) = Rd * temp(i,1) / rg * log(paprs(i,1)/paprs(i,2))
413          zlay(i,1) = dzlay(i,1)/2
414
415          DO kvert=2,klev
416                 IF (kvert.EQ.klev) THEN
417                       dzlay(i,kvert) = 2*(rd * temp(i,kvert) / rg * log(paprs(i,kvert)/pplay(i,kvert)))
418                 ELSE
419                       dzlay(i,kvert) = rd * temp(i,kvert) / rg * log(paprs(i,kvert)/paprs(i,kvert+1))
420                 ENDIF
421                       dzinterf       = rd * temp(i,kvert) / rg * log(pplay(i,kvert-1)/pplay(i,kvert))
422                       zlay(i,kvert)  = zlay(i,kvert-1) + dzinterf
423           ENDDO
424   ENDDO
425   
426   DO i=1,klon
427          k_top = k
428          IF (rneb(i,k) .LE. tresh_cl) THEN
429                 bool_cl = .FALSE.
430          ELSE
431                 bool_cl = .TRUE.
432          ENDIF
433
434          DO WHILE ((bool_cl) .AND. (k_top .LE. klev))
435          ! find cloud top
436                IF (rneb(i,k_top) .GT. tresh_cl) THEN
437                      k_top = k_top + 1
438                ELSE
439                      bool_cl = .FALSE.
440                      k_top   = k_top - 1
441                ENDIF
442          ENDDO
443          k_top=min(k_top,klev)
444
445          !dist to top is dist between current layer and layer of cloud top (from middle to middle) + dist middle to
446          !interf for layer of cloud top
447          distcltop1D(i) = zlay(i,k_top) - zlay(i,k) + dzlay(i,k_top)/2
448          temp_cltop(i)  = temp(i,k_top)
449   ENDDO ! klon
450
451END SUBROUTINE DISTANCE_TO_CLOUD_TOP
452!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
453
454END MODULE lmdz_lscp_tools
455
456
Note: See TracBrowser for help on using the repository browser.