text:
      SUBROUTINE GEO_COND_MAP

c***********************************************************************
c                   subroutine geo_cond_map
c***********************************************************************

c  Program Source:  Naval Ocean Systems Center - Code 542

c  Date:
c     10 May 1996

c  Function:
c     This program provides graphical display of the VLF ground
c     conductivity map.

c     Projections:  azimuthal equidistant
c                   gnomonic
c                   mercator
c                   orthographic
c                   rectangular
c                   stereographic

c     The colors are set as follows:
c        green   3e-2/1e-2
c        blue    3e-3/1e-3
c        purple  3e-5
c        red     1e-5
c        black   3e-4
c        yellow  1e-4

c     Latitudes and longitudes in degrees
c     Sign convention is + for N/E and - for S/W

c  Parameters passed:

c  Parameters returned:

c  Common blocks referenced:
c     geo$data
c     grf$cntr
c     lwpc_lun

c  Functions and subroutines referenced:
c     abs

c     geo_coord
c     grf_cntr_fill
c     grf_color
c     grf_limits
c     open_dat
c     recvr2

c     Requires COND$F.DAT

c  References:

c  Change History:

c     18 Dec 1997   Added stereographic projection.

c     10 Oct 1999   Better test for the antipode inserted.

c*******************!***************************************************

      IMPLICIT NONE

c     LWPC parameters
      include      'lwpc_geo.cmn'
      include      'lwpc_lun.cmn'
      include      'grf_cntr.cmn'

c     Conductivity map data
      integer       mxcoord
      parameter    (mxcoord=4000)

      integer  *  2 cond_lat(mxcoord)
      integer  *  2 cond_lon(mxcoord)
      integer  *  2 ndx
      integer  *  2 nrp

      character*  8 color_list(6)
      character* 40 file_name
      logical       first
      integer       k
      integer       n
      integer       m
      integer       mcode
      integer       mxpass
      integer       nrcoord
      integer       nrpass
      real          crclt
      real          dtr
      real          geo_brng
      real          geo_rng
      real          geo_x
      real          geo_y
      real          lat
      real          latx
      real          lon
      real          lonx
      real          pi
      real          prv_brng
      real          plt_x
      real          prv_x
      real          plt_y
      real          prv_y
      real          rclt
      real          rlng
      real          tst_scale
      real          srclt
cxxcDEBUG
cxx      real          xp(mxcoord)
cxx      real          yp(mxcoord)

      data          first/.true./
      data          file_name/'cond$f.dat'/

      data          color_list/'red','purple','yellow','black',
     &                         'blue','green'/

      data          dtr/.01745329252/
      data          pi/3.14159265359/

      if (first) then

c        Open the data file:
         call OPEN_DAT (lwpcGROUND_lun,file_name,'unformatted')

         first=.false.
      end if

c     Store current color
      call GRF_COLOR ('store')

      nrholes=0

21    read (lwpcGROUND_lun,end=29)
     &      ndx,nrp,
     &     (cond_lon(n),cond_lat(n),n=1,nrp)

      mcode=ndx/1000

      nrcoord=nrp

      SELECT CASE( geo_prjctn )

      CASE( 'A','G','O','S' )

c        Azimuthal, Gnomonic, Orthographic or Stereographic

c        Check if this area is inside of mapped area
         k=0
         n=0
         latx=-999.
         lonx=-999.
         do m=1,nrcoord

            lat=cond_lat(m)/2.
            lon=cond_lon(m)/2.

            if (ABS(lat-latx) .gt. .5 .or.
     &          ABS(lon-lonx) .gt. .5) then

               call GEO_COORD
     &             (lat,lon,
     &              geo_rng,geo_brng,geo_x,geo_y,
     &              plt_x,plt_y)

               k=k+1
               if (geo_rng .lt. geo_rng_max) n=n+1

               latx=lat
               lonx=lon
            end if
         end do
         if (n .gt. 0) then

c           Set up plot resolution
            tst_scale=.0001

            nrpass=0
            mxpass=1
            do while (nrpass .lt. mxpass)

               m=1
               prv_x=99.
               prv_y=99.
               nrcntrp=0
               do while (m .le. nrcoord)

                  lat=cond_lat(m)/2.
                  lon=cond_lon(m)/2.

                  call GEO_COORD
     &                (lat,lon,
     &                 geo_rng,geo_brng,geo_x,geo_y,
     &                 plt_x,plt_y)
cxxcDEBUG
cxx                  write(lwpcLOG_lun,
cxx     &                '(i5,2f7.1,f6.1,f8.1,2f6.1)')
cxx     &                  m,lat,lon,
cxx     &                  geo_brng*57.296,geo_rng*6366.,
cxx     &                  plt_x,plt_y
cxxcDEBUG
cxx                  xp(m)=plt_x
cxx                  yp(m)=plt_y

                  if (geo_rng .gt. geo_rng_max) then

                     if (ABS(geo_center_lat+lat) .le. 0.1
     &                  .and.
     &                  ((geo_center_lon .gt. 0. .and. lon .lt. 0.)
     &                  .or.
     &                   (geo_center_lon .lt. 0. .and. lon .gt. 0.))
     &                  .and.
     &                  (ABS(geo_center_lon-lon) .gt. 179.9
     &                  .and.
     &                   ABS(geo_center_lon-lon) .lt. 181.1)) then
cxxcDEBUG
cxx                        write(lwpcLOG_lun,
cxx     &                      '(i5,2f7.1,f6.1,f8.1,2f6.0,'' ant'')'
cxx     &                        m,lat,lon,
cxx     &                        geo_brng*57.296,geo_rng*6366.,
cxx     &                        geo_center_lat,geo_center_lon

c                       This is the antipode which should give a
c                       bearing of zero
                        if (prv_brng .lt. 1.57 .or.
     &                      prv_brng .gt. 4.71) then

c                          Force bearing to 0
                           geo_brng=0.
                        else

c                          Force bearing to 180
                           geo_brng=pi
                        end if
                     end if

c                    Move the point to the edge of the map
                     geo_rng=geo_rng_max
                     call RECVR2
     &                   (geo_lng0,
     &                    geo_clt0,geo_cclt0,geo_sclt0,
     &                    geo_brng,geo_rng,
     &                    rlng,rclt,crclt,srclt)

                     lat=90.-rclt/dtr
                     lon=rlng/dtr

                     call GEO_COORD
     &                   (lat,lon,
     &                    geo_rng,geo_brng,geo_x,geo_y,
     &                    plt_x,plt_y)
cxxcDEBUG
cxx                     write(lwpcLOG_lun,
cxx     &                   '(i5,2f7.1,f6.1,f8.1,2f6.1,'' mod'')')
cxx     &                     m,lat,lon,
cxx     &                     geo_brng*57.296,geo_rng*6366.,
cxx     &                     plt_x,plt_y
                  end if
                  prv_brng=geo_brng

                  if ((plt_x-prv_x)**2
     &               +(plt_y-prv_y)**2 .ge. tst_scale) then

                     if (nrcntrp .eq. mxcntrp) then

c                       Current plot resolution is too fine;
c                       restart with coarser value
                        tst_scale=tst_scale*5.
                        m=nrcoord
                        nrcntrp=-1
                        mxpass=mxpass+1
                     else

c                       Store this point
                        nrcntrp=nrcntrp+1
                        cntrpx (nrcntrp)=plt_x
                        cntrpy (nrcntrp)=plt_y

                        prv_x=plt_x
                        prv_y=plt_y
                     end if
                  end if
                  m=m+1
               end do
               nrpass=nrpass+1
            end do
cxxcDEBUG
cxx            call GRF_COLOR
cxx     &          ('red')
cxx            call GRF_CURVE
cxx     &          (xp,yp,nrcoord,0.,0.,1.,1.,1,20)

            if (nrcntrp .gt. 0) then

               nrpts(0)=nrcntrp

c              Ensure that all points are inside the boundary
               call GRF_LIMITS
     &             (mxholes,nrholes,nrpts,mxcntrp,cntrpx,cntrpy,
     &              geo_size_x,geo_size_y)

               call GRF_CNTR_FILL
     &             (0.,1.,0.,1.,'100%',color_list(mcode))
            end if
         end if

      CASE( 'M','R' )

c        Mercator or Rectangular

c        Set up plot resolution
         tst_scale=.0001

         nrpass=0
         mxpass=1
         do while (nrpass .lt. mxpass)

            m=1
            prv_x=99.
            prv_y=99.
            nrcntrp=0
            do while (m .lt. nrcoord)

               lat=cond_lat(m)/2.
               lon=cond_lon(m)/2.

               call GEO_COORD
     &             (lat,lon,
     &              geo_rng,geo_brng,geo_x,geo_y,
     &              plt_x,plt_y)

               if (nrpass .eq. 0) then

                  if (lon .gt. geo_corner1_lon) then

c                    This line extends to the left of the map;
c                    set up a second pass
                     mxpass=2
                     plt_x=0.
                  else
     &            if (lon .lt. geo_corner2_lon) then

c                    This line extends to the right of the map;
c                    set up a second pass
                     mxpass=2
                     plt_x=geo_size_x
                  end if
               else

c                 This is the second pass;
c                 exclude what has already been drawn
                  if (lon .ge. geo_corner1_lon) then

                     plt_x=0.
                  else
     &            if (lon .le. geo_corner2_lon) then

                     plt_x=geo_size_x
                  end if
               end if

               if (plt_y .lt.         0.) then

                  plt_y=0.
               else
     &         if (plt_y .gt. geo_size_y) then

                  plt_y=geo_size_y
               end if

               if (plt_x .ge. 0. .and. plt_x .le. geo_size_x .and.
     &             plt_y .ge. 0. .and. plt_y .le. geo_size_y) then

                  if (nrcntrp .eq. 0 .or.
     &               ((plt_x-prv_x)**2
     &               +(plt_y-prv_y)**2 .ge. tst_scale)) then

                     if (nrcntrp .eq. mxcntrp-3) then

c                       Current plot resolution is too fine;
c                       restart with coarser value
                        tst_scale=tst_scale*5.
                        m=nrcoord
                        nrcntrp=-1
                     else

c                       Store this point
                        nrcntrp=nrcntrp+1
                        cntrpx (nrcntrp)=plt_x
                        cntrpy (nrcntrp)=plt_y

                        prv_x=plt_x
                        prv_y=plt_y
                     end if
                  end if
               end if
               m=m+1
            end do

            if (nrcntrp .gt. 0) then

               nrpts(0)=nrcntrp

               call GRF_CNTR_FILL
     &             (0.,1.,0.,1.,'100%',color_list(mcode))
            end if

            if (nrcntrp .ge. 0) then

               nrpass=nrpass+1
               if (nrpass .lt. mxpass) then

c                 Shift this line segment to the left and
c                 do another pass
                  do m=1,nrcoord
                     cond_lon(m)=cond_lon(m)-720
                  end do
               end if
            end if
         end do

      END SELECT
      go to 21
29    REWIND (lwpcGROUND_lun)

c     Reset color
      call GRF_COLOR ('reset')

      RETURN
      END      ! GEO_COND_MAP

max_stars_repo_name: spinorkit/LWPC

max_issues_repo_name: spinorkit/LWPC

max_forks_repo_name: spinorkit/LWPC


    module mod_vectors
    use, intrinsic :: iso_fortran_env, only : sp=>real32, wp => real64, li => int64
    implicit none ( type, external )

    real(wp), parameter :: nan64 =  transfer(-2251799813685248_li, 1._wp)
    real(wp), parameter :: pi = 4d0*atan(1d0), deg = pi/180d0
    reaL(wp), parameter :: o_(3) = [0d0,0d0,0d0]
    reaL(wp), parameter :: i_(3) = [1d0,0d0,0d0]
    reaL(wp), parameter :: j_(3) = [0d0,1d0,0d0]
    reaL(wp), parameter :: k_(3) = [0d0,0d0,1d0]
    real(wp), parameter :: eye3(3,3) = reshape([i_,j_,k_],[3,3])

    interface eye
        procedure a_identity, a_elemental, v_elemental
    end interface
    interface null
        procedure a_null
    end interface

    interface mag
        procedure v_mag
    end interface
    !interface dot
    !    procedure v_dot
    !end interface
    interface cross
        procedure v_cross_op, v_cross
    end interface
    interface operator (.x.)
        procedure v_cross_op, v_cross
    end interface
    interface diag
        procedure v_diag, D_diag
    end interface
    interface inv
        procedure mat3_inv
    end interface
    interface solve
        procedure v_solve, v_solve2
    end interface

    contains

    pure function i_null(n, index) result(non)
    integer, intent(in) :: n, index(:)
    integer, allocatable :: non(:)
    integer i, k, nk
        nk = size(index)
        k = 0
        allocate(non(n-nk))
        do i=1, n
            if ( any(index==i) ) then
            else
                k = k + 1
                non(k) = i
            end if
        end do
    end function

    pure function a_identity(n) result(A)
    integer, intent(in) :: n
    real(wp) :: A(n,n)
    integer :: i
        A = 0d0
        forall(i=1:n)
            A(i,i) = 1d0
        end forall
    end function

    pure function v_elemental(n,i) result(v)
    integer, intent(in) :: n, i
    real(wp) :: v(n)
        v = 0d0
        v(i) = 1d0
    end function

    pure function a_elemental(n,index) result(A)
    integer, intent(in) :: n, index(:)
    real(wp), allocatable :: A(:,:)
    integer :: i, nk
        nk = size(index)
        allocate(A(n, nk), source=0d0)
        do i=1, nk
            A(index(i),i) = 1d0
        end do
    end function

    pure function a_null(n, index) result(A)
    integer, intent(in) :: n, index(:)
    real(wp), allocatable:: A(:,:)
    integer, allocatable :: non(:)
    integer :: i, nk
        nk = size(index)
        non = i_null(n, index)
        allocate(A(n,n-nk), source=0d0)
        do i=1,n-nk
            A(non(i),i) = 1d0
        end do
    end function

    pure function v_cross(v_1,v_2) result(v_3)
    real(wp), intent(in) :: v_1(3), v_2(3)
    real(wp) :: v_3(3)
        v_3 = [ v_1(2)*v_2(3)-v_1(3)*v_2(2), &
                v_1(3)*v_2(1)-v_1(1)*v_2(3), &
                v_1(1)*v_2(2)-v_1(2)*v_2(1) ]
    end function

    pure function v_cross_op(v) result(vx)
    real(wp), intent(in) :: v(3)
    real(wp) :: vx(3,3)
        vx(1,:) = [0d0, -v(3), v(2)]
        vx(2,:) = [v(3), 0d0, -v(1)]
        vx(3,:) = [-v(2), v(1), 0d0]
    end function

    pure function v_mmoi(v) result(mm)
    real(wp), intent(in) :: v(3)
    real(wp) :: mm(3,3)
        mm(1,:) = [v(2)**2+v(3)**2, -v(1)*v(2), -v(1)*v(3)]
        mm(2,:) = [-v(1)*v(2), v(1)**2+v(3)**2, -v(2)*v(3)]
        mm(3,:) = [-v(1)*v(3), -v(2)*v(3), v(1)**2+v(2)**2]
    end function

    pure function v_mag(v) result(m)
    !$omp declare simd(v_mag)
    real(wp), intent(in) :: v(:)
    real(wp) :: m
        m = sqrt( sum( v**2 ))
    end function

    pure function v_dot(v, w) result(m)
    !$omp declare simd(v_dot)
    real(wp), intent(in) :: v(:), w(:)
    real(wp) :: m
        m = sum(v*w)
	end function

    pure function v_diag(v) result(D)
    !$omp declare simd(v_diag)
    real(wp), intent(in) :: v(:)
    real(wp),allocatable :: D(:,:)
    integer :: n, i
        n = size(v)
        allocate(D(n,n))
        D = 0d0
        forall(i=1:n)
            D(i,i) = v(i)
        end forall
	end function

    pure function D_diag(D) result(v)
    !$omp declare simd(D_diag)
    real(wp), intent(in) :: D(:,:)
    real(wp), allocatable :: v(:)
    integer :: n,i
        n = min(size(D,1),size(D,2))
        allocate(v(n))
        do i=1,n
            v(i) = D(i,i)
        end do
    end function

    pure function v_matmul(A, v) result(w)
    ! evaluate w = A*v
    real(wp), intent(in) :: A(:,:), v(:)
	real(wp) :: w(size(A,1))
    integer :: i
        do i=1,size(A,1)
            w(i) = sum(A(i,:)*v(:))
        end do
    end function

    pure function v_matmul2(v, A) result(w)
    ! evaluate w = v*A
    real(wp), intent(in) :: A(:,:), v(:)
	real(wp) :: w(size(A,2))
    integer :: i
        do i=1,size(A,2)
            w(i) = sum(v(:)*A(:,i))
        end do
    end function

    pure function mat3_inv(A) result(B)
    !$omp declare simd(mat3_inv)
    real(wp), intent(in) :: A(3,3)
	real(wp) :: B(3,3), d, t2, t3, t4, t7, t8, t9, t6
        t2 = A(1,1)*A(2,2)*A(3,3)
        t3 = A(1,2)*A(2,3)*A(3,1)
        t4 = A(1,3)*A(2,1)*A(3,2)
        t7 = A(1,1)*A(2,3)*A(3,2)
        t8 = A(1,2)*A(2,1)*A(3,3)
        t9 = A(1,3)*A(2,2)*A(3,1)
        d = t2+t3+t4-t7-t8-t9
        t6 = 1.0D0/d
        B(1,1) = t6*(A(2,2)*A(3,3)-A(2,3)*A(3,2))
        B(1,2) = -t6*(A(1,2)*A(3,3)-A(1,3)*A(3,2))
        B(1,3) = t6*(A(1,2)*A(2,3)-A(1,3)*A(2,2))
        B(2,1) = -t6*(A(2,1)*A(3,3)-A(2,3)*A(3,1))
        B(2,2) = t6*(A(1,1)*A(3,3)-A(1,3)*A(3,1))
        B(2,3) = -t6*(A(1,1)*A(2,3)-A(1,3)*A(2,1))
        B(3,1) = t6*(A(2,1)*A(3,2)-A(2,2)*A(3,1))
        B(3,2) = -t6*(A(1,1)*A(3,2)-A(1,2)*A(3,1))
        B(3,3) = t6*(A(1,1)*A(2,2)-A(1,2)*A(2,1))
    end function

    pure function v_solve(A, b) result(x)
    !$omp declare simd(v_solve)
    ! solve b = A*x for x
    real(wp), intent(in) :: A(3,3), b(3)
	real(wp) :: x(3), d, w(3)

        d = A(1,1)*(A(2,2)*A(3,3)-A(2,3)*A(3,2))+A(1,2)*(A(2,3)*A(3,1)-A(2,1)*A(3,3))+A(1,3)*(A(2,1)*A(3,2)-A(2,2)*A(3,1))
        w = b/d
        x = [ &
            w(1)*(A(2,2)*A(3,3)-A(2,3)*A(3,2))+w(2)*(A(1,3)*A(3,2)-A(1,2)*A(3,3))+w(3)*(A(1,2)*A(2,3)-A(1,3)*A(2,2)), &
            w(1)*(A(2,3)*A(3,1)-A(2,1)*A(3,3))+w(2)*(A(1,1)*A(3,3)-A(1,3)*A(3,1))+w(3)*(A(1,3)*A(2,1)-A(1,1)*A(2,3)), &
            w(1)*(A(2,1)*A(3,2)-A(2,2)*A(3,1))+w(2)*(A(1,2)*A(3,1)-A(1,1)*A(3,2))+w(3)*(A(1,1)*A(2,2)-A(1,2)*A(2,1))]
    end function

    pure function v_solve2(A, B) result(x)
    !$omp declare simd(v_solve2)
    ! solve b = A*x for x
    real(wp), intent(in) :: A(3,3), B(3,3)
	real(wp) :: x(3,3), d, w(3,3)
    integer :: i

        d = A(1,1)*(A(2,2)*A(3,3)-A(2,3)*A(3,2))+A(1,2)*(A(2,3)*A(3,1)-A(2,1)*A(3,3))+A(1,3)*(A(2,1)*A(3,2)-A(2,2)*A(3,1))

        forall (i=1:3)
            w(:,i) = B(:,i)/d
            x(:,i) = [ &
                w(1,i)*(A(2,2)*A(3,3)-A(2,3)*A(3,2))+w(2,i)*(A(1,3)*A(3,2)-A(1,2)*A(3,3))+w(3,i)*(A(1,2)*A(2,3)-A(1,3)*A(2,2)), &
                w(1,i)*(A(2,3)*A(3,1)-A(2,1)*A(3,3))+w(2,i)*(A(1,1)*A(3,3)-A(1,3)*A(3,1))+w(3,i)*(A(1,3)*A(2,1)-A(1,1)*A(2,3)), &
                w(1,i)*(A(2,1)*A(3,2)-A(2,2)*A(3,1))+w(2,i)*(A(1,2)*A(3,1)-A(1,1)*A(3,2))+w(3,i)*(A(1,1)*A(2,2)-A(1,2)*A(2,1))]
        end forall

        !do concurrent (i=1:3)
        !    w = B(:,i)/d
        !    x(:,i) = [ &
        !        w(1)*(A(2,2)*A(3,3)-A(2,3)*A(3,2))+w(2)*(A(1,3)*A(3,2)-A(1,2)*A(3,3))+w(3)*(A(1,2)*A(2,3)-A(1,3)*A(2,2)), &
        !        w(1)*(A(2,3)*A(3,1)-A(2,1)*A(3,3))+w(2)*(A(1,1)*A(3,3)-A(1,3)*A(3,1))+w(3)*(A(1,3)*A(2,1)-A(1,1)*A(2,3)), &
        !        w(1)*(A(2,1)*A(3,2)-A(2,2)*A(3,1))+w(2)*(A(1,2)*A(3,1)-A(1,1)*A(3,2))+w(3)*(A(1,1)*A(2,2)-A(1,2)*A(2,1))]
        !end do
    end function

    function v_rotate_x(v,q) result(u)
    real(wp), intent(in) :: v(3), q
    real(wp) :: u(3), cq, sq
        cq = cos(q)
        sq = sin(q)
        u = [ v(1), &
            cq*v(2)-sq*v(3), &
            sq*v(2)+cq*v(3)]
    end function
    function v_rotate_y(v,q) result(u)
    real(wp), intent(in) :: v(3), q
    real(wp) :: u(3), cq, sq
        cq = cos(q)
        sq = sin(q)
        u = [ cq*v(1)+sq*v(3), &
            v(2), &
            -sq*v(1)+cq*v(3)]
    end function
    function v_rotate_z(v,q) result(u)
    real(wp), intent(in) :: v(3), q
    real(wp) :: u(3), cq, sq
        cq = cos(q)
        sq = sin(q)
        u = [ cq*v(1)-sq*v(1), &
            sq*v(1)+cq*v(2), &
            v(3)]
    end function

    function a_mmul(A,x,b) result(y)
    !tex: Matrix vector calculation
    ! $$ y = A\,x + b $$
    ! or
    ! $$ y = A\,x  $$
    real(wp), intent(in) :: A(:,:), x(:)
    real(wp), intent(in), optional :: b(:)
    real(wp), allocatable :: y(:)
        if( present(b) ) then
            y = matmul(A,x)  + b
        else
            y = matmul(A,x)
        end if
    end function

    end module


text:
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! This file is dedicated to the public domain.  If your jurisdiction requires a        !
! specific license:                                                                    !
!                                                                                !
! Copyright (c) Stephen McDowell, 2017-2019                                            !
! License:      CC0 1.0 Universal                                                      !
! License Text: https://creativecommons.org/publicdomain/zero/1.0/legalcode            !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!---------------------------------------------------------------------------------------
! MODULE: conversions
!> \author Thomas E. Kurtz
!
!> \file
!
!> \brief A fortran module for converting between degrees and radians.
!
!> Conversions may be done using either single or double precision.  This file is mostly
!> an adaptation of the excellent crash-course slides presented in Fall 2010:
!> https://pages.mtu.edu/~shene/COURSES/cs201/NOTES/F90-Subprograms.pdf
!---------------------------------------------------------------------------------------
module conversions
  use iso_c_binding, only : c_float, c_double
  implicit none

  !> The single precision representation of \f$\Pi\f$.
  real(c_float),  parameter :: pi_s  = 3.14159265358979323846
  !> The double precision representation of \f$\Pi\f$.
  real(c_double), parameter :: pi_d  = 3.14159265358979323846264338327950288
  !> The single precision representation of ``180.0``.
  real(c_float),  parameter :: s_180 = 180.0
  !> The double precision representation of ``180.0``.
  real(c_double), parameter :: d_180 = 180.0

contains
  !-------------------------------------------------------------------------------------
  !> \brief
  !>     Convert degrees into radians using single precision floating point numbers.
  !>
  !> \param degrees_s
  !>     The input single precision number representing degrees to be converted.
  !>
  !> \return
  !>     The input ``degrees_s`` converted to radians.  Return type is single precision.
  !-------------------------------------------------------------------------------------
  function degrees_to_radians_s(degrees_s) result(radians_s)
    real(c_float), intent(in) :: degrees_s
    real(c_float)             :: radians_s
    radians_s = degrees_s * pi_s / s_180
  end function degrees_to_radians_s

  !-------------------------------------------------------------------------------------
  !> \brief
  !>     Convert radians into degrees using single precision floating point numbers.
  !>
  !> \param radians_s
  !>     The input single precision number representing radians to be converted.
  !>
  !> \return
  !>     The input ``radians_s`` converted to degrees.  Return type is single precision.
  !-------------------------------------------------------------------------------------
  function radians_to_degrees_s(radians_s) result(degrees_s)
    real(c_float), intent(in) :: radians_s
    real(c_float)             :: degrees_s
    degrees_s = radians_s * s_180 / pi_s
  end function radians_to_degrees_s

  !-------------------------------------------------------------------------------------
  !> \brief
  !>     Convert degrees into radians using double precision floating point numbers.
  !>
  !> \param degrees_d
  !>     The input double precision number representing degrees to be converted.
  !>
  !> \return
  !>     The input ``degrees_d`` converted to radians.  Return type is double precision.
  !-------------------------------------------------------------------------------------
  function degrees_to_radians_d(degrees_d) result(radians_d)
    real(c_double), intent(in) :: degrees_d
    real(c_double)             :: radians_d
    radians_d = degrees_d * pi_d / d_180
  end function degrees_to_radians_d

  !-------------------------------------------------------------------------------------
  !> \brief
  !>     Convert radians into degrees using double precision floating point numbers.
  !>
  !> \param radians_d
  !>     The input double precision number representing radians to be converted.
  !>
  !> \return
  !>     The input ``radians_d`` converted to degrees.  Return type is double precision.
  !-------------------------------------------------------------------------------------
  function radians_to_degrees_d(radians_d) result(degrees_d)
    real(c_double), intent(in) :: radians_d
    real(c_double)             :: degrees_d
    degrees_d = radians_d * d_180 / pi_d
  end function radians_to_degrees_d
end module conversions

max_stars_repo_name: dweindl/exhale

max_issues_repo_name: dweindl/exhale

max_forks_repo_name: dweindl/exhale


text:
!-----------------------------------------------------------------------
!     Subroutine:    chisq_bmin
!     Authors:       S. Lazerson (lazerson@pppl.gov)
!     Date:          05/26/2012
!     Description:   This subroutine targets the difference between
!                    the minimum value of |B| on a surface in Boozer
!                    coordinates
!-----------------------------------------------------------------------
      SUBROUTINE chisq_bmax(target,sigma,niter,iflag)
!-----------------------------------------------------------------------
!     Libraries
!-----------------------------------------------------------------------
      USE stellopt_runtime
      USE stellopt_targets
      USE equil_utils
      USE read_boozer_mod

!-----------------------------------------------------------------------
!     Input/Output Variables
!
!-----------------------------------------------------------------------
      IMPLICIT NONE
      REAL(rprec), INTENT(in)    ::  target(nrad)
      REAL(rprec), INTENT(in)    ::  sigma(nrad)
      INTEGER,     INTENT(in)    ::  niter
      INTEGER,     INTENT(inout) ::  iflag

!-----------------------------------------------------------------------
!     Local Variables
!
!-----------------------------------------------------------------------
      INTEGER ::  u, v, ik, ier, nu_b, nv_b
      REAL(rprec) :: bmax_global
      REAL(rprec), ALLOCATABLE :: xu_booz(:),xv_booz(:), bmax(:)
      REAL(rprec), ALLOCATABLE :: modb_booz(:,:,:)

!----------------------------------------------------------------------
!     BEGIN SUBROUTINE
!----------------------------------------------------------------------
      IF (iflag < 0) RETURN
      ik = count(sigma < bigno)
      IF (iflag == 1) WRITE(iunit_out,'(A,2(2X,I3.3))') 'B_MAX ',ik,4
      IF (iflag == 1) WRITE(iunit_out,'(A)') 'TARGET  SIGMA  VAL  NS'
      IF (niter >= 0) THEN
         ! First transform to s,u,v space
         nu_b = 8 * mboz
         nv_b = 8 * nboz
         ALLOCATE(bmax(nu_b))
         ALLOCATE(modb_booz(nu_b,nv_b,ns_b))
         ALLOCATE(xu_booz(nu_b),xv_booz(nv_b))
         FORALL(u=1:nu_b) xu_booz(u) = REAL(u-1)/REAL(nu_b)
         FORALL(v=1:nv_b) xv_booz(v) = REAL(v-1)/REAL(nv_b)
         CALL mntouv(1,ns_b,mnboz_b,nu_b,nv_b,xu_booz,xv_booz,bmnc_b,&
                     ixm_b,ixn_b/nfp_b,modb_booz,0,1)
         IF (lasym_b) THEN
            CALL mntouv(1,ns_b,mnboz_b,nu_b,nv_b,xu_booz,xv_booz,bmns_b,&
                        ixm_b,ixn_b/nfp_b,modb_booz,1,0)
         END IF
         ! Now evaluate
         DO ik = 1, nsd
            IF (sigma(ik) >= bigno) CYCLE
            !DO u = 1, nu_b
            !   bmax(u) = MAXVAL(modb_booz(u,:,ik))
            !END DO
            bmax_global = MAXVAL(MAXVAL(modb_booz(:,:,ik),2),1)
            mtargets = mtargets + 1
            targets(mtargets) = target(ik)
            sigmas(mtargets)  = sigma(ik)
            vals(mtargets)    = bmax_global
            IF (iflag == 1) WRITE(iunit_out,'(3ES22.12E3,I3.3)') target(ik),sigma(ik),vals(mtargets),ik
         END DO
         ! DEALLOCATE
         DEALLOCATE(modb_booz)
         DEALLOCATE(xu_booz,xv_booz)
         DEALLOCATE(bmax)
      ELSE
         DO ik = 1, nrad
            IF (sigma(ik) >= bigno) CYCLE
               lbooz(ik) = .true.
               mtargets = mtargets + 1
               IF (niter == -2) target_dex(mtargets)=jtarget_bmax
         END DO
      END IF
      RETURN
!----------------------------------------------------------------------
!     END SUBROUTINE
!----------------------------------------------------------------------
      END SUBROUTINE chisq_bmax

max_stars_repo_name: joseluisvelasco/STELLOPT

max_issues_repo_name: joseluisvelasco/STELLOPT

max_forks_repo_name: joseluisvelasco/STELLOPT


text:
program test
    ! defines the C API
    use mnt_cmdlineargparser_capi_mod
    use, intrinsic :: iso_c_binding, only: c_int, c_double, c_char, c_ptr

    implicit none
    type(c_ptr)                                      :: prsr
    integer(c_int)                                   :: ier, nargs, i, n, verbosity, nargs1
    character(len=1), allocatable                    :: args(:)
    character(len=mnt_string_size)                   :: argv_full
    integer(c_int)                                   :: nval, vval
    real(c_double)                                   :: dval
    character(len=1), dimension(mnt_string_size)     :: ival, oval
    character(len=mnt_string_size)                   :: ival_f, oval_f

    nargs = command_argument_count()
    nargs1 = nargs + 1

    ! args must be a contiguous string
    allocate(args(nargs1 * mnt_string_size))
    args(:) = char(0)

    do i = 0, nargs
        ! i = 0 is the executable
        call get_command_argument(i, argv_full)
        ! add termination character, trim...
        call mnt_f2c_string(argv_full, args(i*mnt_string_size + 1:(i+1)*mnt_string_size))
    enddo

    ! define the command line arguments
    ier = mnt_cmdlineargparser_new(prsr)

    ! set
    ier = mnt_cmdlineargparser_setint(prsr, "-n"//char(0), -123, &
                                      "some integer"//char(0))

    ier = mnt_cmdlineargparser_setdouble(prsr, "-d"//char(0), -1.23_c_double, &
                                         "some double"//char(0))

    ier = mnt_cmdlineargparser_setstring(prsr, "-i"//char(0), &
                                         "Hello input!"//char(0), &
                                         "some input string"//char(0))

    ier = mnt_cmdlineargparser_setstring(prsr, "-o"//char(0), &
                                         "Hello output!"//char(0), &
                                         "some output string"//char(0))

    ier = mnt_cmdlineargparser_setbool(prsr, "-v"//char(0), 0, &
                                         "verbose"//char(0))

    ! print help
    ier = mnt_cmdlineargparser_help(prsr)

    ! parse the command line arguments
    ier = mnt_cmdlineargparser_parse(prsr, nargs1, mnt_string_size, args(1))

    ! extract
    ier = mnt_cmdlineargparser_getint(prsr, "-n"//char(0), nval)
    ier = mnt_cmdlineargparser_getstring(prsr, "-i"//char(0), size(ival), ival)
    ier = mnt_cmdlineargparser_getstring(prsr, "-o"//char(0), size(oval), oval)
    ier = mnt_cmdlineargparser_getdouble(prsr, "-d"//char(0), dval)
    ier = mnt_cmdlineargparser_getbool(prsr, "-v"//char(0), verbosity)

    call mnt_c2f_string(ival, ival_f)
    call mnt_c2f_string(oval, oval_f)

    print *, "-n arg is ", nval
    print *, "-i arg is '", trim(ival_f), "'"
    print *, "-o arg is '", trim(oval_f), "'"
    print *, "-d arg is ", dval
    print *, "-v arg is ", verbosity

    ! clean up
    deallocate(args)
    ier = mnt_cmdlineargparser_del(prsr)

end program

max_stars_repo_name: pletzer/mint

max_issues_repo_name: pletzer/mint

max_forks_repo_name: pletzer/mint


text:
*DECK CBRT
      FUNCTION CBRT (X)
C***BEGIN PROLOGUE  CBRT
C***PURPOSE  Compute the cube root.
C***LIBRARY   SLATEC (FNLIB)
C***CATEGORY  C2
C***TYPE      SINGLE PRECISION (CBRT-S, DCBRT-D, CCBRT-C)
C***KEYWORDS  CUBE ROOT, ELEMENTARY FUNCTIONS, FNLIB, ROOTS
C***AUTHOR  Fullerton, W., (LANL)
C***DESCRIPTION
C
C CBRT(X) calculates the cube root of X.
C
C***REFERENCES  (NONE)
C***ROUTINES CALLED  R1MACH, R9PAK, R9UPAK
C***REVISION HISTORY  (YYMMDD)
C   770601  DATE WRITTEN
C   890531  Changed all specific intrinsics to generic.  (WRB)
C   890531  REVISION DATE from Version 3.2
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C***END PROLOGUE  CBRT
      DIMENSION CBRT2(5)
      SAVE CBRT2, NITER
      DATA CBRT2(1) / 0.6299605249 4743658E0 /
      DATA CBRT2(2) / 0.7937005259 8409974E0 /
      DATA CBRT2(3) / 1.0E0 /
      DATA CBRT2(4) / 1.2599210498 9487316E0 /
      DATA CBRT2(5) / 1.5874010519 6819947E0 /
      DATA NITER / 0 /
C***FIRST EXECUTABLE STATEMENT  CBRT
      IF (NITER.EQ.0) NITER = 1.443*LOG(-.106*LOG(0.1*R1MACH(3))) + 1.
C
      CBRT = 0.0
      IF (X.EQ.0.) RETURN
C
      CALL R9UPAK (ABS(X), Y, N)
      IXPNT = N/3
      IREM = N - 3*IXPNT + 3
C
C THE APPROXIMATION BELOW IS A GENERALIZED CHEBYSHEV SERIES CONVERTED
C TO POLYNOMIAL FORM.  THE APPROX IS NEARLY BEST IN THE SENSE OF
C RELATIVE ERROR WITH 4.085 DIGITS ACCURACY.
C
      CBRT = .439581E0 + Y*(.928549E0 + Y*(-.512653E0 + Y*.144586E0))
C
      DO 10 ITER=1,NITER
        CBRTSQ = CBRT*CBRT
        CBRT = CBRT + (Y-CBRT*CBRTSQ)/(3.0*CBRTSQ)
 10   CONTINUE
C
      CBRT = R9PAK (CBRT2(IREM)*SIGN(CBRT,X), IXPNT)
      RETURN
C
      END

ext: f


program area_mandelbrot
  implicit none
!
  integer, parameter :: sp = kind(1.0)
  integer, parameter :: dp = kind(1.0d0)
  integer :: i, j, iter, numoutside
  integer, parameter :: npoints = 4000, maxiter = 4000
  real (kind=dp) :: area, error
  complex (kind=dp) :: c , z
  real :: starttime, endtime
!
! Calculate area of mandelbrot set
!
!    Outer loops runs over npoints, initialize z=c
!
!    Inner loop has the iteration z=z*z+c, and threshold test
!

  numoutside = 0

  call cpu_time(starttime)

  do i = 0,npoints-1
     do j= 0,npoints-1
        c = cmplx(-2.0+(2.5*i)/npoints + 1.0d-07,(1.125*j)/npoints + 1.0d-07)
        z = c
        iter = 0
        do while (iter < maxiter)
           z = z*z + c
           iter = iter + 1
           if (real(z)*real(z)+aimag(z)*aimag(z) > 4) then
              numoutside = numoutside + 1
              exit
           endif
        end do
     end do
  end do

  call cpu_time(endtime)
!
! Output results
!
  area = 2.0*2.5*1.125 * real(npoints*npoints-numoutside)/real(npoints*npoints)
  error = area/real(npoints)
  print *, "Area of Mandelbrot set = ",area," +/- ",error
  print *, "Tome taken for calculation: ",endtime-starttime
!
  stop
end program area_mandelbrot


text:
! { dg-do compile }
! PR 71085
!
! Testcase from PR by Vladimir Fuka <vladimir.fuka@gmail.com>
!
program pr71085
  print *, f()
contains
  function f()
    integer :: f(iargc()*10)
  end
end

max_stars_repo_name: TUDSSL/TICS

max_issues_repo_name: TUDSSL/TICS

max_forks_repo_name: TUDSSL/TICS


text:
Hank Nguyen

Currently attending UC Davis School of Law, class of 2014

Graduated UC Davis in 2010, with Honors
Bachelor in Political Science and Philosophy

Eagle Scout and Merit Badge Counselor

Hank spends way too much time at The Davis Beer Shoppe, Varsity Theater, and the Davis Farmers Market.

Welcome to the wiki I think the term you were looking for was Davisites Davisite ~Daubert

Thanks! ~Hank Nguyen

Oh hey, wassup? ~Users/MichelleWarren

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
20120523 10:30:02 nbsp Do you own the copyrights to the pictures youre using on Executive Cleaners? If not, you cant use them here. Also, I would appreciate it if you stopped deleting photographs of actual work done by the company from the page. Users/WilliamLewis

20120523 10:37:14 nbsp As another note, please understand that you might be doing something on behalf of the owner of Executive Cleaners, but they do not have exclusive ownership of the entry on the wiki about Executive Cleaners.  They dont own it any more than they would own an article in the Enterprise about Executive Cleaners.   Please see Welcome to the Wiki/Business Owner for more tips that might help you understand. Users/JabberWokky Evan JabberWokky Edwards, jw.dw@timewarp.org, (814) 8898845

20120523 10:43:53 nbsp Ok, but if it is a representation of everyone, how come you are the only one making changes, yet mine all get deleted. All i wanted to do was substitute an ugly looking picture of a repair with something else that looks more elegant. Users/ZacharySturges
 What do illegally used stock photos have to do with Executive Cleaners? This is not a glossy flier or a website for the business. Your changes are being deleted because they are inappropriate for the wiki. Users/WilliamLewis

20120523 10:58:09 nbsp We have permission from Selix tuxedos to use the picture, which is from the catalog of their tuxedos, styles which are available for rental, which was taken by their personal photographer. Forget the alterations photo. Users/ZacharySturges
   Do you have permission to release it to the community under Copyrights Creative Commons?  Remember, as it says at the bottom of each page, any content on the wiki can be used by anybody for any purpose, so long as they credit the source. Users/JabberWokky Evan JabberWokky Edwards

20120612 17:08:25 nbsp Please stop deleting the user uploaded photo, showing an actual product from Executive Cleaners.  You may wish to read the welcome to the wiki page (linked above) again. Users/CovertProfessor


text:
	SUBROUTINE FGTBOX  ( ill, jll, box, xx1, yy1, xx2, yy2,
     +                       zgrid, iret )
C************************************************************************
C* FGTBOX								*
C*									*
C* This subroutine gets the corners of the current grid box.  The grid	*
C* itself is stored in common.						*
C*									*
C* FGTBOX  ( ILL, JLL, BOX, XX1, YY1, XX2, YY2, ZGRID, IRET )		*
C*									*
C* Input parameters:							*
C*	ILL		INTEGER		X coordinate			*
C*	JLL		INTEGER		Y coordinate			*
C*      ZGRID (KX*KY)   REAL            Array for grid points           *
C*									*
C* Output parameters:							*
C*	BOX (4)		REAL		Corners of box			*
C*	XX1		REAL		Lower left x coordinate		*
C*	YY1		REAL		Lower left y coordinate		*
C*	XX2		REAL		Upper right x coordinate	*
C*	YY2		REAL		Upper right y coordinate	*
C*	IRET		INTEGER		Return code			*
C*					  0 = normal return		*
C**									*
C* Log:									*
C* M. desJardins/NMC	12/91				 		*
C* S. Gilbert/NCEP      03/07   Added work arrays to call sequence      *
C************************************************************************
	INCLUDE		'GEMPRM.PRM'
	INCLUDE		'CONTUR.CMN'
C*
	REAL		box (4), zgrid (*)
C*
	INCLUDE		'ERMISS.FNC'
C----------------------------------------------------------------------
	iret = 0
C
C*	Get the corners of the box.
C
	box (1) = zgrid ( (jll-1) * isize + ill )
	box (2) = zgrid ( (jll-1) * isize + ill + 1 )
	box (3) = zgrid ( ( jll ) * isize + ill + 1 )
	box (4) = zgrid ( ( jll ) * isize + ill )
C
C*	Check for missing data.
C
	IF  ( ( ERMISS ( box (1) ) ) .or.
     +	      ( ERMISS ( box (2) ) ) .or.
     +	      ( ERMISS ( box (3) ) ) .or.
     +	      ( ERMISS ( box (4) ) ) )  THEN
	    iret = -1
	  ELSE
	    xx1 = ill
	    yy1 = jll
	    xx2 = xx1 + 1.0
	    yy2 = yy1 + 1.0
	END IF
C*
	RETURN
	END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


module ConnectionModule

  use ConstantsModule, only: DZERO

  type :: ConnectionType
    integer :: noden ! node number of cell in parent where ghost node resides
    integer :: nodem ! node number of the connecting cell to which flow occurs from the ghost node
    integer :: kc, ic, jc, kp, ip, jp  ! layer, row, column indices of child and parent nodes
    ! ihc = 0 -- connection is vertical
    ! ihc = 1 -- connection is horizontal
    ! ihc = 2 -- connection is horizontal for vert. staggered grid
    integer :: ihc
    ! cl1: Distance between the center of cell noden and its
    !      shared face with nodem.
    double precision :: cl1
    ! cl2: Distance between the center of cell nodem and its
    !      shared face with noden.
    double precision :: cl2
    double precision :: hwva ! horizontal width or area of vertical connection
    ! for ghost-node correction
    ! NODEJ1 is contributing (adjacent) node in parent grid used for
    ! head interpolation in direction 1 from cell where ghost node resides.
    ! NODEJ2 is contributing (adjacent) node in parent grid used for
    ! head interpolation in direction 2 from cell where ghost node resides.
    integer :: nodej1 = 0
    integer :: k1, i1, j1  ! layer, row, column indices of contributing node nodej1
    integer :: nodej2 = 0
    integer :: k2, i2, j2  ! layer, row, column indices of contributing node nodej2
    integer :: nodej12 = 0 ! for the off-corner cell
    integer :: k12, i12, j12  ! layer, row, column indices of contributing node nodej12
    ! ALPHAJ1 and ALPHAJ2 are contributing factors (portion of head
    ! attributable to contributing node) for contributing cells 1 and 2.
    double precision :: alphaj1 = DZERO
    double precision :: alphaj2 = DZERO
    double precision :: alphaj12 = DZERO ! for the off-corner cell
  contains
    procedure, public :: WriteConnection
    procedure, public :: WriteGhostNodeCorrection
  end type ConnectionType

  contains

  subroutine WriteConnection(this, iu)
    implicit none
    ! dummy
    class(ConnectionType) :: this
    integer, intent(in)   :: iu
    ! format
    20 format(1x,3(1x,i0),1x,3(1x,i0),2x,i0,3(2x,g0))
    !
    write(iu,20)this%kp, this%ip, this%jp, this%kc, this%ic, this%jc, &
                this%ihc, this%cl1, this%cl2, this%hwva
    !
    return
  end subroutine WriteConnection

  subroutine WriteGhostNodeCorrection(this, iu, numalphaj)
    implicit none
    ! dummy
    class(ConnectionType) :: this
    integer, intent(in)   :: iu
    integer, intent(in)   :: numalphaj
    ! format
    10 format(3(1x,3(1x,i0)),2x,g0)
    20 format(4(1x,3(1x,i0)),2(2x,g0))
    30 format(5(1x,3(1x,i0)),3(2x,g0))
    !
    select case (numalphaj)
    case (1)
      write(iu,10)this%kp, this%ip, this%jp, this%kc, this%ic, this%jc, &
                  this%k1, this%i1, this%j1, this%alphaj1
    case (2)
      write(iu,20)this%kp, this%ip, this%jp, this%kc, this%ic, this%jc, &
                  this%k1, this%i1, this%j1, this%k2, this%i2, this%j2, &
                  this%alphaj1, this%alphaj2
    case (3)
      write(iu,30)this%kp, this%ip, this%jp, this%kc, this%ic, this%jc, &
                  this%k1, this%i1, this%j1, this%k2, this%i2, this%j2, &
                  this%k12, this%i12, this%j12, this%alphaj1, this%alphaj2, &
                  this%alphaj12
    end select
    !
    return
  end subroutine WriteGhostNodeCorrection

end module ConnectionModule


PROGRAM TI_Schrod
  USE NumParameters_m
  USE Basis_m
  USE psi_m
  USE op_m
  IMPLICIT NONE

  TYPE (Basis_t), target :: Basis
  TYPE (psi_t)           :: psi,Hpsi
  TYPE (op_t)            :: H

  !====================================================================
  ! read some informations (basis set/grid) : numbers of basis functions, grid points ...
  ! the basis/grid informations have to be put in a module
  CALL Read_Basis(Basis,nio=in_unitp)
  !====================================================================

  write(out_unitp,*) 'Initialization of a real psi'

  CALL init_psi(psi,Basis,cplx=.FALSE.) ! to be changed
  psi%RVec(:) = ONE
  CALL Write_psi(psi)

  write(out_unitp,*) 'Initialization of a complex psi'
  CALL init_psi(psi,Basis,cplx=.TRUE.) ! to be changed
  psi%CVec(:) = CONE
  CALL Write_psi(psi)

  write(out_unitp,*) ' | H | Psi > calculation'
  CALL Set_op(H,Basis) ! to be change
  CALL calc_OpPsi(H,psi,Hpsi)
  CALL Write_psi(Hpsi)

  write(out_unitp,*) 'deallocation'
  CALL dealloc_Op(H)
  CALL dealloc_psi(psi)
  CALL dealloc_psi(Hpsi)

END PROGRAM TI_Schrod


! { dg-do run }
! Checks Fortran 2003's new_line intrinsic function
! PR fortran/28585
program new_line_check
  implicit none
  character(len=10) :: a1
  character(len=10) :: a2(2)
  character(len=10), parameter :: a3 = "1234567890"
  character(len=10), parameter :: a4(2) = "1234567890"
  character(len=10), parameter :: a5(2) = repeat("1234567890",2)

  if(achar(10) /= new_line('a')) STOP 1

  if (iachar(new_line(a1)) /= 10) STOP 2
  if (iachar(new_line(a2)) /= 10) STOP 3
  if (iachar(new_line(a3)) /= 10) STOP 4
  if (iachar(new_line(a4)) /= 10) STOP 5
  if (iachar(new_line(a5)) /= 10) STOP 6

end program new_line_check