diff --git a/src/ocnicepost.fd/arrays_mod.F90 b/src/ocnicepost.fd/arrays_mod.F90 new file mode 100644 index 00000000..b6e58eb6 --- /dev/null +++ b/src/ocnicepost.fd/arrays_mod.F90 @@ -0,0 +1,108 @@ +module arrays_mod + + use init_mod , only : nxt, nyt, nlevs, nxr, nyr + use init_mod , only : debug, logunit + use init_mod , only : vardefs + + implicit none + + integer :: nbilin2d !< the number of 2D fields mapped bilinearly + integer :: nbilin3d !< the number of 3D fields mapped bilinearly + integer :: nconsd2d !< the number of 2D fields mapped conservatively + + ! source arrays + real, allocatable, dimension(:,:) :: bilin2d !< packed 2D source fields for bilinear remap + real, allocatable, dimension(:,:) :: consd2d !< packed 2D source fields for conservative remap + real, allocatable, dimension(:,:,:) :: bilin3d !< packed 3D source fields for bilinear remap + + type(vardefs), allocatable, dimension(:) :: b2d !< variable metadata for 2D source fields bilinear remap + type(vardefs), allocatable, dimension(:) :: c2d !< variable metadata for 2D source fields conservative remap + type(vardefs), allocatable, dimension(:) :: b3d !< variable metadata for 3D source fields bilinear remap + + ! destination arrays + real, allocatable, dimension(:,:) :: rgb2d !< packed 2D fields with bilinear remap + real, allocatable, dimension(:,:) :: rgc2d !< packed 2D fields with conservative remap + real, allocatable, dimension(:,:,:) :: rgb3d !< packed 3D fields with bilinear remap + + real, allocatable, dimension(:,:) :: dstlon !< the destination grid longitudes + real, allocatable, dimension(:,:) :: dstlat !< the destination grid latitudes + + public setup_packing + +contains + + subroutine setup_packing(nvalid, vars) + + type(vardefs), intent(inout) :: vars(:) + integer , intent(in) :: nvalid + + ! local variables + integer :: n,i,j,k + + ! -------------------------------------------------------- + ! count numbers of fields to remapped for each + ! mapping type; these can be remapped as packed arrays + ! -------------------------------------------------------- + + nbilin2d = 0; nbilin3d = 0; nconsd2d = 0 + do n = 1,nvalid + if (trim(vars(n)%var_remapmethod) == 'bilinear') then + if (vars(n)%var_dimen == 2) nbilin2d = nbilin2d + 1 + if (vars(n)%var_dimen == 3) nbilin3d = nbilin3d + 1 + end if + if (trim(vars(n)%var_remapmethod) == 'conserve')nconsd2d = nconsd2d + 1 !no 3d variables w/ conservative mapping + end do + if (debug) write(logunit,'(3(a,i4))')'bilin 2d ',nbilin2d,' bilin 3d ',nbilin3d,' conserv 2d ',nconsd2d + + ! initialization required when compiled with sinit_arrays=nan + if (nbilin2d > 0) then + allocate(bilin2d(nxt*nyt,nbilin2d)); bilin2d = 0.0 + allocate(b2d(1:nbilin2d)) + if (debug) write(logunit,'(a)')'allocate bilin2d fields and types ' + end if + if (nconsd2d > 0) then + allocate(consd2d(nxt*nyt,nconsd2d)); consd2d = 0.0 + allocate(c2d(1:nconsd2d)) + if (debug) write(logunit,'(a)')'allocate consd2d fields and types ' + end if + if (nbilin3d > 0) then + allocate(bilin3d(nxt*nyt,nlevs,nbilin3d)); bilin3d = 0.0 + allocate(b3d(1:nbilin3d)) + if (debug) write(logunit,'(a)')'allocate bilin3d fields and types ' + end if + + ! -------------------------------------------------------- + ! create types for each packed array and fill values + ! -------------------------------------------------------- + + i = 0; j = 0; k = 0 + do n = 1,nvalid + if (trim(vars(n)%var_remapmethod) == 'bilinear') then + if (vars(n)%var_dimen == 2 .and. allocated(b2d)) then + i = i+1; b2d(i) = vars(n) + end if + if (vars(n)%var_dimen == 3 .and. allocated(b3d)) then + j = j+1; b3d(j) = vars(n) + end if + end if + if (trim(vars(n)%var_remapmethod) == 'conserve' .and. allocated(c2d)) then + k = k+1; c2d(k) = vars(n) + end if + end do + + ! -------------------------------------------------------- + ! create arrays for remapped packed fields + ! -------------------------------------------------------- + + if (nbilin2d > 0) then + allocate(rgb2d(nxr*nyr,nbilin2d)); rgb2d = 0.0 + end if + if (nconsd2d > 0) then + allocate(rgc2d(nxr*nyr,nconsd2d)); rgc2d = 0.0 + end if + if (nbilin3d > 0) then + allocate(rgb3d(nxr*nyr,nlevs,nbilin3d)); rgb3d = 0.0 + end if + + end subroutine setup_packing +end module arrays_mod diff --git a/src/ocnicepost.fd/init_mod.F90 b/src/ocnicepost.fd/init_mod.F90 new file mode 100644 index 00000000..39631f59 --- /dev/null +++ b/src/ocnicepost.fd/init_mod.F90 @@ -0,0 +1,158 @@ +module init_mod + + implicit none + + public + + real, parameter :: maxvars = 50 !< The maximum number of fields expected in a source file + + type :: vardefs + character(len= 10) :: var_name !< A variable's variable name + character(len=120) :: long_name !< A variable's long name + character(len= 20) :: units !< A variable's unit + character(len= 10) :: var_remapmethod !< A variable's mapping method + integer :: var_dimen !< A variable's dimensionality + character(len= 4) :: var_grid !< A variable's input grid location; all output locations are on cell centers + character(len= 10) :: var_pair !< A variable's pair + character(len= 4) :: var_pair_grid !< A pair variable grid + real :: var_fillvalue !< A variable's fillvalue + end type vardefs + + type(vardefs) :: outvars(maxvars) !< An empty structure filled by reading a csv file describing the fields + + character(len=10) :: ftype !< The type of tripole grid (ocean or ice) + character(len=10) :: fsrc !< A character string for tripole grid + character(len=10) :: fdst !< A character string for the destination grid + character(len=120) :: wgtsdir !< The directory containing the regridding weights + character(len=20) :: input_file !< The input file name + character(len=10) :: maskvar !< The variable in the source file used to create the interpolation mask + + ! rotation angles + character(len=10) :: sinvar !< The variable in the source file containing the cosine of the rotation angle (ocean only) + character(len=10) :: cosvar !< The variable in the source file containing the sine of the rotation angle (ocean only) + character(len=10) :: angvar !< The variable in the source file containing the rotation angle (ice only) + + integer :: nxt !< The x-dimension of the source tripole grid + integer :: nyt !< The y-dimension of the source tripole grid + integer :: nlevs !< The vertical dimension of the source tripole grid + + integer :: nxr !< The x-dimension of the destination rectilinear grid + integer :: nyr !< The y-dimension of the destination rectilinear grid + + integer :: logunit !< The log unit + logical :: debug !< If true, print debug messages and intermediate files + logical :: do_ocnpost !< If true, the source file is ocean, otherwise ice + +contains + + subroutine readnml + + ! local variable + character(len=40) :: fname + integer :: ierr, iounit + integer :: srcdims(3), dstdims(2) + + namelist /ocnicepost_nml/ ftype, srcdims, wgtsdir, dstdims, maskvar, sinvar, cosvar, & + angvar, debug + + ! -------------------------------------------------------- + ! read the name list + ! -------------------------------------------------------- + + srcdims = 0; dstdims = 0 + sinvar=''; cosvar=''; angvar='' + + fname = 'ocnicepost.nml' + inquire (file=trim(fname), iostat=ierr) + if (ierr /= 0) then + write (6, '(3a)') 'Error: input file "', trim(fname), '" does not exist.' + stop + end if + + ! Open and read namelist file. + open (action='read', file=trim(fname), iostat=ierr, newunit=iounit) + read (nml=ocnicepost_nml, iostat=ierr, unit=iounit) + if (ierr /= 0) then + write (6, '(a)') 'Error: invalid namelist format.' + end if + close (iounit) + nxt = srcdims(1); nyt = srcdims(2) + nxr = dstdims(1); nyr = dstdims(2) + If (srcdims(3) > 0) nlevs = srcdims(3) + + ! initialize the source file type and variables + if (trim(ftype) == 'ocean') then + do_ocnpost = .true. + else + do_ocnpost = .false. + end if + input_file = trim(ftype)//'.nc' + + open(newunit=logunit, file=trim(ftype)//'.post.log',form='formatted') + if (debug) write(logunit, '(a)')'input file: '//trim(input_file) + + ! set grid names + fsrc = '' + if (nxt == 1440 .and. nyt == 1080) fsrc = 'mx025' ! 1/4deg tripole + if (nxt == 720 .and. nyt == 576) fsrc = 'mx050' ! 1/2deg tripole + if (nxt == 360 .and. nyt == 320) fsrc = 'mx100' ! 1deg tripole + if (nxt == 72 .and. nyt == 35) fsrc = 'mx500' ! 5deg tripole + if (len_trim(fsrc) == 0) then + write(logunit,'(a)')'source grid dimensions unknown' + stop + end if + + fdst = '' + if (nxr == 1440 .and. nyr == 721) fdst = '0p25' ! 1/4deg rectilinear + if (nxr == 720 .and. nyr == 361) fdst = '0p5' ! 1/2 deg rectilinear + if (nxr == 360 .and. nyr == 181) fdst = '1p0' ! 1 deg rectilinear + if (len_trim(fdst) == 0) then + write(logunit,'(a)')'destination grid dimensions unknown' + stop + end if + + !TODO: test for consistency of source/destination resolution + end subroutine readnml + + subroutine readcsv(nvalid) + + integer, intent(out) :: nvalid + + character(len= 40) :: fname + character(len=100) :: chead + character(len= 10) :: c1,c3,c4,c5,c6 + integer :: i2 + integer :: nn,n,ierr,iounit + + ! -------------------------------------------------------- + ! Open and read list of variables + ! -------------------------------------------------------- + + fname=trim(ftype)//'.csv' + open(newunit=iounit, file=trim(fname), status='old', iostat=ierr) + if (ierr /= 0) then + write (6, '(3a)') 'Error: input file "', trim(fname), '" does not exist.' + stop + end if + + read(iounit,*)chead + nn=0 + do n = 1,maxvars + read(iounit,*,iostat=ierr)c1,i2,c3,c4,c5,c6 + if (ierr .ne. 0) exit + if (len_trim(c1) > 0) then + nn = nn+1 + outvars(nn)%var_name = trim(c1) + outvars(nn)%var_dimen = i2 + outvars(nn)%var_grid = trim(c3) + outvars(nn)%var_remapmethod = trim(c4) + outvars(nn)%var_pair = trim(c5) + outvars(nn)%var_pair_grid = trim(c6) + end if + end do + close(iounit) + nvalid = nn + + end subroutine readcsv + +end module init_mod diff --git a/src/ocnicepost.fd/masking_mod.F90 b/src/ocnicepost.fd/masking_mod.F90 new file mode 100644 index 00000000..a2999596 --- /dev/null +++ b/src/ocnicepost.fd/masking_mod.F90 @@ -0,0 +1,161 @@ +module masking_mod + + use netcdf + + use init_mod , only : nxt, nyt, nlevs, nxr, nyr, vardefs + use init_mod , only : wgtsdir, ftype, fsrc, fdst, input_file + use init_mod , only : do_ocnpost, debug, logunit, maskvar + use arrays_mod , only : dstlat, dstlon + use utils_mod , only : getfield, dumpnc, remap + + implicit none + + real, allocatable, dimension(:,:) :: mask3d !< the 3D mask of the source fields on Ct grid points + real, allocatable, dimension(:) :: mask2d !< the 2D mask of the source fields on Ct grid points + + real, allocatable, dimension(:,:) :: rgmask3d !< the 3D mask of the destination fields on Ct grid points + real, allocatable, dimension(:) :: rgmask2d !< the 2D mask of the destination fields on Ct grid points + + public remap_masks + +contains + + subroutine remap_masks(vfill) + + real, intent(out) :: vfill + + ! local variables + integer :: rc, ncid, varid, n + character(len=240) :: wgtsfile + real :: minlat = -79.75 + + real, allocatable, dimension(:) :: out1d + + ! -------------------------------------------------------- + ! obtain the destination lat and lon directly from the weights file + ! -------------------------------------------------------- + + wgtsfile = trim(wgtsdir)//'tripole.'//trim(fsrc)//'.Ct.to.rect.'//trim(fdst)//'.bilinear.nc' + + allocate(dstlon(nxr,nyr)); dstlon = 0.0 + allocate(dstlat(nxr,nyr)); dstlat = 0.0 + allocate(out1d(nxr*nyr)); out1d = 0.0 + + rc = nf90_open(trim(wgtsfile), nf90_nowrite, ncid) + rc = nf90_inq_varid(ncid, 'xc_b', varid) + rc = nf90_get_var(ncid, varid, out1d) + dstlon = reshape(out1d,(/nxr,nyr/)) + rc = nf90_inq_varid(ncid, 'yc_b', varid) + rc = nf90_get_var(ncid, varid, out1d) + dstlat = reshape(out1d,(/nxr,nyr/)) + rc = nf90_close(ncid) + + ! -------------------------------------------------------- + ! mask is a 2d (ice) or 3d (ocn) array which contains 1's + ! on land and 0's at valid points. + ! when remapped, any mask value > 0 identifies land values that + ! have crept into the field. remapped model fields are then + ! masked with this interpolation mask + ! -------------------------------------------------------- + + if (do_ocnpost) then + call makemask3d(vfill) + else + call makemask2d(vfill) + end if + + ! -------------------------------------------------------- + ! remap the source grid 2/3D mask to obtain the interpolation mask. + ! -------------------------------------------------------- + + if (do_ocnpost) then + call remap(trim(wgtsfile), dim2=nlevs, src_field=mask3d, dst_field=rgmask3d) + ! set interpolation mask missing on land, 1.0 on ocean on destination grids + where(rgmask3d > 0.0)rgmask3d = vfill + where(rgmask3d /= vfill)rgmask3d = 1.0 + ! out1d contains dstlat + do n = 1,nlevs + where(out1d(:) <= minlat)rgmask3d(:,n) = vfill + end do + + if (debug) then + write(logunit,'(a,2g14.4)')'mask min/max on destination grid ',minval(rgmask3d),maxval(rgmask3d) + call dumpnc(trim(ftype)//'.'//trim(fdst)//'.rgmask3d.nc', 'rgmask3d', dims=(/nxr,nyr,nlevs/), & + field=rgmask3d) + end if + else + call remap(trim(wgtsfile), src_field=mask2d, dst_field=rgmask2d) + ! set interpolation mask missing on land, 1.0 on ocean on destination grids + where(rgmask2d > 0.0)rgmask2d = vfill + where(rgmask2d /= vfill)rgmask2d = 1.0 + ! out1d contains dstlat + where(out1d(:) <= minlat)rgmask2d(:) = vfill + + if (debug) then + write(logunit,'(a,2g14.4)')'mask min/max on destination grid ',minval(rgmask2d),maxval(rgmask2d) + call dumpnc(trim(ftype)//'.'//trim(fdst)//'.rgmask2d.nc', 'rgmask2d', dims=(/nxr,nyr/), & + field=rgmask2d) + end if + end if + + end subroutine remap_masks + + subroutine makemask3d(vfill) + + real, intent(out) :: vfill + + ! local variables + integer :: rc, ncid, varid + real, allocatable, dimension(:,:,:) :: tmp3d + + allocate(tmp3d(nxt,nyt,nlevs)); tmp3d = 0.0 + + rc = nf90_open(trim(input_file), nf90_nowrite, ncid) + ! Obtain maskvar directly from file to set fill value + rc = nf90_inq_varid(ncid, trim(maskvar), varid) + rc = nf90_get_att(ncid, varid, '_FillValue', vfill) + rc = nf90_get_var(ncid, varid, tmp3d) + rc = nf90_close(ncid) + + mask3d = reshape(tmp3d, (/nxt*nyt,nlevs/)) + ! set mask3d to 0 on ocean, 1 on land on source grid + where(mask3d .eq. vfill)mask3d = 1.0 + where(mask3d .ne. 1.0)mask3d = 0.0 + + if (debug) then + write(logunit,'(a,2g14.4)')'mask3d min/max on source grid ',minval(mask3d),maxval(mask3d) + call dumpnc(trim(ftype)//'.mask3d.nc', 'mask3d', dims=(/nxt,nyt,nlevs/), field=mask3d) + end if + + end subroutine makemask3d + + subroutine makemask2d(vfill) + + real, intent(out) :: vfill + + ! local variables + integer :: rc, ncid, varid + real, allocatable, dimension(:,:) :: tmp2d + + allocate(tmp2d(nxt,nyt)); tmp2d = 0.0 + + rc = nf90_open(trim(input_file), nf90_nowrite, ncid) + ! Obtain maskvar directly from file to set fill value + rc = nf90_inq_varid(ncid, trim(maskvar), varid) + rc = nf90_get_att(ncid, varid, '_FillValue', vfill) + rc = nf90_get_var(ncid, varid, tmp2d) + rc = nf90_close(ncid) + + mask2d = reshape(tmp2d, (/nxt*nyt/)) + ! set mask2d to 0 on ocean, 1 on land on source grid + mask2d = mask2d - 1.0 + where(mask2d .eq. -1.0)mask2d = 1.0 + + if (debug) then + write(logunit,'(a,2g14.4)')'mask2d min/max on source grid ',minval(mask2d),maxval(mask2d) + call dumpnc(trim(ftype)//'.mask2d.nc', 'mask2d', dims=(/nxt,nyt/), field=mask2d) + end if + + end subroutine makemask2d + +end module masking_mod diff --git a/src/ocnicepost.fd/ocnicepost.F90 b/src/ocnicepost.fd/ocnicepost.F90 new file mode 100644 index 00000000..5669f5d0 --- /dev/null +++ b/src/ocnicepost.fd/ocnicepost.F90 @@ -0,0 +1,371 @@ +program ocnicepost + + ! This program will remap MOM6 ocean or CICE6 ice output on the tripole grid to a + ! rectilinear grid using pre-computed ESMF weights to remap the chosen fields to the + ! destination grid and write the results to a new netCDF file. The ESMF weights needed + ! are pre-generated using the UFS-UTILS cpld_gridgen utility described at + ! https://ufs-community.github.io/UFS_UTILS/cpld_gridgen/index.html. Weights are currently + ! generated only for mapping a source grid to a similar or lower resolution rectilinear grid. + ! + ! Two control files determine the code behaviour. A name list file ocnicepost.nml is used + ! to set the input file type (ice or ocean), the location of the ESMF weights, the dimensions + ! of the source and destination grids, the source variable to used to create an interpolation + ! mask, the required rotation variables and a flag to produce debugging output. + ! + ! The list of variables to be remapped is expected in either ocean.csv or ice.csv. Each + ! variable is specified by name, dimensionality, and the required mapping method. For vectors, + ! the paired vector and it's grid is also listed. + ! + ! Either a 2D (ice) or 3D (ocean) interpolation mask is used to mask remapped fields so that + ! only valid ocean grid points on the source grid appear in the remapped fields. + ! + ! Source fields are packed by mapping method and the remapped. Vector fields are first remapped + ! to the center (Ct) grid points and rotated from the I-J orientation to E-W before remapping. + ! The interpolation mask is used to masks out land contaminated points prior to writing the + ! output netCDF file. + + use netcdf + use init_mod , only : nxt, nyt, nlevs, nxr, nyr, outvars, readnml, readcsv + use init_mod , only : wgtsdir, ftype, fsrc, fdst, input_file, cosvar, sinvar, angvar + use init_mod , only : do_ocnpost, debug, logunit + use arrays_mod , only : b2d, c2d, b3d, rgb2d, rgc2d, rgb3d, dstlon, dstlat, setup_packing + use arrays_mod , only : nbilin2d, nbilin3d, nconsd2d, bilin2d, bilin3d, consd2d + use masking_mod, only : mask2d, mask3d, rgmask2d, rgmask3d, remap_masks + use utils_mod , only : getfield, packarrays, remap, dumpnc + + implicit none + + character(len=120) :: wgtsfile + character(len=120) :: fout + + ! dimensions, units and variables from source file used in creation of + ! output netcdf + real, allocatable, dimension(:) :: z_l !< the vertical grid center + real, allocatable, dimension(:) :: z_i !< the vertical grid interfaces + real, allocatable, dimension(:) :: cosrot, sinrot !< the cos and sin of the rotation angle at Ct points (ocean) + real, allocatable, dimension(:) :: anglet !< the rotation angle at the Ct points (ice) + + ! work arrays for output netcdf + real, allocatable, dimension(:,:) :: out2d !< 2D destination grid output array + real, allocatable, dimension(:,:,:) :: out3d !< 3D destination grid output array + + real(kind=8) :: timestamp + character(len= 40) :: timeunit, timecal + character(len= 20) :: vname, vunit + character(len=120) :: vlong + + real :: vfill + integer :: nvalid + integer :: n,rc,ncid,varid + integer :: idimid,jdimid,kdimid,edimid,timid + integer :: idx1,idx2,idx3 + + ! -------------------------------------------------------- + ! read the nml file and a file containing the list of + ! variables to be remapped + ! -------------------------------------------------------- + + call readnml + call readcsv(nvalid) + + ! -------------------------------------------------------- + ! read the source file and obtain the units and long name, + ! rotation angles, vertical grid and time axis + ! -------------------------------------------------------- + + rc = nf90_open(trim(input_file), nf90_nowrite, ncid) + do n = 1,nvalid + rc = nf90_inq_varid(ncid, trim(outvars(n)%var_name), varid) + rc = nf90_get_att(ncid, varid, 'long_name', outvars(n)%long_name) + rc = nf90_get_att(ncid, varid, 'units', outvars(n)%units) + rc = nf90_get_att(ncid, varid, '_FillValue', outvars(n)%var_fillvalue) + end do + + ! timestamp + rc = nf90_inq_varid(ncid, 'time', varid) + rc = nf90_get_var(ncid, varid, timestamp) + rc = nf90_get_att(ncid, varid, 'units', timeunit) + rc = nf90_get_att(ncid, varid, 'calendar', timecal) + if (do_ocnpost) then + allocate(z_l(nlevs)) ; z_l = 0.0 + allocate(z_i(0:nlevs)); z_i = 0.0 + allocate(cosrot(nxt*nyt)); cosrot = 0.0 + allocate(sinrot(nxt*nyt)); sinrot = 0.0 + + ! cell centers + rc = nf90_inq_varid(ncid, 'z_l', varid) + rc = nf90_get_var(ncid, varid, z_l) + ! cell edges + rc = nf90_inq_varid(ncid, 'z_i', varid) + rc = nf90_get_var(ncid, varid, z_i) + rc = nf90_close(ncid) + ! rotation angles + call getfield(trim(input_file), trim(cosvar), dims=(/nxt,nyt/), field=cosrot) + call getfield(trim(input_file), trim(sinvar), dims=(/nxt,nyt/), field=sinrot) + else + allocate(anglet(nxt*nyt)); anglet = 0.0 + call getfield(trim(input_file), trim(angvar), dims=(/nxt,nyt/), field=anglet) + cosrot = cos(anglet) + sinrot = -sin(anglet) + end if + + if (debug) then + do n = 1,nvalid + write(logunit,'(a12,i4,a10,3(a6))')trim(outvars(n)%var_name)//', ',outvars(n)%var_dimen, & + ', '//trim(outvars(n)%var_remapmethod),', '//trim(outvars(n)%var_grid), & + ', '//trim(outvars(n)%var_pair),', '//trim(outvars(n)%var_pair_grid) + end do + end if + + ! -------------------------------------------------------- + ! create interpolation masks + ! -------------------------------------------------------- + + if (do_ocnpost) then + allocate(mask3d(nxt*nyt,nlevs)); mask3d = 0.0 + allocate(rgmask3d(nxr*nyr,nlevs)); rgmask3d = 0.0 + else + allocate(mask2d(nxt*nyt)); mask2d = 0.0 + allocate(rgmask2d(nxr*nyr)); rgmask2d = 0.0 + end if + + call remap_masks(vfill) + + ! -------------------------------------------------------- + ! create packed arrays for mapping and remap packed arrays + ! to the destination grid + ! -------------------------------------------------------- + + call setup_packing(nvalid,outvars) + + ! 2D bilin + if (allocated(bilin2d)) then + + wgtsfile = trim(wgtsdir)//'tripole.'//trim(fsrc)//'.Ct.to.rect.'//trim(fdst)//'.bilinear.nc' + call packarrays(trim(input_file), trim(wgtsdir), cosrot, sinrot, b2d, dims=(/nxt,nyt/), & + nflds=nbilin2d, fields=bilin2d) + call remap(trim(wgtsfile), dim2=nbilin2d, src_field=bilin2d, dst_field=rgb2d) + + if (debug) then + write(logunit,'(a)')'remap 2D fields bilinear with '//trim(wgtsfile) + write(logunit,'(a)')'packed min/max values, mapped min/max values' + do n = 1,nbilin2d + write(logunit,'(i4,a10,3(a2,a6),4g14.4)')n,trim(b2d(n)%var_name),' ', & + trim(b2d(n)%var_grid),' ',trim(b2d(n)%var_pair),' ', trim(b2d(n)%var_pair_grid), & + minval(bilin2d(:,n)), maxval(bilin2d(:,n)),minval(rgb2d(:,n)), maxval(rgb2d(:,n)) + end do + call dumpnc(trim(ftype)//'.'//trim(fsrc)//'.bilin2d.nc', 'bilin2d', dims=(/nxt,nyt/), & + nflds=nbilin2d, field=bilin2d) + call dumpnc(trim(ftype)//'.'//trim(fdst)//'.rgbilin2d.nc', 'rgbilin2d', dims=(/nxr,nyr/), & + nflds=nbilin2d, field=rgb2d) + + end if + end if + + ! 2D conserv + if (allocated(consd2d)) then + + wgtsfile = trim(wgtsdir)//'tripole.'//trim(fsrc)//'.Ct.to.rect.'//trim(fdst)//'.conserve.nc' + call packarrays(trim(input_file), trim(wgtsdir), cosrot, sinrot, c2d, dims=(/nxt,nyt/), & + nflds=nconsd2d, fields=consd2d) + call remap(trim(wgtsfile), dim2=nconsd2d, src_field=consd2d, dst_field=rgc2d) + + if (debug) then + write(logunit,'(a)')'remap 2D fields conserv with '//trim(wgtsfile) + write(logunit,'(a)')'packed min/max values, mapped min/max values' + do n = 1,nconsd2d + write(logunit,'(i4,a10,3(a2,a6),4g14.4)')n,trim(c2d(n)%var_name),' ', & + trim(c2d(n)%var_grid),' ',trim(c2d(n)%var_pair),' ', trim(c2d(n)%var_pair_grid), & + minval(consd2d(:,n)), maxval(consd2d(:,n)), minval(rgc2d(:,n)), maxval(rgc2d(:,n)) + end do + call dumpnc(trim(ftype)//'.'//trim(fsrc)//'.consd2d.nc', 'consd2d', dims=(/nxt,nyt/), & + nflds=nconsd2d, field=consd2d) + call dumpnc(trim(ftype)//'.'//trim(fdst)//'.rgconsd2d.nc', 'rgconsd2d', dims=(/nxr,nyr/), & + nflds=nconsd2d, field=rgc2d) + end if + end if + + ! 3D bilin + if (allocated(bilin3d))then + + wgtsfile = trim(wgtsdir)//'tripole.'//trim(fsrc)//'.Ct.to.rect.'//trim(fdst)//'.bilinear.nc' + call packarrays(trim(input_file), trim(wgtsdir), cosrot, sinrot, b3d, dims=(/nxt,nyt,nlevs/), & + nflds=nbilin3d, fields=bilin3d) + call remap(trim(wgtsfile), nk=nlevs, nflds=nbilin3d, src_field=bilin3d, dst_field=rgb3d) + + if (debug) then + write(logunit,'(a)')'remap 3D fields bilinear with '//trim(wgtsfile) + write(logunit,'(a)')'packed min/max values,mapped min/max values' + do n = 1,nbilin3d + write(logunit,'(i4,a10,3(a2,a6),4g14.4)')n,trim(b3d(n)%var_name),' ', & + trim(b3d(n)%var_grid),' ',trim(b3d(n)%var_pair),' ', trim(b3d(n)%var_pair_grid), & + minval(bilin3d(:,:,n)), maxval(bilin3d(:,:,n)),minval(rgb3d(:,:,n)), maxval(rgb3d(:,:,n)) + end do + call dumpnc(trim(ftype)//'.'//trim(fsrc)//'.bilin3d.nc', 'bilin3d', dims=(/nxt,nyt,nlevs/), & + nk=nlevs, nflds=nbilin3d, field=bilin3d) + call dumpnc(trim(ftype)//'.'//trim(fdst)//'.rgbilin3d.nc', 'rgbilin3d', dims=(/nxr,nyr,nlevs/), & + nk=nlevs, nflds=nbilin3d, field=rgb3d) + end if + end if + + ! -------------------------------------------------------- + ! mask the mapped fields + ! -------------------------------------------------------- + + do n = 1,nbilin2d + if (allocated(rgmask3d)) then + where(rgmask3d(:,1) .eq. vfill)rgb2d(:,n) = vfill + end if + if (allocated(rgmask2d))then + where(rgmask2d(:) .eq. vfill)rgb2d(:,n) = vfill + end if + end do + do n = 1,nconsd2d + if (allocated(rgmask3d)) then + where(rgmask3d(:,1) .eq. vfill)rgc2d(:,n) = vfill + end if + if (allocated(rgmask2d))then + where(rgmask2d(:) .eq. vfill)rgc2d(:,n) = vfill + end if + end do + do n = 1,nbilin3d + if (allocated(rgmask3d)) then + where(rgmask3d(:,:) .eq. vfill)rgb3d(:,:,n) = vfill + end if + end do + + ! -------------------------------------------------------- + ! replace model native speed field with a value calculated + ! from remapped ssu,ssv + ! -------------------------------------------------------- + + if (do_ocnpost) then + do n = 1,nbilin2d + if (trim(b2d(n)%var_name) == 'speed')idx1 = n + if (trim(b2d(n)%var_name) == 'SSU')idx2 = n + if (trim(b2d(n)%var_name) == 'SSV')idx3 = n + enddo + where(rgb2d(:,idx1) .ne. vfill)rgb2d(:,idx1) = & + sqrt(rgb2d(:,idx2)**2 + rgb2d(:,idx3)**2) + end if + + ! -------------------------------------------------------- + ! write the mapped fields + ! -------------------------------------------------------- + + allocate(out2d(nxr,nyr)); out2d = 0.0 + allocate(out3d(nxr,nyr,nlevs)); out3d = 0.0 + + fout = trim(ftype)//'.'//trim(fdst)//'.nc' + if (debug) write(logunit, '(a)')'output file: '//trim(fout) + + rc = nf90_create(trim(fout), nf90_clobber, ncid) + rc = nf90_def_dim(ncid, 'longitude', nxr, idimid) + rc = nf90_def_dim(ncid, 'latitude', nyr, jdimid) + rc = nf90_def_dim(ncid, 'time', nf90_unlimited, timid) + + ! define the time variable + rc = nf90_def_var(ncid, 'time', nf90_double, (/timid/), varid) + rc = nf90_put_att(ncid, varid, 'units', trim(timeunit)) + rc= nf90_put_att(ncid, varid, 'calendar', trim(timecal)) + ! spatial grid + rc = nf90_def_var(ncid, 'longitude', nf90_float, (/idimid/), varid) + rc = nf90_put_att(ncid, varid, 'units', 'degrees_east') + rc = nf90_def_var(ncid, 'latitude', nf90_float, (/jdimid/), varid) + rc = nf90_put_att(ncid, varid, 'units', 'degrees_north') + ! vertical grid + if (do_ocnpost) then + rc = nf90_def_dim(ncid, 'z_l', nlevs , kdimid) + rc = nf90_def_dim(ncid, 'z_i', nlevs+1, edimid) + rc = nf90_def_var(ncid, 'z_l', nf90_float, (/kdimid/), varid) + rc = nf90_put_att(ncid, varid, 'units', 'm') + rc = nf90_put_att(ncid, varid, 'positive', 'down') + rc = nf90_def_var(ncid, 'z_i', nf90_float, (/edimid/), varid) + rc = nf90_put_att(ncid, varid, 'units', 'm') + rc = nf90_put_att(ncid, varid, 'positive', 'down') + end if + + if (allocated(b2d)) then + do n = 1,nbilin2d + vname = trim(b2d(n)%var_name) + vunit = trim(b2d(n)%units) + vlong = trim(b2d(n)%long_name) + vfill = b2d(n)%var_fillvalue + rc = nf90_def_var(ncid, vname, nf90_float, (/idimid,jdimid,timid/), varid) + rc = nf90_put_att(ncid, varid, 'units', vunit) + rc = nf90_put_att(ncid, varid, 'long_name', vlong) + rc = nf90_put_att(ncid, varid, '_FillValue', vfill) + enddo + end if + if (allocated(c2d)) then + do n = 1,nconsd2d + vname = trim(c2d(n)%var_name) + vunit = trim(c2d(n)%units) + vlong = trim(c2d(n)%long_name) + vfill = c2d(n)%var_fillvalue + rc = nf90_def_var(ncid, vname, nf90_float, (/idimid,jdimid,timid/), varid) + rc = nf90_put_att(ncid, varid, 'units', vunit) + rc = nf90_put_att(ncid, varid, 'long_name', vlong) + rc = nf90_put_att(ncid, varid, '_FillValue', vfill) + enddo + end if + if (allocated(b3d)) then + do n = 1,nbilin3d + vname = trim(b3d(n)%var_name) + vunit = trim(b3d(n)%units) + vlong = trim(b3d(n)%long_name) + vfill = b3d(n)%var_fillvalue + rc = nf90_def_var(ncid, vname, nf90_float, (/idimid,jdimid,kdimid,timid/), varid) + rc = nf90_put_att(ncid, varid, 'units', vunit) + rc = nf90_put_att(ncid, varid, 'long_name', vlong) + rc = nf90_put_att(ncid, varid, '_FillValue', vfill) + enddo + end if + rc = nf90_enddef(ncid) + + ! dimensions + rc = nf90_inq_varid(ncid, 'longitude', varid) + rc = nf90_put_var(ncid, varid, dstlon(:,1)) + rc = nf90_inq_varid(ncid, 'latitude', varid) + rc = nf90_put_var(ncid, varid, dstlat(1,:)) + ! time + rc = nf90_inq_varid(ncid, 'time', varid) + rc = nf90_put_var(ncid, varid, timestamp) + ! vertical + if (do_ocnpost) then + rc = nf90_inq_varid(ncid, 'z_l', varid) + rc = nf90_put_var(ncid, varid, z_l) + rc = nf90_inq_varid(ncid, 'z_i', varid) + rc = nf90_put_var(ncid, varid, z_i) + end if + if (allocated(rgb2d)) then + do n = 1,nbilin2d + out2d(:,:) = reshape(rgb2d(:,n), (/nxr,nyr/)) + out2d(:,nyr) = vfill + vname = trim(b2d(n)%var_name) + rc = nf90_inq_varid(ncid, vname, varid) + rc = nf90_put_var(ncid, varid, out2d) + end do + end if + if (allocated(rgc2d)) then + do n = 1,nconsd2d + out2d(:,:) = reshape(rgc2d(:,n), (/nxr,nyr/)) + out2d(:,nyr) = vfill + vname = trim(c2d(n)%var_name) + rc = nf90_inq_varid(ncid, vname, varid) + rc = nf90_put_var(ncid, varid, out2d) + end do + end if + if (allocated(rgb3d)) then + do n = 1,nbilin3d + out3d(:,:,:) = reshape(rgb3d(:,:,n), (/nxr,nyr,nlevs/)) + out3d(:,nyr,:) = vfill + vname = trim(b3d(n)%var_name) + rc = nf90_inq_varid(ncid, vname, varid) + rc = nf90_put_var(ncid, varid, out3d) + end do + end if + rc = nf90_close(ncid) + write(logunit,'(a)')trim(fout)//' done' + +end program ocnicepost diff --git a/src/ocnicepost.fd/utils_mod.F90 b/src/ocnicepost.fd/utils_mod.F90 new file mode 100644 index 00000000..82183c5e --- /dev/null +++ b/src/ocnicepost.fd/utils_mod.F90 @@ -0,0 +1,601 @@ +module utils_mod + + use netcdf + use init_mod, only : debug, logunit, vardefs + + implicit none + + private + + interface getfield + module procedure getfield2d + module procedure getfield3d + end interface getfield + + interface packarrays + module procedure packarrays2d + module procedure packarrays3d + end interface packarrays + + interface remap + module procedure remap1d + module procedure remap2d + module procedure remap3d + end interface remap + + interface getvecpair + module procedure getvecpair2d + module procedure getvecpair3d + end interface getvecpair + + interface dumpnc + module procedure dumpnc1d + module procedure dumpnc2d + module procedure dumpnc3d + module procedure dumpnc3dk + end interface dumpnc + + public getfield + public packarrays + public remap + public dumpnc + +contains + + !---------------------------------------------------------- + ! pack 2D fields into arrays by mapping type + !---------------------------------------------------------- + subroutine packarrays2d(filesrc, wgtsdir, cosrot, sinrot, vars, dims, nflds, fields) + + character(len=*), intent(in) :: filesrc,wgtsdir + real, intent(in) :: cosrot(:),sinrot(:) + type(vardefs), intent(in) :: vars(:) + integer, intent(in) :: dims(:) + integer, intent(in) :: nflds + real, intent(out) :: fields(:,:) + + ! local variables + integer :: n, nn + real, allocatable, dimension(:,:) :: vecpair + character(len=20) :: subname = 'packarrays2d' + + fields=0.0 + + if (debug)write(logunit,'(a)')'enter '//trim(subname) + ! obtain vector pairs + do n = 1,nflds + if (trim(vars(n)%var_grid) == 'Cu' .or. trim(vars(n)%var_grid) == 'Bu_x') then + allocate(vecpair(dims(1)*dims(2),2)); vecpair = 0.0 + call getvecpair(trim(filesrc), trim(wgtsdir), cosrot, sinrot, & + trim(vars(n)%var_name), trim(vars(n)%var_grid(1:2)), & + trim(vars(n)%var_pair), trim(vars(n)%var_pair_grid(1:2)), & + dims=(/dims(1),dims(2)/), vecpair=vecpair) + end if + end do + + ! create packed array + nn = 0 + do n = 1,nflds + if (len_trim(vars(n)%var_pair) == 0) then + nn = nn + 1 + call getfield(trim(filesrc), trim(vars(n)%var_name), dims=(/dims(1),dims(2)/), & + field=fields(:,nn)) + else ! fill with vector pairs + nn = nn+1 + ! ocn vectors + if (trim(vars(n)%var_grid) == 'Cu')fields(:,nn) = vecpair(:,1) + if (trim(vars(n)%var_grid) == 'Cv')fields(:,nn) = vecpair(:,2) + ! ice vectors + if (trim(vars(n)%var_grid) == 'Bu_x')fields(:,nn) = vecpair(:,1) + if (trim(vars(n)%var_grid) == 'Bu_y')fields(:,nn) = vecpair(:,2) + end if + end do + + if (debug)write(logunit,'(a)')'exit '//trim(subname) + end subroutine packarrays2d + + !---------------------------------------------------------- + ! pack 3D fields into arrays by mapping type + !---------------------------------------------------------- + subroutine packarrays3d(filesrc, wgtsdir, cosrot, sinrot, vars, dims, nflds, fields) + + character(len=*), intent(in) :: filesrc,wgtsdir + real, intent(in) :: cosrot(:),sinrot(:) + type(vardefs), intent(in) :: vars(:) + integer, intent(in) :: dims(:) + integer, intent(in) :: nflds + real, intent(out) :: fields(:,:,:) + + ! local variables + integer :: n, nn + real, allocatable, dimension(:,:,:) :: vecpair + character(len=20) :: subname = 'packarrays3d' + + fields=0.0 + + if (debug)write(logunit,'(a)')'enter '//trim(subname) + ! obtain vector pairs + do n = 1,dims(3) + if (trim(vars(n)%var_grid) == 'Cu') then + allocate(vecpair(dims(1)*dims(2),dims(3),2)); vecpair = 0.0 + call getvecpair(trim(filesrc), trim(wgtsdir), cosrot, sinrot, & + trim(vars(n)%var_name), trim(vars(n)%var_grid), & + trim(vars(n)%var_pair), trim(vars(n)%var_pair_grid), & + dims=(/dims(1),dims(2),dims(3)/), vecpair=vecpair) + end if + end do + + ! create packed array + nn = 0 + do n = 1,nflds + if (len_trim(vars(n)%var_pair) == 0) then + nn = nn + 1 + call getfield(trim(filesrc), trim(vars(n)%var_name), dims=(/dims(1),dims(2),dims(3)/), & + field=fields(:,:,nn)) + else ! fill with vector pairs + nn = nn+1 + if (trim(vars(n)%var_grid) == 'Cu')fields(:,:,nn) = vecpair(:,:,1) + if (trim(vars(n)%var_grid) == 'Cv')fields(:,:,nn) = vecpair(:,:,2) + end if + end do + + if (debug)write(logunit,'(a)')'exit '//trim(subname) + end subroutine packarrays3d + + !---------------------------------------------------------- + ! obtain 2D vector pairs mapped to Ct and rotated to EW + !---------------------------------------------------------- + subroutine getvecpair2d(fname, wdir, cosrot, sinrot, vname1, vgrid1, & + vname2, vgrid2, dims, vecpair) + + character(len=*), intent(in) :: fname + character(len=*), intent(in) :: wdir + real, intent(in) :: cosrot(:), sinrot(:) + character(len=*), intent(in) :: vname1, vgrid1, vname2, vgrid2 + integer, intent(in) :: dims(:) + real, intent(out) :: vecpair(:,:) + + ! local variables + integer :: ii + real, dimension(dims(1)*dims(2)) :: urot, vrot + character(len=240) :: wgtsfile + character(len=20) :: subname = 'getvecpair2d' + + if (debug)write(logunit,'(a)')'enter '//trim(subname) + + wgtsfile = trim(wdir)//'tripole.mx025.'//vgrid1//'.to.Ct.bilinear.nc' + call getfield(fname, vname1, dims=dims, field=vecpair(:,1), wgts=trim(wgtsfile)) + if (debug)write(logunit,'(a)')'wgtsfile for 2d vector '//trim(vname1)//' '//trim(wgtsfile) + wgtsfile = trim(wdir)//'tripole.mx025.'//vgrid2//'.to.Ct.bilinear.nc' + call getfield(fname, vname2, dims=dims, field=vecpair(:,2), wgts=trim(wgtsfile)) + if (debug)write(logunit,'(a)')'wgtsfile for 2d vector '//trim(vname2)//' '//trim(wgtsfile) + + urot = 0.0; vrot = 0.0 + do ii = 1,dims(1)*dims(2) + urot(ii) = vecpair(ii,1)*cosrot(ii) + vecpair(ii,2)*sinrot(ii) + vrot(ii) = vecpair(ii,2)*cosrot(ii) - vecpair(ii,1)*sinrot(ii) + end do + vecpair(:,1) = urot(:) + vecpair(:,2) = vrot(:) + + if (debug) write(logunit,'(a)')'exit '//trim(subname) + end subroutine getvecpair2d + + !---------------------------------------------------------- + ! obtain 3D vector pairs, mapped to Ct and rotated to EW + !---------------------------------------------------------- + subroutine getvecpair3d(fname, wdir, cosrot, sinrot, vname1, vgrid1, & + vname2, vgrid2, dims, vecpair) + + character(len=*), intent(in) :: fname + character(len=*), intent(in) :: wdir + real, intent(in) :: cosrot(:), sinrot(:) + character(len=*), intent(in) :: vname1, vgrid1, vname2, vgrid2 + integer, intent(in) :: dims(:) + real, intent(out) :: vecpair(:,:,:) + + ! local variables + integer :: ii,k + real, dimension(dims(1)*dims(2)) :: urot, vrot + character(len=240) :: wgtsfile + character(len=20) :: subname = 'getvecpair3d' + + if (debug)write(logunit,'(a)')'enter '//trim(subname) + + wgtsfile = trim(wdir)//'tripole.mx025.'//vgrid1//'.to.Ct.bilinear.nc' + call getfield(fname, vname1, dims=dims, field=vecpair(:,:,1), wgts=trim(wgtsfile)) + wgtsfile = trim(wdir)//'tripole.mx025.'//vgrid2//'.to.Ct.bilinear.nc' + call getfield(fname, vname2, dims=dims, field=vecpair(:,:,2), wgts=trim(wgtsfile)) + + do k = 1,dims(3) + urot = 0.0; vrot = 0.0 + do ii = 1,dims(1)*dims(2) + urot(ii)= vecpair(ii,k,1)*cosrot(ii) + vecpair(ii,k,2)*sinrot(ii) + vrot(ii)= vecpair(ii,k,2)*cosrot(ii) - vecpair(ii,k,1)*sinrot(ii) + end do + vecpair(:,k,1) = urot(:) + vecpair(:,k,2) = vrot(:) + end do + + if (debug) write(logunit,'(a)')'exit '//trim(subname) + end subroutine getvecpair3d + + !---------------------------------------------------------- + ! obtain a 2D field and return a 1-D vector array + !---------------------------------------------------------- + subroutine getfield2d(fname, vname, dims, field, wgts) + + character(len=*), intent(in) :: fname, vname + integer, intent(in) :: dims(:) + real, intent(out) :: field(:) + character(len=*), optional, intent(in) :: wgts + + ! local variable + integer :: ncid, varid, rc + real :: fval + real, allocatable :: a2d(:,:) + real, allocatable :: atmp(:) + character(len=20) :: subname = 'getfield2d' + + if (debug)write(logunit,'(a)')'enter '//trim(subname)//' variable '//vname + + allocate(a2d(dims(1),dims(2))); a2d = 0.0 + allocate(atmp(dims(1)*dims(2))); atmp = 0.0 + + rc = nf90_open(fname, nf90_nowrite, ncid) + call handle_err(rc,' nf90_open '//fname) + rc = nf90_inq_varid(ncid, vname, varid) + call handle_err(rc,' get variable ID '// vname) + rc = nf90_get_var(ncid, varid, a2d) + call handle_err(rc,' get variable'// vname) + rc = nf90_get_att(ncid, varid, '_FillValue', fval) + rc = nf90_close(ncid) + + atmp(:) = reshape(a2d, (/dims(1)*dims(2)/)) + where(atmp .eq. fval)atmp = 0.0 + if(present(wgts)) then + call remap(trim(wgts), src_field=atmp, dst_field=field) + else + field = atmp + end if + + if (debug) write(logunit,'(a)')'exit '//trim(subname)//' variable '//vname + end subroutine getfield2d + + !---------------------------------------------------------- + ! obtain a 3D field and return a 2-D vector array + !---------------------------------------------------------- + subroutine getfield3d(fname, vname, dims, field, wgts) + + character(len=*), intent(in) :: fname, vname + integer, intent(in) :: dims(:) + real, intent(out) :: field(:,:) + character(len=*), optional, intent(in) :: wgts + + ! local variable + integer :: ncid, varid, rc + real :: fval + real, allocatable :: a3d(:,:,:) + real, allocatable :: atmp(:,:) + character(len=20) :: subname = 'getfield3d' + + if (debug)write(logunit,'(a)')'enter '//trim(subname)//' variable '//vname + + allocate(a3d(dims(1),dims(2),dims(3))); a3d = 0.0 + allocate(atmp(dims(1)*dims(2),dims(3))); atmp = 0.0 + + rc = nf90_open(fname, nf90_nowrite, ncid) + call handle_err(rc,' nf90_open '//fname) + rc = nf90_inq_varid(ncid, vname, varid) + call handle_err(rc,' get variable ID '// vname) + rc = nf90_get_var(ncid, varid, a3d) + call handle_err(rc,' get variable'// vname) + rc = nf90_get_att(ncid, varid, '_FillValue', fval) + rc = nf90_close(ncid) + + atmp(:,:) = reshape(a3d, (/dims(1)*dims(2),dims(3)/)) + where(atmp .eq. fval)atmp = 0.0 + if(present(wgts)) then + call remap(trim(wgts), dim2=dims(3), src_field=atmp, dst_field=field) + else + field = atmp + end if + + if (debug) write(logunit,'(a)')'exit '//trim(subname)//' variable '//vname + end subroutine getfield3d + + !---------------------------------------------------------- + ! remap a 1-D vector array + !---------------------------------------------------------- + subroutine remap1d(fname, src_field, dst_field) + + character(len=*), intent(in) :: fname + real, intent(in) :: src_field(:) + real, intent(out) :: dst_field(:) + + ! local variables + integer :: ncid, rc, id + integer :: i,ii,jj + integer :: n_a, n_b, n_s + integer(kind=4), allocatable, dimension(:) :: col, row + real(kind=8), allocatable, dimension(:) :: S + character(len=20) :: subname = 'remap1d' + + if (debug)write(logunit,'(a)')'enter '//trim(subname) + + ! retrieve the weights + rc = nf90_open(trim(fname), nf90_nowrite, ncid) + call handle_err(rc,' nf90_open '//fname) + rc = nf90_inq_dimid(ncid, 'n_s', id) + rc = nf90_inquire_dimension(ncid, id, len=n_s) + rc = nf90_inq_dimid(ncid, 'n_a', id) + rc = nf90_inquire_dimension(ncid, id, len=n_a) + rc = nf90_inq_dimid(ncid, 'n_b', id) + rc = nf90_inquire_dimension(ncid, id, len=n_b) + + allocate(col(1:n_s)); col = 0 + allocate(row(1:n_s)); row = 0 + allocate( S(1:n_s)); S = 0.0 + + rc = nf90_inq_varid(ncid, 'col', id) + rc = nf90_get_var(ncid, id, col) + rc = nf90_inq_varid(ncid, 'row', id) + rc = nf90_get_var(ncid, id, row) + rc = nf90_inq_varid(ncid, 'S', id) + rc = nf90_get_var(ncid, id, S) + rc = nf90_close(ncid) + + dst_field = 0.0 + do i = 1,n_s + ii = row(i); jj = col(i) + dst_field(ii) = dst_field(ii) + S(i)*real(src_field(jj),8) + enddo + + if (debug) write(logunit,'(a)')'exit '//trim(subname) + end subroutine remap1d + + !---------------------------------------------------------- + ! remap a packed field of either nflds or nlevs + !---------------------------------------------------------- + subroutine remap2d(fname, dim2, src_field, dst_field) + + character(len=*), intent(in) :: fname + integer, intent(in) :: dim2 + real, intent(in) :: src_field(:,:) + real, intent(out) :: dst_field(:,:) + + ! local variables + integer :: ncid, rc, id + integer :: i,ii,jj + integer :: n_a, n_b, n_s + integer(kind=4), allocatable, dimension(:) :: col, row + real(kind=8), allocatable, dimension(:) :: S + character(len=20) :: subname = 'remap2d' + + if (debug)write(logunit,'(a)')'enter '//trim(subname)//' weights = '//trim(fname) + + ! retrieve the weights + rc = nf90_open(trim(fname), nf90_nowrite, ncid) + call handle_err(rc,' nf90_open '//fname) + rc = nf90_inq_dimid(ncid, 'n_s', id) + rc = nf90_inquire_dimension(ncid, id, len=n_s) + rc = nf90_inq_dimid(ncid, 'n_a', id) + rc = nf90_inquire_dimension(ncid, id, len=n_a) + rc = nf90_inq_dimid(ncid, 'n_b', id) + rc = nf90_inquire_dimension(ncid, id, len=n_b) + + allocate(col(1:n_s)); col = 0 + allocate(row(1:n_s)); row = 0 + allocate( S(1:n_s)); S = 0.0 + + rc = nf90_inq_varid(ncid, 'col', id) + rc = nf90_get_var(ncid, id, col) + rc = nf90_inq_varid(ncid, 'row', id) + rc = nf90_get_var(ncid, id, row) + rc = nf90_inq_varid(ncid, 'S', id) + rc = nf90_get_var(ncid, id, S) + rc = nf90_close(ncid) + + dst_field = 0.0 + do i = 1,n_s + ii = row(i); jj = col(i) + dst_field(ii,:) = dst_field(ii,:) + S(i)*real(src_field(jj,:),8) + enddo + + if (debug) write(logunit,'(a)')'exit '//trim(subname) + end subroutine remap2d + + !---------------------------------------------------------- + ! remap a field packed array of nk levels and nflds fields + !---------------------------------------------------------- + subroutine remap3d(fname, nk, nflds, src_field, dst_field) + + character(len=*), intent(in) :: fname + integer, intent(in) :: nk, nflds + real, intent(in) :: src_field(:,:,:) + real, intent(out) :: dst_field(:,:,:) + + ! local variables + integer :: ncid, rc, id + integer :: i,ii,jj + integer :: n_a, n_b, n_s + integer(kind=4), allocatable, dimension(:) :: col, row + real(kind=8), allocatable, dimension(:) :: S + character(len=20) :: subname = 'remap3d' + + if (debug)write(logunit,'(a)')'enter '//trim(subname)//' weights = '//trim(fname) + + ! retrieve the weights + rc = nf90_open(trim(fname), nf90_nowrite, ncid) + call handle_err(rc,' nf90_open '//fname) + rc = nf90_inq_dimid(ncid, 'n_s', id) + rc = nf90_inquire_dimension(ncid, id, len=n_s) + rc = nf90_inq_dimid(ncid, 'n_a', id) + rc = nf90_inquire_dimension(ncid, id, len=n_a) + rc = nf90_inq_dimid(ncid, 'n_b', id) + rc = nf90_inquire_dimension(ncid, id, len=n_b) + + allocate(col(1:n_s)); col = 0 + allocate(row(1:n_s)); row = 0 + allocate( S(1:n_s)); S = 0.0 + + rc = nf90_inq_varid(ncid, 'col', id) + rc = nf90_get_var(ncid, id, col) + rc = nf90_inq_varid(ncid, 'row', id) + rc = nf90_get_var(ncid, id, row) + rc = nf90_inq_varid(ncid, 'S', id) + rc = nf90_get_var(ncid, id, S) + rc = nf90_close(ncid) + + dst_field = 0.0 + do i = 1,n_s + ii = row(i); jj = col(i) + dst_field(ii,:,:) = dst_field(ii,:,:) + S(i)*real(src_field(jj,:,:),8) + enddo + + if (debug) write(logunit,'(a)')'exit '//trim(subname) + end subroutine remap3d + + !---------------------------------------------------------- + ! write a bare netcdf file of a 2D packed field + !---------------------------------------------------------- + subroutine dumpnc2d(fname, vname, dims, nflds, field) + + character(len=*), intent(in) :: fname, vname + integer, intent(in) :: dims(:) + integer, intent(in) :: nflds + real, intent(in) :: field(:,:) + + ! local variable + integer :: ncid, varid, rc, idimid, jdimid, fdimid + real, allocatable :: a3d(:,:,:) + character(len=20) :: subname = 'dumpnc2d' + + if (debug)write(logunit,'(a)')'enter '//trim(subname)//' variable '//vname + allocate(a3d(dims(1),dims(2),nflds)); a3d = 0.0 + + rc = nf90_create(trim(fname), nf90_clobber, ncid) + rc = nf90_def_dim(ncid, 'nx', dims(1), idimid) + rc = nf90_def_dim(ncid, 'ny', dims(2), jdimid) + rc = nf90_def_dim(ncid, 'nf', nflds, fdimid) + rc = nf90_def_var(ncid, vname, nf90_float, (/idimid,jdimid,fdimid/), varid) + rc = nf90_enddef(ncid) + + a3d(:,:,:) = reshape(field(1:dims(1)*dims(2),1:nflds), (/dims(1),dims(2),nflds/)) + rc = nf90_put_var(ncid, varid, a3d) + rc = nf90_close(ncid) + + if (debug)write(logunit,'(a)')'exit '//trim(subname)//' variable '//vname + end subroutine dumpnc2d + + !---------------------------------------------------------- + ! write a bare netcdf file of a packed 3D field + !---------------------------------------------------------- + subroutine dumpnc3d(fname, vname, dims, nk, nflds, field) + + character(len=*), intent(in) :: fname, vname + integer, intent(in) :: dims(:) + integer, intent(in) :: nk, nflds + real, intent(in) :: field(:,:,:) + + ! local variable + integer :: n, ncid, varid, rc, idimid, jdimid, kdimid, fdimid + real, allocatable :: a4d(:,:,:,:) + character(len=20) :: subname = 'dumpnc3d' + + if (debug)write(logunit,'(a)')'enter '//trim(subname)//' variable '//vname + allocate(a4d(dims(1),dims(2),dims(3),nflds)); a4d = 0.0 + + rc = nf90_create(trim(fname), nf90_clobber, ncid) + rc = nf90_def_dim(ncid, 'nx', dims(1), idimid) + rc = nf90_def_dim(ncid, 'ny', dims(2), jdimid) + rc = nf90_def_dim(ncid, 'nk', dims(3), kdimid) + rc = nf90_def_dim(ncid, 'nf', nflds, fdimid) + rc = nf90_def_var(ncid, vname, nf90_float, (/idimid,jdimid,kdimid,fdimid/), varid) + rc = nf90_enddef(ncid) + + do n = 1,nflds + a4d(:,:,:,n) = reshape(field(1:dims(1)*dims(2),1:dims(3),n), (/dims(1),dims(2),dims(3)/)) + end do + rc = nf90_put_var(ncid, varid, a4d) + rc = nf90_close(ncid) + + if (debug)write(logunit,'(a)')'exit '//trim(subname)//' variable '//vname + end subroutine dumpnc3d + + !---------------------------------------------------------- + ! write a bare netcdf file of an unpacked 3D field + !---------------------------------------------------------- + subroutine dumpnc3dk(fname, vname, dims, field) + + character(len=*), intent(in) :: fname, vname + integer, intent(in) :: dims(:) + real, intent(in) :: field(:,:) + + ! local variable + integer :: ncid, varid, rc, idimid, jdimid, kdimid + real, allocatable :: a3d(:,:,:) + character(len=20) :: subname = 'dumpnc3dk' + + if (debug)write(logunit,'(a)')'enter '//trim(subname)//' variable '//vname + allocate(a3d(dims(1),dims(2),dims(3))); a3d = 0.0 + + rc = nf90_create(trim(fname), nf90_clobber, ncid) + rc = nf90_def_dim(ncid, 'nx', dims(1), idimid) + rc = nf90_def_dim(ncid, 'ny', dims(2), jdimid) + rc = nf90_def_dim(ncid, 'nk', dims(3), kdimid) + rc = nf90_def_var(ncid, vname, nf90_float, (/idimid,jdimid,kdimid/), varid) + rc = nf90_enddef(ncid) + + a3d(:,:,:) = reshape(field(1:dims(1)*dims(2),1:dims(3)), (/dims(1),dims(2),dims(3)/)) + rc = nf90_put_var(ncid, varid, a3d) + rc = nf90_close(ncid) + + if (debug)write(logunit,'(a)')'exit '//trim(subname)//' variable '//vname + + end subroutine dumpnc3dk + + !---------------------------------------------------------- + ! write a bare netcdf file of an unpacked 2D field + !---------------------------------------------------------- + subroutine dumpnc1d(fname, vname, dims, field) + + character(len=*), intent(in) :: fname, vname + integer, intent(in) :: dims(:) + real, intent(in) :: field(:) + + ! local variable + integer :: ncid, varid, rc, idimid, jdimid + real, allocatable :: a2d(:,:) + character(len=20) :: subname = 'dumpnc1d' + + if (debug)write(logunit,'(a)')'enter '//trim(subname)//' variable '//vname + allocate(a2d(dims(1),dims(2))); a2d = 0.0 + + rc = nf90_create(trim(fname), nf90_clobber, ncid) + rc = nf90_def_dim(ncid, 'nx', dims(1), idimid) + rc = nf90_def_dim(ncid, 'ny', dims(2), jdimid) + rc = nf90_def_var(ncid, vname, nf90_float, (/idimid,jdimid/), varid) + rc = nf90_enddef(ncid) + + a2d(:,:) = reshape(field(1:dims(1)*dims(2)), (/dims(1),dims(2)/)) + rc = nf90_put_var(ncid, varid, a2d) + rc = nf90_close(ncid) + + if (debug)write(logunit,'(a)')'exit '//trim(subname)//' variable '//vname + + end subroutine dumpnc1d + + !---------------------------------------------------------- + ! handle netcdf errors + !---------------------------------------------------------- + subroutine handle_err(ierr,string) + + integer , intent(in) :: ierr + character(len=*), intent(in) :: string + if (ierr /= nf90_noerr) then + write(logunit,'(a)') '*** ERROR ***: '//trim(string)//':'//trim(nf90_strerror(ierr)) + stop + end if + end subroutine handle_err +end module utils_mod