Context Navigation

condense_cloud.F90 @ 537

Last change on this file since 537 was 526, checked in by jleconte, 13 years ago

13/02/2012 == JL + AS

All outputs are now in netCDF format. Even in 1D (No more G1D)
Clean up of the call to callcorrk when CLFvarying=true
Corrects a bug in writediagspecIR/VI. Output are now in W/m2/cm-1 as a function of the wavenumber in cm-1
Enable writediagspecIR/V to work in the CLFvarying=true case (output now done in Physiq after writediagfi)
Add a simple treatment for the supersaturation of CO2 (see forget et al 2012)
corrects a small bug when no clouds are present in aeropacity

Property svn:executable set to *

File size: 17.6 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,cpp3D)
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	! Allow variations in cp with location
86	REAL cpp3D(ngridmx,nlayermx) ! specific heat capacity at const. pressure
87
88
89
90	!-----------------------------------------------------------------------
91	! Local variables
92
93	INTEGER l,ig,icap,ilay,it,iq
94
95	REAL*8 zt(ngridmx,nlayermx)
96	REAL zq(ngridmx,nlayermx,nqmx)
97	REAL zcpi
98	REAL ztcond (ngridmx,nlayermx)
99	REAL ztnuc (ngridmx,nlayermx)
100	REAL ztcondsol(ngridmx)
101	REAL zdiceco2(ngridmx)
102	REAL zcondicea(ngridmx,nlayermx), zcondices(ngridmx)
103	REAL zfallice(ngridmx), Mfallice(ngridmx)
104	REAL zmflux(nlayermx+1)
105	REAL zu(nlayermx),zv(nlayermx)
106	REAL ztsrf(ngridmx)
107	REAL ztc(nlayermx), ztm(nlayermx+1)
108	REAL zum(nlayermx+1) , zvm(nlayermx+1)
109	LOGICAL condsub(ngridmx)
110	REAL subptimestep
111	Integer Ntime
112	real masse (ngridmx,nlayermx), w(ngridmx,nlayermx,nqmx)
113	real wq(ngridmx,nlayermx+1)
114	real vstokes,reff
115
116	! Special diagnostic variables
117	real tconda1(ngridmx,nlayermx)
118	real tconda2(ngridmx,nlayermx)
119	real zdtsig (ngridmx,nlayermx)
120	real zdt (ngridmx,nlayermx)
121
122	!-----------------------------------------------------------------------
123	! Saved local variables
124
125	REAL emisref(ngridmx)
126	REAL latcond
127	REAL ccond
128	REAL cpice
129	SAVE emisref, cpice
130	SAVE latcond,ccond
131
132	LOGICAL firstcall
133	SAVE firstcall
134	REAL SSUM
135	EXTERNAL SSUM
136
137	DATA latcond /5.9e5/
138	DATA cpice /1000./
139	DATA firstcall/.true./
140
141	REAL CBRT
142	EXTERNAL CBRT
143
144	INTEGER,SAVE :: i_co2ice=0 ! co2 ice
145	CHARACTER(LEN=20) :: tracername ! to temporarily store text
146
147	integer igas
148	integer,save :: igasco2=0
149	character(len=3) :: gasname
150
151	real reffradmin, reffradmax
152
153	real ppco2
154
155	!-----------------------------------------------------------------------
156	! Initializations
157
158	call zerophys(ngridnlayernq,pdqc)
159	call zerophys(ngridnlayernq,pdtc)
160	call zerophys(ngridmxnlayermxnqmx,zq)
161	call zerophys(ngridmx*nlayermx,zt)
162
163	!reffradmin=1.e-7
164	!reffradmax=5.e-4
165	!reffradmin=0.5e-7
166	!reffradmax=0.1e-3 ! FF data
167	reffradmin=0.1e-5
168	reffradmax=0.1e-3 ! JB data
169	! improve this in the future...
170
171	! Initialisation
172	IF (firstcall) THEN
173
174	! find CO2 ice tracer
175	do iq=1,nqmx
176	tracername=noms(iq)
177	if (tracername.eq."co2_ice") then
178	i_co2ice=iq
179	endif
180	enddo
181
182	! find CO2 gas
183	do igas=1,ngasmx
184	gasname=gnom(igas)
185	! gasname=noms(igas) ! was a bug
186	if (gasname.eq."CO2") then
187	igasco2=igas
188	endif
189	enddo
190
191	write(,) "condense_cloud: i_co2ice=",i_co2ice
192
193	if((i_co2ice.lt.1))then
194	print*,'In condens_cloud but no CO2 ice tracer, exiting.'
195	print*,'Still need generalisation to arbitrary species!'
196	stop
197	endif
198
199	ccond=cpp/(g*latcond)
200	print*,'In condens_cloud: ccond=',ccond,' latcond=',latcond
201
202	! Prepare special treatment if gas is not pure CO2
203	!if (addn2) then
204	! m_co2 = 44.01E-3 ! CO2 molecular mass (kg/mol)
205	! m_noco2 = 28.02E-3 ! N2 molecular mass (kg/mol)
206	! Compute A and B coefficient use to compute
207	! mean molecular mass Mair defined by
208	! 1/Mair = q(ico2)/m_co2 + (1-q(ico2))/m_noco2
209	! 1/Mair = A*q(ico2) + B
210	! A = (1/m_co2 - 1/m_noco2)
211	! B = 1/m_noco2
212	!endif
213
214	! Minimum CO2 mixing ratio below which mixing occurs with layer above:
215	!qco2min =0.75
216
217	firstcall=.false.
218	ENDIF
219	zcpi=1./cpp
220
221	!-----------------------------------------------------------------------
222	! Calculation of CO2 condensation / sublimation
223	!
224	! Variables used:
225	! piceco2(ngrid) amount of co2 ice on the ground (kg/m2)
226	! zcondicea(ngrid,l) condensation rate in layer l (kg/m2/s)
227	! zcondices(ngrid) condensation rate on the ground (kg/m2/s)
228	! zfallice(ngrid) flux of ice falling on surface (kg/m2/s)
229	! pdtc(ngrid,nlayermx) dT/dt due to phase changes (K/s)
230
231
232	! Tendencies initially set to 0 (except pdtc)
233	DO l=1,nlayer
234	DO ig=1,ngrid
235	zcondicea(ig,l) = 0.
236	pduc(ig,l) = 0
237	pdvc(ig,l) = 0
238	pdqc(ig,l,i_co2ice) = 0
239	END DO
240	END DO
241
242	DO ig=1,ngrid
243	Mfallice(ig) = 0.
244	zfallice(ig) = 0.
245	zcondices(ig) = 0.
246	pdtsrfc(ig) = 0.
247	pdpsrf(ig) = 0.
248	condsub(ig) = .false.
249	zdiceco2(ig) = 0.
250	ENDDO
251
252	!-----------------------------------------------------------------------
253	! Atmospheric condensation
254
255
256	! Compute CO2 Volume mixing ratio
257	! -------------------------------
258	! if (addn2) then
259	! DO l=1,nlayer
260	! DO ig=1,ngrid
261	! qco2=pq(ig,l,ico2)+pdq(ig,l,ico2)*ptimestep
262	!! Mean air molecular mass = 1/(q(ico2)/m_co2 + (1-q(ico2))/m_noco2)
263	! mmean=1/(A*qco2 +B)
264	! vmr_co2(ig,l) = qco2*mmean/m_co2
265	! ENDDO
266	! ENDDO
267	! else
268	! DO l=1,nlayer
269	! DO ig=1,ngrid
270	! vmr_co2(ig,l)=0.5
271	! ENDDO
272	! ENDDO
273	! end if
274
275	! Forecast the atmospheric frost temperature 'ztcond' and nucleation temperature 'ztnuc'
276	DO l=1,nlayer
277	DO ig=1,ngrid
278	ppco2=gfrac(igasco2)*pplay(ig,l)
279	call get_tcond_co2(ppco2,ztcond(ig,l))
280	call get_tnuc_co2(ppco2,ztnuc(ig,l))
281	ENDDO
282	ENDDO
283
284	! Initialize zq and zt at the beginning of the sub-timestep loop
285	DO l=1,nlayer
286	DO ig=1,ngrid
287	zt(ig,l)=pt(ig,l)
288	zq(ig,l,i_co2ice)=pq(ig,l,i_co2ice)
289	IF( zq(ig,l,i_co2ice).lt.-1.e-6 ) THEN
290	print*,'Uh-oh, zq = ',zq(ig,l,i_co2ice),'at ig,l=',ig,l
291	if(l.eq.1)then
292	print*,'Perhaps the atmosphere is collapsing on surface...?'
293	endif
294	END IF
295	ENDDO
296	ENDDO
297
298	! Calculate the mass of each atmospheric layer (kg.m-2)
299	do ilay=1,nlayer
300	do ig=1, ngrid
301	masse(ig,ilay)=(pplev(ig,ilay) - pplev(ig,ilay+1)) /g
302	end do
303	end do
304
305	! -----------------------------------------------
306	! START CONDENSATION/SEDIMENTATION SUB-TIME LOOP
307	! -----------------------------------------------
308	Ntime = 20 ! number of sub-timestep
309	subptimestep = ptimestep/float(Ntime)
310
311	DO it=1,Ntime
312
313	! Add the tendencies from other physical processes at each subtimstep
314	DO l=1,nlayer
315	DO ig=1,ngrid
316	zt(ig,l) = zt(ig,l) + pdt(ig,l) * subptimestep
317	zq(ig,l,i_co2ice) = zq(ig,l,i_co2ice) + pdq(ig,l,i_co2ice) * subptimestep
318	END DO
319	END DO
320
321
322	! Gravitational sedimentation
323
324	! sedimentation computed from radius computed from q
325	! assuming that the ice is splitted in Nmix particle
326	do ilay=1,nlayer
327	do ig=1, ngrid
328
329	reff = CBRT( 3zq(ig,ilay,i_co2ice)/( 4Nmix_co2pirho_co2 ))
330
331	! there should be a more elegant way of doing this...
332	if(reff.lt.1.e-16) reff=1.e-16
333
334	! update reffrad for radiative transfer
335	if(reff.lt.reffradmin)then
336	reffrad(ig,ilay,1)=reffradmin
337	!print*,'reff below optical limit'
338	elseif(reff.gt.reffradmax)then
339	reffrad(ig,ilay,1)=reffradmax
340	!print*,'reff above optical limit'
341	else
342	reffrad(ig,ilay,1)=reff
343	endif
344
345	call stokes &
346	(pplev(ig,ilay),pt(ig,ilay), &
347	reff,vstokes,rho_co2)
348
349	!w(ig,ilay,i_co2ice) = 0.0
350	w(ig,ilay,i_co2ice) = vstokes * subptimestep * &
351	pplev(ig,ilay)/(r*pt(ig,ilay))
352
353	end do
354	end do
355
356	! Computing q after sedimentation
357
358	call vlz_fi(ngrid,zq(1,1,i_co2ice),2.,masse,w(1,1,i_co2ice),wq)
359
360
361	! Progressively accumulating the flux to the ground
362	! Mfallice is the total amount of ice fallen to the ground
363	do ig=1,ngrid
364	Mfallice(ig) = Mfallice(ig) + wq(ig,i_co2ice)
365	end do
366
367
368	! Condensation / sublimation in the atmosphere
369	! --------------------------------------------
370	! (calculation of zcondicea, zfallice and pdtc)
371	! (MODIFICATIONS FOR EARLY MARS: falling heat neglected, condensation
372	! of CO2 into tracer i_co2ice)
373
374	DO l=nlayer , 1, -1
375	DO ig=1,ngrid
376	pdtc(ig,l)=0.
377
378	! tweak ccond if the gas is non-ideal
379	if(nonideal)then
380	ccond=cpp3D(ig,l)/(g*latcond)
381	endif
382
383
384	! ztcond-> ztnuc in test beneath to nucleate only when super saturation occurs(JL 2011)
385	IF ((zt(ig,l).LT.ztnuc(ig,l)).or.(zq(ig,l,i_co2ice).gt.1.E-10)) THEN
386	condsub(ig)=.true.
387	pdtc(ig,l) = (ztcond(ig,l) - zt(ig,l))/subptimestep
388	pdqc(ig,l,i_co2ice) = pdtc(ig,l)ccondg
389
390	! Case when the ice from above sublimes entirely
391	IF ((zq(ig,l,i_co2ice).lt.-pdqc(ig,l,i_co2ice)*subptimestep) &
392	.AND. (zq(ig,l,i_co2ice).gt.0)) THEN
393
394	pdqc(ig,l,i_co2ice) = -zq(ig,l,i_co2ice)/subptimestep
395	pdtc(ig,l) =-zq(ig,l,i_co2ice)/(ccondgsubptimestep)
396
397	END IF
398
399	! Temperature and q after condensation
400	zt(ig,l) = zt(ig,l) + pdtc(ig,l) * subptimestep
401	zq(ig,l,i_co2ice) = zq(ig,l,i_co2ice) + pdqc(ig,l,i_co2ice) * subptimestep
402	END IF
403
404	ENDDO
405	ENDDO
406	ENDDO ! end of subtimestep loop
407
408	! Computing global tendencies after the subtimestep
409	DO l=1,nlayer
410	DO ig=1,ngrid
411	pdtc(ig,l) = &
412	(zt(ig,l) - (pt(ig,l) + pdt(ig,l)*ptimestep))/ptimestep
413	pdqc(ig,l,i_co2ice) = &
414	(zq(ig,l,i_co2ice)-(pq(ig,l,i_co2ice)+pdq(ig,l,i_co2ice)*ptimestep))/ptimestep
415	END DO
416	END DO
417	DO ig=1,ngrid
418	zfallice(ig) = Mfallice(ig)/ptimestep
419	END DO
420
421
422	!-----------------------------------------------------------------------
423	! Condensation/sublimation on the ground
424	! (calculation of zcondices and pdtsrfc)
425
426	! forecast of ground temperature ztsrf and frost temperature ztcondsol
427	DO ig=1,ngrid
428	ppco2=gfrac(igasco2)*pplay(ig,1)
429	call get_tcond_co2(ppco2,ztcondsol(ig))
430
431	ztsrf(ig) = ptsrf(ig)
432
433	if((ztsrf(ig).le.ztcondsol(ig)+2.0).and.(ngrid.eq.1))then
434	print*,'CO2 is condensing on the surface in 1D. This atmosphere is doomed.'
435	print*,'T_surf = ',ztsrf,'K'
436	print*,'T_cond = ',ztcondsol,'K'
437	open(116,file='surf_vals.out')
438	write(116,*) 0.0, pplev(1,1), 0.0, 0.0
439	close(116)
440	call abort
441	endif
442
443	ztsrf(ig) = ptsrf(ig) + pdtsrf(ig)*ptimestep
444
445	ENDDO
446
447	DO ig=1,ngrid
448	IF(ig.GT.ngrid/2+1) THEN
449	icap=2
450	ELSE
451	icap=1
452	ENDIF
453
454	! Loop over where we have condensation / sublimation
455	IF ((ztsrf(ig) .LT. ztcondsol(ig)) .OR. & ! ground condensation
456	(zfallice(ig).NE.0.) .OR. & ! falling snow
457	((ztsrf(ig) .GT. ztcondsol(ig)) .AND. & ! ground sublimation
458	((piceco2(ig)+zfallice(ig)*ptimestep) .NE. 0.))) THEN
459	condsub(ig) = .true.
460
461	! Condensation or partial sublimation of CO2 ice
462	zcondices(ig)=pcapcal(ig)*(ztcondsol(ig)-ztsrf(ig)) &
463	/(latcond*ptimestep)
464	pdtsrfc(ig) = (ztcondsol(ig) - ztsrf(ig))/ptimestep
465
466	! If the entire CO_2 ice layer sublimes
467	! (including what has just condensed in the atmosphere)
468	IF((piceco2(ig)/ptimestep+zfallice(ig)).LE. &
469	-zcondices(ig))THEN
470	zcondices(ig) = -piceco2(ig)/ptimestep - zfallice(ig)
471	pdtsrfc(ig)=(latcond/pcapcal(ig))* &
472	(zcondices(ig))
473	END IF
474
475	! Changing CO2 ice amount and pressure
476
477	zdiceco2(ig) = zcondices(ig) + zfallice(ig)
478	piceco2(ig) = piceco2(ig) + zdiceco2(ig)*ptimestep
479	pdpsrf(ig) = -zdiceco2(ig)*g
480
481	IF(ABS(pdpsrf(ig)*ptimestep).GT.pplev(ig,1)) THEN
482	PRINT*,'STOP in condens'
483	PRINT*,'condensing more than total mass'
484	PRINT*,'Grid point ',ig
485	PRINT*,'Ps = ',pplev(ig,1)
486	PRINT*,'d Ps = ',pdpsrf(ig)
487	STOP
488	ENDIF
489	END IF
490	ENDDO
491
492	! Surface albedo and emissivity of the ground below the snow (emisref)
493	! --------------------------------------------------------------------
494	do ig =1,ngrid
495	IF(ig.GT.ngrid/2+1) THEN
496	icap=2
497	ELSE
498	icap=1
499	ENDIF
500
501	if(.not.piceco2(ig).ge.0.) THEN
502	if(piceco2(ig).le.-1.e-8) print*, &
503	'WARNING in condense_co2cloud: piceco2(',ig,')=', piceco2(ig)
504	piceco2(ig)=0.
505	endif
506	if (piceco2(ig).gt.0) then
507	psolaralb(ig) = albedice(icap)
508	emisref(ig) = emisice(icap)
509	else
510	psolaralb(ig) = albedodat(ig)
511	emisref(ig) = emissiv
512	pemisurf(ig) = emissiv
513	end if
514	end do
515
516	return
517	end subroutine condense_cloud
518
519	!-------------------------------------------------------------------------
520	subroutine get_tcond_co2(p,tcond)
521	! Calculates the condensation temperature for CO2
522
523	implicit none
524
525	#include "callkeys.h"
526
527	real p, peff, tcond
528	real, parameter :: ptriple=518000.0
529
530	peff=p
531
532	if(peff.lt.ptriple)then
533	tcond = (-3167.8)/(log(.01*peff)-23.23) ! Fanale's formula
534	else
535	tcond = 684.2-92.3log(peff)+4.32log(peff)**2
536	! liquid-vapour transition (based on CRC handbook 2003 data)
537	endif
538	return
539
540	end subroutine get_tcond_co2
541
542
543
544
545	!-------------------------------------------------------------------------
546	subroutine get_tnuc_co2(p,tnuc)
547	! Calculates the nucleation temperature for CO2, based on a simple super saturation criterion
548	! (JL 2011)
549
550	implicit none
551
552	#include "callkeys.h"
553
554	real p, peff, tnuc
555	real, parameter :: ptriple=518000.0
556
557	peff=p/co2supsat
558
559	if(peff.lt.ptriple)then
560	tnuc = (-3167.8)/(log(.01*peff)-23.23) ! Fanale's formula
561	else
562	tnuc = 684.2-92.3log(peff)+4.32log(peff)**2
563	! liquid-vapour transition (based on CRC handbook 2003 data)
564	endif
565	return
566
567	end subroutine get_tnuc_co2

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

Download in other formats:

Original Format