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

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
[1216]	9	use gases_h, only: gfrac, igas_co2
[728]	10	use radii_mod, only : co2_reffrad
	11	use aerosol_mod, only : iaero_co2
[1216]	12	USE surfdat_h, only: albedodat, albedice, emisice, emissiv
	13	USE comgeomfi_h, only: lati
	14	USE tracer_h, only: noms, rho_co2
[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	!
[1308]	34	! pdt(ngrid,nlayer) Time derivative before condensation/sublimation of pt
[305]	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	! -------
[1308]	40	! pdpsrf(ngrid) \ Contribution of condensation/sublimation
	41	! pdtc(ngrid,nlayer) / to the time derivatives of Ps, pt, and ptsrf
	42	! pdtsrfc(ngrid) /
[305]	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
[1308]	58	!#include "dimensions.h"
	59	!#include "dimphys.h"
[305]	60	#include "comcstfi.h"
[1308]	61	!#include "comvert.h"
[305]	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
[1308]	101	REAL*8 zt(ngrid,nlayer)
	102	REAL zq(ngrid,nlayer,nq)
[305]	103	REAL zcpi
[1308]	104	REAL ztcond (ngrid,nlayer)
	105	REAL ztnuc (ngrid,nlayer)
[787]	106	REAL ztcondsol(ngrid)
	107	REAL zdiceco2(ngrid)
[1308]	108	REAL zcondicea(ngrid,nlayer), zcondices(ngrid)
[787]	109	REAL zfallice(ngrid), Mfallice(ngrid)
[1308]	110	REAL zmflux(nlayer+1)
	111	REAL zu(nlayer),zv(nlayer)
[787]	112	REAL ztsrf(ngrid)
[1308]	113	REAL ztc(nlayer), ztm(nlayer+1)
	114	REAL zum(nlayer+1) , zvm(nlayer+1)
[787]	115	LOGICAL condsub(ngrid)
[305]	116	REAL subptimestep
	117	Integer Ntime
[1308]	118	real masse (ngrid,nlayer), w(ngrid,nlayer,nq)
	119	real wq(ngrid,nlayer+1)
[305]	120	real vstokes,reff
	121
	122	! Special diagnostic variables
[1308]	123	real tconda1(ngrid,nlayer)
	124	real tconda2(ngrid,nlayer)
	125	real zdtsig (ngrid,nlayer)
	126	real zdt (ngrid,nlayer)
[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
[1315]	135	!$OMP THREADPRIVATE(latcond,ccond,cpice,emisref)
[305]	136
[858]	137	LOGICAL,SAVE :: firstcall=.true.
[1315]	138	!$OMP THREADPRIVATE(firstcall)
[858]	139	REAL,EXTERNAL :: SSUM
[305]	140
[858]	141	REAL,EXTERNAL :: CBRT
[305]	142
	143	INTEGER,SAVE :: i_co2ice=0 ! co2 ice
[1315]	144	!$OMP THREADPRIVATE(i_co2ice)
[305]	145	CHARACTER(LEN=20) :: tracername ! to temporarily store text
	146
	147	integer igas
	148
	149	real ppco2
	150
	151	!-----------------------------------------------------------------------
	152	! Initializations
	153
[858]	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
[305]	158
	159	! Initialisation
	160	IF (firstcall) THEN
	161
[787]	162	ALLOCATE(emisref(ngrid)) !! this should be deallocated in lastcall...
	163
[305]	164	! find CO2 ice tracer
[787]	165	do iq=1,nq
[305]	166	tracername=noms(iq)
	167	if (tracername.eq."co2_ice") then
	168	i_co2ice=iq
	169	endif
	170	enddo
	171
[486]	172	write(,) "condense_cloud: i_co2ice=",i_co2ice
[305]	173
	174	if((i_co2ice.lt.1))then
[486]	175	print*,'In condens_cloud but no CO2 ice tracer, exiting.'
	176	print*,'Still need generalisation to arbitrary species!'
[305]	177	stop
	178	endif
	179
	180	ccond=cpp/(g*latcond)
[486]	181	print*,'In condens_cloud: ccond=',ccond,' latcond=',latcond
[305]	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)
[1308]	210	! pdtc(ngrid,nlayer) dT/dt due to phase changes (K/s)
[305]	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
[526]	256	! Forecast the atmospheric frost temperature 'ztcond' and nucleation temperature 'ztnuc'
[305]	257	DO l=1,nlayer
	258	DO ig=1,ngrid
[869]	259	ppco2=gfrac(igas_CO2)*pplay(ig,l)
[305]	260	call get_tcond_co2(ppco2,ztcond(ig,l))
[526]	261	call get_tnuc_co2(ppco2,ztnuc(ig,l))
[305]	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
[787]	281	DO ig=1,ngrid
[305]	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
[728]	305	! sedimentation computed from radius computed from q in module radii_mod
[1308]	306	call co2_reffrad(ngrid,nlayer,nq,zq,reffrad)
[728]	307
[305]	308	do ilay=1,nlayer
[787]	309	DO ig=1,ngrid
[305]	310
[858]	311	reff = reffrad(ig,ilay)
[305]	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
[1308]	326	call vlz_fi(ngrid,nlayer,zq(1,1,i_co2ice),2.,masse,w(1,1,i_co2ice),wq)
[305]	327
	328
	329	! Progressively accumulating the flux to the ground
	330	! Mfallice is the total amount of ice fallen to the ground
[787]	331	DO ig=1,ngrid
[305]	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
[526]	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
[305]	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
[869]	391	ppco2=gfrac(igas_CO2)*pplay(ig,1)
[305]	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	! --------------------------------------------------------------------
[787]	457	DO ig=1,ngrid
[1216]	458	IF(lati(ig).LT.0.) THEN
	459	icap=2 ! Southern Hemisphere
[305]	460	ELSE
[1216]	461	icap=1 ! Nortnern hemisphere
[305]	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
[526]	493	peff=p
[305]	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
[526]	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