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

Last change on this file since 4285 was 4072, checked in by evignon, 3 years ago

Mise à jour de la routine de nuages lscp
les principaux changements consistent en:

  • des corrections de bug (déclaration et division

par zero) pour que la routine compile avec les modifications d'OB

  • des restructurations pour que les fonctions de lscp_tools_mod

travaillent sur des tableaux et non des scalaires

  • une séparation des coefficients d'évaporation de la neige et de la pluie
File size: 11.2 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
[4072]49        iwcg=MAX(iwc(i)*1000.,1E-3) ! iwc in g/m3. We set a minimum value to prevent from division by 0
[3999]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
[4072]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)
[3999]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        velo(i)=fallv_tun*min(vmsnow,vmice)/100. ! to m/s
96        dvel=0.2
[4072]97        cvel=velo(i)/((iwcg/1000.*rho(i))**dvel)
[3999]98
99    ELSE
100        ! By default, fallspeed velocity of ice crystals according to Heymsfield & Donner 1990
[4072]101        velo(i) = fallv_tun*3.29/2.0*((iwcg/1000.)**0.16)
[3999]102        dvel=0.16
103        cvel=fallv_tun*3.29/2.0*(rho(i)**0.16)
104    ENDIF
105    ENDDO
106
107END SUBROUTINE FALLICE_VELOCITY
108!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
109
110!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
111SUBROUTINE ICEFRAC_LSCP(klon, temp, sig, icefrac, dicefracdT)
112!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
113 
114  ! Compute the ice fraction 1-xliq (see e.g.
115  ! Doutriaux-Boucher & Quaas 2004, section 2.2.)
116  ! as a function of temperature
117  ! see also Fig 3 of Madeleine et al. 2020, JAMES
118!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
119
120
121    USE print_control_mod, ONLY: lunout, prt_level
122
[4059]123    IMPLICIT NONE
[3999]124
125
126    INCLUDE "YOMCST.h"
127    INCLUDE "nuage.h"
128    INCLUDE "clesphys.h"
129
130
131  ! nuage.h contains:
132  ! t_glace_min: if T < Tmin, the cloud is only made of water ice
133  ! t_glace_max: if T > Tmax, the cloud is only made of liquid water
134  ! exposant_glace: controls the sharpness of the transition
135
136    INTEGER, INTENT(IN)                 :: klon       ! number of horizontal grid points
137    REAL, INTENT(IN), DIMENSION(klon)   :: temp       ! temperature
138    REAL, INTENT(IN), DIMENSION(klon)   :: sig
139    REAL, INTENT(OUT), DIMENSION(klon)  :: icefrac
140    REAL, INTENT(OUT), DIMENSION(klon)  :: dicefracdT
141
142
143    INTEGER i
144    REAL    sig0,www,tmin_tmp,liqfrac_tmp
145    REAL    Dv, denomdep,beta,qsi,dqsidt
146    INTEGER exposant_glace_old
147    REAL t_glace_min_old
148    LOGICAL ice_thermo
149
150    sig0=0.8
151    t_glace_min_old = RTT - 15.0
152    ice_thermo = (iflag_ice_thermo .EQ. 1).OR.(iflag_ice_thermo .GE. 3)
153    IF (ice_thermo) THEN
154      exposant_glace_old = 2
155    ELSE
156      exposant_glace_old = 6
157    ENDIF
158
159
160! calculation of icefrac and dicefrac/dT
161
162    DO i=1,klon
163
164    IF (iflag_t_glace.EQ.1) THEN
165            ! Transition to ice close to surface for T<Tmax
166            ! w=1 at the surface and 0 for sig < sig0
167            www=(max(sig(i)-sig0,0.))/(1.-sig0)
168    ELSEIF (iflag_t_glace.GE.2) THEN
169            ! No convertion to ice close to surface
170            www = 0.
171    ENDIF
172     
173    tmin_tmp=www*t_glace_max+(1.-www)*t_glace_min
174    liqfrac_tmp=  (temp(i)-tmin_tmp) / (t_glace_max-tmin_tmp)
175    liqfrac_tmp = MIN(MAX(liqfrac_tmp,0.0),1.0)
176       
177    IF (iflag_t_glace.GE.3) THEN
178            icefrac(i) = 1.0-liqfrac_tmp**exposant_glace
179            IF ((icefrac(i) .GT.0.) .AND. (liqfrac_tmp .GT. 0)) THEN
180                 dicefracdT(i)= exposant_glace * ((liqfrac_tmp)**(exposant_glace-1.)) &
181                           / (t_glace_min - t_glace_max)
182            ELSE
183
184                 dicefracdT(i)=0.
185            ENDIF
186
187    ELSE
188            icefrac(i) = (1.0-liqfrac_tmp)**exposant_glace
189            IF (icefrac(i) .GT.0.) THEN
190                 dicefracdT(i)= exposant_glace * (icefrac(i)**(exposant_glace-1.)) &
191                           / (t_glace_min - t_glace_max)
192            ENDIF
193
194            IF ((icefrac(i).EQ.0).OR.(icefrac(i).EQ.1)) THEN
195                 dicefracdT(i)=0.
196            ENDIF
197
198    ENDIF
199   
200    ENDDO
201
[4072]202
[3999]203    RETURN
204
205END SUBROUTINE ICEFRAC_LSCP
206!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
207
208
209
[4072]210SUBROUTINE CALC_QSAT_ECMWF(klon,temp,qtot,pressure,tref,phase,flagth,qs,dqs)
[3999]211!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
212    ! Calculate qsat following ECMWF method
213!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
214
[4072]215
[4059]216    IMPLICIT NONE
[3999]217
218    include "YOMCST.h"
219    include "YOETHF.h"
220    include "FCTTRE.h"
221
[4072]222    INTEGER, INTENT(IN) :: klon  ! number of horizontal grid points
223    REAL, INTENT(IN), DIMENSION(klon) :: temp     ! temperature in K
224    REAL, INTENT(IN), DIMENSION(klon) :: qtot     ! total specific water in kg/kg
225    REAL, INTENT(IN), DIMENSION(klon) :: pressure ! pressure in Pa
226    REAL, INTENT(IN)                  :: tref     ! reference temperature in K
[3999]227    LOGICAL, INTENT(IN) :: flagth     ! flag for qsat calculation for thermals
228    INTEGER, INTENT(IN) :: phase
229    ! phase: 0=depend on temperature sign (temp>tref -> liquid, temp<tref, solid)
230    !        1=liquid
231    !        2=solid
232
[4072]233    REAL, INTENT(OUT), DIMENSION(klon) :: qs      ! saturation specific humidity [kg/kg]
234    REAL, INTENT(OUT), DIMENSION(klon) :: dqs     ! derivation of saturation specific humidity wrt T
[3999]235
236    REAL delta, cor, cvm5
[4072]237    INTEGER i
238
239    DO i=1,klon
240
[3999]241    IF (phase .EQ. 1) THEN
242        delta=0.
243    ELSEIF (phase .EQ. 2) THEN
244        delta=1.
245    ELSE
[4072]246        delta=MAX(0.,SIGN(1.,tref-temp(i)))
[3999]247    ENDIF
248
249    IF (flagth) THEN
250    cvm5=R5LES*(1.-delta) + R5IES*delta
251    ELSE
252    cvm5 = R5LES*RLVTT*(1.-delta) + R5IES*RLSTT*delta
[4072]253    cvm5 = cvm5 /RCPD/(1.0+RVTMP2*(qtot(i)))
[3999]254    ENDIF
255
[4072]256    qs(i)= R2ES*FOEEW(temp(i),delta)/pressure(i)
257    qs(i)=MIN(0.5,qs(i))
258    cor=1./(1.-RETV*qs(i))
259    qs(i)=qs(i)*cor
260    dqs(i)= FOEDE(temp(i),delta,cvm5,qs(i),cor)
[3999]261
[4072]262    END DO
263
[3999]264END SUBROUTINE CALC_QSAT_ECMWF
265!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
266
267
268!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
[4072]269SUBROUTINE CALC_GAMMASAT(klon,temp,qtot,pressure,gammasat,dgammasatdt)
[3999]270
271!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
272    ! programme that calculates the gammasat parameter that determines the
273    ! homogeneous condensation thresholds for cold (<0oC) clouds
274    ! condensation at q>gammasat*qsat
275    ! Etienne Vignon, March 2021
276!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
277
278
[4059]279    IMPLICIT NONE
[3999]280
281    include "YOMCST.h"
282    include "YOETHF.h"
283    include "FCTTRE.h"
284    include "nuage.h"
285
[4072]286    INTEGER, INTENT(IN) :: klon                       ! number of horizontal grid points
287    REAL, INTENT(IN), DIMENSION(klon) :: temp         ! temperature in K
288    REAL, INTENT(IN), DIMENSION(klon) :: qtot         ! total specific water in kg/kg
[3999]289
[4072]290    REAL, INTENT(IN), DIMENSION(klon) :: pressure     ! pressure in Pa
[3999]291
[4072]292    REAL, INTENT(OUT), DIMENSION(klon) :: gammasat    ! coefficient to multiply qsat with to calculate saturation
293    REAL, INTENT(OUT), DIMENSION(klon) :: dgammasatdt ! derivative of gammasat wrt temperature
[3999]294
[4072]295    REAL, DIMENSION(klon) ::  qsi,qsl,dqsl,dqsi
296    REAL  fcirrus, fac
[3999]297    REAL, PARAMETER :: acirrus=2.349
298    REAL, PARAMETER :: bcirrus=259.0
299
[4072]300    INTEGER i
301   
302        CALL CALC_QSAT_ECMWF(klon,temp,qtot,pressure,RTT,1,.false.,qsl,dqsl)
303        CALL CALC_QSAT_ECMWF(klon,temp,qtot,pressure,RTT,2,.false.,qsi,dqsi)
[3999]304
[4072]305    DO i=1,klon
[3999]306
[4072]307        IF (temp(i) .GE. RTT) THEN
[3999]308            ! warm clouds: condensation at saturation wrt liquid
[4072]309            gammasat(i)=1.
310            dgammasatdt(i)=0.
[3999]311
[4072]312        ELSEIF ((temp(i) .LT. RTT) .AND. (temp(i) .GT. t_glace_min)) THEN
[3999]313           
314            IF (iflag_gammasat .GE. 2) THEN         
[4072]315               gammasat(i)=qsl(i)/qsi(i)
316               dgammasatdt(i)=(dqsl(i)*qsi(i)-dqsi(i)*qsl(i))/qsi(i)/qsi(i)
[3999]317            ELSE
[4072]318               gammasat(i)=1.
319               dgammasatdt(i)=0.
[3999]320            ENDIF
321
322        ELSE
323
324            IF (iflag_gammasat .GE.1) THEN
325            ! homogeneous freezing of aerosols, according to
326            ! Koop, 2000 and Karcher 2008, QJRMS
327            ! 'Cirrus regime'
[4072]328               fcirrus=acirrus-temp(i)/bcirrus
329               IF (fcirrus .LT. qsl(i)/qsi(i)) THEN
330                  gammasat(i)=qsl(i)/qsi(i)
331                  dgammasatdt(i)=(dqsl(i)*qsi(i)-dqsi(i)*qsl(i))/qsi(i)/qsi(i)
[3999]332               ELSE
[4072]333                  gammasat(i)=fcirrus
334                  dgammasatdt(i)=-1.0/bcirrus
[3999]335               ENDIF
336           
337            ELSE
338
[4072]339               gammasat(i)=1.
340               dgammasatdt(i)=0.
[3999]341
342            ENDIF
343
344        ENDIF
345   
[4072]346    END DO
347
348
[3999]349END SUBROUTINE CALC_GAMMASAT
350
351
352!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
353
354END MODULE LSCP_TOOLS_MOD
355
356
Note: See TracBrowser for help on using the repository browser.