PageRenderTime 57ms CodeModel.GetById 21ms RepoModel.GetById 0ms app.codeStats 1ms

/wrfv2_fire/phys/module_bl_boulac.F

http://github.com/jbeezley/wrf-fire
FORTRAN Legacy | 921 lines | 568 code | 152 blank | 201 comment | 14 complexity | 0590a77da7a2c9d54b43c8e91ec7c0bb MD5 | raw file
Possible License(s): AGPL-1.0
  1. MODULE module_bl_boulac
  2. !USE module_model_constants
  3. !------------------------------------------------------------------------
  4. ! Calculation of the tendency due to momentum, heat
  5. ! and moisture turbulent fluxes follwing the approach
  6. ! of Bougeault and Lacarrere, 1989 (MWR, 117, 1872-1890).
  7. ! The scheme computes a prognostic ecuation for TKE and derives
  8. ! dissipation and turbulent coefficients using length scales.
  9. !
  10. ! Subroutine written by Alberto Martilli, CIEMAT, Spain,
  11. ! e-mail:alberto_martilli@ciemat.es
  12. ! August 2006.
  13. !------------------------------------------------------------------------
  14. ! IN THIS VERSION TKE IS NOT ADVECTED!!!!
  15. ! TO BE CHANGED IN THE FUTURE
  16. !
  17. ! -----------------------------------------------------------------------
  18. ! Constant used in the module
  19. ! ck_b=constant used in the compuation of diffusion coefficients
  20. ! ceps_b=constant used inthe computation of dissipation
  21. ! temin= minimum value allowed for TKE
  22. ! vk=von karman constant
  23. ! -----------------------------------------------------------------------
  24. real ck_b,ceps_b,vk,temin ! constant for Bougeault and Lacarrere
  25. parameter(ceps_b=1/1.4,ck_b=0.4,temin=0.0001,vk=0.4) ! impose minimum values for tke similar to those of MYJ
  26. ! -----------------------------------------------------------------------
  27. CONTAINS
  28. subroutine boulac(frc_urb2d,idiff,flag_bep,dz8w,dt,u_phy,v_phy &
  29. ,th_phy,rho,qv_curr,hfx &
  30. ,qfx,ustar,cp,g &
  31. ,rublten,rvblten,rthblten &
  32. ,rqvblten,rqcblten &
  33. ,tke,dlk,wu,wv,wt,wq,exch_h,exch_m,pblh &
  34. ,a_u_bep,a_v_bep,a_t_bep,a_q_bep &
  35. ,a_e_bep,b_u_bep,b_v_bep &
  36. ,b_t_bep,b_q_bep,b_e_bep,dlg_bep &
  37. ,dl_u_bep,sf_bep,vl_bep &
  38. ,ids,ide, jds,jde, kds,kde &
  39. ,ims,ime, jms,jme, kms,kme &
  40. ,its,ite, jts,jte, kts,kte)
  41. implicit none
  42. !-----------------------------------------------------------------------
  43. ! Input
  44. !------------------------------------------------------------------------
  45. INTEGER:: ids,ide, jds,jde, kds,kde, &
  46. ims,ime, jms,jme, kms,kme, &
  47. its,ite, jts,jte, kts,kte
  48. integer, INTENT(IN) :: idiff
  49. REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) :: DZ8W !vertical resolution
  50. REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) :: qv_curr !moisture
  51. REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) :: RHO !air density
  52. REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) :: TH_PHY !potential temperature
  53. REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) :: U_PHY !x-component of wind
  54. REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) :: V_PHY !y-component of wind
  55. REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN ) :: ustar !friction velocity
  56. REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN ) :: hfx !sensible heat flux (W/m2) at surface
  57. REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN ) :: qfx !moisture flux at surface
  58. real, INTENT(IN ) :: g,cp !gravity and Cp
  59. REAL, INTENT(IN ):: DT ! Time step
  60. REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN ) :: FRC_URB2D !fraction cover urban
  61. REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: PBLH !PBL height
  62. !
  63. ! variable added for urban
  64. REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) ::a_u_bep ! Implicit component for the momemtum in X-direction
  65. REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) ::a_v_bep ! Implicit component for the momemtum in Y-direction
  66. REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) ::a_t_bep ! Implicit component for the Pot. Temp.
  67. REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) ::a_q_bep ! Implicit component for Moisture
  68. REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) ::a_e_bep ! Implicit component for the TKE
  69. REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) ::b_u_bep ! Explicit component for the momemtum in X-direction
  70. REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) ::b_v_bep ! Explicit component for the momemtum in Y-direction
  71. REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) ::b_t_bep ! Explicit component for the Pot. Temp.
  72. REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) ::b_q_bep ! Explicit component for Moisture
  73. REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) ::b_e_bep ! Explicit component for the TKE
  74. REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::dlg_bep ! Height above ground (L_ground in formula (24) of the BLM paper).
  75. REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) ::dl_u_bep ! Length scale (lb in formula (22) ofthe BLM paper).
  76. ! urban surface and volumes
  77. REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) ::sf_bep ! surface of the urban grid cells
  78. REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) ::vl_bep ! volume of the urban grid cells
  79. LOGICAL, INTENT(IN) :: flag_bep !flag for BEP
  80. !
  81. !-----------------------------------------------------------------------
  82. ! Local, carried on from one timestep to the other
  83. !------------------------------------------------------------------------
  84. ! real, save, allocatable, dimension (:,:,:)::TKE ! Turbulent kinetic energy
  85. real, dimension (ims:ime, kms:kme, jms:jme) ::th_0 ! reference state for potential temperature
  86. !------------------------------------------------------------------------
  87. ! Output
  88. !------------------------------------------------------------------------
  89. real, dimension( ims:ime, kms:kme, jms:jme ), INTENT(OUT ) :: exch_h ! exchange coefficient for heat
  90. real, dimension( ims:ime, kms:kme, jms:jme ), INTENT(OUT ) :: exch_m ! exchange coefficient for momentum
  91. real, dimension( ims:ime, kms:kme, jms:jme ), INTENT(INOUT ) :: tke ! Turbulence Kinetic Energy
  92. real, dimension( ims:ime, kms:kme, jms:jme ), INTENT(OUT ) :: wu ! Turbulent flux of momentum (x)
  93. real, dimension( ims:ime, kms:kme, jms:jme ), INTENT(OUT ) :: wv ! Turbulent flux of momentum (y)
  94. real, dimension( ims:ime, kms:kme, jms:jme ), INTENT(OUT ) :: wt ! Turbulent flux of temperature
  95. real, dimension( ims:ime, kms:kme, jms:jme ), INTENT(OUT ) :: wq ! Turbulent flux of water vapor
  96. real, dimension( ims:ime, kms:kme, jms:jme ), INTENT(OUT ) :: dlk ! Turbulent flux of water vapor
  97. ! only if idiff not equal 1:
  98. REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(OUT ) :: RUBLTEN !tendency for U_phy
  99. REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(OUT ) :: RVBLTEN !tendency for V_phy
  100. REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(OUT ) :: RTHBLTEN !tendency for TH_phy
  101. REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(OUT ) :: RQVBLTEN !tendency for QV_curr
  102. REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(OUT ) :: RQCBLTEN !tendency for QV_curr
  103. !--------------------------------------------------------------
  104. ! Local
  105. !--------------------------------------------------------------
  106. ! 1D array used for the input and output of the routine boulac1D
  107. real z1D(kms:kme) ! vertical coordinates (faces of the grid)
  108. real dz1D(kms:kme) ! vertical resolution
  109. real u1D(kms:kme) ! wind speed in the x directions
  110. real v1D(kms:kme) ! wind speed in the y directions
  111. real th1D(kms:kme) ! potential temperature
  112. real q1D(kms:kme) ! potential temperature
  113. real rho1D(kms:kme) ! air density
  114. real rhoz1D(kms:kme) ! air density at the faces
  115. real tke1D(kms:kme) ! air pressure
  116. real th01D(kms:kme) ! reference potential temperature
  117. real dlk1D(kms:kme) ! dlk
  118. real dls1D(kms:kme) ! dls
  119. real exch1D(kms:kme) ! exch
  120. real sf1D(kms:kme) ! surface of the grid cells
  121. real vl1D(kms:kme) ! volume of the grid cells
  122. real a_u1D(kms:kme) ! Implicit component of the momentum sources or sinks in the X-direction
  123. real a_v1D(kms:kme) ! Implicit component of the momentum sources or sinks in the Y-direction
  124. real a_t1D(kms:kme) ! Implicit component of the heat sources or sinks
  125. real a_q1D(kms:kme) ! Implicit component of the moisture sources or sinks
  126. real a_e1D(kms:kme) ! Implicit component of the TKE sources or sinks
  127. real b_u1D(kms:kme) ! Explicit component of the momentum sources or sinks in the X-direction
  128. real b_v1D(kms:kme) ! Explicit component of the momentum sources or sinks in the Y-direction
  129. real b_t1D(kms:kme) ! Explicit component of the heat sources or sinks
  130. real b_q1D(kms:kme) ! Explicit component of the moisture sources or sinks
  131. real b_e1D(kms:kme) ! Explicit component of the TKE sources or sinks
  132. real dlg1D(kms:kme) ! Height above ground (L_ground in formula (24) of the BLM paper).
  133. real dl_u1D(kms:kme) ! Length scale (lb in formula (22) ofthe BLM paper)
  134. real sh1D(kms:kme) ! shear
  135. real bu1D(kms:kme) ! buoyancy
  136. real wu1D(kms:kme) ! turbulent flux of momentum (x component)
  137. real wv1D(kms:kme) ! turbulent flux of momentum (y component)
  138. real wt1D(kms:kme) ! turbulent flux of temperature
  139. real wq1D(kms:kme) ! turbulent flux of water vapor
  140. ! local added only for diagnostic output
  141. real a_e(ims:ime,kms:kme,jms:jme) ! implicit term in TKE
  142. real b_e(ims:ime,kms:kme,jms:jme) ! explicit term in TKE
  143. real bu(ims:ime,kms:kme,jms:jme) ! buoyancy term in TKE
  144. real sh(ims:ime,kms:kme,jms:jme) ! shear term in TKE
  145. real wrk(ims:ime) ! working array
  146. integer ix,iy,iz,id,iz_u,iw_u,ig,ir_u,ix1,iy1
  147. real ufrac_int ! urban fraction
  148. real vect,time_tke,hour,zzz
  149. real ustarf
  150. real summ1,summ2,summ3
  151. save time_tke,hour
  152. !
  153. ! store reference state for potential temperature
  154. ! and initialize tke
  155. !
  156. !here I fix the value of the reference state equal to the value of the potnetial temperature
  157. ! the only use of this variable in the code is to compute the paramter BETA = g/th0
  158. ! I fix it to 300K.
  159. do ix=its,ite
  160. do iy=jts,jte
  161. do iz=kts,kte
  162. ! th_0(ix,iz,iy)=th_phy(ix,iz,iy)
  163. th_0(ix,iz,iy)=300.
  164. enddo
  165. enddo
  166. enddo
  167. bu1d=0.
  168. sh1d=0.
  169. b_e1d=0.
  170. b_u1d=0.
  171. b_v1d=0.
  172. b_t1d=0.
  173. b_q1d=0.
  174. a_e1d=0.
  175. a_u1d=0.
  176. a_v1d=0.
  177. a_t1d=0.
  178. a_q1d=0.
  179. z1D=0. ! vertical coordinates (faces of the grid)
  180. dz1D=0. ! vertical resolution
  181. u1D =0. ! wind speed in the x directions
  182. v1D =0. ! wind speed in the y directions
  183. th1D=0. ! potential temperature
  184. q1D=0. ! potential temperature
  185. rho1D=0. ! air density
  186. rhoz1D=0. ! air density at the faces
  187. tke1D =0. ! air pressure
  188. th01D =0. ! reference potential temperature
  189. dlk1D =0. ! dlk
  190. dls1D =0. ! dls
  191. exch1D=0. ! exch
  192. sf1D =1. ! surface of the grid cells
  193. vl1D =1. ! volume of the grid cells
  194. a_u1D =0. ! Implicit component of the momentum sources or sinks in the X-direction
  195. a_v1D =0. ! Implicit component of the momentum sources or sinks in the Y-direction
  196. a_t1D =0. ! Implicit component of the heat sources or sinks
  197. a_q1D =0. ! Implicit component of the moisture sources or sinks
  198. a_e1D =0. ! Implicit component of the TKE sources or sinks
  199. b_u1D =0. ! Explicit component of the momentum sources or sinks in the X-direction
  200. b_v1D =0. ! Explicit component of the momentum sources or sinks in the Y-direction
  201. b_t1D =0. ! Explicit component of the heat sources or sinks
  202. b_q1D =0. ! Explicit component of the moisture sources or sinks
  203. b_e1D =0. ! Explicit component of the TKE sources or sinks
  204. dlg1D =0. ! Height above ground (L_ground in formula (24) of the BLM paper).
  205. dl_u1D=0. ! Length scale (lb in formula (22) ofthe BLM paper)
  206. sh1D =0. ! shear
  207. bu1D =0. ! buoyancy
  208. wu1D =0. ! turbulent flux of momentum (x component)
  209. wv1D =0. ! turbulent flux of momentum (y component)
  210. wt1D =0. ! turbulent flux of temperature
  211. wq1D =0. ! turbulent flux of water vapor
  212. ! local added only for diagnostic output
  213. ! loop over the columns.
  214. ! put variables in 1D temporary arrays
  215. !
  216. do ix=its,ite
  217. do iy=jts,jte
  218. z1d(kts)=0.
  219. do iz= kts,kte
  220. u1D(iz)=u_phy(ix,iz,iy)
  221. v1D(iz)=v_phy(ix,iz,iy)
  222. th1D(iz)=th_phy(ix,iz,iy)
  223. q1D(iz)=qv_curr(ix,iz,iy)
  224. tke1D(iz)=tke(ix,iz,iy)
  225. rho1D(iz)=rho(ix,iz,iy)
  226. th01D(iz)=th_0(ix,iz,iy)
  227. dz1D(iz)=dz8w(ix,iz,iy)
  228. z1D(iz+1)=z1D(iz)+dz1D(iz)
  229. enddo
  230. rhoz1D(kts)=rho1D(kts)
  231. do iz=kts+1,kte
  232. rhoz1D(iz)=(rho1D(iz)*dz1D(iz-1)+rho1D(iz-1)*dz1D(iz))/(dz1D(iz-1)+dz1D(iz))
  233. enddo
  234. rhoz1D(kte+1)=rho1D(kte)
  235. if(flag_bep)then
  236. do iz=kts,kte
  237. a_e1D(iz)=a_e_bep(ix,iz,iy)
  238. b_e1D(iz)=b_e_bep(ix,iz,iy)
  239. dlg1D(iz)=(z1D(iz)+z1D(iz+1))/2.*(1.-frc_urb2d(ix,iy))+dlg_bep(ix,iz,iy)*frc_urb2d(ix,iy)
  240. dl_u1D(iz)=dl_u_bep(ix,iz,iy)
  241. if((1.-frc_urb2d(ix,iy)).lt.1.)dl_u1D(iz)=dl_u1D(iz)/frc_urb2d(ix,iy)
  242. vl1D(iz)=vl_bep(ix,iz,iy)
  243. sf1D(iz)=sf_bep(ix,iz,iy)
  244. enddo
  245. ufrac_int=frc_urb2d(ix,iy)
  246. sf1D(kte+1)=sf_bep(ix,1,iy)
  247. else
  248. do iz=kts,kte
  249. a_e1D(iz)=0.
  250. b_e1D(iz)=0.
  251. dlg1D(iz)=(z1D(iz)+z1D(iz+1))/2.
  252. dl_u1D(iz)=0.
  253. vl1D(iz)=1.
  254. sf1D(iz)=1.
  255. enddo
  256. ufrac_int=0.
  257. sf1D(kte+1)=1.
  258. endif
  259. ! call the routine that will solve the turbulence in 1D for tke
  260. call boulac1D(ix,iy,ufrac_int,kms,kme,kts,kte,dz1d,z1D,dt,u1D,v1D,th1D,rho1D,rhoz1D,q1D,th01D,&
  261. tke1D,ustar(ix,iy),hfx(ix,iy),qfx(ix,iy),cp,g, &
  262. a_e1D,b_e1D, &
  263. dlg1D,dl_u1D,sf1D,vl1D,dlk1D,dls1D,exch1D,sh1D,bu1D)
  264. call pbl_height(kms,kme,kts,kte,dz1d,z1d,th1D,q1D,pblh(ix,iy))
  265. ! store turbulent exchange coefficients, TKE, and other variables
  266. do iz= kts,kte
  267. a_e(ix,iz,iy)=a_e1D(iz)
  268. b_e(ix,iz,iy)=b_e1D(iz)
  269. if(flag_bep)then
  270. dlg_bep(ix,iz,iy)=dlg1D(iz)
  271. endif
  272. tke(ix,iz,iy)=tke1D(iz)
  273. dlk(ix,iz,iy)=dlk1D(iz)
  274. sh(ix,iz,iy)=sh1D(iz)
  275. bu(ix,iz,iy)=bu1D(iz)
  276. exch_h(ix,iz,iy)=exch1D(iz)
  277. exch_m(ix,iz,iy)=exch1D(iz)
  278. enddo
  279. if(idiff.ne.1)then
  280. ! estimate the tendencies
  281. if(flag_bep)then
  282. do iz=kts,kte
  283. a_t1D(iz)=a_t_bep(ix,iz,iy)
  284. b_t1D(iz)=b_t_bep(ix,iz,iy)
  285. a_u1D(iz)=a_u_bep(ix,iz,iy)
  286. b_u1D(iz)=b_u_bep(ix,iz,iy)
  287. a_v1D(iz)=a_v_bep(ix,iz,iy)
  288. b_v1D(iz)=b_v_bep(ix,iz,iy)
  289. a_q1D(iz)=a_q_bep(ix,iz,iy)
  290. b_q1D(iz)=b_q_bep(ix,iz,iy)
  291. enddo
  292. else
  293. do iz=kts,kte
  294. a_t1D(iz)=0.
  295. b_t1D(iz)=0.
  296. a_u1D(iz)=0.
  297. b_u1D(iz)=0.
  298. a_v1D(iz)=0.
  299. b_v1D(iz)=0.
  300. a_q1D(iz)=0.
  301. b_q1D(iz)=0.
  302. enddo
  303. b_t1D(1)=hfx(ix,iy)/dz1D(1)/rho1D(1)/cp
  304. b_q1D(1)=qfx(ix,iy)/dz1D(1)/rho1D(1)
  305. a_u1D(1)=(-ustar(ix,iy)*ustar(ix,iy)/dz1D(1)/((u1D(1)**2.+v1D(1)**2.)**.5))
  306. a_v1D(1)=(-ustar(ix,iy)*ustar(ix,iy)/dz1D(1)/((u1D(1)**2.+v1D(1)**2.)**.5))
  307. endif
  308. ! solve diffusion equation for momentum x component
  309. call diff(kms,kme,kts,kte,1,1,dt,u1D,rho1D,rhoz1D,exch1D,a_u1D,b_u1D,sf1D,vl1D,dz1D,wu1D)
  310. ! solve diffusion equation for momentum y component
  311. call diff(kms,kme,kts,kte,1,1,dt,v1D,rho1D,rhoz1D,exch1D,a_v1D,b_v1D,sf1D,vl1D,dz1D,wv1D)
  312. ! solve diffusion equation for potential temperature
  313. call diff(kms,kme,kts,kte,1,1,dt,th1D,rho1D,rhoz1D,exch1D,a_t1D,b_t1D,sf1D,vl1D,dz1D,wt1D)
  314. ! solve diffusion equation for water vapor mixing ratio
  315. call diff(kms,kme,kts,kte,1,1,dt,q1D,rho1D,rhoz1D,exch1D,a_q1D,b_q1D,sf1D,vl1D,dz1D,wq1D)
  316. ! compute the tendencies
  317. do iz= kts,kte
  318. rthblten(ix,iz,iy)=rthblten(ix,iz,iy)+(th1D(iz)-th_phy(ix,iz,iy))/dt
  319. rqvblten(ix,iz,iy)=rqvblten(ix,iz,iy)+(q1D(iz)-qv_curr(ix,iz,iy))/dt
  320. rublten(ix,iz,iy)=rublten(ix,iz,iy)+(u1D(iz)-u_phy(ix,iz,iy))/dt
  321. rvblten(ix,iz,iy)=rvblten(ix,iz,iy)+(v1D(iz)-v_phy(ix,iz,iy))/dt
  322. wu(ix,iz,iy)=wu1D(iz)
  323. wv(ix,iz,iy)=wv1D(iz)
  324. wt(ix,iz,iy)=wt1D(iz)
  325. wq(ix,iz,iy)=wq1D(iz)
  326. enddo
  327. endif
  328. enddo ! iy
  329. enddo ! ix
  330. time_tke=time_tke+dt
  331. return
  332. end subroutine boulac
  333. ! ===6=8===============================================================72
  334. ! ===6=8===============================================================72
  335. subroutine boulac1D(ix,iy,ufrac_int,kms,kme,kts,kte,dz,z,dt,u,v,th,rho,rhoz,qa,th0,te, &
  336. ustar,hfx,qfx,cp,g, &
  337. a_e,b_e, &
  338. dlg,dl_u,sf,vl,dlk,dls,exch,sh,bu)
  339. ! ----------------------------------------------------------------------
  340. ! 1D resolution of TKE following Bougeault and Lacarrere
  341. ! ----------------------------------------------------------------------
  342. implicit none
  343. integer iz,ix,iy
  344. ! ----------------------------------------------------------------------
  345. ! INPUT:
  346. ! ----------------------------------------------------------------------
  347. integer kms,kme,kts,kte
  348. real z(kms:kme) ! Altitude above the ground of the cell interfaces.
  349. real dz(kms:kme) ! vertical resolution
  350. real u(kms:kme) ! Wind speed in the x direction
  351. real v(kms:kme) ! Wind speed in the y direction
  352. real th(kms:kme) ! Potential temperature
  353. real rho(kms:kme) ! Air density
  354. real g ! gravity
  355. real cp !
  356. real te(kms:kme) ! turbulent kinetic energy
  357. real qa(kms:kme) ! air humidity
  358. real th0(kms:kme) ! Reference potential temperature
  359. real dt ! Time step
  360. real ustar ! ustar
  361. real hfx ! sensbile heat flux
  362. real qfx ! kinematic latent heat flux
  363. real sf(kms:kme) ! surface of the urban grid cells
  364. real vl(kms:kme) ! volume of the urban grid cells
  365. real a_e(kms:kme) ! Implicit component of the TKE sources or sinks
  366. real b_e(kms:kme) ! Explicit component of the TKE sources or sinks
  367. real dlg(kms:kme) ! Height above ground (L_ground in formula (24) of the BLM paper).
  368. real dl_u(kms:kme) ! Length scale (lb in formula (22) ofthe BLM paper)
  369. real ufrac_int ! urban fraction
  370. ! local variables not needed in principle, but that can be used to estimate the budget and turbulent fluxes
  371. real we(kms:kme),dwe(kms:kme)
  372. ! local variables
  373. real sh(kms:kme) ! shear term in TKE eqn.
  374. real bu(kms:kme) ! buoyancy term in TKE eqn.
  375. real td(kms:kme) ! dissipation term in TKE eqn.
  376. real exch(kms:kme) ! turbulent diffusion coefficients (defined at the faces)
  377. real dls(kms:kme) ! dissipation length scale
  378. real dlk(kms:kme) ! length scale used to estimate exch
  379. real dlu(kms:kme) ! l_up
  380. real dld(kms:kme) ! l_down
  381. real rhoz(kms:kme) !air density at the faces of the cell
  382. real tstar ! derived from hfx and ustar
  383. real beta
  384. real summ1,summ2,summ3,summ4
  385. ! interpolate air density at the faces
  386. ! estimation of tstar
  387. tstar=-hfx/rho(1)/cp/ustar
  388. ! first compute values of dlu and dld (length scales up and down).
  389. call dissip_bougeault(ix,iy,g,kms,kme,kts,kte,z,dz,te,dlu,dld,th,th0)
  390. !then average them to obtain dls and dlk (length scales for dissipation and eddy coefficients)
  391. call length_bougeault(ix,iy,kms,kme,kts,kte,dld,dlu,dlg,dl_u,dls,dlk)
  392. ! compute the turbulent diffusion coefficients exch
  393. call cdtur_bougeault(ix,iy,kms,kme,kts,kte,te,z,dz,exch,dlk)
  394. ! compute source and sink terms in the TKE equation (shear, buoyancy and dissipation)
  395. call tke_bougeault(ix,iy,g,kms,kme,kts,kte,z,dz,vl,u,v,th,te,th0,ustar,tstar,exch,dls,td,sh,bu,b_e,a_e,sf,ufrac_int)
  396. ! solve for tke
  397. call diff(kms,kme,kts,kte,1,1,dt,te,rho,rhoz,exch,a_e,b_e,sf,vl,dz,we)
  398. ! avoid negative values for tke
  399. do iz=kts,kte
  400. if(te(iz).lt.temin) te(iz)=temin
  401. enddo
  402. return
  403. end subroutine boulac1d
  404. !
  405. ! ===6=8===============================================================72
  406. ! ===6=8===============================================================72
  407. subroutine dissip_bougeault(ix,iy,g,kms,kme,kts,kte,z,dz,te,dlu,dld,th,th0)
  408. ! compute the length scales up and down
  409. implicit none
  410. integer kms,kme,kts,kte,iz,izz,ix,iy
  411. real dzt,zup,beta,zup_inf,bbb,tl,zdo,zdo_sup,zzz,g
  412. real te(kms:kme),dlu(kms:kme),dld(kms:kme),dz(kms:kme)
  413. real z(kms:kme),th(kms:kme),th0(kms:kme)
  414. do iz=kts,kte
  415. zup=0.
  416. dlu(iz)=z(kte+1)-z(iz)-dz(iz)/2.
  417. zzz=0.
  418. zup_inf=0.
  419. beta=g/th0(iz) !Buoyancy coefficient
  420. do izz=iz,kte-1
  421. dzt=(dz(izz+1)+dz(izz))/2.
  422. zup=zup-beta*th(iz)*dzt
  423. zup=zup+beta*(th(izz+1)+th(izz))*dzt/2.
  424. zzz=zzz+dzt
  425. if(te(iz).lt.zup.and.te(iz).ge.zup_inf)then
  426. bbb=(th(izz+1)-th(izz))/dzt
  427. if(bbb.ne.0)then
  428. tl=(-beta*(th(izz)-th(iz))+sqrt( max(0.,(beta*(th(izz)-th(iz)))**2.+2.*bbb*beta*(te(iz)-zup_inf))))/bbb/beta
  429. else
  430. if(th(izz).ne.th(iz))then
  431. tl=(te(iz)-zup_inf)/(beta*(th(izz)-th(iz)))
  432. else
  433. tl=0.
  434. endif
  435. endif
  436. dlu(iz)=zzz-dzt+tl
  437. endif
  438. zup_inf=zup
  439. enddo
  440. zdo=0.
  441. zdo_sup=0.
  442. dld(iz)=z(iz)+dz(iz)/2.
  443. zzz=0.
  444. do izz=iz,kts+1,-1
  445. dzt=(dz(izz-1)+dz(izz))/2.
  446. zdo=zdo+beta*th(iz)*dzt
  447. zdo=zdo-beta*(th(izz-1)+th(izz))*dzt/2.
  448. zzz=zzz+dzt
  449. if(te(iz).lt.zdo.and.te(iz).ge.zdo_sup)then
  450. bbb=(th(izz)-th(izz-1))/dzt
  451. if(bbb.ne.0.)then
  452. tl=(beta*(th(izz)-th(iz))+sqrt( max(0.,(beta*(th(izz)-th(iz)))**2.+2.*bbb*beta*(te(iz)-zdo_sup))))/bbb/beta
  453. else
  454. if(th(izz).ne.th(iz))then
  455. tl=(te(iz)-zdo_sup)/(beta*(th(izz)-th(iz)))
  456. else
  457. tl=0.
  458. endif
  459. endif
  460. dld(iz)=zzz-dzt+tl
  461. endif
  462. zdo_sup=zdo
  463. enddo
  464. enddo
  465. end subroutine dissip_bougeault
  466. !
  467. ! ===6=8===============================================================72
  468. ! ===6=8===============================================================72
  469. subroutine length_bougeault(ix,iy,kms,kme,kts,kte,dld,dlu,dlg,dl_u,dls,dlk)
  470. ! compute the length scales for dissipation and turbulent coefficients
  471. implicit none
  472. integer kms,kme,kts,kte,iz,ix,iy
  473. real dlu(kms:kme),dld(kms:kme),dl_u(kms:kme)
  474. real dls(kms:kme),dlk(kms:kme),dlg(kms:kme)
  475. do iz=kts,kte
  476. dld(iz)=min(dld(iz),dlg(iz))
  477. dls(iz)=sqrt(dlu(iz)*dld(iz))
  478. dlk(iz)=min(dlu(iz),dld(iz))
  479. if(dl_u(iz).gt.0.)then
  480. dls(iz)=1./(1./dls(iz)+1./dl_u(iz))
  481. dlk(iz)=1./(1./dlk(iz)+1./dl_u(iz))
  482. endif
  483. enddo
  484. return
  485. end subroutine length_bougeault
  486. !
  487. ! ===6=8===============================================================72
  488. ! ===6=8===============================================================72
  489. subroutine cdtur_bougeault(ix,iy,kms,kme,kts,kte,te,z,dz,exch,dlk)
  490. ! compute turbulent coefficients
  491. implicit none
  492. integer iz,kms,kme,kts,kte,ix,iy
  493. real te_m,dlk_m
  494. real te(kms:kme),exch(kms:kme)
  495. real dz(kms:kme),z(kms:kme)
  496. real dlk(kms:kme)
  497. real fact
  498. exch(kts)=0.
  499. ! do iz=2,nz-1
  500. do iz=kts+1,kte
  501. te_m=(te(iz-1)*dz(iz)+te(iz)*dz(iz-1))/(dz(iz)+dz(iz-1))
  502. dlk_m=(dlk(iz-1)*dz(iz)+dlk(iz)*dz(iz-1))/(dz(iz)+dz(iz-1))
  503. exch(iz)=ck_b*dlk_m*sqrt(te_m)
  504. ! exch(iz)=max(exch(iz),0.0001)
  505. exch(iz)=max(exch(iz),0.1)
  506. enddo
  507. exch(kte+1)=0.1
  508. return
  509. end subroutine cdtur_bougeault
  510. ! ===6=8===============================================================72
  511. ! ===6=8===============================================================72
  512. subroutine diff(kms,kme,kts,kte,iz1,izf,dt,co,rho,rhoz,cd,aa,bb,sf,vl,dz,fc)
  513. !------------------------------------------------------------------------
  514. ! Calculation of the diffusion in 1D
  515. !------------------------------------------------------------------------
  516. ! - Input:
  517. ! nz : number of points
  518. ! iz1 : first calculated point
  519. ! co : concentration of the variable of interest
  520. ! dz : vertical levels
  521. ! cd : diffusion coefficients
  522. ! dtext : external time step
  523. ! rho : density of the air at the center
  524. ! rhoz : density of the air at the face
  525. ! itest : if itest eq 1 then update co, else store in a flux array
  526. ! - Output:
  527. ! co :concentration of the variable of interest
  528. ! - Internal:
  529. ! cddz : constant terms in the equations
  530. ! dt : diffusion time step
  531. ! nt : number of the diffusion time steps
  532. ! cstab : ratio of the stability condition for the time step
  533. !---------------------------------------------------------------------
  534. implicit none
  535. integer iz,iz1,izf
  536. integer kms,kme,kts,kte
  537. real dt,dzv
  538. real co(kms:kme),cd(kms:kme),dz(kms:kme)
  539. real rho(kms:kme),rhoz(kms:kme)
  540. real cddz(kms:kme+1),fc(kms:kme),df(kms:kme)
  541. real a(kms:kme,3),c(kms:kme)
  542. real sf(kms:kme),vl(kms:kme)
  543. real aa(kms:kme),bb(kms:kme)
  544. ! Compute cddz=2*cd/dz
  545. cddz(kts)=sf(kts)*rhoz(kts)*cd(kts)/dz(kts)
  546. do iz=kts+1,kte
  547. cddz(iz)=2.*sf(iz)*rhoz(iz)*cd(iz)/(dz(iz)+dz(iz-1))
  548. enddo
  549. cddz(kte+1)=sf(kte+1)*rhoz(kte+1)*cd(kte+1)/dz(kte)
  550. do iz=kts,iz1-1
  551. a(iz,1)=0.
  552. a(iz,2)=1.
  553. a(iz,3)=0.
  554. c(iz)=co(iz)
  555. enddo
  556. do iz=iz1,kte-izf
  557. dzv=vl(iz)*dz(iz)
  558. a(iz,1)=-cddz(iz)*dt/dzv/rho(iz)
  559. a(iz,2)=1+dt*(cddz(iz)+cddz(iz+1))/dzv/rho(iz)-aa(iz)*dt
  560. a(iz,3)=-cddz(iz+1)*dt/dzv/rho(iz)
  561. c(iz)=co(iz)+bb(iz)*dt
  562. enddo
  563. do iz=kte-(izf-1),kte
  564. a(iz,1)=0.
  565. a(iz,2)=1
  566. a(iz,3)=0.
  567. c(iz)=co(iz)
  568. enddo
  569. call invert (kms,kme,kts,kte,a,c,co)
  570. do iz=kts,iz1
  571. fc(iz)=0.
  572. enddo
  573. do iz=iz1+1,kte
  574. fc(iz)=-(cddz(iz)*(co(iz)-co(iz-1)))/rho(iz)
  575. enddo
  576. ! do iz=1,iz1
  577. ! df(iz)=0.
  578. ! enddo
  579. !
  580. ! do iz=iz1+1,nz-izf
  581. ! dzv=vl(iz)*dz(iz)
  582. ! df(iz)=+(co(iz-1)*cddz(iz)-co(iz)*(cddz(iz)+cddz(iz+1))+co(iz+1)*cddz(iz+1))/dzv/rho(iz)
  583. ! enddo
  584. !
  585. ! do iz=nz-izf,nz
  586. ! df(iz)=0.
  587. ! enddo
  588. return
  589. end subroutine diff
  590. ! ===6=8===============================================================72
  591. ! ===6=8===============================================================72
  592. subroutine buoy(ix,iy,g,kms,kme,kts,kte,th,th0,exch,dz,bu,ustar,tstar,ufrac_int)
  593. ! compute buoyancy term
  594. implicit none
  595. integer kms,kme,kts,kte,iz,ix,iy
  596. real dtdz1,dtdz2,cdm,dtmdz,g
  597. real th(kms:kme),exch(kms:kme),dz(kms:kme),bu(kms:kme)
  598. real th0(kms:kme),ustar,tstar,ufrac_int
  599. ! bu(1)=-ustar*tstar*g/th0(1)*(1.-ufrac_int)
  600. bu(kts)=0.
  601. do iz=kts+1,kte-1
  602. dtdz1=2.*(th(iz)-th(iz-1))/(dz(iz-1)+dz(iz))
  603. dtdz2=2.*(th(iz+1)-th(iz))/(dz(iz+1)+dz(iz))
  604. dtmdz=0.5*(dtdz1+dtdz2)
  605. cdm=0.5*(exch(iz+1)+exch(iz))
  606. bu(iz)=-cdm*dtmdz*g/th0(iz)
  607. enddo
  608. !
  609. bu(kte)=0.
  610. return
  611. end subroutine buoy
  612. ! ===6=8===============================================================72
  613. ! ===6=8===============================================================72
  614. subroutine shear(ix,iy,g,kms,kme,kts,kte,u,v,cdua,dz,sh,ustar,tstar,th,ufrac_int)
  615. ! compute shear term
  616. implicit none
  617. integer kms,kme,kts,kte,iz,ix,iy
  618. real dudz1,dudz2,dvdz1,dvdz2,cdm,dumdz,ustar
  619. real tstar,th,al,phim,g
  620. real u(kms:kme),v(kms:kme),cdua(kms:kme),dz(kms:kme),sh(kms:kme)
  621. real u1,u2,v1,v2,ufrac_int
  622. ! al=vk*g*tstar/(th*(ustar**2.))
  623. ! if(al.ge.0.)phim=1.+4.7*dz(1)/2.*al
  624. ! if(al.lt.0.)phim=(1.-15*dz(1)/2.*al)**(-0.25)
  625. !
  626. ! sh(1)=(ustar**3.)/vk/(dz(1)/2.)*(1.-ufrac_int)
  627. sh(kts)=0.
  628. do iz=kts+1,kte-1
  629. u2=(dz(iz+1)*u(iz)+dz(iz)*u(iz+1))/(dz(iz)+dz(iz+1))
  630. u1=(dz(iz)*u(iz-1)+dz(iz-1)*u(iz))/(dz(iz-1)+dz(iz))
  631. v2=(dz(iz+1)*v(iz)+dz(iz)*v(iz+1))/(dz(iz)+dz(iz+1))
  632. v1=(dz(iz)*v(iz-1)+dz(iz-1)*v(iz))/(dz(iz-1)+dz(iz))
  633. cdm=0.5*(cdua(iz)+cdua(iz+1))
  634. dumdz=((u2-u1)/dz(iz))**2.+((v2-v1)/dz(iz))**2.
  635. sh(iz)=cdm*dumdz
  636. enddo
  637. !!!!!!!
  638. sh(kte)=0.
  639. return
  640. end subroutine shear
  641. ! ===6=8===============================================================72
  642. ! ===6=8===============================================================72
  643. subroutine invert(kms,kme,kts,kte,a,c,x)
  644. !ccccccccccccccccccccccccccccccc
  645. ! Aim: Inversion and resolution of a tridiagonal matrix
  646. ! A X = C
  647. ! Input:
  648. ! a(*,1) lower diagonal (Ai,i-1)
  649. ! a(*,2) principal diagonal (Ai,i)
  650. ! a(*,3) upper diagonal (Ai,i+1)
  651. ! c
  652. ! Output
  653. ! x results
  654. !ccccccccccccccccccccccccccccccc
  655. implicit none
  656. integer in
  657. integer kts,kte,kms,kme
  658. real a(kms:kme,3),c(kms:kme),x(kms:kme)
  659. do in=kte-1,kts,-1
  660. c(in)=c(in)-a(in,3)*c(in+1)/a(in+1,2)
  661. a(in,2)=a(in,2)-a(in,3)*a(in+1,1)/a(in+1,2)
  662. enddo
  663. do in=kts+1,kte
  664. c(in)=c(in)-a(in,1)*c(in-1)/a(in-1,2)
  665. enddo
  666. do in=kts,kte
  667. x(in)=c(in)/a(in,2)
  668. enddo
  669. return
  670. end subroutine invert
  671. ! ===6=8===============================================================72
  672. ! ===6=8===============================================================72
  673. subroutine tke_bougeault(ix,iy,g,kms,kme,kts,kte,z,dz,vl,u,v,th,te,th0,ustar,tstar,exch, &
  674. dls,td,sh,bu,b_e,a_e,sf,ufrac_int)
  675. ! in this routine the shear, buoyancy and part of the dissipation terms
  676. ! of the TKE equation are computed
  677. implicit none
  678. integer kms,kme,kts,kte,iz,ix,iy
  679. real g,ustar,tstar,ufrac_int
  680. real z(kms:kme),dz(kms:kme),u(kms:kme),v(kms:kme),th(kms:kme),th0(kms:kme),te(kms:kme)
  681. real exch(kms:kme),dls(kms:kme),td(kms:kme),sh(kms:kme),bu(kms:kme)
  682. real a_e(kms:kme),b_e(kms:kme)
  683. real vl(kms:kme),sf(kms:kme)
  684. real te1,dl1
  685. call shear(ix,iy,g,kms,kme,kts,kte,u,v,exch,dz,sh,ustar,tstar,th(kts),ufrac_int)
  686. call buoy(ix,iy,g,kms,kme,kts,kte,th,th0,exch,dz,bu,ustar,tstar,ufrac_int)
  687. do iz=kts,kte
  688. te1=max(te(iz),temin)
  689. dl1=max(dls(iz),0.1)
  690. td(iz)=-ceps_b*sqrt(te1)/dl1
  691. sh(iz)=sh(iz)*sf(iz)
  692. bu(iz)=bu(iz)*sf(iz)
  693. a_e(iz)=a_e(iz)+td(iz)
  694. b_e(iz)=b_e(iz)+sh(iz)+bu(iz)
  695. enddo
  696. return
  697. end subroutine tke_bougeault
  698. !######################################################################
  699. subroutine pbl_height(kms,kme,kts,kte,dz,z,th,q,pblh)
  700. ! this routine computes the PBL height
  701. ! with an approach similar to MYNN
  702. implicit none
  703. integer kms,kme,kts,kte,iz
  704. real z(kms:kme),dz(kms:kme),th(kms:kme),q(kms:kme)
  705. real pblh
  706. !Local
  707. real thv(kms:kme),zc(kms:kme)
  708. real thsfc
  709. ! compute the height of the center of the grid cells
  710. do iz=kts,kte
  711. zc(iz)=z(iz)+dz(iz)/2.
  712. enddo
  713. ! compute the virtual potential temperature
  714. do iz=kts,kte
  715. thv(iz)=th(iz)*(1.+0.61*q(iz))
  716. enddo
  717. ! now compute the PBL height
  718. pblh=0.
  719. thsfc=thv(kts)+0.5
  720. do iz=kts+1,kte
  721. if(pblh.eq.0.and.thv(iz).gt.thsfc)then
  722. pblh=zc(iz-1)+(thsfc-thv(iz-1))/(max(0.01,thv(iz)-thv(iz-1)))*(zc(iz)-zc(iz-1))
  723. ! pblh=z(iz-1)+(thsfc-thv(iz-1))/(max(0.01,thv(iz)-thv(iz-1)))*(z(iz)-z(iz-1))
  724. endif
  725. enddo
  726. return
  727. end subroutine pbl_height
  728. ! ===6=8===============================================================72
  729. ! ===6=8===============================================================72
  730. SUBROUTINE BOULACINIT(RUBLTEN,RVBLTEN,RTHBLTEN,RQVBLTEN, &
  731. & TKE_PBL,EXCH_H,RESTART,ALLOWED_TO_READ, &
  732. & IDS,IDE,JDS,JDE,KDS,KDE, &
  733. & IMS,IME,JMS,JME,KMS,KME, &
  734. & ITS,ITE,JTS,JTE,KTS,KTE )
  735. !-----------------------------------------------------------------------
  736. IMPLICIT NONE
  737. !-----------------------------------------------------------------------
  738. LOGICAL,INTENT(IN) :: ALLOWED_TO_READ,RESTART
  739. INTEGER,INTENT(IN) :: IDS,IDE,JDS,JDE,KDS,KDE, &
  740. & IMS,IME,JMS,JME,KMS,KME, &
  741. & ITS,ITE,JTS,JTE,KTS,KTE
  742. REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(OUT) :: EXCH_H, &
  743. & RUBLTEN, &
  744. & RVBLTEN, &
  745. & RTHBLTEN, &
  746. & RQVBLTEN, &
  747. & TKE_PBL
  748. INTEGER :: I,J,K,ITF,JTF,KTF
  749. !-----------------------------------------------------------------------
  750. !-----------------------------------------------------------------------
  751. JTF=MIN0(JTE,JDE-1)
  752. KTF=MIN0(KTE,KDE-1)
  753. ITF=MIN0(ITE,IDE-1)
  754. IF(.NOT.RESTART)THEN
  755. DO J=JTS,JTF
  756. DO K=KTS,KTF
  757. DO I=ITS,ITF
  758. TKE_PBL(I,K,J)=0.0001
  759. RUBLTEN(I,K,J)=0.
  760. RVBLTEN(I,K,J)=0.
  761. RTHBLTEN(I,K,J)=0.
  762. RQVBLTEN(I,K,J)=0.
  763. EXCH_H(I,K,J)=0.
  764. ENDDO
  765. ENDDO
  766. ENDDO
  767. ENDIF
  768. END SUBROUTINE BOULACINIT
  769. END MODULE module_bl_boulac