Context Navigation

condense_cloud.F90 @ 1330

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

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

Download in other formats:

Original Format