Context Navigation

ice_sursat_mod.F90 @ 4142

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

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

Download in other formats:

Original Format