From a9170c241363a6c61da51d6f0aab377b56dcac18 Mon Sep 17 00:00:00 2001 From: "denise.worthen" Date: Sat, 14 Oct 2023 13:08:30 -0400 Subject: [PATCH 1/4] add fortran code for ice or ocn post --- src/ocnicepost.fd/ocnicepost.F90 | 577 +++++++++++++++++++++++++++++ src/ocnicepost.fd/ocnicepost.nml | 5 + src/ocnicepost.fd/outputvars.F90 | 146 ++++++++ src/ocnicepost.fd/utils_mod.F90 | 604 +++++++++++++++++++++++++++++++ 4 files changed, 1332 insertions(+) create mode 100644 src/ocnicepost.fd/ocnicepost.F90 create mode 100644 src/ocnicepost.fd/ocnicepost.nml create mode 100644 src/ocnicepost.fd/outputvars.F90 create mode 100644 src/ocnicepost.fd/utils_mod.F90 diff --git a/src/ocnicepost.fd/ocnicepost.F90 b/src/ocnicepost.fd/ocnicepost.F90 new file mode 100644 index 00000000..ca21bf3a --- /dev/null +++ b/src/ocnicepost.fd/ocnicepost.F90 @@ -0,0 +1,577 @@ +program ocnicepost + + use utils_mod, only : debug, logunit, getfield, packarrays, remap, dumpnc + use outputvars + use netcdf + + implicit none + + character(len= 20) :: fnml, ftype + character(len=120) :: wgtsdir + character(len=120) :: input_file, wgtsfile, output_file + + ! source grid, tripole 1/4 deg, 40 vertical levels + integer, parameter :: nxt = 1440, nyt = 1080, nlevs = 40 + + ! destination grids + integer, parameter :: ndest = 3 + integer, parameter, dimension(ndest) :: nxrs = (/1440, 720, 360/) + integer, parameter, dimension(ndest) :: nyrs = (/ 721, 361, 181/) + character(len=4), dimension(ndest) :: dstgrds = (/'0p25', '0p5 ', '1p0 '/) + + ! packed arrays on source grid + real, allocatable, dimension(:,:) :: bilin2d, consd2d + real, allocatable, dimension(:,:,:) :: bilin3d + ! variable types + type(vardefs), allocatable :: b2d(:) + type(vardefs), allocatable :: c2d(:) + type(vardefs), allocatable :: b3d(:) + + ! source grid fields + real, dimension(nxt,nyt,nlevs) :: tmp3d + real, dimension(nxt*nyt,nlevs) :: mask3d + real, dimension(nxt*nyt) :: mask2d + real, dimension(nxt*nyt) :: cosrot, sinrot, anglet + real, dimension(nlevs) :: z_l + real, dimension(0:nlevs) :: z_i + + ! destination grid fields + real, allocatable, dimension(:,:) :: dstlon, dstlat + real, allocatable, dimension(:,:) :: rgmask3d + real, allocatable, dimension(:) :: rgmask2d + ! output fields (x,y,z) + real, allocatable, dimension(:) :: out1d + real, allocatable, dimension(:,:) :: out2d + real, allocatable, dimension(:,:,:) :: out3d + ! packed remapped fields + real, allocatable, dimension(:,:) :: rgb2d, rgc2d + real, allocatable, dimension(:,:,:) :: rgb3d + + real(kind=8) :: timestamp + character(len= 40) :: timeunit, timecal + character(len= 20) :: vname, vunit + character(len=120) :: vlong + character(len=4) :: dstgrid + + real :: vfill + integer :: nd, nxr, nyr + integer :: i,j,k,n,nn,nvalid,iounit + integer :: rc,ncid,varid,dimid + integer :: nbilin2d,nbilin3d,nconsd2d + integer :: idimid,jdimid,kdimid,edimid,timid + integer :: idx1,idx2,idx3 + logical :: do_ocnpost + + namelist /ocnicepost_nml/ ftype, wgtsdir, debug + + ! -------------------------------------------------------- + ! read the name list + ! -------------------------------------------------------- + + fnml = 'ocnicepost.nml' + inquire (file=trim(fnml), iostat=rc) + if (rc /= 0) then + write (6, '(3a)') 'Error: input file "', trim(fnml), '" does not exist.' + stop + end if + + ! Open and read Namelist file. + open (action='read', file=trim(fnml), iostat=rc, newunit=iounit) + read (nml=ocnicepost_nml, iostat=rc, unit=iounit) + if (rc /= 0) then + write (6, '(a)') 'Error: invalid Namelist format.' + end if + close (iounit) + + ! initialize the source file type and variables + if (trim(ftype) == 'ocean') then + do_ocnpost = .true. + call ocnvars_typedefine + else + do_ocnpost = .false. + call icevars_typedefine + 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) + + ! -------------------------------------------------------- + ! read the source file and obtain the units and long name, + ! rotation angles, vertical grid and time axis + ! -------------------------------------------------------- + + nvalid = 0 + rc = nf90_open(trim(input_file), nf90_nowrite, ncid) + do i = 1,maxvars + if (len_trim(outvars(i)%input_var_name) > 0 ) then + rc = nf90_inq_varid(ncid, trim(outvars(i)%input_var_name), varid) + rc = nf90_get_att(ncid, varid, 'long_name', outvars(i)%long_name) + rc = nf90_get_att(ncid, varid, 'units', outvars(i)%units) + rc = nf90_get_att(ncid, varid, '_FillValue', outvars(i)%var_fillvalue) + nvalid = nvalid+1 + if (trim(outvars(i)%input_var_name) == 'temp')vfill = outvars(i)%var_fillvalue + if (trim(outvars(i)%input_var_name) == 'aice_h')vfill = outvars(i)%var_fillvalue + end if + 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 + ! 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), 'cos_rot', dims=(/nxt,nyt/), field=cosrot) + call getfield(trim(input_file), 'sin_rot', dims=(/nxt,nyt/), field=sinrot) + else + call getfield(trim(input_file), 'ANGLET', dims=(/nxt,nyt/), field=anglet) + cosrot = cos(anglet) + sinrot = -sin(anglet) + end if + + ! -------------------------------------------------------- + ! 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 + rc = nf90_open(trim(input_file), nf90_nowrite, ncid) + ! 3D temp to use as mask, obtain directly from file to preserve vfill + rc = nf90_inq_varid(ncid, 'temp', varid) + 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 + else + call getfield(trim(input_file), 'tmask', dims=(/nxt,nyt/), field=mask2d) + ! 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 if + + ! -------------------------------------------------------- + ! 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(outvars(n)%var_remapmethod) == 'bilinear') then + if (outvars(n)%var_dimen == 2) nbilin2d = nbilin2d + 1 + if (outvars(n)%var_dimen == 3) nbilin3d = nbilin3d + 1 + end if + if (trim(outvars(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 + ! -------------------------------------------------------- + + i = 0; j = 0; k = 0 + do n = 1,nvalid + if (trim(outvars(n)%var_remapmethod) == 'bilinear') then + if (outvars(n)%var_dimen == 2 .and. allocated(b2d)) then + i = i+1; b2d(i) = outvars(n) + end if + if (outvars(n)%var_dimen == 3 .and. allocated(b3d)) then + j = j+1; b3d(j) = outvars(n) + end if + end if + if (trim(outvars(n)%var_remapmethod) == 'conserve' .and. allocated(c2d)) then + k = k+1; c2d(k) = outvars(n) + end if + end do + + ! -------------------------------------------------------- + ! create packed arrays for mapping + ! -------------------------------------------------------- + + ! 2D bilin + if (allocated(bilin2d)) then + call packarrays(trim(input_file), trim(wgtsdir), cosrot, sinrot, b2d, dims=(/nxt,nyt/), & + nflds=nbilin2d, fields=bilin2d) + + if (debug) then + write(logunit,'(/,a)')'2D fields mapped bilin, packed field min/max values' + do n = 1,nbilin2d + write(logunit,'(i6,4(a,a),2g14.4)')n,' ',trim(b2d(n)%input_var_name),' ', & + trim(b2d(n)%var_grid),' ',trim(b2d(n)%var_pair),' ', trim(b2d(n)%var_pair_grid), & + minval(bilin2d(:,n)), maxval(bilin2d(:,n)) + end do + call dumpnc(trim(ftype)//'.bilin2d.nc', 'bilin2d', dims=(/nxt,nyt/), nflds=nbilin2d, field=bilin2d) + end if + end if + + ! 2D conserv + if (allocated(consd2d)) then + call packarrays(trim(input_file), trim(wgtsdir), cosrot, sinrot, c2d, dims=(/nxt,nyt/), & + nflds=nconsd2d, fields=consd2d) + + if (debug) then + write(logunit,'(a)')'2D fields mapped conserv, packed field min/max values' + do n = 1,nconsd2d + write(logunit,'(i6,4(a,a),2g14.4)')n,' ',trim(c2d(n)%input_var_name),' ', & + trim(c2d(n)%var_grid),' ',trim(c2d(n)%var_pair),' ', trim(c2d(n)%var_pair_grid), & + minval(consd2d(:,n)), maxval(consd2d(:,n)) + end do + call dumpnc(trim(ftype)//'.consd2d.nc', 'consd2d', dims=(/nxt,nyt/), nflds=nconsd2d, field=consd2d) + end if + end if + + ! 3D bilin + if (allocated(bilin3d))then + call packarrays(trim(input_file), trim(wgtsdir), cosrot, sinrot, b3d, dims=(/nxt,nyt,nlevs/), & + nflds=nbilin3d, fields=bilin3d) + + if (debug) then + write(logunit,'(a)')'3D fields mapped bilin, packed field min/max values' + do n = 1,nbilin3d + write(logunit,'(i6,4(a,a),2g14.4)')n,' ',trim(b3d(n)%input_var_name),' ', & + trim(b3d(n)%var_grid),' ',trim(b3d(n)%var_pair),' ', trim(b3d(n)%var_pair_grid), & + minval(bilin3d(:,:,n)), maxval(bilin3d(:,:,n)) + end do + call dumpnc(trim(ftype)//'.bilin3d.nc', 'bilin3d', dims=(/nxt,nyt,nlevs/), nk=nlevs, & + nflds=nbilin3d, field=bilin3d) + end if + end if + + ! -------------------------------------------------------- + ! remap packed arrays to each destination grid + ! -------------------------------------------------------- + + do nd = 1,ndest + dstgrid = trim(dstgrds(nd)) + nxr = nxrs(nd); nyr = nyrs(nd) + + 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 + allocate(out3d(nxr,nyr,nlevs)); out3d = 0.0 + end if + allocate(dstlon(nxr,nyr)); dstlon = 0.0 + allocate(dstlat(nxr,nyr)); dstlat = 0.0 + allocate(out1d(nxr*nyr)); out1d = 0.0 + allocate(out2d(nxr,nyr)); out2d = 0.0 + + if (do_ocnpost) then + allocate(rgmask3d(nxr*nyr,nlevs)); rgmask3d = 0.0 + else + allocate(rgmask2d(nxr*nyr)); rgmask2d = 0.0 + end if + ! lat,lon of destination grid can be obtained from xc_b,yc_b in wgtsfile + wgtsfile = trim(wgtsdir)//'tripole.mx025.Ct.to.rect.'//trim(dstgrid)//'.bilinear.nc' + 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) + + if (allocated(bilin2d)) then + wgtsfile = trim(wgtsdir)//'tripole.mx025.Ct.to.rect.'//trim(dstgrid)//'.bilinear.nc' + if (debug) write(logunit,'(/,a)')'remapping 2D fields bilinear with '//trim(wgtsfile) + call remap(trim(wgtsfile), dim2=nbilin2d, src_field=bilin2d, dst_field=rgb2d) + + if (debug) then + write(logunit,'(a)')'2D fields mapped bilin, mapped field min/max values' + do n = 1,nbilin2d + write(logunit,'(i6,4(a,a),2g14.4)')n,' ',trim(b2d(n)%input_var_name),' ', & + trim(b2d(n)%var_grid),' ',trim(b2d(n)%var_pair),' ', trim(b2d(n)%var_pair_grid), & + minval(rgb2d(:,n)), maxval(rgb2d(:,n)) + end do + call dumpnc(trim(ftype)//'.rgbilin2d.'//trim(dstgrid)//'.nc', 'rgbilin2d', dims=(/nxr,nyr/), & + nflds=nbilin2d, field=rgb2d) + end if + end if + + if (allocated(consd2d)) then + wgtsfile = trim(wgtsdir)//'tripole.mx025.Ct.to.rect.'//trim(dstgrid)//'.conserve.nc' + if (debug) write(logunit,'(a)')'remapping 2D fields conserv with '//trim(wgtsfile) + call remap(trim(wgtsfile), dim2=nconsd2d, src_field=consd2d, dst_field=rgc2d) + + if (debug) then + write(logunit,'(a)')'2D fields mapped conserv, mapped field min/max values' + do n = 1,nconsd2d + write(logunit,'(i6,4(a,a),2g14.4)')n,' ',trim(c2d(n)%input_var_name),' ', & + trim(c2d(n)%var_grid),' ',trim(c2d(n)%var_pair),' ', trim(c2d(n)%var_pair_grid), & + minval(rgc2d(:,n)), maxval(rgc2d(:,n)) + end do + call dumpnc(trim(ftype)//'.rgconsd2d.'//trim(dstgrid)//'.nc', 'rgconsd2d', dims=(/nxr,nyr/), & + nflds=nconsd2d, field=rgc2d) + end if + end if + + if (allocated(bilin3d)) then + wgtsfile = trim(wgtsdir)//'tripole.mx025.Ct.to.rect.'//trim(dstgrid)//'.bilinear.nc' + if (debug) write(logunit,'(a)')'remapping 3D fields bilinear with '//trim(wgtsfile) + call remap(trim(wgtsfile), nk=nlevs, nflds=nbilin3d, src_field=bilin3d, dst_field=rgb3d) + + if (debug) then + write(logunit,'(a)')'3D fields mapped bilin, mapped field min/max values' + do n = 1,nbilin3d + write(logunit,'(i6,4(a,a),2g14.4)')n,' ',trim(b3d(n)%input_var_name),' ', & + trim(b3d(n)%var_grid),' ',trim(b3d(n)%var_pair),' ', trim(b3d(n)%var_pair_grid), & + minval(rgb3d(:,:,n)), maxval(rgb3d(:,:,n)) + end do + call dumpnc(trim(ftype)//'.rgbilin3d.'//trim(dstgrid)//'.nc', 'rgbilin3d', dims=(/nxr,nyr,nlevs/), & + nk=nlevs, nflds=nbilin3d, field=rgb3d) + end if + end if + + ! -------------------------------------------------------- + ! remap the source grid 2/3D mask to obtain the interpolation mask. + ! -------------------------------------------------------- + + wgtsfile = trim(wgtsdir)//'tripole.mx025.Ct.to.rect.'//trim(dstgrid)//'.bilinear.nc' + 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(:) <= -79.75)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)//'.rgmask3d.'//trim(dstgrid)//'.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(:) <= -79.75)rgmask2d(:) = vfill + + if (debug) then + write(logunit,'(a,2g14.4)')'mask min/max on destination grid ',minval(rgmask2d),maxval(rgmask2d) + call dumpnc(trim(ftype)//'.rgmask2d.'//trim(dstgrid)//'.nc', 'rgmask2d', dims=(/nxr,nyr/), & + field=rgmask2d) + 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)%output_var_name) == 'speed')idx1 = n + if (trim(b2d(n)%output_var_name) == 'SSU')idx2 = n + if (trim(b2d(n)%output_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 + ! -------------------------------------------------------- + + output_file = trim(ftype)//'.'//trim(dstgrid)//'.nc' + if (debug) write(logunit, '(a)')'output file: '//trim(output_file) + + rc = nf90_create(trim(output_file), 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)%output_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)%output_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)%output_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)%output_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)%output_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)%output_var_name) + rc = nf90_inq_varid(ncid, vname, varid) + rc = nf90_put_var(ncid, varid, out3d) + end do + end if + rc = nf90_close(ncid) + + if (allocated(rgb2d)) deallocate(rgb2d) + if (allocated(rgc2d)) deallocate(rgc2d) + if (allocated(rgb3d)) deallocate(rgb3d) + + if (allocated(out1d)) deallocate(out1d) + if (allocated(out2d)) deallocate(out2d) + if (allocated(out3d)) deallocate(out3d) + + if (allocated(dstlon)) deallocate(dstlon) + if (allocated(dstlat)) deallocate(dstlat) + + if (allocated(rgmask2d)) deallocate(rgmask2d) + if (allocated(rgmask3d)) deallocate(rgmask3d) + + end do !nd + write(logunit,'(a)')'all done!' + +end program ocnicepost diff --git a/src/ocnicepost.fd/ocnicepost.nml b/src/ocnicepost.fd/ocnicepost.nml new file mode 100644 index 00000000..abe2439b --- /dev/null +++ b/src/ocnicepost.fd/ocnicepost.nml @@ -0,0 +1,5 @@ +& ocnicepost_nml +ftype='ocean' +wgtsdir='/scratch1/NCEPDEV/stmp4/Denise.Worthen/CPLD_GRIDGEN/' +debug=.false. +/ \ No newline at end of file diff --git a/src/ocnicepost.fd/outputvars.F90 b/src/ocnicepost.fd/outputvars.F90 new file mode 100644 index 00000000..8d79e2fb --- /dev/null +++ b/src/ocnicepost.fd/outputvars.F90 @@ -0,0 +1,146 @@ +module outputvars + + implicit none + + integer, parameter :: maxvars = 40 !< The maximum number of variables written to a file + + type :: vardefs + character(len= 10) :: input_var_name !< A variable's input variable name + character(len= 10) :: output_var_name !< A variable's output 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), public :: outvars(maxvars) !< Attribute definitions for the variables + +contains + + subroutine ocnvars_typedefine + + ! local variables + integer :: ii = 0 + + !set defaults + outvars(:)%input_var_name='' + outvars(:)%var_grid = 'Ct' + outvars(:)%var_remapmethod = 'bilinear' + outvars(:)%var_dimen = 2 + outvars(:)%var_pair = '' + outvars(:)%var_pair_grid = '' + outvars(:)%long_name = '' ! obtained from input file + outvars(:)%units = '' ! obtained from input file + outvars(:)%var_fillvalue = -1.0 ! obtained from input file + + ! 2D states with native grid location on cell centers; remapped bilinearly + ii = ii + 1; outvars(ii)%input_var_name = 'SSH' + ii = ii + 1; outvars(ii)%input_var_name = 'SST' + ii = ii + 1; outvars(ii)%input_var_name = 'SSS' + ii = ii + 1; outvars(ii)%input_var_name = 'speed' + !ii = ii + 1; outvars(ii)%input_var_name = 'mld' + ii = ii + 1; outvars(ii)%input_var_name = 'ePBL' + ii = ii + 1; outvars(ii)%input_var_name = 'MLD_003' + ii = ii + 1; outvars(ii)%input_var_name = 'MLD_0125' + + ! 2D fluxes with native grid location on cell centers; remapped conservatively + ii = ii + 1; outvars(ii)%input_var_name = 'latent' ; outvars(ii)%var_remapmethod = 'conserve' + ii = ii + 1; outvars(ii)%input_var_name = 'sensible' ; outvars(ii)%var_remapmethod = 'conserve' + ii = ii + 1; outvars(ii)%input_var_name = 'SW' ; outvars(ii)%var_remapmethod = 'conserve' + ii = ii + 1; outvars(ii)%input_var_name = 'LW' ; outvars(ii)%var_remapmethod = 'conserve' + ii = ii + 1; outvars(ii)%input_var_name = 'evap' ; outvars(ii)%var_remapmethod = 'conserve' + ii = ii + 1; outvars(ii)%input_var_name = 'lprec' ; outvars(ii)%var_remapmethod = 'conserve' + ii = ii + 1; outvars(ii)%input_var_name = 'fprec' ; outvars(ii)%var_remapmethod = 'conserve' + ii = ii + 1; outvars(ii)%input_var_name = 'LwLatSens' ; outvars(ii)%var_remapmethod = 'conserve' + ii = ii + 1; outvars(ii)%input_var_name = 'Heat_PmE' ; outvars(ii)%var_remapmethod = 'conserve' + + ! 2D vector states on stagger locations; remapped bilinearly + ii = ii + 1; outvars(ii)%input_var_name = 'SSU' + outvars(ii)%var_grid = 'Cu' + outvars(ii)%var_pair = 'SSV' + outvars(ii)%var_pair_grid = 'Cv' + + ii = ii + 1; outvars(ii)%input_var_name = 'SSV' + outvars(ii)%var_grid = 'Cv' + outvars(ii)%var_pair = 'SSU' + outvars(ii)%var_pair_grid = 'Cu' + + ! 2D vector fluxes on stagger locations; remapped conservatively + ii = ii + 1; outvars(ii)%input_var_name = 'taux' + outvars(ii)%var_grid = 'Cu' + outvars(ii)%var_pair = 'tauy' + outvars(ii)%var_pair_grid = 'Cv' + outvars(ii)%var_remapmethod = 'conserve' + + ii = ii + 1; outvars(ii)%input_var_name = 'tauy' + outvars(ii)%var_grid = 'Cv' + outvars(ii)%var_pair = 'taux' + outvars(ii)%var_pair_grid = 'Cu' + outvars(ii)%var_remapmethod = 'conserve' + + ! 3D scalars with native grid location on cell centers; remapped bilinearly + ii = ii + 1; outvars(ii)%input_var_name = 'temp' ; outvars(ii)%var_dimen = 3 + ii = ii + 1; outvars(ii)%input_var_name = 'so' ; outvars(ii)%var_dimen = 3 + + ! 3D vectors on stagger locations; remapped bilinearly + ii = ii + 1; outvars(ii)%input_var_name = 'uo' + outvars(ii)%var_grid = 'Cu' + outvars(ii)%var_pair = 'vo' + outvars(ii)%var_pair_grid = 'Cv' + outvars(ii)%var_dimen = 3 + + ii = ii + 1; outvars(ii)%input_var_name = 'vo' + outvars(ii)%var_grid = 'Cv' + outvars(ii)%var_pair = 'uo' + outvars(ii)%var_pair_grid = 'Cu' + outvars(ii)%var_dimen = 3 + + ! set default output name + outvars(:)%output_var_name = outvars(:)%input_var_name + + end subroutine ocnvars_typedefine + + subroutine icevars_typedefine + + ! local variables + integer :: ii = 0 + + !set defaults + outvars(:)%input_var_name='' + outvars(:)%var_grid = 'Ct' + outvars(:)%var_remapmethod = 'bilinear' + outvars(:)%var_dimen = 2 + outvars(:)%var_pair = '' + outvars(:)%var_pair_grid = '' + outvars(:)%long_name = '' ! obtained from input file + outvars(:)%units = '' ! obtained from input file + outvars(:)%var_fillvalue = -1.0 ! obtained from input file + + ! 2D states with native grid location on cell centers; remapped bilinearly + ii = ii + 1; outvars(ii)%input_var_name = 'hi_h' + ii = ii + 1; outvars(ii)%input_var_name = 'hs_h' + ii = ii + 1; outvars(ii)%input_var_name = 'aice_h' + ii = ii + 1; outvars(ii)%input_var_name = 'sst_h' + ii = ii + 1; outvars(ii)%input_var_name = 'Tsfc_h' + + ! 2D vector states on stagger locations; remapped bilinearly + ii = ii + 1; outvars(ii)%input_var_name = 'uvel_h' + outvars(ii)%var_grid = 'Bu_x' + outvars(ii)%var_pair = 'vvel_h' + outvars(ii)%var_pair_grid = 'Bu_y' + + ii = ii + 1; outvars(ii)%input_var_name = 'vvel_h' + outvars(ii)%var_grid = 'Bu_y' + outvars(ii)%var_pair = 'uvel_h' + outvars(ii)%var_pair_grid = 'Bu_x' + + ! set default output name + outvars(:)%output_var_name = outvars(:)%input_var_name + + end subroutine icevars_typedefine + +end module outputvars diff --git a/src/ocnicepost.fd/utils_mod.F90 b/src/ocnicepost.fd/utils_mod.F90 new file mode 100644 index 00000000..0df02b78 --- /dev/null +++ b/src/ocnicepost.fd/utils_mod.F90 @@ -0,0 +1,604 @@ +module utils_mod + + use netcdf + use outputvars, only : vardefs + + implicit none + + private + + logical, public :: debug + integer, public :: logunit + + 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)%input_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)%input_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)%input_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)%input_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 = 'getfield3d' + + 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 :: n, 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 :: n, 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 From 6518d72b72a8639869dfebb3b7babdca09b89bad Mon Sep 17 00:00:00 2001 From: "denise.worthen" Date: Fri, 20 Oct 2023 10:29:02 -0400 Subject: [PATCH 2/4] restructure code, use configurable files * add csv files to set source field info * restrict destination grid to a single resolution --- src/ocnicepost.fd/arrays_mod.F90 | 108 +++++++ src/ocnicepost.fd/init_mod.F90 | 159 ++++++++++ src/ocnicepost.fd/masking_mod.F90 | 163 ++++++++++ src/ocnicepost.fd/ocnicepost.F90 | 460 +++++++--------------------- src/ocnicepost.fd/remapping_mod.F90 | 6 + src/ocnicepost.fd/utils_mod.F90 | 19 +- 6 files changed, 559 insertions(+), 356 deletions(-) create mode 100644 src/ocnicepost.fd/arrays_mod.F90 create mode 100644 src/ocnicepost.fd/init_mod.F90 create mode 100644 src/ocnicepost.fd/masking_mod.F90 create mode 100644 src/ocnicepost.fd/remapping_mod.F90 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..a98ae68e --- /dev/null +++ b/src/ocnicepost.fd/init_mod.F90 @@ -0,0 +1,159 @@ +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..c61b7a7d --- /dev/null +++ b/src/ocnicepost.fd/masking_mod.F90 @@ -0,0 +1,163 @@ +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) + ! use 3D source variable to use as mask. Obtain directly from + ! file to obtain 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) + ! use 2D source variable to use as mask. Obtain directly from + ! file to obtain 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 index ca21bf3a..1aa478d3 100644 --- a/src/ocnicepost.fd/ocnicepost.F90 +++ b/src/ocnicepost.fd/ocnicepost.F90 @@ -1,118 +1,60 @@ program ocnicepost - use utils_mod, only : debug, logunit, getfield, packarrays, remap, dumpnc - use outputvars 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= 20) :: fnml, ftype - character(len=120) :: wgtsdir - character(len=120) :: input_file, wgtsfile, output_file - - ! source grid, tripole 1/4 deg, 40 vertical levels - integer, parameter :: nxt = 1440, nyt = 1080, nlevs = 40 - - ! destination grids - integer, parameter :: ndest = 3 - integer, parameter, dimension(ndest) :: nxrs = (/1440, 720, 360/) - integer, parameter, dimension(ndest) :: nyrs = (/ 721, 361, 181/) - character(len=4), dimension(ndest) :: dstgrds = (/'0p25', '0p5 ', '1p0 '/) - - ! packed arrays on source grid - real, allocatable, dimension(:,:) :: bilin2d, consd2d - real, allocatable, dimension(:,:,:) :: bilin3d - ! variable types - type(vardefs), allocatable :: b2d(:) - type(vardefs), allocatable :: c2d(:) - type(vardefs), allocatable :: b3d(:) - - ! source grid fields - real, dimension(nxt,nyt,nlevs) :: tmp3d - real, dimension(nxt*nyt,nlevs) :: mask3d - real, dimension(nxt*nyt) :: mask2d - real, dimension(nxt*nyt) :: cosrot, sinrot, anglet - real, dimension(nlevs) :: z_l - real, dimension(0:nlevs) :: z_i - - ! destination grid fields - real, allocatable, dimension(:,:) :: dstlon, dstlat - real, allocatable, dimension(:,:) :: rgmask3d - real, allocatable, dimension(:) :: rgmask2d - ! output fields (x,y,z) - real, allocatable, dimension(:) :: out1d - real, allocatable, dimension(:,:) :: out2d - real, allocatable, dimension(:,:,:) :: out3d - ! packed remapped fields - real, allocatable, dimension(:,:) :: rgb2d, rgc2d - real, allocatable, dimension(:,:,:) :: rgb3d + 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 - character(len=4) :: dstgrid - real :: vfill - integer :: nd, nxr, nyr - integer :: i,j,k,n,nn,nvalid,iounit - integer :: rc,ncid,varid,dimid - integer :: nbilin2d,nbilin3d,nconsd2d + real :: vfill + integer :: nvalid + integer :: n,rc,ncid,varid integer :: idimid,jdimid,kdimid,edimid,timid integer :: idx1,idx2,idx3 - logical :: do_ocnpost - - namelist /ocnicepost_nml/ ftype, wgtsdir, debug ! -------------------------------------------------------- - ! read the name list + ! read the nml file and a file containing the list of + ! variables to be remapped ! -------------------------------------------------------- - fnml = 'ocnicepost.nml' - inquire (file=trim(fnml), iostat=rc) - if (rc /= 0) then - write (6, '(3a)') 'Error: input file "', trim(fnml), '" does not exist.' - stop - end if - - ! Open and read Namelist file. - open (action='read', file=trim(fnml), iostat=rc, newunit=iounit) - read (nml=ocnicepost_nml, iostat=rc, unit=iounit) - if (rc /= 0) then - write (6, '(a)') 'Error: invalid Namelist format.' - end if - close (iounit) - - ! initialize the source file type and variables - if (trim(ftype) == 'ocean') then - do_ocnpost = .true. - call ocnvars_typedefine - else - do_ocnpost = .false. - call icevars_typedefine - 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) + call readnml + call readcsv(nvalid) ! -------------------------------------------------------- ! read the source file and obtain the units and long name, ! rotation angles, vertical grid and time axis ! -------------------------------------------------------- - nvalid = 0 rc = nf90_open(trim(input_file), nf90_nowrite, ncid) - do i = 1,maxvars - if (len_trim(outvars(i)%input_var_name) > 0 ) then - rc = nf90_inq_varid(ncid, trim(outvars(i)%input_var_name), varid) - rc = nf90_get_att(ncid, varid, 'long_name', outvars(i)%long_name) - rc = nf90_get_att(ncid, varid, 'units', outvars(i)%units) - rc = nf90_get_att(ncid, varid, '_FillValue', outvars(i)%var_fillvalue) - nvalid = nvalid+1 - if (trim(outvars(i)%input_var_name) == 'temp')vfill = outvars(i)%var_fillvalue - if (trim(outvars(i)%input_var_name) == 'aice_h')vfill = outvars(i)%var_fillvalue - end if + 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 @@ -121,6 +63,11 @@ program ocnicepost 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) @@ -129,279 +76,114 @@ program ocnicepost rc = nf90_get_var(ncid, varid, z_i) rc = nf90_close(ncid) ! rotation angles - call getfield(trim(input_file), 'cos_rot', dims=(/nxt,nyt/), field=cosrot) - call getfield(trim(input_file), 'sin_rot', dims=(/nxt,nyt/), field=sinrot) + 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 - call getfield(trim(input_file), 'ANGLET', dims=(/nxt,nyt/), field=anglet) + 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 + ! -------------------------------------------------------- - ! 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 + ! create interpolation masks ! -------------------------------------------------------- if (do_ocnpost) then - rc = nf90_open(trim(input_file), nf90_nowrite, ncid) - ! 3D temp to use as mask, obtain directly from file to preserve vfill - rc = nf90_inq_varid(ncid, 'temp', varid) - 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 + allocate(mask3d(nxt*nyt,nlevs)); mask3d = 0.0 + allocate(rgmask3d(nxr*nyr,nlevs)); rgmask3d = 0.0 else - call getfield(trim(input_file), 'tmask', dims=(/nxt,nyt/), field=mask2d) - ! 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 + allocate(mask2d(nxt*nyt)); mask2d = 0.0 + allocate(rgmask2d(nxr*nyr)); rgmask2d = 0.0 end if - ! -------------------------------------------------------- - ! 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(outvars(n)%var_remapmethod) == 'bilinear') then - if (outvars(n)%var_dimen == 2) nbilin2d = nbilin2d + 1 - if (outvars(n)%var_dimen == 3) nbilin3d = nbilin3d + 1 - end if - if (trim(outvars(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 + call remap_masks(vfill) ! -------------------------------------------------------- - ! create types for each packed array + ! create packed arrays for mapping and remap packed arrays + ! to the destination grid ! -------------------------------------------------------- - i = 0; j = 0; k = 0 - do n = 1,nvalid - if (trim(outvars(n)%var_remapmethod) == 'bilinear') then - if (outvars(n)%var_dimen == 2 .and. allocated(b2d)) then - i = i+1; b2d(i) = outvars(n) - end if - if (outvars(n)%var_dimen == 3 .and. allocated(b3d)) then - j = j+1; b3d(j) = outvars(n) - end if - end if - if (trim(outvars(n)%var_remapmethod) == 'conserve' .and. allocated(c2d)) then - k = k+1; c2d(k) = outvars(n) - end if - end do - - ! -------------------------------------------------------- - ! create packed arrays for mapping - ! -------------------------------------------------------- + call setup_packing(nvalid,outvars) ! 2D bilin if (allocated(bilin2d)) then - call packarrays(trim(input_file), trim(wgtsdir), cosrot, sinrot, b2d, dims=(/nxt,nyt/), & + + 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)')'2D fields mapped bilin, packed field min/max values' + 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,'(i6,4(a,a),2g14.4)')n,' ',trim(b2d(n)%input_var_name),' ', & - trim(b2d(n)%var_grid),' ',trim(b2d(n)%var_pair),' ', trim(b2d(n)%var_pair_grid), & - minval(bilin2d(:,n)), maxval(bilin2d(:,n)) + 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)//'.bilin2d.nc', 'bilin2d', dims=(/nxt,nyt/), nflds=nbilin2d, field=bilin2d) + 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 - call packarrays(trim(input_file), trim(wgtsdir), cosrot, sinrot, c2d, dims=(/nxt,nyt/), & + + 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)')'2D fields mapped conserv, packed field min/max values' + 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,'(i6,4(a,a),2g14.4)')n,' ',trim(c2d(n)%input_var_name),' ', & - trim(c2d(n)%var_grid),' ',trim(c2d(n)%var_pair),' ', trim(c2d(n)%var_pair_grid), & - minval(consd2d(:,n)), maxval(consd2d(:,n)) + 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)//'.consd2d.nc', 'consd2d', dims=(/nxt,nyt/), nflds=nconsd2d, field=consd2d) + 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 - call packarrays(trim(input_file), trim(wgtsdir), cosrot, sinrot, b3d, dims=(/nxt,nyt,nlevs/), & + + 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)')'3D fields mapped bilin, packed field min/max values' + 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,'(i6,4(a,a),2g14.4)')n,' ',trim(b3d(n)%input_var_name),' ', & - trim(b3d(n)%var_grid),' ',trim(b3d(n)%var_pair),' ', trim(b3d(n)%var_pair_grid), & - minval(bilin3d(:,:,n)), maxval(bilin3d(:,:,n)) + 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)//'.bilin3d.nc', 'bilin3d', dims=(/nxt,nyt,nlevs/), nk=nlevs, & - nflds=nbilin3d, field=bilin3d) + 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 - ! -------------------------------------------------------- - ! remap packed arrays to each destination grid - ! -------------------------------------------------------- - - do nd = 1,ndest - dstgrid = trim(dstgrds(nd)) - nxr = nxrs(nd); nyr = nyrs(nd) - - 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 - allocate(out3d(nxr,nyr,nlevs)); out3d = 0.0 - end if - allocate(dstlon(nxr,nyr)); dstlon = 0.0 - allocate(dstlat(nxr,nyr)); dstlat = 0.0 - allocate(out1d(nxr*nyr)); out1d = 0.0 - allocate(out2d(nxr,nyr)); out2d = 0.0 - - if (do_ocnpost) then - allocate(rgmask3d(nxr*nyr,nlevs)); rgmask3d = 0.0 - else - allocate(rgmask2d(nxr*nyr)); rgmask2d = 0.0 - end if - ! lat,lon of destination grid can be obtained from xc_b,yc_b in wgtsfile - wgtsfile = trim(wgtsdir)//'tripole.mx025.Ct.to.rect.'//trim(dstgrid)//'.bilinear.nc' - 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) - - if (allocated(bilin2d)) then - wgtsfile = trim(wgtsdir)//'tripole.mx025.Ct.to.rect.'//trim(dstgrid)//'.bilinear.nc' - if (debug) write(logunit,'(/,a)')'remapping 2D fields bilinear with '//trim(wgtsfile) - call remap(trim(wgtsfile), dim2=nbilin2d, src_field=bilin2d, dst_field=rgb2d) - - if (debug) then - write(logunit,'(a)')'2D fields mapped bilin, mapped field min/max values' - do n = 1,nbilin2d - write(logunit,'(i6,4(a,a),2g14.4)')n,' ',trim(b2d(n)%input_var_name),' ', & - trim(b2d(n)%var_grid),' ',trim(b2d(n)%var_pair),' ', trim(b2d(n)%var_pair_grid), & - minval(rgb2d(:,n)), maxval(rgb2d(:,n)) - end do - call dumpnc(trim(ftype)//'.rgbilin2d.'//trim(dstgrid)//'.nc', 'rgbilin2d', dims=(/nxr,nyr/), & - nflds=nbilin2d, field=rgb2d) - end if - end if - - if (allocated(consd2d)) then - wgtsfile = trim(wgtsdir)//'tripole.mx025.Ct.to.rect.'//trim(dstgrid)//'.conserve.nc' - if (debug) write(logunit,'(a)')'remapping 2D fields conserv with '//trim(wgtsfile) - call remap(trim(wgtsfile), dim2=nconsd2d, src_field=consd2d, dst_field=rgc2d) - - if (debug) then - write(logunit,'(a)')'2D fields mapped conserv, mapped field min/max values' - do n = 1,nconsd2d - write(logunit,'(i6,4(a,a),2g14.4)')n,' ',trim(c2d(n)%input_var_name),' ', & - trim(c2d(n)%var_grid),' ',trim(c2d(n)%var_pair),' ', trim(c2d(n)%var_pair_grid), & - minval(rgc2d(:,n)), maxval(rgc2d(:,n)) - end do - call dumpnc(trim(ftype)//'.rgconsd2d.'//trim(dstgrid)//'.nc', 'rgconsd2d', dims=(/nxr,nyr/), & - nflds=nconsd2d, field=rgc2d) - end if - end if - - if (allocated(bilin3d)) then - wgtsfile = trim(wgtsdir)//'tripole.mx025.Ct.to.rect.'//trim(dstgrid)//'.bilinear.nc' - if (debug) write(logunit,'(a)')'remapping 3D fields bilinear with '//trim(wgtsfile) - call remap(trim(wgtsfile), nk=nlevs, nflds=nbilin3d, src_field=bilin3d, dst_field=rgb3d) - - if (debug) then - write(logunit,'(a)')'3D fields mapped bilin, mapped field min/max values' - do n = 1,nbilin3d - write(logunit,'(i6,4(a,a),2g14.4)')n,' ',trim(b3d(n)%input_var_name),' ', & - trim(b3d(n)%var_grid),' ',trim(b3d(n)%var_pair),' ', trim(b3d(n)%var_pair_grid), & - minval(rgb3d(:,:,n)), maxval(rgb3d(:,:,n)) - end do - call dumpnc(trim(ftype)//'.rgbilin3d.'//trim(dstgrid)//'.nc', 'rgbilin3d', dims=(/nxr,nyr,nlevs/), & - nk=nlevs, nflds=nbilin3d, field=rgb3d) - end if - end if - - ! -------------------------------------------------------- - ! remap the source grid 2/3D mask to obtain the interpolation mask. - ! -------------------------------------------------------- - - wgtsfile = trim(wgtsdir)//'tripole.mx025.Ct.to.rect.'//trim(dstgrid)//'.bilinear.nc' - 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(:) <= -79.75)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)//'.rgmask3d.'//trim(dstgrid)//'.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(:) <= -79.75)rgmask2d(:) = vfill - - if (debug) then - write(logunit,'(a,2g14.4)')'mask min/max on destination grid ',minval(rgmask2d),maxval(rgmask2d) - call dumpnc(trim(ftype)//'.rgmask2d.'//trim(dstgrid)//'.nc', 'rgmask2d', dims=(/nxr,nyr/), & - field=rgmask2d) - end if - end if - ! -------------------------------------------------------- ! mask the mapped fields ! -------------------------------------------------------- @@ -427,6 +209,7 @@ program ocnicepost where(rgmask3d(:,:) .eq. vfill)rgb3d(:,:,n) = vfill end if end do + ! -------------------------------------------------------- ! replace model native speed field with a value calculated ! from remapped ssu,ssv @@ -434,9 +217,9 @@ program ocnicepost if (do_ocnpost) then do n = 1,nbilin2d - if (trim(b2d(n)%output_var_name) == 'speed')idx1 = n - if (trim(b2d(n)%output_var_name) == 'SSU')idx2 = n - if (trim(b2d(n)%output_var_name) == 'SSV')idx3 = n + 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) @@ -446,10 +229,13 @@ program ocnicepost ! write the mapped fields ! -------------------------------------------------------- - output_file = trim(ftype)//'.'//trim(dstgrid)//'.nc' - if (debug) write(logunit, '(a)')'output file: '//trim(output_file) + 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(output_file), nf90_clobber, ncid) + 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) @@ -477,7 +263,7 @@ program ocnicepost if (allocated(b2d)) then do n = 1,nbilin2d - vname = trim(b2d(n)%output_var_name) + vname = trim(b2d(n)%var_name) vunit = trim(b2d(n)%units) vlong = trim(b2d(n)%long_name) vfill = b2d(n)%var_fillvalue @@ -489,7 +275,7 @@ program ocnicepost end if if (allocated(c2d)) then do n = 1,nconsd2d - vname = trim(c2d(n)%output_var_name) + vname = trim(c2d(n)%var_name) vunit = trim(c2d(n)%units) vlong = trim(c2d(n)%long_name) vfill = c2d(n)%var_fillvalue @@ -501,7 +287,7 @@ program ocnicepost end if if (allocated(b3d)) then do n = 1,nbilin3d - vname = trim(b3d(n)%output_var_name) + vname = trim(b3d(n)%var_name) vunit = trim(b3d(n)%units) vlong = trim(b3d(n)%long_name) vfill = b3d(n)%var_fillvalue @@ -532,7 +318,7 @@ program ocnicepost do n = 1,nbilin2d out2d(:,:) = reshape(rgb2d(:,n), (/nxr,nyr/)) out2d(:,nyr) = vfill - vname = trim(b2d(n)%output_var_name) + vname = trim(b2d(n)%var_name) rc = nf90_inq_varid(ncid, vname, varid) rc = nf90_put_var(ncid, varid, out2d) end do @@ -541,7 +327,7 @@ program ocnicepost do n = 1,nconsd2d out2d(:,:) = reshape(rgc2d(:,n), (/nxr,nyr/)) out2d(:,nyr) = vfill - vname = trim(c2d(n)%output_var_name) + vname = trim(c2d(n)%var_name) rc = nf90_inq_varid(ncid, vname, varid) rc = nf90_put_var(ncid, varid, out2d) end do @@ -550,28 +336,12 @@ program ocnicepost do n = 1,nbilin3d out3d(:,:,:) = reshape(rgb3d(:,:,n), (/nxr,nyr,nlevs/)) out3d(:,nyr,:) = vfill - vname = trim(b3d(n)%output_var_name) + 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) - - if (allocated(rgb2d)) deallocate(rgb2d) - if (allocated(rgc2d)) deallocate(rgc2d) - if (allocated(rgb3d)) deallocate(rgb3d) - - if (allocated(out1d)) deallocate(out1d) - if (allocated(out2d)) deallocate(out2d) - if (allocated(out3d)) deallocate(out3d) - - if (allocated(dstlon)) deallocate(dstlon) - if (allocated(dstlat)) deallocate(dstlat) - - if (allocated(rgmask2d)) deallocate(rgmask2d) - if (allocated(rgmask3d)) deallocate(rgmask3d) - - end do !nd - write(logunit,'(a)')'all done!' + write(logunit,'(a)')trim(fout)//' done' end program ocnicepost diff --git a/src/ocnicepost.fd/remapping_mod.F90 b/src/ocnicepost.fd/remapping_mod.F90 new file mode 100644 index 00000000..da51b924 --- /dev/null +++ b/src/ocnicepost.fd/remapping_mod.F90 @@ -0,0 +1,6 @@ +module remapping_mod + +contains + + call remap +end module remapping_mod diff --git a/src/ocnicepost.fd/utils_mod.F90 b/src/ocnicepost.fd/utils_mod.F90 index 0df02b78..82183c5e 100644 --- a/src/ocnicepost.fd/utils_mod.F90 +++ b/src/ocnicepost.fd/utils_mod.F90 @@ -1,15 +1,12 @@ module utils_mod use netcdf - use outputvars, only : vardefs + use init_mod, only : debug, logunit, vardefs implicit none private - logical, public :: debug - integer, public :: logunit - interface getfield module procedure getfield2d module procedure getfield3d @@ -70,7 +67,7 @@ subroutine packarrays2d(filesrc, wgtsdir, cosrot, sinrot, vars, dims, nflds, fie 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)%input_var_name), trim(vars(n)%var_grid(1:2)), & + 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 @@ -81,7 +78,7 @@ subroutine packarrays2d(filesrc, wgtsdir, cosrot, sinrot, vars, dims, nflds, fie do n = 1,nflds if (len_trim(vars(n)%var_pair) == 0) then nn = nn + 1 - call getfield(trim(filesrc), trim(vars(n)%input_var_name), dims=(/dims(1),dims(2)/), & + 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 @@ -122,7 +119,7 @@ subroutine packarrays3d(filesrc, wgtsdir, cosrot, sinrot, vars, dims, nflds, fie 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)%input_var_name), trim(vars(n)%var_grid), & + 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 @@ -133,7 +130,7 @@ subroutine packarrays3d(filesrc, wgtsdir, cosrot, sinrot, vars, dims, nflds, fie do n = 1,nflds if (len_trim(vars(n)%var_pair) == 0) then nn = nn + 1 - call getfield(trim(filesrc), trim(vars(n)%input_var_name), dims=(/dims(1),dims(2),dims(3)/), & + 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 @@ -201,7 +198,7 @@ subroutine getvecpair3d(fname, wdir, cosrot, sinrot, vname1, vgrid1, & integer :: ii,k real, dimension(dims(1)*dims(2)) :: urot, vrot character(len=240) :: wgtsfile - character(len=20) :: subname = 'getfield3d' + character(len=20) :: subname = 'getvecpair3d' if (debug)write(logunit,'(a)')'enter '//trim(subname) @@ -536,7 +533,7 @@ subroutine dumpnc3dk(fname, vname, dims, field) real, intent(in) :: field(:,:) ! local variable - integer :: n, ncid, varid, rc, idimid, jdimid, kdimid + integer :: ncid, varid, rc, idimid, jdimid, kdimid real, allocatable :: a3d(:,:,:) character(len=20) :: subname = 'dumpnc3dk' @@ -568,7 +565,7 @@ subroutine dumpnc1d(fname, vname, dims, field) real, intent(in) :: field(:) ! local variable - integer :: n, ncid, varid, rc, idimid, jdimid + integer :: ncid, varid, rc, idimid, jdimid real, allocatable :: a2d(:,:) character(len=20) :: subname = 'dumpnc1d' From 7709ef822e0f57ad0911b8613653152070dd86bd Mon Sep 17 00:00:00 2001 From: "denise.worthen" Date: Fri, 20 Oct 2023 10:42:33 -0400 Subject: [PATCH 3/4] tweaks, remove unused files --- src/ocnicepost.fd/init_mod.F90 | 1 - src/ocnicepost.fd/masking_mod.F90 | 6 +- src/ocnicepost.fd/ocnicepost.F90 | 312 ++++++++++++++-------------- src/ocnicepost.fd/ocnicepost.nml | 5 - src/ocnicepost.fd/outputvars.F90 | 146 ------------- src/ocnicepost.fd/remapping_mod.F90 | 6 - 6 files changed, 158 insertions(+), 318 deletions(-) delete mode 100644 src/ocnicepost.fd/ocnicepost.nml delete mode 100644 src/ocnicepost.fd/outputvars.F90 delete mode 100644 src/ocnicepost.fd/remapping_mod.F90 diff --git a/src/ocnicepost.fd/init_mod.F90 b/src/ocnicepost.fd/init_mod.F90 index a98ae68e..39631f59 100644 --- a/src/ocnicepost.fd/init_mod.F90 +++ b/src/ocnicepost.fd/init_mod.F90 @@ -43,7 +43,6 @@ module init_mod 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 diff --git a/src/ocnicepost.fd/masking_mod.F90 b/src/ocnicepost.fd/masking_mod.F90 index c61b7a7d..a2999596 100644 --- a/src/ocnicepost.fd/masking_mod.F90 +++ b/src/ocnicepost.fd/masking_mod.F90 @@ -111,8 +111,7 @@ subroutine makemask3d(vfill) allocate(tmp3d(nxt,nyt,nlevs)); tmp3d = 0.0 rc = nf90_open(trim(input_file), nf90_nowrite, ncid) - ! use 3D source variable to use as mask. Obtain directly from - ! file to obtain fill value + ! 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) @@ -141,8 +140,7 @@ subroutine makemask2d(vfill) allocate(tmp2d(nxt,nyt)); tmp2d = 0.0 rc = nf90_open(trim(input_file), nf90_nowrite, ncid) - ! use 2D source variable to use as mask. Obtain directly from - ! file to obtain fill value + ! 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) diff --git a/src/ocnicepost.fd/ocnicepost.F90 b/src/ocnicepost.fd/ocnicepost.F90 index 1aa478d3..a7c5f52f 100644 --- a/src/ocnicepost.fd/ocnicepost.F90 +++ b/src/ocnicepost.fd/ocnicepost.F90 @@ -30,7 +30,7 @@ program ocnicepost character(len= 20) :: vname, vunit character(len=120) :: vlong - real :: vfill + real :: vfill integer :: nvalid integer :: n,rc,ncid,varid integer :: idimid,jdimid,kdimid,edimid,timid @@ -126,7 +126,7 @@ program ocnicepost 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),' ', & + 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 @@ -150,7 +150,7 @@ program ocnicepost 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),' ', & + 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 @@ -173,7 +173,7 @@ program ocnicepost 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),' ', & + 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 @@ -184,164 +184,164 @@ program ocnicepost end if end if - ! -------------------------------------------------------- - ! mask the mapped fields - ! -------------------------------------------------------- + ! -------------------------------------------------------- + ! 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 - 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 + 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 - 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 + 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 - if (allocated(rgmask3d)) then - where(rgmask3d(:,:) .eq. vfill)rgb3d(:,:,n) = vfill - end if + 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 - - ! -------------------------------------------------------- - ! 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 if + rc = nf90_close(ncid) + write(logunit,'(a)')trim(fout)//' done' end program ocnicepost diff --git a/src/ocnicepost.fd/ocnicepost.nml b/src/ocnicepost.fd/ocnicepost.nml deleted file mode 100644 index abe2439b..00000000 --- a/src/ocnicepost.fd/ocnicepost.nml +++ /dev/null @@ -1,5 +0,0 @@ -& ocnicepost_nml -ftype='ocean' -wgtsdir='/scratch1/NCEPDEV/stmp4/Denise.Worthen/CPLD_GRIDGEN/' -debug=.false. -/ \ No newline at end of file diff --git a/src/ocnicepost.fd/outputvars.F90 b/src/ocnicepost.fd/outputvars.F90 deleted file mode 100644 index 8d79e2fb..00000000 --- a/src/ocnicepost.fd/outputvars.F90 +++ /dev/null @@ -1,146 +0,0 @@ -module outputvars - - implicit none - - integer, parameter :: maxvars = 40 !< The maximum number of variables written to a file - - type :: vardefs - character(len= 10) :: input_var_name !< A variable's input variable name - character(len= 10) :: output_var_name !< A variable's output 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), public :: outvars(maxvars) !< Attribute definitions for the variables - -contains - - subroutine ocnvars_typedefine - - ! local variables - integer :: ii = 0 - - !set defaults - outvars(:)%input_var_name='' - outvars(:)%var_grid = 'Ct' - outvars(:)%var_remapmethod = 'bilinear' - outvars(:)%var_dimen = 2 - outvars(:)%var_pair = '' - outvars(:)%var_pair_grid = '' - outvars(:)%long_name = '' ! obtained from input file - outvars(:)%units = '' ! obtained from input file - outvars(:)%var_fillvalue = -1.0 ! obtained from input file - - ! 2D states with native grid location on cell centers; remapped bilinearly - ii = ii + 1; outvars(ii)%input_var_name = 'SSH' - ii = ii + 1; outvars(ii)%input_var_name = 'SST' - ii = ii + 1; outvars(ii)%input_var_name = 'SSS' - ii = ii + 1; outvars(ii)%input_var_name = 'speed' - !ii = ii + 1; outvars(ii)%input_var_name = 'mld' - ii = ii + 1; outvars(ii)%input_var_name = 'ePBL' - ii = ii + 1; outvars(ii)%input_var_name = 'MLD_003' - ii = ii + 1; outvars(ii)%input_var_name = 'MLD_0125' - - ! 2D fluxes with native grid location on cell centers; remapped conservatively - ii = ii + 1; outvars(ii)%input_var_name = 'latent' ; outvars(ii)%var_remapmethod = 'conserve' - ii = ii + 1; outvars(ii)%input_var_name = 'sensible' ; outvars(ii)%var_remapmethod = 'conserve' - ii = ii + 1; outvars(ii)%input_var_name = 'SW' ; outvars(ii)%var_remapmethod = 'conserve' - ii = ii + 1; outvars(ii)%input_var_name = 'LW' ; outvars(ii)%var_remapmethod = 'conserve' - ii = ii + 1; outvars(ii)%input_var_name = 'evap' ; outvars(ii)%var_remapmethod = 'conserve' - ii = ii + 1; outvars(ii)%input_var_name = 'lprec' ; outvars(ii)%var_remapmethod = 'conserve' - ii = ii + 1; outvars(ii)%input_var_name = 'fprec' ; outvars(ii)%var_remapmethod = 'conserve' - ii = ii + 1; outvars(ii)%input_var_name = 'LwLatSens' ; outvars(ii)%var_remapmethod = 'conserve' - ii = ii + 1; outvars(ii)%input_var_name = 'Heat_PmE' ; outvars(ii)%var_remapmethod = 'conserve' - - ! 2D vector states on stagger locations; remapped bilinearly - ii = ii + 1; outvars(ii)%input_var_name = 'SSU' - outvars(ii)%var_grid = 'Cu' - outvars(ii)%var_pair = 'SSV' - outvars(ii)%var_pair_grid = 'Cv' - - ii = ii + 1; outvars(ii)%input_var_name = 'SSV' - outvars(ii)%var_grid = 'Cv' - outvars(ii)%var_pair = 'SSU' - outvars(ii)%var_pair_grid = 'Cu' - - ! 2D vector fluxes on stagger locations; remapped conservatively - ii = ii + 1; outvars(ii)%input_var_name = 'taux' - outvars(ii)%var_grid = 'Cu' - outvars(ii)%var_pair = 'tauy' - outvars(ii)%var_pair_grid = 'Cv' - outvars(ii)%var_remapmethod = 'conserve' - - ii = ii + 1; outvars(ii)%input_var_name = 'tauy' - outvars(ii)%var_grid = 'Cv' - outvars(ii)%var_pair = 'taux' - outvars(ii)%var_pair_grid = 'Cu' - outvars(ii)%var_remapmethod = 'conserve' - - ! 3D scalars with native grid location on cell centers; remapped bilinearly - ii = ii + 1; outvars(ii)%input_var_name = 'temp' ; outvars(ii)%var_dimen = 3 - ii = ii + 1; outvars(ii)%input_var_name = 'so' ; outvars(ii)%var_dimen = 3 - - ! 3D vectors on stagger locations; remapped bilinearly - ii = ii + 1; outvars(ii)%input_var_name = 'uo' - outvars(ii)%var_grid = 'Cu' - outvars(ii)%var_pair = 'vo' - outvars(ii)%var_pair_grid = 'Cv' - outvars(ii)%var_dimen = 3 - - ii = ii + 1; outvars(ii)%input_var_name = 'vo' - outvars(ii)%var_grid = 'Cv' - outvars(ii)%var_pair = 'uo' - outvars(ii)%var_pair_grid = 'Cu' - outvars(ii)%var_dimen = 3 - - ! set default output name - outvars(:)%output_var_name = outvars(:)%input_var_name - - end subroutine ocnvars_typedefine - - subroutine icevars_typedefine - - ! local variables - integer :: ii = 0 - - !set defaults - outvars(:)%input_var_name='' - outvars(:)%var_grid = 'Ct' - outvars(:)%var_remapmethod = 'bilinear' - outvars(:)%var_dimen = 2 - outvars(:)%var_pair = '' - outvars(:)%var_pair_grid = '' - outvars(:)%long_name = '' ! obtained from input file - outvars(:)%units = '' ! obtained from input file - outvars(:)%var_fillvalue = -1.0 ! obtained from input file - - ! 2D states with native grid location on cell centers; remapped bilinearly - ii = ii + 1; outvars(ii)%input_var_name = 'hi_h' - ii = ii + 1; outvars(ii)%input_var_name = 'hs_h' - ii = ii + 1; outvars(ii)%input_var_name = 'aice_h' - ii = ii + 1; outvars(ii)%input_var_name = 'sst_h' - ii = ii + 1; outvars(ii)%input_var_name = 'Tsfc_h' - - ! 2D vector states on stagger locations; remapped bilinearly - ii = ii + 1; outvars(ii)%input_var_name = 'uvel_h' - outvars(ii)%var_grid = 'Bu_x' - outvars(ii)%var_pair = 'vvel_h' - outvars(ii)%var_pair_grid = 'Bu_y' - - ii = ii + 1; outvars(ii)%input_var_name = 'vvel_h' - outvars(ii)%var_grid = 'Bu_y' - outvars(ii)%var_pair = 'uvel_h' - outvars(ii)%var_pair_grid = 'Bu_x' - - ! set default output name - outvars(:)%output_var_name = outvars(:)%input_var_name - - end subroutine icevars_typedefine - -end module outputvars diff --git a/src/ocnicepost.fd/remapping_mod.F90 b/src/ocnicepost.fd/remapping_mod.F90 deleted file mode 100644 index da51b924..00000000 --- a/src/ocnicepost.fd/remapping_mod.F90 +++ /dev/null @@ -1,6 +0,0 @@ -module remapping_mod - -contains - - call remap -end module remapping_mod From 0a404feb83cdba3dc52130bbef27275bc3304fef Mon Sep 17 00:00:00 2001 From: "denise.worthen" Date: Fri, 20 Oct 2023 12:36:55 -0400 Subject: [PATCH 4/4] add program description --- src/ocnicepost.fd/ocnicepost.F90 | 24 ++++++++++++++++++++++++ 1 file changed, 24 insertions(+) diff --git a/src/ocnicepost.fd/ocnicepost.F90 b/src/ocnicepost.fd/ocnicepost.F90 index a7c5f52f..5669f5d0 100644 --- a/src/ocnicepost.fd/ocnicepost.F90 +++ b/src/ocnicepost.fd/ocnicepost.F90 @@ -1,5 +1,29 @@ 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