Context Navigation

module_big_step_utilities_em.F @ 2756

Last change on this file since 2756 was 1729, checked in by aslmd, 7 years ago
MESOSCALE VENUS. commented cp=f(T) to easily find where changes are. removed obsolete commented tests or prints.
File size: 163.2 KB

Line
1	!WRF:MODEL_LAYER:DYNAMICS
2	!
3
4	#if (RWORDSIZE == 4)
5	# define VPOWX vspowx
6	# define VPOW vspow
7	#else
8	# define VPOWX vpowx
9	# define VPOW vpow
10	#endif
11
12
13	MODULE module_big_step_utilities_em
14
15	USE module_domain
16	USE module_model_constants
17	USE module_state_description
18	USE module_configure
19	USE module_wrf_error
20
21	CONTAINS
22
23	!-------------------------------------------------------------------------------
24
25	SUBROUTINE calc_mu_uv ( config_flags, &
26	mu, mub, muu, muv, &
27	ids, ide, jds, jde, kds, kde, &
28	ims, ime, jms, jme, kms, kme, &
29	its, ite, jts, jte, kts, kte )
30
31	IMPLICIT NONE
32
33	! Input data
34
35	TYPE(grid_config_rec_type ) , INTENT(IN ) :: config_flags
36
37	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
38	ims, ime, jms, jme, kms, kme, &
39	its, ite, jts, jte, kts, kte
40
41	REAL, DIMENSION( ims:ime , jms:jme ) , INTENT( OUT) :: muu, muv
42	REAL, DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mu, mub
43
44	! local stuff
45
46	INTEGER :: i, j, itf, jtf, im, jm
47
48	!<DESCRIPTION>
49	!
50	! calc_mu_uv calculates the full column dry-air mass at the staggered
51	! horizontal velocity points (u,v) and places the results in muu and muv.
52	! This routine uses the reference state (mub) and perturbation state (mu)
53	!
54	!</DESCRIPTION>
55
56
57	itf=ite
58	jtf=MIN(jte,jde-1)
59
60	IF ( ( its .NE. ids ) .AND. ( ite .NE. ide ) ) THEN
61	DO j=jts,jtf
62	DO i=its,itf
63	muu(i,j) = 0.5*(mu(i,j)+mu(i-1,j)+mub(i,j)+mub(i-1,j))
64	ENDDO
65	ENDDO
66	ELSE IF ( ( its .EQ. ids ) .AND. ( ite .NE. ide ) ) THEN
67	DO j=jts,jtf
68	DO i=its+1,itf
69	muu(i,j) = 0.5*(mu(i,j)+mu(i-1,j)+mub(i,j)+mub(i-1,j))
70	ENDDO
71	ENDDO
72	i=its
73	im = its
74	if(config_flags%periodic_x) im = its-1
75	DO j=jts,jtf
76	! muu(i,j) = mu(i,j) +mub(i,j)
77	! fix for periodic b.c., 13 march 2004, wcs
78	muu(i,j) = 0.5*(mu(i,j)+mu(im,j)+mub(i,j)+mub(im,j))
79	ENDDO
80	ELSE IF ( ( its .NE. ids ) .AND. ( ite .EQ. ide ) ) THEN
81	DO j=jts,jtf
82	DO i=its,itf-1
83	muu(i,j) = 0.5*(mu(i,j)+mu(i-1,j)+mub(i,j)+mub(i-1,j))
84	ENDDO
85	ENDDO
86	i=ite
87	im = ite-1
88	if(config_flags%periodic_x) im = ite
89	DO j=jts,jtf
90	! muu(i,j) = mu(i-1,j) +mub(i-1,j)
91	! fix for periodic b.c., 13 march 2004, wcs
92	muu(i,j) = 0.5*(mu(i-1,j)+mu(im,j)+mub(i-1,j)+mub(im,j))
93	ENDDO
94	ELSE IF ( ( its .EQ. ids ) .AND. ( ite .EQ. ide ) ) THEN
95	DO j=jts,jtf
96	DO i=its+1,itf-1
97	muu(i,j) = 0.5*(mu(i,j)+mu(i-1,j)+mub(i,j)+mub(i-1,j))
98	ENDDO
99	ENDDO
100	i=its
101	im = its
102	if(config_flags%periodic_x) im = its-1
103	DO j=jts,jtf
104	! muu(i,j) = mu(i,j) +mub(i,j)
105	! fix for periodic b.c., 13 march 2004, wcs
106	muu(i,j) = 0.5*(mu(i,j)+mu(im,j)+mub(i,j)+mub(im,j))
107	ENDDO
108	i=ite
109	im = ite-1
110	if(config_flags%periodic_x) im = ite
111	DO j=jts,jtf
112	! muu(i,j) = mu(i-1,j) +mub(i-1,j)
113	! fix for periodic b.c., 13 march 2004, wcs
114	muu(i,j) = 0.5*(mu(i-1,j)+mu(im,j)+mub(i-1,j)+mub(im,j))
115	ENDDO
116	END IF
117
118	itf=MIN(ite,ide-1)
119	jtf=jte
120
121	IF ( ( jts .NE. jds ) .AND. ( jte .NE. jde ) ) THEN
122	DO j=jts,jtf
123	DO i=its,itf
124	muv(i,j) = 0.5*(mu(i,j)+mu(i,j-1)+mub(i,j)+mub(i,j-1))
125	ENDDO
126	ENDDO
127	ELSE IF ( ( jts .EQ. jds ) .AND. ( jte .NE. jde ) ) THEN
128	DO j=jts+1,jtf
129	DO i=its,itf
130	muv(i,j) = 0.5*(mu(i,j)+mu(i,j-1)+mub(i,j)+mub(i,j-1))
131	ENDDO
132	ENDDO
133	j=jts
134	jm = jts
135	if(config_flags%periodic_y) jm = jts-1
136	DO i=its,itf
137	! muv(i,j) = mu(i,j) +mub(i,j)
138	! fix for periodic b.c., 13 march 2004, wcs
139	muv(i,j) = 0.5*(mu(i,j)+mu(i,jm)+mub(i,j)+mub(i,jm))
140	ENDDO
141	ELSE IF ( ( jts .NE. jds ) .AND. ( jte .EQ. jde ) ) THEN
142	DO j=jts,jtf-1
143	DO i=its,itf
144	muv(i,j) = 0.5*(mu(i,j)+mu(i,j-1)+mub(i,j)+mub(i,j-1))
145	ENDDO
146	ENDDO
147	j=jte
148	jm = jte-1
149	if(config_flags%periodic_y) jm = jte
150	DO i=its,itf
151	muv(i,j) = mu(i,j-1) +mub(i,j-1)
152	! fix for periodic b.c., 13 march 2004, wcs
153	muv(i,j) = 0.5*(mu(i,j-1)+mu(i,jm)+mub(i,j-1)+mub(i,jm))
154	ENDDO
155	ELSE IF ( ( jts .EQ. jds ) .AND. ( jte .EQ. jde ) ) THEN
156	DO j=jts+1,jtf-1
157	DO i=its,itf
158	muv(i,j) = 0.5*(mu(i,j)+mu(i,j-1)+mub(i,j)+mub(i,j-1))
159	ENDDO
160	ENDDO
161	j=jts
162	jm = jts
163	if(config_flags%periodic_y) jm = jts-1
164	DO i=its,itf
165	! muv(i,j) = mu(i,j) +mub(i,j)
166	! fix for periodic b.c., 13 march 2004, wcs
167	muv(i,j) = 0.5*(mu(i,j)+mu(i,jm)+mub(i,j)+mub(i,jm))
168	ENDDO
169	j=jte
170	jm = jte-1
171	if(config_flags%periodic_y) jm = jte
172	DO i=its,itf
173	! muv(i,j) = mu(i,j-1) +mub(i,j-1)
174	! fix for periodic b.c., 13 march 2004, wcs
175	muv(i,j) = 0.5*(mu(i,j-1)+mu(i,jm)+mub(i,j-1)+mub(i,jm))
176	ENDDO
177	END IF
178
179	END SUBROUTINE calc_mu_uv
180
181	!-------------------------------------------------------------------------------
182
183	SUBROUTINE calc_mu_uv_1 ( config_flags, &
184	mu, muu, muv, &
185	ids, ide, jds, jde, kds, kde, &
186	ims, ime, jms, jme, kms, kme, &
187	its, ite, jts, jte, kts, kte )
188
189	IMPLICIT NONE
190
191	! Input data
192
193	TYPE(grid_config_rec_type ) , INTENT(IN ) :: config_flags
194
195	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
196	ims, ime, jms, jme, kms, kme, &
197	its, ite, jts, jte, kts, kte
198
199	REAL, DIMENSION( ims:ime , jms:jme ) , INTENT( OUT) :: muu, muv
200	REAL, DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mu
201
202	! local stuff
203
204	INTEGER :: i, j, itf, jtf, im, jm
205
206	!<DESCRIPTION>
207	!
208	! calc_mu_uv calculates the full column dry-air mass at the staggered
209	! horizontal velocity points (u,v) and places the results in muu and muv.
210	! This routine uses the full state (mu)
211	!
212	!</DESCRIPTION>
213
214	itf=ite
215	jtf=MIN(jte,jde-1)
216
217	IF ( ( its .NE. ids ) .AND. ( ite .NE. ide ) ) THEN
218	DO j=jts,jtf
219	DO i=its,itf
220	muu(i,j) = 0.5*(mu(i,j)+mu(i-1,j))
221	ENDDO
222	ENDDO
223	ELSE IF ( ( its .EQ. ids ) .AND. ( ite .NE. ide ) ) THEN
224	DO j=jts,jtf
225	DO i=its+1,itf
226	muu(i,j) = 0.5*(mu(i,j)+mu(i-1,j))
227	ENDDO
228	ENDDO
229	i=its
230	im = its
231	if(config_flags%periodic_x) im = its-1
232	DO j=jts,jtf
233	muu(i,j) = 0.5*(mu(i,j)+mu(im,j))
234	ENDDO
235	ELSE IF ( ( its .NE. ids ) .AND. ( ite .EQ. ide ) ) THEN
236	DO j=jts,jtf
237	DO i=its,itf-1
238	muu(i,j) = 0.5*(mu(i,j)+mu(i-1,j))
239	ENDDO
240	ENDDO
241	i=ite
242	im = ite-1
243	if(config_flags%periodic_x) im = ite
244	DO j=jts,jtf
245	muu(i,j) = 0.5*(mu(i-1,j)+mu(im,j))
246	ENDDO
247	ELSE IF ( ( its .EQ. ids ) .AND. ( ite .EQ. ide ) ) THEN
248	DO j=jts,jtf
249	DO i=its+1,itf-1
250	muu(i,j) = 0.5*(mu(i,j)+mu(i-1,j))
251	ENDDO
252	ENDDO
253	i=its
254	im = its
255	if(config_flags%periodic_x) im = its-1
256	DO j=jts,jtf
257	muu(i,j) = 0.5*(mu(i,j)+mu(im,j))
258	ENDDO
259	i=ite
260	im = ite-1
261	if(config_flags%periodic_x) im = ite
262	DO j=jts,jtf
263	muu(i,j) = 0.5*(mu(i-1,j)+mu(im,j))
264	ENDDO
265	END IF
266
267	itf=MIN(ite,ide-1)
268	jtf=jte
269
270	IF ( ( jts .NE. jds ) .AND. ( jte .NE. jde ) ) THEN
271	DO j=jts,jtf
272	DO i=its,itf
273	muv(i,j) = 0.5*(mu(i,j)+mu(i,j-1))
274	ENDDO
275	ENDDO
276	ELSE IF ( ( jts .EQ. jds ) .AND. ( jte .NE. jde ) ) THEN
277	DO j=jts+1,jtf
278	DO i=its,itf
279	muv(i,j) = 0.5*(mu(i,j)+mu(i,j-1))
280	ENDDO
281	ENDDO
282	j=jts
283	jm = jts
284	if(config_flags%periodic_y) jm = jts-1
285	DO i=its,itf
286	muv(i,j) = 0.5*(mu(i,j)+mu(i,jm))
287	ENDDO
288	ELSE IF ( ( jts .NE. jds ) .AND. ( jte .EQ. jde ) ) THEN
289	DO j=jts,jtf-1
290	DO i=its,itf
291	muv(i,j) = 0.5*(mu(i,j)+mu(i,j-1))
292	ENDDO
293	ENDDO
294	j=jte
295	jm = jte-1
296	if(config_flags%periodic_y) jm = jte
297	DO i=its,itf
298	muv(i,j) = 0.5*(mu(i,j-1)+mu(i,jm))
299	ENDDO
300	ELSE IF ( ( jts .EQ. jds ) .AND. ( jte .EQ. jde ) ) THEN
301	DO j=jts+1,jtf-1
302	DO i=its,itf
303	muv(i,j) = 0.5*(mu(i,j)+mu(i,j-1))
304	ENDDO
305	ENDDO
306	j=jts
307	jm = jts
308	if(config_flags%periodic_y) jm = jts-1
309	DO i=its,itf
310	muv(i,j) = 0.5*(mu(i,j)+mu(i,jm))
311	ENDDO
312	j=jte
313	jm = jte-1
314	if(config_flags%periodic_y) jm = jte
315	DO i=its,itf
316	muv(i,j) = 0.5*(mu(i,j-1)+mu(i,jm))
317	ENDDO
318	END IF
319
320	END SUBROUTINE calc_mu_uv_1
321
322	!-------------------------------------------------------------------------------
323
324	SUBROUTINE couple_momentum ( muu, ru, u, msfu, &
325	muv, rv, v, msfv, &
326	mut, rw, w, msft, &
327	ids, ide, jds, jde, kds, kde, &
328	ims, ime, jms, jme, kms, kme, &
329	its, ite, jts, jte, kts, kte )
330
331	IMPLICIT NONE
332
333	! Input data
334
335	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
336	ims, ime, jms, jme, kms, kme, &
337	its, ite, jts, jte, kts, kte
338
339	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT( OUT) :: ru, rv, rw
340
341	REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: muu, muv, mut
342	REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: msfu, msfv, msft
343
344	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: u, v, w
345
346	! Local data
347
348	INTEGER :: i, j, k, itf, jtf, ktf
349
350	!<DESCRIPTION>
351	!
352	! couple_momentum couples the velocities to the full column mass and
353	! the map factors.
354	!
355	!</DESCRIPTION>
356
357	ktf=MIN(kte,kde-1)
358
359	itf=ite
360	jtf=MIN(jte,jde-1)
361
362	DO j=jts,jtf
363	DO k=kts,ktf
364	DO i=its,itf
365	ru(i,k,j)=u(i,k,j)*muu(i,j)/msfu(i,j)
366	ENDDO
367	ENDDO
368	ENDDO
369
370	itf=MIN(ite,ide-1)
371	jtf=jte
372
373	DO j=jts,jtf
374	DO k=kts,ktf
375	DO i=its,itf
376	rv(i,k,j)=v(i,k,j)*muv(i,j)/msfv(i,j)
377	ENDDO
378	ENDDO
379	ENDDO
380
381	itf=MIN(ite,ide-1)
382	jtf=MIN(jte,jde-1)
383
384	DO j=jts,jtf
385	DO k=kts,kte
386	DO i=its,itf
387	rw(i,k,j)=w(i,k,j)*mut(i,j)/msft(i,j)
388	ENDDO
389	ENDDO
390	ENDDO
391
392	END SUBROUTINE couple_momentum
393
394	!-------------------------------------------------------------------
395
396	SUBROUTINE calc_mu_staggered ( mu, mub, muu, muv, &
397	ids, ide, jds, jde, kds, kde, &
398	ims, ime, jms, jme, kms, kme, &
399	its, ite, jts, jte, kts, kte )
400
401	IMPLICIT NONE
402
403	! Input data
404
405	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
406	ims, ime, jms, jme, kms, kme, &
407	its, ite, jts, jte, kts, kte
408
409	REAL, DIMENSION( ims:ime , jms:jme ) , INTENT( OUT) :: muu, muv
410	REAL, DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mu, mub
411
412	! local stuff
413
414	INTEGER :: i, j, itf, jtf
415
416	!<DESCRIPTION>
417	!
418	! calc_mu_staggered calculates the full dry air mass at the staggered
419	! velocity points (u,v).
420	!
421	!</DESCRIPTION>
422
423	itf=ite
424	jtf=MIN(jte,jde-1)
425
426	IF ( ( its .NE. ids ) .AND. ( ite .NE. ide ) ) THEN
427	DO j=jts,jtf
428	DO i=its,itf
429	muu(i,j) = 0.5*(mu(i,j)+mu(i-1,j)+mub(i,j)+mub(i-1,j))
430	ENDDO
431	ENDDO
432	ELSE IF ( ( its .EQ. ids ) .AND. ( ite .NE. ide ) ) THEN
433	DO j=jts,jtf
434	DO i=its+1,itf
435	muu(i,j) = 0.5*(mu(i,j)+mu(i-1,j)+mub(i,j)+mub(i-1,j))
436	ENDDO
437	ENDDO
438	i=its
439	DO j=jts,jtf
440	muu(i,j) = mu(i,j) +mub(i,j)
441	ENDDO
442	ELSE IF ( ( its .NE. ids ) .AND. ( ite .EQ. ide ) ) THEN
443	DO j=jts,jtf
444	DO i=its,itf-1
445	muu(i,j) = 0.5*(mu(i,j)+mu(i-1,j)+mub(i,j)+mub(i-1,j))
446	ENDDO
447	ENDDO
448	i=ite
449	DO j=jts,jtf
450	muu(i,j) = mu(i-1,j) +mub(i-1,j)
451	ENDDO
452	ELSE IF ( ( its .EQ. ids ) .AND. ( ite .EQ. ide ) ) THEN
453	DO j=jts,jtf
454	DO i=its+1,itf-1
455	muu(i,j) = 0.5*(mu(i,j)+mu(i-1,j)+mub(i,j)+mub(i-1,j))
456	ENDDO
457	ENDDO
458	i=its
459	DO j=jts,jtf
460	muu(i,j) = mu(i,j) +mub(i,j)
461	ENDDO
462	i=ite
463	DO j=jts,jtf
464	muu(i,j) = mu(i-1,j) +mub(i-1,j)
465	ENDDO
466	END IF
467
468	itf=MIN(ite,ide-1)
469	jtf=jte
470
471	IF ( ( jts .NE. jds ) .AND. ( jte .NE. jde ) ) THEN
472	DO j=jts,jtf
473	DO i=its,itf
474	muv(i,j) = 0.5*(mu(i,j)+mu(i,j-1)+mub(i,j)+mub(i,j-1))
475	ENDDO
476	ENDDO
477	ELSE IF ( ( jts .EQ. jds ) .AND. ( jte .NE. jde ) ) THEN
478	DO j=jts+1,jtf
479	DO i=its,itf
480	muv(i,j) = 0.5*(mu(i,j)+mu(i,j-1)+mub(i,j)+mub(i,j-1))
481	ENDDO
482	ENDDO
483	j=jts
484	DO i=its,itf
485	muv(i,j) = mu(i,j) +mub(i,j)
486	ENDDO
487	ELSE IF ( ( jts .NE. jds ) .AND. ( jte .EQ. jde ) ) THEN
488	DO j=jts,jtf-1
489	DO i=its,itf
490	muv(i,j) = 0.5*(mu(i,j)+mu(i,j-1)+mub(i,j)+mub(i,j-1))
491	ENDDO
492	ENDDO
493	j=jte
494	DO i=its,itf
495	muv(i,j) = mu(i,j-1) +mub(i,j-1)
496	ENDDO
497	ELSE IF ( ( jts .EQ. jds ) .AND. ( jte .EQ. jde ) ) THEN
498	DO j=jts+1,jtf-1
499	DO i=its,itf
500	muv(i,j) = 0.5*(mu(i,j)+mu(i,j-1)+mub(i,j)+mub(i,j-1))
501	ENDDO
502	ENDDO
503	j=jts
504	DO i=its,itf
505	muv(i,j) = mu(i,j) +mub(i,j)
506	ENDDO
507	j=jte
508	DO i=its,itf
509	muv(i,j) = mu(i,j-1) +mub(i,j-1)
510	ENDDO
511	END IF
512
513	END SUBROUTINE calc_mu_staggered
514
515	!-------------------------------------------------------------------------------
516
517	SUBROUTINE couple ( mu, mub, rfield, field, name, &
518	msf, &
519	ids, ide, jds, jde, kds, kde, &
520	ims, ime, jms, jme, kms, kme, &
521	its, ite, jts, jte, kts, kte )
522
523	IMPLICIT NONE
524
525	! Input data
526
527	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
528	ims, ime, jms, jme, kms, kme, &
529	its, ite, jts, jte, kts, kte
530
531	CHARACTER(LEN=1) , INTENT(IN ) :: name
532
533	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT( OUT) :: rfield
534
535	REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mu, mub, msf
536
537	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: field
538
539	! Local data
540
541	INTEGER :: i, j, k, itf, jtf, ktf
542	REAL , DIMENSION(ims:ime,jms:jme) :: muu , muv
543
544	!<DESCRIPTION>
545	!
546	! subroutine couple couples the input variable with the dry-air
547	! column mass (mu).
548	!
549	!</DESCRIPTION>
550
551
552	ktf=MIN(kte,kde-1)
553
554	IF (name .EQ. 'u')THEN
555
556	CALL calc_mu_staggered ( mu, mub, muu, muv, &
557	ids, ide, jds, jde, kds, kde, &
558	ims, ime, jms, jme, kms, kme, &
559	its, ite, jts, jte, kts, kte )
560
561	itf=ite
562	jtf=MIN(jte,jde-1)
563
564	DO j=jts,jtf
565	DO k=kts,ktf
566	DO i=its,itf
567	rfield(i,k,j)=field(i,k,j)*muu(i,j)/msf(i,j)
568	ENDDO
569	ENDDO
570	ENDDO
571
572	ELSE IF (name .EQ. 'v')THEN
573
574	CALL calc_mu_staggered ( mu, mub, muu, muv, &
575	ids, ide, jds, jde, kds, kde, &
576	ims, ime, jms, jme, kms, kme, &
577	its, ite, jts, jte, kts, kte )
578
579	itf=ite
580	itf=MIN(ite,ide-1)
581	jtf=jte
582
583	DO j=jts,jtf
584	DO k=kts,ktf
585	DO i=its,itf
586	rfield(i,k,j)=field(i,k,j)*muv(i,j)/msf(i,j)
587	ENDDO
588	ENDDO
589	ENDDO
590
591	ELSE IF (name .EQ. 'w')THEN
592	itf=MIN(ite,ide-1)
593	jtf=MIN(jte,jde-1)
594	DO j=jts,jtf
595	DO k=kts,kte
596	DO i=its,itf
597	rfield(i,k,j)=field(i,k,j)*(mu(i,j)+mub(i,j))/msf(i,j)
598	ENDDO
599	ENDDO
600	ENDDO
601
602	ELSE IF (name .EQ. 'h')THEN
603	itf=MIN(ite,ide-1)
604	jtf=MIN(jte,jde-1)
605	DO j=jts,jtf
606	DO k=kts,kte
607	DO i=its,itf
608	rfield(i,k,j)=field(i,k,j)*(mu(i,j)+mub(i,j))
609	ENDDO
610	ENDDO
611	ENDDO
612
613	ELSE
614	itf=MIN(ite,ide-1)
615	jtf=MIN(jte,jde-1)
616	DO j=jts,jtf
617	DO k=kts,ktf
618	DO i=its,itf
619	rfield(i,k,j)=field(i,k,j)*(mu(i,j)+mub(i,j))
620	ENDDO
621	ENDDO
622	ENDDO
623
624	ENDIF
625
626	END SUBROUTINE couple
627
628	!-----------------------------------------------------------------------
629
630	SUBROUTINE calc_ww ( mu, ru, rv, ww, &
631	rdx, rdy, msft, dnw, &
632	ids, ide, jds, jde, kds, kde, &
633	ims, ime, jms, jme, kms, kme, &
634	its, ite, jts, jte, kts, kte )
635
636	IMPLICIT NONE
637
638	! Input data
639
640
641	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
642	ims, ime, jms, jme, kms, kme, &
643	its, ite, jts, jte, kts, kte
644
645	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: ru, rv
646	REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mu, msft
647	REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: dnw
648
649	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT ) :: ww
650	REAL , INTENT(IN ) :: rdx, rdy
651
652	! Local data
653
654	INTEGER :: i, j, k, itf, jtf, ktf
655	REAL , DIMENSION( its:ite ) :: dmdt
656
657	!<DESCRIPTION>
658	!
659	! calc_ww calculates omega using the mass-coupled velocities muu, muv.
660	! The algorithm integrates the continuity equation through the column
661	! followed by a diagnosis of omega.
662	!
663	!</DESCRIPTION>
664
665
666	jtf=MIN(jte,jde-1)
667	ktf=MIN(kte,kde-1)
668	itf=MIN(ite,ide-1)
669
670	DO j=jts,jtf
671
672	DO i=its,ite
673	dmdt(i) = 0.
674	ww(i,1,j) = 0.
675	ww(i,kte,j) = 0.
676	ENDDO
677
678	!! DO k=kts,ktf+1
679
680	DO k=kts,ktf
681	DO i=its,itf
682
683	dmdt(i) = dmdt(i) + dnw(k)* ( rdx*(ru(i+1,k,j)-ru(i,k,j)) &
684	+rdy*(rv(i,k,j+1)-rv(i,k,j)) )
685
686	ENDDO
687	ENDDO
688
689	! DO K=2,NZ-1
690	! ww(K,I)=ww(K-1,I)-DNW(K-1)*
691	! & (DMDT+RDX( xmu(i )u(K,I )
692	! & -xmu(im1)*u(k,im1)) )
693	! END DO
694
695	DO k=2,ktf
696	DO i=its,itf
697
698	ww(i,k,j)=ww(i,k-1,j) &
699	- dnw(k-1)* ( dmdt(i) &
700	+rdx*(ru(i+1,k-1,j)-ru(i,k-1,j)) &
701	+rdy*(rv(i,k-1,j+1)-rv(i,k-1,j)) )
702	ENDDO
703	ENDDO
704	ENDDO
705
706	END SUBROUTINE calc_ww
707
708
709	!-------------------------------------------------------------------------------
710
711	SUBROUTINE calc_ww_cp ( u, v, mup, mub, ww, &
712	rdx, rdy, msft, msfu, msfv, dnw, &
713	ids, ide, jds, jde, kds, kde, &
714	ims, ime, jms, jme, kms, kme, &
715	its, ite, jts, jte, kts, kte )
716
717	IMPLICIT NONE
718
719	! Input data
720
721
722	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
723	ims, ime, jms, jme, kms, kme, &
724	its, ite, jts, jte, kts, kte
725
726	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: u, v
727	REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mup, mub, &
728	msft, msfu, msfv
729	REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: dnw
730
731	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT ) :: ww
732	REAL , INTENT(IN ) :: rdx, rdy
733
734	! Local data
735
736	INTEGER :: i, j, k, itf, jtf, ktf
737	REAL , DIMENSION( its:ite ) :: dmdt
738	REAL , DIMENSION( its:ite, kts:kte ) :: divv
739	REAL , DIMENSION( its:ite+1, jts:jte+1 ) :: muu, muv
740
741	!<DESCRIPTION>
742	!
743	! calc_ww calculates omega using the velocities (u,v) and the dry-air
744	! column mass (mup+mub).
745	! The algorithm integrates the continuity equation through the column
746	! followed by a diagnosis of omega.
747	!
748	!</DESCRIPTION>
749
750	!<DESCRIPTION>
751	!
752	! calc_ww_cp calculates omega using the velocities (u,v) and the
753	! column mass mu.
754	!
755	!</DESCRIPTION>
756
757	jtf=MIN(jte,jde-1)
758	ktf=MIN(kte,kde-1)
759	itf=MIN(ite,ide-1)
760
761	! mu coupled with the appropriate map factor
762
763	DO j=jts,jtf
764	DO i=its,min(ite+1,ide)
765	muu(i,j) = 0.5*(mup(i,j)+mub(i,j)+mup(i-1,j)+mub(i-1,j))/msfu(i,j)
766	ENDDO
767	ENDDO
768
769	DO j=jts,min(jte+1,jde)
770	DO i=its,itf
771	muv(i,j) = 0.5*(mup(i,j)+mub(i,j)+mup(i,j-1)+mub(i,j-1))/msfv(i,j)
772	ENDDO
773	ENDDO
774
775	DO j=jts,jtf
776
777	DO i=its,ite
778	dmdt(i) = 0.
779	ww(i,1,j) = 0.
780	ww(i,kte,j) = 0.
781	ENDDO
782
783	DO k=kts,ktf
784	DO i=its,itf
785
786	divv(i,k) = msft(i,j)dnw(k)( rdx(muu(i+1,j)u(i+1,k,j)-muu(i,j)*u(i,k,j)) &
787	+rdy(muv(i,j+1)v(i,k,j+1)-muv(i,j)*v(i,k,j)) )
788
789	! dmdt(i) = dmdt(i) + dnw(k)* ( rdx*(ru(i+1,k,j)-ru(i,k,j)) &
790	! +rdy*(rv(i,k,j+1)-rv(i,k,j)) )
791
792	dmdt(i) = dmdt(i) + divv(i,k)
793
794
795	ENDDO
796	ENDDO
797
798	DO k=2,ktf
799	DO i=its,itf
800
801	! ww(i,k,j)=ww(i,k-1,j) &
802	! - dnw(k-1)* ( dmdt(i) &
803	! +rdx*(ru(i+1,k-1,j)-ru(i,k-1,j)) &
804	! +rdy*(rv(i,k-1,j+1)-rv(i,k-1,j)) )
805
806	ww(i,k,j)=ww(i,k-1,j) - dnw(k-1)*dmdt(i) - divv(i,k-1)
807
808	ENDDO
809	ENDDO
810	ENDDO
811
812
813	END SUBROUTINE calc_ww_cp
814
815
816	!-------------------------------------------------------------------------------
817
818	SUBROUTINE calc_cq ( moist, cqu, cqv, cqw, n_moist, &
819	ids, ide, jds, jde, kds, kde, &
820	ims, ime, jms, jme, kms, kme, &
821	its, ite, jts, jte, kts, kte )
822
823	IMPLICIT NONE
824
825	! Input data
826
827	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
828	ims, ime, jms, jme, kms, kme, &
829	its, ite, jts, jte, kts, kte
830
831	INTEGER , INTENT(IN ) :: n_moist
832
833
834	REAL, DIMENSION( ims:ime, kms:kme , jms:jme , n_moist ), INTENT(IN ) :: moist
835
836	REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT( OUT) :: cqu, cqv, cqw
837
838	! Local stuff
839
840	REAL :: qtot
841
842	INTEGER :: i, j, k, itf, jtf, ktf, ispe
843
844	!<DESCRIPTION>
845	!
846	! calc_cq calculates moist coefficients for the momentum equations.
847	!
848	!</DESCRIPTION>
849
850	itf=ite
851	jtf=MIN(jte,jde-1)
852	ktf=MIN(kte,kde-1)
853
854	IF( n_moist >= PARAM_FIRST_SCALAR ) THEN
855
856	DO j=jts,jtf
857	DO k=kts,ktf
858	DO i=its,itf
859	qtot = 0.
860	!DEC$ loop count(3)
861	DO ispe=PARAM_FIRST_SCALAR,n_moist
862	qtot = qtot + moist(i,k,j,ispe) + moist(i-1,k,j,ispe)
863	ENDDO
864	! qtot = 0.5*( moist(i ,k,j,1)+moist(i ,k,j,2)+moist(i ,k,j,3)+ &
865	! & moist(i-1,k,j,1)+moist(i-1,k,j,2)+moist(i-1,k,j,3) )
866	! cqu(i,k,j) = 1./(1.+qtot)
867	cqu(i,k,j) = 1./(1.+0.5*qtot)
868	ENDDO
869	ENDDO
870	ENDDO
871
872	itf=MIN(ite,ide-1)
873	jtf=jte
874
875	DO j=jts,jtf
876	DO k=kts,ktf
877	DO i=its,itf
878	qtot = 0.
879	!DEC$ loop count(3)
880	DO ispe=PARAM_FIRST_SCALAR,n_moist
881	qtot = qtot + moist(i,k,j,ispe) + moist(i,k,j-1,ispe)
882	ENDDO
883	! qtot = 0.5*( moist(i,k,j ,1)+moist(i,k,j ,2)+moist(i,k,j ,3)+ &
884	! & moist(i,k,j-1,1)+moist(i,k,j-1,2)+moist(i,k,j-1,3) )
885	! cqv(i,k,j) = 1./(1.+qtot)
886	cqv(i,k,j) = 1./(1.+0.5*qtot)
887	ENDDO
888	ENDDO
889	ENDDO
890
891	itf=MIN(ite,ide-1)
892	jtf=MIN(jte,jde-1)
893	DO j=jts,jtf
894	DO k=kts+1,ktf
895	DO i=its,itf
896	qtot = 0.
897	!DEC$ loop count(3)
898	DO ispe=PARAM_FIRST_SCALAR,n_moist
899	qtot = qtot + moist(i,k,j,ispe) + moist(i,k-1,j,ispe)
900	ENDDO
901	! qtot = 0.5*( moist(i,k ,j,1)+moist(i,k ,j,2)+moist(i,k-1,j,3)+ &
902	! & moist(i,k-1,j,1)+moist(i,k-1,j,2)+moist(i,k ,j,3) )
903	! cqw(i,k,j) = qtot
904	cqw(i,k,j) = 0.5*qtot
905	ENDDO
906	ENDDO
907	ENDDO
908
909	ELSE
910
911	DO j=jts,jtf
912	DO k=kts,ktf
913	DO i=its,itf
914	cqu(i,k,j) = 1.
915	ENDDO
916	ENDDO
917	ENDDO
918
919	itf=MIN(ite,ide-1)
920	jtf=jte
921
922	DO j=jts,jtf
923	DO k=kts,ktf
924	DO i=its,itf
925	cqv(i,k,j) = 1.
926	ENDDO
927	ENDDO
928	ENDDO
929
930	itf=MIN(ite,ide-1)
931	jtf=MIN(jte,jde-1)
932	DO j=jts,jtf
933	DO k=kts+1,ktf
934	DO i=its,itf
935	cqw(i,k,j) = 0.
936	ENDDO
937	ENDDO
938	ENDDO
939
940	END IF
941
942	END SUBROUTINE calc_cq
943
944	!----------------------------------------------------------------------
945
946	SUBROUTINE calc_alt ( alt, al, alb, &
947	ids, ide, jds, jde, kds, kde, &
948	ims, ime, jms, jme, kms, kme, &
949	its, ite, jts, jte, kts, kte )
950
951	IMPLICIT NONE
952
953	! Input data
954
955	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
956	ims, ime, jms, jme, kms, kme, &
957	its, ite, jts, jte, kts, kte
958
959	REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(IN ) :: alb, al
960	REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT( OUT) :: alt
961
962	! Local stuff
963
964	INTEGER :: i, j, k, itf, jtf, ktf
965
966	!<DESCRIPTION>
967	!
968	! calc_alt computes the full inverse density
969	!
970	!</DESCRIPTION>
971
972	itf=MIN(ite,ide-1)
973	jtf=MIN(jte,jde-1)
974	ktf=MIN(kte,kde-1)
975
976	DO j=jts,jtf
977	DO k=kts,ktf
978	DO i=its,itf
979	alt(i,k,j) = al(i,k,j)+alb(i,k,j)
980	ENDDO
981	ENDDO
982	ENDDO
983
984
985	END SUBROUTINE calc_alt
986
987	!----------------------------------------------------------------------
988
989	SUBROUTINE calc_p_rho_phi ( moist, n_moist, &
990	al, alb, mu, muts, ph, p, pb, &
991	t, p0, t0, znu, dnw, rdnw, &
992	rdn, non_hydrostatic, &
993	ids, ide, jds, jde, kds, kde, &
994	ims, ime, jms, jme, kms, kme, &
995	its, ite, jts, jte, kts, kte )
996
997	IMPLICIT NONE
998
999	! Input data
1000
1001	LOGICAL , INTENT(IN ) :: non_hydrostatic
1002
1003	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
1004	ims, ime, jms, jme, kms, kme, &
1005	its, ite, jts, jte, kts, kte
1006
1007	INTEGER , INTENT(IN ) :: n_moist
1008
1009	REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(IN ) :: alb, &
1010	pb, &
1011	t
1012
1013	REAL, DIMENSION( ims:ime , kms:kme , jms:jme, n_moist ), INTENT(IN ) :: moist
1014
1015	REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT( OUT) :: al, p
1016
1017	REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(INOUT) :: ph
1018
1019	REAL, DIMENSION( ims:ime , jms:jme ), INTENT(IN ) :: mu, muts
1020
1021	REAL, DIMENSION( kms:kme ), INTENT(IN ) :: znu, dnw, rdnw, rdn
1022
1023	REAL, INTENT(IN ) :: t0, p0
1024
1025	! Local stuff
1026
1027	INTEGER :: i, j, k, itf, jtf, ktf, ispe
1028	REAL :: qvf, qtot, qf1, qf2
1029	REAL, DIMENSION( its:ite) :: temp,cpovcv_v
1030
1031
1032	!<DESCRIPTION>
1033	!
1034	! For the nonhydrostatic option, calc_p_rho_phi calculates the
1035	! diagnostic quantities pressure and (inverse) density from the
1036	! prognostic variables using the equation of state.
1037	!
1038	! For the hydrostatic option, calc_p_rho_phi calculates the
1039	! diagnostic quantities (inverse) density and geopotential from the
1040	! prognostic variables using the equation of state and the hydrostatic
1041	! equation.
1042	!
1043	!</DESCRIPTION>
1044
1045	itf=MIN(ite,ide-1)
1046	jtf=MIN(jte,jde-1)
1047	ktf=MIN(kte,kde-1)
1048
1049	#ifndef INTELMKL
1050	cpovcv_v = cpovcv
1051	#endif
1052
1053	IF (non_hydrostatic) THEN
1054
1055	IF (n_moist >= PARAM_FIRST_SCALAR ) THEN
1056
1057	DO j=jts,jtf
1058	DO k=kts,ktf
1059	DO i=its,itf
1060	qvf = 1.+rvovrd*moist(i,k,j,P_QV)
1061	al(i,k,j)=-1./muts(i,j)(alb(i,k,j)mu(i,j) &
1062	+rdnw(k)*(ph(i,k+1,j)-ph(i,k,j)))
1063	temp(i)=(r_d(t0+t(i,k,j))qvf)/ &
1064	(p0*(al(i,k,j)+alb(i,k,j)))
1065	ENDDO
1066	#ifdef INTELMKL
1067	CALL VPOWX ( itf-its+1, temp(its), cpovcv, p(its,k,j) )
1068	#else
1069	! use vector version from libmassv or from compat lib in frame/libmassv.F
1070	CALL VPOW ( p(its,k,j), temp(its), cpovcv_v(its), itf-its+1 )
1071	#endif
1072	DO i=its,itf
1073	p(i,k,j)= p(i,k,j)*p0-pb(i,k,j)
1074	ENDDO
1075	ENDDO
1076	ENDDO
1077
1078	ELSE
1079
1080	DO j=jts,jtf
1081	DO k=kts,ktf
1082	DO i=its,itf
1083	al(i,k,j)=-1./muts(i,j)(alb(i,k,j)mu(i,j) &
1084	+rdnw(k)*(ph(i,k+1,j)-ph(i,k,j)))
1085	p(i,k,j)=p0( (r_d(t0+t(i,k,j)))/ &
1086	(p0(al(i,k,j)+alb(i,k,j))) )*cpovcv &
1087	-pb(i,k,j)
1088	ENDDO
1089	ENDDO
1090	ENDDO
1091
1092	END IF
1093
1094	ELSE
1095
1096	! hydrostatic pressure, al, and ph1 calc; WCS, 5 sept 2001
1097
1098
1099	IF (n_moist >= PARAM_FIRST_SCALAR ) THEN
1100
1101	DO j=jts,jtf
1102
1103	k=ktf ! top layer
1104	DO i=its,itf
1105
1106	qtot = 0.
1107	DO ispe=PARAM_FIRST_SCALAR,n_moist
1108	qtot = qtot + moist(i,k,j,ispe)
1109	ENDDO
1110	qf2 = 1./(1.+qtot)
1111	qf1 = qtot*qf2
1112
1113	p(i,k,j) = - 0.5(mu(i,j)+qf1muts(i,j))/rdnw(k)/qf2
1114	qvf = 1.+rvovrd*moist(i,k,j,P_QV)
1115	al(i,k,j) = (r_d/p1000mb)(t(i,k,j)+t0)qvf* &
1116	(((p(i,k,j)+pb(i,k,j))/p1000mb)**cvpm) - alb(i,k,j)
1117
1118	ENDDO
1119
1120	DO k=ktf-1,kts,-1 ! remaining layers, integrate down
1121	DO i=its,itf
1122
1123	qtot = 0.
1124	DO ispe=PARAM_FIRST_SCALAR,n_moist
1125	qtot = qtot + 0.5*( moist(i,k ,j,ispe) + moist(i,k+1,j,ispe) )
1126	ENDDO
1127	qf2 = 1./(1.+qtot)
1128	qf1 = qtot*qf2
1129
1130	p(i,k,j) = p(i,k+1,j) - (mu(i,j) + qf1*muts(i,j))/qf2/rdn(k+1)
1131	qvf = 1.+rvovrd*moist(i,k,j,P_QV)
1132	al(i,k,j) = (r_d/p1000mb)(t(i,k,j)+t0)qvf* &
1133	(((p(i,k,j)+pb(i,k,j))/p1000mb)**cvpm) - alb(i,k,j)
1134	ENDDO
1135	ENDDO
1136
1137	DO k=2,ktf+1 ! integrate hydrostatic equation for geopotential
1138	DO i=its,itf
1139
1140	! ph(i,k,j) = ph(i,k-1,j) - (1./rdnw(k-1))*( &
1141	! (muts(i,j)+mu(i,j))*al(i,k-1,j)+ &
1142	! mu(i,j)*alb(i,k-1,j) )
1143	ph(i,k,j) = ph(i,k-1,j) - (dnw(k-1))*( &
1144	(muts(i,j))*al(i,k-1,j)+ &
1145	mu(i,j)*alb(i,k-1,j) )
1146
1147
1148	ENDDO
1149	ENDDO
1150
1151	ENDDO
1152
1153	ELSE
1154
1155	DO j=jts,jtf
1156
1157	k=ktf ! top layer
1158	DO i=its,itf
1159
1160	qtot = 0.
1161	qf2 = 1./(1.+qtot)
1162	qf1 = qtot*qf2
1163
1164	p(i,k,j) = - 0.5(mu(i,j)+qf1muts(i,j))/rdnw(k)/qf2
1165	qvf = 1.
1166	al(i,k,j) = (r_d/p1000mb)(t(i,k,j)+t0)qvf* &
1167	(((p(i,k,j)+pb(i,k,j))/p1000mb)**cvpm) - alb(i,k,j)
1168
1169	ENDDO
1170
1171	DO k=ktf-1,kts,-1 ! remaining layers, integrate down
1172	DO i=its,itf
1173
1174	qtot = 0.
1175	qf2 = 1./(1.+qtot)
1176	qf1 = qtot*qf2
1177
1178	p(i,k,j) = p(i,k+1,j) - (mu(i,j) + qf1*muts(i,j))/qf2/rdn(k+1)
1179	qvf = 1.
1180	al(i,k,j) = (r_d/p1000mb)(t(i,k,j)+t0)qvf* &
1181	(((p(i,k,j)+pb(i,k,j))/p1000mb)**cvpm) - alb(i,k,j)
1182	ENDDO
1183	ENDDO
1184
1185	DO k=2,ktf+1 ! integrate hydrostatic equation for geopotential
1186	DO i=its,itf
1187
1188	! ph(i,k,j) = ph(i,k-1,j) - (1./rdnw(k-1))*( &
1189	! (muts(i,j)+mu(i,j))*al(i,k-1,j)+ &
1190	! mu(i,j)*alb(i,k-1,j) )
1191	ph(i,k,j) = ph(i,k-1,j) - (dnw(k-1))*( &
1192	(muts(i,j))*al(i,k-1,j)+ &
1193	mu(i,j)*alb(i,k-1,j) )
1194
1195
1196	ENDDO
1197	ENDDO
1198
1199	ENDDO
1200
1201	END IF
1202
1203	END IF
1204
1205	END SUBROUTINE calc_p_rho_phi
1206
1207	!----------------------------------------------------------------------
1208
1209	SUBROUTINE calc_php ( php, ph, phb, &
1210	ids, ide, jds, jde, kds, kde, &
1211	ims, ime, jms, jme, kms, kme, &
1212	its, ite, jts, jte, kts, kte )
1213
1214	IMPLICIT NONE
1215
1216	! Input data
1217
1218	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
1219	ims, ime, jms, jme, kms, kme, &
1220	its, ite, jts, jte, kts, kte
1221
1222	REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(IN ) :: phb, ph
1223	REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT( OUT) :: php
1224
1225	! Local stuff
1226
1227	INTEGER :: i, j, k, itf, jtf, ktf
1228
1229	!<DESCRIPTION>
1230	!
1231	! calc_php calculates the full geopotential from the reference state
1232	! geopotential and the perturbation geopotential (phb_ph).
1233	!
1234	!</DESCRIPTION>
1235
1236	itf=MIN(ite,ide-1)
1237	jtf=MIN(jte,jde-1)
1238	ktf=MIN(kte,kde-1)
1239
1240	DO j=jts,jtf
1241	DO k=kts,ktf
1242	DO i=its,itf
1243	php(i,k,j) = 0.5*(phb(i,k,j)+phb(i,k+1,j)+ph(i,k,j)+ph(i,k+1,j))
1244	ENDDO
1245	ENDDO
1246	ENDDO
1247
1248	END SUBROUTINE calc_php
1249
1250	!-------------------------------------------------------------------------------
1251
1252	SUBROUTINE diagnose_w( ph_tend, ph_new, ph_old, w, mu, dt, &
1253	u, v, ht, &
1254	cf1, cf2, cf3, rdx, rdy, msft, &
1255	ids, ide, jds, jde, kds, kde, &
1256	ims, ime, jms, jme, kms, kme, &
1257	its, ite, jts, jte, kts, kte )
1258
1259	IMPLICIT NONE
1260
1261	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
1262	ims, ime, jms, jme, kms, kme, &
1263	its, ite, jts, jte, kts, kte
1264
1265	REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(IN ) :: ph_tend, &
1266	ph_new, &
1267	ph_old, &
1268	u, &
1269	v
1270
1271
1272	REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT( OUT) :: w
1273
1274	REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN ) :: mu, ht, msft
1275
1276	REAL, INTENT(IN ) :: dt, cf1, cf2, cf3, rdx, rdy
1277
1278	INTEGER :: i, j, k, itf, jtf
1279
1280	itf=MIN(ite,ide-1)
1281	jtf=MIN(jte,jde-1)
1282
1283	!<DESCRIPTION>
1284	!
1285	! diagnose_w diagnoses the vertical velocity from the geopoential equation.
1286	! Used with the hydrostatic option.
1287	!
1288	!</DESCRIPTION>
1289
1290	DO j = jts, jtf
1291
1292	! lower b.c. on w
1293
1294	DO i = its, itf
1295	w(i,1,j)= msft(i,j)*( &
1296	.5rdy( &
1297	(ht(i,j+1)-ht(i,j )) &
1298	(cf1v(i,1,j+1)+cf2v(i,2,j+1)+cf3v(i,3,j+1)) &
1299	+(ht(i,j )-ht(i,j-1)) &
1300	(cf1v(i,1,j )+cf2v(i,2,j )+cf3v(i,3,j )) ) &
1301	+.5rdx( &
1302	(ht(i+1,j)-ht(i,j )) &
1303	(cf1u(i+1,1,j)+cf2u(i+1,2,j)+cf3u(i+1,3,j)) &
1304	+(ht(i,j )-ht(i-1,j)) &
1305	(cf1u(i ,1,j)+cf2u(i ,2,j)+cf3u(i ,3,j)) ) &
1306	)
1307	ENDDO
1308
1309	! use geopotential equation to diagnose w
1310
1311	DO k = 2, kte
1312	DO i = its, itf
1313	w(i,k,j) = msft(i,j)*( (ph_new(i,k,j)-ph_old(i,k,j))/dt &
1314	- ph_tend(i,k,j)/mu(i,j) )/g
1315
1316	ENDDO
1317	ENDDO
1318
1319	ENDDO
1320
1321	END SUBROUTINE diagnose_w
1322
1323	!-------------------------------------------------------------------------------
1324
1325	SUBROUTINE rhs_ph( ph_tend, u, v, ww, &
1326	ph, ph_old, phb, w, &
1327	mut, muu, muv, &
1328	fnm, fnp, &
1329	rdnw, cfn, cfn1, rdx, rdy, msft, &
1330	non_hydrostatic, &
1331	config_flags, &
1332	ids, ide, jds, jde, kds, kde, &
1333	ims, ime, jms, jme, kms, kme, &
1334	its, ite, jts, jte, kts, kte )
1335	IMPLICIT NONE
1336
1337	TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
1338
1339	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
1340	ims, ime, jms, jme, kms, kme, &
1341	its, ite, jts, jte, kts, kte
1342
1343	REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(IN ) :: &
1344	u, &
1345	v, &
1346	ww, &
1347	ph, &
1348	ph_old, &
1349	phb, &
1350	w
1351
1352	! pjj/cray
1353	! REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT( OUT) :: ph_tend
1354	REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(INOUT) :: ph_tend
1355
1356	REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN ) :: muu, muv, mut, msft
1357
1358	REAL, DIMENSION( kms:kme ), INTENT(IN ) :: rdnw, fnm, fnp
1359
1360	REAL, INTENT(IN ) :: cfn, cfn1, rdx, rdy
1361
1362	LOGICAL, INTENT(IN ) :: non_hydrostatic
1363
1364	! Local stuff
1365
1366	INTEGER :: i, j, k, itf, jtf, ktf, kz, i_start, j_start
1367	REAL :: ur, ul, ub, vr, vl, vb
1368	REAL, DIMENSION(its:ite,kts:kte) :: wdwn
1369
1370	INTEGER :: advective_order
1371
1372	LOGICAL :: specified
1373
1374	!<DESCRIPTION>
1375	!
1376	! rhs_ph calculates the large-timestep tendency terms for the geopotential
1377	! equation. These terms include the advection and "gw". The geopotential
1378	! equation is cast in advective form, so we don't use the flux form advection
1379	! algorithms here.
1380	!
1381	!</DESCRIPTION>
1382
1383	specified = .false.
1384	if(config_flags%specified .or. config_flags%nested) specified = .true.
1385
1386	advective_order = config_flags%h_sca_adv_order
1387	! advective_order = 2 ! original configuration (pre Oct 2001)
1388
1389	itf=MIN(ite,ide-1)
1390	jtf=MIN(jte,jde-1)
1391	ktf=MIN(kte,kde-1)
1392
1393	! advective form for the geopotential equation
1394
1395	DO j = jts, jtf
1396
1397	DO k = 2, kte
1398	DO i = its, itf
1399	wdwn(i,k) = .5(ww(i,k,j)+ww(i,k-1,j))rdnw(k-1) &
1400	*(ph(i,k,j)-ph(i,k-1,j)+phb(i,k,j)-phb(i,k-1,j))
1401	ENDDO
1402	ENDDO
1403
1404	DO k = 2, kte-1
1405	DO i = its, itf
1406	ph_tend(i,k,j) = ph_tend(i,k,j) &
1407	- (fnm(k)wdwn(i,k+1)+fnp(k)wdwn(i,k))
1408	ENDDO
1409	ENDDO
1410
1411	ENDDO
1412
1413	IF (non_hydrostatic) THEN ! add in "gw" term.
1414	DO j = jts, jtf ! in hydrostatic mode, "gw" will be diagnosed
1415	! after the timestep to give us "w"
1416	DO i = its, itf
1417	ph_tend(i,kde,j) = 0.
1418	ENDDO
1419
1420	DO k = 2, kte
1421	DO i = its, itf
1422	ph_tend(i,k,j) = ph_tend(i,k,j) + mut(i,j)gw(i,k,j)/msft(i,j)
1423	ENDDO
1424	ENDDO
1425
1426	ENDDO
1427
1428	END IF
1429
1430	IF (advective_order <= 2) THEN
1431
1432	! y (v) advection
1433
1434	i_start = its
1435	j_start = jts
1436	itf=MIN(ite,ide-1)
1437	jtf=MIN(jte,jde-1)
1438
1439	IF ( (config_flags%open_ys) .and. jts == jds ) j_start = jts+1
1440	IF ( (config_flags%open_ye) .and. jte == jde ) jtf = jtf-1
1441
1442	DO j = j_start, jtf
1443
1444	DO k = 2, kte-1
1445	DO i = i_start, itf
1446	ph_tend(i,k,j)=ph_tend(i,k,j) - .25rdy &
1447	( muv(i,j+1)(v(i,k,j+1)+v(i,k-1,j+1)) &
1448	(phb(i,k,j+1)-phb(i,k,j )+ph(i,k,j+1)-ph(i,k,j )) &
1449	+muv(i,j )(v(i,k,j )+v(i,k-1,j )) &
1450	(phb(i,k,j )-phb(i,k,j-1)+ph(i,k,j )-ph(i,k,j-1)) )
1451	ENDDO
1452	ENDDO
1453
1454	k = kte
1455	DO i = i_start, itf
1456	ph_tend(i,k,j)=ph_tend(i,k,j) - .5rdy &
1457	( muv(i,j+1)(cfnv(i,k-1,j+1)+cfn1v(i,k-2,j+1)) &
1458	(phb(i,k,j+1)-phb(i,k,j )+ph(i,k,j+1)-ph(i,k,j )) &
1459	+muv(i,j )(cfnv(i,k-1,j )+cfn1v(i,k-2,j )) &
1460	(phb(i,k,j )-phb(i,k,j-1)+ph(i,k,j )-ph(i,k,j-1)) )
1461	ENDDO
1462
1463	ENDDO
1464
1465	! x (u) advection
1466
1467	i_start = its
1468	j_start = jts
1469	itf=MIN(ite,ide-1)
1470	jtf=MIN(jte,jde-1)
1471
1472	IF ( (config_flags%open_xs) .and. its == ids ) i_start = its+1
1473	IF ( (config_flags%open_xe) .and. ite == ide ) itf = itf-1
1474
1475	DO j = j_start, jtf
1476
1477	DO k = 2, kte-1
1478	DO i = i_start, itf
1479	ph_tend(i,k,j)=ph_tend(i,k,j) - .25rdx &
1480	( muu(i+1,j)(u(i+1,k,j)+u(i+1,k-1,j)) &
1481	(phb(i+1,k,j)-phb(i ,k,j)+ph(i+1,k,j)-ph(i ,k,j)) &
1482	+muu(i ,j)(u(i ,k,j)+u(i ,k-1,j)) &
1483	(phb(i ,k,j)-phb(i-1,k,j)+ph(i ,k,j)-ph(i-1,k,j)) )
1484	ENDDO
1485	ENDDO
1486
1487	k = kte
1488	DO i = i_start, itf
1489	ph_tend(i,k,j)=ph_tend(i,k,j) - .5rdx &
1490	( muu(i+1,j)(cfnu(i+1,k-1,j)+cfn1u(i+1,k-2,j)) &
1491	(phb(i+1,k,j)-phb(i ,k,j)+ph(i+1,k,j)-ph(i ,k,j)) &
1492	+muu(i ,j)(cfnu(i ,k-1,j)+cfn1u(i ,k-2,j)) &
1493	(phb(i ,k,j)-phb(i-1,k,j)+ph(i ,k,j)-ph(i-1,k,j)) )
1494	ENDDO
1495
1496	ENDDO
1497
1498	ELSE IF (advective_order <= 4) THEN
1499
1500	! y (v) advection
1501
1502	i_start = its
1503	j_start = jts
1504	itf=MIN(ite,ide-1)
1505	jtf=MIN(jte,jde-1)
1506
1507	IF ( (config_flags%open_ys) .and. jts == jds ) j_start = jts+1
1508	IF ( (config_flags%open_ye) .and. jte == jde ) jtf = jtf-1
1509
1510	DO j = j_start, jtf
1511
1512	DO k = 2, kte-1
1513	DO i = i_start, itf
1514	ph_tend(i,k,j)=ph_tend(i,k,j) - .25rdy ( &
1515	( muv(i,j+1)*(v(i,k,j+1)+v(i,k-1,j+1)) &
1516	+muv(i,j )(v(i,k,j )+v(i,k-1,j )) ) (1./12.)*( &
1517	8.*(ph(i,k,j+1)-ph(i,k,j-1)) &
1518	-(ph(i,k,j+2)-ph(i,k,j-2)) &
1519	+8.*(phb(i,k,j+1)-phb(i,k,j-1)) &
1520	-(phb(i,k,j+2)-phb(i,k,j-2)) ) )
1521
1522
1523	ENDDO
1524	ENDDO
1525
1526	k = kte
1527	DO i = i_start, itf
1528	ph_tend(i,k,j)=ph_tend(i,k,j) - .5rdy ( &
1529	( muv(i,j+1)(cfnv(i,k-1,j+1)+cfn1*v(i,k-2,j+1)) &
1530	+muv(i,j )(cfnv(i,k-1,j )+cfn1v(i,k-2,j )) ) (1./12.)*( &
1531	8.*(ph(i,k,j+1)-ph(i,k,j-1)) &
1532	-(ph(i,k,j+2)-ph(i,k,j-2)) &
1533	+8.*(phb(i,k,j+1)-phb(i,k,j-1)) &
1534	-(phb(i,k,j+2)-phb(i,k,j-2)) ) )
1535
1536	ENDDO
1537
1538	ENDDO
1539
1540
1541	! x (u) advection
1542
1543	i_start = its
1544	j_start = jts
1545	itf=MIN(ite,ide-1)
1546	jtf=MIN(jte,jde-1)
1547
1548	IF ( (config_flags%open_xs) .and. its == ids ) i_start = its+1
1549	IF ( (config_flags%open_xe) .and. ite == ide ) itf = itf-1
1550
1551	DO j = j_start, jtf
1552
1553	DO k = 2, kte-1
1554	DO i = i_start, itf
1555	ph_tend(i,k,j)=ph_tend(i,k,j) - .25rdx( &
1556	( muu(i+1,j)*(u(i+1,k,j)+u(i+1,k-1,j)) &
1557	+muu(i,j )(u(i,k,j )+u(i,k-1,j )) ) (1./12.)*( &
1558	8.*(ph(i+1,k,j)-ph(i-1,k,j)) &
1559	-(ph(i+2,k,j)-ph(i-2,k,j)) &
1560	+8.*(phb(i+1,k,j)-phb(i-1,k,j)) &
1561	-(phb(i+2,k,j)-phb(i-2,k,j)) ) )
1562	ENDDO
1563	ENDDO
1564
1565	k = kte
1566	DO i = i_start, itf
1567	ph_tend(i,k,j)=ph_tend(i,k,j) - .5rdx( &
1568	( muu(i+1,j)(cfnu(i+1,k-1,j)+cfn1*u(i+1,k-2,j)) &
1569	+muu(i,j )(cfnu(i ,k-1,j)+cfn1u(i,k-2,j)) ) (1./12.)*( &
1570	8.*(ph(i+1,k,j)-ph(i-1,k,j)) &
1571	-(ph(i+2,k,j)-ph(i-2,k,j)) &
1572	+8.*(phb(i+1,k,j)-phb(i-1,k,j)) &
1573	-(phb(i+2,k,j)-phb(i-2,k,j)) ) )
1574	ENDDO
1575
1576	ENDDO
1577
1578	ELSE IF (advective_order <= 6) THEN
1579
1580	! y (v) advection
1581
1582	i_start = its
1583	j_start = jts
1584	itf=MIN(ite,ide-1)
1585	jtf=MIN(jte,jde-1)
1586
1587	! IF ( (config_flags%open_ys) .and. jts == jds ) j_start = jts+1
1588	! IF ( (config_flags%open_ye) .and. jte == jde ) jtf = jtf-1
1589
1590	IF (config_flags%open_ys .or. specified ) j_start = max(jts,jds+2)
1591	IF (config_flags%open_ye .or. specified ) jtf = min(jtf,jde-3)
1592
1593	DO j = j_start, jtf
1594
1595	DO k = 2, kte-1
1596	DO i = i_start, itf
1597	ph_tend(i,k,j)=ph_tend(i,k,j) - .25rdy ( &
1598	( muv(i,j+1)*(v(i,k,j+1)+v(i,k-1,j+1)) &
1599	+muv(i,j )(v(i,k,j )+v(i,k-1,j )) ) (1./60.)*( &
1600	45.*(ph(i,k,j+1)-ph(i,k,j-1)) &
1601	-9.*(ph(i,k,j+2)-ph(i,k,j-2)) &
1602	+(ph(i,k,j+3)-ph(i,k,j-3)) &
1603	+45.*(phb(i,k,j+1)-phb(i,k,j-1)) &
1604	-9.*(phb(i,k,j+2)-phb(i,k,j-2)) &
1605	+(phb(i,k,j+3)-phb(i,k,j-3)) ) )
1606
1607
1608	ENDDO
1609	ENDDO
1610
1611	k = kte
1612	DO i = i_start, itf
1613	ph_tend(i,k,j)=ph_tend(i,k,j) - .5rdy ( &
1614	( muv(i,j+1)(cfnv(i,k-1,j+1)+cfn1*v(i,k-2,j+1)) &
1615	+muv(i,j )(cfnv(i,k-1,j )+cfn1v(i,k-2,j )) ) (1./60.)*( &
1616	45.*(ph(i,k,j+1)-ph(i,k,j-1)) &
1617	-9.*(ph(i,k,j+2)-ph(i,k,j-2)) &
1618	+(ph(i,k,j+3)-ph(i,k,j-3)) &
1619	+45.*(phb(i,k,j+1)-phb(i,k,j-1)) &
1620	-9.*(phb(i,k,j+2)-phb(i,k,j-2)) &
1621	+(phb(i,k,j+3)-phb(i,k,j-3)) ) )
1622
1623	ENDDO
1624
1625	ENDDO
1626
1627
1628	! pick up near boundary rows using 4th order stencil
1629	! (open bc copy only goes out to jds-1 and jde, hence 4rth is ok but 6th is too big)
1630
1631	IF ( (config_flags%open_ys) .and. jts <= jds+1 ) THEN
1632
1633	j = jds+1
1634	DO k = 2, kte-1
1635	DO i = i_start, itf
1636	ph_tend(i,k,j)=ph_tend(i,k,j) - .25rdy ( &
1637	( muv(i,j+1)*(v(i,k,j+1)+v(i,k-1,j+1)) &
1638	+muv(i,j )(v(i,k,j )+v(i,k-1,j )) ) (1./12.)*( &
1639	8.*(ph(i,k,j+1)-ph(i,k,j-1)) &
1640	-(ph(i,k,j+2)-ph(i,k,j-2)) &
1641	+8.*(phb(i,k,j+1)-phb(i,k,j-1)) &
1642	-(phb(i,k,j+2)-phb(i,k,j-2)) ) )
1643
1644
1645	ENDDO
1646	ENDDO
1647
1648	k = kte
1649	DO i = i_start, itf
1650	ph_tend(i,k,j)=ph_tend(i,k,j) - .5rdy ( &
1651	( muv(i,j+1)(cfnv(i,k-1,j+1)+cfn1*v(i,k-2,j+1)) &
1652	+muv(i,j )(cfnv(i,k-1,j )+cfn1v(i,k-2,j )) ) (1./12.)*( &
1653	8.*(ph(i,k,j+1)-ph(i,k,j-1)) &
1654	-(ph(i,k,j+2)-ph(i,k,j-2)) &
1655	+8.*(phb(i,k,j+1)-phb(i,k,j-1)) &
1656	-(phb(i,k,j+2)-phb(i,k,j-2)) ) )
1657
1658	ENDDO
1659
1660	END IF
1661
1662	IF ( (config_flags%open_ye) .and. jte >= jde-2 ) THEN
1663
1664	j = jde-2
1665	DO k = 2, kte-1
1666	DO i = i_start, itf
1667	ph_tend(i,k,j)=ph_tend(i,k,j) - .25rdy ( &
1668	( muv(i,j+1)*(v(i,k,j+1)+v(i,k-1,j+1)) &
1669	+muv(i,j )(v(i,k,j )+v(i,k-1,j )) ) (1./12.)*( &
1670	8.*(ph(i,k,j+1)-ph(i,k,j-1)) &
1671	-(ph(i,k,j+2)-ph(i,k,j-2)) &
1672	+8.*(phb(i,k,j+1)-phb(i,k,j-1)) &
1673	-(phb(i,k,j+2)-phb(i,k,j-2)) ) )
1674
1675
1676	ENDDO
1677	ENDDO
1678
1679	k = kte
1680	DO i = i_start, itf
1681	ph_tend(i,k,j)=ph_tend(i,k,j) - .5rdy ( &
1682	( muv(i,j+1)(cfnv(i,k-1,j+1)+cfn1*v(i,k-2,j+1)) &
1683	+muv(i,j )(cfnv(i,k-1,j )+cfn1v(i,k-2,j )) ) (1./12.)*( &
1684	8.*(ph(i,k,j+1)-ph(i,k,j-1)) &
1685	-(ph(i,k,j+2)-ph(i,k,j-2)) &
1686	+8.*(phb(i,k,j+1)-phb(i,k,j-1)) &
1687	-(phb(i,k,j+2)-phb(i,k,j-2)) ) )
1688
1689	ENDDO
1690
1691	END IF
1692
1693	! x (u) advection
1694
1695	i_start = its
1696	j_start = jts
1697	itf=MIN(ite,ide-1)
1698	jtf=MIN(jte,jde-1)
1699
1700	IF (config_flags%open_xs .or. specified ) i_start = max(its,ids+2)
1701	IF (config_flags%open_xe .or. specified ) itf = min(itf,ide-3)
1702	IF ( config_flags%periodic_x ) i_start = its
1703	IF ( config_flags%periodic_x ) itf=MIN(ite,ide-1)
1704
1705	DO j = j_start, jtf
1706
1707	DO k = 2, kte-1
1708	DO i = i_start, itf
1709	ph_tend(i,k,j)=ph_tend(i,k,j) - .25rdx( &
1710	( muu(i+1,j)*(u(i+1,k,j)+u(i+1,k-1,j)) &
1711	+muu(i,j )(u(i,k,j )+u(i,k-1,j )) ) (1./60.)*( &
1712	45.*(ph(i+1,k,j)-ph(i-1,k,j)) &
1713	-9.*(ph(i+2,k,j)-ph(i-2,k,j)) &
1714	+(ph(i+3,k,j)-ph(i-3,k,j)) &
1715	+45.*(phb(i+1,k,j)-phb(i-1,k,j)) &
1716	-9.*(phb(i+2,k,j)-phb(i-2,k,j)) &
1717	+(phb(i+3,k,j)-phb(i-3,k,j)) ) )
1718	ENDDO
1719	ENDDO
1720
1721	k = kte
1722	DO i = i_start, itf
1723	ph_tend(i,k,j)=ph_tend(i,k,j) - .5rdx( &
1724	( muu(i+1,j)(cfnu(i+1,k-1,j)+cfn1*u(i+1,k-2,j)) &
1725	+muu(i,j )(cfnu(i ,k-1,j)+cfn1u(i,k-2,j)) ) (1./60.)*( &
1726	45.*(ph(i+1,k,j)-ph(i-1,k,j)) &
1727	-9.*(ph(i+2,k,j)-ph(i-2,k,j)) &
1728	+(ph(i+3,k,j)-ph(i-3,k,j)) &
1729	+45.*(phb(i+1,k,j)-phb(i-1,k,j)) &
1730	-9.*(phb(i+2,k,j)-phb(i-2,k,j)) &
1731	+(phb(i+3,k,j)-phb(i-3,k,j)) ) )
1732	ENDDO
1733
1734	ENDDO
1735
1736	IF ( (config_flags%open_xs) .and. its <= ids+1 ) THEN
1737	i = ids + 1
1738	DO j = j_start, jtf
1739	DO k = 2, kte-1
1740	ph_tend(i,k,j)=ph_tend(i,k,j) - .25rdx( &
1741	( muu(i+1,j)*(u(i+1,k,j)+u(i+1,k-1,j)) &
1742	+muu(i,j )(u(i,k,j )+u(i,k-1,j )) ) (1./12.)*( &
1743	8.*(ph(i+1,k,j)-ph(i-1,k,j)) &
1744	-(ph(i+2,k,j)-ph(i-2,k,j)) &
1745	+8.*(phb(i+1,k,j)-phb(i-1,k,j)) &
1746	-(phb(i+2,k,j)-phb(i-2,k,j)) ) )
1747	ENDDO
1748	k = kte
1749	ph_tend(i,k,j)=ph_tend(i,k,j) - .5rdx( &
1750	( muu(i+1,j)(cfnu(i+1,k-1,j)+cfn1*u(i+1,k-2,j)) &
1751	+muu(i,j )(cfnu(i ,k-1,j)+cfn1u(i,k-2,j)) ) (1./12.)*( &
1752	8.*(ph(i+1,k,j)-ph(i-1,k,j)) &
1753	-(ph(i+2,k,j)-ph(i-2,k,j)) &
1754	+8.*(phb(i+1,k,j)-phb(i-1,k,j)) &
1755	-(phb(i+2,k,j)-phb(i-2,k,j)) ) )
1756
1757	ENDDO
1758	END IF
1759
1760	IF ( (config_flags%open_xe) .and. ite >= ide-2 ) THEN
1761	i = ide-2
1762	DO j = j_start, jtf
1763	DO k = 2, kte-1
1764	ph_tend(i,k,j)=ph_tend(i,k,j) - .25rdx( &
1765	( muu(i+1,j)*(u(i+1,k,j)+u(i+1,k-1,j)) &
1766	+muu(i,j )(u(i,k,j )+u(i,k-1,j )) ) (1./12.)*( &
1767	8.*(ph(i+1,k,j)-ph(i-1,k,j)) &
1768	-(ph(i+2,k,j)-ph(i-2,k,j)) &
1769	+8.*(phb(i+1,k,j)-phb(i-1,k,j)) &
1770	-(phb(i+2,k,j)-phb(i-2,k,j)) ) )
1771	ENDDO
1772	k = kte
1773	ph_tend(i,k,j)=ph_tend(i,k,j) - .5rdx( &
1774	( muu(i+1,j)(cfnu(i+1,k-1,j)+cfn1*u(i+1,k-2,j)) &
1775	+muu(i,j )(cfnu(i ,k-1,j)+cfn1u(i,k-2,j)) ) (1./12.)*( &
1776	8.*(ph(i+1,k,j)-ph(i-1,k,j)) &
1777	-(ph(i+2,k,j)-ph(i-2,k,j)) &
1778	+8.*(phb(i+1,k,j)-phb(i-1,k,j)) &
1779	-(phb(i+2,k,j)-phb(i-2,k,j)) ) )
1780
1781	ENDDO
1782	END IF
1783
1784	END IF
1785
1786	! lateral open boundary conditions,
1787	! start with north and south (y) boundaries
1788
1789	i_start = its
1790	itf=MIN(ite,ide-1)
1791
1792	! south
1793
1794	IF ( (config_flags%open_ys) .and. jts == jds ) THEN
1795
1796	j=jts
1797
1798	DO k=2,kde
1799	kz = min(k,kde-1)
1800	DO i = its,itf
1801	vb =.5( fnm(kz)(v(i,kz ,j+1)+v(i,kz ,j )) &
1802	+fnp(kz)*(v(i,kz-1,j+1)+v(i,kz-1,j )) )
1803	vl=amin1(vb,0.)
1804	ph_tend(i,k,j)=ph_tend(i,k,j)-rdymut(i,j)( &
1805	+vl*(ph_old(i,k,j+1)-ph_old(i,k,j)))
1806	ENDDO
1807	ENDDO
1808
1809	END IF
1810
1811	! north
1812
1813	IF ( (config_flags%open_ye) .and. jte == jde ) THEN
1814
1815	j=jte-1
1816
1817	DO k=2,kde
1818	kz = min(k,kde-1)
1819	DO i = its,itf
1820	vb=.5( fnm(kz)(v(i,kz ,j+1)+v(i,kz ,j)) &
1821	+fnp(kz)*(v(i,kz-1,j+1)+v(i,kz-1,j)) )
1822	vr=amax1(vb,0.)
1823	ph_tend(i,k,j)=ph_tend(i,k,j)-rdymut(i,j)( &
1824	+vr*(ph_old(i,k,j)-ph_old(i,k,j-1)))
1825	ENDDO
1826	ENDDO
1827
1828	END IF
1829
1830	! now the east and west (y) boundaries
1831
1832	j_start = its
1833	jtf=MIN(jte,jde-1)
1834
1835	! west
1836
1837	IF ( (config_flags%open_xs) .and. its == ids ) THEN
1838
1839	i=its
1840
1841	DO j = jts,jtf
1842	DO k=2,kde-1
1843	kz = k
1844	ub =.5( fnm(kz)(u(i+1,kz ,j)+u(i ,kz ,j)) &
1845	+fnp(kz)*(u(i+1,kz-1,j)+u(i ,kz-1,j)) )
1846	ul=amin1(ub,0.)
1847	ph_tend(i,k,j)=ph_tend(i,k,j)-rdxmut(i,j)( &
1848	+ul*(ph_old(i+1,k,j)-ph_old(i,k,j)))
1849	ENDDO
1850
1851	k = kde
1852	kz = k
1853	ub =.5( fnm(kz)(u(i+1,kz ,j)+u(i ,kz ,j)) &
1854	+fnp(kz)*(u(i+1,kz-1,j)+u(i ,kz-1,j)) )
1855	ul=amin1(ub,0.)
1856	ph_tend(i,k,j)=ph_tend(i,k,j)-rdxmut(i,j)( &
1857	+ul*(ph_old(i+1,k,j)-ph_old(i,k,j)))
1858	ENDDO
1859
1860	END IF
1861
1862	! east
1863
1864	IF ( (config_flags%open_xe) .and. ite == ide ) THEN
1865
1866	i = ite-1
1867
1868	DO j = jts,jtf
1869	DO k=2,kde-1
1870	kz = k
1871	ub=.5( fnm(kz)(u(i+1,kz ,j)+u(i,kz ,j)) &
1872	+fnp(kz)*(u(i+1,kz-1,j)+u(i,kz-1,j)) )
1873	ur=amax1(ub,0.)
1874	ph_tend(i,k,j)=ph_tend(i,k,j)-rdxmut(i,j)( &
1875	+ur*(ph_old(i,k,j)-ph_old(i-1,k,j)))
1876	ENDDO
1877
1878	k = kde
1879	kz = k-1
1880	ub=.5( fnm(kz)(u(i+1,kz ,j)+u(i,kz ,j)) &
1881	+fnp(kz)*(u(i+1,kz-1,j)+u(i,kz-1,j)) )
1882	ur=amax1(ub,0.)
1883	ph_tend(i,k,j)=ph_tend(i,k,j)-rdxmut(i,j)( &
1884	+ur*(ph_old(i,k,j)-ph_old(i-1,k,j)))
1885
1886	ENDDO
1887
1888	END IF
1889
1890	END SUBROUTINE rhs_ph
1891
1892	!-------------------------------------------------------------------------------
1893
1894	SUBROUTINE horizontal_pressure_gradient( ru_tend,rv_tend, &
1895	ph,alt,p,pb,al,php,cqu,cqv, &
1896	muu,muv,mu,fnm,fnp,rdnw, &
1897	cf1,cf2,cf3,rdx,rdy,msft, &
1898	config_flags, non_hydrostatic, &
1899	ids, ide, jds, jde, kds, kde, &
1900	ims, ime, jms, jme, kms, kme, &
1901	its, ite, jts, jte, kts, kte )
1902
1903	IMPLICIT NONE
1904
1905	! Input data
1906
1907
1908	TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
1909
1910	LOGICAL, INTENT (IN ) :: non_hydrostatic
1911
1912	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
1913	ims, ime, jms, jme, kms, kme, &
1914	its, ite, jts, jte, kts, kte
1915
1916	REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(IN ) :: &
1917	ph, &
1918	alt, &
1919	al, &
1920	p, &
1921	pb, &
1922	php, &
1923	cqu, &
1924	cqv
1925
1926
1927	REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(INOUT) :: &
1928	ru_tend, &
1929	rv_tend
1930
1931	REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN ) :: muu, muv, mu, msft
1932
1933	REAL, DIMENSION( kms:kme ), INTENT(IN ) :: rdnw, fnm, fnp
1934
1935	REAL, INTENT(IN ) :: rdx, rdy, cf1, cf2, cf3
1936
1937	INTEGER :: i,j,k, itf, jtf, ktf, i_start, j_start
1938	REAL, DIMENSION( ims:ime, kms:kme ) :: dpn
1939	REAL :: dpx, dpy
1940
1941	LOGICAL :: specified
1942
1943	!<DESCRIPTION>
1944	!
1945	! horizontal_pressure_gradient calculates the
1946	! horizontal pressure gradient terms for the large-timestep tendency
1947	! in the horizontal momentum equations (u,v).
1948	!
1949	!</DESCRIPTION>
1950
1951	specified = .false.
1952	if(config_flags%specified .or. config_flags%nested) specified = .true.
1953
1954	! start with the north-south (y) pressure gradient
1955
1956	itf=MIN(ite,ide-1)
1957	jtf=jte
1958	ktf=MIN(kte,kde-1)
1959	i_start = its
1960	j_start = jts
1961	IF ( (config_flags%open_ys .or. specified .or. &
1962	config_flags%nested ) .and. jts == jds ) j_start = jts+1
1963	IF ( (config_flags%open_ye .or. specified .or. &
1964	config_flags%nested ) .and. jte == jde ) jtf = jtf-1
1965
1966	DO j = j_start, jtf
1967
1968	IF ( non_hydrostatic ) THEN
1969
1970	k=1
1971
1972	DO i = i_start, itf
1973	dpn(i,k) = .5( cf1(p(i,k ,j-1)+p(i,k ,j)) &
1974	+cf2*(p(i,k+1,j-1)+p(i,k+1,j)) &
1975	+cf3*(p(i,k+2,j-1)+p(i,k+2,j)) )
1976	dpn(i,kde) = 0.
1977	ENDDO
1978
1979	DO k=2,ktf
1980	DO i = i_start, itf
1981	dpn(i,k) = .5( fnm(k)(p(i,k ,j-1)+p(i,k ,j)) &
1982	+fnp(k)*(p(i,k-1,j-1)+p(i,k-1,j)) )
1983	END DO
1984	END DO
1985
1986	DO K=1,ktf
1987	DO i = i_start, itf
1988	dpy = .5rdymuv(i,j)*( &
1989	(ph (i,k+1,j)-ph (i,k+1,j-1) + ph(i,k,j)-ph(i,k,j-1)) &
1990	+(alt(i,k ,j)+alt(i,k ,j-1))*(p (i,k,j)-p (i,k,j-1)) &
1991	+(al (i,k ,j)+al (i,k ,j-1))*(pb(i,k,j)-pb(i,k,j-1)) )
1992	dpy = dpy + rdy(php(i,k,j)-php(i,k,j-1)) &
1993	(rdnw(k)(dpn(i,k+1)-dpn(i,k))-.5(mu(i,j-1)+mu(i,j)))
1994	rv_tend(i,k,j) = rv_tend(i,k,j)-cqv(i,k,j)*dpy
1995	END DO
1996	END DO
1997
1998	ELSE
1999
2000	DO K=1,ktf
2001	DO i = i_start, itf
2002	dpy = .5rdymuv(i,j)*( &
2003	(ph (i,k+1,j)-ph (i,k+1,j-1) + ph(i,k,j)-ph(i,k,j-1)) &
2004	+(alt(i,k ,j)+alt(i,k ,j-1))*(p (i,k,j)-p (i,k,j-1)) &
2005	+(al (i,k ,j)+al (i,k ,j-1))*(pb(i,k,j)-pb(i,k,j-1)) )
2006	rv_tend(i,k,j) = rv_tend(i,k,j)-cqv(i,k,j)*dpy
2007	END DO
2008	END DO
2009
2010	END IF
2011
2012	ENDDO
2013
2014	! now the east-west (x) pressure gradient
2015
2016	itf=ite
2017	jtf=MIN(jte,jde-1)
2018	ktf=MIN(kte,kde-1)
2019	i_start = its
2020	j_start = jts
2021	IF ( (config_flags%open_xs .or. specified .or. &
2022	config_flags%nested ) .and. its == ids ) i_start = its+1
2023	IF ( (config_flags%open_xe .or. specified .or. &
2024	config_flags%nested ) .and. ite == ide ) itf = itf-1
2025	IF ( config_flags%periodic_x ) i_start = its
2026	IF ( config_flags%periodic_x ) itf=ite
2027
2028	DO j = j_start, jtf
2029
2030	IF ( non_hydrostatic ) THEN
2031
2032	k=1
2033
2034	DO i = i_start, itf
2035	dpn(i,k) = .5( cf1(p(i-1,k ,j)+p(i,k ,j)) &
2036	+cf2*(p(i-1,k+1,j)+p(i,k+1,j)) &
2037	+cf3*(p(i-1,k+2,j)+p(i,k+2,j)) )
2038	dpn(i,kde) = 0.
2039	ENDDO
2040
2041	DO k=2,ktf
2042	DO i = i_start, itf
2043	dpn(i,k) = .5( fnm(k)(p(i-1,k ,j)+p(i,k ,j)) &
2044	+fnp(k)*(p(i-1,k-1,j)+p(i,k-1,j)) )
2045	END DO
2046	END DO
2047
2048	DO K=1,ktf
2049	DO i = i_start, itf
2050	dpx = .5rdxmuu(i,j)*( &
2051	(ph (i,k+1,j)-ph (i-1,k+1,j) + ph(i,k,j)-ph(i-1,k,j)) &
2052	+(alt(i,k ,j)+alt(i-1,k ,j))*(p (i,k,j)-p (i-1,k,j)) &
2053	+(al (i,k ,j)+al (i-1,k ,j))*(pb(i,k,j)-pb(i-1,k,j)) )
2054	dpx = dpx + rdx(php(i,k,j)-php(i-1,k,j)) &
2055	(rdnw(k)(dpn(i,k+1)-dpn(i,k))-.5(mu(i-1,j)+mu(i,j)))
2056	ru_tend(i,k,j) = ru_tend(i,k,j)-cqu(i,k,j)*dpx
2057	END DO
2058	END DO
2059
2060	ELSE
2061
2062	DO K=1,ktf
2063	DO i = i_start, itf
2064	dpx = .5rdxmuu(i,j)*( &
2065	(ph (i,k+1,j)-ph (i-1,k+1,j) + ph(i,k,j)-ph(i-1,k,j)) &
2066	+(alt(i,k ,j)+alt(i-1,k ,j))*(p (i,k,j)-p (i-1,k,j)) &
2067	+(al (i,k ,j)+al (i-1,k ,j))*(pb(i,k,j)-pb(i-1,k,j)) )
2068	ru_tend(i,k,j) = ru_tend(i,k,j)-cqu(i,k,j)*dpx
2069	END DO
2070	END DO
2071
2072	END IF
2073
2074	ENDDO
2075
2076	END SUBROUTINE horizontal_pressure_gradient
2077
2078	!-------------------------------------------------------------------------------
2079
2080	SUBROUTINE pg_buoy_w( rw_tend, p, cqw, mu, mub, &
2081	rdnw, rdn, g, msft, &
2082	ids, ide, jds, jde, kds, kde, &
2083	ims, ime, jms, jme, kms, kme, &
2084	its, ite, jts, jte, kts, kte )
2085
2086	IMPLICIT NONE
2087
2088	! Input data
2089
2090	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
2091	ims, ime, jms, jme, kms, kme, &
2092	its, ite, jts, jte, kts, kte
2093
2094	REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(IN ) :: p
2095	REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(INOUT) :: cqw
2096
2097
2098	REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(INOUT) :: rw_tend
2099
2100	REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN ) :: mub, mu, msft
2101
2102	REAL, DIMENSION( kms:kme ), INTENT(IN ) :: rdnw, rdn
2103
2104	REAL, INTENT(IN ) :: g
2105
2106	INTEGER :: itf, jtf, i, j, k
2107	REAL :: cq1, cq2
2108
2109
2110	!<DESCRIPTION>
2111	!
2112	! pg_buoy_w calculates the
2113	! vertical pressure gradient and buoyancy terms for the large-timestep
2114	! tendency in the vertical momentum equation.
2115	!
2116	!</DESCRIPTION>
2117
2118	! BUOYANCY AND PRESSURE GRADIENT TERM IN W EQUATION AT TIME T
2119
2120	itf=MIN(ite,ide-1)
2121	jtf=MIN(jte,jde-1)
2122
2123	DO j = jts,jtf
2124
2125	k=kde
2126	DO i=its,itf
2127	cq1 = 1./(1.+cqw(i,k-1,j))
2128	cq2 = cqw(i,k-1,j)*cq1
2129	rw_tend(i,k,j) = rw_tend(i,k,j)+(1./msft(i,j))g( &
2130	cq12.rdnw(k-1)*( -p(i,k-1,j)) &
2131	-mu(i,j)-cq2*mub(i,j) )
2132	END DO
2133
2134	DO k = 2, kde-1
2135	DO i = its,itf
2136	cq1 = 1./(1.+cqw(i,k,j))
2137	cq2 = cqw(i,k,j)*cq1
2138	cqw(i,k,j) = cq1
2139	rw_tend(i,k,j) = rw_tend(i,k,j)+(1./msft(i,j))g( &
2140	cq1rdn(k)(p(i,k,j)-p(i,k-1,j)) &
2141	-mu(i,j)-cq2*mub(i,j) )
2142	END DO
2143	ENDDO
2144
2145
2146	ENDDO
2147
2148	END SUBROUTINE pg_buoy_w
2149
2150	!-------------------------------------------------------------------------------
2151
2152	SUBROUTINE w_damp( rw_tend, ww, w, mut, rdnw, dt, &
2153	w_damping, &
2154	ids, ide, jds, jde, kds, kde, &
2155	ims, ime, jms, jme, kms, kme, &
2156	its, ite, jts, jte, kts, kte )
2157
2158	IMPLICIT NONE
2159
2160	! Input data
2161
2162	INTEGER , INTENT(IN ) :: w_damping
2163
2164	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
2165	ims, ime, jms, jme, kms, kme, &
2166	its, ite, jts, jte, kts, kte
2167
2168	REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(IN ) :: ww, w
2169
2170	REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(INOUT) :: rw_tend
2171
2172	REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN ) :: mut
2173
2174	REAL, DIMENSION( kms:kme ), INTENT(IN ) :: rdnw
2175
2176	REAL, INTENT(IN) :: dt
2177	REAL :: cfl, cf_n, cf_d, maxcfl, maxdub, maxdeta
2178
2179	INTEGER :: itf, jtf, i, j, k, maxi, maxj, maxk
2180	INTEGER :: some
2181	CHARACTER*512 :: temp
2182	CHARACTER (LEN=256) :: time_str
2183	CHARACTER (LEN=256) :: grid_str
2184
2185	!<DESCRIPTION>
2186	!
2187	! w_damp computes a damping term for the vertical velocity when the
2188	! vertical Courant number is too large. This was found to be preferable to
2189	! decreasing the timestep or increasing the diffusion in real-data applications
2190	! that produced potentially-unstable large vertical velocities because of
2191	! unphysically large heating rates coming from the cumulus parameterization
2192	! schemes run at moderately high resolutions (dx ~ O(10) km).
2193	!
2194	!</DESCRIPTION>
2195
2196	itf=MIN(ite,ide-1)
2197	jtf=MIN(jte,jde-1)
2198
2199	some = 0
2200	maxcfl = 0.
2201
2202	IF ( w_damping == 1 ) THEN
2203	DO j = jts,jtf
2204
2205	DO k = 2, kde-1
2206	DO i = its,itf
2207	#if 0
2208	cfl = abs(ww(i,k,j)/mut(i,j)rdnw(k)dt)
2209	if(cfl .gt. w_beta)then
2210	#else
2211	! restructure to get rid of divide
2212	cf_n = abs(ww(i,k,j)rdnw(k)dt)
2213	cf_d = abs(mut(i,j))
2214	if(cf_n .gt. cf_d*w_beta )then
2215	#endif
2216	cfl = abs(ww(i,k,j)/mut(i,j)rdnw(k)dt)
2217	IF ( cfl > maxcfl ) THEN
2218	maxcfl = cfl ; maxi = i ; maxj = j ; maxk = k
2219	maxdub = w(i,k,j) ; maxdeta = -1./rdnw(k)
2220	ENDIF
2221	WRITE(temp,*)i,j,k,' cfl,w,d(eta)=',cfl,w(i,k,j),-1./rdnw(k)
2222	CALL wrf_debug ( 100 , TRIM(temp) )
2223	if ( cfl > 2. ) some = some + 1
2224	rw_tend(i,k,j) = rw_tend(i,k,j)-sign(1.,w(i,k,j))w_alpha(cfl-w_beta)*mut(i,j)
2225	endif
2226	END DO
2227	ENDDO
2228	ENDDO
2229	ELSE
2230	! just print
2231	DO j = jts,jtf
2232
2233	DO k = 2, kde-1
2234	DO i = its,itf
2235	cf_n = abs(ww(i,k,j)rdnw(k)dt)
2236	cf_d = abs(mut(i,j))
2237	if(cf_n .gt. cf_d*w_beta )then
2238	cfl = abs(ww(i,k,j)/mut(i,j)rdnw(k)dt)
2239	IF ( cfl > maxcfl ) THEN
2240	maxcfl = cfl ; maxi = i ; maxj = j ; maxk = k
2241	maxdub = w(i,k,j) ; maxdeta = -1./rdnw(k)
2242	ENDIF
2243	WRITE(temp,*)i,j,k,' cfl,w,d(eta)=',cfl,w(i,k,j),-1./rdnw(k)
2244	CALL wrf_debug ( 100 , TRIM(temp) )
2245	if ( cfl > 2. ) some = some + 1
2246	endif
2247	END DO
2248	ENDDO
2249	ENDDO
2250	ENDIF
2251	IF ( some .GT. 0 ) THEN
2252	CALL get_current_time_string( time_str )
2253	CALL get_current_grid_name( grid_str )
2254	WRITE(wrf_err_message,*)some, &
2255	' points exceeded cfl=2 in domain '//TRIM(grid_str)//' at time '//TRIM(time_str)//' hours'
2256	CALL wrf_debug ( 0 , TRIM(wrf_err_message) )
2257	WRITE(wrf_err_message,*)'MAX AT i,j,k: ',maxi,maxj,maxk,' cfl,w,d(eta)=',maxcfl, &
2258	maxdub,maxdeta
2259	CALL wrf_debug ( 0 , TRIM(wrf_err_message) )
2260	ENDIF
2261
2262	END SUBROUTINE w_damp
2263
2264	!-------------------------------------------------------------------------------
2265
2266	SUBROUTINE horizontal_diffusion ( name, field, tendency, mu, &
2267	config_flags, &
2268	msfu, msfv, msft, khdif, xkmhd, rdx, rdy, &
2269	ids, ide, jds, jde, kds, kde, &
2270	ims, ime, jms, jme, kms, kme, &
2271	its, ite, jts, jte, kts, kte )
2272
2273	IMPLICIT NONE
2274
2275	! Input data
2276
2277	TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
2278
2279	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
2280	ims, ime, jms, jme, kms, kme, &
2281	its, ite, jts, jte, kts, kte
2282
2283	CHARACTER(LEN=1) , INTENT(IN ) :: name
2284
2285	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: field, xkmhd
2286
2287	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency
2288
2289	REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mu
2290
2291	REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: msfu, &
2292	msfv, &
2293	msft
2294
2295	REAL , INTENT(IN ) :: rdx, &
2296	rdy, &
2297	khdif
2298
2299	! Local data
2300
2301	INTEGER :: i, j, k, itf, jtf, ktf
2302
2303	INTEGER :: i_start, i_end, j_start, j_end
2304
2305	REAL :: mrdx, mkrdxm, mkrdxp, &
2306	mrdy, mkrdym, mkrdyp
2307	REAL :: pr_inv
2308
2309	LOGICAL :: specified
2310
2311	!<DESCRIPTION>
2312	!
2313	! horizontal_diffusion computes the horizontal diffusion tendency
2314	! on model horizontal coordinate surfaces.
2315	!
2316	!</DESCRIPTION>
2317
2318	pr_inv = 1./prandtl
2319	specified = .false.
2320	if(config_flags%specified .or. config_flags%nested) specified = .true.
2321
2322	ktf=MIN(kte,kde-1)
2323
2324	IF (name .EQ. 'u') THEN
2325
2326	i_start = its
2327	i_end = ite
2328	j_start = jts
2329	j_end = MIN(jte,jde-1)
2330
2331	IF ( config_flags%open_xs .or. specified ) i_start = MAX(ids+1,its)
2332	IF ( config_flags%open_xe .or. specified ) i_end = MIN(ide-1,ite)
2333	IF ( config_flags%open_ys .or. specified ) j_start = MAX(jds+1,jts)
2334	IF ( config_flags%open_ye .or. specified ) j_end = MIN(jde-2,jte)
2335	IF ( config_flags%periodic_x ) i_start = its
2336	IF ( config_flags%periodic_x ) i_end = ite
2337
2338
2339	DO j = j_start, j_end
2340	DO k=kts,ktf
2341	DO i = i_start, i_end
2342
2343	mkrdxm=msft(i-1,j)mu(i-1,j)xkmhd(i-1,k,j)*rdx
2344	mkrdxp=msft(i,j)mu(i,j)xkmhd(i,k,j)*rdx
2345	mrdx=msfu(i,j)*rdx
2346	mkrdym=0.5(msfu(i,j)+msfu(i,j-1)) &
2347	0.25(mu(i,j)+mu(i,j-1)+mu(i-1,j-1)+mu(i-1,j)) &
2348	0.25(xkmhd(i,k,j)+xkmhd(i,k,j-1)+xkmhd(i-1,k,j-1)+xkmhd(i-1,k,j))rdy
2349	mkrdyp=0.5(msfu(i,j)+msfu(i,j+1)) &
2350	0.25(mu(i,j)+mu(i,j+1)+mu(i-1,j+1)+mu(i-1,j)) &
2351	0.25(xkmhd(i,k,j)+xkmhd(i,k,j+1)+xkmhd(i-1,k,j+1)+xkmhd(i-1,k,j))rdy
2352	mrdy=msfu(i,j)*rdy
2353
2354	tendency(i,k,j)=tendency(i,k,j)+( &
2355	mrdx(mkrdxp(field(i+1,k,j)-field(i ,k,j)) &
2356	-mkrdxm*(field(i ,k,j)-field(i-1,k,j))) &
2357	+mrdy(mkrdyp(field(i,k,j+1)-field(i,k,j )) &
2358	-mkrdym*(field(i,k,j )-field(i,k,j-1))))
2359	ENDDO
2360	ENDDO
2361	ENDDO
2362
2363	ELSE IF (name .EQ. 'v')THEN
2364
2365	i_start = its
2366	i_end = MIN(ite,ide-1)
2367	j_start = jts
2368	j_end = jte
2369
2370	IF ( config_flags%open_xs .or. specified ) i_start = MAX(ids+1,its)
2371	IF ( config_flags%open_xe .or. specified ) i_end = MIN(ide-2,ite)
2372	IF ( config_flags%open_ys .or. specified ) j_start = MAX(jds+1,jts)
2373	IF ( config_flags%open_ye .or. specified ) j_end = MIN(jde-1,jte)
2374	IF ( config_flags%periodic_x ) i_start = its
2375	IF ( config_flags%periodic_x ) i_end = MIN(ite,ide-1)
2376
2377	DO j = j_start, j_end
2378	DO k=kts,ktf
2379	DO i = i_start, i_end
2380
2381	mkrdxm=0.5(msfv(i,j)+msfv(i-1,j)) &
2382	0.25(mu(i,j)+mu(i,j-1)+mu(i-1,j-1)+mu(i-1,j)) &
2383	0.25(xkmhd(i,k,j)+xkmhd(i,k,j-1)+xkmhd(i-1,k,j-1)+xkmhd(i-1,k,j))rdx
2384	mkrdxp=0.5(msfv(i,j)+msfv(i+1,j)) &
2385	0.25(mu(i,j)+mu(i,j-1)+mu(i+1,j-1)+mu(i+1,j)) &
2386	0.25(xkmhd(i,k,j)+xkmhd(i,k,j-1)+xkmhd(i+1,k,j-1)+xkmhd(i+1,k,j))rdx
2387	mrdx=msfv(i,j)*rdx
2388	mkrdym=msft(i,j-1)xkmhd(i,k,j-1)rdy
2389	mkrdyp=msft(i,j)xkmhd(i,k,j)rdy
2390	mrdy=msfv(i,j)*rdy
2391
2392	tendency(i,k,j)=tendency(i,k,j)+( &
2393	mrdx(mkrdxp(field(i+1,k,j)-field(i ,k,j)) &
2394	-mkrdxm*(field(i ,k,j)-field(i-1,k,j))) &
2395	+mrdy(mkrdyp(field(i,k,j+1)-field(i,k,j )) &
2396	-mkrdym*(field(i,k,j )-field(i,k,j-1))))
2397	ENDDO
2398	ENDDO
2399	ENDDO
2400
2401	ELSE IF (name .EQ. 'w')THEN
2402
2403	i_start = its
2404	i_end = MIN(ite,ide-1)
2405	j_start = jts
2406	j_end = MIN(jte,jde-1)
2407
2408	IF ( config_flags%open_xs .or. specified ) i_start = MAX(ids+1,its)
2409	IF ( config_flags%open_xe .or. specified ) i_end = MIN(ide-2,ite)
2410	IF ( config_flags%open_ys .or. specified ) j_start = MAX(jds+1,jts)
2411	IF ( config_flags%open_ye .or. specified ) j_end = MIN(jde-2,jte)
2412	IF ( config_flags%periodic_x ) i_start = its
2413	IF ( config_flags%periodic_x ) i_end = MIN(ite,ide-1)
2414
2415	DO j = j_start, j_end
2416	DO k=kts+1,ktf
2417	DO i = i_start, i_end
2418
2419	mkrdxm=msfu(i,j)* &
2420	0.25(mu(i,j)+mu(i-1,j)+mu(i,j)+mu(i-1,j)) &
2421	0.25(xkmhd(i,k,j)+xkmhd(i-1,k,j)+xkmhd(i,k-1,j)+xkmhd(i-1,k-1,j))rdx
2422	mkrdxp=msfu(i+1,j)* &
2423	0.25(mu(i+1,j)+mu(i,j)+mu(i+1,j)+mu(i,j)) &
2424	0.25(xkmhd(i+1,k,j)+xkmhd(i,k,j)+xkmhd(i+1,k-1,j)+xkmhd(i,k-1,j))rdx
2425	mrdx=msft(i,j)*rdx
2426	mkrdym=msfv(i,j)* &
2427	0.25(mu(i,j)+mu(i,j-1)+mu(i,j)+mu(i,j-1)) &
2428	0.25(xkmhd(i,k,j)+xkmhd(i,k,j-1)+xkmhd(i,k-1,j)+xkmhd(i,k-1,j-1))rdy
2429	mkrdyp=msfv(i,j+1)* &
2430	0.25(mu(i,j+1)+mu(i,j)+mu(i,j+1)+mu(i,j)) &
2431	0.25(xkmhd(i,k,j+1)+xkmhd(i,k,j)+xkmhd(i,k-1,j+1)+xkmhd(i,k-1,j))rdy
2432	mrdy=msft(i,j)*rdy
2433
2434	tendency(i,k,j)=tendency(i,k,j)+( &
2435	mrdx(mkrdxp(field(i+1,k,j)-field(i ,k,j)) &
2436	-mkrdxm*(field(i ,k,j)-field(i-1,k,j))) &
2437	+mrdy(mkrdyp(field(i,k,j+1)-field(i,k,j )) &
2438	-mkrdym*(field(i,k,j )-field(i,k,j-1))))
2439	ENDDO
2440	ENDDO
2441	ENDDO
2442
2443	ELSE
2444
2445
2446	i_start = its
2447	i_end = MIN(ite,ide-1)
2448	j_start = jts
2449	j_end = MIN(jte,jde-1)
2450
2451	IF ( config_flags%open_xs .or. specified ) i_start = MAX(ids+1,its)
2452	IF ( config_flags%open_xe .or. specified ) i_end = MIN(ide-2,ite)
2453	IF ( config_flags%open_ys .or. specified ) j_start = MAX(jds+1,jts)
2454	IF ( config_flags%open_ye .or. specified ) j_end = MIN(jde-2,jte)
2455	IF ( config_flags%periodic_x ) i_start = its
2456	IF ( config_flags%periodic_x ) i_end = MIN(ite,ide-1)
2457
2458	DO j = j_start, j_end
2459	DO k=kts,ktf
2460	DO i = i_start, i_end
2461
2462	mkrdxm=msfu(i,j)0.5(xkmhd(i,k,j)+xkmhd(i-1,k,j))0.5(mu(i,j)+mu(i-1,j))rdxpr_inv
2463	mkrdxp=msfu(i+1,j)0.5(xkmhd(i+1,k,j)+xkmhd(i,k,j))0.5(mu(i+1,j)+mu(i,j))rdxpr_inv
2464	mrdx=msft(i,j)*rdx
2465	mkrdym=msfv(i,j)0.5(xkmhd(i,k,j)+xkmhd(i,k,j-1))0.5(mu(i,j)+mu(i,j-1))rdypr_inv
2466	mkrdyp=msfv(i,j+1)0.5(xkmhd(i,k,j+1)+xkmhd(i,k,j))0.5(mu(i,j+1)+mu(i,j))rdypr_inv
2467	mrdy=msft(i,j)*rdy
2468
2469	tendency(i,k,j)=tendency(i,k,j)+( &
2470	mrdx(mkrdxp(field(i+1,k,j)-field(i ,k,j)) &
2471	-mkrdxm*(field(i ,k,j)-field(i-1,k,j))) &
2472	+mrdy(mkrdyp(field(i,k,j+1)-field(i,k,j )) &
2473	-mkrdym*(field(i,k,j )-field(i,k,j-1))))
2474	ENDDO
2475	ENDDO
2476	ENDDO
2477
2478	ENDIF
2479
2480	END SUBROUTINE horizontal_diffusion
2481
2482	!-----------------------------------------------------------------------------------------
2483
2484	SUBROUTINE horizontal_diffusion_3dmp ( name, field, tendency, mu, &
2485	config_flags, base_3d, &
2486	msfu, msfv, msft, khdif, xkmhd, rdx, rdy, &
2487	ids, ide, jds, jde, kds, kde, &
2488	ims, ime, jms, jme, kms, kme, &
2489	its, ite, jts, jte, kts, kte )
2490
2491	IMPLICIT NONE
2492
2493	! Input data
2494
2495	TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
2496
2497	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
2498	ims, ime, jms, jme, kms, kme, &
2499	its, ite, jts, jte, kts, kte
2500
2501	CHARACTER(LEN=1) , INTENT(IN ) :: name
2502
2503	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: field, &
2504	xkmhd, &
2505	base_3d
2506
2507	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency
2508
2509	REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mu
2510
2511	REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: msfu, &
2512	msfv, &
2513	msft
2514
2515	REAL , INTENT(IN ) :: rdx, &
2516	rdy, &
2517	khdif
2518
2519	! Local data
2520
2521	INTEGER :: i, j, k, itf, jtf, ktf
2522
2523	INTEGER :: i_start, i_end, j_start, j_end
2524
2525	REAL :: mrdx, mkrdxm, mkrdxp, &
2526	mrdy, mkrdym, mkrdyp
2527	REAL :: pr_inv
2528
2529	LOGICAL :: specified
2530
2531	!<DESCRIPTION>
2532	!
2533	! horizontal_diffusion_3dmp computes the horizontal diffusion tendency
2534	! on model horizontal coordinate surfaces. This routine computes diffusion
2535	! a perturbation scalar (field-base_3d).
2536	!
2537	!</DESCRIPTION>
2538
2539	pr_inv = 1./prandtl
2540	specified = .false.
2541	if(config_flags%specified .or. config_flags%nested) specified = .true.
2542
2543	ktf=MIN(kte,kde-1)
2544
2545	i_start = its
2546	i_end = MIN(ite,ide-1)
2547	j_start = jts
2548	j_end = MIN(jte,jde-1)
2549
2550	IF ( config_flags%open_xs .or. specified ) i_start = MAX(ids+1,its)
2551	IF ( config_flags%open_xe .or. specified ) i_end = MIN(ide-2,ite)
2552	IF ( config_flags%open_ys .or. specified ) j_start = MAX(jds+1,jts)
2553	IF ( config_flags%open_ye .or. specified ) j_end = MIN(jde-2,jte)
2554	IF ( config_flags%periodic_x ) i_start = its
2555	IF ( config_flags%periodic_x ) i_end = MIN(ite,ide-1)
2556
2557	DO j = j_start, j_end
2558	DO k=kts,ktf
2559	DO i = i_start, i_end
2560
2561	mkrdxm=msfu(i,j)0.5(xkmhd(i,k,j)+xkmhd(i-1,k,j))0.5(mu(i,j)+mu(i-1,j))rdxpr_inv
2562	mkrdxp=msfu(i+1,j)0.5(xkmhd(i+1,k,j)+xkmhd(i,k,j))0.5(mu(i+1,j)+mu(i,j))rdxpr_inv
2563	mrdx=msft(i,j)*rdx
2564	mkrdym=msfv(i,j)0.5(xkmhd(i,k,j)+xkmhd(i,k,j-1))0.5(mu(i,j)+mu(i,j-1))rdypr_inv
2565	mkrdyp=msfv(i,j+1)0.5(xkmhd(i,k,j+1)+xkmhd(i,k,j))0.5(mu(i,j+1)+mu(i,j))rdypr_inv
2566	mrdy=msft(i,j)*rdy
2567
2568	tendency(i,k,j)=tendency(i,k,j)+( &
2569	mrdx( mkrdxp( field(i+1,k,j) -field(i ,k,j) &
2570	-base_3d(i+1,k,j)+base_3d(i ,k,j) ) &
2571	-mkrdxm*( field(i ,k,j) -field(i-1,k,j) &
2572	-base_3d(i ,k,j)+base_3d(i-1,k,j) ) ) &
2573	+mrdy( mkrdyp( field(i,k,j+1) -field(i,k,j ) &
2574	-base_3d(i,k,j+1)+base_3d(i,k,j ) ) &
2575	-mkrdym*( field(i,k,j ) -field(i,k,j-1) &
2576	-base_3d(i,k,j )+base_3d(i,k,j-1) ) ) &
2577	)
2578	ENDDO
2579	ENDDO
2580	ENDDO
2581
2582	END SUBROUTINE horizontal_diffusion_3dmp
2583
2584	!-----------------------------------------------------------------------------------------
2585
2586	SUBROUTINE vertical_diffusion ( name, field, tendency, &
2587	config_flags, &
2588	alt, mut, rdn, rdnw, kvdif, &
2589	ids, ide, jds, jde, kds, kde, &
2590	ims, ime, jms, jme, kms, kme, &
2591	its, ite, jts, jte, kts, kte )
2592
2593
2594	IMPLICIT NONE
2595
2596	! Input data
2597
2598	TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
2599
2600	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
2601	ims, ime, jms, jme, kms, kme, &
2602	its, ite, jts, jte, kts, kte
2603
2604	CHARACTER(LEN=1) , INTENT(IN ) :: name
2605
2606	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , &
2607	INTENT(IN ) :: field, &
2608	alt
2609
2610	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency
2611
2612	REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mut
2613
2614	REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: rdn, rdnw
2615
2616	REAL , INTENT(IN ) :: kvdif
2617
2618	! Local data
2619
2620	INTEGER :: i, j, k, itf, jtf, ktf
2621	INTEGER :: i_start, i_end, j_start, j_end
2622
2623	REAL , DIMENSION(its:ite, jts:jte) :: vfluxm, vfluxp, zz
2624	REAL , DIMENSION(its:ite, 0:kte+1) :: vflux
2625
2626	REAL :: rdz
2627
2628	LOGICAL :: specified
2629
2630	!<DESCRIPTION>
2631	!
2632	! vertical_diffusion
2633	! computes vertical diffusion tendency.
2634	!
2635	!</DESCRIPTION>
2636
2637	specified = .false.
2638	if(config_flags%specified .or. config_flags%nested) specified = .true.
2639
2640	ktf=MIN(kte,kde-1)
2641
2642	IF (name .EQ. 'w')THEN
2643
2644
2645	i_start = its
2646	i_end = MIN(ite,ide-1)
2647	j_start = jts
2648	j_end = MIN(jte,jde-1)
2649
2650	j_loop_w : DO j = j_start, j_end
2651
2652	DO k=kts,ktf-1
2653	DO i = i_start, i_end
2654	vflux(i,k)= (kvdif/alt(i,k,j))rdnw(k)(field(i,k+1,j)-field(i,k,j))
2655	ENDDO
2656	ENDDO
2657
2658	DO i = i_start, i_end
2659	vflux(i,ktf)=0.
2660	ENDDO
2661
2662	DO k=kts+1,ktf
2663	DO i = i_start, i_end
2664	tendency(i,k,j)=tendency(i,k,j) &
2665	+rdn(k)gg/mut(i,j)/(0.5*(alt(i,k,j)+alt(i,k-1,j))) &
2666	*(vflux(i,k)-vflux(i,k-1))
2667	ENDDO
2668	ENDDO
2669
2670	ENDDO j_loop_w
2671
2672	ELSE IF(name .EQ. 'm')THEN
2673
2674	i_start = its
2675	i_end = MIN(ite,ide-1)
2676	j_start = jts
2677	j_end = MIN(jte,jde-1)
2678
2679	j_loop_s : DO j = j_start, j_end
2680
2681	DO k=kts,ktf-1
2682	DO i = i_start, i_end
2683	vflux(i,k)=kvdifrdn(k+1)/(0.5(alt(i,k,j)+alt(i,k+1,j))) &
2684	*(field(i,k+1,j)-field(i,k,j))
2685	ENDDO
2686	ENDDO
2687
2688	DO i = i_start, i_end
2689	vflux(i,0)=vflux(i,1)
2690	ENDDO
2691
2692	DO i = i_start, i_end
2693	vflux(i,ktf)=0.
2694	ENDDO
2695
2696	DO k=kts,ktf
2697	DO i = i_start, i_end
2698	tendency(i,k,j)=tendency(i,k,j)+g*g/mut(i,j)/alt(i,k,j) &
2699	rdnw(k)(vflux(i,k)-vflux(i,k-1))
2700	ENDDO
2701	ENDDO
2702
2703	ENDDO j_loop_s
2704
2705	ENDIF
2706
2707	END SUBROUTINE vertical_diffusion
2708
2709
2710	!-------------------------------------------------------------------------------
2711
2712	SUBROUTINE vertical_diffusion_mp ( field, tendency, config_flags, &
2713	base, &
2714	alt, mut, rdn, rdnw, kvdif, &
2715	ids, ide, jds, jde, kds, kde, &
2716	ims, ime, jms, jme, kms, kme, &
2717	its, ite, jts, jte, kts, kte )
2718
2719
2720	IMPLICIT NONE
2721
2722	! Input data
2723
2724	TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
2725
2726	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
2727	ims, ime, jms, jme, kms, kme, &
2728	its, ite, jts, jte, kts, kte
2729
2730	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , &
2731	INTENT(IN ) :: field, &
2732	alt
2733
2734	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency
2735
2736	REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mut
2737
2738	REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: rdn, &
2739	rdnw, &
2740	base
2741
2742	REAL , INTENT(IN ) :: kvdif
2743
2744	! Local data
2745
2746	INTEGER :: i, j, k, itf, jtf, ktf
2747	INTEGER :: i_start, i_end, j_start, j_end
2748
2749	REAL , DIMENSION(its:ite, 0:kte+1) :: vflux
2750
2751	REAL :: rdz
2752
2753	LOGICAL :: specified
2754
2755	!<DESCRIPTION>
2756	!
2757	! vertical_diffusion_mp
2758	! computes vertical diffusion tendency of a perturbation variable
2759	! (field-base). Note that base as a 1D (k) field.
2760	!
2761	!</DESCRIPTION>
2762
2763	specified = .false.
2764	if(config_flags%specified .or. config_flags%nested) specified = .true.
2765
2766	ktf=MIN(kte,kde-1)
2767
2768	i_start = its
2769	i_end = MIN(ite,ide-1)
2770	j_start = jts
2771	j_end = MIN(jte,jde-1)
2772
2773	j_loop_s : DO j = j_start, j_end
2774
2775	DO k=kts,ktf-1
2776	DO i = i_start, i_end
2777	vflux(i,k)=kvdifrdn(k+1)/(0.5(alt(i,k,j)+alt(i,k+1,j))) &
2778	*(field(i,k+1,j)-field(i,k,j)-base(k+1)+base(k))
2779	ENDDO
2780	ENDDO
2781
2782	DO i = i_start, i_end
2783	vflux(i,0)=vflux(i,1)
2784	ENDDO
2785
2786	DO i = i_start, i_end
2787	vflux(i,ktf)=0.
2788	ENDDO
2789
2790	DO k=kts,ktf
2791	DO i = i_start, i_end
2792	tendency(i,k,j)=tendency(i,k,j)+g*g/mut(i,j)/alt(i,k,j) &
2793	rdnw(k)(vflux(i,k)-vflux(i,k-1))
2794	ENDDO
2795	ENDDO
2796
2797	ENDDO j_loop_s
2798
2799	END SUBROUTINE vertical_diffusion_mp
2800
2801
2802	!-------------------------------------------------------------------------------
2803
2804	SUBROUTINE vertical_diffusion_3dmp ( field, tendency, config_flags, &
2805	base_3d, &
2806	alt, mut, rdn, rdnw, kvdif, &
2807	ids, ide, jds, jde, kds, kde, &
2808	ims, ime, jms, jme, kms, kme, &
2809	its, ite, jts, jte, kts, kte )
2810
2811
2812	IMPLICIT NONE
2813
2814	! Input data
2815
2816	TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
2817
2818	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
2819	ims, ime, jms, jme, kms, kme, &
2820	its, ite, jts, jte, kts, kte
2821
2822	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , &
2823	INTENT(IN ) :: field, &
2824	alt, &
2825	base_3d
2826
2827	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency
2828
2829	REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mut
2830
2831	REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: rdn, &
2832	rdnw
2833
2834	REAL , INTENT(IN ) :: kvdif
2835
2836	! Local data
2837
2838	INTEGER :: i, j, k, itf, jtf, ktf
2839	INTEGER :: i_start, i_end, j_start, j_end
2840
2841	REAL , DIMENSION(its:ite, 0:kte+1) :: vflux
2842
2843	REAL :: rdz
2844
2845	LOGICAL :: specified
2846
2847	!<DESCRIPTION>
2848	!
2849	! vertical_diffusion_3dmp
2850	! computes vertical diffusion tendency of a perturbation variable
2851	! (field-base_3d).
2852	!
2853	!</DESCRIPTION>
2854
2855	specified = .false.
2856	if(config_flags%specified .or. config_flags%nested) specified = .true.
2857
2858	ktf=MIN(kte,kde-1)
2859
2860	i_start = its
2861	i_end = MIN(ite,ide-1)
2862	j_start = jts
2863	j_end = MIN(jte,jde-1)
2864
2865	j_loop_s : DO j = j_start, j_end
2866
2867	DO k=kts,ktf-1
2868	DO i = i_start, i_end
2869	vflux(i,k)=kvdifrdn(k+1)/(0.5(alt(i,k,j)+alt(i,k+1,j))) &
2870	*( field(i,k+1,j) -field(i,k,j) &
2871	-base_3d(i,k+1,j)+base_3d(i,k,j) )
2872	ENDDO
2873	ENDDO
2874
2875	DO i = i_start, i_end
2876	vflux(i,0)=vflux(i,1)
2877	ENDDO
2878
2879	DO i = i_start, i_end
2880	vflux(i,ktf)=0.
2881	ENDDO
2882
2883	DO k=kts,ktf
2884	DO i = i_start, i_end
2885	tendency(i,k,j)=tendency(i,k,j)+g*g/mut(i,j)/alt(i,k,j) &
2886	rdnw(k)(vflux(i,k)-vflux(i,k-1))
2887	ENDDO
2888	ENDDO
2889
2890	ENDDO j_loop_s
2891
2892	END SUBROUTINE vertical_diffusion_3dmp
2893
2894
2895	!-------------------------------------------------------------------------------
2896
2897
2898	SUBROUTINE vertical_diffusion_u ( field, tendency, &
2899	config_flags, u_base, &
2900	alt, muu, rdn, rdnw, kvdif, &
2901	ids, ide, jds, jde, kds, kde, &
2902	ims, ime, jms, jme, kms, kme, &
2903	its, ite, jts, jte, kts, kte )
2904
2905
2906	IMPLICIT NONE
2907
2908	! Input data
2909
2910	TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
2911
2912	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
2913	ims, ime, jms, jme, kms, kme, &
2914	its, ite, jts, jte, kts, kte
2915
2916	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , &
2917	INTENT(IN ) :: field, &
2918	alt
2919
2920	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency
2921
2922	REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: muu
2923
2924	REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: rdn, rdnw, u_base
2925
2926	REAL , INTENT(IN ) :: kvdif
2927
2928	! Local data
2929
2930	INTEGER :: i, j, k, itf, jtf, ktf
2931	INTEGER :: i_start, i_end, j_start, j_end
2932
2933	REAL , DIMENSION(its:ite, 0:kte+1) :: vflux
2934
2935	REAL :: rdz, zz
2936
2937	LOGICAL :: specified
2938
2939	!<DESCRIPTION>
2940	!
2941	! vertical_diffusion_u computes vertical diffusion tendency for
2942	! the u momentum equation. This routine assumes a constant eddy
2943	! viscosity kvdif.
2944	!
2945	!</DESCRIPTION>
2946
2947	specified = .false.
2948	if(config_flags%specified .or. config_flags%nested) specified = .true.
2949
2950	ktf=MIN(kte,kde-1)
2951
2952	i_start = its
2953	i_end = ite
2954	j_start = jts
2955	j_end = MIN(jte,jde-1)
2956
2957	IF ( config_flags%open_xs .or. specified ) i_start = MAX(ids+1,its)
2958	IF ( config_flags%open_xe .or. specified ) i_end = MIN(ide-1,ite)
2959	IF ( config_flags%periodic_x ) i_start = its
2960	IF ( config_flags%periodic_x ) i_end = ite
2961
2962
2963	j_loop_u : DO j = j_start, j_end
2964
2965	DO k=kts,ktf-1
2966	DO i = i_start, i_end
2967	vflux(i,k)=kvdifrdn(k+1)/(0.25( alt(i ,k ,j) &
2968	+alt(i-1,k ,j) &
2969	+alt(i ,k+1,j) &
2970	+alt(i-1,k+1,j) ) ) &
2971	*(field(i,k+1,j)-field(i,k,j) &
2972	-u_base(k+1) +u_base(k) )
2973	ENDDO
2974	ENDDO
2975
2976	DO i = i_start, i_end
2977	vflux(i,0)=vflux(i,1)
2978	ENDDO
2979
2980	DO i = i_start, i_end
2981	vflux(i,ktf)=0.
2982	ENDDO
2983
2984	DO k=kts,ktf-1
2985	DO i = i_start, i_end
2986	tendency(i,k,j)=tendency(i,k,j)+ &
2987	ggrdnw(k)/muu(i,j)/(0.5(alt(i-1,k,j)+alt(i,k,j))) &
2988	(vflux(i,k)-vflux(i,k-1))
2989	ENDDO
2990	ENDDO
2991
2992	ENDDO j_loop_u
2993
2994	END SUBROUTINE vertical_diffusion_u
2995
2996	!-------------------------------------------------------------------------------
2997
2998
2999	SUBROUTINE vertical_diffusion_v ( field, tendency, &
3000	config_flags, v_base, &
3001	alt, muv, rdn, rdnw, kvdif, &
3002	ids, ide, jds, jde, kds, kde, &
3003	ims, ime, jms, jme, kms, kme, &
3004	its, ite, jts, jte, kts, kte )
3005
3006
3007	IMPLICIT NONE
3008
3009	! Input data
3010
3011	TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
3012
3013	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
3014	ims, ime, jms, jme, kms, kme, &
3015	its, ite, jts, jte, kts, kte
3016
3017	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , &
3018	INTENT(IN ) :: field, &
3019	alt
3020	REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: rdn, rdnw, v_base
3021
3022	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency
3023
3024	REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: muv
3025
3026	REAL , INTENT(IN ) :: kvdif
3027
3028	! Local data
3029
3030	INTEGER :: i, j, k, itf, jtf, ktf, jm1
3031	INTEGER :: i_start, i_end, j_start, j_end
3032
3033	REAL , DIMENSION(its:ite, 0:kte+1) :: vflux
3034
3035	REAL :: rdz, zz
3036
3037	LOGICAL :: specified
3038
3039	!<DESCRIPTION>
3040	!
3041	! vertical_diffusion_v computes vertical diffusion tendency for
3042	! the v momentum equation. This routine assumes a constant eddy
3043	! viscosity kvdif.
3044	!
3045	!</DESCRIPTION>
3046
3047	specified = .false.
3048	if(config_flags%specified .or. config_flags%nested) specified = .true.
3049
3050	ktf=MIN(kte,kde-1)
3051
3052	i_start = its
3053	i_end = MIN(ite,ide-1)
3054	j_start = jts
3055	j_end = MIN(jte,jde-1)
3056
3057	IF ( config_flags%open_ys .or. specified ) j_start = MAX(jds+1,jts)
3058	IF ( config_flags%open_ye .or. specified ) j_end = MIN(jde-1,jte)
3059
3060	j_loop_v : DO j = j_start, j_end
3061	! jm1 = max(j-1,1)
3062	jm1 = j-1
3063
3064	DO k=kts,ktf-1
3065	DO i = i_start, i_end
3066	vflux(i,k)=kvdifrdn(k+1)/(0.25( alt(i,k ,j ) &
3067	+alt(i,k ,jm1) &
3068	+alt(i,k+1,j ) &
3069	+alt(i,k+1,jm1) ) ) &
3070	*(field(i,k+1,j)-field(i,k,j) &
3071	-v_base(k+1) +v_base(k) )
3072	ENDDO
3073	ENDDO
3074
3075	DO i = i_start, i_end
3076	vflux(i,0)=vflux(i,1)
3077	ENDDO
3078
3079	DO i = i_start, i_end
3080	vflux(i,ktf)=0.
3081	ENDDO
3082
3083	DO k=kts,ktf-1
3084	DO i = i_start, i_end
3085	tendency(i,k,j)=tendency(i,k,j)+ &
3086	ggrdnw(k)/muv(i,j)/(0.5(alt(i,k,jm1)+alt(i,k,j))) &
3087	(vflux(i,k)-vflux(i,k-1))
3088	ENDDO
3089	ENDDO
3090
3091	ENDDO j_loop_v
3092
3093	END SUBROUTINE vertical_diffusion_v
3094
3095	!*************** end new mass coordinate routines
3096
3097	!-------------------------------------------------------------------------------
3098
3099	SUBROUTINE calculate_full ( rfield, rfieldb, rfieldp, &
3100	ids, ide, jds, jde, kds, kde, &
3101	ims, ime, jms, jme, kms, kme, &
3102	its, ite, jts, jte, kts, kte )
3103
3104	IMPLICIT NONE
3105
3106	! Input data
3107
3108	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
3109	ims, ime, jms, jme, kms, kme, &
3110	its, ite, jts, jte, kts, kte
3111
3112	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: rfieldb, &
3113	rfieldp
3114
3115	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT ) :: rfield
3116
3117	! Local indices.
3118
3119	INTEGER :: i, j, k, itf, jtf, ktf
3120
3121	!<DESCRIPTION>
3122	!
3123	! calculate_full
3124	! calculates full 3D field from pertubation and base field.
3125	!
3126	!</DESCRIPTION>
3127
3128	itf=MIN(ite,ide-1)
3129	jtf=MIN(jte,jde-1)
3130	ktf=MIN(kte,kde-1)
3131
3132	DO j=jts,jtf
3133	DO k=kts,ktf
3134	DO i=its,itf
3135	rfield(i,k,j)=rfieldb(i,k,j)+rfieldp(i,k,j)
3136	ENDDO
3137	ENDDO
3138	ENDDO
3139
3140	END SUBROUTINE calculate_full
3141
3142	!------------------------------------------------------------------------------
3143
3144	SUBROUTINE coriolis ( ru, rv, rw, ru_tend, rv_tend, rw_tend, &
3145	config_flags, &
3146	f, e, sina, cosa, fzm, fzp, &
3147	ids, ide, jds, jde, kds, kde, &
3148	ims, ime, jms, jme, kms, kme, &
3149	its, ite, jts, jte, kts, kte )
3150
3151	IMPLICIT NONE
3152
3153	! Input data
3154
3155	TYPE(grid_config_rec_type) , INTENT(IN ) :: config_flags
3156
3157	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
3158	ims, ime, jms, jme, kms, kme, &
3159	its, ite, jts, jte, kts, kte
3160
3161	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: ru_tend, &
3162	rv_tend, &
3163	rw_tend
3164	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: ru, &
3165	rv, &
3166	rw
3167
3168	REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: f, &
3169	e, &
3170	sina, &
3171	cosa
3172
3173	REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fzm, &
3174	fzp
3175
3176	! Local indices.
3177
3178	INTEGER :: i, j , k, ktf
3179	INTEGER :: i_start, i_end, j_start, j_end
3180
3181	LOGICAL :: specified
3182
3183	!<DESCRIPTION>
3184	!
3185	! coriolis calculates the large timestep tendency terms in the
3186	! u, v, and w momentum equations arise from the coriolis force.
3187	!
3188	!</DESCRIPTION>
3189
3190	specified = .false.
3191	if(config_flags%specified .or. config_flags%nested) specified = .true.
3192
3193	ktf=MIN(kte,kde-1)
3194
3195	! coriolis for u-momentum equation
3196
3197	i_start = its
3198	i_end = ite
3199	IF ( config_flags%open_xs .or. specified .or. &
3200	config_flags%nested) i_start = MAX(ids+1,its)
3201	IF ( config_flags%open_xe .or. specified .or. &
3202	config_flags%nested) i_end = MIN(ide-1,ite)
3203	IF ( config_flags%periodic_x ) i_start = its
3204	IF ( config_flags%periodic_x ) i_end = ite
3205
3206	DO j = jts, MIN(jte,jde-1)
3207
3208	DO k=kts,ktf
3209	DO i = i_start, i_end
3210
3211	ru_tend(i,k,j)=ru_tend(i,k,j) + 0.5*(f(i,j)+f(i-1,j)) &
3212	0.25(rv(i-1,k,j+1)+rv(i,k,j+1)+rv(i-1,k,j)+rv(i,k,j)) &
3213	- 0.5(e(i,j)+e(i-1,j))0.5*(cosa(i,j)+cosa(i-1,j)) &
3214	0.25(rw(i-1,k+1,j)+rw(i-1,k,j)+rw(i,k+1,j)+rw(i,k,j))
3215
3216	ENDDO
3217	ENDDO
3218
3219	IF ( (config_flags%open_xs) .and. (its == ids) ) THEN
3220
3221	DO k=kts,ktf
3222
3223	ru_tend(its,k,j)=ru_tend(its,k,j) + 0.5*(f(its,j)+f(its,j)) &
3224	0.25(rv(its,k,j+1)+rv(its,k,j+1)+rv(its,k,j)+rv(its,k,j)) &
3225	- 0.5(e(its,j)+e(its,j))0.5*(cosa(its,j)+cosa(its,j)) &
3226	0.25(rw(its,k+1,j)+rw(its,k,j)+rw(its,k+1,j)+rw(its,k,j))
3227
3228	ENDDO
3229
3230	ENDIF
3231
3232	IF ( (config_flags%open_xe) .and. (ite == ide) ) THEN
3233
3234	DO k=kts,ktf
3235
3236	ru_tend(ite,k,j)=ru_tend(ite,k,j) + 0.5*(f(ite-1,j)+f(ite-1,j)) &
3237	0.25(rv(ite-1,k,j+1)+rv(ite-1,k,j+1)+rv(ite-1,k,j)+rv(ite-1,k,j)) &
3238	- 0.5(e(ite-1,j)+e(ite-1,j))0.5*(cosa(ite-1,j)+cosa(ite-1,j)) &
3239	0.25(rw(ite-1,k+1,j)+rw(ite-1,k,j)+rw(ite-1,k+1,j)+rw(ite-1,k,j))
3240
3241	ENDDO
3242
3243	ENDIF
3244
3245	ENDDO
3246
3247	! coriolis term for v-momentum equation
3248
3249	j_start = jts
3250	j_end = jte
3251
3252	IF ( config_flags%open_ys .or. specified .or. &
3253	config_flags%nested) j_start = MAX(jds+1,jts)
3254	IF ( config_flags%open_ye .or. specified .or. &
3255	config_flags%nested) j_end = MIN(jde-1,jte)
3256
3257	IF ( (config_flags%open_ys) .and. (jts == jds) ) THEN
3258
3259	DO k=kts,ktf
3260	DO i=its,MIN(ide-1,ite)
3261
3262	rv_tend(i,k,jts)=rv_tend(i,k,jts) - 0.5*(f(i,jts)+f(i,jts)) &
3263	0.25(ru(i,k,jts)+ru(i+1,k,jts)+ru(i,k,jts)+ru(i+1,k,jts)) &
3264	+ 0.5(e(i,jts)+e(i,jts))0.5*(sina(i,jts)+sina(i,jts)) &
3265	0.25(rw(i,k+1,jts)+rw(i,k,jts)+rw(i,k+1,jts)+rw(i,k,jts))
3266
3267	ENDDO
3268	ENDDO
3269
3270	ENDIF
3271
3272	DO j=j_start, j_end
3273	DO k=kts,ktf
3274	DO i=its,MIN(ide-1,ite)
3275
3276	rv_tend(i,k,j)=rv_tend(i,k,j) - 0.5*(f(i,j)+f(i,j-1)) &
3277	0.25(ru(i,k,j)+ru(i+1,k,j)+ru(i,k,j-1)+ru(i+1,k,j-1)) &
3278	+ 0.5(e(i,j)+e(i,j-1))0.5*(sina(i,j)+sina(i,j-1)) &
3279	0.25(rw(i,k+1,j-1)+rw(i,k,j-1)+rw(i,k+1,j)+rw(i,k,j))
3280
3281	ENDDO
3282	ENDDO
3283	ENDDO
3284
3285
3286	IF ( (config_flags%open_ye) .and. (jte == jde) ) THEN
3287
3288	DO k=kts,ktf
3289	DO i=its,MIN(ide-1,ite)
3290
3291	rv_tend(i,k,jte)=rv_tend(i,k,jte) - 0.5*(f(i,jte-1)+f(i,jte-1)) &
3292	0.25(ru(i,k,jte-1)+ru(i+1,k,jte-1)+ru(i,k,jte-1)+ru(i+1,k,jte-1)) &
3293	+ 0.5(e(i,jte-1)+e(i,jte-1))0.5*(sina(i,jte-1)+sina(i,jte-1)) &
3294	0.25(rw(i,k+1,jte-1)+rw(i,k,jte-1)+rw(i,k+1,jte-1)+rw(i,k,jte-1))
3295
3296	ENDDO
3297	ENDDO
3298
3299	ENDIF
3300
3301	! coriolis term for w-mometum
3302
3303	DO j=jts,MIN(jte, jde-1)
3304	DO k=kts+1,ktf
3305	DO i=its,MIN(ite, ide-1)
3306
3307	rw_tend(i,k,j)=rw_tend(i,k,j) + e(i,j)* &
3308	(cosa(i,j)0.5(fzm(k)*(ru(i,k,j)+ru(i+1,k,j)) &
3309	+fzp(k)*(ru(i,k-1,j)+ru(i+1,k-1,j))) &
3310	-sina(i,j)0.5(fzm(k)*(rv(i,k,j)+rv(i,k,j+1)) &
3311	+fzp(k)*(rv(i,k-1,j)+rv(i,k-1,j+1))))
3312
3313	ENDDO
3314	ENDDO
3315	ENDDO
3316
3317	END SUBROUTINE coriolis
3318
3319	!------------------------------------------------------------------------------
3320
3321	SUBROUTINE perturbation_coriolis ( ru_in, rv_in, rw, ru_tend, rv_tend, rw_tend, &
3322	config_flags, &
3323	u_base, v_base, z_base, &
3324	muu, muv, phb, ph, &
3325	f, e, sina, cosa, fzm, fzp, &
3326	ids, ide, jds, jde, kds, kde, &
3327	ims, ime, jms, jme, kms, kme, &
3328	its, ite, jts, jte, kts, kte )
3329
3330	IMPLICIT NONE
3331
3332	! Input data
3333
3334	TYPE(grid_config_rec_type) , INTENT(IN ) :: config_flags
3335
3336	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
3337	ims, ime, jms, jme, kms, kme, &
3338	its, ite, jts, jte, kts, kte
3339
3340	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: ru_tend, &
3341	rv_tend, &
3342	rw_tend
3343	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: ru_in, &
3344	rv_in, &
3345	rw, &
3346	ph, &
3347	phb
3348
3349
3350	REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: f, &
3351	e, &
3352	sina, &
3353	cosa
3354
3355	REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: muu, &
3356	muv
3357
3358
3359	REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fzm, &
3360	fzp
3361
3362	REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: u_base, &
3363	v_base, &
3364	z_base
3365
3366	! Local storage
3367
3368	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) :: ru, &
3369	rv
3370
3371	REAL :: z_at_u, z_at_v, wkp1, wk, wkm1
3372
3373	! Local indices.
3374
3375	INTEGER :: i, j , k, ktf
3376	INTEGER :: i_start, i_end, j_start, j_end
3377
3378	LOGICAL :: specified
3379
3380	!<DESCRIPTION>
3381	!
3382	! perturbation_coriolis calculates the large timestep tendency terms in the
3383	! u, v, and w momentum equations arise from the coriolis force. This version
3384	! subtracts off the horizontal velocities from the initial sounding when
3385	! computing the forcing terms, hence "perturbation" coriolis.
3386	!
3387	!</DESCRIPTION>
3388
3389	specified = .false.
3390	if(config_flags%specified .or. config_flags%nested) specified = .true.
3391
3392	ktf=MIN(kte,kde-1)
3393
3394	! coriolis for u-momentum equation
3395
3396	i_start = its
3397	i_end = ite
3398	IF ( config_flags%open_xs .or. specified .or. &
3399	config_flags%nested) i_start = MAX(ids+1,its)
3400	IF ( config_flags%open_xe .or. specified .or. &
3401	config_flags%nested) i_end = MIN(ide-1,ite)
3402	IF ( config_flags%periodic_x ) i_start = its
3403	IF ( config_flags%periodic_x ) i_end = ite
3404
3405	! compute perturbation mu*v for use in u momentum equation
3406
3407	DO j = jts, MIN(jte,jde-1)+1
3408	DO k=kts+1,ktf-1
3409	DO i = i_start-1, i_end
3410	z_at_v = 0.25*( phb(i,k,j )+phb(i,k+1,j ) &
3411	+phb(i,k,j-1)+phb(i,k+1,j-1) &
3412	+ph(i,k,j )+ph(i,k+1,j ) &
3413	+ph(i,k,j-1)+ph(i,k+1,j-1))/g
3414	wkp1 = min(1.,max(0.,z_at_v-z_base(k))/(z_base(k+1)-z_base(k)))
3415	wkm1 = min(1.,max(0.,z_base(k)-z_at_v)/(z_base(k)-z_base(k-1)))
3416	wk = 1.-wkp1-wkm1
3417	rv(i,k,j) = rv_in(i,k,j) - muv(i,j)*( &
3418	wkm1*v_base(k-1) &
3419	+wk *v_base(k ) &
3420	+wkp1*v_base(k+1) )
3421	ENDDO
3422	ENDDO
3423	ENDDO
3424
3425
3426	! pick up top and bottom v
3427
3428	DO j = jts, MIN(jte,jde-1)+1
3429	DO i = i_start-1, i_end
3430
3431	k = kts
3432	z_at_v = 0.25*( phb(i,k,j )+phb(i,k+1,j ) &
3433	+phb(i,k,j-1)+phb(i,k+1,j-1) &
3434	+ph(i,k,j )+ph(i,k+1,j ) &
3435	+ph(i,k,j-1)+ph(i,k+1,j-1))/g
3436	wkp1 = min(1.,max(0.,z_at_v-z_base(k))/(z_base(k+1)-z_base(k)))
3437	wk = 1.-wkp1
3438	rv(i,k,j) = rv_in(i,k,j) - muv(i,j)*( &
3439	+wk *v_base(k ) &
3440	+wkp1*v_base(k+1) )
3441
3442	k = ktf
3443	z_at_v = 0.25*( phb(i,k,j )+phb(i,k+1,j ) &
3444	+phb(i,k,j-1)+phb(i,k+1,j-1) &
3445	+ph(i,k,j )+ph(i,k+1,j ) &
3446	+ph(i,k,j-1)+ph(i,k+1,j-1))/g
3447	wkm1 = min(1.,max(0.,z_base(k)-z_at_v)/(z_base(k)-z_base(k-1)))
3448	wk = 1.-wkm1
3449	rv(i,k,j) = rv_in(i,k,j) - muv(i,j)*( &
3450	wkm1*v_base(k-1) &
3451	+wk *v_base(k ) )
3452
3453	ENDDO
3454	ENDDO
3455
3456	! compute coriolis forcing for u
3457
3458	DO j = jts, MIN(jte,jde-1)
3459
3460	DO k=kts,ktf
3461	DO i = i_start, i_end
3462	ru_tend(i,k,j)=ru_tend(i,k,j) + 0.5*(f(i,j)+f(i-1,j)) &
3463	0.25(rv(i-1,k,j+1)+rv(i,k,j+1)+rv(i-1,k,j)+rv(i,k,j)) &
3464	- 0.5(e(i,j)+e(i-1,j))0.5*(cosa(i,j)+cosa(i-1,j)) &
3465	0.25(rw(i-1,k+1,j)+rw(i-1,k,j)+rw(i,k+1,j)+rw(i,k,j))
3466	ENDDO
3467	ENDDO
3468
3469	IF ( (config_flags%open_xs) .and. (its == ids) ) THEN
3470
3471	DO k=kts,ktf
3472
3473	ru_tend(its,k,j)=ru_tend(its,k,j) + 0.5*(f(its,j)+f(its,j)) &
3474	0.25(rv(its,k,j+1)+rv(its,k,j+1)+rv(its,k,j)+rv(its,k,j)) &
3475	- 0.5(e(its,j)+e(its,j))0.5*(cosa(its,j)+cosa(its,j)) &
3476	0.25(rw(its,k+1,j)+rw(its,k,j)+rw(its,k+1,j)+rw(its,k,j))
3477
3478	ENDDO
3479
3480	ENDIF
3481
3482	IF ( (config_flags%open_xe) .and. (ite == ide) ) THEN
3483
3484	DO k=kts,ktf
3485
3486	ru_tend(ite,k,j)=ru_tend(ite,k,j) + 0.5*(f(ite-1,j)+f(ite-1,j)) &
3487	0.25(rv(ite-1,k,j+1)+rv(ite-1,k,j+1)+rv(ite-1,k,j)+rv(ite-1,k,j)) &
3488	- 0.5(e(ite-1,j)+e(ite-1,j))0.5*(cosa(ite-1,j)+cosa(ite-1,j)) &
3489	0.25(rw(ite-1,k+1,j)+rw(ite-1,k,j)+rw(ite-1,k+1,j)+rw(ite-1,k,j))
3490
3491	ENDDO
3492
3493	ENDIF
3494
3495	ENDDO
3496
3497	! coriolis term for v-momentum equation
3498
3499	j_start = jts
3500	j_end = jte
3501
3502	IF ( config_flags%open_ys .or. specified .or. &
3503	config_flags%nested) j_start = MAX(jds+1,jts)
3504	IF ( config_flags%open_ye .or. specified .or. &
3505	config_flags%nested) j_end = MIN(jde-1,jte)
3506
3507	! compute perturbation mu*u for use in v momentum equation
3508
3509	DO j = j_start-1,j_end
3510	DO k=kts+1,ktf-1
3511	DO i = its, MIN(ite,ide-1)+1
3512	z_at_u = 0.25*( phb(i ,k,j)+phb(i ,k+1,j) &
3513	+phb(i-1,k,j)+phb(i-1,k+1,j) &
3514	+ph(i ,k,j)+ph(i ,k+1,j) &
3515	+ph(i-1,k,j)+ph(i-1,k+1,j))/g
3516	wkp1 = min(1.,max(0.,z_at_u-z_base(k))/(z_base(k+1)-z_base(k)))
3517	wkm1 = min(1.,max(0.,z_base(k)-z_at_u)/(z_base(k)-z_base(k-1)))
3518	wk = 1.-wkp1-wkm1
3519	ru(i,k,j) = ru_in(i,k,j) - muu(i,j)*( &
3520	wkm1*u_base(k-1) &
3521	+wk *u_base(k ) &
3522	+wkp1*u_base(k+1) )
3523	ENDDO
3524	ENDDO
3525	ENDDO
3526
3527	! pick up top and bottom u
3528
3529	DO j = j_start-1,j_end
3530	DO i = its, MIN(ite,ide-1)+1
3531
3532	k = kts
3533	z_at_u = 0.25*( phb(i ,k,j)+phb(i ,k+1,j) &
3534	+phb(i-1,k,j)+phb(i-1,k+1,j) &
3535	+ph(i ,k,j)+ph(i ,k+1,j) &
3536	+ph(i-1,k,j)+ph(i-1,k+1,j))/g
3537	wkp1 = min(1.,max(0.,z_at_u-z_base(k))/(z_base(k+1)-z_base(k)))
3538	wk = 1.-wkp1
3539	ru(i,k,j) = ru_in(i,k,j) - muu(i,j)*( &
3540	+wk *u_base(k ) &
3541	+wkp1*u_base(k+1) )
3542
3543
3544	k = ktf
3545	z_at_u = 0.25*( phb(i ,k,j)+phb(i ,k+1,j) &
3546	+phb(i-1,k,j)+phb(i-1,k+1,j) &
3547	+ph(i ,k,j)+ph(i ,k+1,j) &
3548	+ph(i-1,k,j)+ph(i-1,k+1,j))/g
3549	wkm1 = min(1.,max(0.,z_base(k)-z_at_u)/(z_base(k)-z_base(k-1)))
3550	wk = 1.-wkm1
3551	ru(i,k,j) = ru_in(i,k,j) - muu(i,j)*( &
3552	wkm1*u_base(k-1) &
3553	+wk *u_base(k ) )
3554
3555	ENDDO
3556	ENDDO
3557
3558	! compute coriolis forcing for v momentum equation
3559
3560	IF ( (config_flags%open_ys) .and. (jts == jds) ) THEN
3561
3562	DO k=kts,ktf
3563	DO i=its,MIN(ide-1,ite)
3564
3565	rv_tend(i,k,jts)=rv_tend(i,k,jts) - 0.5*(f(i,jts)+f(i,jts)) &
3566	0.25(ru(i,k,jts)+ru(i+1,k,jts)+ru(i,k,jts)+ru(i+1,k,jts)) &
3567	+ 0.5(e(i,jts)+e(i,jts))0.5*(sina(i,jts)+sina(i,jts)) &
3568	0.25(rw(i,k+1,jts)+rw(i,k,jts)+rw(i,k+1,jts)+rw(i,k,jts))
3569
3570	ENDDO
3571	ENDDO
3572
3573	ENDIF
3574
3575	DO j=j_start, j_end
3576	DO k=kts,ktf
3577	DO i=its,MIN(ide-1,ite)
3578
3579	rv_tend(i,k,j)=rv_tend(i,k,j) - 0.5*(f(i,j)+f(i,j-1)) &
3580	0.25(ru(i,k,j)+ru(i+1,k,j)+ru(i,k,j-1)+ru(i+1,k,j-1)) &
3581	+ 0.5(e(i,j)+e(i,j-1))0.5*(sina(i,j)+sina(i,j-1)) &
3582	0.25(rw(i,k+1,j-1)+rw(i,k,j-1)+rw(i,k+1,j)+rw(i,k,j))
3583
3584	ENDDO
3585	ENDDO
3586	ENDDO
3587
3588
3589	IF ( (config_flags%open_ye) .and. (jte == jde) ) THEN
3590
3591	DO k=kts,ktf
3592	DO i=its,MIN(ide-1,ite)
3593
3594	rv_tend(i,k,jte)=rv_tend(i,k,jte) - 0.5*(f(i,jte-1)+f(i,jte-1)) &
3595	0.25(ru(i,k,jte-1)+ru(i+1,k,jte-1)+ru(i,k,jte-1)+ru(i+1,k,jte-1)) &
3596	+ 0.5(e(i,jte-1)+e(i,jte-1))0.5*(sina(i,jte-1)+sina(i,jte-1)) &
3597	0.25(rw(i,k+1,jte-1)+rw(i,k,jte-1)+rw(i,k+1,jte-1)+rw(i,k,jte-1))
3598
3599	ENDDO
3600	ENDDO
3601
3602	ENDIF
3603
3604	! coriolis term for w-mometum
3605
3606	DO j=jts,MIN(jte, jde-1)
3607	DO k=kts+1,ktf
3608	DO i=its,MIN(ite, ide-1)
3609
3610	rw_tend(i,k,j)=rw_tend(i,k,j) + e(i,j)* &
3611	(cosa(i,j)0.5(fzm(k)*(ru(i,k,j)+ru(i+1,k,j)) &
3612	+fzp(k)*(ru(i,k-1,j)+ru(i+1,k-1,j))) &
3613	-sina(i,j)0.5(fzm(k)*(rv(i,k,j)+rv(i,k,j+1)) &
3614	+fzp(k)*(rv(i,k-1,j)+rv(i,k-1,j+1))))
3615
3616	ENDDO
3617	ENDDO
3618	ENDDO
3619
3620	END SUBROUTINE perturbation_coriolis
3621
3622	!------------------------------------------------------------------------------
3623
3624	SUBROUTINE curvature ( ru, rv, rw, u, v, w, ru_tend, rv_tend, rw_tend, &
3625	config_flags, &
3626	msfu, msfv, fzm, fzp, rdx, rdy, &
3627	ids, ide, jds, jde, kds, kde, &
3628	ims, ime, jms, jme, kms, kme, &
3629	its, ite, jts, jte, kts, kte )
3630
3631
3632	IMPLICIT NONE
3633
3634	! Input data
3635
3636	TYPE(grid_config_rec_type) , INTENT(IN ) :: config_flags
3637
3638	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
3639	ims, ime, jms, jme, kms, kme, &
3640	its, ite, jts, jte, kts, kte
3641
3642	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , &
3643	INTENT(INOUT) :: ru_tend, &
3644	rv_tend, &
3645	rw_tend
3646
3647	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , &
3648	INTENT(IN ) :: ru, &
3649	rv, &
3650	rw, &
3651	u, &
3652	v, &
3653	w
3654
3655	REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: msfu, &
3656	msfv
3657
3658	REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fzm, &
3659	fzp
3660
3661	REAL , INTENT(IN ) :: rdx, &
3662	rdy
3663
3664	! Local data
3665
3666	! INTEGER :: i, j, k, itf, jtf, ktf, kp1, im, ip, jm, jp
3667	INTEGER :: i, j, k, itf, jtf, ktf
3668	INTEGER :: i_start, i_end, j_start, j_end
3669	! INTEGER :: irmin, irmax, jrmin, jrmax
3670
3671	REAL , DIMENSION( its-1:ite , kts:kte, jts-1:jte ) :: vxgm
3672
3673	LOGICAL :: specified
3674
3675	!<DESCRIPTION>
3676	!
3677	! curvature calculates the large timestep tendency terms in the
3678	! u, v, and w momentum equations arise from the curvature terms.
3679	!
3680	!</DESCRIPTION>
3681
3682	specified = .false.
3683	if(config_flags%specified .or. config_flags%nested) specified = .true.
3684
3685	itf=MIN(ite,ide-1)
3686	jtf=MIN(jte,jde-1)
3687	ktf=MIN(kte,kde-1)
3688
3689	! irmin = ims
3690	! irmax = ime
3691	! jrmin = jms
3692	! jrmax = jme
3693	! IF ( config_flags%open_xs ) irmin = ids
3694	! IF ( config_flags%open_xe ) irmax = ide-1
3695	! IF ( config_flags%open_ys ) jrmin = jds
3696	! IF ( config_flags%open_ye ) jrmax = jde-1
3697
3698	! Define v cross grad m at scalar points - vxgm(i,j)
3699
3700	i_start = its-1
3701	i_end = ite
3702	j_start = jts-1
3703	j_end = jte
3704
3705	IF ( ( config_flags%open_xs .or. specified .or. &
3706	config_flags%nested) .and. (its == ids) ) i_start = its
3707	IF ( ( config_flags%open_xe .or. specified .or. &
3708	config_flags%nested) .and. (ite == ide) ) i_end = ite-1
3709	IF ( ( config_flags%open_ys .or. specified .or. &
3710	config_flags%nested) .and. (jts == jds) ) j_start = jts
3711	IF ( ( config_flags%open_ye .or. specified .or. &
3712	config_flags%nested) .and. (jte == jde) ) j_end = jte-1
3713	IF ( config_flags%periodic_x ) i_start = its-1
3714	IF ( config_flags%periodic_x ) i_end = ite
3715
3716	DO j=j_start, j_end
3717	DO k=kts,ktf
3718	DO i=i_start, i_end
3719	vxgm(i,k,j)=0.5(u(i,k,j)+u(i+1,k,j))(msfv(i,j+1)-msfv(i,j))*rdy - &
3720	0.5(v(i,k,j)+v(i,k,j+1))(msfu(i+1,j)-msfu(i,j))*rdx
3721	ENDDO
3722	ENDDO
3723	ENDDO
3724
3725	! Pick up the boundary rows for open (radiation) lateral b.c.
3726	! Rather crude at present, we are assuming there is no
3727	! variation in this term at the boundary.
3728
3729	IF ( ( config_flags%open_xs .or. (specified .AND. .NOT. config_flags%periodic_x) .or. &
3730	config_flags%nested) .and. (its == ids) ) THEN
3731
3732	DO j = jts-1, jte
3733	DO k = kts, ktf
3734	vxgm(its-1,k,j) = vxgm(its,k,j)
3735	ENDDO
3736	ENDDO
3737
3738	ENDIF
3739
3740	IF ( ( config_flags%open_xe .or. (specified .AND. .NOT. config_flags%periodic_x) .or. &
3741	config_flags%nested) .and. (ite == ide) ) THEN
3742
3743	DO j = jts-1, jte
3744	DO k = kts, ktf
3745	vxgm(ite,k,j) = vxgm(ite-1,k,j)
3746	ENDDO
3747	ENDDO
3748
3749	ENDIF
3750
3751	IF ( ( config_flags%open_ys .or. specified .or. &
3752	config_flags%nested) .and. (jts == jds) ) THEN
3753
3754	DO k = kts, ktf
3755	DO i = its-1, ite
3756	vxgm(i,k,jts-1) = vxgm(i,k,jts)
3757	ENDDO
3758	ENDDO
3759
3760	ENDIF
3761
3762	IF ( ( config_flags%open_ye .or. specified .or. &
3763	config_flags%nested) .and. (jte == jde) ) THEN
3764
3765	DO k = kts, ktf
3766	DO i = its-1, ite
3767	vxgm(i,k,jte) = vxgm(i,k,jte-1)
3768	ENDDO
3769	ENDDO
3770
3771	ENDIF
3772
3773	! curvature term for u momentum eqn.
3774
3775	i_start = its
3776	IF ( config_flags%open_xs .or. specified .or. &
3777	config_flags%nested) i_start = MAX ( ids+1 , its )
3778	IF ( config_flags%open_xe .or. specified .or. &
3779	config_flags%nested) i_end = MIN ( ide-1 , ite )
3780	IF ( config_flags%periodic_x ) i_start = its
3781	IF ( config_flags%periodic_x ) i_end = ite
3782
3783	DO j=jts,MIN(jde-1,jte)
3784	DO k=kts,ktf
3785	DO i=i_start,i_end
3786
3787	ru_tend(i,k,j)=ru_tend(i,k,j) + 0.5*(vxgm(i,k,j)+vxgm(i-1,k,j)) &
3788	0.25(rv(i-1,k,j+1)+rv(i,k,j+1)+rv(i-1,k,j)+rv(i,k,j)) &
3789	- u(i,k,j)*reradius &
3790	0.25(rw(i-1,k+1,j)+rw(i-1,k,j)+rw(i,k+1,j)+rw(i,k,j))
3791
3792	ENDDO
3793	ENDDO
3794	ENDDO
3795
3796	! curvature term for v momentum eqn.
3797
3798	j_start = jts
3799	IF ( config_flags%open_ys .or. specified .or. &
3800	config_flags%nested) j_start = MAX ( jds+1 , jts )
3801	IF ( config_flags%open_ye .or. specified .or. &
3802	config_flags%nested) j_end = MIN ( jde-1 , jte )
3803
3804	DO j=j_start,j_end
3805	DO k=kts,ktf
3806	DO i=its,MIN(ite,ide-1)
3807
3808	rv_tend(i,k,j)=rv_tend(i,k,j) - 0.5*(vxgm(i,k,j)+vxgm(i,k,j-1)) &
3809	0.25(ru(i,k,j)+ru(i+1,k,j)+ru(i,k,j-1)+ru(i+1,k,j-1)) &
3810	- v(i,k,j)*reradius &
3811	0.25(rw(i,k+1,j-1)+rw(i,k,j-1)+rw(i,k+1,j)+rw(i,k,j))
3812
3813	!
3814	! correction version 2.2.1
3815	!
3816
3817
3818	ENDDO
3819	ENDDO
3820	ENDDO
3821
3822	! curvature term for vertical momentum eqn.
3823
3824	DO j=jts,MIN(jte,jde-1)
3825	DO k=MAX(2,kts),ktf
3826	DO i=its,MIN(ite,ide-1)
3827
3828	rw_tend(i,k,j)=rw_tend(i,k,j) + reradius* &
3829	(0.5(fzm(k)(ru(i,k,j)+ru(i+1,k,j))+fzp(k)*(ru(i,k-1,j)+ru(i+1,k-1,j))) &
3830	0.5(fzm(k)( u(i,k,j) +u(i+1,k,j))+fzp(k)( u(i,k-1,j) +u(i+1,k-1,j))) &
3831	+0.5(fzm(k)(rv(i,k,j)+rv(i,k,j+1))+fzp(k)*(rv(i,k-1,j)+rv(i,k-1,j+1))) &
3832	0.5(fzm(k)( v(i,k,j) +v(i,k,j+1))+fzp(k)( v(i,k-1,j) +v(i,k-1,j+1))))
3833
3834	ENDDO
3835	ENDDO
3836	ENDDO
3837
3838	END SUBROUTINE curvature
3839
3840	!------------------------------------------------------------------------------
3841
3842	SUBROUTINE decouple ( rr, rfield, field, name, config_flags, &
3843	fzm, fzp, &
3844	ids, ide, jds, jde, kds, kde, &
3845	ims, ime, jms, jme, kms, kme, &
3846	its, ite, jts, jte, kts, kte )
3847
3848	IMPLICIT NONE
3849
3850	! Input data
3851
3852	TYPE(grid_config_rec_type) , INTENT(IN ) :: config_flags
3853
3854	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
3855	ims, ime, jms, jme, kms, kme, &
3856	its, ite, jts, jte, kts, kte
3857
3858	CHARACTER(LEN=1) , INTENT(IN ) :: name
3859
3860	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: rfield
3861
3862	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: rr
3863
3864	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT ) :: field
3865
3866	REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fzm, fzp
3867
3868	! Local data
3869
3870	INTEGER :: i, j, k, itf, jtf, ktf
3871
3872	!<DESCRIPTION>
3873	!
3874	! decouple decouples a variable from the column dry-air mass.
3875	!
3876	!</DESCRIPTION>
3877
3878	ktf=MIN(kte,kde-1)
3879
3880	IF (name .EQ. 'u')THEN
3881	itf=ite
3882	jtf=MIN(jte,jde-1)
3883
3884	DO j=jts,jtf
3885	DO k=kts,ktf
3886	DO i=its,itf
3887	field(i,k,j)=rfield(i,k,j)/(0.5*(rr(i,k,j)+rr(i-1,k,j)))
3888	ENDDO
3889	ENDDO
3890	ENDDO
3891
3892	ELSE IF (name .EQ. 'v')THEN
3893	itf=MIN(ite,ide-1)
3894	jtf=jte
3895
3896	DO j=jts,jtf
3897	DO k=kts,ktf
3898	DO i=its,itf
3899	field(i,k,j)=rfield(i,k,j)/(0.5*(rr(i,k,j)+rr(i,k,j-1)))
3900	ENDDO
3901	ENDDO
3902	ENDDO
3903
3904	ELSE IF (name .EQ. 'w')THEN
3905	itf=MIN(ite,ide-1)
3906	jtf=MIN(jte,jde-1)
3907	DO j=jts,jtf
3908	DO k=kts+1,ktf
3909	DO i=its,itf
3910	field(i,k,j)=rfield(i,k,j)/(fzm(k)rr(i,k,j)+fzp(k)rr(i,k-1,j))
3911	ENDDO
3912	ENDDO
3913	ENDDO
3914
3915	DO j=jts,jtf
3916	DO i=its,itf
3917	field(i,kte,j) = 0.
3918	ENDDO
3919	ENDDO
3920
3921	ELSE
3922	itf=MIN(ite,ide-1)
3923	jtf=MIN(jte,jde-1)
3924	! For theta we will decouple tb and tp and add them to give t afterwards
3925	DO j=jts,jtf
3926	DO k=kts,ktf
3927	DO i=its,itf
3928	field(i,k,j)=rfield(i,k,j)/rr(i,k,j)
3929	ENDDO
3930	ENDDO
3931	ENDDO
3932
3933	ENDIF
3934
3935	END SUBROUTINE decouple
3936
3937	!-------------------------------------------------------------------------------
3938
3939
3940	SUBROUTINE zero_tend ( tendency, &
3941	ids, ide, jds, jde, kds, kde, &
3942	ims, ime, jms, jme, kms, kme, &
3943	its, ite, jts, jte, kts, kte )
3944
3945
3946	IMPLICIT NONE
3947
3948	! Input data
3949
3950	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
3951	ims, ime, jms, jme, kms, kme, &
3952	its, ite, jts, jte, kts, kte
3953
3954	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency
3955
3956	! Local data
3957
3958	INTEGER :: i, j, k, itf, jtf, ktf
3959
3960	!<DESCRIPTION>
3961	!
3962	! zero_tend sets the input tendency array to zero.
3963	!
3964	!</DESCRIPTION>
3965
3966	DO j = jts, jte
3967	DO k = kts, kte
3968	DO i = its, ite
3969	tendency(i,k,j) = 0.
3970	ENDDO
3971	ENDDO
3972	ENDDO
3973
3974	END SUBROUTINE zero_tend
3975
3976	!======================================================================
3977	! physics prep routines
3978	!======================================================================
3979
3980	SUBROUTINE phy_prep ( config_flags, & ! input
3981	mu, muu, muv, u, v, w, p, pb, alt, ph, & ! input
3982	phb, t, tsk, moist, n_moist, & ! input
3983	mu_3d, rho, th_phy, p_phy , pi_phy , & ! output
3984	u_phy, v_phy, w_phy, p8w, t_phy, t8w, & ! output
3985	z, z_at_w, dz8w, & ! output
3986	fzm, fzp, & ! params
3987	RTHRATEN, &
3988	RTHBLTEN, RUBLTEN, RVBLTEN, &
3989	RQVBLTEN, RQCBLTEN, RQIBLTEN, &
3990	RTHCUTEN, RQVCUTEN, RQCCUTEN, &
3991	RQRCUTEN, RQICUTEN, RQSCUTEN, &
3992	RTHFTEN, RQVFTEN, &
3993	RUNDGDTEN, RVNDGDTEN, RTHNDGDTEN, &
3994	RQVNDGDTEN, RMUNDGDTEN, &
3995	ids, ide, jds, jde, kds, kde, &
3996	ims, ime, jms, jme, kms, kme, &
3997	its, ite, jts, jte, kts, kte )
3998	!----------------------------------------------------------------------
3999	IMPLICIT NONE
4000	!----------------------------------------------------------------------
4001
4002	TYPE(grid_config_rec_type) , INTENT(IN ) :: config_flags
4003
4004	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
4005	ims, ime, jms, jme, kms, kme, &
4006	its, ite, jts, jte, kts, kte
4007	INTEGER , INTENT(IN ) :: n_moist
4008
4009	REAL, DIMENSION( ims:ime, kms:kme , jms:jme , n_moist ), INTENT(IN) :: moist
4010
4011
4012	REAL , DIMENSION( ims:ime, jms:jme ), INTENT(IN ) :: TSK, mu, muu, muv
4013
4014	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , &
4015	INTENT( OUT) :: u_phy, &
4016	v_phy, &
4017	w_phy, &
4018	pi_phy, &
4019	p_phy, &
4020	p8w, &
4021	t_phy, &
4022	th_phy, &
4023	t8w, &
4024	mu_3d, &
4025	rho, &
4026	z, &
4027	dz8w, &
4028	z_at_w
4029
4030	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , &
4031	INTENT(IN ) :: pb, &
4032	p, &
4033	u, &
4034	v, &
4035	w, &
4036	alt, &
4037	ph, &
4038	phb, &
4039	t
4040
4041
4042	REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fzm, &
4043	fzp
4044
4045	REAL, DIMENSION( ims:ime , kms:kme, jms:jme ), &
4046	INTENT(INOUT) :: RTHRATEN
4047
4048	REAL, DIMENSION( ims:ime , kms:kme, jms:jme ), &
4049	INTENT(INOUT) :: RTHCUTEN, &
4050	RQVCUTEN, &
4051	RQCCUTEN, &
4052	RQRCUTEN, &
4053	RQICUTEN, &
4054	RQSCUTEN
4055
4056	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) , &
4057	INTENT(INOUT) :: RUBLTEN, &
4058	RVBLTEN, &
4059	RTHBLTEN, &
4060	RQVBLTEN, &
4061	RQCBLTEN, &
4062	RQIBLTEN
4063
4064	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) , &
4065	INTENT(INOUT) :: RTHFTEN, &
4066	RQVFTEN
4067
4068	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) , &
4069	INTENT(INOUT) :: RUNDGDTEN, &
4070	RVNDGDTEN, &
4071	RTHNDGDTEN, &
4072	RQVNDGDTEN, &
4073	RMUNDGDTEN
4074
4075	INTEGER :: i_start, i_end, j_start, j_end, k_start, k_end, i_startu, j_startv
4076	INTEGER :: i, j, k
4077	REAL :: w1, w2, z0, z1, z2
4078
4079	!-----------------------------------------------------------------------
4080
4081	!<DESCRIPTION>
4082	!
4083	! phys_prep calculates a number of diagnostic quantities needed by
4084	! the physics routines. It also decouples the physics tendencies from
4085	! the column dry-air mass (the physics routines expect to see/update the
4086	! uncoupled tendencies).
4087	!
4088	!</DESCRIPTION>
4089
4090	! set up loop bounds for this grid's boundary conditions
4091
4092	i_start = its
4093	i_end = min( ite,ide-1 )
4094	j_start = jts
4095	j_end = min( jte,jde-1 )
4096
4097	k_start = kts
4098	k_end = min( kte, kde-1 )
4099
4100	! compute thermodynamics and velocities at pressure points
4101
4102	do j = j_start,j_end
4103	do k = k_start, k_end
4104	do i = i_start, i_end
4105
4106	th_phy(i,k,j) = t(i,k,j) + t0
4107	p_phy(i,k,j) = p(i,k,j) + pb(i,k,j)
4108	pi_phy(i,k,j) = (p_phy(i,k,j)/p1000mb)**rcp
4109	!! TAKE INTO ACCOUNT cp=f(T) on Venus
4110	IF (planet .eq. "venus" ) THEN
4111	t_phy(i,k,j)= (th_phy(i,k,j)*nu - nu(TT00*nu)log((p1000mb/p_phy(i,k,j))rcp))(1/nu)
4112	ELSE
4113	t_phy(i,k,j) = th_phy(i,k,j)*pi_phy(i,k,j)
4114	ENDIF
4115	rho(i,k,j) = 1./alt(i,k,j)*(1.+moist(i,k,j,P_QV))
4116	mu_3d(i,k,j) = mu(i,j)
4117	u_phy(i,k,j) = 0.5*(u(i,k,j)+u(i+1,k,j))
4118	v_phy(i,k,j) = 0.5*(v(i,k,j)+v(i,k,j+1))
4119
4120	enddo
4121	enddo
4122	enddo
4123
4124	! compute z at w points
4125
4126	do j = j_start,j_end
4127	do k = k_start, kte
4128	do i = i_start, i_end
4129	z_at_w(i,k,j) = (phb(i,k,j)+ph(i,k,j))/g
4130	enddo
4131	enddo
4132	enddo
4133
4134	do j = j_start,j_end
4135	do k = k_start, kte-1
4136	do i = i_start, i_end
4137	dz8w(i,k,j) = z_at_w(i,k+1,j)-z_at_w(i,k,j)
4138	enddo
4139	enddo
4140	enddo
4141
4142	do j = j_start,j_end
4143	do i = i_start, i_end
4144	dz8w(i,kte,j) = 0.
4145	enddo
4146	enddo
4147
4148	! compute z at p points (average of z at w points)
4149
4150	do j = j_start,j_end
4151	do k = k_start, k_end
4152	do i = i_start, i_end
4153	z(i,k,j) = 0.5*(z_at_w(i,k,j) + z_at_w(i,k+1,j) )
4154	!!!! MARS MARS ajout aymeric (ainsi que les arguments de cette routine)
4155	w_phy(i,k,j) = 0.5*(w(i,k,j) + w(i,k+1,j) )
4156	enddo
4157	enddo
4158	enddo
4159
4160	! interp t and p at w points
4161
4162	do j = j_start,j_end
4163	do k = 2, k_end
4164	do i = i_start, i_end
4165	p8w(i,k,j) = fzm(k)p_phy(i,k,j)+fzp(k)p_phy(i,k-1,j)
4166	t8w(i,k,j) = fzm(k)t_phy(i,k,j)+fzp(k)t_phy(i,k-1,j)
4167	enddo
4168	enddo
4169	enddo
4170
4171	! extrapolate p and t to surface and top.
4172	! we'll use an extrapolation in z for now
4173
4174	do j = j_start,j_end
4175	do i = i_start, i_end
4176
4177	! bottom
4178
4179	z0 = z_at_w(i,1,j)
4180	z1 = z(i,1,j)
4181	z2 = z(i,2,j)
4182	w1 = (z0 - z2)/(z1 - z2)
4183	w2 = 1. - w1
4184	p8w(i,1,j) = w1p_phy(i,1,j)+w2p_phy(i,2,j)
4185	t8w(i,1,j) = w1t_phy(i,1,j)+w2t_phy(i,2,j)
4186
4187	! top
4188
4189	z0 = z_at_w(i,kte,j)
4190	z1 = z(i,k_end,j)
4191	z2 = z(i,k_end-1,j)
4192	w1 = (z0 - z2)/(z1 - z2)
4193	w2 = 1. - w1
4194
4195	! p8w(i,kde,j) = w1p_phy(i,kde-1,j)+w2p_phy(i,kde-2,j)
4196	!!! bug fix extrapolate ln(p) so p is positive definite
4197	p8w(i,kde,j) = exp(w1log(p_phy(i,kde-1,j))+w2log(p_phy(i,kde-2,j)))
4198	t8w(i,kde,j) = w1t_phy(i,kde-1,j)+w2t_phy(i,kde-2,j)
4199
4200	enddo
4201	enddo
4202
4203	! decouple all physics tendencies
4204
4205	IF (config_flags%ra_lw_physics .gt. 0 .or. config_flags%ra_sw_physics .gt. 0) THEN
4206
4207	DO J=j_start,j_end
4208	DO K=k_start,k_end
4209	DO I=i_start,i_end
4210	RTHRATEN(I,K,J)=RTHRATEN(I,K,J)/mu(I,J)
4211	ENDDO
4212	ENDDO
4213	ENDDO
4214
4215	ENDIF
4216
4217	IF (config_flags%cu_physics .gt. 0) THEN
4218
4219	DO J=j_start,j_end
4220	DO I=i_start,i_end
4221	DO K=k_start,k_end
4222	RTHCUTEN(I,K,J)=RTHCUTEN(I,K,J)/mu(I,J)
4223	ENDDO
4224	ENDDO
4225	ENDDO
4226
4227	IF (P_QV .ge. PARAM_FIRST_SCALAR)THEN
4228	DO J=j_start,j_end
4229	DO I=i_start,i_end
4230	DO K=k_start,k_end
4231	RQVCUTEN(I,K,J)=RQVCUTEN(I,K,J)/mu(I,J)
4232	ENDDO
4233	ENDDO
4234	ENDDO
4235	ENDIF
4236
4237	IF (P_QC .ge. PARAM_FIRST_SCALAR)THEN
4238	DO J=j_start,j_end
4239	DO I=i_start,i_end
4240	DO K=k_start,k_end
4241	RQCCUTEN(I,K,J)=RQCCUTEN(I,K,J)/mu(I,J)
4242	ENDDO
4243	ENDDO
4244	ENDDO
4245	ENDIF
4246
4247	IF (P_QR .ge. PARAM_FIRST_SCALAR)THEN
4248	DO J=j_start,j_end
4249	DO I=i_start,i_end
4250	DO K=k_start,k_end
4251	RQRCUTEN(I,K,J)=RQRCUTEN(I,K,J)/mu(I,J)
4252	ENDDO
4253	ENDDO
4254	ENDDO
4255	ENDIF
4256
4257	IF (P_QI .ge. PARAM_FIRST_SCALAR)THEN
4258	DO J=j_start,j_end
4259	DO I=i_start,i_end
4260	DO K=k_start,k_end
4261	RQICUTEN(I,K,J)=RQICUTEN(I,K,J)/mu(I,J)
4262	ENDDO
4263	ENDDO
4264	ENDDO
4265	ENDIF
4266
4267	IF(P_QS .ge. PARAM_FIRST_SCALAR)THEN
4268	DO J=j_start,j_end
4269	DO I=i_start,i_end
4270	DO K=k_start,k_end
4271	RQSCUTEN(I,K,J)=RQSCUTEN(I,K,J)/mu(I,J)
4272	ENDDO
4273	ENDDO
4274	ENDDO
4275	ENDIF
4276
4277	ENDIF
4278
4279	IF ( (config_flags%bl_pbl_physics .gt. 0) &
4280	.OR. (config_flags%modif_wrf) ) THEN
4281	!****MARS
4282	DO J=j_start,j_end
4283	DO K=k_start,k_end
4284	DO I=i_start,i_end
4285	RUBLTEN(I,K,J) =RUBLTEN(I,K,J)/mu(I,J)
4286	RVBLTEN(I,K,J) =RVBLTEN(I,K,J)/mu(I,J)
4287	RTHBLTEN(I,K,J)=RTHBLTEN(I,K,J)/mu(I,J)
4288	ENDDO
4289	ENDDO
4290	ENDDO
4291
4292	IF (P_QV .ge. PARAM_FIRST_SCALAR) THEN
4293	DO J=j_start,j_end
4294	DO K=k_start,k_end
4295	DO I=i_start,i_end
4296	RQVBLTEN(I,K,J)=RQVBLTEN(I,K,J)/mu(I,J)
4297	ENDDO
4298	ENDDO
4299	ENDDO
4300	ENDIF
4301
4302	IF (P_QC .ge. PARAM_FIRST_SCALAR) THEN
4303	DO J=j_start,j_end
4304	DO K=k_start,k_end
4305	DO I=i_start,i_end
4306	RQCBLTEN(I,K,J)=RQCBLTEN(I,K,J)/mu(I,J)
4307	ENDDO
4308	ENDDO
4309	ENDDO
4310	ENDIF
4311
4312	IF (P_QI .ge. PARAM_FIRST_SCALAR) THEN
4313	DO J=j_start,j_end
4314	DO K=k_start,k_end
4315	DO I=i_start,i_end
4316	RQIBLTEN(I,K,J)=RQIBLTEN(I,K,J)/mu(I,J)
4317	ENDDO
4318	ENDDO
4319	ENDDO
4320	ENDIF
4321
4322	ENDIF
4323
4324	! decouple advective forcing required by Grell-Devenyi scheme
4325
4326	if ( config_flags%cu_physics == GDSCHEME ) then
4327
4328	DO J=j_start,j_end
4329	DO I=i_start,i_end
4330	DO K=k_start,k_end
4331	RTHFTEN(I,K,J)=RTHFTEN(I,K,J)/mu(I,J)
4332	ENDDO
4333	ENDDO
4334	ENDDO
4335
4336	IF (P_QV .ge. PARAM_FIRST_SCALAR)THEN
4337	DO J=j_start,j_end
4338	DO I=i_start,i_end
4339	DO K=k_start,k_end
4340	RQVFTEN(I,K,J)=RQVFTEN(I,K,J)/mu(I,J)
4341	ENDDO
4342	ENDDO
4343	ENDDO
4344	ENDIF
4345
4346	END IF
4347
4348	! fdda
4349	! note fdda u and v tendencies are staggered, also only interior points have muu/muv,
4350	! so only decouple those
4351
4352	IF (config_flags%grid_fdda .gt. 0) THEN
4353
4354	i_startu=MAX(its,ids+1)
4355	j_startv=MAX(jts,jds+1)
4356
4357	DO J=j_start,j_end
4358	DO K=k_start,k_end
4359	DO I=i_startu,i_end
4360	RUNDGDTEN(I,K,J) =RUNDGDTEN(I,K,J)/muu(I,J)
4361	ENDDO
4362	ENDDO
4363	ENDDO
4364	DO J=j_startv,j_end
4365	DO K=k_start,k_end
4366	DO I=i_start,i_end
4367	RVNDGDTEN(I,K,J) =RVNDGDTEN(I,K,J)/muv(I,J)
4368	ENDDO
4369	ENDDO
4370	ENDDO
4371	DO J=j_start,j_end
4372	DO K=k_start,k_end
4373	DO I=i_start,i_end
4374	RTHNDGDTEN(I,K,J)=RTHNDGDTEN(I,K,J)/mu(I,J)
4375	! RMUNDGDTEN(I,J) - no coupling
4376	ENDDO
4377	ENDDO
4378	ENDDO
4379	IF (P_QV .ge. PARAM_FIRST_SCALAR) THEN
4380	DO J=j_start,j_end
4381	DO K=k_start,k_end
4382	DO I=i_start,i_end
4383	RQVNDGDTEN(I,K,J)=RQVNDGDTEN(I,K,J)/mu(I,J)
4384	ENDDO
4385	ENDDO
4386	ENDDO
4387	ENDIF
4388
4389	ENDIF
4390
4391	END SUBROUTINE phy_prep
4392
4393	!------------------------------------------------------------
4394
4395	SUBROUTINE moist_physics_prep_em( t_new, t_old, t0, rho, al, alb, &
4396	p, p8w, p0, pb, ph, phb, &
4397	th_phy, pii, pf, &
4398	z, z_at_w, dz8w, &
4399	dt,h_diabatic, &
4400	config_flags,fzm, fzp, &
4401	ids,ide, jds,jde, kds,kde, &
4402	ims,ime, jms,jme, kms,kme, &
4403	its,ite, jts,jte, kts,kte )
4404
4405	IMPLICIT NONE
4406
4407	! Here we construct full fields
4408	! needed by the microphysics
4409
4410	TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
4411
4412	INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde
4413	INTEGER, INTENT(IN ) :: ims,ime, jms,jme, kms,kme
4414	INTEGER, INTENT(IN ) :: its,ite, jts,jte, kts,kte
4415
4416	REAL, INTENT(IN ) :: dt
4417
4418	REAL, DIMENSION( ims:ime , kms:kme, jms:jme ), &
4419	INTENT(IN ) :: al, &
4420	alb, &
4421	p, &
4422	pb, &
4423	ph, &
4424	phb
4425
4426
4427	REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fzm, &
4428	fzp
4429
4430	REAL, DIMENSION( ims:ime , kms:kme, jms:jme ), &
4431	INTENT( OUT) :: rho, &
4432	th_phy, &
4433	pii, &
4434	pf, &
4435	z, &
4436	z_at_w, &
4437	dz8w, &
4438	p8w
4439	! pjj/cray
4440	! p8w, &
4441	! h_diabatic
4442
4443	REAL, DIMENSION( ims:ime , kms:kme, jms:jme ), &
4444	INTENT(INOUT) :: h_diabatic
4445
4446	REAL, DIMENSION( ims:ime , kms:kme, jms:jme ), &
4447	INTENT(INOUT) :: t_new, &
4448	t_old
4449
4450	REAL, INTENT(IN ) :: t0, p0
4451	REAL :: z0,z1,z2,w1,w2
4452
4453	INTEGER :: i_start, i_end, j_start, j_end, k_start, k_end
4454	INTEGER :: i, j, k
4455
4456	!--------------------------------------------------------------------
4457
4458	!<DESCRIPTION>
4459	!
4460	! moist_phys_prep_em calculates a number of diagnostic quantities needed by
4461	! the microphysics routines.
4462	!
4463	!</DESCRIPTION>
4464
4465	! set up loop bounds for this grid's boundary conditions
4466
4467	i_start = its
4468	i_end = min( ite,ide-1 )
4469	j_start = jts
4470	j_end = min( jte,jde-1 )
4471
4472	k_start = kts
4473	k_end = min( kte, kde-1 )
4474
4475	DO j = j_start, j_end
4476	DO k = k_start, kte
4477	DO i = i_start, i_end
4478	z_at_w(i,k,j) = (ph(i,k,j)+phb(i,k,j))/g
4479	ENDDO
4480	ENDDO
4481	ENDDO
4482
4483	do j = j_start,j_end
4484	do k = k_start, kte-1
4485	do i = i_start, i_end
4486	dz8w(i,k,j) = z_at_w(i,k+1,j)-z_at_w(i,k,j)
4487	enddo
4488	enddo
4489	enddo
4490
4491	do j = j_start,j_end
4492	do i = i_start, i_end
4493	dz8w(i,kte,j) = 0.
4494	enddo
4495	enddo
4496
4497
4498	! compute full pii, rho, and z at the new time-level
4499	! (needed for physics).
4500	! convert perturbation theta to full theta (th_phy)
4501	! use h_diabatic to temporarily save pre-microphysics full theta
4502
4503	DO j = j_start, j_end
4504	DO k = k_start, k_end
4505	DO i = i_start, i_end
4506
4507	#ifdef REVERT
4508	t_new(i,k,j) = t_new(i,k,j)-h_diabatic(i,k,j)*dt
4509	#endif
4510	th_phy(i,k,j) = t_new(i,k,j) + t0
4511	h_diabatic(i,k,j) = th_phy(i,k,j)
4512	rho(i,k,j) = 1./(al(i,k,j)+alb(i,k,j))
4513	pii(i,k,j) = ((p(i,k,j)+pb(i,k,j))/p0)**rcp
4514	z(i,k,j) = 0.5*(z_at_w(i,k,j) +z_at_w(i,k+1,j) )
4515	pf(i,k,j) = p(i,k,j)+pb(i,k,j)
4516
4517	ENDDO
4518	ENDDO
4519	ENDDO
4520
4521	! interp t and p at w points
4522
4523	do j = j_start,j_end
4524	do k = 2, k_end
4525	do i = i_start, i_end
4526	p8w(i,k,j) = fzm(k)pf(i,k,j)+fzp(k)pf(i,k-1,j)
4527	enddo
4528	enddo
4529	enddo
4530
4531	! extrapolate p and t to surface and top.
4532	! we'll use an extrapolation in z for now
4533
4534	do j = j_start,j_end
4535	do i = i_start, i_end
4536
4537	! bottom
4538
4539	z0 = z_at_w(i,1,j)
4540	z1 = z(i,1,j)
4541	z2 = z(i,2,j)
4542	w1 = (z0 - z2)/(z1 - z2)
4543	w2 = 1. - w1
4544	p8w(i,1,j) = w1pf(i,1,j)+w2pf(i,2,j)
4545
4546	! top
4547
4548	z0 = z_at_w(i,kte,j)
4549	z1 = z(i,k_end,j)
4550	z2 = z(i,k_end-1,j)
4551	w1 = (z0 - z2)/(z1 - z2)
4552	w2 = 1. - w1
4553	! p8w(i,kde,j) = w1pf(i,kde-1,j)+w2pf(i,kde-2,j)
4554	p8w(i,kde,j) = exp(w1log(pf(i,kde-1,j))+w2log(pf(i,kde-2,j)))
4555
4556	enddo
4557	enddo
4558
4559	END SUBROUTINE moist_physics_prep_em
4560
4561	!------------------------------------------------------------------------------
4562
4563	SUBROUTINE moist_physics_finish_em( t_new, t_old, t0, mut, &
4564	th_phy, h_diabatic, dt, &
4565	config_flags, &
4566	ids,ide, jds,jde, kds,kde, &
4567	ims,ime, jms,jme, kms,kme, &
4568	its,ite, jts,jte, kts,kte )
4569
4570	IMPLICIT NONE
4571
4572	! Here we construct full fields
4573	! needed by the microphysics
4574
4575	TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
4576
4577	INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde
4578	INTEGER, INTENT(IN ) :: ims,ime, jms,jme, kms,kme
4579	INTEGER, INTENT(IN ) :: its,ite, jts,jte, kts,kte
4580
4581	REAL, DIMENSION( ims:ime , kms:kme, jms:jme ), &
4582	INTENT(INOUT) :: t_new, &
4583	t_old, &
4584	th_phy, &
4585	h_diabatic
4586
4587	REAL, DIMENSION( ims:ime , jms:jme ), INTENT(INOUT) :: mut
4588
4589
4590	REAL, INTENT(IN ) :: t0, dt
4591
4592	INTEGER :: i_start, i_end, j_start, j_end, k_start, k_end
4593	INTEGER :: i, j, k
4594
4595	!--------------------------------------------------------------------
4596
4597	!<DESCRIPTION>
4598	!
4599	! moist_phys_finish_em resets theta to its perturbation value and
4600	! computes and stores the microphysics diabatic heating term.
4601	!
4602	!</DESCRIPTION>
4603
4604	! set up loop bounds for this grid's boundary conditions
4605
4606
4607	i_start = its
4608	i_end = min( ite,ide-1 )
4609	j_start = jts
4610	j_end = min( jte,jde-1 )
4611
4612	k_start = kts
4613	k_end = min( kte, kde-1 )
4614
4615	! add microphysics theta diff to perturbation theta, set h_diabatic
4616
4617	DO j = j_start, j_end
4618	DO k = k_start, k_end
4619	DO i = i_start, i_end
4620
4621	t_new(i,k,j) = t_new(i,k,j) + (th_phy(i,k,j)-h_diabatic(i,k,j))
4622	h_diabatic(i,k,j) = (th_phy(i,k,j)-h_diabatic(i,k,j))/dt
4623	! h_diabatic(i,k,j) = 0.
4624
4625	ENDDO
4626	ENDDO
4627	ENDDO
4628
4629	END SUBROUTINE moist_physics_finish_em
4630
4631	!----------------------------------------------------------------
4632
4633
4634	SUBROUTINE init_module_big_step
4635	END SUBROUTINE init_module_big_step
4636
4637	SUBROUTINE set_tend ( field, field_adv_tend, msf, &
4638	ids, ide, jds, jde, kds, kde, &
4639	ims, ime, jms, jme, kms, kme, &
4640	its, ite, jts, jte, kts, kte )
4641
4642	IMPLICIT NONE
4643
4644	! Input data
4645
4646	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
4647	ims, ime, jms, jme, kms, kme, &
4648	its, ite, jts, jte, kts, kte
4649
4650	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT) :: field
4651
4652	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN) :: field_adv_tend
4653
4654	REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN) :: msf
4655
4656	! Local data
4657
4658	INTEGER :: i, j, k, itf, jtf, ktf
4659
4660	!<DESCRIPTION>
4661	!
4662	! set_tend copies the advective tendency array into the tendency array.
4663	!
4664	!</DESCRIPTION>
4665
4666	jtf = MIN(jte,jde-1)
4667	ktf = MIN(kte,kde-1)
4668	itf = MIN(ite,ide-1)
4669	DO j = jts, jtf
4670	DO k = kts, ktf
4671	DO i = its, itf
4672	field(i,k,j) = field_adv_tend(i,k,j)*msf(i,j)
4673	ENDDO
4674	ENDDO
4675	ENDDO
4676
4677	END SUBROUTINE set_tend
4678
4679	!------------------------------------------------------------------------------
4680
4681	SUBROUTINE rk_rayleigh_damp( ru_tendf, rv_tendf, &
4682	rw_tendf, t_tendf, &
4683	u, v, w, t, t_init, &
4684	mut, muu, muv, ph, phb, &
4685	u_base, v_base, t_base, z_base, &
4686	dampcoef, zdamp, &
4687	ids, ide, jds, jde, kds, kde, &
4688	ims, ime, jms, jme, kms, kme, &
4689	its, ite, jts, jte, kts, kte )
4690
4691	! History: Apr 2005 Modifications by George Bryan, NCAR:
4692	! - Generalized the code in a way that allows for
4693	! simulations with steep terrain.
4694	!
4695	! Jul 2004 Modifications by George Bryan, NCAR:
4696	! - Modified the code to use u_base, v_base, and t_base
4697	! arrays for the background state. Removed the hard-wired
4698	! base-state values.
4699	! - Modified the code to use dampcoef, zdamp, and damp_opt,
4700	! i.e., the upper-level damper variables in namelist.input.
4701	! Removed the hard-wired variables in the older version.
4702	! This damper is used when damp_opt = 2.
4703	! - Modified the code to account for the movement of the
4704	! model surfaces with time. The code now obtains a base-
4705	! state value by interpolation using the "_base" arrays.
4706
4707	! Nov 2003 Bug fix by Jason Knievel, NCAR
4708
4709	! Aug 2003 Meridional dimension, some comments, and
4710	! changes in layout of the code added by
4711	! Jason Knievel, NCAR
4712
4713	! Jul 2003 Original code by Bill Skamarock, NCAR
4714
4715	! Purpose: This routine applies Rayleigh damping to a layer at top
4716	! of the model domain.
4717
4718	!-----------------------------------------------------------------------
4719	! Begin declarations.
4720
4721	IMPLICIT NONE
4722
4723	INTEGER, INTENT( IN ) &
4724	:: ids, ide, jds, jde, kds, kde, &
4725	ims, ime, jms, jme, kms, kme, &
4726	its, ite, jts, jte, kts, kte
4727
4728	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( INOUT ) &
4729	:: ru_tendf, rv_tendf, rw_tendf, t_tendf
4730
4731	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( IN ) &
4732	:: u, v, w, t, t_init, ph, phb
4733
4734	REAL, DIMENSION( ims:ime, jms:jme ), INTENT( IN ) &
4735	:: mut, muu, muv
4736
4737	REAL, DIMENSION( kms:kme ) , INTENT(IN ) &
4738	:: u_base, v_base, t_base, z_base
4739
4740	REAL, INTENT(IN ) &
4741	:: dampcoef, zdamp
4742
4743	! Local variables.
4744
4745	INTEGER &
4746	:: i_start, i_end, j_start, j_end, k_start, k_end, i, j, k, ktf, k1, k2
4747
4748	REAL &
4749	:: pii, dcoef, z, ztop
4750
4751	REAL :: wkp1, wk, wkm1
4752
4753	REAL, DIMENSION( kms:kme ) :: z00, u00, v00, t00
4754
4755	! End declarations.
4756	!-----------------------------------------------------------------------
4757
4758	pii = 2.0 * asin(1.0)
4759
4760	ktf = MIN( kte, kde-1 )
4761
4762	!-----------------------------------------------------------------------
4763	! Adjust u to base state.
4764
4765	DO j = jts, MIN( jte, jde-1 )
4766	DO i = its, MIN( ite, ide )
4767
4768	! Get height at top of model
4769	ztop = 0.5*( phb(i ,kde,j)+phb(i-1,kde,j) &
4770	+ph(i ,kde,j)+ph(i-1,kde,j) )/g
4771
4772	! Find bottom of damping layer
4773	k1 = ktf
4774	z = ztop
4775	DO WHILE( z >= (ztop-zdamp) )
4776	z = 0.25*( phb(i ,k1,j)+phb(i ,k1+1,j) &
4777	+phb(i-1,k1,j)+phb(i-1,k1+1,j) &
4778	+ph(i ,k1,j)+ph(i ,k1+1,j) &
4779	+ph(i-1,k1,j)+ph(i-1,k1+1,j))/g
4780	z00(k1) = z
4781	k1 = k1 - 1
4782	ENDDO
4783	k1 = k1 + 2
4784
4785	! Get reference state at model levels
4786	DO k = k1, ktf
4787	k2 = ktf
4788	DO WHILE( z_base(k2) .gt. z00(k) )
4789	k2 = k2 - 1
4790	ENDDO
4791	if(k2+1.gt.ktf)then
4792	u00(k) = u_base(k2) + ( u_base(k2) - u_base(k2-1) ) &
4793	* ( z00(k) - z_base(k2) ) &
4794	/ ( z_base(k2) - z_base(k2-1) )
4795	else
4796	u00(k) = u_base(k2) + ( u_base(k2+1) - u_base(k2) ) &
4797	* ( z00(k) - z_base(k2) ) &
4798	/ ( z_base(k2+1) - z_base(k2) )
4799	endif
4800	ENDDO
4801
4802	! Apply the Rayleigh damper
4803	DO k = k1, ktf
4804	dcoef = 1.0 - MIN( 1.0, ( ztop - z00(k) ) / zdamp )
4805	dcoef = (SIN( 0.5 * pii * dcoef ) )**2
4806	ru_tendf(i,k,j) = ru_tendf(i,k,j) - &
4807	muu(i,j) * ( dcoef * dampcoef ) * &
4808	( u(i,k,j) - u00(k) )
4809	END DO
4810
4811	END DO
4812	END DO
4813
4814	! End adjustment of u.
4815	!-----------------------------------------------------------------------
4816
4817	!-----------------------------------------------------------------------
4818	! Adjust v to base state.
4819
4820	DO j = jts, MIN( jte, jde )
4821	DO i = its, MIN( ite, ide-1 )
4822
4823	! Get height at top of model
4824	ztop = 0.5*( phb(i,kde,j )+phb(i,kde,j-1) &
4825	+ph(i,kde,j )+ph(i,kde,j-1) )/g
4826
4827	! Find bottom of damping layer
4828	k1 = ktf
4829	z = ztop
4830	DO WHILE( z >= (ztop-zdamp) )
4831	z = 0.25*( phb(i,k1,j )+phb(i,k1+1,j ) &
4832	+phb(i,k1,j-1)+phb(i,k1+1,j-1) &
4833	+ph(i,k1,j )+ph(i,k1+1,j ) &
4834	+ph(i,k1,j-1)+ph(i,k1+1,j-1))/g
4835	z00(k1) = z
4836	k1 = k1 - 1
4837	ENDDO
4838	k1 = k1 + 2
4839
4840	! Get reference state at model levels
4841	DO k = k1, ktf
4842	k2 = ktf
4843	DO WHILE( z_base(k2) .gt. z00(k) )
4844	k2 = k2 - 1
4845	ENDDO
4846	if(k2+1.gt.ktf)then
4847	v00(k) = v_base(k2) + ( v_base(k2) - v_base(k2-1) ) &
4848	* ( z00(k) - z_base(k2) ) &
4849	/ ( z_base(k2) - z_base(k2-1) )
4850	else
4851	v00(k) = v_base(k2) + ( v_base(k2+1) - v_base(k2) ) &
4852	* ( z00(k) - z_base(k2) ) &
4853	/ ( z_base(k2+1) - z_base(k2) )
4854	endif
4855	ENDDO
4856
4857	! Apply the Rayleigh damper
4858	DO k = k1, ktf
4859	dcoef = 1.0 - MIN( 1.0, ( ztop - z00(k) ) / zdamp )
4860	dcoef = (SIN( 0.5 * pii * dcoef ) )**2
4861	rv_tendf(i,k,j) = rv_tendf(i,k,j) - &
4862	muv(i,j) * ( dcoef * dampcoef ) * &
4863	( v(i,k,j) - v00(k) )
4864	END DO
4865
4866	END DO
4867	END DO
4868
4869	! End adjustment of v.
4870	!-----------------------------------------------------------------------
4871
4872	!-----------------------------------------------------------------------
4873	! Adjust w to base state.
4874
4875	DO j = jts, MIN( jte, jde-1 )
4876	DO i = its, MIN( ite, ide-1 )
4877	ztop = ( phb(i,kde,j) + ph(i,kde,j) ) / g
4878	DO k = kts, MIN( kte, kde )
4879	z = ( phb(i,k,j) + ph(i,k,j) ) / g
4880	IF ( z >= (ztop-zdamp) ) THEN
4881	dcoef = 1.0 - MIN( 1.0, ( ztop - z ) / zdamp )
4882	dcoef = ( SIN( 0.5 * pii * dcoef ) )**2
4883	rw_tendf(i,k,j) = rw_tendf(i,k,j) - &
4884	mut(i,j) * ( dcoef * dampcoef ) * w(i,k,j)
4885	END IF
4886	END DO
4887	END DO
4888	END DO
4889
4890	! End adjustment of w.
4891	!-----------------------------------------------------------------------
4892
4893	!-----------------------------------------------------------------------
4894	! Adjust potential temperature to base state.
4895
4896	DO j = jts, MIN( jte, jde-1 )
4897	DO i = its, MIN( ite, ide-1 )
4898
4899	! Get height at top of model
4900	ztop = ( phb(i,kde,j) + ph(i,kde,j) ) / g
4901
4902	! Find bottom of damping layer
4903	k1 = ktf
4904	z = ztop
4905	DO WHILE( z >= (ztop-zdamp) )
4906	z = 0.5 * ( phb(i,k1,j) + phb(i,k1+1,j) + &
4907	ph(i,k1,j) + ph(i,k1+1,j) ) / g
4908	z00(k1) = z
4909	k1 = k1 - 1
4910	ENDDO
4911	k1 = k1 + 2
4912
4913	! Get reference state at model levels
4914	DO k = k1, ktf
4915	k2 = ktf
4916	DO WHILE( z_base(k2) .gt. z00(k) )
4917	k2 = k2 - 1
4918	ENDDO
4919	if(k2+1.gt.ktf)then
4920	t00(k) = t_base(k2) + ( t_base(k2) - t_base(k2-1) ) &
4921	* ( z00(k) - z_base(k2) ) &
4922	/ ( z_base(k2) - z_base(k2-1) )
4923	else
4924	t00(k) = t_base(k2) + ( t_base(k2+1) - t_base(k2) ) &
4925	* ( z00(k) - z_base(k2) ) &
4926	/ ( z_base(k2+1) - z_base(k2) )
4927	endif
4928	ENDDO
4929
4930	! Apply the Rayleigh damper
4931	DO k = k1, ktf
4932	dcoef = 1.0 - MIN( 1.0, ( ztop - z00(k) ) / zdamp )
4933	dcoef = (SIN( 0.5 * pii * dcoef ) )**2
4934	t_tendf(i,k,j) = t_tendf(i,k,j) - &
4935	mut(i,j) * ( dcoef * dampcoef ) * &
4936	( t(i,k,j) - t00(k) )
4937	END DO
4938
4939	END DO
4940	END DO
4941
4942	! End adjustment of potential temperature.
4943	!-----------------------------------------------------------------------
4944
4945	END SUBROUTINE rk_rayleigh_damp
4946
4947	!==============================================================================
4948	!==============================================================================
4949
4950	SUBROUTINE sixth_order_diffusion( name, field, tendency, mu, dt, &
4951	config_flags, &
4952	diff_6th_opt, diff_6th_factor, &
4953	ids, ide, jds, jde, kds, kde, &
4954	ims, ime, jms, jme, kms, kme, &
4955	its, ite, jts, jte, kts, kte )
4956
4957	! History: 14 Nov 2006 Name of variable changed by Jason Knievel
4958	! 07 Jun 2006 Revised and generalized by Jason Knievel
4959	! 25 Apr 2005 Original code by Jason Knievel, NCAR
4960
4961	! Purpose: Apply 6th-order, monotonic (flux-limited), numerical
4962	! diffusion to 3-d velocity and to scalars.
4963
4964	! References: Ming Xue (MWR Aug 2000)
4965	! Durran ("Numerical Methods for Wave Equations..." 1999)
4966	! George Bryan (personal communication)
4967
4968	!------------------------------------------------------------------------------
4969	! Begin: Declarations.
4970
4971	IMPLICIT NONE
4972
4973	INTEGER, INTENT(IN) &
4974	:: ids, ide, jds, jde, kds, kde, &
4975	ims, ime, jms, jme, kms, kme, &
4976	its, ite, jts, jte, kts, kte
4977
4978	TYPE(grid_config_rec_type), INTENT(IN) &
4979	:: config_flags
4980
4981	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) &
4982	:: tendency
4983
4984	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN) &
4985	:: field
4986
4987	REAL, DIMENSION( ims:ime , jms:jme ), INTENT(IN) &
4988	:: mu
4989
4990	REAL, INTENT(IN) &
4991	:: dt
4992
4993	REAL, INTENT(IN) &
4994	:: diff_6th_factor
4995
4996	INTEGER, INTENT(IN) &
4997	:: diff_6th_opt
4998
4999	CHARACTER(LEN=1) , INTENT(IN) &
5000	:: name
5001
5002	INTEGER &
5003	:: i, j, k, &
5004	i_start, i_end, &
5005	j_start, j_end, &
5006	k_start, k_end, &
5007	ktf
5008
5009	REAL &
5010	:: dflux_x_p0, dflux_y_p0, &
5011	dflux_x_p1, dflux_y_p1, &
5012	tendency_x, tendency_y, &
5013	mu_avg_p0, mu_avg_p1, &
5014	diff_6th_coef
5015
5016	LOGICAL &
5017	:: specified
5018
5019	! End: Declarations.
5020	!------------------------------------------------------------------------------
5021
5022	!------------------------------------------------------------------------------
5023	! Begin: Translate the diffusion factor into a diffusion coefficient. See
5024	! Durran's text, section 2.4.3, then adjust for sixth-order diffusion (not
5025	! fourth) and for diffusion in two dimensions (not one). For reference, a
5026	! factor of 1.0 would mean complete diffusion of a 2dx wave in one time step,
5027	! although application of the flux limiter reduces somewhat the effects of
5028	! diffusion for a given coefficient.
5029
5030	diff_6th_coef = diff_6th_factor * 0.015625 / ( 2.0 * dt )
5031
5032	! End: Translate diffusion factor.
5033	!------------------------------------------------------------------------------
5034
5035	!------------------------------------------------------------------------------
5036	! Begin: Assign limits of spatial loops depending on variable to be diffused.
5037	! The halo regions are already filled with values by the time this subroutine
5038	! is called, which allows the stencil to extend beyond the domains' edges.
5039
5040	ktf = MIN( kte, kde-1 )
5041
5042	IF ( name .EQ. 'u' ) THEN
5043
5044	i_start = its
5045	i_end = ite
5046	j_start = jts
5047	j_end = MIN(jde-1,jte)
5048	k_start = kts
5049	k_end = ktf
5050
5051	ELSE IF ( name .EQ. 'v' ) THEN
5052
5053	i_start = its
5054	i_end = MIN(ide-1,ite)
5055	j_start = jts
5056	j_end = jte
5057	k_start = kts
5058	k_end = ktf
5059
5060	ELSE IF ( name .EQ. 'w' ) THEN
5061
5062	i_start = its
5063	i_end = MIN(ide-1,ite)
5064	j_start = jts
5065	j_end = MIN(jde-1,jte)
5066	k_start = kts+1
5067	k_end = ktf
5068
5069	ELSE
5070
5071	i_start = its
5072	i_end = MIN(ide-1,ite)
5073	j_start = jts
5074	j_end = MIN(jde-1,jte)
5075	k_start = kts
5076	k_end = ktf
5077
5078	ENDIF
5079
5080	! End: Assignment of limits of spatial loops.
5081	!------------------------------------------------------------------------------
5082
5083	!------------------------------------------------------------------------------
5084	! Begin: Loop across spatial dimensions.
5085
5086	DO j = j_start, j_end
5087	DO k = k_start, k_end
5088	DO i = i_start, i_end
5089
5090	!------------------------------------------------------------------------------
5091	! Begin: Diffusion in x (i index).
5092
5093	! Calculate the diffusive flux in x direction (from Xue's eq. 3).
5094
5095	dflux_x_p0 = ( 10.0 * ( field(i, k,j) - field(i-1,k,j) ) &
5096	- 5.0 * ( field(i+1,k,j) - field(i-2,k,j) ) &
5097	+ ( field(i+2,k,j) - field(i-3,k,j) ) )
5098
5099	dflux_x_p1 = ( 10.0 * ( field(i+1,k,j) - field(i ,k,j) ) &
5100	- 5.0 * ( field(i+2,k,j) - field(i-1,k,j) ) &
5101	+ ( field(i+3,k,j) - field(i-2,k,j) ) )
5102
5103	! If requested in the namelist (diff_6th_opt=2), prohibit up-gradient diffusion
5104	! (variation on Xue's eq. 10).
5105
5106	IF ( diff_6th_opt .EQ. 2 ) THEN
5107
5108	IF ( dflux_x_p0 * ( field(i ,k,j)-field(i-1,k,j) ) .LE. 0.0 ) THEN
5109	dflux_x_p0 = 0.0
5110	END IF
5111
5112	IF ( dflux_x_p1 * ( field(i+1,k,j)-field(i ,k,j) ) .LE. 0.0 ) THEN
5113	dflux_x_p1 = 0.0
5114	END IF
5115
5116	END IF
5117
5118	! Apply 6th-order diffusion in x direction.
5119
5120	IF ( name .EQ. 'u' ) THEN
5121	mu_avg_p0 = mu(i-1,j)
5122	mu_avg_p1 = mu(i ,j)
5123	ELSE IF ( name .EQ. 'v' ) THEN
5124	mu_avg_p0 = 0.25 * ( &
5125	mu(i-1,j-1) + &
5126	mu(i ,j-1) + &
5127	mu(i-1,j ) + &
5128	mu(i ,j ) )
5129	mu_avg_p1 = 0.25 * ( &
5130	mu(i ,j-1) + &
5131	mu(i+1,j-1) + &
5132	mu(i ,j ) + &
5133	mu(i+1,j ) )
5134	ELSE
5135	mu_avg_p0 = 0.5 * ( &
5136	mu(i-1,j) + &
5137	mu(i ,j) )
5138	mu_avg_p1 = 0.5 * ( &
5139	mu(i ,j) + &
5140	mu(i+1,j) )
5141	END IF
5142
5143	tendency_x = diff_6th_coef * &
5144	( ( mu_avg_p1 * dflux_x_p1 ) - ( mu_avg_p0 * dflux_x_p0 ) )
5145
5146	! End: Diffusion in x.
5147	!------------------------------------------------------------------------------
5148
5149	!------------------------------------------------------------------------------
5150	! Begin: Diffusion in y (j index).
5151
5152	! Calculate the diffusive flux in y direction (from Xue's eq. 3).
5153
5154	dflux_y_p0 = ( 10.0 * ( field(i,k,j ) - field(i,k,j-1) ) &
5155	- 5.0 * ( field(i,k,j+1) - field(i,k,j-2) ) &
5156	+ ( field(i,k,j+2) - field(i,k,j-3) ) )
5157
5158	dflux_y_p1 = ( 10.0 * ( field(i,k,j+1) - field(i,k,j ) ) &
5159	- 5.0 * ( field(i,k,j+2) - field(i,k,j-1) ) &
5160	+ ( field(i,k,j+3) - field(i,k,j-2) ) )
5161
5162	! If requested in the namelist (diff_6th_opt=2), prohibit up-gradient diffusion
5163	! (variation on Xue's eq. 10).
5164
5165	IF ( diff_6th_opt .EQ. 2 ) THEN
5166
5167	IF ( dflux_y_p0 * ( field(i,k,j )-field(i,k,j-1) ) .LE. 0.0 ) THEN
5168	dflux_y_p0 = 0.0
5169	END IF
5170
5171	IF ( dflux_y_p1 * ( field(i,k,j+1)-field(i,k,j ) ) .LE. 0.0 ) THEN
5172	dflux_y_p1 = 0.0
5173	END IF
5174
5175	END IF
5176
5177	! Apply 6th-order diffusion in y direction.
5178
5179	IF ( name .EQ. 'u' ) THEN
5180	mu_avg_p0 = 0.25 * ( &
5181	mu(i-1,j-1) + &
5182	mu(i ,j-1) + &
5183	mu(i-1,j ) + &
5184	mu(i ,j ) )
5185	mu_avg_p1 = 0.25 * ( &
5186	mu(i-1,j ) + &
5187	mu(i ,j ) + &
5188	mu(i-1,j+1) + &
5189	mu(i ,j+1) )
5190	ELSE IF ( name .EQ. 'v' ) THEN
5191	mu_avg_p0 = mu(i,j-1)
5192	mu_avg_p1 = mu(i,j )
5193	ELSE
5194	mu_avg_p0 = 0.5 * ( &
5195	mu(i,j-1) + &
5196	mu(i,j ) )
5197	mu_avg_p1 = 0.5 * ( &
5198	mu(i,j ) + &
5199	mu(i,j+1) )
5200	END IF
5201
5202	tendency_y = diff_6th_coef * &
5203	( ( mu_avg_p1 * dflux_y_p1 ) - ( mu_avg_p0 * dflux_y_p0 ) )
5204
5205	! End: Diffusion in y.
5206	!------------------------------------------------------------------------------
5207
5208	!------------------------------------------------------------------------------
5209	! Begin: Combine diffusion in x and y.
5210
5211	tendency(i,k,j) = tendency(i,k,j) + tendency_x + tendency_y
5212
5213	! End: Combine diffusion in x and y.
5214	!------------------------------------------------------------------------------
5215
5216	ENDDO
5217	ENDDO
5218	ENDDO
5219
5220	! End: Loop across spatial dimensions.
5221	!------------------------------------------------------------------------------
5222
5223	END SUBROUTINE sixth_order_diffusion
5224
5225	!==============================================================================
5226	!==============================================================================
5227
5228	END MODULE module_big_step_utilities_em

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

Download in other formats:

Original Format