Context Navigation

vdifc.F @ 513

Last change on this file since 513 was 501, checked in by acolaitis, 14 years ago

More robust gustiness speed computation. Sorry for changing it so often, this time it should be pretty final... This formulation make the 1D model match LES result very well on high vertical resolution. The offset in mixed layer temperature observed on normal resolution is definitely due to vertical resolution (mixing does not occur high enough in the PBL) and not to surface layer underestimating sensible heat flux.

File size: 26.4 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,
[499]	7	$ pdqdif,pdqsdif,wstar,zcdv_true,zcdh_true,hfmax)
[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)
[496]	46	REAL pu(ngrid,nlay),pv(ngrid,nlay),ph(ngrid,nlay),pt(ngrid,nlay)
[38]	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
[499]	62	REAL wstar(ngridmx), hfmax(ngridmx), zi(ngridmx)
[256]	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
[473]	74	REAL ust(ngridmx),tst(ngridmx)
	75
[38]	76	INTEGER ilev,ig,ilay,nlev
	77
	78	REAL z4st,zdplanck(ngridmx)
	79	REAL zkv(ngridmx,nlayermx+1),zkh(ngridmx,nlayermx+1)
[496]	80	REAL zkq(ngridmx,nlayermx+1)
[38]	81	REAL zcdv(ngridmx),zcdh(ngridmx)
[268]	82	REAL zcdv_true(ngridmx),zcdh_true(ngridmx) ! drag coeff are used by the LES to recompute u* and hfx
[38]	83	REAL zu(ngridmx,nlayermx),zv(ngridmx,nlayermx)
	84	REAL zh(ngridmx,nlayermx)
	85	REAL ztsrf2(ngridmx)
	86	REAL z1(ngridmx),z2(ngridmx)
	87	REAL za(ngridmx,nlayermx),zb(ngridmx,nlayermx)
	88	REAL zb0(ngridmx,nlayermx)
	89	REAL zc(ngridmx,nlayermx),zd(ngridmx,nlayermx)
	90	REAL zcst1
[284]	91	REAL zu2(ngridmx)
[38]	92
	93	EXTERNAL SSUM,SCOPY
	94	REAL SSUM
	95	LOGICAL firstcall
	96	SAVE firstcall
	97
	98	c variable added for CO2 condensation:
	99	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	100	REAL hh , zhcond(ngridmx,nlayermx)
	101	REAL latcond,tcond1mb
	102	REAL acond,bcond
	103	SAVE acond,bcond
	104	DATA latcond,tcond1mb/5.9e5,136.27/
	105
	106	c Tracers :
	107	c ~~~~~~~
	108	INTEGER iq
	109	REAL zq(ngridmx,nlayermx,nqmx)
	110	REAL zq1temp(ngridmx)
	111	REAL rho(ngridmx) ! near surface air density
	112	REAL qsat(ngridmx)
	113	DATA firstcall/.true./
	114
	115	REAL kmixmin
	116
[473]	117	c Mass-variation scheme :
	118	c ~~~~~~~
	119
	120	INTEGER j,l
	121	REAL zcondicea(ngrid,nlayermx)
	122	REAL zt(ngridmx,nlayermx),ztcond(ngridmx,nlayermx+1)
	123	REAL betam(ngridmx,nlayermx),dmice(ngridmx,nlayermx)
	124	REAL pdtc(ngrid,nlayermx)
	125	REAL zhs(ngridmx,nlayermx)
[496]	126	REAL ccond
	127	SAVE ccond
[473]	128
	129	c Theta_m formulation for mass-variation scheme :
	130	c ~~~~~~~
	131
	132	INTEGER ico2,llnt(ngridmx)
	133	SAVE ico2
	134	REAL m_co2, m_noco2, A , B
	135	SAVE A, B, m_co2, m_noco2
	136	REAL vmr_co2(ngridmx,nlayermx)
	137	REAL qco2,mmean
	138
[38]	139	c ** un petit test de coherence
	140	c --------------------------
	141
	142	IF (firstcall) THEN
	143	IF(ngrid.NE.ngridmx) THEN
	144	PRINT*,'STOP dans vdifc'
	145	PRINT*,'probleme de dimensions :'
	146	PRINT*,'ngrid =',ngrid
	147	PRINT*,'ngridmx =',ngridmx
	148	STOP
	149	ENDIF
	150	c To compute: Tcond= 1./(bcond-acondlog(.0095p)) (p in pascal)
	151	bcond=1./tcond1mb
	152	acond=r/latcond
[473]	153	ccond=cpp/(g*latcond)
[38]	154	PRINT*,'In vdifc: Tcond(P=1mb)=',tcond1mb,' Lcond=',latcond
[473]	155	PRINT*,' acond,bcond,ccond',acond,bcond,ccond
[38]	156
[473]	157
	158	ico2=0
	159
	160	if (tracer) then
	161	c Prepare Special treatment if one of the tracer is CO2 gas
	162	do iq=1,nqmx
	163	if (noms(iq).eq."co2") then
	164	ico2=iq
	165	m_co2 = 44.01E-3 ! CO2 molecular mass (kg/mol)
	166	m_noco2 = 33.37E-3 ! Non condensible mol mass (kg/mol)
	167	c Compute A and B coefficient use to compute
	168	c mean molecular mass Mair defined by
	169	c 1/Mair = q(ico2)/m_co2 + (1-q(ico2))/m_noco2
	170	c 1/Mair = A*q(ico2) + B
	171	A =(1/m_co2 - 1/m_noco2)
	172	B=1/m_noco2
	173	endif
	174	enddo
	175	end if
	176
[38]	177	firstcall=.false.
	178	ENDIF
	179
	180
	181
	182
	183
	184	c-----------------------------------------------------------------------
	185	c 1. initialisation
	186	c -----------------
	187
	188	nlev=nlay+1
	189
	190	c ** calcul de rhodz et dtrho/dz=dtrho*2 g/dp
	191	c avec rho=p/RT=p/ (R Theta) (p/ps)**kappa
	192	c ----------------------------------------
	193
	194	DO ilay=1,nlay
	195	DO ig=1,ngrid
	196	za(ig,ilay)=(pplev(ig,ilay)-pplev(ig,ilay+1))/g
[473]	197	! Mass variation scheme:
	198	betam(ig,ilay)=-za(ig,ilay)latcond/(cppppopsk(ig,ilay))
[38]	199	ENDDO
	200	ENDDO
	201
	202	zcst1=4.gptimestep/(r*r)
	203	DO ilev=2,nlev-1
	204	DO ig=1,ngrid
	205	zb0(ig,ilev)=pplev(ig,ilev)*
	206	s (pplev(ig,1)/pplev(ig,ilev))**rcp /
	207	s (ph(ig,ilev-1)+ph(ig,ilev))
	208	zb0(ig,ilev)=zcst1zb0(ig,ilev)zb0(ig,ilev)/
	209	s (pplay(ig,ilev-1)-pplay(ig,ilev))
	210	ENDDO
	211	ENDDO
	212	DO ig=1,ngrid
[473]	213	zb0(ig,1)=ptimesteppplev(ig,1)/(rptsrf(ig))
[38]	214	ENDDO
	215
	216	c ** diagnostique pour l'initialisation
	217	c ----------------------------------
	218
	219	IF(lecrit) THEN
	220	ig=ngrid/2+1
	221	PRINT*,'Pression (mbar) ,altitude (km),u,v,theta, rho dz'
	222	DO ilay=1,nlay
	223	WRITE(*,'(6f11.5)')
	224	s .01pplay(ig,ilay),.001pzlay(ig,ilay),
	225	s pu(ig,ilay),pv(ig,ilay),ph(ig,ilay),za(ig,ilay)
	226	ENDDO
	227	PRINT*,'Pression (mbar) ,altitude (km),zb'
	228	DO ilev=1,nlay
	229	WRITE(*,'(3f15.7)')
	230	s .01pplev(ig,ilev),.001pzlev(ig,ilev),
	231	s zb0(ig,ilev)
	232	ENDDO
	233	ENDIF
	234
[473]	235	c -----------------------------------
[38]	236	c Potential Condensation temperature:
	237	c -----------------------------------
	238
[473]	239	c Compute CO2 Volume mixing ratio
	240	c -------------------------------
	241	if (callcond.and.(ico2.ne.0)) then
	242	DO ilev=1,nlay
	243	DO ig=1,ngrid
	244	qco2=pq(ig,ilev,ico2)+pdqfi(ig,ilev,ico2)*ptimestep
	245	c Mean air molecular mass = 1/(q(ico2)/m_co2 + (1-q(ico2))/m_noco2)
	246	mmean=1/(A*qco2 +B)
	247	vmr_co2(ig,ilev) = qco2*mmean/m_co2
	248	ENDDO
	249	ENDDO
	250	else
	251	DO ilev=1,nlay
	252	DO ig=1,ngrid
	253	vmr_co2(ig,ilev)=0.95
	254	ENDDO
	255	ENDDO
	256	end if
[38]	257
[473]	258	c forecast of atmospheric temperature zt and frost temperature ztcond
	259	c --------------------------------------------------------------------
[38]	260
[473]	261	DO ilev=1,nlay
	262	DO ig=1,ngrid
	263	ztcond(ig,ilev)=
	264	& 1./(bcond-acondlog(.01vmr_co2(ig,ilev)*pplay(ig,ilev)))
	265	if (pplay(ig,ilev).lt.1e-4) ztcond(ig,ilev)=0.0 !mars Monica
	266	ENDDO
	267	ENDDO
	268
	269	ztcond(:,nlay+1)=ztcond(:,nlay)
	270
	271	if (callcond) then
	272	DO ilev=1,nlay
	273	DO ig=1,ngrid
	274	! zhcond(ig,ilev) =
	275	! & (1./(bcond-acondlog(.0095pplay(ig,ilev))))/ppopsk(ig,ilev)
	276	zhcond(ig,ilev) = ztcond(ig,ilev)/ppopsk(ig,ilev)
	277	END DO
	278	END DO
	279	else
	280	call zerophys(ngrid*nlay,zhcond)
	281	end if
	282
	283
[38]	284	c-----------------------------------------------------------------------
	285	c 2. ajout des tendances physiques
	286	c -----------------------------
	287
	288	DO ilev=1,nlay
	289	DO ig=1,ngrid
	290	zu(ig,ilev)=pu(ig,ilev)+pdufi(ig,ilev)*ptimestep
	291	zv(ig,ilev)=pv(ig,ilev)+pdvfi(ig,ilev)*ptimestep
	292	zh(ig,ilev)=ph(ig,ilev)+pdhfi(ig,ilev)*ptimestep
[473]	293	! zh(ig,ilev)=max(zh(ig,ilev),zhcond(ig,ilev))
[38]	294	ENDDO
	295	ENDDO
	296	if(tracer) then
	297	DO iq =1, nq
	298	DO ilev=1,nlay
	299	DO ig=1,ngrid
	300	zq(ig,ilev,iq)=pq(ig,ilev,iq)+pdqfi(ig,ilev,iq)*ptimestep
	301	ENDDO
	302	ENDDO
	303	ENDDO
	304	end if
	305
	306	c-----------------------------------------------------------------------
	307	c 3. schema de turbulence
	308	c --------------------
	309
	310	c ** source d'energie cinetique turbulente a la surface
	311	c (condition aux limites du schema de diffusion turbulente
	312	c dans la couche limite
	313	c ---------------------
	314
[499]	315	CALL vdif_cd(ngrid,nlay,pz0,g,pzlay,pu,pv,wstar,ptsrf,ph
[38]	316	& ,zcdv_true,zcdh_true)
[256]	317
[290]	318	zu2(:)=pu(:,1)pu(:,1)+pv(:,1)pv(:,1)
[256]	319
[291]	320	IF (callrichsl) THEN
[499]	321	zcdv(:)=zcdv_true(:)*sqrt(zu2(:)+
[501]	322	& (log(1.+0.7wstar(:) + 2.3wstar(:)2))2)
	323	! & (14.75sqrt(hfmax(:)) - 10.hfmax(:) + 1.7hfmax(:)2)*2) ! subgrid gustiness added by quadratic interpolation with a coeff beta
[499]	324	zcdh(:)=zcdh_true(:)*sqrt(zu2(:)+
[501]	325	& (log(1.+0.7wstar(:) + 2.3wstar(:)2))2)
	326	! & (8.hfmax(:) + 0.hfmax(:)*2 + 0.hfmax(:)4)2) ! LES comparisons. This parameter is linked to the thermals settings)
[496]	327
[499]	328	ust(:)=sqrt(zcdv_true(:))*sqrt(zu2(:)+
[501]	329	& (log(1.+0.7wstar(:) + 2.3wstar(:)2))2)
	330	! & (14.75sqrt(hfmax(:)) - 10.hfmax(:) + 1.7hfmax(:)2)*2)
[496]	331	tst(:)=(ph(:,1)-ptsrf(:))*zcdh(:)/ust(:)
[473]	332	! ust(:)=sqrt(zcdv_true(:))*sqrt(zu2(:))
	333	! tst(:)=(ph(:,1)-ptsrf(:))*zcdh_true(:)/sqrt(zcdv_true(:))
[284]	334	ELSE
	335	zcdv(:)=zcdv_true(:)*sqrt(zu2(:)) ! 1 / bulk aerodynamic momentum conductance
	336	zcdh(:)=zcdh_true(:)*sqrt(zu2(:)) ! 1 / bulk aerodynamic heat conductance
[501]	337	ust(:)=sqrt(zcdv_true(:))*sqrt(zu2(:))
	338	tst(:)=(ph(:,1)-ptsrf(:))*zcdh_true(:)/sqrt(zcdv_true(:))
[290]	339	ENDIF
[38]	340
[290]	341
[256]	342	! Some usefull diagnostics for the new surface layer parametrization :
[501]	343
[268]	344	! call WRITEDIAGFI(ngridmx,'pcdv',
[256]	345	! & 'momentum drag','no units',
	346	! & 2,zcdv_true)
[268]	347	! call WRITEDIAGFI(ngridmx,'pcdh',
[256]	348	! & 'heat drag','no units',
	349	! & 2,zcdh_true)
[496]	350	! call WRITEDIAGFI(ngridmx,'ust',
[473]	351	! & 'friction velocity','m/s',2,ust)
[496]	352	! call WRITEDIAGFI(ngridmx,'tst',
[473]	353	! & 'friction temperature','K',2,tst)
[268]	354	! call WRITEDIAGFI(ngridmx,'rm-1',
	355	! & 'aerodyn momentum conductance','m/s',
[256]	356	! & 2,zcdv)
[268]	357	! call WRITEDIAGFI(ngridmx,'rh-1',
	358	! & 'aerodyn heat conductance','m/s',
[256]	359	! & 2,zcdh)
	360
[38]	361	c ** schema de diffusion turbulente dans la couche limite
	362	c ----------------------------------------------------
[496]	363	!
	364	CALL vdif_kc(ptimestep,g,pzlev,pzlay
[38]	365	& ,pu,pv,ph,zcdv_true
[325]	366	& ,pq2,zkv,zkh,zq)
[496]	367	!
	368	! pt(:,:)=ph(:,:)*ppopsk(:,:)
	369	! CALL yamada4(ngrid,nlay,ptimestep,g,r,pplev,pt
	370	! s ,pzlev,pzlay,pu,pv,ph,zcdv_true,pq2,zkv,zkh,zkq,ust
	371	! s ,8)
[38]	372
	373	if ((doubleq).and.(ngrid.eq.1)) then
	374	kmixmin = 80. !80.! minimum eddy mix coeff in 1D
	375	do ilev=1,nlay
	376	do ig=1,ngrid
	377	zkh(ig,ilev) = max(kmixmin,zkh(ig,ilev))
	378	zkv(ig,ilev) = max(kmixmin,zkv(ig,ilev))
	379	end do
	380	end do
	381	end if
	382
	383	c ** diagnostique pour le schema de turbulence
	384	c -----------------------------------------
	385
	386	IF(lecrit) THEN
	387	PRINT*
	388	PRINT*,'Diagnostic for the vertical turbulent mixing'
	389	PRINT*,'Cd for momentum and potential temperature'
	390
	391	PRINT*,zcdv(ngrid/2+1),zcdh(ngrid/2+1)
	392	PRINT*,'Mixing coefficient for momentum and pot.temp.'
	393	DO ilev=1,nlay
	394	PRINT*,zkv(ngrid/2+1,ilev),zkh(ngrid/2+1,ilev)
	395	ENDDO
	396	ENDIF
	397
	398
	399
	400
	401	c-----------------------------------------------------------------------
	402	c 4. inversion pour l'implicite sur u
	403	c --------------------------------
	404
	405	c ** l'equation est
	406	c u(t+1) = u(t) + dt * {(du/dt)phys}(t) + dt * {(du/dt)difv}(t+1)
	407	c avec
	408	c /zu/ = u(t) + dt * {(du/dt)phys}(t) (voir paragraphe 2.)
	409	c et
	410	c dt * {(du/dt)difv}(t+1) = dt * {(d/dz)[ Ku (du/dz) ]}(t+1)
	411	c donc les entrees sont /zcdv/ pour la condition a la limite sol
	412	c et /zkv/ = Ku
	413
	414	CALL multipl((nlay-1)*ngrid,zkv(1,2),zb0(1,2),zb(1,2))
	415	CALL multipl(ngrid,zcdv,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)zu(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)*zu(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	DO ig=1,ngrid
	434	zu(ig,1)=zc(ig,1)
	435	ENDDO
	436	DO ilay=2,nlay
	437	DO ig=1,ngrid
	438	zu(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*zu(ig,ilay-1)
	439	ENDDO
	440	ENDDO
	441
	442
	443
	444
	445
	446	c-----------------------------------------------------------------------
	447	c 5. inversion pour l'implicite sur v
	448	c --------------------------------
	449
	450	c ** l'equation est
	451	c v(t+1) = v(t) + dt * {(dv/dt)phys}(t) + dt * {(dv/dt)difv}(t+1)
	452	c avec
	453	c /zv/ = v(t) + dt * {(dv/dt)phys}(t) (voir paragraphe 2.)
	454	c et
	455	c dt * {(dv/dt)difv}(t+1) = dt * {(d/dz)[ Kv (dv/dz) ]}(t+1)
	456	c donc les entrees sont /zcdv/ pour la condition a la limite sol
	457	c et /zkv/ = Kv
	458
	459	DO ig=1,ngrid
	460	z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
	461	zc(ig,nlay)=za(ig,nlay)zv(ig,nlay)z1(ig)
	462	zd(ig,nlay)=zb(ig,nlay)*z1(ig)
	463	ENDDO
	464
	465	DO ilay=nlay-1,1,-1
	466	DO ig=1,ngrid
	467	z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
	468	$ zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
	469	zc(ig,ilay)=(za(ig,ilay)*zv(ig,ilay)+
	470	$ zb(ig,ilay+1)zc(ig,ilay+1))z1(ig)
	471	zd(ig,ilay)=zb(ig,ilay)*z1(ig)
	472	ENDDO
	473	ENDDO
	474
	475	DO ig=1,ngrid
	476	zv(ig,1)=zc(ig,1)
	477	ENDDO
	478	DO ilay=2,nlay
	479	DO ig=1,ngrid
	480	zv(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*zv(ig,ilay-1)
	481	ENDDO
	482	ENDDO
	483
	484
	485
	486
	487
	488	c-----------------------------------------------------------------------
	489	c 6. inversion pour l'implicite sur h sans oublier le couplage
	490	c avec le sol (conduction)
	491	c ------------------------
	492
	493	c ** l'equation est
	494	c h(t+1) = h(t) + dt * {(dh/dt)phys}(t) + dt * {(dh/dt)difv}(t+1)
	495	c avec
	496	c /zh/ = h(t) + dt * {(dh/dt)phys}(t) (voir paragraphe 2.)
	497	c et
	498	c dt * {(dh/dt)difv}(t+1) = dt * {(d/dz)[ Kh (dh/dz) ]}(t+1)
	499	c donc les entrees sont /zcdh/ pour la condition de raccord au sol
	500	c et /zkh/ = Kh
	501	c -------------
	502
[473]	503	c Mass variation scheme:
[38]	504	CALL multipl((nlay-1)*ngrid,zkh(1,2),zb0(1,2),zb(1,2))
	505	CALL multipl(ngrid,zcdh,zb0,zb)
	506
[473]	507	c on initialise dm c
	508
	509	pdtc(:,:)=0.
	510	zt(:,:)=0.
	511	dmice(:,:)=0.
[38]	512
	513	c ** calcul de (d Planck / dT) a la temperature d'interface
	514	c ------------------------------------------------------
	515
	516	z4st=4.5.67e-8ptimestep
	517	DO ig=1,ngrid
	518	zdplanck(ig)=z4stpemis(ig)ptsrf(ig)ptsrf(ig)ptsrf(ig)
	519	ENDDO
	520
[473]	521
	522	! calcul de zc et zd pour la couche top en prenant en compte le terme
	523	! de variation de masse (on fait une boucle pour que ça converge)
	524
	525	! Identification des points de grilles qui ont besoin de la correction
	526
	527	llnt(:)=1
	528
	529	DO ig=1,ngrid
	530	DO l=1,nlay
	531	if(zh(ig,l) .lt. zhcond(ig,l)) then
	532	llnt(ig)=300
	533	! 200 and 100 do not go beyond month 9 with normal dissipation
	534	goto 5
	535	endif
	536	ENDDO
	537	5 continue
	538	ENDDO
	539
	540	DO ig=1,ngrid
	541
	542	! Initialization of z1 and zd, which do not depend on dmice
	543
	544	z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
	545	zd(ig,nlay)=zb(ig,nlay)*z1(ig)
	546
	547	DO ilay=nlay-1,1,-1
	548	z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
	549	$ zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
	550	zd(ig,ilay)=zb(ig,ilay)*z1(ig)
	551	ENDDO
	552
	553	! Convergence loop
	554
	555	DO j=1,llnt(ig)
	556
	557	z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
	558	zc(ig,nlay)=za(ig,nlay)*zh(ig,nlay)
	559	& -betam(ig,nlay)*dmice(ig,nlay)
	560	zc(ig,nlay)=zc(ig,nlay)*z1(ig)
	561	! zd(ig,nlay)=zb(ig,nlay)*z1(ig)
	562
	563	! calcul de zc et zd pour les couches du haut vers le bas
	564
	565	DO ilay=nlay-1,1,-1
	566	z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
	567	$ zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
	568	zc(ig,ilay)=(za(ig,ilay)*zh(ig,ilay)+
	569	$ zb(ig,ilay+1)*zc(ig,ilay+1)-
	570	$ betam(ig,ilay)dmice(ig,ilay))z1(ig)
	571	! zd(ig,ilay)=zb(ig,ilay)*z1(ig)
	572	ENDDO
	573
[38]	574	c ** calcul de la temperature_d'interface et de sa tendance.
	575	c on ecrit que la somme des flux est nulle a l'interface
	576	c a t + \delta t,
	577	c c'est a dire le flux radiatif a {t + \delta t}
	578	c + le flux turbulent a {t + \delta t}
	579	c qui s'ecrit K (T1-Tsurf) avec T1 = d1 Tsurf + c1
	580	c (notation K dt = /cpp*b/)
	581	c + le flux dans le sol a t
	582	c + l'evolution du flux dans le sol lorsque la temperature d'interface
	583	c passe de sa valeur a t a sa valeur a {t + \delta t}.
	584	c ----------------------------------------------------
	585
	586	z1(ig)=pcapcal(ig)ptsrf(ig)+cppzb(ig,1)*zc(ig,1)
	587	s +zdplanck(ig)ptsrf(ig)+ pfluxsrf(ig)ptimestep
	588	z2(ig)= pcapcal(ig)+cppzb(ig,1)(1.-zd(ig,1))+zdplanck(ig)
	589	ztsrf2(ig)=z1(ig)/z2(ig)
[473]	590	! pdtsrf(ig)=(ztsrf2(ig)-ptsrf(ig))/ptimestep !incremented outside loop
[38]	591
[473]	592	zhs(ig,1)=zc(ig,1)+zd(ig,1)*ztsrf2(ig)
[38]	593
	594	c ** et a partir de la temperature au sol on remonte
	595	c -----------------------------------------------
	596
[473]	597	DO ilay=2,nlay
	598	zhs(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*zhs(ig,ilay-1)
	599	ENDDO
	600	DO ilay=1,nlay
	601	zt(ig,ilay)=zhs(ig,ilay)*ppopsk(ig,ilay)
	602	ENDDO
	603
	604	c Condensation/sublimation in the atmosphere
	605	c ------------------------------------------
	606	c (computation of zcondicea and dmice)
	607
	608	zcondicea(ig,:)=0.
	609	pdtc(ig,:)=0.
	610
	611	DO l=nlay , 1, -1
	612	IF(zt(ig,l).LT.ztcond(ig,l)) THEN
	613	pdtc(ig,l)=(ztcond(ig,l) - zt(ig,l))/ptimestep
	614	zcondicea(ig,l)=(pplev(ig,l)-pplev(ig,l+1))
	615	& ccondpdtc(ig,l)
	616	dmice(ig,l)= dmice(ig,l) + zcondicea(ig,l)*ptimestep
	617	END IF
[38]	618	ENDDO
	619
[473]	620	ENDDO !of Do j=1,XXX
[38]	621
[473]	622	ENDDO !of Do ig=1,nlay
	623
	624	pdtsrf(:)=(ztsrf2(:)-ptsrf(:))/ptimestep
	625
	626
[38]	627	c-----------------------------------------------------------------------
	628	c TRACERS
	629	c -------
	630
	631	if(tracer) then
	632
	633	c Using the wind modified by friction for lifting and sublimation
	634	c ----------------------------------------------------------------
	635	DO ig=1,ngrid
[284]	636	zu2(ig)=zu(ig,1)zu(ig,1)+zv(ig,1)zv(ig,1)
	637	zcdv(ig)=zcdv_true(ig)*sqrt(zu2(ig))
	638	zcdh(ig)=zcdh_true(ig)*sqrt(zu2(ig))
[38]	639	ENDDO
	640
	641	c Calcul du flux vertical au bas de la premiere couche (dust) :
	642	c -----------------------------------------------------------
	643	do ig=1,ngridmx
	644	rho(ig) = zb0(ig,1) /ptimestep
	645	c zb(ig,1) = 0.
	646	end do
	647	c Dust lifting:
	648	if (lifting) then
[310]	649	#ifndef MESOSCALE
[38]	650	if (doubleq.AND.submicron) then
	651	do ig=1,ngrid
	652	c if(co2ice(ig).lt.1) then
	653	pdqsdif(ig,igcm_dust_mass) =
	654	& -alpha_lift(igcm_dust_mass)
	655	pdqsdif(ig,igcm_dust_number) =
	656	& -alpha_lift(igcm_dust_number)
	657	pdqsdif(ig,igcm_dust_submicron) =
	658	& -alpha_lift(igcm_dust_submicron)
	659	c end if
	660	end do
	661	else if (doubleq) then
	662	do ig=1,ngrid
[77]	663	!!! soulevement constant
[38]	664	pdqsdif(ig,igcm_dust_mass) =
	665	& -alpha_lift(igcm_dust_mass)
	666	pdqsdif(ig,igcm_dust_number) =
	667	& -alpha_lift(igcm_dust_number)
	668	end do
	669	else if (submicron) then
	670	do ig=1,ngrid
	671	pdqsdif(ig,igcm_dust_submicron) =
	672	& -alpha_lift(igcm_dust_submicron)
	673	end do
	674	else
	675	call dustlift(ngrid,nlay,nq,rho,zcdh_true,zcdh,co2ice,
	676	& pdqsdif)
	677	endif !doubleq.AND.submicron
[310]	678	#else
	679	call dustlift(ngrid,nlay,nq,rho,zcdh_true,zcdh,co2ice,
	680	& pdqsdif)
	681	#endif
[38]	682	else
	683	pdqsdif(1:ngrid,1:nq) = 0.
	684	end if
	685
	686	c OU calcul de la valeur de q a la surface (water) :
	687	c ----------------------------------------
	688	if (water) then
	689	call watersat(ngridmx,ptsrf,pplev(1,1),qsat)
	690	end if
	691
	692	c Inversion pour l'implicite sur q
	693	c --------------------------------
	694	do iq=1,nq
	695	CALL multipl((nlay-1)*ngrid,zkh(1,2),zb0(1,2),zb(1,2))
	696
	697	if ((water).and.(iq.eq.igcm_h2o_vap)) then
	698	c This line is required to account for turbulent transport
	699	c from surface (e.g. ice) to mid-layer of atmosphere:
	700	CALL multipl(ngrid,zcdv,zb0,zb(1,1))
	701	CALL multipl(ngrid,dryness,zb(1,1),zb(1,1))
	702	else ! (re)-initialize zb(:,1)
	703	zb(1:ngrid,1)=0
	704	end if
	705
	706	DO ig=1,ngrid
	707	z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
	708	zc(ig,nlay)=za(ig,nlay)zq(ig,nlay,iq)z1(ig)
	709	zd(ig,nlay)=zb(ig,nlay)*z1(ig)
	710	ENDDO
	711
	712	DO ilay=nlay-1,2,-1
	713	DO ig=1,ngrid
	714	z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
	715	$ zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
	716	zc(ig,ilay)=(za(ig,ilay)*zq(ig,ilay,iq)+
	717	$ zb(ig,ilay+1)zc(ig,ilay+1))z1(ig)
	718	zd(ig,ilay)=zb(ig,ilay)*z1(ig)
	719	ENDDO
	720	ENDDO
	721
	722	if (water.and.(iq.eq.igcm_h2o_ice)) then
	723	! special case for water ice tracer: do not include
	724	! h2o ice tracer from surface (which is set when handling
[473]	725	! h2o vapour case (see further down).
[38]	726	DO ig=1,ngrid
	727	z1(ig)=1./(za(ig,1)+zb(ig,1)+
	728	$ zb(ig,2)*(1.-zd(ig,2)))
	729	zc(ig,1)=(za(ig,1)*zq(ig,1,iq)+
	730	$ zb(ig,2)zc(ig,2)) z1(ig)
	731	ENDDO
	732	else ! general case
	733	DO ig=1,ngrid
	734	z1(ig)=1./(za(ig,1)+zb(ig,1)+
	735	$ zb(ig,2)*(1.-zd(ig,2)))
	736	zc(ig,1)=(za(ig,1)*zq(ig,1,iq)+
	737	$ zb(ig,2)*zc(ig,2) +
	738	$ (-pdqsdif(ig,iq)) ptimestep) z1(ig) !tracer flux from surface
	739	ENDDO
	740	endif ! of if (water.and.(iq.eq.igcm_h2o_ice))
	741
	742	IF ((water).and.(iq.eq.igcm_h2o_vap)) then
	743	c Calculation for turbulent exchange with the surface (for ice)
	744	DO ig=1,ngrid
	745	zd(ig,1)=zb(ig,1)*z1(ig)
	746	zq1temp(ig)=zc(ig,1)+ zd(ig,1)*qsat(ig)
	747
	748	pdqsdif(ig,igcm_h2o_ice)=rho(ig)dryness(ig)zcdv(ig)
	749	& *(zq1temp(ig)-qsat(ig))
	750	c write(,)'flux vers le sol=',pdqsdif(ig,nq)
	751	END DO
	752
	753	DO ig=1,ngrid
	754	if(.not.watercaptag(ig)) then
	755	if ((-pdqsdif(ig,igcm_h2o_ice)*ptimestep)
	756	& .gt.pqsurf(ig,igcm_h2o_ice)) then
	757	c write(,)'on sublime plus que qsurf!'
	758	pdqsdif(ig,igcm_h2o_ice)=
	759	& -pqsurf(ig,igcm_h2o_ice)/ptimestep
	760	c write(,)'flux vers le sol=',pdqsdif(ig,nq)
	761	z1(ig)=1./(za(ig,1)+ zb(ig,2)*(1.-zd(ig,2)))
	762	zc(ig,1)=(za(ig,1)*zq(ig,1,igcm_h2o_vap)+
	763	$ zb(ig,2)*zc(ig,2) +
	764	$ (-pdqsdif(ig,igcm_h2o_ice)) ptimestep) z1(ig)
	765	zq1temp(ig)=zc(ig,1)
	766	endif
	767	endif ! if (.not.watercaptag(ig))
	768	c Starting upward calculations for water :
	769	zq(ig,1,igcm_h2o_vap)=zq1temp(ig)
	770	ENDDO ! of DO ig=1,ngrid
	771	ELSE
	772	c Starting upward calculations for simple mixing of tracer (dust)
	773	DO ig=1,ngrid
	774	zq(ig,1,iq)=zc(ig,1)
	775	ENDDO
	776	END IF ! of IF ((water).and.(iq.eq.igcm_h2o_vap))
	777
	778	DO ilay=2,nlay
	779	DO ig=1,ngrid
	780	zq(ig,ilay,iq)=zc(ig,ilay)+zd(ig,ilay)*zq(ig,ilay-1,iq)
	781	ENDDO
	782	ENDDO
	783	enddo ! of do iq=1,nq
	784	end if ! of if(tracer)
	785
	786	c-----------------------------------------------------------------------
	787	c 8. calcul final des tendances de la diffusion verticale
	788	c ----------------------------------------------------
	789
	790	DO ilev = 1, nlay
	791	DO ig=1,ngrid
	792	pdudif(ig,ilev)=( zu(ig,ilev)-
	793	$ (pu(ig,ilev)+pdufi(ig,ilev)*ptimestep) )/ptimestep
	794	pdvdif(ig,ilev)=( zv(ig,ilev)-
	795	$ (pv(ig,ilev)+pdvfi(ig,ilev)*ptimestep) )/ptimestep
[473]	796	hh = ph(ig,ilev)+pdhfi(ig,ilev)*ptimestep
	797	$ + (latcond*dmice(ig,ilev)/cpp)/ppopsk(ig,ilev)
	798	pdhdif(ig,ilev)=( zhs(ig,ilev)- hh )/ptimestep
[38]	799	ENDDO
	800	ENDDO
	801
	802	if (tracer) then
	803	DO iq = 1, nq
	804	DO ilev = 1, nlay
	805	DO ig=1,ngrid
	806	pdqdif(ig,ilev,iq)=(zq(ig,ilev,iq)-
	807	$ (pq(ig,ilev,iq) + pdqfi(ig,ilev,iq)*ptimestep))/ptimestep
	808	ENDDO
	809	ENDDO
	810	ENDDO
	811	end if
	812
	813	c ** diagnostique final
	814	c ------------------
	815
	816	IF(lecrit) THEN
	817	PRINT*,'In vdif'
	818	PRINT*,'Ts (t) and Ts (t+st)'
	819	WRITE(*,'(a10,3a15)')
	820	s 'theta(t)','theta(t+dt)','u(t)','u(t+dt)'
	821	PRINT*,ptsrf(ngrid/2+1),ztsrf2(ngrid/2+1)
	822	DO ilev=1,nlay
	823	WRITE(*,'(4f15.7)')
[473]	824	s ph(ngrid/2+1,ilev),zhs(ngrid/2+1,ilev),
[38]	825	s pu(ngrid/2+1,ilev),zu(ngrid/2+1,ilev)
	826
	827	ENDDO
	828	ENDIF
	829
	830	RETURN
	831	END

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

Download in other formats:

Original Format