Context Navigation

co2condens_mod.F @ 2646

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

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

Download in other formats:

Original Format