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