Context Navigation

vdifc.F @ 541

Last change on this file since 541 was 529, checked in by acolaitis, 13 years ago
Physiq: changed called to pbl_parameter and vdifc. Vdifc: added a flag to deactivate mass-variation scheme for LES computations for two reason: it is costy for nothing and does not work in polar LES.
File size: 26.9 KB

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

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

Download in other formats:

Original Format