Context Navigation

vdifc.F @ 544

Last change on this file since 544 was 544, checked in by acolaitis, 13 years ago

27/02/12 == AC

Continuation of thermals setting, comparisons with mesoscale results on Case C
Added possibility to call gcm (or 1d) with constant prescribed sensible heat flux, in the spirit of direct comparisons with LES

... This is directly comparable to the variable tke_heat_flux in namelist.input
... Requires the use of yamada4.F from terrestrial GCM (mixes more, seem more numerically stable)
... Usually requires high timesteps (>1000) to avoids crashes. Best approach is to compare
height of first model level z1 and teta_1 between LES and 1D, and increase the timesteps until results
between the two models are comparable (might require a slitghly different tke_heat_flux between the two models
due to difference in vertical diffusion schemes and subgrid effects)
... Syntax for use is to add "tke_heat_flux = ..." in callphys.def

Corrected some stuff with tke transport in thermals (which is anyway deactivated but one day maybe...)

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

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

Download in other formats:

Original Format