Context Navigation

vdifc.F @ 706

Last change on this file since 706 was 669, checked in by tnavarro, 13 years ago

Ice effective radius for output is weighted by ice total surface instead of ice total mass in physiq.F

New distribution for permanent ice reservoirs in surfini.F

icelocationmode = 1 allows for realistic ice distribution, no matter what resolution.
icelocationmode = 2 is predefined 32x24 or 64x48 and should be USED BY DEFAULT. It currently overestimates ice.
icelocationmode = 3 is good for quick logical definitions (paleoclimates,stability studies,etc...)

Possibility to change perihelion date in solar longitude as well as in sol in tabfi.F for newstarts
Improved latent heat release when ground ice sublimates: now a implicit scheme for ground temperature

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

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

Download in other formats:

Original Format