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

Last change on this file since 4545 was 4535, checked in by evignon, 18 months ago

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