Context Navigation

co2condens_mod.F @ 2825

Last change on this file since 2825 was 2823, checked in by emillour, 3 years ago
Mars GCM: Remove the "tracer" (logical) flag as we now always run with at least one tracer. EM
File size: 38.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
[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
[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
	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
[2823]	192	c Prepare Special treatment if one of the tracer is CO2 gas
	193	do iq=1,nq
[38]	194	if (noms(iq).eq."co2") then
	195	ico2=iq
	196	m_co2 = 44.01E-3 ! CO2 molecular mass (kg/mol)
	197	m_noco2 = 33.37E-3 ! Non condensible mol mass (kg/mol)
	198	c Compute A and B coefficient use to compute
	199	c mean molecular mass Mair defined by
	200	c 1/Mair = q(ico2)/m_co2 + (1-q(ico2))/m_noco2
	201	c 1/Mair = A*q(ico2) + B
	202	A =(1/m_co2 - 1/m_noco2)
	203	B=1/m_noco2
	204	endif
[2823]	205	enddo
[890]	206	ENDIF ! of IF (firstcall)
[38]	207	zcpi=1./cpp
[1130]	208
[38]	209	c
	210	c======================================================================
	211	c Calcul of CO2 condensation sublimation
	212	c ============================================================
	213	c
	214	c Used variable :
	215	c piceco2(ngrid) : amount of co2 ice on the ground (kg/m2)
	216	c zcondicea(ngrid,l): condensation rate in layer l (kg/m2/s)
	217	c zcondices(ngrid): condensation rate on the ground (kg/m2/s)
	218	c zfallice(ngrid,l):amount of ice falling from layer l (kg/m2/s)
	219	c
[1047]	220	c pdtc(ngrid,nlayer) : dT/dt due to cond/sub
[38]	221	c
	222	c
[2184]	223	c Tendencies set to 0
[38]	224	c -------------------------------------
[2189]	225	zcondicea(1:ngrid,1:nlayer) = 0.
	226	zfallice(1:ngrid,1:nlayer+1) = 0.
	227	pduc(1:ngrid,1:nlayer) = 0
	228	pdvc(1:ngrid,1:nlayer) = 0
	229	pdtc(1:ngrid,1:nlayer) = 0.
	230	pdqsc(1:ngrid,1:nq) = 0
[38]	231
[2189]	232	pdqc(1:ngrid,1:nlayer,1:nq) = 0
[38]	233
[2189]	234	zcondices(1:ngrid) = 0.
	235	pdtsrfc(1:ngrid) = 0.
	236	pdpsrf(1:ngrid) = 0.
	237	condsub(1:ngrid) = .false.
	238	zdiceco2(1:ngrid) = 0.
	239
[38]	240	zfallheat=0
[2184]	241
	242	zdq_scav(:,:,:)=0.
[38]	243
[2184]	244	c Update tendencies from previous processes
	245	c -------------------------------------
	246	DO l=1,nlayer
	247	DO ig=1,ngrid
	248	zt(ig,l)=pt(ig,l)+ pdt(ig,l)*ptimestep
	249	do iq=1,nq
	250	zq(ig,l,iq)=pq(ig,l,iq)+pdq(ig,l,iq)*ptimestep
	251	enddo
	252	ENDDO
	253	ENDDO
	254
[38]	255	c *************************
	256	c ATMOSPHERIC CONDENSATION
	257	c *************************
	258
	259	c Compute CO2 Volume mixing ratio
	260	c -------------------------------
	261	if (improved_ztcond.and.(ico2.ne.0)) then
	262	DO l=1,nlayer
	263	DO ig=1,ngrid
	264	qco2=pq(ig,l,ico2)+pdq(ig,l,ico2)*ptimestep
	265	c Mean air molecular mass = 1/(q(ico2)/m_co2 + (1-q(ico2))/m_noco2)
	266	mmean=1/(A*qco2 +B)
	267	vmr_co2(ig,l) = qco2*mmean/m_co2
	268	ENDDO
	269	ENDDO
	270	else
	271	DO l=1,nlayer
	272	DO ig=1,ngrid
	273	vmr_co2(ig,l)=0.95
	274	ENDDO
	275	ENDDO
[2184]	276	endif
[38]	277
[1996]	278	IF (.NOT. co2clouds) then
[38]	279	c forecast of atmospheric temperature zt and frost temperature ztcond
	280	c --------------------------------------------------------------------
	281
	282	DO l=1,nlayer
	283	DO ig=1,ngrid
	284	! ztcond(ig,l)=1./(bcond-acondlog(.0095pplay(ig,l)))
[1263]	285	if (pplay(ig,l).ge.1e-4) then
	286	ztcond(ig,l)=
[38]	287	& 1./(bcond-acondlog(.01vmr_co2(ig,l)*pplay(ig,l)))
[1263]	288	else
	289	ztcond(ig,l)=0.0 !mars Monica
	290	endif
[38]	291	ENDDO
	292	ENDDO
[327]	293
[328]	294	ztcond(:,nlayer+1)=ztcond(:,nlayer)
[38]	295
	296	c Condensation/sublimation in the atmosphere
	297	c ------------------------------------------
	298	c (calcul of zcondicea , zfallice and pdtc)
	299	c
	300	DO l=nlayer , 1, -1
	301	DO ig=1,ngrid
	302	pdtc(ig,l)=0.
	303	IF((zt(ig,l).LT.ztcond(ig,l)).or.(zfallice(ig,l+1).gt.0))THEN
	304	condsub(ig)=.true.
	305	IF (zfallice(ig,l+1).gt.0) then
	306	zfallheat=zfallice(ig,l+1)*
	307	& (pphi(ig,l+1)-pphi(ig,l) +
	308	& cpice*(ztcond(ig,l+1)-ztcond(ig,l)))/latcond
	309	ELSE
	310	zfallheat=0.
	311	ENDIF
	312	pdtc(ig,l)=(ztcond(ig,l) - zt(ig,l))/ptimestep
	313	zcondicea(ig,l)=(pplev(ig,l)-pplev(ig,l+1))
	314	& ccondpdtc(ig,l)- zfallheat
	315	c Case when the ice from above sublimes entirely
	316	c """""""""""""""""""""""""""""""""""""""""""""""
	317	IF (zfallice(ig,l+1).lt.- zcondicea(ig,l)) then
	318	pdtc(ig,l)=(-zfallice(ig,l+1)+zfallheat)/
	319	& (ccond*(pplev(ig,l)-pplev(ig,l+1)))
	320	zcondicea(ig,l)= -zfallice(ig,l+1)
	321	END IF
	322
	323	zfallice(ig,l) = zcondicea(ig,l)+zfallice(ig,l+1)
	324	END IF
	325	ENDDO
	326	ENDDO
[1996]	327
[2599]	328	condens_layer(:,:) = zcondicea(:,:)
	329	condens_column(:) = 0.
	330
[2184]	331	if (scavco2cond) then
	332	call scavenging_by_co2(ngrid,nlayer,nq,ptimestep,pplev,zq,
	333	& rdust,zcondicea,zfallice,zdq_scav,pdqsc)
	334	endif
[2660]	335	call WRITEdiagfi(ngrid,"zfallice",
	336	& "",
	337	& " ",2,zfallice(ngrid,1))
[2184]	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
[2660]	352	call WRITEdiagfi(ngrid,"zfallice",
	353	& "",
	354	& " ",2,zdqssed_co2) ! otherwise we have not 1 day(1proc) = 1 day (24procs) test
[2184]	355	ENDIF ! end of if co2clouds
[38]	356
[1996]	357	call WRITEdiagfi(ngrid,"pdtc_atm",
	358	& "temperature tendency due to CO2 condensation",
	359	& " ",3,pdtc)
	360
[2599]	361	call WRITEdiagfi(ngrid,"condens_layer",
[1996]	362	& "",
[2599]	363	& " ",3,condens_layer)
[1996]	364
[38]	365	c *************************
	366	c SURFACE CONDENSATION
	367	c *************************
	368
	369	c forecast of ground temperature ztsrf and frost temperature ztcondsol
	370	c --------------------------------------------------------------------
	371	DO ig=1,ngrid
	372	ztcondsol(ig)=
	373	& 1./(bcond-acondlog(.01vmr_co2(ig,1)*pplev(ig,1)))
	374	ztsrf(ig) = ptsrf(ig) + pdtsrf(ig)*ptimestep
	375	ENDDO
	376
	377	c
	378	c Condensation/sublimation on the ground
	379	c --------------------------------------
[1114]	380	c (compute zcondices and pdtsrfc)
[38]	381	c
[2566]	382	c No microphysics of CO2 clouds
[38]	383	DO ig=1,ngrid
[1541]	384	IF(latitude(ig).lt.0) THEN
[1114]	385	! Southern hemisphere
[38]	386	icap=2
	387	ELSE
[1114]	388	! Northern hemisphere
[38]	389	icap=1
	390	ENDIF
	391	c
	392	c Loop on where we have condensation/ sublimation
	393	IF ((ztsrf(ig) .LT. ztcondsol(ig)) .OR. ! ground cond
	394	$ (zfallice(ig,1).NE.0.) .OR. ! falling snow
	395	$ ((ztsrf(ig) .GT. ztcondsol(ig)) .AND. ! ground sublim.
	396	$ ((piceco2(ig)+zfallice(ig,1)*ptimestep) .NE. 0.))) THEN
[2599]	397	condsub(ig) = .true.
[38]	398
[2599]	399	IF (zfallice(ig,1).gt.0 .and. .not. co2clouds) then
[2739]	400	! potential eneregy release due to the impact of the snowfall
	401	! appendix of forget et al. 1999
	402	zfallheat = zfallice(ig,1) * (pphi(ig,1) +
[2599]	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)
[2660]	443	if (co2clouds) then
	444	IF((piceco2(ig)/ptimestep).LE.
	445	& -zcondices(ig))THEN
	446	zcondices(ig) = -piceco2(ig)/ptimestep
	447	pdtsrfc(ig)=(latcond/pcapcal(ig)) *
	448	& (zcondices(ig)+zfallheat)
	449	END IF
	450	else
[38]	451	IF((piceco2(ig)/ptimestep+zfallice(ig,1)).LE.
	452	& -zcondices(ig))THEN
[2599]	453	zcondices(ig) = -piceco2(ig)/ptimestep - zfallice(ig,1)
[2660]	454	pdtsrfc(ig)=(latcond/pcapcal(ig)) *
[38]	455	& (zcondices(ig)+zfallheat)
	456	END IF
[2660]	457	end if
[38]	458
	459	c Changing CO2 ice amount and pressure :
	460	c """"""""""""""""""""""""""""""""""""
[2599]	461	zdiceco2(ig) = zcondices(ig) + zfallice(ig,1)
	462	& + condens_column(ig)
	463	if (co2clouds) then
	464	! add here only direct condensation/sublimation
	465	piceco2(ig) = piceco2(ig) + zcondices(ig)*ptimestep
	466	else
	467	! add here also CO2 ice in the atmosphere
[38]	468	piceco2(ig) = piceco2(ig) + zdiceco2(ig)*ptimestep
[2599]	469	end if
[38]	470	pdpsrf(ig) = -zdiceco2(ig)*g
	471
	472	IF(ABS(pdpsrf(ig)*ptimestep).GT.pplev(ig,1)) THEN
	473	PRINT*,'STOP in condens'
	474	PRINT*,'condensing more than total mass'
	475	PRINT*,'Grid point ',ig
[2124]	476	PRINT*,'Longitude(degrees): ',longitude_deg(ig)
	477	PRINT*,'Latitude(degrees): ',latitude_deg(ig)
[38]	478	PRINT*,'Ps = ',pplev(ig,1)
	479	PRINT*,'d Ps = ',pdpsrf(ig)
[2399]	480	CALL abort_physic('co2condens',
	481	& 'condensing more than total mass', 1)
[38]	482	ENDIF
[2184]	483	END IF ! if there is condensation/sublimation
[38]	484	ENDDO ! of DO ig=1,ngrid
	485
	486	c ********************************************************************
	487	c Surface albedo and emissivity of the surface below the snow (emisref)
	488	c ********************************************************************
	489
	490	! Check that amont of CO2 ice is not problematic
	491	DO ig=1,ngrid
	492	if(.not.piceco2(ig).ge.0.) THEN
	493	if(piceco2(ig).le.-5.e-8) print*,
[2124]	494	$ 'WARNING co2condens piceco2(',ig,')=', piceco2(ig)
[38]	495	piceco2(ig)=0.
	496	endif
	497	ENDDO
	498
	499	! Set albedo and emissivity of the surface
	500	! ----------------------------------------
	501	CALL albedocaps(zls,ngrid,piceco2,psolaralb,emisref)
	502
	503	! set pemisurf() to emissiv when there is bare surface (needed for co2snow)
	504	DO ig=1,ngrid
	505	if (piceco2(ig).eq.0) then
	506	pemisurf(ig)=emissiv
	507	endif
	508	ENDDO
	509
	510	! firstcall2=.false.
	511	c ***************************************************************
	512	c Correction to account for redistribution between sigma or hybrid
	513	c layers when changing surface pressure (and warming/cooling
	514	c of the CO2 currently changing phase).
	515	c *************************************************************
	516
	517	DO ig=1,ngrid
	518	if (condsub(ig)) then
	519	do l=1,nlayer
	520	ztc(l) =zt(ig,l) +pdtc(ig,l) *ptimestep
	521	zu(l) =pu(ig,l) +pdu( ig,l) *ptimestep
	522	zv(l) =pv(ig,l) +pdv( ig,l) *ptimestep
[1036]	523	do iq=1,nq
[2601]	524	zqc(l,iq)=zq(ig,l,iq)+zdq_scav(ig,l,iq)*ptimestep ! zdq_scav=0 if co2clouds=true
[38]	525	enddo
[2184]	526	enddo
[38]	527
	528	c Mass fluxes through the sigma levels (kg.m-2.s-1) (>0 when up)
	529	c """"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
[2599]	530	zmflux(1) = -zcondices(ig) - zdqssed_co2(ig)
[2566]	531	DO l=1,nlayer
[2599]	532	zmflux(l+1) = zmflux(l) - condens_layer(ig,l)
[1438]	533	#ifndef MESOSCALE
[2599]	534	& + (bp(l)-bp(l+1))*(-pdpsrf(ig)/g)
[2566]	535	c zmflux set to 0 if very low to avoid: top layer is disappearing in v1ld
	536	if (abs(zmflux(l+1)).lt.1E-13.OR.bp(l+1).eq.0.) then
	537	zmflux(l+1)=0.
	538	end if
[1438]	539	#else
[2599]	540	zmflux(l+1) = zmflux(l) - condens_layer(ig,l)
[2566]	541	if (abs(zmflux(l+1)).lt.1E-13) zmflux(l+1)=0.
	542	PRINT*, "MESOSCALE: FLUX THROUGH SIGMA LEVELS from"//
	543	& "dPS HAVE TO BE IMPLEMENTED"
[1438]	544	#endif
[2566]	545	END DO
[2155]	546	#ifdef MESOSCALE
	547	print*,"absurd mass set because bp not available"
	548	print*,"TO BE FIXED"
	549	#else
[2124]	550	c Mass of each layer at the end of timestep
	551	c -----------------------------------------
[38]	552	DO l=1,nlayer
[2124]	553	masse(l)=( pplev(ig,l) - pplev(ig,l+1) +
	554	& (bp(l)-bp(l+1))pdpsrf(ig)ptimestep)/g
[38]	555	END DO
[2155]	556	#endif
[38]	557
	558	c Corresponding fluxes for T,U,V,Q
	559	c """"""""""""""""""""""""""""""""
	560
	561	c averaging operator for TRANSPORT
	562	c """"""""""""""""""""""""""""""""
	563	c Value transfert at the surface interface when condensation
	564	c sublimation:
	565	ztm(1) = ztsrf(ig) + pdtsrfc(ig)*ptimestep
	566	zum(1) = 0
	567	zvm(1) = 0
[1036]	568	do iq=1,nq
[38]	569	zqm(1,iq)=0. ! most tracer do not condense !
	570	enddo
	571	c Special case if one of the tracer is CO2 gas
	572	if (ico2.ne.0) zqm(1,ico2)=1. ! flux is 100% CO2
	573
	574	c Van Leer scheme:
	575	DO l=1,nlayer+1
	576	w(l)=-zmflux(l)*ptimestep
	577	END DO
[1269]	578	call vl1d(nlayer,ztc,2.,masse,w,ztm)
	579	call vl1d(nlayer,zu ,2.,masse,w,zum)
	580	call vl1d(nlayer,zv ,2.,masse,w,zvm)
[2322]	581	! MVals: loop over the fathers ("peres")
	582	do iq=1,nqperes
[38]	583	do l=1,nlayer
[2184]	584	zq1(l)=zqc(l,iq)
[38]	585	enddo
	586	zqm1(1)=zqm(1,iq)
[1269]	587	call vl1d(nlayer,zq1,2.,masse,w,zqm1)
[38]	588	do l=2,nlayer
[2184]	589	zqc(l,iq)=zq1(l)
[38]	590	zqm(l,iq)=zqm1(l)
	591	enddo
[2322]	592	! MVals: loop over the sons ("fils")
	593	if (nqfils(iq).gt.0) then
[2332]	594	if (iq.eq.igcm_h2o_ice) then
	595	iq2=igcm_hdo_ice
	596	else if (iq.eq.igcm_h2o_vap) then
	597	iq2=igcm_hdo_vap
	598	else
	599	call abort_physic("co2condens_mod","invalid isotope",1)
	600	endif
[2322]	601	do l=1,nlayer
	602	if (zqc(l,iq).gt.qperemin) then
	603	Ratio1(l)=zqc(l,iq2)/zqc(l,iq)
	604	else
	605	Ratio1(l)=0.
	606	endif
	607	masseq(l)=max(masse(l)*zqc(l,iq),masseqmin)
	608	wq(l)=w(l)*zqm(l,iq)
	609	enddo
	610	Ratiom1(1)=zqm(1,iq2)
	611	call vl1d(nlayer,Ratio1,2.,masseq,wq,Ratiom1)
	612	zqm(1,iq2)=Ratiom1(1)*zqc(1,iq)
	613	do l=2,nlayer
	614	zqm(l,iq2)=Ratiom1(l)*zqm(l,iq)
	615	enddo
	616	endif !if (nqfils(iq).gt.0) then
	617	enddo !iq=1,nqperes
[38]	618
	619	c Surface condensation affects low winds
	620	if (zmflux(1).lt.0) then
	621	zum(1)= zu(1) * (w(1)/masse(1))
	622	zvm(1)= zv(1) * (w(1)/masse(1))
	623	if (w(1).gt.masse(1)) then ! ensure numerical stability
	624	zum(1)= (zu(1)-zum(2))*masse(1)/w(1) + zum(2)
	625	zvm(1)= (zv(1)-zvm(2))*masse(1)/w(1) + zvm(2)
	626	end if
	627	end if
	628
	629	ztm(nlayer+1)= ztc(nlayer) ! should not be used, but...
	630	zum(nlayer+1)= zu(nlayer) ! should not be used, but...
	631	zvm(nlayer+1)= zv(nlayer) ! should not be used, but...
[1036]	632	do iq=1,nq
[2184]	633	zqm(nlayer+1,iq)= zqc(nlayer,iq)
[38]	634	enddo
[86]	635
	636	#ifdef MESOSCALE
	637	!!!! AS: This part must be commented in the mesoscale model
	638	!!!! AS: ... to avoid instabilities.
	639	!!!! AS: you have to compile with -DMESOSCALE to do so
	640	#else
[38]	641	c Tendencies on T, U, V, Q
	642	c """"""""""""""""""""""""
	643	DO l=1,nlayer
[1996]	644	IF(.not. co2clouds) THEN
[38]	645	c Tendencies on T
	646	zdtsig(ig,l) = (1/masse(l)) *
	647	& ( zmflux(l)*(ztm(l) - ztc(l))
	648	& - zmflux(l+1)*(ztm(l+1) - ztc(l))
[2599]	649	& + condens_layer(ig,l)*(ztcond(ig,l)-ztc(l)) )
[1996]	650	ELSE
	651	zdtsig(ig,l) = (1/masse(l)) *
	652	& ( zmflux(l)*(ztm(l) - ztc(l))
	653	& - zmflux(l+1)*(ztm(l+1) - ztc(l)))
	654	ENDIF
	655	c D.BARDET: for diagnotics
	656	zmflux3D(ig,l)=zmflux(l)
	657	ztm3D(ig,l)=ztm(l)
	658	ztc3D(ig,l)=ztc(l)
	659
[38]	660	pdtc(ig,l) = pdtc(ig,l) + zdtsig(ig,l)
	661
	662	c Tendencies on U
	663	pduc(ig,l) = (1/masse(l)) *
	664	& ( zmflux(l)*(zum(l) - zu(l))
	665	& - zmflux(l+1)*(zum(l+1) - zu(l)) )
	666
	667
	668	c Tendencies on V
	669	pdvc(ig,l) = (1/masse(l)) *
	670	& ( zmflux(l)*(zvm(l) - zv(l))
	671	& - zmflux(l+1)*(zvm(l+1) - zv(l)) )
	672
	673	END DO
[1114]	674
[86]	675	#endif
[38]	676
[2566]	677	do iq=1,nq
	678	! if (noms(iq).eq.'co2') then
	679	if (iq.eq.ico2) then
	680	c SPECIAL Case when the tracer IS CO2 :
	681	DO l=1,nlayer
	682	pdqc(ig,l,iq)= (1/masse(l)) *
	683	& ( zmflux(l)*(zqm(l,iq) - zqc(l,iq))
	684	& - zmflux(l+1)*(zqm(l+1,iq) - zqc(l,iq))
[2599]	685	& + condens_layer(ig,l)*(zqc(l,iq)-1.) )
[2566]	686	END DO
	687	else
	688	DO l=1,nlayer
	689	pdqc(ig,l,iq)= (1/masse(l)) *
	690	& ( zmflux(l)*(zqm(l,iq) - zqc(l,iq))
	691	& - zmflux(l+1)*(zqm(l+1,iq) - zqc(l,iq))
[2599]	692	& + condens_layer(ig,l)*zqc(l,iq) )
[2566]	693
[2601]	694	pdqc(ig,l,iq)=pdqc(ig,l,iq)+zdq_scav(ig,l,iq) ! zdq_scav=0 if co2clouds=true
[2566]	695	END DO
	696	end if
	697	enddo
[38]	698
	699	end if ! if (condsub)
	700	END DO ! loop on ig
	701
	702	c ***************************************************************
	703	c CO2 snow / clouds scheme
	704	c ***************************************************************
[2566]	705	call co2snow(ngrid,nlayer,ptimestep,emisref,condsub,pplev,
[2599]	706	& condens_layer,zcondices,zfallice,pemisurf)
[38]	707	c ***************************************************************
	708	c Ecriture des diagnostiques
	709	c ***************************************************************
	710
	711	c DO l=1,nlayer
	712	c DO ig=1,ngrid
	713	c Taux de cond en kg.m-2.pa-1.s-1
	714	c tconda1(ig,l)=zcondicea(ig,l)/(pplev(ig,l)-pplev(ig,l+1))
	715	c Taux de cond en kg.m-3.s-1
	716	c tconda2(ig,l)=tconda1(ig,l)pplay(ig,l)g/(r*pt(ig,l))
	717	c END DO
	718	c END DO
[1047]	719	c call WRITEDIAGFI(ngrid,'tconda1',
[38]	720	c &'Taux de condensation CO2 atmospherique /Pa',
	721	c & 'kg.m-2.Pa-1.s-1',3,tconda1)
[1047]	722	c call WRITEDIAGFI(ngrid,'tconda2',
[38]	723	c &'Taux de condensation CO2 atmospherique /m',
	724	c & 'kg.m-3.s-1',3,tconda2)
	725
	726	! output falling co2 ice in 1st layer:
[1047]	727	! call WRITEDIAGFI(ngrid,'fallice',
[38]	728	! &'Precipitation of co2 ice',
	729	! & 'kg.m-2.s-1',2,zfallice(1,1))
	730
[1114]	731	#ifndef MESOSCALE
	732	! Extra special case for surface temperature tendency pdtsrfc:
	733	! we want to fix the south pole temperature to CO2 condensation temperature
	734	if (caps.and.(obliquit.lt.27.)) then
	735	! check if last grid point is the south pole
[1541]	736	if (abs(latitude(ngrid)-(-pi/2.)).lt.1.e-5) then
[1114]	737	! NB: Updated surface pressure, at grid point 'ngrid', is
	738	! ps(ngrid)=pplev(ngrid,1)+pdpsrf(ngrid)*ptimestep
[2124]	739	! write(,) "co2condens: South pole: latitude(ngrid)=",
[1541]	740	! & latitude(ngrid)
[1114]	741	ztcondsol(ngrid)=
	742	& 1./(bcond-acondlog(.01vmr_co2(ngrid,1)*
	743	& (pplev(ngrid,1)+pdpsrf(ngrid)*ptimestep)))
	744	pdtsrfc(ngrid)=(ztcondsol(ngrid)-ztsrf(ngrid))/ptimestep
	745	endif
	746	endif
	747	#endif
	748
[2009]	749	END SUBROUTINE co2condens
[38]	750
	751
	752
	753	c *****************************************************************
[1269]	754	SUBROUTINE vl1d(nlayer,q,pente_max,masse,w,qm)
[38]	755	c
	756	c
	757	c Operateur de moyenne inter-couche pour calcul de transport type
	758	c Van-Leer " pseudo amont " dans la verticale
	759	c q,w sont des arguments d'entree pour le s-pg ....
	760	c masse : masse de la couche Dp/g
	761	c w : masse d'atm ``transferee'' a chaque pas de temps (kg.m-2)
	762	c pente_max = 2 conseillee
	763	c
	764	c
	765	c --------------------------------------------------------------------
	766	IMPLICIT NONE
	767
	768	c
	769	c
	770	c
	771	c Arguments:
	772	c ----------
[1270]	773	integer nlayer
[1269]	774	real masse(nlayer),pente_max
	775	REAL q(nlayer),qm(nlayer+1)
	776	REAL w(nlayer+1)
[38]	777	c
	778	c Local
	779	c ---------
	780	c
	781	INTEGER l
	782	c
[1269]	783	real dzq(nlayer),dzqw(nlayer),adzqw(nlayer),dzqmax
[38]	784	real sigw, Mtot, MQtot
	785	integer m
	786	c integer ismax,ismin
	787
	788
	789	c On oriente tout dans le sens de la pression
	790	c W > 0 WHEN DOWN !!!!!!!!!!!!!
	791
[1269]	792	do l=2,nlayer
[38]	793	dzqw(l)=q(l-1)-q(l)
	794	adzqw(l)=abs(dzqw(l))
	795	enddo
	796
[1269]	797	do l=2,nlayer-1
[38]	798	if(dzqw(l)*dzqw(l+1).gt.0.) then
	799	dzq(l)=0.5*(dzqw(l)+dzqw(l+1))
	800	else
	801	dzq(l)=0.
	802	endif
	803	dzqmax=pente_max*min(adzqw(l),adzqw(l+1))
	804	dzq(l)=sign(min(abs(dzq(l)),dzqmax),dzq(l))
	805	enddo
	806
	807	dzq(1)=0.
[1269]	808	dzq(nlayer)=0.
[38]	809
[1269]	810	do l = 1,nlayer-1
[38]	811
	812	c Regular scheme (transfered mass < layer mass)
	813	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	814	if(w(l+1).gt.0. .and. w(l+1).le.masse(l+1)) then
	815	sigw=w(l+1)/masse(l+1)
	816	qm(l+1)=(q(l+1)+0.5(1.-sigw)dzq(l+1))
	817	else if(w(l+1).le.0. .and. -w(l+1).le.masse(l)) then
	818	sigw=w(l+1)/masse(l)
	819	qm(l+1)=(q(l)-0.5(1.+sigw)dzq(l))
	820
	821	c Extended scheme (transfered mass > layer mass)
	822	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	823	else if(w(l+1).gt.0.) then
	824	m=l+1
	825	Mtot = masse(m)
	826	MQtot = masse(m)*q(m)
[1269]	827	do while ((m.lt.nlayer).and.(w(l+1).gt.(Mtot+masse(m+1))))
[38]	828	m=m+1
	829	Mtot = Mtot + masse(m)
	830	MQtot = MQtot + masse(m)*q(m)
	831	end do
[1269]	832	if (m.lt.nlayer) then
[38]	833	sigw=(w(l+1)-Mtot)/masse(m+1)
	834	qm(l+1)= (1/w(l+1))(MQtot + (w(l+1)-Mtot)
	835	& (q(m+1)+0.5(1.-sigw)dzq(m+1)) )
	836	else
	837	w(l+1) = Mtot
	838	qm(l+1) = Mqtot / Mtot
[2399]	839	CALL abort_physic('co2condens',
	840	& 'top layer is disapearing !', 1)
[38]	841	end if
	842	else ! if(w(l+1).lt.0)
	843	m = l-1
	844	Mtot = masse(m+1)
	845	MQtot = masse(m+1)*q(m+1)
[120]	846	if (m.gt.0) then ! because some compilers will have problems
	847	! evaluating masse(0)
	848	do while ((m.gt.0).and.(-w(l+1).gt.(Mtot+masse(m))))
[38]	849	m=m-1
	850	Mtot = Mtot + masse(m+1)
	851	MQtot = MQtot + masse(m+1)*q(m+1)
[120]	852	if (m.eq.0) exit
	853	end do
	854	endif
[38]	855	if (m.gt.0) then
	856	sigw=(w(l+1)+Mtot)/masse(m)
	857	qm(l+1)= (-1/w(l+1))(MQtot + (-w(l+1)-Mtot)
	858	& (q(m)-0.5(1.+sigw)dzq(m)) )
	859	else
	860	qm(l+1)= (-1/w(l+1))(MQtot + (-w(l+1)-Mtot)qm(1))
	861	end if
	862	end if
	863	enddo
	864
[2124]	865	c boundary conditions (not used in co2condens !!)
[1269]	866	c qm(nlayer+1)=0.
[38]	867	c if(w(1).gt.0.) then
	868	c qm(1)=q(1)
	869	c else
	870	c qm(1)=0.
	871	c end if
	872
[2009]	873	END SUBROUTINE vl1d
[2184]	874
	875	c *****************************************************************
	876	SUBROUTINE scavenging_by_co2(ngrid,nlayer,nq,ptimestep,pplev,pq,
	877	& rdust,pcondicea,pfallice,pdq_scav,pdqsc)
	878
	879	c
	880	c
	881	c Calcul de la quantite de poussiere lessivee par les nuages de CO2
	882	c
	883	c --------------------------------------------------------------------
	884	use tracer_mod, only: nqmx, igcm_h2o_vap, igcm_h2o_ice,
	885	& igcm_dust_mass, igcm_dust_number,
	886	& igcm_ccn_mass, igcm_ccn_number,
	887	& rho_dust, nuice_sed, nuice_ref,r3n_q
	888	use comcstfi_h, only: g
[2409]	889	use dust_param_mod, only: freedust
[2184]	890	IMPLICIT NONE
[2409]	891	include "callkeys.h" ! for the flags water and microphys
[2184]	892	c
	893	c
	894	c Arguments:
	895	INTEGER,INTENT(IN) :: ngrid ! number of atmospheric columns
	896	INTEGER,INTENT(IN) :: nlayer ! number of vertical layers
	897	INTEGER,INTENT(IN) :: nq ! number of tracers
	898	REAL,INTENT(IN) :: ptimestep ! physics timestep (s)
	899	REAL,INTENT(IN) :: pplev(ngrid,nlayer+1) ! inter-layer pressure (Pa)
	900	REAL,INTENT(IN) :: pq(ngrid,nlayer,nq)
	901	REAL,INTENT(IN) :: rdust(ngrid,nlayer) ! dust effective radius
[2601]	902	REAL,INTENT(IN) :: pcondicea(ngrid,nlayer) ! condensation rate in layer l (kg/m2/s)
[2184]	903	REAL,INTENT(IN) :: pfallice(ngrid,nlayer+1) ! amount of ice falling from layer l (kg/m2/s)
	904
[2601]	905	REAL,INTENT(OUT) :: pdq_scav(ngrid,nlayer,nq) ! tendency due to scavenging by co2
[2184]	906	REAL,INTENT(OUT) :: pdqsc(ngrid,nq) ! tendency on surface tracers
	907
	908	c Locals:
	909	INTEGER l,ig
	910	REAL scav_ratio_dust, scav_ratio_wice ! ratio of the dust/water ice mass mixing ratios in condensing CO2 ice and in air
	911	REAL scav_dust_mass(nlayer+1) ! dust flux (mass) scavenged towards the lower layer (kg/m2/s) (POSITIVE WHEN DOWNWARD)
	912	REAL scav_dust_number(nlayer+1) ! dust flux (number) scavenged towards the lower layer (kg/m2/s) (POSITIVE WHEN DOWNWARD)
	913	REAL scav_ccn_mass(nlayer+1) ! ccn flux (mass) scavenged towards the lower layer
	914	REAL scav_ccn_number(nlayer+1) ! ccn flux (number) scavenged towards the lower layer
	915	REAL scav_h2o_ice(nlayer+1) ! water ice flux (mass) scavenged towards the lower layer
	916	REAL massl ! mass of the layer l at point ig (kg/m2)
	917
	918	c Initialization:
[2669]	919	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
	920	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]	921	pdq_scav(:,:,:)=0.
[2601]	922	pdqsc(:,:)=0.
	923
[2184]	924	DO ig=1,ngrid
	925	scav_dust_mass(nlayer+1)=0.
	926	scav_dust_number(nlayer+1)=0.
	927	scav_ccn_mass(nlayer+1)=0.
	928	scav_ccn_number(nlayer+1)=0.
	929	scav_h2o_ice(nlayer+1)=0.
	930
	931	DO l=nlayer , 1, -1
	932	massl=(pplev(ig,l)-pplev(ig,l+1))/g
	933	IF(pcondicea(ig,l).GT.0.)THEN ! if CO2 condenses and traps dust/water ice
	934	! Calculation of the tendencies
	935	if (freedust) then
	936	pdq_scav(ig,l,igcm_dust_mass)=-pq(ig,l,igcm_dust_mass)
	937	& /ptimestep*(1-exp(
	938	& -scav_ratio_dustpcondicea(ig,l)ptimestep/massl))
	939
	940	pdq_scav(ig,l,igcm_dust_number)=pdq_scav(ig,l,igcm_dust_mass)
	941	& *r3n_q/rdust(ig,l)
	942	endif
	943	if (freedust.AND.microphys) then
	944	pdq_scav(ig,l,igcm_ccn_mass)=-pq(ig,l,igcm_ccn_mass)
	945	& /ptimestep*(1-exp(
	946	& -scav_ratio_wicepcondicea(ig,l)ptimestep/massl))
	947	pdq_scav(ig,l,igcm_ccn_number)=pdq_scav(ig,l,igcm_ccn_mass)
	948	& *r3n_q/rdust(ig,l)
	949	endif
	950	if (water) then
	951	pdq_scav(ig,l,igcm_h2o_ice)=-pq(ig,l,igcm_h2o_ice)
	952	& /ptimestep*(1-exp(
	953	& -scav_ratio_wicepcondicea(ig,l)ptimestep/massl))
	954	endif
	955
	956	ELSE IF(pcondicea(ig,l).LT.0.)THEN ! if CO2 sublimates and releases dust/water ice
	957	! Calculation of the tendencies
	958	if (freedust) then
	959	pdq_scav(ig,l,igcm_dust_mass)=-pcondicea(ig,l)/massl*
	960	& scav_dust_mass(l+1)/pfallice(ig,l+1)
	961
	962	pdq_scav(ig,l,igcm_dust_number)=-pcondicea(ig,l)/massl*
	963	& scav_dust_number(l+1)/pfallice(ig,l+1)
	964	endif
	965	if (freedust.AND.microphys) then
	966	pdq_scav(ig,l,igcm_ccn_mass)=-pcondicea(ig,l)/massl*
	967	& scav_ccn_mass(l+1)/pfallice(ig,l+1)
	968
	969	pdq_scav(ig,l,igcm_ccn_number)=-pcondicea(ig,l)/massl*
	970	& scav_ccn_number(l+1)/pfallice(ig,l+1)
	971	endif
	972	if (water) then
	973	pdq_scav(ig,l,igcm_h2o_ice)=-pcondicea(ig,l)/massl*
	974	& scav_h2o_ice(l+1)/pfallice(ig,l+1)
	975	endif
	976
	977	END IF
	978	! Calculation of the scavenged dust/wice flux towards the lower layers
	979	if (freedust) then
	980	scav_dust_mass(l)=-pdq_scav(ig,l,igcm_dust_mass)*massl
	981	& +scav_dust_mass(l+1)
	982
	983	scav_dust_number(l)=-pdq_scav(ig,l,igcm_dust_number)*massl
	984	& +scav_dust_number(l+1)
	985	endif
	986	if (freedust.AND.microphys) then
	987	scav_ccn_mass(l)=-pdq_scav(ig,l,igcm_ccn_mass)*massl
	988	& +scav_ccn_mass(l+1)
	989
	990	scav_ccn_number(l)=-pdq_scav(ig,l,igcm_ccn_number)*massl
	991	& +scav_dust_number(l+1)
	992	endif
	993	if (water) then
	994	scav_h2o_ice(l)=-pdq_scav(ig,l,igcm_h2o_ice)*massl
	995	& +scav_h2o_ice(l+1)
	996	endif
	997
	998	ENDDO
	999	! Calculation of the surface tendencies
	1000	if (freedust) then
	1001	pdqsc(ig,igcm_dust_mass)=pdqsc(ig,igcm_dust_mass)
	1002	& +scav_dust_mass(1)
	1003	pdqsc(ig,igcm_dust_number)=pdqsc(ig,igcm_dust_number)
	1004	& +scav_dust_number(1)
	1005	endif
	1006	if (freedust.AND.microphys) then
	1007	pdqsc(ig,igcm_dust_mass)=pdqsc(ig,igcm_dust_mass)
	1008	& +scav_ccn_mass(1)
	1009	pdqsc(ig,igcm_dust_number)=pdqsc(ig,igcm_dust_number)
	1010	& +scav_ccn_number(1)
	1011	endif
	1012	if (water) then
	1013	pdqsc(ig,igcm_h2o_ice)=scav_h2o_ice(1)
	1014	endif
[2601]	1015
	1016	ENDDO ! loop on ngrid
[2184]	1017
	1018	END SUBROUTINE scavenging_by_co2
	1019
[2009]	1020	END MODULE co2condens_mod

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

Download in other formats:

Original Format