Context Navigation

co2condens_mod.F @ 2383

Last change on this file since 2383 was 2332, checked in by mvals, 5 years ago

Mars GCM:
Follow-up of the last commit for the transport of the isotopic ratio: simplification of the transmission of variables from the dynamics to the
physics:

libf/dynphy_lonlat/phymars/: iniphysiq_mod.F90: transmission of the content of 2 variables describing the isotopes instead of 4 (nqperes: number of tracers "peres", nqfils:

number of tracers "fils")

libf/phymars/: phys_state_var_init_mod.F90, tracer_mod.F: idem callsedim_mod.F: idem co2condens_mod.F: idem
libf/phymars/dyn1d: testphys1d.F: idem (the reading interface for traceur.def has been completed to fill the variables nqperes and nqfils).

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

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

Download in other formats:

Original Format