Context Navigation

co2condens_mod.F @ 4086

Last change on this file since 4086 was 3901, checked in by emillour, 8 months ago

Mars PCM:
Some code tidying:

turn aeroptproperties, albedocaps, cvmgp and convadj into modules
remove useless check in callradite
clean module vlz_fi (remove obsolete #ifdef CRAY cpp alternatives)
remove function cvmgp since it was only called under #ifdef CRAY in vlz_fi

File size: 44.9 KB

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

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

Download in other formats:

Original Format