utils_esmf_mod.F90

Go to the documentation of this file.

 
module utils_esmf_mod
 
  use esmf
  use netcdf
  use init_mod  , only : do_ocnprep, debug, logunit, vardefs, fsrc, fdst, ftype, nxr, nyr
  use utils_mod , only : dumpnc, remap
  use arrays_mod, only : maskspval, mask3d
 
  implicit none
 
  private
 
  type(esmf_routehandle) :: rh
  type(esmf_dynamicmask) :: dynamiclevmask
  type(esmf_mesh)        :: meshsrc
  type(esmf_mesh)        :: meshdst
 
  integer :: srctermprocessing = 0
 
  interface remaprh
     module procedure remaprh1d
     module procedure remaprh2d
     module procedure remaprh1ddyn
     module procedure remaprh2ddyn
  end interface remaprh
 
  interface rotremap
     module procedure rotremap2d
     module procedure rotremap3d
  end interface rotremap
 
  public :: createrh
  public :: remaprh
  public :: rotremap
  public :: chkerr
 
  character(len=*), parameter :: u_file_u = __file__   
contains
  subroutine createrh(srcmeshfile,dstmeshfile,rc)
 
    character(len=*), intent(in)  :: srcmeshfile
    character(len=*), intent(in)  :: dstmeshfile
    integer,          intent(out) :: rc
 
    ! local variables
    type(esmf_field)             :: fldsrc
    type(esmf_field)             :: flddst
    type(esmf_regridmethod_flag) :: regridmethod
    type(esmf_extrapmethod_flag) :: extrapmethod
    type(esmf_field)             :: dststatusfield
    integer, pointer             :: dststatus(:)
    character(len=20)            :: subname = 'remapRH1d'
    !----------------------------------------------------------------------------
 
    if (debug)write(logunit,'(a)')'enter '//trim(subname)
    rc = esmf_success
 
    ! use nstod for ice to maintain thermodynamically consistent values
    ! w/in layers and categories
    if (do_ocnprep) then
       regridmethod = esmf_regridmethod_bilinear
       extrapmethod = esmf_extrapmethod_nearest_idavg
    else
       regridmethod = esmf_regridmethod_nearest_stod
       extrapmethod = esmf_extrapmethod_nearest_stod
    end if
 
    meshsrc = esmf_meshcreate(filename=trim(srcmeshfile), fileformat=esmf_fileformat_esmfmesh, rc=rc)
    if (chkerr(rc,__line__,u_file_u)) call esmf_finalize(endflag=esmf_end_abort)
    fldsrc = esmf_fieldcreate(meshsrc, esmf_typekind_r8, name='mshsrc', meshloc=esmf_meshloc_element, rc=rc)
    if (chkerr(rc,__line__,u_file_u)) call esmf_finalize(endflag=esmf_end_abort)
 
    meshdst = esmf_meshcreate(filename=trim(dstmeshfile), fileformat=esmf_fileformat_esmfmesh, rc=rc)
    if (chkerr(rc,__line__,u_file_u)) call esmf_finalize(endflag=esmf_end_abort)
    flddst = esmf_fieldcreate(meshdst, esmf_typekind_r8, name='mshdst', meshloc=esmf_meshloc_element, rc=rc)
    if (chkerr(rc,__line__,u_file_u)) call esmf_finalize(endflag=esmf_end_abort)
 
    dststatusfield = esmf_fieldcreate(meshdst, esmf_typekind_i4, meshloc=esmf_meshloc_element, rc=rc)
    if (chkerr(rc,__line__,u_file_u)) return
    call esmf_fieldget(dststatusfield, farrayptr=dststatus, rc=rc)
    if (chkerr(rc,__line__,u_file_u)) return
 
    call esmf_fieldregridstore(fldsrc, flddst, routehandle=rh, &
         srcmaskvalues=(/0/), dstmaskvalues=(/0/),             &
         regridmethod=regridmethod,                            &
         extrapmethod=extrapmethod,                            &
         polemethod=esmf_polemethod_allavg,                    &
         ignoredegenerate=.true.,                              &
         srctermprocessing=srctermprocessing,                  &
         dststatusfield=dststatusfield,                        &
         unmappedaction=esmf_unmappedaction_ignore, rc=rc)
    if (chkerr(rc,__line__,u_file_u)) call esmf_finalize(endflag=esmf_end_abort)
 
    ! Create a dynamic mask object
    call esmf_dynamicmasksetr8r8r8(dynamiclevmask, dynamicsrcmaskvalue=maskspval, &
         dynamicmaskroutine=dynlevmaskproc, rc=rc)
    if (chkerr(rc,__line__,u_file_u)) return
 
    if (debug) then
       call dumpnc(trim(ftype)//'.'//trim(fdst)//'.dststat.nc', 'dststat',        &
            dims=(/nxr,nyr/), field=real(dststatus,8))
    end if
 
    if (debug)write(logunit,'(a)')'exit '//trim(subname)
  end subroutine createrh
 
  subroutine remaprh1d(kk,src_field,dst_field,rc)
 
    integer,      intent(in)  :: kk
    real(kind=8), intent(in)  :: src_field(:)
    real(kind=8), intent(out) :: dst_field(:)
    integer,      intent(out) :: rc
 
    ! local variables
    type(esmf_field)      :: fldsrc
    type(esmf_field)      :: flddst
    real(kind=8), pointer :: srcptr(:), dstptr(:)
    character(len=20)     :: subname = 'remapRH1d'
    !----------------------------------------------------------------------------
 
    if (debug)write(logunit,'(a,i5)')'enter '//trim(subname)//' ',kk
    rc = esmf_success
 
    fldsrc = esmf_fieldcreate(meshsrc, esmf_typekind_r8, meshloc=esmf_meshloc_element,rc=rc)
    if (chkerr(rc,__line__,u_file_u)) return
    flddst = esmf_fieldcreate(meshdst, esmf_typekind_r8, meshloc=esmf_meshloc_element,rc=rc)
    if (chkerr(rc,__line__,u_file_u)) return
 
    call esmf_fieldget(fldsrc, farrayptr=srcptr, rc=rc)
    if (chkerr(rc,__line__,u_file_u)) return
    call esmf_fieldget(flddst, farrayptr=dstptr, rc=rc)
    if (chkerr(rc,__line__,u_file_u)) return
 
    srcptr = src_field
    if (esmf_routehandleiscreated(rh,rc=rc)) then
       call esmf_fieldregrid(fldsrc, flddst, routehandle=rh, rc=rc)
       if (chkerr(rc,__line__,u_file_u)) return
    else
       call esmf_logwrite(trim(subname)//": RH not created ", esmf_logmsg_info)
       rc=esmf_failure
    end if
    dst_field = dstptr
 
    if (debug)write(logunit,'(a,i5)')'exit '//trim(subname)//' ',kk
  end subroutine remaprh1d
 
  subroutine remaprh2d(src_field,dst_field,rc)
 
    real(kind=8), intent(in)  :: src_field(:,:)
    real(kind=8), intent(out) :: dst_field(:,:)
    integer,      intent(out) :: rc
 
    ! local variables
    type(esmf_field)      :: fldsrc
    type(esmf_field)      :: flddst
    real(kind=8), pointer :: srcptr(:,:), dstptr(:,:)
    character(len=20)     :: subname = 'remapRH2d'
    !----------------------------------------------------------------------------
 
    if (debug)write(logunit,'(a)')'enter '//trim(subname)
    rc = esmf_success
 
    fldsrc = esmf_fieldcreate(meshsrc, esmf_typekind_r8, meshloc=esmf_meshloc_element, &
         ungriddedlbound=(/1/), ungriddedubound=(/size(src_field,1)/),                 &
         gridtofieldmap=(/2/), rc=rc)
    if (chkerr(rc,__line__,u_file_u)) return
    flddst = esmf_fieldcreate(meshdst, esmf_typekind_r8, meshloc=esmf_meshloc_element, &
         ungriddedlbound=(/1/), ungriddedubound=(/size(dst_field,1)/),                 &
         gridtofieldmap=(/2/), rc=rc)
    if (chkerr(rc,__line__,u_file_u)) return
 
    call esmf_fieldget(fldsrc, farrayptr=srcptr, rc=rc)
    if (chkerr(rc,__line__,u_file_u)) return
    call esmf_fieldget(flddst, farrayptr=dstptr, rc=rc)
    if (chkerr(rc,__line__,u_file_u)) return
 
    srcptr = src_field
    if (esmf_routehandleiscreated(rh,rc=rc)) then
       call esmf_fieldregrid(fldsrc, flddst, routehandle=rh, rc=rc)
       if (chkerr(rc,__line__,u_file_u)) return
    else
       call esmf_logwrite(trim(subname)//": RH not created ", esmf_logmsg_info)
       rc=esmf_failure
    end if
    dst_field = dstptr
 
    if (debug)write(logunit,'(a)')'exit '//trim(subname)
  end subroutine remaprh2d
 
  subroutine remaprh1ddyn(kk,src_field,dst_field,hmask,rc)
 
    !nflds,nlen
    integer,      intent(in)  :: kk
    real(kind=8), intent(in)  :: src_field(:)
    real(kind=8), intent(in)  :: hmask(:)
    real(kind=8), intent(out) :: dst_field(:)
    integer,      intent(out) :: rc
 
    ! local variables
    type(esmf_field)      :: fldsrc
    type(esmf_field)      :: flddst
    real(kind=8), pointer :: srcptr(:), dstptr(:)
    character(len=20)     :: subname = 'remapRH1ddyn'
    !----------------------------------------------------------------------------
 
    if (debug)write(logunit,'(a,i5)')'enter '//trim(subname)//' ',kk
    rc = esmf_success
 
    fldsrc = esmf_fieldcreate(meshsrc, esmf_typekind_r8, meshloc=esmf_meshloc_element)
    if (chkerr(rc,__line__,u_file_u)) return
    flddst = esmf_fieldcreate(meshdst, esmf_typekind_r8, meshloc=esmf_meshloc_element)
    if (chkerr(rc,__line__,u_file_u)) return
 
    call esmf_fieldget(fldsrc, farrayptr=srcptr, rc=rc)
    if (chkerr(rc,__line__,u_file_u)) return
    call esmf_fieldget(flddst, farrayptr=dstptr, rc=rc)
    if (chkerr(rc,__line__,u_file_u)) return
 
    srcptr = src_field
    where(hmask .eq. maskspval)srcptr = maskspval
 
    if (esmf_routehandleiscreated(rh,rc=rc)) then
       call esmf_fieldregrid(fldsrc, flddst, routehandle=rh, dynamicmask=dynamiclevmask, rc=rc)
       if (chkerr(rc,__line__,u_file_u)) return
    else
       call esmf_logwrite(trim(subname)//": RH not created ", esmf_logmsg_info)
       rc=esmf_failure
    end if
    dst_field = dstptr
 
    if (debug)write(logunit,'(a,i5)')'exit '//trim(subname)//' ',kk
  end subroutine remaprh1ddyn
 
  subroutine remaprh2ddyn(kk,src_field,dst_field,hmask,rc)
 
    integer,      intent(in)  :: kk
    real(kind=8), intent(in)  :: src_field(:,:)
    real(kind=8), intent(in)  :: hmask(:)
    real(kind=8), intent(out) :: dst_field(:,:)
    integer,      intent(out) :: rc
 
    ! local variables
    type(esmf_field)      :: fldsrc
    type(esmf_field)      :: flddst
    integer               :: i,n
    real(kind=8), pointer :: srcptr(:,:), dstptr(:,:)
    character(len=20)     :: subname = 'remapRH2ddyn'
    !----------------------------------------------------------------------------
 
    if (debug)write(logunit,'(a,i5)')'enter '//trim(subname)//' ',kk
    rc = esmf_success
 
    fldsrc = esmf_fieldcreate(meshsrc, esmf_typekind_r8, meshloc=esmf_meshloc_element, &
         ungriddedlbound=(/1/), ungriddedubound=(/size(src_field,1)/),                 &
         gridtofieldmap=(/2/), rc=rc)
    if (chkerr(rc,__line__,u_file_u)) return
    flddst = esmf_fieldcreate(meshdst, esmf_typekind_r8, meshloc=esmf_meshloc_element, &
         ungriddedlbound=(/1/), ungriddedubound=(/size(dst_field,1)/),                 &
         gridtofieldmap=(/2/), rc=rc)
    if (chkerr(rc,__line__,u_file_u)) return
 
    call esmf_fieldget(fldsrc, farrayptr=srcptr, rc=rc)
    if (chkerr(rc,__line__,u_file_u)) return
    call esmf_fieldget(flddst, farrayptr=dstptr, rc=rc)
    if (chkerr(rc,__line__,u_file_u)) return
 
    srcptr = src_field
    do n = 1,ubound(src_field,2)
       do i = 1,ubound(src_field,1)
          if(hmask(n) .eq. maskspval)srcptr(i,n) = maskspval
       end do
    end do
 
    if (esmf_routehandleiscreated(rh,rc=rc)) then
       call esmf_fieldregrid(fldsrc, flddst, routehandle=rh, dynamicmask=dynamiclevmask, rc=rc)
       if (chkerr(rc,__line__,u_file_u)) return
    else
       call esmf_logwrite(trim(subname)//": RH not created ", esmf_logmsg_info)
       rc=esmf_failure
    end if
    dst_field = dstptr
 
    if (debug)write(logunit,'(a,i5)')'exit '//trim(subname)//' ',kk
  end subroutine remaprh2ddyn
 
  subroutine rotremap2d(wdir, vars, cosrot, sinrot, dims, nflds, fields)
 
    character(len=*), intent(in)    :: wdir
    real(kind=8),     intent(in)    :: cosrot(:),sinrot(:)
    type(vardefs),    intent(in)    :: vars(:)
    integer,          intent(in)    :: dims(:)
    integer,          intent(in)    :: nflds
    real(kind=8),     intent(inout) :: fields(:,:)
 
    integer            :: n, idx1, idx2
    real(kind=8), allocatable, dimension(:) :: urot, vrot
    character(len=10)  :: vgrid1, vgrid2
    character(len=240) :: wgtsfile
    character(len=20)  :: subname = 'rotremap2d'
    !----------------------------------------------------------------------------
 
    if (debug)write(logunit,'(a)')'enter '//trim(subname)
 
    idx1 = 0; idx2 = 0
    do n = 1,nflds
       if (len_trim(vars(n)%var_pair) > 0 .and. idx1 .eq. 0) then
          idx1 = n
          idx2 = n+1
       end if
    end do
    if (idx1 .eq. 0)return
 
    vgrid1 = vars(idx1)%var_grid(1:2)
    vgrid2 = vars(idx1)%var_pair_grid(1:2)
 
    allocate(urot(1:dims(1)*dims(2))); urot = 0.0
    allocate(vrot(1:dims(1)*dims(2))); vrot = 0.0
    urot(:) = fields(idx1,:)*cosrot(:) - fields(idx2,:)*sinrot(:)
    vrot(:) = fields(idx2,:)*cosrot(:) + fields(idx1,:)*sinrot(:)
 
    wgtsfile = trim(wdir)//'tripole.'//trim(fdst)//'.Ct.to.'//trim(vgrid1)//'.bilinear.nc'
    call remap(trim(wgtsfile), urot, fields(idx1,:))
    if (debug) write(logunit,'(a)')'use '//trim(wgtsfile)//' to restagger from Ct to '//trim(vgrid1)
 
    wgtsfile = trim(wdir)//'tripole.'//trim(fdst)//'.Ct.to.'//trim(vgrid2)//'.bilinear.nc'
    call remap(trim(wgtsfile), vrot, fields(idx2,:))
    if (debug) write(logunit,'(a)')'use '//trim(wgtsfile)//' to restagger from Ct to '//trim(vgrid2)
 
    if (debug)write(logunit,'(a)')'exit '//trim(subname)
  end subroutine rotremap2d
 
  subroutine rotremap3d(wdir, vars, cosrot, sinrot, dims, nflds, fields)
 
    character(len=*), intent(in)    :: wdir
    real(kind=8),     intent(in)    :: cosrot(:),sinrot(:)
    type(vardefs),    intent(in)    :: vars(:)
    integer,          intent(in)    :: dims(:)
    integer,          intent(in)    :: nflds
    real(kind=8),     intent(inout) :: fields(:,:,:)
 
    integer            :: k, n, idx1, idx2
    real(kind=8), allocatable, dimension(:) :: urot, vrot
    character(len=10)  :: vgrid1, vgrid2
    character(len=240) :: wgtsfile
    character(len=20)  :: subname = 'rotremap3d'
    !----------------------------------------------------------------------------
 
    if (debug)write(logunit,'(a)')'enter '//trim(subname)
 
    idx1 = 0; idx2 = 0
    do n = 1,nflds
       if (len_trim(vars(n)%var_pair) > 0 .and. idx1 .eq. 0) then
          idx1 = n
          idx2 = n+1
       end if
    end do
    if (idx1 .eq. 0)return
 
    vgrid1 = vars(idx1)%var_grid(1:2)
    vgrid2 = vars(idx1)%var_pair_grid(1:2)
    if (debug) write(logunit,'(a)')'restagger from Ct to '//trim(vgrid1)//' and '//trim(vgrid2)
 
    allocate(urot(1:dims(1)*dims(2))); urot = 0.0
    allocate(vrot(1:dims(1)*dims(2))); vrot = 0.0
    do k = 1,dims(3)
       urot(:) = fields(idx1,k,:)*cosrot(:) - fields(idx2,k,:)*sinrot(:)
       vrot(:) = fields(idx2,k,:)*cosrot(:) + fields(idx1,k,:)*sinrot(:)
       wgtsfile = trim(wdir)//'tripole.'//trim(fdst)//'.Ct.to.'//trim(vgrid1)//'.bilinear.nc'
       call remap(trim(wgtsfile), urot, fields(idx1,k,:))
       wgtsfile = trim(wdir)//'tripole.'//trim(fdst)//'.Ct.to.'//trim(vgrid2)//'.bilinear.nc'
       call remap(trim(wgtsfile), vrot, fields(idx2,k,:))
    end do
 
    if (debug)write(logunit,'(a)')'exit '//trim(subname)
  end subroutine rotremap3d
 
  subroutine dynlevmaskproc(dynamicMaskList, dynamicSrcMaskValue, dynamicDstMaskValue, rc)
 
    ! input/output arguments
    type(esmf_dynamicmaskelementr8r8r8) , pointer :: dynamicMaskList(:)
    real(ESMF_KIND_R8), intent(in), optional      :: dynamicSrcMaskValue
    real(ESMF_KIND_R8), intent(in), optional      :: dynamicDstMaskValue
    integer           , intent(out)               :: rc
 
    ! local variables
    integer  :: i, j
    real(ESMF_KIND_R8) :: renorm
    character(len=20)  :: subname = 'dynLevMaskProc'
    !----------------------------------------------------------------------------
 
    if (debug)write(logunit,'(a)')'enter '//trim(subname)
    rc = esmf_success
 
    if (associated(dynamicmasklist)) then
       do i=1, size(dynamicmasklist)
          dynamicmasklist(i)%dstElement = 0.0 ! set to zero
          renorm = 0.0 ! reset
          do j = 1, size(dynamicmasklist(i)%factor)
             if (dynamicsrcmaskvalue /= dynamicmasklist(i)%srcElement(j)) then
                dynamicmasklist(i)%dstElement = dynamicmasklist(i)%dstElement + &
                     (dynamicmasklist(i)%factor(j) * dynamicmasklist(i)%srcElement(j))
                renorm = renorm + dynamicmasklist(i)%factor(j)
             endif
          enddo
          if (renorm > 0.0) then
             dynamicmasklist(i)%dstElement = dynamicmasklist(i)%dstElement / renorm
          else if (present(dynamicsrcmaskvalue)) then
             dynamicmasklist(i)%dstElement = dynamicsrcmaskvalue
          else
             rc = esmf_rc_arg_bad  ! error detected
             return
          endif
       enddo
    endif
    if (debug)write(logunit,'(a)')'exit '//trim(subname)
 
  end subroutine dynlevmaskproc
 
  logical function chkerr(rc, line, file)
    integer, intent(in) :: rc
    integer, intent(in) :: line
    character(len=*), intent(in) :: file
    integer :: lrc
    !----------------------------------------------------------------------------
    chkerr = .false.
    lrc = rc
    if (esmf_logfounderror(rctocheck=lrc, msg=esmf_logerr_passthru, line=line, file=file)) then
       chkerr = .true.
    endif
  end function chkerr
end module utils_esmf_mod