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