text:
      SUBROUTINE BMG3_SymStd_SETUP_interp_OI( &
     &                KGF, KGC, SO, SOC, CI,&
     &                IIF, JJF, KKF, IIC, JJC, KKC, &
     &                NOG, ifd, NStncl, irelax, yo, &
     &                NOGm, IGRD, JPN, halof&
     &                ) BIND(C, NAME='MPI_BMG3_SymStd_SETUP_interp_OI')

! ======================================================================
!  --------------------
!   DESCRIPTION:
!  --------------------
!
!     BMG3_SymStd_SETUP_interp_OI constructs the operator-induced
!     interpolation operator CI from the fine-grid stencil, SO.  The
!     operator CI interpolates a vector from the coarse grid KC, to the
!     fine grid, KF.
!
! ======================================================================
! $license_flag$
! ======================================================================
!  --------------------
!   INPUT:
!  --------------------
!
!
! ======================================================================
!  --------------------
!   INPUT/OUTPUT:
!  --------------------
!
!
! ======================================================================
!  --------------------
!   OUTPUT:
!  --------------------
!
!
!
! ======================================================================
!  --------------------
!   LOCAL:
!  --------------------
!
!
! ======================================================================

      USE ModInterface
      IMPLICIT NONE

! -----------------------------
!     Includes
!
      INCLUDE 'mpif.h'
      INCLUDE 'MSG_f90.h'

      INCLUDE 'BMG_constants_f90.h'
      INCLUDE 'BMG_stencils_f90.h'
      INCLUDE 'BMG_workspace_f90.h'
      INCLUDE 'BMG_parameters_f90.h'

! ---------------------------
!     Argument Declarations:
!
      integer(c_int), value :: NOG, NOGm, NStncl, KGF, KGC, ifd, irelax
      integer(c_int), value :: JPN
      integer(len_t), value :: IIF, JJF, KKF, IIC, JJC, KKC
      integer(len_t) :: IGRD(NOGm,NBMG_pIGRD)
      real(real_t) :: CI(IIC,JJC,KKC,26), SO(IIF+1,JJF+1,KKF+1,NStncl)
      real(real_t) :: SOC(IIC+1,JJC+1,KKC+1,14), YO(IIF,JJF,2,14)
      type(c_ptr) :: halof
! --------------------------
!     Local Declarations:
!
      INTEGER i, j, k, ic, jc, kc, IIC1, JJC1, KKC1, KPZ

      REAL*8  a, b, c, de, dn, dne, dnw, dp, ds, dse, dsw, dw,&
     &        eMACH, ep, sum
      real*8 D1MACH
      LOGICAL LEFTPLANE, BOTTOMPLANE, FRONTPLANE,&
     &        RIGHTPLANE, TOPPLANE, BACKPLANE
      INTEGER ISTART, JSTART, KSTART, ICSTART, JCSTART, KCSTART, &
     &        ISTARTO, JSTARTO, KSTARTO, ICSTARTO, JCSTARTO,  &
     &        KCSTARTO, ICEND, JCEND, KCEND, ICENDO, JCENDO, KCENDO&
     &

      INTEGER ierror, ptrn, iGs_c, jGs_c, kGs_c
      LOGICAL  PERIODIC_X, PERIODIC_Y, PERIODIC_Z

! ======================================================================

! --------------------------------------------
!     Test Periodicity
! --------------------------------------------

      PERIODIC_X = (JPN.EQ.BMG_BCs_def_per_x     .OR.&
                    JPN.EQ.BMG_BCs_def_per_xy    .OR.&
                    JPN.EQ.BMG_BCs_def_per_xz    .OR.&
                    JPN.EQ.BMG_BCs_def_per_xyz   .OR.&
                    JPN.EQ.BMG_BCs_indef_per_x   .OR.&
                    JPN.EQ.BMG_BCs_indef_per_xy  .OR.&
                    JPN.EQ.BMG_BCs_indef_per_xz  .OR.&
                    JPN.EQ.BMG_BCs_indef_per_xyz)

      PERIODIC_Y = (JPN.EQ.BMG_BCs_def_per_y    .OR.&
                    JPN.EQ.BMG_BCs_def_per_xy   .OR.&
                    JPN.EQ.BMG_BCs_def_per_yz   .OR.&
                    JPN.EQ.BMG_BCs_def_per_xyz  .OR.&
                    JPN.EQ.BMG_BCs_indef_per_y  .OR.&
                    JPN.EQ.BMG_BCs_indef_per_xy .OR.&
                    JPN.EQ.BMG_BCs_indef_per_yz .OR.&
                    JPN.EQ.BMG_BCs_indef_per_xyz)

      PERIODIC_Z = (JPN.EQ.BMG_BCs_def_per_z     .OR.&
                    JPN.EQ.BMG_BCs_def_per_xz    .OR.&
                    JPN.EQ.BMG_BCs_def_per_yz    .OR.&
                    JPN.EQ.BMG_BCs_def_per_xyz   .OR.&
                    JPN.EQ.BMG_BCs_indef_per_z   .OR.&
                    JPN.EQ.BMG_BCs_indef_per_xz  .OR.&
                    JPN.EQ.BMG_BCs_indef_per_yz  .OR.&
                    JPN.EQ.BMG_BCs_indef_per_xyz)

      eMACH = D1MACH(3)

      IIC1 = IIC-1
      JJC1 = JJC-1
      KKC1 = KKC-1

      DO KPZ=1,14
         DO k=1,2
            DO j=1,JJF
               DO i=1,IIF
                  yo(i,j,k,KPZ)= rZERO
               ENDDO
            ENDDO
         ENDDO
      ENDDO

! ------------------------------------------------------------------
!
!     Loop boundaries for interpolation setup:
!
!     The grid is assumed to be logically rectangular, thus all
!     indexing is simply tensor product.  This simplifies matters
!     in that we can consider each dimension independently.
!
!
!     Fine-Grid Indexing:
!
!     (i,j,k)   are used to index fine-grid points
!
!     ISTART
!     ISTARTO
!
!     Coarse-Grid Indexing:
!
!     (ic,jc,kc) are used to index coarse-grid points
!
!     ICSTART
!     ICSTARTO
!
!     ICEND       = iic - 1
!     ICENDO      = iicf - 1
!
! ------------------------------------------------------------------

      iGs_c = IGRD(KGC,idL_BMG_ICoord)
      jGs_c = IGRD(KGC,idL_BMG_JCoord)
      kGs_c = IGRD(KGC,idL_BMG_KCoord)

      LEFTPLANE   = (IGRD(KGF,idL_BMG_ICoord).EQ.1 .and. .not. PERIODIC_X)
      BOTTOMPLANE = (IGRD(KGF,idL_BMG_JCoord).EQ.1 .and. .not. PERIODIC_Y)
      FRONTPLANE  = (IGRD(KGF,idL_BMG_KCoord).EQ.1 .and. .not. PERIODIC_Z)

      RIGHTPLANE  = ((IGRD(KGF,idL_BMG_ICoord)+IGRD(KGF,idL_BMG_NLx)-1)&
     &     .EQ.IGRD(KGF,idL_BMG_NGx) .and. .not. PERIODIC_X)
      TOPPLANE    = ((IGRD(KGF,idL_BMG_JCoord)+IGRD(KGF,idL_BMG_NLy)-1)&
     &     .EQ.IGRD(KGF,idL_BMG_NGy) .and. .not. PERIODIC_Y)
      BACKPLANE   = ((IGRD(KGF,idL_BMG_KCoord)+IGRD(KGF,idL_BMG_NLz)-1)&
     &     .EQ.IGRD(KGF,idL_BMG_NGz) .and. .not. PERIODIC_Z)

      IF (LEFTPLANE) THEN
         ISTARTO  = 2
         ICSTARTO = 3
      ELSE
         IF (mod(IGRD(KGF,idL_BMG_ICoord),2).EQ.1) THEN
            ISTARTO = 0
         ELSE
            ISTARTO = 1
         ENDIF

         ICSTARTO = 2
      ENDIF

      IF (RIGHTPLANE.and.&
     &     ( (IGRD(KGC,idL_BMG_ICoord)+IIC-3)*2 .EQ.&
     &     (IGRD(KGF,idL_BMG_ICoord)+IIF-2) ) ) THEN
         ICENDO = IIC1
      ELSE
         ICENDO = IIC1+1
      ENDIF

      IF (mod(IGRD(KGF,idL_BMG_ICoord),2).EQ.1) THEN
         ISTART = 0
      ELSE
         ISTART = 1
      ENDIF
      ICSTART = 2
      ICEND = IIC1

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

      IF (BOTTOMPLANE) THEN
         JSTARTO  = 2
         JCSTARTO = 3
      ELSE
         IF (mod(IGRD(KGF,idL_BMG_JCoord),2).EQ.1) THEN
            JSTARTO = 0
         ELSE
            JSTARTO = 1
         ENDIF

         JCSTARTO = 2
      ENDIF

      IF (TOPPLANE.and.&
     &     ( (IGRD(KGC,idL_BMG_JCoord)+JJC-3)*2 .EQ.&
     &     (IGRD(KGF,idL_BMG_JCoord)+JJF-2) ) ) THEN
         JCENDO = JJC1
      ELSE
         JCENDO = JJC1+1
      ENDIF

      IF (mod(IGRD(KGF,idL_BMG_JCoord),2).EQ.1) THEN
         JSTART = 0
      ELSE
         JSTART = 1
      ENDIF
      JCSTART = 2
      JCEND = JJC1

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

      IF (FRONTPLANE) THEN
         KSTARTO  = 2
         KCSTARTO = 3
      ELSE
         IF (mod(IGRD(KGF,idL_BMG_KCoord),2).EQ.1) THEN
            KSTARTO = 0
         ELSE
            KSTARTO = 1
         ENDIF

         KCSTARTO = 2
      ENDIF

      IF (BACKPLANE.and.&
     &     ( (IGRD(KGC,idL_BMG_KCoord)+KKC-3)*2 .EQ.&
     &     (IGRD(KGF,idL_BMG_KCoord)+KKF-2) ) ) THEN
         KCENDO = KKC1
      ELSE
         KCENDO = KKC1+1
      ENDIF

      IF (mod(IGRD(KGF,idL_BMG_KCoord),2).EQ.1) THEN
         KSTART = 0
      ELSE
         KSTART = 1
      ENDIF
      KCSTART = 2
      KCEND = KKC1

! --------------------------------------------------------------------
! >>>>>>>>>>>>>>>>>> BEGIN:  INTERPOLATION CONSTRUCTION <<<<<<<<<<<<<<
! --------------------------------------------------------------------

      IF ( kgf.LT.NOG .OR. ifd.NE.1 ) THEN

         !
         !  Interpolation based on a 27-point stencil
         !

         !
         !  (1D) Piecewise linear interpolation along coarse grid lines
         !

         !
         !  Interpolate fine-grid points lieing on coarse x-lines
         !
         !  NB: storage viewed as points lieing on fine-only
         !      y-lines of a coarse xy-plane (i.e., coarse k-plane)
         !
         k=KSTART
         do kc=KCSTART,KCEND
            k=k+2
            j=JSTART
            do jc=JCSTART,JCEND
               j=j+2
               i=ISTARTO
               do ic=ICSTARTO,ICENDO
                  i=i+2
                  a=SO(i-1,j+1,k,kpnw)+SO(i-1,j,k,kpw)+SO(i-1,j,k,kpsw)&
     &                 +SO(i-1,j+1,k,kbnw)+SO(i-1,j,k,kbw)&
     &                 +SO(i-1,j,k,kbsw)+SO(i-1,j+1,k+1,kbse)&
     &                 +SO(i-1,j,k+1,kbe)+SO(i-1,j,k+1,kbne)
                  b=SO(i,j+1,k,kpsw)+SO(i,j,k,kpw)+SO(i,j,k,kpnw)&
     &                 +SO(i,j+1,k,kbne)+SO(i,j,k,kbe)+SO(i,j,k,kbse)&
     &                 +SO(i,j+1,k+1,kbsw)+SO(i,j,k+1,kbw)&
     &                 +SO(i,j,k+1,kbnw)
                  c=a+b+SO(i-1,j,k,kps)+SO(i-1,j+1,k,kps)&
     &                 +SO(i-1,j+1,k,kbn)+SO(i-1,j,k,kb)+SO(i-1,j,k,kbs)&
     &                 +SO(i-1,j+1,k+1,kbs)+SO(i-1,j,k+1,kb)&
     &                 +SO(i-1,j,k+1,kbn)
                  ep=MIN(abs(a/so(i-1,j,k,kp)),abs(b/so(i-1,j,k,kp)))
                  c=a+b+(SO(i-1,j,k,kp)-c)*MAX(SO(i-1,j,k,kp)&
     &                 -(rONE+ep)*c,rZERO)/(abs(SO(i-1,j,k,kp)&
     &                 -(rONE+ep)*c)+eMACH)
                  CI(ic,jc,kc,lxyl)=a/c
                  CI(ic,jc,kc,lxyr)=b/c
               enddo
            enddo
         enddo

         !
         !  Interpolate fine-grid points lieing on coarse y-lines
         !
         !  NB: storage viewed as points lieing on fine-only
         !      x-lines of a coarse xy-plane (i.e., coarse k-plane)
         !
         k=KSTART
         do kc=KCSTART,KCEND
            k=k+2
            j=JSTARTO
            do jc=JCSTARTO,JCENDO
               j=j+2
               i=ISTART
               do ic=ICSTART,ICEND
                  i=i+2
                  a=SO(i,j,k,kpnw)+SO(i,j,k,kps)+SO(i+1,j,k,kpsw)&
     &                 +SO(i,j,k,kbnw)+SO(i,j,k,kbn)+SO(i+1,j,k,kbne)&
     &                 +SO(i,j,k+1,kbse)+SO(i,j,k+1,kbs)&
     &                 +SO(i+1,j,k+1,kbsw)
                  b=SO(i,j-1,k,kpsw)+SO(i,j-1,k,kps)+SO(i+1,j-1,k,kpnw)&
     &                 +SO(i,j-1,k,kbsw)+SO(i,j-1,k,kbs)&
     &                 +SO(i+1,j-1,k,kbse)+SO(i,j-1,k+1,kbne)&
     &                 +SO(i,j-1,k+1,kbn)+SO(i+1,j-1,k+1,kbnw)
                  ep=MIN(abs(a/so(i,j-1,k,kp)),abs(b/so(i,j-1,k,kp)))
                  c=a+b+SO(i,j-1,k,kpw)+SO(i+1,j-1,k,kpw)&
     &                 +SO(i,j-1,k,kbw)+SO(i,j-1,k,kb)+SO(i+1,j-1,k,kbe)&
     &                 +SO(i,j-1,k+1,kbe)+SO(i,j-1,k+1,kb)&
     &                 +SO(i+1,j-1,k+1,kbw)
                  c=a+b+(SO(i,j-1,k,kp)-c)*MAX(SO(i,j-1,k,kp)&
     &                 -(rONE+ep)*c,rZERO)/(abs(SO(i,j-1,k,kp)&
     &                 -(rONE+ep)*c)+eMACH)
                  CI(ic,jc,kc,lxya)=a/c
                  CI(ic,jc,kc,lxyb)=b/c
               enddo
            enddo
         enddo

         !
         !  Interpolate fine-grid points lieing on coarse z-lines
         !
         !  NB: storage viewed as points lieing on fine-only
         !      x-lines of a coarse xz-plane (i.e., coarse j-plane)
         !
         k=KSTARTO
         do kc=KCSTARTO,KCENDO
            k=k+2
            j=JSTART
            do jc=JCSTART,JCEND
               j=j+2
               i=ISTART
               do ic=ICSTART,ICEND
                  i=i+2
                  a=SO(i,j+1,k,kbse)+SO(i,j+1,k,kbs)+SO(i+1,j+1,k,kbsw)&
     &                 +SO(i,j,k,kbe)+SO(i,j,k,kb)+SO(i+1,j,k,kbw)&
     &                 +SO(i,j,k,kbne)+SO(i,j,k,kbn)+SO(i+1,j,k,kbnw)
                  b=SO(i,j+1,k-1,kbnw)+SO(i,j+1,k-1,kbn)&
     &                 +SO(i+1,j+1,k-1,kbne)+SO(i,j,k-1,kbw)&
     &                 +SO(i,j,k-1,kb)+SO(i+1,j,k-1,kbe)&
     &                 +SO(i,j,k-1,kbsw)+SO(i,j,k-1,kbs)&
     &                 +SO(i+1,j,k-1,kbse)
                  c=a+b+SO(i,j,k-1,kpw)+SO(i+1,j,k-1,kpw)&
     &                 +SO(i,j+1,k-1,kpnw)+SO(i,j+1,k-1,kps)&
     &                 +SO(i+1,j+1,k-1,kpsw)+SO(i,j,k-1,kpsw)&
     &                 +SO(i,j,k-1,kps)+SO(i+1,j,k-1,kpnw)
                  ep=MIN(abs(a/so(i,j,k-1,kp)),abs(b/so(i,j,k-1,kp)))
                  c=a+b+(SO(i,j,k-1,kp)-c)*MAX(SO(i,j,k-1,kp)&
     &                 -(rONE+ep)*c,rZERO)/(abs(SO(i,j,k-1,kp)&
     &                 -(rONE+ep)*c)+eMACH)
                  CI(ic,jc,kc,lxza)=a/c
                  CI(ic,jc,kc,lxzb)=b/c
               enddo
            enddo
         enddo

         !
         !  Communitcate piecewise linear interpolation coefficients
         !
         !  NB: This assumes a particular ordering of interplation
         !      indexing. We should make this indirect, and hence
         !      independent of the parameter ordering.
         !
         DO I=lxyl, lxzb
            call halo_exchange(KGC, CI(1,1,1,I), halof)
         ENDDO

         !
         !  (2D) Piecewise bilinear interpolation in coarse-grid planes
         !

         !
         !  Interpolation of fine-grid points on coarse (x,y)-planes
         !  that lie at the intersection of fine-only x-lines and
         !  fine-only y-lines.
         !
         k=KSTART
         do kc=KCSTART,KCEND
            k=k+2
            j=JSTARTO
            do jc=JCSTARTO,JCENDO
               j=j+2
               i=ISTARTO
               do ic=ICSTARTO,ICENDO
                  i=i+2
                  dnw=SO(i-1,j,k,kpnw)+SO(i-1,j,k,kbnw)&
     &                 +SO(i-1,j,k+1,kbse)
                  dn=SO(i-1,j,k,kps)+SO(i-1,j,k,kbn)+SO(i-1,j,k+1,kbs)
                  dne=SO(i,j,k,kpsw)+SO(i,j,k,kbne)+SO(i,j,k+1,kbsw)
                  dw=SO(i-1,j-1,k,kpw)+SO(i-1,j-1,k,kbw)&
     &                 +SO(i-1,j-1,k+1,kbe)
                  de=SO(i,j-1,k,kpw)+SO(i,j-1,k,kbe)+SO(i,j-1,k+1,kbw)
                  dsw=SO(i-1,j-1,k,kpsw)+SO(i-1,j-1,k,kbsw)&
     &                 +SO(i-1,j-1,k+1,kbne)
                  ds=SO(i-1,j-1,k,kps)+SO(i-1,j-1,k,kbs)&
     &                 +SO(i-1,j-1,k+1,kbn)
                  dse=SO(i,j-1,k,kpnw)+SO(i,j-1,k,kbse)&
     &                 +SO(i,j-1,k+1,kbnw)
                  ep=MIN(abs((dsw+dw+dnw)/so(i-1,j-1,k,kp)),&
     &                   abs((dnw+dn+dne)/so(i-1,j-1,k,kp)),&
     &                   abs((dne+de+dse)/so(i-1,j-1,k,kp)),&
     &                   abs((dse+ds+dsw)/so(i-1,j-1,k,kp))&
     &                   )
                  dp=dw+dnw+dn+dne+de+dse+ds+dsw
                  sum=SO(i-1,j-1,k,kp)-SO(i-1,j-1,k,kb)&
     &                 -SO(i-1,j-1,k+1,kb)
                  dp=dp+(sum-dp)*MAX(sum-(rONE+ep)*dp,rZERO)&
     &                 /(abs(sum-(rONE+ep)*dp)+eMACH)
                  dp=rONE/dp
                  CI(ic,jc,kc,lxynw)=dp*(dnw+CI(ic-1,jc,kc,lxya)*dw&
     &                 +CI(ic,jc,kc,lxyl)*dn)
                  CI(ic,jc,kc,lxyne)=dp*(dne+CI(ic,jc,kc,lxyr)*dn&
     &                 +CI(ic,jc,kc,lxya)*de)
                  CI(ic,jc,kc,lxyse)=dp*(dse+CI(ic,jc,kc,lxyb)*de&
     &                 +CI(ic,jc-1,kc,lxyr)*ds)
                  CI(ic,jc,kc,lxysw)=dp*(dsw+CI(ic,jc-1,kc,lxyl)*ds&
     &                 +CI(ic-1,jc,kc,lxyb)*dw)
               enddo
            enddo
         enddo

         !
         !  Interpolation of fine-grid points on coarse (x,z)-planes
         !  that lie at the intersection of fine-only x-lines and
         !  fine-only z-lines.
         !
         k=KSTARTO
         do kc=KCSTARTO,KCENDO
            k=k+2
            j=JSTART
            do jc=JCSTART,JCEND
               j=j+2
               i=ISTARTO
               do ic=ICSTARTO,ICENDO
                  i=i+2
                  dnw=SO(i-1,j+1,k,kbse)+SO(i-1,j,k,kbe)&
     &                 +SO(i-1,j,k,kbne)
                  dn=SO(i-1,j+1,k,kbs)+SO(i-1,j,k,kb)+SO(i-1,j,k,kbn)
                  dne=SO(i,j+1,k,kbsw)+SO(i,j,k,kbw)+SO(i,j,k,kbnw)
                  dw=SO(i-1,j+1,k-1,kpnw)+SO(i-1,j,k-1,kpw)&
     &                 +SO(i-1,j,k-1,kpsw)
                  de=SO(i,j+1,k-1,kpsw)+SO(i,j,k-1,kpw)&
     &                 +SO(i,j,k-1,kpnw)
                  dsw=SO(i-1,j+1,k-1,kbnw)+SO(i-1,j,k-1,kbw)&
     &                 +SO(i-1,j,k-1,kbsw)
                  ds=SO(i-1,j+1,k-1,kbn)+SO(i-1,j,k-1,kb)&
     &                 +SO(i-1,j,k-1,kbs)
                  dse=SO(i,j+1,k-1,kbne)+SO(i,j,k-1,kbe)&
     &                 +SO(i,j,k-1,kbse)
                  ep=MIN(abs((dsw+dw+dnw)/so(i-1,j,k-1,kp)),&
     &                 abs((dnw+dn+dne)/so(i-1,j,k-1,kp)),&
     &                 abs((dne+de+dse)/so(i-1,j,k-1,kp)),&
     &                 abs((dse+ds+dsw)/so(i-1,j,k-1,kp)))
                  dp=dw+dnw+dn+dne+de+dse+ds+dsw
                  sum=SO(i-1,j,k-1,kp)-SO(i-1,j+1,k-1,kps)&
     &                 -SO(i-1,j,k-1,kps)
                  dp=dp+(sum-dp)*MAX(sum-(rONE+ep)*dp,rZERO)&
     &                 /(abs(sum-(rONe+ep)*dp)+eMACH)
                  dp=rONE/dp
                  CI(ic,jc,kc,lxznw)=dp*(dnw+CI(ic-1,jc,kc,lxza)*dw&
     &                 +CI(ic,jc,kc,lxyl)*dn)
                  CI(ic,jc,kc,lxzne)=dp*(dne+CI(ic,jc,kc,lxyr)*dn&
     &                 +CI(ic,jc,kc,lxza)*de)
                  CI(ic,jc,kc,lxzse)=dp*(dse+CI(ic,jc,kc,lxzb)*de&
     &                 +CI(ic,jc,kc-1,lxyr)*ds)
                  CI(ic,jc,kc,lxzsw)=dp*(dsw+CI(ic,jc,kc-1,lxyl)*ds&
     &                 +CI(ic-1,jc,kc,lxzb)*dw)
               enddo
            enddo
         enddo

         !
         !  Interpolation of fine-grid points on coarse (y,z)-planes
         !  that lie at the intersection of fine-only y-lines and
         !  fine-only z-lines.
         !
         k=KSTARTO
         do kc=KCSTARTO,KCENDO
            k=k+2
            j=JSTARTO
            do jc=JCSTARTO,JCENDO
               j=j+2
               i=ISTART
               do ic=ICSTART,ICEND
                  i=i+2
                  dnw=SO(i,j,k,kbse)+SO(i,j,k,kbs)+SO(i+1,j,k,kbsw)
                  dn=SO(i,j-1,k,kbe)+SO(i,j-1,k,kb)+SO(i+1,j-1,k,kbw)
                  dne=SO(i,j-1,k,kbne)+SO(i,j-1,k,kbn)&
     &                 +SO(i+1,j-1,k,kbnw)
                  dw=SO(i,j,k-1,kpnw)+SO(i,j,k-1,kps)&
     &                 +SO(i+1,j,k-1,kpsw)
                  de=SO(i,j-1,k-1,kpsw)+SO(i,j-1,k-1,kps)&
     &                 +SO(i+1,j-1,k-1,kpnw)
                  dsw=SO(i,j,k-1,kbnw)+SO(i,j,k-1,kbn)&
     &                 +SO(i+1,j,k-1,kbne)
                  ds=SO(i,j-1,k-1,kbw)+SO(i,j-1,k-1,kb)&
     &                 +SO(i+1,j-1,k-1,kbe)
                  dse=SO(i,j-1,k-1,kbsw)+SO(i,j-1,k-1,kbs)&
     &                 +SO(i+1,j-1,k-1,kbse)
                  ep=MIN(abs((dsw+dw+dnw)/so(i,j-1,k-1,kp)),&
     &                   abs((dnw+dn+dne)/so(i,j-1,k-1,kp)),&
     &                   abs((dne+de+dse)/so(i,j-1,k-1,kp)),&
     &                   abs((dse+ds+dsw)/so(i,j-1,k-1,kp))&
     &                   )
                  dp=dw+dnw+dn+dne+de+dse+ds+dsw
                  sum=SO(i,j-1,k-1,kp)-SO(i,j-1,k-1,kpw)&
     &                 -SO(i+1,j-1,k-1,kpw)
                  dp=dp+(sum-dp)*MAX(sum-(rONE+ep)*dp,rZERO)&
     &                 /(abs(sum-(rONE+ep)*dp)+eMACH)
                  dp=rONE/dp
                  CI(ic,jc,kc,lyznw)=dp*(dnw+CI(ic,jc,kc,lxza)*dw&
     &                 +CI(ic,jc,kc,lxya)*dn)
                  CI(ic,jc,kc,lyzne)=dp*(dne+CI(ic,jc,kc,lxyb)*dn&
     &                 +CI(ic,jc-1,kc,lxza)*de)
                  CI(ic,jc,kc,lyzse)=dp*(dse+CI(ic,jc-1,kc,lxzb)*de&
     &                 +CI(ic,jc,kc-1,lxyb)*ds)
                  CI(ic,jc,kc,lyzsw)=dp*(dsw+CI(ic,jc,kc-1,lxya)*ds&
     &                 +CI(ic,jc,kc,lxzb)*dw)
               enddo
            enddo
         enddo

         !
         !  Communitcate piecewise bilinear interpolation coefficients
         !
         !  NB: This assumes a particular ordering of interplation
         !      indexing. We should make this indirect, and hence
         !      independent of the parameter ordering.
         !
         DO I=lxyne, lyzse
            call halo_exchange(KGC, CI(1,1,1,I), halof)
         ENDDO

         !
         !  (3D) Piecewise trilinear interpolation.
         !

         !
         !  Interpolation for fine grid points that are at the
         !  center of a coarse-grid cell (i.e., points that are
         !  on fine-only lines in x, y, and z).
         !
         k=KSTARTO
         DO kc=KCSTARTO,KCENDO
            k=k+2
            j=JSTARTO
            DO jc=JCSTARTO,JCENDO
               j=j+2
               i=ISTARTO
               DO ic=ICSTARTO,ICENDO
                  i=i+2
                  yo(ic,JJF,2,kp)=SO(i-1,j-1,k-1,kpw)+SO(i-1,j,k-1,kpnw)&
     &                 +SO(i-1,j,k-1,kps)+SO(i,j,k-1,kpsw)&
     &                 +SO(i,j-1,k-1,kpw)&
     &                 +SO(i,j-1,k-1,kpnw)+SO(i-1,j-1,k-1,kps)&
     &                 +SO(i-1,j-1,k-1,kpsw)+SO(i-1,j-1,k-1,kb)&
     &                 +SO(i-1,j-1,k-1,kbw)+SO(i-1,j,k-1,kbnw)&
     &                 +SO(i-1,j,k-1,kbn)+SO(i,j,k-1,kbne)&
     &                 +SO(i,j-1,k-1,kbe)&
     &                 +SO(i,j-1,k-1,kbse)+SO(i-1,j-1,k-1,kbs)&
     &                 +SO(i-1,j-1,k-1,kbsw)+SO(i-1,j-1,k,kb)&
     &                 +SO(i-1,j-1,k,kbe)+SO(i-1,j,k,kbse)&
     &                 +SO(i-1,j,k,kbs)&
     &                 +SO(i,j,k,kbsw)+SO(i,j-1,k,kbw)+SO(i,j-1,k,kbnw)&
     &                 +SO(i-1,j-1,k,kbn)+SO(i-1,j-1,k,kbne)
                  yo(ic,JJF,2,kpw)&
     &                 =MIN(abs(SO(i-1,j-1,k-1,kpw)+SO(i-1,j,k-1,kpnw)&
     &                 +SO(i-1,j,k,kbse)+SO(i-1,j-1,k,kbe)&
     &                 +SO(i-1,j-1,k,kbne)&
     &                 +SO(i-1,j-1,k-1,kpsw)+SO(i-1,j-1,k-1,kbsw)&
     &                 +SO(i-1,j-1,k-1,kbw)+SO(i-1,j,k-1,kbnw))/&
     &                    so(i-1,j-1,k-1,kp),&
     &                 abs(SO(i,j-1,k-1,kpw)+SO(i,j,k-1,kpsw)&
     &                 +SO(i,j,k,kbsw)+SO(i,j-1,k,kbw)&
     &                 +SO(i,j-1,k,kbnw)+SO(i,j-1,k-1,kpnw)&
     &                 +SO(i,j-1,k-1,kbse)&
     &                 +SO(i,j-1,k-1,kbe)+SO(i,j,k-1,kbne))/&
     &                    so(i-1,j-1,k-1,kp),&
     &                 abs(SO(i-1,j,k-1,kps)+SO(i-1,j,k-1,kpnw)&
     &                 +SO(i-1,j,k,kbse)+SO(i-1,j,k,kbs)+SO(i,j,k,kbsw)&
     &                 +SO(i,j,k-1,kpsw)+SO(i,j,k-1,kbne)&
     &                 +SO(i-1,j,k-1,kbn)&
     &                 +SO(i-1,j,k-1,kbnw)),abs(SO(i-1,j-1,k-1,kps)&
     &                 +SO(i-1,j-1,k-1,kpsw)+SO(i-1,j-1,k,kbne)&
     &                 +SO(i-1,j-1,k,kbn)+SO(i,j-1,k,kbnw)&
     &                 +SO(i,j-1,k-1,kpnw)+SO(i,j-1,k-1,kbse)&
     &                 +SO(i-1,j-1,k-1,kbs)+SO(i,j-1,k-1,kbse))/&
     &                    so(i-1,j-1,k-1,kp))
                  yo(ic,JJF,2,kpw)=MIN(yo(ic,JJF,2,kpw),&
     &                 abs(SO(i-1,j-1,k-1,kb)+SO(i-1,j-1,k-1,kbw)&
     &                 +SO(i-1,j,k-1,kbnw)+SO(i-1,j,k-1,kbn)&
     &                 +SO(i,j,k-1,kbne)&
     &                 +SO(i,j-1,k-1,kbe)+SO(i,j-1,k-1,kbse)&
     &                 +SO(i-1,j-1,k-1,kbs)&
     &                 +SO(i-1,j-1,k-1,kbsw)),abs(SO(i-1,j-1,k,kb)&
     &                 +SO(i-1,j-1,k,kbe)+SO(i-1,j,k,kbse)&
     &                 +SO(i-1,j,k,kbs)&
     &                 +SO(i,j,k,kbsw)+SO(i,j-1,k,kbw)+SO(i,j-1,k,kbnw)&
     &                 +SO(i-1,j-1,k,kbn)+SO(i-1,j-1,k,kbne))/&
     &                    so(i-1,j-1,k-1,kp))
                  yo(ic,JJF,2,kp)&
     &                 =yo(ic,JJF,2,kp)+(SO(i-1,j-1,k-1,kp)&
     &                 -yo(ic,JJF,2,kp))*MAX(SO(i-1,j-1,k-1,kp)&
     &                 -(rONE+yo(ic,JJF,2,kpw))*yo(ic,JJF,2,kp),rZERO)&
     &                 /(abs(SO(i-1,j-1,k-1,kp)-(rONE+yo(ic,JJF,2,kpw))&
     &                 *yo(ic,JJF,2,kp))+eMACH)
                  yo(ic,JJF,2,kp)=rONE/yo(ic,JJF,2,kp)
                  CI(ic,jc,kc,ltnw)&
     &                 =yo(ic,JJF,2,kp)*(SO(i-1,j,k,kbse)&
     &                 +CI(ic-1,jc,kc,lyznw)&
     &                 *SO(i-1,j-1,k-1,kpw)+CI(ic-1,jc,kc,lxza)&
     &                 *SO(i-1,j,k-1,kpnw)&
     &                 +CI(ic,jc,kc,lxznw)*SO(i-1,j,k-1,kps)&
     &                 +CI(ic-1,jc,kc,lxya)&
     &                 *SO(i-1,j-1,k,kbe)+CI(ic,jc,kc,lxyl)&
     &                 *SO(i-1,j,k,kbs)&
     &                 +CI(ic,jc,kc,lxynw)*SO(i-1,j-1,k,kb))
                  CI(ic,jc,kc,ltne)&
     &                 =yo(ic,JJF,2,kp)*(SO(i,j,k,kbsw)&
     &                 +CI(ic,jc,kc,lxzne)&
     &                 *SO(i-1,j,k-1,kps)+CI(ic,jc,kc,lxza)&
     &                 *SO(i,j,k-1,kpsw)&
     &                 +CI(ic,jc,kc,lyznw)*SO(i,j-1,k-1,kpw)&
     &                 +CI(ic,jc,kc,lxyr)&
     &                 *SO(i-1,j,k,kbs)+CI(ic,jc,kc,lxya)&
     &                 *SO(i,j-1,k,kbw)&
     &                 +CI(ic,jc,kc,lxyne)*SO(i-1,j-1,k,kb))
                  CI(ic,jc,kc,lbnw)&
     &                 =yo(ic,JJF,2,kp)*(SO(i-1,j,k-1,kbnw)&
     &                 +CI(ic-1,jc,kc-1,lxya)*SO(i-1,j-1,k-1,kbw)&
     &                 +CI(ic,jc,kc-1,lxyl)*SO(i-1,j,k-1,kbn)&
     &                 +CI(ic,jc,kc-1,lxynw)*SO(i-1,j-1,k-1,kb)&
     &                 +CI(ic-1,jc,kc,lyzsw)*SO(i-1,j-1,k-1,kpw)&
     &                 +CI(ic-1,jc,kc,lxzb)*SO(i-1,j,k-1,kpnw)&
     &                 +CI(ic,jc,kc,lxzsw)*SO(i-1,j,k-1,kps))
                  CI(ic,jc,kc,lbne)&
     &                 =yo(ic,JJF,2,kp)*(SO(i,j,k-1,kbne)&
     &                 +CI(ic,jc,kc-1,lxyne)&
     &                 *SO(i-1,j-1,k-1,kb)+CI(ic,jc,kc-1,lxyr)&
     &                 *SO(i-1,j,k-1,kbn)&
     &                 +CI(ic,jc,kc-1,lxya)*SO(i,j-1,k-1,kbe)&
     &                 +CI(ic,jc,kc,lxzse)&
     &                 *SO(i-1,j,k-1,kps)+CI(ic,jc,kc,lxzb)&
     &                 *SO(i,j,k-1,kpsw)&
     &                 +CI(ic,jc,kc,lyzsw)*SO(i,j-1,k-1,kpw))
                  CI(ic,jc,kc,lbsw)&
     &                 =yo(ic,JJF,2,kp)*(SO(i-1,j-1,k-1,kbsw)&
     &                 +CI(ic-1,jc,kc-1,lxyb)*SO(i-1,j-1,k-1,kbw)&
     &                 +CI(ic,jc,kc-1,lxysw)&
     &                 *SO(i-1,j-1,k-1,kb)+CI(ic,jc-1,kc-1,lxyl)&
     &                 *SO(i-1,j-1,k-1,&
     &                 kbs)+CI(ic-1,jc,kc,lyzse)*SO(i-1,j-1,k-1,kpw)&
     &                 +CI(ic,jc-1,kc,lxzsw)*SO(i-1,j-1,k-1,kps)&
     &                 +CI(ic-1,jc-1,kc,lxzb)*SO(i-1,j-1,k-1,kpsw))
                  CI(ic,jc,kc,ltsw)&
     &                 =yo(ic,JJF,2,kp)*(SO(i-1,j-1,k,kbne)&
     &                 +CI(ic-1,jc,kc,lxyb)*SO(i-1,j-1,k,kbe)&
     &                 +CI(ic,jc,kc,lxysw)*SO(i-1,j-1,k,kb)&
     &                 +CI(ic,jc-1,kc,lxyl)*SO(i-1,j-1,k,kbn)&
     &                 +CI(ic-1,jc,kc,lyzne)&
     &                 *SO(i-1,j-1,k-1,kpw)+CI(ic,jc-1,kc,lxznw)&
     &                 *SO(i-1,j-1,k-1,kps)&
     &                 +CI(ic-1,jc-1,kc,lxza)*SO(i-1,j-1,k-1,kpsw))
                  CI(ic,jc,kc,ltse)&
     &                 =yo(ic,JJF,2,kp)*(SO(i,j-1,k,kbnw)&
     &                 +CI(ic,jc-1,kc,lxyr)&
     &                 *SO(i-1,j-1,k,kbn)+CI(ic,jc,kc,lxyse)&
     &                 *SO(i-1,j-1,k,kb)&
     &                 +CI(ic,jc,kc,lxyb)*SO(i,j-1,k,kbw)&
     &                 +CI(ic,jc-1,kc,lxzne)&
     &                 *SO(i-1,j-1,k-1,kps)+CI(ic,jc,kc,lyzne)&
     &                 *SO(i,j-1,k-1,kpw)&
     &                 +CI(ic,jc-1,kc,lxza)*SO(i,j-1,k-1,kpnw))
                  CI(ic,jc,kc,lbse)&
     &                 =yo(ic,JJF,2,kp)*(SO(i,j-1,k-1,kbse)&
     &                 +CI(ic,jc-1,kc-1,lxyr)*SO(i-1,j-1,k-1,kbs)&
     &                 +CI(ic,jc,kc-1,lxyse)*SO(i-1,j-1,k-1,kb)&
     &                 +CI(ic,jc,kc-1,lxyb)*SO(i,j-1,k-1,kbe)&
     &                 +CI(ic,jc-1,kc,lxzse)*SO(i-1,j-1,k-1,kps)&
     &                 +CI(ic,jc,kc,lyzse)*SO(i,j-1,k-1,kpw)&
     &                 +CI(ic,jc-1,kc,lxzb)*SO(i,j-1,k-1,kpnw))
               ENDDO
            ENDDO
         ENDDO

         !
         !  Communitcate piecewise trilinear interpolation coefficients
         !
         !  NB: This assumes a particular ordering of interplation
         !      indexing. We should make this indirect, and hence
         !      independent of the parameter ordering.
         !
         DO I=lbsw, ltse
            call halo_exchange(KGC, CI(1,1,1,I), halof)
         ENDDO

      ELSE ! if kgf.ge.NOG.and.ifd.eq.1

         !
         !  Interpolation based on a 7-point stencil
         !

         !
         !  (1D) Piecewise linear interpolation along coarse grid lines
         !

         !
         !  Interpolate fine-grid points lieing on coarse x-lines
         !
         !  NB: storage viewed as points lieing on fine-only
         !      y-lines of a coarse xy-plane (i.e., coarse k-plane)
         !
         k=KSTART
         DO kc=KCSTART,KCEND
            k=k+2
            j=JSTART
            DO jc=JCSTART,JCEND
               j=j+2
               i=ISTARTO
               DO ic=ICSTARTO,ICENDO
                  i=i+2
                  a=SO(i-1,j,k,kpw)
                  b=SO(i,j,k,kpw)
                  ep=MIN(abs(a/so(i-1,j,k,kp)),abs(b)/so(i-1,j,k,kp))
                  c=a+b+SO(i-1,j,k,kps)+SO(i-1,j+1,k,kps)&
     &                 +SO(i-1,j,k,kb)+SO(i-1,j,k+1,kb)
                  c=a+b+(SO(i-1,j,k,kp)-c)*MAX(SO(i-1,j,k,kp)&
     &                 -(rONE+ep)*c,rZERO)/(abs(SO(i-1,j,k,kp)&
     &                 -(rONE+ep)*c)+eMACH)
                  CI(ic,jc,kc,lxyl)=a/c
                  CI(ic,jc,kc,lxyr)=b/c
               ENDDO
            ENDDO
         ENDDO

         !
         !  Interpolate fine-grid points lieing on coarse y-lines
         !
         !  NB: storage viewed as points lieing on fine-only
         !      x-lines of a coarse xy-plane (i.e., coarse k-plane)
         !
         k=KSTART
         DO kc=KCSTART,KCEND
            k=k+2
            j=JSTARTO
            DO jc=JCSTARTO,JCENDO
               j=j+2
               i=ISTART
               DO ic=ICSTART,ICEND
                  i=i+2
                  a=SO(i,j,k,kps)
                  b=SO(i,j-1,k,kps)
                  c=a+b+SO(i,j-1,k,kpw)+SO(i+1,j-1,k,kpw)&
     &                 +SO(i,j-1,k,kb)+SO(i,j-1,k+1,kb)
                  ep=MIN(abs(a/so(i,j-1,k,kp)),abs(b/so(i,j-1,k,kp)))
                  c=a+b+(SO(i,j-1,k,kp)-c)*MAX(SO(i,j-1,k,kp)&
     &                 -(rONE+ep)*c,rZERO)/(abs(SO(i,j-1,k,kp)&
     &                 -(rONE+ep)*c)+eMACH)
                  CI(ic,jc,kc,lxya)=a/c
                  CI(ic,jc,kc,lxyb)=b/c
               ENDDO
            ENDDO
         ENDDO

         !
         !  Interpolate fine-grid points lieing on coarse z-lines
         !
         !  NB: storage viewed as points lieing on fine-only
         !      x-lines of a coarse xz-plane (i.e., coarse j-plane)
         !
         k=KSTARTO
         DO kc=KCSTARTO,KCENDO
            k=k+2
            j=JSTART
            DO jc=JCSTART,JCEND
               j=j+2
               i=ISTART
               DO ic=ICSTART,ICEND
                  i=i+2
                  a=SO(i,j,k,kb)
                  b=SO(i,j,k-1,kb)
                  c=a+b+SO(i,j,k-1,kpw)+SO(i+1,j,k-1,kpw)&
     &                 +SO(i,j+1,k-1,kps)+SO(i,j,k-1,kps)
                  ep=MIN(abs(a/so(i,j,k-1,kp)),abs(b/so(i,j,k-1,kp)))
                  c=a+b+(SO(i,j,k-1,kp)-c)*MAX(SO(i,j,k-1,kp)&
     &                 -(rONE+ep)*c,rZERO)/(abs(SO(i,j,k-1,kp)&
     &                 -(rONE+ep)*c)+eMACH)
                  CI(ic,jc,kc,lxza)=a/c
                  CI(ic,jc,kc,lxzb)=b/c
               ENDDO
            ENDDO
         ENDDO

         !
         !  Communitcate piecewise linear interpolation coefficients
         !
         !  NB: This assumes a particular ordering of interplation
         !      indexing. We should make this indirect, and hence
         !      independent of the parameter ordering.
         !
         DO I=lxyl, lxzb
            call halo_exchange(KGC, CI(1,1,1,I), halof)
         ENDDO

         !
         !  (2D) Piecewise bilinear interpolation in coarse-grid planes
         !

         !
         !  Interpolation of fine-grid points on coarse (x,y)-planes
         !  that lie at the intersection of fine-only x-lines and
         !  fine-only y-lines.
         !
         i=ISTARTO
         DO ic=ICSTARTO,ICENDO
            i=i+2
            j=JSTARTO
            DO jc=JCSTARTO,JCENDO
               j=j+2
               k=KSTART
               DO kc=KCSTART,KCEND
                  k=k+2
                  dn=SO(i-1,j,k,kps)
                  dw=SO(i-1,j-1,k,kpw)
                  de=SO(i,j-1,k,kpw)
                  ds=SO(i-1,j-1,k,kps)
                  dp=dw+dn+de+ds
                  sum=SO(i-1,j-1,k,kp)-SO(i-1,j-1,k,kb)&
     &                 -SO(i-1,j-1,k+1,kb)
                  ep=MIN(abs(dw/so(i-1,j-1,k,kp)),&
     &                   abs(dn/so(i-1,j-1,k,kp)),&
     &                   abs(de/so(i-1,j-1,k,kp)),&
     &                   abs(ds/so(i-1,j-1,k,kp)) &
     &                   )
                  dp=dp+(sum-dp)*MAX(sum-(rONE+ep)*dp,rZERO)&
     &                 /(abs(sum-(rONE+ep)*dp)+eMACH)
                  dp=rONE/dp
                  CI(ic,jc,kc,lxynw)=dp*(CI(ic-1,jc,kc,lxya)*dw&
     &                 +CI(ic,jc,kc,lxyl)*dn)
                  CI(ic,jc,kc,lxyne)=dp*(CI(ic,jc,kc,lxyr)*dn&
     &                 +CI(ic,jc,kc,lxya)*de)
                  CI(ic,jc,kc,lxyse)=dp*(CI(ic,jc,kc,lxyb)*de&
     &                 +CI(ic,jc-1,kc,lxyr)*ds)
                  CI(ic,jc,kc,lxysw)=dp*(CI(ic,jc-1,kc,lxyl)*ds&
     &                 +CI(ic-1,jc,kc,lxyb)*dw)
               ENDDO
            ENDDO
         ENDDO

         !
         !  Interpolation of fine-grid points on coarse (x,z)-planes
         !  that lie at the intersection of fine-only x-lines and
         !  fine-only z-lines.
         !
         k=KSTARTO
         DO kc=KCSTARTO,KCENDO
            k=k+2
            j=JSTART
            DO jc=JCSTART,JCEND
               j=j+2
               i=ISTARTO
               DO ic=ICSTARTO,ICENDO
                  i=i+2
                  dn=SO(i-1,j,k,kb)
                  dw=SO(i-1,j,k-1,kpw)
                  de=SO(i,j,k-1,kpw)
                  ds=SO(i-1,j,k-1,kb)
                  dp=dw+dn+de+ds
                  sum=SO(i-1,j,k-1,kp)-SO(i-1,j+1,k-1,kps)&
     &                 -SO(i-1,j,k-1,kps)
                  ep=MIN(abs(dw/so(i-1,j,k-1,kp)),&
     &                   abs(dn/so(i-1,j,k-1,kp)),&
     &                   abs(de/so(i-1,j,k-1,kp)),&
     &                   abs(ds/so(i-1,j,k-1,kp)) &
     &                   )
                  dp=dp+(sum-dp)*MAX(sum-(rONE+ep)*dp,rZERO)&
     &                 /(abs(sum-(rONE+ep)*dp)+eMACH)
                  dp=rONE/dp
                  CI(ic,jc,kc,lxznw)=dp*(CI(ic-1,jc,kc,lxza)*dw&
     &                 +CI(ic,jc,kc,lxyl)*dn)
                  CI(ic,jc,kc,lxzne)=dp*(CI(ic,jc,kc,lxyr)*dn&
     &                 +CI(ic,jc,kc,lxza)*de)
                  CI(ic,jc,kc,lxzse)=dp*(CI(ic,jc,kc,lxzb)*de&
     &                 +CI(ic,jc,kc-1,lxyr)*ds)
                  CI(ic,jc,kc,lxzsw)=dp*(CI(ic,jc,kc-1,lxyl)*ds&
     &                 +CI(ic-1,jc,kc,lxzb)*dw)
               ENDDO
            ENDDO
         ENDDO

         !
         !  Interpolation of fine-grid points on coarse (y,z)-planes
         !  that lie at the intersection of fine-only y-lines and
         !  fine-only z-lines.
         !
         k=KSTARTO
         DO kc=KCSTARTO,KCENDO
            k=k+2
            j=JSTARTO
            DO jc=JCSTARTO,JCENDO
               j=j+2
               i=ISTART
               DO ic=ICSTART,ICEND
                  i=i+2
                  dn=SO(i,j-1,k,kb)
                  dw=SO(i,j,k-1,kps)
                  de=SO(i,j-1,k-1,kps)
                  ds=SO(i,j-1,k-1,kb)
                  dp=dw+dn+de+ds
                  sum=SO(i,j-1,k-1,kp)-SO(i,j-1,k-1,kpw)&
     &                 -SO(i+1,j-1,k-1,kpw)
                  ep=MIN(abs(dw/so(i,j-1,k-1,kp)),&
     &                   abs(dn/so(i,j-1,k-1,kp)),&
     &                   abs(de/so(i,j-1,k-1,kp)),&
     &                   abs(ds/so(i,j-1,k-1,kp)) &
     &                   )
                  dp=dp+(sum-dp)*MAX(sum-(rONE+ep)*dp,rZERO)&
     &                 /(abs(sum-(rONE+ep)*dp)+eMACH)
                  dp=rONE/dp
                  CI(ic,jc,kc,lyznw)=dp*(CI(ic,jc,kc,lxza)*dw&
     &                 +CI(ic,jc,kc,lxya)*dn)
                  CI(ic,jc,kc,lyzne)=dp*(CI(ic,jc,kc,lxyb)*dn&
     &                 +CI(ic,jc-1,kc,lxza)*de)
                  CI(ic,jc,kc,lyzse)=dp*(CI(ic,jc-1,kc,lxzb)*de&
     &                 +CI(ic,jc,kc-1,lxyb)*ds)
                  CI(ic,jc,kc,lyzsw)=dp*(CI(ic,jc,kc-1,lxya)*ds&
     &                 +CI(ic,jc,kc,lxzb)*dw)
               ENDDO
            ENDDO
         ENDDO

         !
         !  Communitcate piecewise bilinear interpolation coefficients
         !
         !  NB: This assumes a particular ordering of interplation
         !      indexing. We should make this indirect, and hence
         !      independent of the parameter ordering.
         !
         DO I=lxyne, lyzse
            call halo_exchange(KGC, CI(1,1,1,I), halof)
         ENDDO

         !
         !  (3D) Piecewise trilinear interpolation.
         !

         !
         !  Interpolation for fine grid points that are at the
         !  center of a coarse-grid cell (i.e., points that are
         !  on fine-only lines in x, y, and z).
         !
         k=KSTARTO
         DO kc=KCSTARTO,KCENDO
            k=k+2
            j=JSTARTO
            DO jc=JCSTARTO,JCENDO
               j=j+2
               i=ISTARTO
               DO ic=ICSTARTO,ICENDO
                  i=i+2
                  dp=SO(i-1,j-1,k-1,kpw)+SO(i-1,j,k-1,kps)&
     &                 +SO(i,j-1,k-1,kpw)+SO(i-1,j-1,k-1,kps)&
     &                 +SO(i-1,j-1,k-1,kb)+SO(i-1,j-1,k,kb)
                  ep=MIN(abs(so(i-1,j-1,k-1,kpw)/so(i-1,j-1,k-1,kp)),&
     &                   abs(so(i-1,j,k-1,kps)/so(i-1,j-1,k-1,kp)),&
     &                   abs(so(i,j-1,k-1,kpw)/so(i-1,j-1,k-1,kp)),&
     &                   abs(so(i-1,j-1,k-1,kps)/so(i-1,j-1,k-1,kp)),&
     &                   abs(so(i-1,j-1,k-1,kb)/so(i-1,j-1,k-1,kp)),&
     &                   abs(so(i-1,j-1,k,kb)/so(i-1,j-1,k-1,kp)) &
     &                   )
                  dp=(SO(i-1,j-1,k-1,kp)-dp)*MAX(SO(i-1,j-1,k-1,kp)&
     &                 -(rONE+ep)*dp,rZERO)/(abs(SO(i-1,j-1,k-1,kp)&
     &                 -(rONE+ep)*dp)+eMACH)+dp
                  dp=rONE/dp
                  CI(ic,jc,kc,ltnw)=dp*(CI(ic-1,jc,kc,lyznw)&
     &                 *SO(i-1,j-1,k-1,kpw)&
     &                 +CI(ic,jc,kc,lxznw)*SO(i-1,j,k-1,kps)&
     &                 +CI(ic,jc,kc,lxynw)*SO(i-1,j-1,k,kb))
                  CI(ic,jc,kc,ltne)=dp*(CI(ic,jc,kc,lxzne)&
     &                 *SO(i-1,j,k-1,kps)&
     &                 +CI(ic,jc,kc,lyznw)*SO(i,j-1,k-1,kpw)&
     &                 +CI(ic,jc,kc,lxyne)*SO(i-1,j-1,k,kb))
                  CI(ic,jc,kc,lbnw)=dp*(CI(ic,jc,kc-1,lxynw)&
     &                 *SO(i-1,j-1,k-1,kb)&
     &                 +CI(ic-1,jc,kc,lyzsw)*SO(i-1,j-1,k-1,kpw)&
     &                 +CI(ic,jc,kc,lxzsw)*SO(i-1,j,k-1,kps))
                  CI(ic,jc,kc,lbne)=dp*(CI(ic,jc,kc-1,lxyne)&
     &                 *SO(i-1,j-1,k-1,kb)&
     &                 +CI(ic,jc,kc,lxzse)*SO(i-1,j,k-1,kps)&
     &                 +CI(ic,jc,kc,lyzsw)*SO(i,j-1,k-1,kpw))
                  CI(ic,jc,kc,lbsw)=dp*(CI(ic,jc,kc-1,lxysw)&
     &                 *SO(i-1,j-1,k-1,kb)&
     &                 +CI(ic-1,jc,kc,lyzse)*SO(i-1,j-1,k-1,kpw)&
     &                 +CI(ic,jc-1,kc,lxzsw)*SO(i-1,j-1,k-1,kps))
                  CI(ic,jc,kc,ltsw)=dp*(CI(ic,jc,kc,lxysw)&
     &                 *SO(i-1,j-1,k,kb)&
     &                 +CI(ic-1,jc,kc,lyzne)*SO(i-1,j-1,k-1,kpw)&
     &                 +CI(ic,jc-1,kc,lxznw)*SO(i-1,j-1,k-1,kps))
                  CI(ic,jc,kc,ltse)=dp*(CI(ic,jc,kc,lxyse)&
     &                 *SO(i-1,j-1,k,kb)&
     &                 +CI(ic,jc-1,kc,lxzne)*SO(i-1,j-1,k-1,kps)&
     &                 +CI(ic,jc,kc,lyzne)*SO(i,j-1,k-1,kpw))
                  CI(ic,jc,kc,lbse)=dp*(CI(ic,jc,kc-1,lxyse)&
     &                 *SO(i-1,j-1,k-1,kb)&
     &                 +CI(ic,jc-1,kc,lxzse)*SO(i-1,j-1,k-1,kps)&
     &                 +CI(ic,jc,kc,lyzse)*SO(i,j-1,k-1,kpw))
               ENDDO
            ENDDO
         ENDDO

         !
         !  Communitcate piecewise trilinear interpolation coefficients
         !
         !  NB: This assumes a particular ordering of interplation
         !      indexing. We should make this indirect, and hence
         !      independent of the parameter ordering.
         !
         DO I=lbsw, ltse
            call halo_exchange(KGC, CI(1,1,1,I), halof)
         ENDDO

      ENDIF ! if kgf.ge.NOG.and.ifd.eq.1

! --------------------------------------------------------------------
! >>>>>>>>>>>>>>>>>> END:  INTERPOLATION CONSTRUCTION <<<<<<<<<<<<<<
! --------------------------------------------------------------------

! ======================================================================

      RETURN
      END

max_stars_repo_name: cedar-framework/cedar

max_issues_repo_name: cedar-framework/cedar

max_forks_repo_name: cedar-framework/cedar


text:
!*****************************************************************************
!
! Different implementations of a simple single precision A-X plus Y function
! on a CPU and GPU, to verify this code works.
!
! x = ax + y
! x,y,z: vector
! a: scalar
!*****************************************************************************

!*****************************************************************************
! Setup/driver code
!*****************************************************************************
subroutine test_saxpy_cuda_kernels
  use cudafor
  use saxpy_cuda_kernels_mod

  real, device :: x_d(2**20), y_d(2**20)
  x_d = 1.0
  y_d = 2.0

  ! Perform SAXPY on 1M elements
  call saxpy_gpu_cuda_kernel<<<4096, 256>>>(2.0, x_d, y_d)

  return
end subroutine test_saxpy_cuda_kernels

max_stars_repo_name: bblakeley/FGPU

max_issues_repo_name: bblakeley/FGPU

max_forks_repo_name: bblakeley/FGPU


! ===================================================================!
!   UBMOD_S -- One dimensional mass balance model for unsaturated    !
!              water flow and solute transport. Version 1.00.        !
!                                                                    !
!   Designed by Wei Mao, Yan Zhu and Jinzhong Yang.                  !
!                                                                    !
!   Cited: Mao W, Yang J, Zhu Y, et al. An efficient soil water      !
!          balance model based on hybrid numerical and statistical   !
!          methods[J]. Journal of hydrology, 2018, 559: 721-735,     !
!          doi: 10.1016/j.jhydrol.2018.02.074.                       !
!                                                                    !
!   Feel free to contact us if you have any question.                !
!       Email: weimao@whu.edu.cn, zyan0701@163.com                   !
!                                                                    !
!                            Last modified: Sep, 2019, by Wei Mao.   !
! ===================================================================!
! ====================================================================
!     Input files:
!     1. Essential files.
!         SELECTOR.IN     The basic input information.
!         uz.IN           The discrete information.
!     2. Optional files.
!         cropdat.dat     The simple crop model.
!         01.wea          Meteorological data.
! ====================================================================
! ====================================================================
!     Output files:
!
! ====================================================================
!   storage Routing Method~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!   *************************end of documentation.
    PROGRAM Main_Program
    USE parm
    IMPLICIT NONE
    REAL (KIND=KR) :: t1, t2, Zero
    INTEGER (KIND=4) :: lengthpath, ierr
    CHARACTER (100) :: filename, datapath

    Zero = 0.0_KR
    Terr = 0_KI

    filename = 'LEVEL_01.DIR'
    OPEN(10,file=filename, status='old',err=901)
    READ(10,'(a)',err=904) datapath
    CLOSE(10)
    lengthpath = Len_Trim(datapath)
!---Open the interface files. The relpath is used here.
! ====================================================================
!   input files
    OPEN(33,file=datapath(1:lengthpath)//'/Selector.in', status='old',err=902)
    OPEN(32,file=datapath(1:lengthpath)//'/Profile.in' , status='old',err=902)
!   output files
    OPEN(80,file=datapath(1:lengthpath)//'/Obsnode.out' ,status='unknown',err=903) ! Observation.
    OPEN(81,file=datapath(1:lengthpath)//'/TPrint.out'  ,status='unknown',err=903) ! Profile results.
    OPEN(89,file=datapath(1:lengthpath)//'/A_Level.out' ,status='unknown',err=903) ! Accumulated results.
    OPEN(90,file=datapath(1:lengthpath)//'/Runtime.out' ,status='unknown',err=903) ! Run time.
    OPEN(91,file=datapath(1:lengthpath)//'/Balance.out' ,status='unknown',err=903) ! Mass Balance.
    OPEN(99,file=datapath(1:lengthpath)//'/Error.out'   ,status='unknown',err=903) ! Error message.
! ====================================================================

!-----Begin of the program.
! ====================================================================
!     subroutine about input information.
!     call for basic information.
    CALL Selector_In
    IF (Terr.ne.0) GOTO (905, 906, 907, 908, 920, 921, 922) Terr
!     call for node information in 1D.
    CALL Profile_In
    IF (Terr.ne.0) GOTO 909
! ====================================================================

!-----preparation of the calculation
! ====================================================================
!   CPU time.
    CALL CPU_time (t1)
!   The initial water amount in model.
    CALL Balance_Initial
    IF (Terr.ne.0) GOTO 916
!     Diffusion model.
    IF (lWat) CALL Diffusion_Model
!     Call for reference Evaportranspiration and division of E&T.
    CALL Boundary_Cond(datapath)
    IF (Terr.ne.0) GOTO (911,912,913,918,910,905) Terr
! ====================================================================

!-----Begin time loop.
100 CONTINUE

! ====================================================================
!   The steady water flow.
    IF (.NOT.lWat) THEN
        CALL Steady_Flow
    ELSE

! ====================================================================
!     Four main processes.
!     Firstly, divide infiltration and surface runoff.
!     There is no surface hydrological process by now.
          !if () then
          !    call SCS
          !elseif () then
          !    call Green-Ampt
          !else
          !    qair = Q_infiltration
          !endif
        !   Set upper boundary condition.
        CALL Set_Input

! ====================================================================
!     Secondly, advective movement driven by gravitational potential.
!       "Tipping-bucket" method.
        CALL Water_Redis
        IF (Terr.ne.0) GOTO (930, 933) Terr

! ====================================================================
!     Thirdly, source/sink term.
!     open the Files that stored E&T and the rain, and the writen ETa.
        IF(bup >= Zero) CALL Water_SetET
        IF (Terr.ne.0) GOTO (931) Terr

! ====================================================================
!     Last, Diffusive soil water movement driven by matric potential.
        CALL Water_Diff
        IF (Terr.ne.0) GOTO (932) Terr
    ENDIF

    IF (lChem .and. Cind > 1) THEN
        th = tht
        t = t - dt
        ct = t + dt
        Cind = Cind
        dt = dt/Cind
        index = Cind
        GOTO 100
    ENDIF

!---Solute transport module.
! ====================================================================
    IF (lChem) THEN
!   1st, the solute convection
        CALL Solute_Conv

!   2nd, the solute source/sink term
        CALL Solute_Sour

!   3nd, the exchange between mobile and immobile regions.
        IF (MIM) THEN
            CALL Solute_Exch
        ELSE
            Conc3 = Conc2
        ENDIF

!   4th, the diffusive term driven by concentration difference.
        CALL Solute_DiffD
    ENDIF

! ====================================================================
!     Output control.
!     Output the hydraulic head and soil moisture in 1D model.
    CALL Obs_Out
    IF (Terr.ne.0) GOTO (915) Terr
!   Call for the water balance in 1D and 3D model.
    CALL BalanceT
    IF (Terr.ne.0) GOTO 916
!   call for new time and time step.
    WRITE(*,*)"t=",sngl(t)

!   P-Level information
    IF (abs(TPrint(Plevel)-t) < Tol) THEN
        CALL thOut
        IF (Terr.ne.0) GOTO (917) Terr
        Plevel = Plevel + 1
    ENDIF

    IF(abs(t-Tend) <= Tol) THEN
        GOTO 200
    ENDIF

    IF (dt < dtOld) THEN
        dt = dtOld
    ENDIF
    dtOld = dt
    CALL Tcontrol
    TLevel = TLevel + 1
    t = t+dt

    !IF (bup == 1) THEN
    !    WRITE(150,'(A18,F10.3,2F10.6)')date,t,Epa,Tra
    !    CALL DateAdd(date,1,date)
    !    Epa=Zero
    !    Tra=Zero
    !ENDIF

    GOTO 100

200 CALL CPU_time (t2)
    WRITE(90,*,err=914)'Real time [sec]',t2-t1

    CLOSE(80)
    CLOSE(81)
    CLOSE(89)
    CLOSE(90)
    CLOSE(91)
    CLOSE(99)
    CLOSE(110)
    CLOSE(130)
    CLOSE(150)

    STOP

! Input or output error.
901 ierr=1
    GOTO 999
902 ierr=2
    GOTO 999
903 ierr=3
    GOTO 999
904 ierr=4
    GOTO 999
905 ierr=5
    GOTO 999
906 ierr=6
    GOTO 999
907 ierr=7
    GOTO 999
908 ierr=8
    GOTO 999
909 ierr=9
    GOTO 999
910 ierr=10
    GOTO 999
911 ierr=11
    GOTO 999
912 ierr=12
    GOTO 999
913 ierr=13
    GOTO 999
914 ierr=14
    GOTO 999
915 ierr=15
    GOTO 999
916 ierr=16
    GOTO 999
917 ierr=17
    GOTO 999
918 ierr=18
    GOTO 999
919 ierr=19
    GOTO 999
! Data check error
920 ierr=20
    GOTO 999
921 ierr=21
    GOTO 999
922 ierr=22
    GOTO 999
923 ierr=23
    GOTO 999

930 ierr=30
    GOTO 999
931 ierr=31
    GOTO 999
932 ierr=32
    GOTO 999
933 ierr=33
    GOTO 999

999 CALL Error_Out(ierr)
    PAUSE
    STOP

    END PROGRAM Main_Program

    SUBROUTINE Error_Out(ierr)
    IMPLICIT NONE
    INTEGER (KIND=4) :: ierr
    CHARACTER (LEN=100), DIMENSION(40) :: cErr

    cErr( 1)='Open file error in file LEVEL_01.DIR !'
    cErr( 2)='Open file error for input files !'
    cErr( 3)='Error when opening the output file !'
    cErr( 4)='Error when reading from an input file LEVEL_01.DIR !'
    cErr( 5)='Error when reading from an input file Selector.in Basic Information !'
    cErr( 6)='Error when reading from an input file Selector.in Material Information !'
    cErr( 7)='Error when reading from an input file Selector.in Time Information !'
    cErr( 8)='Error when reading from an input file Selector.in Solute Information !'
    cErr( 9)='Error when reading from an input file Profile.in !'
    cErr(10)='Error when reading from an input file Met.in !'
    cErr(11)='Error when opening or reading from an input file 01.wea !'
    cErr(12)='Error when opening or reading from an input file cropdat.dat !'
    cErr(13)='Error when opening or reading from an input file crp/et0/eti !'
    cErr(14)='Error when writing to the output file Runtime.out !'
    cErr(15)='Error when writing to the output file A_Level.out !'
    cErr(16)='Error when writing to the output file Balance.out !'
    cErr(17)='Error when writing to the output file TPrint.out !'
    cErr(18)='Error when writing to the output file eta.dat !'
    !
    cErr(20)='Input option is wrong !'
    cErr(21)='Dimension of the Array is exceeded !'
    cErr(22)='The soil water hydraulic parameters is wrong !'

    cErr(30)='Mass balance error in Water_Redis module !'
    cErr(31)='Mass balance error in Water_SetET module !'
    cErr(32)='Mass balance error in Water_Diff module !'
    cErr(32)='Too much infiltration water in Water_Redis module !'

    WRITE(*,*) cErr(ierr)
    WRITE(99,*) cErr(ierr)

    END SUBROUTINE Error_Out


text:
c
c
c     ###################################################
c     ##  COPYRIGHT (C)  1993  by  Jay William Ponder  ##
c     ##              All Rights Reserved              ##
c     ###################################################
c
c     ###############################################################
c     ##                                                           ##
c     ##  subroutine eangang2  --  angle-angle Hessian; numerical  ##
c     ##                                                           ##
c     ###############################################################
c
c
c     "eangang2" calculates the angle-angle potential energy
c     second derivatives with respect to Cartesian coordinates
c     using finite difference methods
c
c
      subroutine eangang2 (i)
      implicit none
      include 'sizes.i'
      include 'angang.i'
      include 'angle.i'
      include 'atoms.i'
      include 'deriv.i'
      include 'group.i'
      include 'hessn.i'
      integer i,j,k,iangang
      integer ia,ib,ic,id,ie
      real*8 eps,fgrp
      real*8 old,term
      real*8, allocatable :: d0(:,:)
      logical proceed
      logical twosided
c
c
c     set stepsize for derivatives and default group weight
c
      eps = 1.0d-5
      fgrp = 1.0d0
      twosided = .false.
      if (n .le. 50)  twosided = .true.
c
c     perform dynamic allocation of some local arrays
c
      allocate (d0(3,n))
c
c     compute numerical angle-angle Hessian for current atom
c
      do iangang = 1, nangang
         j = iaa(1,iangang)
         k = iaa(2,iangang)
         ia = iang(1,j)
         ib = iang(2,j)
         ic = iang(3,j)
         id = iang(1,k)
         ie = iang(3,k)
c
c     decide whether to compute the current interaction
c
         proceed = (i.eq.ia .or. i.eq.ib .or. i.eq.ic
     &                 .or. i.eq.id .or. i.eq.ie)
         if (proceed .and. use_group)
     &      call groups (proceed,fgrp,ia,ib,ic,id,ie,0)
c
c     eliminate any duplicate atoms in the pair of angles
c
         if (proceed) then
            if (id.eq.ia .or. id.eq.ic)  then
               id = ie
               ie = 0
            else if (ie.eq.ia .or. ie.eq.ic) then
               ie = 0
            end if
            term = fgrp / eps
c
c     find first derivatives for the base structure
c
            if (.not. twosided) then
               call eangang2a (iangang)
               do j = 1, 3
                  d0(j,ia) = deaa(j,ia)
                  d0(j,ib) = deaa(j,ib)
                  d0(j,ic) = deaa(j,ic)
                  d0(j,id) = deaa(j,id)
                  if (ie .ne. 0)  d0(j,ie) = deaa(j,ie)
               end do
            end if
c
c     find numerical x-components via perturbed structures
c
            old = x(i)
            if (twosided) then
               x(i) = x(i) - 0.5d0*eps
               call eangang2a (iangang)
               do j = 1, 3
                  d0(j,ia) = deaa(j,ia)
                  d0(j,ib) = deaa(j,ib)
                  d0(j,ic) = deaa(j,ic)
                  d0(j,id) = deaa(j,id)
                  if (ie .ne. 0)  d0(j,ie) = deaa(j,ie)
               end do
            end if
            x(i) = x(i) + eps
            call eangang2a (iangang)
            x(i) = old
            do j = 1, 3
               hessx(j,ia) = hessx(j,ia) + term*(deaa(j,ia)-d0(j,ia))
               hessx(j,ib) = hessx(j,ib) + term*(deaa(j,ib)-d0(j,ib))
               hessx(j,ic) = hessx(j,ic) + term*(deaa(j,ic)-d0(j,ic))
               hessx(j,id) = hessx(j,id) + term*(deaa(j,id)-d0(j,id))
               if (ie .ne. 0)
     &            hessx(j,ie) = hessx(j,ie) + term*(deaa(j,ie)-d0(j,ie))
            end do
c
c     find numerical y-components via perturbed structures
c
            old = y(i)
            if (twosided) then
               y(i) = y(i) - 0.5d0*eps
               call eangang2a (iangang)
               do j = 1, 3
                  d0(j,ia) = deaa(j,ia)
                  d0(j,ib) = deaa(j,ib)
                  d0(j,ic) = deaa(j,ic)
                  d0(j,id) = deaa(j,id)
                  if (ie .ne. 0)  d0(j,ie) = deaa(j,ie)
               end do
            end if
            y(i) = y(i) + eps
            call eangang2a (iangang)
            y(i) = old
            do j = 1, 3
               hessy(j,ia) = hessy(j,ia) + term*(deaa(j,ia)-d0(j,ia))
               hessy(j,ib) = hessy(j,ib) + term*(deaa(j,ib)-d0(j,ib))
               hessy(j,ic) = hessy(j,ic) + term*(deaa(j,ic)-d0(j,ic))
               hessy(j,id) = hessy(j,id) + term*(deaa(j,id)-d0(j,id))
               if (ie .ne. 0)
     &            hessy(j,ie) = hessy(j,ie) + term*(deaa(j,ie)-d0(j,ie))
            end do
c
c     find numerical z-components via perturbed structures
c
            old = z(i)
            if (twosided) then
               z(i) = z(i) - 0.5d0*eps
               call eangang2a (iangang)
               do j = 1, 3
                  d0(j,ia) = deaa(j,ia)
                  d0(j,ib) = deaa(j,ib)
                  d0(j,ic) = deaa(j,ic)
                  d0(j,id) = deaa(j,id)
                  if (ie .ne. 0)  d0(j,ie) = deaa(j,ie)
               end do
            end if
            z(i) = z(i) + eps
            call eangang2a (iangang)
            z(i) = old
            do j = 1, 3
               hessz(j,ia) = hessz(j,ia) + term*(deaa(j,ia)-d0(j,ia))
               hessz(j,ib) = hessz(j,ib) + term*(deaa(j,ib)-d0(j,ib))
               hessz(j,ic) = hessz(j,ic) + term*(deaa(j,ic)-d0(j,ic))
               hessz(j,id) = hessz(j,id) + term*(deaa(j,id)-d0(j,id))
               if (ie .ne. 0)
     &            hessz(j,ie) = hessz(j,ie) + term*(deaa(j,ie)-d0(j,ie))
            end do
         end if
      end do
c
c     perform deallocation of some local arrays
c
      deallocate (d0)
      return
      end
c
c
c     ################################################################
c     ##                                                            ##
c     ##  subroutine eangang2a  --  angle-angle interaction derivs  ##
c     ##                                                            ##
c     ################################################################
c
c
c     "eangang2a" calculates the angle-angle first derivatives for
c     a single interaction with respect to Cartesian coordinates;
c     used in computation of finite difference second derivatives
c
c
      subroutine eangang2a (i)
      implicit none
      include 'sizes.i'
      include 'angang.i'
      include 'angle.i'
      include 'angpot.i'
      include 'atoms.i'
      include 'bound.i'
      include 'deriv.i'
      include 'math.i'
      integer i,j,k
      integer ia,ib,ic,id,ie
      real*8 angle,term
      real*8 dot,cosine
      real*8 dt1,deddt1
      real*8 dt2,deddt2
      real*8 xia,yia,zia
      real*8 xib,yib,zib
      real*8 xic,yic,zic
      real*8 xid,yid,zid
      real*8 xie,yie,zie
      real*8 xab,yab,zab
      real*8 xcb,ycb,zcb
      real*8 xdb,ydb,zdb
      real*8 xeb,yeb,zeb
      real*8 rab2,rcb2
      real*8 rdb2,reb2
      real*8 xp,yp,zp,rp
      real*8 xq,yq,zq,rq
      real*8 terma,termc
      real*8 termd,terme
      real*8 dedxia,dedyia,dedzia
      real*8 dedxib,dedyib,dedzib
      real*8 dedxic,dedyic,dedzic
      real*8 dedxid,dedyid,dedzid
      real*8 dedxie,dedyie,dedzie
c
c
c     set the coordinates of the involved atoms
c
      j = iaa(1,i)
      k = iaa(2,i)
      ia = iang(1,j)
      ib = iang(2,j)
      ic = iang(3,j)
      id = iang(1,k)
      ie = iang(3,k)
      xia = x(ia)
      yia = y(ia)
      zia = z(ia)
      xib = x(ib)
      yib = y(ib)
      zib = z(ib)
      xic = x(ic)
      yic = y(ic)
      zic = z(ic)
      xid = x(id)
      yid = y(id)
      zid = z(id)
      xie = x(ie)
      yie = y(ie)
      zie = z(ie)
c
c     zero out the first derivative components
c
      deaa(1,ia) = 0.0d0
      deaa(2,ia) = 0.0d0
      deaa(3,ia) = 0.0d0
      deaa(1,ib) = 0.0d0
      deaa(2,ib) = 0.0d0
      deaa(3,ib) = 0.0d0
      deaa(1,ic) = 0.0d0
      deaa(2,ic) = 0.0d0
      deaa(3,ic) = 0.0d0
      deaa(1,id) = 0.0d0
      deaa(2,id) = 0.0d0
      deaa(3,id) = 0.0d0
      deaa(1,ie) = 0.0d0
      deaa(2,ie) = 0.0d0
      deaa(3,ie) = 0.0d0
c
c     compute the values of the two bond angles
c
      xab = xia - xib
      yab = yia - yib
      zab = zia - zib
      xcb = xic - xib
      ycb = yic - yib
      zcb = zic - zib
      xdb = xid - xib
      ydb = yid - yib
      zdb = zid - zib
      xeb = xie - xib
      yeb = yie - yib
      zeb = zie - zib
      if (use_polymer) then
         call image (xab,yab,zab)
         call image (xcb,ycb,zcb)
         call image (xdb,ydb,zdb)
         call image (xeb,yeb,zeb)
      end if
      rab2 = xab*xab + yab*yab + zab*zab
      rcb2 = xcb*xcb + ycb*ycb + zcb*zcb
      rdb2 = xdb*xdb + ydb*ydb + zdb*zdb
      reb2 = xeb*xeb + yeb*yeb + zeb*zeb
      xp = ycb*zab - zcb*yab
      yp = zcb*xab - xcb*zab
      zp = xcb*yab - ycb*xab
      xq = yeb*zdb - zeb*ydb
      yq = zeb*xdb - xeb*zdb
      zq = xeb*ydb - yeb*xdb
      rp = sqrt(xp*xp + yp*yp + zp*zp)
      rq = sqrt(xq*xq + yq*yq + zq*zq)
      if (rp*rq .ne. 0.0d0) then
         dot = xab*xcb + yab*ycb + zab*zcb
         cosine = dot / sqrt(rab2*rcb2)
         cosine = min(1.0d0,max(-1.0d0,cosine))
         angle = radian * acos(cosine)
         dt1 = angle - anat(j)
         dot = xdb*xeb + ydb*yeb + zdb*zeb
         cosine = dot / sqrt(rdb2*reb2)
         cosine = min(1.0d0,max(-1.0d0,cosine))
         angle = radian * acos(cosine)
         dt2 = angle - anat(k)
c
c     get the energy and master chain rule terms for derivatives
c
         term = radian * aaunit * kaa(i)
         deddt1 = term * dt2
         deddt2 = term * dt1
c
c     find chain rule terms for the first bond angle deviation
c
         terma = -deddt1 / (rab2*rp)
         termc = deddt1 / (rcb2*rp)
         dedxia = terma * (yab*zp-zab*yp)
         dedyia = terma * (zab*xp-xab*zp)
         dedzia = terma * (xab*yp-yab*xp)
         dedxic = termc * (ycb*zp-zcb*yp)
         dedyic = termc * (zcb*xp-xcb*zp)
         dedzic = termc * (xcb*yp-ycb*xp)
c
c     find chain rule terms for the second bond angle deviation
c
         termd = -deddt2 / (rdb2*rq)
         terme = deddt2 / (reb2*rq)
         dedxid = termd * (ydb*zq-zdb*yq)
         dedyid = termd * (zdb*xq-xdb*zq)
         dedzid = termd * (xdb*yq-ydb*xq)
         dedxie = terme * (yeb*zq-zeb*yq)
         dedyie = terme * (zeb*xq-xeb*zq)
         dedzie = terme * (xeb*yq-yeb*xq)
c
c     get the central atom derivative terms by difference
c
         dedxib = -dedxia - dedxic - dedxid - dedxie
         dedyib = -dedyia - dedyic - dedyid - dedyie
         dedzib = -dedzia - dedzic - dedzid - dedzie
c
c     set the angle-angle interaction first derivatives
c
         deaa(1,ia) = deaa(1,ia) + dedxia
         deaa(2,ia) = deaa(2,ia) + dedyia
         deaa(3,ia) = deaa(3,ia) + dedzia
         deaa(1,ib) = deaa(1,ib) + dedxib
         deaa(2,ib) = deaa(2,ib) + dedyib
         deaa(3,ib) = deaa(3,ib) + dedzib
         deaa(1,ic) = deaa(1,ic) + dedxic
         deaa(2,ic) = deaa(2,ic) + dedyic
         deaa(3,ic) = deaa(3,ic) + dedzic
         deaa(1,id) = deaa(1,id) + dedxid
         deaa(2,id) = deaa(2,id) + dedyid
         deaa(3,id) = deaa(3,id) + dedzid
         deaa(1,ie) = deaa(1,ie) + dedxie
         deaa(2,ie) = deaa(2,ie) + dedyie
         deaa(3,ie) = deaa(3,ie) + dedzie
      end if
      return
      end

ext: f

max_stars_repo_name: htna/HCsbLib

max_issues_repo_name: htna/HCsbLib

max_forks_repo_name: htna/HCsbLib


text:
C A program to test the for2py-pp pre-processing code.  This program has
C several continuation lines.
C
C NOTE: The output of for2py-pp.py for this input file contains lines that
C exceed fixed-form Fortran's line-length limits.  To compile the Fortran
C ode output by for2py-pp.py, use the command
C
C       gfortran -ffree-form -ffree-line-length-none <file-name>

        PROGRAM MAIN

        IMPLICIT NONE

        REAL E,Fc,Lai, nb,N,PT,Pg, di,PAR
        REAL rm,dwf,int, TMAX,TMIN, p1, sla
        REAL PD,EMP1,EMP2,Lfmax,dwc, TMN
        REAL dwr,dw,dn,w,wc,wr,wf,tb,intot, dLAI, FL
        INTEGER DOY, endsim, COUNT
        CHARACTER*10 DYN

        REAL SWFAC1, SWFAC2

        OPEN (2,FILE='INFILE.INP',STATUS='UNKNOWN')

        READ(2,10) Lfmax, EMP2,EMP1,PD,nb,rm,fc,tb,intot,n,lai,w,wr,wc
     &     ,p1,sla
   10   FORMAT(17(1X,F7.4))
        CLOSE(2)

        WRITE(*,11)
        WRITE(*,12)
   11   FORMAT('Results of plant growth simulation: ')
   12   FORMAT(/
     &/,'                Accum',
     &/,'       Number    Temp                                    Leaf',
     &/,'  Day      of  during   Plant  Canopy    Root   Fruit    Area',
     &/,'   of    Leaf  Reprod  Weight  Weight  Weight  weight   Index',
     &/,' Year   Nodes    (oC)  (g/m2)  (g/m2)  (g/m2)  (g/m2) (m2/m2)',
     &/,' ----  ------  ------  ------  ------  ------  ------  ------')

        WRITE(*,11)
        WRITE(*,12)

        WRITE(*,10) Lfmax, EMP2,EMP1,PD,nb,rm,fc,tb,intot,n,lai,w,wr,wc
     &     ,p1,sla

        STOP
        END PROGRAM MAIN

ext: for


text:
If you are looking for the fruit trees fruit, you can find them on our local orange trees, or at a supermarket.  Also check out gold, which can often be orangish.

  Orange Court is home to several local businesses.
  You could live at Orange Tree Apartments in Davis.
  Orange Hut is a Chinese Restaurants Chinese restaurant.  Of course, Orange Chicken is a popular dish, and once choice of a standard to compare them, like in the battle between Panda Express vs. Davis Noodle City.
  People drink Orangina in Davis.
  Bakers Square used to have an orange pie.

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


program l
    use utilities
    use percolation
    implicit none

    integer :: num_Ls, num_xs, i, j, num_samples, fileunit
    integer, dimension(:), allocatable :: Ls
    real(kind=dp), dimension(:), allocatable :: xs
    character(len=:), allocatable :: filename
    real(kind=dp) :: tolerance, p_x

    num_samples = 16000
    tolerance = 0.0001d0

    xs = [0.3, 0.8]
    num_xs = size(xs)
    Ls = [25, 50, 100, 200, 400, 800]
    num_Ls = size(Ls)

    do i = 1, num_xs
        associate(x => xs(i))
            filename = "tmp/l_" // stringfromint(i) // ".dat"
            open(newunit=fileunit, file=filename, status="replace")
            do j = 1, num_Ls
                associate(L => Ls(j))
                    p_x = spanning_probability_inverse(x, L, num_samples/L, tolerance)
                    write(unit=fileunit, fmt=*) L, p_x
                end associate
            end do
            close(unit=fileunit)
        end associate
    end do

end program l


text:
      MODULE GWFWELMODULE
        INTEGER,SAVE,POINTER  ::NWELLS,MXWELL,NWELVL,IWELCB,IPRWEL
        INTEGER,SAVE,POINTER  ::NPWEL,IWELPB,NNPWEL
        CHARACTER(LEN=16),SAVE, DIMENSION(:),   POINTER     ::WELAUX
        REAL,             SAVE, DIMENSION(:,:), POINTER     ::WELL
      TYPE GWFWELTYPE
        INTEGER,POINTER  ::NWELLS,MXWELL,NWELVL,IWELCB,IPRWEL
        INTEGER,POINTER  ::NPWEL,IWELPB,NNPWEL
        CHARACTER(LEN=16), DIMENSION(:),   POINTER     ::WELAUX
        REAL,              DIMENSION(:,:), POINTER     ::WELL
      END TYPE
      TYPE(GWFWELTYPE), SAVE:: GWFWELDAT(10)
      END MODULE GWFWELMODULE

      SUBROUTINE GWF2WEL7AR(IN,IGRID)
C     ******************************************************************
C     ALLOCATE ARRAY STORAGE FOR WELL PACKAGE
C     ******************************************************************
C
C        SPECIFICATIONS:
C     ------------------------------------------------------------------
      USE GLOBAL,       ONLY:IOUT,NCOL,NROW,NLAY,IFREFM
      USE GWFWELMODULE, ONLY:NWELLS,MXWELL,NWELVL,IWELCB,IPRWEL,NPWEL,
     1                       IWELPB,NNPWEL,WELAUX,WELL
C
      CHARACTER*200 LINE
C     ------------------------------------------------------------------
      ALLOCATE(NWELLS,MXWELL,NWELVL,IWELCB,IPRWEL)
      ALLOCATE(NPWEL,IWELPB,NNPWEL)
C
C1------IDENTIFY PACKAGE AND INITIALIZE NWELLS.
      WRITE(IOUT,1)IN
    1 FORMAT(1X,/1X,'WEL -- WELL PACKAGE, VERSION 7, 5/2/2005',
     1' INPUT READ FROM UNIT ',I4)
      NWELLS=0
      NNPWEL=0
C
C2------READ MAXIMUM NUMBER OF WELLS AND UNIT OR FLAG FOR
C2------CELL-BY-CELL FLOW TERMS.
      CALL URDCOM(IN,IOUT,LINE)
      CALL UPARLSTAL(IN,IOUT,LINE,NPWEL,MXPW)
      IF(IFREFM.EQ.0) THEN
         READ(LINE,'(2I10)') MXACTW,IWELCB
         LLOC=21
      ELSE
         LLOC=1
         CALL URWORD(LINE,LLOC,ISTART,ISTOP,2,MXACTW,R,IOUT,IN)
         CALL URWORD(LINE,LLOC,ISTART,ISTOP,2,IWELCB,R,IOUT,IN)
      END IF
      WRITE(IOUT,3) MXACTW
    3 FORMAT(1X,'MAXIMUM OF ',I6,' ACTIVE WELLS AT ONE TIME')
      IF(IWELCB.LT.0) WRITE(IOUT,7)
    7 FORMAT(1X,'CELL-BY-CELL FLOWS WILL BE PRINTED WHEN ICBCFL NOT 0')
      IF(IWELCB.GT.0) WRITE(IOUT,8) IWELCB
    8 FORMAT(1X,'CELL-BY-CELL FLOWS WILL BE SAVED ON UNIT ',I4)
C
C3------READ AUXILIARY VARIABLES AND PRINT FLAG.
      ALLOCATE(WELAUX(20))
      NAUX=0
      IPRWEL=1
   10 CALL URWORD(LINE,LLOC,ISTART,ISTOP,1,N,R,IOUT,IN)
      IF(LINE(ISTART:ISTOP).EQ.'AUXILIARY' .OR.
     1        LINE(ISTART:ISTOP).EQ.'AUX') THEN
         CALL URWORD(LINE,LLOC,ISTART,ISTOP,1,N,R,IOUT,IN)
         IF(NAUX.LT.20) THEN
            NAUX=NAUX+1
            WELAUX(NAUX)=LINE(ISTART:ISTOP)
            WRITE(IOUT,12) WELAUX(NAUX)
   12       FORMAT(1X,'AUXILIARY WELL VARIABLE: ',A)
         END IF
         GO TO 10
      ELSE IF(LINE(ISTART:ISTOP).EQ.'NOPRINT') THEN
         WRITE(IOUT,13)
   13    FORMAT(1X,'LISTS OF WELL CELLS WILL NOT BE PRINTED')
         IPRWEL = 0
         GO TO 10
      END IF
C3A-----THERE ARE FOUR INPUT VALUES PLUS ONE LOCATION FOR
C3A-----CELL-BY-CELL FLOW.
      NWELVL=5+NAUX
C
C4------ALLOCATE SPACE FOR THE WELL DATA.
      IWELPB=MXACTW+1
      MXWELL=MXACTW+MXPW
      ALLOCATE (WELL(NWELVL,MXWELL))
C
C5------READ NAMED PARAMETERS.
      WRITE(IOUT,18) NPWEL
   18 FORMAT(1X,//1X,I5,' Well parameters')
      IF(NPWEL.GT.0) THEN
        LSTSUM=IWELPB
        DO 120 K=1,NPWEL
          LSTBEG=LSTSUM
          CALL UPARLSTRP(LSTSUM,MXWELL,IN,IOUT,IP,'WEL','Q',1,
     &                   NUMINST)
          NLST=LSTSUM-LSTBEG
          IF(NUMINST.EQ.0) THEN
C5A-----READ PARAMETER WITHOUT INSTANCES.
            CALL ULSTRD(NLST,WELL,LSTBEG,NWELVL,MXWELL,1,IN,
     &        IOUT,'WELL NO.  LAYER   ROW   COL   STRESS FACTOR',
     &        WELAUX,20,NAUX,IFREFM,NCOL,NROW,NLAY,4,4,IPRWEL)
          ELSE
C5B-----READ INSTANCES.
            NINLST=NLST/NUMINST
            DO 110 I=1,NUMINST
            CALL UINSRP(I,IN,IOUT,IP,IPRWEL)
            CALL ULSTRD(NINLST,WELL,LSTBEG,NWELVL,MXWELL,1,IN,
     &        IOUT,'WELL NO.  LAYER   ROW   COL   STRESS FACTOR',
     &        WELAUX,20,NAUX,IFREFM,NCOL,NROW,NLAY,4,4,IPRWEL)
            LSTBEG=LSTBEG+NINLST
  110       CONTINUE
          END IF
  120   CONTINUE
      END IF
C
C6------RETURN
      CALL SGWF2WEL7PSV(IGRID)
      RETURN
      END
      SUBROUTINE GWF2WEL7RP(IN,IGRID)
C     ******************************************************************
C     READ WELL DATA FOR A STRESS PERIOD
C     ******************************************************************
C
C        SPECIFICATIONS:
C     ------------------------------------------------------------------
      USE GLOBAL,       ONLY:IOUT,NCOL,NROW,NLAY,IFREFM
      USE GWFWELMODULE, ONLY:NWELLS,MXWELL,NWELVL,IPRWEL,NPWEL,
     1                       IWELPB,NNPWEL,WELAUX,WELL
C
      CHARACTER*6 CWELL
C     ------------------------------------------------------------------
      CALL SGWF2WEL7PNT(IGRID)
C
C1----READ NUMBER OF WELLS (OR FLAG SAYING REUSE WELL DATA).
C1----AND NUMBER OF PARAMETERS
      IF(NPWEL.GT.0) THEN
        IF(IFREFM.EQ.0) THEN
           READ(IN,'(2I10)') ITMP,NP
        ELSE
           READ(IN,*) ITMP,NP
        END IF
      ELSE
         NP=0
         IF(IFREFM.EQ.0) THEN
            READ(IN,'(I10)') ITMP
         ELSE
            READ(IN,*) ITMP
         END IF
      END IF
C
C------Calculate some constants.
      NAUX=NWELVL-5
      IOUTU = IOUT
      IF (IPRWEL.EQ.0) IOUTU=-IOUTU
C
C1A-----IF ITMP LESS THAN ZERO REUSE NON-PARAMETER DATA. PRINT MESSAGE.
C1A-----IF ITMP=>0, SET NUMBER OF NON-PARAMETER WELLS EQUAL TO ITMP.
      IF(ITMP.LT.0) THEN
         WRITE(IOUT,6)
    6    FORMAT(1X,/
     1    1X,'REUSING NON-PARAMETER WELLS FROM LAST STRESS PERIOD')
      ELSE
         NNPWEL=ITMP
      END IF
C
C1B-----IF THERE ARE NEW NON-PARAMETER WELLS, READ THEM.
      MXACTW=IWELPB-1
      IF(ITMP.GT.0) THEN
         IF(NNPWEL.GT.MXACTW) THEN
            WRITE(IOUT,99) NNPWEL,MXACTW
   99       FORMAT(1X,/1X,'THE NUMBER OF ACTIVE WELLS (',I6,
     1                     ') IS GREATER THAN MXACTW(',I6,')')
            CALL USTOP(' ')
         END IF
         CALL ULSTRD(NNPWEL,WELL,1,NWELVL,MXWELL,1,IN,IOUT,
     1            'WELL NO.  LAYER   ROW   COL   STRESS RATE',
     2             WELAUX,20,NAUX,IFREFM,NCOL,NROW,NLAY,4,4,IPRWEL)
      END IF
      NWELLS=NNPWEL
C
C1C-----IF THERE ARE ACTIVE WELL PARAMETERS, READ THEM AND SUBSTITUTE
      CALL PRESET('Q')
      NREAD=NWELVL-1
      IF(NP.GT.0) THEN
         DO 30 N=1,NP
         CALL UPARLSTSUB(IN,'WEL',IOUTU,'Q',WELL,NWELVL,MXWELL,NREAD,
     1                MXACTW,NWELLS,4,4,
     2            'WELL NO.  LAYER   ROW   COL   STRESS RATE',
     3            WELAUX,20,NAUX)
   30    CONTINUE
      END IF
C
C3------PRINT NUMBER OF WELLS IN CURRENT STRESS PERIOD.
      CWELL=' WELLS'
      IF(NWELLS.EQ.1) CWELL=' WELL '
      WRITE(IOUT,101) NWELLS,CWELL
  101 FORMAT(1X,/1X,I6,A)
C
C6------RETURN
      RETURN
      END
      SUBROUTINE GWF2WEL7FM(IGRID)
C     ******************************************************************
C     SUBTRACT Q FROM RHS
C     ******************************************************************
C
C        SPECIFICATIONS:
C     ------------------------------------------------------------------
      USE GLOBAL,       ONLY:IBOUND,RHS,HCOF
      USE GWFWELMODULE, ONLY:NWELLS,WELL
C     ------------------------------------------------------------------
      CALL SGWF2WEL7PNT(IGRID)
C
C1------IF NUMBER OF WELLS <= 0 THEN RETURN.
      IF(NWELLS.LE.0) RETURN
C
C2------PROCESS EACH WELL IN THE WELL LIST.
      DO 100 L=1,NWELLS
      IR=WELL(2,L)
      IC=WELL(3,L)
      IL=WELL(1,L)
      Q=WELL(4,L)
C
C2A-----IF THE CELL IS INACTIVE THEN BYPASS PROCESSING.
      IF(IBOUND(IC,IR,IL).LE.0) GO TO 100
C
C2B-----IF THE CELL IS VARIABLE HEAD THEN SUBTRACT Q FROM
C       THE RHS ACCUMULATOR.
      RHS(IC,IR,IL)=RHS(IC,IR,IL)-Q
  100 CONTINUE
C
C3------RETURN
      RETURN
      END
      SUBROUTINE GWF2WEL7BD(KSTP,KPER,IGRID)
C     ******************************************************************
C     CALCULATE VOLUMETRIC BUDGET FOR WELLS
C     ******************************************************************
C
C        SPECIFICATIONS:
C     ------------------------------------------------------------------
      USE GLOBAL,      ONLY:IOUT,NCOL,NROW,NLAY,IBOUND,BUFF
      USE GWFBASMODULE,ONLY:MSUM,ICBCFL,IAUXSV,DELT,PERTIM,TOTIM,
     1                      VBVL,VBNM
      USE GWFWELMODULE,ONLY:NWELLS,IWELCB,WELL,NWELVL,WELAUX
C
      CHARACTER*16 TEXT
      DOUBLE PRECISION RATIN,RATOUT,QQ
      DATA TEXT /'           WELLS'/
C     ------------------------------------------------------------------
      CALL SGWF2WEL7PNT(IGRID)
C
C1------CLEAR RATIN AND RATOUT ACCUMULATORS, AND SET CELL-BY-CELL
C1------BUDGET FLAG.
      ZERO=0.
      RATIN=ZERO
      RATOUT=ZERO
      IBD=0
      IF(IWELCB.LT.0 .AND. ICBCFL.NE.0) IBD=-1
      IF(IWELCB.GT.0) IBD=ICBCFL
      IBDLBL=0
C
C2-----IF CELL-BY-CELL FLOWS WILL BE SAVED AS A LIST, WRITE HEADER.
      IF(IBD.EQ.2) THEN
         NAUX=NWELVL-5
         IF(IAUXSV.EQ.0) NAUX=0
         CALL UBDSV4(KSTP,KPER,TEXT,NAUX,WELAUX,IWELCB,NCOL,NROW,NLAY,
     1          NWELLS,IOUT,DELT,PERTIM,TOTIM,IBOUND)
      END IF
C
C3------CLEAR THE BUFFER.
      DO 50 IL=1,NLAY
      DO 50 IR=1,NROW
      DO 50 IC=1,NCOL
      BUFF(IC,IR,IL)=ZERO
50    CONTINUE
C
C4------IF THERE ARE NO WELLS, DO NOT ACCUMULATE FLOW.
      IF(NWELLS.EQ.0) GO TO 200
C
C5------LOOP THROUGH EACH WELL CALCULATING FLOW.
      DO 100 L=1,NWELLS
C
C5A-----GET LAYER, ROW & COLUMN OF CELL CONTAINING WELL.
      IR=WELL(2,L)
      IC=WELL(3,L)
      IL=WELL(1,L)
      Q=ZERO
C
C5B-----IF THE CELL IS NO-FLOW OR CONSTANT_HEAD, IGNORE IT.
      IF(IBOUND(IC,IR,IL).LE.0)GO TO 99
C
C5C-----GET FLOW RATE FROM WELL LIST.
      Q=WELL(4,L)
      QQ=Q
C
C5D-----PRINT FLOW RATE IF REQUESTED.
      IF(IBD.LT.0) THEN
         IF(IBDLBL.EQ.0) WRITE(IOUT,61) TEXT,KPER,KSTP
   61    FORMAT(1X,/1X,A,'   PERIOD ',I4,'   STEP ',I3)
         WRITE(IOUT,62) L,IL,IR,IC,Q
   62    FORMAT(1X,'WELL ',I6,'   LAYER ',I3,'   ROW ',I5,'   COL ',I5,
     1       '   RATE ',1PG15.6)
         IBDLBL=1
      END IF
C
C5E-----ADD FLOW RATE TO BUFFER.
      BUFF(IC,IR,IL)=BUFF(IC,IR,IL)+Q
C
C5F-----SEE IF FLOW IS POSITIVE OR NEGATIVE.
      IF(Q.GE.ZERO) THEN
C
C5G-----FLOW RATE IS POSITIVE (RECHARGE). ADD IT TO RATIN.
        RATIN=RATIN+QQ
      ELSE
C
C5H-----FLOW RATE IS NEGATIVE (DISCHARGE). ADD IT TO RATOUT.
        RATOUT=RATOUT-QQ
      END IF
C
C5I-----IF SAVING CELL-BY-CELL FLOWS IN A LIST, WRITE FLOW.  ALSO
C5I-----COPY FLOW TO WELL LIST.
   99 IF(IBD.EQ.2) CALL UBDSVB(IWELCB,NCOL,NROW,IC,IR,IL,Q,
     1                  WELL(:,L),NWELVL,NAUX,5,IBOUND,NLAY)
      WELL(NWELVL,L)=Q
  100 CONTINUE
C
C6------IF CELL-BY-CELL FLOWS WILL BE SAVED AS A 3-D ARRAY,
C6------CALL UBUDSV TO SAVE THEM.
      IF(IBD.EQ.1) CALL UBUDSV(KSTP,KPER,TEXT,IWELCB,BUFF,NCOL,NROW,
     1                          NLAY,IOUT)
C
C7------MOVE RATES, VOLUMES & LABELS INTO ARRAYS FOR PRINTING.
  200 RIN=RATIN
      ROUT=RATOUT
      VBVL(3,MSUM)=RIN
      VBVL(4,MSUM)=ROUT
      VBVL(1,MSUM)=VBVL(1,MSUM)+RIN*DELT
      VBVL(2,MSUM)=VBVL(2,MSUM)+ROUT*DELT
      VBNM(MSUM)=TEXT
C
C8------INCREMENT BUDGET TERM COUNTER(MSUM).
      MSUM=MSUM+1
C
C9------RETURN
      RETURN
      END
      SUBROUTINE GWF2WEL7DA(IGRID)
C  Deallocate WEL MEMORY
      USE GWFWELMODULE
C
        CALL SGWF2WEL7PNT(IGRID)
        DEALLOCATE(NWELLS)
        DEALLOCATE(MXWELL)
        DEALLOCATE(NWELVL)
        DEALLOCATE(IWELCB)
        DEALLOCATE(IPRWEL)
        DEALLOCATE(NPWEL)
        DEALLOCATE(IWELPB)
        DEALLOCATE(NNPWEL)
        DEALLOCATE(WELAUX)
        DEALLOCATE(WELL)
C
      RETURN
      END
      SUBROUTINE SGWF2WEL7PNT(IGRID)
C  Change WEL data to a different grid.
      USE GWFWELMODULE
C
        NWELLS=>GWFWELDAT(IGRID)%NWELLS
        MXWELL=>GWFWELDAT(IGRID)%MXWELL
        NWELVL=>GWFWELDAT(IGRID)%NWELVL
        IWELCB=>GWFWELDAT(IGRID)%IWELCB
        IPRWEL=>GWFWELDAT(IGRID)%IPRWEL
        NPWEL=>GWFWELDAT(IGRID)%NPWEL
        IWELPB=>GWFWELDAT(IGRID)%IWELPB
        NNPWEL=>GWFWELDAT(IGRID)%NNPWEL
        WELAUX=>GWFWELDAT(IGRID)%WELAUX
        WELL=>GWFWELDAT(IGRID)%WELL
C
      RETURN
      END
      SUBROUTINE SGWF2WEL7PSV(IGRID)
C  Save WEL data for a grid.
      USE GWFWELMODULE
C
        GWFWELDAT(IGRID)%NWELLS=>NWELLS
        GWFWELDAT(IGRID)%MXWELL=>MXWELL
        GWFWELDAT(IGRID)%NWELVL=>NWELVL
        GWFWELDAT(IGRID)%IWELCB=>IWELCB
        GWFWELDAT(IGRID)%IPRWEL=>IPRWEL
        GWFWELDAT(IGRID)%NPWEL=>NPWEL
        GWFWELDAT(IGRID)%IWELPB=>IWELPB
        GWFWELDAT(IGRID)%NNPWEL=>NNPWEL
        GWFWELDAT(IGRID)%WELAUX=>WELAUX
        GWFWELDAT(IGRID)%WELL=>WELL
C
      RETURN
      END


text:
! { dg-do compile }
!
! PR fortran/34861, in which the test of conformity of the result array bounds
! would barf because they are not known at compile time in this case.
!
! Contributed by Dick Hendrickson <dick.hendrickson@gmail.com>
!
FUNCTION I_IMFUD0 ( IDA2 , NDS4, NDS3) RESULT(I_IMFUDP)
  INTEGER  ::   NDS4, NDS3
  INTEGER  ::   IDA2(5,NDS4,NDS3,2)
  INTEGER  ::   I_IMFUDP(SIZE(IDA2,1), SIZE(IDA2,2), SIZE(IDA2,3), SIZE(IDA2,4))
  ENTRY I_IMFUDX (NDS4, NDS3, IDA2) RESULT(I_IMFUDP)
  ENTRY I_IMFUDY (NDS3, NDS4, IDA2) RESULT(I_IMFUDP)
  ENTRY I_IMFUDZ (NDS3, IDA2, NDS4) RESULT(I_IMFUDP)
  I_IMFUDP = 1-IDA2(:,:,:,::NDS4-NDS3)
END FUNCTION

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


!*************************************************************************
      SUBROUTINE Light_Tahoe(latit)
!***********************************************************************!
!
!     This routine calculates the light extinction coefficient for the
!     Lake Tahoe, development of the Tahoe light model of Swift (2000)
!
!    scatter = scattering contribution of sediments (clay)
!    a_cdom  = absortion of colored dissolved organic matter
!    a_star  = absortivity coefficient of chlorophyll-material
!    b_star  = scattering coefficient of chlorophyll-material
!    gamma   = constant on equation SD
!    SD      = Secchi Depth
!    mu_zero = cosinus dependence
!    Kd      = diffuse irradiance attenuation coeffient
!    PAR     = Photosynthetically Active Radiation
!***********************************************************************!
      USE DLMWQ_VARIABLES
      USE DLMWQ_SERVICES
		IMPLICIT NONE
      SAVE

		INTEGER(4), PARAMETER:: size = 7
		INTEGER(4):: i,j,j_counter,k,l,month
!IF		INTEGER(4):: jday
		INTEGER(4):: exits

		REAL(8):: Kd
		REAL(8):: A_1(maxns),B_1(maxns),C_1(maxns),BSed_1(maxns)
		REAL(8):: A_mean, B_mean, C_mean, Depth_mean,Bsed_mean,Depth_S
		REAL(8):: Chloro_mean,Parts_mean(size),Volumetric
		REAL(8):: Organic_mean,Total_mean
      REAL(8):: Kd_mean
		REAL(8):: b_sed(maxns)
		REAL(8):: Scatter(size)
		REAL(8):: SD_prima
		REAL(8):: SD_measured(size)
		REAL(8):: POP_mean,PON_mean
		REAL(8):: Vol_mean
      REAL(8):: Nitrate_mean
      REAL(8):: Ammonia_mean
      REAL(8):: THP_mean
		REAL(8):: Total, Organic
		REAL(8):: Dec, radian, tau, elevation, mu_zero, latit
		REAL(8), PARAMETER :: pi=3.1415926535897932384626433832795d0
		REAL(8), PARAMETER :: arfac  = 1.0d+6,volfac = 1000.0d0
		REAL(8):: Ch_abs,Ch_sca,Wa_abs,Wa_sca,CDOMab,Pa_sca,Pa_sc1,Pa_sc2, Pa_sc3,Pa_sc4,Pa_sc5,Pa_sc6
!
!		New Set of Parameters (06/09/2002)
!
		REAL(8),PARAMETER:: 		gamma    = 8.7000D0   			! # 1	! constant, Gamma
		REAL(8),PARAMETER:: 		a_water  = 0.0120D0				! # 5	! /m, Absroption Pure Water
		REAL(8),PARAMETER:: 		a_cdom   = 0.0380D0				! # 2	! /m, Absorption CDOM
		REAL(8),PARAMETER:: 		a_star   = 0.0250D0				! # 4	! m2/mg, Absroption Chla !low value = up the line, high = low the line
		REAL(8),PARAMETER:: 		b_water  = 0.0027D0				! # 5	! /m, Scattering Pure Water
		REAL(8),PARAMETER:: 		b_star   = 0.1050D0				! # 6 ! m2/mg^0.62, Sacttering Chla

		Scatter(1)  = 4.287D-12			! # 7-13 ! m2/particle, Particle Scattering Coefficients
		Scatter(2)  = 3.015D-11       ! m2/particle
		Scatter(3)  = 9.939D-11       ! m2/particle,
		Scatter(4)  = 3.757D-10       ! m2/particle,
		Scatter(5)  = 1.459D-09       ! m2/particle,
		Scatter(6)  = 5.831D-09       ! m2/particle,
		Scatter(7)  = 0.0d0           ! m2/particle,
4004  CONTINUE
!
!    Factor (Ino/Total) particles
!
		fraction  = 0.0D0
		Organic   = 0.0D0
		Total     = 0.0D0

		DO i = 1,size
			Total   = Total   + cf(ns,i)
		ENDDO ! i
		DO i = 1,nchl
			 Organic = Organic + DBLE(phyto_part_fac(i)*wqual(ns,i))
		ENDDO ! i

!
		fraction = Total /(Total+ Organic)
		IF (fraction.ge.1.0D0) fraction = 1.0D0
		IF (fraction.lt.0.0D0) fraction = 0.0D0
!
!
!    Inicialize  variables
!
!
! Solar angle
!
		radian	 = 2.0D0 * pi/360.0D0									! Degrees to radians factor
		jday      = jday - int(jday/1000)*1000							! Jday (XXX) this omitted the year like 1999002
		Dec       = 23.45D0 * dsin((284+jday)*2.0D0*pi/365.0D0)	! Solar Declination
		tau       = DBLE(360.0D0*12.0D0/24.0D0)						! Hour. Fix to midday
		elevation = dsin(dsin(radian*latit)*dsin(radian*dec)	&	! Solar Elevation
		  		  + dcos(radian*latit)*dcos(radian*dec)*dcos(radian*tau))
!gbs	mu_zero   = dcos(dasin(dcos(elevation)/1.33D0))				! Angle mu
		mu_zero   = dcos(dasin(dsin(elevation)/1.33D0))				! Angle mu
		exits=0
!xxxxxxxxxxxxxx LOOP Starts from here xxxxxxxxxxxxxxxxxxxxxxxx
5		CONTINUE

		Total_mean    = 0.0d0
		Organic_mean  = 0.0d0
		fraction_mean = 0.0d0
		DO i = 1,ns
			A_1(i)    = 0.0D0
			B_1(i)    = 0.0D0
			b_sed(i)  = 0.0D0
			BSed_1(i) = 0.0D0
			Kd        = 0.0D0
			DO j = 1,nchl
				A_1(i)    = DBLE(A_1(i)   + DBLE(a_star*wqual(i,j)))
				B_1(i)    = DBLE(B_1(i)   + DBLE(b_star*wqual(i,j)**0.62D0))
			ENDDO ! j

			DO j = 1,7
				b_sed(i)  = DBLE(b_sed(i) +DBLE(DBLE(fraction*Scatter(j)*cf(i,j))))
			ENDDO  !j
			BSed_1(i) = b_sed(i)
			A_1(i)    = DBLE(A_1(i) + a_water + a_cdom)
			B_1(i)    = DBLE(B_1(i) + b_water + b_sed(i))
			C_1(i)    = DBLE(A_1(i) + B_1(i))
			Kd = (A_1(i)/mu_zero) * DSQRT(1.0D0+(0.425D0*mu_zero-0.19D0)* DBLE(B_1(i)/A_1(i)))
			et1(i) = 1.00*Kd
		ENDDO ! i
!
!   Estimate Secchi depth. Use Iterative Process
!
      j_counter = 1
		DO 5000 i = ns,1,-1
			A_mean       = 0.0D0
			B_mean       = 0.0D0
			Bsed_mean    = 0.0D0
			C_mean       = 0.0D0
			Chloro_mean  = 0.0D0
			Depth_mean   = 0.0D0
			Kd_mean      = 0.0D0
			POP_mean     = 0.0D0
			PON_mean     = 0.0D0
			SD_prima     = 0.0D0
			Depth_S      = 0.0D0
			Nitrate_mean = 0.0D0
			Ammonia_mean = 0.0D0
			THP_mean     = 0.0D0
			Vol_mean     = 0.0D0
			Total_mean	  = 0.0D0
			Organic_mean = 0.0d0
			fraction_mean = 0.0d0

			DO l = 1,7
				Parts_mean(l) = 0.0D0
			ENDDO

			DO j = ns,ns - j_counter,-1
				IF(j-1.eq.0.or.ns-j_counter.eq.0) THEN
					Volumetric   = vol(1)*volfac
					Depth_S      = depth(1)
					POP_mean     = POP_mean     + wqual(1,17)*Volumetric
					PON_mean     = PON_mean     + wqual(1,23)*Volumetric
					Nitrate_mean = Nitrate_mean + wqual(1,21)*Volumetric
					Ammonia_mean = Ammonia_mean + wqual(1,22)*Volumetric
					THP_mean     = THP_mean     + wqual(1,16)*Volumetric
					DO k = 1,nchl
						Chloro_mean = Chloro_mean  + wqual(1,k)*Volumetric
					ENDDO    ! k
					DO k = 1,7
						Parts_mean(k) = Parts_mean(k) + cf(1,k)*Volumetric
					ENDDO    ! k
					A_mean =     A_mean     + A_1(1)*depth_S   !
					B_mean =     B_mean     + B_1(1)*depth_S   !
					Bsed_mean  = Bsed_mean  + BSed_1(1)*depth_S
					Depth_mean = Depth_mean + Depth_S
					Vol_mean   = Vol_mean   + Volumetric
					GOTO 3456
				ENDIF
				Volumetric   = vol(j)*volfac
				Depth_S      = depth(j)-depth(j-1)
				POP_mean     = POP_mean     + wqual(j,17)*Volumetric
				PON_mean     = PON_mean     + wqual(j,23)*Volumetric
				Nitrate_mean = Nitrate_mean + wqual(j,21)*Volumetric
				Ammonia_mean = Ammonia_mean + wqual(j,22)*Volumetric
				THP_mean     = THP_mean     + wqual(j,16)*Volumetric
				DO k = 1,nchl
					Chloro_mean = Chloro_mean  + wqual(j,k)*Volumetric
				ENDDO    ! k
				DO k = 1,7
					Parts_mean(k) = Parts_mean(k) + cf(j,k)*Volumetric
				ENDDO    ! k
				A_mean     = A_mean     + A_1(j)   *depth_S !Unit of Length
				B_mean     = B_mean     + B_1(j)   *depth_S !Unit of Length
				Bsed_mean  = Bsed_mean  + BSed_1(j)*depth_S !Unit of Length
				Depth_mean = Depth_mean + Depth_S
				Vol_mean   = Vol_mean   + Volumetric
			ENDDO  ! j

 3456		CONTINUE

			A_mean      = DBLE(A_mean/Depth_mean)
			B_mean      = DBLE(B_mean/Depth_mean)
			Bsed_mean   = DBLE(Bsed_mean/Depth_mean)
			C_mean      = A_mean + B_mean
			Kd_mean     = A_mean * (DSQRT(1+(0.425*mu_zero-0.19D0)*		&
							  DBLE(DBLE(B_mean/A_mean))/mu_zero))

			SD_prima    = DBLE(gamma/DBLE(C_mean + Kd_mean))
			Chloro_mean = DBLE(Chloro_mean/Vol_mean)

			DO k = 1,7
				Parts_mean(k)  = DBLE(Parts_mean(k)/Vol_mean)
			ENDDO       ! k
			Nitrate_mean = DBLE(Nitrate_mean/Vol_mean)
			Ammonia_mean = DBLE(Ammonia_mean/Vol_mean)
			THP_mean     = DBLE(THP_mean/Vol_mean)
			POP_mean     = DBLE(POP_mean/Vol_mean)
			PON_mean     = DBLE(PON_mean/Vol_mean)

			IF(SD_prima.gt.(depth(ns)-depth(i))) THEN
				SD = Depth_mean
				j_counter = j_counter + 1
!				WRITE(*,fmt='(3f10.3)')SD_prima,(depth(ns)-depth(i))
			ELSE
				SD           = SD_prima
				A_Tahoe      = A_mean
				B_Tahoe      = B_mean
				Kd_Tahoe     = Kd_mean
				Bsed_Tahoe   = Bsed_mean
				Chloro_Tahoe = Chloro_mean
				DO k = 1,7
					Parts_Tahoe(k)  = Parts_mean(k)
				ENDDO
				POP_Tahoe     = POP_mean
				PON_Tahoe     = PON_mean
				Nitrate_Tahoe = Nitrate_mean
				Ammonia_Tahoe = Ammonia_mean
				THP_Tahoe     = THP_mean
				DO k = 1,7
					Total_mean   = Total_mean   + Parts_mean(k)
				ENDDO ! k
				Organic_mean = Organic_mean + DBLE(phyto_part_fac(1)*Chloro_mean)
				fraction_mean = Total_mean /(Total_mean+ Organic_mean)
		      IF (fraction_mean.ge.1.0D0) fraction_mean = 1.0D0        !gbs
      		IF (fraction_mean.lt.0.0D0) fraction_mean = 0.0D0	     !gbs
!gbs   Becuase the repeation takes place for some values and it never converge
!gbs   the value has changed. ELSE GOTO 5 repeats and never converge.
!gbs	    IF(DABS(fraction-fraction_mean).gt.0.01) THEN
!gbs       WRITE(*,fmt='(2f10.3)')fraction-fraction_mean,SD
  				IF(DABS(fraction-fraction_mean).gt.0.01) THEN
					fraction = fraction_mean
					exits=exits+1
					IF(exits.gt.5) GOTO 4900
					GOTO 5
				ELSEIF (DABS(fraction-fraction_mean).LE.0.01) THEN
				   GOTO 4900
				ENDIF
				GOTO 4900
			ENDIF
5000	ENDDO  ! i
4900  CONTINUE

!
! Contribution from each source
!
      Ch_abs = 0.0d0
      Ch_sca = 0.0d0
      Wa_abs = 0.0d0
      Wa_sca = 0.0d0
      CDOMab = 0.0d0
      Pa_sca = 0.0d0
      Pa_sc1 = 0.0d0
      Pa_sc2 = 0.0d0
      Pa_sc3 = 0.0d0
      Pa_sc4 = 0.0d0
      Pa_sc5 = 0.0d0
      Pa_sc6 = 0.0d0
		DO k = ns-1, i,-1
			DO j = 1,nchl
				Ch_abs    = Ch_abs   + DBLE(a_star*wqual(k,j))*(depth(k+1)-depth(k))
				Ch_sca    = Ch_sca   + DBLE(b_star*wqual(k,j)**0.62D0)*(depth(k+1)-depth(k))
			ENDDO ! j

			DO j = 1,7
				   Pa_sca = Pa_sca +DBLE(Scatter(j)*cf(k,j))*(depth(k+1)-depth(k))
				 IF (J.EQ.1) THEN
				   Pa_sc1 = Pa_sc1 +DBLE(Scatter(j)*cf(k,j))*(depth(k+1)-depth(k))
				 ELSEIF (J.EQ.2) THEN
				   Pa_sc2 = Pa_sc2 +DBLE(Scatter(j)*cf(k,j))*(depth(k+1)-depth(k))
				 ELSEIF (J.EQ.3) THEN
				   Pa_sc3 = Pa_sc3 +DBLE(Scatter(j)*cf(k,j))*(depth(k+1)-depth(k))
				 ELSEIF (J.EQ.4) THEN
				   Pa_sc4 = Pa_sc4 +DBLE(Scatter(j)*cf(k,j))*(depth(k+1)-depth(k))
				 ELSEIF (J.EQ.5) THEN
				   Pa_sc5 = Pa_sc5 +DBLE(Scatter(j)*cf(k,j))*(depth(k+1)-depth(k))
				 ELSEIF (J.EQ.6) THEN
				   Pa_sc6 = Pa_sc6 +DBLE(Scatter(j)*cf(k,j))*(depth(k+1)-depth(k))
				 ENDIF
			ENDDO  !j
			Wa_abs    =  Wa_abs + a_water*(depth(k+1)-depth(k))
			Wa_sca    =  Wa_sca + b_water*(depth(k+1)-depth(k))
			CDOMab    =  CDOMab + a_cdom*(depth(k+1)-depth(k))
    ENDDO ! k

      IF(ISNAN(SD).or.SD<=10.0d0)SD=10.0d0
      WRITE (99,5005)JDAY,ICLOCK,SD, depth(NS)-DEPTH(I),Ch_abs,Ch_sca,Wa_abs,  &
      Wa_sca,CDOMab,Pa_sca,Pa_sc1,Pa_sc2, Pa_sc3,Pa_sc4,Pa_sc5,Pa_sc6
5005  FORMAT(2I8, 2F10.3, 12F10.6)
		RETURN
		END SUBROUTINE Light_Tahoe


text:
      integer nqcdjets,nqcdstart
      common/nqcdjets/nqcdjets,nqcdstart

max_stars_repo_name: tmartini/JHUGen

max_issues_repo_name: tmartini/JHUGen

max_forks_repo_name: tmartini/JHUGen


text:
C  *********************************************************
C  *                                                       *
C  *    TEST NUMBER: 06.02.01/02                           *
C  *    TEST TITLE : Behavior of 3D mapping                *
C  *                                                       *
C  *    PHIGS Validation Tests, produced by NIST           *
C  *                                                       *
C  *********************************************************

      COMMON /GLOBNU/ CTLHND, ERRSIG, ERRFIL, IERRCT, UNERR,
     1        TESTCT, IFLERR, PASSSW, ERRSW, MAXLIN,
     2        CONID, MEMUN, WKID, WTYPE, GLBLUN, INDLUN,
     3        DUMINT, DUMRL
      INTEGER         CTLHND, ERRSIG, ERRFIL, IERRCT, UNERR,
     1        TESTCT, IFLERR, PASSSW, ERRSW, MAXLIN,
     2        CONID, MEMUN, WKID, WTYPE, GLBLUN, INDLUN,
     3        DUMINT(20), ERRIND
      REAL    DUMRL(20)

      COMMON /GLOBCH/ PIDENT,    GLBERR,    TSTMSG,     FUNCID,
     1                DUMCH
      CHARACTER       PIDENT*40, GLBERR*60, TSTMSG*900, FUNCID*80,
     1                DUMCH(20)*20

C  Declare program-specific variables
      REAL       VWWNLM(4), PJVPLM(6), VWMPMT(4,4)
      REAL       PRP(3), VPD,BPD,FPD, VWCENT(3)
      INTEGER    PJTYPE

C viewtype
      INTEGER    PPARL,     PPERS
      PARAMETER (PPARL = 0, PPERS = 1)

      CALL INITGL ('06.02.01/02')

C open PHIGS
      CALL XPOPPH (ERRFIL, MEMUN)

C  Set 3D mapping parameters for PARALLEL projection:
      VWWNLM(1) = -462.27
      VWWNLM(2) = -80.95
      VWWNLM(3) = -4891.85
      VWWNLM(4) = -94.78

      PJVPLM(1) = 0.2640
      PJVPLM(2) = 0.4082
      PJVPLM(3) = 0.01339
      PJVPLM(4) = 0.7457
      PJVPLM(5) = 0.1395
      PJVPLM(6) = 0.30036

      PJTYPE = PPARL

      PRP(1) = -15.044
      PRP(2) = 1.318
      PRP(3) = 34.873
      VPD = -1.213
      BPD = -3.311
      FPD = 15.332

      CALL SETMSG ('1 2 7 8 9', '<Evaluate view mapping matrix 3> '   //
     1             'should map the corners of the view volume to '    //
     2             'the corresponding corners of the projection '     //
     3             'viewport for a PARALLEL projection.')
C subroutine TSTCOR performs pass/fail test
      CALL TSTCOR  (VWWNLM, PJVPLM, PJTYPE, PRP, VPD,BPD,FPD,
     1              VWMPMT, VWCENT)

      CALL SETMSG ('1 2 10', '<Evaluate view mapping matrix 3> '      //
     1             'should map any two points on a given projector '  //
     2             'to the same x and y values in NPC space for a '   //
     3             'PARALLEL projection.')
C subroutine TSTPRJ performs pass/fail test
      CALL TSTPRJ (PRP, VWMPMT, VWCENT)

C  Set 3D mapping parameters for PERSPECTIVE projection:
      VWWNLM(1) = -62.2236
      VWWNLM(2) = -11.39577
      VWWNLM(3) = -60.530
      VWWNLM(4) =  24.255

      PJVPLM(1) = 0.10639
      PJVPLM(2) = 0.22243
      PJVPLM(3) = 0.10012
      PJVPLM(4) = 0.15392
      PJVPLM(5) = 0.24695
      PJVPLM(6) = 0.88262

      PJTYPE = PPERS

      PRP(1) = -3.6054E-02
      PRP(2) = 2.1181E-03
      PRP(3) = 3.0569

      VPD = 0.16487
      BPD = -5.2858E-03
      FPD = 3.0531

      CALL SETMSG ('1 2 7 8 9', '<Evaluate view mapping matrix 3> '   //
     1             'should map the corners of the view volume to '    //
     2             'the corresponding corners of the projection '     //
     3             'viewport for a PERSPECTIVE projection.')
C subroutine TSTCOR performs pass/fail test
      CALL TSTCOR  (VWWNLM, PJVPLM, PJTYPE, PRP, VPD,BPD,FPD,
     1              VWMPMT, VWCENT)

      CALL SETMSG ('1 2 10', '<Evaluate view mapping matrix 3> '      //
     1             'should map any two points on a given projector '  //
     2             'to the same x and y values in NPC space for a '   //
     3             'PERSPECTIVE projection.')
C subroutine TSTPRJ performs pass/fail test
      CALL TSTPRJ (PRP, VWMPMT, VWCENT)

C  Set 3D mapping parameters for PARALLEL projection, with flat
C  viewport (ZMIN = ZMAX):
      VWWNLM(1) = -2.6073
      VWWNLM(2) = 231.673
      VWWNLM(3) = 830.654
      VWWNLM(4) = 4.74075E+4

      PJVPLM(1) = 0.30322
      PJVPLM(2) = 0.54641
      PJVPLM(3) = 0.28099
      PJVPLM(4) = 0.68383
      PJVPLM(5) = 6.0076E-03
      PJVPLM(6) = 6.0076E-03

      PJTYPE = PPARL

      PRP(1) = 3.8497
      PRP(2) = 72.490
      PRP(3) = -111.363
      VPD = -0.38069
      BPD = -0.46296
      FPD = -4.25281E-02

      CALL SETMSG ('1 2 7 8 9', '<Evaluate view mapping matrix 3> '   //
     1             'should map the corners of the view volume to '    //
     2             'the corresponding corners of a flat projection '  //
     3             'viewport for a PARALLEL projection.')
C subroutine TSTCOR performs pass/fail test
      CALL TSTCOR  (VWWNLM, PJVPLM, PJTYPE, PRP, VPD,BPD,FPD,
     1              VWMPMT, VWCENT)

      CALL SETMSG ('1 2 10', '<Evaluate view mapping matrix 3> '      //
     1             'should map any two points on a given projector '  //
     2             'to the same x and y values in NPC space for a '   //
     3             'flat projection viewport and a PARALLEL '         //
     4             'projection.')
C subroutine TSTPRJ performs pass/fail test
      CALL TSTPRJ (PRP, VWMPMT, VWCENT)

C  Set 3D mapping parameters for PERSPECTIVE projection, with a flat
C  viewport (ZMIN = ZMAX) and view volume (BPD = FPD):
      VWWNLM(1) = -13.7
      VWWNLM(2) =  71.55
      VWWNLM(3) = -51.5509
      VWWNLM(4) =  79.861

      PJVPLM(1) = 0.15580
      PJVPLM(2) = 0.537
      PJVPLM(3) = 1.7186E-03
      PJVPLM(4) = 0.3148
      PJVPLM(5) = 0.73940
      PJVPLM(6) = 0.73940

      PJTYPE = PPERS

      PRP(1) = -1.8005
      PRP(2) = -4.19006E-02
      PRP(3) =  6.3462

      VPD = 0.831
      BPD = -2.1753
      FPD = -2.1753

      CALL SETMSG ('1 2 7 8 9', '<Evaluate view mapping matrix 3> '   //
     1             'should map the corners of a flat view volume to ' //
     2             'the corresponding corners of a flat projection '  //
     3             'viewport for a PERSPECTIVE projection.')
C subroutine TSTCOR performs pass/fail test
      CALL TSTCOR  (VWWNLM, PJVPLM, PJTYPE, PRP, VPD,BPD,FPD,
     1              VWMPMT, VWCENT)

      CALL SETMSG ('1 2 10', '<Evaluate view mapping matrix 3> '      //
     1             'should map any two points on a given projector '  //
     2             'to the same x and y values in NPC space for a '   //
     3             'flat view volume, a flat projection viewport, '   //
     4             'and a PERSPECTIVE projection.')
C subroutine TSTPRJ performs pass/fail test
      CALL TSTPRJ (PRP, VWMPMT, VWCENT)

666   CONTINUE
      CALL ENDIT
      END


text:
great_rsd(aa1,dd1).
great_rsd(g1,h1).
great_rsd(aa1,k1).
great_rsd(aa1,l1).
great_rsd(ee1,jj1).
great_rsd(z1,h1).
great_rsd(kk1,z1).
great_rsd(k1,hh1).
great_rsd(a1,g1).
great_rsd(f1,a1).
great_rsd(a1,u1).
great_rsd(a1,jj1).
great_rsd(cc1,h1).
great_rsd(aa1,h1).
great_rsd(ee1,h1).
great_rsd(g1,x1).
great_rsd(w1,t1).
great_rsd(m1,dd1).
great_rsd(e1,x1).
great_rsd(w1,c1).
great_rsd(t1,kk1).
great_rsd(kk1,x1).
great_rsd(m1,g1).
great_rsd(cc1,jj1).
great_rsd(a1,c1).
great_rsd(m1,hh1).
great_rsd(a1,k1).
great_rsd(c1,dd1).
great_rsd(g1,l1).
great_rsd(e1,q1).
great_rsd(ff1,z1).
great_rsd(ff1,dd1).

max_stars_repo_name: manoelfranca/cilppp

max_issues_repo_name: manoelfranca/cilppp

max_forks_repo_name: manoelfranca/cilppp