Context Navigation

vdifc.F @ 662

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

Thermals: updated computation of wstar to take into account vertical turbulent heat flux from diffusion. Without that, wstar is underestimated. We use the sensible heat flux and tendancies computed in vdifc to compute the vertical turbulent heat flux contribution from the diffusion, and then compute the maximum in the PBL. The resulting total heat flux profile and wstar correspond very well to LES results.

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

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

Download in other formats:

Original Format