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

Last change on this file since 5119 was 5112, checked in by abarral, 4 months ago

Rename modules in phy_common from *_mod > lmdz_*
File size: 32.5 KB
Line 
1MODULE ice_sursat_mod
2
3IMPLICIT NONE
4
5!--flight inventories
6
7REAL, SAVE, ALLOCATABLE :: flight_m(:,:)    !--flown distance m s-1 per cell
8!$OMP THREADPRIVATE(flight_m)
9REAL, 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
15REAL, 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
50REAL, SAVE :: chi=1.1, l_turb=50.0, tun_N=1.3, tun_ratqs=3.0
51REAL, 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
58SUBROUTINE 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
84END SUBROUTINE ice_sursat_init
85
86!*******************************************************************
87
88SUBROUTINE 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
260END SUBROUTINE airplane
261
262!********************************************************************
263! simple routine to initialise flight_m and test a flight corridor
264!--Olivier Boucher - 2021
265
266SUBROUTINE 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
285END SUBROUTINE flight_init
286
287!*******************************************************************
288!--Routine to deal with ice supersaturation
289!--Determines the respective fractions of unsaturated clear sky, ice supersaturated clear sky and cloudy sky
290!--Diagnoses regions prone for non-persistent and persistent contrail formation
291
292!--Audran Borella - 2021
293
294SUBROUTINE ice_sursat(pplay, dpaprs, dtime, i, k, t, q, gamma_ss, &
295                      qsat, t_actuel, rneb_seri, ratqs, rneb, qincld,   &
296                      rnebss, qss, Tcontr, qcontr, qcontr2, fcontrN, fcontrP)
297
298  USE dimphy
299  USE lmdz_print_control,    ONLY: prt_level, lunout
300  USE phys_state_var_mod,   ONLY: pbl_tke, t_ancien
301  USE phys_local_var_mod,   ONLY: N1_ss, N2_ss
302  USE phys_local_var_mod,   ONLY: drneb_sub, drneb_con, drneb_tur, drneb_avi
303!!  USE phys_local_var_mod,   ONLY: Tcontr, qcontr, fcontrN, fcontrP
304  USE indice_sol_mod,       ONLY: is_ave
305  USE lmdz_geometry,         ONLY: cell_area
306
307  IMPLICIT NONE
308  INCLUDE "YOMCST.h"
309  INCLUDE "YOETHF.h"
310  INCLUDE "FCTTRE.h"
311  INCLUDE "clesphys.h"
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_cell*q + 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)>=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>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.*qvc**3 + 7.*PI**2/480.*qvc**5 &
433             + 127.*PI**3/40320.*qvc**7 + 4369.*PI**4/5806080.*qvc**9 &
434             + 34807.*PI**5/182476800.*qvc**11
435        qvc = sqrt(PI)/2.*qvc
436        qvc = (qvc-pdf_b)*pdf_k*sqrt(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)>=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(qvc*gamma_prec/q)/(pdf_k*sqrt(2.))
462        pdf_e2 = pdf_a+pdf_b
463        IF (abs(pdf_e2)>=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)<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<=0.) THEN
478           pdf_N1 = 0.
479           IF (pdf_e2>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<qsat) THEN
490        pdf_a = log(qvc/q)/(pdf_k*sqrt(2.))
491        pdf_e1 = pdf_a+pdf_b
492        IF (abs(pdf_e1)>=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)>=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.5*pdf_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_ss*qsat<gamma_prec*qvc) THEN
519     
520        pdf_a = log(gamma_ss*qsat/q)/(pdf_k*sqrt(2.))
521        pdf_e1 = pdf_a+pdf_b
522        IF (abs(pdf_e1)>=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_prec*qvc/q)/(pdf_k*sqrt(2.))
529        pdf_e2 = pdf_a+pdf_b
530        IF (abs(pdf_e2)>=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_ss*qsat/q)/(pdf_k*sqrt(2.))
546        pdf_e1 = pdf_a+pdf_b
547        IF (abs(pdf_e1)>=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)>=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)>=erf_lim) THEN
573        pdf_e2 = sign(1.,pdf_e2)
574     ELSE
575        pdf_e2 = erf(pdf_e2)
576     ENDIF
577     pdf_e2 = 0.5*q*(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)>=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_ss*qsat,gamma_prec*qvc)/q)/(pdf_k*sqrt(2.))
593     pdf_e2 = pdf_a-pdf_b
594     IF (abs(pdf_e2)>=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.5*q*pdf_N1*(pdf_e1-pdf_e2)
601     
602     qcld = pdf_e1
603
604     ! Partie 2 (sous condition)
605     IF (gamma_ss*qsat>gamma_prec*qvc) THEN
606        pdf_a = log(gamma_ss*qsat/q)/(pdf_k*sqrt(2.))
607        pdf_e1 = pdf_a-pdf_b
608        IF (abs(pdf_e1)>=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.5*q*pdf_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.5*q*(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>eps.AND.rneb<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 = (rneb*V_cell/4./PI/N_cld)**(1./3.)
643       ! On verifie que la longeur de melange n'est pas trop grande
644       IF (L_tur>b_tur) THEN
645          L_tur = b_tur
646       ENDIF
647       
648       V_env = N_cld*4.*PI*(3.*(b_tur**2.)*L_tur + L_tur**3. + 3.*b_tur*(L_tur**2.))
649       V_cld = N_cld*4.*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(sig*V_env,rnebss*V_cell)
655       V_clr = MIN((1.-sig)*V_env,rnebclr*V_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_env*qclr/MAX(eps,rnebclr) + V_cld*qcld/MAX(eps,rneb)
665       q_eq = q_eq/(V_env + V_cld)
666
667       IF (q_eq>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) = rnebss*flight_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<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<eps) THEN
714        ! s'il y a une ISSR on met tout dans l'ISSR, sinon dans le
715        ! clear sky
716        IF (rnebss<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_w*Qheat*(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.43*log(Gcontr-0.053)+0.72*(log(Gcontr-0.053))**2 !--K
751     IF (Gcontr > 0.1) THEN
752
753       Tcontr = 226.69+9.43*log(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/RV*qcontr(i,k)/pplay      !--so as to return to something similar to R2ES*FOEEW/pplay
763       qcontr2 = RD/RV*qcontr/pplay      !--so as to return to something similar to R2ES*FOEEW/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_w*qcontr(i,k) / (pplay+eps_w*qcontr(i,k))
774       !qcontr=MAX(1.E-10,qcontr)
775       !qcontr2 = eps_w*qcontr / (pplay+eps_w*qcontr)
776
777       IF (t < 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>=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)>=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)>=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)>=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)>=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)>=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_prec*qvc,gamma_ss*qsat))/q)/(pdf_k*sqrt(2.))
836         pdf_a = log(MIN(qcontr2,MIN(gamma_prec*qvc,gamma_ss*qsat))/q)/(pdf_k*sqrt(2.))
837         pdf_e2 = pdf_a+pdf_b
838         IF (abs(pdf_e2)>=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_prec*qvc/q)/(pdf_k*sqrt(2.))
848         pdf_e1 = pdf_a+pdf_b
849         IF (abs(pdf_e1)>=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_ss*qsat)/q)/(pdf_k*sqrt(2.))
856         pdf_a = log(MIN(qcontr2,gamma_ss*qsat)/q)/(pdf_k*sqrt(2.))
857         pdf_e2 = pdf_a+pdf_b
858         IF (abs(pdf_e2)>=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)>=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)>=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
902END SUBROUTINE ice_sursat
903
904!*******************************************************************
905
906END MODULE ice_sursat_mod
Note: See TracBrowser for help on using the repository browser.