text:
PROGRAM amsr2lr
      IMPLICIT none
!From Vera Gerald program of 18 Dec 2013
!F90 and other modifications R. Grumbine

      INTEGER iunt
      PARAMETER ( iunt = 12 )

      CHARACTER(8) SUBSET
      CHARACTER(80)  HS1, HS2, HS3, Htmp

!      CHARACTER*(*)   FLT
      CHARACTER(20)   FLT
      PARAMETER     ( FLT = 'bufr.tab' )

!      REAL HDR(9), XLOC(2,1), PRDLR(5,12), PRDHR(5,2)
      REAL(8) HDR(9), XLOC(2,1), PRDLR(5,12), PRDHR(5,2)
!      REAL(8) HDR(9), XLOC(2,1), PRDLR(5,12), PRDHR(5,2)
      DATA HS1, HS2, HS3, Htmp &
    / 'SAID SIID YEAR MNTH DAYS HOUR MINU ORBN SLNM               ', &
      'CLATH CLONH                                                ', &
      'SCCF ALFR VIIRSQ ANPO TMBR                                 ', &
      'MTYP                                                       '/
!
      INTEGER icount, iu, ntank, is
      REAL sccf, alfr, virsq, anpo, btmp
      DOUBLE PRECISION slon, slat

      INTEGER itk, lubfi, idate, iret, isad, issid, myr, mon, mda, mhr
      INTEGER minu, iorbn, iscanl, idatc
      INTEGER nchan, err

      INTEGER openout, outunit, creturns
      LOGICAL ok, hfok, lfok, weatherok
!.................................

       OPEN  ( UNIT = iunt, FILE = FLT )

       CALL W3TAGB('DCODAMSR2',0095,0333,0077,'NP11   ')

      outunit = 52
      OPEN (outunit, FORM="UNFORMATTED", STATUS="NEW")
      creturns = openout(outunit)

      icount = 0
      iu = 11
      ntank = 1

      DO ITK = 1, NTANK
!-----------------------------------------------------------------------
!                       DECODE BUFR FILE
!-----------------------------------------------------------------------
!  OPEN THE BUFR INPUT FILE
!  ------------------------
      LUBFI=IU

      CALL OPENBF(LUBFI,'IN',iunt)
      CALL READMG(LUBFI,SUBSET,IDATE,iret)
      IF(iret.NE.0) THEN
        CALL   closbf(lubfi)
        IF (icount.lt.1) then
          WRITE(6, *) 'READMG ERROR W/ BUFR MSG:iret=', iret
          err = 23
          STOP
!          CALL exit(23)
          go to 99
        ENDIF
      ENDIF
!
      PRINT*, 'READING DATA FOR ', IDATE
!
!  READ A SUBSET - READ A MESSAGE WHEN NO SUBSETS - END WHEN NO MESSAGES
!  ---------------------------------------------------------------------
  150 CONTINUE
      CALL READSB(LUBFI, iret)
      IF(iret.NE.0) THEN
         CALL READMG(LUBFI, SUBSET, IDATE, iret)
         IF(iret.NE.0) GOTO 100
         GOTO 150
      ENDIF
!
      if(icount .LE. 5e8) go to 250
      PRINT *, ' ***** NUMBER OF OBS EXCEEDED 5e8 *****'
      STOP

  250 CONTINUE
!  -----------------------------------------------------------------
!  CALL UFBINT, UFBREP TO GET THE DATA
!
!   SAID SIID YEAR MNTH DAYS HOUR MINU ORBN SLNM
!   CLATH CLONH                                                ',
!   SCCF ALFR VIIRSQ ANPO TMBR
!  -----------------------------------------------
      CALL UFBINT(LUBFI,HDR, 9,  1,iret, HS1)
      CALL UFBrep(LUBFI,XLOC,2,  1,iret, HS2)
      CALL UFBrep(LUBFI,PRDLR,5,12,iret, HS3)
!      PRINT *,'iret hs3 = ',iret
      nchan = iret
!
      isad   = NINT(hdr(1))
      issid  = NINT(hdr(2))
      myr    = NINT(hdr(3))
      mon    = NINT(hdr(4))
      mda    = NINT(hdr(5))
      mhr    = NINT(hdr(6))
      minu   = NINT(hdr(7))
      iorbn  = NINT(hdr(8))
      iscanl = NINT(hdr(9))

!
        SLAT = XLOC(1,1)
        SLON = XLOC(2,1)

        if(slon.lt.-180.0.or.slon.gt.180.0) go to 151

          icount = icount + 1
! Add some checking for quality and being in latitudes suitable for sea
! ice.  31 July 2017 Robert Grumbine
          ok = .TRUE.
          IF (nchan .EQ. 2) THEN
            ok = hfok(PRDLR)
          ELSE
            ok = lfok(PRDLR)
          ENDIF
          ok = ok .AND. (slat > 25.0 .OR. slat < -40.0) !

          IF (ok) THEN
            CALL headerout(isad,myr,mon,mda,mhr,minu,slat,slon, nchan)
            DO is = 1,nchan
               sccf  = PRDLR(1,is)
               alfr  = PRDLR(2,is)
               virsq = PRDLR(3,is)
               anpo  = PRDLR(4,is)
               btmp  = PRDLR(5,is)
               CALL spotout(sccf, alfr, virsq, anpo, btmp)
            ENDDO
          ENDIF

         IF (MOD(icount,INT(1e5)) .EQ. 0) THEN
           WRITE(51,*)isad,myr,mon,mda,mhr,minu,slat,slon, nchan
         ENDIF
  50     FORMAT(1x,i3,1x,i4,3i2,i3,1x,2f11.5,I3)
!         WRITE(52+icount/5000000) sccf/1.e9,alfr,virsq,anpo,btmp
!         IF (MOD(icount,INT(1e5)) .EQ. 0) THEN
!           WRITE(51,54) sccf/1.e9,alfr,virsq,anpo,btmp
!         ENDIF
   54 FORMAT(1x,F6.2,1x,4F13.4)

  151             continue

!       get next report
      GO TO 150

  100 CONTINUE
!
         CALL   closbf(lubfi)
         IF (icount.lt.1) then
           WRITE(6,*) 'READMG ERROR W/ BUFR MSG:iret=',iret
           err = 23
!           CALL exit(23)
           STOP
           go to 99
         endif
!
      WRITE(6,*) 'AFTER WOOLENS INTERFACE:  icount= ',icount
!
      ENDDO
!
!.................
      WRITE(6,692) icount
  692 FORMAT(1x,'END OF JOB','count= ',i8)

      CALL W3TAGE('DCODNET3')

   99 CONTINUE
END
!---------------------------------------------------------------
LOGICAL FUNCTION hfok(obs)
      IMPLICIT none
      REAL(8) obs(5,2) ! i = parameters, 5 = tmbr; j = channel
      LOGICAL tmp
      tmp = .TRUE.
      IF (obs(5,1) > 285 .OR. obs(5,2) > 285. ) tmp = .FALSE.
      IF (obs(5,1) > obs(5,2) ) tmp = .FALSE.
      hfok = tmp
      RETURN
      END
      LOGICAL FUNCTION lfok(obs)
      IMPLICIT none
      REAL(8) obs(5,12) ! i = parameters, 5 = tmbr; j = channel
      LOGICAL tmp
      INTEGER i
      tmp = .TRUE.

      DO i = 1, 6
        IF (obs(5,2*i-1) > 285 .OR. obs(5,2*i) > 285 ) tmp = .FALSE.
        IF (obs(5,2*i-1) > obs(5,2*i) ) tmp = .FALSE.
      ENDDO

      lfok = tmp
      RETURN
      END
      LOGICAL FUNCTION weatherok(obs, lat)
! Weather filter.  Should only apply after regressions to AMSRE have
! been performed.
      IMPLICIT none
      REAL AMSR_GR37LIM, AMSR_GR24LIM
      PARAMETER (AMSR_GR37LIM = 0.046)
      PARAMETER (AMSR_GR24LIM = 0.045)
      REAL(8) obs(5,12), amsre(5,12)
      REAL(8) lat
      LOGICAL tmp
      REAL gr3719, gr2419

      amsre = obs
      CALL regress(amsre, lat)
      tmp = .FALSE.
      gr3719 = (amsre(5,12) - amsre(5,8) ) / (amsre(5,12)+amsre(5,8))
      gr2419 = (amsre(5,10) - amsre(5,8) ) / (amsre(5,10)+amsre(5,8))

      IF ((gr3719 < AMSR_GR37LIM) .AND. (gr2419 < AMSR_GR24LIM) ) THEN
        tmp = .TRUE.
      ENDIF
      weatherok = tmp
      RETURN
      END
      SUBROUTINE regress(amsre, lat)
      IMPLICIT none
      REAL(8) amsre(5,12)
      REAL(8) lat
!Perform limited regression -- just those channels that are used in
!weather filter
      IF (lat > 0) THEN
        amsre(5, 8) = amsre(5, 8)*1.031 - 9.710
        amsre(5,10) = amsre(5,10)*0.999 - 1.706
        amsre(5,12) = amsre(5,12)*0.997 - 2.610
      ELSE
        amsre(5, 8) = amsre(5, 8)*1.032 -10.013
        amsre(5,10) = amsre(5,10)*0.993 - 0.987
        amsre(5,12) = amsre(5,12)*0.995 - 2.400
      ENDIF

      RETURN
END

max_stars_repo_name: rgrumbine/seaice-concentration

max_issues_repo_name: rgrumbine/seaice-concentration

max_forks_repo_name: rgrumbine/seaice-concentration


text:
  !> @file
  !! Other fortran file for f_malloc routines
  !! @author
  !!    Copyright (C) 2012-2015 BigDFT group
  !!    This file is distributed under the terms of the
  !!    GNU General Public License, see ~/COPYING file
  !!    or http://www.gnu.org/copyleft/gpl.txt .
  !!    For the list of contributors, see ~/AUTHORS

  !if the src_ptr is nullified the f_malloc should provide nullification
  if (.not. associated(src_ptr)) then
     m%srcdata_add=int(-1,f_address)
     return
  end if

max_stars_repo_name: adegomme/aiida-bigdft-plugin


module NewStrategyM

  use ProcessM

  implicit none

  private
  public :: NewStrategyDT

  type, abstract :: NewStrategyDT
   contains
     procedure(NewStrategy_Procedure), nopass, deferred :: useNewStrategy
  end type NewStrategyDT

  abstract interface
     subroutine NewStrategy_Procedure(this)
       import ProcessDT
       class(ProcessDT), intent(inout) :: this
     end subroutine NewStrategy_Procedure
  end interface

end module NewStrategyM


text:
MODULE test_select
IMPLICIT NONE
!
TYPE :: parent
  INTEGER(4) :: n
END TYPE parent
!
TYPE, EXTENDS(parent) :: child1
  REAL(8) :: a
END TYPE child1
!
TYPE, EXTENDS(parent) :: child2
  COMPLEX(8) :: a
END TYPE child2
CONTAINS
!
SUBROUTINE test_select_sub(self)
CLASS(parent), INTENT(inout) :: self
! Select statement with binding
SELECT TYPE(this=>self)
TYPE IS(child1)
  this%a
CLASS IS(child2)
  this%a
CLASS DEFAULT
  this%n
END SELECT
! Select statement without binding
SELECT TYPE(self)
TYPE IS(child1)
  self%a
END SELECT
END SUBROUTINE test_select_sub
END MODULE test_select

max_stars_repo_name: regit/fortran-language-server

max_issues_repo_name: regit/fortran-language-server

max_forks_repo_name: regit/fortran-language-server


text:
      logical function occ_bound_p(cbound,occ1,occ2)

      implicit none

      include 'opdim.h'

       character, intent(in) ::
     &     cbound*(*)
      integer(8), intent(in) ::
     &     occ1, occ2

      integer(8) ::
     &     scr1, scr2
      integer ::
     &     iel

      occ_bound_p = .false.
      scr1 = occ1
      scr2 = occ2
      select case(trim(cbound))
      case('<','LT','lt')
        do iel = 1, ngastp*2
          if(mod(scr1,pack_base).ge.mod(scr2,pack_base)) exit
          scr1 = scr1 / pack_base
          scr2 = scr2 / pack_base
          occ_bound_p = iel.eq.ngastp*2
        end do
      case('<=','LE','le')
        do iel = 1, ngastp*2
          if(mod(scr1,pack_base).gt.mod(scr2,pack_base)) exit
          scr1 = scr1 / pack_base
          scr2 = scr2 / pack_base
          occ_bound_p = iel.eq.ngastp*2
        end do
      case('>','GT','gt')
        do iel = 1, ngastp*2
          if(mod(scr1,pack_base).le.mod(scr2,pack_base)) exit
          scr1 = scr1 / pack_base
          scr2 = scr2 / pack_base
          occ_bound_p = iel.eq.ngastp*2
        end do
      case('>=','GE','ge')
        do iel = 1, ngastp*2
          if(mod(scr1,pack_base).lt.mod(scr2,pack_base)) exit
          scr1 = scr1 / pack_base
          scr2 = scr2 / pack_base
          occ_bound_p = iel.eq.ngastp*2
        end do
      case default
        call quit(1,'occ_bound_p','Unknown relation: "'//cbound//'"')
      end select

      return
      end

ext: f

max_stars_repo_name: ak-ustutt/GeCCo-public

max_issues_repo_name: ak-ustutt/GeCCo-public

max_forks_repo_name: ak-ustutt/GeCCo-public


text:
!==========================================================================
elemental function gsw_sstar_from_sp (sp, p, long, lat)
!==========================================================================
!
! Calculates Preformed Salinity, Sstar, from Practical Salinity, SP.
!
! sp     : Practical Salinity                              [unitless]
! p      : sea pressure                                    [dbar]
! long   : longitude                                       [deg E]
! lat    : latitude                                        [deg N]
!
! gsw_sstar_from_sp  : Preformed Salinity                  [g/kg]
!--------------------------------------------------------------------------

use gsw_mod_toolbox, only : gsw_saar, gsw_sa_from_sp_baltic

use gsw_mod_error_functions, only : gsw_error_code, gsw_error_limit

use gsw_mod_teos10_constants, only : gsw_ups

use gsw_mod_kinds

implicit none

real (r8), intent(in) :: sp, p, long, lat

real (r8) :: gsw_sstar_from_sp

real (r8) :: saar, sstar_baltic

character (*), parameter :: func_name = "gsw_sstar_from_sp"

!In the Baltic Sea, Sstar = SA.
sstar_baltic = gsw_sa_from_sp_baltic(sp,long,lat)

if (sstar_baltic .lt. 1e10_r8) then

   gsw_sstar_from_sp = sstar_baltic

else

   saar = gsw_saar(p,long,lat)
   if (saar .gt. gsw_error_limit) then
      gsw_sstar_from_sp = gsw_error_code(1,func_name,saar)
   else
      gsw_sstar_from_sp = gsw_ups*sp*(1 - 0.35_r8*saar)
   end if

end if

return
end function

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

max_stars_repo_name: minsukji/ci-debug

max_issues_repo_name: minsukji/ci-debug

max_forks_repo_name: minsukji/ci-debug


text:
      subroutine correctpc(iparent,isurv,iremov,iparsurv,iparremov,
     &   itettyp,iseedtet,itp1,itet,itetoff,jtet,jtetoff,nelts,
     &   ipielts,nef_cmo,mbndry,imt1,isn1,itetclr,lmoved)
C
C ######################################################################
C
C        $Log: correctpc.f,v $
C        Revision 2.00  2007/11/05 19:45:51  spchu
C        Import to CVS
C
CPVCS
CPVCS       Rev 1.21   02 Oct 2007 12:40:28   spchu
CPVCS    original version
C
C ######################################################################
C
      implicit none

      include 'local_element.h'
      include 'chydro.h'

      pointer (ipielts,ielts)
      integer ielts(*)

      integer itp1(*),iparent(*),itettyp(*),itet(*),itetoff(*),
     &   jtet(*),jtetoff(*),mbndry,imt1(*),isn1(*),
     &   itetclr(*)
      integer iparsurv,iparremov,matelt,nodk,ksurv,ichild,nod1,nod2,
     &   ityp,locnod,j,isurv,iremov,iseedtet,nod,ierrw,nelts,nef_cmo,
     &   ielt,ierr,k,itp
      logical lmoved(*),lset1,lset2,lreal
      character*132 cbuf,logmess

      ityp=itettyp(iseedtet)
      locnod=0
      do j=1,nelmnen(ityp)
         nod=itet(j+itetoff(iseedtet))
         if (iparsurv.eq.iparent(nod)) then
            locnod=j
         endif
      enddo
      if (locnod.eq.0) then
         write(logmess,'(a,i9,a,i9,a,i9)')
     &      'Try ',isurv,'<--',iremov,
     &      ':  FATAL ERROR : topological error at node ',isurv
         call writloga('default',0,logmess,0,ierrw)
         call termcode()
      endif

      call get_elements_around_node(iseedtet,locnod,nelts
     &   ,ipielts,itetoff,jtetoff,itet,jtet,itettyp,iparent,
     &   nef_cmo,mbndry)

c.... Check that there is child at ISURV for each material type
c.... appearing in the retrieved tet list.

      do 200 j=1,nelts
         ielt=ielts(j)
         ityp=itettyp(ielt)
         matelt=itetclr(ielt)

         do k=1,nelmnen(ityp)
            nodk=itet(k+itetoff(ielt))
            if (nodk.eq.isurv) then
               ksurv=k
               if (imt1(nodk).eq.matelt) then
                  goto 200
               else
                  goto 300
               endif
            endif
         enddo
         goto 200

 300     continue

         nod=isn1(nodk)
         do k=1,10000
            if (imt1(nod).eq.matelt.and.itp1(nod).ne.ifitpcup)
     &         then
               itet(ksurv+itetoff(ielt))=nod
               goto 200
            endif
            nod=isn1(nod)
            if (nod.eq.nodk) goto 210
         enddo
         write(logmess,'(a,i9,a,i9,a,i9)')
     &      'Try ',isurv,'<--',iremov,
     &      ':  FAIL : topological error at node ',isurv
         call writloga('default',0,logmess,0,ierrw)
         call termcode()

 210     continue

c.... There is no child for material MATELT of element IELT.
c.... This material must appear as a child at IREMOV.  Make this
c.... dudded child into the necessary child at ISURV.

         nod=isn1(iparent(iremov))
         if (nod.eq.0) then
            write(logmess,'(a,i9,a,i9,a,i9)')
     &         'Try ',isurv,'<--',iremov,
     &         ':  FATAL ERROR : topological error at node ',iremov
            call writloga('default',0,logmess,0,ierrw)
            call termcode()
         endif
         do k=1,10000
            if (imt1(nod).eq.matelt.and.itp1(nod).ne.ifitpcup)
     *         goto 220
            nod=isn1(nod)
            if (nod.eq.iparent(iremov)) goto 230
         enddo
         write(logmess,'(a,i9,a,i9,a,i9)')
     &      'Try ',isurv,'<--',iremov,
     &      ':  FATAL ERROR : topological error at node ',iremov
         call writloga('default',0,logmess,0,ierrw)
         call termcode()

 230     continue
         write(logmess,'(a,i9,a,i9,a,i9,a)')
     &      'Try ',isurv,'<--',iremov,
     &      ':  FATAL ERROR : material ',matelt,' unexpected.'
         call writloga('default',0,logmess,0,ierrw)
         call termcode()

 220     continue

c.... Copy data of a child at ISURV to the child at IREMOV.
c.... Then put this former child at IREMOV into the parent/child
c.... system at ISURV.  Die if first available child at ISURV is
c.... not real.

         ichild=isn1(iparsurv)
         itp=itp1(ichild)
         lset1=.false.
         if (itp.ge.ifitpst1.and.itp.le.ifitpen1) lset1=.true.
         lset2=.false.
         if (itp.ge.ifitpst2.and.itp.le.ifitpen2) lset2=.true.
         lreal=lset1.or.lset2
         if (.not.lreal) then
            write(logmess,'(a,i9,a,i9,a,i9,a)')
     &         'Try ',isurv,'<--',iremov,
     &         ':  FATAL ERROR.'
            call writloga('default',0,logmess,0,ierrw)
            call termcode()
         endif

c.... Take NOD out of parent/child system at IREMOV.
         iparent(nod)=iparsurv
         lmoved(nod)=.true.

         nod1=iparremov
         nod2=isn1(nod1)
         do k=1,10000
            if (nod2.eq.nod) then
               isn1(nod1)=isn1(nod)
               goto 370
            endif
            nod1=nod2
            nod2=isn1(nod2)
         enddo
         logmess='FATAL ERROR'
         call writloga('default',0,logmess,0,ierrw)
         call termcode()

 370     continue

c.... Put NOD into parent/child system at ISURV.

         write(cbuf,'(10(a,i9))') 'copypts/',
     &      ichild,',',ichild,', 1 / ',nod,', 1 ; finish'
         call dotaskx3d(cbuf,ierr)

c.... Correct color of NOD

         imt1(nod)=matelt

c
         nod1=iparsurv
         nod2=isn1(nod1)
         do k=1,10000
            if ((matelt.lt.imt1(nod2).and.nod1.eq.iparsurv).or.
     &         (matelt.lt.imt1(nod2).and.matelt.gt.imt1(nod1))
     &         .or.nod2.eq.iparsurv) then
               isn1(nod)=isn1(nod1)
               isn1(nod1)=nod
               goto 310
            endif
            nod1=nod2
            nod2=isn1(nod2)
         enddo
         logmess='FATAL ERROR'
         call writloga('default',0,logmess,0,ierrw)
         call termcode()

 310     continue
         itet(ksurv+itetoff(ielt))=nod
 200  continue

c.... Loop over child points and make sure each child material
c.... is the material of at least one incident tet.  If not,
c.... dud out this child.

      nod1=iparsurv
      nod2=isn1(iparsurv)

c.... We assume that iparsurv IS a parent node (type IFITPCUP)
c.... AND ISN1(IPARSURV).ne.IPARSURV.  If this isn't true, die.

      if (nod1.eq.nod2.or.itp1(nod1).ne.ifitpcup) then
         logmess='FATAL ERROR'
         call writloga('default',0,logmess,0,ierrw)
         call termcode()
      endif

      do k=1,10000
         do j=1,nelts
            if (itetclr(ielts(j)).eq.imt1(nod2)) goto 340
         enddo

c.... Dud out NOD2.

         isn1(nod1)=isn1(nod2)
         itp1(nod2)=ifitpmrg
         goto 350

c.... Advance to next node.

 340     continue
         nod1=nod2

 350     continue
         nod2=isn1(nod2)
         if (nod2.eq.iparsurv) goto 360

      enddo
      logmess='FATAL ERROR'
      call writloga('default',0,logmess,0,ierrw)
      call termcode()

 360  continue

c.... Check we aren't down to just one child.  If so,
c.... change point types and get rid of parent point.

      nod=isn1(iparsurv)
      if (nod.eq.iparsurv) then
         logmess='FATAL ERROR'
         call writloga('default',0,logmess,0,ierrw)
         call termcode()
      endif
      if (isn1(nod).eq.iparsurv) then
         isn1(nod)=0
         iparent(nod)=nod
         itp1(iparsurv)=ifitpmrg
         if (itp1(nod).eq.ifitpini) then
            itp1(nod)=ifitpint
         elseif (itp1(nod).eq.ifitpirb) then
            itp1(nod)=ifitprfl
         elseif (itp1(nod).eq.ifitpvif) then
            itp1(nod)=ifitpvfb
         elseif (itp1(nod).eq.ifitpalb) then
            itp1(nod)=ifitpvrf
         elseif (itp1(nod).eq.ifitpifb) then
            itp1(nod)=ifitpfre
         elseif (itp1(nod).eq.ifitpirf) then
            itp1(nod)=ifitprfb
         elseif (itp1(nod).eq.ifitpvir) then
            itp1(nod)=ifitpvrb
         else
            write(logmess,'(10(a,i9))')
     &         'Try ',isurv,'<--',iremov,
     &         ':  WARNING : node ',nod,' now has '
            call writloga('default',0,logmess,0,ierrw)
            write(logmess,'(10(a,i9))')
     &         'itp=',itp1(nod),' and only 1 material.'
     &
            call writloga('default',0,logmess,0,ierrw)
         endif
      endif

      return
      end

max_forks_repo_name: daniellivingston/LaGriT


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

20100915 01:58:48 nbsp Howdy Omid, and welcome to the wiki!    My names Evan, pleased to meet you!  I saw you uploaded an image but didnt add it to your user profile, so I figured Id try and help you out.  You can edit this (and any other!) entry by clicking the edit icon.  Once again, welcome to the wiki. Users/JabberWokky Evan JabberWokky Edwards

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
	SUBROUTINE GR_RBAN  ( anlblk, deltan, deltax, deltay, gbnds,
     +			      ebnds, dbnds, iextnd, iret )
C************************************************************************
C* GR_RBAN								*
C*									*
C* This subroutine reads a Barnes analysis block.  All the bounds	*
C* are returned in the order:  lower left latitude; lower left		*
C* longitude; upper right latitude; upper right longitude.		*
C*									*
C* GR_RBAN  ( ANLBLK, DELTAN, DELTAX, DELTAY, GBNDS, EBNDS, DBNDS,	*
C*            IEXTND, IRET )						*
C*									*
C* Input parameters:							*
C*	ANLBLK (LLNANL)	REAL		Analysis block			*
C*									*
C* Output parameters:							*
C*	DELTAN		REAL		Station spacing			*
C*	DELTAX		REAL		Grid spacing in x dir		*
C*	DELTAY		REAL		Grid spacing in y dir		*
C*	GBNDS (4)	REAL		Grid area bounds		*
C*	EBNDS (4)	REAL		Extend area bounds		*
C*	DBNDS (4)	REAL		Data area bounds		*
C*	IEXTND (4)	INTEGER		Extend grid points		*
C*	IRET		INTEGER		Return code			*
C*					  0 = normal return		*
C*					-10 = invalid analysis block	*
C**									*
C* Log:									*
C* M. desJardins/GSFC	 2/85						*
C* M. desJardins/GSFC	 9/88						*
C* K. Brill/GSC          4/90	Read for ANLBLK (1) = 2			*
C* M. desJardins/NMC	 4/91	Set DELTAX, DELTAY			*
C************************************************************************
	INCLUDE		'GEMPRM.PRM'
C*
	REAL		gbnds (*), ebnds (*), dbnds (*), anlblk (*)
	INTEGER		iextnd (*)
C------------------------------------------------------------------------
	iret = 0
C
C*	Check that this is a Barnes analysis block.
C
	IF  ( anlblk (1) .eq. 1.0 )  THEN
	    deltan = anlblk (2)
	    deltax = anlblk (3)
	    deltay = anlblk (4)
	    gamma  = anlblk (5)
	    DO  i = 1, 4
	      gbnds (i) = anlblk ( 5+i)
	      ebnds (i) = anlblk ( 9+i)
	      dbnds (i) = anlblk (13+i)
	    END DO
C
C*	    Compute IEXTND.
C
	    iextnd (1) = ( gbnds (2) - ebnds (2) ) / deltax + 0.1
	    iextnd (2) = ( gbnds (1) - ebnds (1) ) / deltay + 0.1
	    iextnd (3) = ( ebnds (4) - gbnds (4) ) / deltax + 0.1
	    iextnd (4) = ( ebnds (3) - gbnds (3) ) / deltay + 0.1
C*
	  ELSE IF  ( anlblk (1) .eq. 2.0 )  THEN
	    deltan = anlblk (2)
	    deltax = RMISSD
	    deltay = RMISSD
	    iextnd (1) = NINT ( anlblk (3) )
	    iextnd (2) = NINT ( anlblk (4) )
	    iextnd (3) = NINT ( anlblk (5) )
	    iextnd (4) = NINT ( anlblk (6) )
            DO i = 1, 4
	      gbnds (i) = anlblk (  6 + i )
	      ebnds (i) = anlblk ( 10 + i )
	      dbnds (i) = anlblk ( 14 + i )
	    END DO
	  ELSE
	    iret = -10
	END IF
C*
	RETURN
	END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


text:
The International Socialist Organization, or ISO, is a leftist organization which is currently (2007) inactive on the UCD campus due to lack of student interest as shown by scarce meeting attendance.

They hold lots of events that organize activists... in order to mobilize opposition to all forms of oppression and exploitation. Some say these topics suggest George W. Bush is evil.  There are a lot of flyers to these events posted on campus, each bearing a tagline Handbill only.  Do not post.

John Green was an officer of the ISO.

The ISO self identifies as a MarxistLeninistTrotskyist organization, which can be seen on their website which states that The International Socialist Organization (ISO) stands in the revolutionary tradition of Karl Marx, Vladimir Lenin and Leon Trotsky. the books they sell, reading groups they have, the content of the news paper, http://socialistworker.org/ , which includes a The Meaning of Marxism column. Though, it should be noted that Marxs theories sparked philosophies that headed in many different directions and Socialism is, and has historically been, different than what some would see as communism, it can be argued that the ISO in fact is a communist organization as seen by these people because it believes that Capitalism must be overthrown, and as Leninists, they push for socialism as the intermediary step to communism.

To clarify the standard meanings of socialism and communism, there are http://www.marxmail.org/faq/socialism_and_communism.htm plenty http://www.geocities.com/commiett/difference.html of http://www.romm.org/soc_com.html places http://www.pww.org/pastweeks1999/Socialism%20and%20communism.htm you http://www.tysknews.com/Articles/democratism_pt1.htm can http://www.newyouth.com/archives/theory/faq/socialism_communism.asp read wiki:WikiPedia:Socialism about wiki:WikiPedia:Communist it, though it appears that the ISO use the term socialist though some would label them as a communist organization.

 Where They Stand
According to http://bad.eserver.org/issues/2004/65/chretien.html Todd Chretien, Northern California Organizer for the International Socialist Organization, ...while we should fight for whatever reforms we can win — be it defending abortion rights, or strengthening the unions, or winning voting rights for undocumented immigrants, or working alongside the Greens to get independent candidates elected — whatever the reforms are, we believe that the point of a socialist organization is to say in the long term we need to build a revolutionary movement... to overthrow capitalism.

More about them can be read http://www.internationalsocialist.org/ here


!==============================================================================!
  subroutine Cgs(sol, x, b, prec, miter, niter, tol, fin_res, norm)
!------------------------------------------------------------------------------!
!   Solves the linear systems of equations by a precond. CGS Method.           !
!------------------------------------------------------------------------------!
!   Allows preconditioning of the system by:                                   !
!     1. Diagonal preconditioning                                              !
!     2. Incomplete Cholesky preconditioning                                   !
!                                                                              !
!   The type of precondtioning is chosen by setting the variable prec to 0     !
!   (for no preconditioning), 1 (for diagonal preconditioning) or 2 (for       !
!   incomplete Cholesky preconditioning)                                       !
!------------------------------------------------------------------------------!
!----------------------------------[Modules]-----------------------------------!
  use Work_Mod, only: p1         => r_cell_01,  &
                      p2         => r_cell_02,  &
                      q1         => r_cell_03,  &
                      q2         => r_cell_04,  &
                      r1         => r_cell_06,  &
                      r2         => r_cell_07,  &
                      u1         => r_cell_08,  &
                      u2         => r_cell_09,  &
                      v2         => r_cell_10,  &
                      u1_plus_q1 => r_cell_11
!------------------------------------------------------------------------------!
  implicit none
!---------------------------------[Arguments]----------------------------------!
  type(Solver_Type), target :: sol
  real              :: x(-sol % pnt_grid % n_bnd_cells :  &
                          sol % pnt_grid % n_cells)
  real              :: b( sol % pnt_grid % n_cells)  ! [A]{x}={b}
  character(len=80) :: prec                          ! preconditioner
  integer           :: miter                         ! max and actual ...
  integer           :: niter                         ! ... num. iterations
  real              :: tol                           ! tolerance
  real              :: fin_res                       ! residual
  real, optional    :: norm                          ! normalization
!-----------------------------------[Locals]-----------------------------------!
  type(Matrix_Type), pointer :: a
  type(Matrix_Type), pointer :: d
  integer                    :: nt, ni, nb
  real                       :: alfa, beta, rho, rho_old, bnrm2, res
  integer                    :: i, j, k, iter
!==============================================================================!

  ! Take some aliases
  a => sol % a
  d => sol % d
  nt = a % pnt_grid % n_cells
  ni = a % pnt_grid % n_cells - a % pnt_grid % comm % n_buff_cells
  nb = a % pnt_grid % n_bnd_cells

  res = 0.0

  !---------------------!
  !   Preconditioning   !
  !---------------------!
  call Prec_Form(ni, a, d, prec)

  !-----------------------------------!
  !    This is quite tricky point.    !
  !   What if bnrm2 is very small ?   !
  !-----------------------------------!
  if(.not. present(norm)) then
    bnrm2 = Normalized_Root_Mean_Square(ni, b(1:nt), a, x(1:nt))
  else
    bnrm2 = Normalized_Root_Mean_Square(ni, b(1:nt), a, x(1:nt), norm)
  end if

  if(bnrm2 < tol) then
    iter=0
    goto 1
  end if

  !-----------------!
  !   r1 = b - Ax   !
  !-----------------!
  call Residual_Vector(ni, r1(1:nt), b(1:nt), a, x(1:nt))

  !--------------------------------!
  !   Calculate initial residual   !
  !--------------------------------!
  res = Normalized_Root_Mean_Square(ni, r1(1:nt), a, x(1:nt))

  if(res < tol) then
    iter=0
    goto 1
  end if

  !-------------!
  !   r2 = r1   !
  !-------------!
  r2(1:ni) = r1(1:ni)

  !---------------!
  !               !
  !   Main loop   !
  !               !
  !---------------!
  do iter = 1, miter

    !-------------------!
    !   rho = (r2,z1)   !
    !-------------------!
    rho = dot_product(r1(1:ni), r2(1:ni))
    call Comm_Mod_Global_Sum_Real(rho)

    if(iter .eq. 1) then
      u1(1:ni) = r1(1:ni)
      u2(1:ni) = u1(1:ni)
    else
      beta = rho / rho_old
      u1(1:ni) = r1(1:ni) + beta *  q1(1:ni)
      u2(1:ni) = u1(1:ni) + beta * (q1(1:ni) + beta*u2(1:ni))
    end if

    !---------------------!
    !   Solve M p2 = u2   !
    !---------------------!
    call Prec_Solve(ni, a, d, p2(1:nt), u2(1:nt), prec)

    !--------------!
    !   v2 = Ap2   !
    !--------------!
    call Grid_Mod_Exchange_Cells_Real(a % pnt_grid, p2(-nb:ni))
    do i = 1, ni
      v2(i) = 0.0
      do j = a % row(i), a % row(i+1)-1
        k = a % col(j)
        v2(i) = v2(i) + a % val(j) * p2(k)
      end do
      alfa = alfa + r2(i) * v2(i)
    end do

    !------------------------!
    !   alfa = rho/(r2,v2)   !
    !------------------------!
    alfa = dot_product(r2(1:ni), v2(1:ni))
    call Comm_Mod_Global_Sum_Real(alfa)
    alfa=rho / alfa

    !-------------------------!
    !   q1 = u1 - alfa * v2   !
    !-------------------------!
    q1(1:ni) = u1(1:ni) - alfa * v2(1:ni)

    !-------------------------!
    !   solve Mp1 = u1 + q1   !
    !-------------------------!
    u1_plus_q1(1:ni) = u1(1:ni) + q1(1:ni)
    call Prec_Solve(ni, a, d, p1(1:nt), u1_plus_q1(1:nt), prec)

    !---------------------!
    !   x = x + alfa p1   !
    !---------------------!
    x(1:ni) = x(1:ni) + alfa * p1(1:ni)

    !---------------!
    !   q2 = A p1   !
    !---------------!
    call Grid_Mod_Exchange_Cells_Real(a % pnt_grid, p1(-nb:ni))
    do i = 1, ni
      q2(i) = 0.0
      do j = a % row(i), a % row(i+1)-1
        k = a % col(j)
        q2(i) = q2(i) + a % val(j) * p1(k)
      end do
    end do

    !---------------------!
    !   r = r - alfa q2   !
    !---------------------!
    r1(1:ni) = r1(1:ni) - alfa * q2(1:ni)

    !-----------------------!
    !   Check convergence   !
    !-----------------------!
    if(.not. present(norm)) then
      res = Normalized_Root_Mean_Square(ni, r1(1:nt), a, x(1:nt))
    else
      res = Normalized_Root_Mean_Square(ni, r1(1:nt), a, x(1:nt), norm)
    end if

    if(res < tol) goto 1

    rho_old=rho

  end do                ! iter

  !----------------------------------!
  !                                  !
  !   Convergence has been reached   !
  !                                  !
  !----------------------------------!
1 continue
  fin_res = res
  niter  = iter

  end subroutine


text:
      subroutine VSOP87 (tdj,ivers,ibody,prec,lu,r,ierr)
*-----------------------------------------------------------------------
*
*     Reference : Bureau des Longitudes - PBGF9502
*
*     Object :
*
*     Substitution of time in VSOP87 solution written on a file.
*     The file corresponds to a version of VSOP87 theory and to a body.
*
*     Input :
*
*     tdj      julian date (real double precision).
*              time scale : dynamical time TDB.
*
*     ivers    version index (integer).
*              0: VSOP87 (initial solution).
*                 elliptic coordinates
*                 dynamical equinox and ecliptic J2000.
*              1: VSOP87A.
*                 rectangular coordinates
*                 heliocentric positions and velocities
*                 dynamical equinox and ecliptic J2000.
*              2: VSOP87A.
*                 spherical coordinates
*                 heliocentric positions and velocities
*                 dynamical equinox and ecliptic J2000.
*              3: VSOP87C.
*                 rectangular coordinates
*                 heliocentric positions and velocities
*                 dynamical equinox and ecliptic of the date.
*              4: VSOP87D.
*                 spherical coordinates
*                 heliocentric positions and velocities
*                 dynamical equinox and ecliptic of the date.
*              5: VSOP87E.
*                 rectangular coordinates
*                 barycentric positions and velocities
*                 dynamical equinox and ecliptic J2000.
*
*     ibody    body index (integer).
*              0: Sun
*              1: Mercury
*              2: Venus
*              3: Earth
*              4: Mars
*              5: Jupiter
*              6: Saturn
*              7: Uranus
*              8: Neptune
*              9: Earth-Moon barycenter
*
*     prec     relative precision (real double precision).
*
*              if prec is equal to 0 then the precision is the precision
*                 p0 of the complete solution VSOP87.
*                 Mercury    p0 =  0.6 10**-8
*                 Venus      p0 =  2.5 10**-8
*                 Earth      p0 =  2.5 10**-8
*                 Mars       p0 = 10.0 10**-8
*                 Jupiter    p0 = 35.0 10**-8
*                 Saturn     p0 = 70.0 10**-8
*                 Uranus     p0 =  8.0 10**-8
*                 Neptune    p0 = 42.0 10**-8
*
*              if prec is not equal to 0, let us say in between p0 and
*              10**-2, the precision is :
*                 for the positions :
*                 - prec*a0 au for the distances.
*                 - prec rd for the other variables.
*                 for the velocities :
*                 - prec*a0 au/day for the distances.
*                 - prec rd/day for the other variables.
*                   a0 is semi-major axis of the body.
*                 Mercury    a0 =  0.3871 ua
*                 Venus      a0 =  0.7233 ua
*                 Earth      a0 =  1.0000 ua
*                 Mars       a0 =  1.5237 ua
*                 Jupiter    a0 =  5.2026 ua
*                 Saturn     a0 =  9.5547 ua
*                 Uranus     a0 = 19.2181 ua
*                 Neptune    a0 = 30.1096 ua
*
*     lu       logical unit index of the file (integer).
*              The file corresponds to a version of VSOP87 theory and
*              a body, and it must be defined and opened before the
*              first call to subroutine VSOP87.
*
*     Output :
*
*     r(6)     array of the results (real double precision).
*
*              for elliptic coordinates :
*                  1: semi-major axis (au)
*                  2: mean longitude (rd)
*                  3: k = e*cos(pi) (rd)
*                  4: h = e*sin(pi) (rd)
*                  5: q = sin(i/2)*cos(omega) (rd)
*                  6: p = sin(i/2)*sin(omega) (rd)
*                     e:     eccentricity
*                     pi:    perihelion longitude
*                     i:     inclination
*                     omega: ascending node longitude
*
*              for rectangular coordinates :
*                  1: position x (au)
*                  2: position y (au)
*                  3: position z (au)
*                  4: velocity x (au/day)
*                  5: velocity y (au/day)
*                  6: velocity z (au/day)
*
*              for spherical coordinates :
*                  1: longitude (rd)
*                  2: latitude (rd)
*                  3: radius (au)
*                  4: longitude velocity (rd/day)
*                  5: latitude velocity (rd/day)
*                  6: radius velocity (au/day)
*
*     ierr     error index (integer).
*                  0: no error.
*                  1: file error (check up ivers index).
*                  2: file error (check up ibody index).
*                  3: precision error (check up prec parameter).
*                  4: reading file error.
*
*-----------------------------------------------------------------------
*
*     --------------------------------
*     Declarations and initializations
*     --------------------------------
*
      implicit double precision (a-h,o-z)
      character*7 bo,body(0:9)
      dimension r(6),t(-1:5),a0(0:9)
      data body/'SUN','MERCURY','VENUS','EARTH','MARS','JUPITER',
     .          'SATURN','URANUS','NEPTUNE','EMB'/
      data a0/0.01d0,0.3871d0,0.7233d0,1.d0,1.5237d0,5.2026d0,
     .        9.5547d0,19.2181d0,30.1096d0,1.d0/
      data dpi/6.283185307179586d0/
      data t/0.d0,1.d0,5*0.d0/
      data t2000/2451545.d0/
      data a1000/365250.d0/
      k=0
      ideb=0
*
      rewind (lu,err=500)
      do i=1,6
         r(i)=0.d0
      enddo
*
      t(1)=(tdj-t2000)/a1000
      do i=2,5
         t(i)=t(1)*t(i-1)
      enddo
*
      ierr=3
      if (prec.lt.0.d0.or.prec.gt.1.d-2) return
      q=dmax1(3.d0,-dlog10(prec+1.d-50))
*
*     -------------------------------------
*     File reading and substitution of time
*     -------------------------------------
*
100   continue
      read (lu,1001,end=400,err=500) iv,bo,ic,it,in
*
      if (ideb.eq.0) then
         ideb=1
         ierr=1
         if (iv.ne.ivers) return
         ierr=2
         if (bo.ne.body(ibody)) return
         ierr=0
         if (iv.eq.0) k=2
         if (iv.eq.2.or.iv.eq.4) k=1
      endif
*
      if (in.eq.0) goto 100
*
      p=prec/10.d0/(q-2)/(dabs(t(it))+it*dabs(t(it-1))*1.d-4+1.d-50)
      if (k.eq.0.or.(k.ne.0.and.ic.eq.5-2*k)) p=p*a0(ibody)
*
      do 200 n=1,in
         nn=n
         read (lu,1002,err=500) a,b,c
         if (dabs(a).lt.p) goto 300
         u=b+c*t(1)
         cu=dcos(u)
         r(ic)=r(ic)+a*cu*t(it)
         if (iv.eq.0) goto 200
         su=dsin(u)
         r(ic+3)=r(ic+3)+t(it-1)*it*a*cu-t(it)*a*c*su
200   continue
*
      goto 100
300   continue
*
      if (nn.eq.in) goto 100
*
      do n=nn+1,in
         read (lu,1002,err=500)
      enddo
      goto 100
*
400   continue
      if (iv.ne.0) then
         do i=4,6
            r(i)=r(i)/a1000
         enddo
      endif
*
      if (k.eq.0) return
*
      r(k)=dmod(r(k),dpi)
      if (r(k).lt.0.d0) r(k)=r(k)+dpi
      return
*
500   continue
      ierr=4
      return
*
*     -------
*     Formats
*     -------
*
1001  format (17x,i1,4x,a7,12x,i1,17x,i1,i7)
1002  format (79x,f18.11,f14.11,f20.11)
*
      end

max_issues_repo_name: benitocm/my-orbits


text:
!     + + ++ + ++ + ++ + ++ + ++ + ++ + ++ + ++ + ++ + ++ + ++ + ++ + ++ + ++ + ++ + ++ + ++ + ++ + ++ + ++ + ++ + ++ + + + +
!                        ²ú Éú(0,1) ¾ùÔÈ·Ö²¼Ëæ»úÊý×Ó³ÌÐò
!     + + ++ + ++ + ++ + ++ + ++ + ++ + ++ + ++ + ++ + ++ + ++ + ++ + ++ + ++ + ++ + ++ + ++ + ++ + ++ + ++ + ++ + ++ + + + +
      function ran1(idum)
      integer idum,ia,im,iq,ir,ntab,ndiv
       real ran1,am,eps,rnmx
      parameter(ia = 16807,im = 2147483647,am = 1./im,iq = 127773,ir = 2836,&
                ntab = 32,ndiv = 1 + (im-1)/ntab,eps = 1.2e-7,rnmx = 1.-eps)
      integer j,k,iv(ntab),iy
      save iv,iy
      data iv/ntab*0/,iy/0/
      if (idum.le.0.or.iy.eq.0) then
          idum = max(-idum,1)
          do j = ntab + 8,1,-1
             k = idum/iq
             idum = ia*(idum-k*iq)-ir*k
             if (idum.lt.0) idum = idum + im
             if (j.le.ntab) iv(j) = idum
          end do
          iy = iv(1)
      end if
      k = idum/iq
        idum = ia*(idum-k*iq)-ir*k
      if (idum.lt.0) idum = idum + im
      j = 1 + iy/ndiv
      iy = iv(j)
      iv(j) = idum
      ran1 = min(am*iy,rnmx)
      end
!     -----------------------------------------------------------------------

max_stars_repo_name: SpakowitzLab/BasicWLC

max_issues_repo_name: riscalab/wlcsim

max_forks_repo_name: riscalab/wlcsim


text:
! MIT License
!
! Copyright (c) 2010-present David A. Kopriva and other contributors: AUTHORS.md
!
! Permission is hereby granted, free of charge, to any person obtaining a copy
! of this software and associated documentation files (the "Software"), to deal
! in the Software without restriction, including without limitation the rights
! to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
! copies of the Software, and to permit persons to whom the Software is
! furnished to do so, subject to the following conditions:
!
! The above copyright notice and this permission notice shall be included in all
! copies or substantial portions of the Software.
!
! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
! IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
! FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
! AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
! LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
! OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
! SOFTWARE.
!
! HOHQMesh contains code that, to the best of our knowledge, has been released as
! public domain software:
! * `b3hs_hash_key_jenkins`: originally by Rich Townsend,
!    https://groups.google.com/forum/#!topic/comp.lang.fortran/RWoHZFt39ng, 2005
! * `fmin`: originally by George Elmer Forsythe, Michael A. Malcolm, Cleve B. Moler,
!    Computer Methods for Mathematical Computations, 1977
! * `spline`: originally by George Elmer Forsythe, Michael A. Malcolm, Cleve B. Moler,
!    Computer Methods for Mathematical Computations, 1977
! * `seval`: originally by George Elmer Forsythe, Michael A. Malcolm, Cleve B. Moler,
!    Computer Methods for Mathematical Computations, 1977
!
! --- End License
!
!////////////////////////////////////////////////////////////////////////
!
!      SimpleSweep.f90
!      Created: March 28, 2013 9:55 AM
!      By: David Kopriva
!
!     Take a quad mesh generated by the SpecMesh2D code and
!     extrude it vertically in the "z" direction to create a
!     3D Hex mesh or rotate it about the x-axis to generate
!     a volume of revolution.
!
!      SIMPLE_EXTRUSION contains the keys
!         direction
!         height
!         subdivisions
!         start surface name
!         end surface name
!
!      SIMPLE_ROTATION contains the keys
!
!         direction
!         rotation angle factor
!         subdivisions
!         start surface name
!         end surface name
!
!////////////////////////////////////////////////////////////////////////
!
      Module SimpleSweepModule
      USE FTValueDictionaryClass
      USE SMConstants
      USE ProgramGlobals
      USE FTExceptionClass
      USE SharedExceptionManagerModule
      USE HexMeshObjectsModule
      USE ErrorTypesModule
      IMPLICIT NONE
!
!---------------------------------------------------------------
! Define methods to generate a 3D mesh from a 2D one by straight
! line extrusion or simple rotation
!---------------------------------------------------------------
!
      CHARACTER(LEN=STRING_CONSTANT_LENGTH), PARAMETER :: SIMPLE_EXTRUSION_ALGORITHM_KEY    = "SIMPLE_EXTRUSION"
      CHARACTER(LEN=STRING_CONSTANT_LENGTH), PARAMETER :: SIMPLE_ROTATION_ALGORITHM_KEY     = "SIMPLE_ROTATION"
      CHARACTER(LEN=STRING_CONSTANT_LENGTH), PARAMETER :: SIMPLE_EXTRUSION_HEIGHT_KEY       = "height"
      CHARACTER(LEN=STRING_CONSTANT_LENGTH), PARAMETER :: SIMPLE_SWEEP_SUBDIVISIONS_KEY     = "subdivisions"
      CHARACTER(LEN=STRING_CONSTANT_LENGTH), PARAMETER :: SIMPLE_SWEEP_STARTNAME_KEY        = "start surface name"
      CHARACTER(LEN=STRING_CONSTANT_LENGTH), PARAMETER :: SIMPLE_SWEEP_ENDNAME_KEY          = "end surface name"
      CHARACTER(LEN=STRING_CONSTANT_LENGTH), PARAMETER :: SIMPLE_SWEEP_DIRECTION_KEY        = "direction"
      CHARACTER(LEN=STRING_CONSTANT_LENGTH), PARAMETER :: SIMPLE_ROTATION_ANGLE_FRAC_KEY    = "rotation angle factor"
      CHARACTER(LEN=STRING_CONSTANT_LENGTH), PARAMETER :: SIMPLE_SWEEP_PERIODIC_KEY         = "periodic"
      CHARACTER(LEN=STRING_CONSTANT_LENGTH), PARAMETER :: SIMPLE_ROTATION_ANGLE_KEY         = "rotation angle"

      INTEGER, PARAMETER :: SWEEP_FLOOR = 1, SWEEP_CEILING = 2
!

      REAL(KIND=RP), ALLOCATABLE :: chebyPoints(:)
!
!     ------------------------------------------------------------
!     Given the global node ID, get the level and 2D node location
!     in the nodes array
!     ------------------------------------------------------------
!
      INTEGER, ALLOCATABLE :: locAndLevelForNodeID(:,:)

!     ========
      CONTAINS
!     ========
!
!////////////////////////////////////////////////////////////////////////
!
      SUBROUTINE CheckSimpleExtrusionBlock( dict )
!
!        Example block is:
!
!            \begin{SIMPLE_EXTRUSION}
!               direction          = 1 = x, 2 = y, 3 = z
!               height             = 10.0
!               subdivisions       = 5
!               start surface name = "bottom"
!               end surface name   = "top"
!            \end{SIMPLE_EXTRUSION}
!
         IMPLICIT NONE
!
!        ---------
!        Arguments
!        ---------
!
         CLASS(FTValueDictionary) :: dict
!
!        ---------------
!        Local variables
!        ---------------
!
         INTEGER      , EXTERNAL                :: GetIntValue
         REAL(KIND=RP), EXTERNAL                :: GetRealValue
         CHARACTER( LEN=LINE_LENGTH ), EXTERNAL :: GetStringValue
!
!        ---------
!        Direction
!        ---------
!
         IF ( .NOT. dict % containsKey(key = SIMPLE_SWEEP_DIRECTION_KEY) )     THEN
            CALL ThrowErrorExceptionOfType(poster = "CheckSimpleExtrusionBlock", &
                                           msg = "key " // TRIM(SIMPLE_SWEEP_DIRECTION_KEY) // " not found in extrusion block", &
                                           typ = FT_ERROR_FATAL)
         END IF
!
!        ------
!        Height
!        ------
!
         IF ( .NOT. dict % containsKey(key = SIMPLE_EXTRUSION_HEIGHT_KEY) )     THEN
            CALL ThrowErrorExceptionOfType(poster = "CheckSimpleExtrusionBlock", &
                                           msg = "key " // TRIM(SIMPLE_EXTRUSION_HEIGHT_KEY) // " not found in extrusion block", &
                                           typ = FT_ERROR_FATAL)
         END IF
!
!        ------------
!        Subdivisions
!        ------------
!
         IF ( .NOT. dict % containsKey(key = SIMPLE_SWEEP_SUBDIVISIONS_KEY) )     THEN
            CALL ThrowErrorExceptionOfType(poster = "CheckSimpleExtrusionBlock", &
                                           msg = "key " // TRIM(SIMPLE_SWEEP_SUBDIVISIONS_KEY) // " not found in extrusion block", &
                                           typ = FT_ERROR_FATAL)
         END IF
!
!        -------------------
!        Bottom surface name
!        -------------------
!
         IF ( .NOT. dict % containsKey(key = SIMPLE_SWEEP_STARTNAME_KEY) )     THEN
            CALL ThrowErrorExceptionOfType(poster = "CheckSimpleExtrusionBlock", &
                                           msg = "key " // TRIM(SIMPLE_SWEEP_STARTNAME_KEY) // " not found in extrusion block", &
                                           typ = FT_ERROR_FATAL)
         END IF
!
!        ----------------
!        Top surface name
!        ----------------
!
         IF ( .NOT. dict % containsKey(key = SIMPLE_SWEEP_ENDNAME_KEY) )     THEN
            CALL ThrowErrorExceptionOfType(poster = "CheckSimpleExtrusionBlock", &
                                           msg = "key " // TRIM(SIMPLE_SWEEP_ENDNAME_KEY) // " not found in extrusion block", &
                                           typ = FT_ERROR_FATAL)
         END IF

      END SUBROUTINE CheckSimpleExtrusionBlock
!
!////////////////////////////////////////////////////////////////////////
!
      SUBROUTINE CheckSimpleRotationBlock( dict )
!
!        Example block is:
!
!            \begin{SIMPLE_ROTATION}
!               direction               = 1 = x, 2 = y, 3 = z
!               rotation angle factor   = 0.5
!               subdivisions            = 5
!               start surface name      = "bottom"
!               end surface name        = "top"
!            \end{SIMPLE_ROTATION}
!
         IMPLICIT NONE
!
!        ---------
!        Arguments
!        ---------
!
         CLASS(FTValueDictionary) :: dict
!
!        ---------------
!        Local variables
!        ---------------
!
         REAL(KIND=RP)                          :: angleFactor
         INTEGER      , EXTERNAL                :: GetIntValue
         REAL(KIND=RP), EXTERNAL                :: GetRealValue
         CHARACTER( LEN=LINE_LENGTH ), EXTERNAL :: GetStringValue
!
!        ---------
!        Direction
!        ---------
!
         IF ( .NOT. dict % containsKey(key = SIMPLE_SWEEP_DIRECTION_KEY) )     THEN
            CALL ThrowErrorExceptionOfType(poster = "CheckSimpleRotationBlock", &
                                           msg = "key " // TRIM(SIMPLE_SWEEP_DIRECTION_KEY) // " not found in rotation block", &
                                           typ = FT_ERROR_FATAL)
         END IF
!
!        ----------------------
!        Angle (Fraction of PI)
!        ----------------------
!
         IF ( dict % containsKey(key = SIMPLE_ROTATION_ANGLE_FRAC_KEY) )     THEN
            angleFactor = dict % doublePrecisionValueForKey(key = SIMPLE_ROTATION_ANGLE_FRAC_KEY)
            CALL dict % addValueForKey(angleFactor,SIMPLE_ROTATION_ANGLE_KEY)
         ELSE
            CALL ThrowErrorExceptionOfType(poster = "CheckSimpleRotationBlock", &
                                           msg = "key " // TRIM(SIMPLE_ROTATION_ANGLE_FRAC_KEY) // " not found in rotation block", &
                                           typ = FT_ERROR_FATAL)
         END IF
!
!        ------------
!        Subdivisions
!        ------------
!
         IF ( .NOT. dict % containsKey(key = SIMPLE_SWEEP_SUBDIVISIONS_KEY) )     THEN
            CALL ThrowErrorExceptionOfType(poster = "CheckSimpleRotationBlock", &
                                           msg = "key " // TRIM(SIMPLE_SWEEP_SUBDIVISIONS_KEY) // " not found in rotation block", &
                                           typ = FT_ERROR_FATAL)
         END IF
!
!        -------------------
!        Bottom surface name
!        -------------------
!
         IF ( .NOT. dict % containsKey(key = SIMPLE_SWEEP_STARTNAME_KEY) )     THEN
            CALL ThrowErrorExceptionOfType(poster = "CheckSimpleRotationBlock", &
                                           msg = "key " // TRIM(SIMPLE_SWEEP_STARTNAME_KEY) // " not found in rotation block", &
                                           typ = FT_ERROR_FATAL)
         END IF
!
!        ----------------
!        Top surface name
!        ----------------
!
         IF ( .NOT. dict % containsKey(key = SIMPLE_SWEEP_ENDNAME_KEY) )     THEN
            CALL ThrowErrorExceptionOfType(poster = "CheckSimpleRotationBlock", &
                                           msg = "key " // TRIM(SIMPLE_SWEEP_ENDNAME_KEY) // " not found in rotation block", &
                                           typ = FT_ERROR_FATAL)
         END IF

      END SUBROUTINE CheckSimpleRotationBlock
!
!////////////////////////////////////////////////////////////////////////
!
      SUBROUTINE PerformSimpleMeshSweep( project, pMutation, h, parametersDictionary )
         USE MeshProjectClass
         IMPLICIT NONE
!
!        ---------
!        Arguments
!        ---------
!
         TYPE ( MeshProject )        :: project
         CLASS( FTValueDictionary )  :: parametersDictionary
         INTEGER                     :: pMutation
         REAL(KIND=RP)               :: h
!
!        ---------------
!        Local Variables
!        ---------------
!
         TYPE (SMMesh)              , POINTER :: quadMesh
         TYPE (FTMutableObjectArray), POINTER :: quadElementsArray

         INTEGER                :: numberOfLayers
         INTEGER                :: numberOf2DNodes, numberOfQuadElements
         INTEGER                :: numberOfNodes
         INTEGER                :: node2DID
         INTEGER                :: N

         TYPE(SMNodePtr)   , DIMENSION(:), ALLOCATABLE :: quadMeshNodes
         CLASS(SMNode)                   , POINTER     :: currentNode
         CLASS(FTObject)                 , POINTER     :: obj
!
!
         quadMesh              => project % mesh
         N                     =  project % runParams % polynomialOrder
         numberOfQuadElements  =  project % hexMesh % numberOfQuadElements
         numberOfLayers        =  project % hexMesh % numberOfLayers
         numberOf2DNodes       =  quadMesh % nodes % count()
!
!        ---------------------------------------------------------------
!        Make sure that the nodes and elements are consecutively ordered
!        and that the edges refer to the correct elements.
!        ---------------------------------------------------------------
!
         CALL quadMesh % renumberObjects(NODES)
         CALL quadMesh % renumberObjects(ELEMENTS)
         CALL quadMesh % renumberObjects(EDGES)
!
!        ---------------------------------------
!        Gather nodes for easy access
!        TODO: The list class now has a function
!        to return an array of the objects.
!        ---------------------------------------
!
         numberOfNodes   = numberOf2DNodes*(numberOfLayers + 1)
         ALLOCATE( quadMeshNodes(numberOf2DNodes) )

         CALL quadMesh % nodesIterator % setToStart()
         DO WHILE( .NOT.quadMesh % nodesIterator % isAtEnd())

            obj => quadMesh % nodesIterator % object()
            CALL castToSMNode(obj,currentNode)
            node2DID = currentNode % id
            quadMeshNodes(node2DID) % node => currentNode

            CALL quadMesh % nodesIterator % moveToNext()
         END DO
!
!        ---------------------------------------------------------------
!        Allocate connections between global ID and local (2D, level) id
!        ---------------------------------------------------------------
!
         ALLOCATE( locAndLevelForNodeID(2, numberOfNodes) )
!
!        ------------------------------
!        Sweep the skeleton of the mesh
!        ------------------------------
!
         CALL sweepNodes( quadMeshNodes, project % hexMesh, h, pMutation )
         CALL sweepElements( quadMesh, project % hexMesh, &
                             numberofLayers, parametersDictionary )
!
!        -------------------------------------
!        Sweep the internal degrees of freedom
!        -------------------------------------
!
         quadElementsArray => quadMesh % elements % allObjects()
         CALL SweepInternalDOFs(hex8Mesh          = project % hexMesh, &
                                quadElementsArray = quadElementsArray, &
                                N                 = N,                 &
                                h                 = h,                &
                                pmutation         = pMutation)

         CALL releaseFTMutableObjectArray(quadElementsArray)
         DEALLOCATE(quadMeshNodes)
!

      END SUBROUTINE PerformSimpleMeshSweep
!
!////////////////////////////////////////////////////////////////////////
!
      SUBROUTINE sweepNodes( quadMeshNodes, hex8Mesh, h, pMutation )
         IMPLICIT NONE
!
!        ---------
!        Arguments
!        ---------
!
         TYPE ( StructuredHexMesh )    :: hex8Mesh
         REAL(KIND=RP)                 :: h
         TYPE(SMNodePtr), DIMENSION(:) :: quadMeshNodes
         INTEGER                       :: pMutation
!
!        ---------------
!        Local Variables
!        ---------------
!
         INTEGER         :: numberOf2DNodes
         INTEGER         :: nodeID
         INTEGER         :: j, k
         INTEGER         :: numberOfLayers
         REAL(KIND=RP)   :: z, xi
!
!        ---------------------------------------
!        Generate the new nodes for the hex mesh
!        layer by layer. Order the new node IDs
!        layer by layer, too.
!        ---------------------------------------
!
         numberOf2DNodes = SIZE(quadMeshNodes)
         numberOfLayers  = hex8Mesh % numberofLayers

         nodeID = 1
         DO j = 0, numberofLayers

            DO k = 1, numberOf2DNodes
               hex8Mesh % nodes(k,j) % globalID = nodeID
!
!              ----------------------------------
!              Interpolate between top and bottom
!              ----------------------------------
!
               xi = DBLE(j)/DBLE(numberOfLayers)
               z  = quadMeshNodes(k) % node % x(pMutation)*(1.0_RP - xi) + h*xi
               hex8Mesh % nodes(k,j) % x  = extrudedNodeLocation(baseLocation = quadMeshNodes(k) % node % x, &
                                                                 delta        = z, &
                                                                 pmutation    = pMutation)
               locAndLevelForNodeID(1,nodeID) = k
               locAndLevelForNodeID(2,nodeID) = j
               nodeID = nodeID + 1
            END DO

         END DO
!
      END SUBROUTINE sweepNodes
!
!////////////////////////////////////////////////////////////////////////
!
      SUBROUTINE SweepInternalDOFs( hex8Mesh, quadElementsArray, N, h, pmutation)
         USE SMMeshClass
         USE FTMutableObjectArrayClass
         IMPLICIT NONE
!
!        ---------
!        Arguments
!        ---------
!
         TYPE (FTMutableObjectArray), POINTER  :: quadElementsArray
         TYPE ( StructuredHexMesh )            :: hex8Mesh
         REAL(KIND=RP)                         :: h
         INTEGER                               :: pMutation
         INTEGER                               :: N
!
!        ---------------
!        Local Variables
!        ---------------
!
         INTEGER                   :: l, m, i, j, k, nLev
         REAL(KIND=RP)             :: x(3), z, y(3)
         REAL(KIND=RP)             :: delta, eta, xi
         CLASS(FTObject) , POINTER :: obj
         CLASS(SMElement), POINTER :: e
!
!        ------------------------------------------
!        Extend the face points on the quad element
!        up through the hex element.
!        ------------------------------------------
!
         nLev = hex8Mesh % numberofLayers

         DO l = 1, hex8Mesh % numberOfQuadElements
            obj => quadElementsArray % objectAtIndex(l)
            CALL castToSMelement(obj,e)
            DO m = 1, hex8Mesh % numberofLayers
               xi = DBLE(m-1)/DBLE(nLev)
               DO k = 0, N
                  eta = 0.5_RP*(1.0 - COS(k*PI/N))
                  DO j = 0, N
                     DO i = 0, N
                        delta = (h - e % xPatch(pMutation,i,j))/nLev
                        z     = xi*h + (1.0_RP - xi)*e % xPatch(pMutation,i,j) + delta*eta
                        x     = permutePosition(x = e % xPatch(:,i,j),pmutation = pMutation)
                        y     = extrudedNodeLocation(baseLocation = x, delta = z, pmutation = pMutation)

                        hex8Mesh % elements(l,m) % x(:,i,j,k) = y
                     END DO
                  END DO
               END DO
            END DO
         END DO

      END SUBROUTINE SweepInternalDOFs
!
!////////////////////////////////////////////////////////////////////////
!
      SUBROUTINE sweepElements( quadMesh, hex8Mesh, numberofLayers, parametersDictionary )
!
!        -------------------------------
!        Call after generating the nodes
!        -------------------------------
!
         USE MeshProjectClass
         IMPLICIT NONE
!
!        ---------
!        Arguments
!        ---------
!
         TYPE ( SMMesh )             :: quadMesh
         TYPE ( StructuredHexMesh )  :: hex8Mesh
         INTEGER                     :: numberOfLayers
         TYPE( FTValueDictionary)    :: parametersDictionary
!
!        ---------------
!        Local Variables
!        ---------------
!
         INTEGER   :: numberOfQuadElements
         INTEGER   :: elementID, nodeID, node2DID, quadElementID
         INTEGER   :: j, k
         INTEGER   :: flagMap(4) = [1,4,2,6]

         CLASS(SMNode)                   , POINTER     :: node
         CLASS(SMElement)                , POINTER     :: currentQuadElement
         CLASS(FTObject)                 , POINTER     :: obj

         numberOfQuadElements = hex8Mesh % numberOfQuadElements
!
!        ---------------------------------
!        Build the elements layer by layer
!        ---------------------------------
!
         elementID = 1
         DO j = 1, numberOfLayers
            quadElementID = 1

            CALL quadMesh % elementsIterator % setToStart()

            DO WHILE( .NOT. quadMesh % elementsIterator % isAtEnd() )
               obj => quadMesh % elementsIterator % object()
               CALL castToSMElement(obj,currentQuadElement)
!
!              -----------------------
!              Set the element nodeIDs
!              -----------------------
!
               DO k = 1, 4
!
!                 -------------
!                 Bottom of hex
!                 -------------
!
                  obj => currentQuadElement % nodes % objectAtIndex(k)
                  CALL cast(obj,node)
                  node2DID = node % id
                  nodeID   = hex8Mesh % nodes(node2DID,j-1) % globalID
                  hex8Mesh % elements(quadElementID,j) % nodeIDs(k) = nodeID
!
!                 ----------
!                 Top of hex
!                 ----------
!
                  nodeID = hex8Mesh % nodes(node2DID,j) % globalID
                  hex8Mesh % elements(quadElementID,j)  % nodeIDs(k+4) = nodeID

               END DO
!
!              ------------------------------------------------------------------
!              Set boundary condition names at the start and end of the extrusion
!              as defined in the control file
!              ------------------------------------------------------------------
!
               IF ( j == 1 )     THEN
                  hex8Mesh % elements(quadElementID,j) % bFaceName(3) = &
                  parametersDictionary % stringValueForKey(key             = SIMPLE_SWEEP_STARTNAME_KEY,&
                                                           requestedLength = LINE_LENGTH)
               END IF
               IF (j == numberOfLayers)     THEN
                  hex8Mesh % elements(quadElementID,j) % bFaceName(5) = &
                  parametersDictionary % stringValueForKey(key             = SIMPLE_SWEEP_ENDNAME_KEY,&
                                                           requestedLength = LINE_LENGTH)
               END IF
!
!              ----------------------------------------------------------------
!              Use edge info of parent quad element to set boundary curve flags
!              and names for the new hex element
!              ----------------------------------------------------------------
!
               DO k = 1, 4
                  IF ( currentQuadElement % boundaryInfo % bCurveFlag(k) == ON )     THEN
                     hex8Mesh % elements(quadElementID,j) % bFaceFlag(flagMap(k)) = ON
                     hex8Mesh % elements(quadElementID,j) % bFaceFlag(3) = ON
                     hex8Mesh % elements(quadElementID,j) % bFaceFlag(5) = ON
                  END IF
                  hex8Mesh % elements(quadElementID,j) % bFaceName(flagMap(k)) &
                                      = currentQuadElement % boundaryInfo % bCurveName(k)
               END DO

               quadElementID  = quadElementID + 1
               elementID      = elementID + 1

               CALL quadMesh % elementsIterator % moveToNext()
            END DO

         END DO
!
      END SUBROUTINE sweepElements
!
!////////////////////////////////////////////////////////////////////////
!
      FUNCTION extrudedNodeLocation(baseLocation,delta,pmutation)  RESULT(x)
         IMPLICIT NONE
         REAL(KIND=RP) :: baseLocation(3), delta
         INTEGER       :: pmutation
         REAL(KIND=RP) :: x(3)

         x              = baseLocation
         x(pmutation)   = delta
      END FUNCTION extrudedNodeLocation
!
!////////////////////////////////////////////////////////////////////////
!
      FUNCTION rotatedNodeLocation(baseLocation,theta,rotAxis)  RESULT(x)
         IMPLICIT NONE
         REAL(KIND=RP) :: baseLocation(3), theta
         INTEGER       :: rotAxis
         REAL(KIND=RP) :: x(3)
         REAL(KIND=RP) :: r

         x              = baseLocation
         SELECT CASE ( rotAxis )
            CASE( 1 ) ! rotation about x-Axis
               r    = baseLocation(2)
               x(2) = r*COS(theta)
               x(3) = r*SIN(theta)
            CASE (2)  ! rotation about y-Axix
               r    = baseLocation(1)
               x(1) = r*COS(theta)
               x(3) = r*SIN(theta)
            CASE (3)  ! rotation about z-Axis
               r    = baseLocation(2)
               x(2) = r*COS(theta)
               x(1) = r*SIN(theta)
            CASE DEFAULT
         END SELECT

      END FUNCTION rotatedNodeLocation
!
!////////////////////////////////////////////////////////////////////////
!
      FUNCTION permutePosition(x, pmutation)  RESULT(y)
         IMPLICIT NONE
         REAL(KIND=RP), DIMENSION(3) :: x, y
         INTEGER                     :: pmutation

         y  = CSHIFT(x, SHIFT = -pmutation)

      END FUNCTION permutePosition
!
!////////////////////////////////////////////////////////////////////////
!
      SUBROUTINE RotateAll(mesh, N, rotAxis)
        IMPLICIT NONE
!
!       ----------
!       Arguments
!       ----------
!
        TYPE(StructuredHexMesh) :: mesh
        INTEGER                 :: rotAxis
        INTEGER                 :: N
!
!       ---------------
!       Local Variables
!       ---------------
!
         INTEGER       :: rotMap(3) = [3, 3, 1]
         INTEGER       :: k, j, i, l, m
         REAL(KIND=RP) :: rotated(3)
!
!        ----------------
!        Rotate the nodes
!        ----------------
!
         DO l = 0, SIZE(mesh % nodes,2)-1
            DO m = 1, SIZE(mesh % nodes,1)
               rotated = rotatedNodeLocation(baseLocation = mesh % nodes(m,l) % x,            &
                                             theta        = mesh % nodes(m,l) % x(rotMap(rotAxis)), &
                                             rotAxis    = rotAxis)
               mesh % nodes(m,l) % x = rotated
            END DO
         END DO
!
!        -----------------------------------
!        Rotate the internal DOFs
!        With rotation, all faces are curved
!        -----------------------------------
!
         DO l = 1, mesh % numberofLayers             ! level
            DO m = 1, mesh % numberOfQuadElements    ! element on original quad mesh
               mesh % elements(m,l) % bFaceFlag = ON
               DO k = 0, N
                  DO j = 0, N
                     DO i = 0, N
                        rotated = rotatedNodeLocation(baseLocation = mesh % elements(m,l) % x(:,i,j,k), &
                                                      theta        = mesh % elements(m,l) % x((rotMap(rotAxis)),i,j,k),  &
                                                      rotAxis      = rotAxis)
                        mesh % elements(m,l) % x(:,i,j,k) = rotated
                     END DO
                  END DO
               END DO
            END DO
         END DO

      END SUBROUTINE RotateAll

      END Module SimpleSweepModule

max_stars_repo_name: HigherOrderMethods/HOHQMesh

max_issues_repo_name: DavidAKopriva/HOHQMesh

max_forks_repo_name: DavidAKopriva/HOHQMesh


! Copyright (c) 2011-2013 Bartosz Szczesny
! LICENSE: The MIT License (MIT)

program task_4
    integer, parameter :: nsites = 5, nhours = 24
    integer :: site, hour, site_min, hour_min
    character( len = nsites ) :: abcde = 'ABCDE'
    real :: temp( nsites, nhours ), humi( nsites, nhours )
    real :: temp_avg, humi_avg

    ! read input data from two files
    open( unit = 22, file = "temp.txt" )
    open( unit = 33, file = "humi.txt" )
    do site = 1, nsites
        do hour = 1, nhours
            read( 22, * ) temp( site, hour )
            read( 33, * ) humi( site, hour )
        enddo
    enddo
    close( 22 )
    close( 33 )
    do hour = 1, nhours
         print *, ( temp( site, hour ), site = 1, nsites )
    enddo

    ! find where minimum humidity occurred
    site_min = 1
    hour_min = 1
    do site = 1, nsites
        do hour = 1, nhours
            if ( humi( site, hour ) < humi( site_min, hour_min ) ) then
                site_min = site
                hour_min = hour
            endif
        enddo
    enddo
    print *, 'min humidity:'
    print *, 'site = ', abcde( site_min : site_min )
    print *, 'hour =', hour_min
    print *, 'temp =', temp( site_min, hour_min )
    print *, 'humi =', humi( site_min, hour_min )

    ! tempi_avg over all sites at hour = 14
    temp_avg = 0.0
    do site = 1, nsites
        temp_avg = temp_avg + temp( site, 14 )
    enddo
    temp_avg = temp_avg / real( nsites )
    print *, '@ hour = 14, temp_avg =', temp_avg

    ! humi_avg at site C over 24 hours
    site = 3 ! site C
    temp_avg = 0.0
    do hour = 1, nhours
        humi_avg = humi_avg + humi( site, hour )
    enddo
    humi_avg = humi_avg / real( nhours )
    print *, '@ site = C, humi_avg =', humi_avg

    stop
end


text:
Describe Users/ToddFreidman here.

20100413 16:37:28 nbsp Welcome to the Wiki. Please dont remove content. Users/JasonAller


text:
The May 2011 Election had only one item on the ballot:  Measure A, which asks, Recognizing cuts in state funding, shall the Davis Joint Unified School District be authorized to levy an emergency special tax for a period of two years not to exceed the annual rate of $20.00 per dwelling unit for multidwelling parcels and $200.00 per parcel for all other parcels, to fund essential school programs including core subjects, elective classes, reduced class sizes, counselors and school site safety?

This was Davis first MailOnly Elections mailonly election.  The ballots were to be mailed to local voters on April 4 and are due back by 8 p.m., May 3.  The measure required 2/3 approval in order to pass.

As of http://www.yoloelections.org/returns/ 9:52 PM, May 3 (unlikely to change):

See also http://twitter.com/#!/yoloelections twitter source.

Bob Dunning had http://www.davisenterprise.com/opinion/dunning/besuretovotebutdontwastetimelookingforpollingplaces/ predicted on an uncluttered ballot with only one choice to be made, Measure A will carry the day with 70.2 percent of the vote … but don’t bet the mortgage on it

Columnist Rich Rifkin had predicted a 67.2% margin.


text:
!!  Copyright (C)  Stichting Deltares, 2012-2016.
!!
!!  This program is free software: you can redistribute it and/or modify
!!  it under the terms of the GNU General Public License version 3,
!!  as published by the Free Software Foundation.
!!
!!  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 registered trademarks
!!  of Stichting Deltares remain the property of Stichting Deltares. All
!!  rights reserved.

      subroutine expand_frc ( lunrep, procesdef, allitems, sfracs)

      ! expand processes per fractions

      use processet
      use timers       !   performance timers

      implicit none

      ! declaration of arguments

      integer                   :: lunrep          ! report file
      type(procespropcoll)      :: procesdef       ! the process definition
      type(itempropcoll)        :: allitems        ! all items of the proces system
      type(sfracsprop)          :: sfracs          ! substance fraction properties

      ! local decalarations

      type(procesprop), pointer :: proc            ! single process
      type(ioitemprop), pointer :: new_item(:)     ! expanded list of items
      type(stochiprop), pointer :: new_stochi(:)   ! expanded list of items
      type(itemprop)            :: item            ! one item
      integer                   :: nproc           ! number of processes
      integer                   :: iproc           ! loop counter processes
      integer                   :: isfrac          ! index substance fractions
      integer                   :: isfrac2         ! index substance fractions
      integer                   :: nfrac           ! number fractions in substance fraction
      character(len=20)         :: basnam          ! base name substance fractions
      character(len=20)         :: fracnam         ! name of substance fraction
      integer                   :: i_item          ! index input item
      integer                   :: i_flux          ! index flux item
      integer                   :: i_stochi        ! index stochi
      integer                   :: indx            ! index in list
      integer                   :: iret            ! index in collection
      integer                   :: ifrac           ! fraction number
      character(len=3)          :: suffix          ! suffix
      integer                   :: ierr_alloc      ! error indication

      logical                   :: l_expand              ! expand item list
      logical                   :: l_frac(sfracs%nsfrac) ! fraction involved in process

      integer                   :: no_input_new          ! new number of items
      integer                   :: no_output_new         ! new number of items
      integer                   :: no_fluxoutput_new     ! new number of items
      integer                   :: no_fluxstochi_new     ! new number of items
      integer                   :: no_velostochi_new     ! new number of items
      integer                   :: no_dispstochi_new     ! new number of items
      integer                   :: i_new                 ! index in new items
      integer                   :: isfrac_found          ! index substance fractions
      integer                   :: i_star                ! index of * in name
      integer                   :: nzoek                 ! nzoek
      integer(4) :: ithndl = 0
      if (timon) call timstrt( "expand_frc", ithndl )

      ! loop over the processes

      nproc = procesdef%cursize
      do iproc = 1, nproc

         proc => procesdef%procesprops(iproc)

         ! check for expansion and for the fractions involved

         l_expand = .false.
         l_frac   = .false.
         do i_item = 1, proc%no_input
            if ( proc%input_item(i_item)%name(1:1) .eq. '*' ) then
               do isfrac = 1, sfracs%nsfrac
                  basnam = sfracs%name(isfrac)
                  call zoek( basnam, 1, proc%input_item(i_item)%name(2:), 19, indx)
                  if ( indx .gt. 0 ) then
                     l_expand       = .true.
                     l_frac(isfrac) = .true.
                     exit
                  endif
               enddo
            endif
         enddo

         if ( l_expand ) then

            ! count the dimensions of the expanded process

            no_input_new = proc%no_input
            do i_item = 1, proc%no_input

               ! check for the fractions itself

               do isfrac = 1, sfracs%nsfrac
                  if ( l_frac(isfrac) ) then
                     call zoek( sfracs%name(isfrac), 1, proc%input_item(i_item)%name, 20, indx)
                     if ( indx .gt. 0 ) then
                        no_input_new = no_input_new + sfracs%nfrac(isfrac) - 1
                        exit
                     endif
                  endif
               enddo

               ! check for input which needs to be expanded

               i_star = index(proc%input_item(i_item)%name,'*')
               if ( i_star .gt. 1 ) then
                  do isfrac = 1, sfracs%nsfrac
                     nzoek = 20-i_star
                     call zoek( sfracs%name(isfrac), 1, proc%input_item(i_item)%name(i_star+1:),  nzoek, indx)
                     if ( indx .gt. 0 ) then
                        no_input_new = no_input_new + sfracs%nfrac(isfrac) - 1
                        exit
                     endif
                  enddo
               endif

            enddo

            no_output_new = proc%no_output
            do i_item = 1, proc%no_output
               i_star = index(proc%output_item(i_item)%name,'*')
               if ( i_star .gt. 0 ) then
                  do isfrac = 1, sfracs%nsfrac
                     nzoek = 20-i_star
                     call zoek( sfracs%name(isfrac), 1, proc%output_item(i_item)%name(i_star+1:),  nzoek, indx)
                     if ( indx .gt. 0 ) then
                        no_output_new = no_output_new + sfracs%nfrac(isfrac) - 1
                        exit
                     endif
                  enddo
               endif
            enddo

            no_fluxoutput_new   = proc%no_fluxoutput
            do i_flux = 1, proc%no_fluxoutput
               i_star = index(proc%fluxoutput(i_flux)%name,'*')
               if ( i_star .gt. 0 ) then
                  do isfrac = 1, sfracs%nsfrac
                     nzoek = 20-i_star
                     call zoek( sfracs%name(isfrac), 1, proc%fluxoutput(i_flux)%name(i_star+1:),  nzoek, indx)
                     if ( indx .gt. 0 ) then
                        no_fluxoutput_new = no_fluxoutput_new + sfracs%nfrac(isfrac) - 1
                        exit
                     endif
                  enddo
               endif
            enddo

            no_fluxstochi_new = proc%no_fluxstochi
            do i_stochi = 1, proc%no_fluxstochi
               i_star = index(proc%fluxstochi(i_stochi)%ioitem,'*')
               if ( i_star .gt. 0 ) then
                  do isfrac = 1, sfracs%nsfrac
                     nfrac  = sfracs%nfrac(isfrac)
                     basnam = sfracs%name(isfrac)
                     nzoek = 20-i_star
                     call zoek( basnam, 1, proc%fluxstochi(i_stochi)%ioitem(i_star+1:),  nzoek, indx)
                     if ( indx .gt. 0 ) then
                        no_fluxstochi_new = no_fluxstochi_new + nfrac - 1
                        exit
                     endif
                  enddo
               endif
            enddo

            no_velostochi_new = proc%no_velostochi
            do i_item = 1, proc%no_velostochi
               i_star = index(proc%velostochi(i_item)%ioitem,'*')
               if ( i_star .gt. 0 ) then
                  do isfrac = 1, sfracs%nsfrac
                     nzoek = 20-i_star
                     call zoek( sfracs%name(isfrac), 1, proc%velostochi(i_item)%ioitem(i_star+1:),  nzoek, indx)
                     if ( indx .gt. 0 ) then
                        no_velostochi_new = no_velostochi_new + sfracs%nfrac(isfrac) - 1
                        exit
                     endif
                  enddo
               endif
            enddo

            no_dispstochi_new = proc%no_dispstochi
            do i_item = 1, proc%no_dispstochi
               i_star = index(proc%dispstochi(i_item)%ioitem,'*')
               if ( i_star .gt. 0 ) then
                  do isfrac = 1, sfracs%nsfrac
                     nzoek = 20-i_star
                     call zoek( sfracs%name(isfrac), 1, proc%dispstochi(i_item)%ioitem(i_star+1:),  nzoek, indx)
                     if ( indx .gt. 0 ) then
                        no_dispstochi_new = no_dispstochi_new + sfracs%nfrac(isfrac) - 1
                        exit
                     endif
                  enddo
               endif
            enddo

            ! then expand

            write(lunrep,2000) ' expanding fraction input for process [',proc%name,']'

            allocate(new_item(no_input_new))
            i_new = 0
            do i_item = 1, proc%no_input

               ! check for fraction (i_star = 0 then the substance fraction itself, i_star > 1 then process parameter)

               isfrac_found = 0
               i_star = index(proc%input_item(i_item)%name,'*')
               do isfrac = 1, sfracs%nsfrac
                  if ( l_frac(isfrac) ) then
                     nzoek = 20-i_star
                     call zoek( sfracs%name(isfrac), 1, proc%input_item(i_item)%name(i_star+1:),  nzoek, indx)
                     if ( indx .gt. 0 ) then
                        isfrac_found = isfrac
                        exit
                     endif
                  endif
               enddo

               if ( i_star .eq. 1 ) then

                  ! number of fractions

                  i_new = i_new + 1
                  new_item(i_new) = proc%input_item(i_item)
                  new_item(i_new)%indx = i_new
                  if ( isfrac_found .gt. 0 ) then
                     new_item(i_new)%actdef = sfracs%nfrac(isfrac_found)
                  endif

               elseif ( isfrac_found .gt. 0 ) then

                  ! expand input

                  nfrac  = sfracs%nfrac(isfrac_found)
                  basnam = sfracs%name(isfrac_found)
                  do ifrac = 1, nfrac
                     if ( ifrac .lt. 100 ) then
                        write(suffix,'(i2.2)') ifrac
                     else
                        write(suffix,'(i3.3)') ifrac
                     endif
                     fracnam             = trim(basnam)//trim(suffix)

                     i_new = i_new + 1
                     new_item(i_new) = proc%input_item(i_item)
                     new_item(i_new)%indx = i_new
                     if ( i_star .eq. 0 ) then
                        new_item(i_new)%name = fracnam
                     else
                        new_item(i_new)%name = trim(proc%input_item(i_item)%name(1:i_star))//fracnam
                     endif
                     item%name = new_item(i_new)%name
                     iret      = itempropcollfind( allitems, item )
                     if ( iret .le. 0 ) then
                        item%text    = proc%input_item(i_item)%item%text
                        item%default = proc%input_item(i_item)%item%default
                        item%waqtype = proc%input_item(i_item)%item%waqtype
                        iret         = itempropcolladd( allitems, item )
                     endif
                     new_item(i_new)%item=>allitems%itemproppnts(iret)%pnt
                  enddo
               else

                  ! just copy input

                  i_new = i_new + 1
                  new_item(i_new) = proc%input_item(i_item)
                  new_item(i_new)%indx = i_new

               endif

            enddo
            if (proc%no_input .gt. 0) deallocate(proc%input_item)
            proc%input_item => new_item
            proc%no_input   =  no_input_new

            ! output

            i_new   = 0
            allocate(new_item(no_output_new))
            do i_item = 1, proc%no_output
               isfrac_found = 0
               i_star = index(proc%output_item(i_item)%name,'*')
               if ( i_star .gt. 0 ) then
                  do isfrac = 1, sfracs%nsfrac
                     nzoek = 20-i_star
                     call zoek( sfracs%name(isfrac), 1, proc%output_item(i_item)%name(i_star+1:),  nzoek, indx)
                     if ( indx .gt. 0 ) then
                        isfrac_found = isfrac
                        exit
                     endif
                  enddo
               endif
               if ( isfrac_found .gt. 0 ) then
                  nfrac  = sfracs%nfrac(isfrac_found)
                  basnam = sfracs%name(isfrac_found)
                  do ifrac = 1, nfrac
                     if ( ifrac .lt. 100 ) then
                        write(suffix,'(i2.2)') ifrac
                     else
                        write(suffix,'(i3.3)') ifrac
                     endif
                     fracnam             = trim(basnam)//trim(suffix)
                     i_new = i_new + 1
                     new_item(i_new)      = proc%output_item(i_item)
                     new_item(i_new)%indx = i_new
                     new_item(i_new)%name = trim(proc%output_item(i_item)%name(1:i_star))//fracnam
                     item%name = new_item(i_new)%name
                     iret      = itempropcollfind( allitems, item )
                     if ( iret .le. 0 ) then
                        item%text    = proc%output_item(i_item)%item%text
                        item%default = proc%output_item(i_item)%item%default
                        item%waqtype = proc%output_item(i_item)%item%waqtype
                        iret         = itempropcolladd( allitems, item )
                     endif
                     new_item(i_new)%item=>allitems%itemproppnts(iret)%pnt
                  enddo
               else
                  i_new = i_new + 1
                  new_item(i_new) = proc%output_item(i_item)
                  new_item(i_new)%indx = i_new
               endif
            enddo
            if (proc%no_output .gt. 0) deallocate(proc%output_item)
            proc%output_item    => new_item
            proc%no_output = no_output_new

            ! fluxes

            i_new   = 0
            allocate(new_item(no_fluxoutput_new))
            do i_item = 1, proc%no_fluxoutput
               isfrac_found = 0
               i_star = index(proc%fluxoutput(i_item)%name,'*')
               if ( i_star .gt. 0 ) then
                  do isfrac = 1, sfracs%nsfrac
                     nzoek = 20-i_star
                     call zoek( sfracs%name(isfrac), 1, proc%fluxoutput(i_item)%name(i_star+1:),  nzoek, indx)
                     if ( indx .gt. 0 ) then
                        isfrac_found = isfrac
                        exit
                     endif
                  enddo
               endif
               if ( isfrac_found .gt. 0 ) then
                  nfrac  = sfracs%nfrac(isfrac_found)
                  basnam = sfracs%name(isfrac_found)
                  do ifrac = 1, nfrac
                     if ( ifrac .lt. 100 ) then
                        write(suffix,'(i2.2)') ifrac
                     else
                        write(suffix,'(i3.3)') ifrac
                     endif
                     fracnam             = trim(basnam)//trim(suffix)
                     i_new = i_new + 1
                     new_item(i_new)      = proc%fluxoutput(i_item)
                     new_item(i_new)%indx = i_new
                     new_item(i_new)%name = trim(proc%fluxoutput(i_item)%name(1:i_star))//fracnam
                     item%name = new_item(i_new)%name
                     iret      = itempropcollfind( allitems, item )
                     if ( iret .le. 0 ) then
                        item%text    = proc%fluxoutput(i_item)%item%text
                        item%default = proc%fluxoutput(i_item)%item%default
                        item%waqtype = proc%fluxoutput(i_item)%item%waqtype
                        iret         = itempropcolladd( allitems, item )
                     endif
                     new_item(i_new)%item=>allitems%itemproppnts(iret)%pnt
                  enddo
               else
                  i_new = i_new + 1
                  new_item(i_new) = proc%fluxoutput(i_item)
                  new_item(i_new)%indx = i_new
               endif
            enddo
            if (proc%no_fluxoutput .gt. 0) deallocate(proc%fluxoutput)
            proc%fluxoutput    => new_item
            proc%no_fluxoutput = no_fluxoutput_new

            ! flux stochi

            i_new   = 0
            allocate(new_stochi(no_fluxstochi_new))
            do i_stochi = 1, proc%no_fluxstochi
               i_star = index(proc%fluxstochi(i_stochi)%ioitem,'*')
               if ( i_star .gt. 0 ) then
                  do isfrac = 1, sfracs%nsfrac
                     nfrac  = sfracs%nfrac(isfrac)
                     basnam = sfracs%name(isfrac)
                     nzoek = 20-i_star
                     call zoek( basnam, 1, proc%fluxstochi(i_stochi)%ioitem(i_star+1:),  nzoek, indx)
                     if ( indx .gt. 0 ) then
                        isfrac_found = isfrac
                        exit
                     endif
                  enddo
                  if ( isfrac_found .gt. 0 ) then
                     do ifrac = 1, nfrac
                        if ( ifrac .lt. 100 ) then
                           write(suffix,'(i2.2)') ifrac
                        else
                           write(suffix,'(i3.3)') ifrac
                        endif
                        fracnam = trim(basnam)//trim(suffix)
                        i_new   = i_new + 1
                        new_stochi(i_new) = proc%fluxstochi(i_stochi)
                        new_stochi(i_new)%ioitem = trim(proc%fluxstochi(i_stochi)%ioitem(1:i_star))//fracnam
                        call zoek (basnam,1,proc%fluxstochi(i_stochi)%substance,20,indx)
                        if ( indx .gt. 0 ) then
                           new_stochi(i_new)%substance = fracnam
                        else
                           ! look for linked substance (if not found the substance name stays the same
                           do isfrac2 = 1, nfrac
                              call zoek (sfracs%name(isfrac2),1,proc%fluxstochi(i_stochi)%substance,20,indx)
                              if ( indx .gt. 0 .and. sfracs%linklist(isfrac,isfrac2) .eq. 1 ) then
                                 new_stochi(i_new)%substance = trim(sfracs%name(isfrac2))//suffix
                                 exit
                              endif
                           enddo
                        endif
                     enddo
                     isfrac_found = 0
                  else
                     i_new = i_new + 1
                     new_stochi(i_new) = proc%fluxstochi(i_stochi)
                  endif
               else
                  i_new = i_new + 1
                  new_stochi(i_new) = proc%fluxstochi(i_stochi)
               endif
            enddo
            if (proc%no_fluxstochi .gt. 0) deallocate(proc%fluxstochi)
            proc%fluxstochi    => new_stochi
            proc%no_fluxstochi = no_fluxstochi_new

            ! disp stochi

            i_new   = 0
            allocate(new_stochi(no_dispstochi_new))
            do i_stochi = 1, proc%no_dispstochi
               i_star = index(proc%dispstochi(i_stochi)%ioitem,'*')
               if ( i_star .gt. 0 ) then
                  do isfrac = 1, sfracs%nsfrac
                     nfrac  = sfracs%nfrac(isfrac)
                     basnam = sfracs%name(isfrac)
                     nzoek = 20-i_star
                     call zoek( basnam, 1, proc%dispstochi(i_stochi)%ioitem(i_star+1:),  nzoek, indx)
                     if ( indx .gt. 0 ) then
                        isfrac_found = isfrac
                        exit
                     endif
                  enddo
                  if ( isfrac_found .gt. 0 ) then
                     do ifrac = 1, nfrac
                        if ( ifrac .lt. 100 ) then
                           write(suffix,'(i2.2)') ifrac
                        else
                           write(suffix,'(i3.3)') ifrac
                        endif
                        fracnam = trim(basnam)//trim(suffix)
                        i_new   = i_new + 1
                        new_stochi(i_new) = proc%dispstochi(i_stochi)
                        new_stochi(i_new)%ioitem = trim(proc%dispstochi(i_stochi)%ioitem(1:i_star))//fracnam
                        call zoek (basnam,1,proc%dispstochi(i_stochi)%substance,20,indx)
                        if ( indx .gt. 0 ) then
                           new_stochi(i_new)%substance = fracnam
                        else
                           ! look for linked substance (if not found the substance name stays the same
                           do isfrac2 = 1, nfrac
                              call zoek (sfracs%name(isfrac2),1,proc%dispstochi(i_stochi)%substance,20,indx)
                              if ( indx .gt. 0 .and. sfracs%linklist(isfrac,isfrac2) .eq. 1 ) then
                                 new_stochi(i_new)%substance = trim(sfracs%name(isfrac2))//suffix
                                 exit
                              endif
                           enddo
                        endif
                     enddo
                     isfrac_found = 0
                  else
                     i_new = i_new + 1
                     new_stochi(i_new) = proc%dispstochi(i_stochi)
                  endif
               else
                  i_new = i_new + 1
                  new_stochi(i_new) = proc%dispstochi(i_stochi)
               endif
            enddo
            if (proc%no_dispstochi .gt. 0) deallocate(proc%dispstochi)
            proc%dispstochi    => new_stochi
            proc%no_dispstochi = no_dispstochi_new

            ! velo stochi

            i_new   = 0
            allocate(new_stochi(no_velostochi_new))
            do i_stochi = 1, proc%no_velostochi
               i_star = index(proc%velostochi(i_stochi)%ioitem,'*')
               if ( i_star .gt. 0 ) then
                  do isfrac = 1, sfracs%nsfrac
                     nfrac  = sfracs%nfrac(isfrac)
                     basnam = sfracs%name(isfrac)
                     nzoek = 20-i_star
                     call zoek( basnam, 1, proc%velostochi(i_stochi)%ioitem(i_star+1:),  nzoek, indx)
                     if ( indx .gt. 0 ) then
                        isfrac_found = isfrac
                        exit
                     endif
                  enddo
                  if ( isfrac_found .gt. 0 ) then
                     do ifrac = 1, nfrac
                        if ( ifrac .lt. 100 ) then
                           write(suffix,'(i2.2)') ifrac
                        else
                           write(suffix,'(i3.3)') ifrac
                        endif
                        fracnam = trim(basnam)//trim(suffix)
                        i_new   = i_new + 1
                        new_stochi(i_new) = proc%velostochi(i_stochi)
                        new_stochi(i_new)%ioitem = trim(proc%velostochi(i_stochi)%ioitem(1:i_star))//fracnam
                        call zoek (basnam,1,proc%velostochi(i_stochi)%substance,20,indx)
                        if ( indx .gt. 0 ) then
                           new_stochi(i_new)%substance = fracnam
                        else
                           ! look for linked substance (if not found the substance name stays the same
                           do isfrac2 = 1, nfrac
                              call zoek (sfracs%name(isfrac2),1,proc%velostochi(i_stochi)%substance,20,indx)
                              if ( indx .gt. 0 .and. sfracs%linklist(isfrac,isfrac2) .eq. 1 ) then
                                 new_stochi(i_new)%substance = trim(sfracs%name(isfrac2))//suffix
                                 exit
                              endif
                           enddo
                        endif
                     enddo
                     isfrac_found = 0
                  else
                     i_new = i_new + 1
                     new_stochi(i_new) = proc%velostochi(i_stochi)
                  endif
               else
                  i_new = i_new + 1
                  new_stochi(i_new) = proc%velostochi(i_stochi)
               endif
            enddo
            if (proc%no_velostochi .gt. 0) deallocate(proc%velostochi)
            proc%velostochi    => new_stochi
            proc%no_velostochi = no_velostochi_new

         endif

      enddo

      if (timon) call timstop( ithndl )
      return
 2000 format ( 3a )
      end

max_stars_repo_name: liujiamingustc/phd

max_issues_repo_name: liujiamingustc/phd

max_forks_repo_name: liujiamingustc/phd