Context Navigation

ice_sursat_mod.F90 @ 5087

Last change on this file since 5087 was 5084, checked in by Laurent Fairhead, 12 months ago
Reverting to r4065. Updating fortran standard broke too much stuff. Will do it by smaller chunks AB, LF
File size: 32.8 KB

Line
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)/RNAVO44.e-3cell_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)/RNAVO18.e-3cell_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_cellq + 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.qvc3 + 7.PI*2/480.qvc**5 &
433	+ 127.PI3/40320.qvc*7 + 4369.PI*4/5806080.qvc**9 &
434	+ 34807.PI5/182476800.qvc**11
435	qvc = sqrt(PI)/2.*qvc
436	qvc = (qvc-pdf_b)pdf_ksqrt(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(qvcgamma_prec/q)/(pdf_ksqrt(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.5pdf_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_ssqsat.LT.gamma_precqvc) THEN
519
520	pdf_a = log(gamma_ssqsat/q)/(pdf_ksqrt(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_precqvc/q)/(pdf_ksqrt(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_ssqsat/q)/(pdf_ksqrt(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.5q(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_ssqsat,gamma_precqvc)/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.5qpdf_N1*(pdf_e1-pdf_e2)
601
602	qcld = pdf_e1
603
604	! Partie 2 (sous condition)
605	IF (gamma_ssqsat.GT.gamma_precqvc) THEN
606	pdf_a = log(gamma_ssqsat/q)/(pdf_ksqrt(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.5qpdf_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.5q(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 = (rnebV_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_cld4.PI(3.(b_tur*2.)L_tur + L_tur*3. + 3.b_tur(L_tur*2.))
649	V_cld = N_cld4.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(sigV_env,rnebssV_cell)
655	V_clr = MIN((1.-sig)V_env,rnebclrV_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_envqclr/MAX(eps,rnebclr) + V_cldqcld/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) = rnebssflight_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_wQheat(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.43log(Gcontr-0.053)+0.72(log(Gcontr-0.053))**2 !--K
751	IF (Gcontr .GT. 0.1) THEN
752	!
753	Tcontr = 226.69+9.43log(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/RVqcontr(i,k)/pplay !--so as to return to something similar to R2ESFOEEW/pplay
763	qcontr2 = RD/RVqcontr/pplay !--so as to return to something similar to R2ESFOEEW/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_wqcontr(i,k) / (pplay+eps_wqcontr(i,k))
774	!qcontr=MAX(1.E-10,qcontr)
775	!qcontr2 = eps_wqcontr / (pplay+eps_wqcontr)
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_precqvc,gamma_ssqsat))/q)/(pdf_k*sqrt(2.))
836	pdf_a = log(MIN(qcontr2,MIN(gamma_precqvc,gamma_ssqsat))/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_precqvc/q)/(pdf_ksqrt(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_ssqsat)/q)/(pdf_ksqrt(2.))
856	pdf_a = log(MIN(qcontr2,gamma_ssqsat)/q)/(pdf_ksqrt(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

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

Download in other formats:

Original Format