Context Navigation

condense_cloud.F90 @ 724

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

Rev	Line
[305]	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, &
[586]	6	pdqc,reffrad)
[305]	7
	8	use radinc_h, only : naerkind
[471]	9	use gases_h
[305]	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)
[526]	48	! Includes simplifed nucleation by J. Leconte (2011)
[305]	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)
[526]	96	REAL ztnuc (ngridmx,nlayermx)
[305]	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
[486]	188	write(,) "condense_cloud: i_co2ice=",i_co2ice
[305]	189
	190	if((i_co2ice.lt.1))then
[486]	191	print*,'In condens_cloud but no CO2 ice tracer, exiting.'
	192	print*,'Still need generalisation to arbitrary species!'
[305]	193	stop
	194	endif
	195
	196	ccond=cpp/(g*latcond)
[486]	197	print*,'In condens_cloud: ccond=',ccond,' latcond=',latcond
[305]	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
[526]	272	! Forecast the atmospheric frost temperature 'ztcond' and nucleation temperature 'ztnuc'
[305]	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))
[526]	277	call get_tnuc_co2(ppco2,ztnuc(ig,l))
[305]	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
[526]	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
[305]	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
[526]	522	peff=p
[305]	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
[526]	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