507 lines
15 KiB
Fortran
Executable File
507 lines
15 KiB
Fortran
Executable File
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
|
!
|
|
! This program creates a remapping grid file for Gaussian lat/lon
|
|
! grids (for spectral transform codes).
|
|
!
|
|
!-----------------------------------------------------------------------
|
|
!
|
|
! CVS:$Id: convertgauss.f,v 1.3 2000/04/19 22:05:57 pwjones Exp $
|
|
!
|
|
! Copyright (c) 1997, 1998 the Regents of the University of
|
|
! California.
|
|
!
|
|
! 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.
|
|
!
|
|
!***********************************************************************
|
|
|
|
program convert_gauss
|
|
|
|
!-----------------------------------------------------------------------
|
|
!
|
|
! This file creates a remapping grid file for a Gaussian grid
|
|
!
|
|
!-----------------------------------------------------------------------
|
|
|
|
use kinds_mod
|
|
use constants
|
|
use iounits
|
|
use netcdf_mod
|
|
|
|
implicit none
|
|
|
|
!-----------------------------------------------------------------------
|
|
!
|
|
! variables that describe the grid
|
|
!
|
|
! T42: nx=128 ny=64
|
|
! T62: nx=192 ny=94
|
|
!
|
|
!-----------------------------------------------------------------------
|
|
|
|
integer (kind=int_kind), parameter ::
|
|
& nx = 192, ny = 94,
|
|
& grid_size = nx*ny,
|
|
& grid_rank = 2,
|
|
& grid_corners = 4
|
|
|
|
character(char_len), parameter ::
|
|
& grid_name = 'T62 Gaussian Grid',
|
|
& grid_file_out = 'remap_grid_T62.nc'
|
|
|
|
integer (kind=int_kind), dimension(grid_rank) ::
|
|
& grid_dims
|
|
|
|
!-----------------------------------------------------------------------
|
|
!
|
|
! grid coordinates and masks
|
|
!
|
|
!-----------------------------------------------------------------------
|
|
|
|
integer (kind=int_kind), dimension(grid_size) ::
|
|
& grid_imask
|
|
|
|
real (kind=dbl_kind), dimension(grid_size) ::
|
|
& grid_center_lat, ! lat/lon coordinates for
|
|
& grid_center_lon ! each grid center in degrees
|
|
|
|
real (kind=dbl_kind), dimension(grid_corners,grid_size) ::
|
|
& grid_corner_lat, ! lat/lon coordinates for
|
|
& grid_corner_lon ! each grid corner in degrees
|
|
|
|
!-----------------------------------------------------------------------
|
|
!
|
|
! other local variables
|
|
!
|
|
!-----------------------------------------------------------------------
|
|
|
|
integer (kind=int_kind) :: i, j, iunit, atm_add
|
|
|
|
integer (kind=int_kind) ::
|
|
& ncstat, ! general netCDF status variable
|
|
& nc_grid_id, ! netCDF grid dataset id
|
|
& nc_gridsize_id, ! netCDF grid size dim id
|
|
& nc_gridcorn_id, ! netCDF grid corner dim id
|
|
& nc_gridrank_id, ! netCDF grid rank dim id
|
|
& nc_griddims_id, ! netCDF grid dimension size id
|
|
& nc_grdcntrlat_id, ! netCDF grid center lat id
|
|
& nc_grdcntrlon_id, ! netCDF grid center lon id
|
|
& nc_grdimask_id, ! netCDF grid mask id
|
|
& nc_grdcrnrlat_id, ! netCDF grid corner lat id
|
|
& nc_grdcrnrlon_id ! netCDF grid corner lon id
|
|
|
|
integer (kind=int_kind), dimension(2) ::
|
|
& nc_dims2_id ! netCDF dim id array for 2-d arrays
|
|
|
|
real (kind=dbl_kind) :: dlon, minlon, maxlon, centerlon,
|
|
& minlat, maxlat, centerlat
|
|
|
|
real (kind=dbl_kind), dimension(ny) :: gauss_root, gauss_wgt
|
|
|
|
!-----------------------------------------------------------------------
|
|
!
|
|
! compute longitudes of cell centers and corners. set up alon
|
|
! array for search routine.
|
|
!
|
|
!-----------------------------------------------------------------------
|
|
|
|
grid_dims(1) = nx
|
|
grid_dims(2) = ny
|
|
|
|
dlon = 360./nx
|
|
|
|
do i=1,nx
|
|
|
|
centerlon = (i-1)*dlon
|
|
minlon = centerlon - half*dlon
|
|
maxlon = centerlon + half*dlon
|
|
|
|
do j=1,ny
|
|
atm_add = (j-1)*nx + i
|
|
|
|
grid_center_lon(atm_add ) = centerlon
|
|
grid_corner_lon(1,atm_add) = minlon
|
|
grid_corner_lon(2,atm_add) = maxlon
|
|
grid_corner_lon(3,atm_add) = maxlon
|
|
grid_corner_lon(4,atm_add) = minlon
|
|
end do
|
|
|
|
end do
|
|
|
|
!-----------------------------------------------------------------------
|
|
!
|
|
! compute Gaussian latitudes and store in gauss_wgt.
|
|
!
|
|
!-----------------------------------------------------------------------
|
|
|
|
call gquad(ny, gauss_root, gauss_wgt)
|
|
do j=1,ny
|
|
gauss_wgt(j) = pih - gauss_root(ny+1-j)
|
|
end do
|
|
|
|
!-----------------------------------------------------------------------
|
|
!
|
|
! compute latitudes at cell centers and corners. set up alat
|
|
! array for search routine.
|
|
!
|
|
!-----------------------------------------------------------------------
|
|
|
|
do j=1,ny
|
|
centerlat = gauss_wgt(j)
|
|
|
|
if (j .eq. 1) then
|
|
minlat = -pih
|
|
else
|
|
minlat = ATAN((COS(gauss_wgt(j-1)) -
|
|
& COS(gauss_wgt(j )))/
|
|
& (SIN(gauss_wgt(j )) -
|
|
& SIN(gauss_wgt(j-1))))
|
|
endif
|
|
|
|
if (j .eq. ny) then
|
|
maxlat = pih
|
|
else
|
|
maxlat = ATAN((COS(gauss_wgt(j )) -
|
|
& COS(gauss_wgt(j+1)))/
|
|
& (SIN(gauss_wgt(j+1)) -
|
|
& SIN(gauss_wgt(j ))))
|
|
endif
|
|
|
|
do i=1,nx
|
|
atm_add = (j-1)*nx + i
|
|
grid_center_lat(atm_add ) = centerlat*360./pi2
|
|
grid_corner_lat(1,atm_add) = minlat*360./pi2
|
|
grid_corner_lat(2,atm_add) = minlat*360./pi2
|
|
grid_corner_lat(3,atm_add) = maxlat*360./pi2
|
|
grid_corner_lat(4,atm_add) = maxlat*360./pi2
|
|
end do
|
|
|
|
end do
|
|
|
|
!-----------------------------------------------------------------------
|
|
!
|
|
! define mask
|
|
!
|
|
!-----------------------------------------------------------------------
|
|
|
|
grid_imask = 1
|
|
|
|
!-----------------------------------------------------------------------
|
|
!
|
|
! set up attributes for netCDF file
|
|
!
|
|
!-----------------------------------------------------------------------
|
|
|
|
!***
|
|
!*** create netCDF dataset for this grid
|
|
!***
|
|
|
|
ncstat = nf_create (grid_file_out, NF_CLOBBER,
|
|
& nc_grid_id)
|
|
call netcdf_error_handler(ncstat)
|
|
|
|
ncstat = nf_put_att_text (nc_grid_id, NF_GLOBAL, 'title',
|
|
& len_trim(grid_name), grid_name)
|
|
call netcdf_error_handler(ncstat)
|
|
|
|
!***
|
|
!*** define grid size dimension
|
|
!***
|
|
|
|
ncstat = nf_def_dim (nc_grid_id, 'grid_size', grid_size,
|
|
& nc_gridsize_id)
|
|
call netcdf_error_handler(ncstat)
|
|
|
|
!***
|
|
!*** define grid corner dimension
|
|
!***
|
|
|
|
ncstat = nf_def_dim (nc_grid_id, 'grid_corners', grid_corners,
|
|
& nc_gridcorn_id)
|
|
call netcdf_error_handler(ncstat)
|
|
|
|
!***
|
|
!*** define grid rank dimension
|
|
!***
|
|
|
|
ncstat = nf_def_dim (nc_grid_id, 'grid_rank', grid_rank,
|
|
& nc_gridrank_id)
|
|
call netcdf_error_handler(ncstat)
|
|
|
|
!***
|
|
!*** define grid dimension size array
|
|
!***
|
|
|
|
ncstat = nf_def_var (nc_grid_id, 'grid_dims', NF_INT,
|
|
& 1, nc_gridrank_id, nc_griddims_id)
|
|
call netcdf_error_handler(ncstat)
|
|
|
|
!***
|
|
!*** define grid center latitude array
|
|
!***
|
|
|
|
ncstat = nf_def_var (nc_grid_id, 'grid_center_lat', NF_DOUBLE,
|
|
& 1, nc_gridsize_id, nc_grdcntrlat_id)
|
|
call netcdf_error_handler(ncstat)
|
|
|
|
ncstat = nf_put_att_text (nc_grid_id, nc_grdcntrlat_id, 'units',
|
|
& 7, 'degrees')
|
|
call netcdf_error_handler(ncstat)
|
|
|
|
!***
|
|
!*** define grid center longitude array
|
|
!***
|
|
|
|
ncstat = nf_def_var (nc_grid_id, 'grid_center_lon', NF_DOUBLE,
|
|
& 1, nc_gridsize_id, nc_grdcntrlon_id)
|
|
call netcdf_error_handler(ncstat)
|
|
|
|
ncstat = nf_put_att_text (nc_grid_id, nc_grdcntrlon_id, 'units',
|
|
& 7, 'degrees')
|
|
call netcdf_error_handler(ncstat)
|
|
|
|
!***
|
|
!*** define grid mask
|
|
!***
|
|
|
|
ncstat = nf_def_var (nc_grid_id, 'grid_imask', NF_INT,
|
|
& 1, nc_gridsize_id, nc_grdimask_id)
|
|
call netcdf_error_handler(ncstat)
|
|
|
|
ncstat = nf_put_att_text (nc_grid_id, nc_grdimask_id, 'units',
|
|
& 8, 'unitless')
|
|
call netcdf_error_handler(ncstat)
|
|
|
|
!***
|
|
!*** define grid corner latitude array
|
|
!***
|
|
|
|
nc_dims2_id(1) = nc_gridcorn_id
|
|
nc_dims2_id(2) = nc_gridsize_id
|
|
|
|
ncstat = nf_def_var (nc_grid_id, 'grid_corner_lat', NF_DOUBLE,
|
|
& 2, nc_dims2_id, nc_grdcrnrlat_id)
|
|
call netcdf_error_handler(ncstat)
|
|
|
|
ncstat = nf_put_att_text (nc_grid_id, nc_grdcrnrlat_id, 'units',
|
|
& 7, 'degrees')
|
|
call netcdf_error_handler(ncstat)
|
|
|
|
!***
|
|
!*** define grid corner longitude array
|
|
!***
|
|
|
|
ncstat = nf_def_var (nc_grid_id, 'grid_corner_lon', NF_DOUBLE,
|
|
& 2, nc_dims2_id, nc_grdcrnrlon_id)
|
|
call netcdf_error_handler(ncstat)
|
|
|
|
ncstat = nf_put_att_text (nc_grid_id, nc_grdcrnrlon_id, 'units',
|
|
& 7, 'degrees')
|
|
call netcdf_error_handler(ncstat)
|
|
|
|
!***
|
|
!*** end definition stage
|
|
!***
|
|
|
|
ncstat = nf_enddef(nc_grid_id)
|
|
call netcdf_error_handler(ncstat)
|
|
|
|
!-----------------------------------------------------------------------
|
|
!
|
|
! write grid data
|
|
!
|
|
!-----------------------------------------------------------------------
|
|
|
|
ncstat = nf_put_var_int(nc_grid_id, nc_griddims_id, grid_dims)
|
|
call netcdf_error_handler(ncstat)
|
|
|
|
ncstat = nf_put_var_int(nc_grid_id, nc_grdimask_id, grid_imask)
|
|
call netcdf_error_handler(ncstat)
|
|
|
|
ncstat = nf_put_var_double(nc_grid_id, nc_grdcntrlat_id,
|
|
& grid_center_lat)
|
|
call netcdf_error_handler(ncstat)
|
|
|
|
ncstat = nf_put_var_double(nc_grid_id, nc_grdcntrlon_id,
|
|
& grid_center_lon)
|
|
call netcdf_error_handler(ncstat)
|
|
|
|
ncstat = nf_put_var_double(nc_grid_id, nc_grdcrnrlat_id,
|
|
& grid_corner_lat)
|
|
call netcdf_error_handler(ncstat)
|
|
|
|
ncstat = nf_put_var_double(nc_grid_id, nc_grdcrnrlon_id,
|
|
& grid_corner_lon)
|
|
call netcdf_error_handler(ncstat)
|
|
|
|
ncstat = nf_close(nc_grid_id)
|
|
call netcdf_error_handler(ncstat)
|
|
|
|
!-----------------------------------------------------------------------
|
|
|
|
end program convert_gauss
|
|
|
|
!***********************************************************************
|
|
|
|
subroutine gquad(l,root,w)
|
|
|
|
!-----------------------------------------------------------------------
|
|
!
|
|
! This subroutine finds the l roots (in theta) and gaussian weights
|
|
! associated with the legendre polynomial of degree l > 1.
|
|
!
|
|
!-----------------------------------------------------------------------
|
|
|
|
use kinds_mod
|
|
use constants
|
|
|
|
implicit none
|
|
|
|
!-----------------------------------------------------------------------
|
|
!
|
|
! intent(in)
|
|
!
|
|
!-----------------------------------------------------------------------
|
|
|
|
integer (kind=int_kind), intent(in) :: l
|
|
|
|
!-----------------------------------------------------------------------
|
|
!
|
|
! intent(out)
|
|
!
|
|
!-----------------------------------------------------------------------
|
|
|
|
real (kind=dbl_kind), dimension(l), intent(out) ::
|
|
& root, w
|
|
|
|
!-----------------------------------------------------------------------
|
|
!
|
|
! local
|
|
!
|
|
!-----------------------------------------------------------------------
|
|
|
|
integer (kind=int_kind) :: l1, l2, l22, l3, k, i, j
|
|
|
|
real (kind=dbl_kind) ::
|
|
& del,co,p1,p2,p3,t1,t2,slope,s,c,pp1,pp2,p00
|
|
|
|
!-----------------------------------------------------------------------
|
|
!
|
|
! Define useful constants.
|
|
!
|
|
!-----------------------------------------------------------------------
|
|
|
|
del= pi/float(4*l)
|
|
l1 = l+1
|
|
co = float(2*l+3)/float(l1**2)
|
|
p2 = 1.0
|
|
t2 = -del
|
|
l2 = l/2
|
|
k = 1
|
|
p00 = one/sqrt(two)
|
|
|
|
!-----------------------------------------------------------------------
|
|
!
|
|
! Start search for each root by looking for crossing point.
|
|
!
|
|
!-----------------------------------------------------------------------
|
|
|
|
do i=1,l2
|
|
10 t1 = t2
|
|
t2 = t1+del
|
|
p1 = p2
|
|
s = sin(t2)
|
|
c = cos(t2)
|
|
pp1 = 1.0
|
|
p3 = p00
|
|
do j=1,l1
|
|
pp2 = pp1
|
|
pp1 = p3
|
|
p3 = 2.0*sqrt((float(j**2)-0.250)/float(j**2))*c*pp1-
|
|
& sqrt(float((2*j+1)*(j-1)*(j-1))/
|
|
& float((2*j-3)*j*j))*pp2
|
|
end do
|
|
p2 = pp1
|
|
if ((k*p2).gt.0) goto 10
|
|
|
|
!-----------------------------------------------------------------------
|
|
!
|
|
! Now converge using Newton-Raphson.
|
|
!
|
|
!-----------------------------------------------------------------------
|
|
|
|
k = -k
|
|
20 continue
|
|
slope = (t2-t1)/(p2-p1)
|
|
t1 = t2
|
|
t2 = t2-slope*p2
|
|
p1 = p2
|
|
s = sin(t2)
|
|
c = cos(t2)
|
|
pp1 = 1.0
|
|
p3 = p00
|
|
do j=1,l1
|
|
pp2 = pp1
|
|
pp1 = p3
|
|
p3 = 2.0*sqrt((float(j**2)-0.250)/float(j**2))*c*pp1-
|
|
& sqrt(float((2*j+1)*(j-1)*(j-1))/
|
|
& float((2*j-3)*j*j))*pp2
|
|
end do
|
|
p2 = pp1
|
|
if (abs(p2).gt.1.e-10) goto 20
|
|
root(i) = t2
|
|
w(i) = co*(sin(t2)/p3)**2
|
|
end do
|
|
|
|
!-----------------------------------------------------------------------
|
|
!
|
|
! If l is odd, take care of odd point.
|
|
!
|
|
!-----------------------------------------------------------------------
|
|
|
|
l22 = 2*l2
|
|
if (l22 .ne. l) then
|
|
l2 = l2+1
|
|
t2 = pi/2.0
|
|
root(l2) = t2
|
|
s = sin(t2)
|
|
c = cos(t2)
|
|
pp1 = 1.0
|
|
p3 = p00
|
|
do j=1,l1
|
|
pp2 = pp1
|
|
pp1 = p3
|
|
p3 = 2.0*sqrt((float(j**2)-0.250)/float(j**2))*c*pp1-
|
|
& sqrt(float((2*j+1)*(j-1)*(j-1))/
|
|
& float((2*j-3)*j*j))*pp2
|
|
end do
|
|
p2 = pp1
|
|
w(l2) = co/p3**2
|
|
endif
|
|
|
|
!-----------------------------------------------------------------------
|
|
!
|
|
! Use symmetry to compute remaining roots and weights.
|
|
!
|
|
!-----------------------------------------------------------------------
|
|
|
|
l3 = l2+1
|
|
do i=l3,l
|
|
root(i) = pi-root(l-i+1)
|
|
w(i) = w(l-i+1)
|
|
end do
|
|
|
|
!-----------------------------------------------------------------------
|
|
|
|
end subroutine gquad
|
|
|
|
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|