PageRenderTime 51ms CodeModel.GetById 14ms RepoModel.GetById 0ms app.codeStats 1ms

/WPS/ungrib/src/rrpr.F

http://github.com/jbeezley/wrf-fire
FORTRAN Legacy | 887 lines | 476 code | 124 blank | 287 comment | 73 complexity | 6aa7e3b2ae5c1204faa3dfffdcbc4686 MD5 | raw file
Possible License(s): AGPL-1.0
  1. subroutine rrpr(hstart, ntimes, interval, nlvl, maxlvl, plvl, debug_level, out_format, prefix)
  2. ! !
  3. ! In case you are wondering, RRPR stands for "Read, ReProcess, and wRite" !
  4. ! !
  5. !*****************************************************************************!
  6. ! !
  7. use filelist
  8. use gridinfo
  9. use storage_module
  10. use table
  11. use module_debug
  12. use misc_definitions_module
  13. use stringutil
  14. implicit none
  15. !------------------------------------------------------------------------------
  16. ! Arguments:
  17. ! HSTART: Starting date of times to process
  18. character (LEN=19) :: hstart
  19. ! NTIMES: Number of time periods to process
  20. integer :: ntimes
  21. ! INTERVAL: Time inteval (seconds) of time periods to process.
  22. integer :: interval
  23. ! NLVL: The number of levels in the stored data.
  24. integer :: nlvl
  25. ! MAXLVL: The parameterized maximum number of levels to allow.
  26. integer :: maxlvl
  27. ! PLVL: Array of pressure levels (Pa) in the dataset
  28. real , dimension(maxlvl) :: plvl
  29. ! DEBUG_LEVEL: Integer level of debug printing (from namelist)
  30. integer :: debug_level
  31. !------------------------------------------------------------------------------
  32. character (LEN=25) :: units
  33. character (LEN=46) :: Desc
  34. real, allocatable, dimension(:,:) :: scr2d, tmp2d
  35. real, pointer, dimension(:,:) :: ptr2d
  36. integer :: k, kk, mm, n, ierr, ifv
  37. integer :: iunit=13
  38. character(LEN=19) :: hdate, hend
  39. character(LEN=24) :: hdate_output
  40. character(LEN=3) :: out_format
  41. character(LEN=MAX_FILENAME_LEN) :: prefix
  42. real :: xfcst, level
  43. character(LEN=9) :: field
  44. integer :: ntime, idts
  45. ! DATELEN: length of date strings to use for our output file names.
  46. integer :: datelen
  47. ! Decide the length of date strings to use for output file names.
  48. ! DATELEN is 13 for hours, 16 for minutes, and 19 for seconds.
  49. if (mod(interval,3600) == 0) then
  50. datelen = 13
  51. else if (mod(interval, 60) == 0) then
  52. datelen = 16
  53. else
  54. datelen = 19
  55. endif
  56. if ( debug_level .gt. 100 ) then
  57. call mprintf(.true.,DEBUG,"Begin rrpr")
  58. call mprintf(.true.,DEBUG,"nfiles = %i , ntimes = %i )",i1=nfiles,i2=ntimes)
  59. do n = 1, nfiles
  60. call mprintf(.true.,DEBUG,"filedates(%i) = %s",i1=n,s1=filedates(n))
  61. enddo
  62. endif
  63. ! Compute the ending time:
  64. call geth_newdate(hend, hstart, interval*ntimes)
  65. call clear_storage
  66. ! We want to do something for each of the requested times:
  67. TIMELOOP : do ntime = 1, ntimes
  68. idts = (ntime-1) * interval
  69. call geth_newdate(hdate, hstart, idts)
  70. call mprintf(.true.,DEBUG, &
  71. "RRPR: hstart = %s , hdate = %s , idts = %i",s1=hstart,s2=hdate,i1=idts)
  72. ! Loop over the output file dates, and do stuff if the file date matches
  73. ! the requested time we are working on now.
  74. FILELOOP : do n = 1, nfiles
  75. if ( debug_level .gt. 100 ) then
  76. call mprintf(.true.,DEBUG, &
  77. "hstart = %s , hend = %s",s1=hstart,s2=hend)
  78. call mprintf(.true.,DEBUG, &
  79. "filedates(n) = %s",s1=filedates(n))
  80. call mprintf(.true.,DEBUG, &
  81. "filedates(n) = %s",s1=filedates(n)(1:datelen))
  82. end if
  83. if (filedates(n)(1:datelen).ne.hdate(1:datelen)) cycle FILELOOP
  84. if (debug_level .gt. 50 ) then
  85. call mprintf(.true.,INFORM, &
  86. "RRPR Processing : %s",s1=filedates(n)(1:datelen))
  87. endif
  88. open(iunit, file=trim(get_path(prefix))//'PFILE:'//filedates(n)(1:datelen), &
  89. form='unformatted',status='old')
  90. ! Read the file:
  91. rdloop: do
  92. read (iunit, iostat=ierr) ifv
  93. if (ierr.ne.0) exit rdloop
  94. if ( ifv .eq. 5) then ! WPS
  95. read (iunit) hdate_output, xfcst, map%source, field, units, Desc, &
  96. level, map%nx, map%ny, map%igrid
  97. hdate = hdate_output(1:19)
  98. select case (map%igrid)
  99. case (0, 4)
  100. read (iunit) map%startloc, map%lat1, map%lon1, map%dy, map%dx, map%r_earth
  101. case (3)
  102. read (iunit) map%startloc, map%lat1, map%lon1, map%dx, map%dy, map%lov, &
  103. map%truelat1, map%truelat2, map%r_earth
  104. case (5)
  105. read (iunit) map%startloc, map%lat1, map%lon1, map%dx, map%dy, map%lov, &
  106. map%truelat1, map%r_earth
  107. case (1)
  108. read (iunit) map%startloc, map%lat1, map%lon1, map%dy, map%dx, &
  109. map%truelat1, map%r_earth
  110. case default
  111. call mprintf(.true.,ERROR, &
  112. "Unrecognized map%%igrid: %i in RRPR 1",i1=map%igrid)
  113. end select
  114. read (iunit) map%grid_wind
  115. else if ( ifv .eq. 4 ) then ! SI
  116. read (iunit) hdate_output, xfcst, map%source, field, units, desc, level, &
  117. map%nx, map%ny, map%igrid
  118. hdate = hdate_output(1:19)
  119. select case (map%igrid)
  120. case (0, 4)
  121. read(iunit) map%startloc, map%lat1, map%lon1, map%dy, map%dx
  122. case (3)
  123. read (iunit) map%startloc, map%lat1, map%lon1, map%dx, map%dy, &
  124. map%lov, map%truelat1, map%truelat2
  125. case (5)
  126. read (iunit) map%startloc, map%lat1, map%lon1, map%dx, map%dy, &
  127. map%lov, map%truelat1
  128. case default
  129. call mprintf(.true.,ERROR, &
  130. "Unrecognized map%%igrid: %i in RRPR 2",i1=map%igrid)
  131. end select
  132. else if ( ifv .eq. 3 ) then ! MM5
  133. read(iunit) hdate_output, xfcst, field, units, desc, level,&
  134. map%nx, map%ny, map%igrid
  135. hdate = hdate_output(1:19)
  136. select case (map%igrid)
  137. case (3) ! lamcon
  138. read (iunit) map%lat1, map%lon1, map%dx, map%dy, map%lov, &
  139. map%truelat1, map%truelat2
  140. case (5) ! Polar Stereographic
  141. read (iunit) map%lat1, map%lon1, map%dx, map%dy, map%lov, &
  142. map%truelat1
  143. case (0, 4) ! lat/lon
  144. read (iunit) map%lat1, map%lon1, map%dy, map%dx
  145. case (1) ! Mercator
  146. read (iunit) map%lat1, map%lon1, map%dy, map%dx, map%truelat1
  147. case default
  148. call mprintf(.true.,ERROR, &
  149. "Unrecognized map%%igrid: %i in RRPR 3",i1=map%igrid)
  150. end select
  151. else
  152. call mprintf(.true.,ERROR, &
  153. "unknown out_format, ifv = %i",i1=ifv)
  154. endif
  155. allocate(ptr2d(map%nx,map%ny))
  156. read (iunit) ptr2d
  157. call refw_storage(nint(level), field, ptr2d, map%nx, map%ny)
  158. nullify (ptr2d)
  159. enddo rdloop
  160. !
  161. ! We have reached the end of file, so time to close it.
  162. !
  163. close(iunit)
  164. if (debug_level .gt. 100 ) call print_storage
  165. !
  166. ! By now the file has been read completely. Now, see if we need to fill in
  167. ! missing fields:
  168. !
  169. ! Retrieve the number of levels in storage:
  170. !
  171. call get_plvls(plvl, maxlvl, nlvl)
  172. !
  173. ! Fill the surface level (code 200100) from higher 200100s, as necessary
  174. !
  175. do k = 1, nlvl
  176. if ((plvl(k).gt.200100) .and. (plvl(k).lt.200200)) then
  177. ! We found a level between 200100 and 200200, now find the field
  178. ! corresponding to that level.
  179. MLOOP : do mm = 1, maxvar
  180. if (is_there(nint(plvl(k)), namvar(mm))) then
  181. INLOOP : do kk = 200101, nint(plvl(k))
  182. if (is_there(kk, namvar(mm))) then
  183. if ( debug_level .gt. 100 ) then
  184. call mprintf(.true.,DEBUG, &
  185. "Copying %s at level %i to level 200100.",s1=namvar(mm),i1=kk)
  186. end if
  187. call get_dims(kk, namvar(mm))
  188. allocate(scr2d(map%nx,map%ny))
  189. call get_storage &
  190. (kk, namvar(mm), scr2d, map%nx, map%ny)
  191. call put_storage &
  192. (200100,namvar(mm), scr2d,map%nx,map%ny)
  193. deallocate(scr2d)
  194. EXIT INLOOP
  195. endif
  196. enddo INLOOP
  197. endif
  198. enddo MLOOP
  199. endif
  200. enddo
  201. !
  202. ! If upper-air U is missing, see if we can interpolate from surrounding levels.
  203. ! This is a simple vertical interpolation, linear in pressure.
  204. ! Currently, this simply fills in one missing level between two present levels.
  205. !
  206. do k = 2, nlvl-1, 1
  207. if (plvl(k-1) .lt. 200000.) then
  208. if ( (.not. is_there(nint(plvl(k)),'UU')) .and. &
  209. ( is_there(nint(plvl(k-1)), 'UU')) .and.&
  210. ( is_there(nint(plvl(k+1)), 'UU')) ) then
  211. call get_dims(nint(plvl(k+1)), 'UU')
  212. call vntrp(plvl, maxlvl, k, "UU ", map%nx, map%ny)
  213. endif
  214. endif
  215. enddo
  216. !
  217. ! If upper-air V is missing, see if we can interpolate from surrounding levels.
  218. ! This is a simple vertical interpolation, linear in pressure.
  219. ! Currently, this simply fills in one missing level between two present levels.
  220. !
  221. do k = 2, nlvl-1, 1
  222. if (plvl(k-1) .lt. 200000.) then
  223. if ( (.not. is_there(nint(plvl(k)),'VV')) .and. &
  224. ( is_there(nint(plvl(k-1)), 'VV')) .and.&
  225. ( is_there(nint(plvl(k+1)), 'VV')) ) then
  226. call get_dims(nint(plvl(k+1)), 'VV')
  227. call vntrp(plvl, maxlvl, k, "VV ", map%nx, map%ny)
  228. endif
  229. endif
  230. enddo
  231. !
  232. ! If upper-air SPECHUMD is missing, see if we can compute SPECHUMD from QVAPOR:
  233. !--- Tanya's change for initializing WRF with RUC
  234. do k = 1, nlvl
  235. if (plvl(k).lt.200000.) then
  236. if (.not. is_there(nint(plvl(k)), 'SPECHUMD').and. &
  237. is_there(nint(plvl(k)), 'QV')) then
  238. call get_dims(nint(plvl(k)), 'QV')
  239. call compute_spechumd_qvapor(map%nx, map%ny, plvl(k))
  240. endif
  241. endif
  242. enddo
  243. !--- Tanya's change for initializing WRF with RUC
  244. ! This allows for the ingestion for RUC isentropic data
  245. !
  246. do k = 1, nlvl
  247. if (plvl(k).lt.200000.) then
  248. if (.not. is_there(nint(plvl(k)), 'TT').and. &
  249. is_there(nint(plvl(k)), 'VPTMP').and. &
  250. is_there(nint(plvl(k)), 'SPECHUMD')) then
  251. call get_dims(nint(plvl(k)), 'VPTMP')
  252. call compute_t_vptmp(map%nx, map%ny, plvl(k))
  253. endif
  254. endif
  255. enddo
  256. !!!
  257. !
  258. ! If upper-air T is missing, see if we can interpolate from surrounding levels.
  259. ! This is a simple vertical interpolation, linear in pressure.
  260. ! Currently, this simply fills in one missing level between two present levels.
  261. !
  262. do k = 2, nlvl-1, 1
  263. if (plvl(k-1) .lt. 200000.) then
  264. if ( (.not. is_there(nint(plvl(k)),'TT')) .and. &
  265. ( is_there(nint(plvl(k-1)), 'TT')) .and.&
  266. ( is_there(nint(plvl(k+1)), 'TT')) ) then
  267. call get_dims(nint(plvl(k+1)), 'TT')
  268. call vntrp(plvl, maxlvl, k, "TT ", map%nx, map%ny)
  269. endif
  270. endif
  271. enddo
  272. !
  273. ! Check to see if we need to fill HGT from GEOPT.
  274. !
  275. do k = 1, nlvl
  276. if (plvl(k).lt.200000.) then
  277. if (.not. is_there(nint(plvl(k)), 'HGT').and. &
  278. is_there(nint(plvl(k)), 'GEOPT')) then
  279. call get_dims(nint(plvl(k)), 'GEOPT')
  280. allocate(scr2d(map%nx,map%ny))
  281. call get_storage(nint(plvl(k)), 'GEOPT', scr2d, map%nx, map%ny)
  282. scr2d = scr2d / 9.81
  283. call put_storage(nint(plvl(k)), 'HGT', scr2d, map%nx, map%ny)
  284. deallocate(scr2d)
  285. endif
  286. endif
  287. enddo
  288. ! If upper-air RH is missing, see if we can compute RH from Specific Humidity:
  289. do k = 1, nlvl
  290. if (plvl(k).lt.200000.) then
  291. if (.not. is_there(nint(plvl(k)), 'RH') .and. &
  292. is_there(nint(plvl(k)), 'TT') .and. &
  293. is_there(nint(plvl(k)), 'SPECHUMD')) then
  294. call get_dims(nint(plvl(k)), 'TT')
  295. call compute_rh_spechumd_upa(map%nx, map%ny, plvl(k))
  296. endif
  297. endif
  298. enddo
  299. ! If upper-air RH is missing, see if we can compute RH from Vapor Pressure:
  300. ! (Thanks to Bob Hart of PSU ESSC -- 1999-05-27.)
  301. do k = 1, nlvl
  302. if (plvl(k).lt.200000.) then
  303. if (.not. is_there(nint(plvl(k)),'RH').and. &
  304. is_there(nint(plvl(k)), 'TT') .and. &
  305. is_there(nint(plvl(k)),'VAPP')) then
  306. call get_dims(nint(plvl(k)),'TT')
  307. call compute_rh_vapp_upa(map%nx, map%ny, plvl(k))
  308. endif
  309. endif
  310. enddo
  311. ! If upper-air RH is missing, see if we can compute RH from Dewpoint Depression:
  312. do k = 1, nlvl
  313. if (plvl(k).lt.200000.) then
  314. if (.not. is_there(nint(plvl(k)),'RH').and. &
  315. is_there(nint(plvl(k)), 'TT') .and. &
  316. is_there(nint(plvl(k)),'DEPR')) then
  317. call get_dims(nint(plvl(k)),'TT')
  318. call compute_rh_depr(map%nx, map%ny, plvl(k))
  319. endif
  320. endif
  321. enddo
  322. !
  323. ! If upper-air RH is missing, see if we can interpolate from surrounding levels.
  324. ! This is a simple vertical interpolation, linear in pressure.
  325. ! Currently, this simply fills in one missing level between two present levels.
  326. ! May expand this in the future to fill in additional levels. May also expand
  327. ! this in the future to vertically interpolate other variables.
  328. !
  329. do k = 2, nlvl-1, 1
  330. if (plvl(k-1) .lt. 200000.) then
  331. if ( (.not. is_there(nint(plvl(k)),'RH')) .and. &
  332. ( is_there(nint(plvl(k-1)), 'RH')) .and.&
  333. ( is_there(nint(plvl(k+1)), 'RH')) ) then
  334. call get_dims(nint(plvl(k+1)), 'RH')
  335. call vntrp(plvl, maxlvl, k, "RH ", map%nx, map%ny)
  336. endif
  337. endif
  338. enddo
  339. ! Repair GFS and ECMWF pressure-level RH
  340. if (index(map%source,'NCEP GFS') .ne. 0 .or. &
  341. index(map%source,'ECMWF') .ne. 0 ) then
  342. call mprintf(.true.,DEBUG, &
  343. "RRPR: Adjusting GFS/ECMWF RH values ")
  344. do k = 1, nlvl
  345. if ( is_there(nint(plvl(k)),'RH') .and. &
  346. is_there(nint(plvl(k)),'TT') ) then
  347. call fix_gfs_rh (map%nx, map%ny, plvl(k))
  348. endif
  349. enddo
  350. endif
  351. !
  352. ! Check to see if we need to fill RH above 300 mb:
  353. !
  354. if (is_there(30000, 'RH')) then
  355. call get_dims(30000, 'RH')
  356. allocate(scr2d(map%nx,map%ny))
  357. do k = 1, nlvl
  358. ! Set missing RH to 5% between 300 and 70 hPa. Set RH to 0 above 70 hPa.
  359. ! The stratospheric RH will be adjusted further in real.
  360. if (plvl(k).le.7000.) then
  361. scr2d = 0.
  362. else if (plvl(k).lt.30000.) then
  363. scr2d = 5.
  364. endif
  365. if (plvl(k).lt.30000. .and. plvl(k) .gt. 10. ) then
  366. ! levels higher than .1 mb are special - do not fill
  367. if (.not. is_there(nint(plvl(k)), 'RH')) then
  368. call put_storage(nint(plvl(k)),'RH',scr2d,map%nx,map%ny)
  369. call mprintf(.true.,DEBUG, &
  370. "RRPR: RH missing at %i hPa, inserting synthetic RH ",i1=nint(plvl(k)/100.))
  371. endif
  372. endif
  373. enddo
  374. deallocate(scr2d)
  375. endif
  376. !
  377. ! If surface RH is missing, see if we can compute RH from Specific Humidity
  378. ! or Dewpoint or Dewpoint depression:
  379. !
  380. if (.not. is_there (200100, 'RH')) then
  381. if (is_there(200100, 'TT').and. &
  382. is_there(200100, 'PSFC' ) .and. &
  383. is_there(200100, 'SPECHUMD')) then
  384. call get_dims(200100, 'TT')
  385. call compute_rh_spechumd(map%nx, map%ny)
  386. call mprintf(.true.,DEBUG, &
  387. "RRPR: SURFACE RH is computed")
  388. elseif (is_there(200100, 'TT' ).and. &
  389. is_there(200100, 'DEWPT')) then
  390. call get_dims(200100, 'TT')
  391. call compute_rh_dewpt(map%nx, map%ny)
  392. elseif (is_there(200100, 'TT').and. &
  393. is_there(200100, 'DEPR')) then
  394. call get_dims(200100, 'TT')
  395. call compute_rh_depr(map%nx, map%ny, 200100.)
  396. endif
  397. endif
  398. !
  399. ! If surface SNOW is missing, see if we can compute SNOW from SNOWRUC
  400. ! (From Wei Wang, 2007 June 21, modified 12/28/2007)
  401. !
  402. if (.not. is_there(200100, 'SNOW') .and. &
  403. is_there(200100, 'SNOWRUC')) then
  404. call get_dims(200100, 'SNOWRUC')
  405. allocate(scr2d(map%nx,map%ny))
  406. call get_storage(200100, 'SNOWRUC', scr2d, map%nx, map%ny)
  407. scr2d = scr2d * 1000.
  408. call put_storage(200100, 'SNOW', scr2d, map%nx, map%ny)
  409. deallocate(scr2d)
  410. endif
  411. !
  412. ! Check to see if we need to fill SOILHGT from SOILGEO.
  413. ! (From Wei Wang, 2007 June 21)
  414. !
  415. if (.not. is_there(200100, 'SOILHGT') .and. &
  416. is_there(200100, 'SOILGEO')) then
  417. call get_dims(200100, 'SOILGEO')
  418. allocate(scr2d(map%nx,map%ny))
  419. call get_storage(200100, 'SOILGEO', scr2d, map%nx, map%ny)
  420. scr2d = scr2d / 9.81
  421. call put_storage(200100, 'SOILHGT', scr2d, map%nx, map%ny)
  422. deallocate(scr2d)
  423. endif
  424. ! For hybrid-level input, soilgeo is in level 1 (e.g. ERA40)
  425. if (.not. is_there(200100, 'SOILHGT') .and. &
  426. is_there(1, 'SOILGEO')) then
  427. call get_dims(1, 'SOILGEO')
  428. allocate(scr2d(map%nx,map%ny))
  429. call get_storage(1, 'SOILGEO', scr2d, map%nx, map%ny)
  430. scr2d = scr2d / 9.81
  431. call put_storage(200100, 'SOILHGT', scr2d, map%nx, map%ny)
  432. deallocate(scr2d)
  433. endif
  434. ! For ECMWF hybrid-level input, may need to move psfc from level 1 to 2001.
  435. if ( index(map%source,'ECMWF') .ne. 0) then
  436. if (.not. is_there(200100, 'PSFC') .and. &
  437. is_there(1, 'PSFCH')) then
  438. call get_dims(1, 'PSFCH')
  439. allocate(scr2d(map%nx,map%ny))
  440. call get_storage(1, 'PSFCH', scr2d, map%nx, map%ny)
  441. call put_storage(200100, 'PSFC', scr2d, map%nx, map%ny)
  442. deallocate(scr2d)
  443. endif
  444. endif
  445. ! ECMWF snow depth in meters of water equivalent (Table 128). Convert to kg/m2
  446. !
  447. if (is_there(200100, 'SNOW_EC')) then
  448. call get_dims(200100, 'SNOW_EC')
  449. allocate(scr2d(map%nx,map%ny))
  450. call get_storage(200100, 'SNOW_EC', scr2d, map%nx, map%ny)
  451. scr2d = scr2d * 1000.
  452. call put_storage(200100, 'SNOW', scr2d, map%nx, map%ny)
  453. deallocate(scr2d)
  454. endif
  455. ! Convert the ECMWF LANDSEA mask from a fraction to a flag
  456. if ( index(map%source,'ECMWF') .ne. 0) then
  457. if (is_there(200100, 'LANDSEA')) then
  458. call get_dims(200100, 'LANDSEA')
  459. call make_zero_or_one(map%nx, map%ny, 'LANDSEA')
  460. endif
  461. endif
  462. ! NCEP GFS weasd is one-half of the NAM value. Increase it for use in WRF.
  463. ! The GFS-based reanalyses values should be OK as is.
  464. if ((index(map%source,'NCEP GFS') .ne. 0 .or. &
  465. index(map%source,'NCEP GEFS') .ne. 0) .and. &
  466. is_there(200100, 'SNOW')) then
  467. call mprintf(.true.,DEBUG, &
  468. "RRPR: Recomputing SNOW for NCEP GFS")
  469. call get_dims(200100, 'SNOW')
  470. allocate(scr2d(map%nx,map%ny))
  471. call get_storage(200100, 'SNOW', scr2d, map%nx, map%ny)
  472. scr2d = scr2d * 2.
  473. call put_storage(200100, 'SNOW', scr2d, map%nx, map%ny)
  474. deallocate(scr2d)
  475. endif
  476. ! compute physical snow depth (SNOWH) for various models
  477. ! As of March 2011, this is done here instead of real because we have model
  478. ! source information.
  479. if (is_there(200100, 'SNOW') .and. .not. is_there(200100, 'SNOWH')) then
  480. call get_dims(200100, 'SNOW')
  481. allocate(scr2d(map%nx,map%ny))
  482. call get_storage(200100, 'SNOW', scr2d, map%nx, map%ny)
  483. call mprintf(.true.,DEBUG, &
  484. "RRPR: Computing SNOWH from SNOW")
  485. if ( index(map%source,'NCEP ') .ne. 0) then
  486. scr2d = scr2d * 0.005 ! Assume 200:1 ratio as used at NCEP and in NOAH
  487. else if (index(map%source,'ECMWF') .ne. 0) then
  488. if (is_there(200100, 'SNOW_DEN')) then ! If we have snow density, use it to compute snowh
  489. call get_dims(200100, 'SNOW_DEN')
  490. allocate(tmp2d(map%nx,map%ny))
  491. call get_storage(200100, 'SNOW_DEN', tmp2d, map%nx, map%ny)
  492. scr2d = scr2d / tmp2d
  493. deallocate(tmp2d)
  494. else
  495. scr2d = scr2d * 0.004 ! otherwise, assume a density of 250 mm/m (i.e. 250:1 ratio).
  496. endif
  497. else ! Other models
  498. scr2d = scr2d * 0.005 ! Use real's default method (200:1)
  499. endif
  500. call put_storage(200100, 'SNOWH', scr2d, map%nx, map%ny)
  501. deallocate(scr2d)
  502. endif
  503. ! As of March 2011, SEAICE can be a flag or a fraction. It will be converted
  504. ! to the appropriate values in real depending on whether or not the polar mods are used.
  505. !! If we've got a SEAICE field, make sure that it is all Zeros and Ones:
  506. ! if (is_there(200100, 'SEAICE')) then
  507. ! call get_dims(200100, 'SEAICE')
  508. ! call make_zero_or_one(map%nx, map%ny, 'SEAICE')
  509. ! endif
  510. ! If we've got an ICEMASK field, re-flag it for output to met_em and real:
  511. ! Field | GRIB In | Out
  512. ! -------------------------
  513. ! water | 0 | 0
  514. ! land | -1 | 1
  515. ! ice | 1 | 0
  516. if (is_there(200100, 'ICEMASK')) then
  517. call get_dims(200100, 'ICEMASK')
  518. call re_flag_ice_mask(map%nx, map%ny)
  519. endif
  520. ! If we have an ICEFRAC field, convert from % to fraction
  521. if (is_there(200100, 'ICEFRAC')) then
  522. call get_dims(200100, 'ICEFRAC')
  523. allocate(scr2d(map%nx,map%ny))
  524. call get_storage(200100, 'ICEFRAC', scr2d, map%nx, map%ny)
  525. scr2d = scr2d / 100.
  526. call put_storage(200100, 'ICEFRAC', scr2d, map%nx, map%ny)
  527. deallocate(scr2d)
  528. endif
  529. call mprintf(.true.,INFORM, &
  530. "RRPR: hdate = %s ",s1=hdate)
  531. call output(hdate, nlvl, maxlvl, plvl, interval, 2, out_format, prefix, debug_level)
  532. call clear_storage
  533. exit FILELOOP
  534. enddo FILELOOP
  535. enddo TIMELOOP
  536. end subroutine rrpr
  537. subroutine make_zero_or_one(ix, jx, infield)
  538. ! Make sure the input field (SEAICE or LANDSEA) is zero or one.
  539. use storage_module
  540. implicit none
  541. integer :: ix, jx
  542. real, dimension(ix,jx) :: seaice
  543. character(len=*) :: infield
  544. call get_storage(200100, infield, seaice, ix, jx)
  545. where(seaice > 0.5)
  546. seaice = 1.0
  547. elsewhere
  548. seaice = 0.0
  549. end where
  550. call put_storage(200100, infield, seaice, ix, jx)
  551. end subroutine make_zero_or_one
  552. subroutine re_flag_ice_mask(ix, jx)
  553. !
  554. ! Change land points from -1 to 1
  555. ! Change ice points from 1 to 0
  556. ! Water points stay at 0
  557. !
  558. use storage_module
  559. implicit none
  560. integer :: ix, jx
  561. real, dimension(ix,jx) :: iceflag
  562. call get_storage(200100, 'ICEMASK',iceflag, ix, jx)
  563. where(iceflag > 0.5) ! Ice points, set to water value
  564. iceflag = 0.0
  565. end where
  566. where(iceflag < -0.5) ! Land points
  567. iceflag = 1.0
  568. end where
  569. call put_storage(200100, 'ICEMASK',iceflag, ix, jx)
  570. end subroutine re_flag_ice_mask
  571. subroutine compute_spechumd_qvapor(ix, jx, plvl)
  572. ! Compute specific humidity from water vapor mixing ratio.
  573. use storage_module
  574. implicit none
  575. integer :: ix, jx
  576. real :: plvl
  577. real, dimension(ix,jx) :: QVAPOR, SPECHUMD
  578. call get_storage(nint(plvl), 'QV', QVAPOR, ix, jx)
  579. SPECHUMD = QVAPOR/(1.+QVAPOR)
  580. call put_storage(nint(plvl), 'SPECHUMD', spechumd, ix, jx)
  581. if(nint(plvl).eq.1) then
  582. call put_storage(200100,'SPECHUMD', spechumd, ix, jx)
  583. endif
  584. end subroutine compute_spechumd_qvapor
  585. subroutine compute_t_vptmp(ix, jx, plvl)
  586. ! Compute temperature from virtual potential temperature
  587. use storage_module
  588. implicit none
  589. integer :: ix, jx
  590. real :: plvl
  591. real, dimension(ix,jx) :: T, VPTMP, P, Q
  592. real, parameter :: rovcp=0.28571
  593. call get_storage(nint(plvl), 'VPTMP', VPTMP, ix, jx)
  594. IF (nint(plvl) .LT. 200) THEN
  595. call get_storage(nint(plvl), 'PRESSURE', P, ix, jx)
  596. ELSE
  597. p = plvl
  598. ENDIF
  599. call get_storage(nint(plvl), 'SPECHUMD', Q, ix, jx)
  600. t=vptmp * (p*1.e-5)**rovcp * (1./(1.+0.6078*Q))
  601. call put_storage(nint(plvl), 'TT', t, ix, jx)
  602. if(nint(plvl).eq.1) then
  603. call put_storage(200100, 'PSFC', p, ix, jx)
  604. endif
  605. end subroutine compute_t_vptmp
  606. subroutine compute_rh_spechumd(ix, jx)
  607. ! Compute relative humidity from specific humidity.
  608. use storage_module
  609. implicit none
  610. integer :: ix, jx
  611. real, dimension(ix,jx) :: T, P, RH, Q
  612. real, parameter :: svp1=611.2
  613. real, parameter :: svp2=17.67
  614. real, parameter :: svp3=29.65
  615. real, parameter :: svpt0=273.15
  616. real, parameter :: eps = 0.622
  617. call get_storage(200100, 'TT', T, ix, jx)
  618. call get_storage(200100, 'PSFC', P, ix, jx)
  619. call get_storage(200100, 'SPECHUMD', Q, ix, jx)
  620. rh = 1.E2 * (p*q/(q*(1.-eps) + eps))/(svp1*exp(svp2*(t-svpt0)/(T-svp3)))
  621. call put_storage(200100, 'RH', rh, ix, jx)
  622. end subroutine compute_rh_spechumd
  623. subroutine compute_rh_spechumd_upa(ix, jx, plvl)
  624. ! Compute relative humidity from specific humidity.
  625. use storage_module
  626. implicit none
  627. integer :: ix, jx
  628. real :: plvl
  629. real, dimension(ix,jx) :: T, P, RH, Q
  630. real, parameter :: svp1=611.2
  631. real, parameter :: svp2=17.67
  632. real, parameter :: svp3=29.65
  633. real, parameter :: svpt0=273.15
  634. real, parameter :: eps = 0.622
  635. IF ( nint(plvl).LT. 200) THEN
  636. if (is_there(nint(plvl), 'PRESSURE')) then
  637. call get_storage(nint(plvl), 'PRESSURE', P, ix, jx)
  638. else
  639. return ! if we don't have pressure on model levels, return
  640. endif
  641. ELSE
  642. P = plvl
  643. ENDIF
  644. call get_storage(nint(plvl), 'TT', T, ix, jx)
  645. call get_storage(nint(plvl), 'SPECHUMD', Q, ix, jx)
  646. Q=MAX(1.E-10,Q)
  647. rh = 1.E2 * (p*q/(q*(1.-eps) + eps))/(svp1*exp(svp2*(t-svpt0)/(T-svp3)))
  648. call put_storage(nint(plvl), 'RH', rh, ix, jx)
  649. end subroutine compute_rh_spechumd_upa
  650. subroutine compute_rh_vapp_upa(ix, jx, plvl)
  651. ! Compute relative humidity from vapor pressure.
  652. ! Thanks to Bob Hart of PSU ESSC -- 1999-05-27.
  653. use storage_module
  654. implicit none
  655. integer :: ix, jx
  656. real :: plvl
  657. real, dimension(ix,jx) :: P, ES
  658. real, pointer, dimension(:,:) :: T, E, RH
  659. real, parameter :: svp1=611.2
  660. real, parameter :: svp2=17.67
  661. real, parameter :: svp3=29.65
  662. real, parameter :: svpt0=273.15
  663. allocate(RH(ix,jx))
  664. IF ( nint(plvl).LT. 200) THEN
  665. if (is_there(nint(plvl), 'PRESSURE')) then
  666. call get_storage(nint(plvl), 'PRESSURE', P, ix, jx)
  667. else
  668. return ! if we don't have pressure on model levels, return
  669. endif
  670. ELSE
  671. P = plvl
  672. ENDIF
  673. call refr_storage(nint(plvl), 'TT', T, ix, jx)
  674. call refr_storage(nint(plvl), 'VAPP', E, ix, jx)
  675. ES=svp1*exp(svp2*(T-svpt0)/(T-svp3))
  676. rh=min(1.E2*(P-ES)*E/((P-E)*ES), 1.E2)
  677. call refw_storage(nint(plvl), 'RH', rh, ix, jx)
  678. nullify(T,E)
  679. end subroutine compute_rh_vapp_upa
  680. subroutine compute_rh_depr(ix, jx, plvl)
  681. ! Compute relative humidity from Dewpoint Depression
  682. use storage_module
  683. implicit none
  684. integer :: ix, jx
  685. real :: plvl
  686. real, dimension(ix,jx) :: t, depr, rh
  687. real, parameter :: Xlv = 2.5e6
  688. real, parameter :: Rv = 461.5
  689. integer :: i, j
  690. call get_storage(nint(plvl), 'TT', T, ix, jx)
  691. call get_storage(nint(plvl), 'DEPR', DEPR, ix, jx)
  692. where(DEPR < 100.)
  693. rh = exp(Xlv/Rv*(1./T - 1./(T-depr))) * 1.E2
  694. elsewhere
  695. rh = 0.0
  696. endwhere
  697. call put_storage(nint(plvl),'RH ', rh, ix, jx)
  698. end subroutine compute_rh_depr
  699. subroutine compute_rh_dewpt(ix,jx)
  700. ! Compute relative humidity from Dewpoint
  701. use storage_module
  702. implicit none
  703. integer :: ix, jx
  704. real, dimension(ix,jx) :: t, dp, rh
  705. real, parameter :: Xlv = 2.5e6
  706. real, parameter :: Rv = 461.5
  707. call get_storage(200100, 'TT ', T, ix, jx)
  708. call get_storage(200100, 'DEWPT ', DP, ix, jx)
  709. rh = exp(Xlv/Rv*(1./T - 1./dp)) * 1.E2
  710. call put_storage(200100,'RH ', rh, ix, jx)
  711. end subroutine compute_rh_dewpt
  712. subroutine vntrp(plvl, maxlvl, k, name, ix, jx)
  713. use storage_module
  714. implicit none
  715. integer :: ix, jx, k, maxlvl
  716. real, dimension(maxlvl) :: plvl
  717. character(len=8) :: name
  718. real, dimension(ix,jx) :: a, b, c
  719. real :: frc
  720. write(*,'("Interpolating to fill in ", A, " at level ", I8)') trim(name), nint(plvl(k))
  721. call get_storage(nint(plvl(k-1)), name, a, ix, jx)
  722. call get_storage(nint(plvl(k+1)), name, c, ix, jx)
  723. frc = (plvl(k) - plvl(k+1)) / ( plvl(k-1)-plvl(k+1))
  724. b = (1.-frc)*a + frc*c
  725. !KWM b = 0.5 * (a + c)
  726. call put_storage(nint(plvl(k)), name, b, ix, jx)
  727. end subroutine vntrp
  728. subroutine fix_gfs_rh (ix, jx, plvl)
  729. ! This routine replaces GFS RH (wrt ice) with RH wrt liquid (which is what is assumed in real.exe).
  730. use storage_module
  731. implicit none
  732. integer :: ix, jx, i, j
  733. real :: plvl, eis, ews, r
  734. real, allocatable, dimension(:,:) :: rh, tt
  735. allocate(rh(ix,jx))
  736. allocate(tt(ix,jx))
  737. call get_storage(nint(plvl), 'RH', rh, ix, jx)
  738. call get_storage(nint(plvl), 'TT', tt, ix, jx)
  739. do j = 1, jx
  740. do i = 1, ix
  741. if ( tt(i,j) .le. 273.15 ) then
  742. ! Murphy and Koop 2005 ice saturation vapor pressure.
  743. ! eis and ews in hPA, tt is in K
  744. eis = .01 * exp (9.550426 - (5723.265 / tt(i,j)) + (3.53068 * alog(tt(i,j))) &
  745. - (0.00728332 * tt(i,j)))
  746. ! Bolton 1980 liquid saturation vapor pressure. For water saturation, most
  747. ! formulae are very similar from 0 to -20, so we don't need a more exact formula.
  748. ews = 6.112 * exp(17.67 * (tt(i,j)-273.15) / ((tt(i,j)-273.15)+243.5))
  749. if ( tt(i,j) .gt. 253.15 ) then
  750. ! A linear approximation to the GFS blending region ( -20 > T < 0 )
  751. r = ((273.15 - tt(i,j)) / 20.)
  752. r = (r * eis) + ((1-r)*ews)
  753. else
  754. r = eis
  755. endif
  756. rh(i,j) = rh(i,j) * (r / ews)
  757. endif
  758. enddo
  759. enddo
  760. call put_storage(nint(plvl), 'RH', rh, ix, jx)
  761. deallocate (rh)
  762. deallocate (tt)
  763. end subroutine fix_gfs_rh