2025-09-05 09:42:21 +02:00

710 lines
24 KiB
Fortran
Executable File

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! this program is a short driver that tests the remappings using
! a simple analytic field. the results are written in netCDF
! format.
!
!-----------------------------------------------------------------------
!
! CVS:$Id: scrip_test_repeat.f,v 1.3 2000/04/19 21:56:26 pwjones Exp $
!
! Copyright (c) 1997, 1998 the Regents of the University of
! California.
!
! This software and ancillary information (herein called software)
! called SCRIP is made available under the terms described here.
! The software has been approved for release with associated
! LA-CC Number 98-45.
!
! Unless otherwise indicated, this software has been authored
! by an employee or employees of the University of California,
! operator of the Los Alamos National Laboratory under Contract
! No. W-7405-ENG-36 with the U.S. Department of Energy. The U.S.
! Government has rights to use, reproduce, and distribute this
! software. The public may copy and use this software without
! charge, provided that this Notice and any statement of authorship
! are reproduced on all copies. Neither the Government nor the
! University makes any warranty, express or implied, or assumes
! any liability or responsibility for the use of this software.
!
! If software is modified to produce derivative works, such modified
! software should be clearly marked, so as not to confuse it with
! the version available from Los Alamos National Laboratory.
!
!***********************************************************************
program remap_test_repeat
!-----------------------------------------------------------------------
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
implicit none
!-----------------------------------------------------------------------
!
! interface for remap routine
!
!-----------------------------------------------------------------------
interface
subroutine remap(dst_array, map_wts, dst_add, src_add,
& src_array, src_grad1, src_grad2, src_grad3)
use kinds_mod
use constants
implicit none
integer (kind=int_kind), dimension(:), intent(in) ::
& dst_add, src_add
real (kind=dbl_kind), dimension(:,:), intent(in) :: map_wts
real (kind=dbl_kind), dimension(:), intent(in) :: src_array
real (kind=dbl_kind), dimension(:), intent(in), optional ::
& src_grad1, src_grad2, src_grad3
real (kind=dbl_kind), dimension(:), intent(inout) :: dst_array
end subroutine remap
end interface
!-----------------------------------------------------------------------
!
! input namelist variables
!
!-----------------------------------------------------------------------
integer (kind=int_kind) ::
& field_choice, ! choice of field to be interpolated
& num_repeats ! number of times to repeat remappings
character (char_len) ::
& interp_file1, ! filename containing remap data (map1)
& interp_file2, ! filename containing remap data (map2)
& output_file ! filename containing output test data
namelist /remap_inputs/ field_choice, num_repeats,
& interp_file1, interp_file2, output_file
!-----------------------------------------------------------------------
!
! local variables
!
!-----------------------------------------------------------------------
character (char_len) ::
& map_name1, ! name for mapping from grid1 to grid2
& map_name2 ! name for mapping from grid2 to grid1
integer (kind=int_kind) :: ! netCDF ids for files and arrays
& n, ncstat, nc_outfile_id,
& nc_srcgrdcntrlat_id, nc_srcgrdcntrlon_id,
& nc_dstgrdcntrlat_id, nc_dstgrdcntrlon_id,
& nc_srcgrdrank_id, nc_dstgrdrank_id,
& nc_srcgrdimask_id, nc_dstgrdimask_id,
& nc_srcgrdarea_id, nc_dstgrdarea_id,
& nc_srcgrdfrac_id, nc_dstgrdfrac_id,
& nc_srcarray_id, nc_srcgradlat_id, nc_srcgradlon_id,
& nc_dstarray1_id, nc_dstarray2_id,
& nc_dsterror1_id, nc_dsterror2_id
integer (kind=int_kind), dimension(:), allocatable ::
& nc_grid1size_id, nc_grid2size_id
!-----------------------------------------------------------------------
character (char_len) ::
& dim_name ! netCDF dimension name
integer (kind=int_kind) :: i,j,n,
& iunit ! unit number for namelist file
integer (kind=int_kind), dimension(:), allocatable ::
& grid1_imask, grid2_imask
real (kind=dbl_kind) ::
& length ! length scale for cosine hill test field
real (kind=dbl_kind), dimension(:), allocatable ::
& grid1_array,
& grid1_tmp,
& grad1_lat,
& grad1_lon,
& grid2_array,
& grid2_err,
& grid2_tmp,
& grad2_lat,
& grad2_lon
!-----------------------------------------------------------------------
!
! read namelist for file and mapping info
!
!-----------------------------------------------------------------------
call get_unit(iunit)
open(iunit, file='repeat_test_in', status='old', form='formatted')
read(iunit, nml=remap_inputs)
call release_unit(iunit)
write(*,nml=remap_inputs)
!-----------------------------------------------------------------------
!
! read remapping data
!
!-----------------------------------------------------------------------
num_maps = 2
call read_remap(map_name1, map_name2, interp_file1, interp_file2)
!-----------------------------------------------------------------------
!
! allocate arrays
!
!-----------------------------------------------------------------------
allocate (grid1_array (grid1_size),
& grid1_tmp (grid1_size),
& grad1_lat (grid1_size),
& grad1_lon (grid1_size),
& grid1_imask (grid1_size),
& grid2_array (grid2_size),
& grid2_err (grid2_size),
& grid2_tmp (grid2_size),
& grad2_lat (grid2_size),
& grad2_lon (grid2_size),
& grid2_imask (grid2_size))
where (grid1_mask)
grid1_imask = 1
elsewhere
grid1_imask = 0
endwhere
where (grid2_mask)
grid2_imask = 1
elsewhere
grid2_imask = 0
endwhere
!-----------------------------------------------------------------------
!
! setup a NetCDF file for output
!
!-----------------------------------------------------------------------
!***
!*** create netCDF dataset
!***
ncstat = nf_create (output_file, NF_CLOBBER, nc_outfile_id)
call netcdf_error_handler(ncstat)
ncstat = nf_put_att_text (nc_outfile_id, NF_GLOBAL, 'title',
& len_trim(map_name1), map_name1)
call netcdf_error_handler(ncstat)
!***
!*** define grid size dimensions
!***
allocate( nc_grid1size_id(grid1_rank),
& nc_grid2size_id(grid2_rank))
do n=1,grid1_rank
write(dim_name,1000) 'grid1_dim',n
ncstat = nf_def_dim (nc_outfile_id, dim_name,
& grid1_dims(n), nc_grid1size_id(n))
call netcdf_error_handler(ncstat)
end do
do n=1,grid2_rank
write(dim_name,1000) 'grid2_dim',n
ncstat = nf_def_dim (nc_outfile_id, dim_name,
& grid2_dims(n), nc_grid2size_id(n))
call netcdf_error_handler(ncstat)
end do
1000 format(a9,i1)
!***
!*** define grid center latitude array
!***
ncstat = nf_def_var (nc_outfile_id, 'src_grid_center_lat',
& NF_DOUBLE, grid1_rank, nc_grid1size_id,
& nc_srcgrdcntrlat_id)
call netcdf_error_handler(ncstat)
ncstat = nf_put_att_text (nc_outfile_id, nc_srcgrdcntrlat_id,
& 'units', 7, 'radians')
call netcdf_error_handler(ncstat)
ncstat = nf_def_var (nc_outfile_id, 'dst_grid_center_lat',
& NF_DOUBLE, grid2_rank, nc_grid2size_id,
& nc_dstgrdcntrlat_id)
call netcdf_error_handler(ncstat)
ncstat = nf_put_att_text (nc_outfile_id, nc_dstgrdcntrlat_id,
& 'units', 7, 'radians')
call netcdf_error_handler(ncstat)
!***
!*** define grid center longitude array
!***
ncstat = nf_def_var (nc_outfile_id, 'src_grid_center_lon',
& NF_DOUBLE, grid1_rank, nc_grid1size_id,
& nc_srcgrdcntrlon_id)
call netcdf_error_handler(ncstat)
ncstat = nf_put_att_text (nc_outfile_id, nc_srcgrdcntrlon_id,
& 'units', 7, 'radians')
call netcdf_error_handler(ncstat)
ncstat = nf_def_var (nc_outfile_id, 'dst_grid_center_lon',
& NF_DOUBLE, grid2_rank, nc_grid2size_id,
& nc_dstgrdcntrlon_id)
call netcdf_error_handler(ncstat)
ncstat = nf_put_att_text (nc_outfile_id, nc_dstgrdcntrlon_id,
& 'units', 7, 'radians')
call netcdf_error_handler(ncstat)
!***
!*** define grid mask
!***
ncstat = nf_def_var (nc_outfile_id, 'src_grid_imask', NF_INT,
& grid1_rank, nc_grid1size_id, nc_srcgrdimask_id)
call netcdf_error_handler(ncstat)
ncstat = nf_put_att_text (nc_outfile_id, nc_srcgrdimask_id,
& 'units', 8, 'unitless')
call netcdf_error_handler(ncstat)
ncstat = nf_def_var (nc_outfile_id, 'dst_grid_imask', NF_INT,
& grid2_rank, nc_grid2size_id, nc_dstgrdimask_id)
call netcdf_error_handler(ncstat)
ncstat = nf_put_att_text (nc_outfile_id, nc_dstgrdimask_id,
& 'units', 8, 'unitless')
call netcdf_error_handler(ncstat)
!***
!*** define grid area arrays
!***
ncstat = nf_def_var (nc_outfile_id, 'src_grid_area',
& NF_DOUBLE, grid1_rank, nc_grid1size_id,
& nc_srcgrdarea_id)
call netcdf_error_handler(ncstat)
ncstat = nf_def_var (nc_outfile_id, 'dst_grid_area',
& NF_DOUBLE, grid2_rank, nc_grid2size_id,
& nc_dstgrdarea_id)
call netcdf_error_handler(ncstat)
!***
!*** define grid fraction arrays
!***
ncstat = nf_def_var (nc_outfile_id, 'src_grid_frac',
& NF_DOUBLE, grid1_rank, nc_grid1size_id,
& nc_srcgrdfrac_id)
call netcdf_error_handler(ncstat)
ncstat = nf_def_var (nc_outfile_id, 'dst_grid_frac',
& NF_DOUBLE, grid2_rank, nc_grid2size_id,
& nc_dstgrdfrac_id)
call netcdf_error_handler(ncstat)
!***
!*** define source array
!***
ncstat = nf_def_var (nc_outfile_id, 'src_array',
& NF_DOUBLE, grid1_rank, nc_grid1size_id,
& nc_srcarray_id)
call netcdf_error_handler(ncstat)
!***
!*** define gradient arrays
!***
ncstat = nf_def_var (nc_outfile_id, 'src_grad_lat',
& NF_DOUBLE, grid1_rank, nc_grid1size_id,
& nc_srcgradlat_id)
call netcdf_error_handler(ncstat)
ncstat = nf_def_var (nc_outfile_id, 'src_grad_lon',
& NF_DOUBLE, grid1_rank, nc_grid1size_id,
& nc_srcgradlon_id)
call netcdf_error_handler(ncstat)
!***
!*** define destination arrays
!***
ncstat = nf_def_var (nc_outfile_id, 'dst_array1',
& NF_DOUBLE, grid2_rank, nc_grid2size_id,
& nc_dstarray1_id)
call netcdf_error_handler(ncstat)
ncstat = nf_def_var (nc_outfile_id, 'dst_array2',
& NF_DOUBLE, grid2_rank, nc_grid2size_id,
& nc_dstarray2_id)
call netcdf_error_handler(ncstat)
!***
!*** define error arrays
!***
ncstat = nf_def_var (nc_outfile_id, 'dst_error1',
& NF_DOUBLE, grid2_rank, nc_grid2size_id,
& nc_dsterror1_id)
call netcdf_error_handler(ncstat)
ncstat = nf_def_var (nc_outfile_id, 'dst_error2',
& NF_DOUBLE, grid2_rank, nc_grid2size_id,
& nc_dsterror2_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_srcgrdcntrlat_id,
& grid1_center_lat)
call netcdf_error_handler(ncstat)
ncstat = nf_put_var_double(nc_outfile_id, nc_dstgrdcntrlat_id,
& grid2_center_lat)
call netcdf_error_handler(ncstat)
!***
!*** write grid center longitude array
!***
ncstat = nf_put_var_double(nc_outfile_id, nc_srcgrdcntrlon_id,
& grid1_center_lon)
call netcdf_error_handler(ncstat)
ncstat = nf_put_var_double(nc_outfile_id, nc_dstgrdcntrlon_id,
& grid2_center_lon)
call netcdf_error_handler(ncstat)
!***
!*** write grid mask
!***
ncstat = nf_put_var_int(nc_outfile_id, nc_srcgrdimask_id,
& grid1_imask)
call netcdf_error_handler(ncstat)
ncstat = nf_put_var_int(nc_outfile_id, nc_dstgrdimask_id,
& grid2_imask)
call netcdf_error_handler(ncstat)
!***
!*** define grid area arrays
!***
ncstat = nf_put_var_double(nc_outfile_id, nc_srcgrdarea_id,
& grid1_area)
call netcdf_error_handler(ncstat)
ncstat = nf_put_var_double(nc_outfile_id, nc_dstgrdarea_id,
& grid2_area)
call netcdf_error_handler(ncstat)
!***
!*** define grid fraction arrays
!***
ncstat = nf_put_var_double(nc_outfile_id, nc_srcgrdfrac_id,
& grid1_frac)
call netcdf_error_handler(ncstat)
ncstat = nf_put_var_double(nc_outfile_id, nc_dstgrdfrac_id,
& grid2_frac)
call netcdf_error_handler(ncstat)
!-----------------------------------------------------------------------
!
! set up fields for test cases based on user choice
!
!-----------------------------------------------------------------------
select case (field_choice)
case(1) !*** cosine hill at lon=pi and lat=0
length = 0.1*pi2
grid1_array = cos(grid1_center_lat)*cos(grid1_center_lon)
grid2_array = cos(grid2_center_lat)*cos(grid2_center_lon)
grid1_tmp = acos(-grid1_array)/length
grid2_tmp = acos(-grid2_array)/length
where (grid1_tmp <= one)
grad1_lat = (pi/length)*sin(pi*grid1_tmp)*
& sin(grid1_center_lat)*cos(grid1_center_lon)/
& sqrt(one-grid1_array**2)
grad1_lon = (pi/length)*sin(pi*grid1_tmp)*
& sin(grid1_center_lon)/
& sqrt(one-grid1_array**2)
grid1_array = two + cos(pi*grid1_tmp)
elsewhere
grid1_array = one
grad1_lat = zero
grad1_lon = zero
endwhere
where (grid2_tmp <= one)
grad2_lat = (pi/length)*sin(pi*grid2_tmp)*
& sin(grid2_center_lat)*cos(grid2_center_lon)/
& sqrt(one-grid2_array**2)
grad2_lon = (pi/length)*sin(pi*grid2_tmp)*
& sin(grid2_center_lon)/
& sqrt(one-grid2_array**2)
grid2_array = two + cos(pi*grid2_tmp)
elsewhere
grid2_array = one
grad2_lat = zero
grad2_lon = zero
endwhere
where (.not. grid1_mask)
grid1_array = zero
grad1_lat = zero
grad1_lon = zero
end where
where (.not. grid2_mask)
grid2_array = zero
grad2_lat = zero
grad2_lon = zero
end where
case(2) !*** pseudo-spherical harmonic l=2,m=2
where (grid1_mask)
grid1_array = two + cos(grid1_center_lat)**2*
& cos(two*grid1_center_lon)
grad1_lat = -sin(two*grid1_center_lat)*
& cos(two*grid1_center_lon)
grad1_lon = -two*cos(grid1_center_lat)*
& sin(two*grid1_center_lon)
elsewhere
grid1_array = zero
grad1_lat = zero
grad1_lon = zero
end where
where (grid2_mask)
grid2_array = two + cos(grid2_center_lat)**2*
& cos(two*grid2_center_lon)
grad2_lat = -sin(two*grid2_center_lat)*
& cos(two*grid2_center_lon)
grad2_lon = -two*cos(grid2_center_lat)*
& sin(two*grid2_center_lon)
elsewhere
grid2_array = zero
grad2_lat = zero
grad2_lon = zero
end where
case(3) !*** pseudo-spherical harmonic l=32, m=16
where (grid1_mask)
grid1_array = two + sin(two*grid1_center_lat)**16*
& cos(16.*grid1_center_lon)
grad1_lat = 32.*sin(two*grid1_center_lat)**15*
& cos(two*grid1_center_lat)*
& cos(16.*grid1_center_lon)
grad1_lon = -32.*sin(two*grid1_center_lat)**15*
& sin(grid1_center_lat)*
& sin(16.*grid1_center_lon)
elsewhere
grid1_array = zero
grad1_lat = zero
grad1_lon = zero
end where
where (grid2_mask)
grid2_array = two + sin(two*grid2_center_lat)**16*
& cos(16.*grid2_center_lon)
grad2_lat = 32.*sin(two*grid2_center_lat)**15*
& cos(two*grid2_center_lat)*
& cos(16.*grid2_center_lon)
grad2_lon = -32.*sin(two*grid2_center_lat)**15*
& sin(grid2_center_lat)*
& sin(16.*grid2_center_lon)
elsewhere
grid2_array = zero
grad2_lat = zero
grad2_lon = zero
end where
case default
stop 'Bad choice for field to interpolate'
end select
!-----------------------------------------------------------------------
!
! test repeated first-order maps between grid1 and grid2
!
!-----------------------------------------------------------------------
call remap(grid2_tmp, wts_map1, grid2_add_map1, grid1_add_map1,
& grid1_array)
do n=1,num_repeats
call remap(grid1_tmp, wts_map2, grid1_add_map2, grid2_add_map2,
& grid2_tmp)
call remap(grid2_tmp, wts_map1, grid2_add_map1, grid1_add_map1,
& grid1_tmp)
end do
where (grid2_frac > .999)
grid2_err = (grid2_tmp - grid2_array)/grid2_array
elsewhere
grid2_err = zero
end where
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_tmp ),maxval(grid2_tmp )
print *,' Err2 min,max: ',minval(grid2_err),maxval(grid2_err)
print *,' Err2 mean: ',sum(abs(grid2_err))/
& count(grid2_frac > .001)
!***
!*** Conservation Test
!***
print *,'Conservation:'
print *,'Grid1 Integral = ',sum(grid1_array*grid1_area*grid1_frac)
print *,'Grid2 Integral = ',sum(grid2_tmp *grid2_area*grid2_frac)
!-----------------------------------------------------------------------
!
! write results to NetCDF file
!
!-----------------------------------------------------------------------
ncstat = nf_put_var_double(nc_outfile_id, nc_srcarray_id,
& grid1_array)
call netcdf_error_handler(ncstat)
ncstat = nf_put_var_double(nc_outfile_id, nc_dstarray1_id,
& grid2_tmp )
call netcdf_error_handler(ncstat)
ncstat = nf_put_var_double(nc_outfile_id, nc_dsterror1_id,
& grid2_err)
call netcdf_error_handler(ncstat)
!-----------------------------------------------------------------------
!
! test repeated second-order mapppings between grid1 and grid2
!
!-----------------------------------------------------------------------
if (num_wts > 1) then
call remap(grid2_tmp , wts_map1, grid2_add_map1, grid1_add_map1,
& grid1_array, src_grad1=grad1_lat,
& src_grad2=grad1_lon)
do n=1,num_repeats
call remap(grid1_tmp, wts_map2, grid1_add_map2, grid2_add_map2,
& grid2_tmp, src_grad1=grad2_lat,
& src_grad2=grad2_lon)
call remap(grid2_tmp, wts_map1, grid2_add_map1, grid1_add_map1,
& grid1_tmp, src_grad1=grad1_lat,
& src_grad2=grad1_lon)
end do
where (grid2_frac > .999)
grid2_err = (grid2_tmp - grid2_array)/grid2_array
elsewhere
grid2_err = zero
end where
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_tmp ),maxval(grid2_tmp )
print *,' Err2 min,max: ',minval(grid2_err),maxval(grid2_err)
print *,' Err2 mean: ',sum(abs(grid2_err))/
& count(grid2_frac > .001)
!***
!*** Conservation Test
!***
print *,'Conservation:'
print *,'Grid1 Integral = ',sum(grid1_array*grid1_area*grid1_frac)
print *,'Grid2 Integral = ',sum(grid2_tmp *grid2_area*grid2_frac)
!-----------------------------------------------------------------------
!
! write results to NetCDF file
!
!-----------------------------------------------------------------------
ncstat = nf_put_var_double(nc_outfile_id, nc_srcgradlon_id,
& grad1_lon)
call netcdf_error_handler(ncstat)
ncstat = nf_put_var_double(nc_outfile_id, nc_dstarray2_id,
& grid2_tmp )
call netcdf_error_handler(ncstat)
ncstat = nf_put_var_double(nc_outfile_id, nc_dsterror2_id,
& grid2_err)
call netcdf_error_handler(ncstat)
endif
!-----------------------------------------------------------------------
!
! close netCDF file
!
!-----------------------------------------------------------------------
ncstat = nf_close(nc_outfile_id)
call netcdf_error_handler(ncstat)
!-----------------------------------------------------------------------
end program remap_test_repeat
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!