Context Navigation

module_small_step_em.F

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

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

Download in other formats:

Original Format