PageRenderTime 516ms CodeModel.GetById 34ms RepoModel.GetById 8ms app.codeStats 0ms

/WPS/geogrid/util/plotgrid/src/plotgrid.F

http://github.com/jbeezley/wrf-fire
FORTRAN Legacy | 406 lines | 226 code | 73 blank | 107 comment | 41 complexity | 21f9bbc23d52f014e8d014a60c181323 MD5 | raw file
Possible License(s): AGPL-1.0
  1. program plotgrid
  2. use input_module
  3. implicit none
  4. external ulpr
  5. integer :: n, i, j, nx, ny
  6. integer :: istatus, start_mem_i, end_mem_i, start_mem_j, end_mem_j, &
  7. start_mem_k, end_mem_k, dyn_opt, &
  8. west_east_dim, south_north_dim, bottom_top_dim, map_proj, is_water, num_land_cat, &
  9. is_ice, is_urban, isoilwater, grid_id, parent_id, i_parent_start, j_parent_start, &
  10. i_parent_end, j_parent_end, parent_grid_ratio, &
  11. we_patch_s, we_patch_e, we_patch_s_stag, we_patch_e_stag, &
  12. sn_patch_s, sn_patch_e, sn_patch_s_stag, sn_patch_e_stag
  13. real :: width, height
  14. real :: dx, dy, cen_lat, moad_cen_lat, cen_lon, stand_lon, truelat1, truelat2, pole_lat, pole_lon
  15. real :: start_r, start_g, start_b, end_r, end_g, end_b
  16. real :: ll_lat, ll_lon, ur_lat, ur_lon
  17. real :: left, right, bottom, top, maxter, minter
  18. real :: rotang
  19. real, dimension(16) :: corner_lats, corner_lons
  20. real, dimension(10000) :: xcs, ycs
  21. integer, dimension(10) :: iai, iag
  22. integer, dimension(400000) :: iam
  23. integer, allocatable, dimension(:,:) :: lu
  24. real, allocatable, dimension(:,:) :: xlat, xlon, ter
  25. real, dimension(122000) :: rwrk
  26. real, pointer, dimension(:,:,:) :: real_array
  27. character (len=3) :: memorder
  28. character (len=25) :: crotang
  29. character (len=25) :: units
  30. character (len=46) :: desc
  31. character (len=128) :: init_date, cname, stagger, cunits, cdesc, title, startdate, grid_type, mminlu
  32. character (len=128), dimension(3) :: dimnames
  33. call getarg(1,crotang)
  34. read (crotang,'(f)') rotang
  35. write(6,*) 'Plotting with rotation angle ',rotang
  36. call opngks
  37. call gopwk(13, 41, 3)
  38. call gscr(1, 0, 1.00, 1.00, 1.00)
  39. call gscr(1, 1, 0.00, 0.00, 0.00)
  40. call gscr(1, 2, 0.25, 0.25, 0.25)
  41. call gscr(1, 3, 1.00, 1.00, 0.50)
  42. call gscr(1, 4, 0.50, 1.00, 0.50)
  43. call gscr(1, 5, 1.00, 1.00, 0.00)
  44. call gscr(1, 6, 1.00, 1.00, 0.00)
  45. call gscr(1, 7, 0.50, 1.00, 0.50)
  46. call gscr(1, 8, 1.00, 1.00, 0.50)
  47. call gscr(1, 9, 0.50, 1.00, 0.50)
  48. call gscr(1,10, 0.50, 1.00, 0.50)
  49. call gscr(1,11, 1.00, 1.00, 0.50)
  50. call gscr(1,12, 0.00, 1.00, 0.00)
  51. call gscr(1,13, 0.00, 0.50, 0.00)
  52. call gscr(1,14, 0.00, 1.00, 0.00)
  53. call gscr(1,15, 0.00, 0.50, 0.00)
  54. call gscr(1,16, 0.00, 1.00, 0.00)
  55. call gscr(1,17, 0.50, 0.50, 1.00)
  56. call gscr(1,18, 0.00, 1.00, 0.00)
  57. call gscr(1,19, 0.00, 1.00, 0.00)
  58. call gscr(1,20, 0.75, 0.75, 0.75)
  59. call gscr(1,21, 0.75, 0.75, 0.75)
  60. call gscr(1,22, 0.00, 0.50, 0.00)
  61. call gscr(1,23, 0.75, 0.75, 0.75)
  62. call gscr(1,24, 0.75, 0.75, 0.75)
  63. call gscr(1,25, 1.00, 1.00, 1.00)
  64. start_r = 0.00
  65. end_r = 0.50
  66. start_g = 1.00
  67. end_g = 0.25
  68. start_b = 0.00
  69. end_b = 0.00
  70. do i=26,76
  71. call gscr(1,i,start_r+((end_r-start_r)/50.)*real(i-26),start_g+((end_g-start_g)/50.)*real(i-26),start_b+((end_b-start_b)/50.)*real(i-26))
  72. end do
  73. start_r = 0.50
  74. end_r = 1.00
  75. start_g = 0.25
  76. end_g = 1.00
  77. start_b = 0.00
  78. end_b = 1.00
  79. do i=77,126
  80. call gscr(1,i,start_r+((end_r-start_r)/50.)*real(i-77),start_g+((end_g-start_g)/50.)*real(i-77),start_b+((end_b-start_b)/50.)*real(i-77))
  81. end do
  82. start_r = 0.80
  83. end_r = 1.00
  84. start_g = 0.80
  85. end_g = 1.00
  86. start_b = 0.80
  87. end_b = 1.00
  88. do i=127,176
  89. call gscr(1,i,start_r+((end_r-start_r)/50.)*real(i-127),start_g+((end_g-start_g)/50.)*real(i-127),start_b+((end_b-start_b)/50.)*real(i-127))
  90. end do
  91. call get_namelist_params()
  92. do n=1,max_dom
  93. call input_init(n, istatus)
  94. if (istatus /= 0) then
  95. write(6,*) ' '
  96. write(6,*) 'Error: Could not open domain01 file.'
  97. write(6,*) ' '
  98. stop
  99. end if
  100. call read_global_attrs(title, init_date, grid_type, dyn_opt, &
  101. west_east_dim, south_north_dim, bottom_top_dim, &
  102. we_patch_s, we_patch_e, we_patch_s_stag, we_patch_e_stag, &
  103. sn_patch_s, sn_patch_e, sn_patch_s_stag, sn_patch_e_stag, &
  104. map_proj, mminlu, num_land_cat, is_water, &
  105. is_ice, is_urban, isoilwater, grid_id, parent_id, i_parent_start, j_parent_start, &
  106. i_parent_end, j_parent_end, dx, dy, cen_lat, moad_cen_lat, cen_lon, &
  107. stand_lon, truelat1, truelat2, pole_lat, pole_lon, parent_grid_ratio, &
  108. corner_lats, corner_lons)
  109. istatus = 0
  110. do while (istatus == 0)
  111. call read_next_field(start_mem_i, end_mem_i, start_mem_j, end_mem_j, &
  112. start_mem_k, end_mem_k, cname, cunits, cdesc, &
  113. memorder, stagger, dimnames, real_array, istatus)
  114. if (istatus == 0) then
  115. if (index(cname, 'XLAT_M') /= 0) then
  116. nx = end_mem_i - start_mem_i + 1
  117. ny = end_mem_j - start_mem_j + 1
  118. allocate(xlat(nx,ny))
  119. xlat = real_array(:,:,1)
  120. else if (index(cname, 'XLONG_M') /= 0) then
  121. nx = end_mem_i - start_mem_i + 1
  122. ny = end_mem_j - start_mem_j + 1
  123. allocate(xlon(nx,ny))
  124. xlon = real_array(:,:,1)
  125. else if (index(cname, 'LU_INDEX') /= 0) then
  126. nx = end_mem_i - start_mem_i + 1
  127. ny = end_mem_j - start_mem_j + 1
  128. allocate(lu(nx,ny))
  129. lu = nint(real_array(:,:,1))
  130. end if
  131. end if
  132. end do
  133. call input_close()
  134. ll_lat = xlat(1,1)
  135. ll_lon = xlon(1,1)
  136. ur_lat = xlat(nx,ny)
  137. ur_lon = xlon(nx,ny)
  138. ! if (ur_lon < 0.) ur_lon = ur_lon + 360.0
  139. if (n == 1) then
  140. left = 0.0
  141. right = 1.0
  142. bottom = 0.0
  143. top = 1.0
  144. call mappos(left,right,bottom,top)
  145. call mapstc('OU','CO')
  146. call maproj('CE', cen_lat, cen_lon, rotang)
  147. ! call maproj('LC', truelat1, stand_lon, truelat2)
  148. ! call maproj('ST', cen_lat, cen_lon, stand_lon)
  149. call mapset('CO', ll_lat, ll_lon, ur_lat, ur_lon)
  150. call mapint()
  151. end if
  152. call mpsetr('GR', 10.0)
  153. call maptrn(ll_lat, ll_lon, left, bottom)
  154. call maptrn(ur_lat, ur_lon, right, top)
  155. width = 1.02*(right-left)/real(nx)
  156. height = 1.02*(top-bottom)/real(ny)
  157. do j=1,ny
  158. do i=1,nx
  159. call map_square(xlat(i,j), xlon(i,j), width, height, lu(i,j)+1)
  160. end do
  161. end do
  162. if (n > 1) then
  163. call gsplci(0)
  164. call lined(left-width/2.,bottom-height/2.,left-width/2.,top+height/2.)
  165. call lined(left-width/2.,top+height/2.,right+width/2.,top+height/2.)
  166. call lined(right+width/2.,top+height/2.,right+width/2.,bottom-height/2.)
  167. call lined(right+width/2.,bottom-height/2.,left-width/2.,bottom-height/2.)
  168. call sflush()
  169. call gsplci(1)
  170. end if
  171. deallocate(xlat)
  172. deallocate(xlon)
  173. deallocate(lu)
  174. end do
  175. call mplndr('Earth..3',4)
  176. call arinam (iam,400000)
  177. call mapbla (iam)
  178. call arpram (iam,0,0,0)
  179. call mapgrm (iam,xcs,ycs,10000,iai,iag,10,ulpr)
  180. call frame()
  181. do n=1,max_dom
  182. call input_init(n, istatus)
  183. if (istatus /= 0) then
  184. write(6,*) ' '
  185. write(6,*) 'Error: Could not open domain01 file.'
  186. write(6,*) ' '
  187. stop
  188. end if
  189. call read_global_attrs(title, init_date, grid_type, dyn_opt, &
  190. west_east_dim, south_north_dim, bottom_top_dim, &
  191. we_patch_s, we_patch_e, we_patch_s_stag, we_patch_e_stag, &
  192. sn_patch_s, sn_patch_e, sn_patch_s_stag, sn_patch_e_stag, &
  193. map_proj, mminlu, num_land_cat, is_water, &
  194. is_ice, is_urban, isoilwater, grid_id, parent_id, i_parent_start, j_parent_start, &
  195. i_parent_end, j_parent_end, dx, dy, cen_lat, moad_cen_lat, cen_lon, &
  196. stand_lon, truelat1, truelat2, pole_lat, pole_lon, parent_grid_ratio, &
  197. corner_lats, corner_lons)
  198. istatus = 0
  199. do while (istatus == 0)
  200. call read_next_field(start_mem_i, end_mem_i, start_mem_j, end_mem_j, &
  201. start_mem_k, end_mem_k, cname, cunits, cdesc, &
  202. memorder, stagger, dimnames, real_array, istatus)
  203. if (istatus == 0) then
  204. if (index(cname, 'XLAT_M') /= 0) then
  205. nx = end_mem_i - start_mem_i + 1
  206. ny = end_mem_j - start_mem_j + 1
  207. allocate(xlat(nx,ny))
  208. xlat = real_array(:,:,1)
  209. else if (index(cname, 'XLONG_M') /= 0) then
  210. nx = end_mem_i - start_mem_i + 1
  211. ny = end_mem_j - start_mem_j + 1
  212. allocate(xlon(nx,ny))
  213. xlon = real_array(:,:,1)
  214. else if (index(cname, 'HGT_M') /= 0) then
  215. nx = end_mem_i - start_mem_i + 1
  216. ny = end_mem_j - start_mem_j + 1
  217. allocate(ter(nx,ny))
  218. ter = real_array(:,:,1)
  219. else if (index(cname, 'LU_INDEX') /= 0) then
  220. nx = end_mem_i - start_mem_i + 1
  221. ny = end_mem_j - start_mem_j + 1
  222. allocate(lu(nx,ny))
  223. lu = nint(real_array(:,:,1))
  224. end if
  225. end if
  226. end do
  227. call input_close()
  228. ll_lat = xlat(1,1)
  229. ll_lon = xlon(1,1)
  230. ur_lat = xlat(nx,ny)
  231. ur_lon = xlon(nx,ny)
  232. if (n == 1) then
  233. left = 0.0
  234. right = 1.0
  235. bottom = 0.0
  236. top = 1.0
  237. call mappos(left,right,bottom,top)
  238. call mapstc('OU','CO')
  239. call maproj('CE', cen_lat, cen_lon, rotang)
  240. ! call maproj('LC', truelat1, stand_lon, truelat2)
  241. ! call maproj('ST', cen_lat, cen_lon, stand_lon)
  242. call mapset('CO', ll_lat, ll_lon, ur_lat, ur_lon)
  243. call mapint()
  244. maxter = -10000.
  245. minter = 10000.
  246. do j=1,ny
  247. do i=1,nx
  248. if (ter(i,j) > maxter) maxter = ter(i,j)
  249. if (ter(i,j) < minter) minter = ter(i,j)
  250. end do
  251. end do
  252. ! maxter = 3348.42
  253. end if
  254. call maptrn(ll_lat, ll_lon, left, bottom)
  255. call maptrn(ur_lat, ur_lon, right, top)
  256. width = 1.02*(right-left)/real(nx)
  257. height = 1.02*(top-bottom)/real(ny)
  258. do j=1,ny
  259. do i=1,nx
  260. if (lu(i,j) == 16) then
  261. ter(i,j) = ((ter(i,j)-minter) * 99.)/(maxter-minter) + 26.
  262. call map_square(xlat(i,j), xlon(i,j), width, height, 17)
  263. else if (lu(i,j) == 1) then
  264. ter(i,j) = ((ter(i,j)-minter) * 99.)/(maxter-minter) + 26.
  265. call map_square(xlat(i,j), xlon(i,j), width, height, 2)
  266. else if (lu(i,j) == 24) then
  267. ter(i,j) = ((ter(i,j)-minter) * 50.)/(3500.0-minter) + 127.
  268. call map_square(xlat(i,j), xlon(i,j), width, height, nint(ter(i,j)))
  269. else
  270. ter(i,j) = ((ter(i,j)-minter) * 99.)/(maxter-minter) + 26.
  271. call map_square(xlat(i,j), xlon(i,j), width, height, nint(ter(i,j)))
  272. end if
  273. end do
  274. end do
  275. if (n > 1) then
  276. call gsplci(0)
  277. call lined(left-width/2.,bottom-height/2.,left-width/2.,top+height/2.)
  278. call lined(left-width/2.,top+height/2.,right+width/2.,top+height/2.)
  279. call lined(right+width/2.,top+height/2.,right+width/2.,bottom-height/2.)
  280. call lined(right+width/2.,bottom-height/2.,left-width/2.,bottom-height/2.)
  281. call sflush()
  282. call gsplci(1)
  283. end if
  284. deallocate(xlat)
  285. deallocate(xlon)
  286. deallocate(ter)
  287. deallocate(lu)
  288. end do
  289. call mplndr('Earth..3',4)
  290. call arinam (iam,400000)
  291. call mapbla (iam)
  292. call arpram (iam,0,0,0)
  293. call mapgrm (iam,xcs,ycs,10000,iai,iag,10,ulpr)
  294. call gclwk(13)
  295. call clsgks
  296. stop
  297. end program plotgrid
  298. subroutine map_square(rlat, rlon, width, height, colr)
  299. implicit none
  300. ! Arguments
  301. real :: rlat, rlon, width, height
  302. integer :: colr
  303. ! Local variables
  304. real :: u, v
  305. real, dimension(4) :: xra, yra
  306. real, dimension(2000) :: dst
  307. integer, dimension(3000) :: ind
  308. call maptrn(rlat, rlon, u, v)
  309. xra(1) = u-(width/2.)
  310. xra(2) = u+(width/2.)
  311. xra(3) = u+(width/2.)
  312. xra(4) = u-(width/2.)
  313. yra(1) = v-(height/2.)
  314. yra(2) = v-(height/2.)
  315. yra(3) = v+(height/2.)
  316. yra(4) = v+(height/2.)
  317. call sfsgfa(xra, yra, 4, dst, 2000, ind, 3000, colr)
  318. end subroutine map_square
  319. subroutine ulpr(xcs,ycs,ncs,iai,iag,nai)
  320. implicit none
  321. integer, external :: mapaci
  322. integer :: ncs, nai
  323. integer, dimension(nai) :: iai, iag
  324. real, dimension(ncs) :: xcs, ycs
  325. integer :: itm
  326. if (iai(1) >= 0 .and.iai(2) >= 0) then
  327. itm = max0(iai(1),iai(2))
  328. if (mapaci(itm) == 1) then
  329. if (ncs.gt.150) print * , 'ulpr - ncs too big - ',ncs
  330. call gpl(ncs,xcs,ycs)
  331. end if
  332. end if
  333. end subroutine ulpr