Context Navigation

condense_cloud.F90 @ 671

Last change on this file since 671 was 586, checked in by jleconte, 13 years ago
removed cpp3D and nonideal stuff.
Property svn:executable set to ``*
File size: 17.3 KB

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

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

Download in other formats:

Original Format