PageRenderTime 60ms CodeModel.GetById 22ms RepoModel.GetById 0ms app.codeStats 1ms

/wrfv2_fire/external/io_grib1/io_grib1.F

http://github.com/jbeezley/wrf-fire
FORTRAN Legacy | 3560 lines | 2077 code | 693 blank | 790 comment | 215 complexity | 711faad64316b7f9509963d31eec1ed4 MD5 | raw file
Possible License(s): AGPL-1.0

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

  1. !*-----------------------------------------------------------------------------
  2. !*
  3. !* Todd Hutchinson
  4. !* WSI
  5. !* 400 Minuteman Road
  6. !* Andover, MA 01810
  7. !* thutchinson@wsi.com
  8. !*
  9. !*-----------------------------------------------------------------------------
  10. !*
  11. !* This io_grib1 API is designed to read WRF input and write WRF output data
  12. !* in grib version 1 format.
  13. !*
  14. module gr1_data_info
  15. !*
  16. !* This module will hold data internal to this I/O implementation.
  17. !* The variables will be accessible by all functions (provided they have a
  18. !* "USE gr1_data_info" line).
  19. !*
  20. integer , parameter :: FATAL = 1
  21. integer , parameter :: DEBUG = 100
  22. integer , parameter :: DateStrLen = 19
  23. integer , parameter :: firstFileHandle = 8
  24. integer , parameter :: maxFileHandles = 30
  25. integer , parameter :: maxLevels = 1000
  26. integer , parameter :: maxSoilLevels = 100
  27. integer , parameter :: maxDomains = 500
  28. logical , dimension(maxFileHandles) :: committed, opened, used
  29. character*128, dimension(maxFileHandles) :: DataFile
  30. integer, dimension(maxFileHandles) :: FileFd
  31. integer, dimension(maxFileHandles) :: FileStatus
  32. REAL, dimension(maxLevels) :: half_eta, full_eta
  33. REAL, dimension(maxSoilLevels) :: soil_depth, soil_thickness
  34. character*24 :: StartDate = ''
  35. character*24 :: InputProgramName = ''
  36. integer :: projection
  37. integer :: wg_grid_id
  38. real :: dx,dy
  39. real :: truelat1, truelat2
  40. real :: center_lat, center_lon
  41. real :: proj_central_lon
  42. real :: timestep
  43. character, dimension(:), pointer :: grib_tables
  44. logical :: table_filled = .FALSE.
  45. character, dimension(:), pointer :: grid_info
  46. integer :: full_xsize, full_ysize
  47. integer, dimension(maxDomains) :: domains = -1
  48. integer :: this_domain = 0
  49. integer :: max_domain = 0
  50. TYPE :: HandleVar
  51. character, dimension(:), pointer :: fileindex(:)
  52. integer :: CurrentTime
  53. integer :: NumberTimes
  54. character (DateStrLen), dimension(:),pointer :: Times(:)
  55. ENDTYPE
  56. TYPE (HandleVar), dimension(maxFileHandles) :: fileinfo
  57. TYPE :: prevdata
  58. integer :: fcst_secs_rainc
  59. integer :: fcst_secs_rainnc
  60. real, dimension(:,:), pointer :: rainc, rainnc
  61. END TYPE prevdata
  62. TYPE (prevdata), DIMENSION(500) :: lastdata
  63. TYPE :: initdata
  64. real, dimension(:,:), pointer :: snod
  65. END TYPE initdata
  66. TYPE (initdata), dimension(maxDomains) :: firstdata
  67. TYPE :: prestype
  68. real, dimension(:,:,:), pointer :: vals
  69. logical :: newtime
  70. character*120 :: lastDateStr
  71. END TYPE prestype
  72. character*120, dimension(maxDomains) :: lastDateStr
  73. TYPE (prestype), dimension(maxDomains) :: pressure
  74. TYPE (prestype), dimension(maxDomains) :: geopotential
  75. integer :: center, subcenter, parmtbl
  76. character(len=15000), dimension(firstFileHandle:maxFileHandles) :: td_output
  77. character(len=15000), dimension(firstFileHandle:maxFileHandles) :: ti_output
  78. logical :: WrfIOnotInitialized = .true.
  79. end module gr1_data_info
  80. subroutine ext_gr1_ioinit(SysDepInfo,Status)
  81. USE gr1_data_info
  82. implicit none
  83. #include "wrf_status_codes.h"
  84. #include "wrf_io_flags.h"
  85. CHARACTER*(*), INTENT(IN) :: SysDepInfo
  86. integer ,intent(out) :: Status
  87. integer :: i
  88. integer :: size, istat
  89. CHARACTER (LEN=300) :: wrf_err_message
  90. call wrf_debug ( DEBUG , 'Entering ext_gr1_ioinit')
  91. do i=firstFileHandle, maxFileHandles
  92. used(i) = .false.
  93. committed(i) = .false.
  94. opened(i) = .false.
  95. td_output(i) = ''
  96. ti_output(i) = ''
  97. enddo
  98. domains(:) = -1
  99. do i = 1, maxDomains
  100. pressure(i)%newtime = .false.
  101. pressure(i)%lastDateStr = ''
  102. geopotential(i)%newtime = .false.
  103. geopotential(i)%lastDateStr = ''
  104. lastDateStr(i) = ''
  105. enddo
  106. lastdata%fcst_secs_rainc = 0
  107. lastdata%fcst_secs_rainnc = 0
  108. FileStatus(1:maxFileHandles) = WRF_FILE_NOT_OPENED
  109. WrfIOnotInitialized = .false.
  110. Status = WRF_NO_ERR
  111. return
  112. end subroutine ext_gr1_ioinit
  113. !*****************************************************************************
  114. subroutine ext_gr1_ioexit(Status)
  115. USE gr1_data_info
  116. implicit none
  117. #include "wrf_status_codes.h"
  118. integer istat
  119. integer ,intent(out) :: Status
  120. call wrf_debug ( DEBUG , 'Entering ext_gr1_ioexit')
  121. if (table_filled) then
  122. CALL free_gribmap(grib_tables)
  123. DEALLOCATE(grib_tables, stat=istat)
  124. table_filled = .FALSE.
  125. endif
  126. IF ( ASSOCIATED ( grid_info ) ) THEN
  127. DEALLOCATE(grid_info, stat=istat)
  128. ENDIF
  129. NULLIFY(grid_info)
  130. Status = WRF_NO_ERR
  131. return
  132. end subroutine ext_gr1_ioexit
  133. !*****************************************************************************
  134. SUBROUTINE ext_gr1_open_for_read_begin ( FileName , Comm_compute, Comm_io, &
  135. SysDepInfo, DataHandle , Status )
  136. USE gr1_data_info
  137. IMPLICIT NONE
  138. #include "wrf_status_codes.h"
  139. #include "wrf_io_flags.h"
  140. CHARACTER*(*) :: FileName
  141. INTEGER , INTENT(IN) :: Comm_compute , Comm_io
  142. CHARACTER*(*) :: SysDepInfo
  143. INTEGER , INTENT(OUT) :: DataHandle
  144. INTEGER , INTENT(OUT) :: Status
  145. integer :: ierr
  146. integer :: size
  147. integer :: idx
  148. integer :: parmid
  149. integer :: dpth_parmid
  150. integer :: thk_parmid
  151. integer :: leveltype
  152. integer , DIMENSION(1000) :: indices
  153. integer :: numindices
  154. real , DIMENSION(1000) :: levels
  155. real :: tmp
  156. integer :: swapped
  157. integer :: etaidx
  158. integer :: grb_index
  159. integer :: level1, level2
  160. integer :: tablenum
  161. integer :: stat
  162. integer :: endchar
  163. integer :: last_grb_index
  164. CHARACTER (LEN=300) :: wrf_err_message
  165. call wrf_debug ( DEBUG , 'Entering ext_gr1_open_for_read_begin')
  166. CALL gr1_get_new_handle(DataHandle)
  167. if (DataHandle .GT. 0) then
  168. CALL open_file(TRIM(FileName), 'r', FileFd(DataHandle), ierr)
  169. if (ierr .ne. 0) then
  170. Status = WRF_ERR_FATAL_BAD_FILE_STATUS
  171. else
  172. opened(DataHandle) = .true.
  173. DataFile(DataHandle) = TRIM(FileName)
  174. FileStatus(DataHandle) = WRF_FILE_OPENED_NOT_COMMITTED
  175. endif
  176. else
  177. Status = WRF_WARN_TOO_MANY_FILES
  178. return
  179. endif
  180. ! Read the grib index file first
  181. if (.NOT. table_filled) then
  182. table_filled = .TRUE.
  183. CALL GET_GRIB1_TABLES_SIZE(size)
  184. ALLOCATE(grib_tables(1:size), STAT=ierr)
  185. CALL LOAD_GRIB1_TABLES ("gribmap.txt", grib_tables, ierr)
  186. if (ierr .ne. 0) then
  187. DEALLOCATE(grib_tables)
  188. WRITE( wrf_err_message , * ) &
  189. 'Could not open file gribmap.txt '
  190. CALL wrf_error_fatal ( TRIM ( wrf_err_message ) )
  191. Status = WRF_ERR_FATAL_BAD_FILE_STATUS
  192. return
  193. endif
  194. endif
  195. ! Begin by indexing file and reading metadata into structure.
  196. CALL GET_FILEINDEX_SIZE(size)
  197. ALLOCATE(fileinfo(DataHandle)%fileindex(1:size), STAT=ierr)
  198. CALL ALLOC_INDEX_FILE(fileinfo(DataHandle)%fileindex(:))
  199. CALL INDEX_FILE(FileFd(DataHandle),fileinfo(DataHandle)%fileindex(:))
  200. ! Get times into Times variable
  201. CALL GET_NUM_TIMES(fileinfo(DataHandle)%fileindex(:), &
  202. fileinfo(DataHandle)%NumberTimes);
  203. ALLOCATE(fileinfo(DataHandle)%Times(1:fileinfo(DataHandle)%NumberTimes), STAT=ierr)
  204. do idx = 1,fileinfo(DataHandle)%NumberTimes
  205. CALL GET_TIME(fileinfo(DataHandle)%fileindex(:),idx, &
  206. fileinfo(DataHandle)%Times(idx))
  207. enddo
  208. ! CurrentTime starts as 0. The first time in the file is 1. So,
  209. ! until set_time or get_next_time is called, the current time
  210. ! is not set.
  211. fileinfo(DataHandle)%CurrentTime = 0
  212. CALL gr1_fill_eta_levels(fileinfo(DataHandle)%fileindex(:), &
  213. FileFd(DataHandle), &
  214. grib_tables, "ZNW", full_eta)
  215. CALL gr1_fill_eta_levels(fileinfo(DataHandle)%fileindex(:), FileFd(DataHandle), &
  216. grib_tables, "ZNU", half_eta)
  217. !
  218. ! Now, get the soil levels
  219. !
  220. CALL GET_GRIB_PARAM(grib_tables, "ZS", center, subcenter, parmtbl, &
  221. tablenum, dpth_parmid)
  222. CALL GET_GRIB_PARAM(grib_tables,"DZS", center, subcenter, parmtbl, &
  223. tablenum, thk_parmid)
  224. if (dpth_parmid == -1) then
  225. call wrf_message ('Error getting grib parameter')
  226. endif
  227. leveltype = 112
  228. CALL GET_GRIB_INDICES(fileinfo(DataHandle)%fileindex(:),center, subcenter, parmtbl, &
  229. dpth_parmid,"*",leveltype, &
  230. -HUGE(1),-HUGE(1), -HUGE(1),-HUGE(1),indices,numindices)
  231. last_grb_index = -1;
  232. do idx = 1,numindices
  233. CALL READ_GRIB(fileinfo(DataHandle)%fileindex(:), FileFd(DataHandle), &
  234. indices(idx), soil_depth(idx))
  235. !
  236. ! Now read the soil thickenesses
  237. !
  238. CALL GET_LEVEL1(fileinfo(DataHandle)%fileindex(:),indices(idx),level1)
  239. CALL GET_LEVEL2(fileinfo(DataHandle)%fileindex(:),indices(idx),level2)
  240. CALL GET_GRIB_INDEX_GUESS(fileinfo(DataHandle)%fileindex(:), &
  241. center, subcenter, parmtbl, thk_parmid,"*",leveltype, &
  242. level1,level2,-HUGE(1),-HUGE(1), last_grb_index+1, grb_index)
  243. CALL READ_GRIB(fileinfo(DataHandle)%fileindex(:),FileFd(DataHandle),grb_index, &
  244. soil_thickness(idx))
  245. last_grb_index = grb_index
  246. enddo
  247. !
  248. ! Fill up any variables that need to be retrieved from Metadata
  249. !
  250. CALL GET_METADATA_VALUE(fileinfo(DataHandle)%fileindex(:), 'PROGRAM_NAME', "none", &
  251. "none", InputProgramName, stat)
  252. if (stat /= 0) then
  253. CALL wrf_debug (DEBUG , "PROGRAM_NAME not found in input METADATA")
  254. else
  255. endchar = SCAN(InputProgramName," ")
  256. InputProgramName = InputProgramName(1:endchar)
  257. endif
  258. call wrf_debug ( DEBUG , 'Exiting ext_gr1_open_for_read_begin')
  259. RETURN
  260. END SUBROUTINE ext_gr1_open_for_read_begin
  261. !*****************************************************************************
  262. SUBROUTINE ext_gr1_open_for_read_commit( DataHandle , Status )
  263. USE gr1_data_info
  264. IMPLICIT NONE
  265. #include "wrf_status_codes.h"
  266. #include "wrf_io_flags.h"
  267. character(len=1000) :: msg
  268. INTEGER , INTENT(IN ) :: DataHandle
  269. INTEGER , INTENT(OUT) :: Status
  270. call wrf_debug ( DEBUG , 'Entering ext_gr1_open_for_read_commit')
  271. Status = WRF_NO_ERR
  272. if(WrfIOnotInitialized) then
  273. Status = WRF_IO_NOT_INITIALIZED
  274. write(msg,*) 'ext_gr1_ioinit was not called ',__FILE__,', line', __LINE__
  275. call wrf_debug ( FATAL , msg)
  276. return
  277. endif
  278. committed(DataHandle) = .true.
  279. FileStatus(DataHandle) = WRF_FILE_OPENED_FOR_READ
  280. Status = WRF_NO_ERR
  281. RETURN
  282. END SUBROUTINE ext_gr1_open_for_read_commit
  283. !*****************************************************************************
  284. SUBROUTINE ext_gr1_open_for_read ( FileName , Comm_compute, Comm_io, &
  285. SysDepInfo, DataHandle , Status )
  286. USE gr1_data_info
  287. IMPLICIT NONE
  288. #include "wrf_status_codes.h"
  289. #include "wrf_io_flags.h"
  290. CHARACTER*(*) :: FileName
  291. INTEGER , INTENT(IN) :: Comm_compute , Comm_io
  292. CHARACTER*(*) :: SysDepInfo
  293. INTEGER , INTENT(OUT) :: DataHandle
  294. INTEGER , INTENT(OUT) :: Status
  295. call wrf_debug ( DEBUG , 'Entering ext_gr1_open_for_read')
  296. DataHandle = 0 ! dummy setting to quiet warning message
  297. CALL ext_gr1_open_for_read_begin( FileName, Comm_compute, Comm_io, &
  298. SysDepInfo, DataHandle, Status )
  299. IF ( Status .EQ. WRF_NO_ERR ) THEN
  300. FileStatus(DataHandle) = WRF_FILE_OPENED_NOT_COMMITTED
  301. CALL ext_gr1_open_for_read_commit( DataHandle, Status )
  302. ENDIF
  303. return
  304. RETURN
  305. END SUBROUTINE ext_gr1_open_for_read
  306. !*****************************************************************************
  307. SUBROUTINE ext_gr1_open_for_write_begin(FileName, Comm, IOComm, SysDepInfo, &
  308. DataHandle, Status)
  309. USE gr1_data_info
  310. implicit none
  311. #include "wrf_status_codes.h"
  312. #include "wrf_io_flags.h"
  313. character*(*) ,intent(in) :: FileName
  314. integer ,intent(in) :: Comm
  315. integer ,intent(in) :: IOComm
  316. character*(*) ,intent(in) :: SysDepInfo
  317. integer ,intent(out) :: DataHandle
  318. integer ,intent(out) :: Status
  319. integer :: ierr
  320. CHARACTER (LEN=300) :: wrf_err_message
  321. integer :: size
  322. call wrf_debug ( DEBUG , 'Entering ext_gr1_open_for_write_begin')
  323. if (.NOT. table_filled) then
  324. table_filled = .TRUE.
  325. CALL GET_GRIB1_TABLES_SIZE(size)
  326. ALLOCATE(grib_tables(1:size), STAT=ierr)
  327. CALL LOAD_GRIB1_TABLES ("gribmap.txt", grib_tables, ierr)
  328. if (ierr .ne. 0) then
  329. DEALLOCATE(grib_tables)
  330. WRITE( wrf_err_message , * ) &
  331. 'Could not open file gribmap.txt '
  332. CALL wrf_error_fatal ( TRIM ( wrf_err_message ) )
  333. Status = WRF_ERR_FATAL_BAD_FILE_STATUS
  334. return
  335. endif
  336. endif
  337. Status = WRF_NO_ERR
  338. CALL gr1_get_new_handle(DataHandle)
  339. if (DataHandle .GT. 0) then
  340. CALL open_file(TRIM(FileName), 'w', FileFd(DataHandle), ierr)
  341. if (ierr .ne. 0) then
  342. Status = WRF_WARN_WRITE_RONLY_FILE
  343. else
  344. opened(DataHandle) = .true.
  345. DataFile(DataHandle) = TRIM(FileName)
  346. FileStatus(DataHandle) = WRF_FILE_OPENED_NOT_COMMITTED
  347. endif
  348. committed(DataHandle) = .false.
  349. td_output(DataHandle) = ''
  350. else
  351. Status = WRF_WARN_TOO_MANY_FILES
  352. endif
  353. RETURN
  354. END SUBROUTINE ext_gr1_open_for_write_begin
  355. !*****************************************************************************
  356. SUBROUTINE ext_gr1_open_for_write_commit( DataHandle , Status )
  357. USE gr1_data_info
  358. IMPLICIT NONE
  359. #include "wrf_status_codes.h"
  360. #include "wrf_io_flags.h"
  361. INTEGER , INTENT(IN ) :: DataHandle
  362. INTEGER , INTENT(OUT) :: Status
  363. call wrf_debug ( DEBUG , 'Entering ext_gr1_open_for_write_commit')
  364. IF ( opened( DataHandle ) ) THEN
  365. IF ( used( DataHandle ) ) THEN
  366. committed(DataHandle) = .true.
  367. FileStatus(DataHandle) = WRF_FILE_OPENED_FOR_WRITE
  368. ENDIF
  369. ENDIF
  370. Status = WRF_NO_ERR
  371. RETURN
  372. END SUBROUTINE ext_gr1_open_for_write_commit
  373. !*****************************************************************************
  374. subroutine ext_gr1_inquiry (Inquiry, Result, Status)
  375. use gr1_data_info
  376. implicit none
  377. #include "wrf_status_codes.h"
  378. character *(*), INTENT(IN) :: Inquiry
  379. character *(*), INTENT(OUT) :: Result
  380. integer ,INTENT(INOUT) :: Status
  381. SELECT CASE (Inquiry)
  382. CASE ("RANDOM_WRITE","RANDOM_READ")
  383. Result='ALLOW'
  384. CASE ("SEQUENTIAL_WRITE","SEQUENTIAL_READ")
  385. Result='NO'
  386. CASE ("OPEN_READ", "OPEN_WRITE", "OPEN_COMMIT_WRITE")
  387. Result='REQUIRE'
  388. CASE ("OPEN_COMMIT_READ","PARALLEL_IO")
  389. Result='NO'
  390. CASE ("SELF_DESCRIBING","SUPPORT_METADATA","SUPPORT_3D_FIELDS")
  391. Result='YES'
  392. CASE ("MEDIUM")
  393. Result ='FILE'
  394. CASE DEFAULT
  395. Result = 'No Result for that inquiry!'
  396. END SELECT
  397. Status=WRF_NO_ERR
  398. return
  399. end subroutine ext_gr1_inquiry
  400. !*****************************************************************************
  401. SUBROUTINE ext_gr1_inquire_opened ( DataHandle, FileName , FileStat, Status )
  402. USE gr1_data_info
  403. IMPLICIT NONE
  404. #include "wrf_status_codes.h"
  405. #include "wrf_io_flags.h"
  406. INTEGER , INTENT(IN) :: DataHandle
  407. CHARACTER*(*) :: FileName
  408. INTEGER , INTENT(OUT) :: FileStat
  409. INTEGER , INTENT(OUT) :: Status
  410. call wrf_debug ( DEBUG , 'Entering ext_gr1_inquire_opened')
  411. FileStat = WRF_NO_ERR
  412. if ((DataHandle .ge. firstFileHandle) .and. &
  413. (DataHandle .le. maxFileHandles)) then
  414. FileStat = FileStatus(DataHandle)
  415. else
  416. FileStat = WRF_FILE_NOT_OPENED
  417. endif
  418. Status = FileStat
  419. RETURN
  420. END SUBROUTINE ext_gr1_inquire_opened
  421. !*****************************************************************************
  422. SUBROUTINE ext_gr1_ioclose ( DataHandle, Status )
  423. USE gr1_data_info
  424. IMPLICIT NONE
  425. #include "wrf_status_codes.h"
  426. INTEGER DataHandle, Status
  427. INTEGER istat
  428. INTEGER ierr
  429. character(len=1000) :: outstring
  430. character :: lf
  431. lf=char(10)
  432. call wrf_debug ( DEBUG , 'Entering ext_gr1_ioclose')
  433. Status = WRF_NO_ERR
  434. CALL write_file(FileFd(DataHandle), lf//'<METADATA>'//lf,ierr)
  435. outstring = &
  436. '<!-- The following are fields that were supplied to the WRF I/O API.'//lf//&
  437. 'Many variables (but not all) are redundant with the variables within '//lf//&
  438. 'the grib headers. They are stored here, as METADATA, so that the '//lf//&
  439. 'WRF I/O API has simple access to these variables.-->'
  440. CALL write_file(FileFd(DataHandle), trim(outstring), ierr)
  441. if (trim(ti_output(DataHandle)) /= '') then
  442. CALL write_file(FileFd(DataHandle), trim(ti_output(DataHandle)), ierr)
  443. CALL write_file(FileFd(DataHandle), lf, ierr)
  444. endif
  445. if (trim(td_output(DataHandle)) /= '') then
  446. CALL write_file(FileFd(DataHandle), trim(td_output(DataHandle)), ierr)
  447. CALL write_file(FileFd(DataHandle), lf, ierr)
  448. endif
  449. CALL write_file(FileFd(DataHandle), '</METADATA>'//lf,ierr)
  450. ti_output(DataHandle) = ''
  451. td_output(DataHandle) = ''
  452. if (ierr .ne. 0) then
  453. Status = WRF_WARN_WRITE_RONLY_FILE
  454. endif
  455. CALL close_file(FileFd(DataHandle))
  456. used(DataHandle) = .false.
  457. RETURN
  458. END SUBROUTINE ext_gr1_ioclose
  459. !*****************************************************************************
  460. SUBROUTINE ext_gr1_write_field( DataHandle , DateStrIn , VarName , &
  461. Field , FieldType , Comm , IOComm, &
  462. DomainDesc , MemoryOrder , Stagger , &
  463. DimNames , &
  464. DomainStart , DomainEnd , &
  465. MemoryStart , MemoryEnd , &
  466. PatchStart , PatchEnd , &
  467. Status )
  468. USE gr1_data_info
  469. IMPLICIT NONE
  470. #include "wrf_status_codes.h"
  471. #include "wrf_io_flags.h"
  472. #include "wrf_projection.h"
  473. INTEGER , INTENT(IN) :: DataHandle
  474. CHARACTER*(*) :: DateStrIn
  475. CHARACTER(DateStrLen) :: DateStr
  476. CHARACTER*(*) :: VarName
  477. CHARACTER*120 :: OutName
  478. CHARACTER(120) :: TmpVarName
  479. integer ,intent(in) :: FieldType
  480. integer ,intent(inout) :: Comm
  481. integer ,intent(inout) :: IOComm
  482. integer ,intent(in) :: DomainDesc
  483. character*(*) ,intent(in) :: MemoryOrder
  484. character*(*) ,intent(in) :: Stagger
  485. character*(*) , dimension (*) ,intent(in) :: DimNames
  486. integer ,dimension(*) ,intent(in) :: DomainStart, DomainEnd
  487. integer ,dimension(*) ,intent(in) :: MemoryStart, MemoryEnd
  488. integer ,dimension(*) ,intent(in) :: PatchStart, PatchEnd
  489. integer ,intent(out) :: Status
  490. integer :: ierror
  491. character (120) :: msg
  492. integer :: xsize, ysize, zsize
  493. integer :: x, y, z
  494. integer :: x_start,x_end,y_start,y_end,z_start,z_end,ndim
  495. integer :: idx
  496. integer :: proj_center_flag
  497. logical :: vert_stag = .false.
  498. integer :: levelnum
  499. real, DIMENSION(:,:), POINTER :: data,tmpdata
  500. integer, DIMENSION(:), POINTER :: mold
  501. integer :: istat
  502. integer :: accum_period
  503. integer :: size
  504. integer, dimension(1000) :: level1, level2
  505. real, DIMENSION( 1:1,MemoryStart(1):MemoryEnd(1), &
  506. MemoryStart(2):MemoryEnd(2), &
  507. MemoryStart(3):MemoryEnd(3) ) :: Field
  508. real :: fcst_secs
  509. logical :: soil_layers, fraction
  510. integer :: vert_unit
  511. integer :: abc(2,2,2)
  512. integer :: def(8)
  513. logical :: output = .true.
  514. integer :: idx1, idx2, idx3
  515. logical :: new_domain
  516. real :: region_center_lat, region_center_lon
  517. integer :: dom_xsize, dom_ysize;
  518. integer :: ierr
  519. logical :: already_have_domain
  520. call wrf_debug ( DEBUG , 'Entering ext_gr1_write_field for parameter'//VarName)
  521. !
  522. ! If DateStr is all 0's, we reset it to StartDate (if StartDate exists).
  523. ! For some reason,
  524. ! in idealized simulations, StartDate is 0001-01-01_00:00:00 while
  525. ! the first DateStr is 0000-00-00_00:00:00.
  526. !
  527. if (DateStrIn .eq. '0000-00-00_00:00:00') then
  528. if (StartDate .ne. '') then
  529. DateStr = TRIM(StartDate)
  530. else
  531. DateStr = '0001-01-01_00:00:00'
  532. endif
  533. else
  534. DateStr = DateStrIn
  535. endif
  536. !
  537. ! Check if this is a domain that we haven't seen yet. If so, add it to
  538. ! the list of domains.
  539. !
  540. new_domain = .false.
  541. already_have_domain = .false.
  542. do idx = 1, max_domain
  543. if (this_domain .eq. domains(idx)) then
  544. already_have_domain = .true.
  545. endif
  546. enddo
  547. if (.NOT. already_have_domain) then
  548. max_domain = max_domain + 1
  549. domains(max_domain) = this_domain
  550. new_domain = .true.
  551. endif
  552. !
  553. ! If the time has changed, we open a new file. This is a kludge to get
  554. ! around slowness in WRF that occurs when opening a new data file the
  555. ! standard way.
  556. !
  557. #ifdef GRIB_ONE_TIME_PER_FILE
  558. if (lastDateStr(this_domain) .ne. DateStr) then
  559. write(DataFile(DataHandle),'(A8,i2.2,A1,A19)') 'wrfout_d',this_domain,'_',DateStr
  560. call ext_gr1_ioclose ( DataHandle, Status )
  561. CALL open_file(TRIM(DataFile(DataHandle)), 'w', FileFd(DataHandle), ierr)
  562. if (ierr .ne. 0) then
  563. print *,'Could not open new file: ',DataFile(DataHandle)
  564. print *,' Appending to old file.'
  565. else
  566. ! Just set used back to .true. here, since ioclose set it to false.
  567. used(DataHandle) = .true.
  568. endif
  569. td_output(DataHandle) = ''
  570. endif
  571. lastDateStr(this_domain) = DateStr
  572. #endif
  573. output = .true.
  574. zsize = 1
  575. xsize = 1
  576. ysize = 1
  577. OutName = VarName
  578. soil_layers = .false.
  579. fraction = .false.
  580. ! First, handle then special cases for the boundary data.
  581. CALL get_dims(MemoryOrder, PatchStart, PatchEnd, ndim, x_start, x_end, &
  582. y_start, y_end,z_start,z_end)
  583. xsize = x_end - x_start + 1
  584. ysize = y_end - y_start + 1
  585. zsize = z_end - z_start + 1
  586. do idx = 1, len(MemoryOrder)
  587. if ((MemoryOrder(idx:idx) .eq. 'Z') .and. &
  588. (DimNames(idx) .eq. 'soil_layers_stag')) then
  589. soil_layers = .true.
  590. else if ((OutName .eq. 'LANDUSEF') .or. (OutName .eq. 'SOILCBOT') .or. &
  591. (OutName .eq. 'SOILCTOP')) then
  592. fraction = .true.
  593. endif
  594. enddo
  595. if (.not. ASSOCIATED(grid_info)) then
  596. CALL get_grid_info_size(size)
  597. ALLOCATE(grid_info(1:size), STAT=istat)
  598. if (istat .eq. -1) then
  599. DEALLOCATE(grid_info)
  600. Status = WRF_ERR_FATAL_BAD_FILE_STATUS
  601. return
  602. endif
  603. endif
  604. if (new_domain) then
  605. ALLOCATE(firstdata(this_domain)%snod(xsize,ysize))
  606. firstdata(this_domain)%snod(:,:) = 0.0
  607. ALLOCATE(lastdata(this_domain)%rainc(xsize,ysize))
  608. lastdata(this_domain)%rainc(:,:) = 0.0
  609. ALLOCATE(lastdata(this_domain)%rainnc(xsize,ysize))
  610. lastdata(this_domain)%rainnc(:,:) = 0.0
  611. endif
  612. if (zsize .eq. 0) then
  613. zsize = 1
  614. endif
  615. ALLOCATE(data(1:xsize,1:ysize), STAT=istat)
  616. ALLOCATE(mold(1:ysize), STAT=istat)
  617. ALLOCATE(tmpdata(1:xsize,1:ysize), STAT=istat)
  618. if (OutName .eq. 'ZNU') then
  619. do idx = 1, zsize
  620. half_eta(idx) = Field(1,idx,1,1)
  621. enddo
  622. endif
  623. if (OutName .eq. 'ZNW') then
  624. do idx = 1, zsize
  625. full_eta(idx) = Field(1,idx,1,1)
  626. enddo
  627. endif
  628. if (OutName .eq. 'ZS') then
  629. do idx = 1, zsize
  630. soil_depth(idx) = Field(1,idx,1,1)
  631. enddo
  632. endif
  633. if (OutName .eq. 'DZS') then
  634. do idx = 1, zsize
  635. soil_thickness(idx) = Field(1,idx,1,1)
  636. enddo
  637. endif
  638. if ((xsize .lt. 1) .or. (ysize .lt. 1)) then
  639. write(msg,*) 'Cannot output field with memory order: ', &
  640. MemoryOrder,Varname
  641. call wrf_message(msg)
  642. return
  643. endif
  644. call get_vert_stag(OutName,Stagger,vert_stag)
  645. do idx = 1, zsize
  646. call gr1_get_levels(OutName, idx, zsize, soil_layers, vert_stag, fraction, &
  647. vert_unit, level1(idx), level2(idx))
  648. enddo
  649. !
  650. ! Get the center lat/lon for the area being output. For some cases (such
  651. ! as for boundary areas, the center of the area is different from the
  652. ! center of the model grid.
  653. !
  654. if (index(Stagger,'X') .le. 0) then
  655. dom_xsize = full_xsize - 1
  656. else
  657. dom_xsize = full_xsize
  658. endif
  659. if (index(Stagger,'Y') .le. 0) then
  660. dom_ysize = full_ysize - 1
  661. else
  662. dom_ysize = full_ysize
  663. endif
  664. !
  665. ! Handle case of polare stereographic centered on pole. In that case,
  666. ! always set center lon to be the projection central longitude.
  667. !
  668. if ((projection .eq. WRF_POLAR_STEREO) .AND. &
  669. (abs(center_lat - 90.0) < 0.01)) then
  670. center_lon = proj_central_lon
  671. endif
  672. CALL get_region_center(MemoryOrder, projection, center_lat, center_lon, &
  673. dom_xsize, dom_ysize, dx, dy, proj_central_lon, proj_center_flag, &
  674. truelat1, truelat2, xsize, ysize, region_center_lat, region_center_lon)
  675. if ( .not. opened(DataHandle)) then
  676. Status = WRF_WARN_FILE_NOT_OPENED
  677. return
  678. endif
  679. if (opened(DataHandle) .and. committed(DataHandle)) then
  680. #ifdef OUTPUT_FULL_PRESSURE
  681. !
  682. ! The following is a kludge to output full pressure instead of the two
  683. ! fields of base-state pressure and pressure perturbation.
  684. !
  685. ! This code can be turned on by adding -DOUTPUT_FULL_PRESSURE to the
  686. ! compile line
  687. !
  688. if ((OutName .eq. 'P') .or. (OutName.eq.'PB')) then
  689. do idx = 1, len(MemoryOrder)
  690. if (MemoryOrder(idx:idx) .eq. 'X') then
  691. idx1=idx
  692. endif
  693. if (MemoryOrder(idx:idx) .eq. 'Y') then
  694. idx2=idx
  695. endif
  696. if (MemoryOrder(idx:idx) .eq. 'Z') then
  697. idx3=idx
  698. endif
  699. enddo
  700. !
  701. ! Allocate space for pressure values (this variable holds
  702. ! base-state pressure or pressure perturbation to be used
  703. ! later to sum base-state and perturbation pressure to get full
  704. ! pressure).
  705. !
  706. if (.not. ASSOCIATED(pressure(this_domain)%vals)) then
  707. ALLOCATE(pressure(this_domain)%vals(MemoryStart(1):MemoryEnd(1), &
  708. MemoryStart(2):MemoryEnd(2),MemoryStart(3):MemoryEnd(3)))
  709. endif
  710. if (DateStr .NE. &
  711. pressure(this_domain)%lastDateStr) then
  712. pressure(this_domain)%newtime = .true.
  713. endif
  714. if (pressure(this_domain)%newtime) then
  715. pressure(this_domain)%vals = Field(1,:,:,:)
  716. pressure(this_domain)%newtime = .false.
  717. output = .false.
  718. else
  719. output = .true.
  720. endif
  721. pressure(this_domain)%lastDateStr=DateStr
  722. endif
  723. #endif
  724. #ifdef OUTPUT_FULL_GEOPOTENTIAL
  725. !
  726. ! The following is a kludge to output full geopotential height instead
  727. ! of the two fields of base-state geopotential and perturbation
  728. ! geopotential.
  729. !
  730. ! This code can be turned on by adding -DOUTPUT_FULL_GEOPOTENTIAL to the
  731. ! compile line
  732. !
  733. if ((OutName .eq. 'PHB') .or. (OutName.eq.'PH')) then
  734. do idx = 1, len(MemoryOrder)
  735. if (MemoryOrder(idx:idx) .eq. 'X') then
  736. idx1=idx
  737. endif
  738. if (MemoryOrder(idx:idx) .eq. 'Y') then
  739. idx2=idx
  740. endif
  741. if (MemoryOrder(idx:idx) .eq. 'Z') then
  742. idx3=idx
  743. endif
  744. enddo
  745. !
  746. ! Allocate space for geopotential values (this variable holds
  747. ! geopotential to be used
  748. ! later to sum base-state and perturbation to get full
  749. ! geopotential).
  750. !
  751. if (.not. ASSOCIATED(geopotential(this_domain)%vals)) then
  752. ALLOCATE(geopotential(this_domain)%vals(MemoryStart(1):MemoryEnd(1), &
  753. MemoryStart(2):MemoryEnd(2),MemoryStart(3):MemoryEnd(3)))
  754. endif
  755. if (DateStr .NE. &
  756. geopotential(this_domain)%lastDateStr) then
  757. geopotential(this_domain)%newtime = .true.
  758. endif
  759. if (geopotential(this_domain)%newtime) then
  760. geopotential(this_domain)%vals = Field(1,:,:,:)
  761. geopotential(this_domain)%newtime = .false.
  762. output = .false.
  763. else
  764. output = .true.
  765. endif
  766. geopotential(this_domain)%lastDateStr=DateStr
  767. endif
  768. #endif
  769. if (output) then
  770. if (StartDate == '') then
  771. StartDate = DateStr
  772. endif
  773. CALL geth_idts(DateStr,StartDate,fcst_secs)
  774. if (center_lat .lt. 0) then
  775. proj_center_flag = 2
  776. else
  777. proj_center_flag = 1
  778. endif
  779. do z = 1, zsize
  780. SELECT CASE (MemoryOrder)
  781. CASE ('XYZ')
  782. data = Field(1,1:xsize,1:ysize,z)
  783. CASE ('XZY')
  784. data = Field(1,1:xsize,z,1:ysize)
  785. CASE ('YXZ')
  786. do x = 1,xsize
  787. do y = 1,ysize
  788. data(x,y) = Field(1,y,x,z)
  789. enddo
  790. enddo
  791. CASE ('YZX')
  792. do x = 1,xsize
  793. do y = 1,ysize
  794. data(x,y) = Field(1,y,z,x)
  795. enddo
  796. enddo
  797. CASE ('ZXY')
  798. data = Field(1,z,1:xsize,1:ysize)
  799. CASE ('ZYX')
  800. do x = 1,xsize
  801. do y = 1,ysize
  802. data(x,y) = Field(1,z,y,x)
  803. enddo
  804. enddo
  805. CASE ('XY')
  806. data = Field(1,1:xsize,1:ysize,1)
  807. CASE ('YX')
  808. do x = 1,xsize
  809. do y = 1,ysize
  810. data(x,y) = Field(1,y,x,1)
  811. enddo
  812. enddo
  813. CASE ('XSZ')
  814. do x = 1,xsize
  815. do y = 1,ysize
  816. data(x,y) = Field(1,y,z,x)
  817. enddo
  818. enddo
  819. CASE ('XEZ')
  820. do x = 1,xsize
  821. do y = 1,ysize
  822. data(x,y) = Field(1,y,z,x)
  823. enddo
  824. enddo
  825. CASE ('YSZ')
  826. do x = 1,xsize
  827. do y = 1,ysize
  828. data(x,y) = Field(1,x,z,y)
  829. enddo
  830. enddo
  831. CASE ('YEZ')
  832. do x = 1,xsize
  833. do y = 1,ysize
  834. data(x,y) = Field(1,x,z,y)
  835. enddo
  836. enddo
  837. CASE ('XS')
  838. do x = 1,xsize
  839. do y = 1,ysize
  840. data(x,y) = Field(1,y,x,1)
  841. enddo
  842. enddo
  843. CASE ('XE')
  844. do x = 1,xsize
  845. do y = 1,ysize
  846. data(x,y) = Field(1,y,x,1)
  847. enddo
  848. enddo
  849. CASE ('YS')
  850. do x = 1,xsize
  851. do y = 1,ysize
  852. data(x,y) = Field(1,x,y,1)
  853. enddo
  854. enddo
  855. CASE ('YE')
  856. do x = 1,xsize
  857. do y = 1,ysize
  858. data(x,y) = Field(1,x,y,1)
  859. enddo
  860. enddo
  861. CASE ('Z')
  862. data(1,1) = Field(1,z,1,1)
  863. CASE ('z')
  864. data(1,1) = Field(1,z,1,1)
  865. CASE ('C')
  866. data = Field(1,1:xsize,1:ysize,z)
  867. CASE ('c')
  868. data = Field(1,1:xsize,1:ysize,z)
  869. CASE ('0')
  870. data(1,1) = Field(1,1,1,1)
  871. END SELECT
  872. !
  873. ! Here, we convert any integer fields to real
  874. !
  875. if (FieldType == WRF_INTEGER) then
  876. mold = 0
  877. do idx=1,xsize
  878. !
  879. ! The parentheses around data(idx,:) are needed in order
  880. ! to fix a bug with transfer with the xlf compiler on NCAR's
  881. ! IBM (bluesky).
  882. !
  883. data(idx,:)=transfer((data(idx,:)),mold)
  884. enddo
  885. endif
  886. !
  887. ! Here, we do any necessary conversions to the data.
  888. !
  889. ! Potential temperature is sometimes passed in as perturbation
  890. ! potential temperature (i.e., POT-300). Other times (i.e., from
  891. ! WRF SI), it is passed in as full potential temperature.
  892. ! Here, we convert to full potential temperature by adding 300
  893. ! only if POT < 200 K.
  894. !
  895. if (OutName == 'T') then
  896. if (data(1,1) < 200) then
  897. data = data + 300
  898. endif
  899. endif
  900. !
  901. ! For precip, we setup the accumulation period, and output a precip
  902. ! rate for time-step precip.
  903. !
  904. if (OutName .eq. 'RAINNCV') then
  905. ! Convert time-step precip to precip rate.
  906. data = data/timestep
  907. accum_period = 0
  908. else
  909. accum_period = 0
  910. endif
  911. #ifdef OUTPUT_FULL_PRESSURE
  912. !
  913. ! Computation of full-pressure off by default since there are
  914. ! uses for base-state and perturbation (i.e., restarts
  915. !
  916. if ((OutName .eq. 'P') .or. (OutName.eq.'PB')) then
  917. if (idx3 .eq. 1) then
  918. data = data + &
  919. pressure(this_domain)%vals(z, &
  920. patchstart(2):patchend(2),patchstart(3):patchend(3))
  921. elseif (idx3 .eq. 2) then
  922. data = data + &
  923. pressure(this_domain)%vals(patchstart(1):patchend(1), &
  924. z,patchstart(3):patchend(3))
  925. elseif (idx3 .eq. 3) then
  926. data = data + &
  927. pressure(this_domain)%vals(patchstart(1):patchend(1), &
  928. patchstart(2):patchend(2),z)
  929. else
  930. call wrf_message ('error in idx3, continuing')
  931. endif
  932. OutName = 'P'
  933. endif
  934. #endif
  935. #ifdef OUTPUT_FULL_GEOPOTENTIAL
  936. !
  937. ! Computation of full-geopotential off by default since there are
  938. ! uses for base-state and perturbation (i.e., restarts
  939. !
  940. if ((OutName .eq. 'PHB') .or. (OutName.eq.'PH')) then
  941. if (idx3 .eq. 1) then
  942. data = data + &
  943. geopotential(this_domain)%vals(z, &
  944. patchstart(2):patchend(2),patchstart(3):patchend(3))
  945. elseif (idx3 .eq. 2) then
  946. data = data + &
  947. geopotential(this_domain)%vals(patchstart(1):patchend(1), &
  948. z,patchstart(3):patchend(3))
  949. elseif (idx3 .eq. 3) then
  950. data = data + &
  951. geopotential(this_domain)%vals(patchstart(1):patchend(1), &
  952. patchstart(2):patchend(2),z)
  953. else
  954. call wrf_message ('error in idx3, continuing')
  955. endif
  956. OutName = 'PHP'
  957. endif
  958. #endif
  959. !
  960. ! Output current level
  961. !
  962. CALL load_grid_info(OutName, StartDate, vert_unit, level1(z), &
  963. level2(z), fcst_secs, accum_period, wg_grid_id, projection, &
  964. xsize, ysize, region_center_lat, region_center_lon, dx, dy, &
  965. proj_central_lon, proj_center_flag, truelat1, truelat2, &
  966. grib_tables, grid_info)
  967. !
  968. ! Here, we copy data to a temporary array. After write_grib,
  969. ! we copy back from the temporary array to the permanent
  970. ! array. write_grib modifies data. For certain fields that
  971. ! we use below, we want the original (unmodified) data
  972. ! values. This kludge assures that we have the original
  973. ! values.
  974. !
  975. if ((OutName .eq. 'RAINC') .or. (OutName .eq. 'RAINNC') .or. &
  976. (OutName .eq. 'SNOWH')) then
  977. tmpdata(:,:) = data(:,:)
  978. endif
  979. CALL write_grib(grid_info, FileFd(DataHandle), data)
  980. if ((OutName .eq. 'RAINC') .or. (OutName .eq. 'RAINNC') .or. &
  981. (OutName .eq. 'SNOWH')) then
  982. data(:,:) = tmpdata(:,:)
  983. endif
  984. CALL free_grid_info(grid_info)
  985. !
  986. ! If this is the total accumulated rain, call write_grib again
  987. ! to output the accumulation since the last output time as well.
  988. ! This is somewhat of a kludge to meet the requirements of PF.
  989. !
  990. if ((OutName .eq. 'RAINC') .or. (OutName .eq. 'RAINNC') .or. &
  991. (OutName .eq. 'SNOWH')) then
  992. if (OutName .eq. 'RAINC') then
  993. data(:,:) = data(:,:) - lastdata(this_domain)%rainc(:,:)
  994. lastdata(this_domain)%rainc(:,:) = tmpdata(:,:)
  995. accum_period = fcst_secs - &
  996. lastdata(this_domain)%fcst_secs_rainc
  997. lastdata(this_domain)%fcst_secs_rainc = fcst_secs
  998. TmpVarName = 'ACPCP'
  999. else if (OutName .eq. 'RAINNC') then
  1000. tmpdata(:,:) = data(:,:)
  1001. data(:,:) = data(:,:) - lastdata(this_domain)%rainnc(:,:)
  1002. lastdata(this_domain)%rainnc(:,:) = tmpdata(:,:)
  1003. accum_period = fcst_secs - &
  1004. lastdata(this_domain)%fcst_secs_rainnc
  1005. lastdata(this_domain)%fcst_secs_rainnc = fcst_secs
  1006. TmpVarName = 'NCPCP'
  1007. else if (OutName .eq. 'SNOWH') then
  1008. if (fcst_secs .eq. 0) then
  1009. firstdata(this_domain)%snod(:,:) = data(:,:)
  1010. endif
  1011. data(:,:) = data(:,:) - firstdata(this_domain)%snod(:,:)
  1012. TmpVarName = 'SNOWCU'
  1013. endif
  1014. CALL load_grid_info(TmpVarName, StartDate, vert_unit, level1(z),&
  1015. level2(z), fcst_secs, accum_period, wg_grid_id, &
  1016. projection, xsize, ysize, region_center_lat, &
  1017. region_center_lon, dx, dy, proj_central_lon, &
  1018. proj_center_flag, truelat1, truelat2, grib_tables, &
  1019. grid_info)
  1020. CALL write_grib(grid_info, FileFd(DataHandle), data)
  1021. CALL free_grid_info(grid_info)
  1022. endif
  1023. enddo
  1024. endif
  1025. endif
  1026. deallocate(data, STAT = istat)
  1027. deallocate(mold, STAT = istat)
  1028. deallocate(tmpdata, STAT = istat)
  1029. Status = WRF_NO_ERR
  1030. call wrf_debug ( DEBUG , 'Leaving ext_gr1_write_field')
  1031. RETURN
  1032. END SUBROUTINE ext_gr1_write_field
  1033. !*****************************************************************************
  1034. SUBROUTINE ext_gr1_read_field ( DataHandle , DateStr , VarName , Field , &
  1035. FieldType , Comm , IOComm, DomainDesc , MemoryOrder , Stagger , &
  1036. DimNames , DomainStart , DomainEnd , MemoryStart , MemoryEnd , &
  1037. PatchStart , PatchEnd , Status )
  1038. USE gr1_data_info
  1039. IMPLICIT NONE
  1040. #include "wrf_status_codes.h"
  1041. #include "wrf_io_flags.h"
  1042. INTEGER , INTENT(IN) :: DataHandle
  1043. CHARACTER*(*) :: DateStr
  1044. CHARACTER*(*) :: VarName
  1045. CHARACTER (len=400) :: msg
  1046. integer ,intent(inout) :: FieldType
  1047. integer ,intent(inout) :: Comm
  1048. integer ,intent(inout) :: IOComm
  1049. integer ,intent(inout) :: DomainDesc
  1050. character*(*) ,intent(inout) :: MemoryOrder
  1051. character*(*) ,intent(inout) :: Stagger
  1052. character*(*) , dimension (*) ,intent(inout) :: DimNames
  1053. integer ,dimension(*) ,intent(inout) :: DomainStart, DomainEnd
  1054. integer ,dimension(*) ,intent(inout) :: MemoryStart, MemoryEnd
  1055. integer ,dimension(*) ,intent(inout) :: PatchStart, PatchEnd
  1056. integer ,intent(out) :: Status
  1057. INTEGER ,intent(out) :: Field(*)
  1058. integer :: ndim,x_start,x_end,y_start,y_end,z_start,z_end
  1059. integer :: zidx
  1060. REAL, DIMENSION(:,:), POINTER :: data
  1061. logical :: vert_stag
  1062. logical :: soil_layers
  1063. integer :: level1,level2
  1064. integer :: parmid
  1065. integer :: vert_unit
  1066. integer :: grb_index
  1067. integer :: numcols, numrows
  1068. integer :: data_allocated
  1069. integer :: istat
  1070. integer :: tablenum
  1071. integer :: di
  1072. integer :: last_grb_index
  1073. call wrf_debug ( DEBUG , 'Entering ext_gr1_read_field')
  1074. !
  1075. ! Get dimensions of data.
  1076. ! Assume that the domain size in the input data is the same as the Domain
  1077. ! Size from the input arguments.
  1078. !
  1079. CALL get_dims(MemoryOrder,DomainStart,DomainEnd,ndim,x_start,x_end,y_start, &
  1080. y_end,z_start,z_end)
  1081. !
  1082. ! Get grib parameter id
  1083. !
  1084. CALL GET_GRIB_PARAM(grib_tables, VarName, center, subcenter, parmtbl, &
  1085. tablenum, parmid)
  1086. !
  1087. ! Setup the vertical unit and levels
  1088. !
  1089. CALL get_vert_stag(VarName,Stagger,vert_stag)
  1090. CALL get_soil_layers(VarName,soil_layers)
  1091. !
  1092. ! Loop over levels, grabbing data from each level, then assembling into a
  1093. ! 3D array.
  1094. !
  1095. data_allocated = 0
  1096. last_grb_index = -1
  1097. do zidx = z_start,z_end
  1098. CALL gr1_get_levels(VarName,zidx,z_end-z_start,soil_layers,vert_stag, &
  1099. .false., vert_unit,level1,level2)
  1100. CALL GET_GRIB_INDEX_VALIDTIME_GUESS(fileinfo(DataHandle)%fileindex(:), center, &
  1101. subcenter, parmtbl, parmid,DateStr,vert_unit,level1, &
  1102. level2, last_grb_index + 1, grb_index)
  1103. if (grb_index < 0) then
  1104. write(msg,*)'Field not found: parmid: ',VarName,parmid,DateStr, &
  1105. vert_unit,level1,level2
  1106. call wrf_debug (DEBUG , msg)
  1107. cycle
  1108. endif
  1109. if (data_allocated .eq. 0) then
  1110. CALL GET_SIZEOF_GRID(fileinfo(DataHandle)%fileindex(:),grb_index,numcols,numrows)
  1111. allocate(data(z_start:z_end,1:numcols*numrows),stat=istat)
  1112. data_allocated = 1
  1113. endif
  1114. CALL READ_GRIB(fileinfo(DataHandle)%fileindex(:), FileFd(DataHandle), grb_index, &
  1115. data(zidx,:))
  1116. !
  1117. ! Transpose data into the order specified by MemoryOrder, setting only
  1118. ! entries within the memory dimensions
  1119. !
  1120. CALL get_dims(MemoryOrder, MemoryStart, MemoryEnd, ndim, x_start, x_end, &
  1121. y_start, y_end,z_start,z_end)
  1122. if(FieldType == WRF_DOUBLE) then
  1123. di = 2
  1124. else
  1125. di = 1
  1126. endif
  1127. !
  1128. ! Here, we do any necessary conversions to the data.
  1129. !
  1130. ! The WRF executable (wrf.exe) expects perturbation potential
  1131. ! temperature. However, real.exe expects full potential T.
  1132. ! So, if the program is WRF, subtract 300 from Potential Temperature
  1133. ! to get perturbation potential temperature.
  1134. !
  1135. if (VarName == 'T') then
  1136. if ( &
  1137. (InputProgramName .eq. 'REAL_EM') .or. &
  1138. (InputProgramName .eq. 'IDEAL') .or. &
  1139. (InputProgramName .eq. 'NDOWN_EM')) then
  1140. data(zidx,:) = data(zidx,:) - 300
  1141. endif
  1142. endif
  1143. CALL Transpose_grib(MemoryOrder, di, FieldType, Field, &
  1144. MemoryStart(1), MemoryEnd(1), MemoryStart(2), MemoryEnd(2), &
  1145. MemoryStart(3), MemoryEnd(3), &
  1146. data(zidx,:), zidx, numrows, numcols)
  1147. if (zidx .eq. z_end) then
  1148. data_allocated = 0
  1149. deallocate(data)
  1150. endif
  1151. last_grb_index = grb_index
  1152. enddo
  1153. Status = WRF_NO_ERR
  1154. if (grb_index < 0) Status = WRF_WARN_VAR_NF
  1155. call wrf_debug ( DEBUG , 'Leaving ext_gr1_read_field')
  1156. RETURN
  1157. END SUBROUTINE ext_gr1_read_field
  1158. !*****************************************************************************
  1159. SUBROUTINE ext_gr1_get_next_var ( DataHandle, VarName, Status )
  1160. USE gr1_data_info
  1161. IMPLICIT NONE
  1162. #include "wrf_status_codes.h"
  1163. INTEGER , INTENT(IN) :: DataHandle
  1164. CHARACTER*(*) :: VarName
  1165. INTEGER , INTENT(OUT) :: Status
  1166. call wrf_debug ( DEBUG , 'Entering ext_gr1_get_next_var')
  1167. call wrf_message ( 'WARNING: ext_gr1_get_next_var is not supported.')
  1168. Status = WRF_WARN_NOOP
  1169. RETURN
  1170. END SUBROUTINE ext_gr1_get_next_var
  1171. !*****************************************************************************
  1172. subroutine ext_gr1_end_of_frame(DataHandle, Status)
  1173. USE gr1_data_info
  1174. implicit none
  1175. #include "wrf_status_codes.h"
  1176. integer ,intent(in) :: DataHandle
  1177. integer ,intent(out) :: Status
  1178. call wrf_debug ( DEBUG , 'Entering ext_gr1_end_of_frame')
  1179. Status = WRF_WARN_NOOP
  1180. return
  1181. end subroutine ext_gr1_end_of_frame
  1182. !*****************************************************************************
  1183. SUBROUTINE ext_gr1_iosync ( DataHandle, Status )
  1184. USE gr1_data_info
  1185. IMPLICIT NONE
  1186. #include "wrf_status_codes.h"
  1187. INTEGER , INTENT(IN) :: DataHandle
  1188. INTEGER , INTENT(OUT) :: Status
  1189. call wrf_debug ( DEBUG , 'Entering ext_gr1_iosync')
  1190. Status = WRF_NO_ERR
  1191. if (DataHandle .GT. 0) then
  1192. CALL flush_file(FileFd(DataHandle))
  1193. else
  1194. Status = WRF_WARN_TOO_MANY_FILES
  1195. endif
  1196. RETURN
  1197. END SUBROUTINE ext_gr1_iosync
  1198. !*****************************************************************************
  1199. SUBROUTINE ext_gr1_inquire_filename ( DataHandle, FileName , FileStat, &
  1200. Status )
  1201. USE gr1_data_info
  1202. IMPLICIT NONE
  1203. #include "wrf_status_codes.h"
  1204. #include "wrf_io_flags.h"
  1205. INTEGER , INTENT(IN) :: DataHandle
  1206. CHARACTER*(*) :: FileName
  1207. INTEGER , INTENT(OUT) :: FileStat
  1208. INTEGER , INTENT(OUT) :: Status
  1209. CHARACTER *80 SysDepInfo
  1210. call wrf_debug ( DEBUG , 'Entering ext_gr1_inquire_filename')
  1211. FileName = DataFile(DataHandle)
  1212. if ((DataHandle .ge. firstFileHandle) .and. &
  1213. (DataHandle .le. maxFileHandles)) then
  1214. FileStat = FileStatus(DataHandle)
  1215. else
  1216. FileStat = WRF_FILE_NOT_OPENED
  1217. endif
  1218. Status = WRF_NO_ERR
  1219. RETURN
  1220. END SUBROUTINE ext_gr1_inquire_filename
  1221. !*****************************************************************************
  1222. SUBROUTINE ext_gr1_get_var_info ( DataHandle , VarName , NDim , &
  1223. MemoryOrder , Stagger , DomainStart , DomainEnd , WrfType, Status )
  1224. USE gr1_data_info
  1225. IMPLICIT NONE
  1226. #include "wrf_status_codes.h"
  1227. integer ,intent(in) :: DataHandle
  1228. character*(*) ,intent(in) :: VarName
  1229. integer ,intent(out) :: NDim
  1230. character*(*) ,intent(out) :: MemoryOrder
  1231. character*(*) ,intent(out) :: Stagger
  1232. integer ,dimension(*) ,intent(out) :: DomainStart, DomainEnd
  1233. integer ,intent(out) :: WrfType
  1234. integer ,intent(out) :: Status
  1235. call wrf_debug ( DEBUG , 'Entering ext_gr1_get_var_info')
  1236. CALL wrf_message('ext_gr1_get_var_info not supported for grib version1 data')
  1237. Status = WRF_NO_ERR
  1238. RETURN
  1239. END SUBROUTINE ext_gr1_get_var_info
  1240. !*****************************************************************************
  1241. SUBROUTINE ext_gr1_set_time ( DataHandle, DateStr, Status )
  1242. USE gr1_data_info
  1243. IMPLICIT NONE
  1244. #include "wrf_status_codes.h"
  1245. INTEGER , INTENT(IN) :: DataHandle
  1246. CHARACTER*(*) :: DateStr
  1247. INTEGER , INTENT(OUT) :: Status
  1248. integer :: found_time
  1249. integer :: idx
  1250. call wrf_debug ( DEBUG , 'Entering ext_gr1_set_time')
  1251. found_time = 0
  1252. do idx = 1,fileinfo(DataHandle)%NumberTimes
  1253. if (fileinfo(DataHandle)%Times(idx) == DateStr) then
  1254. found_time = 1
  1255. fileinfo(DataHandle)%CurrentTime = idx
  1256. endif
  1257. enddo
  1258. if (found_time == 0) then
  1259. Status = WRF_WARN_TIME_NF
  1260. else
  1261. Status = WRF_NO_ERR
  1262. endif
  1263. RETURN
  1264. END SUBROUTINE ext_gr1_set_time
  1265. !*****************************************************************************
  1266. SUBROUTINE ext_gr1_get_next_time ( DataHandle, DateStr, Status )
  1267. USE gr1_data_info
  1268. IMPLICIT NONE
  1269. #include "wrf_status_codes.h"
  1270. INTEGER , INTENT(IN) :: DataHandle
  1271. CHARACTER*(*) , INTENT(OUT) :: DateStr
  1272. INTEGER , INTENT(OUT) :: Status
  1273. call wrf_debug ( DEBUG , 'Entering ext_gr1_get_next_time')
  1274. if (fileinfo(DataHandle)%CurrentTime == fileinfo(DataHandle)%NumberTimes) then
  1275. Status = WRF_WARN_TIME_EOF
  1276. else
  1277. fileinfo(DataHandle)%CurrentTime = fileinfo(DataHandle)%CurrentTime + 1
  1278. DateStr = fileinfo(DataHandle)%Times(fileinfo(DataHandle)%CurrentTime)
  1279. Status = WRF_NO_ERR
  1280. endif
  1281. RETURN
  1282. END SUBROUTINE ext_gr1_get_next_time
  1283. !*****************************************************************************
  1284. SUBROUTINE ext_gr1_get_previous_time ( DataHandle, DateStr, Status )
  1285. USE gr1_data_info
  1286. IMPLICIT NONE
  1287. #include "wrf_status_codes.h"
  1288. INTEGER , INTENT(IN) :: DataHandle
  1289. CHARACTER*(*) :: DateStr
  1290. INTEGER , INTENT(OUT) :: Status
  1291. call wrf_debug ( DEBUG , 'Entering ext_gr1_get_previous_time')
  1292. if (fileinfo(DataHandle)%CurrentTime <= 0) then
  1293. Status = WRF_WARN_TIME_EOF
  1294. else
  1295. fileinfo(DataHandle)%CurrentTime = fileinfo(DataHandle)%CurrentTime - 1
  1296. DateStr = fileinfo(DataHandle)%Times(fileinfo(DataHandle)%CurrentTime)
  1297. Status = WRF_NO_ERR
  1298. endif
  1299. RETURN
  1300. END SUBROUTINE ext_gr1_get_previous_time
  1301. !******************************************************************************
  1302. !* Start of get_var_ti_* routines
  1303. !*************************************************

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