Context Navigation

vdifc.F @ 295

Last change on this file since 295 was 291, checked in by acolaitis, 14 years ago
Sorry, removed the .eq. .true. to suit picky compilers. This is a bad habit from IDL coding...
File size: 21.9 KB

Rev	Line
[38]	1	SUBROUTINE vdifc(ngrid,nlay,nq,co2ice,ppopsk,
	2	$ ptimestep,pcapcal,lecrit,
	3	$ pplay,pplev,pzlay,pzlev,pz0,
	4	$ pu,pv,ph,pq,ptsrf,pemis,pqsurf,
	5	$ pdufi,pdvfi,pdhfi,pdqfi,pfluxsrf,
	6	$ pdudif,pdvdif,pdhdif,pdtsrf,pq2,
[268]	7	$ pdqdif,pdqsdif,wmax,zcdv_true,zcdh_true)
[38]	8	IMPLICIT NONE
	9
	10	c=======================================================================
	11	c
	12	c subject:
	13	c --------
	14	c Turbulent diffusion (mixing) for potential T, U, V and tracer
	15	c
	16	c Shema implicite
	17	c On commence par rajouter au variables x la tendance physique
	18	c et on resoult en fait:
	19	c x(t+1) = x(t) + dt * (dx/dt)phys(t) + dt * (dx/dt)difv(t+1)
	20	c
	21	c author:
	22	c ------
	23	c Hourdin/Forget/Fournier
	24	c=======================================================================
	25
	26	c-----------------------------------------------------------------------
	27	c declarations:
	28	c -------------
	29
	30	#include "dimensions.h"
	31	#include "dimphys.h"
	32	#include "comcstfi.h"
	33	#include "callkeys.h"
	34	#include "surfdat.h"
	35	#include "comgeomfi.h"
	36	#include "tracer.h"
	37
	38	c
	39	c arguments:
	40	c ----------
	41
	42	INTEGER ngrid,nlay
	43	REAL ptimestep
	44	REAL pplay(ngrid,nlay),pplev(ngrid,nlay+1)
	45	REAL pzlay(ngrid,nlay),pzlev(ngrid,nlay+1)
	46	REAL pu(ngrid,nlay),pv(ngrid,nlay),ph(ngrid,nlay)
	47	REAL ptsrf(ngrid),pemis(ngrid)
	48	REAL pdufi(ngrid,nlay),pdvfi(ngrid,nlay),pdhfi(ngrid,nlay)
	49	REAL pfluxsrf(ngrid)
	50	REAL pdudif(ngrid,nlay),pdvdif(ngrid,nlay),pdhdif(ngrid,nlay)
	51	REAL pdtsrf(ngrid),pcapcal(ngrid)
	52	REAL pq2(ngrid,nlay+1)
	53
	54	c Argument added for condensation:
	55	REAL co2ice (ngrid), ppopsk(ngrid,nlay)
	56	logical lecrit
	57
[224]	58	REAL pz0(ngridmx) ! surface roughness length (m)
	59
[256]	60	c Argument added to account for subgrid gustiness :
	61
	62	REAL wmax(ngridmx)
	63
[38]	64	c Traceurs :
	65	integer nq
	66	REAL pqsurf(ngrid,nq)
	67	real pq(ngrid,nlay,nq), pdqfi(ngrid,nlay,nq)
	68	real pdqdif(ngrid,nlay,nq)
	69	real pdqsdif(ngrid,nq)
	70
	71	c local:
	72	c ------
	73
	74	INTEGER ilev,ig,ilay,nlev
	75
	76	REAL z4st,zdplanck(ngridmx)
	77	REAL zkv(ngridmx,nlayermx+1),zkh(ngridmx,nlayermx+1)
	78	REAL zcdv(ngridmx),zcdh(ngridmx)
[268]	79	REAL zcdv_true(ngridmx),zcdh_true(ngridmx) ! drag coeff are used by the LES to recompute u* and hfx
[38]	80	REAL zu(ngridmx,nlayermx),zv(ngridmx,nlayermx)
	81	REAL zh(ngridmx,nlayermx)
	82	REAL ztsrf2(ngridmx)
	83	REAL z1(ngridmx),z2(ngridmx)
	84	REAL za(ngridmx,nlayermx),zb(ngridmx,nlayermx)
	85	REAL zb0(ngridmx,nlayermx)
	86	REAL zc(ngridmx,nlayermx),zd(ngridmx,nlayermx)
	87	REAL zcst1
[284]	88	REAL zu2(ngridmx)
[38]	89
	90	EXTERNAL SSUM,SCOPY
	91	REAL SSUM
	92	LOGICAL firstcall
	93	SAVE firstcall
	94
	95	c variable added for CO2 condensation:
	96	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	97	REAL hh , zhcond(ngridmx,nlayermx)
	98	REAL latcond,tcond1mb
	99	REAL acond,bcond
	100	SAVE acond,bcond
	101	DATA latcond,tcond1mb/5.9e5,136.27/
	102
	103	c Tracers :
	104	c ~~~~~~~
	105	INTEGER iq
	106	REAL zq(ngridmx,nlayermx,nqmx)
	107	REAL zq1temp(ngridmx)
	108	REAL rho(ngridmx) ! near surface air density
	109	REAL qsat(ngridmx)
	110	DATA firstcall/.true./
	111
	112	REAL kmixmin
	113
	114	c ** un petit test de coherence
	115	c --------------------------
	116
	117	IF (firstcall) THEN
	118	IF(ngrid.NE.ngridmx) THEN
	119	PRINT*,'STOP dans vdifc'
	120	PRINT*,'probleme de dimensions :'
	121	PRINT*,'ngrid =',ngrid
	122	PRINT*,'ngridmx =',ngridmx
	123	STOP
	124	ENDIF
	125	c To compute: Tcond= 1./(bcond-acondlog(.0095p)) (p in pascal)
	126	bcond=1./tcond1mb
	127	acond=r/latcond
	128	PRINT*,'In vdifc: Tcond(P=1mb)=',tcond1mb,' Lcond=',latcond
	129	PRINT*,' acond,bcond',acond,bcond
	130
	131	firstcall=.false.
	132	ENDIF
	133
	134
	135
	136
	137
	138	c-----------------------------------------------------------------------
	139	c 1. initialisation
	140	c -----------------
	141
	142	nlev=nlay+1
	143
	144	c ** calcul de rhodz et dtrho/dz=dtrho*2 g/dp
	145	c avec rho=p/RT=p/ (R Theta) (p/ps)**kappa
	146	c ----------------------------------------
	147
	148	DO ilay=1,nlay
	149	DO ig=1,ngrid
	150	za(ig,ilay)=(pplev(ig,ilay)-pplev(ig,ilay+1))/g
	151	ENDDO
	152	ENDDO
	153
	154	zcst1=4.gptimestep/(r*r)
	155	DO ilev=2,nlev-1
	156	DO ig=1,ngrid
	157	zb0(ig,ilev)=pplev(ig,ilev)*
	158	s (pplev(ig,1)/pplev(ig,ilev))**rcp /
	159	s (ph(ig,ilev-1)+ph(ig,ilev))
	160	zb0(ig,ilev)=zcst1zb0(ig,ilev)zb0(ig,ilev)/
	161	s (pplay(ig,ilev-1)-pplay(ig,ilev))
	162	ENDDO
	163	ENDDO
	164	DO ig=1,ngrid
	165	zb0(ig,1)=ptimesteppplev(ig,1)/(rptsrf(ig))
	166	ENDDO
	167
	168	c ** diagnostique pour l'initialisation
	169	c ----------------------------------
	170
	171	IF(lecrit) THEN
	172	ig=ngrid/2+1
	173	PRINT*,'Pression (mbar) ,altitude (km),u,v,theta, rho dz'
	174	DO ilay=1,nlay
	175	WRITE(*,'(6f11.5)')
	176	s .01pplay(ig,ilay),.001pzlay(ig,ilay),
	177	s pu(ig,ilay),pv(ig,ilay),ph(ig,ilay),za(ig,ilay)
	178	ENDDO
	179	PRINT*,'Pression (mbar) ,altitude (km),zb'
	180	DO ilev=1,nlay
	181	WRITE(*,'(3f15.7)')
	182	s .01pplev(ig,ilev),.001pzlev(ig,ilev),
	183	s zb0(ig,ilev)
	184	ENDDO
	185	ENDIF
	186
	187	c Potential Condensation temperature:
	188	c -----------------------------------
	189
	190	c if (callcond) then
	191	c DO ilev=1,nlay
	192	c DO ig=1,ngrid
	193	c zhcond(ig,ilev) =
	194	c & (1./(bcond-acondlog(.0095pplay(ig,ilev))))/ppopsk(ig,ilev)
	195	c END DO
	196	c END DO
	197	c else
	198	call zerophys(ngrid*nlay,zhcond)
	199	c end if
	200
	201
	202	c-----------------------------------------------------------------------
	203	c 2. ajout des tendances physiques
	204	c -----------------------------
	205
	206	DO ilev=1,nlay
	207	DO ig=1,ngrid
	208	zu(ig,ilev)=pu(ig,ilev)+pdufi(ig,ilev)*ptimestep
	209	zv(ig,ilev)=pv(ig,ilev)+pdvfi(ig,ilev)*ptimestep
	210	zh(ig,ilev)=ph(ig,ilev)+pdhfi(ig,ilev)*ptimestep
	211	zh(ig,ilev)=max(zh(ig,ilev),zhcond(ig,ilev))
	212	ENDDO
	213	ENDDO
	214	if(tracer) then
	215	DO iq =1, nq
	216	DO ilev=1,nlay
	217	DO ig=1,ngrid
	218	zq(ig,ilev,iq)=pq(ig,ilev,iq)+pdqfi(ig,ilev,iq)*ptimestep
	219	ENDDO
	220	ENDDO
	221	ENDDO
	222	end if
	223
	224	c-----------------------------------------------------------------------
	225	c 3. schema de turbulence
	226	c --------------------
	227
	228	c ** source d'energie cinetique turbulente a la surface
	229	c (condition aux limites du schema de diffusion turbulente
	230	c dans la couche limite
	231	c ---------------------
	232
[268]	233	CALL vdif_cd(ngrid,nlay,pz0,g,pzlay,pu,pv,wmax,ptsrf,ph
[38]	234	& ,zcdv_true,zcdh_true)
[256]	235
[290]	236	zu2(:)=pu(:,1)pu(:,1)+pv(:,1)pv(:,1)
[256]	237
[291]	238	IF (callrichsl) THEN
[290]	239	zcdv(:)=zcdv_true(:)sqrt(zu2(:)+(0.7wmax(:))**2) ! subgrid gustiness added by quadratic interpolation with a coeff beta
	240	zcdh(:)=zcdh_true(:)sqrt(zu2(:)+(1.2wmax(:))**2) ! LES comparisons. This parameter is linked to the thermals settings)
[284]	241	ELSE
	242	zcdv(:)=zcdv_true(:)*sqrt(zu2(:)) ! 1 / bulk aerodynamic momentum conductance
	243	zcdh(:)=zcdh_true(:)*sqrt(zu2(:)) ! 1 / bulk aerodynamic heat conductance
[290]	244	ENDIF
[38]	245
[290]	246
[256]	247	! Some usefull diagnostics for the new surface layer parametrization :
	248
[268]	249	! call WRITEDIAGFI(ngridmx,'pcdv',
[256]	250	! & 'momentum drag','no units',
	251	! & 2,zcdv_true)
[268]	252	! call WRITEDIAGFI(ngridmx,'pcdh',
[256]	253	! & 'heat drag','no units',
	254	! & 2,zcdh_true)
[268]	255	! call WRITEDIAGFI(ngridmx,'u*',
[284]	256	! & 'friction velocity','m/s',
[290]	257	! & 2,sqrt(zcdv_true(:))sqrt(zu2(:)+(0.7wmax(:))**2))
[268]	258	! call WRITEDIAGFI(ngridmx,'T*',
	259	! & 'friction temperature','K',
[290]	260	! & 2,zcdh_true(:)sqrt(zu2(:)+(1.2wmax(:))*2)
	261	! & -(ptsrf(:)-ph(:,1))
	262	! & /(sqrt(zcdv_true(:))sqrt(zu2(:)+(0.7wmax(:))**2)))
[268]	263	! call WRITEDIAGFI(ngridmx,'rm-1',
	264	! & 'aerodyn momentum conductance','m/s',
[256]	265	! & 2,zcdv)
[268]	266	! call WRITEDIAGFI(ngridmx,'rh-1',
	267	! & 'aerodyn heat conductance','m/s',
[256]	268	! & 2,zcdh)
	269
[38]	270	c ** schema de diffusion turbulente dans la couche limite
	271	c ----------------------------------------------------
	272
	273	CALL vdif_kc(ptimestep,g,pzlev,pzlay
	274	& ,pu,pv,ph,zcdv_true
	275	& ,pq2,zkv,zkh)
	276
	277	if ((doubleq).and.(ngrid.eq.1)) then
	278	kmixmin = 80. !80.! minimum eddy mix coeff in 1D
	279	do ilev=1,nlay
	280	do ig=1,ngrid
	281	zkh(ig,ilev) = max(kmixmin,zkh(ig,ilev))
	282	zkv(ig,ilev) = max(kmixmin,zkv(ig,ilev))
	283	end do
	284	end do
	285	end if
	286
	287	c ** diagnostique pour le schema de turbulence
	288	c -----------------------------------------
	289
	290	IF(lecrit) THEN
	291	PRINT*
	292	PRINT*,'Diagnostic for the vertical turbulent mixing'
	293	PRINT*,'Cd for momentum and potential temperature'
	294
	295	PRINT*,zcdv(ngrid/2+1),zcdh(ngrid/2+1)
	296	PRINT*,'Mixing coefficient for momentum and pot.temp.'
	297	DO ilev=1,nlay
	298	PRINT*,zkv(ngrid/2+1,ilev),zkh(ngrid/2+1,ilev)
	299	ENDDO
	300	ENDIF
	301
	302
	303
	304
	305	c-----------------------------------------------------------------------
	306	c 4. inversion pour l'implicite sur u
	307	c --------------------------------
	308
	309	c ** l'equation est
	310	c u(t+1) = u(t) + dt * {(du/dt)phys}(t) + dt * {(du/dt)difv}(t+1)
	311	c avec
	312	c /zu/ = u(t) + dt * {(du/dt)phys}(t) (voir paragraphe 2.)
	313	c et
	314	c dt * {(du/dt)difv}(t+1) = dt * {(d/dz)[ Ku (du/dz) ]}(t+1)
	315	c donc les entrees sont /zcdv/ pour la condition a la limite sol
	316	c et /zkv/ = Ku
	317
	318	CALL multipl((nlay-1)*ngrid,zkv(1,2),zb0(1,2),zb(1,2))
	319	CALL multipl(ngrid,zcdv,zb0,zb)
	320
	321	DO ig=1,ngrid
	322	z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
	323	zc(ig,nlay)=za(ig,nlay)zu(ig,nlay)z1(ig)
	324	zd(ig,nlay)=zb(ig,nlay)*z1(ig)
	325	ENDDO
	326
	327	DO ilay=nlay-1,1,-1
	328	DO ig=1,ngrid
	329	z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
	330	$ zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
	331	zc(ig,ilay)=(za(ig,ilay)*zu(ig,ilay)+
	332	$ zb(ig,ilay+1)zc(ig,ilay+1))z1(ig)
	333	zd(ig,ilay)=zb(ig,ilay)*z1(ig)
	334	ENDDO
	335	ENDDO
	336
	337	DO ig=1,ngrid
	338	zu(ig,1)=zc(ig,1)
	339	ENDDO
	340	DO ilay=2,nlay
	341	DO ig=1,ngrid
	342	zu(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*zu(ig,ilay-1)
	343	ENDDO
	344	ENDDO
	345
	346
	347
	348
	349
	350	c-----------------------------------------------------------------------
	351	c 5. inversion pour l'implicite sur v
	352	c --------------------------------
	353
	354	c ** l'equation est
	355	c v(t+1) = v(t) + dt * {(dv/dt)phys}(t) + dt * {(dv/dt)difv}(t+1)
	356	c avec
	357	c /zv/ = v(t) + dt * {(dv/dt)phys}(t) (voir paragraphe 2.)
	358	c et
	359	c dt * {(dv/dt)difv}(t+1) = dt * {(d/dz)[ Kv (dv/dz) ]}(t+1)
	360	c donc les entrees sont /zcdv/ pour la condition a la limite sol
	361	c et /zkv/ = Kv
	362
	363	DO ig=1,ngrid
	364	z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
	365	zc(ig,nlay)=za(ig,nlay)zv(ig,nlay)z1(ig)
	366	zd(ig,nlay)=zb(ig,nlay)*z1(ig)
	367	ENDDO
	368
	369	DO ilay=nlay-1,1,-1
	370	DO ig=1,ngrid
	371	z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
	372	$ zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
	373	zc(ig,ilay)=(za(ig,ilay)*zv(ig,ilay)+
	374	$ zb(ig,ilay+1)zc(ig,ilay+1))z1(ig)
	375	zd(ig,ilay)=zb(ig,ilay)*z1(ig)
	376	ENDDO
	377	ENDDO
	378
	379	DO ig=1,ngrid
	380	zv(ig,1)=zc(ig,1)
	381	ENDDO
	382	DO ilay=2,nlay
	383	DO ig=1,ngrid
	384	zv(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*zv(ig,ilay-1)
	385	ENDDO
	386	ENDDO
	387
	388
	389
	390
	391
	392	c-----------------------------------------------------------------------
	393	c 6. inversion pour l'implicite sur h sans oublier le couplage
	394	c avec le sol (conduction)
	395	c ------------------------
	396
	397	c ** l'equation est
	398	c h(t+1) = h(t) + dt * {(dh/dt)phys}(t) + dt * {(dh/dt)difv}(t+1)
	399	c avec
	400	c /zh/ = h(t) + dt * {(dh/dt)phys}(t) (voir paragraphe 2.)
	401	c et
	402	c dt * {(dh/dt)difv}(t+1) = dt * {(d/dz)[ Kh (dh/dz) ]}(t+1)
	403	c donc les entrees sont /zcdh/ pour la condition de raccord au sol
	404	c et /zkh/ = Kh
	405	c -------------
	406
	407	CALL multipl((nlay-1)*ngrid,zkh(1,2),zb0(1,2),zb(1,2))
	408	CALL multipl(ngrid,zcdh,zb0,zb)
	409
	410	DO ig=1,ngrid
	411	z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
	412	zc(ig,nlay)=za(ig,nlay)zh(ig,nlay)z1(ig)
	413	zd(ig,nlay)=zb(ig,nlay)*z1(ig)
	414	ENDDO
	415
	416	DO ilay=nlay-1,1,-1
	417	DO ig=1,ngrid
	418	z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
	419	$ zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
	420	zc(ig,ilay)=(za(ig,ilay)*zh(ig,ilay)+
	421	$ zb(ig,ilay+1)zc(ig,ilay+1))z1(ig)
	422	zd(ig,ilay)=zb(ig,ilay)*z1(ig)
	423	ENDDO
	424	ENDDO
	425
	426	c ** calcul de (d Planck / dT) a la temperature d'interface
	427	c ------------------------------------------------------
	428
	429	z4st=4.5.67e-8ptimestep
	430	DO ig=1,ngrid
	431	zdplanck(ig)=z4stpemis(ig)ptsrf(ig)ptsrf(ig)ptsrf(ig)
	432	ENDDO
	433
	434	c ** calcul de la temperature_d'interface et de sa tendance.
	435	c on ecrit que la somme des flux est nulle a l'interface
	436	c a t + \delta t,
	437	c c'est a dire le flux radiatif a {t + \delta t}
	438	c + le flux turbulent a {t + \delta t}
	439	c qui s'ecrit K (T1-Tsurf) avec T1 = d1 Tsurf + c1
	440	c (notation K dt = /cpp*b/)
	441	c + le flux dans le sol a t
	442	c + l'evolution du flux dans le sol lorsque la temperature d'interface
	443	c passe de sa valeur a t a sa valeur a {t + \delta t}.
	444	c ----------------------------------------------------
	445
	446	DO ig=1,ngrid
	447	z1(ig)=pcapcal(ig)ptsrf(ig)+cppzb(ig,1)*zc(ig,1)
	448	s +zdplanck(ig)ptsrf(ig)+ pfluxsrf(ig)ptimestep
	449	z2(ig)= pcapcal(ig)+cppzb(ig,1)(1.-zd(ig,1))+zdplanck(ig)
	450	ztsrf2(ig)=z1(ig)/z2(ig)
	451	pdtsrf(ig)=(ztsrf2(ig)-ptsrf(ig))/ptimestep
	452
	453	c Modif speciale CO2 condensation:
	454	c tconds = 1./(bcond-acondlog(.0095pplev(ig,1)))
	455	c if ((callcond).and.
	456	c & ((co2ice(ig).ne.0).or.(ztsrf2(ig).lt.tconds)))then
	457	c zh(ig,1)=zc(ig,1)+zd(ig,1)*tconds
	458	c else
	459	zh(ig,1)=zc(ig,1)+zd(ig,1)*ztsrf2(ig)
	460	c end if
	461	ENDDO
	462
	463	c ** et a partir de la temperature au sol on remonte
	464	c -----------------------------------------------
	465
	466	DO ilay=2,nlay
	467	DO ig=1,ngrid
	468	hh = max( zh(ig,ilay-1) , zhcond(ig,ilay-1) ) ! modif co2cond
	469	zh(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*hh
	470	ENDDO
	471	ENDDO
	472
	473
	474	c-----------------------------------------------------------------------
	475	c TRACERS
	476	c -------
	477
	478	if(tracer) then
	479
	480	c Using the wind modified by friction for lifting and sublimation
	481	c ----------------------------------------------------------------
	482	DO ig=1,ngrid
[284]	483	zu2(ig)=zu(ig,1)zu(ig,1)+zv(ig,1)zv(ig,1)
	484	zcdv(ig)=zcdv_true(ig)*sqrt(zu2(ig))
	485	zcdh(ig)=zcdh_true(ig)*sqrt(zu2(ig))
[38]	486	ENDDO
	487
	488	c Calcul du flux vertical au bas de la premiere couche (dust) :
	489	c -----------------------------------------------------------
	490	do ig=1,ngridmx
	491	rho(ig) = zb0(ig,1) /ptimestep
	492	c zb(ig,1) = 0.
	493	end do
	494	c Dust lifting:
	495	if (lifting) then
	496	if (doubleq.AND.submicron) then
	497	do ig=1,ngrid
	498	c if(co2ice(ig).lt.1) then
	499	pdqsdif(ig,igcm_dust_mass) =
	500	& -alpha_lift(igcm_dust_mass)
	501	pdqsdif(ig,igcm_dust_number) =
	502	& -alpha_lift(igcm_dust_number)
	503	pdqsdif(ig,igcm_dust_submicron) =
	504	& -alpha_lift(igcm_dust_submicron)
	505	c end if
	506	end do
	507	else if (doubleq) then
	508	do ig=1,ngrid
[77]	509	!!! soulevement constant
[38]	510	pdqsdif(ig,igcm_dust_mass) =
	511	& -alpha_lift(igcm_dust_mass)
	512	pdqsdif(ig,igcm_dust_number) =
	513	& -alpha_lift(igcm_dust_number)
	514	end do
	515	else if (submicron) then
	516	do ig=1,ngrid
	517	pdqsdif(ig,igcm_dust_submicron) =
	518	& -alpha_lift(igcm_dust_submicron)
	519	end do
	520	else
	521	call dustlift(ngrid,nlay,nq,rho,zcdh_true,zcdh,co2ice,
	522	& pdqsdif)
	523	endif !doubleq.AND.submicron
	524	else
	525	pdqsdif(1:ngrid,1:nq) = 0.
	526	end if
	527
	528	c OU calcul de la valeur de q a la surface (water) :
	529	c ----------------------------------------
	530	if (water) then
	531	call watersat(ngridmx,ptsrf,pplev(1,1),qsat)
	532	end if
	533
	534	c Inversion pour l'implicite sur q
	535	c --------------------------------
	536	do iq=1,nq
	537	CALL multipl((nlay-1)*ngrid,zkh(1,2),zb0(1,2),zb(1,2))
	538
	539	if ((water).and.(iq.eq.igcm_h2o_vap)) then
	540	c This line is required to account for turbulent transport
	541	c from surface (e.g. ice) to mid-layer of atmosphere:
	542	CALL multipl(ngrid,zcdv,zb0,zb(1,1))
	543	CALL multipl(ngrid,dryness,zb(1,1),zb(1,1))
	544	else ! (re)-initialize zb(:,1)
	545	zb(1:ngrid,1)=0
	546	end if
	547
	548	DO ig=1,ngrid
	549	z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
	550	zc(ig,nlay)=za(ig,nlay)zq(ig,nlay,iq)z1(ig)
	551	zd(ig,nlay)=zb(ig,nlay)*z1(ig)
	552	ENDDO
	553
	554	DO ilay=nlay-1,2,-1
	555	DO ig=1,ngrid
	556	z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
	557	$ zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
	558	zc(ig,ilay)=(za(ig,ilay)*zq(ig,ilay,iq)+
	559	$ zb(ig,ilay+1)zc(ig,ilay+1))z1(ig)
	560	zd(ig,ilay)=zb(ig,ilay)*z1(ig)
	561	ENDDO
	562	ENDDO
	563
	564	if (water.and.(iq.eq.igcm_h2o_ice)) then
	565	! special case for water ice tracer: do not include
	566	! h2o ice tracer from surface (which is set when handling
	567	! h2o vapour case (see further down).
	568	DO ig=1,ngrid
	569	z1(ig)=1./(za(ig,1)+zb(ig,1)+
	570	$ zb(ig,2)*(1.-zd(ig,2)))
	571	zc(ig,1)=(za(ig,1)*zq(ig,1,iq)+
	572	$ zb(ig,2)zc(ig,2)) z1(ig)
	573	ENDDO
	574	else ! general case
	575	DO ig=1,ngrid
	576	z1(ig)=1./(za(ig,1)+zb(ig,1)+
	577	$ zb(ig,2)*(1.-zd(ig,2)))
	578	zc(ig,1)=(za(ig,1)*zq(ig,1,iq)+
	579	$ zb(ig,2)*zc(ig,2) +
	580	$ (-pdqsdif(ig,iq)) ptimestep) z1(ig) !tracer flux from surface
	581	ENDDO
	582	endif ! of if (water.and.(iq.eq.igcm_h2o_ice))
	583
	584	IF ((water).and.(iq.eq.igcm_h2o_vap)) then
	585	c Calculation for turbulent exchange with the surface (for ice)
	586	DO ig=1,ngrid
	587	zd(ig,1)=zb(ig,1)*z1(ig)
	588	zq1temp(ig)=zc(ig,1)+ zd(ig,1)*qsat(ig)
	589
	590	pdqsdif(ig,igcm_h2o_ice)=rho(ig)dryness(ig)zcdv(ig)
	591	& *(zq1temp(ig)-qsat(ig))
	592	c write(,)'flux vers le sol=',pdqsdif(ig,nq)
	593	END DO
	594
	595	DO ig=1,ngrid
	596	if(.not.watercaptag(ig)) then
	597	if ((-pdqsdif(ig,igcm_h2o_ice)*ptimestep)
	598	& .gt.pqsurf(ig,igcm_h2o_ice)) then
	599	c write(,)'on sublime plus que qsurf!'
	600	pdqsdif(ig,igcm_h2o_ice)=
	601	& -pqsurf(ig,igcm_h2o_ice)/ptimestep
	602	c write(,)'flux vers le sol=',pdqsdif(ig,nq)
	603	z1(ig)=1./(za(ig,1)+ zb(ig,2)*(1.-zd(ig,2)))
	604	zc(ig,1)=(za(ig,1)*zq(ig,1,igcm_h2o_vap)+
	605	$ zb(ig,2)*zc(ig,2) +
	606	$ (-pdqsdif(ig,igcm_h2o_ice)) ptimestep) z1(ig)
	607	zq1temp(ig)=zc(ig,1)
	608	endif
	609	endif ! if (.not.watercaptag(ig))
	610	c Starting upward calculations for water :
	611	zq(ig,1,igcm_h2o_vap)=zq1temp(ig)
	612	ENDDO ! of DO ig=1,ngrid
	613	ELSE
	614	c Starting upward calculations for simple mixing of tracer (dust)
	615	DO ig=1,ngrid
	616	zq(ig,1,iq)=zc(ig,1)
	617	ENDDO
	618	END IF ! of IF ((water).and.(iq.eq.igcm_h2o_vap))
	619
	620	DO ilay=2,nlay
	621	DO ig=1,ngrid
	622	zq(ig,ilay,iq)=zc(ig,ilay)+zd(ig,ilay)*zq(ig,ilay-1,iq)
	623	ENDDO
	624	ENDDO
	625	enddo ! of do iq=1,nq
	626	end if ! of if(tracer)
	627
	628	c-----------------------------------------------------------------------
	629	c 8. calcul final des tendances de la diffusion verticale
	630	c ----------------------------------------------------
	631
	632	DO ilev = 1, nlay
	633	DO ig=1,ngrid
	634	pdudif(ig,ilev)=( zu(ig,ilev)-
	635	$ (pu(ig,ilev)+pdufi(ig,ilev)*ptimestep) )/ptimestep
	636	pdvdif(ig,ilev)=( zv(ig,ilev)-
	637	$ (pv(ig,ilev)+pdvfi(ig,ilev)*ptimestep) )/ptimestep
	638	hh = max(ph(ig,ilev)+pdhfi(ig,ilev)*ptimestep ,
	639	$ zhcond(ig,ilev)) ! modif co2cond
	640	pdhdif(ig,ilev)=( zh(ig,ilev)- hh )/ptimestep
	641	ENDDO
	642	ENDDO
	643
	644
	645	if (tracer) then
	646	DO iq = 1, nq
	647	DO ilev = 1, nlay
	648	DO ig=1,ngrid
	649	pdqdif(ig,ilev,iq)=(zq(ig,ilev,iq)-
	650	$ (pq(ig,ilev,iq) + pdqfi(ig,ilev,iq)*ptimestep))/ptimestep
	651	ENDDO
	652	ENDDO
	653	ENDDO
	654	end if
	655
	656	c ** diagnostique final
	657	c ------------------
	658
	659	IF(lecrit) THEN
	660	PRINT*,'In vdif'
	661	PRINT*,'Ts (t) and Ts (t+st)'
	662	WRITE(*,'(a10,3a15)')
	663	s 'theta(t)','theta(t+dt)','u(t)','u(t+dt)'
	664	PRINT*,ptsrf(ngrid/2+1),ztsrf2(ngrid/2+1)
	665	DO ilev=1,nlay
	666	WRITE(*,'(4f15.7)')
	667	s ph(ngrid/2+1,ilev),zh(ngrid/2+1,ilev),
	668	s pu(ngrid/2+1,ilev),zu(ngrid/2+1,ilev)
	669
	670	ENDDO
	671	ENDIF
	672
[161]	673	if (calltherm .and. outptherm) then
	674	if (ngrid .eq. 1) then
	675	call WRITEDIAGFI(ngrid,'zh','zh inside vdifc',
	676	& 'SI',1,ph(:,:)+pdhfi(:,:)*ptimestep)
	677	call WRITEDIAGFI(ngrid,'zkh','zkh',
	678	& 'SI',1,zkh)
	679	endif
	680	endif
	681
[38]	682	RETURN
	683	END

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

Download in other formats:

Original Format