Context Navigation

newcondens.F @ 1112

Last change on this file since 1112 was 1047, checked in by emillour, 11 years ago

Mars GCM:

IMPORTANT CHANGE: Removed all reference/use of ngridmx (dimphys.h) in routines (necessary prerequisite to using parallel dynamics); in most cases this just means adding 'ngrid' as routine argument, and making local saved variables allocatable (and allocated at first call). In the process, had to convert many *.h files to equivalent modules: yomaer.h => yomaer_h.F90 , surfdat.h => surfdat_h.F90 , comsaison.h => comsaison_h.F90 , yomlw.h => yomlw_h.F90 , comdiurn.h => comdiurn_h.F90 , dimradmars.h => dimradmars_mod.F90 , comgeomfi.h => comgeomfi_h.F90, comsoil.h => comsoil_h.F90 , slope.h => slope_mod.F90
Also updated EOF routines, everything is now in eofdump_mod.F90
Removed unused routine lectfux.F (in dyn3d)

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

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

Download in other formats:

Original Format