/WPS/geogrid/util/plotgrid/src/plotgrid.F
FORTRAN Legacy | 406 lines | 226 code | 73 blank | 107 comment | 41 complexity | 21f9bbc23d52f014e8d014a60c181323 MD5 | raw file
Possible License(s): AGPL-1.0
- program plotgrid
- use input_module
- implicit none
- external ulpr
- integer :: n, i, j, nx, ny
- integer :: istatus, start_mem_i, end_mem_i, start_mem_j, end_mem_j, &
- start_mem_k, end_mem_k, dyn_opt, &
- west_east_dim, south_north_dim, bottom_top_dim, map_proj, is_water, num_land_cat, &
- is_ice, is_urban, isoilwater, grid_id, parent_id, i_parent_start, j_parent_start, &
- i_parent_end, j_parent_end, parent_grid_ratio, &
- we_patch_s, we_patch_e, we_patch_s_stag, we_patch_e_stag, &
- sn_patch_s, sn_patch_e, sn_patch_s_stag, sn_patch_e_stag
- real :: width, height
- real :: dx, dy, cen_lat, moad_cen_lat, cen_lon, stand_lon, truelat1, truelat2, pole_lat, pole_lon
- real :: start_r, start_g, start_b, end_r, end_g, end_b
- real :: ll_lat, ll_lon, ur_lat, ur_lon
- real :: left, right, bottom, top, maxter, minter
- real :: rotang
- real, dimension(16) :: corner_lats, corner_lons
- real, dimension(10000) :: xcs, ycs
- integer, dimension(10) :: iai, iag
- integer, dimension(400000) :: iam
- integer, allocatable, dimension(:,:) :: lu
- real, allocatable, dimension(:,:) :: xlat, xlon, ter
- real, dimension(122000) :: rwrk
- real, pointer, dimension(:,:,:) :: real_array
- character (len=3) :: memorder
- character (len=25) :: crotang
- character (len=25) :: units
- character (len=46) :: desc
- character (len=128) :: init_date, cname, stagger, cunits, cdesc, title, startdate, grid_type, mminlu
- character (len=128), dimension(3) :: dimnames
- call getarg(1,crotang)
- read (crotang,'(f)') rotang
- write(6,*) 'Plotting with rotation angle ',rotang
- call opngks
- call gopwk(13, 41, 3)
- call gscr(1, 0, 1.00, 1.00, 1.00)
- call gscr(1, 1, 0.00, 0.00, 0.00)
- call gscr(1, 2, 0.25, 0.25, 0.25)
- call gscr(1, 3, 1.00, 1.00, 0.50)
- call gscr(1, 4, 0.50, 1.00, 0.50)
- call gscr(1, 5, 1.00, 1.00, 0.00)
- call gscr(1, 6, 1.00, 1.00, 0.00)
- call gscr(1, 7, 0.50, 1.00, 0.50)
- call gscr(1, 8, 1.00, 1.00, 0.50)
- call gscr(1, 9, 0.50, 1.00, 0.50)
- call gscr(1,10, 0.50, 1.00, 0.50)
- call gscr(1,11, 1.00, 1.00, 0.50)
- call gscr(1,12, 0.00, 1.00, 0.00)
- call gscr(1,13, 0.00, 0.50, 0.00)
- call gscr(1,14, 0.00, 1.00, 0.00)
- call gscr(1,15, 0.00, 0.50, 0.00)
- call gscr(1,16, 0.00, 1.00, 0.00)
- call gscr(1,17, 0.50, 0.50, 1.00)
- call gscr(1,18, 0.00, 1.00, 0.00)
- call gscr(1,19, 0.00, 1.00, 0.00)
- call gscr(1,20, 0.75, 0.75, 0.75)
- call gscr(1,21, 0.75, 0.75, 0.75)
- call gscr(1,22, 0.00, 0.50, 0.00)
- call gscr(1,23, 0.75, 0.75, 0.75)
- call gscr(1,24, 0.75, 0.75, 0.75)
- call gscr(1,25, 1.00, 1.00, 1.00)
- start_r = 0.00
- end_r = 0.50
- start_g = 1.00
- end_g = 0.25
- start_b = 0.00
- end_b = 0.00
- do i=26,76
- 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))
- end do
- start_r = 0.50
- end_r = 1.00
- start_g = 0.25
- end_g = 1.00
- start_b = 0.00
- end_b = 1.00
- do i=77,126
- 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))
- end do
- start_r = 0.80
- end_r = 1.00
- start_g = 0.80
- end_g = 1.00
- start_b = 0.80
- end_b = 1.00
- do i=127,176
- 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))
- end do
- call get_namelist_params()
- do n=1,max_dom
- call input_init(n, istatus)
- if (istatus /= 0) then
- write(6,*) ' '
- write(6,*) 'Error: Could not open domain01 file.'
- write(6,*) ' '
- stop
- end if
- call read_global_attrs(title, init_date, grid_type, dyn_opt, &
- west_east_dim, south_north_dim, bottom_top_dim, &
- we_patch_s, we_patch_e, we_patch_s_stag, we_patch_e_stag, &
- sn_patch_s, sn_patch_e, sn_patch_s_stag, sn_patch_e_stag, &
- map_proj, mminlu, num_land_cat, is_water, &
- is_ice, is_urban, isoilwater, grid_id, parent_id, i_parent_start, j_parent_start, &
- i_parent_end, j_parent_end, dx, dy, cen_lat, moad_cen_lat, cen_lon, &
- stand_lon, truelat1, truelat2, pole_lat, pole_lon, parent_grid_ratio, &
- corner_lats, corner_lons)
- istatus = 0
- do while (istatus == 0)
- call read_next_field(start_mem_i, end_mem_i, start_mem_j, end_mem_j, &
- start_mem_k, end_mem_k, cname, cunits, cdesc, &
- memorder, stagger, dimnames, real_array, istatus)
- if (istatus == 0) then
- if (index(cname, 'XLAT_M') /= 0) then
- nx = end_mem_i - start_mem_i + 1
- ny = end_mem_j - start_mem_j + 1
- allocate(xlat(nx,ny))
- xlat = real_array(:,:,1)
- else if (index(cname, 'XLONG_M') /= 0) then
- nx = end_mem_i - start_mem_i + 1
- ny = end_mem_j - start_mem_j + 1
- allocate(xlon(nx,ny))
- xlon = real_array(:,:,1)
- else if (index(cname, 'LU_INDEX') /= 0) then
- nx = end_mem_i - start_mem_i + 1
- ny = end_mem_j - start_mem_j + 1
- allocate(lu(nx,ny))
- lu = nint(real_array(:,:,1))
- end if
- end if
- end do
- call input_close()
- ll_lat = xlat(1,1)
- ll_lon = xlon(1,1)
- ur_lat = xlat(nx,ny)
- ur_lon = xlon(nx,ny)
- ! if (ur_lon < 0.) ur_lon = ur_lon + 360.0
- if (n == 1) then
- left = 0.0
- right = 1.0
- bottom = 0.0
- top = 1.0
- call mappos(left,right,bottom,top)
- call mapstc('OU','CO')
- call maproj('CE', cen_lat, cen_lon, rotang)
- ! call maproj('LC', truelat1, stand_lon, truelat2)
- ! call maproj('ST', cen_lat, cen_lon, stand_lon)
- call mapset('CO', ll_lat, ll_lon, ur_lat, ur_lon)
- call mapint()
- end if
- call mpsetr('GR', 10.0)
- call maptrn(ll_lat, ll_lon, left, bottom)
- call maptrn(ur_lat, ur_lon, right, top)
- width = 1.02*(right-left)/real(nx)
- height = 1.02*(top-bottom)/real(ny)
- do j=1,ny
- do i=1,nx
- call map_square(xlat(i,j), xlon(i,j), width, height, lu(i,j)+1)
- end do
- end do
- if (n > 1) then
- call gsplci(0)
- call lined(left-width/2.,bottom-height/2.,left-width/2.,top+height/2.)
- call lined(left-width/2.,top+height/2.,right+width/2.,top+height/2.)
- call lined(right+width/2.,top+height/2.,right+width/2.,bottom-height/2.)
- call lined(right+width/2.,bottom-height/2.,left-width/2.,bottom-height/2.)
- call sflush()
- call gsplci(1)
- end if
- deallocate(xlat)
- deallocate(xlon)
- deallocate(lu)
- end do
- call mplndr('Earth..3',4)
- call arinam (iam,400000)
- call mapbla (iam)
- call arpram (iam,0,0,0)
- call mapgrm (iam,xcs,ycs,10000,iai,iag,10,ulpr)
- call frame()
- do n=1,max_dom
- call input_init(n, istatus)
- if (istatus /= 0) then
- write(6,*) ' '
- write(6,*) 'Error: Could not open domain01 file.'
- write(6,*) ' '
- stop
- end if
- call read_global_attrs(title, init_date, grid_type, dyn_opt, &
- west_east_dim, south_north_dim, bottom_top_dim, &
- we_patch_s, we_patch_e, we_patch_s_stag, we_patch_e_stag, &
- sn_patch_s, sn_patch_e, sn_patch_s_stag, sn_patch_e_stag, &
- map_proj, mminlu, num_land_cat, is_water, &
- is_ice, is_urban, isoilwater, grid_id, parent_id, i_parent_start, j_parent_start, &
- i_parent_end, j_parent_end, dx, dy, cen_lat, moad_cen_lat, cen_lon, &
- stand_lon, truelat1, truelat2, pole_lat, pole_lon, parent_grid_ratio, &
- corner_lats, corner_lons)
- istatus = 0
- do while (istatus == 0)
- call read_next_field(start_mem_i, end_mem_i, start_mem_j, end_mem_j, &
- start_mem_k, end_mem_k, cname, cunits, cdesc, &
- memorder, stagger, dimnames, real_array, istatus)
- if (istatus == 0) then
- if (index(cname, 'XLAT_M') /= 0) then
- nx = end_mem_i - start_mem_i + 1
- ny = end_mem_j - start_mem_j + 1
- allocate(xlat(nx,ny))
- xlat = real_array(:,:,1)
- else if (index(cname, 'XLONG_M') /= 0) then
- nx = end_mem_i - start_mem_i + 1
- ny = end_mem_j - start_mem_j + 1
- allocate(xlon(nx,ny))
- xlon = real_array(:,:,1)
- else if (index(cname, 'HGT_M') /= 0) then
- nx = end_mem_i - start_mem_i + 1
- ny = end_mem_j - start_mem_j + 1
- allocate(ter(nx,ny))
- ter = real_array(:,:,1)
- else if (index(cname, 'LU_INDEX') /= 0) then
- nx = end_mem_i - start_mem_i + 1
- ny = end_mem_j - start_mem_j + 1
- allocate(lu(nx,ny))
- lu = nint(real_array(:,:,1))
- end if
- end if
- end do
- call input_close()
- ll_lat = xlat(1,1)
- ll_lon = xlon(1,1)
- ur_lat = xlat(nx,ny)
- ur_lon = xlon(nx,ny)
- if (n == 1) then
- left = 0.0
- right = 1.0
- bottom = 0.0
- top = 1.0
- call mappos(left,right,bottom,top)
- call mapstc('OU','CO')
- call maproj('CE', cen_lat, cen_lon, rotang)
- ! call maproj('LC', truelat1, stand_lon, truelat2)
- ! call maproj('ST', cen_lat, cen_lon, stand_lon)
- call mapset('CO', ll_lat, ll_lon, ur_lat, ur_lon)
- call mapint()
- maxter = -10000.
- minter = 10000.
- do j=1,ny
- do i=1,nx
- if (ter(i,j) > maxter) maxter = ter(i,j)
- if (ter(i,j) < minter) minter = ter(i,j)
- end do
- end do
- ! maxter = 3348.42
- end if
- call maptrn(ll_lat, ll_lon, left, bottom)
- call maptrn(ur_lat, ur_lon, right, top)
- width = 1.02*(right-left)/real(nx)
- height = 1.02*(top-bottom)/real(ny)
- do j=1,ny
- do i=1,nx
- if (lu(i,j) == 16) then
- ter(i,j) = ((ter(i,j)-minter) * 99.)/(maxter-minter) + 26.
- call map_square(xlat(i,j), xlon(i,j), width, height, 17)
- else if (lu(i,j) == 1) then
- ter(i,j) = ((ter(i,j)-minter) * 99.)/(maxter-minter) + 26.
- call map_square(xlat(i,j), xlon(i,j), width, height, 2)
- else if (lu(i,j) == 24) then
- ter(i,j) = ((ter(i,j)-minter) * 50.)/(3500.0-minter) + 127.
- call map_square(xlat(i,j), xlon(i,j), width, height, nint(ter(i,j)))
- else
- ter(i,j) = ((ter(i,j)-minter) * 99.)/(maxter-minter) + 26.
- call map_square(xlat(i,j), xlon(i,j), width, height, nint(ter(i,j)))
- end if
- end do
- end do
- if (n > 1) then
- call gsplci(0)
- call lined(left-width/2.,bottom-height/2.,left-width/2.,top+height/2.)
- call lined(left-width/2.,top+height/2.,right+width/2.,top+height/2.)
- call lined(right+width/2.,top+height/2.,right+width/2.,bottom-height/2.)
- call lined(right+width/2.,bottom-height/2.,left-width/2.,bottom-height/2.)
- call sflush()
- call gsplci(1)
- end if
- deallocate(xlat)
- deallocate(xlon)
- deallocate(ter)
- deallocate(lu)
- end do
- call mplndr('Earth..3',4)
- call arinam (iam,400000)
- call mapbla (iam)
- call arpram (iam,0,0,0)
- call mapgrm (iam,xcs,ycs,10000,iai,iag,10,ulpr)
- call gclwk(13)
- call clsgks
- stop
- end program plotgrid
- subroutine map_square(rlat, rlon, width, height, colr)
- implicit none
- ! Arguments
- real :: rlat, rlon, width, height
- integer :: colr
- ! Local variables
- real :: u, v
- real, dimension(4) :: xra, yra
- real, dimension(2000) :: dst
- integer, dimension(3000) :: ind
- call maptrn(rlat, rlon, u, v)
- xra(1) = u-(width/2.)
- xra(2) = u+(width/2.)
- xra(3) = u+(width/2.)
- xra(4) = u-(width/2.)
- yra(1) = v-(height/2.)
- yra(2) = v-(height/2.)
- yra(3) = v+(height/2.)
- yra(4) = v+(height/2.)
- call sfsgfa(xra, yra, 4, dst, 2000, ind, 3000, colr)
- end subroutine map_square
- subroutine ulpr(xcs,ycs,ncs,iai,iag,nai)
- implicit none
- integer, external :: mapaci
- integer :: ncs, nai
- integer, dimension(nai) :: iai, iag
- real, dimension(ncs) :: xcs, ycs
- integer :: itm
- if (iai(1) >= 0 .and.iai(2) >= 0) then
- itm = max0(iai(1),iai(2))
- if (mapaci(itm) == 1) then
- if (ncs.gt.150) print * , 'ulpr - ncs too big - ',ncs
- call gpl(ncs,xcs,ycs)
- end if
- end if
- end subroutine ulpr