Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_big_step_utilities_em.F @ 3567

Last change on this file since 3567 was 2761, checked in by aslmd, 2 years ago
Applied planetary adaptation changes to WRFV3. job done previously by LMD_LES_MARS_install. Moved Registry.EM.newphys to put it simply in Registry.EM
File size: 182.9 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
1012	itf=MIN(ite,ide-1)
1013	jtf=MIN(jte,jde-1)
1014	ktf=MIN(kte,kde-1)
1015
1016	#ifndef INTELMKL
1017	cpovcv_v = cpovcv
1018	#endif
1019
1020	IF (non_hydrostatic) THEN
1021
1022	IF (n_moist >= PARAM_FIRST_SCALAR ) THEN
1023
1024	DO j=jts,jtf
1025	DO k=kts,ktf
1026	DO i=its,itf
1027	qvf = 1.+rvovrd*moist(i,k,j,P_QV)
1028	al(i,k,j)=-1./muts(i,j)(alb(i,k,j)mu(i,j) &
1029	+rdnw(k)*(ph(i,k+1,j)-ph(i,k,j)))
1030	temp(i)=(r_d(t0+t(i,k,j))qvf)/ &
1031	(p0*(al(i,k,j)+alb(i,k,j)))
1032	ENDDO
1033	#ifdef INTELMKL
1034	CALL VPOWX ( itf-its+1, temp(its), cpovcv, p(its,k,j) )
1035	#else
1036	! use vector version from libmassv or from compat lib in frame/libmassv.F
1037	CALL VPOW ( p(its,k,j), temp(its), cpovcv_v(its), itf-its+1 )
1038	#endif
1039	DO i=its,itf
1040	p(i,k,j)= p(i,k,j)*p0-pb(i,k,j)
1041	ENDDO
1042	ENDDO
1043	ENDDO
1044
1045	ELSE
1046
1047	DO j=jts,jtf
1048	DO k=kts,ktf
1049	DO i=its,itf
1050	al(i,k,j)=-1./muts(i,j)(alb(i,k,j)mu(i,j) &
1051	+rdnw(k)*(ph(i,k+1,j)-ph(i,k,j)))
1052	p(i,k,j)=p0( (r_d(t0+t(i,k,j)))/ &
1053	(p0(al(i,k,j)+alb(i,k,j))) )*cpovcv &
1054	-pb(i,k,j)
1055	ENDDO
1056	ENDDO
1057	ENDDO
1058
1059	END IF
1060
1061	ELSE
1062
1063	! hydrostatic pressure, al, and ph1 calc; WCS, 5 sept 2001
1064
1065
1066	IF (n_moist >= PARAM_FIRST_SCALAR ) THEN
1067
1068	DO j=jts,jtf
1069
1070	k=ktf ! top layer
1071	DO i=its,itf
1072
1073	qtot = 0.
1074	DO ispe=PARAM_FIRST_SCALAR,n_moist
1075	qtot = qtot + moist(i,k,j,ispe)
1076	ENDDO
1077	qf2 = 1./(1.+qtot)
1078	qf1 = qtot*qf2
1079
1080	p(i,k,j) = - 0.5(mu(i,j)+qf1muts(i,j))/rdnw(k)/qf2
1081	qvf = 1.+rvovrd*moist(i,k,j,P_QV)
1082	al(i,k,j) = (r_d/p1000mb)(t(i,k,j)+t0)qvf* &
1083	(((p(i,k,j)+pb(i,k,j))/p1000mb)**cvpm) - alb(i,k,j)
1084
1085	ENDDO
1086
1087	DO k=ktf-1,kts,-1 ! remaining layers, integrate down
1088	DO i=its,itf
1089
1090	qtot = 0.
1091	DO ispe=PARAM_FIRST_SCALAR,n_moist
1092	qtot = qtot + 0.5*( moist(i,k ,j,ispe) + moist(i,k+1,j,ispe) )
1093	ENDDO
1094	qf2 = 1./(1.+qtot)
1095	qf1 = qtot*qf2
1096
1097	p(i,k,j) = p(i,k+1,j) - (mu(i,j) + qf1*muts(i,j))/qf2/rdn(k+1)
1098	qvf = 1.+rvovrd*moist(i,k,j,P_QV)
1099	al(i,k,j) = (r_d/p1000mb)(t(i,k,j)+t0)qvf* &
1100	(((p(i,k,j)+pb(i,k,j))/p1000mb)**cvpm) - alb(i,k,j)
1101	ENDDO
1102	ENDDO
1103
1104	DO k=2,ktf+1 ! integrate hydrostatic equation for geopotential
1105	DO i=its,itf
1106
1107	! ph(i,k,j) = ph(i,k-1,j) - (1./rdnw(k-1))*( &
1108	! (muts(i,j)+mu(i,j))*al(i,k-1,j)+ &
1109	! mu(i,j)*alb(i,k-1,j) )
1110	ph(i,k,j) = ph(i,k-1,j) - (dnw(k-1))*( &
1111	(muts(i,j))*al(i,k-1,j)+ &
1112	mu(i,j)*alb(i,k-1,j) )
1113
1114
1115	ENDDO
1116	ENDDO
1117
1118	ENDDO
1119
1120	ELSE
1121
1122	DO j=jts,jtf
1123
1124	k=ktf ! top layer
1125	DO i=its,itf
1126
1127	qtot = 0.
1128	qf2 = 1./(1.+qtot)
1129	qf1 = qtot*qf2
1130
1131	p(i,k,j) = - 0.5(mu(i,j)+qf1muts(i,j))/rdnw(k)/qf2
1132	qvf = 1.
1133	al(i,k,j) = (r_d/p1000mb)(t(i,k,j)+t0)qvf* &
1134	(((p(i,k,j)+pb(i,k,j))/p1000mb)**cvpm) - alb(i,k,j)
1135
1136	ENDDO
1137
1138	DO k=ktf-1,kts,-1 ! remaining layers, integrate down
1139	DO i=its,itf
1140
1141	qtot = 0.
1142	qf2 = 1./(1.+qtot)
1143	qf1 = qtot*qf2
1144
1145	p(i,k,j) = p(i,k+1,j) - (mu(i,j) + qf1*muts(i,j))/qf2/rdn(k+1)
1146	qvf = 1.
1147	al(i,k,j) = (r_d/p1000mb)(t(i,k,j)+t0)qvf* &
1148	(((p(i,k,j)+pb(i,k,j))/p1000mb)**cvpm) - alb(i,k,j)
1149	ENDDO
1150	ENDDO
1151
1152	DO k=2,ktf+1 ! integrate hydrostatic equation for geopotential
1153	DO i=its,itf
1154
1155	! ph(i,k,j) = ph(i,k-1,j) - (1./rdnw(k-1))*( &
1156	! (muts(i,j)+mu(i,j))*al(i,k-1,j)+ &
1157	! mu(i,j)*alb(i,k-1,j) )
1158	ph(i,k,j) = ph(i,k-1,j) - (dnw(k-1))*( &
1159	(muts(i,j))*al(i,k-1,j)+ &
1160	mu(i,j)*alb(i,k-1,j) )
1161
1162
1163	ENDDO
1164	ENDDO
1165
1166	ENDDO
1167
1168	END IF
1169
1170	END IF
1171
1172	END SUBROUTINE calc_p_rho_phi
1173
1174	!----------------------------------------------------------------------
1175
1176	SUBROUTINE calc_php ( php, ph, phb, &
1177	ids, ide, jds, jde, kds, kde, &
1178	ims, ime, jms, jme, kms, kme, &
1179	its, ite, jts, jte, kts, kte )
1180
1181	IMPLICIT NONE
1182
1183	! Input data
1184
1185	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
1186	ims, ime, jms, jme, kms, kme, &
1187	its, ite, jts, jte, kts, kte
1188
1189	REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(IN ) :: phb, ph
1190	REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT( OUT) :: php
1191
1192	! Local stuff
1193
1194	INTEGER :: i, j, k, itf, jtf, ktf
1195
1196	!<DESCRIPTION>
1197	!
1198	! calc_php calculates the full geopotential from the reference state
1199	! geopotential and the perturbation geopotential (phb_ph).
1200	!
1201	!</DESCRIPTION>
1202
1203	itf=MIN(ite,ide-1)
1204	jtf=MIN(jte,jde-1)
1205	ktf=MIN(kte,kde-1)
1206
1207	DO j=jts,jtf
1208	DO k=kts,ktf
1209	DO i=its,itf
1210	php(i,k,j) = 0.5*(phb(i,k,j)+phb(i,k+1,j)+ph(i,k,j)+ph(i,k+1,j))
1211	ENDDO
1212	ENDDO
1213	ENDDO
1214
1215	END SUBROUTINE calc_php
1216
1217	!-------------------------------------------------------------------------------
1218
1219	SUBROUTINE diagnose_w( ph_tend, ph_new, ph_old, w, mu, dt, &
1220	u, v, ht, &
1221	cf1, cf2, cf3, rdx, rdy, &
1222	msftx, msfty, &
1223	ids, ide, jds, jde, kds, kde, &
1224	ims, ime, jms, jme, kms, kme, &
1225	its, ite, jts, jte, kts, kte )
1226
1227	IMPLICIT NONE
1228
1229	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
1230	ims, ime, jms, jme, kms, kme, &
1231	its, ite, jts, jte, kts, kte
1232
1233	REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(IN ) :: ph_tend, &
1234	ph_new, &
1235	ph_old, &
1236	u, &
1237	v
1238
1239
1240	REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT( OUT) :: w
1241
1242	REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN ) :: mu, ht, msftx, msfty
1243
1244	REAL, INTENT(IN ) :: dt, cf1, cf2, cf3, rdx, rdy
1245
1246	INTEGER :: i, j, k, itf, jtf
1247
1248	itf=MIN(ite,ide-1)
1249	jtf=MIN(jte,jde-1)
1250
1251	!<DESCRIPTION>
1252	!
1253	! diagnose_w diagnoses the vertical velocity from the geopoential equation.
1254	! Used with the hydrostatic option.
1255	!
1256	!</DESCRIPTION>
1257
1258	DO j = jts, jtf
1259
1260	! lower b.c. on w
1261
1262	! Notes on map scale factors:
1263	! Chain rule: if Z=Z(X,Y) [true at the surface] then
1264	! dZ/dt = dZ/dX * dX/dt + dZ/dY * dY/dt, U=dX/dt, V=dY/dt
1265	! Using capitals to denote actual values
1266	! In mapped values, u=U, v=V, z=Z, 1/dX=mx/dx, 1/dY=my/dy
1267	! => w = dz/dt = mx u dz/dx + my v dz/dy
1268	! [where dz/dx is just the surface height change between x
1269	! gridpoints, and dz/dy is the change between y gridpoints]
1270	! [NB: cf1, cf2 and cf3 do vertical weighting of u or v values
1271	! nearest the surface]
1272
1273	! Previously msft multiplied by rdy and rdx terms.
1274	! Now msfty multiplies rdy term, and msftx multiplies msftx term
1275	DO i = its, itf
1276	w(i,1,j)= msfty(i,j).5rdy*( &
1277	(ht(i,j+1)-ht(i,j )) &
1278	(cf1v(i,1,j+1)+cf2v(i,2,j+1)+cf3v(i,3,j+1)) &
1279	+(ht(i,j )-ht(i,j-1)) &
1280	(cf1v(i,1,j )+cf2v(i,2,j )+cf3v(i,3,j )) ) &
1281	+msftx(i,j).5rdx*( &
1282	(ht(i+1,j)-ht(i,j )) &
1283	(cf1u(i+1,1,j)+cf2u(i+1,2,j)+cf3u(i+1,3,j)) &
1284	+(ht(i,j )-ht(i-1,j)) &
1285	(cf1u(i ,1,j)+cf2u(i ,2,j)+cf3u(i ,3,j)) )
1286	ENDDO
1287
1288	! use geopotential equation to diagnose w
1289
1290	! Further notes on map scale factors
1291	! If ph_tend contains: -mx partial d/dx(mu rho u/my)
1292	! -mx partial d/dy(phi mu v/mx)
1293	! -partial d/dz(phi mu w/my)
1294	! then phi eqn is: partial d/dt(mu phi/my) = ph_tend + mu g w/my
1295	! => w = [my/(mug)][partial d/dt(mu phi/my) - ph_tend]
1296
1297	DO k = 2, kte
1298	DO i = its, itf
1299	w(i,k,j) = msfty(i,j)*( (ph_new(i,k,j)-ph_old(i,k,j))/dt &
1300	- ph_tend(i,k,j)/mu(i,j) )/g
1301
1302	ENDDO
1303	ENDDO
1304
1305	ENDDO
1306
1307	END SUBROUTINE diagnose_w
1308
1309	!-------------------------------------------------------------------------------
1310
1311	SUBROUTINE rhs_ph( ph_tend, u, v, ww, &
1312	ph, ph_old, phb, w, &
1313	mut, muu, muv, &
1314	fnm, fnp, &
1315	rdnw, cfn, cfn1, rdx, rdy, &
1316	msfux, msfuy, msfvx, &
1317	msfvx_inv, msfvy, &
1318	msftx, msfty, &
1319	non_hydrostatic, &
1320	config_flags, &
1321	ids, ide, jds, jde, kds, kde, &
1322	ims, ime, jms, jme, kms, kme, &
1323	its, ite, jts, jte, kts, kte )
1324	IMPLICIT NONE
1325
1326	TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
1327
1328	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
1329	ims, ime, jms, jme, kms, kme, &
1330	its, ite, jts, jte, kts, kte
1331
1332	REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(IN ) :: &
1333	u, &
1334	v, &
1335	ww, &
1336	ph, &
1337	ph_old, &
1338	phb, &
1339	w
1340
1341	REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(INOUT) :: ph_tend
1342
1343	REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN ) :: muu, muv, mut, &
1344	msfux, msfuy, &
1345	msfvx, msfvy, &
1346	msftx, msfty, &
1347	msfvx_inv
1348
1349	REAL, DIMENSION( kms:kme ), INTENT(IN ) :: rdnw, fnm, fnp
1350
1351	REAL, INTENT(IN ) :: cfn, cfn1, rdx, rdy
1352
1353	LOGICAL, INTENT(IN ) :: non_hydrostatic
1354
1355	! Local stuff
1356
1357	INTEGER :: i, j, k, itf, jtf, ktf, kz, i_start, j_start
1358	REAL :: ur, ul, ub, vr, vl, vb
1359	REAL, DIMENSION(its:ite,kts:kte) :: wdwn
1360
1361	INTEGER :: advective_order
1362
1363	LOGICAL :: specified
1364
1365	!<DESCRIPTION>
1366	!
1367	! rhs_ph calculates the large-timestep tendency terms for the geopotential
1368	! equation. These terms include the advection and "gw". The geopotential
1369	! equation is cast in advective form, so we don't use the flux form advection
1370	! algorithms here.
1371	!
1372	!</DESCRIPTION>
1373
1374	specified = .false.
1375	if(config_flags%specified .or. config_flags%nested) specified = .true.
1376
1377	advective_order = config_flags%h_sca_adv_order
1378	! advective_order = 2 ! original configuration (pre Oct 2001)
1379
1380	itf=MIN(ite,ide-1)
1381	jtf=MIN(jte,jde-1)
1382	ktf=MIN(kte,kde-1)
1383
1384	! Notes on map scale factors (WCS, 2 march 2008)
1385	! phi equation is: mu/my d/dt(phi) = -(1/my) mx mu u d/dx(phi)
1386	! -(1/my) my mu v d/dy(phi)
1387	! - omega d/d_eta(phi)
1388	! +mu g w/my
1389	!
1390	! A little further explanation...
1391	! The tendency term we are computing here is for mu/my d/dt(phi). It is advective form
1392	! but it is multiplied be mu/my. It will be decoupled from (mu/my) when the implicit w-phi
1393	! solution is computed in subourine advance_w. The formulation dates from the early
1394	! days of the mass coordinate model when we were testing both a flux and an advective formulation
1395	! for the geopotential equation and different forms of the vertical momentum equation and the
1396	! vertically implicit solver.
1397
1398	! advective form for the geopotential equation
1399
1400	DO j = jts, jtf
1401
1402	DO k = 2, kte
1403	DO i = its, itf
1404	wdwn(i,k) = .5(ww(i,k,j)+ww(i,k-1,j))rdnw(k-1) &
1405	*(ph(i,k,j)-ph(i,k-1,j)+phb(i,k,j)-phb(i,k-1,j))
1406	ENDDO
1407	ENDDO
1408
1409	! RHS term 3 is: - omega partial d/dnu(phi)
1410
1411	DO k = 2, kte-1
1412	DO i = its, itf
1413	ph_tend(i,k,j) = ph_tend(i,k,j) &
1414	- (fnm(k)wdwn(i,k+1)+fnp(k)wdwn(i,k))
1415	ENDDO
1416	ENDDO
1417
1418	ENDDO
1419
1420	IF (non_hydrostatic) THEN ! add in "gw" term.
1421	DO j = jts, jtf ! in hydrostatic mode, "gw" will be diagnosed
1422	! after the timestep to give us "w"
1423	DO i = its, itf
1424	ph_tend(i,kde,j) = 0.
1425	ENDDO
1426
1427	DO k = 2, kte
1428	DO i = its, itf
1429	! phi equation RHS term 4
1430	ph_tend(i,k,j) = ph_tend(i,k,j) + mut(i,j)gw(i,k,j)/msfty(i,j)
1431	ENDDO
1432	ENDDO
1433
1434	ENDDO
1435
1436	END IF
1437
1438	! Notes on map scale factors:
1439	! RHS terms 1 and 2 are: -(1/my) mx u mu partial d/dx(phi)
1440	! -(1/my) my v mu partial d/dy(phi)
1441
1442	IF (advective_order <= 2) THEN
1443
1444	! y (v) advection
1445
1446	i_start = its
1447	j_start = jts
1448	itf=MIN(ite,ide-1)
1449	jtf=MIN(jte,jde-1)
1450
1451	IF ( (config_flags%open_ys) .and. jts == jds ) j_start = jts+1
1452	IF ( (config_flags%open_ye) .and. jte == jde ) jtf = jtf-1
1453
1454	DO j = j_start, jtf
1455
1456	DO k = 2, kte-1
1457	DO i = i_start, itf
1458	ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25rdy/msfty(i,j)) &
1459	( muv(i,j+1)(v(i,k,j+1)+v(i,k-1,j+1))msfvy(i,j+1)* &
1460	(phb(i,k,j+1)-phb(i,k,j )+ph(i,k,j+1)-ph(i,k,j )) &
1461	+muv(i,j )(v(i,k,j )+v(i,k-1,j ))msfvy(i,j )* &
1462	(phb(i,k,j )-phb(i,k,j-1)+ph(i,k,j )-ph(i,k,j-1)) )
1463	ENDDO
1464	ENDDO
1465
1466	k = kte
1467	DO i = i_start, itf
1468	ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5rdy/msfty(i,j)) &
1469	( muv(i,j+1)(cfnv(i,k-1,j+1)+cfn1v(i,k-2,j+1))msfvy(i,j+1)* &
1470	(phb(i,k,j+1)-phb(i,k,j )+ph(i,k,j+1)-ph(i,k,j )) &
1471	+muv(i,j )(cfnv(i,k-1,j )+cfn1v(i,k-2,j ))msfvy(i,j )* &
1472	(phb(i,k,j )-phb(i,k,j-1)+ph(i,k,j )-ph(i,k,j-1)) )
1473	ENDDO
1474
1475	ENDDO
1476
1477	! x (u) advection
1478
1479	i_start = its
1480	j_start = jts
1481	itf=MIN(ite,ide-1)
1482	jtf=MIN(jte,jde-1)
1483
1484	IF ( (config_flags%open_xs) .and. its == ids ) i_start = its+1
1485	IF ( (config_flags%open_xe) .and. ite == ide ) itf = itf-1
1486
1487	DO j = j_start, jtf
1488
1489	DO k = 2, kte-1
1490	DO i = i_start, itf
1491	ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25rdx/msfty(i,j)) &
1492	( muu(i+1,j)(u(i+1,k,j)+u(i+1,k-1,j))msfux(i+1,j)* &
1493	(phb(i+1,k,j)-phb(i ,k,j)+ph(i+1,k,j)-ph(i ,k,j)) &
1494	+muu(i ,j)(u(i ,k,j)+u(i ,k-1,j))msfux(i ,j)* &
1495	(phb(i ,k,j)-phb(i-1,k,j)+ph(i ,k,j)-ph(i-1,k,j)) )
1496	ENDDO
1497	ENDDO
1498
1499	k = kte
1500	DO i = i_start, itf
1501	ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5rdx/msfty(i,j)) &
1502	( muu(i+1,j)(cfnu(i+1,k-1,j)+cfn1u(i+1,k-2,j))msfux(i+1,j)* &
1503	(phb(i+1,k,j)-phb(i ,k,j)+ph(i+1,k,j)-ph(i ,k,j)) &
1504	+muu(i ,j)(cfnu(i ,k-1,j)+cfn1u(i ,k-2,j))msfux( i,j)* &
1505	(phb(i ,k,j)-phb(i-1,k,j)+ph(i ,k,j)-ph(i-1,k,j)) )
1506	ENDDO
1507
1508	ENDDO
1509
1510	ELSE IF (advective_order <= 4) THEN
1511
1512	! y (v) advection
1513
1514	i_start = its
1515	j_start = jts
1516	itf=MIN(ite,ide-1)
1517	jtf=MIN(jte,jde-1)
1518
1519	IF ( (config_flags%open_ys) .and. jts == jds ) j_start = jts+1
1520	IF ( (config_flags%open_ye) .and. jte == jde ) jtf = jtf-1
1521
1522	DO j = j_start, jtf
1523
1524	DO k = 2, kte-1
1525	DO i = i_start, itf
1526	ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25rdy/msfty(i,j))( &
1527	( muv(i,j+1)(v(i,k,j+1)+v(i,k-1,j+1))msfvy(i,j+1) &
1528	+muv(i,j )(v(i,k,j )+v(i,k-1,j ))msfvy(i,j ))* (1./12.)*( &
1529	8.*(ph(i,k,j+1)-ph(i,k,j-1)) &
1530	-(ph(i,k,j+2)-ph(i,k,j-2)) &
1531	+8.*(phb(i,k,j+1)-phb(i,k,j-1)) &
1532	-(phb(i,k,j+2)-phb(i,k,j-2)) ) )
1533
1534
1535	ENDDO
1536	ENDDO
1537
1538	k = kte
1539	DO i = i_start, itf
1540	ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5rdy/msfty(i,j))( &
1541	( muv(i,j+1)(cfnv(i,k-1,j+1)+cfn1v(i,k-2,j+1))msfvy(i,j+1) &
1542	+muv(i,j )(cfnv(i,k-1,j )+cfn1v(i,k-2,j ))msfvy(i,j ))* (1./12.)*( &
1543	8.*(ph(i,k,j+1)-ph(i,k,j-1)) &
1544	-(ph(i,k,j+2)-ph(i,k,j-2)) &
1545	+8.*(phb(i,k,j+1)-phb(i,k,j-1)) &
1546	-(phb(i,k,j+2)-phb(i,k,j-2)) ) )
1547
1548	ENDDO
1549
1550	ENDDO
1551
1552
1553	! x (u) advection
1554
1555	i_start = its
1556	j_start = jts
1557	itf=MIN(ite,ide-1)
1558	jtf=MIN(jte,jde-1)
1559
1560	IF ( (config_flags%open_xs) .and. its == ids ) i_start = its+1
1561	IF ( (config_flags%open_xe) .and. ite == ide ) itf = itf-1
1562
1563	DO j = j_start, jtf
1564
1565	DO k = 2, kte-1
1566	DO i = i_start, itf
1567	ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25rdx/msfty(i,j))( &
1568	( muu(i+1,j)(u(i+1,k,j)+u(i+1,k-1,j))msfux(i+1,j) &
1569	+muu(i,j )(u(i ,k,j)+u(i ,k-1,j))msfux(i ,j) )* (1./12.)*( &
1570	8.*(ph(i+1,k,j)-ph(i-1,k,j)) &
1571	-(ph(i+2,k,j)-ph(i-2,k,j)) &
1572	+8.*(phb(i+1,k,j)-phb(i-1,k,j)) &
1573	-(phb(i+2,k,j)-phb(i-2,k,j)) ) )
1574	ENDDO
1575	ENDDO
1576
1577	k = kte
1578	DO i = i_start, itf
1579	ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5rdx/msfty(i,j))( &
1580	( muu(i+1,j)(cfnu(i+1,k-1,j)+cfn1u(i+1,k-2,j))msfux(i+1,j) &
1581	+muu(i,j )(cfnu(i ,k-1,j)+cfn1u(i ,k-2,j))msfux(i ,j) )* (1./12.)*( &
1582	8.*(ph(i+1,k,j)-ph(i-1,k,j)) &
1583	-(ph(i+2,k,j)-ph(i-2,k,j)) &
1584	+8.*(phb(i+1,k,j)-phb(i-1,k,j)) &
1585	-(phb(i+2,k,j)-phb(i-2,k,j)) ) )
1586	ENDDO
1587
1588	ENDDO
1589
1590	ELSE IF (advective_order <= 6) THEN
1591
1592	! y (v) advection
1593
1594	i_start = its
1595	j_start = jts
1596	itf=MIN(ite,ide-1)
1597	jtf=MIN(jte,jde-1)
1598
1599	! IF ( (config_flags%open_ys) .and. jts == jds ) j_start = jts+1
1600	! IF ( (config_flags%open_ye) .and. jte == jde ) jtf = jtf-1
1601
1602	IF (config_flags%open_ys .or. specified ) j_start = max(jts,jds+2)
1603	IF (config_flags%open_ye .or. specified ) jtf = min(jtf,jde-3)
1604
1605	DO j = j_start, jtf
1606
1607	DO k = 2, kte-1
1608	DO i = i_start, itf
1609	ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25rdy/msfty(i,j)) ( &
1610	( muv(i,j+1)(v(i,k,j+1)+v(i,k-1,j+1))msfvy(i,j+1) &
1611	+muv(i,j )(v(i,k,j )+v(i,k-1,j ))msfvy(i,j ) )* (1./60.)*( &
1612	45.*(ph(i,k,j+1)-ph(i,k,j-1)) &
1613	-9.*(ph(i,k,j+2)-ph(i,k,j-2)) &
1614	+(ph(i,k,j+3)-ph(i,k,j-3)) &
1615	+45.*(phb(i,k,j+1)-phb(i,k,j-1)) &
1616	-9.*(phb(i,k,j+2)-phb(i,k,j-2)) &
1617	+(phb(i,k,j+3)-phb(i,k,j-3)) ) )
1618
1619
1620	ENDDO
1621	ENDDO
1622
1623	k = kte
1624	DO i = i_start, itf
1625	ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5rdy/msfty(i,j)) ( &
1626	( muv(i,j+1)(cfnv(i,k-1,j+1)+cfn1v(i,k-2,j+1))msfvy(i,j+1) &
1627	+muv(i,j )(cfnv(i,k-1,j )+cfn1v(i,k-2,j ))msfvy(i,j ) )* (1./60.)*( &
1628	45.*(ph(i,k,j+1)-ph(i,k,j-1)) &
1629	-9.*(ph(i,k,j+2)-ph(i,k,j-2)) &
1630	+(ph(i,k,j+3)-ph(i,k,j-3)) &
1631	+45.*(phb(i,k,j+1)-phb(i,k,j-1)) &
1632	-9.*(phb(i,k,j+2)-phb(i,k,j-2)) &
1633	+(phb(i,k,j+3)-phb(i,k,j-3)) ) )
1634
1635	ENDDO
1636
1637	ENDDO
1638
1639
1640	! pick up near boundary rows using 4th order stencil
1641	! (open bc copy only goes out to jds-1 and jde, hence 4rth is ok but 6th is too big)
1642
1643	IF ( (config_flags%open_ys) .and. jts <= jds+1 ) THEN
1644
1645	j = jds+1
1646	DO k = 2, kte-1
1647	DO i = i_start, itf
1648	ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25rdy/msfty(i,j)) ( &
1649	( muv(i,j+1)(v(i,k,j+1)+v(i,k-1,j+1))msfvy(i,j+1) &
1650	+muv(i,j )(v(i,k,j )+v(i,k-1,j ))msfvy(i,j ) )* (1./12.)*( &
1651	8.*(ph(i,k,j+1)-ph(i,k,j-1)) &
1652	-(ph(i,k,j+2)-ph(i,k,j-2)) &
1653	+8.*(phb(i,k,j+1)-phb(i,k,j-1)) &
1654	-(phb(i,k,j+2)-phb(i,k,j-2)) ) )
1655
1656
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	+muv(i,j )(cfnv(i,k-1,j )+cfn1v(i,k-2,j ))msfvy(i,j) )* (1./12.)*( &
1665	8.*(ph(i,k,j+1)-ph(i,k,j-1)) &
1666	-(ph(i,k,j+2)-ph(i,k,j-2)) &
1667	+8.*(phb(i,k,j+1)-phb(i,k,j-1)) &
1668	-(phb(i,k,j+2)-phb(i,k,j-2)) ) )
1669
1670	ENDDO
1671
1672	END IF
1673
1674	IF ( (config_flags%open_ye) .and. jte >= jde-2 ) THEN
1675
1676	j = jde-2
1677	DO k = 2, kte-1
1678	DO i = i_start, itf
1679	ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25rdy/msfty(i,j)) ( &
1680	( muv(i,j+1)(v(i,k,j+1)+v(i,k-1,j+1))msfvy(i,j+1) &
1681	+muv(i,j )(v(i,k,j )+v(i,k-1,j ))msfvy(i,j) )* (1./12.)*( &
1682	8.*(ph(i,k,j+1)-ph(i,k,j-1)) &
1683	-(ph(i,k,j+2)-ph(i,k,j-2)) &
1684	+8.*(phb(i,k,j+1)-phb(i,k,j-1)) &
1685	-(phb(i,k,j+2)-phb(i,k,j-2)) ) )
1686
1687
1688	ENDDO
1689	ENDDO
1690
1691	k = kte
1692	DO i = i_start, itf
1693	ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5rdy/msfty(i,j)) ( &
1694	( muv(i,j+1)(cfnv(i,k-1,j+1)+cfn1v(i,k-2,j+1))msfvy(i,j+1) &
1695	+muv(i,j )(cfnv(i,k-1,j )+cfn1v(i,k-2,j ))msfvy(i,j) )* (1./12.)*( &
1696	8.*(ph(i,k,j+1)-ph(i,k,j-1)) &
1697	-(ph(i,k,j+2)-ph(i,k,j-2)) &
1698	+8.*(phb(i,k,j+1)-phb(i,k,j-1)) &
1699	-(phb(i,k,j+2)-phb(i,k,j-2)) ) )
1700
1701	ENDDO
1702
1703	END IF
1704
1705	! x (u) advection
1706
1707	i_start = its
1708	j_start = jts
1709	itf=MIN(ite,ide-1)
1710	jtf=MIN(jte,jde-1)
1711
1712	IF (config_flags%open_xs .or. specified ) i_start = max(its,ids+2)
1713	IF (config_flags%open_xe .or. specified ) itf = min(itf,ide-3)
1714	IF ( config_flags%periodic_x ) i_start = its
1715	IF ( config_flags%periodic_x ) itf=MIN(ite,ide-1)
1716
1717	DO j = j_start, jtf
1718
1719	DO k = 2, kte-1
1720	DO i = i_start, itf
1721	ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25rdx/msfty(i,j))( &
1722	( muu(i+1,j)(u(i+1,k,j)+u(i+1,k-1,j))msfux(i+1,j) &
1723	+muu(i,j )(u(i,k,j )+u(i,k-1,j ))msfux(i,j) )* (1./60.)*( &
1724	45.*(ph(i+1,k,j)-ph(i-1,k,j)) &
1725	-9.*(ph(i+2,k,j)-ph(i-2,k,j)) &
1726	+(ph(i+3,k,j)-ph(i-3,k,j)) &
1727	+45.*(phb(i+1,k,j)-phb(i-1,k,j)) &
1728	-9.*(phb(i+2,k,j)-phb(i-2,k,j)) &
1729	+(phb(i+3,k,j)-phb(i-3,k,j)) ) )
1730	ENDDO
1731	ENDDO
1732
1733	k = kte
1734	DO i = i_start, itf
1735	ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5rdx/msfty(i,j))( &
1736	( muu(i+1,j)(cfnu(i+1,k-1,j)+cfn1u(i+1,k-2,j))msfux(i+1,j) &
1737	+muu(i,j )(cfnu(i ,k-1,j)+cfn1u(i,k-2,j))msfux(i,j) )* (1./60.)*( &
1738	45.*(ph(i+1,k,j)-ph(i-1,k,j)) &
1739	-9.*(ph(i+2,k,j)-ph(i-2,k,j)) &
1740	+(ph(i+3,k,j)-ph(i-3,k,j)) &
1741	+45.*(phb(i+1,k,j)-phb(i-1,k,j)) &
1742	-9.*(phb(i+2,k,j)-phb(i-2,k,j)) &
1743	+(phb(i+3,k,j)-phb(i-3,k,j)) ) )
1744	ENDDO
1745
1746	ENDDO
1747
1748	IF ( (config_flags%open_xs) .and. its <= ids+1 ) THEN
1749	i = ids + 1
1750	DO j = j_start, jtf
1751	DO k = 2, kte-1
1752	ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25rdx/msfty(i,j))( &
1753	( muu(i+1,j)(u(i+1,k,j)+u(i+1,k-1,j))msfux(i+1,j) &
1754	+muu(i,j )(u(i,k,j )+u(i,k-1,j ))msfux(i,j) )* (1./12.)*( &
1755	8.*(ph(i+1,k,j)-ph(i-1,k,j)) &
1756	-(ph(i+2,k,j)-ph(i-2,k,j)) &
1757	+8.*(phb(i+1,k,j)-phb(i-1,k,j)) &
1758	-(phb(i+2,k,j)-phb(i-2,k,j)) ) )
1759	ENDDO
1760	k = kte
1761	ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5rdx/msfty(i,j))( &
1762	( muu(i+1,j)(cfnu(i+1,k-1,j)+cfn1u(i+1,k-2,j))msfux(i+1,j) &
1763	+muu(i,j )(cfnu(i ,k-1,j)+cfn1u(i,k-2,j))msfux(i,j) )* (1./12.)*( &
1764	8.*(ph(i+1,k,j)-ph(i-1,k,j)) &
1765	-(ph(i+2,k,j)-ph(i-2,k,j)) &
1766	+8.*(phb(i+1,k,j)-phb(i-1,k,j)) &
1767	-(phb(i+2,k,j)-phb(i-2,k,j)) ) )
1768
1769	ENDDO
1770	END IF
1771
1772	IF ( (config_flags%open_xe) .and. ite >= ide-2 ) THEN
1773	i = ide-2
1774	DO j = j_start, jtf
1775	DO k = 2, kte-1
1776	ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25rdx/msfty(i,j))( &
1777	( muu(i+1,j)(u(i+1,k,j)+u(i+1,k-1,j))msfux(i+1,j) &
1778	+muu(i,j )(u(i,k,j )+u(i,k-1,j ))msfux(i,j) )* (1./12.)*( &
1779	8.*(ph(i+1,k,j)-ph(i-1,k,j)) &
1780	-(ph(i+2,k,j)-ph(i-2,k,j)) &
1781	+8.*(phb(i+1,k,j)-phb(i-1,k,j)) &
1782	-(phb(i+2,k,j)-phb(i-2,k,j)) ) )
1783	ENDDO
1784	k = kte
1785	ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5rdx/msfty(i,j))( &
1786	( muu(i+1,j)(cfnu(i+1,k-1,j)+cfn1u(i+1,k-2,j))msfux(i+1,j) &
1787	+muu(i,j )(cfnu(i ,k-1,j)+cfn1u(i,k-2,j))msfux(i,j) )* (1./12.)*( &
1788	8.*(ph(i+1,k,j)-ph(i-1,k,j)) &
1789	-(ph(i+2,k,j)-ph(i-2,k,j)) &
1790	+8.*(phb(i+1,k,j)-phb(i-1,k,j)) &
1791	-(phb(i+2,k,j)-phb(i-2,k,j)) ) )
1792
1793	ENDDO
1794	END IF
1795
1796	END IF
1797
1798	! lateral open boundary conditions,
1799	! start with north and south (y) boundaries
1800
1801	i_start = its
1802	itf=MIN(ite,ide-1)
1803
1804	! south
1805
1806	IF ( (config_flags%open_ys) .and. jts == jds ) THEN
1807
1808	j=jts
1809
1810	DO k=2,kde
1811	kz = min(k,kde-1)
1812	DO i = its,itf
1813	vb =.5( fnm(kz)(v(i,kz ,j+1)+v(i,kz ,j )) &
1814	+fnp(kz)*(v(i,kz-1,j+1)+v(i,kz-1,j )) )
1815	vl=amin1(vb,0.)
1816	ph_tend(i,k,j)=ph_tend(i,k,j)-rdymut(i,j)( &
1817	+vl*(ph_old(i,k,j+1)-ph_old(i,k,j)))
1818	ENDDO
1819	ENDDO
1820
1821	END IF
1822
1823	! north
1824
1825	IF ( (config_flags%open_ye) .and. jte == jde ) THEN
1826
1827	j=jte-1
1828
1829	DO k=2,kde
1830	kz = min(k,kde-1)
1831	DO i = its,itf
1832	vb=.5( fnm(kz)(v(i,kz ,j+1)+v(i,kz ,j)) &
1833	+fnp(kz)*(v(i,kz-1,j+1)+v(i,kz-1,j)) )
1834	vr=amax1(vb,0.)
1835	ph_tend(i,k,j)=ph_tend(i,k,j)-rdymut(i,j)( &
1836	+vr*(ph_old(i,k,j)-ph_old(i,k,j-1)))
1837	ENDDO
1838	ENDDO
1839
1840	END IF
1841
1842	! now the east and west (y) boundaries
1843
1844	j_start = its
1845	jtf=MIN(jte,jde-1)
1846
1847	! west
1848
1849	IF ( (config_flags%open_xs) .and. its == ids ) THEN
1850
1851	i=its
1852
1853	DO j = jts,jtf
1854	DO k=2,kde-1
1855	kz = k
1856	ub =.5( fnm(kz)(u(i+1,kz ,j)+u(i ,kz ,j)) &
1857	+fnp(kz)*(u(i+1,kz-1,j)+u(i ,kz-1,j)) )
1858	ul=amin1(ub,0.)
1859	ph_tend(i,k,j)=ph_tend(i,k,j)-(msftx(i,j)/msfty(i,j))rdxmut(i,j)*( &
1860	+ul*(ph_old(i+1,k,j)-ph_old(i,k,j)))
1861	ENDDO
1862
1863	k = kde
1864	kz = k
1865	ub =.5( fnm(kz)(u(i+1,kz ,j)+u(i ,kz ,j)) &
1866	+fnp(kz)*(u(i+1,kz-1,j)+u(i ,kz-1,j)) )
1867	ul=amin1(ub,0.)
1868	ph_tend(i,k,j)=ph_tend(i,k,j)-(msftx(i,j)/msfty(i,j))rdxmut(i,j)*( &
1869	+ul*(ph_old(i+1,k,j)-ph_old(i,k,j)))
1870	ENDDO
1871
1872	END IF
1873
1874	! east
1875
1876	IF ( (config_flags%open_xe) .and. ite == ide ) THEN
1877
1878	i = ite-1
1879
1880	DO j = jts,jtf
1881	DO k=2,kde-1
1882	kz = k
1883	ub=.5( fnm(kz)(u(i+1,kz ,j)+u(i,kz ,j)) &
1884	+fnp(kz)*(u(i+1,kz-1,j)+u(i,kz-1,j)) )
1885	ur=amax1(ub,0.)
1886	ph_tend(i,k,j)=ph_tend(i,k,j)-(msftx(i,j)/msfty(i,j))rdxmut(i,j)*( &
1887	+ur*(ph_old(i,k,j)-ph_old(i-1,k,j)))
1888	ENDDO
1889
1890	k = kde
1891	kz = k-1
1892	ub=.5( fnm(kz)(u(i+1,kz ,j)+u(i,kz ,j)) &
1893	+fnp(kz)*(u(i+1,kz-1,j)+u(i,kz-1,j)) )
1894	ur=amax1(ub,0.)
1895	ph_tend(i,k,j)=ph_tend(i,k,j)-(msftx(i,j)/msfty(i,j))rdxmut(i,j)*( &
1896	+ur*(ph_old(i,k,j)-ph_old(i-1,k,j)))
1897
1898	ENDDO
1899
1900	END IF
1901
1902	END SUBROUTINE rhs_ph
1903
1904
1905	!-------------------------------------------------------------------------------
1906
1907	SUBROUTINE horizontal_pressure_gradient( ru_tend,rv_tend, &
1908	ph,alt,p,pb,al,php,cqu,cqv, &
1909	muu,muv,mu,fnm,fnp,rdnw, &
1910	cf1,cf2,cf3,rdx,rdy,msfux,msfuy,&
1911	msfvx,msfvy,msftx,msfty, &
1912	config_flags, non_hydrostatic, &
1913	top_lid, &
1914	ids, ide, jds, jde, kds, kde, &
1915	ims, ime, jms, jme, kms, kme, &
1916	its, ite, jts, jte, kts, kte )
1917
1918	IMPLICIT NONE
1919
1920	! Input data
1921
1922
1923	TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
1924
1925	LOGICAL, INTENT (IN ) :: non_hydrostatic, top_lid
1926
1927	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
1928	ims, ime, jms, jme, kms, kme, &
1929	its, ite, jts, jte, kts, kte
1930
1931	REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(IN ) :: &
1932	ph, &
1933	alt, &
1934	al, &
1935	p, &
1936	pb, &
1937	php, &
1938	cqu, &
1939	cqv
1940
1941
1942	REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(INOUT) :: &
1943	ru_tend, &
1944	rv_tend
1945
1946	REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN ) :: muu, muv, mu, &
1947	msfux, msfuy, &
1948	msfvx, msfvy, &
1949	msftx, msfty
1950
1951	REAL, DIMENSION( kms:kme ), INTENT(IN ) :: rdnw, fnm, fnp
1952
1953	REAL, INTENT(IN ) :: rdx, rdy, cf1, cf2, cf3
1954
1955	INTEGER :: i,j,k, itf, jtf, ktf, i_start, j_start
1956	REAL, DIMENSION( ims:ime, kms:kme ) :: dpn
1957	REAL :: dpx, dpy
1958
1959	LOGICAL :: specified
1960
1961	!<DESCRIPTION>
1962	!
1963	! horizontal_pressure_gradient calculates the
1964	! horizontal pressure gradient terms for the large-timestep tendency
1965	! in the horizontal momentum equations (u,v).
1966	!
1967	!</DESCRIPTION>
1968
1969	specified = .false.
1970	if(config_flags%specified .or. config_flags%nested) specified = .true.
1971
1972	! Notes on map scale factors:
1973	! Calculates the pressure gradient terms in ADT eqns 44 and 45
1974	! With upper rho -> 'mu', these are:
1975	! Eqn 30: -mu(mx/my)(1/rho)*partial dp/dx
1976	! Eqn 31: -mu(my/mx)(1/rho)*partial dp/dy
1977	!
1978	! As we are on nu, rather than height, surfaces:
1979	!
1980	! mu dp/dx = mu alpha partial dp'/dx + (nu mu partial dmubar/dx) alpha'
1981	! + mu partial dphi'/dx + (partial dphi/dx)*(partial dp'/dnu - mu')
1982	!
1983	! mu dp/dy = mu alpha partial dp'/dy + (nu mu partial dmubar/dy) alpha'
1984	! + mu partial dphi'/dy + (partial dphi/dy)*(partial dp'/dnu - mu')
1985
1986	! start with the north-south (y) pressure gradient
1987
1988	itf=MIN(ite,ide-1)
1989	jtf=jte
1990	ktf=MIN(kte,kde-1)
1991	i_start = its
1992	j_start = jts
1993	IF ( (config_flags%open_ys .or. specified .or. &
1994	config_flags%nested .or. config_flags%polar ) .and. jts == jds ) j_start = jts+1
1995	IF ( (config_flags%open_ye .or. specified .or. &
1996	config_flags%nested .or. config_flags%polar ) .and. jte == jde ) jtf = jtf-1
1997
1998	DO j = j_start, jtf
1999
2000	IF ( non_hydrostatic ) THEN
2001
2002	k=1
2003
2004	DO i = i_start, itf
2005	dpn(i,k) = .5( cf1(p(i,k ,j-1)+p(i,k ,j)) &
2006	+cf2*(p(i,k+1,j-1)+p(i,k+1,j)) &
2007	+cf3*(p(i,k+2,j-1)+p(i,k+2,j)) )
2008	dpn(i,kde) = 0.
2009	ENDDO
2010	IF (top_lid) THEN
2011	DO i = i_start, itf
2012	dpn(i,kde) = .5( cf1(p(i,kde-1,j-1)+p(i,kde-1,j)) &
2013	+cf2*(p(i,kde-2,j-1)+p(i,kde-2,j)) &
2014	+cf3*(p(i,kde-3,j-1)+p(i,kde-3,j)) )
2015	ENDDO
2016	ENDIF
2017
2018	DO k=2,ktf
2019	DO i = i_start, itf
2020	dpn(i,k) = .5( fnm(k)(p(i,k ,j-1)+p(i,k ,j)) &
2021	+fnp(k)*(p(i,k-1,j-1)+p(i,k-1,j)) )
2022	END DO
2023	END DO
2024
2025	! ADT eqn 45: -mu(my/mx)(1/rho)*partial dp/dy
2026	! [alt, al are 1/rho terms; muv, mu are NOT coupled]
2027	DO K=1,ktf
2028	DO i = i_start, itf
2029	! Here are mu dp/dy terms 1-3
2030	dpy = (msfvy(i,j)/msfvx(i,j)).5rdymuv(i,j)( &
2031	(ph (i,k+1,j)-ph (i,k+1,j-1) + ph(i,k,j)-ph(i,k,j-1)) &
2032	+(alt(i,k ,j)+alt(i,k ,j-1))*(p (i,k,j)-p (i,k,j-1)) &
2033	+(al (i,k ,j)+al (i,k ,j-1))*(pb(i,k,j)-pb(i,k,j-1)) )
2034	! Here is mu dp/dy term 4
2035	dpy = dpy + (msfvy(i,j)/msfvx(i,j))rdy(php(i,k,j)-php(i,k,j-1))* &
2036	(rdnw(k)(dpn(i,k+1)-dpn(i,k))-.5(mu(i,j-1)+mu(i,j)))
2037	rv_tend(i,k,j) = rv_tend(i,k,j)-cqv(i,k,j)*dpy
2038	END DO
2039	END DO
2040
2041	ELSE
2042
2043	! ADT eqn 45: -mu(my/mx)(1/rho)*partial dp/dy
2044	! [alt, al are 1/rho terms; muv, mu are NOT coupled]
2045	DO K=1,ktf
2046	DO i = i_start, itf
2047	! Here are mu dp/dy terms 1-3; term 4 not needed if hydrostatic
2048	dpy = (msfvy(i,j)/msfvx(i,j)).5rdymuv(i,j)( &
2049	(ph (i,k+1,j)-ph (i,k+1,j-1) + ph(i,k,j)-ph(i,k,j-1)) &
2050	+(alt(i,k ,j)+alt(i,k ,j-1))*(p (i,k,j)-p (i,k,j-1)) &
2051	+(al (i,k ,j)+al (i,k ,j-1))*(pb(i,k,j)-pb(i,k,j-1)) )
2052	rv_tend(i,k,j) = rv_tend(i,k,j)-cqv(i,k,j)*dpy
2053	END DO
2054	END DO
2055
2056	END IF
2057
2058	ENDDO
2059
2060	! now the east-west (x) pressure gradient
2061
2062	itf=ite
2063	jtf=MIN(jte,jde-1)
2064	ktf=MIN(kte,kde-1)
2065	i_start = its
2066	j_start = jts
2067	IF ( (config_flags%open_xs .or. specified .or. &
2068	config_flags%nested ) .and. its == ids ) i_start = its+1
2069	IF ( (config_flags%open_xe .or. specified .or. &
2070	config_flags%nested ) .and. ite == ide ) itf = itf-1
2071	IF ( config_flags%periodic_x ) i_start = its
2072	IF ( config_flags%periodic_x ) itf=ite
2073
2074	DO j = j_start, jtf
2075
2076	IF ( non_hydrostatic ) THEN
2077
2078	k=1
2079
2080	DO i = i_start, itf
2081	dpn(i,k) = .5( cf1(p(i-1,k ,j)+p(i,k ,j)) &
2082	+cf2*(p(i-1,k+1,j)+p(i,k+1,j)) &
2083	+cf3*(p(i-1,k+2,j)+p(i,k+2,j)) )
2084	dpn(i,kde) = 0.
2085	ENDDO
2086	IF (top_lid) THEN
2087	DO i = i_start, itf
2088	dpn(i,kde) = .5( cf1(p(i-1,kde-1,j)+p(i,kde-1,j)) &
2089	+cf2*(p(i-1,kde-2,j)+p(i,kde-2,j)) &
2090	+cf3*(p(i-1,kde-3,j)+p(i,kde-3,j)) )
2091	ENDDO
2092	ENDIF
2093
2094	DO k=2,ktf
2095	DO i = i_start, itf
2096	dpn(i,k) = .5( fnm(k)(p(i-1,k ,j)+p(i,k ,j)) &
2097	+fnp(k)*(p(i-1,k-1,j)+p(i,k-1,j)) )
2098	END DO
2099	END DO
2100
2101	! ADT eqn 44: -mu(mx/my)(1/rho)*partial dp/dx
2102	! [alt, al are 1/rho terms; muu, mu are NOT coupled]
2103	DO K=1,ktf
2104	DO i = i_start, itf
2105	! Here are mu dp/dy terms 1-3
2106	dpx = (msfux(i,j)/msfuy(i,j)).5rdxmuu(i,j)( &
2107	(ph (i,k+1,j)-ph (i-1,k+1,j) + ph(i,k,j)-ph(i-1,k,j)) &
2108	+(alt(i,k ,j)+alt(i-1,k ,j))*(p (i,k,j)-p (i-1,k,j)) &
2109	+(al (i,k ,j)+al (i-1,k ,j))*(pb(i,k,j)-pb(i-1,k,j)) )
2110	! Here is mu dp/dy term 4
2111	dpx = dpx + (msfux(i,j)/msfuy(i,j))rdx(php(i,k,j)-php(i-1,k,j))* &
2112	(rdnw(k)(dpn(i,k+1)-dpn(i,k))-.5(mu(i-1,j)+mu(i,j)))
2113	ru_tend(i,k,j) = ru_tend(i,k,j)-cqu(i,k,j)*dpx
2114	END DO
2115	END DO
2116
2117	ELSE
2118
2119	! ADT eqn 44: -mu(mx/my)(1/rho)*partial dp/dx
2120	! [alt, al are 1/rho terms; muu, mu are NOT coupled]
2121	DO K=1,ktf
2122	DO i = i_start, itf
2123	! Here are mu dp/dy terms 1-3; term 4 not needed if hydrostatic
2124	dpx = (msfux(i,j)/msfuy(i,j)).5rdxmuu(i,j)( &
2125	(ph (i,k+1,j)-ph (i-1,k+1,j) + ph(i,k,j)-ph(i-1,k,j)) &
2126	+(alt(i,k ,j)+alt(i-1,k ,j))*(p (i,k,j)-p (i-1,k,j)) &
2127	+(al (i,k ,j)+al (i-1,k ,j))*(pb(i,k,j)-pb(i-1,k,j)) )
2128	ru_tend(i,k,j) = ru_tend(i,k,j)-cqu(i,k,j)*dpx
2129	END DO
2130	END DO
2131
2132	END IF
2133
2134	ENDDO
2135
2136	END SUBROUTINE horizontal_pressure_gradient
2137
2138	!-------------------------------------------------------------------------------
2139
2140	SUBROUTINE pg_buoy_w( rw_tend, p, cqw, mu, mub, &
2141	rdnw, rdn, g, msftx, msfty, &
2142	ids, ide, jds, jde, kds, kde, &
2143	ims, ime, jms, jme, kms, kme, &
2144	its, ite, jts, jte, kts, kte )
2145
2146	IMPLICIT NONE
2147
2148	! Input data
2149
2150	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
2151	ims, ime, jms, jme, kms, kme, &
2152	its, ite, jts, jte, kts, kte
2153
2154	REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(IN ) :: p
2155	REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(INOUT) :: cqw
2156
2157
2158	REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(INOUT) :: rw_tend
2159
2160	REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN ) :: mub, mu, msftx, msfty
2161
2162	REAL, DIMENSION( kms:kme ), INTENT(IN ) :: rdnw, rdn
2163
2164	REAL, INTENT(IN ) :: g
2165
2166	INTEGER :: itf, jtf, i, j, k
2167	REAL :: cq1, cq2
2168
2169
2170	!<DESCRIPTION>
2171	!
2172	! pg_buoy_w calculates the
2173	! vertical pressure gradient and buoyancy terms for the large-timestep
2174	! tendency in the vertical momentum equation.
2175	!
2176	!</DESCRIPTION>
2177
2178	! BUOYANCY AND PRESSURE GRADIENT TERM IN W EQUATION AT TIME T
2179
2180	! Map scale factor notes
2181	! ADT eqn 46 RHS terms 6 and 7 (where 7 is "-rho g")
2182	! Dividing by my, and using mu and nu (see Klemp et al. eqns 32, 40)
2183	! term 6: +(g/my) partial dp'/dnu
2184	! term 7: -(g/my) mu'
2185	!
2186	! For moisture-free atmosphere, cq1=1, cq2=0
2187	! => (1./msft(i,j)) * g * [rdn(k)*{p(i,k,j)-p(i,k-1,j)}-mu(i,j)]
2188
2189	itf=MIN(ite,ide-1)
2190	jtf=MIN(jte,jde-1)
2191
2192	DO j = jts,jtf
2193
2194	k=kde
2195	DO i=its,itf
2196	cq1 = 1./(1.+cqw(i,k-1,j))
2197	cq2 = cqw(i,k-1,j)*cq1
2198	rw_tend(i,k,j) = rw_tend(i,k,j)+(1./msfty(i,j))g( &
2199	cq12.rdnw(k-1)*( -p(i,k-1,j)) &
2200	-mu(i,j)-cq2*mub(i,j) )
2201	END DO
2202
2203	DO k = 2, kde-1
2204	DO i = its,itf
2205	cq1 = 1./(1.+cqw(i,k,j))
2206	cq2 = cqw(i,k,j)*cq1
2207	cqw(i,k,j) = cq1
2208	rw_tend(i,k,j) = rw_tend(i,k,j)+(1./msfty(i,j))g( &
2209	cq1rdn(k)(p(i,k,j)-p(i,k-1,j)) &
2210	-mu(i,j)-cq2*mub(i,j) )
2211	END DO
2212	ENDDO
2213
2214
2215	ENDDO
2216
2217	END SUBROUTINE pg_buoy_w
2218
2219	!-------------------------------------------------------------------------------
2220
2221	SUBROUTINE w_damp( rw_tend, max_vert_cfl,max_horiz_cfl, &
2222	u, v, ww, w, mut, rdnw, &
2223	rdx, rdy, msfux, msfuy, &
2224	msfvx, msfvy, dt, &
2225	config_flags, &
2226	ids, ide, jds, jde, kds, kde, &
2227	ims, ime, jms, jme, kms, kme, &
2228	its, ite, jts, jte, kts, kte )
2229
2230	USE module_llxy
2231	IMPLICIT NONE
2232
2233	! Input data
2234
2235	TYPE(grid_config_rec_type ) , INTENT(IN ) :: config_flags
2236
2237	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
2238	ims, ime, jms, jme, kms, kme, &
2239	its, ite, jts, jte, kts, kte
2240
2241	REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(IN ) :: u, v, ww, w
2242
2243	REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(INOUT) :: rw_tend
2244
2245	REAL, INTENT(OUT) :: max_vert_cfl
2246	REAL, INTENT(OUT) :: max_horiz_cfl
2247	REAL :: horiz_cfl
2248
2249	REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN ) :: mut
2250
2251	REAL, DIMENSION( kms:kme ), INTENT(IN ) :: rdnw
2252
2253	REAL, INTENT(IN) :: dt
2254	REAL, INTENT(IN) :: rdx, rdy
2255	REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: msfux, msfuy
2256	REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: msfvx, msfvy
2257
2258	REAL :: vert_cfl, cf_n, cf_d, maxdub, maxdeta
2259
2260	INTEGER :: itf, jtf, i, j, k, maxi, maxj, maxk
2261	INTEGER :: some
2262	CHARACTER*512 :: temp
2263
2264	CHARACTER (LEN=256) :: time_str
2265	CHARACTER (LEN=256) :: grid_str
2266
2267	integer :: total
2268	REAL :: msfuxt , msfxffl
2269
2270	!<DESCRIPTION>
2271	!
2272	! w_damp computes a damping term for the vertical velocity when the
2273	! vertical Courant number is too large. This was found to be preferable to
2274	! decreasing the timestep or increasing the diffusion in real-data applications
2275	! that produced potentially-unstable large vertical velocities because of
2276	! unphysically large heating rates coming from the cumulus parameterization
2277	! schemes run at moderately high resolutions (dx ~ O(10) km).
2278	!
2279	! Additionally, w_damp returns the maximum cfl values due to vertical motion and
2280	! horizontal motion. These values are returned via the max_vert_cfl and
2281	! max_horiz_cfl variables. (Added by T. Hutchinson, WSI, 3/5/2007)
2282	!
2283	!</DESCRIPTION>
2284
2285	itf=MIN(ite,ide-1)
2286	jtf=MIN(jte,jde-1)
2287
2288	some = 0
2289	max_vert_cfl = 0.
2290	max_horiz_cfl = 0.
2291	total = 0
2292
2293	IF(config_flags%map_proj == PROJ_CASSINI ) then
2294	msfxffl = 1.0/COS(config_flags%fft_filter_lat*degrad)
2295	END IF
2296
2297	IF ( config_flags%w_damping == 1 ) THEN
2298	DO j = jts,jtf
2299
2300	DO k = 2, kde-1
2301	DO i = its,itf
2302	#if 1
2303	IF(config_flags%map_proj == PROJ_CASSINI ) then
2304	msfuxt = MIN(msfux(i,j), msfxffl)
2305	ELSE
2306	msfuxt = msfux(i,j)
2307	END IF
2308	vert_cfl = abs(ww(i,k,j)/mut(i,j)rdnw(k)dt)
2309
2310	IF ( vert_cfl > max_vert_cfl ) THEN
2311	max_vert_cfl = vert_cfl ; maxi = i ; maxj = j ; maxk = k
2312	maxdub = w(i,k,j) ; maxdeta = -1./rdnw(k)
2313	ENDIF
2314
2315	horiz_cfl = max( abs(u(i,k,j) * rdx * msfuxt * dt), &
2316	abs(v(i,k,j) * rdy * msfvy(i,j) * dt) )
2317	if (horiz_cfl > max_horiz_cfl) then
2318	max_horiz_cfl = horiz_cfl
2319	endif
2320
2321	if(vert_cfl .gt. w_beta)then
2322	#else
2323	! restructure to get rid of divide
2324	!
2325	! This had been used for efficiency, but with the addition of returning the cfl values,
2326	! the old version (above) was reinstated. (T. Hutchinson, 3/5/2007)
2327	!
2328	cf_n = abs(ww(i,k,j)rdnw(k)dt)
2329	cf_d = abs(mut(i,j))
2330	if(cf_n .gt. cf_d*w_beta )then
2331	#endif
2332
2333	WRITE(temp,*)i,j,k,' vert_cfl,w,d(eta)=',vert_cfl,w(i,k,j),-1./rdnw(k)
2334	CALL wrf_debug ( 100 , TRIM(temp) )
2335	if ( vert_cfl > 2. ) some = some + 1
2336	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)
2337	endif
2338	END DO
2339	ENDDO
2340	ENDDO
2341	ELSE
2342	! just print
2343	DO j = jts,jtf
2344
2345	DO k = 2, kde-1
2346	DO i = its,itf
2347
2348	#if 1
2349	IF(config_flags%map_proj == PROJ_CASSINI ) then
2350	msfuxt = MIN(msfux(i,j), msfxffl)
2351	ELSE
2352	msfuxt = msfux(i,j)
2353	END IF
2354	vert_cfl = abs(ww(i,k,j)/mut(i,j)rdnw(k)dt)
2355
2356	IF ( vert_cfl > max_vert_cfl ) THEN
2357	max_vert_cfl = vert_cfl ; maxi = i ; maxj = j ; maxk = k
2358	maxdub = w(i,k,j) ; maxdeta = -1./rdnw(k)
2359	ENDIF
2360
2361	horiz_cfl = max( abs(u(i,k,j) * rdx * msfuxt * dt), &
2362	abs(v(i,k,j) * rdy * msfvy(i,j) * dt) )
2363
2364	if (horiz_cfl > max_horiz_cfl) then
2365	max_horiz_cfl = horiz_cfl
2366	endif
2367
2368	if(vert_cfl .gt. w_beta)then
2369	#else
2370	! restructure to get rid of divide
2371	!
2372	! This had been used for efficiency, but with the addition of returning the cfl values,
2373	! the old version (above) was reinstated. (T. Hutchinson, 3/5/2007)
2374	!
2375	cf_n = abs(ww(i,k,j)rdnw(k)dt)
2376	cf_d = abs(mut(i,j))
2377	if(cf_n .gt. cf_d*w_beta )then
2378	#endif
2379	WRITE(temp,*)i,j,k,' vert_cfl,w,d(eta)=',vert_cfl,w(i,k,j),-1./rdnw(k)
2380	CALL wrf_debug ( 100 , TRIM(temp) )
2381	if ( vert_cfl > 2. ) some = some + 1
2382	endif
2383	END DO
2384	ENDDO
2385	ENDDO
2386	ENDIF
2387	IF ( some .GT. 0 ) THEN
2388	CALL get_current_time_string( time_str )
2389	CALL get_current_grid_name( grid_str )
2390	WRITE(wrf_err_message,*)some, &
2391	' points exceeded cfl=2 in domain '//TRIM(grid_str)//' at time '//TRIM(time_str)//' hours'
2392	CALL wrf_debug ( 0 , TRIM(wrf_err_message) )
2393	WRITE(wrf_err_message,*)'MAX AT i,j,k: ',maxi,maxj,maxk,' vert_cfl,w,d(eta)=',max_vert_cfl, &
2394	maxdub,maxdeta
2395	CALL wrf_debug ( 0 , TRIM(wrf_err_message) )
2396	ENDIF
2397
2398	END SUBROUTINE w_damp
2399
2400	!-------------------------------------------------------------------------------
2401
2402	SUBROUTINE horizontal_diffusion ( name, field, tendency, mu, &
2403	config_flags, &
2404	msfux, msfuy, msfvx, msfvx_inv, &
2405	msfvy, msftx, msfty, &
2406	khdif, xkmhd, rdx, rdy, &
2407	ids, ide, jds, jde, kds, kde, &
2408	ims, ime, jms, jme, kms, kme, &
2409	its, ite, jts, jte, kts, kte )
2410
2411	IMPLICIT NONE
2412
2413	! Input data
2414
2415	TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
2416
2417	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
2418	ims, ime, jms, jme, kms, kme, &
2419	its, ite, jts, jte, kts, kte
2420
2421	CHARACTER(LEN=1) , INTENT(IN ) :: name
2422
2423	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: field, xkmhd
2424
2425	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency
2426
2427	REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mu
2428
2429	REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: msfux, &
2430	msfuy, &
2431	msfvx, &
2432	msfvx_inv, &
2433	msfvy, &
2434	msftx, &
2435	msfty
2436
2437	REAL , INTENT(IN ) :: rdx, &
2438	rdy, &
2439	khdif
2440
2441	! Local data
2442
2443	INTEGER :: i, j, k, itf, jtf, ktf
2444
2445	INTEGER :: i_start, i_end, j_start, j_end
2446
2447	REAL :: mrdx, mkrdxm, mkrdxp, &
2448	mrdy, mkrdym, mkrdyp
2449
2450	LOGICAL :: specified
2451
2452	!<DESCRIPTION>
2453	!
2454	! horizontal_diffusion computes the horizontal diffusion tendency
2455	! on model horizontal coordinate surfaces.
2456	!
2457	!</DESCRIPTION>
2458
2459	specified = .false.
2460	if(config_flags%specified .or. config_flags%nested) specified = .true.
2461
2462	ktf=MIN(kte,kde-1)
2463
2464	IF (name .EQ. 'u') THEN
2465
2466	i_start = its
2467	i_end = ite
2468	j_start = jts
2469	j_end = MIN(jte,jde-1)
2470
2471	IF ( config_flags%open_xs .or. specified ) i_start = MAX(ids+1,its)
2472	IF ( config_flags%open_xe .or. specified ) i_end = MIN(ide-1,ite)
2473	IF ( config_flags%open_ys .or. specified ) j_start = MAX(jds+1,jts)
2474	IF ( config_flags%open_ye .or. specified ) j_end = MIN(jde-2,jte)
2475	IF ( config_flags%periodic_x ) i_start = its
2476	IF ( config_flags%periodic_x ) i_end = ite
2477
2478
2479	DO j = j_start, j_end
2480	DO k=kts,ktf
2481	DO i = i_start, i_end
2482
2483	! The interior is grad: (m_xd/dx), the exterior is div: (m_xm_y*d/dx(/m_y))
2484	! setting up different averagings of m^2 partial d/dX and m^2 partial d/dY
2485
2486	mkrdxm=(msftx(i-1,j)/msfty(i-1,j))mu(i-1,j)xkmhd(i-1,k,j)*rdx
2487	mkrdxp=(msftx(i,j)/msfty(i,j))mu(i,j)xkmhd(i,k,j)*rdx
2488	mrdx=msfux(i,j)msfuy(i,j)rdx
2489	mkrdym=( (msfuy(i,j)+msfuy(i,j-1))/(msfux(i,j)+msfux(i,j-1)) )* &
2490	0.25(mu(i,j)+mu(i,j-1)+mu(i-1,j-1)+mu(i-1,j)) &
2491	0.25(xkmhd(i,k,j)+xkmhd(i,k,j-1)+xkmhd(i-1,k,j-1)+xkmhd(i-1,k,j))rdy
2492	mkrdyp=( (msfuy(i,j)+msfuy(i,j+1))/(msfux(i,j)+msfux(i,j+1)) )* &
2493	0.25(mu(i,j)+mu(i,j+1)+mu(i-1,j+1)+mu(i-1,j)) &
2494	0.25(xkmhd(i,k,j)+xkmhd(i,k,j+1)+xkmhd(i-1,k,j+1)+xkmhd(i-1,k,j))rdy
2495	! need to do four-corners (t) for diffusion coefficient as there are
2496	! no values at u,v points
2497	! msfuy - has to be y as part of d/dY
2498	! has to be u as we're at a u point
2499	mrdy=msfux(i,j)msfuy(i,j)rdy
2500
2501	! correctly averaged version of rho~ * m^2 *
2502	! [partial d/dX(partial du^/dX) + partial d/dY(partial du^/dY)]
2503	tendency(i,k,j)=tendency(i,k,j)+( &
2504	mrdx(mkrdxp(field(i+1,k,j)-field(i ,k,j)) &
2505	-mkrdxm*(field(i ,k,j)-field(i-1,k,j))) &
2506	+mrdy(mkrdyp(field(i,k,j+1)-field(i,k,j )) &
2507	-mkrdym*(field(i,k,j )-field(i,k,j-1))))
2508	ENDDO
2509	ENDDO
2510	ENDDO
2511
2512	ELSE IF (name .EQ. 'v')THEN
2513
2514	i_start = its
2515	i_end = MIN(ite,ide-1)
2516	j_start = jts
2517	j_end = jte
2518
2519	IF ( config_flags%open_xs .or. specified ) i_start = MAX(ids+1,its)
2520	IF ( config_flags%open_xe .or. specified ) i_end = MIN(ide-2,ite)
2521	IF ( config_flags%open_ys .or. specified ) j_start = MAX(jds+1,jts)
2522	IF ( config_flags%open_ye .or. specified ) j_end = MIN(jde-1,jte)
2523	IF ( config_flags%periodic_x ) i_start = its
2524	IF ( config_flags%periodic_x ) i_end = MIN(ite,ide-1)
2525	IF ( config_flags%polar ) j_start = MAX(jds+1,jts)
2526	IF ( config_flags%polar ) j_end = MIN(jde-1,jte)
2527
2528	DO j = j_start, j_end
2529	DO k=kts,ktf
2530	DO i = i_start, i_end
2531
2532	mkrdxm=( (msfvx(i,j)+msfvx(i-1,j))/(msfvy(i,j)+msfvy(i-1,j)) )* &
2533	0.25(mu(i,j)+mu(i,j-1)+mu(i-1,j-1)+mu(i-1,j)) &
2534	0.25(xkmhd(i,k,j)+xkmhd(i,k,j-1)+xkmhd(i-1,k,j-1)+xkmhd(i-1,k,j))rdx
2535	mkrdxp=( (msfvx(i,j)+msfvx(i+1,j))/(msfvy(i,j)+msfvy(i+1,j)) )* &
2536	0.25(mu(i,j)+mu(i,j-1)+mu(i+1,j-1)+mu(i+1,j)) &
2537	0.25(xkmhd(i,k,j)+xkmhd(i,k,j-1)+xkmhd(i+1,k,j-1)+xkmhd(i+1,k,j))rdx
2538	mrdx=msfvx(i,j)msfvy(i,j)rdx
2539	mkrdym=(msfty(i,j-1)/msftx(i,j-1))xkmhd(i,k,j-1)rdy
2540	mkrdyp=(msfty(i,j)/msftx(i,j))xkmhd(i,k,j)rdy
2541	mrdy=msfvx(i,j)msfvy(i,j)rdy
2542
2543	tendency(i,k,j)=tendency(i,k,j)+( &
2544	mrdx(mkrdxp(field(i+1,k,j)-field(i ,k,j)) &
2545	-mkrdxm*(field(i ,k,j)-field(i-1,k,j))) &
2546	+mrdy(mkrdyp(field(i,k,j+1)-field(i,k,j )) &
2547	-mkrdym*(field(i,k,j )-field(i,k,j-1))))
2548	ENDDO
2549	ENDDO
2550	ENDDO
2551
2552	ELSE IF (name .EQ. 'w')THEN
2553
2554	i_start = its
2555	i_end = MIN(ite,ide-1)
2556	j_start = jts
2557	j_end = MIN(jte,jde-1)
2558
2559	IF ( config_flags%open_xs .or. specified ) i_start = MAX(ids+1,its)
2560	IF ( config_flags%open_xe .or. specified ) i_end = MIN(ide-2,ite)
2561	IF ( config_flags%open_ys .or. specified ) j_start = MAX(jds+1,jts)
2562	IF ( config_flags%open_ye .or. specified ) j_end = MIN(jde-2,jte)
2563	IF ( config_flags%periodic_x ) i_start = its
2564	IF ( config_flags%periodic_x ) i_end = MIN(ite,ide-1)
2565
2566	DO j = j_start, j_end
2567	DO k=kts+1,ktf
2568	DO i = i_start, i_end
2569
2570	mkrdxm=(msfux(i,j)/msfuy(i,j))* &
2571	0.25(mu(i,j)+mu(i-1,j)+mu(i,j)+mu(i-1,j)) &
2572	0.25(xkmhd(i,k,j)+xkmhd(i-1,k,j)+xkmhd(i,k-1,j)+xkmhd(i-1,k-1,j))rdx
2573	mkrdxp=(msfux(i+1,j)/msfuy(i+1,j))* &
2574	0.25(mu(i+1,j)+mu(i,j)+mu(i+1,j)+mu(i,j)) &
2575	0.25(xkmhd(i+1,k,j)+xkmhd(i,k,j)+xkmhd(i+1,k-1,j)+xkmhd(i,k-1,j))rdx
2576	mrdx=msftx(i,j)msfty(i,j)rdx
2577	! mkrdym=(msfvy(i,j)/msfvx(i,j))* &
2578	mkrdym=(msfvy(i,j)msfvx_inv(i,j)) &
2579	0.25(mu(i,j)+mu(i,j-1)+mu(i,j)+mu(i,j-1)) &
2580	0.25(xkmhd(i,k,j)+xkmhd(i,k,j-1)+xkmhd(i,k-1,j)+xkmhd(i,k-1,j-1))rdy
2581	! mkrdyp=(msfvy(i,j+1)/msfvx(i,j+1))* &
2582	mkrdyp=(msfvy(i,j+1)msfvx_inv(i,j+1)) &
2583	0.25(mu(i,j+1)+mu(i,j)+mu(i,j+1)+mu(i,j)) &
2584	0.25(xkmhd(i,k,j+1)+xkmhd(i,k,j)+xkmhd(i,k-1,j+1)+xkmhd(i,k-1,j))rdy
2585	mrdy=msftx(i,j)msfty(i,j)rdy
2586
2587	tendency(i,k,j)=tendency(i,k,j)+( &
2588	mrdx(mkrdxp(field(i+1,k,j)-field(i ,k,j)) &
2589	-mkrdxm*(field(i ,k,j)-field(i-1,k,j))) &
2590	+mrdy(mkrdyp(field(i,k,j+1)-field(i,k,j )) &
2591	-mkrdym*(field(i,k,j )-field(i,k,j-1))))
2592	ENDDO
2593	ENDDO
2594	ENDDO
2595
2596	ELSE
2597
2598
2599	i_start = its
2600	i_end = MIN(ite,ide-1)
2601	j_start = jts
2602	j_end = MIN(jte,jde-1)
2603
2604	IF ( config_flags%open_xs .or. specified ) i_start = MAX(ids+1,its)
2605	IF ( config_flags%open_xe .or. specified ) i_end = MIN(ide-2,ite)
2606	IF ( config_flags%open_ys .or. specified ) j_start = MAX(jds+1,jts)
2607	IF ( config_flags%open_ye .or. specified ) j_end = MIN(jde-2,jte)
2608	IF ( config_flags%periodic_x ) i_start = its
2609	IF ( config_flags%periodic_x ) i_end = MIN(ite,ide-1)
2610
2611	DO j = j_start, j_end
2612	DO k=kts,ktf
2613	DO i = i_start, i_end
2614
2615	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
2616	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
2617	mrdx=msftx(i,j)msfty(i,j)rdx
2618	! 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
2619	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
2620	! 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
2621	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
2622	mrdy=msftx(i,j)msfty(i,j)rdy
2623
2624	tendency(i,k,j)=tendency(i,k,j)+( &
2625	mrdx(mkrdxp(field(i+1,k,j)-field(i ,k,j)) &
2626	-mkrdxm*(field(i ,k,j)-field(i-1,k,j))) &
2627	+mrdy(mkrdyp(field(i,k,j+1)-field(i,k,j )) &
2628	-mkrdym*(field(i,k,j )-field(i,k,j-1))))
2629	ENDDO
2630	ENDDO
2631	ENDDO
2632
2633	ENDIF
2634
2635	END SUBROUTINE horizontal_diffusion
2636
2637	!-----------------------------------------------------------------------------------------
2638
2639	SUBROUTINE horizontal_diffusion_3dmp ( name, field, tendency, mu, &
2640	config_flags, base_3d, &
2641	msfux, msfuy, msfvx, msfvx_inv, &
2642	msfvy, msftx, msfty, &
2643	khdif, xkmhd, rdx, rdy, &
2644	ids, ide, jds, jde, kds, kde, &
2645	ims, ime, jms, jme, kms, kme, &
2646	its, ite, jts, jte, kts, kte )
2647
2648	IMPLICIT NONE
2649
2650	! Input data
2651
2652	TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
2653
2654	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
2655	ims, ime, jms, jme, kms, kme, &
2656	its, ite, jts, jte, kts, kte
2657
2658	CHARACTER(LEN=1) , INTENT(IN ) :: name
2659
2660	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: field, &
2661	xkmhd, &
2662	base_3d
2663
2664	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency
2665
2666	REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mu
2667
2668	REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: msfux, &
2669	msfuy, &
2670	msfvx, &
2671	msfvx_inv, &
2672	msfvy, &
2673	msftx, &
2674	msfty
2675
2676	REAL , INTENT(IN ) :: rdx, &
2677	rdy, &
2678	khdif
2679
2680	! Local data
2681
2682	INTEGER :: i, j, k, itf, jtf, ktf
2683
2684	INTEGER :: i_start, i_end, j_start, j_end
2685
2686	REAL :: mrdx, mkrdxm, mkrdxp, &
2687	mrdy, mkrdym, mkrdyp
2688
2689	LOGICAL :: specified
2690
2691	!<DESCRIPTION>
2692	!
2693	! horizontal_diffusion_3dmp computes the horizontal diffusion tendency
2694	! on model horizontal coordinate surfaces. This routine computes diffusion
2695	! a perturbation scalar (field-base_3d).
2696	!
2697	!</DESCRIPTION>
2698
2699	specified = .false.
2700	if(config_flags%specified .or. config_flags%nested) specified = .true.
2701
2702	ktf=MIN(kte,kde-1)
2703
2704	i_start = its
2705	i_end = MIN(ite,ide-1)
2706	j_start = jts
2707	j_end = MIN(jte,jde-1)
2708
2709	IF ( config_flags%open_xs .or. specified ) i_start = MAX(ids+1,its)
2710	IF ( config_flags%open_xe .or. specified ) i_end = MIN(ide-2,ite)
2711	IF ( config_flags%open_ys .or. specified ) j_start = MAX(jds+1,jts)
2712	IF ( config_flags%open_ye .or. specified ) j_end = MIN(jde-2,jte)
2713	IF ( config_flags%periodic_x ) i_start = its
2714	IF ( config_flags%periodic_x ) i_end = MIN(ite,ide-1)
2715
2716	DO j = j_start, j_end
2717	DO k=kts,ktf
2718	DO i = i_start, i_end
2719
2720	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
2721	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
2722	mrdx=msftx(i,j)msfty(i,j)rdx
2723	! 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
2724	! 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
2725	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
2726	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
2727	mrdy=msftx(i,j)msfty(i,j)rdy
2728
2729	tendency(i,k,j)=tendency(i,k,j)+( &
2730	mrdx( mkrdxp( field(i+1,k,j) -field(i ,k,j) &
2731	-base_3d(i+1,k,j)+base_3d(i ,k,j) ) &
2732	-mkrdxm*( field(i ,k,j) -field(i-1,k,j) &
2733	-base_3d(i ,k,j)+base_3d(i-1,k,j) ) ) &
2734	+mrdy( mkrdyp( field(i,k,j+1) -field(i,k,j ) &
2735	-base_3d(i,k,j+1)+base_3d(i,k,j ) ) &
2736	-mkrdym*( field(i,k,j ) -field(i,k,j-1) &
2737	-base_3d(i,k,j )+base_3d(i,k,j-1) ) ) &
2738	)
2739	ENDDO
2740	ENDDO
2741	ENDDO
2742
2743	END SUBROUTINE horizontal_diffusion_3dmp
2744
2745	!-----------------------------------------------------------------------------------------
2746
2747	SUBROUTINE vertical_diffusion ( name, field, tendency, &
2748	config_flags, &
2749	alt, mut, rdn, rdnw, kvdif, &
2750	ids, ide, jds, jde, kds, kde, &
2751	ims, ime, jms, jme, kms, kme, &
2752	its, ite, jts, jte, kts, kte )
2753
2754
2755	IMPLICIT NONE
2756
2757	! Input data
2758
2759	TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
2760
2761	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
2762	ims, ime, jms, jme, kms, kme, &
2763	its, ite, jts, jte, kts, kte
2764
2765	CHARACTER(LEN=1) , INTENT(IN ) :: name
2766
2767	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , &
2768	INTENT(IN ) :: field, &
2769	alt
2770
2771	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency
2772
2773	REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mut
2774
2775	REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: rdn, rdnw
2776
2777	REAL , INTENT(IN ) :: kvdif
2778
2779	! Local data
2780
2781	INTEGER :: i, j, k, itf, jtf, ktf
2782	INTEGER :: i_start, i_end, j_start, j_end
2783
2784	REAL , DIMENSION(its:ite, jts:jte) :: vfluxm, vfluxp, zz
2785	REAL , DIMENSION(its:ite, 0:kte+1) :: vflux
2786
2787	REAL :: rdz
2788
2789	LOGICAL :: specified
2790
2791	!<DESCRIPTION>
2792	!
2793	! vertical_diffusion
2794	! computes vertical diffusion tendency.
2795	!
2796	!</DESCRIPTION>
2797
2798	specified = .false.
2799	if(config_flags%specified .or. config_flags%nested) specified = .true.
2800
2801	ktf=MIN(kte,kde-1)
2802
2803	IF (name .EQ. 'w')THEN
2804
2805
2806	i_start = its
2807	i_end = MIN(ite,ide-1)
2808	j_start = jts
2809	j_end = MIN(jte,jde-1)
2810
2811	j_loop_w : DO j = j_start, j_end
2812
2813	DO k=kts,ktf-1
2814	DO i = i_start, i_end
2815	vflux(i,k)= (kvdif/alt(i,k,j))rdnw(k)(field(i,k+1,j)-field(i,k,j))
2816	ENDDO
2817	ENDDO
2818
2819	DO i = i_start, i_end
2820	vflux(i,ktf)=0.
2821	ENDDO
2822
2823	DO k=kts+1,ktf
2824	DO i = i_start, i_end
2825	tendency(i,k,j)=tendency(i,k,j) &
2826	+rdn(k)gg/mut(i,j)/(0.5*(alt(i,k,j)+alt(i,k-1,j))) &
2827	*(vflux(i,k)-vflux(i,k-1))
2828	ENDDO
2829	ENDDO
2830
2831	ENDDO j_loop_w
2832
2833	ELSE IF(name .EQ. 'm')THEN
2834
2835	i_start = its
2836	i_end = MIN(ite,ide-1)
2837	j_start = jts
2838	j_end = MIN(jte,jde-1)
2839
2840	j_loop_s : DO j = j_start, j_end
2841
2842	DO k=kts,ktf-1
2843	DO i = i_start, i_end
2844	vflux(i,k)=kvdifrdn(k+1)/(0.5(alt(i,k,j)+alt(i,k+1,j))) &
2845	*(field(i,k+1,j)-field(i,k,j))
2846	ENDDO
2847	ENDDO
2848
2849	DO i = i_start, i_end
2850	vflux(i,0)=vflux(i,1)
2851	ENDDO
2852
2853	DO i = i_start, i_end
2854	vflux(i,ktf)=0.
2855	ENDDO
2856
2857	DO k=kts,ktf
2858	DO i = i_start, i_end
2859	tendency(i,k,j)=tendency(i,k,j)+g*g/mut(i,j)/alt(i,k,j) &
2860	rdnw(k)(vflux(i,k)-vflux(i,k-1))
2861	ENDDO
2862	ENDDO
2863
2864	ENDDO j_loop_s
2865
2866	ENDIF
2867
2868	END SUBROUTINE vertical_diffusion
2869
2870
2871	!-------------------------------------------------------------------------------
2872
2873	SUBROUTINE vertical_diffusion_mp ( field, tendency, config_flags, &
2874	base, &
2875	alt, mut, rdn, rdnw, kvdif, &
2876	ids, ide, jds, jde, kds, kde, &
2877	ims, ime, jms, jme, kms, kme, &
2878	its, ite, jts, jte, kts, kte )
2879
2880
2881	IMPLICIT NONE
2882
2883	! Input data
2884
2885	TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
2886
2887	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
2888	ims, ime, jms, jme, kms, kme, &
2889	its, ite, jts, jte, kts, kte
2890
2891	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , &
2892	INTENT(IN ) :: field, &
2893	alt
2894
2895	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency
2896
2897	REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mut
2898
2899	REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: rdn, &
2900	rdnw, &
2901	base
2902
2903	REAL , INTENT(IN ) :: kvdif
2904
2905	! Local data
2906
2907	INTEGER :: i, j, k, itf, jtf, ktf
2908	INTEGER :: i_start, i_end, j_start, j_end
2909
2910	REAL , DIMENSION(its:ite, 0:kte+1) :: vflux
2911
2912	REAL :: rdz
2913
2914	LOGICAL :: specified
2915
2916	!<DESCRIPTION>
2917	!
2918	! vertical_diffusion_mp
2919	! computes vertical diffusion tendency of a perturbation variable
2920	! (field-base). Note that base as a 1D (k) field.
2921	!
2922	!</DESCRIPTION>
2923
2924	specified = .false.
2925	if(config_flags%specified .or. config_flags%nested) specified = .true.
2926
2927	ktf=MIN(kte,kde-1)
2928
2929	i_start = its
2930	i_end = MIN(ite,ide-1)
2931	j_start = jts
2932	j_end = MIN(jte,jde-1)
2933
2934	j_loop_s : DO j = j_start, j_end
2935
2936	DO k=kts,ktf-1
2937	DO i = i_start, i_end
2938	vflux(i,k)=kvdifrdn(k+1)/(0.5(alt(i,k,j)+alt(i,k+1,j))) &
2939	*(field(i,k+1,j)-field(i,k,j)-base(k+1)+base(k))
2940	ENDDO
2941	ENDDO
2942
2943	DO i = i_start, i_end
2944	vflux(i,0)=vflux(i,1)
2945	ENDDO
2946
2947	DO i = i_start, i_end
2948	vflux(i,ktf)=0.
2949	ENDDO
2950
2951	DO k=kts,ktf
2952	DO i = i_start, i_end
2953	tendency(i,k,j)=tendency(i,k,j)+g*g/mut(i,j)/alt(i,k,j) &
2954	rdnw(k)(vflux(i,k)-vflux(i,k-1))
2955	ENDDO
2956	ENDDO
2957
2958	ENDDO j_loop_s
2959
2960	END SUBROUTINE vertical_diffusion_mp
2961
2962
2963	!-------------------------------------------------------------------------------
2964
2965	SUBROUTINE vertical_diffusion_3dmp ( field, tendency, config_flags, &
2966	base_3d, &
2967	alt, mut, rdn, rdnw, kvdif, &
2968	ids, ide, jds, jde, kds, kde, &
2969	ims, ime, jms, jme, kms, kme, &
2970	its, ite, jts, jte, kts, kte )
2971
2972
2973	IMPLICIT NONE
2974
2975	! Input data
2976
2977	TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
2978
2979	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
2980	ims, ime, jms, jme, kms, kme, &
2981	its, ite, jts, jte, kts, kte
2982
2983	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , &
2984	INTENT(IN ) :: field, &
2985	alt, &
2986	base_3d
2987
2988	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency
2989
2990	REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mut
2991
2992	REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: rdn, &
2993	rdnw
2994
2995	REAL , INTENT(IN ) :: kvdif
2996
2997	! Local data
2998
2999	INTEGER :: i, j, k, itf, jtf, ktf
3000	INTEGER :: i_start, i_end, j_start, j_end
3001
3002	REAL , DIMENSION(its:ite, 0:kte+1) :: vflux
3003
3004	REAL :: rdz
3005
3006	LOGICAL :: specified
3007
3008	!<DESCRIPTION>
3009	!
3010	! vertical_diffusion_3dmp
3011	! computes vertical diffusion tendency of a perturbation variable
3012	! (field-base_3d).
3013	!
3014	!</DESCRIPTION>
3015
3016	specified = .false.
3017	if(config_flags%specified .or. config_flags%nested) specified = .true.
3018
3019	ktf=MIN(kte,kde-1)
3020
3021	i_start = its
3022	i_end = MIN(ite,ide-1)
3023	j_start = jts
3024	j_end = MIN(jte,jde-1)
3025
3026	j_loop_s : DO j = j_start, j_end
3027
3028	DO k=kts,ktf-1
3029	DO i = i_start, i_end
3030	vflux(i,k)=kvdifrdn(k+1)/(0.5(alt(i,k,j)+alt(i,k+1,j))) &
3031	*( field(i,k+1,j) -field(i,k,j) &
3032	-base_3d(i,k+1,j)+base_3d(i,k,j) )
3033	ENDDO
3034	ENDDO
3035
3036	DO i = i_start, i_end
3037	vflux(i,0)=vflux(i,1)
3038	ENDDO
3039
3040	DO i = i_start, i_end
3041	vflux(i,ktf)=0.
3042	ENDDO
3043
3044	DO k=kts,ktf
3045	DO i = i_start, i_end
3046	tendency(i,k,j)=tendency(i,k,j)+g*g/mut(i,j)/alt(i,k,j) &
3047	rdnw(k)(vflux(i,k)-vflux(i,k-1))
3048	ENDDO
3049	ENDDO
3050
3051	ENDDO j_loop_s
3052
3053	END SUBROUTINE vertical_diffusion_3dmp
3054
3055
3056	!-------------------------------------------------------------------------------
3057
3058
3059	SUBROUTINE vertical_diffusion_u ( field, tendency, &
3060	config_flags, u_base, &
3061	alt, muu, rdn, rdnw, kvdif, &
3062	ids, ide, jds, jde, kds, kde, &
3063	ims, ime, jms, jme, kms, kme, &
3064	its, ite, jts, jte, kts, kte )
3065
3066
3067	IMPLICIT NONE
3068
3069	! Input data
3070
3071	TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
3072
3073	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
3074	ims, ime, jms, jme, kms, kme, &
3075	its, ite, jts, jte, kts, kte
3076
3077	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , &
3078	INTENT(IN ) :: field, &
3079	alt
3080
3081	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency
3082
3083	REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: muu
3084
3085	REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: rdn, rdnw, u_base
3086
3087	REAL , INTENT(IN ) :: kvdif
3088
3089	! Local data
3090
3091	INTEGER :: i, j, k, itf, jtf, ktf
3092	INTEGER :: i_start, i_end, j_start, j_end
3093
3094	REAL , DIMENSION(its:ite, 0:kte+1) :: vflux
3095
3096	REAL :: rdz, zz
3097
3098	LOGICAL :: specified
3099
3100	!<DESCRIPTION>
3101	!
3102	! vertical_diffusion_u computes vertical diffusion tendency for
3103	! the u momentum equation. This routine assumes a constant eddy
3104	! viscosity kvdif.
3105	!
3106	!</DESCRIPTION>
3107
3108	specified = .false.
3109	if(config_flags%specified .or. config_flags%nested) specified = .true.
3110
3111	ktf=MIN(kte,kde-1)
3112
3113	i_start = its
3114	i_end = ite
3115	j_start = jts
3116	j_end = MIN(jte,jde-1)
3117
3118	IF ( config_flags%open_xs .or. specified ) i_start = MAX(ids+1,its)
3119	IF ( config_flags%open_xe .or. specified ) i_end = MIN(ide-1,ite)
3120	IF ( config_flags%periodic_x ) i_start = its
3121	IF ( config_flags%periodic_x ) i_end = ite
3122
3123
3124	j_loop_u : DO j = j_start, j_end
3125
3126	DO k=kts,ktf-1
3127	DO i = i_start, i_end
3128	vflux(i,k)=kvdifrdn(k+1)/(0.25( alt(i ,k ,j) &
3129	+alt(i-1,k ,j) &
3130	+alt(i ,k+1,j) &
3131	+alt(i-1,k+1,j) ) ) &
3132	*(field(i,k+1,j)-field(i,k,j) &
3133	-u_base(k+1) +u_base(k) )
3134	ENDDO
3135	ENDDO
3136
3137	DO i = i_start, i_end
3138	vflux(i,0)=vflux(i,1)
3139	ENDDO
3140
3141	DO i = i_start, i_end
3142	vflux(i,ktf)=0.
3143	ENDDO
3144
3145	DO k=kts,ktf-1
3146	DO i = i_start, i_end
3147	tendency(i,k,j)=tendency(i,k,j)+ &
3148	ggrdnw(k)/muu(i,j)/(0.5(alt(i-1,k,j)+alt(i,k,j))) &
3149	(vflux(i,k)-vflux(i,k-1))
3150	ENDDO
3151	ENDDO
3152
3153	ENDDO j_loop_u
3154
3155	END SUBROUTINE vertical_diffusion_u
3156
3157	!-------------------------------------------------------------------------------
3158
3159
3160	SUBROUTINE vertical_diffusion_v ( field, tendency, &
3161	config_flags, v_base, &
3162	alt, muv, rdn, rdnw, kvdif, &
3163	ids, ide, jds, jde, kds, kde, &
3164	ims, ime, jms, jme, kms, kme, &
3165	its, ite, jts, jte, kts, kte )
3166
3167
3168	IMPLICIT NONE
3169
3170	! Input data
3171
3172	TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
3173
3174	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
3175	ims, ime, jms, jme, kms, kme, &
3176	its, ite, jts, jte, kts, kte
3177
3178	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , &
3179	INTENT(IN ) :: field, &
3180	alt
3181	REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: rdn, rdnw, v_base
3182
3183	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency
3184
3185	REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: muv
3186
3187	REAL , INTENT(IN ) :: kvdif
3188
3189	! Local data
3190
3191	INTEGER :: i, j, k, itf, jtf, ktf, jm1
3192	INTEGER :: i_start, i_end, j_start, j_end
3193
3194	REAL , DIMENSION(its:ite, 0:kte+1) :: vflux
3195
3196	REAL :: rdz, zz
3197
3198	LOGICAL :: specified
3199
3200	!<DESCRIPTION>
3201	!
3202	! vertical_diffusion_v computes vertical diffusion tendency for
3203	! the v momentum equation. This routine assumes a constant eddy
3204	! viscosity kvdif.
3205	!
3206	!</DESCRIPTION>
3207
3208	specified = .false.
3209	if(config_flags%specified .or. config_flags%nested) specified = .true.
3210
3211	ktf=MIN(kte,kde-1)
3212
3213	i_start = its
3214	i_end = MIN(ite,ide-1)
3215	j_start = jts
3216	j_end = MIN(jte,jde-1)
3217
3218	IF ( config_flags%open_ys .or. specified ) j_start = MAX(jds+1,jts)
3219	IF ( config_flags%open_ye .or. specified ) j_end = MIN(jde-1,jte)
3220
3221	j_loop_v : DO j = j_start, j_end
3222	! jm1 = max(j-1,1)
3223	jm1 = j-1
3224
3225	DO k=kts,ktf-1
3226	DO i = i_start, i_end
3227	vflux(i,k)=kvdifrdn(k+1)/(0.25( alt(i,k ,j ) &
3228	+alt(i,k ,jm1) &
3229	+alt(i,k+1,j ) &
3230	+alt(i,k+1,jm1) ) ) &
3231	*(field(i,k+1,j)-field(i,k,j) &
3232	-v_base(k+1) +v_base(k) )
3233	ENDDO
3234	ENDDO
3235
3236	DO i = i_start, i_end
3237	vflux(i,0)=vflux(i,1)
3238	ENDDO
3239
3240	DO i = i_start, i_end
3241	vflux(i,ktf)=0.
3242	ENDDO
3243
3244	DO k=kts,ktf-1
3245	DO i = i_start, i_end
3246	tendency(i,k,j)=tendency(i,k,j)+ &
3247	ggrdnw(k)/muv(i,j)/(0.5(alt(i,k,jm1)+alt(i,k,j))) &
3248	(vflux(i,k)-vflux(i,k-1))
3249	ENDDO
3250	ENDDO
3251
3252	ENDDO j_loop_v
3253
3254	END SUBROUTINE vertical_diffusion_v
3255
3256	!*************** end new mass coordinate routines
3257
3258	!-------------------------------------------------------------------------------
3259
3260	SUBROUTINE calculate_full ( rfield, rfieldb, rfieldp, &
3261	ids, ide, jds, jde, kds, kde, &
3262	ims, ime, jms, jme, kms, kme, &
3263	its, ite, jts, jte, kts, kte )
3264
3265	IMPLICIT NONE
3266
3267	! Input data
3268
3269	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
3270	ims, ime, jms, jme, kms, kme, &
3271	its, ite, jts, jte, kts, kte
3272
3273	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: rfieldb, &
3274	rfieldp
3275
3276	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT ) :: rfield
3277
3278	! Local indices.
3279
3280	INTEGER :: i, j, k, itf, jtf, ktf
3281
3282	!<DESCRIPTION>
3283	!
3284	! calculate_full
3285	! calculates full 3D field from pertubation and base field.
3286	!
3287	!</DESCRIPTION>
3288
3289	itf=MIN(ite,ide-1)
3290	jtf=MIN(jte,jde-1)
3291	ktf=MIN(kte,kde-1)
3292
3293	DO j=jts,jtf
3294	DO k=kts,ktf
3295	DO i=its,itf
3296	rfield(i,k,j)=rfieldb(i,k,j)+rfieldp(i,k,j)
3297	ENDDO
3298	ENDDO
3299	ENDDO
3300
3301	END SUBROUTINE calculate_full
3302
3303	!------------------------------------------------------------------------------
3304
3305	SUBROUTINE coriolis ( ru, rv, rw, ru_tend, rv_tend, rw_tend, &
3306	config_flags, &
3307	msftx, msfty, msfux, msfuy, &
3308	msfvx, msfvy, &
3309	f, e, sina, cosa, fzm, fzp, &
3310	ids, ide, jds, jde, kds, kde, &
3311	ims, ime, jms, jme, kms, kme, &
3312	its, ite, jts, jte, kts, kte )
3313
3314	IMPLICIT NONE
3315
3316	! Input data
3317
3318	TYPE(grid_config_rec_type) , INTENT(IN ) :: config_flags
3319
3320	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
3321	ims, ime, jms, jme, kms, kme, &
3322	its, ite, jts, jte, kts, kte
3323
3324	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: ru_tend, &
3325	rv_tend, &
3326	rw_tend
3327	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: ru, &
3328	rv, &
3329	rw
3330
3331	REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: msfux, &
3332	msfuy, &
3333	msfvx, &
3334	msfvy, &
3335	msftx, &
3336	msfty
3337
3338	REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: f, &
3339	e, &
3340	sina, &
3341	cosa
3342
3343	REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fzm, &
3344	fzp
3345
3346	! Local indices.
3347
3348	INTEGER :: i, j , k, ktf
3349	INTEGER :: i_start, i_end, j_start, j_end
3350
3351	LOGICAL :: specified
3352
3353	!<DESCRIPTION>
3354	!
3355	! coriolis calculates the large timestep tendency terms in the
3356	! u, v, and w momentum equations arise from the coriolis force.
3357	!
3358	!</DESCRIPTION>
3359
3360	specified = .false.
3361	if(config_flags%specified .or. config_flags%nested) specified = .true.
3362
3363	ktf=MIN(kte,kde-1)
3364
3365	! coriolis for u-momentum equation
3366
3367	! Notes on map scale factor
3368	! cosa, sina are related to rotating the coordinate frame if desired
3369	! generally sina=0, cosa=1
3370	! ADT eqn 44, RHS terms 6 and 7: -2 mu w omega cos(lat)/my
3371	! + 2 mu v omega sin(lat)/my
3372	! Define f=2 omega sin(lat), e=2 omega cos(lat)
3373	! => terms are: -e mu w / my + f mu v / my
3374	! rv = mu v / mx ; rw = mu w / my
3375	! => terms are: -e rw + f rv *mx / my
3376
3377	i_start = its
3378	i_end = ite
3379	IF ( config_flags%open_xs .or. specified .or. &
3380	config_flags%nested) i_start = MAX(ids+1,its)
3381	IF ( config_flags%open_xe .or. specified .or. &
3382	config_flags%nested) i_end = MIN(ide-1,ite)
3383	IF ( config_flags%periodic_x ) i_start = its
3384	IF ( config_flags%periodic_x ) i_end = ite
3385
3386	DO j = jts, MIN(jte,jde-1)
3387
3388	DO k=kts,ktf
3389	DO i = i_start, i_end
3390
3391	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)) &
3392	0.25(rv(i-1,k,j+1)+rv(i,k,j+1)+rv(i-1,k,j)+rv(i,k,j)) &
3393	- 0.5(e(i,j)+e(i-1,j))0.5*(cosa(i,j)+cosa(i-1,j)) &
3394	0.25(rw(i-1,k+1,j)+rw(i-1,k,j)+rw(i,k+1,j)+rw(i,k,j))
3395
3396	ENDDO
3397	ENDDO
3398
3399	IF ( (config_flags%open_xs) .and. (its == ids) ) THEN
3400
3401	DO k=kts,ktf
3402
3403	ru_tend(its,k,j)=ru_tend(its,k,j) + (msfux(its,j)/msfuy(its,j))0.5(f(its,j)+f(its,j)) &
3404	0.25(rv(its,k,j+1)+rv(its,k,j+1)+rv(its,k,j)+rv(its,k,j)) &
3405	- 0.5(e(its,j)+e(its,j))0.5*(cosa(its,j)+cosa(its,j)) &
3406	0.25(rw(its,k+1,j)+rw(its,k,j)+rw(its,k+1,j)+rw(its,k,j))
3407
3408	ENDDO
3409
3410	ENDIF
3411
3412	IF ( (config_flags%open_xe) .and. (ite == ide) ) THEN
3413
3414	DO k=kts,ktf
3415
3416	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)) &
3417	0.25(rv(ite-1,k,j+1)+rv(ite-1,k,j+1)+rv(ite-1,k,j)+rv(ite-1,k,j)) &
3418	- 0.5(e(ite-1,j)+e(ite-1,j))0.5*(cosa(ite-1,j)+cosa(ite-1,j)) &
3419	0.25(rw(ite-1,k+1,j)+rw(ite-1,k,j)+rw(ite-1,k+1,j)+rw(ite-1,k,j))
3420
3421	ENDDO
3422
3423	ENDIF
3424
3425	ENDDO
3426
3427	! coriolis term for v-momentum equation
3428
3429	! Notes on map scale factors
3430	! ADT eqn 45, RHS terms 6 and 6b [0 for sina=0]: -2 mu u omega sin(lat)/mx + ?
3431	! Define f=2 omega sin(lat), e=2 omega cos(lat)
3432	! => terms are: -f mu u / mx
3433	! ru = mu u / my ; rw = mu w / my
3434	! => terms are: -f ru *my / mx + ?
3435
3436	j_start = jts
3437	j_end = jte
3438
3439	IF ( config_flags%open_ys .or. specified .or. &
3440	config_flags%nested .or. config_flags%polar) j_start = MAX(jds+1,jts)
3441	IF ( config_flags%open_ye .or. specified .or. &
3442	config_flags%nested .or. config_flags%polar) j_end = MIN(jde-1,jte)
3443
3444	IF ( (config_flags%open_ys) .and. (jts == jds) ) THEN
3445
3446	DO k=kts,ktf
3447	DO i=its,MIN(ide-1,ite)
3448
3449	rv_tend(i,k,jts)=rv_tend(i,k,jts) - (msfvy(i,jts)/msfvx(i,jts))0.5(f(i,jts)+f(i,jts)) &
3450	0.25(ru(i,k,jts)+ru(i+1,k,jts)+ru(i,k,jts)+ru(i+1,k,jts)) &
3451	+ (msfvy(i,jts)/msfvx(i,jts))0.5(e(i,jts)+e(i,jts))0.5(sina(i,jts)+sina(i,jts)) &
3452	0.25(rw(i,k+1,jts)+rw(i,k,jts)+rw(i,k+1,jts)+rw(i,k,jts))
3453
3454	ENDDO
3455	ENDDO
3456
3457	ENDIF
3458
3459	DO j=j_start, j_end
3460	DO k=kts,ktf
3461	DO i=its,MIN(ide-1,ite)
3462
3463	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)) &
3464	0.25(ru(i,k,j)+ru(i+1,k,j)+ru(i,k,j-1)+ru(i+1,k,j-1)) &
3465	+ (msfvy(i,j)/msfvx(i,j))0.5(e(i,j)+e(i,j-1))0.5(sina(i,j)+sina(i,j-1)) &
3466	0.25(rw(i,k+1,j-1)+rw(i,k,j-1)+rw(i,k+1,j)+rw(i,k,j))
3467
3468	ENDDO
3469	ENDDO
3470	ENDDO
3471
3472
3473	IF ( (config_flags%open_ye) .and. (jte == jde) ) THEN
3474
3475	DO k=kts,ktf
3476	DO i=its,MIN(ide-1,ite)
3477
3478	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)) &
3479	0.25(ru(i,k,jte-1)+ru(i+1,k,jte-1)+ru(i,k,jte-1)+ru(i+1,k,jte-1)) &
3480	+ (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)) &
3481	0.25(rw(i,k+1,jte-1)+rw(i,k,jte-1)+rw(i,k+1,jte-1)+rw(i,k,jte-1))
3482
3483	ENDDO
3484	ENDDO
3485
3486	ENDIF
3487
3488	! coriolis term for w-mometum
3489
3490	! Notes on map scale factors
3491	! ADT eqn 46/my, RHS terms 5 and 5b [0 for sina=0]: 2 mu u omega cos(lat)/my +?
3492	! Define e=2 omega cos(lat)
3493	! => terms are: e mu u / my + ???
3494	! ru = mu u / my ; ru = mu v / mx
3495	! => terms are: e ru + ???
3496
3497	DO j=jts,MIN(jte, jde-1)
3498	DO k=kts+1,ktf
3499	DO i=its,MIN(ite, ide-1)
3500
3501	rw_tend(i,k,j)=rw_tend(i,k,j) + e(i,j)* &
3502	(cosa(i,j)0.5(fzm(k)*(ru(i,k,j)+ru(i+1,k,j)) &
3503	+fzp(k)*(ru(i,k-1,j)+ru(i+1,k-1,j))) &
3504	-(msftx(i,j)/msfty(i,j))* &
3505	sina(i,j)0.5(fzm(k)*(rv(i,k,j)+rv(i,k,j+1)) &
3506	+fzp(k)*(rv(i,k-1,j)+rv(i,k-1,j+1))))
3507
3508	ENDDO
3509	ENDDO
3510	ENDDO
3511
3512	END SUBROUTINE coriolis
3513
3514	!------------------------------------------------------------------------------
3515
3516	SUBROUTINE perturbation_coriolis ( ru_in, rv_in, rw, ru_tend, rv_tend, rw_tend, &
3517	config_flags, &
3518	u_base, v_base, z_base, &
3519	muu, muv, phb, ph, &
3520	msftx, msfty, msfux, msfuy, msfvx, msfvy, &
3521	f, e, sina, cosa, fzm, fzp, &
3522	ids, ide, jds, jde, kds, kde, &
3523	ims, ime, jms, jme, kms, kme, &
3524	its, ite, jts, jte, kts, kte )
3525
3526	IMPLICIT NONE
3527
3528	! Input data
3529
3530	TYPE(grid_config_rec_type) , INTENT(IN ) :: config_flags
3531
3532	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
3533	ims, ime, jms, jme, kms, kme, &
3534	its, ite, jts, jte, kts, kte
3535
3536	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: ru_tend, &
3537	rv_tend, &
3538	rw_tend
3539	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: ru_in, &
3540	rv_in, &
3541	rw, &
3542	ph, &
3543	phb
3544
3545
3546	REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: msfux, &
3547	msfuy, &
3548	msfvx, &
3549	msfvy, &
3550	msftx, &
3551	msfty
3552
3553	REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: f, &
3554	e, &
3555	sina, &
3556	cosa
3557
3558	REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: muu, &
3559	muv
3560
3561
3562	REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fzm, &
3563	fzp
3564
3565	REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: u_base, &
3566	v_base, &
3567	z_base
3568
3569	! Local storage
3570
3571	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) :: ru, &
3572	rv
3573
3574	REAL :: z_at_u, z_at_v, wkp1, wk, wkm1
3575
3576	! Local indices.
3577
3578	INTEGER :: i, j , k, ktf
3579	INTEGER :: i_start, i_end, j_start, j_end
3580
3581	LOGICAL :: specified
3582
3583	!<DESCRIPTION>
3584	!
3585	! perturbation_coriolis calculates the large timestep tendency terms in the
3586	! u, v, and w momentum equations arise from the coriolis force. This version
3587	! subtracts off the horizontal velocities from the initial sounding when
3588	! computing the forcing terms, hence "perturbation" coriolis.
3589	!
3590	!</DESCRIPTION>
3591
3592	specified = .false.
3593	if(config_flags%specified .or. config_flags%nested) specified = .true.
3594
3595	ktf=MIN(kte,kde-1)
3596
3597	! coriolis for u-momentum equation
3598
3599	i_start = its
3600	i_end = ite
3601	IF ( config_flags%open_xs .or. specified .or. &
3602	config_flags%nested) i_start = MAX(ids+1,its)
3603	IF ( config_flags%open_xe .or. specified .or. &
3604	config_flags%nested) i_end = MIN(ide-1,ite)
3605	IF ( config_flags%periodic_x ) i_start = its
3606	IF ( config_flags%periodic_x ) i_end = ite
3607
3608	! compute perturbation mu*v for use in u momentum equation
3609
3610	DO j = jts, MIN(jte,jde-1)+1
3611	DO k=kts+1,ktf-1
3612	DO i = i_start-1, i_end
3613	z_at_v = 0.25*( phb(i,k,j )+phb(i,k+1,j ) &
3614	+phb(i,k,j-1)+phb(i,k+1,j-1) &
3615	+ph(i,k,j )+ph(i,k+1,j ) &
3616	+ph(i,k,j-1)+ph(i,k+1,j-1))/g
3617	wkp1 = min(1.,max(0.,z_at_v-z_base(k))/(z_base(k+1)-z_base(k)))
3618	wkm1 = min(1.,max(0.,z_base(k)-z_at_v)/(z_base(k)-z_base(k-1)))
3619	wk = 1.-wkp1-wkm1
3620	rv(i,k,j) = rv_in(i,k,j) - muv(i,j)*( &
3621	wkm1*v_base(k-1) &
3622	+wk *v_base(k ) &
3623	+wkp1*v_base(k+1) )
3624	ENDDO
3625	ENDDO
3626	ENDDO
3627
3628
3629	! pick up top and bottom v
3630
3631	DO j = jts, MIN(jte,jde-1)+1
3632	DO i = i_start-1, i_end
3633
3634	k = kts
3635	z_at_v = 0.25*( phb(i,k,j )+phb(i,k+1,j ) &
3636	+phb(i,k,j-1)+phb(i,k+1,j-1) &
3637	+ph(i,k,j )+ph(i,k+1,j ) &
3638	+ph(i,k,j-1)+ph(i,k+1,j-1))/g
3639	wkp1 = min(1.,max(0.,z_at_v-z_base(k))/(z_base(k+1)-z_base(k)))
3640	wk = 1.-wkp1
3641	rv(i,k,j) = rv_in(i,k,j) - muv(i,j)*( &
3642	+wk *v_base(k ) &
3643	+wkp1*v_base(k+1) )
3644
3645	k = ktf
3646	z_at_v = 0.25*( phb(i,k,j )+phb(i,k+1,j ) &
3647	+phb(i,k,j-1)+phb(i,k+1,j-1) &
3648	+ph(i,k,j )+ph(i,k+1,j ) &
3649	+ph(i,k,j-1)+ph(i,k+1,j-1))/g
3650	wkm1 = min(1.,max(0.,z_base(k)-z_at_v)/(z_base(k)-z_base(k-1)))
3651	wk = 1.-wkm1
3652	rv(i,k,j) = rv_in(i,k,j) - muv(i,j)*( &
3653	wkm1*v_base(k-1) &
3654	+wk *v_base(k ) )
3655
3656	ENDDO
3657	ENDDO
3658
3659	! compute coriolis forcing for u
3660
3661	! Map scale factors: see comments above for Coriolis
3662
3663	DO j = jts, MIN(jte,jde-1)
3664
3665	DO k=kts,ktf
3666	DO i = i_start, i_end
3667	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)) &
3668	0.25(rv(i-1,k,j+1)+rv(i,k,j+1)+rv(i-1,k,j)+rv(i,k,j)) &
3669	- 0.5(e(i,j)+e(i-1,j))0.5*(cosa(i,j)+cosa(i-1,j)) &
3670	0.25(rw(i-1,k+1,j)+rw(i-1,k,j)+rw(i,k+1,j)+rw(i,k,j))
3671	ENDDO
3672	ENDDO
3673
3674	IF ( (config_flags%open_xs) .and. (its == ids) ) THEN
3675
3676	DO k=kts,ktf
3677
3678	ru_tend(its,k,j)=ru_tend(its,k,j) + (msfux(its,j)/msfuy(its,j))0.5(f(its,j)+f(its,j)) &
3679	0.25(rv(its,k,j+1)+rv(its,k,j+1)+rv(its,k,j)+rv(its,k,j)) &
3680	- 0.5(e(its,j)+e(its,j))0.5*(cosa(its,j)+cosa(its,j)) &
3681	0.25(rw(its,k+1,j)+rw(its,k,j)+rw(its,k+1,j)+rw(its,k,j))
3682
3683	ENDDO
3684
3685	ENDIF
3686
3687	IF ( (config_flags%open_xe) .and. (ite == ide) ) THEN
3688
3689	DO k=kts,ktf
3690
3691	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)) &
3692	0.25(rv(ite-1,k,j+1)+rv(ite-1,k,j+1)+rv(ite-1,k,j)+rv(ite-1,k,j)) &
3693	- 0.5(e(ite-1,j)+e(ite-1,j))0.5*(cosa(ite-1,j)+cosa(ite-1,j)) &
3694	0.25(rw(ite-1,k+1,j)+rw(ite-1,k,j)+rw(ite-1,k+1,j)+rw(ite-1,k,j))
3695
3696	ENDDO
3697
3698	ENDIF
3699
3700	ENDDO
3701
3702	! coriolis term for v-momentum equation
3703	! Map scale factors: see comments above for Coriolis
3704
3705	j_start = jts
3706	j_end = jte
3707
3708	IF ( config_flags%open_ys .or. specified .or. &
3709	config_flags%nested .or. config_flags%polar) j_start = MAX(jds+1,jts)
3710	IF ( config_flags%open_ye .or. specified .or. &
3711	config_flags%nested .or. config_flags%polar) j_end = MIN(jde-1,jte)
3712
3713	! compute perturbation mu*u for use in v momentum equation
3714
3715	DO j = j_start-1,j_end
3716	DO k=kts+1,ktf-1
3717	DO i = its, MIN(ite,ide-1)+1
3718	z_at_u = 0.25*( phb(i ,k,j)+phb(i ,k+1,j) &
3719	+phb(i-1,k,j)+phb(i-1,k+1,j) &
3720	+ph(i ,k,j)+ph(i ,k+1,j) &
3721	+ph(i-1,k,j)+ph(i-1,k+1,j))/g
3722	wkp1 = min(1.,max(0.,z_at_u-z_base(k))/(z_base(k+1)-z_base(k)))
3723	wkm1 = min(1.,max(0.,z_base(k)-z_at_u)/(z_base(k)-z_base(k-1)))
3724	wk = 1.-wkp1-wkm1
3725	ru(i,k,j) = ru_in(i,k,j) - muu(i,j)*( &
3726	wkm1*u_base(k-1) &
3727	+wk *u_base(k ) &
3728	+wkp1*u_base(k+1) )
3729	ENDDO
3730	ENDDO
3731	ENDDO
3732
3733	! pick up top and bottom u
3734
3735	DO j = j_start-1,j_end
3736	DO i = its, MIN(ite,ide-1)+1
3737
3738	k = kts
3739	z_at_u = 0.25*( phb(i ,k,j)+phb(i ,k+1,j) &
3740	+phb(i-1,k,j)+phb(i-1,k+1,j) &
3741	+ph(i ,k,j)+ph(i ,k+1,j) &
3742	+ph(i-1,k,j)+ph(i-1,k+1,j))/g
3743	wkp1 = min(1.,max(0.,z_at_u-z_base(k))/(z_base(k+1)-z_base(k)))
3744	wk = 1.-wkp1
3745	ru(i,k,j) = ru_in(i,k,j) - muu(i,j)*( &
3746	+wk *u_base(k ) &
3747	+wkp1*u_base(k+1) )
3748
3749
3750	k = ktf
3751	z_at_u = 0.25*( phb(i ,k,j)+phb(i ,k+1,j) &
3752	+phb(i-1,k,j)+phb(i-1,k+1,j) &
3753	+ph(i ,k,j)+ph(i ,k+1,j) &
3754	+ph(i-1,k,j)+ph(i-1,k+1,j))/g
3755	wkm1 = min(1.,max(0.,z_base(k)-z_at_u)/(z_base(k)-z_base(k-1)))
3756	wk = 1.-wkm1
3757	ru(i,k,j) = ru_in(i,k,j) - muu(i,j)*( &
3758	wkm1*u_base(k-1) &
3759	+wk *u_base(k ) )
3760
3761	ENDDO
3762	ENDDO
3763
3764	! compute coriolis forcing for v momentum equation
3765	! Map scale factors: see comments above for Coriolis
3766
3767	IF ( (config_flags%open_ys) .and. (jts == jds) ) THEN
3768
3769	DO k=kts,ktf
3770	DO i=its,MIN(ide-1,ite)
3771
3772	rv_tend(i,k,jts)=rv_tend(i,k,jts) - (msfvy(i,jts)/msfvx(i,jts))0.5(f(i,jts)+f(i,jts)) &
3773	0.25(ru(i,k,jts)+ru(i+1,k,jts)+ru(i,k,jts)+ru(i+1,k,jts)) &
3774	+ (msfvy(i,jts)/msfvx(i,jts))0.5(e(i,jts)+e(i,jts))0.5(sina(i,jts)+sina(i,jts)) &
3775	0.25(rw(i,k+1,jts)+rw(i,k,jts)+rw(i,k+1,jts)+rw(i,k,jts))
3776
3777	ENDDO
3778	ENDDO
3779
3780	ENDIF
3781
3782	DO j=j_start, j_end
3783	DO k=kts,ktf
3784	DO i=its,MIN(ide-1,ite)
3785
3786	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)) &
3787	0.25(ru(i,k,j)+ru(i+1,k,j)+ru(i,k,j-1)+ru(i+1,k,j-1)) &
3788	+ (msfvy(i,j)/msfvx(i,j))0.5(e(i,j)+e(i,j-1))0.5(sina(i,j)+sina(i,j-1)) &
3789	0.25(rw(i,k+1,j-1)+rw(i,k,j-1)+rw(i,k+1,j)+rw(i,k,j))
3790
3791	ENDDO
3792	ENDDO
3793	ENDDO
3794
3795
3796	IF ( (config_flags%open_ye) .and. (jte == jde) ) THEN
3797
3798	DO k=kts,ktf
3799	DO i=its,MIN(ide-1,ite)
3800
3801	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)) &
3802	0.25(ru(i,k,jte-1)+ru(i+1,k,jte-1)+ru(i,k,jte-1)+ru(i+1,k,jte-1)) &
3803	+ (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)) &
3804	0.25(rw(i,k+1,jte-1)+rw(i,k,jte-1)+rw(i,k+1,jte-1)+rw(i,k,jte-1))
3805
3806	ENDDO
3807	ENDDO
3808
3809	ENDIF
3810
3811	! coriolis term for w-mometum
3812	! Map scale factors: see comments above for Coriolis
3813
3814	DO j=jts,MIN(jte, jde-1)
3815	DO k=kts+1,ktf
3816	DO i=its,MIN(ite, ide-1)
3817
3818	rw_tend(i,k,j)=rw_tend(i,k,j) + e(i,j)* &
3819	(cosa(i,j)0.5(fzm(k)*(ru(i,k,j)+ru(i+1,k,j)) &
3820	+fzp(k)*(ru(i,k-1,j)+ru(i+1,k-1,j))) &
3821	-(msftx(i,j)/msfty(i,j))sina(i,j)0.5(fzm(k)(rv(i,k,j)+rv(i,k,j+1)) &
3822	+fzp(k)*(rv(i,k-1,j)+rv(i,k-1,j+1))))
3823
3824	ENDDO
3825	ENDDO
3826	ENDDO
3827
3828	END SUBROUTINE perturbation_coriolis
3829
3830	!------------------------------------------------------------------------------
3831
3832	SUBROUTINE curvature ( ru, rv, rw, u, v, w, ru_tend, rv_tend, rw_tend, &
3833	config_flags, &
3834	msfux, msfuy, msfvx, msfvy, msftx, msfty, &
3835	xlat, fzm, fzp, rdx, rdy, &
3836	ids, ide, jds, jde, kds, kde, &
3837	ims, ime, jms, jme, kms, kme, &
3838	its, ite, jts, jte, kts, kte )
3839
3840
3841	IMPLICIT NONE
3842
3843	! Input data
3844
3845	TYPE(grid_config_rec_type) , INTENT(IN ) :: config_flags
3846
3847	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
3848	ims, ime, jms, jme, kms, kme, &
3849	its, ite, jts, jte, kts, kte
3850
3851	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , &
3852	INTENT(INOUT) :: ru_tend, &
3853	rv_tend, &
3854	rw_tend
3855
3856	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , &
3857	INTENT(IN ) :: ru, &
3858	rv, &
3859	rw, &
3860	u, &
3861	v, &
3862	w
3863
3864	REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: msfux, &
3865	msfuy, &
3866	msfvx, &
3867	msfvy, &
3868	msftx, &
3869	msfty, &
3870	xlat
3871
3872	REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fzm, &
3873	fzp
3874
3875	REAL , INTENT(IN ) :: rdx, &
3876	rdy
3877
3878	! Local data
3879
3880	! INTEGER :: i, j, k, itf, jtf, ktf, kp1, im, ip, jm, jp
3881	INTEGER :: i, j, k, itf, jtf, ktf
3882	INTEGER :: i_start, i_end, j_start, j_end
3883	! INTEGER :: irmin, irmax, jrmin, jrmax
3884
3885	REAL , DIMENSION( its-1:ite , kts:kte, jts-1:jte ) :: vxgm
3886
3887	LOGICAL :: specified
3888
3889	!<DESCRIPTION>
3890	!
3891	! curvature calculates the large timestep tendency terms in the
3892	! u, v, and w momentum equations arise from the curvature terms.
3893	!
3894	!</DESCRIPTION>
3895
3896	specified = .false.
3897	if(config_flags%specified .or. config_flags%nested) specified = .true.
3898
3899	itf=MIN(ite,ide-1)
3900	jtf=MIN(jte,jde-1)
3901	ktf=MIN(kte,kde-1)
3902
3903	! irmin = ims
3904	! irmax = ime
3905	! jrmin = jms
3906	! jrmax = jme
3907	! IF ( config_flags%open_xs ) irmin = ids
3908	! IF ( config_flags%open_xe ) irmax = ide-1
3909	! IF ( config_flags%open_ys ) jrmin = jds
3910	! IF ( config_flags%open_ye ) jrmax = jde-1
3911
3912	! Define v cross grad m at scalar points - vxgm(i,j)
3913
3914	i_start = its-1
3915	i_end = ite
3916	j_start = jts-1
3917	j_end = jte
3918
3919	IF ( ( config_flags%open_xs .or. specified .or. &
3920	config_flags%nested) .and. (its == ids) ) i_start = its
3921	IF ( ( config_flags%open_xe .or. specified .or. &
3922	config_flags%nested) .and. (ite == ide) ) i_end = ite-1
3923	IF ( ( config_flags%open_ys .or. specified .or. &
3924	config_flags%nested .or. config_flags%polar) .and. (jts == jds) ) j_start = jts
3925	IF ( ( config_flags%open_ye .or. specified .or. &
3926	config_flags%nested .or. config_flags%polar) .and. (jte == jde) ) j_end = jte-1
3927	IF ( config_flags%periodic_x ) i_start = its-1
3928	IF ( config_flags%periodic_x ) i_end = ite
3929
3930	DO j=j_start, j_end
3931	DO k=kts,ktf
3932	DO i=i_start, i_end
3933	! Map scale factor notes:
3934	! msf...y is constant everywhere for cylindrical map projection
3935	! msf...x varies with y only
3936	! But we know that this is not = 0 for cylindrical,
3937	! therefore use msfvX in 1st line
3938	! which => by symmetry use msfuY in 2nd line - ???
3939	vxgm(i,k,j)=0.5(u(i,k,j)+u(i+1,k,j))(msfvx(i,j+1)-msfvx(i,j))*rdy - &
3940	0.5(v(i,k,j)+v(i,k,j+1))(msfuy(i+1,j)-msfuy(i,j))*rdx
3941	ENDDO
3942	ENDDO
3943	ENDDO
3944
3945	! Pick up the boundary rows for open (radiation) lateral b.c.
3946	! Rather crude at present, we are assuming there is no
3947	! variation in this term at the boundary.
3948
3949	IF ( ( config_flags%open_xs .or. (specified .AND. .NOT. config_flags%periodic_x) .or. &
3950	config_flags%nested) .and. (its == ids) ) THEN
3951
3952	DO j = jts, jte-1
3953	DO k = kts, ktf
3954	vxgm(its-1,k,j) = vxgm(its,k,j)
3955	ENDDO
3956	ENDDO
3957
3958	ENDIF
3959
3960	IF ( ( config_flags%open_xe .or. (specified .AND. .NOT. config_flags%periodic_x) .or. &
3961	config_flags%nested) .and. (ite == ide) ) THEN
3962
3963	DO j = jts, jte-1
3964	DO k = kts, ktf
3965	vxgm(ite,k,j) = vxgm(ite-1,k,j)
3966	ENDDO
3967	ENDDO
3968
3969	ENDIF
3970
3971	! Polar boundary condition:
3972	! The following change is needed in case one tries using the vxgm route with
3973	! polar B.C.'s in the future, but not needed if 'tan' used
3974	IF ( ( config_flags%open_ys .or. specified .or. &
3975	config_flags%nested .or. config_flags%polar) .and. (jts == jds) ) THEN
3976
3977	DO k = kts, ktf
3978	DO i = its-1, ite
3979	vxgm(i,k,jts-1) = vxgm(i,k,jts)
3980	ENDDO
3981	ENDDO
3982
3983	ENDIF
3984
3985	! Polar boundary condition:
3986	! The following change is needed in case one tries using the vxgm route with
3987	! polar B.C.'s in the future, but not needed if 'tan' used
3988	IF ( ( config_flags%open_ye .or. specified .or. &
3989	config_flags%nested .or. config_flags%polar) .and. (jte == jde) ) THEN
3990
3991	DO k = kts, ktf
3992	DO i = its-1, ite
3993	vxgm(i,k,jte) = vxgm(i,k,jte-1)
3994	ENDDO
3995	ENDDO
3996
3997	ENDIF
3998
3999	! curvature term for u momentum eqn.
4000
4001	! Map scale factor notes:
4002	! ADT eqn 44, RHS terms 4 and 5, in cylindrical: mu u v tan(lat)/(a my)
4003	! - mu u w /(a my)
4004	! ru = mu u / my ; rw = mu w / my ; rv = mu v / mx
4005	! => terms are:
4006	! (mx/my)u rv tan(lat) / a - u rw / a = (u/a)[(mx/my) rv tan(lat) - rw]
4007	! ru v tan(lat) / a - u rw / a
4008	! xlat defined with end points half grid space from pole,
4009	! hence are on u latitude points
4010
4011	i_start = its
4012	IF ( config_flags%open_xs .or. specified .or. &
4013	config_flags%nested) i_start = MAX ( ids+1 , its )
4014	IF ( config_flags%open_xe .or. specified .or. &
4015	config_flags%nested) i_end = MIN ( ide-1 , ite )
4016	IF ( config_flags%periodic_x ) i_start = its
4017	IF ( config_flags%periodic_x ) i_end = ite
4018
4019	! Polar boundary condition
4020	IF ((config_flags%map_proj == 6) .OR. (config_flags%polar)) THEN
4021
4022	DO j=jts,MIN(jde-1,jte)
4023	DO k=kts,ktf
4024	DO i=i_start,i_end
4025
4026	ru_tend(i,k,j)=ru_tend(i,k,j) + u(i,k,j)reradius ( &
4027	(msfux(i,j)/msfuy(i,j))0.25(rv(i-1,k,j+1)+rv(i,k,j+1)+ &
4028	rv(i-1,k,j)+rv(i,k,j))tan(xlat(i,j)degrad) &
4029	- 0.25*(rw(i-1,k+1,j)+rw(i-1,k,j)+rw(i,k+1,j)+rw(i,k,j)) )
4030	ENDDO
4031	ENDDO
4032	ENDDO
4033
4034	ELSE ! normal code
4035
4036
4037	DO j=jts,MIN(jde-1,jte)
4038	DO k=kts,ktf
4039	DO i=i_start,i_end
4040
4041	ru_tend(i,k,j)=ru_tend(i,k,j) + 0.5*(vxgm(i,k,j)+vxgm(i-1,k,j)) &
4042	0.25(rv(i-1,k,j+1)+rv(i,k,j+1)+rv(i-1,k,j)+rv(i,k,j)) &
4043	- u(i,k,j)*reradius &
4044	0.25(rw(i-1,k+1,j)+rw(i-1,k,j)+rw(i,k+1,j)+rw(i,k,j))
4045
4046	ENDDO
4047	ENDDO
4048	ENDDO
4049
4050	END IF
4051
4052	! curvature term for v momentum eqn.
4053
4054	! Map scale factor notes
4055	! ADT eqn 45, RHS terms 4 and 5, in cylindrical: mu u*u tan(lat)/(a mx)
4056	! - mu v w /(a mx)
4057	! ru = mu u / my ; rw = mu w / my ; rv = mu v / mx
4058	! terms are:
4059	! (my/mx)u ru tan(lat) / a - (my/mx)v rw / a
4060	! = [my/(mxa)][u ru tan(lat) - v rw]
4061	! (1/a)[(my/mx)u ru tan(lat) - w rv]
4062	! xlat defined with end points half grid space from pole, hence are on
4063	! u latitude points => av here
4064	!
4065	! in original wrf, there was a sign error for the rw contribution
4066
4067	j_start = jts
4068	IF ( config_flags%open_ys .or. specified .or. &
4069	config_flags%nested .or. config_flags%polar) j_start = MAX ( jds+1 , jts )
4070	IF ( config_flags%open_ye .or. specified .or. &
4071	config_flags%nested .or. config_flags%polar) j_end = MIN ( jde-1 , jte )
4072
4073	IF ((config_flags%map_proj == 6) .OR. (config_flags%polar)) THEN
4074
4075	DO j=j_start,j_end
4076	DO k=kts,ktf
4077	DO i=its,MIN(ite,ide-1)
4078	rv_tend(i,k,j)=rv_tend(i,k,j) - (msfvy(i,j)/msfvx(i,j))reradius ( &
4079	0.25(u(i,k,j)+u(i+1,k,j)+u(i,k,j-1)+u(i+1,k,j-1)) &
4080	tan((xlat(i,j)+xlat(i,j-1))0.5degrad)* &
4081	0.25*(ru(i,k,j)+ru(i+1,k,j)+ru(i,k,j-1)+ru(i+1,k,j-1)) &
4082	- v(i,k,j)0.25(rw(i,k+1,j-1)+rw(i,k,j-1)+ &
4083	rw(i,k+1,j)+rw(i,k,j)) )
4084	ENDDO
4085	ENDDO
4086	ENDDO
4087
4088	ELSE ! normal code
4089
4090	DO j=j_start,j_end
4091	DO k=kts,ktf
4092	DO i=its,MIN(ite,ide-1)
4093
4094	rv_tend(i,k,j)=rv_tend(i,k,j) - 0.5*(vxgm(i,k,j)+vxgm(i,k,j-1)) &
4095	0.25(ru(i,k,j)+ru(i+1,k,j)+ru(i,k,j-1)+ru(i+1,k,j-1)) &
4096	- (msfvy(i,j)/msfvx(i,j))v(i,k,j)reradius &
4097	0.25(rw(i,k+1,j-1)+rw(i,k,j-1)+rw(i,k+1,j)+rw(i,k,j))
4098
4099	ENDDO
4100	ENDDO
4101	ENDDO
4102
4103	END IF
4104
4105	! curvature term for vertical momentum eqn.
4106
4107	! Notes on map scale factors:
4108	! ADT eqn 46, RHS term 4: [mu/(a my)][uu + v*v]
4109	! ru = mu u / my ; rw = mu w / my ; rv = mu v / mx
4110	! terms are: u ru / a + (mx/my)v rv / a
4111
4112	DO j=jts,MIN(jte,jde-1)
4113	DO k=MAX(2,kts),ktf
4114	DO i=its,MIN(ite,ide-1)
4115
4116	rw_tend(i,k,j)=rw_tend(i,k,j) + reradius* &
4117	(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))) &
4118	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))) &
4119	+(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))) &
4120	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))))
4121
4122	ENDDO
4123	ENDDO
4124	ENDDO
4125
4126	END SUBROUTINE curvature
4127
4128	!------------------------------------------------------------------------------
4129
4130	SUBROUTINE decouple ( rr, rfield, field, name, config_flags, &
4131	fzm, fzp, &
4132	ids, ide, jds, jde, kds, kde, &
4133	ims, ime, jms, jme, kms, kme, &
4134	its, ite, jts, jte, kts, kte )
4135
4136	IMPLICIT NONE
4137
4138	! Input data
4139
4140	TYPE(grid_config_rec_type) , INTENT(IN ) :: config_flags
4141
4142	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
4143	ims, ime, jms, jme, kms, kme, &
4144	its, ite, jts, jte, kts, kte
4145
4146	CHARACTER(LEN=1) , INTENT(IN ) :: name
4147
4148	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: rfield
4149
4150	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: rr
4151
4152	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT ) :: field
4153
4154	REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fzm, fzp
4155
4156	! Local data
4157
4158	INTEGER :: i, j, k, itf, jtf, ktf
4159
4160	!<DESCRIPTION>
4161	!
4162	! decouple decouples a variable from the column dry-air mass.
4163	!
4164	!</DESCRIPTION>
4165
4166	ktf=MIN(kte,kde-1)
4167
4168	IF (name .EQ. 'u')THEN
4169	itf=ite
4170	jtf=MIN(jte,jde-1)
4171
4172	DO j=jts,jtf
4173	DO k=kts,ktf
4174	DO i=its,itf
4175	field(i,k,j)=rfield(i,k,j)/(0.5*(rr(i,k,j)+rr(i-1,k,j)))
4176	ENDDO
4177	ENDDO
4178	ENDDO
4179
4180	ELSE IF (name .EQ. 'v')THEN
4181	itf=MIN(ite,ide-1)
4182	jtf=jte
4183
4184	DO j=jts,jtf
4185	DO k=kts,ktf
4186	DO i=its,itf
4187	field(i,k,j)=rfield(i,k,j)/(0.5*(rr(i,k,j)+rr(i,k,j-1)))
4188	ENDDO
4189	ENDDO
4190	ENDDO
4191
4192	ELSE IF (name .EQ. 'w')THEN
4193	itf=MIN(ite,ide-1)
4194	jtf=MIN(jte,jde-1)
4195	DO j=jts,jtf
4196	DO k=kts+1,ktf
4197	DO i=its,itf
4198	field(i,k,j)=rfield(i,k,j)/(fzm(k)rr(i,k,j)+fzp(k)rr(i,k-1,j))
4199	ENDDO
4200	ENDDO
4201	ENDDO
4202
4203	DO j=jts,jtf
4204	DO i=its,itf
4205	field(i,kte,j) = 0.
4206	ENDDO
4207	ENDDO
4208
4209	ELSE
4210	itf=MIN(ite,ide-1)
4211	jtf=MIN(jte,jde-1)
4212	! For theta we will decouple tb and tp and add them to give t afterwards
4213	DO j=jts,jtf
4214	DO k=kts,ktf
4215	DO i=its,itf
4216	field(i,k,j)=rfield(i,k,j)/rr(i,k,j)
4217	ENDDO
4218	ENDDO
4219	ENDDO
4220
4221	ENDIF
4222
4223	END SUBROUTINE decouple
4224
4225	!-------------------------------------------------------------------------------
4226
4227
4228	SUBROUTINE zero_tend ( tendency, &
4229	ids, ide, jds, jde, kds, kde, &
4230	ims, ime, jms, jme, kms, kme, &
4231	its, ite, jts, jte, kts, kte )
4232
4233
4234	IMPLICIT NONE
4235
4236	! Input data
4237
4238	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
4239	ims, ime, jms, jme, kms, kme, &
4240	its, ite, jts, jte, kts, kte
4241
4242	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency
4243
4244	! Local data
4245
4246	INTEGER :: i, j, k, itf, jtf, ktf
4247
4248	!<DESCRIPTION>
4249	!
4250	! zero_tend sets the input tendency array to zero.
4251	!
4252	!</DESCRIPTION>
4253
4254	DO j = jts, jte
4255	DO k = kts, kte
4256	DO i = its, ite
4257	tendency(i,k,j) = 0.
4258	ENDDO
4259	ENDDO
4260	ENDDO
4261
4262	END SUBROUTINE zero_tend
4263
4264	!-------------------------------------------------------------------------------
4265	! Sets the an array on the polar v point(s) to zero
4266	SUBROUTINE zero_pole ( field, &
4267	ids, ide, jds, jde, kds, kde, &
4268	ims, ime, jms, jme, kms, kme, &
4269	its, ite, jts, jte, kts, kte )
4270
4271
4272	IMPLICIT NONE
4273
4274	! Input data
4275
4276	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
4277	ims, ime, jms, jme, kms, kme, &
4278	its, ite, jts, jte, kts, kte
4279
4280	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: field
4281
4282	! Local data
4283
4284	INTEGER :: i, k
4285
4286	IF (jts == jds) THEN
4287	DO k = kts, kte
4288	DO i = its-1, ite+1
4289	field(i,k,jts) = 0.
4290	END DO
4291	END DO
4292	END IF
4293	IF (jte == jde) THEN
4294	DO k = kts, kte
4295	DO i = its-1, ite+1
4296	field(i,k,jte) = 0.
4297	END DO
4298	END DO
4299	END IF
4300
4301	END SUBROUTINE zero_pole
4302
4303	!-------------------------------------------------------------------------------
4304	! Sets the an array on the polar v point(s)
4305	SUBROUTINE pole_point_bc ( field, &
4306	ids, ide, jds, jde, kds, kde, &
4307	ims, ime, jms, jme, kms, kme, &
4308	its, ite, jts, jte, kts, kte )
4309
4310
4311	IMPLICIT NONE
4312
4313	! Input data
4314
4315	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
4316	ims, ime, jms, jme, kms, kme, &
4317	its, ite, jts, jte, kts, kte
4318
4319	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: field
4320
4321	! Local data
4322
4323	INTEGER :: i, k
4324
4325	IF (jts == jds) THEN
4326	DO k = kts, kte
4327	DO i = its, ite
4328	! field(i,k,jts) = 2*field(i,k,jts+1) - field(i,k,jts+2)
4329	field(i,k,jts) = field(i,k,jts+1)
4330	END DO
4331	END DO
4332	END IF
4333	IF (jte == jde) THEN
4334	DO k = kts, kte
4335	DO i = its, ite
4336	! field(i,k,jte) = 2*field(i,k,jte-1) - field(i,k,jte-2)
4337	field(i,k,jte) = field(i,k,jte-1)
4338	END DO
4339	END DO
4340	END IF
4341
4342	END SUBROUTINE pole_point_bc
4343
4344	!======================================================================
4345	! physics prep routines
4346	!======================================================================
4347
4348	SUBROUTINE phy_prep ( config_flags, & ! input
4349	mu, muu, muv, u, v, w, p, pb, alt, ph, & ! input
4350	phb, t, tsk, moist, n_moist, & ! input
4351	mu_3d, rho, th_phy, p_phy , pi_phy , & ! output
4352	u_phy, v_phy, w_phy, p8w, t_phy, t8w, & ! output
4353	z, z_at_w, dz8w, & ! output
4354	fzm, fzp, & ! params
4355	RTHRATEN, &
4356	RTHBLTEN, RUBLTEN, RVBLTEN, &
4357	RQVBLTEN, RQCBLTEN, RQIBLTEN, &
4358	RTHCUTEN, RQVCUTEN, RQCCUTEN, &
4359	RQRCUTEN, RQICUTEN, RQSCUTEN, &
4360	RTHFTEN, RQVFTEN, &
4361	RUNDGDTEN, RVNDGDTEN, RTHNDGDTEN, &
4362	RQVNDGDTEN, RMUNDGDTEN, &
4363	ids, ide, jds, jde, kds, kde, &
4364	ims, ime, jms, jme, kms, kme, &
4365	its, ite, jts, jte, kts, kte )
4366	!----------------------------------------------------------------------
4367	IMPLICIT NONE
4368	!----------------------------------------------------------------------
4369
4370	TYPE(grid_config_rec_type) , INTENT(IN ) :: config_flags
4371
4372	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
4373	ims, ime, jms, jme, kms, kme, &
4374	its, ite, jts, jte, kts, kte
4375	INTEGER , INTENT(IN ) :: n_moist
4376
4377	REAL, DIMENSION( ims:ime, kms:kme , jms:jme , n_moist ), INTENT(IN) :: moist
4378
4379
4380	REAL , DIMENSION( ims:ime, jms:jme ), INTENT(IN ) :: TSK, mu, muu, muv
4381
4382	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , &
4383	INTENT( OUT) :: u_phy, &
4384	v_phy, &
4385	w_phy, &
4386	pi_phy, &
4387	p_phy, &
4388	p8w, &
4389	t_phy, &
4390	th_phy, &
4391	t8w, &
4392	mu_3d, &
4393	rho, &
4394	z, &
4395	dz8w, &
4396	z_at_w
4397
4398	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , &
4399	INTENT(IN ) :: pb, &
4400	p, &
4401	u, &
4402	v, &
4403	w, &
4404	alt, &
4405	ph, &
4406	phb, &
4407	t
4408
4409
4410	REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fzm, &
4411	fzp
4412
4413	REAL, DIMENSION( ims:ime , kms:kme, jms:jme ), &
4414	INTENT(INOUT) :: RTHRATEN
4415
4416	REAL, DIMENSION( ims:ime , kms:kme, jms:jme ), &
4417	INTENT(INOUT) :: RTHCUTEN, &
4418	RQVCUTEN, &
4419	RQCCUTEN, &
4420	RQRCUTEN, &
4421	RQICUTEN, &
4422	RQSCUTEN
4423
4424	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) , &
4425	INTENT(INOUT) :: RUBLTEN, &
4426	RVBLTEN, &
4427	RTHBLTEN, &
4428	RQVBLTEN, &
4429	RQCBLTEN, &
4430	RQIBLTEN
4431
4432	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) , &
4433	INTENT(INOUT) :: RTHFTEN, &
4434	RQVFTEN
4435
4436	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) , &
4437	INTENT(INOUT) :: RUNDGDTEN, &
4438	RVNDGDTEN, &
4439	RTHNDGDTEN, &
4440	RQVNDGDTEN, &
4441	RMUNDGDTEN
4442
4443	INTEGER :: i_start, i_end, j_start, j_end, k_start, k_end, i_startu, j_startv
4444	INTEGER :: i, j, k
4445	REAL :: w1, w2, z0, z1, z2
4446
4447	!-----------------------------------------------------------------------
4448
4449	!<DESCRIPTION>
4450	!
4451	! phys_prep calculates a number of diagnostic quantities needed by
4452	! the physics routines. It also decouples the physics tendencies from
4453	! the column dry-air mass (the physics routines expect to see/update the
4454	! uncoupled tendencies).
4455	!
4456	!</DESCRIPTION>
4457
4458	! set up loop bounds for this grid's boundary conditions
4459
4460	i_start = its
4461	i_end = min( ite,ide-1 )
4462	j_start = jts
4463	j_end = min( jte,jde-1 )
4464
4465	k_start = kts
4466	k_end = min( kte, kde-1 )
4467
4468	! compute thermodynamics and velocities at pressure points
4469
4470	do j = j_start,j_end
4471	do k = k_start, k_end
4472	do i = i_start, i_end
4473
4474	th_phy(i,k,j) = t(i,k,j) + t0
4475	p_phy(i,k,j) = p(i,k,j) + pb(i,k,j)
4476	pi_phy(i,k,j) = (p_phy(i,k,j)/p1000mb)**rcp
4477	!! TAKE INTO ACCOUNT cp=f(T) on Venus
4478	IF (planet.eq. "venus" ) THEN
4479	t_phy(i,k,j)= (th_phy(i,k,j)*nu - nu(TT00*nu)log((p1000mb/p_phy(i,k,j))rcp))(1/nu)
4480	ELSE
4481	t_phy(i,k,j) = th_phy(i,k,j)*pi_phy(i,k,j)
4482	ENDIF
4483	rho(i,k,j) = 1./alt(i,k,j)*(1.+moist(i,k,j,P_QV))
4484	mu_3d(i,k,j) = mu(i,j)
4485	u_phy(i,k,j) = 0.5*(u(i,k,j)+u(i+1,k,j))
4486	v_phy(i,k,j) = 0.5*(v(i,k,j)+v(i,k,j+1))
4487
4488	enddo
4489	enddo
4490	enddo
4491
4492	! compute z at w points
4493
4494	do j = j_start,j_end
4495	do k = k_start, kte
4496	do i = i_start, i_end
4497	z_at_w(i,k,j) = (phb(i,k,j)+ph(i,k,j))/g
4498	enddo
4499	enddo
4500	enddo
4501
4502	do j = j_start,j_end
4503	do k = k_start, kte-1
4504	do i = i_start, i_end
4505	dz8w(i,k,j) = z_at_w(i,k+1,j)-z_at_w(i,k,j)
4506	enddo
4507	enddo
4508	enddo
4509
4510	do j = j_start,j_end
4511	do i = i_start, i_end
4512	dz8w(i,kte,j) = 0.
4513	enddo
4514	enddo
4515
4516	! compute z at p points (average of z at w points)
4517
4518	do j = j_start,j_end
4519	do k = k_start, k_end
4520	do i = i_start, i_end
4521	z(i,k,j) = 0.5*(z_at_w(i,k,j) + z_at_w(i,k+1,j) )
4522	!!!! MARS MARS ajout aymeric (ainsi que les arguments de cette routine)
4523	w_phy(i,k,j) = 0.5*(w(i,k,j) + w(i,k+1,j) )
4524	enddo
4525	enddo
4526	enddo
4527
4528	! interp t and p at w points
4529
4530	do j = j_start,j_end
4531	do k = 2, k_end
4532	do i = i_start, i_end
4533	p8w(i,k,j) = fzm(k)p_phy(i,k,j)+fzp(k)p_phy(i,k-1,j)
4534	t8w(i,k,j) = fzm(k)t_phy(i,k,j)+fzp(k)t_phy(i,k-1,j)
4535	enddo
4536	enddo
4537	enddo
4538
4539	! extrapolate p and t to surface and top.
4540	! we'll use an extrapolation in z for now
4541
4542	do j = j_start,j_end
4543	do i = i_start, i_end
4544
4545	! bottom
4546
4547	z0 = z_at_w(i,1,j)
4548	z1 = z(i,1,j)
4549	z2 = z(i,2,j)
4550	w1 = (z0 - z2)/(z1 - z2)
4551	w2 = 1. - w1
4552	p8w(i,1,j) = w1p_phy(i,1,j)+w2p_phy(i,2,j)
4553	t8w(i,1,j) = w1t_phy(i,1,j)+w2t_phy(i,2,j)
4554
4555	! top
4556
4557	z0 = z_at_w(i,kte,j)
4558	z1 = z(i,k_end,j)
4559	z2 = z(i,k_end-1,j)
4560	w1 = (z0 - z2)/(z1 - z2)
4561	w2 = 1. - w1
4562
4563	! p8w(i,kde,j) = w1p_phy(i,kde-1,j)+w2p_phy(i,kde-2,j)
4564	!!! bug fix extrapolate ln(p) so p is positive definite
4565	p8w(i,kde,j) = exp(w1log(p_phy(i,kde-1,j))+w2log(p_phy(i,kde-2,j)))
4566	t8w(i,kde,j) = w1t_phy(i,kde-1,j)+w2t_phy(i,kde-2,j)
4567
4568	enddo
4569	enddo
4570
4571	! decouple all physics tendencies
4572
4573	IF (config_flags%ra_lw_physics .gt. 0 .or. config_flags%ra_sw_physics .gt. 0) THEN
4574
4575	DO J=j_start,j_end
4576	DO K=k_start,k_end
4577	DO I=i_start,i_end
4578	RTHRATEN(I,K,J)=RTHRATEN(I,K,J)/mu(I,J)
4579	ENDDO
4580	ENDDO
4581	ENDDO
4582
4583	ENDIF
4584
4585	IF (config_flags%cu_physics .gt. 0) THEN
4586
4587	DO J=j_start,j_end
4588	DO I=i_start,i_end
4589	DO K=k_start,k_end
4590	RTHCUTEN(I,K,J)=RTHCUTEN(I,K,J)/mu(I,J)
4591	ENDDO
4592	ENDDO
4593	ENDDO
4594
4595	IF (P_QV .ge. PARAM_FIRST_SCALAR)THEN
4596	DO J=j_start,j_end
4597	DO I=i_start,i_end
4598	DO K=k_start,k_end
4599	RQVCUTEN(I,K,J)=RQVCUTEN(I,K,J)/mu(I,J)
4600	ENDDO
4601	ENDDO
4602	ENDDO
4603	ENDIF
4604
4605	IF (P_QC .ge. PARAM_FIRST_SCALAR)THEN
4606	DO J=j_start,j_end
4607	DO I=i_start,i_end
4608	DO K=k_start,k_end
4609	RQCCUTEN(I,K,J)=RQCCUTEN(I,K,J)/mu(I,J)
4610	ENDDO
4611	ENDDO
4612	ENDDO
4613	ENDIF
4614
4615	IF (P_QR .ge. PARAM_FIRST_SCALAR)THEN
4616	DO J=j_start,j_end
4617	DO I=i_start,i_end
4618	DO K=k_start,k_end
4619	RQRCUTEN(I,K,J)=RQRCUTEN(I,K,J)/mu(I,J)
4620	ENDDO
4621	ENDDO
4622	ENDDO
4623	ENDIF
4624
4625	IF (P_QI .ge. PARAM_FIRST_SCALAR)THEN
4626	DO J=j_start,j_end
4627	DO I=i_start,i_end
4628	DO K=k_start,k_end
4629	RQICUTEN(I,K,J)=RQICUTEN(I,K,J)/mu(I,J)
4630	ENDDO
4631	ENDDO
4632	ENDDO
4633	ENDIF
4634
4635	IF(P_QS .ge. PARAM_FIRST_SCALAR)THEN
4636	DO J=j_start,j_end
4637	DO I=i_start,i_end
4638	DO K=k_start,k_end
4639	RQSCUTEN(I,K,J)=RQSCUTEN(I,K,J)/mu(I,J)
4640	ENDDO
4641	ENDDO
4642	ENDDO
4643	ENDIF
4644
4645	ENDIF
4646
4647	!!MARS MARS
4648	! IF (config_flags%bl_pbl_physics .gt. 0) THEN
4649	IF ( (config_flags%bl_pbl_physics .gt. 0) &
4650	.OR. (config_flags%modif_wrf) ) THEN
4651	!****MARS
4652	DO J=j_start,j_end
4653	DO K=k_start,k_end
4654	DO I=i_start,i_end
4655	RUBLTEN(I,K,J) =RUBLTEN(I,K,J)/mu(I,J)
4656	RVBLTEN(I,K,J) =RVBLTEN(I,K,J)/mu(I,J)
4657	RTHBLTEN(I,K,J)=RTHBLTEN(I,K,J)/mu(I,J)
4658	ENDDO
4659	ENDDO
4660	ENDDO
4661
4662	IF (P_QV .ge. PARAM_FIRST_SCALAR) THEN
4663	DO J=j_start,j_end
4664	DO K=k_start,k_end
4665	DO I=i_start,i_end
4666	RQVBLTEN(I,K,J)=RQVBLTEN(I,K,J)/mu(I,J)
4667	ENDDO
4668	ENDDO
4669	ENDDO
4670	ENDIF
4671
4672	IF (P_QC .ge. PARAM_FIRST_SCALAR) THEN
4673	DO J=j_start,j_end
4674	DO K=k_start,k_end
4675	DO I=i_start,i_end
4676	RQCBLTEN(I,K,J)=RQCBLTEN(I,K,J)/mu(I,J)
4677	ENDDO
4678	ENDDO
4679	ENDDO
4680	ENDIF
4681
4682	IF (P_QI .ge. PARAM_FIRST_SCALAR) THEN
4683	DO J=j_start,j_end
4684	DO K=k_start,k_end
4685	DO I=i_start,i_end
4686	RQIBLTEN(I,K,J)=RQIBLTEN(I,K,J)/mu(I,J)
4687	ENDDO
4688	ENDDO
4689	ENDDO
4690	ENDIF
4691
4692	ENDIF
4693
4694	! decouple advective forcing required by Grell-Devenyi scheme
4695
4696	if(( config_flags%cu_physics == GDSCHEME ) .OR. &
4697	( config_flags%cu_physics == G3SCHEME )) then
4698
4699	DO J=j_start,j_end
4700	DO I=i_start,i_end
4701	DO K=k_start,k_end
4702	RTHFTEN(I,K,J)=RTHFTEN(I,K,J)/mu(I,J)
4703	ENDDO
4704	ENDDO
4705	ENDDO
4706
4707	IF (P_QV .ge. PARAM_FIRST_SCALAR)THEN
4708	DO J=j_start,j_end
4709	DO I=i_start,i_end
4710	DO K=k_start,k_end
4711	RQVFTEN(I,K,J)=RQVFTEN(I,K,J)/mu(I,J)
4712	ENDDO
4713	ENDDO
4714	ENDDO
4715	ENDIF
4716
4717	END IF
4718
4719	! fdda
4720	! note fdda u and v tendencies are staggered, also only interior points have muu/muv,
4721	! so only decouple those
4722
4723	IF (config_flags%grid_fdda .gt. 0) THEN
4724
4725	i_startu=MAX(its,ids+1)
4726	j_startv=MAX(jts,jds+1)
4727
4728	DO J=j_start,j_end
4729	DO K=k_start,k_end
4730	DO I=i_startu,i_end
4731	RUNDGDTEN(I,K,J) =RUNDGDTEN(I,K,J)/muu(I,J)
4732	ENDDO
4733	ENDDO
4734	ENDDO
4735	DO J=j_startv,j_end
4736	DO K=k_start,k_end
4737	DO I=i_start,i_end
4738	RVNDGDTEN(I,K,J) =RVNDGDTEN(I,K,J)/muv(I,J)
4739	ENDDO
4740	ENDDO
4741	ENDDO
4742	DO J=j_start,j_end
4743	DO K=k_start,k_end
4744	DO I=i_start,i_end
4745	RTHNDGDTEN(I,K,J)=RTHNDGDTEN(I,K,J)/mu(I,J)
4746	! RMUNDGDTEN(I,J) - no coupling
4747	ENDDO
4748	ENDDO
4749	ENDDO
4750	IF (P_QV .ge. PARAM_FIRST_SCALAR) THEN
4751	DO J=j_start,j_end
4752	DO K=k_start,k_end
4753	DO I=i_start,i_end
4754	RQVNDGDTEN(I,K,J)=RQVNDGDTEN(I,K,J)/mu(I,J)
4755	ENDDO
4756	ENDDO
4757	ENDDO
4758	ENDIF
4759
4760	ENDIF
4761
4762	END SUBROUTINE phy_prep
4763
4764	!------------------------------------------------------------
4765
4766	SUBROUTINE moist_physics_prep_em( t_new, t_old, t0, rho, al, alb, &
4767	p, p8w, p0, pb, ph, phb, &
4768	th_phy, pii, pf, &
4769	z, z_at_w, dz8w, &
4770	dt,h_diabatic, &
4771	config_flags,fzm, fzp, &
4772	ids,ide, jds,jde, kds,kde, &
4773	ims,ime, jms,jme, kms,kme, &
4774	its,ite, jts,jte, kts,kte )
4775
4776	IMPLICIT NONE
4777
4778	! Here we construct full fields
4779	! needed by the microphysics
4780
4781	TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
4782
4783	INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde
4784	INTEGER, INTENT(IN ) :: ims,ime, jms,jme, kms,kme
4785	INTEGER, INTENT(IN ) :: its,ite, jts,jte, kts,kte
4786
4787	REAL, INTENT(IN ) :: dt
4788
4789	REAL, DIMENSION( ims:ime , kms:kme, jms:jme ), &
4790	INTENT(IN ) :: al, &
4791	alb, &
4792	p, &
4793	pb, &
4794	ph, &
4795	phb
4796
4797
4798	REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fzm, &
4799	fzp
4800
4801	REAL, DIMENSION( ims:ime , kms:kme, jms:jme ), &
4802	INTENT( OUT) :: rho, &
4803	th_phy, &
4804	pii, &
4805	pf, &
4806	z, &
4807	z_at_w, &
4808	dz8w, &
4809	p8w
4810
4811	REAL, DIMENSION( ims:ime , kms:kme, jms:jme ), &
4812	INTENT(INOUT) :: h_diabatic
4813
4814	REAL, DIMENSION( ims:ime , kms:kme, jms:jme ), &
4815	INTENT(INOUT) :: t_new, &
4816	t_old
4817
4818	REAL, INTENT(IN ) :: t0, p0
4819	REAL :: z0,z1,z2,w1,w2
4820
4821	INTEGER :: i_start, i_end, j_start, j_end, k_start, k_end
4822	INTEGER :: i, j, k
4823
4824	!--------------------------------------------------------------------
4825
4826	!<DESCRIPTION>
4827	!
4828	! moist_phys_prep_em calculates a number of diagnostic quantities needed by
4829	! the microphysics routines.
4830	!
4831	!</DESCRIPTION>
4832
4833	! set up loop bounds for this grid's boundary conditions
4834
4835	i_start = its
4836	i_end = min( ite,ide-1 )
4837	j_start = jts
4838	j_end = min( jte,jde-1 )
4839
4840	k_start = kts
4841	k_end = min( kte, kde-1 )
4842
4843	DO j = j_start, j_end
4844	DO k = k_start, kte
4845	DO i = i_start, i_end
4846	z_at_w(i,k,j) = (ph(i,k,j)+phb(i,k,j))/g
4847	ENDDO
4848	ENDDO
4849	ENDDO
4850
4851	do j = j_start,j_end
4852	do k = k_start, kte-1
4853	do i = i_start, i_end
4854	dz8w(i,k,j) = z_at_w(i,k+1,j)-z_at_w(i,k,j)
4855	enddo
4856	enddo
4857	enddo
4858
4859	do j = j_start,j_end
4860	do i = i_start, i_end
4861	dz8w(i,kte,j) = 0.
4862	enddo
4863	enddo
4864
4865
4866	! compute full pii, rho, and z at the new time-level
4867	! (needed for physics).
4868	! convert perturbation theta to full theta (th_phy)
4869	! use h_diabatic to temporarily save pre-microphysics full theta
4870
4871	DO j = j_start, j_end
4872	DO k = k_start, k_end
4873	DO i = i_start, i_end
4874
4875	#ifdef REVERT
4876	t_new(i,k,j) = t_new(i,k,j)-h_diabatic(i,k,j)*dt
4877	#endif
4878	th_phy(i,k,j) = t_new(i,k,j) + t0
4879	h_diabatic(i,k,j) = th_phy(i,k,j)
4880	rho(i,k,j) = 1./(al(i,k,j)+alb(i,k,j))
4881	pii(i,k,j) = ((p(i,k,j)+pb(i,k,j))/p0)**rcp
4882	z(i,k,j) = 0.5*(z_at_w(i,k,j) +z_at_w(i,k+1,j) )
4883	pf(i,k,j) = p(i,k,j)+pb(i,k,j)
4884
4885	ENDDO
4886	ENDDO
4887	ENDDO
4888
4889	! interp t and p at w points
4890
4891	do j = j_start,j_end
4892	do k = 2, k_end
4893	do i = i_start, i_end
4894	p8w(i,k,j) = fzm(k)pf(i,k,j)+fzp(k)pf(i,k-1,j)
4895	enddo
4896	enddo
4897	enddo
4898
4899	! extrapolate p and t to surface and top.
4900	! we'll use an extrapolation in z for now
4901
4902	do j = j_start,j_end
4903	do i = i_start, i_end
4904
4905	! bottom
4906
4907	z0 = z_at_w(i,1,j)
4908	z1 = z(i,1,j)
4909	z2 = z(i,2,j)
4910	w1 = (z0 - z2)/(z1 - z2)
4911	w2 = 1. - w1
4912	p8w(i,1,j) = w1pf(i,1,j)+w2pf(i,2,j)
4913
4914	! top
4915
4916	z0 = z_at_w(i,kte,j)
4917	z1 = z(i,k_end,j)
4918	z2 = z(i,k_end-1,j)
4919	w1 = (z0 - z2)/(z1 - z2)
4920	w2 = 1. - w1
4921	! p8w(i,kde,j) = w1pf(i,kde-1,j)+w2pf(i,kde-2,j)
4922	p8w(i,kde,j) = exp(w1log(pf(i,kde-1,j))+w2log(pf(i,kde-2,j)))
4923
4924	enddo
4925	enddo
4926
4927	END SUBROUTINE moist_physics_prep_em
4928
4929	!------------------------------------------------------------------------------
4930
4931	SUBROUTINE moist_physics_finish_em( t_new, t_old, t0, mut, &
4932	th_phy, h_diabatic, dt, &
4933	config_flags, &
4934	ids,ide, jds,jde, kds,kde, &
4935	ims,ime, jms,jme, kms,kme, &
4936	its,ite, jts,jte, kts,kte )
4937
4938	IMPLICIT NONE
4939
4940	! Here we construct full fields
4941	! needed by the microphysics
4942
4943	TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
4944
4945	INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde
4946	INTEGER, INTENT(IN ) :: ims,ime, jms,jme, kms,kme
4947	INTEGER, INTENT(IN ) :: its,ite, jts,jte, kts,kte
4948
4949	REAL, DIMENSION( ims:ime , kms:kme, jms:jme ), &
4950	INTENT(INOUT) :: t_new, &
4951	t_old, &
4952	th_phy, &
4953	h_diabatic
4954
4955	REAL, DIMENSION( ims:ime , jms:jme ), INTENT(INOUT) :: mut
4956
4957
4958	REAL, INTENT(IN ) :: t0, dt
4959
4960	INTEGER :: i_start, i_end, j_start, j_end, k_start, k_end
4961	INTEGER :: i, j, k
4962
4963	!--------------------------------------------------------------------
4964
4965	!<DESCRIPTION>
4966	!
4967	! moist_phys_finish_em resets theta to its perturbation value and
4968	! computes and stores the microphysics diabatic heating term.
4969	!
4970	!</DESCRIPTION>
4971
4972	! set up loop bounds for this grid's boundary conditions
4973
4974
4975	i_start = its
4976	i_end = min( ite,ide-1 )
4977	j_start = jts
4978	j_end = min( jte,jde-1 )
4979
4980	k_start = kts
4981	k_end = min( kte, kde-1 )
4982
4983	! add microphysics theta diff to perturbation theta, set h_diabatic
4984
4985	IF ( config_flags%no_mp_heating .eq. 0 ) THEN
4986	DO j = j_start, j_end
4987	DO k = k_start, k_end
4988	DO i = i_start, i_end
4989	t_new(i,k,j) = t_new(i,k,j) + (th_phy(i,k,j)-h_diabatic(i,k,j))
4990	h_diabatic(i,k,j) = (th_phy(i,k,j)-h_diabatic(i,k,j))/dt
4991	ENDDO
4992	ENDDO
4993	ENDDO
4994
4995	ELSE
4996
4997	DO j = j_start, j_end
4998	DO k = k_start, k_end
4999	DO i = i_start, i_end
5000	! t_new(i,k,j) = t_new(i,k,j)
5001	h_diabatic(i,k,j) = 0.
5002	ENDDO
5003	ENDDO
5004	ENDDO
5005	ENDIF
5006
5007	END SUBROUTINE moist_physics_finish_em
5008
5009	!----------------------------------------------------------------
5010
5011
5012	SUBROUTINE init_module_big_step
5013	END SUBROUTINE init_module_big_step
5014
5015	SUBROUTINE set_tend ( field, field_adv_tend, msf, &
5016	ids, ide, jds, jde, kds, kde, &
5017	ims, ime, jms, jme, kms, kme, &
5018	its, ite, jts, jte, kts, kte )
5019
5020	IMPLICIT NONE
5021
5022	! Input data
5023
5024	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
5025	ims, ime, jms, jme, kms, kme, &
5026	its, ite, jts, jte, kts, kte
5027
5028	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT) :: field
5029
5030	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN) :: field_adv_tend
5031
5032	REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN) :: msf
5033
5034	! Local data
5035
5036	INTEGER :: i, j, k, itf, jtf, ktf
5037
5038	!<DESCRIPTION>
5039	!
5040	! set_tend copies the advective tendency array into the tendency array.
5041	!
5042	!</DESCRIPTION>
5043
5044	jtf = MIN(jte,jde-1)
5045	ktf = MIN(kte,kde-1)
5046	itf = MIN(ite,ide-1)
5047	DO j = jts, jtf
5048	DO k = kts, ktf
5049	DO i = its, itf
5050	field(i,k,j) = field_adv_tend(i,k,j)*msf(i,j)
5051	ENDDO
5052	ENDDO
5053	ENDDO
5054
5055	END SUBROUTINE set_tend
5056
5057	!------------------------------------------------------------------------------
5058
5059	SUBROUTINE rk_rayleigh_damp( ru_tendf, rv_tendf, &
5060	rw_tendf, t_tendf, &
5061	u, v, w, t, t_init, &
5062	mut, muu, muv, ph, phb, &
5063	u_base, v_base, t_base, z_base, &
5064	dampcoef, zdamp, &
5065	ids, ide, jds, jde, kds, kde, &
5066	ims, ime, jms, jme, kms, kme, &
5067	its, ite, jts, jte, kts, kte )
5068
5069	! History: Apr 2005 Modifications by George Bryan, NCAR:
5070	! - Generalized the code in a way that allows for
5071	! simulations with steep terrain.
5072	!
5073	! Jul 2004 Modifications by George Bryan, NCAR:
5074	! - Modified the code to use u_base, v_base, and t_base
5075	! arrays for the background state. Removed the hard-wired
5076	! base-state values.
5077	! - Modified the code to use dampcoef, zdamp, and damp_opt,
5078	! i.e., the upper-level damper variables in namelist.input.
5079	! Removed the hard-wired variables in the older version.
5080	! This damper is used when damp_opt = 2.
5081	! - Modified the code to account for the movement of the
5082	! model surfaces with time. The code now obtains a base-
5083	! state value by interpolation using the "_base" arrays.
5084
5085	! Nov 2003 Bug fix by Jason Knievel, NCAR
5086
5087	! Aug 2003 Meridional dimension, some comments, and
5088	! changes in layout of the code added by
5089	! Jason Knievel, NCAR
5090
5091	! Jul 2003 Original code by Bill Skamarock, NCAR
5092
5093	! Purpose: This routine applies Rayleigh damping to a layer at top
5094	! of the model domain.
5095
5096	!-----------------------------------------------------------------------
5097	! Begin declarations.
5098
5099	IMPLICIT NONE
5100
5101	INTEGER, INTENT( IN ) &
5102	:: ids, ide, jds, jde, kds, kde, &
5103	ims, ime, jms, jme, kms, kme, &
5104	its, ite, jts, jte, kts, kte
5105
5106	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( INOUT ) &
5107	:: ru_tendf, rv_tendf, rw_tendf, t_tendf
5108
5109	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( IN ) &
5110	:: u, v, w, t, t_init, ph, phb
5111
5112	REAL, DIMENSION( ims:ime, jms:jme ), INTENT( IN ) &
5113	:: mut, muu, muv
5114
5115	REAL, DIMENSION( kms:kme ) , INTENT(IN ) &
5116	:: u_base, v_base, t_base, z_base
5117
5118	REAL, INTENT(IN ) &
5119	:: dampcoef, zdamp
5120
5121	! Local variables.
5122
5123	INTEGER &
5124	:: i_start, i_end, j_start, j_end, k_start, k_end, i, j, k, ktf, k1, k2
5125
5126	REAL &
5127	:: pii, dcoef, z, ztop
5128
5129	REAL :: wkp1, wk, wkm1
5130
5131	REAL, DIMENSION( kms:kme ) :: z00, u00, v00, t00
5132
5133	! End declarations.
5134	!-----------------------------------------------------------------------
5135
5136	pii = 2.0 * asin(1.0)
5137
5138	ktf = MIN( kte, kde-1 )
5139
5140	!-----------------------------------------------------------------------
5141	! Adjust u to base state.
5142
5143	DO j = jts, MIN( jte, jde-1 )
5144	DO i = its, MIN( ite, ide )
5145
5146	! Get height at top of model
5147	ztop = 0.5*( phb(i ,kde,j)+phb(i-1,kde,j) &
5148	+ph(i ,kde,j)+ph(i-1,kde,j) )/g
5149
5150	! Find bottom of damping layer
5151	k1 = ktf
5152	z = ztop
5153	DO WHILE( z >= (ztop-zdamp) )
5154	z = 0.25*( phb(i ,k1,j)+phb(i ,k1+1,j) &
5155	+phb(i-1,k1,j)+phb(i-1,k1+1,j) &
5156	+ph(i ,k1,j)+ph(i ,k1+1,j) &
5157	+ph(i-1,k1,j)+ph(i-1,k1+1,j))/g
5158	z00(k1) = z
5159	k1 = k1 - 1
5160	ENDDO
5161	k1 = k1 + 2
5162
5163	! Get reference state at model levels
5164	DO k = k1, ktf
5165	k2 = ktf
5166	DO WHILE( z_base(k2) .gt. z00(k) )
5167	k2 = k2 - 1
5168	ENDDO
5169	if(k2+1.gt.ktf)then
5170	u00(k) = u_base(k2) + ( u_base(k2) - u_base(k2-1) ) &
5171	* ( z00(k) - z_base(k2) ) &
5172	/ ( z_base(k2) - z_base(k2-1) )
5173	else
5174	u00(k) = u_base(k2) + ( u_base(k2+1) - u_base(k2) ) &
5175	* ( z00(k) - z_base(k2) ) &
5176	/ ( z_base(k2+1) - z_base(k2) )
5177	endif
5178	ENDDO
5179
5180	! Apply the Rayleigh damper
5181	DO k = k1, ktf
5182	dcoef = 1.0 - MIN( 1.0, ( ztop - z00(k) ) / zdamp )
5183	dcoef = (SIN( 0.5 * pii * dcoef ) )**2
5184	ru_tendf(i,k,j) = ru_tendf(i,k,j) - &
5185	muu(i,j) * ( dcoef * dampcoef ) * &
5186	( u(i,k,j) - u00(k) )
5187	END DO
5188
5189	END DO
5190	END DO
5191
5192	! End adjustment of u.
5193	!-----------------------------------------------------------------------
5194
5195	!-----------------------------------------------------------------------
5196	! Adjust v to base state.
5197
5198	DO j = jts, MIN( jte, jde )
5199	DO i = its, MIN( ite, ide-1 )
5200
5201	! Get height at top of model
5202	ztop = 0.5*( phb(i,kde,j )+phb(i,kde,j-1) &
5203	+ph(i,kde,j )+ph(i,kde,j-1) )/g
5204
5205	! Find bottom of damping layer
5206	k1 = ktf
5207	z = ztop
5208	DO WHILE( z >= (ztop-zdamp) )
5209	z = 0.25*( phb(i,k1,j )+phb(i,k1+1,j ) &
5210	+phb(i,k1,j-1)+phb(i,k1+1,j-1) &
5211	+ph(i,k1,j )+ph(i,k1+1,j ) &
5212	+ph(i,k1,j-1)+ph(i,k1+1,j-1))/g
5213	z00(k1) = z
5214	k1 = k1 - 1
5215	ENDDO
5216	k1 = k1 + 2
5217
5218	! Get reference state at model levels
5219	DO k = k1, ktf
5220	k2 = ktf
5221	DO WHILE( z_base(k2) .gt. z00(k) )
5222	k2 = k2 - 1
5223	ENDDO
5224	if(k2+1.gt.ktf)then
5225	v00(k) = v_base(k2) + ( v_base(k2) - v_base(k2-1) ) &
5226	* ( z00(k) - z_base(k2) ) &
5227	/ ( z_base(k2) - z_base(k2-1) )
5228	else
5229	v00(k) = v_base(k2) + ( v_base(k2+1) - v_base(k2) ) &
5230	* ( z00(k) - z_base(k2) ) &
5231	/ ( z_base(k2+1) - z_base(k2) )
5232	endif
5233	ENDDO
5234
5235	! Apply the Rayleigh damper
5236	DO k = k1, ktf
5237	dcoef = 1.0 - MIN( 1.0, ( ztop - z00(k) ) / zdamp )
5238	dcoef = (SIN( 0.5 * pii * dcoef ) )**2
5239	rv_tendf(i,k,j) = rv_tendf(i,k,j) - &
5240	muv(i,j) * ( dcoef * dampcoef ) * &
5241	( v(i,k,j) - v00(k) )
5242	END DO
5243
5244	END DO
5245	END DO
5246
5247	! End adjustment of v.
5248	!-----------------------------------------------------------------------
5249
5250	!-----------------------------------------------------------------------
5251	! Adjust w to base state.
5252
5253	DO j = jts, MIN( jte, jde-1 )
5254	DO i = its, MIN( ite, ide-1 )
5255	ztop = ( phb(i,kde,j) + ph(i,kde,j) ) / g
5256	DO k = kts, MIN( kte, kde )
5257	z = ( phb(i,k,j) + ph(i,k,j) ) / g
5258	IF ( z >= (ztop-zdamp) ) THEN
5259	dcoef = 1.0 - MIN( 1.0, ( ztop - z ) / zdamp )
5260	dcoef = ( SIN( 0.5 * pii * dcoef ) )**2
5261	rw_tendf(i,k,j) = rw_tendf(i,k,j) - &
5262	mut(i,j) * ( dcoef * dampcoef ) * w(i,k,j)
5263	END IF
5264	END DO
5265	END DO
5266	END DO
5267
5268	! End adjustment of w.
5269	!-----------------------------------------------------------------------
5270
5271	!-----------------------------------------------------------------------
5272	! Adjust potential temperature to base state.
5273
5274	DO j = jts, MIN( jte, jde-1 )
5275	DO i = its, MIN( ite, ide-1 )
5276
5277	! Get height at top of model
5278	ztop = ( phb(i,kde,j) + ph(i,kde,j) ) / g
5279
5280	! Find bottom of damping layer
5281	k1 = ktf
5282	z = ztop
5283	DO WHILE( z >= (ztop-zdamp) )
5284	z = 0.5 * ( phb(i,k1,j) + phb(i,k1+1,j) + &
5285	ph(i,k1,j) + ph(i,k1+1,j) ) / g
5286	z00(k1) = z
5287	k1 = k1 - 1
5288	ENDDO
5289	k1 = k1 + 2
5290
5291	! Get reference state at model levels
5292	DO k = k1, ktf
5293	k2 = ktf
5294	DO WHILE( z_base(k2) .gt. z00(k) )
5295	k2 = k2 - 1
5296	ENDDO
5297	if(k2+1.gt.ktf)then
5298	t00(k) = t_base(k2) + ( t_base(k2) - t_base(k2-1) ) &
5299	* ( z00(k) - z_base(k2) ) &
5300	/ ( z_base(k2) - z_base(k2-1) )
5301	else
5302	t00(k) = t_base(k2) + ( t_base(k2+1) - t_base(k2) ) &
5303	* ( z00(k) - z_base(k2) ) &
5304	/ ( z_base(k2+1) - z_base(k2) )
5305	endif
5306	ENDDO
5307
5308	! Apply the Rayleigh damper
5309	DO k = k1, ktf
5310	dcoef = 1.0 - MIN( 1.0, ( ztop - z00(k) ) / zdamp )
5311	dcoef = (SIN( 0.5 * pii * dcoef ) )**2
5312	t_tendf(i,k,j) = t_tendf(i,k,j) - &
5313	mut(i,j) * ( dcoef * dampcoef ) * &
5314	( t(i,k,j) - t00(k) )
5315	END DO
5316
5317	END DO
5318	END DO
5319
5320	! End adjustment of potential temperature.
5321	!-----------------------------------------------------------------------
5322
5323	END SUBROUTINE rk_rayleigh_damp
5324
5325	!==============================================================================
5326	!==============================================================================
5327
5328	SUBROUTINE sixth_order_diffusion( name, field, tendency, mu, dt, &
5329	config_flags, &
5330	diff_6th_opt, diff_6th_factor, &
5331	ids, ide, jds, jde, kds, kde, &
5332	ims, ime, jms, jme, kms, kme, &
5333	its, ite, jts, jte, kts, kte )
5334
5335	! History: 14 Nov 2006 Name of variable changed by Jason Knievel
5336	! 07 Jun 2006 Revised and generalized by Jason Knievel
5337	! 25 Apr 2005 Original code by Jason Knievel, NCAR
5338
5339	! Purpose: Apply 6th-order, monotonic (flux-limited), numerical
5340	! diffusion to 3-d velocity and to scalars.
5341
5342	! References: Ming Xue (MWR Aug 2000)
5343	! Durran ("Numerical Methods for Wave Equations..." 1999)
5344	! George Bryan (personal communication)
5345
5346	!------------------------------------------------------------------------------
5347	! Begin: Declarations.
5348
5349	IMPLICIT NONE
5350
5351	INTEGER, INTENT(IN) &
5352	:: ids, ide, jds, jde, kds, kde, &
5353	ims, ime, jms, jme, kms, kme, &
5354	its, ite, jts, jte, kts, kte
5355
5356	TYPE(grid_config_rec_type), INTENT(IN) &
5357	:: config_flags
5358
5359	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) &
5360	:: tendency
5361
5362	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN) &
5363	:: field
5364
5365	REAL, DIMENSION( ims:ime , jms:jme ), INTENT(IN) &
5366	:: mu
5367
5368	REAL, INTENT(IN) &
5369	:: dt
5370
5371	REAL, INTENT(IN) &
5372	:: diff_6th_factor
5373
5374	INTEGER, INTENT(IN) &
5375	:: diff_6th_opt
5376
5377	CHARACTER(LEN=1) , INTENT(IN) &
5378	:: name
5379
5380	INTEGER &
5381	:: i, j, k, &
5382	i_start, i_end, &
5383	j_start, j_end, &
5384	k_start, k_end, &
5385	ktf
5386
5387	REAL &
5388	:: dflux_x_p0, dflux_y_p0, &
5389	dflux_x_p1, dflux_y_p1, &
5390	tendency_x, tendency_y, &
5391	mu_avg_p0, mu_avg_p1, &
5392	diff_6th_coef
5393
5394	LOGICAL &
5395	:: specified
5396
5397	! End: Declarations.
5398	!------------------------------------------------------------------------------
5399
5400	!------------------------------------------------------------------------------
5401	! Begin: Translate the diffusion factor into a diffusion coefficient. See
5402	! Durran's text, section 2.4.3, then adjust for sixth-order diffusion (not
5403	! fourth) and for diffusion in two dimensions (not one). For reference, a
5404	! factor of 1.0 would mean complete diffusion of a 2dx wave in one time step,
5405	! although application of the flux limiter reduces somewhat the effects of
5406	! diffusion for a given coefficient.
5407
5408	diff_6th_coef = diff_6th_factor * 0.015625 / ( 2.0 * dt )
5409
5410	! End: Translate diffusion factor.
5411	!------------------------------------------------------------------------------
5412
5413	!------------------------------------------------------------------------------
5414	! Begin: Assign limits of spatial loops depending on variable to be diffused.
5415	! The halo regions are already filled with values by the time this subroutine
5416	! is called, which allows the stencil to extend beyond the domains' edges.
5417
5418	ktf = MIN( kte, kde-1 )
5419
5420	IF ( name .EQ. 'u' ) THEN
5421
5422	i_start = its
5423	i_end = ite
5424	j_start = jts
5425	j_end = MIN(jde-1,jte)
5426	k_start = kts
5427	k_end = ktf
5428
5429	ELSE IF ( name .EQ. 'v' ) THEN
5430
5431	i_start = its
5432	i_end = MIN(ide-1,ite)
5433	j_start = jts
5434	j_end = jte
5435	k_start = kts
5436	k_end = ktf
5437
5438	ELSE IF ( name .EQ. 'w' ) THEN
5439
5440	i_start = its
5441	i_end = MIN(ide-1,ite)
5442	j_start = jts
5443	j_end = MIN(jde-1,jte)
5444	k_start = kts+1
5445	k_end = ktf
5446
5447	ELSE
5448
5449	i_start = its
5450	i_end = MIN(ide-1,ite)
5451	j_start = jts
5452	j_end = MIN(jde-1,jte)
5453	k_start = kts
5454	k_end = ktf
5455
5456	ENDIF
5457
5458	! End: Assignment of limits of spatial loops.
5459	!------------------------------------------------------------------------------
5460
5461	!------------------------------------------------------------------------------
5462	! Begin: Loop across spatial dimensions.
5463
5464	DO j = j_start, j_end
5465	DO k = k_start, k_end
5466	DO i = i_start, i_end
5467
5468	!------------------------------------------------------------------------------
5469	! Begin: Diffusion in x (i index).
5470
5471	! Calculate the diffusive flux in x direction (from Xue's eq. 3).
5472
5473	dflux_x_p0 = ( 10.0 * ( field(i, k,j) - field(i-1,k,j) ) &
5474	- 5.0 * ( field(i+1,k,j) - field(i-2,k,j) ) &
5475	+ ( field(i+2,k,j) - field(i-3,k,j) ) )
5476
5477	dflux_x_p1 = ( 10.0 * ( field(i+1,k,j) - field(i ,k,j) ) &
5478	- 5.0 * ( field(i+2,k,j) - field(i-1,k,j) ) &
5479	+ ( field(i+3,k,j) - field(i-2,k,j) ) )
5480
5481	! If requested in the namelist (diff_6th_opt=2), prohibit up-gradient diffusion
5482	! (variation on Xue's eq. 10).
5483
5484	IF ( diff_6th_opt .EQ. 2 ) THEN
5485
5486	IF ( dflux_x_p0 * ( field(i ,k,j)-field(i-1,k,j) ) .LE. 0.0 ) THEN
5487	dflux_x_p0 = 0.0
5488	END IF
5489
5490	IF ( dflux_x_p1 * ( field(i+1,k,j)-field(i ,k,j) ) .LE. 0.0 ) THEN
5491	dflux_x_p1 = 0.0
5492	END IF
5493
5494	END IF
5495
5496	! Apply 6th-order diffusion in x direction.
5497
5498	IF ( name .EQ. 'u' ) THEN
5499	mu_avg_p0 = mu(i-1,j)
5500	mu_avg_p1 = mu(i ,j)
5501	ELSE IF ( name .EQ. 'v' ) THEN
5502	mu_avg_p0 = 0.25 * ( &
5503	mu(i-1,j-1) + &
5504	mu(i ,j-1) + &
5505	mu(i-1,j ) + &
5506	mu(i ,j ) )
5507	mu_avg_p1 = 0.25 * ( &
5508	mu(i ,j-1) + &
5509	mu(i+1,j-1) + &
5510	mu(i ,j ) + &
5511	mu(i+1,j ) )
5512	ELSE
5513	mu_avg_p0 = 0.5 * ( &
5514	mu(i-1,j) + &
5515	mu(i ,j) )
5516	mu_avg_p1 = 0.5 * ( &
5517	mu(i ,j) + &
5518	mu(i+1,j) )
5519	END IF
5520
5521	tendency_x = diff_6th_coef * &
5522	( ( mu_avg_p1 * dflux_x_p1 ) - ( mu_avg_p0 * dflux_x_p0 ) )
5523
5524	! End: Diffusion in x.
5525	!------------------------------------------------------------------------------
5526
5527	!------------------------------------------------------------------------------
5528	! Begin: Diffusion in y (j index).
5529
5530	! Calculate the diffusive flux in y direction (from Xue's eq. 3).
5531
5532	dflux_y_p0 = ( 10.0 * ( field(i,k,j ) - field(i,k,j-1) ) &
5533	- 5.0 * ( field(i,k,j+1) - field(i,k,j-2) ) &
5534	+ ( field(i,k,j+2) - field(i,k,j-3) ) )
5535
5536	dflux_y_p1 = ( 10.0 * ( field(i,k,j+1) - field(i,k,j ) ) &
5537	- 5.0 * ( field(i,k,j+2) - field(i,k,j-1) ) &
5538	+ ( field(i,k,j+3) - field(i,k,j-2) ) )
5539
5540	! If requested in the namelist (diff_6th_opt=2), prohibit up-gradient diffusion
5541	! (variation on Xue's eq. 10).
5542
5543	IF ( diff_6th_opt .EQ. 2 ) THEN
5544
5545	IF ( dflux_y_p0 * ( field(i,k,j )-field(i,k,j-1) ) .LE. 0.0 ) THEN
5546	dflux_y_p0 = 0.0
5547	END IF
5548
5549	IF ( dflux_y_p1 * ( field(i,k,j+1)-field(i,k,j ) ) .LE. 0.0 ) THEN
5550	dflux_y_p1 = 0.0
5551	END IF
5552
5553	END IF
5554
5555	! Apply 6th-order diffusion in y direction.
5556
5557	IF ( name .EQ. 'u' ) THEN
5558	mu_avg_p0 = 0.25 * ( &
5559	mu(i-1,j-1) + &
5560	mu(i ,j-1) + &
5561	mu(i-1,j ) + &
5562	mu(i ,j ) )
5563	mu_avg_p1 = 0.25 * ( &
5564	mu(i-1,j ) + &
5565	mu(i ,j ) + &
5566	mu(i-1,j+1) + &
5567	mu(i ,j+1) )
5568	ELSE IF ( name .EQ. 'v' ) THEN
5569	mu_avg_p0 = mu(i,j-1)
5570	mu_avg_p1 = mu(i,j )
5571	ELSE
5572	mu_avg_p0 = 0.5 * ( &
5573	mu(i,j-1) + &
5574	mu(i,j ) )
5575	mu_avg_p1 = 0.5 * ( &
5576	mu(i,j ) + &
5577	mu(i,j+1) )
5578	END IF
5579
5580	tendency_y = diff_6th_coef * &
5581	( ( mu_avg_p1 * dflux_y_p1 ) - ( mu_avg_p0 * dflux_y_p0 ) )
5582
5583	! End: Diffusion in y.
5584	!------------------------------------------------------------------------------
5585
5586	!------------------------------------------------------------------------------
5587	! Begin: Combine diffusion in x and y.
5588
5589	tendency(i,k,j) = tendency(i,k,j) + tendency_x + tendency_y
5590
5591	! End: Combine diffusion in x and y.
5592	!------------------------------------------------------------------------------
5593
5594	ENDDO
5595	ENDDO
5596	ENDDO
5597
5598	! End: Loop across spatial dimensions.
5599	!------------------------------------------------------------------------------
5600
5601	END SUBROUTINE sixth_order_diffusion
5602
5603	!==============================================================================
5604	!==============================================================================
5605
5606	END MODULE module_big_step_utilities_em

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

Download in other formats:

Original Format