Context Navigation

newcondens.F @ 1009

Last change on this file since 1009 was 890, checked in by emillour, 12 years ago
Mars GCM: minor bug correction in newcondens.F and some code tidying (in co2snow.F also). EM
File size: 31.9 KB

Rev	Line
[38]	1	SUBROUTINE newcondens(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,pdqc,
	6	$ fluxsurf_sw,zls)
	7
	8	IMPLICIT NONE
	9	c=======================================================================
	10	c subject:
	11	c --------
	12	c Condensation/sublimation of CO2 ice on the ground and in the
	13	c atmosphere
	14	c (Scheme described in Forget et al., Icarus, 1998)
	15	c
	16	c author: Francois Forget 1994-1996
	17	c ------
	18	c
	19	c input:
	20	c -----
	21	c ngrid nombre de points de verticales
	22	c (toutes les boucles de la physique sont au
	23	c moins vectorisees sur ngrid)
	24	c nlayer nombre de couches
	25	c pplay(ngrid,nlayer) Pressure levels
	26	c pplev(ngrid,nlayer+1) Pressure levels
	27	c pt(ngrid,nlayer) temperature (en K)
	28	c ptsrf(ngrid) temperature de surface
	29	c
	30	c \
	31	c pdt(ngrid,nlayermx) \ derivee temporelle physique avant condensation
	32	c / ou sublimation pour pt,ptsrf
	33	c pdtsrf(ngrid) /
	34	c
	35	c output:
	36	c -------
	37	c
	38	c pdpsrf(ngrid) \ derivee temporelle physique (contribution de
	39	c pdtc(ngrid,nlayermx) / la condensation ou sublimation) pour Ps,pt,ptsrf
	40	c pdtsrfc(ngrid) /
	41	c
	42	c Entree/sortie
	43	c -------------
	44	c
	45	c piceco2(ngrid) : quantite de glace co2 au sol (kg/m2)
	46	c psolaralb(ngrid,2) : albedo au sol
	47	c pemisurf(ngrid) : emissivite du sol
	48
	49	c
	50	c=======================================================================
	51	c
	52	c 0. Declarations :
	53	c ------------------
	54	c
	55	#include "dimensions.h"
	56	#include "dimphys.h"
	57	#include "comcstfi.h"
	58	#include "surfdat.h"
	59	#include "comgeomfi.h"
	60	#include "comvert.h"
	61	#include "paramet.h"
	62	#include "callkeys.h"
	63	#include "tracer.h"
	64
	65	c-----------------------------------------------------------------------
	66	c Arguments :
	67	c ---------
[890]	68	INTEGER,INTENT(IN) :: ngrid ! number of atmospheric columns
	69	INTEGER,INTENT(IN) :: nlayer ! number of vertical layers
	70	INTEGER,INTENT(IN) :: nq ! number of tracers
[38]	71
[890]	72	REAL,INTENT(IN) :: ptimestep ! physics timestep (s)
	73	REAL,INTENT(IN) :: pcapcal(ngrid)
	74	REAL,INTENT(IN) :: pplay(ngrid,nlayer) !mid-layer pressure (Pa)
	75	REAL,INTENT(IN) :: pplev(ngrid,nlayer+1) ! inter-layer pressure (Pa)
	76	REAL,INTENT(IN) :: ptsrf(ngrid) ! surface temperature (K)
	77	REAL,INTENT(IN) :: pt(ngrid,nlayer) ! atmospheric temperature (K)
	78	REAL,INTENT(IN) :: pphi(ngrid,nlayer) ! geopotential (m2.s-2)
	79	REAL,INTENT(IN) :: pdt(ngrid,nlayer) ! tendency on temperature from
	80	! previous physical processes (K/s)
	81	REAL,INTENT(IN) :: pdu(ngrid,nlayer) ! tendency on zonal wind (m/s2)
	82	! from previous physical processes
	83	REAL,INTENT(IN) :: pdv(ngrid,nlayer) ! tendency on meridional wind (m/s2)
	84	! from previous physical processes
	85	REAL,INTENT(IN) :: pdtsrf(ngrid) ! tendency on surface temperature from
	86	! previous physical processes (K/s)
	87	REAL,INTENT(IN) :: pu(ngrid,nlayer) ! zonal wind (m/s)
	88	REAL,INTENT(IN) :: pv(ngrid,nlayer) ! meridional wind (m/s)
	89	REAL,INTENT(IN) :: pq(ngrid,nlayer,nq) ! tracers (../kg_air)
	90	REAL,INTENT(IN) :: pdq(ngrid,nlayer,nq) ! tendency on tracers from
	91	! previous physical processes
	92	REAL,INTENT(INOUT) :: piceco2(ngrid) ! CO2 ice on the surface (kg.m-2)
	93	REAL,INTENT(INOUT) :: psolaralb(ngrid,2) ! albedo of the surface
	94	REAL,INTENT(INOUT) :: pemisurf(ngrid) ! emissivity of the surface
	95
	96	! tendencies due to CO2 condensation/sublimation:
	97	REAL,INTENT(OUT) :: pdtc(ngrid,nlayer) ! tendency on temperature (K/s)
	98	REAL,INTENT(OUT) :: pdtsrfc(ngrid) ! tendency on surface temperature (K/s)
	99	REAL,INTENT(OUT) :: pdpsrf(ngrid) ! tendency on surface pressure (Pa/s)
	100	REAL,INTENT(OUT) :: pduc(ngrid,nlayer) ! tendency on zonal wind (m.s-2)
	101	REAL,INTENT(OUT) :: pdvc(ngrid,nlayer) ! tendency on meridional wind (m.s-2)
	102	REAL,INTENT(OUT) :: pdqc(ngrid,nlayer,nq) ! tendency on tracers
	103
	104	! added to calculate flux dependent albedo:
	105	REAL,intent(in) :: fluxsurf_sw(ngrid,2)
	106	real,intent(in) :: zls ! solar longitude (rad)
[38]	107
	108	c
	109	c Local variables :
	110	c -----------------
	111
	112	c variables used for albedo parametrization
	113	c --------------------------------------------
	114	INTEGER i,j
	115	REAL Fluxmean(jjp1)
	116	INTEGER l,ig,iq,icap,nmix
	117	LOGICAL transparency, fluxdependent
	118	c flag transparency if you want to make the co2ice semi-transparent
	119	PARAMETER(transparency=.true.)
	120	c flag fluxdependent if you want the co2ice albedo to be dependent on
	121	c the incident solar flux
	122	PARAMETER(fluxdependent=.false.)
[890]	123	REAL slopy,alpha,constA,constB,constT,albediceF_new(ngrid)
	124	REAL zt(ngrid,nlayer)
[38]	125	REAL zcpi
[890]	126	REAL ztcond (ngrid,nlayer+1)
	127	REAL ztcondsol(ngrid)
	128	REAL zdiceco2(ngrid)
	129	REAL zcondicea(ngrid,nlayer)
	130	REAL zcondices(ngrid)
	131	REAL zfallice(ngrid,nlayer+1) , zfallheat
	132	REAL zmflux(nlayer+1)
	133	REAL zu(nlayer),zv(nlayer)
	134	REAL zq(nlayer,nq),zq1(nlayer)
	135	REAL ztsrf(ngrid)
	136	REAL ztc(nlayer), ztm(nlayer+1)
	137	REAL zum(nlayer+1) , zvm(nlayer+1)
	138	REAL zqm(nlayer+1,nq),zqm1(nlayer+1)
	139	REAL masse(nlayer),w(nlayer+1)
	140	REAL Sm(nlayer),Smq(nlayer,nq),mixmas,qmix
	141	LOGICAL condsub(ngrid)
[38]	142
	143	c variable speciale diagnostique
[890]	144	real tconda1(ngrid,nlayer)
	145	real tconda2(ngrid,nlayer)
	146	c REAL zdiceco2a(ngrid) ! for diagnostic only
	147	real zdtsig (ngrid,nlayer)
	148	real zdt (ngrid,nlayer)
	149	real vmr_co2(ngrid,nlayer) ! co2 volume mixing ratio
[38]	150	! improved_ztcond flag: If set to .true. (AND running with a 'co2' tracer)
	151	! then condensation temperature is computed using partial pressure of CO2
	152	logical,parameter :: improved_ztcond=.true.
	153	! Bound co2 (tracer) values...
	154	logical,parameter :: bound_qco2=.false.
	155	real,parameter :: qco2max=1.1
	156	real,parameter :: qco2mini=0.1
	157	real :: zqco2
	158
	159	c local saved variables
[890]	160	integer,save :: ico2 ! index of CO2 tracer
	161	real,save :: qco2min,qco2,mmean
	162	real,allocatable,save :: emisref(:)
	163	real,parameter :: latcond=5.9e5 ! (J/kg) Latent heat of solid CO2 ice
	164	real,parameter :: tcond1mb=136.27 ! condensation temperature (K) at 1 mbar
	165	real,parameter :: cpice=1000. ! (J.kg-1.K-1) specific heat of CO2 ice
	166	REAL,SAVE :: acond,bcond,ccond
	167	! REAL,SAVE :: albediceF(ngridmx)
	168	real,save :: m_co2, m_noco2, A , B
[38]	169
[890]	170	LOGICAL,SAVE :: firstcall = .true. !,firstcall2=.true.
[38]	171
	172	integer flag
	173
	174	c----------------------------------------------------------------------
	175
	176	c Initialisation
	177	c --------------
	178	c
	179	IF (firstcall) THEN
[890]	180
	181	allocate(emisref(ngrid))
	182
[38]	183	bcond=1./tcond1mb
	184	ccond=cpp/(g*latcond)
	185	acond=r/latcond
	186
	187	firstcall=.false.
	188	write(,) 'Newcondens: improved_ztcond=',improved_ztcond
	189	write(,) 'Newcondens: bound_qco2=',bound_qco2
	190	PRINT*,'In newcondens:Tcond(P=1mb)=',tcond1mb,' Lcond=',latcond
	191	PRINT*,'acond,bcond,ccond',acond,bcond,ccond
	192
	193	ico2=0
	194
	195	if (tracer) then
	196	c Prepare Special treatment if one of the tracer is CO2 gas
	197	do iq=1,nqmx
	198	if (noms(iq).eq."co2") then
	199	ico2=iq
	200	m_co2 = 44.01E-3 ! CO2 molecular mass (kg/mol)
	201	m_noco2 = 33.37E-3 ! Non condensible mol mass (kg/mol)
	202	c Compute A and B coefficient use to compute
	203	c mean molecular mass Mair defined by
	204	c 1/Mair = q(ico2)/m_co2 + (1-q(ico2))/m_noco2
	205	c 1/Mair = A*q(ico2) + B
	206	A =(1/m_co2 - 1/m_noco2)
	207	B=1/m_noco2
	208	endif
	209	enddo
	210	c minimum CO2 mix. ratio below which mixing occur with layer above:
	211	qco2min =0.75
	212	end if
[890]	213	ENDIF ! of IF (firstcall)
[38]	214	zcpi=1./cpp
	215	c
	216	c======================================================================
	217	c Calcul of CO2 condensation sublimation
	218	c ============================================================
	219	c
	220	c Used variable :
	221	c piceco2(ngrid) : amount of co2 ice on the ground (kg/m2)
	222	c zcondicea(ngrid,l): condensation rate in layer l (kg/m2/s)
	223	c zcondices(ngrid): condensation rate on the ground (kg/m2/s)
	224	c zfallice(ngrid,l):amount of ice falling from layer l (kg/m2/s)
	225	c
	226	c pdtc(ngrid,nlayermx) : dT/dt due to cond/sub
	227	c
	228	c
	229	c Tendencies set to 0 (except pdtc)
	230	c -------------------------------------
	231	DO l=1,nlayer
	232	DO ig=1,ngrid
	233	zcondicea(ig,l) = 0.
	234	zfallice(ig,l) = 0.
	235	pduc(ig,l) = 0
	236	pdvc(ig,l) = 0
	237	END DO
	238	END DO
	239
	240	DO iq=1,nqmx
	241	DO l=1,nlayer
	242	DO ig=1,ngrid
	243	pdqc(ig,l,iq) = 0
	244	END DO
	245	END DO
	246	END DO
	247
	248	DO ig=1,ngrid
	249	zfallice(ig,nlayer+1) = 0.
	250	zcondices(ig) = 0.
	251	pdtsrfc(ig) = 0.
	252	pdpsrf(ig) = 0.
	253	condsub(ig) = .false.
	254	zdiceco2(ig) = 0.
	255	ENDDO
	256	zfallheat=0
	257
	258	c *************************
	259	c ATMOSPHERIC CONDENSATION
	260	c *************************
	261
	262	c Compute CO2 Volume mixing ratio
	263	c -------------------------------
	264	if (improved_ztcond.and.(ico2.ne.0)) then
	265	DO l=1,nlayer
	266	DO ig=1,ngrid
	267	qco2=pq(ig,l,ico2)+pdq(ig,l,ico2)*ptimestep
	268	c Mean air molecular mass = 1/(q(ico2)/m_co2 + (1-q(ico2))/m_noco2)
	269	mmean=1/(A*qco2 +B)
	270	vmr_co2(ig,l) = qco2*mmean/m_co2
	271	ENDDO
	272	ENDDO
	273	else
	274	DO l=1,nlayer
	275	DO ig=1,ngrid
	276	vmr_co2(ig,l)=0.95
	277	ENDDO
	278	ENDDO
	279	end if
	280
	281	c forecast of atmospheric temperature zt and frost temperature ztcond
	282	c --------------------------------------------------------------------
	283
	284	DO l=1,nlayer
	285	DO ig=1,ngrid
	286	zt(ig,l)=pt(ig,l)+ pdt(ig,l)*ptimestep
	287	! ztcond(ig,l)=1./(bcond-acondlog(.0095pplay(ig,l)))
	288	ztcond(ig,l)=
	289	& 1./(bcond-acondlog(.01vmr_co2(ig,l)*pplay(ig,l)))
	290	if (pplay(ig,l).lt.1e-4) ztcond(ig,l)=0.0 !mars Monica
	291	ENDDO
	292	ENDDO
[327]	293
[328]	294	ztcond(:,nlayer+1)=ztcond(:,nlayer)
[38]	295
	296	c Condensation/sublimation in the atmosphere
	297	c ------------------------------------------
	298	c (calcul of zcondicea , zfallice and pdtc)
	299	c
	300	DO l=nlayer , 1, -1
	301	DO ig=1,ngrid
	302	pdtc(ig,l)=0.
	303	IF((zt(ig,l).LT.ztcond(ig,l)).or.(zfallice(ig,l+1).gt.0))THEN
	304	condsub(ig)=.true.
	305	IF (zfallice(ig,l+1).gt.0) then
	306	zfallheat=zfallice(ig,l+1)*
	307	& (pphi(ig,l+1)-pphi(ig,l) +
	308	& cpice*(ztcond(ig,l+1)-ztcond(ig,l)))/latcond
	309	ELSE
	310	zfallheat=0.
	311	ENDIF
	312	pdtc(ig,l)=(ztcond(ig,l) - zt(ig,l))/ptimestep
	313	zcondicea(ig,l)=(pplev(ig,l)-pplev(ig,l+1))
	314	& ccondpdtc(ig,l)- zfallheat
	315	c Case when the ice from above sublimes entirely
	316	c """""""""""""""""""""""""""""""""""""""""""""""
	317	IF (zfallice(ig,l+1).lt.- zcondicea(ig,l)) then
	318	pdtc(ig,l)=(-zfallice(ig,l+1)+zfallheat)/
	319	& (ccond*(pplev(ig,l)-pplev(ig,l+1)))
	320	zcondicea(ig,l)= -zfallice(ig,l+1)
	321	END IF
	322
	323	zfallice(ig,l) = zcondicea(ig,l)+zfallice(ig,l+1)
	324	END IF
	325	ENDDO
	326	ENDDO
	327
	328	c *************************
	329	c SURFACE CONDENSATION
	330	c *************************
	331
	332	c forecast of ground temperature ztsrf and frost temperature ztcondsol
	333	c --------------------------------------------------------------------
	334	DO ig=1,ngrid
	335	ztcondsol(ig)=
	336	& 1./(bcond-acondlog(.01vmr_co2(ig,1)*pplev(ig,1)))
	337	ztsrf(ig) = ptsrf(ig) + pdtsrf(ig)*ptimestep
	338	ENDDO
	339
	340	c
	341	c Condensation/sublimation on the ground
	342	c --------------------------------------
	343	c (calcul of zcondices , pdtsrfc)
	344	c
	345	DO ig=1,ngrid
	346	IF(ig.GT.ngrid/2+1) THEN
	347	icap=2
	348	ELSE
	349	icap=1
	350	ENDIF
	351
	352	c
	353	c Loop on where we have condensation/ sublimation
	354	IF ((ztsrf(ig) .LT. ztcondsol(ig)) .OR. ! ground cond
	355	$ (zfallice(ig,1).NE.0.) .OR. ! falling snow
	356	$ ((ztsrf(ig) .GT. ztcondsol(ig)) .AND. ! ground sublim.
	357	$ ((piceco2(ig)+zfallice(ig,1)*ptimestep) .NE. 0.))) THEN
	358	condsub(ig) = .true.
	359
	360	IF (zfallice(ig,1).gt.0) then
	361	zfallheat=zfallice(ig,1)*
	362	& (pphi(ig,1)- phisfi(ig) +
[890]	363	& cpice*(ztcond(ig,1)-ztcondsol(ig)))/latcond
[38]	364	ELSE
	365	zfallheat=0.
	366	ENDIF
	367
	368	c condensation or partial sublimation of CO2 ice
	369	c """""""""""""""""""""""""""""""""""""""""""""""
	370	zcondices(ig)=pcapcal(ig)*(ztcondsol(ig)-ztsrf(ig))
	371	& /(latcond*ptimestep) - zfallheat
	372	pdtsrfc(ig) = (ztcondsol(ig) - ztsrf(ig))/ptimestep
	373
	374	c If the entire CO_2 ice layer sublimes
	375	c """"""""""""""""""""""""""""""""""""""""""""""""""""
	376	c (including what has just condensed in the atmosphere)
	377
	378	IF((piceco2(ig)/ptimestep+zfallice(ig,1)).LE.
	379	& -zcondices(ig))THEN
	380	zcondices(ig) = -piceco2(ig)/ptimestep - zfallice(ig,1)
	381	pdtsrfc(ig)=(latcond/pcapcal(ig))*
	382	& (zcondices(ig)+zfallheat)
	383	END IF
	384
	385	c Changing CO2 ice amount and pressure :
	386	c """"""""""""""""""""""""""""""""""""
	387
	388	zdiceco2(ig) = zcondices(ig) + zfallice(ig,1)
	389	piceco2(ig) = piceco2(ig) + zdiceco2(ig)*ptimestep
	390	pdpsrf(ig) = -zdiceco2(ig)*g
	391
	392	IF(ABS(pdpsrf(ig)*ptimestep).GT.pplev(ig,1)) THEN
	393	PRINT*,'STOP in condens'
	394	PRINT*,'condensing more than total mass'
	395	PRINT*,'Grid point ',ig
	396	PRINT*,'Ps = ',pplev(ig,1)
	397	PRINT*,'d Ps = ',pdpsrf(ig)
	398	STOP
	399	ENDIF
	400	END IF ! if there is condensation/sublimmation
	401	ENDDO ! of DO ig=1,ngrid
	402
	403	c ********************************************************************
	404	c Surface albedo and emissivity of the surface below the snow (emisref)
	405	c ********************************************************************
	406	c Prepare the case where albedo varies with insolation:
	407	c ----------------------------------------------------
	408	! if (fluxdependent) then
	409	!
	410	c Calcul du flux moyen (zonal mean)
	411	! do j=1,jjp1
	412	! Fluxmean(j)=0
	413	! do i=1,iim
	414	! ig=1+(j-2)*iim +i
	415	! if(j.eq.1) ig=1
	416	! if(j.eq.jjp1) ig=ngrid
	417	! Fluxmean(j)=Fluxmean(j)+fluxsurf_sw(ig,1)
	418	! $ +fluxsurf_sw(ig,2)
	419	! enddo
	420	! Fluxmean(j)=Fluxmean(j)/float(iim)
	421	! enddo
	422	!
	423	c const A and B used to calculate the albedo which depends on solar flux
	424	c albedice=constA+constB*Flux
	425	c constT = time step to calculate the solar flux when flux decreases
	426	! constA=0.26
	427	c constA=0.33
	428	c constA=0.186
	429	! constB=0.00187
	430	! constT=10
	431	! endif ! of if (fluxdependent)
	432
	433	! Check that amont of CO2 ice is not problematic
	434	DO ig=1,ngrid
	435	if(.not.piceco2(ig).ge.0.) THEN
	436	if(piceco2(ig).le.-5.e-8) print*,
	437	$ 'WARNING newcondens piceco2(',ig,')=', piceco2(ig)
	438	piceco2(ig)=0.
	439	endif
	440	ENDDO
	441
	442	! Set albedo and emissivity of the surface
	443	! ----------------------------------------
	444	CALL albedocaps(zls,ngrid,piceco2,psolaralb,emisref)
	445
	446	c Calcul de l'albedo
	447	c ------------------
	448	! do ig =1,ngrid
	449	! IF(ig.GT.ngrid/2+1) THEN
	450	! icap=2
	451	! ELSE
	452	! icap=1
	453	! ENDIF
	454	! IF(firstcall2) THEN
	455	! albediceF(ig)=albedice(icap)
	456	! ENDIF
	457	c if there is still co2ice ccccccccccccccccccccccc
	458	! if (piceco2(ig).gt.0) then
	459	! emisref(ig) = emisice(icap)
	460
	461	c if flux dependent albedo is used
	462	c --------------------------------
	463	! if (fluxdependent) then
	464	! j=INT((ig-2)/iim)+2
	465	! if(ig.eq.1) j=1
	466	! if(ig.eq.ngrid) j=jjp1
	467	c albediceF_new(ig)=MIN(constA+constB*Fluxmean(j),
	468	c $ constA+constB*250)
	469	! albediceF_new(ig)=constA+constB*Fluxmean(j)
	470	! if (albediceF(ig).gt.albediceF_new(ig)) then
	471	! albediceF(ig)=albediceF(ig)+ ptimestep/(daysec*
	472	! $ constT)*(albediceF_new(ig)-albediceF(ig))
	473	! else
	474	! albediceF(ig)=albediceF_new(ig)
	475	! endif
	476	c if part of the ice is transparent
	477	c slopy = pente de la droite: alpha = y*co2ice/1620
	478	c pour une valeur superieur a une epaisseur de glace donnee
	479	c ici, epaisseur limite = 10cm
	480	! if (transparency) then
	481	! slopy=1/(1620*0.10)
	482	! alpha=MIN(slopy*piceco2(ig),1.)
	483	! psolaralb(ig,1) = alpha*albediceF(ig)
	484	! $ +(1-alpha)*albedodat(ig)
	485	! psolaralb(ig,2) = psolaralb(ig,1)
	486	! else
	487	! psolaralb(ig,1) = albediceF(ig)
	488	! psolaralb(ig,2) = psolaralb(ig,1)
	489	! endif
	490	! else
	491	c transparency set to true and fluxdependent set to false
	492	! if (transparency) then
	493	! slopy=1/(1620*0.10)
	494	! alpha=MIN(slopy*piceco2(ig),1.)
	495	! psolaralb(ig,1) = alpha*albedice(icap)
	496	! $ +(1-alpha)*albedodat(ig)
	497	! psolaralb(ig,2) = psolaralb(ig,1)
	498	! else
	499	c simplest case: transparency and flux dependent set to false
	500	! psolaralb(ig,1) = albedice(icap)
	501	! psolaralb(ig,2) = albedice(icap)
	502	! endif
	503	! endif
	504	c no more co2ice, albedo = ground albedo
	505	! else
	506	! psolaralb(ig,1) = albedodat(ig)
	507	! psolaralb(ig,2) = albedodat(ig)
	508	! emisref(ig) = emissiv
	509	! pemisurf(ig) = emissiv
	510	! endif
	511	! end do ! end of the ig loop
	512
	513	! set pemisurf() to emissiv when there is bare surface (needed for co2snow)
	514	DO ig=1,ngrid
	515	if (piceco2(ig).eq.0) then
	516	pemisurf(ig)=emissiv
	517	endif
	518	ENDDO
	519
	520	! firstcall2=.false.
	521	c ***************************************************************
	522	c Correction to account for redistribution between sigma or hybrid
	523	c layers when changing surface pressure (and warming/cooling
	524	c of the CO2 currently changing phase).
	525	c *************************************************************
	526
	527	DO ig=1,ngrid
	528	if (condsub(ig)) then
	529	do l=1,nlayer
	530	ztc(l) =zt(ig,l) +pdtc(ig,l) *ptimestep
	531	zu(l) =pu(ig,l) +pdu( ig,l) *ptimestep
	532	zv(l) =pv(ig,l) +pdv( ig,l) *ptimestep
	533	do iq=1,nqmx
	534	zq(l,iq)=pq(ig,l,iq)+pdq(ig,l,iq)*ptimestep
	535	enddo
	536	end do
	537
	538	c Mass fluxes through the sigma levels (kg.m-2.s-1) (>0 when up)
	539	c """"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
	540
	541	zmflux(1) = -zcondices(ig)
	542	DO l=1,nlayer
	543	zmflux(l+1) = zmflux(l) -zcondicea(ig,l)
	544	& + (bp(l)-bp(l+1))*(zfallice(ig,1)-zmflux(1))
	545	c zmflux set to 0 if very low to avoid: top layer is disappearing in v1ld
	546	if (abs(zmflux(l+1)).lt.1E-13.OR.bp(l+1).eq.0.) zmflux(l+1)=0.
	547	END DO
	548
	549	c Mass of each layer
	550	c ------------------
	551	DO l=1,nlayer
	552	masse(l)=(pplev(ig,l) - pplev(ig,l+1))/g
	553	END DO
	554
	555
	556	c Corresponding fluxes for T,U,V,Q
	557	c """"""""""""""""""""""""""""""""
	558
	559	c averaging operator for TRANSPORT
	560	c """"""""""""""""""""""""""""""""
	561	c Value transfert at the surface interface when condensation
	562	c sublimation:
	563	ztm(1) = ztsrf(ig) + pdtsrfc(ig)*ptimestep
	564	zum(1) = 0
	565	zvm(1) = 0
	566	do iq=1,nqmx
	567	zqm(1,iq)=0. ! most tracer do not condense !
	568	enddo
	569	c Special case if one of the tracer is CO2 gas
	570	if (ico2.ne.0) zqm(1,ico2)=1. ! flux is 100% CO2
	571
	572	c Van Leer scheme:
	573	DO l=1,nlayer+1
	574	w(l)=-zmflux(l)*ptimestep
	575	END DO
	576	call vl1d(ztc,2.,masse,w,ztm)
	577	call vl1d(zu ,2.,masse,w,zum)
	578	call vl1d(zv ,2.,masse,w,zvm)
	579	do iq=1,nqmx
	580	do l=1,nlayer
	581	zq1(l)=zq(l,iq)
	582	enddo
	583	zqm1(1)=zqm(1,iq)
	584	call vl1d(zq1,2.,masse,w,zqm1)
	585	do l=2,nlayer
	586	zq( l,iq)=zq1(l)
	587	zqm(l,iq)=zqm1(l)
	588	enddo
	589	enddo
	590
	591	c Surface condensation affects low winds
	592	if (zmflux(1).lt.0) then
	593	zum(1)= zu(1) * (w(1)/masse(1))
	594	zvm(1)= zv(1) * (w(1)/masse(1))
	595	if (w(1).gt.masse(1)) then ! ensure numerical stability
	596	zum(1)= (zu(1)-zum(2))*masse(1)/w(1) + zum(2)
	597	zvm(1)= (zv(1)-zvm(2))*masse(1)/w(1) + zvm(2)
	598	end if
	599	end if
	600
	601	ztm(nlayer+1)= ztc(nlayer) ! should not be used, but...
	602	zum(nlayer+1)= zu(nlayer) ! should not be used, but...
	603	zvm(nlayer+1)= zv(nlayer) ! should not be used, but...
	604	do iq=1,nqmx
	605	zqm(nlayer+1,iq)= zq(nlayer,iq)
	606	enddo
[86]	607
	608	#ifdef MESOSCALE
	609	!!!! AS: This part must be commented in the mesoscale model
	610	!!!! AS: ... to avoid instabilities.
	611	!!!! AS: you have to compile with -DMESOSCALE to do so
	612	#else
[38]	613	c Tendencies on T, U, V, Q
	614	c """"""""""""""""""""""""
	615	DO l=1,nlayer
	616
	617	c Tendencies on T
	618	zdtsig(ig,l) = (1/masse(l)) *
	619	& ( zmflux(l)*(ztm(l) - ztc(l))
	620	& - zmflux(l+1)*(ztm(l+1) - ztc(l))
	621	& + zcondicea(ig,l)*(ztcond(ig,l)-ztc(l)) )
	622	pdtc(ig,l) = pdtc(ig,l) + zdtsig(ig,l)
	623
	624	c Tendencies on U
	625	pduc(ig,l) = (1/masse(l)) *
	626	& ( zmflux(l)*(zum(l) - zu(l))
	627	& - zmflux(l+1)*(zum(l+1) - zu(l)) )
	628
	629
	630	c Tendencies on V
	631	pdvc(ig,l) = (1/masse(l)) *
	632	& ( zmflux(l)*(zvm(l) - zv(l))
	633	& - zmflux(l+1)*(zvm(l+1) - zv(l)) )
	634
	635	END DO
[86]	636	#endif
[38]	637
	638	c Tendencies on Q
	639	do iq=1,nqmx
	640	! if (noms(iq).eq.'co2') then
	641	if (iq.eq.ico2) then
	642	c SPECIAL Case when the tracer IS CO2 :
	643	DO l=1,nlayer
	644	pdqc(ig,l,iq)= (1/masse(l)) *
	645	& ( zmflux(l)*(zqm(l,iq) - zq(l,iq))
	646	& - zmflux(l+1)*(zqm(l+1,iq) - zq(l,iq))
	647	& + zcondicea(ig,l)*(zq(l,iq)-1.) )
	648	END DO
	649	else
	650	DO l=1,nlayer
	651	pdqc(ig,l,iq)= (1/masse(l)) *
	652	& ( zmflux(l)*(zqm(l,iq) - zq(l,iq))
	653	& - zmflux(l+1)*(zqm(l+1,iq) - zq(l,iq))
	654	& + zcondicea(ig,l)*zq(l,iq) )
	655	END DO
	656	end if
	657	enddo
	658
	659	c --------------------------------------------------------
	660	c Roughly Simulate Molecular mixing when CO2 is too depleted by
	661	c Surface condensation (mixing starts if qco2 < qco2min )
	662	c FF 06/2004
	663	c WARNING : this is now done in convadj, better (FF 02/2005)
	664	c --------------------------------------------------------
	665	flag=0 ! now done in convadj : must be =0
	666	if (flag.eq.1) then
	667	if(ico2.gt.0) then ! relevant only if one tracer is CO2
	668	if(pq(ig,1,ico2)+(pdq(ig,1,ico2)+pdqc(ig,1,ico2))*ptimestep
	669	& .lt.qco2min) then
	670	do iq=1,nqmx
	671	zq(1,iq)=pq(ig,1,iq)
	672	& +(pdq(ig,1,iq)+pdqc(ig,1,iq))*ptimestep
	673	Smq(1,iq) = masse(1)*zq(1,iq)
	674	end do
	675	Sm(1) = masse(1)
	676	do l =2,nlayermx
	677	do iq=1,nqmx
	678	zq(l,iq)=pq(ig,l,iq)
	679	& +(pdq(ig,l,iq)+pdqc(ig,l,iq))*ptimestep
	680	smq(l,iq) = smq(l-1,iq) + masse(l)*zq(l,iq)
	681	end do
	682	sm(l) = sm(l-1) + masse(l)
	683	if(zq(l,ico2).gt.qco2min) then
	684	c mixmas: mass of atmosphere that must be mixed to reach qco2min
	685	mixmas = (sm(l-1)*zq(l,ico2)-Smq(l-1,ico2))
	686	& /(zq(l,ico2)-qco2min)
	687	if((mixmas.le.sm(l)))then
	688	c OK if mixed mass less than mass of layers affected
	689	nmix=l ! number of layer affected by mixing
	690	goto 99
	691	end if
	692	end if
	693	end do
	694	99 continue
	695	do iq=1,nqmx
	696	qmix=zq(nmix,iq)
	697	& +(Smq(nmix-1,iq)-zq(nmix,iq)*Sm(nmix-1))/mixmas
	698	do l=1,nmix-1
	699	pdqc(ig,l,iq)=
	700	& (qmix-pq(ig,l,iq))/ptimestep - pdq(ig,l,iq)
	701	end do
	702	c layer only partly mixed :
	703	pdqc(ig,nmix,iq)=(
	704	& qmix+(Sm(nmix)-mixmas)*(zq(nmix,iq)-qmix)/masse(nmix)
	705	& -pq(ig,nmix,iq))/ptimestep -pdq(ig,nmix,iq)
	706
	707	end do
	708	end if
	709	end if
	710
	711	endif ! (flag.eq.1)
	712	end if ! if (condsub)
	713	END DO ! loop on ig
	714
	715	c ***************************************************************
	716	c CO2 snow / clouds scheme
	717	c ***************************************************************
	718
	719	call co2snow(ngrid,nlayer,ptimestep,emisref,condsub,pplev,
	720	& zcondicea,zcondices,zfallice,pemisurf)
	721
	722	c ***************************************************************
	723	c Ecriture des diagnostiques
	724	c ***************************************************************
	725
	726	c DO l=1,nlayer
	727	c DO ig=1,ngrid
	728	c Taux de cond en kg.m-2.pa-1.s-1
	729	c tconda1(ig,l)=zcondicea(ig,l)/(pplev(ig,l)-pplev(ig,l+1))
	730	c Taux de cond en kg.m-3.s-1
	731	c tconda2(ig,l)=tconda1(ig,l)pplay(ig,l)g/(r*pt(ig,l))
	732	c END DO
	733	c END DO
	734	c call WRITEDIAGFI(ngridmx,'tconda1',
	735	c &'Taux de condensation CO2 atmospherique /Pa',
	736	c & 'kg.m-2.Pa-1.s-1',3,tconda1)
	737	c call WRITEDIAGFI(ngridmx,'tconda2',
	738	c &'Taux de condensation CO2 atmospherique /m',
	739	c & 'kg.m-3.s-1',3,tconda2)
	740
	741	! output falling co2 ice in 1st layer:
	742	! call WRITEDIAGFI(ngridmx,'fallice',
	743	! &'Precipitation of co2 ice',
	744	! & 'kg.m-2.s-1',2,zfallice(1,1))
	745
	746	!! Specific stuff to bound co2 tracer values ....
	747	if (bound_qco2.and.(ico2.ne.0)) then
[890]	748	do ig=1,ngrid
	749	do l=1,nlayer
[38]	750	zqco2=pq(ig,l,ico2)
	751	& +(pdq(ig,l,ico2)+pdqc(ig,l,ico2))*ptimestep
	752	if (zqco2.gt.qco2max) then
	753	! correct pdqc:
	754	pdqc(ig,l,ico2)=((qco2max-pq(ig,l,ico2))/ptimestep)
	755	& -pdq(ig,l,ico2)
	756	write(,) "newcondens: adapting pdqc(ig,l,ico2)",
	757	& " so that co2 conc. does not exceed",qco2max
	758	write(,) " ig:",ig," l:",l
	759	endif ! of if (zqco2.gt.qco2max)
	760	if (zqco2.lt.qco2mini) then
	761	! correct pdqc:
	762	pdqc(ig,l,ico2)=((qco2mini-pq(ig,l,ico2))/ptimestep)
	763	& -pdq(ig,l,ico2)
	764	write(,) "newcondens: adapting pdqc(ig,l,ico2)",
	765	& " so that co2 conc. is not less than",qco2mini
	766	write(,) " ig:",ig," l:",l
	767	endif ! of if (zqco2.lt.qco2mini)
	768	end do
	769	enddo
	770	endif ! of if (bound_qco2.and.(ico2.ne.0)) then
	771
	772	return
	773	end
	774
	775
	776
	777	c *****************************************************************
	778	SUBROUTINE vl1d(q,pente_max,masse,w,qm)
	779	c
	780	c
	781	c Operateur de moyenne inter-couche pour calcul de transport type
	782	c Van-Leer " pseudo amont " dans la verticale
	783	c q,w sont des arguments d'entree pour le s-pg ....
	784	c masse : masse de la couche Dp/g
	785	c w : masse d'atm ``transferee'' a chaque pas de temps (kg.m-2)
	786	c pente_max = 2 conseillee
	787	c
	788	c
	789	c --------------------------------------------------------------------
	790	IMPLICIT NONE
	791
	792	#include "dimensions.h"
	793
	794	c
	795	c
	796	c
	797	c Arguments:
	798	c ----------
	799	real masse(llm),pente_max
	800	REAL q(llm),qm(llm+1)
	801	REAL w(llm+1)
	802	c
	803	c Local
	804	c ---------
	805	c
	806	INTEGER l
	807	c
	808	real dzq(llm),dzqw(llm),adzqw(llm),dzqmax
	809	real sigw, Mtot, MQtot
	810	integer m
	811	c integer ismax,ismin
	812
	813
	814	c On oriente tout dans le sens de la pression
	815	c W > 0 WHEN DOWN !!!!!!!!!!!!!
	816
	817	do l=2,llm
	818	dzqw(l)=q(l-1)-q(l)
	819	adzqw(l)=abs(dzqw(l))
	820	enddo
	821
	822	do l=2,llm-1
	823	if(dzqw(l)*dzqw(l+1).gt.0.) then
	824	dzq(l)=0.5*(dzqw(l)+dzqw(l+1))
	825	else
	826	dzq(l)=0.
	827	endif
	828	dzqmax=pente_max*min(adzqw(l),adzqw(l+1))
	829	dzq(l)=sign(min(abs(dzq(l)),dzqmax),dzq(l))
	830	enddo
	831
	832	dzq(1)=0.
	833	dzq(llm)=0.
	834
	835	do l = 1,llm-1
	836
	837	c Regular scheme (transfered mass < layer mass)
	838	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	839	if(w(l+1).gt.0. .and. w(l+1).le.masse(l+1)) then
	840	sigw=w(l+1)/masse(l+1)
	841	qm(l+1)=(q(l+1)+0.5(1.-sigw)dzq(l+1))
	842	else if(w(l+1).le.0. .and. -w(l+1).le.masse(l)) then
	843	sigw=w(l+1)/masse(l)
	844	qm(l+1)=(q(l)-0.5(1.+sigw)dzq(l))
	845
	846	c Extended scheme (transfered mass > layer mass)
	847	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	848	else if(w(l+1).gt.0.) then
	849	m=l+1
	850	Mtot = masse(m)
	851	MQtot = masse(m)*q(m)
	852	do while ((m.lt.llm).and.(w(l+1).gt.(Mtot+masse(m+1))))
	853	m=m+1
	854	Mtot = Mtot + masse(m)
	855	MQtot = MQtot + masse(m)*q(m)
	856	end do
	857	if (m.lt.llm) then
	858	sigw=(w(l+1)-Mtot)/masse(m+1)
	859	qm(l+1)= (1/w(l+1))(MQtot + (w(l+1)-Mtot)
	860	& (q(m+1)+0.5(1.-sigw)dzq(m+1)) )
	861	else
	862	w(l+1) = Mtot
	863	qm(l+1) = Mqtot / Mtot
	864	write(,) 'top layer is disapearing !'
	865	stop
	866	end if
	867	else ! if(w(l+1).lt.0)
	868	m = l-1
	869	Mtot = masse(m+1)
	870	MQtot = masse(m+1)*q(m+1)
[120]	871	if (m.gt.0) then ! because some compilers will have problems
	872	! evaluating masse(0)
	873	do while ((m.gt.0).and.(-w(l+1).gt.(Mtot+masse(m))))
[38]	874	m=m-1
	875	Mtot = Mtot + masse(m+1)
	876	MQtot = MQtot + masse(m+1)*q(m+1)
[120]	877	if (m.eq.0) exit
	878	end do
	879	endif
[38]	880	if (m.gt.0) then
	881	sigw=(w(l+1)+Mtot)/masse(m)
	882	qm(l+1)= (-1/w(l+1))(MQtot + (-w(l+1)-Mtot)
	883	& (q(m)-0.5(1.+sigw)dzq(m)) )
	884	else
	885	qm(l+1)= (-1/w(l+1))(MQtot + (-w(l+1)-Mtot)qm(1))
	886	end if
	887	end if
	888	enddo
	889
	890	c boundary conditions (not used in newcondens !!)
	891	c qm(llm+1)=0.
	892	c if(w(1).gt.0.) then
	893	c qm(1)=q(1)
	894	c else
	895	c qm(1)=0.
	896	c end if
	897
	898	return
	899	end

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

Download in other formats:

Original Format