Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

condens_n2.F90 @ 3572

Last change on this file since 3572 was 3175, checked in by emillour, 18 months ago
Pluto PCM: Add the old Pluto LMDZ for reference (required prior step to making an LMDZ.PLUTO using the same framework as the other physics packages). TB+EM
File size: 30.9 KB

Rev	Line
[3175]	1	subroutine condens_n2(ngrid,nlayer,nq,ptimestep, &
	2	pcapcal,pplay,pplev,ptsrf,pt, &
	3	pphi,pdt,pdu,pdv,pdtsrf,pu,pv,pq,pdq, &
	4	picen2,psolaralb,pemisurf, &
	5	pdtc,pdtsrfc,pdpsrf,pduc,pdvc, &
	6	pdqc,pdicen2)
	7
	8	use radinc_h, only : naerkind
	9	use comgeomfi_h
	10	implicit none
	11
	12	!==================================================================
	13	! Purpose
	14	! -------
	15	! Condense and/or sublime N2 ice on the ground and in the
	16	! atmosphere, and sediment the ice.
	17	!
	18	! Inputs
	19	! ------
	20	! ngrid Number of vertical columns
	21	! nlayer Number of layers
	22	! pplay(ngrid,nlayer) Pressure layers
	23	! pplev(ngrid,nlayer+1) Pressure levels
	24	! pt(ngrid,nlayer) Temperature (in K)
	25	! ptsrf(ngrid) Surface temperature
	26	!
	27	! pdt(ngrid,nlayermx) Time derivative before condensation/sublimation of pt
	28	! pdtsrf(ngrid) Time derivative before condensation/sublimation of ptsrf
	29	! picen2(ngrid) n2 ice at the surface (kg/m2)
	30	!
	31	! Outputs
	32	! -------
	33	! pdpsrf(ngrid) \ Contribution of condensation/sublimation
	34	! pdtc(ngrid,nlayermx) / to the time derivatives of Ps, pt, and ptsrf
	35	! pdtsrfc(ngrid) /
	36	! pdicen2(ngrid) Tendancy n2 ice at the surface (kg/m2)
	37	!
	38	! Both
	39	! ----
	40	! psolaralb(ngrid) Albedo at the surface
	41	! pemisurf(ngrid) Emissivity of the surface
	42	!
	43	! Authors
	44	! -------
	45	! Francois Forget (1996,2013)
	46	! Converted to Fortran 90 and slightly modified by R. Wordsworth (2009)
	47	!
	48	!
	49	!==================================================================
	50
	51	#include "dimensions.h"
	52	#include "dimphys.h"
	53	#include "comcstfi.h"
	54	#include "surfdat.h"
	55	#include "comvert.h"
	56	#include "callkeys.h"
	57	#include "tracer.h"
	58
	59	!-----------------------------------------------------------------------
	60	! Arguments
	61
	62	INTEGER ngrid, nlayer, nq
	63
	64	REAL ptimestep
	65	REAL pplay(ngrid,nlayer),pplev(ngrid,nlayer+1)
	66	REAL pcapcal(ngrid)
	67	REAL pt(ngrid,nlayer)
	68	REAL ptsrf(ngrid),flu1(ngrid),flu2(ngrid),flu3(ngrid)
	69	REAL pphi(ngrid,nlayer)
	70	REAL pdt(ngrid,nlayer),pdtsrf(ngrid),pdtc(ngrid,nlayer)
	71	REAL pdtsrfc(ngrid),pdpsrf(ngrid)
	72	REAL picen2(ngrid),psolaralb(ngrid),pemisurf(ngrid)
	73
	74
	75	REAL pu(ngrid,nlayer) , pv(ngrid,nlayer)
	76	REAL pdu(ngrid,nlayer) , pdv(ngrid,nlayer)
	77	REAL pduc(ngrid,nlayer) , pdvc(ngrid,nlayer)
	78	REAL pq(ngridmx,nlayer,nq),pdq(ngrid,nlayer,nq)
	79	REAL pdqc(ngrid,nlayer,nq)
	80
	81	!-----------------------------------------------------------------------
	82	! Local variables
	83
	84	INTEGER l,ig,ilay,it,iq,ich4_gas
	85
	86	REAL*8 zt(ngridmx,nlayermx)
	87	REAL tcond_n2
	88	REAL pcond_n2
	89	REAL glob_average2d ! function
	90	REAL zqn2(ngridmx,nlayermx) ! N2 MMR used to compute Tcond/zqn2
	91	REAL ztcond (ngridmx,nlayermx)
	92	REAL ztcondsol(ngridmx),zfallheat
	93	REAL pdicen2(ngridmx)
	94	REAL zcondicea(ngridmx,nlayermx), zcondices(ngridmx)
	95	REAL zfallice(ngridmx,nlayermx+1)
	96	REAL zmflux(nlayermx+1)
	97	REAL zu(nlayermx),zv(nlayermx)
	98	REAL zq(nlayermx,nqmx),zq1(nlayermx)
	99	REAL ztsrf(ngridmx)
	100	REAL ztc(nlayermx), ztm(nlayermx+1)
	101	REAL zum(nlayermx+1) , zvm(nlayermx+1)
	102	REAL zqm(nlayermx+1,nqmx),zqm1(nlayermx+1)
	103	LOGICAL condsub(ngridmx)
	104	REAL subptimestep
	105	Integer Ntime
	106	real masse (nlayermx),w(nlayermx+1)
	107	real wq(ngridmx,nlayermx+1)
	108	real vstokes,reff
	109	real dWtotsn2
	110	real condnconsn2(ngridmx)
	111	real nconsMAXn2
	112	! Special diagnostic variables
	113	real tconda1(ngridmx,nlayermx)
	114	real tconda2(ngridmx,nlayermx)
	115	real zdtsig (ngridmx,nlayermx)
	116	real zdtlatent (ngridmx,nlayermx)
	117	real zdt (ngridmx,nlayermx)
	118	REAL albediceF(ngridmx)
	119	SAVE albediceF
	120	INTEGER nsubtimestep,itsub !number of subtimestep when calling vl1d
	121	REAL subtimestep !ptimestep/nsubtimestep
	122	REAL dtmax
	123
	124	REAL zplevhist(ngridmx)
	125	REAL zplevnew(ngridmx)
	126	REAL zplev(ngridmx)
	127	REAL zpicen2(ngridmx)
	128	REAL ztsrfhist(ngridmx)
	129	REAL zdtsrf(ngridmx)
	130	REAL globzplevnew
	131
	132	REAL vmrn2(ngridmx)
	133	SAVE vmrn2
	134	REAL stephan
	135	DATA stephan/5.67e-08/ ! Stephan Boltzman constant
	136	SAVE stephan
	137	!-----------------------------------------------------------------------
	138	! Saved local variables
	139
	140	! REAL latcond
	141	REAL ccond
	142	REAL cpice ! for atm condensation
	143	SAVE cpice
	144	! SAVE latcond,ccond
	145	SAVE ccond
	146
	147	LOGICAL firstcall
	148	SAVE firstcall
	149	REAL SSUM
	150	EXTERNAL SSUM
	151
	152	! DATA latcond /2.5e5/
	153	! DATA latcond /1.98e5/
	154	DATA cpice /1300./
	155	DATA firstcall/.true./
	156
	157	INTEGER,SAVE :: i_n2ice=0 ! n2 ice
	158	CHARACTER(LEN=20) :: tracername ! to temporarily store text
	159	logical olkin ! option to prevent N2 ice effect in the south
	160	DATA olkin/.false./
	161	save olkin
	162
	163	CHARACTER(len=10) :: tname
	164
	165	!-----------------------------------------------------------------------
	166
	167	! Initialisation
	168	IF (firstcall) THEN
	169	ccond=cpp/(g*lw_n2)
	170	print*,'In condens_n2cloud: ccond=',ccond,' latcond=',lw_n2
	171
	172	! calculate global mean surface pressure for the fast mode
	173	IF (fast) THEN
	174	DO ig=1,ngridmx
	175	kp(ig) = exp(-phisfi(ig)/(r*38.))
	176	ENDDO
	177	p00=glob_average2d(kp) ! mean pres at ref level
	178	ENDIF
	179
	180	vmrn2(:) = 1.
	181	IF (ch4lag) then
	182	DO ig=1,ngridmx
	183	if (lati(ig)*180./pi.ge.latlag) then
	184	vmrn2(ig) = vmrlag
	185	endif
	186	ENDDO
	187	ENDIF
	188	IF (no_n2frost) then
	189	DO ig=1,ngridmx
	190	if (picen2(ig).eq.0.) then
	191	vmrn2(ig) = 1.e-15
	192	endif
	193	ENDDO
	194	ENDIF
	195	firstcall=.false.
	196	ENDIF
	197
	198	!-----------------------------------------------------------------------
	199	! Calculation of n2 condensation / sublimation
	200
	201	! Variables used:
	202	! picen2(ngrid) : amount of n2 ice on the ground (kg/m2)
	203	! zcondicea(ngrid,nlayermx): condensation rate in layer l (kg/m2/s)
	204	! zcondices(ngrid) : condensation rate on the ground (kg/m2/s)
	205	! zfallice(ngrid,nlayermx) : amount of ice falling from layer l (kg/m2/s)
	206	! zdtlatent(ngrid,nlayermx): dT/dt due to phase changes (K/s)
	207
	208	! Tendencies initially set to 0
	209	zcondices(1:ngrid) = 0.
	210	pdtsrfc(1:ngrid) = 0.
	211	pdpsrf(1:ngrid) = 0.
	212	ztsrfhist(1:ngrid) = 0.
	213	condsub(1:ngrid) = .false.
	214	pdicen2(1:ngrid) = 0.
	215	zfallheat=0
	216	pdqc(1:ngrid,1:nlayer,1:nq)=0
	217	pdtc(1:ngrid,1:nlayer)=0
	218	pduc(1:ngrid,1:nlayer)=0
	219	pdvc(1:ngrid,1:nlayer)=0
	220	zfallice(1:ngrid,1:nlayer+1)=0
	221	zcondicea(1:ngrid,1:nlayer)=0
	222	zdtlatent(1:ngrid,1:nlayer)=0
	223	zt(1:ngrid,1:nlayer)=0.
	224
	225	!-----------------------------------------------------------------------
	226	! Atmospheric condensation
	227
	228	! Condensation / sublimation in the atmosphere
	229	! --------------------------------------------
	230	! (calcul of zcondicea , zfallice and pdtc)
	231
	232	zt(1:ngrid,1:nlayer)=pt(1:ngrid,1:nlayer)+ pdt(1:ngrid,1:nlayer)*ptimestep
	233	if (igcm_n2.ne.0) then
	234	zqn2(1:ngrid,1:nlayer) = 1. ! & temporaire
	235	! zqn2(1:ngrid,1:nlayer)= pq(1:ngrid,1:nlayer,igcm_n2) + pdq(1:ngrid,1:nlayer,igcm_n2)*ptimestep
	236	else
	237	zqn2(1:ngrid,1:nlayer) = 1.
	238	end if
	239
	240	if (.not.fast) then
	241	! Forecast the atmospheric frost temperature 'ztcond' with function tcond_n2
	242	DO l=1,nlayer
	243	DO ig=1,ngrid
	244	ztcond (ig,l) = tcond_n2(pplay(ig,l),zqn2(ig,l))
	245	ENDDO
	246	ENDDO
	247
	248	DO l=nlayer,1,-1
	249	DO ig=1,ngrid
	250	pdtc(ig,l)=0. ! final tendancy temperature set to 0
	251
	252	IF((zt(ig,l).LT.ztcond(ig,l)).or.(zfallice(ig,l+1).gt.0))THEN
	253	condsub(ig)=.true. !condensation at level l
	254	IF (zfallice(ig,l+1).gt.0) then
	255	zfallheat=zfallice(ig,l+1)*&
	256	(pphi(ig,l+1)-pphi(ig,l) +&
	257	cpice*(ztcond(ig,l+1)-ztcond(ig,l)))/lw_n2
	258	ELSE
	259	zfallheat=0.
	260	ENDIF
	261	zdtlatent(ig,l)=(ztcond(ig,l) - zt(ig,l))/ptimestep
	262	zcondicea(ig,l)=(pplev(ig,l)-pplev(ig,l+1))&
	263	ccondzdtlatent(ig,l)- zfallheat
	264	! Case when the ice from above sublimes entirely
	265	! """""""""""""""""""""""""""""""""""""""""""""""
	266	IF ((zfallice(ig,l+1).lt.-zcondicea(ig,l)) &
	267	.AND. (zfallice(ig,l+1).gt.0)) THEN
	268
	269	zdtlatent(ig,l)=(-zfallice(ig,l+1)+zfallheat)/&
	270	(ccond*(pplev(ig,l)-pplev(ig,l+1)))
	271	zcondicea(ig,l)= -zfallice(ig,l+1)
	272	END IF
	273
	274	zfallice(ig,l) = zcondicea(ig,l)+zfallice(ig,l+1)
	275
	276	END IF
	277
	278	ENDDO
	279	ENDDO
	280	endif ! not fast
	281
	282	!-----------------------------------------------------------------------
	283	! Condensation/sublimation on the ground
	284	! (calculation of zcondices and pdtsrfc)
	285
	286	! Adding subtimesteps : in fast version, pressures too low lead to
	287	! instabilities.
	288	IF (fast) THEN
	289	IF (pplev(1,1).gt.0.3) THEN
	290	nsubtimestep= 1
	291	ELSE
	292	nsubtimestep= nbsub !max(nint(ptimestep/dtmax),1)
	293	ENDIF
	294	ELSE
	295	nsubtimestep= 1 ! if more then kp and p00 have to be calculated
	296	! for nofast mode
	297	ENDIF
	298	subtimestep=ptimestep/float(nsubtimestep)
	299
	300	do itsub=1,nsubtimestep
	301	! first loop : getting zplev, ztsurf
	302	IF (itsub.eq.1) then
	303	DO ig=1,ngrid
	304	zplev(ig)=pplev(ig,1)
	305	ztsrfhist(ig)=ptsrf(ig) + pdtsrf(ig)*ptimestep
	306	ztsrf(ig)=ptsrf(ig) + pdtsrf(ig)*subtimestep !!
	307	zpicen2(ig)=picen2(ig)
	308	ENDDO
	309	ELSE
	310	! additional loop :
	311	! 1) pressure updated
	312	! 2) direct redistribution of pressure over the globe
	313	! 3) modification pressure for unstable cases
	314	! 4) pressure update to conserve mass
	315	! 5) temperature updated with radiative tendancies
	316	DO ig=1,ngrid
	317	zplevhist(ig)=zplev(ig)
	318	zplevnew(ig)=zplev(ig)+pdpsrf(ig)*subtimestep ! 1)
	319	!IF (zplevnew(ig).le.0.0001) then
	320	! zplevnew(ig)=0.0001*kp(ig)/p00
	321	!ENDIF
	322	ENDDO
	323	! intermediaire de calcul: valeur moyenne de zplevnew (called twice in the code)
	324	globzplevnew=glob_average2d(zplevnew)
	325	DO ig=1,ngrid
	326	zplev(ig)=kp(ig)*globzplevnew/p00 ! 2)
	327	ENDDO
	328	DO ig=1,ngrid ! 3) unstable case condensation and sublimation: zplev=zplevhist
	329	IF (((pdpsrf(ig).lt.0.).and.(tcond_n2(zplev(ig),zqn2(ig,1)).le.ztsrf(ig))).or. &
	330	((pdpsrf(ig).gt.0.).and.(tcond_n2(zplev(ig),zqn2(ig,1)).ge.ztsrf(ig)))) then
	331	zplev(ig)=zplevhist(ig)
	332	ENDIF
	333	zplevhist(ig)=zplev(ig)
	334	ENDDO
	335	zplev=zplev*globzplevnew/glob_average2d(zplevhist) ! 4)
	336	DO ig=1,ngrid ! 5)
	337	zdtsrf(ig)=pdtsrf(ig) + (stephan/pcapcal(ig))(ptsrf(ig)4-ztsrf(ig)*4)
	338	ztsrf(ig)=ztsrf(ig)+pdtsrfc(ig)subtimestep+zdtsrf(ig)subtimestep
	339	zpicen2(ig)=zpicen2(ig)+pdicen2(ig)*subtimestep
	340	ENDDO
	341	ENDIF ! (itsub=1)
	342
	343	DO ig=1,ngrid
	344	! forecast of frost temperature ztcondsol
	345	!ztcondsol(ig) = tcond_n2(zplev(ig),zqn2(ig,1))
	346	ztcondsol(ig) = tcond_n2(zplev(ig),vmrn2(ig))
	347
	348	! Loop over where we have condensation / sublimation
	349	IF ((ztsrf(ig) .LT. ztcondsol(ig)) .OR. & ! ground cond
	350	((ztsrf(ig) .GT. ztcondsol(ig)) .AND. & ! ground sublim
	351	(zpicen2(ig) .GT. 0.))) THEN
	352	condsub(ig) = .true. ! condensation or sublimation
	353
	354	! Condensation or partial sublimation of N2 ice
	355	if (ztsrf(ig) .LT. ztcondsol(ig)) then ! condensation
	356	! Include a correction to account for the cooling of air near
	357	! the surface before condensing:
	358	zcondices(ig)=pcapcal(ig)*(ztcondsol(ig)-ztsrf(ig)) &
	359	/((lw_n2+cpp(zt(ig,1)-ztcondsol(ig)))subtimestep)
	360	else ! sublimation
	361	zcondices(ig)=pcapcal(ig)*(ztcondsol(ig)-ztsrf(ig)) &
	362	/(lw_n2*subtimestep)
	363	end if
	364
	365	pdtsrfc(ig) = (ztcondsol(ig) - ztsrf(ig))/subtimestep
	366
	367	! partial sublimation of N2 ice
	368	! If the entire N_2 ice layer sublimes
	369	! (including what has just condensed in the atmosphere)
	370	IF((zpicen2(ig)/subtimestep).LE. &
	371	-zcondices(ig))THEN
	372	zcondices(ig) = -zpicen2(ig)/subtimestep
	373	pdtsrfc(ig)=(lw_n2/pcapcal(ig))* &
	374	(zcondices(ig))
	375	END IF
	376
	377	! Changing N2 ice amount and pressure
	378
	379	pdicen2(ig) = zcondices(ig)
	380	pdpsrf(ig) = -pdicen2(ig)*g
	381	! pdpsrf(ig) = 0. ! OPTION to check impact N2 sub/cond
	382	IF (zplev(ig)+pdpsrf(ig)*subtimestep.le.0.0000001) then
	383	pdpsrf(ig)=(0.0000001*kp(ig)/p00-zplev(ig))/subtimestep
	384	pdicen2(ig)=-pdpsrf(ig)/g
	385	ENDIF
	386
	387	ELSE ! no condsub
	388	pdpsrf(ig)=0.
	389	pdicen2(ig)=0.
	390	pdtsrfc(ig)=0.
	391	ENDIF
	392	ENDDO ! ig
	393	enddo ! subtimestep
	394
	395	! Updating pressure, temperature and ice reservoir
	396	DO ig=1,ngrid
	397	pdpsrf(ig)=(zplev(ig)+pdpsrf(ig)*subtimestep-pplev(ig,1))/ptimestep
	398	! Two options here : 1 ok, 2 is wrong
	399	pdicen2(ig)=(zpicen2(ig)+pdicen2(ig)*subtimestep-picen2(ig))/ptimestep
	400	!pdicen2(ig)=-pdpsrf(ig)/g
	401
	402	pdtsrfc(ig)=((ztsrf(ig)+pdtsrfc(ig)*subtimestep)-(ztsrfhist(ig)))/ptimestep
	403
	404	! security
	405	if (picen2(ig) + pdicen2(ig)*ptimestep.lt.0.) then
	406	write(,) 'WARNING in condense_n2:'
	407	write(,) picen2(ig),pdicen2(ig)*ptimestep
	408	pdicen2(ig)= -picen2(ig)/ptimestep
	409	pdpsrf(ig)=-pdicen2(ig)*g
	410	endif
	411
	412	if(.not.picen2(ig).ge.0.) THEN
	413	! if(picen2(ig) + pdicen2(ig)*ptimestep.le.-1.e-8) then
	414	print, 'WARNING NEG RESERVOIR in condense_n2: picen2(',ig,')=', picen2(ig) + pdicen2(ig)ptimestep
	415	! pdicen2(ig)= -picen2(ig)/ptimestep
	416	! else
	417	picen2(ig)=0.0
	418	! endif
	419	endif
	420	ENDDO
	421
	422	! ***************************************************************
	423	! Correction to account for redistribution between sigma or hybrid
	424	! layers when changing surface pressure (and warming/cooling
	425	! of the n2 currently changing phase).
	426	! *************************************************************
	427	if (.not.fast) then
	428	DO ig=1,ngrid
	429	if (condsub(ig)) then
	430
	431	! Mass fluxes through the sigma levels (kg.m-2.s-1) (>0 when up)
	432	! """"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
	433	zmflux(1) = -zcondices(ig)
	434	DO l=1,nlayer
	435	zmflux(l+1) = zmflux(l) -zcondicea(ig,l) &
	436	+ (bp(l)-bp(l+1))*(zfallice(ig,1)-zmflux(1))
	437	! zmflux set to 0 if very low to avoid: top layer is disappearing in v1ld
	438	if (abs(zmflux(l+1)).lt.1E-13.OR.bp(l+1).eq.0.) zmflux(l+1)=0.
	439	END DO
	440
	441	! Mass of each layer
	442	! ------------------
	443	DO l=1,nlayer
	444	masse(l)=(pplev(ig,l) - pplev(ig,l+1))/g
	445	END DO
	446
	447
	448	! Corresponding fluxes for T,U,V,Q
	449	! """"""""""""""""""""""""""""""""
	450	! averaging operator for TRANSPORT
	451	! """"""""""""""""""""""""""""""""
	452
	453	! Subtimestep loop to perform the redistribution gently and simultaneously with
	454	! the other tendencies
	455	! Estimation of subtimestep (using only the first layer, the most critical)
	456	dtmax=ptimestep
	457	if (zmflux(1).gt.1.e-20) then
	458	dtmax=min(dtmax,masse(1)*zqn2(ig,1)/abs(zmflux(1)))
	459	endif
	460	nsubtimestep= max(nint(ptimestep/dtmax),nint(2.))
	461	subtimestep=ptimestep/float(nsubtimestep)
	462
	463	! New flux for each subtimestep
	464	do l=1,nlayer+1
	465	w(l)=-zmflux(l)*subtimestep
	466	enddo
	467	! initializing variables that will vary during subtimestep:
	468	do l=1,nlayer
	469	ztc(l) =pt(ig,l)
	470	zu(l) =pu(ig,l)
	471	zv(l) =pv(ig,l)
	472	do iq=1,nqmx
	473	zq(l,iq) = pq(ig,l,iq)
	474	enddo
	475	end do
	476
	477	! loop over nsubtimestep
	478	! """"""""""""""""""""""
	479	do itsub=1,nsubtimestep
	480	! Progressively adding tendancies from other processes.
	481	do l=1,nlayer
	482	ztc(l) =ztc(l) +(pdt(ig,l) + zdtlatent(ig,l))*subtimestep
	483	zu(l) =zu(l) +pdu( ig,l) * subtimestep
	484	zv(l) =zv(l) +pdv( ig,l) * subtimestep
	485	do iq=1,nqmx
	486	zq(l,iq) = zq(l,iq) + pdq(ig,l,iq)* subtimestep
	487	enddo
	488	end do
	489
	490	! Change of mass in each layer
	491	do l=1,nlayer
	492	masse(l)=masse(l)+pdpsrf(ig)subtimestep(pplev(ig,l) - pplev(ig,l+1))&
	493	/(g*pplev(ig,1))
	494	end do
	495
	496	ztm(2:nlayermx+1)=0.
	497	zum(2:nlayermx+1)=0.
	498	zvm(2:nlayermx+1)=0.
	499	zqm1(1:nlayermx+1)=0.
	500
	501	! Van Leer scheme:
	502	call vl1d(nlayer,ztc,2.,masse,w,ztm)
	503	call vl1d(nlayer,zu ,2.,masse,w,zum)
	504	call vl1d(nlayer,zv ,2.,masse,w,zvm)
	505	do iq=1,nqmx
	506	do l=1,nlayer
	507	zq1(l)=zq(l,iq)
	508	enddo
	509	zqm1(1)=zqm(1,iq)
	510	call vl1d(nlayer,zq1,2.,masse,w,zqm1)
	511	do l=2,nlayer
	512	zqm(l,iq)=zqm1(l)
	513	enddo
	514	enddo
	515
	516	! Correction to prevent negative value for qn2
	517	if (igcm_n2.ne.0) then
	518	zqm(1,igcm_n2)=1.
	519	do l=1,nlayer-1
	520	if (w(l)zqm(l,igcm_n2).gt.zq(l,igcm_n2)masse(l)) then
	521	zqm(l+1,igcm_n2)=max(zqm(l+1,igcm_n2), &
	522	(zqm(l,igcm_n2)w(l) -zq(l,igcm_n2)masse(l))/w(l+1) )
	523	else
	524	exit
	525	endif
	526	end do
	527	end if
	528
	529	! Value transfert at the surface interface when condensation sublimation:
	530
	531	if (zmflux(1).lt.0) then
	532	! Surface condensation
	533	zum(1)= zu(1)
	534	zvm(1)= zv(1)
	535	ztm(1) = ztc(1)
	536	else
	537	! Surface sublimation:
	538	ztm(1) = ztsrf(ig) + pdtsrfc(ig)*ptimestep
	539	zum(1) = 0
	540	zvm(1) = 0
	541	end if
	542	do iq=1,nqmx
	543	zqm(1,iq)=0. ! most tracer do not condense !
	544	enddo
	545	! Special case if the tracer is n2 gas
	546	if (igcm_n2.ne.0) zqm(1,igcm_n2)=1.
	547
	548	!!! Source haze: 0.02 pourcent when n2 sublimes
	549	IF (source_haze) THEN
	550	IF (pdicen2(ig).lt.0) THEN
	551	DO iq=1,nq
	552	tname=noms(iq)
	553	if (tname(1:4).eq."haze") then
	554	!zqm(1,iq)=0.02
	555	!zqm(1,iq)=-pdicen2(ig)*0.02
	556	zqm(1,iq)=-pdicen2(ig)ptimestep0.02
	557	!zqm(10,iq)=-pdicen2(ig)ptimestep100.
	558	!zqm(1,iq)=-pdicen2(ig)*1000000.
	559
	560	endif
	561	ENDDO
	562	ENDIF
	563	ENDIF
	564	ztm(nlayer+1)= ztc(nlayer) ! should not be used, but...
	565	zum(nlayer+1)= zu(nlayer) ! should not be used, but...
	566	zvm(nlayer+1)= zv(nlayer) ! should not be used, but...
	567	do iq=1,nqmx
	568	zqm(nlayer+1,iq)= zq(nlayer,iq)
	569	enddo
	570
	571	! Tendencies on T, U, V, Q
	572	! """""""""""""""""""""""
	573	DO l=1,nlayer
	574
	575	! Tendencies on T
	576	zdtsig(ig,l) = (1/masse(l)) * &
	577	( zmflux(l)*(ztm(l) - ztc(l)) &
	578	- zmflux(l+1)*(ztm(l+1) - ztc(l)) &
	579	+ zcondicea(ig,l)*(ztcond(ig,l)-ztc(l)) )
	580
	581	! Tendencies on U
	582	pduc(ig,l) = (1/masse(l)) * &
	583	( zmflux(l)*(zum(l) - zu(l))&
	584	- zmflux(l+1)*(zum(l+1) - zu(l)) )
	585
	586	! Tendencies on V
	587	pdvc(ig,l) = (1/masse(l)) * &
	588	( zmflux(l)*(zvm(l) - zv(l)) &
	589	- zmflux(l+1)*(zvm(l+1) - zv(l)) )
	590
	591	END DO
	592
	593	! Tendencies on Q
	594	do iq=1,nqmx
	595	if (iq.eq.igcm_n2) then
	596	! SPECIAL Case when the tracer IS N2 :
	597	DO l=1,nlayer
	598	pdqc(ig,l,iq)= (1/masse(l)) * &
	599	( zmflux(l)*(zqm(l,iq) - zq(l,iq)) &
	600	- zmflux(l+1)*(zqm(l+1,iq) - zq(l,iq))&
	601	+ zcondicea(ig,l)*(zq(l,iq)-1.) )
	602	END DO
	603	else
	604	DO l=1,nlayer
	605	pdqc(ig,l,iq)= (1/masse(l)) * &
	606	( zmflux(l)*(zqm(l,iq) - zq(l,iq)) &
	607	- zmflux(l+1)*(zqm(l+1,iq) - zq(l,iq)) &
	608	+ zcondicea(ig,l)*zq(l,iq) )
	609	END DO
	610	end if
	611	enddo
	612	! Update variables at the end of each subtimestep.
	613	do l=1,nlayer
	614	ztc(l) =ztc(l) + zdtsig(ig,l) *subtimestep
	615	zu(l) =zu(l) + pduc(ig,l) *subtimestep
	616	zv(l) =zv(l) + pdvc(ig,l) *subtimestep
	617	do iq=1,nqmx
	618	zq(l,iq) = zq(l,iq) + pdqc(ig,l,iq)*subtimestep
	619	enddo
	620	end do
	621	enddo ! loop on nsubtimestep
	622	! Recomputing Total tendencies
	623	do l=1,nlayer
	624	pdtc(ig,l) = (ztc(l) - zt(ig,l) )/ptimestep
	625	pduc(ig,l) = (zu(l) - (pu(ig,l) + pdu(ig,l)*ptimestep))/ptimestep
	626	pdvc(ig,l) = (zv(l) - (pv(ig,l) + pdv(ig,l)*ptimestep))/ptimestep
	627	do iq=1,nqmx
	628	pdqc(ig,l,iq) = (zq(l,iq) - (pq(ig,l,iq) + pdq(ig,l,iq)*ptimestep))/ptimestep
	629
	630
	631	! Correction temporaire
	632	if (iq.eq.igcm_n2) then
	633	if((pq(ig,l,iq) +(pdqc(ig,l,iq)+ pdq(ig,l,iq))*ptimestep) &
	634	.lt.0.01) then ! if n2 < 1 % !
	635	pdqc(ig,l,iq)=(0.01-pq(ig,l,iq))/ptimestep-pdq(ig,l,iq)
	636	end if
	637	end if
	638
	639	enddo
	640	end do
	641	! *****************************TEMPORAIRE ****************
	642	if (ngrid.eq.1) then
	643	write(,) 'nsubtimestep=' ,nsubtimestep
	644	write(,) 'masse avant' , (pplev(ig,1) - pplev(ig,2))/g
	645	write(,) 'masse apres' , masse(1)
	646	write(,) 'zmfluxDT, l=1' , zmflux(1)ptimestep
	647	write(,) 'zmfluxDT, l=2' , zmflux(2)ptimestep
	648	write(,) 'pq, l=1,2,3' , pq(1,1,1), pq(1,2,1),pq(1,3,1)
	649	write(,) 'zq, l=1,2,3' , zq(1,1), zq(2,1),zq(3,1)
	650	write(,) 'dqDt l=1' , pdq(1,1,1)ptimestep
	651	write(,) 'dqcondDt l=1' , pdqc(1,1,1)ptimestep
	652	end if
	653
	654	! ***********************************************************
	655	end if ! if (condsub)
	656	END DO ! loop on ig
	657	endif ! not fast
	658
	659	! ********** Option Olkin to prevent N2 effect in the south******
	660	112 continue
	661	if (olkin) then
	662	DO ig=1,ngrid
	663	if (lati(ig).lt.0) then
	664	pdtsrfc(ig) = max(0.,pdtsrfc(ig))
	665	pdpsrf(ig) = 0.
	666	pdicen2(ig) = 0.
	667	do l=1,nlayer
	668	pdtc(ig,l) = max(0.,zdtlatent(ig,l))
	669	pduc(ig,l) = 0.
	670	pdvc(ig,l) = 0.
	671	do iq=1,nqmx
	672	pdqc(ig,l,iq) = 0
	673	enddo
	674	end do
	675	end if
	676	END DO
	677	end if
	678	! *******************************************************************
	679
	680	! ***************************************************************
	681	! Ecriture des diagnostiques
	682	! ***************************************************************
	683
	684	! DO l=1,nlayer
	685	! DO ig=1,ngrid
	686	! Taux de cond en kg.m-2.pa-1.s-1
	687	! tconda1(ig,l)=zcondicea(ig,l)/(pplev(ig,l)-pplev(ig,l+1))
	688	! Taux de cond en kg.m-3.s-1
	689	! tconda2(ig,l)=tconda1(ig,l)pplay(ig,l)g/(r*pt(ig,l))
	690	! END DO
	691	! END DO
	692	! call WRITEDIAGFI(ngridmx,'tconda1', &
	693	! 'Taux de condensation N2 atmospherique /Pa', &
	694	! 'kg.m-2.Pa-1.s-1',3,tconda1)
	695	! call WRITEDIAGFI(ngridmx,'tconda2', &
	696	! 'Taux de condensation N2 atmospherique /m', &
	697	! 'kg.m-3.s-1',3,tconda2)
	698
	699
	700	return
	701	end subroutine condens_n2
	702
	703	!-------------------------------------------------------------------------
	704
	705	real function tcond_n2(p,vmr)
	706	! Calculates the condensation temperature for N2 at pressure P and vmr
	707	implicit none
	708	real, intent(in):: p,vmr
	709
	710	! tcond_n2 = (1.)/(0.026315-0.0011877log(.7143p*vmr))
	711	IF (p.ge.0.529995) then
	712	! tcond Fray and Schmitt for N2 phase beta (T>35.6 K) FIT TB
	713	! tcond_n2 = (1.)/(1./63.1470-296.925/(2.5e50.98)log(1./(0.1255701.e5)p*vmr))
	714	tcond_n2 = (1.)/(0.01583606505-1.211938776e-3log(7.963685594e-5p*vmr))
	715	ELSE
	716	! tcond Fray and Schmitt for N2 phase alpha(T<35.6 K) FIT BT
	717	! tcond_n2 = (1.)/(1./35.6-296.925/(2.5e51.09)log(1./(0.508059)pvmr))
	718	tcond_n2 = (1.)/(1./35.6-1.089633028e-3log(1.968275338p*vmr))
	719	ENDIF
	720	return
	721	end function tcond_n2
	722
	723	!-------------------------------------------------------------------------
	724
	725	real function pcond_n2(t,vmr)
	726	! Calculates the condensation pressure for N2 at temperature T and vmr
	727	implicit none
	728	real, intent(in):: t,vmr
	729
	730	! tcond_n2 = (1.)/(0.026315-0.0011877log(.7143p*vmr))
	731	IF (t.ge.35.6) then
	732	! tcond Fray and Schmitt for N2 phase beta (T>35.6 K) FIT TB
	733	! pcond_n2 = 0.1255701.e5/vmrexp((2.5e50.98)/296.925(1./63.1470-1./t))
	734	pcond_n2 = 0.125570e5/vmrexp(825.1241896(1./63.147-1./t))
	735	ELSE
	736	! tcond Fray and Schmitt for N2 phase alpha(T<35.6 K) FIT TB
	737	! pcond_n2 = 0.508059/vmrexp((2.5e51.09)/296.925*(1./35.6-1./t))
	738	pcond_n2 = 0.508059/vmrexp(917.7401701(1./35.6-1./t))
	739	ENDIF
	740	return
	741	end function pcond_n2
	742
	743	!-------------------------------------------------------------------------
	744
	745	real function glob_average2d(var)
	746	! Calculates the global average of variable var
	747	use comgeomfi_h
	748	implicit none
	749	#include "dimensions.h"
	750	#include "dimphys.h"
	751
	752	! INTEGER ngrid
	753	REAL var(ngridmx)
	754	INTEGER ig
	755
	756	glob_average2d = 0.
	757	DO ig=2,ngridmx-1
	758	glob_average2d = glob_average2d + var(ig)*area(ig)
	759	END DO
	760	glob_average2d = glob_average2d + var(1)area(1)iim
	761	glob_average2d = glob_average2d + var(ngridmx)area(ngridmx)iim
	762	glob_average2d = glob_average2d/(totarea+(area(1)+area(ngridmx))*(iim-1))
	763
	764	end function glob_average2d
	765
	766	! *****************************************************************
	767
	768	subroutine vl1d(nlayer,q,pente_max,masse,w,qm)
	769	!
	770	! Operateur de moyenne inter-couche pour calcul de transport type
	771	! Van-Leer " pseudo amont " dans la verticale
	772	! q,w sont des arguments d'entree pour le s-pg ....
	773	! masse : masse de la couche Dp/g
	774	! w : masse d'atm ``transferee'' a chaque pas de temps (kg.m-2)
	775	! pente_max = 2 conseillee
	776	! --------------------------------------------------------------------
	777	IMPLICIT NONE
	778
	779	#include "dimensions.h"
	780
	781	! Arguments:
	782	! ----------
	783	integer nlayer
	784	real masse(nlayer),pente_max
	785	REAL q(nlayer),qm(nlayer+1)
	786	REAL w(nlayer+1)
	787	!
	788	! Local
	789	! ---------
	790	!
	791	INTEGER l
	792	!
	793	real dzq(nlayer),dzqw(nlayer),adzqw(nlayer),dzqmax
	794	real sigw, Mtot, MQtot
	795	integer m
	796
	797
	798	! On oriente tout dans le sens de la pression
	799	! W > 0 WHEN DOWN !!!!!!!!!!!!!
	800
	801	do l=2,nlayer
	802	dzqw(l)=q(l-1)-q(l)
	803	adzqw(l)=abs(dzqw(l))
	804	enddo
	805
	806	do l=2,nlayer-1
	807	if(dzqw(l)*dzqw(l+1).gt.0.) then
	808	dzq(l)=0.5*(dzqw(l)+dzqw(l+1))
	809	else
	810	dzq(l)=0.
	811	endif
	812	dzqmax=pente_max*min(adzqw(l),adzqw(l+1))
	813	dzq(l)=sign(min(abs(dzq(l)),dzqmax),dzq(l))
	814	enddo
	815
	816	dzq(1)=0.
	817	dzq(nlayer)=0.
	818
	819	do l = 1,nlayer-1
	820
	821	! Regular scheme (transfered mass < layer mass)
	822	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	823	if(w(l+1).gt.0. .and. w(l+1).le.masse(l+1)) then
	824	sigw=w(l+1)/masse(l+1)
	825	qm(l+1)=(q(l+1)+0.5(1.-sigw)dzq(l+1))
	826	else if(w(l+1).le.0. .and. -w(l+1).le.masse(l)) then
	827	sigw=w(l+1)/masse(l)
	828	qm(l+1)=(q(l)-0.5(1.+sigw)dzq(l))
	829
	830	! Extended scheme (transfered mass > layer mass)
	831	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	832	else if(w(l+1).gt.0.) then
	833	m=l+1
	834	Mtot = masse(m)
	835	MQtot = masse(m)*q(m)
	836	!if (m.lt.nlayer) then ! because some compilers will have problems
	837	! ! evaluating masse(nlayer+1)
	838	do while ((m.lt.nlayer).and.(w(l+1).gt.(Mtot+masse(m+1))))
	839	m=m+1
	840	Mtot = Mtot + masse(m)
	841	MQtot = MQtot + masse(m)*q(m)
	842	! if (m.eq.nlayer) exit
	843	end do
	844	!endif
	845	if (m.lt.nlayer) then
	846	sigw=(w(l+1)-Mtot)/masse(m+1)
	847	qm(l+1)= (1/w(l+1))(MQtot + (w(l+1)-Mtot) &
	848	(q(m+1)+0.5(1.-sigw)dzq(m+1)) )
	849	else
	850	w(l+1) = Mtot
	851	qm(l+1) = Mqtot / Mtot
	852	write(,) 'top layer is disapearing !'
	853	stop
	854	end if
	855	else ! if(w(l+1).lt.0)
	856	m = l-1
	857	Mtot = masse(m+1)
	858	MQtot = masse(m+1)*q(m+1)
	859	if (m.gt.0) then ! because some compilers will have problems
	860	! evaluating masse(0)
	861	do while ((m.gt.0).and.(-w(l+1).gt.(Mtot+masse(m))))
	862	m=m-1
	863	Mtot = Mtot + masse(m+1)
	864	MQtot = MQtot + masse(m+1)*q(m+1)
	865	if (m.eq.0) exit
	866	end do
	867	endif
	868	if (m.gt.0) then
	869	sigw=(w(l+1)+Mtot)/masse(m)
	870	qm(l+1)= (-1/w(l+1))(MQtot + (-w(l+1)-Mtot) &
	871	(q(m)-0.5(1.+sigw)dzq(m)) )
	872	else
	873	qm(l+1)= (-1/w(l+1))(MQtot + (-w(l+1)-Mtot)qm(1))
	874	end if
	875	end if
	876	enddo
	877
	878	return
	879	end subroutine vl1d
	880

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

Download in other formats:

Original Format