PageRenderTime 26ms CodeModel.GetById 21ms RepoModel.GetById 0ms app.codeStats 0ms

/wrfv2_fire/dyn_em/module_force_scm.F

http://github.com/jbeezley/wrf-fire
FORTRAN Legacy | 554 lines | 461 code | 56 blank | 37 comment | 48 complexity | eebcd172279dcf342214e7a7220aec8b MD5 | raw file
Possible License(s): AGPL-1.0
  1. MODULE module_force_scm
  2. ! AUTHOR: Josh Hacker (NCAR/RAL)
  3. ! Forces a single-column (3x3) version of WRF
  4. CONTAINS
  5. SUBROUTINE force_scm(itimestep, dt, scm_force, dx, num_force_layers &
  6. , scm_th_adv, scm_qv_adv &
  7. , scm_ql_adv &
  8. , scm_wind_adv, scm_vert_adv &
  9. , scm_th_t_tend, scm_qv_t_tend &
  10. , scm_soilT_force, scm_soilQ_force &
  11. , scm_force_th_largescale &
  12. , scm_force_qv_largescale &
  13. , scm_force_ql_largescale &
  14. , scm_force_wind_largescale &
  15. , u_base, v_base, z_base &
  16. , z_force, z_force_tend &
  17. , u_g, v_g &
  18. , u_g_tend, v_g_tend &
  19. , w_subs, w_subs_tend &
  20. , th_upstream_x, th_upstream_x_tend &
  21. , th_upstream_y, th_upstream_y_tend &
  22. , qv_upstream_x, qv_upstream_x_tend &
  23. , qv_upstream_y, qv_upstream_y_tend &
  24. , ql_upstream_x, ql_upstream_x_tend &
  25. , ql_upstream_y, ql_upstream_y_tend &
  26. , u_upstream_x, u_upstream_x_tend &
  27. , u_upstream_y, u_upstream_y_tend &
  28. , v_upstream_x, v_upstream_x_tend &
  29. , v_upstream_y, v_upstream_y_tend &
  30. , th_t_tend, qv_t_tend &
  31. , tau_x, tau_x_tend &
  32. , tau_y, tau_y_tend &
  33. ,th_largescale &
  34. ,th_largescale_tend &
  35. ,qv_largescale &
  36. ,qv_largescale_tend &
  37. ,ql_largescale &
  38. ,ql_largescale_tend &
  39. ,u_largescale &
  40. ,u_largescale_tend &
  41. ,v_largescale &
  42. ,v_largescale_tend &
  43. ,tau_largescale &
  44. ,tau_largescale_tend &
  45. , num_force_soil_layers, num_soil_layers &
  46. , soil_depth_force, zs &
  47. , tslb, smois &
  48. , t_soil_forcing_val, t_soil_forcing_tend &
  49. , q_soil_forcing_val, q_soil_forcing_tend &
  50. , tau_soil &
  51. , z, z_at_w, th, qv, ql, u, v &
  52. , thten, qvten, qlten, uten, vten &
  53. , ids, ide, jds, jde, kds, kde &
  54. , ims, ime, jms, jme, kms, kme &
  55. , ips, ipe, jps, jpe, kps, kpe &
  56. , kts, kte &
  57. )
  58. ! adds forcing to bl tendencies and also to base state/geostrophic winds.
  59. USE module_init_utilities, ONLY : interp_0
  60. IMPLICIT NONE
  61. INTEGER, INTENT(IN ) :: itimestep
  62. INTEGER, INTENT(IN ) :: num_force_layers, scm_force
  63. REAL, INTENT(IN ) :: dt,dx
  64. LOGICAL, INTENT(IN ) :: scm_th_adv, &
  65. scm_qv_adv, &
  66. scm_ql_adv, &
  67. scm_wind_adv, &
  68. scm_vert_adv, &
  69. scm_soilT_force, &
  70. scm_soilQ_force, &
  71. scm_force_th_largescale, &
  72. scm_force_qv_largescale, &
  73. scm_force_ql_largescale, &
  74. scm_force_wind_largescale,&
  75. scm_th_t_tend,&
  76. scm_qv_t_tend
  77. REAL, DIMENSION(ims:ime,kms:kme,jms:jme), INTENT(IN ) :: z, th, qv, ql
  78. REAL, DIMENSION(ims:ime,kms:kme,jms:jme), INTENT(IN ) :: u, v
  79. REAL, DIMENSION(ims:ime,kms:kme,jms:jme), INTENT(IN ) :: z_at_w
  80. REAL, DIMENSION(ims:ime,kms:kme,jms:jme), INTENT(INOUT) :: thten, qvten
  81. REAL, DIMENSION(ims:ime,kms:kme,jms:jme), INTENT(INOUT) :: qlten
  82. REAL, DIMENSION(ims:ime,kms:kme,jms:jme), INTENT(INOUT) :: uten, vten
  83. REAL, DIMENSION( kms:kme ), INTENT(INOUT) :: u_base, v_base
  84. REAL, DIMENSION( kms:kme ), INTENT(INOUT) :: z_base
  85. REAL, DIMENSION(num_force_layers), INTENT (INOUT) :: z_force
  86. REAL, DIMENSION(num_force_layers), INTENT (INOUT) :: u_g,v_g
  87. REAL, DIMENSION(num_force_layers), INTENT (IN) :: z_force_tend
  88. REAL, DIMENSION(num_force_layers), INTENT (IN) :: u_g_tend,v_g_tend
  89. REAL, DIMENSION(num_force_layers), INTENT (IN) :: w_subs_tend
  90. REAL, DIMENSION(num_force_layers), INTENT (IN) :: th_upstream_x_tend
  91. REAL, DIMENSION(num_force_layers), INTENT (IN) :: th_upstream_y_tend
  92. REAL, DIMENSION(num_force_layers), INTENT (IN) :: qv_upstream_x_tend
  93. REAL, DIMENSION(num_force_layers), INTENT (IN) :: qv_upstream_y_tend
  94. REAL, DIMENSION(num_force_layers), INTENT (IN) :: ql_upstream_x_tend
  95. REAL, DIMENSION(num_force_layers), INTENT (IN) :: ql_upstream_y_tend
  96. REAL, DIMENSION(num_force_layers), INTENT (IN) :: u_upstream_x_tend
  97. REAL, DIMENSION(num_force_layers), INTENT (IN) :: u_upstream_y_tend
  98. REAL, DIMENSION(num_force_layers), INTENT (IN) :: v_upstream_x_tend
  99. REAL, DIMENSION(num_force_layers), INTENT (IN) :: v_upstream_y_tend
  100. REAL, DIMENSION(num_force_layers), INTENT (IN) :: th_t_tend
  101. REAL, DIMENSION(num_force_layers), INTENT (IN) :: qv_t_tend
  102. REAL, DIMENSION(num_force_layers), INTENT (IN) :: tau_x_tend
  103. REAL, DIMENSION(num_force_layers), INTENT (IN) :: tau_y_tend
  104. REAL, DIMENSION(num_force_layers), INTENT (INOUT) :: th_upstream_x
  105. REAL, DIMENSION(num_force_layers), INTENT (INOUT) :: th_upstream_y
  106. REAL, DIMENSION(num_force_layers), INTENT (INOUT) :: u_upstream_x
  107. REAL, DIMENSION(num_force_layers), INTENT (INOUT) :: u_upstream_y
  108. REAL, DIMENSION(num_force_layers), INTENT (INOUT) :: v_upstream_x
  109. REAL, DIMENSION(num_force_layers), INTENT (INOUT) :: v_upstream_y
  110. REAL, DIMENSION(num_force_layers), INTENT (INOUT) :: qv_upstream_x
  111. REAL, DIMENSION(num_force_layers), INTENT (INOUT) :: qv_upstream_y
  112. REAL, DIMENSION(num_force_layers), INTENT (INOUT) :: ql_upstream_x
  113. REAL, DIMENSION(num_force_layers), INTENT (INOUT) :: ql_upstream_y
  114. REAL, DIMENSION(num_force_layers), INTENT (INOUT) :: w_subs
  115. REAL, DIMENSION(num_force_layers), INTENT (INOUT) :: tau_x
  116. REAL, DIMENSION(num_force_layers), INTENT (INOUT) :: tau_y
  117. ! WA 1/8/10 for large-scale forcing
  118. REAL, DIMENSION(num_force_layers), INTENT (INOUT) :: th_largescale
  119. REAL, DIMENSION(num_force_layers), INTENT (INOUT) :: th_largescale_tend
  120. REAL, DIMENSION(num_force_layers), INTENT (INOUT) :: u_largescale
  121. REAL, DIMENSION(num_force_layers), INTENT (INOUT) :: u_largescale_tend
  122. REAL, DIMENSION(num_force_layers), INTENT (INOUT) :: v_largescale
  123. REAL, DIMENSION(num_force_layers), INTENT (INOUT) :: v_largescale_tend
  124. REAL, DIMENSION(num_force_layers), INTENT (INOUT) :: qv_largescale
  125. REAL, DIMENSION(num_force_layers), INTENT (INOUT) :: qv_largescale_tend
  126. REAL, DIMENSION(num_force_layers), INTENT (INOUT) :: ql_largescale
  127. REAL, DIMENSION(num_force_layers), INTENT (INOUT) :: ql_largescale_tend
  128. REAL, DIMENSION(num_force_layers), INTENT (INOUT) :: tau_largescale
  129. REAL, DIMENSION(num_force_layers), INTENT (INOUT) :: tau_largescale_tend
  130. ! WA 1/3/10 For soil forcing
  131. INTEGER, INTENT(IN ) :: num_force_soil_layers, num_soil_layers
  132. REAL, DIMENSION(ims:ime,num_soil_layers,jms:jme),INTENT(INOUT) :: tslb, smois
  133. REAL, DIMENSION(num_force_soil_layers), INTENT (INOUT) :: t_soil_forcing_val
  134. REAL, DIMENSION(num_force_soil_layers), INTENT (INOUT) :: t_soil_forcing_tend
  135. REAL, DIMENSION(num_force_soil_layers), INTENT (INOUT) :: q_soil_forcing_val
  136. REAL, DIMENSION(num_force_soil_layers), INTENT (INOUT) :: q_soil_forcing_tend
  137. REAL, DIMENSION(num_force_soil_layers), INTENT (INOUT) :: tau_soil
  138. REAL, DIMENSION(num_force_soil_layers), INTENT (IN ) :: soil_depth_force
  139. REAL, DIMENSION(num_soil_layers), INTENT (IN ) :: zs
  140. INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
  141. ims,ime, jms,jme, kms,kme, &
  142. ips,ipe, jps,jpe, kps,kpe, &
  143. kts,kte
  144. ! Local
  145. INTEGER :: i,j,k
  146. LOGICAL :: debug = .false.
  147. REAL :: t_x, t_y, qv_x, qv_y, ql_x, ql_y
  148. REAL :: u_x, u_y, v_x, v_y
  149. REAL, DIMENSION(kms:kme) :: th_adv_tend, qv_adv_tend, ql_adv_tend
  150. REAL, DIMENSION(kms:kme) :: u_adv_tend, v_adv_tend
  151. REAL, DIMENSION(kms:kme) :: th_t_tend_interp, qv_t_tend_interp
  152. REAL, DIMENSION(kms:kme) :: dthdz, dudz, dvdz, dqvdz, dqldz
  153. REAL :: w
  154. REAL, DIMENSION(kms:kme) :: w_dthdz, w_dudz, w_dvdz, w_dqvdz, w_dqldz
  155. REAL, DIMENSION(kms:kme) :: adv_timescale_x, adv_timescale_y
  156. CHARACTER*256 :: message
  157. ! Large-scale forcing WA 1/8/10
  158. REAL :: t_ls, qv_ls, ql_ls
  159. REAL :: u_ls, v_ls
  160. REAL, DIMENSION(kms:kme) :: th_ls_tend, qv_ls_tend, ql_ls_tend
  161. REAL, DIMENSION(kms:kme) :: u_ls_tend, v_ls_tend
  162. REAL, DIMENSION(kms:kme) :: ls_timescale
  163. ! Soil forcing WA 1/3/10
  164. INTEGER :: ks
  165. REAL :: t_soil, q_soil
  166. REAL, DIMENSION(num_soil_layers) :: t_soil_tend, q_soil_tend
  167. REAL, DIMENSION(num_soil_layers) :: timescale_soil
  168. IF ( scm_force .EQ. 0 ) return
  169. ! NOTES
  170. ! z is kts:kte
  171. ! z_at_w is kms:kme
  172. ! this is a good place for checks on the configuration
  173. if ( z_force(1) > z(ids,1,jds) ) then
  174. CALL wrf_message("First forcing level must be lower than first WRF half-level")
  175. WRITE( message , * ) 'z forcing = ',z_force(1), 'z = ',z(ids,1,jds)
  176. ! print*,"z forcing = ",z_force(1), "z = ",z(ids,1,jds)
  177. CALL wrf_error_fatal( message )
  178. endif
  179. z_force = z_force + dt*z_force_tend
  180. u_g = u_g + dt*u_g_tend
  181. v_g = v_g + dt*v_g_tend
  182. tau_x = tau_x + dt*tau_x_tend
  183. tau_y = tau_y + dt*tau_y_tend
  184. tau_largescale = tau_largescale + dt*tau_largescale_tend
  185. if ( scm_th_adv .AND. th_upstream_x(1) > 0.) then
  186. th_upstream_x = th_upstream_x + dt*th_upstream_x_tend
  187. th_upstream_y = th_upstream_y + dt*th_upstream_y_tend
  188. endif
  189. if ( scm_qv_adv .AND. qv_upstream_x(1) > 0.) then
  190. qv_upstream_x = qv_upstream_x + dt*qv_upstream_x_tend
  191. qv_upstream_y = qv_upstream_y + dt*qv_upstream_y_tend
  192. endif
  193. if ( scm_ql_adv .AND. ql_upstream_x(1) > 0.) then
  194. ql_upstream_x = ql_upstream_x + dt*ql_upstream_x_tend
  195. ql_upstream_y = ql_upstream_y + dt*ql_upstream_y_tend
  196. endif
  197. if ( scm_wind_adv .AND. u_upstream_x(1) > -900.) then
  198. u_upstream_x = u_upstream_x + dt*u_upstream_x_tend
  199. u_upstream_y = u_upstream_y + dt*u_upstream_y_tend
  200. v_upstream_x = v_upstream_x + dt*v_upstream_x_tend
  201. v_upstream_y = v_upstream_y + dt*v_upstream_y_tend
  202. endif
  203. if ( scm_vert_adv ) then
  204. w_subs = w_subs + dt*w_subs_tend
  205. endif
  206. if ( scm_force_th_largescale .AND. th_largescale(1) > 0.) then
  207. th_largescale = th_largescale + dt*th_largescale_tend
  208. endif
  209. if ( scm_force_qv_largescale .AND. qv_largescale(1) > 0.) then
  210. qv_largescale = qv_largescale + dt*qv_largescale_tend
  211. endif
  212. if ( scm_force_ql_largescale.AND. ql_largescale(1) > 0.) then
  213. ql_largescale = ql_largescale + dt*ql_largescale_tend
  214. endif
  215. if ( scm_force_wind_largescale .AND. u_largescale(1) > -900.) then
  216. u_largescale = u_largescale + dt*u_largescale_tend
  217. v_largescale = v_largescale + dt*v_largescale_tend
  218. endif
  219. if ( scm_soilT_force ) then
  220. t_soil_forcing_val = t_soil_forcing_val + dt*t_soil_forcing_tend
  221. endif
  222. if ( scm_soilQ_force ) then
  223. q_soil_forcing_val = q_soil_forcing_val + dt*q_soil_forcing_tend
  224. endif
  225. ! 0 everything in case we don't set it later
  226. th_adv_tend = 0.0
  227. qv_adv_tend = 0.0
  228. ql_adv_tend = 0.0
  229. u_adv_tend = 0.0
  230. v_adv_tend = 0.0
  231. th_ls_tend = 0.0
  232. qv_ls_tend = 0.0
  233. ql_ls_tend = 0.0
  234. u_ls_tend = 0.0
  235. v_ls_tend = 0.0
  236. w_dthdz = 0.0
  237. w_dqvdz = 0.0
  238. w_dqldz = 0.0
  239. w_dudz = 0.0
  240. w_dvdz = 0.0
  241. adv_timescale_x = 0.0
  242. adv_timescale_y = 0.0
  243. th_t_tend_interp =0.0
  244. qv_t_tend_interp =0.0
  245. ! now interpolate forcing to model vertical grid
  246. ! if ( debug ) print*,' z u_base v_base '
  247. CALL wrf_debug(100,'k z_base u_base v_base')
  248. do k = kms,kme-1
  249. z_base(k) = z(ids,k,jds)
  250. u_base(k) = interp_0(u_g,z_force,z_base(k),num_force_layers)
  251. v_base(k) = interp_0(v_g,z_force,z_base(k),num_force_layers)
  252. ! if ( debug ) print*,z_base(k),u_base(k),v_base(k)
  253. WRITE( message, '(i4,3f12.4)' ) k,z_base(k),u_base(k),v_base(k)
  254. CALL wrf_debug ( 100, message )
  255. enddo
  256. if ( scm_th_adv .or. scm_qv_adv .or. scm_ql_adv .or. scm_wind_adv ) then
  257. if ( scm_th_adv ) CALL wrf_debug ( 100, 'k tau_x tau_y t_ups_x t_ups_y t_m ' )
  258. do k = kms,kme-1
  259. adv_timescale_x(k) = interp_0(tau_x,z_force,z(ids,k,jds),num_force_layers)
  260. adv_timescale_y(k) = interp_0(tau_y,z_force,z(ids,k,jds),num_force_layers)
  261. enddo
  262. endif
  263. if ( scm_th_adv ) then
  264. if ( th_upstream_x(1) > 0.) then
  265. do k = kms,kme-1
  266. t_x = interp_0(th_upstream_x,z_force,z(ids,k,jds),num_force_layers)
  267. t_y = interp_0(th_upstream_y,z_force,z(ids,k,jds),num_force_layers)
  268. th_adv_tend(k) = (t_x-th(ids,k,jds))/adv_timescale_x(k) + (t_y-th(ids,k,jds))/adv_timescale_y(k)
  269. WRITE( message, '(i4,5f12.4)' ) k,adv_timescale_x(k), adv_timescale_y(k), t_x, t_y, th(ids,k,jds)
  270. CALL wrf_debug ( 100, message )
  271. enddo
  272. else ! WA if upstream is empty, use tendency only not value+tend
  273. do k = kms,kme-1
  274. t_x = interp_0(dt*th_upstream_x_tend,z_force,z(ids,k,jds),num_force_layers)
  275. t_y = interp_0(dt*th_upstream_y_tend,z_force,z(ids,k,jds),num_force_layers)
  276. th_adv_tend(k) = t_x/adv_timescale_x(k) + t_y/adv_timescale_y(k)
  277. WRITE( message, '(i4,5f12.4)' ) k,adv_timescale_x(k), adv_timescale_y(k), t_x, t_y, th(ids,k,jds)
  278. CALL wrf_debug ( 100, message )
  279. enddo
  280. endif
  281. endif
  282. if (minval(tau_x) < 0) then
  283. print*,tau_x
  284. stop 'TAU_X'
  285. endif
  286. if (minval(tau_y) < 0) then
  287. print*,z_force
  288. print*,tau_y
  289. stop 'TAU_Y'
  290. endif
  291. if ( scm_qv_adv ) then
  292. if ( qv_upstream_x(1) > 0.) then
  293. do k = kms,kme-1
  294. qv_x = interp_0(qv_upstream_x,z_force,z(ids,k,jds),num_force_layers)
  295. qv_y = interp_0(qv_upstream_y,z_force,z(ids,k,jds),num_force_layers)
  296. qv_adv_tend(k) = (qv_x-qv(ids,k,jds))/adv_timescale_x(k) + (qv_y-qv(ids,k,jds))/adv_timescale_y(k)
  297. WRITE( message, * ) 'qv_adv_tend branch 1',k,adv_timescale_x(k), qv_upstream_x(k), adv_timescale_y(k), qv_x, qv_y, qv(ids,k,jds), qv_adv_tend(k)
  298. CALL wrf_debug ( 100, message )
  299. enddo
  300. else ! WA if upstream is empty, use tendency only not value+tend
  301. do k = kms,kme-1
  302. qv_x = interp_0(dt*qv_upstream_x_tend,z_force,z(ids,k,jds),num_force_layers)
  303. qv_y = interp_0(dt*qv_upstream_y_tend,z_force,z(ids,k,jds),num_force_layers)
  304. qv_adv_tend(k) = qv_x/adv_timescale_x(k) + qv_y/adv_timescale_y(k)
  305. WRITE( message, * ) 'qv_adv_tend branch 2',k,adv_timescale_x(k), adv_timescale_y(k), qv_upstream_x(k), qv_x, qv_y, qv(ids,k,jds), qv_adv_tend(k)
  306. CALL wrf_debug ( 100, message )
  307. enddo
  308. endif
  309. endif
  310. if ( scm_ql_adv ) then
  311. if ( ql_upstream_x(1) > 0.) then
  312. do k = kms,kme-1
  313. ql_x = interp_0(ql_upstream_x,z_force,z(ids,k,jds),num_force_layers)
  314. ql_y = interp_0(ql_upstream_y,z_force,z(ids,k,jds),num_force_layers)
  315. ql_adv_tend(k) = (ql_x-ql(ids,k,jds))/adv_timescale_x(k) + (ql_y-ql(ids,k,jds))/adv_timescale_y(k)
  316. enddo
  317. else ! WA if upstream is empty, use tendency only not value+tend
  318. do k = kms,kme-1
  319. ql_x = interp_0(dt*ql_upstream_x_tend,z_force,z(ids,k,jds),num_force_layers)
  320. ql_y = interp_0(dt*ql_upstream_y_tend,z_force,z(ids,k,jds),num_force_layers)
  321. ql_adv_tend(k) = ql_x/adv_timescale_x(k) + ql_y/adv_timescale_y(k)
  322. enddo
  323. endif
  324. endif
  325. if ( scm_wind_adv ) then
  326. if ( u_upstream_x(1) > -900.) then
  327. do k = kms,kme-1
  328. u_x = interp_0(u_upstream_x,z_force,z(ids,k,jds),num_force_layers)
  329. u_y = interp_0(u_upstream_y,z_force,z(ids,k,jds),num_force_layers)
  330. v_x = interp_0(v_upstream_x,z_force,z(ids,k,jds),num_force_layers)
  331. v_y = interp_0(v_upstream_y,z_force,z(ids,k,jds),num_force_layers)
  332. u_adv_tend(k) = (u_x-u(ids,k,jds))/adv_timescale_x(k) + (u_y-u(ids,k,jds))/adv_timescale_y(k)
  333. v_adv_tend(k) = (v_x-v(ids,k,jds))/adv_timescale_x(k) + (v_y-v(ids,k,jds))/adv_timescale_y(k)
  334. enddo
  335. else ! WA if upstream is empty, use tendency only not value+tend
  336. do k = kms,kme-1
  337. u_x = interp_0(dt*u_upstream_x_tend,z_force,z(ids,k,jds),num_force_layers)
  338. u_y = interp_0(dt*u_upstream_y_tend,z_force,z(ids,k,jds),num_force_layers)
  339. v_x = interp_0(dt*v_upstream_x_tend,z_force,z(ids,k,jds),num_force_layers)
  340. v_y = interp_0(dt*v_upstream_y_tend,z_force,z(ids,k,jds),num_force_layers)
  341. u_adv_tend(k) = u_x/adv_timescale_x(k) + u_y/adv_timescale_y(k)
  342. v_adv_tend(k) = v_x/adv_timescale_x(k) + v_y/adv_timescale_y(k)
  343. enddo
  344. endif
  345. endif
  346. if ( scm_th_t_tend ) then
  347. do k = kms,kme-1
  348. th_t_tend_interp(k) = interp_0(th_t_tend,z_force,z(ids,k,jds),num_force_layers)
  349. enddo
  350. endif
  351. if ( scm_qv_t_tend ) then
  352. do k = kms,kme-1
  353. qv_t_tend_interp(k) = interp_0(qv_t_tend,z_force,z(ids,k,jds),num_force_layers)
  354. write(*,'(i3, f20.15)') k, qv_t_tend_interp(k)
  355. enddo
  356. endif
  357. ! Large scale forcing starts here 1/8/10 WA
  358. if ( scm_force_th_largescale .or. scm_force_qv_largescale .or. scm_force_ql_largescale .or. scm_force_wind_largescale ) then
  359. do k = kms,kme-1
  360. ls_timescale(k) = interp_0(tau_largescale,z_force,z(ids,k,jds),num_force_layers)
  361. enddo
  362. endif
  363. if ( scm_force_th_largescale ) then
  364. if ( th_largescale(1) > 0.) then
  365. do k = kms,kme-1
  366. t_ls = interp_0(th_largescale,z_force,z(ids,k,jds),num_force_layers)
  367. th_ls_tend(k) = (t_ls-th(ids,k,jds))/ls_timescale(k)
  368. enddo
  369. else ! WA if upstream is empty, use tendency only not value+tend
  370. do k = kms,kme-1
  371. t_ls = interp_0(dt*th_largescale_tend,z_force,z(ids,k,jds),num_force_layers)
  372. th_ls_tend(k) = t_ls/ls_timescale(k)
  373. enddo
  374. endif
  375. endif
  376. if ( scm_force_qv_largescale ) then
  377. if ( qv_largescale(1) > 0.) then
  378. do k = kms,kme-1
  379. qv_ls = interp_0(qv_largescale,z_force,z(ids,k,jds),num_force_layers)
  380. qv_ls_tend(k) = (qv_ls-qv(ids,k,jds))/ls_timescale(k)
  381. enddo
  382. else ! WA if upstream is empty, use tendency only not value+tend
  383. do k = kms,kme-1
  384. qv_ls = interp_0(dt*qv_largescale_tend,z_force,z(ids,k,jds),num_force_layers)
  385. qv_ls_tend(k) = qv_ls/ls_timescale(k)
  386. enddo
  387. endif
  388. endif
  389. if ( scm_force_ql_largescale ) then
  390. if ( ql_largescale(1) > 0.) then
  391. do k = kms,kme-1
  392. ql_ls = interp_0(ql_largescale,z_force,z(ids,k,jds),num_force_layers)
  393. ql_ls_tend(k) = (ql_ls-ql(ids,k,jds))/ls_timescale(k)
  394. enddo
  395. else ! WA if upstream is empty, use tendency only not value+tend
  396. do k = kms,kme-1
  397. ql_ls = interp_0(dt*ql_largescale_tend,z_force,z(ids,k,jds),num_force_layers)
  398. ql_ls_tend(k) = ql_ls/ls_timescale(k)
  399. enddo
  400. endif
  401. endif
  402. if ( scm_force_wind_largescale ) then
  403. if ( u_largescale(1) > -900.) then
  404. do k = kms,kme-1
  405. u_ls = interp_0(u_largescale,z_force,z(ids,k,jds),num_force_layers)
  406. v_ls = interp_0(v_largescale,z_force,z(ids,k,jds),num_force_layers)
  407. u_ls_tend(k) = (u_ls-u(ids,k,jds))/ls_timescale(k)
  408. v_ls_tend(k) = (v_ls-v(ids,k,jds))/ls_timescale(k)
  409. enddo
  410. else ! WA if upstream is empty, use tendency only not value+tend
  411. do k = kms,kme-1
  412. u_ls = interp_0(dt*u_largescale_tend,z_force,z(ids,k,jds),num_force_layers)
  413. v_ls = interp_0(dt*v_largescale_tend,z_force,z(ids,k,jds),num_force_layers)
  414. u_ls_tend(k) = u_ls/ls_timescale(k)
  415. v_ls_tend(k) = v_ls/ls_timescale(k)
  416. enddo
  417. endif
  418. endif
  419. ! Now do vertical advection. Note that no large-scale vertical advection
  420. ! is implemented at this time, may not make sense anyway (WA).
  421. ! loops are set so that the top and bottom (w=0) are handled correctly
  422. ! vertical derivatives
  423. do k = kms+1,kme-1
  424. dthdz(k) = (th(2,k,2)-th(2,k-1,2))/(z(2,k,2)-z(2,k-1,2))
  425. dqvdz(k) = (qv(2,k,2)-qv(2,k-1,2))/(z(2,k,2)-z(2,k-1,2))
  426. dqldz(k) = (ql(2,k,2)-ql(2,k-1,2))/(z(2,k,2)-z(2,k-1,2))
  427. dudz(k) = (u(2,k,2)-u(2,k-1,2))/(z(2,k,2)-z(2,k-1,2))
  428. dvdz(k) = (v(2,k,2)-v(2,k-1,2))/(z(2,k,2)-z(2,k-1,2))
  429. enddo
  430. ! w on full levels, then advect
  431. if ( scm_vert_adv ) then
  432. do k = kms+1,kme-1
  433. w = interp_0(w_subs,z_force,z_at_w(ids,k,jds),num_force_layers)
  434. w_dthdz(k) = -w*dthdz(k)
  435. w_dqvdz(k) = -w*dqvdz(k)
  436. w_dqldz(k) = -w*dqldz(k)
  437. w_dudz(k) = -w*dudz(k)
  438. w_dvdz(k) = -w*dvdz(k)
  439. enddo
  440. endif
  441. ! set tendencies for return
  442. ! vertical advection tendencies need to be interpolated back to half levels
  443. CALL wrf_debug ( 100, 'j, k, th_adv_ten, qv_adv_ten, ql_adv_ten, u_adv_ten, v_adv_ten')
  444. do j = jms,jme
  445. do k = kms,kme-1
  446. if(j==1) WRITE( message, * ) k,th_adv_tend(k),qv_adv_tend(k),ql_adv_tend(k), u_adv_tend(k),v_adv_tend(k)
  447. if(j==1) CALL wrf_debug ( 100, message )
  448. do i = ims,ime
  449. thten(i,k,j) = thten(i,k,j) + th_adv_tend(k) + &
  450. 0.5*(w_dthdz(k) + w_dthdz(k+1)) + th_t_tend_interp(k)&
  451. + th_ls_tend(k)
  452. qvten(i,k,j) = qvten(i,k,j) + qv_adv_tend(k) + &
  453. 0.5*(w_dqvdz(k) + w_dqvdz(k+1)) + qv_t_tend_interp(k)&
  454. + qv_ls_tend(k)
  455. qlten(i,k,j) = qlten(i,k,j) + ql_adv_tend(k) + &
  456. 0.5*(w_dqldz(k) + w_dqldz(k+1)) &
  457. + ql_ls_tend(k)
  458. uten(i,k,j) = uten(i,k,j) + u_adv_tend(k) + &
  459. 0.5*(w_dudz(k) + w_dudz(k+1)) &
  460. + u_ls_tend(k)
  461. vten(i,k,j) = vten(i,k,j) + v_adv_tend(k) + &
  462. 0.5*(w_dvdz(k) + w_dvdz(k+1)) &
  463. + v_ls_tend(k)
  464. enddo
  465. enddo
  466. enddo
  467. ! soil forcing 1/3/10 WA
  468. if ( scm_soilT_force ) then
  469. do ks = 1,num_soil_layers
  470. t_soil = interp_0(t_soil_forcing_val,soil_depth_force,zs(ks),num_force_soil_layers)
  471. timescale_soil(ks) = interp_0(tau_soil,soil_depth_force,zs(ks),num_force_soil_layers)
  472. t_soil_tend(ks) = (t_soil-tslb(ids,ks,jds))/timescale_soil(ks)
  473. enddo
  474. do j = jms,jme
  475. do ks = 1,num_soil_layers
  476. do i = ims,ime
  477. tslb(ids,ks,jds) = tslb(ids,ks,jds) + t_soil_tend(ks)
  478. enddo
  479. enddo
  480. enddo
  481. endif
  482. if ( scm_soilQ_force ) then
  483. do ks = 1,num_soil_layers
  484. q_soil = interp_0(q_soil_forcing_val,soil_depth_force,zs(ks),num_force_soil_layers)
  485. timescale_soil(ks) = interp_0(tau_soil,soil_depth_force,zs(ks),num_force_soil_layers)
  486. q_soil_tend(ks) = (q_soil-smois(ids,ks,jds))/timescale_soil(ks)
  487. enddo
  488. do j = jms,jme
  489. do ks = 1,num_soil_layers
  490. do i = ims,ime
  491. smois(ids,ks,jds) = smois(ids,ks,jds) + q_soil_tend(ks)
  492. enddo
  493. enddo
  494. enddo
  495. endif
  496. RETURN
  497. END SUBROUTINE force_scm
  498. END MODULE module_force_scm