1 |
|
|
MODULE ice_sursat_mod |
2 |
|
|
|
3 |
|
|
IMPLICIT NONE |
4 |
|
|
|
5 |
|
|
!--flight inventories |
6 |
|
|
! |
7 |
|
|
REAL, SAVE, ALLOCATABLE :: flight_m(:,:) !--flown distance m s-1 per cell |
8 |
|
|
!$OMP THREADPRIVATE(flight_m) |
9 |
|
|
REAL, SAVE, ALLOCATABLE :: flight_h2o(:,:) !--emitted kg H2O s-1 per cell |
10 |
|
|
!$OMP THREADPRIVATE(flight_h2o) |
11 |
|
|
! |
12 |
|
|
!--Fixed Parameters |
13 |
|
|
! |
14 |
|
|
!--safety parameters for ERF function |
15 |
|
|
REAL, PARAMETER :: erf_lim = 5., eps = 1.e-10 |
16 |
|
|
! |
17 |
|
|
!--Tuning parameters (and their default values) |
18 |
|
|
! |
19 |
|
|
!--chi gère la répartition statistique de la longueur des frontières |
20 |
|
|
! entre les zones nuages et ISSR/ciel clair sous-saturé. Gamme de valeur : |
21 |
|
|
! chi > 1, je n'ai pas regardé de limite max (pour chi = 1, la longueur de |
22 |
|
|
! la frontière entre ne nuage et l'ISSR est proportionnelle à la |
23 |
|
|
! répartition ISSR/ciel clair sous-sat dans la maille, i.e. il n'y a pas |
24 |
|
|
! de favorisation de la localisation de l'ISSR près de nuage. Pour chi = inf, |
25 |
|
|
! le nuage n'est en contact qu'avec de l'ISSR, quelle que soit la taille |
26 |
|
|
! de l'ISSR dans la maille.) |
27 |
|
|
! |
28 |
|
|
!--l_turb est la longueur de mélange pour la turbulence. |
29 |
|
|
! dans les tests, ça n'a jamais été modifié pour l'instant. |
30 |
|
|
! |
31 |
|
|
!--tun_N est le paramètre qui contrôle l'importance relative de N_2 par rapport à N_1. |
32 |
|
|
! La valeur est comprise entre 1 et 2 (tun_N = 1 => N_1 = N_2) |
33 |
|
|
! |
34 |
|
|
!--tun_ratqs : paramètre qui modifie ratqs en fonction de la valeur de |
35 |
|
|
! alpha_cld selon la formule ratqs_new = ratqs_old / ( 1 + tun_ratqs * |
36 |
|
|
! alpha_cld ). Dans le rapport il est appelé beta. Il varie entre 0 et 5 |
37 |
|
|
! (tun_ratqs = 0 => pas de modification de ratqs). |
38 |
|
|
! |
39 |
|
|
!--gamma0 and Tgamma: define RHcrit limit above which heterogeneous freezing occurs as a function of T |
40 |
|
|
!--Karcher and Lohmann (2002) |
41 |
|
|
!--gamma = 2.583 - t / 207.83 |
42 |
|
|
!--Ren and MacKenzie (2005) reused by Kärcher |
43 |
|
|
!--gamma = 2.349 - t / 259.0 |
44 |
|
|
! |
45 |
|
|
!--N_cld: number of clouds in cell (needs to be parametrized at some point) |
46 |
|
|
! |
47 |
|
|
!--contrail cross section: typical value found in Freudenthaler et al, GRL, 22, 3501-3504, 1995 |
48 |
|
|
!--in m2, 1000x200 = 200 000 m2 after 15 min |
49 |
|
|
! |
50 |
|
|
REAL, SAVE :: chi=1.1, l_turb=50.0, tun_N=1.3, tun_ratqs=3.0 |
51 |
|
|
REAL, SAVE :: gamma0=2.349, Tgamma=259.0, N_cld=100, contrail_cross_section=200000.0 |
52 |
|
|
!$OMP THREADPRIVATE(chi,l_turb,tun_N,tun_ratqs,gamma0,Tgamma,N_cld,contrail_cross_section) |
53 |
|
|
|
54 |
|
|
CONTAINS |
55 |
|
|
|
56 |
|
|
!******************************************************************* |
57 |
|
|
! |
58 |
|
|
SUBROUTINE ice_sursat_init() |
59 |
|
|
|
60 |
|
|
USE print_control_mod, ONLY: lunout |
61 |
|
|
USE ioipsl_getin_p_mod, ONLY : getin_p |
62 |
|
|
|
63 |
|
|
IMPLICIT NONE |
64 |
|
|
|
65 |
|
|
CALL getin_p('flag_chi',chi) |
66 |
|
|
CALL getin_p('flag_l_turb',l_turb) |
67 |
|
|
CALL getin_p('flag_tun_N',tun_N) |
68 |
|
|
CALL getin_p('flag_tun_ratqs',tun_ratqs) |
69 |
|
|
CALL getin_p('gamma0',gamma0) |
70 |
|
|
CALL getin_p('Tgamma',Tgamma) |
71 |
|
|
CALL getin_p('N_cld',N_cld) |
72 |
|
|
CALL getin_p('contrail_cross_section',contrail_cross_section) |
73 |
|
|
|
74 |
|
|
WRITE(lunout,*) 'Parameters for ice_sursat param' |
75 |
|
|
WRITE(lunout,*) 'flag_chi = ', chi |
76 |
|
|
WRITE(lunout,*) 'flag_l_turb = ', l_turb |
77 |
|
|
WRITE(lunout,*) 'flag_tun_N = ', tun_N |
78 |
|
|
WRITE(lunout,*) 'flag_tun_ratqs = ', tun_ratqs |
79 |
|
|
WRITE(lunout,*) 'gamma0 = ', gamma0 |
80 |
|
|
WRITE(lunout,*) 'Tgamma = ', Tgamma |
81 |
|
|
WRITE(lunout,*) 'N_cld = ', N_cld |
82 |
|
|
WRITE(lunout,*) 'contrail_cross_section = ', contrail_cross_section |
83 |
|
|
|
84 |
|
|
END SUBROUTINE ice_sursat_init |
85 |
|
|
|
86 |
|
|
!******************************************************************* |
87 |
|
|
! |
88 |
|
|
SUBROUTINE airplane(debut,pphis,pplay,paprs,t_seri) |
89 |
|
|
|
90 |
|
|
USE dimphy |
91 |
|
|
USE mod_grid_phy_lmdz, ONLY: klon_glo |
92 |
|
|
USE geometry_mod, ONLY: cell_area |
93 |
|
|
USE phys_cal_mod, ONLY : mth_cur |
94 |
|
|
USE mod_phys_lmdz_mpi_data, ONLY: is_mpi_root |
95 |
|
|
USE mod_phys_lmdz_omp_data, ONLY: is_omp_root |
96 |
|
|
USE mod_phys_lmdz_para, ONLY: scatter, bcast |
97 |
|
|
USE print_control_mod, ONLY: lunout |
98 |
|
|
|
99 |
|
|
IMPLICIT NONE |
100 |
|
|
|
101 |
|
|
INCLUDE "YOMCST.h" |
102 |
|
|
INCLUDE 'netcdf.inc' |
103 |
|
|
|
104 |
|
|
!-------------------------------------------------------- |
105 |
|
|
!--input variables |
106 |
|
|
!-------------------------------------------------------- |
107 |
|
|
LOGICAL, INTENT(IN) :: debut |
108 |
|
|
REAL, INTENT(IN) :: pphis(klon), pplay(klon,klev), paprs(klon,klev+1), t_seri(klon,klev) |
109 |
|
|
|
110 |
|
|
!-------------------------------------------------------- |
111 |
|
|
! ... Local variables |
112 |
|
|
!-------------------------------------------------------- |
113 |
|
|
|
114 |
|
|
CHARACTER (LEN=20) :: modname='airplane_mod' |
115 |
|
|
INTEGER :: i, k, kori, iret, varid, error, ncida, klona |
116 |
|
|
INTEGER,SAVE :: nleva, ntimea |
117 |
|
|
!$OMP THREADPRIVATE(nleva,ntimea) |
118 |
|
|
REAL, ALLOCATABLE :: pkm_airpl_glo(:,:,:) !--km/s |
119 |
|
|
REAL, ALLOCATABLE :: ph2o_airpl_glo(:,:,:) !--molec H2O/cm3/s |
120 |
|
|
REAL, ALLOCATABLE, SAVE :: zmida(:), zinta(:) |
121 |
|
|
REAL, ALLOCATABLE, SAVE :: pkm_airpl(:,:,:) |
122 |
|
|
REAL, ALLOCATABLE, SAVE :: ph2o_airpl(:,:,:) |
123 |
|
|
!$OMP THREADPRIVATE(pkm_airpl,ph2o_airpl,zmida,zinta) |
124 |
|
|
REAL :: zalt(klon,klev+1) |
125 |
|
|
REAL :: zrho, zdz(klon,klev), zfrac |
126 |
|
|
|
127 |
|
|
! |
128 |
|
|
IF (debut) THEN |
129 |
|
|
!-------------------------------------------------------------------------------- |
130 |
|
|
! ... Open the file and read airplane emissions |
131 |
|
|
!-------------------------------------------------------------------------------- |
132 |
|
|
! |
133 |
|
|
IF (is_mpi_root .AND. is_omp_root) THEN |
134 |
|
|
! |
135 |
|
|
iret = nf_open('aircraft_phy.nc', 0, ncida) |
136 |
|
|
IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to open aircraft_phy.nc file',1) |
137 |
|
|
! ... Get lengths |
138 |
|
|
iret = nf_inq_dimid(ncida, 'time', varid) |
139 |
|
|
IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to get time dimid in aircraft_phy.nc file',1) |
140 |
|
|
iret = nf_inq_dimlen(ncida, varid, ntimea) |
141 |
|
|
IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to get time dimlen aircraft_phy.nc file',1) |
142 |
|
|
iret = nf_inq_dimid(ncida, 'vector', varid) |
143 |
|
|
IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to get vector dimid aircraft_phy.nc file',1) |
144 |
|
|
iret = nf_inq_dimlen(ncida, varid, klona) |
145 |
|
|
IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to get vector dimlen aircraft_phy.nc file',1) |
146 |
|
|
iret = nf_inq_dimid(ncida, 'lev', varid) |
147 |
|
|
IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to get lev dimid aircraft_phy.nc file',1) |
148 |
|
|
iret = nf_inq_dimlen(ncida, varid, nleva) |
149 |
|
|
IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to get lev dimlen aircraft_phy.nc file',1) |
150 |
|
|
! |
151 |
|
|
IF ( klona /= klon_glo ) THEN |
152 |
|
|
WRITE(lunout,*) 'klona & klon_glo =', klona, klon_glo |
153 |
|
|
CALL abort_physic(modname,'problem klon in aircraft_phy.nc file',1) |
154 |
|
|
ENDIF |
155 |
|
|
! |
156 |
|
|
IF ( ntimea /= 12 ) THEN |
157 |
|
|
WRITE(lunout,*) 'ntimea=', ntimea |
158 |
|
|
CALL abort_physic(modname,'problem ntime<>12 in aircraft_phy.nc file',1) |
159 |
|
|
ENDIF |
160 |
|
|
! |
161 |
|
|
ALLOCATE(zmida(nleva), STAT=error) |
162 |
|
|
IF (error /= 0) CALL abort_physic(modname,'problem to allocate zmida',1) |
163 |
|
|
ALLOCATE(pkm_airpl_glo(klona,nleva,ntimea), STAT=error) |
164 |
|
|
IF (error /= 0) CALL abort_physic(modname,'problem to allocate pkm_airpl_glo',1) |
165 |
|
|
ALLOCATE(ph2o_airpl_glo(klona,nleva,ntimea), STAT=error) |
166 |
|
|
IF (error /= 0) CALL abort_physic(modname,'problem to allocate ph2o_airpl_glo',1) |
167 |
|
|
! |
168 |
|
|
iret = nf_inq_varid(ncida, 'lev', varid) |
169 |
|
|
IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to get lev dimid aircraft_phy.nc file',1) |
170 |
|
|
iret = nf_get_var_double(ncida, varid, zmida) |
171 |
|
|
IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to read zmida file',1) |
172 |
|
|
! |
173 |
|
|
iret = nf_inq_varid(ncida, 'emi_co2_aircraft', varid) !--CO2 as a proxy for m flown - |
174 |
|
|
IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to get emi_distance dimid aircraft_phy.nc file',1) |
175 |
|
|
iret = nf_get_var_double(ncida, varid, pkm_airpl_glo) |
176 |
|
|
IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to read pkm_airpl file',1) |
177 |
|
|
! |
178 |
|
|
iret = nf_inq_varid(ncida, 'emi_h2o_aircraft', varid) |
179 |
|
|
IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to get emi_h2o_aircraft dimid aircraft_phy.nc file',1) |
180 |
|
|
iret = nf_get_var_double(ncida, varid, ph2o_airpl_glo) |
181 |
|
|
IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to read ph2o_airpl file',1) |
182 |
|
|
! |
183 |
|
|
ENDIF !--is_mpi_root and is_omp_root |
184 |
|
|
! |
185 |
|
|
CALL bcast(nleva) |
186 |
|
|
CALL bcast(ntimea) |
187 |
|
|
! |
188 |
|
|
IF (.NOT.ALLOCATED(zmida)) ALLOCATE(zmida(nleva), STAT=error) |
189 |
|
|
IF (.NOT.ALLOCATED(zinta)) ALLOCATE(zinta(nleva+1), STAT=error) |
190 |
|
|
! |
191 |
|
|
ALLOCATE(pkm_airpl(klon,nleva,ntimea)) |
192 |
|
|
ALLOCATE(ph2o_airpl(klon,nleva,ntimea)) |
193 |
|
|
! |
194 |
|
|
ALLOCATE(flight_m(klon,klev)) |
195 |
|
|
ALLOCATE(flight_h2o(klon,klev)) |
196 |
|
|
! |
197 |
|
|
CALL bcast(zmida) |
198 |
|
|
zinta(1)=0.0 !--surface |
199 |
|
|
DO k=2, nleva |
200 |
|
|
zinta(k) = (zmida(k-1)+zmida(k))/2.0*1000.0 !--conversion from km to m |
201 |
|
|
ENDDO |
202 |
|
|
zinta(nleva+1)=zinta(nleva)+(zmida(nleva)-zmida(nleva-1))*1000.0 !--extrapolation for last interface |
203 |
|
|
!print *,'zinta=', zinta |
204 |
|
|
! |
205 |
|
|
CALL scatter(pkm_airpl_glo,pkm_airpl) |
206 |
|
|
CALL scatter(ph2o_airpl_glo,ph2o_airpl) |
207 |
|
|
! |
208 |
|
|
!$OMP MASTER |
209 |
|
|
IF (is_mpi_root .AND. is_omp_root) THEN |
210 |
|
|
DEALLOCATE(pkm_airpl_glo) |
211 |
|
|
DEALLOCATE(ph2o_airpl_glo) |
212 |
|
|
ENDIF !--is_mpi_root |
213 |
|
|
!$OMP END MASTER |
214 |
|
|
|
215 |
|
|
ENDIF !--debut |
216 |
|
|
! |
217 |
|
|
!--compute altitude of model level interfaces |
218 |
|
|
! |
219 |
|
|
DO i = 1, klon |
220 |
|
|
zalt(i,1)=pphis(i)/RG !--in m |
221 |
|
|
ENDDO |
222 |
|
|
! |
223 |
|
|
DO k=1, klev |
224 |
|
|
DO i = 1, klon |
225 |
|
|
zrho=pplay(i,k)/t_seri(i,k)/RD |
226 |
|
|
zdz(i,k)=(paprs(i,k)-paprs(i,k+1))/zrho/RG |
227 |
|
|
zalt(i,k+1)=zalt(i,k)+zdz(i,k) !--in m |
228 |
|
|
ENDDO |
229 |
|
|
ENDDO |
230 |
|
|
! |
231 |
|
|
!--vertical reprojection |
232 |
|
|
! |
233 |
|
|
flight_m(:,:)=0.0 |
234 |
|
|
flight_h2o(:,:)=0.0 |
235 |
|
|
! |
236 |
|
|
DO k=1, klev |
237 |
|
|
DO kori=1, nleva |
238 |
|
|
DO i=1, klon |
239 |
|
|
!--fraction of layer kori included in layer k |
240 |
|
|
zfrac=max(0.0,min(zalt(i,k+1),zinta(kori+1))-max(zalt(i,k),zinta(kori)))/(zinta(kori+1)-zinta(kori)) |
241 |
|
|
!--reproject |
242 |
|
|
flight_m(i,k)=flight_m(i,k) + pkm_airpl(i,kori,mth_cur)*zfrac |
243 |
|
|
!--reproject |
244 |
|
|
flight_h2o(i,k)=flight_h2o(i,k) + ph2o_airpl(i,kori,mth_cur)*zfrac |
245 |
|
|
ENDDO |
246 |
|
|
ENDDO |
247 |
|
|
ENDDO |
248 |
|
|
! |
249 |
|
|
DO k=1, klev |
250 |
|
|
DO i=1, klon |
251 |
|
|
!--molec.cm-3.s-1 / (molec/mol) * kg CO2/mol * m2 * m * cm3/m3 / (kg CO2/m) => m s-1 per cell |
252 |
|
|
flight_m(i,k)=flight_m(i,k)/RNAVO*44.e-3*cell_area(i)*zdz(i,k)*1.e6/16.37e-3 |
253 |
|
|
flight_m(i,k)=flight_m(i,k)*100.0 !--x100 to augment signal to noise |
254 |
|
|
!--molec.cm-3.s-1 / (molec/mol) * kg H2O/mol * m2 * m * cm3/m3 => kg H2O s-1 per cell |
255 |
|
|
flight_h2o(i,k)=flight_h2o(i,k)/RNAVO*18.e-3*cell_area(i)*zdz(i,k)*1.e6 |
256 |
|
|
ENDDO |
257 |
|
|
ENDDO |
258 |
|
|
! |
259 |
|
|
END SUBROUTINE airplane |
260 |
|
|
|
261 |
|
|
!******************************************************************** |
262 |
|
|
! simple routine to initialise flight_m and test a flight corridor |
263 |
|
|
!--Olivier Boucher - 2021 |
264 |
|
|
! |
265 |
|
|
SUBROUTINE flight_init() |
266 |
|
|
USE dimphy |
267 |
|
|
USE geometry_mod, ONLY: cell_area, latitude_deg, longitude_deg |
268 |
|
|
IMPLICIT NONE |
269 |
|
|
INTEGER :: i |
270 |
|
|
|
271 |
|
|
ALLOCATE(flight_m(klon,klev)) |
272 |
|
|
ALLOCATE(flight_h2o(klon,klev)) |
273 |
|
|
! |
274 |
|
|
flight_m(:,:) = 0.0 !--initialisation |
275 |
|
|
flight_h2o(:,:) = 0.0 !--initialisation |
276 |
|
|
! |
277 |
|
|
DO i=1, klon |
278 |
|
|
IF (latitude_deg(i).GE.42.0.AND.latitude_deg(i).LE.48.0) THEN |
279 |
|
|
flight_m(i,38) = 50000.0 !--5000 m of flight/second in grid cell x 10 scaling |
280 |
|
|
ENDIF |
281 |
|
|
ENDDO |
282 |
|
|
|
283 |
|
|
RETURN |
284 |
|
|
END SUBROUTINE flight_init |
285 |
|
|
|
286 |
|
|
!******************************************************************* |
287 |
|
|
!--Routine to deal with ice supersaturation |
288 |
|
|
!--Determines the respective fractions of unsaturated clear sky, ice supersaturated clear sky and cloudy sky |
289 |
|
|
!--Diagnoses regions prone for non-persistent and persistent contrail formation |
290 |
|
|
! |
291 |
|
|
!--Audran Borella - 2021 |
292 |
|
|
! |
293 |
|
|
SUBROUTINE ice_sursat(pplay, dpaprs, dtime, i, k, t, q, gamma_ss, & |
294 |
|
|
qsat, t_actuel, rneb_seri, ratqs, rneb, qincld, & |
295 |
|
|
rnebss, qss, Tcontr, qcontr, qcontr2, fcontrN, fcontrP) |
296 |
|
|
! |
297 |
|
|
USE dimphy |
298 |
|
|
USE print_control_mod, ONLY: prt_level, lunout |
299 |
|
|
USE phys_state_var_mod, ONLY: pbl_tke, t_ancien |
300 |
|
|
USE phys_local_var_mod, ONLY: N1_ss, N2_ss |
301 |
|
|
USE phys_local_var_mod, ONLY: drneb_sub, drneb_con, drneb_tur, drneb_avi |
302 |
|
|
!! USE phys_local_var_mod, ONLY: Tcontr, qcontr, fcontrN, fcontrP |
303 |
|
|
USE indice_sol_mod, ONLY: is_ave |
304 |
|
|
USE geometry_mod, ONLY: cell_area |
305 |
|
|
! |
306 |
|
|
IMPLICIT NONE |
307 |
|
|
INCLUDE "YOMCST.h" |
308 |
|
|
INCLUDE "YOETHF.h" |
309 |
|
|
INCLUDE "FCTTRE.h" |
310 |
|
|
INCLUDE "clesphys.h" |
311 |
|
|
|
312 |
|
|
! |
313 |
|
|
! Input |
314 |
|
|
! Beware: this routine works on a gridpoint! |
315 |
|
|
! |
316 |
|
|
REAL, INTENT(IN) :: pplay ! layer pressure (Pa) |
317 |
|
|
REAL, INTENT(IN) :: dpaprs ! layer delta pressure (Pa) |
318 |
|
|
REAL, INTENT(IN) :: dtime ! intervalle du temps (s) |
319 |
|
|
REAL, INTENT(IN) :: t ! température advectée (K) |
320 |
|
|
REAL, INTENT(IN) :: qsat ! vapeur de saturation |
321 |
|
|
REAL, INTENT(IN) :: t_actuel ! temperature actuelle de la maille (K) |
322 |
|
|
REAL, INTENT(IN) :: rneb_seri ! fraction nuageuse en memoire |
323 |
|
|
INTEGER, INTENT(IN) :: i, k |
324 |
|
|
! |
325 |
|
|
! Input/output |
326 |
|
|
! |
327 |
|
|
REAL, INTENT(INOUT) :: q ! vapeur de la maille (=zq) |
328 |
|
|
REAL, INTENT(INOUT) :: ratqs ! determine la largeur de distribution de vapeur |
329 |
|
|
REAL, INTENT(INOUT) :: Tcontr, qcontr, qcontr2, fcontrN, fcontrP |
330 |
|
|
! |
331 |
|
|
! Output |
332 |
|
|
! |
333 |
|
|
REAL, INTENT(OUT) :: gamma_ss ! |
334 |
|
|
REAL, INTENT(OUT) :: rneb ! cloud fraction |
335 |
|
|
REAL, INTENT(OUT) :: qincld ! in-cloud total water |
336 |
|
|
REAL, INTENT(OUT) :: rnebss ! ISSR fraction |
337 |
|
|
REAL, INTENT(OUT) :: qss ! in-ISSR total water |
338 |
|
|
! |
339 |
|
|
! Local |
340 |
|
|
! |
341 |
|
|
REAL PI |
342 |
|
|
PARAMETER (PI=4.*ATAN(1.)) |
343 |
|
|
REAL rnebclr, gamma_prec |
344 |
|
|
REAL qclr, qvc, qcld, qi |
345 |
|
|
REAL zrho, zdz, zrhodz |
346 |
|
|
REAL pdf_N, pdf_N1, pdf_N2 |
347 |
|
|
REAL pdf_a, pdf_b |
348 |
|
|
REAL pdf_e1, pdf_e2, pdf_k |
349 |
|
|
REAL drnebss, drnebclr, dqss, dqclr, sum_rneb_rnebss, dqss_avi |
350 |
|
|
REAL V_cell !--volume of the cell |
351 |
|
|
REAL M_cell !--dry mass of the cell |
352 |
|
|
REAL tke, sig, L_tur, b_tur, q_eq |
353 |
|
|
REAL V_env, V_cld, V_ss, V_clr |
354 |
|
|
REAL zcor |
355 |
|
|
! |
356 |
|
|
!--more local variables for diagnostics |
357 |
|
|
!--imported from YOMCST.h |
358 |
|
|
!--eps_w = 0.622 = ratio of molecular masses of water and dry air (kg H2O kg air -1) |
359 |
|
|
!--RCPD = 1004 J kg air−1 K−1 = the isobaric heat capacity of air |
360 |
|
|
!--values from Schumann, Meteorol Zeitschrift, 1996 |
361 |
|
|
!--EiH2O = 1.25 / 2.24 / 8.94 kg H2O / kg fuel for kerosene / methane / dihydrogen |
362 |
|
|
!--Qheat = 43. / 50. / 120. MJ / kg fuel for kerosene / methane / dihydrogen |
363 |
|
|
REAL, PARAMETER :: EiH2O=1.25 !--emission index of water vapour for kerosene (kg kg-1) |
364 |
|
|
REAL, PARAMETER :: Qheat=43.E6 !--specific combustion heat for kerosene (J kg-1) |
365 |
|
|
REAL, PARAMETER :: eta=0.3 !--average propulsion efficiency of the aircraft |
366 |
|
|
!--Gcontr is the slope of the mean phase trajectory in the turbulent exhaust field on an absolute |
367 |
|
|
!--temperature versus water vapor partial pressure diagram. G has the unit of Pa K−1. Rap et al JGR 2010. |
368 |
|
|
REAL :: Gcontr |
369 |
|
|
!--Tcontr = critical temperature for contrail formation (T_LM in Schumann 1996, Eq 31 in appendix 2) |
370 |
|
|
!--qsatliqcontr = e_L(T_LM) in Schumann 1996 but expressed in specific humidity (kg kg humid air-1) |
371 |
|
|
REAL :: qsatliqcontr |
372 |
|
|
|
373 |
|
|
! Initialisations |
374 |
|
|
zrho = pplay / t / RD !--dry density kg m-3 |
375 |
|
|
zrhodz = dpaprs / RG !--dry air mass kg m-2 |
376 |
|
|
zdz = zrhodz / zrho !--cell thickness m |
377 |
|
|
V_cell = zdz * cell_area(i) !--cell volume m3 |
378 |
|
|
M_cell = zrhodz * cell_area(i) !--cell dry air mass kg |
379 |
|
|
! |
380 |
|
|
! Recuperation de la memoire sur la couverture nuageuse |
381 |
|
|
rneb = rneb_seri |
382 |
|
|
! |
383 |
|
|
! Ajout des émissions de H2O dues à l'aviation |
384 |
|
|
! q is the specific humidity (kg/kg humid air) hence the complicated equation to update q |
385 |
|
|
! qnew = ( m_humid_air * qold + dm_H2O ) / ( m_humid_air + dm_H2O ) |
386 |
|
|
! = ( m_dry_air * qold + dm_h2O * (1-qold) ) / (m_dry_air + dm_H2O * (1-qold) ) |
387 |
|
|
! The equation is derived by writing m_humid_air = m_dry_air + m_H2O = m_dry_air / (1-q) |
388 |
|
|
! flight_h2O is in kg H2O / s / cell |
389 |
|
|
! |
390 |
|
|
IF (ok_plane_h2o) THEN |
391 |
|
|
q = ( M_cell*q + flight_h2o(i,k)*dtime*(1.-q) ) / (M_cell + flight_h2o(i,k)*dtime*(1.-q) ) |
392 |
|
|
ENDIF |
393 |
|
|
! |
394 |
|
|
!--Estimating gamma |
395 |
|
|
gamma_ss = MAX(1.0, gamma0 - t_actuel/Tgamma) |
396 |
|
|
!gamma_prec = MAX(1.0, gamma0 - t_ancien(i,k)/Tgamma) !--formulation initiale d Audran |
397 |
|
|
gamma_prec = MAX(1.0, gamma0 - t/Tgamma) !--autre formulation possible basée sur le T du pas de temps |
398 |
|
|
! |
399 |
|
|
! Initialisation de qvc : q_sat du pas de temps precedent |
400 |
|
|
!qvc = R2ES*FOEEW(t_ancien(i,k),1.)/pplay !--formulation initiale d Audran |
401 |
|
|
qvc = R2ES*FOEEW(t,1.)/pplay !--autre formulation possible basée sur le T du pas de temps |
402 |
|
|
qvc = min(0.5,qvc) |
403 |
|
|
zcor = 1./(1.-RETV*qvc) |
404 |
|
|
qvc = qvc*zcor |
405 |
|
|
! |
406 |
|
|
! Modification de ratqs selon formule proposee : ksi_new = ksi_old/(1+beta*alpha_cld) |
407 |
|
|
ratqs = ratqs / (tun_ratqs*rneb_seri + 1.) |
408 |
|
|
! |
409 |
|
|
! Calcul de N |
410 |
|
|
pdf_k = -sqrt(log(1.+ratqs**2.)) |
411 |
|
|
pdf_a = log(qvc/q)/(pdf_k*sqrt(2.)) |
412 |
|
|
pdf_b = pdf_k/(2.*sqrt(2.)) |
413 |
|
|
pdf_e1 = pdf_a+pdf_b |
414 |
|
|
IF (abs(pdf_e1).GE.erf_lim) THEN |
415 |
|
|
pdf_e1 = sign(1.,pdf_e1) |
416 |
|
|
pdf_N = max(0.,sign(rneb,pdf_e1)) |
417 |
|
|
ELSE |
418 |
|
|
pdf_e1 = erf(pdf_e1) |
419 |
|
|
pdf_e1 = 0.5*(1.+pdf_e1) |
420 |
|
|
pdf_N = max(0.,rneb/pdf_e1) |
421 |
|
|
ENDIF |
422 |
|
|
! |
423 |
|
|
! On calcule ensuite N1 et N2. Il y a deux cas : N > 1 et N <= 1 |
424 |
|
|
! Cas 1 : N > 1. N'arrive en theorie jamais, c'est une barriere |
425 |
|
|
! On perd la memoire sur la temperature (sur qvc) pour garder |
426 |
|
|
! celle sur alpha_cld |
427 |
|
|
IF (pdf_N.GT.1.) THEN |
428 |
|
|
! On inverse alpha_cld = int_qvc^infty P(q) dq |
429 |
|
|
! pour determiner qvc = f(alpha_cld) |
430 |
|
|
! On approxime en serie entiere erf-1(x) |
431 |
|
|
qvc = 2.*rneb-1. |
432 |
|
|
qvc = qvc + PI/12.*qvc**3 + 7.*PI**2/480.*qvc**5 & |
433 |
|
|
+ 127.*PI**3/40320.*qvc**7 + 4369.*PI**4/5806080.*qvc**9 & |
434 |
|
|
+ 34807.*PI**5/182476800.*qvc**11 |
435 |
|
|
qvc = sqrt(PI)/2.*qvc |
436 |
|
|
qvc = (qvc-pdf_b)*pdf_k*sqrt(2.) |
437 |
|
|
qvc = q*exp(qvc) |
438 |
|
|
|
439 |
|
|
! On met a jour rneb avec la nouvelle valeur de qvc |
440 |
|
|
! La maj est necessaire a cause de la serie entiere approximative |
441 |
|
|
pdf_a = log(qvc/q)/(pdf_k*sqrt(2.)) |
442 |
|
|
pdf_e1 = pdf_a+pdf_b |
443 |
|
|
IF (abs(pdf_e1).GE.erf_lim) THEN |
444 |
|
|
pdf_e1 = sign(1.,pdf_e1) |
445 |
|
|
ELSE |
446 |
|
|
pdf_e1 = erf(pdf_e1) |
447 |
|
|
ENDIF |
448 |
|
|
pdf_e1 = 0.5*(1.+pdf_e1) |
449 |
|
|
rneb = pdf_e1 |
450 |
|
|
|
451 |
|
|
! Si N > 1, N1 et N2 = 1 |
452 |
|
|
pdf_N1 = 1. |
453 |
|
|
pdf_N2 = 1. |
454 |
|
|
|
455 |
|
|
! Cas 2 : N <= 1 |
456 |
|
|
ELSE |
457 |
|
|
! D'abord on calcule N2 avec le tuning |
458 |
|
|
pdf_N2 = min(1.,pdf_N*tun_N) |
459 |
|
|
|
460 |
|
|
! Puis N1 pour assurer la conservation de alpha_cld |
461 |
|
|
pdf_a = log(qvc*gamma_prec/q)/(pdf_k*sqrt(2.)) |
462 |
|
|
pdf_e2 = pdf_a+pdf_b |
463 |
|
|
IF (abs(pdf_e2).GE.erf_lim) THEN |
464 |
|
|
pdf_e2 = sign(1.,pdf_e2) |
465 |
|
|
ELSE |
466 |
|
|
pdf_e2 = erf(pdf_e2) |
467 |
|
|
ENDIF |
468 |
|
|
pdf_e2 = 0.5*(1.+pdf_e2) ! integrale sous P pour q > gamma qsat |
469 |
|
|
|
470 |
|
|
IF (abs(pdf_e1-pdf_e2).LT.eps) THEN |
471 |
|
|
pdf_N1 = pdf_N2 |
472 |
|
|
ELSE |
473 |
|
|
pdf_N1 = (rneb-pdf_N2*pdf_e2)/(pdf_e1-pdf_e2) |
474 |
|
|
ENDIF |
475 |
|
|
|
476 |
|
|
! Barriere qui traite le cas gamma_prec = 1. |
477 |
|
|
IF (pdf_N1.LE.0.) THEN |
478 |
|
|
pdf_N1 = 0. |
479 |
|
|
IF (pdf_e2.GT.eps) THEN |
480 |
|
|
pdf_N2 = rneb/pdf_e2 |
481 |
|
|
ELSE |
482 |
|
|
pdf_N2 = 0. |
483 |
|
|
ENDIF |
484 |
|
|
ENDIF |
485 |
|
|
ENDIF |
486 |
|
|
|
487 |
|
|
! Physique 1 |
488 |
|
|
! Sublimation |
489 |
|
|
IF (qvc.LT.qsat) THEN |
490 |
|
|
pdf_a = log(qvc/q)/(pdf_k*sqrt(2.)) |
491 |
|
|
pdf_e1 = pdf_a+pdf_b |
492 |
|
|
IF (abs(pdf_e1).GE.erf_lim) THEN |
493 |
|
|
pdf_e1 = sign(1.,pdf_e1) |
494 |
|
|
ELSE |
495 |
|
|
pdf_e1 = erf(pdf_e1) |
496 |
|
|
ENDIF |
497 |
|
|
|
498 |
|
|
pdf_a = log(qsat/q)/(pdf_k*sqrt(2.)) |
499 |
|
|
pdf_e2 = pdf_a+pdf_b |
500 |
|
|
IF (abs(pdf_e2).GE.erf_lim) THEN |
501 |
|
|
pdf_e2 = sign(1.,pdf_e2) |
502 |
|
|
ELSE |
503 |
|
|
pdf_e2 = erf(pdf_e2) |
504 |
|
|
ENDIF |
505 |
|
|
|
506 |
|
|
pdf_e1 = 0.5*pdf_N1*(pdf_e1-pdf_e2) |
507 |
|
|
|
508 |
|
|
! Calcul et ajout de la tendance |
509 |
|
|
drneb_sub(i,k) = - pdf_e1/dtime !--unit [s-1] |
510 |
|
|
rneb = rneb + drneb_sub(i,k)*dtime |
511 |
|
|
ELSE |
512 |
|
|
drneb_sub(i,k) = 0. |
513 |
|
|
ENDIF |
514 |
|
|
|
515 |
|
|
! NOTE : verifier si ca marche bien pour gamma proche de 1. |
516 |
|
|
|
517 |
|
|
! Condensation |
518 |
|
|
IF (gamma_ss*qsat.LT.gamma_prec*qvc) THEN |
519 |
|
|
|
520 |
|
|
pdf_a = log(gamma_ss*qsat/q)/(pdf_k*sqrt(2.)) |
521 |
|
|
pdf_e1 = pdf_a+pdf_b |
522 |
|
|
IF (abs(pdf_e1).GE.erf_lim) THEN |
523 |
|
|
pdf_e1 = sign(1.,pdf_e1) |
524 |
|
|
ELSE |
525 |
|
|
pdf_e1 = erf(pdf_e1) |
526 |
|
|
ENDIF |
527 |
|
|
|
528 |
|
|
pdf_a = log(gamma_prec*qvc/q)/(pdf_k*sqrt(2.)) |
529 |
|
|
pdf_e2 = pdf_a+pdf_b |
530 |
|
|
IF (abs(pdf_e2).GE.erf_lim) THEN |
531 |
|
|
pdf_e2 = sign(1.,pdf_e2) |
532 |
|
|
ELSE |
533 |
|
|
pdf_e2 = erf(pdf_e2) |
534 |
|
|
ENDIF |
535 |
|
|
|
536 |
|
|
pdf_e1 = 0.5*(1.-pdf_N1)*(pdf_e1-pdf_e2) |
537 |
|
|
pdf_e2 = 0.5*(1.-pdf_N2)*(1.+pdf_e2) |
538 |
|
|
|
539 |
|
|
! Calcul et ajout de la tendance |
540 |
|
|
drneb_con(i,k) = (pdf_e1 + pdf_e2)/dtime !--unit [s-1] |
541 |
|
|
rneb = rneb + drneb_con(i,k)*dtime |
542 |
|
|
|
543 |
|
|
ELSE |
544 |
|
|
|
545 |
|
|
pdf_a = log(gamma_ss*qsat/q)/(pdf_k*sqrt(2.)) |
546 |
|
|
pdf_e1 = pdf_a+pdf_b |
547 |
|
|
IF (abs(pdf_e1).GE.erf_lim) THEN |
548 |
|
|
pdf_e1 = sign(1.,pdf_e1) |
549 |
|
|
ELSE |
550 |
|
|
pdf_e1 = erf(pdf_e1) |
551 |
|
|
ENDIF |
552 |
|
|
pdf_e1 = 0.5*(1.-pdf_N2)*(1.+pdf_e1) |
553 |
|
|
|
554 |
|
|
! Calcul et ajout de la tendance |
555 |
|
|
drneb_con(i,k) = pdf_e1/dtime !--unit [s-1] |
556 |
|
|
rneb = rneb + drneb_con(i,k)*dtime |
557 |
|
|
|
558 |
|
|
ENDIF |
559 |
|
|
|
560 |
|
|
! Calcul des grandeurs diagnostiques |
561 |
|
|
! Determination des grandeurs ciel clair |
562 |
|
|
pdf_a = log(qsat/q)/(pdf_k*sqrt(2.)) |
563 |
|
|
pdf_e1 = pdf_a+pdf_b |
564 |
|
|
IF (abs(pdf_e1).GE.erf_lim) THEN |
565 |
|
|
pdf_e1 = sign(1.,pdf_e1) |
566 |
|
|
ELSE |
567 |
|
|
pdf_e1 = erf(pdf_e1) |
568 |
|
|
ENDIF |
569 |
|
|
pdf_e1 = 0.5*(1.-pdf_e1) |
570 |
|
|
|
571 |
|
|
pdf_e2 = pdf_a-pdf_b |
572 |
|
|
IF (abs(pdf_e2).GE.erf_lim) THEN |
573 |
|
|
pdf_e2 = sign(1.,pdf_e2) |
574 |
|
|
ELSE |
575 |
|
|
pdf_e2 = erf(pdf_e2) |
576 |
|
|
ENDIF |
577 |
|
|
pdf_e2 = 0.5*q*(1.-pdf_e2) |
578 |
|
|
|
579 |
|
|
rnebclr = pdf_e1 |
580 |
|
|
qclr = pdf_e2 |
581 |
|
|
|
582 |
|
|
! Determination de q_cld |
583 |
|
|
! Partie 1 |
584 |
|
|
pdf_a = log(max(qsat,qvc)/q)/(pdf_k*sqrt(2.)) |
585 |
|
|
pdf_e1 = pdf_a-pdf_b |
586 |
|
|
IF (abs(pdf_e1).GE.erf_lim) THEN |
587 |
|
|
pdf_e1 = sign(1.,pdf_e1) |
588 |
|
|
ELSE |
589 |
|
|
pdf_e1 = erf(pdf_e1) |
590 |
|
|
ENDIF |
591 |
|
|
|
592 |
|
|
pdf_a = log(min(gamma_ss*qsat,gamma_prec*qvc)/q)/(pdf_k*sqrt(2.)) |
593 |
|
|
pdf_e2 = pdf_a-pdf_b |
594 |
|
|
IF (abs(pdf_e2).GE.erf_lim) THEN |
595 |
|
|
pdf_e2 = sign(1.,pdf_e2) |
596 |
|
|
ELSE |
597 |
|
|
pdf_e2 = erf(pdf_e2) |
598 |
|
|
ENDIF |
599 |
|
|
|
600 |
|
|
pdf_e1 = 0.5*q*pdf_N1*(pdf_e1-pdf_e2) |
601 |
|
|
|
602 |
|
|
qcld = pdf_e1 |
603 |
|
|
|
604 |
|
|
! Partie 2 (sous condition) |
605 |
|
|
IF (gamma_ss*qsat.GT.gamma_prec*qvc) THEN |
606 |
|
|
pdf_a = log(gamma_ss*qsat/q)/(pdf_k*sqrt(2.)) |
607 |
|
|
pdf_e1 = pdf_a-pdf_b |
608 |
|
|
IF (abs(pdf_e1).GE.erf_lim) THEN |
609 |
|
|
pdf_e1 = sign(1.,pdf_e1) |
610 |
|
|
ELSE |
611 |
|
|
pdf_e1 = erf(pdf_e1) |
612 |
|
|
ENDIF |
613 |
|
|
|
614 |
|
|
pdf_e2 = 0.5*q*pdf_N2*(pdf_e2-pdf_e1) |
615 |
|
|
|
616 |
|
|
qcld = qcld + pdf_e2 |
617 |
|
|
|
618 |
|
|
pdf_e2 = pdf_e1 |
619 |
|
|
ENDIF |
620 |
|
|
|
621 |
|
|
! Partie 3 |
622 |
|
|
pdf_e2 = 0.5*q*(1.+pdf_e2) |
623 |
|
|
|
624 |
|
|
qcld = qcld + pdf_e2 |
625 |
|
|
|
626 |
|
|
! Fin du calcul de q_cld |
627 |
|
|
|
628 |
|
|
! Determination des grandeurs ISSR via les equations de conservation |
629 |
|
|
rneb=MIN(rneb, 1. - rnebclr - eps) !--ajout OB - barrière |
630 |
|
|
rnebss = MAX(0.0, 1. - rnebclr - rneb) !--ajout OB |
631 |
|
|
qss = MAX(0.0, q - qclr - qcld) !--ajout OB |
632 |
|
|
|
633 |
|
|
! Physique 2 : Turbulence |
634 |
|
|
IF (rneb.GT.eps.AND.rneb.LT.1.-eps) THEN ! rneb != 0 and != 1 |
635 |
|
|
! |
636 |
|
|
tke = pbl_tke(i,k,is_ave) |
637 |
|
|
! A MODIFIER formule a revoir |
638 |
|
|
L_tur = min(l_turb, sqrt(tke)*dtime) |
639 |
|
|
|
640 |
|
|
! On fait pour l'instant l'hypothese a = 3b. V = 4/3 pi a b**2 = alpha_cld |
641 |
|
|
! donc b = alpha_cld/4pi **1/3. |
642 |
|
|
b_tur = (rneb*V_cell/4./PI/N_cld)**(1./3.) |
643 |
|
|
! On verifie que la longeur de melange n'est pas trop grande |
644 |
|
|
IF (L_tur.GT.b_tur) THEN |
645 |
|
|
L_tur = b_tur |
646 |
|
|
ENDIF |
647 |
|
|
|
648 |
|
|
V_env = N_cld*4.*PI*(3.*(b_tur**2.)*L_tur + L_tur**3. + 3.*b_tur*(L_tur**2.)) |
649 |
|
|
V_cld = N_cld*4.*PI*(3.*(b_tur**2.)*L_tur + L_tur**3. - 3.*b_tur*(L_tur**2.)) |
650 |
|
|
V_cld = abs(V_cld) |
651 |
|
|
|
652 |
|
|
! Repartition statistique des zones de contact nuage-ISSR et nuage-ciel clair |
653 |
|
|
sig = rnebss/(chi*rnebclr+rnebss) |
654 |
|
|
V_ss = MIN(sig*V_env,rnebss*V_cell) |
655 |
|
|
V_clr = MIN((1.-sig)*V_env,rnebclr*V_cell) |
656 |
|
|
V_cld = MIN(V_cld,rneb*V_cell) |
657 |
|
|
V_env = V_ss + V_clr |
658 |
|
|
|
659 |
|
|
! ISSR => rneb |
660 |
|
|
drnebss = -1.*V_ss/V_cell |
661 |
|
|
dqss = drnebss*qss/MAX(eps,rnebss) |
662 |
|
|
|
663 |
|
|
! Clear sky <=> rneb |
664 |
|
|
q_eq = V_env*qclr/MAX(eps,rnebclr) + V_cld*qcld/MAX(eps,rneb) |
665 |
|
|
q_eq = q_eq/(V_env + V_cld) |
666 |
|
|
|
667 |
|
|
IF (q_eq.GT.qsat) THEN |
668 |
|
|
drnebclr = - V_clr/V_cell |
669 |
|
|
dqclr = drnebclr*qclr/MAX(eps,rnebclr) |
670 |
|
|
ELSE |
671 |
|
|
drnebclr = V_cld/V_cell |
672 |
|
|
dqclr = drnebclr*qcld/MAX(eps,rneb) |
673 |
|
|
ENDIF |
674 |
|
|
|
675 |
|
|
! Maj des variables avec les tendances |
676 |
|
|
rnebclr = MAX(0.0,rnebclr + drnebclr) !--OB ajout d'un max avec eps (il faudrait modified drnebclr pour le diag) |
677 |
|
|
qclr = MAX(0.0, qclr + dqclr) !--OB ajout d'un max avec 0 |
678 |
|
|
|
679 |
|
|
rneb = rneb - drnebclr - drnebss |
680 |
|
|
qcld = qcld - dqclr - dqss |
681 |
|
|
|
682 |
|
|
rnebss = MAX(0.0,rnebss + drnebss) !--OB ajout d'un max avec eps (il faudrait modifier drnebss pour le diag) |
683 |
|
|
qss = MAX(0.0, qss + dqss) !--OB ajout d'un max avec 0 |
684 |
|
|
|
685 |
|
|
! Tendances pour le diagnostic |
686 |
|
|
drneb_tur(i,k) = (drnebclr + drnebss)/dtime !--unit [s-1] |
687 |
|
|
|
688 |
|
|
ENDIF ! rneb |
689 |
|
|
|
690 |
|
|
!--add a source of cirrus from aviation contrails |
691 |
|
|
IF (ok_plane_contrail) THEN |
692 |
|
|
drneb_avi(i,k) = rnebss*flight_m(i,k)*contrail_cross_section/V_cell !--tendency rneb due to aviation [s-1] |
693 |
|
|
drneb_avi(i,k) = MIN(drneb_avi(i,k), rnebss/dtime) !--majoration |
694 |
|
|
dqss_avi = qss*drneb_avi(i,k)/MAX(eps,rnebss) !--tendency q aviation [kg kg-1 s-1] |
695 |
|
|
rneb = rneb + drneb_avi(i,k)*dtime !--add tendency to rneb |
696 |
|
|
qcld = qcld + dqss_avi*dtime !--add tendency to qcld |
697 |
|
|
rnebss = rnebss - drneb_avi(i,k)*dtime !--add tendency to rnebss |
698 |
|
|
qss = qss - dqss_avi*dtime !--add tendency to qss |
699 |
|
|
ELSE |
700 |
|
|
drneb_avi(i,k)=0.0 |
701 |
|
|
ENDIF |
702 |
|
|
|
703 |
|
|
! Barrieres |
704 |
|
|
! ISSR trop petite |
705 |
|
|
IF (rnebss.LT.eps) THEN |
706 |
|
|
rneb = MIN(rneb + rnebss,1.0-eps) !--ajout OB barriere |
707 |
|
|
qcld = qcld + qss |
708 |
|
|
rnebss = 0. |
709 |
|
|
qss = 0. |
710 |
|
|
ENDIF |
711 |
|
|
|
712 |
|
|
! le nuage est trop petit |
713 |
|
|
IF (rneb.LT.eps) THEN |
714 |
|
|
! s'il y a une ISSR on met tout dans l'ISSR, sinon dans le |
715 |
|
|
! clear sky |
716 |
|
|
IF (rnebss.LT.eps) THEN |
717 |
|
|
rnebclr = 1. |
718 |
|
|
rnebss = 0. !--ajout OB |
719 |
|
|
qclr = q |
720 |
|
|
ELSE |
721 |
|
|
rnebss = MIN(rnebss + rneb,1.0-eps) !--ajout OB barriere |
722 |
|
|
qss = qss + qcld |
723 |
|
|
ENDIF |
724 |
|
|
rneb = 0. |
725 |
|
|
qcld = 0. |
726 |
|
|
qincld = qsat * gamma_ss |
727 |
|
|
ELSE |
728 |
|
|
qincld = qcld / rneb |
729 |
|
|
ENDIF |
730 |
|
|
|
731 |
|
|
!--OB ajout borne superieure |
732 |
|
|
sum_rneb_rnebss=rneb+rnebss |
733 |
|
|
rneb=rneb*MIN(1.-eps,sum_rneb_rnebss)/MAX(eps,sum_rneb_rnebss) |
734 |
|
|
rnebss=rnebss*MIN(1.-eps,sum_rneb_rnebss)/MAX(eps,sum_rneb_rnebss) |
735 |
|
|
|
736 |
|
|
! On ecrit dans la memoire |
737 |
|
|
N1_ss(i,k) = pdf_N1 |
738 |
|
|
N2_ss(i,k) = pdf_N2 |
739 |
|
|
|
740 |
|
|
!--Diagnostics only used from last iteration |
741 |
|
|
!--test |
742 |
|
|
!!Tcontr(i,k)=200. |
743 |
|
|
!!fcontrN(i,k)=1.0 |
744 |
|
|
!!fcontrP(i,k)=0.5 |
745 |
|
|
! |
746 |
|
|
!--slope of dilution line in exhaust |
747 |
|
|
!--kg H2O/kg fuel * J kg air-1 K-1 * Pa / (kg H2O / kg air * J kg fuel-1) |
748 |
|
|
Gcontr = EiH2O * RCPD * pplay / (eps_w*Qheat*(1.-eta)) !--Pa K-1 |
749 |
|
|
!--critical T_LM below which no liquid contrail can form in exhaust |
750 |
|
|
!Tcontr(i,k) = 226.69+9.43*log(Gcontr-0.053)+0.72*(log(Gcontr-0.053))**2 !--K |
751 |
|
|
IF (Gcontr .GT. 0.1) THEN |
752 |
|
|
! |
753 |
|
|
Tcontr = 226.69+9.43*log(Gcontr-0.053)+0.72*(log(Gcontr-0.053))**2 !--K |
754 |
|
|
!print *,'Tcontr=',iter,i,k,eps_w,pplay,Gcontr,Tcontr(i,k) |
755 |
|
|
!--Psat with index 0 in FOEEW to get saturation wrt liquid water corresponding to Tcontr |
756 |
|
|
!qsatliqcontr = RESTT*FOEEW(Tcontr(i,k),0.) !--Pa |
757 |
|
|
qsatliqcontr = RESTT*FOEEW(Tcontr,0.) !--Pa |
758 |
|
|
!--Critical water vapour above which there is contrail formation for ambiant temperature |
759 |
|
|
!qcontr(i,k) = Gcontr*(t-Tcontr(i,k)) + qsatliqcontr !--Pa |
760 |
|
|
qcontr = Gcontr*(t-Tcontr) + qsatliqcontr !--Pa |
761 |
|
|
!--Conversion of qcontr in specific humidity - method 1 |
762 |
|
|
!qcontr(i,k) = RD/RV*qcontr(i,k)/pplay !--so as to return to something similar to R2ES*FOEEW/pplay |
763 |
|
|
qcontr2 = RD/RV*qcontr/pplay !--so as to return to something similar to R2ES*FOEEW/pplay |
764 |
|
|
!qcontr(i,k) = min(0.5,qcontr(i,k)) !--and then we apply the same correction term as for qsat |
765 |
|
|
qcontr2 = min(0.5,qcontr2) !--and then we apply the same correction term as for qsat |
766 |
|
|
!zcor = 1./(1.-RETV*qcontr(i,k)) !--for consistency with qsat but is it correct at all? |
767 |
|
|
zcor = 1./(1.-RETV*qcontr2) !--for consistency with qsat but is it correct at all as p is dry? |
768 |
|
|
!zcor = 1./(1.+qcontr2) !--for consistency with qsat |
769 |
|
|
!qcontr(i,k) = qcontr(i,k)*zcor |
770 |
|
|
qcontr2 = qcontr2*zcor |
771 |
|
|
qcontr2=MAX(1.e-10,qcontr2) !--eliminate negative values due to extrapolation on dilution curve |
772 |
|
|
!--Conversion of qcontr in specific humidity - method 2 |
773 |
|
|
!qcontr(i,k) = eps_w*qcontr(i,k) / (pplay+eps_w*qcontr(i,k)) |
774 |
|
|
!qcontr=MAX(1.E-10,qcontr) |
775 |
|
|
!qcontr2 = eps_w*qcontr / (pplay+eps_w*qcontr) |
776 |
|
|
! |
777 |
|
|
IF (t .LT. Tcontr) THEN !--contrail formation is possible |
778 |
|
|
! |
779 |
|
|
!--compute fractions of persistent (P) and non-persistent(N) contrail potential regions |
780 |
|
|
!!IF (qcontr(i,k).GE.qsat) THEN |
781 |
|
|
IF (qcontr2.GE.qsat) THEN |
782 |
|
|
!--none of the unsaturated clear sky is prone for contrail formation |
783 |
|
|
!!fcontrN(i,k) = 0.0 |
784 |
|
|
fcontrN = 0.0 |
785 |
|
|
! |
786 |
|
|
!--integral of P(q) from qsat to qcontr in ISSR |
787 |
|
|
pdf_a = log(qsat/q)/(pdf_k*sqrt(2.)) |
788 |
|
|
pdf_e1 = pdf_a+pdf_b |
789 |
|
|
IF (abs(pdf_e1).GE.erf_lim) THEN |
790 |
|
|
pdf_e1 = sign(1.,pdf_e1) |
791 |
|
|
ELSE |
792 |
|
|
pdf_e1 = erf(pdf_e1) |
793 |
|
|
ENDIF |
794 |
|
|
! |
795 |
|
|
!!pdf_a = log(MIN(qcontr(i,k),qvc)/q)/(pdf_k*sqrt(2.)) |
796 |
|
|
pdf_a = log(MIN(qcontr2,qvc)/q)/(pdf_k*sqrt(2.)) |
797 |
|
|
pdf_e2 = pdf_a+pdf_b |
798 |
|
|
IF (abs(pdf_e2).GE.erf_lim) THEN |
799 |
|
|
pdf_e2 = sign(1.,pdf_e2) |
800 |
|
|
ELSE |
801 |
|
|
pdf_e2 = erf(pdf_e2) |
802 |
|
|
ENDIF |
803 |
|
|
! |
804 |
|
|
!!fcontrP(i,k) = MAX(0., 0.5*(pdf_e1-pdf_e2)) |
805 |
|
|
fcontrP = MAX(0., 0.5*(pdf_e1-pdf_e2)) |
806 |
|
|
! |
807 |
|
|
pdf_a = log(qsat/q)/(pdf_k*sqrt(2.)) |
808 |
|
|
pdf_e1 = pdf_a+pdf_b |
809 |
|
|
IF (abs(pdf_e1).GE.erf_lim) THEN |
810 |
|
|
pdf_e1 = sign(1.,pdf_e1) |
811 |
|
|
ELSE |
812 |
|
|
pdf_e1 = erf(pdf_e1) |
813 |
|
|
ENDIF |
814 |
|
|
! |
815 |
|
|
!!pdf_a = log(MIN(qcontr(i,k),qvc)/q)/(pdf_k*sqrt(2.)) |
816 |
|
|
pdf_a = log(MIN(qcontr2,qvc)/q)/(pdf_k*sqrt(2.)) |
817 |
|
|
pdf_e2 = pdf_a+pdf_b |
818 |
|
|
IF (abs(pdf_e2).GE.erf_lim) THEN |
819 |
|
|
pdf_e2 = sign(1.,pdf_e2) |
820 |
|
|
ELSE |
821 |
|
|
pdf_e2 = erf(pdf_e2) |
822 |
|
|
ENDIF |
823 |
|
|
! |
824 |
|
|
!!fcontrP(i,k) = MAX(0., 0.5*(pdf_e1-pdf_e2)) |
825 |
|
|
fcontrP = MAX(0., 0.5*(pdf_e1-pdf_e2)) |
826 |
|
|
! |
827 |
|
|
pdf_a = log(MAX(qsat,qvc)/q)/(pdf_k*sqrt(2.)) |
828 |
|
|
pdf_e1 = pdf_a+pdf_b |
829 |
|
|
IF (abs(pdf_e1).GE.erf_lim) THEN |
830 |
|
|
pdf_e1 = sign(1.,pdf_e1) |
831 |
|
|
ELSE |
832 |
|
|
pdf_e1 = erf(pdf_e1) |
833 |
|
|
ENDIF |
834 |
|
|
! |
835 |
|
|
!!pdf_a = log(MIN(qcontr(i,k),MIN(gamma_prec*qvc,gamma_ss*qsat))/q)/(pdf_k*sqrt(2.)) |
836 |
|
|
pdf_a = log(MIN(qcontr2,MIN(gamma_prec*qvc,gamma_ss*qsat))/q)/(pdf_k*sqrt(2.)) |
837 |
|
|
pdf_e2 = pdf_a+pdf_b |
838 |
|
|
IF (abs(pdf_e2).GE.erf_lim) THEN |
839 |
|
|
pdf_e2 = sign(1.,pdf_e2) |
840 |
|
|
ELSE |
841 |
|
|
pdf_e2 = erf(pdf_e2) |
842 |
|
|
ENDIF |
843 |
|
|
! |
844 |
|
|
!!fcontrP(i,k) = fcontrP(i,k) + MAX(0., 0.5*(1-pdf_N1)*(pdf_e1-pdf_e2)) |
845 |
|
|
fcontrP = fcontrP + MAX(0., 0.5*(1-pdf_N1)*(pdf_e1-pdf_e2)) |
846 |
|
|
! |
847 |
|
|
pdf_a = log(gamma_prec*qvc/q)/(pdf_k*sqrt(2.)) |
848 |
|
|
pdf_e1 = pdf_a+pdf_b |
849 |
|
|
IF (abs(pdf_e1).GE.erf_lim) THEN |
850 |
|
|
pdf_e1 = sign(1.,pdf_e1) |
851 |
|
|
ELSE |
852 |
|
|
pdf_e1 = erf(pdf_e1) |
853 |
|
|
ENDIF |
854 |
|
|
! |
855 |
|
|
!!pdf_a = log(MIN(qcontr(i,k),gamma_ss*qsat)/q)/(pdf_k*sqrt(2.)) |
856 |
|
|
pdf_a = log(MIN(qcontr2,gamma_ss*qsat)/q)/(pdf_k*sqrt(2.)) |
857 |
|
|
pdf_e2 = pdf_a+pdf_b |
858 |
|
|
IF (abs(pdf_e2).GE.erf_lim) THEN |
859 |
|
|
pdf_e2 = sign(1.,pdf_e2) |
860 |
|
|
ELSE |
861 |
|
|
pdf_e2 = erf(pdf_e2) |
862 |
|
|
ENDIF |
863 |
|
|
! |
864 |
|
|
!!fcontrP(i,k) = fcontrP(i,k) + MAX(0., 0.5*(1-pdf_N2)*(pdf_e1-pdf_e2)) |
865 |
|
|
fcontrP = fcontrP + MAX(0., 0.5*(1-pdf_N2)*(pdf_e1-pdf_e2)) |
866 |
|
|
! |
867 |
|
|
ELSE !--qcontr LT qsat |
868 |
|
|
! |
869 |
|
|
!--all of ISSR is prone for contrail formation |
870 |
|
|
!!fcontrP(i,k) = rnebss |
871 |
|
|
fcontrP = rnebss |
872 |
|
|
! |
873 |
|
|
!--integral of zq from qcontr to qsat in unsaturated clear-sky region |
874 |
|
|
!!pdf_a = log(qcontr(i,k)/q)/(pdf_k*sqrt(2.)) |
875 |
|
|
pdf_a = log(qcontr2/q)/(pdf_k*sqrt(2.)) |
876 |
|
|
pdf_e1 = pdf_a+pdf_b !--normalement pdf_b est deja defini |
877 |
|
|
IF (abs(pdf_e1).GE.erf_lim) THEN |
878 |
|
|
pdf_e1 = sign(1.,pdf_e1) |
879 |
|
|
ELSE |
880 |
|
|
pdf_e1 = erf(pdf_e1) |
881 |
|
|
ENDIF |
882 |
|
|
! |
883 |
|
|
pdf_a = log(qsat/q)/(pdf_k*sqrt(2.)) |
884 |
|
|
pdf_e2 = pdf_a+pdf_b |
885 |
|
|
IF (abs(pdf_e2).GE.erf_lim) THEN |
886 |
|
|
pdf_e2 = sign(1.,pdf_e2) |
887 |
|
|
ELSE |
888 |
|
|
pdf_e2 = erf(pdf_e2) |
889 |
|
|
ENDIF |
890 |
|
|
! |
891 |
|
|
!!fcontrN(i,k) = 0.5*(pdf_e1-pdf_e2) |
892 |
|
|
fcontrN = 0.5*(pdf_e1-pdf_e2) |
893 |
|
|
!!fcontrN=2.0 |
894 |
|
|
! |
895 |
|
|
ENDIF |
896 |
|
|
! |
897 |
|
|
ENDIF !-- t < Tcontr |
898 |
|
|
! |
899 |
|
|
ENDIF !-- Gcontr > 0.1 |
900 |
|
|
|
901 |
|
|
RETURN |
902 |
|
|
END SUBROUTINE ice_sursat |
903 |
|
|
! |
904 |
|
|
!******************************************************************* |
905 |
|
|
! |
906 |
|
|
END MODULE ice_sursat_mod |