Context Navigation

condens_co2cloud.F90 @ 136

Last change on this file since 136 was 135, checked in by aslmd, 14 years ago
CHANGEMENT ARBORESCENCE ETAPE 2 -- NON COMPLET
File size: 14.7 KB

Line
1	subroutine condens_co2cloud(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,reffrad,cpp3D)
7
8	use radinc_h, only : naerkind
9
10	implicit none
11
12	!==================================================================
13	! Purpose
14	! -------
15	! Condense and/or sublime CO2 ice on the ground and in the
16	! atmosphere, and sediment the ice.
17	!
18	! Inputs
19	! ------
20	! ngrid Number of vertical columns
21	! nlayer Number of layers
22	! pplay(ngrid,nlayer) Pressure layers
23	! pplev(ngrid,nlayer+1) Pressure levels
24	! pt(ngrid,nlayer) Temperature (in K)
25	! ptsrf(ngrid) Surface temperature
26	!
27	! pdt(ngrid,nlayermx) Time derivative before condensation/sublimation of pt
28	! pdtsrf(ngrid) Time derivative before condensation/sublimation of ptsrf
29	! pqsurf(ngrid,nq) Sedimentation flux at the surface (kg.m-2.s-1)
30	!
31	! Outputs
32	! -------
33	! pdpsrf(ngrid) \ Contribution of condensation/sublimation
34	! pdtc(ngrid,nlayermx) / to the time derivatives of Ps, pt, and ptsrf
35	! pdtsrfc(ngrid) /
36	!
37	! Both
38	! ----
39	! piceco2(ngrid) CO2 ice at the surface (kg/m2)
40	! psolaralb(ngrid) Albedo at the surface
41	! pemisurf(ngrid) Emissivity of the surface
42	!
43	! Authors
44	! -------
45	! Francois Forget (1996)
46	! Converted to Fortran 90 and slightly modified by R. Wordsworth (2009)
47	!
48	!==================================================================
49
50	#include "dimensions.h"
51	#include "dimphys.h"
52	#include "comcstfi.h"
53	#include "surfdat.h"
54	#include "comgeomfi.h"
55	#include "comvert.h"
56	#include "callkeys.h"
57	#include "tracer.h"
58
59	!-----------------------------------------------------------------------
60	! Arguments
61
62	INTEGER ngrid, nlayer, nq
63
64	REAL ptimestep
65	REAL pplay(ngrid,nlayer),pplev(ngrid,nlayer+1)
66	REAL pcapcal(ngrid)
67	REAL pt(ngrid,nlayer)
68	REAL ptsrf(ngrid)
69	REAL pphi(ngrid,nlayer)
70	REAL pdt(ngrid,nlayer),pdtsrf(ngrid),pdtc(ngrid,nlayer)
71	REAL pdtsrfc(ngrid),pdpsrf(ngrid)
72	! REAL piceco2(ngrid),psolaralb(ngrid,2),pemisurf(ngrid)
73	REAL piceco2(ngrid),psolaralb(ngrid),pemisurf(ngrid)
74
75	REAL pu(ngrid,nlayer) , pv(ngrid,nlayer)
76	REAL pdu(ngrid,nlayer) , pdv(ngrid,nlayer)
77	REAL pduc(ngrid,nlayer) , pdvc(ngrid,nlayer)
78	REAL pq(ngridmx,nlayer,nq),pdq(ngrid,nlayer,nq)
79	REAL pdqc(ngrid,nlayer,nq), pdqsc(ngrid,nq)
80
81	REAL reffrad(ngrid,nlayer,naerkind)
82
83	! Allow variations in cp with location
84	REAL cpp3D(ngridmx,nlayermx) ! specific heat capacity at const. pressure
85
86
87
88	!-----------------------------------------------------------------------
89	! Local variables
90
91	INTEGER l,ig,icap,ilay,it,iq
92
93	REAL*8 zt(ngridmx,nlayermx)
94	REAL zq(ngridmx,nlayermx,nqmx)
95	REAL zcpi
96	REAL ztcond (ngridmx,nlayermx)
97	REAL ztcondsol(ngridmx)
98	REAL zdiceco2(ngridmx)
99	REAL zcondicea(ngridmx,nlayermx), zcondices(ngridmx)
100	REAL zfallice(ngridmx), Mfallice(ngridmx)
101	REAL zmflux(nlayermx+1)
102	REAL zu(nlayermx),zv(nlayermx)
103	REAL ztsrf(ngridmx)
104	REAL ztc(nlayermx), ztm(nlayermx+1)
105	REAL zum(nlayermx+1) , zvm(nlayermx+1)
106	LOGICAL condsub(ngridmx)
107	REAL subptimestep
108	Integer Ntime
109	real masse (ngridmx,nlayermx), w(ngridmx,nlayermx,nqmx)
110	real wq(ngridmx,nlayermx+1)
111	real vstokes,reff
112
113	! Special diagnostic variables
114	real tconda1(ngridmx,nlayermx)
115	real tconda2(ngridmx,nlayermx)
116	real zdtsig (ngridmx,nlayermx)
117	real zdt (ngridmx,nlayermx)
118
119	!-----------------------------------------------------------------------
120	! Saved local variables
121
122	REAL emisref(ngridmx)
123	REAL latcond
124	REAL ccond
125	REAL cpice
126	SAVE emisref, cpice
127	SAVE latcond,ccond
128
129	LOGICAL firstcall
130	SAVE firstcall
131	REAL SSUM
132	EXTERNAL SSUM
133
134	DATA latcond /5.9e5/
135	DATA cpice /1000./
136	DATA firstcall/.true./
137
138	REAL CBRT
139	EXTERNAL CBRT
140
141	INTEGER,SAVE :: i_co2ice=0 ! co2 ice
142	CHARACTER(LEN=20) :: tracername ! to temporarily store text
143
144	real reffradmin, reffradmax
145
146	!-----------------------------------------------------------------------
147	! Initializations
148
149	call zerophys(ngridnlayernq,pdqc)
150	call zerophys(ngridnlayernq,pdtc)
151	call zerophys(ngridmxnlayermxnqmx,zq)
152	call zerophys(ngridmx*nlayermx,zt)
153
154	!reffradmin=1.e-7
155	!reffradmax=5.e-4
156	reffradmin=0.5e-7
157	reffradmax=0.1e-3 ! FF data
158	!reffradmin=0.1e-5
159	!reffradmax=0.5e-3 ! JB data
160	! improve this in the future...
161
162	! Initialisation
163	IF (firstcall) THEN
164	do iq=1,nqmx
165	tracername=noms(iq)
166	if (tracername.eq."co2_ice") then
167	i_co2ice=iq
168	endif
169	enddo
170
171	write(,) "condense_co2cloud: i_co2ice=",i_co2ice
172
173	if((i_co2ice.lt.1))then
174	print*,'In condens_co2cloud but no CO2 ice tracer, exiting.'
175	stop
176	endif
177
178	ccond=cpp/(g*latcond)
179	print*,'In condens_co2cloud: ccond=',ccond,' latcond=',latcond
180
181	firstcall=.false.
182	ENDIF
183	zcpi=1./cpp
184
185	!-----------------------------------------------------------------------
186	! Calculation of CO2 condensation / sublimation
187	!
188	! Variables used:
189	! piceco2(ngrid) amount of co2 ice on the ground (kg/m2)
190	! zcondicea(ngrid,l) condensation rate in layer l (kg/m2/s)
191	! zcondices(ngrid) condensation rate on the ground (kg/m2/s)
192	! zfallice(ngrid) flux of ice falling on surface (kg/m2/s)
193	! pdtc(ngrid,nlayermx) dT/dt due to phase changes (K/s)
194
195
196	! Tendencies initially set to 0 (except pdtc)
197	DO l=1,nlayer
198	DO ig=1,ngrid
199	zcondicea(ig,l) = 0.
200	pduc(ig,l) = 0
201	pdvc(ig,l) = 0
202	pdqc(ig,l,i_co2ice) = 0
203	END DO
204	END DO
205
206	DO ig=1,ngrid
207	Mfallice(ig) = 0.
208	zfallice(ig) = 0.
209	zcondices(ig) = 0.
210	pdtsrfc(ig) = 0.
211	pdpsrf(ig) = 0.
212	condsub(ig) = .false.
213	zdiceco2(ig) = 0.
214	ENDDO
215
216	!-----------------------------------------------------------------------
217	! Atmospheric condensation
218
219	! Forecast the atmospheric frost temperature 'ztcond'
220	DO l=1,nlayer
221	DO ig=1,ngrid
222	call get_tcond_co2(pplay(ig,l),ztcond(ig,l),1,1)
223	ENDDO
224	ENDDO
225
226	! Initialize zq and zt at the beginning of the sub-timestep loop
227	DO l=1,nlayer
228	DO ig=1,ngrid
229	zt(ig,l)=pt(ig,l)
230	zq(ig,l,i_co2ice)=pq(ig,l,i_co2ice)
231	IF( zq(ig,l,i_co2ice).lt.-1.e-6 ) THEN
232	print*,'Uh-oh, zq = ',zq(ig,l,i_co2ice),'at ig,l=',ig,l
233	if(l.eq.1)then
234	print*,'Perhaps the atmosphere is collapsing on surface...?'
235	endif
236	stop
237
238	END IF
239	ENDDO
240	ENDDO
241
242	! Calculate the mass of each atmospheric layer (kg.m-2)
243	do ilay=1,nlayer
244	do ig=1, ngrid
245	masse(ig,ilay)=(pplev(ig,ilay) - pplev(ig,ilay+1)) /g
246	end do
247	end do
248
249	! -----------------------------------------------
250	! START CONDENSATION/SEDIMENTATION SUB-TIME LOOP
251	! -----------------------------------------------
252	Ntime = 20 ! number of sub-timestep
253	subptimestep = ptimestep/float(Ntime)
254
255	DO it=1,Ntime
256
257	! Add the tendencies from other physical processes at each subtimstep
258	DO l=1,nlayer
259	DO ig=1,ngrid
260	zt(ig,l) = zt(ig,l) + pdt(ig,l) * subptimestep
261	zq(ig,l,i_co2ice) = zq(ig,l,i_co2ice) + pdq(ig,l,i_co2ice) * subptimestep
262	END DO
263	END DO
264
265
266	! Gravitational sedimentation
267
268	! sedimentation computed from radius computed from q
269	! assuming that the ice is splitted in Nmix particle
270	do ilay=1,nlayer
271	do ig=1, ngrid
272
273	reff = CBRT( 3zq(ig,ilay,i_co2ice)/( 4Nmix_co2pirho_co2 ))
274
275	! there should be a more elegant way of doing this...
276	if(reff.lt.1.e-16) reff=1.e-16
277
278
279	! update reffrad for radiative transfer
280	if(reff.lt.reffradmin)then
281	reffrad(ig,ilay,1)=reffradmin
282	!print*,'reff below optical limit'
283	elseif(reff.gt.reffradmax)then
284	reffrad(ig,ilay,1)=reffradmax
285	!print*,'reff above optical limit'
286	else
287	reffrad(ig,ilay,1)=reff
288	endif
289
290	call stokes &
291	(pplev(ig,ilay),pt(ig,ilay), &
292	reff,vstokes,rho_co2)
293
294	w(ig,ilay,i_co2ice) = vstokes * subptimestep * &
295	pplev(ig,ilay)/(r*pt(ig,ilay))
296
297	end do
298	end do
299
300	! Computing q after sedimentation
301
302	call vlz_fi(ngrid,zq(1,1,i_co2ice),2.,masse,w(1,1,i_co2ice),wq)
303
304
305	! Progressively accumulating the flux to the ground
306	! Mfallice is the total amount of ice fallen to the ground
307	do ig=1,ngrid
308	Mfallice(ig) = Mfallice(ig) + wq(ig,i_co2ice)
309	end do
310
311
312	! Condensation / sublimation in the atmosphere
313	! --------------------------------------------
314	! (calculation of zcondicea, zfallice and pdtc)
315	! (MODIFICATIONS FOR EARLY MARS: falling heat neglected, condensation
316	! of CO2 into tracer i_co2ice)
317
318	DO l=nlayer , 1, -1
319	DO ig=1,ngrid
320	pdtc(ig,l)=0.
321
322	! tweak ccond if the gas is non-ideal
323	if(nonideal)then
324	ccond=cpp3D(ig,l)/(g*latcond)
325	endif
326
327
328	IF ((zt(ig,l).LT.ztcond(ig,l)).or.(zq(ig,l,i_co2ice).gt.0)) THEN
329	condsub(ig)=.true.
330	pdtc(ig,l) = (ztcond(ig,l) - zt(ig,l))/subptimestep
331	pdqc(ig,l,i_co2ice) = pdtc(ig,l)ccondg
332
333	! Case when the ice from above sublimes entirely
334	IF ((zq(ig,l,i_co2ice).lt.-pdqc(ig,l,i_co2ice)*subptimestep) &
335	.AND. (zq(ig,l,i_co2ice).gt.0)) THEN
336
337	pdqc(ig,l,i_co2ice) = -zq(ig,l,i_co2ice)/subptimestep
338	pdtc(ig,l) =-zq(ig,l,i_co2ice)/(ccondgsubptimestep)
339
340	END IF
341
342	! Temperature and q after condensation
343	zt(ig,l) = zt(ig,l) + pdtc(ig,l) * subptimestep
344	zq(ig,l,i_co2ice) = zq(ig,l,i_co2ice) + pdqc(ig,l,i_co2ice) * subptimestep
345	END IF
346
347	ENDDO
348	ENDDO
349	ENDDO ! end of subtimestep loop
350
351	! Computing global tendencies after the subtimestep
352	DO l=1,nlayer
353	DO ig=1,ngrid
354	pdtc(ig,l) = &
355	(zt(ig,l) - (pt(ig,l) + pdt(ig,l)*ptimestep))/ptimestep
356	pdqc(ig,l,i_co2ice) = &
357	(zq(ig,l,i_co2ice)-(pq(ig,l,i_co2ice)+pdq(ig,l,i_co2ice)*ptimestep))/ptimestep
358	END DO
359	END DO
360	DO ig=1,ngrid
361	zfallice(ig) = Mfallice(ig)/ptimestep
362	END DO
363
364
365	!-----------------------------------------------------------------------
366	! Condensation/sublimation on the ground
367	! (calculation of zcondices and pdtsrfc)
368
369	! forecast of ground temperature ztsrf and frost temperature ztcondsol
370	DO ig=1,ngrid
371
372	call get_tcond_co2(pplev(ig,1),ztcondsol(ig))
373
374	ztsrf(ig) = ptsrf(ig)
375
376	if((ztsrf(ig).le.ztcondsol(ig)+2.0).and.(ngrid.eq.1))then
377	print*,'CO2 is condensing on the surface in 1D. This atmosphere is doomed.'
378	print*,'T_surf = ',ztsrf,'K'
379	print*,'T_cond = ',ztcondsol,'K'
380	open(116,file='surf_vals.out')
381	write(116,*) 0.0, pplev(1,1), 0.0, 0.0
382	close(116)
383	call abort
384	endif
385
386	ztsrf(ig) = ptsrf(ig) + pdtsrf(ig)*ptimestep
387
388	ENDDO
389
390	DO ig=1,ngrid
391	IF(ig.GT.ngrid/2+1) THEN
392	icap=2
393	ELSE
394	icap=1
395	ENDIF
396
397	! Loop over where we have condensation / sublimation
398	IF ((ztsrf(ig) .LT. ztcondsol(ig)) .OR. & ! ground condensation
399	(zfallice(ig).NE.0.) .OR. & ! falling snow
400	((ztsrf(ig) .GT. ztcondsol(ig)) .AND. & ! ground sublimation
401	((piceco2(ig)+zfallice(ig)*ptimestep) .NE. 0.))) THEN
402	condsub(ig) = .true.
403
404	! Condensation or partial sublimation of CO2 ice
405	zcondices(ig)=pcapcal(ig)*(ztcondsol(ig)-ztsrf(ig)) &
406	/(latcond*ptimestep)
407	pdtsrfc(ig) = (ztcondsol(ig) - ztsrf(ig))/ptimestep
408
409	! If the entire CO_2 ice layer sublimes
410	! (including what has just condensed in the atmosphere)
411	IF((piceco2(ig)/ptimestep+zfallice(ig)).LE. &
412	-zcondices(ig))THEN
413	zcondices(ig) = -piceco2(ig)/ptimestep - zfallice(ig)
414	pdtsrfc(ig)=(latcond/pcapcal(ig))* &
415	(zcondices(ig))
416	END IF
417
418	! Changing CO2 ice amount and pressure
419
420	zdiceco2(ig) = zcondices(ig) + zfallice(ig)
421	piceco2(ig) = piceco2(ig) + zdiceco2(ig)*ptimestep
422	pdpsrf(ig) = -zdiceco2(ig)*g
423
424	IF(ABS(pdpsrf(ig)*ptimestep).GT.pplev(ig,1)) THEN
425	PRINT*,'STOP in condens'
426	PRINT*,'condensing more than total mass'
427	PRINT*,'Grid point ',ig
428	PRINT*,'Ps = ',pplev(ig,1)
429	PRINT*,'d Ps = ',pdpsrf(ig)
430	STOP
431	ENDIF
432	END IF
433	ENDDO
434
435	! Surface albedo and emissivity of the ground below the snow (emisref)
436	! --------------------------------------------------------------------
437	do ig =1,ngrid
438	IF(ig.GT.ngrid/2+1) THEN
439	icap=2
440	ELSE
441	icap=1
442	ENDIF
443
444	if(.not.piceco2(ig).ge.0.) THEN
445	if(piceco2(ig).le.-1.e-8) print*, &
446	'WARNING in condense_co2cloud: piceco2(',ig,')=', piceco2(ig)
447	piceco2(ig)=0.
448	endif
449	if (piceco2(ig).gt.0) then
450	psolaralb(ig) = albedice(icap)
451	emisref(ig) = emisice(icap)
452	else
453	psolaralb(ig) = albedodat(ig)
454	emisref(ig) = emissiv
455	pemisurf(ig) = emissiv
456	end if
457	end do
458
459	return
460	end subroutine condens_co2cloud
461
462	!-------------------------------------------------------------------------
463	subroutine get_tcond_co2(p,tcond)
464	! Calculates the condensation temperature for CO2
465
466	implicit none
467
468	real p, tcond
469	tcond = (-3167.8)/(log(.01*p)-23.23) ! Fanale's formula
470	return
471
472	end subroutine get_tcond_co2

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

Download in other formats:

Original Format