text:
! THIS VERSION: GALAHAD 2.1 - 22/03/2007 AT 09:00 GMT.

!-*-*-*-*-*-  L A N C E L O T  -B-  SCALN  M O D U L E  *-*-*-*-*-*-*-*

!  Nick Gould, for GALAHAD productions
!  Copyright reserved
!  February 1st 1995

   MODULE LANCELOT_SCALN_double

     USE LANCELOT_double, ONLY: LANCELOT_problem_type, LANCELOT_data_type

     USE LANCELOT_OTHERS_double, ONLY: OTHERS_fdgrad_save_type, OTHERS_fdgrad
     USE LANCELOT_INITW_double
     IMPLICIT NONE

     PRIVATE
     PUBLIC :: SCALN_save_type, SCALN_get_scalings

!  Set precision

     INTEGER, PARAMETER :: wp = KIND( 1.0D+0 )

!  Set other parameters

     REAL ( KIND = wp ), PARAMETER :: zero = 0.0_wp
     REAL ( KIND = wp ), PARAMETER :: point1 = 0.1_wp
     REAL ( KIND = wp ), PARAMETER :: one = 1.0_wp
     REAL ( KIND = wp ), PARAMETER :: ten = 10.0_wp

!  The SCALN_save_type derived type

     TYPE :: SCALN_save_type
       INTEGER :: nsets, igetfd, ntype, ntotel, nvrels, nnza
       LOGICAL :: altriv
       TYPE ( OTHERS_fdgrad_save_type ) :: OTHERS
     END TYPE SCALN_save_type

   CONTAINS

!-*-*-*-  L A N C E L O T  -B-  SCALN_get_scalings  S U B R O U T I N E  -*-*

     SUBROUTINE SCALN_get_scalings(                                            &
                      prob, RANGE , data, ncalc, ICALC, FT, GVALS,             &
                      FUVALS, lfuval, stopga, stopca, scaleg, scalev,          &
                      iout, iprint, buffer, fdgrad, status, S, ELFUN , GROUP )

!  To calculate suitable variable and constraint scalings for a nonlinear
!  programming problem input to LANCELOT

!  status must be set to 0 on entry, and will have the value 0 on
!  successful exit. If status is negative on exit, the user must set
!  appropriate values and re-enter. Possible values of status are:
!   -1  compute the values of the element functions and their derivatives
!       at prob%X, and return the values in FUVALS. All other arguments
!       must be unaltered on re-entry
!   -2  compute the values of the first derivatives of the the group
!       functions evaluated at FT, and return the values in GVALS( :, 2 ).
!       All other arguments must be unaltered on re-entry

!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------

     TYPE ( LANCELOT_problem_type ), INTENT( INOUT ) :: prob
     TYPE ( LANCELOT_data_type ), INTENT( INOUT ) :: data
     TYPE ( SCALN_save_type ), INTENT( INOUT ) :: S
     INTEGER, INTENT( IN ) :: lfuval, iprint, iout, buffer
     INTEGER, INTENT( INOUT ) :: status, ncalc
     REAL ( KIND = wp ), INTENT( INOUT ) :: stopca, stopga
     LOGICAL, INTENT( IN ) :: scaleg, scalev, fdgrad
     INTEGER, INTENT( INOUT ), DIMENSION( MAX( prob%ng, prob%nel ) ) :: ICALC
     REAL ( KIND = wp ), INTENT( INOUT ), DIMENSION( prob%ng ) :: FT
     REAL ( KIND = wp ), INTENT( INOUT ), DIMENSION( prob%ng, 3 ) :: GVALS
     REAL ( KIND = wp ), INTENT( INOUT ), DIMENSION( lfuval ) :: FUVALS

!-----------------------------------------------
!   I n t e r f a c e   B l o c k s
!-----------------------------------------------

     INTERFACE

!  Interface block for RANGE

       SUBROUTINE RANGE ( ielemn, transp, W1, W2, nelvar, ninvar, ieltyp,      &
                          lw1, lw2 )
       INTEGER, INTENT( IN ) :: ielemn, nelvar, ninvar, ieltyp, lw1, lw2
       LOGICAL, INTENT( IN ) :: transp
       REAL ( KIND = KIND( 1.0D+0 ) ), INTENT( IN ), DIMENSION ( lw1 ) :: W1
       REAL ( KIND = KIND( 1.0D+0 ) ), INTENT( OUT ), DIMENSION ( lw2 ) :: W2
       END SUBROUTINE RANGE

!  Interface block for ELFUN

       SUBROUTINE ELFUN ( FUVALS, XVALUE, EPVALU, ncalcf, ITYPEE, ISTAEV,      &
                          IELVAR, INTVAR, ISTADH, ISTEPA, ICALCF, ltypee,      &
                          lstaev, lelvar, lntvar, lstadh, lstepa, lcalcf,      &
                          lfuval, lxvalu, lepvlu, ifflag, ifstat )
       INTEGER, INTENT( IN ) :: ncalcf, ifflag, ltypee, lstaev, lelvar, lntvar
       INTEGER, INTENT( IN ) :: lstadh, lstepa, lcalcf, lfuval, lxvalu, lepvlu
       INTEGER, INTENT( OUT ) :: ifstat
       INTEGER, INTENT( IN ) :: ITYPEE(LTYPEE), ISTAEV(LSTAEV), IELVAR(LELVAR)
       INTEGER, INTENT( IN ) :: INTVAR(LNTVAR), ISTADH(LSTADH), ISTEPA(LSTEPA)
       INTEGER, INTENT( IN ) :: ICALCF(LCALCF)
       REAL ( KIND = KIND( 1.0D+0 ) ), INTENT( IN ) :: XVALUE(LXVALU)
       REAL ( KIND = KIND( 1.0D+0 ) ), INTENT( IN ) :: EPVALU(LEPVLU)
       REAL ( KIND = KIND( 1.0D+0 ) ), INTENT( INOUT ) :: FUVALS(LFUVAL)
       END SUBROUTINE ELFUN

!  Interface block for GROUP

       SUBROUTINE GROUP ( GVALUE, lgvalu, FVALUE, GPVALU, ncalcg,              &
                          ITYPEG, ISTGPA, ICALCG, ltypeg, lstgpa,              &
                          lcalcg, lfvalu, lgpvlu, derivs, igstat )
       INTEGER, INTENT( IN ) :: lgvalu, ncalcg
       INTEGER, INTENT( IN ) :: ltypeg, lstgpa, lcalcg, lfvalu, lgpvlu
       INTEGER, INTENT( OUT ) :: igstat
       LOGICAL, INTENT( IN ) :: derivs
       INTEGER, INTENT( IN ), DIMENSION ( ltypeg ) :: ITYPEG
       INTEGER, INTENT( IN ), DIMENSION ( lstgpa ) :: ISTGPA
       INTEGER, INTENT( IN ), DIMENSION ( lcalcg ) :: ICALCG
       REAL ( KIND = KIND( 1.0D+0 ) ), INTENT( IN ),                           &
                                       DIMENSION ( lfvalu ) :: FVALUE
       REAL ( KIND = KIND( 1.0D+0 ) ), INTENT( IN ),                           &
                                       DIMENSION ( lgpvlu ) :: GPVALU
       REAL ( KIND = KIND( 1.0D+0 ) ), INTENT( INOUT ),                        &
                                       DIMENSION ( lgvalu, 3 ) :: GVALUE
       END SUBROUTINE GROUP

     END INTERFACE

!-----------------------------------------------------
!   O p t i o n a l   D u m m y   A r g u m e n t s
!-----------------------------------------------------

     OPTIONAL :: ELFUN , GROUP

!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------

     INTEGER, PARAMETER :: lmin = 10000
     LOGICAL :: internal_el, internal_gr

     internal_el = PRESENT( ELFUN )
     internal_gr = PRESENT( GROUP )

! Initial entry: set up data

     IF ( status == 0 ) THEN

       IF ( prob%n <= 0 .OR. prob%ng <= 0 ) RETURN

!  If the element functions are to be evaluated internally, check that
!  the user has supplied appropriate information

       IF ( internal_el ) THEN
         IF ( ALLOCATED( prob%ISTEPA ) .AND. ALLOCATED( prob%ITYPEE ) .AND.  &
              ALLOCATED( prob%EPVALU ) ) THEN
           IF ( SIZE( prob%ISTEPA ) < prob%nel + 1 .OR.                        &
                SIZE( prob%ITYPEE ) < prob%nel ) THEN
             status = 10 ; RETURN ; END IF
           IF ( SIZE( prob%EPVALU ) < prob%ISTEPA( prob%nel + 1 ) - 1) THEN
             status = 10 ; RETURN ; END IF
         ELSE
           status = 10 ; RETURN
         END IF
       END IF

!  Do the same if the group functions are to be evaluated internally.

       IF ( internal_gr ) THEN
         IF ( ALLOCATED( prob%ISTGPA ) .AND. ALLOCATED( prob%ITYPEG ) .AND.  &
              ALLOCATED( prob%GPVALU ) ) THEN
           IF ( SIZE( prob%ISTGPA ) < prob%ng + 1 .OR.                         &
                SIZE( prob%ITYPEG ) < prob%ng ) THEN
             status = 11 ; RETURN ; END IF
           IF ( SIZE( prob%GPVALU ) < prob%ISTGPA( prob%ng + 1 ) - 1 ) THEN
             status = 11 ; RETURN ; END IF
         ELSE
           status = 11 ; RETURN
         END IF
       END IF

!  Set initial array lengths for EXTEND arrays

       data%S%EXTEND%lirnh = lmin
       data%S%EXTEND%ljcnh = lmin
       data%S%EXTEND%llink_min = lmin
       data%S%EXTEND%lirnh_min = lmin
       data%S%EXTEND%ljcnh_min = lmin
       data%S%EXTEND%lh_min = lmin
       data%S%EXTEND%lwtran_min = lmin
       data%S%EXTEND%litran_min = lmin
       data%S%EXTEND%lh = lmin

     END IF

!  ========================================================
!  Call the main subroutine to perform the bulk of the work
!  ========================================================

!  Both internal element and group evaluations will be performed
!  -------------------------------------------------------------

     IF ( internal_el .AND. internal_gr ) THEN
       CALL SCALN_get_scalings_main(                                           &
           prob%n, prob%ng, prob%nel, lfuval,                                  &
           prob%X , prob%Y, ncalc, ICALC, prob%ICNA,                           &
           prob%ISTADA, prob%IELING, prob%INTVAR, prob%ISTADG, prob%ISTAEV,    &
           prob%IELVAR, prob%ISTADH, prob%A, prob%B, FT, GVALS, FUVALS,        &
           stopga, stopca, prob%GSCALE, prob%ESCALE, prob%VSCALE, scaleg,      &
           scalev, prob%KNDOFG, prob%ITYPEE, prob%GXEQX, prob%INTREP, RANGE ,  &
           iout, iprint, buffer, fdgrad, status, S,                            &
!    workspace
           data%ISYMMD, data%ISYMMH, data%ISWKSP, data%ISTAJC, data%ISTAGV,    &
           data%ISVGRP, data%ISLGRP, data%IGCOLJ, data%IVALJR, data%IUSED ,    &
           data%ITYPER, data%ISSWTR, data%ISSITR, data%ISET  , data%ISVSET,    &
           data%INVSET, data%LIST_elements, data%FUVALS_temp, data%X0,         &
           data%IW_asmbl, data%NZ_comp_w, data%W_ws, data%W_el,                &
           data%W_in, data%H_el, data%H_in, data%WTRANS,                       &
           data%S%EXTEND%lwtran, data%ITRANS, data%S%EXTEND%litran,            &
           data%S%EXTEND%lwtran_min, data%S%EXTEND%litran_min,                 &
           data%LINK_elem_uses_var,                                            &
           data%S%EXTEND%l_link_e_u_v, data%S%EXTEND%llink_min,                &
           ELFUN = ELFUN , ISTEPA = prob%ISTEPA, EPVALU = prob%EPVALU,         &
           GROUP = GROUP , ISTGPA = prob%ISTGPA,                               &
           ITYPEG = prob%ITYPEG, GPVALU = prob%GPVALU )

!  Just internal element evaluations will be performed
!  ---------------------------------------------------

     ELSE IF ( internal_el ) THEN
       CALL SCALN_get_scalings_main(                                           &
           prob%n, prob%ng, prob%nel, lfuval,                                  &
           prob%X , prob%Y, ncalc, ICALC, prob%ICNA,                           &
           prob%ISTADA, prob%IELING, prob%INTVAR, prob%ISTADG, prob%ISTAEV,    &
           prob%IELVAR, prob%ISTADH, prob%A, prob%B, FT, GVALS, FUVALS,        &
           stopga, stopca, prob%GSCALE, prob%ESCALE, prob%VSCALE, scaleg,      &
           scalev, prob%KNDOFG, prob%ITYPEE, prob%GXEQX, prob%INTREP, RANGE ,  &
           iout, iprint, buffer, fdgrad, status, S,                            &
!    workspace
           data%ISYMMD, data%ISYMMH, data%ISWKSP, data%ISTAJC, data%ISTAGV,    &
           data%ISVGRP, data%ISLGRP, data%IGCOLJ, data%IVALJR, data%IUSED ,    &
           data%ITYPER, data%ISSWTR, data%ISSITR, data%ISET  , data%ISVSET,    &
           data%INVSET, data%LIST_elements, data%FUVALS_temp, data%X0,         &
           data%IW_asmbl, data%NZ_comp_w, data%W_ws, data%W_el,                &
           data%W_in, data%H_el, data%H_in, data%WTRANS,                       &
           data%S%EXTEND%lwtran, data%ITRANS, data%S%EXTEND%litran,            &
           data%S%EXTEND%lwtran_min, data%S%EXTEND%litran_min,                 &
           data%LINK_elem_uses_var,                                            &
           data%S%EXTEND%l_link_e_u_v, data%S%EXTEND%llink_min,                &
           ELFUN = ELFUN , ISTEPA = prob%ISTEPA, EPVALU = prob%EPVALU )

!  Just internal group evaluations will be performed
!  -------------------------------------------------

     ELSE IF ( internal_gr ) THEN
       CALL SCALN_get_scalings_main(                                           &
           prob%n, prob%ng, prob%nel, lfuval,                                  &
           prob%X , prob%Y, ncalc, ICALC, prob%ICNA,                           &
           prob%ISTADA, prob%IELING, prob%INTVAR, prob%ISTADG, prob%ISTAEV,    &
           prob%IELVAR, prob%ISTADH, prob%A, prob%B, FT, GVALS, FUVALS,        &
           stopga, stopca, prob%GSCALE, prob%ESCALE, prob%VSCALE, scaleg,      &
           scalev, prob%KNDOFG, prob%ITYPEE, prob%GXEQX, prob%INTREP, RANGE ,  &
           iout, iprint, buffer, fdgrad, status, S,                            &
!    workspace
           data%ISYMMD, data%ISYMMH, data%ISWKSP, data%ISTAJC, data%ISTAGV,    &
           data%ISVGRP, data%ISLGRP, data%IGCOLJ, data%IVALJR, data%IUSED ,    &
           data%ITYPER, data%ISSWTR, data%ISSITR, data%ISET  , data%ISVSET,    &
           data%INVSET, data%LIST_elements, data%FUVALS_temp, data%X0,         &
           data%IW_asmbl, data%NZ_comp_w, data%W_ws, data%W_el,                &
           data%W_in, data%H_el, data%H_in, data%WTRANS,                       &
           data%S%EXTEND%lwtran, data%ITRANS, data%S%EXTEND%litran,            &
           data%S%EXTEND%lwtran_min, data%S%EXTEND%litran_min,                 &
           data%LINK_elem_uses_var,                                            &
           data%S%EXTEND%l_link_e_u_v, data%S%EXTEND%llink_min,                &
           GROUP = GROUP, ISTGPA = prob%ISTGPA,                                &
           ITYPEG = prob%ITYPEG, GPVALU = prob%GPVALU )

!  Element and group evaluations will be performed via reverse communication
!  -------------------------------------------------------------------------

     ELSE
       CALL SCALN_get_scalings_main(                                           &
           prob%n, prob%ng, prob%nel, lfuval,                                  &
           prob%X , prob%Y, ncalc, ICALC, prob%ICNA,                           &
           prob%ISTADA, prob%IELING, prob%INTVAR, prob%ISTADG, prob%ISTAEV,    &
           prob%IELVAR, prob%ISTADH, prob%A, prob%B, FT, GVALS, FUVALS,        &
           stopga, stopca, prob%GSCALE, prob%ESCALE, prob%VSCALE, scaleg,      &
           scalev, prob%KNDOFG, prob%ITYPEE, prob%GXEQX, prob%INTREP, RANGE ,  &
           iout, iprint, buffer, fdgrad, status, S,                            &
!    workspace
           data%ISYMMD, data%ISYMMH, data%ISWKSP, data%ISTAJC, data%ISTAGV,    &
           data%ISVGRP, data%ISLGRP, data%IGCOLJ, data%IVALJR, data%IUSED ,    &
           data%ITYPER, data%ISSWTR, data%ISSITR, data%ISET  , data%ISVSET,    &
           data%INVSET, data%LIST_elements, data%FUVALS_temp, data%X0,         &
           data%IW_asmbl, data%NZ_comp_w, data%W_ws, data%W_el,                &
           data%W_in, data%H_el, data%H_in, data%WTRANS,                       &
           data%S%EXTEND%lwtran, data%ITRANS, data%S%EXTEND%litran,            &
           data%S%EXTEND%lwtran_min, data%S%EXTEND%litran_min,                 &
           data%LINK_elem_uses_var,                                            &
           data%S%EXTEND%l_link_e_u_v, data%S%EXTEND%llink_min )
     END IF

     RETURN

!  End of subroutine SCALN_get_scalings

     END SUBROUTINE SCALN_get_scalings

!-*-*-  L A N C E L O T  -B-  SCALN_get_scalings_main  S U B R O U T I N E  -*-*

     SUBROUTINE SCALN_get_scalings_main(                                       &
                      n, ng, nel, lfuval,                                      &
                      X , Y , ncalc , ICALC , ICNA  , ISTADA,                  &
                      IELING, INTVAR, ISTADG, ISTAEV, IELVAR, ISTADH,          &
                      A , B , FT    , GVALS , FUVALS, stopga, stopca,          &
                      GSCALE, ESCALE, VSCALE, scaleg, scalev,                  &
                      KNDOFG, ITYPEE, GXEQX , INTREP, RANGE , iout  , iprint,  &
                      buffer, fdgrad, status, S,                               &
!  workspace
                      ISYMMD, ISYMMH, ISWKSP, ISTAJC, ISTAGV, ISVGRP,          &
                      ISLGRP, IGCOLJ, IVALJR, IUSED , ITYPER, ISSWTR,          &
                      ISSITR, ISET  , ISVSET, INVSET, LIST_elements,           &
                      FUVALS_temp, X_temp,                                     &
                      IW_asmbl, NZ_comp_w, W_ws, W_el, W_in, H_el, H_in,       &
                      WTRANS, lwtran, ITRANS, litran,                          &
                      lwtran_min, litran_min, LINK_elem_uses_var  ,            &
                      l_link_e_u_v, llink_min,                                 &
                      ELFUN , ISTEPA, EPVALU, GROUP , ISTGPA, ITYPEG, GPVALU )

!  To calculate suitable variable and constraint scalings for a nonlinear
!  programming problem input to LANCELOT

!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------

     INTEGER, INTENT( IN ) :: n, ng, nel, lfuval
     TYPE ( SCALN_save_type ), INTENT( INOUT ) :: S
     INTEGER, INTENT( IN ) :: iprint, iout, buffer
     INTEGER, INTENT( INOUT ) :: status, ncalc
     REAL ( KIND = wp ), INTENT( INOUT ) :: stopca, stopga
     LOGICAL, INTENT( IN ) :: scaleg, scalev, fdgrad
     INTEGER, INTENT( IN ), DIMENSION( ng      ) :: KNDOFG
     INTEGER, INTENT( IN ), DIMENSION( nel ) :: ITYPEE
     INTEGER, INTENT( IN ), DIMENSION( ng  + 1 ) :: ISTADA, ISTADG
     INTEGER, INTENT( INOUT ), DIMENSION( nel + 1 ) :: ISTAEV, INTVAR
     INTEGER, INTENT( IN ), DIMENSION( ISTADG( ng + 1 ) - 1 ) :: IELING
     INTEGER, INTENT( IN ), DIMENSION( ISTAEV( nel + 1 ) - 1 ) :: IELVAR
     INTEGER, INTENT( IN ), DIMENSION( ISTADA( ng  + 1 ) - 1 ) :: ICNA
     INTEGER, INTENT( OUT ), DIMENSION( nel + 1 ) :: ISTADH
     INTEGER, INTENT( INOUT ),                                                 &
              DIMENSION( MAX( ng, nel ) ) :: ICALC
     REAL ( KIND = wp ), INTENT( INOUT ), DIMENSION( n ) :: X
     REAL ( KIND = wp ), INTENT( IN ),                                         &
            DIMENSION( ISTADA( ng + 1 ) - 1  ) :: A
     REAL ( KIND = wp ), INTENT( IN  ), DIMENSION( ng ) :: B
     REAL ( KIND = wp ), INTENT( IN  ),                                        &
            DIMENSION( ISTADG( ng + 1 ) - 1  ) :: ESCALE
     REAL ( KIND = wp ), INTENT( INOUT ), DIMENSION( n ) :: VSCALE
     REAL ( KIND = wp ), INTENT( INOUT ), DIMENSION( ng ) :: FT, Y, GSCALE
     REAL ( KIND = wp ), INTENT( INOUT ), DIMENSION( ng, 3 ) :: GVALS
     REAL ( KIND = wp ), INTENT( INOUT ),                                      &
            DIMENSION( lfuval ) :: FUVALS
     LOGICAL, INTENT( IN    ), DIMENSION( ng  ) :: GXEQX
     LOGICAL, INTENT( IN    ), DIMENSION( nel ) :: INTREP

!--------------------------------------------------------------
!   D u m m y   A r g u m e n t s  f o r   W o r k s p a c e
!--------------------------------------------------------------

     INTEGER, INTENT( INOUT ) :: lwtran, litran, lwtran_min, litran_min
     INTEGER, INTENT( INOUT ) :: l_link_e_u_v, llink_min
     INTEGER, ALLOCATABLE, DIMENSION( : , : ) :: ISYMMH
     INTEGER, ALLOCATABLE, DIMENSION( : ) :: ISYMMD, ISWKSP, ISTAJC
     INTEGER, ALLOCATABLE, DIMENSION( : ) :: ISTAGV, ISVGRP, ISLGRP, IGCOLJ
     INTEGER, ALLOCATABLE, DIMENSION( : ) :: IVALJR, IUSED , ITYPER, ISSWTR
     INTEGER, ALLOCATABLE, DIMENSION( : ) :: ISSITR, ISET  , ISVSET, INVSET
     INTEGER, ALLOCATABLE, DIMENSION( : ) :: LIST_elements
     INTEGER, ALLOCATABLE, DIMENSION( : ) :: IW_asmbl
     INTEGER, ALLOCATABLE, DIMENSION( : ) :: NZ_comp_w
     REAL ( KIND = wp ), ALLOCATABLE, DIMENSION( : ) :: W_ws
     REAL ( KIND = wp ), ALLOCATABLE, DIMENSION( : ) :: W_el
     REAL ( KIND = wp ), ALLOCATABLE, DIMENSION( : ) :: W_in
     REAL ( KIND = wp ), ALLOCATABLE, DIMENSION( : ) :: H_el
     REAL ( KIND = wp ), ALLOCATABLE, DIMENSION( : ) :: H_in
     INTEGER, ALLOCATABLE, DIMENSION( : ) :: ITRANS, LINK_elem_uses_var
     REAL ( KIND = wp ), ALLOCATABLE, DIMENSION( : ) :: FUVALS_temp, X_temp
     REAL ( KIND = wp ), ALLOCATABLE, DIMENSION( : ) :: WTRANS

!-----------------------------------------------
!   I n t e r f a c e   B l o c k s
!-----------------------------------------------

     INTERFACE

       SUBROUTINE RANGE ( ielemn, transp, W1, W2, nelvar, ninvar, ieltyp,      &
                         lw1, lw2 )
       INTEGER, INTENT( IN ) :: ielemn, nelvar, ninvar, ieltyp, lw1, lw2
       LOGICAL, INTENT( IN ) :: transp
       REAL ( KIND = KIND( 1.0D+0 ) ), INTENT( IN ), DIMENSION ( lw1 ) :: W1
       REAL ( KIND = KIND( 1.0D+0 ) ), INTENT( OUT ), DIMENSION ( lw2 ) :: W2
       END SUBROUTINE RANGE

!  Interface block for ELFUN

       SUBROUTINE ELFUN ( FUVALS, XVALUE, EPVALU, ncalcf, ITYPEE, ISTAEV,      &
                          IELVAR, INTVAR, ISTADH, ISTEPA, ICALCF, ltypee,      &
                          lstaev, lelvar, lntvar, lstadh, lstepa, lcalcf,      &
                          lfuval, lxvalu, lepvlu, ifflag, ifstat )
       INTEGER, INTENT( IN ) :: ncalcf, ifflag, ltypee, lstaev, lelvar, lntvar
       INTEGER, INTENT( IN ) :: lstadh, lstepa, lcalcf, lfuval, lxvalu, lepvlu
       INTEGER, INTENT( OUT ) :: ifstat
       INTEGER, INTENT( IN ) :: ITYPEE(LTYPEE), ISTAEV(LSTAEV), IELVAR(LELVAR)
       INTEGER, INTENT( IN ) :: INTVAR(LNTVAR), ISTADH(LSTADH), ISTEPA(LSTEPA)
       INTEGER, INTENT( IN ) :: ICALCF(LCALCF)
       REAL ( KIND = KIND( 1.0D+0 ) ), INTENT( IN ) :: XVALUE(LXVALU)
       REAL ( KIND = KIND( 1.0D+0 ) ), INTENT( IN ) :: EPVALU(LEPVLU)
       REAL ( KIND = KIND( 1.0D+0 ) ), INTENT( INOUT ) :: FUVALS(LFUVAL)
       END SUBROUTINE ELFUN

!  Interface block for GROUP

       SUBROUTINE GROUP ( GVALUE, lgvalu, FVALUE, GPVALU, ncalcg,              &
                          ITYPEG, ISTGPA, ICALCG, ltypeg, lstgpa,              &
                          lcalcg, lfvalu, lgpvlu, derivs, igstat )
       INTEGER, INTENT( IN ) :: lgvalu, ncalcg
       INTEGER, INTENT( IN ) :: ltypeg, lstgpa, lcalcg, lfvalu, lgpvlu
       INTEGER, INTENT( OUT ) :: igstat
       LOGICAL, INTENT( IN ) :: derivs
       INTEGER, INTENT( IN ), DIMENSION ( ltypeg ) :: ITYPEG
       INTEGER, INTENT( IN ), DIMENSION ( lstgpa ) :: ISTGPA
       INTEGER, INTENT( IN ), DIMENSION ( lcalcg ) :: ICALCG
       REAL ( KIND = KIND( 1.0D+0 ) ), INTENT( IN ),                           &
                                       DIMENSION ( lfvalu ) :: FVALUE
       REAL ( KIND = KIND( 1.0D+0 ) ), INTENT( IN ),                           &
                                       DIMENSION ( lgpvlu ) :: GPVALU
       REAL ( KIND = KIND( 1.0D+0 ) ), INTENT( INOUT ),                        &
                                       DIMENSION ( lgvalu, 3 ) :: GVALUE
       END SUBROUTINE GROUP

     END INTERFACE

!-----------------------------------------------------
!   O p t i o n a l   D u m m y   A r g u m e n t s
!-----------------------------------------------------

     INTEGER, INTENT( IN ), OPTIONAL, DIMENSION( nel + 1 ) :: ISTEPA
     INTEGER, INTENT( IN ), OPTIONAL, DIMENSION( ng + 1 ) :: ISTGPA
     INTEGER, INTENT( IN ), OPTIONAL, DIMENSION( ng ) :: ITYPEG
     REAL ( KIND = wp ), INTENT( IN ), OPTIONAL, DIMENSION( : ) :: EPVALU
     REAL ( KIND = wp ), INTENT( IN ), OPTIONAL, DIMENSION( : ) :: GPVALU
     OPTIONAL :: ELFUN , GROUP

!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------

     INTEGER :: i , ig, nnzgrj, lstagv, lendgv, iestev, nvarel, iel, l
     INTEGER :: ig1, j, ii, k , ieintv, nvargp, lfxi  , lgxi
     INTEGER :: lhxi  , lggfx , ldx   , lgrjac, maxsel, nin
     INTEGER :: ntotin, lnguvl, lnhuvl, alloc_status
     REAL ( KIND = wp ) :: gi    , scale , scalee, ftt   , grpmax
     REAL ( KIND = wp ) :: varmin, tol
     LOGICAL :: centrl, reallocate, external_el, external_gr
     CHARACTER ( LEN = 24 ) :: bad_alloc

!-----------------------------------------------
!   A l l o c a t a b l e   A r r a y s
!-----------------------------------------------

     INTEGER, ALLOCATABLE, DIMENSION( : ) :: IRNGRJ, ICNGRJ
     REAL ( KIND = wp ), ALLOCATABLE, DIMENSION( : ) :: GRPSCA, GRJAC
     REAL ( KIND = wp ), ALLOCATABLE, DIMENSION( : ) :: WK_n, WK_t

     external_el = .NOT. PRESENT( ELFUN )
     external_gr = .NOT. PRESENT( GROUP )

!  Branch to the interior of the code if a re-entry is being made

     SELECT CASE ( status )
     CASE ( - 1 ) ; GO TO 110
     CASE ( - 2 ) ; GO TO 160
     END SELECT

!  0. Initialize data

!  Ensure that all element functions are evaluated

     DO i = 1, nel ; ICALC( i ) = i ; END DO

!  If finite-difference gradients are required, set up the necessary
!  data structures (this is a bit wasteful as not all the partitions
!  are used)

     S%ntotel = ISTADG( ng  + 1 ) - 1
     S%nvrels = ISTAEV( nel + 1 ) - 1
     S%nnza   = ISTADA( ng  + 1 ) - 1
     CALL INITW_initialize_workspace(                                          &
         n, ng, nel, S%ntotel, S%nvrels, S%nnza  , n,  &
         nvargp, IELING, ISTADG, IELVAR, ISTAEV, INTVAR, ISTADH, ICNA, ISTADA, &
         ITYPEE, GXEQX , INTREP, S%altriv,.FALSE., fdgrad, lfxi, lgxi, lhxi  , &
         lggfx , ldx   , lnguvl, lnhuvl, ntotin, S%ntype ,                     &
         S%nsets , maxsel, RANGE , iprint, iout, buffer,                       &
! workspace
         lwtran, litran, lwtran_min, litran_min,  l_link_e_u_v, llink_min,     &
         ITRANS, LINK_elem_uses_var, WTRANS, ISYMMD, ISWKSP, ISTAJC, ISTAGV,   &
         ISVGRP, ISLGRP, IGCOLJ, IVALJR, IUSED , ITYPER, ISSWTR, ISSITR,       &
         ISET  , ISVSET, INVSET, LIST_elements, ISYMMH, IW_asmbl, NZ_comp_w,   &
         W_ws, W_el, W_in, H_el, H_in, status, alloc_status, bad_alloc,        &
         .FALSE. )

     IF ( status > 0 ) THEN ; WRITE( iout, 2100 ) ; GO TO 610 ; END IF

     S%igetfd = 0

!  Set up the starting addresses for the element gradients with respect to
!  their internal variables

!  1. Calculate the required element and group derivatives

!  Return to the calling program to obtain the element function and derivative
!  values at the initial point

     ncalc = nel
     status = - 1
     IF ( external_el ) RETURN

!  Evaluate the element function and derivative value

 100 CONTINUE
     CALL ELFUN ( FUVALS, X, EPVALU, ncalc, ITYPEE, ISTAEV, IELVAR, INTVAR,    &
                  ISTADH, ISTEPA, ICALC, nel, nel + 1, ISTAEV( nel + 1 ) - 1,  &
                  nel + 1, nel + 1, nel + 1, nel, lfuval, n,                   &
                  ISTEPA( nel + 1 ) - 1, 1, i )
     IF ( .NOT. fdgrad )                                                       &
       CALL ELFUN ( FUVALS, X, EPVALU, ncalc, ITYPEE, ISTAEV, IELVAR, INTVAR,  &
                    ISTADH, ISTEPA, ICALC, nel, nel + 1, ISTAEV( nel + 1 ) - 1,&
                    nel + 1, nel + 1, nel + 1, nel, lfuval, n,                 &
                    ISTEPA( nel + 1 ) - 1, 2, i )

 110 CONTINUE

!  If finite-difference gradients are used, compute their values

     IF ( fdgrad .AND. nel > 0 ) THEN

!  Store the values of the nonlinear elements for future use

       IF ( S%igetfd == 0 ) THEN

          reallocate = .TRUE.
          IF ( ALLOCATED( X_temp ) ) THEN
            IF ( SIZE( X_temp ) < n ) THEN ; DEALLOCATE( X_temp )
            ELSE ; reallocate = .FALSE. ; END IF
          END IF
          IF ( reallocate ) THEN
            ALLOCATE( X_temp( n ), STAT = alloc_status )
            IF ( alloc_status /= 0 ) THEN ; bad_alloc = 'X_temp'; GO TO 620
            END IF
          END IF

          reallocate = .TRUE.
          IF ( ALLOCATED( FUVALS_temp ) ) THEN
            IF ( SIZE( FUVALS_temp ) < nel ) THEN
               DEALLOCATE( FUVALS_temp )
            ELSE ; reallocate = .FALSE. ; END IF
          END IF
          IF ( reallocate ) THEN
            ALLOCATE( FUVALS_temp( nel ), STAT = alloc_status )
            IF ( alloc_status /= 0 ) THEN ; bad_alloc = 'FUVAL_'; GO TO 620
            END IF
          END IF

          X_temp( : n ) = X( : n )
          FUVALS_temp( : nel ) = FUVALS( : nel )
          centrl = .TRUE.
       END IF

!  Obtain a further set of differences

       CALL OTHERS_fdgrad( n, nel, lfuval,                                     &
                           S%ntotel, S%nvrels, S%nsets , IELVAR,               &
                           ISTAEV, IELING, ICALC, ncalc, INTVAR,               &
                           S%ntype, X_temp, X, FUVALS, centrl, S%igetfd,       &
                           S%OTHERS, ISVSET, ISET, INVSET, ISSWTR,             &
                           ISSITR, ITYPER, LIST_elements,                      &
                           LINK_elem_uses_var, WTRANS, ITRANS )
       IF ( S%igetfd > 0 ) THEN
         IF ( external_el ) THEN ; RETURN ; ELSE ; GO TO 100 ; END IF
       END IF

!  Restore the values of the nonlinear elements at X

       S%igetfd = S%nsets + 1
       FUVALS( : nel ) = FUVALS_temp( : nel )
     END IF

!  Compute the group argument values FT

     DO ig = 1, ng
       ftt = - B( ig )

!  Include the contribution from the linear element

       DO j = ISTADA( ig ), ISTADA( ig + 1 ) - 1
         IF ( ICNA( j ) <= n ) ftt = ftt + A( j ) * X( ICNA( j ) )
       END DO

!  Include the contributions from the nonlinear elements

!      ftt = ftt + SUM( ESCALE( ISTADG( ig ) : ISTADG( ig + 1 ) - 1 )          &
!            * FUVALS( IELING( ISTADG( ig ) : ISTADG( ig + 1 ) - 1 ) ) )
       DO j = ISTADG( ig ), ISTADG( ig + 1 ) - 1
         ftt = ftt + ESCALE( j ) * FUVALS( IELING( j ) )
       END DO
       FT( ig ) = ftt
       ICALC( ig ) = ig
     END DO

!  Compute the group function values

     IF ( S%altriv ) THEN
       GVALS( : ng, 2 ) = FT( : ng )
     ELSE

!  If necessary, return to the calling program to obtain the group function
!  and derivative values at the initial point

       ncalc = ng
       status = - 2
       IF ( external_gr ) RETURN

!  Evaluate the group function derivatives

       CALL GROUP ( GVALS, ng, FT, GPVALU, ncalc, ITYPEG, ISTGPA, ICALC,       &
                    ng, ng + 1, ng, ng, ISTGPA( ng + 1 ) - 1, .TRUE., i )
     END IF

 160 CONTINUE
     IF ( .NOT. S%altriv ) THEN
       WHERE ( GXEQX( : ng ) ) GVALS( : ng, 2 ) = FT( : ng)
     END IF

!  2. Form the Jacobian matrix of the nonlinear function (group( 1 ), .... ,
!     group(NG))(transpose).

     lgrjac = ISTAGV( ng + 1 ) - 1

!    Allocate space to hold the Jacobian in co-ordinate form

     reallocate = .TRUE.
     IF ( ALLOCATED( IRNGRJ ) ) THEN
       IF ( SIZE( IRNGRJ ) < lgrjac ) THEN ; DEALLOCATE( IRNGRJ )
       ELSE ; reallocate = .FALSE. ; END IF
     END IF
     IF ( reallocate ) THEN
       ALLOCATE( IRNGRJ( lgrjac ), STAT = alloc_status )
       IF ( alloc_status /= 0 ) THEN ; bad_alloc = 'IRNGRJ'; GO TO 620 ; END IF
     END IF

     reallocate = .TRUE.
     IF ( ALLOCATED( ICNGRJ ) ) THEN
       IF ( SIZE( ICNGRJ ) < lgrjac ) THEN ; DEALLOCATE( ICNGRJ )
       ELSE ; reallocate = .FALSE. ; END IF
     END IF
     IF ( reallocate ) THEN
       ALLOCATE( ICNGRJ( lgrjac ), STAT = alloc_status )
       IF ( alloc_status /= 0 ) THEN ; bad_alloc = 'ICNGRJ'; GO TO 620 ; END IF
     END IF

     reallocate = .TRUE.
     IF ( ALLOCATED( GRJAC ) ) THEN
       IF ( SIZE( GRJAC ) < lgrjac ) THEN ; DEALLOCATE( GRJAC )
       ELSE ; reallocate = .FALSE. ; END IF
     END IF
     IF ( reallocate ) THEN
       ALLOCATE( GRJAC( lgrjac ), STAT = alloc_status )
       IF ( alloc_status /= 0 ) THEN ; bad_alloc = 'GRJAC'; GO TO 620 ; END IF
     END IF

     reallocate = .TRUE.
     IF ( ALLOCATED( GRPSCA ) ) THEN
       IF ( SIZE( GRPSCA ) < ng ) THEN ; DEALLOCATE( GRPSCA )
       ELSE ; reallocate = .FALSE. ; END IF
     END IF
     IF ( reallocate ) THEN
       ALLOCATE( GRPSCA( ng ), STAT = alloc_status )
       IF ( alloc_status /= 0 ) THEN ; bad_alloc = 'GRPSCA'; GO TO 620 ; END IF
     END IF

     reallocate = .TRUE.
     IF ( ALLOCATED( WK_n ) ) THEN
       IF ( SIZE( WK_n ) < n ) THEN ; DEALLOCATE( WK_n )
       ELSE ; reallocate = .FALSE. ; END IF
     END IF
     IF ( reallocate ) THEN
       ALLOCATE( WK_n( n ), STAT = alloc_status )
       IF ( alloc_status /= 0 ) THEN ; bad_alloc = 'WK_n'; GO TO 620 ; END IF
     END IF

     reallocate = .TRUE.
     IF ( ALLOCATED( WK_t ) ) THEN
       IF ( SIZE( WK_t ) < MAX( ng, n ) ) THEN ; DEALLOCATE( WK_t )
       ELSE ; reallocate = .FALSE. ; END IF
     END IF
     IF ( reallocate ) THEN
       ALLOCATE( WK_t( MAX( ng, n ) ), STAT = alloc_status )
       IF ( alloc_status /= 0 ) THEN ; bad_alloc = 'WK_t'; GO TO 620 ; END IF
     END IF

!  Consider the IG-th group

     nnzgrj = 0
     DO ig = 1, ng
       ig1 = ig + 1
       lstagv = ISTAGV( ig )
       lendgv = ISTAGV( ig1 ) - 1

!  Initialize the group derivative to zero

       WK_n( ISVGRP( lstagv : lendgv ) ) = zero

!  See if the group has any nonlinear elements

       DO ii = ISTADG( ig ), ISTADG( ig1 ) - 1
         iel = IELING( ii )
         ieintv = INTVAR( iel )
         iestev = ISTAEV( iel )
         nvarel = ISTAEV( iel + 1 ) - iestev
         scalee = ESCALE( ii )
         IF ( INTREP( iel ) ) THEN

!  The IEL-th element has an internal representation

           nin = INTVAR( iel + 1 ) - ieintv
           CALL RANGE ( iel, .TRUE., FUVALS( ieintv : ieintv + nin - 1 ),      &
                        WK_t( : nvarel ), nvarel, nin, ITYPEE( iel ),          &
                        nin, nvarel )
!DIR$ IVDEP
           DO i = 1, nvarel
             j = IELVAR( iestev )
             WK_n( j ) = WK_n( j ) + scalee * WK_t( i )
             iestev = iestev + 1
           END DO
         ELSE

!  The IEL-th element has no internal representation

!DIR$ IVDEP
           DO i = 1, nvarel
             j = IELVAR( iestev )
             WK_n( j ) = WK_n( j ) + scalee * FUVALS( ieintv )
             ieintv = ieintv + 1
             iestev = iestev + 1
           END DO
         END IF
       END DO

!  Include the contribution from the linear element

!DIR$ IVDEP
       DO k = ISTADA( ig ), ISTADA( ig1 ) - 1
         IF ( ICNA( k ) <= n ) THEN
           j = ICNA( k )
           WK_n( j ) = WK_n( j ) + A( k )
         END IF
       END DO

!  Store the values of the nonzero entries of the gradient of the IG-th group
!  in GRJAC, along with their row and column indices in the overall Jacobian

       IF ( .NOT. GXEQX( ig ) ) THEN
         gi = GSCALE( ig ) * GVALS( ig , 2 )
       ELSE
         gi = GSCALE( ig )
       END IF
!DIR$ IVDEP
       DO j = lstagv, lendgv
         i = ISVGRP( j )
         nnzgrj = nnzgrj + 1
         IRNGRJ( nnzgrj ) = ig
         ICNGRJ( nnzgrj ) = i
         GRJAC( nnzgrj ) = gi * WK_n( i )
         IF ( iprint >= 100 ) WRITE( iout, 2070 ) ig, i, GRJAC( nnzgrj )
       END DO
     END DO

!  3. Remove all entries which are smaller than TOL times the
!     largest entries in both their rows and columns

!  Find the largest entries, in absolute value, in each row and column

     tol = SQRT( EPSILON( one ) )
     WK_n = zero; WK_t( : ng ) = zero
     DO k = 1, nnzgrj
       i = IRNGRJ( k )
       j = ICNGRJ( k )
       gi = ABS( GRJAC( k ) )
       WK_t( i ) = MAX( WK_t( i ), gi )
       WK_n( j ) = MAX( WK_n( j ), gi )
     END DO

!  Remove small entries

     l = 0
     DO k = 1, nnzgrj
       i = IRNGRJ( k )
       j = ICNGRJ( k )
       gi = ABS( GRJAC( k ) ) / tol
       IF ( gi > WK_t( i ) .OR. gi > WK_n( j ) ) THEN
         l = l + 1
         IRNGRJ( l ) = i
         ICNGRJ( l ) = j
         GRJAC( l ) = GRJAC( k )
       ELSE
         IF ( iprint >= 100 ) WRITE( iout, 2080 ) i, j, GRJAC( k )
       END IF
     END DO
     nnzgrj = l

!  4. Apply the row and column equilibrating scheme of Curtis and Reid to the
!     Jacobian

     CALL SCALN_curtis_reid( ng, n, nnzgrj, GRJAC, IRNGRJ,                     &
                             ICNGRJ, GRPSCA, VSCALE( : n ) )

!  5. Calculate the scale factors

!  Obtain the smallest variable and group scale factor

     grpmax = zero
     varmin = HUGE( one )
     VSCALE( : n ) = EXP( VSCALE( : n ) )
     varmin = MIN( MINVAL( VSCALE( : n ) ), varmin )
!DIR$ IVDEP
     DO ig = 1, ng
       IF ( KNDOFG( ig ) /= 1 ) THEN
         GRPSCA( ig ) = EXP( GRPSCA( ig ) )
         grpmax = MAX( GRPSCA( ig ),grpmax )
       END IF
     END DO

!  Scale the factors relative to their largest members

     VSCALE( : n ) = ten ** ANINT( LOG10( VSCALE( : n ) / varmin ) )
     WHERE ( KNDOFG( : ng ) /= 1 ) GRPSCA( : ng ) =                            &
             ten ** ANINT( LOG10( GRPSCA( : ng ) / grpmax ) )

!  Use the variable scalings within LANCELOT

     IF ( scalev ) THEN
       scale = zero
!DIR$ IVDEP
       VSCALE( : n ) = VSCALE( : n )
       scale = MAX( scale, MAXVAL( VSCALE( : n ) ) )
       stopga = stopga / scale
     END IF

!  Use the group scalings within LANCELOT. Scale the Lagrange multipliers

     IF ( scaleg ) THEN
       scale = one
!DIR$ IVDEP
       DO ig = 1, ng
         IF ( KNDOFG( ig ) /= 1 ) THEN
           GSCALE( ig ) = GSCALE( ig ) * GRPSCA( ig )
           Y( ig ) = Y( ig ) / GRPSCA( ig )
           scale = MIN( scale, GSCALE( ig ) )
         END IF
       END DO
       stopca = stopca * scale
     END IF

!  Successful conclusion to the calculation

     IF ( iprint > 0 ) THEN
       IF ( scaleg ) THEN
         WRITE( iout, 2040 )
         j = 1
         DO i = 1, ng
           IF ( KNDOFG( i ) == 1 ) THEN
             IF ( i - 1 >= j )                                                 &
               WRITE( iout, 2060 ) ( ig, GRPSCA( ig ), ig = j, i - 1 )
              j = i + 1
           END IF
         END DO
         IF ( ng >= j ) WRITE( iout, 2060 ) ( ig, GRPSCA( ig ), ig = j, ng )
       END IF
       IF ( scalev ) WRITE( iout, 2050 ) ( i, VSCALE( i ),i = 1, n )
       WRITE( iout, 2090 ) stopga, stopca
     END IF
     status = 0

!  De-allocate workspace arrays

     DEALLOCATE( IRNGRJ, ICNGRJ, GRJAC, GRPSCA, WK_n, WK_t )
     IF ( fdgrad .AND. nel > 0 ) DEALLOCATE( X_temp, FUVALS_temp )

 610 CONTINUE
     DEALLOCATE( ISTAGV, ISVGRP )

     RETURN

!  Unsuccessful returns

 620 CONTINUE
     status = 1000 + alloc_status
     WRITE( iout, 2600 ) bad_alloc, alloc_status
     RETURN

!  Non-executable statements

 2040  FORMAT( /,' Multiply groups by the factors:' )
 2050  FORMAT( /,' Divide variables by the factors:', / 4( I6, ES12.4 ) )
 2060  FORMAT( 4( I6, ES12.4 ) )
 2070  FORMAT( ' Row ', I6, ' column ', I6, ' value ', ES12.4 )
 2080  FORMAT( ' Row ', I6, ' column ', I6, ' value ', ES12.4,' removed' )
 2090  FORMAT( /,' Scaled projected gradient tolerance = ', ES12.4, /,         &
                 ' Scaled constraint tolerance         = ', ES12.4 )
 2100  FORMAT( ' ** Return from SCALN. Insufficient workspace ' )
 2600  FORMAT( ' ** Message from -SCALN_get_scalings-', /,                     &
               ' Allocation error, for ', A6, ', status = ', I6 )

!  End of subroutine SCALN_get_scalings_main

     END SUBROUTINE SCALN_get_scalings_main

!-*-*-*-  L A N C E L O T  -B-  SCALN_curtis_reid  S U B R O U T I N E  -*-*

     SUBROUTINE SCALN_curtis_reid(  m , n , ne, A , IRNA  , ICNA  , R , C )

!   Compute row and column scalings using the algorithm of Curtis and Reid
!   (J.I.M.A. 10 (1972) 118-124) by approximately minimizing the function

!        sum (nonzero A) ( log(|a_ij|) - r_i - c_j)^2

!   The required scalings are then e^-r and e^-c respectively

!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------

     INTEGER, INTENT( IN ) :: m , n , ne
     INTEGER, INTENT( IN  ), DIMENSION( ne ) :: IRNA, ICNA
     REAL ( KIND = wp ), INTENT( IN  ), DIMENSION( ne ) :: A
     REAL ( KIND = wp ), INTENT( OUT ), DIMENSION( m  ) :: R
     REAL ( KIND = wp ), INTENT( OUT ), DIMENSION( n  ) :: C

!-----------------------------------------------
!   L o c a l   P a r a m e t e r s
!-----------------------------------------------

     INTEGER, PARAMETER :: itmax = 100

!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------

     INTEGER :: i, j, k, iter
     REAL ( KIND = wp ) :: e, e_old, eprod, q, q_old
     REAL ( KIND = wp ) :: ss, ss_old, stop_tol, loga, val, absa

!-----------------------------------------------
!   A u t o m a t i c   A r r a y s
!-----------------------------------------------

     REAL ( KIND = wp ), DIMENSION( m ) :: ROW_COUNT, M_INV_SIG
     REAL ( KIND = wp ), DIMENSION( n ) :: COL_COUNT, COL_RHS, SHIFT_C

     IF ( m < 1 .OR. n < 1 ) RETURN

!  Set the stopping tolerance

     stop_tol = ne * point1

!  Initialise for accumulation of sums and products

     R = zero ; C = zero
     ROW_COUNT = zero ; COL_COUNT = zero ; COL_RHS = zero

!  Count non-zeros in the rows, and compute r.h.s. vectors; use R to store
!  the row r.h.s. (sigma in Curtis+Reid)

     DO k = 1, ne
       absa = ABS( A( k ) )
       IF ( absa /= zero ) THEN
         i = IRNA( k ) ; j = ICNA( k )
         loga = LOG( absa )
         ROW_COUNT( i ) = ROW_COUNT( i ) + one
         COL_COUNT( j ) = COL_COUNT( j ) + one
         R ( i )      = R ( i ) + loga
         COL_RHS( j ) = COL_RHS( j ) + loga
       END IF
     END DO

!  Account for structural singularity

     WHERE ( ROW_COUNT == zero ) ROW_COUNT = one
     WHERE ( COL_COUNT == zero ) COL_COUNT = one

!  Use Reid's special purpose method for matrices with property "A". Comments
!  refer to equation numbers in Curtis and Reid's paper

!  Form M^-1 sigma and N^-1 tau (in C+R's notation)

     M_INV_SIG = R / ROW_COUNT
     COL_RHS = COL_RHS / COL_COUNT

!  Compute initial residual vector

     R = M_INV_SIG
     DO k = 1, ne
       IF ( A( k ) /= zero ) THEN
         i = IRNA( k ) ; j = ICNA( k )
         R( i ) = R( i ) - COL_RHS( j ) / ROW_COUNT( i )  ! (4.3)
       END IF
     END DO

!  Set initial values

     e = zero ; q = one ;
!    ss = DOT_PRODUCT( ROW_COUNT, R ** 2 )
     ss = zero ; DO k = 1, m ; ss = ss + ROW_COUNT( k ) * R( k ) ** 2 ; END DO
     IF ( ss > stop_tol ) THEN

!  Iteration loop

       SHIFT_C = zero
       DO iter = 1, itmax

!  Update column residual vector

         DO k = 1, ne
           IF ( A( k ) /= zero ) THEN
             i = IRNA( k ) ; j = ICNA( k )
             C( j ) = C( j ) + R( i )
           END IF
         END DO

!  Rescale column residual

         ss_old = ss ; ss = zero
         DO j = 1, n
           val = - C( j ) / q
           C( j ) = val / COL_COUNT( j )     ! (4.4a)
           ss = ss + val * C( j )            ! (4.5a)
         END DO

!  Test for termination

          IF ( ss <= stop_tol ) THEN
            R = zero ; GO TO 100
          END IF
          e_old = e ; e = q * ss / ss_old   ! (4.6)
          q = one - e                       ! (4.7)
          R = R * e * ROW_COUNT

!  Rescale row residual vector

          eprod = e * e_old

!  Update row residual

          DO k = 1, ne
            IF ( A( k ) /= zero ) THEN
              i = IRNA( k ) ; j = ICNA( k )
              R( i ) = R( i ) + C( j )
            END IF
          END DO

!  Again, rescale row residual

          ss_old = ss ; ss = zero
          DO i = 1, m
            val = - R( i ) / q
            R( i ) = val / ROW_COUNT( i )    ! (4.4b)
            ss = ss + val * R( i )           ! (4.5b)
          END DO

!  Test for termination

          IF ( ss <= stop_tol ) THEN
            COL_RHS = COL_RHS + ( C + eprod * SHIFT_C ) / q ; EXIT
          END IF
          e_old = e ; e = q * ss / ss_old   ! (4.6)
          q_old = q ; q = one - e           ! (4.7)

!  Rescale column residual vector

          SHIFT_C = ( C + eprod * SHIFT_C ) / ( q * q_old )
          COL_RHS = COL_RHS + SHIFT_C

!  Update column scaling factors

          C = C * e * COL_COUNT
        END DO
      END IF

!  Rescale row scaling factors

     R = R * ROW_COUNT

!  Rescale row scaling factors

 100 CONTINUE
     DO k = 1, ne
       IF ( A( k ) /= zero ) THEN
         i = IRNA( k )
         j = ICNA( k )
         R( i ) = R( i ) + COL_RHS( j )
       END IF
     END DO

!  Final rescaling of row and column scaling factors

     R = R / ROW_COUNT - M_INV_SIG ; C = - COL_RHS
     RETURN

!  End of subroutine SCALN_curtis_reid

     END SUBROUTINE SCALN_curtis_reid

!  End of module LANCELOT_SCALN

   END MODULE LANCELOT_SCALN_double

max_forks_repo_name: dynamics-of-stellar-systems/dynamite


text:
      program nga02

C     Extraction de conditions de non-debordement: you want N<=33 before
C     the print statement. It works with the derivative fix-point

      real a(33)

      read *, n

      if(n.lt.5) stop

      do i = 1, n
         if(i.gt.33) stop
         a(i) = 0.
      enddo

      print *, a, n

      end

ext: f


program modif

USE netcdf

IMPLICIT NONE

INTEGER :: fidWOA, status, dimID_nbounds, dimID_z, dimID_y, dimID_x, mnbounds, mz, my, mx, thermal_forcing_ID, z_bnds_ID, z_ID, y_ID, x_ID, fidDFT, my2, mx2, thickness_ID, surface_ID, fidM, ii, jj, kk, kinf, ksup

CHARACTER(LEN=150) :: file_WOA, file_out, file_DFT

REAL*8,ALLOCATABLE,DIMENSION(:) :: x, y, z

REAL*8,ALLOCATABLE,DIMENSION(:,:) :: thickness, surface, ice_draft, z_bnds, TF_draft

REAL*8,ALLOCATABLE,DIMENSION(:,:,:) :: thermal_forcing

LOGICAL :: ll_z_upward

file_DFT  = 'bedmap2_8km.nc'
!file_WOA  = 'woa_meop_thermal_forcing_1955-2012_8km_x_60m.nc'
file_WOA  = 'obs_thermal_forcing_1995-2017_8km_x_60m.nc'
file_out  = 'obs_thermal_forcing_on_bedmap2_ice_draft.nc'

!---------------------------------------
! Read 3D thermal forcing :

write(*,*) 'Reading ', TRIM(file_WOA)

status = NF90_OPEN(TRIM(file_WOA),0,fidWOA); call erreur(status,.TRUE.,"read")

status = NF90_INQ_DIMID(fidWOA,"nbounds",dimID_nbounds); call erreur(status,.TRUE.,"inq_dimID_nbounds")
status = NF90_INQ_DIMID(fidWOA,"z",dimID_z); call erreur(status,.TRUE.,"inq_dimID_z")
status = NF90_INQ_DIMID(fidWOA,"y",dimID_y); call erreur(status,.TRUE.,"inq_dimID_y")
status = NF90_INQ_DIMID(fidWOA,"x",dimID_x); call erreur(status,.TRUE.,"inq_dimID_x")

status = NF90_INQUIRE_DIMENSION(fidWOA,dimID_nbounds,len=mnbounds); call erreur(status,.TRUE.,"inq_dim_nbounds")
status = NF90_INQUIRE_DIMENSION(fidWOA,dimID_z,len=mz); call erreur(status,.TRUE.,"inq_dim_z")
status = NF90_INQUIRE_DIMENSION(fidWOA,dimID_y,len=my); call erreur(status,.TRUE.,"inq_dim_y")
status = NF90_INQUIRE_DIMENSION(fidWOA,dimID_x,len=mx); call erreur(status,.TRUE.,"inq_dim_x")

ALLOCATE(  thermal_forcing(my,mx,mz)  )
ALLOCATE(  TF_draft(my,mx)  )
ALLOCATE(  z_bnds(mnbounds,mz)  )
ALLOCATE(  z(mz)  )
ALLOCATE(  y(my)  )
ALLOCATE(  x(mx)  )

status = NF90_INQ_VARID(fidWOA,"thermal_forcing",thermal_forcing_ID); call erreur(status,.TRUE.,"inq_thermal_forcing_ID")
status = NF90_INQ_VARID(fidWOA,"z_bnds",z_bnds_ID); call erreur(status,.TRUE.,"inq_z_bnds_ID")
status = NF90_INQ_VARID(fidWOA,"z",z_ID); call erreur(status,.TRUE.,"inq_z_ID")
status = NF90_INQ_VARID(fidWOA,"y",y_ID); call erreur(status,.TRUE.,"inq_y_ID")
status = NF90_INQ_VARID(fidWOA,"x",x_ID); call erreur(status,.TRUE.,"inq_x_ID")

status = NF90_GET_VAR(fidWOA,thermal_forcing_ID,thermal_forcing); call erreur(status,.TRUE.,"getvar_thermal_forcing")
status = NF90_GET_VAR(fidWOA,z_bnds_ID,z_bnds); call erreur(status,.TRUE.,"getvar_z_bnds")
status = NF90_GET_VAR(fidWOA,z_ID,z); call erreur(status,.TRUE.,"getvar_z")
status = NF90_GET_VAR(fidWOA,y_ID,y); call erreur(status,.TRUE.,"getvar_y")
status = NF90_GET_VAR(fidWOA,x_ID,x); call erreur(status,.TRUE.,"getvar_x")

status = NF90_CLOSE(fidWOA); call erreur(status,.TRUE.,"close_file")

!---------------------------------------
! Read BEDMAP2's ice draft :

write(*,*) 'Reading ', TRIM(file_DFT)

status = NF90_OPEN(TRIM(file_DFT),0,fidDFT); call erreur(status,.TRUE.,"read")

status = NF90_INQ_DIMID(fidDFT,"y",dimID_y); call erreur(status,.TRUE.,"inq_dimID_y")
status = NF90_INQ_DIMID(fidDFT,"x",dimID_x); call erreur(status,.TRUE.,"inq_dimID_x")

status = NF90_INQUIRE_DIMENSION(fidDFT,dimID_y,len=my2); call erreur(status,.TRUE.,"inq_dim_y")
status = NF90_INQUIRE_DIMENSION(fidDFT,dimID_x,len=mx2); call erreur(status,.TRUE.,"inq_dim_x")

if ( mx .ne. mx2 .or. my .ne. my2 ) then
  write(*,*) '~!@#$%^* ERROR : dimension mismatch >>>>> stop'
  stop
endif

ALLOCATE(  thickness(mx,my)  )
ALLOCATE(  surface(mx,my)  )
ALLOCATE(  ice_draft(mx,my)  )

status = NF90_INQ_VARID(fidDFT,"thickness",thickness_ID); call erreur(status,.TRUE.,"inq_thickness_ID")
status = NF90_INQ_VARID(fidDFT,"surface",surface_ID); call erreur(status,.TRUE.,"inq_surface_ID")

status = NF90_GET_VAR(fidDFT,thickness_ID,thickness); call erreur(status,.TRUE.,"getvar_thickness")
status = NF90_GET_VAR(fidDFT,surface_ID,surface); call erreur(status,.TRUE.,"getvar_surface")

status = NF90_CLOSE(fidDFT); call erreur(status,.TRUE.,"close_file")

ice_draft(:,:) = surface(:,:) - thickness(:,:)
DEALLOCATE( surface, thickness )

!----------------------------------------------------------------------
! Linear interpolation of the thermal forcing on the ice draft depth :

if ( z(2) .gt. z(1) ) then
  ll_z_upward = .true.
else
  ll_z_upward = .false.
endif

do ii=1,mx
do jj=1,my

  if ( ice_draft(ii,jj) .lt. 0.d0 ) then

    if ( ll_z_upward ) then
      kinf=1
      do kk=2,mz-1
        if ( z(kk) .le. ice_draft(ii,jj) ) kinf = kk
      enddo
      ksup = kinf + 1
    else
      ksup=mz
      do kk=mz-1,2,-1
        if ( z(kk) .le. ice_draft(ii,jj) ) ksup = kk
      enddo
      kinf = ksup - 1
    endif
    write(*,*) 'zinf, ice_draft, zsup = ', z(kinf), ice_draft(ii,jj), z(ksup)

    TF_draft(ii,jj) = (   (z(ksup)-ice_draft(ii,jj)) * thermal_forcing(ii,jj,kinf)   &
    &                   + (ice_draft(ii,jj)-z(kinf)) * thermal_forcing(ii,jj,ksup) ) / (z(ksup)-z(kinf))

    write(*,555) ii, jj, kinf, ksup, z(kinf), z(ksup), ice_draft(ii,jj), thermal_forcing(ii,jj,kinf), thermal_forcing(ii,jj,ksup), TF_draft(ii,jj)

  else

    TF_draft(ii,jj) = 0.d0 ! will ensure zero melt in case of unexpected overlap
    !write(*,*) ii, jj, '........'

  endif

enddo
enddo

555 FORMAT(i4,i4,i3,i3,' ',f7.1,' ',f7.1,' ',f7.1,' ',f5.2,' ',f5.2,' ',f5.2)

!---------------------------------------
! Writing new netcdf file :

write(*,*) 'Creating ', TRIM(file_out)

status = NF90_CREATE(TRIM(file_out),NF90_NOCLOBBER,fidM); call erreur(status,.TRUE.,'create')

status = NF90_DEF_DIM(fidM,"y",my,dimID_y); call erreur(status,.TRUE.,"def_dimID_y")
status = NF90_DEF_DIM(fidM,"x",mx,dimID_x); call erreur(status,.TRUE.,"def_dimID_x")

status = NF90_DEF_VAR(fidM,"TF_draft",NF90_DOUBLE,(/dimID_y,dimID_x/),thermal_forcing_ID); call erreur(status,.TRUE.,"def_var_thermal_forcing_ID")
status = NF90_DEF_VAR(fidM,"y",NF90_DOUBLE,(/dimID_y/),y_ID); call erreur(status,.TRUE.,"def_var_y_ID")
status = NF90_DEF_VAR(fidM,"x",NF90_DOUBLE,(/dimID_x/),x_ID); call erreur(status,.TRUE.,"def_var_x_ID")

status = NF90_PUT_ATT(fidM,thermal_forcing_ID,"long_name","thermal forcing on the ice draft"); call erreur(status,.TRUE.,"put_att_thermal_forcing_ID")
status = NF90_PUT_ATT(fidM,thermal_forcing_ID,"units","degrees_celsius"); call erreur(status,.TRUE.,"put_att_thermal_forcing_ID")
!status = NF90_PUT_ATT(fidM,thermal_forcing_ID,"_FillValue",0.d0); call erreur(status,.TRUE.,"put_att_thermal_forcing_ID")
status = NF90_PUT_ATT(fidM,y_ID,"units","meters"); call erreur(status,.TRUE.,"put_att_y_ID")
!status = NF90_PUT_ATT(fidM,y_ID,"_FillValue",0.d0); call erreur(status,.TRUE.,"put_att_y_ID")
status = NF90_PUT_ATT(fidM,x_ID,"units","meters"); call erreur(status,.TRUE.,"put_att_x_ID")
!status = NF90_PUT_ATT(fidM,x_ID,"_FillValue",0.d0); call erreur(status,.TRUE.,"put_att_x_ID")

status = NF90_PUT_ATT(fidM,NF90_GLOBAL,"history","Created using interpolate_TF_on_bedmap2_ice_draft.f90"); call erreur(status,.TRUE.,"put_att_GLOBAL_ID")

status = NF90_ENDDEF(fidM); call erreur(status,.TRUE.,"fin_definition")

status = NF90_PUT_VAR(fidM,thermal_forcing_ID,TF_draft); call erreur(status,.TRUE.,"var_thermal_forcing_ID")
status = NF90_PUT_VAR(fidM,y_ID,y); call erreur(status,.TRUE.,"var_y_ID")
status = NF90_PUT_VAR(fidM,x_ID,x); call erreur(status,.TRUE.,"var_x_ID")

status = NF90_CLOSE(fidM); call erreur(status,.TRUE.,"final")

end program modif

SUBROUTINE erreur(iret, lstop, chaine)
  ! pour les messages d'erreur
  USE netcdf
  INTEGER, INTENT(in)                     :: iret
  LOGICAL, INTENT(in)                     :: lstop
  CHARACTER(LEN=*), INTENT(in)            :: chaine
  !
  CHARACTER(LEN=80)                       :: message
  !
  IF ( iret .NE. 0 ) THEN
    WRITE(*,*) 'ROUTINE: ', TRIM(chaine)
    WRITE(*,*) 'ERREUR: ', iret
    message=NF90_STRERROR(iret)
    WRITE(*,*) 'CA VEUT DIRE:',TRIM(message)
    IF ( lstop ) STOP
  ENDIF
  !
END SUBROUTINE erreur


text:
great_ne(q1,p1).
great_ne(k1,c1).
great_ne(jj1,w1).
great_ne(q1,jj1).
great_ne(cc1,k1).
great_ne(p1,aa1).
great_ne(jj1,hh1).
great_ne(z1,p1).
great_ne(t1,dd1).
great_ne(t1,u1).
great_ne(m1,y1).
great_ne(w1,b1).
great_ne(bb1,ii1).
great_ne(u1,hh1).
great_ne(dd1,m1).
great_ne(w1,c1).
great_ne(b1,h1).
great_ne(hh1,k1).
great_ne(q1,z1).
great_ne(z1,ii1).
great_ne(z1,aa1).
great_ne(q1,m1).
great_ne(t1,ii1).
great_ne(z1,w1).
great_ne(p1,jj1).
great_ne(q1,ff1).
great_ne(z1,x1).
great_ne(ii1,i1).
great_ne(z1,t1).
great_ne(y1,w1).
great_ne(hh1,ff1).
great_ne(ee1,b1).
great_ne(aa1,k1).
great_ne(q1,w1).

max_stars_repo_name: manoelfranca/cilppp

max_issues_repo_name: manoelfranca/cilppp

max_forks_repo_name: manoelfranca/cilppp


!
!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
subroutine dstuart (nu, xo, tridlc, ss,ds,op, u,t)

!     Calls Stuart's angle subroutine to get triangle dislocation results.

!     Comninou and Dundurs use a 1,2,3 = NED coordinate system.
!       Their simple angular dislocation is in the 1,3 plane
!       with its tip under the origin and its arm pointing toward +x1.
  implicit none
  integer :: i,k,iflag
  real(8) :: nu,ss,ds,op
  real(8) :: xo(3), tridlc(9), burg(3), u(3), t(9)
  real(8) :: utmp(3), ttmp(9)
  real(8) :: tri(3,4)
  integer :: nwarn
  common/warn0/ nwarn

!      nu         = poissons ratio
!      xo(3)      = observation point in global NED=123 coordinate system
!      tridlc(9)  = x2,x1,x3 coordinates for the three triangle corners.
!      ss,ds,op   = burg(3)    = burgers vector
!                     (ss = + for LL; ds = + for normal; op = + for expansion)
!      u(3)       = displacements = u1,u2,u3
!      t(9)       = tilts = t11,t12,t13,t21,t22,t23,t31,t32,t33
!                     (where tij = d(ui)/dxj)

  real(8), parameter :: pi=3.14159265358979323846d0
  real(8), parameter :: deg2rad=pi/180.d0
  real(8), parameter :: rad2deg=180.d0/pi

!.....Initialize displacement and tilt arrays.
  u = 0.d0
  t = 0.d0

!.....Fill array with triangle corner coordinates.
!.......(Corner 4 is corner 1 repeated.)
  do i=1,4
    do  k=1,3
      tri(k,i) = tridlc(k+3*mod(i-1,3))
    end do
  end do

!.....Express burger's vector in global 123=NED coordinates.
      call burger2global (ss,ds,op, tri, burg)
!        print '(a,6f8.3)', 'ss,op,ds,burg =(stu) ', ss,op,ds,burg

!.....For each of 3 sides, calculate fields.
  do k = 1,3
    call twoangles (nu,xo,tri(1,k),tri(1,k+1),burg,utmp,ttmp)
    u = u + utmp
    t = t + ttmp
  end do

!.....Apply additional displacement if pt is under the triangle.
  call undertriangle (xo, tri, iflag)

  if (iflag .eq. 1) then
!        Point lies under interior of triangle.
    u = u - burg
  else if (iflag .ne. 0) then
!        Point lies on triangle or on edge of prism below triangle.
    nwarn = nwarn + 1
    if (nwarn.lt.10) then
      print *, ' *** Warning, possible singular location'
      print *, '      for point = ', xo
    else if (nwarn.eq.11) then
      print *, ' *** Warning, possible singular location etc...'
    endif
    u = u - 0.5*burg
  endif

  return
end subroutine dstuart

subroutine burger2global (ss,ds,op, tri, burg)

!.....Express burger's vector in global 123=NED coordinates.
!     Note special definition of ss and ds needed for horizontal dislocations
!         and of ds for vertical dislocations.
  implicit none
  real(8) :: ss,ds,op
  real(8) :: tri(3,3), burg(3)
  real(8) :: side12(3), side13(3)
  real(8) :: perp(3), hor(3), dwndip(3), b(3)
  real(8), parameter :: pi=3.14159265358979323846d0
  real(8), parameter :: deg2rad=pi/180.d0
  real(8), parameter :: rad2deg=180.d0/pi
  real(8) :: strikerad,strike,dip

!     tri(3,i) = coordinates of ith corner in 123 coordinate system.
!                   (i=4 returns to i=1)

!.....For a horizontal triangle define strike- and dip-directions by fiat.
  if (tri(3,1).eq.tri(3,2) .and. tri(3,1).eq.tri(3,3)) then
!        Let 1 = ss-axis, 2 = ds-axis, 3 = op-axis
    burg(1) = ss
    burg(2) = ds
    burg(3) = -op
    return
  endif

!.....Otherwise, construct a coordinate system attached to the triangle.

!.....Calculate a normal vector.
  side12(:) = tri(:,2) - tri(:,1)
  side13(:) = tri(:,3) - tri(:,1)
  call cross (side12,side13,perp)
  call unit (perp)

!.....Vertical triangle.
!       By definition, for a vertical triangle, side from which corners
!          appear clockwise is considered top side.

!.....For dipping triangle...
  if (perp(3).lt.0.d0) then
!       Dipping triangle has CCW order to corners viewed from above.
!        Switch corners 2 and 3 and reverse perp.
    call switch3(tri(1,2),tri(1,3))
    perp(:) = -1.d0*perp(:)
!       Perp should now be a unit normal vector pointing down,
!        (away from top), unless triangle is vertical.
  endif

!.....Otherwise define a horizontal vector along strike.
!      (Strike is to left hand as one stands looking down dip
!         -- perp points opposite down-dip.)
  strikerad = atan2(perp(2),perp(1)) + pi*0.5d0
  hor(1) = cos(strikerad)
  hor(2) = sin(strikerad)
  hor(3) = 0.d0

!.....Next a vector in the plane of the triangle pointing in the
!       down-dip direction.
  call cross (perp,hor,dwndip)

!.....Express Burgers vector in this coordinate system.
  burg(:) = ss*hor(:) + ds*dwndip(:) - op*perp(:)

!.....Do a check.
  strike = strikerad * rad2deg
  dip = acos(-perp(3)) * rad2deg
  call global2triburger (burg, strike, dip, b)

  return
end subroutine burger2global

subroutine global2triburger (burg, strike, dip, b)

!.....Express burger's vector as b="ss,op,ds" (i.e., in local coordinates
!      attached to dlc.)  Order is here set to match Comninou & Dunders
!      conventions.
  implicit none
  real(8) :: burg(3), b(3)
  real(8) :: perp(3), hor(3), dwndip(3)
  real(8), parameter :: pi=3.14159265358979323846d0
  real(8), parameter :: deg2rad=pi/180.d0
  real(8), parameter :: rad2deg=180.d0/pi
  real(8) :: strike,dip,strikerad,diprad,dipdirrad

!.....Construct a coordinates system attached to the triangle.

!.....For a horizontal triangle define strike- and dip-directions by fiat.
  if (dip .eq. 0.d0) then
!        Let 1 = ss-axis, 2 = op-axis, 3 = ds-axis
    strike = 0.d0
  endif

!.....Define a horizontal vector along strike.
  strikerad = strike * deg2rad
  hor(1) = cos(strikerad)
  hor(2) = sin(strikerad)
  hor(3) = 0.d0

!.....Next a vector in the plane of the triangle pointing in the
!       down-dip dipdirection.
  diprad = dip * deg2rad
  dipdirrad = strikerad + pi*0.5
  dwndip(1) = cos (diprad) * cos(dipdirrad)
  dwndip(2) = cos (diprad) * sin(dipdirrad)
  dwndip(3) = sin (diprad)

!.....Finally a perpendicular vector
  call cross (hor,dwndip,perp)

!.....Get local burgers vector.
  b(1) =   dot_product(burg,hor)
  b(2) = - dot_product(burg,perp)
  b(3) =   dot_product(burg,dwndip)

  return
end subroutine global2triburger

subroutine twoangles (nu,xo,xcorna,xcornb, burg, u,t)

!     Calculates the effect of a vertical side formed by two angles
!       located at points xcorna and xcornb.
  implicit none
  real(8) :: nu, xo(3), xcorna(3), xcornb(3), burg(3), u(3), t(9)
  real(8) :: xop(3), xca(3), xcb(3), b(3)
  real(8) :: v(3,3), dv(3,3,3)
  real(8) :: va(3,3), dva(3,3,3)
  real(8) :: vb(3,3), dvb(3,3,3)

!     nu     = Poissons ratio
!     xo(3)  = observation point in x1,x2,x3 = NED coordinates.
!     xca(3) = x1,x2,x3 coordinates of corner A.
!     xcb(3) = x1,x2,x3 coordinates of corner B.
  real(8), parameter :: pi=3.14159265358979323846d0
  real(8), parameter :: deg2rad=pi/180.d0
  real(8), parameter :: rad2deg=180.d0/pi
  integer :: iswitch,j,k,iret
  real(8) :: theta,strike,dip,rot,deptha,depthb,xlen,betarad,x1,x2,x3

!.....Copy corners.
  xca(:) = xcorna(:)
  xcb(:) = xcornb(:)

!.....See which is shallower, A or B, and switch if necessary.
  if (xca(3).gt.xcb(3)) then
!       B is shallower than A
    call switch3 (xca, xcb)
    iswitch = 1
  else
    iswitch = 0
  endif

!.....Calculate angle from A to B.
  theta = atan2(xcb(2)-xca(2), xcb(1)-xca(1)) * rad2deg

!.....Set parameters.
  strike = theta
  deptha = xca(3)
  depthb = xcb(3)
  xlen = sqrt((xcb(1)-xca(1))**2 + (xcb(2)-xca(2))**2)
  betarad =  atan2 (xlen, xcb(3)-xca(3))

!.....Locate the observation point relative to the shallower of A or B.
  xop(1) =  xo(1) - xca(1)
  xop(2) =  xo(2) - xca(2)
  xop(3) =  xo(3)
!     If strike .ne. 0, need to rotate obs points and output fields.
  rot = -strike
  call vector_rot3(xop,-rot)
  x1 = xop(1)
  x2 = xop(2)
  x3 = xop(3)

!.....Compute the displacement and tilts.
  call comdun2 ( nu, x1,x2,x3, deptha, betarad, va, dva, iret )
  call comdun2 ( nu, x1-xlen,x2,x3, depthb, betarad, vb, dvb, iret )

!.....Combine the results from both angles.
  v(:,:) = va(:,:) - vb(:,:)
  dv(:,:,:) = dva(:,:,:) - dvb(:,:,:)

!.....Calculate displacements(u) and tilts(t)
  dip = 90.d0
  call global2triburger (burg, strike, dip, b)

  if (iswitch.eq.1)  b(:) = -1d0*b(:)

  u = matmul(b,v)
  do j=1,3
    do k=1,3
      t(k + 3*(j-1)) = dot_product(b(:),dv(:,j,k))
    end do
  end do

!.....Rotate fields back to global NED=123 coordinate system.
  call vector_rot3 (u,rot)
  call tensor_rot3 (t,rot)

  return
end subroutine twoangles

subroutine comdun2 ( nu, x1, x2, x3, a, beta, v, dv, iret )
!     Performs some checks before calling comdun
  implicit none
  real(8) :: nu,x1,x2,x3,a,beta
  integer :: iret,nwarn
  real(8) :: v(3,3), dv(3,3,3)
  real(8) :: vp(3,3), dvp(3,3,3)
  real(8) :: vm(3,3), dvm(3,3,3)
  real(8) :: betatol = 1.d-4
  real(8) :: tol = 1.d-4
  real(8) :: x2p,x2m,tan2angle

  common/nwarn00/ nwarn

!.....Check for angle beta - if close to zero, return zero result.
  if (beta .lt. betatol) then
    v(:,:) = 0.d0
    dv(:,:,:) = 0.d0
    return
  endif

!.....Check for distance from x1-x3 plane.  If too close, average.
  if (abs(x2).lt. tol) then
    x2p = abs(x2) + tol
    x2m = -x2p
    call comdun ( nu, x1, x2p, x3, a, beta, vp, dvp, iret )
    call comdun ( nu, x1, x2m, x3, a, beta, vm, dvm, iret )
    v(:,:) = 0.5d0*(vp(:,:)+vm(:,:))
    dv(:,:,:) = 0.5d0*(dvp(:,:,:)+dvm(:,:,:))
  else
!        Normal call to comdun
    call comdun ( nu, x1, x2, x3, a, beta, v, dv, iret )
  endif

!c.....Check to see if point is inside angle.
  tan2angle = atan2(x1,x3-a)
  if (abs(x2).lt.tol .and. (tan2angle.le.beta.and.tan2angle.ge.0.d0)) then
!        Close to the displacement singularity.
    nwarn = nwarn + 1
    if (nwarn.lt.10) then
      print *, '  *** Displacement singularity...'
    else if (nwarn.eq.11) then
      print *, '  *** Displacement singularity...etc,etc'
    endif
    iret = 1
  endif

  return
end subroutine comdun2

subroutine undertriangle (xo, tri, iflag)

!     Returns a flag = 1 if the point xo lies under the triangle interior.
!     Returns a flag = 2 if the point xo lies under a triangle edge or corner.
!     Returns a flag = 3 if the point xo lies on a triangle edge or corner.
!     Returns a flag = 0 otherwise.
  implicit none
  integer :: iflag
  real(8) :: xo(3), tri(3,4)
  real(8) :: xc(3), yc(3)
  real(8) :: side12(3), side13(3), perp(3)
  real(8) :: tol = 1.d-6
  real(8) :: xpt,ypt,zpt,zontri,ht
  integer :: inside

  iflag=0
!     Make arrays of x,y coords of triangle corners.
  xc(1:3) = tri(1,1:3)
  yc(1:3) = tri(2,1:3)

  xpt = xo(1)
  ypt = xo(2)
  zpt = xo(3)
  iflag = inside (xpt,ypt, xc,yc,3)

!     Point is outside vertical prism....
  if (iflag.eq.0)  return

!     Point is inside vertical prism.... is it below triangle?

!     Form perpendicular to the triangle.
  side12(:) = tri(:,2) - tri(:,1)
  side13(:) = tri(:,3) - tri(:,1)
  call cross (side12,side13,perp)
  call unit (perp)

  if (perp(3).eq.0.0) then
!        Triangle is vertical, return edge flag.
    iflag = 2
    return
  endif

!     Calculate height (+) below or (-) above triangle surface.
  zontri = (dot_product(tri(:,1),perp(:)) - xpt*perp(1) - ypt*perp(2))/perp(3)
  ht = zpt - zontri
  if (abs(ht).lt.tol) then
    iflag = 3
  else if (ht.lt.-tol) then
    iflag = 0
  else
    iflag = abs(iflag)
  endif

  return
end subroutine undertriangle

integer function inside (x0,y0, px,py,n)

!     From J.Berger, et al. BSSA v. 74, p. 1849-1862, October 1984.

!     x0,y0 = point to test
!     px,py,n = corners of polygon of n sides.
!     Return value = 0 if point is outside
!                  = +/-1 if point is inside
!                  = 2 if point is on an edge or vertex.
  integer :: i
  integer :: n
  real(8) :: x0,y0,px(n), py(n)
  inside = 0
  do i=1,n-1
    isicr=ksicr(px(i)-x0,py(i)-y0,px(i+1)-x0,py(i+1)-y0)
    if (isicr.eq.4) then
      inside = 2
      return
    endif
    inside = inside + isicr
  end do

  isicr=ksicr(px(n)-x0,py(n)-y0,px(1)-x0,py(1)-y0)

  if(isicr.eq.4) then
    inside = 2
    return
  endif

  inside = (inside + isicr)/2.
  return
end function inside

function ksicr (x1,y1,x2,y2)
  implicit none
  integer :: ksicr
  real(8) :: x1,x2,y1,y2
  if (y1*y2 .gt. 0.0) goto 600

  if (x1*y2 .ne. x2*y1 .or. x1*x2 .gt. 0.0) goto 100

  ksicr = 4
  return

100 if (y1*y2 .lt. 0.0) goto 300
  if (y2.eq.0.0) goto 200
  if (x1.gt.0.0) goto 600
  if (y2.gt.0.0) goto 700
  goto 800

200 if (y1.eq.0.0 .or. x2.gt.0.0) goto 600
  if (y1.gt.0.0) goto 800
  goto 700

300 if (y1.gt.0.0) goto 400
  if (x1*y2 .ge. y1*x2) goto 600
  ksicr = 2
  return

400 if (y1*x2.ge.x1*y2) goto 600
  ksicr = -2
  return

600 ksicr = 0
  return

700 ksicr = 1
  return

800 ksicr = -1
  return

end function ksicr

subroutine vector_rot3 (u,phi)

!     Modified from Laurie Erickson's version of the DIS3D program.

!     Rotates vector u around the x3 axis in a CW sense
!        by angle phi in the x1-x2 plane.
  implicit none
  real(8) :: phi
  real(8) :: u(3), utmp(3)
  real(8), parameter :: pi=3.14159265358979323846d0
  real(8), parameter :: deg2rad=pi/180.d0
  real(8), parameter :: rad2deg=180.d0/pi

  real(8) :: phirad,cp,sp

  if (phi.eq.0.0) return

  phirad = phi*deg2rad
  cp = cos(phirad)
  sp = sin(phirad)

  utmp(1:2) = u(1:2)

  u(1) =  cp*utmp(1) + sp*utmp(2)
  u(2) = -sp*utmp(1) + cp*utmp(2)

  return
end subroutine vector_rot3

subroutine tensor_rot3 (t,phi)

!     Modified from Laurie Erickson's version of the DIS3D program.

!     Rotates a general tensor t by the angle phi around the
!       x3 axis (within the x1-x2 plane).
  real(8) :: phi
  real(8) :: t(9), ttmp(9)
  real(8), parameter :: pi=3.14159265358979323846d0
  real(8), parameter :: deg2rad=pi/180.d0
  real(8), parameter :: rad2deg=180.d0/pi
  real(8) :: phirad,cp,sp

  if (phi.eq.0.d0) return

  phirad = phi*deg2rad
  cp = cos(phirad)
  sp = sin(phirad)

  ttmp(1:8) = t(1:8)

  t(1) =  (cp**2)*ttmp(1)+cp*sp*ttmp(2)+cp*sp*ttmp(4)+(sp**2)*ttmp(5)
  t(2) = -cp*sp*ttmp(1)+(cp**2)*ttmp(2)-(sp**2)*ttmp(4)+sp*cp*ttmp(5)
  t(3) =  cp*ttmp(3) + sp*ttmp(6)
  t(4) = -sp*cp*ttmp(1)-(sp**2)*ttmp(2)+(cp**2)*ttmp(4)+cp*sp*ttmp(5)
  t(5) =  (sp**2)*ttmp(1)-sp*cp*ttmp(2)-sp*cp*ttmp(4)+(cp**2)*ttmp(5)
  t(6) = -sp*ttmp(3) + cp*ttmp(6)
  t(7) =  cp*ttmp(7) + sp*ttmp(8)
  t(8) = -sp*ttmp(7) + cp*ttmp(8)

  return
end subroutine tensor_rot3

subroutine switch (a, b)
!     Switches real numbers.
  implicit none
  real(8) :: a,b,temp
  temp = a
  a = b
  b = temp

  return
end subroutine switch

subroutine switch3 (xca, xcb)
!     Interchanges two vectors
  implicit none
  real(8) :: xca(3), xcb(3), xtmp(3)

  xtmp(:) = xca(:)
  xca(:) = xcb(:)
  xcb(:) = xtmp(:)

  return
end subroutine switch3

subroutine cross (a,b,acrsb)

!     Calculates the cross product, acrsb, of vectors a and b.
  implicit none
  real(8) :: a(3), b(3), acrsb(3)

  acrsb(1)= a(2)*b(3) - a(3)*b(2)
  acrsb(2)= a(3)*b(1) - a(1)*b(3)
  acrsb(3)= a(1)*b(2) - a(2)*b(1)

  return
end subroutine cross

subroutine unit(x)
!     Makes a unit vector out of x.
  implicit none
  real(8) :: x(3)
  real(8) :: d
  d = sqrt( dot_product(x,x))

 if(d.gt.0.d0) then
   x(:) = x(:)/d
  else
   print *, ' ** Warning... zero unit vector in sub unit.'
   stop
  endif

  return
end subroutine unit


MODULE gridmod
  !-------------------------------------------------------------------
  !  Created     : Mateo Duque Villegas
  !  Last updated: 2-Feb-2018
  !-------------------------------------------------------------------
  IMPLICIT NONE

CONTAINS

    SUBROUTINE latlon(blat,blon,ydegr,xdegr,nlat,mlon,ylat,xlon)

    IMPLICIT NONE

    !! Global namespace
    INTEGER(KIND=4), INTENT(in) :: nlat
    INTEGER(KIND=4), INTENT(in) :: mlon
    REAL(KIND=4),    INTENT(in) :: ydegr
    REAL(KIND=4),    INTENT(in) :: xdegr
    REAL(KIND=4),    INTENT(in) :: blat
    REAL(KIND=4),    INTENT(in) :: blon
    REAL(KIND=8),    DIMENSION(nlat), INTENT(out)  :: ylat
    REAL(KIND=8),    DIMENSION(mlon), INTENT(out)  :: xlon

    !! Local namespace
    INTEGER(KIND=4)                                :: ilon
    INTEGER(KIND=4)                                :: ilat

    !! Longitudes going from -180 to 180
    DO ilon = 1,mlon
       xlon(ilon) = -blon + (ilon-1)*xdegr
    END DO

    !! Latitudes going from 90 to -90
    DO ilat = 1,nlat
       ylat(ilat) = blat - ilat*ydegr
    END DO

    RETURN

  END SUBROUTINE latlon

END MODULE gridmod


text:
	!----------------------------------------------------------------------------
    	subroutine grainsize(coords,ldistance,rdistance,lm,
     1		 slpdir1,slpnor1,nslptl,feature,oriensh,noel,
     2		 mask,mgrains,deadpoint)

	!subroutine containing the length scale calculation
	!INPUTS :
	!coords             - voxel coordinates
	!nodeout            - columns in nodex, nodey, nodez,boundfeat
	!totalfeat          - total number of grains in the RVE
	!nodex,nodey,nodez  - x, y, and z coordinates for each boundary node
	!slpdir1             - slip directions in local (crystal) coordinate system
	!slpnor1             - slip plane normals in local (crystal) coordinate system
	!elcent             - array containing the centroids of all voxels
	!feature            - grain number
	!boundgrain         - boundary grain ids
	!oriensh            - orientations in euler angles in each grain
	!noel               - voxel number

	!OUTPUTS:
	!rdistance          - distance from the voxel to grain boundary
	!ldistance          - distance along slip length
	!lm                 - Luster-Morris parameter

	parameter(nd=150)
	INCLUDE 'param_array.inc'

	integer*8:: i, arraysize,tester,grainb,
     1  arraysizeb,valb,NSLPTL,j,MININDEXSING,CHARCTERLENGT,
     2  intotalfeat, mgrains(2),minidexinsing
	integer :: val,noel

	real*8, allocatable ::  mindist(:),vect(:,:),
     1  sortmindist(:),sortedbgrains(:),uniqu(:),uniindex(:),
     2  uniqgraindist(:),zrot(:,:),slpdirrotate(:,:),
     3  grainboundindex(:) ,featureboundsnodes(:,:),
     4  bxvalue(:),byvalue(:),bzvalue(:),btotal(:,:),bangle(:,:),
     5  normarray(:),normslp(:),minangleval(:),minangleval180(:),
     6  minangleactual(:),rnodes(:,:),vectr(:,:),
     7  nodesnot(:,:),lxtotal(:,:,:),lxnormin(:,:),lxnorm(:),
     8  vectrnorm(:),lxangle(:,:),minxval0(:,:),minxval180(:,:),
     9  minxval0act(:,:),minxval180act(:,:),minxangle(:),
     &  minxangleact(:),lxnodes(:,:),lxnodesindex(:),
     1  slpdirrotatenorm(:),veca(:,:),vecb(:,:),checkangle(:),
     2  vanorm(:),vbnorm(:),slpplanrotate(:,:),
     3  slpdirrotatea(:,:),slpplanrotatea(:,:),slpplanrotatenorm(:),
     4  slpdirrotateanorm(:),slpplanrotateanorm(:),cosphi(:),cosk(:)

	integer, allocatable :: minindex(:),minang0(:,:),
     1 minang180(:,:),
     2  minangleactloc(:),boundnodegrainb(:),nodesnotindex(:)

	real*8 :: coords(3),vectotal(3),
     1 euler1,totalrot(3,3),groundbundid,pi,charctlenght,
     2 twovoxdist(3),centrodes(3),totalrota(3,3),
     3 nodex(totalfeat,nodeout), boundgrain(totalfeat,nodeout),
     4 nodey(totalfeat,nodeout), nodez(totalfeat,nodeout),
     5 elcent(totalels,3),oriensh(totalfeat,3),sortmindistsing,
     6 uniqgraindistsing,FEATURE,testera

	real*8 :: lm(nd),ldistance(nd),rdistance(nd)

	common nodex,nodey,nodez,boundgrain,elcent,intotalfeat

	real*8 :: SLPDIR1(3,nd), SLPNOR1(3,nd),SLPDIR(nd,3),
     1            slpplane(nd,3)

	LOGICAL, allocatable :: mk(:)

	logical :: locgrain= .False.,mask,deadpoint

	PI=4.D0*DATAN(1.D0)

	!adjust the slip plane normal and slip direction to removing
	 !any trailing zeros

	do i=1,nd
		do j=1,3
			slpdir(i,j)=slpdir1(j,i)
			slpplane(i,j)=slpnor1(j,i)
		end do
	end do

	! check to remove values of interest ** this for loop can be
	! replaced with 'findloc'**

	checkerloop: DO i=1,int(nodeout)

		if(int(nodey(int(feature),i)) == 9999999)then

			arraysize=i-1

			exit checkerloop
		else

			arraysize=i

		end if
	end do checkerloop

	!find the closest boundary.  This will form 'Grain B' from which
	!the information on orientation is extracted to rotate the slip
	!systems

	!create a vector from the current voxel position to boundary nodes
		allocate (vect(3,arraysize))

	vect(1,:)=nodex(int(feature),1:arraysize)-coords(1)

	vect(2,:)=nodey(int(feature),1:arraysize)-coords(2)

	vect(3,:)=nodez(int(feature),1:arraysize)-coords(3)

	vectotal=(/vect(1,:),vect(2,:),vect(3,:)/)

	allocate (mindist(arraysize))
	allocate (sortmindist(arraysize))
	allocate (mk(arraysize))
	allocate (minindex(arraysize))
	allocate (sortedbgrains(arraysize))

	! calculate the distance of the current location to the grain boundary and
	! multiplying to bring distances to microns
		mindist=(((vect(1,:)**2.D0) + (vect(2,:)**2.D0)+
     1              (vect(3,:)**2.D0))**0.5D0)

	! sort distance by the minimum size and and keep indices of ordered minimum

	minindexsing=minloc(mindist,1)
	sortmindistsing=minval(mindist)

	!using the index of unique grains, the corresponding distance value can be extracted
	allocate (uniqgraindist(size(uniindex)))
	uniqgraindistsing=int(boundgrain(int(feature),minindexsing))

	allocate( slpdirrotate(size(slpdir,1),3))
	allocate( slpplanrotate(size(slpplane,1),3))

	grainb=uniqgraindistsing

	!rotate the slip systems using the angles from each identified grain
	call eulercosmatrix(oriensh(grainb,:),totalrot)
	do j =1, size(slpdir,1)
		slpdirrotate(j,:)= matmul(totalrot, slpdir(j,:))
		slpplanrotate(j,:)=matmul(totalrot, slpplane(j,:))
	end do

	allocate( slpdirrotatea(size(slpdir,1),3))
	allocate( slpplanrotatea(size(slpplane,1),3))
	!rotate slip systems for the current grain
	call eulercosmatrix(oriensh(int(feature),:),totalrota)

	do j =1, size(slpdir,1)
		slpdirrotatea(j,:)= matmul(totalrota, slpdir(j,:))
		slpplanrotatea(j,:)=matmul(totalrota, slpplane(j,:))
	end do

	!using the identified closest grain, construct vectors extending from the
	!current voxel position to the grain boundary.

	!firstly extract only the nodes at the shared boundary between grain A and
	!grain B.

	!firtly find the size of the array at the indentified grain B
	! check to remove values of interest **this for loop can be replaced with
	! 'findloc' if it is available**

	checkerloopb: DO i=1,int(nodeout)
		if(nodex(int(grainb),i) ==9999999)then
			arraysizeb=i-1
			exit checkerloopb
		else
			arraysizeb=i

		 end if
	end do checkerloopb

	!************************************************************
	!find index of array locations which are true.
	!For loops are used here to make it more generically applicable
	!but you can also use 'all' first to create a boolean array
	! and the 'findloc' to find .TRUE. values.

	val=1
	do i=1,arraysizeb
		if(int(boundgrain(grainb,i)) .eq. int(feature))then
			val=val+1
		end if
	end do

	allocate (grainboundindex(val-1))
	val=1
	do i=1,arraysizeb
		if(int(boundgrain(grainb,i)) .eq. int(feature))then
			grainboundindex(val)=i
			val=val+1
		end if
	end do

	   !**********************************************************

	   !find the coordinates of the identified shared boundary nodes
	   allocate (featureboundsnodes(val-1,3))
	   featureboundsnodes=reshape((/nodex(grainb,grainboundindex),
     1   	 nodey(grainb,grainboundindex),
     2    	nodez(grainb,grainboundindex)/),(/val-1,3/))

	   !calculate the vector from the boundary nodes to the current
	   !voxel
	   allocate (bxvalue(val-1))
	   allocate (byvalue(val-1))
	   allocate (bzvalue(val-1))
	   bxvalue=featureboundsnodes(:,1)-coords(1)
	   byvalue=featureboundsnodes(:,2)-coords(2)
	   bzvalue=featureboundsnodes(:,3)-coords(3)

	   allocate (btotal(val-1,3))
	   btotal=reshape((/bxvalue,byvalue,bzvalue/),(/val-1,3/))

	   !find the closest angle between the vectors from the interface
	   !boundary to the voxel and the slip direction in grain B.

	   allocate (bangle(size(slpdirrotate,1),val-1))

	   !these subroutine calls can just be replaced with the
	   !intrisic function norm2*********
	   allocate (normarray(val-1))
	   allocate (normslp(size(slpdirrotate,1)))

	   normarray=norm2(btotal,dim=2)
	   !call enorm(btotal,val-1,normarray)

	   !call enorm(slpdirrotate,size(slpdirrotate,1),normslp)
	   normslp=norm2(slpdirrotate,dim=2)

	   !****************************************************
	   !determine if the current voxel is on the boundary.
	   !If this is true, this voxel should form the grain boundary voxel
	   !To do this, determine if any euclidean distances are less than
	   !the elements characteristic length

	   !calculate characteristic length of voxel

	   twovoxdist=(elcent(3,:)-elcent(2,:))/2.D0
	   charcterlengt=(twovoxdist(1)**2.D0 + twovoxdist(2)**2.D0 +
     1   twovoxdist(3)**2.D0)**0.5D0

	  !this can be done with findloc if available
	  do i=1,size(normarray,1)
		if(normarray(i) .lt. charcterlengt)then
			locgrain=.True.
		 end if

	  end do

	   allocate (minang0(size(slpdirrotate,1),1))
	   allocate (minangleval(size(slpdirrotate,1)))
	   allocate (minang180(size(slpdirrotate,1),1))
	   allocate (minangleval180(size(slpdirrotate,1)))
	   allocate (minangleactual(size(slpdirrotate,1)))
	   allocate (minangleactloc(size(slpdirrotate,1)))
	   allocate (boundnodegrainb(size(slpdirrotate,1)))
	   allocate (rnodes(size(slpdirrotate,1),3))

	   if (locgrain == .FALSE.) then

	   !voxel not on the boundary

		   do i=1,size(slpdirrotate,1)
		   !calculate angle between the slipdirection in grian B and the vector
		   !created between the current voxel and the boundary nodes.
				bangle(i,:)=dacos(matmul(btotal,slpdirrotate(i,:))/
     1       			(normarray*normslp(i)))

		   !find minimum angle and the location index
		   !of the array for each slip direction
				minang0(i,:)=minloc(bangle(i,:))

				minangleval(i)=bangle(i,minang0(i,1))

			!find the minimums at angles of 180 degrees
				minang180(i,:)=minloc(abs(bangle(i,:)-pi))
				minangleval180(i)=abs(bangle(i,minang180(i,1))-pi)
			!find the smaller of the 0 and 180 degrees to determine
			!the true mininum value
				if (minangleval(i) .lt. minangleval180(i)) then
					minangleactual(i)=minangleval(i)
					minangleactloc(i)=minang0(i,1)
				else
					minangleactual(i)=minangleval180(i)
					minangleactloc(i)=minang180(i,1)
				end if
		   end do

		   !using the minimum angles, the nodes at the boundary  which connect to the voxel location
		   !can be determined

		   boundnodegrainb=grainboundindex(minangleactloc)
		   rnodes=reshape((/nodex(grainb,int(boundnodegrainb)),
     1       		nodey(grainb,int(boundnodegrainb)),
     2      		 nodez(grainb,int(boundnodegrainb))/),
     3      		(/size(slpdirrotate,1),3/))

		   ! the corresponding vector for each slip direction
		   allocate (vectr(size(slpdirrotate,1),3))

		   vectr=reshape((/nodex(grainb,int(boundnodegrainb))-coords(1),
     1              	nodey(grainb,int(boundnodegrainb))-coords(2),
     2              	nodez(grainb,int(boundnodegrainb))-coords(3)
     3              	/),(/size(slpdirrotate,1),3/))

		   ! calculate the distance from the voxel centroid to the boundary nodes
		   rdistance=norm2(vectr,dim=2)
		   !call enorm(vectr,size(slpdirrotate,1),rdistance)

	   else

	   !since the voxel is on the boundary, the distance to the boundary is taken
	   !as half the voxel characteristic length

		   !rdistance(1:size(slpdirrotate,1))=charcterlengt*0.5D0
		   rdistance=0.5D0

		   do i=1,3
			  rnodes(1:size(slpdirrotate,1),i)=elcent(int(noel),i)
		   end do
	   end if

	   ! calculate the distance from the boundary to boundary in grain B along
	   ! the slip direction

	   !****************************************************************
	   !This is to be done by excluding the nodes at the boundary.
	   !For loops are used here to make it more generically applicable
	   !but you can also use 'all' first to create a boolean array
	   ! and the 'findloc' to find .TRUE. values.

		valb=1
		do i=1,arraysizeb
			if(int(boundgrain(grainb,i)) .ne. int(feature))then
				valb=valb+1
			end if
		end do

		allocate (nodesnotindex(valb-1))
		valb=1
		do i=1,arraysizeb
			if(int(boundgrain(grainb,i)) .ne. int(feature))then
				nodesnotindex(valb)=i
				valb=valb+1
			end if
		end do

		!**********************************************************
		!node coordinates not at the shared boundary.
		allocate (nodesnot(valb-1,3))

		nodesnot=reshape((/nodex(grainb,int(nodesnotindex)),
     1            	nodey(grainb,int(nodesnotindex)),
     2           	 nodez(grainb,int(nodesnotindex))
     3            	/),(/valb-1,3/))

		allocate (lxtotal(size(slpdirrotate,1),valb-1,3))
		allocate (lxnorm(valb-1))
		allocate (lxnormin(valb-1,3))
		allocate (lxangle(size(slpdirrotate,1),valb-1))
		allocate (minxval0(size(slpdirrotate,1),1))
		allocate (minxval180(size(slpdirrotate,1),1))
		allocate (minxval180act(size(slpdirrotate,1),1))
		allocate (minxval0act(size(slpdirrotate,1),1))
		allocate (minxangle(size(slpdirrotate,1)))
		allocate (minxangleact(size(slpdirrotate,1)))

		!calculate vector from the boundary to the nodes in Grain B not
		!at the boundary

		! check to see if the current voxel is at the boundary or not
		! If not, progress in the following
		allocate (vectrnorm(size(slpdirrotate,1)))
		allocate (slpdirrotatenorm(size(slpdirrotate,1)))
		if (locgrain == .FALSE.)then

			!call enorm(vectr,size(slpdirrotate,1),vectrnorm)
			vectrnorm=norm2(vectr,dim=2)
			do i=1,size(slpdirrotate,1)
				lxtotal(i,:,:)=reshape((/nodex(grainb,int(nodesnotindex))
     1            			-rnodes(i,1),
     2            			nodey(grainb,int(nodesnotindex))-rnodes(i,2),
     3           			 nodez(grainb,int(nodesnotindex))-rnodes(i,3)
     4            			/),(/valb-1,3/))
				lxnormin=lxtotal(i,:,:)
				!******can be replaced with norm2 function
				!call enorm(lxnormin,valb-1,lxnorm)
				lxnorm=norm2(lxnormin,dim=2)

				!find the angle between the vectors from the shared boundary to
				!boundary nodes in grain b and the vector developed from the voxel
				!to the shared boundary
				lxangle(i,:)= dacos(matmul(lxtotal(i,:,:),vectr(i,:))/(lxnorm*
     1         				vectrnorm(i)))
				minxval0(i,:)=minloc(lxangle(i,:))
				minxval0act(i,:)=lxangle(i,int(minxval0(i,:)))
				minxval180(i,:)=minloc(abs(lxangle(i,:)-pi))
				minxval180act(i,:)=abs(pi-lxangle(i,int(minxval180(i,:))))

				minxangle(i)=minxval0(i,1)
				minxangleact(i)=minxval0act(i,1)

			 end do
		 else

			!call enorm(slpdirrotate,size(slpdirrotate,1),slpdirrotatenorm)
			slpdirrotatenorm=norm2(slpdirrotate,dim=2)

			do i=1,size(slpdirrotate,1)
				lxtotal(i,:,:)=reshape((/nodex(grainb,int(nodesnotindex))
     1         				 -rnodes(i,1),
     2          			 nodey(grainb,int(nodesnotindex))-rnodes(i,2),
     3         				 nodez(grainb,int(nodesnotindex))-rnodes(i,3)
     4          			/),(/valb-1,3/))
				lxnormin=lxtotal(i,:,:)
				!******can be replaced with norm2 function
				!call enorm(lxnormin,valb-1,lxnorm)
				lxnorm=norm2(lxnormin,dim=2)
				!find the angle between the vectors from the shared boundary to
				!boundary nodes in grain b and the vector developed from the voxel
				!to the shared boundary
				lxangle(i,:)= dacos(matmul(lxtotal(i,:,:),slpdirrotate(i,:))/(lxnorm*
     1       				      slpdirrotatenorm(i)))
				minxval0(i,:)=minloc(lxangle(i,:))
				minxval0act(i,:)=lxangle(i,int(minxval0(i,:)))
				minxval180(i,:)=minloc(abs(lxangle(i,:)-pi))
				minxval180act(i,:)=abs(pi-lxangle(i,int(minxval180(i,:))))

				!since the vectors in this case case be parallel but in opposite directions
				!an 'if' loop is used all sort through the 0deg and 180deg possible
				!minimums
				if(minxval0act(i,1)<minxval180act(i,1))then
					minxangle(i)=minxval0(i,1)
					minxangleact(i)=minxval0act(i,1)
				else
					minxangle(i)=minxval180(i,1)
					minxangleact(i)=minxval180act(i,1)
				end if

			 end do
		 end if

		 allocate (lxnodesindex(size(slpdirrotate,1)))
		 allocate (lxnodes(size(slpdirrotate,1),3))
		 !using the location of the closest vector in grain B to the developed vector
		 !along grain B's slip systems in grain A, the vector and corresponding nodes
		 !can be determined.
		 lxnodesindex=nodesnotindex(minxangle)
		 lxnodes=reshape((/nodex(grainb,lxnodesindex),
     1  		nodey(grainb,lxnodesindex),
     2  		nodez(grainb,lxnodesindex)/),(/size(slpdirrotate,1),3/))

		 !need to check the angle formed by the lxnodes vector and the slip vector
		 !in grain A.  If the difference is large (greater than a specified minimum
		 !allowable angle), the slip distance in grian B should be ignored since it
		 !is unreliable
		 allocate (vecb(size(slpdirrotate,1),3))
		 allocate (veca(size(slpdirrotate,1),3))
		 if (locgrain == .FALSE.)then
			!vector from shared bounday nodes to
			vecb=lxnodes-rnodes
			!vector from voxel to shared boundary
			do i=1,size(slpdirrotate,1)
				veca(i,:)=elcent(int(noel),:)-rnodes(i,:)
			end do
		  else
			!vector from shared bounday nodes to
			vecb=lxnodes-rnodes
			!voxel at the bounary so use slip direction in grain B
			veca=slpdirrotate
		  end if

		  !calculate the angle between vector a and b
		  allocate (vanorm(size(slpdirrotate,1)))
		  allocate (vbnorm(size(slpdirrotate,1)))
		  allocate (checkangle(size(slpdirrotate,1)))
		  !call enorm(veca,size(slpdirrotate,1),vanorm)
		  vanorm=norm2(veca,dim=2)
		  !call enorm(vecb,size(slpdirrotate,1),vbnorm)
		   vbnorm=norm2(vecb,dim=2)
                  ldistance=vbnorm
!        -------------------------------------------------------------------------
!                 This section is added to check the angle to see if it is
!                 within a tolerance. This can be updated further by the user
!        -------------------------------------------------------------------------
!		  do i=1,size(slpdirrotate,1)
!			checkangle(i)=dacos(dot_product(veca(i,:),vecb(i,:))
!     1      			 /(vanorm(i)*vbnorm(i)))
		  !check to see if the angle is below the recommended minimum
!			if(checkangle(i)<abs(pi-(40.D0*pi/180.D0)))then
			   !lxnodes(i,:)=rnodes(i,:)
		   !distance is taken only as half the characteristic voxel length

			   !ldistance(i)=charcterlengt*0.5D0
!			 else
!			   ldistance(i)=vbnorm(i)
!			 end if
!		  end do
!-----------------------------------------------------------------------------------

		  !calculate the Luster-Morris parameter
		  allocate (slpplanrotatenorm(size(slpplanrotate,1)))
		  allocate (slpdirrotateanorm(size(slpdir,1)))
		  allocate (slpplanrotateanorm(size(slpplane,1)))

		  !call enorm(slpplanrotate,size(slpplanrotate,1),slpplanrotatenorm)
		  slpplanrotatenorm=norm2(slpplanrotate,dim=2)
		  !call enorm(slpdirrotatea,size(slpdir,1),slpdirrotateanorm)
		   slpdirrotateanorm=norm2(slpdirrotatea,dim=2)
		  !call enorm(slpplanrotatea,size(slpplane,1),slpplanrotateanorm)
		   slpplanrotateanorm=norm2(slpplanrotatea,dim=2)

		  allocate (cosphi(size(slpplane,1)))
		  allocate (cosk(size(slpplane,1)))

		  do i=1,size(slpplane,1)
			cosphi(i)=dot_product(slpplanrotate(i,:),slpplanrotatea(i,:))/
     1        			(slpplanrotatenorm(i)*slpplanrotateanorm(i))

			cosk(i)=dot_product(slpdirrotatea(i,:),slpdirrotate(i,:))/
     1       			 (slpdirrotateanorm(i)*normslp(i))
			lm(i)=abs(cosphi(i))*abs(cosk(i))
		  end do

		  !need to make sure to deallocate arrays
		 deallocate (vect,mindist  ,sortmindist, mk, minindex,sortedbgrains,
     1       		uniqgraindist,slpdirrotate,slpplanrotate,slpdirrotatea,
     2      		slpplanrotatea,grainboundindex,featureboundsnodes,
     3       		bxvalue,byvalue,bzvalue,btotal,bangle,normarray,normslp,
     4       		minangleval,minang180,minangleval180,minang0,
     5        		minangleactual,minangleactloc,boundnodegrainb,rnodes,
     6        		nodesnotindex,nodesnot,lxtotal,lxnorm,lxangle,lxnormin,
     7        		minxval0,minxval180,minxval180act,minxval0act,minxangle,
     8        		minxangleact,vectrnorm,slpdirrotatenorm,lxnodesindex,
     9        		lxnodes,vecb,veca,vanorm,vbnorm,checkangle,
     &        		slpplanrotatenorm,slpdirrotateanorm,slpplanrotateanorm,
     1        		cosphi,cosk)

	end subroutine grainsize

!-------------------------------------------------------------------
!        Alternative subourinte for norm2
!        ==================================
	subroutine enorm(array,val,normarray)
	real*8, allocatable :: array(:,:)
	integer*8, intent(in) :: val
	real*8 :: normarray(val)

	normarray=(array(:,1)**2.0 + array(:,2)**2.0 + array(:,3)**2.0)**0.50

	end subroutine enorm

!----------------------------------------------------------------------
!        Construct the rotation matrix based on Euler angles
!        ===================================================
	subroutine eulercosmatrix(orient,totalrot)
		Real*8 :: euler1,euler2,euler3,pi,zrot(3,3),xrot(3,3),zrot2(3,3),
     1  	 orient(1,3)
		real*8 :: totalrot(3,3)
		integer :: i, l,j

		PI=4.D0*DATAN(1.D0)
		euler1=orient(1,1)*PI/180.D0
		euler2=orient(1,2)*PI/180.D0
		euler3=orient(1,3)*PI/180.D0

		zrot=reshape((/cosd(euler1),-sind(euler1),0.d0,
     1   		sind(euler1),cosd(euler1),0.d0,
     2      		 0.d0,0.d0,1.d0/),(/3,3/))

		xrot=reshape((/1.d0,0.d0,0.d0,
     1   		0.d0,cosd(euler2),-sind(euler2),
     2  		 0.d0, sind(euler2),cosd(euler2)/),(/3,3/))

		zrot2=reshape((/cosd(euler3),-sind(euler3),0.d0,
     1    		sind(euler3),cosd(euler3),0.d0,
     2          	 0.d0,0.d0,1.d0/),(/3,3/))

		!total rotation matrix
		totalrot=transpose(matmul(matmul(zrot2,xrot),zrot))

		end subroutine eulercosmatrix

!----------------------------------------------------------------------

ext: for

max_stars_repo_path: Additional Files/Length scale subroutines/grainsize_subroutine.for

max_stars_repo_name: DylanAgius/LengMorph

max_issues_repo_name: DylanAgius/LengMorph

max_forks_repo_name: DylanAgius/LengMorph


text:
      SUBROUTINE sprstp(sa,ija,sb,ijb)
      INTEGER ija(*),ijb(*)
      REAL sa(*),sb(*)
CU    USES iindexx
      INTEGER j,jl,jm,jp,ju,k,m,n2,noff,inc,iv
      REAL v
      n2=ija(1)
      do 11 j=1,n2-2
        sb(j)=sa(j)
11    continue
      call iindexx(ija(n2-1)-ija(1),ija(n2),ijb(n2))
      jp=0
      do 13 k=ija(1),ija(n2-1)-1
        m=ijb(k)+n2-1
        sb(k)=sa(m)
        do 12 j=jp+1,ija(m)
          ijb(j)=k
12      continue
        jp=ija(m)
        jl=1
        ju=n2-1
5       if (ju-jl.gt.1) then
          jm=(ju+jl)/2
          if(ija(jm).gt.m)then
            ju=jm
          else
            jl=jm
          endif
          goto 5
        endif
        ijb(k)=jl
13    continue
      do 14 j=jp+1,n2-1
        ijb(j)=ija(n2-1)
14    continue
      do 16 j=1,n2-2
        jl=ijb(j+1)-ijb(j)
        noff=ijb(j)-1
        inc=1
1       inc=3*inc+1
        if(inc.le.jl)goto 1
2       continue
          inc=inc/3
          do 15 k=noff+inc+1,noff+jl
            iv=ijb(k)
            v=sb(k)
            m=k
3           if(ijb(m-inc).gt.iv)then
              ijb(m)=ijb(m-inc)
              sb(m)=sb(m-inc)
              m=m-inc
              if(m-noff.le.inc)goto 4
            goto 3
            endif
4           ijb(m)=iv
            sb(m)=v
15        continue
        if(inc.gt.1)goto 2
16    continue
      return
      END

max_stars_repo_path: NR-Functions/Numerical Recipes- Example & Functions/Functions/sprstp.for