PageRenderTime 108ms CodeModel.GetById 22ms 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
  1. !wrf:MODEL_LAYER:DYNAMICS
  2. !
  3. #if (RWORDSIZE == 4)
  4. # define VPOWX vspowx
  5. # define VPOW vspow
  6. #else
  7. # define VPOWX vpowx
  8. # define VPOW vpow
  9. #endif
  10. MODULE module_big_step_utilities_em
  11. USE module_model_constants
  12. USE module_state_description, only: p_qg, p_qs, p_qi, gdscheme, tiedtkescheme, kfetascheme, g3scheme, &
  13. p_qv, param_first_scalar, p_qr, p_qc
  14. USE module_configure, ONLY : grid_config_rec_type
  15. CONTAINS
  16. !-------------------------------------------------------------------------------
  17. SUBROUTINE calc_mu_uv ( config_flags, &
  18. mu, mub, muu, muv, &
  19. ids, ide, jds, jde, kds, kde, &
  20. ims, ime, jms, jme, kms, kme, &
  21. its, ite, jts, jte, kts, kte )
  22. IMPLICIT NONE
  23. ! Input data
  24. TYPE(grid_config_rec_type ) , INTENT(IN ) :: config_flags
  25. INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
  26. ims, ime, jms, jme, kms, kme, &
  27. its, ite, jts, jte, kts, kte
  28. REAL, DIMENSION( ims:ime , jms:jme ) , INTENT( OUT) :: muu, muv
  29. REAL, DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mu, mub
  30. ! local stuff
  31. INTEGER :: i, j, itf, jtf, im, jm
  32. !<DESCRIPTION>
  33. !
  34. ! calc_mu_uv calculates the full column dry-air mass at the staggered
  35. ! horizontal velocity points (u,v) and places the results in muu and muv.
  36. ! This routine uses the reference state (mub) and perturbation state (mu)
  37. !
  38. !</DESCRIPTION>
  39. itf=ite
  40. jtf=MIN(jte,jde-1)
  41. IF ( ( its .NE. ids ) .AND. ( ite .NE. ide ) ) THEN
  42. DO j=jts,jtf
  43. DO i=its,itf
  44. muu(i,j) = 0.5*(mu(i,j)+mu(i-1,j)+mub(i,j)+mub(i-1,j))
  45. ENDDO
  46. ENDDO
  47. ELSE IF ( ( its .EQ. ids ) .AND. ( ite .NE. ide ) ) THEN
  48. DO j=jts,jtf
  49. DO i=its+1,itf
  50. muu(i,j) = 0.5*(mu(i,j)+mu(i-1,j)+mub(i,j)+mub(i-1,j))
  51. ENDDO
  52. ENDDO
  53. i=its
  54. im = its
  55. if(config_flags%periodic_x) im = its-1
  56. DO j=jts,jtf
  57. ! muu(i,j) = mu(i,j) +mub(i,j)
  58. ! fix for periodic b.c., 13 march 2004, wcs
  59. muu(i,j) = 0.5*(mu(i,j)+mu(im,j)+mub(i,j)+mub(im,j))
  60. ENDDO
  61. ELSE IF ( ( its .NE. ids ) .AND. ( ite .EQ. ide ) ) THEN
  62. DO j=jts,jtf
  63. DO i=its,itf-1
  64. muu(i,j) = 0.5*(mu(i,j)+mu(i-1,j)+mub(i,j)+mub(i-1,j))
  65. ENDDO
  66. ENDDO
  67. i=ite
  68. im = ite-1
  69. if(config_flags%periodic_x) im = ite
  70. DO j=jts,jtf
  71. ! muu(i,j) = mu(i-1,j) +mub(i-1,j)
  72. ! fix for periodic b.c., 13 march 2004, wcs
  73. muu(i,j) = 0.5*(mu(i-1,j)+mu(im,j)+mub(i-1,j)+mub(im,j))
  74. ENDDO
  75. ELSE IF ( ( its .EQ. ids ) .AND. ( ite .EQ. ide ) ) THEN
  76. DO j=jts,jtf
  77. DO i=its+1,itf-1
  78. muu(i,j) = 0.5*(mu(i,j)+mu(i-1,j)+mub(i,j)+mub(i-1,j))
  79. ENDDO
  80. ENDDO
  81. i=its
  82. im = its
  83. if(config_flags%periodic_x) im = its-1
  84. DO j=jts,jtf
  85. ! muu(i,j) = mu(i,j) +mub(i,j)
  86. ! fix for periodic b.c., 13 march 2004, wcs
  87. muu(i,j) = 0.5*(mu(i,j)+mu(im,j)+mub(i,j)+mub(im,j))
  88. ENDDO
  89. i=ite
  90. im = ite-1
  91. if(config_flags%periodic_x) im = ite
  92. DO j=jts,jtf
  93. ! muu(i,j) = mu(i-1,j) +mub(i-1,j)
  94. ! fix for periodic b.c., 13 march 2004, wcs
  95. muu(i,j) = 0.5*(mu(i-1,j)+mu(im,j)+mub(i-1,j)+mub(im,j))
  96. ENDDO
  97. END IF
  98. itf=MIN(ite,ide-1)
  99. jtf=jte
  100. IF ( ( jts .NE. jds ) .AND. ( jte .NE. jde ) ) THEN
  101. DO j=jts,jtf
  102. DO i=its,itf
  103. muv(i,j) = 0.5*(mu(i,j)+mu(i,j-1)+mub(i,j)+mub(i,j-1))
  104. ENDDO
  105. ENDDO
  106. ELSE IF ( ( jts .EQ. jds ) .AND. ( jte .NE. jde ) ) THEN
  107. DO j=jts+1,jtf
  108. DO i=its,itf
  109. muv(i,j) = 0.5*(mu(i,j)+mu(i,j-1)+mub(i,j)+mub(i,j-1))
  110. ENDDO
  111. ENDDO
  112. j=jts
  113. jm = jts
  114. if(config_flags%periodic_y) jm = jts-1
  115. DO i=its,itf
  116. ! muv(i,j) = mu(i,j) +mub(i,j)
  117. ! fix for periodic b.c., 13 march 2004, wcs
  118. muv(i,j) = 0.5*(mu(i,j)+mu(i,jm)+mub(i,j)+mub(i,jm))
  119. ENDDO
  120. ELSE IF ( ( jts .NE. jds ) .AND. ( jte .EQ. jde ) ) THEN
  121. DO j=jts,jtf-1
  122. DO i=its,itf
  123. muv(i,j) = 0.5*(mu(i,j)+mu(i,j-1)+mub(i,j)+mub(i,j-1))
  124. ENDDO
  125. ENDDO
  126. j=jte
  127. jm = jte-1
  128. if(config_flags%periodic_y) jm = jte
  129. DO i=its,itf
  130. muv(i,j) = mu(i,j-1) +mub(i,j-1)
  131. ! fix for periodic b.c., 13 march 2004, wcs
  132. muv(i,j) = 0.5*(mu(i,j-1)+mu(i,jm)+mub(i,j-1)+mub(i,jm))
  133. ENDDO
  134. ELSE IF ( ( jts .EQ. jds ) .AND. ( jte .EQ. jde ) ) THEN
  135. DO j=jts+1,jtf-1
  136. DO i=its,itf
  137. muv(i,j) = 0.5*(mu(i,j)+mu(i,j-1)+mub(i,j)+mub(i,j-1))
  138. ENDDO
  139. ENDDO
  140. j=jts
  141. jm = jts
  142. if(config_flags%periodic_y) jm = jts-1
  143. DO i=its,itf
  144. ! muv(i,j) = mu(i,j) +mub(i,j)
  145. ! fix for periodic b.c., 13 march 2004, wcs
  146. muv(i,j) = 0.5*(mu(i,j)+mu(i,jm)+mub(i,j)+mub(i,jm))
  147. ENDDO
  148. j=jte
  149. jm = jte-1
  150. if(config_flags%periodic_y) jm = jte
  151. DO i=its,itf
  152. ! muv(i,j) = mu(i,j-1) +mub(i,j-1)
  153. ! fix for periodic b.c., 13 march 2004, wcs
  154. muv(i,j) = 0.5*(mu(i,j-1)+mu(i,jm)+mub(i,j-1)+mub(i,jm))
  155. ENDDO
  156. END IF
  157. END SUBROUTINE calc_mu_uv
  158. !-------------------------------------------------------------------------------
  159. SUBROUTINE calc_mu_uv_1 ( config_flags, &
  160. mu, muu, muv, &
  161. ids, ide, jds, jde, kds, kde, &
  162. ims, ime, jms, jme, kms, kme, &
  163. its, ite, jts, jte, kts, kte )
  164. IMPLICIT NONE
  165. ! Input data
  166. TYPE(grid_config_rec_type ) , INTENT(IN ) :: config_flags
  167. INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
  168. ims, ime, jms, jme, kms, kme, &
  169. its, ite, jts, jte, kts, kte
  170. REAL, DIMENSION( ims:ime , jms:jme ) , INTENT( OUT) :: muu, muv
  171. REAL, DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mu
  172. ! local stuff
  173. INTEGER :: i, j, itf, jtf, im, jm
  174. !<DESCRIPTION>
  175. !
  176. ! calc_mu_uv calculates the full column dry-air mass at the staggered
  177. ! horizontal velocity points (u,v) and places the results in muu and muv.
  178. ! This routine uses the full state (mu)
  179. !
  180. !</DESCRIPTION>
  181. itf=ite
  182. jtf=MIN(jte,jde-1)
  183. IF ( ( its .NE. ids ) .AND. ( ite .NE. ide ) ) THEN
  184. DO j=jts,jtf
  185. DO i=its,itf
  186. muu(i,j) = 0.5*(mu(i,j)+mu(i-1,j))
  187. ENDDO
  188. ENDDO
  189. ELSE IF ( ( its .EQ. ids ) .AND. ( ite .NE. ide ) ) THEN
  190. DO j=jts,jtf
  191. DO i=its+1,itf
  192. muu(i,j) = 0.5*(mu(i,j)+mu(i-1,j))
  193. ENDDO
  194. ENDDO
  195. i=its
  196. im = its
  197. if(config_flags%periodic_x) im = its-1
  198. DO j=jts,jtf
  199. muu(i,j) = 0.5*(mu(i,j)+mu(im,j))
  200. ENDDO
  201. ELSE IF ( ( its .NE. ids ) .AND. ( ite .EQ. ide ) ) THEN
  202. DO j=jts,jtf
  203. DO i=its,itf-1
  204. muu(i,j) = 0.5*(mu(i,j)+mu(i-1,j))
  205. ENDDO
  206. ENDDO
  207. i=ite
  208. im = ite-1
  209. if(config_flags%periodic_x) im = ite
  210. DO j=jts,jtf
  211. muu(i,j) = 0.5*(mu(i-1,j)+mu(im,j))
  212. ENDDO
  213. ELSE IF ( ( its .EQ. ids ) .AND. ( ite .EQ. ide ) ) THEN
  214. DO j=jts,jtf
  215. DO i=its+1,itf-1
  216. muu(i,j) = 0.5*(mu(i,j)+mu(i-1,j))
  217. ENDDO
  218. ENDDO
  219. i=its
  220. im = its
  221. if(config_flags%periodic_x) im = its-1
  222. DO j=jts,jtf
  223. muu(i,j) = 0.5*(mu(i,j)+mu(im,j))
  224. ENDDO
  225. i=ite
  226. im = ite-1
  227. if(config_flags%periodic_x) im = ite
  228. DO j=jts,jtf
  229. muu(i,j) = 0.5*(mu(i-1,j)+mu(im,j))
  230. ENDDO
  231. END IF
  232. itf=MIN(ite,ide-1)
  233. jtf=jte
  234. IF ( ( jts .NE. jds ) .AND. ( jte .NE. jde ) ) THEN
  235. DO j=jts,jtf
  236. DO i=its,itf
  237. muv(i,j) = 0.5*(mu(i,j)+mu(i,j-1))
  238. ENDDO
  239. ENDDO
  240. ELSE IF ( ( jts .EQ. jds ) .AND. ( jte .NE. jde ) ) THEN
  241. DO j=jts+1,jtf
  242. DO i=its,itf
  243. muv(i,j) = 0.5*(mu(i,j)+mu(i,j-1))
  244. ENDDO
  245. ENDDO
  246. j=jts
  247. jm = jts
  248. if(config_flags%periodic_y) jm = jts-1
  249. DO i=its,itf
  250. muv(i,j) = 0.5*(mu(i,j)+mu(i,jm))
  251. ENDDO
  252. ELSE IF ( ( jts .NE. jds ) .AND. ( jte .EQ. jde ) ) THEN
  253. DO j=jts,jtf-1
  254. DO i=its,itf
  255. muv(i,j) = 0.5*(mu(i,j)+mu(i,j-1))
  256. ENDDO
  257. ENDDO
  258. j=jte
  259. jm = jte-1
  260. if(config_flags%periodic_y) jm = jte
  261. DO i=its,itf
  262. muv(i,j) = 0.5*(mu(i,j-1)+mu(i,jm))
  263. ENDDO
  264. ELSE IF ( ( jts .EQ. jds ) .AND. ( jte .EQ. jde ) ) THEN
  265. DO j=jts+1,jtf-1
  266. DO i=its,itf
  267. muv(i,j) = 0.5*(mu(i,j)+mu(i,j-1))
  268. ENDDO
  269. ENDDO
  270. j=jts
  271. jm = jts
  272. if(config_flags%periodic_y) jm = jts-1
  273. DO i=its,itf
  274. muv(i,j) = 0.5*(mu(i,j)+mu(i,jm))
  275. ENDDO
  276. j=jte
  277. jm = jte-1
  278. if(config_flags%periodic_y) jm = jte
  279. DO i=its,itf
  280. muv(i,j) = 0.5*(mu(i,j-1)+mu(i,jm))
  281. ENDDO
  282. END IF
  283. END SUBROUTINE calc_mu_uv_1
  284. !-------------------------------------------------------------------------------
  285. ! Map scale factor comments for this routine:
  286. ! Locally not changed, but sent the correct map scale factors
  287. ! from module_em (msfuy, msfvx, msfty)
  288. SUBROUTINE couple_momentum ( muu, ru, u, msfu, &
  289. muv, rv, v, msfv, msfv_inv, &
  290. mut, rw, w, msft, &
  291. ids, ide, jds, jde, kds, kde, &
  292. ims, ime, jms, jme, kms, kme, &
  293. its, ite, jts, jte, kts, kte )
  294. IMPLICIT NONE
  295. ! Input data
  296. INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
  297. ims, ime, jms, jme, kms, kme, &
  298. its, ite, jts, jte, kts, kte
  299. REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT( OUT) :: ru, rv, rw
  300. REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: muu, muv, mut
  301. REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: msfu, msfv, msft, msfv_inv
  302. REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: u, v, w
  303. ! Local data
  304. INTEGER :: i, j, k, itf, jtf, ktf
  305. !<DESCRIPTION>
  306. !
  307. ! couple_momentum couples the velocities to the full column mass and
  308. ! the map factors.
  309. !
  310. !</DESCRIPTION>
  311. ktf=MIN(kte,kde-1)
  312. itf=ite
  313. jtf=MIN(jte,jde-1)
  314. DO j=jts,jtf
  315. DO k=kts,ktf
  316. DO i=its,itf
  317. ru(i,k,j)=u(i,k,j)*muu(i,j)/msfu(i,j)
  318. ENDDO
  319. ENDDO
  320. ENDDO
  321. itf=MIN(ite,ide-1)
  322. jtf=jte
  323. DO j=jts,jtf
  324. DO k=kts,ktf
  325. DO i=its,itf
  326. rv(i,k,j)=v(i,k,j)*muv(i,j)*msfv_inv(i,j)
  327. ! rv(i,k,j)=v(i,k,j)*muv(i,j)/msfv(i,j)
  328. ENDDO
  329. ENDDO
  330. ENDDO
  331. itf=MIN(ite,ide-1)
  332. jtf=MIN(jte,jde-1)
  333. DO j=jts,jtf
  334. DO k=kts,kte
  335. DO i=its,itf
  336. rw(i,k,j)=w(i,k,j)*mut(i,j)/msft(i,j)
  337. ENDDO
  338. ENDDO
  339. ENDDO
  340. END SUBROUTINE couple_momentum
  341. !-------------------------------------------------------------------
  342. SUBROUTINE calc_mu_staggered ( mu, mub, muu, muv, &
  343. ids, ide, jds, jde, kds, kde, &
  344. ims, ime, jms, jme, kms, kme, &
  345. its, ite, jts, jte, kts, kte )
  346. IMPLICIT NONE
  347. ! Input data
  348. INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
  349. ims, ime, jms, jme, kms, kme, &
  350. its, ite, jts, jte, kts, kte
  351. REAL, DIMENSION( ims:ime , jms:jme ) , INTENT( OUT) :: muu, muv
  352. REAL, DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mu, mub
  353. ! local stuff
  354. INTEGER :: i, j, itf, jtf
  355. !<DESCRIPTION>
  356. !
  357. ! calc_mu_staggered calculates the full dry air mass at the staggered
  358. ! velocity points (u,v).
  359. !
  360. !</DESCRIPTION>
  361. itf=ite
  362. jtf=MIN(jte,jde-1)
  363. IF ( ( its .NE. ids ) .AND. ( ite .NE. ide ) ) THEN
  364. DO j=jts,jtf
  365. DO i=its,itf
  366. muu(i,j) = 0.5*(mu(i,j)+mu(i-1,j)+mub(i,j)+mub(i-1,j))
  367. ENDDO
  368. ENDDO
  369. ELSE IF ( ( its .EQ. ids ) .AND. ( ite .NE. ide ) ) THEN
  370. DO j=jts,jtf
  371. DO i=its+1,itf
  372. muu(i,j) = 0.5*(mu(i,j)+mu(i-1,j)+mub(i,j)+mub(i-1,j))
  373. ENDDO
  374. ENDDO
  375. i=its
  376. DO j=jts,jtf
  377. muu(i,j) = mu(i,j) +mub(i,j)
  378. ENDDO
  379. ELSE IF ( ( its .NE. ids ) .AND. ( ite .EQ. ide ) ) THEN
  380. DO j=jts,jtf
  381. DO i=its,itf-1
  382. muu(i,j) = 0.5*(mu(i,j)+mu(i-1,j)+mub(i,j)+mub(i-1,j))
  383. ENDDO
  384. ENDDO
  385. i=ite
  386. DO j=jts,jtf
  387. muu(i,j) = mu(i-1,j) +mub(i-1,j)
  388. ENDDO
  389. ELSE IF ( ( its .EQ. ids ) .AND. ( ite .EQ. ide ) ) THEN
  390. DO j=jts,jtf
  391. DO i=its+1,itf-1
  392. muu(i,j) = 0.5*(mu(i,j)+mu(i-1,j)+mub(i,j)+mub(i-1,j))
  393. ENDDO
  394. ENDDO
  395. i=its
  396. DO j=jts,jtf
  397. muu(i,j) = mu(i,j) +mub(i,j)
  398. ENDDO
  399. i=ite
  400. DO j=jts,jtf
  401. muu(i,j) = mu(i-1,j) +mub(i-1,j)
  402. ENDDO
  403. END IF
  404. itf=MIN(ite,ide-1)
  405. jtf=jte
  406. IF ( ( jts .NE. jds ) .AND. ( jte .NE. jde ) ) THEN
  407. DO j=jts,jtf
  408. DO i=its,itf
  409. muv(i,j) = 0.5*(mu(i,j)+mu(i,j-1)+mub(i,j)+mub(i,j-1))
  410. ENDDO
  411. ENDDO
  412. ELSE IF ( ( jts .EQ. jds ) .AND. ( jte .NE. jde ) ) THEN
  413. DO j=jts+1,jtf
  414. DO i=its,itf
  415. muv(i,j) = 0.5*(mu(i,j)+mu(i,j-1)+mub(i,j)+mub(i,j-1))
  416. ENDDO
  417. ENDDO
  418. j=jts
  419. DO i=its,itf
  420. muv(i,j) = mu(i,j) +mub(i,j)
  421. ENDDO
  422. ELSE IF ( ( jts .NE. jds ) .AND. ( jte .EQ. jde ) ) THEN
  423. DO j=jts,jtf-1
  424. DO i=its,itf
  425. muv(i,j) = 0.5*(mu(i,j)+mu(i,j-1)+mub(i,j)+mub(i,j-1))
  426. ENDDO
  427. ENDDO
  428. j=jte
  429. DO i=its,itf
  430. muv(i,j) = mu(i,j-1) +mub(i,j-1)
  431. ENDDO
  432. ELSE IF ( ( jts .EQ. jds ) .AND. ( jte .EQ. jde ) ) THEN
  433. DO j=jts+1,jtf-1
  434. DO i=its,itf
  435. muv(i,j) = 0.5*(mu(i,j)+mu(i,j-1)+mub(i,j)+mub(i,j-1))
  436. ENDDO
  437. ENDDO
  438. j=jts
  439. DO i=its,itf
  440. muv(i,j) = mu(i,j) +mub(i,j)
  441. ENDDO
  442. j=jte
  443. DO i=its,itf
  444. muv(i,j) = mu(i,j-1) +mub(i,j-1)
  445. ENDDO
  446. END IF
  447. END SUBROUTINE calc_mu_staggered
  448. !-------------------------------------------------------------------------------
  449. SUBROUTINE couple ( mu, mub, rfield, field, name, &
  450. msf, &
  451. ids, ide, jds, jde, kds, kde, &
  452. ims, ime, jms, jme, kms, kme, &
  453. its, ite, jts, jte, kts, kte )
  454. IMPLICIT NONE
  455. ! Input data
  456. INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
  457. ims, ime, jms, jme, kms, kme, &
  458. its, ite, jts, jte, kts, kte
  459. CHARACTER(LEN=1) , INTENT(IN ) :: name
  460. REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT( OUT) :: rfield
  461. REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mu, mub, msf
  462. REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: field
  463. ! Local data
  464. INTEGER :: i, j, k, itf, jtf, ktf
  465. REAL , DIMENSION(ims:ime,jms:jme) :: muu , muv
  466. !<DESCRIPTION>
  467. !
  468. ! subroutine couple couples the input variable with the dry-air
  469. ! column mass (mu).
  470. !
  471. !</DESCRIPTION>
  472. ktf=MIN(kte,kde-1)
  473. IF (name .EQ. 'u')THEN
  474. CALL calc_mu_staggered ( mu, mub, muu, muv, &
  475. ids, ide, jds, jde, kds, kde, &
  476. ims, ime, jms, jme, kms, kme, &
  477. its, ite, jts, jte, kts, kte )
  478. itf=ite
  479. jtf=MIN(jte,jde-1)
  480. DO j=jts,jtf
  481. DO k=kts,ktf
  482. DO i=its,itf
  483. rfield(i,k,j)=field(i,k,j)*muu(i,j)/msf(i,j)
  484. ENDDO
  485. ENDDO
  486. ENDDO
  487. ELSE IF (name .EQ. 'v')THEN
  488. CALL calc_mu_staggered ( mu, mub, muu, muv, &
  489. ids, ide, jds, jde, kds, kde, &
  490. ims, ime, jms, jme, kms, kme, &
  491. its, ite, jts, jte, kts, kte )
  492. itf=ite
  493. itf=MIN(ite,ide-1)
  494. jtf=jte
  495. DO j=jts,jtf
  496. DO k=kts,ktf
  497. DO i=its,itf
  498. rfield(i,k,j)=field(i,k,j)*muv(i,j)/msf(i,j)
  499. ENDDO
  500. ENDDO
  501. ENDDO
  502. ELSE IF (name .EQ. 'w')THEN
  503. itf=MIN(ite,ide-1)
  504. jtf=MIN(jte,jde-1)
  505. DO j=jts,jtf
  506. DO k=kts,kte
  507. DO i=its,itf
  508. rfield(i,k,j)=field(i,k,j)*(mu(i,j)+mub(i,j))/msf(i,j)
  509. ENDDO
  510. ENDDO
  511. ENDDO
  512. ELSE IF (name .EQ. 'h')THEN
  513. itf=MIN(ite,ide-1)
  514. jtf=MIN(jte,jde-1)
  515. DO j=jts,jtf
  516. DO k=kts,kte
  517. DO i=its,itf
  518. rfield(i,k,j)=field(i,k,j)*(mu(i,j)+mub(i,j))
  519. ENDDO
  520. ENDDO
  521. ENDDO
  522. ELSE
  523. itf=MIN(ite,ide-1)
  524. jtf=MIN(jte,jde-1)
  525. DO j=jts,jtf
  526. DO k=kts,ktf
  527. DO i=its,itf
  528. rfield(i,k,j)=field(i,k,j)*(mu(i,j)+mub(i,j))
  529. ENDDO
  530. ENDDO
  531. ENDDO
  532. ENDIF
  533. END SUBROUTINE couple
  534. !-------------------------------------------------------------------------------
  535. SUBROUTINE calc_ww_cp ( u, v, mup, mub, ww, &
  536. rdx, rdy, msftx, msfty, &
  537. msfux, msfuy, msfvx, msfvx_inv, &
  538. msfvy, dnw, &
  539. ids, ide, jds, jde, kds, kde, &
  540. ims, ime, jms, jme, kms, kme, &
  541. its, ite, jts, jte, kts, kte )
  542. IMPLICIT NONE
  543. ! Input data
  544. INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
  545. ims, ime, jms, jme, kms, kme, &
  546. its, ite, jts, jte, kts, kte
  547. REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: u, v
  548. REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mup, mub, &
  549. msftx, msfty, &
  550. msfux, msfuy, &
  551. msfvx, msfvy, &
  552. msfvx_inv
  553. REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: dnw
  554. REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT ) :: ww
  555. REAL , INTENT(IN ) :: rdx, rdy
  556. ! Local data
  557. INTEGER :: i, j, k, itf, jtf, ktf
  558. REAL , DIMENSION( its:ite ) :: dmdt
  559. REAL , DIMENSION( its:ite, kts:kte ) :: divv
  560. REAL , DIMENSION( its:ite+1, jts:jte+1 ) :: muu, muv
  561. !<DESCRIPTION>
  562. !
  563. ! calc_ww calculates omega using the velocities (u,v) and the dry-air
  564. ! column mass (mup+mub).
  565. ! The algorithm integrates the continuity equation through the column
  566. ! followed by a diagnosis of omega.
  567. !
  568. !</DESCRIPTION>
  569. !<DESCRIPTION>
  570. !
  571. ! calc_ww_cp calculates omega using the velocities (u,v) and the
  572. ! column mass mu.
  573. !
  574. !</DESCRIPTION>
  575. jtf=MIN(jte,jde-1)
  576. ktf=MIN(kte,kde-1)
  577. itf=MIN(ite,ide-1)
  578. ! mu coupled with the appropriate map factor
  579. DO j=jts,jtf
  580. DO i=its,min(ite+1,ide)
  581. ! u is always coupled with my
  582. muu(i,j) = 0.5*(mup(i,j)+mub(i,j)+mup(i-1,j)+mub(i-1,j))/msfuy(i,j)
  583. ENDDO
  584. ENDDO
  585. DO j=jts,min(jte+1,jde)
  586. DO i=its,itf
  587. ! v is always coupled with mx
  588. ! muv(i,j) = 0.5*(mup(i,j)+mub(i,j)+mup(i,j-1)+mub(i,j-1))/msfvx(i,j)
  589. muv(i,j) = 0.5*(mup(i,j)+mub(i,j)+mup(i,j-1)+mub(i,j-1))*msfvx_inv(i,j)
  590. ENDDO
  591. ENDDO
  592. DO j=jts,jtf
  593. DO i=its,ite
  594. dmdt(i) = 0.
  595. ww(i,1,j) = 0.
  596. ww(i,kte,j) = 0.
  597. ENDDO
  598. ! Comments on the modifications for map scale factors
  599. ! ADT eqn 47 / my (putting rho -> 'mu') is:
  600. ! (1/my) partial d mu/dt = -mx partial d/dx(mu u/my)
  601. ! -mx partial d/dy(mu v/mx)
  602. ! -partial d/dz(mu w/my)
  603. !
  604. ! Using nu instead of z the last term becomes:
  605. ! -partial d/dnu(mu (dnu/dt)/my)
  606. !
  607. ! Integrating with respect to nu over ALL levels, with dnu/dt=0 at top
  608. ! and bottom, the last term becomes = 0
  609. !
  610. ! Integral|bot->top[(1/my) partial d mu/dt]dnu =
  611. ! Integral|bot->top[-mx partial d/dx(mu u/my)
  612. ! -mx partial d/dy(mu v/mx)]dnu
  613. !
  614. ! muu='mu'[on u]/my, muv='mu'[on v]/mx
  615. ! (1/my) partial d mu/dt is independent of nu
  616. ! => LHS = Integral|bot->top[con]dnu = conservation*(-1) = -dmdt
  617. !
  618. ! => dmdt = mx*Integral|bot->top[partial d/dx(mu u/my) +
  619. ! partial d/dy(mu v/mx)]dnu
  620. ! => dmdt = sum_bot->top[divv]
  621. ! where
  622. ! divv=mx*[partial d/dx(mu u/my) + partial d/dy(mu v/mx)]*delta nu
  623. DO k=kts,ktf
  624. DO i=its,itf
  625. 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)) &
  626. +rdy*(muv(i,j+1)*v(i,k,j+1)-muv(i,j)*v(i,k,j)) )
  627. ! dmdt(i) = dmdt(i) + dnw(k)* ( rdx*(ru(i+1,k,j)-ru(i,k,j)) &
  628. ! +rdy*(rv(i,k,j+1)-rv(i,k,j)) )
  629. dmdt(i) = dmdt(i) + divv(i,k)
  630. ENDDO
  631. ENDDO
  632. ! Further map scale factor notes:
  633. ! Now integrate from bottom to top, level by level:
  634. ! mu dnu/dt/my [k+1] = mu dnu/dt/my [k] + [-(1/my) partial d mu/dt
  635. ! -mx partial d/dx(mu u/my)
  636. ! -mx partial d/dy(mu v/mx)]*dnu[k->k+1]
  637. ! ww [k+1] = ww [k] -(1/my) partial d mu/dt * dnu[k->k+1] - divv[k]
  638. ! = ww [k] -dmdt * dnw[k] - divv[k]
  639. DO k=2,ktf
  640. DO i=its,itf
  641. ! ww(i,k,j)=ww(i,k-1,j) &
  642. ! - dnw(k-1)* ( dmdt(i) &
  643. ! +rdx*(ru(i+1,k-1,j)-ru(i,k-1,j)) &
  644. ! +rdy*(rv(i,k-1,j+1)-rv(i,k-1,j)) )
  645. ww(i,k,j)=ww(i,k-1,j) - dnw(k-1)*dmdt(i) - divv(i,k-1)
  646. ENDDO
  647. ENDDO
  648. ENDDO
  649. END SUBROUTINE calc_ww_cp
  650. !-------------------------------------------------------------------------------
  651. SUBROUTINE calc_cq ( moist, cqu, cqv, cqw, n_moist, &
  652. ids, ide, jds, jde, kds, kde, &
  653. ims, ime, jms, jme, kms, kme, &
  654. its, ite, jts, jte, kts, kte )
  655. IMPLICIT NONE
  656. ! Input data
  657. INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
  658. ims, ime, jms, jme, kms, kme, &
  659. its, ite, jts, jte, kts, kte
  660. INTEGER , INTENT(IN ) :: n_moist
  661. REAL, DIMENSION( ims:ime, kms:kme , jms:jme , n_moist ), INTENT(IN ) :: moist
  662. REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT( OUT) :: cqu, cqv, cqw
  663. ! Local stuff
  664. REAL :: qtot
  665. INTEGER :: i, j, k, itf, jtf, ktf, ispe
  666. !<DESCRIPTION>
  667. !
  668. ! calc_cq calculates moist coefficients for the momentum equations.
  669. !
  670. !</DESCRIPTION>
  671. itf=ite
  672. jtf=MIN(jte,jde-1)
  673. ktf=MIN(kte,kde-1)
  674. IF( n_moist >= PARAM_FIRST_SCALAR ) THEN
  675. DO j=jts,jtf
  676. DO k=kts,ktf
  677. DO i=its,itf
  678. qtot = 0.
  679. !DEC$ loop count(3)
  680. DO ispe=PARAM_FIRST_SCALAR,n_moist
  681. qtot = qtot + moist(i,k,j,ispe) + moist(i-1,k,j,ispe)
  682. ENDDO
  683. ! qtot = 0.5*( moist(i ,k,j,1)+moist(i ,k,j,2)+moist(i ,k,j,3)+ &
  684. ! & moist(i-1,k,j,1)+moist(i-1,k,j,2)+moist(i-1,k,j,3) )
  685. ! cqu(i,k,j) = 1./(1.+qtot)
  686. cqu(i,k,j) = 1./(1.+0.5*qtot)
  687. ENDDO
  688. ENDDO
  689. ENDDO
  690. itf=MIN(ite,ide-1)
  691. jtf=jte
  692. DO j=jts,jtf
  693. DO k=kts,ktf
  694. DO i=its,itf
  695. qtot = 0.
  696. !DEC$ loop count(3)
  697. DO ispe=PARAM_FIRST_SCALAR,n_moist
  698. qtot = qtot + moist(i,k,j,ispe) + moist(i,k,j-1,ispe)
  699. ENDDO
  700. ! qtot = 0.5*( moist(i,k,j ,1)+moist(i,k,j ,2)+moist(i,k,j ,3)+ &
  701. ! & moist(i,k,j-1,1)+moist(i,k,j-1,2)+moist(i,k,j-1,3) )
  702. ! cqv(i,k,j) = 1./(1.+qtot)
  703. cqv(i,k,j) = 1./(1.+0.5*qtot)
  704. ENDDO
  705. ENDDO
  706. ENDDO
  707. itf=MIN(ite,ide-1)
  708. jtf=MIN(jte,jde-1)
  709. DO j=jts,jtf
  710. DO k=kts+1,ktf
  711. DO i=its,itf
  712. qtot = 0.
  713. !DEC$ loop count(3)
  714. DO ispe=PARAM_FIRST_SCALAR,n_moist
  715. qtot = qtot + moist(i,k,j,ispe) + moist(i,k-1,j,ispe)
  716. ENDDO
  717. ! qtot = 0.5*( moist(i,k ,j,1)+moist(i,k ,j,2)+moist(i,k-1,j,3)+ &
  718. ! & moist(i,k-1,j,1)+moist(i,k-1,j,2)+moist(i,k ,j,3) )
  719. ! cqw(i,k,j) = qtot
  720. cqw(i,k,j) = 0.5*qtot
  721. ENDDO
  722. ENDDO
  723. ENDDO
  724. ELSE
  725. DO j=jts,jtf
  726. DO k=kts,ktf
  727. DO i=its,itf
  728. cqu(i,k,j) = 1.
  729. ENDDO
  730. ENDDO
  731. ENDDO
  732. itf=MIN(ite,ide-1)
  733. jtf=jte
  734. DO j=jts,jtf
  735. DO k=kts,ktf
  736. DO i=its,itf
  737. cqv(i,k,j) = 1.
  738. ENDDO
  739. ENDDO
  740. ENDDO
  741. itf=MIN(ite,ide-1)
  742. jtf=MIN(jte,jde-1)
  743. DO j=jts,jtf
  744. DO k=kts+1,ktf
  745. DO i=its,itf
  746. cqw(i,k,j) = 0.
  747. ENDDO
  748. ENDDO
  749. ENDDO
  750. END IF
  751. END SUBROUTINE calc_cq
  752. !----------------------------------------------------------------------
  753. SUBROUTINE calc_alt ( alt, al, alb, &
  754. ids, ide, jds, jde, kds, kde, &
  755. ims, ime, jms, jme, kms, kme, &
  756. its, ite, jts, jte, kts, kte )
  757. IMPLICIT NONE
  758. ! Input data
  759. INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
  760. ims, ime, jms, jme, kms, kme, &
  761. its, ite, jts, jte, kts, kte
  762. REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(IN ) :: alb, al
  763. REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT( OUT) :: alt
  764. ! Local stuff
  765. INTEGER :: i, j, k, itf, jtf, ktf
  766. !<DESCRIPTION>
  767. !
  768. ! calc_alt computes the full inverse density
  769. !
  770. !</DESCRIPTION>
  771. itf=MIN(ite,ide-1)
  772. jtf=MIN(jte,jde-1)
  773. ktf=MIN(kte,kde-1)
  774. DO j=jts,jtf
  775. DO k=kts,ktf
  776. DO i=its,itf
  777. alt(i,k,j) = al(i,k,j)+alb(i,k,j)
  778. ENDDO
  779. ENDDO
  780. ENDDO
  781. END SUBROUTINE calc_alt
  782. !----------------------------------------------------------------------
  783. SUBROUTINE calc_p_rho_phi ( moist, n_moist, hypsometric_opt, &
  784. al, alb, mu, muts, ph, phb, p, pb, &
  785. t, p0, t0, ptop, znu, znw, dnw, rdnw, &
  786. rdn, non_hydrostatic, &
  787. ids, ide, jds, jde, kds, kde, &
  788. ims, ime, jms, jme, kms, kme, &
  789. its, ite, jts, jte, kts, kte )
  790. IMPLICIT NONE
  791. ! Input data
  792. LOGICAL , INTENT(IN ) :: non_hydrostatic
  793. INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
  794. ims, ime, jms, jme, kms, kme, &
  795. its, ite, jts, jte, kts, kte
  796. INTEGER , INTENT(IN ) :: n_moist
  797. INTEGER , INTENT(IN ) :: hypsometric_opt
  798. REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(IN ) :: alb, &
  799. pb, &
  800. t
  801. REAL, DIMENSION( ims:ime , kms:kme , jms:jme, n_moist ), INTENT(IN ) :: moist
  802. REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT( OUT) :: al, p
  803. REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(INOUT) :: ph, phb
  804. REAL, DIMENSION( ims:ime , jms:jme ), INTENT(IN ) :: mu, muts
  805. REAL, DIMENSION( kms:kme ), INTENT(IN ) :: znu, znw, dnw, rdnw, rdn
  806. REAL, INTENT(IN ) :: t0, p0, ptop
  807. ! Local stuff
  808. INTEGER :: i, j, k, itf, jtf, ktf, ispe
  809. REAL :: qvf, qtot, qf1, qf2
  810. REAL, DIMENSION( its:ite) :: temp,cpovcv_v
  811. REAL :: pfu, phm, pfd
  812. !<DESCRIPTION>
  813. !
  814. ! For the nonhydrostatic option, calc_p_rho_phi calculates the
  815. ! diagnostic quantities pressure and (inverse) density from the
  816. ! prognostic variables using the equation of state.
  817. !
  818. ! For the hydrostatic option, calc_p_rho_phi calculates the
  819. ! diagnostic quantities (inverse) density and geopotential from the
  820. ! prognostic variables using the equation of state and the hydrostatic
  821. ! equation.
  822. !
  823. !</DESCRIPTION>
  824. itf=MIN(ite,ide-1)
  825. jtf=MIN(jte,jde-1)
  826. ktf=MIN(kte,kde-1)
  827. #ifndef INTELMKL
  828. cpovcv_v = cpovcv
  829. #endif
  830. IF (non_hydrostatic) THEN
  831. IF (hypsometric_opt == 1) THEN
  832. DO j=jts,jtf
  833. DO k=kts,ktf
  834. DO i=its,itf
  835. 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)))
  836. END DO
  837. END DO
  838. END DO
  839. ELSE IF (hypsometric_opt == 2) THEN
  840. ! The relation used to get specific volume, al, is: al = -dZ/dp,
  841. ! where dp = mut * d(eta). The pressure depth, dp, is replaced with
  842. ! p*(dp/p) ~ p*LOG((p+0.5dp)/(p-0.5dp)). Difference between dp and p*dLOG(p)
  843. ! is as follows: p*dLOG(p) - dp = 1/12*(dp/p)**3 + 1/90*(dp/p)**5 + ...
  844. ! Therefore, p*dLOG(p) is always larger than dp and the difference is
  845. ! in proportion to dp/p. TKW, 02/16/2010
  846. DO j=jts,jtf
  847. DO k=kts,ktf
  848. DO i=its,itf
  849. pfu = muts(i,j)*znw(k+1)+ptop
  850. pfd = muts(i,j)*znw(k) +ptop
  851. phm = muts(i,j)*znu(k) +ptop
  852. 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)
  853. END DO
  854. END DO
  855. END DO
  856. ELSE
  857. CALL wrf_error_fatal ( 'calc_p_rho_phi: hypsometric_opt should be 1 or 2' )
  858. END IF
  859. IF (n_moist >= PARAM_FIRST_SCALAR ) THEN
  860. DO j=jts,jtf
  861. DO k=kts,ktf
  862. DO i=its,itf
  863. qvf = 1.+rvovrd*moist(i,k,j,P_QV)
  864. temp(i)=(r_d*(t0+t(i,k,j))*qvf)/(p0*(al(i,k,j)+alb(i,k,j)))
  865. ENDDO
  866. #ifdef INTELMKL
  867. CALL VPOWX ( itf-its+1, temp(its), cpovcv, p(its,k,j) )
  868. #else
  869. ! use vector version from libmassv or from compat lib in frame/libmassv.F
  870. CALL VPOW ( p(its,k,j), temp(its), cpovcv_v(its), itf-its+1 )
  871. #endif
  872. DO i=its,itf
  873. p(i,k,j)= p(i,k,j)*p0-pb(i,k,j)
  874. ENDDO
  875. ENDDO
  876. ENDDO
  877. ELSE
  878. DO j=jts,jtf
  879. DO k=kts,ktf
  880. DO i=its,itf
  881. p(i,k,j)=p0*( (r_d*(t0+t(i,k,j)))/ &
  882. (p0*(al(i,k,j)+alb(i,k,j))) )**cpovcv &
  883. -pb(i,k,j)
  884. ENDDO
  885. ENDDO
  886. ENDDO
  887. END IF
  888. ELSE
  889. ! hydrostatic pressure, al, and ph1 calc; WCS, 5 sept 2001
  890. IF (n_moist >= PARAM_FIRST_SCALAR ) THEN
  891. DO j=jts,jtf
  892. k=ktf ! top layer
  893. DO i=its,itf
  894. qtot = 0.
  895. DO ispe=PARAM_FIRST_SCALAR,n_moist
  896. qtot = qtot + moist(i,k,j,ispe)
  897. ENDDO
  898. qf2 = 1.
  899. qf1 = qtot*qf2
  900. p(i,k,j) = - 0.5*(mu(i,j)+qf1*muts(i,j))/rdnw(k)/qf2
  901. qvf = 1.+rvovrd*moist(i,k,j,P_QV)
  902. al(i,k,j) = (r_d/p1000mb)*(t(i,k,j)+t0)*qvf* &
  903. (((p(i,k,j)+pb(i,k,j))/p1000mb)**cvpm) - alb(i,k,j)
  904. ENDDO
  905. DO k=ktf-1,kts,-1 ! remaining layers, integrate down
  906. DO i=its,itf
  907. qtot = 0.
  908. DO ispe=PARAM_FIRST_SCALAR,n_moist
  909. qtot = qtot + 0.5*( moist(i,k ,j,ispe) + moist(i,k+1,j,ispe) )
  910. ENDDO
  911. qf2 = 1.
  912. qf1 = qtot*qf2
  913. p(i,k,j) = p(i,k+1,j) - (mu(i,j) + qf1*muts(i,j))/qf2/rdn(k+1)
  914. qvf = 1.+rvovrd*moist(i,k,j,P_QV)
  915. al(i,k,j) = (r_d/p1000mb)*(t(i,k,j)+t0)*qvf* &
  916. (((p(i,k,j)+pb(i,k,j))/p1000mb)**cvpm) - alb(i,k,j)
  917. ENDDO
  918. ENDDO
  919. ENDDO
  920. ELSE
  921. DO j=jts,jtf
  922. k=ktf ! top layer
  923. DO i=its,itf
  924. qtot = 0.
  925. qf2 = 1.
  926. qf1 = qtot*qf2
  927. p(i,k,j) = - 0.5*(mu(i,j)+qf1*muts(i,j))/rdnw(k)/qf2
  928. qvf = 1.
  929. al(i,k,j) = (r_d/p1000mb)*(t(i,k,j)+t0)*qvf* &
  930. (((p(i,k,j)+pb(i,k,j))/p1000mb)**cvpm) - alb(i,k,j)
  931. ENDDO
  932. DO k=ktf-1,kts,-1 ! remaining layers, integrate down
  933. DO i=its,itf
  934. qtot = 0.
  935. qf2 = 1.
  936. qf1 = qtot*qf2
  937. p(i,k,j) = p(i,k+1,j) - (mu(i,j) + qf1*muts(i,j))/qf2/rdn(k+1)
  938. qvf = 1.
  939. al(i,k,j) = (r_d/p1000mb)*(t(i,k,j)+t0)*qvf* &
  940. (((p(i,k,j)+pb(i,k,j))/p1000mb)**cvpm) - alb(i,k,j)
  941. ENDDO
  942. ENDDO
  943. ENDDO
  944. END IF
  945. IF (hypsometric_opt == 1) THEN
  946. DO j=jts,jtf
  947. DO k=2,ktf+1 ! integrate hydrostatic equation for geopotential
  948. DO i=its,itf
  949. ph(i,k,j) = ph(i,k-1,j) - (dnw(k-1))*( &
  950. (muts(i,j))*al(i,k-1,j)+ &
  951. mu(i,j)*alb(i,k-1,j) )
  952. ENDDO
  953. ENDDO
  954. ENDDO
  955. ELSE
  956. ! Revised hypsometric eq.: dZ=-al*p*dLOG(p), where p is dry pressure
  957. DO j=jts,jtf
  958. DO i=its,itf
  959. ph(i,kts,j) = phb(i,kts,j)
  960. END DO
  961. DO k=kts+1,ktf+1
  962. DO i=its,itf
  963. pfu = muts(i,j)*znw(k) +ptop
  964. pfd = muts(i,j)*znw(k-1)+ptop
  965. phm = muts(i,j)*znu(k-1)+ptop
  966. ph(i,k,j) = ph(i,k-1,j) + (al(i,k-1,j)+alb(i,k-1,j))*phm*LOG(pfd/pfu)
  967. ENDDO
  968. ENDDO
  969. DO k=kts,ktf+1
  970. DO i=its,itf
  971. ph(i,k,j) = ph(i,k,j) - phb(i,k,j)
  972. END DO
  973. END DO
  974. END DO
  975. END IF
  976. END IF
  977. END SUBROUTINE calc_p_rho_phi
  978. !----------------------------------------------------------------------
  979. SUBROUTINE calc_php ( php, ph, phb, &
  980. ids, ide, jds, jde, kds, kde, &
  981. ims, ime, jms, jme, kms, kme, &
  982. its, ite, jts, jte, kts, kte )
  983. IMPLICIT NONE
  984. ! Input data
  985. INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
  986. ims, ime, jms, jme, kms, kme, &
  987. its, ite, jts, jte, kts, kte
  988. REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(IN ) :: phb, ph
  989. REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT( OUT) :: php
  990. ! Local stuff
  991. INTEGER :: i, j, k, itf, jtf, ktf
  992. !<DESCRIPTION>
  993. !
  994. ! calc_php calculates the full geopotential from the reference state
  995. ! geopotential and the perturbation geopotential (phb_ph).
  996. !
  997. !</DESCRIPTION>
  998. itf=MIN(ite,ide-1)
  999. jtf=MIN(jte,jde-1)
  1000. ktf=MIN(kte,kde-1)
  1001. DO j=jts,jtf
  1002. DO k=kts,ktf
  1003. DO i=its,itf
  1004. php(i,k,j) = 0.5*(phb(i,k,j)+phb(i,k+1,j)+ph(i,k,j)+ph(i,k+1,j))
  1005. ENDDO
  1006. ENDDO
  1007. ENDDO
  1008. END SUBROUTINE calc_php
  1009. !-------------------------------------------------------------------------------
  1010. SUBROUTINE diagnose_w( ph_tend, ph_new, ph_old, w, mu, dt, &
  1011. u, v, ht, &
  1012. cf1, cf2, cf3, rdx, rdy, &
  1013. msftx, msfty, &
  1014. ids, ide, jds, jde, kds, kde, &
  1015. ims, ime, jms, jme, kms, kme, &
  1016. its, ite, jts, jte, kts, kte )
  1017. IMPLICIT NONE
  1018. INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
  1019. ims, ime, jms, jme, kms, kme, &
  1020. its, ite, jts, jte, kts, kte
  1021. REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(IN ) :: ph_tend, &
  1022. ph_new, &
  1023. ph_old, &
  1024. u, &
  1025. v
  1026. REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT( OUT) :: w
  1027. REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN ) :: mu, ht, msftx, msfty
  1028. REAL, INTENT(IN ) :: dt, cf1, cf2, cf3, rdx, rdy
  1029. INTEGER :: i, j, k, itf, jtf
  1030. itf=MIN(ite,ide-1)
  1031. jtf=MIN(jte,jde-1)
  1032. !<DESCRIPTION>
  1033. !
  1034. ! diagnose_w diagnoses the vertical velocity from the geopoential equation.
  1035. ! Used with the hydrostatic option.
  1036. !
  1037. !</DESCRIPTION>
  1038. DO j = jts, jtf
  1039. ! lower b.c. on w
  1040. ! Notes on map scale factors:
  1041. ! Chain rule: if Z=Z(X,Y) [true at the surface] then
  1042. ! dZ/dt = dZ/dX * dX/dt + dZ/dY * dY/dt, U=dX/dt, V=dY/dt
  1043. ! Using capitals to denote actual values
  1044. ! In mapped values, u=U, v=V, z=Z, 1/dX=mx/dx, 1/dY=my/dy
  1045. ! => w = dz/dt = mx u dz/dx + my v dz/dy
  1046. ! [where dz/dx is just the surface height change between x
  1047. ! gridpoints, and dz/dy is the change between y gridpoints]
  1048. ! [NB: cf1, cf2 and cf3 do vertical weighting of u or v values
  1049. ! nearest the surface]
  1050. ! Previously msft multiplied by rdy and rdx terms.
  1051. ! Now msfty multiplies rdy term, and msftx multiplies msftx term
  1052. DO i = its, itf
  1053. w(i,1,j)= msfty(i,j)*.5*rdy*( &
  1054. (ht(i,j+1)-ht(i,j )) &
  1055. *(cf1*v(i,1,j+1)+cf2*v(i,2,j+1)+cf3*v(i,3,j+1)) &
  1056. +(ht(i,j )-ht(i,j-1)) &
  1057. *(cf1*v(i,1,j )+cf2*v(i,2,j )+cf3*v(i,3,j )) ) &
  1058. +msftx(i,j)*.5*rdx*( &
  1059. (ht(i+1,j)-ht(i,j )) &
  1060. *(cf1*u(i+1,1,j)+cf2*u(i+1,2,j)+cf3*u(i+1,3,j)) &
  1061. +(ht(i,j )-ht(i-1,j)) &
  1062. *(cf1*u(i ,1,j)+cf2*u(i ,2,j)+cf3*u(i ,3,j)) )
  1063. ENDDO
  1064. ! use geopotential equation to diagnose w
  1065. ! Further notes on map scale factors
  1066. ! If ph_tend contains: -mx partial d/dx(mu rho u/my)
  1067. ! -mx partial d/dy(phi mu v/mx)
  1068. ! -partial d/dz(phi mu w/my)
  1069. ! then phi eqn is: partial d/dt(mu phi/my) = ph_tend + mu g w/my
  1070. ! => w = [my/(mu*g)]*[partial d/dt(mu phi/my) - ph_tend]
  1071. DO k = 2, kte
  1072. DO i = its, itf
  1073. w(i,k,j) = msfty(i,j)*( (ph_new(i,k,j)-ph_old(i,k,j))/dt &
  1074. - ph_tend(i,k,j)/mu(i,j) )/g
  1075. ENDDO
  1076. ENDDO
  1077. ENDDO
  1078. END SUBROUTINE diagnose_w
  1079. !-------------------------------------------------------------------------------
  1080. SUBROUTINE rhs_ph( ph_tend, u, v, ww, &
  1081. ph, ph_old, phb, w, &
  1082. mut, muu, muv, &
  1083. fnm, fnp, &
  1084. rdnw, cfn, cfn1, rdx, rdy, &
  1085. msfux, msfuy, msfvx, &
  1086. msfvx_inv, msfvy, &
  1087. msftx, msfty, &
  1088. non_hydrostatic, &
  1089. config_flags, &
  1090. ids, ide, jds, jde, kds, kde, &
  1091. ims, ime, jms, jme, kms, kme, &
  1092. its, ite, jts, jte, kts, kte )
  1093. IMPLICIT NONE
  1094. TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
  1095. INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
  1096. ims, ime, jms, jme, kms, kme, &
  1097. its, ite, jts, jte, kts, kte
  1098. REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(IN ) :: &
  1099. u, &
  1100. v, &
  1101. ww, &
  1102. ph, &
  1103. ph_old, &
  1104. phb, &
  1105. w
  1106. REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(INOUT) :: ph_tend
  1107. REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN ) :: muu, muv, mut, &
  1108. msfux, msfuy, &
  1109. msfvx, msfvy, &
  1110. msftx, msfty, &
  1111. msfvx_inv
  1112. REAL, DIMENSION( kms:kme ), INTENT(IN ) :: rdnw, fnm, fnp
  1113. REAL, INTENT(IN ) :: cfn, cfn1, rdx, rdy
  1114. LOGICAL, INTENT(IN ) :: non_hydrostatic
  1115. ! Local stuff
  1116. INTEGER :: i, j, k, itf, jtf, ktf, kz, i_start, j_start
  1117. REAL :: ur, ul, ub, vr, vl, vb
  1118. REAL, DIMENSION(its:ite,kts:kte) :: wdwn
  1119. INTEGER :: advective_order
  1120. LOGICAL :: specified
  1121. !<DESCRIPTION>
  1122. !
  1123. ! rhs_ph calculates the large-timestep tendency terms for the geopotential
  1124. ! equation. These terms include the advection and "gw". The geopotential
  1125. ! equation is cast in advective form, so we don't use the flux form advection
  1126. ! algorithms here.
  1127. !
  1128. !</DESCRIPTION>
  1129. specified = .false.
  1130. if(config_flags%specified .or. config_flags%nested) specified = .true.
  1131. advective_order = config_flags%h_sca_adv_order
  1132. itf=MIN(ite,ide-1)
  1133. jtf=MIN(jte,jde-1)
  1134. ktf=MIN(kte,kde-1)
  1135. ! Notes on map scale factors (WCS, 2 march 2008)
  1136. ! phi equation is: mu/my d/dt(phi) = -(1/my) mx mu u d/dx(phi)
  1137. ! -(1/my) my mu v d/dy(phi)
  1138. ! - omega d/d_eta(phi)
  1139. ! +mu g w/my
  1140. !
  1141. ! A little further explanation...
  1142. ! The tendency term we are computing here is for mu/my d/dt(phi). It is advective form
  1143. ! but it is multiplied be mu/my. It will be decoupled from (mu/my) when the implicit w-phi
  1144. ! solution is computed in subourine advance_w. The formulation dates from the early
  1145. ! days of the mass coordinate model when we were testing both a flux and an advective formulation
  1146. ! for the geopotential equation and different forms of the vertical momentum equation and the
  1147. ! vertically implicit solver.
  1148. ! advective form for the geopotential equation
  1149. DO j = jts, jtf
  1150. DO k = 2, kte
  1151. DO i = its, itf
  1152. wdwn(i,k) = .5*(ww(i,k,j)+ww(i,k-1,j))*rdnw(k-1) &
  1153. *(ph(i,k,j)-ph(i,k-1,j)+phb(i,k,j)-phb(i,k-1,j))
  1154. ENDDO
  1155. ENDDO
  1156. ! RHS term 3 is: - omega partial d/dnu(phi)
  1157. DO k = 2, kte-1
  1158. DO i = its, itf
  1159. ph_tend(i,k,j) = ph_tend(i,k,j) &
  1160. - (fnm(k)*wdwn(i,k+1)+fnp(k)*wdwn(i,k))
  1161. ENDDO
  1162. ENDDO
  1163. ENDDO
  1164. IF (non_hydrostatic) THEN ! add in "gw" term.
  1165. DO j = jts, jtf ! in hydrostatic mode, "gw" will be diagnosed
  1166. ! after the timestep to give us "w"
  1167. DO i = its, itf
  1168. ph_tend(i,kde,j) = 0.
  1169. ENDDO
  1170. DO k = 2, kte
  1171. DO i = its, itf
  1172. ! phi equation RHS term 4
  1173. ph_tend(i,k,j) = ph_tend(i,k,j) + mut(i,j)*g*w(i,k,j)/msfty(i,j)
  1174. ENDDO
  1175. ENDDO
  1176. ENDDO
  1177. END IF
  1178. ! Notes on map scale factors:
  1179. ! RHS terms 1 and 2 are: -(1/my) mx u mu partial d/dx(phi)
  1180. ! -(1/my) my v mu partial d/dy(phi)
  1181. IF (advective_order <= 2) THEN
  1182. ! y (v) advection
  1183. i_start = its
  1184. j_start = jts
  1185. itf=MIN(ite,ide-1)
  1186. jtf=MIN(jte,jde-1)
  1187. IF ( (config_flags%open_ys .or. specified) .and. jts == jds ) j_start = jts+1
  1188. IF ( (config_flags%open_ye .or. specified) .and. jte == jde ) jtf = jtf-2
  1189. DO j = j_start, jtf
  1190. DO k = 2, kte-1
  1191. DO i = i_start, itf
  1192. ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdy/msfty(i,j))* &
  1193. ( muv(i,j+1)*(v(i,k,j+1)+v(i,k-1,j+1))*msfvy(i,j+1)* &
  1194. (phb(i,k,j+1)-phb(i,k,j )+ph(i,k,j+1)-ph(i,k,j )) &
  1195. +muv(i,j )*(v(i,k,j )+v(i,k-1,j ))*msfvy(i,j )* &
  1196. (phb(i,k,j )-phb(i,k,j-1)+ph(i,k,j )-ph(i,k,j-1)) )
  1197. ENDDO
  1198. ENDDO
  1199. k = kte
  1200. DO i = i_start, itf
  1201. ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdy/msfty(i,j))* &
  1202. ( muv(i,j+1)*(cfn*v(i,k-1,j+1)+cfn1*v(i,k-2,j+1))*msfvy(i,j+1)* &
  1203. (phb(i,k,j+1)-phb(i,k,j )+ph(i,k,j+1)-ph(i,k,j )) &
  1204. +muv(i,j )*(cfn*v(i,k-1,j )+cfn1*v(i,k-2,j ))*msfvy(i,j )* &
  1205. (phb(i,k,j )-phb(i,k,j-1)+ph(i,k,j )-ph(i,k,j-1)) )
  1206. ENDDO
  1207. ENDDO
  1208. ! x (u) advection
  1209. i_start = its
  1210. j_start = jts
  1211. itf=MIN(ite,ide-1)
  1212. jtf=MIN(jte,jde-1)
  1213. IF ( (config_flags%open_xs .or. specified) .and. its == ids ) i_start = its+1
  1214. IF ( (config_flags%open_xe .or. specified) .and. ite == ide ) itf = itf-2
  1215. DO j = j_start, jtf
  1216. DO k = 2, kte-1
  1217. DO i = i_start, itf
  1218. ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdx/msfty(i,j))* &
  1219. ( muu(i+1,j)*(u(i+1,k,j)+u(i+1,k-1,j))*msfux(i+1,j)* &
  1220. (phb(i+1,k,j)-phb(i ,k,j)+ph(i+1,k,j)-ph(i ,k,j)) &
  1221. +muu(i ,j)*(u(i ,k,j)+u(i ,k-1,j))*msfux(i ,j)* &
  1222. (phb(i ,k,j)-phb(i-1,k,j)+ph(i ,k,j)-ph(i-1,k,j)) )
  1223. ENDDO
  1224. ENDDO
  1225. k = kte
  1226. DO i = i_start, itf
  1227. ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdx/msfty(i,j))* &
  1228. ( muu(i+1,j)*(cfn*u(i+1,k-1,j)+cfn1*u(i+1,k-2,j))*msfux(i+1,j)* &
  1229. (phb(i+1,k,j)-phb(i ,k,j)+ph(i+1,k,j)-ph(i ,k,j)) &
  1230. +muu(i ,j)*(cfn*u(i ,k-1,j)+cfn1*u(i ,k-2,j))*msfux( i,j)* &
  1231. (phb(i ,k,j)-phb(i-1,k,j)+ph(i ,k,j)-ph(i-1,k,j)) )
  1232. ENDDO
  1233. ENDDO
  1234. ELSE IF (advective_order <= 4) THEN
  1235. ! y (v) advection
  1236. i_start = its
  1237. j_start = jts
  1238. itf=MIN(ite,ide-1)
  1239. jtf=MIN(jte,jde-1)
  1240. IF ( (config_flags%open_ys .or. specified) .and. jts == jds ) j_start = jts+2
  1241. IF ( (config_flags%open_ye .or. specified) .and. jte == jde ) jtf = jtf-3
  1242. DO j = j_start, jtf
  1243. DO k = 2, kte-1
  1244. DO i = i_start, itf
  1245. ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdy/msfty(i,j))*( &
  1246. ( muv(i,j+1)*(v(i,k,j+1)+v(i,k-1,j+1))*msfvy(i,j+1) &
  1247. +muv(i,j )*(v(i,k,j )+v(i,k-1,j ))*msfvy(i,j ))* (1./12.)*( &
  1248. 8.*(ph(i,k,j+1)-ph(i,k,j-1)) &
  1249. -(ph(i,k,j+2)-ph(i,k,j-2)) &
  1250. +8.*(phb(i,k,j+1)-phb(i,k,j-1)) &
  1251. -(phb(i,k,j+2)-phb(i,k,j-2)) ) )
  1252. ENDDO
  1253. ENDDO
  1254. k = kte
  1255. DO i = i_start, itf
  1256. ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdy/msfty(i,j))*( &
  1257. ( muv(i,j+1)*(cfn*v(i,k-1,j+1)+cfn1*v(i,k-2,j+1))*msfvy(i,j+1) &
  1258. +muv(i,j )*(cfn*v(i,k-1,j )+cfn1*v(i,k-2,j ))*msfvy(i,j ))* (1./12.)*( &
  1259. 8.*(ph(i,k,j+1)-ph(i,k,j-1)) &
  1260. -(ph(i,k,j+2)-ph(i,k,j-2)) &
  1261. +8.*(phb(i,k,j+1)-phb(i,k,j-1)) &
  1262. -(phb(i,k,j+2)-phb(i,k,j-2)) ) )
  1263. ENDDO
  1264. ENDDO
  1265. ! pick up near boundary rows using 2nd order stencil for open and specified conditions
  1266. IF ( (config_flags%open_ys .or. specified) .and. jts <= jds+1 ) THEN
  1267. j = jds+1
  1268. DO k = 2, kte-1
  1269. DO i = i_start, itf
  1270. ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdy/msfty(i,j))* &
  1271. ( muv(i,j+1)*(v(i,k,j+1)+v(i,k-1,j+1))*msfvy(i,j+1)* &
  1272. (phb(i,k,j+1)-phb(i,k,j )+ph(i,k,j+1)-ph(i,k,j )) &
  1273. +muv(i,j )*(v(i,k,j )+v(i,k-1,j ))*msfvy(i,j )* &
  1274. (phb(i,k,j )-phb(i,k,j-1)+ph(i,k,j )-ph(i,k,j-1)) )
  1275. ENDDO
  1276. ENDDO
  1277. k = kte
  1278. DO i = i_start, itf
  1279. ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdy/msfty(i,j))* &
  1280. ( muv(i,j+1)*(cfn*v(i,k-1,j+1)+cfn1*v(i,k-2,j+1))*msfvy(i,j+1)* &
  1281. (phb(i,k,j+1)-phb(i,k,j )+ph(i,k,j+1)-ph(i,k,j )) &
  1282. +muv(i,j )*(cfn*v(i,k-1,j )+cfn1*v(i,k-2,j ))*msfvy(i,j )* &
  1283. (phb(i,k,j )-phb(i,k,j-1)+ph(i,k,j )-ph(i,k,j-1)) )
  1284. ENDDO
  1285. END IF
  1286. IF ( (config_flags%open_ye .or. specified) .and. jte >= jde-2 ) THEN
  1287. j = jde-2
  1288. DO k = 2, kte-1
  1289. DO i = i_start, itf
  1290. ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdy/msfty(i,j))* &
  1291. ( muv(i,j+1)*(v(i,k,j+1)+v(i,k-1,j+1))*msfvy(i,j+1)* &
  1292. (phb(i,k,j+1)-phb(i,k,j )+ph(i,k,j+1)-ph(i,k,j )) &
  1293. +muv(i,j )*(v(i,k,j )+v(i,k-1,j ))*msfvy(i,j )* &
  1294. (phb(i,k,j )-phb(i,k,j-1)+ph(i,k,j )-ph(i,k,j-1)) )
  1295. ENDDO
  1296. ENDDO
  1297. k = kte
  1298. DO i = i_start, itf
  1299. ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdy/msfty(i,j))* &
  1300. ( muv(i,j+1)*(cfn*v(i,k-1,j+1)+cfn1*v(i,k-2,j+1))*msfvy(i,j+1)* &
  1301. (phb(i,k,j+1)-phb(i,k,j )+ph(i,k,j+1)-ph(i,k,j )) &
  1302. +muv(i,j )*(cfn*v(i,k-1,j )+cfn1*v(i,k-2,j ))*msfvy(i,j )* &
  1303. (phb(i,k,j )-phb(i,k,j-1)+ph(i,k,j )-ph(i,k,j-1)) )
  1304. ENDDO
  1305. END IF
  1306. ! x (u) advection
  1307. i_start = its
  1308. j_start = jts
  1309. itf=MIN(ite,ide-1)
  1310. jtf=MIN(jte,jde-1)
  1311. IF ( (config_flags%open_xs) .and. its == ids ) i_start = its+2
  1312. IF ( (config_flags%open_xe) .and. ite == ide ) itf = itf-3
  1313. DO j = j_start, jtf
  1314. DO k = 2, kte-1
  1315. DO i = i_start, itf
  1316. ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdx/msfty(i,j))*( &
  1317. ( muu(i+1,j)*(u(i+1,k,j)+u(i+1,k-1,j))*msfux(i+1,j) &
  1318. +muu(i,j )*(u(i ,k,j)+u(i ,k-1,j))*msfux(i ,j) )* (1./12.)*( &
  1319. 8.*(ph(i+1,k,j)-ph(i-1,k,j)) &
  1320. -(ph(i+2,k,j)-ph(i-2,k,j)) &
  1321. +8.*(phb(i+1,k,j)-phb(i-1,k,j)) &
  1322. -(phb(i+2,k,j)-phb(i-2,k,j)) ) )
  1323. ENDDO
  1324. ENDDO
  1325. k = kte
  1326. DO i = i_start, itf
  1327. ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdx/msfty(i,j))*( &
  1328. ( muu(i+1,j)*(cfn*u(i+1,k-1,j)+cfn1*u(i+1,k-2,j))*msfux(i+1,j) &
  1329. +muu(i,j )*(cfn*u(i ,k-1,j)+cfn1*u(i ,k-2,j))*msfux(i ,j) )* (1./12.)*( &
  1330. 8.*(ph(i+1,k,j)-ph(i-1,k,j)) &
  1331. -(ph(i+2,k,j)-ph(i-2,k,j)) &
  1332. +8.*(phb(i+1,k,j)-phb(i-1,k,j)) &
  1333. -(phb(i+2,k,j)-phb(i-2,k,j)) ) )
  1334. ENDDO
  1335. ENDDO
  1336. ! pick up near boundary rows using 2nd order stencil for open and specified conditions
  1337. IF ( (config_flags%open_xs .or. specified) .and. its <= ids+1 ) THEN
  1338. i = ids + 1
  1339. DO j = j_start, jtf
  1340. DO k = 2, kte-1
  1341. ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdx/msfty(i,j))* &
  1342. ( muu(i+1,j)*(u(i+1,k,j)+u(i+1,k-1,j))*msfux(i+1,j)* &
  1343. (phb(i+1,k,j)-phb(i ,k,j)+ph(i+1,k,j)-ph(i ,k,j)) &
  1344. +muu(i ,j)*(u(i ,k,j)+u(i ,k-1,j))*msfux(i ,j)* &
  1345. (phb(i ,k,j)-phb(i-1,k,j)+ph(i ,k,j)-ph(i-1,k,j)) )
  1346. ENDDO
  1347. ENDDO
  1348. k = kte
  1349. DO j = j_start, jtf
  1350. ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdx/msfty(i,j))* &
  1351. ( muu(i+1,j)*(cfn*u(i+1,k-1,j)+cfn1*u(i+1,k-2,j))*msfux(i+1,j)* &
  1352. (phb(i+1,k,j)-phb(i ,k,j)+ph(i+1,k,j)-ph(i ,k,j)) &
  1353. +muu(i ,j)*(cfn*u(i ,k-1,j)+cfn1*u(i ,k-2,j))*msfux( i,j)* &
  1354. (phb(i ,k,j)-phb(i-1,k,j)+ph(i ,k,j)-ph(i-1,k,j)) )
  1355. ENDDO
  1356. END IF
  1357. IF ( (config_flags%open_xe .or. specified) .and. ite >= ide-2 ) THEN
  1358. i = ide-2
  1359. DO j = j_start, jtf
  1360. DO k = 2, kte-1
  1361. ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdx/msfty(i,j))* &
  1362. ( muu(i+1,j)*(u(i+1,k,j)+u(i+1,k-1,j))*msfux(i+1,j)* &
  1363. (phb(i+1,k,j)-phb(i ,k,j)+ph(i+1,k,j)-ph(i ,k,j)) &
  1364. +muu(i ,j)*(u(i ,k,j)+u(i ,k-1,j))*msfux(i ,j)* &
  1365. (phb(i ,k,j)-phb(i-1,k,j)+ph(i ,k,j)-ph(i-1,k,j)) )
  1366. ENDDO
  1367. ENDDO
  1368. k = kte
  1369. DO j = j_start, jtf
  1370. ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdx/msfty(i,j))* &
  1371. ( muu(i+1,j)*(cfn*u(i+1,k-1,j)+cfn1*u(i+1,k-2,j))*msfux(i+1,j)* &
  1372. (phb(i+1,k,j)-phb(i ,k,j)+ph(i+1,k,j)-ph(i ,k,j)) &
  1373. +muu(i ,j)*(cfn*u(i ,k-1,j)+cfn1*u(i ,k-2,j))*msfux( i,j)* &
  1374. (phb(i ,k,j)-phb(i-1,k,j)+ph(i ,k,j)-ph(i-1,k,j)) )
  1375. ENDDO
  1376. END IF
  1377. !--------------------------
  1378. ELSE IF (advective_order <= 6) THEN
  1379. !!! NOTE: this branch has been looked at and fixed with changes for overdecomposition
  1380. !!! the branches covering the other advective_order cases have not. 20090923. JM
  1381. ! y (v) advection
  1382. i_start = its
  1383. j_start = jts
  1384. itf=MIN(ite,ide-1)
  1385. jtf=MIN(jte,jde-1)
  1386. IF (config_flags%open_ys .or. specified ) j_start = max(jts,jds+3)
  1387. IF (config_flags%open_ye .or. specified ) jtf = min(jtf,jde-4)
  1388. DO j = j_start, jtf
  1389. DO k = 2, kte-1
  1390. DO i = i_start, itf
  1391. ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdy/msfty(i,j))* ( &
  1392. ( muv(i,j+1)*(v(i,k,j+1)+v(i,k-1,j+1))*msfvy(i,j+1) &
  1393. +muv(i,j )*(v(i,k,j )+v(i,k-1,j ))*msfvy(i,j ) )* (1./60.)*( &
  1394. 45.*(ph(i,k,j+1)-ph(i,k,j-1)) &
  1395. -9.*(ph(i,k,j+2)-ph(i,k,j-2)) &
  1396. +(ph(i,k,j+3)-ph(i,k,j-3)) &
  1397. +45.*(phb(i,k,j+1)-phb(i,k,j-1)) &
  1398. -9.*(phb(i,k,j+2)-phb(i,k,j-2)) &
  1399. +(phb(i,k,j+3)-phb(i,k,j-3)) ) )
  1400. ENDDO
  1401. ENDDO
  1402. k = kte
  1403. DO i = i_start, itf
  1404. ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdy/msfty(i,j))* ( &
  1405. ( muv(i,j+1)*(cfn*v(i,k-1,j+1)+cfn1*v(i,k-2,j+1))*msfvy(i,j+1) &
  1406. +muv(i,j )*(cfn*v(i,k-1,j )+cfn1*v(i,k-2,j ))*msfvy(i,j ) )* (1./60.)*( &
  1407. 45.*(ph(i,k,j+1)-ph(i,k,j-1)) &
  1408. -9.*(ph(i,k,j+2)-ph(i,k,j-2)) &
  1409. +(ph(i,k,j+3)-ph(i,k,j-3)) &
  1410. +45.*(phb(i,k,j+1)-phb(i,k,j-1)) &
  1411. -9.*(phb(i,k,j+2)-phb(i,k,j-2)) &
  1412. +(phb(i,k,j+3)-phb(i,k,j-3)) ) )
  1413. ENDDO
  1414. ENDDO
  1415. ! 4th order stencil for open or specified conditions two in form the boundary
  1416. IF ( (config_flags%open_ys .or. specified) .and. jts <= jds+2 .and. jds+2 <= jte ) THEN
  1417. j = jds+2
  1418. DO k = 2, kte-1
  1419. DO i = i_start, itf
  1420. ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdy/msfty(i,j))* ( &
  1421. ( muv(i,j+1)*(v(i,k,j+1)+v(i,k-1,j+1))*msfvy(i,j+1) &
  1422. +muv(i,j )*(v(i,k,j )+v(i,k-1,j ))*msfvy(i,j ) )* (1./12.)*( &
  1423. 8.*(ph(i,k,j+1)-ph(i,k,j-1)) &
  1424. -(ph(i,k,j+2)-ph(i,k,j-2)) &
  1425. +8.*(phb(i,k,j+1)-phb(i,k,j-1)) &
  1426. -(phb(i,k,j+2)-phb(i,k,j-2)) ) )
  1427. ENDDO
  1428. ENDDO
  1429. k = kte
  1430. DO i = i_start, itf
  1431. ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdy/msfty(i,j))* ( &
  1432. ( muv(i,j+1)*(cfn*v(i,k-1,j+1)+cfn1*v(i,k-2,j+1))*msfvy(i,j+1) &
  1433. +muv(i,j )*(cfn*v(i,k-1,j )+cfn1*v(i,k-2,j ))*msfvy(i,j) )* (1./12.)*( &
  1434. 8.*(ph(i,k,j+1)-ph(i,k,j-1)) &
  1435. -(ph(i,k,j+2)-ph(i,k,j-2)) &
  1436. +8.*(phb(i,k,j+1)-phb(i,k,j-1)) &
  1437. -(phb(i,k,j+2)-phb(i,k,j-2)) ) )
  1438. ENDDO
  1439. END IF
  1440. IF ( (config_flags%open_ye .or. specified) .and. jts <= jde-3 .and. jde-3 <= jte ) THEN
  1441. j = jde-3
  1442. DO k = 2, kte-1
  1443. DO i = i_start, itf
  1444. ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdy/msfty(i,j))* ( &
  1445. ( muv(i,j+1)*(v(i,k,j+1)+v(i,k-1,j+1))*msfvy(i,j+1) &
  1446. +muv(i,j )*(v(i,k,j )+v(i,k-1,j ))*msfvy(i,j) )* (1./12.)*( &
  1447. 8.*(ph(i,k,j+1)-ph(i,k,j-1)) &
  1448. -(ph(i,k,j+2)-ph(i,k,j-2)) &
  1449. +8.*(phb(i,k,j+1)-phb(i,k,j-1)) &
  1450. -(phb(i,k,j+2)-phb(i,k,j-2)) ) )
  1451. ENDDO
  1452. ENDDO
  1453. k = kte
  1454. DO i = i_start, itf
  1455. ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdy/msfty(i,j))* ( &
  1456. ( muv(i,j+1)*(cfn*v(i,k-1,j+1)+cfn1*v(i,k-2,j+1))*msfvy(i,j+1) &
  1457. +muv(i,j )*(cfn*v(i,k-1,j )+cfn1*v(i,k-2,j ))*msfvy(i,j) )* (1./12.)*( &
  1458. 8.*(ph(i,k,j+1)-ph(i,k,j-1)) &
  1459. -(ph(i,k,j+2)-ph(i,k,j-2)) &
  1460. +8.*(phb(i,k,j+1)-phb(i,k,j-1)) &
  1461. -(phb(i,k,j+2)-phb(i,k,j-2)) ) )
  1462. ENDDO
  1463. END IF
  1464. ! 2nd order stencil for open and specified conditions one row in from boundary
  1465. IF ( (config_flags%open_ys .or. specified) .and. jts <= jds+1 .and. jds+1 <= jte ) THEN
  1466. j = jds+1
  1467. DO k = 2, kte-1
  1468. DO i = i_start, itf
  1469. ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdy/msfty(i,j))* &
  1470. ( muv(i,j+1)*(v(i,k,j+1)+v(i,k-1,j+1))*msfvy(i,j+1)* &
  1471. (phb(i,k,j+1)-phb(i,k,j )+ph(i,k,j+1)-ph(i,k,j )) &
  1472. +muv(i,j )*(v(i,k,j )+v(i,k-1,j ))*msfvy(i,j )* &
  1473. (phb(i,k,j )-phb(i,k,j-1)+ph(i,k,j )-ph(i,k,j-1)) )
  1474. ENDDO
  1475. ENDDO
  1476. k = kte
  1477. DO i = i_start, itf
  1478. ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdy/msfty(i,j))* &
  1479. ( muv(i,j+1)*(cfn*v(i,k-1,j+1)+cfn1*v(i,k-2,j+1))*msfvy(i,j+1)* &
  1480. (phb(i,k,j+1)-phb(i,k,j )+ph(i,k,j+1)-ph(i,k,j )) &
  1481. +muv(i,j )*(cfn*v(i,k-1,j )+cfn1*v(i,k-2,j ))*msfvy(i,j )* &
  1482. (phb(i,k,j )-phb(i,k,j-1)+ph(i,k,j )-ph(i,k,j-1)) )
  1483. ENDDO
  1484. END IF
  1485. IF ( (config_flags%open_ye .or. specified) .and. jts <= jde-2 .and. jde-2 <= jte ) THEN
  1486. j = jde-2
  1487. DO k = 2, kte-1
  1488. DO i = i_start, itf
  1489. ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdy/msfty(i,j))* &
  1490. ( muv(i,j+1)*(v(i,k,j+1)+v(i,k-1,j+1))*msfvy(i,j+1)* &
  1491. (phb(i,k,j+1)-phb(i,k,j )+ph(i,k,j+1)-ph(i,k,j )) &
  1492. +muv(i,j )*(v(i,k,j )+v(i,k-1,j ))*msfvy(i,j )* &
  1493. (phb(i,k,j )-phb(i,k,j-1)+ph(i,k,j )-ph(i,k,j-1)) )
  1494. ENDDO
  1495. ENDDO
  1496. k = kte
  1497. DO i = i_start, itf
  1498. ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdy/msfty(i,j))* &
  1499. ( muv(i,j+1)*(cfn*v(i,k-1,j+1)+cfn1*v(i,k-2,j+1))*msfvy(i,j+1)* &
  1500. (phb(i,k,j+1)-phb(i,k,j )+ph(i,k,j+1)-ph(i,k,j )) &
  1501. +muv(i,j )*(cfn*v(i,k-1,j )+cfn1*v(i,k-2,j ))*msfvy(i,j )* &
  1502. (phb(i,k,j )-phb(i,k,j-1)+ph(i,k,j )-ph(i,k,j-1)) )
  1503. ENDDO
  1504. END IF
  1505. ! x (u) advection
  1506. i_start = its
  1507. j_start = jts
  1508. itf=MIN(ite,ide-1)
  1509. jtf=MIN(jte,jde-1)
  1510. IF (config_flags%open_xs .or. specified ) i_start = max(its,ids+3)
  1511. IF (config_flags%open_xe .or. specified ) itf = min(itf,ide-4)
  1512. DO j = j_start, jtf
  1513. DO k = 2, kte-1
  1514. DO i = i_start, itf
  1515. ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdx/msfty(i,j))*( &
  1516. ( muu(i+1,j)*(u(i+1,k,j)+u(i+1,k-1,j))*msfux(i+1,j) &
  1517. +muu(i,j )*(u(i,k,j )+u(i,k-1,j ))*msfux(i,j) )* (1./60.)*( &
  1518. 45.*(ph(i+1,k,j)-ph(i-1,k,j)) &
  1519. -9.*(ph(i+2,k,j)-ph(i-2,k,j)) &
  1520. +(ph(i+3,k,j)-ph(i-3,k,j)) &
  1521. +45.*(phb(i+1,k,j)-phb(i-1,k,j)) &
  1522. -9.*(phb(i+2,k,j)-phb(i-2,k,j)) &
  1523. +(phb(i+3,k,j)-phb(i-3,k,j)) ) )
  1524. ENDDO
  1525. ENDDO
  1526. k = kte
  1527. DO i = i_start, itf
  1528. ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdx/msfty(i,j))*( &
  1529. ( muu(i+1,j)*(cfn*u(i+1,k-1,j)+cfn1*u(i+1,k-2,j))*msfux(i+1,j) &
  1530. +muu(i,j )*(cfn*u(i ,k-1,j)+cfn1*u(i,k-2,j))*msfux(i,j) )* (1./60.)*( &
  1531. 45.*(ph(i+1,k,j)-ph(i-1,k,j)) &
  1532. -9.*(ph(i+2,k,j)-ph(i-2,k,j)) &
  1533. +(ph(i+3,k,j)-ph(i-3,k,j)) &
  1534. +45.*(phb(i+1,k,j)-phb(i-1,k,j)) &
  1535. -9.*(phb(i+2,k,j)-phb(i-2,k,j)) &
  1536. +(phb(i+3,k,j)-phb(i-3,k,j)) ) )
  1537. ENDDO
  1538. ENDDO
  1539. ! 4th order stencil two in from the boundary for open and specified conditions
  1540. IF ( (config_flags%open_xs) .and. its <= ids+2 .and. ids+2 <= ite ) THEN
  1541. i = ids + 2
  1542. DO j = j_start, jtf
  1543. DO k = 2, kte-1
  1544. ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdx/msfty(i,j))*( &
  1545. ( muu(i+1,j)*(u(i+1,k,j)+u(i+1,k-1,j))*msfux(i+1,j) &
  1546. +muu(i,j )*(u(i,k,j )+u(i,k-1,j ))*msfux(i,j) )* (1./12.)*( &
  1547. 8.*(ph(i+1,k,j)-ph(i-1,k,j)) &
  1548. -(ph(i+2,k,j)-ph(i-2,k,j)) &
  1549. +8.*(phb(i+1,k,j)-phb(i-1,k,j)) &
  1550. -(phb(i+2,k,j)-phb(i-2,k,j)) ) )
  1551. ENDDO
  1552. k = kte
  1553. ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdx/msfty(i,j))*( &
  1554. ( muu(i+1,j)*(cfn*u(i+1,k-1,j)+cfn1*u(i+1,k-2,j))*msfux(i+1,j) &
  1555. +muu(i,j )*(cfn*u(i ,k-1,j)+cfn1*u(i,k-2,j))*msfux(i,j) )* (1./12.)*( &
  1556. 8.*(ph(i+1,k,j)-ph(i-1,k,j)) &
  1557. -(ph(i+2,k,j)-ph(i-2,k,j)) &
  1558. +8.*(phb(i+1,k,j)-phb(i-1,k,j)) &
  1559. -(phb(i+2,k,j)-phb(i-2,k,j)) ) )
  1560. ENDDO
  1561. END IF
  1562. IF ( (config_flags%open_xe) .and. its <= ide-3 .and. ide-3 <= ite ) THEN
  1563. i = ide-3
  1564. DO j = j_start, jtf
  1565. DO k = 2, kte-1
  1566. ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdx/msfty(i,j))*( &
  1567. ( muu(i+1,j)*(u(i+1,k,j)+u(i+1,k-1,j))*msfux(i+1,j) &
  1568. +muu(i,j )*(u(i,k,j )+u(i,k-1,j ))*msfux(i,j) )* (1./12.)*( &
  1569. 8.*(ph(i+1,k,j)-ph(i-1,k,j)) &
  1570. -(ph(i+2,k,j)-ph(i-2,k,j)) &
  1571. +8.*(phb(i+1,k,j)-phb(i-1,k,j)) &
  1572. -(phb(i+2,k,j)-phb(i-2,k,j)) ) )
  1573. ENDDO
  1574. k = kte
  1575. ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdx/msfty(i,j))*( &
  1576. ( muu(i+1,j)*(cfn*u(i+1,k-1,j)+cfn1*u(i+1,k-2,j))*msfux(i+1,j) &
  1577. +muu(i,j )*(cfn*u(i ,k-1,j)+cfn1*u(i,k-2,j))*msfux(i,j) )* (1./12.)*( &
  1578. 8.*(ph(i+1,k,j)-ph(i-1,k,j)) &
  1579. -(ph(i+2,k,j)-ph(i-2,k,j)) &
  1580. +8.*(phb(i+1,k,j)-phb(i-1,k,j)) &
  1581. -(phb(i+2,k,j)-phb(i-2,k,j)) ) )
  1582. ENDDO
  1583. END IF
  1584. ! 2nd order stencil for open and specified conditions one in from the boundary
  1585. IF ( (config_flags%open_xs .or. specified) .and. its <= ids+1 .and. ids+1 <= ite ) THEN
  1586. i = ids + 1
  1587. DO j = j_start, jtf
  1588. DO k = 2, kte-1
  1589. ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdx/msfty(i,j))* &
  1590. ( muu(i+1,j)*(u(i+1,k,j)+u(i+1,k-1,j))*msfux(i+1,j)* &
  1591. (phb(i+1,k,j)-phb(i ,k,j)+ph(i+1,k,j)-ph(i ,k,j)) &
  1592. +muu(i ,j)*(u(i ,k,j)+u(i ,k-1,j))*msfux(i ,j)* &
  1593. (phb(i ,k,j)-phb(i-1,k,j)+ph(i ,k,j)-ph(i-1,k,j)) )
  1594. ENDDO
  1595. ENDDO
  1596. k = kte
  1597. DO j = j_start, jtf
  1598. ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdx/msfty(i,j))* &
  1599. ( muu(i+1,j)*(cfn*u(i+1,k-1,j)+cfn1*u(i+1,k-2,j))*msfux(i+1,j)* &
  1600. (phb(i+1,k,j)-phb(i ,k,j)+ph(i+1,k,j)-ph(i ,k,j)) &
  1601. +muu(i ,j)*(cfn*u(i ,k-1,j)+cfn1*u(i ,k-2,j))*msfux( i,j)* &
  1602. (phb(i ,k,j)-phb(i-1,k,j)+ph(i ,k,j)-ph(i-1,k,j)) )
  1603. ENDDO
  1604. END IF
  1605. IF ( (config_flags%open_xe .or. specified) .and. its <= ide-2 .and. ide-2 <= ite ) THEN
  1606. i = ide-2
  1607. DO j = j_start, jtf
  1608. DO k = 2, kte-1
  1609. ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdx/msfty(i,j))* &
  1610. ( muu(i+1,j)*(u(i+1,k,j)+u(i+1,k-1,j))*msfux(i+1,j)* &
  1611. (phb(i+1,k,j)-phb(i ,k,j)+ph(i+1,k,j)-ph(i ,k,j)) &
  1612. +muu(i ,j)*(u(i ,k,j)+u(i ,k-1,j))*msfux(i ,j)* &
  1613. (phb(i ,k,j)-phb(i-1,k,j)+ph(i ,k,j)-ph(i-1,k,j)) )
  1614. ENDDO
  1615. ENDDO
  1616. k = kte
  1617. DO j = j_start, jtf
  1618. ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdx/msfty(i,j))* &
  1619. ( muu(i+1,j)*(cfn*u(i+1,k-1,j)+cfn1*u(i+1,k-2,j))*msfux(i+1,j)* &
  1620. (phb(i+1,k,j)-phb(i ,k,j)+ph(i+1,k,j)-ph(i ,k,j)) &
  1621. +muu(i ,j)*(cfn*u(i ,k-1,j)+cfn1*u(i ,k-2,j))*msfux( i,j)* &
  1622. (phb(i ,k,j)-phb(i-1,k,j)+ph(i ,k,j)-ph(i-1,k,j)) )
  1623. ENDDO
  1624. END IF
  1625. END IF ! 6th order advection
  1626. ! lateral open boundary conditions,
  1627. ! start with north and south (y) boundaries
  1628. i_start = its
  1629. itf=MIN(ite,ide-1)
  1630. ! south
  1631. IF ( (config_flags%open_ys) .and. jts == jds ) THEN
  1632. j=jts
  1633. DO k=2,kde
  1634. kz = min(k,kde-1)
  1635. DO i = its,itf
  1636. vb =.5*( fnm(kz)*(v(i,kz ,j+1)+v(i,kz ,j )) &
  1637. +fnp(kz)*(v(i,kz-1,j+1)+v(i,kz-1,j )) )
  1638. vl=amin1(vb,0.)
  1639. ph_tend(i,k,j)=ph_tend(i,k,j)-rdy*mut(i,j)*( &
  1640. +vl*(ph_old(i,k,j+1)-ph_old(i,k,j)))
  1641. ENDDO
  1642. ENDDO
  1643. END IF
  1644. ! north
  1645. IF ( (config_flags%open_ye) .and. jte == jde ) THEN
  1646. j=jte-1
  1647. DO k=2,kde
  1648. kz = min(k,kde-1)
  1649. DO i = its,itf
  1650. vb=.5*( fnm(kz)*(v(i,kz ,j+1)+v(i,kz ,j)) &
  1651. +fnp(kz)*(v(i,kz-1,j+1)+v(i,kz-1,j)) )
  1652. vr=amax1(vb,0.)
  1653. ph_tend(i,k,j)=ph_tend(i,k,j)-rdy*mut(i,j)*( &
  1654. +vr*(ph_old(i,k,j)-ph_old(i,k,j-1)))
  1655. ENDDO
  1656. ENDDO
  1657. END IF
  1658. ! now the east and west (y) boundaries
  1659. j_start = its
  1660. jtf=MIN(jte,jde-1)
  1661. ! west
  1662. IF ( (config_flags%open_xs) .and. its == ids ) THEN
  1663. i=its
  1664. DO j = jts,jtf
  1665. DO k=2,kde-1
  1666. kz = k
  1667. ub =.5*( fnm(kz)*(u(i+1,kz ,j)+u(i ,kz ,j)) &
  1668. +fnp(kz)*(u(i+1,kz-1,j)+u(i ,kz-1,j)) )
  1669. ul=amin1(ub,0.)
  1670. ph_tend(i,k,j)=ph_tend(i,k,j)-(msftx(i,j)/msfty(i,j))*rdx*mut(i,j)*( &
  1671. +ul*(ph_old(i+1,k,j)-ph_old(i,k,j)))
  1672. ENDDO
  1673. k = kde
  1674. kz = k
  1675. ub =.5*( fnm(kz)*(u(i+1,kz ,j)+u(i ,kz ,j)) &
  1676. +fnp(kz)*(u(i+1,kz-1,j)+u(i ,kz-1,j)) )
  1677. ul=amin1(ub,0.)
  1678. ph_tend(i,k,j)=ph_tend(i,k,j)-(msftx(i,j)/msfty(i,j))*rdx*mut(i,j)*( &
  1679. +ul*(ph_old(i+1,k,j)-ph_old(i,k,j)))
  1680. ENDDO
  1681. END IF
  1682. ! east
  1683. IF ( (config_flags%open_xe) .and. ite == ide ) THEN
  1684. i = ite-1
  1685. DO j = jts,jtf
  1686. DO k=2,kde-1
  1687. kz = k
  1688. ub=.5*( fnm(kz)*(u(i+1,kz ,j)+u(i,kz ,j)) &
  1689. +fnp(kz)*(u(i+1,kz-1,j)+u(i,kz-1,j)) )
  1690. ur=amax1(ub,0.)
  1691. ph_tend(i,k,j)=ph_tend(i,k,j)-(msftx(i,j)/msfty(i,j))*rdx*mut(i,j)*( &
  1692. +ur*(ph_old(i,k,j)-ph_old(i-1,k,j)))
  1693. ENDDO
  1694. k = kde
  1695. kz = k-1
  1696. ub=.5*( fnm(kz)*(u(i+1,kz ,j)+u(i,kz ,j)) &
  1697. +fnp(kz)*(u(i+1,kz-1,j)+u(i,kz-1,j)) )
  1698. ur=amax1(ub,0.)
  1699. ph_tend(i,k,j)=ph_tend(i,k,j)-(msftx(i,j)/msfty(i,j))*rdx*mut(i,j)*( &
  1700. +ur*(ph_old(i,k,j)-ph_old(i-1,k,j)))
  1701. ENDDO
  1702. END IF
  1703. END SUBROUTINE rhs_ph
  1704. !-------------------------------------------------------------------------------
  1705. SUBROUTINE horizontal_pressure_gradient( ru_tend,rv_tend, &
  1706. ph,alt,p,pb,al,php,cqu,cqv, &
  1707. muu,muv,mu,fnm,fnp,rdnw, &
  1708. cf1,cf2,cf3,rdx,rdy,msfux,msfuy,&
  1709. msfvx,msfvy,msftx,msfty, &
  1710. config_flags, non_hydrostatic, &
  1711. top_lid, &
  1712. ids, ide, jds, jde, kds, kde, &
  1713. ims, ime, jms, jme, kms, kme, &
  1714. its, ite, jts, jte, kts, kte )
  1715. IMPLICIT NONE
  1716. ! Input data
  1717. TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
  1718. LOGICAL, INTENT (IN ) :: non_hydrostatic, top_lid
  1719. INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
  1720. ims, ime, jms, jme, kms, kme, &
  1721. its, ite, jts, jte, kts, kte
  1722. REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(IN ) :: &
  1723. ph, &
  1724. alt, &
  1725. al, &
  1726. p, &
  1727. pb, &
  1728. php, &
  1729. cqu, &
  1730. cqv
  1731. REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(INOUT) :: &
  1732. ru_tend, &
  1733. rv_tend
  1734. REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN ) :: muu, muv, mu, &
  1735. msfux, msfuy, &
  1736. msfvx, msfvy, &
  1737. msftx, msfty
  1738. REAL, DIMENSION( kms:kme ), INTENT(IN ) :: rdnw, fnm, fnp
  1739. REAL, INTENT(IN ) :: rdx, rdy, cf1, cf2, cf3
  1740. INTEGER :: i,j,k, itf, jtf, ktf, i_start, j_start
  1741. REAL, DIMENSION( ims:ime, kms:kme ) :: dpn
  1742. REAL :: dpx, dpy
  1743. LOGICAL :: specified
  1744. !<DESCRIPTION>
  1745. !
  1746. ! horizontal_pressure_gradient calculates the
  1747. ! horizontal pressure gradient terms for the large-timestep tendency
  1748. ! in the horizontal momentum equations (u,v).
  1749. !
  1750. !</DESCRIPTION>
  1751. specified = .false.
  1752. if(config_flags%specified .or. config_flags%nested) specified = .true.
  1753. ! Notes on map scale factors:
  1754. ! Calculates the pressure gradient terms in ADT eqns 44 and 45
  1755. ! With upper rho -> 'mu', these are:
  1756. ! Eqn 30: -mu*(mx/my)*(1/rho)*partial dp/dx
  1757. ! Eqn 31: -mu*(my/mx)*(1/rho)*partial dp/dy
  1758. !
  1759. ! As we are on nu, rather than height, surfaces:
  1760. !
  1761. ! mu dp/dx = mu alpha partial dp'/dx + (nu mu partial dmubar/dx) alpha'
  1762. ! + mu partial dphi'/dx + (partial dphi/dx)*(partial dp'/dnu - mu')
  1763. !
  1764. ! mu dp/dy = mu alpha partial dp'/dy + (nu mu partial dmubar/dy) alpha'
  1765. ! + mu partial dphi'/dy + (partial dphi/dy)*(partial dp'/dnu - mu')
  1766. ! start with the north-south (y) pressure gradient
  1767. itf=MIN(ite,ide-1)
  1768. jtf=jte
  1769. ktf=MIN(kte,kde-1)
  1770. i_start = its
  1771. j_start = jts
  1772. IF ( (config_flags%open_ys .or. specified .or. &
  1773. config_flags%nested .or. config_flags%polar ) .and. jts == jds ) j_start = jts+1
  1774. IF ( (config_flags%open_ye .or. specified .or. &
  1775. config_flags%nested .or. config_flags%polar ) .and. jte == jde ) jtf = jtf-1
  1776. DO j = j_start, jtf
  1777. IF ( non_hydrostatic ) THEN
  1778. k=1
  1779. DO i = i_start, itf
  1780. dpn(i,k) = .5*( cf1*(p(i,k ,j-1)+p(i,k ,j)) &
  1781. +cf2*(p(i,k+1,j-1)+p(i,k+1,j)) &
  1782. +cf3*(p(i,k+2,j-1)+p(i,k+2,j)) )
  1783. dpn(i,kde) = 0.
  1784. ENDDO
  1785. IF (top_lid) THEN
  1786. DO i = i_start, itf
  1787. dpn(i,kde) = .5*( cf1*(p(i,kde-1,j-1)+p(i,kde-1,j)) &
  1788. +cf2*(p(i,kde-2,j-1)+p(i,kde-2,j)) &
  1789. +cf3*(p(i,kde-3,j-1)+p(i,kde-3,j)) )
  1790. ENDDO
  1791. ENDIF
  1792. DO k=2,ktf
  1793. DO i = i_start, itf
  1794. dpn(i,k) = .5*( fnm(k)*(p(i,k ,j-1)+p(i,k ,j)) &
  1795. +fnp(k)*(p(i,k-1,j-1)+p(i,k-1,j)) )
  1796. END DO
  1797. END DO
  1798. ! ADT eqn 45: -mu*(my/mx)*(1/rho)*partial dp/dy
  1799. ! [alt, al are 1/rho terms; muv, mu are NOT coupled]
  1800. DO K=1,ktf
  1801. DO i = i_start, itf
  1802. ! Here are mu dp/dy terms 1-3
  1803. dpy = (msfvy(i,j)/msfvx(i,j))*.5*rdy*muv(i,j)*( &
  1804. (ph (i,k+1,j)-ph (i,k+1,j-1) + ph(i,k,j)-ph(i,k,j-1)) &
  1805. +(alt(i,k ,j)+alt(i,k ,j-1))*(p (i,k,j)-p (i,k,j-1)) &
  1806. +(al (i,k ,j)+al (i,k ,j-1))*(pb(i,k,j)-pb(i,k,j-1)) )
  1807. ! Here is mu dp/dy term 4
  1808. dpy = dpy + (msfvy(i,j)/msfvx(i,j))*rdy*(php(i,k,j)-php(i,k,j-1))* &
  1809. (rdnw(k)*(dpn(i,k+1)-dpn(i,k))-.5*(mu(i,j-1)+mu(i,j)))
  1810. rv_tend(i,k,j) = rv_tend(i,k,j)-cqv(i,k,j)*dpy
  1811. END DO
  1812. END DO
  1813. ELSE
  1814. ! ADT eqn 45: -mu*(my/mx)*(1/rho)*partial dp/dy
  1815. ! [alt, al are 1/rho terms; muv, mu are NOT coupled]
  1816. DO K=1,ktf
  1817. DO i = i_start, itf
  1818. ! Here are mu dp/dy terms 1-3; term 4 not needed if hydrostatic
  1819. dpy = (msfvy(i,j)/msfvx(i,j))*.5*rdy*muv(i,j)*( &
  1820. (ph (i,k+1,j)-ph (i,k+1,j-1) + ph(i,k,j)-ph(i,k,j-1)) &
  1821. +(alt(i,k ,j)+alt(i,k ,j-1))*(p (i,k,j)-p (i,k,j-1)) &
  1822. +(al (i,k ,j)+al (i,k ,j-1))*(pb(i,k,j)-pb(i,k,j-1)) )
  1823. rv_tend(i,k,j) = rv_tend(i,k,j)-cqv(i,k,j)*dpy
  1824. END DO
  1825. END DO
  1826. END IF
  1827. ENDDO
  1828. ! now the east-west (x) pressure gradient
  1829. itf=ite
  1830. jtf=MIN(jte,jde-1)
  1831. ktf=MIN(kte,kde-1)
  1832. i_start = its
  1833. j_start = jts
  1834. IF ( (config_flags%open_xs .or. specified .or. &
  1835. config_flags%nested ) .and. its == ids ) i_start = its+1
  1836. IF ( (config_flags%open_xe .or. specified .or. &
  1837. config_flags%nested ) .and. ite == ide ) itf = itf-1
  1838. IF ( config_flags%periodic_x ) i_start = its
  1839. IF ( config_flags%periodic_x ) itf=ite
  1840. DO j = j_start, jtf
  1841. IF ( non_hydrostatic ) THEN
  1842. k=1
  1843. DO i = i_start, itf
  1844. dpn(i,k) = .5*( cf1*(p(i-1,k ,j)+p(i,k ,j)) &
  1845. +cf2*(p(i-1,k+1,j)+p(i,k+1,j)) &
  1846. +cf3*(p(i-1,k+2,j)+p(i,k+2,j)) )
  1847. dpn(i,kde) = 0.
  1848. ENDDO
  1849. IF (top_lid) THEN
  1850. DO i = i_start, itf
  1851. dpn(i,kde) = .5*( cf1*(p(i-1,kde-1,j)+p(i,kde-1,j)) &
  1852. +cf2*(p(i-1,kde-2,j)+p(i,kde-2,j)) &
  1853. +cf3*(p(i-1,kde-3,j)+p(i,kde-3,j)) )
  1854. ENDDO
  1855. ENDIF
  1856. DO k=2,ktf
  1857. DO i = i_start, itf
  1858. dpn(i,k) = .5*( fnm(k)*(p(i-1,k ,j)+p(i,k ,j)) &
  1859. +fnp(k)*(p(i-1,k-1,j)+p(i,k-1,j)) )
  1860. END DO
  1861. END DO
  1862. ! ADT eqn 44: -mu*(mx/my)*(1/rho)*partial dp/dx
  1863. ! [alt, al are 1/rho terms; muu, mu are NOT coupled]
  1864. DO K=1,ktf
  1865. DO i = i_start, itf
  1866. ! Here are mu dp/dy terms 1-3
  1867. dpx = (msfux(i,j)/msfuy(i,j))*.5*rdx*muu(i,j)*( &
  1868. (ph (i,k+1,j)-ph (i-1,k+1,j) + ph(i,k,j)-ph(i-1,k,j)) &
  1869. +(alt(i,k ,j)+alt(i-1,k ,j))*(p (i,k,j)-p (i-1,k,j)) &
  1870. +(al (i,k ,j)+al (i-1,k ,j))*(pb(i,k,j)-pb(i-1,k,j)) )
  1871. ! Here is mu dp/dy term 4
  1872. dpx = dpx + (msfux(i,j)/msfuy(i,j))*rdx*(php(i,k,j)-php(i-1,k,j))* &
  1873. (rdnw(k)*(dpn(i,k+1)-dpn(i,k))-.5*(mu(i-1,j)+mu(i,j)))
  1874. ru_tend(i,k,j) = ru_tend(i,k,j)-cqu(i,k,j)*dpx
  1875. END DO
  1876. END DO
  1877. ELSE
  1878. ! ADT eqn 44: -mu*(mx/my)*(1/rho)*partial dp/dx
  1879. ! [alt, al are 1/rho terms; muu, mu are NOT coupled]
  1880. DO K=1,ktf
  1881. DO i = i_start, itf
  1882. ! Here are mu dp/dy terms 1-3; term 4 not needed if hydrostatic
  1883. dpx = (msfux(i,j)/msfuy(i,j))*.5*rdx*muu(i,j)*( &
  1884. (ph (i,k+1,j)-ph (i-1,k+1,j) + ph(i,k,j)-ph(i-1,k,j)) &
  1885. +(alt(i,k ,j)+alt(i-1,k ,j))*(p (i,k,j)-p (i-1,k,j)) &
  1886. +(al (i,k ,j)+al (i-1,k ,j))*(pb(i,k,j)-pb(i-1,k,j)) )
  1887. ru_tend(i,k,j) = ru_tend(i,k,j)-cqu(i,k,j)*dpx
  1888. END DO
  1889. END DO
  1890. END IF
  1891. ENDDO
  1892. END SUBROUTINE horizontal_pressure_gradient
  1893. !-------------------------------------------------------------------------------
  1894. SUBROUTINE pg_buoy_w( rw_tend, p, cqw, mu, mub, &
  1895. rdnw, rdn, g, msftx, msfty, &
  1896. ids, ide, jds, jde, kds, kde, &
  1897. ims, ime, jms, jme, kms, kme, &
  1898. its, ite, jts, jte, kts, kte )
  1899. IMPLICIT NONE
  1900. ! Input data
  1901. INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
  1902. ims, ime, jms, jme, kms, kme, &
  1903. its, ite, jts, jte, kts, kte
  1904. REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(IN ) :: p
  1905. REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(INOUT) :: cqw
  1906. REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(INOUT) :: rw_tend
  1907. REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN ) :: mub, mu, msftx, msfty
  1908. REAL, DIMENSION( kms:kme ), INTENT(IN ) :: rdnw, rdn
  1909. REAL, INTENT(IN ) :: g
  1910. INTEGER :: itf, jtf, i, j, k
  1911. REAL :: cq1, cq2
  1912. !<DESCRIPTION>
  1913. !
  1914. ! pg_buoy_w calculates the
  1915. ! vertical pressure gradient and buoyancy terms for the large-timestep
  1916. ! tendency in the vertical momentum equation.
  1917. !
  1918. !</DESCRIPTION>
  1919. ! BUOYANCY AND PRESSURE GRADIENT TERM IN W EQUATION AT TIME T
  1920. ! Map scale factor notes
  1921. ! ADT eqn 46 RHS terms 6 and 7 (where 7 is "-rho g")
  1922. ! Dividing by my, and using mu and nu (see Klemp et al. eqns 32, 40)
  1923. ! term 6: +(g/my) partial dp'/dnu
  1924. ! term 7: -(g/my) mu'
  1925. !
  1926. ! For moisture-free atmosphere, cq1=1, cq2=0
  1927. ! => (1./msft(i,j)) * g * [rdn(k)*{p(i,k,j)-p(i,k-1,j)}-mu(i,j)]
  1928. itf=MIN(ite,ide-1)
  1929. jtf=MIN(jte,jde-1)
  1930. DO j = jts,jtf
  1931. k=kde
  1932. DO i=its,itf
  1933. cq1 = 1./(1.+cqw(i,k-1,j))
  1934. cq2 = cqw(i,k-1,j)*cq1
  1935. rw_tend(i,k,j) = rw_tend(i,k,j)+(1./msfty(i,j))*g*( &
  1936. cq1*2.*rdnw(k-1)*( -p(i,k-1,j)) &
  1937. -mu(i,j)-cq2*mub(i,j) )
  1938. END DO
  1939. DO k = 2, kde-1
  1940. DO i = its,itf
  1941. cq1 = 1./(1.+cqw(i,k,j))
  1942. cq2 = cqw(i,k,j)*cq1
  1943. cqw(i,k,j) = cq1
  1944. rw_tend(i,k,j) = rw_tend(i,k,j)+(1./msfty(i,j))*g*( &
  1945. cq1*rdn(k)*(p(i,k,j)-p(i,k-1,j)) &
  1946. -mu(i,j)-cq2*mub(i,j) )
  1947. END DO
  1948. ENDDO
  1949. ENDDO
  1950. END SUBROUTINE pg_buoy_w
  1951. !-------------------------------------------------------------------------------
  1952. SUBROUTINE w_damp( rw_tend, max_vert_cfl,max_horiz_cfl, &
  1953. u, v, ww, w, mut, rdnw, &
  1954. rdx, rdy, msfux, msfuy, &
  1955. msfvx, msfvy, dt, &
  1956. config_flags, &
  1957. ids, ide, jds, jde, kds, kde, &
  1958. ims, ime, jms, jme, kms, kme, &
  1959. its, ite, jts, jte, kts, kte )
  1960. USE module_llxy
  1961. IMPLICIT NONE
  1962. ! Input data
  1963. TYPE(grid_config_rec_type ) , INTENT(IN ) :: config_flags
  1964. INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
  1965. ims, ime, jms, jme, kms, kme, &
  1966. its, ite, jts, jte, kts, kte
  1967. REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(IN ) :: u, v, ww, w
  1968. REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(INOUT) :: rw_tend
  1969. REAL, INTENT(OUT) :: max_vert_cfl
  1970. REAL, INTENT(OUT) :: max_horiz_cfl
  1971. REAL :: horiz_cfl
  1972. REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN ) :: mut
  1973. REAL, DIMENSION( kms:kme ), INTENT(IN ) :: rdnw
  1974. REAL, INTENT(IN) :: dt
  1975. REAL, INTENT(IN) :: rdx, rdy
  1976. REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: msfux, msfuy
  1977. REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: msfvx, msfvy
  1978. REAL :: vert_cfl, cf_n, cf_d, maxdub, maxdeta
  1979. INTEGER :: itf, jtf, i, j, k, maxi, maxj, maxk
  1980. INTEGER :: some
  1981. CHARACTER*512 :: temp
  1982. CHARACTER (LEN=256) :: time_str
  1983. CHARACTER (LEN=256) :: grid_str
  1984. integer :: total
  1985. REAL :: msfuxt , msfxffl
  1986. !<DESCRIPTION>
  1987. !
  1988. ! w_damp computes a damping term for the vertical velocity when the
  1989. ! vertical Courant number is too large. This was found to be preferable to
  1990. ! decreasing the timestep or increasing the diffusion in real-data applications
  1991. ! that produced potentially-unstable large vertical velocities because of
  1992. ! unphysically large heating rates coming from the cumulus parameterization
  1993. ! schemes run at moderately high resolutions (dx ~ O(10) km).
  1994. !
  1995. ! Additionally, w_damp returns the maximum cfl values due to vertical motion and
  1996. ! horizontal motion. These values are returned via the max_vert_cfl and
  1997. ! max_horiz_cfl variables. (Added by T. Hutchinson, WSI, 3/5/2007)
  1998. !
  1999. !</DESCRIPTION>
  2000. itf=MIN(ite,ide-1)
  2001. jtf=MIN(jte,jde-1)
  2002. some = 0
  2003. max_vert_cfl = 0.
  2004. max_horiz_cfl = 0.
  2005. total = 0
  2006. IF(config_flags%map_proj == PROJ_CASSINI ) then
  2007. msfxffl = 1.0/COS(config_flags%fft_filter_lat*degrad)
  2008. END IF
  2009. IF ( config_flags%w_damping == 1 ) THEN
  2010. DO j = jts,jtf
  2011. DO k = 2, kde-1
  2012. DO i = its,itf
  2013. #if 1
  2014. IF(config_flags%map_proj == PROJ_CASSINI ) then
  2015. msfuxt = MIN(msfux(i,j), msfxffl)
  2016. ELSE
  2017. msfuxt = msfux(i,j)
  2018. END IF
  2019. vert_cfl = abs(ww(i,k,j)/mut(i,j)*rdnw(k)*dt)
  2020. IF ( vert_cfl > max_vert_cfl ) THEN
  2021. max_vert_cfl = vert_cfl ; maxi = i ; maxj = j ; maxk = k
  2022. maxdub = w(i,k,j) ; maxdeta = -1./rdnw(k)
  2023. ENDIF
  2024. horiz_cfl = max( abs(u(i,k,j) * rdx * msfuxt * dt), &
  2025. abs(v(i,k,j) * rdy * msfvy(i,j) * dt) )
  2026. if (horiz_cfl > max_horiz_cfl) then
  2027. max_horiz_cfl = horiz_cfl
  2028. endif
  2029. if(vert_cfl .gt. w_beta)then
  2030. #else
  2031. ! restructure to get rid of divide
  2032. !
  2033. ! This had been used for efficiency, but with the addition of returning the cfl values,
  2034. ! the old version (above) was reinstated. (T. Hutchinson, 3/5/2007)
  2035. !
  2036. cf_n = abs(ww(i,k,j)*rdnw(k)*dt)
  2037. cf_d = abs(mut(i,j))
  2038. if(cf_n .gt. cf_d*w_beta )then
  2039. #endif
  2040. !$OMP CRITICAL(big_step_utilities)
  2041. WRITE(temp,*)i,j,k,' vert_cfl,w,d(eta)=',vert_cfl,w(i,k,j),-1./rdnw(k)
  2042. CALL wrf_debug ( 100 , TRIM(temp) )
  2043. !$OMP END CRITICAL(big_step_utilities)
  2044. if ( vert_cfl > 2. ) some = some + 1
  2045. rw_tend(i,k,j) = rw_tend(i,k,j)-sign(1.,w(i,k,j))*w_alpha*(vert_cfl-w_beta)*mut(i,j)
  2046. endif
  2047. END DO
  2048. ENDDO
  2049. ENDDO
  2050. ELSE
  2051. ! just print
  2052. DO j = jts,jtf
  2053. DO k = 2, kde-1
  2054. DO i = its,itf
  2055. #if 1
  2056. IF(config_flags%map_proj == PROJ_CASSINI ) then
  2057. msfuxt = MIN(msfux(i,j), msfxffl)
  2058. ELSE
  2059. msfuxt = msfux(i,j)
  2060. END IF
  2061. vert_cfl = abs(ww(i,k,j)/mut(i,j)*rdnw(k)*dt)
  2062. IF ( vert_cfl > max_vert_cfl ) THEN
  2063. max_vert_cfl = vert_cfl ; maxi = i ; maxj = j ; maxk = k
  2064. maxdub = w(i,k,j) ; maxdeta = -1./rdnw(k)
  2065. ENDIF
  2066. horiz_cfl = max( abs(u(i,k,j) * rdx * msfuxt * dt), &
  2067. abs(v(i,k,j) * rdy * msfvy(i,j) * dt) )
  2068. if (horiz_cfl > max_horiz_cfl) then
  2069. max_horiz_cfl = horiz_cfl
  2070. endif
  2071. if(vert_cfl .gt. w_beta)then
  2072. #else
  2073. ! restructure to get rid of divide
  2074. !
  2075. ! This had been used for efficiency, but with the addition of returning the cfl values,
  2076. ! the old version (above) was reinstated. (T. Hutchinson, 3/5/2007)
  2077. !
  2078. cf_n = abs(ww(i,k,j)*rdnw(k)*dt)
  2079. cf_d = abs(mut(i,j))
  2080. if(cf_n .gt. cf_d*w_beta )then
  2081. #endif
  2082. !$OMP CRITICAL(big_step_utilities)
  2083. WRITE(temp,*)i,j,k,' vert_cfl,w,d(eta)=',vert_cfl,w(i,k,j),-1./rdnw(k)
  2084. CALL wrf_debug ( 100 , TRIM(temp) )
  2085. !$OMP END CRITICAL(big_step_utilities)
  2086. if ( vert_cfl > 2. ) some = some + 1
  2087. endif
  2088. END DO
  2089. ENDDO
  2090. ENDDO
  2091. ENDIF
  2092. IF ( some .GT. 0 ) THEN
  2093. CALL get_current_time_string( time_str )
  2094. CALL get_current_grid_name( grid_str )
  2095. !$OMP CRITICAL(big_step_utilities)
  2096. WRITE(wrf_err_message,*)some, &
  2097. ' points exceeded cfl=2 in domain '//TRIM(grid_str)//' at time '//TRIM(time_str)//' hours'
  2098. CALL wrf_debug ( 0 , TRIM(wrf_err_message) )
  2099. WRITE(wrf_err_message,*)'MAX AT i,j,k: ',maxi,maxj,maxk,' vert_cfl,w,d(eta)=',max_vert_cfl, &
  2100. maxdub,maxdeta
  2101. CALL wrf_debug ( 0 , TRIM(wrf_err_message) )
  2102. !$OMP END CRITICAL(big_step_utilities)
  2103. ENDIF
  2104. END SUBROUTINE w_damp
  2105. !-------------------------------------------------------------------------------
  2106. SUBROUTINE horizontal_diffusion ( name, field, tendency, mu, &
  2107. config_flags, &
  2108. msfux, msfuy, msfvx, msfvx_inv, &
  2109. msfvy, msftx, msfty, &
  2110. khdif, xkmhd, rdx, rdy, &
  2111. ids, ide, jds, jde, kds, kde, &
  2112. ims, ime, jms, jme, kms, kme, &
  2113. its, ite, jts, jte, kts, kte )
  2114. IMPLICIT NONE
  2115. ! Input data
  2116. TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
  2117. INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
  2118. ims, ime, jms, jme, kms, kme, &
  2119. its, ite, jts, jte, kts, kte
  2120. CHARACTER(LEN=1) , INTENT(IN ) :: name
  2121. REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: field, xkmhd
  2122. REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency
  2123. REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mu
  2124. REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: msfux, &
  2125. msfuy, &
  2126. msfvx, &
  2127. msfvx_inv, &
  2128. msfvy, &
  2129. msftx, &
  2130. msfty
  2131. REAL , INTENT(IN ) :: rdx, &
  2132. rdy, &
  2133. khdif
  2134. ! Local data
  2135. INTEGER :: i, j, k, itf, jtf, ktf
  2136. INTEGER :: i_start, i_end, j_start, j_end
  2137. REAL :: mrdx, mkrdxm, mkrdxp, &
  2138. mrdy, mkrdym, mkrdyp
  2139. LOGICAL :: specified
  2140. !<DESCRIPTION>
  2141. !
  2142. ! horizontal_diffusion computes the horizontal diffusion tendency
  2143. ! on model horizontal coordinate surfaces.
  2144. !
  2145. !</DESCRIPTION>
  2146. specified = .false.
  2147. if(config_flags%specified .or. config_flags%nested) specified = .true.
  2148. ktf=MIN(kte,kde-1)
  2149. IF (name .EQ. 'u') THEN
  2150. i_start = its
  2151. i_end = ite
  2152. j_start = jts
  2153. j_end = MIN(jte,jde-1)
  2154. IF ( config_flags%open_xs .or. specified ) i_start = MAX(ids+1,its)
  2155. IF ( config_flags%open_xe .or. specified ) i_end = MIN(ide-1,ite)
  2156. IF ( config_flags%open_ys .or. specified ) j_start = MAX(jds+1,jts)
  2157. IF ( config_flags%open_ye .or. specified ) j_end = MIN(jde-2,jte)
  2158. IF ( config_flags%periodic_x ) i_start = its
  2159. IF ( config_flags%periodic_x ) i_end = ite
  2160. DO j = j_start, j_end
  2161. DO k=kts,ktf
  2162. DO i = i_start, i_end
  2163. ! The interior is grad: (m_x*d/dx), the exterior is div: (m_x*m_y*d/dx(/m_y))
  2164. ! setting up different averagings of m^2 partial d/dX and m^2 partial d/dY
  2165. mkrdxm=(msftx(i-1,j)/msfty(i-1,j))*mu(i-1,j)*xkmhd(i-1,k,j)*rdx
  2166. mkrdxp=(msftx(i,j)/msfty(i,j))*mu(i,j)*xkmhd(i,k,j)*rdx
  2167. mrdx=msfux(i,j)*msfuy(i,j)*rdx
  2168. mkrdym=( (msfuy(i,j)+msfuy(i,j-1))/(msfux(i,j)+msfux(i,j-1)) )* &
  2169. 0.25*(mu(i,j)+mu(i,j-1)+mu(i-1,j-1)+mu(i-1,j))* &
  2170. 0.25*(xkmhd(i,k,j)+xkmhd(i,k,j-1)+xkmhd(i-1,k,j-1)+xkmhd(i-1,k,j))*rdy
  2171. mkrdyp=( (msfuy(i,j)+msfuy(i,j+1))/(msfux(i,j)+msfux(i,j+1)) )* &
  2172. 0.25*(mu(i,j)+mu(i,j+1)+mu(i-1,j+1)+mu(i-1,j))* &
  2173. 0.25*(xkmhd(i,k,j)+xkmhd(i,k,j+1)+xkmhd(i-1,k,j+1)+xkmhd(i-1,k,j))*rdy
  2174. ! need to do four-corners (t) for diffusion coefficient as there are
  2175. ! no values at u,v points
  2176. ! msfuy - has to be y as part of d/dY
  2177. ! has to be u as we're at a u point
  2178. mrdy=msfux(i,j)*msfuy(i,j)*rdy
  2179. ! correctly averaged version of rho~ * m^2 *
  2180. ! [partial d/dX(partial du^/dX) + partial d/dY(partial du^/dY)]
  2181. tendency(i,k,j)=tendency(i,k,j)+( &
  2182. mrdx*(mkrdxp*(field(i+1,k,j)-field(i ,k,j)) &
  2183. -mkrdxm*(field(i ,k,j)-field(i-1,k,j))) &
  2184. +mrdy*(mkrdyp*(field(i,k,j+1)-field(i,k,j )) &
  2185. -mkrdym*(field(i,k,j )-field(i,k,j-1))))
  2186. ENDDO
  2187. ENDDO
  2188. ENDDO
  2189. ELSE IF (name .EQ. 'v')THEN
  2190. i_start = its
  2191. i_end = MIN(ite,ide-1)
  2192. j_start = jts
  2193. j_end = jte
  2194. IF ( config_flags%open_xs .or. specified ) i_start = MAX(ids+1,its)
  2195. IF ( config_flags%open_xe .or. specified ) i_end = MIN(ide-2,ite)
  2196. IF ( config_flags%open_ys .or. specified ) j_start = MAX(jds+1,jts)
  2197. IF ( config_flags%open_ye .or. specified ) j_end = MIN(jde-1,jte)
  2198. IF ( config_flags%periodic_x ) i_start = its
  2199. IF ( config_flags%periodic_x ) i_end = MIN(ite,ide-1)
  2200. IF ( config_flags%polar ) j_start = MAX(jds+1,jts)
  2201. IF ( config_flags%polar ) j_end = MIN(jde-1,jte)
  2202. DO j = j_start, j_end
  2203. DO k=kts,ktf
  2204. DO i = i_start, i_end
  2205. mkrdxm=( (msfvx(i,j)+msfvx(i-1,j))/(msfvy(i,j)+msfvy(i-1,j)) )* &
  2206. 0.25*(mu(i,j)+mu(i,j-1)+mu(i-1,j-1)+mu(i-1,j))* &
  2207. 0.25*(xkmhd(i,k,j)+xkmhd(i,k,j-1)+xkmhd(i-1,k,j-1)+xkmhd(i-1,k,j))*rdx
  2208. mkrdxp=( (msfvx(i,j)+msfvx(i+1,j))/(msfvy(i,j)+msfvy(i+1,j)) )* &
  2209. 0.25*(mu(i,j)+mu(i,j-1)+mu(i+1,j-1)+mu(i+1,j))* &
  2210. 0.25*(xkmhd(i,k,j)+xkmhd(i,k,j-1)+xkmhd(i+1,k,j-1)+xkmhd(i+1,k,j))*rdx
  2211. mrdx=msfvx(i,j)*msfvy(i,j)*rdx
  2212. mkrdym=(msfty(i,j-1)/msftx(i,j-1))*xkmhd(i,k,j-1)*rdy
  2213. mkrdyp=(msfty(i,j)/msftx(i,j))*xkmhd(i,k,j)*rdy
  2214. mrdy=msfvx(i,j)*msfvy(i,j)*rdy
  2215. tendency(i,k,j)=tendency(i,k,j)+( &
  2216. mrdx*(mkrdxp*(field(i+1,k,j)-field(i ,k,j)) &
  2217. -mkrdxm*(field(i ,k,j)-field(i-1,k,j))) &
  2218. +mrdy*(mkrdyp*(field(i,k,j+1)-field(i,k,j )) &
  2219. -mkrdym*(field(i,k,j )-field(i,k,j-1))))
  2220. ENDDO
  2221. ENDDO
  2222. ENDDO
  2223. ELSE IF (name .EQ. 'w')THEN
  2224. i_start = its
  2225. i_end = MIN(ite,ide-1)
  2226. j_start = jts
  2227. j_end = MIN(jte,jde-1)
  2228. IF ( config_flags%open_xs .or. specified ) i_start = MAX(ids+1,its)
  2229. IF ( config_flags%open_xe .or. specified ) i_end = MIN(ide-2,ite)
  2230. IF ( config_flags%open_ys .or. specified ) j_start = MAX(jds+1,jts)
  2231. IF ( config_flags%open_ye .or. specified ) j_end = MIN(jde-2,jte)
  2232. IF ( config_flags%periodic_x ) i_start = its
  2233. IF ( config_flags%periodic_x ) i_end = MIN(ite,ide-1)
  2234. DO j = j_start, j_end
  2235. DO k=kts+1,ktf
  2236. DO i = i_start, i_end
  2237. mkrdxm=(msfux(i,j)/msfuy(i,j))* &
  2238. 0.25*(mu(i,j)+mu(i-1,j)+mu(i,j)+mu(i-1,j))* &
  2239. 0.25*(xkmhd(i,k,j)+xkmhd(i-1,k,j)+xkmhd(i,k-1,j)+xkmhd(i-1,k-1,j))*rdx
  2240. mkrdxp=(msfux(i+1,j)/msfuy(i+1,j))* &
  2241. 0.25*(mu(i+1,j)+mu(i,j)+mu(i+1,j)+mu(i,j))* &
  2242. 0.25*(xkmhd(i+1,k,j)+xkmhd(i,k,j)+xkmhd(i+1,k-1,j)+xkmhd(i,k-1,j))*rdx
  2243. mrdx=msftx(i,j)*msfty(i,j)*rdx
  2244. ! mkrdym=(msfvy(i,j)/msfvx(i,j))* &
  2245. mkrdym=(msfvy(i,j)*msfvx_inv(i,j))* &
  2246. 0.25*(mu(i,j)+mu(i,j-1)+mu(i,j)+mu(i,j-1))* &
  2247. 0.25*(xkmhd(i,k,j)+xkmhd(i,k,j-1)+xkmhd(i,k-1,j)+xkmhd(i,k-1,j-1))*rdy
  2248. ! mkrdyp=(msfvy(i,j+1)/msfvx(i,j+1))* &
  2249. mkrdyp=(msfvy(i,j+1)*msfvx_inv(i,j+1))* &
  2250. 0.25*(mu(i,j+1)+mu(i,j)+mu(i,j+1)+mu(i,j))* &
  2251. 0.25*(xkmhd(i,k,j+1)+xkmhd(i,k,j)+xkmhd(i,k-1,j+1)+xkmhd(i,k-1,j))*rdy
  2252. mrdy=msftx(i,j)*msfty(i,j)*rdy
  2253. tendency(i,k,j)=tendency(i,k,j)+( &
  2254. mrdx*(mkrdxp*(field(i+1,k,j)-field(i ,k,j)) &
  2255. -mkrdxm*(field(i ,k,j)-field(i-1,k,j))) &
  2256. +mrdy*(mkrdyp*(field(i,k,j+1)-field(i,k,j )) &
  2257. -mkrdym*(field(i,k,j )-field(i,k,j-1))))
  2258. ENDDO
  2259. ENDDO
  2260. ENDDO
  2261. ELSE
  2262. i_start = its
  2263. i_end = MIN(ite,ide-1)
  2264. j_start = jts
  2265. j_end = MIN(jte,jde-1)
  2266. IF ( config_flags%open_xs .or. specified ) i_start = MAX(ids+1,its)
  2267. IF ( config_flags%open_xe .or. specified ) i_end = MIN(ide-2,ite)
  2268. IF ( config_flags%open_ys .or. specified ) j_start = MAX(jds+1,jts)
  2269. IF ( config_flags%open_ye .or. specified ) j_end = MIN(jde-2,jte)
  2270. IF ( config_flags%periodic_x ) i_start = its
  2271. IF ( config_flags%periodic_x ) i_end = MIN(ite,ide-1)
  2272. DO j = j_start, j_end
  2273. DO k=kts,ktf
  2274. DO i = i_start, i_end
  2275. mkrdxm=(msfux(i,j)/msfuy(i,j))*0.5*(xkmhd(i,k,j)+xkmhd(i-1,k,j))*0.5*(mu(i,j)+mu(i-1,j))*rdx
  2276. mkrdxp=(msfux(i+1,j)/msfuy(i+1,j))*0.5*(xkmhd(i+1,k,j)+xkmhd(i,k,j))*0.5*(mu(i+1,j)+mu(i,j))*rdx
  2277. mrdx=msftx(i,j)*msfty(i,j)*rdx
  2278. ! mkrdym=(msfvy(i,j)/msfvx(i,j))*0.5*(xkmhd(i,k,j)+xkmhd(i,k,j-1))*0.5*(mu(i,j)+mu(i,j-1))*rdy
  2279. mkrdym=(msfvy(i,j)*msfvx_inv(i,j))*0.5*(xkmhd(i,k,j)+xkmhd(i,k,j-1))*0.5*(mu(i,j)+mu(i,j-1))*rdy
  2280. ! mkrdyp=(msfvy(i,j+1)/msfvx(i,j+1))*0.5*(xkmhd(i,k,j+1)+xkmhd(i,k,j))*0.5*(mu(i,j+1)+mu(i,j))*rdy
  2281. mkrdyp=(msfvy(i,j+1)*msfvx_inv(i,j+1))*0.5*(xkmhd(i,k,j+1)+xkmhd(i,k,j))*0.5*(mu(i,j+1)+mu(i,j))*rdy
  2282. mrdy=msftx(i,j)*msfty(i,j)*rdy
  2283. tendency(i,k,j)=tendency(i,k,j)+( &
  2284. mrdx*(mkrdxp*(field(i+1,k,j)-field(i ,k,j)) &
  2285. -mkrdxm*(field(i ,k,j)-field(i-1,k,j))) &
  2286. +mrdy*(mkrdyp*(field(i,k,j+1)-field(i,k,j )) &
  2287. -mkrdym*(field(i,k,j )-field(i,k,j-1))))
  2288. ENDDO
  2289. ENDDO
  2290. ENDDO
  2291. ENDIF
  2292. END SUBROUTINE horizontal_diffusion
  2293. !-----------------------------------------------------------------------------------------
  2294. SUBROUTINE horizontal_diffusion_3dmp ( name, field, tendency, mu, &
  2295. config_flags, base_3d, &
  2296. msfux, msfuy, msfvx, msfvx_inv, &
  2297. msfvy, msftx, msfty, &
  2298. khdif, xkmhd, rdx, rdy, &
  2299. ids, ide, jds, jde, kds, kde, &
  2300. ims, ime, jms, jme, kms, kme, &
  2301. its, ite, jts, jte, kts, kte )
  2302. IMPLICIT NONE
  2303. ! Input data
  2304. TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
  2305. INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
  2306. ims, ime, jms, jme, kms, kme, &
  2307. its, ite, jts, jte, kts, kte
  2308. CHARACTER(LEN=1) , INTENT(IN ) :: name
  2309. REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: field, &
  2310. xkmhd, &
  2311. base_3d
  2312. REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency
  2313. REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mu
  2314. REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: msfux, &
  2315. msfuy, &
  2316. msfvx, &
  2317. msfvx_inv, &
  2318. msfvy, &
  2319. msftx, &
  2320. msfty
  2321. REAL , INTENT(IN ) :: rdx, &
  2322. rdy, &
  2323. khdif
  2324. ! Local data
  2325. INTEGER :: i, j, k, itf, jtf, ktf
  2326. INTEGER :: i_start, i_end, j_start, j_end
  2327. REAL :: mrdx, mkrdxm, mkrdxp, &
  2328. mrdy, mkrdym, mkrdyp
  2329. LOGICAL :: specified
  2330. !<DESCRIPTION>
  2331. !
  2332. ! horizontal_diffusion_3dmp computes the horizontal diffusion tendency
  2333. ! on model horizontal coordinate surfaces. This routine computes diffusion
  2334. ! a perturbation scalar (field-base_3d).
  2335. !
  2336. !</DESCRIPTION>
  2337. specified = .false.
  2338. if(config_flags%specified .or. config_flags%nested) specified = .true.
  2339. ktf=MIN(kte,kde-1)
  2340. i_start = its
  2341. i_end = MIN(ite,ide-1)
  2342. j_start = jts
  2343. j_end = MIN(jte,jde-1)
  2344. IF ( config_flags%open_xs .or. specified ) i_start = MAX(ids+1,its)
  2345. IF ( config_flags%open_xe .or. specified ) i_end = MIN(ide-2,ite)
  2346. IF ( config_flags%open_ys .or. specified ) j_start = MAX(jds+1,jts)
  2347. IF ( config_flags%open_ye .or. specified ) j_end = MIN(jde-2,jte)
  2348. IF ( config_flags%periodic_x ) i_start = its
  2349. IF ( config_flags%periodic_x ) i_end = MIN(ite,ide-1)
  2350. DO j = j_start, j_end
  2351. DO k=kts,ktf
  2352. DO i = i_start, i_end
  2353. mkrdxm=(msfux(i,j)/msfuy(i,j))*0.5*(xkmhd(i,k,j)+xkmhd(i-1,k,j))*0.5*(mu(i,j)+mu(i-1,j))*rdx
  2354. mkrdxp=(msfux(i+1,j)/msfuy(i+1,j))*0.5*(xkmhd(i+1,k,j)+xkmhd(i,k,j))*0.5*(mu(i+1,j)+mu(i,j))*rdx
  2355. mrdx=msftx(i,j)*msfty(i,j)*rdx
  2356. ! mkrdym=(msfvy(i,j)/msfvx(i,j))*0.5*(xkmhd(i,k,j)+xkmhd(i,k,j-1))*0.5*(mu(i,j)+mu(i,j-1))*rdy
  2357. ! mkrdyp=(msfvy(i,j+1)/msfvx(i,j+1))*0.5*(xkmhd(i,k,j+1)+xkmhd(i,k,j))*0.5*(mu(i,j+1)+mu(i,j))*rdy
  2358. mkrdym=(msfvy(i,j)*msfvx_inv(i,j))*0.5*(xkmhd(i,k,j)+xkmhd(i,k,j-1))*0.5*(mu(i,j)+mu(i,j-1))*rdy
  2359. mkrdyp=(msfvy(i,j+1)*msfvx_inv(i,j+1))*0.5*(xkmhd(i,k,j+1)+xkmhd(i,k,j))*0.5*(mu(i,j+1)+mu(i,j))*rdy
  2360. mrdy=msftx(i,j)*msfty(i,j)*rdy
  2361. tendency(i,k,j)=tendency(i,k,j)+( &
  2362. mrdx*( mkrdxp*( field(i+1,k,j) -field(i ,k,j) &
  2363. -base_3d(i+1,k,j)+base_3d(i ,k,j) ) &
  2364. -mkrdxm*( field(i ,k,j) -field(i-1,k,j) &
  2365. -base_3d(i ,k,j)+base_3d(i-1,k,j) ) ) &
  2366. +mrdy*( mkrdyp*( field(i,k,j+1) -field(i,k,j ) &
  2367. -base_3d(i,k,j+1)+base_3d(i,k,j ) ) &
  2368. -mkrdym*( field(i,k,j ) -field(i,k,j-1) &
  2369. -base_3d(i,k,j )+base_3d(i,k,j-1) ) ) &
  2370. )
  2371. ENDDO
  2372. ENDDO
  2373. ENDDO
  2374. END SUBROUTINE horizontal_diffusion_3dmp
  2375. !-----------------------------------------------------------------------------------------
  2376. SUBROUTINE vertical_diffusion ( name, field, tendency, &
  2377. config_flags, &
  2378. alt, mut, rdn, rdnw, kvdif, &
  2379. ids, ide, jds, jde, kds, kde, &
  2380. ims, ime, jms, jme, kms, kme, &
  2381. its, ite, jts, jte, kts, kte )
  2382. IMPLICIT NONE
  2383. ! Input data
  2384. TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
  2385. INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
  2386. ims, ime, jms, jme, kms, kme, &
  2387. its, ite, jts, jte, kts, kte
  2388. CHARACTER(LEN=1) , INTENT(IN ) :: name
  2389. REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , &
  2390. INTENT(IN ) :: field, &
  2391. alt
  2392. REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency
  2393. REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mut
  2394. REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: rdn, rdnw
  2395. REAL , INTENT(IN ) :: kvdif
  2396. ! Local data
  2397. INTEGER :: i, j, k, itf, jtf, ktf
  2398. INTEGER :: i_start, i_end, j_start, j_end
  2399. REAL , DIMENSION(its:ite, jts:jte) :: vfluxm, vfluxp, zz
  2400. REAL , DIMENSION(its:ite, 0:kte+1) :: vflux
  2401. REAL :: rdz
  2402. LOGICAL :: specified
  2403. !<DESCRIPTION>
  2404. !
  2405. ! vertical_diffusion
  2406. ! computes vertical diffusion tendency.
  2407. !
  2408. !</DESCRIPTION>
  2409. specified = .false.
  2410. if(config_flags%specified .or. config_flags%nested) specified = .true.
  2411. ktf=MIN(kte,kde-1)
  2412. IF (name .EQ. 'w')THEN
  2413. i_start = its
  2414. i_end = MIN(ite,ide-1)
  2415. j_start = jts
  2416. j_end = MIN(jte,jde-1)
  2417. j_loop_w : DO j = j_start, j_end
  2418. DO k=kts,ktf-1
  2419. DO i = i_start, i_end
  2420. vflux(i,k)= (kvdif/alt(i,k,j))*rdnw(k)*(field(i,k+1,j)-field(i,k,j))
  2421. ENDDO
  2422. ENDDO
  2423. DO i = i_start, i_end
  2424. vflux(i,ktf)=0.
  2425. ENDDO
  2426. DO k=kts+1,ktf
  2427. DO i = i_start, i_end
  2428. tendency(i,k,j)=tendency(i,k,j) &
  2429. +rdn(k)*g*g/mut(i,j)/(0.5*(alt(i,k,j)+alt(i,k-1,j))) &
  2430. *(vflux(i,k)-vflux(i,k-1))
  2431. ENDDO
  2432. ENDDO
  2433. ENDDO j_loop_w
  2434. ELSE IF(name .EQ. 'm')THEN
  2435. i_start = its
  2436. i_end = MIN(ite,ide-1)
  2437. j_start = jts
  2438. j_end = MIN(jte,jde-1)
  2439. j_loop_s : DO j = j_start, j_end
  2440. DO k=kts,ktf-1
  2441. DO i = i_start, i_end
  2442. vflux(i,k)=kvdif*rdn(k+1)/(0.5*(alt(i,k,j)+alt(i,k+1,j))) &
  2443. *(field(i,k+1,j)-field(i,k,j))
  2444. ENDDO
  2445. ENDDO
  2446. DO i = i_start, i_end
  2447. vflux(i,0)=vflux(i,1)
  2448. ENDDO
  2449. DO i = i_start, i_end
  2450. vflux(i,ktf)=0.
  2451. ENDDO
  2452. DO k=kts,ktf
  2453. DO i = i_start, i_end
  2454. tendency(i,k,j)=tendency(i,k,j)+g*g/mut(i,j)/alt(i,k,j) &
  2455. *rdnw(k)*(vflux(i,k)-vflux(i,k-1))
  2456. ENDDO
  2457. ENDDO
  2458. ENDDO j_loop_s
  2459. ENDIF
  2460. END SUBROUTINE vertical_diffusion
  2461. !-------------------------------------------------------------------------------
  2462. SUBROUTINE vertical_diffusion_mp ( field, tendency, config_flags, &
  2463. base, &
  2464. alt, mut, rdn, rdnw, kvdif, &
  2465. ids, ide, jds, jde, kds, kde, &
  2466. ims, ime, jms, jme, kms, kme, &
  2467. its, ite, jts, jte, kts, kte )
  2468. IMPLICIT NONE
  2469. ! Input data
  2470. TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
  2471. INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
  2472. ims, ime, jms, jme, kms, kme, &
  2473. its, ite, jts, jte, kts, kte
  2474. REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , &
  2475. INTENT(IN ) :: field, &
  2476. alt
  2477. REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency
  2478. REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mut
  2479. REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: rdn, &
  2480. rdnw, &
  2481. base
  2482. REAL , INTENT(IN ) :: kvdif
  2483. ! Local data
  2484. INTEGER :: i, j, k, itf, jtf, ktf
  2485. INTEGER :: i_start, i_end, j_start, j_end
  2486. REAL , DIMENSION(its:ite, 0:kte+1) :: vflux
  2487. REAL :: rdz
  2488. LOGICAL :: specified
  2489. !<DESCRIPTION>
  2490. !
  2491. ! vertical_diffusion_mp
  2492. ! computes vertical diffusion tendency of a perturbation variable
  2493. ! (field-base). Note that base as a 1D (k) field.
  2494. !
  2495. !</DESCRIPTION>
  2496. specified = .false.
  2497. if(config_flags%specified .or. config_flags%nested) specified = .true.
  2498. ktf=MIN(kte,kde-1)
  2499. i_start = its
  2500. i_end = MIN(ite,ide-1)
  2501. j_start = jts
  2502. j_end = MIN(jte,jde-1)
  2503. j_loop_s : DO j = j_start, j_end
  2504. DO k=kts,ktf-1
  2505. DO i = i_start, i_end
  2506. vflux(i,k)=kvdif*rdn(k+1)/(0.5*(alt(i,k,j)+alt(i,k+1,j))) &
  2507. *(field(i,k+1,j)-field(i,k,j)-base(k+1)+base(k))
  2508. ENDDO
  2509. ENDDO
  2510. DO i = i_start, i_end
  2511. vflux(i,0)=vflux(i,1)
  2512. ENDDO
  2513. DO i = i_start, i_end
  2514. vflux(i,ktf)=0.
  2515. ENDDO
  2516. DO k=kts,ktf
  2517. DO i = i_start, i_end
  2518. tendency(i,k,j)=tendency(i,k,j)+g*g/mut(i,j)/alt(i,k,j) &
  2519. *rdnw(k)*(vflux(i,k)-vflux(i,k-1))
  2520. ENDDO
  2521. ENDDO
  2522. ENDDO j_loop_s
  2523. END SUBROUTINE vertical_diffusion_mp
  2524. !-------------------------------------------------------------------------------
  2525. SUBROUTINE vertical_diffusion_3dmp ( field, tendency, config_flags, &
  2526. base_3d, &
  2527. alt, mut, rdn, rdnw, kvdif, &
  2528. ids, ide, jds, jde, kds, kde, &
  2529. ims, ime, jms, jme, kms, kme, &
  2530. its, ite, jts, jte, kts, kte )
  2531. IMPLICIT NONE
  2532. ! Input data
  2533. TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
  2534. INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
  2535. ims, ime, jms, jme, kms, kme, &
  2536. its, ite, jts, jte, kts, kte
  2537. REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , &
  2538. INTENT(IN ) :: field, &
  2539. alt, &
  2540. base_3d
  2541. REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency
  2542. REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mut
  2543. REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: rdn, &
  2544. rdnw
  2545. REAL , INTENT(IN ) :: kvdif
  2546. ! Local data
  2547. INTEGER :: i, j, k, itf, jtf, ktf
  2548. INTEGER :: i_start, i_end, j_start, j_end
  2549. REAL , DIMENSION(its:ite, 0:kte+1) :: vflux
  2550. REAL :: rdz
  2551. LOGICAL :: specified
  2552. !<DESCRIPTION>
  2553. !
  2554. ! vertical_diffusion_3dmp
  2555. ! computes vertical diffusion tendency of a perturbation variable
  2556. ! (field-base_3d).
  2557. !
  2558. !</DESCRIPTION>
  2559. specified = .false.
  2560. if(config_flags%specified .or. config_flags%nested) specified = .true.
  2561. ktf=MIN(kte,kde-1)
  2562. i_start = its
  2563. i_end = MIN(ite,ide-1)
  2564. j_start = jts
  2565. j_end = MIN(jte,jde-1)
  2566. j_loop_s : DO j = j_start, j_end
  2567. DO k=kts,ktf-1
  2568. DO i = i_start, i_end
  2569. vflux(i,k)=kvdif*rdn(k+1)/(0.5*(alt(i,k,j)+alt(i,k+1,j))) &
  2570. *( field(i,k+1,j) -field(i,k,j) &
  2571. -base_3d(i,k+1,j)+base_3d(i,k,j) )
  2572. ENDDO
  2573. ENDDO
  2574. DO i = i_start, i_end
  2575. vflux(i,0)=vflux(i,1)
  2576. ENDDO
  2577. DO i = i_start, i_end
  2578. vflux(i,ktf)=0.
  2579. ENDDO
  2580. DO k=kts,ktf
  2581. DO i = i_start, i_end
  2582. tendency(i,k,j)=tendency(i,k,j)+g*g/mut(i,j)/alt(i,k,j) &
  2583. *rdnw(k)*(vflux(i,k)-vflux(i,k-1))
  2584. ENDDO
  2585. ENDDO
  2586. ENDDO j_loop_s
  2587. END SUBROUTINE vertical_diffusion_3dmp
  2588. !-------------------------------------------------------------------------------
  2589. SUBROUTINE vertical_diffusion_u ( field, tendency, &
  2590. config_flags, u_base, &
  2591. alt, muu, rdn, rdnw, kvdif, &
  2592. ids, ide, jds, jde, kds, kde, &
  2593. ims, ime, jms, jme, kms, kme, &
  2594. its, ite, jts, jte, kts, kte )
  2595. IMPLICIT NONE
  2596. ! Input data
  2597. TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
  2598. INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
  2599. ims, ime, jms, jme, kms, kme, &
  2600. its, ite, jts, jte, kts, kte
  2601. REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , &
  2602. INTENT(IN ) :: field, &
  2603. alt
  2604. REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency
  2605. REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: muu
  2606. REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: rdn, rdnw, u_base
  2607. REAL , INTENT(IN ) :: kvdif
  2608. ! Local data
  2609. INTEGER :: i, j, k, itf, jtf, ktf
  2610. INTEGER :: i_start, i_end, j_start, j_end
  2611. REAL , DIMENSION(its:ite, 0:kte+1) :: vflux
  2612. REAL :: rdz, zz
  2613. LOGICAL :: specified
  2614. !<DESCRIPTION>
  2615. !
  2616. ! vertical_diffusion_u computes vertical diffusion tendency for
  2617. ! the u momentum equation. This routine assumes a constant eddy
  2618. ! viscosity kvdif.
  2619. !
  2620. !</DESCRIPTION>
  2621. specified = .false.
  2622. if(config_flags%specified .or. config_flags%nested) specified = .true.
  2623. ktf=MIN(kte,kde-1)
  2624. i_start = its
  2625. i_end = ite
  2626. j_start = jts
  2627. j_end = MIN(jte,jde-1)
  2628. IF ( config_flags%open_xs .or. specified ) i_start = MAX(ids+1,its)
  2629. IF ( config_flags%open_xe .or. specified ) i_end = MIN(ide-1,ite)
  2630. IF ( config_flags%periodic_x ) i_start = its
  2631. IF ( config_flags%periodic_x ) i_end = ite
  2632. j_loop_u : DO j = j_start, j_end
  2633. DO k=kts,ktf-1
  2634. DO i = i_start, i_end
  2635. vflux(i,k)=kvdif*rdn(k+1)/(0.25*( alt(i ,k ,j) &
  2636. +alt(i-1,k ,j) &
  2637. +alt(i ,k+1,j) &
  2638. +alt(i-1,k+1,j) ) ) &
  2639. *(field(i,k+1,j)-field(i,k,j) &
  2640. -u_base(k+1) +u_base(k) )
  2641. ENDDO
  2642. ENDDO
  2643. DO i = i_start, i_end
  2644. vflux(i,0)=vflux(i,1)
  2645. ENDDO
  2646. DO i = i_start, i_end
  2647. vflux(i,ktf)=0.
  2648. ENDDO
  2649. DO k=kts,ktf-1
  2650. DO i = i_start, i_end
  2651. tendency(i,k,j)=tendency(i,k,j)+ &
  2652. g*g*rdnw(k)/muu(i,j)/(0.5*(alt(i-1,k,j)+alt(i,k,j)))* &
  2653. (vflux(i,k)-vflux(i,k-1))
  2654. ENDDO
  2655. ENDDO
  2656. ENDDO j_loop_u
  2657. END SUBROUTINE vertical_diffusion_u
  2658. !-------------------------------------------------------------------------------
  2659. SUBROUTINE vertical_diffusion_v ( field, tendency, &
  2660. config_flags, v_base, &
  2661. alt, muv, rdn, rdnw, kvdif, &
  2662. ids, ide, jds, jde, kds, kde, &
  2663. ims, ime, jms, jme, kms, kme, &
  2664. its, ite, jts, jte, kts, kte )
  2665. IMPLICIT NONE
  2666. ! Input data
  2667. TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
  2668. INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
  2669. ims, ime, jms, jme, kms, kme, &
  2670. its, ite, jts, jte, kts, kte
  2671. REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , &
  2672. INTENT(IN ) :: field, &
  2673. alt
  2674. REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: rdn, rdnw, v_base
  2675. REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency
  2676. REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: muv
  2677. REAL , INTENT(IN ) :: kvdif
  2678. ! Local data
  2679. INTEGER :: i, j, k, itf, jtf, ktf, jm1
  2680. INTEGER :: i_start, i_end, j_start, j_end
  2681. REAL , DIMENSION(its:ite, 0:kte+1) :: vflux
  2682. REAL :: rdz, zz
  2683. LOGICAL :: specified
  2684. !<DESCRIPTION>
  2685. !
  2686. ! vertical_diffusion_v computes vertical diffusion tendency for
  2687. ! the v momentum equation. This routine assumes a constant eddy
  2688. ! viscosity kvdif.
  2689. !
  2690. !</DESCRIPTION>
  2691. specified = .false.
  2692. if(config_flags%specified .or. config_flags%nested) specified = .true.
  2693. ktf=MIN(kte,kde-1)
  2694. i_start = its
  2695. i_end = MIN(ite,ide-1)
  2696. j_start = jts
  2697. j_end = MIN(jte,jde-1)
  2698. IF ( config_flags%open_ys .or. specified ) j_start = MAX(jds+1,jts)
  2699. IF ( config_flags%open_ye .or. specified ) j_end = MIN(jde-1,jte)
  2700. j_loop_v : DO j = j_start, j_end
  2701. ! jm1 = max(j-1,1)
  2702. jm1 = j-1
  2703. DO k=kts,ktf-1
  2704. DO i = i_start, i_end
  2705. vflux(i,k)=kvdif*rdn(k+1)/(0.25*( alt(i,k ,j ) &
  2706. +alt(i,k ,jm1) &
  2707. +alt(i,k+1,j ) &
  2708. +alt(i,k+1,jm1) ) ) &
  2709. *(field(i,k+1,j)-field(i,k,j) &
  2710. -v_base(k+1) +v_base(k) )
  2711. ENDDO
  2712. ENDDO
  2713. DO i = i_start, i_end
  2714. vflux(i,0)=vflux(i,1)
  2715. ENDDO
  2716. DO i = i_start, i_end
  2717. vflux(i,ktf)=0.
  2718. ENDDO
  2719. DO k=kts,ktf-1
  2720. DO i = i_start, i_end
  2721. tendency(i,k,j)=tendency(i,k,j)+ &
  2722. g*g*rdnw(k)/muv(i,j)/(0.5*(alt(i,k,jm1)+alt(i,k,j)))* &
  2723. (vflux(i,k)-vflux(i,k-1))
  2724. ENDDO
  2725. ENDDO
  2726. ENDDO j_loop_v
  2727. END SUBROUTINE vertical_diffusion_v
  2728. !*************** end new mass coordinate routines
  2729. !-------------------------------------------------------------------------------
  2730. SUBROUTINE calculate_full ( rfield, rfieldb, rfieldp, &
  2731. ids, ide, jds, jde, kds, kde, &
  2732. ims, ime, jms, jme, kms, kme, &
  2733. its, ite, jts, jte, kts, kte )
  2734. IMPLICIT NONE
  2735. ! Input data
  2736. INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
  2737. ims, ime, jms, jme, kms, kme, &
  2738. its, ite, jts, jte, kts, kte
  2739. REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: rfieldb, &
  2740. rfieldp
  2741. REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT ) :: rfield
  2742. ! Local indices.
  2743. INTEGER :: i, j, k, itf, jtf, ktf
  2744. !<DESCRIPTION>
  2745. !
  2746. ! calculate_full
  2747. ! calculates full 3D field from pertubation and base field.
  2748. !
  2749. !</DESCRIPTION>
  2750. itf=MIN(ite,ide-1)
  2751. jtf=MIN(jte,jde-1)
  2752. ktf=MIN(kte,kde-1)
  2753. DO j=jts,jtf
  2754. DO k=kts,ktf
  2755. DO i=its,itf
  2756. rfield(i,k,j)=rfieldb(i,k,j)+rfieldp(i,k,j)
  2757. ENDDO
  2758. ENDDO
  2759. ENDDO
  2760. END SUBROUTINE calculate_full
  2761. !------------------------------------------------------------------------------
  2762. SUBROUTINE coriolis ( ru, rv, rw, ru_tend, rv_tend, rw_tend, &
  2763. config_flags, &
  2764. msftx, msfty, msfux, msfuy, &
  2765. msfvx, msfvy, &
  2766. f, e, sina, cosa, fzm, fzp, &
  2767. ids, ide, jds, jde, kds, kde, &
  2768. ims, ime, jms, jme, kms, kme, &
  2769. its, ite, jts, jte, kts, kte )
  2770. IMPLICIT NONE
  2771. ! Input data
  2772. TYPE(grid_config_rec_type) , INTENT(IN ) :: config_flags
  2773. INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
  2774. ims, ime, jms, jme, kms, kme, &
  2775. its, ite, jts, jte, kts, kte
  2776. REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: ru_tend, &
  2777. rv_tend, &
  2778. rw_tend
  2779. REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: ru, &
  2780. rv, &
  2781. rw
  2782. REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: msfux, &
  2783. msfuy, &
  2784. msfvx, &
  2785. msfvy, &
  2786. msftx, &
  2787. msfty
  2788. REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: f, &
  2789. e, &
  2790. sina, &
  2791. cosa
  2792. REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fzm, &
  2793. fzp
  2794. ! Local indices.
  2795. INTEGER :: i, j , k, ktf
  2796. INTEGER :: i_start, i_end, j_start, j_end
  2797. LOGICAL :: specified
  2798. !<DESCRIPTION>
  2799. !
  2800. ! coriolis calculates the large timestep tendency terms in the
  2801. ! u, v, and w momentum equations arise from the coriolis force.
  2802. !
  2803. !</DESCRIPTION>
  2804. specified = .false.
  2805. if(config_flags%specified .or. config_flags%nested) specified = .true.
  2806. ktf=MIN(kte,kde-1)
  2807. ! coriolis for u-momentum equation
  2808. ! Notes on map scale factor
  2809. ! cosa, sina are related to rotating the coordinate frame if desired
  2810. ! generally sina=0, cosa=1
  2811. ! ADT eqn 44, RHS terms 6 and 7: -2 mu w omega cos(lat)/my
  2812. ! + 2 mu v omega sin(lat)/my
  2813. ! Define f=2 omega sin(lat), e=2 omega cos(lat)
  2814. ! => terms are: -e mu w / my + f mu v / my
  2815. ! rv = mu v / mx ; rw = mu w / my
  2816. ! => terms are: -e rw + f rv *mx / my
  2817. i_start = its
  2818. i_end = ite
  2819. IF ( config_flags%open_xs .or. specified .or. &
  2820. config_flags%nested) i_start = MAX(ids+1,its)
  2821. IF ( config_flags%open_xe .or. specified .or. &
  2822. config_flags%nested) i_end = MIN(ide-1,ite)
  2823. IF ( config_flags%periodic_x ) i_start = its
  2824. IF ( config_flags%periodic_x ) i_end = ite
  2825. DO j = jts, MIN(jte,jde-1)
  2826. DO k=kts,ktf
  2827. DO i = i_start, i_end
  2828. ru_tend(i,k,j)=ru_tend(i,k,j) + (msfux(i,j)/msfuy(i,j))*0.5*(f(i,j)+f(i-1,j)) &
  2829. *0.25*(rv(i-1,k,j+1)+rv(i,k,j+1)+rv(i-1,k,j)+rv(i,k,j)) &
  2830. - 0.5*(e(i,j)+e(i-1,j))*0.5*(cosa(i,j)+cosa(i-1,j)) &
  2831. *0.25*(rw(i-1,k+1,j)+rw(i-1,k,j)+rw(i,k+1,j)+rw(i,k,j))
  2832. ENDDO
  2833. ENDDO
  2834. ! boundary loops for coriolis not needed for open bdy (commented out 20100611 JD)
  2835. ! IF ( (config_flags%open_xs) .and. (its == ids) ) THEN
  2836. ! DO k=kts,ktf
  2837. !
  2838. ! ru_tend(its,k,j)=ru_tend(its,k,j) + (msfux(its,j)/msfuy(its,j))*0.5*(f(its,j)+f(its,j)) &
  2839. ! *0.25*(rv(its,k,j+1)+rv(its,k,j+1)+rv(its,k,j)+rv(its,k,j)) &
  2840. ! - 0.5*(e(its,j)+e(its,j))*0.5*(cosa(its,j)+cosa(its,j)) &
  2841. ! *0.25*(rw(its,k+1,j)+rw(its,k,j)+rw(its,k+1,j)+rw(its,k,j))
  2842. ! ENDDO
  2843. ! ENDIF
  2844. ! IF ( (config_flags%open_xe) .and. (ite == ide) ) THEN
  2845. ! DO k=kts,ktf
  2846. !
  2847. ! ru_tend(ite,k,j)=ru_tend(ite,k,j) + (msfux(ite,j)/msfuy(ite,j))*0.5*(f(ite-1,j)+f(ite-1,j)) &
  2848. ! *0.25*(rv(ite-1,k,j+1)+rv(ite-1,k,j+1)+rv(ite-1,k,j)+rv(ite-1,k,j)) &
  2849. ! - 0.5*(e(ite-1,j)+e(ite-1,j))*0.5*(cosa(ite-1,j)+cosa(ite-1,j)) &
  2850. ! *0.25*(rw(ite-1,k+1,j)+rw(ite-1,k,j)+rw(ite-1,k+1,j)+rw(ite-1,k,j))
  2851. ! ENDDO
  2852. ! ENDIF
  2853. ENDDO
  2854. ! coriolis term for v-momentum equation
  2855. ! Notes on map scale factors
  2856. ! ADT eqn 45, RHS terms 6 and 6b [0 for sina=0]: -2 mu u omega sin(lat)/mx + ?
  2857. ! Define f=2 omega sin(lat), e=2 omega cos(lat)
  2858. ! => terms are: -f mu u / mx
  2859. ! ru = mu u / my ; rw = mu w / my
  2860. ! => terms are: -f ru *my / mx + ?
  2861. j_start = jts
  2862. j_end = jte
  2863. IF ( config_flags%open_ys .or. specified .or. &
  2864. config_flags%nested .or. config_flags%polar) j_start = MAX(jds+1,jts)
  2865. IF ( config_flags%open_ye .or. specified .or. &
  2866. config_flags%nested .or. config_flags%polar) j_end = MIN(jde-1,jte)
  2867. ! boundary loops for coriolis not needed for open bdy (commented out 20100611 JD)
  2868. ! IF ( (config_flags%open_ys) .and. (jts == jds) ) THEN
  2869. ! DO k=kts,ktf
  2870. ! DO i=its,MIN(ide-1,ite)
  2871. !
  2872. ! rv_tend(i,k,jts)=rv_tend(i,k,jts) - (msfvy(i,jts)/msfvx(i,jts))*0.5*(f(i,jts)+f(i,jts)) &
  2873. ! *0.25*(ru(i,k,jts)+ru(i+1,k,jts)+ru(i,k,jts)+ru(i+1,k,jts)) &
  2874. ! + (msfvy(i,jts)/msfvx(i,jts))*0.5*(e(i,jts)+e(i,jts))*0.5*(sina(i,jts)+sina(i,jts)) &
  2875. ! *0.25*(rw(i,k+1,jts)+rw(i,k,jts)+rw(i,k+1,jts)+rw(i,k,jts))
  2876. ! ENDDO
  2877. ! ENDDO
  2878. ! ENDIF
  2879. DO j=j_start, j_end
  2880. DO k=kts,ktf
  2881. DO i=its,MIN(ide-1,ite)
  2882. rv_tend(i,k,j)=rv_tend(i,k,j) - (msfvy(i,j)/msfvx(i,j))*0.5*(f(i,j)+f(i,j-1)) &
  2883. *0.25*(ru(i,k,j)+ru(i+1,k,j)+ru(i,k,j-1)+ru(i+1,k,j-1)) &
  2884. + (msfvy(i,j)/msfvx(i,j))*0.5*(e(i,j)+e(i,j-1))*0.5*(sina(i,j)+sina(i,j-1)) &
  2885. *0.25*(rw(i,k+1,j-1)+rw(i,k,j-1)+rw(i,k+1,j)+rw(i,k,j))
  2886. ENDDO
  2887. ENDDO
  2888. ENDDO
  2889. ! boundary loops for coriolis not needed for open bdy (commented out 20100611 JD)
  2890. ! IF ( (config_flags%open_ye) .and. (jte == jde) ) THEN
  2891. ! DO k=kts,ktf
  2892. ! DO i=its,MIN(ide-1,ite)
  2893. !
  2894. ! rv_tend(i,k,jte)=rv_tend(i,k,jte) - (msfvy(i,jte)/msfvx(i,jte))*0.5*(f(i,jte-1)+f(i,jte-1)) &
  2895. ! *0.25*(ru(i,k,jte-1)+ru(i+1,k,jte-1)+ru(i,k,jte-1)+ru(i+1,k,jte-1)) &
  2896. ! + (msfvy(i,jte)/msfvx(i,jte))*0.5*(e(i,jte-1)+e(i,jte-1))*0.5*(sina(i,jte-1)+sina(i,jte-1)) &
  2897. ! *0.25*(rw(i,k+1,jte-1)+rw(i,k,jte-1)+rw(i,k+1,jte-1)+rw(i,k,jte-1))
  2898. ! ENDDO
  2899. ! ENDDO
  2900. ! ENDIF
  2901. ! coriolis term for w-mometum
  2902. ! Notes on map scale factors
  2903. ! ADT eqn 46/my, RHS terms 5 and 5b [0 for sina=0]: 2 mu u omega cos(lat)/my +?
  2904. ! Define e=2 omega cos(lat)
  2905. ! => terms are: e mu u / my + ???
  2906. ! ru = mu u / my ; ru = mu v / mx
  2907. ! => terms are: e ru + ???
  2908. DO j=jts,MIN(jte, jde-1)
  2909. DO k=kts+1,ktf
  2910. DO i=its,MIN(ite, ide-1)
  2911. rw_tend(i,k,j)=rw_tend(i,k,j) + e(i,j)* &
  2912. (cosa(i,j)*0.5*(fzm(k)*(ru(i,k,j)+ru(i+1,k,j)) &
  2913. +fzp(k)*(ru(i,k-1,j)+ru(i+1,k-1,j))) &
  2914. -(msftx(i,j)/msfty(i,j))* &
  2915. sina(i,j)*0.5*(fzm(k)*(rv(i,k,j)+rv(i,k,j+1)) &
  2916. +fzp(k)*(rv(i,k-1,j)+rv(i,k-1,j+1))))
  2917. ENDDO
  2918. ENDDO
  2919. ENDDO
  2920. END SUBROUTINE coriolis
  2921. !------------------------------------------------------------------------------
  2922. SUBROUTINE perturbation_coriolis ( ru_in, rv_in, rw, ru_tend, rv_tend, rw_tend, &
  2923. config_flags, &
  2924. u_base, v_base, z_base, &
  2925. muu, muv, phb, ph, &
  2926. msftx, msfty, msfux, msfuy, msfvx, msfvy, &
  2927. f, e, sina, cosa, fzm, fzp, &
  2928. ids, ide, jds, jde, kds, kde, &
  2929. ims, ime, jms, jme, kms, kme, &
  2930. its, ite, jts, jte, kts, kte )
  2931. IMPLICIT NONE
  2932. ! Input data
  2933. TYPE(grid_config_rec_type) , INTENT(IN ) :: config_flags
  2934. INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
  2935. ims, ime, jms, jme, kms, kme, &
  2936. its, ite, jts, jte, kts, kte
  2937. REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: ru_tend, &
  2938. rv_tend, &
  2939. rw_tend
  2940. REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: ru_in, &
  2941. rv_in, &
  2942. rw, &
  2943. ph, &
  2944. phb
  2945. REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: msfux, &
  2946. msfuy, &
  2947. msfvx, &
  2948. msfvy, &
  2949. msftx, &
  2950. msfty
  2951. REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: f, &
  2952. e, &
  2953. sina, &
  2954. cosa
  2955. REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: muu, &
  2956. muv
  2957. REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fzm, &
  2958. fzp
  2959. REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: u_base, &
  2960. v_base, &
  2961. z_base
  2962. ! Local storage
  2963. REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) :: ru, &
  2964. rv
  2965. REAL :: z_at_u, z_at_v, wkp1, wk, wkm1
  2966. ! Local indices.
  2967. INTEGER :: i, j , k, ktf
  2968. INTEGER :: i_start, i_end, j_start, j_end
  2969. LOGICAL :: specified
  2970. !<DESCRIPTION>
  2971. !
  2972. ! perturbation_coriolis calculates the large timestep tendency terms in the
  2973. ! u, v, and w momentum equations arise from the coriolis force. This version
  2974. ! subtracts off the horizontal velocities from the initial sounding when
  2975. ! computing the forcing terms, hence "perturbation" coriolis.
  2976. !
  2977. !</DESCRIPTION>
  2978. specified = .false.
  2979. if(config_flags%specified .or. config_flags%nested) specified = .true.
  2980. ktf=MIN(kte,kde-1)
  2981. ! coriolis for u-momentum equation
  2982. i_start = its
  2983. i_end = ite
  2984. IF ( config_flags%open_xs .or. specified .or. &
  2985. config_flags%nested) i_start = MAX(ids+1,its)
  2986. IF ( config_flags%open_xe .or. specified .or. &
  2987. config_flags%nested) i_end = MIN(ide-1,ite)
  2988. IF ( config_flags%periodic_x ) i_start = its
  2989. IF ( config_flags%periodic_x ) i_end = ite
  2990. ! compute perturbation mu*v for use in u momentum equation
  2991. DO j = jts, MIN(jte,jde-1)+1
  2992. DO k=kts+1,ktf-1
  2993. DO i = i_start-1, i_end
  2994. z_at_v = 0.25*( phb(i,k,j )+phb(i,k+1,j ) &
  2995. +phb(i,k,j-1)+phb(i,k+1,j-1) &
  2996. +ph(i,k,j )+ph(i,k+1,j ) &
  2997. +ph(i,k,j-1)+ph(i,k+1,j-1))/g
  2998. wkp1 = min(1.,max(0.,z_at_v-z_base(k))/(z_base(k+1)-z_base(k)))
  2999. wkm1 = min(1.,max(0.,z_base(k)-z_at_v)/(z_base(k)-z_base(k-1)))
  3000. wk = 1.-wkp1-wkm1
  3001. rv(i,k,j) = rv_in(i,k,j) - muv(i,j)*( &
  3002. wkm1*v_base(k-1) &
  3003. +wk *v_base(k ) &
  3004. +wkp1*v_base(k+1) )
  3005. ENDDO
  3006. ENDDO
  3007. ENDDO
  3008. ! pick up top and bottom v
  3009. DO j = jts, MIN(jte,jde-1)+1
  3010. DO i = i_start-1, i_end
  3011. k = kts
  3012. z_at_v = 0.25*( phb(i,k,j )+phb(i,k+1,j ) &
  3013. +phb(i,k,j-1)+phb(i,k+1,j-1) &
  3014. +ph(i,k,j )+ph(i,k+1,j ) &
  3015. +ph(i,k,j-1)+ph(i,k+1,j-1))/g
  3016. wkp1 = min(1.,max(0.,z_at_v-z_base(k))/(z_base(k+1)-z_base(k)))
  3017. wk = 1.-wkp1
  3018. rv(i,k,j) = rv_in(i,k,j) - muv(i,j)*( &
  3019. +wk *v_base(k ) &
  3020. +wkp1*v_base(k+1) )
  3021. k = ktf
  3022. z_at_v = 0.25*( phb(i,k,j )+phb(i,k+1,j ) &
  3023. +phb(i,k,j-1)+phb(i,k+1,j-1) &
  3024. +ph(i,k,j )+ph(i,k+1,j ) &
  3025. +ph(i,k,j-1)+ph(i,k+1,j-1))/g
  3026. wkm1 = min(1.,max(0.,z_base(k)-z_at_v)/(z_base(k)-z_base(k-1)))
  3027. wk = 1.-wkm1
  3028. rv(i,k,j) = rv_in(i,k,j) - muv(i,j)*( &
  3029. wkm1*v_base(k-1) &
  3030. +wk *v_base(k ) )
  3031. ENDDO
  3032. ENDDO
  3033. ! compute coriolis forcing for u
  3034. ! Map scale factors: see comments above for Coriolis
  3035. DO j = jts, MIN(jte,jde-1)
  3036. DO k=kts,ktf
  3037. DO i = i_start, i_end
  3038. ru_tend(i,k,j)=ru_tend(i,k,j) + (msfux(i,j)/msfuy(i,j))*0.5*(f(i,j)+f(i-1,j)) &
  3039. *0.25*(rv(i-1,k,j+1)+rv(i,k,j+1)+rv(i-1,k,j)+rv(i,k,j)) &
  3040. - 0.5*(e(i,j)+e(i-1,j))*0.5*(cosa(i,j)+cosa(i-1,j)) &
  3041. *0.25*(rw(i-1,k+1,j)+rw(i-1,k,j)+rw(i,k+1,j)+rw(i,k,j))
  3042. ENDDO
  3043. ENDDO
  3044. ! boundary loops for perturbation coriolis is needed for open bdy (20110225 JD)
  3045. IF ( (config_flags%open_xs) .and. (its == ids) ) THEN
  3046. DO k=kts,ktf
  3047. ru_tend(its,k,j)=ru_tend(its,k,j) + (msfux(its,j)/msfuy(its,j))*0.5*(f(its,j)+f(its,j)) &
  3048. *0.25*(rv(its,k,j+1)+rv(its,k,j+1)+rv(its,k,j)+rv(its,k,j)) &
  3049. - 0.5*(e(its,j)+e(its,j))*0.5*(cosa(its,j)+cosa(its,j)) &
  3050. *0.25*(rw(its,k+1,j)+rw(its,k,j)+rw(its,k+1,j)+rw(its,k,j))
  3051. ENDDO
  3052. ENDIF
  3053. IF ( (config_flags%open_xe) .and. (ite == ide) ) THEN
  3054. DO k=kts,ktf
  3055. ru_tend(ite,k,j)=ru_tend(ite,k,j) + (msfux(ite,j)/msfuy(ite,j))*0.5*(f(ite-1,j)+f(ite-1,j)) &
  3056. *0.25*(rv(ite-1,k,j+1)+rv(ite-1,k,j+1)+rv(ite-1,k,j)+rv(ite-1,k,j)) &
  3057. - 0.5*(e(ite-1,j)+e(ite-1,j))*0.5*(cosa(ite-1,j)+cosa(ite-1,j)) &
  3058. *0.25*(rw(ite-1,k+1,j)+rw(ite-1,k,j)+rw(ite-1,k+1,j)+rw(ite-1,k,j))
  3059. ENDDO
  3060. ENDIF
  3061. ENDDO
  3062. ! coriolis term for v-momentum equation
  3063. ! Map scale factors: see comments above for Coriolis
  3064. j_start = jts
  3065. j_end = jte
  3066. IF ( config_flags%open_ys .or. specified .or. &
  3067. config_flags%nested .or. config_flags%polar) j_start = MAX(jds+1,jts)
  3068. IF ( config_flags%open_ye .or. specified .or. &
  3069. config_flags%nested .or. config_flags%polar) j_end = MIN(jde-1,jte)
  3070. ! compute perturbation mu*u for use in v momentum equation
  3071. DO j = j_start-1,j_end
  3072. DO k=kts+1,ktf-1
  3073. DO i = its, MIN(ite,ide-1)+1
  3074. z_at_u = 0.25*( phb(i ,k,j)+phb(i ,k+1,j) &
  3075. +phb(i-1,k,j)+phb(i-1,k+1,j) &
  3076. +ph(i ,k,j)+ph(i ,k+1,j) &
  3077. +ph(i-1,k,j)+ph(i-1,k+1,j))/g
  3078. wkp1 = min(1.,max(0.,z_at_u-z_base(k))/(z_base(k+1)-z_base(k)))
  3079. wkm1 = min(1.,max(0.,z_base(k)-z_at_u)/(z_base(k)-z_base(k-1)))
  3080. wk = 1.-wkp1-wkm1
  3081. ru(i,k,j) = ru_in(i,k,j) - muu(i,j)*( &
  3082. wkm1*u_base(k-1) &
  3083. +wk *u_base(k ) &
  3084. +wkp1*u_base(k+1) )
  3085. ENDDO
  3086. ENDDO
  3087. ENDDO
  3088. ! pick up top and bottom u
  3089. DO j = j_start-1,j_end
  3090. DO i = its, MIN(ite,ide-1)+1
  3091. k = kts
  3092. z_at_u = 0.25*( phb(i ,k,j)+phb(i ,k+1,j) &
  3093. +phb(i-1,k,j)+phb(i-1,k+1,j) &
  3094. +ph(i ,k,j)+ph(i ,k+1,j) &
  3095. +ph(i-1,k,j)+ph(i-1,k+1,j))/g
  3096. wkp1 = min(1.,max(0.,z_at_u-z_base(k))/(z_base(k+1)-z_base(k)))
  3097. wk = 1.-wkp1
  3098. ru(i,k,j) = ru_in(i,k,j) - muu(i,j)*( &
  3099. +wk *u_base(k ) &
  3100. +wkp1*u_base(k+1) )
  3101. k = ktf
  3102. z_at_u = 0.25*( phb(i ,k,j)+phb(i ,k+1,j) &
  3103. +phb(i-1,k,j)+phb(i-1,k+1,j) &
  3104. +ph(i ,k,j)+ph(i ,k+1,j) &
  3105. +ph(i-1,k,j)+ph(i-1,k+1,j))/g
  3106. wkm1 = min(1.,max(0.,z_base(k)-z_at_u)/(z_base(k)-z_base(k-1)))
  3107. wk = 1.-wkm1
  3108. ru(i,k,j) = ru_in(i,k,j) - muu(i,j)*( &
  3109. wkm1*u_base(k-1) &
  3110. +wk *u_base(k ) )
  3111. ENDDO
  3112. ENDDO
  3113. ! compute coriolis forcing for v momentum equation
  3114. ! Map scale factors: see comments above for Coriolis
  3115. ! boundary loops for perturbation coriolis is needed for open bdy (20110225 JD)
  3116. IF ( (config_flags%open_ys) .and. (jts == jds) ) THEN
  3117. DO k=kts,ktf
  3118. DO i=its,MIN(ide-1,ite)
  3119. rv_tend(i,k,jts)=rv_tend(i,k,jts) - (msfvy(i,jts)/msfvx(i,jts))*0.5*(f(i,jts)+f(i,jts)) &
  3120. *0.25*(ru(i,k,jts)+ru(i+1,k,jts)+ru(i,k,jts)+ru(i+1,k,jts)) &
  3121. + (msfvy(i,jts)/msfvx(i,jts))*0.5*(e(i,jts)+e(i,jts))*0.5*(sina(i,jts)+sina(i,jts)) &
  3122. *0.25*(rw(i,k+1,jts)+rw(i,k,jts)+rw(i,k+1,jts)+rw(i,k,jts))
  3123. ENDDO
  3124. ENDDO
  3125. ENDIF
  3126. DO j=j_start, j_end
  3127. DO k=kts,ktf
  3128. DO i=its,MIN(ide-1,ite)
  3129. rv_tend(i,k,j)=rv_tend(i,k,j) - (msfvy(i,j)/msfvx(i,j))*0.5*(f(i,j)+f(i,j-1)) &
  3130. *0.25*(ru(i,k,j)+ru(i+1,k,j)+ru(i,k,j-1)+ru(i+1,k,j-1)) &
  3131. + (msfvy(i,j)/msfvx(i,j))*0.5*(e(i,j)+e(i,j-1))*0.5*(sina(i,j)+sina(i,j-1)) &
  3132. *0.25*(rw(i,k+1,j-1)+rw(i,k,j-1)+rw(i,k+1,j)+rw(i,k,j))
  3133. ENDDO
  3134. ENDDO
  3135. ENDDO
  3136. ! boundary loops for perturbation coriolis is needed for open bdy (20110225 JD)
  3137. IF ( (config_flags%open_ye) .and. (jte == jde) ) THEN
  3138. DO k=kts,ktf
  3139. DO i=its,MIN(ide-1,ite)
  3140. rv_tend(i,k,jte)=rv_tend(i,k,jte) - (msfvy(i,jte)/msfvx(i,jte))*0.5*(f(i,jte-1)+f(i,jte-1)) &
  3141. *0.25*(ru(i,k,jte-1)+ru(i+1,k,jte-1)+ru(i,k,jte-1)+ru(i+1,k,jte-1)) &
  3142. + (msfvy(i,jte)/msfvx(i,jte))*0.5*(e(i,jte-1)+e(i,jte-1))*0.5*(sina(i,jte-1)+sina(i,jte-1)) &
  3143. *0.25*(rw(i,k+1,jte-1)+rw(i,k,jte-1)+rw(i,k+1,jte-1)+rw(i,k,jte-1))
  3144. ENDDO
  3145. ENDDO
  3146. ENDIF
  3147. ! coriolis term for w-mometum
  3148. ! Map scale factors: see comments above for Coriolis
  3149. DO j=jts,MIN(jte, jde-1)
  3150. DO k=kts+1,ktf
  3151. DO i=its,MIN(ite, ide-1)
  3152. rw_tend(i,k,j)=rw_tend(i,k,j) + e(i,j)* &
  3153. (cosa(i,j)*0.5*(fzm(k)*(ru(i,k,j)+ru(i+1,k,j)) &
  3154. +fzp(k)*(ru(i,k-1,j)+ru(i+1,k-1,j))) &
  3155. -(msftx(i,j)/msfty(i,j))*sina(i,j)*0.5*(fzm(k)*(rv(i,k,j)+rv(i,k,j+1)) &
  3156. +fzp(k)*(rv(i,k-1,j)+rv(i,k-1,j+1))))
  3157. ENDDO
  3158. ENDDO
  3159. ENDDO
  3160. END SUBROUTINE perturbation_coriolis
  3161. !------------------------------------------------------------------------------
  3162. SUBROUTINE curvature ( ru, rv, rw, u, v, w, ru_tend, rv_tend, rw_tend, &
  3163. config_flags, &
  3164. msfux, msfuy, msfvx, msfvy, msftx, msfty, &
  3165. xlat, fzm, fzp, rdx, rdy, &
  3166. ids, ide, jds, jde, kds, kde, &
  3167. ims, ime, jms, jme, kms, kme, &
  3168. its, ite, jts, jte, kts, kte )
  3169. IMPLICIT NONE
  3170. ! Input data
  3171. TYPE(grid_config_rec_type) , INTENT(IN ) :: config_flags
  3172. INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
  3173. ims, ime, jms, jme, kms, kme, &
  3174. its, ite, jts, jte, kts, kte
  3175. REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , &
  3176. INTENT(INOUT) :: ru_tend, &
  3177. rv_tend, &
  3178. rw_tend
  3179. REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , &
  3180. INTENT(IN ) :: ru, &
  3181. rv, &
  3182. rw, &
  3183. u, &
  3184. v, &
  3185. w
  3186. REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: msfux, &
  3187. msfuy, &
  3188. msfvx, &
  3189. msfvy, &
  3190. msftx, &
  3191. msfty, &
  3192. xlat
  3193. REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fzm, &
  3194. fzp
  3195. REAL , INTENT(IN ) :: rdx, &
  3196. rdy
  3197. ! Local data
  3198. ! INTEGER :: i, j, k, itf, jtf, ktf, kp1, im, ip, jm, jp
  3199. INTEGER :: i, j, k, itf, jtf, ktf
  3200. INTEGER :: i_start, i_end, j_start, j_end
  3201. ! INTEGER :: irmin, irmax, jrmin, jrmax
  3202. REAL , DIMENSION( its-1:ite , kts:kte, jts-1:jte ) :: vxgm
  3203. LOGICAL :: specified
  3204. !<DESCRIPTION>
  3205. !
  3206. ! curvature calculates the large timestep tendency terms in the
  3207. ! u, v, and w momentum equations arise from the curvature terms.
  3208. !
  3209. !</DESCRIPTION>
  3210. specified = .false.
  3211. if(config_flags%specified .or. config_flags%nested) specified = .true.
  3212. itf=MIN(ite,ide-1)
  3213. jtf=MIN(jte,jde-1)
  3214. ktf=MIN(kte,kde-1)
  3215. ! irmin = ims
  3216. ! irmax = ime
  3217. ! jrmin = jms
  3218. ! jrmax = jme
  3219. ! IF ( config_flags%open_xs ) irmin = ids
  3220. ! IF ( config_flags%open_xe ) irmax = ide-1
  3221. ! IF ( config_flags%open_ys ) jrmin = jds
  3222. ! IF ( config_flags%open_ye ) jrmax = jde-1
  3223. ! Define v cross grad m at scalar points - vxgm(i,j)
  3224. i_start = its-1
  3225. i_end = ite
  3226. j_start = jts-1
  3227. j_end = jte
  3228. IF ( ( config_flags%open_xs .or. specified .or. &
  3229. config_flags%nested) .and. (its == ids) ) i_start = its
  3230. IF ( ( config_flags%open_xe .or. specified .or. &
  3231. config_flags%nested) .and. (ite == ide) ) i_end = ite-1
  3232. IF ( ( config_flags%open_ys .or. specified .or. &
  3233. config_flags%nested .or. config_flags%polar) .and. (jts == jds) ) j_start = jts
  3234. IF ( ( config_flags%open_ye .or. specified .or. &
  3235. config_flags%nested .or. config_flags%polar) .and. (jte == jde) ) j_end = jte-1
  3236. IF ( config_flags%periodic_x ) i_start = its-1
  3237. IF ( config_flags%periodic_x ) i_end = ite
  3238. DO j=j_start, j_end
  3239. DO k=kts,ktf
  3240. DO i=i_start, i_end
  3241. ! Map scale factor notes:
  3242. ! msf...y is constant everywhere for cylindrical map projection
  3243. ! msf...x varies with y only
  3244. ! But we know that this is not = 0 for cylindrical,
  3245. ! therefore use msfvX in 1st line
  3246. ! which => by symmetry use msfuY in 2nd line - ???
  3247. vxgm(i,k,j)=0.5*(u(i,k,j)+u(i+1,k,j))*(msfvx(i,j+1)-msfvx(i,j))*rdy - &
  3248. 0.5*(v(i,k,j)+v(i,k,j+1))*(msfuy(i+1,j)-msfuy(i,j))*rdx
  3249. ENDDO
  3250. ENDDO
  3251. ENDDO
  3252. ! Pick up the boundary rows for open (radiation) lateral b.c.
  3253. ! Rather crude at present, we are assuming there is no
  3254. ! variation in this term at the boundary.
  3255. IF ( ( config_flags%open_xs .or. (specified .AND. .NOT. config_flags%periodic_x) .or. &
  3256. config_flags%nested) .and. (its == ids) ) THEN
  3257. DO j = jts, jte-1
  3258. DO k = kts, ktf
  3259. vxgm(its-1,k,j) = vxgm(its,k,j)
  3260. ENDDO
  3261. ENDDO
  3262. ENDIF
  3263. IF ( ( config_flags%open_xe .or. (specified .AND. .NOT. config_flags%periodic_x) .or. &
  3264. config_flags%nested) .and. (ite == ide) ) THEN
  3265. DO j = jts, jte-1
  3266. DO k = kts, ktf
  3267. vxgm(ite,k,j) = vxgm(ite-1,k,j)
  3268. ENDDO
  3269. ENDDO
  3270. ENDIF
  3271. ! Polar boundary condition:
  3272. ! The following change is needed in case one tries using the vxgm route with
  3273. ! polar B.C.'s in the future, but not needed if 'tan' used
  3274. IF ( ( config_flags%open_ys .or. specified .or. &
  3275. config_flags%nested .or. config_flags%polar) .and. (jts == jds) ) THEN
  3276. DO k = kts, ktf
  3277. DO i = its-1, ite
  3278. vxgm(i,k,jts-1) = vxgm(i,k,jts)
  3279. ENDDO
  3280. ENDDO
  3281. ENDIF
  3282. ! Polar boundary condition:
  3283. ! The following change is needed in case one tries using the vxgm route with
  3284. ! polar B.C.'s in the future, but not needed if 'tan' used
  3285. IF ( ( config_flags%open_ye .or. specified .or. &
  3286. config_flags%nested .or. config_flags%polar) .and. (jte == jde) ) THEN
  3287. DO k = kts, ktf
  3288. DO i = its-1, ite
  3289. vxgm(i,k,jte) = vxgm(i,k,jte-1)
  3290. ENDDO
  3291. ENDDO
  3292. ENDIF
  3293. ! curvature term for u momentum eqn.
  3294. ! Map scale factor notes:
  3295. ! ADT eqn 44, RHS terms 4 and 5, in cylindrical: mu u v tan(lat)/(a my)
  3296. ! - mu u w /(a my)
  3297. ! ru = mu u / my ; rw = mu w / my ; rv = mu v / mx
  3298. ! => terms are:
  3299. ! (mx/my)*u rv tan(lat) / a - u rw / a = (u/a)*[(mx/my) rv tan(lat) - rw]
  3300. ! ru v tan(lat) / a - u rw / a
  3301. ! xlat defined with end points half grid space from pole,
  3302. ! hence are on u latitude points
  3303. i_start = its
  3304. IF ( config_flags%open_xs .or. specified .or. &
  3305. config_flags%nested) i_start = MAX ( ids+1 , its )
  3306. IF ( config_flags%open_xe .or. specified .or. &
  3307. config_flags%nested) i_end = MIN ( ide-1 , ite )
  3308. IF ( config_flags%periodic_x ) i_start = its
  3309. IF ( config_flags%periodic_x ) i_end = ite
  3310. ! Polar boundary condition
  3311. IF ((config_flags%map_proj == 6) .OR. (config_flags%polar)) THEN
  3312. DO j=jts,MIN(jde-1,jte)
  3313. DO k=kts,ktf
  3314. DO i=i_start,i_end
  3315. ru_tend(i,k,j)=ru_tend(i,k,j) + u(i,k,j)*reradius* ( &
  3316. (msfux(i,j)/msfuy(i,j))*0.25*(rv(i-1,k,j+1)+rv(i,k,j+1)+ &
  3317. rv(i-1,k,j)+rv(i,k,j))*tan(xlat(i,j)*degrad) &
  3318. - 0.25*(rw(i-1,k+1,j)+rw(i-1,k,j)+rw(i,k+1,j)+rw(i,k,j)) )
  3319. ENDDO
  3320. ENDDO
  3321. ENDDO
  3322. ELSE ! normal code
  3323. DO j=jts,MIN(jde-1,jte)
  3324. DO k=kts,ktf
  3325. DO i=i_start,i_end
  3326. ru_tend(i,k,j)=ru_tend(i,k,j) + 0.5*(vxgm(i,k,j)+vxgm(i-1,k,j)) &
  3327. *0.25*(rv(i-1,k,j+1)+rv(i,k,j+1)+rv(i-1,k,j)+rv(i,k,j)) &
  3328. - u(i,k,j)*reradius &
  3329. *0.25*(rw(i-1,k+1,j)+rw(i-1,k,j)+rw(i,k+1,j)+rw(i,k,j))
  3330. ENDDO
  3331. ENDDO
  3332. ENDDO
  3333. END IF
  3334. ! curvature term for v momentum eqn.
  3335. ! Map scale factor notes
  3336. ! ADT eqn 45, RHS terms 4 and 5, in cylindrical: - mu u*u tan(lat)/(a mx)
  3337. ! - mu v w /(a mx)
  3338. ! ru = mu u / my ; rw = mu w / my ; rv = mu v / mx
  3339. ! terms are:
  3340. ! - (my/mx)*u ru tan(lat) / a - (my/mx)*v rw / a
  3341. ! = - [my/(mx*a)]*[u ru tan(lat) + v rw]
  3342. ! - (1/a)*[(my/mx)*u ru tan(lat) + w rv]
  3343. ! xlat defined with end points half grid space from pole, hence are on
  3344. ! u latitude points => av here
  3345. !
  3346. ! in original wrf, there was a sign error for the rw contribution
  3347. j_start = jts
  3348. IF ( config_flags%open_ys .or. specified .or. &
  3349. config_flags%nested .or. config_flags%polar) j_start = MAX ( jds+1 , jts )
  3350. IF ( config_flags%open_ye .or. specified .or. &
  3351. config_flags%nested .or. config_flags%polar) j_end = MIN ( jde-1 , jte )
  3352. IF ((config_flags%map_proj == 6) .OR. (config_flags%polar)) THEN
  3353. DO j=j_start,j_end
  3354. DO k=kts,ktf
  3355. DO i=its,MIN(ite,ide-1)
  3356. rv_tend(i,k,j)=rv_tend(i,k,j) - (msfvy(i,j)/msfvx(i,j))*reradius* ( &
  3357. 0.25*(u(i,k,j)+u(i+1,k,j)+u(i,k,j-1)+u(i+1,k,j-1))* &
  3358. tan((xlat(i,j)+xlat(i,j-1))*0.5*degrad)* &
  3359. 0.25*(ru(i,k,j)+ru(i+1,k,j)+ru(i,k,j-1)+ru(i+1,k,j-1)) &
  3360. + v(i,k,j)*0.25*(rw(i,k+1,j-1)+rw(i,k,j-1)+ &
  3361. rw(i,k+1,j)+rw(i,k,j)) )
  3362. ENDDO
  3363. ENDDO
  3364. ENDDO
  3365. ELSE ! normal code
  3366. DO j=j_start,j_end
  3367. DO k=kts,ktf
  3368. DO i=its,MIN(ite,ide-1)
  3369. rv_tend(i,k,j)=rv_tend(i,k,j) - 0.5*(vxgm(i,k,j)+vxgm(i,k,j-1)) &
  3370. *0.25*(ru(i,k,j)+ru(i+1,k,j)+ru(i,k,j-1)+ru(i+1,k,j-1)) &
  3371. - (msfvy(i,j)/msfvx(i,j))*v(i,k,j)*reradius &
  3372. *0.25*(rw(i,k+1,j-1)+rw(i,k,j-1)+rw(i,k+1,j)+rw(i,k,j))
  3373. ENDDO
  3374. ENDDO
  3375. ENDDO
  3376. END IF
  3377. ! curvature term for vertical momentum eqn.
  3378. ! Notes on map scale factors:
  3379. ! ADT eqn 46, RHS term 4: [mu/(a my)]*[u*u + v*v]
  3380. ! ru = mu u / my ; rw = mu w / my ; rv = mu v / mx
  3381. ! terms are: u ru / a + (mx/my)v rv / a
  3382. DO j=jts,MIN(jte,jde-1)
  3383. DO k=MAX(2,kts),ktf
  3384. DO i=its,MIN(ite,ide-1)
  3385. rw_tend(i,k,j)=rw_tend(i,k,j) + reradius* &
  3386. (0.5*(fzm(k)*(ru(i,k,j)+ru(i+1,k,j))+fzp(k)*(ru(i,k-1,j)+ru(i+1,k-1,j))) &
  3387. *0.5*(fzm(k)*( u(i,k,j) +u(i+1,k,j))+fzp(k)*( u(i,k-1,j) +u(i+1,k-1,j))) &
  3388. +(msftx(i,j)/msfty(i,j))*0.5*(fzm(k)*(rv(i,k,j)+rv(i,k,j+1))+fzp(k)*(rv(i,k-1,j)+rv(i,k-1,j+1))) &
  3389. *0.5*(fzm(k)*( v(i,k,j) +v(i,k,j+1))+fzp(k)*( v(i,k-1,j) +v(i,k-1,j+1))))
  3390. ENDDO
  3391. ENDDO
  3392. ENDDO
  3393. END SUBROUTINE curvature
  3394. !------------------------------------------------------------------------------
  3395. SUBROUTINE decouple ( rr, rfield, field, name, config_flags, &
  3396. fzm, fzp, &
  3397. ids, ide, jds, jde, kds, kde, &
  3398. ims, ime, jms, jme, kms, kme, &
  3399. its, ite, jts, jte, kts, kte )
  3400. IMPLICIT NONE
  3401. ! Input data
  3402. TYPE(grid_config_rec_type) , INTENT(IN ) :: config_flags
  3403. INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
  3404. ims, ime, jms, jme, kms, kme, &
  3405. its, ite, jts, jte, kts, kte
  3406. CHARACTER(LEN=1) , INTENT(IN ) :: name
  3407. REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: rfield
  3408. REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: rr
  3409. REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT ) :: field
  3410. REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fzm, fzp
  3411. ! Local data
  3412. INTEGER :: i, j, k, itf, jtf, ktf
  3413. !<DESCRIPTION>
  3414. !
  3415. ! decouple decouples a variable from the column dry-air mass.
  3416. !
  3417. !</DESCRIPTION>
  3418. ktf=MIN(kte,kde-1)
  3419. IF (name .EQ. 'u')THEN
  3420. itf=ite
  3421. jtf=MIN(jte,jde-1)
  3422. DO j=jts,jtf
  3423. DO k=kts,ktf
  3424. DO i=its,itf
  3425. field(i,k,j)=rfield(i,k,j)/(0.5*(rr(i,k,j)+rr(i-1,k,j)))
  3426. ENDDO
  3427. ENDDO
  3428. ENDDO
  3429. ELSE IF (name .EQ. 'v')THEN
  3430. itf=MIN(ite,ide-1)
  3431. jtf=jte
  3432. DO j=jts,jtf
  3433. DO k=kts,ktf
  3434. DO i=its,itf
  3435. field(i,k,j)=rfield(i,k,j)/(0.5*(rr(i,k,j)+rr(i,k,j-1)))
  3436. ENDDO
  3437. ENDDO
  3438. ENDDO
  3439. ELSE IF (name .EQ. 'w')THEN
  3440. itf=MIN(ite,ide-1)
  3441. jtf=MIN(jte,jde-1)
  3442. DO j=jts,jtf
  3443. DO k=kts+1,ktf
  3444. DO i=its,itf
  3445. field(i,k,j)=rfield(i,k,j)/(fzm(k)*rr(i,k,j)+fzp(k)*rr(i,k-1,j))
  3446. ENDDO
  3447. ENDDO
  3448. ENDDO
  3449. DO j=jts,jtf
  3450. DO i=its,itf
  3451. field(i,kte,j) = 0.
  3452. ENDDO
  3453. ENDDO
  3454. ELSE
  3455. itf=MIN(ite,ide-1)
  3456. jtf=MIN(jte,jde-1)
  3457. ! For theta we will decouple tb and tp and add them to give t afterwards
  3458. DO j=jts,jtf
  3459. DO k=kts,ktf
  3460. DO i=its,itf
  3461. field(i,k,j)=rfield(i,k,j)/rr(i,k,j)
  3462. ENDDO
  3463. ENDDO
  3464. ENDDO
  3465. ENDIF
  3466. END SUBROUTINE decouple
  3467. !-------------------------------------------------------------------------------
  3468. SUBROUTINE zero_tend ( tendency, &
  3469. ids, ide, jds, jde, kds, kde, &
  3470. ims, ime, jms, jme, kms, kme, &
  3471. its, ite, jts, jte, kts, kte )
  3472. IMPLICIT NONE
  3473. ! Input data
  3474. INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
  3475. ims, ime, jms, jme, kms, kme, &
  3476. its, ite, jts, jte, kts, kte
  3477. REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency
  3478. ! Local data
  3479. INTEGER :: i, j, k, itf, jtf, ktf
  3480. !<DESCRIPTION>
  3481. !
  3482. ! zero_tend sets the input tendency array to zero.
  3483. !
  3484. !</DESCRIPTION>
  3485. DO j = jts, jte
  3486. DO k = kts, kte
  3487. DO i = its, ite
  3488. tendency(i,k,j) = 0.
  3489. ENDDO
  3490. ENDDO
  3491. ENDDO
  3492. END SUBROUTINE zero_tend
  3493. !-------------------------------------------------------------------------------
  3494. ! Sets the an array on the polar v point(s) to zero
  3495. SUBROUTINE zero_pole ( field, &
  3496. ids, ide, jds, jde, kds, kde, &
  3497. ims, ime, jms, jme, kms, kme, &
  3498. its, ite, jts, jte, kts, kte )
  3499. IMPLICIT NONE
  3500. ! Input data
  3501. INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
  3502. ims, ime, jms, jme, kms, kme, &
  3503. its, ite, jts, jte, kts, kte
  3504. REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: field
  3505. ! Local data
  3506. INTEGER :: i, k
  3507. IF (jts == jds) THEN
  3508. DO k = kts, kte
  3509. DO i = its-1, ite+1
  3510. field(i,k,jts) = 0.
  3511. END DO
  3512. END DO
  3513. END IF
  3514. IF (jte == jde) THEN
  3515. DO k = kts, kte
  3516. DO i = its-1, ite+1
  3517. field(i,k,jte) = 0.
  3518. END DO
  3519. END DO
  3520. END IF
  3521. END SUBROUTINE zero_pole
  3522. !-------------------------------------------------------------------------------
  3523. ! Sets the an array on the polar v point(s)
  3524. SUBROUTINE pole_point_bc ( field, &
  3525. ids, ide, jds, jde, kds, kde, &
  3526. ims, ime, jms, jme, kms, kme, &
  3527. its, ite, jts, jte, kts, kte )
  3528. IMPLICIT NONE
  3529. ! Input data
  3530. INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
  3531. ims, ime, jms, jme, kms, kme, &
  3532. its, ite, jts, jte, kts, kte
  3533. REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: field
  3534. ! Local data
  3535. INTEGER :: i, k
  3536. IF (jts == jds) THEN
  3537. DO k = kts, kte
  3538. DO i = its, ite
  3539. ! field(i,k,jts) = 2*field(i,k,jts+1) - field(i,k,jts+2)
  3540. field(i,k,jts) = field(i,k,jts+1)
  3541. END DO
  3542. END DO
  3543. END IF
  3544. IF (jte == jde) THEN
  3545. DO k = kts, kte
  3546. DO i = its, ite
  3547. ! field(i,k,jte) = 2*field(i,k,jte-1) - field(i,k,jte-2)
  3548. field(i,k,jte) = field(i,k,jte-1)
  3549. END DO
  3550. END DO
  3551. END IF
  3552. END SUBROUTINE pole_point_bc
  3553. !======================================================================
  3554. ! physics prep routines
  3555. !======================================================================
  3556. SUBROUTINE phy_prep ( config_flags, & ! input
  3557. mu, muu, muv, u, v, p, pb, alt, ph, & ! input
  3558. phb, t, tsk, moist, n_moist, & ! input
  3559. rho, th_phy, p_phy , pi_phy , & ! output
  3560. u_phy, v_phy, p8w, t_phy, t8w, & ! output
  3561. z, z_at_w, dz8w, & ! output
  3562. p_hyd, p_hyd_w, dnw, & ! output
  3563. fzm, fzp, znw, p_top, & ! params
  3564. RTHRATEN, &
  3565. RTHBLTEN, RUBLTEN, RVBLTEN, &
  3566. RQVBLTEN, RQCBLTEN, RQIBLTEN, &
  3567. RUCUTEN, RVCUTEN, RTHCUTEN, &
  3568. RQVCUTEN, RQCCUTEN, RQRCUTEN, &
  3569. RQICUTEN, RQSCUTEN, &
  3570. RUSHTEN, RVSHTEN, RTHSHTEN, &
  3571. RQVSHTEN, RQCSHTEN, RQRSHTEN, &
  3572. RQISHTEN, RQSSHTEN, RQGSHTEN, &
  3573. RTHFTEN, RQVFTEN, &
  3574. RUNDGDTEN, RVNDGDTEN, RTHNDGDTEN, &
  3575. RPHNDGDTEN,RQVNDGDTEN, RMUNDGDTEN, &
  3576. ids, ide, jds, jde, kds, kde, &
  3577. ims, ime, jms, jme, kms, kme, &
  3578. its, ite, jts, jte, kts, kte )
  3579. !----------------------------------------------------------------------
  3580. IMPLICIT NONE
  3581. !----------------------------------------------------------------------
  3582. TYPE(grid_config_rec_type) , INTENT(IN ) :: config_flags
  3583. INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
  3584. ims, ime, jms, jme, kms, kme, &
  3585. its, ite, jts, jte, kts, kte
  3586. INTEGER , INTENT(IN ) :: n_moist
  3587. REAL, DIMENSION( ims:ime, kms:kme , jms:jme , n_moist ), INTENT(IN) :: moist
  3588. REAL , DIMENSION( ims:ime, jms:jme ), INTENT(IN ) :: TSK, mu, muu, muv
  3589. REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , &
  3590. INTENT( OUT) :: u_phy, &
  3591. v_phy, &
  3592. pi_phy, &
  3593. p_phy, &
  3594. p8w, &
  3595. t_phy, &
  3596. th_phy, &
  3597. t8w, &
  3598. rho, &
  3599. z, &
  3600. dz8w, &
  3601. z_at_w
  3602. REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , &
  3603. INTENT( OUT) :: p_hyd, &
  3604. p_hyd_w
  3605. REAL , INTENT(IN ) :: p_top
  3606. REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , &
  3607. INTENT(IN ) :: pb, &
  3608. p, &
  3609. u, &
  3610. v, &
  3611. alt, &
  3612. ph, &
  3613. phb, &
  3614. t
  3615. REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fzm, &
  3616. fzp
  3617. REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: znw, &
  3618. dnw
  3619. REAL, DIMENSION( ims:ime , kms:kme, jms:jme ), &
  3620. INTENT(INOUT) :: RTHRATEN
  3621. REAL, DIMENSION( ims:ime , kms:kme, jms:jme ), &
  3622. INTENT(INOUT) :: RUCUTEN, &
  3623. RVCUTEN, &
  3624. RTHCUTEN, &
  3625. RQVCUTEN, &
  3626. RQCCUTEN, &
  3627. RQRCUTEN, &
  3628. RQICUTEN, &
  3629. RQSCUTEN, &
  3630. RUSHTEN, &
  3631. RVSHTEN, &
  3632. RTHSHTEN, &
  3633. RQVSHTEN, &
  3634. RQCSHTEN, &
  3635. RQRSHTEN, &
  3636. RQISHTEN, &
  3637. RQSSHTEN, &
  3638. RQGSHTEN
  3639. REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) , &
  3640. INTENT(INOUT) :: RUBLTEN, &
  3641. RVBLTEN, &
  3642. RTHBLTEN, &
  3643. RQVBLTEN, &
  3644. RQCBLTEN, &
  3645. RQIBLTEN
  3646. REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) , &
  3647. INTENT(INOUT) :: RTHFTEN, &
  3648. RQVFTEN
  3649. REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) , &
  3650. INTENT(INOUT) :: RUNDGDTEN, &
  3651. RVNDGDTEN, &
  3652. RTHNDGDTEN, &
  3653. RPHNDGDTEN, &
  3654. RQVNDGDTEN
  3655. REAL, DIMENSION( ims:ime, jms:jme ) , &
  3656. INTENT(INOUT) :: RMUNDGDTEN
  3657. INTEGER :: i_start, i_end, j_start, j_end, k_start, k_end, i_startu, j_startv
  3658. INTEGER :: i, j, k
  3659. REAL :: w1, w2, z0, z1, z2
  3660. REAL :: qtot
  3661. INTEGER :: n
  3662. !-----------------------------------------------------------------------
  3663. !<DESCRIPTION>
  3664. !
  3665. ! phys_prep calculates a number of diagnostic quantities needed by
  3666. ! the physics routines. It also decouples the physics tendencies from
  3667. ! the column dry-air mass (the physics routines expect to see/update the
  3668. ! uncoupled tendencies).
  3669. !
  3670. !</DESCRIPTION>
  3671. ! set up loop bounds for this grid's boundary conditions
  3672. i_start = its
  3673. i_end = min( ite,ide-1 )
  3674. j_start = jts
  3675. j_end = min( jte,jde-1 )
  3676. k_start = kts
  3677. k_end = min( kte, kde-1 )
  3678. ! compute thermodynamics and velocities at pressure points (or half levels)
  3679. do j = j_start,j_end
  3680. do k = k_start, k_end
  3681. do i = i_start, i_end
  3682. th_phy(i,k,j) = t(i,k,j) + t0
  3683. p_phy(i,k,j) = p(i,k,j) + pb(i,k,j)
  3684. pi_phy(i,k,j) = (p_phy(i,k,j)/p1000mb)**rcp
  3685. t_phy(i,k,j) = th_phy(i,k,j)*pi_phy(i,k,j)
  3686. rho(i,k,j) = 1./alt(i,k,j)*(1.+moist(i,k,j,P_QV))
  3687. u_phy(i,k,j) = 0.5*(u(i,k,j)+u(i+1,k,j))
  3688. v_phy(i,k,j) = 0.5*(v(i,k,j)+v(i,k,j+1))
  3689. enddo
  3690. enddo
  3691. enddo
  3692. ! compute z at w points
  3693. do j = j_start,j_end
  3694. do k = k_start, kte
  3695. do i = i_start, i_end
  3696. z_at_w(i,k,j) = (phb(i,k,j)+ph(i,k,j))/g
  3697. enddo
  3698. enddo
  3699. enddo
  3700. do j = j_start,j_end
  3701. do k = k_start, kte-1
  3702. do i = i_start, i_end
  3703. dz8w(i,k,j) = z_at_w(i,k+1,j)-z_at_w(i,k,j)
  3704. enddo
  3705. enddo
  3706. enddo
  3707. do j = j_start,j_end
  3708. do i = i_start, i_end
  3709. dz8w(i,kte,j) = 0.
  3710. enddo
  3711. enddo
  3712. ! compute z at p points or half levels (average of z at full levels)
  3713. do j = j_start,j_end
  3714. do k = k_start, k_end
  3715. do i = i_start, i_end
  3716. z(i,k,j) = 0.5*(z_at_w(i,k,j) +z_at_w(i,k+1,j) )
  3717. enddo
  3718. enddo
  3719. enddo
  3720. ! interp t and p to full levels
  3721. do j = j_start,j_end
  3722. do k = 2, k_end
  3723. do i = i_start, i_end
  3724. p8w(i,k,j) = fzm(k)*p_phy(i,k,j)+fzp(k)*p_phy(i,k-1,j)
  3725. t8w(i,k,j) = fzm(k)*t_phy(i,k,j)+fzp(k)*t_phy(i,k-1,j)
  3726. enddo
  3727. enddo
  3728. enddo
  3729. ! extrapolate p and t to surface and top.
  3730. ! we'll use an extrapolation in z for now
  3731. do j = j_start,j_end
  3732. do i = i_start, i_end
  3733. ! bottom
  3734. z0 = z_at_w(i,1,j)
  3735. z1 = z(i,1,j)
  3736. z2 = z(i,2,j)
  3737. w1 = (z0 - z2)/(z1 - z2)
  3738. w2 = 1. - w1
  3739. p8w(i,1,j) = w1*p_phy(i,1,j)+w2*p_phy(i,2,j)
  3740. t8w(i,1,j) = w1*t_phy(i,1,j)+w2*t_phy(i,2,j)
  3741. ! top
  3742. z0 = z_at_w(i,kte,j)
  3743. z1 = z(i,k_end,j)
  3744. z2 = z(i,k_end-1,j)
  3745. w1 = (z0 - z2)/(z1 - z2)
  3746. w2 = 1. - w1
  3747. ! p8w(i,kde,j) = w1*p_phy(i,kde-1,j)+w2*p_phy(i,kde-2,j)
  3748. !!! bug fix extrapolate ln(p) so p is positive definite
  3749. p8w(i,kde,j) = exp(w1*log(p_phy(i,kde-1,j))+w2*log(p_phy(i,kde-2,j)))
  3750. t8w(i,kde,j) = w1*t_phy(i,kde-1,j)+w2*t_phy(i,kde-2,j)
  3751. enddo
  3752. enddo
  3753. ! calculate hydrostatic pressure at both full and half levels
  3754. ! first, full level p: assuming dry over model top
  3755. do j = j_start,j_end
  3756. do i = i_start, i_end
  3757. p_hyd_w(i,kte,j) = p_top
  3758. enddo
  3759. enddo
  3760. do j = j_start,j_end
  3761. do k = kte-1, k_start, -1
  3762. do i = i_start, i_end
  3763. qtot = 0.
  3764. do n = PARAM_FIRST_SCALAR,n_moist
  3765. qtot = qtot + moist(i,k,j,n)
  3766. enddo
  3767. p_hyd_w(i,k,j) = p_hyd_w(i,k+1,j) - (1.+qtot)*mu(i,j)*dnw(k)
  3768. ! p_hyd_w(i,k,j) = p_hyd_w(i,k+1,j)+1./alt(i,k,j)*(1.+moist(i,k,j,P_QV))*g*dz8w(i,k,j)
  3769. enddo
  3770. enddo
  3771. enddo
  3772. ! now calculate hydrostatic pressure at half levels
  3773. do j = j_start,j_end
  3774. do k = k_start, k_end
  3775. do i = i_start, i_end
  3776. p_hyd(i,k,j) = 0.5*(p_hyd_w(i,k,j)+p_hyd_w(i,k+1,j))
  3777. enddo
  3778. enddo
  3779. enddo
  3780. ! decouple all physics tendencies
  3781. IF (config_flags%ra_lw_physics .gt. 0 .or. config_flags%ra_sw_physics .gt. 0) THEN
  3782. DO J=j_start,j_end
  3783. DO K=k_start,k_end
  3784. DO I=i_start,i_end
  3785. RTHRATEN(I,K,J)=RTHRATEN(I,K,J)/mu(I,J)
  3786. ENDDO
  3787. ENDDO
  3788. ENDDO
  3789. ENDIF
  3790. IF (config_flags%cu_physics .gt. 0) THEN
  3791. DO J=j_start,j_end
  3792. DO I=i_start,i_end
  3793. DO K=k_start,k_end
  3794. RUCUTEN(I,K,J) =RUCUTEN(I,K,J)/mu(I,J)
  3795. RVCUTEN(I,K,J) =RVCUTEN(I,K,J)/mu(I,J)
  3796. RTHCUTEN(I,K,J)=RTHCUTEN(I,K,J)/mu(I,J)
  3797. ENDDO
  3798. ENDDO
  3799. ENDDO
  3800. IF (P_QV .ge. PARAM_FIRST_SCALAR)THEN
  3801. DO J=j_start,j_end
  3802. DO I=i_start,i_end
  3803. DO K=k_start,k_end
  3804. RQVCUTEN(I,K,J)=RQVCUTEN(I,K,J)/mu(I,J)
  3805. ENDDO
  3806. ENDDO
  3807. ENDDO
  3808. ENDIF
  3809. IF (P_QC .ge. PARAM_FIRST_SCALAR)THEN
  3810. DO J=j_start,j_end
  3811. DO I=i_start,i_end
  3812. DO K=k_start,k_end
  3813. RQCCUTEN(I,K,J)=RQCCUTEN(I,K,J)/mu(I,J)
  3814. ENDDO
  3815. ENDDO
  3816. ENDDO
  3817. ENDIF
  3818. IF (P_QR .ge. PARAM_FIRST_SCALAR)THEN
  3819. DO J=j_start,j_end
  3820. DO I=i_start,i_end
  3821. DO K=k_start,k_end
  3822. RQRCUTEN(I,K,J)=RQRCUTEN(I,K,J)/mu(I,J)
  3823. ENDDO
  3824. ENDDO
  3825. ENDDO
  3826. ENDIF
  3827. IF (P_QI .ge. PARAM_FIRST_SCALAR)THEN
  3828. DO J=j_start,j_end
  3829. DO I=i_start,i_end
  3830. DO K=k_start,k_end
  3831. RQICUTEN(I,K,J)=RQICUTEN(I,K,J)/mu(I,J)
  3832. ENDDO
  3833. ENDDO
  3834. ENDDO
  3835. ENDIF
  3836. IF(P_QS .ge. PARAM_FIRST_SCALAR)THEN
  3837. DO J=j_start,j_end
  3838. DO I=i_start,i_end
  3839. DO K=k_start,k_end
  3840. RQSCUTEN(I,K,J)=RQSCUTEN(I,K,J)/mu(I,J)
  3841. ENDDO
  3842. ENDDO
  3843. ENDDO
  3844. ENDIF
  3845. ENDIF
  3846. IF (config_flags%shcu_physics .gt. 0) THEN
  3847. DO J=j_start,j_end
  3848. DO I=i_start,i_end
  3849. DO K=k_start,k_end
  3850. RUSHTEN(I,K,J) =RUSHTEN(I,K,J)/mu(I,J)
  3851. RVSHTEN(I,K,J) =RVSHTEN(I,K,J)/mu(I,J)
  3852. RTHSHTEN(I,K,J)=RTHSHTEN(I,K,J)/mu(I,J)
  3853. ENDDO
  3854. ENDDO
  3855. ENDDO
  3856. IF (P_QV .ge. PARAM_FIRST_SCALAR)THEN
  3857. DO J=j_start,j_end
  3858. DO I=i_start,i_end
  3859. DO K=k_start,k_end
  3860. RQVSHTEN(I,K,J)=RQVSHTEN(I,K,J)/mu(I,J)
  3861. ENDDO
  3862. ENDDO
  3863. ENDDO
  3864. ENDIF
  3865. IF (P_QC .ge. PARAM_FIRST_SCALAR)THEN
  3866. DO J=j_start,j_end
  3867. DO I=i_start,i_end
  3868. DO K=k_start,k_end
  3869. RQCSHTEN(I,K,J)=RQCSHTEN(I,K,J)/mu(I,J)
  3870. ENDDO
  3871. ENDDO
  3872. ENDDO
  3873. ENDIF
  3874. IF (P_QR .ge. PARAM_FIRST_SCALAR)THEN
  3875. DO J=j_start,j_end
  3876. DO I=i_start,i_end
  3877. DO K=k_start,k_end
  3878. RQRSHTEN(I,K,J)=RQRSHTEN(I,K,J)/mu(I,J)
  3879. ENDDO
  3880. ENDDO
  3881. ENDDO
  3882. ENDIF
  3883. IF (P_QI .ge. PARAM_FIRST_SCALAR)THEN
  3884. DO J=j_start,j_end
  3885. DO I=i_start,i_end
  3886. DO K=k_start,k_end
  3887. RQISHTEN(I,K,J)=RQISHTEN(I,K,J)/mu(I,J)
  3888. ENDDO
  3889. ENDDO
  3890. ENDDO
  3891. ENDIF
  3892. IF(P_QS .ge. PARAM_FIRST_SCALAR)THEN
  3893. DO J=j_start,j_end
  3894. DO I=i_start,i_end
  3895. DO K=k_start,k_end
  3896. RQSSHTEN(I,K,J)=RQSSHTEN(I,K,J)/mu(I,J)
  3897. ENDDO
  3898. ENDDO
  3899. ENDDO
  3900. ENDIF
  3901. IF(P_QG .ge. PARAM_FIRST_SCALAR)THEN
  3902. DO J=j_start,j_end
  3903. DO I=i_start,i_end
  3904. DO K=k_start,k_end
  3905. RQGSHTEN(I,K,J)=RQGSHTEN(I,K,J)/mu(I,J)
  3906. ENDDO
  3907. ENDDO
  3908. ENDDO
  3909. ENDIF
  3910. ENDIF
  3911. IF (config_flags%bl_pbl_physics .gt. 0) THEN
  3912. DO J=j_start,j_end
  3913. DO K=k_start,k_end
  3914. DO I=i_start,i_end
  3915. RUBLTEN(I,K,J) =RUBLTEN(I,K,J)/mu(I,J)
  3916. RVBLTEN(I,K,J) =RVBLTEN(I,K,J)/mu(I,J)
  3917. RTHBLTEN(I,K,J)=RTHBLTEN(I,K,J)/mu(I,J)
  3918. ENDDO
  3919. ENDDO
  3920. ENDDO
  3921. IF (P_QV .ge. PARAM_FIRST_SCALAR) THEN
  3922. DO J=j_start,j_end
  3923. DO K=k_start,k_end
  3924. DO I=i_start,i_end
  3925. RQVBLTEN(I,K,J)=RQVBLTEN(I,K,J)/mu(I,J)
  3926. ENDDO
  3927. ENDDO
  3928. ENDDO
  3929. ENDIF
  3930. IF (P_QC .ge. PARAM_FIRST_SCALAR) THEN
  3931. DO J=j_start,j_end
  3932. DO K=k_start,k_end
  3933. DO I=i_start,i_end
  3934. RQCBLTEN(I,K,J)=RQCBLTEN(I,K,J)/mu(I,J)
  3935. ENDDO
  3936. ENDDO
  3937. ENDDO
  3938. ENDIF
  3939. IF (P_QI .ge. PARAM_FIRST_SCALAR) THEN
  3940. DO J=j_start,j_end
  3941. DO K=k_start,k_end
  3942. DO I=i_start,i_end
  3943. RQIBLTEN(I,K,J)=RQIBLTEN(I,K,J)/mu(I,J)
  3944. ENDDO
  3945. ENDDO
  3946. ENDDO
  3947. ENDIF
  3948. ENDIF
  3949. ! decouple advective forcing required by Grell-Devenyi scheme
  3950. if(( config_flags%cu_physics == GDSCHEME ) .OR. &
  3951. ( config_flags%cu_physics == G3SCHEME ) .OR. &
  3952. ( config_flags%cu_physics == KFETASCHEME ) .OR. &
  3953. ( config_flags%cu_physics == TIEDTKESCHEME ) ) then ! Tiedtke ZCX&YQW
  3954. DO J=j_start,j_end
  3955. DO I=i_start,i_end
  3956. DO K=k_start,k_end
  3957. RTHFTEN(I,K,J)=RTHFTEN(I,K,J)/mu(I,J)
  3958. ENDDO
  3959. ENDDO
  3960. ENDDO
  3961. IF (P_QV .ge. PARAM_FIRST_SCALAR)THEN
  3962. DO J=j_start,j_end
  3963. DO I=i_start,i_end
  3964. DO K=k_start,k_end
  3965. RQVFTEN(I,K,J)=RQVFTEN(I,K,J)/mu(I,J)
  3966. ENDDO
  3967. ENDDO
  3968. ENDDO
  3969. ENDIF
  3970. END IF
  3971. ! fdda
  3972. ! note fdda u and v tendencies are staggered, also only interior points have muu/muv,
  3973. ! so only decouple those
  3974. IF (config_flags%grid_fdda .gt. 0) THEN
  3975. i_startu=MAX(its,ids+1)
  3976. j_startv=MAX(jts,jds+1)
  3977. DO J=j_start,j_end
  3978. DO K=k_start,k_end
  3979. DO I=i_startu,i_end
  3980. RUNDGDTEN(I,K,J) =RUNDGDTEN(I,K,J)/muu(I,J)
  3981. ENDDO
  3982. ENDDO
  3983. ENDDO
  3984. DO J=j_startv,j_end
  3985. DO K=k_start,k_end
  3986. DO I=i_start,i_end
  3987. RVNDGDTEN(I,K,J) =RVNDGDTEN(I,K,J)/muv(I,J)
  3988. ENDDO
  3989. ENDDO
  3990. ENDDO
  3991. DO J=j_start,j_end
  3992. DO K=k_start,k_end
  3993. DO I=i_start,i_end
  3994. RTHNDGDTEN(I,K,J)=RTHNDGDTEN(I,K,J)/mu(I,J)
  3995. ! RMUNDGDTEN(I,J) - no coupling
  3996. ENDDO
  3997. ENDDO
  3998. ENDDO
  3999. IF (config_flags%grid_fdda .EQ. 2) THEN
  4000. DO J=j_start,j_end
  4001. DO K=k_start,kte
  4002. DO I=i_start,i_end
  4003. RPHNDGDTEN(I,K,J)=RPHNDGDTEN(I,K,J)/mu(I,J)
  4004. ENDDO
  4005. ENDDO
  4006. ENDDO
  4007. ELSE IF (config_flags%grid_fdda .EQ. 1) THEN
  4008. IF (P_QV .ge. PARAM_FIRST_SCALAR) THEN
  4009. DO J=j_start,j_end
  4010. DO K=k_start,k_end
  4011. DO I=i_start,i_end
  4012. RQVNDGDTEN(I,K,J)=RQVNDGDTEN(I,K,J)/mu(I,J)
  4013. ENDDO
  4014. ENDDO
  4015. ENDDO
  4016. ENDIF
  4017. ENDIF
  4018. ENDIF
  4019. END SUBROUTINE phy_prep
  4020. !------------------------------------------------------------
  4021. SUBROUTINE moist_physics_prep_em( t_new, t_old, t0, rho, al, alb, &
  4022. p, p8w, p0, pb, ph, phb, &
  4023. th_phy, pii, pf, &
  4024. z, z_at_w, dz8w, &
  4025. dt,h_diabatic, &
  4026. config_flags,fzm, fzp, &
  4027. ids,ide, jds,jde, kds,kde, &
  4028. ims,ime, jms,jme, kms,kme, &
  4029. its,ite, jts,jte, kts,kte )
  4030. IMPLICIT NONE
  4031. ! Here we construct full fields
  4032. ! needed by the microphysics
  4033. TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
  4034. INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde
  4035. INTEGER, INTENT(IN ) :: ims,ime, jms,jme, kms,kme
  4036. INTEGER, INTENT(IN ) :: its,ite, jts,jte, kts,kte
  4037. REAL, INTENT(IN ) :: dt
  4038. REAL, DIMENSION( ims:ime , kms:kme, jms:jme ), &
  4039. INTENT(IN ) :: al, &
  4040. alb, &
  4041. p, &
  4042. pb, &
  4043. ph, &
  4044. phb
  4045. REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fzm, &
  4046. fzp
  4047. REAL, DIMENSION( ims:ime , kms:kme, jms:jme ), &
  4048. INTENT( OUT) :: rho, &
  4049. th_phy, &
  4050. pii, &
  4051. pf, &
  4052. z, &
  4053. z_at_w, &
  4054. dz8w, &
  4055. p8w
  4056. REAL, DIMENSION( ims:ime , kms:kme, jms:jme ), &
  4057. INTENT(INOUT) :: h_diabatic
  4058. REAL, DIMENSION( ims:ime , kms:kme, jms:jme ), &
  4059. INTENT(INOUT) :: t_new, &
  4060. t_old
  4061. REAL, INTENT(IN ) :: t0, p0
  4062. REAL :: z0,z1,z2,w1,w2
  4063. INTEGER :: i_start, i_end, j_start, j_end, k_start, k_end
  4064. INTEGER :: i, j, k
  4065. !--------------------------------------------------------------------
  4066. !<DESCRIPTION>
  4067. !
  4068. ! moist_phys_prep_em calculates a number of diagnostic quantities needed by
  4069. ! the microphysics routines.
  4070. !
  4071. !</DESCRIPTION>
  4072. ! set up loop bounds for this grid's boundary conditions
  4073. i_start = its
  4074. i_end = min( ite,ide-1 )
  4075. j_start = jts
  4076. j_end = min( jte,jde-1 )
  4077. k_start = kts
  4078. k_end = min( kte, kde-1 )
  4079. DO j = j_start, j_end
  4080. DO k = k_start, kte
  4081. DO i = i_start, i_end
  4082. z_at_w(i,k,j) = (ph(i,k,j)+phb(i,k,j))/g
  4083. ENDDO
  4084. ENDDO
  4085. ENDDO
  4086. do j = j_start,j_end
  4087. do k = k_start, kte-1
  4088. do i = i_start, i_end
  4089. dz8w(i,k,j) = z_at_w(i,k+1,j)-z_at_w(i,k,j)
  4090. enddo
  4091. enddo
  4092. enddo
  4093. do j = j_start,j_end
  4094. do i = i_start, i_end
  4095. dz8w(i,kte,j) = 0.
  4096. enddo
  4097. enddo
  4098. ! compute full pii, rho, and z at the new time-level
  4099. ! (needed for physics).
  4100. ! convert perturbation theta to full theta (th_phy)
  4101. ! use h_diabatic to temporarily save pre-microphysics full theta
  4102. DO j = j_start, j_end
  4103. DO k = k_start, k_end
  4104. DO i = i_start, i_end
  4105. #ifdef REVERT
  4106. t_new(i,k,j) = t_new(i,k,j)-h_diabatic(i,k,j)*dt
  4107. #endif
  4108. th_phy(i,k,j) = t_new(i,k,j) + t0
  4109. h_diabatic(i,k,j) = th_phy(i,k,j)
  4110. rho(i,k,j) = 1./(al(i,k,j)+alb(i,k,j))
  4111. pii(i,k,j) = ((p(i,k,j)+pb(i,k,j))/p0)**rcp
  4112. z(i,k,j) = 0.5*(z_at_w(i,k,j) +z_at_w(i,k+1,j) )
  4113. pf(i,k,j) = p(i,k,j)+pb(i,k,j)
  4114. ENDDO
  4115. ENDDO
  4116. ENDDO
  4117. ! interp t and p at w points
  4118. do j = j_start,j_end
  4119. do k = 2, k_end
  4120. do i = i_start, i_end
  4121. p8w(i,k,j) = fzm(k)*pf(i,k,j)+fzp(k)*pf(i,k-1,j)
  4122. enddo
  4123. enddo
  4124. enddo
  4125. ! extrapolate p and t to surface and top.
  4126. ! we'll use an extrapolation in z for now
  4127. do j = j_start,j_end
  4128. do i = i_start, i_end
  4129. ! bottom
  4130. z0 = z_at_w(i,1,j)
  4131. z1 = z(i,1,j)
  4132. z2 = z(i,2,j)
  4133. w1 = (z0 - z2)/(z1 - z2)
  4134. w2 = 1. - w1
  4135. p8w(i,1,j) = w1*pf(i,1,j)+w2*pf(i,2,j)
  4136. ! top
  4137. z0 = z_at_w(i,kte,j)
  4138. z1 = z(i,k_end,j)
  4139. z2 = z(i,k_end-1,j)
  4140. w1 = (z0 - z2)/(z1 - z2)
  4141. w2 = 1. - w1
  4142. ! p8w(i,kde,j) = w1*pf(i,kde-1,j)+w2*pf(i,kde-2,j)
  4143. p8w(i,kde,j) = exp(w1*log(pf(i,kde-1,j))+w2*log(pf(i,kde-2,j)))
  4144. enddo
  4145. enddo
  4146. END SUBROUTINE moist_physics_prep_em
  4147. !------------------------------------------------------------------------------
  4148. SUBROUTINE moist_physics_finish_em( t_new, t_old, t0, mut, &
  4149. th_phy, h_diabatic, dt, &
  4150. config_flags, &
  4151. #if ( WRF_DFI_RADAR == 1 )
  4152. dfi_tten_rad,dfi_stage, &
  4153. #endif
  4154. ids,ide, jds,jde, kds,kde, &
  4155. ims,ime, jms,jme, kms,kme, &
  4156. its,ite, jts,jte, kts,kte )
  4157. IMPLICIT NONE
  4158. ! Here we construct full fields
  4159. ! needed by the microphysics
  4160. TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
  4161. INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde
  4162. INTEGER, INTENT(IN ) :: ims,ime, jms,jme, kms,kme
  4163. INTEGER, INTENT(IN ) :: its,ite, jts,jte, kts,kte
  4164. REAL, DIMENSION( ims:ime , kms:kme, jms:jme ), &
  4165. INTENT(INOUT) :: t_new, &
  4166. t_old, &
  4167. th_phy, &
  4168. h_diabatic
  4169. #if ( WRF_DFI_RADAR == 1 )
  4170. REAL, DIMENSION( ims:ime , kms:kme, jms:jme ), &
  4171. INTENT(IN), OPTIONAL :: dfi_tten_rad
  4172. INTEGER, INTENT(IN ) ,OPTIONAL :: dfi_stage
  4173. REAL :: dfi_tten_max, old_max
  4174. #endif
  4175. REAL mpten, mptenmax, mptenmin
  4176. REAL, DIMENSION( ims:ime , jms:jme ), INTENT(INOUT) :: mut
  4177. REAL, INTENT(IN ) :: t0, dt
  4178. INTEGER :: i_start, i_end, j_start, j_end, k_start, k_end
  4179. INTEGER :: i, j, k, imax, jmax, imin, jmin
  4180. !--------------------------------------------------------------------
  4181. !<DESCRIPTION>
  4182. !
  4183. ! moist_phys_finish_em resets theta to its perturbation value and
  4184. ! computes and stores the microphysics diabatic heating term.
  4185. !
  4186. !</DESCRIPTION>
  4187. ! set up loop bounds for this grid's boundary conditions
  4188. i_start = its
  4189. i_end = min( ite,ide-1 )
  4190. j_start = jts
  4191. j_end = min( jte,jde-1 )
  4192. ! i_start=max(its,ids+4)
  4193. ! i_end=min(ite,ide-5)
  4194. ! j_start=max(jts,jds+4)
  4195. ! j_end=min(jte,jde-5)
  4196. k_start = kts
  4197. k_end = min( kte, kde-1 )
  4198. #if ( WRF_DFI_RADAR == 1 )
  4199. IF ( PRESENT(dfi_stage) .and. PRESENT(dfi_tten_rad) ) THEN
  4200. IF ( dfi_stage ==DFI_FWD ) THEN
  4201. !$OMP CRITICAL(big_step_utilities)
  4202. WRITE(wrf_err_message,*)'Add radar tendency: i_start,j_start: ', i_start, j_start
  4203. CALL wrf_debug ( 100 , TRIM(wrf_err_message) )
  4204. !$OMP END CRITICAL(big_step_utilities)
  4205. ENDIF
  4206. ENDIF
  4207. dfi_tten_max=-999
  4208. old_max=-999
  4209. #endif
  4210. ! add microphysics theta diff to perturbation theta, set h_diabatic
  4211. IF ( config_flags%no_mp_heating .eq. 0 ) THEN
  4212. mptenmax = 0.
  4213. mptenmin = 999.
  4214. DO j = j_start, j_end
  4215. DO k = k_start, k_end
  4216. DO i = i_start, i_end
  4217. mpten = th_phy(i,k,j)-h_diabatic(i,k,j)
  4218. #if ( WRF_DFI_RADAR == 1 )
  4219. if(mpten.gt.mptenmax) then
  4220. mptenmax=mpten
  4221. imax=i
  4222. jmax=j
  4223. endif
  4224. if(mpten.lt.mptenmin) then
  4225. mptenmin=mpten
  4226. imin=i
  4227. jmin=j
  4228. endif
  4229. mpten=min(config_flags%mp_tend_lim*dt, mpten)
  4230. mpten=max(-config_flags%mp_tend_lim*dt, mpten)
  4231. if(k < k_end ) then
  4232. if(dfi_tten_max < dfi_tten_rad(i,k,j) ) dfi_tten_max = dfi_tten_rad(i,k,j)
  4233. if(old_max < (th_phy(i,k,j)-h_diabatic(i,k,j)) ) old_max=th_phy(i,k,j)-h_diabatic(i,k,j)
  4234. endif
  4235. IF ( PRESENT(dfi_stage) .and. PRESENT(dfi_tten_rad) ) THEN
  4236. IF ( dfi_stage == DFI_FWD .and. dfi_tten_rad(i,k,j) >= -0.1 .and. &
  4237. dfi_tten_rad(i,k,j) <= 0.1 .and. k < k_end ) THEN
  4238. ! add radar temp tendency
  4239. ! there is radar coverage
  4240. t_new(i,k,j) = t_new(i,k,j) + (dfi_tten_rad(i,k,j))*dt
  4241. ELSE
  4242. ! no radar coverage
  4243. t_new(i,k,j) = t_new(i,k,j) + mpten
  4244. ENDIF
  4245. ENDIF
  4246. #else
  4247. t_new(i,k,j) = t_new(i,k,j) + mpten
  4248. #endif
  4249. h_diabatic(i,k,j) = mpten/dt
  4250. ENDDO
  4251. ENDDO
  4252. ENDDO
  4253. ELSE
  4254. DO j = j_start, j_end
  4255. DO k = k_start, k_end
  4256. DO i = i_start, i_end
  4257. ! t_new(i,k,j) = t_new(i,k,j)
  4258. h_diabatic(i,k,j) = 0.
  4259. ENDDO
  4260. ENDDO
  4261. ENDDO
  4262. ENDIF
  4263. END SUBROUTINE moist_physics_finish_em
  4264. !----------------------------------------------------------------
  4265. SUBROUTINE init_module_big_step
  4266. END SUBROUTINE init_module_big_step
  4267. SUBROUTINE set_tend ( field, field_adv_tend, msf, &
  4268. ids, ide, jds, jde, kds, kde, &
  4269. ims, ime, jms, jme, kms, kme, &
  4270. its, ite, jts, jte, kts, kte )
  4271. IMPLICIT NONE
  4272. ! Input data
  4273. INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
  4274. ims, ime, jms, jme, kms, kme, &
  4275. its, ite, jts, jte, kts, kte
  4276. REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT) :: field
  4277. REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN) :: field_adv_tend
  4278. REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN) :: msf
  4279. ! Local data
  4280. INTEGER :: i, j, k, itf, jtf, ktf
  4281. !<DESCRIPTION>
  4282. !
  4283. ! set_tend copies the advective tendency array into the tendency array.
  4284. !
  4285. !</DESCRIPTION>
  4286. jtf = MIN(jte,jde-1)
  4287. ktf = MIN(kte,kde-1)
  4288. itf = MIN(ite,ide-1)
  4289. DO j = jts, jtf
  4290. DO k = kts, ktf
  4291. DO i = its, itf
  4292. field(i,k,j) = field_adv_tend(i,k,j)*msf(i,j)
  4293. ENDDO
  4294. ENDDO
  4295. ENDDO
  4296. END SUBROUTINE set_tend
  4297. !------------------------------------------------------------------------------
  4298. SUBROUTINE rk_rayleigh_damp( ru_tendf, rv_tendf, &
  4299. rw_tendf, t_tendf, &
  4300. u, v, w, t, t_init, &
  4301. mut, muu, muv, ph, phb, &
  4302. u_base, v_base, t_base, z_base, &
  4303. dampcoef, zdamp, &
  4304. ids, ide, jds, jde, kds, kde, &
  4305. ims, ime, jms, jme, kms, kme, &
  4306. its, ite, jts, jte, kts, kte )
  4307. ! History: Apr 2005 Modifications by George Bryan, NCAR:
  4308. ! - Generalized the code in a way that allows for
  4309. ! simulations with steep terrain.
  4310. !
  4311. ! Jul 2004 Modifications by George Bryan, NCAR:
  4312. ! - Modified the code to use u_base, v_base, and t_base
  4313. ! arrays for the background state. Removed the hard-wired
  4314. ! base-state values.
  4315. ! - Modified the code to use dampcoef, zdamp, and damp_opt,
  4316. ! i.e., the upper-level damper variables in namelist.input.
  4317. ! Removed the hard-wired variables in the older version.
  4318. ! This damper is used when damp_opt = 2.
  4319. ! - Modified the code to account for the movement of the
  4320. ! model surfaces with time. The code now obtains a base-
  4321. ! state value by interpolation using the "_base" arrays.
  4322. ! Nov 2003 Bug fix by Jason Knievel, NCAR
  4323. ! Aug 2003 Meridional dimension, some comments, and
  4324. ! changes in layout of the code added by
  4325. ! Jason Knievel, NCAR
  4326. ! Jul 2003 Original code by Bill Skamarock, NCAR
  4327. ! Purpose: This routine applies Rayleigh damping to a layer at top
  4328. ! of the model domain.
  4329. !-----------------------------------------------------------------------
  4330. ! Begin declarations.
  4331. IMPLICIT NONE
  4332. INTEGER, INTENT( IN ) &
  4333. :: ids, ide, jds, jde, kds, kde, &
  4334. ims, ime, jms, jme, kms, kme, &
  4335. its, ite, jts, jte, kts, kte
  4336. REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( INOUT ) &
  4337. :: ru_tendf, rv_tendf, rw_tendf, t_tendf
  4338. REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( IN ) &
  4339. :: u, v, w, t, t_init, ph, phb
  4340. REAL, DIMENSION( ims:ime, jms:jme ), INTENT( IN ) &
  4341. :: mut, muu, muv
  4342. REAL, DIMENSION( kms:kme ) , INTENT(IN ) &
  4343. :: u_base, v_base, t_base, z_base
  4344. REAL, INTENT(IN ) &
  4345. :: dampcoef, zdamp
  4346. ! Local variables.
  4347. INTEGER &
  4348. :: i_start, i_end, j_start, j_end, k_start, k_end, i, j, k, ktf, k1, k2
  4349. REAL &
  4350. :: pii, dcoef, z, ztop
  4351. REAL :: wkp1, wk, wkm1
  4352. REAL, DIMENSION( kms:kme ) :: z00, u00, v00, t00
  4353. ! End declarations.
  4354. !-----------------------------------------------------------------------
  4355. pii = 2.0 * asin(1.0)
  4356. ktf = MIN( kte, kde-1 )
  4357. !-----------------------------------------------------------------------
  4358. ! Adjust u to base state.
  4359. DO j = jts, MIN( jte, jde-1 )
  4360. DO i = its, MIN( ite, ide )
  4361. ! Get height at top of model
  4362. ztop = 0.5*( phb(i ,kde,j)+phb(i-1,kde,j) &
  4363. +ph(i ,kde,j)+ph(i-1,kde,j) )/g
  4364. ! Find bottom of damping layer
  4365. k1 = ktf
  4366. z = ztop
  4367. DO WHILE( z >= (ztop-zdamp) )
  4368. z = 0.25*( phb(i ,k1,j)+phb(i ,k1+1,j) &
  4369. +phb(i-1,k1,j)+phb(i-1,k1+1,j) &
  4370. +ph(i ,k1,j)+ph(i ,k1+1,j) &
  4371. +ph(i-1,k1,j)+ph(i-1,k1+1,j))/g
  4372. z00(k1) = z
  4373. k1 = k1 - 1
  4374. ENDDO
  4375. k1 = k1 + 2
  4376. ! Get reference state at model levels
  4377. DO k = k1, ktf
  4378. k2 = ktf
  4379. DO WHILE( z_base(k2) .gt. z00(k) )
  4380. k2 = k2 - 1
  4381. ENDDO
  4382. if(k2+1.gt.ktf)then
  4383. u00(k) = u_base(k2) + ( u_base(k2) - u_base(k2-1) ) &
  4384. * ( z00(k) - z_base(k2) ) &
  4385. / ( z_base(k2) - z_base(k2-1) )
  4386. else
  4387. u00(k) = u_base(k2) + ( u_base(k2+1) - u_base(k2) ) &
  4388. * ( z00(k) - z_base(k2) ) &
  4389. / ( z_base(k2+1) - z_base(k2) )
  4390. endif
  4391. ENDDO
  4392. ! Apply the Rayleigh damper
  4393. DO k = k1, ktf
  4394. dcoef = 1.0 - MIN( 1.0, ( ztop - z00(k) ) / zdamp )
  4395. dcoef = (SIN( 0.5 * pii * dcoef ) )**2
  4396. ru_tendf(i,k,j) = ru_tendf(i,k,j) - &
  4397. muu(i,j) * ( dcoef * dampcoef ) * &
  4398. ( u(i,k,j) - u00(k) )
  4399. END DO
  4400. END DO
  4401. END DO
  4402. ! End adjustment of u.
  4403. !-----------------------------------------------------------------------
  4404. !-----------------------------------------------------------------------
  4405. ! Adjust v to base state.
  4406. DO j = jts, MIN( jte, jde )
  4407. DO i = its, MIN( ite, ide-1 )
  4408. ! Get height at top of model
  4409. ztop = 0.5*( phb(i,kde,j )+phb(i,kde,j-1) &
  4410. +ph(i,kde,j )+ph(i,kde,j-1) )/g
  4411. ! Find bottom of damping layer
  4412. k1 = ktf
  4413. z = ztop
  4414. DO WHILE( z >= (ztop-zdamp) )
  4415. z = 0.25*( phb(i,k1,j )+phb(i,k1+1,j ) &
  4416. +phb(i,k1,j-1)+phb(i,k1+1,j-1) &
  4417. +ph(i,k1,j )+ph(i,k1+1,j ) &
  4418. +ph(i,k1,j-1)+ph(i,k1+1,j-1))/g
  4419. z00(k1) = z
  4420. k1 = k1 - 1
  4421. ENDDO
  4422. k1 = k1 + 2
  4423. ! Get reference state at model levels
  4424. DO k = k1, ktf
  4425. k2 = ktf
  4426. DO WHILE( z_base(k2) .gt. z00(k) )
  4427. k2 = k2 - 1
  4428. ENDDO
  4429. if(k2+1.gt.ktf)then
  4430. v00(k) = v_base(k2) + ( v_base(k2) - v_base(k2-1) ) &
  4431. * ( z00(k) - z_base(k2) ) &
  4432. / ( z_base(k2) - z_base(k2-1) )
  4433. else
  4434. v00(k) = v_base(k2) + ( v_base(k2+1) - v_base(k2) ) &
  4435. * ( z00(k) - z_base(k2) ) &
  4436. / ( z_base(k2+1) - z_base(k2) )
  4437. endif
  4438. ENDDO
  4439. ! Apply the Rayleigh damper
  4440. DO k = k1, ktf
  4441. dcoef = 1.0 - MIN( 1.0, ( ztop - z00(k) ) / zdamp )
  4442. dcoef = (SIN( 0.5 * pii * dcoef ) )**2
  4443. rv_tendf(i,k,j) = rv_tendf(i,k,j) - &
  4444. muv(i,j) * ( dcoef * dampcoef ) * &
  4445. ( v(i,k,j) - v00(k) )
  4446. END DO
  4447. END DO
  4448. END DO
  4449. ! End adjustment of v.
  4450. !-----------------------------------------------------------------------
  4451. !-----------------------------------------------------------------------
  4452. ! Adjust w to base state.
  4453. DO j = jts, MIN( jte, jde-1 )
  4454. DO i = its, MIN( ite, ide-1 )
  4455. ztop = ( phb(i,kde,j) + ph(i,kde,j) ) / g
  4456. DO k = kts, MIN( kte, kde )
  4457. z = ( phb(i,k,j) + ph(i,k,j) ) / g
  4458. IF ( z >= (ztop-zdamp) ) THEN
  4459. dcoef = 1.0 - MIN( 1.0, ( ztop - z ) / zdamp )
  4460. dcoef = ( SIN( 0.5 * pii * dcoef ) )**2
  4461. rw_tendf(i,k,j) = rw_tendf(i,k,j) - &
  4462. mut(i,j) * ( dcoef * dampcoef ) * w(i,k,j)
  4463. END IF
  4464. END DO
  4465. END DO
  4466. END DO
  4467. ! End adjustment of w.
  4468. !-----------------------------------------------------------------------
  4469. !-----------------------------------------------------------------------
  4470. ! Adjust potential temperature to base state.
  4471. DO j = jts, MIN( jte, jde-1 )
  4472. DO i = its, MIN( ite, ide-1 )
  4473. ! Get height at top of model
  4474. ztop = ( phb(i,kde,j) + ph(i,kde,j) ) / g
  4475. ! Find bottom of damping layer
  4476. k1 = ktf
  4477. z = ztop
  4478. DO WHILE( z >= (ztop-zdamp) )
  4479. z = 0.5 * ( phb(i,k1,j) + phb(i,k1+1,j) + &
  4480. ph(i,k1,j) + ph(i,k1+1,j) ) / g
  4481. z00(k1) = z
  4482. k1 = k1 - 1
  4483. ENDDO
  4484. k1 = k1 + 2
  4485. ! Get reference state at model levels
  4486. DO k = k1, ktf
  4487. k2 = ktf
  4488. DO WHILE( z_base(k2) .gt. z00(k) )
  4489. k2 = k2 - 1
  4490. ENDDO
  4491. if(k2+1.gt.ktf)then
  4492. t00(k) = t_base(k2) + ( t_base(k2) - t_base(k2-1) ) &
  4493. * ( z00(k) - z_base(k2) ) &
  4494. / ( z_base(k2) - z_base(k2-1) )
  4495. else
  4496. t00(k) = t_base(k2) + ( t_base(k2+1) - t_base(k2) ) &
  4497. * ( z00(k) - z_base(k2) ) &
  4498. / ( z_base(k2+1) - z_base(k2) )
  4499. endif
  4500. ENDDO
  4501. ! Apply the Rayleigh damper
  4502. DO k = k1, ktf
  4503. dcoef = 1.0 - MIN( 1.0, ( ztop - z00(k) ) / zdamp )
  4504. dcoef = (SIN( 0.5 * pii * dcoef ) )**2
  4505. t_tendf(i,k,j) = t_tendf(i,k,j) - &
  4506. mut(i,j) * ( dcoef * dampcoef ) * &
  4507. ( t(i,k,j) - t00(k) )
  4508. END DO
  4509. END DO
  4510. END DO
  4511. ! End adjustment of potential temperature.
  4512. !-----------------------------------------------------------------------
  4513. END SUBROUTINE rk_rayleigh_damp
  4514. !==============================================================================
  4515. !==============================================================================
  4516. SUBROUTINE theta_relaxation( t_tendf, t, t_init, &
  4517. mut, ph, phb, &
  4518. t_base, z_base, &
  4519. ids, ide, jds, jde, kds, kde, &
  4520. ims, ime, jms, jme, kms, kme, &
  4521. its, ite, jts, jte, kts, kte )
  4522. ! Purpose: Newtonian relaxation on potential temperature. Serves two
  4523. ! purposes: 1) to mimic atmospheric radiation in a simple
  4524. ! manner, and 2) to keep the vertical profile of temperature
  4525. ! close to the initial (base-state) profile, which is useful
  4526. ! for certain idealized applications.
  4527. ! Reference: Rotunno and Emanuel, 1987, JAS, p. 546
  4528. !-----------------------------------------------------------------------
  4529. ! Begin declarations.
  4530. IMPLICIT NONE
  4531. INTEGER, INTENT( IN ) &
  4532. :: ids, ide, jds, jde, kds, kde, &
  4533. ims, ime, jms, jme, kms, kme, &
  4534. its, ite, jts, jte, kts, kte
  4535. REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( INOUT ) &
  4536. :: t_tendf
  4537. REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( IN ) &
  4538. :: t, t_init, ph, phb
  4539. REAL, DIMENSION( ims:ime, jms:jme ), INTENT( IN ) &
  4540. :: mut
  4541. REAL, DIMENSION( kms:kme ) , INTENT(IN ) &
  4542. :: t_base, z_base
  4543. ! Local variables.
  4544. INTEGER :: i, j, k, ktf, k2
  4545. REAL :: tau_r , rmax , rmin , inv_tau_r , inv_g , rterm
  4546. REAL, DIMENSION( kms:kme ) :: z00,t00
  4547. ! End declarations.
  4548. !-----------------------------------------------------------------------
  4549. ! set tau_r to 12 h, following RE87
  4550. tau_r = 12.0*3600.0
  4551. ! limit rterm to +/- 2 K/day
  4552. rmax = 2.0/86400.0
  4553. rmin = -rmax
  4554. ktf = MIN( kte, kde-1 )
  4555. inv_tau_r = 1.0/tau_r
  4556. inv_g = 1.0/g
  4557. !-----------------------------------------------------------------------
  4558. ! Adjust potential temperature to base state.
  4559. DO j = jts, MIN( jte, jde-1 )
  4560. DO i = its, MIN( ite, ide-1 )
  4561. ! Get height of model levels:
  4562. DO k = kts, ktf
  4563. z00(k) = 0.5 * ( phb(i,k,j) + phb(i,k+1,j) + &
  4564. ph(i,k,j) + ph(i,k+1,j) ) * inv_g
  4565. ENDDO
  4566. ! Get reference state:
  4567. DO k = kts, ktf
  4568. k2 = ktf
  4569. DO WHILE( z_base(k2) .gt. z00(k) .and. k2 .gt. 1 )
  4570. k2 = k2 - 1
  4571. ENDDO
  4572. if(k2+1.gt.ktf)then
  4573. t00(k) = t_base(k2) + ( t_base(k2) - t_base(k2-1) ) &
  4574. * ( z00(k) - z_base(k2) ) &
  4575. / ( z_base(k2) - z_base(k2-1) )
  4576. else
  4577. t00(k) = t_base(k2) + ( t_base(k2+1) - t_base(k2) ) &
  4578. * ( z00(k) - z_base(k2) ) &
  4579. / ( z_base(k2+1) - z_base(k2) )
  4580. endif
  4581. ENDDO
  4582. ! Apply the RE87 R term:
  4583. DO k = kts, ktf
  4584. rterm = -( t(i,k,j) - t00(k) )*inv_tau_r
  4585. ! limit rterm:
  4586. rterm = min( rterm , rmax )
  4587. rterm = max( rterm , rmin )
  4588. t_tendf(i,k,j) = t_tendf(i,k,j) + mut(i,j)*rterm
  4589. END DO
  4590. END DO
  4591. END DO
  4592. END SUBROUTINE theta_relaxation
  4593. !==============================================================================
  4594. !==============================================================================
  4595. SUBROUTINE sixth_order_diffusion( name, field, tendency, mu, dt, &
  4596. config_flags, &
  4597. diff_6th_opt, diff_6th_factor, &
  4598. ids, ide, jds, jde, kds, kde, &
  4599. ims, ime, jms, jme, kms, kme, &
  4600. its, ite, jts, jte, kts, kte )
  4601. ! History: 14 Nov 2006 Name of variable changed by Jason Knievel
  4602. ! 07 Jun 2006 Revised and generalized by Jason Knievel
  4603. ! 25 Apr 2005 Original code by Jason Knievel, NCAR
  4604. ! Purpose: Apply 6th-order, monotonic (flux-limited), numerical
  4605. ! diffusion to 3-d velocity and to scalars.
  4606. ! References: Ming Xue (MWR Aug 2000)
  4607. ! Durran ("Numerical Methods for Wave Equations..." 1999)
  4608. ! George Bryan (personal communication)
  4609. !------------------------------------------------------------------------------
  4610. ! Begin: Declarations.
  4611. IMPLICIT NONE
  4612. INTEGER, INTENT(IN) &
  4613. :: ids, ide, jds, jde, kds, kde, &
  4614. ims, ime, jms, jme, kms, kme, &
  4615. its, ite, jts, jte, kts, kte
  4616. TYPE(grid_config_rec_type), INTENT(IN) &
  4617. :: config_flags
  4618. REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) &
  4619. :: tendency
  4620. REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN) &
  4621. :: field
  4622. REAL, DIMENSION( ims:ime , jms:jme ), INTENT(IN) &
  4623. :: mu
  4624. REAL, INTENT(IN) &
  4625. :: dt
  4626. REAL, INTENT(IN) &
  4627. :: diff_6th_factor
  4628. INTEGER, INTENT(IN) &
  4629. :: diff_6th_opt
  4630. CHARACTER(LEN=1) , INTENT(IN) &
  4631. :: name
  4632. INTEGER &
  4633. :: i, j, k, &
  4634. i_start, i_end, &
  4635. j_start, j_end, &
  4636. k_start, k_end, &
  4637. ktf
  4638. REAL &
  4639. :: dflux_x_p0, dflux_y_p0, &
  4640. dflux_x_p1, dflux_y_p1, &
  4641. tendency_x, tendency_y, &
  4642. mu_avg_p0, mu_avg_p1, &
  4643. diff_6th_coef
  4644. LOGICAL &
  4645. :: specified
  4646. ! End: Declarations.
  4647. !------------------------------------------------------------------------------
  4648. !------------------------------------------------------------------------------
  4649. ! Begin: Translate the diffusion factor into a diffusion coefficient. See
  4650. ! Durran's text, section 2.4.3, then adjust for sixth-order diffusion (not
  4651. ! fourth) and for diffusion in two dimensions (not one). For reference, a
  4652. ! factor of 1.0 would mean complete diffusion of a 2dx wave in one time step,
  4653. ! although application of the flux limiter reduces somewhat the effects of
  4654. ! diffusion for a given coefficient.
  4655. diff_6th_coef = diff_6th_factor * 0.015625 / ( 2.0 * dt )
  4656. ! End: Translate diffusion factor.
  4657. !------------------------------------------------------------------------------
  4658. !------------------------------------------------------------------------------
  4659. ! Begin: Assign limits of spatial loops depending on variable to be diffused.
  4660. ! The halo regions are already filled with values by the time this subroutine
  4661. ! is called, which allows the stencil to extend beyond the domains' edges.
  4662. ktf = MIN( kte, kde-1 )
  4663. IF ( name .EQ. 'u' ) THEN
  4664. i_start = its
  4665. i_end = ite
  4666. j_start = jts
  4667. j_end = MIN(jde-1,jte)
  4668. k_start = kts
  4669. k_end = ktf
  4670. ELSE IF ( name .EQ. 'v' ) THEN
  4671. i_start = its
  4672. i_end = MIN(ide-1,ite)
  4673. j_start = jts
  4674. j_end = jte
  4675. k_start = kts
  4676. k_end = ktf
  4677. ELSE IF ( name .EQ. 'w' ) THEN
  4678. i_start = its
  4679. i_end = MIN(ide-1,ite)
  4680. j_start = jts
  4681. j_end = MIN(jde-1,jte)
  4682. k_start = kts+1
  4683. k_end = ktf
  4684. ELSE
  4685. i_start = its
  4686. i_end = MIN(ide-1,ite)
  4687. j_start = jts
  4688. j_end = MIN(jde-1,jte)
  4689. k_start = kts
  4690. k_end = ktf
  4691. ENDIF
  4692. ! End: Assignment of limits of spatial loops.
  4693. !------------------------------------------------------------------------------
  4694. !------------------------------------------------------------------------------
  4695. ! Begin: Loop across spatial dimensions.
  4696. DO j = j_start, j_end
  4697. DO k = k_start, k_end
  4698. DO i = i_start, i_end
  4699. !------------------------------------------------------------------------------
  4700. ! Begin: Diffusion in x (i index).
  4701. ! Calculate the diffusive flux in x direction (from Xue's eq. 3).
  4702. dflux_x_p0 = ( 10.0 * ( field(i, k,j) - field(i-1,k,j) ) &
  4703. - 5.0 * ( field(i+1,k,j) - field(i-2,k,j) ) &
  4704. + ( field(i+2,k,j) - field(i-3,k,j) ) )
  4705. dflux_x_p1 = ( 10.0 * ( field(i+1,k,j) - field(i ,k,j) ) &
  4706. - 5.0 * ( field(i+2,k,j) - field(i-1,k,j) ) &
  4707. + ( field(i+3,k,j) - field(i-2,k,j) ) )
  4708. ! If requested in the namelist (diff_6th_opt=2), prohibit up-gradient diffusion
  4709. ! (variation on Xue's eq. 10).
  4710. IF ( diff_6th_opt .EQ. 2 ) THEN
  4711. IF ( dflux_x_p0 * ( field(i ,k,j)-field(i-1,k,j) ) .LE. 0.0 ) THEN
  4712. dflux_x_p0 = 0.0
  4713. END IF
  4714. IF ( dflux_x_p1 * ( field(i+1,k,j)-field(i ,k,j) ) .LE. 0.0 ) THEN
  4715. dflux_x_p1 = 0.0
  4716. END IF
  4717. END IF
  4718. ! Apply 6th-order diffusion in x direction.
  4719. IF ( name .EQ. 'u' ) THEN
  4720. mu_avg_p0 = mu(i-1,j)
  4721. mu_avg_p1 = mu(i ,j)
  4722. ELSE IF ( name .EQ. 'v' ) THEN
  4723. mu_avg_p0 = 0.25 * ( &
  4724. mu(i-1,j-1) + &
  4725. mu(i ,j-1) + &
  4726. mu(i-1,j ) + &
  4727. mu(i ,j ) )
  4728. mu_avg_p1 = 0.25 * ( &
  4729. mu(i ,j-1) + &
  4730. mu(i+1,j-1) + &
  4731. mu(i ,j ) + &
  4732. mu(i+1,j ) )
  4733. ELSE
  4734. mu_avg_p0 = 0.5 * ( &
  4735. mu(i-1,j) + &
  4736. mu(i ,j) )
  4737. mu_avg_p1 = 0.5 * ( &
  4738. mu(i ,j) + &
  4739. mu(i+1,j) )
  4740. END IF
  4741. tendency_x = diff_6th_coef * &
  4742. ( ( mu_avg_p1 * dflux_x_p1 ) - ( mu_avg_p0 * dflux_x_p0 ) )
  4743. ! End: Diffusion in x.
  4744. !------------------------------------------------------------------------------
  4745. !------------------------------------------------------------------------------
  4746. ! Begin: Diffusion in y (j index).
  4747. ! Calculate the diffusive flux in y direction (from Xue's eq. 3).
  4748. dflux_y_p0 = ( 10.0 * ( field(i,k,j ) - field(i,k,j-1) ) &
  4749. - 5.0 * ( field(i,k,j+1) - field(i,k,j-2) ) &
  4750. + ( field(i,k,j+2) - field(i,k,j-3) ) )
  4751. dflux_y_p1 = ( 10.0 * ( field(i,k,j+1) - field(i,k,j ) ) &
  4752. - 5.0 * ( field(i,k,j+2) - field(i,k,j-1) ) &
  4753. + ( field(i,k,j+3) - field(i,k,j-2) ) )
  4754. ! If requested in the namelist (diff_6th_opt=2), prohibit up-gradient diffusion
  4755. ! (variation on Xue's eq. 10).
  4756. IF ( diff_6th_opt .EQ. 2 ) THEN
  4757. IF ( dflux_y_p0 * ( field(i,k,j )-field(i,k,j-1) ) .LE. 0.0 ) THEN
  4758. dflux_y_p0 = 0.0
  4759. END IF
  4760. IF ( dflux_y_p1 * ( field(i,k,j+1)-field(i,k,j ) ) .LE. 0.0 ) THEN
  4761. dflux_y_p1 = 0.0
  4762. END IF
  4763. END IF
  4764. ! Apply 6th-order diffusion in y direction.
  4765. IF ( name .EQ. 'u' ) THEN
  4766. mu_avg_p0 = 0.25 * ( &
  4767. mu(i-1,j-1) + &
  4768. mu(i ,j-1) + &
  4769. mu(i-1,j ) + &
  4770. mu(i ,j ) )
  4771. mu_avg_p1 = 0.25 * ( &
  4772. mu(i-1,j ) + &
  4773. mu(i ,j ) + &
  4774. mu(i-1,j+1) + &
  4775. mu(i ,j+1) )
  4776. ELSE IF ( name .EQ. 'v' ) THEN
  4777. mu_avg_p0 = mu(i,j-1)
  4778. mu_avg_p1 = mu(i,j )
  4779. ELSE
  4780. mu_avg_p0 = 0.5 * ( &
  4781. mu(i,j-1) + &
  4782. mu(i,j ) )
  4783. mu_avg_p1 = 0.5 * ( &
  4784. mu(i,j ) + &
  4785. mu(i,j+1) )
  4786. END IF
  4787. tendency_y = diff_6th_coef * &
  4788. ( ( mu_avg_p1 * dflux_y_p1 ) - ( mu_avg_p0 * dflux_y_p0 ) )
  4789. ! End: Diffusion in y.
  4790. !------------------------------------------------------------------------------
  4791. !------------------------------------------------------------------------------
  4792. ! Begin: Combine diffusion in x and y.
  4793. tendency(i,k,j) = tendency(i,k,j) + tendency_x + tendency_y
  4794. ! End: Combine diffusion in x and y.
  4795. !------------------------------------------------------------------------------
  4796. ENDDO
  4797. ENDDO
  4798. ENDDO
  4799. ! End: Loop across spatial dimensions.
  4800. !------------------------------------------------------------------------------
  4801. END SUBROUTINE sixth_order_diffusion
  4802. !==============================================================================
  4803. !==============================================================================
  4804. END MODULE module_big_step_utilities_em