text:
! (C) Copyright 2005- ECMWF.
!
! This software is licensed under the terms of the Apache Licence Version 2.0
! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
!
! In applying this licence, ECMWF does not waive the privileges and immunities granted to it by
! virtue of its status as an intergovernmental organisation nor does it submit to any jurisdiction.
!
!
!  Description: How to create a new GRIB message from a sample.
!
!
program sample
   use eccodes
   implicit none
   integer  :: err
   integer  :: outfile, datafile
   integer  :: igribsample, igribclone, igribdata, size1
   integer  :: date1, startStep, endStep, table2Version, indicatorOfParameter
   integer  :: decimalPrecision
   character(len=10) stepType
   real(kind=8), dimension(:), allocatable   :: v1, v2, v

   date1 = 20080104
   startStep = 0
   endStep = 12
   stepType = 'accum'
   table2Version = 2
   indicatorOfParameter = 61
   decimalPrecision = 2

   ! a new GRIB message is loaded from an existing sample.
   ! Samples are searched in a default sample path (use codes_info
   ! to see where that is). The default sample path can be changed by
   ! setting the environment variable ECCODES_SAMPLES_PATH
   call codes_grib_new_from_samples(igribsample, "regular_latlon_surface.grib1")

   call codes_open_file(outfile, 'f_out.samples.grib1', 'w')
   call codes_open_file(datafile, '../../data/tp_ecmwf.grib', 'r')

   call codes_grib_new_from_file(datafile, igribdata, err)

   call codes_get_size(igribdata, 'values', size1)
   allocate (v(size1))
   allocate (v1(size1))
   allocate (v2(size1))

   call codes_get(igribdata, 'values', v)

   v = v*1000.0 ! different units for the output grib
   v1 = v

   do while (err /= CODES_END_OF_FILE)

      call codes_clone(igribsample, igribclone) ! clone sample before modifying it

      call codes_set(igribclone, 'dataDate', date1)
      call codes_set(igribclone, 'table2Version', table2Version)
      call codes_set(igribclone, 'indicatorOfParameter', indicatorOfParameter)

      call codes_set(igribclone, 'stepType', stepType)
      call codes_set(igribclone, 'startStep', startStep)
      call codes_set(igribclone, 'endStep', endStep)

      call codes_set(igribclone, 'decimalPrecision', decimalPrecision)

      call codes_set(igribclone, 'values', v)

      call codes_write(igribclone, outfile)

      call codes_grib_new_from_file(datafile, igribdata, err)

      if (err == 0) then
         call codes_get(igribdata, 'values', v2)

         v2 = v2*1000.0 ! different units for the output grib

         v = v2 - v1 ! accumulation from startStep to endStep

         v1 = v2 ! save previous step field

         startStep = startStep + 12
         endStep = endStep + 12

      end if

   end do

   call codes_release(igribsample)
   deallocate (v)
   deallocate (v1)
   deallocate (v2)

   call codes_close_file(outfile)

end program sample

max_stars_repo_name: mtasaka/eccodes

max_issues_repo_name: mtasaka/eccodes

max_forks_repo_name: mtasaka/eccodes


text:
      module tfmessage
      use tfstk

      contains
      type (sad_descriptor) function kxmessagename(mess)
      implicit none
      type (sad_descriptor) k1
      type (sad_symdef), pointer :: symd
      integer*4 i,l,isp1
      character*(*) ,intent(in):: mess
      i=index(mess,'::')
      if(i .le. 0)then
        write(*,*)'itfmessagename implementation error: ',mess
        call abort
      endif
      isp1=isp+1
      ktastk(isp1)=ktfoper+mtfmessagename
      k1=kxsymbolz(mess(1:i-1),i-1,symd)
      if(ktfoperq(symd%value))then
        dtastk(isp1+1)=symd%value
      else
        dtastk(isp1+1)=k1
      endif
      l=len(mess)
      isp=isp1+2
      dtastk(isp)=kxsymbolz(mess(i+2:l),l-i-1)
      kxmessagename=kxcomposev(isp1)
      isp=isp1-1
      return
      end function

      end module

      integer*4 function itfmessage(level,mess0,arg0)
      use tfstk
      use tfmessage
      use tfcsi
      use eeval
      implicit none
      type (sad_descriptor) dm,ks,mn
      type (sad_dlist), pointer :: klx,klm,klhm,klhms
      integer*4 ,intent(in):: level
      integer*4 irtc,ltr0,l,la,isp0
      character*(*) ,intent(in):: mess0,arg0
      character*256 mess,arg
      logical*4 ,save::iter=.false.
      data mn%k /0/
      if(iter .or. level .lt. ierrorth .or.
     $     rlist(iaximmediate) .lt. 0.d0)then
        itfmessage=-1
        return
      endif
      if(mn%k .eq. 0)then
        mn=kxsymbolz('MessageString',13)
      endif
      iter=.true.
      ltr0=ltrace
      ltrace=0
      l=len_trim(mess0)
      mess=mess0(:l)
      la=len_trim(arg0)
      arg=arg0(:la)
      dm=dtfcopy1(kxmessagename(mess(:l)))
      call descr_sad(dm,klm)
      ks=tfleval(klm,.true.,irtc)
      if(irtc .ne. 0)then
        if(irtc .gt. 0)then
          call tfaddmessage(' ',0,icslfno())
        endif
        itfmessage=-1
        call tflocal1(dm%k)
      elseif(ktfstringq(ks))then
        call tflocal(kerror)
        kerror=ktflist+ktadaloc(0,4,klx)
        klx%rbody(1)=dble(level)
        klx%dbody(2)=kxadaloc(0,1,klhm)
        klx%dbody(3)=kxsalocb(0,mess,l)
        klx%dbody(4)=kxadaloc(0,1,klhms)
        klhm%head%k=ktfoper+mtfhold
        klhm%dbody(1)=dm
        klhms%head%k=ktfoper+mtfhold
        isp0=isp
        isp=isp+1
        dtastk(isp)=mn
        isp=isp+1
        dtastk(isp)=ks
        call tfstringliststk(arg(:la))
        klhms%dbody(1)=dtfcopy1(kxcomposev(isp0+1))
        isp=isp0
        rlist(ierrorgen)=rlist(ierrorgen)+1.d0
        ierrorprint=sad_loc(klx%head%k)
        itfmessage=1
      else
        call tflocal1(dm%k)
        itfmessage=-1
      endif
      iter=.false.
      ltrace=ltr0
      return
      end

c     tfstringliststk destroy `string' to decode backslash escape
      subroutine tfstringliststk(string)
      use tfstk
      implicit none
      character*(*) ,intent(inout):: string
      integer*4 l,is,is1,in
      l=len_trim(string)
c      write(*,*)'tfstringliststk ',l,' ',string(1:l)
      is=1
      do while(is .le. l)
c     Search '"' from string(is:l) with backslash escape
        is1=is
        do while(is .le. l .and. string(is1:is1) .ne. '"')
           if(string(is1:is1) .eq. '\\') then
              string(is1:l-1)=string(is1+1:l)
              l=l-1
           endif
           is1=is1+1
        enddo
        if(is1 .gt. l) then
           is1=0
        endif
        if(is1 .le. 0)then
          isp=isp+1
          dtastk(isp)=kxsalocb(-1,string(is:),l-is+1)
          return
        endif
        is=is1+1
c     Search '"' from string(is+1:l) with backslash escape
        in=is
        do while(in .le. l .and. string(in:in) .ne. '"')
           if(string(in:in) .eq. '\\') then
              string(in:l-1)=string(in+1:l)
              l=l-1
           endif
           in=in+1
        enddo
        in=in-is
        if(in .gt. l) then
           in=0
        endif
        if(in .le. 0)then
          isp=isp+1
          dtastk(isp)=kxsalocb(-1,string(is:),l-is+1)
          return
        endif
        isp=isp+1
        dtastk(isp)=kxsalocb(-1,string(is:),in)
        is=is+in+1
      enddo
      return
      end

      integer*4 function itfmessageexp(level,mess,k)
      use tfstk
      implicit none
      type (sad_descriptor) ,intent(in):: k
      integer*4 ,intent(in):: level
      integer*4 nc,itfmessage
      character*(*) ,intent(in):: mess
      character*256 buf
      call tfconvstrs(buf(2:),k,nc,.true.,' ')
      nc=max(0,min(254,nc))
      buf(1:1)='"'
      buf(nc+2:nc+2)='"'
      itfmessageexp=itfmessage(level,mess,buf(:nc+2))
      return
      end

      integer*4 function itfmessagestr(level,mess,str)
      use tfstk
      use strbuf
      implicit none
      integer*4 ,parameter :: lbuf=128
      type (sad_strbuf) , pointer :: strb
      character*(*) ,intent(in):: str
      integer*4 ,intent(in):: level
      character*(*) ,intent(in):: mess
      integer*4 irtc,isp0,itfmessage
      isp0=isp
      call getstringbuf(strb,lbuf*2,.true.)
      call tfquotestring(strb,str,min(lbuf-1,len(str)),0,irtc)
      itfmessagestr=itfmessage(level,mess,strb%str(1:strb%nch))
      isp=isp0
      return
      end

      subroutine tferrorhandle(kx,irtc)
      use tfstk
      use tfcsi
      implicit none
      type (sad_descriptor) ,intent(in):: kx
      integer*4 ,intent(in):: irtc
      if(ierrorprint .ne. 0)then
        ierrorf=kx%k
        if(irtc .gt. 0)then
          if(tflistq(kerror))then
            call tfaddmessage(' ',0,icslfno())
          endif
        endif
        ierrorprint=0
      endif
      return
      end

      subroutine tferrorhandles(kx,str,irtc)
      use tfstk
      use tfcsi
      implicit none
      type (sad_descriptor) ,intent(in):: kx
      character*(*) ,intent(in):: str
      integer*4 ,intent(in):: irtc
      if(ierrorprint .ne. 0)then
        ierrorf=kx%k
        if(irtc .gt. 0)then
          if(tflistq(kerror))then
            call tfaddmessage(str(1:len_trim(str)),
     $           len_trim(str),icslfno())
          endif
        endif
        ierrorprint=0
      endif
      return
      end

      subroutine tfaddmessage(str,ip,lfno)
      use tfstk
      use efun
      implicit none
      type (sad_descriptor) kx,kxaddmess
      type (sad_dlist), pointer :: kle
      integer*4 ,intent(in):: ip,lfno
      integer*4 irtc,isp1,ltr0
      character*(*) ,intent(in):: str
      logical*4 iter
      data kxaddmess%k,iter /0,.false./
      if(iter)then
        call tfdebugprint(kerror,'???Error in error handling.',10)
        return
      endif
      if(.not. tflistq(kerror,kle))then
        return
      endif
      iter=.true.
      ltr0=ltrace
      ltrace=0
      if(kxaddmess%k .eq. 0)then
        kxaddmess=kxsymbolz('Add$Message',11)
      endif
      isp1=isp+1
      isp=isp1+1
      dtastk(isp1)=kxaddmess
      ktastk(isp)=kerror
c      call tfdebugprint(kerror,'addmessage',1)
      kx=tfefunref(isp1,.false.,irtc)
c      call tfdebugprint(kx,'addmessage',1)
      isp=isp1-1
      if(irtc .eq. 0)then
        if(ktflistq(kx))then
          call tfemes(kx,str,ip,lfno)
        else
          call tflocal(kerror)
          kerror=0
          call tclrfpe
        endif
      endif
      ierrorf=0
      ierrorprint=0
      iter=.false.
      ltrace=ltr0
      return
      end

      subroutine tfemes(k,string,ip,lfno)
      use tfstk
      implicit none
      type (sad_descriptor) ,intent(in):: k
      type (sad_dlist), pointer :: kl
      character*(*) ,intent(in):: string
      integer*4 ,intent(in):: ip,lfno
      integer*4 ls,ifchar,is,i,lb,nc,itfgetrecl,nc1
      character*10 autofg
      character*512 buff,tfconvstr,tfgetstr
      if(ktfaddr(k) .eq. 0 .or. ierrorprint .eq. 0)then
        ierrorf=0
        return
      elseif(tflistq(k,kl))then
        buff(1:3)='???'
        buff(4:)=tfgetstr(kl%dbody(3),nc)
        buff(nc+4:nc+6)=': '
        buff(nc+7:)=tfgetstr(kl%dbody(4),nc1)
        call capita(buff(nc+7:nc+7))
      else
        buff='???Unexpected error: '//autofg(dble(k%k),'S10')
      endif
      ierrorprint=0
      lb=min(len_trim(buff),len(buff)-1)
      if(ierrorf .ne. 0)then
        buff(lb+1:)=' in '//tfconvstr(ierrorf,nc,'*')
        if(lb+nc+4 .gt. len(buff))then
          buff(len(buff)-3:)=' ...'
        endif
        lb=min(lb+nc+4,len(buff))
      else
        lb=lb+1
        buff(lb:lb)=':'
      endif
      call tftruncprint(buff(1:lb),itfgetrecl()-1,' ,})]'//char(10),
     $     .false.,lfno)
      if(lfno .ne. 6 .and. lfno .ne. 0)then
        call tftruncprint(buff(1:lb),itfgetrecl()-1,' ,})]'//char(10),
     $       .false.,6)
      endif
      if(ip .gt. 0)then
        do i=ip-2,1,-1
          if(string(i:i) .eq. char(10))then
            is=i+1
            go to 1
          endif
        enddo
        is=max(1,ip-32)
 1      ls=ifchar(string(1:ip),char(10),is)-1
        if(ls .le. 0)then
          ls=len_trim(string(1:ip))-1
        endif
        if(ls .ge. is)then
          write(lfno,'(a)')string(is:ls)
          if(lfno .ne. 6)then
            write(6,'(a)')string(is:ls)
          endif
          if(ierrorf .eq. 0)then
            ls=ip-is
            if(ls .le. 0)then
              ls=1
            endif
            buff(1:ls)=' '
            buff(ls:ls)='^'
            write(lfno,'(a)')buff(1:ls)
            if(lfno .ne. 6)then
              write(6,'(a)')buff(1:ls)
            endif
          endif
        endif
      endif
      ierrorf=0
      call tflocal(kerror)
      kerror=0
      call tclrfpe
      return
      end

      subroutine tfreseterror
      use tfstk
      implicit none
      if(kerror .ne. 0)then
        call tflocal(kerror)
        kerror=0
      endif
      if(ierrorprint .ne. 0)then
        ierrorprint=0
      endif
      return
      end

      subroutine tffserrorhandle(l,irtc)
      use tfstk
      implicit none
      integer*4 l,irtc
      call tferrorhandle(sad_descr(dble(l)),irtc)
      call tfreseterror
      return
      end

ext: f

max_stars_repo_name: noboruatkek/SAD

max_issues_repo_name: noboruatkek/SAD

max_forks_repo_name: noboruatkek/SAD


text:
Hay (aka Straw) is a multipurpose material grown locally. It can be used for pet/horse bedding, building material like the Straw Bale Bench, landscaping mulch, bicycle race crash barriers, outdoor seating and all sorts of other stuff.

Where to buy?
  http://www.higbyfeed.com/ Higbys Feed and Farm Supply, Dixon
  Windmill Feed off County Road 95 in Woodland.
  House on 8th street between B & F.  To find this place, just drive your truck down the street yelling, Hay, hay! out the window until the farmer comes out to sell some to you. you cant miss it its 45 houses in on 8th street just West of F, a big stack of hay with a sign that says HAY for SALE on it...

  Small amounts (for small mammal pets) can be obtained at PetCo.  (Maybe Ace?)

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


program gaussseli
    ! This program simply takes in an nxn matrix, and an nx1 matrix, converts into a reduced row echelon form using gauss jordan pivoting,
    ! ,then calculates the determinant of nxn, if determinant is non-singular the program proceeds to calculate inverse
    ! and solution of (nXn) x = (nX1) form equation (i.e. Ax=b)
    !Various parts of the code are highlighted by comments as to what they do.
    implicit none
    integer :: n,i,k,j
    real , allocatable :: A(:,:), temp(:)
    real :: q, ep=10**(-7)
    print*, "Give the dimension of matrix A in form AX=b"
    read*, n
    allocate(A(n,n+1+n), temp(n+n+1))   !Allocate nxn for original matrix, nx1 for b, nxn for I all concatenated
    !!!!!!!!!!!!!!!!!!!!!!!!!!!! READING MATRIX FROM USER !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    print*, "Give the rows of matrix A"
    do i=1,n
        read*, A(i,1:n)
    enddo
    print*,"Give values of b"
    do i=1,n
        read*, A(i,n+1)
    enddo
    !!!!!!!!!!!!! MATRIX FORM IS A(augmented) = A(read)|b(read)|I CONCATENATED AND PROCESSED TOGETHER THROUGHOUT !!!!!!!!!!!!!!
    do i=1,n
        A(i,n+i+1) =1   !getting the identity of the third part in matrix
    enddo
    print*, "MATRIX = A | b | Identity"
    do i=1,n
        print*,A(i,:)
    enddo
    print*, "================================================"
    !!!!!!!!!!!!!!!! CONVERTING THE MATRIX TO UPPER TRIANGULAR REDUCED ECHELON !!!!!!!!!!!!!!!!!!!
    do k = 1,n
        do i=k,n
            if (abs(A(i,k)) .gt. abs(A(k,k))) then
                temp = A(k,:) ; A(k,:) = A (i,:); A(i,:) = temp    !swapping higher terms to higher position in matrix
            endif
        enddo
        do j=k+1,n
            q=(A(j,k))/A(k,k)
            do i=k,n+n+1
                A(j,i) = A(j,i) - q* A(k,i)       !finally making the lower triangular terms zero
            enddo
        enddo
    enddo
    !!!!!!!!!!!!!!!!!!!!!!!!! CONVERSION DONE !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    print*, "The modified/augmented matrix:"
    do i=1,n
        print*,A(i,:)
    enddo
    print*, "====================================================="
    !!!!!!!!!! CHECK DETERMINANT AFTER REDUCING IT TO UPPER TRIANGULAR !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    print*, "Checking determinant"
    q = 1    ! determinant
    do i=1,n
        q = q * A(i,i)      !determinant is simply the product of diagonal elements in upper diagonal
    enddo
    print*, "The determinant is:", q
    If (abs(q) .le. ep) then
        print*, "Determinant is SINGULAR, unsolvable"
        STOP
    endif
    print*, "====================================================="
    !!!!!!!!!!!!!!!!!!!!!!!!!!!CALCULATE FINAL RESULT!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    temp =0
    temp(n) = (A(n,n+1))/A(n,n)
    do i = n-1,1,-1
        q=0
        do j=i+1,n
            q = q + A(i,j)*temp(j)
        enddo
        temp(i)= (A(i,n+1) -q)/A(i,i)
    enddo
    print*, "The solutions are", temp(1:n)
    print*, "========================================================="
    !!!!!!!!!!!!!!!!!!!!!!!!!!!! COMPUTE INVERSE NOW !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    print*, "The inverse of A is"
    do i = n,1,-1
        A(i,:) = A(i,:)/A(i,i)
    enddo
    do i = n,1,-1
        do j=1,i-1
            A(j,:) = A(j,:) - A(i,:)*A(j,i)
        enddo
    enddo
    do i=1,n
        print*,A(i,n+2:)
    enddo
    print*, "========================================================="
end program gaussseli


NOT USED NOW
module subtools

    use constants_module
    use arrays_module
    use var_module
    use arrays_section_module
    use xsec_attribute_module

contains

 !   subroutine calc_q_sk_multi(i,currentQ,q_sk_multi)
 !       implicit none
 !
 !       integer,intent(in) :: i
 !       real,intent(in) :: currentQ
 !       real,intent(out) :: q_sk_multi
 !       integer :: pp, tableLength
 !       real :: r_interpo_nn

 !       q_sk_multi = 1.0
 !       do pp = 1, noQSKtable
 !           if (  ( eachQSKtableNodeRange(1,pp) - i) * ( eachQSKtableNodeRange(2,pp) - i) .le. 0 ) then
 !               tableLength = Q_sk_tableEntry(pp)
 !               q_sk_multi = r_interpo_nn(Q_Sk_Table(1,1:tableLength,pp), &
 !               Q_Sk_Table(2,1:tableLength,pp),tableLength,currentQ)
 !           end if
 !       end do

 !   end subroutine

!+++-------------------------------------------------------------------
!+ Compute Q from a given value of Y using Manning's equation, which
!+ becomes uniform flow eqns when slp is So and nonuniform flow eqns
!+ when slp is Sf.
!+++-------------------------------------------------------------------
    subroutine ManningEq_QforY(y_mn, bo_mn, slp, n_mn, q_mn)

        implicit none

        real,intent(in) :: y_mn, bo_mn, slp, n_mn
        real,intent(out) :: q_mn
        real :: ar, peri, hydr
        integer :: chshp

        !+++----------------------------------------------------+
        !+ Rectangular channels
        !+++----------------------------------------------------+
        if (chshp.eq.1) then
            !* area
            ar=y_mn*bo_mn
            !* perimeter
            peri=2.0*y_mn+bo_mn
            !* hydraulic radius
            hydr=ar/peri
            !* Note that n_mn=1/Mannin's N.
            !q_mn=n_mn*ar*(hydr**(2.0/3.0))*sqrt(slp)
            q_mn=n_mn*ar*(hydr**(2.0/3.0))*(slp**0.5)
        end if

    end subroutine ManningEq_QforY

!+++-------------------------------------------------------------------------------
!+ Compute Y and thus area from a given value of Q using Manning's equation, which
!+ becomes uniform flow eqns when slp is So and nonuniform flow eqns
!+ when slp is Sf. The numerical method used in this subroutine is described in
!+ p.88, Chaudhry book.
!+++-------------------------------------------------------------------------------
    !subroutine ManningEq_YforQ()!(y_mn, bo_mn, slp, n_mn, q_mn)

    !end subroutine ManningEq_YforQ

!+++-------------------------------------------------------------------------------
!+ computation of normal depth in regular/trapezoidal x-section using
!+ Newton-Raphson method. Refer to Appendix C-2, Chaudhary and p71,RM1_MESH
!+++-------------------------------------------------------------------------------
    subroutine normal_crit_y(i, q_sk_multi, So, dsc, y_norm, y_crit, area_n, area_c)

        implicit none

        integer, intent(in) :: i
        real, intent(in) :: q_sk_multi, So, dsc
        real, intent(out) :: y_norm, y_crit, area_n, area_c
        real :: area_0, width_0, errorY, hydR_0, r_interpol!, fro
        integer :: trapnm_app, recnm_app, iter

        !c1=manN*dsc/(c0*sqrt(So))

        !y_norm=1.0    !*initial estimate

            elevTable = xsec_tab(1,:,i)
            areaTable = xsec_tab(2,:,i)
            rediTable = xsec_tab(4,:,i)
            topwTable = xsec_tab(6,:,i)
            area_0 = r_interpol(elevTable,areaTable,nel,oldY(i))
            width_0= r_interpol(elevTable,topwTable,nel,oldY(i))

            area_c=area_0
            errorY = 100.
            !print*, 'inside function', i, q_sk_multi, So, dsc, area_0, width_0
            !pause
            do while (errorY .gt. 0.00001)

                hydR_0 = r_interpol(areaTable,rediTable,nel,area_0)
                area_n = dsc/sk(i)/q_sk_multi/ hydR_0 ** (2./3.) / sqrt(So)

                errorY = abs(area_n - area_0) / area_n
                area_0 = area_n
                area_c = (dsc * dsc * width_0 / grav) ** (1./3.)
                width_0  = r_interpol(areaTable,topwTable,nel,area_c)
                !fro=abs(dsc)/sqrt(grav*area_c**3.0/width_0)
                !print*, 'check point -1', area_0,  area_c, fro

            enddo
                !print*, 'nazmul j', j, oldY(1,j)

            y_norm = r_interpol(areaTable,elevTable,nel,area_0)
            y_crit = r_interpol(areaTable,elevTable,nel,area_c)

            !print*, 'check point 0', So, S_0,dsc,width_0, y_norm-z(i), y_crit-z(i), area_n, area_norm, area_c, fro
            !print*, 'check point -1', dsc, y_norm-z(1),  y_crit-z(1)
            !pause 500
    end subroutine normal_crit_y

 end module


text:
C     Last change:  ERB   6 Sep 2002    3:09 pm
C=======================================================================
      SUBROUTINE SEN1HFB6FM(BOTM,DELC,DELR,HNEW,HFB,IIPP,MXACTFB,MXHFB,
     &                      NBOTM,NCOL,NLAY,NROW,RHS,LAYHDT,NHFBNP,
     &                      IBOUND)
C-----VERSION 20010222 ERB
C     ******************************************************************
C     FOR PARAMETERIZED HORIZONTAL FLOW BARRIERS: FORMULATE RIGHT-HAND
C     SIDE TO CALCULATE SENSITIVITIES FOR HFB PARAMETERS
C     ******************************************************************
C        SPECIFICATIONS:
C     ------------------------------------------------------------------
      REAL BOTM, C, DE, DELC, DELR, DU, DV, HCDW, HFB, HO, RHS, TDW,
     &     U, V, ZERO
      INTEGER IBOUND, I1, I2, II, J1, J2, K, NBOTM, NCOL, NLAY,
     &        NROW
      DOUBLE PRECISION HNEW(NCOL,NROW,NLAY)
      DIMENSION BOTM(NCOL,NROW,0:NBOTM), DELC(NROW), DELR(NCOL),
     &          HFB(7,MXHFB), RHS(NCOL,NROW,NLAY), LAYHDT(NLAY),
     &          IBOUND(NCOL,NROW,NLAY)
      COMMON /DISCOM/LBOTM(999),LAYCBD(999)
      INCLUDE 'param.inc'
C     ------------------------------------------------------------------
C
      ZERO = 0.0
C
C     COPY UNMODIFIED CONDUCTANCE FROM ACTIVE PART OF HFB
C     ARRAY TO PART WHERE PARAMETER INFORMATION IS STORED
      CALL SSEN1HFB6FC(HFB,IIPP,MXACTFB,MXHFB,NHFBNP)
C
C     LOOP THROUGH BARRIERS CONTROLLED BY THIS PARAMETER
      DO 50 II = IPLOC(1,IIPP),IPLOC(2,IIPP)
        K = HFB(1,II)
        I1 = HFB(2,II)
        J1 = HFB(3,II)
        I2 = HFB(4,II)
        J2 = HFB(5,II)
C
C-------IF CELL ON EITHER SIDE OF BARRIER IS INACTIVE, CONTRIBUTION TO
C       SENSITIVITY IS ZERO -- ERB 8/25/03
        IF (IBOUND(J1,I1,K).EQ.0 .OR. IBOUND(J2,I2,K).EQ.0) GOTO 50
C
C-------GET THE HYDRAULIC CHARACTERISTIC FOR THIS BARRIER
C-------AND CALCULATE THO, TH1, HO, AND HP
        HCDW = HFB(6,II)*B(IIPP)
        TH0 = BOTM(J1,I1,LBOTM(K)-1) - BOTM(J1,I1,LBOTM(K))
        TH1 = BOTM(J2,I2,LBOTM(K)-1) - BOTM(J2,I2,LBOTM(K))
        HO = HNEW(J1,I1,K)   ! Changed 9/6/02 - ERB
        HP = HNEW(J2,I2,K)   ! Changed 9/6/02 - ERB
C-------FOR VARIABLE THICKNESS LAYERS, CALCULATE SATURATED THICKNESSES
        IF (LAYHDT(K).GT.0) THEN
          IF (HO.LT.BOTM(J1,I1,LBOTM(K)-1)) TH0 =
     &                                      HO - BOTM(J1,I1,LBOTM(K))
          IF (HP.LT.BOTM(J2,I2,LBOTM(K)-1)) TH1 =                       ! Changed 9/6/02 - ERB
     &                                      HP - BOTM(J2,I2,LBOTM(K))   ! Changed 9/6/02 - ERB
        ENDIF
        TDW = HCDW*(TH0+TH1)/2.0
C-------VARIABLES THAT DEPEND ON WHETHER THE BARRIER IS
C-------BETWEEN ROWS OR COLUMNS
c        IF (I1.EQ.I2) THEN
        IF (J1.EQ.J2) THEN   ! Changed 9/6/02 - ERB
          DE = DELR(J1)
        ELSE
          DE = DELC(I1)
        ENDIF
        HH = HO - HP
C-------U AND V ARE THE NUMERATOR AND DENOMINATOR, RESPECTIVELY, OF THE
C----------CONDUCTANCE INCLUDING THE HORIZONTAL FLOW BARRIER TERM.
C-------C IS THE CONDUCTANCE WITHOUT THE HORIZONTAL FLOW BARRIER TERM.
C-------DU IS THE DERIVATIVE OF THE NUMERATOR WITH RESPECT TO THE
C----------HFB PARAMETER
C-------DV IS THE DERIVATIVE OF THE DENOMINATOR WITH RESPECT TO THE
C----------HFB PARAMETER
C-------USE THE BASIC EQUATION FOR THE DERIVATIVE OF U/V
C----------THAT IS,(V DU - U DV)/V**2
        C = HFB(7,II)
        U = C*TDW*DE
        V = C+TDW*DE
        DU = C*DE*HFB(6,II)*(TH0+TH1)/2.0
        DV = DE*HFB(6,II)*(TH0+TH1)/2.0
        CO = (V*DU - U*DV)/V**2
        IF (CO.NE.ZERO) THEN
C---------CALCULATE RHS
          RHS(J2,I2,K) = RHS(J2,I2,K) - CO*HH
          RHS(J1,I1,K) = RHS(J1,I1,K) + CO*HH
        ENDIF
50    CONTINUE
      RETURN
      END
C=======================================================================
      SUBROUTINE SSEN1HFB6MD(C,CHAR,CO,DELC,DELR,HFB,I,J,K,MXACTFB,
     &                       NCOL,NHFB,NROW,TH0,TH1)
C-----VERSION 19990507 ERB
C     ******************************************************************
C     FOR HORIZONTAL FLOW BARRIERS: MODIFY CONDUCTANCE DERIVATIVES FROM
C     SSENLPF1CH (OR SSENHUF1CH) AS NEEDED.
C     ******************************************************************
C        SPECIFICATIONS:
C     ------------------------------------------------------------------
      REAL CO, DE, DELC, DELR, DU, DV, HFB, TDW, TH0, TH1, U, V
      INTEGER I, I1, I2, J, J1, J2, K, NCOL, NHFB, NROW
      CHARACTER*2 CHAR
      DIMENSION HFB(7,MXACTFB), DELR(NCOL), DELC(NROW)
C     ------------------------------------------------------------------
C
C-----DETERMINE IF THE CONDUCTANCE INVOLVES A HORIZONTAL FLOW BARRIER
C
C-----CHECK THE ROW AND COLUMN BY LOOPING THROUGH ALL HORIZONTAL FLOW
C-----BOUNDARY CELLS
C
      DO 40 NB = 1,NHFB
        K1 = HFB(1,NB)
        IF (K1.EQ.K) THEN
          I1 = HFB(2,NB)
          J1 = HFB(3,NB)
          I2 = HFB(4,NB)
          J2 = HFB(5,NB)
          IF (I1.EQ.I.AND.J1.EQ.J .OR. I2.EQ.I.AND.J2.EQ.J) THEN
            IF (I1.EQ.I.AND.J1.EQ.J .AND.
     &         (I2.EQ.I+1 .OR. J2.EQ.J+1)) THEN
              IF (I2.EQ.I+1 .AND. CHAR.EQ.'CC') GOTO 50
              IF (J2.EQ.J+1 .AND. CHAR.EQ.'CR') GOTO 50
            ELSEIF (I2.EQ.I.AND.J2.EQ.J .AND.
     &             (I1.EQ.I+1. .OR. J1.EQ.J+1)) THEN
              IF (I1.EQ.I+1 .AND. CHAR.EQ.'CC') GOTO 50
              IF (J1.EQ.J+1 .AND. CHAR.EQ.'CR') GOTO 50
            ENDIF
          ENDIF
        ENDIF
40    CONTINUE
      RETURN
C-----THE CONDUCTANCE INVOLVES A HORIZONTAL FLOW BARRIER
50    CONTINUE
      TDW = HFB(6,NB)*(TH0+TH1)/2.0
      IF (I1.EQ.I2) THEN
        DE = DELC(I1)
      ELSE
        DE = DELR(J1)
      ENDIF
C-----U AND V ARE THE NUMERATOR AND DENOMINATOR, RESPECTIVELY, OF THE
C--------CONDUCTANCE INCLUDING THE HORIZONTAL FLOW BARRIER TERMS.
C-----C IS THE CONDUCTANCE WITHOUT THE HORIZONTAL FLOW BARRIER TERMS.
C-----DU IS THE DERIVATIVE OF THE NUMERATOR WITH RESPECT TO THE
C--------T PARAMETER
C-----DV IS THE DERIVATIVE OF THE DENOMINATOR WITH RESPECT TO THE
C--------T PARAMETER
C-----USE THE BASIC EQUATION FOR THE DERIVATIVE OF U/V
C--------THAT IS,(V DU - U DV)/V**2
      U = C*TDW*DE
      V = C+TDW*DE
      DU = TDW*DE*CO
      DV = CO
      CO = (V*DU - U*DV)/V**2
      RETURN
      END
C=======================================================================
      SUBROUTINE SSEN1HFB6FC(HFB,IIPP,MXACTFB,MXHFB,NHFBNP)
C
C-----VERSION 19990524 ERB
C     ******************************************************************
C     COPY UNMODIFIED CONDUCTANCE FOR PARAMETERIZED HORIZONTAL FLOW
C     BARRIERS FROM ACTIVE PART OF HFB ARRAY TO PARAMETER PART
C     ******************************************************************
C        SPECIFICATIONS:
C     ------------------------------------------------------------------
      REAL HFB(7,MXHFB)
      INCLUDE 'param.inc'
C     ------------------------------------------------------------------
      IFIRST = NHFBNP + 1
C     LOOP THROUGH BARRIERS CONTROLLED BY PARAMETER IIPP
      DO 50 IBP = IPLOC(1,IIPP),IPLOC(2,IIPP)
        IBPK = HFB(1,IBP)
        IBPI1 = HFB(2,IBP)
        IBPJ1 = HFB(3,IBP)
        IBPI2 = HFB(4,IBP)
        IBPJ2 = HFB(5,IBP)
C       SEARCH FOR MATCHING BARRIER IN ACTIVE PART OF HFB ARRAY
        DO 30 IBA = IFIRST,MXACTFB
          IBAK = HFB(1,IBA)
          IBAI1 = HFB(2,IBA)
          IBAJ1 = HFB(3,IBA)
          IBAI2 = HFB(4,IBA)
          IBAJ2 = HFB(5,IBA)
          IF (IBPK.EQ.IBAK .AND. IBPI1.EQ.IBAI1 .AND. IBPJ1.EQ.IBAJ1
     &        .AND. IBPI2.EQ.IBAI2 .AND. IBPJ2.EQ.IBAJ2) THEN
C           COPY UNMODIFIED CONDUCTANCE
            HFB(7,IBP) = HFB(7,IBA)
            GOTO 40
          ENDIF
   30   CONTINUE
   40   CONTINUE
   50 CONTINUE
C
      RETURN
      END

max_stars_repo_name: quanpands/wflow

max_issues_repo_name: quanpands/wflow

max_forks_repo_name: quanpands/wflow


text:
module type_function_info
    use mod_kinds,  only: ik
    use type_seed,  only: seed_t

    !> Container for flux information.
    !!
    !!
    !!  @author Nathan A. Wukie
    !!  @date   2/1/2016
    !!
    !!
    !---------------------------------------------------------------------
    type, public :: function_info_t

        integer(ik)     :: type     !< Function type. ex: 'boundary_advective_flux', 'boundary_diffusive_flux', etc.
        integer(ik)     :: ifcn     !< Index of the function of the give type being computed.
        integer(ik)     :: idepend  !< Dependency index of a given element to the function. Chimera could be > 1.
        integer(ik)     :: idiff    !< Index of the direction being linearized. XI_MIN, ETA_MAX, DIAG, etc.
        type(seed_t)    :: seed     !< Indices of the element being linearized

    end type function_info_t
    !*********************************************************************

end module type_function_info

max_stars_repo_name: wanglican/ChiDG

max_issues_repo_name: haohb/ChiDG

max_forks_repo_name: haohb/ChiDG


text:
*
*     ..........................................................   :
*                                                                  :
*          Block                                                   :
*                 S U B M A T R I X    E L E M E N T S             :
*                                                                  :
*     For Calculation Angular Momentum Coefficients for            :
*     Relativistic operators                                       :
*                                                                  :
*     Written by G. Gaigalas,                                      :
*                  Department  of  Computer Science,               :
*                  Vanderbilt University,  Nashville               :
*                                                   October 1996   :
*                                                                  :
*     ..........................................................   :
*        i)  Spin - orbit                                          :
*     ..........................................................   :
*
*     -------------------------------------------------------------
*      S P I N O R
*     -------------------------------------------------------------
*
*     CALCULATE ONE-ELECTRON NUCLEAR SPIN-ORBIT OPERATOR           *
*                                                                  *
*     Written by G. Gaigalas,                                      *
*     Universite Libre de Bruxelles, Belgium         October 1995  *
*
      SUBROUTINE SPINOR(L1,L2,I,A)
      IMPLICIT DOUBLEPRECISION (A-H,O-Z)
      COMMON/CONSTS/ZERO,TENTH,HALF,ONE,TWO,THREE,FOUR,SEVEN,ELEVEN,EPS
      COMMON/MEDEFN/IHSH,NJ(16),LJ(16),NOSH1(16),NOSH2(16),J1QN1(31,3),
     :     J1QN2(31,3),IJFUL(16)
      A=ZERO
      I=5
      IF(L1.NE.L2)RETURN
      IF(L1.EQ.0)RETURN
      A=DSQRT(DBLE(6*L1*(L1+1)*(2*L1+1)))
      KK=IHSH+IHSH-1
      A=-A*DSQRT(DBLE(J1QN1(KK,2)*J1QN1(KK,3)))
      RETURN
      END


text:
          SUBROUTINE get_grid_dims(nproc,nprow,npcol)
! ...     This subroutine factorizes the number of processors (NPROC)
! ...     into NPROW and NPCOL,  that are the sizes of the 2D processors mesh.
            integer nproc,nprow,npcol
            integer sqrtnp,i
            sqrtnp = INT( SQRT( DBLE(nproc) ) + 1 )
            DO i=1,sqrtnp
              IF(MOD(nproc,i).EQ.0) nprow = i
            END DO
            npcol = nproc/nprow
            RETURN
          END SUBROUTINE

max_stars_repo_name: sscerr/eTF

max_issues_repo_name: sscerr/eTF

max_forks_repo_name: sscerr/eTF


text:
!====================================================================
!-----------------------------------------------------------------------
!       NASA/GSFC, Data Assimilation Office, Code 910.3, GEOS/DAS      !
!-----------------------------------------------------------------------
!BOP
! !MODULE: m_RadDataF -- Utilities for extracting raob data from data files
!
! !DESCRIPTION:
!
! !INTERFACE:
!
      module m_RadDataF
      use m_RadData, only : radcor_profiles
!ods      use m_RadODS,  only : ods_info
      use m_RadBufr, only : prep_info, prep_options
      implicit    NONE
      private   ! except

      public ::
     .   Inq_Raobs,          ! Inquire routine for obtaining obs
     .   Get_Raobs,          ! Get raob data
     .   Put_Raobs,          ! Put raob data to buffer file
     .   Clean_Raobs,        ! Clean up data structures
     .   Set_Options,        ! Set the options in file_options.
     .   file_info,          ! File information
     .   file_options,       ! Options for reading/writing to file
     .   radcor_profiles     ! Radcor data structure

      interface Clean_Raobs
         module procedure
     .      Clean_Raobs_
      end interface
      interface Inq_Raobs
         module procedure
     .      Inq_Raobs_
      end interface
      interface Get_Raobs
         module procedure
     .      Get_Raobs_
      end interface
      interface Put_Raobs
         module procedure
     .      Put_Raobs_
      end interface
      interface Set_Options
         module procedure
     .       Init_Options,
     .       Copy_Options
      end interface
!
! !REVISION HISTORY:
!     16Apr2003  C. Redder  Original code
!     13Nov2003  C. Redder  Replaced the data structure radcor with
!                           radcor_profiles.
!     20Apr2004  C. Redder  Added generic interface, SetOptions, and
!                           data structure, file_options
!EOP
!-----------------------------------------------------------------

      character (len=*), parameter :: MyModule = 'm_RadDataF'

      type file_info
!ods         type (  ods_info    ) :: ODS
         type ( prep_info    ) :: Prep
      end type

      type file_options
         type ( prep_options ) :: Prep
      end type

      contains
!....................................................................

!-------------------------------------------------------------------------
!         Nasa/GSFC, Data Assimilation Office, Code 910.3, GEOS/DAS      !
!-------------------------------------------------------------------------
! BOP
!
! !ROUTINE: Init_Options () -- initialize option flags for reading/writing to data files
!
! !DESCRIPTION:
!     The routine initializes the contents of the input data structure.
!
! !INTERFACE:
      subroutine Init_Options ( this )
      use m_RadBufr, only : Set_PrepOpt => Set_Options
      implicit   NONE
!
! !OUTPUT PARAMETERS:
      type ( file_options ), intent   (out), optional ::
     .   this              ! Output data structure.
!
! !SEE ALSO:
!
! !REVISION HISTORY:
!     20Apr2004  C. Redder   Origional code
!
! EOP
!-------------------------------------------------------------------------

      call Set_PrepOpt ( this % Prep )

      return
      end subroutine Init_Options

!....................................................................

!-------------------------------------------------------------------------
!         Nasa/GSFC, Data Assimilation Office, Code 910.3, GEOS/DAS      !
!-------------------------------------------------------------------------
! BOP
!
! !ROUTINE: Copy_Options () --- Copy option flags for reading/writing to data files
!
! !DESCRIPTION:
!     The routine copies the contents from the input data structure to
!     the output data structure if it is present
!
! !INTERFACE:
      subroutine Copy_Options ( this_in, this_out )
      use m_RadBufr, only : Set_PrepOpt => Set_Options
!
! !INPUT PARAMETERS:
      implicit   NONE
      type ( file_options ), intent (in)  ::
     .   this_in
!
! !OUTPUT PARAMETERS:
      type ( file_options ), intent (out) ::
     .   this_out            ! Output data structure.
!
! !SEE ALSO:
!
! !REVISION HISTORY:
!     20Apr2004  C. Redder   Origional code
!
! EOP
!-------------------------------------------------------------------------

      call Set_PrepOpt ( this_in % Prep, this_out % Prep )

      return
      end subroutine Copy_Options

!.................................................................

!-------------------------------------------------------------------------
!         Nasa/GSFC, Data Assimilation Office, Code 910.3, GEOS/DAS      !
!-------------------------------------------------------------------------
! BOP
!
! !ROUTINE: Inq_Raobs_ () --- Inquire data from data file (ODS or Prep)
!
! !DESCRIPTION:
!
! !INTERFACE:
      subroutine Inq_Raobs_ ( File, Type, stat, NSyn, DTBeg, DTEnd )
!ods      use m_RadODS,   only : ODS_Inq,  Get_ODS  => Get_Raobs
      use m_RadBufr,  only : Bufr_Inq, Get_Prep => Get_Raobs
      use m_AdvError, only : ErrStat,  Alloc, WPErr, ItoDate
      use m_TextUtil, only : LowerCase
!
! !INPUT PARAMETERS:
      implicit   NONE
      character (len=*),  intent (in)  ::
     .   File,            ! Data filename
     .   Type             ! File type (i.e. 'ODS' or 'Prep',
     .                    !   case insensitive)
!
! !OUTPUT PARAMETERS:
      integer, optional,  intent (out) ::
     .   NSyn,            ! Number of synoptic times per day
     .   DTBeg, DTEnd     ! Date/time tag (in YYYYMMHHHH format) for the
     .                    ! ... latest synoptic Julian hour on file
      integer,            intent (out) ::
     .   stat             ! Returned status code
!
! !SEE ALSO:
!
! !REVISION HISTORY:
!     03Jul2003  C. Redder   Original code
!     19Feb2004  C. Redder   Made file type case insensitive
! EOP
!-------------------------------------------------------------------------
      character (len=*), parameter ::
     .   MyName  = MyModule // '::Inq_Raobs_'

      integer :: DTBeg_, DTEnd_, NSyn_
      character ( len = 500 ) :: Tail
      character ( len =  80 ) :: Prep_Type, ODS_Type = "!?*", Type_

      Tail = '\n   file  = \C'   // File

!     Get parameters
!     --------------
      call Bufr_Inq ( Prep_Type = Prep_Type )
!ods      call ODS_Inq  ( ODS_Type  = ODS_Type  )
      Prep_Type = LowerCase ( Prep_Type )
!ods      ODS_Type  = LowerCase (  ODS_Type )

!     Get data from the ODS file
!     --------------------------
      Type_ = LowerCase ( Type )
      if      ( Type_ .eq. ODS_Type  ) then
!ods         call ODS_Inq  ( File, stat, NSyn  = NSyn_,
!ods     .                               DTBeg = DTBeg_,
!ods     .                               DTEnd = DTEnd_ )
         if ( stat .ne. 0 ) then
            call WPErr ( MyName, ErrStat ( 'ODS_Inq',  stat ))
            return
         end if

!     ... or "Prep" data file
!     -----------------------
      else if ( Type_ .eq. Prep_Type ) then
         call Bufr_Inq ( File, stat, SynTag = DTBeg_ )
         if ( stat .ne. 0 ) then
            call WPErr ( MyName, ErrStat ( 'Bufr_Inq', stat ))
            return
         end if
         DTEnd_ = DTBeg_
         NSyn_  = 1

!     ... and other type results in an error
!    ---------------------------------------
      else
         call WPErr ( MyName, 'Invalid file type, ' // trim ( Type )
     .                      // Tail )
         stat = 1
         return
      end if

      if ( present ( NSyn  ))  NSyn  = NSyn_
      if ( present ( DTBeg ))  DTBeg = DTBeg_
      if ( present ( DTEnd ))  DTEnd = DTEnd_

      return
      end subroutine Inq_Raobs_
!.................................................................

!-------------------------------------------------------------------------
!         Nasa/GSFC, Data Assimilation Office, Code 910.3, GEOS/DAS      !
!-------------------------------------------------------------------------
! BOP
!
! !ROUTINE: Clean_Raobs_ () --- Clean up data structure (ODS or Prep)
!
! !DESCRIPTION:
!
! !INTERFACE:
      subroutine Clean_Raobs_ ( FInfo, Obs, stat )
!ods      use m_RadODS,   only : ODS_Inq,  Clean_ODS  => Clean_Raobs
      use m_RadBufr,  only : Bufr_Inq, Clean_Prep => Clean_Raobs
      use m_RadData,  only : radcor_profiles
      use m_AdvError, only : ErrStat,  Alloc, WPErr, ItoDate
      use m_TextUtil, only : LowerCase
!
! !INPUT/OUTPUT PARAMETERS:
      type ( file_info ),       intent (inout) ::
     .   FInfo                  ! Information about the data file which is
                                !   used to write the data to the input file
      type ( radcor_profiles ), intent (inout) ::
     .   Obs                    ! Raob data structure
! !OUTPUT PARAMETERS:
      integer,                  intent (out)   ::
     .   stat                   ! Returned status code
!
! !SEE ALSO:
!
! !REVISION HISTORY:
!     03Jul2003  C. Redder   Original code
!     13Nov2003  C. Redder   Replaced the data structure radcor with
!                            radcor_profiles.
!     19Feb2004  C. Redder   Made file type case insensitive
! EOP
!-------------------------------------------------------------------------
      character (len=*), parameter ::
     .   MyName  = MyModule // '::Clean_Raobs_'
      character ( len = 500 ) :: Tail
      character ( len =  80 ) :: Prep_Type, ODS_Type = "!?*", Type

      Tail = '\n   Input file  = \C'    // trim ( Obs % ObsFile )

!     Get parameters
!     --------------
      call Bufr_Inq ( Prep_Type = Prep_Type )
!ods      call ODS_Inq  ( ODS_Type  = ODS_Type  )
      Prep_Type = LowerCase ( Prep_Type )
!ods      ODS_Type  = LowerCase (  ODS_Type )

!     Get data from the ODS file
!     --------------------------
      Type = LowerCase ( Obs % FileType )
      if      ( Type .eq. ODS_Type  ) then
!ods         call Clean_ODS  ( FInfo % ODS, Obs, stat )
         if ( stat .ne. 0 ) then
            call WPErr ( MyName, ErrStat ( 'Clean_ODS',  stat ))
            return
         end if

!     ... or "Prep" data file
!     -----------------------
      else if ( Type .eq. Prep_Type ) then
         call Clean_Prep ( FInfo % Prep, Obs, stat )
         if ( stat .ne. 0 ) then
            call WPErr ( MyName, ErrStat ( 'Clean_Prep', stat ))
            return
         end if

!     ... and other type results in an error
!    ---------------------------------------
      else
         call WPErr ( MyName, 'Invalid file type, ' // trim ( Type )
     .                      // Tail )
         stat = 1
         return
      end if

      return
      end subroutine Clean_Raobs_
!.................................................................

!-------------------------------------------------------------------------
!         Nasa/GSFC, Data Assimilation Office, Code 910.3, GEOS/DAS      !
!-------------------------------------------------------------------------
! BOP
!
! !ROUTINE: Get_Raobs_ () --- Read data from input file (ODS or Prep)
!
! !DESCRIPTION:
!
! !INTERFACE:
      subroutine Get_Raobs_ ( File, Type, NYMDH, FInfo, Obs, stat,
     .                        Options )
!ods      use m_RadODS,    only : ODS_Inq,  Get_ODS  => Get_Raobs
      use m_RadBufr,   only : Bufr_Inq, Get_Prep => Get_Raobs
      use m_RadData,   only : radcor_profiles
      use m_AdvError,  only : ErrStat,  Alloc, WPErr, ItoDate
      use m_TextUtil,  only : LowerCase
!
! !INPUT PARAMETERS:
      implicit   NONE
      character (len=*),        intent (in)    ::
     .   File,                  ! Input filename
     .   Type                   ! File type (i.e. 'ODS' or 'Prep',
     .                          !   case insensitive)
      integer,                  intent (in)    ::
     .   NYMDH                  ! Synoptic date and time for the desired
                                !   data, in YYYYMMDDHH format.
      type ( file_options),     intent (in), optional ::
     .   Options                ! Options for reading/writing.
!
! !OUTPUT PARAMETERS:
      type ( file_info ),       intent (inout) ::
     .   FInfo                  ! Information about the data file which
                                !   is used to write the data to the
                                !   input file
      type ( radcor_profiles ), intent (inout) ::
     .   Obs                    ! Raob data structure
      integer,                  intent (out)   ::
     .   stat                   ! Returned status code
!
! !SEE ALSO:
!
! !REVISION HISTORY:
!     03Jul2003  C. Redder   Original code
!     13Nov2003  C. Redder   Replaced the data structure radcor with
!                            radcor_profiles.
!     19Feb2004  C. Redder   Made file type case insensitive
!     04Mar2004  C. Redder   Added optional arguments, PName and
!                            updated_QC.
!     01Apr2004  C. Redder   Removed the optional arguments, PName and
!                            updated_QM and added the optional argument,
!                            Options.
!     01Jun2004  C. Redder   Fixed bug in implementing options
!
! EOP
!-------------------------------------------------------------------------
      character (len=*), parameter ::
     .   MyName  = MyModule // '::Get_Raobs_'
      type ( file_options )   :: Opt
      integer :: NYMD, NHMS
      character ( len = 500 ) :: Tail
      character ( len =  80 ) :: Prep_Type, ODS_Type = "!?*", Type_

      NYMD = NYMDH / 100
      NHMS = mod ( NYMDH, 100 ) * 10000

      Tail = '\n   Date/time = \C' // ItoDate ( NYMD, NHMS )
     .    // '\n   File      = '   // File

!     Get parameters
!     --------------
      call Bufr_Inq ( Prep_Type = Prep_Type )
!ods      call ODS_Inq  ( ODS_Type  = ODS_Type  )
      Prep_Type = LowerCase ( Prep_Type )
!ods      ODS_Type  = LowerCase (  ODS_Type )

!     Implement options
!     -----------------
      call Init_Options ( Opt )
      if ( present ( Options )) call Copy_Options ( Options, Opt )

!     Get data from the ODS file
!     --------------------------
      Type_     = LowerCase ( Type )
      if      ( Type_ .eq. ODS_Type  ) then
!ods         call Get_ODS  ( File, NYMD, NHMS, FInfo % ODS,  Obs, stat )
         if ( stat .ne. 0 ) then
            call WPErr ( MyName, ErrStat ( 'Get_ODS',  stat ))
            return
         end if

!     ... or "Prep" data file
!     -----------------------
      else if ( Type_ .eq. Prep_Type ) then
         call Get_Prep ( File, NYMD, NHMS, FInfo % Prep, Obs, stat,
     .                   Options = Opt % Prep )
c         call Get_Prep ( File, NYMDH, FInfo % Prep, Obs, stat )
         if ( stat .ne. 0 ) then
            call WPErr ( MyName, ErrStat ( 'Get_Prep', stat ))
            return
         end if

!     ... and other type results in an error
!    ---------------------------------------
      else
         call WPErr ( MyName, 'Invalid file type, ' // trim ( Type )
     .                      // Tail )
         stat = 1
         return
      end if

      return
      end subroutine Get_Raobs_
!.................................................................

!-------------------------------------------------------------------------
!         Nasa/GSFC, Data Assimilation Office, Code 910.3, GEOS/DAS      !
!-------------------------------------------------------------------------
! BOP
!
! !ROUTINE: Put_Raobs_ () --- Write data to output file (ODS or Prep)
!
! !DESCRIPTION:
!
! !INTERFACE:
      subroutine Put_Raobs_ ( Tplate, NYMDH, FInfo, Obs, stat,
     .                        ExpID,  cleanup, Options )
!ods      use m_RadODS,    only : ODS_Inq,  Put_ODS  => Put_Raobs
      use m_RadBufr,   only : Bufr_Inq, Put_Prep => Put_Raobs
      use m_RadData,   only : radcor_profiles
      use m_AdvError,  only : ErrStat,  Alloc, WPErr, ItoDate
      use m_RadTplate, only : StrTemplate
      use m_TextUtil,  only : LowerCase
!
! !INPUT PARAMETERS:
      implicit   NONE
      integer,                  intent (in)    ::
     .   NYMDH                  ! Synoptic date and time for the
!                               !    desired data, in YYYYMMDDHH format.
      character ( len = * ),    intent (in), optional ::
     .   ExpID
      logical, optional,        intent (in)    ::
     .   cleanup                ! = .true. to deallocate and reset
                                !    input data structure.  Default:
                                !    cleanup = .false.
      type ( file_options ),    intent (in), optional ::
     .   Options                ! Other file options for writing
! !INPUT/OUTPUT PARAMETERS:
      character (len=*),        intent (inout) ::
     .   Tplate                 ! On input,  template of output filename.
                                ! On output, generated output filename
      type ( file_info ),       intent (inout) ::
     .   FInfo                  ! Information about the data file which
                                !   is used to write the data to the
                                !   output file
      type ( radcor_profiles ), intent (inout) ::
     .   Obs                    ! Raob data structure
                                !   On output, both data structures are
                                !   reset if the optional arguent,
!                               !   cleanup = .true.
! !OUTPUT PARAMETERS:
      integer,                  intent (out)   ::
     .   stat                   ! Returned status code
!
! !SEE ALSO:
!
! !REVISION HISTORY:
!     14Jul2003  C. Redder   Original code
!     13Nov2003  C. Redder   Replaced the data structure radcor with
!                            radcor_profiles.
!     19Feb2004  C. Redder   Made file type case insensitive
!     01Jun2004  C. Redder   Fixed bug in implementing options
!
! EOP
!-------------------------------------------------------------------------

      character (len=*), parameter ::
     .   MyName  = MyModule // '::Put_Raobs_'
      type ( file_options ) :: Opt
      integer :: NYMD, NHMS, iBeg, Len, iEnd
      logical :: set_id
      character ( len = 255 ) :: InFile, OutFile, ExpID_
      character ( len = 500 ) :: Tail
      character ( len =  80 ) :: Prep_Type, ODS_Type = "!?*", Type

      NYMD   = NYMDH / 100
      NHMS   = mod ( NYMDH, 100 ) * 10000
      InFile  = Obs % ObsFile
      Tail   = '\n   Date/time     = \C'  // ItoDate ( NYMD, NHMS )
     .      // '\n   Input file    = '    // trim ( InFile )
     .      // '\n   File template = '    // trim ( Tplate )

!     Implement options
!     -----------------
      call Init_Options ( Opt )
      if ( present ( Options )) call Copy_Options ( Options, Opt )

!     Get the ID tag
!     --------------
      if ( present ( ExpID )) then
         ExpID_ = ExpID
      else
         iBeg   = scan ( InFile, '/', back = .true. ) + 1
          Len   = scan ( InFile ( iBeg : ), '.' ) - 1
         if ( Len .eq. 0 )
     .      Len = len_trim ( InFile ( iBeg : ) )
         iEnd   = iBeg + Len - 1
         ExpID_ = InFile ( iBeg : iEnd )
      end if

!     Generate file name from template
!     --------------------------------
      call StrTemplate  ( OutFile, TPlate,
     .                    xid  = trim ( ExpID_ ),
     .                    nymd = NYMD,
     .                    nhms = NHMS, stat = stat )
      if ( stat .ne. 0 ) then
          call WPErr ( MyName, ErrStat ( 'StrTemplate', stat )
     .                      // Tail )
          return
      end if

      Tail   = '\n   Date/time   = \C'  // ItoDate ( NYMD, NHMS )
     .      // '\n   Input file  = '    // trim ( InFile  )
     .      // '\n   Output file = '    // trim ( OutFile )

!     Get parameters
!     --------------
      call Bufr_Inq ( Prep_Type = Prep_Type )
!ods      call ODS_Inq  ( ODS_Type  = ODS_Type  )
      Prep_Type = LowerCase ( Prep_Type )
!ods      ODS_Type  = LowerCase (  ODS_Type )

!     Write data to the ODS file
!     --------------------------
      Type = LowerCase ( Obs % FileType )
      if      ( Type .eq. ODS_Type  ) then
!ods         call Put_ODS  ( OutFile, NYMD, NHMS, FInfo % ODS,  Obs, stat,
!ods     .                   cleanup = cleanup )
         if ( stat .ne. 0 ) then
            call WPErr ( MyName, ErrStat ( 'Put_ODS',  stat ))
            return
         end if

!     ... or "Prep" data file
!     -----------------------
      else if ( Type .eq. Prep_Type ) then
         call Put_Prep ( OutFile, NYMD, NHMS, FInfo % Prep, Obs, stat,
     .                   Options = Opt % Prep,
     .                   cleanup = cleanup )
         if ( stat .ne. 0 ) then
            call WPErr ( MyName, ErrStat ( 'Put_Prep', stat ))
            return
         end if

!     ... and other type results in an error
!    ---------------------------------------
      else
         call WPErr ( MyName, 'Invalid file type, ' // trim ( Type )
     .                      // Tail )
         stat = 1
         return
      end if

      stat   = 0
      Tplate = OutFile
      return
      end subroutine Put_Raobs_
!.................................................................
      end module m_RadDataF
!====================================================================

max_stars_repo_name: GEOS-ESM/GEOSadas

max_issues_repo_name: GEOS-ESM/GEOSadas

max_forks_repo_name: GEOS-ESM/GEOSadas


text: Sacramento Sac gal working at UCD


text:
MODULE problem_setup_module

  USE ic_module
  USE shared_data
  USE timer
  USE fields
  USE strings

  IMPLICIT NONE

  PRIVATE
  PUBLIC :: problem_setup

CONTAINS

  SUBROUTINE problem_setup(deck_state, problem)

    INTEGER, INTENT(IN) :: deck_state
    CHARACTER(LEN=c_max_string_length), INTENT(IN) :: problem

    IF (deck_state == c_ds_first) THEN
      CALL control_initialise
      CALL boundary_initialise(problem)
      CALL custom_problem_setup(deck_state, problem)
    ELSE
      CALL fields_initialise(problem)
      CALL custom_problem_setup(deck_state, problem)
      CALL laser_finalise
      CALL species_finalise
    ENDIF

  END SUBROUTINE problem_setup

  SUBROUTINE boundary_initialise(problem)

    CHARACTER(LEN=c_max_string_length), INTENT(IN) :: problem
    INTEGER :: i

    ! Could add problem specific boundary conditions here
    ! Default values
    DO i = c_bd_x_min, c_bd_z_max
      bc_field(i) = c_bc_periodic
      bc_particle(i) = c_bc_periodic
    ENDDO

  END SUBROUTINE boundary_initialise

  SUBROUTINE control_initialise

    ! Default values
    allow_cpu_reduce = .TRUE.
    timer_collect = .TRUE.
    nx_global = -1
    ny_global = -1
    nz_global = -1
    x_min =  c_largest_number
    x_max = -c_largest_number
    y_min =  c_largest_number
    y_max = -c_largest_number
    z_min =  c_largest_number
    z_max = -c_largest_number
    nprocx = 0
    nprocy = 0
    nprocz = 0
    npart_global = 0
    nsteps = -1
    t_end = c_largest_number
    dt_multiplier = 0.95_num
    dlb_threshold = 1.0_num
    use_balance = .FALSE.
    CALL set_field_order(2)
    stdout_frequency = 10
    use_random_seed = .FALSE.
    particle_push_start_time = 0.0_num
    n_species = 0
    fixed_fields = .FALSE.

  END SUBROUTINE control_initialise

  SUBROUTINE fields_initialise(problem)

    CHARACTER(LEN=c_max_string_length), INTENT(IN) :: problem
    INTEGER :: i,j,k
    REAL(num) :: xp
    REAL(num), PARAMETER :: em_wave_mag = 1e6_num

    ! Could include problem specific fields here.
    ! Default to zero
    ex = 0.0_num
    ey = 0.0_num
    ez = 0.0_num
    bx = 0.0_num
    by = 0.0_num
    bz = 0.0_num
    SELECT CASE (TRIM(problem))

    CASE ('one_stream')
       IF (.false.) THEN
          bx = 1.0_num
          ey = 1.0_num

          !To do this right, need to get the offsets to the grid positions.
          do i=1-ng,nx+ng
             do j=1-ng,ny+ng
                do k=1-ng,nz+ng
                   !Correct offset for bz, is there a variable
                   ! storing offset per field/direction?
                   xp = x_grid_min_local+(i-0.5_num)*dx
                   bz(i,j,k) = cos(xp)
                end do
             end do
          end do
       END IF

    CASE ('drift_kin_default')
       ex = 0.0_num
       ey = 0.0_num
       ez = 0.0_num
       bx = 0.0_num
       by = 0.0_num
       bz = 1.0_num

    CASE ('em_wave')
      ex = 0.0_num
      ez = 0.0_num
      bx = 0.0_num
      by = 0.0_num

      DO i = 1-ng, nx+ng
        xp = x(i)
        ey(i,:,:) = em_wave_mag * COS(xp)
        xp = xp + 0.5_num * dx
        bz(i,:,:) = em_wave_mag * COS(xp) / c
      END DO

    CASE default

    END SELECT

  END SUBROUTINE fields_initialise

  SUBROUTINE laser_finalise

    INTEGER :: errcode, ierr
    TYPE(laser_block), POINTER :: current
    INTEGER :: error, io, iu

    errcode = c_err_none

    error = 0
    current => laser_x_min
    DO WHILE(ASSOCIATED(current))
      IF (current%omega < 0.0_num) error = IOR(error, 1)
      IF (current%amp < 0.0_num) error = IOR(error, 2)
      current => current%next
    ENDDO

    current => laser_x_max
    DO WHILE(ASSOCIATED(current))
      IF (current%omega < 0.0_num) error = IOR(error, 1)
      IF (current%amp < 0.0_num) error = IOR(error, 2)
      current => current%next
    ENDDO

    current => laser_y_min
    DO WHILE(ASSOCIATED(current))
      IF (current%omega < 0.0_num) error = IOR(error, 1)
      IF (current%amp < 0.0_num) error = IOR(error, 2)
      current => current%next
    ENDDO

    current => laser_y_max
    DO WHILE(ASSOCIATED(current))
      IF (current%omega < 0.0_num) error = IOR(error, 1)
      IF (current%amp < 0.0_num) error = IOR(error, 2)
      current => current%next
    ENDDO

    current => laser_z_min
    DO WHILE(ASSOCIATED(current))
      IF (current%omega < 0.0_num) error = IOR(error, 1)
      IF (current%amp < 0.0_num) error = IOR(error, 2)
      current => current%next
    ENDDO

    current => laser_z_max
    DO WHILE(ASSOCIATED(current))
      IF (current%omega < 0.0_num) error = IOR(error, 1)
      IF (current%amp < 0.0_num) error = IOR(error, 2)
      current => current%next
    ENDDO

    IF (IAND(error,1) /= 0) THEN
      IF (rank == 0) THEN
        DO iu = 1, nio_units ! Print to stdout and to file
          io = io_units(iu)
          WRITE(io,*) '*** ERROR ***'
          WRITE(io,*) 'Must define a "lambda" or "omega" for every laser.'
        ENDDO
      ENDIF
      errcode = c_err_missing_elements
      CALL MPI_ABORT(MPI_COMM_WORLD, errcode, ierr)
    ENDIF

    IF (IAND(error,2) /= 0) THEN
      IF (rank == 0) THEN
        DO iu = 1, nio_units ! Print to stdout and to file
          io = io_units(iu)
          WRITE(io,*) '*** ERROR ***'
          WRITE(io,*) 'Must define an "amp" or "irradiance" for every laser.'
        ENDDO
      ENDIF
      errcode = c_err_missing_elements
      CALL MPI_ABORT(MPI_COMM_WORLD, errcode, ierr)
    ENDIF

  END SUBROUTINE laser_finalise

  SUBROUTINE species_finalise

    INTEGER :: i, io, iu
    CHARACTER(LEN=8) :: string
    INTEGER :: errcode, ierr
    TYPE(particle_species), POINTER :: species, next_species

    ! Sanity check
    next_species => species_list
    DO i = 1, n_species
      species => next_species
      next_species => species%next

      IF (rank == 0) THEN
        CALL integer_as_string(i, string)
        PRINT*, 'Name of species ', TRIM(ADJUSTL(string)), ' is ', &
            TRIM(species%name)
      ENDIF

      IF (species%mass < 0) THEN
        IF (rank == 0) THEN
          DO iu = 1, nio_units ! Print to stdout and to file
            io = io_units(iu)
            WRITE(io,*) '*** ERROR ***'
            WRITE(io,*) 'Particle species "', TRIM(species%name), &
                '" has negative mass'
          ENDDO
        ENDIF
        CALL MPI_ABORT(MPI_COMM_WORLD, errcode, ierr)
      ENDIF

      IF (.NOT. (species%charge < HUGE(1.0_num))) THEN
        IF (rank == 0) THEN
          DO iu = 1, nio_units ! Print to stdout and to file
            io = io_units(iu)
            WRITE(io,*) '*** ERROR ***'
            WRITE(io,*) 'No charge specified for particle species "', &
                TRIM(species%name), '"'
          ENDDO
        ENDIF
        CALL MPI_ABORT(MPI_COMM_WORLD, errcode, ierr)
      ENDIF

      IF (species%npart_per_cell >= 0) THEN
        IF (species%count >= 0 .AND. rank == 0) THEN
          DO iu = 1, nio_units ! Print to stdout and to file
            io = io_units(iu)
            WRITE(io,*) '*** WARNING ***'
            WRITE(io,*) 'Two forms of npart used for particle species "', &
                TRIM(species%name), '"'
            WRITE(io,*) 'Just using "npart_per_cell".'
          ENDDO
        ENDIF
        species%count = INT(species%npart_per_cell, i8)
      ENDIF
    ENDDO

  END SUBROUTINE species_finalise

END MODULE problem_setup_module

max_stars_repo_name: TomGoffrey/minepoch

max_issues_repo_name: TomGoffrey/minepoch

max_forks_repo_name: TomGoffrey/minepoch


text:
************************************************************************
*
*     Function GETINT              Called by: AMLLOAD, LADLOAD, MLELOAD,
*                                             SETXL2
*
*     convert all or part of a character variable to an integer.
*
*     This function seems trivial, but it allows one to access
*     individual characters within the string variable when the
*     string variable is a vector in the calling program (i.e. it
*     doens't seem possible to use the ':' notation when the character
*     variable is a vector - we get around this by passing a single
*     element of the vector to this function).
*
************************************************************************
      INTEGER*4 FUNCTION GETINT(STRVAR,ISTART,IEND)
*
*     function arguments
*
      CHARACTER*(*) STRVAR
      INTEGER*4 ISTART,IEND
      READ(STRVAR(ISTART:IEND),*) GETINT

      RETURN
      END

max_stars_repo_name: lthiamodelers/baseflow-coefficients

max_issues_repo_name: lthiamodelers/baseflow-coefficients

max_forks_repo_name: lthiamodelers/baseflow-coefficients


text:
Describe Users/dianacoss here.

Hi Diana, is the University Farm Circle no longer doing the Moms Groups walks? I reverted the change because the page used to have several tidbits of information and links to it, but I see you removed mention of it at the UFC page too.  users/edwins ES

20100218 22:44:38 nbsp Hi, Yes, could you please delete this UFC moms group, it no longer exists, I tried to do it months ago, and someone brought to my attention that it was still up, unfortunately, she waited at Borders to go on the walk, so I dont want that to happen to anyone else. Users/dianacoss

20100218 22:57:55 nbsp You could mark it as no longer existing by putting {{{Include(NoLongerExists)}}} on the page. Users/JasonAller


subroutine VEL (MONCE,IONCE,StabCriteria,YCOUNT,Obst_Hgt,dual_regime, &
  zo_patch,VEGVELinit_vel,PSLCALINIT)
  use simsphere_mod
  implicit none

!  Subroutine VEL computes the Monin Obukhov Length, the Friction
!  Velocity and the Integral of Heat Diffusivity.

! Code Altered 5th May 1992 ... Code transferred to Bri.for
! Alterations 16th July to account for different Roughness Lenghts
! associated with partial vegetation calculations.

  real :: KMM,KW,KS,KX
  real(kind=4) :: Obst_Hgt,zo_patch
  integer(kind=1) :: StabCriteria
  logical :: dual_regime
  integer :: MONCE, IONCE
  real :: REFLEV, ZTEN, CMH, CMW, YCOUNT
  real :: USCRN, RZAZO, U_Patch, Ustar_Patch, RZA_Obst, RObst_Patch
  real :: PTMP20, BOWEN, RTRANS, RTRANO3
  real :: SA, SRFCZO, SO10M, STEN, FM, FTEN, USTAR1
  real :: CHIO, CHIA, CHI20, T_ft, FT20, SObst_Hgt, FObst_Hgt
  real :: CH_Obst_Hgt, FT_Obst_Hgt, RZA_Obst_Hgt
  real :: PSIHNEW, REKUST
  real :: Rtrans_patch, RtransW_patch

  integer(kind=1) :: VEGVELinit_vel
  integer :: PSLCALINIT

!      INCLUDE 'modvars.h'

!  ZO Roughness height, ZA Top of surface layer (50m)

!  ZTEN - height at 10 m, REFLEV - "Screen or Anemometer Height".

  DATA REFLEV, ZTEN, KS, KW / 2., 10., 2.49951e-2, 2.97681e-2 /
  DATA CMH,CMW / 2*1. /

!  The model assumes neutral conditions at the start of the run where
!  HEAT = 0.  Therefore calc surface wind profile and resistances for
!  the surface layer on the basis log wind profile law.

  IF (StabCriteria .eq. 0 .or. (StabCriteria .eq. 1 .and.               &
     Heat <= 0.0 )) THEN ! Neutral Profile (Rare)

!         ^ Allow for negative heat flux

    USTAR = You_Star (AWIND,ZA,ZO,0.0,KARMAN)      ! Friction Velocity

    USCRN = WindF (USTAR,Reflev,ZO,0.0,KARMAN)      ! Wind Speeds
    UTEN =  WindF (USTAR,ZTEN,ZO,0.0,KARMAN)

    RZAZO = R_ohms (USTAR,ZA,ZO,0.0,KARMAN)        ! Resistances
    RZASCR = R_ohms (USTAR,ZA,REFLEV,0.0,KARMAN)

  if (dual_regime) then

     U_Patch = WindF (USTAR,Obst_Hgt,ZO,0.0,KARMAN)
     Ustar_Patch = You_Star (U_Patch,Obst_Hgt,zo_patch,0.0,KARMAN)
     RZA_Obst = R_ohms (USTAR,ZA,Obst_Hgt,0.0,KARMAN)
     RObst_Patch = R_ohms (Ustar_Patch,Obst_Hgt,zo_patch,0.0,KARMAN)

  end if

!  Potential, actual temp, specific humidity at 2 m. Pass
!  to vegetation component.

  PTMP20 = APTEMP + ( HEAT * RZASCR / ( DENS * CP ) )
  TA = PTMP20 - Tdif_s
  QA = AHUM + ( EVAP * RZASCR / ( DENS * LE ) )

!  Unstable case ... use the nondimensional functions of Z/L in
!  wind profiles for momentum and heat (FM,FT) to determine the
!  resistance term over the surface layer.

  else if ( StabCriteria .eq. 1 .and. heat > 0) then ! Unstable
    IF (EVAP < 0.000001) EVAP = 0.000001
    BOWEN = HEAT / EVAP
    IF ( ABS( BOWEN ) > 10 ) BOWEN = 10 * BOWEN / ABS ( BOWEN )
    MOL =  - ( ( USTAR**3 ) * aptemp * CP * DENS ) /                    &
           ( KARMAN * GRAV * HEAT * ( 1 + 0.07 / ( ABS ( BOWEN ) ) ) )

!  Dimensionless wind shear

    SA = Stab (ZA,MOL)
    SRFCZO = Stab (ZO,MOL)
    SO10M = Stab (REFLEV,MOL)
    STEN = Stab (ZTEN,MOL)

!  Stability Correction for momentum. Benoit solution.

    FM = FstabM (SRFCZO, SA)
    FTEN = FstabM (SRFCZO, STEN)

    USTAR1 = You_Star (AWIND,ZA,ZO,FM,KARMAN) ! Friction Velocity
    USTAR = ( USTAR + USTAR1 ) / 2 ! Smooth

! Dimensionless temperature gradient again using the integrated form.

    CHIO = StabH (ZO,MOL)
    CHIA = StabH (ZA,MOL)
    CHI20 = StabH (REFLEV,MOL)

!  Stability Correction for heat. Benoit solution.

    T_ft = FstabH (CHIO,CHIA)
    FT20 = FstabH (CHI20,CHIA)

    UTEN = WindF (USTAR,ZTEN,ZO,FTEN,KARMAN)       ! Surface Winds at
    USCRN = WindF (USTAR,REFLEV,ZO,0.0,KARMAN)     ! Screen Level & 10 metres

    RZASCR = R_ohms (USTAR,ZA,REFLEV,FT20,KARMAN) ! Resistances in
    RZAZO = R_ohms (USTAR,ZA,ZO,T_ft,KARMAN)  ! Surface Layer

    if (dual_regime) then                   ! As above

      Sobst_Hgt = Stab (Obst_Hgt,MOL)
      Fobst_Hgt = FstabM (SRFCZO, Sobst_Hgt)
      U_Patch = WindF (USTAR,Obst_Hgt,ZO,0.0,KARMAN)
      Ustar_Patch = You_Star(U_Patch,Obst_Hgt,zo_patch,Fobst_Hgt,KARMAN)

      CH_Obst_Hgt = StabH (Obst_Hgt,MOL)
      FT_Obst_Hgt = FstabH (CH_Obst_Hgt,CHIA)
      RZA_Obst_Hgt = R_ohms (USTAR,ZA,Obst_Hgt,FT_Obst_Hgt,KARMAN)
      RObst_Patch = R_ohms (Ustar_Patch,Obst_Hgt,zo_patch,0.0,KARMAN)

    end if

!   As in neutral case calculate TA and QA.

    PTMP20 = APTEMP + ( HEAT * RZASCR / ( DENS * CP ) )
    TA = PTMP20 - Tdif_s
    QA = AHUM + ( EVAP * RZASCR / ( DENS * LE ) )

  else

!  Call the nightime routine.  Stable case

    CALL BRI (MONCE,PSIHNEW,YCOUNT,ZTEN) ! Stable

    REKUST = 1 / (KARMAN * USTAR)
    RZAZO =  REKUST * ( ALOG( ZA / ZO ) + psihnew )

!       RZASCR = REKUST * ( ALOG( ZA / REFLEV ) + psihref )

!  note that stability correction function corresponds to
!  old method described in Panofsky and Dutton.  New function
!  employed in Bri.for.  Differences between two are
!  slight during the short period following radiation sunset
!  when the above functions are employed.

    TA = T1
    QA = AHUM

  END IF

!  Set USTAR equal to some non-zero value if small.

  IF ( USTAR < 0.01 ) USTAR = 0.01
! TJC Work around a problem with this being unitialized.
! TJC Related to dual_regime and zo_patch
  ustar_patch = ustar

!  Calc the diffusivities and resistances to heat and water for
!  the transition layer.  KS & KW are the molecular conductivities
!  for heat and water.  CMH and CMW are the scaling factors .. see
!  manual.

  IF (IONCE == 0) THEN
    KMM = CMH * KS / ( DENS * CP )
    KX =  CMW * KW / ( DENS * CP )
  END IF

  RTRANS = ResTrn (USTAR,ZO,KMM,KARMAN) ! Resistance Transition Layer
  RTRANW = ResTrn (USTAR,ZO,KX,KARMAN)
  RTRANO3 = ResTrn (USTAR,ZO,KX/1.32,KARMAN)

!  Finally compute the total series resistance over both layers.

  GBL_sum = RZAZO + RTRANS
  SUMW = RZAZO + RTRANW
  sumo3 = rzazo + rtrano3

  If( dual_regime) then
    Rtrans_Patch = ResTrn (ustar_patch,ZO,KMM,KARMAN)
    RtransW_Patch = ResTrn (ustar_patch,ZO,KX,KARMAN)

    GBL_sum = RZA_Obst_Hgt + RObst_Patch + Rtrans_Patch
    SUMW = RZA_Obst_Hgt + RObst_Patch + RtransW_Patch
  End if

!  If vegetation included, call vegetation component.

  IF (SWAVE > 0 .AND. RNET > 0 .AND. FRVEG > 0) THEN
    CALL VEGVEL(VEGVELinit_vel, PSLCALINIT)
  END IF

  return

end subroutine


text:
Davis Tree Works is one of the local businesses offering tree services.

20131012 15:32:09 nbsp Daniel was recommended by JAGCO Roofing & Painting for removal of a tree that unfortunately grew in a bad location. He did an excellent job including stump grinding and was very reasonably priced. I was impressed enough to create this wiki page. Previously I had used Kiwi Tree prior to Don leaving for a cooler climate, but Im very happy to have been referred to Daniel. Users/Dr.Baker

20131225 18:03:22 nbsp I needed a tree in my backyard cut back from my house and one removed. Daniel came out next morning, gave me a very competitive estimate and did the work within a few days. He did a great job. I can easily recommend him. Users/VeronikaMonell


text:
subroutine comp_emcoefs(r1_0,r2_0,r1_2,r2_2,j1,j2,th,nu0,Aji,nsI_tot,nsQ_tot,es_I,es_Q)
use fwigxjpf
!subroutine to compute the emission coefficients
!as of now, only implementation for the ratio nq/ni and eq/ei
!theta angle is required
!tested on a matlab routine

implicit none

integer j1,j2

double precision th,tki,tkq
double precision wjj2_1,wjj2_2
double precision fas,fac1,fac2,w6_1,w6_2
double precision r1_0,r2_0,r1_2,r2_2
double precision g1,g2,nu0,b_to_a,a_to_b
double precision nsI_tot,nsQ_tot
double precision na_I,na_Q,ns_I,ns_Q
double precision es_I,es_Q
double precision pQ,eQ,Aji,Bji,hv_fac

!input
!r1_0,2		-- lower level rank0,2 population
!r2_0,2		-- upper level rank0,2 population
!j1/j2		-- lower/upper level angular momentum
!th		-- propagation-magnetic field angle
!nu0		-- transition frequency

!output
!oQ		-- \eta_Q/\eta_I -factor
!eQ		-- \epsilon_Q / \epsilon_I - factor

!compute the w_{jj'}^k-factors
!only required for k=0 and k=2
!Eq.(38) Landi-84

fas = (-1.d0)**(1.d0 + j1 + j2)
fac1 = dsqrt(3.d0 * (2.d0 * j1 + 1.d0))
fac2 = dsqrt(3.d0 * (2.d0 * j2 + 1.d0))

w6_1 = fwig6jj(2*  1 , 2*  1 , 2*  2 ,&
     &         2* j1,  2*  j1, 2*  j2)

w6_2 = fwig6jj(2*  1 , 2*  1 , 2*  2 ,&
     &         2* j2,  2*  j2, 2*  j1)

wjj2_1 = fas * fac1 * w6_1
wjj2_2 = fas * fac2 * w6_2

!degeneracy
g1 = 2.d0 * j1 + 1.d0
g2 = 2.d0 * j2 + 1.d0

!
!T^K_Q(i,\Omega) rotation elements from Bommier 1997b
tki = (3.d0 * dcos(th)**2.d0 - 1.d0 )/(2.d0 * dsqrt(2.d0))
tkq = -3.d0 * dsin(th)**2.d0 / (2.d0 * dsqrt(2.d0))

!write(*,*)
!write(*,*)tki,tkq
!write(*,*)

!Eqs.(40)-(42) Landi-84
na_I = 1.d0 + wjj2_1 * (r1_2/r1_0) * tki
na_Q = wjj2_1 * (r1_2/r1_0) * tkq

ns_I = 1.d0 + wjj2_2 * (r2_2/r2_0) * tki
ns_Q = wjj2_2 * (r2_2/r2_0) * tkq

a_to_b = 6.781962D49    !best to use L&L relation here
                        !B21 = [c^2 / (2*h * v^3)] * A21
                        !a_to_b = c^2 / (2*h)

!write(*,*)na_I,na_Q
!write(*,*)ns_I,ns_Q

b_to_a = nu0**3.d0 / a_to_b

Bji = Aji / b_to_a
!write(*,*)
!write(*,*)Aji,Bji

es_I = ns_I
es_Q = ns_Q

!now properly substract the absorption and stimulated emission
!so that the proportionalities stay in-tact

!na ~ B_12 * N1
!ns ~ B_21 * N2 = B_12 * g1 * N2 / g2
!N1 ~ sqrt(g1) * r1_0
!N2 ~ sqrt(g2) * r2_0 		--- Eq.(43) Landi-84

nsI_tot = dsqrt(g1)*r1_0 * na_I - dsqrt(g2)*r2_0*ns_I*g1/g2
nsQ_tot = dsqrt(g1)*r1_0 * na_Q - dsqrt(g2)*r2_0*ns_Q*g1/g2

hv_fac = 5.273D-35 * nu0            !(h/4*pi) * nu0

!write(*,*)
!write(*,*)hv_fac

es_I = hv_fac * Aji * es_I * dsqrt(g2) * r2_0
es_Q = hv_fac * Aji * es_Q * dsqrt(g2) * r2_0

nsI_tot = hv_fac * (g2/g1)*Bji * nsI_tot        !g2/g1 * B21 = B12
nsQ_tot = hv_fac * (g2/g1)*Bji * nsQ_tot

!write(*,*)nsI_tot,nsQ_tot
!write(*,*)es_I,es_Q
!stop
end subroutine

subroutine get_rhovec(c_dim,count_array,nlev,j_lev,ntrans,trans,AC_trans,t_int,kmax,Temp,D0,rhovec,extra_info,rates)
!subroutine that obtains the rho-vector of a certain node
!by setting up the SEE and subsequently solving them
implicit none
integer nlev,ntrans,t_int,c_dim
integer i,j,k,jc,ic,INFO,n
integer kmax,extra_info

integer, dimension(nlev) :: j_lev,count_array
integer, dimension(2,ntrans) :: trans
integer, dimension(c_dim) :: IPIV

double precision Temp

double precision, dimension(5,ntrans) :: AC_trans
double precision, dimension(nlev,2:kmax) :: D0
double precision, dimension(c_dim,c_dim) :: Amat
double precision, dimension(c_dim+1) :: yvec
double precision, dimension(c_dim+1,c_dim) :: Bmat
double precision, dimension(c_dim+1) :: dvec
double precision, dimension(c_dim) :: WORK,rhovec
double precision, dimension(c_dim,c_dim) :: C_mat
double precision, dimension(c_dim,c_dim) :: B_mat
double precision, dimension(c_dim,c_dim) :: R_mat
double precision, dimension(3,c_dim) :: rates

!input:
!c_dim          -- dimension of the SEE matrix
!count_array    -- bookkeeping array for the (polarized) levels
!nlev           -- number of levels
!j_lev          -- associated J-angular momentum of the levels
!ntrans         -- number of transitions
!trans          -- associated levels with these transitions
!AC_trans       -- associated rates with these transtions
!t_int          --

call setup_stateq(c_dim,count_array,nlev,j_lev,trans,ntrans,kmax,&
&AC_trans,D0,Temp,Amat,R_mat,B_mat,C_mat)

!like in Hazel, replace the first row with the level constraint
!that is set at 1d7 --- densities we are interested in are
!probably lower, so we'll set 1.d3

Amat(1,:) = 0.d0
do i = 1,nlev
  ic = count_array(i) + 1
  jc = j_lev(i)

  Amat(1,ic) = dsqrt(2.d0 * jc + 1.d0)
end do

!Amat * rhovec = dvec
!dvec(1) = 1.d3 and rest = 0.d0

call DGETRF(c_dim,c_dim,Amat,c_dim,IPIV,INFO)
if (INFO.ne.0) then
  write(*,*)'Density-matrix calculation: failed'
  stop
else
  call DGETRI(c_dim,Amat,c_dim,IPIV,WORK,c_dim,INFO)
  if (INFO.ne.0) then
    write(*,*)'Density-matrix calculation: failed'
    stop
  endif
endif

!see beneath to get all the populations, but we are actually
!interested only in the populations for our levels
!dvec(1) = 1.d3 and rest = 0.d0
do i = 1,c_dim
  rhovec(i) = Amat(i,1)
end do

if (extra_info.eq.1) then
  rates(:,:) = 0.d0
  do i = 1,c_dim
    do j = 1,c_dim
      rates(1,i) = rates(1,i) + dabs(R_mat(i,j))*rhovec(j)
      rates(2,i) = rates(2,i) + dabs(B_mat(i,j))*rhovec(j)
      rates(3,i) = rates(3,i) + dabs(C_mat(i,j))*rhovec(j)
!      write(*,*)j,dabs(C_mat(i,j)),rhovec(j)
    end do
!    stop
  end do
endif

end subroutine

max_stars_repo_name: blankhaar/PORTAL

max_issues_repo_name: blankhaar/PORTAL

max_forks_repo_name: blankhaar/PORTAL


text:
The UC Davis Bowling Club is a student organizations student club on campus at UC Davis. They practice in the MU Games Area.
The mens club sport team is led by captain Jesse Anzelc and finished their 20112012 season ranked 35th in the nation, out of about 150 teams.  The team earned two collegiate titles this season, the Camarillo Open (November 2011) and the Silicon Valley Classic (February 2012).


program exercise_3_4

!сумма элементов массива с индексами с меняющимся шагом = 3, начиная от 5го элемента
   implicit none

   character(*), parameter     :: input_file = "../data/input.txt", output_file = "output.txt"
   integer                     ::In = 0, Out = 0, N=0
   real(4)                     :: SumA=0
   real(4), allocatable        ::  A (:)

open (file=input_file, newunit=In)
   read (In, *) N  ! число -  размерность массива
   allocate (A(N))
   read (In, *) A  !читаем N  следующиx элементов из файла
close (In)

   SumA = Sum(A(5:N:3)) !Правильно - экономия памяти. в кеш памяти хранится только каждый 3й элемент. перебор от 5го элемента до  элемента N c шагом 3

   open (file=output_file, newunit=Out)
      write (Out, '(F0.2)') A  !оставить 2 знака после запятой
      write (Out, *) SumA
   close (Out)

end program exercise_3_4


          program demo_chomp

          use M_strings, only : chomp
          implicit none
          character(len=100)            :: inline
          character(len=:),allocatable  :: token
          character(len=*),parameter    :: delimiters=' ;,'
          integer                       :: ios
          integer                       :: icount
          integer                       :: itoken
             icount=0
             do        ! read lines from stdin until end-of-file or error
                read (unit=*,fmt="(a)",iostat=ios) inline
                if(ios.ne.0)stop
                icount=icount+1
                itoken=0
                write(*,*)'INLINE ',trim(inline)
                do while ( chomp(inline,token,delimiters).ge. 0)
                   itoken=itoken+1
                   print *, itoken,'TOKEN=['//trim(token)//']'
                enddo
             enddo

          end program demo_chomp


text:
!##############################################################################
!# Tutorial 005a: Basic math
!##############################################################################

module tutorial005a

  ! Include basic Feat-2 modules
  use fsystem
  use genoutput
  use mprimitives

  implicit none
  private

  public :: start_tutorial005a

contains

  ! ***************************************************************************

  subroutine start_tutorial005a

    ! Declare some variables
    real(DP) :: dval1, dval2, dval3
    real(DP), dimension(4) :: Dcoeff

    ! Print a message
    call output_lbrk()
    call output_separator (OU_SEP_STAR)
    call output_line ("This is FEAT-2. Tutorial 005a")
    call output_separator (OU_SEP_MINUS)

    ! =================================
    ! Kronecker Symbol
    ! =================================

    call output_line ("kronecker(3,5) = " // trim(sys_siL(mprim_kronecker(3,5),10)) )
    call output_line ("kronecker(5,5) = " // trim(sys_siL(mprim_kronecker(5,5),10)) )

    ! =================================
    ! Signum
    ! =================================

    call output_line ("signum(-5) = " // trim(sys_siL(mprim_signum(-5_I32),10)) )
    call output_line ("signum( 0) = " // trim(sys_siL(mprim_signum( 0_I32),10)) )
    call output_line ("signum( 2) = " // trim(sys_siL(mprim_signum( 2_I32),10)) )

    ! =================================
    ! Linear rescaling
    ! =================================

    ! Linear interpolation of [0,1] to [5,7]
    call mprim_linearRescale(0.00_DP,0.0_DP,1.0_DP,5.0_DP,7.0_DP,dval1)
    call mprim_linearRescale(0.25_DP,0.0_DP,1.0_DP,5.0_DP,7.0_DP,dval2)
    call mprim_linearRescale(1.00_DP,0.0_DP,1.0_DP,5.0_DP,7.0_DP,dval3)

    call output_line ("linrescale(0.00,[0,1],[5,7]) = " // trim(sys_sdL(dval1,10)) )
    call output_line ("linrescale(0.25,[0,1],[5,7]) = " // trim(sys_sdL(dval2,10)) )
    call output_line ("linrescale(1.00,[0,1],[5,7]) = " // trim(sys_sdL(dval3,10)) )

    ! =================================
    ! Quadratic interpolation
    ! =================================

    ! Quadratic interpolation, i.e. evaluate a quadratic
    ! polynomial p with p(-1)=2, p(0)=0, p(1)=2
    call mprim_quadraticInterpolation (-1.0_DP,2.0_DP,0.0_DP,2.0_DP,dval1)
    call mprim_quadraticInterpolation ( 0.0_DP,2.0_DP,0.0_DP,2.0_DP,dval2)
    call mprim_quadraticInterpolation ( 0.5_DP,2.0_DP,0.0_DP,2.0_DP,dval3)

    call output_line ("quadinterpol(-1.0,[2,0,2]) = " // trim(sys_sdL(dval1,10)) )
    call output_line ("quadinterpol( 0.0,[2,0,2]) = " // trim(sys_sdL(dval2,10)) )
    call output_line ("quadinterpol( 0.5,[2,0,2]) = " // trim(sys_sdL(dval3,10)) )

    ! =================================
    ! Parabolic profile
    ! =================================

    ! This is a quadratic polynomial p with p(0)=p(len)=0, p(mid)=max.
    dval1 = mprim_getParabolicProfile (0.5_DP,2.0_DP,10.0_DP)
    dval2 = mprim_getParabolicProfile (1.0_DP,2.0_DP,10.0_DP)
    dval3 = mprim_getParabolicProfile (1.5_DP,2.0_DP,10.0_DP)

    call output_line ("parprofile(0.5,len=2,max=10) = " // trim(sys_sdL(dval1,10)) )
    call output_line ("parprofile(1.0,len=2,max=10) = " // trim(sys_sdL(dval2,10)) )
    call output_line ("parprofile(1.5,len=2,max=10) = " // trim(sys_sdL(dval3,10)) )

    ! =================================
    ! Conversion of polar coordinates
    ! in cartesian coordinates
    ! =================================

    call mprim_polarToCartesian(0.5_DP,0.0_DP,dval1,dval2)
    call output_line ("polar(0.5,0.0    )  = " // trim(sys_sdL(dval1,10)) // "," &
                                               // trim(sys_sdL(dval2,10)))

    call mprim_polarToCartesian(0.5_DP,0.25_DP*SYS_PI,dval1,dval2)
    call output_line ("polar(0.5,0.25 Pi)  = " // trim(sys_sdL(dval1,10)) // "," &
                                               // trim(sys_sdL(dval2,10)))

    call mprim_polarToCartesian(0.5_DP,-0.75*SYS_PI,dval1,dval2)
    call output_line ("polar(0.5,-0.75 Pi) = " // trim(sys_sdL(dval1,10)) // "," &
                                               // trim(sys_sdL(dval2,10)))

    ! =================================
    ! Conversion of cartesian coordinates
    ! in polar coordinates
    ! =================================

    call mprim_cartesianToPolar(0.5_DP,0.0_DP,dval1,dval2)
    call output_line ("cartesian(0.5,0.0)         = " // trim(sys_sdL(dval1,10)) // "," &
                                                      // trim(sys_sdL(dval2,10)))

    call mprim_cartesianToPolar(0.25_DP*sqrt(2.0_DP),0.25_DP*sqrt(2.0_DP),dval1,dval2)
    call output_line ("cartesian(0.3535,0.3535)   = " // trim(sys_sdL(dval1,10)) // "," &
                                                      // trim(sys_sdL(dval2,10)))

    call mprim_cartesianToPolar(-0.25_DP*sqrt(2.0_DP),-0.25_DP*sqrt(2.0_DP),dval1,dval2)
    call output_line ("cartesian(-0.3535,-0.3535) = " // trim(sys_sdL(dval1,10)) // "," &
                                                      // trim(sys_sdL(dval2,10)))

    ! =================================
    ! Horner scheme
    ! =================================

    ! Initialise a cubic polynomial: p(x,y) = x^3 - 2x^2 - 3*x + 4
    Dcoeff(:) = (/ 4.0_DP, 3.0_DP, -2.0_DP, 1.0_DP /)

    ! Get the function value as well as 1st and 2nd derivative in x=2.
    dval1 = mprim_horner   (2.0_DP,Dcoeff)
    dval2 = mprim_hornerd1 (2.0_DP,Dcoeff)
    dval3 = mprim_hornerd2 (2.0_DP,Dcoeff)

    call output_line ("p  (2.0) = " // trim(sys_sdL(dval1,10)) )
    call output_line ("p' (2.0) = " // trim(sys_sdL(dval2,10)) )
    call output_line ("p''(2.0) = " // trim(sys_sdL(dval3,10)) )

  end subroutine

end module

max_stars_repo_licenses: Intel Unlicense

max_issues_repo_licenses: Intel Unlicense

max_forks_repo_licenses: Intel Unlicense


text:
      program main
c---- Целые числа от 2 до 100: решето Эратосфена
      integer prime (100)
c---- Обнулим prime
      do 10 n=2, 100
 10       prime (n) = 1
c---- Вычислим prime
      do 20 n=2, 100
          if (prime(n) .eq. 0) goto 20
              maxk = 100 / n
	      if (maxk .lt. 2) goto 20
              do 30 k=2, maxk
 30               prime (n*k) = 0
 20       continue
c---- Печать prime
      do 40 n=2, 100
 40       if (prime(n) .ne. 0) print 1000, n
      stop
 1000 format (1x, i4)
      end


text:
c
c
c     ###################################################
c     ##  COPYRIGHT (C)  1993  by  Jay William Ponder  ##
c     ###################################################
c
c     #############################################################
c     ##                                                         ##
c     ##  subroutine eangang  --  angle-angle cross term energy  ##
c     ##                                                         ##
c     #############################################################
c
c
c     "eangang" calculates the angle-angle potential energy
c
c
      subroutine eangang
      implicit none
      include 'sizes.i'
      include 'angang.i'
      include 'angle.i'
      include 'angpot.i'
      include 'atoms.i'
      include 'bound.i'
      include 'energi.i'
      include 'group.i'
      include 'math.i'
      include 'usage.i'
      integer i,k,iangang
      integer ia,ib,ic,id,ie
      real*8 e,dt1,dt2,fgrp
      real*8 angle,dot,cosine
      real*8 xia,yia,zia
      real*8 xib,yib,zib
      real*8 xic,yic,zic
      real*8 xid,yid,zid
      real*8 xie,yie,zie
      real*8 xab,yab,zab
      real*8 xcb,ycb,zcb
      real*8 xdb,ydb,zdb
      real*8 xeb,yeb,zeb
      real*8 rab2,rcb2
      real*8 rdb2,reb2
      logical proceed
c
c
c     zero out the angle-angle cross term energy
c
      eaa = 0.0d0
c
c     calculate the angle-angle interaction energy term
c
      do iangang = 1, nangang
         i = iaa(1,iangang)
         k = iaa(2,iangang)
         ia = iang(1,i)
         ib = iang(2,i)
         ic = iang(3,i)
         id = iang(1,k)
         ie = iang(3,k)
c
c     decide whether to compute the current interaction
c
         proceed = .true.
         if (use_group)  call groups (proceed,fgrp,ia,ib,ic,id,ie,0)
         if (proceed)  proceed = (use(ia) .or. use(ib) .or. use(ic)
     &                               .or. use(id) .or. use(ie))
c
c     get the coordinates of the atoms in the angle
c
         if (proceed) then
            xia = x(ia)
            yia = y(ia)
            zia = z(ia)
            xib = x(ib)
            yib = y(ib)
            zib = z(ib)
            xic = x(ic)
            yic = y(ic)
            zic = z(ic)
            xid = x(id)
            yid = y(id)
            zid = z(id)
            xie = x(ie)
            yie = y(ie)
            zie = z(ie)
c
c     compute the values of the two bond angles
c
            xab = xia - xib
            yab = yia - yib
            zab = zia - zib
            xcb = xic - xib
            ycb = yic - yib
            zcb = zic - zib
            xdb = xid - xib
            ydb = yid - yib
            zdb = zid - zib
            xeb = xie - xib
            yeb = yie - yib
            zeb = zie - zib
            if (use_polymer) then
               call image (xab,yab,zab)
               call image (xcb,ycb,zcb)
               call image (xdb,ydb,zdb)
               call image (xeb,yeb,zeb)
            end if
            rab2 = xab*xab + yab*yab + zab*zab
            rcb2 = xcb*xcb + ycb*ycb + zcb*zcb
            rdb2 = xdb*xdb + ydb*ydb + zdb*zdb
            reb2 = xeb*xeb + yeb*yeb + zeb*zeb
            if (rab2*rcb2*rdb2*reb2 .ne. 0.0d0) then
               dot = xab*xcb + yab*ycb + zab*zcb
               cosine = dot / sqrt(rab2*rcb2)
               cosine = min(1.0d0,max(-1.0d0,cosine))
               angle = radian * acos(cosine)
               dt1 = angle - anat(i)
               dot = xdb*xeb + ydb*yeb + zdb*zeb
               cosine = dot / sqrt(rdb2*reb2)
               cosine = min(1.0d0,max(-1.0d0,cosine))
               angle = radian * acos(cosine)
               dt2 = angle - anat(k)
c
c     get the angle-angle interaction energy
c
               e = aaunit * kaa(iangang) * dt1 * dt2
c
c     scale the interaction based on its group membership
c
               if (use_group)  e = e * fgrp
c
c     increment the total angle-angle energy
c
               eaa = eaa + e
            end if
         end if
      end do
      return
      end

max_stars_repo_name: htna/HCsbLib

max_issues_repo_name: htna/HCsbLib

max_forks_repo_name: htna/HCsbLib