Context Navigation

mass_redistribution.F90 @ 1529

Last change on this file since 1529 was 1529, checked in by emillour, 9 years ago

Generic GCM: Towards a cleaner separation between physics and dynamics

Got rid of references to "dimensions.h" from physics packages: use nbp_lon (=iim), nbp_lat (=jjp1) and nbp_lev from module mod_grid_phy_lmdz (in phy_common) instead.
Removed module "comhdiff_mod.F90", as it is only used by module surf_heat_transp_mod.F90, moved module variables there.
Addedin "surf_heat_transp_mod" local versions of some arrays and routines (from dyn3d) required to compute gradient, divergence, etc. on the global dynamics grid. As before, the slab ocean only works in serial.

File size: 13.4 KB

Line
1	SUBROUTINE mass_redistribution(ngrid,nlayer,nq,ptimestep, &
2	rnat,pcapcal,pplay,pplev,pt,ptsrf,pq,pqs, &
3	pu,pv,pdt,pdtsrf,pdq,pdu,pdv,pdmassmr, &
4	pdtmr,pdtsrfmr,pdpsrfmr,pdumr,pdvmr,pdqmr,pdqsmr)
5
6	USE watercommon_h, Only: Tsat_water,RLVTT
7	use surfdat_h
8	use radcommon_h, only: glat
9	USE comgeomfi_h
10	USE tracer_h
11	USE planete_mod, only: bp
12	use comcstfi_mod, only: g
13	USE callkeys_mod, ONLY: water
14
15	IMPLICIT NONE
16	!=======================================================================
17	! subject:
18	! --------
19	! Mass and momentum fluxes through sigma levels as the surface pressure is modified are also taken into account
20	!
21	! author: Jeremy Leconte 2012 (from F.Forget 1998)
22	! ------
23	!
24	! input:
25	! -----
26	! ngrid nombre de points de verticales
27	! (toutes les boucles de la physique sont au
28	! moins vectorisees sur ngrid)
29	! nlayer nombre de couches
30	! pplay(ngrid,nlayer) Pressure levels
31	! pplev(ngrid,nlayer+1) Pressure levels
32	! nq Number of tracers
33	!
34	! pt(ngrid,nlayer) temperature (en K)
35	! pq(ngrid,nlayer,nq) tracer specific concentration (kg/kg of air)
36	! pu,pv (ngrid,nlayer) wind velocity (m/s)
37	!
38	!
39	! pdX physical tendency of X before mass redistribution
40	!
41	! pdmassmr air Mass added to the atmosphere in each layer (kg/m2/s)
42	!
43	! output:
44	! -------
45	!
46	! pdXmr(ngrid) physical tendency of X after mass redistribution
47	!
48	!
49	!
50	!=======================================================================
51	!
52	! 0. Declarations :
53	! ------------------
54
55	!-----------------------------------------------------------------------
56	! Arguments :
57	! ---------
58	INTEGER ngrid, nlayer, nq
59	REAL ptimestep
60	REAL pcapcal(ngrid)
61	INTEGER rnat(ngrid)
62	REAL pplay(ngrid,nlayer),pplev(ngrid,nlayer+1)
63	REAL pt(ngrid,nlayer),pdt(ngrid,nlayer)
64	REAL ptsrf(ngrid),pdtsrf(ngrid)
65	REAL pdtmr(ngrid,nlayer)
66	REAL pu(ngrid,nlayer) , pv(ngrid,nlayer)
67	REAL pdu(ngrid,nlayer) , pdv(ngrid,nlayer)
68	REAL pdmassmr(ngrid,nlayer)
69	REAL pdumr(ngrid,nlayer) , pdvmr(ngrid,nlayer)
70	REAL pq(ngrid,nlayer,nq),pdq(ngrid,nlayer,nq)
71	REAL pqs(ngrid,nq)
72	REAL pdqmr(ngrid,nlayer,nq),pdqsmr(ngrid,nq)
73	REAL pdpsrfmr(ngrid),pdtsrfmr(ngrid)
74	!
75	! Local variables :
76	! -----------------
77
78	! Boiling/sublimation
79	REAL Tsat(ngrid),zmassboil(ngrid)
80
81	! vertical reorganization of sigma levels
82	REAL zzu(nlayer),zzv(nlayer)
83	REAL zzq(nlayer,nq),zzt(nlayer)
84	! Dummy variables
85	INTEGER n,l,ig,iq
86	REAL zdtsig(ngrid,nlayer)
87	REAL zmass(ngrid,nlayer),zzmass(nlayer),w(nlayer+1)
88	REAL zdmass_sum(ngrid,nlayer+1)
89	REAL zmflux(nlayer+1)
90	REAL zq1(nlayer)
91	REAL ztsrf(ngrid)
92	REAL ztm(nlayer+1)
93	REAL zum(nlayer+1) , zvm(nlayer+1)
94	REAL zqm(nlayer+1,nq),zqm1(nlayer+1)
95	REAL sigma(nlayer+1)
96
97	! local saved variables
98	LOGICAL, SAVE :: firstcall=.true.
99	!$OMP THREADPRIVATE(firstcall)
100
101	!----------------------------------------------------------------------
102
103	! Initialisation
104	! --------------
105	!
106	IF (firstcall) THEN
107	firstcall=.false.
108	ENDIF
109	!
110	!======================================================================
111	! Calcul of h2o condensation
112	! ============================================================
113	!
114	! Used variable :
115	! pdmassmr : air Mass added to the atmosphere in each layer per unit time (kg/m2/s)
116	! zdmass_sum(ngrid,l) : total air mass added to the atm above layer l per unit time (kg/m2/s)
117	!
118	!
119	! Surface tracer Tendencies set to 0
120	! -------------------------------------
121	pdqsmr(1:ngrid,1:nq)=0.
122
123	ztsrf(1:ngrid) = ptsrf(1:ngrid) + pdtsrf(1:ngrid)*ptimestep
124
125
126	DO ig=1,ngrid
127	zdmass_sum(ig,nlayer+1)=0.
128	DO l = nlayer, 1, -1
129	zmass(ig,l) = (pplev(ig,l)-pplev(ig,l+1))/glat(ig)
130	zdmass_sum(ig,l)= zdmass_sum(ig,l+1)+pdmassmr(ig,l)
131	END DO
132	END DO
133
134
135	if (water) then
136	do ig=1,ngrid
137	call Tsat_water(pplev(ig,1)+zdmass_sum(ig,1)gptimestep,Tsat(ig))
138	enddo
139	call writediagfi(ngrid,'Tsat','saturation temperature at surface','',2,Tsat)
140
141	do ig=1,ngrid
142	if (ztsrf(ig).gt.Tsat(ig)) then
143	zmassboil(ig)=(ptsrf(ig)-Tsat(ig))*pcapcal(ig)/RLVTT/ptimestep
144	if ((zmassboil(ig)*ptimestep.gt.pqs(ig,igcm_h2o_vap)).and.(rnat(ig).eq.1)) then
145	zmassboil(ig)=pqs(ig,igcm_h2o_vap)/ptimestep
146	endif
147	zmassboil(ig)=zmassboil(ig)*0.0 !momentary, should be 1. JL12
148	pdqsmr(ig,igcm_h2o_vap)=-zmassboil(ig)
149	pdtsrfmr(ig)=-zmassboil(ig)*RLVTT/pcapcal(ig)
150	ztsrf(ig)=ptsrf(ig)+pdtsrfmr(ig)*ptimestep
151	else
152	zmassboil(ig)=0.
153	pdtsrfmr(ig)=0.
154	endif
155	enddo
156	endif
157
158	! *************************
159	! UPDATE SURFACE
160	! *************************
161	! Changing pressure at the surface:
162	! """"""""""""""""""""""""""""""""""""
163
164	pdpsrfmr(1:ngrid) = (zdmass_sum(1:ngrid,1)+zmassboil(1:ngrid))*g
165
166	do ig = 1, ngrid
167	IF(ABS(pdpsrfmr(ig)*ptimestep).GT.pplev(ig,1)) THEN
168	PRINT*,'STOP in condens in mass_redistribution'
169	PRINT*,'condensing more than total mass'
170	PRINT*,'Grid point ',ig
171	PRINT*,'Ps = ',pplev(ig,1)
172	PRINT,'d Ps = ',pdpsrfmr(ig)ptimestep
173	STOP
174	ENDIF
175	enddo ! of DO ig=1,ngrid
176
177
178	! ***************************************************************
179	! Correction to account for redistribution between sigma or hybrid
180	! layers when changing surface pressure
181	! zzx quantities have dimension (nlayer) to speed up calculation
182	! *************************************************************
183
184	DO ig=1,ngrid
185	zzt(1:nlayer) = pt(ig,1:nlayer) + pdt(ig,1:nlayer) * ptimestep
186	zzu(1:nlayer) = pu(ig,1:nlayer) + pdu(ig,1:nlayer) * ptimestep
187	zzv(1:nlayer) = pv(ig,1:nlayer) + pdv(ig,1:nlayer) * ptimestep
188	zzq(1:nlayer,1:nq)=pq(ig,1:nlayer,1:nq)+pdq(ig,1:nlayer,1:nq)*ptimestep ! must add the water that has fallen???
189
190	! Mass fluxes of air through the sigma levels (kg.m-2.s-1) (>0 when up)
191	! """"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
192	zmflux(1) = zmassboil(ig)
193	sigma(1)=1
194	DO l=1,nlayer
195	! Ehouarn: shouldn't we rather compute sigma levels than use bp()?
196	! sigma(l+1)=pplev(ig,l+1)/pplev(ig,1)
197	! zmflux(l+1) = zmflux(l) + pdmassmr(ig,l) - &
198	! (sigma(l)-sigma(l+1))*(zdmass_sum(ig,1)+zmflux(1))
199	! if (abs(zmflux(l+1)).lt.1E-13.OR.sigma(l+1).eq.0.) zmflux(l+1)=0.
200	! Ehouarn: but for now leave things as before
201	zmflux(l+1) = zmflux(l) + pdmassmr(ig,l) - (bp(l)-bp(l+1))*(zdmass_sum(ig,1)+zmflux(1))
202	! zmflux set to 0 if very low to avoid: top layer is disappearing in v1ld
203	if (abs(zmflux(l+1)).lt.1E-13.OR.bp(l+1).eq.0.) zmflux(l+1)=0.
204	END DO
205
206	! Mass of each layer
207	! ------------------
208	zzmass(1:nlayer)=zmass(ig,1:nlayer)(1.+pdpsrfmr(ig)ptimestep/pplev(ig,1))
209
210
211	! Corresponding fluxes for T,U,V,Q
212	! """"""""""""""""""""""""""""""""
213
214	! averaging operator for TRANSPORT
215	! """"""""""""""""""""""""""""""""
216	! Value transfert at the surface interface when condensation
217	! sublimation:
218	ztm(1) = ztsrf(ig)
219	zum(1) = 0.
220	zvm(1) = 0.
221	zqm(1,1:nq)=0. ! most tracer do not condense !
222	if (water) zqm(1,igcm_h2o_vap)=1. ! flux is 100% h2o at surface
223
224	! Van Leer scheme:
225	w(1:nlayer+1)=-zmflux(1:nlayer+1)*ptimestep
226	call vl1d(nlayer,zzt,2.,zzmass,w,ztm)
227	call vl1d(nlayer,zzu,2.,zzmass,w,zum)
228	call vl1d(nlayer,zzv,2.,zzmass,w,zvm)
229	do iq=1,nq
230	zq1(1:nlayer)=zzq(1:nlayer,iq)
231	zqm1(1)=zqm(1,iq)
232	! print*,iq
233	! print*,zq1
234	call vl1d(nlayer,zq1,2.,zzmass,w,zqm1)
235	do l=2,nlayer
236	zzq(l,iq)=zq1(l)
237	zqm(l,iq)=zqm1(l)
238	enddo
239	enddo
240
241	! Surface condensation affects low winds
242	if (zmflux(1).lt.0) then
243	zum(1)= zzu(1) * (w(1)/zzmass(1))
244	zvm(1)= zzv(1) * (w(1)/zzmass(1))
245	if (w(1).gt.zzmass(1)) then ! ensure numerical stability
246	zum(1)= (zzu(1)-zum(2))*zzmass(1)/w(1) + zum(2)
247	zvm(1)= (zzv(1)-zvm(2))*zzmass(1)/w(1) + zvm(2)
248	end if
249	end if
250
251	ztm(nlayer+1)= zzt(nlayer) ! should not be used, but...
252	zum(nlayer+1)= zzu(nlayer) ! should not be used, but...
253	zvm(nlayer+1)= zzv(nlayer) ! should not be used, but...
254	zqm(nlayer+1,1:nq)= zzq(nlayer,1:nq)
255
256	! Tendencies on T, U, V, Q
257	! """"""""""""""""""""""""
258	DO l=1,nlayer
259
260	! Tendencies on T
261	pdtmr(ig,l) = (1/zzmass(l)) * &
262	(zmflux(l)(ztm(l) - zzt(l))-zmflux(l+1)(ztm(l+1)-zzt(l)))
263	!JL12 the last term in Newcondens has been set to zero because we are only dealing with redistribution here
264
265	! Tendencies on U
266	pdumr(ig,l) = (1/zzmass(l)) ( zmflux(l)(zum(l) - zzu(l)) - zmflux(l+1)*(zum(l+1) - zzu(l)) )
267
268	! Tendencies on V
269	pdvmr(ig,l) = (1/zzmass(l)) ( zmflux(l)(zvm(l) - zzv(l)) - zmflux(l+1)*(zvm(l+1) - zzv(l)) )
270
271	END DO
272
273	! Tendencies on Q
274	do iq=1,nq
275	DO l=1,nlayer
276	pdqmr(ig,l,iq)= (1/zzmass(l)) * &
277	(zmflux(l)(zqm(l,iq)-zzq(l,iq))- zmflux(l+1)(zqm(l+1,iq)-zzq(l,iq)) - pdmassmr(ig,l)*zzq(l,iq))
278	END DO
279	enddo
280
281	END DO ! loop on ig
282
283	CONTAINS
284
285	! *****************************************************************
286	SUBROUTINE vl1d(llm,q,pente_max,zzmass,w,qm)
287	!
288	!
289	! Operateur de moyenne inter-couche pour calcul de transport type
290	! Van-Leer " pseudo amont " dans la verticale
291	! q,w sont des arguments d'entree pour le s-pg ....
292	! masse : masse de la couche Dp/g
293	! w : masse d'atm ``transferee'' a chaque pas de temps (kg.m-2)
294	! pente_max = 2 conseillee
295	!
296	!
297	! --------------------------------------------------------------------
298
299	IMPLICIT NONE
300
301	! Arguments:
302	! ----------
303	integer,intent(in) :: llm
304	real zzmass(llm),pente_max
305	REAL q(llm),qm(llm+1)
306	REAL w(llm+1)
307	!
308	! Local
309	! ---------
310	!
311	INTEGER l
312	!
313	real dzq(llm),dzqw(llm),adzqw(llm),dzqmax
314	real sigw, Mtot, MQtot
315	integer m
316	! integer ismax,ismin
317
318
319	! On oriente tout dans le sens de la pression
320	! W > 0 WHEN DOWN !!!!!!!!!!!!!
321
322	do l=2,llm
323	dzqw(l)=q(l-1)-q(l)
324	adzqw(l)=abs(dzqw(l))
325	enddo
326
327	do l=2,llm-1
328	if(dzqw(l)*dzqw(l+1).gt.0.) then
329	dzq(l)=0.5*(dzqw(l)+dzqw(l+1))
330	else
331	dzq(l)=0.
332	endif
333	dzqmax=pente_max*min(adzqw(l),adzqw(l+1))
334	dzq(l)=sign(min(abs(dzq(l)),dzqmax),dzq(l))
335	enddo
336
337	dzq(1)=0.
338	dzq(llm)=0.
339
340	do l = 1,llm-1
341
342	! Regular scheme (transfered mass < layer mass)
343	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
344	if(w(l+1).gt.0. .and. w(l+1).le.zzmass(l+1)) then
345	sigw=w(l+1)/zzmass(l+1)
346	qm(l+1)=(q(l+1)+0.5(1.-sigw)dzq(l+1))
347	else if(w(l+1).le.0. .and. -w(l+1).le.zzmass(l)) then
348	sigw=w(l+1)/zzmass(l)
349	qm(l+1)=(q(l)-0.5(1.+sigw)dzq(l))
350
351	! Extended scheme (transfered mass > layer mass)
352	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
353	else if(w(l+1).gt.0.) then
354	m=l+1
355	Mtot = zzmass(m)
356	MQtot = zzmass(m)*q(m)
357	do while ((m.lt.llm).and.(w(l+1).gt.(Mtot+zzmass(m+1))))
358	m=m+1
359	Mtot = Mtot + zzmass(m)
360	MQtot = MQtot + zzmass(m)*q(m)
361	end do
362	if (m.lt.llm) then
363	sigw=(w(l+1)-Mtot)/zzmass(m+1)
364	qm(l+1)= (1/w(l+1))(MQtot + (w(l+1)-Mtot)(q(m+1)+0.5(1.-sigw)dzq(m+1)) )
365	else
366	! w(l+1) = Mtot
367	! qm(l+1) = Mqtot / Mtot
368	write(,) 'top layer is disappearing !',l,Mtot,w(l+1),qm(l+1)
369	print*,zzmass
370	print*,w
371	print*,q
372	print*,qm
373	stop
374	end if
375	else ! if(w(l+1).lt.0)
376	m = l-1
377	Mtot = zzmass(m+1)
378	MQtot = zzmass(m+1)*q(m+1)
379	if (m.gt.0) then ! because some compilers will have problems
380	! evaluating zzmass(0)
381	do while ((m.gt.0).and.(-w(l+1).gt.(Mtot+zzmass(m))))
382	m=m-1
383	Mtot = Mtot + zzmass(m+1)
384	MQtot = MQtot + zzmass(m+1)*q(m+1)
385	if (m.eq.0) exit
386	end do
387	endif
388	if (m.gt.0) then
389	sigw=(w(l+1)+Mtot)/zzmass(m)
390	qm(l+1)= (-1/w(l+1))(MQtot + (-w(l+1)-Mtot)(q(m)-0.5(1.+sigw)dzq(m)) )
391	else
392	qm(l+1)= (-1/w(l+1))(MQtot + (-w(l+1)-Mtot)qm(1))
393	end if
394	end if
395	enddo
396
397	! boundary conditions (not used in newcondens !!)
398	! qm(llm+1)=0.
399	! if(w(1).gt.0.) then
400	! qm(1)=q(1)
401	! else
402	! qm(1)=0.
403	! end if
404
405	END SUBROUTINE vl1d
406
407	END SUBROUTINE mass_redistribution

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

Download in other formats:

Original Format