! ``mappm.f90`` contains the subroutine ``mappm`` and its dependencies ``cs_profile``, ``cs_limiters``,
! ``ppm_profile``, and ``ppm_limiters``. These are all copied from the FV3GFS module ``fv_mapz.F90``.
! The only code changes that have been made to ``mappm.f90`` compared to the FV3GFS code are:
!   - removed pre-processor ifdef statements
!   - added definitions as needed for the constants r3, r23, r12 within each subroutine
!     (these were defined as global variables in ``fv_mapz.F90``)

!>@brief The subroutine 'mappm' is a general-purpose routine for remapping
!! one set of vertical levels to another.
 subroutine mappm(km, pe1, q1, kn, pe2, q2, i1, i2, iv, kord, ptop)

! IV = 0: constituents
! IV = 1: potential temp
! IV =-1: winds

! Mass flux preserving mapping: q1(im,km) -> q2(im,kn)

! pe1: pressure at layer edges (from model top to bottom surface)
!      in the original vertical coordinate
! pe2: pressure at layer edges (from model top to bottom surface)
!      in the new vertical coordinate

 integer, intent(in):: i1, i2, km, kn, kord, iv
 real, intent(in ):: pe1(i1:i2,km+1), pe2(i1:i2,kn+1) !< pe1: pressure at layer edges from model top to bottom
                                                      !!      surface in the ORIGINAL vertical coordinate
                                                      !< pe2: pressure at layer edges from model top to bottom
                                                      !!      surface in the NEW vertical coordinate
! Mass flux preserving mapping: q1(im,km) -> q2(im,kn)
 real, intent(in )::  q1(i1:i2,km)
 real, intent(out)::  q2(i1:i2,kn)
 real, intent(IN) :: ptop
! local
      real  qs(i1:i2)
      real dp1(i1:i2,km)
      real a4(4,i1:i2,km)
      integer i, k, l
      integer k0, k1
      real pl, pr, tt, delp, qsum, dpsum, esl
      real, parameter:: r3 = 1./3., r23 = 2./3.

      do k=1,km
         do i=i1,i2
             dp1(i,k) = pe1(i,k+1) - pe1(i,k)
            a4(1,i,k) = q1(i,k)
         enddo
      enddo

      if ( kord >7 ) then
           call  cs_profile( qs, a4, dp1, km, i1, i2, iv, kord )
      else
           call ppm_profile( a4, dp1, km, i1, i2, iv, kord )
      endif

!------------------------------------
! Lowest layer: constant distribution
!------------------------------------

      do 5555 i=i1,i2
         k0 = 1
      do 555 k=1,kn

         if(pe2(i,k) .le. pe1(i,1)) then
! above old ptop
            q2(i,k) = q1(i,1)
         elseif(pe2(i,k) .ge. pe1(i,km+1)) then
! Entire grid below old ps
            q2(i,k) = q1(i,km)
         else

         do 45 L=k0,km
! locate the top edge at pe2(i,k)
         if( pe2(i,k) .ge. pe1(i,L) .and.        &
             pe2(i,k) .le. pe1(i,L+1)    ) then
             k0 = L
             PL = (pe2(i,k)-pe1(i,L)) / dp1(i,L)
             if(pe2(i,k+1) .le. pe1(i,L+1)) then

! entire new grid is within the original grid
               PR = (pe2(i,k+1)-pe1(i,L)) / dp1(i,L)
               TT = r3*(PR*(PR+PL)+PL**2)
               q2(i,k) = a4(2,i,L) + 0.5*(a4(4,i,L)+a4(3,i,L)  &
                       - a4(2,i,L))*(PR+PL) - a4(4,i,L)*TT
              goto 555
             else
! Fractional area...
              delp = pe1(i,L+1) - pe2(i,k)
              TT   = r3*(1.+PL*(1.+PL))
              qsum = delp*(a4(2,i,L)+0.5*(a4(4,i,L)+            &
                     a4(3,i,L)-a4(2,i,L))*(1.+PL)-a4(4,i,L)*TT)
              dpsum = delp
              k1 = L + 1
             goto 111
             endif
         endif
45       continue

111      continue
         do 55 L=k1,km
         if( pe2(i,k+1) .gt. pe1(i,L+1) ) then

! Whole layer..

            qsum  =  qsum + dp1(i,L)*q1(i,L)
            dpsum = dpsum + dp1(i,L)
         else
           delp = pe2(i,k+1)-pe1(i,L)
           esl  = delp / dp1(i,L)
           qsum = qsum + delp * (a4(2,i,L)+0.5*esl*            &
                 (a4(3,i,L)-a4(2,i,L)+a4(4,i,L)*(1.-r23*esl)) )
          dpsum = dpsum + delp
           k0 = L
           goto 123
         endif
55       continue
        delp = pe2(i,k+1) - pe1(i,km+1)
        if(delp > 0.) then
! Extended below old ps
           qsum = qsum + delp * q1(i,km)
          dpsum = dpsum + delp
        endif
123     q2(i,k) = qsum / dpsum
      endif
555   continue
5555  continue

 end subroutine mappm

!>@brief The subroutine 'cs_profile' performs the optimized vertical profile reconstruction:
!>@date April 2008
!>@author S. J. Lin, NOAA/GFDL
 subroutine cs_profile(qs, a4, delp, km, i1, i2, iv, kord)
! Optimized vertical profile reconstruction:
! Latest: Apr 2008 S.-J. Lin, NOAA/GFDL
 integer, intent(in):: i1, i2
 integer, intent(in):: km      !< vertical dimension
 integer, intent(in):: iv      !< iv =-1: winds
                               !< iv = 0: positive definite scalars
                               !< iv = 1: others
 integer, intent(in):: kord
 real, intent(in)   ::   qs(i1:i2)
 real, intent(in)   :: delp(i1:i2,km)     !< layer pressure thickness
 real, intent(inout):: a4(4,i1:i2,km)     !< Interpolated values
!-----------------------------------------------------------------------
 logical, dimension(i1:i2,km):: extm, ext5, ext6
 real  gam(i1:i2,km)
 real    q(i1:i2,km+1)
 real   d4(i1:i2)
 real   bet, a_bot, grat
 real   pmp_1, lac_1, pmp_2, lac_2, x0, x1
 integer i, k, im

 if ( iv .eq. -2 ) then
      do i=i1,i2
         gam(i,2) = 0.5
           q(i,1) = 1.5*a4(1,i,1)
      enddo
      do k=2,km-1
         do i=i1, i2
                  grat = delp(i,k-1) / delp(i,k)
                   bet =  2. + grat + grat - gam(i,k)
                q(i,k) = (3.*(a4(1,i,k-1)+a4(1,i,k)) - q(i,k-1))/bet
            gam(i,k+1) = grat / bet
         enddo
      enddo
      do i=i1,i2
            grat = delp(i,km-1) / delp(i,km)
         q(i,km) = (3.*(a4(1,i,km-1)+a4(1,i,km)) - grat*qs(i) - q(i,km-1)) /  &
                   (2. + grat + grat - gam(i,km))
         q(i,km+1) = qs(i)
      enddo
      do k=km-1,1,-1
        do i=i1,i2
           q(i,k) = q(i,k) - gam(i,k+1)*q(i,k+1)
        enddo
      enddo
 else
  do i=i1,i2
         grat = delp(i,2) / delp(i,1)   ! grid ratio
          bet = grat*(grat+0.5)
       q(i,1) = ( (grat+grat)*(grat+1.)*a4(1,i,1) + a4(1,i,2) ) / bet
     gam(i,1) = ( 1. + grat*(grat+1.5) ) / bet
  enddo

  do k=2,km
     do i=i1,i2
           d4(i) = delp(i,k-1) / delp(i,k)
             bet =  2. + d4(i) + d4(i) - gam(i,k-1)
          q(i,k) = ( 3.*(a4(1,i,k-1)+d4(i)*a4(1,i,k)) - q(i,k-1) )/bet
        gam(i,k) = d4(i) / bet
     enddo
  enddo

  do i=i1,i2
         a_bot = 1. + d4(i)*(d4(i)+1.5)
     q(i,km+1) = (2.*d4(i)*(d4(i)+1.)*a4(1,i,km)+a4(1,i,km-1)-a_bot*q(i,km))  &
               / ( d4(i)*(d4(i)+0.5) - a_bot*gam(i,km) )
  enddo

  do k=km,1,-1
     do i=i1,i2
        q(i,k) = q(i,k) - gam(i,k)*q(i,k+1)
     enddo
  enddo
 endif
!----- Perfectly linear scheme --------------------------------
 if ( abs(kord) > 16 ) then
  do k=1,km
     do i=i1,i2
        a4(2,i,k) = q(i,k  )
        a4(3,i,k) = q(i,k+1)
        a4(4,i,k) = 3.*(2.*a4(1,i,k) - (a4(2,i,k)+a4(3,i,k)))
     enddo
  enddo
  return
 endif
!----- Perfectly linear scheme --------------------------------

!------------------
! Apply constraints
!------------------
  im = i2 - i1 + 1

! Apply *large-scale* constraints
  do i=i1,i2
     q(i,2) = min( q(i,2), max(a4(1,i,1), a4(1,i,2)) )
     q(i,2) = max( q(i,2), min(a4(1,i,1), a4(1,i,2)) )
  enddo

  do k=2,km
     do i=i1,i2
        gam(i,k) = a4(1,i,k) - a4(1,i,k-1)
     enddo
  enddo

! Interior:
  do k=3,km-1
     do i=i1,i2
        if ( gam(i,k-1)*gam(i,k+1)>0. ) then
! Apply large-scale constraint to ALL fields if not local max/min
             q(i,k) = min( q(i,k), max(a4(1,i,k-1),a4(1,i,k)) )
             q(i,k) = max( q(i,k), min(a4(1,i,k-1),a4(1,i,k)) )
        else
          if ( gam(i,k-1) > 0. ) then
! There exists a local max
               q(i,k) = max(q(i,k), min(a4(1,i,k-1),a4(1,i,k)))
          else
! There exists a local min
                 q(i,k) = min(q(i,k), max(a4(1,i,k-1),a4(1,i,k)))
               if ( iv==0 ) q(i,k) = max(0., q(i,k))
          endif
        endif
     enddo
  enddo

! Bottom:
  do i=i1,i2
     q(i,km) = min( q(i,km), max(a4(1,i,km-1), a4(1,i,km)) )
     q(i,km) = max( q(i,km), min(a4(1,i,km-1), a4(1,i,km)) )
  enddo

  do k=1,km
     do i=i1,i2
        a4(2,i,k) = q(i,k  )
        a4(3,i,k) = q(i,k+1)
     enddo
  enddo

  do k=1,km
     if ( k==1 .or. k==km ) then
       do i=i1,i2
          extm(i,k) = (a4(2,i,k)-a4(1,i,k)) * (a4(3,i,k)-a4(1,i,k)) > 0.
       enddo
     else
       do i=i1,i2
          extm(i,k) = gam(i,k)*gam(i,k+1) < 0.
       enddo
     endif
     if ( abs(kord) > 9 ) then
       do i=i1,i2
          x0 = 2.*a4(1,i,k) - (a4(2,i,k)+a4(3,i,k))
          x1 = abs(a4(2,i,k)-a4(3,i,k))
          a4(4,i,k) = 3.*x0
          ext5(i,k) = abs(x0) > x1
          ext6(i,k) = abs(a4(4,i,k)) > x1
       enddo
     endif
  enddo

!---------------------------
! Apply subgrid constraints:
!---------------------------
! f(s) = AL + s*[(AR-AL) + A6*(1-s)]         ( 0 <= s  <= 1 )
! Top 2 and bottom 2 layers always use monotonic mapping

  if ( iv==0 ) then
     do i=i1,i2
        a4(2,i,1) = max(0., a4(2,i,1))
     enddo
  elseif ( iv==-1 ) then
      do i=i1,i2
         if ( a4(2,i,1)*a4(1,i,1) <= 0. ) a4(2,i,1) = 0.
      enddo
  elseif ( iv==2 ) then
     do i=i1,i2
        a4(2,i,1) = a4(1,i,1)
        a4(3,i,1) = a4(1,i,1)
        a4(4,i,1) = 0.
     enddo
  endif

  if ( iv/=2 ) then
     do i=i1,i2
        a4(4,i,1) = 3.*(2.*a4(1,i,1) - (a4(2,i,1)+a4(3,i,1)))
     enddo
     call cs_limiters(im, extm(i1,1), a4(1,i1,1), 1)
  endif

! k=2
   do i=i1,i2
      a4(4,i,2) = 3.*(2.*a4(1,i,2) - (a4(2,i,2)+a4(3,i,2)))
   enddo
   call cs_limiters(im, extm(i1,2), a4(1,i1,2), 2)

!-------------------------------------
! Huynh's 2nd constraint for interior:
!-------------------------------------
  do k=3,km-2
     if ( abs(kord)<9 ) then
       do i=i1,i2
! Left  edges
          pmp_1 = a4(1,i,k) - 2.*gam(i,k+1)
          lac_1 = pmp_1 + 1.5*gam(i,k+2)
          a4(2,i,k) = min(max(a4(2,i,k), min(a4(1,i,k), pmp_1, lac_1)),   &
                                         max(a4(1,i,k), pmp_1, lac_1) )
! Right edges
          pmp_2 = a4(1,i,k) + 2.*gam(i,k)
          lac_2 = pmp_2 - 1.5*gam(i,k-1)
          a4(3,i,k) = min(max(a4(3,i,k), min(a4(1,i,k), pmp_2, lac_2)),    &
                                         max(a4(1,i,k), pmp_2, lac_2) )

          a4(4,i,k) = 3.*(2.*a4(1,i,k) - (a4(2,i,k)+a4(3,i,k)))
       enddo

     elseif ( abs(kord)==9 ) then
       do i=i1,i2
          if ( extm(i,k) .and. extm(i,k-1) ) then  ! c90_mp122
! grid-scale 2-delta-z wave detected
               a4(2,i,k) = a4(1,i,k)
               a4(3,i,k) = a4(1,i,k)
               a4(4,i,k) = 0.
          else if ( extm(i,k) .and. extm(i,k+1) ) then  ! c90_mp122
! grid-scale 2-delta-z wave detected
               a4(2,i,k) = a4(1,i,k)
               a4(3,i,k) = a4(1,i,k)
               a4(4,i,k) = 0.
          else
            a4(4,i,k) = 6.*a4(1,i,k) - 3.*(a4(2,i,k)+a4(3,i,k))
! Check within the smooth region if subgrid profile is non-monotonic
            if( abs(a4(4,i,k)) > abs(a4(2,i,k)-a4(3,i,k)) ) then
                  pmp_1 = a4(1,i,k) - 2.*gam(i,k+1)
                  lac_1 = pmp_1 + 1.5*gam(i,k+2)
              a4(2,i,k) = min(max(a4(2,i,k), min(a4(1,i,k), pmp_1, lac_1)),  &
                                             max(a4(1,i,k), pmp_1, lac_1) )
                  pmp_2 = a4(1,i,k) + 2.*gam(i,k)
                  lac_2 = pmp_2 - 1.5*gam(i,k-1)
              a4(3,i,k) = min(max(a4(3,i,k), min(a4(1,i,k), pmp_2, lac_2)),  &
                                             max(a4(1,i,k), pmp_2, lac_2) )
              a4(4,i,k) = 6.*a4(1,i,k) - 3.*(a4(2,i,k)+a4(3,i,k))
            endif
          endif
       enddo
     elseif ( abs(kord)==10 ) then
       do i=i1,i2
          if( ext5(i,k) ) then
              if( ext5(i,k-1) .or. ext5(i,k+1) ) then
                   a4(2,i,k) = a4(1,i,k)
                   a4(3,i,k) = a4(1,i,k)
              elseif ( ext6(i,k-1) .or. ext6(i,k+1) ) then
                   pmp_1 = a4(1,i,k) - 2.*gam(i,k+1)
                   lac_1 = pmp_1 + 1.5*gam(i,k+2)
                   a4(2,i,k) = min(max(a4(2,i,k), min(a4(1,i,k), pmp_1, lac_1)),  &
                                                  max(a4(1,i,k), pmp_1, lac_1) )
                   pmp_2 = a4(1,i,k) + 2.*gam(i,k)
                   lac_2 = pmp_2 - 1.5*gam(i,k-1)
                   a4(3,i,k) = min(max(a4(3,i,k), min(a4(1,i,k), pmp_2, lac_2)),  &
                                                  max(a4(1,i,k), pmp_2, lac_2) )
              endif
          elseif( ext6(i,k) ) then
              if( ext5(i,k-1) .or. ext5(i,k+1) ) then
                  pmp_1 = a4(1,i,k) - 2.*gam(i,k+1)
                  lac_1 = pmp_1 + 1.5*gam(i,k+2)
                  a4(2,i,k) = min(max(a4(2,i,k), min(a4(1,i,k), pmp_1, lac_1)),  &
                                                 max(a4(1,i,k), pmp_1, lac_1) )
                  pmp_2 = a4(1,i,k) + 2.*gam(i,k)
                  lac_2 = pmp_2 - 1.5*gam(i,k-1)
                  a4(3,i,k) = min(max(a4(3,i,k), min(a4(1,i,k), pmp_2, lac_2)),  &
                                                 max(a4(1,i,k), pmp_2, lac_2) )
              endif
          endif
       enddo
       do i=i1,i2
          a4(4,i,k) = 3.*(2.*a4(1,i,k) - (a4(2,i,k)+a4(3,i,k)))
       enddo
     elseif ( abs(kord)==12 ) then
       do i=i1,i2
          if( extm(i,k) ) then
! grid-scale 2-delta-z wave detected
              a4(2,i,k) = a4(1,i,k)
              a4(3,i,k) = a4(1,i,k)
              a4(4,i,k) = 0.
          else        ! not a local extremum
            a4(4,i,k) = 6.*a4(1,i,k) - 3.*(a4(2,i,k)+a4(3,i,k))
! Check within the smooth region if subgrid profile is non-monotonic
            if( abs(a4(4,i,k)) > abs(a4(2,i,k)-a4(3,i,k)) ) then
                  pmp_1 = a4(1,i,k) - 2.*gam(i,k+1)
                  lac_1 = pmp_1 + 1.5*gam(i,k+2)
              a4(2,i,k) = min(max(a4(2,i,k), min(a4(1,i,k), pmp_1, lac_1)),  &
                                             max(a4(1,i,k), pmp_1, lac_1) )
                  pmp_2 = a4(1,i,k) + 2.*gam(i,k)
                  lac_2 = pmp_2 - 1.5*gam(i,k-1)
              a4(3,i,k) = min(max(a4(3,i,k), min(a4(1,i,k), pmp_2, lac_2)),  &
                                             max(a4(1,i,k), pmp_2, lac_2) )
              a4(4,i,k) = 6.*a4(1,i,k) - 3.*(a4(2,i,k)+a4(3,i,k))
            endif
          endif
       enddo
     elseif ( abs(kord)==13 ) then
       do i=i1,i2
          if( ext6(i,k) ) then
             if ( ext6(i,k-1) .and. ext6(i,k+1) ) then
! grid-scale 2-delta-z wave detected
                 a4(2,i,k) = a4(1,i,k)
                 a4(3,i,k) = a4(1,i,k)
             endif
          endif
       enddo
       do i=i1,i2
          a4(4,i,k) = 3.*(2.*a4(1,i,k) - (a4(2,i,k)+a4(3,i,k)))
       enddo
     elseif ( abs(kord)==14 ) then

       do i=i1,i2
          a4(4,i,k) = 3.*(2.*a4(1,i,k) - (a4(2,i,k)+a4(3,i,k)))
       enddo

     elseif ( abs(kord)==15 ) then   ! revised kord=9 scehem
       do i=i1,i2
          if ( ext5(i,k) ) then  ! c90_mp122
             if ( ext5(i,k-1) .or. ext5(i,k+1) ) then  ! c90_mp122
! grid-scale 2-delta-z wave detected
                  a4(2,i,k) = a4(1,i,k)
                  a4(3,i,k) = a4(1,i,k)
             endif
          elseif( ext6(i,k) ) then
! Check within the smooth region if subgrid profile is non-monotonic
                  pmp_1 = a4(1,i,k) - 2.*gam(i,k+1)
                  lac_1 = pmp_1 + 1.5*gam(i,k+2)
              a4(2,i,k) = min(max(a4(2,i,k), min(a4(1,i,k), pmp_1, lac_1)),  &
                                             max(a4(1,i,k), pmp_1, lac_1) )
                  pmp_2 = a4(1,i,k) + 2.*gam(i,k)
                  lac_2 = pmp_2 - 1.5*gam(i,k-1)
              a4(3,i,k) = min(max(a4(3,i,k), min(a4(1,i,k), pmp_2, lac_2)),  &
                                             max(a4(1,i,k), pmp_2, lac_2) )
          endif
       enddo
       do i=i1,i2
          a4(4,i,k) = 3.*(2.*a4(1,i,k) - (a4(2,i,k)+a4(3,i,k)))
       enddo
     elseif ( abs(kord)==16 ) then
       do i=i1,i2
          if( ext5(i,k) ) then
             if ( ext5(i,k-1) .or. ext5(i,k+1) ) then
                 a4(2,i,k) = a4(1,i,k)
                 a4(3,i,k) = a4(1,i,k)
             elseif ( ext6(i,k-1) .or. ext6(i,k+1) ) then
                 ! Left  edges
                 pmp_1 = a4(1,i,k) - 2.*gam(i,k+1)
                 lac_1 = pmp_1 + 1.5*gam(i,k+2)
                 a4(2,i,k) = min(max(a4(2,i,k), min(a4(1,i,k), pmp_1, lac_1)),   &
                                     max(a4(1,i,k), pmp_1, lac_1) )
                 ! Right edges
                 pmp_2 = a4(1,i,k) + 2.*gam(i,k)
                 lac_2 = pmp_2 - 1.5*gam(i,k-1)
                 a4(3,i,k) = min(max(a4(3,i,k), min(a4(1,i,k), pmp_2, lac_2)),    &
                                     max(a4(1,i,k), pmp_2, lac_2) )
             endif
          endif
       enddo
       do i=i1,i2
          a4(4,i,k) = 3.*(2.*a4(1,i,k) - (a4(2,i,k)+a4(3,i,k)))
       enddo
     else      ! kord = 11
       do i=i1,i2
         if ( ext5(i,k) .and. (ext5(i,k-1) .or. ext5(i,k+1)) ) then
! Noisy region:
              a4(2,i,k) = a4(1,i,k)
              a4(3,i,k) = a4(1,i,k)
              a4(4,i,k) = 0.
         else
              a4(4,i,k) = 3.*(2.*a4(1,i,k) - (a4(2,i,k)+a4(3,i,k)))
         endif
       enddo
     endif

! Additional constraint to ensure positivity
     if ( iv==0 ) call cs_limiters(im, extm(i1,k), a4(1,i1,k), 0)

  enddo      ! k-loop

!----------------------------------
! Bottom layer subgrid constraints:
!----------------------------------
  if ( iv==0 ) then
     do i=i1,i2
        a4(3,i,km) = max(0., a4(3,i,km))
     enddo
  elseif ( iv .eq. -1 ) then
      do i=i1,i2
         if ( a4(3,i,km)*a4(1,i,km) <= 0. )  a4(3,i,km) = 0.
      enddo
  endif

  do k=km-1,km
     do i=i1,i2
        a4(4,i,k) = 3.*(2.*a4(1,i,k) - (a4(2,i,k)+a4(3,i,k)))
     enddo
     if(k==(km-1)) call cs_limiters(im, extm(i1,k), a4(1,i1,k), 2)
     if(k== km   ) call cs_limiters(im, extm(i1,k), a4(1,i1,k), 1)
  enddo

 end subroutine cs_profile

 subroutine cs_limiters(im, extm, a4, iv)
 integer, intent(in) :: im
 integer, intent(in) :: iv
 logical, intent(in) :: extm(im)
 real , intent(inout) :: a4(4,im)   !< PPM array
! LOCAL VARIABLES:
 real  da1, da2, a6da
 integer i
 real, parameter :: r12 = 1./12.

 if ( iv==0 ) then
! Positive definite constraint
    do i=1,im
    if( a4(1,i)<=0.) then
        a4(2,i) = a4(1,i)
        a4(3,i) = a4(1,i)
        a4(4,i) = 0.
    else
      if( abs(a4(3,i)-a4(2,i)) < -a4(4,i) ) then
         if( (a4(1,i)+0.25*(a4(3,i)-a4(2,i))**2/a4(4,i)+a4(4,i)*r12) < 0. ) then
! local minimum is negative
             if( a4(1,i)<a4(3,i) .and. a4(1,i)<a4(2,i) ) then
                 a4(3,i) = a4(1,i)
                 a4(2,i) = a4(1,i)
                 a4(4,i) = 0.
             elseif( a4(3,i) > a4(2,i) ) then
                 a4(4,i) = 3.*(a4(2,i)-a4(1,i))
                 a4(3,i) = a4(2,i) - a4(4,i)
             else
                 a4(4,i) = 3.*(a4(3,i)-a4(1,i))
                 a4(2,i) = a4(3,i) - a4(4,i)
             endif
         endif
      endif
    endif
    enddo
 elseif ( iv==1 ) then
    do i=1,im
      if( (a4(1,i)-a4(2,i))*(a4(1,i)-a4(3,i))>=0. ) then
         a4(2,i) = a4(1,i)
         a4(3,i) = a4(1,i)
         a4(4,i) = 0.
      else
         da1  = a4(3,i) - a4(2,i)
         da2  = da1**2
         a6da = a4(4,i)*da1
         if(a6da < -da2) then
            a4(4,i) = 3.*(a4(2,i)-a4(1,i))
            a4(3,i) = a4(2,i) - a4(4,i)
         elseif(a6da > da2) then
            a4(4,i) = 3.*(a4(3,i)-a4(1,i))
            a4(2,i) = a4(3,i) - a4(4,i)
         endif
      endif
    enddo
 else
! Standard PPM constraint
    do i=1,im
      if( extm(i) ) then
         a4(2,i) = a4(1,i)
         a4(3,i) = a4(1,i)
         a4(4,i) = 0.
      else
         da1  = a4(3,i) - a4(2,i)
         da2  = da1**2
         a6da = a4(4,i)*da1
         if(a6da < -da2) then
            a4(4,i) = 3.*(a4(2,i)-a4(1,i))
            a4(3,i) = a4(2,i) - a4(4,i)
         elseif(a6da > da2) then
            a4(4,i) = 3.*(a4(3,i)-a4(1,i))
            a4(2,i) = a4(3,i) - a4(4,i)
         endif
      endif
    enddo
 endif
 end subroutine cs_limiters

 subroutine ppm_profile(a4, delp, km, i1, i2, iv, kord)

! INPUT PARAMETERS:
 integer, intent(in):: iv      !< iv =-1: winds iv = 0: positive definite scalars iv = 1: others iv = 2: temp (if remap_t) and w (iv=-2)
 integer, intent(in):: i1      !< Starting longitude
 integer, intent(in):: i2      !< Finishing longitude
 integer, intent(in):: km      !< Vertical dimension
 integer, intent(in):: kord    !< Order (or more accurately method no.):
                               !
 real , intent(in):: delp(i1:i2,km)     !< Layer pressure thickness

! !INPUT/OUTPUT PARAMETERS:
 real , intent(inout):: a4(4,i1:i2,km)  !< Interpolated values

! DESCRIPTION:
!
!   Perform the piecewise parabolic reconstruction
!
! !REVISION HISTORY:
! S.-J. Lin   revised at GFDL 2007
!-----------------------------------------------------------------------
! local arrays:
      real    dc(i1:i2,km)
      real    h2(i1:i2,km)
      real  delq(i1:i2,km)
      real   df2(i1:i2,km)
      real    d4(i1:i2,km)

! local scalars:
      integer i, k, km1, lmt, it
      real  fac
      real  a1, a2, c1, c2, c3, d1, d2
      real  qm, dq, lac, qmp, pmp

      km1 = km - 1
       it = i2 - i1 + 1

      do k=2,km
         do i=i1,i2
            delq(i,k-1) =   a4(1,i,k) - a4(1,i,k-1)
              d4(i,k  ) = delp(i,k-1) + delp(i,k)
         enddo
      enddo

      do k=2,km1
         do i=i1,i2
                 c1  = (delp(i,k-1)+0.5*delp(i,k))/d4(i,k+1)
                 c2  = (delp(i,k+1)+0.5*delp(i,k))/d4(i,k)
            df2(i,k) = delp(i,k)*(c1*delq(i,k) + c2*delq(i,k-1)) /      &
                                    (d4(i,k)+delp(i,k+1))
            dc(i,k) = sign( min(abs(df2(i,k)),              &
                            max(a4(1,i,k-1),a4(1,i,k),a4(1,i,k+1))-a4(1,i,k),  &
                  a4(1,i,k)-min(a4(1,i,k-1),a4(1,i,k),a4(1,i,k+1))), df2(i,k) )
         enddo
      enddo

!-----------------------------------------------------------
! 4th order interpolation of the provisional cell edge value
!-----------------------------------------------------------

      do k=3,km1
         do i=i1,i2
            c1 = delq(i,k-1)*delp(i,k-1) / d4(i,k)
            a1 = d4(i,k-1) / (d4(i,k) + delp(i,k-1))
            a2 = d4(i,k+1) / (d4(i,k) + delp(i,k))
            a4(2,i,k) = a4(1,i,k-1) + c1 + 2./(d4(i,k-1)+d4(i,k+1)) *    &
                      ( delp(i,k)*(c1*(a1 - a2)+a2*dc(i,k-1)) -          &
                        delp(i,k-1)*a1*dc(i,k  ) )
         enddo
      enddo

!     if(km>8 .and. kord>4) call steepz(i1, i2, km, a4, df2, dc, delq, delp, d4)

! Area preserving cubic with 2nd deriv. = 0 at the boundaries
! Top
      do i=i1,i2
         d1 = delp(i,1)
         d2 = delp(i,2)
         qm = (d2*a4(1,i,1)+d1*a4(1,i,2)) / (d1+d2)
         dq = 2.*(a4(1,i,2)-a4(1,i,1)) / (d1+d2)
         c1 = 4.*(a4(2,i,3)-qm-d2*dq) / ( d2*(2.*d2*d2+d1*(d2+3.*d1)) )
         c3 = dq - 0.5*c1*(d2*(5.*d1+d2)-3.*d1*d1)
         a4(2,i,2) = qm - 0.25*c1*d1*d2*(d2+3.*d1)
! Top edge:
!-------------------------------------------------------
         a4(2,i,1) = d1*(2.*c1*d1**2-c3) + a4(2,i,2)
!-------------------------------------------------------
!        a4(2,i,1) = (12./7.)*a4(1,i,1)-(13./14.)*a4(1,i,2)+(3./14.)*a4(1,i,3)
!-------------------------------------------------------
! No over- and undershoot condition
         a4(2,i,2) = max( a4(2,i,2), min(a4(1,i,1), a4(1,i,2)) )
         a4(2,i,2) = min( a4(2,i,2), max(a4(1,i,1), a4(1,i,2)) )
         dc(i,1) =  0.5*(a4(2,i,2) - a4(1,i,1))
      enddo

! Enforce monotonicity  within the top layer

      if( iv==0 ) then
         do i=i1,i2
            a4(2,i,1) = max(0., a4(2,i,1))
            a4(2,i,2) = max(0., a4(2,i,2))
         enddo
      elseif( iv==-1 ) then
         do i=i1,i2
            if ( a4(2,i,1)*a4(1,i,1) <= 0. ) a4(2,i,1) = 0.
         enddo
      elseif( abs(iv)==2 ) then
         do i=i1,i2
            a4(2,i,1) = a4(1,i,1)
            a4(3,i,1) = a4(1,i,1)
         enddo
      endif

! Bottom
! Area preserving cubic with 2nd deriv. = 0 at the surface
      do i=i1,i2
         d1 = delp(i,km)
         d2 = delp(i,km1)
         qm = (d2*a4(1,i,km)+d1*a4(1,i,km1)) / (d1+d2)
         dq = 2.*(a4(1,i,km1)-a4(1,i,km)) / (d1+d2)
         c1 = (a4(2,i,km1)-qm-d2*dq) / (d2*(2.*d2*d2+d1*(d2+3.*d1)))
         c3 = dq - 2.0*c1*(d2*(5.*d1+d2)-3.*d1*d1)
         a4(2,i,km) = qm - c1*d1*d2*(d2+3.*d1)
! Bottom edge:
!-----------------------------------------------------
         a4(3,i,km) = d1*(8.*c1*d1**2-c3) + a4(2,i,km)
!        dc(i,km) = 0.5*(a4(3,i,km) - a4(1,i,km))
!-----------------------------------------------------
!        a4(3,i,km) = (12./7.)*a4(1,i,km)-(13./14.)*a4(1,i,km-1)+(3./14.)*a4(1,i,km-2)
! No over- and under-shoot condition
         a4(2,i,km) = max( a4(2,i,km), min(a4(1,i,km), a4(1,i,km1)) )
         a4(2,i,km) = min( a4(2,i,km), max(a4(1,i,km), a4(1,i,km1)) )
         dc(i,km) = 0.5*(a4(1,i,km) - a4(2,i,km))
      enddo

! Enforce constraint on the "slope" at the surface

      if( iv==0 ) then
          do i=i1,i2
             a4(2,i,km) = max(0.,a4(2,i,km))
             a4(3,i,km) = max(0.,a4(3,i,km))
          enddo
      elseif( iv<0 ) then
          do i=i1,i2
             if( a4(1,i,km)*a4(3,i,km) <= 0. )  a4(3,i,km) = 0.
          enddo
      endif

   do k=1,km1
      do i=i1,i2
         a4(3,i,k) = a4(2,i,k+1)
      enddo
   enddo

!-----------------------------------------------------------
! f(s) = AL + s*[(AR-AL) + A6*(1-s)]         ( 0 <= s  <= 1 )
!-----------------------------------------------------------
! Top 2 and bottom 2 layers always use monotonic mapping
      do k=1,2
         do i=i1,i2
            a4(4,i,k) = 3.*(2.*a4(1,i,k) - (a4(2,i,k)+a4(3,i,k)))
         enddo
         call ppm_limiters(dc(i1,k), a4(1,i1,k), it, 0)
      enddo

      if(kord >= 7) then
!-----------------------
! Huynh's 2nd constraint
!-----------------------
      do k=2,km1
         do i=i1,i2
! Method#1
!           h2(i,k) = delq(i,k) - delq(i,k-1)
! Method#2 - better
            h2(i,k) = 2.*(dc(i,k+1)/delp(i,k+1) - dc(i,k-1)/delp(i,k-1))  &
                     / ( delp(i,k)+0.5*(delp(i,k-1)+delp(i,k+1)) )        &
                     * delp(i,k)**2
! Method#3
!!!            h2(i,k) = dc(i,k+1) - dc(i,k-1)
         enddo
      enddo

      fac = 1.5           ! original quasi-monotone

      do k=3,km-2
        do i=i1,i2
! Right edges
!        qmp   = a4(1,i,k) + 2.0*delq(i,k-1)
!        lac   = a4(1,i,k) + fac*h2(i,k-1) + 0.5*delq(i,k-1)
!
         pmp   = 2.*dc(i,k)
         qmp   = a4(1,i,k) + pmp
         lac   = a4(1,i,k) + fac*h2(i,k-1) + dc(i,k)
         a4(3,i,k) = min(max(a4(3,i,k), min(a4(1,i,k), qmp, lac)),    &
                                        max(a4(1,i,k), qmp, lac) )
! Left  edges
!        qmp   = a4(1,i,k) - 2.0*delq(i,k)
!        lac   = a4(1,i,k) + fac*h2(i,k+1) - 0.5*delq(i,k)
!
         qmp   = a4(1,i,k) - pmp
         lac   = a4(1,i,k) + fac*h2(i,k+1) - dc(i,k)
         a4(2,i,k) = min(max(a4(2,i,k),  min(a4(1,i,k), qmp, lac)),   &
                     max(a4(1,i,k), qmp, lac))
!-------------
! Recompute A6
!-------------
         a4(4,i,k) = 3.*(2.*a4(1,i,k) - (a4(2,i,k)+a4(3,i,k)))
        enddo
! Additional constraint to ensure positivity when kord=7
         if (iv == 0 .and. kord >= 6 )                      &
             call ppm_limiters(dc(i1,k), a4(1,i1,k), it, 2)
      enddo

      else

         lmt = kord - 3
         lmt = max(0, lmt)
         if (iv == 0) lmt = min(2, lmt)

         do k=3,km-2
            if( kord /= 4) then
              do i=i1,i2
                 a4(4,i,k) = 3.*(2.*a4(1,i,k) - (a4(2,i,k)+a4(3,i,k)))
              enddo
            endif
            if(kord/=6) call ppm_limiters(dc(i1,k), a4(1,i1,k), it, lmt)
         enddo
      endif

      do k=km1,km
         do i=i1,i2
            a4(4,i,k) = 3.*(2.*a4(1,i,k) - (a4(2,i,k)+a4(3,i,k)))
         enddo
         call ppm_limiters(dc(i1,k), a4(1,i1,k), it, 0)
      enddo

 end subroutine ppm_profile

 subroutine ppm_limiters(dm, a4, itot, lmt)

! INPUT PARAMETERS:
      real , intent(in):: dm(*)     !< Linear slope
      integer, intent(in) :: itot      !< Total Longitudes
      integer, intent(in) :: lmt       !< 0: Standard PPM constraint 1: Improved full monotonicity constraint
                                       !< (Lin) 2: Positive definite constraint
                                       !< 3: do nothing (return immediately)
! INPUT/OUTPUT PARAMETERS:
      real , intent(inout) :: a4(4,*)   !< PPM array AA <-- a4(1,i) AL <-- a4(2,i) AR <-- a4(3,i) A6 <-- a4(4,i)
! LOCAL VARIABLES:
      real  qmp
      real  da1, da2, a6da
      real  fmin
      integer i
      real, parameter :: r12 = 1./12.

! Developer: S.-J. Lin

      if ( lmt == 3 ) return

      if(lmt == 0) then
! Standard PPM constraint
      do i=1,itot
      if(dm(i) == 0.) then
         a4(2,i) = a4(1,i)
         a4(3,i) = a4(1,i)
         a4(4,i) = 0.
      else
         da1  = a4(3,i) - a4(2,i)
         da2  = da1**2
         a6da = a4(4,i)*da1
         if(a6da < -da2) then
            a4(4,i) = 3.*(a4(2,i)-a4(1,i))
            a4(3,i) = a4(2,i) - a4(4,i)
         elseif(a6da > da2) then
            a4(4,i) = 3.*(a4(3,i)-a4(1,i))
            a4(2,i) = a4(3,i) - a4(4,i)
         endif
      endif
      enddo

      elseif (lmt == 1) then

! Improved full monotonicity constraint (Lin 2004)
! Note: no need to provide first guess of A6 <-- a4(4,i)
      do i=1, itot
           qmp = 2.*dm(i)
         a4(2,i) = a4(1,i)-sign(min(abs(qmp),abs(a4(2,i)-a4(1,i))), qmp)
         a4(3,i) = a4(1,i)+sign(min(abs(qmp),abs(a4(3,i)-a4(1,i))), qmp)
         a4(4,i) = 3.*( 2.*a4(1,i) - (a4(2,i)+a4(3,i)) )
      enddo

      elseif (lmt == 2) then

! Positive definite constraint
      do i=1,itot
      if( abs(a4(3,i)-a4(2,i)) < -a4(4,i) ) then
      fmin = a4(1,i)+0.25*(a4(3,i)-a4(2,i))**2/a4(4,i)+a4(4,i)*r12
         if( fmin < 0. ) then
         if(a4(1,i)<a4(3,i) .and. a4(1,i)<a4(2,i)) then
            a4(3,i) = a4(1,i)
            a4(2,i) = a4(1,i)
            a4(4,i) = 0.
         elseif(a4(3,i) > a4(2,i)) then
            a4(4,i) = 3.*(a4(2,i)-a4(1,i))
            a4(3,i) = a4(2,i) - a4(4,i)
         else
            a4(4,i) = 3.*(a4(3,i)-a4(1,i))
            a4(2,i) = a4(3,i) - a4(4,i)
         endif
         endif
      endif
      enddo

      endif

 end subroutine ppm_limiters


SUBROUTINE apr_srp (GNSSid, BLKTYP, X_SIDE, Z_SIDE, A_SOLAR, F0)

! ----------------------------------------------------------------------
! SUBROUTINE: apr_srp
! ----------------------------------------------------------------------
! Purpose:
! TBD
! ----------------------------------------------------------------------
! Input arguments:
! - GNSSid       : id of satellite constellation
! - blktyp       : satellite block type
!
! Output arguments:
! - X_SIDE       : area of X-sdie in body-fixed frame
! - Z_SIDE       : area of Z-side in body-fixed frame
! - A_SOLAR      : area of solar panel
! - F0           : net force (computed from cannoball) in D direction
! ----------------------------------------------------------------------
! Author :	John Donovan
!
! Created:	04 Feb 2020
!
! Copyright:  GEOSCIENCE AUSTRALIA
! ----------------------------------------------------------------------

      USE mdl_precision
      IMPLICIT NONE

! ----------------------------------------------------------------------
! Dummy arguments declaration
! ----------------------------------------------------------------------
      CHARACTER (LEN=1), INTENT(IN) :: GNSSid
      CHARACTER (LEN=10), INTENT(IN) :: BLKTYP
! ----------------------------------------------------------------------
! Local variables declaration
! ----------------------------------------------------------------------
! ----------------------------------------------------------------------
! Satellite physical informaiton
! ----------------------------------------------------------------------
      REAL (KIND = prec_q), INTENT(OUT) :: X_SIDE, Z_SIDE
      REAL (KIND = prec_q), INTENT(OUT) :: A_SOLAR, F0
! ----------------------------------------------------------------------

! default init
Z_SIDE = 0.d0
X_SIDE = 0.d0
A_SOLAR = 0.d0
F0 = 0.d0

IF (GNSSid == 'G') THEN
! GPS constellation
! -----------------
! I
         if (TRIM(BLKTYP)=='GPS-I')then
         Z_SIDE = 3.020D0
         X_SIDE = 1.728D0
         A_SOLAR= 6.053D0
         F0 = 4.54d-5

! II and IIA
         else if (TRIM(BLKTYP)=='GPS-II' .or. TRIM(BLKTYP)=='GPS-IIA') then
         Z_SIDE = 2.881D0
         X_SIDE = 2.893D0
         A_SOLAR= 11.871D0
         F0 = 8.695d-5
! IIF
         else if(TRIM(BLKTYP)=='GPS-IIF') then
         Z_SIDE = 5.05D0
         X_SIDE = 4.55D0
         A_SOLAR= 22.25D0
         F0 = 16.7d-5
! IIR
         else if (TRIM(BLKTYP)=='GPS-IIR' .or. TRIM(BLKTYP)=='GPS-IIR-A' .or. &
                  TRIM(BLKTYP)=='GPS-IIR-B'.or.TRIM(BLKTYP)=='GPS-IIR-M') then
         Z_SIDE = 4.25D0
         X_SIDE = 4.11D0
         A_SOLAR= 13.92D0
         F0 = 11.15d-5
! III
         else if (TRIM(BLKTYP)=='GPS-IIIA') then
         Z_SIDE = 4.38D0
         X_SIDE = 6.05D0
         A_SOLAR= 22.25D0
         F0 = 11.0d-5

         else
                 print *,'apr_srp - Unknown block type: ',BLKTYP, ", GNSSid = 'G'"
                 STOP
         end if

ELSE IF (GNSSid == 'R') THEN
! GLONASS constellation
! ---------------------
         if(TRIM(BLKTYP)=='GLO'.or.TRIM(BLKTYP)=='GLO-M'.or.TRIM(BLKTYP)=='GLO-M+'.or.&
            TRIM(BLKTYP)=='GLO-K1A'.or.TRIM(BLKTYP)=='GLO-K1B')then
         Z_SIDE = 1.6620D0
         X_SIDE = 4.200D0
         A_SOLAR= 23.616D0
! GLONASS-K
         if(TRIM(BLKTYP)=='GLO-K1A'.or.TRIM(BLKTYP)=='GLO-K1B') F0 = 10.0d-5
! GLONASS-M
         if(TRIM(BLKTYP)=='GLO'.or.TRIM(BLKTYP)=='GLO-M'.or.TRIM(BLKTYP)=='GLO-M+') F0 = 20.9d-5

         else
                 print *,'apr_srp - Unknown block type: ',BLKTYP, ", GNSSid = 'R'"
                 STOP
         end if

ELSE IF (GNSSid == 'E') THEN
! GALILEO constellation
! ---------------------
         if (TRIM(BLKTYP)=='GAL-1'.or.TRIM(BLKTYP)=='GAL-2') then
         Z_SIDE = 3.002D0
         X_SIDE = 1.323D0
         A_SOLAR= 11.0D0
         F0 = 8.35d-5

         else
                 print *,'apr_srp - Unknown block type: ',BLKTYP, ", GNSSid = 'E'"
                 STOP
         end if

ELSE IF (GNSSid == 'C') THEN
! BDS constellation
! -----------------
         Z_SIDE = 3.96D0
         X_SIDE = 4.5D0
         A_SOLAR= 22.44D0
! BDS GEO
         if(TRIM(BLKTYP)=='BDS-2G') F0 = 50.1d-5
! BDS MEO
         if(TRIM(BLKTYP)=='BDS-2M') F0 = 8.35d-5
! BDS IGSO
         if(TRIM(BLKTYP)=='BDS-2I') F0 = 17.0d-5

         if (BLKTYP(1:3)/='BDS') then
                 print *,'apr_srp - Unknown block type: ',BLKTYP, ", GNSSid = 'C'"
                 STOP
         end if

ELSE IF (GNSSid == 'J') THEN
! QZSS constellation
! ------------------
! QZSS-1
         if (TRIM(BLKTYP)=='QZS-1') then
         Z_SIDE = 5.6D0
         X_SIDE = 9.0D0
         A_SOLAR= 45.0D0
         F0 = 35.0d-5
! QZSS-2I
         else if (TRIM(BLKTYP)=='QZS-2I') then
         Z_SIDE = 5.6D0
         X_SIDE = 10.1D0
         A_SOLAR= 29.8D0
         F0 = 25.0d-5
! QZSS-2G
         else if (TRIM(BLKTYP)=='QZS-2G') then
         Z_SIDE = 5.6D0
         X_SIDE = 10.1D0
         A_SOLAR= 29.8D0
         F0 = 25.0d-5

         else
                 print *,'apr_srp - Unknown block type: ',BLKTYP, ", GNSSid = 'J'"
                 STOP
         end if

END IF

END SUBROUTINE


text:
      SUBROUTINE MESGBC(LUNIN,MESGTYP,ICOMP)

C$$$  SUBPROGRAM DOCUMENTATION BLOCK
C
C SUBPROGRAM:    MESGBC
C   PRGMMR: KEYSER           ORG: NP22       DATE: 2004-06-29
C
C ABSTRACT: THIS SUBROUTINE EXAMINES A BUFR MESSAGE AND RETURNS BOTH
C  THE MESSAGE TYPE FROM SECTION 1 AND A MESSAGE COMPRESSION INDICATOR
C  UNPACKED FROM SECTION 3.  IT OBTAINS THE BUFR MESSAGE VIA TWO
C  DIFFERENT METHODS, BASED UPON THE SIGN OF LUNIN.
C     IF LUNIN IS GREATER THAN ZERO, THIS SUBROUTINE READS AND EXAMINES
C  SECTION 1 OF MESSAGES IN A BUFR FILE IN SEQUENCE UNTIL IT FINDS THE
C  FIRST MESSAGE THAT ACTUALLY CONTAINS REPORT DATA {I.E., BEYOND THE
C  BUFR TABLE (DICTIONARY) MESSAGES AT THE TOP AND, FOR DUMP FILES,
C  BEYOND THE TWO DUMMY MESSAGES CONTAINING THE CENTER TIME AND THE
C  DUMP TIME}.  IT THEN RETURNS THE MESSAGE TYPE AND COMPRESSION
C  INDICATOR FOR THIS FIRST DATA MESSAGE.  IN THIS CASE, THE BUFR FILE
C  SHOULD NOT BE OPENED VIA BUFR ARCHIVE LIBRARY SUBROUTINE OPENBF
C  PRIOR TO CALLING THIS SUBROUTINE.  HOWEVER, THE BUFR FILE MUST BE
C  CONNECTED TO UNIT ABS(LUNIN).  WHEN USED THIS WAY, THIS SUBROUTINE
C  IS IDENTICAL TO BUFR ARCHIVE LIBRARY SUBROUTINE MESGBF EXCEPT MESGBF
C  DOES NOT RETURN ANY INFORMATION ABOUT COMPRESSION AND MESGBF READS
C  UNTIL IT FINDS THE FIRST NON-DICTIONARY MESSAGE REGARDLESS OF
C  WHETHER OR NOT IT CONTAINS ANY REPORTS (I.E., IT WOULD STOP AT THE
C  DUMMY MESSAGE CONTAINING THE CENTER TIME FOR DUMP FILES).
C     THE SECOND METHOD IN WHICH THIS SUBROUTINE CAN BE USED OCCURS
C  WHEN LUNIN IS PASSED IN WITH A VALUE LESS THAN ZERO.  IN THIS CASE,
C  IT SIMPLY RETURNS THE MESSAGE TYPE AND COMPRESSION INDICATOR FOR THE
C  BUFR MESSAGE CURRENTLY STORED IN THE INTERNAL MESSAGE BUFFER (ARRAY
C  MBAY IN COMMON BLOCK /BITBUF/).  IN THIS CASE, THE BUFR FILE
C  CONNECTED TO ABS(LUNIN) MUST HAVE BEEN PREVIOUSLY OPENED FOR INPUT
C  OPERATIONS BY BUFR ARCHIVE LIBRARY SUBROUTINE OPENBF, AND THE BUFR
C  MESSAGE MUST HAVE BEEN READ INTO MEMORY BY BUFR ARCHIVE LIBRARY
C  ROUTINE READMG OR EQUIVALENT.
C
C PROGRAM HISTORY LOG:
C 2003-11-04  D. KEYSER  -- ORIGINAL AUTHOR
C 2004-06-29  D. KEYSER  -- ADDED NEW OPTION TO RETURN MESSAGE TYPE AND
C                           COMPRESSION INDICATOR FOR BUFR MESSAGE
C                           CURRENTLY STORED IN MEMORY (TRIGGERED BY
C                           INPUT ARGUMENT LUNIN LESS THAN ZERO)
C 2004-08-09  J. ATOR    -- MAXIMUM MESSAGE LENGTH INCREASED FROM
C                           20,000 TO 50,000 BYTES
C 2005-11-29  J. ATOR    -- USE IUPBS01, GETLENS AND RDMSGW
C DART $Id$
C
C USAGE:    CALL MESGBC (LUNIN, MESGTYP, ICOMP)
C   INPUT ARGUMENT LIST:
C     LUNIN    - INTEGER: ABSOLUTE VALUE IS FORTRAN LOGICAL UNIT NUMBER
C                FOR BUFR FILE
C                  - IF LUNIN IS GREATER THAN ZERO, THIS SUBROUTINE
C                    READS THROUGH ALL BUFR MESSAGES FROM BEGINNING OF
C                    FILE UNTIL IT FINDS THE FIRST MESSAGE CONTAINING
C                    REPORT DATA
C                  - IF LUNIN IS LESS THAN ZERO, THIS SUBROUTINE
C                    OPERATES ON THE BUFR MESSAGE CURRENTLY STORED IN
C                    MEMORY
C
C   OUTPUT ARGUMENT LIST:
C     MESGTYP  - INTEGER: BUFR MESSAGE TYPE FOR EITHER THE FIRST
C                MESSAGE IN FILE CONTAINING REPORT DATA (IF LUNIN > 0),
C                OR FOR THE MESSAGE CURRENTLY IN MEMORY (IF LUNIN < 0)
C                    -256 = for LUNIN > 0 case only: no messages read
C                           or error reading file
C                     < 0 = for LUNIN > 0 case only: none of the
C                           messages read contain reports; this is the
C                           negative of the message type the last
C                           message read (i.e., -11 indicates the BUFR
C                           file contains only BUFR table messages)
C     ICOMP    - INTEGER: BUFR MESSAGE COMPRESSION SWITCH:
C                      -3 = for LUNIN > 0 case only: BUFR file does not
C                           exist
C                      -2 = for LUNIN > 0 case only: BUFR file does not
C                           contain any report messages
C                      -1 = for LUNIN > 0 case only: cannot determine
C                           if first BUFR message containing report
C                           data is compressed due to error reading
C                           file
C                       0 = BUFR message (either first containing
C                           report data if LUNIN > 0, or that currently
C                           in memory if LUNIN < 0) is NOT compressed
C                       1 = BUFR message (either first containing
C                           report data if LUNIN > 0, or that currently
C                           in memory if LUNIN < 0) IS compressed
C
C   INPUT FILES:
C     UNIT ABS(LUNIN) - BUFR FILE
C
C REMARKS:
C    THIS ROUTINE CALLS:        GETLENS  IUPB     IUPBS01  RDMSGW
C                               STATUS   WRDLEN
C    THIS ROUTINE IS CALLED BY: COPYSB   UFBTAB
C                               Also called by application programs.
C
C ATTRIBUTES:
C   LANGUAGE: FORTRAN 77
C   MACHINE:  PORTABLE TO ALL PLATFORMS
C
C$$$

      INCLUDE 'bufrlib.prm'

      COMMON /BITBUF/ MAXBYT,IBIT,IBAY(MXMSGLD4),MBYT(NFILES),
     .                MBAY(MXMSGLD4,NFILES)

      DIMENSION    MSGS(MXMSGLD4)

C-----------------------------------------------------------------------
C-----------------------------------------------------------------------

      LUNIT = ABS(LUNIN)

C  DETERMINE METHOD OF OPERATION BASED ON SIGN OF LUNIN
C   LUNIN > 0 - REWIND AND LOOK FOR FIRST DATA MESSAGE (ITYPE = 0)
C   LUNIN < 0 - LOOK AT MESSAGE CURRENLY IN MEMORY     (ITYPE = 1)
C  ---------------------------------------------------------------

      ITYPE = 0
      IF(LUNIT.NE.LUNIN) ITYPE = 1

      ICOMP   =   -1
      MESGTYP = -256

      IF(ITYPE.EQ.0) THEN

         IREC    =    0
         NSUB    =    0

C  SINCE OPENBF HAS NOT YET BEEN CALLED, MUST CALL WRDLEN TO GET
C  MACHINE INFO NEEDED LATER
C  -------------------------------------------------------------

         CALL WRDLEN

C  READ PAST ANY BUFR TABLES AND RETURN THE FIRST MESSAGE TYPE FOUND
C  -----------------------------------------------------------------

         REWIND LUNIT

1        CALL RDMSGW(LUNIT,MSGS,IER)
         IF(IER.EQ.-1) GOTO 900
         IF(IER.EQ.-2) GOTO 901

         IREC = IREC + 1

         MESGTYP = IUPBS01(MSGS,'MTYP')

         CALL GETLENS(MSGS,2,LEN0,LEN1,LEN2,L3,L4,L5)
         IPT = LEN0+LEN1+LEN2

         NSUB = IUPB(MSGS,IPT+5,16)

         IF(MESGTYP.EQ.11 .OR. NSUB.EQ.0) GO TO 1

      ELSE

C  RETURN MESSAGE TYPE FOR MESSAGE CURRENTLY STORED IN MEMORY
C  ----------------------------------------------------------

         CALL STATUS(LUNIT,LUN,IL,IM)

         DO I=1,12
           MSGS(I) = MBAY(I,LUN)
         ENDDO

         MESGTYP = IUPBS01(MSGS,'MTYP')

         CALL GETLENS(MSGS,2,LEN0,LEN1,LEN2,L3,L4,L5)
         IPT = LEN0+LEN1+LEN2

      END IF

      ICOMP = 0
      NCMP = IUPB(MSGS,IPT+7,8)

C  BUFR MESSAGES ARE COMPRESSED IF BELOW IS TRUE
C  ---------------------------------------------

      IF(IAND(NCMP,64).GT.0) ICOMP = 1
      if(itype.eq.0)  REWIND LUNIT
      GOTO 100

C  CAN ONLY GET TO STATEMENTS 900 OR 901 WHEN ITYPE = 0
C  ----------------------------------------------------

900   IF(IREC.EQ.0) THEN
         MESGTYP = -256
         ICOMP =     -3
      ELSE
         IF(MESGTYP.GE.0) MESGTYP = -MESGTYP
         ICOMP  = -2
         REWIND LUNIT
      ENDIF
      GOTO 100

901   MESGTYP = -256
      ICOMP =     -1
      REWIND LUNIT

C  EXIT
C  ----

100   RETURN
      END

ext: f

max_stars_repo_name: hkershaw-brown/feature-preprocess

max_issues_repo_name: hkershaw-brown/feature-preprocess

max_forks_repo_name: hkershaw-brown/feature-preprocess


text:
! { dg-do compile }
! Checks the fix for PR34008, in which INTENT(INOUT) was disallowed
! for the first argument of assign_m, whereas both INOUT and OUT
! should be allowed.
!
! Contributed by Harald Anlauf <anlauf@gmx.de>
!
module mo_memory
  implicit none
  type t_mi
     logical       :: alloc = .false.
  end type t_mi
  type t_m
     type(t_mi)    :: i                         ! meta data
     real, pointer :: ptr (:,:,:,:) => NULL ()
  end type t_m

  interface assignment (=)
     module  procedure assign_m
  end interface
contains
  elemental subroutine assign_m (y, x)
    !---------------------------------------
    ! overwrite intrinsic assignment routine
    !---------------------------------------
    type (t_m), intent(inout) :: y
    type (t_m), intent(in)    :: x
    y% i = x% i
    if (y% i% alloc) y% ptr = x% ptr
  end subroutine assign_m
end module mo_memory

module gfcbug74
  use mo_memory, only: t_m, assignment (=)
  implicit none
  type t_atm
     type(t_m) :: m(42)
  end type t_atm
contains
  subroutine assign_atm_to_atm (y, x)
    type (t_atm), intent(inout) :: y
    type (t_atm), intent(in)    :: x
    integer :: i
!   do i=1,42; y% m(i) = x% m(i); end do    ! Works
    y% m = x% m                             ! ICE
  end subroutine assign_atm_to_atm
end module gfcbug74
! { dg-final { cleanup-modules "mo_memory gfcbug74" } }

max_stars_repo_name: maurizioabba/rose

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


text:
c
c-----------------------------------------------------------------------
c
      subroutine p19old(xcen,ycen,rho,u,v,p)
      implicit double precision(a-h,o-z)
c
c     # For the point (x,y), compute the true solution at this point
c
       refrho   = 1.d0
       refmach  = 2.25d0
       refc     = 1.d0

       radius = dsqrt(xcen*xcen + ycen*ycen)
       if (radius .lt. 0.725d0) radius = .725d0
       ratio = (1.d0/radius)**2
       theta = datan2(ycen,xcen)
       rho  =  refrho*
     &         (1.d0+.4d0*.5d0*(refmach**2)*(1.-ratio))**(1./.4d0)
       vel = refc*refmach*(1.d0/radius)
       u = vel*dsin(theta)
       v = -vel*dcos(theta)
       p = (1.d0/1.4d0)*rho**1.4d0

       return
       end


text:
c*************************************************************************
c                            RMVS_STEP_IN_MVPL.F
c*************************************************************************
c Rearrange the order of the planets for planocentric encounter. Also
c moves the arraw from xpltb to xplte
c
c             Input:
c                 ipl            ==>  The planet in the center
c                 nbod           ==>  number of massive bodies (int scalar)
c                 mass           ==>  Heliocentric mass of bodies (real array)
c                 xpl,ypl,zpl    ==>  Heliocentric position of planet wrt time
c                                       (real arrays)
c             Output:
c                 masst          ==>  rearranged masses
c              xpltb,ypltb,zpltb ==>  rearranged pl x's at beginning of step
c                                       (1d real arrays)
c              xplte,yplte,zplte ==>  rearranged pl x's at end of step
c                                       (1d real arrays)
c
c
c Remarks: Adopted from hal's wiscl_fk.f
c Authors:  Hal Levison
c Date:    2/19/93
c Last revision:

      subroutine rmvs_step_in_mvpl(ipl,nbod,mass,xpl,ypl,zpl,
     &             masst,xpltb,ypltb,zpltb,xplte,yplte,zplte)

      include '../swift.inc'
      include 'rmvs.inc'

c...  Inputs Only:
      integer ipl,nbod
      real*8 mass(NPLMAX)
      real*8 xpl(NPLMAX),ypl(NPLMAX)
      real*8 zpl(NPLMAX)

c...  Outputs:
      real*8 xpltb(NPLMAX),ypltb(NPLMAX),zpltb(NPLMAX)
      real*8 xplte(NPLMAX),yplte(NPLMAX),zplte(NPLMAX)
      real*8 masst(NPLMAX)

c...  Internals
      integer i

c----
c...  Executable code

c...  Move things over
      do i = 1,nbod
         xpltb(i) = xplte(i)
         ypltb(i) = yplte(i)
         zpltb(i) = zplte(i)
      enddo

c...  first just move the planets
      do i = 2,nbod
         masst(i) = mass(i)
         xplte(i) = xpl(i) - xpl(ipl)
         yplte(i) = ypl(i) - ypl(ipl)
         zplte(i) = zpl(i) - zpl(ipl)
      enddo

c...  now switch things
      masst(1) = mass(ipl)
      xplte(1) = 0.0
      yplte(1) = 0.0
      zplte(1) = 0.0

      masst(ipl) = mass(1)
      xplte(ipl) = -xpl(ipl)
      yplte(ipl) = -ypl(ipl)
      zplte(ipl) = -zpl(ipl)

      return
      end   ! rmvs_step_in_mvpl
c---------------------------------------------------------------------------

max_stars_repo_name: Simske/exostriker


! { dg-do  run }
! { dg-shouldfail "Dim argument incorrect in SUM intrinsic: is 5, should be between 1 and 2" }

program summation

  integer, parameter :: n1=5, n2=7
  integer, dimension(1:n1,1:n2) :: arr
  integer, dimension(1:n1) :: r2
  integer, dimension(1:n2) :: r1
  integer :: i,j
  character(len=80) :: c1, c2
  character(len=50) :: fmt = '(10I5)'
  do j=1,n2
     do i=1,n1
        arr(i,j) = i+j*10
     end do
  end do

  r1 = sum(arr,dim=1,mask=.true.)
  write (unit=c2, fmt=fmt) r1
  call print_sum(1,c1)
  if (c1 /= c2) STOP 1
  r2 = sum(arr,dim=2,mask=.true.)
  write (unit=c2, fmt=fmt) r2
  call print_sum(2,c1)
  if (c1 /= c2) STOP 2
  call print_sum(5,c1)

contains

  subroutine print_sum(d, c)
    integer, intent(in) :: d
    character(len=80), intent(out) :: c
    write (unit=c, fmt=fmt) sum(arr,dim=d,mask=.true.)
   end subroutine

end


text:
C Copyright(C) 1999-2020 National Technology & Engineering Solutions
C of Sandia, LLC (NTESS).  Under the terms of Contract DE-NA0003525 with
C NTESS, the U.S. Government retains certain rights in this software.
C
C See packages/seacas/LICENSE for details

C=======================================================================
      SUBROUTINE SPMAIN (A, NEUTRL, NAMES, NPTIMS, IPTIMS,
     *  TIMES, NENUM,
     &   XN, YN, ZN, XE, YE, ZE,
     &   IE2ELB, ISEVOK, LIDSP, BLKCOL, IDELB,
     *  MAPEL, MAPND)
C=======================================================================

C   --*** SPMAIN *** (SPLOT) SPLOT main plot routine
C   --   Modified by John Glick - 11/9/88
C   --   Written by Amy Gilkey - revised 03/07/88
C   --
C   --SPMAIN calculates the distances, reads the variables, and displays
C   --the distance-versus-variable plots.
C   --
C   --Parameters:
C   --   A          - IN  - the dynamic memory base array
C   --   NEUTRL     - IN  - the type of neutral file to write.
C   --   NAMES      - IN  - the variable names
C   --   NPTIMS     - IN  - the number of selected time steps
C   --   IPTIMS     - IN  - the selected time steps
C   --   TIMES      - IN  - the database times
C   --   NENUM      - IN  - the selected node/element numbers
C   --   XN,YN,ZN   - IN  - the nodal mesh coordinates
C   --   XE,YE,ZE   - IN  - the element centroid mesh coordinates
C   --   IE2ELB     - IN  - the element block for each element
C   --   ISEVOK     - IN  - the element block variable truth table;
C   --                      variable i of block j exists iff ISEVOK(j,i)
C   --   LIDSP(0:*) - IN  - the indices of the selected variables
C   --                      whose values will be displayed on the plot legend.
C   --                      LIDSP(0) = the number of variables in the list.
C   --                      LIDSP(i) identifies the ith variable in the list.
C   --                      LIDSP(i) > 0, LIDSP(i) is the id of a history var.
C   --                      LIDSP(i) < 0, -LIDSP(i) is the id of a global var.
C   --                      LIDSP(i) = 0, TIME is displayed on the plot legend.
C   --   BLKCOL     - I/O - the user selected colors of the element blocks.
C   --                      BLKCOL(0) = 1 if the user defined material
C   --                                  colors should be used in mesh plots.
C   --                                = -1 if program selected colors should
C   --                                  be used.
C   --                      BLKCOL(i) = the user selected color of element
C   --                                  block i:
C   --                                  -2 - no color selected by user.
C   --                                  -1 - black
C   --                                   0 - white
C   --                                   1 - red
C   --                                   2 - green
C   --                                   3 - yellow
C   --                                   4 - blue
C   --                                   5 - cyan
C   --                                   6 - magenta
C   --
C   --Common Variables:
C   --   Uses NNENUM of /SELNE/
C   --   Uses NSPVAR of /SPVARS/

      include 'params.blk'

      PARAMETER (KLFT=1, KRGT=2, KBOT=3, KTOP=4)

      include 'dbnums.blk'
      include 'selne.blk'
      include 'spvars.blk'

      DIMENSION A(*)
      CHARACTER*(*) NAMES(*)
      INTEGER IPTIMS(NPTIMS)
      REAL TIMES(*)
      INTEGER NENUM(NNENUM)
      REAL XN(*), YN(*), ZN(*)
      REAL XE(*), YE(*), ZE(*)
      INTEGER IE2ELB(*)
      LOGICAL ISEVOK(*)
      INTEGER LIDSP(0:*)
      INTEGER BLKCOL(0:NELBLK)
      INTEGER IDELB(*)
      INTEGER MAPEL(*), MAPND(*)

C   --Use the selected color table
      CALL GRCOLU ('ALTERNATE')

C   --Reserve storage for plot variables

      LPLVAL = NNENUM * NSPVAR * NPTIMS
      CALL MDRSRV ('PLTVAL', KPLVAL, LPLVAL)
      CALL MDSTAT (NERR, MEM)
      IF (NERR .GT. 0) GOTO 100

C   --Read plot variables

      CALL SPREAD (A, NPTIMS, IPTIMS, NENUM, A(KPLVAL))

      CALL MDSTAT (NERR, MEM)
      IF (NERR .GT. 0) GOTO 100

C   --Identify and index element variables with undefined values

      CALL MDRSRV ('IXSEGV', KXSEGV, NSPVAR)
      CALL MDSTAT (NERR, MEM)
      IF (NERR .GT. 0) GOTO 100

      CALL SPDSEG (A, NENUM, IE2ELB, ISEVOK,
     &   NSEGV, A(KXSEGV), KSEGEL)

      CALL MDSTAT (NERR, MEM)
      IF (NERR .GT. 0) GOTO 100

C   --Calculate distances

      CALL MDRSRV ('DIST', KDIST, NNENUM)
      CALL MDSTAT (NERR, MEM)
      IF (NERR .GT. 0) GOTO 100

      IF (NODVAR) THEN
         CALL CALDIS (NNENUM, NENUM, XN, YN, ZN, A(KDIST))
      ELSE
         CALL CALDIS (NNENUM, NENUM, XE, YE, ZE, A(KDIST))
      END IF

C   --Plot requested curves

      IF (NSEGV .GT. 0) THEN
         CALL MDRSRV ('XDIST', KXDIST, NNENUM)
         CALL MDSTAT (NERR, MEM)
         IF (NERR .GT. 0) GOTO 100
      END IF

      CALL SPPLOT (A, NEUTRL, NPTIMS, IPTIMS, TIMES, NENUM,
     &   A(KDIST), A(KPLVAL), NAMES,
     &   A(KXSEGV), A(KSEGEL), A(KXDIST), LIDSP, BLKCOL,
     *  MAPEL, MAPND)

      CALL MDDEL ('PLTVAL')
      CALL MDDEL ('DIST')
      IF (NSEGV .GT. 0) THEN
         CALL MDDEL ('XDIST')
      END IF
      CALL MDDEL ('IXSEGV')
      IF (NSEGV .GT. 0) THEN
         CALL MDDEL ('ISEGEL')
      END IF

  100 CONTINUE
      RETURN
      END

max_stars_repo_name: jschueller/seacas

max_issues_repo_name: jschueller/seacas

max_forks_repo_name: jschueller/seacas


text:
!  Copyright (c) December 5 2006 - , CSC - IT Center for Science Ltd., Finland
!
!  This program is free software; you can redistribute it and/or
!  modify it under the terms of the GNU General Public License
!  as published by the Free Software Foundation; either version 2
!  of the License, or (at your option) any later version.
!
!  This program is distributed in the hope that it will be useful,
!  but WITHOUT ANY WARRANTY; without even the implied warranty of
!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!  GNU General Public License for more details.
!
!  You should have received a copy of the GNU General Public License
!  along with this program (in file fem/GPL-2); if not, write to the
!  Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
!  Boston, MA 02110-1301, USA.
!
!
! Wrappers for a bunch of C functions defined in binio.c to read/write
! from/to a binary stream in a portable (independent of endianness) manner.
!
! The OPTIONAL 'Status' argument, if PRESENT, will be set to > 0 for error (use
! StrErrorF to create meaningful error messages), 0 for success and (for the
! BinRead* procedures) -1 for end-of-file.  If .NOT.PRESENT(Status), then, on
! error, an error message will be written and the program terminated.
!
MODULE BinIO

    USE Kinds

    IMPLICIT NONE

    ! TODO: Use some kind of autoconf magic to set these to the C parameters
    ! SEEK_SET etc. That way we could pass the directly onwards to fseeko in the
    ! C code.
    INTEGER, PARAMETER :: BIN_SEEK_SET = 0, BIN_SEEK_CUR = 1, BIN_SEEK_END = 2

    INTERFACE
        FUNCTION BinFTell( Unit ) BIND(C, NAME="binftell_c")
            USE Kinds
            USE, INTRINSIC :: ISO_C_BINDING
            INTEGER(C_INT), INTENT(IN) :: Unit
            INTEGER(IntOff_k) :: BinFTell
        END FUNCTION BinFTell

        SUBROUTINE BinFSeek( Unit, Offset, Whence ) BIND(C, NAME="binfseek_c")
            USE Kinds
            USE, INTRINSIC :: ISO_C_BINDING
            INTEGER(C_INT), INTENT(IN) :: Unit
            INTEGER(IntOff_k), INTENT(IN) :: Offset
            INTEGER, INTENT(IN) :: Whence
        END SUBROUTINE BinFSeek

        SUBROUTINE BinEndianess(e) BIND(C, NAME="binendianess_c")
          USE, INTRINSIC :: ISO_C_BINDING
          CHARACTER(KIND=C_CHAR), INTENT(OUT) :: e(*)
        END SUBROUTINE BinEndianess

        SUBROUTINE BinSetInputEndianess( Unit, e ) BIND(C, NAME="binsetinputendianess_c")
          USE, INTRINSIC :: ISO_C_BINDING
          INTEGER(C_INT) :: Unit
          CHARACTER(KIND=C_CHAR) :: e(*)
        END SUBROUTINE BinSetInputEndianess

        SUBROUTINE BinOpen_C(Unit,File,FileLen,Action,Status) BIND(C, NAME='binopen_c')
          USE, INTRINSIC :: ISO_C_BINDING
          CHARACTER(KIND=C_CHAR) :: File(*), Action(*)
          INTEGER(C_INT) :: Unit,FileLen,Status
        END SUBROUTINE BinOpen_C

        SUBROUTINE BinClose_C(unit, stat) BIND(C, NAME="binclose_c")
            USE, INTRINSIC :: ISO_C_BINDING
            INTEGER(C_INT) :: unit, stat
        END SUBROUTINE BinClose_C

        SUBROUTINE BinWriteChar_C( Unit, c, Status ) BIND(C, NAME='binwritechar_c')
          USE, INTRINSIC :: ISO_C_BINDING
          CHARACTER(C_CHAR) :: c(*)
          INTEGER(C_INT) :: Unit, Status
        END SUBROUTINE BinWriteChar_C

        SUBROUTINE BinReadInt4_C(Unit, a, Status_) BIND(C, NAME="binreadint4_c")
            USE, INTRINSIC :: ISO_C_BINDING
            INTEGER(C_INT) :: Unit, a, Status_
        END SUBROUTINE BinReadInt4_C

        SUBROUTINE BinWriteInt4_C( Unit, a, Status_ ) BIND(C, NAME="binwriteint4_c")
            USE, INTRINSIC :: ISO_C_BINDING
            INTEGER(C_INT) :: Unit, Status_, a ! Note that this assumes LP64!
        END SUBROUTINE BinWriteInt4_c

        SUBROUTINE BinReadInt8_C(unit, a, stat) BIND(C, NAME="binreadint8_c")
            USE, INTRINSIC :: ISO_C_BINDING
            INTEGER(C_INT) :: unit, stat
            INTEGER(C_INT_LEAST64_T) :: a
        END SUBROUTINE BinReadInt8_C

        SUBROUTINE BinWriteInt8_C(unit, a, stat) BIND(C, NAME="binwriteint8_c")
            USE, INTRINSIC :: ISO_C_BINDING
            INTEGER(C_INT) :: unit, stat
            INTEGER(C_INT_LEAST64_T) :: a
        END SUBROUTINE BinWriteInt8_C

        SUBROUTINE BinReadDouble_C(unit, a, stat) BIND(C, NAME="binreaddouble_c")
            USE, INTRINSIC :: ISO_C_BINDING
            INTEGER(C_INT) :: unit, stat
            REAL(C_DOUBLE) :: a
        END SUBROUTINE BinReadDouble_C

        SUBROUTINE BinWriteDouble_C(unit, a, stat) BIND(C, NAME="binwritedouble_c")
            USE, INTRINSIC :: ISO_C_BINDING
            INTEGER(C_INT) :: unit, stat
            REAL(C_DOUBLE) :: a
        END SUBROUTINE BinWriteDouble_C

        SUBROUTINE BinWriteString_C( Unit, s, len, Status ) BIND(C, NAME='binwritestring_c')
          USE, INTRINSIC :: ISO_C_BINDING
          CHARACTER(C_CHAR) :: s(*)
          INTEGER(C_INT) :: Unit, len, Status
        END SUBROUTINE BinWriteString_C

        SUBROUTINE BinReadString_C( Unit, s, len, Status ) BIND(C, NAME='binreadstring_c')
          USE, INTRINSIC :: ISO_C_BINDING
          CHARACTER(C_CHAR) :: s(*)
          INTEGER(C_INT) :: Unit, len, Status
        END SUBROUTINE BinReadString_C

        SUBROUTINE StrErrorF_C( e,s,len) BIND(C, NAME="strerrorf_c")
          USE, INTRINSIC :: ISO_C_BINDING
          CHARACTER(C_CHAR) :: s(*)
          INTEGER(C_INT) :: e,len
        END SUBROUTINE StrErrorF_C
    END INTERFACE

    PRIVATE :: HandleStatus

CONTAINS

    SUBROUTINE HandleStatus( Status, Status_, MsgPrefix )
        INTEGER, OPTIONAL, INTENT(OUT) :: Status
        INTEGER, INTENT(IN) :: Status_
        CHARACTER(LEN=*), INTENT(IN) :: MsgPrefix
        CHARACTER(LEN=100) :: Msg
        INTEGER, PARAMETER :: STDERR = 0

        IF ( PRESENT( Status ) ) THEN
            Status = Status_
        ELSE
            IF ( Status_ > 0 ) THEN
                CALL StrErrorF( Status_, Msg )
                WRITE( STDERR, * ) TRIM(MsgPrefix) // ": " // TRIM(Msg)
                STOP
            END IF
        END IF
    END SUBROUTINE HandleStatus

    SUBROUTINE BinOpen( Unit, File, Action, Status )
        USE, INTRINSIC :: ISO_C_BINDING
        INTEGER, INTENT(IN) :: Unit
        CHARACTER(*), INTENT(IN) :: File
        CHARACTER(*), INTENT(IN) :: Action ! "write", "append" or "read"
        INTEGER, OPTIONAL, INTENT(OUT) :: Status
        INTEGER :: Status_

        CALL BinOpen_C( Unit, TRIM(File)//C_NULL_CHAR, LEN(TRIM(File))+1, Action, Status_ )
        CALL HandleStatus( Status, Status_, "BINIO: Can't open file " &
                                            // TRIM(File) )
    END SUBROUTINE BinOpen

    SUBROUTINE BinClose( Unit, Status )
        INTEGER, INTENT(IN) :: Unit
        INTEGER, OPTIONAL, INTENT(OUT) :: Status
        INTEGER :: Status_

        CALL BinClose_C( Unit, Status_ )
        CALL HandleStatus( Status, Status_, "BINIO: Can't close file" )
    END SUBROUTINE BinClose

    SUBROUTINE BinWriteInt4( Unit, a, Status )
        INTEGER, INTENT(IN) :: Unit
        INTEGER(Int4_k), INTENT(IN) :: a
        INTEGER, OPTIONAL, INTENT(OUT) :: Status
        INTEGER :: Status_

        CALL BinWriteInt4_C( Unit, a, Status_ )
        CALL HandleStatus( Status, Status_, "BINIO: Error writing Int4" )
    END SUBROUTINE BinWriteInt4

    SUBROUTINE BinReadInt4( Unit, a, Status )
        INTEGER, INTENT(IN) :: Unit
        INTEGER(Int4_k), INTENT(OUT) :: a
        INTEGER, OPTIONAL, INTENT(OUT) :: Status
        INTEGER :: Status_

        CALL BinReadInt4_C( Unit, a, Status_ )
        CALL HandleStatus( Status, Status_, "BINIO: Error reading Int4" )
    END SUBROUTINE BinReadInt4

    SUBROUTINE BinWriteInt8( Unit, a, Status )
        INTEGER, INTENT(IN) :: Unit
        INTEGER(Int8_k), INTENT(IN) :: a
        INTEGER, OPTIONAL, INTENT(OUT) :: Status
        INTEGER :: Status_

        CALL BinWriteInt8_C( Unit, a, Status_ )
        CALL HandleStatus( Status, Status_, "BINIO: Error writing Int8" )
    END SUBROUTINE BinWriteInt8

    SUBROUTINE BinReadInt8( Unit, a, Status )
        INTEGER, INTENT(IN) :: Unit
        INTEGER(Int8_k), INTENT(OUT) :: a
        INTEGER, OPTIONAL, INTENT(OUT) :: Status
        INTEGER :: Status_

        CALL BinReadInt8_C( Unit, a, Status_ )
        CALL HandleStatus( Status, Status_, "BINIO: Error reading Int8" )
    END SUBROUTINE BinReadInt8

    SUBROUTINE BinWriteDouble( Unit, a, Status )
        INTEGER, INTENT(IN) :: Unit
        DOUBLE PRECISION, INTENT(IN) :: a
        INTEGER, OPTIONAL, INTENT(OUT) :: Status
        INTEGER :: Status_

        CALL BinWriteDouble_C( Unit, a, Status_ )
        CALL HandleStatus( Status, Status_, "BINIO: Error writing Double" )
    END SUBROUTINE BinWriteDouble

    SUBROUTINE BinReadDouble( Unit, a, Status )
        INTEGER, INTENT(IN) :: Unit
        DOUBLE PRECISION, INTENT(OUT) :: a
        INTEGER, OPTIONAL, INTENT(OUT) :: Status
        INTEGER :: Status_

        CALL BinReadDouble_C( Unit, a, Status_ )
        CALL HandleStatus( Status, Status_, "BINIO: Error reading Double" )
    END SUBROUTINE BinReadDouble

    ! Write a CHARACTER(1).  (Note: to read a CHARACTER(1), with no '\0' at the
    ! end, just use BinReadString with a CHARACTER(1) as argument.
    SUBROUTINE BinWriteChar( UNIT, c, Status )
        INTEGER, INTENT(IN) :: Unit
        CHARACTER, INTENT(IN) :: c
        INTEGER, OPTIONAL, INTENT(OUT) :: Status
        INTEGER :: Status_

        CALL BinWriteChar_C( Unit, c, Status_ )
        CALL HandleStatus( Status, Status_, "BINIO: Error writing char" )
    END SUBROUTINE BinWriteChar

    ! Write 's' to file pointed to by 'unit', and append a '\0'.
    SUBROUTINE BinWriteString( UNIT, s, Status )
        INTEGER, INTENT(IN) :: Unit
        CHARACTER(*), INTENT(IN) :: s
        INTEGER, OPTIONAL, INTENT(OUT) :: Status
        INTEGER :: Status_

        CALL BinWriteString_C( Unit, s, LEN(s), Status_ )
        CALL HandleStatus( Status, Status_, "BINIO: Error writing string" )
    END SUBROUTINE BinWriteString

    ! Read bytes up to next NULL byte (but no more than len(s)) into 's'
    ! from file pointed to by 'unit', and pad with spaces. The NULL byte
    ! will be removed from the stream, but not added to 's'.
    SUBROUTINE BinReadString( Unit, s, Status )
        INTEGER, INTENT(IN) :: Unit
        CHARACTER(*), INTENT(OUT) :: s
        INTEGER, OPTIONAL, INTENT(OUT) :: Status
        INTEGER :: Status_

        CALL BinReadString_C( Unit, s, LEN(s), Status_ )
        CALL HandleStatus( Status, Status_, "BINIO: Error reading string" )
    END SUBROUTINE BinReadString

    ! Return a string representation of the error code 'e' (as returned in a
    ! status variable) in 's'. (Actually just a wrapper to the C function
    ! strerror().)
    SUBROUTINE StrErrorF( e, s )
        INTEGER, INTENT(IN) :: e
        CHARACTER(*), INTENT(OUT) :: s

        CALL StrErrorF_C( e, s, LEN(s) )
    END SUBROUTINE StrErrorF

END MODULE BinIO

max_stars_repo_name: readex-eu/readex-apps

max_issues_repo_name: jcmcmurry/pipelining

max_forks_repo_name: jcmcmurry/pipelining


!
!MIT License
!
!Copyright (c) 2020 flavianowilliams
!
!Permission is hereby granted, free of charge, to any person obtaining a copy
!of this software and associated documentation files (the "Software"), to deal
!in the Software without restriction, including without limitation the rights
!to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
!copies of the Software, and to permit persons to whom the Software is
!furnished to do so, subject to the following conditions:
!
!THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
!IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
!FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
!AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
!LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
!OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
!SOFTWARE.
!
module structure_module
  !*******************************************************************************************
  !*******************************************************************************************

  use constants_module

  implicit none

  integer i

  private
  public :: structure

  type, extends(constants) :: structure
     integer, private           :: nmol
     integer, private           :: natom
     integer, private           :: nfree
     character(9), private      :: gsym
     integer, allocatable       :: ntmol(:)
     integer, allocatable       :: nxmol(:)
     character(10), allocatable :: namemol(:)
     real(8), private           :: volume
     real(8), private           :: a
     real(8), private           :: b
     real(8), private           :: c
     real(8), private           :: alpha
     real(8), private           :: beta
     real(8), private           :: gamma
     real(8), allocatable       :: xa(:)
     real(8), allocatable       :: ya(:)
     real(8), allocatable       :: za(:)
     real(8)                    :: v(3,3)
     real(8)                    :: sys_shift(3)
   contains
     procedure          :: sites
     procedure          :: translate
     procedure          :: set_lattice_constants
     procedure          :: set_lattice_angles
     procedure          :: set_symmetry
     procedure          :: set_nmol
     procedure          :: get_nmol
     procedure          :: set_natom
     procedure          :: get_natom
     procedure          :: set_namemol
     procedure          :: set_ntmol
     procedure          :: set_nxmol
     procedure          :: set_volume
     procedure          :: get_volume
     procedure          :: set_volume2
     procedure          :: set_nfree
     procedure          :: get_nfree
     procedure          :: set_a
     procedure          :: get_a
     procedure          :: set_b
     procedure          :: get_b
     procedure          :: set_c
     procedure          :: get_c
     procedure          :: get_alpha
     procedure          :: get_beta
     procedure          :: get_gamma
     procedure          :: get_gsym
     procedure          :: ccp
     procedure          :: mic
     procedure          :: random_coordinates
  end type structure

contains

  subroutine translate(this)
    class(structure), intent(inout) :: this
    do i=1,this%get_natom()
       this%xa(i)=this%xa(i)+1.0d0*this%sys_shift(1)*this%a
       this%ya(i)=this%ya(i)+1.0d0*this%sys_shift(2)*this%b
       this%za(i)=this%za(i)+1.0d0*this%sys_shift(3)*this%c
    end do
  end subroutine translate

  subroutine set_nmol(this,nmol)
    class(structure), intent(inout) :: this
    integer, intent(in)             :: nmol
    this%nmol=nmol
  end subroutine set_nmol

  integer function get_nmol(this)
    class(structure), intent(in) :: this
    get_nmol=this%nmol
  end function get_nmol

  subroutine set_namemol(this,namemol)
    class(structure), intent(inout) :: this
    integer                         :: i
    character(10), intent(in)       :: namemol(:)
    do i=1,this%get_nmol()
       this%namemol(i)=namemol(i)
    end do
  end subroutine set_namemol

  subroutine set_ntmol(this,ntmol)
    class(structure), intent(inout) :: this
    integer                         :: i
    integer, intent(in)             :: ntmol(:)
    do i=1,this%get_nmol()
       this%ntmol(i)=ntmol(i)
    end do
  end subroutine set_ntmol

  subroutine set_nxmol(this,nxmol)
    class(structure), intent(inout) :: this
    integer                         :: i
    integer, intent(in)             :: nxmol(:)
    do i=1,this%get_nmol()
       this%nxmol(i)=nxmol(i)
    end do
  end subroutine set_nxmol

  subroutine set_volume(this)
    implicit none
    class(structure), intent(inout) :: this
    real(8)                         :: volume,vl(3)
    vl(1)=(this%v(1,2)*this%v(2,3)-this%v(1,3)*this%v(2,2))
    vl(2)=(this%v(1,3)*this%v(2,1)-this%v(1,1)*this%v(2,3))
    vl(3)=(this%v(1,1)*this%v(2,2)-this%v(1,2)*this%v(2,1))
    volume=0.d0
    do i=1,3
       volume=volume+this%v(3,i)*vl(i)
    end do
    this%volume=abs(volume)
  end subroutine set_volume

  subroutine set_volume2(this,volume)
    implicit none
    class(structure), intent(inout) :: this
    real(8), intent(in)             :: volume
    this%volume=volume
  end subroutine set_volume2

  double precision function get_volume(this)
    class(structure), intent(in) :: this
    get_volume=this%volume
  end function get_volume

  subroutine set_a(this,a)
    implicit none
    class(structure), intent(inout) :: this
    real(8), intent(in)             :: a
    this%a=a
  end subroutine set_a

  double precision function get_a(this)
    class(structure), intent(in) :: this
    get_a=this%a
  end function get_a

  subroutine set_b(this,b)
    implicit none
    class(structure), intent(inout) :: this
    real(8), intent(in)             :: b
    this%b=b
  end subroutine set_b

  double precision function get_b(this)
    class(structure), intent(in) :: this
    get_b=this%b
  end function get_b

  subroutine set_c(this,c)
    implicit none
    class(structure), intent(inout) :: this
    real(8), intent(in)             :: c
    this%c=c
  end subroutine set_c

  double precision function get_c(this)
    class(structure), intent(in) :: this
    get_c=this%c
  end function get_c

  double precision function get_alpha(this)
    class(structure), intent(in) :: this
    get_alpha=this%alpha
  end function get_alpha

  double precision function get_beta(this)
    class(structure), intent(in) :: this
    get_beta=this%beta
  end function get_beta

  double precision function get_gamma(this)
    class(structure), intent(in) :: this
    get_gamma=this%gamma
  end function get_gamma

  character(9) function get_gsym(this)
    class(structure), intent(in) :: this
    get_gsym=this%gsym
  end function get_gsym

  subroutine set_natom(this)
    class(structure), intent(inout) :: this
    integer                      :: nx
    nx=0
    do i=1,this%nmol
       nx=nx+this%ntmol(i)*this%nxmol(i)
    end do
    this%natom=nx
  end subroutine set_natom

  function get_natom(this)
    class(structure), intent(inout) :: this
    integer                      :: get_natom
    get_natom=this%natom
  end function get_natom

  subroutine set_nfree(this)
    class(structure), intent(inout) :: this
    this%nfree=3*this%natom-3
  end subroutine set_nfree

  integer function get_nfree(this)
    class(structure), intent(inout) :: this
    get_nfree=this%nfree
  end function get_nfree

  subroutine sites(this)
    class(structure), intent(inout) :: this
    integer                        :: nx
    character(2)                   :: at
    nx=this%get_natom()
    allocate(this%xa(nx),this%ya(nx),this%za(nx))
    read(9,*)nx
    read(9,*)
    if(nx.ne.this%get_natom())goto 1
    do i=1,this%get_natom()
       read(9,*)at,this%xa(i),this%ya(i),this%za(i)
    end do
    return
1   write(6,*)&
         'ERROR: The number of atoms in HICOLM.xyz does not match with that defined in the INPUT file!'
    write(6,*)'Hint: Change the HICOLM.xyz or check the input in the &SYSTEM section.'
    stop
  end subroutine sites

  subroutine set_lattice_constants(this)
    class(structure), intent(inout) :: this
    this%a=sqrt(this%v(1,1)**2+this%v(1,2)**2+this%v(1,3)**2)
    this%b=sqrt(this%v(2,1)**2+this%v(2,2)**2+this%v(2,3)**2)
    this%c=sqrt(this%v(3,1)**2+this%v(3,2)**2+this%v(3,3)**2)
  end subroutine set_lattice_constants

  subroutine set_lattice_angles(this)
    class(structure), intent(inout) :: this
    real(8)                      :: sum
    sum=0.d0
    do i=1,3
       sum=sum+this%v(2,i)*this%v(3,i)
    end do
    this%alpha=acos(sum/(this%b*this%c))
    sum=0.d0
    do i=1,3
       sum=sum+this%v(1,i)*this%v(3,i)
    end do
    this%beta=acos(sum/(this%a*this%c))
    sum=0.d0
    do i=1,3
       sum=sum+this%v(1,i)*this%v(2,i)
    end do
    this%gamma=acos(sum/(this%a*this%b))
  end subroutine set_lattice_angles

  subroutine set_symmetry(this)
    class(structure), intent(inout) :: this
    real(8)                      :: prec
    character(9)                 :: cvar
    prec=1.d-2
    !tryclinic
    cvar='Tryclinic'
    !-cubic
    if(abs(2.d0*this%alpha-this%get_pi()).le.prec)then
       if(abs(2.d0*this%beta-this%get_pi()).le.prec)then
          if(abs(2.d0*this%gamma-this%get_pi()).le.prec)then
             if(abs(this%a-this%b).le.prec)then
                if(abs(this%a-this%c).le.prec)cvar='Cubic'
             end if
          end if
       end if
    end if
    !-Hexagonal
    if(abs(2.d0*this%alpha-acos(-1.d0)).le.prec)then
       if(abs(2.d0*this%beta-acos(-1.d0)).le.prec)then
          if(abs(this%a-this%b).le.prec)then
             if(abs(2.d0*this%gamma-this%get_pi()/6.d0).le.prec.or.&
                  abs(2.d0*this%gamma-this%get_pi()/3.d0).le.prec)cvar='Hexagonal'
          end if
       end if
    end if
    this%gsym=cvar
  end subroutine set_symmetry

  subroutine ccp(this)
    implicit none
    class(structure), intent(inout) :: this
    integer                         :: i
    real(8)                         :: xvz,yvz,zvz
    do i=1,this%get_natom()
       xvz=this%v(1,1)+this%v(2,1)+this%v(3,1)
       yvz=this%v(1,2)+this%v(2,2)+this%v(3,2)
       zvz=this%v(1,3)+this%v(2,3)+this%v(3,3)
       this%xa(i)=this%xa(i)-xvz*nint(this%xa(i)/xvz)
       this%ya(i)=this%ya(i)-yvz*nint(this%ya(i)/yvz)
       this%za(i)=this%za(i)-zvz*nint(this%za(i)/zvz)
    end do
  end subroutine ccp

  subroutine mic(this,i1,i2,xvz,yvz,zvz)
    implicit none
    class(structure), intent(inout) :: this
    integer, intent(in)             :: i1,i2
    real(8), intent(out)            :: xvz,yvz,zvz
    real(8)                         :: xx,yy,zz
    xx=this%v(1,1)+this%v(2,1)+this%v(3,1)
    yy=this%v(1,2)+this%v(2,2)+this%v(3,2)
    zz=this%v(1,3)+this%v(2,3)+this%v(3,3)
    xvz=(this%xa(i2)-this%xa(i1))-xx*int(2.d0*(this%xa(i2)-this%xa(i1))/xx)
    yvz=(this%ya(i2)-this%ya(i1))-yy*int(2.d0*(this%ya(i2)-this%ya(i1))/yy)
    zvz=(this%za(i2)-this%za(i1))-zz*int(2.d0*(this%za(i2)-this%za(i1))/zz)
  end subroutine mic

  subroutine random_coordinates(this)
    implicit none
    class(structure), intent(inout) :: this
    integer                         :: i,j,k,nx
    real(8)                         :: srand,vx,vy,vz,vr,dr
    dr=0.1d0
    nx=1
    do i=1,this%get_nmol()
       do j=1,this%ntmol(i)
          call random_number(srand)
          vx=(2.0d0*srand-1.0d0)*dr
          call random_number(srand)
          vy=(2.0d0*srand-1.0d0)*dr
          call random_number(srand)
          vz=(2.0d0*srand-1.0d0)*dr
          vr=sqrt(vx**2+vy**2+vz**2)
          vx=vx*dr/vr
          vy=vy*dr/vr
          vz=vz*dr/vr
          do k=1,this%nxmol(i)
             this%xa(nx)=this%xa(nx)+vx
             this%ya(nx)=this%ya(nx)+vy
             this%za(nx)=this%za(nx)+vz
             nx=nx+1
          end do
       end do
    end do
  end subroutine random_coordinates

end module structure_module


text:
!     test case insensitive for symobls
! i should be treated as I, not another function i
! Liao, 5/26/2009
      SUBROUTINE A28
      IMPLICIT NONE
      INTEGER I
      i = 1
      END SUBROUTINE A28


! ###############################################################################################################################
! Begin MIT license text.
! _______________________________________________________________________________________________________

! Copyright 2019 Dr William R Case, Jr (dbcase29@gmail.com)

! Permission is hereby granted, free of charge, to any person obtaining a copy of this software and
! associated documentation files (the "Software"), to deal in the Software without restriction, including
! without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
! copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to
! the following conditions:

! The above copyright notice and this permission notice shall be included in all copies or substantial
! portions of the Software and documentation.

! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
! OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
! FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
! AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
! LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
! OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
! THE SOFTWARE.
! _______________________________________________________________________________________________________

! End MIT license text.

   MODULE LINK4_Interface

   INTERFACE

      SUBROUTINE LINK4

      USE PENTIUM_II_KIND, ONLY       :  BYTE, LONG, DOUBLE
      USE IOUNT1, ONLY                :  WRT_BUG, WRT_ERR, WRT_LOG, ERR, ERRSTAT, F04, F06, L1M, L3A, SC1
      USE IOUNT1, ONLY                :  LINK1M,  LINK2I,  LINK3A, L1M_MSG, L3A_MSG
      USE SCONTR, ONLY                :  BLNK_SUB_NAM, COMM, FATAL_ERR, LINKNO, MBUG, NDOFL,                                       &
                                         NTERM_KLL, NTERM_KLLD, NTERM_KLLDn,                                                       &
                                         NTERM_MLL, NTERM_MLLn,                                                                    &
                                         NVEC, NUM_EIGENS, NUM_KLLD_DIAG_ZEROS, NUM_MLL_DIAG_ZEROS, SOL_NAME, WARN_ERR
      USE TIMDAT, ONLY                :  YEAR, MONTH, DAY, HOUR, MINUTE, SEC, SFRAC, STIME, TSEC
      USE CONSTANTS_1, ONLY           :  ZERO, ONE
      USE PARAMS, ONLY                :  EPSIL, LANCMETH, SPARSTOR, SUPINFO
      USE MODEL_STUF, ONLY            :  EIG_COMP, EIG_CRIT, EIG_FRQ1, EIG_FRQ2, EIG_GRID, EIG_METH, EIG_MSGLVL, EIG_LAP_MAT_TYPE, &
                                         EIG_MODE, EIG_N1, EIG_N2, EIG_NCVFACL, EIG_NORM, EIG_SID, EIG_SIGMA, EIG_VECS, MAXMIJ,    &
                                         MIJ_COL, MIJ_ROW, NUM_FAIL_CRIT

      USE SPARSE_MATRICES, ONLY       :  I_KLL, J_KLL, KLL, I_KLLD, J_KLLD, KLLD, I_KLLDn, J_KLLDn, KLLDn,                         &
                                         I_MLL, J_MLL, MLL, I_MLLn, J_MLLn, MLLn
      USE OUTPUT4_MATRICES, ONLY      :  NUM_OU4_REQUESTS
      USE EIGEN_MATRICES_1, ONLY      :  GEN_MASS, MODE_NUM, EIGEN_VAL, EIGEN_VEC
      USE LAPACK_DPB_MATRICES, ONLY   :  ABAND, BBAND
      USE DEBUG_PARAMETERS, ONLY      :  DEBUG

      IMPLICIT NONE

      CHARACTER, PARAMETER            :: CR13 = CHAR(13)   ! This causes a carriage return simulating the "+" action in a FORMAT

      INTEGER(LONG), PARAMETER        :: P_LINKNO = 2        ! Prior LINK no's that should have run before this LINK can execute.

      END SUBROUTINE LINK4

   END INTERFACE

   END MODULE LINK4_Interface


text:
!
! Copyright (c) Stanford University, The Regents of the University of
!               California, and others.
!
! All Rights Reserved.
!
! See Copyright-SimVascular.txt for additional details.
!
! Permission is hereby granted, free of charge, to any person obtaining
! a copy of this software and associated documentation files (the
! "Software"), to deal in the Software without restriction, including
! without limitation the rights to use, copy, modify, merge, publish,
! distribute, sublicense, and/or sell copies of the Software, and to
! permit persons to whom the Software is furnished to do so, subject
! to the following conditions:
!
! The above copyright notice and this permission notice shall be included
! in all copies or substantial portions of the Software.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
! IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
! PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
! OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
! EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
! PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
! PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
! LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
! NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
! SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
!
!--------------------------------------------------------------------
!
!     This is for solving heat equation in a solid (simple diffusion/
!     Laplace equation)
!
!--------------------------------------------------------------------

      SUBROUTINE CONSTRUCT_HEATS(lM, Ag, Yg)
      USE COMMOD
      USE ALLFUN
      IMPLICIT NONE
      TYPE(mshType), INTENT(IN) :: lM
      REAL(KIND=RKIND), INTENT(IN) :: Ag(tDof,tnNo), Yg(tDof,tnNo)

      INTEGER(KIND=IKIND) a, e, g, Ac, eNoN, cPhys
      REAL(KIND=RKIND) w, Jac, ksix(nsd,nsd)

      INTEGER(KIND=IKIND), ALLOCATABLE :: ptr(:)
      REAL(KIND=RKIND), ALLOCATABLE :: xl(:,:), al(:,:), yl(:,:),
     2   N(:), Nx(:,:), lR(:,:), lK(:,:,:)

      eNoN = lM%eNoN

!     HEATS: dof = 1
      ALLOCATE(ptr(eNoN), xl(nsd,eNoN), al(tDof,eNoN), yl(tDof,eNoN),
     2   N(eNoN), Nx(nsd,eNoN), lR(dof,eNoN), lK(dof*dof,eNoN,eNoN))

!     Loop over all elements of mesh
      DO e=1, lM%nEl
!        Update domain and proceed if domain phys and eqn phys match
         cDmn  = DOMAIN(lM, cEq, e)
         cPhys = eq(cEq)%dmn(cDmn)%phys
         IF (cPhys .NE. phys_heatS) CYCLE

!        Update shape functions for NURBS
         IF (lM%eType .EQ. eType_NRB) CALL NRBNNX(lM, e)

!        Create local copies
         DO a=1, eNoN
            Ac = lM%IEN(a,e)
            ptr(a)  = Ac
            xl(:,a) = x(:,Ac)
            al(:,a) = Ag(:,Ac)
            yl(:,a) = Yg(:,Ac)
         END DO

!        Gauss integration
         lR = 0._RKIND
         lK = 0._RKIND
         DO g=1, lM%nG
            IF (g.EQ.1 .OR. .NOT.lM%lShpF) THEN
               CALL GNN(eNoN, nsd, lM%Nx(:,:,g), xl, Nx, Jac, ksix)
               IF (ISZERO(Jac)) err = "Jac < 0 @ element "//e
            END IF
            w = lM%w(g) * Jac
            N = lM%N(:,g)

            IF (nsd .EQ. 3) THEN
               CALL HEATS3D(eNoN, w, N, Nx, al, yl, lR, lK)

            ELSE IF (nsd .EQ. 2) THEN
               CALL HEATS2D(eNoN, w, N, Nx, al, yl, lR, lK)

            END IF
         END DO ! g: loop

!        Assembly
#ifdef WITH_TRILINOS
         IF (eq(cEq)%assmTLS) THEN
            CALL TRILINOS_DOASSEM(eNoN, ptr, lK, lR)
         ELSE
#endif
            CALL DOASSEM(eNoN, ptr, lK, lR)
#ifdef WITH_TRILINOS
         END IF
#endif
      END DO ! e: loop

      DEALLOCATE(ptr, xl, al, yl, N, Nx, lR, lK)

      RETURN
      END SUBROUTINE CONSTRUCT_HEATS
!####################################################################
!     This is for solving the heat equation in a solid
      PURE SUBROUTINE HEATS3D (eNoN, w, N, Nx, al, yl, lR, lK)
      USE COMMOD
      IMPLICIT NONE
      INTEGER(KIND=IKIND), INTENT(IN) :: eNoN
      REAL(KIND=RKIND), INTENT(IN) :: w, N(eNoN), Nx(nsd,eNoN),
     2   al(tDof,eNoN), yl(tDof,eNoN)
      REAL(KIND=RKIND), INTENT(INOUT) :: lR(1,eNoN), lK(1,eNoN,eNoN)

      INTEGER(KIND=IKIND) i, a, b
      REAL(KIND=RKIND) nu, T1, amd, wl, Td, Tx(nsd), s, rho

      nu  = eq(cEq)%dmn(cDmn)%prop(conductivity)
      s   = eq(cEq)%dmn(cDmn)%prop(source_term)
      rho = eq(cEq)%dmn(cDmn)%prop(solid_density)

      T1  = eq(cEq)%af*eq(cEq)%gam*dt
      amd = eq(cEq)%am * rho/T1
      i   = eq(cEq)%s

      wl = w*T1

      Td = -s
      Tx = 0._RKIND
      DO a=1, eNoN
         Td = Td + N(a)*al(i,a)

         Tx(1) = Tx(1) + Nx(1,a)*yl(i,a)
         Tx(2) = Tx(2) + Nx(2,a)*yl(i,a)
         Tx(3) = Tx(3) + Nx(3,a)*yl(i,a)
      END DO
      Td = Td * rho

      DO a=1,eNoN
         lR(1,a) = lR(1,a) + w*(N(a)*Td
     2      + (Nx(1,a)*Tx(1) + Nx(2,a)*Tx(2) + Nx(3,a)*Tx(3))*nu)

         DO b=1,eNoN
            lK(1,a,b) = lK(1,a,b) + wl*(N(a)*N(b)*amd
     2         + nu*(Nx(1,a)*Nx(1,b) +Nx(2,a)*Nx(2,b) +Nx(3,a)*Nx(3,b)))
         END DO
      END DO

      RETURN
      END SUBROUTINE HEATS3D
!--------------------------------------------------------------------
!     This is for solving the heat equation in a solid
      PURE SUBROUTINE HEATS2D (eNoN, w, N, Nx, al, yl, lR, lK)
      USE COMMOD
      IMPLICIT NONE
      INTEGER(KIND=IKIND), INTENT(IN) :: eNoN
      REAL(KIND=RKIND), INTENT(IN) :: w, N(eNoN), Nx(nsd,eNoN),
     2   al(tDof,eNoN), yl(tDof,eNoN)
      REAL(KIND=RKIND), INTENT(INOUT) :: lR(1,eNoN), lK(1,eNoN,eNoN)

      INTEGER(KIND=IKIND) i, a, b
      REAL(KIND=RKIND) nu, T1, amd, wl, Td, Tx(nsd), s, rho

      nu  = eq(cEq)%dmn(cDmn)%prop(conductivity)
      s   = eq(cEq)%dmn(cDmn)%prop(source_term)
      rho = eq(cEq)%dmn(cDmn)%prop(solid_density)

      T1  = eq(cEq)%af*eq(cEq)%gam*dt
      amd = eq(cEq)%am * rho/T1
      i   = eq(cEq)%s

      wl = w*T1

      Td = -s
      Tx = 0._RKIND
      DO a=1, eNoN
         Td = Td + N(a)*al(i,a)

         Tx(1) = Tx(1) + Nx(1,a)*yl(i,a)
         Tx(2) = Tx(2) + Nx(2,a)*yl(i,a)
      END DO
      Td = Td * rho

      DO a=1, eNoN
         lR(1,a) = lR(1,a) + w*(N(a)*Td + (Nx(1,a)*Tx(1)
     2      + Nx(2,a)*Tx(2))*nu)

         DO b=1, eNoN
            lK(1,a,b) = lK(1,a,b) + wl*(N(a)*N(b)*amd
     2         + nu*(Nx(1,a)*Nx(1,b) + Nx(2,a)*Nx(2,b)))
         END DO
      END DO

      RETURN
      END SUBROUTINE HEATS2D
!--------------------------------------------------------------------
      PURE SUBROUTINE BHEATS (eNoN, w, N, h, lR)
      USE COMMOD
      IMPLICIT NONE
      INTEGER(KIND=IKIND), INTENT(IN) :: eNoN
      REAL(KIND=RKIND), INTENT(IN) :: w, N(eNoN), h
      REAL(KIND=RKIND), INTENT(INOUT) :: lR(dof,eNoN)

      INTEGER(KIND=IKIND) a

      DO a=1, eNoN
         lR(1,a) = lR(1,a) + w*N(a)*h
      END DO

      RETURN
      END SUBROUTINE BHEATS
!####################################################################

max_stars_repo_name: ktbolt/svFSI

max_issues_repo_name: ktbolt/svFSI

max_forks_repo_name: ktbolt/svFSI


program bigmat
use m_matrix, only : lala
   call lala(4000) ! small memory footprint initialization
   call lala( 'a=<1 2 3 4; 5 6 7 8>;quit')
   call lala( 'save(''file1'');quit')
   call lala( 'help;quit')
   call lala( ' ')
end program bigmat