text:
subroutine write_swan_datafile (var1  , var2       , mmax   , nmax, covered, &
                              & filnam, extr_var1, extr_var2, &
                              & sumvars, positiveonly)
!----- GPL ---------------------------------------------------------------------
!
!  Copyright (C)  Stichting Deltares, 2011-2016.
!
!  This program is free software: you can redistribute it and/or modify
!  it under the terms of the GNU General Public License as published by
!  the Free Software Foundation version 3.
!
!  This program is distributed in the hope that it will be useful,
!  but WITHOUT ANY WARRANTY; without even the implied warranty of
!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!  GNU General Public License for more details.
!
!  You should have received a copy of the GNU General Public License
!  along with this program.  If not, see <http://www.gnu.org/licenses/>.
!
!  contact: delft3d.support@deltares.nl
!  Stichting Deltares
!  P.O. Box 177
!  2600 MH Delft, The Netherlands
!
!  All indications and logos of, and references to, "Delft3D" and "Deltares"
!  are registered trademarks of Stichting Deltares, and remain the property of
!  Stichting Deltares. All rights reserved.
!
!-------------------------------------------------------------------------------
!  $Id: write_swan_datafile.f90 5717 2016-01-12 11:35:24Z mourits $
!  $HeadURL: https://svn.oss.deltares.nl/repos/delft3d/tags/6686/src/engines_gpl/wave/packages/io/src/write_swan_datafile.f90 $
!!--description-----------------------------------------------------------------
! NONE
!!--pseudo code and references--------------------------------------------------
! NONE
!!--declarations----------------------------------------------------------------
    implicit none
!
! Local parameters
!
    integer     , parameter :: iindex = 1
    integer     , parameter :: jindex = 2
!
! Global variables
!
    integer                       , intent(in)  :: mmax
    integer                       , intent(in)  :: nmax
    integer, dimension(mmax, nmax), intent(in)  :: covered
    real   , dimension(mmax, nmax)              :: var1
    real   , dimension(mmax, nmax)              :: var2
    logical                       , intent(in)  :: extr_var1
    logical                       , intent(in)  :: extr_var2
    logical                       , intent(in)  :: sumvars
    logical                       , intent(in)  :: positiveonly
    character(*)                  , intent(in)  :: filnam
!
! Local variables
!
    integer, dimension(:,:,:), allocatable :: closestPoint
    integer                                :: dm
    integer                                :: i
    integer                                :: j
    integer                                :: i1
    integer                                :: j1
    integer                                :: ierr
    integer                                :: lunfil
    integer, external                      :: new_lun
    integer                                :: sweep
    integer                                :: sweepEnd
    integer                                :: sweepStart
    integer                                :: sweepStep
    real   , dimension(:,:)  , allocatable :: cpDistance
    real                                   :: distance
!
!! executable statements -------------------------------------------------------
!
    if (extr_var1 .or. extr_var2) then
       allocate(closestPoint(mmax, nmax, max(iindex,jindex)))
       allocate(cpDistance  (mmax, nmax))
       do j=1,nmax
          do i=1,mmax
             if (covered(i, j) == 0) then
                !
                ! point is outside all FLOW grids:
                ! - closest point is not defined (zero)
                ! - distance to closest point is infinity
                !
                closestPoint(i, j, iindex) = 0
                closestPoint(i, j, jindex) = 0
                cpDistance  (i, j)         = 1.0e20
             else
                !
                ! point is inside a FLOW grid:
                ! - closest point is the current point itself
                ! - distance to closest point is zero
                !
                closestPoint(i, j, iindex) = i
                closestPoint(i, j, jindex) = j
                cpDistance  (i, j)         = 0.0
             endif
          enddo
       enddo
       !
       ! Detect closest covered point for all not covered points
       !
       do sweep = 1, 2
          if (sweep == 1) then
             !
             ! From top to bottom
             !
             sweepStart =  1
             sweepEnd   = nmax
             sweepStep  =  1
          else
             !
             ! From bottom to top
             !
             sweepStart = nmax
             sweepEnd   =  1
             sweepStep  = -1
          endif
          do j = sweepStart, sweepEnd, sweepStep
             do i = 1, mmax
                if (covered(i, j) == 0) then
                   !
                   ! point is outside all FLOW grids
                   !
                   do dm = -1, 1
                      i1 = i + dm
                      j1 = j - sweepStep
                      if (      i1>=1 .and. i1<=mmax &
                        & .and. j1>=1 .and. j1<=nmax  ) then
                         if (closestPoint(i1, j1, 1)>0) then
                            distance =   (i - closestPoint(i1, j1, iindex))**2 &
                                     & + (j - closestPoint(i1, j1, jindex))**2
                            if (distance < cpDistance(i, j)) then
                               closestPoint(i, j, iindex) = closestPoint(i1, j1, iindex)
                               closestPoint(i, j, jindex) = closestPoint(i1, j1, jindex)
                               cpDistance  (i, j)         = distance
                            endif
                         endif
                      endif
                   enddo
                   i1 = i - 1
                   j1 = j
                   if (      i1>=1 .and. i1<=mmax &
                     & .and. j1>=1 .and. j1<=nmax  ) then
                      if (closestPoint(i1, j1, 1)>0) then
                         distance =   (i - closestPoint(i1, j1, iindex))**2 &
                                  & + (j - closestPoint(i1, j1, jindex))**2
                         if (distance < cpDistance(i, j)) then
                            closestPoint(i, j, iindex) = closestPoint(i1, j1, iindex)
                            closestPoint(i, j, jindex) = closestPoint(i1, j1, jindex)
                            cpDistance  (i, j)         = distance
                         endif
                      endif
                   endif
                endif
             enddo
             do i = mmax, 1, -1
                if (covered(i, j) == 0) then
                   !
                   ! point is outside all FLOW grids
                   !
                   i1 = i + 1
                   j1 = j
                   if (      i1>=1 .and. i1<=mmax &
                     & .and. j1>=1 .and. j1<=nmax  ) then
                      if (closestPoint(i1, j1, 1)>0) then
                         distance =   (i - closestPoint(i1, j1, iindex))**2 &
                                  & + (j - closestPoint(i1, j1, jindex))**2
                         if (distance < cpDistance(i, j)) then
                            closestPoint(i, j, iindex) = closestPoint(i1, j1, iindex)
                            closestPoint(i, j, jindex) = closestPoint(i1, j1, jindex)
                            cpDistance  (i, j)         = distance
                         endif
                      endif
                   endif
                endif
             enddo
          enddo
       enddo
       if (extr_var1) then
          do j=1,nmax
             do i=1,mmax
                !
                ! Only extrapolate from closestPoints, covered by valid points from source grid
                ! If the closestPoint is covered by invalid source grid points
                ! then covered == -1 and the value at i,j is not changed
                !
                if (covered(i, j) == 0 .and. &
                  & covered(closestPoint(i, j, iindex), closestPoint(i, j, jindex)) == 1) then
                   var1(i, j) = var1(closestPoint(i, j, iindex), closestPoint(i, j, jindex))
                endif
             enddo
          enddo
       endif
       if (extr_var2) then
          do j=1,nmax
             do i=1,mmax
                !
                ! Only extrapolate from closestPoints, covered by valid points from source grid
                ! If the closestPoint is covered by invalid source grid points
                ! then covered == -1 and the value at i,j is not changed
                !
                if (covered(i, j) == 0 .and. &
                  & covered(closestPoint(i, j, iindex), closestPoint(i, j, jindex)) == 1) then
                   var2(i, j) = var2(closestPoint(i, j, iindex), closestPoint(i, j, jindex))
                endif
             enddo
          enddo
       endif
       deallocate(closestPoint, stat=ierr)
       deallocate(cpDistance  , stat=ierr)
    endif
    !
    ! write var1 and var2 to file for SWAN
    !
    lunfil = new_lun()
    open (lunfil, file = filnam, status = 'unknown')
    !
    ! Up to now, SWAN data files are produced for writing
    ! var1 and var2 (current, wind) or var1+var2 (bottom)
    ! logical sumvars is used to distinct these two cases
    !
    if (sumvars) then
       if (positiveonly) then
          write (lunfil,'(4(3X,E13.6))') ( ( max(var1(i,j)+var2(i,j),0.0), i=1,mmax), j=1,nmax )
       else
          write (lunfil,'(4(3X,E13.6))') ( (var1(i,j)+var2(i,j), i=1,mmax), j=1,nmax )
       endif
    else
       if (positiveonly) then
          write (lunfil,'(4(3X,E13.6))') ( ( max(var1(i,j),0.0), i=1,mmax), j=1,nmax )
          write (lunfil,'(4(3X,E13.6))') ( ( max(var2(i,j),0.0), i=1,mmax), j=1,nmax )
       else
          write (lunfil,'(4(3X,E13.6))') ( (var1(i,j), i=1,mmax), j=1,nmax )
          write (lunfil,'(4(3X,E13.6))') ( (var2(i,j), i=1,mmax), j=1,nmax )
       endif
    endif
    close (lunfil)
end subroutine write_swan_datafile


text:
Tal Topf is a candidate for ASUCD ASUCD Senate Senate in the Fall 2012 ASUCD Election running on the NOW slate.

She was awarded a trip to Washington DC for her opposition to Israeli divestment from The David Project.

 Candidate Statement

Hey Everyone!

My name is TAL TOPF and I am running for ASUCD Senate with the NOW slate! I am a Communications and Psychology double major from Los Angeles. It’s THAT time of year when we get to elect the group of people who will make a difference on our campus. There are many reasons I think I should be elected an ASUCD Senator, but by far the most important is that I LOVE DAVIS! I enjoy every moment I spend on campus and this is greatly due to the amazing ASUCD resources available to me. I work at the information center at the Memorial Union and I am in a sorority on campus. Being a student here is a great experience but I think I can improve campus life even more.

If elected, here are just a few of the changes I would like to see at UC Davis:

EXPAND RESEARCH AND INTERNSHIP OPPORTUNITIES
Internship and research opportunities are important to a wellrounded college experience. Research can often guide our decisions on what we want to do after college, allow us to use what we learn in class, and get references for graduate school or careers. There are so many research and internship opportunities, but there isn’t an easy mechanism for students to find these experiences. I want to create a website that will connect students with professors looking for undergraduate interns. This website will not only provide students with more opportunities, but will also continue to make UC Davis a top research institution.

MAKE CLUBS AND ORGANIZATIONS MORE ACCESSIBLE TO STUDENTS
We are a huge and diverse campus! With so many opportunities, it becomes hard to keep up with when clubs have meetings, fundraisers, and other events. I want to assure that, even if you are not a member of a club, you will know what is going on and what kind of events there are with all different campus groups. I will do this by working with Center for Student Involvement so that it can better fit students needs. I want us to celebrate and support our entire Davis community and increase Aggie Pride.

NOW is the time to VOTE TAL TOPF #1 and the rest of the NOW slate #2#5.

VOTE NOVEMBER 1316 @ http://elections.ucdavis.edu

20121115 09:44:16 nbsp Will you support a bill that would divest ASUCD funds in companies that support the nation state Israel due to their crimes against humanity and the illegal occupation of Palestine? Users/HannaFletcher

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
!
!     CalculiX - A 3-dimensional finite element program
!              Copyright (C) 1998-2019 Guido Dhondt
!
!     This program is free software; you can redistribute it and/or
!     modify it under the terms of the GNU General Public License as
!     published by the Free Software Foundation(version 2);
!
!
!     This program is distributed in the hope that it will be useful,
!     but WITHOUT ANY WARRANTY; without even the implied warranty of
!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
!     GNU General Public License for more details.
!
!     You should have received a copy of the GNU General Public License
!     along with this program; if not, write to the Free Software
!     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
!
      subroutine contactdampings(inpc,textpart,elcon,nelcon,
     &  nmat,ntmat_,ncmat_,irstrt,istep,istat,n,iline,ipol,inl,ipoinp,
     &  inp,ipoinpc,imat,ier)
!
!     reading the input deck: *CONTACT DAMPING
!
      implicit none
!
      character*1 inpc(*)
      character*132 textpart(16)
!
      integer nelcon(2,*),nmat,ntmat_,istep,istat,ipoinpc(0:*),
     &  n,key,i,ncmat_,irstrt(*),iline,ipol,inl,ipoinp(2,*),inp(3,*),
     &  imat,ier
!
      real*8 elcon(0:ncmat_,ntmat_,*)
!
      if((istep.gt.0).and.(irstrt(1).ge.0)) then
         write(*,*) '*ERROR reading *CONTACT DAMPING:'
         write(*,*) '       *CONTACT DAMPING should be placed'
         write(*,*) '       before all step definitions'
         ier=1
         return
      endif
!
      if(imat.eq.0) then
         write(*,*) '*ERROR reading *CONTACT DAMPING:'
         write(*,*) '       *CONTACT DAMPING should be preceded'
         write(*,*) '       by a *SURFACE INTERACTION card'
         ier=1
         return
      endif
!
!     default: no tangential damping
!
      elcon(8,1,imat)=0.d0
!
      do i=2,n
         if(textpart(i)(1:16).eq.'TANGENTFRACTION=') then
            read(textpart(i)(17:36),'(f20.0)',iostat=istat)
     &               elcon(8,1,imat)
            if(istat.gt.0) then
               call inputerror(inpc,ipoinpc,iline,
     &              "*CONTACT DAMPING%",ier)
               return
            endif
         else
            write(*,*)
     &   '*WARNING reading *CONTACT DAMPING: parameter not recognized:'
            write(*,*) '         ',
     &                 textpart(i)(1:index(textpart(i),' ')-1)
            call inputwarning(inpc,ipoinpc,iline,
     &"*CONTACT DAMPING%")
         endif
      enddo
!
      nelcon(1,imat)=max(nelcon(1,imat),8)
      nelcon(2,imat)=1
!
!     no temperature dependence allowed; last line is decisive
!
      do
         call getnewline(inpc,textpart,istat,n,key,iline,ipol,inl,
     &        ipoinp,inp,ipoinpc)
         if((istat.lt.0).or.(key.eq.1)) return
         read(textpart(1)(1:20),'(f20.0)',iostat=istat)
     &        elcon(5,1,imat)
         if(istat.gt.0) then
            call inputerror(inpc,ipoinpc,iline,
     &        "*CONTACT DAMPING%",ier)
            return
         endif
         elcon(0,1,imat)=0.d0
      enddo
!
      return
      end

max_stars_repo_name: alleindrach/calculix-desktop

max_issues_repo_name: alleindrach/calculix-desktop

max_forks_repo_name: alleindrach/calculix-desktop


text:
module trace_module

  implicit none

  private

  public trace

contains

  subroutine trace(q,c,flatn,qd_l1,qd_l2,qd_h1,qd_h2, &
                   dloga,dloga_l1,dloga_l2,dloga_h1,dloga_h2, &
                   qxm,qxp,qym,qyp,qpd_l1,qpd_l2,qpd_h1,qpd_h2, &
                   src,src_l1,src_l2,src_h1,src_h2, &
                   ilo1,ilo2,ihi1,ihi2,dx,dy,dt)

    use meth_params_module, only : plm_iorder, QVAR, QRHO, QU, QV, &
                                   QREINT, QPRES, &
                                   npassive, qpass_map, small_dens, small_pres, ppm_type, use_pslope
    use slope_module, only : uslope, pslope, multid_slope
    use amrex_constants_module

    implicit none

    integer ilo1,ilo2,ihi1,ihi2
    integer qd_l1,qd_l2,qd_h1,qd_h2
    integer dloga_l1,dloga_l2,dloga_h1,dloga_h2
    integer qpd_l1,qpd_l2,qpd_h1,qpd_h2
    integer src_l1,src_l2,src_h1,src_h2

    double precision dx, dy, dt
    double precision     q(qd_l1:qd_h1,qd_l2:qd_h2,QVAR)
    double precision     c(qd_l1:qd_h1,qd_l2:qd_h2)
    double precision flatn(qd_l1:qd_h1,qd_l2:qd_h2)
    double precision dloga(dloga_l1:dloga_h1,dloga_l2:dloga_h2)
    double precision qxm(qpd_l1:qpd_h1,qpd_l2:qpd_h2,QVAR)
    double precision qxp(qpd_l1:qpd_h1,qpd_l2:qpd_h2,QVAR)
    double precision qym(qpd_l1:qpd_h1,qpd_l2:qpd_h2,QVAR)
    double precision qyp(qpd_l1:qpd_h1,qpd_l2:qpd_h2,QVAR)
    double precision src(src_l1:src_h1,src_l2:src_h2,QVAR)

    double precision, allocatable :: dqx(:,:,:), dqy(:,:,:)

    ! Local variables
    integer i, j
    integer n, ipassive

    double precision dtdx, dtdy
    double precision cc, csq, rho, u, v, p, rhoe
    double precision drho, du, dv, dp, drhoe

    double precision enth, alpham, alphap, alpha0r, alpha0e
    double precision alpha0u, alpha0v
    double precision spzero
    double precision apright, amright, azrright, azeright
    double precision azu1rght, azv1rght
    double precision apleft, amleft, azrleft, azeleft
    double precision azu1left, azv1left
    double precision acmprght, acmpleft, acmpbot, acmptop
    double precision sourcr,sourcp,source,courn,eta,dlogatmp

    double precision :: rho_ref, u_ref, v_ref, p_ref, rhoe_ref
    double precision :: e(3)

    if (ppm_type .ne. 0) then
       print *,'Oops -- shouldnt be in trace with ppm_type != 0'
       call bl_error("Error:: trace_2d.f90")
    end if

    dtdx = dt/dx
    dtdy = dt/dy

    allocate(dqx(qpd_l1:qpd_h1,qpd_l2:qpd_h2,QVAR))
    allocate(dqy(qpd_l1:qpd_h1,qpd_l2:qpd_h2,QVAR))

    ! Compute slopes
    if (plm_iorder == 1) then
       dqx(ilo1-1:ihi1+1,ilo2-1:ihi2+1,1:QVAR) = ZERO
       dqy(ilo1-1:ihi1+1,ilo2-1:ihi2+1,1:QVAR) = ZERO

    elseif (plm_iorder == 2) then
       ! these are piecewise linear slopes.  The limiter is a 4th order
       ! limiter, but the overall method will be second order.
       call uslope(q, &
                   flatn, qd_l1, qd_l2, qd_h1, qd_h2, &
                   dqx  ,qpd_l1,qpd_l2,qpd_h1,qpd_h2, &
                   ilo1,ilo2,ihi1,ihi2,QVAR,1)

       call uslope(q,flatn,qd_l1,qd_l2,qd_h1,qd_h2, &
                   dqy,qpd_l1,qpd_l2,qpd_h1,qpd_h2, &
                   ilo1,ilo2,ihi1,ihi2,QVAR,2)

       if (use_pslope .eq. 1) then
          call pslope(q(:,:,QPRES),q(:,:,QRHO),  &
                      flatn         , qd_l1, qd_l2, qd_h1, qd_h2, &
                      dqx(:,:,QPRES),qpd_l1,qpd_l2,qpd_h1,qpd_h2, &
                      src           ,src_l1,src_l2,src_h1,src_h2, &
                      ilo1,ilo2,ihi1,ihi2,dx,dy,1)

          call pslope(q(:,:,QPRES),q(:,:,QRHO),  &
                      flatn         , qd_l1, qd_l2, qd_h1, qd_h2, &
                      dqy(:,:,QPRES),qpd_l1,qpd_l2,qpd_h1,qpd_h2, &
                      src           ,src_l1,src_l2,src_h1,src_h2, &
                      ilo1,ilo2,ihi1,ihi2,dx,dy,2)
       endif

    elseif (plm_iorder == -2) then
       ! these are also piecewise linear, but it uses a multidimensional
       ! reconstruction based on the BDS advection method to construct
       ! the x- and y-slopes together
       do n = 1, QVAR
          call multid_slope(q(:,:,n), flatn, &
                            qd_l1, qd_l2, qd_h1, qd_h2, &
                            dqx(:,:,n), dqy(:,:,n), &
                            qpd_l1, qpd_l2, qpd_h1, qpd_h2, &
                            dx, dy, &
                            ilo1, ilo2, ihi1, ihi2)
       enddo

    else
       call bl_error("ERROR: invalid value of islope")

    endif

    !-------------------------------------------------------------------------
    ! x-direction
    !-------------------------------------------------------------------------

    ! Compute left and right traced states

    do j = ilo2-1, ihi2+1
       do i = ilo1-1, ihi1+1

          cc = c(i,j)
          csq = cc**2
          rho = q(i,j,QRHO)
          u = q(i,j,QU)
          v = q(i,j,QV)
          p = q(i,j,QPRES)
          rhoe = q(i,j,QREINT)
          enth = ( (rhoe+p)/rho )/csq

          drho = dqx(i,j,QRHO)
          du = dqx(i,j,QU)
          dv = dqx(i,j,QV)
          dp = dqx(i,j,QPRES)
          drhoe = dqx(i,j,QREINT)

          alpham = HALF*(dp/(rho*cc) - du)*rho/cc
          alphap = HALF*(dp/(rho*cc) + du)*rho/cc
          alpha0r = drho - dp/csq
          alpha0e = drhoe - dp*enth
          alpha0v = dv

          e(1) = u-cc
          e(2) = u
          e(3) = u+cc

          ! construct the right state on the i-1/2 interface

          rho_ref = rho - HALF*(ONE + dtdx*min(e(1),ZERO))*drho
          u_ref = u - HALF*(ONE + dtdx*min(e(1),ZERO))*du
          v_ref = v - HALF*(ONE + dtdx*min(e(1),ZERO))*dv
          p_ref = p - HALF*(ONE + dtdx*min(e(1),ZERO))*dp
          rhoe_ref = rhoe - HALF*(ONE + dtdx*min(e(1),ZERO))*drhoe

          ! this is -(1/2) ( 1 + dt/dx lambda) (l . dq) r
          apright = 0.25d0*dtdx*(e(1) - e(3))*(ONE - sign(ONE,e(3)))*alphap
          amright = 0.25d0*dtdx*(e(1) - e(1))*(ONE - sign(ONE,e(1)))*alpham

          azrright = 0.25e0*dtdx*(e(1)-e(2))*(ONE - sign(ONE,e(2)))*alpha0r
          azeright = 0.25e0*dtdx*(e(1)-e(2))*(ONE - sign(ONE,e(2)))*alpha0e
          azv1rght = 0.25e0*dtdx*(e(1)-e(2))*(ONE - sign(ONE,e(2)))*alpha0v

          if (i .ge. ilo1) then
             qxp(i,j,QRHO) = rho_ref + apright + amright + azrright
             qxp(i,j,QRHO) = max(small_dens,qxp(i,j,QRHO))
             qxp(i,j,QU) = u_ref + (apright - amright)*cc/rho
             qxp(i,j,QV) = v_ref + azv1rght
             qxp(i,j,QPRES) = p_ref + (apright + amright)*csq
             qxp(i,j,QPRES) = max(qxp(i,j,QPRES),small_pres)
             qxp(i,j,QREINT) = rhoe_ref + (apright + amright)*enth*csq + azeright
          end if

          ! construct the left state on the i+1/2 interface

          rho_ref = rho + HALF*(ONE - dtdx*max(e(3),ZERO))*drho
          u_ref = u + HALF*(ONE - dtdx*max(e(3),ZERO))*du
          v_ref = v + HALF*(ONE - dtdx*max(e(3),ZERO))*dv
          p_ref = p + HALF*(ONE - dtdx*max(e(3),ZERO))*dp
          rhoe_ref = rhoe + HALF*(ONE - dtdx*max(e(3),ZERO))*drhoe

          apleft = 0.25d0*dtdx*(e(3) - e(3))*(ONE + sign(ONE,e(3)))*alphap
          amleft = 0.25d0*dtdx*(e(3) - e(1))*(ONE + sign(ONE,e(1)))*alpham

          azrleft = 0.25d0*dtdx*(e(3) - e(2))*(ONE + sign(ONE,e(2)))*alpha0r
          azeleft = 0.25d0*dtdx*(e(3) - e(2))*(ONE + sign(ONE,e(2)))*alpha0e
          azv1left = 0.25d0*dtdx*(e(3) - e(2))*(ONE + sign(ONE,e(2)))*alpha0v

          if (i .le. ihi1) then
             qxm(i+1,j,QRHO) = rho_ref + apleft + amleft + azrleft
             qxm(i+1,j,QRHO) = max(qxm(i+1,j,QRHO),small_dens)
             qxm(i+1,j,QU) = u_ref + (apleft - amleft)*cc/rho
             qxm(i+1,j,QV) = v_ref + azv1left
             qxm(i+1,j,QPRES) = p_ref + (apleft + amleft)*csq
             qxm(i+1,j,QPRES) = max(qxm(i+1,j,QPRES), small_pres)
             qxm(i+1,j,QREINT) = rhoe_ref + (apleft + amleft)*enth*csq + azeleft
          end if

          ! geometry source terms
          if(dloga(i,j).ne.0)then
             courn = dtdx*(cc+abs(u))
             eta = (ONE-courn)/(cc*dt*abs(dloga(i,j)))
             dlogatmp = min(eta,ONE)*dloga(i,j)
             sourcr = -HALF*dt*rho*dlogatmp*u
             sourcp = sourcr*csq
             source = sourcp*enth
             if (i .le. ihi1) then
                qxm(i+1,j,QRHO) = qxm(i+1,j,QRHO) + sourcr
                qxm(i+1,j,QRHO) = max(qxm(i+1,j,QRHO),small_dens)
                qxm(i+1,j,QPRES) = qxm(i+1,j,QPRES) + sourcp
                qxm(i+1,j,QREINT) = qxm(i+1,j,QREINT) + source
             end if
             if (i .ge. ilo1) then
                qxp(i,j,QRHO) = qxp(i,j,QRHO) + sourcr
                qxp(i,j,QRHO) = max(qxp(i,j,QRHO),small_dens)
                qxp(i,j,QPRES) = qxp(i,j,QPRES) + sourcp
                qxp(i,j,QREINT) = qxp(i,j,QREINT) + source
             end if
          endif

       enddo
    enddo

    ! We do all passively advected quantities in one loop
    do ipassive = 1, npassive
       n = qpass_map(ipassive)
       do j = ilo2-1, ihi2+1

          ! Right state
          do i = ilo1, ihi1+1
             u = q(i,j,QU)
             if (u .gt. ZERO) then
                spzero = -ONE
             else
                spzero = u*dtdx
             endif
             acmprght = HALF*(-ONE - spzero )*dqx(i,j,n)
             qxp(i,j,n) = q(i,j,n) + acmprght
          enddo

          ! Left state
          do i = ilo1-1, ihi1
             u = q(i,j,QU)
             if (u .ge. ZERO) then
                spzero = u*dtdx
             else
                spzero = ONE
             endif
             acmpleft = HALF*(ONE - spzero )*dqx(i,j,n)
             qxm(i+1,j,n) = q(i,j,n) + acmpleft
          enddo

       enddo
    enddo

    !-------------------------------------------------------------------------
    ! y-direction
    !-------------------------------------------------------------------------

    ! Compute left and right traced states
    do j = ilo2-1, ihi2+1
       do i = ilo1-1, ihi1+1

          cc = c(i,j)
          csq = cc**2
          rho = q(i,j,QRHO)
          u = q(i,j,QU)
          v = q(i,j,QV)
          p = q(i,j,QPRES)
          rhoe = q(i,j,QREINT)
          enth = ( (rhoe+p)/rho )/csq

          drho = dqy(i,j,QRHO)
          du = dqy(i,j,QU)
          dv = dqy(i,j,QV)
          dp = dqy(i,j,QPRES)
          drhoe = dqy(i,j,QREINT)

          alpham = HALF*(dp/(rho*cc) - dv)*rho/cc
          alphap = HALF*(dp/(rho*cc) + dv)*rho/cc
          alpha0r = drho - dp/csq
          alpha0e = drhoe - dp*enth
          alpha0u = du

          e(1) = v-cc
          e(2) = v
          e(3) = v+cc

          ! construct the right state on the j-1/2 interface
          rho_ref = rho - HALF*(ONE + dtdy*min(e(1),ZERO))*drho
          u_ref = u - HALF*(ONE + dtdy*min(e(1),ZERO))*du
          v_ref = v - HALF*(ONE + dtdy*min(e(1),ZERO))*dv
          p_ref = p - HALF*(ONE + dtdy*min(e(1),ZERO))*dp
          rhoe_ref = rhoe - HALF*(ONE + dtdy*min(e(1),ZERO))*drhoe

          apright = 0.25d0*dtdy*(e(1) - e(3))*(ONE - sign(ONE,e(3)))*alphap
          amright = 0.25d0*dtdy*(e(1) - e(1))*(ONE - sign(ONE,e(1)))*alpham

          azrright = 0.25e0*dtdy*(e(1)-e(2))*(ONE - sign(ONE,e(2)))*alpha0r
          azeright = 0.25e0*dtdy*(e(1)-e(2))*(ONE - sign(ONE,e(2)))*alpha0e
          azu1rght = 0.25e0*dtdy*(e(1)-e(2))*(ONE - sign(ONE,e(2)))*alpha0u

          if (j .ge. ilo2) then
             qyp(i,j,QRHO) = rho_ref + apright + amright + azrright
             qyp(i,j,QRHO) = max(small_dens, qyp(i,j,QRHO))
             qyp(i,j,QV) = v_ref + (apright - amright)*cc/rho
             qyp(i,j,QU) = u_ref + azu1rght
             qyp(i,j,QPRES) = p_ref + (apright + amright)*csq
             qyp(i,j,QPRES) = max(qyp(i,j,QPRES), small_pres)
             qyp(i,j,QREINT) = rhoe_ref + (apright + amright)*enth*csq + azeright
          end if

          ! construct the left state on the j+1/2 interface

          rho_ref = rho + HALF*(ONE - dtdy*max(e(3),ZERO))*drho
          u_ref = u + HALF*(ONE - dtdy*max(e(3),ZERO))*du
          v_ref = v + HALF*(ONE - dtdy*max(e(3),ZERO))*dv
          p_ref = p + HALF*(ONE - dtdy*max(e(3),ZERO))*dp
          rhoe_ref = rhoe + HALF*(ONE - dtdy*max(e(3),ZERO))*drhoe

          apleft = 0.25d0*dtdy*(e(3) - e(3))*(ONE + sign(ONE,e(3)))*alphap
          amleft = 0.25d0*dtdy*(e(3) - e(1))*(ONE + sign(ONE,e(1)))*alpham

          azrleft = 0.25d0*dtdy*(e(3) - e(2))*(ONE + sign(ONE,e(2)))*alpha0r
          azeleft = 0.25d0*dtdy*(e(3) - e(2))*(ONE + sign(ONE,e(2)))*alpha0e
          azu1left = 0.25d0*dtdy*(e(3) - e(2))*(ONE + sign(ONE,e(2)))*alpha0u

          if (j .le. ihi2) then
             qym(i,j+1,QRHO) = rho_ref + apleft + amleft + azrleft
             qym(i,j+1,QRHO) = max(small_dens, qym(i,j+1,QRHO))
             qym(i,j+1,QV) = v_ref + (apleft - amleft)*cc/rho
             qym(i,j+1,QU) = u_ref + azu1left
             qym(i,j+1,QPRES) = p_ref + (apleft + amleft)*csq
             qym(i,j+1,QPRES) = max(qym(i,j+1,QPRES), small_pres)
             qym(i,j+1,QREINT) = rhoe_ref + (apleft + amleft)*enth*csq + azeleft
          end if

       enddo
    enddo

    do ipassive = 1, npassive
       n = qpass_map(ipassive)
       do i = ilo1-1, ihi1+1

          ! Top state
          do j = ilo2, ihi2+1
             v = q(i,j,QV)
             if (v .gt. ZERO) then
                spzero = -ONE
             else
                spzero = v*dtdy
             endif
             acmptop = HALF*(-ONE - spzero )*dqy(i,j,n)
             qyp(i,j,n) = q(i,j,n) + acmptop
          enddo

          ! Bottom state
          do j = ilo2-1, ihi2
             v = q(i,j,QV)
             if (v .ge. ZERO) then
                spzero = v*dtdy
             else
                spzero = ONE
             endif
             acmpbot = HALF*(ONE - spzero )*dqy(i,j,n)
             qym(i,j+1,n) = q(i,j,n) + acmpbot
          enddo

       enddo
    enddo

    deallocate(dqx,dqy)

  end subroutine trace

end module trace_module

max_stars_repo_name: ndeak/PeleC

max_issues_repo_name: ndeak/PeleC

max_forks_repo_name: ndeak/PeleC


text:
!****************************************************************
!*                                                            *
!*   module MO_SALSA_PROPERTIES                                 *
!*                                                              *
!*   Contains subroutines and functions that are used           *
!*   to calculate particle properties during simulation         *
!*                                                              *
!****************************************************************
MODULE mo_salsa_properties

CONTAINS

  ! fxm: should sea salt form a solid particle when prh is very low
  !  (even though it could be mixed with e.g. sulphate)?
  ! fxm: crashes if no sulphate or sea salt
  ! fxm: do we really need to consider Kelvin effect for regime 2
  !********************************************************************
  !
  ! subroutine WETSIZE()
  !
  !********************************************************************
  !
  ! Purpose:
  ! --------
  ! Calculates ambient sizes of particles by equilibrating
  !  soluble fraction of particles with water
  !
  !
  ! Method:
  ! -------
  ! Following chemical components are assumed water-soluble
  ! - (ammonium) sulphate (100%)
  ! - sea salt (100 %)
  ! - organic carbon (epsoc * 100%)
  !
  ! Exact thermodynamic considerations neglected
  ! - If particles contain no sea salt, calculation according to
  !  sulphate properties
  ! - If contain sea salt but no sulphate, calculation according to
  !  sea salt properties
  ! - If contain both sulphate and sea salt
  !  -> the molar fraction of these compounds determines
  !     which one of them is used as the basis of calculation
  !
  ! If sulphate and sea salt coexist in a particle,
  !   it is assumed that the Cl is replaced by sulphate;
  !   thus only either sulphate + organics or sea salt + organics
  !   is included in the calculation of soluble fraction.
  !
  ! Molality parameterizations taken from table 1 of
  !  Tang: Mixed-salt aerosols of atmospheric importance,
  !   JGR, 102 (D2), 1883-1893 (1997)
  !
  !
  ! Interface:
  ! ----------
  ! Called from main aerosol model
  !
  !
  ! Coded by:
  ! ---------
  ! Hannele Korhonen (FMI) 2005
  ! Harri Kokkola (FMI) 2006
  ! Matti Niskanen(FMI) 2012
  ! Anton Laakso  (FMI) 2013
  !
  ! Modified for the new aerosol datatype,
  ! Juha Tonttila (FMI) 2014
  !---------------------------------------------------------------------

  SUBROUTINE equilibration(kproma, kbdim, klev,    &
                           prh, ptemp, paero, init )

    USE mo_submctl, ONLY : &
         t_section,    &
         pi6,          & ! pi/6
         in1a, fn1a,   &
         in2a,    &
         fn2b,   &
         boltz,        & ! Boltzmann constant [J/K]
         nlim,         & ! lowest possible particle conc. in a bin [#/m3]

                         ! molar masses [kg/mol]
         mwa,          & ! water
                         ! molecular volumes [m3]
         mvsu,         & ! sulphate
                         ! density [kg/m3]
         rhowa,        & ! water

         surfw0,       & ! surface tension of water [J/m2]
         epsoc,        & ! fxm

         rhosu, msu,   & ! properties of compounds
         rhooc, moc,   &
         rhoss, mss,   &
         rhono, mno

    IMPLICIT NONE

    !-- input variables -------------
    INTEGER, INTENT(in) ::          &
         kproma,                    & ! number of horiz. grid kproma
         kbdim,                     & ! dimension for arrays
         klev                         ! number of vertical levels

    REAL, INTENT(in) ::        &
         prh(kbdim,klev),          & ! relative humidity [0-1]
         ptemp(kbdim,klev)           ! temperature [K]

    LOGICAL, INTENT(in) :: init  ! TRUE: Initialization call
                                 ! FALSE: Normal runtime: update water content only for 1a

    !-- output variables -------------
    TYPE(t_section), INTENT(inout) :: paero(kbdim,klev,fn2b)

    !-- local variables --------------
    INTEGER :: ii, jj, kk             ! loop indices
    INTEGER :: count

    REAL ::      &
         zbinmol(7), &   ! binary molality of individual components [mol/kg]
         zvpart(7),  &   ! volume of chem. compounds in one particle [fxm]
         zke,        &   ! Kelvin term
         zaw,        &   ! water activity [0-1]
         zlwc,       &   ! liquid water content [kg/m3-air]
         zdold,      &   !
         zrh
    REAL :: zcore,  &
                zdwet

    zcore = 0.
    zvpart(:) = 0.
    zlwc = 0.
    zdwet = 0.

    !----------------------------------------------------------------------
    !-- 1) Regime 1: sulphate and partly water-soluble OC -----------------
    !                This is done for every CALL
    zke = 1.001
    DO kk = in1a,fn1a      ! size bin
       DO jj = 1,klev      ! vertical grid
          DO ii = 1,kbdim ! horizontal grid

             !-- initialize
             zbinmol = 0.
             zdold = 1.

             IF ((paero(ii,jj,kk)%numc > nlim)) THEN

                !-- volume of sulphate and OC in one particle [fxm]
                !   + NO/NH

                zvpart(1:2) = paero(ii,jj,kk)%volc(1:2)/paero(ii,jj,kk)%numc
                zvpart(6:7) = paero(ii,jj,kk)%volc(6:7)/paero(ii,jj,kk)%numc ! NO + NH

                !-- total volume of one dry particle [fxm]
                zcore = sum(zvpart(1:2))
                zdwet = paero(ii,jj,kk)%dwet

                ! Relative Humidity:
                zrh = prh(ii,jj)
                zrh = MAX(zrh , 0.05)
                zrh = MIN(zrh , 0.98)

                count = 0
                DO WHILE(abs(zdwet/zdold-1.) > 1.e-2)
                   zdold = max(zdwet,1.e-20)

                   zaw = zrh/zke

                   !-- binary molalities [mol/kg]
                   zbinmol(1) =                      & ! sulphate
                        + 1.1065495e+2            &
                        - 3.6759197e+2 * zaw      &
                        + 5.0462934e+2 * zaw**2   &
                        - 3.1543839e+2 * zaw**3   &
                        + 6.770824e+1  * zaw**4

                   zbinmol(2) = 1./(zaw*mwa)-1./mwa ! organic carbon

                   zbinmol(6) =                      & ! nitric acid
                        + 2.306844303e+1          &
                        - 3.563608869e+1 * zaw    &
                        - 6.210577919e+1 * zaw**2 &
                        + 5.510176187e+2 * zaw**3 &
                        - 1.460055286e+3 * zaw**4 &
                        + 1.894467542e+3 * zaw**5 &
                        - 1.220611402e+3 * zaw**6 &
                        + 3.098597737e+2 * zaw**7

                   !
                   ! Calculate the liquid water content (kg/m3-air) using ZSR
                   ! (see e.g. equation (9.98) in Seinfeld and Pandis (1998))
                   !
                   zlwc = (paero(ii,jj,kk)%volc(1)*(rhosu/msu))/zbinmol(1)       + &
                          epsoc * paero(ii,jj,kk)%volc(2)*(rhooc/moc)/zbinmol(2) + &
                          (paero(ii,jj,kk)%volc(6)*(rhono/mno))/zbinmol(6)

                   !-- particle wet radius [m]
                   zdwet = (zlwc/paero(ii,jj,kk)%numc/rhowa/pi6 + &
                        (SUM(zvpart(6:7))/pi6)**(1./3.) + &
                        zcore/pi6)**(1./3.)

                   zke = exp(2.*surfw0*mvsu/(boltz*ptemp(ii,jj)*zdwet))
                   !-- Kelvin effect

                   count = count + 1
                   IF (count > 1000) STOP 'SALSA equilibration (regime 1): no convergence!!'

                END DO

                ! Instead of lwc, use the volume concentration of water from now on
                ! (easy to convert...)
                paero(ii,jj,kk)%volc(8) = zlwc/rhowa

                ! If this is initialization, update the core and wet diameter
                IF (init) THEN
                   paero(ii,jj,kk)%dwet = zdwet
                   paero(ii,jj,kk)%core = zcore
                END IF

             ELSE
                ! If initialization
                !-- 1.2) empty bins given bin average values -----------------
                IF (init) THEN
                   paero(ii,jj,kk)%dwet = paero(ii,jj,kk)%dmid
                   paero(ii,jj,kk)%core = pi6*paero(ii,jj,kk)%dmid**3
                END IF
             END IF
          END DO
       END DO
    END DO

    !-- 2) Regime 2a: sulphate, OC, BC and sea salt ----------------------------
    !                 This is done only for initialization call, otherwise the
    !                 water contents are computed via condensation

    IF (init) THEN

       ! loops over:
       DO kk = in2a,fn2b      ! size bin
          DO jj = 1,klev      ! vertical grid
             DO ii = 1,kbdim ! horizontal grid

                zke = 1.02

                !-- initialize
                zbinmol = 0.
                zdold = 1.

                !-- 1) particle properties calculated for non-empty bins ---------
                IF ((paero(ii,jj,kk)%numc > nlim)) THEN

                   !-- volume in one particle [fxm]
                   zvpart = 0.
                   zvpart(1:7) = paero(ii,jj,kk)%volc(1:7)/paero(ii,jj,kk)%numc

                   !-- total volume of one dry particle [fxm]
                   zcore = sum(zvpart(1:5))
                   zdwet = paero(ii,jj,kk)%dwet

                   ! Relative Humidity:
                   zrh = prh(ii,jj)
                   zrh = MAX(zrh , 0.37)
                   zrh = MIN(zrh , 0.98)

                   count = 0
                   DO WHILE(abs(zdwet/zdold-1.) > 1.e-12)
                      zdold = max(zdwet,1.e-20)

                      zaw = zrh/zke

                      !-- binary molalities [mol/kg]
                      zbinmol(1) =                     &  ! sulphate
                           + 1.1065495e+2           &
                           - 3.6759197e+2 * zaw     &
                           + 5.0462934e+2 * zaw**2  &
                           - 3.1543839e+2 * zaw**3  &
                           + 6.770824e+1  * zaw**4

                      zbinmol(2) = 1./(zaw*mwa)-1./mwa ! organic carbon

                      zbinmol(6) =                      & ! nitric acid
                           + 2.306844303e+1          &
                           - 3.563608869e+1 * zaw    &
                           - 6.210577919e+1 * zaw**2 &
                           + 5.510176187e+2 * zaw**3 &
                           - 1.460055286e+3 * zaw**4 &
                           + 1.894467542e+3 * zaw**5 &
                           - 1.220611402e+3 * zaw**6 &
                           + 3.098597737e+2 * zaw**7

                      zbinmol(5) =                     &  ! sea salt (NaCl)
                           + 5.875248e+1            &  !
                           - 1.8781997e+2 * zaw     &
                           + 2.7211377e+2 * zaw**2  &
                           - 1.8458287e+2 * zaw**3  &
                           + 4.153689e+1  * zaw**4

                      !-- calculate the liquid water content (kg/m3-air)
                      zlwc = (paero(ii,jj,kk)%volc(1)*(rhosu/msu))/zbinmol(1) +                 &
                           epsoc * (paero(ii,jj,kk)%volc(2)*(rhooc/moc))/zbinmol(2) +         &
                           (paero(ii,jj,kk)%volc(6)*(rhono/mno))/zbinmol(6)         +         &
                           (paero(ii,jj,kk)%volc(5)*(rhoss/mss))/zbinmol(5)

                      !-- particle wet radius [m]
                      zdwet = (zlwc/paero(ii,jj,kk)%numc/rhowa/pi6 +  &
                           (SUM(zvpart(6:7))/pi6)**(1./3.) + &
                           zcore/pi6)**(1./3.)

                      !-- Kelvin effect
                      zke = exp(2.*surfw0*mvsu/(boltz*ptemp(ii,jj)*zdwet))

                      count = count + 1
                      IF (count > 1000) STOP 'SALSA equilibration (regime 2): no convergence!!'

                   END DO

                   ! Liquid water content; instead of LWC use the volume concentration
                   !plwc(ii,jj,kk)=zlwc
                   paero(ii,jj,kk)%volc(8) = zlwc/rhowa
                   paero(ii,jj,kk)%dwet = zdwet
                   paero(ii,jj,kk)%core = zcore

                ELSE
                   !-- 2.2) empty bins given bin average values -------------------------
                   paero(ii,jj,kk)%dwet = paero(ii,jj,kk)%dmid
                   paero(ii,jj,kk)%core = pi6*paero(ii,jj,kk)%dmid**3
                END IF

             END DO
          END DO
       END DO

    END IF

  END SUBROUTINE equilibration

  ! Juha: It should not be necessary to do this since cloud water content is Always calculated via condensation equations
  !               - This is done for initialization call when non-zero cloud or ice number concentration is initialized
  !                 manually. Otherwise water contents are computed via condensation
  !               - Not an equilibrium, but fixed droplet/ice diameter
  !               - Regimes 2a & 2b for ice and cloud droplets
  ! ----------------------------------------------------------------------------------------------------------------------
  SUBROUTINE equilibration_cloud(kproma, kbdim, klev,    &
                           prh, ptemp, pcloud, pice )

    USE mo_submctl, ONLY : &
         t_section,    &
         pi6,          & ! pi/6
         ica, fcb,     &
         iia, fib,     &
         ncld,nice,  &
         nlim, prlim

    IMPLICIT NONE

    !-- input variables -------------
    INTEGER, INTENT(in) ::          &
         kproma,                    & ! number of horiz. grid kproma
         kbdim,                     & ! dimension for arrays
         klev                         ! number of vertical levels

    REAL, INTENT(in) ::        &
         prh(kbdim,klev),          & ! relative humidity [0-1]
         ptemp(kbdim,klev)           ! temperature [K]

    !-- output variables -------------
    TYPE(t_section), INTENT(inout) :: pcloud(kbdim,klev,ncld), &
                                      pice(kbdim,klev,nice)

    !-- local variables --------------
    INTEGER :: ii, jj, kk             ! loop indices

    DO jj = 1,klev      ! vertical grid
        DO ii = 1,kbdim ! horizontal grid
            ! Cloud droplets
            DO kk = ica%cur,fcb%cur      ! size bin
                IF ((pcloud(ii,jj,kk)%numc > nlim)) THEN
                   !-- 1) particle properties calculated for non-empty bins ---------
                   pcloud(ii,jj,kk)%volc(8) = pi6*12.0E-6**3*pcloud(ii,jj,kk)%numc ! set water content according to 12um diameter
                   pcloud(ii,jj,kk)%dwet = 12.0E-6
                END IF
            ENDDO

            ! Ice
            DO kk = iia%cur,fib%cur      ! size bin
                IF ((pice(ii,jj,kk)%numc > prlim)) THEN
                    !-- 1) particle properties calculated for non-empty bins ---------
                   pice(ii,jj,kk)%volc(8) = pi6*30.0E-6**3*pice(ii,jj,kk)%numc  ! set water content according to 30um diameter
                   pice(ii,jj,kk)%dwet = 30.0E-6! sqrt(3*B/A)
                END IF
            END DO
        END DO
    END DO

  END SUBROUTINE equilibration_cloud

END MODULE mo_salsa_properties

max_stars_repo_name: huangynj/UCLALES-SALSA

max_issues_repo_name: huangynj/UCLALES-SALSA

max_forks_repo_name: huangynj/UCLALES-SALSA


text:
module test_stdlib_bitset_64
    use testdrive, only : new_unittest, unittest_type, error_type, check
    use :: stdlib_kinds, only : int8, int16, int32, int64
    use stdlib_bitsets
    implicit none
    private

    public :: collect_stdlib_bitset_64

    character(*), parameter :: &
        bitstring_0   = '000000000000000000000000000000000', &
        bitstring_33  = '100000000000000000000000000000000', &
        bitstring_all = '111111111111111111111111111111111'

contains

    !> Collect all exported unit tests
    subroutine collect_stdlib_bitset_64(testsuite)
        !> Collection of tests
        type(unittest_type), allocatable, intent(out) :: testsuite(:)

        testsuite = [ &
            new_unittest("string-operations-0", test_string_operations_0), &
            new_unittest("string-operations-1", test_string_operations_1), &
            new_unittest("string-operations-3", test_string_operations_3), &
            new_unittest("string-operations-4", test_string_operations_4), &
            new_unittest("io", test_io), &
            new_unittest("initialization", test_initialization), &
            new_unittest("bitset-inquiry", test_bitset_inquiry), &
            new_unittest("bit-operations", test_bit_operations), &
            new_unittest("bitset-comparisons", test_bitset_comparisons), &
            new_unittest("bitset-operations-and", test_bitset_operations_and), &
            new_unittest("bitset-operations-nand", test_bitset_operations_nand), &
            new_unittest("bitset-operations-or", test_bitset_operations_or), &
            new_unittest("bitset-operations-xor", test_bitset_operations_xor) &
            ]

    end subroutine collect_stdlib_bitset_64

    subroutine test_string_operations_0(error)
        !> Error handling
        type(error_type), allocatable, intent(out) :: error

        type(bitset_64) :: set
        character(:), allocatable :: string0

        call set%from_string(bitstring_0)

        call check(error, bits(set), 33)
        if (allocated(error)) return

        call check(error, set%none())
        if (allocated(error)) return

        call check(error, .not.set%any())
        if (allocated(error)) return

        call set%to_string(string0)

        call check(error, string0, bitstring_0)
        if (allocated(error)) return

        call set%write_bitset(string0)

        call check(error, string0, ('S33B' // bitstring_0))
        if (allocated(error)) return

    end subroutine test_string_operations_0

    subroutine test_string_operations_1(error)
        !> Error handling
        type(error_type), allocatable, intent(out) :: error

        type(bitset_64) :: set
        character(:), allocatable :: string0

        call set%from_string(bitstring_all)

        call check(error, bits(set), 33)
        if (allocated(error)) return

        call check(error, .not.set%none() )
        if (allocated(error)) return

        call check(error, set%any() )
        if (allocated(error)) return

        call check(error, set%all() )
        if (allocated(error)) return

        call set%to_string(string0)

        call check(error, string0, bitstring_all)
        if (allocated(error)) return

        call set%write_bitset(string0)

        call check(error, string0, ('S33B' // bitstring_all))
        if (allocated(error)) return

    end subroutine test_string_operations_1

    subroutine test_string_operations_3(error)
        !> Error handling
        type(error_type), allocatable, intent(out) :: error

        type(bitset_64) :: set
        integer         :: status

        call set%read_bitset(bitstring_0, status)

        call check(error, status /= success)
        if (allocated(error)) return

        call set%read_bitset('s33b' // bitstring_0, status)

        call check(error, bits(set), 33)
        if (allocated(error)) return

        call check(error, set%none())
        if (allocated(error)) return

    end subroutine test_string_operations_3

    subroutine test_string_operations_4(error)
        !> Error handling
        type(error_type), allocatable, intent(out) :: error

        type(bitset_64) :: set

        call set%read_bitset('s33b' // bitstring_all )

        call check(error, bits(set), 33)
        if (allocated(error)) return

        call check(error, .not.set%none())
        if (allocated(error)) return

        call check(error, set%any())
        if (allocated(error)) return

        call check(error, set%all())
        if (allocated(error)) return

    end subroutine test_string_operations_4

    subroutine test_io(error)
        !> Error handling
        type(error_type), allocatable, intent(out) :: error

        integer :: unit
        type(bitset_64) :: set0, set1, set2, set3, set4, set5

        call set0%from_string(bitstring_0)
        call set1%from_string(bitstring_all)
        call set2%from_string(bitstring_33)

        open( newunit=unit, status='scratch', form='formatted', &
              action='readwrite' )
        call set2%write_bitset(unit)
        call set1%write_bitset(unit)
        call set0%write_bitset(unit)
        rewind( unit )
        call set3%read_bitset(unit)
        call set5%read_bitset(unit)
        call set4%read_bitset(unit)

        call check(error, set4 == set0 .and. set5 == set1 .and. set3 == set2, &
            'transfer to and from units using bitset literals failed.')
        if (.not.allocated(error)) then
            rewind( unit )

            call set2%write_bitset(unit, advance='no')
            call set1%write_bitset(unit, advance='no')
            call set0%write_bitset(unit)
            rewind( unit )
            call set3%read_bitset(unit, advance='no')
            call set4%read_bitset(unit, advance='no')
            call set5%read_bitset(unit)

            call check(error, set5 == set0 .and. set4 == set1 .and. set3 == set2, &
                'transfer to and from units using bitset literals with advance="no" failed.')
        end if
        close(unit)
        if (allocated(error)) return

        open( newunit=unit, form='unformatted', status='scratch', &
              action='readwrite' )
        call set2%output(unit)
        call set1%output(unit)
        call set0%output(unit)
        rewind( unit )
        call set5%input(unit)
        call set4%input(unit)
        call set3%input(unit)
        close( unit )

        call check(error, set3 == set0 .and. set4 == set1 .and. set5 == set2, &
            'transfer to and from units using output and input failed.')
        if (allocated(error)) return

        open( newunit=unit, form='unformatted', access='stream', &
              status='scratch', action='readwrite' )
        call set2%output(unit)
        call set1%output(unit)
        call set0%output(unit)
        rewind( unit )
        call set5%input(unit)
        call set4%input(unit)
        call set3%input(unit)
        close( unit )

        call check(error, set3 == set0 .and. set4 == set1 .and. set5 == set2, &
            'transfer to and from units using stream output and input failed.')
        if (allocated(error)) return

    end subroutine test_io

    subroutine test_initialization(error)
        !> Error handling
        type(error_type), allocatable, intent(out) :: error

        logical(int8)  :: log1(64) = .true.
        logical(int16) :: log2(31) = .false.
        logical(int32) :: log3(15) = .true.
        logical(int64) :: log4(33) = .false.
        logical(int8), allocatable  :: log5(:)
        logical(int16), allocatable :: log6(:)
        logical(int32), allocatable :: log7(:)
        logical(int64), allocatable :: log8(:)

        type(bitset_64) :: set4, set5

        set5 = log1
        call check(error, set5%bits(), 64, &
            'initialization with logical(int8) failed to set the right size.')
        if (allocated(error)) return
        call check(error, set5%all(), &
            'initialization with logical(int8) failed to set the right values.')
        if (allocated(error)) return

        set5 = log2
        call check(error, set5%bits(), 31, &
            'initialization with logical(int16) failed to set the right size.')
        if (allocated(error)) return
        call check(error, set5%none(), &
            'initialization with logical(int16) failed to set the right values.')
        if (allocated(error)) return

        set5 = log3
        call check(error, set5%bits(), 15, &
            'initialization with logical(int32) failed to set the right size.')
        if (allocated(error)) return
        call check(error, set5%all(), &
            'initialization with logical(int32) failed to set the right values.')
        if (allocated(error)) return

        set5 = log4
        call check(error, set5%bits(), 33, &
            'initialization with logical(int64) failed to set the right size.')
        if (allocated(error)) return
        call check(error, set5%none(), &
            'initialization with logical(int64) failed to set the right values.')
        if (allocated(error)) return

        set5 = log1
        call extract( set4, set5, 1_bits_kind, 33_bits_kind )
        call check(error, set4%bits(), 33, &
            'initialization with extract failed to set the right size.')
        if (allocated(error)) return
        call check(error, set4%all(), &
            'initialization with extract failed to set the right values.')
        if (allocated(error)) return

        set4 = set5
        call check(error, set4%bits(), 64, &
            'initialization with simple assignment failed to set the right size.')
        if (allocated(error)) return
        call check(error, set4%all(), &
            'initialization with simple assignment failed to set the right values.')
        if (allocated(error)) return

        log5 = set5
        call check(error, size(log5), 64, &
            'initialization of logical(int8) with assignment failed to set the right size.')
        if (allocated(error)) return
        call check(error, all(log5) .eqv. .true., &  ! FIXME
            'initialization of logical(int8) with assignment failed to set the right values.')
        if (allocated(error)) return

        log6 = set5
        call check(error, size(log6), 64, &
            'initialization of logical(int16) with assignment failed to set the right size.')
        if (allocated(error)) return
        call check(error, all(log6) .eqv. .true., &  ! FIXME
            'initialization of logical(int16) with assignment failed to set the right values.')
        if (allocated(error)) return

        log7 = set5
        call check(error, size(log7), 64, &
            'initialization of logical(int32) with assignment failed to set the right size.')
        if (allocated(error)) return
        call check(error, all(log7), &
            'initialization of logical(int32) with assignment failed to set the right values.')
        if (allocated(error)) return

        log8 = set5
        call check(error, size(log8), 64, &
            'initialization of logical(int64) with assignment failed to set the right size.')
        if (allocated(error)) return
        call check(error, merge(.true., .false., all(log8)), &  ! FIXME
            'initialization of logical(int64) with assignment failed to set the right values.')
        if (allocated(error)) return

    end subroutine test_initialization

    subroutine test_bitset_inquiry(error)
        !> Error handling
        type(error_type), allocatable, intent(out) :: error

        integer(bits_kind) :: i

        type(bitset_64) :: set0, set1

        call set0%from_string(bitstring_0)
        call set1%from_string(bitstring_all)

        call check(error, set0%none(), 'set0 did not have none set which ' // &
            'was unexpected')
        if (allocated(error)) return
        call check(error, .not. set0%any(), 'set0 had some bits set which ' // &
            'was unexpected.')
        if (allocated(error)) return

        call set0%not()
        call check(error, set0%all(), 'set0 did not have all bits set ' // &
            'which was unexpected')
        if (allocated(error)) return
        call check(error, set0%any(), 'set0 had no bits set which ' // &
            'was unexpected.')
        if (allocated(error)) return

        call check(error, set1%any(), 'set1 had no bits set ' // &
            'which was unexpected')
        if (allocated(error)) return
        call check(error, set1%all(), 'set1 did not have all bits set ' // &
            'which was unexpected.')
        if (allocated(error)) return

        call set0%not()
        do i=0, set0%bits() - 1
            call check(error, .not. set0%test(i), &
                'against expectations set0 has at least 1 bit set.')
        end do

        do i=0, set1%bits() - 1
            call check(error, set1%test(i), &
                'against expectations set1 has at least 1 bit unset.')
        end do

        do i=0, set0%bits() - 1
            call check(error, .not.( set0%value(i) /= 0), &
                'against expectations set0 has at least 1 bit set.')
        end do

        do i=0, set1%bits() - 1
            call check(error, .not.( set1%value(i) /= 1), &
                'against expectations set1 has at least 1 bit unset.')
        end do

        call check(error, set0%bits() == 33, 'et0 unexpectedly does not have 33 bits.')
        if (allocated(error)) return

    end subroutine test_bitset_inquiry

    subroutine test_bit_operations(error)
        !> Error handling
        type(error_type), allocatable, intent(out) :: error

        type(bitset_64) :: set1

        call set1%from_string(bitstring_all)

        call check(error, set1%all(), 'set1 is not all set.')
        if (allocated(error)) return

        call set1%clear(0_bits_kind)
        call check(error, .not. set1%test(0_bits_kind), 'did not clear the first bit in set1.')
        if (allocated(error)) return
        call check(error, set1%test(1_bits_kind), 'cleared more than one bit in set1.')
        if (allocated(error)) return

        call set1%clear(1_bits_kind, 32_bits_kind)
        call check(error, set1%none(), 'did not clear remaining bits in set1.')
        if (allocated(error)) return

        call set1%flip(0_bits_kind)
        call check(error, set1%test(0_bits_kind), 'did not flip the first bit in set1.')
        if (allocated(error)) return
        call check(error, .not. set1%test(1_bits_kind), 'flipped more than one bit in set1.')
        if (allocated(error)) return

        call set1%flip(1_bits_kind, 32_bits_kind)
        call check(error, set1%all(), 'did not flip remaining bits in set1.')
        if (allocated(error)) return

        call set1%not()
        call check(error, set1%none(), 'did not unset bits in set1.')
        if (allocated(error)) return

        call set1%set(0_bits_kind)
        call check(error, set1%test(0_bits_kind), 'did not set the first bit in set1.')
        if (allocated(error)) return
        call check(error, .not. set1%test(1_bits_kind), 'set more than one bit in set1.')
        if (allocated(error)) return

        call set1%set(1_bits_kind, 32_bits_kind)
        call check(error, set1%all(), 'did not set the remaining bits in set1.')
        if (allocated(error)) return

    end subroutine test_bit_operations

    subroutine test_bitset_comparisons(error)
        !> Error handling
        type(error_type), allocatable, intent(out) :: error

        type(bitset_64) :: set0, set1, set2

        call set0%from_string(bitstring_0)
        call set1%from_string(bitstring_all)
        call set2%from_string(bitstring_33)

        call check(error, set0 == set0 .and. set1 == set1 .and. set2 == set2 .and. &
            .not. set0 == set1 .and. .not. set0 == set2 .and. .not.   &
            set1 == set2, 'failed 64 bit equality tests.')
        if (allocated(error)) return

        call check(error, set0 /= set1 .and. set1 /= set2 .and. set0 /= set2 .and. &
            .not. set0 /= set0 .and. .not. set1 /= set1 .and. .not.   &
            set2 /= set2, 'failed 64 bit inequality tests.')
        if (allocated(error)) return

        call check(error, set1 > set0 .and. set2 > set0 .and. set1 > set2 .and. &
            .not. set0 > set1 .and. .not. set1 > set1 .and. .not.  &
            set2 > set1, 'failed 64 bit greater than tests.')
        if (allocated(error)) return

        call check(error, set1 >= set0 .and. set1 >= set2 .and. set2 >= set2 .and. &
            .not. set0 >= set1 .and. .not. set0 >= set1 .and. .not.  &
            set2 >= set1, 'failed 64 bit greater than or equal tests.')
        if (allocated(error)) return

        call check(error, set0 < set1 .and. set0 < set1 .and. set2 < set1 .and. &
            .not. set1 < set0 .and. .not. set0 < set0 .and. .not.  &
            set1 < set2, 'failed 64 bit less than tests.')
        if (allocated(error)) return

        call check(error, set0 <= set1 .and. set2 <= set1 .and. set2 <= set2 .and. &
            .not. set1 <= set0 .and. .not. set2 <= set0 .and. .not.  &
            set1 <= set2, 'failed 64 bit less than or equal tests.')
        if (allocated(error)) return

    end subroutine test_bitset_comparisons

    subroutine test_bitset_operations_and(error)
        !> Error handling
        type(error_type), allocatable, intent(out) :: error

        type(bitset_64) :: set3, set4, set0

        call set0%from_string( bitstring_all )
        call set4%from_string( bitstring_all )
        call and( set0, set4 ) ! all all
        call check(error, set0%all(), 'first test of AND failed.')
        if (allocated(error)) return

        call set4%from_string( bitstring_0 )
        call set3%from_string( bitstring_all )
        call and( set3, set4 ) ! all none
        call check(error, set3%none(), 'second test of AND failed.')
        if (allocated(error)) return

        call set3%from_string( bitstring_all )
        call set4%from_string( bitstring_0 )
        call and( set4, set3 ) ! none all
        call check(error, set4%none(), 'third test of AND failed.')
        if (allocated(error)) return

        call set3%from_string( bitstring_0 )
        call and( set4, set3 ) ! none none
        call check(error, set4%none(), 'fourth test of AND failed.')
        if (allocated(error)) return

    end subroutine test_bitset_operations_and

    subroutine test_bitset_operations_nand(error)
        !> Error handling
        type(error_type), allocatable, intent(out) :: error

        type(bitset_64) :: set3, set4

        call set3%from_string( bitstring_all )
        call set4%from_string( bitstring_all )
        call and_not( set4, set3 ) ! all all
        call check(error, set4%none(), 'first test of AND_NOT failed.')
        if (allocated(error)) return

        call set4%from_string( bitstring_0 )
        call and_not( set4, set3 ) ! none all
        call check(error, set4%none(), 'second test of AND_NOT failed.')
        if (allocated(error)) return

        call set3%from_string( bitstring_all )
        call set4%from_string( bitstring_0 )
        call and_not( set3, set4 ) ! all none
        call check(error, set3%all(), 'third  test of AND_NOT failed.')
        if (allocated(error)) return

        call set3%from_string( bitstring_0 )
        call set4%from_string( bitstring_0 )
        call and_not( set3, set4 ) ! none none
        call check(error, set3%none(), 'fourth  test of AND_NOT failed.')
        if (allocated(error)) return

    end subroutine test_bitset_operations_nand

    subroutine test_bitset_operations_or(error)
        !> Error handling
        type(error_type), allocatable, intent(out) :: error

        type(bitset_64) :: set3, set4

        call set3%from_string( bitstring_all )
        call set4%from_string( bitstring_all )
        call or( set3, set4 ) ! all all
        call check(error, set3%all(), 'first test of OR failed.')
        if (allocated(error)) return

        call set3%from_string( bitstring_0 )
        call or( set4, set3 ) ! all none
        call check(error, set4%all(), 'second test of OR failed.')
        if (allocated(error)) return

        call or( set3, set4 ) ! none all
        call check(error, set3%all(), 'third test of OR failed.')
        if (allocated(error)) return

        call set3%from_string( bitstring_0 )
        call set4%from_string( bitstring_0 )
        call or( set4, set3 ) !none none
        call check(error, set4%none(), 'fourth test of OR failed.')
        if (allocated(error)) return

    end subroutine test_bitset_operations_or

    subroutine test_bitset_operations_xor(error)
        !> Error handling
        type(error_type), allocatable, intent(out) :: error

        type(bitset_64) :: set3, set4

        call set3%from_string( bitstring_0 )
        call set4%from_string( bitstring_0 )
        call xor( set3, set4 ) ! none none
        call check(error, set3%none(), 'first test of XOR failed.')
        if (allocated(error)) return

        call set4%from_string( bitstring_all )
        call xor( set3, set4 ) ! none all
        call check(error, set3%all(), 'second test of XOR failed.')
        if (allocated(error)) return

        call set4%from_string( bitstring_0 )
        call xor( set3, set4 ) ! all none
        call check(error, set3%all(), 'third test of XOR failed.')
        if (allocated(error)) return

        call set4%from_string( bitstring_all )
        call xor( set3, set4 ) ! all all
        call check(error, set3%none(), 'fourth test of XOR failed.')
        if (allocated(error)) return

    end subroutine test_bitset_operations_xor

end module test_stdlib_bitset_64

program tester
    use, intrinsic :: iso_fortran_env, only : error_unit
    use testdrive, only : run_testsuite, new_testsuite, testsuite_type
    use test_stdlib_bitset_64, only : collect_stdlib_bitset_64
    implicit none
    integer :: stat, is
    type(testsuite_type), allocatable :: testsuites(:)
    character(len=*), parameter :: fmt = '("#", *(1x, a))'

    stat = 0

    testsuites = [ &
        new_testsuite("stdlib-bitset-64", collect_stdlib_bitset_64) &
        ]

    do is = 1, size(testsuites)
        write(error_unit, fmt) "Testing:", testsuites(is)%name
        call run_testsuite(testsuites(is)%collect, error_unit, stat)
    end do

    if (stat > 0) then
        write(error_unit, '(i0, 1x, a)') stat, "test(s) failed!"
        error stop
    end if
end program

max_stars_repo_name: freevryheid/stdlib

max_issues_repo_name: seanwallawalla-forks/stdlib

max_forks_repo_name: seanwallawalla-forks/stdlib


text:
PROGRAM bufr_readsst
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
!  Original program by Bert Katz X/X/2005
!  Modified 11 Dec 2008 by Robert Grumbine
!                                                                      C
!     IN THE SCRIPT THAT RUNS THIS PROGRAM, ONE MUST ISSUE THE         C
!     COMMANDS:                                                        C
!                                                                      C
!     export XLFRTEOPTS="unit_vars=yes"                                C
!     export XLFUNIT_11="NAME OF INPUT FILE (INCLUDING FULL PATH)"     C
!                                                                      C
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      !RG IMPLICIT NONE

      CHARACTER(80)  IDST1,IDST2,IDCHN

      DATA IDST1  /'YEAR MNTH DAYS HOUR MINU SECO SAID CLATH CLONH'/
      DATA IDST2  /'SAZA SOZA FOVN CLAVR '/
      DATA IDCHN  /'INCN ALBD TMBR '/

      REAL*8 XDATA(18),XCHAN(3,5)
      DOUBLE PRECISION dst(18), chan(3,5)
      CHARACTER(8)  SUBSET
      INTEGER outunit, creturns
      INTEGER lat
      PARAMETER (lat = 8)

      DATA IUNT /11/

      outunit = 51
      creturns = openout(outunit)

      CALL OPENBF(IUNT,'IN',IUNT)

      CALL READMG(IUNT,SUBSET,IDATE,IERR)
      IF(IERR.NE.0) GO TO 8888

      IREC = 0

10    CONTINUE

      CALL READSB(IUNT,IERR)
      IF(IERR.NE.0) THEN
        CALL READMG(IUNT,SUBSET,IDATE,IERR)
        IF(IERR.NE.0) GO TO 8888
      ELSE
        CALL UFBINT(IUNT,XDATA( 1),9,1,IRET,IDST1)
        CALL UFBINT(IUNT,XDATA(10),4,1,IRET,IDST2)
        CALL UFBREP(IUNT,XCHAN    ,3,5,IRET,IDCHN)
        IREC = IREC + 1
        do lll = 1, 9
          dst(lll)   = XDATA(lll)
          dst(lll+9) = XDATA(lll+9)
        enddo
        chan = XCHAN

        IF (MOD(IREC,100000) .EQ. 0) THEN
          PRINT *,'irec = ',IREC/100000,'e5'
        ENDIF
!Since the concern is polar, trim off the low-latitude data
        IF (ABS(XDATA(lat)) .GT. 30) THEN
          CALL gacout(dst, chan)
        ENDIF

      ENDIF
      GO TO 10

!  I/O ERROR EXIT
!  --------------

 8888 CONTINUE
      PRINT *,'IERR = ',IERR
      WRITE(6,1001) IREC
 1001 FORMAT(' EOF ENCOUNTERED IN SST INPUT FILE AFTER RECORD NO. ',I8)
      STOP

!  NORMAL EXIT
!  -----------

END

max_stars_repo_name: NOAA-EMC/seaice-concentration

max_issues_repo_name: NOAA-EMC/seaice-concentration

max_forks_repo_name: rgrumbine/concentration


#include "f_defs.h"

module realspace_common_m
  use global_m
  implicit none
  private
  public :: realspace_t

  type realspace_t
     character*20 :: filename         !< either epsmat.h5 or chimat.h5
     real(DP) :: epsinvhead           !< need to be set for filename = epsmat.h5
     integer :: nrq                   !< number of qpoints in RBZ
     real(DP), pointer :: rq(:,:)     !< qpoints in RBZ
     integer :: nfq                   !< number of qpoints in FBZ
     real(DP), pointer :: fq(:,:)     !< qpoints in FBZ
     integer, pointer :: nmtx(:)      !< number of Gvectors for each RBZ qpoint
     integer :: nmtx_max              !< nmtx_max = MAXVAL(nmtx(:))
     integer :: nmtx_used             !< nmtx_used <= MINVAL(nmtx(:))
     real(DP) :: ecut                 !< nmtx_used is determined from ecut
     integer :: qgrid(3)              !< qgrid
     integer :: ispin                 !< which spin component to plot
     integer :: nr2                   !< number of r2 points
     real(DP), allocatable :: r2(:,:) !< coordinate of r2 points, with respect to lattice vectors
     logical :: unfold                !< use symmetries to unfold grid
     integer :: downsample(3)         !< downsample each dimention of the real-space grid by this factor. Defaults to 2.
     integer :: nsuper(3)             !< number of unit-cell repetitions
     integer, pointer :: isrtx(:,:)   !< gvec%components(:,isrtx(ig)) is the ig-th Gvector in epsmat
     integer, pointer :: isrtxi(:,:)  !< ig_gvec = isrtx(ig), isrtxi(ig_gvec) = ig
     logical :: realpart              !< output real part of scfft
     logical :: imagpart              !< output imaginary part of scfft
     logical :: abs                   !< output abs(scfft)
     logical :: abs2                  !< output abs(scfft)**2

     !> Frequency dependence
     integer :: freq_dep        !< frequency dependence of the inverse dielectric matrix
                                ! 0: static calculation 2: full frequency 3: two imaginary frequencies
     integer :: nfreq           !< number of frequencies used in full frequency calculation
     integer :: nfreq_imag      !< number of imaginary freqs for CD (also 1 for GN GPP)
     complex(DP) :: freq_target
     integer :: FFTgrid(3), FFTfactor
     logical :: high_resolution
     logical :: low_comm, mid_comm
     integer :: rq_blocksize, nrq_loc, nrq_loc_max
  end type realspace_t

end module realspace_common_m


text:
      program t
      integer,dimension(4)::i=(/0,(j*2,j=1,2),6/)
      print *,i
      end program t

max_stars_repo_name: maddenp/ppp

max_issues_repo_name: maddenp/ppp

max_forks_repo_name: maddenp/ppp


!***********************************/
!	Name:Gnuplot用の出力結果を1つにまとめてGif作成用のファイルを作るプログラム
!	Alias:MergeGnuPlot
!	Description:
!	Type:Geom,iStep
!	Input:Density/Density*.txt,Velocity/Velocity*.txt,Pressure/Pressure*.txt
!	Output:
!	Note:これ単体のQMFというプログラムが存在する(Gnuplotスクリプトを吐き出すのでそっちの方が有能)
!	Author:Akitaka Toyota
!	Date:2017.10.26
!	Update:2017.11.09
!	Other:
!***********************************/
subroutine MergeGnuPlot(Geom,iStep)
    use StructVar_Mod
    use LoopVar_Mod
    implicit none

    type(Geometry),intent(in) :: Geom
    integer,intent(in) :: iStep

    integer :: iProgress
    double precision :: Tmp
    double precision, allocatable :: TmpQuantity(:)
    character(len=64), target :: cStep,cDensity,cVelocity,cPressure
!    call Conserve2Primitive(Geom,CC)

    allocate(TmpQuantity(Geom%CellNumber(1)))
    write(cStep,*) iStep
    iProgress = int(dble(iStep)/10.0d0)

    write(6,*) "Density.txt Generating ..."
    open(unit=2,file="Density.txt",status='unknown')
        do iLoop = 0, iStep
            write(cDensity, '("Density/Density",i6.6,".txt")') iLoop

            open(unit=1,file=cDensity,status='unknown')
                do iCenterX=1, Geom%CellNumber(1)
                    read(1,*) Tmp,TmpQuantity(iCenterX)
                end do
            close(1)

            do iCenterX=1, Geom%CellNumber(1)
                write(2,*) Geom%Bound(1,1) + 2.0d0 * Geom%Width(1,1,1) * (dble(iCenterX)-0.5d0),TmpQuantity(iCenterX)
            end do
            write(2,*) ""
            write(2,*) ""

            if(mod(iLoop,iProgress) == 0) then
                write(6,"((1x,i3),'% Complete')") int(dble(iLoop)/dble(iProgress)*10.0d0)
            end if
        end do
    close(2)

    write(6,*) "Velocity.txt Generating..."
    open(unit=2,file='Velocity.txt',status='unknown')
        do iLoop = 0, iStep
            write(cVelocity, '("Velocity/Velocity",i6.6,".txt")') iLoop

            open(unit=1,file=cVelocity,status='unknown')
                do iCenterX=1, Geom%CellNumber(1)
                    read(1,*) Tmp,TmpQuantity(iCenterX)
                end do
            close(1)

            do iCenterX=1, Geom%CellNumber(1)
                write(2,*) Geom%Bound(1,1) + 2.0d0*Geom%Width(1,1,1)*(dble(iCenterX)-0.5d0), TmpQuantity(iCenterX)
            end do
            write(2,*) ""
            write(2,*) ""
            if(mod(iLoop,iProgress) == 0) then
                write(6,"((1x,i3),'% Complete')") int(dble(iLoop)/dble(iProgress)*10.0d0)
            end if
        end do
    close(2)

    write(6,*) "Pressure.txt Generating...."
    open(unit=2,file='Pressure.txt',status='unknown')
        do iLoop = 0, iStep
            write(cPressure, '("Pressure/Pressure",i6.6,".txt")') iLoop

            open(unit=1,file=cPressure,status='unknown')
                do iCenterX=1, Geom%CellNumber(1)
                    read(1,*) Tmp,TmpQuantity(iCenterX)
                end do
            close(1)

            do iCenterX=1, Geom%CellNumber(1)
                write(2,*) Geom%Bound(1,1) + 2.0d0*Geom%Width(1,1,1)*(dble(iCenterX)-0.5d0), TmpQuantity(iCenterX)
            end do

            write(2,*) ""
            write(2,*) ""
            if(mod(iLoop,iProgress) == 0) then
                write(6,"((1x,i3),'% Complete')") int(dble(iLoop)/dble(iProgress)*10.0d0)
            end if
        end do
    close(2)

return
end subroutine MergeGnuPlot


text:
subroutine shafga(xieta,ndime,ltopo,ngaus,shaga)

!-----------------------------------------------------------------------
!
! This routine evaluates shape functions associated to gauss points
! for linear and quadratic isoparametric elements
!
!-----------------------------------------------------------------------
  use      typre
  implicit none
  integer(ip), intent(in)  :: ndime,ltopo,ngaus
  real(rp),    intent(in)  :: xieta(ndime)
  real(rp),    intent(out) :: shaga(ngaus)

  select case(ndime)
  case(1)
     call shaga1(xieta(1),ngaus,shaga)
  case(2)
     call shaga2(xieta(1),xieta(ndime),ltopo,ngaus,shaga)
  case(3)
     call shaga3(xieta(1),xieta(2),xieta(ndime),ltopo,ngaus,shaga)
  end select

end subroutine shafga

max_stars_repo_name: ciaid-colombia/InsFEM

max_issues_repo_name: ciaid-colombia/InsFEM

max_forks_repo_name: ciaid-colombia/InsFEM


text:
| Example of GCD algorithm:

: gcd ( a b -- gcd )
	0; tuck mod gcd ;

: .gcd ( a b -- )
	." GCD of "
	2dup . ." and " . ." is:"
	gcd . cr ;

100 80 .gcd
23 17 .gcd
64 12 .gcd
5 0 .gcd
bye

max_stars_repo_name: ronaaron/reva

max_issues_repo_name: ronaaron/reva

max_forks_repo_name: ronaaron/reva


text:
      SUBROUTINE CCHOL(N,ND,A,FAIL)
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION A(ND,ND)
      INTEGER FAIL
C
C
C SUBROUTINE CCHOL COMPUTES CHOLESKI
C DECOMPOSITION OF GIVEN COMPLEX POSITIVE DEFINITE
C MATRIX A.
C
C INPUT DATA
C
C  N --- ORDER OF MATRIX
C  ND -- DIMENSION OF ARRAY A (IT CAN BE
C        GREATER THAN OR EQUAL TO N)
C  A --- GIVEN MATRIX
C        IT IS SUPPOSED TO BE STORED IN THE
C        FOLLOWING WAY.  DIAGONAL ELEMENTS,
C        BEING REAL, ARE STORED ON THE DIAGONAL,
C        REAL PARTS OF OFFDIAGONAL ELEMENTS
C        ARE STORED IN THE LOWER TRIANGLE OF A,
C        IMAG PARTS OF THE LOWER TRIANGLE ARE
C        STORED IN THE UPPER TRIANGLE OF A.
C
C     EXIT INFORMATION
C
C  A --- COMPLEX ELEMENTS OF MATRIX L, DEFINED BY
C     A=L*L(H)
C        ARE STORED IN THE SAME WAY AS ORIGINAL
C        ELEMENTS OF A, THAT IS REAL PARTS OF THE
C        LOWER TRIANGLE OF L IN THE LOWER TRIANGLE
C        OF A AND THE CORRESPONDING IMAG PARTS IN
C        THE UPPER TRIANGLE OF A.
C  FAIL --- IS SET TO ZERO IF THE DECOMPOSITION WAS
C           SUCCESFUL AND TO NONZERO IF
C           THE MATRIX WAS NOT POSITIVE DEFINITE.
C
C
C    PROGRAMMED BY E. ZAKRAJSEK
C     JUNE 20, 1974
C
C
C
C     SUPPOSE DECOMPOSITION WILL FAIL
C
      FAIL=1
C
      DO 13 I=1,N
C
C     TEST FOR POSITIVE DEFINITNESS
C
      IF(A(I,I).LE.0.0D0) RETURN
C
C      COMPUTE COLUMN I
C
      A(I,I)=DSQRT(A(I,I))
      IF(I.EQ.N) GO TO 13
      IA=I+1
      DO 10 J=IA,N
      A(J,I)=A(J,I)/A(I,I)
   10 A(I,J)=A(I,J)/A(I,I)
C
C     REDUCE REMAINING COLUMNS
C
      DO 12 K=IA,N
      A(K,K)=A(K,K)-A(K,I)*A(K,I)-A(I,K)*A(I,K)
      IF(K.EQ.N) GO TO 12
      KA=K+1
      DO 11 J=KA,N
      A(J,K)=A(J,K)-A(J,I)*A(K,I)-A(I,J)*A(I,K)
   11 A(K,J)=A(K,J)-A(I,J)*A(K,I)+A(J,I)*A(I,K)
   12 CONTINUE
   13 CONTINUE
      FAIL=0
      RETURN
      END
      SUBROUTINE CTRED2(N,ND,A,B,D,E,F)
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION A(ND,ND),B(ND,ND),D(ND),E(ND),F(ND)
C
C
C     SUBROUTINE CTRED2 REDUCES GIVEN COMPLEX
C     HERMITIAN MATRIX TO A TRIDIAGONAL FORM
C
C     PARAMETERS
C
C     N    --- ORDER OF THE MATRIX
C     ND   --- DIMENSION OF ARRAYS A AND B
C     A    --- GIVEN MATRIX, REPLACED BY REAL PART
C              OF THE TRANSFORMATION MATRIX
C     B    --- IMAG PART OF TRANSFORMATION MATRIX
C     D    --- DIAGONAL PART OF THE TRIADIAGONAL MATRIX
C     E    --- REAL PART OF THE CODIAGONAL OF THE
C              TRIDIAGONAL MATRIX
C              (LAST N-1 LOCATIONS)
C     F    --- IMAG PARTS OF THE LOWER CODIAGONAL.
C
C     THE GIVEN MATRIX SHOULD BE STORED IN THE
C     FOLLOWING WAY
C
C          --- DIAGONAL ELEMENTS IN THE DIAGONAL
C          --- REAL PART OF THE LOWER TRIANGLE IN THE
C              LOWER TRIANGLE
C          --- IMAG PARTS OF THE LOWER TRIANGLE
C              IN THE UPPER TRIANGLE
C
C
C     PROGRAMMED BY E. ZAKRAJSEK
C     JUNE 20,1974
C
C
C
      CHEP=2.0D0**(-56)
      D(1)=A(1,1)
      IF(N.EQ.1) GO TO 31
C
C MAIN K LOOP
C
      DO 30 K=2,N
      L=K-1
C
C     COMPUTE NORM
C
      ALL=0.D0
      DO 10 I=K,N
   10 ALL=ALL+A(I,L)*A(I,L) + A(L,I)*A(L,I)
      ALL=DSQRT(ALL)
C
C     COMPUTE CONSTANTS
C
      C=1.0D0
      S=0.D0
      R = DSQRT(A(K,L)*A(K,L) + A(L,K)*A(L,K))
      IF(DABS(R).LT.1.D-50) R=0.D0
      IF(R.EQ.0.0D0) GOTO 11
      C=A(K,L)/R
      S=A(L,K)/R
   11 ALR=ALL*C
      ALI=ALL*S
      A(L,L)=0.0D0
C
C     TEST FOR SUPERFLUOUS TRANSFORMATION
C
      SM=ALL*(ALL+R)
      IF(DABS(SM).LT.1.D-50) SM=0.D0
      IF(SM.EQ.0.D0) GO TO 20
      G=1.0D0/SM
      A(L,L)=G
C
      A(K,L)=A(K,L)+ALR
      A(L,K) =A(L,K)+ALI
C
C     NOW COMPUTE U=A*W
C     AND STORE INTO (E,F)
C
      T=0.0D0
      DO 16 I=K,N
      C=A(I,I)*A(I,L)
      S=A(I,I)*A(L,I)
      IF(I.EQ.K)GOTO 13
      IA=I-1
      DO 12 J=K,IA
      C=C+A(I,J)*A(J,L)-A(J,I)*A(L,J)
   12 S=S+A(I,J)*A(L,J)+A(J,I)*A(J,L)
   13 IF(I.EQ.N) GOTO 15
      IA=I+1
      DO 14 J=IA,N
      C=C+A(J,I)*A(J,L)+A(I,J)*A(L,J)
   14 S=S+A(J,I)*A(L,J)-A(I,J)*A(J,L)
   15 E(I)=G*C
      F(I)=G*S
   16 T=T+A(I,L)*C+A(L,I) *S
      T=T*(G*G)
C
C    TRANSFORM  MATRIX
C
      DO 18 I=K,N
      A(I,I)=A(I,I)-2.0D0*(A(I,L)*E(I)+A(L,I)*F(I))+
     *      T*(A(I,L)*A(I,L)+A(L,I)*A(L,I))
      IF(I.EQ.K) GOTO 18
      IA=I-1
      DO 17 J=K,IA
      A(I,J)=A(I,J)-A(I,L)*E(J)-A(L,I)*F(J)
     *      -A(J,L)*E(I)-A(L,J)*F(I)
     *      +T*(A(I,L)*A(J,L)+A(L,I)*A(L,J))
      A(J,I)=A(J,I)-A(L,I)*E(J)+A(I,L)*F(J)
     *      +A(L,J)*E(I)-A(J,L)*F(I)
     *      +T*(A(L,I)*A(J,L)-A(I,L)*A(L,J))
   17 CONTINUE
   18 CONTINUE
C
C     STORE DIAGONAL AND CODIAGONAL ELEMENTS
C
   20 D(K)=A(K,K)
      E(K)=-ALR
      F(K)=-ALI
   30 CONTINUE
C
C     NOW ACCUMULATE TRANSFORMATIONS
C
   31 A(N,N)=1.D0
      B(N,N)=0.D0
      IF(N.EQ.1) RETURN
      DO 40 KK=2,N
      K=N-KK+2
      L=K-1
C
C     SKIP TRANSFORMATION IF UNIT
C
      IF(DABS(A(L,L)).LT.1.D-50) A(L,L)=0.D0
      IF(A(L,L).EQ.0.0D0) GO TO 36
C
C     COMPUTE PRODUCT
C
      DO 35 J=K,N
      C=0.D0
      S=0.D0
      DO 33 I=K,N
      C=C+A(I,L)*A(I,J)+A(L,I)*B(I,J)
   33 S=S+A(I,L)*B(I,J)-A(L,I)*A(I,J)
      C=C*A(L,L)
      S=S*A(L,L)
      DO 34 I=K,N
      A(I,J)=A(I,J)-C*A(I,L)+S*A(L,I)
   34 B(I,J)=B(I,J)-C*A(L,I)-S*A(I,L)
   35 CONTINUE
C
C     MAKE NEW LINE
C
   36 DO 37 I=K,N
      A(I,L)=0.D0
      A(L,I)=0.D0
      B(I,L)=0.D0
   37 B(L,I)=0.D0
      A(L,L)=1.D0
      B(L,L)=0.D0
   40 CONTINUE
      RETURN
      END
      SUBROUTINE CTQL2(N,ND,D,E,F,A,B,FAIL)
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION A(ND,ND), B(ND,ND), D(ND), E(ND), F(ND)
      INTEGER FAIL
C
C
C     SUBROUTINE CTQL2 COMPUTES THE EIGENVALUES AND
C     EIGENVECTORS OF A COMPLEX HERMITIAN TRIDIAGONAL
C     MATRIX
C
C     PARAMETERS
C
C     N    --- ORDER OF MATRIX
C     ND   --- DIMENSION OF A AND B
C     D    --- DIAGONAL GIVEN
C     E    --- REAL PART OF CODIAGONAL GIVEN
C              (LAST N-1 LOCATIONS)
C     F    --- IMAG PART OF THE LOWER CODIAGONAL
C     A    --- REAL PART OF EIGENVECTORS
C     B    --- IMAG PART OF EIGENVECTORS
C     FAIL --- RECEIVES VALUE OF 1 INSTEAD OF ZERO
C              IF SOME EIGENVALUE TAKES MORE THAN 30
C              ITERATIONS.
C
C
C     EIGENVALUES ARE OBTAINED IN INCREASING OF
C     MAGNITUDE IN VECTOR D, EIGENVECTORS ARE STORED
C     BY COLUMNS.  ARRAYS A AND B SHOULD BE PRESET TO
C     SOME UNITARY MATRIX SUCH AS THE IDENTITY MATRIX
C     OR THE MATRIX PRODUCED BY CTRED2.
C
C
C     PROGRAMMED BY E.  ZAKRAJSEK
C     JUNE 21, 1974
C
C
C
C
C     ***************************************
C     *                                     *
C     * NEXT LINE OF PROGRAM DEFINES        *
C     * MACHINE DEPENDENT CONSTANT CHEP     *
C     * DEFINED AS THE SMALLEST REAL        *
C     * NUMBER FOR WHICH                    *
C     *                                     *
C     *        1.0+CHEP .GT. 1.0            *
C     *                                     *
C     ***************************************
C
      CHEP=2.0D0**(-56)
C
C     FIRST MAKE REAL CODIAGONAL MOVED DOWN
C     TO FIRST LOCATION
C
      IF(N.EQ.1) GO TO 12
      DO 11 K=2,N
      R=DSQRT(E(K)*E(K)+F(K)*F(K))
      IF(DABS(R).LT.1.D-50) R=0.D0
      IF(R.EQ.0.0D0) GO TO 11
      C=E(K)/R
      S=F(K)/R
C
C     ACCUMULATE ROTATION
C
      DO 10 I=1,N
      P=A(I,K)*C-B(I,K)*S
      B(I,K)=A(I,K)*S+B(I,K)*C
   10 A(I,K)=P
C
C     TRANSFORM NEXT E
C
      IF(K.EQ.N) GO TO 11
      L=K+1
      P=E(L)*C-F(L)*S
      F(L)=E(L)*S+F(L)*C
      E(L)=P
   11 E(K-1)=R
   12 E(N)=0.D0
C
C     INITIALIZE
C
      BB=0.D0
      FF=0.D0
      FAIL=1
C
C     MAIN LOOP
C
      DO 32 L=1,N
      J=0
      H=CHEP*(DABS(D(L))+DABS(E(L)))
      IF(BB.LT.H) BB=H
C
C     LOOK FOR SMALL SUBDIAGONAL ELEMENT
C
      DO 20 M=L,N
      IF(DABS(E(M)).LE.BB) GO TO 21
   20 CONTINUE
   21 IF(M.EQ.L) GO TO 31
C
C     NEXT ITERATION
C
   24 IF(J.EQ.30) RETURN
      J=J+1
C
C     FORM SHIFT
C
      P=(D(L+1)-D(L))/(2.0D0*E(L))
      R=DSQRT(1.0D0+P*P)
      H=D(L)-E(L)/(P+DSIGN(R,P))
C
      DO 25 I=L,N
   25 D(I)=D(I)-H
      FF=FF+H
C
C     QL TRANSFORMATION
C
      P=D(M)
      C=1.0D0
      S=0.0D0
      MA=M-1
C
      DO 30 IA=L,MA
      I=MA-IA+L
      I1=I+1
      G=C*E(I)
      H=C*P
      IF(DABS(P).LT.DABS(E(I))) GO TO 26
C
      C=E(I)/P
      R=DSQRT(C*C+1.0D0)
      E(I1 )=S*P*R
      S=C/R
      C=1.0D0/R
      GO TO 27
C
   26 C=P/E(I)
      R=DSQRT(C*C+1. D0)
      E(I1 )=S*E(I)*R
      S=1.0D0/R
      C=C/R
C
   27 P=C*D(I)-S*G
      D(I1 )=H+S*(C*G+S*D(I))
C
C     FORM VECTOR
C
      DO 28 K=1,N
      HR=A(K,I1 )
      HI=B(K,I1 )
      A(K,I1 )=S*A(K,I)+C*HR
      B(K,I1 )=S*B(K,I)+C*HI
      A(K,I)=C*A(K,I)-S*HR
   28 B(K,I)=C*B(K,I)-S*HI
C
   30 CONTINUE
C
      E(L)=S*P
      D(L)=C*P
      IF(DABS(E(L)).GT.BB) GO TO 24
C
C     ROOT FOUND
C
   31 D(L)=D(L)+FF
   32 CONTINUE
C
C     ORDER EIGENVALUES AND EIGENVECTORS
C
      DO 42 I=1,N
      K=I
      DO 40 J=I,N
      IF(D(J).LT.D(K)) K=J
   40 CONTINUE
C
C
      IF(K.EQ.I) GO TO 42
      P=D(I)
      D(I)=D(K)
      D(K)=P
C
      DO 41 J=  1,N
      P=A(J,I)
      A(J,I)=A(J,K)
      A(J,K)=P
C
      P=B(J,I)
      B(J,I)=B(J,K)
   41 B(J,K)=P
   42 CONTINUE
      FAIL=0
      RETURN
      END

max_stars_repo_name: richardjgowers/yaehmop

max_issues_repo_name: richardjgowers/yaehmop

max_forks_repo_name: richardjgowers/yaehmop


text:
	SUBROUTINE IM_NEX2  ( imgfil, iret )
C************************************************************************
C* IM_NEX2								*
C*									*
C*									*
C* IM_NEX2  ( IMGFIL, IRET )						*
C*									*
C* Input parameters:							*
C*	IMGFIL		CHAR*		Image file			*
C*									*
C* Output parameters:							*
C*	IRET		INTEGER		Return code			*
C*					  0 = normal return		*
C*					 -2 = Error opening/reading file*
C**									*
C* Log:									*
C* S. Chiswell	Unidata	 4/04	Modified to use GSATMG			*
C************************************************************************
	INCLUDE		'GEMPRM.PRM'
	INCLUDE		'IMGDEF.CMN'
	INCLUDE		'AREAFL.CMN'
	INCLUDE		'imcmn.cmn'
C*
	CHARACTER*(*)	imgfil
C*
	CHARACTER	filnam*160, proj*4
	CHARACTER	radar_field*8, sweep_num*8
	CHARACTER	lev2tmp*20, tmpl*80, radid*4
C------------------------------------------------------------------------
	iret = 0
C
C*	Add a NULL to the end of the file name.
C
	CALL ST_NULL ( imgfil, filnam, lenf, ier )
C
	IF ( ICHAR(im_sweep_num(1:1)) .eq. 0 ) THEN
		sweep_num = '0'
	ELSE
		sweep_num = im_sweep_num
	END IF
C
	IF ( ICHAR(im_radar_field(1:1)) .eq. 0 ) THEN
                radar_field = 'dz'
        ELSE
                radar_field = im_radar_field
        END IF
C
	CALL ST_NULL ( sweep_num, sweep_num, lenf, ief )
	CALL ST_NULL ( radar_field, radar_field, lenf, ief )
C
C*	see if Radar ID can be extracted from file name (as a failsafe)
C
	radid = '    '
	lev2tmp = 'NEXRII'
	CALL ST_NULL (lev2tmp, lev2tmp, lenf, ief)
	CALL ctb_dttmpl ( lev2tmp, tmpl, ief )
	IF ( ief .eq. 0 ) THEN
	    CALL ST_RNUL ( tmpl, tmpl, ilen1, ief)
	    CALL ST_RPST ( tmpl, '%SITE%', '****', ipos, tmpl, ief)
	    IF ( ( ief .eq. 0 ) .and. ( ipos .gt. 0 ) ) THEN
		CALL ST_LSTR ( tmpl, ilen1, ief )
		CALL ST_LSTR ( imgfil, ilen2, ief )
		i = ilen2
		DO WHILE ( ( i .gt. 0 ) .and. (imgfil (i:i) .ne. '/' ) )
		   i = i - 1
		END DO
		IF ( ilen1 .eq. (ilen2 - i) ) THEN
		    ipos = ipos + i
		    radid = imgfil(ipos:ipos+3)
		END IF
	    END IF
	END IF
C
C*	Read the Level II Archive file header and get necessary information.
C
	CALL IM_RNEX2 ( filnam, radid, sweep_num, radar_field, imnpix,
     +			imnlin, imsorc, imtype, imdate, imtime, range,
     +			belv, iswpnum, clat, clon,
     +			xll, xur, yll, yur, vel_res, ier )
	IF  ( ier .ne. 0 )  THEN
	    iret = -2
	    RETURN
	END IF
C
C*	Set the remaining image definition common variables.
C
	imldat = imnpix * imnlin
C
C*	kilometers per pixel of image
C
	rmxres = ( 2. * range ) / float(imnpix - 1)
	rmyres = (2. * range) / float(imnlin - 1)
	rmxysc = rmxres / rmyres
	imdpth = 1
	imradf = 1
	imleft = 1
	imtop  = 1
	imrght = imnpix
	imbot  = imnlin
	rmbelv = belv

	cmstyp = radar_field
	imnchl = iswpnum
C	CALL ST_NUMB ( sweep_num, imnchl, ier )
C
C*	Set the map projection.
C
	iadir (6)  = 1
	iadir (7)  = 1
	iadir (12) = 1
	iadir (13) = 1
	CALL ST_CTOI ( 'RADR', 1, ianav (1), ier )
        ianav (2)  = imnlin / 2
        ianav (3)  = imnpix / 2
C
C*	Set the station lat/lon in navigation block
C
c	ianav (4) = clat_dms
C
C*      Set west lon (neg) to match MCIDAS west lon (pos)
C
c	ianav (5) = -1. * clon_dms
C
	deg = ABS ( clat )
        ideg = INT ( deg )
	fms = deg - float ( ideg )
	imin = INT ( fms * 60 )
	fsec = fms * 60. - float ( imin )
	isec = INT ( fsec * 60. + . 5 )
	IF ( isec .gt. 59 ) THEN
	   isec = 0
	   imin = imin + 1
	END IF
	IF ( imin .gt. 59 ) THEN
	   imin = 0
	   ideg = ideg + 1
	END IF
	ianav (4) = ( ( ideg * 100 ) + imin ) * 100 + isec
	IF ( clat .lt. 0 ) ianav (4) = -1 * ianav (4)

	deg = ABS ( clon )
        ideg = INT ( deg )
	fms = deg - float ( ideg )
	imin = INT ( fms * 60. )
	fsec = fms * 60 - float ( imin )
	isec = INT ( fsec * 60. + . 5 )
	IF ( isec .gt. 59 ) THEN
	   isec = 0
	   imin = imin + 1
	END IF
	IF ( imin .gt. 59 ) THEN
	   imin = 0
	   ideg = ideg + 1
	END IF
	ianav (5) = ( ( ideg * 100 ) + imin ) * 100 + isec
C
C*          Set west lon (neg) to match MCIDAS west lon (pos)
C
	IF ( clon .gt. 0 ) ianav (5) = -1 * ianav (5)
C
C	resolution
	ianav (6) = 1000 * rmxres
	ianav (7) = 0
	ianav (8) = 0

        CALL GSATMG ( imgfil, iadir, ianav, imleft, imtop, imrght,
     +                    imbot, ier )

	CALL IM_CALBAR_INIT ( ier )
C
        cmcalb = 'PROD'
	imcalbar = 1
C
C*	Read in min and max pixel and corresponding data
C*	values. Scale value is for data points.
C
	if ( imtype .eq. 225 ) then
           cmbunt = 'dBZ'
	   CALL IM_CALBAR ( 0., 107., -32., 75., ier )
	   cmblev (1) = 'ND'
	   cmblev (108) = '>=75'
	   imndlv = 108
	else if ( imtype .eq. 226 ) then
           cmbunt = 'm s**-1'
	   if (vel_res .eq. .5 ) then
	      CALL IM_CALBAR ( 0., 255., -63.5, 63.5, ier )
	   else
	      iminval = -127.
	      imaxval = 127
	      CALL IM_CALBAR ( 0., 255., -127., 127., ier )
	   end if
	   cmblev (1) = 'ND'
	   cmblev (256) = 'RF'
	   imndlv = 256
	else
           cmbunt = 'm s**-1'
	   CALL IM_CALBAR ( 0., 255., -64., 64., ier )
	   cmblev (1) = 'ND'
	   cmblev (256) = 'RF'
	   imndlv = 256
	endif

	RETURN
	END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak