    program dyndot
    use m_dyn
    use m_ioutil, only : luavail
    use m_die, only: die
    implicit none
    character(len=*), parameter :: myname='dyndot'
    integer,parameter :: dyntype=5
    integer,parameter :: nfiles=2
    integer nymd, nhms, lu, n, freq, vectype, prec, ier, nstep
    integer i, nf, iarg, argc, ndim2, ndim3, intarg, iargc
    character(len=255) :: dynfile(nfiles)
    character(len=255) argv
    type(dyn_vect) w_1
    type(dyn_vect) w_2
    logical use_ps
    real :: total, dotp(5)
    integer im,jm,km,lm,k,rc
    integer im1,jm1,km1,lm1
    integer im2,jm2,km2,lm2
    logical adm(nfiles)
    logical pncf

    prec = 0
    vectype = 5
    use_ps=.false. ! default is to use delp
    adm=.false.
    pncf=.false.

    argc = iargc()
    if ( argc < 1 ) call usage_()
    nf=0; iarg=0
    do i = 1, 32767
       iarg = iarg + 1
       if ( iarg .gt. argc ) exit
       call GetArg ( iarg, argv )
       select case (argv)
          case ("-use_ps")
             use_ps = .true.
          case ("-adm1")
             adm(1) = .true.
          case ("-adm2")
             adm(2) = .true.
          case ("-pncf")
             pncf = .true.
!         case ("-adm")
!            if ( iarg+1 .gt. argc ) call usage_()
!            iarg = iarg + 1
!            call GetArg ( iarg, argv )
!            read(argv,*) intarg
!            if(intarg>0) adm(1)=.true.
!            iarg = iarg + 1
!            call GetArg ( iarg, argv )
!            read(argv,*) intarg
!            if(intarg>0) adm(2)=.true.
          case default
             nf = nf + 1
             if ( nf .gt. nfiles ) call die(myname,'too many eta files')
             dynfile(nf) = trim(argv)
       end select
    enddo

    do nf=1,nfiles
       write(6,'(a,i3,2a)') "Dyn Vector ", nf , ": ", trim(dynfile(nf))
    enddo

    call dyn_getdim ( trim(dynfile(1)), im1, jm1, km1, lm1, rc )
    call dyn_getdim ( trim(dynfile(2)), im2, jm2, km2, lm2, rc )
    if(km1/=km2) then ! ignore diff in lm for now
       print *, trim(myname), ': km/lm file 1 = ', km1,lm1
       print *, trim(myname), ': km/lm file 2 = ', km2,lm2
       call die(myname,'inconsistent km/lm')
    endif
    im=min(im1,im2)
    jm=min(jm1,jm2)
    km=min(km1,km2)
    lm=min(lm1,lm2)
    ndim2 = im*jm

!   check dimensions and take proper action
    n=1
    if ( im1==im2 .and. jm1==jm2 ) then
       call dyn_get ( trim(dynfile(1)), nymd, nhms, w_1, ier, timidx=n, &
                                      freq=freq, nstep=nstep, vectype=vectype, pncf=pncf )
       call dyn_get ( trim(dynfile(2)), nymd, nhms, w_2, ier, timidx=n, &
                                      freq=freq, nstep=nstep, vectype=vectype, pncf=pncf )
    else

       if (im1<im2 .and. jm1<jm2 ) then

          call dyn_get ( trim(dynfile(1)), nymd, nhms, w_1, ier, timidx=n, &
                                         freq=freq, nstep=nstep, vectype=vectype, pncf=pncf )

!         Initialize dimension of output (interpolated) vector
!         ----------------------------------------------------
          call dyn_init ( im, jm, km, lm, w_2, rc, &
                          w_1%grid%ptop, w_1%grid%ks, w_1%grid%ak, w_1%grid%bk, vectype=dyntype )
            if ( rc/=0 ) then
                 call die (myname, ': Error initializing dyn vector(w_2)')
            endif

          call getvector_ ( dynfile(2), w_2, adm(1) )

!         Recalculate ps for consistency w/ delp
!         --------------------------------------
          w_2%ps = 0.0 ! perturbations
          do k=1,km
             w_2%ps = w_2%ps + w_2%delp(:,:,k)
          enddo

       else

          call dyn_get ( trim(dynfile(2)), nymd, nhms, w_2, ier, timidx=n, &
                                         freq=freq, nstep=nstep, vectype=vectype, pncf=pncf )

!         Initialize dimension of output (interpolated) vector
!         ----------------------------------------------------
          call dyn_init ( im, jm, km, lm, w_1, rc, &
                          w_2%grid%ptop, w_2%grid%ks, w_2%grid%ak, w_2%grid%bk, vectype=dyntype )
            if ( rc/=0 ) then
                 call die (myname, ': Error initializing dyn vector(w_1)')
            endif

          call getvector_ ( dynfile(1), w_1, adm(2) )

       endif

    endif
    ndim3 = ndim2 * km
    call dotall_ (w_1,w_2)

    contains

    subroutine getvector_ (fname,x,adm)
    use m_maph_pert, only: h_map_pert
    implicit none
    character(len=*) fname
    type(dyn_vect) x
    type(dyn_vect) www
    logical adm

          call dyn_get ( trim(fname), nymd, nhms, www, rc, timidx=n, &
                         freq=freq, vectype=dyntype, pncf=pncf )

!         Interpolate to required resolution
!         ----------------------------------
          call h_map_pert ( www, x, rc, adm=adm )
               if ( rc/=0 ) then
                    call dyn_clean ( x )
                    call dyn_clean ( www )
                    print *, 'h_map error code = ', rc
                    call die(myname,' failed in h_map')
               else
                    print*, myname, ': interpolated perturbation: ', trim(fname)
                    print*, myname, ': from ', www%grid%im, 'x', www%grid%jm, ' to ', x%grid%im, 'x', x%grid%jm
                    call dyn_clean ( www )
               endif

    end subroutine getvector_

    subroutine dotall_ (x,y)
    type(dyn_vect) x,y
    real ddot_

    dotp(1) = ddot_(ndim3,x%u ,1,y%u ,1)
    dotp(2) = ddot_(ndim3,x%v ,1,y%v ,1)
    dotp(3) = ddot_(ndim3,x%pt,1,y%pt,1)
    if ( use_ps ) then
       dotp(4) = ddot_(ndim2,x%ps  ,1,y%ps  ,1)
    else
       dotp(4) = ddot_(ndim3,x%delp,1,y%delp,1)
    endif
    dotp(5) = ddot_(ndim3,x%q(:,:,:,1) ,1,y%q(:,:,:,1) ,1)  ! only 1st entry in q for now
    total   = sum(dotp)

    lu= luavail()
    open(lu, file='dyndot.txt', form='formatted')
    write (lu, '(1p,6e16.6)') total, dotp(1:5)
    close(lu)

    end subroutine dotall_

    subroutine usage_
       print *, "Purpose: Calculate dot product between two dyn vectors"
       print *, "         Needless to say user should know what he/she "
       print *, "         is calculation or else this is meaningless."
       print *, "   "
       print *, "Usage: dyndot.x [options] dynfile1 dynfile2 "
       print *, "   "
       print *, " Optional arguments:  "
       print *, "   "
       print *, " -use_ps    Specify to use PS instead of DELP in dot product "
       print *, " -adm1      Specify if 1st file is a sensitivity vector"
       print *, " -adm2      Specify if 2nd file is a sensitivity vector"
       print *, "   "
       print *, "Output written to file dyndot.txt"
       stop
    end subroutine usage_

    end program dyndot

!
! For some reason, this BLAS function cannot be found by the
! intel compiler on some Linux platforms. (Arlindo)
!

      real FUNCTION DDOT_(N,DX,INCX,DY,INCY)
!     .. Scalar Arguments ..
      INTEGER INCX,INCY,N
!     ..
!     .. Array Arguments ..
      real DX(*),DY(*)
!     ..
!
!  Purpose
!  =======
!
!     forms the dot product of two vectors.
!     uses unrolled loops for increments equal to one.
!     jack dongarra, linpack, 3/11/78.
!     modified 12/3/93, array(1) declarations changed to array(!)
!
!
!     .. Local Scalars ..
      real*8 DTEMP
      INTEGER I,IX,IY,M,MP1
!     ..
!     .. Intrinsic Functions ..
      INTRINSIC MOD
!     ..
      DDOT_ = 0.0d0
      DTEMP = 0.0d0
      IF (N.LE.0) RETURN
      IF (INCX.EQ.1 .AND. INCY.EQ.1) GO TO 20
!
!        code for unequal increments or equal increments
!          not equal to 1
!
      IX = 1
      IY = 1
      IF (INCX.LT.0) IX = (-N+1)*INCX + 1
      IF (INCY.LT.0) IY = (-N+1)*INCY + 1
      DO 10 I = 1,N
          DTEMP = DTEMP + DX(IX)*DY(IY)
          IX = IX + INCX
          IY = IY + INCY
   10 CONTINUE
      DDOT_ = DTEMP
      RETURN
!
!        code for both increments equal to 1
!
!
!        clean-up loop
!
   20 M = MOD(N,5)
      IF (M.EQ.0) GO TO 40
      DO 30 I = 1,M
          DTEMP = DTEMP + DX(I)*DY(I)
   30 CONTINUE
      IF (N.LT.5) GO TO 60
   40 MP1 = M + 1
      DO 50 I = MP1,N,5
          DTEMP = DTEMP + DX(I)*DY(I) + DX(I+1)*DY(I+1) + &
                  DX(I+2)*DY(I+2) + DX(I+3)*DY(I+3) + DX(I+4)*DY(I+4)
   50 CONTINUE
   60 DDOT_ = DTEMP
      RETURN
      END


text:
SUBROUTINE fc2d_geoclaw_update_gauge (mx,my,mbc,meqn,xlower,ylower, &
    dx,dy,q,maux,aux,xc,yc,qvar,avar)

    use geoclaw_module, only: dry_tolerance

    implicit none

    integer :: mx, my, mbc, meqn, maux
    real(kind=8) :: xlower, ylower, dx, dy, xc, yc
    real(kind=8) :: q(meqn, 1-mbc:mx+mbc, 1-mbc:my+mbc)
    real(kind=8) :: aux(maux, 1-mbc:mx+mbc, 1-mbc:my+mbc)
    real(kind=8) :: qvar(meqn), avar(maux)

    !! local variables:
    real(kind=8) :: xcent,ycent,xoff,yoff
    integer :: ioff,joff,iindex,jindex,mq
    real(kind=8) :: h(4),drytol2,topo

    iindex =  int((xc-xlower)/dx) + 1
    jindex =  int((yc-ylower)/dy) + 1

    xcent  = xlower + (iindex-.5d0)*dx
    ycent  = ylower + (jindex-.5d0)*dy
    xoff   = (xc-xcent)/dx
    yoff   = (yc-ycent)/dy

    drytol2 = 0.1d0 * dry_tolerance

    h(1) = q(1,iindex,jindex)
    h(2) = q(1,iindex+1,jindex)
    h(3) = q(1,iindex,jindex+1)
    h(4) = q(1,iindex+1,jindex+1)

    if ((h(1) < drytol2) .or.  &
        (h(2) < drytol2) .or.  &
        (h(3) < drytol2) .or.  &
        (h(4) < drytol2)) then
        !! One of the cells is dry, so just use value from grid cell
        !! that contains gauge rather than interpolating

        !! icell = int(1.d0 + (xc - xlower) / hx)
        !! jcell = int(1.d0 + (yc - ylower) / hy)
        do mq=1,3
            qvar(mq) = q(mq,iindex,jindex)
        enddo
        !! This is the bottom layer and we should figure out the
        !! topography
        topo = aux(1,iindex,jindex)
    else
        !! Linear interpolation between four cells
        do mq=1,3
            qvar(mq) = (1.d0 - xoff) * (1.d0 - yoff) &
                * q(mq,iindex,jindex)  &
                + xoff*(1.d0 - yoff) * q(mq,iindex+1,jindex)  &
                + (1.d0 - xoff) * yoff * q(mq,iindex,jindex+1)  &
                + xoff * yoff * q(mq,iindex+1,jindex+1)
        enddo
        topo = (1.d0 - xoff) * (1.d0 - yoff)  &
                  * aux(1,iindex,jindex)  &
                  + xoff * (1.d0 - yoff) * aux(1,iindex+1,jindex)  &
                  + (1.d0 - xoff) * yoff * aux(1,iindex,jindex+1)  &
                  + xoff * yoff * aux(1,iindex+1,jindex+1)
    endif

    ! Extract surfaces
    !!eta = qvar(1) + topo  !! Compute this when printing buffers
    avar(1) = topo

    ! Zero out tiny values to prevent later problems reading data,
    ! as done in valout.f
    do mq = 1,3
        if (abs(qvar(mq)) < 1d-99) qvar(mq) = 0.d0
    end do
    !! if (abs(eta) < 1d-99) eta = 0.d0

END SUBROUTINE fc2d_geoclaw_update_gauge

max_stars_repo_name: ECLAIRWaveS/ForestClaw

max_issues_repo_name: ECLAIRWaveS/ForestClaw

max_forks_repo_name: ECLAIRWaveS/ForestClaw


text:
program multi_svd
    use utils_mod, only : get_file_name, get_dataset_name
    use matrix_mod, only : read_matrix, print_matrix, compute_svd, &
                           reconstruct_matrix, compute_sigma
    use, intrinsic :: iso_fortran_env, only : dp => REAL64, &
                                              i8 => INT64, &
                                              error_unit
    use :: omp_lib
    implicit none
    logical, parameter :: verbose = .false.
    real(kind=dp), allocatable, dimension(:, :) :: A, U, VT, Sigma, orig_A
    real(kind=dp), allocatable, dimension(:) :: S
    character(len=1024) :: base_file_name, file_name, dataset_name
    integer :: nr_rows, nr_cols, error, thread_id

    call get_file_name(base_file_name)
    call get_dataset_name(dataset_name)
    !$omp parallel default(none) &
    !$omp          private(file_name, A, U, VT, Sigma, orig_A, S, &
    !$omp                  nr_rows, nr_cols, error, thread_id) &
    !$omp          shared(base_file_name, dataset_name)
    thread_id = omp_get_thread_num()
    write (file_name, '(A, I0, A)') trim(base_file_name), thread_id + 1, '.h5'
    if (verbose) then
        print '(A, I0, A, A)', 'thread ', thread_id, ' reading ', trim(file_name)
    end if
    !$omp critical
    call read_matrix(file_name, dataset_name, A)
    !$omp end critical
    nr_rows = size(A, 1)
    nr_cols = size(A, 2)

    ! compute SVD
    call compute_svd(A, U, S, VT)

    ! compute Sigma, which is nicer than S, and which is needed to
    ! reconstruct the matrix
    call compute_sigma(S, nr_rows, nr_cols, sigma)

    ! print results if verbose
    if (verbose) then
        !$omp critical
        print "(A)", "U"
        call print_matrix(U)
        print "(A)", "Sigma"
        call print_matrix(Sigma)
        print "(A)", "V^t"
        call print_matrix(VT)
        !$omp end critical
    end if

    ! copy original matrix
    allocate(orig_A, source=A, stat=error)
    if (error /= 0) then
    end if
    orig_A = A
    deallocate(A)

    call reconstruct_matrix(U, S, VT, A)

    ! print reconstructted and original matrix if verbose
    if (verbose) then
        !$omp critical
        print "(A)", "A"
        call print_matrix(A)
        print "(A)", "orig. A"
        call print_matrix(orig_A)
        !$omp end critical
    end if

    ! compute and print relative error
    print "(A, E10.3)", 'relative error = ', compute_error(orig_A, A)

    deallocate(A, orig_A, U, S, VT, Sigma)

    !$omp end parallel

contains

    function compute_error(orig_A, A) result(rel_err)
        use, intrinsic :: iso_fortran_env, only : dp => REAL64
        implicit none
        real(kind=dp), dimension(:, :), intent(in) :: orig_A, A
        real(kind=dp) :: rel_err
        integer :: i, j
        real(kind=dp) :: norm
        rel_err = 0.0_dp
        do j = 1, size(A, 2)
            do i = 1, size(A, 1)
                norm = abs(orig_A(i, j))
                rel_err = rel_err + abs(A(i, j) - orig_A(i, j))/norm
            end do
        end do
        rel_err = rel_err/(size(A, 1)*size(A, 2))
    end function compute_error

end program multi_svd

max_stars_repo_name: Gjacquenot/training-material

max_issues_repo_name: Gjacquenot/training-material

max_forks_repo_name: Gjacquenot/training-material


text:
       subroutine apomak
*++
*
* APOMASK V1.0 - create an apodising mask for image restoration
*
* This is an IRAF application written using the F77/VOS.
*
* Purpose:
*
* Create an image which is zero everywhere except for a central circular region
* where it is one. The region outside the central circle fades off as a gaussian.
* The user supplies the radius of the central circle and the sigma of the fade-off.
*
* History:
*
* First version, Richard Hook, ST-ECF, Jan 92
* Added limb darkening, Richard Hook, ST-ECF, Apr 94
* Renamed APOMASK, Richard Hook, Mar 99
* Given version number, changed parameter names to all lowercase,
*  Richard Hook, February 2000
*
*--

* Iraf global space
*
      real Memr(1)
      integer Memi(1)
      equivalence (Memr,Memi)
      common /Mem/Memr

* Local variables
      integer nx,ny,dims(2),id
      character*80 image
      real rad,sigma,xpos,ypos,k

* Get the name of the output image
      call uclgst('output',image,istat)

* and its dimensions
      call uclgsi('xdim',nx,istat)
      call uclgsi('ydim',ny,istat)

      dims(1)=nx
      dims(2)=ny

* and the size and fade-off parameters
      call uclgsr('radius',rad,istat)
      call uclgsr('sigma',sigma,istat)
      call uclgsr('k',k,istat)

* and the position for the centre of the mask
      call uclgsr('xpos',xpos,istat)
      call uclgsr('ypos',ypos,istat)

* Create the image
      call uimcre(image,6,2,dims,id,istat)
      if(istat.ne.0) then
         call umsput('Unable to create output image',
     :               1,0,istat)
         go to 99
      endif

* Allocate memory for the image
      call udmget(nx*ny,6,ipnt,istat)
      if(istat.ne.0) then
         call umsput('Unable to allocate memory for data',
     :               1,0,istat)
         go to 99
      endif

* Fill the image
      call circle(Memr(ipnt),nx,ny,rad,sigma,k,xpos,ypos)

* Write out the result
      call uips2r(id,1,nx,1,ny,Memr(ipnt),istat)

 99   continue
      call uimclo(id,istat)
      call udmfre(ipnt,6,istat)

      end

      subroutine circle(data,nx,ny,rad,sigma,k,xcen,ycen)
*
* Fill in a circular region with 1s and the rest zero except for
* a fade off around the circle with a specified sigma.
*
      implicit none

      integer nx,ny
      real data(nx,ny),rad,sigma,dx,dy,xcen,ycen,k,r
      integer i,j

      do j=1,ny
          dy=real(j)-ycen
          do i=1,nx
             dx=real(i)-xcen
             r=sqrt(dx*dx+dy*dy)
             if(r.le.rad) then
                data(i,j)=(1.0-(r/rad)**2)**(k-0.5)
             else if(r.le.rad+5*sigma) then
                data(i,j)=exp(-((r-rad)**2.0/(2.0*sigma*sigma)))
             else
                data(i,j)=0.0
             endif
          enddo
      enddo

      end

ext: f

max_stars_repo_name: iraf-community/stsdas


text:
      program main
        real x(3,3)
        data x/1, 0, 0, 0, 1, 0, 0, 0, 1/
        print *, x
      end program main

max_stars_repo_name: camfort/fortran-vars

max_issues_repo_name: camfort/fortran-vars

max_forks_repo_name: camfort/fortran-vars


text:
great(b1,n1).
great(l1,ll1).
great(kk1,w1).
great(ff1,i1).
great(a1,z1).
great(a1,kk1).
great(aa1,b1).
great(bb1,aa1).
great(d1,s1).
great(z1,i1).
great(r1,n1).
great(k1,cc1).
great(y1,w1).
great(ee1,aa1).
great(f1,w1).
great(g1,b1).
great(r1,x1).
great(ii1,n1).
great(aa1,o1).
great(hh1,z1).
great(u1,o1).
great(a1,g1).
great(ii1,z1).
great(g1,t1).
great(d1,l1).
great(u1,dd1).
great(j1,cc1).
great(z1,r1).
great(j1,ff1).
great(a1,aa1).
great(u1,p1).
great(d1,w1).
great(u1,aa1).
great(u1,s1).
great(d1,i1).
great(ff1,s1).
great(l1,b1).
great(k1,i1).
great(ii1,l1).
great(q1,r1).
great(l1,g1).
great(e1,d1).
great(ii1,ee1).
great(j1,q1).
great(d1,p1).
great(e1,n1).
great(h1,r1).
great(bb1,t1).
great(z1,dd1).
great(jj1,x1).
great(l1,bb1).
great(z1,j1).
great(ii1,q1).
great(dd1,w1).
great(dd1,t1).
great(d1,ff1).
great(f1,hh1).
great(q1,m1).
great(v1,m1).
great(x1,p1).
great(hh1,kk1).
great(c1,n1).
great(j1,s1).
great(p1,m1).
great(jj1,b1).
great(ii1,kk1).

max_stars_repo_name: manoelfranca/cilppp

max_issues_repo_name: manoelfranca/cilppp

max_forks_repo_name: manoelfranca/cilppp


! This file created from test/mpi/f77/rma/winattrf.f with f77tof90
! -*- Mode: Fortran; -*-
!
!  (C) 2003 by Argonne National Laboratory.
!      See COPYRIGHT in top-level directory.
!
      program main
      use mpi
      integer errs, ierr
      integer (kind=MPI_ADDRESS_KIND) extrastate, valin, valout, val

      integer comm, win, buf(10)
      integer curcount, keyval
      logical flag
      external mycopyfn, mydelfn
      integer callcount, delcount
      common /myattr/ callcount, delcount
!
! The only difference between the MPI-2 and MPI-1 attribute caching
! routines in Fortran is that the take an address-sized integer
! instead of a simple integer.  These still are not pointers,
! so the values are still just integers.
!
      errs      = 0
      callcount = 0
      delcount  = 0
      call mtest_init( ierr )
      call mpi_comm_dup( MPI_COMM_WORLD, comm, ierr )
! Create a new window; use val for an address-sized int
      val = 10
      call mpi_win_create( buf, val, 1, &
      &                        MPI_INFO_NULL, comm, win, ierr )
!
      extrastate = 1001
      call mpi_win_create_keyval( mycopyfn, mydelfn, keyval,  &
      &                             extrastate, ierr )
      flag = .true.
      call mpi_win_get_attr( win, keyval, valout, flag, ierr )
      if (flag) then
         errs = errs + 1
         print *, ' get attr returned true when no attr set'
      endif

      valin = 2003
      call mpi_win_set_attr( win, keyval, valin, ierr )
      flag = .false.
      valout = -1
      call mpi_win_get_attr( win, keyval, valout, flag, ierr )
      if (valout .ne. 2003) then
         errs = errs + 1
         print *, 'Unexpected value (should be 2003)', valout,  &
      &            ' from attr'
      endif

      valin = 2001
      call mpi_win_set_attr( win, keyval, valin, ierr )
      flag = .false.
      valout = -1
      call mpi_win_get_attr( win, keyval, valout, flag, ierr )
      if (valout .ne. 2001) then
         errs = errs + 1
         print *, 'Unexpected value (should be 2001)', valout,  &
      &            ' from attr'
      endif
!
! Test the attr delete function
      delcount   = 0
      call mpi_win_delete_attr( win, keyval, ierr )
      if (delcount .ne. 1) then
         errs = errs + 1
         print *, ' Delete_attr did not call delete function'
      endif
      flag = .true.
      call mpi_win_get_attr( win, keyval, valout, flag, ierr )
      if (flag) then
         errs = errs + 1
         print *, ' Delete_attr did not delete attribute'
      endif

! Test the delete function on window free
      valin = 2001
      call mpi_win_set_attr( win, keyval, valin, ierr )
      curcount = delcount
      call mpi_win_free( win, ierr )
      if (delcount .ne. curcount + 1) then
         errs = errs + 1
         print *, ' did not get expected value of delcount ',  &
      &          delcount, curcount + 1
      endif

      ierr = -1
      call mpi_win_free_keyval( keyval, ierr )
      if (ierr .ne. MPI_SUCCESS) then
         errs = errs + 1
         call mtestprinterror( ierr )
      endif
!
! The MPI standard defines null copy and duplicate functions.
! However, are only used when an object is duplicated.  Since
! MPI_Win objects cannot be duplicated, so under normal circumstances,
! these will not be called.  Since they are defined, they should behave
! as defined.  To test them, we simply call them here
      flag   = .false.
      valin  = 7001
      valout = -1
      ierr   = -1
      call MPI_WIN_DUP_FN( win, keyval, extrastate, valin, valout, &
      &     flag, ierr )
      if (.not. flag) then
         errs = errs + 1
         print *, " Flag was false after MPI_WIN_DUP_FN"
      else if (valout .ne. 7001) then
         errs = errs + 1
         if (valout .eq. -1 ) then
          print *, " output attr value was not copied in MPI_WIN_DUP_FN"
         endif
         print *, " value was ", valout, " but expected 7001"
      else if (ierr .ne. MPI_SUCCESS) then
         errs = errs + 1
         print *, " MPI_WIN_DUP_FN did not return MPI_SUCCESS"
      endif

      flag   = .true.
      valin  = 7001
      valout = -1
      ierr   = -1
      call MPI_WIN_NULL_COPY_FN( win, keyval, extrastate, valin, valout &
      &     ,flag, ierr )
      if (flag) then
         errs = errs + 1
         print *, " Flag was true after MPI_WIN_NULL_COPY_FN"
      else if (valout .ne. -1) then
         errs = errs + 1
         print *, &
      &        " output attr value was copied in MPI_WIN_NULL_COPY_FN"
      else if (ierr .ne. MPI_SUCCESS) then
         errs = errs + 1
         print *, " MPI_WIN_NULL_COPY_FN did not return MPI_SUCCESS"
      endif
!
      call mpi_comm_free( comm, ierr )
      call mtest_finalize( errs )
      call mpi_finalize( ierr )
      end
!
! Note that the copyfn is unused for MPI windows, since there is
! (and because of alias rules, can be) no MPI_Win_dup function
      subroutine mycopyfn( oldwin, keyval, extrastate, valin, valout, &
      &                     flag, ierr )
      use mpi
      integer oldwin, keyval, ierr
      integer (kind=MPI_ADDRESS_KIND) extrastate, valin, valout, val

      logical flag
      integer callcount, delcount
      common /myattr/ callcount, delcount
! increment the attribute by 2
      valout = valin + 2
      callcount = callcount + 1
!
! Since we should *never* call this, indicate an error
      print *, ' Unexpected use of mycopyfn'
      flag = .false.
      ierr = MPI_ERR_OTHER
      end
!
      subroutine mydelfn( win, keyval, val, extrastate, ierr )
      use mpi
      integer win, keyval, ierr
      integer (kind=MPI_ADDRESS_KIND) extrastate, valin, valout, val

      integer callcount, delcount
      common /myattr/ callcount, delcount
      delcount = delcount + 1
      if (extrastate .eq. 1001) then
         ierr = MPI_SUCCESS
      else
         print *, ' Unexpected value of extrastate = ', extrastate
         ierr = MPI_ERR_OTHER
      endif
      end


text:
program main
  implicit none
  integer, parameter :: max_greeting_length=64
  integer image
  character(len=max_greeting_length), allocatable :: greeting(:)[:]
  allocate(greeting(num_images())[*])
  write(greeting(this_image()),*) "Hello from image ",this_image()," of ", num_images()
  greeting(this_image())[1] = greeting (this_image())
  sync all
  if (this_image()==1) then
    do concurrent(image=1:num_images())
      print *,greeting(image)
    end do
  end if
end program

max_stars_repo_name: sourceryinstitute/miniFAVOR

max_issues_repo_name: sourceryinstitute/miniFAVOR

max_forks_repo_name: sourceryinstitute/miniFAVOR


text:
C
C $Id: slgetr.f,v 1.7 2008-07-27 00:17:26 haley Exp $
C
C                Copyright (C)  2000
C        University Corporation for Atmospheric Research
C                All Rights Reserved
C
C The use of this Software is governed by a License Agreement.
C
      SUBROUTINE SLGETR (PNAM,RVAL)
C
C Get in RVAL the real value of the STITLE parameter named PNAM.
C
        CHARACTER*(*) PNAM
C
C The common block SLCOMN holds all of the internal parameters of
C the package STITLE except for color-table parameters.
C
        COMMON /SLCOMN/ GPSZ,IBGC,IBGF,ICOP,IDOT,IFGC,IFGF,IJMP,IMAP,
     +                  INCU,IWLU,IWRK,IWWI,IXND,IXST,OORV,PCSZ,RNFS,
     +                  RVPB,RVPL,RVPR,RVPT,TFIN,TFOU,TGP1,TGP2,TGP3
        SAVE   /SLCOMN/
C
C Define a temporary variable in which to put the first three characters
C of PNAM.
C
        CHARACTER*3 CTMP
C
C Define a character variable in which messages may be formed.
C
        CHARACTER*29 CMSG
C
C Do a call forcing a BLOCKDATA to be loaded from a binary library.
C
        CALL SLBLDA
C
C Check for an uncleared prior error.
C
        IF (ICFELL('SLGETR - UNCLEARED PRIOR ERROR',1).NE.0) RETURN
C
C Extract the first three characters of the parameter name.
C
        CTMP=PNAM
C
C If the parameter name has less than three characters, log an error.
C
        IF (LEN(PNAM).LT.3) GO TO 901
C
C See what the parameter name is ...
C
C ... the flag controlling the output of the alignment frames, ...
C
        IF      (CTMP.EQ.'ALN'.OR.CTMP.EQ.'aln') THEN
C
          RVAL=REAL(IDOT)
C
C ... the blue component of the background color, ...
C
        ELSE IF (CTMP.EQ.'BGB'.OR.CTMP.EQ.'bgb') THEN
C
          CALL SLGCLR (IBGC,DUMI,DUMI,RVAL)
          IF (ICFELL('SLGETR',2).NE.0) RETURN
C
C ... the default background color index ...
C
        ELSE IF (CTMP.EQ.'BGC'.OR.CTMP.EQ.'bgc') THEN
C
          RVAL=REAL(IBGC)
C
C ... the background color fade flag ...
C
        ELSE IF (CTMP.EQ.'BGF'.OR.CTMP.EQ.'bgf') THEN
C
          RVAL=REAL(IBGF)
C
C ... the green component of the background color, ...
C
        ELSE IF (CTMP.EQ.'BGG'.OR.CTMP.EQ.'bgg') THEN
C
          CALL SLGCLR (IBGC,DUMI,RVAL,DUMI)
          IF (ICFELL('SLGETR',3).NE.0) RETURN
C
C ... the red component of the background color, ...
C
        ELSE IF (CTMP.EQ.'BGR'.OR.CTMP.EQ.'bgr') THEN
C
          CALL SLGCLR (IBGC,RVAL,DUMI,DUMI)
          IF (ICFELL('SLGETR',4).NE.0) RETURN
C
C ... the blue component of the foreground color, ...
C
        ELSE IF (CTMP.EQ.'FGB'.OR.CTMP.EQ.'fgb') THEN
C
          CALL SLGPAI (PNAM,4,IPAI)
          IF (IPAI.LE.0) IPAI=IFGC
          CALL SLGCLR (IPAI,DUMI,DUMI,RVAL)
          IF (ICFELL('SLGETR',5).NE.0) RETURN
C
C ... the default foreground color index ...
C
        ELSE IF (CTMP.EQ.'FGC'.OR.CTMP.EQ.'fgc') THEN
C
          RVAL=REAL(IFGC)
C
C ... the foreground color fade flag ...
C
        ELSE IF (CTMP.EQ.'FGF'.OR.CTMP.EQ.'fgf') THEN
C
          RVAL=REAL(IFGF)
C
C ... the green component of the foreground color, ...
C
        ELSE IF (CTMP.EQ.'FGG'.OR.CTMP.EQ.'fgg') THEN
C
          CALL SLGPAI (PNAM,4,IPAI)
          IF (IPAI.LE.0) IPAI=IFGC
          CALL SLGCLR (IPAI,DUMI,RVAL,DUMI)
          IF (ICFELL('SLGETR',6).NE.0) RETURN
C
C ... the red component of the foreground color, ...
C
        ELSE IF (CTMP.EQ.'FGR'.OR.CTMP.EQ.'fgr') THEN
C
          CALL SLGPAI (PNAM,4,IPAI)
          IF (IPAI.LE.0) IPAI=IFGC
          CALL SLGCLR (IPAI,RVAL,DUMI,DUMI)
          IF (ICFELL('SLGETR',7).NE.0) RETURN
C
C ... the fade-in time, ...
C
        ELSE IF (CTMP.EQ.'FIN'.OR.CTMP.EQ.'fin') THEN
C
          RVAL=TFIN
C
C ... the fade-out time, ...
C
        ELSE IF (CTMP.EQ.'FOU'.OR.CTMP.EQ.'fou') THEN
C
          RVAL=TFOU
C
C ... the interline spacing ("gap size"), ...
C
        ELSE IF (CTMP.EQ.'GSZ'.OR.CTMP.EQ.'gsz') THEN
C
          RVAL=GPSZ
C
C ... the centering parameter, ...
C
        ELSE IF (CTMP.EQ.'ICO'.OR.CTMP.EQ.'ico') THEN
C
          RVAL=REAL(ICOP)
C
C ... the FORTRAN logical unit number for "card" input, ...
C
        ELSE IF (CTMP.EQ.'ICU'.OR.CTMP.EQ.'icu') THEN
C
          RVAL=REAL(INCU)
C
C ... the interline spacing for practice runs, ...
C
        ELSE IF (CTMP.EQ.'INC'.OR.CTMP.EQ.'inc') THEN
C
          RVAL=REAL(IJMP)
C
C ... the FORTRAN logical unit number for WISS, ...
C
        ELSE IF (CTMP.EQ.'LOG'.OR.CTMP.EQ.'log') THEN
C
          RVAL=REAL(IWLU)
C
C ... the plotter X coordinate of the left edge of the viewport, ...
C
        ELSE IF (CTMP.EQ.'LX1'.OR.CTMP.EQ.'lx1') THEN
C
          RVAL=32767.*RVPL
C
C ... the plotter X coordinate of the right edge of the viewport, ...
C
        ELSE IF (CTMP.EQ.'LX2'.OR.CTMP.EQ.'lx2') THEN
C
          RVAL=32767.*RVPR
C
C ... the plotter Y coordinate of the bottom edge of the viewport, ...
C
        ELSE IF (CTMP.EQ.'LY1'.OR.CTMP.EQ.'ly1') THEN
C
          RVAL=32767.*RVPB
C
C ... the plotter Y coordinate of the top edge of the viewport, ...
C
        ELSE IF (CTMP.EQ.'LY2'.OR.CTMP.EQ.'ly2') THEN
C
          RVAL=32767.*RVPT
C
C ... the PLOTCHAR mapping flag, ...
C
        ELSE IF (CTMP.EQ.'MAP'.OR.CTMP.EQ.'map') THEN
C
          RVAL=REAL(IMAP)
C
C ... the number of frames per second, ...
C
        ELSE IF (CTMP.EQ.'NFS'.OR.CTMP.EQ.'nfs') THEN
C
          RVAL=RNFS
C
C ... the horizontal scroll end coordinate, ...
C
        ELSE IF (CTMP.EQ.'NXE'.OR.CTMP.EQ.'nxe') THEN
C
          RVAL=REAL(IXND)
C
C ... the horizontal scroll start coordinate, ...
C
        ELSE IF (CTMP.EQ.'NXS'.OR.CTMP.EQ.'nxs') THEN
C
          RVAL=REAL(IXST)
C
C ... the PLOTCHAR out-of-range value, ...
C
        ELSE IF (CTMP.EQ.'ORV'.OR.CTMP.EQ.'orv') THEN
C
          RVAL=OORV
C
C ... the character size, ...
C
        ELSE IF (CTMP.EQ.'PSZ'.OR.CTMP.EQ.'psz') THEN
C
          RVAL=PCSZ
C
C ... the flag controlling suppression of background color during a
C fade-in or fade-out, ...
C
        ELSE IF (CTMP.EQ.'SBK'.OR.CTMP.EQ.'sbk') THEN
C
          IF (IBGF.NE.-1) THEN
            RVAL=0.
          ELSE
            RVAL=1.
          END IF
C
C ... the flag controlling suppression of foreground color during a
C fade-in or fade-out, ...
C
        ELSE IF (CTMP.EQ.'SFG'.OR.CTMP.EQ.'sfg') THEN
C
          IF (IFGF.NE.-1) THEN
            RVAL=0.
          ELSE
            RVAL=1.
          END IF
C
C ... the initial blank-frame count (for FTITLE), ...
C
        ELSE IF (CTMP.EQ.'TM1'.OR.CTMP.EQ.'tm1') THEN
C
          RVAL=TGP1
C
C ... the intermediate blank-frame count (for FTITLE), ...
C
        ELSE IF (CTMP.EQ.'TM2'.OR.CTMP.EQ.'tm2') THEN
C
          RVAL=TGP2
C
C ... the terminal extra blank-frame count (for FTITLE), ...
C
        ELSE IF (CTMP.EQ.'TM3'.OR.CTMP.EQ.'tm3') THEN
C
          RVAL=TGP2
C
C ... the NDC Y coordinate of the bottom edge of the viewport, ...
C
        ELSE IF (CTMP.EQ.'VPB'.OR.CTMP.EQ.'vpb') THEN
C
          RVAL=RVPB
C
C ... the NDC X coordinate of the left edge of the viewport, ...
C
        ELSE IF (CTMP.EQ.'VPL'.OR.CTMP.EQ.'vpl') THEN
C
          RVAL=RVPL
C
C ... the NDC X coordinate of the right edge of the viewport, ...
C
        ELSE IF (CTMP.EQ.'VPR'.OR.CTMP.EQ.'vpr') THEN
C
          RVAL=RVPR
C
C ... the NDC Y coordinate of the top edge of the viewport, ...
C
        ELSE IF (CTMP.EQ.'VPT'.OR.CTMP.EQ.'vpt') THEN
C
          RVAL=RVPT
C
C ... or the workstation identifier for WISS.
C
        ELSE IF (CTMP.EQ.'WID'.OR.CTMP.EQ.'wid') THEN
C
          RVAL=REAL(IWWI)
C
        ELSE
C
C Otherwise, the parameter name is not recognized.
C
          GO TO 901
C
        END IF
C
        RETURN
C
C Error exit.
C
  901   CMSG(1:29)='SLGETR - INVALID KEYWORD: '//CTMP
        CALL SETER (CMSG(1:29),8,1)
        RETURN
C
      END

max_stars_repo_name: tenomoto/ncl

max_issues_repo_name: tenomoto/ncl

max_forks_repo_name: tenomoto/ncl


text:
c
c -------------------------------------------------------------
c
      subroutine putnod (mptr)
c
      implicit double precision (a-h,o-z)

      include  "call.i"

c :::::::::::::::::::::::::::::: PUTNOD :::::::::::::::::::::;
c
c  return mptr node to the linked list kept in node array.
c
c ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::;
c
      node(nextfree, mptr) = ndfree
      ndfree               = mptr
c
      return
      end

max_issues_repo_name: che-wenchao/D-Claw

max_forks_repo_name: che-wenchao/D-Claw


text:
C Example of Fortran II program:
C This program has one data card input, with simple zero-data edit check.
C If one of the input values is zero, then the program will end with an
C error code of "STOP 1" in the job control card listing following the
C execution of the program. Normal output will be one line printed with
C A, B, C, and AREA on the "report" page following the compilation with
C no specific units are stated; and with a normal "STOP" in the job card
C listing.

C AREA OF A TRIANGLE WITH A STANDARD SQUARE ROOT FUNCTION
C INPUT - CARD READER UNIT 5, INTEGER INPUT
C OUTPUT - LINE PRINTER UNIT 6, REAL OUTPUT
C INPUT ERROR DISPAY ERROR OUTPUT CODE 1 IN JOB CONTROL LISTING
      INTEGER A,B,C
      READ(5,501) A,B,C
  501 FORMAT(3I5)
      IF(A.EQ.0 .OR. B.EQ.0 .OR. C.EQ.0) STOP 1
      S = (A + B + C) / 2.0
      AREA = SQRT( S * (S - A) * (S - B) * (S - C) / 2.0)
      WRITE(6,601) A,B,C,AREA
  601 FORMAT(4H A=,I5,5H B=,I5,5H C= ,I5,8H AREA=,F10.2,12HSQUARE)
      STOP
      END

max_stars_repo_name: maurizioabba/rose

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


text:
      subroutine tmpadj(hour,halfdy)
c
c     +++PURPOSE+++
c     This function adjusts the surface temperature on an hourly
c     basis by using the daily minumum and maximum temperatures.
c     It also uses a similar method for computing the hourly
c     dewpoint based on the daily average.
c
c     Authors(s):  John Witte, UofMn WCES @ USDA-ARS-NCSRL
c     Date: 04/01/93

c     Verified and tested by Reza Savabi, USDA-ARS, NSERL 317-494-5051
c                  August 1994
cd    The program is modiffied by S. Dun Aug 14, 2008 to direcly use SI units
c
c     +++ARGUMENT DECLARATIONS+++
      real    halfdy
      integer hour
c
c     +++ARGUMENT DEFINITIONS+++
c     halfdy - The amount of time from sunrise until noon (hrs).
c     hour   - The hour of the day that we are calculating.
c
c     +++PARAMETERS+++
      include 'pmxpln.inc'
      include 'pmxnsl.inc'
      include 'pmxhil.inc'
      include 'pmxtls.inc'
      include 'pmxtil.inc'
      include 'pmxcrp.inc'
      include 'pmxelm.inc'
c
c     +++COMMON BLOCKS+++
      include  'cclim.inc'
c       read:  tave,vwind,hradmj,hrtemp,radmj
c
      include  'cwint.inc'
c       read:  snodpy(mxplane),tfrdp(mxplan),tthawd(mxplan),surtmp
c              frdp(mxplan),thdp(mxplan),densg(mxplan)
c
      include  'cwater.inc'
c       read:  salb(mxplan),
c
      include  'ccover.inc'
c       read:  canhgt(mxplan)
c
      include  'ccrpout.inc'
c       read:  rrc(mxplan)
c
      include  'cstruc.inc'
c       read:  iplane
cd      Added by S. Dun, Nov. 21, 2006
      include 'cangie.inc'
c      read : solef (solar energy adjust fator for sloped surface)
cd      End adding
c
c     +++LOCAL VARIABLES+++
c     save
      real    alb,rads,lytomj,clouds,ktemp,aemiss,displ,wsrgh,etrgh,
     1        convht,vkcons,denair,hcpair,htwind,radcof,semiss,sbcons,
     2        netrad,grdp,gutdp,sysdep,ksnow,kres,kftill,kfutil,
     3        numer,denom,surfcn,effk,bakup, gtdp, hrtave,factor
c    2        netrad,radadj,grdp,gutdp,sysdep,ksnow,kres,kftill,kfutil,
c
c     +++LOCAL DEFINITIONS+++
c     alb    - Albedo of the top layer for the system (Dec %).
c     rads   - Radiation on a sloping surface (W/m^2).
c     lytomj - Conversn factor of solar radiation (Langly's to MJ/m^2).
c     clouds - Estimated cloud cover (Dec %).
c     ktemp  - Average temperature in degrees Kelvin units.
c     aemiss - Atmospheric emissivity.
c     displ  - Zero-plane displacment, a function of canopy cover (m).
c     wsrgh  - Roughness parameter for wind speed (m).
c     etrgh  - Roughness parameter for energy transfer (m).
c     convht - Convective heat transfer coeff (J/m^3-K).
c     vkcons - Von Karman constant.
c     denair - Approximate density of air (Kg/m^3).
c     hcpair - Heat capacity of air (J/Kg-K).
c     htwind - Height at which wind velocity is measured (m).
c     radcof - Long_wave radiation transfer coefficient (W/m^2-K).
c     semiss - Surface emissivity for all layers.
c     sbcons - Stefan-Bolzman constant (W/m^2-K^4).
c     netrad - Calculated net radiation (W/m^2).
c     radadj - Daily radiation adjusted to hourly (MJ/m^2).
c     grdp   - Depth from soil's surface to first 0 degree isotherm (m).
c     gutdp  - Temperature gradient depth of untilled layer (m).
c     sysdep - Depth of the frozen layer system from surface to
c              first 0 degree isotherm (m).
c     ksnow  -   These are temporary variables representing the thermal
c     kres   --_ conductivity of the snow, residue, frozen tilled and
c     kftill --  frozen untilled layers.  Reason for temps is because
c     kfutil -   if no such layer exists, the values are set to zero.
c     numer  - Temporary numerator.
c     denom  - Temporary denominator.
c     surfcn - Surface conductivity of top layers (W/m-C).
c     effk   - Effective thermal conductivity of the entire layered
c              system to 0 degree isotherm (W/m-C).
c     bakup  - Temporary variable used in simplifying equations.
c      hrtave - hourly air temperature to aviod messing up with daily average temperature stored
c
c     +++DATA INITIALIZATIONS+++
c     data   sbcons/8.1247e-11/
c     data   vkcons/0.4/
c     data   hcpair/101.0/
c     data   denair/1.0/
c     data   htwind/2.0/
c     data   lytomj/0.04184/
c     data   kres/0.0232/
c
c     +++END SPECIFICATIONS+++
cd      sbcons = 8.1247e-11
      sbcons = 5.6697e-8
      vkcons = 0.4
      hcpair = 1012.0
      denair = 1.2
      htwind = 2.0
      lytomj = 0.04184
c
c XXX  WHY is the value for KRES here set to 0.0232, while in
c      subroutine FROST, kres = 0.168 ??????   dcf  2/23/94
c
cd      kres = 0.0232
        kres = 0.05
c       kres = 0.168
        kres = kres*kresf

      grdp = 0.0
      gtdp = 0.0
      gutdp = 0.0

      alb = salb(iplane)
      if (snodpy(iplane) .gt. 0.01) alb = 0.5

c -- Convert from MJ/m^2/hr to L/min.

cd      rads = (hradmj / 60) * (1 / lytomj)
c -- Convert from MJ/m^2/hr to W/m^2
      rads = (hradmj / 3600) *1E+6

c      Modified by S. Dun, Jan 12, 2007
      clouds = CloudC
c
cd      if (rpoth .lt. 0.0001) then
cd        clouds = 1.0
cd      else
cd        clouds = (1.0 - radmj /rpoth) / 0.7
cd      endif

cd      if (clouds .gt. 0.999) then
cd        clouds = 1.0
cd      endif
cd      End modifying

      call hrtmp(hour,halfdy)

      hrtave = hrtemp
      ktemp = hrtave + 273.16

      semiss = 1.0
      aemiss = (1.0 - 0.84 * clouds) * (1.0 - 0.261 *
     1         exp(7.77e-4  * (hrtave * hrtave))) + (0.84 * clouds)

      displ = 0.77 * canhgt(iplane)
c reza added next line 7/27/94
c      if(displ.gt.htwind)displ=0.77*htwind
c reza added the above line 7/27/94
c fix from John Laflen to prevent divide by zero below
       if(displ.ge.htwind)displ=0.77*htwind
c fix from John Laflen above - 10/2/2001
      if((snodpy(iplane) .lt. 0.01).and.(canhgt(iplane) .gt. 0.0)) then
        wsrgh = 0.13 * canhgt(iplane)
      elseif (snodpy(iplane) .lt. 0.01) then
        wsrgh = rrc(iplane)
      elseif (snodpy(iplane) .gt. canhgt(iplane)) then
        wsrgh = 0.0002
      else
        wsrgh = 0.13 * (canhgt(iplane) - snodpy(iplane))
      endif

      if (wsrgh .lt. 0.001) wsrgh = 0.001
c Reza added 7/27/94
      if(wsrgh.gt..26)wsrgh=0.26
c end Reza 7/27/94
      etrgh = 0.2 * wsrgh

      convht =((vkcons * vkcons) * denair * hcpair) /
     1         (log((htwind - displ + etrgh) / etrgh) *
     2          log((htwind - displ + wsrgh) / wsrgh))

      radcof = 4.0 * semiss * sbcons * (ktemp**3)
c
c XXX  WHY is RADADJ set equal to RADCOF here, then never used???????
c      Is the correct variable being used below (RADCOF) or should
c      RADADJ be used???  Can RADADJ be deleted????   dcf   2/22/94
c
c     radadj = radcof

cd      netrad = (1.0 - alb)*rads + (aemiss - semiss)*sbcons*hrtave**4
      netrad = (1.0 - alb)*rads + (aemiss - semiss)*sbcons*ktemp**4

c -- Now we calculate the gradient depth for both top frost and none...

c -- If average temperature is below zero.

      if (hrtave .lt. 0.0) then

c -- If top frost is present...

        if (tfrdp(iplane) .gt. 0.001) then
          if (thdp(iplane) .gt. 0.001) then
            grdp = 0.0
          else
            grdp = tfrdp(iplane)
          endif
        else

c -- If no top frost is present...

          if (frdp(iplane) .gt. 0.001) then
            if (thdp(iplane) .gt. 0.001) then
              grdp = 0.0
            else
              grdp = frdp(iplane)
            endif
          endif
        endif

c -- The above assumption comes directly out of Flershinger's notes
c -- and assumes frost exists all the way from the frost depth to
c -- the surface.

c -- Now, we consider the gradient depth is avg temperature is above 0.

      else

c -- If top frost is present...

        if (tfrdp(iplane) .gt. 0.001) then
            grdp = thdp(iplane)
        else

c -- If no top frost is present...

          if (frdp(iplane) .gt. 0.001) then
            if (thdp(iplane) .gt. 0.001) then
              grdp = thdp(iplane)
            else
              grdp = 0.0
            endif
          endif
        endif
      endif

      if ((grdp .le. tilld(iplane)) .or. (tilld(iplane) .gt. 1.0)) then
        gtdp = grdp
        gutdp = 0.0
      else
        gtdp = tilld(iplane)
        gutdp = grdp - tilld(iplane)
      endif

      if (gtdp .lt. 0.001)  gtdp = 0.0
      if (gutdp .lt. 0.001) gutdp = 0.0
      if (grdp .lt. 0.001)  grdp = 0.0

      sysdep = snodpy(iplane) + resdep(iplane) + grdp

c -- Now we make sure that each layer exists, if not we simply
c -- set the "tc" value of that layer to 1.0 so that multiplying
c -- by the "tc" will leave the final value unaffected by that
c -- absent layer.

      if (snodpy(iplane) .lt. 0.0001) then
        ksnow = 1.0
      else
c         Considering thermal conductivity of the snow pack
c         model from Sturm et al., 1997
          if (densg(iplane).lt. 156) then
                ksnow = 0.023 + 0.234 * (densg(iplane)/1000.)
          else
                ksnow = 0.138 - 1.01*(densg(iplane)/1000.)
     1                  + 3.233* (densg(iplane)/1000.)** 2
          endif
cd          ksnow = 0.046 + 0.0451 * (densg(iplane)/1000) + 0.3447 *
cd     1            ((densg(iplane)/1000) * (densg(iplane)/1000))
c
          ksnow = ksnow * ksnowf
      endif

      if (resdep(iplane) .lt. 0.0001) kres = 1.0
c
c XXX IS the following right???  Values for KFTILL when gtdp.ge.0.0001
c     were not defined.  I got the value of 1.75 from subroutine FROST
c     REZA ???
      if (gtdp .lt. 0.0001) then
        kftill = 1.0
      else
        kftill = 1.75
      endif
c
c XXX IS the following right???  Values for KFUTIL when gutdp.ge.0.0001
c     were not defined.  I got the value of 2.1 from subroutine FROST
c     REZA ???
      if (gutdp .lt. 0.0001) then
        kfutil = 1.0
      else
        kfutil = 2.1
      endif
cd   Added by S. Dun, May 25, 2009
c    When there is no frost layers in the soil, we assume 0c isotherm is at 0.1 mm
      if (sysdep .lt. 0.0001) then
         kftill = 1.75
         grdp = 0.001
         gtdp = 0.001
         sysdep = snodpy(iplane) + resdep(iplane) + grdp
      endif
cd   End adding

c -- Now we use the harmonic mean method to calculate the surface
c -- temperature for the hour.  NOTE that the final calculation of
c -- numerator/denominator has been broken up into two parts... the
c -- "numer" and the "denom".

      numer = (ksnow * kres * kftill * kfutil) *
     1         (snodpy(iplane) + resdep(iplane) + grdp)

      denom = (ksnow * kres * kftill * gutdp) +
     1         (ksnow * kres * gtdp * kfutil) +
     2          (ksnow * resdep(iplane) * kftill * kfutil) +
     3           (snodpy(iplane) * kres * kftill * kfutil)

c -- Now we must make sure we have no floating point errors...

      if ((denom .lt. -0.0001) .or. (denom .gt. 0.0001)) then
        surfcn = numer / denom

      else
        surfcn = 0.0
      endif

c -- If the second case is chosen, an error has occurred!

c -- The following line converts from MJoules to Langly's

cd      effk = surfcn / 0.0697333
c       Directly use the SI units
        effk = surfcn

      if (sysdep .gt. 0.0) then
cd         surtmp(hour) = (netrad+(radcof+convht*(vwind*100))*hrtave+0.0)/
cd     1                   (radcof+convht*(vwind*100)+(effk/sysdep))
         surtmp(hour) = (netrad+(radcof+convht*(vwind))*hrtave)/
     1                   (radcof+convht*(vwind)+(effk/sysdep))
      else
cd         surtmp(hour) = (netrad+(radcof+convht*(vwind*100))*hrtave+0.0)/
cd     1                  (radcof+convht*(vwind*100))
         surtmp(hour) = (netrad+(radcof+convht*(vwind))*hrtave)/
     1                  (radcof+convht*(vwind))
      endif
cd   Added by S. Dun, May 25, 2009
c    When there is snow, we assume surface temperture is not greater than 0C
      if ((surtmp(hour).gt. 0.0).and.(snodpy(iplane) .gt. 0.001)) then
          surtmp(hour) = 0.0
      endif
cd   End adding
      return
      end

ext: for

max_stars_repo_name: jarad/dep

max_issues_repo_name: jarad/dep

max_forks_repo_name: akrherz/idep


text: #redirect Governors Ball

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


! { dg-do compile }
! { dg-options "-fimplicit-none" }
!
! PR 45748: [4.5/4.6 Regression] -fimplicit-none failures when using intrinsic MAX
!
! Contributed by Themos Tsikas <themos.tsikas@gmail.com>

SUBROUTINE BUG(WORK)
  INTRINSIC MAX
  DOUBLE PRECISION WORK(MAX(2,3))
END


text:
      SUBROUTINE DIJKL1 (C,N,NATI,W,CIJ,WCIJ,CKL)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      INCLUDE 'SIZES'
      DIMENSION C(N,*), W(*)
      DIMENSION CIJ(10*MAXORB), WCIJ(10*MAXORB), CKL(10*MAXORB)
************************************************************************
*
*   DIJKL1 IS SIMILAR TO IJKL.  THE MAIN DIFFERENCES ARE THAT
*   THE ARRAY W CONTAINS THE TWO ELECTRON INTEGRALS BETWEEN
*   ONE ATOM (NATI) AND ALL THE OTHER ATOMS IN THE SYSTEM.
*
*           ON EXIT
*
*   THE ARRAY XY IS FILLED WITH THE DIFFERENTIALS OF THE
*   TWO-ELECTRON INTEGRALS OVER ACTIVE-SPACE M.O.S W.R.T. MOTION
*   OF THE ATOM NATI.
************************************************************************
      COMMON /MOLKST/ NUMAT,NAT(NUMATM),NFIRST(NUMATM),NMIDLE(NUMATM),
     1                NLAST(NUMATM), NORBS, NELECS,NALPHA,NBETA,
     2                NCLOSE,NOPEN,NDUMY,FRACT
      COMMON /CIBITS/ NMOS,LAB,NELEC, NBO(3)
      COMMON /XYIJKL/ XY(NMECI,NMECI,NMECI,NMECI)
      DIMENSION NB(0:8)
      DATA NB /1,0,0,10,0,0,0,0,45/
      NA=NMOS
      DO 110 I=1,NA
         DO 110 J=1,I
            IPQ=0
            DO 20 II=1,NUMAT
               IF(II.EQ.NATI) GOTO 20
               DO 10 IP=NFIRST(II),NLAST(II)
                  DO 10 IQ=NFIRST(II),IP
                     IPQ=IPQ+1
                     CIJ(IPQ)=C(IP,I)*C(IQ,J)+C(IP,J)*C(IQ,I)
   10          CONTINUE
   20       CONTINUE
            I77=IPQ+1
            DO 30 IP=NFIRST(NATI),NLAST(NATI)
               DO 30 IQ=NFIRST(NATI),IP
                  IPQ=IPQ+1
                  CIJ(IPQ)=C(IP,I)*C(IQ,J)+C(IP,J)*C(IQ,I)
   30       CONTINUE
            DO 40 II=1,IPQ
   40       WCIJ(II)=0.D0
            KR=1
            JS=1
            NBJ=NB(NLAST(NATI)-NFIRST(NATI))
            DO 50 II=1,NUMAT
               IF (II.EQ.NATI) GOTO 50
               NBI=NB(NLAST(II)-NFIRST(II))
               CALL FORMXY
     1(W(KR), KR, WCIJ(I77), WCIJ(JS), CIJ(I77), NBJ, CIJ(JS), NBI)
               JS=JS+NBI
   50       CONTINUE
            DO 100 K=1,I
               IF(K.EQ.I) THEN
                  LL=J
               ELSE
                  LL=K
               ENDIF
               DO 100 L=1,LL
                  IPQ=0
                  DO 70 II=1,NUMAT
                     IF(II.EQ.NATI) GOTO 70
                     DO 60 IP=NFIRST(II),NLAST(II)
                        DO 60 IQ=NFIRST(II),IP
                           IPQ=IPQ+1
                           CKL(IPQ)=C(IP,K)*C(IQ,L)+C(IP,L)*C(IQ,K)
   60                CONTINUE
   70             CONTINUE
                  DO 80 IP=NFIRST(NATI),NLAST(NATI)
                     DO 80 IQ=NFIRST(NATI),IP
                        IPQ=IPQ+1
                        CKL(IPQ)=C(IP,K)*C(IQ,L)+C(IP,L)*C(IQ,K)
   80             CONTINUE
                  SUM=0.D0
                  DO 90 II=1,IPQ
   90             SUM=SUM+CKL(II)*WCIJ(II)
                  XY(I,J,K,L)=SUM
                  XY(I,J,L,K)=SUM
                  XY(J,I,K,L)=SUM
                  XY(J,I,L,K)=SUM
                  XY(K,L,I,J)=SUM
                  XY(K,L,J,I)=SUM
                  XY(L,K,I,J)=SUM
                  XY(L,K,J,I)=SUM
  100       CONTINUE
  110 CONTINUE
      RETURN
      END

max_stars_repo_name: openmopac/MOPAC-archive

max_issues_repo_name: openmopac/MOPAC-archive

max_forks_repo_name: openmopac/MOPAC-archive