!----------------------------------------------------------------------- ! ! 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 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!