text:
C
C Basic image manipulation routines for images
C
C
C
       subroutine img_apply( def_in, def_out, din, op, va, s )
C      -------------------------------------------------------
C
C Apply a specified operation to the supplied (sub-)image
C
C Given:
C   input definition record
       integer    def_in(*)
C   input (sub-)image
       integer    def_out(*)
C   input data
       real*4     din(*)
C   operation to apply
       integer    op
C   control data values
       real*4     va(*)
C
C Updated:
C   status
       integer    s
C
C-

       include 'ic.inc'
       integer  img_index
       integer  i, j, k, t, v, ni, m, it
       real*4   b, xt

       if (s.ne.0) return

       do m=1,len_defn
         image_defn(m) = def_in(m)
       enddo
       b = blank
       do m=1,len_defn
         image_defn(m) = def_out(m)
       enddo
       if (op.eq.1) then
       do t = t1, t2
          do k = z1, z2
             do j = y1, y2
                do i = x1, x2
                   do v = v1,v2
                      ni = img_index( def_in, i, j, k, t, v )
                      if (ni.gt.0) then
                         if (din(ni).gt.0.0) then
                            din(ni) = log(din(ni))
                         else
                            din(ni) = blank
                         endif
                      endif
                   enddo
                enddo
             enddo
          enddo
       enddo

       elseif (op.eq.2) then
       do t = t1, t2
          do k = z1, z2
             do j = y1, y2
                do i = x1, x2
                   do v = v1,v2
                      ni = img_index( def_in, i, j, k, t, v )
                      if (ni.gt.0) then
                         if (din(ni).gt.0.0) then
                            din(ni) = log10(din(ni))
                         else
                            din(ni) = blank
                         endif
                      endif
                   enddo
                enddo
             enddo
          enddo
       enddo

       elseif (op.eq.3) then
       do t = t1, t2
          do k = z1, z2
             do j = y1, y2
                do i = x1, x2
                   do v = v1,v2
                      ni = img_index( def_in, i, j, k, t, v )
                      if (ni.gt.0) then
                         if (din(ni).ne.b) then
                            din(ni) = exp(din(ni))
                         endif
                      endif
                   enddo
                enddo
             enddo
          enddo
       enddo

       elseif (op.eq.5) then
       if (op.eq.4) va(1) = 10.0
       do t = t1, t2
          do k = z1, z2
             do j = y1, y2
                do i = x1, x2
                   do v = v1,v2
                      ni = img_index( def_in, i, j, k, t, v )
                      if (ni.gt.0) then
                         if (din(ni).ne.b) then
                            din(ni) = va(1)**(din(ni))
                         endif
                      endif
                   enddo
                enddo
             enddo
          enddo
       enddo

       elseif (op.eq.6) then
       do t = t1, t2
          do k = z1, z2
             do j = y1, y2
                do i = x1, x2
                   do v = v1,v2
                      ni = img_index( def_in, i, j, k, t, v )
                      if (ni.gt.0) then
                         if (din(ni).ne.b) then
                            din(ni) = din(ni)**va(1)
                         endif
                      endif
                   enddo
                enddo
             enddo
          enddo
       enddo

       elseif (op.eq.7) then
       do t = t1, t2
          do k = z1, z2
             do j = y1, y2
                do i = x1, x2
                   do v = v1,v2
                      ni = img_index( def_in, i, j, k, t, v )
                      if (ni.gt.0) then
                         if (din(ni).ne.b) then
                            din(ni) = sin(din(ni))
                         endif
                      endif
                   enddo
                enddo
             enddo
          enddo
       enddo

       elseif (op.eq.8) then
       do t = t1, t2
          do k = z1, z2
             do j = y1, y2
                do i = x1, x2
                   do v = v1,v2
                      ni = img_index( def_in, i, j, k, t, v )
                      if (ni.gt.0) then
                         if (din(ni).ne.b) then
                            din(ni) = cos(din(ni))
                         endif
                      endif
                   enddo
                enddo
             enddo
          enddo
       enddo

       elseif (op.eq.9) then
       do t = t1, t2
          do k = z1, z2
             do j = y1, y2
                do i = x1, x2
                   do v = v1,v2
                      ni = img_index( def_in, i, j, k, t, v )
                      if (ni.gt.0) then
                         if (din(ni).ne.b) then
                            din(ni) = tan(din(ni))
                         endif
                      endif
                   enddo
                enddo
             enddo
          enddo
       enddo

       elseif (op.eq.10) then
       do t = t1, t2
          do k = z1, z2
             do j = y1, y2
                do i = x1, x2
                   do v = v1,v2
                      ni = img_index( def_in, i, j, k, t, v )
                      if (ni.gt.0) then
                         if (din(ni).ne.b) then
                            din(ni) = asin(din(ni))
                         endif
                      endif
                   enddo
                enddo
             enddo
          enddo
       enddo

       elseif (op.eq.11) then
       do t = t1, t2
          do k = z1, z2
             do j = y1, y2
                do i = x1, x2
                   do v = v1,v2
                      ni = img_index( def_in, i, j, k, t, v )
                      if (ni.gt.0) then
                         if (din(ni).ne.b) then
                            din(ni) = acos(din(ni))
                         endif
                      endif
                   enddo
                enddo
             enddo
          enddo
       enddo

       elseif (op.eq.12) then
       do t = t1, t2
          do k = z1, z2
             do j = y1, y2
                do i = x1, x2
                   do v = v1,v2
                      ni = img_index( def_in, i, j, k, t, v )
                      if (ni.gt.0) then
                         if (din(ni).ne.b) then
                            din(ni) = atan(din(ni))
                         endif
                      endif
                   enddo
                enddo
             enddo
          enddo
       enddo

       elseif (op.eq.13) then
       do t = t1, t2
          do k = z1, z2
             do j = y1, y2
                do i = x1, x2
                   do v = v1,v2
                      ni = img_index( def_in, i, j, k, t, v )
                      if (ni.gt.0) then
                         if (din(ni).ne.b) then
                            din(ni) = din(ni) * va(1)
                         endif
                      endif
                   enddo
                enddo
             enddo
          enddo
       enddo

       elseif (op.eq.14) then
       do t = t1, t2
          do k = z1, z2
             do j = y1, y2
                do i = x1, x2
                   do v = v1,v2
                      ni = img_index( def_in, i, j, k, t, v )
                      if (ni.gt.0) then
                         if (din(ni).ne.b) then
                            din(ni) = din(ni) + va(1)
                         endif
                      endif
                   enddo
                enddo
             enddo
          enddo
       enddo

       elseif (op.eq.15) then
       do t = t1, t2
          do k = z1, z2
             do j = y1, y2
                do i = x1, x2
                   do v = v1,v2
                      ni = img_index( def_in, i, j, k, t, v )
                      if (ni.gt.0) then
                         if (din(ni).ne.b) then
                            if (din(ni).lt.va(1)) then
                               din(ni) = 0.0
                            else
                               din(ni) = 1.0
                            endif
                         endif
                      endif
                   enddo
                enddo
             enddo
          enddo
       enddo

       elseif (op.eq.16) then
       do t = t1, t2
          do k = z1, z2
             do j = y1, y2
                do i = x1, x2
                   do v = v1,v2
                      ni = img_index( def_in, i, j, k, t, v )
                      if (ni.gt.0) then
                         if (din(ni).ne.b) then
                            din(ni) = va(1)
                         endif
                      endif
                   enddo
                enddo
             enddo
          enddo
       enddo

       elseif (op.eq.17) then
       do t = t1, t2
          do k = z1, z2
             do j = y1, y2
                do i = x1, x2
                   do v = v1,v2
                      ni = img_index( def_in, i, j, k, t, v )
                      if (ni.gt.0) then
                         if (din(ni).ne.b) then
                            if (din(ni).lt.va(1) .or.
     *                          din(ni).gt.va(2) ) then
                               din(ni) = va(3)
                            endif
                         endif
                      endif
                   enddo
                enddo
             enddo
          enddo
       enddo

       elseif (op.eq.18) then
       do t = t1, t2
          do k = z1, z2
             do j = y1, y2
                do i = x1, x2
                   do v = v1,v2
                      ni = img_index( def_in, i, j, k, t, v )
                      if (ni.gt.0) then
                         din(ni) = blank
                      endif
                   enddo
                enddo
             enddo
          enddo
       enddo

       elseif (op.eq.19) then
       do t = t1, t2
          do k = z1, z2
             do j = y1, y2
                do i = x1, x2
                   do v = v1,v2
                      ni = img_index( def_in, i, j, k, t, v )
                      if (ni.gt.0) then
                         if (din(ni).eq.b) then
                            din(ni) = va(1)
                         endif
                      endif
                   enddo
                enddo
             enddo
          enddo
       enddo

       elseif (op.eq.20) then
       do t = t1, t2
          do k = z1, z2
             do j = y1, y2
                do i = x1, x2
                   do v = v1,v2
                      ni = img_index( def_in, i, j, k, t, v )
                      if (ni.gt.0) then
                         if (din(ni).ne.b .and.din(ni).ne.0.0) then
                            din(ni) = 1.0/din(ni)
                         endif
                      endif
                   enddo
                enddo
             enddo
          enddo
       enddo

       elseif (op.eq.21) then
       do t = t1, t2
          do k = z1, z2
             do j = y1, y2
                do i = x1, x2
                   do v = v1,v2
                      ni = img_index( def_in, i, j, k, t, v )
                      if (ni.gt.0) then
                         if (din(ni).ne.b .and.din(ni).ge.0.0) then
                            din(ni) = sqrt(din(ni))
                         endif
                      endif
                   enddo
                enddo
             enddo
          enddo
       enddo

       elseif (op.eq.22) then
       do t = t1, t2
          do k = z1, z2
             do j = y1, y2
                do i = x1, x2
                   do v = v1,v2
                      ni = img_index( def_in, i, j, k, t, v )
                      if (ni.gt.0) then
                         if (din(ni).ne.b) then
                            xt = din(ni)
                            if (xt.gt.va(2)) xt = va(2)
                            it = (xt - va(1))/va(3)
                            din(ni) = va(1) + it*va(3)
                         endif
                      endif
                   enddo
                enddo
             enddo
          enddo
       enddo

       elseif (op.eq.30) then
       do t = t1, t2
          do k = z1, z2
             do j = y1, y2
                do i = x1, x2
                   do v = v1,v2
                      ni = img_index( def_in, i, j, k, t, v )
                      if (ni.gt.0) then
                         if (din(ni).ne.b) then
                            if (din(ni).ge.va(1) .and.
     *                          din(ni).le.va(2) ) then
                               din(ni) = va(3)
                            endif
                         endif
                      endif
                   enddo
                enddo
             enddo
          enddo
       enddo
       elseif (op.eq.31) then
       do t = t1, t2
          do k = z1, z2
             do j = y1, y2
                do i = x1, x2
                   do v = v1,v2
                      ni = img_index( def_in, i, j, k, t, v )
                      if (ni.gt.0) then
                         if (din(ni).ne.b) then
                            if (din(ni).ge.va(1) .and.
     *                          din(ni).le.va(2) ) then
                               din(ni) = din(ni) + va(3)
                            endif
                         endif
                      endif
                   enddo
                enddo
             enddo
          enddo
       enddo
       elseif (op.eq.32) then
       do t = t1, t2
          do k = z1, z2
             do j = y1, y2
                do i = x1, x2
                   do v = v1,v2
                      ni = img_index( def_in, i, j, k, t, v )
                      if (ni.gt.0) then
                         if (din(ni).ne.b) then
                            if (din(ni).ge.va(1) .and.
     *                          din(ni).le.va(2) ) then
                               din(ni) = din(ni) - va(3)
                            endif
                         endif
                      endif
                   enddo
                enddo
             enddo
          enddo
       enddo
       elseif (op.eq.33) then
       do t = t1, t2
          do k = z1, z2
             do j = y1, y2
                do i = x1, x2
                   do v = v1,v2
                      ni = img_index( def_in, i, j, k, t, v )
                      if (ni.gt.0) then
                         if (din(ni).ne.b) then
                            if (din(ni).ge.va(1) .and.
     *                          din(ni).le.va(2) ) then
                               din(ni) = din(ni) * va(3)
                            endif
                         endif
                      endif
                   enddo
                enddo
             enddo
          enddo
       enddo
       elseif (op.eq.34) then
       do t = t1, t2
          do k = z1, z2
             do j = y1, y2
                do i = x1, x2
                   do v = v1,v2
                      ni = img_index( def_in, i, j, k, t, v )
                      if (ni.gt.0) then
                         if (din(ni).ne.b) then
                            if (din(ni).ge.va(1) .and.
     *                          din(ni).le.va(2) ) then
                               din(ni) = din(ni) / va(3)
                            endif
                         endif
                      endif
                   enddo
                enddo
             enddo
          enddo
       enddo

       endif
       end
C
C
       subroutine img_combine( def1, def_out, d1, op, def2, d2, s )
C      ------------------------------------------------------------
C
C Combine the intersection of the images def1/def_out and def2
C
C Given:
C   definition record image 1
       integer    def1(*)
C   subimage definition record
       integer    def_out(*)
C
C Updated:
C   input data image 1 (and output)
       real*4     d1(*)
C
C Given:
C   operator code
       integer    op
C
C   definition record image 2
       integer    def2(*)
C   input data image 2
       real*4     d2(*)
C
C Updated:
C   status
       integer    s
C
C-

       include 'ic.inc'
       integer  img_index
       integer  i, j, k, t, v, n1, n2, m
       real*4   b1, b2

       if (s.ne.0) return

       do m=1,len_defn
         image_defn(m) = def1(m)
       enddo
       b1 = blank
       do m=1,len_defn
         image_defn(m) = def2(m)
       enddo
       b2 = blank
       do m=1,len_defn
         image_defn(m) = def_out(m)
       enddo
       n2 = 0
       if (op.eq.1) then
       do t = t1, t2
          do k = z1, z2
             do j = y1, y2
                do i = x1, x2
                   do v = v1,v2
                      n1 = img_index( def1, i, j, k, t, v )
                      n2 = img_index( def2, i, j, k, t, v )
                      if (n1.gt.0 .and. n2.gt.0) then
                         d1(n1) = d2(n2)
                      endif
                   enddo
                enddo
             enddo
          enddo
       enddo

       elseif (op.eq.2) then
       do t = t1, t2
          do k = z1, z2
             do j = y1, y2
                do i = x1, x2
                   do v = v1,v2
                      n1 = img_index( def1, i, j, k, t, v )
                      n2 = img_index( def2, i, j, k, t, v )
                      if (n1.gt.0 .and. n2.gt.0) then
                         if (d1(n1).ne.b1  .and. d2(n2).ne.b2) then
                            d1(n1) = d1(n1) + d2(n2)
                         endif
                      endif
                   enddo
                enddo
             enddo
          enddo
       enddo

       elseif (op.eq.3) then
       do t = t1, t2
          do k = z1, z2
             do j = y1, y2
                do i = x1, x2
                   do v = v1,v2
                      n1 = img_index( def1, i, j, k, t, v )
                      n2 = img_index( def2, i, j, k, t, v )
                      if (n1.gt.0 .and. n2.gt.0) then
                         if (d1(n1).ne.b1  .and. d2(n2).ne.b2) then
                            d1(n1) = d1(n1) - d2(n2)
                         endif
                      endif
                   enddo
                enddo
             enddo
          enddo
       enddo

       elseif (op.eq.4) then
       do t = t1, t2
          do k = z1, z2
             do j = y1, y2
                do i = x1, x2
                   do v = v1,v2
                      n1 = img_index( def1, i, j, k, t, v )
                      n2 = img_index( def2, i, j, k, t, v )
                      if (n1.gt.0 .and. n2.gt.0) then
                         if (d1(n1).ne.b1  .and. d2(n2).ne.b2) then
                            d1(n1) = d1(n1) * d2(n2)
                         endif
                      endif
                   enddo
                enddo
             enddo
          enddo
       enddo

       elseif (op.eq.5) then
       do t = t1, t2
          do k = z1, z2
             do j = y1, y2
                do i = x1, x2
                   do v = v1,v2
                      n1 = img_index( def1, i, j, k, t, v )
                      n2 = img_index( def2, i, j, k, t, v )
                      if (n1.gt.0 .and. n2.gt.0) then
                         if (d1(n1).ne.b1  .and. d2(n2).ne.b2) then
                            if (d2(n2).ne.0.0) then
                              d1(n1) = d1(n1) / d2(n2)
                            else
                              d1(n1) = blank
                            endif
                         endif
                      endif
                   enddo
                enddo
             enddo
          enddo
       enddo

       else if (op.eq.-1) then
       do t = t1, t2
          do k = z1, z2
             do j = y1, y2
                do i = x1, x2
                   do v = v1,v2
                      n1 = img_index( def1, i, j, k, t, v )
                      n2 = n2 + 1
                      if (n1.gt.0 .and. n2.le.def2(2)) then
                         d1(n1) = d2(n2)
                      endif
                   enddo
                enddo
             enddo
          enddo
       enddo

       elseif (op.eq.-2) then
       do t = t1, t2
          do k = z1, z2
             do j = y1, y2
                do i = x1, x2
                   do v = v1,v2
                      n1 = img_index( def1, i, j, k, t, v )
                      n2 = n2 + 1
                      if (n1.gt.0 .and. n2.le.def2(2)) then
                         if (d1(n1).ne.b1  .and. d2(n2).ne.b2) then
                            d1(n1) = d1(n1) + d2(n2)
                         endif
                      endif
                   enddo
                enddo
             enddo
          enddo
       enddo

       elseif (op.eq.-3) then
       do t = t1, t2
          do k = z1, z2
             do j = y1, y2
                do i = x1, x2
                   do v = v1,v2
                      n1 = img_index( def1, i, j, k, t, v )
                      n2 = n2 + 1
                      if (n1.gt.0 .and. n2.le.def2(2)) then
                         if (d1(n1).ne.b1  .and. d2(n2).ne.b2) then
                            d1(n1) = d1(n1) - d2(n2)
                         endif
                      endif
                   enddo
                enddo
             enddo
          enddo
       enddo

       elseif (op.eq.-4) then
       do t = t1, t2
          do k = z1, z2
             do j = y1, y2
                do i = x1, x2
                   do v = v1,v2
                      n1 = img_index( def1, i, j, k, t, v )
                      n2 = n2 + 1
                      if (n1.gt.0 .and. n2.le.def2(2)) then
                         if (d1(n1).ne.b1  .and. d2(n2).ne.b2) then
                            d1(n1) = d1(n1) * d2(n2)
                         endif
                      endif
                   enddo
                enddo
             enddo
          enddo
       enddo

       elseif (op.eq.-5) then
       do t = t1, t2
          do k = z1, z2
             do j = y1, y2
                do i = x1, x2
                   do v = v1,v2
                      n1 = img_index( def1, i, j, k, t, v )
                      n2 = n2 + 1
                      if (n1.gt.0 .and. n2.le.def2(2)) then
                         if (d1(n1).ne.b1  .and. d2(n2).ne.b2) then
                            if (d2(n2).ne.0.0) then
                              d1(n1) = d1(n1) / d2(n2)
                            else
                              d1(n1) = blank
                            endif
                         endif
                      endif
                   enddo
                enddo
             enddo
          enddo
       enddo

       endif
       end
C
C
       subroutine img_combine2( def1, def_out, d1, op, def2, d2,
     *                          def3, d3, s )
C      ---------------------------------------------------------
C
C Combine the intersection of the images def1/def_out and def2
C
C Given:
C   definition record image 1
       integer    def1(*)
C   subimage definition record
       integer    def_out(*)
C
C Updated:
C   input data image 1 (and output)
       real*4     d1(*)
C
C Given:
C   operator code
       integer    op
C
C   definition record image 2
       integer    def2(*)
C   input data image 2
       real*4     d2(*)
C
C   definition record image 3
       integer    def3(*)
C   input data image 3
       real*4     d3(*)
C
C Updated:
C   status
       integer    s
C
C-

       include 'ic.inc'
       integer  img_index
       integer  i, j, k, t, v, n1, n2, n3, m
       real*4   b1, b2, b3

       if (s.ne.0) return

       do m=1,len_defn
         image_defn(m) = def1(m)
       enddo
       b1 = blank
       do m=1,len_defn
         image_defn(m) = def2(m)
       enddo
       b2 = blank
       do m=1,len_defn
         image_defn(m) = def3(m)
       enddo
       b3 = blank
       do m=1,len_defn
         image_defn(m) = def_out(m)
       enddo
       n2 = 0
       if (op.eq.1) then
       do t = t1, t2
          do k = z1, z2
             do j = y1, y2
                do i = x1, x2
                   do v = v1,v2
                      n1 = img_index( def1, i, j, k, t, v )
                      n2 = img_index( def2, i, j, k, t, v )
                      n3 = img_index( def3, i, j, k, t, v )
                      if (n1.gt.0 .and. n2.gt.0 .and. n3.gt.0) then
                        if (d2(n2).ne.b2 .and. d3(n3).ne.b3) then
                           d1(n1) = atan2( d2(n2), d3(n3) )
                        else
                           d1(n1) = blank
                        endif
                      else
                         d1(n1) = blank
                      endif
                   enddo
                enddo
             enddo
          enddo
       enddo
       endif
       end

max_stars_repo_name: CavendishAstrophysics/anmap

max_issues_repo_name: CavendishAstrophysics/anmap

max_forks_repo_name: CavendishAstrophysics/anmap


text:
c$pragma aux Sys_color "*_" parm (value)

      SUBROUTINE GRF_COLOR
     &          (clr_rgb)

c***********************************************************************

c  Function:
c     Determines the color index and then calls the appropriate device
c     color routine

c  Parameters passed:
c     clr_rgb       [s] color name or RGB as "RRR,GGG,BBB"

c  Parameters returned:

c  Common blocks referenced:
c     graphics

c  Functions and subroutines referenced:
c     index

c     decode_lst_int
c     str_length
c     str_lower
c     str_trim_char

c     hpgl_color
c     sys_color

c  Entry:

c  References:

c  Change History:

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

      include      'graphics.cmn'

      character*(*) clr_rgb

      character* 12 xclr_rgb
      character* 12 oldColor(2)
      character*200 error_msg
      logical       match
      integer       ColorIndex,rgb(3),
     &              oldColorNdx

      integer       sys_color,
     &              str_length

c     Default color
      data          oldColorNdx/-9/

c     Determine if color name or rgb value specified
      xclr_rgb=clr_rgb
      nl=STR_LENGTH(xclr_rgb)
      if (nl .eq. 0) then

c        End of graph
         ColorIndex=0
      else

c        Trim leading spaces from the color name
         call STR_TRIM_CHAR (xclr_rgb,' ',nl,'Leading_Blanks')

c        Convert to lower case
         call STR_LOWER (xclr_rgb,1,nl)

         if (xclr_rgb .eq. 'store') then

c           Store current color and exit
            oldColorNdx=ColorIndex
            oldColor(1)=graphics_Color(1)
            oldColor(2)=graphics_Color(2)
            RETURN
         else
     &   if (xclr_rgb .eq. 'reset') then

            if (oldColorNdx .eq. -9) then

c              No stored color; use first one in the list
               ColorIndex       =1
               graphics_Color(1)=graphics_ColorList(1,1)
               graphics_Color(2)=graphics_ColorList(2,1)
            else

c              Restore previous color
               ColorIndex       =oldColorNdx
               graphics_Color(1)=oldColor(1)
               graphics_Color(2)=oldColor(2)
            end if
         else
     &   if (INDEX(xclr_rgb,',') .eq. 0) then

c           Name of color specified

c           If color is same as current, do nothing.
            nl=STR_LENGTH (graphics_Color(1))
            if (xclr_rgb(:nl) .eq. graphics_Color(1)(:nl)) RETURN

c           Set requested color
            nc=1
            match=.false.
            do while (nc .le. graphics_NrColors .and. .not.match)
               nl=STR_LENGTH (graphics_ColorList(1,nc))
               if (xclr_rgb(:nl) .eq.
     &            graphics_ColorList(1,nc)(:nl)) then

                  graphics_Color(1)=graphics_ColorList(1,nc)
                  graphics_Color(2)=graphics_ColorList(2,nc)
                  if (STR_LENGTH(graphics_device) .eq. 3 .and.
     &                graphics_device(:3) .eq. 'sys') then

c                    Output to system device;
c                    set index in the color list
                     ColorIndex=nc
                  else
     &            if (graphics_device(:4) .eq. 'sys-' .or.
     &                graphics_device(:4) .eq. 'sys_') then

c                    Output to system device;
c                    set index in the color list
                     ColorIndex=nc
                  else

c                    Output is to a file;
c                    set equivalent pen number
                     ColorIndex=graphics_ColorPen(nc)
                  end if
                  match=.true.
               else
                  nc=nc+1
               end if
            end do
            if (.not.match) then

               call STR_UPPER (clr_rgb,0,0)
               write(error_msg,
     &             '(''[GRF_COLOR]: '',
     &               ''Requested color ('',a,'') '',
     &               ''not found in the color list'')')
     &                 clr_rgb(:STR_LENGTH(clr_rgb))
               call LWPC_ERROR ('WARNING',error_msg)

               ColorIndex=1
               graphics_Color(1)=graphics_ColorList(1,1)
               graphics_Color(2)=graphics_ColorList(2,1)
            end if
         else

c           RGB specified

c           If color is same as current, do nothing.
            nl=STR_LENGTH (graphics_Color(2))
            if (xclr_rgb(:nl) .eq. graphics_Color(2)(:nl)) RETURN

c           Set requested color
            graphics_Color(1)='rgb'
            graphics_Color(2)=xclr_rgb
            if (graphics_device(:3) .ne. 'sys') then

c              Set first pen; the RGB color will be ignored
               ColorIndex=1
            else

c              Set index to the last color/pen
               ColorIndex=graphics_NrColors
            end if
         end if
      end if

c     Change color
      if (graphics_device(:3) .eq. 'sys') then

         call DECODE_LIST_INT (graphics_Color(2),3,nl,rgb)
         iresult=SYS_COLOR (rgb(1),rgb(2),rgb(3))
      else

         call HPGL_COLOR (ColorIndex)
      end if

      RETURN
      END      ! GRF_COLOR

ext: for

max_stars_repo_name: spinorkit/LWPC

max_issues_repo_name: spinorkit/LWPC

max_forks_repo_name: spinorkit/LWPC


text:
Jump to Timeline #Navigation Navigation

In 1983, the Davis originating band True West released their selftitled first EP, The Davis Shoppers Guide started appearing in peoples mailboxes, and Davis High School grad Denise Curry won the silver medal for basketball in the Pan American games.

 May

   4^th^: Thong H. Huynh stabbed to death at Davis Senior High School  prompts increased awareness of racial tension

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
	program ftr
	implicit doubleprecision (a-h,o-z)
	dimension w(20000),t(20000)
	complex*16 ac(20000),ft(20000),ci,te
	real*4 a1,a2,a3,a4,a5
	character*32 fil1,fil2,fil3
	ci=(0.d0,1.d0)
	open(20,file='ftr_in',status='old')
	read(20,*) nm,beta
	read(20,*) fil1,fil2
        open(25,file=fil1,status='old')
        open(26,file=fil2)
c       energy  units
	eu=27.2d0/1836.15
	emx=50d0
	write(*,*) 'energy unit',eu

	do 1 i=1,nm
	  read(25,1000) x1,ac(i),x2
	 t(i)=x1
1	continue
	ac(1)=ac(1)/2.d0
	ac(nm)=ac(nm)/2.d0
c-------- damp the signal ----------------
        beta=-log(beta/x2)/t(nm)**2
        do 21 i=1,nm
21       ac(i)=ac(i)*exp(-beta*t(i)**2)

	pi=4.d0*atan(1.d0)
	dt=abs(t(2)-t(1))
	write(*,*) dt,nm
	h=0.5d0/nm/dt*pi
	nw=emx/h
	write(*,*) h,nw

	do 3 i=1,nw
	  w(i)=h*i
3	continue

	do 2 k=1,nw
	  ft(k)=(0.d0,0.d0)
	  do 4 i=1,nm
	    te=exp(ci*w(k)*t(i))
4	    ft(k)=ac(i)*te+ft(k)
c	  fin=dreal(ft(k))*dt
	  fin=abs(ft(k))*dt
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
          write(26,1000) w(k),fin,ft(k)*dt
2	  continue

1001    format(20(e14.5,1x))
1000    format(20(e14.7,1x))
	stop
	end


text:
subroutine tem_coupling_outerr(a)
   use Mod_Temperature
   implicit none
   class(TemperatureProblem) :: a

   !Check Advection from outerr Physical problem has not been set
   if(a%kfl_advec==1.and.(a%kfl_ExternalAdvection==0)) &
      call runend('EXTERNAL CONVECTION IN TEMPER IS IMPOSSIBLE IF EXTERNAL VELOCITY IS NOT SET')

end subroutine

max_stars_repo_name: ciaid-colombia/InsFEM

max_issues_repo_name: ciaid-colombia/InsFEM

max_forks_repo_name: ciaid-colombia/InsFEM


text:
************************************************************************
*
*     Subroutine INPUTH4                               Called by: INIT3
*
*     input UCFLAG, ULFLAG, CNAME
*
************************************************************************
      SUBROUTINE INPUTH4(CFACTOR,CNAME,CUNITSTR,LFACTOR,LUNITSTR,NCONST)
*
*     logical devices
*
      INCLUDE 'lda.inc'
*
*     subroutine arguments
*
      INTEGER*4 NCONST
      DOUBLE PRECISION CFACTOR(*),LFACTOR(*)
      CHARACTER*4 CUNITSTR(*)
      CHARACTER*9 LUNITSTR(*)
      CHARACTER*45 CNAME(*)
*
*     local vars
*
      INTEGER*4 J,UCFLAG,ULFLAG
      CHARACTER*8 CUNIT(2) / 'mg/L (1)','ug/L (2)' /
      CHARACTER*9 LUNIT(4) /'kg/d (1)','g/d (2)','lb/d (3)','ton/d (4)'/
      CHARACTER*500 BUFFER
*
*     input format statements
*
 1000 FORMAT(A45,I5,I5)
*
*     output format statements
*
 2000 FORMAT(/,T47,'Conc.',T64,'Load',/,1X,'Constituent',34X,
     &       'Units (Flag)',T64,'Units (Flag)',/,1X,74('-'))
 2100 FORMAT(1X,A45,A8,T64,A9)
*
*     Read name and unit flags for each constituent.  For valid unit
*     flags, set conversion factors (CFACTOR and LFACTOR) and unit
*     strings (CUNITSTR and LUNITSTR).  Call error routines if unit
*     flags are invalid.
*
      WRITE(LDECHO,2000)
      DO 10 J=1,NCONST
         CALL GETLINE(LDHEAD,BUFFER)
         READ (BUFFER,1000) CNAME(J),UCFLAG,ULFLAG
         WRITE (LDECHO,2100) CNAME(J),CUNIT(UCFLAG),LUNIT(ULFLAG)
         IF (UCFLAG .EQ. 1) THEN
            CFACTOR(J) = 86400.D-3 * 0.3048D0 * 0.3048D0 * 0.3048D0
            CUNITSTR(J) = 'MG/L'
         ELSEIF (UCFLAG .EQ. 2) THEN
            CFACTOR(J) = 86400.D-6 * 0.3048D0 * 0.3048D0 * 0.3048D0
            CUNITSTR(J) = 'UG/L'
         ELSE
            CALL ERROR1(6,UCFLAG)
         ENDIF
         IF (ULFLAG .EQ. 1) THEN
            LFACTOR(J) = 1.D0
            LUNITSTR(J) =  '[KG/DAY]'
         ELSEIF (ULFLAG .EQ. 2) THEN
            LFACTOR(J) = 1000.D0
            LUNITSTR(J) = '[G/DAY]'
         ELSEIF (ULFLAG .EQ. 3) THEN
            LFACTOR(J) = 2.2046D0
            LUNITSTR(J) = '[LBS/DAY]'
         ELSEIF (ULFLAG .EQ. 4) THEN
            LFACTOR(J) = 2.2046D0/2000.D0
            LUNITSTR(J) = '[TON/DAY]'
         ELSE
            CALL ERROR1(7,ULFLAG)
         ENDIF
 10   CONTINUE

      RETURN
      END

max_stars_repo_name: lthiamodelers/baseflow-coefficients

max_issues_repo_name: lthiamodelers/baseflow-coefficients

max_forks_repo_name: lthiamodelers/baseflow-coefficients


text:
      subroutine hruyr

!!    ~ ~ ~ PURPOSE ~ ~ ~
!!    this subroutine writes annual HRU output to the output.hru file

!!    ~ ~ ~ INCOMING VARIABLES ~ ~ ~
!!    name          |units         |definition
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
!!    bio_yrms(:)   |metric tons/ha|annual biomass (dry weight) in the HRU
!!    cpnm(:)       |NA            |four character code to represent crop name
!!    deepst(:)     |mm H2O        |depth of water in deep aquifer
!!    hru_km(:)     |km^2          |area of HRU in square kilometers
!!    hru_sub(:)    |none          |subbasin in which HRU is located
!!    hrugis(:)     |none          |GIS code printed to output files(output.hru,.rch)
!!    hruyro(1,:)   |mm H2O        |precipitation in HRU during year
!!    hruyro(2,:)   |mm H2O        |amount of precipitation falling as freezing
!!                                 |rain/snow in HRU during year
!!    hruyro(3,:)   |mm H2O        |amount of snow melt in HRU during year
!!    hruyro(4,:)   |mm H2O        |amount of surface runoff to main channel
!!                                 |from HRU during year (ignores impact of
!!                                 |transmission losses)
!!    hruyro(5,:)   |mm H2O        |amount of lateral flow contribution to main
!!                                 |channel from HRU during year
!!    hruyro(6,:)   |mm H2O        |amount of groundwater flow contribution to
!!                                 |main channel from HRU during year
!!    hruyro(7,:)   |mm H2O        |amount of water moving from shallow aquifer
!!                                 |to plants or soil profile in HRU during year
!!    hruyro(8,:)   |mm H2O        |amount of water recharging deep aquifer in
!!                                 |HRU during year
!!    hruyro(9,:)   |mm H2O        |total amount of water entering both aquifers
!!                                 |from HRU during year
!!    hruyro(10,:)  |mm H2O        |water yield (total amount of water entering
!!                                 |main channel) from HRU during year
!!    hruyro(11,:)  |mm H2O        |amount of water percolating out of the soil
!!                                 |profile and into the vadose zone in HRU
!!                                 |during year
!!    hruyro(12,:)  |mm H2O        |actual evapotranspiration in HRU during year
!!    hruyro(13,:)  |mm H2O        |amount of transmission losses from tributary
!!                                 |channels in HRU for year
!!    hruyro(14,:)  |metric tons/ha|sediment yield from HRU for year
!!    hruyro(17,:)  |kg N/ha       |amount of nitrogen applied in continuous
!!                                 |fertilizer operation during year in HRU
!!    hruyro(18,:)  |kg P/ha       |amount of phosphorus applied in continuous
!!                                 |fertilizer operation during year in HRU
!!    hruyro(23,:)  |mm H2O        |amount of water removed from shallow aquifer
!!                                 |in HRU for irrigation during year
!!    hruyro(24,:)  |mm H2O        |amount of water removed from deep aquifer
!!                                 |in HRU for irrigation during year
!!    hruyro(25,:)  |mm H2O        |potential evapotranspiration in HRU during
!!                                 |year
!!    hruyro(26,:)  |kg N/ha       |annual amount of N (organic & mineral)
!!                                 |applied in HRU during grazing
!!    hruyro(27,:)  |kg P/ha       |annual amount of P (organic & mineral)
!!                                 |applied in HRU during grazing
!!    hruyro(28,:)  |kg N/ha       |annual amount of N (organic & mineral)
!!                                 |auto-applied in HRU
!!    hruyro(29,:)  |kg P/ha       |annual amount of P (organic & mineral)
!!                                 |auto-applied in HRU
!!    hruyro(31,:)  |stress days   |water stress days in HRU during year
!!    hruyro(32,:)  |stress days   |temperature stress days in HRU during year
!!    hruyro(33,:)  |stress days   |nitrogen stress days in HRU during year
!!    hruyro(34,:)  |stress days   |phosphorus stress days in HRU during year
!!    hruyro(35,:)  |kg N/ha       |organic nitrogen in surface runoff in HRU
!!                                 |during year
!!    hruyro(36,:)  |kg P/ha       |organic phosphorus in surface runoff in HRU
!!                                 |during year
!!    hruyro(37,:)  |kg N/ha       |nitrate in surface runoff in HRU during year
!!    hruyro(38,:)  |kg N/ha       |nitrate in lateral flow in HRU during year
!!    hruyro(39,:)  |kg P/ha       |soluble phosphorus in surface runoff in HRU
!!                                 |during year
!!    hruyro(40,:)  |kg N/ha       |amount of nitrogen removed from soil by plant
!!                                 |uptake in HRU during year
!!    hruyro(41,:)  |kg N/ha       |nitrate percolating past bottom of soil
!!                                 |profile in HRU during year
!!    hruyro(42,:)  |kg P/ha       |amount of phosphorus removed from soil by
!!                                 |plant uptake in HRU during year
!!    hruyro(43,:)  |kg P/ha       |amount of phosphorus moving from labile
!!                                 |mineral to active mineral pool in HRU during
!!                                 |year
!!    hruyro(44,:)  |kg P/ha       |amount of phosphorus moving from active
!!                                 |mineral to stable mineral pool in HRU during
!!                                 |year
!!    hruyro(45,:)  |kg N/ha       |amount of nitrogen applied to HRU in
!!                                 |fertilizer and grazing operations during
!!                                 |year
!!    hruyro(46,:)  |kg P/ha       |amount of phosphorus applied to HRU in
!!                                 |fertilizer and grazing operations during
!!                                 |year
!!    hruyro(47,:)  |kg N/ha       |amount of nitrogen added to soil by fixation
!!                                 |in HRU during year
!!    hruyro(48,:)  |kg N/ha       |amount of nitrogen lost by denitrification
!!                                 |in HRU during year
!!    hruyro(49,:)  |kg N/ha       |amount of nitrogen moving from active organic
!!                                 |to nitrate pool in HRU during year
!!    hruyro(50,:)  |kg N/ha       |amount of nitrogen moving from active organic
!!                                 |to stable organic pool in HRU during year
!!    hruyro(51,:)  |kg P/ha       |amount of phosphorus moving from organic to
!!                                 |labile mineral pool in HRU during year
!!    hruyro(52,:)  |kg N/ha       |amount of nitrogen moving from fresh organic
!!                                 |to nitrate and active organic pools in HRU
!!                                 |during year
!!    hruyro(53,:)  |kg P/ha       |amount of phosphorus moving from fresh
!!                                 |organic to the labile mineral and organic
!!                                 |pools in HRU during year
!!    hruyro(54,:)  |kg N/ha       |amount of nitrogen added to soil in rain
!!                                 |during year
!!    hruyro(61,:)  |metric tons/ha|daily soil loss predicted with USLE equation
!!    hruyro(63,:)  |# bacteria/ha |less persistent bacteria transported to main
!!                                 |channel from HRU during year
!!    hruyro(64,:)  |# bacteria/ha |persistent bacteria transported to main
!!                                 |channel from HRU during year
!!    hruyro(65,:)  |kg N/ha       |nitrate loading from groundwater in HRU to
!!                                 |main channel during year
!!    hruyro(66,:)  |kg P/ha       |soluble P loading from groundwater in HRU to
!!                                 |main channel during year
!!    hruyro(67,:)  |kg P/ha       |loading of mineral P attached to sediment
!!                                 |in HRU to main channel during year
!!    icr(:)        |none          |sequence number of crop grown within the
!!                                 |current year
!!    idplt(:)      |none          |land cover code from crop.dat
!!    ipdvas(:)     |none          |output variable codes for output.hru file
!!    isproj        |none          |special project code:
!!                                 |1 test rewind (run simulation twice)
!!    itots         |none          |number of output variables printed (output.hru)
!!    iyr           |year          |current year of simulation (eg 1980)
!!    lai_yrmx(:)   |none          |maximum leaf area index for the year in the
!!                                 |HRU
!!    mhruo         |none          |maximum number of variables written to
!!                                 |HRU output file (output.hru)
!!    nhru          |none          |number of HRUs in watershed
!!    nmgt(:)       |none          |management code (for GIS output only)
!!    nro(:)        |none          |sequence number of year in rotation
!!    shallst(:)    |mm H2O        |depth of water in shallow aquifer
!!    sol_sw(:)     |mm H2O        |amount of water stored in the soil profile
!!                                 |on any given day
!!    yldanu(:)     |metric tons/ha|annual yield (dry weight) in the HRU
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~

!!    ~ ~ ~ LOCAL DEFINITIONS ~ ~ ~
!!    name        |units         |definition
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
!!    ii          |none          |counter
!!    j           |none          |HRU number
!!    pdvas(:)    |varies        |array to hold HRU output values
!!    pdvs(:)     |varies        |array to hold selected HRU output values
!!                               |when user doesn't want to print all
!!    sb          |none          |subbasin number
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~

!!    ~ ~ ~ ~ ~ ~ END SPECIFICATIONS ~ ~ ~ ~ ~ ~

      use parm

      integer :: j, sb, ii, iflag
      real, dimension (mhruo) :: pdvas, pdvs
      character (len=4) :: cropname

      do j = 1, nhru
        iflag = 0
	    sb = hru_sub(j)

        do ii = 1, itoth
          if (ipdhru(ii) == j) iflag = 1
        end do

        if (iflag == 1) then

        pdvas = 0.
        pdvs = 0.

        pdvas(1) = hruyro(1,j)
        pdvas(2) = hruyro(2,j)
        pdvas(3) = hruyro(3,j)
        pdvas(4) = hruyro(22,j)
        pdvas(5) = hruyro(25,j)
        pdvas(6) = hruyro(12,j)
        pdvas(7) = hruyro(21,j) / Real(366 - leapyr)
        pdvas(8) = sol_sw(j)
        pdvas(9) = hruyro(11,j)
        pdvas(10) = hruyro(9,j)
        pdvas(11) = hruyro(8,j)
        pdvas(12) = hruyro(7,j)
        pdvas(13) = hruyro(23,j)
        pdvas(14) = hruyro(24,j)
        pdvas(15) = shallst(j)
        pdvas(16) = deepst(j)
        pdvas(17) = hruyro(19,j)
        pdvas(18) = hruyro(4,j)
        pdvas(19) = hruyro(13,j)
        pdvas(20) = hruyro(5,j)
        pdvas(21) = hruyro(6,j)
        pdvas(22) = hruyro(10,j)
        pdvas(23) = hruyro(20,j) / Real(366 - leapyr)
        pdvas(24) = hruyro(57,j) / Real(366 - leapyr)
        pdvas(25) = hruyro(55,j) / Real(366 - leapyr)
        pdvas(26) = hruyro(56,j) / Real(366 - leapyr)
        pdvas(27) = hruyro(30,j) / Real(366 - leapyr)
        pdvas(28) = hruyro(58,j) / Real(366 - leapyr)
        pdvas(29) = hruyro(14,j)
        pdvas(30) = hruyro(61,j)
        pdvas(31) = hruyro(45,j)
        pdvas(32) = hruyro(46,j)
        pdvas(33) = hruyro(28,j)
        pdvas(34) = hruyro(29,j)
        pdvas(35) = hruyro(26,j)
        pdvas(36) = hruyro(27,j)
        pdvas(37) = hruyro(17,j)
        pdvas(38) = hruyro(18,j)
        pdvas(39) = hruyro(54,j)
        pdvas(40) = hruyro(47,j)
        pdvas(41) = hruyro(52,j)
        pdvas(42) = hruyro(49,j)
        pdvas(43) = hruyro(50,j)
        pdvas(44) = hruyro(53,j)
        pdvas(45) = hruyro(51,j)
        pdvas(46) = hruyro(43,j)
        pdvas(47) = hruyro(44,j)
        pdvas(48) = hruyro(48,j)
        pdvas(49) = hruyro(40,j)
        pdvas(50) = hruyro(42,j)
        pdvas(51) = hruyro(35,j)
        pdvas(52) = hruyro(36,j)
        pdvas(53) = hruyro(67,j)
        pdvas(54) = hruyro(37,j)
        pdvas(55) = hruyro(38,j)
        pdvas(56) = hruyro(41,j)
        pdvas(57) = hruyro(65,j)
        pdvas(58) = hruyro(39,j)
        pdvas(59) = hruyro(66,j)
        pdvas(60) = hruyro(31,j)
        pdvas(61) = hruyro(32,j)
        pdvas(62) = hruyro(33,j)
        pdvas(63) = hruyro(34,j)
        pdvas(64) = bio_yrms(j)
        pdvas(65) = lai_yrmx(j)
        pdvas(66) = yldanu(j)
        pdvas(67) = hruyro(63,j)
        pdvas(68) = hruyro(64,j)
        pdvas(69) = wtab(j)  !! based on 30 day antecedent climate(mm) (prec,et)
!       pdvas(70) = wtabelo   !! based on depth from soil surface (mm)
        pdvas(70) = wat_tbl(j)   !! based on depth from soil surface (mm): Dmoriasi 4/08/2014
!!      added current snow content in the hru (not summed)
        pdvas(71) = sno_hru(j)

!!      added current soil carbon for first layer
        pdvas(72) = cmup_kgh(j)    !! first soil layer only
!!      added current soil carbon integrated - ggreagating all soil layers
        pdvas(73) = cmtot_kgh(j)

!!    adding qtile to output.hru write 3/2/2010 gsm
        pdvas(74) = hruyro(62,j)
!!    tileno3 - output.hru
        pdvas(75) = hruyro(68,j)
!!    latno3 - output.hru
        pdvas(76) = hruyro(69,j)
!!    gwq deep
        pdvas(77) = hruyro(70,j)
!!    latq contribution
        pdvas(78) = hruyro(71,j)
!!      phos due to crack flow (tvap)
        pdvas(79) = hruyro(72,j)

# GWP edits
        pdvas(80) = hruyro(73,j) / Real(366 - leapyr)
        pdvas(81) = hruyro(74,j) / Real(366 - leapyr)
        pdvas(82) = hruyro(75,j) / Real(366 - leapyr)
        pdvas(83) = hruyro(76,j) / Real(366 - leapyr)
        pdvas(84) = hruyro(77,j) / Real(366 - leapyr)
        pdvas(85) = hruyro(78,j) / Real(366 - leapyr)
        pdvas(86) = hruyro(79,j) / Real(366 - leapyr)
        pdvas(87) = hruyro(80,j) / Real(366 - leapyr)
        pdvas(88) = hruyro(81,j) / Real(366 - leapyr)
        pdvas(89) = hruyro(82,j) / Real(366 - leapyr)

# end GWP edits

        if (ipdvas(1) > 0) then
          do ii = 1, itots
            pdvs(ii) = pdvas(ipdvas(ii))
          end do

        idplant = idplt(j)
        if (idplant > 0) then
          cropname = cpnm(idplant)
        else
          cropname = "NOCR"
        endif

          if (iscen == 1 .and. isproj == 0) then
          write (28,1000) cropname, j, subnum(j), hruno(j), sb,
     &             nmgt(j), iyr, hru_km(j), (pdvs(ii), ii = 1, itots)
          else if (isproj == 1) then
          write (21,1000) cropname, j, subnum(j), hruno(j),
     &            sb, nmgt(j), iyr, hru_km(j), (pdvs(ii), ii = 1, itots)
          else if (iscen == 1 .and. isproj == 2) then
          write (28,2000) cropname, j, subnum(j), hruno(j), sb,
     &    nmgt(j), iyr, hru_km(j), (pdvs(ii), ii = 1, itots), iyr
          endif
        else
          if (iscen == 1 .and. isproj == 0) then
          write (28,1001) cropname, j, subnum(j), hruno(j), sb,
     &            nmgt(j), iyr, hru_km(j), (pdvas(ii), ii = 1, mhruo)
          else if (isproj == 1) then
          write (21,1001) cropname, j, subnum(j), hruno(j),
     &           sb, nmgt(j), iyr, hru_km(j), (pdvas(ii), ii = 1, mhruo)
          else if (iscen == 1 .and. isproj == 2) then
          write (28,1001) cropname, j, subnum(j), hruno(j), sb,
     &    nmgt(j), iyr, hru_km(j), (pdvas(ii), ii = 1, mhruo), iyr
          endif
        end if
        end if
      end do

      return

 1000 format (a4,i5,1x,a5,a4,i5,1x,i4,1x,i4,e10.5,66f10.3,1x,
     *e10.5,1x,e10.5,8e10.3,6f10.3)
 2000 format (a4,i5,1x,a5,a4,i5,1x,i4,1x,i4,e10.5,66f10.3,1x,
     *e10.5,1x,e10.5,5e10.3,6f10.3,1x,i4)
 1001 format (a4,i7,1x,a5,a4,i5,1x,i4,1x,i4,e10.5,66f10.3,1x,
     *e10.5,1x,e10.5,3e10.3,6f10.3,1x,i4)
      end


! Copyright (c) 2019, NVIDIA CORPORATION.  All rights reserved.
!
! Licensed under the Apache License, Version 2.0 (the "License");
! you may not use this file except in compliance with the License.
! You may obtain a copy of the License at
!
!     http://www.apache.org/licenses/LICENSE-2.0
!
! Unless required by applicable law or agreed to in writing, software
! distributed under the License is distributed on an "AS IS" BASIS,
! WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
! See the License for the specific language governing permissions and
! limitations under the License.

! Verify miscellaneous bugs

! The function result must be declared after the dummy arguments
module m1
contains
  function f1(x) result(y)
    integer :: x(:)
    integer :: y(size(x))
  end
  function f2(x)
    integer :: x(:)
    integer :: f2(size(x))
  end
end

!Expect: m1.mod
!module m1
!contains
! function f1(x) result(y)
!  integer(4)::x(:)
!  integer(4)::y(1_8:1_8*size(x,dim=1))
! end
! function f2(x)
!  integer(4)::x(:)
!  integer(4)::f2(1_8:1_8*size(x,dim=1))
! end
!end

! Order of names in PUBLIC statement shouldn't affect .mod file.
module m2
  public :: a
  type t
  end type
  type(t), parameter :: a = t()
end

!Expect: m2.mod
!module m2
! type::t
! end type
! type(t),parameter::a=t()
!end

! Don't write out intrinsics
module m3a
  integer, parameter :: i4 = selected_int_kind(9)
end
module m3b
  use m3a
  integer(i4) :: j
end

!Expect: m3a.mod
!module m3a
! integer(4),parameter::i4=4_4
!end

!Expect: m3b.mod
!module m3b
! use m3a,only:i4
! integer(4)::j
!end

! Test that character literals written with backslash escapes are read correctly.
module m4a
  character(1), parameter :: a = achar(1)
end
module m4b
  use m4a
  character(1), parameter :: b = a
end

!Expect: m4a.mod
!module m4a
! character(1_4,1),parameter::a=1_"\001"
!end

!Expect: m4b.mod
!module m4b
! use m4a,only:a
! character(1_4,1),parameter::b=1_"\001"
!end


text:
C
C
C +-----------------------------------------------------------------+
C |                                                                 |
C |                Copyright (C) 1986 by UCAR                       |
C |        University Corporation for Atmospheric Research          |
C |                    All Rights Reserved                          |
C |                                                                 |
C |                 NCARGRAPHICS  Version 1.00                      |
C |                                                                 |
C +-----------------------------------------------------------------+
C
C
c + noao: block data gridt changed to run time initialization
c       BLOCK DATA GRIDT
        subroutine gridt
C
C
        COMMON /LAB/ SIZX,SIZY,XDEC,YDEC,IXORI
        COMMON /CLAB/ XFMT, YFMT
        COMMON /TICK/ MAJX, MINX, MAJY, MINY
        COMMON /GRIINT/ IGRIMJ, IGRIMN, IGRITX
        CHARACTER*8 XFMT,YFMT
        REAL MAJX,MINX,MAJY,MINY
C
c +noao: following flag added to prevent initializing more than once
        logical first
        SAVE
        data first /.true./
        if (.not. first) then
            return
        endif
        first = .false.
C
c       DATA XFMT,YFMT /'(E10.3) ','(E10.3) '/
        XFMT = '(E10.3) '
        YFMT = '(E10.3) '
c
c       DATA SIZX,SIZY / 0.01, 0.01 /
        SIZX = 0.01
        SIZY = 0.01
c
c       DATA XDEC,YDEC / 0., 0. /
        XDEC = 0.
        YDEC = 0.
c
c       DATA IXORI / 0 /
        IXORI = 0
c
c       DATA MAJX,MINX,MAJY,MINY / 0., 0., 0., 0./
        MAJX = 0.
        MINX = 0.
        MAJY = 0.
        MINY = 0.
c
c       DATA IGRIMJ,IGRIMN,IGRITX / 1, 1, 1/
c+noao: These values changed so major axes and labels are bold
        IGRIMJ = 2
        IGRIMN = 1
        IGRITX = 2
C - noao
        END
C REVISION HISTORY---------------
C----------------------------------------------------------

max_stars_repo_name: ysBach/irafdocgen

max_issues_repo_name: kirxkirx/iraf

max_forks_repo_name: kirxkirx/iraf


!! module for hilbert transform
module m_hilbert_transform
    use constants, only : CUSTOM_REAL
    implicit none

    complex :: CI=(0.0,1.0)

contains

    subroutine hilbert(trace,nt)
        !! return hilbert transform of real signal trace(nt)
        !! i.e. imaginary part of analytic signal
        !! a = cmplx(trace, hilbert(trace))

        integer,                  intent(in)    :: nt
        real(kind=CUSTOM_REAL),         intent(inout) :: trace(nt)

        complex, allocatable, dimension(:) :: C
        integer :: NPT,IMID

        ! extend nt to a power of 2
        IF ( nt <= 0 ) STOP 'FATAL ERROR in HILBERT: nt must be positive'

        NPT = 2**( INT( LOG10( REAL( nt ) ) / 0.30104 ) + 1 )
        ! IF ( NPT /= nt) print*,'pad trace from length ', nt, ' to ',NPT
        IF (NPT > 16784) STOP 'FATAL ERROR in HILBERT: nt(NPT) exceeds 16784 '

        allocate(C(NPT))
        C=cmplx(0.,0.)
        C(1:nt)=cmplx(trace(1:nt),0.0)

        ! Fourier transform
        call CFFT(C,NPT,1)
        ! scaling
        C=C/NPT

        !  Multiply by i * sgn( f )
        IMID = NPT / 2
        C( 1:IMID-1 ) = -CI * C( 1:IMID-1 )   ! pos. spectrum (-i)
        C( IMID     ) = 0.0                   ! d.c. component
        C(IMID+1:NPT) = CI * C( IMID+1:NPT )   ! neg. spectrum (i)

        ! inverse Fourier transform
        call  CFFT(C,NPT,-1)

        ! output
        trace(1:nt)=real(C(1:nt))

        deallocate(C)

    end subroutine hilbert

    SUBROUTINE CFFT(trace,N,iforw)
        !!! complex FFT
        IMPLICIT NONE
        integer,         intent(in)    :: N, iforw
        complex,         intent(inout) :: trace(N)
        INTEGER :: I1, I2A, I2B, I3, I3Rev, IP1, IP2, ISign
        REAL    ::  theta, sinth
        COMPLEX :: TEMP, W, WSTP

        ISIGN = -IFORW
        I3REV = 1

        DO I3 = 1, N
        IF ( I3 < I3REV ) THEN   ! switch values
            TEMP          = trace( I3 )
            trace( I3 )    = trace( I3REV )
            trace( I3REV ) = TEMP
        ENDIF
        ! following loop is just to compute I3REV
        IP1 = N / 2
        DO WHILE (  I3REV > IP1 )
        IF ( IP1 <= 1 ) EXIT
        I3REV = I3REV - IP1
        IP1   = IP1 / 2
        END DO
        I3REV = I3REV + IP1
        END DO

        IP1 = 1

        DO WHILE  ( IP1 < N )
        IP2   = IP1 * 2
        THETA = 6.283185307 / FLOAT( ISIGN * IP2 )
        SINTH = SIN( THETA / 2.)
        WSTP  = CMPLX( -2. * SINTH * SINTH, SIN( THETA ) )
        W     = 1.

        DO I1 = 1, IP1
        DO I3 = I1, N, IP2
        I2A = I3
        I2B = I2A + IP1
        TEMP      = W * trace( I2B )
        trace( I2B ) = trace( I2A ) - TEMP
        trace( I2A ) = trace( I2A ) + TEMP
        END DO

        W = W * WSTP + W
        END DO

        IP1 = IP2
        END DO

        RETURN
    END SUBROUTINE CFFT

end module m_hilbert_transform


text:
This editor can edit this entry and tell us a bit about themselves by clicking the Edit icon.

20080821 14:04:04 nbsp Welcome to the Wiki Howdy, Mr. or Ms. Spears, and Welcome to the Wiki!     My names Evan, pleased to meet you!  Users/JabberWokky Evan JabberWokky Edwards

20080821 12:05:03 nbsp What up with you boss! Hit me up next time your in Dtown! Users/StevenDaubert


text:
c     # ------------------------------------------------------------------
c     # MAPC2M_PILLOWSPHERE
c     # ------------------------------------------------------------------
c     #
c     # Maps a logically rectangular Cartesian grid in [-1,1]x[-1,1] to
c     # either the upper hemisphere (blockno == 0) or the lower hemisphere
c     # (blockno == 1).
c     #
c     # ------------------------------------------------------------------
      subroutine mapc2m_pillowsphere(blockno, xc1,yc1,xp,yp,zp)
      implicit none

      external mapc2m_cart
      double precision xc1,yc1, xp, yp, zp

      double precision x1, y1, d, rp2, xc, yc, zc
      integer blockno
      logical l1, l2
      logical ispillowsphere

c     # Map to [-1,1]x[-1,1]
      call mapc2m_cart(blockno,xc1,yc1,xc,yc,zc)

c     # Map xc and yc from ghost cells to interior
      d = max(xc - 1.0,0.d0) + max(-1 - xc,0.d0)
      x1 = ((1-d)/(1+d))*xc

      d = max(yc - 1,0.d0) + max(-1 - yc,0.d0)
      y1 = ((1-d)/(1+d))*yc

c     # Get circle of radius sqrt(2.d0).  Cluster radial
c     # direction towards boundary
      call mapc2p_circle_sp(x1,y1,xp,yp)

c     # Set z value
      rp2 = xp**2 + yp**2
      if (abs(rp2 - 2.d0) .lt. 1d-10) then
         zp = 0.d0
      else
         zp = sqrt(2.d0 - rp2)
      endif

c     # Values that are outside of [-1,1]x[-1,1] (in ghost cell regions)
c     # to the lower hemisphere
c     # Juqueen complained about the 'xor', because it is not in the
c     # Fortran standard.  Use .neqv. instead.
      l1 = abs(xc) .gt. 1
      l2 = abs(yc) .gt. 1
      if (l1 .neqv. l2) then
         zp = -zp
      endif

c     # This maps everything to the unit sphere
      xp = xp/sqrt(2.d0)
      yp = yp/sqrt(2.d0)
      zp = zp/sqrt(2.d0)

c     # Set lower hemisphere
      if (ispillowsphere()) then
         if (blockno .eq. 1) then
            zp = -zp
         endif
      endif

      return
      end
c --------------------------------- MAPC2M -------------------------------

c     # Map single grid to the disk.  Since this disk will be used for the
c     #
      subroutine mapc2p_circle_sp(x1,y1,xp,yp)
      implicit none

      double precision xc,yc,xp,yp, x1,y1
      double precision xi,eta,x,y, minxy,maxxy
      double precision xit, etat, dd

      double precision pi, pi2
      common /compi/ pi, pi2

      xc = x1
      yc = y1

      if (x1 .lt. -1.d0) then
         xc = -(xc + 2)
      endif

      xi = min(abs(xc),abs(yc))
      eta = max(abs(xc),abs(yc))
      eta = max(eta,1.d-10)

      dd = sin(pi*eta/2.d0)
c      dd = eta*(2-eta)

      xit = (xi/eta)*dd
      etat = dd

      call map_north_sector_sp(xit,etat,x,y)

      minxy = min(abs(x),abs(y))
      maxxy = max(abs(x),abs(y))

      if (abs(xc) .le. abs(yc))  then
c        # In NS sectors
         xp = sign(1.d0,xc)*minxy
         yp = sign(1.d0,yc)*maxxy
      else
c        # In EW sectors
         xp = sign(1.d0,xc)*maxxy
         yp = sign(1.d0,yc)*minxy
      endif

      end

      subroutine map_north_sector_sp(xi,eta,x,y)
      implicit none

      double precision xi,eta,x,y

      x = xi
      y = sqrt(2 - xi**2) - sqrt(2 - eta**2) + eta

      end

max_stars_repo_name: ECLAIRWaveS/ForestClaw

max_issues_repo_name: ECLAIRWaveS/ForestClaw

max_forks_repo_name: ECLAIRWaveS/ForestClaw


text:
C    Copyright(C) 1999-2021 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
C
C     C* FILE: [.QMESH]PICKM1.FOR
C     C* MODIFIED BY: TED BLACKER
C     C* MODIFICATION DATE: 7/6/90
C     C* MODIFICATION: COMPLETED HEADER INFORMATION
C
      SUBROUTINE PICKM1 (N, X, Y, ANGLE, M1, IFIRST, REAL)
C***********************************************************************
C
C     SUBROUTINE PICKM1 = DETERMINES A REASONABLE SHAPE FOR A LOGICAL
C     RECTANGLE WITH PERIMETER GIVEN IN X AND Y
C
C***********************************************************************
C
      DIMENSION X (N), Y (N), ANGLE (N)
      DIMENSION SMANG (7), INDEX (7)
C
      LOGICAL REAL
C
C     FORM THE LIST OF SMALLEST ANGLES
C
      NSA = 6
      DO 100 I = 1, NSA
         SMANG (I) = 10.
         INDEX (I) = 0
 100  CONTINUE
C
      PI = ATAN2(0.0, -1.0)
      TWOPI = PI + PI
      AGOLD = ATAN2 (Y (1) - Y (N), X (1) - X (N))
C
      DO 130 J = 1, N
C
C     GET THE ANGLE FORMED BY THIS SET OF POINTS
C
         NEXT = J + 1
         IF  (NEXT .GT. N) NEXT = 1
         AGNEW = ATAN2 (Y (NEXT) - Y (J) ,  X (NEXT) - X (J))
         DIFF = AGNEW - AGOLD
         IF (DIFF .GT. PI)DIFF = DIFF - TWOPI
         IF (DIFF .LT.  - PI)DIFF = DIFF + TWOPI
         ANGLE (J) = PI - DIFF
         AGOLD = AGNEW
C
C     SORT THIS ANGLE AGAINST PREVIOUS ANGLES TO SEE IF IT IS ONE OF
C     THE SMALLEST
C
         SMANG (NSA + 1) = ANGLE (J)
         INDEX (NSA + 1) = J
         DO II = 1, NSA
            I = NSA + 1 - II
            IF  (SMANG (I + 1) .GE. SMANG (I)) GO TO 120
            TEMP = SMANG (I)
            ITEMP = INDEX (I)
            SMANG (I) = SMANG (I + 1)
            INDEX (I) = INDEX (I + 1)
            SMANG (I + 1) = TEMP
            INDEX (I + 1) = ITEMP
         end do
 120     CONTINUE
C
 130  CONTINUE
C
C     DETERMINE OPTIMUM ORIGIN / SHAPE COMBINATION
C
      ATOL = PI * 150. / 180.
      IFIRST = 1
      M1 = N / 4
      M2 = N / 2 - M1
      I2 = 1 + M1
      I3 = I2 + M2
      I4 = I3 + M1
      GBEST = ANGLE (1) + ANGLE (I2) + ANGLE (I3) + ANGLE (I4)
      BADANG = AMAX1 (ANGLE (1), ANGLE (I2), ANGLE (I3), ANGLE (I4))
C
      MMAX = N / 2 - 1
      AMAXEL = DBLE(N / 4) * DBLE( (N + 2) / 4)
      DO 150 ISA = 1, NSA
         IF (SMANG (ISA) .LE. ATOL) THEN
            I1 = INDEX (ISA)
            DO 140 M = 1, MMAX
               M2 = N / 2 - M
               I2 = I1 + M
               IF  (I2 .GT. N) I2 = I2 - N
               I3 = I2 + M2
               IF  (I3 .GT. N) I3 = I3 - N
               I4 = I3 + M
               IF  (I4 .GT. N) I4 = I4 - N
               AFAC = ANGLE (I1) + ANGLE (I2) + ANGLE (I3) + ANGLE (I4)
               ERAT = AMIN1 (AMAXEL / DBLE(M * M2) ,  5.)
               EFAC =  (ERAT + 15.) / 16.
               GVAL = AFAC * EFAC
               IF (GVAL .LT. GBEST) THEN
                  BADANG = AMAX1 (ANGLE (I1), ANGLE (I2), ANGLE (I3),
     &                 ANGLE (I4))
                  IFIRST = I1
                  M1 = M
                  GBEST = GVAL
               ENDIF
 140        CONTINUE
         ENDIF
 150  CONTINUE
      IF ( (REAL) .AND. (BADANG .GT. 2.62)) THEN
         CALL MESSAGE(' **  WARNING: CORNER (S) OF THE REGION HAVE  **')
         CALL MESSAGE(' **           LARGE ANGLES  (> 150 DEGREES.) **')
         CALL MESSAGE(' **           POORLY FORMED MESH MAY RESULT  **')
      ENDIF
C
      RETURN
      END

max_stars_repo_name: tjotaha/seacas

max_issues_repo_name: tjotaha/seacas

max_forks_repo_name: tjotaha/seacas