!----------------------------------------------------------------------- program shave_nc !----------------------------------------------------------------------- use netcdf !----------------------------------------------------------------------- implicit none !----------------------------------------------------------------------- ! !*** The grid driver step in FV3 preprocessing generates a grid_tile !*** file and an oro_tile file for the regional domain. The final !*** size of these files' domains must include the halo surrounding !*** the computational domain. However the original size of these !*** domains must exceed the domain size plus haloes so that the !*** topography filtering program can produce correct values over !*** the halo region. Then before the files go into the chgres !*** job their domains must be shaved down to only the computational !*** interior and the halo which is what this code does. ! !----------------------------------------------------------------------- ! integer,parameter :: kdbl=selected_real_kind(p=13,r=200) ! character(len=255) :: filename_full,filename_shaved integer :: idim_compute,jdim_compute,halo integer :: i_count_compute,j_count_compute & ,i_count_super,j_count_super integer :: i_start,j_start,i_count,j_count & ,n_count,n_shave,n_start integer :: n,na,natts,ncid_in,ncid_out,nctype,nd,ndims,ngatts & ,nvars,unlimdimid 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_att,name_dim,name_xdim,name_ydim & ,name_var,xdim,ydim 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 ! !-----------------------------------------------------------------------