1 |
|
|
MODULE LSCP_TOOLS_MOD |
2 |
|
|
|
3 |
|
|
IMPLICIT NONE |
4 |
|
|
|
5 |
|
|
CONTAINS |
6 |
|
|
|
7 |
|
|
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
8 |
|
56160 |
SUBROUTINE 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 |
✓✓ |
55879200 |
DO i=1,klon |
40 |
|
|
|
41 |
✓✓ |
55823040 |
IF (ptconv(i)) THEN |
42 |
|
2488970 |
fallv_tun=ffallv_con |
43 |
|
|
ELSE |
44 |
|
53334070 |
fallv_tun=ffallv_lsc |
45 |
|
|
ENDIF |
46 |
|
|
|
47 |
|
55823040 |
tempc=temp(i)-273.15 ! celcius temp |
48 |
|
55823040 |
iwcg=MAX(iwc(i)*1000.,1E-3) ! iwc in g/m3. We set a minimum value to prevent from division by 0 |
49 |
|
55823040 |
phpa=pres(i)/100. ! pressure in hPa |
50 |
|
|
|
51 |
✗✓ |
55879200 |
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 |
✗✓ |
55823040 |
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 |
|
55823040 |
velo(i) = fallv_tun*cice_velo*((iwcg/1000.)**dice_velo) |
97 |
|
|
ENDIF |
98 |
|
|
ENDDO |
99 |
|
|
|
100 |
|
56160 |
END SUBROUTINE FALLICE_VELOCITY |
101 |
|
|
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
102 |
|
|
|
103 |
|
|
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
104 |
|
70200 |
SUBROUTINE 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 |
✗✓ |
70200 |
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 |
✗✓ |
70200 |
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 |
✓✓ |
69849000 |
DO i=1,klon |
149 |
|
|
|
150 |
|
|
! old function with sole dependence upon temperature |
151 |
✗✓ |
69778800 |
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 |
✓✗ |
69778800 |
IF (iflag_t_glace .EQ. 3) THEN |
168 |
|
69778800 |
liqfrac_tmp = (temp(i)-t_glace_min) / (t_glace_max-t_glace_min) |
169 |
|
69778800 |
liqfrac_tmp = MIN(MAX(liqfrac_tmp,0.0),1.0) |
170 |
|
69778800 |
icefrac(i) = 1.0-liqfrac_tmp**exposant_glace |
171 |
✓✓✓✓
|
69778800 |
IF ((icefrac(i) .GT.0.) .AND. (liqfrac_tmp .GT. 0.)) THEN |
172 |
|
|
dicefracdT(i)= exposant_glace * ((liqfrac_tmp)**(exposant_glace-1.)) & |
173 |
|
16075885 |
/ (t_glace_min - t_glace_max) |
174 |
|
|
ELSE |
175 |
|
53702915 |
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 |
✗✓ |
69778800 |
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 |
✗✓ |
69778800 |
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 |
✗✓ |
69849000 |
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 |
|
70200 |
RETURN |
231 |
|
|
|
232 |
|
|
END SUBROUTINE ICEFRAC_LSCP |
233 |
|
|
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
234 |
|
|
|
235 |
|
|
|
236 |
|
|
|
237 |
|
235872 |
SUBROUTINE 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 |
✓✓ |
234692640 |
DO i=1,klon |
267 |
|
|
|
268 |
✓✓ |
234456768 |
IF (phase .EQ. 1) THEN |
269 |
|
|
delta=0. |
270 |
✓✓ |
156304512 |
ELSEIF (phase .EQ. 2) THEN |
271 |
|
|
delta=1. |
272 |
|
|
ELSE |
273 |
|
78152256 |
delta=MAX(0.,SIGN(1.,tref-temp(i))) |
274 |
|
|
ENDIF |
275 |
|
|
|
276 |
✗✓ |
234456768 |
IF (flagth) THEN |
277 |
|
|
cvm5=R5LES*(1.-delta) + R5IES*delta |
278 |
|
|
ELSE |
279 |
|
234456768 |
cvm5 = R5LES*RLVTT*(1.-delta) + R5IES*RLSTT*delta |
280 |
|
234456768 |
cvm5 = cvm5 /RCPD/(1.0+RVTMP2*(qtot(i))) |
281 |
|
|
ENDIF |
282 |
|
|
|
283 |
|
234456768 |
qs(i)= R2ES*FOEEW(temp(i),delta)/pressure(i) |
284 |
|
234456768 |
qs(i)=MIN(0.5,qs(i)) |
285 |
|
234456768 |
cor=1./(1.-RETV*qs(i)) |
286 |
|
234456768 |
qs(i)=qs(i)*cor |
287 |
|
235872 |
dqs(i)= FOEDE(temp(i),delta,cvm5,qs(i),cor) |
288 |
|
|
|
289 |
|
|
END DO |
290 |
|
|
|
291 |
|
235872 |
END SUBROUTINE CALC_QSAT_ECMWF |
292 |
|
|
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
293 |
|
|
|
294 |
|
|
|
295 |
|
|
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
296 |
|
56160 |
SUBROUTINE 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 |
|
234456768 |
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 |
|
112320 |
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 |
|
56160 |
CALL CALC_QSAT_ECMWF(klon,temp,qtot,pressure,RTT,1,.false.,qsl,dqsl) |
327 |
|
56160 |
CALL CALC_QSAT_ECMWF(klon,temp,qtot,pressure,RTT,2,.false.,qsi,dqsi) |
328 |
|
|
|
329 |
✓✓ |
55879200 |
DO i=1,klon |
330 |
|
|
|
331 |
✓✓ |
55879200 |
IF (temp(i) .GE. RTT) THEN |
332 |
|
|
! warm clouds: condensation at saturation wrt liquid |
333 |
|
8162357 |
gammasat(i)=1. |
334 |
|
8162357 |
dgammasatdt(i)=0. |
335 |
|
|
|
336 |
✓✗✓✓
|
47660683 |
ELSEIF ((temp(i) .LT. RTT) .AND. (temp(i) .GT. t_glace_min)) THEN |
337 |
|
|
|
338 |
✗✓ |
12833884 |
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 |
|
12833884 |
gammasat(i)=1. |
343 |
|
12833884 |
dgammasatdt(i)=0. |
344 |
|
|
ENDIF |
345 |
|
|
|
346 |
|
|
ELSE |
347 |
|
|
|
348 |
✗✓ |
34826799 |
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 |
|
34826799 |
gammasat(i)=1. |
364 |
|
34826799 |
dgammasatdt(i)=0. |
365 |
|
|
|
366 |
|
|
ENDIF |
367 |
|
|
|
368 |
|
|
ENDIF |
369 |
|
|
|
370 |
|
|
END DO |
371 |
|
|
|
372 |
|
|
|
373 |
|
56160 |
END SUBROUTINE CALC_GAMMASAT |
374 |
|
|
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
375 |
|
|
|
376 |
|
|
|
377 |
|
|
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
378 |
|
|
SUBROUTINE 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 |
|
|
|
441 |
|
|
END SUBROUTINE DISTANCE_TO_CLOUD_TOP |
442 |
|
|
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
443 |
|
|
|
444 |
|
|
END MODULE LSCP_TOOLS_MOD |
445 |
|
|
|
446 |
|
|
|