source: LMDZ6/branches/Amaury_dev/libf/phylmd/ice_sursat_mod.F90 @ 5143

Last change on this file since 5143 was 5143, checked in by abarral, 3 months ago

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