Context Navigation

newcondens.F @ 3574

Last change on this file since 3574 was 1867, checked in by aslmd, 7 years ago
in mesoscale old physics, corrected runtime errors with picky compilers such as gfortran
File size: 26.4 KB

Rev	Line
[57]	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)
	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	#include "fisice.h"
	66
	67	c-----------------------------------------------------------------------
	68	c Arguments :
	69	c ---------
	70	INTEGER ngrid, nlayer, nq
	71
	72	REAL ptimestep
	73	REAL pplay(ngrid,nlayer),pplev(ngrid,nlayer+1)
	74	REAL pcapcal(ngrid)
	75	REAL pt(ngrid,nlayer)
	76	REAL ptsrf(ngrid)
	77	REAL pphi(ngrid,nlayer)
	78	REAL pdt(ngrid,nlayer),pdtsrf(ngrid),pdtc(ngrid,nlayer)
	79	REAL pdtsrfc(ngrid),pdpsrf(ngrid)
	80	REAL piceco2(ngrid),psolaralb(ngrid,2),pemisurf(ngrid)
	81
	82	REAL pu(ngrid,nlayer) , pv(ngrid,nlayer)
	83	REAL pdu(ngrid,nlayer) , pdv(ngrid,nlayer)
	84	REAL pduc(ngrid,nlayer) , pdvc(ngrid,nlayer)
	85	REAL pq(ngridmx,nlayer,nq),pdq(ngrid,nlayer,nq)
	86	REAL pdqc(ngrid,nlayer,nq)
	87	REAL fluxsurf_sw(ngrid,2) ! added to calculate flux dependent albedo
	88	c
	89	c Local variables :
	90	c -----------------
	91
	92	c variables used for albedo parametrization
	93	c --------------------------------------------
	94	INTEGER i,j
	95	REAL Fluxmean(jjp1)
	96	INTEGER l,ig,iq,icap,nmix
	97	LOGICAL transparency, fluxdependent
	98	c flag transparency if you want to make the co2ice semi-transparent
	99	PARAMETER(transparency=.true.)
	100	c flag fluxdependent if you want the co2ice albedo to be dependent on
	101	c the incident solar flux
	102	PARAMETER(fluxdependent=.false.)
	103	REAL slopy,alpha,constA,constB,constT,albediceF_new(ngridmx)
	104	REAL zt(ngridmx,nlayermx)
	105	REAL zcpi
	106	REAL ztcond (ngridmx,nlayermx)
	107	REAL ztcondsol(ngridmx)
	108	REAL zdiceco2(ngridmx)
	109	c REAL zcondicea(ngridmx,nlayermx) ! Already defined in fisice.h
	110	REAL zcondices(ngridmx)
	111	REAL zfallice(ngridmx,nlayermx+1) , zfallheat
	112	REAL zmflux(nlayermx+1)
	113	REAL zu(nlayermx),zv(nlayermx)
	114	REAL zq(nlayermx,nqmx),zq1(nlayermx)
	115	REAL ztsrf(ngridmx)
	116	REAL ztc(nlayermx), ztm(nlayermx+1)
	117	REAL zum(nlayermx+1) , zvm(nlayermx+1)
	118	REAL zqm(nlayermx+1,nqmx),zqm1(nlayermx+1)
	119	REAL masse(nlayermx),w(nlayermx+1)
	120	REAL Sm(nlayermx),Smq(nlayermx,nqmx),mixmas,qmix
	121	LOGICAL condsub(ngridmx)
	122
	123	c variable speciale diagnostique
	124	real tconda1(ngridmx,nlayermx)
	125	real tconda2(ngridmx,nlayermx)
	126	c REAL zdiceco2a(ngridmx) ! for diagnostic only
	127	real zdtsig (ngridmx,nlayermx)
	128	real zdt (ngridmx,nlayermx)
	129
	130
	131	c local saved variables
	132	integer ico2
	133	real qco2min
	134	save ico2,qco2min
	135	REAL emisref(ngridmx)
	136	REAL latcond,tcond1mb
	137	REAL acond,bcond,ccond,cpice
	138	REAL albediceF(ngridmx)
	139	SAVE emisref
	140	SAVE latcond,acond,bcond,ccond,cpice
	141	SAVE albediceF
	142
	143
	144	LOGICAL firstcall,firstcall2
	145	SAVE firstcall,firstcall2
	146	REAL SSUM
	147	EXTERNAL SSUM
	148
	149	c common/scratch/zt,ztcond,zcondicea,zfallice,tconda1
	150	c , ,tconda2,zdtsig,zdt,zu,zv
	151
	152
	153	DATA latcond,tcond1mb/5.9e5,136.27/
	154	DATA cpice /1000./
	155	DATA firstcall/.true./
	156	DATA firstcall2/.true./
	157
	158	integer flag
	159
	160	c----------------------------------------------------------------------
	161
	162	c Initialisation
	163	c --------------
	164	c
	165	IF (firstcall) THEN
	166	bcond=1./tcond1mb
	167	ccond=cpp/(g*latcond)
	168	acond=r/latcond
	169
	170	firstcall=.false.
	171	PRINT*,'In newcondens:Tcond(P=1mb)=',tcond1mb,' Lcond=',latcond
	172	PRINT*,'acond,bcond,ccond',acond,bcond,acond
	173
	174	ico2=0
	175	if (tracer) then
	176	c Prepare Special treatment if one of the tracer is CO2 gas
	177	do iq=1,nqmx
	178	if (noms(iq).eq."co2") ico2=iq
	179	enddo
	180	c minimum CO2 mix. ratio below which mixing occur with layer above:
	181	qco2min =0.75
	182	end if
	183	ENDIF
	184	zcpi=1./cpp
	185	c
	186	c======================================================================
	187	c Calcul of CO2 condensation sublimation
	188	c ============================================================
	189	c
	190	c Used variable :
	191	c piceco2(ngrid) : amount of co2 ice on the ground (kg/m2)
	192	c zcondicea(ngrid,l): condensation rate in layer l (kg/m2/s)
	193	c zcondices(ngrid): condensation rate on the ground (kg/m2/s)
	194	c zfallice(ngrid,l):amount of ice falling from layer l (kg/m2/s)
	195	c
	196	c pdtc(ngrid,nlayermx) : dT/dt due to cond/sub
	197	c
	198	c
	199	c Tendencies set to 0 (except pdtc)
	200	c -------------------------------------
	201	DO l=1,nlayer
	202	DO ig=1,ngrid
	203	zcondicea(ig,l) = 0.
	204	zfallice(ig,l) = 0.
	205	pduc(ig,l) = 0
	206	pdvc(ig,l) = 0
	207	END DO
	208	END DO
	209
	210	DO iq=1,nqmx
	211	DO l=1,nlayer
	212	DO ig=1,ngrid
	213	pdqc(ig,l,iq) = 0
	214	END DO
	215	END DO
	216	END DO
	217
	218	DO ig=1,ngrid
	219	zfallice(ig,nlayer+1) = 0.
	220	zcondices(ig) = 0.
	221	pdtsrfc(ig) = 0.
	222	pdpsrf(ig) = 0.
	223	condsub(ig) = .false.
	224	zdiceco2(ig) = 0.
	225	ENDDO
	226	zfallheat=0
	227
	228	c *************************
	229	c ATMOSPHERIC CONDENSATION
	230	c *************************
	231
	232	c forecast of atmospheric temperature zt and frost temperature ztcond
	233	c --------------------------------------------------------------------
	234
	235	DO l=1,nlayer
	236	DO ig=1,ngrid
	237	zt(ig,l)=pt(ig,l)+ pdt(ig,l)*ptimestep
	238	ztcond(ig,l)=1./(bcond-acondlog(.0095pplay(ig,l)))
	239	if (pplay(ig,l).lt.1e-4) ztcond(ig,l)=0.0 !mars Monica
	240	ENDDO
	241	ENDDO
	242
	243	c Condensation/sublimation in the atmosphere
	244	c ------------------------------------------
	245	c (calcul of zcondicea , zfallice and pdtc)
	246	c
	247	DO l=nlayer , 1, -1
	248	DO ig=1,ngrid
	249	pdtc(ig,l)=0.
	250	IF((zt(ig,l).LT.ztcond(ig,l)).or.(zfallice(ig,l+1).gt.0))THEN
	251	condsub(ig)=.true.
	252	IF (zfallice(ig,l+1).gt.0) then
	253	zfallheat=zfallice(ig,l+1)*
	254	& (pphi(ig,l+1)-pphi(ig,l) +
	255	& cpice*(ztcond(ig,l+1)-ztcond(ig,l)))/latcond
	256	ELSE
	257	zfallheat=0.
	258	ENDIF
	259	pdtc(ig,l)=(ztcond(ig,l) - zt(ig,l))/ptimestep
	260	zcondicea(ig,l)=(pplev(ig,l)-pplev(ig,l+1))
	261	& ccondpdtc(ig,l)- zfallheat
	262	c Case when the ice from above sublimes entirely
	263	c """""""""""""""""""""""""""""""""""""""""""""""
	264	IF (zfallice(ig,l+1).lt.- zcondicea(ig,l)) then
	265	pdtc(ig,l)=(-zfallice(ig,l+1)+zfallheat)/
	266	& (ccond*(pplev(ig,l)-pplev(ig,l+1)))
	267	zcondicea(ig,l)= -zfallice(ig,l+1)
	268	END IF
	269
	270	zfallice(ig,l) = zcondicea(ig,l)+zfallice(ig,l+1)
	271	END IF
	272	ENDDO
	273	ENDDO
	274
	275	c *************************
	276	c SURFACE CONDENSATION
	277	c *************************
	278
	279	c forecast of ground temperature ztsrf and frost temperature ztcondsol
	280	c --------------------------------------------------------------------
	281	DO ig=1,ngrid
	282	ztcondsol(ig)=1./(bcond-acondlog(.0095pplev(ig,1)))
	283	ztsrf(ig) = ptsrf(ig) + pdtsrf(ig)*ptimestep
	284	ENDDO
	285
	286	c
	287	c Condensation/sublimation on the ground
	288	c --------------------------------------
	289	c (calcul of zcondices , pdtsrfc)
	290	c
	291	DO ig=1,ngrid
	292	IF(ig.GT.ngrid/2+1) THEN
	293	icap=2
	294	ELSE
	295	icap=1
	296	ENDIF
	297
	298	c
	299	c Loop on where we have condensation/ sublimation
	300	IF ((ztsrf(ig) .LT. ztcondsol(ig)) .OR. ! ground cond
	301	$ (zfallice(ig,1).NE.0.) .OR. ! falling snow
	302	$ ((ztsrf(ig) .GT. ztcondsol(ig)) .AND. ! ground sublim.
	303	$ ((piceco2(ig)+zfallice(ig,1)*ptimestep) .NE. 0.))) THEN
	304	condsub(ig) = .true.
	305
	306	IF (zfallice(ig,1).gt.0) then
	307	zfallheat=zfallice(ig,1)*
	308	& (pphi(ig,1)- phisfi(ig) +
	309	& cpice(ztcond(ig,1)-ztcondsol(ig)))/(latcondptimestep)
	310	ELSE
	311	zfallheat=0.
	312	ENDIF
	313
	314	c condensation or partial sublimation of CO2 ice
	315	c """""""""""""""""""""""""""""""""""""""""""""""
	316	zcondices(ig)=pcapcal(ig)*(ztcondsol(ig)-ztsrf(ig))
	317	& /(latcond*ptimestep) - zfallheat
	318	pdtsrfc(ig) = (ztcondsol(ig) - ztsrf(ig))/ptimestep
	319
	320	c If the entire CO_2 ice layer sublimes
	321	c """"""""""""""""""""""""""""""""""""""""""""""""""""
	322	c (including what has just condensed in the atmosphere)
	323
	324	IF((piceco2(ig)/ptimestep+zfallice(ig,1)).LE.
	325	& -zcondices(ig))THEN
	326	zcondices(ig) = -piceco2(ig)/ptimestep - zfallice(ig,1)
	327	pdtsrfc(ig)=(latcond/pcapcal(ig))*
	328	& (zcondices(ig)+zfallheat)
	329	END IF
	330
	331	c Changing CO2 ice amount and pressure :
	332	c """"""""""""""""""""""""""""""""""""
	333
	334	zdiceco2(ig) = zcondices(ig) + zfallice(ig,1)
	335	piceco2(ig) = piceco2(ig) + zdiceco2(ig)*ptimestep
	336	pdpsrf(ig) = -zdiceco2(ig)*g
	337
	338	IF(ABS(pdpsrf(ig)*ptimestep).GT.pplev(ig,1)) THEN
	339	PRINT*,'STOP in condens'
	340	PRINT*,'condensing more than total mass'
	341	PRINT*,'Grid point ',ig
	342	PRINT*,'Ps = ',pplev(ig,1)
	343	PRINT*,'d Ps = ',pdpsrf(ig)
	344	STOP
	345	ENDIF
	346	END IF
	347	ENDDO
	348
	349	c ********************************************************************
	350	c Surface albedo and emissivity of the surface below the snow (emisref)
	351	c ********************************************************************
	352	c Prepare the case where albedo varies with insolation:
	353	c ----------------------------------------------------
	354	if (fluxdependent) then
	355
	356	c Calcul du flux moyen (zonal mean)
	357	do j=1,jjp1
	358	Fluxmean(j)=0
	359	do i=1,iim
	360	ig=1+(j-2)*iim +i
	361	if(j.eq.1) ig=1
	362	if(j.eq.jjp1) ig=ngrid
	363	Fluxmean(j)=Fluxmean(j)+fluxsurf_sw(ig,1)
	364	$ +fluxsurf_sw(ig,2)
	365	enddo
	366	Fluxmean(j)=Fluxmean(j)/float(iim)
	367	enddo
	368
	369	c const A and B used to calculate the albedo which depends on solar flux
	370	c albedice=constA+constB*Flux
	371	c constT = time step to calculate the solar flux when flux decreases
	372	constA=0.26
	373	c constA=0.33
	374	c constA=0.186
	375	constB=0.00187
	376	constT=10
	377	endif
	378
	379	c Calcul de l'albedo
	380	c ------------------
	381	do ig =1,ngrid
	382	IF(ig.GT.ngrid/2+1) THEN
	383	icap=2
	384	ELSE
	385	icap=1
	386	ENDIF
	387	IF(firstcall2) THEN
	388	albediceF(ig)=albedice(icap)
	389	ENDIF
	390	c test on the existence of co2ice ccccccccccccccccccccccc
	391	if(.not.piceco2(ig).ge.0.) THEN
	392	if(piceco2(ig).le.-5.e-8) print*,
	393	$ 'WARNING newcondens piceco2(',ig,')=', piceco2(ig)
	394	piceco2(ig)=0.
	395	endif
	396	c if there is still co2ice ccccccccccccccccccccccc
	397	if (piceco2(ig).gt.0) then
	398	emisref(ig) = emisice(icap)
	399
	400	c if flux dependent albedo is used
	401	c --------------------------------
	402	if (fluxdependent) then
	403	j=INT((ig-2)/iim)+2
	404	if(ig.eq.1) j=1
	405	if(ig.eq.ngrid) j=jjp1
	406	c albediceF_new(ig)=MIN(constA+constB*Fluxmean(j),
	407	c $ constA+constB*250)
	408	albediceF_new(ig)=constA+constB*Fluxmean(j)
	409	if (albediceF(ig).gt.albediceF_new(ig)) then
	410	albediceF(ig)=albediceF(ig)+ ptimestep/(daysec*
	411	$ constT)*(albediceF_new(ig)-albediceF(ig))
	412	else
	413	albediceF(ig)=albediceF_new(ig)
	414	endif
	415	c if part of the ice is transparent
	416	c slopy = pente de la droite: alpha = y*co2ice/1620
	417	c pour une valeur superieur a une epaisseur de glace donnee
	418	c ici, epaisseur limite = 10cm
	419	if (transparency) then
	420	slopy=1/(1620*0.10)
	421	alpha=MIN(slopy*piceco2(ig),1.)
	422	psolaralb(ig,1) = alpha*albediceF(ig)
	423	$ +(1-alpha)*albedodat(ig)
	424	psolaralb(ig,2) = psolaralb(ig,1)
	425	else
	426	psolaralb(ig,1) = albediceF(ig)
	427	psolaralb(ig,2) = psolaralb(ig,1)
	428	endif
	429	else
	430	c transparency set to true and fluxdependent set to false
	431	if (transparency) then
	432	slopy=1/(1620*0.10)
	433	alpha=MIN(slopy*piceco2(ig),1.)
	434	psolaralb(ig,1) = alpha*albedice(icap)
	435	$ +(1-alpha)*albedodat(ig)
	436	psolaralb(ig,2) = psolaralb(ig,1)
	437	else
	438	c simplest case: transparency and flux dependent set to false
	439	psolaralb(ig,1) = albedice(icap)
	440	psolaralb(ig,2) = albedice(icap)
	441	endif
	442	endif
	443	c no more co2ice, albedo = ground albedo
	444	else
	445	psolaralb(ig,1) = albedodat(ig)
	446	psolaralb(ig,2) = albedodat(ig)
	447	emisref(ig) = emissiv
	448	pemisurf(ig) = emissiv
	449	endif
	450	end do ! end of the ig loop
	451
	452	firstcall2=.false.
	453	c ***************************************************************
	454	c Correction to account for redistribution between sigma or hybrid
	455	c layers when changing surface pressure (and warming/cooling
	456	c of the CO2 currently changing phase).
	457	c *************************************************************
	458
	459	DO ig=1,ngrid
	460	if (condsub(ig)) then
	461	do l=1,nlayer
	462	ztc(l) =zt(ig,l) +pdtc(ig,l) *ptimestep
	463	zu(l) =pu(ig,l) +pdu( ig,l) *ptimestep
	464	zv(l) =pv(ig,l) +pdv( ig,l) *ptimestep
	465	do iq=1,nqmx
	466	zq(l,iq)=pq(ig,l,iq)+pdq(ig,l,iq)*ptimestep
	467	enddo
	468	end do
	469
	470	c Mass fluxes through the sigma levels (kg.m-2.s-1) (>0 when up)
	471	c """"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
	472
	473	zmflux(1) = -zcondices(ig)
	474	DO l=1,nlayer
	475	zmflux(l+1) = zmflux(l) -zcondicea(ig,l)
	476	& + (bp(l)-bp(l+1))*(zfallice(ig,1)-zmflux(1))
	477	c zmflux set to 0 if very low to avoid: top layer is disappearing in v1ld
	478	if (abs(zmflux(l+1)).lt.1E-13.OR.bp(l+1).eq.0.) zmflux(l+1)=0.
	479	END DO
	480
	481	c Mass of each layer
	482	c ------------------
	483	DO l=1,nlayer
	484	masse(l)=(pplev(ig,l) - pplev(ig,l+1))/g
	485	END DO
	486
	487
	488	c Corresponding fluxes for T,U,V,Q
	489	c """"""""""""""""""""""""""""""""
	490
	491	c averaging operator for TRANSPORT
	492	c """"""""""""""""""""""""""""""""
	493	c Value transfert at the surface interface when condensation
	494	c sublimation:
	495	ztm(1) = ztsrf(ig) + pdtsrfc(ig)*ptimestep
	496	zum(1) = 0
	497	zvm(1) = 0
	498	do iq=1,nqmx
	499	zqm(1,iq)=0. ! most tracer do not condense !
	500	enddo
	501	c Special case if one of the tracer is CO2 gas
	502	if (ico2.ne.0) zqm(1,ico2)=1. ! flux is 100% CO2
	503
	504	c Van Leer scheme:
	505	DO l=1,nlayer+1
	506	w(l)=-zmflux(l)*ptimestep
	507	END DO
	508	call vl1d(ztc,2.,masse,w,ztm)
	509	call vl1d(zu ,2.,masse,w,zum)
	510	call vl1d(zv ,2.,masse,w,zvm)
	511	do iq=1,nqmx
	512	do l=1,nlayer
	513	zq1(l)=zq(l,iq)
	514	enddo
	515	zqm1(1)=zqm(1,iq)
	516	call vl1d(zq1,2.,masse,w,zqm1)
	517	do l=2,nlayer
	518	zq( l,iq)=zq1(l)
	519	zqm(l,iq)=zqm1(l)
	520	enddo
	521	enddo
	522
	523	c Surface condensation affects low winds
	524	if (zmflux(1).lt.0) then
	525	zum(1)= zu(1) * (w(1)/masse(1))
	526	zvm(1)= zv(1) * (w(1)/masse(1))
	527	if (w(1).gt.masse(1)) then ! ensure numerical stability
	528	zum(1)= (zu(1)-zum(2))*masse(1)/w(1) + zum(2)
	529	zvm(1)= (zv(1)-zvm(2))*masse(1)/w(1) + zvm(2)
	530	end if
	531	end if
	532
	533	ztm(nlayer+1)= ztc(nlayer) ! should not be used, but...
	534	zum(nlayer+1)= zu(nlayer) ! should not be used, but...
	535	zvm(nlayer+1)= zv(nlayer) ! should not be used, but...
	536	do iq=1,nqmx
	537	zqm(nlayer+1,iq)= zq(nlayer,iq)
	538	enddo
	539
	540	c Tendencies on T, U, V, Q
	541	c """"""""""""""""""""""""
	542	DO l=1,nlayer
	543
	544	c Tendencies on T
	545	zdtsig(ig,l) = (1/masse(l)) *
	546	& ( zmflux(l)*(ztm(l) - ztc(l))
	547	& - zmflux(l+1)*(ztm(l+1) - ztc(l))
	548	& + zcondicea(ig,l)*(ztcond(ig,l)-ztc(l)) )
	549	pdtc(ig,l) = pdtc(ig,l) + zdtsig(ig,l)
	550
	551	c Tendencies on U
	552	pduc(ig,l) = (1/masse(l)) *
	553	& ( zmflux(l)*(zum(l) - zu(l))
	554	& - zmflux(l+1)*(zum(l+1) - zu(l)) )
	555
	556
	557	c Tendencies on V
	558	pdvc(ig,l) = (1/masse(l)) *
	559	& ( zmflux(l)*(zvm(l) - zv(l))
	560	& - zmflux(l+1)*(zvm(l+1) - zv(l)) )
	561
	562	END DO
	563
	564	c Tendencies on Q
	565	do iq=1,nqmx
	566	if (noms(iq).eq.'co2') then
	567	c SPECIAL Case when the tracer IS CO2 :
	568	DO l=1,nlayer
	569	pdqc(ig,l,iq)= (1/masse(l)) *
	570	& ( zmflux(l)*(zqm(l,iq) - zq(l,iq))
	571	& - zmflux(l+1)*(zqm(l+1,iq) - zq(l,iq))
	572	& + zcondicea(ig,l)*(zq(l,iq)-1.) )
	573	END DO
	574	else
	575	DO l=1,nlayer
	576	pdqc(ig,l,iq)= (1/masse(l)) *
	577	& ( zmflux(l)*(zqm(l,iq) - zq(l,iq))
	578	& - zmflux(l+1)*(zqm(l+1,iq) - zq(l,iq))
	579	& + zcondicea(ig,l)*zq(l,iq) )
	580	END DO
	581	end if
	582	enddo
	583
	584	c --------------------------------------------------------
	585	c Roughly Simulate Molecular mixing when CO2 is too depleted by
	586	c Surface condensation (mixing starts if qco2 < qco2min )
	587	c FF 06/2004
	588	c WARNING : this is now done in convadj, better (FF 02/2005)
	589	c --------------------------------------------------------
	590	flag=0 ! now done in convadj : must be =0
	591	if (flag.eq.1) then
	592	if(ico2.gt.0) then ! relevant only if one tracer is CO2
	593	if(pq(ig,1,ico2)+(pdq(ig,1,ico2)+pdqc(ig,1,ico2))*ptimestep
	594	& .lt.qco2min) then
	595	do iq=1,nqmx
	596	zq(1,iq)=pq(ig,1,iq)
	597	& +(pdq(ig,1,iq)+pdqc(ig,1,iq))*ptimestep
	598	Smq(1,iq) = masse(1)*zq(1,iq)
	599	end do
	600	Sm(1) = masse(1)
	601	do l =2,nlayermx
	602	do iq=1,nqmx
	603	zq(l,iq)=pq(ig,l,iq)
	604	& +(pdq(ig,l,iq)+pdqc(ig,l,iq))*ptimestep
	605	smq(l,iq) = smq(l-1,iq) + masse(l)*zq(l,iq)
	606	end do
	607	sm(l) = sm(l-1) + masse(l)
	608	if(zq(l,ico2).gt.qco2min) then
	609	c mixmas: mass of atmosphere that must be mixed to reach qco2min
	610	mixmas = (sm(l-1)*zq(l,ico2)-Smq(l-1,ico2))
	611	& /(zq(l,ico2)-qco2min)
	612	if((mixmas.le.sm(l)))then
	613	c OK if mixed mass less than mass of layers affected
	614	nmix=l ! number of layer affected by mixing
	615	goto 99
	616	end if
	617	end if
	618	end do
	619	99 continue
	620	do iq=1,nqmx
	621	qmix=zq(nmix,iq)
	622	& +(Smq(nmix-1,iq)-zq(nmix,iq)*Sm(nmix-1))/mixmas
	623	do l=1,nmix-1
	624	pdqc(ig,l,iq)=
	625	& (qmix-pq(ig,l,iq))/ptimestep - pdq(ig,l,iq)
	626	end do
	627	c layer only partly mixed :
	628	pdqc(ig,nmix,iq)=(
	629	& qmix+(Sm(nmix)-mixmas)*(zq(nmix,iq)-qmix)/masse(nmix)
	630	& -pq(ig,nmix,iq))/ptimestep -pdq(ig,nmix,iq)
	631
	632	end do
	633	end if
	634	end if
	635
	636	endif ! (flag.eq.1)
	637	end if ! if (condsub)
	638	END DO ! loop on ig
	639
	640	c ***************************************************************
	641	c CO2 snow / clouds scheme
	642	c ***************************************************************
	643
	644	call co2snow(ngrid,nlayer,ptimestep,emisref,condsub,pplev,
	645	& zcondicea,zcondices,zfallice,pemisurf)
	646
	647	c ***************************************************************
	648	c Ecriture des diagnostiques
	649	c ***************************************************************
	650
	651	c DO l=1,nlayer
	652	c DO ig=1,ngrid
	653	c Taux de cond en kg.m-2.pa-1.s-1
	654	c tconda1(ig,l)=zcondicea(ig,l)/(pplev(ig,l)-pplev(ig,l+1))
	655	c Taux de cond en kg.m-3.s-1
	656	c tconda2(ig,l)=tconda1(ig,l)pplay(ig,l)g/(r*pt(ig,l))
	657	c END DO
	658	c END DO
	659	c call WRITEDIAGFI(ngridmx,'tconda1',
	660	c &'Taux de condensation CO2 atmospherique /Pa',
	661	c & 'kg.m-2.Pa-1.s-1',3,tconda1)
	662	c call WRITEDIAGFI(ngridmx,'tconda2',
	663	c &'Taux de condensation CO2 atmospherique /m',
	664	c & 'kg.m-3.s-1',3,tconda2)
	665
	666	return
	667	end
	668
	669
	670
	671	c *****************************************************************
	672	SUBROUTINE vl1d(q,pente_max,masse,w,qm)
	673	c
	674	c
	675	c Operateur de moyenne inter-couche pour calcul de transport type
	676	c Van-Leer " pseudo amont " dans la verticale
	677	c q,w sont des arguments d'entree pour le s-pg ....
	678	c masse : masse de la couche Dp/g
	679	c w : masse d'atm ``transferee'' a chaque pas de temps (kg.m-2)
	680	c pente_max = 2 conseillee
	681	c
	682	c
	683	c --------------------------------------------------------------------
	684	IMPLICIT NONE
	685
	686	#include "dimensions.h"
	687
	688	c
	689	c
	690	c
	691	c Arguments:
	692	c ----------
	693	real masse(llm),pente_max
	694	REAL q(llm),qm(llm+1)
	695	REAL w(llm+1)
	696	c
	697	c Local
	698	c ---------
	699	c
	700	INTEGER l
	701	c
	702	real dzq(llm),dzqw(llm),adzqw(llm),dzqmax
	703	real sigw, Mtot, MQtot
	704	integer m
	705	c integer ismax,ismin
	706
	707
	708	c On oriente tout dans le sens de la pression
	709	c W > 0 WHEN DOWN !!!!!!!!!!!!!
	710
	711	do l=2,llm
	712	dzqw(l)=q(l-1)-q(l)
	713	adzqw(l)=abs(dzqw(l))
	714	enddo
	715
	716	do l=2,llm-1
	717	if(dzqw(l)*dzqw(l+1).gt.0.) then
	718	dzq(l)=0.5*(dzqw(l)+dzqw(l+1))
	719	else
	720	dzq(l)=0.
	721	endif
	722	dzqmax=pente_max*min(adzqw(l),adzqw(l+1))
	723	dzq(l)=sign(min(abs(dzq(l)),dzqmax),dzq(l))
	724	enddo
	725
	726	dzq(1)=0.
	727	dzq(llm)=0.
	728
	729	do l = 1,llm-1
	730
	731	c Regular scheme (transfered mass < layer mass)
	732	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	733	if(w(l+1).gt.0. .and. w(l+1).le.masse(l+1)) then
	734	sigw=w(l+1)/masse(l+1)
	735	qm(l+1)=(q(l+1)+0.5(1.-sigw)dzq(l+1))
	736	else if(w(l+1).le.0. .and. -w(l+1).le.masse(l)) then
	737	sigw=w(l+1)/masse(l)
	738	qm(l+1)=(q(l)-0.5(1.+sigw)dzq(l))
	739
	740	c Extended scheme (transfered mass > layer mass)
	741	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	742	else if(w(l+1).gt.0.) then
	743	m=l+1
	744	Mtot = masse(m)
	745	MQtot = masse(m)*q(m)
	746	do while ((m.lt.llm).and.(w(l+1).gt.(Mtot+masse(m+1))))
	747	m=m+1
	748	Mtot = Mtot + masse(m)
	749	MQtot = MQtot + masse(m)*q(m)
	750	end do
	751	if (m.lt.llm) then
	752	sigw=(w(l+1)-Mtot)/masse(m+1)
	753	qm(l+1)= (1/w(l+1))(MQtot + (w(l+1)-Mtot)
	754	& (q(m+1)+0.5(1.-sigw)dzq(m+1)) )
	755	else
	756	w(l+1) = Mtot
	757	qm(l+1) = Mqtot / Mtot
	758	write(,) 'top layer is disapearing !'
	759	stop
	760	end if
	761	else ! if(w(l+1).lt.0)
	762	m = l-1
	763	Mtot = masse(m+1)
	764	MQtot = masse(m+1)*q(m+1)
[1867]	765	do while (m.gt.0)
	766	do while (-w(l+1).gt.(Mtot+masse(m)))
[57]	767	m=m-1
	768	Mtot = Mtot + masse(m+1)
	769	MQtot = MQtot + masse(m+1)*q(m+1)
[1867]	770	end do
[57]	771	end do
	772	if (m.gt.0) then
	773	sigw=(w(l+1)+Mtot)/masse(m)
	774	qm(l+1)= (-1/w(l+1))(MQtot + (-w(l+1)-Mtot)
	775	& (q(m)-0.5(1.+sigw)dzq(m)) )
	776	else
	777	qm(l+1)= (-1/w(l+1))(MQtot + (-w(l+1)-Mtot)qm(1))
	778	end if
	779	end if
	780	enddo
	781
	782	c boundary conditions (not used in newcondens !!)
	783	c qm(llm+1)=0.
	784	c if(w(1).gt.0.) then
	785	c qm(1)=q(1)
	786	c else
	787	c qm(1)=0.
	788	c end if
	789
	790	return
	791	end

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

Download in other formats:

Original Format