model_grid.F90

Go to the documentation of this file.

 module model_grid

 use esmf
 use esmf_logpublicmod

 implicit none

 private

 character(len=5), allocatable, public  :: tiles_target_grid(:)
                                           
 character(len=50), public              :: input_grid_type = "latlon"
                                           

 ! Made lsoil_target non-parameter to allow for RAP land surface initiation
 integer, public                        :: lsoil_target = 4 ! # soil layers
                                           
 integer, public                        :: i_input
                                           
 integer, public                        :: j_input
                                           
 integer, public                        :: ip1_input
                                           
 integer, public                        :: jp1_input
                                           
 integer, public                        :: i_target
                                           
 integer, public                        :: j_target
                                           
 integer, public                        :: ip1_target
                                           
 integer, public                        :: jp1_target
                                           
 integer, public                        :: num_tiles_input_grid
                                           
 integer, public                        :: num_tiles_target_grid
                                           

 type(esmf_grid),  public               :: input_grid
                                           
 type(esmf_grid),  public               :: target_grid
                                           

 type(esmf_field),  public              :: latitude_input_grid
                                           
 type(esmf_field),  public              :: longitude_input_grid
                                           
 type(esmf_field),  public              :: latitude_s_input_grid
                                           
 type(esmf_field),  public              :: longitude_s_input_grid
                                           
 type(esmf_field),  public              :: latitude_w_input_grid
                                           
 type(esmf_field),  public              :: longitude_w_input_grid
                                           

 type(esmf_field),  public              :: landmask_target_grid
                                           
 type(esmf_field),  public              :: latitude_target_grid
                                           
 type(esmf_field),  public              :: latitude_s_target_grid
                                           
 type(esmf_field),  public              :: latitude_w_target_grid
                                           
 type(esmf_field),  public              :: longitude_target_grid
                                           
 type(esmf_field),  public              :: longitude_s_target_grid
                                           
 type(esmf_field),  public              :: longitude_w_target_grid
                                           
 type(esmf_field),  public              :: seamask_target_grid
                                           
 type(esmf_field),  public              :: terrain_target_grid
                                           

 public :: define_target_grid
 public :: define_input_grid
 public :: cleanup_input_target_grid_data

 contains

 subroutine define_input_grid(localpet, npets)

 use program_setup, only       : input_type

 implicit none

 integer, intent(in)          :: localpet, npets

 if (trim(input_type) == "gaussian_nemsio" .or. &
     trim(input_type) == "gfs_gaussian_nemsio" .or. &
     trim(input_type) == "gfs_sigio" .or. &
     trim(input_type) == "gaussian_netcdf") then
   call define_input_grid_gaussian(localpet, npets)
 elseif (trim(input_type) == "grib2") then
   call define_input_grid_grib2(localpet,npets)
 else
   call define_input_grid_mosaic(localpet, npets)
 endif

 end subroutine define_input_grid

 subroutine define_input_grid_gaussian(localpet, npets)

 use nemsio_module

 use program_setup, only       : data_dir_input_grid, &
                                 atm_files_input_grid, &
                                 sfc_files_input_grid, &
                                 input_type, &
                                 convert_atm, convert_sfc

 use sfcio_module
 use sigio_module
 use netcdf

 implicit none

 integer, intent(in)              :: localpet, npets

 character(len=250)               :: the_file

 integer                          :: i, j, rc, clb(2), cub(2), ncid, id_grid
 integer(sfcio_intkind)           :: rc2
 integer(sigio_intkind)           :: rc3

 real(esmf_kind_r8), allocatable  :: latitude(:,:)
 real(esmf_kind_r8), allocatable  :: longitude(:,:)
 real(esmf_kind_r8), pointer      :: lat_src_ptr(:,:)
 real(esmf_kind_r8), pointer      :: lon_src_ptr(:,:)
 real(esmf_kind_r8), pointer      :: lat_corner_src_ptr(:,:)
 real(esmf_kind_r8), pointer      :: lon_corner_src_ptr(:,:)
 real(esmf_kind_r8)               :: deltalon
 real(esmf_kind_r8), allocatable  :: slat(:), wlat(:)

 type(nemsio_gfile)               :: gfile
 type(esmf_polekind_flag)         :: polekindflag(2)
 type(sfcio_head)                 :: sfchead
 type(sigio_head)                 :: sighead

 print*,"- DEFINE INPUT GRID OBJECT FOR GAUSSIAN DATA."

 num_tiles_input_grid = 1

 if (convert_sfc) then
   the_file=trim(data_dir_input_grid) // "/" // trim(sfc_files_input_grid(1))
 elseif (convert_atm) then
   the_file=trim(data_dir_input_grid) // "/" // trim(atm_files_input_grid(1))
 endif

 if (trim(input_type) == "gfs_sigio") then  ! sigio/sfcio format, used by
                                               ! spectral gfs prior to 7/19/2017.

   if (convert_sfc) then   ! sfcio format
     print*,"- OPEN AND READ ", trim(the_file)
     call sfcio_sropen(21, trim(the_file), rc2)
     if (rc2 /= 0) call error_handler("OPENING FILE", rc2)
     call sfcio_srhead(21, sfchead, rc2)
     if (rc2 /= 0) call error_handler("READING FILE", rc2)
     call sfcio_sclose(21, rc2)
     i_input = sfchead%lonb
     j_input = sfchead%latb
   elseif (convert_atm) then ! sigio format
     print*,"- OPEN AND READ ", trim(the_file)
     call sigio_sropen(21, trim(the_file), rc3)
     if (rc3 /= 0) call error_handler("OPENING FILE", rc3)
     call sigio_srhead(21, sighead, rc3)
     if (rc3 /= 0) call error_handler("READING FILE", rc3)
     call sigio_sclose(21, rc3)
     i_input = sighead%lonb
     j_input = sighead%latb
   endif

 elseif (trim(input_type) == "gaussian_netcdf") then

   print*,'- OPEN AND READ: ',trim(the_file)
   rc=nf90_open(trim(the_file),nf90_nowrite,ncid)
   call netcdf_err(rc, 'opening file')

   print*,"- READ grid_xt"
   rc=nf90_inq_dimid(ncid, 'grid_xt', id_grid)
   call netcdf_err(rc, 'reading grid_xt id')
   rc=nf90_inquire_dimension(ncid,id_grid,len=i_input)
   call netcdf_err(rc, 'reading grid_xt')

   print*,"- READ grid_yt"
   rc=nf90_inq_dimid(ncid, 'grid_yt', id_grid)
   call netcdf_err(rc, 'reading grid_yt id')
   rc=nf90_inquire_dimension(ncid,id_grid,len=j_input)
   call netcdf_err(rc, 'reading grid_yt')

   rc = nf90_close(ncid)

 else ! nemsio format

   call nemsio_init(iret=rc)

   print*,"- OPEN AND READ ", trim(the_file)
   call nemsio_open(gfile, the_file, "read", iret=rc)
   if (rc /= 0) call error_handler("OPENING FILE", rc)

   call nemsio_getfilehead(gfile, iret=rc, dimx=i_input, dimy=j_input)
   if (rc /= 0) call error_handler("READING FILE", rc)

   call nemsio_close(gfile)
 
 endif

 ip1_input = i_input + 1
 jp1_input = j_input + 1

 polekindflag(1:2) = esmf_polekind_monopole

 print*,"- CALL GridCreate1PeriDim FOR INPUT GRID."
 input_grid = esmf_gridcreate1peridim(minindex=(/1,1/), &
                                    maxindex=(/i_input,j_input/), &
                                    polekindflag=polekindflag, &
                                    periodicdim=1, &
                                    poledim=2,  &
                                    coordsys=esmf_coordsys_sph_deg, &
                                    regdecomp=(/1,npets/),  &
                                    indexflag=esmf_index_global, rc=rc)
 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
   call error_handler("IN GridCreate1PeriDim", rc)

 print*,"- CALL FieldCreate FOR INPUT GRID LATITUDE."
 latitude_input_grid = esmf_fieldcreate(input_grid, &
                                   typekind=esmf_typekind_r8, &
                                   staggerloc=esmf_staggerloc_center, &
                                   name="input_grid_latitude", rc=rc)
 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
   call error_handler("IN FieldCreate", rc)

 print*,"- CALL FieldCreate FOR INPUT GRID LONGITUDE."
 longitude_input_grid = esmf_fieldcreate(input_grid, &
                                   typekind=esmf_typekind_r8, &
                                   staggerloc=esmf_staggerloc_center, &
                                   name="input_grid_longitude", rc=rc)
 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
   call error_handler("IN FieldCreate", rc)

 allocate(longitude(i_input,j_input))
 allocate(latitude(i_input,j_input))

 deltalon = 360.0_esmf_kind_r8 / real(i_input,kind=esmf_kind_r8)
 do i = 1, i_input
   longitude(i,:) = real((i-1),kind=esmf_kind_r8) * deltalon
 enddo

 allocate(slat(j_input))
 allocate(wlat(j_input))
 call splat(4, j_input, slat, wlat)

 do i = 1, j_input
   latitude(:,i) = 90.0_esmf_kind_r8 - (acos(slat(i))* 180.0_esmf_kind_r8 / &
                  (4.0_esmf_kind_r8*atan(1.0_esmf_kind_r8)))
 enddo

 deallocate(slat, wlat)

 print*,"- CALL FieldScatter FOR INPUT GRID LONGITUDE."
 call esmf_fieldscatter(longitude_input_grid, longitude, rootpet=0, rc=rc)
 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
    call error_handler("IN FieldScatter", rc)

 print*,"- CALL FieldScatter FOR INPUT GRID LATITUDE."
 call esmf_fieldscatter(latitude_input_grid, latitude, rootpet=0, rc=rc)
 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
    call error_handler("IN FieldScatter", rc)

 print*,"- CALL GridAddCoord FOR INPUT GRID."
 call esmf_gridaddcoord(input_grid, &
                        staggerloc=esmf_staggerloc_center, rc=rc)
 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
    call error_handler("IN GridAddCoord", rc)

 print*,"- CALL GridGetCoord FOR INPUT GRID X-COORD."
 nullify(lon_src_ptr)
 call esmf_gridgetcoord(input_grid, &
                        staggerloc=esmf_staggerloc_center, &
                        coorddim=1, &
                        farrayptr=lon_src_ptr, rc=rc)
 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
    call error_handler("IN GridGetCoord", rc)

 print*,"- CALL GridGetCoord FOR INPUT GRID Y-COORD."
 nullify(lat_src_ptr)
 call esmf_gridgetcoord(input_grid, &
                        staggerloc=esmf_staggerloc_center, &
                        coorddim=2, &
                        computationallbound=clb, &
                        computationalubound=cub, &
                        farrayptr=lat_src_ptr, rc=rc)
 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
    call error_handler("IN GridGetCoord", rc)

 do j = clb(2), cub(2)
   do i = clb(1), cub(1)
     lon_src_ptr(i,j) = longitude(i,j)
     if (lon_src_ptr(i,j) > 360.0_esmf_kind_r8) lon_src_ptr(i,j) = lon_src_ptr(i,j) - 360.0_esmf_kind_r8
     lat_src_ptr(i,j) = latitude(i,j)
   enddo
 enddo

 print*,"- CALL GridAddCoord FOR INPUT GRID."
 call esmf_gridaddcoord(input_grid, &
                        staggerloc=esmf_staggerloc_corner, rc=rc)
 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
    call error_handler("IN GridAddCoord", rc)

 print*,"- CALL GridGetCoord FOR INPUT GRID X-COORD."
 nullify(lon_corner_src_ptr)
 call esmf_gridgetcoord(input_grid, &
                        staggerloc=esmf_staggerloc_corner, &
                        coorddim=1, &
                        farrayptr=lon_corner_src_ptr, rc=rc)
 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
    call error_handler("IN GridGetCoord", rc)

 print*,"- CALL GridGetCoord FOR INPUT GRID Y-COORD."
 nullify(lat_corner_src_ptr)
 call esmf_gridgetcoord(input_grid, &
                        staggerloc=esmf_staggerloc_corner, &
                        coorddim=2, &
                        computationallbound=clb, &
                        computationalubound=cub, &
                        farrayptr=lat_corner_src_ptr, rc=rc)
 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
    call error_handler("IN GridGetCoord", rc)

 do j = clb(2), cub(2)
   do i = clb(1), cub(1)
     lon_corner_src_ptr(i,j) = longitude(i,1) - (0.5_esmf_kind_r8*deltalon)
     if (lon_corner_src_ptr(i,j) > 360.0_esmf_kind_r8) lon_corner_src_ptr(i,j) = lon_corner_src_ptr(i,j) - 360.0_esmf_kind_r8
     if (j == 1) then 
       lat_corner_src_ptr(i,j) = 90.0_esmf_kind_r8
       cycle
     endif
     if (j == jp1_input) then
       lat_corner_src_ptr(i,j) = -90.0_esmf_kind_r8
       cycle
     endif
     lat_corner_src_ptr(i,j) = 0.5_esmf_kind_r8 * (latitude(i,j-1)+ latitude(i,j))
   enddo
 enddo

 deallocate(latitude,longitude)

 end subroutine define_input_grid_gaussian

 subroutine define_input_grid_mosaic(localpet, npets)

 use netcdf
 use program_setup, only       : mosaic_file_input_grid,  &
                                 orog_dir_input_grid, &
                                 orog_files_input_grid

 implicit none

 character(len=500)           :: the_file

 integer, intent(in)          :: localpet, npets

 integer                      :: id_tiles, id_dim, tile
 integer                      :: extra, error, ncid
 integer, allocatable         :: decomptile(:,:)

 integer(esmf_kind_i8), allocatable    :: landmask_one_tile(:,:)

 real(esmf_kind_r8), allocatable       :: latitude_one_tile(:,:)
 real(esmf_kind_r8), allocatable       :: latitude_s_one_tile(:,:)
 real(esmf_kind_r8), allocatable       :: latitude_w_one_tile(:,:)
 real(esmf_kind_r8), allocatable       :: longitude_one_tile(:,:)
 real(esmf_kind_r8), allocatable       :: longitude_s_one_tile(:,:)
 real(esmf_kind_r8), allocatable       :: longitude_w_one_tile(:,:)

 print*,'- OPEN INPUT GRID MOSAIC FILE: ',trim(mosaic_file_input_grid)
 error=nf90_open(trim(mosaic_file_input_grid),nf90_nowrite,ncid)
 call netcdf_err(error, 'opening grid mosaic file')

 print*,"- READ NUMBER OF TILES"
 error=nf90_inq_dimid(ncid, 'ntiles', id_tiles)
 call netcdf_err(error, 'reading ntiles id')
 error=nf90_inquire_dimension(ncid,id_tiles,len=num_tiles_input_grid)
 call netcdf_err(error, 'reading ntiles')

 error = nf90_close(ncid)

 print*,'- NUMBER OF TILES, INPUT MODEL GRID IS ', num_tiles_input_grid

 if (mod(npets,num_tiles_input_grid) /= 0) then
   call error_handler("MUST RUN WITH A TASK COUNT THAT IS A MULTIPLE OF 6.", 1)
 endif

!-----------------------------------------------------------------------
! Create ESMF grid object for the model grid.
!-----------------------------------------------------------------------

 extra = npets / num_tiles_input_grid

 allocate(decomptile(2,num_tiles_input_grid))

 do tile = 1, num_tiles_input_grid
   decomptile(:,tile)=(/1,extra/)
 enddo

 print*,"- CALL GridCreateMosaic FOR INPUT MODEL GRID"
 input_grid = esmf_gridcreatemosaic(filename=trim(mosaic_file_input_grid), &
                                  regdecompptile=decomptile, &
                                  staggerloclist=(/esmf_staggerloc_center, esmf_staggerloc_corner, &
                                                   esmf_staggerloc_edge1, esmf_staggerloc_edge2/), &
                                  indexflag=esmf_index_global, &
                                  tilefilepath=trim(orog_dir_input_grid), &
                                  rc=error)
 if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
    call error_handler("IN GridCreateMosaic", error)

!-----------------------------------------------------------------------
! Read the mask and lat/lons.
!-----------------------------------------------------------------------

 print*,"- CALL FieldCreate FOR INPUT GRID LATITUDE."
 latitude_input_grid = esmf_fieldcreate(input_grid, &
                                   typekind=esmf_typekind_r8, &
                                   staggerloc=esmf_staggerloc_center, &
                                   name="input_grid_latitude", &
                                   rc=error)
 if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
    call error_handler("IN FieldCreate", error)

 print*,"- CALL FieldCreate FOR INPUT GRID LONGITUDE."
 longitude_input_grid = esmf_fieldcreate(input_grid, &
                                   typekind=esmf_typekind_r8, &
                                   staggerloc=esmf_staggerloc_center, &
                                   name="input_grid_longitude", &
                                   rc=error)
 if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
    call error_handler("IN FieldCreate", error)

 print*,"- CALL FieldCreate FOR INPUT GRID LATITUDE_S."
 latitude_s_input_grid = esmf_fieldcreate(input_grid, &
                                   typekind=esmf_typekind_r8, &
                                   staggerloc=esmf_staggerloc_edge2, &
                                   name="input_grid_latitude_s", &
                                   rc=error)
 if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
    call error_handler("IN FieldCreate", error)

 print*,"- CALL FieldCreate FOR INPUT GRID LONGITUDE_S."
 longitude_s_input_grid = esmf_fieldcreate(input_grid, &
                                   typekind=esmf_typekind_r8, &
                                   staggerloc=esmf_staggerloc_edge2, &
                                   name="input_grid_longitude_s", &
                                   rc=error)
 if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
    call error_handler("IN FieldCreate", error)

 print*,"- CALL FieldCreate FOR INPUT GRID LATITUDE_W."
 latitude_w_input_grid = esmf_fieldcreate(input_grid, &
                                   typekind=esmf_typekind_r8, &
                                   staggerloc=esmf_staggerloc_edge1, &
                                   name="input_grid_latitude_w", &
                                   rc=error)
 if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
    call error_handler("IN FieldCreate", error)

 print*,"- CALL FieldCreate FOR INPUT GRID LONGITUDE_W."
 longitude_w_input_grid = esmf_fieldcreate(input_grid, &
                                   typekind=esmf_typekind_r8, &
                                   staggerloc=esmf_staggerloc_edge1, &
                                   name="input_grid_longitude_w", &
                                   rc=error)
 if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
    call error_handler("IN FieldCreate", error)

 the_file = trim(orog_dir_input_grid) // trim(orog_files_input_grid(1))

 print*,'- OPEN FIRST INPUT GRID OROGRAPHY FILE: ',trim(the_file)
 error=nf90_open(trim(the_file),nf90_nowrite,ncid)
 call netcdf_err(error, 'opening ororgraphy file')
 print*,"- READ GRID DIMENSIONS"
 error=nf90_inq_dimid(ncid, 'lon', id_dim)
 call netcdf_err(error, 'reading lon id')
 error=nf90_inquire_dimension(ncid,id_dim,len=i_input)
 call netcdf_err(error, 'reading lon')
 error=nf90_inq_dimid(ncid, 'lat', id_dim)
 call netcdf_err(error, 'reading lat id')
 error=nf90_inquire_dimension(ncid,id_dim,len=j_input)
 call netcdf_err(error, 'reading lat')
 error = nf90_close(ncid)

 print*,"- I/J DIMENSIONS OF THE INPUT GRID TILES ", i_input, j_input

 ip1_input = i_input + 1
 jp1_input = j_input + 1

 if (localpet == 0) then
   allocate(longitude_one_tile(i_input,j_input))
   allocate(longitude_s_one_tile(i_input,jp1_input))
   allocate(longitude_w_one_tile(ip1_input,j_input))
   allocate(latitude_one_tile(i_input,j_input))
   allocate(latitude_s_one_tile(i_input,jp1_input))
   allocate(latitude_w_one_tile(ip1_input,j_input))
   allocate(landmask_one_tile(i_input,j_input))
 else
   allocate(longitude_one_tile(0,0))
   allocate(longitude_s_one_tile(0,0))
   allocate(longitude_w_one_tile(0,0))
   allocate(latitude_one_tile(0,0))
   allocate(latitude_s_one_tile(0,0))
   allocate(latitude_w_one_tile(0,0))
   allocate(landmask_one_tile(0,0))
 endif

 do tile = 1, num_tiles_input_grid
   if (localpet == 0) then
     call get_model_latlons(mosaic_file_input_grid, orog_dir_input_grid, num_tiles_input_grid, tile, &
                            i_input, j_input, ip1_input, jp1_input, latitude_one_tile, &
                            latitude_s_one_tile, latitude_w_one_tile, longitude_one_tile, &
                            longitude_s_one_tile, longitude_w_one_tile)
   endif
   print*,"- CALL FieldScatter FOR INPUT GRID LATITUDE. TILE IS: ", tile
   call esmf_fieldscatter(latitude_input_grid, latitude_one_tile, rootpet=0, tile=tile, rc=error)
   if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
     call error_handler("IN FieldScatter", error)
   print*,"- CALL FieldScatter FOR INPUT GRID LONGITUDE. TILE IS: ", tile
   call esmf_fieldscatter(longitude_input_grid, longitude_one_tile, rootpet=0, tile=tile, rc=error)
   if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
     call error_handler("IN FieldScatter", error)
   print*,"- CALL FieldScatter FOR INPUT GRID LATITUDE_S. TILE IS: ", tile
   call esmf_fieldscatter(latitude_s_input_grid, latitude_s_one_tile, rootpet=0, tile=tile, rc=error)
   if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
     call error_handler("IN FieldScatter", error)
   print*,"- CALL FieldScatter FOR INPUT GRID LONGITUDE_S. TILE IS: ", tile
   call esmf_fieldscatter(longitude_s_input_grid, longitude_s_one_tile, rootpet=0, tile=tile, rc=error)
   if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
     call error_handler("IN FieldScatter", error)
   print*,"- CALL FieldScatter FOR INPUT GRID LATITUDE_W. TILE IS: ", tile
   call esmf_fieldscatter(latitude_w_input_grid, latitude_w_one_tile, rootpet=0, tile=tile, rc=error)
   if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
     call error_handler("IN FieldScatter", error)
   print*,"- CALL FieldScatter FOR INPUT GRID LONGITUDE_W. TILE IS: ", tile
   call esmf_fieldscatter(longitude_w_input_grid, longitude_w_one_tile, rootpet=0, tile=tile, rc=error)
   if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
     call error_handler("IN FieldScatter", error)
 enddo

 deallocate(longitude_one_tile)
 deallocate(longitude_s_one_tile)
 deallocate(longitude_w_one_tile)
 deallocate(latitude_one_tile)
 deallocate(latitude_s_one_tile)
 deallocate(latitude_w_one_tile)
 deallocate(landmask_one_tile)

 end subroutine define_input_grid_mosaic

 subroutine define_input_grid_grib2(localpet,npets)

 use grib_mod
 use gdswzd_mod
 use program_setup, only       : grib2_file_input_grid, data_dir_input_grid

 implicit none

 integer, intent(in)              :: localpet, npets

 character(len=500)               :: the_file

 integer                          :: i, j, k, jdisc, jgdtn, jpdtn, lugb, lugi
 integer                          :: jids(200), jgdt(200), jpdt(200), rc
 integer                          :: kgds(200), nret, clb(2), cub(2)

 logical                          :: unpack

 real                             :: res
 real, allocatable                :: rlon(:,:),rlat(:,:),xpts(:,:),ypts(:,:)
 real, allocatable                :: rlon_corner(:,:),rlat_corner(:,:)
 real, allocatable                :: rlon_diff(:,:),rlat_diff(:,:)
 real, allocatable                :: xpts_corner(:,:),ypts_corner(:,:)
 real(esmf_kind_r8), allocatable  :: latitude(:,:)
 real(esmf_kind_r8), allocatable  :: longitude(:,:)
 real(esmf_kind_r8), allocatable  :: latitude_corner(:,:)
 real(esmf_kind_r8), allocatable  :: longitude_corner(:,:)
 real(esmf_kind_r8), pointer      :: lat_src_ptr(:,:)
 real(esmf_kind_r8), pointer      :: lat_corner_src_ptr(:,:)
 real(esmf_kind_r8), pointer      :: lon_src_ptr(:,:)
 real(esmf_kind_r8), pointer      :: lon_corner_src_ptr(:,:)

 type(esmf_polekind_flag)         :: polekindflag(2)

 type(gribfield)                  :: gfld

 the_file = trim(data_dir_input_grid) // "/" // grib2_file_input_grid

 lugb=12

 print*,"- OPEN AND READ INPUT DATA GRIB2 FILE: ", trim(the_file)
 call baopenr(lugb,the_file,rc)
 if (rc /= 0) call error_handler("OPENING FILE", rc)

! Read the first record and get the grid definition template.

 j       = 0      ! Search at beginning of file
 lugi    = 0      ! No grib index file
 jdisc   = -1     ! Search for any discipline
 jpdtn   = -1     ! Search for any product definition template number
 jgdtn   = -1     ! Search for any grid definition template number
 jids    = -9999  ! Array of values in identification section, set to wildcard.
 jgdt    = -9999  ! Array of values in grid definition template, set to wildcard.
 jpdt    = -9999  ! Array of values in product definition template, set to wildcard.
 unpack  = .false. ! unpack data
   
 call getgb2(lugb, lugi, j, jdisc, jids, jpdtn, jpdt, jgdtn, jgdt, &
             unpack, k, gfld, rc)
 if (rc /= 0) call error_handler("DEGRIBBING INPUT FILE.", rc)

 call baclose(lugb,rc)

 if (gfld%igdtnum == 0) then
   print*,"- INPUT DATA ON LAT/LON GRID."
   input_grid_type = 'latlon'
 elseif (gfld%igdtnum == 30) then
   print*,"- INPUT DATA ON LAMBERT CONFORMAL GRID."
   input_grid_type = 'lambert'
 elseif (gfld%igdtnum == 32769) then
   print*,"- INPUT DATA ON ROTATED LAT/LON GRID."
   input_grid_type = 'rotated_latlon'
 else
   call error_handler("INPUT GRID TEMPLATE NOT SUPPORTED.", 2)
 endif

 kgds = 0
 call gdt_to_gds(gfld%igdtnum, gfld%igdtmpl, gfld%igdtlen, kgds, i_input, j_input, res)

 ip1_input = i_input + 1
 jp1_input = j_input + 1

 allocate(rlat(i_input,j_input))
 allocate(rlon(i_input,j_input))
 allocate(rlat_diff(i_input,j_input))
 allocate(rlon_diff(i_input,j_input))
 allocate(xpts(i_input,j_input))
 allocate(ypts(i_input,j_input))
 allocate(rlat_corner(ip1_input,jp1_input))
 allocate(rlon_corner(ip1_input,jp1_input))
 allocate(xpts_corner(ip1_input,jp1_input))
 allocate(ypts_corner(ip1_input,jp1_input))

 do j = 1, j_input
 do i = 1, i_input
   xpts(i,j) = float(i)
   ypts(i,j) = float(j)
 enddo
 enddo

 print*,"- COMPUTE GRID CELL CENTER COORDINATES."
 call gdswzd(kgds,1,(i_input*j_input),-9999.,xpts,ypts,rlon,rlat,nret)

 if (nret /= (i_input*j_input)) then
   call error_handler("GDSWZD RETURNED WRONG NUMBER OF POINTS.", 2)
 endif

 deallocate(xpts, ypts)

 do j = 1, jp1_input
 do i = 1, ip1_input
   xpts_corner(i,j) = float(i) - 0.5
   ypts_corner(i,j) = float(j) - 0.5
 enddo
 enddo

 print*,"- COMPUTE GRID CELL CORNER COORDINATES."
 call gdswzd(kgds,1,(ip1_input*jp1_input),-9999.,xpts_corner,ypts_corner,rlon_corner,rlat_corner,nret)

 if (nret /= (ip1_input*jp1_input)) then
   call error_handler("GDSWZD RETURNED WRONG NUMBER OF POINTS.", 2)
 endif

 deallocate(xpts_corner, ypts_corner)

 if (gfld%igdtnum == 0) then ! gfs lat/lon data

   print*,"- CALL GridCreate1PeriDim FOR INPUT GRID."

   polekindflag(1:2) = esmf_polekind_monopole

   input_grid = esmf_gridcreate1peridim(minindex=(/1,1/), &
                                    maxindex=(/i_input,j_input/), &
                                    polekindflag=polekindflag, &
                                    periodicdim=1, &
                                    poledim=2,  &
                                    coordsys=esmf_coordsys_sph_deg, &
                                    regdecomp=(/1,npets/),  &
                                    indexflag=esmf_index_global, rc=rc)

 else

   print*,"- CALL GridCreateNoPeriDim FOR INPUT GRID."

   input_grid = esmf_gridcreatenoperidim(maxindex=(/i_input,j_input/), & 
                                  indexflag=esmf_index_global, &
                                  rc=rc)

 endif

 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
   call error_handler("IN GridCreate1PeriDim", rc)

 print*,"- CALL FieldCreate FOR INPUT GRID LATITUDE."
 latitude_input_grid = esmf_fieldcreate(input_grid, &
                                   typekind=esmf_typekind_r8, &
                                   staggerloc=esmf_staggerloc_center, &
                                   name="input_grid_latitude", rc=rc)
 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
   call error_handler("IN FieldCreate", rc)

 print*,"- CALL FieldCreate FOR INPUT GRID LONGITUDE."
 longitude_input_grid = esmf_fieldcreate(input_grid, &
                                   typekind=esmf_typekind_r8, &
                                   staggerloc=esmf_staggerloc_center, &
                                   name="input_grid_longitude", rc=rc)
 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
   call error_handler("IN FieldCreate", rc)

 allocate(latitude(i_input,j_input))
 allocate(longitude(i_input,j_input))

 latitude = rlat
 longitude = rlon

 deallocate (rlat, rlon)

 print*,"- CALL FieldScatter FOR INPUT GRID LONGITUDE."
 call esmf_fieldscatter(longitude_input_grid, longitude, rootpet=0, rc=rc)
 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
    call error_handler("IN FieldScatter", rc)

 print*,"- CALL FieldScatter FOR INPUT GRID LATITUDE."
 call esmf_fieldscatter(latitude_input_grid, latitude, rootpet=0, rc=rc)
 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
    call error_handler("IN FieldScatter", rc)

 print*,"- CALL GridAddCoord FOR INPUT GRID."
 call esmf_gridaddcoord(input_grid, &
                        staggerloc=esmf_staggerloc_center, rc=rc)
 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
    call error_handler("IN GridAddCoord", rc)

 print*,"- CALL GridGetCoord FOR INPUT GRID X-COORD."
 nullify(lon_src_ptr)
 call esmf_gridgetcoord(input_grid, &
                        staggerloc=esmf_staggerloc_center, &
                        coorddim=1, &
                        farrayptr=lon_src_ptr, rc=rc)
 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
    call error_handler("IN GridGetCoord", rc)

 print*,"- CALL GridGetCoord FOR INPUT GRID Y-COORD."
 nullify(lat_src_ptr)
 call esmf_gridgetcoord(input_grid, &
                        staggerloc=esmf_staggerloc_center, &
                        coorddim=2, &
                        computationallbound=clb, &
                        computationalubound=cub, &
                        farrayptr=lat_src_ptr, rc=rc)
 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
    call error_handler("IN GridGetCoord", rc)

 do j = clb(2), cub(2)
   do i = clb(1), cub(1)
     lon_src_ptr(i,j) = longitude(i,j)
     if (lon_src_ptr(i,j) > 360.0_esmf_kind_r8) lon_src_ptr(i,j) = lon_src_ptr(i,j) - 360.0_esmf_kind_r8
     lat_src_ptr(i,j) = latitude(i,j)
   enddo
 enddo

 deallocate(latitude, longitude)

 allocate(latitude_corner(ip1_input,jp1_input))
 allocate(longitude_corner(ip1_input,jp1_input))

 latitude_corner = rlat_corner
 longitude_corner = rlon_corner

 deallocate (rlat_corner, rlon_corner)

 print*,"- CALL GridAddCoord FOR INPUT GRID."
 call esmf_gridaddcoord(input_grid, &
                        staggerloc=esmf_staggerloc_corner, rc=rc)
 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
    call error_handler("IN GridAddCoord", rc)

 print*,"- CALL GridGetCoord FOR INPUT GRID X-COORD."
 nullify(lon_corner_src_ptr)
 call esmf_gridgetcoord(input_grid, &
                        staggerloc=esmf_staggerloc_corner, &
                        coorddim=1, &
                        farrayptr=lon_corner_src_ptr, rc=rc)
 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
    call error_handler("IN GridGetCoord", rc)

 print*,"- CALL GridGetCoord FOR INPUT GRID Y-COORD."
 nullify(lat_corner_src_ptr)
 call esmf_gridgetcoord(input_grid, &
                        staggerloc=esmf_staggerloc_corner, &
                        coorddim=2, &
                        computationallbound=clb, &
                        computationalubound=cub, &
                        farrayptr=lat_corner_src_ptr, rc=rc)
 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
    call error_handler("IN GridGetCoord", rc)

 do j = clb(2), cub(2)
   do i = clb(1), cub(1)
     lon_corner_src_ptr(i,j) = longitude_corner(i,j)
     if (lon_corner_src_ptr(i,j) > 360.0_esmf_kind_r8) lon_corner_src_ptr(i,j) = lon_corner_src_ptr(i,j) - 360.0_esmf_kind_r8
     lat_corner_src_ptr(i,j) = latitude_corner(i,j)
   enddo
 enddo

 deallocate(latitude_corner, longitude_corner)

 end subroutine define_input_grid_grib2

 subroutine define_target_grid(localpet, npets)

 use netcdf
 use program_setup, only       : mosaic_file_target_grid, &
                                 orog_dir_target_grid,    &
                                 orog_files_target_grid,  &
                                 nsoill_out

 implicit none

 integer, intent(in)                   :: localpet, npets

 character(len=500)                    :: the_file

 integer                               :: error, ncid, extra
 integer                               :: id_tiles
 integer                               :: id_dim, id_grid_tiles
 integer                               :: tile
 integer, allocatable                  :: decomptile(:,:)
 integer(esmf_kind_i8), allocatable    :: landmask_one_tile(:,:)
 integer(esmf_kind_i8), allocatable    :: seamask_one_tile(:,:)
 
 real(esmf_kind_r8), allocatable       :: latitude_one_tile(:,:)
 real(esmf_kind_r8), allocatable       :: latitude_s_one_tile(:,:)
 real(esmf_kind_r8), allocatable       :: latitude_w_one_tile(:,:)
 real(esmf_kind_r8), allocatable       :: longitude_one_tile(:,:)
 real(esmf_kind_r8), allocatable       :: longitude_s_one_tile(:,:)
 real(esmf_kind_r8), allocatable       :: longitude_w_one_tile(:,:)
 real(esmf_kind_r8), allocatable       :: terrain_one_tile(:,:)

 lsoil_target = nsoill_out
 
 print*,'- OPEN TARGET GRID MOSAIC FILE: ',trim(mosaic_file_target_grid)
 error=nf90_open(trim(mosaic_file_target_grid),nf90_nowrite,ncid)
 call netcdf_err(error, 'opening grid mosaic file')

 print*,"- READ NUMBER OF TILES"
 error=nf90_inq_dimid(ncid, 'ntiles', id_tiles)
 call netcdf_err(error, 'reading ntile id')
 error=nf90_inquire_dimension(ncid,id_tiles,len=num_tiles_target_grid)
 call netcdf_err(error, 'reading ntiles')
 error=nf90_inq_varid(ncid, 'gridtiles', id_grid_tiles)
 call netcdf_err(error, 'reading gridtiles id')
 allocate(tiles_target_grid(num_tiles_target_grid))
 tiles_target_grid="NULL"
 print*,"- READ TILE NAMES"
 error=nf90_get_var(ncid, id_grid_tiles, tiles_target_grid)
 call netcdf_err(error, 'reading gridtiles')

 error = nf90_close(ncid)

 print*,'- NUMBER OF TILES, TARGET MODEL GRID IS ', num_tiles_target_grid

 if (mod(npets,num_tiles_target_grid) /= 0) then
   call error_handler("MUST RUN WITH TASK COUNT THAT IS A MULTIPLE OF # OF TILES.", 1)
 endif

!-----------------------------------------------------------------------
! Get the model grid specs and land mask from the orography files.
!-----------------------------------------------------------------------

 the_file = trim(orog_dir_target_grid) // trim(orog_files_target_grid(1))

 print*,'- OPEN FIRST TARGET GRID OROGRAPHY FILE: ',trim(the_file)
 error=nf90_open(trim(the_file),nf90_nowrite,ncid)
 call netcdf_err(error, 'opening orography file')
 print*,"- READ GRID DIMENSIONS"
 error=nf90_inq_dimid(ncid, 'lon', id_dim)
 call netcdf_err(error, 'reading lon id')
 error=nf90_inquire_dimension(ncid,id_dim,len=i_target)
 call netcdf_err(error, 'reading lon')
 error=nf90_inq_dimid(ncid, 'lat', id_dim)
 call netcdf_err(error, 'reading lat id')
 error=nf90_inquire_dimension(ncid,id_dim,len=j_target)
 call netcdf_err(error, 'reading lat')
 error = nf90_close(ncid)

 print*,"- I/J DIMENSIONS OF THE TARGET GRID TILES ", i_target, j_target

 ip1_target = i_target + 1
 jp1_target = j_target + 1

!-----------------------------------------------------------------------
! Create ESMF grid object for the model grid.
!-----------------------------------------------------------------------

 extra = npets / num_tiles_target_grid

 allocate(decomptile(2,num_tiles_target_grid))

 do tile = 1, num_tiles_target_grid
   decomptile(:,tile)=(/1,extra/)
 enddo

 print*,"- CALL GridCreateMosaic FOR TARGET GRID"
 target_grid = esmf_gridcreatemosaic(filename=trim(mosaic_file_target_grid), &
                                  regdecompptile=decomptile, &
                                  staggerloclist=(/esmf_staggerloc_center, esmf_staggerloc_corner, &
                                                   esmf_staggerloc_edge1, esmf_staggerloc_edge2/), &
                                  indexflag=esmf_index_global, &
                                  tilefilepath=trim(orog_dir_target_grid), rc=error)
 if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
    call error_handler("IN GridCreateMosaic", error)

!-----------------------------------------------------------------------
! Set target model landmask (1 - land, 0 - not land) and 
! seamask (1 - non-land, 0 -land).  Read lat/lon on target grid.
!-----------------------------------------------------------------------

 print*,"- CALL FieldCreate FOR TARGET GRID LANDMASK."
 landmask_target_grid = esmf_fieldcreate(target_grid, &
                                   typekind=esmf_typekind_i8, &
                                   staggerloc=esmf_staggerloc_center, &
                                   name="target_grid_landmask", rc=error)
 if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
    call error_handler("IN FieldCreate", error)

 print*,"- CALL FieldCreate FOR TARGET GRID SEAMASK."
 seamask_target_grid = esmf_fieldcreate(target_grid, &
                                   typekind=esmf_typekind_i8, &
                                   staggerloc=esmf_staggerloc_center, &
                                   name="target_grid_seamask", rc=error)
 if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
    call error_handler("IN FieldCreate", error)

 print*,"- CALL FieldCreate FOR TARGET GRID LATITUDE."
 latitude_target_grid = esmf_fieldcreate(target_grid, &
                                   typekind=esmf_typekind_r8, &
                                   staggerloc=esmf_staggerloc_center, &
                                   name="target_grid_latitude", rc=error)
 if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
    call error_handler("IN FieldCreate", error)

 print*,"- CALL FieldCreate FOR TARGET GRID LATITUDE_S."
 latitude_s_target_grid = esmf_fieldcreate(target_grid, &
                                   typekind=esmf_typekind_r8, &
                                   staggerloc=esmf_staggerloc_edge2, &
                                   name="target_grid_latitude_s", rc=error)
 if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
    call error_handler("IN FieldCreate", error)

 print*,"- CALL FieldCreate FOR TARGET GRID LATITUDE_W."
 latitude_w_target_grid = esmf_fieldcreate(target_grid, &
                                   typekind=esmf_typekind_r8, &
                                   staggerloc=esmf_staggerloc_edge1, &
                                   name="target_grid_latitude_w", &
                                   rc=error)
 if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
    call error_handler("IN FieldCreate", error)

 print*,"- CALL FieldCreate FOR TARGET GRID LONGITUDE."
 longitude_target_grid = esmf_fieldcreate(target_grid, &
                                   typekind=esmf_typekind_r8, &
                                   staggerloc=esmf_staggerloc_center, &
                                   name="target_grid_longitude", &
                                   rc=error)
 if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
    call error_handler("IN FieldCreate", error)

 print*,"- CALL FieldCreate FOR TARGET GRID LONGITUDE_S."
 longitude_s_target_grid = esmf_fieldcreate(target_grid, &
                                   typekind=esmf_typekind_r8, &
                                   staggerloc=esmf_staggerloc_edge2, &
                                   name="target_grid_longitude_s", &
                                   rc=error)
 if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
    call error_handler("IN FieldCreate", error)

 print*,"- CALL FieldCreate FOR TARGET GRID LONGITUDE_W."
 longitude_w_target_grid = esmf_fieldcreate(target_grid, &
                                   typekind=esmf_typekind_r8, &
                                   staggerloc=esmf_staggerloc_edge1, &
                                   name="target_grid_longitude_w", &
                                   rc=error)
 if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
    call error_handler("IN FieldCreate", error)

 print*,"- CALL FieldCreate FOR TARGET GRID TERRAIN."
 terrain_target_grid = esmf_fieldcreate(target_grid, &
                                   typekind=esmf_typekind_r8, &
                                   staggerloc=esmf_staggerloc_center, &
                                   name="target_grid_terrain", &
                                   rc=error)
 if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
    call error_handler("IN FieldCreate", error)

 if (localpet == 0) then
   allocate(landmask_one_tile(i_target,j_target))
   allocate(seamask_one_tile(i_target,j_target))
   allocate(latitude_one_tile(i_target,j_target))
   allocate(latitude_s_one_tile(i_target,jp1_target))
   allocate(latitude_w_one_tile(ip1_target,j_target))
   allocate(longitude_one_tile(i_target,j_target))
   allocate(longitude_s_one_tile(i_target,jp1_target))
   allocate(longitude_w_one_tile(ip1_target,j_target))
   allocate(terrain_one_tile(i_target,j_target))
 else
   allocate(landmask_one_tile(0,0))
   allocate(seamask_one_tile(0,0))
   allocate(longitude_one_tile(0,0))
   allocate(longitude_s_one_tile(0,0))
   allocate(longitude_w_one_tile(0,0))
   allocate(latitude_one_tile(0,0))
   allocate(latitude_s_one_tile(0,0))
   allocate(latitude_w_one_tile(0,0))
   allocate(terrain_one_tile(0,0))
 endif

 do tile = 1, num_tiles_target_grid
   if (localpet == 0) then
     the_file = trim(orog_dir_target_grid) // trim(orog_files_target_grid(tile))
     call get_model_mask_terrain(trim(the_file), i_target, j_target, landmask_one_tile, &
                                 terrain_one_tile)
     seamask_one_tile = 0
     where(landmask_one_tile == 0) seamask_one_tile = 1
     call get_model_latlons(mosaic_file_target_grid, orog_dir_target_grid, num_tiles_target_grid, tile, &
                            i_target, j_target, ip1_target, jp1_target, latitude_one_tile, &
                            latitude_s_one_tile, latitude_w_one_tile, longitude_one_tile, &
                            longitude_s_one_tile, longitude_w_one_tile)
   endif
   print*,"- CALL FieldScatter FOR TARGET GRID LANDMASK. TILE IS: ", tile
   call esmf_fieldscatter(landmask_target_grid, landmask_one_tile, rootpet=0, tile=tile, rc=error)
   if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
      call error_handler("IN FieldScatter", error)
   print*,"- CALL FieldScatter FOR TARGET GRID SEAMASK. TILE IS: ", tile
   call esmf_fieldscatter(seamask_target_grid, seamask_one_tile, rootpet=0, tile=tile, rc=error)
   if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
      call error_handler("IN FieldScatter", error)
   print*,"- CALL FieldScatter FOR TARGET GRID LONGITUDE. TILE IS: ", tile
   call esmf_fieldscatter(longitude_target_grid, longitude_one_tile, rootpet=0, tile=tile, rc=error)
   if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
      call error_handler("IN FieldScatter", error)
   print*,"- CALL FieldScatter FOR TARGET GRID LONGITUDE_S. TILE IS: ", tile
   call esmf_fieldscatter(longitude_s_target_grid, longitude_s_one_tile, rootpet=0, tile=tile, rc=error)
   if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
      call error_handler("IN FieldScatter", error)
   print*,"- CALL FieldScatter FOR TARGET GRID LONGITUDE_W. TILE IS: ", tile
   call esmf_fieldscatter(longitude_w_target_grid, longitude_w_one_tile, rootpet=0, tile=tile, rc=error)
   if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
      call error_handler("IN FieldScatter", error)
   print*,"- CALL FieldScatter FOR TARGET GRID LATITUDE. TILE IS: ", tile
   call esmf_fieldscatter(latitude_target_grid, latitude_one_tile, rootpet=0, tile=tile, rc=error)
   if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
      call error_handler("IN FieldScatter", error)
   print*,"- CALL FieldScatter FOR TARGET GRID LATITUDE_S. TILE IS: ", tile
   call esmf_fieldscatter(latitude_s_target_grid, latitude_s_one_tile, rootpet=0, tile=tile, rc=error)
   if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
      call error_handler("IN FieldScatter", error)
   print*,"- CALL FieldScatter FOR TARGET GRID LATITUDE_W. TILE IS: ", tile
   call esmf_fieldscatter(latitude_w_target_grid, latitude_w_one_tile, rootpet=0, tile=tile, rc=error)
   if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
      call error_handler("IN FieldScatter", error)
   print*,"- CALL FieldScatter FOR TARGET GRID TERRAIN. TILE IS: ", tile
   call esmf_fieldscatter(terrain_target_grid, terrain_one_tile, rootpet=0, tile=tile, rc=error)
   if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
      call error_handler("IN FieldScatter", error)
 enddo

 deallocate(landmask_one_tile)
 deallocate(seamask_one_tile)
 deallocate(longitude_one_tile)
 deallocate(longitude_s_one_tile)
 deallocate(longitude_w_one_tile)
 deallocate(latitude_one_tile)
 deallocate(latitude_s_one_tile)
 deallocate(latitude_w_one_tile)
 deallocate(terrain_one_tile)

 end subroutine define_target_grid

 subroutine get_model_latlons(mosaic_file, orog_dir, num_tiles, tile, &
                              i_tile, j_tile, ip1_tile, jp1_tile,  &
                              latitude, latitude_s, latitude_w, &
                              longitude, longitude_s, longitude_w)

 use netcdf

 implicit none

 character(len=*), intent(in)      :: mosaic_file, orog_dir

 integer, intent(in)               :: num_tiles, tile
 integer, intent(in)               :: i_tile, j_tile
 integer, intent(in)               :: ip1_tile, jp1_tile

 real(esmf_kind_r8), intent(out)   :: latitude(i_tile, j_tile)
 real(esmf_kind_r8), intent(out)   :: latitude_s(i_tile, jp1_tile)
 real(esmf_kind_r8), intent(out)   :: latitude_w(ip1_tile, j_tile)
 real(esmf_kind_r8), intent(out)   :: longitude(i_tile, j_tile)
 real(esmf_kind_r8), intent(out)   :: longitude_s(i_tile, jp1_tile)
 real(esmf_kind_r8), intent(out)   :: longitude_w(ip1_tile, j_tile)

 character(len=25)                 :: grid_files(num_tiles)
 character(len=255)                :: grid_file

 integer                           :: error, id_var, ncid
 integer                           :: id_dim, nxp, nyp, i, j, ii, jj

 real(esmf_kind_r8), allocatable   :: tmpvar(:,:)

 print*,"- READ MODEL GRID FILE"

 print*,'- OPEN MOSAIC FILE: ', trim(mosaic_file)
 error=nf90_open(trim(mosaic_file), nf90_nowrite, ncid)
 call netcdf_err(error, 'opening mosaic file')

 print*,"- READ GRID FILE NAMES"
 error=nf90_inq_varid(ncid, 'gridfiles', id_var)
 call netcdf_err(error, 'reading gridfiles id')
 error=nf90_get_var(ncid, id_var, grid_files)
 call netcdf_err(error, 'reading gridfiles')

 error = nf90_close(ncid)

 grid_file = trim(orog_dir) // trim(grid_files(tile))

 print*,'- OPEN GRID FILE: ', trim(grid_file)
 error=nf90_open(trim(grid_file), nf90_nowrite, ncid)
 call netcdf_err(error, 'opening grid file')

 print*,'- READ NXP ID'
 error=nf90_inq_dimid(ncid, 'nxp', id_dim)
 call netcdf_err(error, 'reading nxp id')

 print*,'- READ NXP'
 error=nf90_inquire_dimension(ncid,id_dim,len=nxp)
 call netcdf_err(error, 'reading nxp')

 print*,'- READ NYP ID'
 error=nf90_inq_dimid(ncid, 'nyp', id_dim)
 call netcdf_err(error, 'reading nyp id')

 print*,'- READ NYP'
 error=nf90_inquire_dimension(ncid,id_dim,len=nyp)
 call netcdf_err(error, 'reading nyp')

 if ((nxp/2 /= i_tile) .or. (nyp/2 /= j_tile)) then
   call error_handler("DIMENSION MISMATCH IN GRID FILE.", 1)
 endif

 allocate(tmpvar(nxp,nyp))

 print*,'- READ LONGITUDE ID'
 error=nf90_inq_varid(ncid, 'x', id_var)
 call netcdf_err(error, 'reading longitude id')

 print*,'- READ LONGITUDE'
 error=nf90_get_var(ncid, id_var, tmpvar)
 call netcdf_err(error, 'reading longitude')

 do j = 1, j_tile
 do i = 1, i_tile
   ii = 2*i
   jj = 2*j
   longitude(i,j) = tmpvar(ii,jj)
 enddo
 enddo

 do j = 1, jp1_tile
 do i = 1, i_tile
   ii = 2*i
   jj = (2*j) - 1
   longitude_s(i,j) = tmpvar(ii,jj)
 enddo
 enddo

 do j = 1, j_tile
 do i = 1, ip1_tile
   ii = (2*i) - 1
   jj = 2*j
   longitude_w(i,j) = tmpvar(ii,jj)
 enddo
 enddo

 print*,'- READ LATITUDE ID'
 error=nf90_inq_varid(ncid, 'y', id_var)
 call netcdf_err(error, 'reading latitude id')

 print*,'- READ LATIITUDE'
 error=nf90_get_var(ncid, id_var, tmpvar)
 call netcdf_err(error, 'reading latitude')

 do j = 1, j_tile
 do i = 1, i_tile
   ii = 2*i
   jj = 2*j
   latitude(i,j) = tmpvar(ii,jj)
 enddo
 enddo

 do j = 1, jp1_tile
 do i = 1, i_tile
   ii = 2*i
   jj = (2*j) - 1
   latitude_s(i,j) = tmpvar(ii,jj)
 enddo
 enddo

 do j = 1, j_tile
 do i = 1, ip1_tile
   ii = (2*i) - 1
   jj = 2*j
   latitude_w(i,j) = tmpvar(ii,jj)
 enddo
 enddo

 deallocate(tmpvar)

 error = nf90_close(ncid)

 end subroutine get_model_latlons
 
 subroutine get_model_mask_terrain(orog_file, idim, jdim, mask, terrain)

 use netcdf

 implicit none

 character(len=*), intent(in)       :: orog_file

 integer, intent(in)                :: idim, jdim
 integer(esmf_kind_i8), intent(out) :: mask(idim,jdim)

 real(esmf_kind_i8), intent(out)    :: terrain(idim,jdim)

 integer                            :: error, lat, lon
 integer                            :: ncid, id_dim, id_var

 real(kind=4), allocatable          :: dummy(:,:)

 print*,"- READ MODEL LAND MASK FILE"

 print*,'- OPEN LAND MASK FILE: ', orog_file
 error=nf90_open(orog_file,nf90_nowrite,ncid)
 call netcdf_err(error, 'opening land mask file')

 print*,"- READ I-DIMENSION"
 error=nf90_inq_dimid(ncid, 'lon', id_dim)
 call netcdf_err(error, 'reading idim id')
 error=nf90_inquire_dimension(ncid,id_dim,len=lon)
 call netcdf_err(error, 'reading idim')

 print*,"- READ J-DIMENSION"
 error=nf90_inq_dimid(ncid, 'lat', id_dim)
 call netcdf_err(error, 'reading jdim id')
 error=nf90_inquire_dimension(ncid,id_dim,len=lat)
 call netcdf_err(error, 'reading jdim')

 print*,"- I/J DIMENSIONS: ", lon, lat

 if ((lon /= idim) .or. (lat /= jdim)) then
   call error_handler("MISMATCH IN DIMENSIONS.", 1)
 endif

 allocate(dummy(idim,jdim))

 print*,"- READ LAND MASK"
 error=nf90_inq_varid(ncid, 'slmsk', id_var)
 call netcdf_err(error, 'reading slmsk id')
 error=nf90_get_var(ncid, id_var, dummy)
 call netcdf_err(error, 'reading slmsk')
 mask = nint(dummy)

 print*,"- READ RAW OROGRAPHY."
 error=nf90_inq_varid(ncid, 'orog_raw', id_var)
 call netcdf_err(error, 'reading orog_raw id')
 error=nf90_get_var(ncid, id_var, dummy)
 call netcdf_err(error, 'reading orog_raw')
 terrain = dummy

 error = nf90_close(ncid)

 deallocate (dummy)

 end subroutine get_model_mask_terrain

 subroutine cleanup_input_target_grid_data

 implicit none

 integer                                :: rc

 print*,"- DESTROY MODEL DATA."
 
 call esmf_fielddestroy(latitude_input_grid,rc=rc)
 call esmf_fielddestroy(longitude_input_grid,rc=rc)
 if (esmf_fieldiscreated(latitude_s_input_grid)) then
   call esmf_fielddestroy(latitude_s_input_grid, rc=rc)
 endif
 if (esmf_fieldiscreated(latitude_w_input_grid)) then
   call esmf_fielddestroy(latitude_w_input_grid, rc=rc)
 endif
 if (esmf_fieldiscreated(longitude_s_input_grid)) then
   call esmf_fielddestroy(longitude_s_input_grid, rc=rc)
 endif
 if (esmf_fieldiscreated(longitude_w_input_grid)) then
   call esmf_fielddestroy(longitude_w_input_grid, rc=rc)
 endif
 call esmf_fielddestroy(landmask_target_grid, rc=rc)
 call esmf_fielddestroy(latitude_target_grid, rc=rc)
 if (esmf_fieldiscreated(latitude_s_target_grid)) then
   call esmf_fielddestroy(latitude_s_target_grid, rc=rc)
 endif
 if (esmf_fieldiscreated(latitude_w_target_grid)) then
   call esmf_fielddestroy(latitude_w_target_grid, rc=rc)
 endif
 call esmf_fielddestroy(longitude_target_grid, rc=rc)
 if (esmf_fieldiscreated(longitude_s_target_grid)) then
   call esmf_fielddestroy(longitude_s_target_grid, rc=rc)
 endif
 if (esmf_fieldiscreated(longitude_w_target_grid)) then
   call esmf_fielddestroy(longitude_w_target_grid, rc=rc)
 endif
 call esmf_fielddestroy(seamask_target_grid, rc=rc)
 call esmf_fielddestroy(terrain_target_grid, rc=rc)
 call esmf_griddestroy(input_grid, rc=rc)
 call esmf_griddestroy(target_grid, rc=rc)

 end subroutine cleanup_input_target_grid_data

 subroutine gdt_to_gds(igdtnum, igdstmpl, igdtlen, kgds, ni, nj, res)

 implicit none

 integer, intent(in   )  :: igdtnum, igdtlen, igdstmpl(igdtlen)
 integer, intent(  out)  :: kgds(200), ni, nj
 integer                 :: iscale

 real,    intent(  out)  :: res

 kgds=0

 if (igdtnum.eq.32769) then        ! rot lat/lon b grid

     iscale=igdstmpl(10)*igdstmpl(11)
     if (iscale == 0) iscale = 1e6
     kgds(1)=205                    ! oct 6,     rotated lat/lon for Non-E
                                    !            Stagger grid
     kgds(2)=igdstmpl(8)            ! octs 7-8,  Ni
     ni = kgds(2)
     kgds(3)=igdstmpl(9)            ! octs 9-10, Nj
     nj = kgds(3)
     kgds(4)=nint(float(igdstmpl(12))/float(iscale)*1000.)  ! octs 11-13, Lat of
                                                            ! 1st grid point
     kgds(5)=nint(float(igdstmpl(13))/float(iscale)*1000.)  ! octs 14-16, Lon of
                                                            ! 1st grid point

     kgds(6)=0                      ! oct 17, resolution and component flags
     if (igdstmpl(1)==2 ) kgds(6)=64
     if ( btest(igdstmpl(14),4).OR.btest(igdstmpl(14),5) )  kgds(6)=kgds(6)+128
     if ( btest(igdstmpl(14),3) ) kgds(6)=kgds(6)+8

     kgds(7)=nint(float(igdstmpl(15))/float(iscale)*1000.)       ! octs 18-20,
                                                                 ! Lat of cent of rotation
     kgds(8)=nint(float(igdstmpl(16))/float(iscale)*1000.)       ! octs 21-23,
                                                                 ! Lon of cent of rotation
     kgds(9)=nint(float(igdstmpl(17))/float(iscale)*1000.)       ! octs 24-25,
                                                                 ! Di
     kgds(10)=nint(float(igdstmpl(18))/float(iscale)*1000.)      ! octs 26-27,
                                                                 ! Dj

     kgds(11) = 0                   ! oct 28, scan mode
     if (btest(igdstmpl(19),7)) kgds(11) = 128
     if (btest(igdstmpl(19),6)) kgds(11) = kgds(11) +  64
     if (btest(igdstmpl(19),5)) kgds(11) = kgds(11) +  32

     kgds(12)=nint(float(igdstmpl(20))/float(iscale)*1000.) ! octs 29-31, Lat of
                                                            ! last grid point
     kgds(13)=nint(float(igdstmpl(21))/float(iscale)*1000.) ! octs 32-34, Lon of
                                                            ! last grid point

     kgds(19)=0    ! oct 4, # vert coordinate parameters
     kgds(20)=255  ! oct 5, used for thinned grids, set to 255

     res = ((float(kgds(9)) / 1000.0) + (float(kgds(10)) / 1000.0)) &
             * 0.5 * 111.0

   elseif(igdtnum==30) then

     kgds(1)=3                      ! oct 6,     lambert conformal
     kgds(2)=igdstmpl(8)            ! octs 7-8,  Ni
     ni = kgds(2)
     kgds(3)=igdstmpl(9)            ! octs 9-10, Nj
     nj = kgds(3)

     iscale = 1e6
     kgds(4) = nint(float(igdstmpl(10))/1000.0)
     kgds(5) = nint(float(igdstmpl(11))/1000.0)

     kgds(6)=0                      ! oct 17, resolution and component flags
     if (igdstmpl(1)==2 ) kgds(6)=64
     if ( btest(igdstmpl(12),4).OR.btest(igdstmpl(12),5) )  kgds(6)=kgds(6)+128
     if ( btest(igdstmpl(12),3) ) kgds(6)=kgds(6)+8

     kgds(7) = nint(float(igdstmpl(14))/1000.0)
     kgds(8) = nint(float(igdstmpl(15))/1000.0)
     kgds(9) = nint(float(igdstmpl(16))/1000.0)
     kgds(10) = 0 

     kgds(11) = 0                   ! oct 28, scan mode
     if (btest(igdstmpl(18),7)) kgds(11) = 128
     if (btest(igdstmpl(18),6)) kgds(11) = kgds(11) +  64
     if (btest(igdstmpl(18),5)) kgds(11) = kgds(11) +  32

     kgds(12) = nint(float(igdstmpl(19))/1000.0)
     kgds(13) = nint(float(igdstmpl(20))/1000.0)
     kgds(14) = -90
     kgds(15) = 0

   elseif(igdtnum==0) then  ! lat/lon grid

     iscale=igdstmpl(10)*igdstmpl(11)
     if (iscale == 0) iscale = 1e6
     kgds(1)=0                   ! oct 6, data representation type.
     kgds(2)=igdstmpl(8)         ! octs 7-8, Ni
     ni = kgds(2)
     kgds(3)=igdstmpl(9)         ! octs 9-10, Nj
     nj = kgds(3)
     kgds(4)=nint(float(igdstmpl(12))/float(iscale)*1000.)  ! octs 11-13, Lat of 1st grid point
     kgds(5)=nint(float(igdstmpl(13))/float(iscale)*1000.)  ! octs 14-16, Lon of 1st grid point

     kgds(6)=0                   ! oct 17, resolution and component flags
     if (igdstmpl(1)==2 ) kgds(6)=64
     if ( btest(igdstmpl(14),4).OR.btest(igdstmpl(14),5) ) kgds(6)=kgds(6)+128
     if ( btest(igdstmpl(14),3) ) kgds(6)=kgds(6)+8

     kgds(7)=nint(float(igdstmpl(15))/float(iscale)*1000.)  ! octs 18-20, Lat of last grid point
     kgds(8)=nint(float(igdstmpl(16))/float(iscale)*1000.)  ! octs 21-23, Lon of last grid point
     kgds(9)=nint(float(igdstmpl(17))/float(iscale)*1000.)  ! octs 24-25, "i" resolution.
     kgds(10)=nint(float(igdstmpl(18))/float(iscale)*1000.) ! octs 26-27, "j" resolution.

     kgds(11) = 0              ! oct 28, scan mode
     if (btest(igdstmpl(19),7)) kgds(11) = 128
     if (btest(igdstmpl(19),6)) kgds(11) = kgds(11) +  64
     if (btest(igdstmpl(19),5)) kgds(11) = kgds(11) +  32

     kgds(12)=0      ! octs 29-32, reserved
     kgds(19)=0      ! oct 4, # vert coordinate parameters
     kgds(20)=255    ! oct 5, used for thinned grids, set to 255

   else
      
     call error_handler("UNRECOGNIZED INPUT GRID TYPE ", 1)

 endif

 end subroutine gdt_to_gds

 end module model_grid