shave_nc.F90

Go to the documentation of this file.

      program shave_nc
      use netcdf
      implicit none
      integer,parameter :: kdbl=selected_real_kind(p=13,r=200)
      character(len=255) :: filename_full ! the netcdf input file
      character(len=255) :: filename_shaved ! the NetCDF file for shaved data
      integer :: idim_compute,jdim_compute,halo
      integer :: i_count_compute,j_count_compute                        &
                ,i_count_super,j_count_super
      integer :: i_start ! starting i of 2-D data with halo rows
      integer :: j_start ! starting j of 2-D data with halo rows
      integer :: i_count ! i extent of 2-D data with halo rows
      integer :: j_count ! j extent of 2-D data with halo rows
      integer :: n_count
      integer :: n_shave ! number of rows to shave off full data to leave halo rows
      integer :: n_start ! location where Start reading incoming character data

      integer :: n,na,ncid_in,ncid_out,nd,ndims,ngatts     &
                ,nvars,unlimdimid
      integer :: natts ! # of attributes
      integer :: nctype ! type of the nth variable
      integer :: dim_id,len_dim,len_x,len_y,var_id,xdim_id,xdim_id_out  &
                ,ydim_id,ydim_id_out
      integer :: istat
      integer,dimension(1:2) :: dimids=(/0,0/)
      real,dimension(:)  ,allocatable :: var_1d_with_halo
      real,dimension(:,:),allocatable :: var_2d_with_halo
      real(kind=kdbl),dimension(:,:),allocatable :: var_2d_dbl_with_halo
!     real*8,dimension(:,:),allocatable :: var_2d_dbl_with_halo
      character(len=50) :: file,name_dim,name_xdim,name_ydim   &
                          ,xdim,ydim
      character(len=50) :: name_att ! the attribute's name
      character(len=50) :: name_var ! the name of the nth variable
      character(len=255) :: att=' '
      character(len=255),dimension(:),allocatable :: var_1d_char

!-----------------------------------------------------------------------
!***  Read in the required compute dimensions, halo size, and filenames.
!-----------------------------------------------------------------------
      read(5,*)idim_compute,jdim_compute,halo,filename_full,filename_shaved
      write(6,*)' id ',idim_compute,' jd ',jdim_compute,' halo ',halo
      write(6,*)' fn_f ',trim(filename_full)
      write(6,*)' fn_s ',trim(filename_shaved)
      i_count_compute=idim_compute+2*halo
      j_count_compute=jdim_compute+2*halo
      i_count_super  =2*i_count_compute
      j_count_super  =2*j_count_compute

!-----------------------------------------------------------------------
!***  Open the netcdf file with the incoming data to be shaved.
!-----------------------------------------------------------------------
      call check(nf90_open(filename_full,nf90_nowrite,ncid_in))            ! Open the netcdf file; get the file ID.
      call check(nf90_inquire(ncid_in,ndims,nvars,ngatts,unlimdimid))      ! Find the number of variables in the file.

!-----------------------------------------------------------------------
!***  Create the NetCDF file that the shaved data will be written into.
!***  Match the NetCDF format of the input file.
!-----------------------------------------------------------------------
      call check(nf90_create(filename_shaved                            &  !
                            ,or(nf90_classic_model,nf90_netcdf4)        &  ! Create NetCDF file for shaved data.
                            ,ncid_out))                                    !

!-----------------------------------------------------------------------
!***  Replicate the dimensions from the input to the output file
!***  but change values as needed to account for the shaving.  We
!***  know the grid file and orog file have given names for their
!***  x and y dimensions so use those to adjust to shaved values.
!***  NOTE:  Gridpoints in the grid file are on the domain's
!***         supergrid while points in the orog file are on the
!***         compute grid.
!-----------------------------------------------------------------------
      do nd=1,ndims
        call check(nf90_inquire_dimension(ncid_in,nd,name_dim,len_dim))    ! Get this dimension's name and value.
        select case (name_dim)
          case ('nx')                      !--- 
            len_dim=i_count_super               !
            xdim=name_dim                       !
            file='grid_file'                    !
          case ('nxp')                          ! Used by the
            len_dim=i_count_super+1             ! grid file.
          case ('ny')                           !
            len_dim=j_count_super               !
            ydim=name_dim                       !
          case ('nyp')                          !
            len_dim=j_count_super+1        !--- 
          case ('lon')                     !--- 
            len_dim=i_count_compute             !
            xdim=name_dim                       ! Used by the
            file='orog_file'                    ! orog file.
          case ('lat')                          !
            len_dim=j_count_compute             !
            ydim=name_dim                  !---
        end select
        dim_id=nd
        call check(nf90_def_dim(ncid_out,name_dim,len_dim,dim_id))         !-- Insert dimension into the output file.
      enddo

!-----------------------------------------------------------------------
!***  The output file's variables must be defined while that file
!***  is still in define mode.  Loop through the variables in the 
!***  input file and define each of them in the output file.
!-----------------------------------------------------------------------
      do n=1,nvars
        var_id=n
        call check(nf90_inquire_variable(ncid_in,var_id,name_var,nctype &  !-- name and type of this variable
                  ,ndims,dimids,natts))                                    !-- # of dimensions, ID, and attributes in this variable
        if(ndims==1)then
          call check(nf90_def_var(ncid_out,name_var,nctype,dimids(1),var_id)) !-- Define this 1-D variable in the output file.
        elseif(ndims==2)then
          call check(nf90_def_var(ncid_out,name_var,nctype,dimids,var_id))    !-- Define this 2-D variable in the output file.
        endif

!-----------------------------------------------------------------------
!***  Copy this variable's attributes to the output file's variable.
!-----------------------------------------------------------------------
        if(natts>0)then
          do na=1,natts
            call check(nf90_inq_attname(ncid_in,var_id,na,name_att))       !-- Get the attribute's name and ID from input file.
            call check(nf90_copy_att(ncid_in,var_id,name_att,ncid_out,var_id))  !-- Copy to output file.
          enddo
        endif
      enddo

!-----------------------------------------------------------------------
!***  Copy the global attributes to the output file.
!-----------------------------------------------------------------------
      do n=1,ngatts
        call check(nf90_inq_attname(ncid_in,nf90_global,n,name_att))
        call check(nf90_copy_att(ncid_in,nf90_global,name_att,ncid_out,nf90_global))
      enddo
      call check(nf90_enddef(ncid_out))                                    !-- Put the output file into data mode.

!-----------------------------------------------------------------------
!***  Get the x and y extents of the incoming grid with extra rows 
!***  so we can find determine how many rows to shave off. 
!-----------------------------------------------------------------------
      call check(nf90_inq_dimid(ncid_in,xdim,xdim_id))                     !-- Find the ID of the x dimension.
      call check(nf90_inq_dimid(ncid_in,ydim,ydim_id))                     !-- Find the ID of the y dimension.
      call check(nf90_inquire_dimension(ncid_in,xdim_id,name_xdim,len_x))  !-- Length of x dimension of vars in incoming file.
      call check(nf90_inquire_dimension(ncid_in,ydim_id,name_ydim,len_y))  !-- Length of y dimension of vars in incoming file.
      if(trim(file)=='orog_file')then
        i_start=(len_x-idim_compute)/2-halo+1                              !-- Starting i of 2-D data with halo rows on compute grid.
        j_start=(len_y-jdim_compute)/2-halo+1                              !-- Starting j of 2-D data with halo rows on compute grid.
 
      elseif(trim(file)=='grid_file')then
        i_start=(len_x-2*idim_compute)/2-2*halo+1                          !-- Starting i of 2-D data with halo rows on supergrid.
        j_start=(len_y-2*jdim_compute)/2-2*halo+1                          !-- Starting j of 2-D data with halo rows on supergrid.
      endif

!-----------------------------------------------------------------------
!***  We assume the # of extra rows on the incoming data is the same
!***  in both x and y so the # of rows to shave off to leave the
!***  halo rows is also the same in x and y.  So consider only the 
!***  values from the x dimension.
!-----------------------------------------------------------------------
      n_shave=i_start-1                                                    !-- # of rows to shave off full data to leave halo rows.

!-----------------------------------------------------------------------
!***  Now loop through all the variables in the input netcdf file,
!***  read in the data excluding all extra rows except for halo rows,
!***  and then write that out to the output file.
!-----------------------------------------------------------------------
      var_loop: do n=1,nvars
        var_id=n
        call check(nf90_inquire_variable(ncid_in,var_id,name_var,nctype &  !-- The name and type of the nth variable
                  ,ndims,dimids,natts))                                    !-- The dimensions, ID, and attributes in the nth variable
        call check(nf90_inquire_dimension(ncid_in,dimids(1),name_xdim,len_x)) !-- Get the length of the input 1st dimension.
        if(ndims==2)then
          call check(nf90_inquire_dimension(ncid_in,dimids(2),name_ydim,len_y)) !-- Get the length of the input y dimension.
        endif

!-------------------
!***  1-D variables
!-------------------
        if(ndims==1)then

!---------------
!***  Character
!---------------
          if(nctype==nf90_char)then
            n_start=1                                                      !-- Start reading incoming character data at this location.
            n_count=len_x                                                  !-- Character data is not gridded so not shaved.
            allocate(var_1d_char(1:n_count),stat=istat)
            call check(nf90_get_var(ncid_in,var_id,var_1d_char(:) &        !-- Fill the 1-D character variable.
                                   ,start=(/n_start/) &
                                   ,count=(/n_count/)))
            call check(nf90_put_var(ncid_out,var_id,var_1d_char))          !-- Write out the 1-D character variable.
            deallocate(var_1d_char)

!---------------
!***  Numerical
!---------------
          else
            n_start=n_shave+1                                              !-- Start reading incoming data at this location.
            n_count=len_dim-2*n_shave                                      !-- # of datapoints to fill in the shaved 1-D variable.
            allocate(var_1d_with_halo(1:n_count),stat=istat)
            call check(nf90_get_var(ncid_in,var_id,var_1d_with_halo(:) &   !-- Fill the shaved 1-D variable.
                                   ,start=(/n_start/) &
                                   ,count=(/n_count/)))
            call check(nf90_put_var(ncid_out,var_id,var_1d_with_halo))     !-- Write out the shaved 1-D variable.
            deallocate(var_1d_with_halo)
          endif

!-------------------
!***  2-D variables
!-------------------
        elseif(ndims==2)then
          if(trim(file)=='orog_file')then
            i_start=(len_x-idim_compute)/2-halo+1                          !-- Starting i of 2-D data with halo rows on compute grid.
            j_start=(len_y-jdim_compute)/2-halo+1                          !-- Starting j of 2-D data with halo rows on compute grid.
            i_count=i_count_compute                                        !-- i extent of 2-D data with halo rows on compute grid.
            j_count=j_count_compute                                        !-- j extent of 2-D data with halo rows on compute grid.
          elseif(trim(file)=='grid_file')then
            i_start=(len_x-2*idim_compute)/2-2*halo+1                      !-- Starting i of 2-D data with halo rows on supergrid.
            j_start=(len_y-2*jdim_compute)/2-2*halo+1                      !-- Starting j of 2-D data with halo rows on supergrid.
            i_count=i_count_super                                          !-- i extent of 2-D data with halo rows on supergrid.
            j_count=j_count_super                                          !-- j extent of 2-D data with halo rows on supergrid.
            if(trim(name_xdim)=='nxp')then
              i_count=i_count+1                                            !-- nxp is # of cell corners in x, not centers.
            endif
            if(trim(name_ydim)=='nyp')then
              j_count=j_count+1                                            !-- nyp is # of cell corners in y, not centers.
            endif
          endif
          if(nctype==nf90_float)then                                       !-- Single precision real variables
            allocate(var_2d_with_halo(i_count,j_count),stat=istat)
            call check(nf90_get_var(ncid_in,var_id,var_2d_with_halo(:,:) & !-- Fill array with compute data plus halo rows.
                                   ,start=(/i_start,j_start/) &
                                   ,count=(/i_count,j_count/)))
            call check(nf90_put_var(ncid_out,var_id,var_2d_with_halo))     !-- Write out the shaved 2-D single precision variable.
            deallocate(var_2d_with_halo)
          elseif(nctype==nf90_double)then                                  !-- Double precision real variables
            allocate(var_2d_dbl_with_halo(i_count,j_count),stat=istat)
            call check(nf90_get_var(ncid_in,var_id,var_2d_dbl_with_halo(:,:) &  !-- Fill array with compute data plus halo rows.
                                   ,start=(/i_start,j_start/) &
                                   ,count=(/i_count,j_count/)))
            call check(nf90_put_var(ncid_out,var_id,var_2d_dbl_with_halo)) !-- Write out the shaved 2-D double precision variable.
            deallocate(var_2d_dbl_with_halo)
          endif
        endif
      enddo var_loop
      call check(nf90_close(ncid_out))
      call check(nf90_close(ncid_in))

      contains

      subroutine check(status)
      integer,intent(in) :: status
      if(status /= nf90_noerr) then
        print *, trim(nf90_strerror(status))
        stop "Stopped"
      end if
      end subroutine check
      end program shave_nc