Context Navigation

vdifc_pluto_mod.F90 @ 3501

Last change on this file since 3501 was 3258, checked in by afalco, 16 months ago
Pluto PCM: Methane/CO taken into account in callcorrk AF
File size: 24.4 KB

Rev	Line
[3247]	1	module vdifc_pluto_mod
	2
	3	implicit none
	4
	5	contains
	6
	7	SUBROUTINE vdifc_pluto(ngrid,nlay,nq,ppopsk, &
	8	ptimestep,pcapcal,lecrit, &
	9	pplay,pplev,pzlay,pzlev,pz0, &
	10	pu,pv,ph,pq,pt,ptsrf,pemis,pqsurf, &
	11	pdufi,pdvfi,pdhfi,pdqfi,pdtfi,pfluxsrf, &
	12	pdudif,pdvdif,pdhdif,pdtsrf,pq2, &
	13	pdqdif,pdqsdif,qsat_ch4,qsat_ch4_l1)
	14
	15	use comgeomfi_h
	16	use callkeys_mod, only: carbox, methane, condcosurf, condensn2, condmetsurf,&
	17	kmix_proffix, vertdiff, tracer, kmixmin
	18	use datafile_mod, only: datadir
	19	use surfdat_h, only: phisfi
[3258]	20	use comcstfi_mod, only: g, r, rcp, cpp
[3247]	21	USE tracer_h, only: igcm_ch4_gas, igcm_ch4_ice, igcm_co_gas, igcm_co_ice,&
	22	igcm_n2, lw_ch4, lw_co, lw_n2
	23
	24	implicit none
	25
	26	!=======================================================================
	27	!
	28	! subject:
	29	! --------
	30	! Turbulent diffusion (mixing) for potential T, U, V and tracer
	31	!
	32	! Shema implicite
	33	! On commence par rajouter au variables x la tendance physique
	34	! et on resoult en fait:
	35	! x(t+1) = x(t) + dt * (dx/dt)phys(t) + dt * (dx/dt)difv(t+1)
	36	!
	37	! author:
	38	! ------
	39	! Hourdin/Forget/Fournier
	40	!=======================================================================
	41
	42	!-----------------------------------------------------------------------
	43	! declarations:
	44	! -------------
	45
	46	#include "dimensions.h"
	47
	48	!
	49	! arguments:
	50	! ----------
	51
	52	INTEGER ngrid,nlay
	53	REAL ptimestep
	54	REAL pplay(ngrid,nlay),pplev(ngrid,nlay+1)
	55	REAL pzlay(ngrid,nlay),pzlev(ngrid,nlay+1)
	56	REAL pu(ngrid,nlay),pv(ngrid,nlay),ph(ngrid,nlay)
	57	REAL ptsrf(ngrid),pemis(ngrid)
	58	REAL pdufi(ngrid,nlay),pdvfi(ngrid,nlay),pdhfi(ngrid,nlay)
	59	REAL pdtfi(ngrid,nlay)
	60	REAL pt(ngrid,nlay)
	61	REAL pfluxsrf(ngrid)
	62	REAL pdudif(ngrid,nlay),pdvdif(ngrid,nlay),pdhdif(ngrid,nlay)
	63	REAL pdtsrf(ngrid),pcapcal(ngrid)
	64	REAL pq2(ngrid,nlay+1)
	65	REAL qsat_ch4(ngrid)
	66	REAL qsat_co(ngrid)
	67	REAL qsat_ch4_l1(ngrid)
	68	REAL zq1temp_ch4(ngrid)
	69
	70	! Argument added for condensation:
	71	! REAL n2ice (ngrid)
	72	REAL ppopsk(ngrid,nlay)
	73	logical lecrit
	74	REAL pz0
	75
	76	! Traceurs :
	77	integer nq
	78	REAL pqsurf(ngrid,nq)
	79	real pq(ngrid,nlay,nq), pdqfi(ngrid,nlay,nq)
	80	real pdqdif(ngrid,nlay,nq)
	81	real pdqdifeddy(ngrid,nlay,nq)
	82	real pdqsdif(ngrid,nq),pdqsdif1(ngrid,nq)
	83
	84	! local:
	85	! ------
	86
	87	INTEGER ilev,ig,ilay,nlev
	88
	89	REAL z4st,zdplanck(ngrid)
	90	REAL zkv(ngrid,nlay+1),zkh(ngrid,nlay+1)
	91	REAL zcdv(ngrid),zcdh(ngrid),sat2(ngrid)
	92	REAL zcdv_true(ngrid),zcdh_true(ngrid)
	93	REAL zu(ngrid,nlay),zv(ngrid,nlay)
	94	REAL zh(ngrid,nlay),zt(ngrid,nlay)
	95	REAL ztsrf2(ngrid),sat(ngrid),sat1(ngrid)
	96	REAL z1(ngrid),z2(ngrid)
	97	REAL za(ngrid,nlay),zb(ngrid,nlay)
	98	REAL zb0(ngrid,nlay)
	99	REAL zc(ngrid,nlay),zd(ngrid,nlay)
	100	REAL zcst1
	101	REAL zu2
	102
	103	EXTERNAL SSUM,SCOPY
	104	REAL SSUM
	105	LOGICAL firstcall
	106	SAVE firstcall
	107
	108	!!read fixed profile for kmix
	109	integer Nfine,ifine
	110	parameter(Nfine=701)
	111	character(len=100) :: file_path
	112	real,save :: levdat(Nfine),kmixdat(Nfine)
	113	real :: kmix_z(nlay) ! kmix from kmix_proffix
	114
	115	! variable added for N2 condensation:
	116	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	117	REAL hh , zhcond(ngrid,nlay)
	118	! REAL latcond,tcond1p4Pa
	119	REAL tcond1p4Pa
	120	REAL acond,bcond
	121	SAVE acond,bcond
	122	! DATA latcond,tcond1p4Pa/2.5e5,38/
	123	DATA tcond1p4Pa/38/
	124
	125	! Tracers :
	126	! ~~~~~~~
	127	INTEGER iq
	128	REAL zq(ngrid,nlay,nq)
	129	REAL zq1temp_co(ngrid)
	130	REAL rho(ngrid) ! near surface air density
	131	DATA firstcall/.true./
	132
	133	! ** un petit test de coherence
	134	! --------------------------
	135
	136	IF (firstcall) THEN
	137	IF(ngrid.NE.ngrid) THEN
	138	write(,) 'STOP dans vdifc'
	139	write(,) 'probleme de dimensions :'
	140	write(,) 'ngrid =',ngrid
	141	write(,) 'ngrid =',ngrid
	142	STOP
	143	ENDIF
	144	! To compute: Tcond= 1./(bcond-acondlog(.7143p)) (p in pascal)
	145	write(,) 'In vdifc: Tcond(P=1.4Pa)=',tcond1p4Pa,' Lcond=',lw_n2
	146	bcond=1./tcond1p4Pa
	147	acond=r/lw_n2
	148	write(,) ' acond,bcond',acond,bcond
	149
	150	firstcall=.false.
	151
	152	! If fixed profile of kmix
	153	IF (kmix_proffix) then
	154	file_path=trim(datadir)//'/gas_prop/kmix.txt'
	155	open(114,file=file_path,form='formatted')
	156	do ifine=1,Nfine
	157	read(114,*) levdat(ifine), kmixdat(ifine)
	158	enddo
	159	close(114)
	160	ENDIF
	161	ENDIF
	162
	163	!-----------------------------------------------------------------------
	164	! 1. initialisation
	165	! -----------------
	166
	167	nlev=nlay+1
	168
	169	! ** calcul de rhodz et dtrho/dz=dtrho*2 g/dp
	170	! avec rho=p/RT=p/ (R Theta) (p/ps)**kappa
	171	! ----------------------------------------
	172
	173	DO ilay=1,nlay
	174	DO ig=1,ngrid
	175	za(ig,ilay)=(pplev(ig,ilay)-pplev(ig,ilay+1))/g
	176	ENDDO
	177	ENDDO
	178
	179	zcst1=4.gptimestep/(r*r)
	180	DO ilev=2,nlev-1
	181	DO ig=1,ngrid
	182	zb0(ig,ilev)=pplev(ig,ilev)* &
	183	(pplev(ig,1)/pplev(ig,ilev))**rcp / &
	184	(ph(ig,ilev-1)+ph(ig,ilev))
	185	zb0(ig,ilev)=zcst1zb0(ig,ilev)zb0(ig,ilev)/ &
	186	(pplay(ig,ilev-1)-pplay(ig,ilev))
	187	! write(300,*)'zb0',zb0(ig,ilev)
	188	ENDDO
	189	ENDDO
	190	DO ig=1,ngrid
	191	zb0(ig,1)=ptimesteppplev(ig,1)/(rptsrf(ig))
	192	ENDDO
	193
	194	! ** diagnostique pour l'initialisation
	195	! ----------------------------------
	196
	197	! IF(lecrit) THEN
	198	! ig=ngrid/2+1
	199	! write(,) 'Pression (mbar) ,altitude (km),u,v,theta, rho dz'
	200	! DO ilay=1,nlay
	201	! WRITE(*,'(6f11.5)') &
	202	! .01pplay(ig,ilay),.001pzlay(ig,ilay), &
	203	! pu(ig,ilay),pv(ig,ilay),ph(ig,ilay),za(ig,ilay)
	204	! ENDDO
	205	! write(,) 'Pression (mbar) ,altitude (km),zb'
	206	! DO ilev=1,nlay
	207	! WRITE(*,'(3f15.7)') &
	208	! .01pplev(ig,ilev),.001pzlev(ig,ilev), &
	209	! zb0(ig,ilev)
	210	! ENDDO
	211	! ENDIF
	212
	213	! Potential Condensation temperature:
	214	! -----------------------------------
	215
	216	if (condensn2) then
	217	DO ilev=1,nlay
	218	DO ig=1,ngrid
	219	zhcond(ig,ilev) = &
	220	(1./(bcond-acondlog(.7143pplay(ig,ilev))))/ppopsk(ig,ilev)
	221	END DO
	222	END DO
	223	else
	224	DO ilev=1,nlay
	225	DO ig=1,ngrid
	226	zhcond(ig,ilev) = 0.
	227	END DO
	228	END DO
	229	end if
	230
	231
	232	!-----------------------------------------------------------------------
	233	! 2. ajout des tendances physiques
	234	! -----------------------------
	235
	236	DO ilev=1,nlay
	237	DO ig=1,ngrid
	238	zt(ig,ilev)=pt(ig,ilev)+pdtfi(ig,ilev)*ptimestep
	239	zu(ig,ilev)=pu(ig,ilev)+pdufi(ig,ilev)*ptimestep
	240	zv(ig,ilev)=pv(ig,ilev)+pdvfi(ig,ilev)*ptimestep
	241	zh(ig,ilev)=ph(ig,ilev)+pdhfi(ig,ilev)*ptimestep
	242	zh(ig,ilev)=max(zh(ig,ilev),zhcond(ig,ilev))
	243	ENDDO
	244	ENDDO
	245	if(tracer) then
	246	DO iq =1, nq
	247	DO ilev=1,nlay
	248	DO ig=1,ngrid
	249	zq(ig,ilev,iq)=pq(ig,ilev,iq)+pdqfi(ig,ilev,iq)*ptimestep
	250	ENDDO
	251	ENDDO
	252	ENDDO
	253	end if
	254
	255	!-----------------------------------------------------------------------
	256	! 3. schema de turbulence
	257	! --------------------
	258
	259	! ** source d'energie cinetique turbulente a la surface
	260	! (condition aux limites du schema de diffusion turbulente
	261	! dans la couche limite
	262	! ---------------------
	263	CALL vdif_cd( ngrid,nlay,pz0,g,pzlay,pu,pv,ptsrf,ph &
	264	,zcdv_true,zcdh_true)
	265	DO ig=1,ngrid
	266	zu2=pu(ig,1)pu(ig,1)+pv(ig,1)pv(ig,1)
	267	!TB16: GCM wind for flat hemisphere
	268	IF (phisfi(ig).eq.0.) zu2=max(2.,zu2)
	269
	270	zcdv(ig)=zcdv_true(ig)*sqrt(zu2)
	271	zcdh(ig)=zcdh_true(ig)*sqrt(zu2)
	272	ENDDO
	273
	274	! ** schema de diffusion turbulente dans la couche limite
	275	! ----------------------------------------------------
	276
	277	CALL vdif_kc(ngrid,nlay,ptimestep,g,pzlev,pzlay &
	278	,pu,pv,ph,zcdv_true &
	279	,pq2,zkv,zkh)
	280
	281
	282	! Adding eddy mixing to mimic 3D general circulation in 1D
	283	! RW FF 2010
	284	if ((ngrid.eq.1)) then
	285	!kmixmin is the minimum eddy mix coeff in 1D
	286
	287	! If fixed profile of kmix
	288	IF (kmix_proffix) then
	289	!! Interpolate on the model vertical grid
	290	CALL interp_line(levdat,kmixdat,Nfine,&
	291	pzlay(1,:)/1000.,kmix_z(:),nlay)
	292
	293	do ilev=1,nlay
	294	zkh(1,ilev) = max(kmix_z(ilev),zkh(1,ilev))
	295	zkv(1,ilev) = max(kmix_z(ilev),zkv(1,ilev))
	296	!zkh(1,ilev) = kmixmin
	297	!zkv(1,ilev) = kmixmin
	298	end do
	299	ELSE
	300	do ilev=1,nlay
	301	zkh(1,ilev) = max(kmixmin,zkh(1,ilev))
	302	zkv(1,ilev) = max(kmixmin,zkv(1,ilev))
	303	!zkh(1,ilev) = kmixmin
	304	!zkv(1,ilev) = kmixmin
	305	end do
	306	ENDIF
	307	endif ! ngrid.eq.1
	308
	309	!! Temporary:
	310	! zkh = zkh*0.1
	311	! zkv = zkv*0.1
	312
	313	! ** diagnostique pour le schema de turbulence
	314	! -----------------------------------------
	315
	316	! IF(lecrit) THEN
	317	! write(,) ! write(,) 'Diagnostic for the vertical turbulent mixing'
	318	! write(,) 'Cd for momentum and potential temperature'
	319
	320	! write(,) zcdv(ngrid/2+1),zcdh(ngrid/2+1)
	321	! write(,) 'Mixing coefficient for momentum and pot.temp.'
	322	! DO ilev=1,nlay
	323	! write(,) zkv(ngrid/2+1,ilev),zkh(ngrid/2+1,ilev)
	324	! ENDDO
	325	! ENDIF
	326
	327	!-----------------------------------------------------------------------
	328	! 4. inversion pour l'implicite sur u
	329	! --------------------------------
	330
	331	! ** l'equation est
	332	! u(t+1) = u(t) + dt * {(du/dt)phys}(t) + dt * {(du/dt)difv}(t+1)
	333	! avec
	334	! /zu/ = u(t) + dt * {(du/dt)phys}(t) (voir paragraphe 2.)
	335	! et
	336	! dt * {(du/dt)difv}(t+1) = dt * {(d/dz)[ Ku (du/dz) ]}(t+1)
	337	! donc les entrees sont /zcdv/ pour la condition a la limite sol
	338	! et /zkv/ = Ku
	339
	340	CALL multipl((nlay-1)*ngrid,zkv(1,2),zb0(1,2),zb(1,2))
	341	CALL multipl(ngrid,zcdv,zb0,zb)
	342
	343	DO ig=1,ngrid
	344	z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
	345	zc(ig,nlay)=za(ig,nlay)zu(ig,nlay)z1(ig)
	346	zd(ig,nlay)=zb(ig,nlay)*z1(ig)
	347	ENDDO
	348
	349	DO ilay=nlay-1,1,-1
	350	DO ig=1,ngrid
	351	z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+ &
	352	zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
	353	zc(ig,ilay)=(za(ig,ilay)*zu(ig,ilay)+ &
	354	zb(ig,ilay+1)zc(ig,ilay+1))z1(ig)
	355	zd(ig,ilay)=zb(ig,ilay)*z1(ig)
	356	ENDDO
	357	ENDDO
	358
	359	DO ig=1,ngrid
	360	zu(ig,1)=zc(ig,1)
	361	ENDDO
	362	DO ilay=2,nlay
	363	DO ig=1,ngrid
	364	zu(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*zu(ig,ilay-1)
	365	ENDDO
	366	ENDDO
	367
	368	!-----------------------------------------------------------------------
	369	! 5. inversion pour l'implicite sur v
	370	! --------------------------------
	371
	372	! ** l'equation est
	373	! v(t+1) = v(t) + dt * {(dv/dt)phys}(t) + dt * {(dv/dt)difv}(t+1)
	374	! avec
	375	! /zv/ = v(t) + dt * {(dv/dt)phys}(t) (voir paragraphe 2.)
	376	! et
	377	! dt * {(dv/dt)difv}(t+1) = dt * {(d/dz)[ Kv (dv/dz) ]}(t+1)
	378	! donc les entrees sont /zcdv/ pour la condition a la limite sol
	379	! et /zkv/ = Kv
	380
	381	DO ig=1,ngrid
	382	z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
	383	zc(ig,nlay)=za(ig,nlay)zv(ig,nlay)z1(ig)
	384	zd(ig,nlay)=zb(ig,nlay)*z1(ig)
	385	ENDDO
	386
	387	DO ilay=nlay-1,1,-1
	388	DO ig=1,ngrid
	389	z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+ &
	390	zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
	391	zc(ig,ilay)=(za(ig,ilay)*zv(ig,ilay)+ &
	392	zb(ig,ilay+1)zc(ig,ilay+1))z1(ig)
	393	zd(ig,ilay)=zb(ig,ilay)*z1(ig)
	394	ENDDO
	395	ENDDO
	396
	397	DO ig=1,ngrid
	398	zv(ig,1)=zc(ig,1)
	399	ENDDO
	400	DO ilay=2,nlay
	401	DO ig=1,ngrid
	402	zv(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*zv(ig,ilay-1)
	403	ENDDO
	404	ENDDO
	405
	406	!-----------------------------------------------------------------------
	407	! 6. inversion pour l'implicite sur h sans oublier le couplage
	408	! avec le sol (conduction)
	409	! ------------------------
	410
	411	! ** l'equation est
	412	! h(t+1) = h(t) + dt * {(dh/dt)phys}(t) + dt * {(dh/dt)difv}(t+1)
	413	! avec
	414	! /zh/ = h(t) + dt * {(dh/dt)phys}(t) (voir paragraphe 2.)
	415	! et
	416	! dt * {(dh/dt)difv}(t+1) = dt * {(d/dz)[ Kh (dh/dz) ]}(t+1)
	417	! donc les entrees sont /zcdh/ pour la condition de raccord au sol
	418	! et /zkh/ = Kh
	419	! -------------
	420
	421	CALL multipl((nlay-1)*ngrid,zkh(1,2),zb0(1,2),zb(1,2))
	422	CALL multipl(ngrid,zcdh,zb0,zb)
	423
	424	DO ig=1,ngrid
	425	z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
	426	zc(ig,nlay)=za(ig,nlay)zh(ig,nlay)z1(ig)
	427	zd(ig,nlay)=zb(ig,nlay)*z1(ig)
	428	ENDDO
	429
	430	DO ilay=nlay-1,1,-1
	431	DO ig=1,ngrid
	432	z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+ &
	433	zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
	434	zc(ig,ilay)=(za(ig,ilay)*zh(ig,ilay)+ &
	435	zb(ig,ilay+1)zc(ig,ilay+1))z1(ig)
	436	zd(ig,ilay)=zb(ig,ilay)*z1(ig)
	437	ENDDO
	438	ENDDO
	439
	440	! ** calcul de (d Planck / dT) a la temperature d'interface
	441	! ------------------------------------------------------
	442
	443	z4st=4.5.67e-8ptimestep
	444	DO ig=1,ngrid
	445	zdplanck(ig)=z4stpemis(ig)ptsrf(ig)ptsrf(ig)ptsrf(ig)
	446	ENDDO
	447
	448	! ** calcul de la temperature_d'interface et de sa tendance.
	449	! on ecrit que la somme des flux est nulle a l'interface
	450	! a t + \delta t,
	451	! !'est a dire le flux radiatif a {t + \delta t}
	452	! + le flux turbulent a {t + \delta t} &
	453	! qui 'ecrit K (T1-Tsurf) avec T1 = d1 Tsurf + c1
	454	! (notation K dt = /cpp*b/)
	455	! + le flux dans le sol a t
	456	! + l'evolution du flux dans le sol lorsque la temperature d'interface
	457	! passe de sa valeur a t a sa valeur a {t + \delta t}.
	458	! ----------------------------------------------------
	459
	460	DO ig=1,ngrid
	461	z1(ig)=pcapcal(ig)ptsrf(ig)+cppzb(ig,1)*zc(ig,1) &
	462	+zdplanck(ig)ptsrf(ig)+ pfluxsrf(ig)ptimestep
	463	z2(ig)= pcapcal(ig)+cppzb(ig,1)(1.-zd(ig,1))+zdplanck(ig)
	464	ztsrf2(ig)=z1(ig)/z2(ig)
	465	pdtsrf(ig)=(ztsrf2(ig)-ptsrf(ig))/ptimestep
	466
	467	! Modif speciale N2 condensation:
	468	! tconds = 1./(bcond-acondlog(.7143pplev(ig,1)))
	469	! if ((condensn2).and. &
	470	! ((n2ice(ig).ne.0).or.(ztsrf2(ig).lt.tconds)))then
	471	! zh(ig,1)=zc(ig,1)+zd(ig,1)*tconds
	472	! else
	473	zh(ig,1)=zc(ig,1)+zd(ig,1)*ztsrf2(ig)
	474	! end if
	475	ENDDO
	476
	477	! ** et a partir de la temperature au sol on remonte
	478	! -----------------------------------------------
	479
	480	DO ilay=2,nlay
	481	DO ig=1,ngrid
	482	hh = max( zh(ig,ilay-1) , zhcond(ig,ilay-1) ) ! modif n2cond
	483	zh(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*hh
	484	ENDDO
	485	ENDDO
	486
	487
	488	!-----------------------------------------------------------------------
	489	! TRACERS
	490	! -------
	491
	492	if(tracer) then
	493
	494	! Using the wind modified by friction for lifting and sublimation
	495	! ----------------------------------------------------------------
	496	! This is computed above
	497
	498	! DO ig=1,ngrid
	499	! zu2=zu(ig,1)zu(ig,1)+zv(ig,1)zv(ig,1)
	500	! zcdv(ig)=zcdv_true(ig)*sqrt(zu2)
	501	! zcdh(ig)=zcdh_true(ig)*sqrt(zu2)
	502	! ENDDO
	503
	504	! Calcul flux vertical au bas de la premiere couche (cf dust on Mars)
	505	! -----------------------------------------------------------
	506	do ig=1,ngrid
	507	rho(ig) = zb0(ig,1) /ptimestep
	508	! zb(ig,1) = 0.
	509	end do
	510
	511	pdqsdif(:,:) = 0
	512	pdqdif(:,:,:)=0.
	513
	514
	515	! TB: Eddy lifting of tracers :
	516	! ****************************************************************
	517	! DO ig=1,ngrid
	518	!! option : injection only on an equatorial band
	519	!! if (abs(lati(ig))*180./pi.le.25.) then
	520	! pdqsdif(ig,igcm_eddy1e6) =-1.e-12
	521	! pdqsdif(ig,igcm_eddy1e7) =-1.e-12
	522	! pdqsdif(ig,igcm_eddy5e7) =-1.e-12
	523	! pdqsdif(ig,igcm_eddy1e8) =-1.e-12
	524	! pdqsdif(ig,igcm_eddy5e8) =-1.e-12
	525	! endif
	526	! ENDDO
	527
	528	! pdqdifeddy(:,:,:)=0.
	529	! injection a 50 km
	530	! DO ig=1,ngrid
	531	! pdqdifeddy(ig,17,igcm_eddy1e6)=1e-12
	532	! pdqdifeddy(ig,17,igcm_eddy1e7)=1e-12
	533	! pdqdifeddy(ig,17,igcm_eddy5e7)=1e-12
	534	! pdqdifeddy(ig,17,igcm_eddy1e8)=1e-12
	535	! pdqdifeddy(ig,17,igcm_eddy5e8)=1e-12
	536	! ENDDO
	537
	538	! Inversion pour l'implicite sur q
	539	! --------------------------------
	540	do iq=1,nq
	541	CALL multipl((nlay-1)*ngrid,zkh(1,2),zb0(1,2),zb(1,2))
	542
	543	if ((methane).and.(iq.eq.igcm_ch4_gas)) then
	544	! This line is required to account for turbulent transport
	545	! from surface (e.g. ice) to mid-layer of atmosphere:
	546	CALL multipl(ngrid,zcdv,zb0,zb(1,1))
	547	else if ((carbox).and.(iq.eq.igcm_co_gas)) then
	548	CALL multipl(ngrid,zcdv,zb0,zb(1,1))
	549	else ! (re)-initialize zb(:,1)
	550	zb(1:ngrid,1)=0
	551	end if
	552
	553	DO ig=1,ngrid
	554	z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
	555	zc(ig,nlay)=za(ig,nlay)zq(ig,nlay,iq)z1(ig)
	556	zd(ig,nlay)=zb(ig,nlay)*z1(ig)
	557	ENDDO
	558
	559	DO ilay=nlay-1,2,-1
	560	DO ig=1,ngrid
	561	z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+ &
	562	zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
	563	zc(ig,ilay)=(za(ig,ilay)*zq(ig,ilay,iq)+ &
	564	zb(ig,ilay+1)zc(ig,ilay+1))z1(ig)
	565	zd(ig,ilay)=zb(ig,ilay)*z1(ig)
	566	ENDDO
	567	ENDDO
	568
	569	! special case for methane and CO ice tracer: do not include
	570	! ice tracer from surface (which is set when handling
	571	! vapour case (see further down).
	572	if (methane.and.(iq.eq.igcm_ch4_ice)) then
	573	DO ig=1,ngrid
	574	z1(ig)=1./(za(ig,1)+zb(ig,1)+ &
	575	zb(ig,2)*(1.-zd(ig,2)))
	576	zc(ig,1)=(za(ig,1)*zq(ig,1,iq)+ &
	577	zb(ig,2)zc(ig,2)) z1(ig)
	578	ENDDO
	579	else if (carbox.and.(iq.eq.igcm_co_ice)) then
	580	DO ig=1,ngrid
	581	z1(ig)=1./(za(ig,1)+zb(ig,1)+ &
	582	zb(ig,2)*(1.-zd(ig,2)))
	583	zc(ig,1)=(za(ig,1)*zq(ig,1,iq)+ &
	584	zb(ig,2)zc(ig,2)) z1(ig)
	585	ENDDO
	586
	587	else ! general case
	588	DO ig=1,ngrid
	589	z1(ig)=1./(za(ig,1)+zb(ig,1)+ &
	590	zb(ig,2)*(1.-zd(ig,2)))
	591	zc(ig,1)=(za(ig,1)*zq(ig,1,iq)+ &
	592	zb(ig,2)*zc(ig,2) + &
	593	(-pdqsdif(ig,iq)) ptimestep) z1(ig) !tracer flux from surface
	594	ENDDO
	595	endif ! of if (methane.and.(iq.eq.igcm_ch4_ice))
	596
	597	! Calculation for turbulent exchange with the surface (for ice)
	598	IF (methane.and.(iq.eq.igcm_ch4_gas)) then
	599
	600	!! calcul de la valeur de q a la surface :
	601	call methanesat(ngrid,ptsrf,pplev(1,1), &
	602	qsat_ch4(:),pqsurf(:,igcm_n2))
	603
	604	!! For output:
	605	call methanesat(ngrid,zt(:,1),pplev(1,1), &
	606	qsat_ch4_l1(:),pqsurf(:,igcm_n2))
	607
	608	!! Prevent CH4 condensation at the surface
	609	if (.not.condmetsurf) then
	610	qsat_ch4=qsat_ch4*1.e6
	611	endif
	612
	613	DO ig=1,ngrid
	614	zd(ig,1)=zb(ig,1)*z1(ig)
	615	zq1temp_ch4(ig)=zc(ig,1)+ zd(ig,1)*qsat_ch4(ig)
	616	pdqsdif(ig,igcm_ch4_ice)=rho(ig)*zcdv(ig) &
	617	*(zq1temp_ch4(ig)-qsat_ch4(ig))
	618	END DO
	619
	620	DO ig=1,ngrid
	621	if ((-pdqsdif(ig,igcm_ch4_ice)*ptimestep) &
	622	.gt.(pqsurf(ig,igcm_ch4_ice))) then
	623
	624	pdqsdif(ig,igcm_ch4_ice)= &
	625	-pqsurf(ig,igcm_ch4_ice)/ptimestep
	626
	627	z1(ig)=1./(za(ig,1)+ zb(ig,2)*(1.-zd(ig,2)))
	628
	629	zc(ig,1)=(za(ig,1)*zq(ig,1,igcm_ch4_gas)+ &
	630	zb(ig,2)*zc(ig,2) + &
	631	(-pdqsdif(ig,igcm_ch4_ice)) ptimestep) z1(ig)
	632
	633	zq1temp_ch4(ig)=zc(ig,1)
	634	endif
	635
	636	zq(ig,1,igcm_ch4_gas)=zq1temp_ch4(ig)
	637
	638	! TB: MODIF speciale pour CH4
	639	pdtsrf(ig)=pdtsrf(ig)+ &
	640	lw_ch4*pdqsdif(ig,igcm_ch4_ice)/pcapcal(ig)
	641
	642
	643	ENDDO ! of DO ig=1,ngrid
	644
	645	ELSE IF (carbox.and.(iq.eq.igcm_co_gas)) then
	646
	647	! calcul de la valeur de q a la surface :
	648	call cosat(ngrid,ptsrf,pplev(1,1),qsat_co, &
	649	pqsurf(:,igcm_n2),pqsurf(:,igcm_ch4_ice))
	650
	651	!! Prevent CO condensation at the surface
	652	if (.not.condcosurf) then
	653	qsat_co=qsat_co*1.e6
	654	endif
	655
	656	DO ig=1,ngrid
	657	zd(ig,1)=zb(ig,1)*z1(ig)
	658	zq1temp_co(ig)=zc(ig,1)+ zd(ig,1)*qsat_co(ig)
	659	pdqsdif(ig,igcm_co_ice)=rho(ig)*zcdv(ig) &
	660	*(zq1temp_co(ig)-qsat_co(ig))
	661	END DO
	662
	663
	664	DO ig=1,ngrid
	665	! -------------------------------------------------------------
	666	! TEMPORAIRE : pour initialiser CO si glacier N2
	667	! meme si il n'y a pas de CO disponible
	668	! if (pqsurf(ig,igcm_n2).le.10.) then
	669	! -------------------------------------------------------------
	670	!
	671	if ((-pdqsdif(ig,igcm_co_ice)*ptimestep) &
	672	.gt.(pqsurf(ig,igcm_co_ice))) then
	673	pdqsdif(ig,igcm_co_ice)= &
	674	-pqsurf(ig,igcm_co_ice)/ptimestep
	675	z1(ig)=1./(za(ig,1)+ zb(ig,2)*(1.-zd(ig,2)))
	676	zc(ig,1)=(za(ig,1)*zq(ig,1,igcm_co_gas)+ &
	677	zb(ig,2)*zc(ig,2) + &
	678	(-pdqsdif(ig,igcm_co_ice)) ptimestep) z1(ig)
	679	zq1temp_co(ig)=zc(ig,1)
	680	endif
	681	! endif
	682
	683	zq(ig,1,igcm_co_gas)=zq1temp_co(ig)
	684
	685	! MODIF speciale pour CO / corrected by FF 2016
	686	pdtsrf(ig)=pdtsrf(ig)+ &
	687	lw_co*pdqsdif(ig,igcm_co_ice)/pcapcal(ig)
	688
	689	ENDDO ! of DO ig=1,ngrid
	690
	691	ELSE ! if (methane)
	692
	693	DO ig=1,ngrid
	694	zq(ig,1,iq)=zc(ig,1)
	695	ENDDO
	696
[3258]	697	END IF ! of IF ((methane).and.(iq.eq.igcm_ch4_gas))
[3247]	698
	699	!! Diffusion verticale : shut down vertical transport of vertdiff = false
	700	if (vertdiff) then
	701	DO ilay=2,nlay
	702	DO ig=1,ngrid
	703	zq(ig,ilay,iq)=zc(ig,ilay)+zd(ig,ilay)*zq(ig,ilay-1,iq)
	704	ENDDO
	705	ENDDO
	706	endif
	707
	708	enddo ! of do iq=1,nq
	709	end if ! of if(tracer)
	710
	711	!-----------------------------------------------------------------------
	712	! 8. calcul final des tendances de la diffusion verticale
	713	! ----------------------------------------------------
	714	DO ilev = 1, nlay
	715	DO ig=1,ngrid
	716	pdudif(ig,ilev)=( zu(ig,ilev)- &
	717	(pu(ig,ilev)+pdufi(ig,ilev)*ptimestep) )/ptimestep
	718	pdvdif(ig,ilev)=( zv(ig,ilev)- &
	719	(pv(ig,ilev)+pdvfi(ig,ilev)*ptimestep) )/ptimestep
	720	hh = max(ph(ig,ilev)+pdhfi(ig,ilev)*ptimestep , &
	721	zhcond(ig,ilev)) ! modif n2cond
	722	pdhdif(ig,ilev)=( zh(ig,ilev)- hh )/ptimestep
	723	ENDDO
	724	ENDDO
	725
	726	if (tracer) then
	727	DO iq = 1, nq
	728	DO ilev = 1, nlay
	729	DO ig=1,ngrid
	730	pdqdif(ig,ilev,iq)=(zq(ig,ilev,iq)- &
	731	(pq(ig,ilev,iq) + pdqfi(ig,ilev,iq)*ptimestep))/ptimestep
	732	! pdqdif(ig,ilev,iq)=pdqdifeddy(ig,ilev,iq)+(zq(ig,ilev,iq)- &
	733	! (pq(ig,ilev,iq) + pdqfi(ig,ilev,iq)*ptimestep))/ptimestep
	734	ENDDO
	735	ENDDO
	736	ENDDO
	737	end if
	738
	739	! ** diagnostique final
	740	! ------------------
	741
	742	IF(lecrit) THEN
	743	write(,) 'In vdif'
	744	write(,) 'Ts (t) and Ts (t+st)'
	745	WRITE(*,'(a10,3a15)') &
	746	'theta(t)','theta(t+dt)','u(t)','u(t+dt)'
	747	write(,) ptsrf(ngrid/2+1),ztsrf2(ngrid/2+1)
	748	DO ilev=1,nlay
	749	WRITE(*,'(4f15.7)') &
	750	ph(ngrid/2+1,ilev),zh(ngrid/2+1,ilev), &
	751	pu(ngrid/2+1,ilev),zu(ngrid/2+1,ilev)
	752
	753	ENDDO
	754	ENDIF
	755
	756	RETURN
	757	! END
	758
	759	end subroutine vdifc_pluto
	760
	761	end module vdifc_pluto_mod

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

Download in other formats:

Original Format