forcing-perturbation/interpolation/scrip_sources/scrip_use_extrap.f

!-----------------------------------------------------------------------
!
! This script interpolates a (time, latitude, longitude) netcdf file 
! using the interpolation weights computed by SCRIP :
! Spherical Coordinate Remapping and Interpolation Package
! and then extrapolates using the nearest neighbour method
!
! The arguments are passed through a namelist named scrip_use_in:
!&remap_inputs
!    remap_wgt   = 'Weights from SCRIP'
!    infile      = 'input netcdf file'
!    invertlat   = TRUE/FALSE : should the latitudes be reverted ?
!    var         = 'netcdf variable name'
!    fromregular = TRUE/FALSE : is the input grid regular ?
!    outfile     = 'output netcdf file'
!/
!
! History : Virginie Guemas - Initial version                   - 2011 
!-----------------------------------------------------------------------
      program scrip_use_extrap

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

      use kinds_mod    ! defines common data types
      use constants    ! defines common constants
      use iounits      ! I/O unit manager
      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
      use read_input_file ! routines to read file to interpolate

      implicit none

!-----------------------------------------------------------------------
!
!     input namelist variables
!
!-----------------------------------------------------------------------

      character (char_len) :: 
     &  remap_wgt,     ! filename containing remap data (map1)
     &  infile,        ! filename containing input field
     &  var,           ! var name in infile
     &  outfile        ! filename to output interpolated field
    
      logical :: fromregular

      namelist /remap_inputs/ remap_wgt, infile, invertlat, var,
     &                        fromregular, outfile

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

      character (char_len) :: 
     &        map_name      ! name for mapping from grid1 to grid2

      integer (kind=int_kind) ::    ! netCDF ids for files and arrays
     &        ncstat, nc_outfile_id, 
     &        nc_dstgrdcntrlat_id, nc_dstgrdcntrlon_id,
     &        nc_dstarray1_id, nc_dstarray1a_id, nc_dstarray2_id,
     &        nc_vartime_id, nc_srcarray_id

      integer (kind=int_kind), dimension(:), allocatable ::
     &        nc_grid2size_id, nc_grid1size_id

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

      character (char_len) :: 
     &          dim_name, attname    ! netCDF dimension name

      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,jatt,
     &    iunit,jtime            ! unit number for namelist file

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

      real (kind=dbl_kind) ::
     &    delew, delns,     ! variables for computing bicub gradients
     &    length,           ! length scale for cosine hill test field
     &    dist_min,distance ! Extrapolation

      real (kind=dbl_kind), dimension(:), allocatable ::
     &    grid1_tmp, 
     &    grad1_lat, 
     &    grad1_lon, 
     &    grad1_latlon, 
     &    grad1_lat_zero, 
     &    grad1_lon_zero, 
     &    grid2_array, 
     &    grid2_err,
     &    grid2_tmp,
     &    grid1_mask_grid2,
     &    dstcoslon,
     &    dstcoslat,
     &    srccoslon,
     &    srccoslat,
     &    dstsinlon,
     &    dstsinlat,
     &    srcsinlon,
     &    srcsinlat

!-----------------------------------------------------------------------
!
!     read namelist for file and mapping info
!
!-----------------------------------------------------------------------

      call get_unit(iunit)
      open(iunit, file='scrip_use_in', status='old', form='formatted')
      read(iunit, nml=remap_inputs)
      call release_unit(iunit)
      write(*,nml=remap_inputs)

!-----------------------------------------------------------------------
!
!     read remapping data
!
!-----------------------------------------------------------------------

      call read_remap(map_name, remap_wgt)

!-----------------------------------------------------------------------
!
!     read input file
!
!-----------------------------------------------------------------------

      call read_input(infile, var)

!-----------------------------------------------------------------------
!
!     allocate arrays
!
!-----------------------------------------------------------------------

      allocate (grid1_tmp       (grid1_size),
     &          grad1_lat       (grid1_size), 
     &          grad1_lon       (grid1_size), 
     &          grad1_latlon    (grid1_size),
     &          grad1_lat_zero  (grid1_size), 
     &          grad1_lon_zero  (grid1_size), 
     &          grid1_imask     (grid1_size),
     &          grid2_array     (grid2_size), 
     &          grid2_err       (grid2_size),
     &          grid2_tmp       (grid2_size),
     &          grid2_imask     (grid2_size),
     &          grid2_count     (grid2_size),
     &          grid1_mask_grid2(grid2_size))

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

      allocate (srccoslon       (grid1_size),
     &          srcsinlon       (grid1_size),
     &          srccoslat       (grid1_size),
     &          srcsinlat       (grid1_size),
     &          dstcoslon       (grid2_size),
     &          dstsinlon       (grid2_size),
     &          dstcoslat       (grid2_size),
     &          dstsinlat       (grid2_size),
     &          neighbour       (grid2_size)) 
 
      dstcoslon = cos(grid2_center_lon*pi/180)
      dstsinlon = sin(grid2_center_lon*pi/180)
      srccoslon = cos(grid1_center_lon*pi/180)
      srcsinlon = sin(grid1_center_lon*pi/180)

      dstcoslat = cos(grid2_center_lat*pi/180)
      dstsinlat = sin(grid2_center_lat*pi/180)
      srccoslat = cos(grid1_center_lat*pi/180)
      srcsinlat = sin(grid1_center_lat*pi/180)

!-----------------------------------------------------------------------
!
!     setup a NetCDF file for output
!
!-----------------------------------------------------------------------

      !***
      !*** create netCDF dataset 
      !***

      ncstat = nf_create (outfile, NF_CLOBBER, nc_outfile_id)
      call netcdf_error_handler(ncstat)

      ncstat = nf_put_att_text (nc_outfile_id, NF_GLOBAL, 'title',
     &                          len_trim(map_name), map_name)
      call netcdf_error_handler(ncstat)

      !***
      !*** define grid size dimensions
      !***

      allocate( nc_grid2size_id(grid2_rank+1))
      allocate( nc_grid1size_id(grid1_rank+1))
    
      ncstat = nf_def_dim (nc_outfile_id, 'x',
     &                       grid2_dims(1), nc_grid2size_id(1))
      call netcdf_error_handler(ncstat)

      ncstat = nf_def_dim (nc_outfile_id, 'y',
     &                       grid2_dims(2), nc_grid2size_id(2))
      call netcdf_error_handler(ncstat)

      ncstat = nf_def_dim (nc_outfile_id, 'inx',
     &                       grid1_dims(1), nc_grid1size_id(1))
      call netcdf_error_handler(ncstat)

      ncstat = nf_def_dim (nc_outfile_id, 'iny',
     &                       grid1_dims(2), nc_grid1size_id(2))
      call netcdf_error_handler(ncstat)

      !***
      !*** Create time dimension
      !***

      ncstat = nf_def_dim (nc_outfile_id, 'time',
     &                       NF_UNLIMITED,  nc_grid2size_id(3))
      call netcdf_error_handler(ncstat)
      nc_grid1size_id(3)=nc_grid2size_id(3)
      
      ncstat = nf_def_var (nc_outfile_id, 'time',
     &           NF_DOUBLE, 1, nc_grid2size_id(3), nc_vartime_id)
      call netcdf_error_handler(ncstat)

      if ( nc_time_id > -1 ) then
        if ( natts >= 1 ) then 
          do jatt = 1,natts
            ncstat = nf_inq_attname(nc_infile_id, nc_time_id, jatt,
     &                              attname)
            call netcdf_error_handler(ncstat)

            ncstat = nf_copy_att(nc_infile_id, nc_time_id, attname,
     &                       nc_outfile_id, nc_vartime_id)
          enddo
       endif
      endif

      !***
      !*** define grid center latitude array
      !***

      ncstat = nf_def_var (nc_outfile_id, 'latitude', 
     &                     NF_DOUBLE, grid2_rank, nc_grid2size_id 
     &                     (1:grid2_rank), nc_dstgrdcntrlat_id)
      call netcdf_error_handler(ncstat)

      ncstat = nf_put_att_text (nc_outfile_id, nc_dstgrdcntrlat_id, 
     &                          'units', 7, 'degrees')
      call netcdf_error_handler(ncstat)

      !***
      !*** define grid center longitude array
      !***

      ncstat = nf_def_var (nc_outfile_id, 'longitude', 
     &                     NF_DOUBLE, grid2_rank, nc_grid2size_id 
     &                     (1:grid2_rank), nc_dstgrdcntrlon_id)
      call netcdf_error_handler(ncstat)

      ncstat = nf_put_att_text (nc_outfile_id, nc_dstgrdcntrlon_id, 
     &                          'units', 7, 'degrees')
      call netcdf_error_handler(ncstat)

      !***
      !*** define source array
      !***

      ncstat = nf_def_var (nc_outfile_id, 'input',
     &                     NF_DOUBLE, (grid1_rank+1), nc_grid1size_id,
     &                     nc_srcarray_id)
      call netcdf_error_handler(ncstat)
      ncstat = nf_put_att_double (nc_outfile_id, nc_srcarray_id,
     &                     'missing_value', NF_DOUBLE,1,dble(1e20))
      call netcdf_error_handler(ncstat)

      !***
      !*** define destination arrays
      !***

      ncstat = nf_def_var (nc_outfile_id, var, 
     &                     NF_DOUBLE, ( grid2_rank + 1 ), 
     &                     nc_grid2size_id, nc_dstarray1_id )
      call netcdf_error_handler(ncstat)
      ncstat = nf_put_att_double (nc_outfile_id, nc_dstarray1_id,
     &                     'missing_value', NF_DOUBLE,1,dble(1e20))
      call netcdf_error_handler(ncstat)

      if ( nvaratts >= 1 ) then
        do jatt = 1,nvaratts
          ncstat = nf_inq_attname(nc_infile_id, nc_invar_id, jatt,
     &                            attname)
          call netcdf_error_handler(ncstat)
          
          if ((attname .ne. '_FillValue') .and. (attname .ne.
     &        'missing_value') ) then 
            ncstat = nf_copy_att(nc_infile_id, nc_invar_id, attname,
     &                         nc_outfile_id, nc_dstarray1_id)
            call netcdf_error_handler(ncstat)
          endif
        enddo
      endif

      if ( nglobatts >= 1 ) then
        do jatt = 1,nglobatts
          ncstat = nf_inq_attname(nc_infile_id, NF_GLOBAL, jatt,
     &                            attname)
          call netcdf_error_handler(ncstat)

          ncstat = nf_copy_att(nc_infile_id, NF_GLOBAL, attname,
     &                       nc_outfile_id, NF_GLOBAL)
          call netcdf_error_handler(ncstat)
        enddo
      endif


      ncstat = nf_close(nc_infile_id)
      call netcdf_error_handler(ncstat)

      !***
      !*** end definition stage
      !***

      ncstat = nf_enddef(nc_outfile_id)
      call netcdf_error_handler(ncstat)

!-----------------------------------------------------------------------
!
!     write some grid info
!
!-----------------------------------------------------------------------

      !***
      !*** write grid center latitude array
      !***

      ncstat = nf_put_var_double(nc_outfile_id, nc_dstgrdcntrlat_id,
     &                           grid2_center_lat*180./pi)
      call netcdf_error_handler(ncstat)

      !***
      !*** write grid center longitude array
      !***

      ncstat = nf_put_var_double(nc_outfile_id, nc_dstgrdcntrlon_id,
     &                           grid2_center_lon*180./pi)
      call netcdf_error_handler(ncstat)

!-----------------------------------------------------------------------
!
!     Interpolate the input mask
!
!-----------------------------------------------------------------------

      call remap(grid1_mask_grid2, wts_map1, grid2_add_map1,
     &           grid1_add_map1, real(grid1_imask,kind=dbl_kind))

!-----------------------------------------------------------------------
!
!     Prepare extrapolation 
!
!-----------------------------------------------------------------------

      neighbour = 0
      do n = 1,grid2_size
        if ( grid2_imask(n)<0.5 .or. grid1_mask_grid2(n) == 0. ) then
          dist_min     = 1e20

          do j = 1,grid1_size

            if ( grid1_imask(j) == 1 ) then
              distance  = acos(dstcoslat(n)*srccoslat(j)*
     &          (dstcoslon(n)*srccoslon(j)+dstsinlon(n)*srcsinlon(j))
     &          +dstsinlat(n)*srcsinlat(j))
              if ( distance < dist_min ) then
                dist_min     = distance
                neighbour(n) = j
              end if
            end if

          end do
 
        end if
      end do

!-----------------------------------------------------------------------
!
!     Write time dimension
!
!-----------------------------------------------------------------------

      do jtime = 1,ntime

        ncstat = nf_put_vara_double(nc_outfile_id, nc_vartime_id,
     &                              jtime, 1, time(jtime))
        call netcdf_error_handler(ncstat)
        
!-----------------------------------------------------------------------
!
!     if bicubic or 2nd-order conservative, 3 gradients needed in space
!
!-----------------------------------------------------------------------

      if ( fromregular .and. (map_type == map_type_bicubic .or.  
     &    map_type == map_type_conserv) ) then

        nx = grid1_dims(1)
        ny = grid1_dims(2)

        do n=1,grid1_size

          grad1_lat(n) = zero
          grad1_lon(n) = zero
          grad1_latlon(n) = zero

          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,jtime) - 
     &                            grid1_array(iw,jtime))

            !*** 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,jtime) -
     &                            grid1_array(is,jtime))

            !*** 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,jtime) -
     &                                 grid1_array(inw,jtime))

            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,jtime) -
     &                                 grid1_array(isw,jtime))

            grad1_latlon(n) = delns*(grad1_lat_zero(n) -
     &                               grad1_lon_zero(n))

          endif
        enddo
      endif

!-----------------------------------------------------------------------
!
!     remapping from grid1 to grid2
!
!-----------------------------------------------------------------------

      grad1_lat_zero = zero
      grad1_lon_zero = zero

      if (map_type == map_type_bicubic) then   
        if (fromregular) then 
          call remap(grid2_tmp, wts_map1, grid2_add_map1, 
     &               grid1_add_map1, grid1_array(:,jtime), 
     &                          src_grad1=grad1_lat,
     &                          src_grad2=grad1_lon,
     &                          src_grad3=grad1_latlon)
        else 
          print*,"Input grid is not regular, bicubic interpolation "
          stop"              is not possible : We stop "
        endif 
      elseif (map_type == map_type_conserv .and. fromregular ) then
        call remap(grid2_tmp,wts_map1,grid2_add_map1,grid1_add_map1,
     &             grid1_array(:,jtime), src_grad1=grad1_lat,
     &                          src_grad2=grad1_lon)
      else    
        call remap(grid2_tmp, wts_map1, grid2_add_map1, grid1_add_map1,
     &             grid1_array(:,jtime))
      endif

      if (map_type == map_type_conserv) then
        select case (norm_opt)
        case (norm_opt_none)
          grid2_err = grid2_frac*grid2_area
          where (grid2_err /= zero)
            grid2_tmp = grid2_tmp/grid2_err
          else where
            grid2_tmp = zero
          end where
        case (norm_opt_frcarea)
        case (norm_opt_dstarea)
          where (grid2_frac /= zero)
            grid2_tmp = grid2_tmp/grid2_frac
          else where
            grid2_tmp = zero
          end where
        end select
      end if

!-----------------------------------------------------------------------
!
!     extrapolation
!
!-----------------------------------------------------------------------
!
      do n = 1,grid2_size
        if ( grid2_imask(n)<0.5 .or. grid1_mask_grid2(n) == 0. ) then

          if ( neighbour(n) > 0 ) then
            grid2_tmp(n) = grid1_array(neighbour(n),jtime)
          else
            grid2_tmp(n) = 1e20
          end if
 
        end if
      end do

!-----------------------------------------------------------------------
!
!     write results to NetCDF file
!
!-----------------------------------------------------------------------

      ncstat = nf_put_vara_double(nc_outfile_id, nc_dstarray1_id,
     &         (/1, 1, jtime/), (/grid2_dims, 1/), grid2_tmp  )
      call netcdf_error_handler(ncstat)

      where (grid1_imask<0.5)
        grid1_array(:,jtime)=1e20
      end where

      ncstat = nf_put_vara_double(nc_outfile_id, nc_srcarray_id,
     &         (/1, 1, jtime/), (/grid1_dims, 1/), grid1_array(:,jtime))
      call netcdf_error_handler(ncstat)

      enddo

!-----------------------------------------------------------------------
!
!     close netCDF file
!
!-----------------------------------------------------------------------

      ncstat = nf_close(nc_outfile_id)
      call netcdf_error_handler(ncstat)

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

      end program scrip_use_extrap

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!