Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_small_step_em.F @ 1

Last change on this file since 1 was 1, checked in by lfita, 10 years ago

-- --- Opening of the WRF+LMDZ coupling repository --- -- -

WRF: version v3.3
LMDZ: version v1818

More details in:

File size: 61.7 KB

Line
1	!WRF:MODEL_LAYER:DYNAMICS
2	!
3
4	! SMALL_STEP code for the geometric height coordinate model
5	!
6	!---------------------------------------------------------------------------
7
8	MODULE module_small_step_em
9
10	USE module_configure
11	USE module_model_constants
12
13	CONTAINS
14
15	!----------------------------------------------------------------------
16
17	SUBROUTINE small_step_prep( u_1, u_2, v_1, v_2, w_1, w_2, &
18	t_1, t_2, ph_1, ph_2, &
19	mub, mu_1, mu_2, &
20	muu, muus, muv, muvs, &
21	mut, muts, mudf, &
22	u_save, v_save, w_save, &
23	t_save, ph_save, mu_save, &
24	ww, ww_save, &
25	dnw, c2a, pb, p, alt, &
26	msfux, msfuy, msfvx, &
27	msfvx_inv, &
28	msfvy, msftx, msfty, &
29	rdx, rdy, &
30	rk_step, &
31	ids,ide, jds,jde, kds,kde, &
32	ims,ime, jms,jme, kms,kme, &
33	its,ite, jts,jte, kts,kte )
34
35	IMPLICIT NONE ! religion first
36
37	! declarations for the stuff coming in
38
39	INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde
40	INTEGER, INTENT(IN ) :: ims,ime, jms,jme, kms,kme
41	INTEGER, INTENT(IN ) :: its,ite, jts,jte, kts,kte
42
43	INTEGER, INTENT(IN ) :: rk_step
44
45	REAL, DIMENSION(ims:ime, kms:kme, jms:jme),INTENT(INOUT) :: u_1, &
46	v_1, &
47	w_1, &
48	t_1, &
49	ph_1
50
51	REAL, DIMENSION(ims:ime, kms:kme, jms:jme),INTENT( OUT) :: u_save, &
52	v_save, &
53	w_save, &
54	t_save, &
55	ph_save
56
57	REAL, DIMENSION(ims:ime, kms:kme, jms:jme),INTENT(INOUT) :: u_2, &
58	v_2, &
59	w_2, &
60	t_2, &
61	ph_2
62
63	REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT( OUT) :: c2a, &
64	ww_save
65
66	REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN ) :: pb, &
67	p, &
68	alt, &
69	ww
70
71	REAL, DIMENSION(ims:ime, jms:jme) , INTENT(INOUT) :: mu_1,mu_2
72
73	REAL, DIMENSION(ims:ime, jms:jme) , INTENT(IN ) :: mub, &
74	muu, &
75	muv, &
76	mut, &
77	msfux,&
78	msfuy,&
79	msfvx,&
80	msfvx_inv,&
81	msfvy,&
82	msftx,&
83	msfty
84
85	REAL, DIMENSION(ims:ime, jms:jme) , INTENT( OUT) :: muus, &
86	muvs, &
87	muts, &
88	mudf
89
90	REAL, DIMENSION(ims:ime, jms:jme) , INTENT( OUT) :: mu_save
91
92	REAL, DIMENSION(kms:kme, jms:jme) , INTENT(IN ) :: dnw
93
94	REAL, INTENT(IN) :: rdx,rdy
95
96	! local variables
97
98	INTEGER :: i, j, k
99	INTEGER :: i_start, i_end, j_start, j_end, k_start, k_end
100	INTEGER :: i_endu, j_endv
101
102	!<DESCRIPTION>
103	!
104	! small_step_prep prepares the prognostic variables for the small timestep.
105	! This includes switching time-levels in the arrays and computing coupled
106	! perturbation variables for the small timestep
107	! (i.e. muu" = mu(t)u(t)-mu()u(); muu" is advanced during the small
108	! timesteps
109	!
110	!</DESCRIPTION>
111
112	i_start = its
113	i_end = min(ite,ide-1)
114	j_start = jts
115	j_end = min(jte,jde-1)
116	k_start = kts
117	k_end = min(kte,kde-1)
118
119	i_endu = ite
120	j_endv = jte
121
122	! if this is the first RK step, reset _1 to _2
123	! (we are replacing the t-dt fields with the time t fields)
124
125	IF ((rk_step == 1) ) THEN
126
127	DO j=j_start, j_end
128	DO i=i_start, i_end
129	mu_1(i,j)=mu_2(i,j)
130	ww_save(i,kde,j) = 0.
131	ww_save(i,1,j) = 0.
132	mudf(i,j) = 0. ! initialize external mode div damp to zero
133	ENDDO
134	ENDDO
135
136	DO j=j_start, j_end
137	DO k=k_start, k_end
138	DO i=i_start, i_endu
139	u_1(i,k,j) = u_2(i,k,j)
140	ENDDO
141	ENDDO
142	ENDDO
143
144	DO j=j_start, j_endv
145	DO k=k_start, k_end
146	DO i=i_start, i_end
147	v_1(i,k,j) = v_2(i,k,j)
148	ENDDO
149	ENDDO
150	ENDDO
151
152	DO j=j_start, j_end
153	DO k=k_start, k_end
154	DO i=i_start, i_end
155	t_1(i,k,j) = t_2(i,k,j)
156	ENDDO
157	ENDDO
158	ENDDO
159
160	DO j=j_start, j_end
161	DO k=k_start, min(kde,kte)
162	DO i=i_start, i_end
163	w_1(i,k,j) = w_2(i,k,j)
164	ph_1(i,k,j) = ph_2(i,k,j)
165	ENDDO
166	ENDDO
167	ENDDO
168
169	DO j=j_start, j_end
170	DO i=i_start, i_end
171	muts(i,j)=mub(i,j)+mu_2(i,j)
172	ENDDO
173	DO i=i_start, i_endu
174	! rk_step==1, WCS fix for tiling
175	! muus(i,j)=0.5*(mub(i,j)+mu_2(i,j)+mub(i-1,j)+mu_2(i-1,j))
176	muus(i,j) = muu(i,j)
177	ENDDO
178	ENDDO
179
180	DO j=j_start, j_endv
181	DO i=i_start, i_end
182	! rk_step==1, WCS fix for tiling
183	! muvs(i,j)=0.5*(mub(i,j)+mu_2(i,j)+mub(i,j-1)+mu_2(i,j-1))
184	muvs(i,j) = muv(i,j)
185	ENDDO
186	ENDDO
187
188	DO j=j_start, j_end
189	DO i=i_start, i_end
190	mu_save(i,j)=mu_2(i,j)
191	mu_2(i,j)=0.
192	! mu_2(i,j)=mu_2(i,j)-mu_2(i,j)
193	ENDDO
194	ENDDO
195
196	ELSE
197
198	DO j=j_start, j_end
199	DO i=i_start, i_end
200	muts(i,j)=mub(i,j)+mu_1(i,j)
201	ENDDO
202	DO i=i_start, i_endu
203	muus(i,j)=0.5*(mub(i,j)+mu_1(i,j)+mub(i-1,j)+mu_1(i-1,j))
204	ENDDO
205	ENDDO
206
207	DO j=j_start, j_endv
208	DO i=i_start, i_end
209	muvs(i,j)=0.5*(mub(i,j)+mu_1(i,j)+mub(i,j-1)+mu_1(i,j-1))
210	ENDDO
211	ENDDO
212
213	DO j=j_start, j_end
214	DO i=i_start, i_end
215	mu_save(i,j)=mu_2(i,j)
216	mu_2(i,j)=mu_1(i,j)-mu_2(i,j)
217	ENDDO
218	ENDDO
219
220
221	END IF
222
223	! set up the small timestep variables
224
225	DO j=j_start, j_end
226	DO i=i_start, i_end
227	ww_save(i,kde,j) = 0.
228	ww_save(i,1,j) = 0.
229	ENDDO
230	ENDDO
231
232	DO j=j_start, j_end
233	DO k=k_start, k_end
234	DO i=i_start, i_end
235	c2a(i,k,j) = cpovcv*(pb(i,k,j)+p(i,k,j))/alt(i,k,j)
236	ENDDO
237	ENDDO
238	ENDDO
239
240	DO j=j_start, j_end
241	DO k=k_start, k_end
242	DO i=i_start, i_endu
243	u_save(i,k,j) = u_2(i,k,j)
244	! u coupled with my
245	u_2(i,k,j) = (muus(i,j)u_1(i,k,j)-muu(i,j)u_2(i,k,j))/msfuy(i,j)
246	ENDDO
247	ENDDO
248	ENDDO
249
250	DO j=j_start, j_endv
251	DO k=k_start, k_end
252	DO i=i_start, i_end
253	v_save(i,k,j) = v_2(i,k,j)
254	! v coupled with mx
255	! v_2(i,k,j) = (muvs(i,j)v_1(i,k,j)-muv(i,j)v_2(i,k,j))/msfvx(i,j)
256	v_2(i,k,j) = (muvs(i,j)v_1(i,k,j)-muv(i,j)v_2(i,k,j))*msfvx_inv(i,j)
257	ENDDO
258	ENDDO
259	ENDDO
260
261	DO j=j_start, j_end
262	DO k=k_start, k_end
263	DO i=i_start, i_end
264	t_save(i,k,j) = t_2(i,k,j)
265	t_2(i,k,j) = muts(i,j)t_1(i,k,j)-mut(i,j)t_2(i,k,j)
266	ENDDO
267	ENDDO
268	ENDDO
269
270	DO j=j_start, j_end
271	! DO k=k_start, min(kde,kte)
272	DO k=k_start, kde
273	DO i=i_start, i_end
274	w_save(i,k,j) = w_2(i,k,j)
275	! w coupled with my
276	w_2(i,k,j) = (muts(i,j)* w_1(i,k,j)-mut(i,j)* w_2(i,k,j))/msfty(i,j)
277	ph_save(i,k,j) = ph_2(i,k,j)
278	ph_2(i,k,j) = ph_1(i,k,j)-ph_2(i,k,j)
279	ENDDO
280	ENDDO
281	ENDDO
282
283	DO j=j_start, j_end
284	! DO k=k_start, min(kde,kte)
285	DO k=k_start, kde
286	DO i=i_start, i_end
287	ww_save(i,k,j) = ww(i,k,j)
288	ENDDO
289	ENDDO
290	ENDDO
291
292	END SUBROUTINE small_step_prep
293
294	!-------------------------------------------------------------------------
295
296
297	SUBROUTINE small_step_finish( u_2, u_1, v_2, v_1, w_2, w_1, &
298	t_2, t_1, ph_2, ph_1, ww, ww1, &
299	mu_2, mu_1, &
300	mut, muts, muu, muus, muv, muvs, &
301	u_save, v_save, w_save, &
302	t_save, ph_save, mu_save, &
303	msfux, msfuy, msfvx, msfvy, &
304	msftx, msfty, &
305	h_diabatic, &
306	number_of_small_timesteps,dts, &
307	rk_step, rk_order, &
308	ids,ide, jds,jde, kds,kde, &
309	ims,ime, jms,jme, kms,kme, &
310	its,ite, jts,jte, kts,kte )
311
312
313	IMPLICIT NONE ! religion first
314
315	! stuff passed in
316
317	INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde
318	INTEGER, INTENT(IN ) :: ims,ime, jms,jme, kms,kme
319	INTEGER, INTENT(IN ) :: its,ite, jts,jte, kts,kte
320	INTEGER, INTENT(IN ) :: number_of_small_timesteps
321	INTEGER, INTENT(IN ) :: rk_step, rk_order
322	REAL, INTENT(IN ) :: dts
323
324	REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN ) :: u_1, &
325	v_1, &
326	w_1, &
327	t_1, &
328	ww1, &
329	ph_1
330
331	REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: u_2, &
332	v_2, &
333	w_2, &
334	t_2, &
335	ww, &
336	ph_2
337
338	REAL, DIMENSION(ims:ime, kms:kme, jms:jme),INTENT(IN ) :: u_save, &
339	v_save, &
340	w_save, &
341	t_save, &
342	ph_save, &
343	h_diabatic
344
345	REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: muus, muvs
346	REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: mu_2, mu_1
347	REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: mut, muts, &
348	muu, muv, mu_save
349	REAL, DIMENSION(ims:ime, jms:jme), INTENT(IN ) :: msfux, msfuy, &
350	msfvx, msfvy, &
351	msftx, msfty
352
353
354	! local stuff
355
356	INTEGER :: i,j,k
357	INTEGER :: i_start, i_end, j_start, j_end, i_endu, j_endv
358
359	!<DESCRIPTION>
360	!
361	! small_step_finish reconstructs the full uncoupled prognostic variables
362	! from the coupled perturbation variables used in the small timesteps.
363	!
364	!</DESCRIPTION>
365
366	i_start = its
367	i_end = min(ite,ide-1)
368	j_start = jts
369	j_end = min(jte,jde-1)
370
371	i_endu = ite
372	j_endv = jte
373
374	! addition of time level t back into variables
375
376	DO j = j_start, j_endv
377	DO k = kds, kde-1
378	DO i = i_start, i_end
379	! v coupled with mx
380	v_2(i,k,j) = (msfvx(i,j)v_2(i,k,j) + v_save(i,k,j)muv(i,j))/muvs(i,j)
381	ENDDO
382	ENDDO
383	ENDDO
384
385	DO j = j_start, j_end
386	DO k = kds, kde-1
387	DO i = i_start, i_endu
388	! u coupled with my
389	u_2(i,k,j) = (msfuy(i,j)u_2(i,k,j) + u_save(i,k,j)muu(i,j))/muus(i,j)
390	ENDDO
391	ENDDO
392	ENDDO
393
394	DO j = j_start, j_end
395	DO k = kds, kde
396	DO i = i_start, i_end
397	! w coupled with my
398	w_2(i,k,j) = (msfty(i,j)w_2(i,k,j) + w_save(i,k,j)mut(i,j))/muts(i,j)
399	ph_2(i,k,j) = ph_2(i,k,j) + ph_save(i,k,j)
400	ww(i,k,j) = ww(i,k,j) + ww1(i,k,j)
401	ENDDO
402	ENDDO
403	ENDDO
404
405	#ifdef REVERT
406	DO j = j_start, j_end
407	DO k = kds, kde-1
408	DO i = i_start, i_end
409	t_2(i,k,j) = (t_2(i,k,j) + t_save(i,k,j)*mut(i,j))/muts(i,j)
410	ENDDO
411	ENDDO
412	ENDDO
413	#else
414	IF ( rk_step < rk_order ) THEN
415	DO j = j_start, j_end
416	DO k = kds, kde-1
417	DO i = i_start, i_end
418	t_2(i,k,j) = (t_2(i,k,j) + t_save(i,k,j)*mut(i,j))/muts(i,j)
419	ENDDO
420	ENDDO
421	ENDDO
422	ELSE
423
424	DO j = j_start, j_end
425	DO k = kds, kde-1
426	DO i = i_start, i_end
427	t_2(i,k,j) = (t_2(i,k,j) - dtsnumber_of_small_timestepsmut(i,j)*h_diabatic(i,k,j) &
428	+ t_save(i,k,j)*mut(i,j))/muts(i,j)
429	ENDDO
430	ENDDO
431	ENDDO
432	ENDIF
433	#endif
434
435	DO j = j_start, j_end
436	DO i = i_start, i_end
437	mu_2(i,j) = mu_2(i,j) + mu_save(i,j)
438	ENDDO
439	ENDDO
440
441	END SUBROUTINE small_step_finish
442
443	!-----------------------------------------------------------------------
444
445	SUBROUTINE calc_p_rho( al, p, ph, &
446	alt, t_2, t_1, c2a, pm1, &
447	mu, muts, znu, t0, &
448	rdnw, dnw, smdiv, &
449	non_hydrostatic, step, &
450	ids, ide, jds, jde, kds, kde, &
451	ims, ime, jms, jme, kms, kme, &
452	its,ite, jts,jte, kts,kte )
453
454	IMPLICIT NONE ! religion first
455
456	! declarations for the stuff coming in
457
458	INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde
459	INTEGER, INTENT(IN ) :: ims,ime, jms,jme, kms,kme
460	INTEGER, INTENT(IN ) :: its,ite, jts,jte, kts,kte
461
462	INTEGER, INTENT(IN ) :: step
463
464	REAL, DIMENSION(ims:ime, kms:kme, jms:jme),INTENT( OUT) :: al, &
465	p
466
467	REAL, DIMENSION(ims:ime, kms:kme, jms:jme),INTENT(IN ) :: alt, &
468	t_2, &
469	t_1, &
470	c2a
471
472	REAL, DIMENSION(ims:ime, kms:kme, jms:jme),INTENT(INOUT) :: ph, pm1
473
474	REAL, DIMENSION(ims:ime, jms:jme) , INTENT(IN ) :: mu, &
475	muts
476
477	REAL, DIMENSION(kms:kme) , INTENT(IN ) :: dnw, &
478	rdnw, &
479	znu
480
481	REAL, INTENT(IN ) :: t0, smdiv
482
483	LOGICAL, INTENT(IN ) :: non_hydrostatic
484
485	! local variables
486
487	INTEGER :: i, j, k
488	INTEGER :: i_start, i_end, j_start, j_end, k_start, k_end
489	REAL :: ptmp
490
491	!<DESCRIPTION>
492	!
493	! For the nonhydrostatic option,
494	! calc_p_rho computes the perturbation inverse density and
495	! perturbation pressure from the hydrostatic relation and the
496	! linearized equation of state, respectively.
497	!
498	! For the hydrostatic option,
499	! calc_p_rho computes the perturbation pressure, perturbation density,
500	! and perturbation geopotential
501	! from the vertical coordinate definition, linearized equation of state
502	! and the hydrostatic relation, respectively.
503	!
504	! forward weighting of the pressure (divergence damping) is also
505	! computed here.
506	!
507	!</DESCRIPTION>
508	#ifdef LMDZ1
509	INTEGER :: im2,km2,jm2
510
511	im2=24
512	jm2=21
513	km2=38
514
515	PRINT *,' calc_p_rho: inside: ',im2,km2,jm2
516	PRINT *,' al: ',al(im2,km2,jm2), ' p: ', p(im2,1,jm2),' ph: ',ph(im2,1,jm2)
517	#endif
518
519
520	i_start = its
521	i_end = min(ite,ide-1)
522	j_start = jts
523	j_end = min(jte,jde-1)
524	k_start = kts
525	k_end = min(kte,kde-1)
526
527	IF (non_hydrostatic) THEN
528	DO j=j_start, j_end
529	DO k=k_start, k_end
530	DO i=i_start, i_end
531
532	! al computation is all dry, so ok with moisture
533
534	al(i,k,j)=-1./muts(i,j)(alt(i,k,j)mu(i,j) &
535	+rdnw(k)*(ph(i,k+1,j)-ph(i,k,j)))
536
537	! this is temporally linearized p, no moisture correction needed
538
539	p(i,k,j)=c2a(i,k,j)(alt(i,k,j)(t_2(i,k,j)-mu(i,j)*t_1(i,k,j)) &
540	/(muts(i,j)*(t0+t_1(i,k,j)))-al (i,k,j))
541
542	ENDDO
543	ENDDO
544	ENDDO
545	#ifdef LMDZ1
546	PRINT *,' calc_p_rho: after nonhydrostatic'
547	PRINT *,' al: ',al(im2,km2,jm2), ' p: ', p(im2,1,jm2),' ph: ',ph(im2,1,jm2)
548	PRINT *,' muts: ',muts(im2,jm2), 'alt: ', alt(im2,1,jm2),' mu: ',mu(im2,jm2), &
549	'rdnw: ', rdnw(1),' ph 1: ',ph(im2,2,jm2)
550	PRINT *,' c2a: ',c2a(im2,1,jm2), 't_2: ', t_2(im2,1,jm2),' t_1: ',t_1(im2,1,jm2), &
551	't0: ', t0,' pm1: ',pm1(im2,1,jm2)
552	#endif
553
554	ELSE ! hydrostatic calculation
555
556	DO j=j_start, j_end
557	DO k=k_start, k_end
558	DO i=i_start, i_end
559	p(i,k,j)=mu(i,j)*znu(k)
560	al(i,k,j)=alt(i,k,j)(t_2(i,k,j)-mu(i,j)t_1(i,k,j)) &
561	/(muts(i,j)*(t0+t_1(i,k,j)))-p(i,k,j)/c2a(i,k,j)
562	ph(i,k+1,j)=ph(i,k,j)-dnw(k)(muts(i,j)al (i,k,j) &
563	+mu(i,j)*alt(i,k,j))
564	ENDDO
565	ENDDO
566	ENDDO
567
568	END IF
569
570	! divergence damping setup
571
572	IF (step == 0) then ! we're initializing small timesteps
573	DO j=j_start, j_end
574	DO k=k_start, k_end
575	DO i=i_start, i_end
576	pm1(i,k,j)=p(i,k,j)
577	ENDDO
578	ENDDO
579	ENDDO
580	ELSE ! we're in the small timesteps
581	DO j=j_start, j_end ! and adding div damping component
582	DO k=k_start, k_end
583	DO i=i_start, i_end
584	ptmp = p(i,k,j)
585	p(i,k,j) = p(i,k,j) + smdiv*(p(i,k,j)-pm1(i,k,j))
586	pm1(i,k,j) = ptmp
587	ENDDO
588	ENDDO
589	ENDDO
590	END IF
591	#ifdef LMDZ1
592	PRINT *,' calc_p_rho: after divergence damping'
593	PRINT *,' pm1: ',pm1(im2,km2,jm2), ' smdiv: ', smdiv,' ptmp: ',ptmp
594	#endif
595
596	END SUBROUTINE calc_p_rho
597
598	!----------------------------------------------------------------------
599
600	SUBROUTINE calc_coef_w( a,alpha,gamma, &
601	mut, cqw, &
602	rdn, rdnw, c2a, &
603	dts, g, epssm, top_lid, &
604	ids,ide, jds,jde, kds,kde, & ! domain dims
605	ims,ime, jms,jme, kms,kme, & ! memory dims
606	its,ite, jts,jte, kts,kte ) ! tile dims
607
608	IMPLICIT NONE ! religion first
609
610	! passed in through the call
611
612	INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde
613	INTEGER, INTENT(IN ) :: ims,ime, jms,jme, kms,kme
614	INTEGER, INTENT(IN ) :: its,ite, jts,jte, kts,kte
615
616	LOGICAL, INTENT(IN ) :: top_lid
617
618	REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN ) :: c2a, &
619	cqw
620
621	REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: alpha, &
622	gamma, &
623	a
624
625	REAL, DIMENSION(ims:ime, jms:jme), INTENT(IN ) :: mut
626
627	REAL, DIMENSION(kms:kme), INTENT(IN ) :: rdn, &
628	rdnw
629
630	REAL, INTENT(IN ) :: epssm, &
631	dts, &
632	g
633
634	! Local stack data.
635
636	REAL, DIMENSION(ims:ime) :: cof
637	REAL :: b, c
638
639	INTEGER :: i, j, k, i_start, i_end, j_start, j_end, k_start, k_end
640	INTEGER :: ij, ijp, ijm, lid_flag
641
642	!<DESCRIPTION>
643	!
644	! calc_coef_w calculates the coefficients needed for the
645	! implicit solution of the vertical momentum and geopotential equations.
646	! This requires solution of a tri-diagonal equation.
647	!
648	!</DESCRIPTION>
649
650	i_start = its
651	i_end = min(ite,ide-1)
652	j_start = jts
653	j_end = min(jte,jde-1)
654	k_start = kts
655	k_end = kte-1
656
657	lid_flag=1
658	IF(top_lid)lid_flag=0
659
660	outer_j_loop: DO j = j_start, j_end
661
662	DO i = i_start, i_end
663	cof(i) = (.5dtsg(1.+epssm)/mut(i,j))*2
664	a(i, 2 ,j) = 0.
665	a(i,kde,j) =-2.cof(i)rdnw(kde-1)*2c2a(i,kde-1,j)*lid_flag
666	gamma(i,1 ,j) = 0.
667	ENDDO
668
669	DO k=3,kde-1
670	DO i=i_start, i_end
671	a(i,k,j) = -cqw(i,k,j)cof(i)rdn(k)* rdnw(k-1)*c2a(i,k-1,j)
672	ENDDO
673	ENDDO
674
675
676	DO k=2,kde-1
677	DO i=i_start, i_end
678	b = 1.+cqw(i,k,j)cof(i)rdn(k)(rdnw(k )c2a(i,k,j ) &
679	+rdnw(k-1)*c2a(i,k-1,j))
680	c = -cqw(i,k,j)cof(i)rdn(k)rdnw(k )c2a(i,k,j )
681	alpha(i,k,j) = 1./(b-a(i,k,j)*gamma(i,k-1,j))
682	gamma(i,k,j) = c*alpha(i,k,j)
683	ENDDO
684	ENDDO
685
686	DO i=i_start, i_end
687	b = 1.+2.cof(i)rdnw(kde-1)*2c2a(i,kde-1,j)
688	c = 0.
689	alpha(i,kde,j) = 1./(b-a(i,kde,j)*gamma(i,kde-1,j))
690	gamma(i,kde,j) = c*alpha(i,kde,j)
691	ENDDO
692
693	ENDDO outer_j_loop
694
695	END SUBROUTINE calc_coef_w
696
697	!----------------------------------------------------------------------
698
699	SUBROUTINE advance_uv ( u, ru_tend, v, rv_tend, &
700	p, pb, &
701	ph, php, alt, al, mu, &
702	muu, cqu, muv, cqv, mudf, &
703	msfux, msfuy, msfvx, &
704	msfvx_inv, msfvy, &
705	rdx, rdy, dts, &
706	cf1, cf2, cf3, fnm, fnp, &
707	emdiv, &
708	rdnw, config_flags, spec_zone, &
709	non_hydrostatic, top_lid, &
710	ids, ide, jds, jde, kds, kde, &
711	ims, ime, jms, jme, kms, kme, &
712	its, ite, jts, jte, kts, kte )
713
714
715
716	IMPLICIT NONE ! religion first
717
718	! stuff coming in
719
720	TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
721
722	LOGICAL, INTENT(IN ) :: non_hydrostatic, top_lid
723
724	INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde
725	INTEGER, INTENT(IN ) :: ims,ime, jms,jme, kms,kme
726	INTEGER, INTENT(IN ) :: its,ite, jts,jte, kts,kte
727	INTEGER, INTENT(IN ) :: spec_zone
728
729	REAL, DIMENSION( ims:ime , kms:kme, jms:jme ), &
730	INTENT(INOUT) :: &
731	u, &
732	v
733
734	REAL, DIMENSION( ims:ime , kms:kme, jms:jme ), &
735	INTENT(IN ) :: &
736	ru_tend, &
737	rv_tend, &
738	ph, &
739	php, &
740	p, &
741	pb, &
742	alt, &
743	al, &
744	cqu, &
745	cqv
746
747
748	REAL, DIMENSION( ims:ime , jms:jme ), INTENT(IN ) :: muu, &
749	muv, &
750	mu, &
751	mudf
752
753
754	REAL, DIMENSION( kms:kme ), INTENT(IN ) :: fnm, &
755	fnp , &
756	rdnw
757
758	REAL, DIMENSION( ims:ime , jms:jme ), INTENT(IN ) :: msfux, &
759	msfuy, &
760	msfvx, &
761	msfvy, &
762	msfvx_inv
763
764	REAL, INTENT(IN ) :: rdx, &
765	rdy, &
766	dts, &
767	cf1, &
768	cf2, &
769	cf3, &
770	emdiv
771
772
773	! Local 3d array from the stack (note tile size)
774
775	REAL, DIMENSION (its:ite, kts:kte) :: dpn, dpxy
776	REAL, DIMENSION (its:ite) :: mudf_xy
777	REAL :: dx, dy
778
779	INTEGER :: i,j,k, i_start, i_end, j_start, j_end, k_start, k_end
780	INTEGER :: i_endu, j_endv, k_endw
781	INTEGER :: i_start_up, i_end_up, j_start_up, j_end_up
782	INTEGER :: i_start_vp, i_end_vp, j_start_vp, j_end_vp
783
784	INTEGER :: i_start_u_tend, i_end_u_tend, j_start_v_tend, j_end_v_tend
785
786	#ifdef LMDZ1
787	INTEGER :: im2,km2,jm2
788
789	im2=24
790	jm2=21
791	km2=1
792
793	PRINT *,' advance_uv: inside'
794	PRINT *,' al: ',al(im2,km2,jm2), ' p: ', p(im2,1,jm2),' ph: ',ph(im2,1,jm2)
795	#endif
796
797	!<DESCRIPTION>
798	!
799	! advance_uv advances the explicit perturbation horizontal momentum
800	! equations (u,v) by adding in the large-timestep tendency along with
801	! the small timestep pressure gradient tendency.
802	!
803	!</DESCRIPTION>
804
805	! now, the real work.
806	! set the loop bounds taking into account boundary conditions.
807
808	IF( config_flags%nested .or. config_flags%specified ) THEN
809	i_start = max( its,ids+spec_zone )
810	i_end = min( ite,ide-spec_zone-1 )
811	j_start = max( jts,jds+spec_zone )
812	j_end = min( jte,jde-spec_zone-1 )
813	k_start = kts
814	k_end = min( kte, kde-1 )
815
816	i_endu = min( ite,ide-spec_zone )
817	j_endv = min( jte,jde-spec_zone )
818	k_endw = k_end
819
820	IF( config_flags%periodic_x) THEN
821	i_start = its
822	i_end = min(ite,ide-1)
823	i_endu = ite
824	ENDIF
825	ELSE
826	i_start = its
827	i_end = min(ite,ide-1)
828	j_start = jts
829	j_end = min(jte,jde-1)
830	k_start = kts
831	k_end = kte-1
832
833	i_endu = ite
834	j_endv = jte
835	k_endw = k_end
836	ENDIF
837
838	i_start_up = i_start
839	i_end_up = i_endu
840	j_start_up = j_start
841	j_end_up = j_end
842
843	i_start_vp = i_start
844	i_end_vp = i_end
845	j_start_vp = j_start
846	j_end_vp = j_endv
847
848	IF ( (config_flags%open_xs .or. &
849	config_flags%symmetric_xs ) &
850	.and. (its == ids) ) &
851	i_start_up = i_start_up + 1
852
853	IF ( (config_flags%open_xe .or. &
854	config_flags%symmetric_xe ) &
855	.and. (ite == ide) ) &
856	i_end_up = i_end_up - 1
857
858	IF ( (config_flags%open_ys .or. &
859	config_flags%symmetric_ys .or. &
860	config_flags%polar ) &
861	.and. (jts == jds) ) &
862	j_start_vp = j_start_vp + 1
863
864	IF ( (config_flags%open_ye .or. &
865	config_flags%symmetric_ye .or. &
866	config_flags%polar ) &
867	.and. (jte == jde) ) &
868	j_end_vp = j_end_vp - 1
869
870	i_start_u_tend = i_start
871	i_end_u_tend = i_endu
872	j_start_v_tend = j_start
873	j_end_v_tend = j_endv
874
875	IF ( config_flags%symmetric_xs .and. (its == ids) ) &
876	i_start_u_tend = i_start_u_tend+1
877	IF ( config_flags%symmetric_xe .and. (ite == ide) ) &
878	i_end_u_tend = i_end_u_tend-1
879	IF ( config_flags%symmetric_ys .and. (jts == jds) ) &
880	j_start_v_tend = j_start_v_tend+1
881	IF ( config_flags%symmetric_ye .and. (jte == jde) ) &
882	j_end_v_tend = j_end_v_tend-1
883
884	dx = 1./rdx
885	dy = 1./rdy
886
887	! start real calculations.
888	! first, u
889
890	u_outer_j_loop: DO j = j_start, j_end
891
892	DO k = k_start, k_end
893	DO i = i_start_u_tend, i_end_u_tend
894	u(i,k,j) = u(i,k,j) + dts*ru_tend(i,k,j)
895	ENDDO
896	ENDDO
897
898	DO i = i_start_up, i_end_up
899	mudf_xy(i)= -emdivdx(mudf(i,j)-mudf(i-1,j))/msfuy(i,j)
900	ENDDO
901
902	DO k = k_start, k_end
903	DO i = i_start_up, i_end_up
904
905	! Comments on map scale factors:
906	! x pressure gradient: ADT eqn 44, penultimate term on RHS
907	! = -(mx/my)(mu/rho)partial dp/dx
908	! [i.e., first rho->mu; 2nd still rho; alpha=1/rho]
909	! Klemp et al. splits into 2 terms:
910	! mu alpha partial dp/dx + partial dp/dnu * partial dphi/dx
911	! then into 4 terms:
912	! mu alpha partial dp'/dx + nu mu alpha' partial dmubar/dx +
913	! + mu partial dphi/dx + partial dphi'/dx * (partial dp'/dnu - mu')
914	!
915	! first 3 terms:
916	! ph, alt, p, al, pb not coupled
917	! since we want tendency to fit ADT eqn 44 (coupled) we need to
918	! multiply by (mx/my):
919	!
920	dpxy(i,k)= (msfux(i,j)/msfuy(i,j)).5rdxmuu(i,j)( &
921	((ph (i,k+1,j)-ph (i-1,k+1,j))+(ph (i,k,j)-ph (i-1,k,j))) &
922	+(alt(i,k ,j)+alt(i-1,k ,j))*(p (i,k,j)-p (i-1,k,j)) &
923	+(al (i,k ,j)+al (i-1,k ,j))*(pb (i,k,j)-pb (i-1,k,j)) )
924
925	ENDDO
926	ENDDO
927
928	IF (non_hydrostatic) THEN
929
930	DO i = i_start_up, i_end_up
931	dpn(i,1) = .5( cf1(p(i,1,j)+p(i-1,1,j)) &
932	+cf2*(p(i,2,j)+p(i-1,2,j)) &
933	+cf3*(p(i,3,j)+p(i-1,3,j)) )
934	dpn(i,kde) = 0.
935	ENDDO
936	IF (top_lid) THEN
937	DO i = i_start_up, i_end_up
938	dpn(i,kde) =.5( cf1(p(i-1,kde-1,j)+p(i,kde-1,j)) &
939	+cf2*(p(i-1,kde-2,j)+p(i,kde-2,j)) &
940	+cf3*(p(i-1,kde-3,j)+p(i,kde-3,j)) )
941	ENDDO
942	ENDIF
943
944	DO k = k_start+1, k_end
945	DO i = i_start_up, i_end_up
946	dpn(i,k) = .5( fnm(k)(p(i,k ,j)+p(i-1,k ,j)) &
947	+fnp(k)*(p(i,k-1,j)+p(i-1,k-1,j)) )
948	ENDDO
949	ENDDO
950
951	! Comments on map scale factors:
952	! 4th term:
953	! php, dpn, mu not coupled
954	! since we want tendency to fit ADT eqn 44 (coupled) we need to
955	! multiply by (mx/my):
956
957	DO k = k_start, k_end
958	DO i = i_start_up, i_end_up
959	dpxy(i,k)=dpxy(i,k) + (msfux(i,j)/msfuy(i,j))rdx(php(i,k,j)-php(i-1,k,j))* &
960	(rdnw(k)(dpn(i,k+1)-dpn(i,k))-.5(mu(i-1,j)+mu(i,j)))
961	ENDDO
962	ENDDO
963
964
965	END IF
966
967
968	DO k = k_start, k_end
969	DO i = i_start_up, i_end_up
970	u(i,k,j)=u(i,k,j)-dtscqu(i,k,j)dpxy(i,k)+mudf_xy(i)
971	ENDDO
972	ENDDO
973
974	ENDDO u_outer_j_loop
975
976	! now v
977
978	v_outer_j_loop: DO j = j_start_v_tend, j_end_v_tend
979
980
981	DO k = k_start, k_end
982	DO i = i_start, i_end
983	v(i,k,j) = v(i,k,j) + dts*rv_tend(i,k,j)
984	ENDDO
985	ENDDO
986
987	DO i = i_start, i_end
988	mudf_xy(i)= -emdivdy(mudf(i,j)-mudf(i,j-1))*msfvx_inv(i,j)
989	ENDDO
990
991	IF ( ( j >= j_start_vp) &
992	.and.( j <= j_end_vp ) ) THEN
993
994	DO k = k_start, k_end
995	DO i = i_start, i_end
996
997	! Comments on map scale factors:
998	! y pressure gradient: ADT eqn 45, penultimate term on RHS
999	! = -(my/mx)(mu/rho)partial dp/dy
1000	! [i.e., first rho->mu; 2nd still rho; alpha=1/rho]
1001	! Klemp et al. splits into 2 terms:
1002	! mu alpha partial dp/dy + partial dp/dnu * partial dphi/dy
1003	! then into 4 terms:
1004	! mu alpha partial dp'/dy + nu mu alpha' partial dmubar/dy +
1005	! + mu partial dphi/dy + partial dphi'/dy * (partial dp'/dnu - mu')
1006	!
1007	! first 3 terms:
1008	! ph, alt, p, al, pb not coupled
1009	! since we want tendency to fit ADT eqn 45 (coupled) we need to
1010	! multiply by (my/mx):
1011	! mudf_xy is NOT a map scale factor coupling
1012	! it is some sort of divergence damping
1013
1014	dpxy(i,k)= (msfvy(i,j)/msfvx(i,j)).5rdymuv(i,j)( &
1015	((ph(i,k+1,j)-ph(i,k+1,j-1))+(ph (i,k,j)-ph (i,k,j-1))) &
1016	+(alt(i,k ,j)+alt(i,k ,j-1))*(p (i,k,j)-p (i,k,j-1)) &
1017	+(al (i,k ,j)+al (i,k ,j-1))*(pb (i,k,j)-pb (i,k,j-1)) )
1018	ENDDO
1019	ENDDO
1020
1021
1022	IF (non_hydrostatic) THEN
1023
1024	DO i = i_start, i_end
1025	dpn(i,1) = .5( cf1(p(i,1,j)+p(i,1,j-1)) &
1026	+cf2*(p(i,2,j)+p(i,2,j-1)) &
1027	+cf3*(p(i,3,j)+p(i,3,j-1)) )
1028	dpn(i,kde) = 0.
1029	ENDDO
1030	IF (top_lid) THEN
1031	DO i = i_start, i_end
1032	dpn(i,kde) =.5( cf1(p(i,kde-1,j-1)+p(i,kde-1,j)) &
1033	+cf2*(p(i,kde-2,j-1)+p(i,kde-2,j)) &
1034	+cf3*(p(i,kde-3,j-1)+p(i,kde-3,j)) )
1035	ENDDO
1036	ENDIF
1037
1038	DO k = k_start+1, k_end
1039	DO i = i_start, i_end
1040	dpn(i,k) = .5( fnm(k)(p(i,k ,j)+p(i,k ,j-1)) &
1041	+fnp(k)*(p(i,k-1,j)+p(i,k-1,j-1)) )
1042	ENDDO
1043	ENDDO
1044
1045	! Comments on map scale factors:
1046	! 4th term:
1047	! php, dpn, mu not coupled
1048	! since we want tendency to fit ADT eqn 45 (coupled) we need to
1049	! multiply by (my/mx):
1050
1051	DO k = k_start, k_end
1052	DO i = i_start, i_end
1053	dpxy(i,k)=dpxy(i,k) + (msfvy(i,j)/msfvx(i,j))rdy(php(i,k,j)-php(i,k,j-1))* &
1054	(rdnw(k)(dpn(i,k+1)-dpn(i,k))-.5(mu(i,j-1)+mu(i,j)))
1055	ENDDO
1056	ENDDO
1057
1058
1059	END IF
1060
1061
1062	DO k = k_start, k_end
1063	DO i = i_start, i_end
1064	v(i,k,j)=v(i,k,j)-dtscqv(i,k,j)dpxy(i,k)+mudf_xy(i)
1065	ENDDO
1066	ENDDO
1067	END IF
1068
1069	ENDDO v_outer_j_loop
1070
1071	! The check for j_start_vp and j_end_vp makes sure that the edges in v
1072	! are not updated. Let's add a polar boundary condition check here for
1073	! safety to ensure that v at the poles never gets to be non-zero in the
1074	! small time steps.
1075	IF (config_flags%polar) THEN
1076	IF (jts == jds) THEN
1077	DO k = k_start, k_end
1078	DO i = i_start, i_end
1079	v(i,k,jds) = 0.
1080	ENDDO
1081	ENDDO
1082	END IF
1083	IF (jte == jde) THEN
1084	DO k = k_start, k_end
1085	DO i = i_start, i_end
1086	v(i,k,jde) = 0.
1087	ENDDO
1088	ENDDO
1089	END IF
1090	END IF
1091
1092
1093	END SUBROUTINE advance_uv
1094
1095	!---------------------------------------------------------------------
1096
1097	SUBROUTINE advance_mu_t( ww, ww_1, u, u_1, v, v_1, &
1098	mu, mut, muave, muts, muu, muv, &
1099	mudf, uam, vam, wwam, t, t_1, &
1100	t_ave, ft, mu_tend, &
1101	rdx, rdy, dts, epssm, &
1102	dnw, fnm, fnp, rdnw, &
1103	msfux, msfuy, msfvx, msfvx_inv, &
1104	msfvy, msftx, msfty, &
1105	step, config_flags, &
1106	ids, ide, jds, jde, kds, kde, &
1107	ims, ime, jms, jme, kms, kme, &
1108	its, ite, jts, jte, kts, kte )
1109
1110	IMPLICIT NONE ! religion first
1111
1112	! stuff coming in
1113
1114	TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
1115
1116	INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde
1117	INTEGER, INTENT(IN ) :: ims,ime, jms,jme, kms,kme
1118	INTEGER, INTENT(IN ) :: its,ite, jts,jte, kts,kte
1119
1120	INTEGER, INTENT(IN ) :: step
1121
1122	REAL, DIMENSION( ims:ime , kms:kme, jms:jme ), &
1123	INTENT(IN ) :: &
1124	u, &
1125	v, &
1126	u_1, &
1127	v_1, &
1128	t_1, &
1129	ft
1130
1131	REAL, DIMENSION( ims:ime , kms:kme, jms:jme ), &
1132	INTENT(INOUT) :: &
1133	ww, &
1134	ww_1, &
1135	t, &
1136	t_ave, &
1137	uam, &
1138	vam, &
1139	wwam
1140
1141	REAL, DIMENSION( ims:ime , jms:jme ), INTENT(IN ) :: muu, &
1142	muv, &
1143	mut, &
1144	msfux,&
1145	msfuy,&
1146	msfvx,&
1147	msfvx_inv,&
1148	msfvy,&
1149	msftx,&
1150	msfty,&
1151	mu_tend
1152
1153	REAL, DIMENSION( ims:ime , jms:jme ), INTENT( OUT) :: muave, &
1154	muts, &
1155	mudf
1156
1157	REAL, DIMENSION( ims:ime , jms:jme ), INTENT(INOUT) :: mu
1158
1159	REAL, DIMENSION( kms:kme ), INTENT(IN ) :: fnm, &
1160	fnp, &
1161	dnw, &
1162	rdnw
1163
1164
1165	REAL, INTENT(IN ) :: rdx, &
1166	rdy, &
1167	dts, &
1168	epssm
1169
1170	! Local arrays from the stack (note tile size)
1171
1172	REAL, DIMENSION (its:ite, kts:kte) :: wdtn, dvdxi
1173	REAL, DIMENSION (its:ite) :: dmdt
1174
1175	INTEGER :: i,j,k, i_start, i_end, j_start, j_end, k_start, k_end
1176	INTEGER :: i_endu, j_endv
1177	REAL :: acc
1178
1179	!<DESCRIPTION>
1180	!
1181	! advance_mu_t advances the explicit perturbation theta equation and the mass
1182	! conservation equation. In addition, the small timestep omega is updated,
1183	! and some quantities needed in other places are squirrelled away.
1184	!
1185	!</DESCRIPTION>
1186
1187	! now, the real work.
1188	! set the loop bounds taking into account boundary conditions.
1189	#ifdef LMDZ1
1190	INTEGER :: im2,km2,jm2
1191
1192	im2 = config_flags%i_check_point
1193	jm2 = config_flags%j_check_point
1194	km2 = config_flags%k_check_point
1195
1196	PRINT *,' advance_mu_t: inside'
1197	PRINT *,' ww: ',ww(im2,km2,jm2),' ww_1: ', ww_1(im2,km2,jm2), ' u: ',u(im2,km2,jm2), &
1198	' u_1: ',u_1(im2,km2,jm2), ' v: ',v(im2,km2,jm2), ' v_1: ',v_1(im2,km2,jm2), ' mu: ', mu(im2,jm2), &
1199	' mu_t: ',mut(im2,jm2), ' muave: ',muave(im2,jm2), ' muts: ',muts(im2,jm2), &
1200	' muu: ',muu(im2,jm2), ' muv: ',muv(im2,jm2), ' mudf: ', mudf(im2,jm2), &
1201	' uam: ',uam(im2,km2,jm2), ' vam: ',vam(im2,km2,jm2), ' wwam: ',wwam(im2,km2,jm2), ' t: ',t(im2,km2,jm2), &
1202	' t_1: ',t_1(im2,km2,jm2), ' t_ave: ', t_ave(im2,km2,jm2), ' ft: ',ft(im2,km2,jm2), &
1203	' mu_tend: ', mu_tend(im2,jm2), ' rdx: ', rdx, ' rdy: ',rdy, &
1204	' dts: ',dts, ' epssm: ',epssm, ' dnw: ', dnw(km2), &
1205	' fnm: ',fnm(km2), ' fnp: ',fnp(km2), ' rndw: ',rdnw(km2), &
1206	' msfux: ',msfux(im2,jm2), ' msfuy: ',msfuy(im2,jm2), 'msfvx: ',msfvx(im2,jm2), &
1207	' msfvx_inv: ',msfvx_inv(im2,jm2), 'msfvy: ',msfvy(im2,jm2), &
1208	' msftx: ',msftx(im2,jm2), 'msfty: ',msfty(im2,jm2), ' step: ', step
1209	#endif
1210
1211	i_start = its
1212	i_end = min(ite,ide-1)
1213	j_start = jts
1214	j_end = min(jte,jde-1)
1215	k_start = kts
1216	k_end = kte-1
1217	IF ( .NOT. config_flags%periodic_x )THEN
1218	IF ( config_flags%specified .or. config_flags%nested ) then
1219	i_start = max(its,ids+1)
1220	i_end = min(ite,ide-2)
1221	ENDIF
1222	ENDIF
1223	IF ( config_flags%specified .or. config_flags%nested ) then
1224	j_start = max(jts,jds+1)
1225	j_end = min(jte,jde-2)
1226	ENDIF
1227
1228	i_endu = ite
1229	j_endv = jte
1230
1231	! CALCULATION OF WW (dETA/dt)
1232	DO j = j_start, j_end
1233
1234	DO i=i_start, i_end
1235	dmdt(i) = 0.
1236	ENDDO
1237	! NOTE: mu is not coupled with the map scale factor.
1238	! ww (omega) IS coupled with the map scale factor.
1239	! Being coupled with the map scale factor means
1240	! multiplication by (1/msft) in this case.
1241
1242	! Comments on map scale factors
1243	! ADT eqn 47:
1244	! partial drho/dt = -mx*my[partial d/dx(rho u/my) + partial d/dy(rho v/mx)]
1245	! -partial d/dz(rho w)
1246	! with rho -> mu, dividing by my, and with partial d/dnu(rho nu/my [=ww])
1247	! as the final term (because we're looking for d_nu_/dt)
1248	!
1249	! begin by integrating with respect to nu from bottom to top
1250	! BCs are ww=0 at both
1251	! final term gives 0
1252	! first term gives Integral([1/my]partial d mu/dt) over total column = dm/dt
1253	! RHS remaining is Integral(-mx[partial d/dx(mu u/my) +
1254	! partial d/dy(mu v/mx)]) over column
1255	! lines below find RHS terms at each level then set dmdt = sum over all levels
1256	!
1257	! [don't divide the below by msfty until find ww, since dmdt is used in
1258	! the meantime]
1259
1260	#ifdef LMDZ1
1261	IF (j == jm2) PRINT *,' Lluis advance_mu_t k msftx msfty rdy v/m muv v_1 ' // &
1262	'msfvx_inv rdx u/m muu u_1 msfuy dvdxi dnw ______ '
1263	#endif
1264	DO k=k_start, k_end
1265	DO i=i_start, i_end
1266	dvdxi(i,k) = msftx(i,j)msfty(i,j)( &
1267	rdy( (v(i,k,j+1)+muv(i,j+1)v_1(i,k,j+1)*msfvx_inv(i,j+1)) &
1268	-(v(i,k,j )+muv(i,j )v_1(i,k,j )msfvx_inv(i,j )) ) &
1269	+rdx( (u(i+1,k,j)+muu(i+1,j)u_1(i+1,k,j)/msfuy(i+1,j)) &
1270	-(u(i,k,j )+muu(i ,j)*u_1(i,k,j )/msfuy(i ,j)) ))
1271	dmdt(i) = dmdt(i) + dnw(k)*dvdxi(i,k)
1272	#ifdef LMDZ1
1273	IF (i == im2 .AND. j == jm2) PRINT *,k, msftx(im2,jm2), msfty(im2,jm2), rdy, v(im2,k,jm2), muv(im2,jm2), &
1274	v_1(im2,k,jm2), msfvx_inv(im2,jm2), rdx, u(im2,k,jm2), muu(im2,jm2), u_1(im2,k,jm2), &
1275	msfuy(im2,jm2), dvdxi(im2,k), dnw(k)
1276	#endif
1277	ENDDO
1278	ENDDO
1279	#ifdef LMDZ1
1280	IF (j == jm2) THEN
1281	PRINT *,' advance_mu_t: after dvdx/dmdt'
1282	PRINT *,' dvdxi: ',dvdxi(im2,km2),' dmdt: ', dmdt(im2)
1283	END IF
1284	#endif
1285
1286	DO i=i_start, i_end
1287	muave(i,j) = mu(i,j)
1288	mu(i,j) = mu(i,j)+dts*(dmdt(i)+mu_tend(i,j))
1289	mudf(i,j) = (dmdt(i)+mu_tend(i,j)) ! save tendency for div damp filter
1290	muts(i,j) = mut(i,j)+mu(i,j)
1291	muave(i,j) =.5((1.+epssm)mu(i,j)+(1.-epssm)*muave(i,j))
1292	ENDDO
1293	#ifdef LMDZ1
1294	IF (j == jm2) THEN
1295	PRINT *,' advance_mu_t: after muave'
1296	PRINT *,' muave: ',muave(im2,jm2),' mu: ', mu(im2,jm2), ' mudf: ',mudf(im2,jm2), &
1297	' muts: ',muts(im2,jm2), ' muave: ',muave(im2,jm2)
1298	END IF
1299	#endif
1300
1301	DO k=2,k_end
1302	DO i=i_start, i_end
1303	ww(i,k,j)=ww(i,k-1,j)-dnw(k-1)*(dmdt(i)+dvdxi(i,k-1)+mu_tend(i,j))/msfty(i,j)
1304	ENDDO
1305	END DO
1306	#ifdef LMDZ1
1307	IF (j == jm2) THEN
1308	PRINT *,' advance_mu_t: after w'
1309	PRINT *,' ww: ',ww(im2,km2,jm2)
1310	END IF
1311	#endif
1312
1313	! NOTE: ww_1 (large timestep ww) is already coupled with the
1314	! map scale factor
1315
1316	DO k=1,k_end
1317	DO i=i_start, i_end
1318	ww(i,k,j)=ww(i,k,j)-ww_1(i,k,j)
1319	END DO
1320	END DO
1321	#ifdef LMDZ1
1322	IF (j == jm2) THEN
1323	PRINT *,' advance_mu_t: after w 1'
1324	PRINT *,' ww: ',ww(im2,km2,jm2)
1325	END IF
1326	#endif
1327	ENDDO
1328
1329	! CALCULATION OF THETA
1330
1331	! NOTE: theta'' is not coupled with the map-scale factor,
1332	! while the theta'' tendency is coupled (i.e., mult by 1/msft)
1333
1334	! Comments on map scale factors
1335	! BUT NOTE THAT both are mass coupled
1336	! in flux form equations (Klemp et al.) Theta = mu*theta
1337	!
1338	! scalar eqn: partial d/dt(rho q/my) = -mx[partial d/dx(q rho u/my) +
1339	! partial d/dy(q rho v/mx)]
1340	! - partial d/dz(q rho w/my)
1341	! with rho -> mu, and with partial d/dnu(q rho nu/my) as the final term
1342	!
1343	! adding previous tendency contribution which was map scale factor coupled
1344	! (had been divided by msfty)
1345	! need to uncouple before updating uncoupled Theta (by adding)
1346
1347	DO j=j_start, j_end
1348	DO k=1,k_end
1349	DO i=i_start, i_end
1350	t_ave(i,k,j) = t(i,k,j)
1351	t (i,k,j) = t(i,k,j) + msfty(i,j)dtsft(i,k,j)
1352	END DO
1353	END DO
1354	ENDDO
1355
1356	DO j=j_start, j_end
1357
1358	DO i=i_start, i_end
1359	wdtn(i,1 )=0.
1360	wdtn(i,kde)=0.
1361	ENDDO
1362
1363	DO k=2,k_end
1364	DO i=i_start, i_end
1365	! for scalar eqn RHS term 3
1366	wdtn(i,k)= ww(i,k,j)(fnm(k)t_1(i,k ,j)+fnp(k)*t_1(i,k-1,j))
1367	ENDDO
1368	ENDDO
1369
1370	! scalar eqn, RHS terms 1, 2 and 3
1371	! multiply by msfty to uncouple result for Theta from map scale factor
1372
1373	DO k=1,k_end
1374	DO i=i_start, i_end
1375	! multiplication by msfty uncouples result for Theta
1376	t(i,k,j) = t(i,k,j) - dtsmsfty(i,j)( &
1377	! multiplication by mx needed for RHS terms 1 & 2
1378	msftx(i,j)*( &
1379	.5rdy &
1380	( v(i,k,j+1)*(t_1(i,k,j+1)+t_1(i,k, j )) &
1381	-v(i,k,j )*(t_1(i,k, j )+t_1(i,k,j-1)) ) &
1382	+ .5rdx &
1383	( u(i+1,k,j)*(t_1(i+1,k,j)+t_1(i ,k,j)) &
1384	-u(i ,k,j)*(t_1(i ,k,j)+t_1(i-1,k,j)) ) ) &
1385	+ rdnw(k)*( wdtn(i,k+1)-wdtn(i,k) ) )
1386	ENDDO
1387	ENDDO
1388
1389	ENDDO
1390
1391	END SUBROUTINE advance_mu_t
1392
1393
1394
1395	!------------------------------------------------------------
1396
1397	SUBROUTINE advance_w( w, rw_tend, ww, w_save, u, v, &
1398	mu1, mut, muave, muts, &
1399	t_2ave, t_2, t_1, &
1400	ph, ph_1, phb, ph_tend, &
1401	ht, c2a, cqw, alt, alb, &
1402	a, alpha, gamma, &
1403	rdx, rdy, dts, t0, epssm, &
1404	dnw, fnm, fnp, rdnw, rdn, &
1405	cf1, cf2, cf3, msftx, msfty,&
1406	config_flags, top_lid, &
1407	ids,ide, jds,jde, kds,kde, & ! domain dims
1408	ims,ime, jms,jme, kms,kme, & ! memory dims
1409	its,ite, jts,jte, kts,kte ) ! tile dims
1410
1411	! We have used msfty for msft_inv but have not thought through w equation
1412	! pieces properly yet, so we will have to hope that it is okay
1413	! We think we have found a slight error in surface w calculation
1414
1415	IMPLICIT NONE ! religion first
1416
1417	! stuff coming in
1418
1419
1420	TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
1421
1422	INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde
1423	INTEGER, INTENT(IN ) :: ims,ime, jms,jme, kms,kme
1424	INTEGER, INTENT(IN ) :: its,ite, jts,jte, kts,kte
1425
1426	LOGICAL, INTENT(IN ) :: top_lid
1427
1428	REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
1429	INTENT(INOUT) :: &
1430	t_2ave, &
1431	w, &
1432	ph
1433
1434	REAL, DIMENSION( ims:ime , kms:kme, jms:jme ), &
1435	INTENT(IN ) :: &
1436	rw_tend, &
1437	ww, &
1438	w_save, &
1439	u, &
1440	v, &
1441	t_2, &
1442	t_1, &
1443	ph_1, &
1444	phb, &
1445	ph_tend, &
1446	alpha, &
1447	gamma, &
1448	a, &
1449	c2a, &
1450	cqw, &
1451	alb, &
1452	alt
1453
1454	REAL, DIMENSION( ims:ime , jms:jme ), &
1455	INTENT(IN ) :: &
1456	mu1, &
1457	mut, &
1458	muave, &
1459	muts, &
1460	ht, &
1461	msftx, &
1462	msfty
1463
1464	REAL, DIMENSION( kms:kme ), INTENT(IN ) :: fnp, &
1465	fnm, &
1466	rdnw, &
1467	rdn, &
1468	dnw
1469
1470	REAL, INTENT(IN ) :: rdx, &
1471	rdy, &
1472	dts, &
1473	cf1, &
1474	cf2, &
1475	cf3, &
1476	t0, &
1477	epssm
1478
1479	! Stack based 3d data, tile size.
1480
1481	REAL, DIMENSION( its:ite ) :: mut_inv, msft_inv
1482	REAL, DIMENSION( its:ite, kts:kte ) :: rhs, wdwn
1483	INTEGER :: i,j,k, i_start, i_end, j_start, j_end, k_start, k_end
1484	REAL, DIMENSION (kts:kte) :: dampwt
1485	real :: htop,hbot,hdepth,hk
1486	real :: pi,dampmag
1487
1488	!<DESCRIPTION>
1489	!
1490	! advance_w advances the implicit w and geopotential equations.
1491	!
1492	!</DESCRIPTION>
1493
1494	! set loop limits.
1495	! Currently set for periodic boundary conditions
1496
1497	i_start = its
1498	i_end = min(ite,ide-1)
1499	j_start = jts
1500	j_end = min(jte,jde-1)
1501	k_start = kts
1502	k_end = kte-1
1503	IF ( .NOT. config_flags%periodic_x )THEN
1504	IF ( config_flags%specified .or. config_flags%nested ) then
1505	i_start = max(its,ids+1)
1506	i_end = min(ite,ide-2)
1507	ENDIF
1508	ENDIF
1509	IF ( config_flags%specified .or. config_flags%nested ) then
1510	j_start = max(jts,jds+1)
1511	j_end = min(jte,jde-2)
1512	ENDIF
1513
1514	pi = 4.*atan(1.)
1515	dampmag = dts*config_flags%dampcoef
1516	hdepth=config_flags%zdamp
1517
1518	! calculation of phi and w equations
1519
1520	! Comments on map scale factors:
1521	! phi equation is:
1522	! partial d/dt(rho phi/my) = -mx partial d/dx(phi rho u/my)
1523	! -mx partial d/dy(phi rho v/mx)
1524	! - partial d/dz(phi rho w/my) + rho g w/my
1525	! as with scalar equation, use uncoupled value (here phi) to find the
1526	! coupled tendency (rho phi/my)
1527	! here as usual rho -> ~'mu'
1528	!
1529	! w eqn [divided by my] is:
1530	! partial d/dt(rho w/my) = -mx partial d/dx(w rho u/my)
1531	! -mx partial d/dy(v rho v/mx)
1532	! - partial d/dz(w rho w/my)
1533	! +rho[(uu+vv)/a + 2 u omega cos(lat) -
1534	! (1/rho) partial dp/dz - g + Fz]/my
1535	! here as usual rho -> ~'mu'
1536	!
1537	! 'u,v,w' sent here must be coupled variables (= rho w/my etc.) by their usage
1538
1539
1540	DO i=i_start, i_end
1541	rhs(i,1) = 0.
1542	ENDDO
1543
1544	j_loop_w: DO j = j_start, j_end
1545	DO i=i_start, i_end
1546	mut_inv(i) = 1./mut(i,j)
1547	msft_inv(i) = 1./msfty(i,j)
1548	ENDDO
1549
1550	DO k=1, k_end
1551	DO i=i_start, i_end
1552	t_2ave(i,k,j)=.5((1.+epssm)t_2(i,k,j) &
1553	+(1.-epssm)*t_2ave(i,k,j))
1554	t_2ave(i,k,j)=(t_2ave(i,k,j) + muave(i,j)*t0) &
1555	/(muts(i,j)*(t0+t_1(i,k,j)))
1556	wdwn(i,k+1)=.5(ww(i,k+1,j)+ww(i,k,j))rdnw(k) &
1557	*(ph_1(i,k+1,j)-ph_1(i,k,j)+phb(i,k+1,j)-phb(i,k,j))
1558	rhs(i,k+1) = dts(ph_tend(i,k+1,j) + .5g(1.-epssm)w(i,k+1,j))
1559
1560	ENDDO
1561	ENDDO
1562
1563	! building up RHS of phi equation
1564	! ph_tend contains terms 1 and 2; now adding 3 and 4 in stages:
1565	! here rhs = delta t [ph_tend + ~gw/2 - ~ww partial d phi/dz]
1566	DO k=2,k_end
1567	DO i=i_start, i_end
1568	rhs(i,k) = rhs(i,k)-dts( fnm(k)wdwn(i,k+1) &
1569	+fnp(k)*wdwn(i,k ) )
1570	ENDDO
1571	ENDDO
1572
1573	! NOTE: phi'' is not coupled with the map-scale factor (1/m),
1574	! but it's tendency is, so must multiply by msft here
1575
1576	! Comments on map scale factors:
1577	! building up RHS of phi equation
1578	! ph_tend contains terms 1 and 2; now adding 3 and 4 in stages:
1579	! here rhs = ph_previous + (msft/mu)[rhs_previous - ~ww delta t *
1580	! partial d phi/dz]
1581	! = ph_previous + (msftdelta t/mu)[ph_tend + ~gw/2 - ~ww
1582	! partial d phi/dz]
1583	DO k=2,k_end+1
1584	DO i=i_start, i_end
1585	rhs(i,k) = ph(i,k,j) + msfty(i,j)rhs(i,k)mut_inv(i)
1586	if(top_lid .and. k.eq.k_end+1)rhs(i,k)=0.
1587	ENDDO
1588	ENDDO
1589
1590	! lower boundary condition on w
1591
1592	! Comments on map scale factors:
1593	! Chain rule: if Z=Z(X,Y) [true at the surface] then
1594	! dZ/dt = dZ/dX * dX/dt + dZ/dY * dY/dt, U=dX/dt, V=dY/dt
1595	! using capitals to denote actual values
1596	! In mapped values, u=U, v=V, z=Z, 1/dX=mx/dx, 1/dY=my/dy
1597	! w = dz/dt = mx u dz/dx + my v dz/dy
1598	! [where dz/dx is just the surface height change between x
1599	! gridpoints, and dz/dy is the change between y gridpoints]
1600	! [cf1, cf2 and cf3 do vertical weighting of u or v values nearest
1601	! the surface]
1602
1603	DO i=i_start, i_end
1604	w(i,1,j)= &
1605
1606	msfty(i,j).5rdy*( &
1607	(ht(i,j+1)-ht(i,j )) &
1608	(cf1v(i,1,j+1)+cf2v(i,2,j+1)+cf3v(i,3,j+1)) &
1609	+(ht(i,j )-ht(i,j-1)) &
1610	(cf1v(i,1,j )+cf2v(i,2,j )+cf3v(i,3,j )) ) &
1611
1612	+msftx(i,j).5rdx*( &
1613	(ht(i+1,j)-ht(i,j )) &
1614	(cf1u(i+1,1,j)+cf2u(i+1,2,j)+cf3u(i+1,3,j)) &
1615	+(ht(i,j )-ht(i-1,j)) &
1616	(cf1u(i ,1,j)+cf2u(i ,2,j)+cf3u(i ,3,j)) )
1617
1618	ENDDO
1619	!
1620	! Jammed 3 doubly nested loops over k/i into 1 for slight improvement
1621	! in efficiency. No change in results (bit-for-bit). JM 20040514
1622	! (left a blank line where the other two k/i-loops were)
1623	!
1624	! above surface, begin by adding delta t * previous (coupled) w tendency
1625	DO k=2,k_end
1626	DO i=i_start, i_end
1627	w(i,k,j)=w(i,k,j)+dts*rw_tend(i,k,j) &
1628	+ msft_inv(i)cqw(i,k,j)( &
1629	+.5dtsgmut_inv(i)rdn(k)* &
1630	(c2a(i,k ,j)*rdnw(k ) &
1631	((1.+epssm)(rhs(i,k+1 )-rhs(i,k )) &
1632	+(1.-epssm)*(ph(i,k+1,j)-ph(i,k ,j))) &
1633	-c2a(i,k-1,j)*rdnw(k-1) &
1634	((1.+epssm)(rhs(i,k )-rhs(i,k-1 )) &
1635	+(1.-epssm)*(ph(i,k ,j)-ph(i,k-1,j))))) &
1636
1637	+dtsgmsft_inv(i)(rdn(k) &
1638	(c2a(i,k ,j)alt(i,k ,j)t_2ave(i,k ,j) &
1639	-c2a(i,k-1,j)alt(i,k-1,j)t_2ave(i,k-1,j)) &
1640	- muave(i,j))
1641	ENDDO
1642	ENDDO
1643
1644	K=k_end+1
1645
1646	DO i=i_start, i_end
1647	w(i,k,j)=w(i,k,j)+dts*rw_tend(i,k,j) &
1648	+msft_inv(i)*( &
1649	-.5dtsgmut_inv(i)rdnw(k-1)*22.*c2a(i,k-1,j) &
1650	((1.+epssm)(rhs(i,k )-rhs(i,k-1 )) &
1651	+(1.-epssm)*(ph(i,k,j)-ph(i,k-1,j))) &
1652	-dtsg(2.rdnw(k-1) &
1653	c2a(i,k-1,j)alt(i,k-1,j)t_2ave(i,k-1,j) &
1654	+ muave(i,j)) )
1655	if(top_lid)w(i,k,j) = 0.
1656	ENDDO
1657
1658	DO k=2,k_end+1
1659	DO i=i_start, i_end
1660	w(i,k,j)=(w(i,k,j)-a(i,k,j)w(i,k-1,j))alpha(i,k,j)
1661	ENDDO
1662	ENDDO
1663
1664	DO k=k_end,2,-1
1665	DO i=i_start, i_end
1666	w (i,k,j)=w (i,k,j)-gamma(i,k,j)*w(i,k+1,j)
1667	ENDDO
1668	ENDDO
1669
1670	IF (config_flags%damp_opt .eq. 3) THEN
1671	DO k=k_end+1,2,-1
1672	DO i=i_start, i_end
1673	htop=(ph_1(i,k_end+1,j)+phb(i,k_end+1,j))/g
1674	hk=(ph_1(i,k,j)+phb(i,k,j))/g
1675	hbot=htop-hdepth
1676	dampwt(k) = 0.
1677	if(hk .ge. hbot)then
1678	dampwt(k) = dampmagsin(0.5pi(hk-hbot)/hdepth)sin(0.5pi(hk-hbot)/hdepth)
1679	endif
1680	w(i,k,j) = (w(i,k,j) - dampwt(k)mut(i,j)w_save(i,k,j))/(1.+dampwt(k))
1681	ENDDO
1682	ENDDO
1683	ENDIF
1684
1685	DO k=k_end+1,2,-1
1686	DO i=i_start, i_end
1687	ph(i,k,j) = rhs(i,k)+msfty(i,j).5dtsg(1.+epssm) &
1688	*w(i,k,j)/muts(i,j)
1689	ENDDO
1690	ENDDO
1691
1692	ENDDO j_loop_w
1693
1694	END SUBROUTINE advance_w
1695
1696	!---------------------------------------------------------------------
1697
1698	SUBROUTINE sumflux ( ru, rv, ww, &
1699	u_lin, v_lin, ww_lin, &
1700	muu, muv, &
1701	ru_m, rv_m, ww_m, epssm, &
1702	msfux, msfuy, msfvx, msfvx_inv, msfvy, &
1703	iteration , number_of_small_timesteps, &
1704	ids,ide, jds,jde, kds,kde, &
1705	ims,ime, jms,jme, kms,kme, &
1706	its,ite, jts,jte, kts,kte )
1707
1708
1709	IMPLICIT NONE ! religion first
1710
1711	! declarations for the stuff coming in
1712
1713	INTEGER, INTENT(IN ) :: number_of_small_timesteps
1714	INTEGER, INTENT(IN ) :: iteration
1715	INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde
1716	INTEGER, INTENT(IN ) :: ims,ime, jms,jme, kms,kme
1717	INTEGER, INTENT(IN ) :: its,ite, jts,jte, kts,kte
1718
1719	REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN ) :: ru, &
1720	rv, &
1721	ww, &
1722	u_lin, &
1723	v_lin, &
1724	ww_lin
1725
1726
1727	REAL, DIMENSION(ims:ime, kms:kme, jms:jme) , INTENT(INOUT) :: ru_m, &
1728	rv_m, &
1729	ww_m
1730	REAL, DIMENSION(ims:ime, jms:jme) , INTENT(IN ) :: muu, muv, &
1731	msfux, msfuy, &
1732	msfvx, msfvy, msfvx_inv
1733
1734	INTEGER :: mini, minj, mink
1735
1736
1737	REAL, INTENT(IN ) :: epssm
1738	INTEGER :: i,j,k
1739
1740
1741	!<DESCRIPTION>
1742	!
1743	! update the small-timestep time-averaged mass fluxes; these
1744	! are needed for consistent mass-conserving scalar advection.
1745	!
1746	!</DESCRIPTION>
1747
1748	IF (iteration == 1 )THEN
1749	DO j = jts, jte
1750	DO k = kts, kte
1751	DO i = its, ite
1752	ru_m(i,k,j) = 0.
1753	rv_m(i,k,j) = 0.
1754	ww_m(i,k,j) = 0.
1755	ENDDO
1756	ENDDO
1757	ENDDO
1758	ENDIF
1759
1760	mini = min(ide-1,ite)
1761	minj = min(jde-1,jte)
1762	mink = min(kde-1,kte)
1763
1764
1765	DO j = jts, minj
1766	DO k = kts, mink
1767	DO i = its, mini
1768	ru_m(i,k,j) = ru_m(i,k,j) + ru(i,k,j)
1769	rv_m(i,k,j) = rv_m(i,k,j) + rv(i,k,j)
1770	ww_m(i,k,j) = ww_m(i,k,j) + ww(i,k,j)
1771	ENDDO
1772	ENDDO
1773	ENDDO
1774
1775	IF (ite .GT. mini) THEN
1776	DO j = jts, minj
1777	DO k = kts, mink
1778	DO i = mini+1, ite
1779	ru_m(i,k,j) = ru_m(i,k,j) + ru(i,k,j)
1780	ENDDO
1781	ENDDO
1782	ENDDO
1783	END IF
1784	IF (jte .GT. minj) THEN
1785	DO j = minj+1, jte
1786	DO k = kts, mink
1787	DO i = its, mini
1788	rv_m(i,k,j) = rv_m(i,k,j) + rv(i,k,j)
1789	ENDDO
1790	ENDDO
1791	ENDDO
1792	END IF
1793	IF ( kte .GT. mink) THEN
1794	DO j = jts, minj
1795	DO k = mink+1, kte
1796	DO i = its, mini
1797	ww_m(i,k,j) = ww_m(i,k,j) + ww(i,k,j)
1798	ENDDO
1799	ENDDO
1800	ENDDO
1801	END IF
1802
1803	IF (iteration == number_of_small_timesteps) THEN
1804
1805	DO j = jts, minj
1806	DO k = kts, mink
1807	DO i = its, mini
1808	ru_m(i,k,j) = ru_m(i,k,j) / number_of_small_timesteps &
1809	+ muu(i,j)*u_lin(i,k,j)/msfuy(i,j)
1810	rv_m(i,k,j) = rv_m(i,k,j) / number_of_small_timesteps &
1811	+ muv(i,j)v_lin(i,k,j)msfvx_inv(i,j)
1812	ww_m(i,k,j) = ww_m(i,k,j) / number_of_small_timesteps &
1813	+ ww_lin(i,k,j)
1814	ENDDO
1815	ENDDO
1816	ENDDO
1817
1818
1819	IF (ite .GT. mini) THEN
1820	DO j = jts, minj
1821	DO k = kts, mink
1822	DO i = mini+1, ite
1823	ru_m(i,k,j) = ru_m(i,k,j) / number_of_small_timesteps &
1824	+ muu(i,j)*u_lin(i,k,j)/msfuy(i,j)
1825	ENDDO
1826	ENDDO
1827	ENDDO
1828	END IF
1829	IF (jte .GT. minj) THEN
1830	DO j = minj+1, jte
1831	DO k = kts, mink
1832	DO i = its, mini
1833	rv_m(i,k,j) = rv_m(i,k,j) / number_of_small_timesteps &
1834	+ muv(i,j)v_lin(i,k,j)msfvx_inv(i,j)
1835	ENDDO
1836	ENDDO
1837	ENDDO
1838	END IF
1839	IF ( kte .GT. mink) THEN
1840	DO j = jts, minj
1841	DO k = mink+1, kte
1842	DO i = its, mini
1843	ww_m(i,k,j) = ww_m(i,k,j) / number_of_small_timesteps &
1844	+ ww_lin(i,k,j)
1845	ENDDO
1846	ENDDO
1847	ENDDO
1848	END IF
1849
1850	ENDIF
1851
1852
1853	END SUBROUTINE sumflux
1854
1855	!---------------------------------------------------------------------
1856
1857	SUBROUTINE init_module_small_step
1858	END SUBROUTINE init_module_small_step
1859
1860	END MODULE module_small_step_em

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

Download in other formats:

Original Format