Context Navigation

lmdz_lscp_tools.f90 @ 5763

Last change on this file since 5763 was 5727, checked in by evignon, 2 weeks ago
legere modification de lscp pour enlever un warning obsolete sur la parameterisation de phase
File size: 38.8 KB

Rev	Line
[4665]	1	MODULE lmdz_lscp_tools
[3999]	2
	3	IMPLICIT NONE
	4
	5	CONTAINS
	6
	7	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
	8	SUBROUTINE FALLICE_VELOCITY(klon,iwc,temp,rho,pres,ptconv,velo)
	9
	10	! Ref:
	11	! Stubenrauch, C. J., Bonazzola, M.,
	12	! Protopapadaki, S. E., & Musat, I. (2019).
	13	! New cloud system metrics to assess bulk
	14	! ice cloud schemes in a GCM. Journal of
	15	! Advances in Modeling Earth Systems, 11,
	16	! 3212–3234. https://doi.org/10.1029/2019MS001642
	17
[4665]	18	use lmdz_lscp_ini, only: iflag_vice, ffallv_con, ffallv_lsc
	19	use lmdz_lscp_ini, only: cice_velo, dice_velo
[4535]	20
[3999]	21	IMPLICIT NONE
	22
	23	INTEGER, INTENT(IN) :: klon
	24	REAL, INTENT(IN), DIMENSION(klon) :: iwc ! specific ice water content [kg/m3]
	25	REAL, INTENT(IN), DIMENSION(klon) :: temp ! temperature [K]
	26	REAL, INTENT(IN), DIMENSION(klon) :: rho ! dry air density [kg/m3]
	27	REAL, INTENT(IN), DIMENSION(klon) :: pres ! air pressure [Pa]
	28	LOGICAL, INTENT(IN), DIMENSION(klon) :: ptconv ! convective point [-]
	29
	30	REAL, INTENT(OUT), DIMENSION(klon) :: velo ! fallspeed velocity of crystals [m/s]
	31
	32
	33	INTEGER i
[4559]	34	REAL logvm,iwcg,tempc,phpa,fallv_tun
[3999]	35	REAL m2ice, m2snow, vmice, vmsnow
	36	REAL aice, bice, asnow, bsnow
	37
	38
	39	DO i=1,klon
	40
	41	IF (ptconv(i)) THEN
	42	fallv_tun=ffallv_con
	43	ELSE
	44	fallv_tun=ffallv_lsc
	45	ENDIF
	46
	47	tempc=temp(i)-273.15 ! celcius temp
[4072]	48	iwcg=MAX(iwc(i)*1000.,1E-3) ! iwc in g/m3. We set a minimum value to prevent from division by 0
[3999]	49	phpa=pres(i)/100. ! pressure in hPa
	50
	51	IF (iflag_vice .EQ. 1) THEN
	52	! so-called 'empirical parameterization' in Stubenrauch et al. 2019
	53	if (tempc .GE. -60.0) then
	54
	55	logvm= -0.0000414122tempctempclog(iwcg)-0.00538922tempc*log(iwcg) &
	56	-0.0516344log(iwcg)+0.00216078tempc + 1.9714
	57	velo(i)=exp(logvm)
	58	else
	59	velo(i)=65.0(iwcg0.2)(150./phpa)**0.15
	60	endif
	61
	62	velo(i)=fallv_tun*velo(i)/100.0 ! from cm/s to m/s
	63
	64	ELSE IF (iflag_vice .EQ. 2) THEN
	65	! so called PSDM empirical coherent bulk ice scheme in Stubenrauch et al. 2019
	66	aice=0.587
	67	bice=2.45
	68	asnow=0.0444
	69	bsnow=2.1
	70
	71	m2ice=((iwcg0.001/aice)/(exp(13.6-bice7.76+0.479bice2) &
	72	exp((-0.0361+bice0.0151+0.00149bice*2)tempc))) &
	73	*(1./(0.807+bice0.00581+0.0457bice*2))
	74
[4072]	75	vmice=100.1042.4exp(13.6-(bice+1)7.76+0.479(bice+1.)*2)exp((-0.0361+ &
	76	(bice+1.)0.0151+0.00149(bice+1.)*2)tempc) &
	77	(m2ice(0.807+(bice+1.)0.00581+0.0457(bice+1.)2))/(iwcg0.001/aice)
[3999]	78
	79
	80	vmice=vmice((1000./phpa)*0.2)
	81
	82	m2snow=((iwcg0.001/asnow)/(exp(13.6-bsnow7.76+0.479bsnow2) &
	83	exp((-0.0361+bsnow0.0151+0.00149bsnow*2)tempc))) &
	84	*(1./(0.807+bsnow0.00581+0.0457bsnow*2))
	85
	86
	87	vmsnow=100.14.3exp(13.6-(bsnow+.416)7.76+0.479(bsnow+.416)**2)&
	88	exp((-0.0361+(bsnow+.416)0.0151+0.00149(bsnow+.416)2)tempc)&
	89	(m2snow(0.807+(bsnow+.416)0.00581+0.0457(bsnow+.416)2))/(iwcg0.001/asnow)
	90
	91	vmsnow=vmsnow((1000./phpa)*0.35)
	92	velo(i)=fallv_tun*min(vmsnow,vmice)/100. ! to m/s
	93
	94	ELSE
	95	! By default, fallspeed velocity of ice crystals according to Heymsfield & Donner 1990
[4559]	96	velo(i) = fallv_tuncice_velo((iwcg/1000.)**dice_velo)
[3999]	97	ENDIF
	98	ENDDO
	99
	100	END SUBROUTINE FALLICE_VELOCITY
	101	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
	102
	103	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
[4639]	104	SUBROUTINE ICEFRAC_LSCP(klon, temp, iflag_ice_thermo, distcltop, temp_cltop, icefrac, dicefracdT)
[3999]	105	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
	106
	107	! Compute the ice fraction 1-xliq (see e.g.
	108	! Doutriaux-Boucher & Quaas 2004, section 2.2.)
	109	! as a function of temperature
	110	! see also Fig 3 of Madeleine et al. 2020, JAMES
	111	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
	112
	113
	114	USE print_control_mod, ONLY: lunout, prt_level
[4665]	115	USE lmdz_lscp_ini, ONLY: t_glace_min, t_glace_max, exposant_glace, iflag_t_glace
	116	USE lmdz_lscp_ini, ONLY : RTT, dist_liq, temp_nowater
[3999]	117
[4059]	118	IMPLICIT NONE
[3999]	119
	120
[4639]	121	INTEGER, INTENT(IN) :: klon ! number of horizontal grid points
	122	REAL, INTENT(IN), DIMENSION(klon) :: temp ! temperature
	123	REAL, INTENT(IN), DIMENSION(klon) :: distcltop ! distance to cloud top
	124	REAL, INTENT(IN), DIMENSION(klon) :: temp_cltop ! temperature of cloud top
[4535]	125	INTEGER, INTENT(IN) :: iflag_ice_thermo
[3999]	126	REAL, INTENT(OUT), DIMENSION(klon) :: icefrac
	127	REAL, INTENT(OUT), DIMENSION(klon) :: dicefracdT
	128
	129
	130	INTEGER i
[4562]	131	REAL liqfrac_tmp, dicefrac_tmp
[3999]	132	REAL Dv, denomdep,beta,qsi,dqsidt
	133	LOGICAL ice_thermo
	134
[4562]	135	CHARACTER (len = 20) :: modname = 'lscp_tools'
	136	CHARACTER (len = 80) :: abort_message
	137
[5007]	138	IF ((iflag_t_glace.LT.2)) THEN !.OR. (iflag_t_glace.GT.6)) THEN
[4562]	139	abort_message = 'lscp cannot be used if iflag_t_glace<2 or >6'
	140	CALL abort_physic(modname,abort_message,1)
[3999]	141	ENDIF
	142
[4562]	143	IF (.NOT.((iflag_ice_thermo .EQ. 1).OR.(iflag_ice_thermo .GE. 3))) THEN
	144	abort_message = 'lscp cannot be used without ice thermodynamics'
	145	CALL abort_physic(modname,abort_message,1)
	146	ENDIF
[3999]	147
	148
	149	DO i=1,klon
[4562]	150
	151	! old function with sole dependence upon temperature
	152	IF (iflag_t_glace .EQ. 2) THEN
	153	liqfrac_tmp = (temp(i)-t_glace_min) / (t_glace_max-t_glace_min)
	154	liqfrac_tmp = MIN(MAX(liqfrac_tmp,0.0),1.0)
[3999]	155	icefrac(i) = (1.0-liqfrac_tmp)**exposant_glace
	156	IF (icefrac(i) .GT.0.) THEN
	157	dicefracdT(i)= exposant_glace * (icefrac(i)**(exposant_glace-1.)) &
	158	/ (t_glace_min - t_glace_max)
	159	ENDIF
	160
	161	IF ((icefrac(i).EQ.0).OR.(icefrac(i).EQ.1)) THEN
	162	dicefracdT(i)=0.
	163	ENDIF
	164
[4562]	165	ENDIF
[3999]	166
[4562]	167	! function of temperature used in CMIP6 physics
	168	IF (iflag_t_glace .EQ. 3) THEN
	169	liqfrac_tmp = (temp(i)-t_glace_min) / (t_glace_max-t_glace_min)
	170	liqfrac_tmp = MIN(MAX(liqfrac_tmp,0.0),1.0)
	171	icefrac(i) = 1.0-liqfrac_tmp**exposant_glace
	172	IF ((icefrac(i) .GT.0.) .AND. (liqfrac_tmp .GT. 0.)) THEN
	173	dicefracdT(i)= exposant_glace * ((liqfrac_tmp)**(exposant_glace-1.)) &
	174	/ (t_glace_min - t_glace_max)
	175	ELSE
	176	dicefracdT(i)=0.
	177	ENDIF
	178	ENDIF
[4072]	179
[4562]	180	! for iflag_t_glace .GE. 4, the liquid fraction depends upon temperature at cloud top
	181	! and then decreases with decreasing height
[3999]	182
[4562]	183	!with linear function of temperature at cloud top
	184	IF (iflag_t_glace .EQ. 4) THEN
	185	liqfrac_tmp = (temp(i)-t_glace_min) / (t_glace_max-t_glace_min)
	186	liqfrac_tmp = MIN(MAX(liqfrac_tmp,0.0),1.0)
	187	icefrac(i) = MAX(MIN(1.,1.0 - liqfrac_tmp*exp(-distcltop(i)/dist_liq)),0.)
	188	dicefrac_tmp = - temp(i)/(t_glace_max-t_glace_min)
	189	dicefracdT(i) = dicefrac_tmp*exp(-distcltop(i)/dist_liq)
	190	IF ((liqfrac_tmp .LE.0) .OR. (liqfrac_tmp .GE. 1)) THEN
	191	dicefracdT(i) = 0.
	192	ENDIF
	193	ENDIF
	194
	195	! with CMIP6 function of temperature at cloud top
[5007]	196	IF ((iflag_t_glace .EQ. 5) .OR. (iflag_t_glace .EQ. 7)) THEN
[4562]	197	liqfrac_tmp = (temp(i)-t_glace_min) / (t_glace_max-t_glace_min)
	198	liqfrac_tmp = MIN(MAX(liqfrac_tmp,0.0),1.0)
	199	liqfrac_tmp = liqfrac_tmp**exposant_glace
	200	icefrac(i) = MAX(MIN(1.,1.0 - liqfrac_tmp*exp(-distcltop(i)/dist_liq)),0.)
	201	IF ((liqfrac_tmp .LE.0) .OR. (liqfrac_tmp .GE. 1)) THEN
	202	dicefracdT(i) = 0.
	203	ELSE
	204	dicefracdT(i) = exposant_glace((liqfrac_tmp)*(exposant_glace-1.))/(t_glace_min- t_glace_max) &
	205	*exp(-distcltop(i)/dist_liq)
	206	ENDIF
	207	ENDIF
	208
	209	! with modified function of temperature at cloud top
	210	! to get largere values around 260 K, works well with t_glace_min = 241K
	211	IF (iflag_t_glace .EQ. 6) THEN
	212	IF (temp(i) .GT. t_glace_max) THEN
	213	liqfrac_tmp = 1.
	214	ELSE
	215	liqfrac_tmp = -((temp(i)-t_glace_max) / (t_glace_max-t_glace_min))**2+1.
	216	ENDIF
	217	liqfrac_tmp = MIN(MAX(liqfrac_tmp,0.0),1.0)
	218	icefrac(i) = MAX(MIN(1.,1.0 - liqfrac_tmp*exp(-distcltop(i)/dist_liq)),0.)
	219	IF ((liqfrac_tmp .LE.0) .OR. (liqfrac_tmp .GE. 1)) THEN
	220	dicefracdT(i) = 0.
	221	ELSE
	222	dicefracdT(i) = 2*((temp(i)-t_glace_max) / (t_glace_max-t_glace_min))/(t_glace_max-t_glace_min) &
	223	*exp(-distcltop(i)/dist_liq)
	224	ENDIF
	225	ENDIF
	226
[5614]	227	! if temperature or temperature of cloud top <-40°C,
[4639]	228	IF (iflag_t_glace .GE. 4) THEN
	229	IF ((temp_cltop(i) .LE. temp_nowater) .AND. (temp(i) .LE. t_glace_max)) THEN
	230	icefrac(i) = 1.
	231	dicefracdT(i) = 0.
	232	ENDIF
	233	ENDIF
[5007]	234
[4562]	235
	236	ENDDO ! klon
	237	RETURN
	238
[3999]	239	END SUBROUTINE ICEFRAC_LSCP
	240	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
	241
[5383]	242
[5727]	243	SUBROUTINE ICEFRAC_LSCP_TURB(klon, dtime, pticefracturb, temp, pplay, paprsdn, paprsup, wvel, qice_ini, snowcld, qtot_incl, cldfra, tke, &
[5614]	244	tke_dissip, sursat_e, invtau_e, qliq, qvap_cld, qice, icefrac, dicefracdT, cldfraliq, sigma2_icefracturb, mean_icefracturb)
[5007]	245	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
	246	! Compute the liquid, ice and vapour content (+ice fraction) based
[5170]	247	! on turbulence (see Fields 2014, Furtado 2016, Raillard 2025)
[5383]	248	! L.Raillard (23/09/24)
[5614]	249	! E.Vignon (03/2025) : additional elements for treatment of convective
	250	! boundary layer clouds
[5007]	251	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
[3999]	252
	253
[5007]	254	USE lmdz_lscp_ini, ONLY : prt_level, lunout
	255	USE lmdz_lscp_ini, ONLY : RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT, RD, RG, RV, RPI
	256	USE lmdz_lscp_ini, ONLY : seuil_neb, temp_nowater
[5614]	257	USE lmdz_lscp_ini, ONLY : naero5, gamma_snwretro, gamma_taud, capa_crystal, rho_ice
[5633]	258	USE lmdz_lscp_ini, ONLY : eps, snow_fallspeed
[5007]	259
	260	IMPLICIT NONE
	261
[5383]	262	INTEGER, INTENT(IN) :: klon !--number of horizontal grid points
	263	REAL, INTENT(IN) :: dtime !--time step [s]
[5614]	264	LOGICAL, INTENT(IN), DIMENSION(klon) :: pticefracturb !--grid points concerned by this routine
[5383]	265	REAL, INTENT(IN), DIMENSION(klon) :: temp !--temperature
	266	REAL, INTENT(IN), DIMENSION(klon) :: pplay !--pressure in the middle of the layer [Pa]
	267	REAL, INTENT(IN), DIMENSION(klon) :: paprsdn !--pressure at the bottom interface of the layer [Pa]
	268	REAL, INTENT(IN), DIMENSION(klon) :: paprsup !--pressure at the top interface of the layer [Pa]
[5614]	269	REAL, INTENT(IN), DIMENSION(klon) :: wvel !--vertical velocity [m/s]
[5383]	270	REAL, INTENT(IN), DIMENSION(klon) :: qtot_incl !--specific total cloud water in-cloud content [kg/kg]
	271	REAL, INTENT(IN), DIMENSION(klon) :: cldfra !--cloud fraction in gridbox [-]
	272	REAL, INTENT(IN), DIMENSION(klon) :: tke !--turbulent kinetic energy [m2/s2]
	273	REAL, INTENT(IN), DIMENSION(klon) :: tke_dissip !--TKE dissipation [m2/s3]
[5007]	274
[5383]	275	REAL, INTENT(IN), DIMENSION(klon) :: qice_ini !--initial specific ice content gridbox-mean [kg/kg]
[5614]	276	REAL, INTENT(IN), DIMENSION(klon) :: snowcld !--in-cloud snowfall flux [kg/m2/s]
	277	REAL, INTENT(IN), DIMENSION(klon) :: sursat_e !--environment supersaturation [-]
	278	REAL, INTENT(IN), DIMENSION(klon) :: invtau_e !--inverse time-scale of mixing with environment [s-1]
[5383]	279	REAL, INTENT(OUT), DIMENSION(klon) :: qliq !--specific liquid content gridbox-mean [kg/kg]
	280	REAL, INTENT(OUT), DIMENSION(klon) :: qvap_cld !--specific cloud vapor content, gridbox-mean [kg/kg]
	281	REAL, INTENT(OUT), DIMENSION(klon) :: qice !--specific ice content gridbox-mean [kg/kg]
[5007]	282
[5614]	283	REAL, INTENT(INOUT), DIMENSION(klon) :: icefrac !--fraction of ice in condensed water [-]
	284	REAL, INTENT(INOUT), DIMENSION(klon) :: dicefracdT
	285
[5383]	286	REAL, INTENT(OUT), DIMENSION(klon) :: cldfraliq !--fraction of cldfra where liquid [-]
	287	REAL, INTENT(OUT), DIMENSION(klon) :: sigma2_icefracturb !--Sigma2 of the ice supersaturation PDF [-]
	288	REAL, INTENT(OUT), DIMENSION(klon) :: mean_icefracturb !--Mean of the ice supersaturation PDF [-]
[5007]	289
[5383]	290	REAL, DIMENSION(klon) :: qzero, qsatl, dqsatl, qsati, dqsati !--specific humidity saturation values
[5007]	291	INTEGER :: i
	292
[5383]	293	REAL :: qvap_incl, qice_incl, qliq_incl, qiceini_incl !--In-cloud specific quantities [kg/kg]
	294	REAL :: water_vapor_diff !--Water-vapour diffusion coeff in air (f(T,P)) [m2/s]
	295	REAL :: air_thermal_conduct !--Thermal conductivity of air (f(T)) [J/m/K/s]
	296	REAL :: C0 !--Lagrangian structure function [-]
[5007]	297	REAL :: tau_dissipturb
[5614]	298	REAL :: invtau_phaserelax
	299	REAL :: sigma2_pdf
[5007]	300	REAL :: ai, bi, B0
	301	REAL :: sursat_iceliq
[5383]	302	REAL :: sursat_equ
[5007]	303	REAL :: liqfra_max
	304	REAL :: sursat_iceext
[5383]	305	REAL :: nb_crystals !--number concentration of ice crystals [#/m3]
	306	REAL :: moment1_PSD !--1st moment of ice PSD
	307	REAL :: N0_PSD, lambda_PSD !--parameters of the exponential PSD
[5007]	308
	309	REAL :: cldfra1D
[5383]	310	REAL :: rho_air
	311	REAL :: psati !--saturation vapor pressure wrt ice [Pa]
[5614]	312
[5383]	313
[5614]	314	REAL :: tempvig1, tempvig2
	315
	316	tempvig1 = -21.06 + RTT
	317	tempvig2 = -30.35 + RTT
[5383]	318	C0 = 10. !--value assumed in Field2014
[5007]	319	sursat_iceext = -0.1
	320	qzero(:) = 0.
	321	cldfraliq(:) = 0.
[5614]	322	qliq(:) = 0.
	323	qice(:) = 0.
	324	qvap_cld(:) = 0.
[5007]	325	sigma2_icefracturb(:) = 0.
[5383]	326	mean_icefracturb(:) = 0.
[5007]	327
[5383]	328	!--wrt liquid
[5007]	329	CALL calc_qsat_ecmwf(klon,temp(:),qzero(:),pplay(:),RTT,1,.false.,qsatl(:),dqsatl(:))
	330	!--wrt ice
	331	CALL calc_qsat_ecmwf(klon,temp(:),qzero(:),pplay(:),RTT,2,.false.,qsati(:),dqsati(:))
	332
	333
	334	DO i=1,klon
[5383]	335	rho_air = pplay(i) / temp(i) / RD
[5007]	336	! because cldfra is intent in, but can be locally modified due to test
	337	cldfra1D = cldfra(i)
[5614]	338	! activate param for concerned grid points and for cloudy conditions
	339	IF ((pticefracturb(i)) .AND. (cldfra(i) .GT. 0.)) THEN
[5007]	340	IF (cldfra(i) .GE. 1.0) THEN
	341	cldfra1D = 1.0
	342	END IF
	343
	344	! T>0°C, no ice allowed
	345	IF ( temp(i) .GE. RTT ) THEN
	346	qvap_cld(i) = qsatl(i) * cldfra1D
	347	qliq(i) = MAX(0.0,qtot_incl(i)-qsatl(i)) * cldfra1D
	348	qice(i) = 0.
	349	cldfraliq(i) = 1.
	350	icefrac(i) = 0.
	351	dicefracdT(i) = 0.
	352
	353	! T<-38°C, no liquid allowed
	354	ELSE IF ( temp(i) .LE. temp_nowater) THEN
	355	qvap_cld(i) = qsati(i) * cldfra1D
	356	qliq(i) = 0.
	357	qice(i) = MAX(0.0,qtot_incl(i)-qsati(i)) * cldfra1D
	358	cldfraliq(i) = 0.
	359	icefrac(i) = 1.
	360	dicefracdT(i) = 0.
	361
[5614]	362
[5383]	363	!---------------------------------------------------------
	364	!-- MIXED PHASE TEMPERATURE REGIME
	365	!---------------------------------------------------------
	366	!--In the mixed phase regime (-38°C< T <0°C) we distinguish
	367	!--3 possible subcases.
	368	!--1. No pre-existing ice
	369	!--2A. Pre-existing ice and no turbulence
	370	!--2B. Pre-existing ice and turbulence
[5007]	371	ELSE
[5383]	372
[5633]	373	! gamma_snwretro controls the contribution of snowflakes to the negative feedback
	374	! note that for reasons related to inetarctions with the condensation iteration in lscp_main
	375	! we consider here the mean snowflake concentration in the mesh (not the in-cloud concentration)
	376	! when poprecip is active, it will be worth testing considering the incloud fraction, dividing
[5727]	377	! by znebprecipcld
[5633]	378	! qiceini_incl = qice_ini(i) / cldfra1D + &
	379	! gamma_snwretro * snowcld(i) * RG * dtime / ( paprsdn(i) - paprsup(i) )
	380	! assuming constant snowfall velocity
	381	qiceini_incl = qice_ini(i) / cldfra1D + gamma_snwretro * snowcld(i) / pplay(i) * RD * temp(i) / snow_fallspeed
[5383]	382
[5614]	383	!--1. No preexisting ice and no mixing with environment: if vertical motion, fully liquid
[5383]	384	!--cloud else fully iced cloud
[5614]	385	IF ( (qiceini_incl .LT. eps) .AND. (invtau_e(i) .LT. eps) ) THEN
	386	IF ( (wvel(i) .GT. eps) .OR. (tke(i) .GT. eps) ) THEN
[5383]	387	qvap_cld(i) = qsatl(i) * cldfra1D
	388	qliq(i) = MAX(0.,qtot_incl(i)-qsatl(i)) * cldfra1D
	389	qice(i) = 0.
	390	cldfraliq(i) = 1.
	391	icefrac(i) = 0.
	392	dicefracdT(i) = 0.
	393	ELSE
	394	qvap_cld(i) = qsati(i) * cldfra1D
	395	qliq(i) = 0.
	396	qice(i) = MAX(0.,qtot_incl(i)-qsati(i)) * cldfra1D
	397	cldfraliq(i) = 0.
	398	icefrac(i) = 1.
	399	dicefracdT(i) = 0.
	400	ENDIF
[5007]	401
	402
[5614]	403	!--2. Pre-existing ice and/or mixing with environment:computation of ice properties for
[5383]	404	!--feedback
	405	ELSE
	406
[5007]	407	sursat_iceliq = qsatl(i)/qsati(i) - 1.
	408	psati = qsati(i) * pplay(i) / (RD/RV)
[5383]	409
[5007]	410	!--We assume an exponential ice PSD whose parameters
	411	!--are computed following Morrison&Gettelman 2008
	412	!--Ice number density is assumed equals to INP density
[5614]	413	!--which is for naero5>0 a function of temperature (DeMott 2010)
[5007]	414	!--bi and B0 are microphysical function characterizing
	415	!--vapor/ice interactions
	416	!--tau_phase_relax is the typical time of vapor deposition
	417	!--onto ice crystals
	418
[5614]	419	!--For naero5<=0 INP density is derived from the empirical fit
	420	!--from MARCUS campaign from Vignon 2021
	421	!--/!\ Note that option is very specific and should be use for
	422	!--the Southern Ocean and the Antarctic
	423
	424	IF (naero5 .LE. 0) THEN
	425	IF ( temp(i) .GT. tempvig1 ) THEN
	426	nb_crystals = 1.e3 * 10*(-0.14(temp(i)-tempvig1) - 2.88)
	427	ELSE IF ( temp(i) .GT. tempvig2 ) THEN
	428	nb_crystals = 1.e3 * 10*(-0.31(temp(i)-tempvig1) - 2.88)
	429	ELSE
	430	nb_crystals = 1.e3 * 10**(0.)
	431	ENDIF
	432	ELSE
	433	nb_crystals = 1.e3 * 5.94e-5 * ( RTT - temp(i) )*3.33 naero5*(0.0264(RTT-temp(i))+0.0033)
	434	ENDIF
	435	lambda_PSD = ( (RPIrho_icenb_crystals) / (rho_air * MAX(qiceini_incl , eps) ) ) ** (1./3.)
[5383]	436	N0_PSD = nb_crystals * lambda_PSD
	437	moment1_PSD = N0_PSD/lambda_PSD**2
	438
[5007]	439	!--Formulae for air thermal conductivity and water vapor diffusivity
	440	!--comes respectively from Beard and Pruppacher (1971)
	441	!--and Hall and Pruppacher (1976)
	442
	443	air_thermal_conduct = ( 5.69 + 0.017 * ( temp(i) - RTT ) ) * 1.e-3 * 4.184
	444	water_vapor_diff = 2.111e-5 ( temp(i) / RTT )*1.94 ( 101325 / pplay(i) )
	445
	446	bi = 1./((qsati(i)+qsatl(i))/2.) + RLSTT2 / RCPD / RV / temp(i)2
	447	B0 = 4. * RPI * capa_crystal * 1. / ( RLSTT2 / air_thermal_conduct / RV / temp(i)2 &
	448	+ RV * temp(i) / psati / water_vapor_diff )
[5614]	449	invtau_phaserelax = bi * B0 * moment1_PSD
[5383]	450
	451	ai = RG / RD / temp(i) * ( RD * RLSTT / RCPD / RV / temp(i) - 1. )
[5614]	452	sursat_equ = (ai * wvel(i) + sursat_e(i)*invtau_e(i)) / (invtau_phaserelax + invtau_e(i))
	453	! as sursaturation is by definition lower than -1 and
	454	! because local supersaturation > 1 are never found in the atmosphere
[5007]	455
[5614]	456	!--2A. No TKE : stationnary binary solution depending on vertical velocity and mixing with env.
[5383]	457	! If Sequ > Siw liquid cloud, else ice cloud
	458	IF ( tke_dissip(i) .LE. eps ) THEN
[5614]	459	sigma2_icefracturb(i)= 0.
	460	mean_icefracturb(i) = sursat_equ
[5383]	461	IF (sursat_equ .GT. sursat_iceliq) THEN
	462	qvap_cld(i) = qsatl(i) * cldfra1D
	463	qliq(i) = MAX(0.,qtot_incl(i)-qsatl(i)) * cldfra1D
	464	qice(i) = 0.
	465	cldfraliq(i) = 1.
	466	icefrac(i) = 0.
	467	dicefracdT(i) = 0.
	468	ELSE
	469	qvap_cld(i) = qsati(i) * cldfra1D
	470	qliq(i) = 0.
	471	qice(i) = MAX(0.,qtot_incl(i)-qsati(i)) * cldfra1D
	472	cldfraliq(i) = 0.
	473	icefrac(i) = 1.
	474	dicefracdT(i) = 0.
	475	ENDIF
	476
	477	!--2B. TKE and ice : ice supersaturation PDF
	478	!--we compute the cloud liquid properties with a Gaussian PDF
	479	!--of ice supersaturation F(Si) (Field2014, Furtado2016).
	480	!--Parameters of the PDF are function of turbulence and
	481	!--microphysics/existing ice.
	482	ELSE
	483
	484	!--Tau_dissipturb is the time needed for turbulence to decay
	485	!--due to viscosity
	486	tau_dissipturb = gamma_taud * 2. * 2./3. * tke(i) / tke_dissip(i) / C0
[5007]	487
[5383]	488	!--------------------- PDF COMPUTATIONS ---------------------
	489	!--Formulae for sigma2_pdf (variance), mean of PDF in Raillard2025
	490	!--cloud liquid fraction and in-cloud liquid content are given
	491	!--by integrating resp. F(Si) and Si*F(Si)
	492	!--Liquid is limited by the available water vapor trough a
	493	!--maximal liquid fraction
	494	!--qice_ini(i) / cldfra1D = qiceincld without precip
[5007]	495
[5383]	496	liqfra_max = MAX(0., (MIN (1.,( qtot_incl(i) - (qice_ini(i) / cldfra1D) - qsati(i) * (1 + sursat_iceext ) ) / ( qsatl(i) - qsati(i) ) ) ) )
[5633]	497	sigma2_pdf = 1./2. * ( ai*2 ) 2./3. * tke(i) * tau_dissipturb / (invtau_phaserelax + invtau_e(i))
[5614]	498	! sursat ranges between -1 and 1, so we prevent sigma2 so exceed 1
	499	cldfraliq(i) = 0.5 * (1. - erf( ( sursat_iceliq - sursat_equ) / (SQRT(2.* sigma2_pdf) ) ) )
[5383]	500	IF (cldfraliq(i) .GT. liqfra_max) THEN
	501	cldfraliq(i) = liqfra_max
	502	ENDIF
	503
[5614]	504	qliq_incl = qsati(i) * SQRT(sigma2_pdf) / SQRT(2.RPI) EXP( -1.(sursat_iceliq - sursat_equ)2. / (2.sigma2_pdf) ) &
	505	- qsati(i) * cldfraliq(i) * (sursat_iceliq - sursat_equ )
[5383]	506
	507	sigma2_icefracturb(i)= sigma2_pdf
[5614]	508	mean_icefracturb(i) = sursat_equ
[5170]	509
[5383]	510	!------------ SPECIFIC VAPOR CONTENT AND WATER CONSERVATION ------------
[5007]	511
[5383]	512	IF ( (qliq_incl .LE. eps) .OR. (cldfraliq(i) .LE. eps) ) THEN
	513	qliq_incl = 0.
	514	cldfraliq(i) = 0.
	515	END IF
	516
	517	!--Specific humidity is the max between qsati and the weighted mean between
	518	!--qv in MPC patches and qv in ice-only parts. We assume that MPC parts are
	519	!--always at qsatl and ice-only parts slightly subsaturated (qsati*sursat_iceext+1)
	520	!--The whole cloud can therefore be supersaturated but never subsaturated.
[5170]	521
[5383]	522	qvap_incl = MAX(qsati(i), ( 1. - cldfraliq(i) ) * (sursat_iceext + 1.) * qsati(i) + cldfraliq(i) * qsatl(i) )
[5170]	523
[5383]	524	IF ( qvap_incl .GE. qtot_incl(i) ) THEN
	525	qvap_incl = qsati(i)
[5614]	526	qliq_incl = MAX(0.0,qtot_incl(i) - qvap_incl)
[5383]	527	qice_incl = 0.
[5170]	528
[5383]	529	ELSEIF ( (qvap_incl + qliq_incl) .GE. qtot_incl(i) ) THEN
	530	qliq_incl = MAX(0.0,qtot_incl(i) - qvap_incl)
	531	qice_incl = 0.
	532	ELSE
	533	qice_incl = qtot_incl(i) - qvap_incl - qliq_incl
	534	END IF
[5007]	535
[5383]	536	qvap_cld(i) = qvap_incl * cldfra1D
	537	qliq(i) = qliq_incl * cldfra1D
	538	qice(i) = qice_incl * cldfra1D
[5614]	539	IF ((qice(i)+qliq(i)) .GT. 0.) THEN
	540	icefrac(i) = qice(i) / ( qice(i) + qliq(i) )
	541	ELSE
	542	icefrac(i) = 1. ! to keep computation of qsat wrt ice in condensation loop in lmdz_lscp_main
	543	ENDIF
[5383]	544	dicefracdT(i) = 0.
[5007]	545
[5383]	546	END IF ! Enough TKE
	547
	548	END IF ! End qini
[5007]	549
	550	END IF ! ! MPC temperature
	551
[5614]	552	END IF ! pticefracturb and cldfra
[5007]	553
	554	ENDDO ! klon
	555	END SUBROUTINE ICEFRAC_LSCP_TURB
[5383]	556	!
[5007]	557	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
	558
	559
[4072]	560	SUBROUTINE CALC_QSAT_ECMWF(klon,temp,qtot,pressure,tref,phase,flagth,qs,dqs)
[3999]	561	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
	562	! Calculate qsat following ECMWF method
	563	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
	564
[4072]	565
[5284]	566	USE yoethf_mod_h
[5285]	567	USE yomcst_mod_h
[5274]	568	IMPLICIT NONE
[3999]	569
[5274]	570
[3999]	571	include "FCTTRE.h"
	572
[4072]	573	INTEGER, INTENT(IN) :: klon ! number of horizontal grid points
	574	REAL, INTENT(IN), DIMENSION(klon) :: temp ! temperature in K
	575	REAL, INTENT(IN), DIMENSION(klon) :: qtot ! total specific water in kg/kg
	576	REAL, INTENT(IN), DIMENSION(klon) :: pressure ! pressure in Pa
	577	REAL, INTENT(IN) :: tref ! reference temperature in K
[3999]	578	LOGICAL, INTENT(IN) :: flagth ! flag for qsat calculation for thermals
	579	INTEGER, INTENT(IN) :: phase
	580	! phase: 0=depend on temperature sign (temp>tref -> liquid, temp<tref, solid)
	581	! 1=liquid
	582	! 2=solid
	583
[4072]	584	REAL, INTENT(OUT), DIMENSION(klon) :: qs ! saturation specific humidity [kg/kg]
	585	REAL, INTENT(OUT), DIMENSION(klon) :: dqs ! derivation of saturation specific humidity wrt T
[3999]	586
	587	REAL delta, cor, cvm5
[4072]	588	INTEGER i
	589
	590	DO i=1,klon
	591
[3999]	592	IF (phase .EQ. 1) THEN
	593	delta=0.
	594	ELSEIF (phase .EQ. 2) THEN
	595	delta=1.
	596	ELSE
[4072]	597	delta=MAX(0.,SIGN(1.,tref-temp(i)))
[3999]	598	ENDIF
	599
	600	IF (flagth) THEN
	601	cvm5=R5LES(1.-delta) + R5IESdelta
	602	ELSE
	603	cvm5 = R5LESRLVTT(1.-delta) + R5IESRLSTTdelta
[4072]	604	cvm5 = cvm5 /RCPD/(1.0+RVTMP2*(qtot(i)))
[3999]	605	ENDIF
	606
[4072]	607	qs(i)= R2ES*FOEEW(temp(i),delta)/pressure(i)
	608	qs(i)=MIN(0.5,qs(i))
	609	cor=1./(1.-RETV*qs(i))
	610	qs(i)=qs(i)*cor
	611	dqs(i)= FOEDE(temp(i),delta,cvm5,qs(i),cor)
[3999]	612
[4072]	613	END DO
	614
[3999]	615	END SUBROUTINE CALC_QSAT_ECMWF
	616	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
	617
	618
	619	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
[4072]	620	SUBROUTINE CALC_GAMMASAT(klon,temp,qtot,pressure,gammasat,dgammasatdt)
[3999]	621
	622	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
	623	! programme that calculates the gammasat parameter that determines the
	624	! homogeneous condensation thresholds for cold (<0oC) clouds
	625	! condensation at q>gammasat*qsat
	626	! Etienne Vignon, March 2021
	627	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
	628
[5204]	629	use lmdz_lscp_ini, only: iflag_gammasat, temp_nowater, RTT
	630	use lmdz_lscp_ini, only: a_homofreez, b_homofreez, delta_hetfreez
[3999]	631
[4059]	632	IMPLICIT NONE
[3999]	633
	634
[4072]	635	INTEGER, INTENT(IN) :: klon ! number of horizontal grid points
	636	REAL, INTENT(IN), DIMENSION(klon) :: temp ! temperature in K
	637	REAL, INTENT(IN), DIMENSION(klon) :: qtot ! total specific water in kg/kg
[3999]	638
[4072]	639	REAL, INTENT(IN), DIMENSION(klon) :: pressure ! pressure in Pa
[3999]	640
[4072]	641	REAL, INTENT(OUT), DIMENSION(klon) :: gammasat ! coefficient to multiply qsat with to calculate saturation
	642	REAL, INTENT(OUT), DIMENSION(klon) :: dgammasatdt ! derivative of gammasat wrt temperature
[3999]	643
[4072]	644	REAL, DIMENSION(klon) :: qsi,qsl,dqsl,dqsi
[5204]	645	REAL f_homofreez, fac
[3999]	646
[4072]	647	INTEGER i
	648
	649	CALL CALC_QSAT_ECMWF(klon,temp,qtot,pressure,RTT,1,.false.,qsl,dqsl)
	650	CALL CALC_QSAT_ECMWF(klon,temp,qtot,pressure,RTT,2,.false.,qsi,dqsi)
[3999]	651
[5204]	652	DO i = 1, klon
[3999]	653
[5204]	654	IF ( temp(i) .GE. RTT ) THEN
[3999]	655	! warm clouds: condensation at saturation wrt liquid
[5204]	656	gammasat(i) = 1.
	657	dgammasatdt(i) = 0.
[3999]	658
[5204]	659	ELSE
	660	! cold clouds: qsi > qsl
[3999]	661
[5204]	662	! homogeneous freezing of aerosols, according to
	663	! Koop, 2000 and Ren and MacKenzie, 2005 (QJRMS)
	664	! 'Cirrus regime'
	665	! if f_homofreez > qsl / qsi, liquid nucleation
	666	! if f_homofreez < qsl / qsi, homogeneous freezing of aerosols
	667	! Note: f_homofreez = qsl / qsi for temp ~= -38degC
	668	f_homofreez = a_homofreez - temp(i) / b_homofreez
	669
	670	IF ( iflag_gammasat .GE. 3 ) THEN
	671	! condensation at homogeneous freezing threshold for temp < -38 degC
	672	! condensation at liquid saturation for temp > -38 degC
	673	IF ( f_homofreez .LE. qsl(i) / qsi(i) ) THEN
	674	gammasat(i) = f_homofreez
	675	dgammasatdt(i) = - 1. / b_homofreez
	676	ELSE
	677	gammasat(i) = qsl(i) / qsi(i)
	678	dgammasatdt(i) = ( dqsl(i) * qsi(i) - dqsi(i) * qsl(i) ) / qsi(i) / qsi(i)
	679	ENDIF
[3999]	680
[5204]	681	ELSEIF ( iflag_gammasat .EQ. 2 ) THEN
	682	! condensation at homogeneous freezing threshold for temp < -38 degC
	683	! condensation at a threshold linearly decreasing between homogeneous
	684	! freezing and ice saturation for -38 degC < temp < temp_nowater
	685	! condensation at ice saturation for temp > temp_nowater
	686	! If temp_nowater = 235.15 K, this is equivalent to iflag_gammasat = 1
	687	IF ( f_homofreez .LE. qsl(i) / qsi(i) ) THEN
	688	gammasat(i) = f_homofreez
	689	dgammasatdt(i) = - 1. / b_homofreez
	690	ELSEIF ( temp(i) .LE. temp_nowater ) THEN
	691	! Here, we assume that f_homofreez = qsl / qsi for temp = -38 degC = 235.15 K
	692	dgammasatdt(i) = ( a_homofreez - 235.15 / b_homofreez - 1. ) &
	693	/ ( 235.15 - temp_nowater )
	694	gammasat(i) = dgammasatdt(i) * ( temp(i) - temp_nowater ) + 1.
	695	ELSE
	696	gammasat(i) = 1.
	697	dgammasatdt(i) = 0.
	698	ENDIF
[3999]	699
[5204]	700	ELSEIF ( iflag_gammasat .EQ. 1 ) THEN
	701	! condensation at homogeneous freezing threshold for temp < -38 degC
	702	! condensation at ice saturation for temp > -38 degC
	703	IF ( f_homofreez .LE. qsl(i) / qsi(i) ) THEN
	704	gammasat(i) = f_homofreez
	705	dgammasatdt(i) = - 1. / b_homofreez
	706	ELSE
	707	gammasat(i) = 1.
	708	dgammasatdt(i) = 0.
	709	ENDIF
	710
[3999]	711	ELSE
[5204]	712	! condensation at ice saturation for temp < -38 degC
	713	! condensation at ice saturation for temp > -38 degC
	714	gammasat(i) = 1.
	715	dgammasatdt(i) = 0.
[3999]	716
	717	ENDIF
	718
[5204]	719	! Note that the delta_hetfreez parameter allows to linearly decrease the
	720	! condensation threshold between the calculated threshold and the ice saturation
	721	! for delta_hetfreez = 1, the threshold is the calculated condensation threshold
	722	! for delta_hetfreez = 0, the threshold is the ice saturation
	723	gammasat(i) = ( 1. - delta_hetfreez ) + delta_hetfreez * gammasat(i)
	724	dgammasatdt(i) = delta_hetfreez * dgammasatdt(i)
	725
[3999]	726	ENDIF
	727
[4072]	728	END DO
	729
	730
[3999]	731	END SUBROUTINE CALC_GAMMASAT
[4562]	732	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
[3999]	733
	734
[4562]	735	!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
[4639]	736	SUBROUTINE DISTANCE_TO_CLOUD_TOP(klon,klev,k,temp,pplay,paprs,rneb,distcltop1D,temp_cltop)
[4562]	737	!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
	738
[4665]	739	USE lmdz_lscp_ini, ONLY : rd,rg,tresh_cl
[4562]	740
	741	IMPLICIT NONE
	742
	743	INTEGER, INTENT(IN) :: klon,klev !number of horizontal and vertical grid points
	744	INTEGER, INTENT(IN) :: k ! vertical index
	745	REAL, INTENT(IN), DIMENSION(klon,klev) :: temp ! temperature in K
	746	REAL, INTENT(IN), DIMENSION(klon,klev) :: pplay ! pressure middle layer in Pa
	747	REAL, INTENT(IN), DIMENSION(klon,klev+1) :: paprs ! pressure interfaces in Pa
	748	REAL, INTENT(IN), DIMENSION(klon,klev) :: rneb ! cloud fraction
	749
	750	REAL, INTENT(OUT), DIMENSION(klon) :: distcltop1D ! distance from cloud top
[4639]	751	REAL, INTENT(OUT), DIMENSION(klon) :: temp_cltop ! temperature of cloud top
	752
[4562]	753	REAL dzlay(klon,klev)
	754	REAL zlay(klon,klev)
	755	REAL dzinterf
	756	INTEGER i,k_top, kvert
	757	LOGICAL bool_cl
	758
	759
	760	DO i=1,klon
	761	! Initialization height middle of first layer
	762	dzlay(i,1) = Rd * temp(i,1) / rg * log(paprs(i,1)/paprs(i,2))
	763	zlay(i,1) = dzlay(i,1)/2
	764
	765	DO kvert=2,klev
	766	IF (kvert.EQ.klev) THEN
	767	dzlay(i,kvert) = 2(rd temp(i,kvert) / rg * log(paprs(i,kvert)/pplay(i,kvert)))
	768	ELSE
	769	dzlay(i,kvert) = rd * temp(i,kvert) / rg * log(paprs(i,kvert)/paprs(i,kvert+1))
	770	ENDIF
	771	dzinterf = rd * temp(i,kvert) / rg * log(pplay(i,kvert-1)/pplay(i,kvert))
	772	zlay(i,kvert) = zlay(i,kvert-1) + dzinterf
	773	ENDDO
	774	ENDDO
	775
	776	DO i=1,klon
[4639]	777	k_top = k
[4562]	778	IF (rneb(i,k) .LE. tresh_cl) THEN
	779	bool_cl = .FALSE.
	780	ELSE
	781	bool_cl = .TRUE.
	782	ENDIF
	783
	784	DO WHILE ((bool_cl) .AND. (k_top .LE. klev))
	785	! find cloud top
	786	IF (rneb(i,k_top) .GT. tresh_cl) THEN
	787	k_top = k_top + 1
	788	ELSE
	789	bool_cl = .FALSE.
	790	k_top = k_top - 1
	791	ENDIF
	792	ENDDO
	793	k_top=min(k_top,klev)
	794
	795	!dist to top is dist between current layer and layer of cloud top (from middle to middle) + dist middle to
	796	!interf for layer of cloud top
	797	distcltop1D(i) = zlay(i,k_top) - zlay(i,k) + dzlay(i,k_top)/2
[4639]	798	temp_cltop(i) = temp(i,k_top)
[4562]	799	ENDDO ! klon
	800
	801	END SUBROUTINE DISTANCE_TO_CLOUD_TOP
[3999]	802	!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
	803
[5204]	804	!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
	805	FUNCTION GAMMAINC ( p, x )
	806
	807	!*****************************************************************************80
	808	!
	809	!! GAMMAINC computes the regularized lower incomplete Gamma Integral
	810	!
	811	! Modified:
	812	!
	813	! 20 January 2008
	814	!
	815	! Author:
	816	!
	817	! Original FORTRAN77 version by B Shea.
	818	! FORTRAN90 version by John Burkardt.
	819	!
	820	! Reference:
	821	!
	822	! B Shea,
	823	! Algorithm AS 239:
	824	! Chi-squared and Incomplete Gamma Integral,
	825	! Applied Statistics,
	826	! Volume 37, Number 3, 1988, pages 466-473.
	827	!
	828	! Parameters:
	829	!
	830	! Input, real X, P, the parameters of the incomplete
	831	! gamma ratio. 0 <= X, and 0 < P.
	832	!
	833	! Output, real GAMMAINC, the value of the incomplete
	834	! Gamma integral.
	835	!
	836	IMPLICIT NONE
	837
	838	REAL A
	839	REAL AN
	840	REAL ARG
	841	REAL B
	842	REAL C
	843	REAL, PARAMETER :: ELIMIT = - 88.0E+00
	844	REAL GAMMAINC
	845	REAL, PARAMETER :: OFLO = 1.0E+37
	846	REAL P
	847	REAL, PARAMETER :: PLIMIT = 1000.0E+00
	848	REAL PN1
	849	REAL PN2
	850	REAL PN3
	851	REAL PN4
	852	REAL PN5
	853	REAL PN6
	854	REAL RN
	855	REAL, PARAMETER :: TOL = 1.0E-14
	856	REAL X
	857	REAL, PARAMETER :: XBIG = 1.0E+08
	858
	859	GAMMAINC = 0.0E+00
	860
	861	IF ( X == 0.0E+00 ) THEN
	862	GAMMAINC = 0.0E+00
	863	RETURN
	864	END IF
	865	!
	866	! IF P IS LARGE, USE A NORMAL APPROXIMATION.
	867	!
	868	IF ( PLIMIT < P ) THEN
	869
	870	PN1 = 3.0E+00 * SQRT ( P ) * ( ( X / P )**( 1.0E+00 / 3.0E+00 ) &
	871	+ 1.0E+00 / ( 9.0E+00 * P ) - 1.0E+00 )
	872
	873	GAMMAINC = 0.5E+00 * ( 1. + ERF ( PN1 ) )
	874	RETURN
	875
	876	END IF
	877	!
	878	! IF X IS LARGE SET GAMMAD = 1.
	879	!
	880	IF ( XBIG < X ) THEN
	881	GAMMAINC = 1.0E+00
	882	RETURN
	883	END IF
	884	!
	885	! USE PEARSON'S SERIES EXPANSION.
	886	! (NOTE THAT P IS NOT LARGE ENOUGH TO FORCE OVERFLOW IN ALOGAM).
	887	!
	888	IF ( X <= 1.0E+00 .OR. X < P ) THEN
	889
	890	ARG = P * LOG ( X ) - X - LOG_GAMMA ( P + 1.0E+00 )
	891	C = 1.0E+00
	892	GAMMAINC = 1.0E+00
	893	A = P
	894
	895	DO
	896
	897	A = A + 1.0E+00
	898	C = C * X / A
	899	GAMMAINC = GAMMAINC + C
	900
	901	IF ( C <= TOL ) THEN
	902	EXIT
	903	END IF
	904
	905	END DO
	906
	907	ARG = ARG + LOG ( GAMMAINC )
	908
	909	IF ( ELIMIT <= ARG ) THEN
	910	GAMMAINC = EXP ( ARG )
	911	ELSE
	912	GAMMAINC = 0.0E+00
	913	END IF
	914	!
	915	! USE A CONTINUED FRACTION EXPANSION.
	916	!
	917	ELSE
	918
	919	ARG = P * LOG ( X ) - X - LOG_GAMMA ( P )
	920	A = 1.0E+00 - P
	921	B = A + X + 1.0E+00
	922	C = 0.0E+00
	923	PN1 = 1.0E+00
	924	PN2 = X
	925	PN3 = X + 1.0E+00
	926	PN4 = X * B
	927	GAMMAINC = PN3 / PN4
	928
	929	DO
	930
	931	A = A + 1.0E+00
	932	B = B + 2.0E+00
	933	C = C + 1.0E+00
	934	AN = A * C
	935	PN5 = B * PN3 - AN * PN1
	936	PN6 = B * PN4 - AN * PN2
	937
	938	IF ( PN6 /= 0.0E+00 ) THEN
	939
	940	RN = PN5 / PN6
	941
	942	IF ( ABS ( GAMMAINC - RN ) <= MIN ( TOL, TOL * RN ) ) THEN
	943	EXIT
	944	END IF
	945
	946	GAMMAINC = RN
	947
	948	END IF
	949
	950	PN1 = PN3
	951	PN2 = PN4
	952	PN3 = PN5
	953	PN4 = PN6
	954	!
	955	! RE-SCALE TERMS IN CONTINUED FRACTION IF TERMS ARE LARGE.
	956	!
	957	IF ( OFLO <= ABS ( PN5 ) ) THEN
	958	PN1 = PN1 / OFLO
	959	PN2 = PN2 / OFLO
	960	PN3 = PN3 / OFLO
	961	PN4 = PN4 / OFLO
	962	END IF
	963
	964	END DO
	965
	966	ARG = ARG + LOG ( GAMMAINC )
	967
	968	IF ( ELIMIT <= ARG ) THEN
	969	GAMMAINC = 1.0E+00 - EXP ( ARG )
	970	ELSE
	971	GAMMAINC = 1.0E+00
	972	END IF
	973
	974	END IF
	975
	976	RETURN
	977	END FUNCTION GAMMAINC
	978	!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
	979
[4665]	980	END MODULE lmdz_lscp_tools
[3999]	981
	982

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

Download in other formats:

Original Format