Context Navigation

vdifc.F @ 1242

Last change on this file since 1242 was 1242, checked in by aslmd, 11 years ago
LMDZ.MARS. MESOSCALE. Corrected problems affecting LES runs with new physics (sensibFlux calculated in vdifc is just used as an input to thermals, it does not have the right sign and apparently is not W/m2). Further reduction or factorisation of mesoscale includes. Also made ustar and tstar available in turb_mod from the calculations in vdifc.
File size: 28.3 KB

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

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

Download in other formats:

Original Format