Context Navigation

condense_cloud.F90 @ 470

Last change on this file since 470 was 305, checked in by rwordsworth, 13 years ago
Several new files added as part of the climate evolution model (main program kcm.F90). Some general cleanup in physiq.F90 and callcorrk.F90. Bugs in dust radiative transfer and H2 Rayleigh scattering corrected.
Property svn:executable set to ``*
File size: 16.6 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, &
	6	pdqc,reffrad,cpp3D)
	7
	8	use radinc_h, only : naerkind
	9
	10	implicit none
	11
	12	!==================================================================
	13	! Purpose
	14	! -------
	15	! Condense and/or sublime CO2 ice on the ground and in the
	16	! atmosphere, and sediment the ice.
	17	!
	18	! Inputs
	19	! ------
	20	! ngrid Number of vertical columns
	21	! nlayer Number of layers
	22	! pplay(ngrid,nlayer) Pressure layers
	23	! pplev(ngrid,nlayer+1) Pressure levels
	24	! pt(ngrid,nlayer) Temperature (in K)
	25	! ptsrf(ngrid) Surface temperature
	26	!
	27	! pdt(ngrid,nlayermx) Time derivative before condensation/sublimation of pt
	28	! pdtsrf(ngrid) Time derivative before condensation/sublimation of ptsrf
	29	! pqsurf(ngrid,nq) Sedimentation flux at the surface (kg.m-2.s-1)
	30	!
	31	! Outputs
	32	! -------
	33	! pdpsrf(ngrid) \ Contribution of condensation/sublimation
	34	! pdtc(ngrid,nlayermx) / to the time derivatives of Ps, pt, and ptsrf
	35	! pdtsrfc(ngrid) /
	36	!
	37	! Both
	38	! ----
	39	! piceco2(ngrid) CO2 ice at the surface (kg/m2)
	40	! psolaralb(ngrid) Albedo at the surface
	41	! pemisurf(ngrid) Emissivity of the surface
	42	!
	43	! Authors
	44	! -------
	45	! Francois Forget (1996)
	46	! Converted to Fortran 90 and slightly modified by R. Wordsworth (2009)
	47	!
	48	!==================================================================
	49
	50	#include "dimensions.h"
	51	#include "dimphys.h"
	52	#include "comcstfi.h"
	53	#include "surfdat.h"
	54	#include "comgeomfi.h"
	55	#include "comvert.h"
	56	#include "callkeys.h"
	57	#include "tracer.h"
	58	#include "gases.h"
	59
	60	!-----------------------------------------------------------------------
	61	! Arguments
	62
	63	INTEGER ngrid, nlayer, nq
	64
	65	REAL ptimestep
	66	REAL pplay(ngrid,nlayer),pplev(ngrid,nlayer+1)
	67	REAL pcapcal(ngrid)
	68	REAL pt(ngrid,nlayer)
	69	REAL ptsrf(ngrid)
	70	REAL pphi(ngrid,nlayer)
	71	REAL pdt(ngrid,nlayer),pdtsrf(ngrid),pdtc(ngrid,nlayer)
	72	REAL pdtsrfc(ngrid),pdpsrf(ngrid)
	73	! REAL piceco2(ngrid),psolaralb(ngrid,2),pemisurf(ngrid)
	74	REAL piceco2(ngrid),psolaralb(ngrid),pemisurf(ngrid)
	75
	76	REAL pu(ngrid,nlayer) , pv(ngrid,nlayer)
	77	REAL pdu(ngrid,nlayer) , pdv(ngrid,nlayer)
	78	REAL pduc(ngrid,nlayer) , pdvc(ngrid,nlayer)
	79	REAL pq(ngridmx,nlayer,nq),pdq(ngrid,nlayer,nq)
	80	REAL pdqc(ngrid,nlayer,nq), pdqsc(ngrid,nq)
	81
	82	REAL reffrad(ngrid,nlayer,naerkind)
	83
	84	! Allow variations in cp with location
	85	REAL cpp3D(ngridmx,nlayermx) ! specific heat capacity at const. pressure
	86
	87
	88
	89	!-----------------------------------------------------------------------
	90	! Local variables
	91
	92	INTEGER l,ig,icap,ilay,it,iq
	93
	94	REAL*8 zt(ngridmx,nlayermx)
	95	REAL zq(ngridmx,nlayermx,nqmx)
	96	REAL zcpi
	97	REAL ztcond (ngridmx,nlayermx)
	98	REAL ztcondsol(ngridmx)
	99	REAL zdiceco2(ngridmx)
	100	REAL zcondicea(ngridmx,nlayermx), zcondices(ngridmx)
	101	REAL zfallice(ngridmx), Mfallice(ngridmx)
	102	REAL zmflux(nlayermx+1)
	103	REAL zu(nlayermx),zv(nlayermx)
	104	REAL ztsrf(ngridmx)
	105	REAL ztc(nlayermx), ztm(nlayermx+1)
	106	REAL zum(nlayermx+1) , zvm(nlayermx+1)
	107	LOGICAL condsub(ngridmx)
	108	REAL subptimestep
	109	Integer Ntime
	110	real masse (ngridmx,nlayermx), w(ngridmx,nlayermx,nqmx)
	111	real wq(ngridmx,nlayermx+1)
	112	real vstokes,reff
	113
	114	! Special diagnostic variables
	115	real tconda1(ngridmx,nlayermx)
	116	real tconda2(ngridmx,nlayermx)
	117	real zdtsig (ngridmx,nlayermx)
	118	real zdt (ngridmx,nlayermx)
	119
	120	!-----------------------------------------------------------------------
	121	! Saved local variables
	122
	123	REAL emisref(ngridmx)
	124	REAL latcond
	125	REAL ccond
	126	REAL cpice
	127	SAVE emisref, cpice
	128	SAVE latcond,ccond
	129
	130	LOGICAL firstcall
	131	SAVE firstcall
	132	REAL SSUM
	133	EXTERNAL SSUM
	134
	135	DATA latcond /5.9e5/
	136	DATA cpice /1000./
	137	DATA firstcall/.true./
	138
	139	REAL CBRT
	140	EXTERNAL CBRT
	141
	142	INTEGER,SAVE :: i_co2ice=0 ! co2 ice
	143	CHARACTER(LEN=20) :: tracername ! to temporarily store text
	144
	145	integer igas
	146	integer,save :: igasco2=0
	147	character(len=3) :: gasname
	148
	149	real reffradmin, reffradmax
	150
	151	real ppco2
	152
	153	!-----------------------------------------------------------------------
	154	! Initializations
	155
	156	call zerophys(ngridnlayernq,pdqc)
	157	call zerophys(ngridnlayernq,pdtc)
	158	call zerophys(ngridmxnlayermxnqmx,zq)
	159	call zerophys(ngridmx*nlayermx,zt)
	160
	161	!reffradmin=1.e-7
	162	!reffradmax=5.e-4
	163	!reffradmin=0.5e-7
	164	!reffradmax=0.1e-3 ! FF data
	165	reffradmin=0.1e-5
	166	reffradmax=0.1e-3 ! JB data
	167	! improve this in the future...
	168
	169	! Initialisation
	170	IF (firstcall) THEN
	171
	172	! find CO2 ice tracer
	173	do iq=1,nqmx
	174	tracername=noms(iq)
	175	if (tracername.eq."co2_ice") then
	176	i_co2ice=iq
	177	endif
	178	enddo
	179
	180	! find CO2 gas
	181	do igas=1,ngasmx
	182	gasname=gnom(igas)
	183	! gasname=noms(igas) ! was a bug
	184	if (gasname.eq."CO2") then
	185	igasco2=igas
	186	endif
	187	enddo
	188
	189	write(,) "condense_co2cloud: i_co2ice=",i_co2ice
	190
	191	if((i_co2ice.lt.1))then
	192	print*,'In condens_co2cloud but no CO2 ice tracer, exiting.'
	193	stop
	194	endif
	195
	196	ccond=cpp/(g*latcond)
	197	print*,'In condens_co2cloud: 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'
	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	ENDDO
	278	ENDDO
	279
	280	! Initialize zq and zt at the beginning of the sub-timestep loop
	281	DO l=1,nlayer
	282	DO ig=1,ngrid
	283	zt(ig,l)=pt(ig,l)
	284	zq(ig,l,i_co2ice)=pq(ig,l,i_co2ice)
	285	IF( zq(ig,l,i_co2ice).lt.-1.e-6 ) THEN
	286	print*,'Uh-oh, zq = ',zq(ig,l,i_co2ice),'at ig,l=',ig,l
	287	if(l.eq.1)then
	288	print*,'Perhaps the atmosphere is collapsing on surface...?'
	289	endif
	290	END IF
	291	ENDDO
	292	ENDDO
	293
	294	! Calculate the mass of each atmospheric layer (kg.m-2)
	295	do ilay=1,nlayer
	296	do ig=1, ngrid
	297	masse(ig,ilay)=(pplev(ig,ilay) - pplev(ig,ilay+1)) /g
	298	end do
	299	end do
	300
	301	! -----------------------------------------------
	302	! START CONDENSATION/SEDIMENTATION SUB-TIME LOOP
	303	! -----------------------------------------------
	304	Ntime = 20 ! number of sub-timestep
	305	subptimestep = ptimestep/float(Ntime)
	306
	307	DO it=1,Ntime
	308
	309	! Add the tendencies from other physical processes at each subtimstep
	310	DO l=1,nlayer
	311	DO ig=1,ngrid
	312	zt(ig,l) = zt(ig,l) + pdt(ig,l) * subptimestep
	313	zq(ig,l,i_co2ice) = zq(ig,l,i_co2ice) + pdq(ig,l,i_co2ice) * subptimestep
	314	END DO
	315	END DO
	316
	317
	318	! Gravitational sedimentation
	319
	320	! sedimentation computed from radius computed from q
	321	! assuming that the ice is splitted in Nmix particle
	322	do ilay=1,nlayer
	323	do ig=1, ngrid
	324
	325	reff = CBRT( 3zq(ig,ilay,i_co2ice)/( 4Nmix_co2pirho_co2 ))
	326
	327	! there should be a more elegant way of doing this...
	328	if(reff.lt.1.e-16) reff=1.e-16
	329
	330	! update reffrad for radiative transfer
	331	if(reff.lt.reffradmin)then
	332	reffrad(ig,ilay,1)=reffradmin
	333	!print*,'reff below optical limit'
	334	elseif(reff.gt.reffradmax)then
	335	reffrad(ig,ilay,1)=reffradmax
	336	!print*,'reff above optical limit'
	337	else
	338	reffrad(ig,ilay,1)=reff
	339	endif
	340
	341	call stokes &
	342	(pplev(ig,ilay),pt(ig,ilay), &
	343	reff,vstokes,rho_co2)
	344
	345	!w(ig,ilay,i_co2ice) = 0.0
	346	w(ig,ilay,i_co2ice) = vstokes * subptimestep * &
	347	pplev(ig,ilay)/(r*pt(ig,ilay))
	348
	349	end do
	350	end do
	351
	352	! Computing q after sedimentation
	353
	354	call vlz_fi(ngrid,zq(1,1,i_co2ice),2.,masse,w(1,1,i_co2ice),wq)
	355
	356
	357	! Progressively accumulating the flux to the ground
	358	! Mfallice is the total amount of ice fallen to the ground
	359	do ig=1,ngrid
	360	Mfallice(ig) = Mfallice(ig) + wq(ig,i_co2ice)
	361	end do
	362
	363
	364	! Condensation / sublimation in the atmosphere
	365	! --------------------------------------------
	366	! (calculation of zcondicea, zfallice and pdtc)
	367	! (MODIFICATIONS FOR EARLY MARS: falling heat neglected, condensation
	368	! of CO2 into tracer i_co2ice)
	369
	370	DO l=nlayer , 1, -1
	371	DO ig=1,ngrid
	372	pdtc(ig,l)=0.
	373
	374	! tweak ccond if the gas is non-ideal
	375	if(nonideal)then
	376	ccond=cpp3D(ig,l)/(g*latcond)
	377	endif
	378
	379
	380	IF ((zt(ig,l).LT.ztcond(ig,l)).or.(zq(ig,l,i_co2ice).gt.0)) THEN
	381	condsub(ig)=.true.
	382	pdtc(ig,l) = (ztcond(ig,l) - zt(ig,l))/subptimestep
	383	pdqc(ig,l,i_co2ice) = pdtc(ig,l)ccondg
	384
	385	! Case when the ice from above sublimes entirely
	386	IF ((zq(ig,l,i_co2ice).lt.-pdqc(ig,l,i_co2ice)*subptimestep) &
	387	.AND. (zq(ig,l,i_co2ice).gt.0)) THEN
	388
	389	pdqc(ig,l,i_co2ice) = -zq(ig,l,i_co2ice)/subptimestep
	390	pdtc(ig,l) =-zq(ig,l,i_co2ice)/(ccondgsubptimestep)
	391
	392	END IF
	393
	394	! Temperature and q after condensation
	395	zt(ig,l) = zt(ig,l) + pdtc(ig,l) * subptimestep
	396	zq(ig,l,i_co2ice) = zq(ig,l,i_co2ice) + pdqc(ig,l,i_co2ice) * subptimestep
	397	END IF
	398
	399	ENDDO
	400	ENDDO
	401	ENDDO ! end of subtimestep loop
	402
	403	! Computing global tendencies after the subtimestep
	404	DO l=1,nlayer
	405	DO ig=1,ngrid
	406	pdtc(ig,l) = &
	407	(zt(ig,l) - (pt(ig,l) + pdt(ig,l)*ptimestep))/ptimestep
	408	pdqc(ig,l,i_co2ice) = &
	409	(zq(ig,l,i_co2ice)-(pq(ig,l,i_co2ice)+pdq(ig,l,i_co2ice)*ptimestep))/ptimestep
	410	END DO
	411	END DO
	412	DO ig=1,ngrid
	413	zfallice(ig) = Mfallice(ig)/ptimestep
	414	END DO
	415
	416
	417	!-----------------------------------------------------------------------
	418	! Condensation/sublimation on the ground
	419	! (calculation of zcondices and pdtsrfc)
	420
	421	! forecast of ground temperature ztsrf and frost temperature ztcondsol
	422	DO ig=1,ngrid
	423	ppco2=gfrac(igasco2)*pplay(ig,1)
	424	call get_tcond_co2(ppco2,ztcondsol(ig))
	425
	426	ztsrf(ig) = ptsrf(ig)
	427
	428	if((ztsrf(ig).le.ztcondsol(ig)+2.0).and.(ngrid.eq.1))then
	429	print*,'CO2 is condensing on the surface in 1D. This atmosphere is doomed.'
	430	print*,'T_surf = ',ztsrf,'K'
	431	print*,'T_cond = ',ztcondsol,'K'
	432	open(116,file='surf_vals.out')
	433	write(116,*) 0.0, pplev(1,1), 0.0, 0.0
	434	close(116)
	435	call abort
	436	endif
	437
	438	ztsrf(ig) = ptsrf(ig) + pdtsrf(ig)*ptimestep
	439
	440	ENDDO
	441
	442	DO ig=1,ngrid
	443	IF(ig.GT.ngrid/2+1) THEN
	444	icap=2
	445	ELSE
	446	icap=1
	447	ENDIF
	448
	449	! Loop over where we have condensation / sublimation
	450	IF ((ztsrf(ig) .LT. ztcondsol(ig)) .OR. & ! ground condensation
	451	(zfallice(ig).NE.0.) .OR. & ! falling snow
	452	((ztsrf(ig) .GT. ztcondsol(ig)) .AND. & ! ground sublimation
	453	((piceco2(ig)+zfallice(ig)*ptimestep) .NE. 0.))) THEN
	454	condsub(ig) = .true.
	455
	456	! Condensation or partial sublimation of CO2 ice
	457	zcondices(ig)=pcapcal(ig)*(ztcondsol(ig)-ztsrf(ig)) &
	458	/(latcond*ptimestep)
	459	pdtsrfc(ig) = (ztcondsol(ig) - ztsrf(ig))/ptimestep
	460
	461	! If the entire CO_2 ice layer sublimes
	462	! (including what has just condensed in the atmosphere)
	463	IF((piceco2(ig)/ptimestep+zfallice(ig)).LE. &
	464	-zcondices(ig))THEN
	465	zcondices(ig) = -piceco2(ig)/ptimestep - zfallice(ig)
	466	pdtsrfc(ig)=(latcond/pcapcal(ig))* &
	467	(zcondices(ig))
	468	END IF
	469
	470	! Changing CO2 ice amount and pressure
	471
	472	zdiceco2(ig) = zcondices(ig) + zfallice(ig)
	473	piceco2(ig) = piceco2(ig) + zdiceco2(ig)*ptimestep
	474	pdpsrf(ig) = -zdiceco2(ig)*g
	475
	476	IF(ABS(pdpsrf(ig)*ptimestep).GT.pplev(ig,1)) THEN
	477	PRINT*,'STOP in condens'
	478	PRINT*,'condensing more than total mass'
	479	PRINT*,'Grid point ',ig
	480	PRINT*,'Ps = ',pplev(ig,1)
	481	PRINT*,'d Ps = ',pdpsrf(ig)
	482	STOP
	483	ENDIF
	484	END IF
	485	ENDDO
	486
	487	! Surface albedo and emissivity of the ground below the snow (emisref)
	488	! --------------------------------------------------------------------
	489	do ig =1,ngrid
	490	IF(ig.GT.ngrid/2+1) THEN
	491	icap=2
	492	ELSE
	493	icap=1
	494	ENDIF
	495
	496	if(.not.piceco2(ig).ge.0.) THEN
	497	if(piceco2(ig).le.-1.e-8) print*, &
	498	'WARNING in condense_co2cloud: piceco2(',ig,')=', piceco2(ig)
	499	piceco2(ig)=0.
	500	endif
	501	if (piceco2(ig).gt.0) then
	502	psolaralb(ig) = albedice(icap)
	503	emisref(ig) = emisice(icap)
	504	else
	505	psolaralb(ig) = albedodat(ig)
	506	emisref(ig) = emissiv
	507	pemisurf(ig) = emissiv
	508	end if
	509	end do
	510
	511	return
	512	end subroutine condense_cloud
	513
	514	!-------------------------------------------------------------------------
	515	subroutine get_tcond_co2(p,tcond)
	516	! Calculates the condensation temperature for CO2
	517
	518	implicit none
	519
	520	#include "callkeys.h"
	521
	522	real p, peff, tcond
	523	real, parameter :: ptriple=518000.0
	524
	525	peff=p!/co2supsat
	526
	527	if(peff.lt.ptriple)then
	528	tcond = (-3167.8)/(log(.01*peff)-23.23) ! Fanale's formula
	529	else
	530	tcond = 684.2-92.3log(peff)+4.32log(peff)**2
	531	! liquid-vapour transition (based on CRC handbook 2003 data)
	532	endif
	533	return
	534
	535	end subroutine get_tcond_co2

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

Download in other formats:

Original Format