Context Navigation

co2condens_mod.F @ 2437

Last change on this file since 2437 was 2409, checked in by emillour, 5 years ago
Mars GCM: Code tidying : make a "dust_param_mod" module to store dust cycle relevant flags and variables (and remove them from callkeys.h) EM
File size: 36.4 KB

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

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

Download in other formats:

Original Format