Context Navigation

vdifc.F @ 498

Last change on this file since 498 was 496, checked in by acolaitis, 14 years ago

M 489 libf/phymars/thermcell_main_mars.F90
^{--------------------> New parameters for thermals, following latest version of the relevant article}

M 489 libf/phymars/physiq.F
^{--------------------> Minor modifications}

D 489 libf/phymars/surflayer_interpol.F
A 0 libf/phymars/pbl_parameters.F
^{--------------------> Replaced surflayer_interpol.F by a cleaner and richer version called pbl_parameters.F}

This routine is the base of what will later be implemented in the MCD and as a tool to
compute surface properties from output fields (so that it will ultimately disappear from libf
and might become an utility)

M 489 libf/phymars/vdif_cd.F
^{--------------------> Minor modification to the Richardson computation}

M 489 libf/phymars/vdifc.F
^{--------------------> Minor modification to the aerodynamic conductances computation, replacing wmax by hfmax, which}

is a better proxy for convective activity. Might be replaced by w_star later.

File size: 26.4 KB

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

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

Download in other formats:

Original Format