source: LMDZ6/trunk/libf/phylmd/ice_sursat_mod.F90 @ 4142

Last change on this file since 4142 was 4099, checked in by oboucher, 3 years ago

update to ice_sursat module

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