mpas_increment_mod.F90

Go to the documentation of this file.

! (C) Copyright 2017 UCAR
! 
! This software is licensed under the terms of the Apache Licence Version 2.0
! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. 
 
module mpas_increment_mod
 
use atlas_module, only: atlas_geometry, atlas_field, atlas_fieldset, atlas_real
use fckit_configuration_module, only: fckit_configuration
use fckit_log_module, only: fckit_log
 
!oops
use datetime_mod
use kinds, only: kind_real
use oops_variables_mod
 
!ufo
use ufo_geovals_mod
use ufo_vars_mod
 
!MPAS-Model
use mpas_constants
use mpas_dmpar
use mpas_derived_types
use mpas_kind_types, only: strkind
use mpas_pool_routines
 
!mpas-jedi
use mpas_constants_mod
use mpas_geom_mod
use mpas_fields_mod
use mpas2ufo_vars_mod
use mpas4da_mod
 
implicit none
 
private
 
public :: diff_incr, dirac, &
        & set_atlas, to_atlas, from_atlas
 
integer, parameter    :: max_string=8000
character(max_string) :: message
 
! ------------------------------------------------------------------------------
 
contains
 
! ------------------------------------------------------------------------------
 
subroutine diff_incr(lhs,x1,x2)
 
   implicit none
   class(mpas_fields), intent(inout) :: lhs
   class(mpas_fields), intent(in)    :: x1
   class(mpas_fields), intent(in)    :: x2
   character(len=StrKIND) :: kind_op
 
   call lhs%zeros()
   if (x1%geom%nCells==x2%geom%nCells .and. x1%geom%nVertLevels==x2%geom%nVertLevels) then
     if (lhs%geom%nCells==x1%geom%nCells .and. lhs%geom%nVertLevels==x1%geom%nVertLevels) then
        kind_op = 'sub'
        call da_operator(trim(kind_op), lhs % subFields, x1 % subFields, x2 % subFields, fld_select = lhs % fldnames_ci)
     else
       call abor1_ftn("mpas_increment:diff_incr: dimension mismatch between the two variables.")
     endif
   else
     call abor1_ftn("mpas_increment:diff_incr: states not at same resolution")
   endif
 
   return
 
end subroutine diff_incr
 
! ------------------------------------------------------------------------------
 
subroutine dirac(self, f_conf)
 
   implicit none
   class(mpas_fields),        intent(inout) :: self
   type(fckit_configuration), intent(in)    :: f_conf   !< Configuration
   character(len=:), allocatable :: str
   integer                :: ndir, idir, ildir, ndirlocal
   character(len=3)       :: idirchar
   character(len=StrKIND) :: dirvar
   type (mpas_pool_iterator_type) :: poolitr
   real (kind=kind_real), dimension(:,:), pointer :: r2d_ptr_a
   real (kind=kind_real), dimension(:), pointer   :: r1d_ptr_a
   integer :: nearestcell
   integer, allocatable, dimension(:) :: dirowned, dirownedglobal
   real (kind=kind_real), allocatable, dimension(:) :: dirlats
   real (kind=kind_real), allocatable, dimension(:) :: dirlons
   integer, allocatable, dimension(:) :: dircells
 
   ! Get number and positions of Diracs
   call f_conf%get_or_die("ndir",ndir)
 
   allocate( dirowned(ndir) )
   allocate( dircells(ndir) )
 
   call f_conf%get_or_die("dirLats",dirlats)
   call f_conf%get_or_die("dirLons",dirlons)
 
   call f_conf%get_or_die("ildir",ildir)
   call f_conf%get_or_die("dirvar",str)
   dirvar = str
 
   !Test if dir is owned and find the nearest local cell
   ! (repurposed from MPAS-Release/src/core_atmosphere/diagnostics/soundings.F)
 
   ndirlocal = 0
   do idir=1,ndir
      nearestcell = self % geom % nCellsSolve
      nearestcell = nearest_cell( (dirlats(idir) * mpas_jedi_deg2rad_kr), &
                                  (dirlons(idir) * mpas_jedi_deg2rad_kr), &
                                  nearestcell, self % geom % nCells, self % geom % maxEdges, &
                                  self % geom % nEdgesOnCell, self % geom % cellsOnCell, &
                                  self % geom % latCell, self % geom % lonCell )
 
      if (nearestcell <= self % geom % nCellsSolve) then
          dirowned(idir) = 1
          dircells(idir) = nearestcell
          ndirlocal = ndirlocal + 1
      else
          dirowned(idir) = 0
          dircells(idir) = self % geom % nCells + 1
      end if
   end do
 
   write(message,*) 'This processor owns ',ndirlocal,' dirac forcing locations'
   call fckit_log%debug(message)
 
   ! Check
   if (ndir<1) call abor1_ftn("mpas_increment:dirac non-positive ndir")
 
   allocate( dirownedglobal(ndir) )
   call mpas_dmpar_max_int_array( self % geom % domain % dminfo, ndir, dirowned, dirownedglobal)
   if ( any(dirownedglobal.lt.1) ) then
         call abor1_ftn("mpas_increment:dirac invalid Lat/Lon")
   end if
 
   call mpas_dmpar_sum_int_array( self % geom % domain % dminfo, ndir, dirowned, dirownedglobal)
   if ( any(dirownedglobal.gt.1) ) then
         call abor1_ftn("mpas_increment:duplicated dirac on >1 processors")
   end if
   deallocate( dirownedglobal )
 
  if ((ildir < 1) .or. (ildir > self % geom % nVertLevels)) then
      call abor1_ftn("mpas_increment:dirac invalid ildir")
   endif
 
   ! Setup Diracs
   call self%zeros()
 
   call mpas_pool_begin_iteration(self % subFields)
 
   do while ( mpas_pool_get_next_member(self % subFields, poolitr) )
      ! Pools may in general contain dimensions, namelist options, fields, or other pools,
      ! so we select only those members of the pool that are fields
      if (poolitr % memberType == mpas_pool_field) then
         ! Fields can be integer, logical, or real. Here, we operate only on real-valued fields
         write(message,*) 'poolItr % nDims , poolItr % memberName =', poolitr % nDims , poolitr % memberName
         call fckit_log%debug(message)
 
         if (poolitr % dataType == mpas_pool_real) then
            ! Depending on the dimensionality of the field, we need to set pointers of
            ! the correct type
            if( trim(dirvar) .eq. trim(poolitr % memberName) ) then
               ndirlocal = 0
               do idir=1, ndir
                  if ( dirowned(idir).eq.1 ) then
                     if (poolitr % nDims == 1) then
                        call mpas_pool_get_array(self % subFields, trim(poolitr % memberName), r1d_ptr_a)
                        r1d_ptr_a( dircells(idir) ) = mpas_jedi_one_kr
                     else if (poolitr % nDims == 2) then
                        call mpas_pool_get_array(self % subFields, trim(poolitr % memberName), r2d_ptr_a)
                        r2d_ptr_a( ildir, dircells(idir) ) = mpas_jedi_one_kr
                     else if (poolitr % nDims == 3) then
                        call fckit_log%info ('Not implemented yet')
                     end if
                     ndirlocal = ndirlocal + 1
                  end if
               end do
            end if
         end if
      end if
   end do
 
   deallocate( dirowned )
   deallocate( dirlats )
   deallocate( dirlons )
   deallocate( dircells )
 
end subroutine dirac
 
! ------------------------------------------------------------------------------
 
!!TODO: Alternatively could make the function nearest_cell public in MPAS
!!      and then define an interface to it within this module (cleaner)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Finds the MPAS grid cell nearest to (target_lat, target_lon)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
integer function nearest_cell(target_lat, target_lon, start_cell, nCells, maxEdges, &
                              nEdgesOnCell, cellsOnCell, latCell, lonCell)
 
   implicit none
 
   real (kind=kind_real), intent(in) :: target_lat, target_lon
   integer, intent(in) :: start_cell
   integer, intent(in) :: ncells, maxedges
   integer, dimension(nCells), intent(in) :: nedgesoncell
   integer, dimension(maxEdges,nCells), intent(in) :: cellsoncell
   real (kind=kind_real), dimension(nCells), intent(in) :: latcell, loncell
 
   integer :: i
   integer :: icell
   integer :: current_cell
   real (kind=kind_real) :: current_distance, d
   real (kind=kind_real) :: nearest_distance
   type (atlas_geometry) :: ageometry
 
   nearest_cell = start_cell
   current_cell = -1
 
   ageometry = atlas_geometry("UnitSphere")
 
   do while (nearest_cell /= current_cell)
      current_cell = nearest_cell
      current_distance = ageometry%distance(loncell(current_cell), latcell(current_cell),  target_lon, target_lat)
      nearest_cell = current_cell
      nearest_distance = current_distance
      do i = 1, nedgesoncell(current_cell)
         icell = cellsoncell(i,current_cell)
         if (icell <= ncells) then
            d = ageometry%distance(loncell(icell), latcell(icell), target_lon, target_lat)
            if (d < nearest_distance) then
               nearest_cell = icell
               nearest_distance = d
            end if
         end if
      end do
   end do
end function nearest_cell
 
! ------------------------------------------------------------------------------
 
subroutine set_atlas(self, geom, vars, afieldset)
 
   implicit none
 
   type(mpas_fields),    intent(in)    :: self
   type(mpas_geom),      intent(in)    :: geom
   type(oops_variables), intent(in)    :: vars
   type(atlas_fieldset), intent(inout) :: afieldset
 
   integer :: jvar, nlevels
   logical :: var_found
   type(atlas_field) :: afield
   type(mpas_pool_iterator_type) :: poolitr
 
   do jvar = 1,vars%nvars()
      var_found = .false.
      call mpas_pool_begin_iteration(self%subFields)
      do while (mpas_pool_get_next_member(self%subFields, poolitr))
         if (trim(vars%variable(jvar))==trim(poolitr%memberName)) then
            if (.not.afieldset%has_field(vars%variable(jvar))) then
               ! Create field
               if (poolitr%nDims==1) then
                  nlevels = 0
               else if (poolitr%nDims==2) then
                  nlevels = geom%nVertLevels
               else if (poolitr%nDims==3) then
                  call abor1_ftn('not implemented yet')
               end if
               afield = geom%afunctionspace%create_field(name=vars%variable(jvar),kind=atlas_real(kind_real),levels=nlevels)
 
               ! Add field
               call afieldset%add(afield)
 
               ! Release pointer
               call afield%final()
            end if
 
            ! Set flag
            var_found = .true.
            exit
         end if
      end do
      if (.not.var_found) call abor1_ftn('variable '//trim(vars%variable(jvar))//' not found in increment')
   end do
 
end subroutine set_atlas
 
! ------------------------------------------------------------------------------
 
subroutine to_atlas(self, geom, vars, afieldset)
 
   implicit none
 
   type(mpas_fields),    intent(in)    :: self
   type(mpas_geom),      intent(in)    :: geom
   type(oops_variables), intent(in)    :: vars
   type(atlas_fieldset), intent(inout) :: afieldset
 
   integer :: jvar, nlevels
   real(kind=kind_real), pointer :: real_ptr_1(:), real_ptr_2(:,:)
   real(kind=kind_real), pointer :: r1d_ptr_a(:), r2d_ptr_a(:,:)
   logical :: var_found
   type(atlas_field) :: afield
   type(mpas_pool_iterator_type) :: poolitr
 
   do jvar = 1,vars%nvars()
      var_found = .false.
      call mpas_pool_begin_iteration(self%subFields)
      do while (mpas_pool_get_next_member(self%subFields, poolitr))
         if (trim(vars%variable(jvar))==trim(poolitr%memberName)) then
            if (afieldset%has_field(vars%variable(jvar))) then
               ! Get field
               afield = afieldset%field(vars%variable(jvar))
            else
               ! Create field
               if (poolitr%nDims==1) then
                  nlevels = 0
               else if (poolitr%nDims==2) then
                  nlevels = geom%nVertLevels
               else if (poolitr%nDims==3) then
                  call abor1_ftn('not implemented yet')
               end if
               afield = geom%afunctionspace%create_field(name=vars%variable(jvar),kind=atlas_real(kind_real),levels=nlevels)
 
               ! Add field
               call afieldset%add(afield)
            end if
 
            ! Copy data
            if (poolitr%nDims==1) then
               call afield%data(real_ptr_1)
               call mpas_pool_get_array(self%subFields, trim(poolitr%memberName), r1d_ptr_a)
               real_ptr_1 = r1d_ptr_a(1:geom%nCellsSolve)
            else if (poolitr%nDims==2) then
               call afield%data(real_ptr_2)
               call mpas_pool_get_array(self%subFields, trim(poolitr%memberName), r2d_ptr_a)
               real_ptr_2 = r2d_ptr_a(1:geom%nVertLevels,1:geom%nCellsSolve)
            end if
 
            ! Release pointer
            call afield%final()
 
            ! Set flag
            var_found = .true.
            exit
         end if
      end do
      if (.not.var_found) call abor1_ftn('variable '//trim(vars%variable(jvar))//' not found in increment')
   end do
 
end subroutine to_atlas
 
! ------------------------------------------------------------------------------
 
subroutine from_atlas(self, geom, vars, afieldset)
 
   implicit none
 
   type(mpas_fields),    intent(inout) :: self
   type(mpas_geom),      intent(in)    :: geom
   type(oops_variables), intent(in)    :: vars
   type(atlas_fieldset), intent(in)    :: afieldset
 
   integer :: jvar
   real(kind=kind_real), pointer :: real_ptr_1(:), real_ptr_2(:,:)
   real(kind=kind_real), pointer :: r1d_ptr_a(:), r2d_ptr_a(:,:)
   logical :: var_found
   type(atlas_field) :: afield
   type(mpas_pool_iterator_type) :: poolitr
 
   do jvar = 1,vars%nvars()
      var_found = .false.
      call mpas_pool_begin_iteration(self%subFields)
      do while (mpas_pool_get_next_member(self%subFields, poolitr))
         if (trim(vars%variable(jvar))==trim(poolitr%memberName)) then
            ! Get field
            afield = afieldset%field(vars%variable(jvar))
 
            ! Copy data
            if (poolitr%nDims==1) then
               call afield%data(real_ptr_1)
               call mpas_pool_get_array(self%subFields, trim(poolitr%memberName), r1d_ptr_a)
               r1d_ptr_a(1:geom%nCellsSolve) = real_ptr_1
            else if (poolitr%nDims==2) then
               call afield%data(real_ptr_2)
               call mpas_pool_get_array(self%subFields, trim(poolitr%memberName), r2d_ptr_a)
               r2d_ptr_a(1:geom%nVertLevels,1:geom%nCellsSolve) = real_ptr_2
            else if (poolitr%nDims==3) then
               call abor1_ftn('not implemented yet')
            end if
 
            ! Release pointer
            call afield%final()
 
            ! Set flag
            var_found = .true.
            exit
         end if
      end do
      if (.not.var_found) call abor1_ftn('variable '//trim(vars%variable(jvar))//' not found in increment')
   end do
 
   ! TODO: Since only local locations are updated/transferred from ug, 
   !       need MPAS HALO comms before using these fields in MPAS
 
end subroutine from_atlas
 
! ------------------------------------------------------------------------------
 
end module mpas_increment_mod