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

Line
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