Context Navigation

ice_sursat_mod.F90 @ 4461

Last change on this file since 4461 was 4260, checked in by lguez, 2 years ago

Bug fix: split too long lines

The Fortran 2003 standard says a line may contain no more than 132
characters.

File size: 32.9 KB

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

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

Download in other formats:

Original Format