PageRenderTime 52ms CodeModel.GetById 16ms RepoModel.GetById 0ms app.codeStats 0ms

/wrfv2_fire/dyn_em/module_small_step_em.F

https://github.com/jbeezley/wrf-fire
FORTRAN Legacy | 1772 lines | 1148 code | 325 blank | 299 comment | 13 complexity | d599c6f45b76074d79b180a438c884fd MD5 | raw file
Possible License(s): AGPL-1.0

Large files files are truncated, but you can click here to view the full file

    

  1. !WRF:MODEL_LAYER:DYNAMICS
    2. 

  2. !
    3. 

  3. 

  4. !  SMALL_STEP code for the geometric height coordinate model
    5. 

  5. !
    6. 

  6. !---------------------------------------------------------------------------
    7. 

  7. 

  8. MODULE module_small_step_em
    9. 

  9. 

  10.    USE module_configure
    11. 

  11.    USE module_model_constants
    12. 

  12. 

  13. CONTAINS
    14. 

  14. 

  15. !----------------------------------------------------------------------
    16. 

  16. 

  17. SUBROUTINE small_step_prep( u_1, u_2, v_1, v_2, w_1, w_2, &
    18. 

  18.                             t_1, t_2, ph_1, ph_2,         &
    19. 

  19.                             mub, mu_1, mu_2,              &
    20. 

  20.                             muu, muus, muv, muvs,         &
    21. 

  21.                             mut, muts, mudf,              &
    22. 

  22.                             u_save, v_save, w_save,       &
    23. 

  23.                             t_save, ph_save, mu_save,     &
    24. 

  24.                             ww, ww_save,                  &
    25. 

  25.                             dnw, c2a, pb, p, alt,         &
    26. 

  26.                             msfux, msfuy, msfvx,          &
    27. 

  27.                             msfvx_inv,                    &
    28. 

  28.                             msfvy, msftx, msfty,          &
    29. 

  29.                             rdx, rdy,                     &
    30. 

  30.                             rk_step,                      &
    31. 

  31.                             ids,ide, jds,jde, kds,kde,    & 
    32. 

  32.                             ims,ime, jms,jme, kms,kme,    &
    33. 

  33.                             its,ite, jts,jte, kts,kte    )
    34. 

  34. 

  35.   IMPLICIT NONE  ! religion first
    36. 

  36. 

  37. ! declarations for the stuff coming in
    38. 

  38. 

  39.   INTEGER,      INTENT(IN   )    :: ids,ide, jds,jde, kds,kde
    40. 

  40.   INTEGER,      INTENT(IN   )    :: ims,ime, jms,jme, kms,kme
    41. 

  41.   INTEGER,      INTENT(IN   )    :: its,ite, jts,jte, kts,kte
    42. 

  42. 

  43.   INTEGER,      INTENT(IN   )    :: rk_step
    44. 

  44. 

  45.   REAL, DIMENSION(ims:ime, kms:kme, jms:jme),INTENT(INOUT) :: u_1,   &
    46. 

  46.                                                               v_1,   &
    47. 

  47.                                                               w_1,   &
    48. 

  48.                                                               t_1,   &
    49. 

  49.                                                               ph_1
    50. 

  50. 

  51.   REAL, DIMENSION(ims:ime, kms:kme, jms:jme),INTENT(  OUT) :: u_save,   &
    52. 

  52.                                                               v_save,   &
    53. 

  53.                                                               w_save,   &
    54. 

  54.                                                               t_save,   &
    55. 

  55.                                                               ph_save
    56. 

  56. 

  57.   REAL, DIMENSION(ims:ime, kms:kme, jms:jme),INTENT(INOUT) :: u_2,   &
    58. 

  58.                                                               v_2,   &
    59. 

  59.                                                               w_2,   &
    60. 

  60.                                                               t_2,   &
    61. 

  61.                                                               ph_2
    62. 

  62. 

  63.   REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(  OUT) :: c2a, &
    64. 

  64.                                                                ww_save
    65. 

  65. 

  66.   REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN   ) ::  pb,  &
    67. 

  67.                                                                 p,   &
    68. 

  68.                                                                 alt, &
    69. 

  69.                                                                 ww
    70. 

  70. 

  71.   REAL, DIMENSION(ims:ime, jms:jme)         , INTENT(INOUT) :: mu_1,mu_2
    72. 

  72. 

  73.   REAL, DIMENSION(ims:ime, jms:jme)         , INTENT(IN   ) :: mub,  &
    74. 

  74.                                                                muu,  &
    75. 

  75.                                                                muv,  &
    76. 

  76.                                                                mut,  &
    77. 

  77.                                                                msfux,&
    78. 

  78.                                                                msfuy,&
    79. 

  79.                                                                msfvx,&
    80. 

  80.                                                                msfvx_inv,&
    81. 

  81.                                                                msfvy,&
    82. 

  82.                                                                msftx,&
    83. 

  83.                                                                msfty
    84. 

  84. 

  85.   REAL, DIMENSION(ims:ime, jms:jme)         , INTENT(  OUT) :: muus, &
    86. 

  86.                                                                muvs, &
    87. 

  87.                                                                muts, &
    88. 

  88.                                                                mudf
    89. 

  89. 

  90.   REAL, DIMENSION(ims:ime, jms:jme)         , INTENT(  OUT) :: mu_save
    91. 

  91. 

  92.   REAL, DIMENSION(kms:kme, jms:jme)         , INTENT(IN   ) :: dnw
    93. 

  93. 

  94.   REAL, INTENT(IN) :: rdx,rdy
    95. 

  95. 

  96. ! local variables
    97. 

  97. 

  98.   INTEGER :: i, j, k
    99. 

  99.   INTEGER :: i_start, i_end, j_start, j_end, k_start, k_end
    100. 

  100.   INTEGER :: i_endu, j_endv
    101. 

  101. 

  102. 

  103. !<DESCRIPTION>
    104. 

  104. !
    105. 

  105. !  small_step_prep prepares the prognostic variables for the small timestep.
    106. 

  106. !  This includes switching time-levels in the arrays and computing coupled 
    107. 

  107. !  perturbation variables for the small timestep 
    108. 

  108. !  (i.e. mu*u" = mu(t)*u(t)-mu(*)*u(*); mu*u" is advanced during the small 
    109. 

  109. !  timesteps
    110. 

  110. !
    111. 

  111. !</DESCRIPTION>
    112. 

  112. 

  113.     i_start = its
    114. 

  114.     i_end   = min(ite,ide-1)
    115. 

  115.     j_start = jts
    116. 

  116.     j_end   = min(jte,jde-1)
    117. 

  117.     k_start = kts
    118. 

  118.     k_end = min(kte,kde-1)
    119. 

  119. 

  120.     i_endu = ite
    121. 

  121.     j_endv = jte
    122. 

  122. 

  123.     !  if this is the first RK step, reset *_1 to *_2
    124. 

  124.     !  (we are replacing the t-dt fields with the time t fields)
    125. 

  125. 

  126.     IF ((rk_step == 1) ) THEN
    127. 

  127. 

  128.       DO j=j_start, j_end
    129. 

  129.       DO i=i_start, i_end
    130. 

  130.         mu_1(i,j)=mu_2(i,j)
    131. 

  131.         ww_save(i,kde,j) = 0.
    132. 

  132.         ww_save(i,1,j) = 0.
    133. 

  133.         mudf(i,j) = 0.  !  initialize external mode div damp to zero
    134. 

  134.       ENDDO
    135. 

  135.       ENDDO
    136. 

  136. 

  137.       DO j=j_start, j_end
    138. 

  138.       DO k=k_start, k_end
    139. 

  139.       DO i=i_start, i_endu
    140. 

  140.         u_1(i,k,j) = u_2(i,k,j)
    141. 

  141.       ENDDO
    142. 

  142.       ENDDO
    143. 

  143.       ENDDO
    144. 

  144. 

  145.       DO j=j_start, j_endv
    146. 

  146.       DO k=k_start, k_end
    147. 

  147.       DO i=i_start, i_end
    148. 

  148.         v_1(i,k,j) = v_2(i,k,j)
    149. 

  149.       ENDDO
    150. 

  150.       ENDDO
    151. 

  151.       ENDDO
    152. 

  152. 

  153.       DO j=j_start, j_end
    154. 

  154.       DO k=k_start, k_end
    155. 

  155.       DO i=i_start, i_end
    156. 

  156.         t_1(i,k,j) = t_2(i,k,j)
    157. 

  157.       ENDDO
    158. 

  158.       ENDDO
    159. 

  159.       ENDDO
    160. 

  160. 

  161.       DO j=j_start, j_end
    162. 

  162.       DO k=k_start, min(kde,kte)
    163. 

  163.       DO i=i_start, i_end
    164. 

  164.         w_1(i,k,j)  = w_2(i,k,j)
    165. 

  165.         ph_1(i,k,j) = ph_2(i,k,j)
    166. 

  166.       ENDDO
    167. 

  167.       ENDDO
    168. 

  168.       ENDDO
    169. 

  169. 

  170.     DO j=j_start, j_end
    171. 

  171.       DO i=i_start, i_end
    172. 

  172.         muts(i,j)=mub(i,j)+mu_2(i,j)
    173. 

  173.       ENDDO
    174. 

  174.       DO i=i_start, i_endu
    175. 

  175. !  rk_step==1, WCS fix for tiling
    176. 

  176. !        muus(i,j)=0.5*(mub(i,j)+mu_2(i,j)+mub(i-1,j)+mu_2(i-1,j))
    177. 

  177.         muus(i,j) = muu(i,j)
    178. 

  178.       ENDDO
    179. 

  179.     ENDDO
    180. 

  180. 

  181.     DO j=j_start, j_endv
    182. 

  182.     DO i=i_start, i_end
    183. 

  183. !  rk_step==1, WCS fix for tiling
    184. 

  184. !      muvs(i,j)=0.5*(mub(i,j)+mu_2(i,j)+mub(i,j-1)+mu_2(i,j-1))
    185. 

  185.         muvs(i,j) = muv(i,j)
    186. 

  186.     ENDDO
    187. 

  187.     ENDDO
    188. 

  188. 

  189.     DO j=j_start, j_end
    190. 

  190.       DO i=i_start, i_end
    191. 

  191.         mu_save(i,j)=mu_2(i,j)
    192. 

  192.         mu_2(i,j)=0.
    193. 

  193. !       mu_2(i,j)=mu_2(i,j)-mu_2(i,j)
    194. 

  194.       ENDDO
    195. 

  195.     ENDDO
    196. 

  196. 

  197.     ELSE
    198. 

  198. 

  199.     DO j=j_start, j_end
    200. 

  200.       DO i=i_start, i_end
    201. 

  201.         muts(i,j)=mub(i,j)+mu_1(i,j)
    202. 

  202.       ENDDO
    203. 

  203.       DO i=i_start, i_endu
    204. 

  204.         muus(i,j)=0.5*(mub(i,j)+mu_1(i,j)+mub(i-1,j)+mu_1(i-1,j))
    205. 

  205.       ENDDO
    206. 

  206.     ENDDO
    207. 

  207. 

  208.     DO j=j_start, j_endv
    209. 

  209.     DO i=i_start, i_end
    210. 

  210.       muvs(i,j)=0.5*(mub(i,j)+mu_1(i,j)+mub(i,j-1)+mu_1(i,j-1))
    211. 

  211.     ENDDO
    212. 

  212.     ENDDO
    213. 

  213. 

  214.     DO j=j_start, j_end
    215. 

  215.       DO i=i_start, i_end
    216. 

  216.         mu_save(i,j)=mu_2(i,j)
    217. 

  217.         mu_2(i,j)=mu_1(i,j)-mu_2(i,j)
    218. 

  218.       ENDDO
    219. 

  219.     ENDDO
    220. 

  220. 

  221. 

  222.     END IF
    223. 

  223. 

  224.     ! set up the small timestep variables
    225. 

  225. 

  226.       DO j=j_start, j_end
    227. 

  227.       DO i=i_start, i_end
    228. 

  228.         ww_save(i,kde,j) = 0.
    229. 

  229.         ww_save(i,1,j) = 0.
    230. 

  230.       ENDDO
    231. 

  231.       ENDDO
    232. 

  232. 

  233.     DO j=j_start, j_end
    234. 

  234.     DO k=k_start, k_end
    235. 

  235.     DO i=i_start, i_end
    236. 

  236.       c2a(i,k,j) = cpovcv*(pb(i,k,j)+p(i,k,j))/alt(i,k,j)
    237. 

  237.     ENDDO
    238. 

  238.     ENDDO
    239. 

  239.     ENDDO
    240. 

  240. 

  241.     DO j=j_start, j_end
    242. 

  242.     DO k=k_start, k_end
    243. 

  243.     DO i=i_start, i_endu
    244. 

  244.       u_save(i,k,j) = u_2(i,k,j)
    245. 

  245.       ! u coupled with my
    246. 

  246.       u_2(i,k,j) = (muus(i,j)*u_1(i,k,j)-muu(i,j)*u_2(i,k,j))/msfuy(i,j)
    247. 

  247.     ENDDO
    248. 

  248.     ENDDO
    249. 

  249.     ENDDO
    250. 

  250. 

  251.     DO j=j_start, j_endv
    252. 

  252.     DO k=k_start, k_end
    253. 

  253.     DO i=i_start, i_end
    254. 

  254.       v_save(i,k,j) = v_2(i,k,j)
    255. 

  255.       ! v coupled with mx
    256. 

  256. !      v_2(i,k,j) = (muvs(i,j)*v_1(i,k,j)-muv(i,j)*v_2(i,k,j))/msfvx(i,j)
    257. 

  257.       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)
    258. 

  258.     ENDDO
    259. 

  259.     ENDDO
    260. 

  260.     ENDDO
    261. 

  261. 

  262.     DO j=j_start, j_end
    263. 

  263.     DO k=k_start, k_end
    264. 

  264.     DO i=i_start, i_end
    265. 

  265.       t_save(i,k,j) = t_2(i,k,j)
    266. 

  266.       t_2(i,k,j) = muts(i,j)*t_1(i,k,j)-mut(i,j)*t_2(i,k,j)
    267. 

  267.     ENDDO
    268. 

  268.     ENDDO
    269. 

  269.     ENDDO
    270. 

  270. 

  271.     DO j=j_start, j_end
    272. 

  272. !    DO k=k_start, min(kde,kte)
    273. 

  273.     DO k=k_start, kde
    274. 

  274.     DO i=i_start, i_end
    275. 

  275.       w_save(i,k,j) = w_2(i,k,j)
    276. 

  276.       ! w coupled with my
    277. 

  277.       w_2(i,k,j)  = (muts(i,j)* w_1(i,k,j)-mut(i,j)* w_2(i,k,j))/msfty(i,j)
    278. 

  278.       ph_save(i,k,j) = ph_2(i,k,j)
    279. 

  279.       ph_2(i,k,j) = ph_1(i,k,j)-ph_2(i,k,j)
    280. 

  280.     ENDDO
    281. 

  281.     ENDDO
    282. 

  282.     ENDDO
    283. 

  283. 

  284.       DO j=j_start, j_end
    285. 

  285. !      DO k=k_start, min(kde,kte)
    286. 

  286.       DO k=k_start, kde
    287. 

  287.       DO i=i_start, i_end
    288. 

  288.         ww_save(i,k,j) = ww(i,k,j)
    289. 

  289.       ENDDO
    290. 

  290.       ENDDO
    291. 

  291.       ENDDO
    292. 

  292. 

  293. END SUBROUTINE small_step_prep
    294. 

  294. 

  295. !-------------------------------------------------------------------------
    296. 

  296. 

  297. 

  298. SUBROUTINE small_step_finish( u_2, u_1, v_2, v_1, w_2, w_1,    &
    299. 

  299.                               t_2, t_1, ph_2, ph_1, ww, ww1,   &
    300. 

  300.                               mu_2, mu_1,                      &
    301. 

  301.                               mut, muts, muu, muus, muv, muvs, &
    302. 

  302.                               u_save, v_save, w_save,          &
    303. 

  303.                               t_save, ph_save, mu_save,        &
    304. 

  304.                               msfux, msfuy, msfvx, msfvy,      &
    305. 

  305.                               msftx, msfty,                    &
    306. 

  306.                               h_diabatic,                      &
    307. 

  307.                               number_of_small_timesteps,dts,   &
    308. 

  308.                               rk_step, rk_order,               &
    309. 

  309.                               ids,ide, jds,jde, kds,kde,       &
    310. 

  310.                               ims,ime, jms,jme, kms,kme,       &
    311. 

  311.                               its,ite, jts,jte, kts,kte       )
    312. 

  312. 

  313. 

  314.   IMPLICIT NONE  ! religion first
    315. 

  315. 

  316. !  stuff passed in
    317. 

  317. 

  318.   INTEGER,                  INTENT(IN   ) :: ids,ide, jds,jde, kds,kde
    319. 

  319.   INTEGER,                  INTENT(IN   ) :: ims,ime, jms,jme, kms,kme
    320. 

  320.   INTEGER,                  INTENT(IN   ) :: its,ite, jts,jte, kts,kte
    321. 

  321.   INTEGER,                  INTENT(IN   ) :: number_of_small_timesteps
    322. 

  322.   INTEGER,                  INTENT(IN   ) :: rk_step, rk_order
    323. 

  323.   REAL,                     INTENT(IN   ) :: dts
    324. 

  324. 

  325.   REAL,   DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN   ) :: u_1, &
    326. 

  326.                                                                  v_1, &
    327. 

  327.                                                                  w_1, &
    328. 

  328.                                                                  t_1, &
    329. 

  329.                                                                  ww1, &
    330. 

  330.                                                                  ph_1
    331. 

  331. 

  332.   REAL,   DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: u_2, &
    333. 

  333.                                                                  v_2, &
    334. 

  334.                                                                  w_2, &
    335. 

  335.                                                                  t_2, &
    336. 

  336.                                                                  ww,  &
    337. 

  337.                                                                  ph_2
    338. 

  338. 

  339.   REAL, DIMENSION(ims:ime, kms:kme, jms:jme),INTENT(IN   ) :: u_save,   &
    340. 

  340.                                                               v_save,   &
    341. 

  341.                                                               w_save,   &
    342. 

  342.                                                               t_save,   &
    343. 

  343.                                                               ph_save,  &
    344. 

  344.                                                              h_diabatic
    345. 

  345. 

  346.   REAL,   DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: muus, muvs
    347. 

  347.   REAL,   DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: mu_2, mu_1
    348. 

  348.   REAL,   DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: mut, muts, &
    349. 

  349.                                                         muu, muv, mu_save
    350. 

  350.   REAL,   DIMENSION(ims:ime, jms:jme), INTENT(IN   ) :: msfux, msfuy, &
    351. 

  351.                                                         msfvx, msfvy, &
    352. 

  352.                                                         msftx, msfty
    353. 

  353. 

  354. 

  355. ! local stuff
    356. 

  356. 

  357.   INTEGER         :: i,j,k
    358. 

  358.   INTEGER :: i_start, i_end, j_start, j_end, i_endu, j_endv
    359. 

  359. 

  360. !<DESCRIPTION>
    361. 

  361. !
    362. 

  362. !  small_step_finish reconstructs the full uncoupled prognostic variables
    363. 

  363. !  from the coupled perturbation variables used in the small timesteps.
    364. 

  364. !
    365. 

  365. !</DESCRIPTION>
    366. 

  366. 

  367.   i_start = its
    368. 

  368.   i_end   = min(ite,ide-1)
    369. 

  369.   j_start = jts
    370. 

  370.   j_end   = min(jte,jde-1)
    371. 

  371. 

  372.   i_endu = ite
    373. 

  373.   j_endv = jte
    374. 

  374. 

  375. !    addition of time level t back into variables
    376. 

  376. 

  377.   DO j = j_start, j_endv
    378. 

  378.   DO k = kds, kde-1
    379. 

  379.   DO i = i_start, i_end
    380. 

  380.     ! v coupled with mx
    381. 

  381.     v_2(i,k,j) = (msfvx(i,j)*v_2(i,k,j) + v_save(i,k,j)*muv(i,j))/muvs(i,j)
    382. 

  382.   ENDDO
    383. 

  383.   ENDDO
    384. 

  384.   ENDDO
    385. 

  385. 

  386.   DO j = j_start, j_end
    387. 

  387.   DO k = kds, kde-1
    388. 

  388.   DO i = i_start, i_endu
    389. 

  389.     ! u coupled with my
    390. 

  390.     u_2(i,k,j) = (msfuy(i,j)*u_2(i,k,j) + u_save(i,k,j)*muu(i,j))/muus(i,j)
    391. 

  391.   ENDDO
    392. 

  392.   ENDDO
    393. 

  393.   ENDDO
    394. 

  394. 

  395.   DO j = j_start, j_end
    396. 

  396.   DO k = kds, kde
    397. 

  397.   DO i = i_start, i_end
    398. 

  398.     ! w coupled with my
    399. 

  399.     w_2(i,k,j) = (msfty(i,j)*w_2(i,k,j) + w_save(i,k,j)*mut(i,j))/muts(i,j)
    400. 

  400.     ph_2(i,k,j) = ph_2(i,k,j) + ph_save(i,k,j)
    401. 

  401.     ww(i,k,j) = ww(i,k,j) + ww1(i,k,j)
    402. 

  402.   ENDDO
    403. 

  403.   ENDDO
    404. 

  404.   ENDDO
    405. 

  405. 

  406. #ifdef REVERT
    407. 

  407.   DO j = j_start, j_end
    408. 

  408.   DO k = kds, kde-1
    409. 

  409.   DO i = i_start, i_end
    410. 

  410.     t_2(i,k,j) = (t_2(i,k,j) + t_save(i,k,j)*mut(i,j))/muts(i,j)
    411. 

  411.   ENDDO
    412. 

  412.   ENDDO
    413. 

  413.   ENDDO
    414. 

  414. #else
    415. 

  415.   IF ( rk_step < rk_order ) THEN
    416. 

  416.     DO j = j_start, j_end
    417. 

  417.     DO k = kds, kde-1
    418. 

  418.     DO i = i_start, i_end
    419. 

  419.       t_2(i,k,j) = (t_2(i,k,j) + t_save(i,k,j)*mut(i,j))/muts(i,j)
    420. 

  420.     ENDDO
    421. 

  421.     ENDDO
    422. 

  422.     ENDDO
    423. 

  423.   ELSE
    424. 

  424. 

  425.     DO j = j_start, j_end
    426. 

  426.     DO k = kds, kde-1
    427. 

  427.     DO i = i_start, i_end
    428. 

  428.       t_2(i,k,j) = (t_2(i,k,j) - dts*number_of_small_timesteps*mut(i,j)*h_diabatic(i,k,j) &
    429. 

  429.                     + t_save(i,k,j)*mut(i,j))/muts(i,j)
    430. 

  430.     ENDDO
    431. 

  431.     ENDDO
    432. 

  432.     ENDDO
    433. 

  433.   ENDIF
    434. 

  434. #endif
    435. 

  435. 

  436.   DO j = j_start, j_end
    437. 

  437.   DO i = i_start, i_end
    438. 

  438.     mu_2(i,j) = mu_2(i,j) + mu_save(i,j)
    439. 

  439.   ENDDO
    440. 

  440.   ENDDO
    441. 

  441. 

  442. END SUBROUTINE small_step_finish
    443. 

  443. 

  444. !-----------------------------------------------------------------------
    445. 

  445. 

  446. SUBROUTINE calc_p_rho( al, p, ph,                    &
    447. 

  447.                        alt, t_2, t_1, c2a, pm1,      &
    448. 

  448.                        mu, muts, znu, t0,            &
    449. 

  449.                        rdnw, dnw, smdiv,             &
    450. 

  450.                        non_hydrostatic, step,        &
    451. 

  451.                        ids, ide, jds, jde, kds, kde, &
    452. 

  452.                        ims, ime, jms, jme, kms, kme, &
    453. 

  453.                        its,ite, jts,jte, kts,kte    )
    454. 

  454. 

  455.   IMPLICIT NONE  ! religion first
    456. 

  456. 

  457. ! declarations for the stuff coming in
    458. 

  458. 

  459.   INTEGER,      INTENT(IN   )    :: ids,ide, jds,jde, kds,kde
    460. 

  460.   INTEGER,      INTENT(IN   )    :: ims,ime, jms,jme, kms,kme
    461. 

  461.   INTEGER,      INTENT(IN   )    :: its,ite, jts,jte, kts,kte
    462. 

  462. 

  463.   INTEGER,      INTENT(IN   )    :: step
    464. 

  464. 

  465.   REAL, DIMENSION(ims:ime, kms:kme, jms:jme),INTENT(  OUT) :: al,   &
    466. 

  466.                                                                p
    467. 

  467. 

  468.   REAL, DIMENSION(ims:ime, kms:kme, jms:jme),INTENT(IN   ) :: alt,   &
    469. 

  469.                                                               t_2,   &
    470. 

  470.                                                               t_1,   &
    471. 

  471.                                                               c2a
    472. 

  472. 

  473.   REAL, DIMENSION(ims:ime, kms:kme, jms:jme),INTENT(INOUT) :: ph, pm1
    474. 

  474. 

  475.   REAL, DIMENSION(ims:ime, jms:jme)         , INTENT(IN   ) :: mu,   &
    476. 

  476.                                                                muts
    477. 

  477. 

  478.   REAL, DIMENSION(kms:kme)         , INTENT(IN   ) :: dnw,  &
    479. 

  479.                                                       rdnw, &
    480. 

  480.                                                       znu
    481. 

  481. 

  482.   REAL,                                       INTENT(IN   ) :: t0, smdiv
    483. 

  483. 

  484.   LOGICAL, INTENT(IN   )  :: non_hydrostatic
    485. 

  485. 

  486. ! local variables
    487. 

  487. 

  488.   INTEGER :: i, j, k
    489. 

  489.   INTEGER :: i_start, i_end, j_start, j_end, k_start, k_end
    490. 

  490.   REAL    :: ptmp
    491. 

  491. 

  492. !<DESCRIPTION>
    493. 

  493. !
    494. 

  494. !  For the nonhydrostatic option,
    495. 

  495. !  calc_p_rho computes the perturbation inverse density and 
    496. 

  496. !  perturbation pressure from the hydrostatic relation and the 
    497. 

  497. !  linearized equation of state, respectively.
    498. 

  498. !
    499. 

  499. !  For the hydrostatic option,
    500. 

  500. !  calc_p_rho computes the perturbation pressure, perturbation density,
    501. 

  501. !  and perturbation geopotential
    502. 

  502. !  from the vertical coordinate definition, linearized equation of state 
    503. 

  503. !  and the hydrostatic relation, respectively.
    504. 

  504. !
    505. 

  505. !  forward weighting of the pressure (divergence damping) is also
    506. 

  506. !  computed here.
    507. 

  507. !
    508. 

  508. !</DESCRIPTION>
    509. 

  509. 

  510.    i_start = its
    511. 

  511.    i_end   = min(ite,ide-1)
    512. 

  512.    j_start = jts
    513. 

  513.    j_end   = min(jte,jde-1)
    514. 

  514.    k_start = kts
    515. 

  515.    k_end = min(kte,kde-1)
    516. 

  516. 

  517.    IF (non_hydrostatic) THEN
    518. 

  518.      DO j=j_start, j_end
    519. 

  519.      DO k=k_start, k_end
    520. 

  520.      DO i=i_start, i_end
    521. 

  521. 

  522. !  al computation is all dry, so ok with moisture
    523. 

  523. 

  524.       al(i,k,j)=-1./muts(i,j)*(alt(i,k,j)*mu(i,j)  &
    525. 

  525.              +rdnw(k)*(ph(i,k+1,j)-ph(i,k,j)))
    526. 

  526. 

  527. !  this is temporally linearized p, no moisture correction needed
    528. 

  528. 

  529.       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))  &
    530. 

  530.                        /(muts(i,j)*(t0+t_1(i,k,j)))-al (i,k,j))
    531. 

  531. 

  532.      ENDDO
    533. 

  533.      ENDDO
    534. 

  534.      ENDDO
    535. 

  535. 

  536.    ELSE  ! hydrostatic calculation
    537. 

  537. 

  538.        DO j=j_start, j_end
    539. 

  539.        DO k=k_start, k_end
    540. 

  540.        DO i=i_start, i_end
    541. 

  541.          p(i,k,j)=mu(i,j)*znu(k)
    542. 

  542.          al(i,k,j)=alt(i,k,j)*(t_2(i,k,j)-mu(i,j)*t_1(i,k,j))            &
    543. 

  543.                       /(muts(i,j)*(t0+t_1(i,k,j)))-p(i,k,j)/c2a(i,k,j)
    544. 

  544.          ph(i,k+1,j)=ph(i,k,j)-dnw(k)*(muts(i,j)*al (i,k,j)              &
    545. 

  545.                           +mu(i,j)*alt(i,k,j))
    546. 

  546.        ENDDO
    547. 

  547.        ENDDO
    548. 

  548.        ENDDO
    549. 

  549. 

  550.    END IF
    551. 

  551. 

  552. !  divergence damping setup
    553. 

  553.  
    554. 

  554.      IF (step == 0) then   ! we're initializing small timesteps
    555. 

  555.        DO j=j_start, j_end
    556. 

  556.        DO k=k_start, k_end
    557. 

  557.        DO i=i_start, i_end
    558. 

  558.          pm1(i,k,j)=p(i,k,j)
    559. 

  559.        ENDDO
    560. 

  560.        ENDDO
    561. 

  561.        ENDDO 
    562. 

  562.      ELSE                     ! we're in the small timesteps 
    563. 

  563.        DO j=j_start, j_end    ! and adding div damping component
    564. 

  564.        DO k=k_start, k_end
    565. 

  565.        DO i=i_start, i_end
    566. 

  566.          ptmp = p(i,k,j)
    567. 

  567.          p(i,k,j) = p(i,k,j) + smdiv*(p(i,k,j)-pm1(i,k,j))
    568. 

  568.          pm1(i,k,j) = ptmp
    569. 

  569.        ENDDO
    570. 

  570.        ENDDO
    571. 

  571.        ENDDO
    572. 

  572.      END IF
    573. 

  573. 

  574. END SUBROUTINE calc_p_rho
    575. 

  575. 

  576. !----------------------------------------------------------------------
    577. 

  577. 

  578. SUBROUTINE calc_coef_w( a,alpha,gamma,              &
    579. 

  579.                         mut, cqw,                   &
    580. 

  580.                         rdn, rdnw, c2a,             &
    581. 

  581.                         dts, g, epssm, top_lid,     &
    582. 

  582.                         ids,ide, jds,jde, kds,kde,  & ! domain dims
    583. 

  583.                         ims,ime, jms,jme, kms,kme,  & ! memory dims
    584. 

  584.                         its,ite, jts,jte, kts,kte  )  ! tile   dims
    585. 

  585.                                                    
    586. 

  586.       IMPLICIT NONE  ! religion first
    587. 

  587. 

  588. !  passed in through the call
    589. 

  589. 

  590.   INTEGER,      INTENT(IN   )    :: ids,ide, jds,jde, kds,kde
    591. 

  591.   INTEGER,      INTENT(IN   )    :: ims,ime, jms,jme, kms,kme
    592. 

  592.   INTEGER,      INTENT(IN   )    :: its,ite, jts,jte, kts,kte
    593. 

  593. 

  594.   LOGICAL,      INTENT(IN   )    :: top_lid
    595. 

  595. 

  596.   REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN   ) :: c2a,  &
    597. 

  597.                                                                cqw
    598. 

  598. 

  599.   REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: alpha, &
    600. 

  600.                                                                gamma, &
    601. 

  601.                                                                a
    602. 

  602. 

  603.   REAL, DIMENSION(ims:ime, jms:jme),          INTENT(IN   ) :: mut
    604. 

  604. 

  605.   REAL, DIMENSION(kms:kme),                   INTENT(IN   ) :: rdn,   &
    606. 

  606.                                                                rdnw
    607. 

  607. 

  608.   REAL,                                       INTENT(IN   ) :: epssm, &
    609. 

  609.                                                                dts,   &
    610. 

  610.                                                                g
    611. 

  611. 

  612. !  Local stack data.
    613. 

  613. 

  614.   REAL, DIMENSION(ims:ime)                         :: cof
    615. 

  615.   REAL  :: b, c
    616. 

  616. 

  617.   INTEGER :: i, j, k, i_start, i_end, j_start, j_end, k_start, k_end
    618. 

  618.   INTEGER :: ij, ijp, ijm, lid_flag
    619. 

  619. 

  620. !<DESCRIPTION>
    621. 

  621. !
    622. 

  622. !  calc_coef_w calculates the coefficients needed for the 
    623. 

  623. !  implicit solution of the vertical momentum and geopotential equations.
    624. 

  624. !  This requires solution of a tri-diagonal equation.
    625. 

  625. !
    626. 

  626. !</DESCRIPTION>
    627. 

  627. 

  628.       i_start = its
    629. 

  629.       i_end   = min(ite,ide-1)
    630. 

  630.       j_start = jts
    631. 

  631.       j_end   = min(jte,jde-1)
    632. 

  632.       k_start = kts
    633. 

  633.       k_end   = kte-1
    634. 

  634. 

  635.       lid_flag=1
    636. 

  636.       IF(top_lid)lid_flag=0
    637. 

  637. 

  638.      outer_j_loop:  DO j = j_start, j_end
    639. 

  639. 

  640.         DO i = i_start, i_end
    641. 

  641.           cof(i)  = (.5*dts*g*(1.+epssm)/mut(i,j))**2
    642. 

  642.           a(i, 2 ,j) = 0.
    643. 

  643.           a(i,kde,j) =-2.*cof(i)*rdnw(kde-1)**2*c2a(i,kde-1,j)*lid_flag
    644. 

  644.           gamma(i,1 ,j) = 0.
    645. 

  645.         ENDDO
    646. 

  646. 

  647.         DO k=3,kde-1
    648. 

  648.         DO i=i_start, i_end
    649. 

  649.           a(i,k,j) =   -cqw(i,k,j)*cof(i)*rdn(k)* rdnw(k-1)*c2a(i,k-1,j)    
    650. 

  650.         ENDDO
    651. 

  651.         ENDDO
    652. 

  652. 

  653. 

  654.         DO k=2,kde-1
    655. 

  655.         DO i=i_start, i_end
    656. 

  656.           b = 1.+cqw(i,k,j)*cof(i)*rdn(k)*(rdnw(k  )*c2a(i,k,j  )  &
    657. 

  657.                                        +rdnw(k-1)*c2a(i,k-1,j))
    658. 

  658.           c =   -cqw(i,k,j)*cof(i)*rdn(k)*rdnw(k  )*c2a(i,k,j  )
    659. 

  659.           alpha(i,k,j) = 1./(b-a(i,k,j)*gamma(i,k-1,j))
    660. 

  660.           gamma(i,k,j) = c*alpha(i,k,j)
    661. 

  661.         ENDDO
    662. 

  662.         ENDDO
    663. 

  663. 

  664.         DO i=i_start, i_end
    665. 

  665.           b = 1.+2.*cof(i)*rdnw(kde-1)**2*c2a(i,kde-1,j)
    666. 

  666.           c = 0.
    667. 

  667.           alpha(i,kde,j) = 1./(b-a(i,kde,j)*gamma(i,kde-1,j))
    668. 

  668.           gamma(i,kde,j) = c*alpha(i,kde,j)
    669. 

  669.         ENDDO
    670. 

  670. 

  671.       ENDDO outer_j_loop
    672. 

  672. 

  673. END SUBROUTINE calc_coef_w
    674. 

  674. 

  675. !----------------------------------------------------------------------
    676. 

  676. 

  677. SUBROUTINE advance_uv ( u, ru_tend, v, rv_tend,        &
    678. 

  678.                         p, pb,                         &
    679. 

  679.                         ph, php, alt, al, mu,          &
    680. 

  680.                         muu, cqu, muv, cqv, mudf,      &
    681. 

  681.                         msfux, msfuy, msfvx,           &
    682. 

  682.                         msfvx_inv, msfvy,              &
    683. 

  683.                         rdx, rdy, dts,                 &
    684. 

  684.                         cf1, cf2, cf3, fnm, fnp,       &
    685. 

  685.                         emdiv,                         &
    686. 

  686.                         rdnw, config_flags, spec_zone, &
    687. 

  687.                         non_hydrostatic, top_lid,      &
    688. 

  688.                         ids, ide, jds, jde, kds, kde,  &
    689. 

  689.                         ims, ime, jms, jme, kms, kme,  &
    690. 

  690.                         its, ite, jts, jte, kts, kte  )
    691. 

  691. 

  692. 

  693. 

  694.       IMPLICIT NONE  ! religion first
    695. 

  695. 

  696. ! stuff coming in
    697. 

  697. 

  698.       TYPE(grid_config_rec_type), INTENT(IN   ) :: config_flags
    699. 

  699. 

  700.       LOGICAL, INTENT(IN   ) :: non_hydrostatic, top_lid
    701. 

  701. 

  702.       INTEGER,      INTENT(IN   )    :: ids,ide, jds,jde, kds,kde
    703. 

  703.       INTEGER,      INTENT(IN   )    :: ims,ime, jms,jme, kms,kme
    704. 

  704.       INTEGER,      INTENT(IN   )    :: its,ite, jts,jte, kts,kte
    705. 

  705.       INTEGER,      INTENT(IN   )    :: spec_zone
    706. 

  706. 

  707.       REAL, DIMENSION( ims:ime , kms:kme, jms:jme ),  &
    708. 

  708.             INTENT(INOUT) ::                          &
    709. 

  709.                                                   u,  &
    710. 

  710.                                                   v
    711. 

  711. 

  712.       REAL, DIMENSION( ims:ime , kms:kme, jms:jme ), &
    713. 

  713.             INTENT(IN   ) ::                          &
    714. 

  714.                                              ru_tend, &
    715. 

  715.                                              rv_tend, &
    716. 

  716.                                              ph,      &
    717. 

  717.                                              php,     &
    718. 

  718.                                              p,       &
    719. 

  719.                                              pb,      &
    720. 

  720.                                              alt,     &
    721. 

  721.                                              al,      &
    722. 

  722.                                              cqu,     &
    723. 

  723.                                              cqv
    724. 

  724. 

  725. 

  726.       REAL, DIMENSION( ims:ime , jms:jme ),    INTENT(IN   ) :: muu,  &
    727. 

  727.                                                                 muv,  &
    728. 

  728.                                                                 mu,   &
    729. 

  729.                                                                 mudf
    730. 

  730. 

  731. 

  732.       REAL, DIMENSION( kms:kme ),              INTENT(IN   ) :: fnm,    &
    733. 

  733.                                                                 fnp ,   &
    734. 

  734.                                                                 rdnw
    735. 

  735. 

  736.       REAL, DIMENSION( ims:ime , jms:jme ),    INTENT(IN   ) :: msfux,  &
    737. 

  737.                                                                 msfuy,  &
    738. 

  738.                                                                 msfvx,  &
    739. 

  739.                                                                 msfvy,  &
    740. 

  740.                                                                 msfvx_inv
    741. 

  741. 

  742.       REAL,                                    INTENT(IN   ) :: rdx,    &
    743. 

  743.                                                                 rdy,    &
    744. 

  744.                                                                 dts,    &
    745. 

  745.                                                                 cf1,    &
    746. 

  746.                                                                 cf2,    &
    747. 

  747.                                                                 cf3,    &
    748. 

  748.                                                                 emdiv
    749. 

  749.     
    750. 

  750. 

  751. !  Local 3d array from the stack (note tile size)
    752. 

  752. 

  753.       REAL, DIMENSION (its:ite, kts:kte) :: dpn, dpxy
    754. 

  754.       REAL, DIMENSION (its:ite) :: mudf_xy
    755. 

  755.       REAL                      :: dx, dy
    756. 

  756. 

  757.       INTEGER :: i,j,k, i_start, i_end, j_start, j_end, k_start, k_end
    758. 

  758.       INTEGER :: i_endu, j_endv, k_endw
    759. 

  759.       INTEGER :: i_start_up, i_end_up, j_start_up, j_end_up
    760. 

  760.       INTEGER :: i_start_vp, i_end_vp, j_start_vp, j_end_vp
    761. 

  761. 

  762.       INTEGER :: i_start_u_tend, i_end_u_tend, j_start_v_tend, j_end_v_tend
    763. 

  763. 

  764. !<DESCRIPTION>
    765. 

  765. !
    766. 

  766. !  advance_uv advances the explicit perturbation horizontal momentum 
    767. 

  767. !  equations (u,v) by adding in the large-timestep tendency along with 
    768. 

  768. !  the small timestep pressure gradient tendency.
    769. 

  769. !
    770. 

  770. !</DESCRIPTION>
    771. 

  771. 

  772. !  now, the real work.
    773. 

  773. !  set the loop bounds taking into account boundary conditions.
    774. 

  774. 

  775.     IF( config_flags%nested .or. config_flags%specified ) THEN
    776. 

  776.       i_start = max( its,ids+spec_zone )
    777. 

  777.       i_end   = min( ite,ide-spec_zone-1 )
    778. 

  778.       j_start = max( jts,jds+spec_zone )
    779. 

  779.       j_end   = min( jte,jde-spec_zone-1 )
    780. 

  780.       k_start = kts
    781. 

  781.       k_end   = min( kte, kde-1 )
    782. 

  782. 

  783.       i_endu = min( ite,ide-spec_zone )
    784. 

  784.       j_endv = min( jte,jde-spec_zone )
    785. 

  785.       k_endw = k_end
    786. 

  786. 

  787.       IF( config_flags%periodic_x) THEN
    788. 

  788.         i_start = its
    789. 

  789.         i_end   = min(ite,ide-1)
    790. 

  790.         i_endu = ite
    791. 

  791.       ENDIF
    792. 

  792.     ELSE
    793. 

  793.       i_start = its
    794. 

  794.       i_end   = min(ite,ide-1)
    795. 

  795.       j_start = jts
    796. 

  796.       j_end   = min(jte,jde-1)
    797. 

  797.       k_start = kts
    798. 

  798.       k_end   = kte-1
    799. 

  799. 

  800.       i_endu = ite
    801. 

  801.       j_endv = jte
    802. 

  802.       k_endw = k_end
    803. 

  803.     ENDIF
    804. 

  804. 

  805.       i_start_up = i_start
    806. 

  806.       i_end_up   = i_endu
    807. 

  807.       j_start_up = j_start
    808. 

  808.       j_end_up   = j_end
    809. 

  809. 

  810.       i_start_vp = i_start
    811. 

  811.       i_end_vp   = i_end
    812. 

  812.       j_start_vp = j_start
    813. 

  813.       j_end_vp   = j_endv
    814. 

  814. 

  815.       IF ( (config_flags%open_xs   .or.     &
    816. 

  816.             config_flags%symmetric_xs   )   &
    817. 

  817.             .and. (its == ids) )            &
    818. 

  818.                  i_start_up = i_start_up + 1
    819. 

  819. 

  820.       IF ( (config_flags%open_xe    .or.  &
    821. 

  821.             config_flags%symmetric_xe   ) &
    822. 

  822.              .and. (ite == ide) )         &
    823. 

  823.                  i_end_up   = i_end_up - 1
    824. 

  824. 

  825.       IF ( (config_flags%open_ys    .or.   &
    826. 

  826.             config_flags%symmetric_ys  .or.   &
    827. 

  827.             config_flags%polar   )  &
    828. 

  828.                      .and. (jts == jds) )  &
    829. 

  829.                  j_start_vp = j_start_vp + 1
    830. 

  830. 

  831.       IF ( (config_flags%open_ye     .or. &
    832. 

  832.             config_flags%symmetric_ye  .or.   &
    833. 

  833.             config_flags%polar   )  &
    834. 

  834.             .and. (jte == jde) )          &
    835. 

  835.                  j_end_vp   = j_end_vp - 1
    836. 

  836. 

  837.       i_start_u_tend = i_start
    838. 

  838.       i_end_u_tend   = i_endu
    839. 

  839.       j_start_v_tend = j_start
    840. 

  840.       j_end_v_tend   = j_endv
    841. 

  841. 

  842.       IF ( config_flags%symmetric_xs .and. (its == ids) ) &
    843. 

  843.                      i_start_u_tend = i_start_u_tend+1
    844. 

  844.       IF ( config_flags%symmetric_xe .and. (ite == ide) ) &
    845. 

  845.                      i_end_u_tend = i_end_u_tend-1
    846. 

  846.       IF ( config_flags%symmetric_ys .and. (jts == jds) ) &
    847. 

  847.                      j_start_v_tend = j_start_v_tend+1
    848. 

  848.       IF ( config_flags%symmetric_ye .and. (jte == jde) ) &
    849. 

  849.                      j_end_v_tend = j_end_v_tend-1
    850. 

  850. 

  851.    dx = 1./rdx
    852. 

  852.    dy = 1./rdy
    853. 

  853. 

  854. !  start real calculations.
    855. 

  855. !  first, u
    856. 

  856. 

  857.   u_outer_j_loop: DO j = j_start, j_end
    858. 

  858. 

  859.    DO k = k_start, k_end
    860. 

  860.    DO i = i_start_u_tend, i_end_u_tend
    861. 

  861.      u(i,k,j) = u(i,k,j) + dts*ru_tend(i,k,j)
    862. 

  862.    ENDDO
    863. 

  863.    ENDDO
    864. 

  864. 

  865.    DO i = i_start_up, i_end_up
    866. 

  866.      mudf_xy(i)= -emdiv*dx*(mudf(i,j)-mudf(i-1,j))/msfuy(i,j)
    867. 

  867.    ENDDO
    868. 

  868. 

  869.    DO k = k_start, k_end
    870. 

  870.    DO i = i_start_up, i_end_up
    871. 

  871. 

  872. !     Comments on map scale factors:   
    873. 

  873. !     x pressure gradient: ADT eqn 44, penultimate term on RHS
    874. 

  874. !        = -(mx/my)*(mu/rho)*partial dp/dx
    875. 

  875. !     [i.e., first rho->mu; 2nd still rho; alpha=1/rho]
    876. 

  876. !     Klemp et al. splits into 2 terms:
    877. 

  877. !        mu alpha partial dp/dx + partial dp/dnu * partial dphi/dx 
    878. 

  878. !     then into 4 terms:
    879. 

  879. !        mu alpha partial dp'/dx + nu mu alpha' partial dmubar/dx +
    880. 

  880. !      + mu partial dphi/dx + partial dphi'/dx * (partial dp'/dnu - mu')
    881. 

  881. !
    882. 

  882. !     first 3 terms:
    883. 

  883. !     ph, alt, p, al, pb not coupled
    884. 

  884. !     since we want tendency to fit ADT eqn 44 (coupled) we need to
    885. 

  885. !     multiply by (mx/my):
    886. 

  886. !
    887. 

  887.      dpxy(i,k)= (msfux(i,j)/msfuy(i,j))*.5*rdx*muu(i,j)*(         &
    888. 

  888.        ((ph (i,k+1,j)-ph (i-1,k+1,j))+(ph (i,k,j)-ph (i-1,k,j)))  &
    889. 

  889.       +(alt(i,k  ,j)+alt(i-1,k  ,j))*(p  (i,k,j)-p  (i-1,k,j))  &
    890. 

  890.       +(al (i,k  ,j)+al (i-1,k  ,j))*(pb (i,k,j)-pb (i-1,k,j)) ) 
    891. 

  891. 

  892.    ENDDO
    893. 

  893.    ENDDO
    894. 

  894. 

  895.    IF (non_hydrostatic) THEN
    896. 

  896. 

  897.      DO i = i_start_up, i_end_up
    898. 

  898.        dpn(i,1) = .5*( cf1*(p(i,1,j)+p(i-1,1,j))  &
    899. 

  899.                       +cf2*(p(i,2,j)+p(i-1,2,j))  &
    900. 

  900.                       +cf3*(p(i,3,j)+p(i-1,3,j)) )
    901. 

  901.        dpn(i,kde) = 0.
    902. 

  902.      ENDDO
    903. 

  903.      IF (top_lid) THEN
    904. 

  904.        DO i = i_start_up, i_end_up
    905. 

  905.          dpn(i,kde) =.5*( cf1*(p(i-1,kde-1,j)+p(i,kde-1,j))   &
    906. 

  906.                          +cf2*(p(i-1,kde-2,j)+p(i,kde-2,j))   &
    907. 

  907.                          +cf3*(p(i-1,kde-3,j)+p(i,kde-3,j))  )
    908. 

  908.        ENDDO
    909. 

  909.      ENDIF
    910. 

  910. 

  911.      DO k = k_start+1, k_end
    912. 

  912.        DO i = i_start_up, i_end_up
    913. 

  913.          dpn(i,k) = .5*( fnm(k)*(p(i,k  ,j)+p(i-1,k  ,j))   &
    914. 

  914.                         +fnp(k)*(p(i,k-1,j)+p(i-1,k-1,j)) )
    915. 

  915.        ENDDO
    916. 

  916.      ENDDO
    917. 

  917. 

  918. !    Comments on map scale factors:
    919. 

  919. !    4th term:
    920. 

  920. !    php, dpn, mu not coupled
    921. 

  921. !    since we want tendency to fit ADT eqn 44 (coupled) we need to
    922. 

  922. !    multiply by (mx/my):
    923. 

  923. 

  924.      DO k = k_start, k_end
    925. 

  925.        DO i = i_start_up, i_end_up
    926. 

  926.          dpxy(i,k)=dpxy(i,k) + (msfux(i,j)/msfuy(i,j))*rdx*(php(i,k,j)-php(i-1,k,j))*        &
    927. 

  927.            (rdnw(k)*(dpn(i,k+1)-dpn(i,k))-.5*(mu(i-1,j)+mu(i,j)))
    928. 

  928.        ENDDO
    929. 

  929.      ENDDO
    930. 

  930. 

  931.  
    932. 

  932.    END IF
    933. 

  933. 

  934. 

  935.    DO k = k_start, k_end
    936. 

  936.      DO i = i_start_up, i_end_up
    937. 

  937.        u(i,k,j)=u(i,k,j)-dts*cqu(i,k,j)*dpxy(i,k)+mudf_xy(i) 
    938. 

  938.      ENDDO
    939. 

  939.    ENDDO
    940. 

  940. 

  941.    ENDDO u_outer_j_loop
    942. 

  942. 

  943. ! now v
    944. 

  944. 

  945.   v_outer_j_loop: DO j = j_start_v_tend, j_end_v_tend
    946. 

  946. 

  947. 

  948.    DO k = k_start, k_end
    949. 

  949.    DO i = i_start, i_end
    950. 

  950.      v(i,k,j) = v(i,k,j) + dts*rv_tend(i,k,j)
    951. 

  951.    ENDDO
    952. 

  952.    ENDDO
    953. 

  953. 

  954.    DO i = i_start, i_end
    955. 

  955.      mudf_xy(i)= -emdiv*dy*(mudf(i,j)-mudf(i,j-1))*msfvx_inv(i,j)
    956. 

  956.    ENDDO
    957. 

  957. 

  958.      IF (     ( j >= j_start_vp)  &
    959. 

  959.        .and.( j <= j_end_vp  ) )  THEN
    960. 

  960. 

  961.          DO k = k_start, k_end
    962. 

  962. 	 DO i = i_start, i_end
    963. 

  963. 

  964. !      Comments on map scale factors:
    965. 

  965. !      y pressure gradient: ADT eqn 45, penultimate term on RHS
    966. 

  966. !         = -(my/mx)*(mu/rho)*partial dp/dy
    967. 

  967. !      [i.e., first rho->mu; 2nd still rho; alpha=1/rho]
    968. 

  968. !      Klemp et al. splits into 2 terms:
    969. 

  969. !         mu alpha partial dp/dy + partial dp/dnu * partial dphi/dy 
    970. 

  970. !      then into 4 terms:
    971. 

  971. !         mu alpha partial dp'/dy + nu mu alpha' partial dmubar/dy +
    972. 

  972. !       + mu partial dphi/dy + partial dphi'/dy * (partial dp'/dnu - mu')
    973. 

  973. !
    974. 

  974. !      first 3 terms:
    975. 

  975. !      ph, alt, p, al, pb not coupled
    976. 

  976. !      since we want tendency to fit ADT eqn 45 (coupled) we need to
    977. 

  977. !      multiply by (my/mx):
    978. 

  978. !      mudf_xy is NOT a map scale factor coupling
    979. 

  979. !      it is some sort of divergence damping
    980. 

  980. 

  981.        dpxy(i,k)= (msfvy(i,j)/msfvx(i,j))*.5*rdy*muv(i,j)*(       &
    982. 

  982.         ((ph(i,k+1,j)-ph(i,k+1,j-1))+(ph (i,k,j)-ph (i,k,j-1)))   &
    983. 

  983.         +(alt(i,k  ,j)+alt(i,k  ,j-1))*(p  (i,k,j)-p  (i,k,j-1))  &
    984. 

  984.         +(al (i,k  ,j)+al (i,k  ,j-1))*(pb (i,k,j)-pb (i,k,j-1)) ) 
    985. 

  985. 	    ENDDO
    986. 

  986. 	    ENDDO
    987. 

  987. 

  988. 

  989.      IF (non_hydrostatic) THEN
    990. 

  990. 

  991.        DO i = i_start, i_end     
    992. 

  992.          dpn(i,1) = .5*( cf1*(p(i,1,j)+p(i,1,j-1))  &
    993. 

  993.                         +cf2*(p(i,2,j)+p(i,2,j-1))  &
    994. 

  994.                         +cf3*(p(i,3,j)+p(i,3,j-1)) )
    995. 

  995.          dpn(i,kde) = 0.
    996. 

  996.        ENDDO
    997. 

  997.      IF (top_lid) THEN
    998. 

  998.        DO i = i_start, i_end
    999. 

  999.          dpn(i,kde) =.5*( cf1*(p(i,kde-1,j-1)+p(i,kde-1,j))   &
    1000. 

  1000.                          +cf2*(p(i,kde-2,j-1)+p(i,kde-2,j))   &
    1001. 

  1001.                          +cf3*(p(i,kde-3,j-1)+p(i,kde-3,j))  )
    1002. 

  1002.        ENDDO
    1003. 

  1003.      ENDIF
    1004. 

  1004. 

  1005.        DO k = k_start+1, k_end
    1006. 

  1006.          DO i = i_start, i_end
    1007. 

  1007.            dpn(i,k) = .5*( fnm(k)*(p(i,k  ,j)+p(i,k  ,j-1))  &
    1008. 

  1008.                           +fnp(k)*(p(i,k-1,j)+p(i,k-1,j-1)) )
    1009. 

  1009.          ENDDO
    1010. 

  1010.        ENDDO
    1011. 

  1011. 

  1012. !      Comments on map scale factors:
    1013. 

  1013. !      4th term:
    1014. 

  1014. !      php, dpn, mu not coupled
    1015. 

  1015. !      since we want tendency to fit ADT eqn 45 (coupled) we need to
    1016. 

  1016. !      multiply by (my/mx):
    1017. 

  1017. 

  1018.        DO k = k_start, k_end
    1019. 

  1019.          DO i = i_start, i_end
    1020. 

  1020.            dpxy(i,k)=dpxy(i,k) + (msfvy(i,j)/msfvx(i,j))*rdy*(php(i,k,j)-php(i,k,j-1))*    &
    1021. 

  1021.              (rdnw(k)*(dpn(i,k+1)-dpn(i,k))-.5*(mu(i,j-1)+mu(i,j)))
    1022. 

  1022.          ENDDO
    1023. 

  1023.        ENDDO
    1024. 

  1024. 

  1025. 

  1026.      END IF
    1027. 

  1027. 

  1028. 

  1029.      DO k = k_start, k_end
    1030. 

  1030.        DO i = i_start, i_end
    1031. 

  1031.          v(i,k,j)=v(i,k,j)-dts*cqv(i,k,j)*dpxy(i,k)+mudf_xy(i) 
    1032. 

  1032.        ENDDO
    1033. 

  1033.      ENDDO
    1034. 

  1034.    END IF
    1035. 

  1035. 

  1036.   ENDDO  v_outer_j_loop
    1037. 

  1037. 

  1038.   ! The check for j_start_vp and j_end_vp makes sure that the edges in v
    1039. 

  1039.   ! are not updated.  Let's add a polar boundary condition check here for
    1040. 

  1040.   ! safety to ensure that v at the poles never gets to be non-zero in the
    1041. 

  1041.   ! small time steps.
    1042. 

  1042.   IF (config_flags%polar) THEN
    1043. 

  1043.      IF (jts == jds) THEN
    1044. 

  1044.         DO k = k_start, k_end
    1045. 

  1045.         DO i = i_start, i_end
    1046. 

  1046.            v(i,k,jds) = 0.
    1047. 

  1047.         ENDDO
    1048. 

  1048.         ENDDO
    1049. 

  1049.      END IF
    1050. 

  1050.      IF (jte == jde) THEN
    1051. 

  1051.         DO k = k_start, k_end
    1052. 

  1052.         DO i = i_start, i_end
    1053. 

  1053.            v(i,k,jde) = 0.
    1054. 

  1054.         ENDDO
    1055. 

  1055.         ENDDO
    1056. 

  1056.      END IF
    1057. 

  1057.   END IF
    1058. 

  1058. 

  1059. 

  1060. END SUBROUTINE advance_uv
    1061. 

  1061. 

  1062. !---------------------------------------------------------------------
    1063. 

  1063. 

  1064. SUBROUTINE advance_mu_t( ww, ww_1, u, u_1, v, v_1,            &
    1065. 

  1065.                          mu, mut, muave, muts, muu, muv,      &
    1066. 

  1066.                          mudf, uam, vam, wwam, t, t_1,        &
    1067. 

  1067.                          t_ave, ft, mu_tend,                  &
    1068. 

  1068.                          rdx, rdy, dts, epssm,                &
    1069. 

  1069.                          dnw, fnm, fnp, rdnw,                 &
    1070. 

  1070.                          msfux, msfuy, msfvx, msfvx_inv,      &
    1071. 

  1071.                          msfvy, msftx, msfty,                 &
    1072. 

  1072.                          step, config_flags,                  &
    1073. 

  1073.                          ids, ide, jds, jde, kds, kde,        &
    1074. 

  1074.                          ims, ime, jms, jme, kms, kme,        &
    1075. 

  1075.                          its, ite, jts, jte, kts, kte        )
    1076. 

  1076. 

  1077.   IMPLICIT NONE  ! religion first
    1078. 

  1078. 

  1079. ! stuff coming in
    1080. 

  1080. 

  1081.   TYPE(grid_config_rec_type), INTENT(IN   ) :: config_flags
    1082. 

  1082. 

  1083.   INTEGER,      INTENT(IN   )    :: ids,ide, jds,jde, kds,kde
    1084. 

  1084.   INTEGER,      INTENT(IN   )    :: ims,ime, jms,jme, kms,kme
    1085. 

  1085.   INTEGER,      INTENT(IN   )    :: its,ite, jts,jte, kts,kte
    1086. 

  1086. 

  1087.   INTEGER,      INTENT(IN   )    :: step
    1088. 

  1088. 

  1089.   REAL, DIMENSION( ims:ime , kms:kme, jms:jme ),   &
    1090. 

  1090.             INTENT(IN   ) ::                       &
    1091. 

  1091.                                               u,   &
    1092. 

  1092.                                               v,   &
    1093. 

  1093.                                               u_1, &
    1094. 

  1094.                                               v_1, &
    1095. 

  1095.                                               t_1, &
    1096. 

  1096.                                               ft
    1097. 

  1097. 

  1098.   REAL, DIMENSION( ims:ime , kms:kme, jms:jme ),      &
    1099. 

  1099.             INTENT(INOUT) ::                          &
    1100. 

  1100.                                               ww,     &
    1101. 

  1101.                                               ww_1,   &
    1102. 

  1102.                                               t,      &
    1103. 

  1103.                                               t_ave,  &
    1104. 

  1104.                                               uam,    &
    1105. 

  1105.                                               vam,    &
    1106. 

  1106.                                               wwam
    1107. 

  1107.                                               
    1108. 

  1108.   REAL, DIMENSION( ims:ime , jms:jme ),    INTENT(IN   ) :: muu,  &
    1109. 

  1109.                                                             muv,  &
    1110. 

  1110.                                                             mut,  &
    1111. 

  1111.                                                             msfux,&
    1112. 

  1112.                                                             msfuy,&
    1113. 

  1113.                                                             msfvx,&
    1114. 

  1114.                                                             msfvx_inv,&
    1115. 

  1115.                                                             msfvy,&
    1116. 

  1116.                                                             msftx,&
    1117. 

  1117.                                                             msfty,&
    1118. 

  1118.                                                             mu_tend
    1119. 

  1119. 

  1120.   REAL, DIMENSION( ims:ime , jms:jme ),    INTENT(  OUT) :: muave, &
    1121. 

  1121.                                                             muts,  &
    1122. 

  1122.                                                             mudf
    1123. 

  1123. 

  1124.   REAL, DIMENSION( ims:ime , jms:jme ),    INTENT(INOUT) :: mu
    1125. 

  1125. 

  1126.   REAL, DIMENSION( kms:kme ),              INTENT(IN   ) :: fnm,    &
    1127. 

  1127.                                                             fnp,    &
    1128. 

  1128.                                                             dnw,    &
    1129. 

  1129.                                                             rdnw
    1130. 

  1130. 

  1131. 

  1132.   REAL,                                    INTENT(IN   ) :: rdx,    &
    1133. 

  1133.                                                             rdy,    &
    1134. 

  1134.                                                             dts,    &
    1135. 

  1135.                                                             epssm
    1136. 

  1136. 

  1137. !  Local arrays from the stack (note tile size)
    1138. 

  1138. 

  1139.   REAL, DIMENSION (its:ite, kts:kte) :: wdtn, dvdxi
    1140. 

  1140.   REAL, DIMENSION (its:ite) :: dmdt
    1141. 

  1141. 

  1142.   INTEGER :: i,j,k, i_start, i_end, j_start, j_end, k_start, k_end
    1143. 

  1143.   INTEGER :: i_endu, j_endv
    1144. 

  1144.   REAL    :: acc
    1145. 

  1145. 

  1146. !<DESCRIPTION>
    1147. 

  1147. !
    1148. 

  1148. !  advance_mu_t advances the explicit perturbation theta equation and the mass
    1149. 

  1149. !  conservation equation.  In addition, the small timestep omega is updated,
    1150. 

  1150. !  and some quantities needed in other places are squirrelled away.
    1151. 

  1151. !
    1152. 

  1152. !</DESCRIPTION>
    1153. 

  1153. 

  1154. !  now, the real work.
    1155. 

  1155. !  set the loop bounds taking into account boundary conditions.
    1156. 

  1156. 

  1157.   i_start = its
    1158. 

  1158.   i_end   = min(ite,ide-1)
    1159. 

  1159.   j_start = jts
    1160. 

  1160.   j_end   = min(jte,jde-1)
    1161. 

  1161.   k_start = kts
    1162. 

  1162.   k_end   = kte-1
    1163. 

  1163.   IF ( .NOT. config_flags%periodic_x )THEN
    1164. 

  1164.      IF ( config_flags%specified .or. config_flags%nested ) then
    1165. 

  1165.        i_start = max(its,ids+1)
    1166. 

  1166.        i_end   = min(ite,ide-2)
    1167. 

  1167.      ENDIF
    1168. 

  1168.   ENDIF
    1169. 

  1169.   IF ( config_flags%specified .or. config_flags%nested ) then
    1170. 

  1170.      j_start = max(jts,jds+1)
    1171. 

  1171.      j_end   = min(jte,jde-2)
    1172. 

  1172.   ENDIF
    1173. 

  1173. 

  1174.   i_endu = ite
    1175. 

  1175.   j_endv = jte
    1176. 

  1176. 

  1177. !        CALCULATION OF WW (dETA/dt)
    1178. 

  1178.    DO j = j_start, j_end
    1179. 

  1179. 

  1180.      DO i=i_start, i_end
    1181. 

  1181.             dmdt(i) = 0.
    1182. 

  1182.      ENDDO
    1183. 

  1183. !  NOTE:  mu is not coupled with the map scale factor.
    1184. 

  1184. !         ww (omega) IS coupled with the map scale factor.
    1185. 

  1185. !         Being coupled with the map scale factor means 
    1186. 

  1186. !           multiplication by (1/msft) in this case.
    1187. 

  1187. 

  1188. !  Comments on map scale factors
    1189. 

  1189. !  ADT eqn 47: 
    1190. 

  1190. !  partial drho/dt = -mx*my[partial d/dx(rho u/my) + partial d/dy(rho v/mx)]
    1191. 

  1191. !                    -partial d/dz(rho w)
    1192. 

  1192. !  with rho -> mu, dividing by my, and with partial d/dnu(rho nu/my [=ww])
    1193. 

  1193. !  as the final term (because we're looking for d_nu_/dt)
    1194. 

  1194. !
    1195. 

  1195. !  begin by integrating with respect to nu from bottom to top
    1196. 

  1196. !  BCs are ww=0 at both
    1197. 

  1197. !  final term gives 0
    1198. 

  1198. !  first term gives Integral([1/my]partial d mu/dt) over total column = dm/dt
    1199. 

  1199. !  RHS remaining is Integral(-mx[partial d/dx(mu u/my) + 
    1200. 

  1200. !                                partial d/dy(mu v/mx)]) over column
    1201. 

  1201. !  lines below find RHS terms at each level then set dmdt = sum over all levels
    1202. 

  1202. !
    1203. 

  1203. !  [don't divide the below by msfty until find ww, since dmdt is used in
    1204. 

  1204. !   the meantime]
    1205. 

  1205. 

  1206.      DO k=k_start, k_end
    1207. 

  1207.      DO i=i_start, i_end
    1208. 

  1208.          dvdxi(i,k) = msftx(i,j)*msfty(i,j)*(                                      &
    1209. 

  1209.                      rdy*( (v(i,k,j+1)+muv(i,j+1)*v_1(i,k,j+1)*msfvx_inv(i,j+1))   &
    1210. 

  1210.                           -(v(i,k,j  )+muv(i,j  )*v_1(i,k,j  )*msfvx_inv(i,j  )) ) &
    1211. 

  1211.                     +rdx*( (u(i+1,k,j)+muu(i+1,j)*u_1(i+1,k,j)/msfuy(i+1,j))       &
    1212. 

  1212.                           -(u(i,k,j  )+muu(i  ,j)*u_1(i,k,j  )/msfuy(i  ,j)) ))
    1213. 

  1213.         dmdt(i)    = dmdt(i) + dnw(k)*dvdxi(i,k)
    1214. 

  1214.      ENDDO
    1215. 

  1215.      ENDDO
    1216. 

  1216.      DO i=i_start, i_end
    1217. 

  1217.        muave(i,j) = mu(i,j)
    1218. 

  1218.        mu(i,j) = mu(i,j)+dts*(dmdt(i)+mu_tend(i,j))
    1219. 

  1219.        mudf(i,j) = (dmdt(i)+mu_tend(i,j)) ! save tendency for div damp filter
    1220. 

  1220.        muts(i,j) = mut(i,j)+mu(i,j)
    1221. 

  1221.        muave(i,j) =.5*((1.+epssm)*mu(i,j)+(1.-epssm)*muave(i,j))
    1222. 

  1222.      ENDDO
    1223. 

  1223. 

  1224.      DO k=2,k_end
    1225. 

  1225.      DO i=i_start, i_end
    1226. 

  1226.        ww(i,k,j)=ww(i,k-1,j)-dnw(k-1)*(dmdt(i)+dvdxi(i,k-1)+mu_tend(i,j))/msfty(i,j)
    1227. 

  1227.      ENDDO
    1228. 

  1228.      END DO
    1229. 

  1229. 

  1230. !  NOTE:  ww_1 (large timestep ww) is already coupled with the 
    1231. 

  1231. !         map scale factor
    1232. 

  1232. 

  1233.      DO k=1,k_end
    1234. 

  1234.      DO i=i_start, i_end
    1235. 

  1235.        ww(i,k,j)=ww(i,k,j)-ww_1(i,k,j)
    1236. 

  1236.      END DO
    1237. 

  1237.      END DO
    1238. 

  1238. 

  1239.    ENDDO
    1240. 

  1240. 

  1241. ! CALCULATION OF THETA
    1242. 

  1242. 

  1243. ! NOTE: theta'' is not coupled with the map-scale factor, 
    1244. 

  1244. !       while the theta'' tendency is coupled (i.e., mult by 1/msft)
    1245. 

  1245. 

  1246. ! Comments on map scale factors
    1247. 

  1247. ! BUT NOTE THAT both are mass coupled
    1248. 

  1248. ! in flux form equations (Klemp et al.) Theta = mu*theta
    1249. 

  1249. !
    1250. 

  1250. ! scalar eqn: partial d/dt(rho q/my) = -mx[partial d/dx(q rho u/my) + 
    1251. 

  1251. !                                          partial d/dy(q rho v/mx)]
    1252. 

  1252. !                                      - partial d/dz(q rho w/my)
    1253. 

  1253. ! with rho -> mu, and with partial d/dnu(q rho nu/my) as the final term
    1254. 

  1254. !
    1255. 

  1255. ! adding previous tendency contribution which was map scale factor coupled
    1256. 

  1256. ! (had been divided by msfty)
    1257. 

  1257. ! need to uncouple before updating uncoupled Theta (by adding)
    1258. 

  1258. 

  1259.    DO j=j_start, j_end
    1260. 

  1260.      DO k=1,k_end
    1261. 

  1261.      DO i=i_start, i_end
    1262. 

  1262.        t_ave(i,k,j) = t(i,k,j)
    1263. 

  1263.        t   (i,k,j) = t(i,k,j) + msfty(i,j)*dts*ft(i,k,j)
    1264. 

  1264.      END DO
    1265. 

  1265.      END DO
    1266. 

  1266.    ENDDO   
    1267. 

  1267. 

  1268.    DO j=j_start, j_end
    1269. 

  1269. 

  1270.      DO i=i_start, i_end
    1271. 

  1271.        wdtn(i,1  )=0.
    1272. 

  1272.        wdtn(i,kde)=0.
    1273. 

  1273.      ENDDO
    1274. 

  1274. 

  1275.      DO k=2,k_end
    1276. 

  1276.      DO i=i_start, i_end
    1277. 

  1277.         ! for scalar eqn RHS term 3
    1278. 

  1278.         wdtn(i,k)= ww(i,k,j)*(fnm(k)*t_1(i,k  ,j)+fnp(k)*t_1(i,k-1,j))
    1279. 

  1279.      ENDDO
    1280. 

  1280.      ENDDO
    1281. 

  1281. 

  1282. ! scalar eqn, RHS terms 1, 2 and 3
    1283. 

  1283. ! multiply by msfty to uncouple result for Theta from map scale factor
    1284. 

  1284. 

  1285.      DO k=1,k_end
    1286. 

  1286.      DO i=i_start, i_end
    1287. 

  1287.        ! multiplication by msfty uncouples result for Theta
    1288. 

  1288.        t(i,k,j) = t(i,k,j) - dts*msfty(i,j)*(              &
    1289. 

  1289.                           ! multiplication by mx needed for RHS terms 1 & 2
    1290. 

  1290.                           msftx(i,j)*(                     &
    1291. 

  1291.                .5*rdy*                                     &
    1292. 

  1292.               ( v(i,k,j+1)*(t_1(i,k,j+1)+t_1(i,k, j ))     &
    1293. 

  1293.                -v(i,k,j  )*(t_1(i,k, j )+t_1(i,k,j-1)) )   &
    1294. 

  1294.              + .5*rdx*                                     &
    1295. 

  1295.               ( u(i+1,k,j)*(t_1(i+1,k,j)+t_1(i  ,k,j))     &
    1296. 

  1296.                -u(i  ,k,j)*(t_1(i  ,k,j)+t_1(i-1,k,j)) ) ) &
    1297. 

  1297.              + rdnw(k)*( wdtn(i,k+1)-wdtn(i,k) ) )       
    1298. 

  1298.      ENDDO
    1299. 

  1299.      ENDDO
    1300. 

  1300. 

  1301.    ENDDO
    1302. 

  1302. 

  1303. END SUBROUTINE advance_mu_t
    1304. 

  1304.           
    1305. 

  1305. 

  1306. 

  1307. !------------------------------------------------------------
    1308. 

  1308. 

  1309. SUBROUTINE advance_w( w, rw_tend, ww, w_save, u, v,  &
    1310. 

  1310.                       mu1, mut, muave, muts,      &
    1311. 

  1311.                       t_2ave, t_2, t_1,           &
    1312. 

  1312.                       ph, ph_1, phb, ph_tend,     &
    1313. 

  1313.                       ht, c2a, cqw, alt, alb,     &
    1314. 

  1314.                       a, alpha, gamma,            &
    1315. 

  1315.                       rdx, rdy, dts, t0, epssm,   &
    1316. 

  1316.                       dnw, fnm, fnp, rdnw, rdn,   &
    1317. 

  1317.                       cf1, cf2, cf3, msftx, msfty,&
    1318. 

  1318.                       config_flags, top_lid,      &
    1319. 

  1319.                       ids,ide, jds,jde, kds,kde,  & ! domain dims
    1320. 

  1320.                       ims,ime, jms,jme, kms,kme,  & ! memory dims
    1321. 

  1321.                       its,ite, jts,jte, kts,kte  )  ! tile   dims
    1322. 

  1322. 

  1323. ! We have used msfty for msft_inv but have not thought through w equation
    1324. 

  1324. ! pieces properly yet, so we will have to hope that it is okay
    1325. 

  1325. ! We think we have found a slight error in surface w calculation
    1326. 

  1326. 

  1327.   IMPLICIT NONE ! religion first
    1328. 

  1328.       
    1329. 

  1329. ! stuff coming in
    1330. 

  1330. 

  1331. 

  1332.   TYPE(grid_config_rec_type), INTENT(IN   ) :: config_flags
    1333. 

  1333. 

  1334.   INTEGER,      INTENT(IN   )    :: ids,ide, jds,jde, kds,kde
    1335. 

  1335.   INTEGER,      INTENT(IN   )    :: ims,ime, jms,jme, kms,kme
    1336. 

  1336.   INTEGER,      INTENT(IN   )    :: its,ite, jts,jte, kts,kte
    1337. 

  1337. 

  1338.   LOGICAL,      INTENT(IN   )    :: top_lid
    1339. 

  1339. 

  1340.       REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
    1341. 

  1341.             INTENT(INOUT) ::                          &
    1342. 

  1342.                                              t_2ave,  &
    1343. 

  1343.                                              w,       &
    1344. 

  1344.                                              ph
    1345. 

  1345. 

  1346.       REAL, DIMENSION(  ims:ime , kms:kme, jms:jme ), &
    1347. 

  1347.             INTENT(IN   ) ::                          &
    1348. 

  1348.                                              rw_tend, &
    1349. 

  1349.                                              ww,      &
    1350. 

  1350.                                              w_save,  &
    1351. 

  1351.                                              u,       &
    1352. 

  1352.                                              v,       &
    1353. 

  1353.                                              t_2,     &
    1354. 

  1354.                                              t_1,     &
    1355. 

  1355.                                              ph_1,    &
    1356. 

  1356.                                              phb,     &
    1357. 

  1357.                                              ph_tend, &
    1358. 

  1358.                                              alpha,   &
    1359. 

  1359.                                              gamma,   &
    1360. 

  1360.                                              a,       &
    1361. 

  1361.                                              c2a,     &
    1362. 

  1362.                                              cqw,     &
    1363. 

  1363.                                              alb,     &
    1364. 

  1364.                                              alt
    1365. 

  1365. 

  1366.       REAL, DIMENSION( ims:ime , jms:jme ), &
    1367. 

  1367.             INTENT(IN   )  ::               &
    1368. 

  1368.                                    mu1,     &
    1369. 

  1369.                                    mut,     &
    1370. 

  1370.                                    muave,   &
    1371. 

  1371.                                    muts,    &
    1372. 

  1372.                                    ht,      &
    1373. 

  1373.                                    msftx,   &
    1374. 

  1374.                                    msfty
    1375. 

  1375. 

  1376.       REAL, DIMENSION( kms:kme ),  INTENT(IN   )  :: fnp,     &
    1377. 

  1377.                                                      fnm,     &
    1378. 

  1378.                                                      rdnw,    &
    1379. 

  1379.                                                      rdn,     &
    1380. 

  1380.                                                      dnw
    1381. 

  1381. 

  1382.       REAL,   INTENT(IN   )  :: rdx,     &
    1383. 

  1383.                                 rdy,     &
    1384. 

  1384.                                 dts,     &
    1385. 

  1385.                                 cf1,     &
    1386. 

  1386.                                 cf2,     &
    1387. 

  1387.                                 cf3,     &
    1388. 

  1388.                                 t0,      &
    1389. 

  1389.                                 epssm
    1390. 

  1390. 

  1391. !  Stack based 3d data, tile size.
    1392. 

  1392. 

  1393.       REAL, DIMENSION( its:ite ) :: mut_inv, msft_inv
    1394. 

  1394.       REAL, DIMENSION( its:ite, kts:kte ) :: rhs, wdwn
    1395. 

  1395.       INTEGER :: i,j,k, i_start, i_end, j_start, j_end, k_start, k_end
    1396. 

  1396.       REAL, DIMENSION (kts:kte) :: dampwt
    1397. 

  1397.       real :: htop,hbot,hdepth,hk
    1398. 

  1398.       real    :: pi,dampmag
    1399. 

  1399. 

  1400. !<DESCRIPTION>
    1401. 

  1401. !
    1402. 

  1402. !  advance_w advances the implicit w and geopotential equations.
    1403. 

  1403. !
    1404. 

  1404. !</DESCRIPTION>
    1405. 

  1405. 

  1406. !  set loop limits.
    1407. 

  1407. !  Currently set for periodic boundary conditions
    1408. 

  1408. 

  1409.       i_start = its
    1410. 

  1410.       i_end   = min(ite,ide-1)
    1411. 

  1411.       j_start = jts
    1412. 

  1412.       j_end   = min(jte,jde-1)
    1413. 

  1413.       k_start = kts
    1414. 

  1414.       k_end   = kte-1
    1415. 

  1415.       IF ( .NOT. config_flags%periodic_x )THEN
    1416. 

  1416.          IF ( config_flags%specified .or. config_flags%nested ) then
    1417. 

  1417.            i_start = max(its,ids+1)
    1418. 

  1418.            i_end   = min(ite,ide-2)
    1419. 

  1419.          ENDIF
    1420. 

  1420.       ENDIF
    1421. 

  1421.       IF ( config_flags%specified .or. config_flags%nested ) then
    1422. 

  1422.          j_start = max(jts,jds+1)
    1423. 

  1423.          j_end   = min(jte,jde-2)
    1424. 

  1424.       ENDIF
    1425. 

  1425. 

  1426.    pi = 4.*atan(1.)
    1427. 

  1427.    dampmag = dts*config_flags%dampcoef
    1428. 

  1428.    hdepth=config_flags%zdamp
    1429. 

  1429. 

  1430. ! calculation of phi and w equations
    1431. 

  1431. 

  1432. !   Comments on map scale factors:
    1433. 

  1433. !   phi equation is:
    1434. 

  1434. !    partial d/dt(rho phi/my) = -mx partial d/dx(phi rho u/my) 
    1435. 

  1435. !                               -mx partial d/dy(phi rho v/mx)
    1436. 

  1436. !                               - partial d/dz(phi rho w/my) + rho g w/my
    1437. 

  1437. !   as with scalar equation, use uncoupled value (here phi) to find the
    1438. 

  1438. !   coupled tendency (rho phi/my)
    1439. 

  1439. !   here as usual rho -> ~'mu'
    1440. 

  1440. !
    1441. 

  1441. !   w eqn [divided by my] is:
    1442. 

  1442. !     partial d/dt(rho w/my) = -mx partial d/dx(w rho u/my)
    1443. 

  1443. !                              -mx partial d/dy(v rho v/mx)
    1444. 

  1444. !                              - partial d/dz(w rho w/my)
    1445. 

  1445. !                              +rho[(u*u+v*v)/a + 2 u omega cos(lat) -
    1446. 

  1446. !                                   (1/rho) partial dp/dz - g + Fz]/my
    1447. 

  1447. !   here as usual rho -> ~'mu'
    1448. 

  1448. !
    1449. 

  1449. !  'u,v,w' sent here must be coupled variables (= rho w/my etc.) by their usage
    1450. 

  1450. 

  1451. 

  1452.     DO i=i_start, i_end
    1453. 

  1453.       rhs(i,1) = 0.
    1454. 

  1454.     ENDDO
    1455. 

  1455. 

  1456.   j_loop_w:  DO j = j_start, j_end
    1457. 

  1457.     DO i=i_start, i_end
    1458. 

  1458.        mut_inv(i) = 1./mut(i,j)
    1459. 

  1459.        msft_inv(i) = 1./msfty(i,j)
    1460. 

  1460.     ENDDO
    1461. 

  1461. 

  1462.     DO k=1, k_end
    1463. 

  1463.     DO i=i_start, i_end
    1464. 

  1464.       t_2ave(i,k,j)=.5*((1.+epssm)*t_2(i,k,j)       &
    1465. 

  1465.                     +(1.-epssm)*t_2ave(i,k,j))
    1466. 

  1466.       t_2ave(i,k,j)=(t_2ave(i,k,j) + muave(i,j)*t0) &
    1467. 

  1467.                     /(muts(i,j)*(t0+t_1(i,k,j)))
    1468. 

  1468.       wdwn(i,k+1)=.5*(ww(i,k+1,j)+ww(i,k,j))*rdnw(k)    &
    1469. 

  1469.            *(ph_1(i,k+1,j)-ph_1(i,k,j)+phb(i,k+1,j)-phb(i,k,j))
    1470. 

  1470.       rhs(i,k+1) = dts*(ph_tend(i,k+1,j) + .5*g*(1.-epssm)*w(i,k+1,j))
    1471. 

  1471. 

  1472.     ENDDO
    1473. 

  1473.     ENDDO
    1474. 

  1474. 

  1475. ! …
    1476.

Large files files are truncated, but you can click here to view the full file