Merge pull request #25 from DeniseWorthen/feature/addfortranpost

Add fortran code for ice or ocn post

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

src/ocnicepost.fd/init_mod.F90

@@ -0,0 +1,158 @@
+module init_mod
+
+  implicit none
+
+  public
+
+  real, parameter :: maxvars = 50              !< The maximum number of fields expected in a source file
+
+  type :: vardefs
+     character(len= 10)   :: var_name          !< A variable's variable name
+     character(len=120)   :: long_name         !< A variable's long name
+     character(len= 20)   :: units             !< A variable's unit
+     character(len= 10)   :: var_remapmethod   !< A variable's mapping method
+     integer              :: var_dimen         !< A variable's dimensionality
+     character(len=  4)   :: var_grid          !< A variable's input grid location; all output locations are on cell centers
+     character(len= 10)   :: var_pair          !< A variable's pair
+     character(len=  4)   :: var_pair_grid     !< A pair variable grid
+     real                 :: var_fillvalue     !< A variable's fillvalue
+  end type vardefs
+
+  type(vardefs) :: outvars(maxvars)            !< An empty structure filled by reading a csv file describing the fields
+
+  character(len=10)  :: ftype      !< The type of tripole grid (ocean or ice)
+  character(len=10)  :: fsrc       !< A character string for tripole grid
+  character(len=10)  :: fdst       !< A character string for the destination grid
+  character(len=120) :: wgtsdir    !< The directory containing the regridding weights
+  character(len=20)  :: input_file !< The input file name
+  character(len=10)  :: maskvar    !< The variable in the source file used to create the interpolation mask
+
+  ! rotation angles
+  character(len=10)  :: sinvar    !< The variable in the source file containing the cosine of the rotation angle (ocean only)
+  character(len=10)  :: cosvar    !< The variable in the source file containing the sine of the rotation angle (ocean only)
+  character(len=10)  :: angvar    !< The variable in the source file containing the rotation angle (ice only)
+
+  integer :: nxt        !< The x-dimension of the source tripole grid
+  integer :: nyt        !< The y-dimension of the source tripole grid
+  integer :: nlevs      !< The vertical dimension of the source tripole grid
+
+  integer :: nxr        !< The x-dimension of the destination rectilinear grid
+  integer :: nyr        !< The y-dimension of the destination rectilinear grid
+
+  integer :: logunit    !< The log unit
+  logical :: debug      !< If true, print debug messages and intermediate files
+  logical :: do_ocnpost !< If true, the source file is ocean, otherwise ice
+
+contains
+
+  subroutine readnml
+
+    ! local variable
+    character(len=40) :: fname
+    integer :: ierr, iounit
+    integer :: srcdims(3), dstdims(2)
+
+    namelist /ocnicepost_nml/ ftype, srcdims, wgtsdir, dstdims, maskvar, sinvar, cosvar, &
+         angvar, debug
+
+    ! --------------------------------------------------------
+    ! read the name list
+    ! --------------------------------------------------------
+
+    srcdims = 0; dstdims = 0
+    sinvar=''; cosvar=''; angvar=''
+
+    fname = 'ocnicepost.nml'
+    inquire (file=trim(fname), iostat=ierr)
+    if (ierr /= 0) then
+       write (6, '(3a)') 'Error: input file "', trim(fname), '" does not exist.'
+       stop
+    end if
+
+    ! Open and read namelist file.
+    open (action='read', file=trim(fname), iostat=ierr, newunit=iounit)
+    read (nml=ocnicepost_nml, iostat=ierr, unit=iounit)
+    if (ierr /= 0) then
+       write (6, '(a)') 'Error: invalid namelist format.'
+    end if
+    close (iounit)
+    nxt = srcdims(1); nyt = srcdims(2)
+    nxr = dstdims(1); nyr = dstdims(2)
+    If (srcdims(3) > 0) nlevs = srcdims(3)
+
+    ! initialize the source file type and variables
+    if (trim(ftype) == 'ocean') then
+       do_ocnpost = .true.
+       else
+       do_ocnpost = .false.
+    end if
+    input_file = trim(ftype)//'.nc'
+
+    open(newunit=logunit, file=trim(ftype)//'.post.log',form='formatted')
+    if (debug) write(logunit, '(a)')'input file: '//trim(input_file)
+
+    ! set grid names
+    fsrc = ''
+    if (nxt == 1440 .and. nyt == 1080) fsrc = 'mx025'    ! 1/4deg tripole
+    if (nxt == 720  .and. nyt == 576) fsrc = 'mx050'     ! 1/2deg tripole
+    if (nxt == 360  .and. nyt == 320) fsrc = 'mx100'     ! 1deg tripole
+    if (nxt == 72   .and. nyt == 35) fsrc = 'mx500'      ! 5deg tripole
+    if (len_trim(fsrc) == 0) then
+       write(logunit,'(a)')'source grid dimensions unknown'
+       stop
+    end if
+
+    fdst = ''
+    if (nxr == 1440 .and. nyr == 721) fdst = '0p25'      ! 1/4deg rectilinear
+    if (nxr == 720  .and. nyr == 361) fdst = '0p5'       ! 1/2 deg rectilinear
+    if (nxr == 360  .and. nyr == 181) fdst = '1p0'       ! 1 deg rectilinear
+    if (len_trim(fdst) == 0) then
+       write(logunit,'(a)')'destination grid dimensions unknown'
+       stop
+    end if
+
+    !TODO: test for consistency of source/destination resolution
+  end subroutine readnml
+
+  subroutine readcsv(nvalid)
+
+    integer, intent(out) :: nvalid
+
+    character(len= 40) :: fname
+    character(len=100) :: chead
+    character(len= 10) :: c1,c3,c4,c5,c6
+    integer :: i2
+    integer :: nn,n,ierr,iounit
+
+    ! --------------------------------------------------------
+    ! Open and read list of variables
+    ! --------------------------------------------------------
+
+    fname=trim(ftype)//'.csv'
+    open(newunit=iounit, file=trim(fname), status='old', iostat=ierr)
+    if (ierr /= 0) then
+       write (6, '(3a)') 'Error: input file "', trim(fname), '" does not exist.'
+       stop
+    end if
+
+    read(iounit,*)chead
+    nn=0
+    do n = 1,maxvars
+       read(iounit,*,iostat=ierr)c1,i2,c3,c4,c5,c6
+       if (ierr .ne. 0) exit
+       if (len_trim(c1) > 0) then
+          nn = nn+1
+          outvars(nn)%var_name = trim(c1)
+          outvars(nn)%var_dimen = i2
+          outvars(nn)%var_grid = trim(c3)
+          outvars(nn)%var_remapmethod = trim(c4)
+          outvars(nn)%var_pair = trim(c5)
+          outvars(nn)%var_pair_grid = trim(c6)
+       end if
+    end do
+    close(iounit)
+    nvalid = nn
+
+  end subroutine readcsv
+
+end module init_mod