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