Context Navigation

co2condens_mod.F @ 2184

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

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

Download in other formats:

Original Format