Context Navigation

thermcell_dqupdown.F90 @ 171

Last change on this file since 171 was 161, checked in by acolaitis, 13 years ago

================================================
======== IMPLEMENTATION OF THERMALS ============
================================================

Author: A. Colaitis (2011-06-16)

The main goal of this revision is to start including the thermals into the model
for development purposes. Users should not use the thermals yet, as
several major configuration changes still need to be done.

This version includes :

updraft and downdraft parametrizations
velocity in the thermal, including drag
plume height analysis
closure equation
updraft transport of heat, tracers and momentum
downdraft transport of heat

This model should not be used without upcoming developments, namely :

downdraft transport of tracers and momentum
updraft & downdraft transport of q2 (tke)
revision of vdif_kc to compute q2 for non-stratified cases

Thermals could also include in a later revision :

momentum loss during transport (horizontal drag)

Compilation of the thermals has been successfully tested on ifort, gfortran and pgf90

================================================
================================================

M libf/phymars/callkeys.h
M libf/phymars/inifis.F

Added new control flags to call the thermals :

calltherm (false by default) <- to call thermals
outptherm (false by default) <- to output thermal-related diagnostics (for dev purposes)

================================================

M libf/phymars/vdifc.F
^{------> added a temporary output for thermal-related diagnostics}

M libf/phymars/testphys1d.F
^{------> added treatment for a initialization from a profile of neutral gas (ar)}

-> will be transformed in a decaying tracer for thermal diagnostics

M libf/phymars/physiq.F
^{------> added a section to call the thermals}

-> changed the call to convadj
-> added thermal-related outputs for diagnostics

M libf/phymars/convadj.F
^{------> takes now into account the height of thermals to execute convective adjustment}

=> note : convective adjustment needs to be activated when using thermals, in case of a

second instable layer above the thermals

================================================

A libf/phymars/calltherm_interface.F90
^{------> Interface between physiq.F and the thermals}

A libf/phymars/calltherm_mars.F90
^{------> Routine running the sub-timestep of the thermals}

A libf/phymars/thermcell_main_mars.F90
^{------> Main thermals routine specific to Martian physics}

A libf/phymars/thermcell_dqupdown.F90
^{------> Thermals subroutine computing transport of quantities by updrafts and downdrafts}

A libf/phymars/thermcell.F90
^{------> Module including parameters from the Earth to Mars importation. Will disappear in future dev}

================================================
================================================

Property svn:executable set to *

File size: 9.5 KB

Rev	Line
[161]	1	subroutine thermcell_dqupdown(ngrid,nlay,ptimestep,fm0,entr0, &
	2	& detr0,masse0,q_therm,dq_therm,ztvd,fm_down,ztv,charvar,lmax)
	3	implicit none
	4
	5	! #include "iniprint.h"
	6	!=======================================================================
	7	!
	8	! Calcul du transport verticale dans la couche limite en presence
	9	! de "thermiques" explicitement representes
	10	! calcul du dq/dt une fois qu'on connait les ascendances
	11	! Version modifiee pour prendre les downdrafts a la place de la
	12	! subsidence compensatoire
	13	!=======================================================================
	14
	15	! ============================ INPUTS ============================
	16
	17	INTEGER, INTENT(IN) :: ngrid,nlay
	18	REAL, INTENT(IN) :: ptimestep
	19	REAL, INTENT(IN) :: fm0(ngrid,nlay+1)
	20	REAL, INTENT(IN) :: entr0(ngrid,nlay),detr0(ngrid,nlay)
	21	REAL, INTENT(IN) :: q_therm(ngrid,nlay)
	22	REAL, INTENT(IN) :: fm_down(ngrid,nlay+1)
	23	REAL, INTENT(IN) :: ztvd(ngrid,nlay)
	24	REAL, INTENT(IN) :: ztv(ngrid,nlay)
	25	CHARACTER (LEN=20), INTENT(IN) :: charvar
	26	REAL, INTENT(IN) :: masse0(ngrid,nlay)
	27	INTEGER, INTENT(IN) :: lmax(ngrid)
	28
	29	! ============================ OUTPUTS ===========================
	30
	31	REAL, INTENT(OUT) :: dq_therm(ngrid,nlay)
	32
	33	! ============================ LOCAL =============================
	34
	35	! REAL detr0(ngrid,nlay)
	36	REAL detrd(ngrid,nlay)
	37	REAL entrd(ngrid,nlay)
	38	REAL fmd(ngrid,nlay+1)
	39	REAL q(ngrid,nlay)
	40	REAL qa(ngrid,nlay)
	41	REAL qd(ngrid,nlay)
	42	INTEGER ig,k
	43	LOGICAL active(ngrid,nlay)
	44	INTEGER lmax_down(ngrid),lmin_down(ngrid)
	45	INTEGER ncorec
	46
	47	! =========== Init ==============================================
	48
	49	entrd(:,:)=0.
	50	detrd(:,:)=0.
	51	qa(:,:)=q_therm(:,:)
	52	q(:,:)=q_therm(:,:)
	53	qd(:,:)=q_therm(:,:)
	54	active(:,:)=.false.
	55
	56	! previous calculation of zdthl_down uses the divergence of fmd
	57	! so it can be negative without problem. Here we include the sign
	58	! of fmd in the equations, so it has to be positive
	59	!
	60	! fmd(:,:)=-fm_down(:,:)
	61	!
	62	!! ========== Entrainment, Detrainement and Mass =================
	63	!
	64	!! ========== DOWNDRAFT TRANSPORT DISABLED FOR NOW ===============
	65	!
	66	! do ig=1,ngrid
	67	! if (ztv(ig,nlay)-ztvd(ig,nlay) .gt. 0.5) then
	68	! print*,"downdraft non nul derniere couche !!! (dqupdown)"
	69	! endif
	70	! detrd(ig,nlay)=0.
	71	! entrd(ig,nlay)=0.
	72	! enddo
	73	!
	74	! do k=nlay-1,1,-1
	75	! do ig=1,ngrid
	76	!
	77	! if (ztv(ig,k)-ztvd(ig,k) .gt. 0.0001) then
	78	! detrd(ig,k)=MAX(0.,(fmd(ig,k+1)*(ztv(ig,k)-ztvd(ig,k+1))) &
	79	! & /(ztv(ig,k)-ztvd(ig,k)) - fmd(ig,k))
	80	! entrd(ig,k)=MAX(0.,(fmd(ig,k+1)*(ztvd(ig,k)-ztvd(ig,k+1))) &
	81	! & /(ztv(ig,k)-ztvd(ig,k)))
	82	!
	83	! endif
	84	! enddo
	85	! enddo
	86	!
	87	!! ======= We have computed entrainment and detrainment from a prescribed
	88	!! mass flux and potential temp profile. Due to the way downdraft are parametrized,
	89	!! this can yield negative entr and detr. We force it to be positive, but in
	90	!! order to conserve tracers, we need to recompute an adequate mass flux
	91	!! and modify interface rates, to preserve consistency.
	92	! lmax_down(:)=1
	93	! lmin_down(:)=1
	94	!
	95	! do k=1,nlay
	96	! do ig=1,ngrid
	97	! if ((entrd(ig,k).gt.0.) .or. (detrd(ig,k).gt.0.)) then
	98	!! if (entrd(ig,k).gt.detrd(ig,k)) then
	99	! lmax_down(ig)=min(k,lmax(ig))
	100	! endif
	101	! enddo
	102	! enddo
	103	! do k=nlay,1,-1
	104	! do ig=1,ngrid
	105	! if ((entrd(ig,k).gt.0.) .or. (detrd(ig,k).gt.0.)) then
	106	!! if (detrd(ig,k).gt.entrd(ig,k)) then
	107	! lmin_down(ig)=k
	108	! endif
	109	! enddo
	110	! enddo
	111	!
	112	!
	113	! fmd(:,:)=0.
	114	!
	115	! do ig=1,ngrid
	116	! if ((lmax_down(ig).gt.1) .and. ((lmax_down(ig)-lmin_down(ig)).gt.1)) then
	117	!! fmd(ig,lmax_down(ig))=0.
	118	!! entrd(ig,lmax_down(ig))=detrd(ig,lmax_down(ig))
	119	!! detrd(ig,lmax_down(ig))=0.
	120	!! print*,lmin_down(ig),lmax_down(ig),lmax(ig)
	121	!
	122	! fmd(ig,lmax_down(ig)+1)=0.
	123	!
	124	! do k=lmax_down(ig),lmin_down(ig)+1,-1
	125	! fmd(ig,k)=fmd(ig,k+1)+entrd(ig,k)-detrd(ig,k)
	126	! enddo
	127	!
	128	! fmd(ig,lmin_down(ig))=0.
	129	! detrd(ig,lmin_down(ig))=fmd(ig,lmin_down(ig)+1)+entrd(ig,lmin_down(ig))
	130	!
	131	! else
	132	! entrd(ig,:)=0.
	133	! detrd(ig,:)=0.
	134	! active(ig,:)=.false.
	135	! endif
	136	!
	137	! enddo
	138	! ncorec=0
	139	! do k=nlay,2,-1
	140	! do ig=1,ngrid
	141	! if (fmd(ig,k).lt.0.) then
	142	!! detrd(ig,k)=max(0.,detrd(ig,k)+fmd(ig,k-1))
	143	!! fmd(ig,k-1)=0.
	144	!! entrd(ig,k-1)=0.
	145	!! detrd(ig,k-1)=0.
	146	!! lmin_down(ig)=k-1
	147	! fmd(ig,k)=fmd(ig,k+1)
	148	! detrd(ig,k)=entrd(ig,k)
	149	! ncorec=ncorec+1
	150	!! fmd(ig,k)=0.
	151	!! detrd(ig,k)=entrd(ig,k)+fmd(ig,k+1)
	152	!
	153	! endif
	154	! enddo
	155	! enddo
	156	!
	157	! if (ncorec .ne. 0) then
	158	! print*, 'corrections for negative downward mass flux :',ncorec
	159	! endif
	160	! print*, lmin_down(:),lmax_down(:)
	161	!
	162	! do k=2,nlay
	163	! do ig=1,ngrid
	164	! active(ig,k)=(k.ge.lmin_down(ig)).and.(k.le.lmax_down(ig)) &
	165	! & .and.(((fmd(ig,k)+detrd(ig,k))ptimestep).gt.1.e-6masse0(ig,k))
	166	! enddo
	167	! enddo
	168	!
	169	!
	170	! do ig=1,ngrid
	171	! do k=lmin_down(ig),lmax_down(ig)
	172	! if(.not.active(ig,k)) then
	173	! active(ig,:)=.false.
	174	! endif
	175	! enddo
	176	! enddo
	177	!
	178	! if(charvar .eq. 'tke') then
	179	! active(:,:)=.false.
	180	! endif
	181	!
	182	!! do ig=1,ngrid
	183	!! active(ig,lmax_down(ig))=(((fmd(ig,lmax_down(ig))+detrd(ig,lmax_down(ig)))* &
	184	!! & ptimestep).gt.1.e-8*masse0(ig,lmax_down(ig)))
	185	!! enddo
	186	!!
	187	!! do ig=1,ngrid
	188	!! if (lmax_down(ig).gt.1) then
	189	!! do k=lmax_down(ig)-1,lmin_down(ig),-1
	190	!! active(ig,k)=(((fmd(ig,k)+detrd(ig,k))* &
	191	!! & ptimestep).gt.1.e-8*masse0(ig,k)) &
	192	!! & .and. active(ig,k+1)
	193	!! enddo
	194	!! else
	195	!! active(ig,:)=.false.
	196	!! endif
	197	!! enddo
	198	!!
	199	!! ========== qa : q in updraft ==================================
	200	!
	201	do k=2,nlay
	202	do ig=1,ngrid
	203	if ((fm0(ig,k+1)+detr0(ig,k))*ptimestep.gt. &
	204	& 1.e-5*masse0(ig,k)) then
	205	qa(ig,k)=(fm0(ig,k)qa(ig,k-1)+entr0(ig,k)q(ig,k)) &
	206	& /(fm0(ig,k+1)+detr0(ig,k))
	207
	208	if ((qa(ig,k).lt.0.) .and. (charvar .ne. 'momentum')) then
	209	print*,'qa<0!!!',charvar,ig,k,fm0(ig,k),qa(ig,k-1),entr0(ig,k),q(ig,k),fm0(ig,k+1),detr0(ig,k)
	210	print*,'---------> Cancelling qa'
	211	qa(ig,k)=q(ig,k)
	212	endif
	213	else
	214	qa(ig,k)=q(ig,k)
	215	endif
	216	enddo
	217	enddo
	218
	219	! ========== qd : q in downdraft =================================
	220	!
	221	! do k=nlay-1,1,-1
	222	! do ig=1,ngrid
	223	! if (active(ig,k)) then
	224	! qd(ig,k)=(fmd(ig,k+1)qd(ig,k+1)+entrd(ig,k)q(ig,k)) &
	225	! & /(fmd(ig,k)+detrd(ig,k))
	226	!
	227	! if ((qd(ig,k).lt.0.) .and. (charvar .ne. 'momentum')) then
	228	! print*,'qd<0!!!',charvar,ig,k,fmd(ig,k),qd(ig,k),entrd(ig,k),q(ig,k),fmd(ig,k+1),detrd(ig,k),lmin_down(ig),lmax_down(ig)
	229	! print*, '---------> cancelling qd, no downdraft for this gridpoint'
	230	! qd(ig,k)=q(ig,k)
	231	! active(ig,:)=.false.
	232	! endif
	233	! else
	234	! qd(ig,k)=q(ig,k)
	235	! endif
	236	!! print*,'active,k,entr,detr,q,qd (down) :',active(ig,k),k,entrd(ig,k),detrd(ig,k),q(ig,k),qd(ig,k)
	237	! enddo
	238	! enddo
	239	!
	240	!
	241	! ====== dq ======================================================
	242
	243	do ig=1,ngrid
	244	if(active(ig,1)) then
	245
	246	dq_therm(ig,1)=(detr0(ig,1)qa(ig,1)+detrd(ig,1)qd(ig,1) &
	247	& +fm0(ig,2)*q(ig,2) &
	248	& -entr0(ig,1)q(ig,1)-entrd(ig,1)q(ig,1) &
	249	& -fmd(ig,2)*q(ig,1)) &
	250	& *ptimestep/masse0(ig,1)
	251
	252	else
	253	dq_therm(ig,1)=(detr0(ig,1)qa(ig,1)+fm0(ig,2)q(ig,2) &
	254	& -entr0(ig,1)*q(ig,1)) &
	255	& *ptimestep/masse0(ig,1)
	256
	257	endif
	258	enddo
	259
	260	do k=2,nlay-1
	261	do ig=1, ngrid
	262
	263	if(active(ig,k)) then
	264
	265	dq_therm(ig,k)=(detr0(ig,k)qa(ig,k)+detrd(ig,k)qd(ig,k) &
	266	& +fm0(ig,k+1)q(ig,k+1)+fmd(ig,k)q(ig,k-1) &
	267	& -entr0(ig,k)q(ig,k)-entrd(ig,k)q(ig,k) &
	268	& -fm0(ig,k)q(ig,k)-fmd(ig,k+1)q(ig,k)) &
	269	& *ptimestep/masse0(ig,k)
	270
	271
	272	else
	273	dq_therm(ig,k)=(detr0(ig,k)qa(ig,k)+fm0(ig,k+1)q(ig,k+1) &
	274	& -entr0(ig,k)q(ig,k)-fm0(ig,k)q(ig,k)) &
	275	& *ptimestep/masse0(ig,k)
	276
	277
	278	endif
	279
	280	enddo
	281	enddo
	282
	283	do ig=1, ngrid
	284
	285	if(active(ig,nlay)) then
	286
	287	dq_therm(ig,nlay)=(detr0(ig,nlay)qa(ig,nlay)+detrd(ig,nlay)qd(ig,nlay) &
	288	& +fmd(ig,nlay)*q(ig,nlay-1) &
	289	& -entr0(ig,nlay)q(ig,nlay)-entrd(ig,nlay)q(ig,nlay) &
	290	& -fm0(ig,nlay)*q(ig,nlay)) &
	291	& *ptimestep/masse0(ig,nlay)
	292
	293	else
	294	dq_therm(ig,nlay)=(detr0(ig,nlay)*qa(ig,nlay) &
	295	& -entr0(ig,nlay)q(ig,nlay)-fm0(ig,nlay)q(ig,nlay)) &
	296	& *ptimestep/masse0(ig,nlay)
	297	endif
	298
	299	enddo
	300	return
	301	end

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

Download in other formats:

Original Format