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

/WPS/util/src/calc_ecmwf_p.F

http://github.com/jbeezley/wrf-fire
FORTRAN Legacy | 446 lines | 237 code | 104 blank | 105 comment | 58 complexity | f10b4012f397ee3620d75f1454b4c7f6 MD5 | raw file
Possible License(s): AGPL-1.0
  1. module coefficients
  2. integer :: n_levels
  3. real, allocatable, dimension(:) :: a, b
  4. contains
  5. !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  6. ! Name: read_coeffs
  7. !
  8. ! Notes: Obtain table of coefficients for input by this routine from
  9. ! http://www.ecmwf.int/products/data/technical/model_levels/index.html
  10. !
  11. !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  12. subroutine read_coeffs()
  13. implicit none
  14. integer :: i, nlvl, istatus
  15. open(21,file='ecmwf_coeffs',form='formatted',status='old',iostat=istatus)
  16. n_levels = 0
  17. if (istatus /= 0) then
  18. write(6,*) 'ERROR: Error opening ecmwf_coeffs'
  19. return
  20. end if
  21. read(21,*,iostat=istatus) nlvl
  22. do while (istatus == 0)
  23. n_levels = n_levels + 1
  24. read(21,*,iostat=istatus) nlvl
  25. end do
  26. rewind(21)
  27. n_levels = n_levels - 1
  28. allocate(a(0:n_levels))
  29. allocate(b(0:n_levels))
  30. write(6,*) ' '
  31. write(6,*) 'Coefficients for each level:',n_levels
  32. do i=0,n_levels
  33. read(21,*,iostat=istatus) nlvl, a(i), b(i)
  34. write(6,'(i5,5x,f12.6,2x,f12.10)') nlvl, a(i), b(i)
  35. end do
  36. write(6,*) ' '
  37. close(21)
  38. end subroutine read_coeffs
  39. !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  40. ! Name: cleanup_coeffs
  41. !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  42. subroutine cleanup_coeffs()
  43. implicit none
  44. n_levels = 0
  45. if (allocated(a)) deallocate(a)
  46. if (allocated(b)) deallocate(b)
  47. end subroutine cleanup_coeffs
  48. end module coefficients
  49. !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  50. ! calc_ecmwf_p
  51. !
  52. ! The purpose of this program is to compute a 3d pressure field for ECMWF
  53. ! model-level data sets; the code works in the WPS intermediate file format,
  54. ! reading a PSFC field from intermediate files, the A and B coefficients
  55. ! from a text file, ecmwf_coeffs, and writes the pressure data to an
  56. ! intermediate file.
  57. !
  58. ! November 2008
  59. ! Note: This program now also computes height for each level, this is needed by real
  60. ! modified by: Daniel van Dijke, MeteoConsult B.V., The Netherlands
  61. ! Chiel van Heerwaarden, Wageningen University, The Netherlands
  62. !
  63. !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  64. program calc_ecmwf_p
  65. use coefficients
  66. use date_pack
  67. use gridinfo_module
  68. use misc_definitions_module
  69. use module_debug
  70. use read_met_module
  71. use stringutil
  72. use write_met_module
  73. implicit none
  74. ! Local variables
  75. integer :: i, idiff, n_times, t, istatus, fg_idx, counter
  76. real :: a_full, b_full
  77. character (len=19) :: valid_date, temp_date
  78. character (len=128) :: input_name
  79. real, allocatable, dimension(:,:) :: psfc, hgtsfc, hgtprev, pstart, pend
  80. real, allocatable, dimension(:,:,:) :: tt, qv, hgt_3ddata
  81. logical :: is_psfc
  82. type (met_data) :: ecmwf_data, p_data, rh_data, hgt_data
  83. !
  84. ! Setup (read namelist and check on time range)
  85. !
  86. call get_namelist_params()
  87. call set_debug_level(WARN)
  88. ! Compute number of times that we will process
  89. call geth_idts(end_date(1), start_date(1), idiff)
  90. call mprintf((idiff < 0),ERROR,'Ending date is earlier than starting date in namelist for domain %i.', i1=1)
  91. n_times = idiff / interval_seconds
  92. ! Check that the interval evenly divides the range of times to process
  93. call mprintf((mod(idiff, interval_seconds) /= 0),WARN, &
  94. 'In namelist, interval_seconds does not evenly divide '// &
  95. '(end_date - start_date) for domain %i. Only %i time periods '// &
  96. 'will be processed.', i1=1, i2=n_times)
  97. fg_idx = 1
  98. input_name = fg_name(fg_idx)
  99. !
  100. ! Get coefficients for model level pressures
  101. !
  102. call read_coeffs()
  103. !
  104. ! Loop over all prefixes listed in namelist for fg_name
  105. !
  106. do while (input_name /= '*')
  107. !
  108. ! Loop over all times to be processed for this domain
  109. !
  110. do t=0,n_times
  111. ! Get the date string for the current time
  112. call geth_newdate(valid_date, trim(start_date(1)), t*interval_seconds)
  113. temp_date = ' '
  114. write(temp_date,'(a19)') valid_date(1:10)//'_'//valid_date(12:19)
  115. ! Initialize the module for reading in the met fields
  116. call read_met_init(trim(input_name), .false., temp_date(1:13), istatus)
  117. is_psfc = .false.
  118. if (istatus == 0) then
  119. call mprintf(.true.,STDOUT,'Reading from %s at time %s', s1=input_name, s2=temp_date(1:13))
  120. ! Process all fields and levels from the current file; read_next_met_field()
  121. ! will return a non-zero status when there are no more fields to be read.
  122. do while (istatus == 0)
  123. call read_next_met_field(ecmwf_data, istatus)
  124. if (istatus == 0) then
  125. ! Have we found either PSFC or LOGSFP?
  126. if ((trim(ecmwf_data%field) == 'PSFC' .and. ecmwf_data%xlvl == 200100.) &
  127. .or. trim(ecmwf_data%field) == 'LOGSFP') then
  128. p_data = ecmwf_data
  129. p_data%field = 'PRESSURE '
  130. p_data%desc = 'Pressure'
  131. rh_data = ecmwf_data
  132. rh_data%field = 'RH '
  133. rh_data%units = '%'
  134. rh_data%desc = 'Relative humidity'
  135. if (.not. allocated(psfc)) then
  136. allocate(psfc(ecmwf_data%nx,ecmwf_data%ny))
  137. end if
  138. call mprintf(.true.,STDOUT,'Found %s field in %s:%s', &
  139. s1=ecmwf_data%field, s2=input_name, s3=temp_date(1:13))
  140. is_psfc = .true.
  141. if (trim(ecmwf_data%field) == 'LOGSFP') then
  142. psfc(:,:) = exp(ecmwf_data%slab(:,:))
  143. else
  144. psfc(:,:) = ecmwf_data%slab(:,:)
  145. end if
  146. !CvH_DvD: - Store geopotential height
  147. else if (trim(ecmwf_data%field) == 'SOILHGT' .and. ecmwf_data%xlvl == 200100.) then
  148. hgt_data = ecmwf_data
  149. hgt_data%field = 'GHT '
  150. hgt_data%units = 'm'
  151. hgt_data%desc = 'Height'
  152. if (.not. allocated(hgtsfc)) then
  153. allocate(hgtsfc(ecmwf_data%nx,ecmwf_data%ny))
  154. end if
  155. call mprintf(.true.,STDOUT,'Found %s field in %s:%s', &
  156. s1=ecmwf_data%field, s2=input_name, s3=temp_date(1:13))
  157. hgtsfc(:,:) = ecmwf_data%slab(:,:)
  158. ! Have we found surface geopotential?
  159. else if (trim(ecmwf_data%field) == 'SOILGEO' .and. ecmwf_data%xlvl == 1.) then
  160. hgt_data = ecmwf_data
  161. hgt_data%field = 'GHT '
  162. hgt_data%units = 'm' ! units on output after conversion to height
  163. hgt_data%desc = 'Height'
  164. if (.not. allocated(hgtsfc)) then
  165. allocate(hgtsfc(ecmwf_data%nx,ecmwf_data%ny))
  166. end if
  167. call mprintf(.true.,STDOUT,'Found %s field in %s:%s', &
  168. s1=ecmwf_data%field, s2=input_name, s3=temp_date(1:13))
  169. hgtsfc(:,:) = ecmwf_data%slab(:,:)/9.81
  170. ! Have we found temperature?
  171. else if (trim(ecmwf_data%field) == 'TT') then
  172. if (.not. allocated(tt)) then
  173. allocate(tt(ecmwf_data%nx,ecmwf_data%ny,n_levels+1)) ! Extra level is for surface
  174. end if
  175. if (nint(ecmwf_data%xlvl) >= 1 .and. &
  176. nint(ecmwf_data%xlvl) <= n_levels) then
  177. tt(:,:,nint(ecmwf_data%xlvl)) = ecmwf_data%slab
  178. else if (nint(ecmwf_data%xlvl) == 200100) then
  179. tt(:,:,n_levels+1) = ecmwf_data%slab
  180. end if
  181. ! Have we found specific humidity?
  182. else if (trim(ecmwf_data%field) == 'SPECHUMD') then
  183. if (.not. allocated(qv)) then
  184. allocate(qv(ecmwf_data%nx,ecmwf_data%ny,n_levels+1)) ! Extra level is for surface
  185. end if
  186. if (nint(ecmwf_data%xlvl) >= 1 .and. &
  187. nint(ecmwf_data%xlvl) <= n_levels) then
  188. qv(:,:,nint(ecmwf_data%xlvl)) = ecmwf_data%slab
  189. else if (nint(ecmwf_data%xlvl) == 200100) then
  190. qv(:,:,n_levels+1) = ecmwf_data%slab
  191. end if
  192. end if
  193. if (associated(ecmwf_data%slab)) deallocate(ecmwf_data%slab)
  194. end if
  195. end do
  196. call read_met_close()
  197. ! Now write out, for each level, the pressure field
  198. if (is_psfc) then
  199. allocate(p_data%slab(p_data%nx,p_data%ny))
  200. allocate(rh_data%slab(rh_data%nx,rh_data%ny))
  201. !CvH_DvD: add HGT variable
  202. allocate(hgt_data%slab(hgt_data%nx,hgt_data%ny))
  203. call write_met_init(trim(get_path(input_name))//'PRES', .false., temp_date(1:13), istatus)
  204. if (allocated(tt) .and. allocated(qv)) then
  205. p_data%xlvl = 200100.
  206. p_data%slab = psfc
  207. ! Surface RH should be computed from surface DEWPT by ungrib
  208. ! rh_data%xlvl = 200100.
  209. ! call calc_rh(tt(:,:,n_levels+1), qv(:,:,n_levels+1), psfc, rh_data%slab, rh_data%nx, rh_data%ny)
  210. ! call write_next_met_field(rh_data, istatus)
  211. call write_next_met_field(p_data, istatus)
  212. else
  213. call mprintf(.true.,WARN,'Either TT or SPECHUMD not found. No RH will be computed.')
  214. end if
  215. !CvH_DvD: if tt, qv and hgtsfc are available compute hgt
  216. if (allocated(tt) .and. allocated(qv) .and. allocated(hgtsfc)) then
  217. if (.not. allocated(hgtprev)) then
  218. allocate(hgtprev(ecmwf_data%nx,ecmwf_data%ny))
  219. end if
  220. if (.not. allocated(pstart)) then
  221. allocate(pstart(ecmwf_data%nx,ecmwf_data%ny))
  222. end if
  223. if (.not. allocated(pend)) then
  224. allocate(pend(ecmwf_data%nx,ecmwf_data%ny))
  225. end if
  226. if (.not. allocated(hgt_3ddata)) then
  227. allocate(hgt_3ddata(ecmwf_data%nx,ecmwf_data%ny, n_levels))
  228. end if
  229. ! CvH_DvD: interpolate Q and T if they are not availalbe at all levels
  230. do i = n_levels, 1, -1
  231. if (tt(1,1,i) .eq. 0) then ! CvH_DvD: If TT is zero, level is unknown, so moisture is zero too
  232. if (i .ne. n_levels) then
  233. counter=1
  234. ! CvH_DvD: Find first available level
  235. do while ((tt(1,1,i-counter) .eq. 0) .and. (i-counter .gt. 1))
  236. counter=counter+1
  237. end do
  238. if (tt(1,1,i-counter) .gt. 0) then
  239. ! CvH_DvD: If TT is zero, we're at the top level which is missing, so use previous level,
  240. ! mod level will remove it anyway
  241. tt(:,:,i)=tt(:,:,i+1)+(tt(:,:,i-counter)-tt(:,:,i+1))/(counter+1)
  242. qv(:,:,i)=qv(:,:,i+1)+(qv(:,:,i-counter)-qv(:,:,i+1))/(counter+1)
  243. else
  244. tt(:,:,i)=tt(:,:,i+1)
  245. qv(:,:,i)=qv(:,:,i+1)
  246. end if
  247. else
  248. write(6,*) 'WARNING First level is missing!' ! CvH_DvD: First level should be there!
  249. end if
  250. end if
  251. end do
  252. ! CvH_DvD: first previous hgt is surface hgt/soilgeo
  253. hgtprev = hgtsfc
  254. ! CvH_DvD: - Loop from surface to top, note: 1=top & n_levels=surface
  255. do i = n_levels, 1, -1
  256. ! CvH_DvD: compute half level pressure at current half-level
  257. pend = 0.5 * (a(i-1) + a(i)) + psfc * 0.5 * (b(i-1) + b(i))
  258. if (i .eq. n_levels) then
  259. ! CvH_DvD: use Tv not T, use Psfc and T_halflevel_1=AveT lowest level
  260. hgt_3ddata(:,:,i) = hgtprev + 287.05 * tt(:,:,i) * (1. + 0.61 * qv(:,:,i)) * log(psfc/pend) / 9.81
  261. else
  262. ! CvH_DvD: compute half level pressure beneath current half-level
  263. pstart = 0.5 * (a(i+1) + a(i)) + psfc * 0.5 * (b(i+1) + b(i))
  264. ! CvH_DvD: use Tv not T, create T at full leve beneath current half level
  265. hgt_3ddata(:,:,i) = hgtprev + 287.05 * (tt(:,:,i) + tt(:,:,i+1))/2. * &
  266. (1. + 0.61 * (qv(:,:,i)+qv(:,:,i+1))/2.) * log(pstart/pend) / 9.81
  267. end if
  268. hgtprev = hgt_3ddata(:,:,i)
  269. end do
  270. end if
  271. do i = 1, n_levels
  272. a_full = 0.5 * (a(i-1) + a(i)) ! A and B are dimensioned (0:n_levels)
  273. b_full = 0.5 * (b(i-1) + b(i))
  274. p_data%xlvl = real(i)
  275. p_data%slab = a_full + psfc * b_full
  276. if (allocated(tt) .and. allocated(qv)) then
  277. rh_data%xlvl = real(i)
  278. call calc_rh(tt(:,:,i), qv(:,:,i), p_data%slab, rh_data%slab, rh_data%nx, rh_data%ny)
  279. call write_next_met_field(rh_data, istatus)
  280. if (allocated(hgtsfc)) then
  281. ! CvH_DvD: put hgt_3ddata into hgt_data object
  282. hgt_data%xlvl = real(i)
  283. hgt_data%slab = hgt_3ddata(:,:,i)
  284. call write_next_met_field(hgt_data, istatus)
  285. else
  286. call mprintf(.true., WARN, &
  287. 'Either SOILHGT or SOILGEO are required to create 3-d GHT field, which is required '// &
  288. 'for a correct vertical interpolation in real.')
  289. end if
  290. end if
  291. call write_next_met_field(p_data, istatus)
  292. end do
  293. call write_met_close()
  294. deallocate(p_data%slab)
  295. deallocate(rh_data%slab)
  296. ! CvH_DvD: deallocate stuff
  297. deallocate(hgt_data%slab)
  298. end if
  299. if (allocated(tt)) deallocate(tt)
  300. if (allocated(qv)) deallocate(qv)
  301. if (allocated(psfc)) deallocate(psfc)
  302. ! CvH_DvD: deallocate stuff
  303. if (allocated(hgt_3ddata)) deallocate(hgt_3ddata)
  304. if (allocated(pstart)) deallocate(pstart)
  305. if (allocated(pend)) deallocate(pend)
  306. if (allocated(hgtprev)) deallocate(hgtprev)
  307. else
  308. call mprintf(.true.,ERROR,'Problem opening %s at time %s', s1=input_name, s2=temp_date(1:13))
  309. end if
  310. end do
  311. ! CvH_DvD: only now deallocate hgtsfc,
  312. ! because oper EC does not have hgtsfc @ fps
  313. if (allocated(hgtsfc)) deallocate(hgtsfc)
  314. fg_idx = fg_idx + 1
  315. input_name = fg_name(fg_idx)
  316. end do
  317. call cleanup_coeffs()
  318. stop
  319. end program calc_ecmwf_p
  320. subroutine calc_rh(t, qv, p, rh, nx, ny)
  321. implicit none
  322. ! Arguments
  323. integer, intent(in) :: nx, ny
  324. real, dimension(nx, ny), intent(in) :: t, qv, p
  325. real, dimension(nx, ny), intent(out) :: rh
  326. ! Constants
  327. real, parameter :: svp1=0.6112
  328. real, parameter :: svp2=17.67
  329. real, parameter :: svp3=29.65
  330. real, parameter :: svpt0=273.15
  331. real, parameter :: eps = 0.622
  332. ! Local variables
  333. integer :: i, j
  334. real :: es, e
  335. do j=1,ny
  336. do i=1,nx
  337. es=svp1*10.*exp(svp2*(t(i,j)-svpt0)/(t(i,j)-svp3))
  338. e=qv(i,j)*p(i,j)/100./(eps+qv(i,j)*(1.-eps)) ! qv is specific humidity
  339. rh(i,j) = 100.0 * e/es
  340. end do
  341. end do
  342. end subroutine calc_rh