! This is a Ferret external function.  Given an input scalar or field A,
! and an input scalar or field P, computes A-FLOOR(A/P)*P.
! Uses the f90 MODULO function.
!
! In this subroutine we provide information about
! the function.  The user configurable information
! consists of the following:
!
! descr              Text description of the function
!
! num_args           Required number of arguments
!
! axis_inheritance   Type of axis for the result
!                       ( CUSTOM, IMPLIED_BY_ARGS, NORMAL, ABSTRACT )
!                       CUSTOM          - user defined axis
!                       IMPLIED_BY_ARGS - same axis as the incoming argument
!                       NORMAL          - the result is normal to this axis
!                       ABSTRACT        - an axis which only has index values
!
! piecemeal_ok       For memory optimization:
!                       axes where calculation may be performed piecemeal
!                       ( YES, NO )
!
!
! For each argument we provide the following information:
!
! name               Text name for an argument
!
! unit               Text units for an argument
!
! desc               Text description of an argument
!
! axis_influence     Are this argument's axes the same as the result grid?
!                       ( YES, NO )
!
! axis_extend       How much does Ferret need to extend arg limits relative to result
!


      subroutine modulo_init(id)
      implicit none

! use glenn carver's f90 version of ef includes

      include 'ferret_cmn/EF_Util_f90.inc'

      integer ::  id, arg

      CALL ef_version_test(ef_version)

!**********************************************************************
!                                           USER CONFIGURABLE PORTION |
!                                                                     |
!                                                                     V

      call ef_set_desc(id, 'A-FLOOR(A/P)*P')

      call ef_set_num_args(id, 2)
      call ef_set_has_vari_args(id, NO)

      call ef_set_axis_inheritance(id, IMPLIED_BY_ARGS, &
          IMPLIED_BY_ARGS,IMPLIED_BY_ARGS,IMPLIED_BY_ARGS)
      call ef_set_piecemeal_ok(id, YES, YES, YES, YES)

      arg = 1
      call ef_set_arg_name(id, arg, 'A')
      call ef_set_arg_desc(id, arg, 'Input data')
      call ef_set_axis_influence(id, arg, YES, YES, YES, YES)
      arg = 2
      call ef_set_arg_name(id, arg, 'P')
      call ef_set_arg_desc(id, arg, 'Modulo length')
      call ef_set_axis_influence(id, arg, YES, YES, YES, YES)

!                                                                     ^
!                                                                     |
!                                           USER CONFIGURABLE PORTION |
!***********************************************************************

    end subroutine modulo_init


!
! In this subroutine we compute the result
!
      subroutine modulo_compute(id, arg_1, arg_2, result)

      implicit none

      include 'ferret_cmn/EF_Util_f90.inc'
      include 'ferret_cmn/EF_mem_subsc_f90.inc'

      real :: bad_flag(EF_MAX_ARGS), bad_flag_result
      real :: arg_1(mem1lox:mem1hix, mem1loy:mem1hiy, mem1loz:mem1hiz, mem1lot:mem1hit)
      real :: arg_2(mem2lox:mem2hix, mem2loy:mem2hiy, mem2loz:mem2hiz, mem2lot:mem2hit)
      real :: result(memreslox:memreshix, memresloy:memreshiy, memresloz:memreshiz, memreslot:memreshit)

! After initialization, the 'res_' arrays contain indexing information
! for the result axes.  The 'arg_' arrays will contain the indexing
! information for each variable's axes.

      integer :: res_lo_ss(4), res_hi_ss(4), res_incr(4)
      integer :: arg_lo_ss(4,EF_MAX_ARGS), arg_hi_ss(4,EF_MAX_ARGS),&
                 arg_incr(4,EF_MAX_ARGS)

!**********************************************************************
!                                           USER CONFIGURABLE PORTION |
!                                                                     |
!                                                                     V
      integer :: id, i, j, k, l
      integer :: i1, j1, k1, l1
      integer :: i2, j2, k2, l2

      CALL ef_get_res_subscripts(id, res_lo_ss, res_hi_ss, res_incr)
      CALL ef_get_arg_subscripts(id, arg_lo_ss, arg_hi_ss, arg_incr)
      CALL ef_get_bad_flags(id, bad_flag, bad_flag_result)

      i1 = arg_lo_ss(X_AXIS,ARG1)
      i2 = arg_lo_ss(X_AXIS,ARG2)
      do i=res_lo_ss(X_AXIS), res_hi_ss(X_AXIS)

         j1 = arg_lo_ss(Y_AXIS,ARG1)
         j2 = arg_lo_ss(Y_AXIS,ARG2)
         do j=res_lo_ss(Y_AXIS), res_hi_ss(Y_AXIS)

            k1 = arg_lo_ss(Z_AXIS,ARG1)
            k2 = arg_lo_ss(Z_AXIS,ARG2)
            do k=res_lo_ss(Z_AXIS), res_hi_ss(Z_AXIS)

               l1 = arg_lo_ss(T_AXIS,ARG1)
               l2 = arg_lo_ss(T_AXIS,ARG2)
               do l=res_lo_ss(T_AXIS), res_hi_ss(T_AXIS)

                  if ( arg_1(i1,j1,k1,l1) == bad_flag(1) .or. &
                       arg_2(i2,j2,k2,l2) == bad_flag(2) .or. &
                       arg_2(i2,j2,k2,l2) == 0.) then

                     result(i,j,k,l) = bad_flag_result

                  else

                     result(i,j,k,l) = modulo(arg_1(i1,j1,k1,l1),arg_2(i2,j2,k2,l2))

                  endif

                  l1 = l1 + arg_incr(T_AXIS,ARG1)
                  l2 = l2 + arg_incr(T_AXIS,ARG2)
               enddo

               k1 = k1 + arg_incr(Z_AXIS,ARG1)
               k2 = k2 + arg_incr(Z_AXIS,ARG2)
            enddo
            j1 = j1 + arg_incr(Y_AXIS,ARG1)
            j2 = j2 + arg_incr(Y_AXIS,ARG2)
         enddo
         i1 = i1 + arg_incr(X_AXIS,ARG1)
         i2 = i2 + arg_incr(X_AXIS,ARG2)
      enddo
!                                                                     |
!                                           USER CONFIGURABLE PORTION |
!**********************************************************************

!---------------------------------------------------------------------
      end subroutine modulo_compute