Context Navigation

ice_sursat_mod.F90 @ 4074

Last change on this file since 4074 was 4074, checked in by oboucher, 3 years ago
only compute contrail diagnostics in the right range of Gcontr
File size: 30.6 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format