Context Navigation

condense_co2.F90 @ 1484

Last change on this file since 1484 was 1484, checked in by mturbet, 9 years ago
Harmonization of surface tracer tendencies calculation
Property svn:executable set to ``*
File size: 17.7 KB

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

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

Download in other formats:

Original Format