program plotgrids use map_utils implicit none ! Parameters integer, parameter :: MAX_DOMAINS = 21 ! Variables integer :: iproj_type, n_domains, io_form_output, dyn_opt integer :: i, j, max_dom, funit, io_form_geogrid integer :: interval_seconds integer, dimension(MAX_DOMAINS) :: parent_grid_ratio, parent_id, ixdim, jydim integer, dimension(MAX_DOMAINS) :: i_parent_start, j_parent_start, & s_we, e_we, s_sn, e_sn, & start_year, start_month, start_day, start_hour, & end_year, end_month, end_day, end_hour real :: known_lat, known_lon, stand_lon, truelat1, truelat2, known_x, known_y, & dxkm, dykm, phi, lambda, ref_lat, ref_lon, ref_x, ref_y real :: dx, dy real :: ri, rj, rlats, rlons, rlate, rlone real :: polat , rot real :: rparent_gridpts real :: xa,xb,ya,yb,xxa,xxy,yya,yyb real :: xs, xe, ys, ye integer :: jproj, jgrid, jlts, iusout, idot, ier integer :: ltype , idom real, dimension(MAX_DOMAINS) :: parent_ll_x, parent_ll_y, parent_ur_x, parent_ur_y, & dom_cen_lat, dom_cen_lon character (len=128) :: geog_data_path, opt_output_from_geogrid_path, opt_geogrid_tbl_path character (len=128), dimension(MAX_DOMAINS) :: geog_data_res character (len=128) :: map_proj,ncep_proc_path, ncep_proc_prefix, & ncep_proc_domain_type character (len=128), dimension(MAX_DOMAINS) :: start_date, end_date character (len=3) :: wrf_core character (len=1) :: gridtype logical :: do_tiled_output integer :: debug_level logical :: is_used logical :: ncep_processing, do_gwd, just_last, use_igbp, ncep_proc_grib2 logical :: movable_nests type (proj_info) :: map_projection namelist /share/ wrf_core, max_dom, start_date, end_date, & start_year, end_year, start_month, end_month, & start_day, end_day, start_hour, end_hour, & interval_seconds, & io_form_geogrid, opt_output_from_geogrid_path, debug_level namelist /geogrid/ parent_id, parent_grid_ratio, & i_parent_start, j_parent_start, s_we, e_we, s_sn, e_sn, & map_proj, ref_x, ref_y, ref_lat, ref_lon, dom_cen_lat, & dom_cen_lon, & truelat1, truelat2, stand_lon, dx, dy, & geog_data_res, geog_data_path, opt_geogrid_tbl_path, & ncep_processing, ncep_proc_path, ncep_proc_prefix, & ncep_proc_domain_type, do_gwd, just_last, use_igbp, & ncep_proc_grib2,movable_nests ! Set defaults for namelist variables debug_level = 0 io_form_geogrid = 2 wrf_core = 'ARW' max_dom = 1 geog_data_path = 'NOT_SPECIFIED' ref_x = NAN ref_y = NAN ref_lat = NAN ref_lon = NAN dx = 10000. dy = 10000. map_proj = 'Lambert' truelat1 = NAN truelat2 = NAN stand_lon = NAN do i=1,MAX_DOMAINS geog_data_res(i) = 'default' parent_id(i) = 1 parent_grid_ratio(i) = INVALID s_we(i) = 1 e_we(i) = INVALID s_sn(i) = 1 e_sn(i) = INVALID start_year(i) = 0 start_month(i) = 0 start_day(i) = 0 start_hour(i) = 0 end_year(i) = 0 end_month(i) = 0 end_day(i) = 0 end_hour(i) = 0 start_date(i) = '0000-00-00_00:00:00' end_date(i) = '0000-00-00_00:00:00' end do opt_output_from_geogrid_path = './' opt_geogrid_tbl_path = 'geogrid/' interval_seconds = INVALID ! Read parameters from Fortran namelist do funit=10,100 inquire(unit=funit, opened=is_used) if (.not. is_used) exit end do open(funit,file='namelist.nps',status='old',form='formatted',err=1000) read(funit,share) read(funit,geogrid) close(funit) dxkm = dx dykm = dy known_lat = ref_lat known_lon = ref_lon known_x = ref_x known_y = ref_y ! Convert wrf_core to uppercase letters do i=1,3 if (ichar(wrf_core(i:i)) >= 97) wrf_core(i:i) = char(ichar(wrf_core(i:i))-32) end do ! Before doing anything else, we must have a valid grid type gridtype = ' ' if (wrf_core == 'ARW') then gridtype = 'C' dyn_opt = 2 else if (wrf_core == 'NMM') then gridtype = 'E' dyn_opt = 4 CALL plot_e_grid ( e_we(1)-1 , e_sn(1)-1 , ref_lat , -1. * ref_lon , dy , dx ) stop else if (wrf_core == 'NMB') then gridtype = 'B' dyn_opt = 4 CALL plot_b_grid ( e_we , e_sn , ref_lat , -1. * ref_lon , dy , dx , & debug_level,dom_cen_lat, dom_cen_lon, parent_grid_ratio, & parent_id, max_dom ) stop end if if (gridtype /= 'C' .and. gridtype /= 'E') then write(6,*) 'A valid wrf_core must be specified in the namelist. '// & 'Currently, only "ARW" and "NMM" are supported.' stop end if if (max_dom > MAX_DOMAINS) then write(6,*) 'In namelist, max_dom must be <= ',MAX_DOMAINS,'. To run with more'// & ' than ',MAX_DOMAINS,' domains, increase the MAX_DOMAINS parameter.' stop end if ! Every domain must have a valid parent id do i=2,max_dom if (parent_id(i) <= 0 .or. parent_id(i) >= i) then write(6,*) 'In namelist, the parent_id of domain ',i,' must be in '// & 'the range 1 to ',i-1 stop end if end do ! Convert map_proj to uppercase letters do i=1,len(map_proj) if (ichar(map_proj(i:i)) >= 97) map_proj(i:i) = char(ichar(map_proj(i:i))-32) end do ! Assign parameters to module variables if ((index(map_proj, 'LAMBERT') /= 0) .and. & (len_trim(map_proj) == len('LAMBERT'))) then iproj_type = PROJ_LC rot=truelat1 polat=truelat2 jproj = 3 else if ((index(map_proj, 'MERCATOR') /= 0) .and. & (len_trim(map_proj) == len('MERCATOR'))) then iproj_type = PROJ_MERC rot=0. polat=0. jproj = 9 else if ((index(map_proj, 'POLAR') /= 0) .and. & (len_trim(map_proj) == len('POLAR'))) then iproj_type = PROJ_PS rot=0. polat=SIGN(90., ref_lat) jproj = 1 else if ((index(map_proj, 'ROTATED_LL') /= 0) .and. & (len_trim(map_proj) == len('ROTATED_LL'))) then iproj_type = PROJ_ROTLL else write(6,*) 'In namelist, invalid map_proj specified. Valid '// & 'projections are "lambert", "mercator", "polar", '// & 'and "rotated_ll".' end if n_domains = max_dom do i=1,n_domains ixdim(i) = e_we(i) - s_we(i) + 1 jydim(i) = e_sn(i) - s_sn(i) + 1 end do if (gridtype == 'E') then phi = dykm*real(jydim(1)-1)/2. lambda = dxkm*real(ixdim(1)-1) end if ! If the user hasn't supplied a known_x and known_y, assume the center of domain 1 if (known_x == NAN) known_x = ixdim(1) / 2. if (known_y == NAN) known_y = jydim(1) / 2. ! Checks specific to C grid if (gridtype == 'C') then ! C grid does not support the rotated lat/lon projection if (iproj_type == PROJ_ROTLL) then write(6,*) 'Rotated lat/lon projection is not supported for the ARW core. '// & 'Valid projecitons are "lambert", "mercator", and "polar".' stop end if ! Check that nests have an acceptable number of grid points in each dimension do i=2,n_domains rparent_gridpts = real(ixdim(i)-1)/real(parent_grid_ratio(i)) if (floor(rparent_gridpts) /= ceiling(rparent_gridpts)) then write(6,*) 'For nest ',i,' (e_we-s_we+1) must be one greater than an '// & 'integer multiple of the parent_grid_ratio.' stop end if rparent_gridpts = real(jydim(i)-1)/real(parent_grid_ratio(i)) if (floor(rparent_gridpts) /= ceiling(rparent_gridpts)) then write(6,*) 'For nest ',i,' (e_sn-s_sn+1) must be one greater than an '// & 'integer multiple of the parent_grid_ratio.' stop end if end do end if do i=1,n_domains parent_ll_x(i) = real(i_parent_start(i)) parent_ll_y(i) = real(j_parent_start(i)) parent_ur_x(i) = real(i_parent_start(i))+real(ixdim(i))/real(parent_grid_ratio(i))-1. parent_ur_y(i) = real(j_parent_start(i))+real(jydim(i))/real(parent_grid_ratio(i))-1. end do call map_init(map_projection) call map_set(iproj_type, map_projection, & lat1=known_lat, & lon1=known_lon, & knowni=known_x, & knownj=known_y, & dx=dx, & stdlon=stand_lon, & truelat1=truelat1, & truelat2=truelat2, & ixdim=ixdim(1), & jydim=jydim(1)) call ij_to_latlon(map_projection, 1., 1., rlats, rlons) call ij_to_latlon(map_projection, real(e_we(1)), real(e_sn(1)) , rlate, rlone) call opngks ! Set some colors call gscr(1, 0, 1.00, 1.00, 1.00) call gscr(1, 1, 0.00, 0.00, 0.00) ! Do not grind them with details jgrid=10 jlts=-2 iusout=1 idot=0 call supmap(jproj,polat,stand_lon,rot,& rlats,rlons,rlate,rlone, & jlts,jgrid,iusout,idot,ier) call setusv('LW',1000) call perim(e_we(1)-1,1,e_sn(1)-1,1) call getset(xa,xb,ya,yb,xxa,xxy,yya,yyb,ltype) call set (xa,xb,ya,yb, & 1.,real(e_we(1)),1.,real(e_sn(1)),ltype) do idom = 2 , max_dom call getxy ( xs, xe, ys, ye, & idom , max_dom , & e_we , e_sn , & parent_id , parent_grid_ratio , & i_parent_start , j_parent_start ) call line ( xs , ys , xe , ys ) call line ( xe , ys , xe , ye ) call line ( xe , ye , xs , ye ) call line ( xs , ye , xs , ys ) end do call frame call clsgks stop 1000 write(6,*) 'Error opening namelist.wps' stop end program plotgrids subroutine getxy ( xs, xe, ys, ye, & dom_id , num_domains , & e_we , e_sn , & parent_id , parent_grid_ratio , & i_parent_start , j_parent_start ) implicit none integer , intent(in) :: dom_id integer , intent(in) :: num_domains integer , intent(in) , dimension(num_domains):: e_we , e_sn , & parent_id , parent_grid_ratio , & i_parent_start , j_parent_start real , intent(out) :: xs, xe, ys, ye ! local vars integer :: idom xs = 0. xe = e_we(dom_id) -1 ys = 0. ye = e_sn(dom_id) -1 idom = dom_id compute_xy : DO xs = (i_parent_start(idom) + xs -1 ) / & real(parent_grid_ratio(parent_id(idom))) xe = xe / real(parent_grid_ratio(idom)) ys = (j_parent_start(idom) + ys -1 ) / & real(parent_grid_ratio(parent_id(idom))) ye = ye / real(parent_grid_ratio(idom)) idom = parent_id(idom) if ( idom .EQ. 1 ) then exit compute_xy end if END DO compute_xy xs = xs + 1 xe = xs + xe ys = ys + 1 ye = ys + ye end subroutine getxy SUBROUTINE plot_b_grid ( im , jm , rlat0d , rlon0d , dphd , dlmd , debug_level, & dom_cen_lat, dom_cen_lon, parent_grid_ratio, parent_id, max_dom ) ! This routine generates a gmeta file of the area covered by an Arakawa b-grid. ! We assume that NCAR Graphics has not been called yet (and will be closed ! upon exit). The required input fields are as from the WPS namelist file. IMPLICIT NONE ! Input map parameters for E/B grid. INTEGER , INTENT(IN) :: im(21) , & ! number of H points in odd rows jm(21) , & ! number of rows max_dom , & ! total number of domains debug_level ! debug REAL , INTENT(IN) :: rlat0d , & ! latitude of grid center (degrees) rlon0d ! longitude of grid center (degrees, times -1) REAL , INTENT(IN) :: dphd , & ! angular distance between rows (degrees) dlmd ! angular distance between adjacent H and V points (degrees) REAL, INTENT(IN) :: dom_cen_lat(21),dom_cen_lon(21) INTEGER, INTENT(IN) :: parent_grid_ratio(21),parent_id(21) ! Some local vars REAL :: rlat1d , & rlon1d INTEGER :: ngpwe , & ngpsn , & ilowl , & jlowl INTEGER :: imt , imtjm REAL :: latlft,lonlft,latrgt,lonrgt imt=2*im(1)-1 imtjm=imt*jm(1) rlat1d=rlat0d rlon1d=rlon0d ngpwe=im(1) ngpsn=jm(1) write(0,*) 'ngpwe , ngpsn: ' , ngpwe , ngpsn ! Get lat and lon of left and right points. CALL corners_bgrid ( rlat1d,rlon1d,im(1),jm(1),rlat0d,rlon0d,dphd,dlmd,& ngpwe,ngpsn,ilowl,jlowl,latlft,lonlft,latrgt,lonrgt) write(0,*) 'latlft: ', latlft write(0,*) 'lonlft: ', lonlft write(0,*) 'latrgt: ', latrgt write(0,*) 'lonrgt: ', lonrgt ! With corner points, make map background. CALL mapbkg_bgrid ( im,jm,ilowl,jlowl,ngpwe,ngpsn,& rlat0d,rlon0d,latlft,lonlft,latrgt,lonrgt,& dlmd,dphd,debug_level, & dom_cen_lat, dom_cen_lon, parent_grid_ratio, parent_id, max_dom ) END SUBROUTINE plot_b_grid !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!! E GRID MAP INFO BELOW !!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE plot_e_grid ( im , jm , rlat0d , rlon0d , dphd , dlmd ) ! This routine generates a gmeta file of the area covered by an Arakawa e-grid. ! We assume that NCAR Graphics has not been called yet (and will be closed ! upon exit). The required input fields are as from the WPS namelist file. IMPLICIT NONE ! 15 April 2005 NCEP/EMC ! The Code and some instructions are provided by Tom BLACK to ! NCAR/DTC Meral Demirtas ! 4 May 2005 NCAR/DTC Meral DEMIRTAS ! - An include file (plot_inc) is added to get ! Domain size: IM,JM ! Central latitute and longnitute: RLAT0D,RLON0D ! Horizontal resolution: DPHD, DLMD ! Feb 2007 NCAR/MMM ! Turn into f90 ! Add implicit none ! Remove non-mapping portions ! Make part of WPS domain plotting utility ! Input map parameters for E grid. INTEGER , INTENT(IN) :: im , & ! number of H points in odd rows jm ! number of rows REAL , INTENT(IN) :: rlat0d , & ! latitude of grid center (degrees) rlon0d ! longitude of grid center (degrees, times -1) REAL , INTENT(IN) :: dphd , & ! angular distance between rows (degrees) dlmd ! angular distance between adjacent H and V points (degrees) ! Some local vars REAL :: rlat1d , & rlon1d INTEGER :: ngpwe , & ngpsn , & ilowl , & jlowl INTEGER :: imt , imtjm REAL :: latlft,lonlft,latrgt,lonrgt imt=2*im-1 imtjm=imt*jm rlat1d=rlat0d rlon1d=rlon0d ngpwe=2*im-1 ngpsn=jm ! Get lat and lon of left and right points. CALL corners ( rlat1d,rlon1d,im,jm,rlat0d,rlon0d,dphd,dlmd,& ngpwe,ngpsn,ilowl,jlowl,latlft,lonlft,latrgt,lonrgt) ! With corner points, make map background. CALL mapbkg_egrid ( imt,jm,ilowl,jlowl,ngpwe,ngpsn,& rlat0d,rlon0d,latlft,lonlft,latrgt,lonrgt,& dlmd,dphd) END SUBROUTINE plot_e_grid !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE corners_bgrid ( glatd,glond,im,jm,tph0d,tlm0d,dphd,dlmd,& ngpwe,ngpsn,ilowl,jlowl,glatl,glonl,glatr,glonr) IMPLICIT NONE REAL , INTENT(IN) :: glatd,glond,tph0d,tlm0d,dphd,dlmd,& glatl,glonl,glatr,glonr INTEGER , INTENT(IN) :: im,jm,ngpwe,ngpsn INTEGER , INTENT(OUT) :: ilowl,jlowl ! Local vars REAL , PARAMETER :: d2r = 1.74532925E-2 , r2D = 1./D2R REAL :: glat , glon , dph , dlm , tph0 , tlm0 REAL :: x , y , z , tlat , tlon , tlat1 , tlat2 , tlon1 , tlon2 REAL :: row , col REAL :: dlm1 , dlm2 , d1 , d2 , d3 , d4 , dmin INTEGER :: jmt , ii , jj , iuppr , juppr, imt INTEGER :: nrow , ncol jmt = jm/2+1 imt = im/2+1 ! Convert from geodetic to transformed coordinates (degrees). glat = glatd * d2r glon = glond * d2r dph = dphd * d2r dlm = dlmd * d2r tph0 = tph0d * d2r tlm0 = tlm0d * d2r write(0,*) 'glat, glon: ', glat, glon write(0,*) 'dph, dlm: ', dph, dlm write(0,*) 'tph0, tlm0: ', tph0, tlm0 x = COS(tph0) * COS(glat) * COS(glon-tlm0)+SIN(tph0) * SIN(glat) y = -COS(glat) * SIN(glon-tlm0) z = COS(tph0) * SIN(glat)-SIN(tph0) * COS(glat) * COS(glon-tlm0) write(0,*) 'x, y, z: ', x, y, z tlat = r2d * ATAN(z/SQRT(x*x + y*y)) tlon = r2d * ATAN(y/x) ! Find the real (non-integer) row and column of the input location on ! the filled e-grid. write(0,*) 'tlat, tlon: ', tlat, tlon row = tlat/dphd+jmt col = tlon/dlmd+imt nrow = INT(row) ncol = INT(col) write(0,*) 'ncol, nrow: ', ncol, nrow tlat = tlat * d2r tlon = tlon * d2r ! E2 E3 ! ! ! X ! E1 E4 tlat1 = (nrow-jmt) * dph tlat2 = tlat1+dph tlon1 = (ncol-imt) * dlm tlon2 = tlon1+dlm write(0,*) 'tlat1, tlat2: ', tlat1, tlat2 write(0,*) 'tlon1, tlon2: ', tlon1, tlon2 dlm1 = tlon-tlon1 dlm2 = tlon-tlon2 d1 = ACOS(COS(tlat) * COS(tlat1) * COS(dlm1)+SIN(tlat) * SIN(tlat1)) d2 = ACOS(COS(tlat) * COS(tlat2) * COS(dlm1)+SIN(tlat) * SIN(tlat2)) d3 = ACOS(COS(tlat) * COS(tlat2) * COS(dlm2)+SIN(tlat) * SIN(tlat2)) d4 = ACOS(COS(tlat) * COS(tlat1) * COS(dlm2)+SIN(tlat) * SIN(tlat1)) dmin = MIN(d1,d2,d3,d4) IF ( ABS(dmin-d1) .LT. 1.e-6 ) THEN ii = ncol jj = nrow ELSE IF ( ABS(dmin-d2) .LT. 1.e-6 ) THEN ii = ncol jj = nrow+1 ELSE IF ( ABS(dmin-d3) .LT. 1.e-6 ) THEN ii = ncol+1 jj = nrow+1 ELSE IF ( ABS(dmin-d4) .LT. 1.e-6 ) THEN ii = ncol+1 jj = nrow END IF ! write(0,*) 'ii,jj final: ', ii, jj ! Now find the i and j of the lower left corner of the desired grid ! region and of the upper right. ilowl = ii-ngpwe/2 jlowl = jj-ngpsn/2 iuppr = ii+ngpwe/2 juppr = jj+ngpsn/2 ! write(0,*) 'ilowl, jlowl: ', ilowl, jlowl ! write(0,*) 'iuppr, juppr: ', iuppr,juppr iuppr=im juppr=jm ! Find their geodetic coordinates. write(0,*) 'ilowl, jlowl, im, jm: ', ilowl, jlowl, im, jm write(0,*) 'tph0d, tlm0d, dphd, dlmd: ', tph0d, tlm0d, dphd, dlmd CALL b2t2g(ilowl,jlowl,im,jm,tph0d,tlm0d,dphd,dlmd,glatl,glonl) write(0,*) 'glatl,glonl: ', glatl,glonl CALL b2t2g(iuppr,juppr,im,jm,tph0d,tlm0d,dphd,dlmd,glatr,glonr) write(0,*) 'glatr,glonr: ', glatr,glonr ! CALL b2t2g(ilowl-10,jlowl-10,im,jm,tph0d,tlm0d,dphd,dlmd,glatl,glonl) ! write(0,*) 'glatl,glonl: ', glatl,glonl ! CALL b2t2g(iuppr+10,juppr+10,im,jm,tph0d,tlm0d,dphd,dlmd,glatr,glonr) ! write(0,*) 'glatr,glonr: ', glatr,glonr END SUBROUTINE corners_bgrid SUBROUTINE corners ( glatd,glond,im,jm,tph0d,tlm0d,dphd,dlmd,& ngpwe,ngpsn,ilowl,jlowl,glatl,glonl,glatr,glonr) IMPLICIT NONE REAL , INTENT(IN) :: glatd,glond,tph0d,tlm0d,dphd,dlmd,& glatl,glonl,glatr,glonr INTEGER , INTENT(IN) :: im,jm,ngpwe,ngpsn INTEGER , INTENT(OUT) :: ilowl,jlowl ! Local vars REAL , PARAMETER :: d2r = 1.74532925E-2 , r2D = 1./D2R REAL :: glat , glon , dph , dlm , tph0 , tlm0 REAL :: x , y , z , tlat , tlon , tlat1 , tlat2 , tlon1 , tlon2 REAL :: row , col REAL :: dlm1 , dlm2 , d1 , d2 , d3 , d4 , dmin INTEGER :: jmt , ii , jj , iuppr , juppr INTEGER :: nrow , ncol jmt = jm/2+1 ! Convert from geodetic to transformed coordinates (degrees). glat = glatd * d2r glon = glond * d2r dph = dphd * d2r dlm = dlmd * d2r tph0 = tph0d * d2r tlm0 = tlm0d * d2r x = COS(tph0) * COS(glat) * COS(glon-tlm0)+SIN(tph0) * SIN(glat) y = -COS(glat) * SIN(glon-tlm0) z = COS(tph0) * SIN(glat)-SIN(tph0) * COS(glat) * COS(glon-tlm0) tlat = r2d * ATAN(z/SQRT(x*x + y*y)) tlon = r2d * ATAN(y/x) ! Find the real (non-integer) row and column of the input location on ! the filled e-grid. row = tlat/dphd+jmt col = tlon/dlmd+im nrow = INT(row) ncol = INT(col) tlat = tlat * d2r tlon = tlon * d2r ! E2 E3 ! ! ! X ! E1 E4 tlat1 = (nrow-jmt) * dph tlat2 = tlat1+dph tlon1 = (ncol-im) * dlm tlon2 = tlon1+dlm dlm1 = tlon-tlon1 dlm2 = tlon-tlon2 d1 = ACOS(COS(tlat) * COS(tlat1) * COS(dlm1)+SIN(tlat) * SIN(tlat1)) d2 = ACOS(COS(tlat) * COS(tlat2) * COS(dlm1)+SIN(tlat) * SIN(tlat2)) d3 = ACOS(COS(tlat) * COS(tlat2) * COS(dlm2)+SIN(tlat) * SIN(tlat2)) d4 = ACOS(COS(tlat) * COS(tlat1) * COS(dlm2)+SIN(tlat) * SIN(tlat1)) dmin = MIN(d1,d2,d3,d4) IF ( ABS(dmin-d1) .LT. 1.e-6 ) THEN ii = ncol jj = nrow ELSE IF ( ABS(dmin-d2) .LT. 1.e-6 ) THEN ii = ncol jj = nrow+1 ELSE IF ( ABS(dmin-d3) .LT. 1.e-6 ) THEN ii = ncol+1 jj = nrow+1 ELSE IF ( ABS(dmin-d4) .LT. 1.e-6 ) THEN ii = ncol+1 jj = nrow END IF ! Now find the i and j of the lower left corner of the desired grid ! region and of the upper right. ilowl = ii-ngpwe/2 jlowl = jj-ngpsn/2 iuppr = ii+ngpwe/2 juppr = jj+ngpsn/2 ! Find their geodetic coordinates. CALL e2t2g(ilowl,jlowl,im,jm,tph0d,tlm0d,dphd,dlmd,glatl,glonl) CALL e2t2g(iuppr,juppr,im,jm,tph0d,tlm0d,dphd,dlmd,glatr,glonr) END SUBROUTINE corners !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE mapbkg_bgrid ( im,jm,ilowl,jlowl,ngpwe,ngpsn,& rlat0d,rlon0d,glatl,glonl,glatr,glonr,& dlmd,dphd, debug_level, & dom_cen_lat, dom_cen_lon, parent_grid_ratio, parent_id, max_dom ) ! IMPLICIT NONE ! Some local vars real, dimension(21) :: dom_cen_lat, dom_cen_lon, dx, dy integer, dimension(21) :: parent_grid_ratio, parent_id INTEGER :: max_dom, I, im(21),jm(21),debug_level CHARACTER (LEN=97) :: string REAL:: wbd_guess, sbd_guess, wbd_loc, sbd_loc REAL:: xs, xe, ys, ye REAL:: wbd(21),sbd(21) ! Yet more center lon messing around, hoo boy. ! clonx=180.-rlon0d clonx=-rlon0d ! Open up NCAR Graphics CALL opngks CALL gopwk(8,9,3) CALL gsclip(0) write(0,*) 'clonx: ', clonx ! Make the background white, and the foreground black. CALL gscr ( 1 , 0 , 1., 1., 1. ) CALL gscr ( 1 , 1 , 0., 0., 0. ) ! Line width twice default thickness. CALL setusv('LW',1800) ! Make map outline a solid line, not dots. CALL mapsti('MV',8) CALL mapsti('DO',0) ! Map outlines are political and states. CALL mapstc('OU','PS') ! Cylindrical equidistant. CALL maproj('CE',rlat0d,clonx,0.) ! Specify corner points. CALL mapset('CO',glatl,glonl,glatr,glonr) ! Lat lon lines every 5 degrees. CALL mapsti('GR',5) ! Map takes up this much real estate. CALL mappos( 0.03 , 0.97 , 0.03 , 0.97 ) ! Initialize and draw map. CALL mapint ! CALL mapdrw ! Add approx grid point tick marks ! CALL setusv('LW',1000) ! CALL perim(im,1,jm,1) !!! these x y values are in degrees of latitude longitude if (debug_level > 0) then write(0,*) 'start new stuff' write(0,*) 'max_dom: ', max_dom endif dx(1)=dlmd dy(1)=dphd if (debug_level > 0) then write(0,*) 'dx(1), dy(1): ', dx(1), dy(1) write(0,*) 'im(1), jm(1): ', im(1), jm(1) endif wbd(1)=-((im(1)-1)/2.)*dlmd sbd(1)=-((jm(1)-1)/2.)*dphd if (debug_level > 0) then write(0,*) 'im(1), wbd(1): ', im(1), wbd(1) write(0,*) 'jm(1), sbd(1): ', jm(1), sbd(1) write(0,*) 'dom_cen_lat/lon(1): ', dom_cen_lat(1), dom_cen_lon(1) write(0,*) 'parent_id,ratio(1): ', parent_id(1), parent_grid_ratio(1) endif write(0,*) 'loop for 2 to max_dom: ', max_dom do I=2,max_dom if (debug_level > 0) then write(0,*) ' ' write(0,*) '-----------------------------------' write(0,*) ' ' write(0,*) 'dom_cen_lat(i): ', dom_cen_lat(i) write(0,*) 'dom_cen_lon(i): ', dom_cen_lon(i) write(0,*) 'parent_id,ratio(i): ', parent_id(i), parent_grid_ratio(i) endif call tll (dom_cen_lon(i),dom_cen_lat(i),xcenter, ycenter, dom_cen_lat(1),dom_cen_lon(1)) if (debug_level > 0) then write(0,*) 'im,jm(i): ', im(i),jm(i) write(0,*) 'xcenter, ycenter: ', xcenter, ycenter write(0,*) 'dlmd, dphd: ', dlmd, dphd endif wbd_guess = xcenter - & (im(i)-1)*0.5*dx(parent_id(i))/parent_grid_ratio(i) sbd_guess = ycenter - & (jm(i)-1)*0.5*dy(parent_id(i))/parent_grid_ratio(i) if (debug_level > 0) then write(0,*) 'wbd_guess, sbd_guess: ', wbd_guess, sbd_guess write(0,*) 'i_parent_start parens: ', & ( wbd_guess - wbd(parent_id(i)))/dx(parent_id(i)) write(0,*) 'j_parent_start parens: ', & ( sbd_guess - sbd(parent_id(i)))/dy(parent_id(i)) endif i_parent_start = nint (1 + ( wbd_guess - wbd(parent_id(i)))/dx(parent_id(i)) ) j_parent_start = nint (1 + ( sbd_guess - sbd(parent_id(i)))/dy(parent_id(i)) ) if (debug_level > 0) then write(0,*) 'wbd(parent_id(i)), sbd(parent_id(i)): ', & wbd(parent_id(i)), sbd(parent_id(i)) write(0,*) 'dx,dy(parent_id(i)): ', dx(parent_id(i)), dy(parent_id(i)) write(0,*) 'i/j_parent_start: ', i_parent_start, j_parent_start endif i_parent_end = i_parent_start + (im(i)-1)/parent_grid_ratio(i) j_parent_end = j_parent_start + (jm(i)-1)/parent_grid_ratio(i) if (debug_level > 0) then write(0,*) 'i/j_parent_end: ', i_parent_end, j_parent_end endif wbd(i)= wbd(parent_id(i)) + (i_parent_start-1)*dx(parent_id(i)) sbd(i)= sbd(parent_id(i)) + (j_parent_start-1)*dy(parent_id(i)) dx(i)=dx(parent_id(i))/parent_grid_ratio(i) dy(i)=dy(parent_id(i))/parent_grid_ratio(i) if (debug_level > 0) then write(0,*) 'dx(i), dy(i): ', dx(i), dy(i) endif wbd_loc=-((im(i)-1)/2)*dx(i) sbd_loc=-((jm(i)-1)/2)*dy(i) if (debug_level > 0) then write(0,*) 'wbd_loc, sbd_loc: ', wbd_loc, sbd_loc endif xs=wbd(i) xe=wbd(i)-2*wbd_loc ys=sbd(i) ye=sbd(i)-2*sbd_loc if (debug_level > 0) then write(0,*) 'xs, xe: ', wbd(i), wbd(i)-2*wbd_loc write(0,*) 'ys, ye: ', sbd(i), sbd(i)-2*sbd_loc write(0,*) '+ + + + + + + +' endif write(0,*) 'here a' call line( xs , ys , xe , ys ) call line( xe , ys , xe , ye ) call line( xe , ye , xs , ye ) call line( xs , ye , xs , ys ) enddo write(0,*) 'here b' CALL mapdrw ! Put on a quicky description. WRITE ( string , FMT = '("B-GRID E_WE = ",I4,", E_SN = ",I4 , & &", DX = ",F6.4,", DY = ",F6.4 , & &", REF_LAT = ",F8.3,", REF_LON = ",F8.3)') & im(1),jm(1),dlmd,dphd,rlat0d,-1.*rlon0d CALL getset(xa,xb,ya,yb,xxa,xxy,yya,yyb,ltype) ! CALL set (xa,xb,ya,yb,0.,1.,0.,1.,lytpe) ! CALL gsclip(0) ! CALL pchiqu ( 0.5 , -5. , string , 8. , 0., 0. ) CALL pchiqu ( xxa , yya - (yyb-yya)/20., string , 8. , 0., -1. ) CALL frame ! Close workstation and NCAR Grpahics. CALL gclwk(8) CALL clsgks END SUBROUTINE mapbkg_bgrid SUBROUTINE mapbkg_egrid ( imt,jm,ilowl,jlowl,ngpwe,ngpsn,& rlat0d,rlon0d,glatl,glonl,glatr,glonr,& dlmd,dphd) ! IMPLICIT NONE ! Some local vars CHARACTER (LEN=97) :: string ! Yet more center lon messing around, hoo boy. ! clonx=180.-rlon0d clonx=-rlon0d ! Open up NCAR Graphics CALL opngks CALL gopwk(8,9,3) CALL gsclip(0) ! Make the background white, and the foreground black. CALL gscr ( 1 , 0 , 1., 1., 1. ) CALL gscr ( 1 , 1 , 0., 0., 0. ) ! Line width twice default thickness. CALL setusv('LW',2000) ! Make map outline a solid line, not dots. CALL mapsti('MV',8) CALL mapsti('DO',0) ! Map outlines are political and states. CALL mapstc('OU','PS') ! Cylindrical equidistant. CALL maproj('CE',rlat0d,clonx,0.) ! Specify corner points. CALL mapset('CO',glatl,glonl,glatr,glonr) ! Lat lon lines every 5 degrees. CALL mapsti('GR',5) ! Map takes up this much real estate. CALL mappos( 0.05 , 0.95 , 0.05 , 0.95 ) ! Initialize and draw map. CALL mapint CALL mapdrw ! Add approx grid point tick marks CALL setusv('LW',1000) CALL perim(((imt+3)/2)-1,1,jm-1,1) ! Put on a quicky description. WRITE ( string , FMT = '("E-GRID E_WE = ",I4,", E_SN = ",I4 , & &", DX = ",F6.4,", DY = ",F6.4 , & &", REF_LAT = ",F8.3,", REF_LON = ",F8.3)') & (imt+3)/2,jm+1,dlmd,dphd,rlat0d,-1.*rlon0d CALL getset(xa,xb,ya,yb,xxa,xxy,yya,yyb,ltype) ! CALL set (xa,xb,ya,yb,0.,1.,0.,1.,lytpe) ! CALL gsclip(0) ! CALL pchiqu ( 0.5 , -5. , string , 8. , 0., 0. ) CALL pchiqu ( xxa , yya - (yyb-yya)/20., string , 8. , 0., -1. ) CALL frame ! Close workstation and NCAR Grpahics. CALL gclwk(8) CALL clsgks END SUBROUTINE mapbkg_egrid !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE e2t2g ( ncol,nrow,im,jm,tph0d,tlm0d,dphd,dlmd,glatd,glond) ! IMPLICIT NONE REAL , PARAMETER :: D2R=1.74532925E-2 , R2D=1./D2R DPH=DPHD*D2R DLM=DLMD*D2R TPH0=TPH0D*D2R TLM0=TLM0D*D2R !*** FIND THE TRANSFORMED LAT (POSITIVE NORTH) AND LON (POSITIVE EAST) TLATD=(NROW-(JM+1)/2)*DPHD TLOND=(NCOL-IM)*DLMD !*** NOW CONVERT TO GEODETIC LAT (POSITIVE NORTH) AND LON (POSITIVE EAST) TLATR=TLATD*D2R TLONR=TLOND*D2R ARG1=SIN(TLATR)*COS(TPH0)+COS(TLATR)*SIN(TPH0)*COS(TLONR) GLATR=ASIN(ARG1) GLATD=GLATR*R2D ARG2=COS(TLATR)*COS(TLONR)/(COS(GLATR)*COS(TPH0))- & TAN(GLATR)*TAN(TPH0) IF(ABS(ARG2).GT.1.)ARG2=ABS(ARG2)/ARG2 FCTR=1. IF(TLOND.GT.0.)FCTR=-1. GLOND=TLM0D+FCTR*ACOS(ARG2)*R2D GLOND=-GLOND END SUBROUTINE e2t2g SUBROUTINE b2t2g ( ncol,nrow,im,jm,tph0d,tlm0d,dphd,dlmd,glatd,glond) ! IMPLICIT NONE INTEGER, INTENT(IN) :: NCOL, NROW, IM, JM REAL, INTENT(IN):: TPH0D, TLM0D, DPHD, DLMD REAL, INTENT(OUT):: GLATD, GLOND REAL , PARAMETER :: D2R=1.74532925E-2 , R2D=1./D2R DPH=DPHD*D2R DLM=DLMD*D2R TPH0=TPH0D*D2R TLM0=TLM0D*D2R !*** FIND THE TRANSFORMED LAT (POSITIVE NORTH) AND LON (POSITIVE EAST) TLATD=(NROW-(JM+1)/2)*DPHD TLOND=(NCOL-(IM+1)/2)*DLMD ! write(0,*) 'TLATD, TLOND: ', TLATD, TLOND !*** NOW CONVERT TO GEODETIC LAT (POSITIVE NORTH) AND LON (POSITIVE EAST) TLATR=TLATD*D2R TLONR=TLOND*D2R ARG1=SIN(TLATR)*COS(TPH0)+COS(TLATR)*SIN(TPH0)*COS(TLONR) GLATR=ASIN(ARG1) GLATD=GLATR*R2D write(0,*) 'GLATD into ARG2 calc: ', GLATD ARG2=COS(TLATR)*COS(TLONR)/(COS(GLATR)*COS(TPH0))- & TAN(GLATR)*TAN(TPH0) write(0,*) 'ARG2: ', arg2 IF(ABS(ARG2).GT.1.)ARG2=ABS(ARG2)/ARG2 write(0,*) 'now ARG2: ', arg2 FCTR=1. IF(TLOND.GT.0.)FCTR=-1. write(0,*) 'TLM0D, FCTR, ACOS(ARG2)*R2D: ', TLM0D, FCTR, ACOS(ARG2)*R2D GLOND=TLM0D+FCTR*ACOS(ARG2)*R2D GLOND=-GLOND write(0,*) 'exit with GLATD, GLOND: ', GLATD, GLOND END SUBROUTINE b2t2g subroutine tll(almd,aphd,tlmd,tphd,tph0d,tlm0d) !------------------------------------------------------------------------------- real, intent(in) :: almd, aphd real, intent(out) :: tlmd, tphd real, intent(in) :: tph0d, tlm0d !------------------------------------------------------------------------------- real, parameter :: pi=3.141592654 real, parameter :: dtr=pi/180.0 ! real :: tph0, ctph0, stph0, relm, srlm, crlm real :: aph, sph, cph, cc, anum, denom !------------------------------------------------------------------------------- ! if (tlm0d==0.0.and.tph0d==0.0) then tlmd=almd tphd=aphd else tph0=tph0d*dtr ctph0=cos(tph0) stph0=sin(tph0) ! relm=(almd-tlm0d)*dtr srlm=sin(relm) crlm=cos(relm) aph=aphd*dtr sph=sin(aph) cph=cos(aph) cc=cph*crlm anum=cph*srlm denom=ctph0*cc+stph0*sph ! tlmd=atan2(anum,denom)/dtr tphd=asin(ctph0*sph-stph0*cc)/dtr end if ! return ! end subroutine tll