Context Navigation

module_big_step_utilities_em.F

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

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

WRF: version v3.3
LMDZ: version v1818

More details in:

File size: 205.4 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, ONLY : domain
16	USE module_model_constants
17	USE module_state_description
18	USE module_configure, ONLY : grid_config_rec_type
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	! Map scale factor comments for this routine:
325	! Locally not changed, but sent the correct map scale factors
326	! from module_em (msfuy, msfvx, msfty)
327
328	SUBROUTINE couple_momentum ( muu, ru, u, msfu, &
329	muv, rv, v, msfv, msfv_inv, &
330	mut, rw, w, msft, &
331	ids, ide, jds, jde, kds, kde, &
332	ims, ime, jms, jme, kms, kme, &
333	its, ite, jts, jte, kts, kte )
334
335	IMPLICIT NONE
336
337	! Input data
338
339	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
340	ims, ime, jms, jme, kms, kme, &
341	its, ite, jts, jte, kts, kte
342
343	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT( OUT) :: ru, rv, rw
344
345	REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: muu, muv, mut
346	REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: msfu, msfv, msft, msfv_inv
347
348	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: u, v, w
349
350	! Local data
351
352	INTEGER :: i, j, k, itf, jtf, ktf
353
354	!<DESCRIPTION>
355	!
356	! couple_momentum couples the velocities to the full column mass and
357	! the map factors.
358	!
359	!</DESCRIPTION>
360
361	ktf=MIN(kte,kde-1)
362
363	itf=ite
364	jtf=MIN(jte,jde-1)
365
366	DO j=jts,jtf
367	DO k=kts,ktf
368	DO i=its,itf
369	ru(i,k,j)=u(i,k,j)*muu(i,j)/msfu(i,j)
370	ENDDO
371	ENDDO
372	ENDDO
373
374	itf=MIN(ite,ide-1)
375	jtf=jte
376
377	DO j=jts,jtf
378	DO k=kts,ktf
379	DO i=its,itf
380	rv(i,k,j)=v(i,k,j)muv(i,j)msfv_inv(i,j)
381	! rv(i,k,j)=v(i,k,j)*muv(i,j)/msfv(i,j)
382	ENDDO
383	ENDDO
384	ENDDO
385
386	itf=MIN(ite,ide-1)
387	jtf=MIN(jte,jde-1)
388
389	DO j=jts,jtf
390	DO k=kts,kte
391	DO i=its,itf
392	rw(i,k,j)=w(i,k,j)*mut(i,j)/msft(i,j)
393	ENDDO
394	ENDDO
395	ENDDO
396
397	END SUBROUTINE couple_momentum
398
399	!-------------------------------------------------------------------
400
401	SUBROUTINE calc_mu_staggered ( mu, mub, muu, muv, &
402	ids, ide, jds, jde, kds, kde, &
403	ims, ime, jms, jme, kms, kme, &
404	its, ite, jts, jte, kts, kte )
405
406	IMPLICIT NONE
407
408	! Input data
409
410	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
411	ims, ime, jms, jme, kms, kme, &
412	its, ite, jts, jte, kts, kte
413
414	REAL, DIMENSION( ims:ime , jms:jme ) , INTENT( OUT) :: muu, muv
415	REAL, DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mu, mub
416
417	! local stuff
418
419	INTEGER :: i, j, itf, jtf
420
421	!<DESCRIPTION>
422	!
423	! calc_mu_staggered calculates the full dry air mass at the staggered
424	! velocity points (u,v).
425	!
426	!</DESCRIPTION>
427
428	itf=ite
429	jtf=MIN(jte,jde-1)
430
431	IF ( ( its .NE. ids ) .AND. ( ite .NE. ide ) ) THEN
432	DO j=jts,jtf
433	DO i=its,itf
434	muu(i,j) = 0.5*(mu(i,j)+mu(i-1,j)+mub(i,j)+mub(i-1,j))
435	ENDDO
436	ENDDO
437	ELSE IF ( ( its .EQ. ids ) .AND. ( ite .NE. ide ) ) THEN
438	DO j=jts,jtf
439	DO i=its+1,itf
440	muu(i,j) = 0.5*(mu(i,j)+mu(i-1,j)+mub(i,j)+mub(i-1,j))
441	ENDDO
442	ENDDO
443	i=its
444	DO j=jts,jtf
445	muu(i,j) = mu(i,j) +mub(i,j)
446	ENDDO
447	ELSE IF ( ( its .NE. ids ) .AND. ( ite .EQ. ide ) ) THEN
448	DO j=jts,jtf
449	DO i=its,itf-1
450	muu(i,j) = 0.5*(mu(i,j)+mu(i-1,j)+mub(i,j)+mub(i-1,j))
451	ENDDO
452	ENDDO
453	i=ite
454	DO j=jts,jtf
455	muu(i,j) = mu(i-1,j) +mub(i-1,j)
456	ENDDO
457	ELSE IF ( ( its .EQ. ids ) .AND. ( ite .EQ. ide ) ) THEN
458	DO j=jts,jtf
459	DO i=its+1,itf-1
460	muu(i,j) = 0.5*(mu(i,j)+mu(i-1,j)+mub(i,j)+mub(i-1,j))
461	ENDDO
462	ENDDO
463	i=its
464	DO j=jts,jtf
465	muu(i,j) = mu(i,j) +mub(i,j)
466	ENDDO
467	i=ite
468	DO j=jts,jtf
469	muu(i,j) = mu(i-1,j) +mub(i-1,j)
470	ENDDO
471	END IF
472
473	itf=MIN(ite,ide-1)
474	jtf=jte
475
476	IF ( ( jts .NE. jds ) .AND. ( jte .NE. jde ) ) THEN
477	DO j=jts,jtf
478	DO i=its,itf
479	muv(i,j) = 0.5*(mu(i,j)+mu(i,j-1)+mub(i,j)+mub(i,j-1))
480	ENDDO
481	ENDDO
482	ELSE IF ( ( jts .EQ. jds ) .AND. ( jte .NE. jde ) ) THEN
483	DO j=jts+1,jtf
484	DO i=its,itf
485	muv(i,j) = 0.5*(mu(i,j)+mu(i,j-1)+mub(i,j)+mub(i,j-1))
486	ENDDO
487	ENDDO
488	j=jts
489	DO i=its,itf
490	muv(i,j) = mu(i,j) +mub(i,j)
491	ENDDO
492	ELSE IF ( ( jts .NE. jds ) .AND. ( jte .EQ. jde ) ) THEN
493	DO j=jts,jtf-1
494	DO i=its,itf
495	muv(i,j) = 0.5*(mu(i,j)+mu(i,j-1)+mub(i,j)+mub(i,j-1))
496	ENDDO
497	ENDDO
498	j=jte
499	DO i=its,itf
500	muv(i,j) = mu(i,j-1) +mub(i,j-1)
501	ENDDO
502	ELSE IF ( ( jts .EQ. jds ) .AND. ( jte .EQ. jde ) ) THEN
503	DO j=jts+1,jtf-1
504	DO i=its,itf
505	muv(i,j) = 0.5*(mu(i,j)+mu(i,j-1)+mub(i,j)+mub(i,j-1))
506	ENDDO
507	ENDDO
508	j=jts
509	DO i=its,itf
510	muv(i,j) = mu(i,j) +mub(i,j)
511	ENDDO
512	j=jte
513	DO i=its,itf
514	muv(i,j) = mu(i,j-1) +mub(i,j-1)
515	ENDDO
516	END IF
517
518	END SUBROUTINE calc_mu_staggered
519
520	!-------------------------------------------------------------------------------
521
522	SUBROUTINE couple ( mu, mub, rfield, field, name, &
523	msf, &
524	ids, ide, jds, jde, kds, kde, &
525	ims, ime, jms, jme, kms, kme, &
526	its, ite, jts, jte, kts, kte )
527
528	IMPLICIT NONE
529
530	! Input data
531
532	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
533	ims, ime, jms, jme, kms, kme, &
534	its, ite, jts, jte, kts, kte
535
536	CHARACTER(LEN=1) , INTENT(IN ) :: name
537
538	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT( OUT) :: rfield
539
540	REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mu, mub, msf
541
542	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: field
543
544	! Local data
545
546	INTEGER :: i, j, k, itf, jtf, ktf
547	REAL , DIMENSION(ims:ime,jms:jme) :: muu , muv
548
549	!<DESCRIPTION>
550	!
551	! subroutine couple couples the input variable with the dry-air
552	! column mass (mu).
553	!
554	!</DESCRIPTION>
555
556
557	ktf=MIN(kte,kde-1)
558
559	IF (name .EQ. 'u')THEN
560
561	CALL calc_mu_staggered ( mu, mub, muu, muv, &
562	ids, ide, jds, jde, kds, kde, &
563	ims, ime, jms, jme, kms, kme, &
564	its, ite, jts, jte, kts, kte )
565
566	itf=ite
567	jtf=MIN(jte,jde-1)
568
569	DO j=jts,jtf
570	DO k=kts,ktf
571	DO i=its,itf
572	rfield(i,k,j)=field(i,k,j)*muu(i,j)/msf(i,j)
573	ENDDO
574	ENDDO
575	ENDDO
576
577	ELSE IF (name .EQ. 'v')THEN
578
579	CALL calc_mu_staggered ( mu, mub, muu, muv, &
580	ids, ide, jds, jde, kds, kde, &
581	ims, ime, jms, jme, kms, kme, &
582	its, ite, jts, jte, kts, kte )
583
584	itf=ite
585	itf=MIN(ite,ide-1)
586	jtf=jte
587
588	DO j=jts,jtf
589	DO k=kts,ktf
590	DO i=its,itf
591	rfield(i,k,j)=field(i,k,j)*muv(i,j)/msf(i,j)
592	ENDDO
593	ENDDO
594	ENDDO
595
596	ELSE IF (name .EQ. 'w')THEN
597	itf=MIN(ite,ide-1)
598	jtf=MIN(jte,jde-1)
599	DO j=jts,jtf
600	DO k=kts,kte
601	DO i=its,itf
602	rfield(i,k,j)=field(i,k,j)*(mu(i,j)+mub(i,j))/msf(i,j)
603	ENDDO
604	ENDDO
605	ENDDO
606
607	ELSE IF (name .EQ. 'h')THEN
608	itf=MIN(ite,ide-1)
609	jtf=MIN(jte,jde-1)
610	DO j=jts,jtf
611	DO k=kts,kte
612	DO i=its,itf
613	rfield(i,k,j)=field(i,k,j)*(mu(i,j)+mub(i,j))
614	ENDDO
615	ENDDO
616	ENDDO
617
618	ELSE
619	itf=MIN(ite,ide-1)
620	jtf=MIN(jte,jde-1)
621	DO j=jts,jtf
622	DO k=kts,ktf
623	DO i=its,itf
624	rfield(i,k,j)=field(i,k,j)*(mu(i,j)+mub(i,j))
625	ENDDO
626	ENDDO
627	ENDDO
628
629	ENDIF
630
631	END SUBROUTINE couple
632
633
634	!-------------------------------------------------------------------------------
635
636	SUBROUTINE calc_ww_cp ( u, v, mup, mub, ww, &
637	rdx, rdy, msftx, msfty, &
638	msfux, msfuy, msfvx, msfvx_inv, &
639	msfvy, dnw, &
640	ids, ide, jds, jde, kds, kde, &
641	ims, ime, jms, jme, kms, kme, &
642	its, ite, jts, jte, kts, kte )
643
644	IMPLICIT NONE
645
646	! Input data
647
648
649	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
650	ims, ime, jms, jme, kms, kme, &
651	its, ite, jts, jte, kts, kte
652
653	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: u, v
654	REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mup, mub, &
655	msftx, msfty, &
656	msfux, msfuy, &
657	msfvx, msfvy, &
658	msfvx_inv
659	REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: dnw
660
661	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT ) :: ww
662	REAL , INTENT(IN ) :: rdx, rdy
663
664	! Local data
665
666	INTEGER :: i, j, k, itf, jtf, ktf
667	REAL , DIMENSION( its:ite ) :: dmdt
668	REAL , DIMENSION( its:ite, kts:kte ) :: divv
669	REAL , DIMENSION( its:ite+1, jts:jte+1 ) :: muu, muv
670
671	!<DESCRIPTION>
672	!
673	! calc_ww calculates omega using the velocities (u,v) and the dry-air
674	! column mass (mup+mub).
675	! The algorithm integrates the continuity equation through the column
676	! followed by a diagnosis of omega.
677	!
678	!</DESCRIPTION>
679
680	!<DESCRIPTION>
681	!
682	! calc_ww_cp calculates omega using the velocities (u,v) and the
683	! column mass mu.
684	!
685	!</DESCRIPTION>
686
687	jtf=MIN(jte,jde-1)
688	ktf=MIN(kte,kde-1)
689	itf=MIN(ite,ide-1)
690
691	! mu coupled with the appropriate map factor
692
693	DO j=jts,jtf
694	DO i=its,min(ite+1,ide)
695	! u is always coupled with my
696	muu(i,j) = 0.5*(mup(i,j)+mub(i,j)+mup(i-1,j)+mub(i-1,j))/msfuy(i,j)
697	ENDDO
698	ENDDO
699
700	DO j=jts,min(jte+1,jde)
701	DO i=its,itf
702	! v is always coupled with mx
703	! muv(i,j) = 0.5*(mup(i,j)+mub(i,j)+mup(i,j-1)+mub(i,j-1))/msfvx(i,j)
704	muv(i,j) = 0.5(mup(i,j)+mub(i,j)+mup(i,j-1)+mub(i,j-1))msfvx_inv(i,j)
705	ENDDO
706	ENDDO
707
708	DO j=jts,jtf
709
710	DO i=its,ite
711	dmdt(i) = 0.
712	ww(i,1,j) = 0.
713	ww(i,kte,j) = 0.
714	ENDDO
715
716	! Comments on the modifications for map scale factors
717	! ADT eqn 47 / my (putting rho -> 'mu') is:
718	! (1/my) partial d mu/dt = -mx partial d/dx(mu u/my)
719	! -mx partial d/dy(mu v/mx)
720	! -partial d/dz(mu w/my)
721	!
722	! Using nu instead of z the last term becomes:
723	! -partial d/dnu(mu (dnu/dt)/my)
724	!
725	! Integrating with respect to nu over ALL levels, with dnu/dt=0 at top
726	! and bottom, the last term becomes = 0
727	!
728	! Integral\|bot->top[(1/my) partial d mu/dt]dnu =
729	! Integral\|bot->top[-mx partial d/dx(mu u/my)
730	! -mx partial d/dy(mu v/mx)]dnu
731	!
732	! muu='mu'[on u]/my, muv='mu'[on v]/mx
733	! (1/my) partial d mu/dt is independent of nu
734	! => LHS = Integral\|bot->top[con]dnu = conservation*(-1) = -dmdt
735	!
736	! => dmdt = mx*Integral\|bot->top[partial d/dx(mu u/my) +
737	! partial d/dy(mu v/mx)]dnu
738	! => dmdt = sum_bot->top[divv]
739	! where
740	! divv=mx[partial d/dx(mu u/my) + partial d/dy(mu v/mx)]delta nu
741
742	DO k=kts,ktf
743	DO i=its,itf
744
745	divv(i,k) = msftx(i,j)dnw(k)( rdx(muu(i+1,j)u(i+1,k,j)-muu(i,j)*u(i,k,j)) &
746	+rdy(muv(i,j+1)v(i,k,j+1)-muv(i,j)*v(i,k,j)) )
747
748	! dmdt(i) = dmdt(i) + dnw(k)* ( rdx*(ru(i+1,k,j)-ru(i,k,j)) &
749	! +rdy*(rv(i,k,j+1)-rv(i,k,j)) )
750
751	dmdt(i) = dmdt(i) + divv(i,k)
752
753
754	ENDDO
755	ENDDO
756
757	! Further map scale factor notes:
758	! Now integrate from bottom to top, level by level:
759	! mu dnu/dt/my [k+1] = mu dnu/dt/my [k] + [-(1/my) partial d mu/dt
760	! -mx partial d/dx(mu u/my)
761	! -mx partial d/dy(mu v/mx)]*dnu[k->k+1]
762	! ww [k+1] = ww [k] -(1/my) partial d mu/dt * dnu[k->k+1] - divv[k]
763	! = ww [k] -dmdt * dnw[k] - divv[k]
764
765	DO k=2,ktf
766	DO i=its,itf
767
768	! ww(i,k,j)=ww(i,k-1,j) &
769	! - dnw(k-1)* ( dmdt(i) &
770	! +rdx*(ru(i+1,k-1,j)-ru(i,k-1,j)) &
771	! +rdy*(rv(i,k-1,j+1)-rv(i,k-1,j)) )
772
773	ww(i,k,j)=ww(i,k-1,j) - dnw(k-1)*dmdt(i) - divv(i,k-1)
774
775	ENDDO
776	ENDDO
777	ENDDO
778
779
780	END SUBROUTINE calc_ww_cp
781
782
783	!-------------------------------------------------------------------------------
784
785	SUBROUTINE calc_cq ( moist, cqu, cqv, cqw, n_moist, &
786	ids, ide, jds, jde, kds, kde, &
787	ims, ime, jms, jme, kms, kme, &
788	its, ite, jts, jte, kts, kte )
789
790	IMPLICIT NONE
791
792	! Input data
793
794	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
795	ims, ime, jms, jme, kms, kme, &
796	its, ite, jts, jte, kts, kte
797
798	INTEGER , INTENT(IN ) :: n_moist
799
800
801	REAL, DIMENSION( ims:ime, kms:kme , jms:jme , n_moist ), INTENT(IN ) :: moist
802
803	REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT( OUT) :: cqu, cqv, cqw
804
805	! Local stuff
806
807	REAL :: qtot
808
809	INTEGER :: i, j, k, itf, jtf, ktf, ispe
810
811	!<DESCRIPTION>
812	!
813	! calc_cq calculates moist coefficients for the momentum equations.
814	!
815	!</DESCRIPTION>
816
817	itf=ite
818	jtf=MIN(jte,jde-1)
819	ktf=MIN(kte,kde-1)
820
821	IF( n_moist >= PARAM_FIRST_SCALAR ) THEN
822
823	DO j=jts,jtf
824	DO k=kts,ktf
825	DO i=its,itf
826	qtot = 0.
827	!DEC$ loop count(3)
828	DO ispe=PARAM_FIRST_SCALAR,n_moist
829	qtot = qtot + moist(i,k,j,ispe) + moist(i-1,k,j,ispe)
830	ENDDO
831	! qtot = 0.5*( moist(i ,k,j,1)+moist(i ,k,j,2)+moist(i ,k,j,3)+ &
832	! & moist(i-1,k,j,1)+moist(i-1,k,j,2)+moist(i-1,k,j,3) )
833	! cqu(i,k,j) = 1./(1.+qtot)
834	cqu(i,k,j) = 1./(1.+0.5*qtot)
835	ENDDO
836	ENDDO
837	ENDDO
838
839	itf=MIN(ite,ide-1)
840	jtf=jte
841
842	DO j=jts,jtf
843	DO k=kts,ktf
844	DO i=its,itf
845	qtot = 0.
846	!DEC$ loop count(3)
847	DO ispe=PARAM_FIRST_SCALAR,n_moist
848	qtot = qtot + moist(i,k,j,ispe) + moist(i,k,j-1,ispe)
849	ENDDO
850	! qtot = 0.5*( moist(i,k,j ,1)+moist(i,k,j ,2)+moist(i,k,j ,3)+ &
851	! & moist(i,k,j-1,1)+moist(i,k,j-1,2)+moist(i,k,j-1,3) )
852	! cqv(i,k,j) = 1./(1.+qtot)
853	cqv(i,k,j) = 1./(1.+0.5*qtot)
854	ENDDO
855	ENDDO
856	ENDDO
857
858	itf=MIN(ite,ide-1)
859	jtf=MIN(jte,jde-1)
860	DO j=jts,jtf
861	DO k=kts+1,ktf
862	DO i=its,itf
863	qtot = 0.
864	!DEC$ loop count(3)
865	DO ispe=PARAM_FIRST_SCALAR,n_moist
866	qtot = qtot + moist(i,k,j,ispe) + moist(i,k-1,j,ispe)
867	ENDDO
868	! qtot = 0.5*( moist(i,k ,j,1)+moist(i,k ,j,2)+moist(i,k-1,j,3)+ &
869	! & moist(i,k-1,j,1)+moist(i,k-1,j,2)+moist(i,k ,j,3) )
870	! cqw(i,k,j) = qtot
871	cqw(i,k,j) = 0.5*qtot
872	ENDDO
873	ENDDO
874	ENDDO
875
876	ELSE
877
878	DO j=jts,jtf
879	DO k=kts,ktf
880	DO i=its,itf
881	cqu(i,k,j) = 1.
882	ENDDO
883	ENDDO
884	ENDDO
885
886	itf=MIN(ite,ide-1)
887	jtf=jte
888
889	DO j=jts,jtf
890	DO k=kts,ktf
891	DO i=its,itf
892	cqv(i,k,j) = 1.
893	ENDDO
894	ENDDO
895	ENDDO
896
897	itf=MIN(ite,ide-1)
898	jtf=MIN(jte,jde-1)
899	DO j=jts,jtf
900	DO k=kts+1,ktf
901	DO i=its,itf
902	cqw(i,k,j) = 0.
903	ENDDO
904	ENDDO
905	ENDDO
906
907	END IF
908
909	END SUBROUTINE calc_cq
910
911	!----------------------------------------------------------------------
912
913	SUBROUTINE calc_alt ( alt, al, alb, &
914	ids, ide, jds, jde, kds, kde, &
915	ims, ime, jms, jme, kms, kme, &
916	its, ite, jts, jte, kts, kte )
917
918	IMPLICIT NONE
919
920	! Input data
921
922	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
923	ims, ime, jms, jme, kms, kme, &
924	its, ite, jts, jte, kts, kte
925
926	REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(IN ) :: alb, al
927	REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT( OUT) :: alt
928
929	! Local stuff
930
931	INTEGER :: i, j, k, itf, jtf, ktf
932
933	!<DESCRIPTION>
934	!
935	! calc_alt computes the full inverse density
936	!
937	!</DESCRIPTION>
938
939	itf=MIN(ite,ide-1)
940	jtf=MIN(jte,jde-1)
941	ktf=MIN(kte,kde-1)
942
943	DO j=jts,jtf
944	DO k=kts,ktf
945	DO i=its,itf
946	alt(i,k,j) = al(i,k,j)+alb(i,k,j)
947	ENDDO
948	ENDDO
949	ENDDO
950
951
952	END SUBROUTINE calc_alt
953
954	!----------------------------------------------------------------------
955
956	SUBROUTINE calc_p_rho_phi ( moist, n_moist, &
957	al, alb, mu, muts, ph, p, pb, &
958	t, p0, t0, znu, dnw, rdnw, &
959	rdn, non_hydrostatic, &
960	ids, ide, jds, jde, kds, kde, &
961	ims, ime, jms, jme, kms, kme, &
962	its, ite, jts, jte, kts, kte )
963
964	IMPLICIT NONE
965
966	! Input data
967
968	LOGICAL , INTENT(IN ) :: non_hydrostatic
969
970	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
971	ims, ime, jms, jme, kms, kme, &
972	its, ite, jts, jte, kts, kte
973
974	INTEGER , INTENT(IN ) :: n_moist
975
976	REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(IN ) :: alb, &
977	pb, &
978	t
979
980	REAL, DIMENSION( ims:ime , kms:kme , jms:jme, n_moist ), INTENT(IN ) :: moist
981
982	REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT( OUT) :: al, p
983
984	REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(INOUT) :: ph
985
986	REAL, DIMENSION( ims:ime , jms:jme ), INTENT(IN ) :: mu, muts
987
988	REAL, DIMENSION( kms:kme ), INTENT(IN ) :: znu, dnw, rdnw, rdn
989
990	REAL, INTENT(IN ) :: t0, p0
991
992	! Local stuff
993
994	INTEGER :: i, j, k, itf, jtf, ktf, ispe
995	REAL :: qvf, qtot, qf1, qf2
996	REAL, DIMENSION( its:ite) :: temp,cpovcv_v
997
998
999	!<DESCRIPTION>
1000	!
1001	! For the nonhydrostatic option, calc_p_rho_phi calculates the
1002	! diagnostic quantities pressure and (inverse) density from the
1003	! prognostic variables using the equation of state.
1004	!
1005	! For the hydrostatic option, calc_p_rho_phi calculates the
1006	! diagnostic quantities (inverse) density and geopotential from the
1007	! prognostic variables using the equation of state and the hydrostatic
1008	! equation.
1009	!
1010	!</DESCRIPTION>
1011	#ifdef LMDZ1
1012	INTEGER :: ix,im2,km2,jm2
1013
1014	im2=24
1015	jm2=21
1016	km2=38
1017
1018	PRINT *,' calc_p_rho_phi inside: ',im2,km2,jm2
1019	PRINT *,' al: ',al(im2,km2,jm2), ' p: ', p(im2,km2,jm2), &
1020	' ph: ',ph(im2,km2,jm2), ' mu: ', mu(im2,jm2), ' muts: ', muts(im2,jm2), &
1021	' t: ',t(im2,km2,jm2), ' p0: ',p0, ' t0: ', t0, 'znu: ',znu(km2), &
1022	' dnw: ',dnw(km2), ' rdnw: ',rdnw(km2), ' rdn: ', rdn(km2), ' moist: ', &
1023	moist(im2,km2,jm2,:)
1024	#endif
1025
1026	itf=MIN(ite,ide-1)
1027	jtf=MIN(jte,jde-1)
1028	ktf=MIN(kte,kde-1)
1029
1030	#ifndef INTELMKL
1031	cpovcv_v = cpovcv
1032	#endif
1033
1034	IF (non_hydrostatic) THEN
1035
1036	IF (n_moist >= PARAM_FIRST_SCALAR ) THEN
1037
1038	DO j=jts,jtf
1039	DO k=kts,ktf
1040	DO i=its,itf
1041	qvf = 1.+rvovrd*moist(i,k,j,P_QV)
1042	al(i,k,j)=-1./muts(i,j)(alb(i,k,j)mu(i,j) &
1043	+rdnw(k)*(ph(i,k+1,j)-ph(i,k,j)))
1044	temp(i)=(r_d(t0+t(i,k,j))qvf)/ &
1045	(p0*(al(i,k,j)+alb(i,k,j)))
1046	#ifdef LMDZ1
1047	IF (i == im2 .AND. k==km2 .AND. j==jm2) THEN
1048	PRINT *,' calc_p_rho_phi after non_hydrostatic: ',im2,km2,jm2
1049	PRINT *,' al: ',al(im2,km2,jm2), ' p: ', p(im2,km2,jm2), &
1050	' ph: ',ph(im2,km2,jm2), ' mu: ', mu(im2,jm2), ' muts: ', muts(im2,jm2), &
1051	' t: ',t(im2,km2,jm2), ' p0: ',p0, ' t0: ', t0, 'znu: ',znu(km2), &
1052	' dnw: ',dnw(km2), ' rdnw: ',rdnw(km2), ' rdn: ', rdn(km2), ' qvf: ', &
1053	moist(im2,km2,jm2,p_qv)
1054	END IF
1055	#endif
1056	ENDDO
1057	#ifdef LMDZ1
1058	IF (k==km2 .AND. j==jm2) THEN
1059	PRINT *,' calc_p_rho_phi before VPOW: ',im2,km2,jm2
1060	PRINT *,' p sfc: ', p(its,km2,jm2), ' temp: ',temp(its), ' cpovcv_v: ',&
1061	cpovcv_v(its),' p: ', p(im2,km2,jm2)
1062	END IF
1063	#endif
1064	#ifdef INTELMKL
1065	CALL VPOWX ( itf-its+1, temp(its), cpovcv, p(its,k,j) )
1066	#else
1067	! use vector version from libmassv or from compat lib in frame/libmassv.F
1068	! L. Fita, LMD April 2014. On the writting of a wrfout, this gives a NaN value at a
1069	! certain vertical level not at the its !!!!
1070	CALL VPOW ( p(its,k,j), temp(its), cpovcv_v(its), itf-its+1 )
1071	! DO ix=its,itf-its+1
1072	! p(ix,k,j) = temp(ix)**cpovcv_v(ix)
1073	! IF (k==km2 .AND. j==jm2) PRINT *,ix, ': ', p(ix,k,j), temp(ix), cpovcv_v(ix)
1074	! END DO
1075	#endif
1076	#ifdef LMDZ1
1077	IF (k==km2 .AND. j==jm2) THEN
1078	PRINT *,' calc_p_rho_phi after VPOW: ',im2,km2,jm2
1079	PRINT *,' p sfc: ', p(its,km2,jm2), ' temp: ',temp(its), ' cpovcv_v: ', &
1080	cpovcv_v(its), ' pow: ',itf-its+1,' p: ', p(im2,km2,jm2)
1081	PRINT *,' calc_p_rho_phi before p_compute: ',im2,km2,jm2
1082	PRINT *,' p: ', p(im2,km2,jm2), ' p0: ',p0, ' pb: ', pb(im2,km2,jm2)
1083	END IF
1084	#endif
1085	DO i=its,itf
1086	p(i,k,j)= p(i,k,j)*p0-pb(i,k,j)
1087	#ifdef LMDZ1
1088	IF (i == im2 .AND. k==km2 .AND. j==jm2) THEN
1089	PRINT *,' calc_p_rho_phi after p_compute: ',im2,km2,jm2
1090	PRINT *,' p: ', p(im2,km2,jm2), ' p0: ',p0, ' pb: ', pb(im2,km2,jm2)
1091	END IF
1092	#endif
1093	ENDDO
1094	ENDDO
1095	ENDDO
1096
1097	ELSE
1098
1099	DO j=jts,jtf
1100	DO k=kts,ktf
1101	DO i=its,itf
1102	al(i,k,j)=-1./muts(i,j)(alb(i,k,j)mu(i,j) &
1103	+rdnw(k)*(ph(i,k+1,j)-ph(i,k,j)))
1104	p(i,k,j)=p0( (r_d(t0+t(i,k,j)))/ &
1105	(p0(al(i,k,j)+alb(i,k,j))) )*cpovcv &
1106	-pb(i,k,j)
1107	#ifdef LMDZ1
1108	IF (i == im2 .AND. k==km2 .AND. j==jm2) THEN
1109	PRINT *,' calc_p_rho_phi after al_p_compute 2: ',im2,km2,jm2
1110	PRINT *,' p: ', p(im2,km2,jm2), ' al: ',al(im2,km2,jm2)
1111	END IF
1112	#endif
1113	ENDDO
1114	ENDDO
1115	ENDDO
1116
1117	END IF
1118
1119	ELSE
1120
1121	! hydrostatic pressure, al, and ph1 calc; WCS, 5 sept 2001
1122
1123
1124	IF (n_moist >= PARAM_FIRST_SCALAR ) THEN
1125
1126	DO j=jts,jtf
1127
1128	k=ktf ! top layer
1129	DO i=its,itf
1130
1131	qtot = 0.
1132	DO ispe=PARAM_FIRST_SCALAR,n_moist
1133	qtot = qtot + moist(i,k,j,ispe)
1134	ENDDO
1135	qf2 = 1.
1136	qf1 = qtot*qf2
1137
1138	p(i,k,j) = - 0.5(mu(i,j)+qf1muts(i,j))/rdnw(k)/qf2
1139	qvf = 1.+rvovrd*moist(i,k,j,P_QV)
1140	al(i,k,j) = (r_d/p1000mb)(t(i,k,j)+t0)qvf* &
1141	(((p(i,k,j)+pb(i,k,j))/p1000mb)**cvpm) - alb(i,k,j)
1142
1143	ENDDO
1144
1145	DO k=ktf-1,kts,-1 ! remaining layers, integrate down
1146	DO i=its,itf
1147
1148	qtot = 0.
1149	DO ispe=PARAM_FIRST_SCALAR,n_moist
1150	qtot = qtot + 0.5*( moist(i,k ,j,ispe) + moist(i,k+1,j,ispe) )
1151	ENDDO
1152	qf2 = 1.
1153	qf1 = qtot*qf2
1154
1155	p(i,k,j) = p(i,k+1,j) - (mu(i,j) + qf1*muts(i,j))/qf2/rdn(k+1)
1156	qvf = 1.+rvovrd*moist(i,k,j,P_QV)
1157	al(i,k,j) = (r_d/p1000mb)(t(i,k,j)+t0)qvf* &
1158	(((p(i,k,j)+pb(i,k,j))/p1000mb)**cvpm) - alb(i,k,j)
1159	ENDDO
1160	ENDDO
1161
1162	DO k=2,ktf+1 ! integrate hydrostatic equation for geopotential
1163	DO i=its,itf
1164
1165	! ph(i,k,j) = ph(i,k-1,j) - (1./rdnw(k-1))*( &
1166	! (muts(i,j)+mu(i,j))*al(i,k-1,j)+ &
1167	! mu(i,j)*alb(i,k-1,j) )
1168	ph(i,k,j) = ph(i,k-1,j) - (dnw(k-1))*( &
1169	(muts(i,j))*al(i,k-1,j)+ &
1170	mu(i,j)*alb(i,k-1,j) )
1171
1172
1173	ENDDO
1174	ENDDO
1175
1176	ENDDO
1177
1178	ELSE
1179
1180	DO j=jts,jtf
1181
1182	k=ktf ! top layer
1183	DO i=its,itf
1184
1185	qtot = 0.
1186	qf2 = 1.
1187	qf1 = qtot*qf2
1188
1189	p(i,k,j) = - 0.5(mu(i,j)+qf1muts(i,j))/rdnw(k)/qf2
1190	qvf = 1.
1191	al(i,k,j) = (r_d/p1000mb)(t(i,k,j)+t0)qvf* &
1192	(((p(i,k,j)+pb(i,k,j))/p1000mb)**cvpm) - alb(i,k,j)
1193
1194	#ifdef LMDZ1
1195	IF (i == im2 .AND. j==jm2) THEN
1196	PRINT *,' calc_p_rho_phi after top_layer: ',im2,jm2
1197	PRINT *,' p: ', p(im2,k,jm2), ' al: ',al(im2,k,jm2)
1198	END IF
1199	#endif
1200	ENDDO
1201
1202	DO k=ktf-1,kts,-1 ! remaining layers, integrate down
1203	DO i=its,itf
1204
1205	qtot = 0.
1206	qf2 = 1.
1207	qf1 = qtot*qf2
1208
1209	p(i,k,j) = p(i,k+1,j) - (mu(i,j) + qf1*muts(i,j))/qf2/rdn(k+1)
1210	qvf = 1.
1211	al(i,k,j) = (r_d/p1000mb)(t(i,k,j)+t0)qvf* &
1212	(((p(i,k,j)+pb(i,k,j))/p1000mb)**cvpm) - alb(i,k,j)
1213	#ifdef LMDZ1
1214	IF (i == im2 .AND. k==km2 .AND. j==jm2) THEN
1215	PRINT *,' calc_p_rho_phi after other_layers: ',im2,km2,jm2
1216	PRINT *,' p: ', p(im2,km2,jm2), ' al: ',al(im2,km2,jm2)
1217	END IF
1218	#endif
1219	ENDDO
1220	ENDDO
1221
1222	DO k=2,ktf+1 ! integrate hydrostatic equation for geopotential
1223	DO i=its,itf
1224
1225	! ph(i,k,j) = ph(i,k-1,j) - (1./rdnw(k-1))*( &
1226	! (muts(i,j)+mu(i,j))*al(i,k-1,j)+ &
1227	! mu(i,j)*alb(i,k-1,j) )
1228	ph(i,k,j) = ph(i,k-1,j) - (dnw(k-1))*( &
1229	(muts(i,j))*al(i,k-1,j)+ &
1230	mu(i,j)*alb(i,k-1,j) )
1231
1232
1233	ENDDO
1234	ENDDO
1235
1236	ENDDO
1237
1238	END IF
1239
1240	END IF
1241
1242	END SUBROUTINE calc_p_rho_phi
1243
1244	!----------------------------------------------------------------------
1245
1246	SUBROUTINE calc_php ( php, ph, phb, &
1247	ids, ide, jds, jde, kds, kde, &
1248	ims, ime, jms, jme, kms, kme, &
1249	its, ite, jts, jte, kts, kte )
1250
1251	IMPLICIT NONE
1252
1253	! Input data
1254
1255	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
1256	ims, ime, jms, jme, kms, kme, &
1257	its, ite, jts, jte, kts, kte
1258
1259	REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(IN ) :: phb, ph
1260	REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT( OUT) :: php
1261
1262	! Local stuff
1263
1264	INTEGER :: i, j, k, itf, jtf, ktf
1265
1266	!<DESCRIPTION>
1267	!
1268	! calc_php calculates the full geopotential from the reference state
1269	! geopotential and the perturbation geopotential (phb_ph).
1270	!
1271	!</DESCRIPTION>
1272
1273	itf=MIN(ite,ide-1)
1274	jtf=MIN(jte,jde-1)
1275	ktf=MIN(kte,kde-1)
1276
1277	DO j=jts,jtf
1278	DO k=kts,ktf
1279	DO i=its,itf
1280	php(i,k,j) = 0.5*(phb(i,k,j)+phb(i,k+1,j)+ph(i,k,j)+ph(i,k+1,j))
1281	ENDDO
1282	ENDDO
1283	ENDDO
1284
1285	END SUBROUTINE calc_php
1286
1287	!-------------------------------------------------------------------------------
1288
1289	SUBROUTINE diagnose_w( ph_tend, ph_new, ph_old, w, mu, dt, &
1290	u, v, ht, &
1291	cf1, cf2, cf3, rdx, rdy, &
1292	msftx, msfty, &
1293	ids, ide, jds, jde, kds, kde, &
1294	ims, ime, jms, jme, kms, kme, &
1295	its, ite, jts, jte, kts, kte )
1296
1297	IMPLICIT NONE
1298
1299	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
1300	ims, ime, jms, jme, kms, kme, &
1301	its, ite, jts, jte, kts, kte
1302
1303	REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(IN ) :: ph_tend, &
1304	ph_new, &
1305	ph_old, &
1306	u, &
1307	v
1308
1309
1310	REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT( OUT) :: w
1311
1312	REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN ) :: mu, ht, msftx, msfty
1313
1314	REAL, INTENT(IN ) :: dt, cf1, cf2, cf3, rdx, rdy
1315
1316	INTEGER :: i, j, k, itf, jtf
1317
1318	itf=MIN(ite,ide-1)
1319	jtf=MIN(jte,jde-1)
1320
1321	!<DESCRIPTION>
1322	!
1323	! diagnose_w diagnoses the vertical velocity from the geopoential equation.
1324	! Used with the hydrostatic option.
1325	!
1326	!</DESCRIPTION>
1327
1328	DO j = jts, jtf
1329
1330	! lower b.c. on w
1331
1332	! Notes on map scale factors:
1333	! Chain rule: if Z=Z(X,Y) [true at the surface] then
1334	! dZ/dt = dZ/dX * dX/dt + dZ/dY * dY/dt, U=dX/dt, V=dY/dt
1335	! Using capitals to denote actual values
1336	! In mapped values, u=U, v=V, z=Z, 1/dX=mx/dx, 1/dY=my/dy
1337	! => w = dz/dt = mx u dz/dx + my v dz/dy
1338	! [where dz/dx is just the surface height change between x
1339	! gridpoints, and dz/dy is the change between y gridpoints]
1340	! [NB: cf1, cf2 and cf3 do vertical weighting of u or v values
1341	! nearest the surface]
1342
1343	! Previously msft multiplied by rdy and rdx terms.
1344	! Now msfty multiplies rdy term, and msftx multiplies msftx term
1345	DO i = its, itf
1346	w(i,1,j)= msfty(i,j).5rdy*( &
1347	(ht(i,j+1)-ht(i,j )) &
1348	(cf1v(i,1,j+1)+cf2v(i,2,j+1)+cf3v(i,3,j+1)) &
1349	+(ht(i,j )-ht(i,j-1)) &
1350	(cf1v(i,1,j )+cf2v(i,2,j )+cf3v(i,3,j )) ) &
1351	+msftx(i,j).5rdx*( &
1352	(ht(i+1,j)-ht(i,j )) &
1353	(cf1u(i+1,1,j)+cf2u(i+1,2,j)+cf3u(i+1,3,j)) &
1354	+(ht(i,j )-ht(i-1,j)) &
1355	(cf1u(i ,1,j)+cf2u(i ,2,j)+cf3u(i ,3,j)) )
1356	ENDDO
1357
1358	! use geopotential equation to diagnose w
1359
1360	! Further notes on map scale factors
1361	! If ph_tend contains: -mx partial d/dx(mu rho u/my)
1362	! -mx partial d/dy(phi mu v/mx)
1363	! -partial d/dz(phi mu w/my)
1364	! then phi eqn is: partial d/dt(mu phi/my) = ph_tend + mu g w/my
1365	! => w = [my/(mug)][partial d/dt(mu phi/my) - ph_tend]
1366
1367	DO k = 2, kte
1368	DO i = its, itf
1369	w(i,k,j) = msfty(i,j)*( (ph_new(i,k,j)-ph_old(i,k,j))/dt &
1370	- ph_tend(i,k,j)/mu(i,j) )/g
1371
1372	ENDDO
1373	ENDDO
1374
1375	ENDDO
1376
1377	END SUBROUTINE diagnose_w
1378
1379	!-------------------------------------------------------------------------------
1380
1381	SUBROUTINE rhs_ph( ph_tend, u, v, ww, &
1382	ph, ph_old, phb, w, &
1383	mut, muu, muv, &
1384	fnm, fnp, &
1385	rdnw, cfn, cfn1, rdx, rdy, &
1386	msfux, msfuy, msfvx, &
1387	msfvx_inv, msfvy, &
1388	msftx, msfty, &
1389	non_hydrostatic, &
1390	config_flags, &
1391	ids, ide, jds, jde, kds, kde, &
1392	ims, ime, jms, jme, kms, kme, &
1393	its, ite, jts, jte, kts, kte )
1394	IMPLICIT NONE
1395
1396	TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
1397
1398	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
1399	ims, ime, jms, jme, kms, kme, &
1400	its, ite, jts, jte, kts, kte
1401
1402	REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(IN ) :: &
1403	u, &
1404	v, &
1405	ww, &
1406	ph, &
1407	ph_old, &
1408	phb, &
1409	w
1410
1411	REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(INOUT) :: ph_tend
1412
1413	REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN ) :: muu, muv, mut, &
1414	msfux, msfuy, &
1415	msfvx, msfvy, &
1416	msftx, msfty, &
1417	msfvx_inv
1418
1419	REAL, DIMENSION( kms:kme ), INTENT(IN ) :: rdnw, fnm, fnp
1420
1421	REAL, INTENT(IN ) :: cfn, cfn1, rdx, rdy
1422
1423	LOGICAL, INTENT(IN ) :: non_hydrostatic
1424
1425	! Local stuff
1426
1427	INTEGER :: i, j, k, itf, jtf, ktf, kz, i_start, j_start
1428	REAL :: ur, ul, ub, vr, vl, vb
1429	REAL, DIMENSION(its:ite,kts:kte) :: wdwn
1430
1431	INTEGER :: advective_order
1432
1433	LOGICAL :: specified
1434
1435	!<DESCRIPTION>
1436	!
1437	! rhs_ph calculates the large-timestep tendency terms for the geopotential
1438	! equation. These terms include the advection and "gw". The geopotential
1439	! equation is cast in advective form, so we don't use the flux form advection
1440	! algorithms here.
1441	!
1442	!</DESCRIPTION>
1443
1444	specified = .false.
1445	if(config_flags%specified .or. config_flags%nested) specified = .true.
1446
1447	advective_order = config_flags%h_sca_adv_order
1448
1449	itf=MIN(ite,ide-1)
1450	jtf=MIN(jte,jde-1)
1451	ktf=MIN(kte,kde-1)
1452
1453	! Notes on map scale factors (WCS, 2 march 2008)
1454	! phi equation is: mu/my d/dt(phi) = -(1/my) mx mu u d/dx(phi)
1455	! -(1/my) my mu v d/dy(phi)
1456	! - omega d/d_eta(phi)
1457	! +mu g w/my
1458	!
1459	! A little further explanation...
1460	! The tendency term we are computing here is for mu/my d/dt(phi). It is advective form
1461	! but it is multiplied be mu/my. It will be decoupled from (mu/my) when the implicit w-phi
1462	! solution is computed in subourine advance_w. The formulation dates from the early
1463	! days of the mass coordinate model when we were testing both a flux and an advective formulation
1464	! for the geopotential equation and different forms of the vertical momentum equation and the
1465	! vertically implicit solver.
1466
1467	! advective form for the geopotential equation
1468
1469	DO j = jts, jtf
1470
1471	DO k = 2, kte
1472	DO i = its, itf
1473	wdwn(i,k) = .5(ww(i,k,j)+ww(i,k-1,j))rdnw(k-1) &
1474	*(ph(i,k,j)-ph(i,k-1,j)+phb(i,k,j)-phb(i,k-1,j))
1475	ENDDO
1476	ENDDO
1477
1478	! RHS term 3 is: - omega partial d/dnu(phi)
1479
1480	DO k = 2, kte-1
1481	DO i = its, itf
1482	ph_tend(i,k,j) = ph_tend(i,k,j) &
1483	- (fnm(k)wdwn(i,k+1)+fnp(k)wdwn(i,k))
1484	ENDDO
1485	ENDDO
1486
1487	ENDDO
1488
1489	IF (non_hydrostatic) THEN ! add in "gw" term.
1490	DO j = jts, jtf ! in hydrostatic mode, "gw" will be diagnosed
1491	! after the timestep to give us "w"
1492	DO i = its, itf
1493	ph_tend(i,kde,j) = 0.
1494	ENDDO
1495
1496	DO k = 2, kte
1497	DO i = its, itf
1498	! phi equation RHS term 4
1499	ph_tend(i,k,j) = ph_tend(i,k,j) + mut(i,j)gw(i,k,j)/msfty(i,j)
1500	ENDDO
1501	ENDDO
1502
1503	ENDDO
1504
1505	END IF
1506
1507	! Notes on map scale factors:
1508	! RHS terms 1 and 2 are: -(1/my) mx u mu partial d/dx(phi)
1509	! -(1/my) my v mu partial d/dy(phi)
1510
1511	IF (advective_order <= 2) THEN
1512
1513	! y (v) advection
1514
1515	i_start = its
1516	j_start = jts
1517	itf=MIN(ite,ide-1)
1518	jtf=MIN(jte,jde-1)
1519
1520	IF ( (config_flags%open_ys .or. specified) .and. jts == jds ) j_start = jts+1
1521	IF ( (config_flags%open_ye .or. specified) .and. jte == jde ) jtf = jtf-2
1522
1523	DO j = j_start, jtf
1524
1525	DO k = 2, kte-1
1526	DO i = i_start, itf
1527	ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25rdy/msfty(i,j)) &
1528	( muv(i,j+1)(v(i,k,j+1)+v(i,k-1,j+1))msfvy(i,j+1)* &
1529	(phb(i,k,j+1)-phb(i,k,j )+ph(i,k,j+1)-ph(i,k,j )) &
1530	+muv(i,j )(v(i,k,j )+v(i,k-1,j ))msfvy(i,j )* &
1531	(phb(i,k,j )-phb(i,k,j-1)+ph(i,k,j )-ph(i,k,j-1)) )
1532	ENDDO
1533	ENDDO
1534
1535	k = kte
1536	DO i = i_start, itf
1537	ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5rdy/msfty(i,j)) &
1538	( muv(i,j+1)(cfnv(i,k-1,j+1)+cfn1v(i,k-2,j+1))msfvy(i,j+1)* &
1539	(phb(i,k,j+1)-phb(i,k,j )+ph(i,k,j+1)-ph(i,k,j )) &
1540	+muv(i,j )(cfnv(i,k-1,j )+cfn1v(i,k-2,j ))msfvy(i,j )* &
1541	(phb(i,k,j )-phb(i,k,j-1)+ph(i,k,j )-ph(i,k,j-1)) )
1542	ENDDO
1543
1544	ENDDO
1545
1546	! x (u) advection
1547
1548	i_start = its
1549	j_start = jts
1550	itf=MIN(ite,ide-1)
1551	jtf=MIN(jte,jde-1)
1552
1553	IF ( (config_flags%open_xs .or. specified) .and. its == ids ) i_start = its+1
1554	IF ( (config_flags%open_xe .or. specified) .and. ite == ide ) itf = itf-2
1555
1556	DO j = j_start, jtf
1557
1558	DO k = 2, kte-1
1559	DO i = i_start, itf
1560	ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25rdx/msfty(i,j)) &
1561	( muu(i+1,j)(u(i+1,k,j)+u(i+1,k-1,j))msfux(i+1,j)* &
1562	(phb(i+1,k,j)-phb(i ,k,j)+ph(i+1,k,j)-ph(i ,k,j)) &
1563	+muu(i ,j)(u(i ,k,j)+u(i ,k-1,j))msfux(i ,j)* &
1564	(phb(i ,k,j)-phb(i-1,k,j)+ph(i ,k,j)-ph(i-1,k,j)) )
1565	ENDDO
1566	ENDDO
1567
1568	k = kte
1569	DO i = i_start, itf
1570	ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5rdx/msfty(i,j)) &
1571	( muu(i+1,j)(cfnu(i+1,k-1,j)+cfn1u(i+1,k-2,j))msfux(i+1,j)* &
1572	(phb(i+1,k,j)-phb(i ,k,j)+ph(i+1,k,j)-ph(i ,k,j)) &
1573	+muu(i ,j)(cfnu(i ,k-1,j)+cfn1u(i ,k-2,j))msfux( i,j)* &
1574	(phb(i ,k,j)-phb(i-1,k,j)+ph(i ,k,j)-ph(i-1,k,j)) )
1575	ENDDO
1576
1577	ENDDO
1578
1579	ELSE IF (advective_order <= 4) THEN
1580
1581	! y (v) advection
1582
1583	i_start = its
1584	j_start = jts
1585	itf=MIN(ite,ide-1)
1586	jtf=MIN(jte,jde-1)
1587
1588	IF ( (config_flags%open_ys .or. specified) .and. jts == jds ) j_start = jts+2
1589	IF ( (config_flags%open_ye .or. specified) .and. jte == jde ) jtf = jtf-3
1590
1591	DO j = j_start, jtf
1592
1593	DO k = 2, kte-1
1594	DO i = i_start, itf
1595	ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25rdy/msfty(i,j))( &
1596	( muv(i,j+1)(v(i,k,j+1)+v(i,k-1,j+1))msfvy(i,j+1) &
1597	+muv(i,j )(v(i,k,j )+v(i,k-1,j ))msfvy(i,j ))* (1./12.)*( &
1598	8.*(ph(i,k,j+1)-ph(i,k,j-1)) &
1599	-(ph(i,k,j+2)-ph(i,k,j-2)) &
1600	+8.*(phb(i,k,j+1)-phb(i,k,j-1)) &
1601	-(phb(i,k,j+2)-phb(i,k,j-2)) ) )
1602
1603
1604	ENDDO
1605	ENDDO
1606
1607	k = kte
1608	DO i = i_start, itf
1609	ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5rdy/msfty(i,j))( &
1610	( muv(i,j+1)(cfnv(i,k-1,j+1)+cfn1v(i,k-2,j+1))msfvy(i,j+1) &
1611	+muv(i,j )(cfnv(i,k-1,j )+cfn1v(i,k-2,j ))msfvy(i,j ))* (1./12.)*( &
1612	8.*(ph(i,k,j+1)-ph(i,k,j-1)) &
1613	-(ph(i,k,j+2)-ph(i,k,j-2)) &
1614	+8.*(phb(i,k,j+1)-phb(i,k,j-1)) &
1615	-(phb(i,k,j+2)-phb(i,k,j-2)) ) )
1616
1617	ENDDO
1618
1619	ENDDO
1620
1621	! pick up near boundary rows using 2nd order stencil for open and specified conditions
1622
1623	IF ( (config_flags%open_ys .or. specified) .and. jts <= jds+1 ) THEN
1624
1625	j = jds+1
1626	DO k = 2, kte-1
1627	DO i = i_start, itf
1628	ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25rdy/msfty(i,j)) &
1629	( muv(i,j+1)(v(i,k,j+1)+v(i,k-1,j+1))msfvy(i,j+1)* &
1630	(phb(i,k,j+1)-phb(i,k,j )+ph(i,k,j+1)-ph(i,k,j )) &
1631	+muv(i,j )(v(i,k,j )+v(i,k-1,j ))msfvy(i,j )* &
1632	(phb(i,k,j )-phb(i,k,j-1)+ph(i,k,j )-ph(i,k,j-1)) )
1633	ENDDO
1634	ENDDO
1635
1636	k = kte
1637	DO i = i_start, itf
1638	ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5rdy/msfty(i,j)) &
1639	( muv(i,j+1)(cfnv(i,k-1,j+1)+cfn1v(i,k-2,j+1))msfvy(i,j+1)* &
1640	(phb(i,k,j+1)-phb(i,k,j )+ph(i,k,j+1)-ph(i,k,j )) &
1641	+muv(i,j )(cfnv(i,k-1,j )+cfn1v(i,k-2,j ))msfvy(i,j )* &
1642	(phb(i,k,j )-phb(i,k,j-1)+ph(i,k,j )-ph(i,k,j-1)) )
1643	ENDDO
1644
1645	END IF
1646
1647	IF ( (config_flags%open_ye .or. specified) .and. jte >= jde-2 ) THEN
1648
1649	j = jde-2
1650	DO k = 2, kte-1
1651	DO i = i_start, itf
1652	ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25rdy/msfty(i,j)) &
1653	( muv(i,j+1)(v(i,k,j+1)+v(i,k-1,j+1))msfvy(i,j+1)* &
1654	(phb(i,k,j+1)-phb(i,k,j )+ph(i,k,j+1)-ph(i,k,j )) &
1655	+muv(i,j )(v(i,k,j )+v(i,k-1,j ))msfvy(i,j )* &
1656	(phb(i,k,j )-phb(i,k,j-1)+ph(i,k,j )-ph(i,k,j-1)) )
1657	ENDDO
1658	ENDDO
1659
1660	k = kte
1661	DO i = i_start, itf
1662	ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5rdy/msfty(i,j)) &
1663	( muv(i,j+1)(cfnv(i,k-1,j+1)+cfn1v(i,k-2,j+1))msfvy(i,j+1)* &
1664	(phb(i,k,j+1)-phb(i,k,j )+ph(i,k,j+1)-ph(i,k,j )) &
1665	+muv(i,j )(cfnv(i,k-1,j )+cfn1v(i,k-2,j ))msfvy(i,j )* &
1666	(phb(i,k,j )-phb(i,k,j-1)+ph(i,k,j )-ph(i,k,j-1)) )
1667	ENDDO
1668
1669	END IF
1670
1671	! x (u) advection
1672
1673	i_start = its
1674	j_start = jts
1675	itf=MIN(ite,ide-1)
1676	jtf=MIN(jte,jde-1)
1677
1678	IF ( (config_flags%open_xs) .and. its == ids ) i_start = its+2
1679	IF ( (config_flags%open_xe) .and. ite == ide ) itf = itf-3
1680
1681	DO j = j_start, jtf
1682
1683	DO k = 2, kte-1
1684	DO i = i_start, itf
1685	ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25rdx/msfty(i,j))( &
1686	( muu(i+1,j)(u(i+1,k,j)+u(i+1,k-1,j))msfux(i+1,j) &
1687	+muu(i,j )(u(i ,k,j)+u(i ,k-1,j))msfux(i ,j) )* (1./12.)*( &
1688	8.*(ph(i+1,k,j)-ph(i-1,k,j)) &
1689	-(ph(i+2,k,j)-ph(i-2,k,j)) &
1690	+8.*(phb(i+1,k,j)-phb(i-1,k,j)) &
1691	-(phb(i+2,k,j)-phb(i-2,k,j)) ) )
1692	ENDDO
1693	ENDDO
1694
1695	k = kte
1696	DO i = i_start, itf
1697	ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5rdx/msfty(i,j))( &
1698	( muu(i+1,j)(cfnu(i+1,k-1,j)+cfn1u(i+1,k-2,j))msfux(i+1,j) &
1699	+muu(i,j )(cfnu(i ,k-1,j)+cfn1u(i ,k-2,j))msfux(i ,j) )* (1./12.)*( &
1700	8.*(ph(i+1,k,j)-ph(i-1,k,j)) &
1701	-(ph(i+2,k,j)-ph(i-2,k,j)) &
1702	+8.*(phb(i+1,k,j)-phb(i-1,k,j)) &
1703	-(phb(i+2,k,j)-phb(i-2,k,j)) ) )
1704	ENDDO
1705
1706	ENDDO
1707
1708	! pick up near boundary rows using 2nd order stencil for open and specified conditions
1709
1710	IF ( (config_flags%open_xs .or. specified) .and. its <= ids+1 ) THEN
1711
1712	i = ids + 1
1713
1714	DO j = j_start, jtf
1715	DO k = 2, kte-1
1716	ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25rdx/msfty(i,j)) &
1717	( muu(i+1,j)(u(i+1,k,j)+u(i+1,k-1,j))msfux(i+1,j)* &
1718	(phb(i+1,k,j)-phb(i ,k,j)+ph(i+1,k,j)-ph(i ,k,j)) &
1719	+muu(i ,j)(u(i ,k,j)+u(i ,k-1,j))msfux(i ,j)* &
1720	(phb(i ,k,j)-phb(i-1,k,j)+ph(i ,k,j)-ph(i-1,k,j)) )
1721	ENDDO
1722	ENDDO
1723
1724	k = kte
1725	DO j = j_start, jtf
1726	ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5rdx/msfty(i,j)) &
1727	( muu(i+1,j)(cfnu(i+1,k-1,j)+cfn1u(i+1,k-2,j))msfux(i+1,j)* &
1728	(phb(i+1,k,j)-phb(i ,k,j)+ph(i+1,k,j)-ph(i ,k,j)) &
1729	+muu(i ,j)(cfnu(i ,k-1,j)+cfn1u(i ,k-2,j))msfux( i,j)* &
1730	(phb(i ,k,j)-phb(i-1,k,j)+ph(i ,k,j)-ph(i-1,k,j)) )
1731	ENDDO
1732
1733	END IF
1734
1735	IF ( (config_flags%open_xe .or. specified) .and. ite >= ide-2 ) THEN
1736
1737	i = ide-2
1738	DO j = j_start, jtf
1739	DO k = 2, kte-1
1740	ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25rdx/msfty(i,j)) &
1741	( muu(i+1,j)(u(i+1,k,j)+u(i+1,k-1,j))msfux(i+1,j)* &
1742	(phb(i+1,k,j)-phb(i ,k,j)+ph(i+1,k,j)-ph(i ,k,j)) &
1743	+muu(i ,j)(u(i ,k,j)+u(i ,k-1,j))msfux(i ,j)* &
1744	(phb(i ,k,j)-phb(i-1,k,j)+ph(i ,k,j)-ph(i-1,k,j)) )
1745	ENDDO
1746	ENDDO
1747
1748	k = kte
1749	DO j = j_start, jtf
1750	ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5rdx/msfty(i,j)) &
1751	( muu(i+1,j)(cfnu(i+1,k-1,j)+cfn1u(i+1,k-2,j))msfux(i+1,j)* &
1752	(phb(i+1,k,j)-phb(i ,k,j)+ph(i+1,k,j)-ph(i ,k,j)) &
1753	+muu(i ,j)(cfnu(i ,k-1,j)+cfn1u(i ,k-2,j))msfux( i,j)* &
1754	(phb(i ,k,j)-phb(i-1,k,j)+ph(i ,k,j)-ph(i-1,k,j)) )
1755	ENDDO
1756
1757	END IF
1758
1759	!--------------------------
1760
1761	ELSE IF (advective_order <= 6) THEN
1762
1763	!!! NOTE: this branch has been looked at and fixed with changes for overdecomposition
1764	!!! the branches covering the other advective_order cases have not. 20090923. JM
1765
1766	! y (v) advection
1767
1768	i_start = its
1769	j_start = jts
1770	itf=MIN(ite,ide-1)
1771	jtf=MIN(jte,jde-1)
1772
1773	IF (config_flags%open_ys .or. specified ) j_start = max(jts,jds+3)
1774	IF (config_flags%open_ye .or. specified ) jtf = min(jtf,jde-4)
1775
1776	DO j = j_start, jtf
1777
1778	DO k = 2, kte-1
1779	DO i = i_start, itf
1780	ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25rdy/msfty(i,j)) ( &
1781	( muv(i,j+1)(v(i,k,j+1)+v(i,k-1,j+1))msfvy(i,j+1) &
1782	+muv(i,j )(v(i,k,j )+v(i,k-1,j ))msfvy(i,j ) )* (1./60.)*( &
1783	45.*(ph(i,k,j+1)-ph(i,k,j-1)) &
1784	-9.*(ph(i,k,j+2)-ph(i,k,j-2)) &
1785	+(ph(i,k,j+3)-ph(i,k,j-3)) &
1786	+45.*(phb(i,k,j+1)-phb(i,k,j-1)) &
1787	-9.*(phb(i,k,j+2)-phb(i,k,j-2)) &
1788	+(phb(i,k,j+3)-phb(i,k,j-3)) ) )
1789
1790
1791	ENDDO
1792	ENDDO
1793
1794	k = kte
1795	DO i = i_start, itf
1796	ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5rdy/msfty(i,j)) ( &
1797	( muv(i,j+1)(cfnv(i,k-1,j+1)+cfn1v(i,k-2,j+1))msfvy(i,j+1) &
1798	+muv(i,j )(cfnv(i,k-1,j )+cfn1v(i,k-2,j ))msfvy(i,j ) )* (1./60.)*( &
1799	45.*(ph(i,k,j+1)-ph(i,k,j-1)) &
1800	-9.*(ph(i,k,j+2)-ph(i,k,j-2)) &
1801	+(ph(i,k,j+3)-ph(i,k,j-3)) &
1802	+45.*(phb(i,k,j+1)-phb(i,k,j-1)) &
1803	-9.*(phb(i,k,j+2)-phb(i,k,j-2)) &
1804	+(phb(i,k,j+3)-phb(i,k,j-3)) ) )
1805
1806	ENDDO
1807
1808	ENDDO
1809
1810	! 4th order stencil for open or specified conditions two in form the boundary
1811
1812	IF ( (config_flags%open_ys .or. specified) .and. jts <= jds+2 .and. jds+2 <= jte ) THEN
1813
1814	j = jds+2
1815	DO k = 2, kte-1
1816	DO i = i_start, itf
1817	ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25rdy/msfty(i,j)) ( &
1818	( muv(i,j+1)(v(i,k,j+1)+v(i,k-1,j+1))msfvy(i,j+1) &
1819	+muv(i,j )(v(i,k,j )+v(i,k-1,j ))msfvy(i,j ) )* (1./12.)*( &
1820	8.*(ph(i,k,j+1)-ph(i,k,j-1)) &
1821	-(ph(i,k,j+2)-ph(i,k,j-2)) &
1822	+8.*(phb(i,k,j+1)-phb(i,k,j-1)) &
1823	-(phb(i,k,j+2)-phb(i,k,j-2)) ) )
1824
1825
1826	ENDDO
1827	ENDDO
1828
1829	k = kte
1830	DO i = i_start, itf
1831	ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5rdy/msfty(i,j)) ( &
1832	( muv(i,j+1)(cfnv(i,k-1,j+1)+cfn1v(i,k-2,j+1))msfvy(i,j+1) &
1833	+muv(i,j )(cfnv(i,k-1,j )+cfn1v(i,k-2,j ))msfvy(i,j) )* (1./12.)*( &
1834	8.*(ph(i,k,j+1)-ph(i,k,j-1)) &
1835	-(ph(i,k,j+2)-ph(i,k,j-2)) &
1836	+8.*(phb(i,k,j+1)-phb(i,k,j-1)) &
1837	-(phb(i,k,j+2)-phb(i,k,j-2)) ) )
1838
1839	ENDDO
1840
1841	END IF
1842
1843	IF ( (config_flags%open_ye .or. specified) .and. jts <= jde-3 .and. jde-3 <= jte ) THEN
1844	j = jde-3
1845	DO k = 2, kte-1
1846	DO i = i_start, itf
1847	ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25rdy/msfty(i,j)) ( &
1848	( muv(i,j+1)(v(i,k,j+1)+v(i,k-1,j+1))msfvy(i,j+1) &
1849	+muv(i,j )(v(i,k,j )+v(i,k-1,j ))msfvy(i,j) )* (1./12.)*( &
1850	8.*(ph(i,k,j+1)-ph(i,k,j-1)) &
1851	-(ph(i,k,j+2)-ph(i,k,j-2)) &
1852	+8.*(phb(i,k,j+1)-phb(i,k,j-1)) &
1853	-(phb(i,k,j+2)-phb(i,k,j-2)) ) )
1854
1855
1856	ENDDO
1857	ENDDO
1858
1859	k = kte
1860	DO i = i_start, itf
1861	ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5rdy/msfty(i,j)) ( &
1862	( muv(i,j+1)(cfnv(i,k-1,j+1)+cfn1v(i,k-2,j+1))msfvy(i,j+1) &
1863	+muv(i,j )(cfnv(i,k-1,j )+cfn1v(i,k-2,j ))msfvy(i,j) )* (1./12.)*( &
1864	8.*(ph(i,k,j+1)-ph(i,k,j-1)) &
1865	-(ph(i,k,j+2)-ph(i,k,j-2)) &
1866	+8.*(phb(i,k,j+1)-phb(i,k,j-1)) &
1867	-(phb(i,k,j+2)-phb(i,k,j-2)) ) )
1868
1869	ENDDO
1870
1871	END IF
1872
1873	! 2nd order stencil for open and specified conditions one row in from boundary
1874
1875	IF ( (config_flags%open_ys .or. specified) .and. jts <= jds+1 .and. jds+1 <= jte ) THEN
1876
1877	j = jds+1
1878	DO k = 2, kte-1
1879	DO i = i_start, itf
1880	ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25rdy/msfty(i,j)) &
1881	( muv(i,j+1)(v(i,k,j+1)+v(i,k-1,j+1))msfvy(i,j+1)* &
1882	(phb(i,k,j+1)-phb(i,k,j )+ph(i,k,j+1)-ph(i,k,j )) &
1883	+muv(i,j )(v(i,k,j )+v(i,k-1,j ))msfvy(i,j )* &
1884	(phb(i,k,j )-phb(i,k,j-1)+ph(i,k,j )-ph(i,k,j-1)) )
1885	ENDDO
1886	ENDDO
1887
1888	k = kte
1889	DO i = i_start, itf
1890	ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5rdy/msfty(i,j)) &
1891	( muv(i,j+1)(cfnv(i,k-1,j+1)+cfn1v(i,k-2,j+1))msfvy(i,j+1)* &
1892	(phb(i,k,j+1)-phb(i,k,j )+ph(i,k,j+1)-ph(i,k,j )) &
1893	+muv(i,j )(cfnv(i,k-1,j )+cfn1v(i,k-2,j ))msfvy(i,j )* &
1894	(phb(i,k,j )-phb(i,k,j-1)+ph(i,k,j )-ph(i,k,j-1)) )
1895	ENDDO
1896
1897	END IF
1898
1899	IF ( (config_flags%open_ye .or. specified) .and. jts <= jde-2 .and. jde-2 <= jte ) THEN
1900
1901	j = jde-2
1902	DO k = 2, kte-1
1903	DO i = i_start, itf
1904	ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25rdy/msfty(i,j)) &
1905	( muv(i,j+1)(v(i,k,j+1)+v(i,k-1,j+1))msfvy(i,j+1)* &
1906	(phb(i,k,j+1)-phb(i,k,j )+ph(i,k,j+1)-ph(i,k,j )) &
1907	+muv(i,j )(v(i,k,j )+v(i,k-1,j ))msfvy(i,j )* &
1908	(phb(i,k,j )-phb(i,k,j-1)+ph(i,k,j )-ph(i,k,j-1)) )
1909	ENDDO
1910	ENDDO
1911
1912	k = kte
1913	DO i = i_start, itf
1914	ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5rdy/msfty(i,j)) &
1915	( muv(i,j+1)(cfnv(i,k-1,j+1)+cfn1v(i,k-2,j+1))msfvy(i,j+1)* &
1916	(phb(i,k,j+1)-phb(i,k,j )+ph(i,k,j+1)-ph(i,k,j )) &
1917	+muv(i,j )(cfnv(i,k-1,j )+cfn1v(i,k-2,j ))msfvy(i,j )* &
1918	(phb(i,k,j )-phb(i,k,j-1)+ph(i,k,j )-ph(i,k,j-1)) )
1919	ENDDO
1920
1921	END IF
1922
1923	! x (u) advection
1924
1925	i_start = its
1926	j_start = jts
1927	itf=MIN(ite,ide-1)
1928	jtf=MIN(jte,jde-1)
1929
1930	IF (config_flags%open_xs .or. specified ) i_start = max(its,ids+3)
1931	IF (config_flags%open_xe .or. specified ) itf = min(itf,ide-4)
1932
1933	DO j = j_start, jtf
1934
1935	DO k = 2, kte-1
1936	DO i = i_start, itf
1937	ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25rdx/msfty(i,j))( &
1938	( muu(i+1,j)(u(i+1,k,j)+u(i+1,k-1,j))msfux(i+1,j) &
1939	+muu(i,j )(u(i,k,j )+u(i,k-1,j ))msfux(i,j) )* (1./60.)*( &
1940	45.*(ph(i+1,k,j)-ph(i-1,k,j)) &
1941	-9.*(ph(i+2,k,j)-ph(i-2,k,j)) &
1942	+(ph(i+3,k,j)-ph(i-3,k,j)) &
1943	+45.*(phb(i+1,k,j)-phb(i-1,k,j)) &
1944	-9.*(phb(i+2,k,j)-phb(i-2,k,j)) &
1945	+(phb(i+3,k,j)-phb(i-3,k,j)) ) )
1946	ENDDO
1947	ENDDO
1948
1949	k = kte
1950	DO i = i_start, itf
1951	ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5rdx/msfty(i,j))( &
1952	( muu(i+1,j)(cfnu(i+1,k-1,j)+cfn1u(i+1,k-2,j))msfux(i+1,j) &
1953	+muu(i,j )(cfnu(i ,k-1,j)+cfn1u(i,k-2,j))msfux(i,j) )* (1./60.)*( &
1954	45.*(ph(i+1,k,j)-ph(i-1,k,j)) &
1955	-9.*(ph(i+2,k,j)-ph(i-2,k,j)) &
1956	+(ph(i+3,k,j)-ph(i-3,k,j)) &
1957	+45.*(phb(i+1,k,j)-phb(i-1,k,j)) &
1958	-9.*(phb(i+2,k,j)-phb(i-2,k,j)) &
1959	+(phb(i+3,k,j)-phb(i-3,k,j)) ) )
1960	ENDDO
1961
1962	ENDDO
1963
1964	! 4th order stencil two in from the boundary for open and specified conditions
1965
1966	IF ( (config_flags%open_xs) .and. its <= ids+2 .and. ids+2 <= ite ) THEN
1967	i = ids + 2
1968	DO j = j_start, jtf
1969	DO k = 2, kte-1
1970	ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25rdx/msfty(i,j))( &
1971	( muu(i+1,j)(u(i+1,k,j)+u(i+1,k-1,j))msfux(i+1,j) &
1972	+muu(i,j )(u(i,k,j )+u(i,k-1,j ))msfux(i,j) )* (1./12.)*( &
1973	8.*(ph(i+1,k,j)-ph(i-1,k,j)) &
1974	-(ph(i+2,k,j)-ph(i-2,k,j)) &
1975	+8.*(phb(i+1,k,j)-phb(i-1,k,j)) &
1976	-(phb(i+2,k,j)-phb(i-2,k,j)) ) )
1977	ENDDO
1978	k = kte
1979	ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5rdx/msfty(i,j))( &
1980	( muu(i+1,j)(cfnu(i+1,k-1,j)+cfn1u(i+1,k-2,j))msfux(i+1,j) &
1981	+muu(i,j )(cfnu(i ,k-1,j)+cfn1u(i,k-2,j))msfux(i,j) )* (1./12.)*( &
1982	8.*(ph(i+1,k,j)-ph(i-1,k,j)) &
1983	-(ph(i+2,k,j)-ph(i-2,k,j)) &
1984	+8.*(phb(i+1,k,j)-phb(i-1,k,j)) &
1985	-(phb(i+2,k,j)-phb(i-2,k,j)) ) )
1986
1987	ENDDO
1988	END IF
1989
1990	IF ( (config_flags%open_xe) .and. its <= ide-3 .and. ide-3 <= ite ) THEN
1991	i = ide-3
1992	DO j = j_start, jtf
1993	DO k = 2, kte-1
1994	ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25rdx/msfty(i,j))( &
1995	( muu(i+1,j)(u(i+1,k,j)+u(i+1,k-1,j))msfux(i+1,j) &
1996	+muu(i,j )(u(i,k,j )+u(i,k-1,j ))msfux(i,j) )* (1./12.)*( &
1997	8.*(ph(i+1,k,j)-ph(i-1,k,j)) &
1998	-(ph(i+2,k,j)-ph(i-2,k,j)) &
1999	+8.*(phb(i+1,k,j)-phb(i-1,k,j)) &
2000	-(phb(i+2,k,j)-phb(i-2,k,j)) ) )
2001	ENDDO
2002	k = kte
2003	ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5rdx/msfty(i,j))( &
2004	( muu(i+1,j)(cfnu(i+1,k-1,j)+cfn1u(i+1,k-2,j))msfux(i+1,j) &
2005	+muu(i,j )(cfnu(i ,k-1,j)+cfn1u(i,k-2,j))msfux(i,j) )* (1./12.)*( &
2006	8.*(ph(i+1,k,j)-ph(i-1,k,j)) &
2007	-(ph(i+2,k,j)-ph(i-2,k,j)) &
2008	+8.*(phb(i+1,k,j)-phb(i-1,k,j)) &
2009	-(phb(i+2,k,j)-phb(i-2,k,j)) ) )
2010
2011	ENDDO
2012	END IF
2013
2014	! 2nd order stencil for open and specified conditions one in from the boundary
2015
2016	IF ( (config_flags%open_xs .or. specified) .and. its <= ids+1 .and. ids+1 <= ite ) THEN
2017
2018	i = ids + 1
2019
2020	DO j = j_start, jtf
2021	DO k = 2, kte-1
2022	ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25rdx/msfty(i,j)) &
2023	( muu(i+1,j)(u(i+1,k,j)+u(i+1,k-1,j))msfux(i+1,j)* &
2024	(phb(i+1,k,j)-phb(i ,k,j)+ph(i+1,k,j)-ph(i ,k,j)) &
2025	+muu(i ,j)(u(i ,k,j)+u(i ,k-1,j))msfux(i ,j)* &
2026	(phb(i ,k,j)-phb(i-1,k,j)+ph(i ,k,j)-ph(i-1,k,j)) )
2027	ENDDO
2028	ENDDO
2029
2030	k = kte
2031	DO j = j_start, jtf
2032	ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5rdx/msfty(i,j)) &
2033	( muu(i+1,j)(cfnu(i+1,k-1,j)+cfn1u(i+1,k-2,j))msfux(i+1,j)* &
2034	(phb(i+1,k,j)-phb(i ,k,j)+ph(i+1,k,j)-ph(i ,k,j)) &
2035	+muu(i ,j)(cfnu(i ,k-1,j)+cfn1u(i ,k-2,j))msfux( i,j)* &
2036	(phb(i ,k,j)-phb(i-1,k,j)+ph(i ,k,j)-ph(i-1,k,j)) )
2037	ENDDO
2038
2039	END IF
2040
2041	IF ( (config_flags%open_xe .or. specified) .and. its <= ide-2 .and. ide-2 <= ite ) THEN
2042
2043	i = ide-2
2044	DO j = j_start, jtf
2045	DO k = 2, kte-1
2046	ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25rdx/msfty(i,j)) &
2047	( muu(i+1,j)(u(i+1,k,j)+u(i+1,k-1,j))msfux(i+1,j)* &
2048	(phb(i+1,k,j)-phb(i ,k,j)+ph(i+1,k,j)-ph(i ,k,j)) &
2049	+muu(i ,j)(u(i ,k,j)+u(i ,k-1,j))msfux(i ,j)* &
2050	(phb(i ,k,j)-phb(i-1,k,j)+ph(i ,k,j)-ph(i-1,k,j)) )
2051	ENDDO
2052	ENDDO
2053
2054	k = kte
2055	DO j = j_start, jtf
2056	ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5rdx/msfty(i,j)) &
2057	( muu(i+1,j)(cfnu(i+1,k-1,j)+cfn1u(i+1,k-2,j))msfux(i+1,j)* &
2058	(phb(i+1,k,j)-phb(i ,k,j)+ph(i+1,k,j)-ph(i ,k,j)) &
2059	+muu(i ,j)(cfnu(i ,k-1,j)+cfn1u(i ,k-2,j))msfux( i,j)* &
2060	(phb(i ,k,j)-phb(i-1,k,j)+ph(i ,k,j)-ph(i-1,k,j)) )
2061	ENDDO
2062
2063	END IF
2064
2065	END IF ! 6th order advection
2066
2067	! lateral open boundary conditions,
2068	! start with north and south (y) boundaries
2069
2070	i_start = its
2071	itf=MIN(ite,ide-1)
2072
2073	! south
2074
2075	IF ( (config_flags%open_ys) .and. jts == jds ) THEN
2076
2077	j=jts
2078
2079	DO k=2,kde
2080	kz = min(k,kde-1)
2081	DO i = its,itf
2082	vb =.5( fnm(kz)(v(i,kz ,j+1)+v(i,kz ,j )) &
2083	+fnp(kz)*(v(i,kz-1,j+1)+v(i,kz-1,j )) )
2084	vl=amin1(vb,0.)
2085	ph_tend(i,k,j)=ph_tend(i,k,j)-rdymut(i,j)( &
2086	+vl*(ph_old(i,k,j+1)-ph_old(i,k,j)))
2087	ENDDO
2088	ENDDO
2089
2090	END IF
2091
2092	! north
2093
2094	IF ( (config_flags%open_ye) .and. jte == jde ) THEN
2095
2096	j=jte-1
2097
2098	DO k=2,kde
2099	kz = min(k,kde-1)
2100	DO i = its,itf
2101	vb=.5( fnm(kz)(v(i,kz ,j+1)+v(i,kz ,j)) &
2102	+fnp(kz)*(v(i,kz-1,j+1)+v(i,kz-1,j)) )
2103	vr=amax1(vb,0.)
2104	ph_tend(i,k,j)=ph_tend(i,k,j)-rdymut(i,j)( &
2105	+vr*(ph_old(i,k,j)-ph_old(i,k,j-1)))
2106	ENDDO
2107	ENDDO
2108
2109	END IF
2110
2111	! now the east and west (y) boundaries
2112
2113	j_start = its
2114	jtf=MIN(jte,jde-1)
2115
2116	! west
2117
2118	IF ( (config_flags%open_xs) .and. its == ids ) THEN
2119
2120	i=its
2121
2122	DO j = jts,jtf
2123	DO k=2,kde-1
2124	kz = k
2125	ub =.5( fnm(kz)(u(i+1,kz ,j)+u(i ,kz ,j)) &
2126	+fnp(kz)*(u(i+1,kz-1,j)+u(i ,kz-1,j)) )
2127	ul=amin1(ub,0.)
2128	ph_tend(i,k,j)=ph_tend(i,k,j)-(msftx(i,j)/msfty(i,j))rdxmut(i,j)*( &
2129	+ul*(ph_old(i+1,k,j)-ph_old(i,k,j)))
2130	ENDDO
2131
2132	k = kde
2133	kz = k
2134	ub =.5( fnm(kz)(u(i+1,kz ,j)+u(i ,kz ,j)) &
2135	+fnp(kz)*(u(i+1,kz-1,j)+u(i ,kz-1,j)) )
2136	ul=amin1(ub,0.)
2137	ph_tend(i,k,j)=ph_tend(i,k,j)-(msftx(i,j)/msfty(i,j))rdxmut(i,j)*( &
2138	+ul*(ph_old(i+1,k,j)-ph_old(i,k,j)))
2139	ENDDO
2140
2141	END IF
2142
2143	! east
2144
2145	IF ( (config_flags%open_xe) .and. ite == ide ) THEN
2146
2147	i = ite-1
2148
2149	DO j = jts,jtf
2150	DO k=2,kde-1
2151	kz = k
2152	ub=.5( fnm(kz)(u(i+1,kz ,j)+u(i,kz ,j)) &
2153	+fnp(kz)*(u(i+1,kz-1,j)+u(i,kz-1,j)) )
2154	ur=amax1(ub,0.)
2155	ph_tend(i,k,j)=ph_tend(i,k,j)-(msftx(i,j)/msfty(i,j))rdxmut(i,j)*( &
2156	+ur*(ph_old(i,k,j)-ph_old(i-1,k,j)))
2157	ENDDO
2158
2159	k = kde
2160	kz = k-1
2161	ub=.5( fnm(kz)(u(i+1,kz ,j)+u(i,kz ,j)) &
2162	+fnp(kz)*(u(i+1,kz-1,j)+u(i,kz-1,j)) )
2163	ur=amax1(ub,0.)
2164	ph_tend(i,k,j)=ph_tend(i,k,j)-(msftx(i,j)/msfty(i,j))rdxmut(i,j)*( &
2165	+ur*(ph_old(i,k,j)-ph_old(i-1,k,j)))
2166
2167	ENDDO
2168
2169	END IF
2170
2171	END SUBROUTINE rhs_ph
2172
2173
2174	!-------------------------------------------------------------------------------
2175
2176	SUBROUTINE horizontal_pressure_gradient( ru_tend,rv_tend, &
2177	ph,alt,p,pb,al,php,cqu,cqv, &
2178	muu,muv,mu,fnm,fnp,rdnw, &
2179	cf1,cf2,cf3,rdx,rdy,msfux,msfuy,&
2180	msfvx,msfvy,msftx,msfty, &
2181	config_flags, non_hydrostatic, &
2182	top_lid, &
2183	ids, ide, jds, jde, kds, kde, &
2184	ims, ime, jms, jme, kms, kme, &
2185	its, ite, jts, jte, kts, kte )
2186
2187	IMPLICIT NONE
2188
2189	! Input data
2190
2191
2192	TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
2193
2194	LOGICAL, INTENT (IN ) :: non_hydrostatic, top_lid
2195
2196	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
2197	ims, ime, jms, jme, kms, kme, &
2198	its, ite, jts, jte, kts, kte
2199
2200	REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(IN ) :: &
2201	ph, &
2202	alt, &
2203	al, &
2204	p, &
2205	pb, &
2206	php, &
2207	cqu, &
2208	cqv
2209
2210
2211	REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(INOUT) :: &
2212	ru_tend, &
2213	rv_tend
2214
2215	REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN ) :: muu, muv, mu, &
2216	msfux, msfuy, &
2217	msfvx, msfvy, &
2218	msftx, msfty
2219
2220	REAL, DIMENSION( kms:kme ), INTENT(IN ) :: rdnw, fnm, fnp
2221
2222	REAL, INTENT(IN ) :: rdx, rdy, cf1, cf2, cf3
2223
2224	INTEGER :: i,j,k, itf, jtf, ktf, i_start, j_start
2225	REAL, DIMENSION( ims:ime, kms:kme ) :: dpn
2226	REAL :: dpx, dpy
2227
2228	LOGICAL :: specified
2229
2230	!<DESCRIPTION>
2231	!
2232	! horizontal_pressure_gradient calculates the
2233	! horizontal pressure gradient terms for the large-timestep tendency
2234	! in the horizontal momentum equations (u,v).
2235	!
2236	!</DESCRIPTION>
2237
2238	specified = .false.
2239	if(config_flags%specified .or. config_flags%nested) specified = .true.
2240
2241	! Notes on map scale factors:
2242	! Calculates the pressure gradient terms in ADT eqns 44 and 45
2243	! With upper rho -> 'mu', these are:
2244	! Eqn 30: -mu(mx/my)(1/rho)*partial dp/dx
2245	! Eqn 31: -mu(my/mx)(1/rho)*partial dp/dy
2246	!
2247	! As we are on nu, rather than height, surfaces:
2248	!
2249	! mu dp/dx = mu alpha partial dp'/dx + (nu mu partial dmubar/dx) alpha'
2250	! + mu partial dphi'/dx + (partial dphi/dx)*(partial dp'/dnu - mu')
2251	!
2252	! mu dp/dy = mu alpha partial dp'/dy + (nu mu partial dmubar/dy) alpha'
2253	! + mu partial dphi'/dy + (partial dphi/dy)*(partial dp'/dnu - mu')
2254
2255	! start with the north-south (y) pressure gradient
2256
2257	itf=MIN(ite,ide-1)
2258	jtf=jte
2259	ktf=MIN(kte,kde-1)
2260	i_start = its
2261	j_start = jts
2262	IF ( (config_flags%open_ys .or. specified .or. &
2263	config_flags%nested .or. config_flags%polar ) .and. jts == jds ) j_start = jts+1
2264	IF ( (config_flags%open_ye .or. specified .or. &
2265	config_flags%nested .or. config_flags%polar ) .and. jte == jde ) jtf = jtf-1
2266
2267	DO j = j_start, jtf
2268
2269	IF ( non_hydrostatic ) THEN
2270
2271	k=1
2272
2273	DO i = i_start, itf
2274	dpn(i,k) = .5( cf1(p(i,k ,j-1)+p(i,k ,j)) &
2275	+cf2*(p(i,k+1,j-1)+p(i,k+1,j)) &
2276	+cf3*(p(i,k+2,j-1)+p(i,k+2,j)) )
2277	dpn(i,kde) = 0.
2278	ENDDO
2279	IF (top_lid) THEN
2280	DO i = i_start, itf
2281	dpn(i,kde) = .5( cf1(p(i,kde-1,j-1)+p(i,kde-1,j)) &
2282	+cf2*(p(i,kde-2,j-1)+p(i,kde-2,j)) &
2283	+cf3*(p(i,kde-3,j-1)+p(i,kde-3,j)) )
2284	ENDDO
2285	ENDIF
2286
2287	DO k=2,ktf
2288	DO i = i_start, itf
2289	dpn(i,k) = .5( fnm(k)(p(i,k ,j-1)+p(i,k ,j)) &
2290	+fnp(k)*(p(i,k-1,j-1)+p(i,k-1,j)) )
2291	END DO
2292	END DO
2293
2294	! ADT eqn 45: -mu(my/mx)(1/rho)*partial dp/dy
2295	! [alt, al are 1/rho terms; muv, mu are NOT coupled]
2296	DO K=1,ktf
2297	DO i = i_start, itf
2298	! Here are mu dp/dy terms 1-3
2299	dpy = (msfvy(i,j)/msfvx(i,j)).5rdymuv(i,j)( &
2300	(ph (i,k+1,j)-ph (i,k+1,j-1) + ph(i,k,j)-ph(i,k,j-1)) &
2301	+(alt(i,k ,j)+alt(i,k ,j-1))*(p (i,k,j)-p (i,k,j-1)) &
2302	+(al (i,k ,j)+al (i,k ,j-1))*(pb(i,k,j)-pb(i,k,j-1)) )
2303	! Here is mu dp/dy term 4
2304	dpy = dpy + (msfvy(i,j)/msfvx(i,j))rdy(php(i,k,j)-php(i,k,j-1))* &
2305	(rdnw(k)(dpn(i,k+1)-dpn(i,k))-.5(mu(i,j-1)+mu(i,j)))
2306	rv_tend(i,k,j) = rv_tend(i,k,j)-cqv(i,k,j)*dpy
2307	END DO
2308	END DO
2309
2310	ELSE
2311
2312	! ADT eqn 45: -mu(my/mx)(1/rho)*partial dp/dy
2313	! [alt, al are 1/rho terms; muv, mu are NOT coupled]
2314	DO K=1,ktf
2315	DO i = i_start, itf
2316	! Here are mu dp/dy terms 1-3; term 4 not needed if hydrostatic
2317	dpy = (msfvy(i,j)/msfvx(i,j)).5rdymuv(i,j)( &
2318	(ph (i,k+1,j)-ph (i,k+1,j-1) + ph(i,k,j)-ph(i,k,j-1)) &
2319	+(alt(i,k ,j)+alt(i,k ,j-1))*(p (i,k,j)-p (i,k,j-1)) &
2320	+(al (i,k ,j)+al (i,k ,j-1))*(pb(i,k,j)-pb(i,k,j-1)) )
2321	rv_tend(i,k,j) = rv_tend(i,k,j)-cqv(i,k,j)*dpy
2322	END DO
2323	END DO
2324
2325	END IF
2326
2327	ENDDO
2328
2329	! now the east-west (x) pressure gradient
2330
2331	itf=ite
2332	jtf=MIN(jte,jde-1)
2333	ktf=MIN(kte,kde-1)
2334	i_start = its
2335	j_start = jts
2336	IF ( (config_flags%open_xs .or. specified .or. &
2337	config_flags%nested ) .and. its == ids ) i_start = its+1
2338	IF ( (config_flags%open_xe .or. specified .or. &
2339	config_flags%nested ) .and. ite == ide ) itf = itf-1
2340	IF ( config_flags%periodic_x ) i_start = its
2341	IF ( config_flags%periodic_x ) itf=ite
2342
2343	DO j = j_start, jtf
2344
2345	IF ( non_hydrostatic ) THEN
2346
2347	k=1
2348
2349	DO i = i_start, itf
2350	dpn(i,k) = .5( cf1(p(i-1,k ,j)+p(i,k ,j)) &
2351	+cf2*(p(i-1,k+1,j)+p(i,k+1,j)) &
2352	+cf3*(p(i-1,k+2,j)+p(i,k+2,j)) )
2353	dpn(i,kde) = 0.
2354	ENDDO
2355	IF (top_lid) THEN
2356	DO i = i_start, itf
2357	dpn(i,kde) = .5( cf1(p(i-1,kde-1,j)+p(i,kde-1,j)) &
2358	+cf2*(p(i-1,kde-2,j)+p(i,kde-2,j)) &
2359	+cf3*(p(i-1,kde-3,j)+p(i,kde-3,j)) )
2360	ENDDO
2361	ENDIF
2362
2363	DO k=2,ktf
2364	DO i = i_start, itf
2365	dpn(i,k) = .5( fnm(k)(p(i-1,k ,j)+p(i,k ,j)) &
2366	+fnp(k)*(p(i-1,k-1,j)+p(i,k-1,j)) )
2367	END DO
2368	END DO
2369
2370	! ADT eqn 44: -mu(mx/my)(1/rho)*partial dp/dx
2371	! [alt, al are 1/rho terms; muu, mu are NOT coupled]
2372	DO K=1,ktf
2373	DO i = i_start, itf
2374	! Here are mu dp/dy terms 1-3
2375	dpx = (msfux(i,j)/msfuy(i,j)).5rdxmuu(i,j)( &
2376	(ph (i,k+1,j)-ph (i-1,k+1,j) + ph(i,k,j)-ph(i-1,k,j)) &
2377	+(alt(i,k ,j)+alt(i-1,k ,j))*(p (i,k,j)-p (i-1,k,j)) &
2378	+(al (i,k ,j)+al (i-1,k ,j))*(pb(i,k,j)-pb(i-1,k,j)) )
2379	! Here is mu dp/dy term 4
2380	dpx = dpx + (msfux(i,j)/msfuy(i,j))rdx(php(i,k,j)-php(i-1,k,j))* &
2381	(rdnw(k)(dpn(i,k+1)-dpn(i,k))-.5(mu(i-1,j)+mu(i,j)))
2382	ru_tend(i,k,j) = ru_tend(i,k,j)-cqu(i,k,j)*dpx
2383	END DO
2384	END DO
2385
2386	ELSE
2387
2388	! ADT eqn 44: -mu(mx/my)(1/rho)*partial dp/dx
2389	! [alt, al are 1/rho terms; muu, mu are NOT coupled]
2390	DO K=1,ktf
2391	DO i = i_start, itf
2392	! Here are mu dp/dy terms 1-3; term 4 not needed if hydrostatic
2393	dpx = (msfux(i,j)/msfuy(i,j)).5rdxmuu(i,j)( &
2394	(ph (i,k+1,j)-ph (i-1,k+1,j) + ph(i,k,j)-ph(i-1,k,j)) &
2395	+(alt(i,k ,j)+alt(i-1,k ,j))*(p (i,k,j)-p (i-1,k,j)) &
2396	+(al (i,k ,j)+al (i-1,k ,j))*(pb(i,k,j)-pb(i-1,k,j)) )
2397	ru_tend(i,k,j) = ru_tend(i,k,j)-cqu(i,k,j)*dpx
2398	END DO
2399	END DO
2400
2401	END IF
2402
2403	ENDDO
2404
2405	END SUBROUTINE horizontal_pressure_gradient
2406
2407	!-------------------------------------------------------------------------------
2408
2409	SUBROUTINE pg_buoy_w( rw_tend, p, cqw, mu, mub, &
2410	rdnw, rdn, g, msftx, msfty, &
2411	ids, ide, jds, jde, kds, kde, &
2412	ims, ime, jms, jme, kms, kme, &
2413	its, ite, jts, jte, kts, kte )
2414
2415	IMPLICIT NONE
2416
2417	! Input data
2418
2419	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
2420	ims, ime, jms, jme, kms, kme, &
2421	its, ite, jts, jte, kts, kte
2422
2423	REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(IN ) :: p
2424	REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(INOUT) :: cqw
2425
2426
2427	REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(INOUT) :: rw_tend
2428
2429	REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN ) :: mub, mu, msftx, msfty
2430
2431	REAL, DIMENSION( kms:kme ), INTENT(IN ) :: rdnw, rdn
2432
2433	REAL, INTENT(IN ) :: g
2434
2435	INTEGER :: itf, jtf, i, j, k
2436	REAL :: cq1, cq2
2437
2438
2439	!<DESCRIPTION>
2440	!
2441	! pg_buoy_w calculates the
2442	! vertical pressure gradient and buoyancy terms for the large-timestep
2443	! tendency in the vertical momentum equation.
2444	!
2445	!</DESCRIPTION>
2446
2447	! BUOYANCY AND PRESSURE GRADIENT TERM IN W EQUATION AT TIME T
2448
2449	! Map scale factor notes
2450	! ADT eqn 46 RHS terms 6 and 7 (where 7 is "-rho g")
2451	! Dividing by my, and using mu and nu (see Klemp et al. eqns 32, 40)
2452	! term 6: +(g/my) partial dp'/dnu
2453	! term 7: -(g/my) mu'
2454	!
2455	! For moisture-free atmosphere, cq1=1, cq2=0
2456	! => (1./msft(i,j)) * g * [rdn(k)*{p(i,k,j)-p(i,k-1,j)}-mu(i,j)]
2457
2458	itf=MIN(ite,ide-1)
2459	jtf=MIN(jte,jde-1)
2460
2461	DO j = jts,jtf
2462
2463	k=kde
2464	DO i=its,itf
2465	cq1 = 1./(1.+cqw(i,k-1,j))
2466	cq2 = cqw(i,k-1,j)*cq1
2467	rw_tend(i,k,j) = rw_tend(i,k,j)+(1./msfty(i,j))g( &
2468	cq12.rdnw(k-1)*( -p(i,k-1,j)) &
2469	-mu(i,j)-cq2*mub(i,j) )
2470	END DO
2471
2472	DO k = 2, kde-1
2473	DO i = its,itf
2474	cq1 = 1./(1.+cqw(i,k,j))
2475	cq2 = cqw(i,k,j)*cq1
2476	cqw(i,k,j) = cq1
2477	rw_tend(i,k,j) = rw_tend(i,k,j)+(1./msfty(i,j))g( &
2478	cq1rdn(k)(p(i,k,j)-p(i,k-1,j)) &
2479	-mu(i,j)-cq2*mub(i,j) )
2480	END DO
2481	ENDDO
2482
2483
2484	ENDDO
2485
2486	END SUBROUTINE pg_buoy_w
2487
2488	!-------------------------------------------------------------------------------
2489
2490	SUBROUTINE w_damp( rw_tend, max_vert_cfl,max_horiz_cfl, &
2491	u, v, ww, w, mut, rdnw, &
2492	rdx, rdy, msfux, msfuy, &
2493	msfvx, msfvy, dt, &
2494	config_flags, &
2495	ids, ide, jds, jde, kds, kde, &
2496	ims, ime, jms, jme, kms, kme, &
2497	its, ite, jts, jte, kts, kte )
2498
2499	USE module_llxy
2500	IMPLICIT NONE
2501
2502	! Input data
2503
2504	TYPE(grid_config_rec_type ) , INTENT(IN ) :: config_flags
2505
2506	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
2507	ims, ime, jms, jme, kms, kme, &
2508	its, ite, jts, jte, kts, kte
2509
2510	REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(IN ) :: u, v, ww, w
2511
2512	REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(INOUT) :: rw_tend
2513
2514	REAL, INTENT(OUT) :: max_vert_cfl
2515	REAL, INTENT(OUT) :: max_horiz_cfl
2516	REAL :: horiz_cfl
2517
2518	REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN ) :: mut
2519
2520	REAL, DIMENSION( kms:kme ), INTENT(IN ) :: rdnw
2521
2522	REAL, INTENT(IN) :: dt
2523	REAL, INTENT(IN) :: rdx, rdy
2524	REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: msfux, msfuy
2525	REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: msfvx, msfvy
2526
2527	REAL :: vert_cfl, cf_n, cf_d, maxdub, maxdeta
2528
2529	INTEGER :: itf, jtf, i, j, k, maxi, maxj, maxk
2530	INTEGER :: some
2531	CHARACTER*512 :: temp
2532
2533	CHARACTER (LEN=256) :: time_str
2534	CHARACTER (LEN=256) :: grid_str
2535
2536	integer :: total
2537	REAL :: msfuxt , msfxffl
2538
2539	!<DESCRIPTION>
2540	!
2541	! w_damp computes a damping term for the vertical velocity when the
2542	! vertical Courant number is too large. This was found to be preferable to
2543	! decreasing the timestep or increasing the diffusion in real-data applications
2544	! that produced potentially-unstable large vertical velocities because of
2545	! unphysically large heating rates coming from the cumulus parameterization
2546	! schemes run at moderately high resolutions (dx ~ O(10) km).
2547	!
2548	! Additionally, w_damp returns the maximum cfl values due to vertical motion and
2549	! horizontal motion. These values are returned via the max_vert_cfl and
2550	! max_horiz_cfl variables. (Added by T. Hutchinson, WSI, 3/5/2007)
2551	!
2552	!</DESCRIPTION>
2553
2554	itf=MIN(ite,ide-1)
2555	jtf=MIN(jte,jde-1)
2556
2557	some = 0
2558	max_vert_cfl = 0.
2559	max_horiz_cfl = 0.
2560	total = 0
2561
2562	IF(config_flags%map_proj == PROJ_CASSINI ) then
2563	msfxffl = 1.0/COS(config_flags%fft_filter_lat*degrad)
2564	END IF
2565
2566	IF ( config_flags%w_damping == 1 ) THEN
2567	DO j = jts,jtf
2568
2569	DO k = 2, kde-1
2570	DO i = its,itf
2571	#if 1
2572	IF(config_flags%map_proj == PROJ_CASSINI ) then
2573	msfuxt = MIN(msfux(i,j), msfxffl)
2574	ELSE
2575	msfuxt = msfux(i,j)
2576	END IF
2577	vert_cfl = abs(ww(i,k,j)/mut(i,j)rdnw(k)dt)
2578
2579	IF ( vert_cfl > max_vert_cfl ) THEN
2580	max_vert_cfl = vert_cfl ; maxi = i ; maxj = j ; maxk = k
2581	maxdub = w(i,k,j) ; maxdeta = -1./rdnw(k)
2582	ENDIF
2583
2584	horiz_cfl = max( abs(u(i,k,j) * rdx * msfuxt * dt), &
2585	abs(v(i,k,j) * rdy * msfvy(i,j) * dt) )
2586	if (horiz_cfl > max_horiz_cfl) then
2587	max_horiz_cfl = horiz_cfl
2588	endif
2589
2590	if(vert_cfl .gt. w_beta)then
2591	#else
2592	! restructure to get rid of divide
2593	!
2594	! This had been used for efficiency, but with the addition of returning the cfl values,
2595	! the old version (above) was reinstated. (T. Hutchinson, 3/5/2007)
2596	!
2597	cf_n = abs(ww(i,k,j)rdnw(k)dt)
2598	cf_d = abs(mut(i,j))
2599	if(cf_n .gt. cf_d*w_beta )then
2600	#endif
2601
2602	WRITE(temp,*)i,j,k,' vert_cfl,w,d(eta)=',vert_cfl,w(i,k,j),-1./rdnw(k)
2603	CALL wrf_debug ( 100 , TRIM(temp) )
2604	if ( vert_cfl > 2. ) some = some + 1
2605	rw_tend(i,k,j) = rw_tend(i,k,j)-sign(1.,w(i,k,j))w_alpha(vert_cfl-w_beta)*mut(i,j)
2606	endif
2607	END DO
2608	ENDDO
2609	ENDDO
2610	ELSE
2611	! just print
2612	DO j = jts,jtf
2613
2614	DO k = 2, kde-1
2615	DO i = its,itf
2616
2617	#if 1
2618	IF(config_flags%map_proj == PROJ_CASSINI ) then
2619	msfuxt = MIN(msfux(i,j), msfxffl)
2620	ELSE
2621	msfuxt = msfux(i,j)
2622	END IF
2623	vert_cfl = abs(ww(i,k,j)/mut(i,j)rdnw(k)dt)
2624
2625	IF ( vert_cfl > max_vert_cfl ) THEN
2626	max_vert_cfl = vert_cfl ; maxi = i ; maxj = j ; maxk = k
2627	maxdub = w(i,k,j) ; maxdeta = -1./rdnw(k)
2628	ENDIF
2629
2630	horiz_cfl = max( abs(u(i,k,j) * rdx * msfuxt * dt), &
2631	abs(v(i,k,j) * rdy * msfvy(i,j) * dt) )
2632
2633	if (horiz_cfl > max_horiz_cfl) then
2634	max_horiz_cfl = horiz_cfl
2635	endif
2636
2637	if(vert_cfl .gt. w_beta)then
2638	#else
2639	! restructure to get rid of divide
2640	!
2641	! This had been used for efficiency, but with the addition of returning the cfl values,
2642	! the old version (above) was reinstated. (T. Hutchinson, 3/5/2007)
2643	!
2644	cf_n = abs(ww(i,k,j)rdnw(k)dt)
2645	cf_d = abs(mut(i,j))
2646	if(cf_n .gt. cf_d*w_beta )then
2647	#endif
2648	WRITE(temp,*)i,j,k,' vert_cfl,w,d(eta)=',vert_cfl,w(i,k,j),-1./rdnw(k)
2649	CALL wrf_debug ( 100 , TRIM(temp) )
2650	if ( vert_cfl > 2. ) some = some + 1
2651	endif
2652	END DO
2653	ENDDO
2654	ENDDO
2655	ENDIF
2656	IF ( some .GT. 0 ) THEN
2657	CALL get_current_time_string( time_str )
2658	CALL get_current_grid_name( grid_str )
2659	WRITE(temp,*)some, &
2660	' points exceeded cfl=2 in domain '//TRIM(grid_str)//' at time '//TRIM(time_str)//' hours'
2661	CALL wrf_debug ( 0 , TRIM(temp) )
2662	WRITE(temp,*)'MAX AT i,j,k: ',maxi,maxj,maxk,' vert_cfl,w,d(eta)=',max_vert_cfl, &
2663	maxdub,maxdeta
2664	CALL wrf_debug ( 0 , TRIM(temp) )
2665	ENDIF
2666
2667	END SUBROUTINE w_damp
2668
2669	!-------------------------------------------------------------------------------
2670
2671	SUBROUTINE horizontal_diffusion ( name, field, tendency, mu, &
2672	config_flags, &
2673	msfux, msfuy, msfvx, msfvx_inv, &
2674	msfvy, msftx, msfty, &
2675	khdif, xkmhd, rdx, rdy, &
2676	ids, ide, jds, jde, kds, kde, &
2677	ims, ime, jms, jme, kms, kme, &
2678	its, ite, jts, jte, kts, kte )
2679
2680	IMPLICIT NONE
2681
2682	! Input data
2683
2684	TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
2685
2686	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
2687	ims, ime, jms, jme, kms, kme, &
2688	its, ite, jts, jte, kts, kte
2689
2690	CHARACTER(LEN=1) , INTENT(IN ) :: name
2691
2692	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: field, xkmhd
2693
2694	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency
2695
2696	REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mu
2697
2698	REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: msfux, &
2699	msfuy, &
2700	msfvx, &
2701	msfvx_inv, &
2702	msfvy, &
2703	msftx, &
2704	msfty
2705
2706	REAL , INTENT(IN ) :: rdx, &
2707	rdy, &
2708	khdif
2709
2710	! Local data
2711
2712	INTEGER :: i, j, k, itf, jtf, ktf
2713
2714	INTEGER :: i_start, i_end, j_start, j_end
2715
2716	REAL :: mrdx, mkrdxm, mkrdxp, &
2717	mrdy, mkrdym, mkrdyp
2718
2719	LOGICAL :: specified
2720
2721	!<DESCRIPTION>
2722	!
2723	! horizontal_diffusion computes the horizontal diffusion tendency
2724	! on model horizontal coordinate surfaces.
2725	!
2726	!</DESCRIPTION>
2727
2728	specified = .false.
2729	if(config_flags%specified .or. config_flags%nested) specified = .true.
2730
2731	ktf=MIN(kte,kde-1)
2732
2733	IF (name .EQ. 'u') THEN
2734
2735	i_start = its
2736	i_end = ite
2737	j_start = jts
2738	j_end = MIN(jte,jde-1)
2739
2740	IF ( config_flags%open_xs .or. specified ) i_start = MAX(ids+1,its)
2741	IF ( config_flags%open_xe .or. specified ) i_end = MIN(ide-1,ite)
2742	IF ( config_flags%open_ys .or. specified ) j_start = MAX(jds+1,jts)
2743	IF ( config_flags%open_ye .or. specified ) j_end = MIN(jde-2,jte)
2744	IF ( config_flags%periodic_x ) i_start = its
2745	IF ( config_flags%periodic_x ) i_end = ite
2746
2747
2748	DO j = j_start, j_end
2749	DO k=kts,ktf
2750	DO i = i_start, i_end
2751
2752	! The interior is grad: (m_xd/dx), the exterior is div: (m_xm_y*d/dx(/m_y))
2753	! setting up different averagings of m^2 partial d/dX and m^2 partial d/dY
2754
2755	mkrdxm=(msftx(i-1,j)/msfty(i-1,j))mu(i-1,j)xkmhd(i-1,k,j)*rdx
2756	mkrdxp=(msftx(i,j)/msfty(i,j))mu(i,j)xkmhd(i,k,j)*rdx
2757	mrdx=msfux(i,j)msfuy(i,j)rdx
2758	mkrdym=( (msfuy(i,j)+msfuy(i,j-1))/(msfux(i,j)+msfux(i,j-1)) )* &
2759	0.25(mu(i,j)+mu(i,j-1)+mu(i-1,j-1)+mu(i-1,j)) &
2760	0.25(xkmhd(i,k,j)+xkmhd(i,k,j-1)+xkmhd(i-1,k,j-1)+xkmhd(i-1,k,j))rdy
2761	mkrdyp=( (msfuy(i,j)+msfuy(i,j+1))/(msfux(i,j)+msfux(i,j+1)) )* &
2762	0.25(mu(i,j)+mu(i,j+1)+mu(i-1,j+1)+mu(i-1,j)) &
2763	0.25(xkmhd(i,k,j)+xkmhd(i,k,j+1)+xkmhd(i-1,k,j+1)+xkmhd(i-1,k,j))rdy
2764	! need to do four-corners (t) for diffusion coefficient as there are
2765	! no values at u,v points
2766	! msfuy - has to be y as part of d/dY
2767	! has to be u as we're at a u point
2768	mrdy=msfux(i,j)msfuy(i,j)rdy
2769
2770	! correctly averaged version of rho~ * m^2 *
2771	! [partial d/dX(partial du^/dX) + partial d/dY(partial du^/dY)]
2772	tendency(i,k,j)=tendency(i,k,j)+( &
2773	mrdx(mkrdxp(field(i+1,k,j)-field(i ,k,j)) &
2774	-mkrdxm*(field(i ,k,j)-field(i-1,k,j))) &
2775	+mrdy(mkrdyp(field(i,k,j+1)-field(i,k,j )) &
2776	-mkrdym*(field(i,k,j )-field(i,k,j-1))))
2777	ENDDO
2778	ENDDO
2779	ENDDO
2780
2781	ELSE IF (name .EQ. 'v')THEN
2782
2783	i_start = its
2784	i_end = MIN(ite,ide-1)
2785	j_start = jts
2786	j_end = jte
2787
2788	IF ( config_flags%open_xs .or. specified ) i_start = MAX(ids+1,its)
2789	IF ( config_flags%open_xe .or. specified ) i_end = MIN(ide-2,ite)
2790	IF ( config_flags%open_ys .or. specified ) j_start = MAX(jds+1,jts)
2791	IF ( config_flags%open_ye .or. specified ) j_end = MIN(jde-1,jte)
2792	IF ( config_flags%periodic_x ) i_start = its
2793	IF ( config_flags%periodic_x ) i_end = MIN(ite,ide-1)
2794	IF ( config_flags%polar ) j_start = MAX(jds+1,jts)
2795	IF ( config_flags%polar ) j_end = MIN(jde-1,jte)
2796
2797	DO j = j_start, j_end
2798	DO k=kts,ktf
2799	DO i = i_start, i_end
2800
2801	mkrdxm=( (msfvx(i,j)+msfvx(i-1,j))/(msfvy(i,j)+msfvy(i-1,j)) )* &
2802	0.25(mu(i,j)+mu(i,j-1)+mu(i-1,j-1)+mu(i-1,j)) &
2803	0.25(xkmhd(i,k,j)+xkmhd(i,k,j-1)+xkmhd(i-1,k,j-1)+xkmhd(i-1,k,j))rdx
2804	mkrdxp=( (msfvx(i,j)+msfvx(i+1,j))/(msfvy(i,j)+msfvy(i+1,j)) )* &
2805	0.25(mu(i,j)+mu(i,j-1)+mu(i+1,j-1)+mu(i+1,j)) &
2806	0.25(xkmhd(i,k,j)+xkmhd(i,k,j-1)+xkmhd(i+1,k,j-1)+xkmhd(i+1,k,j))rdx
2807	mrdx=msfvx(i,j)msfvy(i,j)rdx
2808	mkrdym=(msfty(i,j-1)/msftx(i,j-1))xkmhd(i,k,j-1)rdy
2809	mkrdyp=(msfty(i,j)/msftx(i,j))xkmhd(i,k,j)rdy
2810	mrdy=msfvx(i,j)msfvy(i,j)rdy
2811
2812	tendency(i,k,j)=tendency(i,k,j)+( &
2813	mrdx(mkrdxp(field(i+1,k,j)-field(i ,k,j)) &
2814	-mkrdxm*(field(i ,k,j)-field(i-1,k,j))) &
2815	+mrdy(mkrdyp(field(i,k,j+1)-field(i,k,j )) &
2816	-mkrdym*(field(i,k,j )-field(i,k,j-1))))
2817	ENDDO
2818	ENDDO
2819	ENDDO
2820
2821	ELSE IF (name .EQ. 'w')THEN
2822
2823	i_start = its
2824	i_end = MIN(ite,ide-1)
2825	j_start = jts
2826	j_end = MIN(jte,jde-1)
2827
2828	IF ( config_flags%open_xs .or. specified ) i_start = MAX(ids+1,its)
2829	IF ( config_flags%open_xe .or. specified ) i_end = MIN(ide-2,ite)
2830	IF ( config_flags%open_ys .or. specified ) j_start = MAX(jds+1,jts)
2831	IF ( config_flags%open_ye .or. specified ) j_end = MIN(jde-2,jte)
2832	IF ( config_flags%periodic_x ) i_start = its
2833	IF ( config_flags%periodic_x ) i_end = MIN(ite,ide-1)
2834
2835	DO j = j_start, j_end
2836	DO k=kts+1,ktf
2837	DO i = i_start, i_end
2838
2839	mkrdxm=(msfux(i,j)/msfuy(i,j))* &
2840	0.25(mu(i,j)+mu(i-1,j)+mu(i,j)+mu(i-1,j)) &
2841	0.25(xkmhd(i,k,j)+xkmhd(i-1,k,j)+xkmhd(i,k-1,j)+xkmhd(i-1,k-1,j))rdx
2842	mkrdxp=(msfux(i+1,j)/msfuy(i+1,j))* &
2843	0.25(mu(i+1,j)+mu(i,j)+mu(i+1,j)+mu(i,j)) &
2844	0.25(xkmhd(i+1,k,j)+xkmhd(i,k,j)+xkmhd(i+1,k-1,j)+xkmhd(i,k-1,j))rdx
2845	mrdx=msftx(i,j)msfty(i,j)rdx
2846	! mkrdym=(msfvy(i,j)/msfvx(i,j))* &
2847	mkrdym=(msfvy(i,j)msfvx_inv(i,j)) &
2848	0.25(mu(i,j)+mu(i,j-1)+mu(i,j)+mu(i,j-1)) &
2849	0.25(xkmhd(i,k,j)+xkmhd(i,k,j-1)+xkmhd(i,k-1,j)+xkmhd(i,k-1,j-1))rdy
2850	! mkrdyp=(msfvy(i,j+1)/msfvx(i,j+1))* &
2851	mkrdyp=(msfvy(i,j+1)msfvx_inv(i,j+1)) &
2852	0.25(mu(i,j+1)+mu(i,j)+mu(i,j+1)+mu(i,j)) &
2853	0.25(xkmhd(i,k,j+1)+xkmhd(i,k,j)+xkmhd(i,k-1,j+1)+xkmhd(i,k-1,j))rdy
2854	mrdy=msftx(i,j)msfty(i,j)rdy
2855
2856	tendency(i,k,j)=tendency(i,k,j)+( &
2857	mrdx(mkrdxp(field(i+1,k,j)-field(i ,k,j)) &
2858	-mkrdxm*(field(i ,k,j)-field(i-1,k,j))) &
2859	+mrdy(mkrdyp(field(i,k,j+1)-field(i,k,j )) &
2860	-mkrdym*(field(i,k,j )-field(i,k,j-1))))
2861	ENDDO
2862	ENDDO
2863	ENDDO
2864
2865	ELSE
2866
2867
2868	i_start = its
2869	i_end = MIN(ite,ide-1)
2870	j_start = jts
2871	j_end = MIN(jte,jde-1)
2872
2873	IF ( config_flags%open_xs .or. specified ) i_start = MAX(ids+1,its)
2874	IF ( config_flags%open_xe .or. specified ) i_end = MIN(ide-2,ite)
2875	IF ( config_flags%open_ys .or. specified ) j_start = MAX(jds+1,jts)
2876	IF ( config_flags%open_ye .or. specified ) j_end = MIN(jde-2,jte)
2877	IF ( config_flags%periodic_x ) i_start = its
2878	IF ( config_flags%periodic_x ) i_end = MIN(ite,ide-1)
2879
2880	DO j = j_start, j_end
2881	DO k=kts,ktf
2882	DO i = i_start, i_end
2883
2884	mkrdxm=(msfux(i,j)/msfuy(i,j))0.5(xkmhd(i,k,j)+xkmhd(i-1,k,j))0.5(mu(i,j)+mu(i-1,j))*rdx
2885	mkrdxp=(msfux(i+1,j)/msfuy(i+1,j))0.5(xkmhd(i+1,k,j)+xkmhd(i,k,j))0.5(mu(i+1,j)+mu(i,j))*rdx
2886	mrdx=msftx(i,j)msfty(i,j)rdx
2887	! mkrdym=(msfvy(i,j)/msfvx(i,j))0.5(xkmhd(i,k,j)+xkmhd(i,k,j-1))0.5(mu(i,j)+mu(i,j-1))*rdy
2888	mkrdym=(msfvy(i,j)msfvx_inv(i,j))0.5(xkmhd(i,k,j)+xkmhd(i,k,j-1))0.5(mu(i,j)+mu(i,j-1))rdy
2889	! mkrdyp=(msfvy(i,j+1)/msfvx(i,j+1))0.5(xkmhd(i,k,j+1)+xkmhd(i,k,j))0.5(mu(i,j+1)+mu(i,j))*rdy
2890	mkrdyp=(msfvy(i,j+1)msfvx_inv(i,j+1))0.5(xkmhd(i,k,j+1)+xkmhd(i,k,j))0.5(mu(i,j+1)+mu(i,j))rdy
2891	mrdy=msftx(i,j)msfty(i,j)rdy
2892
2893	tendency(i,k,j)=tendency(i,k,j)+( &
2894	mrdx(mkrdxp(field(i+1,k,j)-field(i ,k,j)) &
2895	-mkrdxm*(field(i ,k,j)-field(i-1,k,j))) &
2896	+mrdy(mkrdyp(field(i,k,j+1)-field(i,k,j )) &
2897	-mkrdym*(field(i,k,j )-field(i,k,j-1))))
2898	ENDDO
2899	ENDDO
2900	ENDDO
2901
2902	ENDIF
2903
2904	END SUBROUTINE horizontal_diffusion
2905
2906	!-----------------------------------------------------------------------------------------
2907
2908	SUBROUTINE horizontal_diffusion_3dmp ( name, field, tendency, mu, &
2909	config_flags, base_3d, &
2910	msfux, msfuy, msfvx, msfvx_inv, &
2911	msfvy, msftx, msfty, &
2912	khdif, xkmhd, rdx, rdy, &
2913	ids, ide, jds, jde, kds, kde, &
2914	ims, ime, jms, jme, kms, kme, &
2915	its, ite, jts, jte, kts, kte )
2916
2917	IMPLICIT NONE
2918
2919	! Input data
2920
2921	TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
2922
2923	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
2924	ims, ime, jms, jme, kms, kme, &
2925	its, ite, jts, jte, kts, kte
2926
2927	CHARACTER(LEN=1) , INTENT(IN ) :: name
2928
2929	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: field, &
2930	xkmhd, &
2931	base_3d
2932
2933	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency
2934
2935	REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mu
2936
2937	REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: msfux, &
2938	msfuy, &
2939	msfvx, &
2940	msfvx_inv, &
2941	msfvy, &
2942	msftx, &
2943	msfty
2944
2945	REAL , INTENT(IN ) :: rdx, &
2946	rdy, &
2947	khdif
2948
2949	! Local data
2950
2951	INTEGER :: i, j, k, itf, jtf, ktf
2952
2953	INTEGER :: i_start, i_end, j_start, j_end
2954
2955	REAL :: mrdx, mkrdxm, mkrdxp, &
2956	mrdy, mkrdym, mkrdyp
2957
2958	LOGICAL :: specified
2959
2960	!<DESCRIPTION>
2961	!
2962	! horizontal_diffusion_3dmp computes the horizontal diffusion tendency
2963	! on model horizontal coordinate surfaces. This routine computes diffusion
2964	! a perturbation scalar (field-base_3d).
2965	!
2966	!</DESCRIPTION>
2967
2968	specified = .false.
2969	if(config_flags%specified .or. config_flags%nested) specified = .true.
2970
2971	ktf=MIN(kte,kde-1)
2972
2973	i_start = its
2974	i_end = MIN(ite,ide-1)
2975	j_start = jts
2976	j_end = MIN(jte,jde-1)
2977
2978	IF ( config_flags%open_xs .or. specified ) i_start = MAX(ids+1,its)
2979	IF ( config_flags%open_xe .or. specified ) i_end = MIN(ide-2,ite)
2980	IF ( config_flags%open_ys .or. specified ) j_start = MAX(jds+1,jts)
2981	IF ( config_flags%open_ye .or. specified ) j_end = MIN(jde-2,jte)
2982	IF ( config_flags%periodic_x ) i_start = its
2983	IF ( config_flags%periodic_x ) i_end = MIN(ite,ide-1)
2984
2985	DO j = j_start, j_end
2986	DO k=kts,ktf
2987	DO i = i_start, i_end
2988
2989	mkrdxm=(msfux(i,j)/msfuy(i,j))0.5(xkmhd(i,k,j)+xkmhd(i-1,k,j))0.5(mu(i,j)+mu(i-1,j))*rdx
2990	mkrdxp=(msfux(i+1,j)/msfuy(i+1,j))0.5(xkmhd(i+1,k,j)+xkmhd(i,k,j))0.5(mu(i+1,j)+mu(i,j))*rdx
2991	mrdx=msftx(i,j)msfty(i,j)rdx
2992	! mkrdym=(msfvy(i,j)/msfvx(i,j))0.5(xkmhd(i,k,j)+xkmhd(i,k,j-1))0.5(mu(i,j)+mu(i,j-1))*rdy
2993	! mkrdyp=(msfvy(i,j+1)/msfvx(i,j+1))0.5(xkmhd(i,k,j+1)+xkmhd(i,k,j))0.5(mu(i,j+1)+mu(i,j))*rdy
2994	mkrdym=(msfvy(i,j)msfvx_inv(i,j))0.5(xkmhd(i,k,j)+xkmhd(i,k,j-1))0.5(mu(i,j)+mu(i,j-1))rdy
2995	mkrdyp=(msfvy(i,j+1)msfvx_inv(i,j+1))0.5(xkmhd(i,k,j+1)+xkmhd(i,k,j))0.5(mu(i,j+1)+mu(i,j))rdy
2996	mrdy=msftx(i,j)msfty(i,j)rdy
2997
2998	tendency(i,k,j)=tendency(i,k,j)+( &
2999	mrdx( mkrdxp( field(i+1,k,j) -field(i ,k,j) &
3000	-base_3d(i+1,k,j)+base_3d(i ,k,j) ) &
3001	-mkrdxm*( field(i ,k,j) -field(i-1,k,j) &
3002	-base_3d(i ,k,j)+base_3d(i-1,k,j) ) ) &
3003	+mrdy( mkrdyp( field(i,k,j+1) -field(i,k,j ) &
3004	-base_3d(i,k,j+1)+base_3d(i,k,j ) ) &
3005	-mkrdym*( field(i,k,j ) -field(i,k,j-1) &
3006	-base_3d(i,k,j )+base_3d(i,k,j-1) ) ) &
3007	)
3008	ENDDO
3009	ENDDO
3010	ENDDO
3011
3012	END SUBROUTINE horizontal_diffusion_3dmp
3013
3014	!-----------------------------------------------------------------------------------------
3015
3016	SUBROUTINE vertical_diffusion ( name, field, tendency, &
3017	config_flags, &
3018	alt, mut, rdn, rdnw, kvdif, &
3019	ids, ide, jds, jde, kds, kde, &
3020	ims, ime, jms, jme, kms, kme, &
3021	its, ite, jts, jte, kts, kte )
3022
3023
3024	IMPLICIT NONE
3025
3026	! Input data
3027
3028	TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
3029
3030	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
3031	ims, ime, jms, jme, kms, kme, &
3032	its, ite, jts, jte, kts, kte
3033
3034	CHARACTER(LEN=1) , INTENT(IN ) :: name
3035
3036	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , &
3037	INTENT(IN ) :: field, &
3038	alt
3039
3040	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency
3041
3042	REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mut
3043
3044	REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: rdn, rdnw
3045
3046	REAL , INTENT(IN ) :: kvdif
3047
3048	! Local data
3049
3050	INTEGER :: i, j, k, itf, jtf, ktf
3051	INTEGER :: i_start, i_end, j_start, j_end
3052
3053	REAL , DIMENSION(its:ite, jts:jte) :: vfluxm, vfluxp, zz
3054	REAL , DIMENSION(its:ite, 0:kte+1) :: vflux
3055
3056	REAL :: rdz
3057
3058	LOGICAL :: specified
3059
3060	!<DESCRIPTION>
3061	!
3062	! vertical_diffusion
3063	! computes vertical diffusion tendency.
3064	!
3065	!</DESCRIPTION>
3066
3067	specified = .false.
3068	if(config_flags%specified .or. config_flags%nested) specified = .true.
3069
3070	ktf=MIN(kte,kde-1)
3071
3072	IF (name .EQ. 'w')THEN
3073
3074
3075	i_start = its
3076	i_end = MIN(ite,ide-1)
3077	j_start = jts
3078	j_end = MIN(jte,jde-1)
3079
3080	j_loop_w : DO j = j_start, j_end
3081
3082	DO k=kts,ktf-1
3083	DO i = i_start, i_end
3084	vflux(i,k)= (kvdif/alt(i,k,j))rdnw(k)(field(i,k+1,j)-field(i,k,j))
3085	ENDDO
3086	ENDDO
3087
3088	DO i = i_start, i_end
3089	vflux(i,ktf)=0.
3090	ENDDO
3091
3092	DO k=kts+1,ktf
3093	DO i = i_start, i_end
3094	tendency(i,k,j)=tendency(i,k,j) &
3095	+rdn(k)gg/mut(i,j)/(0.5*(alt(i,k,j)+alt(i,k-1,j))) &
3096	*(vflux(i,k)-vflux(i,k-1))
3097	ENDDO
3098	ENDDO
3099
3100	ENDDO j_loop_w
3101
3102	ELSE IF(name .EQ. 'm')THEN
3103
3104	i_start = its
3105	i_end = MIN(ite,ide-1)
3106	j_start = jts
3107	j_end = MIN(jte,jde-1)
3108
3109	j_loop_s : DO j = j_start, j_end
3110
3111	DO k=kts,ktf-1
3112	DO i = i_start, i_end
3113	vflux(i,k)=kvdifrdn(k+1)/(0.5(alt(i,k,j)+alt(i,k+1,j))) &
3114	*(field(i,k+1,j)-field(i,k,j))
3115	ENDDO
3116	ENDDO
3117
3118	DO i = i_start, i_end
3119	vflux(i,0)=vflux(i,1)
3120	ENDDO
3121
3122	DO i = i_start, i_end
3123	vflux(i,ktf)=0.
3124	ENDDO
3125
3126	DO k=kts,ktf
3127	DO i = i_start, i_end
3128	tendency(i,k,j)=tendency(i,k,j)+g*g/mut(i,j)/alt(i,k,j) &
3129	rdnw(k)(vflux(i,k)-vflux(i,k-1))
3130	ENDDO
3131	ENDDO
3132
3133	ENDDO j_loop_s
3134
3135	ENDIF
3136
3137	END SUBROUTINE vertical_diffusion
3138
3139
3140	!-------------------------------------------------------------------------------
3141
3142	SUBROUTINE vertical_diffusion_mp ( field, tendency, config_flags, &
3143	base, &
3144	alt, mut, rdn, rdnw, kvdif, &
3145	ids, ide, jds, jde, kds, kde, &
3146	ims, ime, jms, jme, kms, kme, &
3147	its, ite, jts, jte, kts, kte )
3148
3149
3150	IMPLICIT NONE
3151
3152	! Input data
3153
3154	TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
3155
3156	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
3157	ims, ime, jms, jme, kms, kme, &
3158	its, ite, jts, jte, kts, kte
3159
3160	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , &
3161	INTENT(IN ) :: field, &
3162	alt
3163
3164	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency
3165
3166	REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mut
3167
3168	REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: rdn, &
3169	rdnw, &
3170	base
3171
3172	REAL , INTENT(IN ) :: kvdif
3173
3174	! Local data
3175
3176	INTEGER :: i, j, k, itf, jtf, ktf
3177	INTEGER :: i_start, i_end, j_start, j_end
3178
3179	REAL , DIMENSION(its:ite, 0:kte+1) :: vflux
3180
3181	REAL :: rdz
3182
3183	LOGICAL :: specified
3184
3185	!<DESCRIPTION>
3186	!
3187	! vertical_diffusion_mp
3188	! computes vertical diffusion tendency of a perturbation variable
3189	! (field-base). Note that base as a 1D (k) field.
3190	!
3191	!</DESCRIPTION>
3192
3193	specified = .false.
3194	if(config_flags%specified .or. config_flags%nested) specified = .true.
3195
3196	ktf=MIN(kte,kde-1)
3197
3198	i_start = its
3199	i_end = MIN(ite,ide-1)
3200	j_start = jts
3201	j_end = MIN(jte,jde-1)
3202
3203	j_loop_s : DO j = j_start, j_end
3204
3205	DO k=kts,ktf-1
3206	DO i = i_start, i_end
3207	vflux(i,k)=kvdifrdn(k+1)/(0.5(alt(i,k,j)+alt(i,k+1,j))) &
3208	*(field(i,k+1,j)-field(i,k,j)-base(k+1)+base(k))
3209	ENDDO
3210	ENDDO
3211
3212	DO i = i_start, i_end
3213	vflux(i,0)=vflux(i,1)
3214	ENDDO
3215
3216	DO i = i_start, i_end
3217	vflux(i,ktf)=0.
3218	ENDDO
3219
3220	DO k=kts,ktf
3221	DO i = i_start, i_end
3222	tendency(i,k,j)=tendency(i,k,j)+g*g/mut(i,j)/alt(i,k,j) &
3223	rdnw(k)(vflux(i,k)-vflux(i,k-1))
3224	ENDDO
3225	ENDDO
3226
3227	ENDDO j_loop_s
3228
3229	END SUBROUTINE vertical_diffusion_mp
3230
3231
3232	!-------------------------------------------------------------------------------
3233
3234	SUBROUTINE vertical_diffusion_3dmp ( field, tendency, config_flags, &
3235	base_3d, &
3236	alt, mut, rdn, rdnw, kvdif, &
3237	ids, ide, jds, jde, kds, kde, &
3238	ims, ime, jms, jme, kms, kme, &
3239	its, ite, jts, jte, kts, kte )
3240
3241
3242	IMPLICIT NONE
3243
3244	! Input data
3245
3246	TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
3247
3248	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
3249	ims, ime, jms, jme, kms, kme, &
3250	its, ite, jts, jte, kts, kte
3251
3252	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , &
3253	INTENT(IN ) :: field, &
3254	alt, &
3255	base_3d
3256
3257	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency
3258
3259	REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mut
3260
3261	REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: rdn, &
3262	rdnw
3263
3264	REAL , INTENT(IN ) :: kvdif
3265
3266	! Local data
3267
3268	INTEGER :: i, j, k, itf, jtf, ktf
3269	INTEGER :: i_start, i_end, j_start, j_end
3270
3271	REAL , DIMENSION(its:ite, 0:kte+1) :: vflux
3272
3273	REAL :: rdz
3274
3275	LOGICAL :: specified
3276
3277	!<DESCRIPTION>
3278	!
3279	! vertical_diffusion_3dmp
3280	! computes vertical diffusion tendency of a perturbation variable
3281	! (field-base_3d).
3282	!
3283	!</DESCRIPTION>
3284
3285	specified = .false.
3286	if(config_flags%specified .or. config_flags%nested) specified = .true.
3287
3288	ktf=MIN(kte,kde-1)
3289
3290	i_start = its
3291	i_end = MIN(ite,ide-1)
3292	j_start = jts
3293	j_end = MIN(jte,jde-1)
3294
3295	j_loop_s : DO j = j_start, j_end
3296
3297	DO k=kts,ktf-1
3298	DO i = i_start, i_end
3299	vflux(i,k)=kvdifrdn(k+1)/(0.5(alt(i,k,j)+alt(i,k+1,j))) &
3300	*( field(i,k+1,j) -field(i,k,j) &
3301	-base_3d(i,k+1,j)+base_3d(i,k,j) )
3302	ENDDO
3303	ENDDO
3304
3305	DO i = i_start, i_end
3306	vflux(i,0)=vflux(i,1)
3307	ENDDO
3308
3309	DO i = i_start, i_end
3310	vflux(i,ktf)=0.
3311	ENDDO
3312
3313	DO k=kts,ktf
3314	DO i = i_start, i_end
3315	tendency(i,k,j)=tendency(i,k,j)+g*g/mut(i,j)/alt(i,k,j) &
3316	rdnw(k)(vflux(i,k)-vflux(i,k-1))
3317	ENDDO
3318	ENDDO
3319
3320	ENDDO j_loop_s
3321
3322	END SUBROUTINE vertical_diffusion_3dmp
3323
3324
3325	!-------------------------------------------------------------------------------
3326
3327
3328	SUBROUTINE vertical_diffusion_u ( field, tendency, &
3329	config_flags, u_base, &
3330	alt, muu, rdn, rdnw, kvdif, &
3331	ids, ide, jds, jde, kds, kde, &
3332	ims, ime, jms, jme, kms, kme, &
3333	its, ite, jts, jte, kts, kte )
3334
3335
3336	IMPLICIT NONE
3337
3338	! Input data
3339
3340	TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
3341
3342	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
3343	ims, ime, jms, jme, kms, kme, &
3344	its, ite, jts, jte, kts, kte
3345
3346	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , &
3347	INTENT(IN ) :: field, &
3348	alt
3349
3350	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency
3351
3352	REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: muu
3353
3354	REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: rdn, rdnw, u_base
3355
3356	REAL , INTENT(IN ) :: kvdif
3357
3358	! Local data
3359
3360	INTEGER :: i, j, k, itf, jtf, ktf
3361	INTEGER :: i_start, i_end, j_start, j_end
3362
3363	REAL , DIMENSION(its:ite, 0:kte+1) :: vflux
3364
3365	REAL :: rdz, zz
3366
3367	LOGICAL :: specified
3368
3369	!<DESCRIPTION>
3370	!
3371	! vertical_diffusion_u computes vertical diffusion tendency for
3372	! the u momentum equation. This routine assumes a constant eddy
3373	! viscosity kvdif.
3374	!
3375	!</DESCRIPTION>
3376
3377	specified = .false.
3378	if(config_flags%specified .or. config_flags%nested) specified = .true.
3379
3380	ktf=MIN(kte,kde-1)
3381
3382	i_start = its
3383	i_end = ite
3384	j_start = jts
3385	j_end = MIN(jte,jde-1)
3386
3387	IF ( config_flags%open_xs .or. specified ) i_start = MAX(ids+1,its)
3388	IF ( config_flags%open_xe .or. specified ) i_end = MIN(ide-1,ite)
3389	IF ( config_flags%periodic_x ) i_start = its
3390	IF ( config_flags%periodic_x ) i_end = ite
3391
3392
3393	j_loop_u : DO j = j_start, j_end
3394
3395	DO k=kts,ktf-1
3396	DO i = i_start, i_end
3397	vflux(i,k)=kvdifrdn(k+1)/(0.25( alt(i ,k ,j) &
3398	+alt(i-1,k ,j) &
3399	+alt(i ,k+1,j) &
3400	+alt(i-1,k+1,j) ) ) &
3401	*(field(i,k+1,j)-field(i,k,j) &
3402	-u_base(k+1) +u_base(k) )
3403	ENDDO
3404	ENDDO
3405
3406	DO i = i_start, i_end
3407	vflux(i,0)=vflux(i,1)
3408	ENDDO
3409
3410	DO i = i_start, i_end
3411	vflux(i,ktf)=0.
3412	ENDDO
3413
3414	DO k=kts,ktf-1
3415	DO i = i_start, i_end
3416	tendency(i,k,j)=tendency(i,k,j)+ &
3417	ggrdnw(k)/muu(i,j)/(0.5(alt(i-1,k,j)+alt(i,k,j))) &
3418	(vflux(i,k)-vflux(i,k-1))
3419	ENDDO
3420	ENDDO
3421
3422	ENDDO j_loop_u
3423
3424	END SUBROUTINE vertical_diffusion_u
3425
3426	!-------------------------------------------------------------------------------
3427
3428
3429	SUBROUTINE vertical_diffusion_v ( field, tendency, &
3430	config_flags, v_base, &
3431	alt, muv, rdn, rdnw, kvdif, &
3432	ids, ide, jds, jde, kds, kde, &
3433	ims, ime, jms, jme, kms, kme, &
3434	its, ite, jts, jte, kts, kte )
3435
3436
3437	IMPLICIT NONE
3438
3439	! Input data
3440
3441	TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
3442
3443	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
3444	ims, ime, jms, jme, kms, kme, &
3445	its, ite, jts, jte, kts, kte
3446
3447	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , &
3448	INTENT(IN ) :: field, &
3449	alt
3450	REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: rdn, rdnw, v_base
3451
3452	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency
3453
3454	REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: muv
3455
3456	REAL , INTENT(IN ) :: kvdif
3457
3458	! Local data
3459
3460	INTEGER :: i, j, k, itf, jtf, ktf, jm1
3461	INTEGER :: i_start, i_end, j_start, j_end
3462
3463	REAL , DIMENSION(its:ite, 0:kte+1) :: vflux
3464
3465	REAL :: rdz, zz
3466
3467	LOGICAL :: specified
3468
3469	!<DESCRIPTION>
3470	!
3471	! vertical_diffusion_v computes vertical diffusion tendency for
3472	! the v momentum equation. This routine assumes a constant eddy
3473	! viscosity kvdif.
3474	!
3475	!</DESCRIPTION>
3476
3477	specified = .false.
3478	if(config_flags%specified .or. config_flags%nested) specified = .true.
3479
3480	ktf=MIN(kte,kde-1)
3481
3482	i_start = its
3483	i_end = MIN(ite,ide-1)
3484	j_start = jts
3485	j_end = MIN(jte,jde-1)
3486
3487	IF ( config_flags%open_ys .or. specified ) j_start = MAX(jds+1,jts)
3488	IF ( config_flags%open_ye .or. specified ) j_end = MIN(jde-1,jte)
3489
3490	j_loop_v : DO j = j_start, j_end
3491	! jm1 = max(j-1,1)
3492	jm1 = j-1
3493
3494	DO k=kts,ktf-1
3495	DO i = i_start, i_end
3496	vflux(i,k)=kvdifrdn(k+1)/(0.25( alt(i,k ,j ) &
3497	+alt(i,k ,jm1) &
3498	+alt(i,k+1,j ) &
3499	+alt(i,k+1,jm1) ) ) &
3500	*(field(i,k+1,j)-field(i,k,j) &
3501	-v_base(k+1) +v_base(k) )
3502	ENDDO
3503	ENDDO
3504
3505	DO i = i_start, i_end
3506	vflux(i,0)=vflux(i,1)
3507	ENDDO
3508
3509	DO i = i_start, i_end
3510	vflux(i,ktf)=0.
3511	ENDDO
3512
3513	DO k=kts,ktf-1
3514	DO i = i_start, i_end
3515	tendency(i,k,j)=tendency(i,k,j)+ &
3516	ggrdnw(k)/muv(i,j)/(0.5(alt(i,k,jm1)+alt(i,k,j))) &
3517	(vflux(i,k)-vflux(i,k-1))
3518	ENDDO
3519	ENDDO
3520
3521	ENDDO j_loop_v
3522
3523	END SUBROUTINE vertical_diffusion_v
3524
3525	!*************** end new mass coordinate routines
3526
3527	!-------------------------------------------------------------------------------
3528
3529	SUBROUTINE calculate_full ( rfield, rfieldb, rfieldp, &
3530	ids, ide, jds, jde, kds, kde, &
3531	ims, ime, jms, jme, kms, kme, &
3532	its, ite, jts, jte, kts, kte )
3533
3534	IMPLICIT NONE
3535
3536	! Input data
3537
3538	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
3539	ims, ime, jms, jme, kms, kme, &
3540	its, ite, jts, jte, kts, kte
3541
3542	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: rfieldb, &
3543	rfieldp
3544
3545	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT ) :: rfield
3546
3547	! Local indices.
3548
3549	INTEGER :: i, j, k, itf, jtf, ktf
3550
3551	!<DESCRIPTION>
3552	!
3553	! calculate_full
3554	! calculates full 3D field from pertubation and base field.
3555	!
3556	!</DESCRIPTION>
3557
3558	itf=MIN(ite,ide-1)
3559	jtf=MIN(jte,jde-1)
3560	ktf=MIN(kte,kde-1)
3561
3562	DO j=jts,jtf
3563	DO k=kts,ktf
3564	DO i=its,itf
3565	rfield(i,k,j)=rfieldb(i,k,j)+rfieldp(i,k,j)
3566	ENDDO
3567	ENDDO
3568	ENDDO
3569
3570	END SUBROUTINE calculate_full
3571
3572	!------------------------------------------------------------------------------
3573
3574	SUBROUTINE coriolis ( ru, rv, rw, ru_tend, rv_tend, rw_tend, &
3575	config_flags, &
3576	msftx, msfty, msfux, msfuy, &
3577	msfvx, msfvy, &
3578	f, e, sina, cosa, fzm, fzp, &
3579	ids, ide, jds, jde, kds, kde, &
3580	ims, ime, jms, jme, kms, kme, &
3581	its, ite, jts, jte, kts, kte )
3582
3583	IMPLICIT NONE
3584
3585	! Input data
3586
3587	TYPE(grid_config_rec_type) , INTENT(IN ) :: config_flags
3588
3589	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
3590	ims, ime, jms, jme, kms, kme, &
3591	its, ite, jts, jte, kts, kte
3592
3593	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: ru_tend, &
3594	rv_tend, &
3595	rw_tend
3596	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: ru, &
3597	rv, &
3598	rw
3599
3600	REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: msfux, &
3601	msfuy, &
3602	msfvx, &
3603	msfvy, &
3604	msftx, &
3605	msfty
3606
3607	REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: f, &
3608	e, &
3609	sina, &
3610	cosa
3611
3612	REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fzm, &
3613	fzp
3614
3615	! Local indices.
3616
3617	INTEGER :: i, j , k, ktf
3618	INTEGER :: i_start, i_end, j_start, j_end
3619
3620	LOGICAL :: specified
3621
3622	!<DESCRIPTION>
3623	!
3624	! coriolis calculates the large timestep tendency terms in the
3625	! u, v, and w momentum equations arise from the coriolis force.
3626	!
3627	!</DESCRIPTION>
3628
3629	specified = .false.
3630	if(config_flags%specified .or. config_flags%nested) specified = .true.
3631
3632	ktf=MIN(kte,kde-1)
3633
3634	! coriolis for u-momentum equation
3635
3636	! Notes on map scale factor
3637	! cosa, sina are related to rotating the coordinate frame if desired
3638	! generally sina=0, cosa=1
3639	! ADT eqn 44, RHS terms 6 and 7: -2 mu w omega cos(lat)/my
3640	! + 2 mu v omega sin(lat)/my
3641	! Define f=2 omega sin(lat), e=2 omega cos(lat)
3642	! => terms are: -e mu w / my + f mu v / my
3643	! rv = mu v / mx ; rw = mu w / my
3644	! => terms are: -e rw + f rv *mx / my
3645
3646	i_start = its
3647	i_end = ite
3648	IF ( config_flags%open_xs .or. specified .or. &
3649	config_flags%nested) i_start = MAX(ids+1,its)
3650	IF ( config_flags%open_xe .or. specified .or. &
3651	config_flags%nested) i_end = MIN(ide-1,ite)
3652	IF ( config_flags%periodic_x ) i_start = its
3653	IF ( config_flags%periodic_x ) i_end = ite
3654
3655	DO j = jts, MIN(jte,jde-1)
3656
3657	DO k=kts,ktf
3658	DO i = i_start, i_end
3659
3660	ru_tend(i,k,j)=ru_tend(i,k,j) + (msfux(i,j)/msfuy(i,j))0.5(f(i,j)+f(i-1,j)) &
3661	0.25(rv(i-1,k,j+1)+rv(i,k,j+1)+rv(i-1,k,j)+rv(i,k,j)) &
3662	- 0.5(e(i,j)+e(i-1,j))0.5*(cosa(i,j)+cosa(i-1,j)) &
3663	0.25(rw(i-1,k+1,j)+rw(i-1,k,j)+rw(i,k+1,j)+rw(i,k,j))
3664
3665	ENDDO
3666	ENDDO
3667
3668	! boundary loops for coriolis not needed for open bdy (commented out 20100611 JD)
3669	! IF ( (config_flags%open_xs) .and. (its == ids) ) THEN
3670
3671	! DO k=kts,ktf
3672	!
3673	! ru_tend(its,k,j)=ru_tend(its,k,j) + (msfux(its,j)/msfuy(its,j))0.5(f(its,j)+f(its,j)) &
3674	! 0.25(rv(its,k,j+1)+rv(its,k,j+1)+rv(its,k,j)+rv(its,k,j)) &
3675	! - 0.5(e(its,j)+e(its,j))0.5*(cosa(its,j)+cosa(its,j)) &
3676	! 0.25(rw(its,k+1,j)+rw(its,k,j)+rw(its,k+1,j)+rw(its,k,j))
3677
3678	! ENDDO
3679
3680	! ENDIF
3681
3682	! IF ( (config_flags%open_xe) .and. (ite == ide) ) THEN
3683
3684	! DO k=kts,ktf
3685	!
3686	! ru_tend(ite,k,j)=ru_tend(ite,k,j) + (msfux(ite,j)/msfuy(ite,j))0.5(f(ite-1,j)+f(ite-1,j)) &
3687	! 0.25(rv(ite-1,k,j+1)+rv(ite-1,k,j+1)+rv(ite-1,k,j)+rv(ite-1,k,j)) &
3688	! - 0.5(e(ite-1,j)+e(ite-1,j))0.5*(cosa(ite-1,j)+cosa(ite-1,j)) &
3689	! 0.25(rw(ite-1,k+1,j)+rw(ite-1,k,j)+rw(ite-1,k+1,j)+rw(ite-1,k,j))
3690
3691	! ENDDO
3692
3693	! ENDIF
3694
3695	ENDDO
3696
3697	! coriolis term for v-momentum equation
3698
3699	! Notes on map scale factors
3700	! ADT eqn 45, RHS terms 6 and 6b [0 for sina=0]: -2 mu u omega sin(lat)/mx + ?
3701	! Define f=2 omega sin(lat), e=2 omega cos(lat)
3702	! => terms are: -f mu u / mx
3703	! ru = mu u / my ; rw = mu w / my
3704	! => terms are: -f ru *my / mx + ?
3705
3706	j_start = jts
3707	j_end = jte
3708
3709	IF ( config_flags%open_ys .or. specified .or. &
3710	config_flags%nested .or. config_flags%polar) j_start = MAX(jds+1,jts)
3711	IF ( config_flags%open_ye .or. specified .or. &
3712	config_flags%nested .or. config_flags%polar) j_end = MIN(jde-1,jte)
3713
3714	! boundary loops for coriolis not needed for open bdy (commented out 20100611 JD)
3715	! IF ( (config_flags%open_ys) .and. (jts == jds) ) THEN
3716
3717	! DO k=kts,ktf
3718	! DO i=its,MIN(ide-1,ite)
3719	!
3720	! rv_tend(i,k,jts)=rv_tend(i,k,jts) - (msfvy(i,jts)/msfvx(i,jts))0.5(f(i,jts)+f(i,jts)) &
3721	! 0.25(ru(i,k,jts)+ru(i+1,k,jts)+ru(i,k,jts)+ru(i+1,k,jts)) &
3722	! + (msfvy(i,jts)/msfvx(i,jts))0.5(e(i,jts)+e(i,jts))0.5(sina(i,jts)+sina(i,jts)) &
3723	! 0.25(rw(i,k+1,jts)+rw(i,k,jts)+rw(i,k+1,jts)+rw(i,k,jts))
3724
3725	! ENDDO
3726	! ENDDO
3727
3728	! ENDIF
3729
3730	DO j=j_start, j_end
3731	DO k=kts,ktf
3732	DO i=its,MIN(ide-1,ite)
3733
3734	rv_tend(i,k,j)=rv_tend(i,k,j) - (msfvy(i,j)/msfvx(i,j))0.5(f(i,j)+f(i,j-1)) &
3735	0.25(ru(i,k,j)+ru(i+1,k,j)+ru(i,k,j-1)+ru(i+1,k,j-1)) &
3736	+ (msfvy(i,j)/msfvx(i,j))0.5(e(i,j)+e(i,j-1))0.5(sina(i,j)+sina(i,j-1)) &
3737	0.25(rw(i,k+1,j-1)+rw(i,k,j-1)+rw(i,k+1,j)+rw(i,k,j))
3738
3739	ENDDO
3740	ENDDO
3741	ENDDO
3742
3743
3744	! boundary loops for coriolis not needed for open bdy (commented out 20100611 JD)
3745	! IF ( (config_flags%open_ye) .and. (jte == jde) ) THEN
3746
3747	! DO k=kts,ktf
3748	! DO i=its,MIN(ide-1,ite)
3749	!
3750	! rv_tend(i,k,jte)=rv_tend(i,k,jte) - (msfvy(i,jte)/msfvx(i,jte))0.5(f(i,jte-1)+f(i,jte-1)) &
3751	! 0.25(ru(i,k,jte-1)+ru(i+1,k,jte-1)+ru(i,k,jte-1)+ru(i+1,k,jte-1)) &
3752	! + (msfvy(i,jte)/msfvx(i,jte))0.5(e(i,jte-1)+e(i,jte-1))0.5(sina(i,jte-1)+sina(i,jte-1)) &
3753	! 0.25(rw(i,k+1,jte-1)+rw(i,k,jte-1)+rw(i,k+1,jte-1)+rw(i,k,jte-1))
3754
3755	! ENDDO
3756	! ENDDO
3757
3758	! ENDIF
3759
3760	! coriolis term for w-mometum
3761
3762	! Notes on map scale factors
3763	! ADT eqn 46/my, RHS terms 5 and 5b [0 for sina=0]: 2 mu u omega cos(lat)/my +?
3764	! Define e=2 omega cos(lat)
3765	! => terms are: e mu u / my + ???
3766	! ru = mu u / my ; ru = mu v / mx
3767	! => terms are: e ru + ???
3768
3769	DO j=jts,MIN(jte, jde-1)
3770	DO k=kts+1,ktf
3771	DO i=its,MIN(ite, ide-1)
3772
3773	rw_tend(i,k,j)=rw_tend(i,k,j) + e(i,j)* &
3774	(cosa(i,j)0.5(fzm(k)*(ru(i,k,j)+ru(i+1,k,j)) &
3775	+fzp(k)*(ru(i,k-1,j)+ru(i+1,k-1,j))) &
3776	-(msftx(i,j)/msfty(i,j))* &
3777	sina(i,j)0.5(fzm(k)*(rv(i,k,j)+rv(i,k,j+1)) &
3778	+fzp(k)*(rv(i,k-1,j)+rv(i,k-1,j+1))))
3779
3780	ENDDO
3781	ENDDO
3782	ENDDO
3783
3784	END SUBROUTINE coriolis
3785
3786	!------------------------------------------------------------------------------
3787
3788	SUBROUTINE perturbation_coriolis ( ru_in, rv_in, rw, ru_tend, rv_tend, rw_tend, &
3789	config_flags, &
3790	u_base, v_base, z_base, &
3791	muu, muv, phb, ph, &
3792	msftx, msfty, msfux, msfuy, msfvx, msfvy, &
3793	f, e, sina, cosa, fzm, fzp, &
3794	ids, ide, jds, jde, kds, kde, &
3795	ims, ime, jms, jme, kms, kme, &
3796	its, ite, jts, jte, kts, kte )
3797
3798	IMPLICIT NONE
3799
3800	! Input data
3801
3802	TYPE(grid_config_rec_type) , INTENT(IN ) :: config_flags
3803
3804	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
3805	ims, ime, jms, jme, kms, kme, &
3806	its, ite, jts, jte, kts, kte
3807
3808	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: ru_tend, &
3809	rv_tend, &
3810	rw_tend
3811	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: ru_in, &
3812	rv_in, &
3813	rw, &
3814	ph, &
3815	phb
3816
3817
3818	REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: msfux, &
3819	msfuy, &
3820	msfvx, &
3821	msfvy, &
3822	msftx, &
3823	msfty
3824
3825	REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: f, &
3826	e, &
3827	sina, &
3828	cosa
3829
3830	REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: muu, &
3831	muv
3832
3833
3834	REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fzm, &
3835	fzp
3836
3837	REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: u_base, &
3838	v_base, &
3839	z_base
3840
3841	! Local storage
3842
3843	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) :: ru, &
3844	rv
3845
3846	REAL :: z_at_u, z_at_v, wkp1, wk, wkm1
3847
3848	! Local indices.
3849
3850	INTEGER :: i, j , k, ktf
3851	INTEGER :: i_start, i_end, j_start, j_end
3852
3853	LOGICAL :: specified
3854
3855	!<DESCRIPTION>
3856	!
3857	! perturbation_coriolis calculates the large timestep tendency terms in the
3858	! u, v, and w momentum equations arise from the coriolis force. This version
3859	! subtracts off the horizontal velocities from the initial sounding when
3860	! computing the forcing terms, hence "perturbation" coriolis.
3861	!
3862	!</DESCRIPTION>
3863
3864	specified = .false.
3865	if(config_flags%specified .or. config_flags%nested) specified = .true.
3866
3867	ktf=MIN(kte,kde-1)
3868
3869	! coriolis for u-momentum equation
3870
3871	i_start = its
3872	i_end = ite
3873	IF ( config_flags%open_xs .or. specified .or. &
3874	config_flags%nested) i_start = MAX(ids+1,its)
3875	IF ( config_flags%open_xe .or. specified .or. &
3876	config_flags%nested) i_end = MIN(ide-1,ite)
3877	IF ( config_flags%periodic_x ) i_start = its
3878	IF ( config_flags%periodic_x ) i_end = ite
3879
3880	! compute perturbation mu*v for use in u momentum equation
3881
3882	DO j = jts, MIN(jte,jde-1)+1
3883	DO k=kts+1,ktf-1
3884	DO i = i_start-1, i_end
3885	z_at_v = 0.25*( phb(i,k,j )+phb(i,k+1,j ) &
3886	+phb(i,k,j-1)+phb(i,k+1,j-1) &
3887	+ph(i,k,j )+ph(i,k+1,j ) &
3888	+ph(i,k,j-1)+ph(i,k+1,j-1))/g
3889	wkp1 = min(1.,max(0.,z_at_v-z_base(k))/(z_base(k+1)-z_base(k)))
3890	wkm1 = min(1.,max(0.,z_base(k)-z_at_v)/(z_base(k)-z_base(k-1)))
3891	wk = 1.-wkp1-wkm1
3892	rv(i,k,j) = rv_in(i,k,j) - muv(i,j)*( &
3893	wkm1*v_base(k-1) &
3894	+wk *v_base(k ) &
3895	+wkp1*v_base(k+1) )
3896	ENDDO
3897	ENDDO
3898	ENDDO
3899
3900
3901	! pick up top and bottom v
3902
3903	DO j = jts, MIN(jte,jde-1)+1
3904	DO i = i_start-1, i_end
3905
3906	k = kts
3907	z_at_v = 0.25*( phb(i,k,j )+phb(i,k+1,j ) &
3908	+phb(i,k,j-1)+phb(i,k+1,j-1) &
3909	+ph(i,k,j )+ph(i,k+1,j ) &
3910	+ph(i,k,j-1)+ph(i,k+1,j-1))/g
3911	wkp1 = min(1.,max(0.,z_at_v-z_base(k))/(z_base(k+1)-z_base(k)))
3912	wk = 1.-wkp1
3913	rv(i,k,j) = rv_in(i,k,j) - muv(i,j)*( &
3914	+wk *v_base(k ) &
3915	+wkp1*v_base(k+1) )
3916
3917	k = ktf
3918	z_at_v = 0.25*( phb(i,k,j )+phb(i,k+1,j ) &
3919	+phb(i,k,j-1)+phb(i,k+1,j-1) &
3920	+ph(i,k,j )+ph(i,k+1,j ) &
3921	+ph(i,k,j-1)+ph(i,k+1,j-1))/g
3922	wkm1 = min(1.,max(0.,z_base(k)-z_at_v)/(z_base(k)-z_base(k-1)))
3923	wk = 1.-wkm1
3924	rv(i,k,j) = rv_in(i,k,j) - muv(i,j)*( &
3925	wkm1*v_base(k-1) &
3926	+wk *v_base(k ) )
3927
3928	ENDDO
3929	ENDDO
3930
3931	! compute coriolis forcing for u
3932
3933	! Map scale factors: see comments above for Coriolis
3934
3935	DO j = jts, MIN(jte,jde-1)
3936
3937	DO k=kts,ktf
3938	DO i = i_start, i_end
3939	ru_tend(i,k,j)=ru_tend(i,k,j) + (msfux(i,j)/msfuy(i,j))0.5(f(i,j)+f(i-1,j)) &
3940	0.25(rv(i-1,k,j+1)+rv(i,k,j+1)+rv(i-1,k,j)+rv(i,k,j)) &
3941	- 0.5(e(i,j)+e(i-1,j))0.5*(cosa(i,j)+cosa(i-1,j)) &
3942	0.25(rw(i-1,k+1,j)+rw(i-1,k,j)+rw(i,k+1,j)+rw(i,k,j))
3943	ENDDO
3944	ENDDO
3945
3946	! boundary loops for perturbation coriolis is needed for open bdy (20110225 JD)
3947	IF ( (config_flags%open_xs) .and. (its == ids) ) THEN
3948
3949	DO k=kts,ktf
3950
3951	ru_tend(its,k,j)=ru_tend(its,k,j) + (msfux(its,j)/msfuy(its,j))0.5(f(its,j)+f(its,j)) &
3952	0.25(rv(its,k,j+1)+rv(its,k,j+1)+rv(its,k,j)+rv(its,k,j)) &
3953	- 0.5(e(its,j)+e(its,j))0.5*(cosa(its,j)+cosa(its,j)) &
3954	0.25(rw(its,k+1,j)+rw(its,k,j)+rw(its,k+1,j)+rw(its,k,j))
3955
3956	ENDDO
3957
3958	ENDIF
3959
3960	IF ( (config_flags%open_xe) .and. (ite == ide) ) THEN
3961
3962	DO k=kts,ktf
3963
3964	ru_tend(ite,k,j)=ru_tend(ite,k,j) + (msfux(ite,j)/msfuy(ite,j))0.5(f(ite-1,j)+f(ite-1,j)) &
3965	0.25(rv(ite-1,k,j+1)+rv(ite-1,k,j+1)+rv(ite-1,k,j)+rv(ite-1,k,j)) &
3966	- 0.5(e(ite-1,j)+e(ite-1,j))0.5*(cosa(ite-1,j)+cosa(ite-1,j)) &
3967	0.25(rw(ite-1,k+1,j)+rw(ite-1,k,j)+rw(ite-1,k+1,j)+rw(ite-1,k,j))
3968
3969	ENDDO
3970
3971	ENDIF
3972
3973	ENDDO
3974
3975	! coriolis term for v-momentum equation
3976	! Map scale factors: see comments above for Coriolis
3977
3978	j_start = jts
3979	j_end = jte
3980
3981	IF ( config_flags%open_ys .or. specified .or. &
3982	config_flags%nested .or. config_flags%polar) j_start = MAX(jds+1,jts)
3983	IF ( config_flags%open_ye .or. specified .or. &
3984	config_flags%nested .or. config_flags%polar) j_end = MIN(jde-1,jte)
3985
3986	! compute perturbation mu*u for use in v momentum equation
3987
3988	DO j = j_start-1,j_end
3989	DO k=kts+1,ktf-1
3990	DO i = its, MIN(ite,ide-1)+1
3991	z_at_u = 0.25*( phb(i ,k,j)+phb(i ,k+1,j) &
3992	+phb(i-1,k,j)+phb(i-1,k+1,j) &
3993	+ph(i ,k,j)+ph(i ,k+1,j) &
3994	+ph(i-1,k,j)+ph(i-1,k+1,j))/g
3995	wkp1 = min(1.,max(0.,z_at_u-z_base(k))/(z_base(k+1)-z_base(k)))
3996	wkm1 = min(1.,max(0.,z_base(k)-z_at_u)/(z_base(k)-z_base(k-1)))
3997	wk = 1.-wkp1-wkm1
3998	ru(i,k,j) = ru_in(i,k,j) - muu(i,j)*( &
3999	wkm1*u_base(k-1) &
4000	+wk *u_base(k ) &
4001	+wkp1*u_base(k+1) )
4002	ENDDO
4003	ENDDO
4004	ENDDO
4005
4006	! pick up top and bottom u
4007
4008	DO j = j_start-1,j_end
4009	DO i = its, MIN(ite,ide-1)+1
4010
4011	k = kts
4012	z_at_u = 0.25*( phb(i ,k,j)+phb(i ,k+1,j) &
4013	+phb(i-1,k,j)+phb(i-1,k+1,j) &
4014	+ph(i ,k,j)+ph(i ,k+1,j) &
4015	+ph(i-1,k,j)+ph(i-1,k+1,j))/g
4016	wkp1 = min(1.,max(0.,z_at_u-z_base(k))/(z_base(k+1)-z_base(k)))
4017	wk = 1.-wkp1
4018	ru(i,k,j) = ru_in(i,k,j) - muu(i,j)*( &
4019	+wk *u_base(k ) &
4020	+wkp1*u_base(k+1) )
4021
4022
4023	k = ktf
4024	z_at_u = 0.25*( phb(i ,k,j)+phb(i ,k+1,j) &
4025	+phb(i-1,k,j)+phb(i-1,k+1,j) &
4026	+ph(i ,k,j)+ph(i ,k+1,j) &
4027	+ph(i-1,k,j)+ph(i-1,k+1,j))/g
4028	wkm1 = min(1.,max(0.,z_base(k)-z_at_u)/(z_base(k)-z_base(k-1)))
4029	wk = 1.-wkm1
4030	ru(i,k,j) = ru_in(i,k,j) - muu(i,j)*( &
4031	wkm1*u_base(k-1) &
4032	+wk *u_base(k ) )
4033
4034	ENDDO
4035	ENDDO
4036
4037	! compute coriolis forcing for v momentum equation
4038	! Map scale factors: see comments above for Coriolis
4039
4040	! boundary loops for perturbation coriolis is needed for open bdy (20110225 JD)
4041	IF ( (config_flags%open_ys) .and. (jts == jds) ) THEN
4042
4043	DO k=kts,ktf
4044	DO i=its,MIN(ide-1,ite)
4045
4046	rv_tend(i,k,jts)=rv_tend(i,k,jts) - (msfvy(i,jts)/msfvx(i,jts))0.5(f(i,jts)+f(i,jts)) &
4047	0.25(ru(i,k,jts)+ru(i+1,k,jts)+ru(i,k,jts)+ru(i+1,k,jts)) &
4048	+ (msfvy(i,jts)/msfvx(i,jts))0.5(e(i,jts)+e(i,jts))0.5(sina(i,jts)+sina(i,jts)) &
4049	0.25(rw(i,k+1,jts)+rw(i,k,jts)+rw(i,k+1,jts)+rw(i,k,jts))
4050
4051	ENDDO
4052	ENDDO
4053
4054	ENDIF
4055
4056	DO j=j_start, j_end
4057	DO k=kts,ktf
4058	DO i=its,MIN(ide-1,ite)
4059
4060	rv_tend(i,k,j)=rv_tend(i,k,j) - (msfvy(i,j)/msfvx(i,j))0.5(f(i,j)+f(i,j-1)) &
4061	0.25(ru(i,k,j)+ru(i+1,k,j)+ru(i,k,j-1)+ru(i+1,k,j-1)) &
4062	+ (msfvy(i,j)/msfvx(i,j))0.5(e(i,j)+e(i,j-1))0.5(sina(i,j)+sina(i,j-1)) &
4063	0.25(rw(i,k+1,j-1)+rw(i,k,j-1)+rw(i,k+1,j)+rw(i,k,j))
4064
4065	ENDDO
4066	ENDDO
4067	ENDDO
4068
4069
4070	! boundary loops for perturbation coriolis is needed for open bdy (20110225 JD)
4071	IF ( (config_flags%open_ye) .and. (jte == jde) ) THEN
4072
4073	DO k=kts,ktf
4074	DO i=its,MIN(ide-1,ite)
4075
4076	rv_tend(i,k,jte)=rv_tend(i,k,jte) - (msfvy(i,jte)/msfvx(i,jte))0.5(f(i,jte-1)+f(i,jte-1)) &
4077	0.25(ru(i,k,jte-1)+ru(i+1,k,jte-1)+ru(i,k,jte-1)+ru(i+1,k,jte-1)) &
4078	+ (msfvy(i,jte)/msfvx(i,jte))0.5(e(i,jte-1)+e(i,jte-1))0.5(sina(i,jte-1)+sina(i,jte-1)) &
4079	0.25(rw(i,k+1,jte-1)+rw(i,k,jte-1)+rw(i,k+1,jte-1)+rw(i,k,jte-1))
4080
4081	ENDDO
4082	ENDDO
4083
4084	ENDIF
4085
4086	! coriolis term for w-mometum
4087	! Map scale factors: see comments above for Coriolis
4088
4089	DO j=jts,MIN(jte, jde-1)
4090	DO k=kts+1,ktf
4091	DO i=its,MIN(ite, ide-1)
4092
4093	rw_tend(i,k,j)=rw_tend(i,k,j) + e(i,j)* &
4094	(cosa(i,j)0.5(fzm(k)*(ru(i,k,j)+ru(i+1,k,j)) &
4095	+fzp(k)*(ru(i,k-1,j)+ru(i+1,k-1,j))) &
4096	-(msftx(i,j)/msfty(i,j))sina(i,j)0.5(fzm(k)(rv(i,k,j)+rv(i,k,j+1)) &
4097	+fzp(k)*(rv(i,k-1,j)+rv(i,k-1,j+1))))
4098
4099	ENDDO
4100	ENDDO
4101	ENDDO
4102
4103	END SUBROUTINE perturbation_coriolis
4104
4105	!------------------------------------------------------------------------------
4106
4107	SUBROUTINE curvature ( ru, rv, rw, u, v, w, ru_tend, rv_tend, rw_tend, &
4108	config_flags, &
4109	msfux, msfuy, msfvx, msfvy, msftx, msfty, &
4110	xlat, fzm, fzp, rdx, rdy, &
4111	ids, ide, jds, jde, kds, kde, &
4112	ims, ime, jms, jme, kms, kme, &
4113	its, ite, jts, jte, kts, kte )
4114
4115
4116	IMPLICIT NONE
4117
4118	! Input data
4119
4120	TYPE(grid_config_rec_type) , INTENT(IN ) :: config_flags
4121
4122	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
4123	ims, ime, jms, jme, kms, kme, &
4124	its, ite, jts, jte, kts, kte
4125
4126	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , &
4127	INTENT(INOUT) :: ru_tend, &
4128	rv_tend, &
4129	rw_tend
4130
4131	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , &
4132	INTENT(IN ) :: ru, &
4133	rv, &
4134	rw, &
4135	u, &
4136	v, &
4137	w
4138
4139	REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: msfux, &
4140	msfuy, &
4141	msfvx, &
4142	msfvy, &
4143	msftx, &
4144	msfty, &
4145	xlat
4146
4147	REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fzm, &
4148	fzp
4149
4150	REAL , INTENT(IN ) :: rdx, &
4151	rdy
4152
4153	! Local data
4154
4155	! INTEGER :: i, j, k, itf, jtf, ktf, kp1, im, ip, jm, jp
4156	INTEGER :: i, j, k, itf, jtf, ktf
4157	INTEGER :: i_start, i_end, j_start, j_end
4158	! INTEGER :: irmin, irmax, jrmin, jrmax
4159
4160	REAL , DIMENSION( its-1:ite , kts:kte, jts-1:jte ) :: vxgm
4161
4162	LOGICAL :: specified
4163
4164	!<DESCRIPTION>
4165	!
4166	! curvature calculates the large timestep tendency terms in the
4167	! u, v, and w momentum equations arise from the curvature terms.
4168	!
4169	!</DESCRIPTION>
4170
4171	specified = .false.
4172	if(config_flags%specified .or. config_flags%nested) specified = .true.
4173
4174	itf=MIN(ite,ide-1)
4175	jtf=MIN(jte,jde-1)
4176	ktf=MIN(kte,kde-1)
4177
4178	! irmin = ims
4179	! irmax = ime
4180	! jrmin = jms
4181	! jrmax = jme
4182	! IF ( config_flags%open_xs ) irmin = ids
4183	! IF ( config_flags%open_xe ) irmax = ide-1
4184	! IF ( config_flags%open_ys ) jrmin = jds
4185	! IF ( config_flags%open_ye ) jrmax = jde-1
4186
4187	! Define v cross grad m at scalar points - vxgm(i,j)
4188
4189	i_start = its-1
4190	i_end = ite
4191	j_start = jts-1
4192	j_end = jte
4193
4194	IF ( ( config_flags%open_xs .or. specified .or. &
4195	config_flags%nested) .and. (its == ids) ) i_start = its
4196	IF ( ( config_flags%open_xe .or. specified .or. &
4197	config_flags%nested) .and. (ite == ide) ) i_end = ite-1
4198	IF ( ( config_flags%open_ys .or. specified .or. &
4199	config_flags%nested .or. config_flags%polar) .and. (jts == jds) ) j_start = jts
4200	IF ( ( config_flags%open_ye .or. specified .or. &
4201	config_flags%nested .or. config_flags%polar) .and. (jte == jde) ) j_end = jte-1
4202	IF ( config_flags%periodic_x ) i_start = its-1
4203	IF ( config_flags%periodic_x ) i_end = ite
4204
4205	DO j=j_start, j_end
4206	DO k=kts,ktf
4207	DO i=i_start, i_end
4208	! Map scale factor notes:
4209	! msf...y is constant everywhere for cylindrical map projection
4210	! msf...x varies with y only
4211	! But we know that this is not = 0 for cylindrical,
4212	! therefore use msfvX in 1st line
4213	! which => by symmetry use msfuY in 2nd line - ???
4214	vxgm(i,k,j)=0.5(u(i,k,j)+u(i+1,k,j))(msfvx(i,j+1)-msfvx(i,j))*rdy - &
4215	0.5(v(i,k,j)+v(i,k,j+1))(msfuy(i+1,j)-msfuy(i,j))*rdx
4216	ENDDO
4217	ENDDO
4218	ENDDO
4219
4220	! Pick up the boundary rows for open (radiation) lateral b.c.
4221	! Rather crude at present, we are assuming there is no
4222	! variation in this term at the boundary.
4223
4224	IF ( ( config_flags%open_xs .or. (specified .AND. .NOT. config_flags%periodic_x) .or. &
4225	config_flags%nested) .and. (its == ids) ) THEN
4226
4227	DO j = jts, jte-1
4228	DO k = kts, ktf
4229	vxgm(its-1,k,j) = vxgm(its,k,j)
4230	ENDDO
4231	ENDDO
4232
4233	ENDIF
4234
4235	IF ( ( config_flags%open_xe .or. (specified .AND. .NOT. config_flags%periodic_x) .or. &
4236	config_flags%nested) .and. (ite == ide) ) THEN
4237
4238	DO j = jts, jte-1
4239	DO k = kts, ktf
4240	vxgm(ite,k,j) = vxgm(ite-1,k,j)
4241	ENDDO
4242	ENDDO
4243
4244	ENDIF
4245
4246	! Polar boundary condition:
4247	! The following change is needed in case one tries using the vxgm route with
4248	! polar B.C.'s in the future, but not needed if 'tan' used
4249	IF ( ( config_flags%open_ys .or. specified .or. &
4250	config_flags%nested .or. config_flags%polar) .and. (jts == jds) ) THEN
4251
4252	DO k = kts, ktf
4253	DO i = its-1, ite
4254	vxgm(i,k,jts-1) = vxgm(i,k,jts)
4255	ENDDO
4256	ENDDO
4257
4258	ENDIF
4259
4260	! Polar boundary condition:
4261	! The following change is needed in case one tries using the vxgm route with
4262	! polar B.C.'s in the future, but not needed if 'tan' used
4263	IF ( ( config_flags%open_ye .or. specified .or. &
4264	config_flags%nested .or. config_flags%polar) .and. (jte == jde) ) THEN
4265
4266	DO k = kts, ktf
4267	DO i = its-1, ite
4268	vxgm(i,k,jte) = vxgm(i,k,jte-1)
4269	ENDDO
4270	ENDDO
4271
4272	ENDIF
4273
4274	! curvature term for u momentum eqn.
4275
4276	! Map scale factor notes:
4277	! ADT eqn 44, RHS terms 4 and 5, in cylindrical: mu u v tan(lat)/(a my)
4278	! - mu u w /(a my)
4279	! ru = mu u / my ; rw = mu w / my ; rv = mu v / mx
4280	! => terms are:
4281	! (mx/my)u rv tan(lat) / a - u rw / a = (u/a)[(mx/my) rv tan(lat) - rw]
4282	! ru v tan(lat) / a - u rw / a
4283	! xlat defined with end points half grid space from pole,
4284	! hence are on u latitude points
4285
4286	i_start = its
4287	IF ( config_flags%open_xs .or. specified .or. &
4288	config_flags%nested) i_start = MAX ( ids+1 , its )
4289	IF ( config_flags%open_xe .or. specified .or. &
4290	config_flags%nested) i_end = MIN ( ide-1 , ite )
4291	IF ( config_flags%periodic_x ) i_start = its
4292	IF ( config_flags%periodic_x ) i_end = ite
4293
4294	! Polar boundary condition
4295	IF ((config_flags%map_proj == 6) .OR. (config_flags%polar)) THEN
4296
4297	DO j=jts,MIN(jde-1,jte)
4298	DO k=kts,ktf
4299	DO i=i_start,i_end
4300
4301	ru_tend(i,k,j)=ru_tend(i,k,j) + u(i,k,j)reradius ( &
4302	(msfux(i,j)/msfuy(i,j))0.25(rv(i-1,k,j+1)+rv(i,k,j+1)+ &
4303	rv(i-1,k,j)+rv(i,k,j))tan(xlat(i,j)degrad) &
4304	- 0.25*(rw(i-1,k+1,j)+rw(i-1,k,j)+rw(i,k+1,j)+rw(i,k,j)) )
4305	ENDDO
4306	ENDDO
4307	ENDDO
4308
4309	ELSE ! normal code
4310
4311
4312	DO j=jts,MIN(jde-1,jte)
4313	DO k=kts,ktf
4314	DO i=i_start,i_end
4315
4316	ru_tend(i,k,j)=ru_tend(i,k,j) + 0.5*(vxgm(i,k,j)+vxgm(i-1,k,j)) &
4317	0.25(rv(i-1,k,j+1)+rv(i,k,j+1)+rv(i-1,k,j)+rv(i,k,j)) &
4318	- u(i,k,j)*reradius &
4319	0.25(rw(i-1,k+1,j)+rw(i-1,k,j)+rw(i,k+1,j)+rw(i,k,j))
4320
4321	ENDDO
4322	ENDDO
4323	ENDDO
4324
4325	END IF
4326
4327	! curvature term for v momentum eqn.
4328
4329	! Map scale factor notes
4330	! ADT eqn 45, RHS terms 4 and 5, in cylindrical: - mu u*u tan(lat)/(a mx)
4331	! - mu v w /(a mx)
4332	! ru = mu u / my ; rw = mu w / my ; rv = mu v / mx
4333	! terms are:
4334	! - (my/mx)u ru tan(lat) / a - (my/mx)v rw / a
4335	! = - [my/(mxa)][u ru tan(lat) + v rw]
4336	! - (1/a)[(my/mx)u ru tan(lat) + w rv]
4337	! xlat defined with end points half grid space from pole, hence are on
4338	! u latitude points => av here
4339	!
4340	! in original wrf, there was a sign error for the rw contribution
4341
4342	j_start = jts
4343	IF ( config_flags%open_ys .or. specified .or. &
4344	config_flags%nested .or. config_flags%polar) j_start = MAX ( jds+1 , jts )
4345	IF ( config_flags%open_ye .or. specified .or. &
4346	config_flags%nested .or. config_flags%polar) j_end = MIN ( jde-1 , jte )
4347
4348	IF ((config_flags%map_proj == 6) .OR. (config_flags%polar)) THEN
4349
4350	DO j=j_start,j_end
4351	DO k=kts,ktf
4352	DO i=its,MIN(ite,ide-1)
4353	rv_tend(i,k,j)=rv_tend(i,k,j) - (msfvy(i,j)/msfvx(i,j))reradius ( &
4354	0.25(u(i,k,j)+u(i+1,k,j)+u(i,k,j-1)+u(i+1,k,j-1)) &
4355	tan((xlat(i,j)+xlat(i,j-1))0.5degrad)* &
4356	0.25*(ru(i,k,j)+ru(i+1,k,j)+ru(i,k,j-1)+ru(i+1,k,j-1)) &
4357	+ v(i,k,j)0.25(rw(i,k+1,j-1)+rw(i,k,j-1)+ &
4358	rw(i,k+1,j)+rw(i,k,j)) )
4359	ENDDO
4360	ENDDO
4361	ENDDO
4362
4363	ELSE ! normal code
4364
4365	DO j=j_start,j_end
4366	DO k=kts,ktf
4367	DO i=its,MIN(ite,ide-1)
4368
4369	rv_tend(i,k,j)=rv_tend(i,k,j) - 0.5*(vxgm(i,k,j)+vxgm(i,k,j-1)) &
4370	0.25(ru(i,k,j)+ru(i+1,k,j)+ru(i,k,j-1)+ru(i+1,k,j-1)) &
4371	- (msfvy(i,j)/msfvx(i,j))v(i,k,j)reradius &
4372	0.25(rw(i,k+1,j-1)+rw(i,k,j-1)+rw(i,k+1,j)+rw(i,k,j))
4373
4374	ENDDO
4375	ENDDO
4376	ENDDO
4377
4378	END IF
4379
4380	! curvature term for vertical momentum eqn.
4381
4382	! Notes on map scale factors:
4383	! ADT eqn 46, RHS term 4: [mu/(a my)][uu + v*v]
4384	! ru = mu u / my ; rw = mu w / my ; rv = mu v / mx
4385	! terms are: u ru / a + (mx/my)v rv / a
4386
4387	DO j=jts,MIN(jte,jde-1)
4388	DO k=MAX(2,kts),ktf
4389	DO i=its,MIN(ite,ide-1)
4390
4391	rw_tend(i,k,j)=rw_tend(i,k,j) + reradius* &
4392	(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))) &
4393	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))) &
4394	+(msftx(i,j)/msfty(i,j))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))) &
4395	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))))
4396
4397	ENDDO
4398	ENDDO
4399	ENDDO
4400
4401	END SUBROUTINE curvature
4402
4403	!------------------------------------------------------------------------------
4404
4405	SUBROUTINE decouple ( rr, rfield, field, name, config_flags, &
4406	fzm, fzp, &
4407	ids, ide, jds, jde, kds, kde, &
4408	ims, ime, jms, jme, kms, kme, &
4409	its, ite, jts, jte, kts, kte )
4410
4411	IMPLICIT NONE
4412
4413	! Input data
4414
4415	TYPE(grid_config_rec_type) , INTENT(IN ) :: config_flags
4416
4417	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
4418	ims, ime, jms, jme, kms, kme, &
4419	its, ite, jts, jte, kts, kte
4420
4421	CHARACTER(LEN=1) , INTENT(IN ) :: name
4422
4423	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: rfield
4424
4425	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: rr
4426
4427	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT ) :: field
4428
4429	REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fzm, fzp
4430
4431	! Local data
4432
4433	INTEGER :: i, j, k, itf, jtf, ktf
4434
4435	!<DESCRIPTION>
4436	!
4437	! decouple decouples a variable from the column dry-air mass.
4438	!
4439	!</DESCRIPTION>
4440
4441	ktf=MIN(kte,kde-1)
4442
4443	IF (name .EQ. 'u')THEN
4444	itf=ite
4445	jtf=MIN(jte,jde-1)
4446
4447	DO j=jts,jtf
4448	DO k=kts,ktf
4449	DO i=its,itf
4450	field(i,k,j)=rfield(i,k,j)/(0.5*(rr(i,k,j)+rr(i-1,k,j)))
4451	ENDDO
4452	ENDDO
4453	ENDDO
4454
4455	ELSE IF (name .EQ. 'v')THEN
4456	itf=MIN(ite,ide-1)
4457	jtf=jte
4458
4459	DO j=jts,jtf
4460	DO k=kts,ktf
4461	DO i=its,itf
4462	field(i,k,j)=rfield(i,k,j)/(0.5*(rr(i,k,j)+rr(i,k,j-1)))
4463	ENDDO
4464	ENDDO
4465	ENDDO
4466
4467	ELSE IF (name .EQ. 'w')THEN
4468	itf=MIN(ite,ide-1)
4469	jtf=MIN(jte,jde-1)
4470	DO j=jts,jtf
4471	DO k=kts+1,ktf
4472	DO i=its,itf
4473	field(i,k,j)=rfield(i,k,j)/(fzm(k)rr(i,k,j)+fzp(k)rr(i,k-1,j))
4474	ENDDO
4475	ENDDO
4476	ENDDO
4477
4478	DO j=jts,jtf
4479	DO i=its,itf
4480	field(i,kte,j) = 0.
4481	ENDDO
4482	ENDDO
4483
4484	ELSE
4485	itf=MIN(ite,ide-1)
4486	jtf=MIN(jte,jde-1)
4487	! For theta we will decouple tb and tp and add them to give t afterwards
4488	DO j=jts,jtf
4489	DO k=kts,ktf
4490	DO i=its,itf
4491	field(i,k,j)=rfield(i,k,j)/rr(i,k,j)
4492	ENDDO
4493	ENDDO
4494	ENDDO
4495
4496	ENDIF
4497
4498	END SUBROUTINE decouple
4499
4500	!-------------------------------------------------------------------------------
4501
4502
4503	SUBROUTINE zero_tend ( tendency, &
4504	ids, ide, jds, jde, kds, kde, &
4505	ims, ime, jms, jme, kms, kme, &
4506	its, ite, jts, jte, kts, kte )
4507
4508
4509	IMPLICIT NONE
4510
4511	! Input data
4512
4513	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
4514	ims, ime, jms, jme, kms, kme, &
4515	its, ite, jts, jte, kts, kte
4516
4517	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency
4518
4519	! Local data
4520
4521	INTEGER :: i, j, k, itf, jtf, ktf
4522
4523	!<DESCRIPTION>
4524	!
4525	! zero_tend sets the input tendency array to zero.
4526	!
4527	!</DESCRIPTION>
4528
4529	DO j = jts, jte
4530	DO k = kts, kte
4531	DO i = its, ite
4532	tendency(i,k,j) = 0.
4533	ENDDO
4534	ENDDO
4535	ENDDO
4536
4537	END SUBROUTINE zero_tend
4538
4539	!-------------------------------------------------------------------------------
4540	! Sets the an array on the polar v point(s) to zero
4541	SUBROUTINE zero_pole ( field, &
4542	ids, ide, jds, jde, kds, kde, &
4543	ims, ime, jms, jme, kms, kme, &
4544	its, ite, jts, jte, kts, kte )
4545
4546
4547	IMPLICIT NONE
4548
4549	! Input data
4550
4551	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
4552	ims, ime, jms, jme, kms, kme, &
4553	its, ite, jts, jte, kts, kte
4554
4555	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: field
4556
4557	! Local data
4558
4559	INTEGER :: i, k
4560
4561	IF (jts == jds) THEN
4562	DO k = kts, kte
4563	DO i = its-1, ite+1
4564	field(i,k,jts) = 0.
4565	END DO
4566	END DO
4567	END IF
4568	IF (jte == jde) THEN
4569	DO k = kts, kte
4570	DO i = its-1, ite+1
4571	field(i,k,jte) = 0.
4572	END DO
4573	END DO
4574	END IF
4575
4576	END SUBROUTINE zero_pole
4577
4578	!-------------------------------------------------------------------------------
4579	! Sets the an array on the polar v point(s)
4580	SUBROUTINE pole_point_bc ( field, &
4581	ids, ide, jds, jde, kds, kde, &
4582	ims, ime, jms, jme, kms, kme, &
4583	its, ite, jts, jte, kts, kte )
4584
4585
4586	IMPLICIT NONE
4587
4588	! Input data
4589
4590	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
4591	ims, ime, jms, jme, kms, kme, &
4592	its, ite, jts, jte, kts, kte
4593
4594	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: field
4595
4596	! Local data
4597
4598	INTEGER :: i, k
4599
4600	IF (jts == jds) THEN
4601	DO k = kts, kte
4602	DO i = its, ite
4603	! field(i,k,jts) = 2*field(i,k,jts+1) - field(i,k,jts+2)
4604	field(i,k,jts) = field(i,k,jts+1)
4605	END DO
4606	END DO
4607	END IF
4608	IF (jte == jde) THEN
4609	DO k = kts, kte
4610	DO i = its, ite
4611	! field(i,k,jte) = 2*field(i,k,jte-1) - field(i,k,jte-2)
4612	field(i,k,jte) = field(i,k,jte-1)
4613	END DO
4614	END DO
4615	END IF
4616
4617	END SUBROUTINE pole_point_bc
4618
4619	!======================================================================
4620	! physics prep routines
4621	!======================================================================
4622
4623	SUBROUTINE phy_prep ( config_flags, & ! input
4624	mu, muu, muv, u, v, p, pb, alt, ph, & ! input
4625	phb, t, tsk, moist, n_moist, & ! input
4626	rho, th_phy, p_phy , pi_phy , & ! output
4627	u_phy, v_phy, p8w, t_phy, t8w, & ! output
4628	z, z_at_w, dz8w, & ! output
4629	p_hyd, p_hyd_w, & ! output
4630	fzm, fzp, znw, p_top, & ! params
4631	RTHRATEN, &
4632	RTHBLTEN, RUBLTEN, RVBLTEN, &
4633	RQVBLTEN, RQCBLTEN, RQIBLTEN, &
4634	RUCUTEN, RVCUTEN, RTHCUTEN, &
4635	RQVCUTEN, RQCCUTEN, RQRCUTEN, &
4636	RQICUTEN, RQSCUTEN, &
4637	RUSHTEN, RVSHTEN, RTHSHTEN, &
4638	RQVSHTEN, RQCSHTEN, RQRSHTEN, &
4639	RQISHTEN, RQSSHTEN, RQGSHTEN, &
4640	RTHFTEN, RQVFTEN, &
4641	RUNDGDTEN, RVNDGDTEN, RTHNDGDTEN, &
4642	RPHNDGDTEN,RQVNDGDTEN, RMUNDGDTEN, &
4643	ids, ide, jds, jde, kds, kde, &
4644	ims, ime, jms, jme, kms, kme, &
4645	its, ite, jts, jte, kts, kte )
4646	!----------------------------------------------------------------------
4647	IMPLICIT NONE
4648	!----------------------------------------------------------------------
4649
4650	TYPE(grid_config_rec_type) , INTENT(IN ) :: config_flags
4651
4652	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
4653	ims, ime, jms, jme, kms, kme, &
4654	its, ite, jts, jte, kts, kte
4655	INTEGER , INTENT(IN ) :: n_moist
4656
4657	REAL, DIMENSION( ims:ime, kms:kme , jms:jme , n_moist ), INTENT(IN) :: moist
4658
4659
4660	REAL , DIMENSION( ims:ime, jms:jme ), INTENT(IN ) :: TSK, mu, muu, muv
4661
4662	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , &
4663	INTENT( OUT) :: u_phy, &
4664	v_phy, &
4665	pi_phy, &
4666	p_phy, &
4667	p8w, &
4668	t_phy, &
4669	th_phy, &
4670	t8w, &
4671	rho, &
4672	z, &
4673	dz8w, &
4674	z_at_w
4675
4676	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , &
4677	INTENT( OUT) :: p_hyd, &
4678	p_hyd_w
4679
4680	REAL , INTENT(IN ) :: p_top
4681
4682	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , &
4683	INTENT(IN ) :: pb, &
4684	p, &
4685	u, &
4686	v, &
4687	alt, &
4688	ph, &
4689	phb, &
4690	t
4691
4692
4693	REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fzm, &
4694	fzp
4695
4696	REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: znw
4697
4698	REAL, DIMENSION( ims:ime , kms:kme, jms:jme ), &
4699	INTENT(INOUT) :: RTHRATEN
4700
4701	REAL, DIMENSION( ims:ime , kms:kme, jms:jme ), &
4702	INTENT(INOUT) :: RUCUTEN, &
4703	RVCUTEN, &
4704	RTHCUTEN, &
4705	RQVCUTEN, &
4706	RQCCUTEN, &
4707	RQRCUTEN, &
4708	RQICUTEN, &
4709	RQSCUTEN, &
4710	RUSHTEN, &
4711	RVSHTEN, &
4712	RTHSHTEN, &
4713	RQVSHTEN, &
4714	RQCSHTEN, &
4715	RQRSHTEN, &
4716	RQISHTEN, &
4717	RQSSHTEN, &
4718	RQGSHTEN
4719
4720	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) , &
4721	INTENT(INOUT) :: RUBLTEN, &
4722	RVBLTEN, &
4723	RTHBLTEN, &
4724	RQVBLTEN, &
4725	RQCBLTEN, &
4726	RQIBLTEN
4727
4728	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) , &
4729	INTENT(INOUT) :: RTHFTEN, &
4730	RQVFTEN
4731
4732	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) , &
4733	INTENT(INOUT) :: RUNDGDTEN, &
4734	RVNDGDTEN, &
4735	RTHNDGDTEN, &
4736	RPHNDGDTEN, &
4737	RQVNDGDTEN
4738
4739	REAL, DIMENSION( ims:ime, jms:jme ) , &
4740	INTENT(INOUT) :: RMUNDGDTEN
4741
4742	INTEGER :: i_start, i_end, j_start, j_end, k_start, k_end, i_startu, j_startv
4743	INTEGER :: i, j, k
4744	REAL :: w1, w2, z0, z1, z2
4745	REAL :: e_vapor
4746
4747	#ifdef LMDZ1
4748	CHARACTER*256 :: message
4749	INTEGER :: im2, jm2, km2
4750	INTEGER :: ix,iy,iz
4751	CHARACTER(LEN=50) :: errmsg
4752
4753	errmsg = 'ERROR -- error -- ERROR -- error'
4754
4755	im2 = config_flags%i_check_point
4756	jm2 = config_flags%j_check_point
4757	km2 = config_flags%k_check_point
4758	#endif
4759
4760	#ifdef LMDZ1
4761	WRITE(message, *)' phy_prep: inside'
4762	CALL wrf_debug(200, message)
4763	WRITE(message,*)' psfc_tend: 0.00000000 p sfc: ', &
4764	p8w(im2,kms,jm2)
4765	CALL wrf_debug(200, message)
4766	#endif
4767
4768	!-----------------------------------------------------------------------
4769
4770	!<DESCRIPTION>
4771	!
4772	! phys_prep calculates a number of diagnostic quantities needed by
4773	! the physics routines. It also decouples the physics tendencies from
4774	! the column dry-air mass (the physics routines expect to see/update the
4775	! uncoupled tendencies).
4776	!
4777	!</DESCRIPTION>
4778
4779	! set up loop bounds for this grid's boundary conditions
4780
4781	i_start = its
4782	i_end = min( ite,ide-1 )
4783	j_start = jts
4784	j_end = min( jte,jde-1 )
4785
4786	k_start = kts
4787	k_end = min( kte, kde-1 )
4788
4789	! compute thermodynamics and velocities at pressure points (or half levels)
4790
4791	do j = j_start,j_end
4792	do k = k_start, k_end
4793	do i = i_start, i_end
4794
4795	th_phy(i,k,j) = t(i,k,j) + t0
4796	p_phy(i,k,j) = p(i,k,j) + pb(i,k,j)
4797	pi_phy(i,k,j) = (p_phy(i,k,j)/p1000mb)**rcp
4798	t_phy(i,k,j) = th_phy(i,k,j)*pi_phy(i,k,j)
4799	rho(i,k,j) = 1./alt(i,k,j)*(1.+moist(i,k,j,P_QV))
4800	u_phy(i,k,j) = 0.5*(u(i,k,j)+u(i+1,k,j))
4801	v_phy(i,k,j) = 0.5*(v(i,k,j)+v(i,k,j+1))
4802
4803	enddo
4804	enddo
4805	enddo
4806
4807	! compute z at w points
4808
4809	do j = j_start,j_end
4810	do k = k_start, kte
4811	do i = i_start, i_end
4812	z_at_w(i,k,j) = (phb(i,k,j)+ph(i,k,j))/g
4813	enddo
4814	enddo
4815	enddo
4816
4817	do j = j_start,j_end
4818	do k = k_start, kte-1
4819	do i = i_start, i_end
4820	dz8w(i,k,j) = z_at_w(i,k+1,j)-z_at_w(i,k,j)
4821	enddo
4822	enddo
4823	enddo
4824
4825	do j = j_start,j_end
4826	do i = i_start, i_end
4827	dz8w(i,kte,j) = 0.
4828	enddo
4829	enddo
4830
4831	! compute z at p points or half levels (average of z at full levels)
4832
4833	do j = j_start,j_end
4834	do k = k_start, k_end
4835	do i = i_start, i_end
4836	z(i,k,j) = 0.5*(z_at_w(i,k,j) +z_at_w(i,k+1,j) )
4837	enddo
4838	enddo
4839	enddo
4840
4841	! interp t and p to full levels
4842
4843	do j = j_start,j_end
4844	do k = 2, k_end
4845	do i = i_start, i_end
4846	p8w(i,k,j) = fzm(k)p_phy(i,k,j)+fzp(k)p_phy(i,k-1,j)
4847	t8w(i,k,j) = fzm(k)t_phy(i,k,j)+fzp(k)t_phy(i,k-1,j)
4848	enddo
4849	enddo
4850	enddo
4851
4852	#ifdef LMDZ1
4853	WRITE(message, *)' phy_prep: z-levels'
4854	CALL wrf_debug(200, message)
4855	WRITE(message,*)' psfc_tend: 0.00000000 p sfc: ', &
4856	p8w(im2,kms,jm2)
4857	CALL wrf_debug(200, message)
4858	#endif
4859
4860	! extrapolate p and t to surface and top.
4861	! we'll use an extrapolation in z for now
4862
4863	do j = j_start,j_end
4864	do i = i_start, i_end
4865
4866	! bottom
4867
4868	z0 = z_at_w(i,1,j)
4869	z1 = z(i,1,j)
4870	z2 = z(i,2,j)
4871	w1 = (z0 - z2)/(z1 - z2)
4872	w2 = 1. - w1
4873	p8w(i,1,j) = w1p_phy(i,1,j)+w2p_phy(i,2,j)
4874	t8w(i,1,j) = w1t_phy(i,1,j)+w2t_phy(i,2,j)
4875
4876	! top
4877
4878	z0 = z_at_w(i,kte,j)
4879	z1 = z(i,k_end,j)
4880	z2 = z(i,k_end-1,j)
4881	w1 = (z0 - z2)/(z1 - z2)
4882	w2 = 1. - w1
4883
4884	! p8w(i,kde,j) = w1p_phy(i,kde-1,j)+w2p_phy(i,kde-2,j)
4885	!!! bug fix extrapolate ln(p) so p is positive definite
4886	p8w(i,kde,j) = exp(w1log(p_phy(i,kde-1,j))+w2log(p_phy(i,kde-2,j)))
4887	t8w(i,kde,j) = w1t_phy(i,kde-1,j)+w2t_phy(i,kde-2,j)
4888
4889	enddo
4890	enddo
4891
4892	#ifdef LMDZ1
4893	WRITE(message, *)' phy_prep: after p-top/bottom'
4894	CALL wrf_debug(200, message)
4895	WRITE(message,*)' psfc_tend: 0.00000000 p sfc: ', &
4896	p8w(im2,kde,jm2)
4897	CALL wrf_debug(200, message)
4898	z0 = z_at_w(im2,1,jm2)
4899	z1 = z(im2,1,jm2)
4900	z2 = z(im2,2,jm2)
4901	w1 = (z0 - z2)/(z1 - z2)
4902	w2 = 1. - w1
4903	WRITE(message,*)' z0: ',z0,' z1: ',z1,' z2: ',z2,' w1: ',w1,' w2: ',w2
4904	CALL wrf_debug(200, message)
4905	WRITE(message,*)' phb: ', phb(im2,1,jm2),' ph: ',ph(im2,1,jm2),' phb 2: ', &
4906	phb(im2,2,jm2),' ph 2: ',ph(im2,2,jm2),' pb: ', pb(im2,1,jm2),' p: ',p(im2,1,jm2),' pb 2: ', &
4907	pb(im2,2,jm2),' p 2: ',p(im2,2,jm2)
4908	CALL wrf_debug(200, message)
4909	WRITE(message,*)' p_phy: ',p_phy(im2,1,jm2),' p_phy 2: ',p_phy(im2,2,jm2)
4910	CALL wrf_debug(200, message)
4911	#endif
4912
4913	! calculate hydrostatic pressure at both full and half levels
4914	! first, full level p: assuming dry over model top
4915
4916	do j = j_start,j_end
4917	do i = i_start, i_end
4918	p_hyd_w(i,kte,j) = p_top
4919	enddo
4920	enddo
4921
4922	e_vapor = 0.
4923	do j = j_start,j_end
4924	do k = kte-1, k_start, -1
4925	do i = i_start, i_end
4926	! e_vapor = 1./alt(i,k,j)moist(i,k,j,P_QV)g*dz8w(i,k,j)
4927	! p_hyd_w(i,k,j) = p_hyd_w(i,k+1,j)+e_vapor
4928	p_hyd_w(i,k,j) = p_hyd_w(i,k+1,j)+1./alt(i,k,j)(1.+moist(i,k,j,P_QV))g*dz8w(i,k,j)
4929	enddo
4930	enddo
4931	enddo
4932
4933	! now calculate hydrostatic pressure at half levels
4934
4935	do j = j_start,j_end
4936	do k = k_start, k_end
4937	do i = i_start, i_end
4938	p_hyd(i,k,j) = 0.5*(p_hyd_w(i,k,j)+p_hyd_w(i,k+1,j))
4939	enddo
4940	enddo
4941	enddo
4942
4943	! decouple all physics tendencies
4944
4945	IF (config_flags%ra_lw_physics .gt. 0 .or. config_flags%ra_sw_physics .gt. 0) THEN
4946
4947	DO J=j_start,j_end
4948	DO K=k_start,k_end
4949	DO I=i_start,i_end
4950	RTHRATEN(I,K,J)=RTHRATEN(I,K,J)/mu(I,J)
4951	ENDDO
4952	ENDDO
4953	ENDDO
4954
4955	ENDIF
4956
4957	IF (config_flags%cu_physics .gt. 0) THEN
4958
4959	DO J=j_start,j_end
4960	DO I=i_start,i_end
4961	DO K=k_start,k_end
4962	RUCUTEN(I,K,J) =RUCUTEN(I,K,J)/mu(I,J)
4963	RVCUTEN(I,K,J) =RVCUTEN(I,K,J)/mu(I,J)
4964	RTHCUTEN(I,K,J)=RTHCUTEN(I,K,J)/mu(I,J)
4965	ENDDO
4966	ENDDO
4967	ENDDO
4968
4969	IF (P_QV .ge. PARAM_FIRST_SCALAR)THEN
4970	DO J=j_start,j_end
4971	DO I=i_start,i_end
4972	DO K=k_start,k_end
4973	RQVCUTEN(I,K,J)=RQVCUTEN(I,K,J)/mu(I,J)
4974	ENDDO
4975	ENDDO
4976	ENDDO
4977	ENDIF
4978
4979	IF (P_QC .ge. PARAM_FIRST_SCALAR)THEN
4980	DO J=j_start,j_end
4981	DO I=i_start,i_end
4982	DO K=k_start,k_end
4983	RQCCUTEN(I,K,J)=RQCCUTEN(I,K,J)/mu(I,J)
4984	ENDDO
4985	ENDDO
4986	ENDDO
4987	ENDIF
4988
4989	IF (P_QR .ge. PARAM_FIRST_SCALAR)THEN
4990	DO J=j_start,j_end
4991	DO I=i_start,i_end
4992	DO K=k_start,k_end
4993	RQRCUTEN(I,K,J)=RQRCUTEN(I,K,J)/mu(I,J)
4994	ENDDO
4995	ENDDO
4996	ENDDO
4997	ENDIF
4998
4999	IF (P_QI .ge. PARAM_FIRST_SCALAR)THEN
5000	DO J=j_start,j_end
5001	DO I=i_start,i_end
5002	DO K=k_start,k_end
5003	RQICUTEN(I,K,J)=RQICUTEN(I,K,J)/mu(I,J)
5004	ENDDO
5005	ENDDO
5006	ENDDO
5007	ENDIF
5008
5009	IF(P_QS .ge. PARAM_FIRST_SCALAR)THEN
5010	DO J=j_start,j_end
5011	DO I=i_start,i_end
5012	DO K=k_start,k_end
5013	RQSCUTEN(I,K,J)=RQSCUTEN(I,K,J)/mu(I,J)
5014	ENDDO
5015	ENDDO
5016	ENDDO
5017	ENDIF
5018
5019	ENDIF
5020
5021	IF (config_flags%shcu_physics .gt. 0) THEN
5022
5023	DO J=j_start,j_end
5024	DO I=i_start,i_end
5025	DO K=k_start,k_end
5026	RUSHTEN(I,K,J) =RUSHTEN(I,K,J)/mu(I,J)
5027	RVSHTEN(I,K,J) =RVSHTEN(I,K,J)/mu(I,J)
5028	RTHSHTEN(I,K,J)=RTHSHTEN(I,K,J)/mu(I,J)
5029	ENDDO
5030	ENDDO
5031	ENDDO
5032
5033	IF (P_QV .ge. PARAM_FIRST_SCALAR)THEN
5034	DO J=j_start,j_end
5035	DO I=i_start,i_end
5036	DO K=k_start,k_end
5037	RQVSHTEN(I,K,J)=RQVSHTEN(I,K,J)/mu(I,J)
5038	ENDDO
5039	ENDDO
5040	ENDDO
5041	ENDIF
5042
5043	IF (P_QC .ge. PARAM_FIRST_SCALAR)THEN
5044	DO J=j_start,j_end
5045	DO I=i_start,i_end
5046	DO K=k_start,k_end
5047	RQCSHTEN(I,K,J)=RQCSHTEN(I,K,J)/mu(I,J)
5048	ENDDO
5049	ENDDO
5050	ENDDO
5051	ENDIF
5052
5053	IF (P_QR .ge. PARAM_FIRST_SCALAR)THEN
5054	DO J=j_start,j_end
5055	DO I=i_start,i_end
5056	DO K=k_start,k_end
5057	RQRSHTEN(I,K,J)=RQRSHTEN(I,K,J)/mu(I,J)
5058	ENDDO
5059	ENDDO
5060	ENDDO
5061	ENDIF
5062
5063	IF (P_QI .ge. PARAM_FIRST_SCALAR)THEN
5064	DO J=j_start,j_end
5065	DO I=i_start,i_end
5066	DO K=k_start,k_end
5067	RQISHTEN(I,K,J)=RQISHTEN(I,K,J)/mu(I,J)
5068	ENDDO
5069	ENDDO
5070	ENDDO
5071	ENDIF
5072
5073	IF(P_QS .ge. PARAM_FIRST_SCALAR)THEN
5074	DO J=j_start,j_end
5075	DO I=i_start,i_end
5076	DO K=k_start,k_end
5077	RQSSHTEN(I,K,J)=RQSSHTEN(I,K,J)/mu(I,J)
5078	ENDDO
5079	ENDDO
5080	ENDDO
5081	ENDIF
5082
5083	IF(P_QG .ge. PARAM_FIRST_SCALAR)THEN
5084	DO J=j_start,j_end
5085	DO I=i_start,i_end
5086	DO K=k_start,k_end
5087	RQGSHTEN(I,K,J)=RQGSHTEN(I,K,J)/mu(I,J)
5088	ENDDO
5089	ENDDO
5090	ENDDO
5091	ENDIF
5092
5093	ENDIF
5094
5095	IF (config_flags%bl_pbl_physics .gt. 0) THEN
5096
5097	DO J=j_start,j_end
5098	DO K=k_start,k_end
5099	DO I=i_start,i_end
5100	RUBLTEN(I,K,J) =RUBLTEN(I,K,J)/mu(I,J)
5101	RVBLTEN(I,K,J) =RVBLTEN(I,K,J)/mu(I,J)
5102	RTHBLTEN(I,K,J)=RTHBLTEN(I,K,J)/mu(I,J)
5103	ENDDO
5104	ENDDO
5105	ENDDO
5106
5107	IF (P_QV .ge. PARAM_FIRST_SCALAR) THEN
5108	DO J=j_start,j_end
5109	DO K=k_start,k_end
5110	DO I=i_start,i_end
5111	RQVBLTEN(I,K,J)=RQVBLTEN(I,K,J)/mu(I,J)
5112	ENDDO
5113	ENDDO
5114	ENDDO
5115	ENDIF
5116
5117	IF (P_QC .ge. PARAM_FIRST_SCALAR) THEN
5118	DO J=j_start,j_end
5119	DO K=k_start,k_end
5120	DO I=i_start,i_end
5121	RQCBLTEN(I,K,J)=RQCBLTEN(I,K,J)/mu(I,J)
5122	ENDDO
5123	ENDDO
5124	ENDDO
5125	ENDIF
5126
5127	IF (P_QI .ge. PARAM_FIRST_SCALAR) THEN
5128	DO J=j_start,j_end
5129	DO K=k_start,k_end
5130	DO I=i_start,i_end
5131	RQIBLTEN(I,K,J)=RQIBLTEN(I,K,J)/mu(I,J)
5132	ENDDO
5133	ENDDO
5134	ENDDO
5135	ENDIF
5136
5137	ENDIF
5138
5139	! decouple advective forcing required by Grell-Devenyi scheme
5140
5141	if(( config_flags%cu_physics == GDSCHEME ) .OR. &
5142	( config_flags%cu_physics == G3SCHEME ) .OR. &
5143	( config_flags%cu_physics == KFETASCHEME ) .OR. &
5144	( config_flags%cu_physics == TIEDTKESCHEME ) ) then ! Tiedtke ZCX&YQW
5145
5146	DO J=j_start,j_end
5147	DO I=i_start,i_end
5148	DO K=k_start,k_end
5149	RTHFTEN(I,K,J)=RTHFTEN(I,K,J)/mu(I,J)
5150	ENDDO
5151	ENDDO
5152	ENDDO
5153
5154	IF (P_QV .ge. PARAM_FIRST_SCALAR)THEN
5155	DO J=j_start,j_end
5156	DO I=i_start,i_end
5157	DO K=k_start,k_end
5158	RQVFTEN(I,K,J)=RQVFTEN(I,K,J)/mu(I,J)
5159	ENDDO
5160	ENDDO
5161	ENDDO
5162	ENDIF
5163
5164	END IF
5165
5166	! fdda
5167	! note fdda u and v tendencies are staggered, also only interior points have muu/muv,
5168	! so only decouple those
5169
5170	IF (config_flags%grid_fdda .gt. 0) THEN
5171
5172	i_startu=MAX(its,ids+1)
5173	j_startv=MAX(jts,jds+1)
5174
5175	DO J=j_start,j_end
5176	DO K=k_start,k_end
5177	DO I=i_startu,i_end
5178	RUNDGDTEN(I,K,J) =RUNDGDTEN(I,K,J)/muu(I,J)
5179	ENDDO
5180	ENDDO
5181	ENDDO
5182	DO J=j_startv,j_end
5183	DO K=k_start,k_end
5184	DO I=i_start,i_end
5185	RVNDGDTEN(I,K,J) =RVNDGDTEN(I,K,J)/muv(I,J)
5186	ENDDO
5187	ENDDO
5188	ENDDO
5189	DO J=j_start,j_end
5190	DO K=k_start,k_end
5191	DO I=i_start,i_end
5192	RTHNDGDTEN(I,K,J)=RTHNDGDTEN(I,K,J)/mu(I,J)
5193	! RMUNDGDTEN(I,J) - no coupling
5194	ENDDO
5195	ENDDO
5196	ENDDO
5197
5198	IF (config_flags%grid_fdda .EQ. 2) THEN
5199	DO J=j_start,j_end
5200	DO K=k_start,kte
5201	DO I=i_start,i_end
5202	RPHNDGDTEN(I,K,J)=RPHNDGDTEN(I,K,J)/mu(I,J)
5203	ENDDO
5204	ENDDO
5205	ENDDO
5206
5207	ELSE IF (config_flags%grid_fdda .EQ. 1) THEN
5208	IF (P_QV .ge. PARAM_FIRST_SCALAR) THEN
5209	DO J=j_start,j_end
5210	DO K=k_start,k_end
5211	DO I=i_start,i_end
5212	RQVNDGDTEN(I,K,J)=RQVNDGDTEN(I,K,J)/mu(I,J)
5213	ENDDO
5214	ENDDO
5215	ENDDO
5216	ENDIF
5217	ENDIF
5218
5219	ENDIF
5220
5221	END SUBROUTINE phy_prep
5222
5223	!------------------------------------------------------------
5224
5225	SUBROUTINE moist_physics_prep_em( t_new, t_old, t0, rho, al, alb, &
5226	p, p8w, p0, pb, ph, phb, &
5227	th_phy, pii, pf, &
5228	z, z_at_w, dz8w, &
5229	dt,h_diabatic, &
5230	config_flags,fzm, fzp, &
5231	ids,ide, jds,jde, kds,kde, &
5232	ims,ime, jms,jme, kms,kme, &
5233	its,ite, jts,jte, kts,kte )
5234
5235	IMPLICIT NONE
5236
5237	! Here we construct full fields
5238	! needed by the microphysics
5239
5240	TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
5241
5242	INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde
5243	INTEGER, INTENT(IN ) :: ims,ime, jms,jme, kms,kme
5244	INTEGER, INTENT(IN ) :: its,ite, jts,jte, kts,kte
5245
5246	REAL, INTENT(IN ) :: dt
5247
5248	REAL, DIMENSION( ims:ime , kms:kme, jms:jme ), &
5249	INTENT(IN ) :: al, &
5250	alb, &
5251	p, &
5252	pb, &
5253	ph, &
5254	phb
5255
5256
5257	REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fzm, &
5258	fzp
5259
5260	REAL, DIMENSION( ims:ime , kms:kme, jms:jme ), &
5261	INTENT( OUT) :: rho, &
5262	th_phy, &
5263	pii, &
5264	pf, &
5265	z, &
5266	z_at_w, &
5267	dz8w, &
5268	p8w
5269
5270	REAL, DIMENSION( ims:ime , kms:kme, jms:jme ), &
5271	INTENT(INOUT) :: h_diabatic
5272
5273	REAL, DIMENSION( ims:ime , kms:kme, jms:jme ), &
5274	INTENT(INOUT) :: t_new, &
5275	t_old
5276
5277	REAL, INTENT(IN ) :: t0, p0
5278	REAL :: z0,z1,z2,w1,w2
5279
5280	INTEGER :: i_start, i_end, j_start, j_end, k_start, k_end
5281	INTEGER :: i, j, k
5282
5283	!--------------------------------------------------------------------
5284
5285	!<DESCRIPTION>
5286	!
5287	! moist_phys_prep_em calculates a number of diagnostic quantities needed by
5288	! the microphysics routines.
5289	!
5290	!</DESCRIPTION>
5291
5292	! set up loop bounds for this grid's boundary conditions
5293
5294	i_start = its
5295	i_end = min( ite,ide-1 )
5296	j_start = jts
5297	j_end = min( jte,jde-1 )
5298
5299	k_start = kts
5300	k_end = min( kte, kde-1 )
5301
5302	DO j = j_start, j_end
5303	DO k = k_start, kte
5304	DO i = i_start, i_end
5305	z_at_w(i,k,j) = (ph(i,k,j)+phb(i,k,j))/g
5306	ENDDO
5307	ENDDO
5308	ENDDO
5309
5310	do j = j_start,j_end
5311	do k = k_start, kte-1
5312	do i = i_start, i_end
5313	dz8w(i,k,j) = z_at_w(i,k+1,j)-z_at_w(i,k,j)
5314	enddo
5315	enddo
5316	enddo
5317
5318	do j = j_start,j_end
5319	do i = i_start, i_end
5320	dz8w(i,kte,j) = 0.
5321	enddo
5322	enddo
5323
5324
5325	! compute full pii, rho, and z at the new time-level
5326	! (needed for physics).
5327	! convert perturbation theta to full theta (th_phy)
5328	! use h_diabatic to temporarily save pre-microphysics full theta
5329
5330	DO j = j_start, j_end
5331	DO k = k_start, k_end
5332	DO i = i_start, i_end
5333
5334	#ifdef REVERT
5335	t_new(i,k,j) = t_new(i,k,j)-h_diabatic(i,k,j)*dt
5336	#endif
5337	th_phy(i,k,j) = t_new(i,k,j) + t0
5338	h_diabatic(i,k,j) = th_phy(i,k,j)
5339	rho(i,k,j) = 1./(al(i,k,j)+alb(i,k,j))
5340	pii(i,k,j) = ((p(i,k,j)+pb(i,k,j))/p0)**rcp
5341	z(i,k,j) = 0.5*(z_at_w(i,k,j) +z_at_w(i,k+1,j) )
5342	pf(i,k,j) = p(i,k,j)+pb(i,k,j)
5343
5344	ENDDO
5345	ENDDO
5346	ENDDO
5347
5348	! interp t and p at w points
5349
5350	do j = j_start,j_end
5351	do k = 2, k_end
5352	do i = i_start, i_end
5353	p8w(i,k,j) = fzm(k)pf(i,k,j)+fzp(k)pf(i,k-1,j)
5354	enddo
5355	enddo
5356	enddo
5357
5358	! extrapolate p and t to surface and top.
5359	! we'll use an extrapolation in z for now
5360
5361	do j = j_start,j_end
5362	do i = i_start, i_end
5363
5364	! bottom
5365
5366	z0 = z_at_w(i,1,j)
5367	z1 = z(i,1,j)
5368	z2 = z(i,2,j)
5369	w1 = (z0 - z2)/(z1 - z2)
5370	w2 = 1. - w1
5371	p8w(i,1,j) = w1pf(i,1,j)+w2pf(i,2,j)
5372
5373	! top
5374
5375	z0 = z_at_w(i,kte,j)
5376	z1 = z(i,k_end,j)
5377	z2 = z(i,k_end-1,j)
5378	w1 = (z0 - z2)/(z1 - z2)
5379	w2 = 1. - w1
5380	! p8w(i,kde,j) = w1pf(i,kde-1,j)+w2pf(i,kde-2,j)
5381	p8w(i,kde,j) = exp(w1log(pf(i,kde-1,j))+w2log(pf(i,kde-2,j)))
5382
5383	enddo
5384	enddo
5385
5386	END SUBROUTINE moist_physics_prep_em
5387
5388	!------------------------------------------------------------------------------
5389
5390	SUBROUTINE moist_physics_finish_em( t_new, t_old, t0, mut, &
5391	th_phy, h_diabatic, dt, &
5392	config_flags, &
5393	#if ( WRF_DFI_RADAR == 1 )
5394	dfi_tten_rad,dfi_stage, &
5395	#endif
5396	ids,ide, jds,jde, kds,kde, &
5397	ims,ime, jms,jme, kms,kme, &
5398	its,ite, jts,jte, kts,kte )
5399
5400	IMPLICIT NONE
5401
5402	! Here we construct full fields
5403	! needed by the microphysics
5404
5405	TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
5406
5407	INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde
5408	INTEGER, INTENT(IN ) :: ims,ime, jms,jme, kms,kme
5409	INTEGER, INTENT(IN ) :: its,ite, jts,jte, kts,kte
5410
5411	REAL, DIMENSION( ims:ime , kms:kme, jms:jme ), &
5412	INTENT(INOUT) :: t_new, &
5413	t_old, &
5414	th_phy, &
5415	h_diabatic
5416	#if ( WRF_DFI_RADAR == 1 )
5417	REAL, DIMENSION( ims:ime , kms:kme, jms:jme ), &
5418	INTENT(IN), OPTIONAL :: dfi_tten_rad
5419	INTEGER, INTENT(IN ) ,OPTIONAL :: dfi_stage
5420	REAL :: dfi_tten_max, old_max
5421	#endif
5422
5423	REAL mpten, mptenmax, mptenmin
5424
5425	REAL, DIMENSION( ims:ime , jms:jme ), INTENT(INOUT) :: mut
5426
5427
5428	REAL, INTENT(IN ) :: t0, dt
5429
5430	INTEGER :: i_start, i_end, j_start, j_end, k_start, k_end
5431	INTEGER :: i, j, k, imax, jmax, imin, jmin
5432
5433	!--------------------------------------------------------------------
5434
5435	!<DESCRIPTION>
5436	!
5437	! moist_phys_finish_em resets theta to its perturbation value and
5438	! computes and stores the microphysics diabatic heating term.
5439	!
5440	!</DESCRIPTION>
5441
5442	! set up loop bounds for this grid's boundary conditions
5443
5444
5445	i_start = its
5446	i_end = min( ite,ide-1 )
5447	j_start = jts
5448	j_end = min( jte,jde-1 )
5449	! i_start=max(its,ids+4)
5450	! i_end=min(ite,ide-5)
5451	! j_start=max(jts,jds+4)
5452	! j_end=min(jte,jde-5)
5453
5454	k_start = kts
5455	k_end = min( kte, kde-1 )
5456
5457	#if ( WRF_DFI_RADAR == 1 )
5458	IF ( PRESENT(dfi_stage) .and. PRESENT(dfi_tten_rad) ) THEN
5459	IF ( dfi_stage ==DFI_FWD ) THEN
5460	WRITE(wrf_err_message,*)'Add radar tendency: i_start,j_start: ', i_start, j_start
5461	CALL wrf_debug ( 100 , TRIM(wrf_err_message) )
5462	ENDIF
5463	ENDIF
5464	dfi_tten_max=-999
5465	old_max=-999
5466	#endif
5467
5468	! add microphysics theta diff to perturbation theta, set h_diabatic
5469
5470	IF ( config_flags%no_mp_heating .eq. 0 ) THEN
5471	mptenmax = 0.
5472	mptenmin = 999.
5473	DO j = j_start, j_end
5474	DO k = k_start, k_end
5475	DO i = i_start, i_end
5476	mpten = th_phy(i,k,j)-h_diabatic(i,k,j)
5477	#if ( WRF_DFI_RADAR == 1 )
5478	if(mpten.gt.mptenmax) then
5479	mptenmax=mpten
5480	imax=i
5481	jmax=j
5482	endif
5483	if(mpten.lt.mptenmin) then
5484	mptenmin=mpten
5485	imin=i
5486	jmin=j
5487	endif
5488	mpten=min(config_flags%mp_tend_lim*dt, mpten)
5489	mpten=max(-config_flags%mp_tend_lim*dt, mpten)
5490
5491	if(k < k_end ) then
5492	if(dfi_tten_max < dfi_tten_rad(i,k,j) ) dfi_tten_max = dfi_tten_rad(i,k,j)
5493	if(old_max < (th_phy(i,k,j)-h_diabatic(i,k,j)) ) old_max=th_phy(i,k,j)-h_diabatic(i,k,j)
5494	endif
5495
5496	IF ( PRESENT(dfi_stage) .and. PRESENT(dfi_tten_rad) ) THEN
5497	IF ( dfi_stage == DFI_FWD .and. dfi_tten_rad(i,k,j) >= -0.1 .and. &
5498	dfi_tten_rad(i,k,j) <= 0.1 .and. k < k_end ) THEN
5499	! add radar temp tendency
5500	! there is radar coverage
5501	t_new(i,k,j) = t_new(i,k,j) + (dfi_tten_rad(i,k,j))*dt
5502	ELSE
5503	! no radar coverage
5504	t_new(i,k,j) = t_new(i,k,j) + mpten
5505	ENDIF
5506	ENDIF
5507	#else
5508	t_new(i,k,j) = t_new(i,k,j) + mpten
5509	#endif
5510	h_diabatic(i,k,j) = mpten/dt
5511	ENDDO
5512	ENDDO
5513	ENDDO
5514
5515	ELSE
5516
5517	DO j = j_start, j_end
5518	DO k = k_start, k_end
5519	DO i = i_start, i_end
5520	! t_new(i,k,j) = t_new(i,k,j)
5521	h_diabatic(i,k,j) = 0.
5522	ENDDO
5523	ENDDO
5524	ENDDO
5525	ENDIF
5526
5527	END SUBROUTINE moist_physics_finish_em
5528
5529	!----------------------------------------------------------------
5530
5531
5532	SUBROUTINE init_module_big_step
5533	END SUBROUTINE init_module_big_step
5534
5535	SUBROUTINE set_tend ( field, field_adv_tend, msf, &
5536	ids, ide, jds, jde, kds, kde, &
5537	ims, ime, jms, jme, kms, kme, &
5538	its, ite, jts, jte, kts, kte )
5539
5540	IMPLICIT NONE
5541
5542	! Input data
5543
5544	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
5545	ims, ime, jms, jme, kms, kme, &
5546	its, ite, jts, jte, kts, kte
5547
5548	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT) :: field
5549
5550	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN) :: field_adv_tend
5551
5552	REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN) :: msf
5553
5554	! Local data
5555
5556	INTEGER :: i, j, k, itf, jtf, ktf
5557
5558	!<DESCRIPTION>
5559	!
5560	! set_tend copies the advective tendency array into the tendency array.
5561	!
5562	!</DESCRIPTION>
5563
5564	jtf = MIN(jte,jde-1)
5565	ktf = MIN(kte,kde-1)
5566	itf = MIN(ite,ide-1)
5567	DO j = jts, jtf
5568	DO k = kts, ktf
5569	DO i = its, itf
5570	field(i,k,j) = field_adv_tend(i,k,j)*msf(i,j)
5571	ENDDO
5572	ENDDO
5573	ENDDO
5574
5575	END SUBROUTINE set_tend
5576
5577	!------------------------------------------------------------------------------
5578
5579	SUBROUTINE rk_rayleigh_damp( ru_tendf, rv_tendf, &
5580	rw_tendf, t_tendf, &
5581	u, v, w, t, t_init, &
5582	mut, muu, muv, ph, phb, &
5583	u_base, v_base, t_base, z_base, &
5584	dampcoef, zdamp, &
5585	ids, ide, jds, jde, kds, kde, &
5586	ims, ime, jms, jme, kms, kme, &
5587	its, ite, jts, jte, kts, kte )
5588
5589	! History: Apr 2005 Modifications by George Bryan, NCAR:
5590	! - Generalized the code in a way that allows for
5591	! simulations with steep terrain.
5592	!
5593	! Jul 2004 Modifications by George Bryan, NCAR:
5594	! - Modified the code to use u_base, v_base, and t_base
5595	! arrays for the background state. Removed the hard-wired
5596	! base-state values.
5597	! - Modified the code to use dampcoef, zdamp, and damp_opt,
5598	! i.e., the upper-level damper variables in namelist.input.
5599	! Removed the hard-wired variables in the older version.
5600	! This damper is used when damp_opt = 2.
5601	! - Modified the code to account for the movement of the
5602	! model surfaces with time. The code now obtains a base-
5603	! state value by interpolation using the "_base" arrays.
5604
5605	! Nov 2003 Bug fix by Jason Knievel, NCAR
5606
5607	! Aug 2003 Meridional dimension, some comments, and
5608	! changes in layout of the code added by
5609	! Jason Knievel, NCAR
5610
5611	! Jul 2003 Original code by Bill Skamarock, NCAR
5612
5613	! Purpose: This routine applies Rayleigh damping to a layer at top
5614	! of the model domain.
5615
5616	!-----------------------------------------------------------------------
5617	! Begin declarations.
5618
5619	IMPLICIT NONE
5620
5621	INTEGER, INTENT( IN ) &
5622	:: ids, ide, jds, jde, kds, kde, &
5623	ims, ime, jms, jme, kms, kme, &
5624	its, ite, jts, jte, kts, kte
5625
5626	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( INOUT ) &
5627	:: ru_tendf, rv_tendf, rw_tendf, t_tendf
5628
5629	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( IN ) &
5630	:: u, v, w, t, t_init, ph, phb
5631
5632	REAL, DIMENSION( ims:ime, jms:jme ), INTENT( IN ) &
5633	:: mut, muu, muv
5634
5635	REAL, DIMENSION( kms:kme ) , INTENT(IN ) &
5636	:: u_base, v_base, t_base, z_base
5637
5638	REAL, INTENT(IN ) &
5639	:: dampcoef, zdamp
5640
5641	! Local variables.
5642
5643	INTEGER &
5644	:: i_start, i_end, j_start, j_end, k_start, k_end, i, j, k, ktf, k1, k2
5645
5646	REAL &
5647	:: pii, dcoef, z, ztop
5648
5649	REAL :: wkp1, wk, wkm1
5650
5651	REAL, DIMENSION( kms:kme ) :: z00, u00, v00, t00
5652
5653	! End declarations.
5654	!-----------------------------------------------------------------------
5655
5656	pii = 2.0 * asin(1.0)
5657
5658	ktf = MIN( kte, kde-1 )
5659
5660	!-----------------------------------------------------------------------
5661	! Adjust u to base state.
5662
5663	DO j = jts, MIN( jte, jde-1 )
5664	DO i = its, MIN( ite, ide )
5665
5666	! Get height at top of model
5667	ztop = 0.5*( phb(i ,kde,j)+phb(i-1,kde,j) &
5668	+ph(i ,kde,j)+ph(i-1,kde,j) )/g
5669
5670	! Find bottom of damping layer
5671	k1 = ktf
5672	z = ztop
5673	DO WHILE( z >= (ztop-zdamp) )
5674	z = 0.25*( phb(i ,k1,j)+phb(i ,k1+1,j) &
5675	+phb(i-1,k1,j)+phb(i-1,k1+1,j) &
5676	+ph(i ,k1,j)+ph(i ,k1+1,j) &
5677	+ph(i-1,k1,j)+ph(i-1,k1+1,j))/g
5678	z00(k1) = z
5679	k1 = k1 - 1
5680	ENDDO
5681	k1 = k1 + 2
5682
5683	! Get reference state at model levels
5684	DO k = k1, ktf
5685	k2 = ktf
5686	DO WHILE( z_base(k2) .gt. z00(k) )
5687	k2 = k2 - 1
5688	ENDDO
5689	if(k2+1.gt.ktf)then
5690	u00(k) = u_base(k2) + ( u_base(k2) - u_base(k2-1) ) &
5691	* ( z00(k) - z_base(k2) ) &
5692	/ ( z_base(k2) - z_base(k2-1) )
5693	else
5694	u00(k) = u_base(k2) + ( u_base(k2+1) - u_base(k2) ) &
5695	* ( z00(k) - z_base(k2) ) &
5696	/ ( z_base(k2+1) - z_base(k2) )
5697	endif
5698	ENDDO
5699
5700	! Apply the Rayleigh damper
5701	DO k = k1, ktf
5702	dcoef = 1.0 - MIN( 1.0, ( ztop - z00(k) ) / zdamp )
5703	dcoef = (SIN( 0.5 * pii * dcoef ) )**2
5704	ru_tendf(i,k,j) = ru_tendf(i,k,j) - &
5705	muu(i,j) * ( dcoef * dampcoef ) * &
5706	( u(i,k,j) - u00(k) )
5707	END DO
5708
5709	END DO
5710	END DO
5711
5712	! End adjustment of u.
5713	!-----------------------------------------------------------------------
5714
5715	!-----------------------------------------------------------------------
5716	! Adjust v to base state.
5717
5718	DO j = jts, MIN( jte, jde )
5719	DO i = its, MIN( ite, ide-1 )
5720
5721	! Get height at top of model
5722	ztop = 0.5*( phb(i,kde,j )+phb(i,kde,j-1) &
5723	+ph(i,kde,j )+ph(i,kde,j-1) )/g
5724
5725	! Find bottom of damping layer
5726	k1 = ktf
5727	z = ztop
5728	DO WHILE( z >= (ztop-zdamp) )
5729	z = 0.25*( phb(i,k1,j )+phb(i,k1+1,j ) &
5730	+phb(i,k1,j-1)+phb(i,k1+1,j-1) &
5731	+ph(i,k1,j )+ph(i,k1+1,j ) &
5732	+ph(i,k1,j-1)+ph(i,k1+1,j-1))/g
5733	z00(k1) = z
5734	k1 = k1 - 1
5735	ENDDO
5736	k1 = k1 + 2
5737
5738	! Get reference state at model levels
5739	DO k = k1, ktf
5740	k2 = ktf
5741	DO WHILE( z_base(k2) .gt. z00(k) )
5742	k2 = k2 - 1
5743	ENDDO
5744	if(k2+1.gt.ktf)then
5745	v00(k) = v_base(k2) + ( v_base(k2) - v_base(k2-1) ) &
5746	* ( z00(k) - z_base(k2) ) &
5747	/ ( z_base(k2) - z_base(k2-1) )
5748	else
5749	v00(k) = v_base(k2) + ( v_base(k2+1) - v_base(k2) ) &
5750	* ( z00(k) - z_base(k2) ) &
5751	/ ( z_base(k2+1) - z_base(k2) )
5752	endif
5753	ENDDO
5754
5755	! Apply the Rayleigh damper
5756	DO k = k1, ktf
5757	dcoef = 1.0 - MIN( 1.0, ( ztop - z00(k) ) / zdamp )
5758	dcoef = (SIN( 0.5 * pii * dcoef ) )**2
5759	rv_tendf(i,k,j) = rv_tendf(i,k,j) - &
5760	muv(i,j) * ( dcoef * dampcoef ) * &
5761	( v(i,k,j) - v00(k) )
5762	END DO
5763
5764	END DO
5765	END DO
5766
5767	! End adjustment of v.
5768	!-----------------------------------------------------------------------
5769
5770	!-----------------------------------------------------------------------
5771	! Adjust w to base state.
5772
5773	DO j = jts, MIN( jte, jde-1 )
5774	DO i = its, MIN( ite, ide-1 )
5775	ztop = ( phb(i,kde,j) + ph(i,kde,j) ) / g
5776	DO k = kts, MIN( kte, kde )
5777	z = ( phb(i,k,j) + ph(i,k,j) ) / g
5778	IF ( z >= (ztop-zdamp) ) THEN
5779	dcoef = 1.0 - MIN( 1.0, ( ztop - z ) / zdamp )
5780	dcoef = ( SIN( 0.5 * pii * dcoef ) )**2
5781	rw_tendf(i,k,j) = rw_tendf(i,k,j) - &
5782	mut(i,j) * ( dcoef * dampcoef ) * w(i,k,j)
5783	END IF
5784	END DO
5785	END DO
5786	END DO
5787
5788	! End adjustment of w.
5789	!-----------------------------------------------------------------------
5790
5791	!-----------------------------------------------------------------------
5792	! Adjust potential temperature to base state.
5793
5794	DO j = jts, MIN( jte, jde-1 )
5795	DO i = its, MIN( ite, ide-1 )
5796
5797	! Get height at top of model
5798	ztop = ( phb(i,kde,j) + ph(i,kde,j) ) / g
5799
5800	! Find bottom of damping layer
5801	k1 = ktf
5802	z = ztop
5803	DO WHILE( z >= (ztop-zdamp) )
5804	z = 0.5 * ( phb(i,k1,j) + phb(i,k1+1,j) + &
5805	ph(i,k1,j) + ph(i,k1+1,j) ) / g
5806	z00(k1) = z
5807	k1 = k1 - 1
5808	ENDDO
5809	k1 = k1 + 2
5810
5811	! Get reference state at model levels
5812	DO k = k1, ktf
5813	k2 = ktf
5814	DO WHILE( z_base(k2) .gt. z00(k) )
5815	k2 = k2 - 1
5816	ENDDO
5817	if(k2+1.gt.ktf)then
5818	t00(k) = t_base(k2) + ( t_base(k2) - t_base(k2-1) ) &
5819	* ( z00(k) - z_base(k2) ) &
5820	/ ( z_base(k2) - z_base(k2-1) )
5821	else
5822	t00(k) = t_base(k2) + ( t_base(k2+1) - t_base(k2) ) &
5823	* ( z00(k) - z_base(k2) ) &
5824	/ ( z_base(k2+1) - z_base(k2) )
5825	endif
5826	ENDDO
5827
5828	! Apply the Rayleigh damper
5829	DO k = k1, ktf
5830	dcoef = 1.0 - MIN( 1.0, ( ztop - z00(k) ) / zdamp )
5831	dcoef = (SIN( 0.5 * pii * dcoef ) )**2
5832	t_tendf(i,k,j) = t_tendf(i,k,j) - &
5833	mut(i,j) * ( dcoef * dampcoef ) * &
5834	( t(i,k,j) - t00(k) )
5835	END DO
5836
5837	END DO
5838	END DO
5839
5840	! End adjustment of potential temperature.
5841	!-----------------------------------------------------------------------
5842
5843	END SUBROUTINE rk_rayleigh_damp
5844
5845	!==============================================================================
5846	!==============================================================================
5847
5848	SUBROUTINE theta_relaxation( t_tendf, t, t_init, &
5849	mut, ph, phb, &
5850	t_base, z_base, &
5851	ids, ide, jds, jde, kds, kde, &
5852	ims, ime, jms, jme, kms, kme, &
5853	its, ite, jts, jte, kts, kte )
5854
5855	! Purpose: Newtonian relaxation on potential temperature. Serves two
5856	! purposes: 1) to mimic atmospheric radiation in a simple
5857	! manner, and 2) to keep the vertical profile of temperature
5858	! close to the initial (base-state) profile, which is useful
5859	! for certain idealized applications.
5860
5861	! Reference: Rotunno and Emanuel, 1987, JAS, p. 546
5862
5863	!-----------------------------------------------------------------------
5864	! Begin declarations.
5865
5866	IMPLICIT NONE
5867
5868	INTEGER, INTENT( IN ) &
5869	:: ids, ide, jds, jde, kds, kde, &
5870	ims, ime, jms, jme, kms, kme, &
5871	its, ite, jts, jte, kts, kte
5872
5873	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( INOUT ) &
5874	:: t_tendf
5875
5876	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( IN ) &
5877	:: t, t_init, ph, phb
5878
5879	REAL, DIMENSION( ims:ime, jms:jme ), INTENT( IN ) &
5880	:: mut
5881
5882	REAL, DIMENSION( kms:kme ) , INTENT(IN ) &
5883	:: t_base, z_base
5884
5885	! Local variables.
5886
5887	INTEGER :: i, j, k, ktf, k2
5888	REAL :: tau_r , rmax , rmin , inv_tau_r , inv_g , rterm
5889	REAL, DIMENSION( kms:kme ) :: z00,t00
5890
5891	! End declarations.
5892	!-----------------------------------------------------------------------
5893
5894	! set tau_r to 12 h, following RE87
5895	tau_r = 12.0*3600.0
5896
5897	! limit rterm to +/- 2 K/day
5898	rmax = 2.0/86400.0
5899	rmin = -rmax
5900
5901	ktf = MIN( kte, kde-1 )
5902	inv_tau_r = 1.0/tau_r
5903	inv_g = 1.0/g
5904
5905	!-----------------------------------------------------------------------
5906	! Adjust potential temperature to base state.
5907
5908	DO j = jts, MIN( jte, jde-1 )
5909	DO i = its, MIN( ite, ide-1 )
5910
5911	! Get height of model levels:
5912	DO k = kts, ktf
5913	z00(k) = 0.5 * ( phb(i,k,j) + phb(i,k+1,j) + &
5914	ph(i,k,j) + ph(i,k+1,j) ) * inv_g
5915	ENDDO
5916
5917	! Get reference state:
5918	DO k = kts, ktf
5919	k2 = ktf
5920	DO WHILE( z_base(k2) .gt. z00(k) .and. k2 .gt. 1 )
5921	k2 = k2 - 1
5922	ENDDO
5923	if(k2+1.gt.ktf)then
5924	t00(k) = t_base(k2) + ( t_base(k2) - t_base(k2-1) ) &
5925	* ( z00(k) - z_base(k2) ) &
5926	/ ( z_base(k2) - z_base(k2-1) )
5927	else
5928	t00(k) = t_base(k2) + ( t_base(k2+1) - t_base(k2) ) &
5929	* ( z00(k) - z_base(k2) ) &
5930	/ ( z_base(k2+1) - z_base(k2) )
5931	endif
5932	ENDDO
5933
5934	! Apply the RE87 R term:
5935	DO k = kts, ktf
5936	rterm = -( t(i,k,j) - t00(k) )*inv_tau_r
5937	! limit rterm:
5938	rterm = min( rterm , rmax )
5939	rterm = max( rterm , rmin )
5940	t_tendf(i,k,j) = t_tendf(i,k,j) + mut(i,j)*rterm
5941	END DO
5942
5943	END DO
5944	END DO
5945
5946	END SUBROUTINE theta_relaxation
5947
5948	!==============================================================================
5949	!==============================================================================
5950
5951	SUBROUTINE sixth_order_diffusion( name, field, tendency, mu, dt, &
5952	config_flags, &
5953	diff_6th_opt, diff_6th_factor, &
5954	ids, ide, jds, jde, kds, kde, &
5955	ims, ime, jms, jme, kms, kme, &
5956	its, ite, jts, jte, kts, kte )
5957
5958	! History: 14 Nov 2006 Name of variable changed by Jason Knievel
5959	! 07 Jun 2006 Revised and generalized by Jason Knievel
5960	! 25 Apr 2005 Original code by Jason Knievel, NCAR
5961
5962	! Purpose: Apply 6th-order, monotonic (flux-limited), numerical
5963	! diffusion to 3-d velocity and to scalars.
5964
5965	! References: Ming Xue (MWR Aug 2000)
5966	! Durran ("Numerical Methods for Wave Equations..." 1999)
5967	! George Bryan (personal communication)
5968
5969	!------------------------------------------------------------------------------
5970	! Begin: Declarations.
5971
5972	IMPLICIT NONE
5973
5974	INTEGER, INTENT(IN) &
5975	:: ids, ide, jds, jde, kds, kde, &
5976	ims, ime, jms, jme, kms, kme, &
5977	its, ite, jts, jte, kts, kte
5978
5979	TYPE(grid_config_rec_type), INTENT(IN) &
5980	:: config_flags
5981
5982	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) &
5983	:: tendency
5984
5985	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN) &
5986	:: field
5987
5988	REAL, DIMENSION( ims:ime , jms:jme ), INTENT(IN) &
5989	:: mu
5990
5991	REAL, INTENT(IN) &
5992	:: dt
5993
5994	REAL, INTENT(IN) &
5995	:: diff_6th_factor
5996
5997	INTEGER, INTENT(IN) &
5998	:: diff_6th_opt
5999
6000	CHARACTER(LEN=1) , INTENT(IN) &
6001	:: name
6002
6003	INTEGER &
6004	:: i, j, k, &
6005	i_start, i_end, &
6006	j_start, j_end, &
6007	k_start, k_end, &
6008	ktf
6009
6010	REAL &
6011	:: dflux_x_p0, dflux_y_p0, &
6012	dflux_x_p1, dflux_y_p1, &
6013	tendency_x, tendency_y, &
6014	mu_avg_p0, mu_avg_p1, &
6015	diff_6th_coef
6016
6017	LOGICAL &
6018	:: specified
6019
6020	! End: Declarations.
6021	!------------------------------------------------------------------------------
6022
6023	!------------------------------------------------------------------------------
6024	! Begin: Translate the diffusion factor into a diffusion coefficient. See
6025	! Durran's text, section 2.4.3, then adjust for sixth-order diffusion (not
6026	! fourth) and for diffusion in two dimensions (not one). For reference, a
6027	! factor of 1.0 would mean complete diffusion of a 2dx wave in one time step,
6028	! although application of the flux limiter reduces somewhat the effects of
6029	! diffusion for a given coefficient.
6030
6031	diff_6th_coef = diff_6th_factor * 0.015625 / ( 2.0 * dt )
6032
6033	! End: Translate diffusion factor.
6034	!------------------------------------------------------------------------------
6035
6036	!------------------------------------------------------------------------------
6037	! Begin: Assign limits of spatial loops depending on variable to be diffused.
6038	! The halo regions are already filled with values by the time this subroutine
6039	! is called, which allows the stencil to extend beyond the domains' edges.
6040
6041	ktf = MIN( kte, kde-1 )
6042
6043	IF ( name .EQ. 'u' ) THEN
6044
6045	i_start = its
6046	i_end = ite
6047	j_start = jts
6048	j_end = MIN(jde-1,jte)
6049	k_start = kts
6050	k_end = ktf
6051
6052	ELSE IF ( name .EQ. 'v' ) THEN
6053
6054	i_start = its
6055	i_end = MIN(ide-1,ite)
6056	j_start = jts
6057	j_end = jte
6058	k_start = kts
6059	k_end = ktf
6060
6061	ELSE IF ( name .EQ. 'w' ) THEN
6062
6063	i_start = its
6064	i_end = MIN(ide-1,ite)
6065	j_start = jts
6066	j_end = MIN(jde-1,jte)
6067	k_start = kts+1
6068	k_end = ktf
6069
6070	ELSE
6071
6072	i_start = its
6073	i_end = MIN(ide-1,ite)
6074	j_start = jts
6075	j_end = MIN(jde-1,jte)
6076	k_start = kts
6077	k_end = ktf
6078
6079	ENDIF
6080
6081	! End: Assignment of limits of spatial loops.
6082	!------------------------------------------------------------------------------
6083
6084	!------------------------------------------------------------------------------
6085	! Begin: Loop across spatial dimensions.
6086
6087	DO j = j_start, j_end
6088	DO k = k_start, k_end
6089	DO i = i_start, i_end
6090
6091	!------------------------------------------------------------------------------
6092	! Begin: Diffusion in x (i index).
6093
6094	! Calculate the diffusive flux in x direction (from Xue's eq. 3).
6095
6096	dflux_x_p0 = ( 10.0 * ( field(i, k,j) - field(i-1,k,j) ) &
6097	- 5.0 * ( field(i+1,k,j) - field(i-2,k,j) ) &
6098	+ ( field(i+2,k,j) - field(i-3,k,j) ) )
6099
6100	dflux_x_p1 = ( 10.0 * ( field(i+1,k,j) - field(i ,k,j) ) &
6101	- 5.0 * ( field(i+2,k,j) - field(i-1,k,j) ) &
6102	+ ( field(i+3,k,j) - field(i-2,k,j) ) )
6103
6104	! If requested in the namelist (diff_6th_opt=2), prohibit up-gradient diffusion
6105	! (variation on Xue's eq. 10).
6106
6107	IF ( diff_6th_opt .EQ. 2 ) THEN
6108
6109	IF ( dflux_x_p0 * ( field(i ,k,j)-field(i-1,k,j) ) .LE. 0.0 ) THEN
6110	dflux_x_p0 = 0.0
6111	END IF
6112
6113	IF ( dflux_x_p1 * ( field(i+1,k,j)-field(i ,k,j) ) .LE. 0.0 ) THEN
6114	dflux_x_p1 = 0.0
6115	END IF
6116
6117	END IF
6118
6119	! Apply 6th-order diffusion in x direction.
6120
6121	IF ( name .EQ. 'u' ) THEN
6122	mu_avg_p0 = mu(i-1,j)
6123	mu_avg_p1 = mu(i ,j)
6124	ELSE IF ( name .EQ. 'v' ) THEN
6125	mu_avg_p0 = 0.25 * ( &
6126	mu(i-1,j-1) + &
6127	mu(i ,j-1) + &
6128	mu(i-1,j ) + &
6129	mu(i ,j ) )
6130	mu_avg_p1 = 0.25 * ( &
6131	mu(i ,j-1) + &
6132	mu(i+1,j-1) + &
6133	mu(i ,j ) + &
6134	mu(i+1,j ) )
6135	ELSE
6136	mu_avg_p0 = 0.5 * ( &
6137	mu(i-1,j) + &
6138	mu(i ,j) )
6139	mu_avg_p1 = 0.5 * ( &
6140	mu(i ,j) + &
6141	mu(i+1,j) )
6142	END IF
6143
6144	tendency_x = diff_6th_coef * &
6145	( ( mu_avg_p1 * dflux_x_p1 ) - ( mu_avg_p0 * dflux_x_p0 ) )
6146
6147	! End: Diffusion in x.
6148	!------------------------------------------------------------------------------
6149
6150	!------------------------------------------------------------------------------
6151	! Begin: Diffusion in y (j index).
6152
6153	! Calculate the diffusive flux in y direction (from Xue's eq. 3).
6154
6155	dflux_y_p0 = ( 10.0 * ( field(i,k,j ) - field(i,k,j-1) ) &
6156	- 5.0 * ( field(i,k,j+1) - field(i,k,j-2) ) &
6157	+ ( field(i,k,j+2) - field(i,k,j-3) ) )
6158
6159	dflux_y_p1 = ( 10.0 * ( field(i,k,j+1) - field(i,k,j ) ) &
6160	- 5.0 * ( field(i,k,j+2) - field(i,k,j-1) ) &
6161	+ ( field(i,k,j+3) - field(i,k,j-2) ) )
6162
6163	! If requested in the namelist (diff_6th_opt=2), prohibit up-gradient diffusion
6164	! (variation on Xue's eq. 10).
6165
6166	IF ( diff_6th_opt .EQ. 2 ) THEN
6167
6168	IF ( dflux_y_p0 * ( field(i,k,j )-field(i,k,j-1) ) .LE. 0.0 ) THEN
6169	dflux_y_p0 = 0.0
6170	END IF
6171
6172	IF ( dflux_y_p1 * ( field(i,k,j+1)-field(i,k,j ) ) .LE. 0.0 ) THEN
6173	dflux_y_p1 = 0.0
6174	END IF
6175
6176	END IF
6177
6178	! Apply 6th-order diffusion in y direction.
6179
6180	IF ( name .EQ. 'u' ) THEN
6181	mu_avg_p0 = 0.25 * ( &
6182	mu(i-1,j-1) + &
6183	mu(i ,j-1) + &
6184	mu(i-1,j ) + &
6185	mu(i ,j ) )
6186	mu_avg_p1 = 0.25 * ( &
6187	mu(i-1,j ) + &
6188	mu(i ,j ) + &
6189	mu(i-1,j+1) + &
6190	mu(i ,j+1) )
6191	ELSE IF ( name .EQ. 'v' ) THEN
6192	mu_avg_p0 = mu(i,j-1)
6193	mu_avg_p1 = mu(i,j )
6194	ELSE
6195	mu_avg_p0 = 0.5 * ( &
6196	mu(i,j-1) + &
6197	mu(i,j ) )
6198	mu_avg_p1 = 0.5 * ( &
6199	mu(i,j ) + &
6200	mu(i,j+1) )
6201	END IF
6202
6203	tendency_y = diff_6th_coef * &
6204	( ( mu_avg_p1 * dflux_y_p1 ) - ( mu_avg_p0 * dflux_y_p0 ) )
6205
6206	! End: Diffusion in y.
6207	!------------------------------------------------------------------------------
6208
6209	!------------------------------------------------------------------------------
6210	! Begin: Combine diffusion in x and y.
6211
6212	tendency(i,k,j) = tendency(i,k,j) + tendency_x + tendency_y
6213
6214	! End: Combine diffusion in x and y.
6215	!------------------------------------------------------------------------------
6216
6217	ENDDO
6218	ENDDO
6219	ENDDO
6220
6221	! End: Loop across spatial dimensions.
6222	!------------------------------------------------------------------------------
6223
6224	END SUBROUTINE sixth_order_diffusion
6225
6226	!==============================================================================
6227	!==============================================================================
6228
6229	END MODULE module_big_step_utilities_em

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

Download in other formats:

Original Format