Context Navigation

newcondens.F @ 1242

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

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

Download in other formats:

Original Format