Context Navigation

condense_co2.F90 @ 1542

Last change on this file since 1542 was 1542, checked in by emillour, 9 years ago

Generic GCM:

fix for 1D in writediagfi to enable writing at "ecritphy" rate.
move iniprint.h to "misc"
Some code cleanup in anticipation of future updates:
- changed variable names in comgeomphy.F90: give them more explicit names: rlond => longitude , rlatd => latitude, airephy => cell_area, cuphy => dx , cvphy => dy
- removed long(), lati() and area() from comgeomfi_h.F90, use longitude(), latitude() and cell_are() from comgeomphy.F90 instead

Property svn:executable set to *

File size: 18.4 KB

Line
1	subroutine condense_co2(ngrid,nlayer,nq,ptimestep, &
2	pcapcal,pplay,pplev,ptsrf,pt, &
3	pdt,pdtsrf,pq,pdq, &
4	pqsurf,pdqsurfc,albedo,pemisurf, &
5	albedo_bareground,albedo_co2_ice_SPECTV, &
6	pdtc,pdtsrfc,pdpsrfc,pdqc)
7
8	use radinc_h, only : L_NSPECTV, naerkind
9	use gases_h, only: gfrac, igas_co2
10	use radii_mod, only : co2_reffrad
11	use aerosol_mod, only : iaero_co2
12	USE surfdat_h, only: emisice, emissiv
13	USE comgeomphy, only: latitude ! in radians
14	USE tracer_h, only: noms, rho_co2
15	use comcstfi_mod, only: g, r, cpp
16
17	implicit none
18
19	!==================================================================
20	! Purpose
21	! -------
22	! Condense and/or sublime CO2 ice on the ground and in the atmosphere, and sediment the ice.
23	!
24	! Inputs
25	! ------
26	! ngrid Number of vertical columns.
27	! nlayer Number of vertical layers.
28	! nq Number of tracers.
29	! ptimestep Duration of the physical timestep (s).
30	! pplay(ngrid,nlayer) Pressure layers (Pa).
31	! pplev(ngrid,nlayer+1) Pressure levels (Pa).
32	! pt(ngrid,nlayer) Atmospheric Temperatures (K).
33	! ptsrf(ngrid) Surface temperatures (K).
34	! pq(ngrid,nlayer,nq) Atmospheric tracers mixing ratios (kg/kg of air).
35	! pqsurf(ngrid,nq) Surface tracers (kg/m2).
36	!
37	! pdt(ngrid,nlayer) Time derivative before condensation/sublimation of pt.
38	! pdtsrf(ngrid) Time derivative before condensation/sublimation of ptsrf.
39	! pdq(ngrid,nlayer,nq) Time derivative before condensation/sublimation of
40	!
41	! albedo_bareground(ngrid) Albedo of the bare ground.
42	! albedo_co2_ice_SPECTV(L_NSPECTV) Spectral albedo of CO2 ice.
43	!
44	! Outputs
45	! -------
46	! pdpsrfc(ngrid) \ Contribution of condensation/sublimation
47	! pdtc(ngrid,nlayer) \ to the time derivatives of
48	! pdtsrfc(ngrid) / Surface Pressure, Atmospheric Temperatures,
49	! pdqsurfc(ngrid) / Surface Temperatures, Surface Tracers,
50	! pdqc(ngrid,nlayer,nq) / and Atmospheric Tracers.*
51	!
52	! pemisurf(ngrid) Emissivity of the surface.
53	!
54	! Both
55	! ----
56	! albedo(ngrid,L_NSPECTV) Spectral albedo of the surface.
57	!
58	! Authors
59	! -------
60	! Francois Forget (1996)
61	! Converted to Fortran 90 and slightly modified by R. Wordsworth (2009)
62	! Includes simplifed nucleation by J. Leconte (2011)
63	!
64	!==================================================================
65
66	!--------------------------
67	! Arguments
68	!--------------------------
69
70
71	INTEGER,INTENT(IN) :: ngrid
72	INTEGER,INTENT(IN) :: nlayer
73	INTEGER,INTENT(IN) :: nq
74	REAL,INTENT(IN) :: ptimestep
75	REAL,INTENT(IN) :: pcapcal(ngrid)
76	REAL,INTENT(IN) :: pplay(ngrid,nlayer)
77	REAL,INTENT(IN) :: pplev(ngrid,nlayer+1)
78	REAL,INTENT(IN) :: ptsrf(ngrid)
79	REAL,INTENT(IN) :: pt(ngrid,nlayer)
80	REAL,INTENT(IN) :: pdt(ngrid,nlayer)
81	REAL,INTENT(IN) :: pdtsrf(ngrid)
82	REAL,INTENT(IN) :: pq(ngrid,nlayer,nq)
83	REAL,INTENT(IN) :: pqsurf(ngrid,nq)
84	REAL,INTENT(IN) :: pdq(ngrid,nlayer,nq)
85	REAL,INTENT(IN) :: albedo_bareground(ngrid)
86	REAL,INTENT(IN) :: albedo_co2_ice_SPECTV(L_NSPECTV)
87	REAL,INTENT(INOUT) :: albedo(ngrid,L_NSPECTV)
88	REAL,INTENT(OUT) :: pemisurf(ngrid)
89	REAL,INTENT(OUT) :: pdtc(ngrid,nlayer)
90	REAL,INTENT(OUT) :: pdtsrfc(ngrid)
91	REAL,INTENT(OUT) :: pdpsrfc(ngrid)
92	REAL,INTENT(OUT) :: pdqc(ngrid,nlayer,nq)
93	REAL,INTENT(OUT) :: pdqsurfc(ngrid)
94
95	!------------------------------
96	! Local variables
97	!------------------------------
98
99	INTEGER l,ig,icap,ilay,iq,nw,igas,it
100
101	REAL reffrad(ngrid,nlayer) ! Radius (m) of the CO2 ice particles.
102	REAL*8 zt(ngrid,nlayer) ! Updated Atmospheric Temperatures (K).
103	REAL ztsrf(ngrid) ! Updated Surface Temperatures (K).
104	REAL zq(ngrid,nlayer,nq) ! Updated Atmospheric tracers mixing ratios (kg/kg of air).
105	REAL piceco2(ngrid) ! Updated Surface Tracer (kg/m2).
106	REAL ztcond (ngrid,nlayer) ! Atmospheric Temperatures of condensation of CO2.
107	REAL ztnuc (ngrid,nlayer) ! Atmospheric Nucleation Temperatures.
108	REAL ztcondsol(ngrid) ! Temperatures of condensation of CO2 at the surface.
109	REAL zcondices(ngrid) ! Condensation rate on the ground (kg/m2/s).
110	REAL zfallice(ngrid) ! Flux of ice falling on the surface (kg/m2/s).
111	REAL Mfallice(ngrid) ! Total amount of ice fallen to the ground during the timestep (kg/m2).
112	REAL wq(ngrid,nlayer+1) ! Total amount of ice fallen to the ground during the timestep (kg/m2).
113	REAL subptimestep ! Duration of the subtimestep (s) for the sedimentation.
114	Integer Ntime ! Number of subtimesteps.
115	REAL masse (ngrid,nlayer) ! Mass of atmospheric layers (kg/m2)
116	REAL w(ngrid,nlayer,nq) !
117	REAL vstokes,reff !
118	REAL ppco2 !
119
120
121	!------------------------------------------
122	! Saved local variables
123	!------------------------------------------
124
125
126	REAL,SAVE :: latcond=5.9e5
127	REAL,SAVE :: ccond
128	REAL,SAVE,ALLOCATABLE,DIMENSION(:) :: emisref
129	!$OMP THREADPRIVATE(latcond,ccond,emisref)
130	LOGICAL,SAVE :: firstcall=.true.
131	!$OMP THREADPRIVATE(firstcall)
132	INTEGER,SAVE :: i_co2ice=0 ! co2 ice
133	!$OMP THREADPRIVATE(i_co2ice)
134	CHARACTER(LEN=20) :: tracername ! to temporarily store text
135
136
137	!------------------------------------------------
138	! Initialization at the first call
139	!------------------------------------------------
140
141
142	IF (firstcall) THEN
143
144	ALLOCATE(emisref(ngrid))
145	! Find CO2 ice tracer.
146	do iq=1,nq
147	tracername=noms(iq)
148	if (tracername.eq."co2_ice") then
149	i_co2ice=iq
150	endif
151	enddo
152
153	write(,) "condense_co2: i_co2ice=",i_co2ice
154
155	if((i_co2ice.lt.1))then
156	print*,'In condens_cloud but no CO2 ice tracer, exiting.'
157	print*,'Still need generalisation to arbitrary species!'
158	stop
159	endif
160
161	ccond=cpp/(g*latcond)
162	print*,'In condens_cloud: ccond=',ccond,' latcond=',latcond
163
164	! Prepare special treatment if gas is not pure CO2
165	! if (addn2) then
166	! m_co2 = 44.01E-3 ! CO2 molecular mass (kg/mol)
167	! m_noco2 = 28.02E-3 ! N2 molecular mass (kg/mol)
168	! Compute A and B coefficient use to compute
169	! mean molecular mass Mair defined by
170	! 1/Mair = q(ico2)/m_co2 + (1-q(ico2))/m_noco2
171	! 1/Mair = A*q(ico2) + B
172	! A = (1/m_co2 - 1/m_noco2)
173	! B = 1/m_noco2
174	! endif
175
176	! Minimum CO2 mixing ratio below which mixing occurs with layer above : qco2min =0.75
177
178	firstcall=.false.
179	ENDIF
180
181
182	!------------------------------------------------
183	! Tendencies initially set to 0
184	!------------------------------------------------
185
186
187	pdqc(1:ngrid,1:nlayer,1:nq) = 0.
188	pdtc(1:ngrid,1:nlayer) = 0.
189	zq(1:ngrid,1:nlayer,1:nq) = 0.
190	zt(1:ngrid,1:nlayer) = 0.
191	Mfallice(1:ngrid) = 0.
192	zfallice(1:ngrid) = 0.
193	zcondices(1:ngrid) = 0.
194	pdtsrfc(1:ngrid) = 0.
195	pdpsrfc(1:ngrid) = 0.
196	pdqsurfc(1:ngrid) = 0.
197
198
199	!----------------------------------
200	! Atmospheric condensation
201	!----------------------------------
202
203
204	! Compute CO2 Volume mixing ratio
205	! -------------------------------
206	! if (addn2) then
207	! DO l=1,nlayer
208	! DO ig=1,ngrid
209	! qco2=pq(ig,l,ico2)+pdq(ig,l,ico2)*ptimestep
210	! Mean air molecular mass = 1/(q(ico2)/m_co2 + (1-q(ico2))/m_noco2)
211	! mmean=1/(A*qco2 +B)
212	! vmr_co2(ig,l) = qco2*mmean/m_co2
213	! ENDDO
214	! ENDDO
215	! else
216	! DO l=1,nlayer
217	! DO ig=1,ngrid
218	! vmr_co2(ig,l)=0.5
219	! ENDDO
220	! ENDDO
221	! end if
222
223
224	! Forecast the atmospheric frost temperature 'ztcond' and nucleation temperature 'ztnuc'.
225	DO l=1,nlayer
226	DO ig=1,ngrid
227	ppco2=gfrac(igas_CO2)*pplay(ig,l)
228	call get_tcond_co2(ppco2,ztcond(ig,l))
229	call get_tnuc_co2(ppco2,ztnuc(ig,l))
230	ENDDO
231	ENDDO
232
233	! Initialize zq and zt at the beginning of the sub-timestep loop and qsurf.
234	DO ig=1,ngrid
235	piceco2(ig)=pqsurf(ig,i_co2ice)
236	DO l=1,nlayer
237	zt(ig,l)=pt(ig,l)
238	zq(ig,l,i_co2ice)=pq(ig,l,i_co2ice)
239	IF( zq(ig,l,i_co2ice).lt.-1.e-6 ) THEN
240	print*,'Uh-oh, zq = ',zq(ig,l,i_co2ice),'at ig,l=',ig,l
241	if(l.eq.1)then
242	print*,'Perhaps the atmosphere is collapsing on surface...?'
243	endif
244	END IF
245	ENDDO
246	ENDDO
247
248	! Calculate the mass of each atmospheric layer (kg.m-2)
249	do ilay=1,nlayer
250	DO ig=1,ngrid
251	masse(ig,ilay)=(pplev(ig,ilay) - pplev(ig,ilay+1)) /g
252	end do
253	end do
254
255
256	!-----------------------------------------------------------
257	! START CONDENSATION/SEDIMENTATION SUB-TIME LOOP
258	!-----------------------------------------------------------
259
260
261	Ntime = 20 ! number of sub-timestep
262	subptimestep = ptimestep/float(Ntime)
263
264	! Add the tendencies from other physical processes at each subtimstep.
265	DO it=1,Ntime
266	DO l=1,nlayer
267	DO ig=1,ngrid
268	zt(ig,l) = zt(ig,l) + pdt(ig,l) * subptimestep
269	zq(ig,l,i_co2ice) = zq(ig,l,i_co2ice) + pdq(ig,l,i_co2ice) * subptimestep
270	END DO
271	END DO
272
273	! Gravitational sedimentation starts.
274
275	! Sedimentation computed from radius computed from q in module radii_mod.
276	call co2_reffrad(ngrid,nlayer,nq,zq,reffrad)
277
278	DO ilay=1,nlayer
279	DO ig=1,ngrid
280
281	reff = reffrad(ig,ilay)
282
283	call stokes &
284	(pplev(ig,ilay),pt(ig,ilay), &
285	reff,vstokes,rho_co2)
286
287	!w(ig,ilay,i_co2ice) = 0.0
288	w(ig,ilay,i_co2ice) = vstokes * subptimestep * &
289	pplev(ig,ilay)/(r*pt(ig,ilay))
290
291	END DO
292	END DO
293
294	! Computing q after sedimentation
295	call vlz_fi(ngrid,nlayer,zq(1,1,i_co2ice),2.,masse,w(1,1,i_co2ice),wq)
296
297
298	! Progressively accumulating the flux to the ground.
299	! Mfallice is the total amount of ice fallen to the ground.
300	DO ig=1,ngrid
301	Mfallice(ig) = Mfallice(ig) + wq(ig,i_co2ice)
302	END DO
303
304	!----------------------------------------------------------
305	! Condensation / sublimation in the atmosphere
306	!----------------------------------------------------------
307	! (MODIFICATIONS FOR EARLY MARS: falling heat neglected, condensation of CO2 into tracer i_co2ice)
308
309
310	DO l=nlayer , 1, -1
311	DO ig=1,ngrid
312	pdtc(ig,l)=0.
313
314	! ztcond-> ztnuc in test beneath to nucleate only when super saturation occurs(JL 2011)
315	IF ((zt(ig,l).LT.ztnuc(ig,l)).or.(zq(ig,l,i_co2ice).gt.1.E-10)) THEN
316	pdtc(ig,l) = (ztcond(ig,l) - zt(ig,l))/subptimestep
317	pdqc(ig,l,i_co2ice) = pdtc(ig,l)ccondg
318
319	! Case when the ice from above sublimes entirely
320	IF ((zq(ig,l,i_co2ice).lt.-pdqc(ig,l,i_co2ice)*subptimestep) &
321	.AND. (zq(ig,l,i_co2ice).gt.0)) THEN
322
323	pdqc(ig,l,i_co2ice) = -zq(ig,l,i_co2ice)/subptimestep
324	pdtc(ig,l) =-zq(ig,l,i_co2ice)/(ccondgsubptimestep)
325
326	END IF
327
328	! Temperature and q after condensation
329	zt(ig,l) = zt(ig,l) + pdtc(ig,l) * subptimestep
330	zq(ig,l,i_co2ice) = zq(ig,l,i_co2ice) + pdqc(ig,l,i_co2ice) * subptimestep
331	END IF
332
333	ENDDO
334	ENDDO
335
336	ENDDO! end of subtimestep loop.
337
338	! Computing global tendencies after the subtimestep.
339	DO l=1,nlayer
340	DO ig=1,ngrid
341	pdtc(ig,l) = &
342	(zt(ig,l) - (pt(ig,l) + pdt(ig,l)*ptimestep))/ptimestep
343	pdqc(ig,l,i_co2ice) = &
344	(zq(ig,l,i_co2ice)-(pq(ig,l,i_co2ice)+pdq(ig,l,i_co2ice)*ptimestep))/ptimestep
345	END DO
346	END DO
347	DO ig=1,ngrid
348	zfallice(ig) = Mfallice(ig)/ptimestep
349	END DO
350
351
352	!-----------------------------------------------------------------------
353	! Condensation/sublimation on the ground
354	!-----------------------------------------------------------------------
355
356
357	! Forecast of ground temperature ztsrf and frost temperature ztcondsol.
358	DO ig=1,ngrid
359	ppco2=gfrac(igas_CO2)*pplay(ig,1)
360	call get_tcond_co2(ppco2,ztcondsol(ig))
361
362	ztsrf(ig) = ptsrf(ig)
363
364	if((ztsrf(ig).le.ztcondsol(ig)+2.0).and.(ngrid.eq.1))then
365	print*,'CO2 is condensing on the surface in 1D. This atmosphere is doomed.'
366	print*,'T_surf = ',ztsrf,'K'
367	print*,'T_cond = ',ztcondsol,'K'
368	open(116,file='surf_vals.out')
369	write(116,*) 0.0, pplev(1,1), 0.0, 0.0
370	close(116)
371	call abort
372	endif
373
374	ztsrf(ig) = ptsrf(ig) + pdtsrf(ig)*ptimestep
375
376	ENDDO
377
378	DO ig=1,ngrid
379
380	IF(ig.GT.ngrid/2+1) THEN
381	icap=2
382	ELSE
383	icap=1
384	ENDIF
385
386	! Loop over where we have condensation / sublimation
387	IF ((ztsrf(ig) .LT. ztcondsol(ig)) .OR. & ! ground condensation
388	(zfallice(ig).NE.0.) .OR. & ! falling snow
389	((ztsrf(ig) .GT. ztcondsol(ig)) .AND. & ! ground sublimation
390	((piceco2(ig)+zfallice(ig)*ptimestep) .NE. 0.))) THEN
391
392
393	! Condensation or partial sublimation of CO2 ice
394	zcondices(ig)=pcapcal(ig)*(ztcondsol(ig)-ztsrf(ig)) &
395	/(latcond*ptimestep)
396	pdtsrfc(ig) = (ztcondsol(ig) - ztsrf(ig))/ptimestep
397
398	! If the entire CO_2 ice layer sublimes
399	! (including what has just condensed in the atmosphere)
400	IF((piceco2(ig)/ptimestep+zfallice(ig)).LE. &
401	-zcondices(ig))THEN
402	zcondices(ig) = -piceco2(ig)/ptimestep - zfallice(ig)
403	pdtsrfc(ig)=(latcond/pcapcal(ig))* &
404	(zcondices(ig))
405	END IF
406
407	! Changing CO2 ice amount and pressure
408	piceco2(ig) = piceco2(ig) + pdqsurfc(ig)*ptimestep
409	pdqsurfc(ig) = zcondices(ig) + zfallice(ig)
410	pdpsrfc(ig) = -pdqsurfc(ig)*g
411
412	IF(ABS(pdpsrfc(ig)*ptimestep).GT.pplev(ig,1)) THEN
413	PRINT*,'STOP in condens in condense_co2'
414	PRINT*,'condensing more than total mass'
415	PRINT*,'Grid point ',ig
416	PRINT*,'Ps = ',pplev(ig,1)
417	PRINT*,'d Ps = ',pdpsrfc(ig)
418	STOP
419	ENDIF
420	END IF
421
422	ENDDO ! end of ngrid loop.
423
424
425	!---------------------------------------------------------------------------------------------
426	! Surface albedo and emissivity of the ground below the snow (emisref)
427	!---------------------------------------------------------------------------------------------
428
429
430	DO ig=1,ngrid
431
432	IF(latitude(ig).LT.0.) THEN
433	icap=2 ! Southern Hemisphere
434	ELSE
435	icap=1 ! Nortnern hemisphere
436	ENDIF
437
438	if(.not.piceco2(ig).ge.0.) THEN
439	if(piceco2(ig).le.-1.e-8) print*, &
440	'WARNING : in condense_co2cloud: piceco2(',ig,')=', piceco2(ig)
441	piceco2(ig)=0.
442	endif
443	if (piceco2(ig) .gt. 1.) then ! CO2 Albedo condition changed to ~1 mm coverage. Change by MT2015.
444	DO nw=1,L_NSPECTV
445	albedo(ig,nw) = albedo_co2_ice_SPECTV(nw)
446	ENDDO
447	emisref(ig) = emisice(icap)
448	else
449	DO nw=1,L_NSPECTV
450	albedo(ig,nw) = albedo_bareground(ig) ! Note : If you have some water, it will be taken into account in the "hydrol" routine.
451	ENDDO
452	emisref(ig) = emissiv
453	pemisurf(ig) = emissiv
454	end if
455
456	END DO
457
458	return
459
460	end subroutine condense_co2
461
462
463
464	!-------------------------------------------------------------------------
465	!-------------------------------------------------------------------------
466	!-------------------------------------------------------------------------
467	!-------------------------------------------------------------------------
468	!-------------------------------------------------------------------------
469	!-------------------------------------------------------------------------
470
471
472
473	subroutine get_tcond_co2(p,tcond) ! Calculates the condensation temperature for CO2
474
475
476	implicit none
477
478	real p, peff, tcond
479	real, parameter :: ptriple=518000.0
480
481	peff=p
482
483	if(peff.lt.ptriple) then
484	tcond = (-3167.8)/(log(.01*peff)-23.23) ! Fanale's formula.
485	else
486	tcond = 684.2-92.3log(peff)+4.32log(peff)**2 ! liquid-vapour transition (based on CRC handbook 2003 data)
487	endif
488	return
489
490	end subroutine get_tcond_co2
491
492
493
494	!-------------------------------------------------------------------------
495	!-------------------------------------------------------------------------
496	!-------------------------------------------------------------------------
497	!-------------------------------------------------------------------------
498	!-------------------------------------------------------------------------
499	!-------------------------------------------------------------------------
500
501
502
503	subroutine get_tnuc_co2(p,tnuc)
504	! Calculates the nucleation temperature for CO2, based on a simple super saturation criterion. JL 2011.
505
506	use callkeys_mod, only: co2supsat
507
508	implicit none
509
510	real p, peff, tnuc
511	real, parameter :: ptriple=518000.0
512
513	peff=p/co2supsat
514
515	if(peff.lt.ptriple) then
516	tnuc = (-3167.8)/(log(.01*peff)-23.23) ! Fanale's formula
517	else
518	tnuc = 684.2-92.3log(peff)+4.32log(peff)**2
519	! liquid-vapour transition (based on CRC handbook 2003 data)
520	endif
521
522	return
523
524	end subroutine get_tnuc_co2

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

Download in other formats:

Original Format