Context Navigation

co2condens_mod.F @ 2156

Last change on this file since 2156 was 2155, checked in by aslmd, 5 years ago
another thing done in the CO2 physics that is not suitable for mesoscale, to be corrected
File size: 27.8 KB

Rev	Line
[2009]	1	MODULE co2condens_mod
	2
	3	IMPLICIT NONE
	4
	5	CONTAINS
	6
	7	SUBROUTINE co2condens(ngrid,nlayer,nq,ptimestep,
[38]	8	$ pcapcal,pplay,pplev,ptsrf,pt,
	9	$ pphi,pdt,pdu,pdv,pdtsrf,pu,pv,pq,pdq,
	10	$ piceco2,psolaralb,pemisurf,
	11	$ pdtc,pdtsrfc,pdpsrf,pduc,pdvc,pdqc,
[1996]	12	$ fluxsurf_sw,zls,
	13	$ zdqssed_co2,pcondicea_co2microp,
	14	$ zdtcloudco2)
[38]	15
[1036]	16	use tracer_mod, only: noms
[1047]	17	use surfdat_h, only: emissiv, phisfi
[2124]	18	use geometry_mod, only: latitude, ! grid point latitudes (rad)
	19	& longitude_deg, latitude_deg
[2009]	20	use planete_h, only: obliquit
	21	use comcstfi_h, only: cpp, g, r, pi
[1432]	22	#ifndef MESOSCALE
[2124]	23	USE vertical_layers_mod, ONLY: ap, bp
[1432]	24	#endif
[38]	25	IMPLICIT NONE
	26	c=======================================================================
	27	c subject:
	28	c --------
	29	c Condensation/sublimation of CO2 ice on the ground and in the
	30	c atmosphere
	31	c (Scheme described in Forget et al., Icarus, 1998)
	32	c
[2009]	33	c author: Francois Forget 1994-1996 ; updated 1996 -- 2018
[38]	34	c ------
[1996]	35	c adapted to external CO2 ice clouds scheme by Deborah Bardet (2018) '
[38]	36	c
	37	c=======================================================================
	38	c
	39	c 0. Declarations :
	40	c ------------------
	41	c
[1528]	42	include "callkeys.h"
[38]	43
	44	c-----------------------------------------------------------------------
	45	c Arguments :
	46	c ---------
[890]	47	INTEGER,INTENT(IN) :: ngrid ! number of atmospheric columns
	48	INTEGER,INTENT(IN) :: nlayer ! number of vertical layers
	49	INTEGER,INTENT(IN) :: nq ! number of tracers
[38]	50
[890]	51	REAL,INTENT(IN) :: ptimestep ! physics timestep (s)
	52	REAL,INTENT(IN) :: pcapcal(ngrid)
	53	REAL,INTENT(IN) :: pplay(ngrid,nlayer) !mid-layer pressure (Pa)
	54	REAL,INTENT(IN) :: pplev(ngrid,nlayer+1) ! inter-layer pressure (Pa)
	55	REAL,INTENT(IN) :: ptsrf(ngrid) ! surface temperature (K)
	56	REAL,INTENT(IN) :: pt(ngrid,nlayer) ! atmospheric temperature (K)
	57	REAL,INTENT(IN) :: pphi(ngrid,nlayer) ! geopotential (m2.s-2)
	58	REAL,INTENT(IN) :: pdt(ngrid,nlayer) ! tendency on temperature from
	59	! previous physical processes (K/s)
	60	REAL,INTENT(IN) :: pdu(ngrid,nlayer) ! tendency on zonal wind (m/s2)
	61	! from previous physical processes
	62	REAL,INTENT(IN) :: pdv(ngrid,nlayer) ! tendency on meridional wind (m/s2)
	63	! from previous physical processes
	64	REAL,INTENT(IN) :: pdtsrf(ngrid) ! tendency on surface temperature from
	65	! previous physical processes (K/s)
	66	REAL,INTENT(IN) :: pu(ngrid,nlayer) ! zonal wind (m/s)
	67	REAL,INTENT(IN) :: pv(ngrid,nlayer) ! meridional wind (m/s)
	68	REAL,INTENT(IN) :: pq(ngrid,nlayer,nq) ! tracers (../kg_air)
	69	REAL,INTENT(IN) :: pdq(ngrid,nlayer,nq) ! tendency on tracers from
	70	! previous physical processes
[1996]	71
	72	REAL,INTENT(IN) :: zdqssed_co2(ngrid) ! CO2 flux at the surface (kg.m-2.s-1)
	73	REAL,INTENT(IN) :: pcondicea_co2microp(ngrid,nlayer)! tendency due to CO2 condensation (kg/kg.s-1)
	74	REAL,INTENT(IN) :: zdtcloudco2(ngrid,nlayer) ! tendency on temperature due to latent heat
	75
[890]	76	REAL,INTENT(INOUT) :: piceco2(ngrid) ! CO2 ice on the surface (kg.m-2)
	77	REAL,INTENT(INOUT) :: psolaralb(ngrid,2) ! albedo of the surface
	78	REAL,INTENT(INOUT) :: pemisurf(ngrid) ! emissivity of the surface
	79
	80	! tendencies due to CO2 condensation/sublimation:
	81	REAL,INTENT(OUT) :: pdtc(ngrid,nlayer) ! tendency on temperature (K/s)
	82	REAL,INTENT(OUT) :: pdtsrfc(ngrid) ! tendency on surface temperature (K/s)
	83	REAL,INTENT(OUT) :: pdpsrf(ngrid) ! tendency on surface pressure (Pa/s)
	84	REAL,INTENT(OUT) :: pduc(ngrid,nlayer) ! tendency on zonal wind (m.s-2)
	85	REAL,INTENT(OUT) :: pdvc(ngrid,nlayer) ! tendency on meridional wind (m.s-2)
	86	REAL,INTENT(OUT) :: pdqc(ngrid,nlayer,nq) ! tendency on tracers
	87
	88	! added to calculate flux dependent albedo:
	89	REAL,intent(in) :: fluxsurf_sw(ngrid,2)
	90	real,intent(in) :: zls ! solar longitude (rad)
[38]	91
	92	c
	93	c Local variables :
	94	c -----------------
	95
	96	INTEGER i,j
[2009]	97	INTEGER l,ig,iq,icap
[890]	98	REAL zt(ngrid,nlayer)
[38]	99	REAL zcpi
[1114]	100	REAL ztcond (ngrid,nlayer+1) ! CO2 condensation temperature (atm)
	101	REAL ztcondsol(ngrid) ! CO2 condensation temperature (surface)
[890]	102	REAL zdiceco2(ngrid)
[2009]	103	REAL zcondicea(ngrid,nlayer) ! condensation rate in layer l (kg/m2/s)
	104	REAL zcondices(ngrid) ! condensation rate on the ground (kg/m2/s)
	105	REAL zfallice(ngrid,nlayer+1) ! amount of ice falling from layer l (kg/m2/s)
	106	REAL zfallheat
[890]	107	REAL zmflux(nlayer+1)
	108	REAL zu(nlayer),zv(nlayer)
	109	REAL zq(nlayer,nq),zq1(nlayer)
	110	REAL ztsrf(ngrid)
	111	REAL ztc(nlayer), ztm(nlayer+1)
	112	REAL zum(nlayer+1) , zvm(nlayer+1)
	113	REAL zqm(nlayer+1,nq),zqm1(nlayer+1)
	114	REAL masse(nlayer),w(nlayer+1)
	115	REAL Sm(nlayer),Smq(nlayer,nq),mixmas,qmix
[2124]	116	REAL availco2
[890]	117	LOGICAL condsub(ngrid)
[38]	118
[1224]	119	real :: emisref(ngrid)
	120
[38]	121	c variable speciale diagnostique
[890]	122	real tconda1(ngrid,nlayer)
	123	real tconda2(ngrid,nlayer)
	124	c REAL zdiceco2a(ngrid) ! for diagnostic only
	125	real zdtsig (ngrid,nlayer)
	126	real zdt (ngrid,nlayer)
	127	real vmr_co2(ngrid,nlayer) ! co2 volume mixing ratio
[38]	128	! improved_ztcond flag: If set to .true. (AND running with a 'co2' tracer)
	129	! then condensation temperature is computed using partial pressure of CO2
	130	logical,parameter :: improved_ztcond=.true.
	131
	132	c local saved variables
[890]	133	integer,save :: ico2 ! index of CO2 tracer
[2124]	134	real,save :: qco2,mmean
[890]	135	real,parameter :: latcond=5.9e5 ! (J/kg) Latent heat of solid CO2 ice
	136	real,parameter :: tcond1mb=136.27 ! condensation temperature (K) at 1 mbar
	137	real,parameter :: cpice=1000. ! (J.kg-1.K-1) specific heat of CO2 ice
	138	REAL,SAVE :: acond,bcond,ccond
	139	real,save :: m_co2, m_noco2, A , B
[38]	140
[890]	141	LOGICAL,SAVE :: firstcall = .true. !,firstcall2=.true.
[38]	142
[2124]	143	c D.BARDET: for debug
[1996]	144	real ztc3D(ngrid,nlayer)
	145	REAL ztm3D(ngrid,nlayer)
	146	REAL zmflux3D(ngrid,nlayer)
[38]	147	c----------------------------------------------------------------------
	148
	149	c Initialisation
	150	c --------------
	151	c
[1779]	152	! AS: firstcall OK absolute
[38]	153	IF (firstcall) THEN
[890]	154
[38]	155	bcond=1./tcond1mb
	156	ccond=cpp/(g*latcond)
	157	acond=r/latcond
	158
	159	firstcall=.false.
[2124]	160	write(,) 'CO2condens: improved_ztcond=',improved_ztcond
	161	PRINT*,'In co2condens:Tcond(P=1mb)=',tcond1mb,' Lcond=',latcond
[38]	162	PRINT*,'acond,bcond,ccond',acond,bcond,ccond
	163
	164	ico2=0
	165
	166	if (tracer) then
	167	c Prepare Special treatment if one of the tracer is CO2 gas
[1036]	168	do iq=1,nq
[38]	169	if (noms(iq).eq."co2") then
	170	ico2=iq
	171	m_co2 = 44.01E-3 ! CO2 molecular mass (kg/mol)
	172	m_noco2 = 33.37E-3 ! Non condensible mol mass (kg/mol)
	173	c Compute A and B coefficient use to compute
	174	c mean molecular mass Mair defined by
	175	c 1/Mair = q(ico2)/m_co2 + (1-q(ico2))/m_noco2
	176	c 1/Mair = A*q(ico2) + B
	177	A =(1/m_co2 - 1/m_noco2)
	178	B=1/m_noco2
	179	endif
	180	enddo
	181	end if
[890]	182	ENDIF ! of IF (firstcall)
[38]	183	zcpi=1./cpp
[1130]	184
[38]	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
[1047]	196	c pdtc(ngrid,nlayer) : dT/dt due to cond/sub
[38]	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
[1036]	210	DO iq=1,nq
[38]	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 Compute CO2 Volume mixing ratio
	233	c -------------------------------
	234	if (improved_ztcond.and.(ico2.ne.0)) then
	235	DO l=1,nlayer
	236	DO ig=1,ngrid
	237	qco2=pq(ig,l,ico2)+pdq(ig,l,ico2)*ptimestep
	238	c Mean air molecular mass = 1/(q(ico2)/m_co2 + (1-q(ico2))/m_noco2)
	239	mmean=1/(A*qco2 +B)
	240	vmr_co2(ig,l) = qco2*mmean/m_co2
	241	ENDDO
	242	ENDDO
	243	else
	244	DO l=1,nlayer
	245	DO ig=1,ngrid
	246	vmr_co2(ig,l)=0.95
	247	ENDDO
	248	ENDDO
	249	end if
	250
[1996]	251	IF (.NOT. co2clouds) then
[38]	252	c forecast of atmospheric temperature zt and frost temperature ztcond
	253	c --------------------------------------------------------------------
	254
	255	DO l=1,nlayer
	256	DO ig=1,ngrid
	257	zt(ig,l)=pt(ig,l)+ pdt(ig,l)*ptimestep
	258	! ztcond(ig,l)=1./(bcond-acondlog(.0095pplay(ig,l)))
[1263]	259	if (pplay(ig,l).ge.1e-4) then
	260	ztcond(ig,l)=
[38]	261	& 1./(bcond-acondlog(.01vmr_co2(ig,l)*pplay(ig,l)))
[1263]	262	else
	263	ztcond(ig,l)=0.0 !mars Monica
	264	endif
[38]	265	ENDDO
	266	ENDDO
[327]	267
[328]	268	ztcond(:,nlayer+1)=ztcond(:,nlayer)
[38]	269
	270	c Condensation/sublimation in the atmosphere
	271	c ------------------------------------------
	272	c (calcul of zcondicea , zfallice and pdtc)
	273	c
	274	DO l=nlayer , 1, -1
	275	DO ig=1,ngrid
	276	pdtc(ig,l)=0.
	277	IF((zt(ig,l).LT.ztcond(ig,l)).or.(zfallice(ig,l+1).gt.0))THEN
	278	condsub(ig)=.true.
	279	IF (zfallice(ig,l+1).gt.0) then
	280	zfallheat=zfallice(ig,l+1)*
	281	& (pphi(ig,l+1)-pphi(ig,l) +
	282	& cpice*(ztcond(ig,l+1)-ztcond(ig,l)))/latcond
	283	ELSE
	284	zfallheat=0.
	285	ENDIF
	286	pdtc(ig,l)=(ztcond(ig,l) - zt(ig,l))/ptimestep
	287	zcondicea(ig,l)=(pplev(ig,l)-pplev(ig,l+1))
	288	& ccondpdtc(ig,l)- zfallheat
	289	c Case when the ice from above sublimes entirely
	290	c """""""""""""""""""""""""""""""""""""""""""""""
	291	IF (zfallice(ig,l+1).lt.- zcondicea(ig,l)) then
	292	pdtc(ig,l)=(-zfallice(ig,l+1)+zfallheat)/
	293	& (ccond*(pplev(ig,l)-pplev(ig,l+1)))
	294	zcondicea(ig,l)= -zfallice(ig,l+1)
	295	END IF
	296
	297	zfallice(ig,l) = zcondicea(ig,l)+zfallice(ig,l+1)
	298	END IF
	299	ENDDO
	300	ENDDO
[1996]	301
	302	ELSE
	303
	304	DO ig=1,ngrid
	305	zfallice(ig,1) = zdqssed_co2(ig)
	306	ENDDO
	307	DO l=nlayer , 1, -1
	308	DO ig=1,ngrid
	309	zcondicea(ig,l) = pcondicea_co2microp(ig,l)*
	310	& (pplev(ig,l) - pplev(ig,l+1))/g
	311	ENDDO
	312	ENDDO
	313	DO l=nlayer, 1, -1
	314	DO ig=1, ngrid
	315	zt(ig,l)=pt(ig,l)+ pdt(ig,l)*ptimestep
	316	pdtc(ig,l)=0.
	317	ENDDO
	318	ENDDO
	319	ENDIF ! if not co2clouds
[38]	320
[1996]	321	call WRITEdiagfi(ngrid,"pdtc_atm",
	322	& "temperature tendency due to CO2 condensation",
	323	& " ",3,pdtc)
	324
	325	call WRITEdiagfi(ngrid,"zcondicea",
	326	& "",
	327	& " ",3,zcondicea)
	328
	329	call WRITEdiagfi(ngrid,"zfallice",
	330	& "",
	331	& " ",2,zfallice(ngrid,1))
	332
[38]	333	c *************************
	334	c SURFACE CONDENSATION
	335	c *************************
	336
	337	c forecast of ground temperature ztsrf and frost temperature ztcondsol
	338	c --------------------------------------------------------------------
	339	DO ig=1,ngrid
	340	ztcondsol(ig)=
	341	& 1./(bcond-acondlog(.01vmr_co2(ig,1)*pplev(ig,1)))
	342	ztsrf(ig) = ptsrf(ig) + pdtsrf(ig)*ptimestep
	343	ENDDO
	344
	345	c
	346	c Condensation/sublimation on the ground
	347	c --------------------------------------
[1114]	348	c (compute zcondices and pdtsrfc)
[38]	349	c
	350	DO ig=1,ngrid
[1541]	351	IF(latitude(ig).lt.0) THEN
[1114]	352	! Southern hemisphere
[38]	353	icap=2
	354	ELSE
[1114]	355	! Northern hemisphere
[38]	356	icap=1
	357	ENDIF
	358
	359	c
	360	c Loop on where we have condensation/ sublimation
	361	IF ((ztsrf(ig) .LT. ztcondsol(ig)) .OR. ! ground cond
	362	$ (zfallice(ig,1).NE.0.) .OR. ! falling snow
	363	$ ((ztsrf(ig) .GT. ztcondsol(ig)) .AND. ! ground sublim.
	364	$ ((piceco2(ig)+zfallice(ig,1)*ptimestep) .NE. 0.))) THEN
	365	condsub(ig) = .true.
	366
	367	IF (zfallice(ig,1).gt.0) then
	368	zfallheat=zfallice(ig,1)*
	369	& (pphi(ig,1)- phisfi(ig) +
[890]	370	& cpice*(ztcond(ig,1)-ztcondsol(ig)))/latcond
[38]	371	ELSE
	372	zfallheat=0.
	373	ENDIF
	374
	375	c condensation or partial sublimation of CO2 ice
	376	c """""""""""""""""""""""""""""""""""""""""""""""
	377	zcondices(ig)=pcapcal(ig)*(ztcondsol(ig)-ztsrf(ig))
	378	& /(latcond*ptimestep) - zfallheat
	379	pdtsrfc(ig) = (ztcondsol(ig) - ztsrf(ig))/ptimestep
[2124]	380	zdiceco2(ig) = zcondices(ig) + zfallice(ig,1)
	381
[2153]	382	#ifdef MESOSCALE
	383	print*, "not enough CO2 tracer in 1st layer to condense"
	384	print*, ">>> to be implemented in the mesoscale case"
	385	print*, "because this uses ap levels..."
	386	#else
[2124]	387	c If there is not enough CO2 tracer in 1st layer to condense
	388	c """"""""""""""""""""""""""""""""""""""""""""""""""""""
	389	IF(ico2.ne.0) then
	390	c Available CO2 tracer in layer 1 at end of timestep (kg/m2)
	391	availco2= pq(ig,1,ico2)*((ap(1)-ap(2))+
	392	& (bp(1)-bp(2))(pplev(ig,1)/g-zdiceco2(ig)ptimestep))
	393
	394	IF ((zcondices(ig) + zcondicea(ig,1))*ptimestep
	395	& .gt.availco2) then
	396	zcondices(ig) = availco2/ptimestep -zcondicea(ig,1)
	397	zdiceco2(ig) = zcondices(ig) + zfallice(ig,1)
	398	pdtsrfc(ig)=(latcond/pcapcal(ig))*
	399	& (zcondices(ig)+zfallheat)
	400	ENDIF
	401	ENDIF
[2153]	402	#endif
[2124]	403
[38]	404	c If the entire CO_2 ice layer sublimes
	405	c """"""""""""""""""""""""""""""""""""""""""""""""""""
	406	c (including what has just condensed in the atmosphere)
	407
	408	IF((piceco2(ig)/ptimestep+zfallice(ig,1)).LE.
	409	& -zcondices(ig))THEN
	410	zcondices(ig) = -piceco2(ig)/ptimestep - zfallice(ig,1)
	411	pdtsrfc(ig)=(latcond/pcapcal(ig))*
	412	& (zcondices(ig)+zfallheat)
[2124]	413	zdiceco2(ig) = zcondices(ig) + zfallice(ig,1)
[38]	414	END IF
	415
	416	c Changing CO2 ice amount and pressure :
	417	c """"""""""""""""""""""""""""""""""""
	418
	419	piceco2(ig) = piceco2(ig) + zdiceco2(ig)*ptimestep
	420	pdpsrf(ig) = -zdiceco2(ig)*g
	421
	422	IF(ABS(pdpsrf(ig)*ptimestep).GT.pplev(ig,1)) THEN
	423	PRINT*,'STOP in condens'
	424	PRINT*,'condensing more than total mass'
	425	PRINT*,'Grid point ',ig
[2124]	426	PRINT*,'Longitude(degrees): ',longitude_deg(ig)
	427	PRINT*,'Latitude(degrees): ',latitude_deg(ig)
[38]	428	PRINT*,'Ps = ',pplev(ig,1)
	429	PRINT*,'d Ps = ',pdpsrf(ig)
	430	STOP
	431	ENDIF
	432	END IF ! if there is condensation/sublimmation
	433	ENDDO ! of DO ig=1,ngrid
	434
	435	c ********************************************************************
	436	c Surface albedo and emissivity of the surface below the snow (emisref)
	437	c ********************************************************************
	438
	439	! Check that amont of CO2 ice is not problematic
	440	DO ig=1,ngrid
	441	if(.not.piceco2(ig).ge.0.) THEN
	442	if(piceco2(ig).le.-5.e-8) print*,
[2124]	443	$ 'WARNING co2condens piceco2(',ig,')=', piceco2(ig)
[38]	444	piceco2(ig)=0.
	445	endif
	446	ENDDO
	447
	448	! Set albedo and emissivity of the surface
	449	! ----------------------------------------
	450	CALL albedocaps(zls,ngrid,piceco2,psolaralb,emisref)
	451
	452	! set pemisurf() to emissiv when there is bare surface (needed for co2snow)
	453	DO ig=1,ngrid
	454	if (piceco2(ig).eq.0) then
	455	pemisurf(ig)=emissiv
	456	endif
	457	ENDDO
	458
	459	! firstcall2=.false.
	460	c ***************************************************************
	461	c Correction to account for redistribution between sigma or hybrid
	462	c layers when changing surface pressure (and warming/cooling
	463	c of the CO2 currently changing phase).
	464	c *************************************************************
	465
	466	DO ig=1,ngrid
	467	if (condsub(ig)) then
	468	do l=1,nlayer
	469	ztc(l) =zt(ig,l) +pdtc(ig,l) *ptimestep
	470	zu(l) =pu(ig,l) +pdu( ig,l) *ptimestep
	471	zv(l) =pv(ig,l) +pdv( ig,l) *ptimestep
[1036]	472	do iq=1,nq
[38]	473	zq(l,iq)=pq(ig,l,iq)+pdq(ig,l,iq)*ptimestep
	474	enddo
	475	end do
	476
	477	c Mass fluxes through the sigma levels (kg.m-2.s-1) (>0 when up)
	478	c """"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
	479
	480	zmflux(1) = -zcondices(ig)
	481	DO l=1,nlayer
	482	zmflux(l+1) = zmflux(l) -zcondicea(ig,l)
[1438]	483	#ifndef MESOSCALE
[38]	484	& + (bp(l)-bp(l+1))*(zfallice(ig,1)-zmflux(1))
	485	c zmflux set to 0 if very low to avoid: top layer is disappearing in v1ld
	486	if (abs(zmflux(l+1)).lt.1E-13.OR.bp(l+1).eq.0.) zmflux(l+1)=0.
[1438]	487	#else
	488	if (abs(zmflux(l+1)).lt.1E-13) zmflux(l+1)=0.
	489	#endif
[38]	490	END DO
	491
[2155]	492	#ifdef MESOSCALE
	493	print*,"absurd mass set because bp not available"
	494	print*,"TO BE FIXED"
	495	#else
[2124]	496	c Mass of each layer at the end of timestep
	497	c -----------------------------------------
[38]	498	DO l=1,nlayer
[2124]	499	masse(l)=( pplev(ig,l) - pplev(ig,l+1) +
	500	& (bp(l)-bp(l+1))pdpsrf(ig)ptimestep)/g
[38]	501	END DO
[2155]	502	#endif
[38]	503
	504	c Corresponding fluxes for T,U,V,Q
	505	c """"""""""""""""""""""""""""""""
	506
	507	c averaging operator for TRANSPORT
	508	c """"""""""""""""""""""""""""""""
	509	c Value transfert at the surface interface when condensation
	510	c sublimation:
	511	ztm(1) = ztsrf(ig) + pdtsrfc(ig)*ptimestep
	512	zum(1) = 0
	513	zvm(1) = 0
[1036]	514	do iq=1,nq
[38]	515	zqm(1,iq)=0. ! most tracer do not condense !
	516	enddo
	517	c Special case if one of the tracer is CO2 gas
	518	if (ico2.ne.0) zqm(1,ico2)=1. ! flux is 100% CO2
	519
	520	c Van Leer scheme:
	521	DO l=1,nlayer+1
	522	w(l)=-zmflux(l)*ptimestep
	523	END DO
[1269]	524	call vl1d(nlayer,ztc,2.,masse,w,ztm)
	525	call vl1d(nlayer,zu ,2.,masse,w,zum)
	526	call vl1d(nlayer,zv ,2.,masse,w,zvm)
[1036]	527	do iq=1,nq
[38]	528	do l=1,nlayer
	529	zq1(l)=zq(l,iq)
	530	enddo
	531	zqm1(1)=zqm(1,iq)
[1269]	532	call vl1d(nlayer,zq1,2.,masse,w,zqm1)
[38]	533	do l=2,nlayer
	534	zq( l,iq)=zq1(l)
	535	zqm(l,iq)=zqm1(l)
	536	enddo
	537	enddo
	538
	539	c Surface condensation affects low winds
	540	if (zmflux(1).lt.0) then
	541	zum(1)= zu(1) * (w(1)/masse(1))
	542	zvm(1)= zv(1) * (w(1)/masse(1))
	543	if (w(1).gt.masse(1)) then ! ensure numerical stability
	544	zum(1)= (zu(1)-zum(2))*masse(1)/w(1) + zum(2)
	545	zvm(1)= (zv(1)-zvm(2))*masse(1)/w(1) + zvm(2)
	546	end if
	547	end if
	548
	549	ztm(nlayer+1)= ztc(nlayer) ! should not be used, but...
	550	zum(nlayer+1)= zu(nlayer) ! should not be used, but...
	551	zvm(nlayer+1)= zv(nlayer) ! should not be used, but...
[1036]	552	do iq=1,nq
[38]	553	zqm(nlayer+1,iq)= zq(nlayer,iq)
	554	enddo
[86]	555
	556	#ifdef MESOSCALE
	557	!!!! AS: This part must be commented in the mesoscale model
	558	!!!! AS: ... to avoid instabilities.
	559	!!!! AS: you have to compile with -DMESOSCALE to do so
	560	#else
[38]	561	c Tendencies on T, U, V, Q
	562	c """"""""""""""""""""""""
	563	DO l=1,nlayer
[1996]	564	IF(.not. co2clouds) THEN
[38]	565	c Tendencies on T
	566	zdtsig(ig,l) = (1/masse(l)) *
	567	& ( zmflux(l)*(ztm(l) - ztc(l))
	568	& - zmflux(l+1)*(ztm(l+1) - ztc(l))
	569	& + zcondicea(ig,l)*(ztcond(ig,l)-ztc(l)) )
[1996]	570	ELSE
	571	zdtsig(ig,l) = (1/masse(l)) *
	572	& ( zmflux(l)*(ztm(l) - ztc(l))
	573	& - zmflux(l+1)*(ztm(l+1) - ztc(l)))
	574	ENDIF
	575	c D.BARDET: for diagnotics
	576	zmflux3D(ig,l)=zmflux(l)
	577	ztm3D(ig,l)=ztm(l)
	578	ztc3D(ig,l)=ztc(l)
	579
[38]	580	pdtc(ig,l) = pdtc(ig,l) + zdtsig(ig,l)
	581
	582	c Tendencies on U
	583	pduc(ig,l) = (1/masse(l)) *
	584	& ( zmflux(l)*(zum(l) - zu(l))
	585	& - zmflux(l+1)*(zum(l+1) - zu(l)) )
	586
	587
	588	c Tendencies on V
	589	pdvc(ig,l) = (1/masse(l)) *
	590	& ( zmflux(l)*(zvm(l) - zv(l))
	591	& - zmflux(l+1)*(zvm(l+1) - zv(l)) )
	592
	593	END DO
[1114]	594
[86]	595	#endif
[38]	596
	597	c Tendencies on Q
[1036]	598	do iq=1,nq
[38]	599	! if (noms(iq).eq.'co2') then
	600	if (iq.eq.ico2) then
	601	c SPECIAL Case when the tracer IS CO2 :
	602	DO l=1,nlayer
	603	pdqc(ig,l,iq)= (1/masse(l)) *
	604	& ( zmflux(l)*(zqm(l,iq) - zq(l,iq))
	605	& - zmflux(l+1)*(zqm(l+1,iq) - zq(l,iq))
	606	& + zcondicea(ig,l)*(zq(l,iq)-1.) )
	607	END DO
	608	else
	609	DO l=1,nlayer
	610	pdqc(ig,l,iq)= (1/masse(l)) *
	611	& ( zmflux(l)*(zqm(l,iq) - zq(l,iq))
	612	& - zmflux(l+1)*(zqm(l+1,iq) - zq(l,iq))
	613	& + zcondicea(ig,l)*zq(l,iq) )
	614	END DO
	615	end if
	616	enddo
	617
	618	end if ! if (condsub)
	619	END DO ! loop on ig
	620
	621	c ***************************************************************
	622	c CO2 snow / clouds scheme
	623	c ***************************************************************
	624
	625	call co2snow(ngrid,nlayer,ptimestep,emisref,condsub,pplev,
	626	& zcondicea,zcondices,zfallice,pemisurf)
	627
	628	c ***************************************************************
	629	c Ecriture des diagnostiques
	630	c ***************************************************************
	631
	632	c DO l=1,nlayer
	633	c DO ig=1,ngrid
	634	c Taux de cond en kg.m-2.pa-1.s-1
	635	c tconda1(ig,l)=zcondicea(ig,l)/(pplev(ig,l)-pplev(ig,l+1))
	636	c Taux de cond en kg.m-3.s-1
	637	c tconda2(ig,l)=tconda1(ig,l)pplay(ig,l)g/(r*pt(ig,l))
	638	c END DO
	639	c END DO
[1047]	640	c call WRITEDIAGFI(ngrid,'tconda1',
[38]	641	c &'Taux de condensation CO2 atmospherique /Pa',
	642	c & 'kg.m-2.Pa-1.s-1',3,tconda1)
[1047]	643	c call WRITEDIAGFI(ngrid,'tconda2',
[38]	644	c &'Taux de condensation CO2 atmospherique /m',
	645	c & 'kg.m-3.s-1',3,tconda2)
	646
	647	! output falling co2 ice in 1st layer:
[1047]	648	! call WRITEDIAGFI(ngrid,'fallice',
[38]	649	! &'Precipitation of co2 ice',
	650	! & 'kg.m-2.s-1',2,zfallice(1,1))
	651
[1114]	652	#ifndef MESOSCALE
	653	! Extra special case for surface temperature tendency pdtsrfc:
	654	! we want to fix the south pole temperature to CO2 condensation temperature
	655	if (caps.and.(obliquit.lt.27.)) then
	656	! check if last grid point is the south pole
[1541]	657	if (abs(latitude(ngrid)-(-pi/2.)).lt.1.e-5) then
[1114]	658	! NB: Updated surface pressure, at grid point 'ngrid', is
	659	! ps(ngrid)=pplev(ngrid,1)+pdpsrf(ngrid)*ptimestep
[2124]	660	! write(,) "co2condens: South pole: latitude(ngrid)=",
[1541]	661	! & latitude(ngrid)
[1114]	662	ztcondsol(ngrid)=
	663	& 1./(bcond-acondlog(.01vmr_co2(ngrid,1)*
	664	& (pplev(ngrid,1)+pdpsrf(ngrid)*ptimestep)))
	665	pdtsrfc(ngrid)=(ztcondsol(ngrid)-ztsrf(ngrid))/ptimestep
	666	endif
	667	endif
	668	#endif
	669
[2009]	670	END SUBROUTINE co2condens
[38]	671
	672
	673
	674	c *****************************************************************
[1269]	675	SUBROUTINE vl1d(nlayer,q,pente_max,masse,w,qm)
[38]	676	c
	677	c
	678	c Operateur de moyenne inter-couche pour calcul de transport type
	679	c Van-Leer " pseudo amont " dans la verticale
	680	c q,w sont des arguments d'entree pour le s-pg ....
	681	c masse : masse de la couche Dp/g
	682	c w : masse d'atm ``transferee'' a chaque pas de temps (kg.m-2)
	683	c pente_max = 2 conseillee
	684	c
	685	c
	686	c --------------------------------------------------------------------
	687	IMPLICIT NONE
	688
	689	c
	690	c
	691	c
	692	c Arguments:
	693	c ----------
[1270]	694	integer nlayer
[1269]	695	real masse(nlayer),pente_max
	696	REAL q(nlayer),qm(nlayer+1)
	697	REAL w(nlayer+1)
[38]	698	c
	699	c Local
	700	c ---------
	701	c
	702	INTEGER l
	703	c
[1269]	704	real dzq(nlayer),dzqw(nlayer),adzqw(nlayer),dzqmax
[38]	705	real sigw, Mtot, MQtot
	706	integer m
	707	c integer ismax,ismin
	708
	709
	710	c On oriente tout dans le sens de la pression
	711	c W > 0 WHEN DOWN !!!!!!!!!!!!!
	712
[1269]	713	do l=2,nlayer
[38]	714	dzqw(l)=q(l-1)-q(l)
	715	adzqw(l)=abs(dzqw(l))
	716	enddo
	717
[1269]	718	do l=2,nlayer-1
[38]	719	if(dzqw(l)*dzqw(l+1).gt.0.) then
	720	dzq(l)=0.5*(dzqw(l)+dzqw(l+1))
	721	else
	722	dzq(l)=0.
	723	endif
	724	dzqmax=pente_max*min(adzqw(l),adzqw(l+1))
	725	dzq(l)=sign(min(abs(dzq(l)),dzqmax),dzq(l))
	726	enddo
	727
	728	dzq(1)=0.
[1269]	729	dzq(nlayer)=0.
[38]	730
[1269]	731	do l = 1,nlayer-1
[38]	732
	733	c Regular scheme (transfered mass < layer mass)
	734	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	735	if(w(l+1).gt.0. .and. w(l+1).le.masse(l+1)) then
	736	sigw=w(l+1)/masse(l+1)
	737	qm(l+1)=(q(l+1)+0.5(1.-sigw)dzq(l+1))
	738	else if(w(l+1).le.0. .and. -w(l+1).le.masse(l)) then
	739	sigw=w(l+1)/masse(l)
	740	qm(l+1)=(q(l)-0.5(1.+sigw)dzq(l))
	741
	742	c Extended scheme (transfered mass > layer mass)
	743	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	744	else if(w(l+1).gt.0.) then
	745	m=l+1
	746	Mtot = masse(m)
	747	MQtot = masse(m)*q(m)
[1269]	748	do while ((m.lt.nlayer).and.(w(l+1).gt.(Mtot+masse(m+1))))
[38]	749	m=m+1
	750	Mtot = Mtot + masse(m)
	751	MQtot = MQtot + masse(m)*q(m)
	752	end do
[1269]	753	if (m.lt.nlayer) then
[38]	754	sigw=(w(l+1)-Mtot)/masse(m+1)
	755	qm(l+1)= (1/w(l+1))(MQtot + (w(l+1)-Mtot)
	756	& (q(m+1)+0.5(1.-sigw)dzq(m+1)) )
	757	else
	758	w(l+1) = Mtot
	759	qm(l+1) = Mqtot / Mtot
	760	write(,) 'top layer is disapearing !'
	761	stop
	762	end if
	763	else ! if(w(l+1).lt.0)
	764	m = l-1
	765	Mtot = masse(m+1)
	766	MQtot = masse(m+1)*q(m+1)
[120]	767	if (m.gt.0) then ! because some compilers will have problems
	768	! evaluating masse(0)
	769	do while ((m.gt.0).and.(-w(l+1).gt.(Mtot+masse(m))))
[38]	770	m=m-1
	771	Mtot = Mtot + masse(m+1)
	772	MQtot = MQtot + masse(m+1)*q(m+1)
[120]	773	if (m.eq.0) exit
	774	end do
	775	endif
[38]	776	if (m.gt.0) then
	777	sigw=(w(l+1)+Mtot)/masse(m)
	778	qm(l+1)= (-1/w(l+1))(MQtot + (-w(l+1)-Mtot)
	779	& (q(m)-0.5(1.+sigw)dzq(m)) )
	780	else
	781	qm(l+1)= (-1/w(l+1))(MQtot + (-w(l+1)-Mtot)qm(1))
	782	end if
	783	end if
	784	enddo
	785
[2124]	786	c boundary conditions (not used in co2condens !!)
[1269]	787	c qm(nlayer+1)=0.
[38]	788	c if(w(1).gt.0.) then
	789	c qm(1)=q(1)
	790	c else
	791	c qm(1)=0.
	792	c end if
	793
[2009]	794	END SUBROUTINE vl1d
	795
	796	END MODULE co2condens_mod

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

Download in other formats:

Original Format