text:
      subroutine set_mptr(itmpin,failed_nodes,ith,istp,ithp,
     2     idaughter)
!***********************************************************************
!  Copyright, 2004,  The  Regents  of the  University of California.
!  This program was prepared by the Regents of the University of
!  California at Los Alamos National Laboratory (the University) under
!  contract No. W-7405-ENG-36 with the U.S. Department of Energy (DOE).
!  All rights in the program are reserved by the DOE and the University.
!  Permission is granted to the public to copy and use this software
!  without charge, provided that this Notice and any statement of
!  authorship are reproduced on all copies. Neither the U.S. Government
!  nor the University makes any warranty, express or implied, or
!  assumes any liability or responsibility for the use of this software.
C***********************************************************************
CD1
CD1  PURPOSE
CD1
CD1  This subroutine sets up the starting nodes and times for particle
CD1  tracking, activated with macro mptr.  It is called only once, the
CD1  first time the particle tracker is invoked.
CD1
C***********************************************************************
CD2
CD2  REVISION HISTORY
CD2
CD2 Revision                    ECD
CD2 Date         Programmer     Number  Comments
CD2
CD2 10-JAN-96    S. Henderson   22      Add prolog.
CD2              S. Henderson           Initial implementation.
CD2
CD2 $Log:   /pvcs.config/fehm90/src/set_mptr.f_a  $
!D2
!D2    Rev 2.5   06 Jan 2004 10:43:52   pvcs
!D2 FEHM Version 2.21, STN 10086-2.21-00, Qualified October 2003
!D2
!D2    Rev 2.4   29 Jan 2003 09:15:08   pvcs
!D2 FEHM Version 2.20, STN 10086-2.20-00
!D2
!D2    Rev 2.3   14 Nov 2001 13:04:00   pvcs
!D2 FEHM Version 2.12, STN 10086-2.12-00
!D2 Update the GoldSim / FEHM interface
!D2
!D2    Rev 2.2   06 Jun 2001 13:26:28   pvcs
!D2 FEHM Version 2.11, STN 10086-2.11-00
!D2
!D2    Rev 2.1   30 Nov 2000 12:07:22   pvcs
!D2 FEHM Version 2.10, STN 10086-2.10-00
!D2
!D2    Rev 2.0   Fri May 07 14:45:14 1999   pvcs
!D2 FEHM Version 2.0, SC-194 (Fortran 90)
CD2
CD2 modified to accommodate the using of variable confactors in
CD2 rip-fehm runs.                                  2/6/98.
CD2 Add changes to accommodate the ingrowth model.
CD2                                                 July 19, 1997 cli
CD2
CD2 Made changes for multiple sepcies simulations and transient source
CD2 terms.                Apr. 18, 1997 cli
CD2
CD2    Rev 1.8   Wed Mar 27 11:35:48 1996   hend
CD2 Updated for ptrk output options
CD2
CD2    Rev 1.7   Tue Mar 05 16:20:20 1996   hend
CD2 Added initialization for p array
CD2
CD2    Rev 1.6   Fri Jan 12 17:56:26 1996   llt
CD2 changed mmgetblk agruments
CD2
CD2    Rev 1.5   Thu Jan 11 10:25:56 1996   hend
CD2 Added Prolog
CD2
CD2    Rev 1.4   08/07/95 11:39:40   awolf
CD2 Fixed for DPDP (n0 instead of neq). Also changed for
CD2 breakthrough output
CD2
CD2    Rev 1.3   05/08/95 12:49:10   robinson
CD2 Calls plotc1 to write coordinate info to .trc file
CD2
CD2    Rev 1.2   03/15/95 17:05:34   robinson
CD2 Added diffusion and dispersion to particle tracking model
CD2
CD2    Rev 1.1   02/02/95 15:22:50   llt
CD2 added pvcs log info
CD2
C***********************************************************************
CD3
CD3  REQUIREMENTS TRACEABILITY
CD3
CD3  2.3.5 Cell-based particle-tracking module
CD3
C***********************************************************************
CD4
CD4  SPECIAL COMMENTS AND REFERENCES
CD4
CD4 Requirements from SDN: 10086-RD-2.20-00
CD4   SOFTWARE REQUIREMENTS DOCUMENT (RD) for the
CD4   FEHM Application Version 2.20
CD4
C***********************************************************************
c
c Date: Summer '94
c Author: Stephen Henderson, EES-5
c         Contact: Bruce Robinson, EES-4
c                  George Zyvoloski, EES-5
      use comdti
      use comai
      use combi
      use comci
      use comdi
      use comflow
      use compart
      use davidi
      implicit none

      integer:: i,i2,ithp,num,split,itmp,ith,istp,j,k,lns2
      integer:: iold,iold1,jj,lns,lns1,kin1,itmpin,kin2,kin2_1,rloc
      integer:: idaughter,lstart,nps, failed_nodes !cli_added failed_nodes
      real:: ran_sp
      real*8:: inc,inflow,tmim, denom

      real, parameter:: dayinsecs=86400.

 !cli_added 8/11/1999

c Set the initial node position and starting times for every particle. The
c manner in which these are set is dependent on the type of problem being
c run. Pcnsk must either be all positive or all negative. Positive
c indicates injecting according to source flow, while negative indicates
c evenly distributing the particles over the nodes. Pcnsk gives the
c relative factor to distribute -- i.e. a 2 indicates twice as many
c particles in this node or set of nodes.
c make sure pcnsk is either always positive or always negative

 !set up failed package nodes for radionuclide release.
 !failed_nodes<0,for single node release; otherwise,group of nodes
      iold=ithp
      if(failed_nodes < 0 )then
         kin1=-failed_nodes
         kin2=-failed_nodes
         kin2_1=1
      else
         kin1=insnode(itmpin-1)
         kin2=kin1+failed_nodes
         kin1=kin1+1
         kin2_1=failed_nodes
      endif

 !inject particles into the system

      if(pcnsk(kin1) > 0.) then
   !constant concentration release
         inflow=0
         if((trak_type(ith).eq.1).and.(ico2.lt.0)) then
  !Liquid only and isothermal air-water
            do i=kin1,kin2
               j=ptindex(i)
               if(sk(j).lt.0.) then
                  inflow=inflow+abs(sk(j))*pcnsk(i)
               endif
            enddo
            do k=1,istp
               tmim=t2sk(k)
               itmp=-1
               do i=kin1,kin2
                  j=ptindex(i)
                  if(sk(j).lt.0.)then
                     num=int(dum_p(k)*abs(sk(j))*pcnsk(i)/inflow)
                     if(num.ne.0)then
                        inc=(t2sk(k)-t1sk(k))/num
                        itmp=1
                        if(inc.lt.tmim)tmim=inc
                        do i2=1,num
                           ithp=ithp+1
                           box(ithp,ith)=j
                           frac_done(ithp,ith)=0.
                           theta(ithp,ith)=-99999.
                           start_time(ithp,ith)=
     2                          dayinsecs*(t1sk(k)+inc*(i2-1))
                        enddo
                     endif
                  endif
               enddo
               if(itmp.eq.1)then
	            !converte 2-D arrays into 1-D arrays 1/26/05
                  nsegs(ith)=nsegs(ith)+1
	            lns1=aidex(ith)+nsegs(ith)
                  lns=lns1-1
	            nsport(lns1)=ithp
	            lstart=nsport(lns)+1
	            lsport(lns)=lstart
	            nps=nsport(lns1)-nsport(lns)
                  dum_p(k)=nps
                  if(idaughter.ne.0)then
                     lns2=lns+lns
                     nprevd(lns)=0
                     tmsport(lns2-1)=t1sk(k)*dayinsecs
                     tmsport(lns2)=t2sk(k)*dayinsecs
                     ivdt(lns)=nps/(t2sk(k)-t1sk(k))
                     call sortranl(nps,start_time(lstart,ith),
     2                    box(lstart,ith),tmim,rseed)
                  endif
               endif
            enddo
c     else if vapor only and isothermal air-water
         else if((trak_type(ith).eq.2).and.(ico2.lt.0)) then
            do i=kin1,kin2
               j=ptindex(i)
               if(qh(j).lt.0.) then
                  inflow=inflow+abs(qh(j))*pcnsk(i)
               endif
            enddo
            do k=1,istp
               tmim=t2sk(k)
               itmp=-1
               do i=kin1,kin2
                  j=ptindex(i)
                  if(qh(j).lt.0.)then
                     num=int(dum_p(k)*abs(qh(j))*pcnsk(i)/inflow)
                     if(num.ne.0)then
                        inc=(t2sk(k)-t1sk(k))/num
                        itmp=1
                        if(inc.lt.tmim)tmim=inc
                        do i2=1,num
                           ithp=ithp+1
                           box(ithp,ith)=j
                           frac_done(ithp,ith)=0.
                           theta(ithp,ith)=-99999.
                           start_time(ithp,ith)=dayinsecs*(t1sk(k)+
     2                          inc*(i2-1))
                        enddo
                     endif
                  endif
               enddo
               if(itmp.eq.1)then
	            !converte 2-D arrays into 1-D arrays 1/26/05
                  nsegs(ith)=nsegs(ith)+1
	            lns1=aidex(ith)+nsegs(ith)
                  lns=lns1-1
	            nsport(lns1)=ithp
	            lstart=nsport(lns)+1
	            lsport(lns)=lstart
                  nps=nsport(lns1)-nsport(lns)
	            dum_p(k)=nps
                  if(idaughter.ne.0)then
                     lns2=lns+lns
                     nprevd(lns)=0
                     tmsport(lns2-1)=t1sk(k)*dayinsecs
                     tmsport(lns2)=t2sk(k)*dayinsecs
                     ivdt(lns)=nps/(t2sk(k)-t1sk(k))
                     call sortranl(nps,start_time(lstart,ith),
     2                    box(lstart,ith),tmim,rseed)
                  endif
               endif
            enddo
c     else if liquid only and not isothermal air-water
         else if ((trak_type(ith).eq.1).and.(ico2.ge.0)) then
            do i=kin1,kin2
               j=ptindex(i)
               if((sk(j)*s(j).lt.0.)) then
                  inflow=inflow+abs(sk(j))*pcnsk(i)
               endif
            enddo
            do k=i,istp
               tmim=t2sk(k)
               itmp=-1
               do i=kin1,kin2
                  j=ptindex(i)
                  if(sk(j)*s(j).lt.0)then
                     num=int(dum_p(k)*abs(sk(j))*s(j)*pcnsk(i)/inflow)
                     if(num.ne.0)then
                        inc=(t2sk(k)-t1sk(k))/num
                        itmp=1
                        if(inc.lt.tmim)tmim=inc
                        do i2=1,num
                           ithp=ithp+1
                           box(ithp,ith)=j
                           frac_done(ithp,ith)=0.
                           theta(ithp,ith)=-99999.
                           start_time(ithp,ith)=dayinsecs*(t1sk(k)+
     2                          inc*(i2-1))
                        enddo
                     endif
                  endif
               enddo
               if(itmp.eq.1)then
	            !converte 2-D arrays into 1-D arrays 1/26/05
                  nsegs(ith)=nsegs(ith)+1
	            lns1=aidex(ith)+nsegs(ith)
                  lns=lns1-1
	            nsport(lns1)=ithp
	            lstart=nsport(lns)+1
	            lsport(lns)=lstart
                  nps=nsport(lns1)-nsport(lns)
	            dum_p(k)=nps
                  if(idaughter.ne.0)then
                     lns2=lns+lns
                     nprevd(lns)=0
                     tmsport(lns2-1)=t1sk(k)*dayinsecs
                     tmsport(lns2)=t2sk(k)*dayinsecs
                     ivdt(lns)=nps/(t2sk(k)-t1sk(k))
                     call sortranl(nps,start_time(lstart,ith),
     2                    box(lstart,ith),tmim,rseed)
                  endif
               endif
            enddo
c     else if vapor only and not isothermal air-water
         else if((trak_type(ith).eq.2).and.(ico2.ge.0)) then
            do i=kin1,kin2
               j=ptindex(i)
               if(sk(j)*(1-s(j)).lt.0.) then
                  inflow=inflow+abs(sk(j))*pcnsk(i)
               endif
            enddo
            do k=1,istp
               tmim=t2sk(k)
               itmp=-1
               do i=kin1,kin2
                  j=ptindex(i)
                  if(sk(j)*(1-s(j)).lt.0.)then
                     num=int(dum_p(k)*abs(sk(j))*(1-s(j))*pcnsk(i)/
     2                    inflow)
                     if(num.ne.0)then
                        inc=(t2sk(k)-t1sk(k))/num
                        itmp=1
                        if(inc.lt.tmim)tmim=inc
                        do i2=1,num
                           ithp=ithp+1
                           box(ithp,ith)=j
                           frac_done(ithp,ith)=0.
                           theta(ithp,ith)=-99999.
                           start_time(ithp,ith)=dayinsecs*(t1sk(k)+
     2                          inc*(i2-1))
                        enddo
                     endif
                  endif
               enddo
               if(itmp.eq.1)then
	            !converte 2-D arrays into 1-D arrays 1/26/05
                  nsegs(ith)=nsegs(ith)+1
	            lns1=aidex(ith)+nsegs(ith)
                  lns=lns1-1
	            nsport(lns1)=ithp
	            lstart=nsport(lns)+1
	            lsport(lns)=lstart
                  nps=nsport(lns1)-nsport(lns)
	            dum_p(k)=nps
                  if(idaughter.ne.0)then
                     lns2=lns+lns
                     nprevd(lns)=0
                     tmsport(lns2-1)=t1sk(k)*dayinsecs
                     tmsport(lns2)=t2sk(k)*dayinsecs
                     ivdt(lns)=nps/(t2sk(k)-t1sk(k))
                     call sortranl(nps,start_time(lstart,ith),
     2                    box(lstart,ith),tmim,rseed)
                  endif
               endif
            enddo

         endif

c     Take care of negative pcnsk
      else

         split=0
         do i=kin1,kin2
            split=split-pcnsk(i)
         enddo
         do k=1,istp
            tmim=t2sk(k)
            itmp=-1
            if(kin2_1.gt.dum_p(k))then
               do i=1,dum_p(k)
                  rloc=int(kin2_1*ran_sp(rseed))+kin1
                  ithp=ithp+1

CHari for fracture-matrix stuff
c                  write(6,*)'in set_mptr'
c                  write(6,*)'species','bin','fracflow'
c	            write(6,*)ith, itmpin-1,fracflow(ith,itmpin-1)
                  if(ripfehm.ne.1)then
                     j=ptindex(rloc)
                  else
                     if(ran_sp(rseed).le. fracflow(ith, itmpin-1))then
                        j=ptindex(rloc)
c                       write(6,*)'put particle in fracture',j
                     else
                        j=ptindex(rloc)+neq
c	                write(6,*)'put particle in matrix',j
                    endif
                 endif
                 box(ithp,ith)=j
                 frac_done(ithp,ith)=0.
                 theta(ithp,ith)=-99999.
                 start_time(ithp,ith)=dayinsecs*((t2sk(k)-t1sk(k))*
     2                ran_sp(rseed)+t1sk(k))
              enddo
              itmp=1
           else
              do i=kin1,kin2
CHari for fracture-matrix stuff
                 if(ripfehm.ne.1)then
                    j=ptindex(i)
                 else
                    if(ran_sp(rseed).le. fracflow(ith, itmpin-1))then
                       j=ptindex(i)
c	               write(6,*)'put particle in fracture',j
                    else
                       j=ptindex(i)+neq
c	               write(6,*)'put particle in matrix',j
                    endif
                 endif
                 num=int(dum_p(k)*(-pcnsk(i))/split+0.5)
                 if(num.ne.0)then
                    inc=(t2sk(k)-t1sk(k))/num
                    itmp=1
                    if(inc.lt.tmim)tmim=inc
                    do i2=1,num
                       ithp=ithp+1
                       box(ithp,ith)=j
                       frac_done(ithp,ith)=0.
                       theta(ithp,ith)=-99999.
                       start_time(ithp,ith)=dayinsecs*(t1sk(k)+
     2                      inc*(i2-1))
                    enddo
                 endif
              enddo
           endif
            if(itmp.eq.1)then
	        !converte 2-D arrays into 1-D arrays 1/26/05
              nsegs(ith)=nsegs(ith)+1
	        lns1=aidex(ith)+nsegs(ith)
              lns=lns1-1
	        nsport(lns1)=ithp
	        lstart=nsport(lns)+1
	        lsport(lns)=lstart
              nps=nsport(lns1)-nsport(lns)
	        dum_p(k)=nps
              if(idaughter.ne.0)then
                lns2=lns+lns
                nprevd(lns)=0
                tmsport(lns2-1)=t1sk(k)*dayinsecs
                tmsport(lns2)=t2sk(k)*dayinsecs
                ivdt(lns)=nps/(t2sk(k)-t1sk(k))
                call sortranl(nps,start_time(lstart,ith),
     2               box(lstart,ith),tmim,rseed)
              endif
            endif
         enddo

      endif

      iold1=iold+1
      num_particles(ith)=ithp

      if(pout.eq.0) then
         if(trak_type(ith).eq.1) then
            do i=iold1,ithp
               j=box(i,ith)
               jj=j+npt(ith)
               if (irdof .ne. 13 .or. ifree .ne. 0) then
                  denom = (rolf(j)*ps(j)*s(j)*sx1(j))
               else
                  denom = (rolf(j)*ps(j)*sx1(j))
               end if
               if(denom .ne. 0.) then
                  an(jj)=an(jj)+1/denom
               else
                  an(jj) = 0.
               end if
               anv(j)=an(j)
            enddo
         else
            do i=iold1,ithp
               j=box(i,ith)
               jj=j+npt(ith)
               denom = (rovf(j)*ps(j)*(1-s(j))*sx1(j))
               if(denom .ne. 0.) then
                  an(jj)=an(jj)+1/denom
               else
                  an(jj) = 0.
               end if
               anv(jj)=an(jj)
            enddo
         endif
      endif

      return
      end

max_stars_repo_name: satkarra/FEHM

max_issues_repo_name: satkarra/FEHM

max_forks_repo_name: satkarra/FEHM


text:
SUBROUTINE ZGESVD_F95( A, S, U, VT, WW, JOB, INFO )
!
!  -- LAPACK95 interface driver routine (version 3.0) --
!     UNI-C, Denmark; Univ. of Tennessee, USA; NAG Ltd., UK
!     September, 2000
!
!  .. USE STATEMENTS ..
   USE LA_PRECISION, ONLY: WP => DP
   USE LA_AUXMOD, ONLY: ERINFO, LSAME
   USE F77_LAPACK, ONLY: GESVD_F77 => LA_GESVD
!  .. IMPLICIT STATEMENT ..
   IMPLICIT NONE
!  .. SCALAR ARGUMENTS ..
   CHARACTER(LEN=1), OPTIONAL, INTENT(IN) :: JOB
   INTEGER, INTENT(OUT), OPTIONAL :: INFO
!  .. ARRAY ARGUMENTS ..
   COMPLEX(WP), INTENT(INOUT) :: A(:,:)
   REAL(WP), INTENT(OUT) :: S(:)
   REAL(WP), INTENT(OUT), OPTIONAL :: WW(:)
   COMPLEX(WP), INTENT(OUT), OPTIONAL, TARGET :: U(:,:), VT(:,:)
!----------------------------------------------------------------------
!
!  Purpose
!  =======
!
!       LA_GESVD and LA_GESDD compute the singular values and,
! optionally, the left and/or right singular vectors from the singular
! value decomposition  (SVD) of a real or complex m by n matrix A. The
! SVD of A is written
!                    A = U * SIGMA * V^H
! where SIGMA is an  m by n matrix which is zero except for its
! min(m, n) diagonal elements, U is an m by m orthogonal (unitary)
! matrix, and V is an n by n orthogonal (unitary) matrix. The diagonal
! elements of SIGMA , i.e., the values
!
!      sigma(i)= SIGMA(i,i), i = 1, 2,..., min(m, n)
! are the singular values of A; they are real and non-negative, and are
! returned in descending order. The first min(m, n) columns of U and V
! are the left and right singular vectors of A, respectively.
! LA_GESDD solves the same problem as LA_GESVD but uses a divide and
! conquer method if singular vectors are desired. For large matrices it
! is usually much faster than LA_GESVD when singular vectors are
! desired, but uses more workspace.
!
! Note: The routine returns V^H , not V .
!
! ========
!
!    SUBROUTINE LA_GESVD / LA_GESDD( A, S, U=u, VT=vt, &
!              WW=ww, JOB=job, INFO=info )
!      <type>(<wp>), INTENT(INOUT) :: A(:,:)
!      REAL(<wp>), INTENT(OUT) :: S(:)
!      <type>(<wp>), INTENT(OUT), OPTIONAL :: U(:,:), VT(:,:)
!      REAL(<wp>), INTENT(OUT), OPTIONAL :: WW(:)
!      CHARACTER(LEN=1), INTENT(IN), OPTIONAL :: JOB
!      INTEGER, INTENT(OUT), OPTIONAL :: INFO
!      where
!      <type> ::= REAL | COMPLEX
!      <wp>   ::= KIND(1.0) | KIND(1.0D0)
!
! Arguments
! =========
!
! A      (input/output) REAL or COMPLEX array, shape (:, :) with
!        size(A, 1) = m and size(A, 2) = n.
!        On entry, the matrix A.
!        On exit, if JOB = 'U' and U is not present, then A is
!        overwritten with the first min(m, n) columns of U (the left
!        singular vectors, stored columnwise).
!        If JOB = 'V' and VT is not present, then A is overwritten with
!        the first min(m, n) rows of V^H (the right singular vectors,
!        stored rowwise).
!        In all cases the original contents of A are destroyed.
! S      (output) REAL array, shape (:) with size(S) = min(m, n).
!        The singular values of A, sorted so that S(i) >= S(i+1).
! U      Optional (output) REAL or COMPLEX array, shape (:, :) with
!        size(U, 1) = m  and size(U, 2) = m or min(m, n).
!        If size(U, 2) = m, U contains the m by m matrix U .
!        If size(U; 2) = min(m, n), U contains the first min(m, n)
!        columns of U (the left singular vectors, stored columnwise).
! VT     Optional (output) REAL or COMPLEX array, shape (:, :) with
!        size(VT, 1) = n or min(m, n) and size(VT, 2) = n.
!        If size(VT, 1) = n , VT contains the n by n matrix V^H .
!        If size(VT, 1) = min(m, n), VT contains the first min(m, n)
!        rows of V^H (the right singular vectors, stored rowwise).
! WW     Optional (output) REAL array, shape (:) with size(WW) =
!        min(m, n) - 1
!        If INFO > 0, WW contains the unconverged superdiagonal elements
!        of an upper bidiagonal matrix B whose diagonal is in SIGMA (not
!        necessarily sorted). B has the same singular values as A.
!        Note: WW is a dummy argument for LA_GESDD.
! JOB    Optional (input) CHARACTER(LEN=1).
!        = 'N': neither columns of U nor rows of V^H are returned in
!          array A.
!        = 'U': if U is not present, the first min(m, n) columns of U
!          (the left singular vectors) are returned in array A;
!        = 'V': if VT is not present, the first min(m, n) rows of V^H
!          (the right singular vectors) are returned in array A;
!        Default value: 'N'.
! INFO   Optional (output) INTEGER.
!        = 0: successful exit.
!        < 0: if INFO = -i, the i-th argument had an illegal value.
!        > 0: The algorithm did not converge.
!        If INFO is not present and an error occurs, then the program is
!        terminated with an error message.
!----------------------------------------------------------------------
!  .. LOCAL PARAMETERS ..
   CHARACTER(LEN=8), PARAMETER :: SRNAME = 'LA_GESVD'
!  .. LOCAL SCALARS ..
   CHARACTER(LEN=1) :: LJOB
   CHARACTER(LEN=1) :: LJOBU, LJOBVT
   INTEGER, SAVE :: LWORK = 0
   INTEGER :: N, M, LINFO, LD, ISTAT, ISTAT1, S1U, S2U, S1VT, S2VT, &
              NN, MN, SWW
!  .. LOCAL ARRAYS ..
   COMPLEX(WP), TARGET :: LLU(1,1), LLVT(1,1)
   COMPLEX(WP), POINTER :: WORK(:)
   REAL(WP), POINTER :: RWORK(:)
   COMPLEX(WP) :: WORKMIN(1)
!  .. INTRINSIC FUNCTIONS ..
   INTRINSIC MIN, MAX, PRESENT, SIZE
!  .. EXECUTABLE STATEMENTS ..
   LINFO = 0; ISTAT = 0; M = SIZE(A,1); N = SIZE(A,2)
   LD = MAX(1,M); MN = MIN(M,N)
   IF( PRESENT(JOB) )THEN; LJOB = JOB; ELSE; LJOB = 'N'; ENDIF
   IF( PRESENT(U) )THEN; S1U = SIZE(U,1); S2U = SIZE(U,2)
   ELSE; S1U = 1; S2U = 1; END IF
   IF( PRESENT(VT) )THEN; S1VT = SIZE(VT,1); S2VT = SIZE(VT,2)
   ELSE; S1VT = 1; S2VT = 1; END IF
   IF( PRESENT(WW) )THEN; SWW = SIZE(WW); ELSE; SWW = MN-1; ENDIF
!  .. TEST THE ARGUMENTS
   IF( M < 0 .OR. N < 0 )THEN; LINFO = -1
   ELSE IF( SIZE( S ) /= MN )THEN; LINFO = -2
   ELSE IF( PRESENT(U) .AND. ( S1U /= M .OR. &
            ( S2U /= M .AND. S2U /= MN ) ) )THEN; LINFO = -3
   ELSE IF( PRESENT(VT) .AND. ( ( S1VT /= N .AND. S1VT /= MN ) &
            .OR. S2VT /= N ) )THEN; LINFO = -4
   ELSE IF( SWW /= MN-1 .AND. MN > 0 ) THEN; LINFO = -5
   ELSE IF( PRESENT(JOB) .AND. ( .NOT. ( LSAME(LJOB,'U') .OR. &
            LSAME(LJOB,'V') .OR. LSAME(LJOB,'N') ) .OR. &
            LSAME(LJOB,'U') .AND. PRESENT(U) .OR. &
            LSAME(LJOB,'V') .AND. PRESENT(VT)) )THEN; LINFO = -6
   ELSE
      IF( PRESENT(U) )THEN
         IF( S2U == M )THEN; LJOBU = 'A'; ELSE; LJOBU = 'S'; ENDIF
      ELSE; IF( LSAME(LJOB,'U') ) THEN; LJOBU = 'O'
         ELSE; LJOBU = 'N'; ENDIF
      ENDIF
      IF( PRESENT(VT) )THEN
         IF( S1VT == N )THEN; LJOBVT = 'A'; ELSE; LJOBVT = 'S'; ENDIF
      ELSE; IF( LSAME(LJOB,'V') )THEN; LJOBVT = 'O'
         ELSE; LJOBVT = 'N'; ENDIF
      ENDIF
        ALLOCATE( RWORK( MAX(1, 5*MIN(M,N))), STAT=ISTAT )
      IF( ISTAT == 0 )THEN
        LWORK = -1
	IF( PRESENT(U) ) THEN
	IF ( PRESENT(VT) )THEN
	CALL GESVD_F77( LJOBU, LJOBVT, M, N, A, LD, S, U, MAX(1,S1U), &
&             VT, MAX(1,S1VT), WORKMIN, LWORK, RWORK, LINFO )
           ELSE
             CALL GESVD_F77( LJOBU, LJOBVT, M, N, A, LD, S, U, MAX(1,S1U), &
 &             LLVT, MAX(1,S1VT), WORKMIN, LWORK, RWORK, LINFO )
           ENDIF
         ELSE
           IF ( PRESENT(VT) )THEN
             CALL GESVD_F77( LJOBU, LJOBVT, M, N, A, LD, S, LLU, MAX(1,S1U), &
 &             VT, MAX(1,S1VT), WORKMIN, LWORK, RWORK, LINFO )
           ELSE
             CALL GESVD_F77( LJOBU, LJOBVT, M, N, A, LD, S, LLU, MAX(1,S1U), &
 &             LLVT, MAX(1,S1VT), WORKMIN, LWORK, RWORK, LINFO )
           ENDIF
         ENDIF
         LWORK = WORKMIN(1)
! THE NEXT LINE SHOULD BE REMOVED
         LWORK = LWORK + 1
         ALLOCATE(WORK(LWORK), STAT = ISTAT)
         IF( ISTAT /= 0 )THEN
           DEALLOCATE(WORK,STAT=ISTAT1)
	 NN = 2*MIN(M,N)+MAX(M,N)
         LWORK = MAX( 1, NN); ALLOCATE(WORK(LWORK), STAT=ISTAT)
            IF( ISTAT == 0) CALL ERINFO( -200, SRNAME, LINFO )
         END IF
      END IF
      IF( ISTAT == 0 ) THEN
          IF( PRESENT(U) ) THEN
  	    IF ( PRESENT(VT) )THEN
	      CALL GESVD_F77( LJOBU, LJOBVT, M, N, A, LD, S, U, MAX(1,S1U), &
	      VT, MAX(1,S1VT), WORK, LWORK, RWORK, LINFO )
	    ELSE
	      CALL GESVD_F77( LJOBU, LJOBVT, M, N, A, LD, S, U, MAX(1,S1U), &
&             LLVT, MAX(1,S1VT), WORK, LWORK, RWORK, LINFO )
            ENDIF
	  ELSE
	   IF ( PRESENT(VT) )THEN
	     CALL GESVD_F77( LJOBU, LJOBVT, M, N, A, LD, S, LLU, MAX(1,S1U), &
&             VT, MAX(1,S1VT), WORK, LWORK, RWORK, LINFO )
           ELSE
	     CALL GESVD_F77( LJOBU, LJOBVT, M, N, A, LD, S, LLU, MAX(1,S1U), &
&             LLVT, MAX(1,S1VT), WORK, LWORK, RWORK, LINFO )
           ENDIF
	 ENDIF
         LWORK = INT(WORK(1)+1)
         IF( LINFO > 0 .AND. PRESENT(WW) ) WW(1:MN-1) = RWORK(1:MN-1)
      ELSE; LINFO = -100; ENDIF
      DEALLOCATE(WORK, RWORK, STAT=ISTAT1)
   ENDIF
   CALL ERINFO(LINFO,SRNAME,INFO,ISTAT)
END SUBROUTINE ZGESVD_F95

max_stars_repo_name: IRFM/ALOHA

max_issues_repo_name: IRFM/ALOHA

max_forks_repo_name: IRFM/ALOHA


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

C=======================================================================
      SUBROUTINE SNAP(ndbin, A, IA, isnp, X, Y, Z, NDIM, numnp,
     *  numess, idess, neess, ixeess,
     *  nelb, idelb, numelb, numlnk)
C=======================================================================

      DIMENSION A(*), IA(*)
      REAL X(*), Y(*), Z(*)
      INTEGER NUMESS
      INTEGER IDESS(*)
      INTEGER NEESS(*)
      INTEGER IXEESS(*)
      INTEGER IDELB(NELB)
      INTEGER NUMELB(NELB)
      INTEGER NUMLNK(NELB)

      include 'exodusII.inc'
      INCLUDE 'gp_snap.blk'

C   --   X,Y,Z -- REAL - IN/OUT - Coordinates of nodes
C   --   NDIM   - IN - the number of dimensions
C   --   NUMESS - IN - the number of side sets
C   --   IDESS - OUT - the side set ID for each set
C   --   NEESS - OUT - the number of elements for each set
C   --   IXEESS - OUT - the index of the first element for each set

      call mdrsrv('ISCR',  KISCR,  NUMNP)
      call mdrsrv('VNORM', KVNORM, NUMNP*3)
      call mdrsrv('PLANE', KPLANE, 0)

      indsl = 0
      indma = 0

C ... Find master index
      do 100 i=1, numess
        if (idess(i) .eq. idssma(isnp)) then
          indma = i
          go to 110
        endif
 100  continue
 110  continue

C ... Find slave index
      do 120 i=1, numess
        if (idess(i) .eq. idsssl(isnp)) then
          indsl = i
          go to 130
        endif
 120  continue
 130  continue

C------------------------------------------------------------------------
C ... Get the number of nodes on the slave surface
      call exgsp(ndbin, idess(indsl), nssess, nsdess, ierr)
      if (nsdess .eq. 0) then
C...     Distribution factors not stored, estimate max size of node list
C        based on maximum of 9 nodes/face.
         nsdess = nssess * 9
      end if
C...  Allocate storage for sideset node and count list.
      call mdrsrv('NODSLV', islvnd, nsdess)
      call mdrsrv('NNDSLV', islvnn, nssess)
C...  Get the sideset nodes...
      call exgssn(ndbin, idess(indsl), ia(islvnn), ia(islvnd), ierr)

C ... Get the number of nodes on the master surface
      call exgsp(ndbin, idess(indma), nmsess, nmdess, ierr)
      if (nmdess .eq. 0) then
C...     Distribution factors not stored, estimate max size of node list
C        based on maximum of 9 nodes/face.
         nmdess = nmsess * 9
      end if
C...  Allocate storage for sideset node and count list.
      call mdrsrv('NODMAS', imasnd, nmdess)
      call mdrsrv('NNDMAS', imasnn, nmsess)
C...  Get the sideset nodes...
      call exgssn(ndbin, idess(indma), ia(imasnn), ia(imasnd), ierr)
C------------------------------------------------------------------------
C ... Allocate temporary storage
      nummaf = nmsess
      call mdlong('PLANE', KPLANE, nummaf*4)

      call inirea(3*numnp,  0.0, a(kvnorm))
      call inirea(4*nummaf, 0.0, a(kplane))

C ... Calculate the slave surface normals
      if (ndim .eq. 3) then
        call setnor(.FALSE., numnp, x, y, z, 3, a(kvnorm),
     *    usnorm(isnp), vector(1, isnp),
     *    nssess, nsdess, ia(islvnd))

C ... Calculate the master surface planes
        call setnor(.TRUE., numnp, x, y, z, 4, a(kplane),
     *    0, vector(1, isnp),
     *    nmsess, nmdess, ia(imasnd))

        call snpnod(numnp, ndim, x, y, z,
     *    a(kvnorm), nssess, nsdess, ia(islvnd),
     *    a(kplane), nmsess, nmdess, ia(imasnd),
     *    snptol(isnp), delmax(isnp))
      else
        call setnr2 (.FALSE., numnp, x, y, 3, a(kvnorm),
     *    usnorm(isnp), vector(1, isnp),
     *    nssess, nsdess, ia(islvnd))

        call snpnd2(numnp, ndim, x, y, a(kvnorm),
     *    nmsess, nmdess, ia(imasnd),
     *    snptol(isnp), delmax(isnp))
      end if

      call mddel('NNDMAS')
      call mddel('NODMAS')
      call mddel('NNDSLV')
      call mddel('NODSLV')

      call mddel('PLANE')
      call mddel('VNORM')
      call mddel('ISCR')

      RETURN
      END

ext: f

max_stars_repo_name: jschueller/seacas

max_issues_repo_name: jschueller/seacas

max_forks_repo_name: jschueller/seacas


text:
The Davis area has a few food blogs...here are links.  There was an article featuring several of these blogs in the January 28, 2007 issue of The Davis Enterprise (section C1) if theres an archive out there.

   http://farmophile.wordpress.com. Ranked one of the Top 50 Farm Blogs on the web by Seametrics, farmophile chronicles a Davis familys adventures with local farms, farmers and food.

   http://vanillagarlic.blogspot.com/ Vanilla Garlic by Users/ProtoGarrett.  This Davis blog posts very unique and yummy cupcake recipes as well as restaurant reviews and hilarious essays about the authors cat, EatBeast.

   http://www.somethinginseason.com/ Something in Season by Brendon.  This Davis blog posts glutenfree recipes focusing on seasonal ingredients.

   http://chucrutecomsalsicha.com/ Chucrute com Salsicha  by Fernanda Guimaraes Rosa.  This Davis blog is written in Portuguese so I cant summarize (dont speak Portuguese), beautiful pictures, though!

   http://crumpeted.blogspot.com/ Crumpet Attack! This Davis blog is written by a high school student who uses many ingredients from her own garden.

   http://whatdidyoueat.typepad.com/what_did_you_eat/ What Did You Eat? by Sher.  This Davis blog discusses...well, what did you eat?

   http://theshamelessfoodies.blogspot.com/ The Shameless Foodies by Users/KristenMartin Kristen and Users/BrooksHayes Brooks.  This Davis blog posts restaurant reviews and other culinary adventures.

   http://sacatomato.com/ Sacatomato by Ann Martin Rolke and Jennifer Cliff.  This blog posts all food things Sacramento, lots of events listed.

   http://elise.com/recipes/ Simply Recipes by Elise Bauer.  This Sacramentoarea food blog was recently voted Best Food Blog Overall in the 2006 Food Blog Awards.

   http://pearsandplums.com/recipes/ Pears and Plums by Laura. This blog by a local Davis resident posts about traditional foods (organic foods, grassfed beef, raw milk, pasturedpoultry, and more )

  http://thequarteracrefarm.com/ The QuarterAcre Farm Growing Food and Eating Well  How I Kept the Patio, Lost the Lawn, and Fed My Family for A year

  http://rollingsreliable.com/ RollingsReliable by Users/OliviaY  A baking/cooking blog by local Daviswikian who tries out both her own recipes and recipes from a vast collection of cookbooks.

  http://toxicfoodie.org/ Toxic Foodie by a UC Davis food toxicology undergrad! A food blog with a twist! Giving you the 411 on whats in the food you eat!

If anyone can pretty up this list, please do.  Please add more links to the list as you discover them!

There is even a Local Food Week.

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
!
! LBLRTM_r7p1_Module
!
! Module containing procedures for LBLRTM Record r7p1.
!
!
! CREATION HISTORY:
!       Written by:   Paul van Delst, 24-Dec-2012
!                     paul.vandelst@noaa.gov
!

MODULE LBLRTM_r7p1_Module

  ! -----------------
  ! Environment setup
  ! -----------------
  ! Module usage
  USE Type_Kinds     , ONLY: fp
  USE File_Utility   , ONLY: File_Open
  USE Message_Handler, ONLY: SUCCESS, FAILURE, INFORMATION, Display_Message
  ! Line-by-line model parameters
  USE LBL_Parameters
  ! Disable implicit typing
  IMPLICIT NONE

  ! ----------
  ! Visibility
  ! ----------
  ! Everything private by default
  PRIVATE
  ! Datatypes
  PUBLIC :: LBLRTM_r7p1_type
  ! Procedures
  PUBLIC :: LBLRTM_r7p1_Write

  ! -----------------
  ! Module parameters
  ! -----------------
  CHARACTER(*), PARAMETER :: MODULE_VERSION_ID = &
  ! Message string length
  INTEGER, PARAMETER :: ML = 256
  ! The record I/O format
  CHARACTER(*), PARAMETER :: LBLRTM_R7P1_FMT = '(f10.3,f10.4,3i5,10x,a35)'

  ! -------------
  ! Derived types
  ! -------------
  TYPE :: LBLRTM_r7p1_type
    REAL(fp)      :: v1f    = 0.0_fp  !  Wavenumber of initial filter value
    REAL(fp)      :: dvf    = 0.0_fp  !  Wavenumber increment between filter values
    INTEGER       :: npts   = 0       !  Number of filter values
    INTEGER       :: jemit  = 0       !  Integration flag. 0==transmittance; 1==radiance
    INTEGER       :: nnfile = 0       !  Unit number for filtered sequential output
    CHARACTER(35) :: heddr  = ''      !  User identification
  END TYPE LBLRTM_r7p1_type

CONTAINS

  FUNCTION LBLRTM_r7p1_Write(r7p1,fid) RESULT(err_stat)

    ! Arguments
    TYPE(LBLRTM_r7p1_type), INTENT(IN) :: r7p1
    INTEGER               , INTENT(IN) :: fid
    ! Function result
    INTEGER :: err_stat
    ! Function parameters
    CHARACTER(*), PARAMETER :: ROUTINE_NAME = 'LBLRTM_r7p1_Write'
    ! Function variables
    CHARACTER(ML) :: msg
    CHARACTER(ML) :: io_msg
    INTEGER :: io_stat

    ! Setup
    err_stat = SUCCESS
    ! ...Check unit is open
    IF ( .NOT. File_Open(fid) ) THEN
      msg = 'File unit is not connected'
      CALL Cleanup(); RETURN
    END IF

    ! Write the record
    WRITE( fid,FMT=LBLRTM_R7P1_FMT,IOSTAT=io_stat,IOMSG=io_msg) r7p1
    IF ( io_stat /= 0 ) THEN
      msg = 'Error writing record - '//TRIM(io_msg)
      CALL Cleanup(); RETURN
    END IF

  CONTAINS

    SUBROUTINE CleanUp()
      err_stat = FAILURE
      CALL Display_Message( ROUTINE_NAME,msg,err_stat )
    END SUBROUTINE CleanUp

  END FUNCTION LBLRTM_r7p1_Write

END MODULE LBLRTM_r7p1_Module

max_stars_repo_name: hsbadr/crtm

max_issues_repo_name: hsbadr/crtm

max_forks_repo_name: hsbadr/crtm


text:
      program prog
      write(10, end=20) num
   20 continue
      end


text:
#include "rundeck_opts.h"

      MODULE ODIAG
!@sum  ODIAG ocean diagnostic arrays (incl. dynamic sea ice)
!@auth Gary Russell/Gavin Schmidt
!@ver  1.0
#ifdef TRACERS_OCEAN
      USE TRACER_COM, only : ntm
#endif
      USE OCEAN, only : im,jm,lmo
      USE STRAITS, only : nmst
      USE DAGCOM, only : npts  ! needed for conservation diags
      IMPLICIT NONE
      SAVE
      INTEGER, PARAMETER :: KOIJ=6,KOIJL=22,KOL=6,KOLNST=8,
     *     KACCO=IM*JM*KOIJ + IM*JM*LMO*KOIJL + LMO*KOL + LMO*NMST
     *     *KOLNST
!@var OIJ   lat-lon ocean diagnostics (on ocean grid)
!@var OIJL  3-dimensional ocean diagnostics
!@var OL    vertical ocean diagnostics
!@var OLNST strait diagnostics
      REAL*8, DIMENSION(IM,JM,KOIJ)  :: OIJ
      REAL*8, DIMENSION(IM,JM,LMO,KOIJL) :: OIJL
      REAL*8, DIMENSION(LMO,KOL)   :: OL
      REAL*8, DIMENSION(LMO,NMST,KOLNST):: OLNST
!@var IJ_xxx Names for OIJ diagnostics
      INTEGER IJ_HBL,IJ_BO,IJ_BOSOL,IJ_USTAR,IJ_SSH,IJ_PB
!@var lname_oij Long names for OIJ diagnostics
      CHARACTER*50, DIMENSION(KOIJ) :: LNAME_OIJ
!@var sname_oij Short names for OIJ diagnostics
      CHARACTER*30, DIMENSION(KOIJ) :: SNAME_OIJ
!@var units_oij Units for OIJ diagnostics
      CHARACTER*50, DIMENSION(KOIJ) :: UNITS_OIJ
!@var ia_oij IDACC numbers for OIJ diagnostics
      INTEGER, DIMENSION(KOIJ) :: IA_OIJ
!@var scale_oij scales for OIJ diagnostics
      REAL*8, DIMENSION(KOIJ) :: SCALE_OIJ
!@var ijgrid_oij Grid descriptor for OIJ diagnostics
       INTEGER, DIMENSION(KOIJ) :: IJGRID_OIJ

!@var IJL_xxx Names for OIJL diagnostics
      INTEGER IJL_MO,IJL_G0M,IJL_S0M,IJL_GFLX,IJL_SFLX,IJL_MFU,IJL_MFV
     *     ,IJL_MFW,IJL_GGMFL,IJL_SGMFL,IJL_KVM,IJL_KVG,IJL_WGFL
     *     ,IJL_WSFL
!@var LN_xxx Names for OLNST diagnostics
      INTEGER LN_KVM,LN_KVG,LN_WGFL,LN_WSFL,LN_MFLX,LN_GFLX,LN_SFLX
     *     ,LN_ICFL
!@var L_xxx Names for OL diagnostics
      INTEGER L_RHO,L_TEMP,L_SALT
!@var icon_xx indexes for conservation quantities
      INTEGER icon_OCE,icon_OKE,icon_OAM,icon_OMS,icon_OSL
!@var kbasin integer index of which basin a particular ocean point is in
      INTEGER, DIMENSION(IM,JM) :: KBASIN
!@var XLB label for diagnostic titles
      CHARACTER XLB*30
!@var FLAT latitude values on primary and secondary ocean grids
      REAL*8, DIMENSION(JM,2) :: FLAT
!@var FLON longitude values on primary and secondary ocean grids
      REAL*8, DIMENSION(IM,2) :: FLON
!@var iu_otj unit number for ascii output of ocean transports
      INTEGER iu_otj
!@var BASIN names of ocean basins for diag output
      CHARACTER*16, DIMENSION(4) :: BASIN=
     *     (/"Atlantic","Pacific ","Indian  ","Global  "/)
C****

#ifdef TRACERS_OCEAN
!@var KTOIJL number of 3-dimensional ocean tracer diagnostics
      INTEGER, PARAMETER :: KTOIJL=10
!@var TOIJL  3-dimensional ocean tracer diagnostics
      REAL*8, DIMENSION(IM,JM,LMO,KTOIJL,NTM)  :: TOIJL
!@var toijl_xxx indices for TOIJL diags
      INTEGER, PARAMETER :: toijl_conc=1,toijl_tflx=2,toijl_gmfl=6
!@var TLNST strait diagnostics
      REAL*8, DIMENSION(LMO,NMST,KOLNST,NTM):: TLNST
#endif
      END MODULE ODIAG

      SUBROUTINE io_ocdiag(kunit,it,iaction,ioerr)
!@sum  io_ocdiag reads and writes ocean diagnostic arrays to file
!@auth Gavin Schmidt
!@ver  1.0
      USE MODEL_COM, only : ioread,iowrite,iowrite_mon,iowrite_single
     *     ,irsfic,irerun,irsficno,ioread_single,lhead
      USE ODIAG
      IMPLICIT NONE

      INTEGER kunit   !@var kunit unit number of read/write
      INTEGER iaction !@var iaction flag for reading or writing to file
!@var IOERR 1 (or -1) if there is (or is not) an error in i/o
      INTEGER, INTENT(INOUT) :: IOERR
!@var HEADER Character string label for individual records
      CHARACTER*80 :: HEADER, MODULE_HEADER = "OCDIAG01"
!@var it input/ouput value of hour
      INTEGER, INTENT(INOUT) :: it
!@var OIJ4,OIJL4,OLNST4,OL4 dummy arrays for reading diag. files
      REAL*4, DIMENSION(IM,JM,KOIJ)  :: OIJ4
      REAL*4, DIMENSION(IM,JM,LMO,KOIJL) :: OIJL4
      REAL*4, DIMENSION(LMO,KOL)   :: OL4
      REAL*4, DIMENSION(LMO,NMST,KOLNST):: OLNST4
#ifdef TRACERS_OCEAN
      REAL*4, DIMENSION(IM,JM,LMO,KTOIJL,NTM)  :: TOIJL4
      REAL*4, DIMENSION(LMO,NMST,KOLNST,NTM):: TLNST4
!@var TR_HEADER Character string label for individual tracer records
      CHARACTER*80 :: TR_HEADER, TR_MODULE_HEADER = "TROCDIAG01"

      write(TR_MODULE_HEADER(lhead+1:80),'(a19,i2,a1,i2,a9,i4,a4)')
     *     'R8 Toijl(im,jm,lmo,',ktoijl,',',ntm,'), Tlnst(',
     *     LMO*NMST*KOLNST*NTM,'),it'
#endif
      write(MODULE_HEADER(lhead+1:80),'(a13,i2,a13,i2,a1,  i2,a5,i2,
     *  a1,i2,a8,i4,a)') 'R8 Oij(im,jm,',koij,'),Oijl(im,jm,',lmo,',',
     *  koijl,'),Ol(',lmo,   ',',kol,'),OLNST(',LMO*NMST*KOLNST,'),it'

      SELECT CASE (IACTION)
      CASE (IOWRITE,IOWRITE_MON)  ! output to standard restart file
        WRITE (kunit,err=10) MODULE_HEADER,OIJ,OIJL,OL,OLNST,it
#ifdef TRACERS_OCEAN
        WRITE (kunit,err=10) TR_MODULE_HEADER,TOIJL,TLNST,it
#endif
      CASE (IOWRITE_SINGLE)    ! output to acc file
        MODULE_HEADER(LHEAD+1:LHEAD+2) = 'R4'
        WRITE (kunit,err=10) MODULE_HEADER,REAL(OIJ,KIND=4),
     *        REAL(OIJL,KIND=4),REAL(OL,KIND=4),REAL(OLNST,KIND=4),it
#ifdef TRACERS_OCEAN
        TR_MODULE_HEADER(LHEAD+1:LHEAD+2) = 'R4'
        WRITE (kunit,err=10) TR_MODULE_HEADER,REAL(TOIJL,KIND=4),
     *       REAL(TLNST,KIND=4),it
#endif
      CASE (IOREAD:)            ! input from restart file
        SELECT CASE (IACTION)
        CASE (IRSFICNO)         ! initial conditions
        CASE (ioread_single)    ! accumulate diagnostic files
          READ (kunit,err=10) HEADER,OIJ4,OIJL4,OL4,OLNST4,it
C**** accumulate diagnostics
          OIJ=OIJ+OIJ4
          OIJL=OIJL+OIJL4
          OL=OL+OL4
          OLNST=OLNST+OLNST4
          IF (HEADER(1:LHEAD).NE.MODULE_HEADER(1:LHEAD)) THEN
            PRINT*,"Discrepancy in module version ",HEADER
     *           ,MODULE_HEADER
            GO TO 10
          END IF
#ifdef TRACERS_OCEAN
          READ (kunit,err=10) TR_HEADER,TOIJL4,TLNST4,it
C**** accumulate diagnostics
          TOIJL=TOIJL+TOIJL4
          TLNST=TLNST+TLNST4
          IF (TR_HEADER(1:LHEAD).NE.TR_MODULE_HEADER(1:LHEAD)) THEN
            PRINT*,"Discrepancy in module version ",TR_HEADER
     *           ,TR_MODULE_HEADER
            GO TO 10
          END IF
#endif
        CASE (ioread,irerun)    ! restarts
          READ (kunit,err=10) HEADER,OIJ,OIJL,OL,OLNST,it
          IF (HEADER(1:LHEAD).NE.MODULE_HEADER(1:LHEAD)) THEN
            PRINT*,"Discrepancy in module version ",HEADER
     *           ,MODULE_HEADER
            GO TO 10
          END IF
#ifdef TRACERS_OCEAN
          READ (kunit,err=10) TR_HEADER,TOIJL,TLNST,it
          IF (TR_HEADER(1:LHEAD).NE.TR_MODULE_HEADER(1:LHEAD)) THEN
            PRINT*,"Discrepancy in module version ",TR_HEADER
     *           ,TR_MODULE_HEADER
            GO TO 10
          END IF
#endif
        END SELECT
      END SELECT

      RETURN
 10   IOERR=1
      RETURN
C****
      END SUBROUTINE io_ocdiag

      SUBROUTINE DIAGCO (M)
!@sum  DIAGCO Keeps track of the ocean conservation properties
!@auth Gary Russell/Gavin Schmidt
!@ver  1.0
      USE ODIAG, only : icon_OCE,icon_OKE,icon_OMS,icon_OSL,icon_OAM
      IMPLICIT NONE
!@var M index denoting from where DIAGCO is called (see DIAGCA)
      INTEGER, INTENT(IN) :: M
      REAL*8, EXTERNAL :: conserv_OCE,conserv_OKE,conserv_OMS
     *     ,conserv_OSL,conserv_OAM

C**** OCEAN MASS
      CALL conserv_DIAG(M,conserv_OMS,icon_OMS)

C**** OCEAN ANGULAR MOMENTUM
      CALL conserv_DIAG(M,conserv_OAM,icon_OAM)

C**** OCEAN KINETIC ENERGY
      CALL conserv_DIAG(M,conserv_OKE,icon_OKE)

C**** OCEAN POTENTIAL ENTHALPY
      CALL conserv_DIAG(M,conserv_OCE,icon_OCE)

C**** OCEAN SALT
      CALL conserv_DIAG(M,conserv_OSL,icon_OSL)

C****
      RETURN
      END SUBROUTINE DIAGCO

      SUBROUTINE init_ODIAG
!@sum  init_ODIAG initialises ocean diagnostics
!@auth Gavin Schmidt
!@ver  1.0
      USE CONSTANT, only : bygrav
      USE MODEL_COM, only : dtsrc
      USE OCEAN, only : ze
      USE DAGCOM, only : ia_src,conpt0,zoc,zoc1
      USE ODIAG
      IMPLICIT NONE
      LOGICAL :: QCON(NPTS), T = .TRUE. , F = .FALSE.
      CHARACTER CONPT(NPTS)*10
      INTEGER k

C**** Set names for OIJL diagnostics
      IJL_MO=1  ; IJL_G0M=2  ; IJL_S0M=3   ; IJL_MFU=4    ; IJL_MFV=5
      IJL_MFW=6 ; IJL_KVM=13; IJL_KVG=14 ; IJL_WGFL=15 ; IJL_WSFL=16
C**** These flux diagnostics need 3 spots each
      IJL_GFLX=7 ; IJL_SFLX=10 ; IJL_GGMFL=17 ; IJL_SGMFL=20
C**** Set names for OLNST diagnostics
      LN_KVM=1  ; LN_KVG=2  ; LN_WGFL=3 ; LN_WSFL=4
      LN_MFLX=5 ; LN_GFLX=6 ; LN_SFLX=7 ; LN_ICFL=8
C**** Set names for OL diagnostics
      L_RHO=1   ; L_TEMP=2  ; L_SALT=3

C**** set properties for OIJ diagnostics
      k=0
c
      k=k+1
      IJ_HBL=k
      lname_oij(k)="Ocean Boundary layer depth (KPP) x PO4"
      sname_oij(k)="oij_hbl"
      units_oij(k)="m"
      ia_oij(k)=ia_src
      scale_oij(k)=0.25
      ijgrid_oij(k)=1

      k=k+1
      IJ_BO=k
      lname_oij(k)="Surface buoyancy forcing (KPP) x PO4"
      sname_oij(k)="oij_bo"
      units_oij(k)="10^-7 m^2/s^3"
      ia_oij(k)=ia_src
      scale_oij(k)=0.25*1d7
      ijgrid_oij(k)=1

      k=k+1
      IJ_BOSOL=k
      lname_oij(k)="Surface solar buoyancy flux x PO4"
      sname_oij(k)="oij_bosol"
      units_oij(k)="10^-7 m^2/s^3"
      ia_oij(k)=ia_src
      scale_oij(k)=0.25*1d7
      ijgrid_oij(k)=1

      k=k+1
      IJ_USTAR=k
      lname_oij(k)="Surface friction speed x PO4"
      sname_oij(k)="oij_ustar"
      units_oij(k)="m/s"
      ia_oij(k)=ia_src
      scale_oij(k)=0.25
      ijgrid_oij(k)=1

      k=k+1
      IJ_SSH=k
      lname_oij(k)="Ocean surface height"
      sname_oij(k)="oij_ssh"
      units_oij(k)="m"
      ia_oij(k)=ia_src
      scale_oij(k)=bygrav
      ijgrid_oij(k)=1

      k=k+1
      IJ_PB=k
      lname_oij(k)="Ocean bottom pressure anomaly"
      sname_oij(k)="oij_pb"
      units_oij(k)="Pa"
      ia_oij(k)=ia_src
      scale_oij(k)=1.
      ijgrid_oij(k)=1

      if (k.gt.KOIJ) then
        write(6,*) "Too many OIJ diagnostics: increase KOIJ to at least"
     *       ,k
        call stop_model("OIJ diagnostic error",255)
      end if

C**** Set up oceanic component conservation diagnostics
C**** Oceanic mass
      CONPT=CONPT0
      CONPT(8)="OCN PHYS"
      QCON=(/ F, F, F, T, F, F, F, T, T, T, T/)
      CALL SET_CON(QCON,CONPT,"OCN MASS","(10**2 KG/M^2) ",
     *     "(10**-8 KG/SM^2)",1d-2,1d8,icon_OMS)
C**** Oceanic angular momentum
      QCON=(/ F, F, F, T, F, F, F, T, T, T, T/)
      CALL SET_CON(QCON,CONPT,"OCN AM  ","(10**12 JS/M^2)",
     *     "(10**2 J/M^2)   ",1d-12,1d-2,icon_OAM)
C**** Oceanic kinetic energy
      QCON=(/ F, F, F, T, F, F, F, T, T, T, T/)
      CALL SET_CON(QCON,CONPT,"OCEAN KE","(J/M^2)        ",
     *     "(10**-6 W/M^2)  ",1d0,1d6,icon_OKE)
C**** Oceanic potential enthalpy (heat)
      QCON=(/ F, F, F, T, F, F, F, T, T, T, T/)
      CALL SET_CON(QCON,CONPT,"OCN HEAT","(10**6 J/M^2)  ",
     *     "(10**-2 W/M^2)  ",1d-6,1d2,icon_OCE)
C**** Oceanic salt mass
      QCON=(/ F, F, F, T, F, F, F, T, T, T, T/)
      CALL SET_CON(QCON,CONPT,"OCN SALT","(10 KG/M^2)    ",
     *     "(10**-9 KG/SM^2)",1d-1,1d9,icon_OSL)
C**** Initialise ocean basins
      CALL OBASIN

C**** Define ocean depths for diagnostic output
      ZOC1(1:LMO+1) = ZE(0:LMO)
      ZOC(1:LMO) = 0.5*(ZE(1:LMO)+ZE(0:LMO-1))

C****
      RETURN
      END SUBROUTINE init_ODIAG

      SUBROUTINE reset_odiag(isum)
!@sum  reset_odiag zeros out ocean diagnostics if needed
!@auth G. Schmidt
!@ver  1.0
      USE ODIAG, only : oij,oijl,ol,olnst
#ifdef TRACERS_OCEAN
     *     ,toijl,tlnst
#endif
      IMPLICIT NONE
      INTEGER, INTENT(IN) :: isum  ! needed for plug-play compatibility

      OIJ=0. ; OIJL=0. ; OL=0. ; OLNST=0.

#ifdef TRACERS_OCEAN
      TOIJL=0. ; TLNST = 0.
#endif

      return
      END SUBROUTINE reset_odiag


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

         function foounfr(i)
           integer foounfr,i,j
           open(12,file='unffile1.dat', form='unformatted')
           read(12) j
!            print *, "i, j:", i,j
           foounfr=i
         end

         subroutine unfread
         integer a,b,foounfr
         integer c(3)
         character*20 writetome
         c(1) = 1
         c(2) = 2
         c(3) = 3
         a=3
         b=4
         open(10,file='unffile2.dat',form='unformatted')
         read(10) a,c(foounfr(2)),b
!          print *, "a c b", a, c, b
         if  (a .eq. 1 .and. c(1) .eq. 1 .and. c(2) .eq. 1 .and. &
             c(3) .eq. 3 .and. b .eq. 3) then
             print *, "PASS"
         else
             print *, "FAIL"
         endif
         end

         function foow2(i)
         integer i,foow2
         write(10) i,i+2,i+3
!            print *, "write to 10", i, i+2, i+3
         foow2 = i+1
         end

         function foow(i)
           integer foow,i,foow2
           open(12,file='unffile1.dat', form='unformatted')
           write(unit=12) i+7, foow2(1), i
!            print *, "write to 12", i+7, i+1, i
           foow=i+3
           close(12)
         end

         program mm
         integer a,b,foow
         character*20 writetome
         a=1
         b=4
         open(10,file='unffile2.dat', form='unformatted')
!          print *, "write to 10", a
         write(10) a,foow(1),b
!          print *, "write to 10", b
         close(10)

         call unfread()
         end

max_stars_repo_name: kammerdienerb/flang

max_issues_repo_name: kammerdienerb/flang

max_forks_repo_name: kammerdienerb/flang


text:
!-----------------------------------------------------------------------
!     Module:        fieldlines_init_errorfield
!     Authors:       S. Lazerson (samuel.lazerson@ipp.mpg.de)
!     Date:          05/14/2021
!     Description:   This subroutine adds and errorfield to the existing
!                    background field.  We assume a perturbation to the
!                    radial component of the magnetic field
!                       B_R=B_0*cos(n*phi)
!                    then from div(B)=0 and B_Z=0 we get
!                       B_phi=-B_0*sin(n*phi)/n
!-----------------------------------------------------------------------
      SUBROUTINE fieldlines_init_errorfield
!-----------------------------------------------------------------------
!     Libraries
!-----------------------------------------------------------------------
      USE stel_kinds, ONLY: rprec
      USE fieldlines_runtime
      USE fieldlines_grid, ONLY: raxis,phiaxis,zaxis, nr, nphi, nz, &
                                 rmin, rmax, zmin, zmax, phimin, &
                                 phimax, B_R, B_Z, B_PHI
      USE mpi_params
      USE mpi_inc
!-----------------------------------------------------------------------
!     Local Variables
!          ier            Error Flag
!          iunit          File ID Number
!-----------------------------------------------------------------------
      IMPLICIT NONE
      INTEGER, PARAMETER :: BYTE_8 = SELECTED_INT_KIND (8)
      INTEGER(KIND=BYTE_8),ALLOCATABLE :: mnum(:), moffsets(:)
      INTEGER :: numprocs_local, mylocalid, mylocalmaster
      INTEGER :: MPI_COMM_LOCAL
      INTEGER(KIND=BYTE_8) :: chunk
      INTEGER :: ier, iunit, s, i, j, mystart, myend, k, ig
!-----------------------------------------------------------------------
!     Begin Subroutine
!-----------------------------------------------------------------------

      ! Divide up Work
      numprocs_local = 1; mylocalid = master
#if defined(MPI_OPT)
      CALL MPI_COMM_DUP( MPI_COMM_SHARMEM, MPI_COMM_LOCAL, ierr_mpi)
      CALL MPI_COMM_RANK( MPI_COMM_LOCAL, mylocalid, ierr_mpi )              ! MPI
      CALL MPI_COMM_SIZE( MPI_COMM_LOCAL, numprocs_local, ierr_mpi )          ! MPI
#endif
      mylocalmaster = master

      IF (lverb) THEN
         WRITE(6,'(A)')   '----- ERROR FIELD Information -----'
         DO ig = 1, 20
            IF (errorfield_amp(ig) .eq. 0) CYCLE
            WRITE(6,'(A,I2.2,A,E11.4,A,F5.3,A)') 'n = ',ig,'; B = ',errorfield_amp(ig),&
                          ' [T]; phase = ',errorfield_phase(ig), ' [rad]'
         END DO
         CALL FLUSH(6)
      END IF

      ! Break up the Work
      CALL MPI_CALC_MYRANGE(MPI_COMM_LOCAL,1, nr*nphi*nz, mystart, myend)

      ! Loop
      ! Here's the logic the index defines n
      DO s = mystart, myend
         i = MOD(s-1,nr)+1
         j = MOD(s-1,nr*nphi)
         j = FLOOR(REAL(j) / REAL(nr))+1
         k = CEILING(REAL(s) / REAL(nr*nphi))
         DO ig = 1, 20
            B_R(i,j,k) = B_R(i,j,k) + errorfield_amp(ig)*COS(ig*phiaxis(j)-errorfield_phase(ig))
            B_PHI(i,j,k) = B_PHI(i,j,k) - errorfield_amp(ig)*SIN(ig*phiaxis(j)-errorfield_phase(ig))/ig
         END DO
      END DO

#if defined(MPI_OPT)
      CALL MPI_BARRIER(MPI_COMM_LOCAL,ierr_mpi)
      IF (ierr_mpi /=0) CALL handle_err(MPI_BARRIER_ERR,'fieldlines_init_coil1',ierr_mpi)
      CALL MPI_COMM_FREE(MPI_COMM_LOCAL,ierr_mpi)
      IF (ierr_mpi /= MPI_SUCCESS) CALL handle_err(MPI_ERR,'fieldlines_init_coil: MPI_COMM_LOCAL',ierr_mpi)
      CALL MPI_BARRIER(MPI_COMM_FIELDLINES,ierr_mpi)
      IF (ierr_mpi /=0) CALL handle_err(MPI_BARRIER_ERR,'fieldlines_init_coil',ierr_mpi)
#endif

      RETURN
!-----------------------------------------------------------------------
!     End Subroutine
!-----------------------------------------------------------------------
      END SUBROUTINE fieldlines_init_errorfield

max_stars_repo_name: joseluisvelasco/STELLOPT

max_issues_repo_name: joseluisvelasco/STELLOPT

max_forks_repo_name: joseluisvelasco/STELLOPT


text:
 Ms. Day can use this space to introduce herself to the community.

20120208 16:09:14 nbsp A belated Welcome to the Wiki, HLDAY! Are you related professionally with the Sterling Pointe Apartments? Users/PeteB

20120208 16:36:24 nbsp Hi Users/StevenDaubert

20120208 19:21:30 nbsp Let me add my welcome.  You might also find it useful to read Welcome to the wiki/business owner.  It will explain why the comments were restored to the Sterling Point Apts page. Users/CovertProfessor

20120213 09:40:22 nbsp Please dont delete user comments. Users/TomGarberson

20120215 17:33:25 nbsp Please stop deleting users comments or you will be put up for a ban vote. This is the 3rd time youve done this. Users/PeteB

20120223 14:11:08 nbsp I might be beating a dead horse here... But please do not delete user comments.  Users/LoriOrf

20120224 11:11:11 nbsp HLDAYThere is a standing rule that we dont delete other peoples comments. We have repeatedly been trying to talk to you about this issue on this page Users/PeteB

20120224 11:13:05 nbsp Please, please, please read the Welcome to the Wiki/Business Owner page.

Regarding the comments, one of the ones youre deleting was posted just 15 days ago:
 20120209 12:28:04 nbsp This place has changed so much for the worse. Not only has there been landscaping work going on for months (including on Thanksgiving day), but it will continue to go on for months to come. The workers are extremely rude to the residents and it has become an inconvenience to even live here. The management does not give out notices to all residents for things like parking lot work which means that your car gets towed (at your own expense) for parking at your own apartment. I used to love this place but I do not recommend signing a lease here because the new owners are more concerned about making money than they are about the people who live here. Users/LMT

When did management change? Users/TomGarberson

20120224 11:29:19 nbsp Its not your page. Sorry. Users/MeggoWaffle

20120224 12:19:23 nbsp Im sorry you feel youve been treated unfairly.  Im not sure how to resolve your concerns without removing community input entirely, though.  Of course comments from users are going to reflect their opinions; thats exactly what they are.  The wiki is used to gather information about everything in Davis, including both raw fact and peoples viewpoints.  Thats part of why its so great.  You can learn about restaurants, apartments, events, streets, activities, and just about anything else in Davis from many different peoples perspectives.  You can also learn about customer service from the way business owners/representatives respond to negative feedback.  That intro page has some examples of people responding extremely well.

Comments from before a change in management are often archived  removed from the page and put onto a separate comment archive page.  The change of management is absolutely good information, and it should be reflected on the page.  Hence my question: when did management change? Users/TomGarberson

20121213 12:03:27 nbsp HLDAYplease stop trying to make the inappropriate changes on Brush Creek Apartments again. Please tell us what concerns you have and we can work together to repair what needs to be updated Users/PeteB

20121213 12:36:29 nbsp Hilary email me at my last @ gmail Users/StevenDaubert

20130129 16:11:40 nbsp HLDAYWe have repeatedly asked you not to delete other peoples comments from Brush Creek Apartments. We have offered to work with you with absolutely no response from you. If you continue to do this, I will request that the page be locked so only admins can make changed to it. Users/PeteB

20130311 15:05:49 nbsp Please dont remove relevant historical information from Oakshade Commons Apartments. Thank you Users/PeteB


text:
      REAL             FUNCTION SBESJ1 (X)
c Copyright (c) 1996 California Institute of Technology, Pasadena, CA.
c ALL RIGHTS RESERVED.
c Based on Government Sponsored Research NAS7-03001.
c>> 1996-03-30 SBESJ1 Krogh  Added external statement.
C>> 1995-11-13 SBESJ1 Krogh  Changed data statment for C converstion.
C>> 1995-11-03 SBESJ1 Krogh  Removed blanks in numbers for C conversion.
C>> 1994-11-11 SBESJ1 Krogh   Declared all vars.
C>> 1994-10-20 SBESJ1 Krogh  Changes to use M77CON
C>> 1992-02-07 SBESJ1 WV Snyder Correct sign for X < 0.
C>> 1990-11-29 SBESJ1 CLL
C>> 1985-08-02 SBESJ1 Lawson  Initial code.
C JUNE 1978 EDITION.  W. FULLERTON, C3, LOS ALAMOS SCIENTIFIC LAB.
C C.L.LAWSON & S.CHAN, JPL, 1984 FEB ADAPTED TO JPL MATH77 LIBRARY.
c     ------------------------------------------------------------------
c--S replaces "?": ?BESJ1, ?BMP1, ?INITS, ?CSEVL
c     ------------------------------------------------------------------
      EXTERNAL R1MACH, SCSEVL
      INTEGER NTJ1
      REAL             X, BJ1CS(19), AMPL, THETA, XSML, XMIN, Y,
     1  R1MACH, SCSEVL
C
C SERIES FOR BJ1        ON THE INTERVAL  0.          TO  1.60000E+01
C                                        WITH WEIGHTED ERROR   1.16E-33
C                                         LOG WEIGHTED ERROR  32.93
C                               SIGNIFICANT FIGURES REQUIRED  32.36
C                                    DECIMAL PLACES REQUIRED  33.57
C
      SAVE NTJ1, XSML, XMIN
C
      DATA BJ1CS / -.117261415133327865606240574524003E+0,
     *  -.253615218307906395623030884554698E+0,
     *  +.501270809844695685053656363203743E-1,
     *  -.463151480962508191842619728789772E-2,
     *  +.247996229415914024539124064592364E-3,
     *  -.867894868627882584521246435176416E-5,
     *  +.214293917143793691502766250991292E-6,
     *  -.393609307918317979229322764073061E-8,
     *  +.559118231794688004018248059864032E-10,
     *  -.632761640466139302477695274014880E-12,
     *  +.584099161085724700326945563268266E-14,
     *  -.448253381870125819039135059199999E-16,
     *  +.290538449262502466306018688000000E-18,
     *  -.161173219784144165412118186666666E-20,
     *  +.773947881939274637298346666666666E-23,
     *  -.324869378211199841143466666666666E-25,
     *  +.120223767722741022720000000000000E-27,
     *  -.395201221265134933333333333333333E-30,
     *  +.116167808226645333333333333333333E-32 /
C
      DATA NTJ1, XSML, XMIN / 0, 2*0.E0 /
C     ------------------------------------------------------------------
      IF (NTJ1 .eq. 0) then
         call SINITS (BJ1CS, 19, 0.1E0*R1MACH(3), NTJ1)
         XSML = SQRT (4.0E0*R1MACH(3))
         XMIN = 2.0E0*R1MACH(1)
      endif
C
      Y = ABS(X)
C
      IF (Y .LE. XMIN) THEN
        SBESJ1 = 0.E0
      ELSE IF (Y .LE. XSML) THEN
        SBESJ1 = .5E0 * X
      ELSE IF (Y .LE. 4.E0) THEN
        SBESJ1 = X * ( .25E0 + SCSEVL (.125E0*Y*Y-1.E0,BJ1CS, NTJ1) )
      ELSE
        CALL SBMP1 (Y, AMPL, THETA)
        SBESJ1 = AMPL * COS(THETA)
        IF (X .LT. 0.0e0) SBESJ1 = -SBESJ1
      END IF
C
      RETURN
      END


!-----------------------------------------------------------------------------
!
! QFLAG2IJ
!
! The program QFLAG2IJ reads the flag in dtm_q_output file and store its
! coordinates in cellsID file.
!
!-----------------------------------------------------------------------------

program qflag2ij

use mpar

implicit none

integer(kind=4) i,j,i_basin,cont,ht
integer(kind=4),allocatable :: q_flag(:,:)

call rparfile('hap.in')

allocate(q_flag(N,M))

!call load_dtm("basin_b","basin_i")

open (15,file='dtm_q_output',access='sequential',status='unknown')
open (16,file='cellsID',access='sequential',status='unknown')

write(6,*)
write(6,'(1x,a)') 'Select the header type:'
write(6,'(1x,a)') '0) None'
write(6,'(1x,a)') '1) ESRI ascii file'
write(6,'(1x,a)') '2) GRASS ascii file'
write(6,'(1x,a)') '(Ctrl C to exit)'
write(6,*)
write(6,'(1x,a,$)') '-> '
read(5,*) ht

if (ht /= 0) then
   do i=1,6
      read(15,*)
   end do
end if

do j=M,1,-1
   read(15,*) (q_flag(i,j),i=1,N)
end do

cont=0
do i=1,N
   do j=1,M
      if (q_flag(i,j) == 1) then
         cont=cont+1
         write(16,*) i,j
	 end if
   end do
end do

write(6,*)
write(6,'(1x,i5,1x,a)') cont,'cells selected'
write(6,*)

end program qflag2ij


text:
program problem85
  implicit none
  integer :: besti,bestj,best,temp,i,j

  best=huge(best)
  besti=0
  bestj=0

  do i=1,100
     do j=1,100
        temp=score(i,j)
        if (abs(temp-2000000)<best) then
           besti=i
           bestj=j
           best=abs(temp-2000000)
        end if
     end do
  end do

  print *, best,besti*bestj

contains

  pure function score(x,y)
    implicit none
    integer, intent(in) :: x,y
    integer             :: u,v,i,j,score

    score=0
    do i=1,x
       u=(x+1-i)
       do j=1,y
          v=(y+1-j)
          score=score+(u*v)
       end do
    end do

  end function score

end program problem85

max_stars_repo_name: jamesmcclain/Dioskouroi

max_issues_repo_name: jamesmcclain/Dioskouroi

max_forks_repo_name: jamesmcclain/Dioskouroi


text:
program type_conversion
    use, intrinsic :: iso_fortran_env, only : &
        SP => REAL32, DP => REAL64, &
        I32 => INT32, I64 => INT64
    implicit none
    real(kind=SP) :: x_sp
    real(kind=DP) :: x_dp
    integer(kind=i32) :: i_I32
    integer(kind=i64) :: i_I64

    x_dp = 1.0e100_DP
    x_sp = real(x_dp, kind=SP)
    print *, x_sp

    i_i64 = 2**40
    i_I32 = int(i_I64, kind=I32)
    print *, i_I32
end program type_conversion

max_stars_repo_name: caguerra/Fortran-MOOC

max_issues_repo_name: caguerra/Fortran-MOOC

max_forks_repo_name: caguerra/Fortran-MOOC


text:
subroutine wrd_main(lundia    ,error     ,ndro      ,itdroc    ,runtxt    , &
                  & trifil    ,dtsec     ,gdp       )
!----- GPL ---------------------------------------------------------------------
!
!  Copyright (C)  Stichting Deltares, 2011-2016.
!
!  This program is free software: you can redistribute it and/or modify
!  it under the terms of the GNU General Public License as published by
!  the Free Software Foundation version 3.
!
!  This program is distributed in the hope that it will be useful,
!  but WITHOUT ANY WARRANTY; without even the implied warranty of
!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!  GNU General Public License for more details.
!
!  You should have received a copy of the GNU General Public License
!  along with this program.  If not, see <http://www.gnu.org/licenses/>.
!
!  contact: delft3d.support@deltares.nl
!  Stichting Deltares
!  P.O. Box 177
!  2600 MH Delft, The Netherlands
!
!  All indications and logos of, and references to, "Delft3D" and "Deltares"
!  are registered trademarks of Stichting Deltares, and remain the property of
!  Stichting Deltares. All rights reserved.
!
!-------------------------------------------------------------------------------
!  $Id: wrd_main.f90 6528 2016-09-08 09:35:26Z ye $
!  $HeadURL: https://svn.oss.deltares.nl/repos/delft3d/tags/6686/src/engines_gpl/flow2d3d/packages/io/src/output/wrd_main.f90 $
!!--description-----------------------------------------------------------------
!
!    Function: Main routine for writing the FLOW HIS file.
!
!!--pseudo code and references--------------------------------------------------
! NONE
!!--declarations----------------------------------------------------------------
    use precision
    use dfparall, only: inode, master, parll
    use datagroups
    use netcdf
    !
    use globaldata
    !
    implicit none
    !
    include 'fsm.i'
    include 'tri-dyn.igd'
    !
    type(globdat),target :: gdp
    !
    ! The following list of pointer parameters is used to point inside the gdp structure
    !
    integer(pntrsize)                    , pointer :: namdro
    integer(pntrsize)                    , pointer :: mndro
    integer(pntrsize)                    , pointer :: itdro
    integer(pntrsize)                    , pointer :: ibuff
    integer(pntrsize)                    , pointer :: dxydro
    integer(pntrsize)                    , pointer :: xydro
    integer                              , pointer :: itdrof
    integer                              , pointer :: itdroi
    integer                              , pointer :: itdate
    real(fp)                             , pointer :: dt
    real(fp)                             , pointer :: tunit
    real(fp)                             , pointer :: tzone
    logical                              , pointer :: sferic
    !
    type (datagroup)                     , pointer :: group
    logical                              , pointer :: first
!
! Global variables
!
    integer                                                             , intent(in)  :: ndro   !  Description and declaration in dimens.igs
    integer                                                             , intent(in)  :: itdroc !  Current time counter for the drogue data file
    integer                                                                           :: lundia !  Description and declaration in inout.igs
    logical                                                                           :: error  !!  Flag=TRUE if an error is encountered
    real(fp)                                                                          :: dtsec  !!  Integration time step [in seconds]
    character(30) , dimension(10)                                       , intent(in)  :: runtxt !!  Textual description of model input
    character(*)                                                        , intent(in)  :: trifil !  File name for FLOW NEFIS output
                                                                                !  files (tri"h/m"-"casl""labl".dat/def)
!
! Local variables
!
    integer                                           :: irequest
    integer                                           :: ierror
    integer                                           :: filetype
    integer                                           :: fds
    integer                                , external :: open_datdef
    integer                                , external :: clsnef
    character(256)                                    :: filename
    !
    character(256)                                    :: version_full
    character(8)                                      :: cdate
    character(10)                                     :: ctime
    character(5)                                      :: czone
    character(1024)                                   :: error_string
    !
    character(16) :: simdat  ! Simulation date representing the flow condition at this date
    character(20) :: rundat  ! Execution date of the simulation
    character(6)  :: ftype   ! String containing to which output file version group should be written

!
!! executable statements -------------------------------------------------------
!
    namdro              => gdp%gdr_i_ch%namdro
    mndro               => gdp%gdr_i_ch%mndro
    itdro               => gdp%gdr_i_ch%itdro
    ibuff               => gdp%gdr_i_ch%ibuff
    dxydro              => gdp%gdr_i_ch%dxydro
    xydro               => gdp%gdr_i_ch%xydro
    itdrof              => gdp%gdinttim%itdrof
    itdroi              => gdp%gdinttim%itdroi
    itdate              => gdp%gdexttim%itdate
    dt                  => gdp%gdexttim%dt
    tunit               => gdp%gdexttim%tunit
    tzone               => gdp%gdexttim%tzone
    sferic              => gdp%gdtricom%sferic
    !
    call getdatagroup(gdp, FILOUT_DRO, 'dro-const', group)
    first               => group%first
    !
    call dattim(rundat    )
    simdat(1:16)  = 'yyyymmdd  hhmmss'
    simdat(1:4)   = rundat(1:4)
    simdat(5:6)   = rundat(6:7)
    simdat(7:8)   = rundat(9:10)
    simdat(11:12) = rundat(12:13)
    simdat(13:14) = rundat(15:16)
    simdat(15:16) = rundat(18:19)
    !
    call getfullversionstring_flow2d3d(version_full)
    call date_and_time(cdate, ctime, czone)
    !
    filename = trifil(1:3) // 'd' // trifil(5:)
    filetype = getfiletype(gdp, FILOUT_DRO)
    if (filetype == FTYPE_NETCDF) filename = trim(filename)//'.nc'
    !
    ierror = 0
    do irequest = REQUESTTYPE_DEFINE, REQUESTTYPE_WRITE
       !
       ! request REQUESTTYPE_DEFINE: define all groups, dimensions, and elements
       !         REQUESTTYPE_WRITE : write the data
       !
       if (irequest == REQUESTTYPE_DEFINE) then
          if (.not.first) cycle
          if (inode /= master) cycle
          call delnef(filename,gdp       )
       endif
       !
       ! create or open the file
       !
       if (inode /= master) then
          ! only the master needs to open the file
       elseif (first .and. irequest == REQUESTTYPE_WRITE) then
          ! the file has already been opened in step 1
       elseif (filetype == FTYPE_NEFIS) then
          if (first .and. irequest == REQUESTTYPE_DEFINE) then
             write(lundia,*) 'Creating new '//trim(filename)//'.dat'
             write(lundia,*) 'Creating new '//trim(filename)//'.def'
          endif
          ierror = open_datdef(filename ,fds      , .false.)
          if (ierror /= 0) then
             write(error_string,'(2a)') 'While trying to open dat/def-file ',trim(filename)
             call prtnefiserr(trim(error_string), gdp)
          endif
       elseif (filetype == FTYPE_NETCDF) then
          if (first .and. irequest == REQUESTTYPE_DEFINE) then
             write(lundia,*) 'Creating new '//trim(filename)
             ierror = nf90_create(filename, or(nf90_clobber,nf90_64bit_offset), fds); call nc_check_err(lundia, ierror, "creating file", filename)
             !
             ! global attributes
             !
             ierror = nf90_put_att(fds, nf90_global,  'Conventions', 'CF-1.6'); call nc_check_err(lundia, ierror, "put_att global Conventions", filename)
             ierror = nf90_put_att(fds, nf90_global,  'institution', trim('Deltares')); call nc_check_err(lundia, ierror, "put_att global institution", filename)
             ierror = nf90_put_att(fds, nf90_global,  'references', trim('www.deltares.nl')); call nc_check_err(lundia, ierror, "put_att global references", filename)
             ierror = nf90_put_att(fds, nf90_global,  'source', trim(version_full)); call nc_check_err(lundia, ierror, "put_att global source", filename)
             ierror = nf90_put_att(fds, nf90_global,  'history', &
                    'This file is created on '//cdate(1:4)//'-'//cdate(5:6)//'-'//cdate(7:8)//'T'//ctime(1:2)//':'//ctime(3:4)//':'//ctime(5:6)//czone(1:5)// &
                    ', '//trim('Delft3D')); call nc_check_err(lundia, ierror, "put_att global history", filename)
          else
             ierror = nf90_open(filename, NF90_WRITE, fds); call nc_check_err(lundia, ierror, "opening file", filename)
          endif
       endif
       if (ierror/=0) goto 9999
       !
       ! time independent data
       !
       if (first) then
          call wridro(lundia    ,error     ,trifil    ,ndro      ,itdrof    , &
                    & itdroi    ,simdat    ,itdate    ,tunit     ,dt        , &
                    & ch(namdro),i(mndro)  ,i(itdro)  ,i(ibuff)  ,r(dxydro) , &
                    & irequest  ,fds       ,gdp       )
          if (error) goto 9999
       endif
       !
       ! data per time step
       !
       if (irequest==REQUESTTYPE_DEFINE .or. .not. first) then
          call wrtdro(lundia    ,error     ,trifil    ,itdroc    ,itdrof    , &
                    & itdroi    ,ndro      ,r(xydro)  ,sferic    ,irequest  , &
                    & fds       ,itdate    ,dtsec     ,gdp       )
          if (error) goto 9999
       endif
       !
       if (irequest == REQUESTTYPE_DEFINE) then
          !
          ! upon first request: define all groups, dimensions, and elements on file
          !
          call defnewgrp(fds, FILOUT_DRO, gdp, filename)
          !
          if (filetype == FTYPE_NETCDF) then
             !
             ierror = nf90_enddef(fds); call nc_check_err(lundia, ierror, "enddef", filename)
             if (ierror/=0) goto 9999
             !
          endif
       endif
    enddo
    !
9999 continue
    if (ierror/= 0) error = .true.
    if (inode == master) then
       if (filetype == FTYPE_NETCDF) then
          ierror = nf90_sync(fds); call nc_check_err(lundia, ierror, "sync file", filename)
          ierror = nf90_close(fds); call nc_check_err(lundia, ierror, "closing file", filename)
       elseif (filetype == FTYPE_NEFIS) then
          ierror = clsnef(fds)
          !
          ! wridoc needs to be called AFTER defnewgrp because this routine will
          ! add the dro-version not just to the file, but also to the group
          ! administration in memory and calling defnewgrp afterwards would thus
          ! duplicate definition of the dro-version group.
          !
          if (first) then
             ftype = 'dro'
             call wridoc(error     ,trifil    ,ftype     ,simdat    ,runtxt    , &
                       & .false.   ,''        ,gdp       )
          endif
       endif
       if (ierror/= 0) error = .true.
    endif
    first = .false.
end subroutine wrd_main

max_stars_repo_name: liujiamingustc/phd

max_issues_repo_name: liujiamingustc/phd

max_forks_repo_name: liujiamingustc/phd


module fboxlib

  use iso_c_binding
  use blobjects
  implicit none

contains

  subroutine pybl_cstring(cstr, slen, fstr)
    integer(c_int),    intent(in) :: slen
    character(c_char), intent(in) :: cstr(slen)
    character(len=slen), intent(out) :: fstr

    integer :: i

    do i = 1, slen
       fstr(i:i) = cstr(i)
    end do
  end subroutine pybl_cstring

  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  ! open, close, set_comm

  subroutine pybl_hello() bind(c)
    print *, 'Hello from Fortran, with FBoxLib!'
  end subroutine pybl_hello

  subroutine pybl_open() bind(c)
    use boxlib
    call boxlib_initialize()
  end subroutine pybl_open

  subroutine pybl_close() bind(c)
    use boxlib
    call boxlib_finalize()
  end subroutine pybl_close

  subroutine pybl_mpi_rank(r) bind(c)
    use parallel
    integer(c_int), intent(out) :: r
    r = parallel_myproc()
  end subroutine pybl_mpi_rank

  subroutine pybl_mpi_size(r) bind(c)
    use parallel
    integer(c_int), intent(out) :: r
    r = parallel_nprocs()
  end subroutine pybl_mpi_size

  subroutine pybl_mpi_reduce_max(r) bind(c)
    use parallel
    real(c_double), intent(inout) :: r
    real(c_double) :: s
    s = r
    call parallel_reduce(r, s, MPI_MAX)
  end subroutine pybl_mpi_reduce_max

  ! subroutine pybl_set_comm(comm)
  !   use parallel
  !   implicit none
  !   integer, intent(in) :: comm
  !   call parallel_set_m_comm(comm)
  ! end subroutine pybl_set_comm

  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  ! boxarray routines

  subroutine pybl_boxarray_create_from_boxes(boxes,nboxes,dim,cptr) bind(c)
    integer(c_int), intent(in   ), value :: dim, nboxes
    integer(c_int), intent(in   )        :: boxes(dim,2,nboxes)
    type(c_ptr),    intent(  out)        :: cptr

    type(boxarray), pointer :: ba
    integer   :: i
    type(box) :: bs(nboxes)

    do i=1,nboxes
       bs(i) = make_box(boxes(:,1,i), boxes(:,2,i))
    end do

    call pybl_boxarray_new(cptr,ba)
    call build(ba, bs)
  end subroutine pybl_boxarray_create_from_boxes

  subroutine pybl_boxarray_print(cptr) bind(c)
    type(c_ptr), intent(in   ), value :: cptr
    type(boxarray), pointer :: ba
    call pybl_boxarray_get(cptr, ba)
    call print(ba)
  end subroutine pybl_boxarray_print

  subroutine pybl_boxarray_nboxes(cptr,nb) bind(c)
    type(c_ptr),    intent(in   ), value :: cptr
    integer(c_int), intent(  out)        :: nb
    type(boxarray), pointer :: ba
    call pybl_boxarray_get(cptr, ba)
    nb = nboxes(ba)
  end subroutine pybl_boxarray_nboxes

  subroutine pybl_boxarray_dim(cptr,dim) bind(c)
    type(c_ptr),    intent(in   ), value :: cptr
    integer(c_int), intent(  out)        :: dim
    type(boxarray), pointer :: ba
    call pybl_boxarray_get(cptr, ba)
    dim = ba%dim
  end subroutine pybl_boxarray_dim

  subroutine pybl_boxarray_maxsize(cptr,maxsize,ndim) bind(c)
    type(c_ptr),    intent(in   ), value :: cptr
    integer(c_int), intent(in   ), value :: ndim
    integer(c_int), intent(in   )        :: maxsize(ndim)
    type(boxarray), pointer :: ba
    call pybl_boxarray_get(cptr, ba)
    call boxarray_maxsize(ba, maxsize)
  end subroutine pybl_boxarray_maxsize

  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  ! layout routines

  subroutine pybl_layout_create_from_boxarray(bacptr,dim,pmask,cptr) bind(c)
    type(c_ptr),    intent(in   ),  value :: bacptr
    type(c_ptr),    intent(  out)         :: cptr
    integer(c_int), intent(in   ),  value :: dim
    integer(c_int), intent(in   )         :: pmask(dim)

    type(boxarray), pointer :: ba
    type(layout), pointer   :: la
    logical                 :: lpmask(dim)

    call pybl_boxarray_get(bacptr, ba)
    call pybl_layout_new(cptr, la)

    lpmask = pmask == 1
    call build(la, ba, boxarray_bbox(ba), pmask=lpmask)
  end subroutine pybl_layout_create_from_boxarray

  subroutine pybl_layout_print(cptr) bind(c)
    type(c_ptr), intent(in   ), value :: cptr
    type(layout), pointer :: la
    call pybl_layout_get(cptr, la)
    call print(la)
  end subroutine pybl_layout_print

  subroutine pybl_layout_nboxes(cptr, boxes) bind(c)
    type(c_ptr),    intent(in   ), value :: cptr
    integer(c_int), intent(  out)        :: boxes
    type(layout), pointer :: la
    call pybl_layout_get(cptr,la)
    boxes = nboxes(la)
  end subroutine pybl_layout_nboxes

  subroutine pybl_layout_get_box(cptr, intbuf, i) bind(c)
    type(c_ptr),    intent(in   ), value :: cptr, intbuf
    integer(c_int), intent(in   ), value :: i
    type(layout), pointer :: la
    integer, pointer :: buf(:)
    type(box) :: lbx
    call c_f_pointer(intbuf, buf, [ 7 ])
    call pybl_layout_get(cptr,la)
    lbx = get_box(la, i)
    buf(1) = lbx%dim
    buf(2:4) = lbx%lo
    buf(5:7) = lbx%hi
  end subroutine pybl_layout_get_box

  subroutine pybl_layout_local(cptr, nbox, islocal) bind(c)
    type(c_ptr),    intent(in   ), value :: cptr
    integer(c_int), intent(in   ), value :: nbox
    integer(c_int), intent(  out)        :: islocal
    type(layout), pointer :: la

    call pybl_layout_get(cptr,la)
    islocal = 0
    if (local(la, nbox)) then
       islocal = 1
    end if
  end subroutine pybl_layout_local

  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  ! multifab routines

  subroutine pybl_lmultifab_create_from_layout(la_cptr,cptr) bind(c)
    type(c_ptr),    intent(in   ), value :: la_cptr
    type(c_ptr),    intent(  out)        :: cptr

    type(layout), pointer   :: la
    type(lmultifab), pointer :: mfab

    call pybl_layout_get(la_cptr,la)
    call pybl_lmultifab_new(cptr,mfab)

    call build(mfab, la)
    call setval(mfab, .false.)
  end subroutine pybl_lmultifab_create_from_layout

  subroutine pybl_multifab_create_from_layout(la_cptr,nc,ng,cptr) bind(c)
    type(c_ptr),    intent(in   ), value :: la_cptr
    integer(c_int), intent(in   ), value :: nc, ng
    type(c_ptr),    intent(  out)        :: cptr

    type(layout), pointer   :: la
    type(multifab), pointer :: mfab

    call pybl_layout_get(la_cptr,la)
    call pybl_multifab_new(cptr,mfab)

    call build(mfab, la, nc=nc, ng=ng) !, interleave=interleave)
    call setval(mfab, 0.0d0)
  end subroutine pybl_multifab_create_from_layout

  subroutine pybl_multifab_print(cptr) bind(c)
    type(c_ptr), intent(in   ), value :: cptr
    type(multifab), pointer :: mfab
    call pybl_multifab_get(cptr, mfab)
    call print(mfab)
  end subroutine pybl_multifab_print

  subroutine pybl_multifab_destroy(cptr) bind(c)
    type(c_ptr), intent(in   ), value :: cptr
    type(multifab), pointer :: mfab
    call pybl_multifab_get(cptr, mfab)
    call destroy(mfab)
  end subroutine pybl_multifab_destroy

  subroutine pybl_multifab_fill_boundary(cptr) bind(c)
    type(c_ptr), intent(in   ), value :: cptr
    type(multifab), pointer :: mfab
    call pybl_multifab_get(cptr, mfab)
    call fill_boundary(mfab)
  end subroutine pybl_multifab_fill_boundary

  subroutine pybl_multifab_write(cptr, dirname, dlen, header, hlen) bind(c)
    type(c_ptr),       intent(in   ), value :: cptr
    integer(c_int),    intent(in   ), value :: dlen, hlen
    character(c_char), intent(in   )        :: dirname(dlen)
    character(c_char), intent(in   )        :: header(hlen)

    type(multifab), pointer :: mfab
    character(len=dlen) :: dname
    character(len=hlen) :: hname

    call pybl_cstring(dirname, dlen, dname)
    call pybl_cstring(header, hlen, hname)

    call pybl_multifab_get(cptr, mfab)
    call fabio_write(mfab, dname, hname)
  end subroutine pybl_multifab_write

  subroutine pybl_multifab_read(dirname, dlen, header, hlen, cptr) bind(c)
    integer(c_int),    intent(in   ), value :: dlen, hlen
    character(c_char), intent(in   )  :: dirname(dlen)
    character(c_char), intent(in   )  :: header(hlen)
    type(c_ptr),       intent(  out) :: cptr

    type(multifab), pointer :: mfab
    character(len=dlen) :: dname
    character(len=hlen) :: hname

    call pybl_cstring(dirname, dlen, dname)
    call pybl_cstring(header, hlen, hname)
    call pybl_multifab_new(cptr, mfab)
    call fabio_multifab_read_d(mfab, dname, hname)
  end subroutine pybl_multifab_read

  subroutine pybl_multifab_copy(dcptr, scptr) bind(c)
    type(c_ptr), intent(in   ), value :: dcptr, scptr
    type(multifab), pointer :: dmfab, smfab
    call pybl_multifab_get(dcptr, dmfab)
    call pybl_multifab_get(scptr, smfab)
    call copy(dmfab, 1, smfab, 1, ncomp(smfab))
  end subroutine pybl_multifab_copy

  subroutine pybl_multifab_info(cptr, dim, nboxes, nc, ng) bind(c)
    type(c_ptr),    intent(in   ), value :: cptr
    integer(c_int), intent(  out)        :: dim, nboxes, nc, ng
    type(multifab), pointer :: mfab
    call pybl_multifab_get(cptr,mfab)
    dim = mfab%dim
    nboxes = nfabs(mfab)
    nc = mfab%nc
    ng = mfab%ng
  end subroutine pybl_multifab_info

  subroutine pybl_multifab_layout(cptr, laptr) bind(c)
    type(c_ptr), intent(in   ), value :: cptr
    type(c_ptr), intent(  out)        :: laptr
    type(multifab), pointer :: mfab
    type(layout), pointer :: la
    call pybl_multifab_get(cptr,mfab)
    la => mfab%la
    laptr = c_loc(la)
  end subroutine pybl_multifab_layout

  subroutine pybl_lmultifab_info(cptr, dim, nboxes, nc, ng) bind(c)
    type(c_ptr),    intent(in   ), value :: cptr
    integer(c_int), intent(  out)        :: dim, nboxes, nc, ng
    type(lmultifab), pointer :: mfab
    call pybl_lmultifab_get(cptr,mfab)
    dim = mfab%dim
    nboxes = nfabs(mfab)
    nc = mfab%nc
    ng = mfab%ng
  end subroutine pybl_lmultifab_info

  subroutine pybl_fab_info(cptr, nbox, dim, nc, bx_lo, bx_hi, pbx_lo, pbx_hi, ibx_lo, ibx_hi) bind(c)
    type(c_ptr),            intent(in   ), value        :: cptr
    integer(c_int),         intent(in   ), value        :: nbox
    integer(c_int),         intent(  out)               :: dim, nc
    integer(c_int),         intent(  out), dimension(3) :: bx_lo, bx_hi, pbx_lo, pbx_hi, ibx_lo, ibx_hi

    type(multifab), pointer :: mfab
    type(box) :: bx

    call pybl_multifab_get(cptr,mfab)
    dim = get_dim(mfab%fbs(nbox))
    nc = ncomp(mfab%fbs(nbox))

    bx = get_box(mfab%fbs(nbox))
    bx_lo = bx%lo
    bx_hi = bx%hi

    bx = get_pbox(mfab%fbs(nbox))
    pbx_lo = bx%lo
    pbx_hi = bx%hi

    bx = get_ibox(mfab%fbs(nbox))
    ibx_lo = bx%lo
    ibx_hi = bx%hi
  end subroutine pybl_fab_info

  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  ! regrid

  subroutine pybl_regrid(lacptrs, mfcptrs, nlevs, max_levs, tag_boxes_cb) bind(c)
    use bc_module
    use define_bc_module
    use regrid_module

    integer(c_int), intent(in   ), value :: max_levs
    integer(c_int), intent(inout)        :: nlevs
    type(c_ptr),    intent(inout)        :: lacptrs(max_levs), mfcptrs(max_levs)
    type(c_ptr),    intent(in   ), value :: tag_boxes_cb

    type(multifab), pointer :: mf
    type(layout), pointer   :: la
    type(multifab), target  :: mfs(max_levs)

    integer           :: i, n, dim, new_nlevs, phys_bc_in(3,2)
    type(ml_boxarray) :: mba
    type(ml_layout)   :: mla
    type(bc_tower)    :: bct
    real(dp_t)        :: dx(max_levs)

    ! build multilevel layout
    call pybl_layout_get(lacptrs(1),la)
    dim = get_dim(la)

    call build(mba, max_levs, dim)
    mba%rr = 2
    dx = 1
    do i=1,max_levs-1
       dx(i+1) = dx(i)/2
    end do

    mba%pd(1) = get_box(la, 1)
    do n=2,max_levs
       mba%pd(n) = refine(mba%pd(n-1),mba%rr((n-1),:))
    enddo

    do i=1,nlevs
       call pybl_layout_get(lacptrs(i),la)
       call pybl_multifab_get(mfcptrs(i),mf)
       mba%bas(i) = get_boxarray(la)
       call boxarray_verify_dim(mba%bas(i))
       mfs(i) = mf
    end do

    call ml_layout_restricted_build(mla,mba,nlevs)
    ! call destroy(mba)

    phys_bc_in = PERIODIC
    call bc_tower_init(bct, max_levs, dim, phys_bc_in)
    do n = 1,nlevs
       call bc_tower_level_build(bct,n,mla%la(n))
    end do

    call regrid(mla, mfs, nlevs, max_levs, dx, bct, 2, 64, tag_boxes_cb)

    do i=1,nlevs
       call pybl_multifab_new(mfcptrs(i),mf)
       mf = mfs(i)
    end do

  end subroutine pybl_regrid

end module fboxlib