Skip to content

Instantly share code, notes, and snippets.

@andersy005

andersy005/zonal_avg_mod.F90

Created December 20, 2019 19:49
Show Gist options
  • Save andersy005/44e063e56bed328e8651eb41633c8230 to your computer and use it in GitHub Desktop.
Save andersy005/44e063e56bed328e8651eb41633c8230 to your computer and use it in GitHub Desktop.
Download ZIP
compute zonal mean on a POP grid by @klindsay28
Raw
zonal_avg_mod.F90
module zonal_avg_mod
!-----------------------------------------------------------------------------
! This module computes basin zonal averages on POP grids. The steps are
! 1) set up the destination grid
! 2) compute averaging weights for each grid cell
! 3) compute normalizing weights for each basin (if required)
! 4) compute basin zonal averages
!
! The averages can be masked in longitude or latitude.
!
! CVS:$Id: zonal_avg_mod.F90,v 1.7 2002/02/11 05:04:40 klindsay Exp $
! CVS:$Name: $
! CVS:$Log: zonal_avg_mod.F90,v $
! CVS:Revision 1.7 2002/02/11 05:04:40 klindsay
! CVS:add options -v and -x
! CVS:
! CVS:Revision 1.6 2002/02/11 03:18:13 klindsay
! CVS:add -N option, which produces area-weighted sums, as opposed to averages
! CVS: together with running sums on latitude, is useful for
! CVS: making inferred transports from fluxes
! CVS:use NF_FILL_DOUBLE for missing_value and _FillValue if variable doesn't
! CVS: already have one
! CVS:
! CVS:Revision 1.5 2002/02/07 13:33:00 klindsay
! CVS:scale line integral areas for partial cells so that they sum to model
! CVS:generated area
! CVS:
! CVS:Revision 1.4 2002/02/07 00:08:05 klindsay
! CVS:fix bug (array length mismatch) in computation of dest_lat_axis
! CVS:also use TLAT for dest_lat_axis where grid has constant latitude
! CVS:continue to use averaging for TLAT where grid has variable latitude
! CVS:
! CVS:Revision 1.3 2001/11/21 22:56:44 klindsay
! CVS:Changed out_dimind to in_dimind, as it really specifies
! CVS:a dimension index in the input variable.
! CVS:
! CVS:Revision 1.2 2001/11/21 20:08:31 klindsay
! CVS:add additional valid depth axis names
! CVS:
! CVS:Revision 1.1 2001/10/15 17:37:54 klindsay
! CVS:initial checkin
! CVS:
!-----------------------------------------------------------------------------
use kinds_mod
use constants
use nf_wrap
use POP_grid_mod
use sphere_area_mod
implicit none
private
public :: &
gen_dest_lat_axis, &
def_dest_dimsnvars, &
put_dest_dim_vars, &
gen_weights, &
compute_zon_avg, &
varlist_include, &
parse_varlist, &
free_zon_avg
integer(kind=int_kind), save :: &
dest_lat_len, & ! length of generated latitude axis
init_unlimdim_len ! initial length of unlimited dimension in outfile
real(kind=r8), dimension(:), allocatable, save :: &
dest_lat_axis, & ! cell midpoints of destination axis
dest_lat_axis_edges ! cell edges of destination axis
integer(kind=int_kind), dimension(:,:), allocatable, save :: &
wght_cnt, & ! no. of weights for (i,j) T cell
lat_axis_ind0, & ! first lat_axis cell (i,j) T cell intersects with
Twght_ind0 ! first index in Twght array
real(kind=r8), dimension(:), allocatable, save :: &
Twght ! area common to tracer cells and latitude bands
real(kind=r8), dimension(:,:), allocatable, save :: &
var_in_lev, & ! data for 1 level on source grid
var_out_lev ! basin zonal averages for 1 level
real(kind=r8), dimension(:,:,:), allocatable, save :: &
Twght_basin_sum_kmt ! sum of Twght across each basin for kmt value
logical(kind=log_kind), save :: &
varlist_include = .true. ! is varlist an inclusion list or exclusion list
character(len=long_char_len), save :: &
varlist_pad = 'empty' ! comma padded version of varlist
contains
!-----------------------------------------------------------------------------
subroutine gen_dest_lat_axis(grid_var, latmin, latmax)
!--------------------------------------------------------------------------
! Generate destination axis for zonal averaging. For the genral case,
! 1) generate axis_edges ignoring clipping
! 2) generate axis
! 3) clip where needed
!--------------------------------------------------------------------------
!--------------------------------------------------------------------------
! arguments
!--------------------------------------------------------------------------
type(POP_grid), intent(in) :: &
grid_var ! grid to generate destination axis from
real(kind=r8), intent(in) :: &
latmin, latmax ! latitude bounds of averaging
!--------------------------------------------------------------------------
! local variables
!--------------------------------------------------------------------------
real(kind=r8), dimension(:), allocatable :: &
tmp_dest_lat_axis_edges ! cell edges of destination axis before clipping
real(kind=r8) :: &
dlat ! change in lat across first T cell
integer(kind=int_kind) :: &
j, & ! loop index
max_j_single, & ! maximum j value for which latitude on U grid is constant
tmp_len, & ! length of tmp_dest_lat_axis_edges - 1
jlo, jhi ! lo & hi j indices of tmp_dest_lat_axis_edges that
! get copied to dest_lat_axis_edges
!--------------------------------------------------------------------------
! compute max_j_reg
!--------------------------------------------------------------------------
max_j_single = 0
do j = 1,grid_var%nlat
if (maxval(grid_var%ulat(:,j)) - minval(grid_var%ulat(:,j)) < 1.0e-10_r8) then
max_j_single = j
else
exit
end if
end do
!--------------------------------------------------------------------------
! generate destination axis ignoring clipping
! first generate edges, then generate cell midpoints
!--------------------------------------------------------------------------
if (max_j_single == grid_var%nlat) then
!-----------------------------------------------------------------------
! case 1 : all constant j lines have constant latitude
!-----------------------------------------------------------------------
tmp_len = grid_var%nlat
allocate(tmp_dest_lat_axis_edges(0:tmp_len))
tmp_dest_lat_axis_edges(0) = grid_var%tvert_latmin
tmp_dest_lat_axis_edges(1:tmp_len) = grid_var%ulat(1,:)
else if (max_j_single == 0) then
!-----------------------------------------------------------------------
! case 2 : no constant j lines have constant latitude
!-----------------------------------------------------------------------
tmp_len = grid_var%nlat
allocate(tmp_dest_lat_axis_edges(0:tmp_len))
do j = 0,tmp_len
tmp_dest_lat_axis_edges(j) = grid_var%tvert_latmin + &
(grid_var%tvert_latmax - grid_var%tvert_latmin) * &
real(j,r8) / real(tmp_len,r8)
end do
else if (grid_var%ulat(1,max_j_single) == c0) then
!-----------------------------------------------------------------------
! case 3 : constant j lines have constant latitude go to equator
!-----------------------------------------------------------------------
tmp_len = max_j_single
if (-grid_var%ulat(1,1) > grid_var%tvert_latmax) then
do j = max_j_single-1,1,-1
tmp_len = tmp_len + 1
if (-grid_var%ulat(1,j) > grid_var%tvert_latmax) exit
end do
allocate(tmp_dest_lat_axis_edges(0:tmp_len))
tmp_dest_lat_axis_edges(0) = grid_var%tvert_latmin
tmp_dest_lat_axis_edges(1:max_j_single) = &
grid_var%ulat(1,1:max_j_single)
do j = max_j_single+1,tmp_len
tmp_dest_lat_axis_edges(j) = &
-tmp_dest_lat_axis_edges(max_j_single-(j-max_j_single))
end do
else
tmp_len = tmp_len + (max_j_single-1)
dlat = grid_var%ulat(1,1) - grid_var%tvert_latmin
tmp_len = tmp_len + &
ceiling((grid_var%tvert_latmax + grid_var%ulat(1,1)) / dlat)
allocate(tmp_dest_lat_axis_edges(0:tmp_len))
tmp_dest_lat_axis_edges(0) = grid_var%tvert_latmin
tmp_dest_lat_axis_edges(1:max_j_single) = &
grid_var%ulat(1,1:max_j_single)
do j = max_j_single+1,2*max_j_single-1
tmp_dest_lat_axis_edges(j) = &
-tmp_dest_lat_axis_edges(max_j_single-(j-max_j_single))
end do
do j = 2*max_j_single,tmp_len
tmp_dest_lat_axis_edges(j) = tmp_dest_lat_axis_edges(j-1) + dlat
end do
end if
else
!-----------------------------------------------------------------------
! case 4 : constant j lines do not stop at equator
!-----------------------------------------------------------------------
tmp_len = grid_var%nlat
allocate(tmp_dest_lat_axis_edges(0:tmp_len))
tmp_dest_lat_axis_edges(0) = grid_var%tvert_latmin
tmp_dest_lat_axis_edges(1:max_j_single) = &
grid_var%ulat(1,1:max_j_single)
do j = max_j_single+1,tmp_len
tmp_dest_lat_axis_edges(j) = grid_var%ulat(1,max_j_single) + &
(grid_var%tvert_latmax - grid_var%ulat(1,max_j_single)) * &
real(j-max_j_single,r8) / real(tmp_len-max_j_single,r8)
end do
end if
!--------------------------------------------------------------------------
! generate actual returned axes, taking clipping into account
!--------------------------------------------------------------------------
if (tmp_dest_lat_axis_edges(0) < latmin) then
do j = 1,tmp_len
if (tmp_dest_lat_axis_edges(j) > latmin) then
jlo = j-1
exit
end if
end do
else
jlo = 0
end if
if (tmp_dest_lat_axis_edges(tmp_len) > latmax) then
do j = tmp_len-1,0,-1
if (tmp_dest_lat_axis_edges(j) < latmax) then
jhi = j+1
exit
end if
end do
else
jhi = tmp_len
end if
dest_lat_len = jhi - jlo
allocate(dest_lat_axis(dest_lat_len), dest_lat_axis_edges(0:dest_lat_len))
dest_lat_axis_edges = tmp_dest_lat_axis_edges(jlo:jhi)
!--------------------------------------------------------------------------
! compute midpoints of cells for tracer points
!--------------------------------------------------------------------------
dest_lat_axis = dest_lat_axis_edges(0:dest_lat_len-1) + p5 * &
(dest_lat_axis_edges(1:dest_lat_len) - dest_lat_axis_edges(0:dest_lat_len-1))
!--------------------------------------------------------------------------
! in between where the source grid edges has constant latitude edges,
! use the source grid tracer point latitudes
!
! check first generated point for reasonableness, as some versions of POP
! produced unreasonable TLAT values for j=1
!--------------------------------------------------------------------------
if (max_j_single > jlo) then
dest_lat_axis(1:max_j_single-jlo) = grid_var%tlat(1,jlo+1:max_j_single)
if (dest_lat_axis(1) > dest_lat_axis_edges(1) .or. &
dest_lat_axis(1) < dest_lat_axis_edges(0)) then
dest_lat_axis(1) = dest_lat_axis_edges(0) + p5 * &
(dest_lat_axis_edges(1) - dest_lat_axis_edges(0))
end if
end if
!--------------------------------------------------------------------------
! clip latitude axes
!--------------------------------------------------------------------------
if (dest_lat_axis(1) < latmin) &
dest_lat_axis(1) = latmin
if (dest_lat_axis(dest_lat_len) > latmax) &
dest_lat_axis(dest_lat_len) = latmax
if (dest_lat_axis_edges(0) < latmin) &
dest_lat_axis_edges(0) = latmin
if (dest_lat_axis_edges(dest_lat_len) > latmax) &
dest_lat_axis_edges(dest_lat_len) = latmax
!--------------------------------------------------------------------------
! free tmp_dest_lat_axis_edges
!--------------------------------------------------------------------------
deallocate(tmp_dest_lat_axis_edges)
end subroutine gen_dest_lat_axis
!-----------------------------------------------------------------------------
subroutine def_dest_dimsnvars(outfilename, infilename, basin_cnt, time_const, lopen)
!--------------------------------------------------------------------------
! Define in outfile the variables from infile that are to be averaged.
! Ensure that the required dimensions and their associated variables are
! defined as well. Also copy global attributes from infile to outfile.
!--------------------------------------------------------------------------
!--------------------------------------------------------------------------
! arguments
!--------------------------------------------------------------------------
character(len=*), intent(in) :: &
outfilename, & ! name of file where dims & vars are to be defined
infilename ! name of file where dims & vars come from
integer(kind=int_kind), intent(in) :: &
basin_cnt ! number of basins in rmask, including global
logical(kind=log_kind), intent(in) :: &
time_const, & ! are time constant variables to be averaged?
lopen ! should file be opened vs. created
!--------------------------------------------------------------------------
! local variables
!--------------------------------------------------------------------------
character(len=*), parameter :: &
subname = 'zonal_avg_mod:def_dest_dims'
integer(kind=int_kind) :: &
i, & ! loop index
nlon_dimid, & ! dimid on "nlon" in infile
nlat_dimid, & ! dimid on "nlat" in infile
stat, & ! status from NetCDF calls
in_ncid, & ! file ID for infile
out_ncid, & ! file ID for outfile
dimlen, & ! dimension length
basin_dimid, & ! basin dimension ID in outfile
lat_dimid, & ! latitude dimension ID in outfile
in_nvars, & ! number of variables in infile
in_unlimdim, & ! unlimited dimension from infile
in_varid, & ! variable ID from infile
in_ndims, & ! ndims of in_varid
out_dimid, & ! dimid in outfile
out_varid, & ! varid in outfile
vartype, & ! NetCDF type of variable
dimvarid, & ! varid of a coordinate variable
dimvarndims, & ! number of dimensions of a coordinate variable
dimvardimids(1), & ! dimids of a coordinate variable
dimvartype, & ! NetCDF type of a coordinate variable
natts, & ! number of attributes a variable has
attnum, & ! number of a specific attribute
swap_dimid ! temporary dimid used for swapping
integer(kind=int_kind), dimension(NF_MAX_VAR_DIMS) :: &
in_dimids, & ! dimids of in_varid
out_dimids ! dimids of out_varid
character(len=char_len) :: &
att_name, & ! name of an attribute
att_text_val, & ! attribute value
varname, & ! variable name
dimname ! dimension name
real(kind=r8) :: &
msv ! value to be used where there is no data
!--------------------------------------------------------------------------
call nf_open_wrap(infilename, 0, in_ncid, subname)
if (lopen) then
call nf_open_wrap(outfilename, NF_WRITE, out_ncid, subname)
call nf_redef_wrap(out_ncid, subname)
else
call nf_create_wrap(outfilename, 0, out_ncid, subname)
end if
!--------------------------------------------------------------------------
! Define basin & latitude axes in outfile.
!--------------------------------------------------------------------------
call nf_inq_dimid_wrap(out_ncid, 'basins', basin_dimid, subname, &
ALLOW=NF_EBADDIM, stat_out=stat)
if (stat == NF_EBADDIM) then
call nf_def_dim_wrap(out_ncid, 'basins', basin_cnt, basin_dimid, subname)
end if
call nf_inq_dimid_wrap(out_ncid, 'lat_t', lat_dimid, subname, &
ALLOW=NF_EBADDIM, stat_out=stat)
if (stat == NF_EBADDIM) then
dimlen = dest_lat_len
call nf_def_dim_wrap(out_ncid, 'lat_t', dimlen, lat_dimid, subname)
end if
call nf_inq_varid_wrap(out_ncid, 'lat_t', out_varid, subname, &
ALLOW=NF_ENOTVAR, stat_out=stat)
if (stat == NF_ENOTVAR) then
call nf_def_var_wrap(out_ncid, 'lat_t', NF_FLOAT, 1, (/ lat_dimid /), &
out_varid, subname)
att_text_val = 'Latitude Cell Centers'
call nf_put_att_wrap(out_ncid, out_varid, 'long_name', &
len_trim(att_text_val), trim(att_text_val), subname)
att_text_val = 'degrees_north'
call nf_put_att_wrap(out_ncid, out_varid, 'units', &
len_trim(att_text_val), trim(att_text_val), subname)
att_text_val = 'lat_t_edges'
call nf_put_att_wrap(out_ncid, out_varid, 'edges', &
len_trim(att_text_val), trim(att_text_val), subname)
end if
call nf_inq_dimid_wrap(out_ncid, 'lat_t_edges', out_dimid, subname, &
ALLOW=NF_EBADDIM, stat_out=stat)
if (stat == NF_EBADDIM) then
dimlen = dest_lat_len+1
call nf_def_dim_wrap(out_ncid, 'lat_t_edges', dimlen, out_dimid, subname)
end if
call nf_inq_varid_wrap(out_ncid, 'lat_t_edges', out_varid, subname, &
ALLOW=NF_ENOTVAR, stat_out=stat)
if (stat == NF_ENOTVAR) then
call nf_def_var_wrap(out_ncid, 'lat_t_edges', NF_FLOAT, 1, &
(/ out_dimid /), out_varid, subname)
end if
!--------------------------------------------------------------------------
! Loop over variables in infile.
! 1) If variable is not to be averaged, go to next variable.
! 2) If variable is already defined, go to next variable.
! This does check that previously defined variable has same dimensions.
! 3) For each variable dimension not defined in outfile, except nlon & nlat
! 3a) define the dimension in outfile
! 3b) if it has a coordinate variable, define it, copying attributes
! 4) Define variable, changing (nlon,nlat) to (basins,lat_t).
! 5) Copy non-coordinate attributes from infile to outfile.
!--------------------------------------------------------------------------
!--------------------------------------------------------------------------
! obtain dimids of unlimdim, "nlon", and "nlat" from infile
! they are treated specially
!--------------------------------------------------------------------------
call nf_inq_unlimdim_wrap(in_ncid, in_unlimdim, subname)
call nf_inq_varid_wrap(in_ncid, TLONG_name, in_varid, subname)
call nf_inq_vardimid_wrap(in_ncid, in_varid, in_dimids, subname)
nlon_dimid = in_dimids(1)
nlat_dimid = in_dimids(2)
call nf_inq_nvars_wrap(in_ncid, in_nvars, subname)
do in_varid = 1,in_nvars
if (.not. to_comp_avg(in_ncid, in_varid, time_const)) cycle
!-----------------------------------------------------------------------
! Is the variable already defined in outfile?
!-----------------------------------------------------------------------
varname = ''
call nf_inq_varname_wrap(in_ncid, in_varid, varname, subname)
call nf_inq_varid_wrap(out_ncid, trim(varname), out_varid, subname, &
ALLOW=NF_ENOTVAR, stat_out=stat)
if (stat == NF_NOERR) cycle
call nf_inq_vartype_wrap(in_ncid, in_varid, vartype, subname)
call nf_inq_varndims_wrap(in_ncid, in_varid, in_ndims, subname)
call nf_inq_vardimid_wrap(in_ncid, in_varid, in_dimids, subname)
!-----------------------------------------------------------------------
! Construct dimids array
! Initial construction has nlon -> lat, nlat -> basin
! basin is shifted to end of dimid array later
!-----------------------------------------------------------------------
do i = 1,in_ndims
if (in_dimids(i) == nlon_dimid) then
out_dimids(i) = lat_dimid
cycle
end if
if (in_dimids(i) == nlat_dimid) then
out_dimids(i) = basin_dimid
cycle
end if
!--------------------------------------------------------------------
! Is the dimension already defined in outfile?
!--------------------------------------------------------------------
dimname = ''
call nf_inq_dimname_wrap(in_ncid, in_dimids(i), dimname, subname)
call nf_inq_dimid_wrap(out_ncid, trim(dimname), out_dimids(i), subname, &
ALLOW=NF_EBADDIM, stat_out=stat)
if (stat == NF_NOERR) cycle
!--------------------------------------------------------------------
! What is the length of the dimension? Treat the unlimited
! dimension correctly.
!--------------------------------------------------------------------
if (in_dimids(i) == in_unlimdim) then
dimlen = NF_UNLIMITED
else
call nf_inq_dimlen_wrap(in_ncid, in_dimids(i), dimlen, subname)
end if
!--------------------------------------------------------------------
! Define dimension in outfile
!--------------------------------------------------------------------
call nf_def_dim_wrap(out_ncid, trim(dimname), dimlen, out_dimids(i), subname)
!--------------------------------------------------------------------
! Does this dimension have a corresponding coordinate variable?
! I.e. a variable w/ the same name, defined on a single dimension
! that is the dimension in question
!--------------------------------------------------------------------
call nf_inq_varid_wrap(in_ncid, trim(dimname), dimvarid, subname, &
ALLOW=NF_ENOTVAR, stat_out=stat)
if (stat == NF_ENOTVAR) cycle
call nf_inq_varndims_wrap(in_ncid, dimvarid, dimvarndims, subname)
if (dimvarndims /= 1) cycle
call nf_inq_vardimid_wrap(in_ncid, dimvarid, dimvardimids, subname)
if (dimvardimids(1) /= in_dimids(i)) cycle
!--------------------------------------------------------------------
! Define coordinate variable in outfile
!--------------------------------------------------------------------
call nf_inq_vartype_wrap(in_ncid, dimvarid, dimvartype, subname)
call nf_def_var_wrap(out_ncid, trim(dimname), dimvartype, 1, &
(/ out_dimids(i) /), out_varid, subname)
!--------------------------------------------------------------------
! Copy non-coordinates attributes
!--------------------------------------------------------------------
call nf_inq_varnatts_wrap(in_ncid, dimvarid, natts, subname)
do attnum = 1,natts
att_name = ''
call nf_inq_attname_wrap(in_ncid, dimvarid, attnum, att_name, subname)
call nf_copy_att_wrap(in_ncid, dimvarid, trim(att_name), &
out_ncid, out_varid, subname)
end do
end do
!-----------------------------------------------------------------------
! Shift basin_dimid to end of list, but not beyond unlimited dimension
!-----------------------------------------------------------------------
do i = 1,in_ndims-1
if (out_dimids(i) == basin_dimid .and. in_dimids(i+1) /= in_unlimdim) then
swap_dimid = out_dimids(i+1)
out_dimids(i+1) = out_dimids(i)
out_dimids(i) = swap_dimid
end if
end do
!-----------------------------------------------------------------------
! Define variable in outfile
!-----------------------------------------------------------------------
call nf_def_var_wrap(out_ncid, trim(varname), vartype, in_ndims, &
out_dimids, out_varid, subname)
!-----------------------------------------------------------------------
! Copy non-coordinates attributes
!-----------------------------------------------------------------------
call nf_inq_varnatts_wrap(in_ncid, in_varid, natts, subname)
do attnum = 1,natts
att_name = ''
call nf_inq_attname_wrap(in_ncid, in_varid, attnum, att_name, subname)
if (att_name == 'coordinates') cycle
call nf_copy_att_wrap(in_ncid, in_varid, trim(att_name), &
out_ncid, out_varid, subname)
end do
!-----------------------------------------------------------------------
! If variable does not have _FillValue or missing_value attribute,
! put default value in for both.
!-----------------------------------------------------------------------
call nf_get_att_wrap(out_ncid, out_varid, '_FillValue', msv, subname, &
ALLOW=NF_ENOTATT, stat_out=stat)
if (stat == NF_ENOTATT) then
call nf_get_att_wrap(out_ncid, out_varid, 'missing_value', msv, subname, &
ALLOW=NF_ENOTATT, stat_out=stat)
if (stat == NF_ENOTATT) then
call nf_put_att_wrap(out_ncid, out_varid, '_FillValue', &
vartype, 1, NF_FILL_DOUBLE, subname)
call nf_put_att_wrap(out_ncid, out_varid, 'missing_value', &
vartype, 1, NF_FILL_DOUBLE, subname)
end if
end if
end do
!--------------------------------------------------------------------------
! Copy global attributes from infile to outfile
!--------------------------------------------------------------------------
call nf_inq_varnatts_wrap(in_ncid, NF_GLOBAL, natts, subname)
do attnum = 1,natts
att_name = ''
call nf_inq_attname_wrap(in_ncid, NF_GLOBAL, attnum, att_name, subname)
call nf_copy_att_wrap(in_ncid, NF_GLOBAL, trim(att_name), &
out_ncid, NF_GLOBAL, subname)
end do
call nf_close_wrap(out_ncid, subname)
call nf_close_wrap(in_ncid, subname)
end subroutine def_dest_dimsnvars
!-----------------------------------------------------------------------------
subroutine parse_varlist(varlist)
!--------------------------------------------------------------------------
! arguments
!--------------------------------------------------------------------------
character(len=long_char_len), intent(in) :: varlist
!--------------------------------------------------------------------------
varlist_pad = ',' // trim(varlist) // ','
end subroutine parse_varlist
!-----------------------------------------------------------------------------
logical(kind=log_kind) function to_comp_avg(ncid, varid, time_const)
!--------------------------------------------------------------------------
! arguments
!--------------------------------------------------------------------------
integer(kind=int_kind) :: &
ncid, & ! file ID
varid ! variable ID
logical(kind=log_kind) :: &
time_const ! are time constant variables to be averaged?
!--------------------------------------------------------------------------
! local variables
!--------------------------------------------------------------------------
character(len=*), parameter :: &
subname = 'zonal_avg_mod:to_comp_avg'
character(len=char_len) :: &
varname, & ! name of variable
att_text_val ! attribute value
integer(kind=int_kind) :: &
stat, & ! status from NetCDF calls
vartype, & ! variable type
timedim, & ! time dimension in ncid
ndims ! ndims of in_varid
integer(kind=int_kind), dimension(NF_MAX_VAR_DIMS) :: &
dimids ! dimids of varid
logical(kind=log_kind) :: &
var_in_list ! is variable in varlist
!--------------------------------------------------------------------------
to_comp_avg = .false.
call nf_inq_vartype_wrap(ncid, varid, vartype, subname)
if (vartype /= NF_DOUBLE .and. vartype /= NF_FLOAT) return
att_text_val = ''
call nf_get_att_wrap(ncid, varid, 'coordinates', att_text_val, &
ALLOW=NF_ENOTATT, stat_out=stat)
if (stat == NF_ENOTATT) return
if (index(att_text_val, TLONG_name) == 0) return
if (index(att_text_val, TLAT_name) == 0) return
if (varlist_pad == 'empty') then
if (.not. time_const) then
call nf_inq_dimid_wrap(ncid, 'time', timedim, subname, &
ALLOW=NF_EBADDIM, stat_out=stat)
if (stat /= NF_EBADDIM) then
call nf_inq_varndims_wrap(ncid, varid, ndims, subname)
call nf_inq_vardimid_wrap(ncid, varid, dimids, subname)
to_comp_avg = any(dimids(1:ndims) == timedim)
else
to_comp_avg = .false.
end if
else
to_comp_avg = .true.
end if
else
call nf_inq_varname_wrap(ncid, varid, varname, subname)
var_in_list = index(varlist_pad, ',' // trim(varname) // ',') > 0
to_comp_avg = var_in_list .eqv. varlist_include
end if
end function to_comp_avg
!-----------------------------------------------------------------------------
subroutine put_dest_dim_vars(outfilename, infilename)
!--------------------------------------------------------------------------
! put latitude axes and coordinate variables from infile into outfile
!--------------------------------------------------------------------------
!--------------------------------------------------------------------------
! arguments
!--------------------------------------------------------------------------
character(len=*), intent(in) :: &
outfilename, & ! file where dims are to be put
infilename ! file where dims come from
!--------------------------------------------------------------------------
! local variables
!--------------------------------------------------------------------------
character(len=*), parameter :: &
subname = 'zonal_avg_mod:put_dest_dim_vars'
integer(kind=int_kind) :: &
stat, & ! status from NetCDF calls
in_ncid, & ! file ID for infile
out_ncid, & ! file ID for outfile
out_varid, & ! varid of coordinate variable in outfile
out_unlimdim, & ! unlimited dimension from infile
out_ndims, & ! number of dimensions in outfile
out_dimid, & ! dimid in outfile
out_varndims, & ! number of dimensions of coordinate variable in outfile
out_vardimids(1),& ! dimids of coordinate variable in outfile
in_dimid, & ! dimid in infile
in_varid, & ! varid of coordinate variable in infile
in_varndims, & ! number of dimensions of coordinate variable in infile
in_vardimids(1), & ! dimids of coordinate variable in infile
vartype, & ! NetCDF type of coordinate variable
dimlen, & ! length of dimension
out_start(1), & ! start index of coordinate variable in outfile
out_count(1) ! length of coordinate variable in outfile
real(kind=r8), dimension(:), allocatable :: &
double_space ! space for double coordinate variable
real(kind=r4), dimension(:), allocatable :: &
float_space ! space for float coordinate variable
integer(kind=int_kind), dimension(:), allocatable :: &
int_space ! space for integer coordinate variable
character(len=char_len) :: &
dimname ! dimension name
!--------------------------------------------------------------------------
call nf_open_wrap(infilename, 0, in_ncid, subname)
call nf_open_wrap(outfilename, NF_WRITE, out_ncid, subname)
call nf_inq_varid_wrap(out_ncid, 'lat_t', out_varid, subname)
call nf_put_var_wrap(out_ncid, out_varid, dest_lat_axis, subname)
call nf_inq_varid_wrap(out_ncid, 'lat_t_edges', out_varid, subname)
call nf_put_var_wrap(out_ncid, out_varid, dest_lat_axis_edges, subname)
!--------------------------------------------------------------------------
! Loop over dimensions in outfile. For each dimension that has a
! coordinate variable, check to see if there is a corresponding
! dimension with coordinate variable in infile. If there is, copy
! the coordinate variable from infile to outfile.
!--------------------------------------------------------------------------
call nf_inq_unlimdim_wrap(out_ncid, out_unlimdim, subname)
call nf_inq_ndims_wrap(out_ncid, out_ndims, subname)
do out_dimid = 1,out_ndims
dimname = ''
call nf_inq_dimname_wrap(out_ncid, out_dimid, dimname, subname)
!-----------------------------------------------------------------------
! Does dimname have a corresponding coordinate variable in outfile?
!-----------------------------------------------------------------------
call nf_inq_varid_wrap(out_ncid, trim(dimname), out_varid, subname, &
ALLOW=NF_ENOTVAR, stat_out=stat)
if (stat == NF_ENOTVAR) cycle
call nf_inq_varndims_wrap(out_ncid, out_varid, out_varndims, subname)
if (out_varndims /= 1) cycle
call nf_inq_vardimid_wrap(out_ncid, out_varid, out_vardimids, subname)
if (out_vardimids(1) /= out_dimid) cycle
!-----------------------------------------------------------------------
! Does dimname occur in infile?
!-----------------------------------------------------------------------
call nf_inq_dimid_wrap(in_ncid, trim(dimname), in_dimid, subname, &
ALLOW=NF_EBADDIM, stat_out=stat)
if (stat == NF_EBADDIM) cycle
!-----------------------------------------------------------------------
! Does dimname have a corresponding coordinate variable in infile?
!-----------------------------------------------------------------------
call nf_inq_varid_wrap(in_ncid, trim(dimname), in_varid, subname, &
ALLOW=NF_ENOTVAR, stat_out=stat)
if (stat == NF_ENOTVAR) cycle
call nf_inq_varndims_wrap(in_ncid, in_varid, in_varndims, subname)
if (in_varndims /= 1) cycle
call nf_inq_vardimid_wrap(in_ncid, in_varid, in_vardimids, subname)
if (in_vardimids(1) /= in_dimid) cycle
!-----------------------------------------------------------------------
! allocate space for reading in variable
! read in variable
! write out variable
! deallocate space
!
! Use put_vara when writing out so that same call sequence can be
! used for unlimited dimension as other dimensions.
!-----------------------------------------------------------------------
call nf_inq_dimlen_wrap(in_ncid, in_dimid, dimlen, subname)
if (out_dimid /= out_unlimdim) then
out_start(1) = 1
else
call nf_inq_dimlen_wrap(out_ncid, out_dimid, init_unlimdim_len, subname)
out_start(1) = init_unlimdim_len + 1
end if
out_count(1) = dimlen
call nf_inq_vartype_wrap(in_ncid, in_varid, vartype, subname)
select case (vartype)
case (NF_DOUBLE)
allocate(double_space(dimlen))
call nf_get_var_wrap(in_ncid, in_varid, double_space, subname)
call nf_put_vara_wrap(out_ncid, out_varid, out_start, out_count, &
double_space, subname)
deallocate(double_space)
case (NF_FLOAT)
allocate(float_space(dimlen))
call nf_get_var_wrap(in_ncid, in_varid, float_space, subname)
call nf_put_vara_wrap(out_ncid, out_varid, out_start, out_count, &
float_space, subname)
deallocate(float_space)
case (NF_INT, NF_SHORT)
allocate(int_space(dimlen))
call nf_get_var_wrap(in_ncid, in_varid, int_space, subname)
call nf_put_vara_wrap(out_ncid, out_varid, out_start, out_count, &
int_space, subname)
deallocate(int_space)
end select
end do
call nf_close_wrap(out_ncid, subname)
call nf_close_wrap(in_ncid, subname)
end subroutine put_dest_dim_vars
!-----------------------------------------------------------------------------
subroutine gen_weights(grid_var, normalizing, rmask, &
latmin, latmax, lonmin, lonmax)
!--------------------------------------------------------------------------
! arguments
!--------------------------------------------------------------------------
type(POP_grid), intent(in) :: &
grid_var ! source grid
logical(kind=log_kind), intent(in) :: &
normalizing ! compute area averages (.T.) or area sums (.F.)
integer(kind=int_kind), dimension(:,:), intent(in) :: &
rmask ! basin mask
real(kind=r8), intent(in) :: &
latmin, latmax, & ! latitude bounds of averaging
lonmin, lonmax ! longitude bounds of averaging
!--------------------------------------------------------------------------
! local variables
!--------------------------------------------------------------------------
integer(kind=int_kind) :: &
i, j, k, & ! loop indices
j2, & ! loop index into weights for a T cell
j3, & ! loop index into dest_lat_axis
j4, & ! loop index into Twght
jlo, jhi, & ! lo & hi j indices of a T cell intersection with
! latitude bands from dest_lat_axis
Twght_ind0_prev, & ! Twght_ind0 from 'previous' T cell
wght_cnt_prev, & ! wght_cnt from 'previous' T cell
basin_cnt, & ! no. of basins (=maxval(abs(rmask)))
level_cnt, & ! no. of vertical levels (=maxval(kmt))
basin_ind ! basin index of current T cell
real(kind=r8) :: &
tmp_lat, & ! temporary latitude value
scale, & ! scale factor to convert line integral area
! to model generated area
LI_area_band ! line integral area of T cell in latitude band
real(kind=r8), dimension(4) :: &
Tvert_lat, & ! latitude of tracer cell vertices (degrees_north)
Tvert_long ! longitude of tracer cell vertices (degrees_east)
!--------------------------------------------------------------------------
allocate( &
wght_cnt(grid_var%nlon, grid_var%nlat), &
lat_axis_ind0(grid_var%nlon, grid_var%nlat), &
Twght_ind0(grid_var%nlon, grid_var%nlat))
do j = 1,grid_var%nlat
do i = 1,grid_var%nlon
if (grid_var%kmt(i,j) == 0 .or. rmask(i,j) == 0) then
wght_cnt(i,j) = 0
lat_axis_ind0(i,j) = 1
cycle
end if
call gen_Tcell_vertices(grid_var, i, j, Tvert_lat, Tvert_long)
if (i == grid_var%pole_i .and. j == grid_var%pole_j) then
tmp_lat = c90
else
tmp_lat = maxval(Tvert_lat)
end if
if ((tmp_lat <= latmin) .or. &
(minval(Tvert_lat) >= latmax) .or. &
(maxval(Tvert_long) <= lonmin) .or. &
(minval(Tvert_long) >= lonmax)) then
wght_cnt(i, j) = 0
cycle
end if
do jhi = dest_lat_len,1,-1
if (tmp_lat > dest_lat_axis_edges(jhi-1)) exit
end do
tmp_lat = minval(Tvert_lat)
do jlo = 1,dest_lat_len
if (tmp_lat < dest_lat_axis_edges(jlo)) exit
end do
wght_cnt(i,j) = jhi - jlo + 1
lat_axis_ind0(i,j) = jlo
end do
end do
allocate(Twght(sum(wght_cnt)))
Twght_ind0_prev = 1
wght_cnt_prev = 0
do j = 1,grid_var%nlat
do i = 1,grid_var%nlon
Twght_ind0(i,j) = Twght_ind0_prev + wght_cnt_prev
if (wght_cnt(i,j) == 0 .or. grid_var%kmt(i,j) == 0) cycle
call gen_Tcell_vertices(grid_var, i, j, Tvert_lat, Tvert_long)
!--------------------------------------------------------------------
! normalize weights so that the weights for a full cell would
! sum to the model generated area
!--------------------------------------------------------------------
scale = grid_var%tarea(i,j) / clip_sphere_area(4, Tvert_lat, &
Tvert_long, -90.0_r8, 90.0_r8, -720.0_r8, 720.0_r8, &
grid_var%lat_threshold)
do j2 = 0,wght_cnt(i,j)-1
j3 = lat_axis_ind0(i,j) + j2
j4 = Twght_ind0(i,j) + j2
LI_area_band = clip_sphere_area(4, Tvert_lat, Tvert_long, &
dest_lat_axis_edges(j3-1), dest_lat_axis_edges(j3), &
lonmin, lonmax, grid_var%lat_threshold)
Twght(j4) = LI_area_band * scale
end do
Twght_ind0_prev = Twght_ind0(i,j)
wght_cnt_prev = wght_cnt(i,j)
end do
end do
level_cnt = maxval(grid_var%kmt)
basin_cnt = maxval(abs(rmask))
allocate( &
var_in_lev(grid_var%nlon,grid_var%nlat), &
var_out_lev(dest_lat_len, 0:basin_cnt))
allocate(Twght_basin_sum_kmt(dest_lat_len, 0:basin_cnt, level_cnt))
Twght_basin_sum_kmt = c0
do k = 1,level_cnt
do j = 1,grid_var%nlat
do i = 1,grid_var%nlon
if (k > grid_var%kmt(i,j) .or. rmask(i,j) == 0) cycle
basin_ind = abs(rmask(i,j))
do j2 = 0,wght_cnt(i,j)-1
j3 = lat_axis_ind0(i,j) + j2
j4 = Twght_ind0(i,j) + j2
Twght_basin_sum_kmt(j3, basin_ind, k) = &
Twght_basin_sum_kmt(j3, basin_ind, k) + Twght(j4)
end do
end do
end do
end do
do k = 1,level_cnt
do basin_ind = 1,basin_cnt
do j3 = 1,dest_lat_len
Twght_basin_sum_kmt(j3, 0, k) = &
Twght_basin_sum_kmt(j3, 0, k) + &
Twght_basin_sum_kmt(j3, basin_ind, k)
end do
end do
end do
end subroutine gen_weights
!-----------------------------------------------------------------------------
subroutine compute_zon_avg(outfilename, grid_var, infilename, normalizing, &
rmask, time_const)
!--------------------------------------------------------------------------
! arguments
!--------------------------------------------------------------------------
character(len=*), intent(in) :: &
outfilename, & ! name of file where dims & vars are to be defined
infilename ! name of file where dims & vars come from
type(POP_grid), intent(in) :: &
grid_var ! grid to generate destination axis from
logical(kind=log_kind), intent(in) :: &
normalizing ! compute area averages (.T.) or area sums (.F.)
integer(kind=int_kind), dimension(:,:), intent(in) :: &
rmask ! basin mask
logical(kind=log_kind), intent(in) :: &
time_const ! are time constant variables to be averaged?
!--------------------------------------------------------------------------
! local variables
!--------------------------------------------------------------------------
character(len=*), parameter :: &
subname = 'zonal_avg_mod:compute_zon_avg'
integer(kind=int_kind) :: &
stat, & ! status from NetCDF calls
in_ncid, & ! file ID for infile
out_ncid, & ! file ID for outfile
in_nvars, & ! number of variables in infile
in_varid, & ! varid of coordinate variable in infile
out_varid, & ! varid of coordinate variable in outfile
in_varndims, & ! number of dimensions of variable in infile
nlon_dimid, & ! dimid on "nlon" in infile
nlat_dimid, & ! dimid on "nlat" in infile
in_unlimdim, & ! unlimited dimension from infile
dimind, & ! loop index
depthind, & ! which dimension corresponds to depth
swap_int ! temporary for swapping integers
logical(kind=log_kind) :: &
carry ! carry flag for incrementing in_start
integer(kind=int_kind), dimension(NF_MAX_VAR_DIMS) :: &
in_dimids, & ! dimids of in_varid
in_dimlens, & ! dimlens of in_dimids
in_start, & ! start index for reading from infile
in_count, & ! count for reading from infile
out_start, & ! start index for writing to outfile
out_count, & ! count for writing to outfile
in_dimind ! maps dimensions from infile to outfile
real(kind=r8) :: &
msv ! value to be used where there is no data
character(len=char_len) :: &
varname, & ! variable name
dimname ! dimension name
!--------------------------------------------------------------------------
call nf_open_wrap(infilename, 0, in_ncid, subname)
call nf_open_wrap(outfilename, NF_WRITE, out_ncid, subname)
call nf_inq_unlimdim_wrap(in_ncid, in_unlimdim, subname)
call nf_inq_varid_wrap(in_ncid, TLONG_name, in_varid, subname)
call nf_inq_vardimid_wrap(in_ncid, in_varid, in_dimids, subname)
nlon_dimid = in_dimids(1)
nlat_dimid = in_dimids(2)
call nf_inq_nvars_wrap(in_ncid, in_nvars, subname)
do in_varid = 1,in_nvars
if (.not. to_comp_avg(in_ncid, in_varid, time_const)) cycle
varname = ''
call nf_inq_varname_wrap(in_ncid, in_varid, varname, subname)
call nf_inq_varid_wrap(out_ncid, trim(varname), out_varid, subname)
!-----------------------------------------------------------------------
! get _FillValue or missing_value value
!-----------------------------------------------------------------------
call nf_get_att_wrap(out_ncid, out_varid, '_FillValue', msv, subname, &
ALLOW=NF_ENOTATT, stat_out=stat)
if (stat == NF_ENOTATT) then
call nf_get_att_wrap(out_ncid, out_varid, 'missing_value', msv, subname, &
ALLOW=NF_ENOTATT, stat_out=stat)
if (stat == NF_ENOTATT) then
msv = NF_FILL_DOUBLE
end if
end if
call nf_inq_varndims_wrap(in_ncid, in_varid, in_varndims, subname)
call nf_inq_vardimid_wrap(in_ncid, in_varid, in_dimids, subname)
!-----------------------------------------------------------------------
! attempt to detect depth dimension, if not found, use surface mask
! get dimension lengths and initialize in_start & in_count
!-----------------------------------------------------------------------
depthind = -1
do dimind = 1,in_varndims
dimname = ''
call nf_inq_dimname_wrap(in_ncid, in_dimids(dimind), dimname, subname)
if (dimname(1:3) == 'z_t') depthind = dimind
if (dimname(1:3) == 'z_w') depthind = dimind
if (dimname == 'depth') depthind = dimind
if (dimname == 'Depth') depthind = dimind
if (dimname == 'DEPTH') depthind = dimind
call nf_inq_dimlen_wrap(in_ncid, in_dimids(dimind), in_dimlens(dimind), &
subname)
in_start(dimind) = 1
if (in_dimids(dimind) == nlon_dimid) then
in_count(dimind) = in_dimlens(dimind)
else if (in_dimids(dimind) == nlat_dimid) then
in_count(dimind) = in_dimlens(dimind)
else
in_count(dimind) = 1
end if
end do
!-----------------------------------------------------------------------
! generate out_count & in_dimind, first mapping nlat -> basin
! then shift basin to end of array
!-----------------------------------------------------------------------
do dimind = 1,in_varndims
if (in_dimids(dimind) == nlon_dimid) then
out_count(dimind) = dest_lat_len
else if (in_dimids(dimind) == nlat_dimid) then
out_count(dimind) = maxval(abs(rmask)) + 1
else
out_count(dimind) = 1
end if
in_dimind(dimind) = dimind
end do
do dimind = 1,in_varndims-1
if (out_count(dimind) == maxval(abs(rmask)) + 1 .and. &
in_dimids(dimind+1) /= in_unlimdim) then
swap_int = out_count(dimind+1)
out_count(dimind+1) = out_count(dimind)
out_count(dimind) = swap_int
swap_int = in_dimind(dimind+1)
in_dimind(dimind+1) = in_dimind(dimind)
in_dimind(dimind) = swap_int
end if
end do
!-----------------------------------------------------------------------
! loop over non-nlon,nlat dimensions doing the following
! 1) read a level
! 2) compute zonal average of the level
! 3) set out_start from in_start
! 4) write out the average
! 5) advance in_start counters
!-----------------------------------------------------------------------
do
call nf_get_vara_wrap(in_ncid, in_varid, in_start, in_count, &
var_in_lev, subname)
if (depthind == -1) then
call compute_zon_avg_lev(grid_var, normalizing, rmask, 1, msv)
else
call compute_zon_avg_lev(grid_var, normalizing, rmask, &
in_start(depthind), msv)
end if
do dimind = 1,in_varndims
if (in_dimids(in_dimind(dimind)) /= in_unlimdim) then
out_start(dimind) = in_start(in_dimind(dimind))
else
out_start(dimind) = in_start(in_dimind(dimind)) + init_unlimdim_len
end if
end do
call nf_put_vara_wrap(out_ncid, out_varid, out_start, out_count, &
var_out_lev, subname)
carry = .true.
do dimind = 1,in_varndims
if (in_dimids(dimind) /= nlon_dimid .and. &
in_dimids(dimind) /= nlat_dimid .and. carry) then
in_start(dimind) = in_start(dimind) + 1
carry = (in_start(dimind) > in_dimlens(dimind))
if (carry) in_start(dimind) = 1
end if
end do
if (carry) exit
end do
end do
call nf_close_wrap(out_ncid, subname)
call nf_close_wrap(in_ncid, subname)
end subroutine compute_zon_avg
!-----------------------------------------------------------------------------
subroutine compute_zon_avg_lev(grid_var, normalizing, rmask, k, msv)
!--------------------------------------------------------------------------
! arguments
!--------------------------------------------------------------------------
type(POP_grid), intent(in) :: &
grid_var ! source grid
logical(kind=log_kind), intent(in) :: &
normalizing ! compute area averages (.T.) or area sums (.F.)
integer(kind=int_kind), dimension(:,:), intent(in) :: &
rmask ! basin mask
integer(kind=int_kind), intent(in) :: &
k ! depth level
real(kind=r8), intent(in) :: &
msv ! value to use where there is no data
!--------------------------------------------------------------------------
! local variables
!--------------------------------------------------------------------------
integer(kind=int_kind) :: &
i, j, & ! loop indices
j2, & ! loop index into weights for a T cell
j3, & ! loop index into dest_lat_axis
j4, & ! loop index into Twght
basin_cnt, & ! no. of basins (=maxval(abs(rmask)))
basin_ind ! basin index of current T cell
!--------------------------------------------------------------------------
var_out_lev = c0
basin_cnt = maxval(abs(rmask))
do j = 1,grid_var%nlat
do i = 1,grid_var%nlon
if (k > grid_var%kmt(i,j) .or. rmask(i,j) == 0) cycle
basin_ind = abs(rmask(i,j))
do j2 = 0,wght_cnt(i,j)-1
j3 = lat_axis_ind0(i,j) + j2
j4 = Twght_ind0(i,j) + j2
var_out_lev(j3, basin_ind) = var_out_lev(j3, basin_ind) + &
Twght(j4) * var_in_lev(i,j)
end do
end do
end do
do basin_ind = 1,basin_cnt
do j3 = 1,dest_lat_len
var_out_lev(j3, 0) = var_out_lev(j3, 0) + var_out_lev(j3, basin_ind)
end do
end do
if (normalizing) then
do basin_ind = 0,basin_cnt
do j3 = 1,dest_lat_len
if (Twght_basin_sum_kmt(j3, basin_ind, k) > c0) then
var_out_lev(j3, basin_ind) = var_out_lev(j3, basin_ind) / &
Twght_basin_sum_kmt(j3, basin_ind, k)
else
var_out_lev(j3, basin_ind) = msv
end if
end do
end do
else
do basin_ind = 0,basin_cnt
do j3 = 1,dest_lat_len
if (Twght_basin_sum_kmt(j3, basin_ind, k) == c0) then
var_out_lev(j3, basin_ind) = msv
end if
end do
end do
end if
end subroutine compute_zon_avg_lev
!-----------------------------------------------------------------------------
subroutine free_zon_avg
deallocate( &
dest_lat_axis, &
dest_lat_axis_edges, &
wght_cnt, &
lat_axis_ind0, &
Twght_ind0, &
Twght, &
var_in_lev, &
var_out_lev )
if (allocated(Twght_basin_sum_kmt)) deallocate(Twght_basin_sum_kmt)
end subroutine free_zon_avg
!-----------------------------------------------------------------------------
end module zonal_avg_mod
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment