Context Navigation

vdifc.F @ 1917

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

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

Download in other formats:

Original Format