Context Navigation

vdifc.F @ 3436

Last change on this file since 3436 was 3237, checked in by afalco, 9 months ago
Pluto PCM: Import methane and ch4 clouds from pluto.old AF
File size: 23.9 KB

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

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

Download in other formats:

Original Format