Context Navigation

condense_cloud.F90 @ 868

Last change on this file since 868 was 858, checked in by emillour, 13 years ago

Generic GCM:

Fixed sedimentation issue: ensure in callsedim that the correct radii are provided to newsedim and also that the updated temperature and tracer mixing ratios are used to compute sedimentation.
Updated the way aerosol radii are considered and used; routines in radii_mod (h2o_reffrad, co2_reffrad, etc.) only handle a single aerosol. The idea here is that these can be called from anywhere and that the caller doesn't need to have the full (naerkind size) array of aerosol radii.
cleanup (addition of intent(..) to routine arguments) in various routines

Property svn:executable set to *

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

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

Download in other formats:

Original Format