Context Navigation

vdifc_mod.F @ 2252

Last change on this file since 2252 was 2218, checked in by emillour, 5 years ago
Mars GCM: Change "latentheat" flag to a more descriptive "latentheat_surfwater" and set its default value to .true. EM
File size: 31.5 KB

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

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

Download in other formats:

Original Format