Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_small_step_em.F @ 2759

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

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

Download in other formats:

Original Format