PageRenderTime 67ms CodeModel.GetById 24ms RepoModel.GetById 1ms app.codeStats 1ms

/wrfv2_fire/dyn_em/module_big_step_utilities_em.F

https://github.com/jbeezley/wrf-fire
FORTRAN Legacy | 6155 lines | 3822 code | 1351 blank | 982 comment | 56 complexity | 28eae140e72dd9b7b229817c61222a96 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. #if (RWORDSIZE == 4)
    5. 

  5. #   define VPOWX vspowx
    6. 

  6. #   define VPOW  vspow
    7. 

  7. #else
    8. 

  8. #   define VPOWX vpowx
    9. 

  9. #   define VPOW  vpow
    10. 

  10. #endif
    11. 

  11. 

  12. 

  13. MODULE module_big_step_utilities_em
    14. 

  14. 

  15.    USE module_model_constants
    16. 

  16.    USE module_state_description, only: p_qg, p_qs, p_qi, gdscheme, tiedtkescheme, kfetascheme, g3scheme, &
    17. 

  17.    p_qv, param_first_scalar, p_qr, p_qc
    18. 

  18.    USE module_configure, ONLY : grid_config_rec_type
    19. 

  19. 

  20. CONTAINS
    21. 

  21. 

  22. !-------------------------------------------------------------------------------
    23. 

  23. 

  24. SUBROUTINE calc_mu_uv ( config_flags,                 &
    25. 

  25.                         mu, mub, muu, muv,            &
    26. 

  26.                         ids, ide, jds, jde, kds, kde, &
    27. 

  27.                         ims, ime, jms, jme, kms, kme, &
    28. 

  28.                         its, ite, jts, jte, kts, kte )
    29. 

  29. 

  30.    IMPLICIT NONE
    31. 

  31.    
    32. 

  32.    ! Input data
    33. 

  33. 

  34.    TYPE(grid_config_rec_type   ) ,   INTENT(IN   ) :: config_flags
    35. 

  35. 

  36.    INTEGER ,          INTENT(IN   ) :: ids, ide, jds, jde, kds, kde, &
    37. 

  37.                                        ims, ime, jms, jme, kms, kme, &
    38. 

  38.                                        its, ite, jts, jte, kts, kte 
    39. 

  39. 

  40.    REAL, DIMENSION( ims:ime , jms:jme ) , INTENT(  OUT) :: muu, muv
    41. 

  41.    REAL, DIMENSION( ims:ime , jms:jme ) , INTENT(IN   ) :: mu, mub
    42. 

  42. 

  43.    !  local stuff
    44. 

  44. 

  45.    INTEGER :: i, j, itf, jtf, im, jm
    46. 

  46. 

  47. !<DESCRIPTION>
    48. 

  48. !
    49. 

  49. !  calc_mu_uv calculates the full column dry-air mass at the staggered
    50. 

  50. !  horizontal velocity points (u,v) and places the results in muu and muv.
    51. 

  51. !  This routine uses the reference state (mub) and perturbation state (mu)
    52. 

  52. !
    53. 

  53. !</DESCRIPTION>
    54. 

  54. 

  55. 

  56.       itf=ite
    57. 

  57.       jtf=MIN(jte,jde-1)
    58. 

  58. 

  59.       IF      ( ( its .NE. ids ) .AND. ( ite .NE. ide ) ) THEN
    60. 

  60.          DO j=jts,jtf
    61. 

  61.          DO i=its,itf
    62. 

  62.             muu(i,j) = 0.5*(mu(i,j)+mu(i-1,j)+mub(i,j)+mub(i-1,j))
    63. 

  63.          ENDDO
    64. 

  64.          ENDDO
    65. 

  65.       ELSE IF ( ( its .EQ. ids ) .AND. ( ite .NE. ide ) ) THEN
    66. 

  66.          DO j=jts,jtf
    67. 

  67.          DO i=its+1,itf
    68. 

  68.             muu(i,j) = 0.5*(mu(i,j)+mu(i-1,j)+mub(i,j)+mub(i-1,j))
    69. 

  69.          ENDDO
    70. 

  70.          ENDDO
    71. 

  71.          i=its
    72. 

  72.          im = its
    73. 

  73.          if(config_flags%periodic_x) im = its-1
    74. 

  74.          DO j=jts,jtf
    75. 

  75. !            muu(i,j) =      mu(i,j)          +mub(i,j)
    76. 

  76. !  fix for periodic b.c., 13 march 2004, wcs
    77. 

  77.             muu(i,j) = 0.5*(mu(i,j)+mu(im,j)+mub(i,j)+mub(im,j))
    78. 

  78.          ENDDO
    79. 

  79.       ELSE IF ( ( its .NE. ids ) .AND. ( ite .EQ. ide ) ) THEN
    80. 

  80.          DO j=jts,jtf
    81. 

  81.          DO i=its,itf-1
    82. 

  82.             muu(i,j) = 0.5*(mu(i,j)+mu(i-1,j)+mub(i,j)+mub(i-1,j))
    83. 

  83.          ENDDO
    84. 

  84.          ENDDO
    85. 

  85.          i=ite
    86. 

  86.          im = ite-1
    87. 

  87.          if(config_flags%periodic_x) im = ite
    88. 

  88.          DO j=jts,jtf
    89. 

  89. !            muu(i,j) =      mu(i-1,j)        +mub(i-1,j)
    90. 

  90. !  fix for periodic b.c., 13 march 2004, wcs
    91. 

  91.             muu(i,j) = 0.5*(mu(i-1,j)+mu(im,j)+mub(i-1,j)+mub(im,j))
    92. 

  92.          ENDDO
    93. 

  93.       ELSE IF ( ( its .EQ. ids ) .AND. ( ite .EQ. ide ) ) THEN
    94. 

  94.          DO j=jts,jtf
    95. 

  95.          DO i=its+1,itf-1
    96. 

  96.             muu(i,j) = 0.5*(mu(i,j)+mu(i-1,j)+mub(i,j)+mub(i-1,j))
    97. 

  97.          ENDDO
    98. 

  98.          ENDDO
    99. 

  99.          i=its
    100. 

  100.          im = its
    101. 

  101.          if(config_flags%periodic_x) im = its-1
    102. 

  102.          DO j=jts,jtf
    103. 

  103. !            muu(i,j) =      mu(i,j)          +mub(i,j)
    104. 

  104. !  fix for periodic b.c., 13 march 2004, wcs
    105. 

  105.             muu(i,j) = 0.5*(mu(i,j)+mu(im,j)+mub(i,j)+mub(im,j))
    106. 

  106.          ENDDO
    107. 

  107.          i=ite
    108. 

  108.          im = ite-1
    109. 

  109.          if(config_flags%periodic_x) im = ite
    110. 

  110.          DO j=jts,jtf
    111. 

  111. !            muu(i,j) =      mu(i-1,j)        +mub(i-1,j)
    112. 

  112. !  fix for periodic b.c., 13 march 2004, wcs
    113. 

  113.             muu(i,j) = 0.5*(mu(i-1,j)+mu(im,j)+mub(i-1,j)+mub(im,j))
    114. 

  114.          ENDDO
    115. 

  115.       END IF
    116. 

  116. 

  117.       itf=MIN(ite,ide-1)
    118. 

  118.       jtf=jte
    119. 

  119. 

  120.       IF      ( ( jts .NE. jds ) .AND. ( jte .NE. jde ) ) THEN
    121. 

  121.          DO j=jts,jtf
    122. 

  122.          DO i=its,itf
    123. 

  123.              muv(i,j) = 0.5*(mu(i,j)+mu(i,j-1)+mub(i,j)+mub(i,j-1))
    124. 

  124.          ENDDO
    125. 

  125.          ENDDO
    126. 

  126.       ELSE IF ( ( jts .EQ. jds ) .AND. ( jte .NE. jde ) ) THEN
    127. 

  127.          DO j=jts+1,jtf
    128. 

  128.          DO i=its,itf
    129. 

  129.              muv(i,j) = 0.5*(mu(i,j)+mu(i,j-1)+mub(i,j)+mub(i,j-1))
    130. 

  130.          ENDDO
    131. 

  131.          ENDDO
    132. 

  132.          j=jts
    133. 

  133.          jm = jts
    134. 

  134.          if(config_flags%periodic_y) jm = jts-1
    135. 

  135.          DO i=its,itf
    136. 

  136. !             muv(i,j) =      mu(i,j)          +mub(i,j)
    137. 

  137. !  fix for periodic b.c., 13 march 2004, wcs
    138. 

  138.              muv(i,j) = 0.5*(mu(i,j)+mu(i,jm)+mub(i,j)+mub(i,jm))
    139. 

  139.          ENDDO
    140. 

  140.       ELSE IF ( ( jts .NE. jds ) .AND. ( jte .EQ. jde ) ) THEN
    141. 

  141.          DO j=jts,jtf-1
    142. 

  142.          DO i=its,itf
    143. 

  143.              muv(i,j) = 0.5*(mu(i,j)+mu(i,j-1)+mub(i,j)+mub(i,j-1))
    144. 

  144.          ENDDO
    145. 

  145.          ENDDO
    146. 

  146.          j=jte
    147. 

  147.          jm = jte-1
    148. 

  148.          if(config_flags%periodic_y) jm = jte
    149. 

  149.          DO i=its,itf
    150. 

  150.              muv(i,j) =      mu(i,j-1)        +mub(i,j-1)
    151. 

  151. !  fix for periodic b.c., 13 march 2004, wcs
    152. 

  152.              muv(i,j) = 0.5*(mu(i,j-1)+mu(i,jm)+mub(i,j-1)+mub(i,jm))
    153. 

  153.          ENDDO
    154. 

  154.       ELSE IF ( ( jts .EQ. jds ) .AND. ( jte .EQ. jde ) ) THEN
    155. 

  155.          DO j=jts+1,jtf-1
    156. 

  156.          DO i=its,itf
    157. 

  157.              muv(i,j) = 0.5*(mu(i,j)+mu(i,j-1)+mub(i,j)+mub(i,j-1))
    158. 

  158.          ENDDO
    159. 

  159.          ENDDO
    160. 

  160.          j=jts
    161. 

  161.          jm = jts
    162. 

  162.          if(config_flags%periodic_y) jm = jts-1
    163. 

  163.          DO i=its,itf
    164. 

  164. !             muv(i,j) =      mu(i,j)          +mub(i,j)
    165. 

  165. !  fix for periodic b.c., 13 march 2004, wcs
    166. 

  166.              muv(i,j) = 0.5*(mu(i,j)+mu(i,jm)+mub(i,j)+mub(i,jm))
    167. 

  167.          ENDDO
    168. 

  168.          j=jte
    169. 

  169.          jm = jte-1
    170. 

  170.          if(config_flags%periodic_y) jm = jte
    171. 

  171.          DO i=its,itf
    172. 

  172. !             muv(i,j) =      mu(i,j-1)        +mub(i,j-1)
    173. 

  173. !  fix for periodic b.c., 13 march 2004, wcs
    174. 

  174.              muv(i,j) = 0.5*(mu(i,j-1)+mu(i,jm)+mub(i,j-1)+mub(i,jm))
    175. 

  175.          ENDDO
    176. 

  176.       END IF
    177. 

  177. 

  178. END SUBROUTINE calc_mu_uv
    179. 

  179. 

  180. !-------------------------------------------------------------------------------
    181. 

  181. 

  182. SUBROUTINE calc_mu_uv_1 ( config_flags,                 &
    183. 

  183.                           mu, muu, muv,                 &
    184. 

  184.                           ids, ide, jds, jde, kds, kde, &
    185. 

  185.                           ims, ime, jms, jme, kms, kme, &
    186. 

  186.                           its, ite, jts, jte, kts, kte )
    187. 

  187. 

  188.    IMPLICIT NONE
    189. 

  189.    
    190. 

  190.    ! Input data
    191. 

  191. 

  192.    TYPE(grid_config_rec_type   ) ,   INTENT(IN   ) :: config_flags
    193. 

  193. 

  194.    INTEGER ,          INTENT(IN   ) :: ids, ide, jds, jde, kds, kde, &
    195. 

  195.                                        ims, ime, jms, jme, kms, kme, &
    196. 

  196.                                        its, ite, jts, jte, kts, kte 
    197. 

  197. 

  198.    REAL, DIMENSION( ims:ime , jms:jme ) , INTENT(  OUT) :: muu, muv
    199. 

  199.    REAL, DIMENSION( ims:ime , jms:jme ) , INTENT(IN   ) :: mu
    200. 

  200. 

  201.    !  local stuff
    202. 

  202. 

  203.    INTEGER :: i, j, itf, jtf, im, jm
    204. 

  204. 

  205. !<DESCRIPTION>
    206. 

  206. !
    207. 

  207. !  calc_mu_uv calculates the full column dry-air mass at the staggered
    208. 

  208. !  horizontal velocity points (u,v) and places the results in muu and muv.
    209. 

  209. !  This routine uses the full state (mu)
    210. 

  210. !
    211. 

  211. !</DESCRIPTION>
    212. 

  212.    
    213. 

  213.       itf=ite
    214. 

  214.       jtf=MIN(jte,jde-1)
    215. 

  215. 

  216.       IF      ( ( its .NE. ids ) .AND. ( ite .NE. ide ) ) THEN
    217. 

  217.          DO j=jts,jtf
    218. 

  218.          DO i=its,itf
    219. 

  219.             muu(i,j) = 0.5*(mu(i,j)+mu(i-1,j))
    220. 

  220.          ENDDO
    221. 

  221.          ENDDO
    222. 

  222.       ELSE IF ( ( its .EQ. ids ) .AND. ( ite .NE. ide ) ) THEN
    223. 

  223.          DO j=jts,jtf
    224. 

  224.          DO i=its+1,itf
    225. 

  225.             muu(i,j) = 0.5*(mu(i,j)+mu(i-1,j))
    226. 

  226.          ENDDO
    227. 

  227.          ENDDO
    228. 

  228.          i=its
    229. 

  229.          im = its
    230. 

  230.          if(config_flags%periodic_x) im = its-1
    231. 

  231.          DO j=jts,jtf
    232. 

  232.             muu(i,j) = 0.5*(mu(i,j)+mu(im,j))
    233. 

  233.          ENDDO
    234. 

  234.       ELSE IF ( ( its .NE. ids ) .AND. ( ite .EQ. ide ) ) THEN
    235. 

  235.          DO j=jts,jtf
    236. 

  236.          DO i=its,itf-1
    237. 

  237.             muu(i,j) = 0.5*(mu(i,j)+mu(i-1,j))
    238. 

  238.          ENDDO
    239. 

  239.          ENDDO
    240. 

  240.          i=ite
    241. 

  241.          im = ite-1
    242. 

  242.          if(config_flags%periodic_x) im = ite
    243. 

  243.          DO j=jts,jtf
    244. 

  244.             muu(i,j) = 0.5*(mu(i-1,j)+mu(im,j))
    245. 

  245.          ENDDO
    246. 

  246.       ELSE IF ( ( its .EQ. ids ) .AND. ( ite .EQ. ide ) ) THEN
    247. 

  247.          DO j=jts,jtf
    248. 

  248.          DO i=its+1,itf-1
    249. 

  249.             muu(i,j) = 0.5*(mu(i,j)+mu(i-1,j))
    250. 

  250.          ENDDO
    251. 

  251.          ENDDO
    252. 

  252.          i=its
    253. 

  253.          im = its
    254. 

  254.          if(config_flags%periodic_x) im = its-1
    255. 

  255.          DO j=jts,jtf
    256. 

  256.             muu(i,j) = 0.5*(mu(i,j)+mu(im,j))
    257. 

  257.          ENDDO
    258. 

  258.          i=ite
    259. 

  259.          im = ite-1
    260. 

  260.          if(config_flags%periodic_x) im = ite
    261. 

  261.          DO j=jts,jtf
    262. 

  262.             muu(i,j) = 0.5*(mu(i-1,j)+mu(im,j))
    263. 

  263.          ENDDO
    264. 

  264.       END IF
    265. 

  265. 

  266.       itf=MIN(ite,ide-1)
    267. 

  267.       jtf=jte
    268. 

  268. 

  269.       IF      ( ( jts .NE. jds ) .AND. ( jte .NE. jde ) ) THEN
    270. 

  270.          DO j=jts,jtf
    271. 

  271.          DO i=its,itf
    272. 

  272.              muv(i,j) = 0.5*(mu(i,j)+mu(i,j-1))
    273. 

  273.          ENDDO
    274. 

  274.          ENDDO
    275. 

  275.       ELSE IF ( ( jts .EQ. jds ) .AND. ( jte .NE. jde ) ) THEN
    276. 

  276.          DO j=jts+1,jtf
    277. 

  277.          DO i=its,itf
    278. 

  278.              muv(i,j) = 0.5*(mu(i,j)+mu(i,j-1))
    279. 

  279.          ENDDO
    280. 

  280.          ENDDO
    281. 

  281.          j=jts
    282. 

  282.          jm = jts
    283. 

  283.          if(config_flags%periodic_y) jm = jts-1
    284. 

  284.          DO i=its,itf
    285. 

  285.              muv(i,j) = 0.5*(mu(i,j)+mu(i,jm))
    286. 

  286.          ENDDO
    287. 

  287.       ELSE IF ( ( jts .NE. jds ) .AND. ( jte .EQ. jde ) ) THEN
    288. 

  288.          DO j=jts,jtf-1
    289. 

  289.          DO i=its,itf
    290. 

  290.              muv(i,j) = 0.5*(mu(i,j)+mu(i,j-1))
    291. 

  291.          ENDDO
    292. 

  292.          ENDDO
    293. 

  293.          j=jte
    294. 

  294.          jm = jte-1
    295. 

  295.          if(config_flags%periodic_y) jm = jte
    296. 

  296.          DO i=its,itf
    297. 

  297.              muv(i,j) = 0.5*(mu(i,j-1)+mu(i,jm))
    298. 

  298.          ENDDO
    299. 

  299.       ELSE IF ( ( jts .EQ. jds ) .AND. ( jte .EQ. jde ) ) THEN
    300. 

  300.          DO j=jts+1,jtf-1
    301. 

  301.          DO i=its,itf
    302. 

  302.              muv(i,j) = 0.5*(mu(i,j)+mu(i,j-1))
    303. 

  303.          ENDDO
    304. 

  304.          ENDDO
    305. 

  305.          j=jts
    306. 

  306.          jm = jts
    307. 

  307.          if(config_flags%periodic_y) jm = jts-1
    308. 

  308.          DO i=its,itf
    309. 

  309.              muv(i,j) = 0.5*(mu(i,j)+mu(i,jm))
    310. 

  310.          ENDDO
    311. 

  311.          j=jte
    312. 

  312.          jm = jte-1
    313. 

  313.          if(config_flags%periodic_y) jm = jte
    314. 

  314.          DO i=its,itf
    315. 

  315.              muv(i,j) = 0.5*(mu(i,j-1)+mu(i,jm))
    316. 

  316.          ENDDO
    317. 

  317.       END IF
    318. 

  318. 

  319. END SUBROUTINE calc_mu_uv_1
    320. 

  320. 

  321. !-------------------------------------------------------------------------------
    322. 

  322. 

  323. ! Map scale factor comments for this routine:
    324. 

  324. ! Locally not changed, but sent the correct map scale factors
    325. 

  325. ! from module_em (msfuy, msfvx, msfty)
    326. 

  326. 

  327. SUBROUTINE couple_momentum ( muu, ru, u, msfu,              &
    328. 

  328.                              muv, rv, v, msfv, msfv_inv,    &
    329. 

  329.                              mut, rw, w, msft,              &
    330. 

  330.                              ids, ide, jds, jde, kds, kde,  &
    331. 

  331.                              ims, ime, jms, jme, kms, kme,  &
    332. 

  332.                              its, ite, jts, jte, kts, kte  )
    333. 

  333. 

  334.    IMPLICIT NONE
    335. 

  335. 

  336.    ! Input data
    337. 

  337. 

  338.    INTEGER ,             INTENT(IN   ) :: ids, ide, jds, jde, kds, kde, &
    339. 

  339.                                           ims, ime, jms, jme, kms, kme, &
    340. 

  340.                                           its, ite, jts, jte, kts, kte
    341. 

  341. 

  342.    REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(  OUT) :: ru, rv, rw
    343. 

  343. 

  344.    REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN   ) :: muu, muv, mut
    345. 

  345.    REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN   ) :: msfu, msfv, msft, msfv_inv
    346. 

  346.    
    347. 

  347.    REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN   ) :: u, v, w
    348. 

  348.    
    349. 

  349.    ! Local data
    350. 

  350.    
    351. 

  351.    INTEGER :: i, j, k, itf, jtf, ktf
    352. 

  352.    
    353. 

  353. !<DESCRIPTION>
    354. 

  354. !
    355. 

  355. ! couple_momentum couples the velocities to the full column mass and
    356. 

  356. ! the map factors.
    357. 

  357. !
    358. 

  358. !</DESCRIPTION>
    359. 

  359. 

  360.    ktf=MIN(kte,kde-1)
    361. 

  361.    
    362. 

  362.       itf=ite
    363. 

  363.       jtf=MIN(jte,jde-1)
    364. 

  364. 

  365.       DO j=jts,jtf
    366. 

  366.       DO k=kts,ktf
    367. 

  367.       DO i=its,itf
    368. 

  368.          ru(i,k,j)=u(i,k,j)*muu(i,j)/msfu(i,j)
    369. 

  369.       ENDDO
    370. 

  370.       ENDDO
    371. 

  371.       ENDDO
    372. 

  372. 

  373.       itf=MIN(ite,ide-1)
    374. 

  374.       jtf=jte
    375. 

  375. 

  376.       DO j=jts,jtf
    377. 

  377.       DO k=kts,ktf
    378. 

  378.       DO i=its,itf
    379. 

  379.            rv(i,k,j)=v(i,k,j)*muv(i,j)*msfv_inv(i,j)
    380. 

  380. !           rv(i,k,j)=v(i,k,j)*muv(i,j)/msfv(i,j)
    381. 

  381.       ENDDO
    382. 

  382.       ENDDO
    383. 

  383.       ENDDO
    384. 

  384. 

  385.       itf=MIN(ite,ide-1)
    386. 

  386.       jtf=MIN(jte,jde-1)
    387. 

  387. 

  388.       DO j=jts,jtf
    389. 

  389.       DO k=kts,kte
    390. 

  390.       DO i=its,itf
    391. 

  391.          rw(i,k,j)=w(i,k,j)*mut(i,j)/msft(i,j)
    392. 

  392.       ENDDO
    393. 

  393.       ENDDO
    394. 

  394.       ENDDO
    395. 

  395. 

  396. END SUBROUTINE couple_momentum
    397. 

  397. 

  398. !-------------------------------------------------------------------
    399. 

  399. 

  400. SUBROUTINE calc_mu_staggered ( mu, mub, muu, muv,            &
    401. 

  401.                                   ids, ide, jds, jde, kds, kde, &
    402. 

  402.                                   ims, ime, jms, jme, kms, kme, &
    403. 

  403.                                   its, ite, jts, jte, kts, kte )
    404. 

  404. 

  405.    IMPLICIT NONE
    406. 

  406.    
    407. 

  407.    ! Input data
    408. 

  408. 

  409.    INTEGER ,          INTENT(IN   ) :: ids, ide, jds, jde, kds, kde, &
    410. 

  410.                                        ims, ime, jms, jme, kms, kme, &
    411. 

  411.                                        its, ite, jts, jte, kts, kte 
    412. 

  412. 

  413.    REAL, DIMENSION( ims:ime , jms:jme ) , INTENT(  OUT) :: muu, muv
    414. 

  414.    REAL, DIMENSION( ims:ime , jms:jme ) , INTENT(IN   ) :: mu, mub
    415. 

  415. 

  416.    !  local stuff
    417. 

  417. 

  418.    INTEGER :: i, j, itf, jtf
    419. 

  419. 

  420. !<DESCRIPTION>
    421. 

  421. !
    422. 

  422. ! calc_mu_staggered calculates the full dry air mass at the staggered
    423. 

  423. ! velocity points (u,v).
    424. 

  424. !
    425. 

  425. !</DESCRIPTION>
    426. 

  426.    
    427. 

  427.       itf=ite
    428. 

  428.       jtf=MIN(jte,jde-1)
    429. 

  429. 

  430.       IF      ( ( its .NE. ids ) .AND. ( ite .NE. ide ) ) THEN
    431. 

  431.          DO j=jts,jtf
    432. 

  432.          DO i=its,itf
    433. 

  433.             muu(i,j) = 0.5*(mu(i,j)+mu(i-1,j)+mub(i,j)+mub(i-1,j))
    434. 

  434.          ENDDO
    435. 

  435.          ENDDO
    436. 

  436.       ELSE IF ( ( its .EQ. ids ) .AND. ( ite .NE. ide ) ) THEN
    437. 

  437.          DO j=jts,jtf
    438. 

  438.          DO i=its+1,itf
    439. 

  439.             muu(i,j) = 0.5*(mu(i,j)+mu(i-1,j)+mub(i,j)+mub(i-1,j))
    440. 

  440.          ENDDO
    441. 

  441.          ENDDO
    442. 

  442.          i=its
    443. 

  443.          DO j=jts,jtf
    444. 

  444.             muu(i,j) =      mu(i,j)          +mub(i,j)
    445. 

  445.          ENDDO
    446. 

  446.       ELSE IF ( ( its .NE. ids ) .AND. ( ite .EQ. ide ) ) THEN
    447. 

  447.          DO j=jts,jtf
    448. 

  448.          DO i=its,itf-1
    449. 

  449.             muu(i,j) = 0.5*(mu(i,j)+mu(i-1,j)+mub(i,j)+mub(i-1,j))
    450. 

  450.          ENDDO
    451. 

  451.          ENDDO
    452. 

  452.          i=ite
    453. 

  453.          DO j=jts,jtf
    454. 

  454.             muu(i,j) =      mu(i-1,j)        +mub(i-1,j)
    455. 

  455.          ENDDO
    456. 

  456.       ELSE IF ( ( its .EQ. ids ) .AND. ( ite .EQ. ide ) ) THEN
    457. 

  457.          DO j=jts,jtf
    458. 

  458.          DO i=its+1,itf-1
    459. 

  459.             muu(i,j) = 0.5*(mu(i,j)+mu(i-1,j)+mub(i,j)+mub(i-1,j))
    460. 

  460.          ENDDO
    461. 

  461.          ENDDO
    462. 

  462.          i=its
    463. 

  463.          DO j=jts,jtf
    464. 

  464.             muu(i,j) =      mu(i,j)          +mub(i,j)
    465. 

  465.          ENDDO
    466. 

  466.          i=ite
    467. 

  467.          DO j=jts,jtf
    468. 

  468.             muu(i,j) =      mu(i-1,j)        +mub(i-1,j)
    469. 

  469.          ENDDO
    470. 

  470.       END IF
    471. 

  471. 

  472.       itf=MIN(ite,ide-1)
    473. 

  473.       jtf=jte
    474. 

  474. 

  475.       IF      ( ( jts .NE. jds ) .AND. ( jte .NE. jde ) ) THEN
    476. 

  476.          DO j=jts,jtf
    477. 

  477.          DO i=its,itf
    478. 

  478.              muv(i,j) = 0.5*(mu(i,j)+mu(i,j-1)+mub(i,j)+mub(i,j-1))
    479. 

  479.          ENDDO
    480. 

  480.          ENDDO
    481. 

  481.       ELSE IF ( ( jts .EQ. jds ) .AND. ( jte .NE. jde ) ) THEN
    482. 

  482.          DO j=jts+1,jtf
    483. 

  483.          DO i=its,itf
    484. 

  484.              muv(i,j) = 0.5*(mu(i,j)+mu(i,j-1)+mub(i,j)+mub(i,j-1))
    485. 

  485.          ENDDO
    486. 

  486.          ENDDO
    487. 

  487.          j=jts
    488. 

  488.          DO i=its,itf
    489. 

  489.              muv(i,j) =      mu(i,j)          +mub(i,j)
    490. 

  490.          ENDDO
    491. 

  491.       ELSE IF ( ( jts .NE. jds ) .AND. ( jte .EQ. jde ) ) THEN
    492. 

  492.          DO j=jts,jtf-1
    493. 

  493.          DO i=its,itf
    494. 

  494.              muv(i,j) = 0.5*(mu(i,j)+mu(i,j-1)+mub(i,j)+mub(i,j-1))
    495. 

  495.          ENDDO
    496. 

  496.          ENDDO
    497. 

  497.          j=jte
    498. 

  498.          DO i=its,itf
    499. 

  499.              muv(i,j) =      mu(i,j-1)        +mub(i,j-1)
    500. 

  500.          ENDDO
    501. 

  501.       ELSE IF ( ( jts .EQ. jds ) .AND. ( jte .EQ. jde ) ) THEN
    502. 

  502.          DO j=jts+1,jtf-1
    503. 

  503.          DO i=its,itf
    504. 

  504.              muv(i,j) = 0.5*(mu(i,j)+mu(i,j-1)+mub(i,j)+mub(i,j-1))
    505. 

  505.          ENDDO
    506. 

  506.          ENDDO
    507. 

  507.          j=jts
    508. 

  508.          DO i=its,itf
    509. 

  509.              muv(i,j) =      mu(i,j)          +mub(i,j)
    510. 

  510.          ENDDO
    511. 

  511.          j=jte
    512. 

  512.          DO i=its,itf
    513. 

  513.              muv(i,j) =      mu(i,j-1)        +mub(i,j-1)
    514. 

  514.          ENDDO
    515. 

  515.       END IF
    516. 

  516. 

  517. END SUBROUTINE calc_mu_staggered
    518. 

  518. 

  519. !-------------------------------------------------------------------------------
    520. 

  520. 

  521. SUBROUTINE couple ( mu, mub, rfield, field, name, &
    522. 

  522.                     msf,                          &
    523. 

  523.                     ids, ide, jds, jde, kds, kde, &
    524. 

  524.                     ims, ime, jms, jme, kms, kme, &
    525. 

  525.                     its, ite, jts, jte, kts, kte )
    526. 

  526. 

  527.    IMPLICIT NONE
    528. 

  528. 

  529.    ! Input data
    530. 

  530. 

  531.    INTEGER ,             INTENT(IN   ) :: ids, ide, jds, jde, kds, kde, &
    532. 

  532.                                           ims, ime, jms, jme, kms, kme, &
    533. 

  533.                                           its, ite, jts, jte, kts, kte
    534. 

  534. 

  535.    CHARACTER(LEN=1) ,     INTENT(IN   ) :: name
    536. 

  536. 

  537.    REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(  OUT) :: rfield
    538. 

  538. 

  539.    REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN   ) :: mu, mub, msf
    540. 

  540.    
    541. 

  541.    REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN   ) :: field
    542. 

  542.    
    543. 

  543.    ! Local data
    544. 

  544.    
    545. 

  545.    INTEGER :: i, j, k, itf, jtf, ktf
    546. 

  546.    REAL , DIMENSION(ims:ime,jms:jme) :: muu , muv
    547. 

  547. 

  548. !<DESCRIPTION>
    549. 

  549. !
    550. 

  550. ! subroutine couple couples the input variable with the dry-air 
    551. 

  551. ! column mass (mu).  
    552. 

  552. !
    553. 

  553. !</DESCRIPTION>
    554. 

  554. 

  555.    
    556. 

  556.    ktf=MIN(kte,kde-1)
    557. 

  557.    
    558. 

  558.    IF (name .EQ. 'u')THEN
    559. 

  559. 

  560.       CALL calc_mu_staggered ( mu, mub, muu, muv,            &
    561. 

  561.                                   ids, ide, jds, jde, kds, kde, &
    562. 

  562.                                   ims, ime, jms, jme, kms, kme, &
    563. 

  563.                                   its, ite, jts, jte, kts, kte )
    564. 

  564. 

  565.       itf=ite
    566. 

  566.       jtf=MIN(jte,jde-1)
    567. 

  567. 

  568.       DO j=jts,jtf
    569. 

  569.       DO k=kts,ktf
    570. 

  570.       DO i=its,itf
    571. 

  571.          rfield(i,k,j)=field(i,k,j)*muu(i,j)/msf(i,j)
    572. 

  572.       ENDDO
    573. 

  573.       ENDDO
    574. 

  574.       ENDDO
    575. 

  575. 

  576.    ELSE IF (name .EQ. 'v')THEN
    577. 

  577. 

  578.       CALL calc_mu_staggered ( mu, mub, muu, muv,            &
    579. 

  579.                                ids, ide, jds, jde, kds, kde, &
    580. 

  580.                                ims, ime, jms, jme, kms, kme, &
    581. 

  581.                                its, ite, jts, jte, kts, kte )
    582. 

  582. 

  583.       itf=ite
    584. 

  584.       itf=MIN(ite,ide-1)
    585. 

  585.       jtf=jte
    586. 

  586. 

  587.       DO j=jts,jtf
    588. 

  588.       DO k=kts,ktf
    589. 

  589.       DO i=its,itf
    590. 

  590.            rfield(i,k,j)=field(i,k,j)*muv(i,j)/msf(i,j)
    591. 

  591.       ENDDO
    592. 

  592.       ENDDO
    593. 

  593.       ENDDO
    594. 

  594. 

  595.    ELSE IF (name .EQ. 'w')THEN
    596. 

  596.       itf=MIN(ite,ide-1)
    597. 

  597.       jtf=MIN(jte,jde-1)
    598. 

  598.       DO j=jts,jtf
    599. 

  599.       DO k=kts,kte
    600. 

  600.       DO i=its,itf
    601. 

  601.          rfield(i,k,j)=field(i,k,j)*(mu(i,j)+mub(i,j))/msf(i,j)
    602. 

  602.       ENDDO
    603. 

  603.       ENDDO
    604. 

  604.       ENDDO
    605. 

  605. 

  606.    ELSE IF (name .EQ. 'h')THEN
    607. 

  607.       itf=MIN(ite,ide-1)
    608. 

  608.       jtf=MIN(jte,jde-1)
    609. 

  609.       DO j=jts,jtf
    610. 

  610.       DO k=kts,kte
    611. 

  611.       DO i=its,itf
    612. 

  612.          rfield(i,k,j)=field(i,k,j)*(mu(i,j)+mub(i,j))
    613. 

  613.       ENDDO
    614. 

  614.       ENDDO
    615. 

  615.       ENDDO
    616. 

  616. 

  617.    ELSE 
    618. 

  618.       itf=MIN(ite,ide-1)
    619. 

  619.       jtf=MIN(jte,jde-1)
    620. 

  620.       DO j=jts,jtf
    621. 

  621.       DO k=kts,ktf
    622. 

  622.       DO i=its,itf
    623. 

  623.          rfield(i,k,j)=field(i,k,j)*(mu(i,j)+mub(i,j))
    624. 

  624.       ENDDO
    625. 

  625.       ENDDO
    626. 

  626.       ENDDO
    627. 

  627.    
    628. 

  628.    ENDIF
    629. 

  629. 

  630. END SUBROUTINE couple
    631. 

  631. 

  632. 

  633. !-------------------------------------------------------------------------------
    634. 

  634. 

  635. SUBROUTINE calc_ww_cp ( u, v, mup, mub, ww,              &
    636. 

  636.                         rdx, rdy, msftx, msfty,          &
    637. 

  637.                         msfux, msfuy, msfvx, msfvx_inv,  &
    638. 

  638.                         msfvy, dnw,                      &
    639. 

  639.                         ids, ide, jds, jde, kds, kde,    &
    640. 

  640.                         ims, ime, jms, jme, kms, kme,    &
    641. 

  641.                         its, ite, jts, jte, kts, kte    )
    642. 

  642. 

  643.    IMPLICIT NONE
    644. 

  644. 

  645.    ! Input data
    646. 

  646. 

  647. 

  648.    INTEGER ,    INTENT(IN   ) :: ids, ide, jds, jde, kds, kde, &
    649. 

  649.                                  ims, ime, jms, jme, kms, kme, &
    650. 

  650.                                  its, ite, jts, jte, kts, kte
    651. 

  651. 

  652.    REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN   ) :: u, v
    653. 

  653.    REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN   ) :: mup, mub, &
    654. 

  654.                                                             msftx, msfty, &
    655. 

  655.                                                             msfux, msfuy, &
    656. 

  656.                                                             msfvx, msfvy, &
    657. 

  657.                                                             msfvx_inv
    658. 

  658.    REAL , DIMENSION( kms:kme ) , INTENT(IN   ) :: dnw
    659. 

  659.    
    660. 

  660.    REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT  ) :: ww
    661. 

  661.    REAL , INTENT(IN   )  :: rdx, rdy
    662. 

  662.    
    663. 

  663.    ! Local data
    664. 

  664.    
    665. 

  665.    INTEGER :: i, j, k, itf, jtf, ktf
    666. 

  666.    REAL , DIMENSION( its:ite ) :: dmdt
    667. 

  667.    REAL , DIMENSION( its:ite, kts:kte ) :: divv
    668. 

  668.    REAL , DIMENSION( its:ite+1, jts:jte+1 ) :: muu, muv
    669. 

  669. 

  670. !<DESCRIPTION>
    671. 

  671. !
    672. 

  672. !  calc_ww calculates omega using the velocities (u,v) and the dry-air 
    673. 

  673. !  column mass (mup+mub).
    674. 

  674. !  The algorithm integrates the continuity equation through the column
    675. 

  675. !  followed by a diagnosis of omega.
    676. 

  676. !
    677. 

  677. !</DESCRIPTION>
    678. 

  678. 

  679. !<DESCRIPTION>
    680. 

  680. !
    681. 

  681. !  calc_ww_cp calculates omega using the velocities (u,v) and the 
    682. 

  682. !  column mass mu.
    683. 

  683. !
    684. 

  684. !</DESCRIPTION>
    685. 

  685. 

  686.     jtf=MIN(jte,jde-1)
    687. 

  687.     ktf=MIN(kte,kde-1)  
    688. 

  688.     itf=MIN(ite,ide-1)
    689. 

  689. 

  690. !  mu coupled with the appropriate map factor
    691. 

  691. 

  692.       DO j=jts,jtf
    693. 

  693.       DO i=its,min(ite+1,ide)
    694. 

  694.         ! u is always coupled with my
    695. 

  695.         muu(i,j) = 0.5*(mup(i,j)+mub(i,j)+mup(i-1,j)+mub(i-1,j))/msfuy(i,j)
    696. 

  696.       ENDDO
    697. 

  697.       ENDDO
    698. 

  698. 

  699.       DO j=jts,min(jte+1,jde)
    700. 

  700.       DO i=its,itf
    701. 

  701.        ! v is always coupled with mx
    702. 

  702. !        muv(i,j) = 0.5*(mup(i,j)+mub(i,j)+mup(i,j-1)+mub(i,j-1))/msfvx(i,j)
    703. 

  703.         muv(i,j) = 0.5*(mup(i,j)+mub(i,j)+mup(i,j-1)+mub(i,j-1))*msfvx_inv(i,j)
    704. 

  704.       ENDDO
    705. 

  705.       ENDDO
    706. 

  706. 

  707.       DO j=jts,jtf
    708. 

  708. 

  709.         DO i=its,ite
    710. 

  710.           dmdt(i) = 0.
    711. 

  711.           ww(i,1,j) = 0.
    712. 

  712.           ww(i,kte,j) = 0.
    713. 

  713.         ENDDO
    714. 

  714. 

  715. !       Comments on the modifications for map scale factors
    716. 

  716. !       ADT eqn 47 / my (putting rho -> 'mu') is:
    717. 

  717. !       (1/my) partial d mu/dt = -mx partial d/dx(mu u/my) 
    718. 

  718. !                                -mx partial d/dy(mu v/mx)
    719. 

  719. !                                -partial d/dz(mu w/my)
    720. 

  720. !
    721. 

  721. !       Using nu instead of z the last term becomes:
    722. 

  722. !                                -partial d/dnu(mu (dnu/dt)/my)
    723. 

  723. !
    724. 

  724. !       Integrating with respect to nu over ALL levels, with dnu/dt=0 at top
    725. 

  725. !       and bottom, the last term becomes = 0
    726. 

  726. !
    727. 

  727. !       Integral|bot->top[(1/my) partial d mu/dt]dnu = 
    728. 

  728. !       Integral|bot->top[-mx partial d/dx(mu u/my) 
    729. 

  729. !                         -mx partial d/dy(mu v/mx)]dnu
    730. 

  730. !
    731. 

  731. !       muu='mu'[on u]/my, muv='mu'[on v]/mx
    732. 

  732. !       (1/my) partial d mu/dt is independent of nu
    733. 

  733. !         => LHS = Integral|bot->top[con]dnu = conservation*(-1) = -dmdt
    734. 

  734. !
    735. 

  735. !         => dmdt = mx*Integral|bot->top[partial d/dx(mu u/my) +
    736. 

  736. !                                        partial d/dy(mu v/mx)]dnu
    737. 

  737. !         => dmdt = sum_bot->top[divv]
    738. 

  738. !       where
    739. 

  739. !         divv=mx*[partial d/dx(mu u/my) + partial d/dy(mu v/mx)]*delta nu
    740. 

  740. 

  741.         DO k=kts,ktf
    742. 

  742.         DO i=its,itf
    743. 

  743. 

  744.           divv(i,k) = msftx(i,j)*dnw(k)*( rdx*(muu(i+1,j)*u(i+1,k,j)-muu(i,j)*u(i,k,j))  &
    745. 

  745.                                         +rdy*(muv(i,j+1)*v(i,k,j+1)-muv(i,j)*v(i,k,j))   )
    746. 

  746. 

  747. !          dmdt(i) = dmdt(i) + dnw(k)* ( rdx*(ru(i+1,k,j)-ru(i,k,j))  &
    748. 

  748. !                                       +rdy*(rv(i,k,j+1)-rv(i,k,j))   )
    749. 

  749. 

  750.           dmdt(i) = dmdt(i) + divv(i,k)
    751. 

  751. 

  752. 

  753.         ENDDO
    754. 

  754.         ENDDO
    755. 

  755. 

  756. !       Further map scale factor notes:
    757. 

  757. !       Now integrate from bottom to top, level by level:
    758. 

  758. !       mu dnu/dt/my [k+1] = mu dnu/dt/my [k] + [-(1/my) partial d mu/dt 
    759. 

  759. !                           -mx partial d/dx(mu u/my)
    760. 

  760. !                           -mx partial d/dy(mu v/mx)]*dnu[k->k+1]
    761. 

  761. !       ww [k+1] = ww [k] -(1/my) partial d mu/dt * dnu[k->k+1] - divv[k]
    762. 

  762. !                = ww [k] -dmdt * dnw[k] - divv[k]
    763. 

  763. 

  764.         DO k=2,ktf
    765. 

  765.         DO i=its,itf
    766. 

  766. 

  767. !           ww(i,k,j)=ww(i,k-1,j)                                       &
    768. 

  768. !                        - dnw(k-1)* ( dmdt(i)                          &
    769. 

  769. !                                     +rdx*(ru(i+1,k-1,j)-ru(i,k-1,j))  &
    770. 

  770. !                                     +rdy*(rv(i,k-1,j+1)-rv(i,k-1,j)) )
    771. 

  771. 

  772.            ww(i,k,j)=ww(i,k-1,j) - dnw(k-1)*dmdt(i) - divv(i,k-1)
    773. 

  773. 

  774.         ENDDO
    775. 

  775.         ENDDO
    776. 

  776.      ENDDO
    777. 

  777. 

  778. 

  779. END SUBROUTINE calc_ww_cp
    780. 

  780. 

  781. 

  782. !-------------------------------------------------------------------------------
    783. 

  783.  
    784. 

  784. SUBROUTINE calc_cq ( moist, cqu, cqv, cqw, n_moist, &
    785. 

  785.                      ids, ide, jds, jde, kds, kde,  &
    786. 

  786.                      ims, ime, jms, jme, kms, kme,  &
    787. 

  787.                      its, ite, jts, jte, kts, kte  )
    788. 

  788. 

  789.    IMPLICIT NONE
    790. 

  790.    
    791. 

  791.    ! Input data
    792. 

  792. 

  793.    INTEGER ,          INTENT(IN   ) :: ids, ide, jds, jde, kds, kde, &
    794. 

  794.                                        ims, ime, jms, jme, kms, kme, &
    795. 

  795.                                        its, ite, jts, jte, kts, kte 
    796. 

  796. 

  797.    INTEGER ,          INTENT(IN   ) :: n_moist
    798. 

  798.    
    799. 

  799. 

  800.    REAL, DIMENSION( ims:ime, kms:kme , jms:jme , n_moist ), INTENT(IN   ) :: moist
    801. 

  801.                                               
    802. 

  802.    REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(  OUT) :: cqu, cqv, cqw
    803. 

  803. 

  804.    ! Local stuff
    805. 

  805. 

  806.    REAL :: qtot
    807. 

  807.    
    808. 

  808.    INTEGER :: i, j, k, itf, jtf, ktf, ispe
    809. 

  809. 

  810. !<DESCRIPTION>
    811. 

  811. !
    812. 

  812. !  calc_cq calculates moist coefficients for the momentum equations.
    813. 

  813. !
    814. 

  814. !</DESCRIPTION>
    815. 

  815. 

  816.       itf=ite
    817. 

  817.       jtf=MIN(jte,jde-1)
    818. 

  818.       ktf=MIN(kte,kde-1)
    819. 

  819. 

  820.       IF(  n_moist >= PARAM_FIRST_SCALAR ) THEN
    821. 

  821. 

  822.         DO j=jts,jtf
    823. 

  823.         DO k=kts,ktf
    824. 

  824.         DO i=its,itf
    825. 

  825.           qtot = 0.
    826. 

  826. !DEC$ loop count(3)
    827. 

  827.           DO ispe=PARAM_FIRST_SCALAR,n_moist
    828. 

  828.             qtot = qtot + moist(i,k,j,ispe) + moist(i-1,k,j,ispe)
    829. 

  829.           ENDDO
    830. 

  830. !           qtot = 0.5*( moist(i  ,k,j,1)+moist(i  ,k,j,2)+moist(i  ,k,j,3)+  &
    831. 

  831. !     &                  moist(i-1,k,j,1)+moist(i-1,k,j,2)+moist(i-1,k,j,3) )
    832. 

  832. !           cqu(i,k,j) = 1./(1.+qtot)
    833. 

  833.            cqu(i,k,j) = 1./(1.+0.5*qtot)
    834. 

  834.         ENDDO
    835. 

  835.         ENDDO
    836. 

  836.         ENDDO
    837. 

  837. 

  838.         itf=MIN(ite,ide-1)
    839. 

  839.         jtf=jte
    840. 

  840. 

  841.         DO j=jts,jtf
    842. 

  842.         DO k=kts,ktf
    843. 

  843.         DO i=its,itf
    844. 

  844.           qtot = 0.
    845. 

  845. !DEC$ loop count(3)
    846. 

  846.           DO ispe=PARAM_FIRST_SCALAR,n_moist
    847. 

  847.             qtot = qtot + moist(i,k,j,ispe) + moist(i,k,j-1,ispe)
    848. 

  848.           ENDDO
    849. 

  849. !           qtot = 0.5*( moist(i,k,j  ,1)+moist(i,k,j  ,2)+moist(i,k,j  ,3)+  &
    850. 

  850. !     &                  moist(i,k,j-1,1)+moist(i,k,j-1,2)+moist(i,k,j-1,3) )
    851. 

  851. !           cqv(i,k,j) = 1./(1.+qtot)
    852. 

  852.            cqv(i,k,j) = 1./(1.+0.5*qtot)
    853. 

  853.         ENDDO
    854. 

  854.         ENDDO
    855. 

  855.         ENDDO
    856. 

  856. 

  857.         itf=MIN(ite,ide-1)
    858. 

  858.         jtf=MIN(jte,jde-1)
    859. 

  859.         DO j=jts,jtf
    860. 

  860.         DO k=kts+1,ktf
    861. 

  861.         DO i=its,itf
    862. 

  862.           qtot = 0.
    863. 

  863. !DEC$ loop count(3)
    864. 

  864.           DO ispe=PARAM_FIRST_SCALAR,n_moist
    865. 

  865.             qtot = qtot + moist(i,k,j,ispe) + moist(i,k-1,j,ispe)
    866. 

  866.           ENDDO
    867. 

  867. !           qtot = 0.5*( moist(i,k  ,j,1)+moist(i,k  ,j,2)+moist(i,k-1,j,3)+  &
    868. 

  868. !     &                  moist(i,k-1,j,1)+moist(i,k-1,j,2)+moist(i,k  ,j,3) )
    869. 

  869. !           cqw(i,k,j) = qtot
    870. 

  870.            cqw(i,k,j) = 0.5*qtot
    871. 

  871.         ENDDO
    872. 

  872.         ENDDO
    873. 

  873.         ENDDO
    874. 

  874. 

  875.       ELSE
    876. 

  876. 

  877.         DO j=jts,jtf
    878. 

  878.         DO k=kts,ktf
    879. 

  879.         DO i=its,itf
    880. 

  880.            cqu(i,k,j) = 1.
    881. 

  881.         ENDDO
    882. 

  882.         ENDDO
    883. 

  883.         ENDDO
    884. 

  884. 

  885.         itf=MIN(ite,ide-1)
    886. 

  886.         jtf=jte
    887. 

  887. 

  888.         DO j=jts,jtf
    889. 

  889.         DO k=kts,ktf
    890. 

  890.         DO i=its,itf
    891. 

  891.            cqv(i,k,j) = 1.
    892. 

  892.         ENDDO
    893. 

  893.         ENDDO
    894. 

  894.         ENDDO
    895. 

  895. 

  896.         itf=MIN(ite,ide-1)
    897. 

  897.         jtf=MIN(jte,jde-1)
    898. 

  898.         DO j=jts,jtf
    899. 

  899.         DO k=kts+1,ktf
    900. 

  900.         DO i=its,itf
    901. 

  901.            cqw(i,k,j) = 0.
    902. 

  902.         ENDDO
    903. 

  903.         ENDDO
    904. 

  904.         ENDDO
    905. 

  905. 

  906.       END IF
    907. 

  907. 

  908. END SUBROUTINE calc_cq
    909. 

  909. 

  910. !----------------------------------------------------------------------
    911. 

  911. 

  912. SUBROUTINE calc_alt ( alt, al, alb,                  &
    913. 

  913.                       ids, ide, jds, jde, kds, kde,  &
    914. 

  914.                       ims, ime, jms, jme, kms, kme,  &
    915. 

  915.                       its, ite, jts, jte, kts, kte  )
    916. 

  916. 

  917.    IMPLICIT NONE
    918. 

  918.    
    919. 

  919.    ! Input data
    920. 

  920. 

  921.    INTEGER ,          INTENT(IN   ) :: ids, ide, jds, jde, kds, kde, &
    922. 

  922.                                        ims, ime, jms, jme, kms, kme, &
    923. 

  923.                                        its, ite, jts, jte, kts, kte 
    924. 

  924. 

  925.    REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(IN   ) :: alb, al
    926. 

  926.    REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(  OUT) :: alt
    927. 

  927. 

  928.    ! Local stuff
    929. 

  929. 

  930.    INTEGER :: i, j, k, itf, jtf, ktf
    931. 

  931. 

  932. !<DESCRIPTION>
    933. 

  933. !
    934. 

  934. ! calc_alt computes the full inverse density
    935. 

  935. !
    936. 

  936. !</DESCRIPTION>
    937. 

  937. 

  938.       itf=MIN(ite,ide-1)
    939. 

  939.       jtf=MIN(jte,jde-1)
    940. 

  940.       ktf=MIN(kte,kde-1)
    941. 

  941. 

  942.       DO j=jts,jtf
    943. 

  943.       DO k=kts,ktf
    944. 

  944.       DO i=its,itf
    945. 

  945.         alt(i,k,j) = al(i,k,j)+alb(i,k,j)
    946. 

  946.       ENDDO
    947. 

  947.       ENDDO
    948. 

  948.       ENDDO
    949. 

  949. 

  950. 

  951. END SUBROUTINE calc_alt
    952. 

  952. 

  953. !----------------------------------------------------------------------
    954. 

  954. 

  955. SUBROUTINE calc_p_rho_phi ( moist, n_moist, hypsometric_opt,      &
    956. 

  956.                             al, alb, mu, muts, ph, phb, p, pb,     &
    957. 

  957.                             t, p0, t0, ptop, znu, znw, dnw, rdnw,  &
    958. 

  958.                             rdn, non_hydrostatic,          &
    959. 

  959.                             ids, ide, jds, jde, kds, kde,  &
    960. 

  960.                             ims, ime, jms, jme, kms, kme,  &
    961. 

  961.                             its, ite, jts, jte, kts, kte  )
    962. 

  962. 

  963.   IMPLICIT NONE
    964. 

  964.    
    965. 

  965.    ! Input data
    966. 

  966. 

  967.   LOGICAL ,          INTENT(IN   ) :: non_hydrostatic
    968. 

  968. 

  969.   INTEGER ,          INTENT(IN   ) :: ids, ide, jds, jde, kds, kde, &
    970. 

  970.                                       ims, ime, jms, jme, kms, kme, &
    971. 

  971.                                       its, ite, jts, jte, kts, kte 
    972. 

  972. 

  973.   INTEGER ,          INTENT(IN   ) :: n_moist
    974. 

  974.   INTEGER ,          INTENT(IN   ) :: hypsometric_opt
    975. 

  975. 

  976.   REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(IN   ) :: alb,  &
    977. 

  977.                                                                    pb,   &
    978. 

  978.                                                                    t
    979. 

  979. 

  980.   REAL, DIMENSION( ims:ime , kms:kme , jms:jme, n_moist ), INTENT(IN   ) :: moist
    981. 

  981. 

  982.   REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(  OUT) :: al, p
    983. 

  983. 

  984.   REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(INOUT) :: ph, phb
    985. 

  985. 

  986.   REAL, DIMENSION( ims:ime , jms:jme ), INTENT(IN   ) :: mu, muts
    987. 

  987. 

  988.   REAL, DIMENSION( kms:kme ), INTENT(IN   ) :: znu, znw, dnw, rdnw, rdn
    989. 

  989. 

  990.   REAL,   INTENT(IN   ) :: t0, p0, ptop
    991. 

  991. 

  992.   ! Local stuff
    993. 

  993. 

  994.   INTEGER :: i, j, k, itf, jtf, ktf, ispe
    995. 

  995.   REAL    :: qvf, qtot, qf1, qf2
    996. 

  996.   REAL, DIMENSION( its:ite) :: temp,cpovcv_v
    997. 

  997.   REAL    :: pfu, phm, pfd
    998. 

  998. 

  999. !<DESCRIPTION>
    1000. 

  1000. !
    1001. 

  1001. ! For the nonhydrostatic option, calc_p_rho_phi calculates the
    1002. 

  1002. ! diagnostic quantities pressure and (inverse) density from the
    1003. 

  1003. ! prognostic variables using the equation of state.
    1004. 

  1004. !
    1005. 

  1005. ! For the hydrostatic option, calc_p_rho_phi calculates the
    1006. 

  1006. ! diagnostic quantities (inverse) density and geopotential from the
    1007. 

  1007. ! prognostic variables using the equation of state and the hydrostatic 
    1008. 

  1008. ! equation.
    1009. 

  1009. !
    1010. 

  1010. !</DESCRIPTION>
    1011. 

  1011. 

  1012.   itf=MIN(ite,ide-1)
    1013. 

  1013.   jtf=MIN(jte,jde-1)
    1014. 

  1014.   ktf=MIN(kte,kde-1)
    1015. 

  1015. 

  1016. #ifndef INTELMKL
    1017. 

  1017.   cpovcv_v = cpovcv
    1018. 

  1018. #endif
    1019. 

  1019. 

  1020.   IF (non_hydrostatic) THEN
    1021. 

  1021. 

  1022.       IF (hypsometric_opt == 1) THEN
    1023. 

  1023.         DO j=jts,jtf
    1024. 

  1024.         DO k=kts,ktf
    1025. 

  1025.         DO i=its,itf
    1026. 

  1026.           al(i,k,j)=-1./muts(i,j)*(alb(i,k,j)*mu(i,j) + rdnw(k)*(ph(i,k+1,j)-ph(i,k,j)))
    1027. 

  1027.         END DO
    1028. 

  1028.         END DO
    1029. 

  1029.         END DO
    1030. 

  1030.       ELSE IF (hypsometric_opt == 2) THEN
    1031. 

  1031. 

  1032.         ! The relation used to get specific volume, al, is: al = -dZ/dp,
    1033. 

  1033.         ! where dp = mut * d(eta). The pressure depth, dp, is replaced with
    1034. 

  1034.         ! p*(dp/p) ~ p*LOG((p+0.5dp)/(p-0.5dp)). Difference between dp and p*dLOG(p)
    1035. 

  1035.         ! is as follows: p*dLOG(p) - dp = 1/12*(dp/p)**3 + 1/90*(dp/p)**5 + ...
    1036. 

  1036.         ! Therefore, p*dLOG(p) is always larger than dp and the difference is
    1037. 

  1037.         ! in proportion to dp/p. TKW, 02/16/2010
    1038. 

  1038. 

  1039.         DO j=jts,jtf
    1040. 

  1040.         DO k=kts,ktf
    1041. 

  1041.         DO i=its,itf
    1042. 

  1042.           pfu = muts(i,j)*znw(k+1)+ptop
    1043. 

  1043.           pfd = muts(i,j)*znw(k)  +ptop
    1044. 

  1044.           phm = muts(i,j)*znu(k)  +ptop
    1045. 

  1045.           al(i,k,j) = (ph(i,k+1,j)-ph(i,k,j)+phb(i,k+1,j)-phb(i,k,j))/phm/LOG(pfd/pfu)-alb(i,k,j)
    1046. 

  1046.         END DO
    1047. 

  1047.         END DO
    1048. 

  1048.         END DO
    1049. 

  1049.       ELSE
    1050. 

  1050.         CALL wrf_error_fatal ( 'calc_p_rho_phi: hypsometric_opt should be 1 or 2' )
    1051. 

  1051.       END IF
    1052. 

  1052. 

  1053.       IF (n_moist >= PARAM_FIRST_SCALAR ) THEN  
    1054. 

  1054. 

  1055.         DO j=jts,jtf
    1056. 

  1056.         DO k=kts,ktf
    1057. 

  1057.         DO i=its,itf
    1058. 

  1058.           qvf = 1.+rvovrd*moist(i,k,j,P_QV)
    1059. 

  1059.           temp(i)=(r_d*(t0+t(i,k,j))*qvf)/(p0*(al(i,k,j)+alb(i,k,j)))
    1060. 

  1060.         ENDDO
    1061. 

  1061. #ifdef INTELMKL
    1062. 

  1062.         CALL VPOWX ( itf-its+1, temp(its), cpovcv, p(its,k,j) )
    1063. 

  1063. #else
    1064. 

  1064. ! use vector version from libmassv or from compat lib in frame/libmassv.F
    1065. 

  1065.         CALL VPOW  ( p(its,k,j), temp(its), cpovcv_v(its), itf-its+1 )
    1066. 

  1066. #endif
    1067. 

  1067.         DO i=its,itf
    1068. 

  1068.            p(i,k,j)= p(i,k,j)*p0-pb(i,k,j)
    1069. 

  1069.         ENDDO
    1070. 

  1070.         ENDDO
    1071. 

  1071.         ENDDO
    1072. 

  1072. 

  1073.       ELSE
    1074. 

  1074. 

  1075.         DO j=jts,jtf
    1076. 

  1076.         DO k=kts,ktf
    1077. 

  1077.         DO i=its,itf
    1078. 

  1078.           p(i,k,j)=p0*( (r_d*(t0+t(i,k,j)))/                     &
    1079. 

  1079.                         (p0*(al(i,k,j)+alb(i,k,j))) )**cpovcv  &
    1080. 

  1080.                            -pb(i,k,j)
    1081. 

  1081.         ENDDO
    1082. 

  1082.         ENDDO
    1083. 

  1083.         ENDDO
    1084. 

  1084. 

  1085.       END IF
    1086. 

  1086. 

  1087.    ELSE
    1088. 

  1088. 

  1089. !  hydrostatic pressure, al, and ph1 calc; WCS, 5 sept 2001
    1090. 

  1090. 

  1091. 

  1092.       IF (n_moist >= PARAM_FIRST_SCALAR ) THEN  
    1093. 

  1093. 

  1094.         DO j=jts,jtf
    1095. 

  1095. 

  1096.           k=ktf          ! top layer
    1097. 

  1097.           DO i=its,itf
    1098. 

  1098. 

  1099.             qtot = 0.
    1100. 

  1100.             DO ispe=PARAM_FIRST_SCALAR,n_moist
    1101. 

  1101.               qtot = qtot + moist(i,k,j,ispe)
    1102. 

  1102.             ENDDO
    1103. 

  1103.             qf2 = 1.
    1104. 

  1104.             qf1 = qtot*qf2
    1105. 

  1105. 

  1106.             p(i,k,j) = - 0.5*(mu(i,j)+qf1*muts(i,j))/rdnw(k)/qf2
    1107. 

  1107.             qvf = 1.+rvovrd*moist(i,k,j,P_QV)
    1108. 

  1108.             al(i,k,j) = (r_d/p1000mb)*(t(i,k,j)+t0)*qvf* &
    1109. 

  1109.                 (((p(i,k,j)+pb(i,k,j))/p1000mb)**cvpm) - alb(i,k,j)
    1110. 

  1110. 

  1111.           ENDDO
    1112. 

  1112. 

  1113.           DO k=ktf-1,kts,-1  ! remaining layers, integrate down
    1114. 

  1114.             DO i=its,itf
    1115. 

  1115. 

  1116.             qtot = 0.
    1117. 

  1117.             DO ispe=PARAM_FIRST_SCALAR,n_moist
    1118. 

  1118.               qtot = qtot + 0.5*(  moist(i,k  ,j,ispe) + moist(i,k+1,j,ispe) )
    1119. 

  1119.             ENDDO
    1120. 

  1120.             qf2 = 1.
    1121. 

  1121.             qf1 = qtot*qf2
    1122. 

  1122. 

  1123.             p(i,k,j) = p(i,k+1,j) - (mu(i,j) + qf1*muts(i,j))/qf2/rdn(k+1)
    1124. 

  1124.             qvf = 1.+rvovrd*moist(i,k,j,P_QV)
    1125. 

  1125.             al(i,k,j) = (r_d/p1000mb)*(t(i,k,j)+t0)*qvf* &
    1126. 

  1126.                         (((p(i,k,j)+pb(i,k,j))/p1000mb)**cvpm) - alb(i,k,j)
    1127. 

  1127.             ENDDO
    1128. 

  1128.           ENDDO
    1129. 

  1129. 

  1130.         ENDDO
    1131. 

  1131. 

  1132.       ELSE
    1133. 

  1133. 

  1134.         DO j=jts,jtf
    1135. 

  1135. 

  1136.           k=ktf          ! top layer
    1137. 

  1137.           DO i=its,itf
    1138. 

  1138. 

  1139.             qtot = 0.
    1140. 

  1140.             qf2 = 1.
    1141. 

  1141.             qf1 = qtot*qf2
    1142. 

  1142. 

  1143.             p(i,k,j) = - 0.5*(mu(i,j)+qf1*muts(i,j))/rdnw(k)/qf2
    1144. 

  1144.             qvf = 1.
    1145. 

  1145.             al(i,k,j) = (r_d/p1000mb)*(t(i,k,j)+t0)*qvf* &
    1146. 

  1146.                 (((p(i,k,j)+pb(i,k,j))/p1000mb)**cvpm) - alb(i,k,j)
    1147. 

  1147. 

  1148.           ENDDO
    1149. 

  1149. 

  1150.           DO k=ktf-1,kts,-1  ! remaining layers, integrate down
    1151. 

  1151.             DO i=its,itf
    1152. 

  1152. 

  1153.             qtot = 0.
    1154. 

  1154.             qf2 = 1.
    1155. 

  1155.             qf1 = qtot*qf2
    1156. 

  1156. 

  1157.             p(i,k,j) = p(i,k+1,j) - (mu(i,j) + qf1*muts(i,j))/qf2/rdn(k+1)
    1158. 

  1158.             qvf = 1.
    1159. 

  1159.             al(i,k,j) = (r_d/p1000mb)*(t(i,k,j)+t0)*qvf* &
    1160. 

  1160.                         (((p(i,k,j)+pb(i,k,j))/p1000mb)**cvpm) - alb(i,k,j)
    1161. 

  1161.             ENDDO
    1162. 

  1162.           ENDDO
    1163. 

  1163. 

  1164.         ENDDO
    1165. 

  1165. 

  1166.      END IF
    1167. 

  1167. 

  1168.      IF (hypsometric_opt == 1) THEN
    1169. 

  1169.         DO j=jts,jtf
    1170. 

  1170.           DO k=2,ktf+1  ! integrate hydrostatic equation for geopotential
    1171. 

  1171.             DO i=its,itf
    1172. 

  1172.               ph(i,k,j) = ph(i,k-1,j) - (dnw(k-1))*(           &
    1173. 

  1173.                            (muts(i,j))*al(i,k-1,j)+            &
    1174. 

  1174.                             mu(i,j)*alb(i,k-1,j)  )
    1175. 

  1175.             ENDDO
    1176. 

  1176.           ENDDO
    1177. 

  1177.         ENDDO
    1178. 

  1178.      ELSE 
    1179. 

  1179. 

  1180.      ! Revised hypsometric eq.: dZ=-al*p*dLOG(p), where p is dry pressure
    1181. 

  1181. 

  1182.       DO j=jts,jtf
    1183. 

  1183.         DO i=its,itf
    1184. 

  1184.            ph(i,kts,j) = phb(i,kts,j)
    1185. 

  1185.         END DO
    1186. 

  1186. 

  1187.         DO k=kts+1,ktf+1
    1188. 

  1188.           DO i=its,itf
    1189. 

  1189.             pfu = muts(i,j)*znw(k)  +ptop
    1190. 

  1190.             pfd = muts(i,j)*znw(k-1)+ptop
    1191. 

  1191.             phm = muts(i,j)*znu(k-1)+ptop
    1192. 

  1192.             ph(i,k,j) = ph(i,k-1,j) + (al(i,k-1,j)+alb(i,k-1,j))*phm*LOG(pfd/pfu)
    1193. 

  1193.           ENDDO
    1194. 

  1194.         ENDDO
    1195. 

  1195. 

  1196.         DO k=kts,ktf+1
    1197. 

  1197.           DO i=its,itf
    1198. 

  1198.              ph(i,k,j) = ph(i,k,j) - phb(i,k,j)
    1199. 

  1199.           END DO
    1200. 

  1200.         END DO
    1201. 

  1201.       END DO
    1202. 

  1202. 

  1203.      END IF
    1204. 

  1204. 

  1205.    END IF
    1206. 

  1206. 

  1207. END SUBROUTINE calc_p_rho_phi
    1208. 

  1208. 

  1209. !----------------------------------------------------------------------
    1210. 

  1210. 

  1211. SUBROUTINE calc_php ( php, ph, phb,                  &
    1212. 

  1212.                       ids, ide, jds, jde, kds, kde,  &
    1213. 

  1213.                       ims, ime, jms, jme, kms, kme,  &
    1214. 

  1214.                       its, ite, jts, jte, kts, kte  )
    1215. 

  1215. 

  1216.    IMPLICIT NONE
    1217. 

  1217.    
    1218. 

  1218.    ! Input data
    1219. 

  1219. 

  1220.    INTEGER ,          INTENT(IN   ) :: ids, ide, jds, jde, kds, kde, &
    1221. 

  1221.                                        ims, ime, jms, jme, kms, kme, &
    1222. 

  1222.                                        its, ite, jts, jte, kts, kte 
    1223. 

  1223. 

  1224.    REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(IN   ) :: phb, ph
    1225. 

  1225.    REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(  OUT) :: php
    1226. 

  1226. 

  1227.    ! Local stuff
    1228. 

  1228. 

  1229.    INTEGER :: i, j, k, itf, jtf, ktf
    1230. 

  1230. 

  1231. !<DESCRIPTION>
    1232. 

  1232. !
    1233. 

  1233. !  calc_php calculates the full geopotential from the reference state
    1234. 

  1234. !  geopotential and the perturbation geopotential (phb_ph).
    1235. 

  1235. !
    1236. 

  1236. !</DESCRIPTION>
    1237. 

  1237. 

  1238.       itf=MIN(ite,ide-1)
    1239. 

  1239.       jtf=MIN(jte,jde-1)
    1240. 

  1240.       ktf=MIN(kte,kde-1)
    1241. 

  1241. 

  1242.       DO j=jts,jtf
    1243. 

  1243.       DO k=kts,ktf
    1244. 

  1244.       DO i=its,itf
    1245. 

  1245.         php(i,k,j) = 0.5*(phb(i,k,j)+phb(i,k+1,j)+ph(i,k,j)+ph(i,k+1,j))
    1246. 

  1246.       ENDDO
    1247. 

  1247.       ENDDO
    1248. 

  1248.       ENDDO
    1249. 

  1249. 

  1250. END SUBROUTINE calc_php
    1251. 

  1251. 

  1252. !-------------------------------------------------------------------------------
    1253. 

  1253. 

  1254. SUBROUTINE diagnose_w( ph_tend, ph_new, ph_old, w, mu, dt,  &
    1255. 

  1255.                        u, v, ht,                            &
    1256. 

  1256.                        cf1, cf2, cf3, rdx, rdy,             &
    1257. 

  1257.                        msftx, msfty,                        &
    1258. 

  1258.                        ids, ide, jds, jde, kds, kde,        &
    1259. 

  1259.                        ims, ime, jms, jme, kms, kme,        &
    1260. 

  1260.                        its, ite, jts, jte, kts, kte        )
    1261. 

  1261. 

  1262.    IMPLICIT NONE
    1263. 

  1263. 

  1264.    INTEGER ,          INTENT(IN   ) :: ids, ide, jds, jde, kds, kde, &
    1265. 

  1265.                                        ims, ime, jms, jme, kms, kme, &
    1266. 

  1266.                                        its, ite, jts, jte, kts, kte 
    1267. 

  1267. 

  1268.    REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(IN   ) ::   ph_tend, &
    1269. 

  1269.                                                                      ph_new,  &
    1270. 

  1270.                                                                      ph_old,  &
    1271. 

  1271.                                                                      u,       &
    1272. 

  1272.                                                                      v
    1273. 

  1273. 

  1274. 

  1275.    REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(  OUT) :: w
    1276. 

  1276. 

  1277.    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN   ) :: mu, ht, msftx, msfty
    1278. 

  1278. 

  1279.    REAL, INTENT(IN   ) :: dt, cf1, cf2, cf3, rdx, rdy
    1280. 

  1280. 

  1281.    INTEGER :: i, j, k, itf, jtf
    1282. 

  1282. 

  1283.    itf=MIN(ite,ide-1)
    1284. 

  1284.    jtf=MIN(jte,jde-1)
    1285. 

  1285. 

  1286. !<DESCRIPTION>
    1287. 

  1287. !
    1288. 

  1288. ! diagnose_w diagnoses the vertical velocity from the geopoential equation.
    1289. 

  1289. ! Used with the hydrostatic option.
    1290. 

  1290. !
    1291. 

  1291. !</DESCRIPTION>
    1292. 

  1292. 

  1293.    DO j = jts, jtf
    1294. 

  1294. 

  1295. !  lower b.c. on w
    1296. 

  1296. 

  1297. !  Notes on map scale factors:
    1298. 

  1298. !  Chain rule: if Z=Z(X,Y) [true at the surface] then
    1299. 

  1299. !  dZ/dt = dZ/dX * dX/dt + dZ/dY * dY/dt, U=dX/dt, V=dY/dt
    1300. 

  1300. !  Using capitals to denote actual values
    1301. 

  1301. !  In mapped values, u=U, v=V, z=Z, 1/dX=mx/dx, 1/dY=my/dy
    1302. 

  1302. !    => w = dz/dt = mx u dz/dx + my v dz/dy
    1303. 

  1303. !  [where dz/dx is just the surface height change between x
    1304. 

  1304. !   gridpoints, and dz/dy is the change between y gridpoints]
    1305. 

  1305. !  [NB: cf1, cf2 and cf3 do vertical weighting of u or v values
    1306. 

  1306. !   nearest the surface]
    1307. 

  1307. 

  1308. !  Previously msft multiplied by rdy and rdx terms.
    1309. 

  1309. !  Now msfty multiplies rdy term, and msftx multiplies msftx term
    1310. 

  1310.      DO i = its, itf
    1311. 

  1311.          w(i,1,j)=  msfty(i,j)*.5*rdy*(                      &
    1312. 

  1312.                            (ht(i,j+1)-ht(i,j  ))             &
    1313. 

  1313.           *(cf1*v(i,1,j+1)+cf2*v(i,2,j+1)+cf3*v(i,3,j+1))    &
    1314. 

  1314.                           +(ht(i,j  )-ht(i,j-1))             &
    1315. 

  1315.           *(cf1*v(i,1,j  )+cf2*v(i,2,j  )+cf3*v(i,3,j  ))  ) &
    1316. 

  1316.                  +msftx(i,j)*.5*rdx*(                        &
    1317. 

  1317.                            (ht(i+1,j)-ht(i,j  ))             &
    1318. 

  1318.           *(cf1*u(i+1,1,j)+cf2*u(i+1,2,j)+cf3*u(i+1,3,j))    &
    1319. 

  1319.                           +(ht(i,j  )-ht(i-1,j))             &
    1320. 

  1320.           *(cf1*u(i  ,1,j)+cf2*u(i  ,2,j)+cf3*u(i  ,3,j))  )
    1321. 

  1321.      ENDDO
    1322. 

  1322. 

  1323. !  use geopotential equation to diagnose w
    1324. 

  1324. 

  1325. !  Further notes on map scale factors
    1326. 

  1326. !  If ph_tend contains:  -mx partial d/dx(mu rho u/my) 
    1327. 

  1327. !                        -mx partial d/dy(phi mu v/mx)
    1328. 

  1328. !                        -partial d/dz(phi mu w/my)
    1329. 

  1329. !  then phi eqn is: partial d/dt(mu phi/my) = ph_tend + mu g w/my
    1330. 

  1330. !    => w = [my/(mu*g)]*[partial d/dt(mu phi/my) - ph_tend]
    1331. 

  1331. 

  1332.      DO k = 2, kte
    1333. 

  1333.      DO i = its, itf
    1334. 

  1334.        w(i,k,j) =  msfty(i,j)*(  (ph_new(i,k,j)-ph_old(i,k,j))/dt       &
    1335. 

  1335.                                - ph_tend(i,k,j)/mu(i,j)        )/g 
    1336. 

  1336. 

  1337.      ENDDO
    1338. 

  1338.      ENDDO
    1339. 

  1339. 

  1340.    ENDDO
    1341. 

  1341. 

  1342. END SUBROUTINE diagnose_w
    1343. 

  1343. 

  1344. !-------------------------------------------------------------------------------
    1345. 

  1345. 

  1346. SUBROUTINE rhs_ph( ph_tend, u, v, ww,               &
    1347. 

  1347.                    ph, ph_old, phb, w,              &
    1348. 

  1348.                    mut, muu, muv,                   &
    1349. 

  1349.                    fnm, fnp,                        &
    1350. 

  1350.                    rdnw, cfn, cfn1, rdx, rdy,       &
    1351. 

  1351.                    msfux, msfuy, msfvx,             &
    1352. 

  1352.                    msfvx_inv, msfvy,                &
    1353. 

  1353.                    msftx, msfty,                    &
    1354. 

  1354.                    non_hydrostatic,                 &
    1355. 

  1355.                    config_flags,                    &
    1356. 

  1356.                    ids, ide, jds, jde, kds, kde,    &
    1357. 

  1357.                    ims, ime, jms, jme, kms, kme,    &
    1358. 

  1358.                    its, ite, jts, jte, kts, kte    )
    1359. 

  1359.    IMPLICIT NONE
    1360. 

  1360. 

  1361.    TYPE(grid_config_rec_type), INTENT(IN   ) :: config_flags
    1362. 

  1362. 

  1363.    INTEGER ,          INTENT(IN   ) :: ids, ide, jds, jde, kds, kde, &
    1364. 

  1364.                                        ims, ime, jms, jme, kms, kme, &
    1365. 

  1365.                                        its, ite, jts, jte, kts, kte 
    1366. 

  1366. 

  1367.    REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(IN   ) ::        &
    1368. 

  1368.                                                                      u,   &
    1369. 

  1369.                                                                      v,   &
    1370. 

  1370.                                                                      ww,  &
    1371. 

  1371.                                                                      ph,  &
    1372. 

  1372.                                                                      ph_old, &
    1373. 

  1373.                                                                      phb, & 
    1374. 

  1374.                                                                     w
    1375. 

  1375. 

  1376.    REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(INOUT) :: ph_tend
    1377. 

  1377. 

  1378.    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN   ) :: muu, muv, mut,   &
    1379. 

  1379.                                                             msfux, msfuy, &
    1380. 

  1380.                                                             msfvx, msfvy, &
    1381. 

  1381.                                                             msftx, msfty, &
    1382. 

  1382.                                                             msfvx_inv
    1383. 

  1383. 

  1384.    REAL, DIMENSION( kms:kme ), INTENT(IN   ) :: rdnw, fnm, fnp
    1385. 

  1385. 

  1386.    REAL,  INTENT(IN   ) :: cfn, cfn1, rdx, rdy
    1387. 

  1387. 

  1388.    LOGICAL,  INTENT(IN   )  ::  non_hydrostatic
    1389. 

  1389. 

  1390.    ! Local stuff
    1391. 

  1391. 

  1392.    INTEGER :: i, j, k, itf, jtf, ktf, kz, i_start, j_start
    1393. 

  1393.    REAL    :: ur, ul, ub, vr, vl, vb
    1394. 

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

  1395. 

  1396.    INTEGER :: advective_order
    1397. 

  1397. 

  1398.    LOGICAL :: specified
    1399. 

  1399. 

  1400. !<DESCRIPTION>
    1401. 

  1401. !
    1402. 

  1402. ! rhs_ph calculates the large-timestep tendency terms for the geopotential
    1403. 

  1403. ! equation.  These terms include the advection and "gw".  The geopotential
    1404. 

  1404. ! equation is cast in advective form, so we don't use the flux form advection
    1405. 

  1405. ! algorithms here.
    1406. 

  1406. !
    1407. 

  1407. !</DESCRIPTION>
    1408. 

  1408. 

  1409.    specified = .false.
    1410. 

  1410.    if(config_flags%specified .or. config_flags%nested) specified = .true.
    1411. 

  1411. 

  1412.    advective_order = config_flags%h_sca_adv_order 
    1413. 

  1413. 

  1414.    itf=MIN(ite,ide-1)
    1415. 

  1415.    jtf=MIN(jte,jde-1)
    1416. 

  1416.    ktf=MIN(kte,kde-1)
    1417. 

  1417. 

  1418. !  Notes on map scale factors (WCS, 2 march 2008)
    1419. 

  1419. !  phi equation is:   mu/my d/dt(phi) = -(1/my) mx mu u  d/dx(phi)
    1420. 

  1420. !                                       -(1/my) my mu v d/dy(phi)
    1421. 

  1421. !                                       - omega d/d_eta(phi)
    1422. 

  1422. !                                               +mu g w/my
    1423. 

  1423. !
    1424. 

  1424. !  A little further explanation...
    1425. 

  1425. !  The tendency term we are computing here is for mu/my d/dt(phi).  It is advective form 
    1426. 

  1426. !  but it is multiplied be mu/my.  It will be decoupled from (mu/my) when the implicit w-phi
    1427. 

  1427. !  solution is computed in subourine advance_w.  The formulation dates from the early 
    1428. 

  1428. !  days of the mass coordinate model when we were testing both a flux and an advective formulation
    1429. 

  1429. !  for the geopotential equation and different forms of the vertical momentum equation and the 
    1430. 

  1430. !  vertically implicit solver.
    1431. 

  1431. 

  1432. ! advective form for the geopotential equation
    1433. 

  1433. 

  1434.    DO j = jts, jtf
    1435. 

  1435. 

  1436.      DO k = 2, kte
    1437. 

  1437.      DO i = its, itf
    1438. 

  1438.           wdwn(i,k) = .5*(ww(i,k,j)+ww(i,k-1,j))*rdnw(k-1)               &
    1439. 

  1439.                         *(ph(i,k,j)-ph(i,k-1,j)+phb(i,k,j)-phb(i,k-1,j))
    1440. 

  1440.      ENDDO
    1441. 

  1441.      ENDDO
    1442. 

  1442. 

  1443. !  RHS term 3 is: - omega partial d/dnu(phi)
    1444. 

  1444. 

  1445.      DO k = 2, kte-1
    1446. 

  1446.      DO i = its, itf
    1447. 

  1447.            ph_tend(i,k,j) = ph_tend(i,k,j)                           &
    1448. 

  1448.                              - (fnm(k)*wdwn(i,k+1)+fnp(k)*wdwn(i,k))
    1449. 

  1449.      ENDDO
    1450. 

  1450.      ENDDO
    1451. 

  1451. 

  1452.    ENDDO
    1453. 

  1453. 

  1454.    IF (non_hydrostatic) THEN  ! add in "gw" term.
    1455. 

  1455.    DO j = jts, jtf            ! in hydrostatic mode, "gw" will be diagnosed
    1456. 

  1456.                               ! after the timestep to give us "w"
    1457. 

  1457.      DO i = its, itf
    1458. 

  1458.         ph_tend(i,kde,j) = 0.
    1459. 

  1459.      ENDDO
    1460. 

  1460. 

  1461.      DO k = 2, kte
    1462. 

  1462.      DO i = its, itf
    1463. 

  1463.         ! phi equation RHS term 4
    1464. 

  1464.         ph_tend(i,k,j) = ph_tend(i,k,j) + mut(i,j)*g*w(i,k,j)/msfty(i,j)
    1465. 

  1465.      ENDDO
    1466. 

  1466.      ENDDO
    1467. 

  1467. 

  1468.    ENDDO
    1469. 

  1469. 

  1470.    END IF
    1471. 

  1471. 

  1472. !  Notes on map scale factors:
    1473. 

  1473. !  RHS terms 1 and 2 are: -(1/my) mx u mu partial d/dx(phi)
    1474. 

  1474. !                         -(1/my) my v mu partial d/dy(phi)
    1475. 

  1475. 

  1476.    IF (advective_order <= 2) THEN
    1477. 

  1477. 

  1478. !  y (v) advection
    1479. 

  1479. 

  1480.    i_start = its
    1481. 

  1481.    j_start = jts
    1482. 

  1482.    itf=MIN(ite,ide-1)
    1483. 

  1483.    jtf=MIN(jte,jde-1)
    1484. 

  1484. 

  1485.    IF ( (config_flags%open_ys .or. specified) .and. jts == jds ) j_start = jts+1
    1486. 

  1486.    IF ( (config_flags%open_ye .or. specified) .and. jte == jde ) jtf = jtf-2
    1487. 

  1487. 

  1488.    DO j = j_start, jtf
    1489. 

  1489. 

  1490.      DO k = 2, kte-1
    1491. 

  1491.      DO i = i_start, itf
    1492. 

  1492.         ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdy/msfty(i,j))*          &
    1493. 

  1493.                  ( muv(i,j+1)*(v(i,k,j+1)+v(i,k-1,j+1))*msfvy(i,j+1)*   &
    1494. 

  1494.                   (phb(i,k,j+1)-phb(i,k,j  )+ph(i,k,j+1)-ph(i,k,j  ))   &
    1495. 

  1495.                   +muv(i,j  )*(v(i,k,j  )+v(i,k-1,j  ))*msfvy(i,j  )*   &
    1496. 

  1496.                   (phb(i,k,j  )-phb(i,k,j-1)+ph(i,k,j  )-ph(i,k,j-1)) )
    1497. 

  1497.      ENDDO
    1498. 

  1498.      ENDDO
    1499. 

  1499. 

  1500.      k = kte
    1501. 

  1501.      DO i = i_start, itf
    1502. 

  1502.         ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdy/msfty(i,j))*                        &
    1503. 

  1503.                   ( muv(i,j+1)*(cfn*v(i,k-1,j+1)+cfn1*v(i,k-2,j+1))*msfvy(i,j+1)*    &
    1504. 

  1504.                    (phb(i,k,j+1)-phb(i,k,j  )+ph(i,k,j+1)-ph(i,k,j  ))               &
    1505. 

  1505.                    +muv(i,j  )*(cfn*v(i,k-1,j  )+cfn1*v(i,k-2,j  ))*msfvy(i,j  )*    &
    1506. 

  1506.                    (phb(i,k,j  )-phb(i,k,j-1)+ph(i,k,j  )-ph(i,k,j-1))              )
    1507. 

  1507.      ENDDO
    1508. 

  1508. 

  1509.    ENDDO
    1510. 

  1510. 

  1511. !  x (u) advection
    1512. 

  1512. 

  1513.    i_start = its
    1514. 

  1514.    j_start = jts
    1515. 

  1515.    itf=MIN(ite,ide-1)
    1516. 

  1516.    jtf=MIN(jte,jde-1)
    1517. 

  1517. 

  1518.    IF ( (config_flags%open_xs .or. specified) .and. its == ids ) i_start = its+1
    1519. 

  1519.    IF ( (config_flags%open_xe .or. specified) .and. ite == ide ) itf = itf-2
    1520. 

  1520. 

  1521.    DO j = j_start, jtf
    1522. 

  1522. 

  1523.      DO k = 2, kte-1
    1524. 

  1524.      DO i = i_start, itf
    1525. 

  1525.         ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdx/msfty(i,j))*         &
    1526. 

  1526.                  ( muu(i+1,j)*(u(i+1,k,j)+u(i+1,k-1,j))*msfux(i+1,j)*  &
    1527. 

  1527.                   (phb(i+1,k,j)-phb(i  ,k,j)+ph(i+1,k,j)-ph(i  ,k,j))  &
    1528. 

  1528.                   +muu(i  ,j)*(u(i  ,k,j)+u(i  ,k-1,j))*msfux(i  ,j)*  &
    1529. 

  1529.                   (phb(i  ,k,j)-phb(i-1,k,j)+ph(i  ,k,j)-ph(i-1,k,j))  )
    1530. 

  1530.      ENDDO
    1531. 

  1531.      ENDDO
    1532. 

  1532.  
    1533. 

  1533.      k = kte
    1534. 

  1534.      DO i = i_start, itf
    1535. 

  1535.         ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdx/msfty(i,j))*                        &
    1536. 

  1536.                   ( muu(i+1,j)*(cfn*u(i+1,k-1,j)+cfn1*u(i+1,k-2,j))*msfux(i+1,j)*    &
    1537. 

  1537.                    (phb(i+1,k,j)-phb(i  ,k,j)+ph(i+1,k,j)-ph(i  ,k,j))               &
    1538. 

  1538.                    +muu(i  ,j)*(cfn*u(i  ,k-1,j)+cfn1*u(i  ,k-2,j))*msfux(  i,j)*    &
    1539. 

  1539.                    (phb(i  ,k,j)-phb(i-1,k,j)+ph(i  ,k,j)-ph(i-1,k,j))             )
    1540. 

  1540.      ENDDO
    1541. 

  1541. 

  1542.    ENDDO
    1543. 

  1543. 

  1544.    ELSE IF (advective_order <= 4) THEN
    1545. 

  1545. 

  1546. !  y (v) advection
    1547. 

  1547. 

  1548.    i_start = its
    1549. 

  1549.    j_start = jts
    1550. 

  1550.    itf=MIN(ite,ide-1)
    1551. 

  1551.    jtf=MIN(jte,jde-1)
    1552. 

  1552. 

  1553.    IF ( (config_flags%open_ys .or. specified) .and. jts == jds ) j_start = jts+2
    1554. 

  1554.    IF ( (config_flags%open_ye .or. specified) .and. jte == jde ) jtf = jtf-3
    1555. 

  1555. 

  1556.    DO j = j_start, jtf
    1557. 

  1557. 

  1558.      DO k = 2, kte-1
    1559. 

  1559.      DO i = i_start, itf
    1560. 

  1560.         ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdy/msfty(i,j))*(                     &
    1561. 

  1561.                  ( muv(i,j+1)*(v(i,k,j+1)+v(i,k-1,j+1))*msfvy(i,j+1)                &
    1562. 

  1562.                   +muv(i,j  )*(v(i,k,j  )+v(i,k-1,j  ))*msfvy(i,j  ))* (1./12.)*(   &
    1563. 

  1563.                     8.*(ph(i,k,j+1)-ph(i,k,j-1))                                    &
    1564. 

  1564.                       -(ph(i,k,j+2)-ph(i,k,j-2))                                    &
    1565. 

  1565.                    +8.*(phb(i,k,j+1)-phb(i,k,j-1))                                  &
    1566. 

  1566.                       -(phb(i,k,j+2)-phb(i,k,j-2))  )   )                
    1567. 

  1567. 

  1568. 

  1569.      ENDDO
    1570. 

  1570.      ENDDO
    1571. 

  1571. 

  1572.      k = kte
    1573. 

  1573.      DO i = i_start, itf
    1574. 

  1574.         ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdy/msfty(i,j))*(                                 &
    1575. 

  1575.                  ( muv(i,j+1)*(cfn*v(i,k-1,j+1)+cfn1*v(i,k-2,j+1))*msfvy(i,j+1)                &
    1576. 

  1576.                   +muv(i,j  )*(cfn*v(i,k-1,j  )+cfn1*v(i,k-2,j  ))*msfvy(i,j  ))* (1./12.)*(   &
    1577. 

  1577.                     8.*(ph(i,k,j+1)-ph(i,k,j-1))                                               &
    1578. 

  1578.                       -(ph(i,k,j+2)-ph(i,k,j-2))                                               &
    1579. 

  1579.                    +8.*(phb(i,k,j+1)-phb(i,k,j-1))                                             &
    1580. 

  1580.                       -(phb(i,k,j+2)-phb(i,k,j-2))  )   )                
    1581. 

  1581. 

  1582.      ENDDO
    1583. 

  1583. 

  1584.    ENDDO
    1585. 

  1585. 

  1586. !  pick up near boundary rows using 2nd order stencil for open and specified conditions
    1587. 

  1587. 

  1588.    IF ( (config_flags%open_ys .or. specified) .and. jts <= jds+1 )  THEN
    1589. 

  1589. 

  1590.      j = jds+1
    1591. 

  1591.      DO k = 2, kte-1
    1592. 

  1592.      DO i = i_start, itf
    1593. 

  1593.         ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdy/msfty(i,j))*          &
    1594. 

  1594.                  ( muv(i,j+1)*(v(i,k,j+1)+v(i,k-1,j+1))*msfvy(i,j+1)*   &
    1595. 

  1595.                   (phb(i,k,j+1)-phb(i,k,j  )+ph(i,k,j+1)-ph(i,k,j  ))   &
    1596. 

  1596.                   +muv(i,j  )*(v(i,k,j  )+v(i,k-1,j  ))*msfvy(i,j  )*   &
    1597. 

  1597.                   (phb(i,k,j  )-phb(i,k,j-1)+ph(i,k,j  )-ph(i,k,j-1)) )
    1598. 

  1598.      ENDDO
    1599. 

  1599.      ENDDO
    1600. 

  1600. 

  1601.      k = kte
    1602. 

  1602.      DO i = i_start, itf
    1603. 

  1603.         ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdy/msfty(i,j))*                        &
    1604. 

  1604.                   ( muv(i,j+1)*(cfn*v(i,k-1,j+1)+cfn1*v(i,k-2,j+1))*msfvy(i,j+1)*    &
    1605. 

  1605.                    (phb(i,k,j+1)-phb(i,k,j  )+ph(i,k,j+1)-ph(i,k,j  ))               &
    1606. 

  1606.                    +muv(i,j  )*(cfn*v(i,k-1,j  )+cfn1*v(i,k-2,j  ))*msfvy(i,j  )*    &
    1607. 

  1607.                    (phb(i,k,j  )-phb(i,k,j-1)+ph(i,k,j  )-ph(i,k,j-1))              )
    1608. 

  1608.      ENDDO
    1609. 

  1609. 

  1610.    END IF
    1611. 

  1611. 

  1612.    IF ( (config_fla…
    1613.

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