Context Navigation

co2condens_mod.F @ 2997

Last change on this file since 2997 was 2977, checked in by llange, 19 months ago
Mars PCM Change the latent heat when computing the amount of CO2 that condenses on a slope: a warm air must be cooled off before condensing on the slope. Adapted from pluton PCM (see Eq.12 from Forget et al., 2017) LL
File size: 41.9 KB

Rev	Line
[2009]	1	MODULE co2condens_mod
	2
	3	IMPLICIT NONE
	4
[2184]	5	logical, save :: scavco2cond = .false. ! flag for using scavenging_by_co2
[2616]	6	!$OMP THREADPRIVATE(scavco2cond)
[2184]	7
[2009]	8	CONTAINS
	9
[2953]	10	SUBROUTINE co2condens(ngrid,nlayer,nq,nslope,ptimestep,
[38]	11	$ pcapcal,pplay,pplev,ptsrf,pt,
	12	$ pphi,pdt,pdu,pdv,pdtsrf,pu,pv,pq,pdq,
[2184]	13	$ piceco2,psolaralb,pemisurf,rdust,
[38]	14	$ pdtc,pdtsrfc,pdpsrf,pduc,pdvc,pdqc,
[1996]	15	$ fluxsurf_sw,zls,
	16	$ zdqssed_co2,pcondicea_co2microp,
[2524]	17	$ pdqsc)
[38]	18
[2332]	19	use tracer_mod, only: noms, igcm_h2o_ice, igcm_h2o_vap,
[2184]	20	& igcm_dust_mass, igcm_dust_number,
[2322]	21	& igcm_ccn_mass, igcm_ccn_number,
[2332]	22	& igcm_hdo_ice, igcm_hdo_vap,
	23	& nqperes,nqfils, ! MVals: variables isotopes
[2566]	24	& qperemin,masseqmin,
	25	& igcm_co2
[2739]	26	use surfdat_h, only: emissiv
[2124]	27	use geometry_mod, only: latitude, ! grid point latitudes (rad)
	28	& longitude_deg, latitude_deg
[2009]	29	use planete_h, only: obliquit
	30	use comcstfi_h, only: cpp, g, r, pi
[2409]	31	use dust_param_mod, only: freedust
[2932]	32	use write_output_mod, only: write_output
[1432]	33	#ifndef MESOSCALE
[2124]	34	USE vertical_layers_mod, ONLY: ap, bp
[1432]	35	#endif
[2953]	36	use comslope_mod, ONLY: subslope_dist,def_slope_mean
[38]	37	IMPLICIT NONE
	38	c=======================================================================
	39	c subject:
	40	c --------
	41	c Condensation/sublimation of CO2 ice on the ground and in the
	42	c atmosphere
	43	c (Scheme described in Forget et al., Icarus, 1998)
	44	c
[2009]	45	c author: Francois Forget 1994-1996 ; updated 1996 -- 2018
[38]	46	c ------
[1996]	47	c adapted to external CO2 ice clouds scheme by Deborah Bardet (2018) '
[38]	48	c
	49	c=======================================================================
	50	c
	51	c 0. Declarations :
	52	c ------------------
	53	c
[1528]	54	include "callkeys.h"
[38]	55
	56	c-----------------------------------------------------------------------
	57	c Arguments :
	58	c ---------
[890]	59	INTEGER,INTENT(IN) :: ngrid ! number of atmospheric columns
	60	INTEGER,INTENT(IN) :: nlayer ! number of vertical layers
	61	INTEGER,INTENT(IN) :: nq ! number of tracers
[2953]	62	INTEGER,INTENT(IN) :: nslope ! number of subslope
[38]	63
[890]	64	REAL,INTENT(IN) :: ptimestep ! physics timestep (s)
[2953]	65	REAL,INTENT(IN) :: pcapcal(ngrid,nslope)
[890]	66	REAL,INTENT(IN) :: pplay(ngrid,nlayer) !mid-layer pressure (Pa)
	67	REAL,INTENT(IN) :: pplev(ngrid,nlayer+1) ! inter-layer pressure (Pa)
[2953]	68	REAL,INTENT(IN) :: ptsrf(ngrid,nslope) ! surface temperature (K)
[890]	69	REAL,INTENT(IN) :: pt(ngrid,nlayer) ! atmospheric temperature (K)
	70	REAL,INTENT(IN) :: pphi(ngrid,nlayer) ! geopotential (m2.s-2)
	71	REAL,INTENT(IN) :: pdt(ngrid,nlayer) ! tendency on temperature from
	72	! previous physical processes (K/s)
	73	REAL,INTENT(IN) :: pdu(ngrid,nlayer) ! tendency on zonal wind (m/s2)
	74	! from previous physical processes
	75	REAL,INTENT(IN) :: pdv(ngrid,nlayer) ! tendency on meridional wind (m/s2)
	76	! from previous physical processes
[2953]	77	REAL,INTENT(IN) :: pdtsrf(ngrid,nslope) ! tendency on surface temperature from
[890]	78	! previous physical processes (K/s)
	79	REAL,INTENT(IN) :: pu(ngrid,nlayer) ! zonal wind (m/s)
	80	REAL,INTENT(IN) :: pv(ngrid,nlayer) ! meridional wind (m/s)
	81	REAL,INTENT(IN) :: pq(ngrid,nlayer,nq) ! tracers (../kg_air)
	82	REAL,INTENT(IN) :: pdq(ngrid,nlayer,nq) ! tendency on tracers from
	83	! previous physical processes
[1996]	84
	85	REAL,INTENT(IN) :: zdqssed_co2(ngrid) ! CO2 flux at the surface (kg.m-2.s-1)
	86	REAL,INTENT(IN) :: pcondicea_co2microp(ngrid,nlayer)! tendency due to CO2 condensation (kg/kg.s-1)
	87
[2953]	88	REAL,INTENT(INOUT) :: piceco2(ngrid,nslope) ! CO2 ice on the surface (kg.m-2)
	89	REAL,INTENT(INOUT) :: psolaralb(ngrid,2,nslope) ! albedo of the surface
	90	REAL,INTENT(INOUT) :: pemisurf(ngrid,nslope) ! emissivity of the surface
[2184]	91	REAL,INTENT(IN) :: rdust(ngrid,nlayer) ! dust effective radius
[890]	92
	93	! tendencies due to CO2 condensation/sublimation:
	94	REAL,INTENT(OUT) :: pdtc(ngrid,nlayer) ! tendency on temperature (K/s)
[2953]	95	REAL,INTENT(OUT) :: pdtsrfc(ngrid,nslope) ! tendency on surface temperature (K/s)
[890]	96	REAL,INTENT(OUT) :: pdpsrf(ngrid) ! tendency on surface pressure (Pa/s)
	97	REAL,INTENT(OUT) :: pduc(ngrid,nlayer) ! tendency on zonal wind (m.s-2)
	98	REAL,INTENT(OUT) :: pdvc(ngrid,nlayer) ! tendency on meridional wind (m.s-2)
	99	REAL,INTENT(OUT) :: pdqc(ngrid,nlayer,nq) ! tendency on tracers
[2184]	100	REAL,INTENT(OUT) :: pdqsc(ngrid,nq) ! tendency on surface tracers
[890]	101
	102	! added to calculate flux dependent albedo:
	103	REAL,intent(in) :: fluxsurf_sw(ngrid,2)
	104	real,intent(in) :: zls ! solar longitude (rad)
[38]	105
	106	c
	107	c Local variables :
	108	c -----------------
	109
	110	INTEGER i,j
[2009]	111	INTEGER l,ig,iq,icap
[890]	112	REAL zt(ngrid,nlayer)
[38]	113	REAL zcpi
[1114]	114	REAL ztcond (ngrid,nlayer+1) ! CO2 condensation temperature (atm)
	115	REAL ztcondsol(ngrid) ! CO2 condensation temperature (surface)
[2953]	116	REAL zdiceco2(ngrid,nslope)
	117	REAL zdiceco2_mesh_avg(ngrid)
[2009]	118	REAL zcondicea(ngrid,nlayer) ! condensation rate in layer l (kg/m2/s)
[2953]	119	REAL zcondices(ngrid,nslope) ! condensation rate on the ground (kg/m2/s)
	120	REAL zcondices_mesh_avg(ngrid) ! condensation rate on the ground (kg/m2/s)
[2009]	121	REAL zfallice(ngrid,nlayer+1) ! amount of ice falling from layer l (kg/m2/s)
[2566]	122	REAL condens_layer(ngrid,nlayer) ! co2clouds: condensation rate in layer l (kg/m2/s)
	123	REAL condens_column(ngrid) ! co2clouds: sum(condens_layer(ig,:)) (kg/m2/s)
[2009]	124	REAL zfallheat
[890]	125	REAL zmflux(nlayer+1)
	126	REAL zu(nlayer),zv(nlayer)
[2184]	127	REAL zqc(nlayer,nq),zq1(nlayer)
[2953]	128	REAL ztsrf(ngrid,nslope)
[890]	129	REAL ztc(nlayer), ztm(nlayer+1)
	130	REAL zum(nlayer+1) , zvm(nlayer+1)
	131	REAL zqm(nlayer+1,nq),zqm1(nlayer+1)
	132	REAL masse(nlayer),w(nlayer+1)
	133	REAL Sm(nlayer),Smq(nlayer,nq),mixmas,qmix
[2124]	134	REAL availco2
[2953]	135	LOGICAL condsub(ngrid,nslope)
[38]	136
[2953]	137	real :: emisref(ngrid,nslope)
[2184]	138
[2601]	139	REAL zdq_scav(ngrid,nlayer,nq) ! tendency due to scavenging by co2
[2184]	140	REAL zq(ngrid,nlayer,nq)
[1224]	141
[38]	142	c variable speciale diagnostique
[890]	143	real tconda1(ngrid,nlayer)
	144	real tconda2(ngrid,nlayer)
	145	c REAL zdiceco2a(ngrid) ! for diagnostic only
	146	real zdtsig (ngrid,nlayer)
	147	real zdt (ngrid,nlayer)
	148	real vmr_co2(ngrid,nlayer) ! co2 volume mixing ratio
[38]	149	! improved_ztcond flag: If set to .true. (AND running with a 'co2' tracer)
	150	! then condensation temperature is computed using partial pressure of CO2
	151	logical,parameter :: improved_ztcond=.true.
	152
	153	c local saved variables
[890]	154	integer,save :: ico2 ! index of CO2 tracer
[2124]	155	real,save :: qco2,mmean
[890]	156	real,parameter :: latcond=5.9e5 ! (J/kg) Latent heat of solid CO2 ice
	157	real,parameter :: tcond1mb=136.27 ! condensation temperature (K) at 1 mbar
	158	real,parameter :: cpice=1000. ! (J.kg-1.K-1) specific heat of CO2 ice
	159	REAL,SAVE :: acond,bcond,ccond
	160	real,save :: m_co2, m_noco2, A , B
[38]	161
[890]	162	LOGICAL,SAVE :: firstcall = .true. !,firstcall2=.true.
[2616]	163
	164	!$OMP THREADPRIVATE(ico2,qco2,mmean,acond,bcond,ccond,m_co2,m_noco2)
	165	!$OMP THREADPRIVATE(A,B,firstcall)
[38]	166
[2124]	167	c D.BARDET: for debug
[1996]	168	real ztc3D(ngrid,nlayer)
	169	REAL ztm3D(ngrid,nlayer)
	170	REAL zmflux3D(ngrid,nlayer)
[2322]	171
	172	c MVals: variables isotopes
	173	REAL Ratio1(nlayer),Ratiom1(nlayer+1)
	174	REAL masseq(nlayer),wq(nlayer+1)
	175	INTEGER ifils,iq2
[2953]	176
	177	c Subslope:
	178
	179	REAL :: emisref_tmp(ngrid)
	180	REAL :: alb_tmp(ngrid,2) ! local
	181	REAL :: zcondices_tmp(ngrid) ! local
	182	REAL :: piceco2_tmp(ngrid) ! local
	183	REAL :: pemisurf_tmp(ngrid)! local
	184	LOGICAL :: condsub_tmp(ngrid) !local
	185	REAL :: zfallice_tmp(ngrid,nlayer+1)
	186	REAL :: condens_layer_tmp(ngrid,nlayer) ! co2clouds: condensation rate in layer l (kg/m2/s)
	187	INTEGER :: islope
[38]	188	c----------------------------------------------------------------------
	189
	190	c Initialisation
	191	c --------------
	192	c
[1779]	193	! AS: firstcall OK absolute
[38]	194	IF (firstcall) THEN
[890]	195
[38]	196	bcond=1./tcond1mb
	197	ccond=cpp/(g*latcond)
	198	acond=r/latcond
	199
	200	firstcall=.false.
[2124]	201	write(,) 'CO2condens: improved_ztcond=',improved_ztcond
	202	PRINT*,'In co2condens:Tcond(P=1mb)=',tcond1mb,' Lcond=',latcond
[38]	203	PRINT*,'acond,bcond,ccond',acond,bcond,ccond
	204
	205	ico2=0
	206
[2823]	207	c Prepare Special treatment if one of the tracer is CO2 gas
	208	do iq=1,nq
[38]	209	if (noms(iq).eq."co2") then
	210	ico2=iq
	211	m_co2 = 44.01E-3 ! CO2 molecular mass (kg/mol)
	212	m_noco2 = 33.37E-3 ! Non condensible mol mass (kg/mol)
	213	c Compute A and B coefficient use to compute
	214	c mean molecular mass Mair defined by
	215	c 1/Mair = q(ico2)/m_co2 + (1-q(ico2))/m_noco2
	216	c 1/Mair = A*q(ico2) + B
	217	A =(1/m_co2 - 1/m_noco2)
	218	B=1/m_noco2
	219	endif
[2823]	220	enddo
[890]	221	ENDIF ! of IF (firstcall)
[38]	222	zcpi=1./cpp
[1130]	223
[38]	224	c
	225	c======================================================================
	226	c Calcul of CO2 condensation sublimation
	227	c ============================================================
	228	c
	229	c Used variable :
	230	c piceco2(ngrid) : amount of co2 ice on the ground (kg/m2)
	231	c zcondicea(ngrid,l): condensation rate in layer l (kg/m2/s)
	232	c zcondices(ngrid): condensation rate on the ground (kg/m2/s)
	233	c zfallice(ngrid,l):amount of ice falling from layer l (kg/m2/s)
	234	c
[1047]	235	c pdtc(ngrid,nlayer) : dT/dt due to cond/sub
[38]	236	c
	237	c
[2184]	238	c Tendencies set to 0
[38]	239	c -------------------------------------
[2189]	240	zcondicea(1:ngrid,1:nlayer) = 0.
	241	zfallice(1:ngrid,1:nlayer+1) = 0.
	242	pduc(1:ngrid,1:nlayer) = 0
	243	pdvc(1:ngrid,1:nlayer) = 0
	244	pdtc(1:ngrid,1:nlayer) = 0.
	245	pdqsc(1:ngrid,1:nq) = 0
[38]	246
[2189]	247	pdqc(1:ngrid,1:nlayer,1:nq) = 0
[38]	248
[2953]	249	zcondices(1:ngrid,1:nslope) = 0.
	250	zcondices_mesh_avg(1:ngrid)=0.
	251	pdtsrfc(1:ngrid,1:nslope) = 0.
[2189]	252	pdpsrf(1:ngrid) = 0.
[2953]	253	condsub(1:ngrid,1:nslope) = .false.
	254	zdiceco2(1:ngrid,1:nslope) = 0.
	255	zdiceco2_mesh_avg(1:ngrid)=0.
[2189]	256
[38]	257	zfallheat=0
[2184]	258
	259	zdq_scav(:,:,:)=0.
[38]	260
[2184]	261	c Update tendencies from previous processes
	262	c -------------------------------------
	263	DO l=1,nlayer
	264	DO ig=1,ngrid
	265	zt(ig,l)=pt(ig,l)+ pdt(ig,l)*ptimestep
	266	do iq=1,nq
	267	zq(ig,l,iq)=pq(ig,l,iq)+pdq(ig,l,iq)*ptimestep
	268	enddo
	269	ENDDO
	270	ENDDO
	271
[38]	272	c *************************
	273	c ATMOSPHERIC CONDENSATION
	274	c *************************
	275
	276	c Compute CO2 Volume mixing ratio
	277	c -------------------------------
	278	if (improved_ztcond.and.(ico2.ne.0)) then
	279	DO l=1,nlayer
	280	DO ig=1,ngrid
	281	qco2=pq(ig,l,ico2)+pdq(ig,l,ico2)*ptimestep
	282	c Mean air molecular mass = 1/(q(ico2)/m_co2 + (1-q(ico2))/m_noco2)
	283	mmean=1/(A*qco2 +B)
	284	vmr_co2(ig,l) = qco2*mmean/m_co2
	285	ENDDO
	286	ENDDO
	287	else
	288	DO l=1,nlayer
	289	DO ig=1,ngrid
	290	vmr_co2(ig,l)=0.95
	291	ENDDO
	292	ENDDO
[2184]	293	endif
[38]	294
[1996]	295	IF (.NOT. co2clouds) then
[38]	296	c forecast of atmospheric temperature zt and frost temperature ztcond
	297	c --------------------------------------------------------------------
	298
	299	DO l=1,nlayer
	300	DO ig=1,ngrid
	301	! ztcond(ig,l)=1./(bcond-acondlog(.0095pplay(ig,l)))
[1263]	302	if (pplay(ig,l).ge.1e-4) then
	303	ztcond(ig,l)=
[38]	304	& 1./(bcond-acondlog(.01vmr_co2(ig,l)*pplay(ig,l)))
[1263]	305	else
	306	ztcond(ig,l)=0.0 !mars Monica
	307	endif
[38]	308	ENDDO
	309	ENDDO
[327]	310
[328]	311	ztcond(:,nlayer+1)=ztcond(:,nlayer)
[38]	312
	313	c Condensation/sublimation in the atmosphere
	314	c ------------------------------------------
	315	c (calcul of zcondicea , zfallice and pdtc)
	316	c
	317	DO l=nlayer , 1, -1
	318	DO ig=1,ngrid
	319	pdtc(ig,l)=0.
	320	IF((zt(ig,l).LT.ztcond(ig,l)).or.(zfallice(ig,l+1).gt.0))THEN
[2953]	321	condsub(ig,:)=.true.
[38]	322	IF (zfallice(ig,l+1).gt.0) then
	323	zfallheat=zfallice(ig,l+1)*
	324	& (pphi(ig,l+1)-pphi(ig,l) +
	325	& cpice*(ztcond(ig,l+1)-ztcond(ig,l)))/latcond
	326	ELSE
	327	zfallheat=0.
	328	ENDIF
	329	pdtc(ig,l)=(ztcond(ig,l) - zt(ig,l))/ptimestep
	330	zcondicea(ig,l)=(pplev(ig,l)-pplev(ig,l+1))
	331	& ccondpdtc(ig,l)- zfallheat
	332	c Case when the ice from above sublimes entirely
	333	c """""""""""""""""""""""""""""""""""""""""""""""
	334	IF (zfallice(ig,l+1).lt.- zcondicea(ig,l)) then
	335	pdtc(ig,l)=(-zfallice(ig,l+1)+zfallheat)/
	336	& (ccond*(pplev(ig,l)-pplev(ig,l+1)))
	337	zcondicea(ig,l)= -zfallice(ig,l+1)
	338	END IF
	339
	340	zfallice(ig,l) = zcondicea(ig,l)+zfallice(ig,l+1)
	341	END IF
	342	ENDDO
	343	ENDDO
[1996]	344
[2599]	345	condens_layer(:,:) = zcondicea(:,:)
	346	condens_column(:) = 0.
	347
[2184]	348	if (scavco2cond) then
	349	call scavenging_by_co2(ngrid,nlayer,nq,ptimestep,pplev,zq,
	350	& rdust,zcondicea,zfallice,zdq_scav,pdqsc)
	351	endif
[2934]	352	call write_output("co2condens_zfallice",
	353	& "CO2 ice tendency on the surface",
	354	& "kg.m-2.s-1",zfallice(:,1))
[2184]	355	ELSE ! if co2 clouds
[2566]	356	condens_layer(:,:) = 0.
	357	condens_column(:) = 0.
[1996]	358	DO l=nlayer , 1, -1
	359	DO ig=1,ngrid
[2566]	360	condens_layer(ig,l) = pcondicea_co2microp(ig,l)*
	361	& (pplev(ig,l) - pplev(ig,l+1))/g
[1996]	362	ENDDO
	363	ENDDO
[2566]	364	DO ig=1,ngrid
	365	condens_column(ig) = sum(condens_layer(ig,:))
[2599]	366	zfallice(ig,1) = zdqssed_co2(ig)
[2953]	367	DO islope = 1,nslope
	368	piceco2(ig,islope) = piceco2(ig,islope) +
	369	& zdqssed_co2(ig)ptimestep
	370	& cos(pi*def_slope_mean(islope)/180.)
	371	ENDDO
[2566]	372	ENDDO
[2934]	373	call write_output("co2condens_zfallice",
	374	& "CO2 ice tendency on the surface",
	375	& "kg.m-2.s-1",zdqssed_co2(:)) ! otherwise we have not 1 day(1proc) = 1 day (24procs) test
[2184]	376	ENDIF ! end of if co2clouds
[38]	377
[2934]	378	call write_output("co2condens_pdtc",
	379	& "Temperature tendency due to CO2 condensation",
	380	& "K.s-1",pdtc(:,:))
[1996]	381
[2934]	382	! call write_output("condens_layer",
	383	! & "",
	384	! & " ",condens_layer(:,:))
[1996]	385
[38]	386	c *************************
	387	c SURFACE CONDENSATION
	388	c *************************
	389
	390	c forecast of ground temperature ztsrf and frost temperature ztcondsol
	391	c --------------------------------------------------------------------
	392	DO ig=1,ngrid
	393	ztcondsol(ig)=
	394	& 1./(bcond-acondlog(.01vmr_co2(ig,1)*pplev(ig,1)))
[2953]	395	DO islope = 1,nslope
	396	ztsrf(ig,islope) = ptsrf(ig,islope) +
	397	& pdtsrf(ig,islope)*ptimestep
	398	ENDDO
[38]	399	ENDDO
	400
	401	c
	402	c Condensation/sublimation on the ground
	403	c --------------------------------------
[1114]	404	c (compute zcondices and pdtsrfc)
[38]	405	c
[2566]	406	c No microphysics of CO2 clouds
[38]	407	DO ig=1,ngrid
[1541]	408	IF(latitude(ig).lt.0) THEN
[1114]	409	! Southern hemisphere
[38]	410	icap=2
	411	ELSE
[1114]	412	! Northern hemisphere
[38]	413	icap=1
	414	ENDIF
[2953]	415
	416	DO islope = 1,nslope
	417	c Need first to put piceco2_slope(ig,islope) in kg/m^2 flat
	418
	419	piceco2(ig,islope) = piceco2(ig,islope)
	420	& /cos(pi*def_slope_mean(islope)/180.)
	421
[38]	422	c
	423	c Loop on where we have condensation/ sublimation
[2953]	424	IF ((ztsrf(ig,islope) .LT. ztcondsol(ig)) .OR. ! ground cond
[38]	425	$ (zfallice(ig,1).NE.0.) .OR. ! falling snow
[2953]	426	$ ((ztsrf(ig,islope) .GT. ztcondsol(ig)) .AND. ! ground sublim.
	427	$ ((piceco2(ig,islope)+zfallice(ig,1)*ptimestep)
	428	$ .NE. 0.))) THEN
	429	condsub(ig,islope) = .true.
[38]	430
[2599]	431	IF (zfallice(ig,1).gt.0 .and. .not. co2clouds) then
[2739]	432	! potential eneregy release due to the impact of the snowfall
	433	! appendix of forget et al. 1999
	434	zfallheat = zfallice(ig,1) * (pphi(ig,1) +
[2599]	435	& cpice*(ztcond(ig,1)-ztcondsol(ig)))/latcond
[38]	436	ELSE
[2599]	437	zfallheat = 0.
[38]	438	ENDIF
	439	c condensation or partial sublimation of CO2 ice
	440	c """""""""""""""""""""""""""""""""""""""""""""""
[2953]	441	if(ztsrf(ig,islope).LT. ztcondsol(ig)) then
	442	c condensation
	443	zcondices(ig,islope)=pcapcal(ig,islope)
	444	& *(ztcondsol(ig)-ztsrf(ig,islope))
[2977]	445	& /((latcond+cpp(zt(ig,1)-ztcondsol(ig)))ptimestep)
[2953]	446	& - zfallheat
	447	else
	448	c sublimation
	449	zcondices(ig,islope)=pcapcal(ig,islope)
	450	& *(ztcondsol(ig)-ztsrf(ig,islope))
	451	& /(latcond*ptimestep)
	452	& - zfallheat
	453	endif
	454	pdtsrfc(ig,islope) = (ztcondsol(ig) - ztsrf(ig,islope))
	455	& /ptimestep
[2153]	456	#ifdef MESOSCALE
	457	print*, "not enough CO2 tracer in 1st layer to condense"
	458	print*, ">>> to be implemented in the mesoscale case"
	459	print*, "because this uses ap levels..."
	460	#else
[2124]	461	c If there is not enough CO2 tracer in 1st layer to condense
	462	c """"""""""""""""""""""""""""""""""""""""""""""""""""""
	463	IF(ico2.ne.0) then
[2599]	464	c Available CO2 tracer in layer 1 at end of timestep (kg/m2)
[2566]	465	#ifndef MESOSCALE
[2599]	466	availco2 = pq(ig,1,ico2)*((ap(1)-ap(2))+
	467	& (bp(1)-bp(2))*(pplev(ig,1)/g -
[2953]	468	& (zcondices(ig,islope) + zfallice(ig,1))
	469	& *ptimestep))
	470	if ((zcondices(ig,islope) + condens_layer(ig,1))
	471	& *ptimestep
[2599]	472	& .gt.availco2) then
[2953]	473	zcondices(ig,islope) = availco2/ptimestep -
[2599]	474	& condens_layer(ig,1)
[2953]	475	pdtsrfc(ig,islope) = (latcond/pcapcal(ig,islope))*
	476	& (zcondices(ig,islope)+zfallheat)
[2599]	477	end if
[2566]	478	#else
	479	availco2 = pq(ig,1,igcm_co2)
	480	PRINT*, "MESOSCALE: CO2 tracer AT FIRST LEVEL IS NOT'
	481	& ' CORRECTED FROM SIGMA LEVELS"
	482	#endif
	483	ENDIF
[2153]	484	#endif
[2124]	485
[2953]	486	c If the entire CO2 ice layer sublimes on the slope
[38]	487	c """"""""""""""""""""""""""""""""""""""""""""""""""""
	488	c (including what has just condensed in the atmosphere)
[2660]	489	if (co2clouds) then
[2953]	490	IF((piceco2(ig,islope)/ptimestep).LE.
	491	& -zcondices(ig,islope))THEN
	492	zcondices(ig,islope) = -piceco2(ig,islope)/ptimestep
	493	pdtsrfc(ig,islope)=(latcond/pcapcal(ig,islope)) *
	494	& (zcondices(ig,islope)+zfallheat)
[2660]	495	END IF
	496	else
[2953]	497	IF((piceco2(ig,islope)/ptimestep+zfallice(ig,1)).LE.
	498	& -zcondices(ig,islope))THEN
	499	zcondices(ig,islope) = -piceco2(ig,islope)
	500	& /ptimestep - zfallice(ig,1)
	501	pdtsrfc(ig,islope)=(latcond/pcapcal(ig,islope)) *
	502	& (zcondices(ig,islope)+zfallheat)
[38]	503	END IF
[2660]	504	end if
[38]	505
[2953]	506	c Changing CO2 ice amount and pressure per slope:
[38]	507	c """"""""""""""""""""""""""""""""""""
[2953]	508	zdiceco2(ig,islope) = zcondices(ig,islope)+zfallice(ig,1)
[2599]	509	& + condens_column(ig)
	510	if (co2clouds) then
	511	! add here only direct condensation/sublimation
[2953]	512	piceco2(ig,islope) = piceco2(ig,islope) +
	513	& zcondices(ig,islope)*ptimestep
[2599]	514	else
	515	! add here also CO2 ice in the atmosphere
[2953]	516	piceco2(ig,islope) = piceco2(ig,islope) +
	517	& zdiceco2(ig,islope)*ptimestep
[2599]	518	end if
[38]	519
[2953]	520	zcondices_mesh_avg(ig) = zcondices_mesh_avg(ig) +
	521	& zcondices(ig,islope)* subslope_dist(ig,islope)
	522
	523	zdiceco2_mesh_avg(ig) = zdiceco2_mesh_avg(ig) +
	524	& zdiceco2(ig,islope)* subslope_dist(ig,islope)
	525
	526	END IF ! if there is condensation/sublimation
	527
	528	piceco2(ig,islope) = piceco2(ig,islope)
	529	& * cos(pi*def_slope_mean(islope)/180.)
	530
	531	ENDDO !islope
	532
	533	pdpsrf(ig) = -zdiceco2_mesh_avg(ig)*g
	534
[38]	535	IF(ABS(pdpsrf(ig)*ptimestep).GT.pplev(ig,1)) THEN
	536	PRINT*,'STOP in condens'
	537	PRINT*,'condensing more than total mass'
	538	PRINT*,'Grid point ',ig
[2124]	539	PRINT*,'Longitude(degrees): ',longitude_deg(ig)
	540	PRINT*,'Latitude(degrees): ',latitude_deg(ig)
[38]	541	PRINT*,'Ps = ',pplev(ig,1)
	542	PRINT*,'d Ps = ',pdpsrf(ig)
[2399]	543	CALL abort_physic('co2condens',
	544	& 'condensing more than total mass', 1)
[38]	545	ENDIF
[2953]	546
[38]	547	ENDDO ! of DO ig=1,ngrid
[2953]	548
[38]	549
	550	c ********************************************************************
	551	c Surface albedo and emissivity of the surface below the snow (emisref)
	552	c ********************************************************************
	553
	554	! Check that amont of CO2 ice is not problematic
	555	DO ig=1,ngrid
[2953]	556	DO islope = 1,nslope
	557	if(.not.piceco2(ig,islope).ge.0.) THEN
	558	if(piceco2(ig,islope).le.-5.e-8) print*,
	559	$ 'WARNING co2condens piceco2(',ig,')=', piceco2(ig,islope)
	560	piceco2(ig,islope)=0.
[38]	561	endif
[2953]	562	ENDDO
[38]	563	ENDDO
	564
	565	! Set albedo and emissivity of the surface
	566	! ----------------------------------------
[2953]	567	DO islope = 1,nslope
	568	piceco2_tmp(:) = piceco2(:,islope)
	569	alb_tmp(:,:) = psolaralb(:,:,islope)
	570	emisref_tmp(:) = 0.
	571	CALL albedocaps(zls,ngrid,piceco2_tmp,alb_tmp,emisref_tmp)
	572	psolaralb(:,1,islope) = alb_tmp(:,1)
	573	psolaralb(:,2,islope) = alb_tmp(:,2)
	574	emisref(:,islope) = emisref_tmp(:)
	575	ENDDO
[38]	576
	577	! set pemisurf() to emissiv when there is bare surface (needed for co2snow)
	578	DO ig=1,ngrid
[2953]	579	DO islope = 1,nslope
	580	if (piceco2(ig,islope).eq.0) then
	581	pemisurf(ig,islope)=emissiv
	582	endif
	583	ENDDO
[38]	584	ENDDO
	585
	586	! firstcall2=.false.
	587	c ***************************************************************
	588	c Correction to account for redistribution between sigma or hybrid
	589	c layers when changing surface pressure (and warming/cooling
	590	c of the CO2 currently changing phase).
	591	c *************************************************************
	592
	593	DO ig=1,ngrid
[2953]	594	if (any(condsub(ig,:))) then
[38]	595	do l=1,nlayer
	596	ztc(l) =zt(ig,l) +pdtc(ig,l) *ptimestep
	597	zu(l) =pu(ig,l) +pdu( ig,l) *ptimestep
	598	zv(l) =pv(ig,l) +pdv( ig,l) *ptimestep
[1036]	599	do iq=1,nq
[2601]	600	zqc(l,iq)=zq(ig,l,iq)+zdq_scav(ig,l,iq)*ptimestep ! zdq_scav=0 if co2clouds=true
[38]	601	enddo
[2184]	602	enddo
[38]	603
	604	c Mass fluxes through the sigma levels (kg.m-2.s-1) (>0 when up)
	605	c """"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
[2953]	606	zmflux(1) = -zcondices_mesh_avg(ig) - zdqssed_co2(ig)
[2566]	607	DO l=1,nlayer
[2599]	608	zmflux(l+1) = zmflux(l) - condens_layer(ig,l)
[1438]	609	#ifndef MESOSCALE
[2599]	610	& + (bp(l)-bp(l+1))*(-pdpsrf(ig)/g)
[2566]	611	c zmflux set to 0 if very low to avoid: top layer is disappearing in v1ld
	612	if (abs(zmflux(l+1)).lt.1E-13.OR.bp(l+1).eq.0.) then
	613	zmflux(l+1)=0.
	614	end if
[1438]	615	#else
[2599]	616	zmflux(l+1) = zmflux(l) - condens_layer(ig,l)
[2566]	617	if (abs(zmflux(l+1)).lt.1E-13) zmflux(l+1)=0.
	618	PRINT*, "MESOSCALE: FLUX THROUGH SIGMA LEVELS from"//
	619	& "dPS HAVE TO BE IMPLEMENTED"
[1438]	620	#endif
[2566]	621	END DO
[2155]	622	#ifdef MESOSCALE
	623	print*,"absurd mass set because bp not available"
	624	print*,"TO BE FIXED"
	625	#else
[2124]	626	c Mass of each layer at the end of timestep
	627	c -----------------------------------------
[38]	628	DO l=1,nlayer
[2124]	629	masse(l)=( pplev(ig,l) - pplev(ig,l+1) +
	630	& (bp(l)-bp(l+1))pdpsrf(ig)ptimestep)/g
[38]	631	END DO
[2155]	632	#endif
[38]	633
	634	c Corresponding fluxes for T,U,V,Q
	635	c """"""""""""""""""""""""""""""""
	636
	637	c averaging operator for TRANSPORT
	638	c """"""""""""""""""""""""""""""""
	639	c Value transfert at the surface interface when condensation
	640	c sublimation:
[2903]	641	ztm(1) = ztcondsol(ig)
[38]	642	zum(1) = 0
	643	zvm(1) = 0
[1036]	644	do iq=1,nq
[38]	645	zqm(1,iq)=0. ! most tracer do not condense !
	646	enddo
	647	c Special case if one of the tracer is CO2 gas
	648	if (ico2.ne.0) zqm(1,ico2)=1. ! flux is 100% CO2
	649
	650	c Van Leer scheme:
	651	DO l=1,nlayer+1
	652	w(l)=-zmflux(l)*ptimestep
	653	END DO
[1269]	654	call vl1d(nlayer,ztc,2.,masse,w,ztm)
	655	call vl1d(nlayer,zu ,2.,masse,w,zum)
	656	call vl1d(nlayer,zv ,2.,masse,w,zvm)
[2322]	657	! MVals: loop over the fathers ("peres")
	658	do iq=1,nqperes
[38]	659	do l=1,nlayer
[2184]	660	zq1(l)=zqc(l,iq)
[38]	661	enddo
	662	zqm1(1)=zqm(1,iq)
[1269]	663	call vl1d(nlayer,zq1,2.,masse,w,zqm1)
[38]	664	do l=2,nlayer
[2184]	665	zqc(l,iq)=zq1(l)
[38]	666	zqm(l,iq)=zqm1(l)
	667	enddo
[2322]	668	! MVals: loop over the sons ("fils")
	669	if (nqfils(iq).gt.0) then
[2332]	670	if (iq.eq.igcm_h2o_ice) then
	671	iq2=igcm_hdo_ice
	672	else if (iq.eq.igcm_h2o_vap) then
	673	iq2=igcm_hdo_vap
	674	else
	675	call abort_physic("co2condens_mod","invalid isotope",1)
	676	endif
[2322]	677	do l=1,nlayer
	678	if (zqc(l,iq).gt.qperemin) then
	679	Ratio1(l)=zqc(l,iq2)/zqc(l,iq)
	680	else
	681	Ratio1(l)=0.
	682	endif
	683	masseq(l)=max(masse(l)*zqc(l,iq),masseqmin)
	684	wq(l)=w(l)*zqm(l,iq)
	685	enddo
	686	Ratiom1(1)=zqm(1,iq2)
	687	call vl1d(nlayer,Ratio1,2.,masseq,wq,Ratiom1)
	688	zqm(1,iq2)=Ratiom1(1)*zqc(1,iq)
	689	do l=2,nlayer
	690	zqm(l,iq2)=Ratiom1(l)*zqm(l,iq)
	691	enddo
	692	endif !if (nqfils(iq).gt.0) then
	693	enddo !iq=1,nqperes
[38]	694
	695	c Surface condensation affects low winds
	696	if (zmflux(1).lt.0) then
	697	zum(1)= zu(1) * (w(1)/masse(1))
	698	zvm(1)= zv(1) * (w(1)/masse(1))
	699	if (w(1).gt.masse(1)) then ! ensure numerical stability
	700	zum(1)= (zu(1)-zum(2))*masse(1)/w(1) + zum(2)
	701	zvm(1)= (zv(1)-zvm(2))*masse(1)/w(1) + zvm(2)
	702	end if
	703	end if
	704
	705	ztm(nlayer+1)= ztc(nlayer) ! should not be used, but...
	706	zum(nlayer+1)= zu(nlayer) ! should not be used, but...
	707	zvm(nlayer+1)= zv(nlayer) ! should not be used, but...
[1036]	708	do iq=1,nq
[2184]	709	zqm(nlayer+1,iq)= zqc(nlayer,iq)
[38]	710	enddo
[86]	711
	712	#ifdef MESOSCALE
	713	!!!! AS: This part must be commented in the mesoscale model
	714	!!!! AS: ... to avoid instabilities.
	715	!!!! AS: you have to compile with -DMESOSCALE to do so
	716	#else
[38]	717	c Tendencies on T, U, V, Q
	718	c """"""""""""""""""""""""
	719	DO l=1,nlayer
[1996]	720	IF(.not. co2clouds) THEN
[38]	721	c Tendencies on T
	722	zdtsig(ig,l) = (1/masse(l)) *
	723	& ( zmflux(l)*(ztm(l) - ztc(l))
	724	& - zmflux(l+1)*(ztm(l+1) - ztc(l))
[2599]	725	& + condens_layer(ig,l)*(ztcond(ig,l)-ztc(l)) )
[1996]	726	ELSE
	727	zdtsig(ig,l) = (1/masse(l)) *
	728	& ( zmflux(l)*(ztm(l) - ztc(l))
	729	& - zmflux(l+1)*(ztm(l+1) - ztc(l)))
	730	ENDIF
	731	c D.BARDET: for diagnotics
	732	zmflux3D(ig,l)=zmflux(l)
	733	ztm3D(ig,l)=ztm(l)
	734	ztc3D(ig,l)=ztc(l)
	735
[38]	736	pdtc(ig,l) = pdtc(ig,l) + zdtsig(ig,l)
	737
	738	c Tendencies on U
	739	pduc(ig,l) = (1/masse(l)) *
	740	& ( zmflux(l)*(zum(l) - zu(l))
	741	& - zmflux(l+1)*(zum(l+1) - zu(l)) )
	742
	743
	744	c Tendencies on V
	745	pdvc(ig,l) = (1/masse(l)) *
	746	& ( zmflux(l)*(zvm(l) - zv(l))
	747	& - zmflux(l+1)*(zvm(l+1) - zv(l)) )
	748
	749	END DO
[1114]	750
[86]	751	#endif
[38]	752
[2566]	753	do iq=1,nq
	754	! if (noms(iq).eq.'co2') then
	755	if (iq.eq.ico2) then
	756	c SPECIAL Case when the tracer IS CO2 :
	757	DO l=1,nlayer
	758	pdqc(ig,l,iq)= (1/masse(l)) *
	759	& ( zmflux(l)*(zqm(l,iq) - zqc(l,iq))
	760	& - zmflux(l+1)*(zqm(l+1,iq) - zqc(l,iq))
[2599]	761	& + condens_layer(ig,l)*(zqc(l,iq)-1.) )
[2566]	762	END DO
	763	else
	764	DO l=1,nlayer
	765	pdqc(ig,l,iq)= (1/masse(l)) *
	766	& ( zmflux(l)*(zqm(l,iq) - zqc(l,iq))
	767	& - zmflux(l+1)*(zqm(l+1,iq) - zqc(l,iq))
[2599]	768	& + condens_layer(ig,l)*zqc(l,iq) )
[2566]	769
[2601]	770	pdqc(ig,l,iq)=pdqc(ig,l,iq)+zdq_scav(ig,l,iq) ! zdq_scav=0 if co2clouds=true
[2566]	771	END DO
	772	end if
	773	enddo
[38]	774
	775	end if ! if (condsub)
	776	END DO ! loop on ig
	777
	778	c ***************************************************************
	779	c CO2 snow / clouds scheme
	780	c ***************************************************************
[2953]	781	DO islope = 1,nslope
	782	emisref_tmp(:) = emisref(:,islope)
	783	condsub_tmp(:) = condsub(:,islope)
	784	condens_layer_tmp(:,:) = condens_layer(:,:)*
	785	& cos(pi*def_slope_mean(islope)/180.)
	786	zcondices_tmp(:) = zcondices(:,islope)*
	787	& cos(pi*def_slope_mean(islope)/180.)
	788	zfallice_tmp(:,:) = zfallice(:,:)*
	789	& cos(pi*def_slope_mean(islope)/180.)
	790	pemisurf_tmp(:) = pemisurf(:,islope)
	791
	792	call co2snow(ngrid,nlayer,ptimestep,emisref_tmp,condsub_tmp,
	793	& pplev,condens_layer_tmp,zcondices_tmp,zfallice_tmp,
	794	& pemisurf_tmp)
	795	pemisurf(:,islope) = pemisurf_tmp(:)
	796
	797	ENDDO
[38]	798	c ***************************************************************
	799	c Ecriture des diagnostiques
	800	c ***************************************************************
	801
	802	c DO l=1,nlayer
	803	c DO ig=1,ngrid
	804	c Taux de cond en kg.m-2.pa-1.s-1
	805	c tconda1(ig,l)=zcondicea(ig,l)/(pplev(ig,l)-pplev(ig,l+1))
	806	c Taux de cond en kg.m-3.s-1
	807	c tconda2(ig,l)=tconda1(ig,l)pplay(ig,l)g/(r*pt(ig,l))
	808	c END DO
	809	c END DO
[2932]	810	c call write_output('tconda1',
[38]	811	c &'Taux de condensation CO2 atmospherique /Pa',
[2932]	812	c & 'kg.m-2.Pa-1.s-1',tconda1)
	813	c call write_output('tconda2',
[38]	814	c &'Taux de condensation CO2 atmospherique /m',
[2932]	815	c & 'kg.m-3.s-1',tconda2)
[38]	816
	817	! output falling co2 ice in 1st layer:
[2932]	818	! call write_output('fallice',
[38]	819	! &'Precipitation of co2 ice',
[2932]	820	! & 'kg.m-2.s-1',zfallice(1,1))
[38]	821
[1114]	822	#ifndef MESOSCALE
	823	! Extra special case for surface temperature tendency pdtsrfc:
	824	! we want to fix the south pole temperature to CO2 condensation temperature
	825	if (caps.and.(obliquit.lt.27.)) then
	826	! check if last grid point is the south pole
[1541]	827	if (abs(latitude(ngrid)-(-pi/2.)).lt.1.e-5) then
[1114]	828	! NB: Updated surface pressure, at grid point 'ngrid', is
	829	! ps(ngrid)=pplev(ngrid,1)+pdpsrf(ngrid)*ptimestep
[2903]	830	ztcondsol(ngrid)=
	831	& 1./(bcond-acondlog(.01vmr_co2(ngrid,1)*
	832	& (pplev(ngrid,1)+pdpsrf(ngrid)*ptimestep)))
[2953]	833	DO islope = 1,nslope
	834	pdtsrfc(ngrid,islope)=(ztcondsol(ngrid)-
	835	& ztsrf(ngrid,islope))/ptimestep
	836	ENDDO ! islope = 1,nslope
[1114]	837	endif
	838	endif
	839	#endif
	840
[2009]	841	END SUBROUTINE co2condens
[38]	842
	843
	844
	845	c *****************************************************************
[1269]	846	SUBROUTINE vl1d(nlayer,q,pente_max,masse,w,qm)
[38]	847	c
	848	c
	849	c Operateur de moyenne inter-couche pour calcul de transport type
	850	c Van-Leer " pseudo amont " dans la verticale
	851	c q,w sont des arguments d'entree pour le s-pg ....
	852	c masse : masse de la couche Dp/g
	853	c w : masse d'atm ``transferee'' a chaque pas de temps (kg.m-2)
	854	c pente_max = 2 conseillee
	855	c
	856	c
	857	c --------------------------------------------------------------------
	858	IMPLICIT NONE
	859
	860	c
	861	c
	862	c
	863	c Arguments:
	864	c ----------
[1270]	865	integer nlayer
[1269]	866	real masse(nlayer),pente_max
	867	REAL q(nlayer),qm(nlayer+1)
	868	REAL w(nlayer+1)
[38]	869	c
	870	c Local
	871	c ---------
	872	c
	873	INTEGER l
	874	c
[1269]	875	real dzq(nlayer),dzqw(nlayer),adzqw(nlayer),dzqmax
[38]	876	real sigw, Mtot, MQtot
	877	integer m
	878	c integer ismax,ismin
	879
	880
	881	c On oriente tout dans le sens de la pression
	882	c W > 0 WHEN DOWN !!!!!!!!!!!!!
	883
[1269]	884	do l=2,nlayer
[38]	885	dzqw(l)=q(l-1)-q(l)
	886	adzqw(l)=abs(dzqw(l))
	887	enddo
	888
[1269]	889	do l=2,nlayer-1
[38]	890	if(dzqw(l)*dzqw(l+1).gt.0.) then
	891	dzq(l)=0.5*(dzqw(l)+dzqw(l+1))
	892	else
	893	dzq(l)=0.
	894	endif
	895	dzqmax=pente_max*min(adzqw(l),adzqw(l+1))
	896	dzq(l)=sign(min(abs(dzq(l)),dzqmax),dzq(l))
	897	enddo
	898
	899	dzq(1)=0.
[1269]	900	dzq(nlayer)=0.
[38]	901
[1269]	902	do l = 1,nlayer-1
[38]	903
	904	c Regular scheme (transfered mass < layer mass)
	905	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	906	if(w(l+1).gt.0. .and. w(l+1).le.masse(l+1)) then
	907	sigw=w(l+1)/masse(l+1)
	908	qm(l+1)=(q(l+1)+0.5(1.-sigw)dzq(l+1))
	909	else if(w(l+1).le.0. .and. -w(l+1).le.masse(l)) then
	910	sigw=w(l+1)/masse(l)
	911	qm(l+1)=(q(l)-0.5(1.+sigw)dzq(l))
	912
	913	c Extended scheme (transfered mass > layer mass)
	914	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	915	else if(w(l+1).gt.0.) then
	916	m=l+1
	917	Mtot = masse(m)
	918	MQtot = masse(m)*q(m)
[1269]	919	do while ((m.lt.nlayer).and.(w(l+1).gt.(Mtot+masse(m+1))))
[38]	920	m=m+1
	921	Mtot = Mtot + masse(m)
	922	MQtot = MQtot + masse(m)*q(m)
	923	end do
[1269]	924	if (m.lt.nlayer) then
[38]	925	sigw=(w(l+1)-Mtot)/masse(m+1)
	926	qm(l+1)= (1/w(l+1))(MQtot + (w(l+1)-Mtot)
	927	& (q(m+1)+0.5(1.-sigw)dzq(m+1)) )
	928	else
	929	w(l+1) = Mtot
	930	qm(l+1) = Mqtot / Mtot
[2399]	931	CALL abort_physic('co2condens',
	932	& 'top layer is disapearing !', 1)
[38]	933	end if
	934	else ! if(w(l+1).lt.0)
	935	m = l-1
	936	Mtot = masse(m+1)
	937	MQtot = masse(m+1)*q(m+1)
[120]	938	if (m.gt.0) then ! because some compilers will have problems
	939	! evaluating masse(0)
	940	do while ((m.gt.0).and.(-w(l+1).gt.(Mtot+masse(m))))
[38]	941	m=m-1
	942	Mtot = Mtot + masse(m+1)
	943	MQtot = MQtot + masse(m+1)*q(m+1)
[120]	944	if (m.eq.0) exit
	945	end do
	946	endif
[38]	947	if (m.gt.0) then
	948	sigw=(w(l+1)+Mtot)/masse(m)
	949	qm(l+1)= (-1/w(l+1))(MQtot + (-w(l+1)-Mtot)
	950	& (q(m)-0.5(1.+sigw)dzq(m)) )
	951	else
	952	qm(l+1)= (-1/w(l+1))(MQtot + (-w(l+1)-Mtot)qm(1))
	953	end if
	954	end if
	955	enddo
	956
[2124]	957	c boundary conditions (not used in co2condens !!)
[1269]	958	c qm(nlayer+1)=0.
[38]	959	c if(w(1).gt.0.) then
	960	c qm(1)=q(1)
	961	c else
	962	c qm(1)=0.
	963	c end if
	964
[2009]	965	END SUBROUTINE vl1d
[2184]	966
	967	c *****************************************************************
	968	SUBROUTINE scavenging_by_co2(ngrid,nlayer,nq,ptimestep,pplev,pq,
	969	& rdust,pcondicea,pfallice,pdq_scav,pdqsc)
	970
	971	c
	972	c
	973	c Calcul de la quantite de poussiere lessivee par les nuages de CO2
	974	c
	975	c --------------------------------------------------------------------
	976	use tracer_mod, only: nqmx, igcm_h2o_vap, igcm_h2o_ice,
	977	& igcm_dust_mass, igcm_dust_number,
	978	& igcm_ccn_mass, igcm_ccn_number,
	979	& rho_dust, nuice_sed, nuice_ref,r3n_q
	980	use comcstfi_h, only: g
[2409]	981	use dust_param_mod, only: freedust
[2184]	982	IMPLICIT NONE
[2409]	983	include "callkeys.h" ! for the flags water and microphys
[2184]	984	c
	985	c
	986	c Arguments:
	987	INTEGER,INTENT(IN) :: ngrid ! number of atmospheric columns
	988	INTEGER,INTENT(IN) :: nlayer ! number of vertical layers
	989	INTEGER,INTENT(IN) :: nq ! number of tracers
	990	REAL,INTENT(IN) :: ptimestep ! physics timestep (s)
	991	REAL,INTENT(IN) :: pplev(ngrid,nlayer+1) ! inter-layer pressure (Pa)
	992	REAL,INTENT(IN) :: pq(ngrid,nlayer,nq)
	993	REAL,INTENT(IN) :: rdust(ngrid,nlayer) ! dust effective radius
[2601]	994	REAL,INTENT(IN) :: pcondicea(ngrid,nlayer) ! condensation rate in layer l (kg/m2/s)
[2184]	995	REAL,INTENT(IN) :: pfallice(ngrid,nlayer+1) ! amount of ice falling from layer l (kg/m2/s)
	996
[2601]	997	REAL,INTENT(OUT) :: pdq_scav(ngrid,nlayer,nq) ! tendency due to scavenging by co2
[2184]	998	REAL,INTENT(OUT) :: pdqsc(ngrid,nq) ! tendency on surface tracers
	999
	1000	c Locals:
	1001	INTEGER l,ig
	1002	REAL scav_ratio_dust, scav_ratio_wice ! ratio of the dust/water ice mass mixing ratios in condensing CO2 ice and in air
	1003	REAL scav_dust_mass(nlayer+1) ! dust flux (mass) scavenged towards the lower layer (kg/m2/s) (POSITIVE WHEN DOWNWARD)
	1004	REAL scav_dust_number(nlayer+1) ! dust flux (number) scavenged towards the lower layer (kg/m2/s) (POSITIVE WHEN DOWNWARD)
	1005	REAL scav_ccn_mass(nlayer+1) ! ccn flux (mass) scavenged towards the lower layer
	1006	REAL scav_ccn_number(nlayer+1) ! ccn flux (number) scavenged towards the lower layer
	1007	REAL scav_h2o_ice(nlayer+1) ! water ice flux (mass) scavenged towards the lower layer
	1008	REAL massl ! mass of the layer l at point ig (kg/m2)
	1009
	1010	c Initialization:
[2669]	1011	scav_ratio_dust = 20 !1 !10 !100 !1000 ! the scavenging ratio value of 20 is a good compromise to remove the dust in the polar night
	1012	scav_ratio_wice = scav_ratio_dust ! while not drying up the water cycle (which occurs at scav_ratio_wice values above 50 at least)
[2184]	1013	pdq_scav(:,:,:)=0.
[2601]	1014	pdqsc(:,:)=0.
	1015
[2184]	1016	DO ig=1,ngrid
	1017	scav_dust_mass(nlayer+1)=0.
	1018	scav_dust_number(nlayer+1)=0.
	1019	scav_ccn_mass(nlayer+1)=0.
	1020	scav_ccn_number(nlayer+1)=0.
	1021	scav_h2o_ice(nlayer+1)=0.
	1022
	1023	DO l=nlayer , 1, -1
	1024	massl=(pplev(ig,l)-pplev(ig,l+1))/g
	1025	IF(pcondicea(ig,l).GT.0.)THEN ! if CO2 condenses and traps dust/water ice
	1026	! Calculation of the tendencies
	1027	if (freedust) then
	1028	pdq_scav(ig,l,igcm_dust_mass)=-pq(ig,l,igcm_dust_mass)
	1029	& /ptimestep*(1-exp(
	1030	& -scav_ratio_dustpcondicea(ig,l)ptimestep/massl))
	1031
	1032	pdq_scav(ig,l,igcm_dust_number)=pdq_scav(ig,l,igcm_dust_mass)
	1033	& *r3n_q/rdust(ig,l)
	1034	endif
	1035	if (freedust.AND.microphys) then
	1036	pdq_scav(ig,l,igcm_ccn_mass)=-pq(ig,l,igcm_ccn_mass)
	1037	& /ptimestep*(1-exp(
	1038	& -scav_ratio_wicepcondicea(ig,l)ptimestep/massl))
	1039	pdq_scav(ig,l,igcm_ccn_number)=pdq_scav(ig,l,igcm_ccn_mass)
	1040	& *r3n_q/rdust(ig,l)
	1041	endif
	1042	if (water) then
	1043	pdq_scav(ig,l,igcm_h2o_ice)=-pq(ig,l,igcm_h2o_ice)
	1044	& /ptimestep*(1-exp(
	1045	& -scav_ratio_wicepcondicea(ig,l)ptimestep/massl))
	1046	endif
	1047
	1048	ELSE IF(pcondicea(ig,l).LT.0.)THEN ! if CO2 sublimates and releases dust/water ice
	1049	! Calculation of the tendencies
	1050	if (freedust) then
	1051	pdq_scav(ig,l,igcm_dust_mass)=-pcondicea(ig,l)/massl*
	1052	& scav_dust_mass(l+1)/pfallice(ig,l+1)
	1053
	1054	pdq_scav(ig,l,igcm_dust_number)=-pcondicea(ig,l)/massl*
	1055	& scav_dust_number(l+1)/pfallice(ig,l+1)
	1056	endif
	1057	if (freedust.AND.microphys) then
	1058	pdq_scav(ig,l,igcm_ccn_mass)=-pcondicea(ig,l)/massl*
	1059	& scav_ccn_mass(l+1)/pfallice(ig,l+1)
	1060
	1061	pdq_scav(ig,l,igcm_ccn_number)=-pcondicea(ig,l)/massl*
	1062	& scav_ccn_number(l+1)/pfallice(ig,l+1)
	1063	endif
	1064	if (water) then
	1065	pdq_scav(ig,l,igcm_h2o_ice)=-pcondicea(ig,l)/massl*
	1066	& scav_h2o_ice(l+1)/pfallice(ig,l+1)
	1067	endif
	1068
	1069	END IF
	1070	! Calculation of the scavenged dust/wice flux towards the lower layers
	1071	if (freedust) then
	1072	scav_dust_mass(l)=-pdq_scav(ig,l,igcm_dust_mass)*massl
	1073	& +scav_dust_mass(l+1)
	1074
	1075	scav_dust_number(l)=-pdq_scav(ig,l,igcm_dust_number)*massl
	1076	& +scav_dust_number(l+1)
	1077	endif
	1078	if (freedust.AND.microphys) then
	1079	scav_ccn_mass(l)=-pdq_scav(ig,l,igcm_ccn_mass)*massl
	1080	& +scav_ccn_mass(l+1)
	1081
	1082	scav_ccn_number(l)=-pdq_scav(ig,l,igcm_ccn_number)*massl
	1083	& +scav_dust_number(l+1)
	1084	endif
	1085	if (water) then
	1086	scav_h2o_ice(l)=-pdq_scav(ig,l,igcm_h2o_ice)*massl
	1087	& +scav_h2o_ice(l+1)
	1088	endif
	1089
	1090	ENDDO
	1091	! Calculation of the surface tendencies
	1092	if (freedust) then
	1093	pdqsc(ig,igcm_dust_mass)=pdqsc(ig,igcm_dust_mass)
	1094	& +scav_dust_mass(1)
	1095	pdqsc(ig,igcm_dust_number)=pdqsc(ig,igcm_dust_number)
	1096	& +scav_dust_number(1)
	1097	endif
	1098	if (freedust.AND.microphys) then
	1099	pdqsc(ig,igcm_dust_mass)=pdqsc(ig,igcm_dust_mass)
	1100	& +scav_ccn_mass(1)
	1101	pdqsc(ig,igcm_dust_number)=pdqsc(ig,igcm_dust_number)
	1102	& +scav_ccn_number(1)
	1103	endif
	1104	if (water) then
	1105	pdqsc(ig,igcm_h2o_ice)=scav_h2o_ice(1)
	1106	endif
[2601]	1107
	1108	ENDDO ! loop on ngrid
[2184]	1109
	1110	END SUBROUTINE scavenging_by_co2
	1111
[2009]	1112	END MODULE co2condens_mod

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

Download in other formats:

Original Format