ecearth3/sources/nemo-3.6/TOOLS/WEIGHTS/nocsutil/scripinterp_mod.F90

! **************************************************************************

  module scripinterp_mod

! ==========================================================================

  use kinds_mod    ! defines common data types
  use constants    ! defines common constants
  use iounits      ! I/O unit manager
  use netcdf
  use netcdf_mod   ! netcdf I/O stuff
  use grids        ! module containing grid info
  use remap_vars   ! module containing remapping info
  use remap_mod    ! module containing remapping routines
  use remap_read   ! routines for reading remap files

  implicit none

  real(kind=dbl_kind), dimension(:), allocatable ::   &
    grid_out

  integer(kind=int_kind) ::                           &  ! netCDF ids for files and arrays
          ncstat, nc_outfile_id, nc_infile_id
  integer (kind=int_kind), dimension(4) ::    &
      input_dims, input_dimids, input_count
  real (kind=dbl_kind), dimension(:), allocatable ::   &
      scale
  integer (kind=int_kind), dimension(:), allocatable ::        &
      nc_xysize_id, nc_gridsize_id, nc_gridsize,           &
      nc_variable_id
  integer :: nc_lon_id, nc_lat_id, nc_array_id

  character (char_len) ::  &
          map_name      ! name for mapping from grid1 to grid2
  character (1), parameter ::  &
          separator = '|'
  
  ! - input namelist variables
  
  character (char_len) ::  &
    input_file,            &  ! filename containing input fields
    interp_file,           &  ! filename containing remap data (map1)
    input_name                ! name of variable to grid
  integer (kind=int_kind), dimension(4) :: &
    input_stride,  &          ! how much of input array to process
    input_start,   &          ! where to start
    input_stop                ! and where to stop
  character (char_len), dimension(4) ::  &
    input_vars                ! input variables to be copied
  
  ! - output namelist variables
  
  character (char_len) ::  &
    output_file,           &  ! filename for test results
    output_mode,           &  ! 'create' or 'append'
    output_name,           &  ! name of new grid variable
    output_lat,            &  ! as it says
    output_lon,            &  ! as it says
    output_ydir               ! choose to invert output arrays in y direction
  character (char_len), dimension(4) ::  &
    output_dims,           &    ! name of new grid variable
    output_vars                 ! variables to copy
  character (char_len), dimension(10) ::  &
    output_attributes,     &  ! attributes of stuff in output file
    output_scaling             ! scaling factor to apply to input data to get
                               ! correct units in output file
contains

! ==========================================================================

  subroutine process_grid(nm_in)
  
    !-----------------------------------------------------------------------
    ! - dummy variables

    character (char_len) ::  &
      nm_in                    ! name of input namelist file
  
    !-----------------------------------------------------------------------
    ! - local variables

    integer (kind=int_kind), dimension(4) ::    &
        astart, acount, plus_one
    integer (kind=int_kind), dimension(3) ::    &
        write_dims
    integer (kind=int_kind) ::                  &
      i1, i2, jdim, n, nx, ny, nloop,           &
      nc_input_id, nc_input_rank,               &
      vstart, vstride, numv

    real (kind=dbl_kind), dimension(:), allocatable  ::    &
      grid1_array
    real (kind=dbl_kind), dimension(:,:), allocatable  ::    &
      var_out

    plus_one(:) = 1

    !-----------------------------------------------------------------------

    call read_mappings(nm_in)

    !-----------------------------------------------------------------------
    ! - read input grid
    ! - WARNING - lots of assumptions here at the moment
  
    ncstat = nf90_open( input_file, NF90_NOWRITE, nc_infile_id )
    call netcdf_error_handler(ncstat,"open")
  
    ncstat = nf90_inq_varid( nc_infile_id, input_name, nc_input_id )
    call netcdf_error_handler(ncstat,"inq_varid")
  
    input_dims(:) = 0
    ncstat = nf90_inquire_variable( nc_infile_id, nc_input_id,   &
                  ndims=nc_input_rank, dimids=input_dimids(:) )
    call netcdf_error_handler(ncstat,"inquire_variable")

    do jdim = 1,nc_input_rank
      ncstat = nf90_inquire_dimension(nc_infile_id,              &
                    input_dimids(jdim), len=input_dims(jdim) )
      call netcdf_error_handler(ncstat,"inquire_dimension")
    enddo
  
    ! - dimids seem to be returned in storage order so the outer dimension of 
    ! - the array as described by the netcdf file becomes the first dimension
    ! - as far as f90 is concerned (confused? you will be!)

    do jdim = 1,nc_input_rank
      if (input_stop(jdim) == 0) then
        input_stop(jdim) = input_dims(jdim)
      endif
      input_count(jdim) = input_stop(jdim) - input_start(jdim) + 1
    enddo

    ! - rashly we assume x followed by y
    nx = input_dims(1)
    ny = input_dims(2)
    write(*,fmt='("Input grid dimensions are:",2i6)') nx, ny
    if (nx*ny /= grid1_size) then
      write(6,*) "mismatch between input grid and remap data"
      stop
    endif

    ! - calculate number of horizontal slices to process
    ! - at the moment this is not very general and will only work with 3 dimensions

    acount(1:nc_input_rank) = &
             (input_stop(1:nc_input_rank)-input_start(1:nc_input_rank)+1) / &
             input_stride(1:nc_input_rank)
    nloop = 1
    do jdim = 1,nc_input_rank
      nloop = nloop*acount(jdim)
    enddo
    nloop = nloop/grid1_size
    write(6,*) "total slices requested: ",nloop

    vstart = input_start(nc_input_rank)            ! ie extra var has outer dimension
    vstride = input_stride(nc_input_rank)

    ! - in general we cant read in the whole array so do it slice by slice
    ! - slow but sure

    write(6,*) "allocating input and output grids"
    allocate( grid1_array(grid1_size))
    allocate( grid_out(grid2_size) )

    numv = 0
    do n = 1,4
      if (trim(input_vars(n)) /= '-' .and. &
          trim(output_vars(n)) /= '-') numv = numv + 1
    enddo

    write_dims(1) = grid2_dims(1)
    write_dims(2) = grid2_dims(2)
    write_dims(3) = nloop
    call define_grid(write_dims(1:3) , 2+numv)

    astart(:) = input_start(:)
    astart(3) = astart(3) - input_stride(3)
    acount(:) = 1
    acount(1) = nx
    acount(2) = ny

    do n = 1,nloop

      write(6,*) "processing slice: ",n
      astart(3) = astart(3) + input_stride(3)
      ncstat = nf90_get_var(nc_infile_id, nc_input_id, grid1_array,  &
                            start=astart(1:nc_input_rank),              &
                            count=acount(1:nc_input_rank))
      call netcdf_error_handler(ncstat,"get_var")

      call calculate_grid(grid1_array, grid_out)

      call write_grid(grid_out, n, write_dims(1:2) , 2)

    enddo

    ! ---------------------------------------------------------------------
    ! - now for any extra variables to copy

    if (numv > 0) then

      write(6,*) "reading ",numv," extra variables"
      allocate( var_out(nloop,numv) )
  
      do n = 1,numv
        write(6,*) "looking for variable: ",trim(input_vars(n))
        ncstat = nf90_inq_varid( nc_infile_id, trim(input_vars(n)), nc_input_id )
        call netcdf_error_handler(ncstat,"inq_varid")
      
        input_dims(:) = 0
        ncstat = nf90_inquire_variable( nc_infile_id, nc_input_id,   &
                      ndims=nc_input_rank, dimids=input_dimids(:) )
        call netcdf_error_handler(ncstat,"inquire_variable")
  
        if (nc_input_rank /= 1) then
          write(6,*) 'sorry, only rank 1 variables can be copied'
          cycle
        endif
        ncstat = nf90_inquire_dimension(nc_infile_id,              &
                      input_dimids(1), len=input_dims(1) )
        call netcdf_error_handler(ncstat,"inquire_dimension")
  
        ncstat = nf90_get_var(nc_infile_id, nc_input_id, var_out(1:nloop,n),  &
                              start=(/ vstart /), stride=(/ vstride /))
        call netcdf_error_handler(ncstat,"get_var")
      enddo
  
      call write_extra(var_out, numv+2)
      deallocate(var_out)

    endif

    ncstat = nf90_close(nc_outfile_id)
    call netcdf_error_handler(ncstat,"out close")
    ncstat = nf90_close(nc_infile_id)
    call netcdf_error_handler(ncstat,"in close")

    ! ---------------------------------------------------------------------

    deallocate( grid1_array, grid_out)

    ! ---------------------------------------------------------------------

  end subroutine process_grid

  ! ==========================================================================

  subroutine define_grid(thedims, therank)

    !-----------------------------------------------------------------------
    ! - dummy variables
    
    integer (kind=int_kind) :: &
      therank
    integer (kind=int_kind), dimension(therank) :: &
      thedims

    !-----------------------------------------------------------------------
    ! - local variables

    integer ::               &
      k, n, ilon, ilat, icolon, i1, i2, natt, nvar, id, jd, kd, nd
    character (char_len) ::  &
      aname, vname, att
    real (kind=dbl_kind) ::   s

    ! - netcdf variables

    integer :: xtype

    !-----------------------------------------------------------------------
    ! -  define grid size dimensions
  
    allocate(nc_xysize_id(grid2_rank))
    allocate(nc_gridsize_id(therank))
    allocate(nc_gridsize(therank))
    allocate(nc_variable_id(therank-2))
    
    !-----------------------------------------------------------------------
    ! - setup a NetCDF file for output
  
    xtype = NF90_FLOAT

    write(6,*) 'creating output file'
    ncstat = nf90_create (output_file, NF90_CLOBBER, nc_outfile_id)
    call netcdf_error_handler(ncstat,"create")
  
    write(6,*) 'setting global attributes'
    ncstat = nf90_put_att(nc_outfile_id, NF90_GLOBAL, 'title', map_name)
    call netcdf_error_handler(ncstat,"put_att")
  
    write(6,*) 'setting dimensions'
    do n=1,therank
      if (n .eq. therank .and. therank .gt. 2) then
        write(6,*) '    unlimited dim ',trim(output_dims(n)),' size: ',thedims(n)
        ncstat = nf90_def_dim (nc_outfile_id, output_dims(n), NF90_UNLIMITED, &
                               nc_gridsize_id(n))
      else
        write(6,*) '    dim ',trim(output_dims(n)),' size: ',thedims(n)
        ncstat = nf90_def_dim (nc_outfile_id, output_dims(n), thedims(n), &
                               nc_gridsize_id(n))
      endif
      call netcdf_error_handler(ncstat,"def_dim")
    end do
    nc_gridsize(:) = thedims(1:therank)

    ! - at the moment there is an assumption here that the ordering is (lon,lat)

    ilon = 1
    ilat = 2
    nc_xysize_id(1) = nc_gridsize_id(ilon)
    nc_xysize_id(2) = nc_gridsize_id(ilat)

    ! ----------------------------------------------------------------
    ! -  define grid center longitude array
  
    write(6,*) 'defining longitude variable'
    ncstat = nf90_def_var (nc_outfile_id, output_lon,   &
                         xtype, nc_xysize_id,   &
                         nc_lon_id)
    call netcdf_error_handler(ncstat,"def_var")
  
    ncstat = nf90_put_att (nc_outfile_id, nc_lon_id, 'units', 'degrees')
    call netcdf_error_handler(ncstat,"put_att")
  
    ! ----------------------------------------------------------------
    ! -  define grid center latitude array
  
    write(6,*) 'defining latitude variable'
    ncstat = nf90_def_var (nc_outfile_id, output_lat,     &
                         xtype, nc_xysize_id,   &
                         nc_lat_id)
    call netcdf_error_handler(ncstat,"def_var")
  
    ncstat = nf90_put_att (nc_outfile_id, nc_lat_id, 'units', 'degrees')
    call netcdf_error_handler(ncstat,"put_att")
  
    ! ----------------------------------------------------------------
    ! -  define copy variables array
  
    write(6,*) 'defining copy variables'
    do n = 3,therank
      ncstat = nf90_def_var (nc_outfile_id, output_vars(n-2),   &
                           xtype, nc_gridsize_id(n),    &
                           nc_variable_id(n-2))
      call netcdf_error_handler(ncstat,"def_var")
    enddo
  
    ! ----------------------------------------------------------------
    ! -  define output array
  
    write(6,*) 'defining grid variable'
    ncstat = nf90_def_var (nc_outfile_id, output_name,               &
                         xtype, nc_gridsize_id,    &
                         nc_array_id)
    call netcdf_error_handler(ncstat,"def_var")
  
    ! ----------------------------------------------------------------
    ! - output attributes has to come after all other definitions
    ! - this code currently a bit murky, needs a rewrite
  
    ncstat = nf90_inquire (nc_outfile_id, nVariables=nvar)
    call netcdf_error_handler(ncstat,"inquire")
    do n = 1,10
      att = trim(output_attributes(n))
      natt = len(att)
      if (att /= '-') then
        i1 = index(att,separator)
        aname = att(1:i1-1)
        do k = 1,nvar
          ncstat = nf90_inquire_variable(nc_outfile_id, k, vname)
          call netcdf_error_handler(ncstat,"inquire_variable")
          if (vname == aname) then
            i2 = index(att,separator,.true.)
            ncstat = nf90_put_att (nc_outfile_id, k, &
                                   att(i1+1:i2-1), att(i2+1:natt))
            call netcdf_error_handler(ncstat,"put_att")
            exit      ! from innermost do
          endif
        enddo
      endif
    enddo
  
    ! output scaling

    allocate (scale(nvar))
    scale(:) = 1.0

    do n = 1,10
      att = trim(output_scaling(n))
      natt = len(att)
      if (att /= '-') then
        i1 = index(att,separator)
        aname = att(1:i1-1)
        do k = 1,nvar
          ncstat = nf90_inquire_variable(nc_outfile_id, k, vname)
          call netcdf_error_handler(ncstat,"inquire_variable")
          if (vname == aname) then
            i2 = index(att,separator,.true.)
            read(att(i2+1:natt),*) scale(k)
            call netcdf_error_handler(ncstat,"put_att")
            exit      ! from innermost do
          endif
        enddo
      endif
    enddo
  
    ! ----------------------------------------------------------------
    ! -  end definition stage
  
    ncstat = nf90_enddef(nc_outfile_id)
    call netcdf_error_handler(ncstat,"enddef")
 
  end subroutine define_grid

  ! ==========================================================================

  subroutine write_grid(thegrid, thelevel, thedims, therank)

    !-----------------------------------------------------------------------
    ! - dummy variables
    
    integer (kind=int_kind), intent(in) :: &
      therank, thelevel
    real (kind=dbl_kind), dimension(:), intent(in) ::   &
      thegrid
    integer (kind=int_kind), dimension(therank) :: &
      thedims

    !-----------------------------------------------------------------------
    ! - local variables

    integer ::               &
      k, n, ilon, ilat, icolon, j1, j2, dj, natt, nvar, id, jd, kd, nd
    character (char_len) ::  &
      aname, vname, att
    real (kind=dbl_kind), dimension(:,:,:), allocatable ::   &
      data
    real (kind=dbl_kind) ::   s
    real (kind=dbl_kind), parameter :: todeg = 57.295779513082323
    integer (kind=int_kind), dimension(3) ::   &
      start

    ! - netcdf variables

    integer :: xtype

    !-----------------------------------------------------------------------
    ! - write results to NetCDF file
  
    allocate (data(thedims(1),thedims(2),1))
    if (output_ydir .eq. 'invert') then
      j1 = thedims(2)
      j2 = 1
      dj = -1
    else
      j1 = 1
      j2 = thedims(2)
      dj = 1
    endif

    if (thelevel .eq. 1) then

      ! -  grid center latitude array
  
      write(6,*) 'writing latitude variable'
      s = scale(nc_lat_id)
      nd = 0
      do jd = j1,j2,dj
        do id =1,thedims(1)
          nd = nd + 1
          data(id,jd,1) = s*todeg*grid2_center_lat(nd)
        enddo
      enddo
      ncstat = nf90_put_var(nc_outfile_id, nc_lat_id, data(:,:,1))
      call netcdf_error_handler(ncstat,"put_var")
  
      ! -  grid center longitude array
  
      write(6,*) 'writing longitude variable'
      s = scale(nc_lon_id)
      nd = 0
      do jd = j1,j2,dj
        do id =1,thedims(1)
          nd = nd + 1
          data(id,jd,1) = s*todeg*grid2_center_lon(nd)
        enddo
      enddo
      ncstat = nf90_put_var(nc_outfile_id, nc_lon_id, data(:,:,1))
      call netcdf_error_handler(ncstat,"put_var")

    endif

    !-----------------------------------------------------------------------
    ! -  new grid

    write(6,*) 'writing grid variable'
    n = therank
    s = scale(nc_array_id)
    nd = 0
    do jd = j1,j2,dj
      do id =1,thedims(1)
        nd = nd + 1
        data(id,jd,1) = thegrid(nd)
      enddo
    enddo
    write(6,*) 'scaling data '
    data(:,:,1) = s*data(:,:,1)
    start(:) = (/ 1, 1, thelevel /)
    ncstat = nf90_put_var(nc_outfile_id, nc_array_id, data, start)
    call netcdf_error_handler(ncstat,"put_var")
    deallocate(data)

  end subroutine write_grid

  ! ==========================================================================

  subroutine write_extra(thevars, therank)

    real (kind=dbl_kind), dimension(:,:), intent(in) ::   &
      thevars
    real (kind=dbl_kind), dimension(:), allocatable ::   &
      thedata
    integer (kind=int_kind), intent(in) :: &
      therank
    real (kind=dbl_kind) ::   s
    integer :: n

    allocate( thedata(size(thevars,1)) )

    ! -  copy variable arrays

    write(6,*) 'writing copy variables'
    do n = 3,therank
      s = scale(nc_variable_id(n-2))
      thedata(:) = s*thevars(:,n-2)
      ncstat = nf90_put_var(nc_outfile_id, nc_variable_id(n-2), thedata)
      call netcdf_error_handler(ncstat,"put_var")
    enddo

    deallocate( thedata )

  end subroutine write_extra

  ! ==========================================================================

  subroutine close_grid()

  !     close netCDF file
  
    write(6,*) 'closing file'
    ncstat = nf90_close(nc_outfile_id)
    call netcdf_error_handler(ncstat,"close")
  
  end subroutine close_grid

  ! ==========================================================================

  subroutine read_mappings(nm_in)

    !-----------------------------------------------------------------------
    ! - dummy variables

    character (char_len) ::  &
      nm_in                    ! name of input namelist file
  
    !-----------------------------------------------------------------------
    ! - local variables
  
    character (char_len) ::   &
              dim_name    ! netCDF dimension name
  
    integer (kind=int_kind) ::     &
        iunit                                                ! unit number for namelist file
  
    !-----------------------------------------------------------------------
    ! - namelist block

    namelist /interp_inputs/ input_file, interp_file, input_name, &
                             input_stride, input_start, input_stop, &
                             input_vars
    namelist /interp_outputs/ output_dims, output_file, output_mode, output_name, &
                              output_lat, output_lon, output_ydir,  &
                              output_scaling, output_vars, output_attributes
  
    !-----------------------------------------------------------------------
    ! - read namelist for file and mapping info
  
    input_stride(:) = 1
    input_start(:) = 1
    input_stop(:) = 0
    output_scaling(:) = '-'
    input_vars(:) = '-'
    output_lon = '-'
    output_lat = '-'
    output_vars(:) = '-'
    output_ydir = 'none'
    output_attributes(:) = '-'

    call get_unit(iunit)
    open(iunit, file=nm_in, status='old', form='formatted')
    read(iunit, nml=interp_inputs)
    read(iunit, nml=interp_outputs)
    call release_unit(iunit)
    write(*,nml=interp_inputs)
    write(*,nml=interp_outputs)
    if (trim(output_mode) == "create") then
      if (trim(output_lon) == '-' .or. trim(output_lat) == '-') then
        write(6,*) 'if creating, need to supply lon and lat names'
        stop
      endif
    endif

    !-----------------------------------------------------------------------
    ! - read remapping data
    ! - via the scrip package this sets variables:
    !    grid1_size, grid2_size:     sizes of input and output grids
    !    grid1_mask, grid2_mask:     masks
    !    grid1_rank, grid2_rank:     ranks
  
    call read_remap(map_name, interp_file)
  
  end subroutine read_mappings

  ! ==========================================================================

  subroutine calculate_grid(grid1_array, grid2_array)

    !-----------------------------------------------------------------------
    ! - dummy variables

    real (kind=dbl_kind), intent(in), dimension(:) ::    &
      grid1_array
    real (kind=dbl_kind), intent(out), dimension(:) ::    &
      grid2_array

    !-----------------------------------------------------------------------
    ! - local variables

    integer (kind=int_kind), dimension(:), allocatable :: &
        grid1_imask, grid2_imask, grid2_count

    real (kind=dbl_kind), dimension(:), allocatable :: &
        grid1_tmp,           &
        grad1_lat,           &
        grad1_lon,           &
        grad1_latlon,        &
        grad1_lat_zero,      &
        grad1_lon_zero,      &
        grid2_tmp1,           &
        grid2_tmp2

    real (kind=dbl_kind) ::  &
        delew, delns            ! variables for computing bicub gradients

    integer (kind=int_kind) ::                          &
        i,j,n,imin,imax,idiff,                          &
        ip1,im1,jp1,jm1,nx,ny,                          &  ! for computing bicub gradients
        in,is,ie,iw,ine,inw,ise,isw

    logical, parameter :: lat_gradient = .false.
  
    write(6,*) 'starting'

    !-----------------------------------------------------------------------
    ! - allocate arrays
  
    allocate (grid1_tmp      (grid1_size),   &
              grad1_lat      (grid1_size),    &
              grad1_lon      (grid1_size),    &
              grad1_lat_zero (grid1_size),    &
              grad1_lon_zero (grid1_size),    &
              grid1_imask    (grid1_size),   &
              grid2_tmp1     (grid2_size),   &
              grid2_tmp2     (grid2_size),   &
              grid2_imask    (grid2_size),   &
              grid2_count    (grid2_size))
  
    write(6,*) 'allocated'
    write(6,*) grid1_size,grid2_size

    grid1_imask(:) = 1
    grid2_imask(:) = 1
    where (grid1_mask)
      grid1_imask = 1
    elsewhere
      grid1_imask = 0
    endwhere
    where (grid2_mask)
      grid2_imask = 1
    elsewhere
      grid2_imask = 0
    endwhere

    write(6,*) 'masked'

    grad1_lat_zero = zero
    grad1_lon_zero = zero
    nx = input_dims(1)
    ny = input_dims(2)
    write(6,*) nx,ny
  
    !-----------------------------------------------------------------------
    ! - if bicubic, we need 3 gradients in logical space
  
    if (map_type == map_type_bicubic) then
  
      write(6,*) 'bicubic'
      write(6,*) grid1_size

      allocate (grad1_latlon (grid1_size)) 
  
      do n=1,grid1_size
  
        grad1_lat(n) = zero
        grad1_lon(n) = zero
        grad1_latlon(n) = zero
  
!       if (n.ge.8000) write(6,*) 0,grid1_mask(n),nx
        if (grid1_mask(n)) then
  
          delew = half
          delns = half
  
          j = (n-1)/nx + 1
          i = n - (j-1)*nx
  
          ip1 = i+1
          im1 = i-1
          jp1 = j+1
          jm1 = j-1
  
          if (ip1 > nx) ip1 = ip1 - nx
          if (im1 < 1 ) im1 = nx
          if (jp1 > ny) then
            jp1 = j
            delns = one
          endif
          if (jm1 < 1 ) then
            jm1 = j
            delns = one
          endif
  
          in  = (jp1-1)*nx + i
          is  = (jm1-1)*nx + i
          ie  = (j  -1)*nx + ip1
          iw  = (j  -1)*nx + im1
  
          ine = (jp1-1)*nx + ip1
          inw = (jp1-1)*nx + im1
          ise = (jm1-1)*nx + ip1
          isw = (jm1-1)*nx + im1
  
          ! -  compute i-gradient
  
          if (.not. grid1_mask(ie)) then
            ie = n
            delew = one
          endif
          if (.not. grid1_mask(iw)) then
            iw = n
            delew = one
          endif
   
          grad1_lat(n) = delew*(grid1_array(ie) - grid1_array(iw))
!         if (n.ge.8000) write(6,*) 1,grad1_lat(n)
  
          ! -  compute j-gradient
  
          if (.not. grid1_mask(in)) then
            in = n
            delns = one
          endif
          if (.not. grid1_mask(is)) then
            is = n
            delns = one
          endif
   
          grad1_lon(n) = delns*(grid1_array(in) - grid1_array(is))
!         if (n.ge.8000) write(6,*) 2,grad1_lon(n)
  
          ! -  compute ij-gradient
  
          delew = half
          if (jp1 == j .or. jm1 == j) then
            delns = one
          else 
            delns = half
          endif
  
          if (.not. grid1_mask(ine)) then
            if (in /= n) then
              ine = in
              delew = one
            else if (ie /= n) then
              ine = ie
              inw = iw
              if (inw == n) delew = one
              delns = one
            else
              ine = n
              inw = iw
              delew = one
              delns = one
            endif
          endif
  
          if (.not. grid1_mask(inw)) then
            if (in /= n) then
              inw = in
              delew = one
            else if (iw /= n) then
              inw = iw
              ine = ie
              if (ie == n) delew = one
              delns = one
            else
              inw = n
              ine = ie
              delew = one
              delns = one
            endif
          endif
  
          grad1_lat_zero(n) = delew*(grid1_array(ine) - grid1_array(inw))
!         if (n.ge.8000) write(6,*) 3,grad1_lat_zero(n)
  
          if (.not. grid1_mask(ise)) then
            if (is /= n) then
              ise = is
              delew = one
            else if (ie /= n) then
              ise = ie
              isw = iw
              if (isw == n) delew = one
              delns = one
            else
              ise = n
              isw = iw
              delew = one
              delns = one
            endif
          endif
  
          if (.not. grid1_mask(isw)) then
            if (is /= n) then
              isw = is
              delew = one
            else if (iw /= n) then
              isw = iw
              ise = ie
              if (ie == n) delew = one
              delns = one
            else
              isw = n
              ise = ie
              delew = one
              delns = one
            endif
          endif
  
          grad1_lon_zero(n) = delew*(grid1_array(ise) - grid1_array(isw))
          grad1_latlon(n) = delns*(grad1_lat_zero(n) - grad1_lon_zero(n))
!         if (n.ge.8000) write(6,*) 4,grad1_lon_zero(n),grad1_latlon(n)
  
        endif
      enddo
  
      write(6,*) 'remapping'
      call remap(grid2_array, wts_map1, grid2_add_map1, grid1_add_map1, &
                 grid1_array, src_grad1=grad1_lat,                    &
                 src_grad2=grad1_lon, src_grad3=grad1_latlon)
  
      print *,'Third order mapping from grid1 to grid2:'
      print *,'----------------------------------------'
      print *,'Grid1 min,max: ',minval(grid1_array),maxval(grid1_array)
      print *,'Grid2 min,max: ',minval(grid2_array  ),maxval(grid2_array  )
    
      ! -  Conservation Test
    
      print *,'Conservation:'
      print *,'Grid1 Integral = ',sum(grid1_array*grid1_area*grid1_frac)
      print *,'Grid2 Integral = ',sum(grid2_array  *grid2_area*grid2_frac)
    
    !-----------------------------------------------------------------------
    ! - a first-order map from grid1 to grid2
  
    else if (map_type /= map_type_conserv .AND.map_type /= map_type_bicubic) then

      write(6,*) 'bilinear or conservative'
  
      call remap(grid2_array, wts_map1, grid2_add_map1, grid1_add_map1,grid1_array)
  
      print *,'First order mapping from grid1 to grid2:'
      print *,'----------------------------------------'
      print *,'Grid1 min,max: ',minval(grid1_array),maxval(grid1_array)
      print *,'Grid2 min,max: ',minval(grid2_array  ),maxval(grid2_array  )
    
      ! -  Conservation Test
    
      print *,'Conservation:'
      print *,'Grid1 Integral = ',sum(grid1_array*grid1_area*grid1_frac)
      print *,'Grid2 Integral = ',sum(grid2_array  *grid2_area*grid2_frac)
    
    !-----------------------------------------------------------------------
    ! - conservative mappings:
    ! - a second-order map from grid1 to grid2 with only lat grads
  
    else if (map_type == map_type_conserv .AND. lat_gradient) then
  
      call remap(grid2_array, wts_map1, grid2_add_map1, grid1_add_map1,  &
                 grid1_array, src_grad1=grad1_lat,src_grad2=grad1_lon_zero) 
  
      select case (norm_opt)
      case (norm_opt_none)
        grid2_tmp2 = grid2_frac*grid2_area
        where (grid2_tmp2 /= zero)
          grid2_array = grid2_array/grid2_tmp2
        elsewhere
          grid2_array = zero
        end where
      case (norm_opt_frcarea)
      case (norm_opt_dstarea)
        where (grid2_frac /= zero)
          grid2_array = grid2_array/grid2_frac
        elsewhere
          grid2_array = zero
        end where
      end select
  
      print *,'Second order mapping from grid1 to grid2 (lat only):'
      print *,'----------------------------------------'
      print *,'Grid1 min,max: ',minval(grid1_array),maxval(grid1_array)
      print *,'Grid2 min,max: ',minval(grid2_array  ),maxval(grid2_array  )
  
      ! -  Conservation Test
  
      print *,'Conservation:'
      print *,'Grid1 Integral = ',sum(grid1_array*grid1_area*grid1_frac)
      print *,'Grid2 Integral = ',sum(grid2_array*grid2_area*grid2_frac)
  
    !-----------------------------------------------------------------------
    ! - conservative mappings:
    ! - a second-order map from grid1 to grid2 both gradients
  
    else if (map_type == map_type_conserv .AND..NOT. lat_gradient) then
  
      call remap(grid2_array,wts_map1,grid2_add_map1,grid1_add_map1,  &
                 grid1_array, src_grad1=grad1_lat,src_grad2=grad1_lon) 
  
      select case (norm_opt)
      case (norm_opt_none)
        grid2_tmp2 = grid2_frac*grid2_area
        where (grid2_tmp2 /= zero)
          grid2_array = grid2_array/grid2_tmp2
        elsewhere
          grid2_array = zero
        end where
      case (norm_opt_frcarea)
      case (norm_opt_dstarea)
        where (grid2_frac /= zero)
          grid2_array = grid2_array/grid2_frac
        elsewhere
          grid2_array = zero
        end where
      end select
  
      print *,'Second order mapping from grid1 to grid2:'
      print *,'-----------------------------------------'
      print *,'Grid1 min,max: ',minval(grid1_array),maxval(grid1_array)
      print *,'Grid2 min,max: ',minval(grid2_array  ),maxval(grid2_array  )
  
      ! -  Conservation Test
  
      print *,'Conservation:'
      print *,'Grid1 Integral = ',sum(grid1_array*grid1_area*grid1_frac)
      print *,'Grid2 Integral = ',sum(grid2_array*grid2_area*grid2_frac)
  
    endif
  
    !-----------------------------------------------------------------------
    !     calculate some statistics
  
    grid2_count = zero
    grid2_tmp1  = zero
    grid2_tmp2  = zero
  
    print *,'number of sparse matrix entries ',num_links_map1
    do n=1,num_links_map1
      grid2_count(grid2_add_map1(n)) = grid2_count(grid2_add_map1(n)) + 1
      if (wts_map1(1,n) > one .or. wts_map1(1,n) < zero) then
        grid2_tmp1(grid2_add_map1(n)) = grid2_tmp1(grid2_add_map1(n)) + 1
        grid2_tmp2(grid2_add_map1(n)) = max(abs(wts_map1(1,n)),grid2_tmp2(grid2_add_map1(n)) )
      endif
    end do
  
    do n=1,grid2_size
      if (grid2_tmp1(n) > zero) print *,n,grid2_tmp2(n)
    end do
  
    imin = minval(grid2_count, mask=(grid2_count > 0))
    imax = maxval(grid2_count)
    idiff =  (imax - imin)/10 + 1
    print *,'total number of dest cells ',grid2_size
    print *,'number of cells participating in remap ',count(grid2_count > zero)
    print *,'min no of entries/row = ',imin
    print *,'max no of entries/row = ',imax
  
    imax = imin + idiff
    do n=1,10
      print *,'num of rows with entries between ',imin,' - ',imax-1,   &
               count(grid2_count >= imin .and. grid2_count < imax)
      imin = imin + idiff
      imax = imax + idiff
    end do
  
    !-----------------------------------------------------------------------
    ! - deallocate arrays
  
    deallocate (grid1_tmp, grad1_lat, grad1_lon,    &
                grad1_lat_zero, grad1_lon_zero, grid1_imask,     &
                grid2_tmp1, grid2_tmp2,   &
                grid2_imask, grid2_count)
  
  end subroutine calculate_grid

  ! ==========================================================================

end module scripinterp_mod