Context Navigation

condense_cloud.F90 @ 1145

Last change on this file since 1145 was 869, checked in by aslmd, 12 years ago

24/01/2013 == AS + JL

A more robust way to refer to gas type.

Gas names with an arbitrary number of characters (<20) can be used This is good for C2H2, C2H6, H2SO4, C17H21NO4, etc... !!! Remember this must be compliant with Q.dat in corrk_data !!!
igas_... labels are assigned once for all in su_gases Then using igas_... everywhere instead of gnom (except for kcm stuff)
Users can still use e.g. H2_ but H2 also works
Simplified condense_cloud so that igas_CO2 is used directly

Property svn:executable set to *

File size: 16.7 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)
7
8	use radinc_h, only : naerkind
9	use gases_h
10	use radii_mod, only : co2_reffrad
11	use aerosol_mod, only : iaero_co2
12	USE surfdat_h
13	USE comgeomfi_h
14	USE tracer_h
15
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,nlayermx) 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,nlayermx) / 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	#include "dimensions.h"
59	#include "dimphys.h"
60	#include "comcstfi.h"
61	#include "comvert.h"
62	#include "callkeys.h"
63
64	!-----------------------------------------------------------------------
65	! Arguments
66
67	INTEGER,INTENT(IN) :: ngrid
68	INTEGER,INTENT(IN) :: nlayer
69	INTEGER,INTENT(IN) :: nq
70	REAL,INTENT(IN) :: ptimestep
71	REAL,INTENT(IN) :: pcapcal(ngrid)
72	REAL,INTENT(IN) :: pplay(ngrid,nlayer)
73	REAL,INTENT(IN) :: pplev(ngrid,nlayer+1)
74	REAL,INTENT(IN) :: ptsrf(ngrid)
75	REAL,INTENT(IN) :: pt(ngrid,nlayer)
76	REAL,INTENT(IN) :: pphi(ngrid,nlayer)
77	REAL,INTENT(IN) :: pdt(ngrid,nlayer)
78	REAL,INTENT(IN) :: pdu(ngrid,nlayer)
79	REAL,INTENT(IN) :: pdv(ngrid,nlayer)
80	REAL,INTENT(IN) :: pdtsrf(ngrid)
81	REAL,INTENT(IN) :: pu(ngrid,nlayer)
82	REAL,INTENT(IN) :: pv(ngrid,nlayer)
83	REAL,INTENT(IN) :: pq(ngrid,nlayer,nq)
84	REAL,INTENT(IN) :: pdq(ngrid,nlayer,nq)
85	REAL,INTENT(INOUT) :: piceco2(ngrid)
86	REAL,INTENT(OUT) :: psolaralb(ngrid)
87	REAL,INTENT(OUT) :: pemisurf(ngrid)
88	REAL,INTENT(OUT) :: pdtc(ngrid,nlayer)
89	REAL,INTENT(OUT) :: pdtsrfc(ngrid)
90	REAL,INTENT(OUT) :: pdpsrf(ngrid)
91	REAL,INTENT(OUT) :: pduc(ngrid,nlayer)
92	REAL,INTENT(OUT) :: pdvc(ngrid,nlayer)
93	REAL,INTENT(OUT) :: pdqc(ngrid,nlayer,nq)
94
95	!-----------------------------------------------------------------------
96	! Local variables
97
98	INTEGER l,ig,icap,ilay,it,iq
99
100	REAL reffrad(ngrid,nlayer) ! radius (m) of the co2 ice particles
101	REAL*8 zt(ngrid,nlayermx)
102	REAL zq(ngrid,nlayermx,nq)
103	REAL zcpi
104	REAL ztcond (ngrid,nlayermx)
105	REAL ztnuc (ngrid,nlayermx)
106	REAL ztcondsol(ngrid)
107	REAL zdiceco2(ngrid)
108	REAL zcondicea(ngrid,nlayermx), zcondices(ngrid)
109	REAL zfallice(ngrid), Mfallice(ngrid)
110	REAL zmflux(nlayermx+1)
111	REAL zu(nlayermx),zv(nlayermx)
112	REAL ztsrf(ngrid)
113	REAL ztc(nlayermx), ztm(nlayermx+1)
114	REAL zum(nlayermx+1) , zvm(nlayermx+1)
115	LOGICAL condsub(ngrid)
116	REAL subptimestep
117	Integer Ntime
118	real masse (ngrid,nlayermx), w(ngrid,nlayermx,nq)
119	real wq(ngrid,nlayermx+1)
120	real vstokes,reff
121
122	! Special diagnostic variables
123	real tconda1(ngrid,nlayermx)
124	real tconda2(ngrid,nlayermx)
125	real zdtsig (ngrid,nlayermx)
126	real zdt (ngrid,nlayermx)
127
128	!-----------------------------------------------------------------------
129	! Saved local variables
130
131	REAL,SAVE :: latcond=5.9e5
132	REAL,SAVE :: ccond
133	REAL,SAVE :: cpice=1000.
134	REAL,SAVE,ALLOCATABLE,DIMENSION(:) :: emisref
135
136	LOGICAL,SAVE :: firstcall=.true.
137	REAL,EXTERNAL :: SSUM
138
139	REAL,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
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_cloud: 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,nlayermx) 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	zdiceco2(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,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,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	zdiceco2(ig) = zcondices(ig) + zfallice(ig)
438	piceco2(ig) = piceco2(ig) + zdiceco2(ig)*ptimestep
439	pdpsrf(ig) = -zdiceco2(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(ig.GT.ngrid/2+1) THEN
456	icap=2
457	ELSE
458	icap=1
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.0) then
467	psolaralb(ig) = albedice(icap)
468	emisref(ig) = emisice(icap)
469	else
470	psolaralb(ig) = albedodat(ig)
471	emisref(ig) = emissiv
472	pemisurf(ig) = emissiv
473	end if
474	end do
475
476	return
477	end subroutine condense_cloud
478
479	!-------------------------------------------------------------------------
480	subroutine get_tcond_co2(p,tcond)
481	! Calculates the condensation temperature for CO2
482
483	implicit none
484
485	#include "callkeys.h"
486
487	real p, peff, tcond
488	real, parameter :: ptriple=518000.0
489
490	peff=p
491
492	if(peff.lt.ptriple)then
493	tcond = (-3167.8)/(log(.01*peff)-23.23) ! Fanale's formula
494	else
495	tcond = 684.2-92.3log(peff)+4.32log(peff)**2
496	! liquid-vapour transition (based on CRC handbook 2003 data)
497	endif
498	return
499
500	end subroutine get_tcond_co2
501
502
503
504
505	!-------------------------------------------------------------------------
506	subroutine get_tnuc_co2(p,tnuc)
507	! Calculates the nucleation temperature for CO2, based on a simple super saturation criterion
508	! (JL 2011)
509
510	implicit none
511
512	#include "callkeys.h"
513
514	real p, peff, tnuc
515	real, parameter :: ptriple=518000.0
516
517	peff=p/co2supsat
518
519	if(peff.lt.ptriple)then
520	tnuc = (-3167.8)/(log(.01*peff)-23.23) ! Fanale's formula
521	else
522	tnuc = 684.2-92.3log(peff)+4.32log(peff)**2
523	! liquid-vapour transition (based on CRC handbook 2003 data)
524	endif
525	return
526
527	end subroutine get_tnuc_co2

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

Download in other formats:

Original Format