Context Navigation

condense_co2.F90 @ 1477

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

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

Download in other formats:

Original Format