Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

ice_sursat_mod.F90 @ 4582

Last change on this file since 4582 was 4368, checked in by lguez, 2 years ago
Sync latest trunk changes to Ocean_skin
File size: 32.9 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 "fisrtilp.h"
311	INCLUDE "clesphys.h"
312
313	!
314	! Input
315	! Beware: this routine works on a gridpoint!
316	!
317	REAL, INTENT(IN) :: pplay ! layer pressure (Pa)
318	REAL, INTENT(IN) :: dpaprs ! layer delta pressure (Pa)
319	REAL, INTENT(IN) :: dtime ! intervalle du temps (s)
320	REAL, INTENT(IN) :: t ! température advectée (K)
321	REAL, INTENT(IN) :: qsat ! vapeur de saturation
322	REAL, INTENT(IN) :: t_actuel ! temperature actuelle de la maille (K)
323	REAL, INTENT(IN) :: rneb_seri ! fraction nuageuse en memoire
324	INTEGER, INTENT(IN) :: i, k
325	!
326	! Input/output
327	!
328	REAL, INTENT(INOUT) :: q ! vapeur de la maille (=zq)
329	REAL, INTENT(INOUT) :: ratqs ! determine la largeur de distribution de vapeur
330	REAL, INTENT(INOUT) :: Tcontr, qcontr, qcontr2, fcontrN, fcontrP
331	!
332	! Output
333	!
334	REAL, INTENT(OUT) :: gamma_ss !
335	REAL, INTENT(OUT) :: rneb ! cloud fraction
336	REAL, INTENT(OUT) :: qincld ! in-cloud total water
337	REAL, INTENT(OUT) :: rnebss ! ISSR fraction
338	REAL, INTENT(OUT) :: qss ! in-ISSR total water
339	!
340	! Local
341	!
342	REAL PI
343	PARAMETER (PI=4.*ATAN(1.))
344	REAL rnebclr, gamma_prec
345	REAL qclr, qvc, qcld, qi
346	REAL zrho, zdz, zrhodz
347	REAL pdf_N, pdf_N1, pdf_N2
348	REAL pdf_a, pdf_b
349	REAL pdf_e1, pdf_e2, pdf_k
350	REAL drnebss, drnebclr, dqss, dqclr, sum_rneb_rnebss, dqss_avi
351	REAL V_cell !--volume of the cell
352	REAL M_cell !--dry mass of the cell
353	REAL tke, sig, L_tur, b_tur, q_eq
354	REAL V_env, V_cld, V_ss, V_clr
355	REAL zcor
356	!
357	!--more local variables for diagnostics
358	!--imported from YOMCST.h
359	!--eps_w = 0.622 = ratio of molecular masses of water and dry air (kg H2O kg air -1)
360	!--RCPD = 1004 J kg air−1 K−1 = the isobaric heat capacity of air
361	!--values from Schumann, Meteorol Zeitschrift, 1996
362	!--EiH2O = 1.25 / 2.24 / 8.94 kg H2O / kg fuel for kerosene / methane / dihydrogen
363	!--Qheat = 43. / 50. / 120. MJ / kg fuel for kerosene / methane / dihydrogen
364	REAL, PARAMETER :: EiH2O=1.25 !--emission index of water vapour for kerosene (kg kg-1)
365	REAL, PARAMETER :: Qheat=43.E6 !--specific combustion heat for kerosene (J kg-1)
366	REAL, PARAMETER :: eta=0.3 !--average propulsion efficiency of the aircraft
367	!--Gcontr is the slope of the mean phase trajectory in the turbulent exhaust field on an absolute
368	!--temperature versus water vapor partial pressure diagram. G has the unit of Pa K−1. Rap et al JGR 2010.
369	REAL :: Gcontr
370	!--Tcontr = critical temperature for contrail formation (T_LM in Schumann 1996, Eq 31 in appendix 2)
371	!--qsatliqcontr = e_L(T_LM) in Schumann 1996 but expressed in specific humidity (kg kg humid air-1)
372	REAL :: qsatliqcontr
373
374	! Initialisations
375	zrho = pplay / t / RD !--dry density kg m-3
376	zrhodz = dpaprs / RG !--dry air mass kg m-2
377	zdz = zrhodz / zrho !--cell thickness m
378	V_cell = zdz * cell_area(i) !--cell volume m3
379	M_cell = zrhodz * cell_area(i) !--cell dry air mass kg
380	!
381	! Recuperation de la memoire sur la couverture nuageuse
382	rneb = rneb_seri
383	!
384	! Ajout des émissions de H2O dues à l'aviation
385	! q is the specific humidity (kg/kg humid air) hence the complicated equation to update q
386	! qnew = ( m_humid_air * qold + dm_H2O ) / ( m_humid_air + dm_H2O )
387	! = ( m_dry_air * qold + dm_h2O * (1-qold) ) / (m_dry_air + dm_H2O * (1-qold) )
388	! The equation is derived by writing m_humid_air = m_dry_air + m_H2O = m_dry_air / (1-q)
389	! flight_h2O is in kg H2O / s / cell
390	!
391	IF (ok_plane_h2o) THEN
392	q = ( M_cellq + flight_h2o(i,k)dtime(1.-q) ) / (M_cell + flight_h2o(i,k)dtime*(1.-q) )
393	ENDIF
394	!
395	!--Estimating gamma
396	gamma_ss = MAX(1.0, gamma0 - t_actuel/Tgamma)
397	!gamma_prec = MAX(1.0, gamma0 - t_ancien(i,k)/Tgamma) !--formulation initiale d Audran
398	gamma_prec = MAX(1.0, gamma0 - t/Tgamma) !--autre formulation possible basée sur le T du pas de temps
399	!
400	! Initialisation de qvc : q_sat du pas de temps precedent
401	!qvc = R2ES*FOEEW(t_ancien(i,k),1.)/pplay !--formulation initiale d Audran
402	qvc = R2ES*FOEEW(t,1.)/pplay !--autre formulation possible basée sur le T du pas de temps
403	qvc = min(0.5,qvc)
404	zcor = 1./(1.-RETV*qvc)
405	qvc = qvc*zcor
406	!
407	! Modification de ratqs selon formule proposee : ksi_new = ksi_old/(1+beta*alpha_cld)
408	ratqs = ratqs / (tun_ratqs*rneb_seri + 1.)
409	!
410	! Calcul de N
411	pdf_k = -sqrt(log(1.+ratqs**2.))
412	pdf_a = log(qvc/q)/(pdf_k*sqrt(2.))
413	pdf_b = pdf_k/(2.*sqrt(2.))
414	pdf_e1 = pdf_a+pdf_b
415	IF (abs(pdf_e1).GE.erf_lim) THEN
416	pdf_e1 = sign(1.,pdf_e1)
417	pdf_N = max(0.,sign(rneb,pdf_e1))
418	ELSE
419	pdf_e1 = erf(pdf_e1)
420	pdf_e1 = 0.5*(1.+pdf_e1)
421	pdf_N = max(0.,rneb/pdf_e1)
422	ENDIF
423	!
424	! On calcule ensuite N1 et N2. Il y a deux cas : N > 1 et N <= 1
425	! Cas 1 : N > 1. N'arrive en theorie jamais, c'est une barriere
426	! On perd la memoire sur la temperature (sur qvc) pour garder
427	! celle sur alpha_cld
428	IF (pdf_N.GT.1.) THEN
429	! On inverse alpha_cld = int_qvc^infty P(q) dq
430	! pour determiner qvc = f(alpha_cld)
431	! On approxime en serie entiere erf-1(x)
432	qvc = 2.*rneb-1.
433	qvc = qvc + PI/12.qvc3 + 7.PI*2/480.qvc**5 &
434	+ 127.PI3/40320.qvc*7 + 4369.PI*4/5806080.qvc**9 &
435	+ 34807.PI5/182476800.qvc**11
436	qvc = sqrt(PI)/2.*qvc
437	qvc = (qvc-pdf_b)pdf_ksqrt(2.)
438	qvc = q*exp(qvc)
439
440	! On met a jour rneb avec la nouvelle valeur de qvc
441	! La maj est necessaire a cause de la serie entiere approximative
442	pdf_a = log(qvc/q)/(pdf_k*sqrt(2.))
443	pdf_e1 = pdf_a+pdf_b
444	IF (abs(pdf_e1).GE.erf_lim) THEN
445	pdf_e1 = sign(1.,pdf_e1)
446	ELSE
447	pdf_e1 = erf(pdf_e1)
448	ENDIF
449	pdf_e1 = 0.5*(1.+pdf_e1)
450	rneb = pdf_e1
451
452	! Si N > 1, N1 et N2 = 1
453	pdf_N1 = 1.
454	pdf_N2 = 1.
455
456	! Cas 2 : N <= 1
457	ELSE
458	! D'abord on calcule N2 avec le tuning
459	pdf_N2 = min(1.,pdf_N*tun_N)
460
461	! Puis N1 pour assurer la conservation de alpha_cld
462	pdf_a = log(qvcgamma_prec/q)/(pdf_ksqrt(2.))
463	pdf_e2 = pdf_a+pdf_b
464	IF (abs(pdf_e2).GE.erf_lim) THEN
465	pdf_e2 = sign(1.,pdf_e2)
466	ELSE
467	pdf_e2 = erf(pdf_e2)
468	ENDIF
469	pdf_e2 = 0.5*(1.+pdf_e2) ! integrale sous P pour q > gamma qsat
470
471	IF (abs(pdf_e1-pdf_e2).LT.eps) THEN
472	pdf_N1 = pdf_N2
473	ELSE
474	pdf_N1 = (rneb-pdf_N2*pdf_e2)/(pdf_e1-pdf_e2)
475	ENDIF
476
477	! Barriere qui traite le cas gamma_prec = 1.
478	IF (pdf_N1.LE.0.) THEN
479	pdf_N1 = 0.
480	IF (pdf_e2.GT.eps) THEN
481	pdf_N2 = rneb/pdf_e2
482	ELSE
483	pdf_N2 = 0.
484	ENDIF
485	ENDIF
486	ENDIF
487
488	! Physique 1
489	! Sublimation
490	IF (qvc.LT.qsat) THEN
491	pdf_a = log(qvc/q)/(pdf_k*sqrt(2.))
492	pdf_e1 = pdf_a+pdf_b
493	IF (abs(pdf_e1).GE.erf_lim) THEN
494	pdf_e1 = sign(1.,pdf_e1)
495	ELSE
496	pdf_e1 = erf(pdf_e1)
497	ENDIF
498
499	pdf_a = log(qsat/q)/(pdf_k*sqrt(2.))
500	pdf_e2 = pdf_a+pdf_b
501	IF (abs(pdf_e2).GE.erf_lim) THEN
502	pdf_e2 = sign(1.,pdf_e2)
503	ELSE
504	pdf_e2 = erf(pdf_e2)
505	ENDIF
506
507	pdf_e1 = 0.5pdf_N1(pdf_e1-pdf_e2)
508
509	! Calcul et ajout de la tendance
510	drneb_sub(i,k) = - pdf_e1/dtime !--unit [s-1]
511	rneb = rneb + drneb_sub(i,k)*dtime
512	ELSE
513	drneb_sub(i,k) = 0.
514	ENDIF
515
516	! NOTE : verifier si ca marche bien pour gamma proche de 1.
517
518	! Condensation
519	IF (gamma_ssqsat.LT.gamma_precqvc) THEN
520
521	pdf_a = log(gamma_ssqsat/q)/(pdf_ksqrt(2.))
522	pdf_e1 = pdf_a+pdf_b
523	IF (abs(pdf_e1).GE.erf_lim) THEN
524	pdf_e1 = sign(1.,pdf_e1)
525	ELSE
526	pdf_e1 = erf(pdf_e1)
527	ENDIF
528
529	pdf_a = log(gamma_precqvc/q)/(pdf_ksqrt(2.))
530	pdf_e2 = pdf_a+pdf_b
531	IF (abs(pdf_e2).GE.erf_lim) THEN
532	pdf_e2 = sign(1.,pdf_e2)
533	ELSE
534	pdf_e2 = erf(pdf_e2)
535	ENDIF
536
537	pdf_e1 = 0.5(1.-pdf_N1)(pdf_e1-pdf_e2)
538	pdf_e2 = 0.5(1.-pdf_N2)(1.+pdf_e2)
539
540	! Calcul et ajout de la tendance
541	drneb_con(i,k) = (pdf_e1 + pdf_e2)/dtime !--unit [s-1]
542	rneb = rneb + drneb_con(i,k)*dtime
543
544	ELSE
545
546	pdf_a = log(gamma_ssqsat/q)/(pdf_ksqrt(2.))
547	pdf_e1 = pdf_a+pdf_b
548	IF (abs(pdf_e1).GE.erf_lim) THEN
549	pdf_e1 = sign(1.,pdf_e1)
550	ELSE
551	pdf_e1 = erf(pdf_e1)
552	ENDIF
553	pdf_e1 = 0.5(1.-pdf_N2)(1.+pdf_e1)
554
555	! Calcul et ajout de la tendance
556	drneb_con(i,k) = pdf_e1/dtime !--unit [s-1]
557	rneb = rneb + drneb_con(i,k)*dtime
558
559	ENDIF
560
561	! Calcul des grandeurs diagnostiques
562	! Determination des grandeurs ciel clair
563	pdf_a = log(qsat/q)/(pdf_k*sqrt(2.))
564	pdf_e1 = pdf_a+pdf_b
565	IF (abs(pdf_e1).GE.erf_lim) THEN
566	pdf_e1 = sign(1.,pdf_e1)
567	ELSE
568	pdf_e1 = erf(pdf_e1)
569	ENDIF
570	pdf_e1 = 0.5*(1.-pdf_e1)
571
572	pdf_e2 = pdf_a-pdf_b
573	IF (abs(pdf_e2).GE.erf_lim) THEN
574	pdf_e2 = sign(1.,pdf_e2)
575	ELSE
576	pdf_e2 = erf(pdf_e2)
577	ENDIF
578	pdf_e2 = 0.5q(1.-pdf_e2)
579
580	rnebclr = pdf_e1
581	qclr = pdf_e2
582
583	! Determination de q_cld
584	! Partie 1
585	pdf_a = log(max(qsat,qvc)/q)/(pdf_k*sqrt(2.))
586	pdf_e1 = pdf_a-pdf_b
587	IF (abs(pdf_e1).GE.erf_lim) THEN
588	pdf_e1 = sign(1.,pdf_e1)
589	ELSE
590	pdf_e1 = erf(pdf_e1)
591	ENDIF
592
593	pdf_a = log(min(gamma_ssqsat,gamma_precqvc)/q)/(pdf_k*sqrt(2.))
594	pdf_e2 = pdf_a-pdf_b
595	IF (abs(pdf_e2).GE.erf_lim) THEN
596	pdf_e2 = sign(1.,pdf_e2)
597	ELSE
598	pdf_e2 = erf(pdf_e2)
599	ENDIF
600
601	pdf_e1 = 0.5qpdf_N1*(pdf_e1-pdf_e2)
602
603	qcld = pdf_e1
604
605	! Partie 2 (sous condition)
606	IF (gamma_ssqsat.GT.gamma_precqvc) THEN
607	pdf_a = log(gamma_ssqsat/q)/(pdf_ksqrt(2.))
608	pdf_e1 = pdf_a-pdf_b
609	IF (abs(pdf_e1).GE.erf_lim) THEN
610	pdf_e1 = sign(1.,pdf_e1)
611	ELSE
612	pdf_e1 = erf(pdf_e1)
613	ENDIF
614
615	pdf_e2 = 0.5qpdf_N2*(pdf_e2-pdf_e1)
616
617	qcld = qcld + pdf_e2
618
619	pdf_e2 = pdf_e1
620	ENDIF
621
622	! Partie 3
623	pdf_e2 = 0.5q(1.+pdf_e2)
624
625	qcld = qcld + pdf_e2
626
627	! Fin du calcul de q_cld
628
629	! Determination des grandeurs ISSR via les equations de conservation
630	rneb=MIN(rneb, 1. - rnebclr - eps) !--ajout OB - barrière
631	rnebss = MAX(0.0, 1. - rnebclr - rneb) !--ajout OB
632	qss = MAX(0.0, q - qclr - qcld) !--ajout OB
633
634	! Physique 2 : Turbulence
635	IF (rneb.GT.eps.AND.rneb.LT.1.-eps) THEN ! rneb != 0 and != 1
636	!
637	tke = pbl_tke(i,k,is_ave)
638	! A MODIFIER formule a revoir
639	L_tur = min(l_turb, sqrt(tke)*dtime)
640
641	! On fait pour l'instant l'hypothese a = 3b. V = 4/3 pi a b**2 = alpha_cld
642	! donc b = alpha_cld/4pi **1/3.
643	b_tur = (rnebV_cell/4./PI/N_cld)*(1./3.)
644	! On verifie que la longeur de melange n'est pas trop grande
645	IF (L_tur.GT.b_tur) THEN
646	L_tur = b_tur
647	ENDIF
648
649	V_env = N_cld4.PI(3.(b_tur*2.)L_tur + L_tur*3. + 3.b_tur(L_tur*2.))
650	V_cld = N_cld4.PI(3.(b_tur*2.)L_tur + L_tur*3. - 3.b_tur(L_tur*2.))
651	V_cld = abs(V_cld)
652
653	! Repartition statistique des zones de contact nuage-ISSR et nuage-ciel clair
654	sig = rnebss/(chi*rnebclr+rnebss)
655	V_ss = MIN(sigV_env,rnebssV_cell)
656	V_clr = MIN((1.-sig)V_env,rnebclrV_cell)
657	V_cld = MIN(V_cld,rneb*V_cell)
658	V_env = V_ss + V_clr
659
660	! ISSR => rneb
661	drnebss = -1.*V_ss/V_cell
662	dqss = drnebss*qss/MAX(eps,rnebss)
663
664	! Clear sky <=> rneb
665	q_eq = V_envqclr/MAX(eps,rnebclr) + V_cldqcld/MAX(eps,rneb)
666	q_eq = q_eq/(V_env + V_cld)
667
668	IF (q_eq.GT.qsat) THEN
669	drnebclr = - V_clr/V_cell
670	dqclr = drnebclr*qclr/MAX(eps,rnebclr)
671	ELSE
672	drnebclr = V_cld/V_cell
673	dqclr = drnebclr*qcld/MAX(eps,rneb)
674	ENDIF
675
676	! Maj des variables avec les tendances
677	rnebclr = MAX(0.0,rnebclr + drnebclr) !--OB ajout d'un max avec eps (il faudrait modified drnebclr pour le diag)
678	qclr = MAX(0.0, qclr + dqclr) !--OB ajout d'un max avec 0
679
680	rneb = rneb - drnebclr - drnebss
681	qcld = qcld - dqclr - dqss
682
683	rnebss = MAX(0.0,rnebss + drnebss) !--OB ajout d'un max avec eps (il faudrait modifier drnebss pour le diag)
684	qss = MAX(0.0, qss + dqss) !--OB ajout d'un max avec 0
685
686	! Tendances pour le diagnostic
687	drneb_tur(i,k) = (drnebclr + drnebss)/dtime !--unit [s-1]
688
689	ENDIF ! rneb
690
691	!--add a source of cirrus from aviation contrails
692	IF (ok_plane_contrail) THEN
693	drneb_avi(i,k) = rnebssflight_m(i,k)contrail_cross_section/V_cell !--tendency rneb due to aviation [s-1]
694	drneb_avi(i,k) = MIN(drneb_avi(i,k), rnebss/dtime) !--majoration
695	dqss_avi = qss*drneb_avi(i,k)/MAX(eps,rnebss) !--tendency q aviation [kg kg-1 s-1]
696	rneb = rneb + drneb_avi(i,k)*dtime !--add tendency to rneb
697	qcld = qcld + dqss_avi*dtime !--add tendency to qcld
698	rnebss = rnebss - drneb_avi(i,k)*dtime !--add tendency to rnebss
699	qss = qss - dqss_avi*dtime !--add tendency to qss
700	ELSE
701	drneb_avi(i,k)=0.0
702	ENDIF
703
704	! Barrieres
705	! ISSR trop petite
706	IF (rnebss.LT.eps) THEN
707	rneb = MIN(rneb + rnebss,1.0-eps) !--ajout OB barriere
708	qcld = qcld + qss
709	rnebss = 0.
710	qss = 0.
711	ENDIF
712
713	! le nuage est trop petit
714	IF (rneb.LT.eps) THEN
715	! s'il y a une ISSR on met tout dans l'ISSR, sinon dans le
716	! clear sky
717	IF (rnebss.LT.eps) THEN
718	rnebclr = 1.
719	rnebss = 0. !--ajout OB
720	qclr = q
721	ELSE
722	rnebss = MIN(rnebss + rneb,1.0-eps) !--ajout OB barriere
723	qss = qss + qcld
724	ENDIF
725	rneb = 0.
726	qcld = 0.
727	qincld = qsat * gamma_ss
728	ELSE
729	qincld = qcld / rneb
730	ENDIF
731
732	!--OB ajout borne superieure
733	sum_rneb_rnebss=rneb+rnebss
734	rneb=rneb*MIN(1.-eps,sum_rneb_rnebss)/MAX(eps,sum_rneb_rnebss)
735	rnebss=rnebss*MIN(1.-eps,sum_rneb_rnebss)/MAX(eps,sum_rneb_rnebss)
736
737	! On ecrit dans la memoire
738	N1_ss(i,k) = pdf_N1
739	N2_ss(i,k) = pdf_N2
740
741	!--Diagnostics only used from last iteration
742	!--test
743	!!Tcontr(i,k)=200.
744	!!fcontrN(i,k)=1.0
745	!!fcontrP(i,k)=0.5
746	!
747	!--slope of dilution line in exhaust
748	!--kg H2O/kg fuel * J kg air-1 K-1 * Pa / (kg H2O / kg air * J kg fuel-1)
749	Gcontr = EiH2O * RCPD * pplay / (eps_wQheat(1.-eta)) !--Pa K-1
750	!--critical T_LM below which no liquid contrail can form in exhaust
751	!Tcontr(i,k) = 226.69+9.43log(Gcontr-0.053)+0.72(log(Gcontr-0.053))**2 !--K
752	IF (Gcontr .GT. 0.1) THEN
753	!
754	Tcontr = 226.69+9.43log(Gcontr-0.053)+0.72(log(Gcontr-0.053))**2 !--K
755	!print *,'Tcontr=',iter,i,k,eps_w,pplay,Gcontr,Tcontr(i,k)
756	!--Psat with index 0 in FOEEW to get saturation wrt liquid water corresponding to Tcontr
757	!qsatliqcontr = RESTT*FOEEW(Tcontr(i,k),0.) !--Pa
758	qsatliqcontr = RESTT*FOEEW(Tcontr,0.) !--Pa
759	!--Critical water vapour above which there is contrail formation for ambiant temperature
760	!qcontr(i,k) = Gcontr*(t-Tcontr(i,k)) + qsatliqcontr !--Pa
761	qcontr = Gcontr*(t-Tcontr) + qsatliqcontr !--Pa
762	!--Conversion of qcontr in specific humidity - method 1
763	!qcontr(i,k) = RD/RVqcontr(i,k)/pplay !--so as to return to something similar to R2ESFOEEW/pplay
764	qcontr2 = RD/RVqcontr/pplay !--so as to return to something similar to R2ESFOEEW/pplay
765	!qcontr(i,k) = min(0.5,qcontr(i,k)) !--and then we apply the same correction term as for qsat
766	qcontr2 = min(0.5,qcontr2) !--and then we apply the same correction term as for qsat
767	!zcor = 1./(1.-RETV*qcontr(i,k)) !--for consistency with qsat but is it correct at all?
768	zcor = 1./(1.-RETV*qcontr2) !--for consistency with qsat but is it correct at all as p is dry?
769	!zcor = 1./(1.+qcontr2) !--for consistency with qsat
770	!qcontr(i,k) = qcontr(i,k)*zcor
771	qcontr2 = qcontr2*zcor
772	qcontr2=MAX(1.e-10,qcontr2) !--eliminate negative values due to extrapolation on dilution curve
773	!--Conversion of qcontr in specific humidity - method 2
774	!qcontr(i,k) = eps_wqcontr(i,k) / (pplay+eps_wqcontr(i,k))
775	!qcontr=MAX(1.E-10,qcontr)
776	!qcontr2 = eps_wqcontr / (pplay+eps_wqcontr)
777	!
778	IF (t .LT. Tcontr) THEN !--contrail formation is possible
779	!
780	!--compute fractions of persistent (P) and non-persistent(N) contrail potential regions
781	!!IF (qcontr(i,k).GE.qsat) THEN
782	IF (qcontr2.GE.qsat) THEN
783	!--none of the unsaturated clear sky is prone for contrail formation
784	!!fcontrN(i,k) = 0.0
785	fcontrN = 0.0
786	!
787	!--integral of P(q) from qsat to qcontr in ISSR
788	pdf_a = log(qsat/q)/(pdf_k*sqrt(2.))
789	pdf_e1 = pdf_a+pdf_b
790	IF (abs(pdf_e1).GE.erf_lim) THEN
791	pdf_e1 = sign(1.,pdf_e1)
792	ELSE
793	pdf_e1 = erf(pdf_e1)
794	ENDIF
795	!
796	!!pdf_a = log(MIN(qcontr(i,k),qvc)/q)/(pdf_k*sqrt(2.))
797	pdf_a = log(MIN(qcontr2,qvc)/q)/(pdf_k*sqrt(2.))
798	pdf_e2 = pdf_a+pdf_b
799	IF (abs(pdf_e2).GE.erf_lim) THEN
800	pdf_e2 = sign(1.,pdf_e2)
801	ELSE
802	pdf_e2 = erf(pdf_e2)
803	ENDIF
804	!
805	!!fcontrP(i,k) = MAX(0., 0.5*(pdf_e1-pdf_e2))
806	fcontrP = MAX(0., 0.5*(pdf_e1-pdf_e2))
807	!
808	pdf_a = log(qsat/q)/(pdf_k*sqrt(2.))
809	pdf_e1 = pdf_a+pdf_b
810	IF (abs(pdf_e1).GE.erf_lim) THEN
811	pdf_e1 = sign(1.,pdf_e1)
812	ELSE
813	pdf_e1 = erf(pdf_e1)
814	ENDIF
815	!
816	!!pdf_a = log(MIN(qcontr(i,k),qvc)/q)/(pdf_k*sqrt(2.))
817	pdf_a = log(MIN(qcontr2,qvc)/q)/(pdf_k*sqrt(2.))
818	pdf_e2 = pdf_a+pdf_b
819	IF (abs(pdf_e2).GE.erf_lim) THEN
820	pdf_e2 = sign(1.,pdf_e2)
821	ELSE
822	pdf_e2 = erf(pdf_e2)
823	ENDIF
824	!
825	!!fcontrP(i,k) = MAX(0., 0.5*(pdf_e1-pdf_e2))
826	fcontrP = MAX(0., 0.5*(pdf_e1-pdf_e2))
827	!
828	pdf_a = log(MAX(qsat,qvc)/q)/(pdf_k*sqrt(2.))
829	pdf_e1 = pdf_a+pdf_b
830	IF (abs(pdf_e1).GE.erf_lim) THEN
831	pdf_e1 = sign(1.,pdf_e1)
832	ELSE
833	pdf_e1 = erf(pdf_e1)
834	ENDIF
835	!
836	!!pdf_a = log(MIN(qcontr(i,k),MIN(gamma_precqvc,gamma_ssqsat))/q)/(pdf_k*sqrt(2.))
837	pdf_a = log(MIN(qcontr2,MIN(gamma_precqvc,gamma_ssqsat))/q)/(pdf_k*sqrt(2.))
838	pdf_e2 = pdf_a+pdf_b
839	IF (abs(pdf_e2).GE.erf_lim) THEN
840	pdf_e2 = sign(1.,pdf_e2)
841	ELSE
842	pdf_e2 = erf(pdf_e2)
843	ENDIF
844	!
845	!!fcontrP(i,k) = fcontrP(i,k) + MAX(0., 0.5(1-pdf_N1)(pdf_e1-pdf_e2))
846	fcontrP = fcontrP + MAX(0., 0.5(1-pdf_N1)(pdf_e1-pdf_e2))
847	!
848	pdf_a = log(gamma_precqvc/q)/(pdf_ksqrt(2.))
849	pdf_e1 = pdf_a+pdf_b
850	IF (abs(pdf_e1).GE.erf_lim) THEN
851	pdf_e1 = sign(1.,pdf_e1)
852	ELSE
853	pdf_e1 = erf(pdf_e1)
854	ENDIF
855	!
856	!!pdf_a = log(MIN(qcontr(i,k),gamma_ssqsat)/q)/(pdf_ksqrt(2.))
857	pdf_a = log(MIN(qcontr2,gamma_ssqsat)/q)/(pdf_ksqrt(2.))
858	pdf_e2 = pdf_a+pdf_b
859	IF (abs(pdf_e2).GE.erf_lim) THEN
860	pdf_e2 = sign(1.,pdf_e2)
861	ELSE
862	pdf_e2 = erf(pdf_e2)
863	ENDIF
864	!
865	!!fcontrP(i,k) = fcontrP(i,k) + MAX(0., 0.5(1-pdf_N2)(pdf_e1-pdf_e2))
866	fcontrP = fcontrP + MAX(0., 0.5(1-pdf_N2)(pdf_e1-pdf_e2))
867	!
868	ELSE !--qcontr LT qsat
869	!
870	!--all of ISSR is prone for contrail formation
871	!!fcontrP(i,k) = rnebss
872	fcontrP = rnebss
873	!
874	!--integral of zq from qcontr to qsat in unsaturated clear-sky region
875	!!pdf_a = log(qcontr(i,k)/q)/(pdf_k*sqrt(2.))
876	pdf_a = log(qcontr2/q)/(pdf_k*sqrt(2.))
877	pdf_e1 = pdf_a+pdf_b !--normalement pdf_b est deja defini
878	IF (abs(pdf_e1).GE.erf_lim) THEN
879	pdf_e1 = sign(1.,pdf_e1)
880	ELSE
881	pdf_e1 = erf(pdf_e1)
882	ENDIF
883	!
884	pdf_a = log(qsat/q)/(pdf_k*sqrt(2.))
885	pdf_e2 = pdf_a+pdf_b
886	IF (abs(pdf_e2).GE.erf_lim) THEN
887	pdf_e2 = sign(1.,pdf_e2)
888	ELSE
889	pdf_e2 = erf(pdf_e2)
890	ENDIF
891	!
892	!!fcontrN(i,k) = 0.5*(pdf_e1-pdf_e2)
893	fcontrN = 0.5*(pdf_e1-pdf_e2)
894	!!fcontrN=2.0
895	!
896	ENDIF
897	!
898	ENDIF !-- t < Tcontr
899	!
900	ENDIF !-- Gcontr > 0.1
901
902	RETURN
903	END SUBROUTINE ice_sursat
904	!
905	!*******************************************************************
906	!
907	END MODULE ice_sursat_mod

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

Download in other formats:

Original Format