text:
      FUNCTION expdev(idum)
      INTEGER idum
      REAL expdev
CU    USES ran1
      REAL dum,ran1
1     dum=ran1(idum)
      if(dum.eq.0.)goto 1
      expdev=-log(dum)
      return
      END

ext: for

max_stars_repo_path: NR-Functions/Numerical Recipes- Example & Functions/Functions/expdev.for


text:
!
!  Written by Leandro Martínez, 2009-2011.
!  Copyright (c) 2009-2018, Leandro Martínez, Jose Mario Martinez,
!  Ernesto G. Birgin.
!
!
! Subroutine comprest: Compute the function value relative to
!                      to the restrictions for one atom
!

subroutine comprest(icart,f)

  use sizes
  use compute_data, only : xcart, restpars, scale, scale2, nratom, ityperest, iratom

  implicit none
  integer :: iratcount, irest, icart
  double precision :: xmin, ymin, zmin, clength, a1, a2, a3, a4, w, b1, b2, b3, d, a5, a6
  double precision :: f
  double precision :: xmax, ymax, zmax
  double precision :: v1, v2, v3
  double precision :: vnorm

  f = 0.d0
  do iratcount = 1, nratom(icart)
    irest = iratom(icart,iratcount)
    if(ityperest(irest).eq.2) then
      clength = restpars(irest,4)
      xmin = restpars(irest,1)
      ymin = restpars(irest,2)
      zmin = restpars(irest,3)
      xmax = restpars(irest,1) + clength
      ymax = restpars(irest,2) + clength
      zmax = restpars(irest,3) + clength
      a1 = dmin1(xcart(icart,1) - xmin, 0.d0)
      a2 = dmin1(xcart(icart,2) - ymin, 0.d0)
      a3 = dmin1(xcart(icart,3) - zmin, 0.d0)
      f = f + scale*(a1 * a1 + a2 * a2 + a3 * a3)
      a1 = dmax1(xcart(icart,1) - xmax, 0.d0)
      a2 = dmax1(xcart(icart,2) - ymax, 0.d0)
      a3 = dmax1(xcart(icart,3) - zmax, 0.d0)
      f = f + scale*(a1 * a1 + a2 * a2 + a3 * a3)
    else if(ityperest(irest).eq.3) then
      xmin = restpars(irest,1)
      ymin = restpars(irest,2)
      zmin = restpars(irest,3)
      xmax = restpars(irest,4)
      ymax = restpars(irest,5)
      zmax = restpars(irest,6)
      a1 = dmin1(xcart(icart,1) - xmin, 0.d0)
      a2 = dmin1(xcart(icart,2) - ymin, 0.d0)
      a3 = dmin1(xcart(icart,3) - zmin, 0.d0)
      f = f + scale*(a1 * a1 + a2 * a2 + a3 * a3)
      a1 = dmax1(xcart(icart,1) - xmax, 0.d0)
      a2 = dmax1(xcart(icart,2) - ymax, 0.d0)
      a3 = dmax1(xcart(icart,3) - zmax, 0.d0)
      f = f + scale*(a1 * a1 + a2 * a2 + a3 * a3)
    else if(ityperest(irest).eq.4) then
      w = (xcart(icart,1)-restpars(irest,1))**2 + &
          (xcart(icart,2)-restpars(irest,2))**2 + &
          (xcart(icart,3)-restpars(irest,3))**2 - &
          restpars(irest,4)**2
      a1 = dmax1(w,0.d0)
      f = f + scale2*a1*a1
    else if(ityperest(irest).eq.5) then
      a1 = (xcart(icart,1)-restpars(irest,1))**2 / restpars(irest,4)**2
      a2 = (xcart(icart,2)-restpars(irest,2))**2 / restpars(irest,5)**2
      a3 = (xcart(icart,3)-restpars(irest,3))**2 / restpars(irest,6)**2
      a4 = restpars(irest,7)**2
      w = a1 + a2 + a3 - a4
      a1 = dmax1(w,0.d0)
      f = f + scale2*a1*a1
    else if(ityperest(irest).eq.6) then
      xmin = restpars(irest,1)
      ymin = restpars(irest,2)
      zmin = restpars(irest,3)
      xmax = restpars(irest,1) + restpars(irest,4)
      ymax = restpars(irest,2) + restpars(irest,4)
      zmax = restpars(irest,3) + restpars(irest,4)
      a1 = dmax1(xcart(icart,1) - xmin,0.d0)
      a2 = dmax1(xcart(icart,2) - ymin,0.d0)
      a3 = dmax1(xcart(icart,3) - zmin,0.d0)
      a4 = dmax1(xmax - xcart(icart,1),0.d0)
      a5 = dmax1(ymax - xcart(icart,2),0.d0)
      a6 = dmax1(zmax - xcart(icart,3),0.d0)
      f = f + a1*a2*a3*a4*a5*a6
    else if(ityperest(irest).eq.7) then
      xmin = restpars(irest,1)
      ymin = restpars(irest,2)
      zmin = restpars(irest,3)
      xmax = restpars(irest,4)
      ymax = restpars(irest,5)
      zmax = restpars(irest,6)
      a1 = dmax1(xcart(icart,1) - xmin,0.d0)
      a2 = dmax1(xcart(icart,2) - ymin,0.d0)
      a3 = dmax1(xcart(icart,3) - zmin,0.d0)
      a4 = dmax1(xmax - xcart(icart,1),0.d0)
      a5 = dmax1(ymax - xcart(icart,2),0.d0)
      a6 = dmax1(zmax - xcart(icart,3),0.d0)
      f = f + a1*a2*a3*a4*a5*a6
    else if(ityperest(irest).eq.8) then
      w = (xcart(icart,1)-restpars(irest,1))**2 + &
          (xcart(icart,2)-restpars(irest,2))**2 + &
          (xcart(icart,3)-restpars(irest,3))**2 - &
          restpars(irest,4)**2
      a1 = dmin1(w,0.d0)
      f = f + scale2*a1*a1
    else if(ityperest(irest).eq.9) then
      a1 = (xcart(icart,1)-restpars(irest,1))**2 / restpars(irest,4)**2
      a2 = (xcart(icart,2)-restpars(irest,2))**2 / restpars(irest,5)**2
      a3 = (xcart(icart,3)-restpars(irest,3))**2 / restpars(irest,6)**2
      a4 = restpars(irest,7)**2
      w = a1 + a2 + a3 - a4
      a1 = dmin1(w,0.d0)
      f = f + a1*a1
    else if(ityperest(irest).eq.10) then
      w = restpars(irest,1)*xcart(icart,1) + &
          restpars(irest,2)*xcart(icart,2) + &
          restpars(irest,3)*xcart(icart,3) - &
          restpars(irest,4)
      a1 = dmin1(w,0.d0)
      f = f + scale * a1*a1
    else if(ityperest(irest).eq.11) then
      w = restpars(irest,1)*xcart(icart,1) + &
          restpars(irest,2)*xcart(icart,2) + &
          restpars(irest,3)*xcart(icart,3) - &
          restpars(irest,4)
      a1 = dmax1(w,0.d0)
      f = f + scale * a1*a1
    else if(ityperest(irest).eq.12) then
      a1 = xcart(icart,1) - restpars(irest,1)
      a2 = xcart(icart,2) - restpars(irest,2)
      a3 = xcart(icart,3) - restpars(irest,3)
      vnorm = sqrt(restpars(irest,4)**2 + restpars(irest,5)**2 + restpars(irest,6)**2)
      v1 = restpars(irest,4)/vnorm
      v2 = restpars(irest,5)/vnorm
      v3 = restpars(irest,6)/vnorm
      b1 = v1 * a1
      b2 = v2 * a2
      b3 = v3 * a3
      w = b1 + b2 + b3
      d = ( a1 - v1*w )**2 + ( a2 - v2*w )**2 + ( a3 - v3*w )**2
      f = f + scale2 * ( &
           dmax1(-w , 0.d0)**2 + &
           dmax1(w - restpars(irest,9), 0.d0)**2 + &
           dmax1(d - restpars(irest,7)**2 , 0.d0 )**2 )
    else if(ityperest(irest).eq.13) then
      a1 = xcart(icart,1) - restpars(irest,1)
      a2 = xcart(icart,2) - restpars(irest,2)
      a3 = xcart(icart,3) - restpars(irest,3)
      vnorm = sqrt(restpars(irest,4)**2 + restpars(irest,5)**2 + restpars(irest,6)**2)
      v1 = restpars(irest,4)/vnorm
      v2 = restpars(irest,5)/vnorm
      v3 = restpars(irest,6)/vnorm
      b1 = v1 * a1
      b2 = v2 * a2
      b3 = v3 * a3
      w = b1 + b2 + b3
      d = ( a1 - v1*w )**2 +( a2 - v2*w )**2 + ( a3 - v3*w )**2
      f = f + scale2 * ( &
           dmin1(-w , 0.d0)**2 * &
           dmin1(w - restpars(irest,9), 0.d0)**2 * &
           dmin1(d - restpars(irest,7)**2 , 0.d0 )**2 )
    else if(ityperest(irest).eq.14) then
      w = restpars(irest,6)*exp( &
          -(xcart(icart,1) - restpars(irest,1))**2 &
          /(2*restpars(irest,3)**2) &
          -(xcart(icart,2) - restpars(irest,2))**2 &
          /(2*restpars(irest,4)**2)) &
          -(xcart(icart,3)-restpars(irest,5))
      a1 = dmax1(w,0.d0)
      f = f + scale * a1*a1
    else if(ityperest(irest).eq.15) then
      w = restpars(irest,6)*exp( &
          -(xcart(icart,1) - restpars(irest,1))**2 &
          /(2*restpars(irest,3)**2) &
          -(xcart(icart,2) - restpars(irest,2))**2 &
          /(2*restpars(irest,4)**2)) &
          -(xcart(icart,3)-restpars(irest,5))
      a1 = dmin1(w,0.d0)
      f = f + scale * a1*a1
    end if
  end do
  return
end subroutine comprest

max_stars_repo_name: SSchott/packmol

max_issues_repo_name: SSchott/packmol

max_forks_repo_name: SSchott/packmol


text:
!   -*- Mode: Fortran; -*-
!
!   (C) 2014 by Argonne National Laboratory.
!   See COPYRIGHT in top-level directory.
!
module mpi_f08_link_constants
use, intrinsic :: iso_c_binding, only: c_ptr, c_int, c_char, c_loc
use :: mpi_f08_types, only : MPI_Status
use :: mpi_c_interface_types, only : c_Status

implicit none

! Named link time constants

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! MPI_STATUS_IGNORE, MPI_STATUSES_IGNORE
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

! In MPICH C code, MPI_STATUS_IGNORE, MPI_STATUSES_IGNORE are pointers of bad address

! Fortran programmers use these MPI_Status variables for INOUT status arguments
type(MPI_Status), bind(C, name="MPIR_F08_MPI_STATUS_IGNORE_OBJ"), target :: MPI_STATUS_IGNORE
type(MPI_Status), dimension(1), bind(C, name="MPIR_F08_MPI_STATUSES_IGNORE_OBJ"), target :: MPI_STATUSES_IGNORE

type(c_ptr), bind(C, name="MPIR_C_MPI_STATUS_IGNORE") :: MPIR_C_MPI_STATUS_IGNORE
type(c_ptr), bind(C, name="MPIR_C_MPI_STATUSES_IGNORE") :: MPIR_C_MPI_STATUSES_IGNORE

! Though these two variables are required by MPI-3 Standard, they are not used in MPICH
type(c_ptr), bind(C, name="MPI_F08_STATUS_IGNORE")   :: MPI_F08_STATUS_IGNORE   ! Point to MPI_STATUS_IGNORE
type(c_ptr), bind(C, name="MPI_F08_STATUSES_IGNORE") :: MPI_F08_STATUSES_IGNORE ! Point to MPI_STATUSES_IGNORE

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!  MPI_ARGV_NULL, MPI_ARGVS_NULL
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! In MPICH C code, MPI_ARGV_NULL and MPI_ARGVS_NULL are NULL pointers
character(len=1), dimension(1), target :: MPI_ARGV_NULL
character(len=1), dimension(1,1), target :: MPI_ARGVS_NULL

type(c_ptr), bind(C, name="MPIR_C_MPI_ARGV_NULL") :: MPIR_C_MPI_ARGV_NULL
type(c_ptr), bind(C, name="MPIR_C_MPI_ARGVS_NULL") :: MPIR_C_MPI_ARGVS_NULL

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!  MPI_ERRCODES_IGNORE
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
integer, dimension(1), target :: MPI_ERRCODES_IGNORE
type(c_ptr), bind(C, name="MPIR_C_MPI_ERRCODES_IGNORE") :: MPIR_C_MPI_ERRCODES_IGNORE

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!  MPI_UNWEIGHTED, MPI_WEIGHTS_EMPTY
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Type of MPI_UNWEIGHTED and MPI_WEIGHTS_EMPTY is not integer(c_int) since we only care their address

integer, dimension(1), target :: MPI_UNWEIGHTED
integer, dimension(1), target :: MPI_WEIGHTS_EMPTY

type(c_ptr), protected, bind(C, name="MPIR_C_MPI_UNWEIGHTED")    :: MPIR_C_MPI_UNWEIGHTED
type(c_ptr), protected, bind(C, name="MPIR_C_MPI_WEIGHTS_EMPTY") :: MPIR_C_MPI_WEIGHTS_EMPTY

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!  MPI_IN_PLACE
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
integer(c_int), bind(C, name="MPIR_F08_MPI_IN_PLACE"), target :: MPI_IN_PLACE

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!  MPI_BOTTOM
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Buffer Address Constants
! A.1.1 p. 663
integer(c_int), bind(C, name="MPIR_F08_MPI_BOTTOM"), target :: MPI_BOTTOM

end module mpi_f08_link_constants

max_stars_repo_name: humairakamal/FG-MPI

max_issues_repo_name: ucd-plse/mpi-error-prop

max_forks_repo_name: ucd-plse/mpi-error-prop


text:
!======================================================================!
  subroutine UserProbe1D_kv2f(grid, name_out)
!----------------------------------------------------------------------!
! Reads the ".1D" file created by the "Generator" and averages the     !
! results in the planes defined by coordinates in it. Then averages    !
! the values of Umean, Vmean, Wmean, uu, vv, ww, uv, uw and vw and     !
! writes them into file ".1Dr".                                        !
!----------------------------------------------------------------------!
  use all_mod
  use allp_mod
  use les_mod
  use pro_mod
  use par_mod
  use rans_mod
!----------------------------------------------------------------------!
  implicit none
!-----------------------------[Arguments]------------------------------!
  real :: y(-NbC:NC)
!------------------------------[Calling]-------------------------------!
  interface
    logical function Approx(A,B,tol)
      real           :: A,B
      real, optional :: tol
    end function Approx
  end interface
!-------------------------------[Locals]-------------------------------!
  integer             :: Nprob, pl, c, dummy
  character           :: namCoo*80, namPro*80, answer*80
  character, optional :: name_out*(*)
  real,allocatable    :: z_p(:), Ump(:), Vmp(:), Wmp(:), &
                                 uup(:), vvp(:), wwp(:), &
                                 uvp(:), uwp(:), vwp(:), &
                                 Tmp(:), TTp(:),         &
                                 uTp(:), vTp(:), wTp(:), &
                                 Ksgsp(:), VISt_m(:), VISt_m1(:) ,VISt_m2(:),&
                                 Kin_mp(:), Eps_mp(:), v_2_mp(:), f22_mp(:), &
                                 var1(:), var2(:)
  integer,allocatable :: Np(:)
!======================================================================!

!  if(this  < 2)  &
!    write(*,*) '# Input probe file name [skip cancels]:'
!  call ReadC(7,inp,tn,ts,te)
!  read(inp(ts(1):te(1)),'(A80)') namPro

  if(PRESENT(name_out)) then
    write(*,*) name_out
!---- save the name
    namPro = name_out
  else
    namPro = 'linija'
    answer=namPro
    call ToUppr(answer)
    if(answer == 'SKIP') return
  end if

  call wait

!>>>>>>>>>>>>>>>>>>>>>>!
!     read 1D file     !
!>>>>>>>>>>>>>>>>>>>>>>!
  namCoo = name
  namCoo(len_trim(name)+1:len_trim(name)+3) = '.1D'
  write(6, *) '# Now reading the file:', namCoo
  open(9, file = namCoo)

!---- write the number of probes
  read(9,'(I8)') Nprob
  allocate(z_p(Nprob))

!---- write the probe coordinates out
  do pl=1,Nprob
    read(9,'(I8,1PE17.8)') dummy, z_p(pl)
  end do

  close(9)

  allocate(Np(Nprob));    Np=0
  allocate(Ump(Nprob));   Ump=0.0
  allocate(Vmp(Nprob));   Vmp=0.0
  allocate(Wmp(Nprob));   Wmp=0.0
  allocate(uup(Nprob));   uup=0.0
  allocate(vvp(Nprob));   vvp=0.0
  allocate(wwp(Nprob));   wwp=0.0
  allocate(uvp(Nprob));   uvp=0.0
  allocate(uwp(Nprob));   uwp=0.0
  allocate(vwp(Nprob));   vwp=0.0
  allocate(Ksgsp(Nprob)); Ksgsp=0.0
  allocate(VISt_m(Nprob)); VISt_m=0.0
  allocate(VISt_m1(Nprob)); VISt_m1=0.0
  allocate(VISt_m2(Nprob)); VISt_m2=0.0
  allocate(Kin_mp(Nprob)); Kin_mp=0.0
  allocate(Eps_mp(Nprob)); Eps_mp=0.0
  allocate(v_2_mp(Nprob)); v_2_mp=0.0
  allocate(f22_mp(Nprob)); f22_mp=0.0
  allocate(var1(Nprob)); var1=0.0
  allocate(var2(Nprob)); var2=0.0
  if(HOT==YES) then
    allocate(Tmp(Nprob));   Tmp=0.0
    allocate(TTp(Nprob));   TTp=0.0
    allocate(uTp(Nprob));   uTp=0.0
    allocate(vTp(Nprob));   vTp=0.0
    allocate(wTp(Nprob));   wTp=0.0
  end if

  call CalcShear(U % n, V % n, W % n, Shear)
  call GraPhi(U % n, 1, Ux,.TRUE.)    ! dU/dx
  call GraPhi(U % n, 2, Uy,.TRUE.)    ! dU/dy
  call GraPhi(U % mean, 3, Uz,.TRUE.)    ! dU/dz
  call GraPhi(V % n, 1, Vx,.TRUE.)    ! dV/dx
  call GraPhi(V % n, 2, Vy,.TRUE.)    ! dV/dy
  call GraPhi(V % n, 3, Vz,.TRUE.)    ! dV/dz
  call GraPhi(W % mean, 1, Wx,.TRUE.)    ! dW/dx
  call GraPhi(W % n, 2, Wy,.TRUE.)    ! dW/dy
  call GraPhi(W % mean, 3, Wz,.TRUE.)    ! dW/dz

!+++++++++++++++++++++++++++++!
!     average the results     !
!+++++++++++++++++++++++++++++!
  do pl=1,Nprob
    do c=1,NC
      if( Approx( zc(c),z_p(pl) ) ) then
        Np(pl)    = Np(pl) + 1
        Ump(pl)   = Ump(pl) + U % mean(c)
        Vmp(pl)   = Vmp(pl) + V % mean(c)
        Wmp(pl)   = Wmp(pl) + W % mean(c)
        uup(pl)   = uup(pl) + (uu % mean(c)-U % mean(c) * U % mean(c))
        vvp(pl)   = vvp(pl) + (vv % mean(c)-V % mean(c) * V % mean(c))
        wwp(pl)   = wwp(pl) + (ww % mean(c)-W % mean(c) * W % mean(c))
        uvp(pl)   = uvp(pl) + (uv % mean(c)-U % mean(c) * V % mean(c))
        uwp(pl)   = uwp(pl) + (uw % mean(c)-U % mean(c) * W % mean(c))
        vwp(pl)   = vwp(pl) + VISt(c)*abs(Uz(c))
!        vwp(pl)   = vwp(pl) + Cdyn(c)
!        vwp(pl)   = vwp(pl) + Kin % mean(c)
        VISt_m(pl)= VISt_m(pl) + Uf(c)
        VISt_m1(pl)= VISt_m1(pl) + VISt_r(c)
!        VISt_m2(pl)= VISt_m2(pl) + abs(uw % mean(c)-U % mean(c) * W % mean(c))/(abs(Uz(c))+tiny)
        VISt_m2(pl)= VISt_m2(pl) + VISt(c)
        Kin_mp(pl)= Kin_mp(pl) + Kin % n(c)
        Eps_mp(pl)= Eps_mp(pl) + Eps % n(c)
        v_2_mp(pl)= v_2_mp(pl) + v_2 % n(c)
        f22_mp(pl)= f22_mp(pl) + f22 % n(c)
        Ksgsp(pl) = Ksgsp(pl) + Ksgs(c)
        var2(pl)  = var2(pl) + Kin%mean(c)**1.5/(volume(c)**0.0333*0.01)
        if( Approx( xc(c),0.0490625 ) ) then
          if( Approx( yc(c),0.0490625 ) ) then
            var1(pl)  = var1(pl) + U % n(c)
          end if
        end if
        if(HOT==YES) then
          Tmp(pl)   = Tmp(pl) + T % mean(c)
          TTp(pl)   = TTp(pl) + (TT % mean(c) - T % mean(c) * T % mean(c))
          uTp(pl)   = uTp(pl) + (uT % mean(c) - u % mean(c) * T % mean(c))
          vTp(pl)   = vTp(pl) + (vT % mean(c) - v % mean(c) * T % mean(c))
          wTp(pl)   = wTp(pl) + (wT % mean(c) - w % mean(c) * T % mean(c))
        end if
      end if
    end do
  end do

!---- average over all processors
  do pl=1,Nprob
    call IGlSum(Np(pl))

    call GloSum(Ump(pl))
    call GloSum(Vmp(pl))
    call GloSum(Wmp(pl))

    call GloSum(var2(pl))

    call GloSum(uup(pl))
    call GloSum(vvp(pl))
    call GloSum(wwp(pl))

    call GloSum(uvp(pl))
    call GloSum(uwp(pl))
    call GloSum(vwp(pl))

    call GloSum(VISt_m(pl))
    call GloSum(VISt_m1(pl))
    call GloSum(VISt_m2(pl))

    call GloSum(Ksgsp(pl))

    call GloSum(Kin_mp(pl))
    call GloSum(Eps_mp(pl))
    call GloSum(v_2_mp(pl))
    call GloSum(f22_mp(pl))

    if(HOT==YES) then
      call GloSum(Tmp(pl))
      call GloSum(TTp(pl))
      call GloSum(uTp(pl))
      call GloSum(vTp(pl))
      call GloSum(wTp(pl))
    end if
  end do

  do pl=1,Nprob
    Ump(pl)   = Ump(pl) / Np(pl)
    Vmp(pl)   = Vmp(pl) / Np(pl)
    Wmp(pl)   = Wmp(pl) / Np(pl)
    var2(pl)  = var2(pl) / Np(pl)
    uup(pl)   = uup(pl) / Np(pl)
    vvp(pl)   = vvp(pl) / Np(pl)
    wwp(pl)   = wwp(pl) / Np(pl)
    uvp(pl)   = uvp(pl) / Np(pl)
    uwp(pl)   = uwp(pl) / Np(pl)
    vwp(pl)   = vwp(pl) / Np(pl)
    Ksgsp(pl) = Ksgsp(pl) / Np(pl)
    VISt_m(pl) = VISt_m(pl) / Np(pl)
    VISt_m1(pl) = VISt_m1(pl) / Np(pl)
    VISt_m2(pl) = VISt_m2(pl) / Np(pl)
    Kin_mp(pl) = Kin_mp(pl) / Np(pl)
    Eps_mp(pl) = Eps_mp(pl) / Np(pl)
    v_2_mp(pl) = v_2_mp(pl) / Np(pl)
    f22_mp(pl) = f22_mp(pl) / Np(pl)
    if(HOT==YES) then
      Tmp(pl) = Tmp(pl) / Np(pl)
      TTp(pl) = TTp(pl) / Np(pl)
      uTp(pl) = uTp(pl) / Np(pl)
      vTp(pl) = vTp(pl) / Np(pl)
      wTp(pl) = wTp(pl) / Np(pl)
    end if
  end do

!<<<<<<<<<<<<<<<<<<<<<<<<<!
!     create 1Dr file     !
!<<<<<<<<<<<<<<<<<<<<<<<<<!
  if(This < 2) then
    namPro(len_trim(namPro)+1:len_trim(namPro)+4) = '.1Dr'
    write(6, *) '# Now creating the file:', namPro
    open(9, file=namPro)

!---- write the number of probes
    write(9,'(A1,I8)') '#', Nprob
    write(9,'(A80)')  '# 1-z  2-Umean  3-Vmean  4-Wmean  5-uu  6-vv  7-ww  8-uv  9-uw  10-vw  11-Ksgs'
    if(HOT==YES) then
      write(9,'(A38)') '# 12-Tmean  13-TT  14-ut  15-vt  16-wt'
    end if

!---- write the values at the probes
    do pl=1,Nprob
      if(HOT==YES) then
        write(9,'(16E17.7,2I8)')                    &
                       z_p(pl),                     &
                       Ump(pl), Vmp(pl), Wmp(pl),   &
                       uup(pl), vvp(pl), wwp(pl),   &
                       uvp(pl), uwp(pl), vwp(pl),   &
                       Ksgsp(pl), Tmp(pl), TTp(pl), &
                       uTp(pl), vTp(pl), wTp(pl),   &
                       pl, Np(pl)
      else
        write(9,'(20E17.7)')                   &
                       z_p(pl),                    &
                       Ump(pl), Vmp(pl), Wmp(pl),  &
                       uup(pl), vvp(pl), wwp(pl),  &
                       uvp(pl), uwp(pl), vwp(pl),  &
                       VISt_m(pl), VISt_m1(pl),VISt_m2(pl), &
                       Kin_mp(pl), Eps_mp(pl), Ksgsp(pl), &
                       v_2_mp(pl), f22_mp(pl), var1(pl), var2(pl)
      end if
    end do

    close(9)
  end if

  deallocate(Np)
  deallocate(Ump)
  deallocate(Vmp)
  deallocate(Wmp)
  deallocate(var1)
  deallocate(var2)
  deallocate(uup)
  deallocate(vvp)
  deallocate(wwp)
  deallocate(uvp)
  deallocate(uwp)
  deallocate(vwp)
  deallocate(Ksgsp)
  deallocate(VISt_m)
  deallocate(VISt_m1)
  deallocate(VISt_m2)
  deallocate(Kin_mp)
  deallocate(Eps_mp)
  if(HOT==YES) then
    deallocate(Tmp)
    deallocate(TTp)
    deallocate(uTp)
    deallocate(vTp)
    deallocate(wTp)
  end if
  end subroutine UserProbe1D_kv2f

max_stars_repo_name: mhadziabdic/T-FlowS

max_issues_repo_name: mhadziabdic/T-FlowS

max_forks_repo_name: mhadziabdic/T-FlowS


MODULE generic_procedure_module
!
!  Purpose:
!    To define the derived data type for 2D vectors,
!    plus two generic bound procedures.
!
!  Record of revisions:
!       Date       Programmer          Description of change
!       ====       ==========          =====================
!     12/27/06    S. J. Chapman        Original code
!
IMPLICIT NONE

! Declare type vector
TYPE :: vector
   REAL :: x                         ! X value
   REAL :: y                         ! Y value
CONTAINS
   GENERIC :: add => vector_plus_vector, vector_plus_scalar
   PROCEDURE,PASS :: vector_plus_vector
   PROCEDURE,PASS :: vector_plus_scalar
END TYPE vector

! Add procedures
CONTAINS

   TYPE (vector) FUNCTION vector_plus_vector ( this, v2 )
   !
   !  Purpose:
   !    To add two vectors.
   !
   !  Record of revisions:
   !       Date       Programmer          Description of change
   !       ====       ==========          =====================
   !     12/27/06    S. J. Chapman        Original code
   !
   IMPLICIT NONE

   ! Data dictionary: declare calling parameter types & definitions
   CLASS(vector),INTENT(IN) :: this         ! First vector
   CLASS(vector),INTENT(IN) :: v2           ! Second vector

   ! Add the vectors
   vector_plus_vector%x = this%x + v2%x
   vector_plus_vector%y = this%y + v2%y

   END FUNCTION vector_plus_vector

   TYPE (vector) FUNCTION vector_plus_scalar ( this, s )
   !
   !  Purpose:
   !    To add a vector and a scalar.
   !
   !  Record of revisions:
   !       Date       Programmer          Description of change
   !       ====       ==========          =====================
   !     12/27/06    S. J. Chapman        Original code
   !
   IMPLICIT NONE

   ! Data dictionary: declare calling parameter types & definitions
   CLASS(vector),INTENT(IN) :: this         ! First vector
   REAL,INTENT(IN) :: s                     ! Scalar

   ! Add the points
   vector_plus_scalar%x = this%x + s
   vector_plus_scalar%y = this%y + s

   END FUNCTION vector_plus_scalar

END MODULE generic_procedure_module


text: #redirect Sheila the Dog

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
A Departed businesses departed business, EP Copier was a tiny Copy Shops copy shop that had four copy machines. For under 50 copies they charged $.06 +tax per sheet. Over 50 copies the price dropped to a cheap $.05 +tax. They also sold a small selection of Chinese media newspapers ($.50.75), office supplies and womens purses. They may have also been located at 140 B Street at some point.

The closed their doors as of September 1, 2005.

A business called HELP Copier now exists in the same place, providing the same services. Did they change owners, or just names? Users/JessicaLuedtke
  Maybe someone else knows better since I havent seen the new copy shop. But the old EP Copier had completely vacated their store when they closed. Users/JackHaskel


text:
	SUBROUTINE UT_MDQI  ( idxstn, idxvar,
     +			      obs, qcd, iqca, iqcr, mxstns,
     +			      rval, cvalqd, rvalqa, rvalqr, iret )
C************************************************************************
C* UT_MDQI								*
C*									*
C* This subroutine retrieves an interface value (and all associated QC	*
C* values!) for a specified variable at a specified station from the	*
C* MADIS software observation and QC arrays.				*
C*									*
C* UT_MDQI  ( IDXSTN, IDXVAR,						*
C*	      OBS, QCD, IQCA, IQCR, MXSTNS,				*
C*	      RVAL, CVALQD, RVALQA, RVALQR, IRET )			*
C*									*
C* Input parameters:							*
C*	IDXSTN		INTEGER		Station index into MADIS arrays	*
C*	IDXVAR		INTEGER		Variable index into MADIS arrays*
C*	OBS(*,*)	REAL		MADIS observations array	*
C*	QCD(*,*)	CHARACTER	QCD values associated with OBS	*
C*	IQCA(*,*)	INTEGER		QCA values associated with OBS	*
C*	IQCR(*,*)	INTEGER		QCR values associated with OBS	*
C*	MXSTNS		INTEGER		Maximum number of stations	*
C*					within MADIS arrays		*
C*									*
C* Output parameters:							*
C*	RVAL		REAL		Interface value			*
C*	CVALQD		CHARACTER	QCD value associated with RVAL	*
C*	RVALQA		REAL		QCA value associated with RVAL	*
C*	RVALQR		REAL		QCA value associated with RVAL	*
C*	IRET		INTEGER		Return code 			*
C*					  0 = normal return 		*
C*									*
C**									*
C* Log:									*
C* J. Ator/NCEP		06/01						*
C* J. Ator/NCEP		09/01	Modify to use FSL MADIS software	*
C************************************************************************
	REAL		obs ( mxstns, * )
C*
	CHARACTER	qcd ( mxstns, * ), cvalqd
C*
	INTEGER		iqca ( mxstns, * ), iqcr ( mxstns, * )
C-----------------------------------------------------------------------
	iret = 0
C
	rval = UT_MDRI ( obs ( idxstn, idxvar ) )
C
	cvalqd = qcd ( idxstn, idxvar )
	rvalqa = UT_MDRI ( FLOAT ( iqca ( idxstn, idxvar ) ) )
	rvalqr = UT_MDRI ( FLOAT ( iqcr ( idxstn, idxvar ) ) )
C*
	RETURN
	END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


!
! Distributed under the OSI-approved Apache License, Version 2.0.  See
!  accompanying file Copyright.txt for details.
!
!  adios2_io_mod.f90 : ADIOS2 Fortran bindings for IO class
!
!   Created on: Mar 13, 2017
!       Author: William F Godoy godoywf@ornl.gov
!

module adios2_io_mod
    use adios2_io_open_mod
    use adios2_io_define_variable_mod
    use adios2_io_define_attribute_mod
    implicit none

contains

    subroutine adios2_set_engine(io, engine_type, ierr)
        type(adios2_io), intent(in) :: io
        character*(*), intent(in) :: engine_type
        integer, intent(out) :: ierr

        call adios2_set_engine_f2c(io%f2c, &
                                   TRIM(ADJUSTL(engine_type))//char(0), &
                                   ierr)
    end subroutine

    subroutine adios2_set_parameter(io, key, value, ierr)
        type(adios2_io), intent(in) :: io
        character*(*), intent(in) :: key
        character*(*), intent(in) :: value
        integer, intent(out) :: ierr

        call adios2_set_parameter_f2c(io%f2c, TRIM(ADJUSTL(key))//char(0), &
                                      TRIM(ADJUSTL(value))//char(0), ierr)
    end subroutine

    subroutine adios2_add_transport(transport_index, io, type, ierr)
        integer, intent(out):: transport_index
        type(adios2_io), intent(in) :: io
        character*(*), intent(in) :: type
        integer, intent(out) :: ierr

        call adios2_add_transport_f2c(transport_index, io%f2c, &
                                      TRIM(ADJUSTL(type))//char(0), ierr)

    end subroutine

    subroutine adios2_set_transport_parameter(io, transport_index, key, value, &
                                              ierr)
        type(adios2_io), intent(in):: io
        integer, intent(in):: transport_index
        character*(*), intent(in) :: key
        character*(*), intent(in) :: value
        integer, intent(out):: ierr

        call adios2_set_transport_parameter_f2c(io%f2c, transport_index, &
                                                TRIM(ADJUSTL(key))//char(0), &
                                                TRIM(ADJUSTL(value))//char(0), &
                                                ierr)
    end subroutine

    subroutine adios2_inquire_variable(variable, io, name, ierr)
        type(adios2_variable), intent(out) :: variable
        type(adios2_io), intent(in) :: io
        character*(*), intent(in) :: name
        integer, intent(out) :: ierr

        call adios2_inquire_variable_f2c(variable%f2c, io%f2c, &
                                         TRIM(ADJUSTL(name))//char(0), ierr)
    end subroutine

    subroutine adios2_remove_variable(io, name, ierr)
        type(adios2_io), intent(in) :: io
        character*(*), intent(in) :: name
        integer, intent(out) :: ierr
        ! Local
        type(adios2_variable):: variable

        call adios2_inquire_variable(variable, io, name, ierr)
        if(ierr == adios2_found) then
            variable%valid = .false.
        end if

        call adios2_remove_variable_f2c(io%f2c, TRIM(ADJUSTL(name))//char(0), ierr)

    end subroutine

    subroutine adios2_remove_all_variables(io, ierr)
        type(adios2_io), intent(in) :: io
        integer, intent(out) :: ierr

        call adios2_remove_all_variables_f2c(io%f2c, ierr)

    end subroutine

    subroutine adios2_inquire_attribute(attribute, io, name, ierr)
        type(adios2_attribute), intent(out) :: attribute
        type(adios2_io), intent(in) :: io
        character*(*), intent(in) :: name
        integer, intent(out) :: ierr

        call adios2_inquire_attribute_f2c(attribute%f2c, io%f2c, &
                                          TRIM(ADJUSTL(name))//char(0), ierr)
    end subroutine

    subroutine adios2_remove_attribute(io, name, ierr)
        type(adios2_io), intent(in) :: io
        character*(*), intent(in) :: name
        integer, intent(out) :: ierr

        ! Local
        type(adios2_attribute):: attribute

        call adios2_inquire_attribute(attribute, io, name, ierr)
        if(ierr == adios2_found) then
            attribute%valid = .false.
        end if

        call adios2_remove_attribute_f2c(io%f2c, TRIM(ADJUSTL(name))//char(0), &
                                         ierr)
    end subroutine

    subroutine adios2_remove_all_attributes(io, ierr)
        type(adios2_io), intent(in) :: io
        integer, intent(out) :: ierr

        call adios2_remove_all_attributes_f2c(io%f2c, ierr)

    end subroutine

    subroutine adios2_flush_all_engines(io, ierr)
        type(adios2_io), intent(in) :: io
        integer, intent(out) :: ierr

        call adios2_flush_all_engines_f2c(io%f2c, ierr)

    end subroutine

    subroutine adios2_io_engine_type(io, engine_type, ierr)
        type(adios2_io), intent(in) :: io
        character(len=:), allocatable, intent(out) :: engine_type
        integer, intent(out) :: ierr

        character(len=32) :: c_engine_type
        integer :: length

        call adios2_io_engine_type_f2c(io%f2c, c_engine_type, length, ierr)
        call adios2_StringC2F(c_engine_type, length, engine_type)

    end subroutine

end module


text:
      subroutine tflifetrack(isp1,kx,irtc)
      use tfstk
      implicit none
      integer*8 kx
      integer*4 isp1,irtc,itfmessage,nip,i,lstr
      parameter (lstr=512)
      real*8 gw
      character*(lstr) file,comm
      if(isp .ne. isp1+13)then
        irtc=itfmessage(9,'General::narg','"13"')
        return
      endif
      if(ktfnonrealq(ktastk(isp1+10)))then
        go to 9000
      endif
      nip=int(rtastk(isp1+10))
      if(nip .le. 0)then
        irtc=itfmessage(9,'General::wrongval',
     $       '" # of IP must be positive"')
        return
      endif
      do i=1,9
        if(.not. tfreallistq(ktastk(isp1+i)))then
          go to 9000
        endif
      enddo
      if(ilist(2,ktfaddr(ktastk(isp1+1))-1) .ne. 3)then
        irtc=itfmessage(9,'General::wrongval',
     $       '"EmittanceStrong={ex,ey,ez}"')
        return
      endif
      if(ilist(2,ktfaddr(ktastk(isp1+2))-1) .ne. 36*nip)then
        irtc=itfmessage(9,'General::wrongval',
     $       '"Dim[BeamMatrixStrong]={#IP,6,6}"')
        return
      endif
      if(ilist(2,ktfaddr(ktastk(isp1+3))-1) .ne. 3)then
        irtc=itfmessage(9,'General::wrongval',
     $       '"EmittanceWeak={ex,ey,ez}"')
        return
      endif
      if(ilist(2,ktfaddr(ktastk(isp1+4))-1) .ne. 36*nip)then
        irtc=itfmessage(9,'General::wrongval',
     $       '"Dim[BeamMatrixWeak]={#IP,6,6}"')
        return
      endif
      if(ilist(2,ktfaddr(ktastk(isp1+5))-1) .ne. 6*nip)then
        irtc=itfmessage(9,'General::wrongval',
     $       '"Dim[OrbitWeak]={#IP,6}"')
        return
      endif
      if(ilist(2,ktfaddr(ktastk(isp1+6))-1) .ne. 6*nip)then
        irtc=itfmessage(9,'General::wrongval',
     $       '"Dim[Separations]={#IP,6}"')
        return
      endif
      if(ilist(2,ktfaddr(ktastk(isp1+7))-1) .ne. 3)then
        irtc=itfmessage(9,'General::wrongval',
     $       '"DampingRate={dx,dy,dz}"')
        return
      endif
      if(ilist(2,ktfaddr(ktastk(isp1+8))-1) .ne. 3)then
        irtc=itfmessage(9,'General::wrongval',
     $       '"Apertures={ax,ay,az}"')
        return
      endif
      if(ilist(2,ktfaddr(ktastk(isp1+9))-1) .ne. nip)then
        irtc=itfmessage(9,'General::wrongval',
     $       '"Dim[Particles]={#IP}"')
        return
      endif
      if(ktfnonrealq(ktastk(isp1+11)))then
        irtc=itfmessage(9,'General::wrongval',
     $       '"GammaWeak=Real"')
        return
      endif
      if(ktfnonstringq(ktastk(isp1+12)))then
        irtc=itfmessage(9,'General::wrongval',
     $       '"OutputFile=String"')
        return
      endif
      if(ktfnonstringq(ktastk(isp1+13)))then
        irtc=itfmessage(9,'General::wrongval',
     $       '"Comment=String"')
        return
      endif
      gw=rtastk(isp1+11)
      file=""
      call tmovb(rlist(ktfaddr(ktastk(isp1+12))+1),file,
     $     min(ilist(1,ktfaddr(ktastk(isp1+12))),lstr))
      comm=""
      call tmovb(rlist(ktfaddr(ktastk(isp1+12))+1),comm,
     $     min(ilist(1,ktfaddr(ktastk(isp1+12))),lstr))
      call sad_lifetrack()
c$$$     $     rlist(ktfaddr(ktastk(isp1+1))),
c$$$     $     rlist(ktfaddr(ktastk(isp1+2))),
c$$$     $     rlist(ktfaddr(ktastk(isp1+3))),
c$$$     $     rlist(ktfaddr(ktastk(isp1+4))),
c$$$     $     rlist(ktfaddr(ktastk(isp1+5))),
c$$$     $     rlist(ktfaddr(ktastk(isp1+6))),
c$$$     $     rlist(ktfaddr(ktastk(isp1+7))),
c$$$     $     rlist(ktfaddr(ktastk(isp1+8))),
c$$$     $     rlist(ktfaddr(ktastk(isp1+9))),
c$$$     $     nip,gw,file,comm)
      irtc=0
      kx=ktfoper+mtfnull
      return
 9000 irtc=itfmessage(9,'General::wrongtype',
     $       '"9*NumList,nip,GammaWeak,Out_String,Comment_String"')
      return
      end

      subroutine sad_lifetrack
      return
      end

max_stars_repo_name: noboruatkek/SAD

max_issues_repo_name: noboruatkek/SAD

max_forks_repo_name: noboruatkek/SAD


!
!                Parallel Sparse BLAS  version 3.5
!      (C) Copyright 2006-2018
!        Salvatore Filippone
!        Alfredo Buttari
!
!    Redistribution and use in source and binary forms, with or without
!    modification, are permitted provided that the following conditions
!    are met:
!      1. Redistributions of source code must retain the above copyright
!         notice, this list of conditions and the following disclaimer.
!      2. Redistributions in binary form must reproduce the above copyright
!         notice, this list of conditions, and the following disclaimer in the
!         documentation and/or other materials provided with the distribution.
!      3. The name of the PSBLAS group or the names of its contributors may
!         not be used to endorse or promote products derived from this
!         software without specific written permission.
!
!    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
!    ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
!    TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
!    PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS
!    BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
!    CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
!    SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
!    INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
!    CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
!    ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
!    POSSIBILITY OF SUCH DAMAGE.
!
!
!
!  Sorting routines
!  References:
!  D. Knuth
!  The Art of Computer Programming, vol. 3
!  Addison-Wesley
!
!  Aho, Hopcroft, Ullman
!  Data Structures and Algorithms
!  Addison-Wesley
!
module psb_d_sort_mod
  use psb_const_mod

  @INTE@

  interface psb_msort_unique
    subroutine psb_dmsort_u(x,nout,dir)
      import
      real(psb_dpk_), intent(inout)           :: x(:)
      integer(psb_ipk_), intent(out)             :: nout
      integer(psb_ipk_), optional, intent(in)    :: dir
    end subroutine psb_dmsort_u
  end interface psb_msort_unique

  type psb_d_heap
    integer(psb_ipk_) :: last, dir
    real(psb_dpk_), allocatable    :: keys(:)
  contains
    procedure, pass(heap) :: init       => psb_d_init_heap
    procedure, pass(heap) :: howmany    => psb_d_howmany
    procedure, pass(heap) :: insert     => psb_d_insert_heap
    procedure, pass(heap) :: get_first  => psb_d_heap_get_first
    procedure, pass(heap) :: dump       => psb_d_dump_heap
    procedure, pass(heap) :: free       => psb_d_free_heap
  end type psb_d_heap

  type psb_d_idx_heap
    integer(psb_ipk_) :: last, dir
    real(psb_dpk_), allocatable    :: keys(:)
    integer(psb_ipk_), allocatable :: idxs(:)
  contains
    procedure, pass(heap) :: init       => psb_d_idx_init_heap
    procedure, pass(heap) :: howmany    => psb_d_idx_howmany
    procedure, pass(heap) :: insert     => psb_d_idx_insert_heap
    procedure, pass(heap) :: get_first  => psb_d_idx_heap_get_first
    procedure, pass(heap) :: dump       => psb_d_idx_dump_heap
    procedure, pass(heap) :: free       => psb_d_idx_free_heap
  end type psb_d_idx_heap

  interface psb_msort
    subroutine psb_dmsort(x,ix,dir,flag)
      import
      real(psb_dpk_), intent(inout)  :: x(:)
      integer(psb_ipk_), optional, intent(in)    :: dir, flag
      integer(psb_ipk_), optional, intent(inout) :: ix(:)
    end subroutine psb_dmsort
  end interface psb_msort

  interface psb_bsrch
    function  psb_dbsrch(key,n,v) result(ipos)
      import
      integer(psb_ipk_) :: ipos, n
      real(psb_dpk_) :: key
      real(psb_dpk_) :: v(:)
    end function psb_dbsrch
  end interface psb_bsrch

  interface psb_ssrch
    function psb_dssrch(key,n,v) result(ipos)
      import
      implicit none
      integer(psb_ipk_) :: ipos, n
      real(psb_dpk_) :: key
      real(psb_dpk_) :: v(:)
    end function psb_dssrch
  end interface psb_ssrch

  interface
    subroutine psi_d_msort_up(n,k,l,iret)
      import
      implicit none
      integer(psb_ipk_) :: n, iret
      real(psb_dpk_)  ::  k(n)
      integer(psb_ipk_) :: l(0:n+1)
    end subroutine psi_d_msort_up
    subroutine psi_d_msort_dw(n,k,l,iret)
      import
      implicit none
      integer(psb_ipk_) :: n, iret
      real(psb_dpk_)  ::  k(n)
      integer(psb_ipk_) :: l(0:n+1)
    end subroutine psi_d_msort_dw
  end interface
  interface
    subroutine psi_d_amsort_up(n,k,l,iret)
      import
      implicit none
      integer(psb_ipk_) :: n, iret
      real(psb_dpk_)  ::  k(n)
      integer(psb_ipk_) :: l(0:n+1)
    end subroutine psi_d_amsort_up
    subroutine psi_d_amsort_dw(n,k,l,iret)
      import
      implicit none
      integer(psb_ipk_) :: n, iret
      real(psb_dpk_)  ::  k(n)
      integer(psb_ipk_) :: l(0:n+1)
    end subroutine psi_d_amsort_dw
  end interface

  interface psb_qsort
    subroutine psb_dqsort(x,ix,dir,flag)
      import
      real(psb_dpk_), intent(inout)  :: x(:)
      integer(psb_ipk_), optional, intent(in)    :: dir, flag
      integer(psb_ipk_), optional, intent(inout) :: ix(:)
    end subroutine psb_dqsort
  end interface psb_qsort

  interface psb_isort
    subroutine psb_disort(x,ix,dir,flag)
      import
      real(psb_dpk_), intent(inout)  :: x(:)
      integer(psb_ipk_), optional, intent(in)    :: dir, flag
      integer(psb_ipk_), optional, intent(inout) :: ix(:)
    end subroutine psb_disort
  end interface psb_isort

  interface psb_hsort
    subroutine psb_dhsort(x,ix,dir,flag)
      import
      real(psb_dpk_), intent(inout)  :: x(:)
      integer(psb_ipk_), optional, intent(in)    :: dir, flag
      integer(psb_ipk_), optional, intent(inout) :: ix(:)
    end subroutine psb_dhsort
  end interface psb_hsort

  interface
    subroutine psi_d_insert_heap(key,last,heap,dir,info)
      import
      implicit none

      !
      ! Input:
      !   key:  the new value
      !   last: pointer to the last occupied element in heap
      !   heap: the heap
      !   dir:  sorting direction

      real(psb_dpk_), intent(in)     :: key
      real(psb_dpk_), intent(inout)  :: heap(:)
      integer(psb_ipk_), intent(in)     :: dir
      integer(psb_ipk_), intent(inout)  :: last
      integer(psb_ipk_), intent(out)    :: info
    end subroutine psi_d_insert_heap
  end interface

  interface
    subroutine psi_d_idx_insert_heap(key,index,last,heap,idxs,dir,info)
      import
      implicit none

      !
      ! Input:
      !   key:  the new value
      !   last: pointer to the last occupied element in heap
      !   heap: the heap
      !   dir:  sorting direction

      real(psb_dpk_), intent(in)     :: key
      real(psb_dpk_), intent(inout)  :: heap(:)
      integer(psb_ipk_), intent(in)     :: index
      integer(psb_ipk_), intent(in)     :: dir
      integer(psb_ipk_), intent(inout)  :: idxs(:)
      integer(psb_ipk_), intent(inout)  :: last
      integer(psb_ipk_), intent(out)    :: info
    end subroutine psi_d_idx_insert_heap
  end interface

  interface
    subroutine psi_d_heap_get_first(key,last,heap,dir,info)
      import
      implicit none
      real(psb_dpk_), intent(inout)  :: key
      integer(psb_ipk_), intent(inout)  :: last
      integer(psb_ipk_), intent(in)     :: dir
      real(psb_dpk_), intent(inout)  :: heap(:)
      integer(psb_ipk_), intent(out)    :: info
    end subroutine psi_d_heap_get_first
  end interface

  interface
    subroutine psi_d_idx_heap_get_first(key,index,last,heap,idxs,dir,info)
      import
      real(psb_dpk_), intent(inout)    :: key
      integer(psb_ipk_), intent(out)    :: index
      real(psb_dpk_), intent(inout)    :: heap(:)
      integer(psb_ipk_), intent(in)     :: dir
      integer(psb_ipk_), intent(inout)  :: last
      integer(psb_ipk_), intent(inout)  :: idxs(:)
      integer(psb_ipk_), intent(out)    :: info
    end subroutine psi_d_idx_heap_get_first
  end interface

  interface
    subroutine psi_disrx_up(n,x,ix)
      import
      real(psb_dpk_), intent(inout)  :: x(:)
      integer(psb_ipk_), intent(inout) :: ix(:)
      integer(psb_ipk_), intent(in)   :: n
    end subroutine psi_disrx_up
    subroutine psi_disrx_dw(n,x,ix)
      import
      real(psb_dpk_), intent(inout)  :: x(:)
      integer(psb_ipk_), intent(inout) :: ix(:)
      integer(psb_ipk_), intent(in)   :: n
    end subroutine psi_disrx_dw
    subroutine psi_disr_up(n,x)
      import
      real(psb_dpk_), intent(inout)  :: x(:)
      integer(psb_ipk_), intent(in)   :: n
    end subroutine psi_disr_up
    subroutine psi_disr_dw(n,x)
      import
      real(psb_dpk_), intent(inout)  :: x(:)
      integer(psb_ipk_), intent(in)   :: n
    end subroutine psi_disr_dw
    subroutine psi_daisrx_up(n,x,ix)
      import
      real(psb_dpk_), intent(inout)  :: x(:)
      integer(psb_ipk_), intent(inout) :: ix(:)
      integer(psb_ipk_), intent(in)   :: n
    end subroutine psi_daisrx_up
    subroutine psi_daisrx_dw(n,x,ix)
      import
      real(psb_dpk_), intent(inout)  :: x(:)
      integer(psb_ipk_), intent(inout) :: ix(:)
      integer(psb_ipk_), intent(in)   :: n
    end subroutine psi_daisrx_dw
    subroutine psi_daisr_up(n,x)
      import
      real(psb_dpk_), intent(inout)  :: x(:)
      integer(psb_ipk_), intent(in)   :: n
    end subroutine psi_daisr_up
    subroutine psi_daisr_dw(n,x)
      import
      real(psb_dpk_), intent(inout)  :: x(:)
      integer(psb_ipk_), intent(in)   :: n
    end subroutine psi_daisr_dw
  end interface

  interface
    subroutine psi_dqsrx_up(n,x,ix)
      import
      real(psb_dpk_), intent(inout)  :: x(:)
      integer(psb_ipk_), intent(inout) :: ix(:)
      integer(psb_ipk_), intent(in)   :: n
    end subroutine psi_dqsrx_up
    subroutine psi_dqsrx_dw(n,x,ix)
      import
      real(psb_dpk_), intent(inout)  :: x(:)
      integer(psb_ipk_), intent(inout) :: ix(:)
      integer(psb_ipk_), intent(in)   :: n
    end subroutine psi_dqsrx_dw
    subroutine psi_dqsr_up(n,x)
      import
      real(psb_dpk_), intent(inout)  :: x(:)
      integer(psb_ipk_), intent(in)   :: n
    end subroutine psi_dqsr_up
    subroutine psi_dqsr_dw(n,x)
      import
      real(psb_dpk_), intent(inout)  :: x(:)
      integer(psb_ipk_), intent(in)   :: n
    end subroutine psi_dqsr_dw
    subroutine psi_daqsrx_up(n,x,ix)
      import
      real(psb_dpk_), intent(inout)  :: x(:)
      integer(psb_ipk_), intent(inout) :: ix(:)
      integer(psb_ipk_), intent(in)   :: n
    end subroutine psi_daqsrx_up
    subroutine psi_daqsrx_dw(n,x,ix)
      import
      real(psb_dpk_), intent(inout)  :: x(:)
      integer(psb_ipk_), intent(inout) :: ix(:)
      integer(psb_ipk_), intent(in)   :: n
    end subroutine psi_daqsrx_dw
    subroutine psi_daqsr_up(n,x)
      import
      real(psb_dpk_), intent(inout)  :: x(:)
      integer(psb_ipk_), intent(in)   :: n
    end subroutine psi_daqsr_up
    subroutine psi_daqsr_dw(n,x)
      import
      real(psb_dpk_), intent(inout)  :: x(:)
      integer(psb_ipk_), intent(in)   :: n
    end subroutine psi_daqsr_dw
  end interface

contains

  subroutine psb_d_init_heap(heap,info,dir)
    use psb_realloc_mod, only : psb_ensure_size
    implicit none
    class(psb_d_heap), intent(inout) :: heap
    integer(psb_ipk_), intent(out)            :: info
    integer(psb_ipk_), intent(in), optional   :: dir

    info = psb_success_
    heap%last=0
    if (present(dir)) then
      heap%dir = dir
    else
      heap%dir = psb_sort_up_
    endif
    select case(heap%dir)
    case (psb_sort_up_,psb_sort_down_,psb_asort_up_,psb_asort_down_)
      ! ok, do nothing
    case default
      write(psb_err_unit,*) 'Invalid direction, defaulting to psb_sort_up_'
      heap%dir = psb_sort_up_
    end select
    call psb_ensure_size(psb_heap_resize,heap%keys,info)

    return
  end subroutine psb_d_init_heap

  function psb_d_howmany(heap) result(res)
    implicit none
    class(psb_d_heap), intent(in) :: heap
    integer(psb_ipk_) :: res
    res  = heap%last
  end function psb_d_howmany

  subroutine psb_d_insert_heap(key,heap,info)
    use psb_realloc_mod, only : psb_ensure_size
    implicit none

    real(psb_dpk_), intent(in)              :: key
    class(psb_d_heap), intent(inout) :: heap
    integer(psb_ipk_), intent(out)                       :: info

    info = psb_success_
    if (heap%last < 0) then
      write(psb_err_unit,*) 'Invalid last in heap ',heap%last
      info = heap%last
      return
    endif

    call psb_ensure_size(heap%last+1,heap%keys,info,addsz=psb_heap_resize)
    if (info /= psb_success_) then
      write(psb_err_unit,*) 'Memory allocation failure in heap_insert'
      info = -5
      return
    end if
    call psi_d_insert_heap(key,&
         & heap%last,heap%keys,heap%dir,info)

    return
  end subroutine psb_d_insert_heap

  subroutine psb_d_heap_get_first(key,heap,info)
    implicit none

    class(psb_d_heap), intent(inout) :: heap
    integer(psb_ipk_), intent(out)     :: info
    real(psb_dpk_), intent(out)       :: key

    info = psb_success_

    call psi_d_heap_get_first(key,&
         & heap%last,heap%keys,heap%dir,info)

    return
  end subroutine psb_d_heap_get_first

  subroutine psb_d_dump_heap(iout,heap,info)

    implicit none
    class(psb_d_heap), intent(in) :: heap
    integer(psb_ipk_), intent(out)    :: info
    integer(psb_ipk_), intent(in)     :: iout

    info = psb_success_
    if (iout < 0) then
      write(psb_err_unit,*) 'Invalid file '
      info =-1
      return
    end if

    write(iout,*) 'Heap direction ',heap%dir
    write(iout,*) 'Heap size      ',heap%last
    if ((heap%last > 0).and.((.not.allocated(heap%keys)).or.&
         & (size(heap%keys)<heap%last))) then
      write(iout,*) 'Inconsistent size/allocation status!!'
    else
      write(iout,*) heap%keys(1:heap%last)
    end if
  end subroutine psb_d_dump_heap

  subroutine psb_d_free_heap(heap,info)
    implicit none
    class(psb_d_heap), intent(inout) :: heap
    integer(psb_ipk_), intent(out)           :: info

    info=psb_success_
    if (allocated(heap%keys)) deallocate(heap%keys,stat=info)

  end subroutine psb_d_free_heap

  subroutine psb_d_idx_init_heap(heap,info,dir)
    use psb_realloc_mod, only : psb_ensure_size
    implicit none
    class(psb_d_idx_heap), intent(inout) :: heap
    integer(psb_ipk_), intent(out)            :: info
    integer(psb_ipk_), intent(in), optional   :: dir

    info = psb_success_
    heap%last=0
    if (present(dir)) then
      heap%dir = dir
    else
      heap%dir = psb_sort_up_
    endif
    select case(heap%dir)
    case (psb_sort_up_,psb_sort_down_,psb_asort_up_,psb_asort_down_)
      ! ok, do nothing
    case default
      write(psb_err_unit,*) 'Invalid direction, defaulting to psb_sort_up_'
      heap%dir = psb_sort_up_
    end select

    call psb_ensure_size(psb_heap_resize,heap%keys,info)
    call psb_ensure_size(psb_heap_resize,heap%idxs,info)
    return
  end subroutine psb_d_idx_init_heap

  function psb_d_idx_howmany(heap) result(res)
    implicit none
    class(psb_d_idx_heap), intent(in) :: heap
    integer(psb_ipk_) :: res
    res  = heap%last
  end function psb_d_idx_howmany

  subroutine psb_d_idx_insert_heap(key,index,heap,info)
    use psb_realloc_mod, only : psb_ensure_size
    implicit none

    real(psb_dpk_), intent(in)              :: key
    integer(psb_ipk_), intent(in)                        :: index
    class(psb_d_idx_heap), intent(inout) :: heap
    integer(psb_ipk_), intent(out)                       :: info

    info = psb_success_
    if (heap%last < 0) then
      write(psb_err_unit,*) 'Invalid last in heap ',heap%last
      info = heap%last
      return
    endif

    call psb_ensure_size(heap%last+1,heap%keys,info,addsz=psb_heap_resize)
    if (info == psb_success_) &
         & call psb_ensure_size(heap%last+1,heap%idxs,info,addsz=psb_heap_resize)
    if (info /= psb_success_) then
      write(psb_err_unit,*) 'Memory allocation failure in heap_insert'
      info = -5
      return
    end if
    call psi_d_idx_insert_heap(key,index,&
         & heap%last,heap%keys,heap%idxs,heap%dir,info)

    return
  end subroutine psb_d_idx_insert_heap

  subroutine psb_d_idx_heap_get_first(key,index,heap,info)
    implicit none

    class(psb_d_idx_heap), intent(inout) :: heap
    integer(psb_ipk_), intent(out)       :: index
    integer(psb_ipk_), intent(out)       :: info
    real(psb_dpk_), intent(out)           :: key

    info = psb_success_

    call psi_d_idx_heap_get_first(key,index,&
         & heap%last,heap%keys,heap%idxs,heap%dir,info)

    return
  end subroutine psb_d_idx_heap_get_first

  subroutine psb_d_idx_dump_heap(iout,heap,info)

    implicit none
    class(psb_d_idx_heap), intent(in) :: heap
    integer(psb_ipk_), intent(out)    :: info
    integer(psb_ipk_), intent(in)     :: iout

    info = psb_success_
    if (iout < 0) then
      write(psb_err_unit,*) 'Invalid file '
      info =-1
      return
    end if

    write(iout,*) 'Heap direction ',heap%dir
    write(iout,*) 'Heap size      ',heap%last
    if ((heap%last > 0).and.((.not.allocated(heap%keys)).or.&
         & (size(heap%keys)<heap%last))) then
      write(iout,*) 'Inconsistent size/allocation status!!'
    else    if ((heap%last > 0).and.((.not.allocated(heap%idxs)).or.&
         & (size(heap%idxs)<heap%last))) then
      write(iout,*) 'Inconsistent size/allocation status!!'
    else
      write(iout,*) heap%keys(1:heap%last)
      write(iout,*) heap%idxs(1:heap%last)
    end if
  end subroutine psb_d_idx_dump_heap

  subroutine psb_d_idx_free_heap(heap,info)
    implicit none
    class(psb_d_idx_heap), intent(inout) :: heap
    integer(psb_ipk_), intent(out)           :: info

    info=psb_success_
    if (allocated(heap%keys)) deallocate(heap%keys,stat=info)
    if ((info == psb_success_).and.(allocated(heap%idxs))) deallocate(heap%idxs,stat=info)

  end subroutine psb_d_idx_free_heap

end module psb_d_sort_mod


text:
program Example
  implicit none

  real :: ra(5,10)
  integer :: ia(5,10)
  integer :: i, j

  call random_number(ra)
  ia = int(ra * 20.0) + 1

outer: do i = 1, size(ia, 1)
         do j = 1, size(ia, 2)
           write(*, "(i3)", advance="no") ia(i,j)
           if (ia(i,j) == 20) exit outer
         end do
         write(*,*)
       end do outer

end program Example

max_issues_repo_name: seanwallawalla-forks/RosettaCodeData

max_forks_repo_name: seanwallawalla-forks/RosettaCodeData


text:
According to some, cow tipping is just an Urban Legends urban legend and is not actually physically feasible. A http://emntserver.unl.edu/MechanicsPages/MattSemke/The%20Statics%20of%20Cow%20Tipping.htm loose mathematical interpretation says it is possible if at least two people are doing the pushing. http://www.timesonline.co.uk/printFriendly/0,,1218582462,00.html This article about the physics of cow tipping also confirms that at least two people are needed to exert enough force.

Just because it may be possible to do it, however, doesnt mean it will happen. Cows do not sleep standing up, nor do their knees lock like a horses would. They do however doze. While traditional parts of the Cow Tipping folklore may not be true, there is anecdotal evidence suggesting that you can actually push a cow over  er, go cow tipping  with a small group of people if you are quiet about it. It physically should be possible if you are fast enough to do so before the cow reacts (requiring more people), or if the cow stands there rigidly and lets you.

UC Davis Chancellor Larry Vanderhoef was emailed about the rumored antitipping regulation by Users/KrisFricke on 200612 and received the following reply: UCD Student Affairs Student Affairs UCD Administrative Officers Vice Chancellor Judy Sakaki will be following up on this and will get back to you.  In the meantime, so far as I know there is no cow tipping policy or rule on campus.  At least I have not run across it during my twenty years here. An Urban Legends urban legend...?

Users/JabberWokky has seen it done in Belle Glade or western Lake Worth in Florida.  It also involved a cow that was awake, heavy rope and a pickup truck.  It was done by a group of people specifically to prove that it can be done (they did it a few times).  The cow went over, got up pissed, and seemed unharmed.  Despite the apparent lack of physical harm, it was certainly abuse.   Despite the people who insist that there is proof that the official way to do it is scientifically impossible (is there an official way to eat an apple or carry a book?), there is an easy way to do it that is actually practiced in rural Florida.  It may be touted impossible as a theoretical activity, but empirical evidence trumps even the best math.  Users/MichaelGiardina, on the other hand, has a video of former roommates who were less successful. Lets just say the bulls didnt like the idea and chased his friends all about.

Users/AllisonEriksen may have the best statement on the issue with her feeling that cows should be Tipping tipped like anyone else would be: with 15% of the cost of the service provided.

Cow tipping is supposedly fatal to the cow upon which this cruel activity is rendered, as the cows internal organs are not designed to withstand the sudden impact of a sideways fall. Cows are also the third most deadly animal in the United States, with only horses and deer killing more people each year.  Proceed with caution, or better yet, dont proceed at all with this tipping business...unless it involves the aforementioned Eriksen suggestion, minus the actual physical violence.

Yes, I know, Im as shocked as you are  I feel like my whole life has been a lie. Also Id heard that cow tipping is specifically an expellable offense here at UC Davis, from a source I thought was reliable (maybe my official UC Davis tour as a prospective student four years ago?)  now I want to confirm whether or not this rule is true in light of the impossibility of its realisation. (I think Im gonna email the chancellor on this) Users/KrisFricke Kris Fricke

Cowtipping is widely considered uncool in this town. It is detrimental to the cows and to the humans involved. I think my very first night in the dorms as a freshman, someone asked the RA about how prevalent cowtipping was here, and the gave them a firm admonishment. I have never heard of it happening or considered doing it since then. Users/KayRide

I think Cow Tipping being an urban legend is an urban legend.  As for it being cause for expulsion at Davis, I hear so is believing that the word gulliable is in the dictionary.  Users/JonerikStorm

Hate to ruin the fun, but cow tipping is an urban myth/legend. You cant tip a 1300 pound animal over, especially since shes going to move before she tips. :( Users/JessicaRuygt ImR
       Contrary (or in addition) to what Users/JabberWokky said above, the times that Ive heard about it being done, in the Davis Midwest Connection midwest or what have you, is when the animal is asleep (and thus not conscious enough to run), and when you have a whole gang of people to push the cow at one time.  Of course, this is complete heresay.  Anyhow, I would think it is feasibly possible, although still totally uncool.  Users/KalenRidenour

It should also be noted that unlike horses, cows knees do not lock when they sleep. Consequently, sleeping cows are not the immobilized hulks commonly pictured in cow tipping folklore. Users/KrisFricke

In the act known as cow tipping it is understood that the cow is asleep.  Cows do NOT sleep standing up and there is no way one is going to just tip over when it can step away when being pushed. It is an urban legend, made up so that people who know better can watch more gullible people get chased around in a pen full of cow pies. I grew up in dairy country and everyone knew that this was bunk.  Users/AllisonEriksen

Is there any real reason the cow has to be asleep? Why couldnt you just ram it with your (hopefully very inexpensive) car? As long as youre practicing immorality (how would you like it if a cow tipped YOU over, especially if you weighed 3/4 of a ton?) you might as well go whole hog. Im no physicist, but it seems to me that an oldsmobile traveling 35 mph could topple a cow. Im not endorsing this, by the way, just mentioning the possibility.  Users/CadeGrunst
    Thus crossing the line between tipping and senselessly ramming into an animal with your car.  Unless your car is very high off the ground, youll only injure the cow and get charged with animal cruelty.   Users/AllisonEriksen

Why is there such obsession with tipping over a cow? Why not a moose or elephant? Perhaps buffalo tipping? Walrus rumping? eh... Users/MichaelGiardina

20051205 14:41:58 nbsp Yes, it can be done. My professors have done it. Users/MichaelGiardina

20060216 22:04:41 nbsp Ive tried to tip my friends cows and its pretty hard as they tend to wake up/hear you/get angry we were chased and had to jump over a fence.  Besides, theres always the danger the cow will decide to run you down like that one mother did in Point Reyes!! Death by Cow sounds, other than highly ironic, rather painful. Users/LucyFitzGibbon

20070124 05:11:35 nbsp Cow tipping is pretty awesome :D Users/GreggAlexander

20070308 16:25:31 nbsp Hilarious, I grew up in Texas amid a cornicopia of cowtipping stories told by people with mental capacities lower than their bovine targets.  It s so funny how in spite of common sense and psyhics people still say it can happen cause their friend, told them they did it once Users/JarrettNoble

20070308 17:25:59 nbsp Cows and horses actually make a lot of sense as causes of mortality, but deer?  Surely not MOREso than cows and horses? Users/KrisFricke

20070308 17:41:26 nbsp deer through the windshield? Either hitting them and having them go through the windshield, or people swerving to avoid them maybe.  Deer are srs bizness. Users/Edwins ES

20070308 17:45:57 nbsp Yeah I almost swerved to miss a deer once.  After I ran it over I felt bad.  But not as bad as I would of if I had crashed my car. Users/JarrettNoble

20071011 02:57:20 nbsp When I first came to Davis, and took the campus tour, I asked about cow tipping. Sadly the guide only had heard of such things, from a friend of a friend, it saddened me, but I am glad I came here none the less. Anyone want to try? Users/DavidPoole

20071107 11:38:57 nbsp On the news right now is a cow who tipped itself off a cliff in Washington, landing on a family filled minivan. Nobody was killed, but the cow was euthanized on the scene. Users/GregFurnish
  Skywalker Ranch had a problem with that for awhile... a pasture had an overhang by a road, and they had some narrow misses when some of the cows started falling over the side of the cliff.  Im not sure it counts as tipping so much as plummeting.  You do have to wonder what a cow thinks when it falls on a family and their minivan: Mmmm... grass... boy do I love grass... its pretty much everything I love.  Tasty tasty grass... mmm... I... (slip) wait... theres no grass here... just a sort of whistling sound and things below me that are getting larger very quic... ARUGH!!! FOR THE LOVE OF GHOD AND ALL THAT IS HOLY, THAT HURT!!!  Gah!  Even worse, theres nothing but broken glass (not grass), twisted metal (not grass) and some seriously bewildered people (still not grass).  This sucks.  Hey, look over there... some more grass.  (meander over) Mmmm... grass... I love grass... grass tastes soooo awesome....  Users/JabberWokky

20080815 20:24:33 nbsp If you think being a plummeting cow is bad, imagine being a plummeting sperm whale: http://www.youtube.com/watch?vSzHuBYEqhqE

Ive never tried to tip, or otherwise displace a cow,  But I have attempted to move horses sideways on many occasions (not to tip them  I tip them with carrots!)  One of two things happens when you try to do this:

1) The horse stays right where it is, and possibly gives you a funny look

2) The horse steps over an inch or two and continues doing what it was doing.  It may also give you a funny look. Users/IDoNotExist
  My comment is pretty much entirely a reference to the plummeting sperm whale in HHGttG; the book reads much like the passage I wrote including (if I recall correctly) the whistling sound reference.  Only in my version, the cow survives and is (mentally) unscathed.  :) jw

20100118 14:16:06 nbsp Cow tipping is not a myth.  I am from rural Florida.  You take a group, obviously.  The cows are not asleep, but they are dozing.  You choose cows close to the fence or you will be chased down and injured.  This has never been described as a solitary action, but a group event.  Even after you push the cows over, you are in danger of injury.  You must run fast to get away.. apparently some of you are missing the point that that is the fun of cow tipping.  (Well, as it is considered by the teenagers who do this sort of thing.)  Ive always heard that it is not good for the cows, and when the troublemakers come to school bragging about it, other students would admonish them.   Users/Cerebral

  So it is like a bull run then? Users/hankim

20100729 17:54:32 nbsp Regardless of whether cow tipping is actually possible or not, if the University is so against it, why do they sell those obnoxious UC Davis Cow Tipping Team shirts and mugs in the bookstore? Doesnt really help our image... Users/Shellesen


text:
! KGEN-generated Fortran source file
!
! Filename    : rrtmg_sw_spcvmc.f90
! Generated at: 2015-07-07 00:48:25
! KGEN version: 0.4.13

    MODULE rrtmg_sw_spcvmc
        USE kgen_utils_mod, ONLY : kgen_dp, check_t, kgen_init_check, kgen_print_check
        !  --------------------------------------------------------------------------
        ! |                                                                          |
        ! |  Copyright 2002-2007, Atmospheric & Environmental Research, Inc. (AER).  |
        ! |  This software may be used, copied, or redistributed as long as it is    |
        ! |  not sold and this copyright notice is reproduced on each copy made.     |
        ! |  This model is provided as is without any express or implied warranties. |
        ! |                       (http://www.rtweb.aer.com/)                        |
        ! |                                                                          |
        !  --------------------------------------------------------------------------
        ! ------- Modules -------
        USE shr_kind_mod, ONLY: r8 => shr_kind_r8
        !      use parkind, only : jpim, jprb
        USE parrrsw, ONLY: ngptsw
        USE rrsw_tbl, ONLY: od_lo
        USE rrsw_tbl, ONLY: bpade
        USE rrsw_tbl, ONLY: tblint
        USE rrsw_tbl, ONLY: exp_tbl
        USE rrsw_wvn, ONLY: ngc
        USE rrsw_wvn, ONLY: ngs
        USE rrtmg_sw_reftra, ONLY: reftra_sw
        USE rrtmg_sw_taumol, ONLY: taumol_sw
        USE rrtmg_sw_vrtqdr, ONLY: vrtqdr_sw
        IMPLICIT NONE
        CONTAINS

        ! write subroutines
        ! No subroutines
        ! No module extern variables
        ! ---------------------------------------------------------------------------

        SUBROUTINE spcvmc_sw(lchnk, ncol, nlayers, istart, iend, icpr, idelm, iout, pavel, tavel, pz, tz, tbound, palbd, palbp, &
        pcldfmc, ptaucmc, pasycmc, pomgcmc, ptaormc, ptaua, pasya, pomga, prmu0, coldry, wkl, adjflux, laytrop, layswtch, laylow, &
        jp, jt, jt1, co2mult, colch4, colco2, colh2o, colmol, coln2o, colo2, colo3, fac00, fac01, fac10, fac11, selffac, selffrac,&
         indself, forfac, forfrac, indfor, pbbfd, pbbfu, pbbcd, pbbcu, puvfd, puvcd, pnifd, pnicd, pnifu, pnicu, pbbfddir, &
        pbbcddir, puvfddir, puvcddir, pnifddir, pnicddir, pbbfsu, pbbfsd)
            ! ---------------------------------------------------------------------------
            !
            ! Purpose: Contains spectral loop to compute the shortwave radiative fluxes,
            !          using the two-stream method of H. Barker and McICA, the Monte-Carlo
            !          Independent Column Approximation, for the representation of
            !          sub-grid cloud variability (i.e. cloud overlap).
            !
            ! Interface:  *spcvmc_sw* is called from *rrtmg_sw.F90* or rrtmg_sw.1col.F90*
            !
            ! Method:
            !    Adapted from two-stream model of H. Barker;
            !    Two-stream model options (selected with kmodts in rrtmg_sw_reftra.F90):
            !        1: Eddington, 2: PIFM, Zdunkowski et al., 3: discret ordinates
            !
            ! Modifications:
            !
            ! Original: H. Barker
            ! Revision: Merge with RRTMG_SW: J.-J.Morcrette, ECMWF, Feb 2003
            ! Revision: Add adjustment for Earth/Sun distance : MJIacono, AER, Oct 2003
            ! Revision: Bug fix for use of PALBP and PALBD: MJIacono, AER, Nov 2003
            ! Revision: Bug fix to apply delta scaling to clear sky: AER, Dec 2004
            ! Revision: Code modified so that delta scaling is not done in cloudy profiles
            !           if routine cldprop is used; delta scaling can be applied by swithcing
            !           code below if cldprop is not used to get cloud properties.
            !           AER, Jan 2005
            ! Revision: Modified to use McICA: MJIacono, AER, Nov 2005
            ! Revision: Uniform formatting for RRTMG: MJIacono, AER, Jul 2006
            ! Revision: Use exponential lookup table for transmittance: MJIacono, AER,
            !           Aug 2007
            !
            ! ------------------------------------------------------------------
            ! ------- Declarations ------
            ! ------- Input -------
            INTEGER, intent(in) :: lchnk
            INTEGER, intent(in) :: nlayers
            INTEGER, intent(in) :: istart
            INTEGER, intent(in) :: iend
            INTEGER, intent(in) :: icpr
            INTEGER, intent(in) :: idelm ! delta-m scaling flag
            ! [0 = direct and diffuse fluxes are unscaled]
            ! [1 = direct and diffuse fluxes are scaled]
            INTEGER, intent(in) :: iout
            INTEGER, intent(in) :: ncol ! column loop index
            INTEGER, intent(in) :: laytrop(ncol)
            INTEGER, intent(in) :: layswtch(ncol)
            INTEGER, intent(in) :: laylow(ncol)
            INTEGER, intent(in) :: indfor(:,:)
            !   Dimensions: (ncol,nlayers)
            INTEGER, intent(in) :: indself(:,:)
            !   Dimensions: (ncol,nlayers)
            INTEGER, intent(in) :: jp(:,:)
            !   Dimensions: (ncol,nlayers)
            INTEGER, intent(in) :: jt(:,:)
            !   Dimensions: (ncol,nlayers)
            INTEGER, intent(in) :: jt1(:,:)
            !   Dimensions: (ncol,nlayers)
            REAL(KIND=r8), intent(in) :: pavel(:,:) ! layer pressure (hPa, mb)
            !   Dimensions: (ncol,nlayers)
            REAL(KIND=r8), intent(in) :: tavel(:,:) ! layer temperature (K)
            !   Dimensions: (ncol,nlayers)
            REAL(KIND=r8), intent(in) :: pz(:,0:) ! level (interface) pressure (hPa, mb)
            !   Dimensions: (ncol,0:nlayers)
            REAL(KIND=r8), intent(in) :: tz(:,0:) ! level temperatures (hPa, mb)
            !   Dimensions: (ncol,0:nlayers)
            REAL(KIND=r8), intent(in) :: tbound(ncol) ! surface temperature (K)
            REAL(KIND=r8), intent(in) :: wkl(:,:,:) ! molecular amounts (mol/cm2)
            !   Dimensions: (ncol,mxmol,nlayers)
            REAL(KIND=r8), intent(in) :: coldry(:,:) ! dry air column density (mol/cm2)
            !   Dimensions: (ncol,nlayers)
            REAL(KIND=r8), intent(in) :: colmol(:,:)
            !   Dimensions: (ncol,nlayers)
            REAL(KIND=r8), intent(in) :: adjflux(:,:) ! Earth/Sun distance adjustment
            !   Dimensions: (ncol,jpband)
            REAL(KIND=r8), intent(in) :: palbd(:,:) ! surface albedo (diffuse)
            !   Dimensions: (ncol,nbndsw)
            REAL(KIND=r8), intent(in) :: palbp(:,:) ! surface albedo (direct)
            !   Dimensions: (ncol, nbndsw)
            REAL(KIND=r8), intent(in) :: prmu0(ncol) ! cosine of solar zenith angle
            REAL(KIND=r8), intent(in) :: pcldfmc(:,:,:) ! cloud fraction [mcica]
            !   Dimensions: (ncol,nlayers,ngptsw)
            REAL(KIND=r8), intent(in) :: ptaucmc(:,:,:) ! cloud optical depth [mcica]
            !   Dimensions: (ncol,nlayers,ngptsw)
            REAL(KIND=r8), intent(in) :: pasycmc(:,:,:) ! cloud asymmetry parameter [mcica]
            !   Dimensions: (ncol,nlayers,ngptsw)
            REAL(KIND=r8), intent(in) :: pomgcmc(:,:,:) ! cloud single scattering albedo [mcica]
            !   Dimensions: (ncol,nlayers,ngptsw)
            REAL(KIND=r8), intent(in) :: ptaormc(:,:,:) ! cloud optical depth, non-delta scaled [mcica]
            !   Dimensions: (ncol,nlayers,ngptsw)
            REAL(KIND=r8), intent(in) :: ptaua(:,:,:) ! aerosol optical depth
            !   Dimensions: (ncol,nlayers,nbndsw)
            REAL(KIND=r8), intent(in) :: pasya(:,:,:) ! aerosol asymmetry parameter
            !   Dimensions: (ncol,nlayers,nbndsw)
            REAL(KIND=r8), intent(in) :: pomga(:,:,:) ! aerosol single scattering albedo
            !   Dimensions: (ncol,nlayers,nbndsw)
            REAL(KIND=r8), intent(in) :: colh2o(:,:)
            !   Dimensions: (ncol,nlayers)
            REAL(KIND=r8), intent(in) :: colco2(:,:)
            !   Dimensions: (ncol,nlayers)
            REAL(KIND=r8), intent(in) :: colch4(:,:)
            !   Dimensions: (ncol,nlayers)
            REAL(KIND=r8), intent(in) :: co2mult(:,:)
            !   Dimensions: (ncol,nlayers)
            REAL(KIND=r8), intent(in) :: colo3(:,:)
            !   Dimensions: (ncol,nlayers)
            REAL(KIND=r8), intent(in) :: colo2(:,:)
            !   Dimensions: (ncol,nlayers)
            REAL(KIND=r8), intent(in) :: coln2o(:,:)
            !   Dimensions: (ncol,nlayers)
            REAL(KIND=r8), intent(in) :: forfac(:,:)
            !   Dimensions: (ncol,nlayers)
            REAL(KIND=r8), intent(in) :: forfrac(:,:)
            !   Dimensions: (ncol,nlayers)
            REAL(KIND=r8), intent(in) :: selffac(:,:)
            !   Dimensions: (ncol,nlayers)
            REAL(KIND=r8), intent(in) :: selffrac(:,:)
            !   Dimensions: (ncol,nlayers)
            REAL(KIND=r8), intent(in) :: fac00(:,:)
            !   Dimensions: (ncol,nlayers)
            REAL(KIND=r8), intent(in) :: fac01(:,:)
            !   Dimensions: (ncol,nlayers)
            REAL(KIND=r8), intent(in) :: fac10(:,:)
            !   Dimensions: (ncol,nlayers)
            REAL(KIND=r8), intent(in) :: fac11(:,:)
            !   Dimensions: (ncol,nlayers)
            ! ------- Output -------
            !   All Dimensions: (nlayers+1)
            REAL(KIND=r8), intent(out) :: pbbcd(:,:)
            REAL(KIND=r8), intent(out) :: pbbcu(:,:)
            REAL(KIND=r8), intent(out) :: pbbfd(:,:)
            REAL(KIND=r8), intent(out) :: pbbfu(:,:)
            REAL(KIND=r8), intent(out) :: pbbfddir(:,:)
            REAL(KIND=r8), intent(out) :: pbbcddir(:,:)
            REAL(KIND=r8), intent(out) :: puvcd(:,:)
            REAL(KIND=r8), intent(out) :: puvfd(:,:)
            REAL(KIND=r8), intent(out) :: puvcddir(:,:)
            REAL(KIND=r8), intent(out) :: puvfddir(:,:)
            REAL(KIND=r8), intent(out) :: pnicd(:,:)
            REAL(KIND=r8), intent(out) :: pnifd(:,:)
            REAL(KIND=r8), intent(out) :: pnicddir(:,:)
            REAL(KIND=r8), intent(out) :: pnifddir(:,:)
            ! Added for net near-IR flux diagnostic
            REAL(KIND=r8), intent(out) :: pnicu(:,:)
            REAL(KIND=r8), intent(out) :: pnifu(:,:)
            ! Output - inactive                                              !   All Dimensions: (nlayers+1)
            !      real(kind=r8), intent(out) :: puvcu(:)
            !      real(kind=r8), intent(out) :: puvfu(:)
            !      real(kind=r8), intent(out) :: pvscd(:)
            !      real(kind=r8), intent(out) :: pvscu(:)
            !      real(kind=r8), intent(out) :: pvsfd(:)
            !      real(kind=r8), intent(out) :: pvsfu(:)
            REAL(KIND=r8), intent(out) :: pbbfsu(:,:,:) ! shortwave spectral flux up (nswbands,nlayers+1)
            REAL(KIND=r8), intent(out) :: pbbfsd(:,:,:) ! shortwave spectral flux down (nswbands,nlayers+1)
            ! ------- Local -------
            LOGICAL :: lrtchkclr(ncol,nlayers)
            LOGICAL :: lrtchkcld(ncol,nlayers)
            INTEGER :: klev
            INTEGER :: ib1
            INTEGER :: ib2
            INTEGER :: ibm
            INTEGER :: igt
            INTEGER :: ikl
            INTEGER :: iw(ncol)
            INTEGER :: jk
            INTEGER :: jb
            INTEGER :: jg, iplon
            !      integer, parameter :: nuv = ??
            !      integer, parameter :: nvs = ??
            INTEGER :: itind(ncol)
            REAL(KIND=r8) :: ze1(ncol)
            REAL(KIND=r8) :: tblind(ncol)
            REAL(KIND=r8) :: zclear(ncol)
            REAL(KIND=r8) :: zcloud(ncol)
            REAL(KIND=r8) :: zdbt(ncol,nlayers+1)
            REAL(KIND=r8) :: zdbt_nodel(ncol,nlayers+1)
            REAL(KIND=r8) :: zgcc(ncol,nlayers)
            REAL(KIND=r8) :: zgco(ncol,nlayers)
            REAL(KIND=r8) :: zomcc(ncol,nlayers)
            REAL(KIND=r8) :: zomco(ncol,nlayers)
            REAL(KIND=r8) :: zrdndc(ncol,nlayers+1)
            REAL(KIND=r8) :: zrdnd(ncol,nlayers+1)
            REAL(KIND=r8) :: zrefc(ncol,nlayers+1)
            REAL(KIND=r8) :: zrefo(ncol,nlayers+1)
            REAL(KIND=r8) :: zref(	ncol,nlayers+1)
            REAL(KIND=r8) :: zrefdc(ncol,nlayers+1)
            REAL(KIND=r8) :: zrefdo(ncol,nlayers+1)
            REAL(KIND=r8) :: zrefd(ncol,nlayers+1)
            REAL(KIND=r8) :: zrup(ncol,nlayers+1)
            REAL(KIND=r8) :: zrupd(ncol,nlayers+1)
            REAL(KIND=r8) :: zrupc(ncol,nlayers+1)
            REAL(KIND=r8) :: zrupdc(ncol,nlayers+1)
            REAL(KIND=r8) :: ztauc(ncol,nlayers)
            REAL(KIND=r8) :: ztauo(ncol,nlayers)
            REAL(KIND=r8) :: ztdbt(ncol,nlayers+1)
            REAL(KIND=r8) :: ztrac(ncol,nlayers+1)
            REAL(KIND=r8) :: ztrao(ncol,nlayers+1)
            REAL(KIND=r8) :: ztra(ncol,nlayers+1)
            REAL(KIND=r8) :: ztradc(ncol,nlayers+1)
            REAL(KIND=r8) :: ztrado(ncol,nlayers+1)
            REAL(KIND=r8) :: ztrad(ncol,nlayers+1)
            REAL(KIND=r8) :: ztdbtc(ncol,nlayers+1)
            REAL(KIND=r8) :: zdbtc(ncol,nlayers+1)
            REAL(KIND=r8) :: zincflx(ncol,ngptsw)
            REAL(KIND=r8) :: zdbtc_nodel(ncol,nlayers+1)
            REAL(KIND=r8) :: ztdbtc_nodel(ncol,nlayers+1)
            REAL(KIND=r8) :: ztdbt_nodel(ncol,nlayers+1)
            REAL(KIND=r8) :: zdbtmc(ncol)
            REAL(KIND=r8) :: zdbtmo(ncol)
            REAL(KIND=r8) :: zf
            REAL(KIND=r8) :: repclc(ncol)
            REAL(KIND=r8) :: tauorig(ncol)
            REAL(KIND=r8) :: zwf
            !     real(kind=r8) :: zincflux                                   ! inactive
            ! Arrays from rrtmg_sw_taumoln routines
            !      real(kind=r8) :: ztaug(nlayers,16), ztaur(nlayers,16)
            !      real(kind=r8) :: zsflxzen(16)
            REAL(KIND=r8) :: ztaug(ncol,nlayers,ngptsw)
            REAL(KIND=r8) :: ztaur(ncol,nlayers,ngptsw)
            REAL(KIND=r8) :: zsflxzen(ncol,ngptsw)
            ! Arrays from rrtmg_sw_vrtqdr routine
            REAL(KIND=r8) :: zcd(ncol,nlayers+1,ngptsw)
            REAL(KIND=r8) :: zcu(ncol,nlayers+1,ngptsw)
            REAL(KIND=r8) :: zfd(ncol,nlayers+1,ngptsw)
            REAL(KIND=r8) :: zfu(ncol,nlayers+1,ngptsw)
            ! Inactive arrays
            !     real(kind=r8) :: zbbcd(nlayers+1), zbbcu(nlayers+1)
            !     real(kind=r8) :: zbbfd(nlayers+1), zbbfu(nlayers+1)
            !     real(kind=r8) :: zbbfddir(nlayers+1), zbbcddir(nlayers+1)
            ! ------------------------------------------------------------------
            ! Initializations
      ib1 = istart
      ib2 = iend
      klev = nlayers
      !djp repclc(iplon) = 1.e-12_r8
      repclc(:) = 1.e-12_r8
            !      zincflux = 0.0_r8
        do iplon=1,ncol
      do jk=1,klev+1
         pbbcd(iplon,jk)=0._r8
         pbbcu(iplon,jk)=0._r8
         pbbfd(iplon,jk)=0._r8
         pbbfu(iplon,jk)=0._r8
         pbbcddir(iplon,jk)=0._r8
         pbbfddir(iplon,jk)=0._r8
         puvcd(iplon,jk)=0._r8
         puvfd(iplon,jk)=0._r8
         puvcddir(iplon,jk)=0._r8
         puvfddir(iplon,jk)=0._r8
         pnicd(iplon,jk)=0._r8
         pnifd(iplon,jk)=0._r8
         pnicddir(iplon,jk)=0._r8
         pnifddir(iplon,jk)=0._r8
         pnicu(iplon,jk)=0._r8
         pnifu(iplon,jk)=0._r8
      enddo
        end do
        call taumol_sw(ncol,klev, &
                     colh2o, colco2, colch4, colo2, colo3, colmol, &
                     laytrop, jp, jt, jt1, &
                     fac00, fac01, fac10, fac11, &
                     selffac, selffrac, indself, forfac, forfrac,indfor, &
                     zsflxzen, ztaug, ztaur)

      jb = ib1-1                  ! ??? ! ???
			do iplon=1,ncol
					iw(iplon) =0
			end do
         do jb = ib1, ib2
			           ibm = jb-15
         igt = ngc(ibm)
                ! Reinitialize g-point counter for each band if output for each band is requested.
                !        do jk=1,klev+1
                !           zbbcd(jk)=0.0_r8
                !           zbbcu(jk)=0.0_r8
                !           zbbfd(jk)=0.0_r8
                !           zbbfu(jk)=0.0_r8
                !        enddo
                ! Top of g-point interval loop within each band (iw(iplon) is cumulative counter)

					DO IPLON=1,ncol
						if (iout.gt.0.and.ibm.ge.2) iw(iplon)= ngs(ibm-1)
					END do
         do jg = 1,igt
				  do iplon=1,ncol

	            iw(iplon) = iw(iplon)+1
                    ! Apply adjustment for correct Earth/Sun distance and zenith angle to incoming solar flux
            zincflx(iplon,iw(iplon)) = adjflux(iplon,jb) * zsflxzen(iplon,iw(iplon)) * prmu0(iplon)
                    !             zincflux = zincflux + adjflux(jb) * zsflxzen(iw(iplon)) * prmu0           ! inactive
                    ! Compute layer reflectances and transmittances for direct and diffuse sources,
                    ! first clear then cloudy
                    ! zrefc(iplon,jk)  direct albedo for clear
                    ! zrefo(iplon,jk)  direct albedo for cloud
                    ! zrefdc(iplon,jk) diffuse albedo for clear
                    ! zrefdo(iplon,jk) diffuse albedo for cloud
                    ! ztrac(iplon,jk)  direct transmittance for clear
                    ! ztrao(iplon,jk)  direct transmittance for cloudy
                    ! ztradc(iplon,jk) diffuse transmittance for clear
                    ! ztrado(iplon,jk) diffuse transmittance for cloudy
                    !
                    ! zref(iplon,jk)   direct reflectance
                    ! zrefd(iplon,jk)  diffuse reflectance
                    ! ztra(iplon,jk)   direct transmittance
                    ! ztrad(iplon,jk)  diffuse transmittance
                    !
                    ! zdbtc(iplon,jk)  clear direct beam transmittance
                    ! zdbto(jk)  cloudy direct beam transmittance
                    ! zdbt(iplon,jk)   layer mean direct beam transmittance
                    ! ztdbt(iplon,jk)  total direct beam transmittance at levels
                    ! Clear-sky
                    !   TOA direct beam
            ztdbtc(iplon,1)=1.0_r8
            ztdbtc_nodel(iplon,1)=1.0_r8
                    !   Surface values
            zdbtc(iplon,klev+1) =0.0_r8
            ztrac(iplon,klev+1) =0.0_r8
            ztradc(iplon,klev+1)=0.0_r8
            zrefc(iplon,klev+1) =palbp(iplon,ibm)
            zrefdc(iplon,klev+1)=palbd(iplon,ibm)
            zrupc(iplon,klev+1) =palbp(iplon,ibm)
            zrupdc(iplon,klev+1)=palbd(iplon,ibm)
                    ! Cloudy-sky
                    !   Surface values
            ztrao(iplon,klev+1) =0.0_r8
            ztrado(iplon,klev+1)=0.0_r8
            zrefo(iplon,klev+1) =palbp(iplon,ibm)
            zrefdo(iplon,klev+1)=palbd(iplon,ibm)
                    ! Total sky
                    !   TOA direct beam
            ztdbt(iplon,1)=1.0_r8
            ztdbt_nodel(iplon,1)=1.0_r8
                    !   Surface values
            zdbt(iplon,klev+1) =0.0_r8
            ztra(iplon,klev+1) =0.0_r8
            ztrad(iplon,klev+1)=0.0_r8
            zref(iplon,klev+1) =palbp(iplon,ibm)
            zrefd(iplon,klev+1)=palbd(iplon,ibm)
            zrup(iplon,klev+1) =palbp(iplon,ibm)
            zrupd(iplon,klev+1)=palbd(iplon,ibm)
                    ! Top of layer loop
            do jk=1,klev
                        ! Note: two-stream calculations proceed from top to bottom;
                        !   RRTMG_SW quantities are given bottom to top and are reversed here
               ikl=klev+1-jk
                        ! Set logical flag to do REFTRA calculation
                        !   Do REFTRA for all clear layers
               lrtchkclr(iplon,jk)=.true.
                        !   Do REFTRA only for cloudy layers in profile, since already done for clear layers
               lrtchkcld(iplon,jk)=.false.
               lrtchkcld(iplon,jk)=(pcldfmc(iplon,ikl,iw(iplon)) > repclc(iplon))
                        ! Clear-sky optical parameters - this section inactive
                        !   Original
                        !               ztauc(iplon,jk) = ztaur(ikl,iw(iplon)) + ztaug(ikl,iw(iplon))
                        !               zomcc(iplon,jk) = ztaur(ikl,iw(iplon)) / ztauc(iplon,jk)
                        !               zgcc(iplon,jk) = 0.0001_r8
                        !   Total sky optical parameters
                        !               ztauo(iplon,jk) = ztaur(ikl,iw(iplon)) + ztaug(ikl,iw(iplon)) + ptaucmc(ikl,iw(iplon))
                        !               zomco(iplon,jk) = ptaucmc(ikl,iw(iplon)) * pomgcmc(ikl,iw(iplon)) + ztaur(ikl,iw(iplon))
                        !               zgco (jk) = (ptaucmc(ikl,iw(iplon)) * pomgcmc(ikl,iw(iplon)) * pasycmc(ikl,iw(iplon)) + &
                        !                           ztaur(ikl,iw(iplon)) * 0.0001_r8) / zomco(iplon,jk)
                        !               zomco(iplon,jk) = zomco(iplon,jk) / ztauo(iplon,jk)
                        ! Clear-sky optical parameters including aerosols
               ztauc(iplon,jk) = ztaur(iplon,ikl,iw(iplon)) + ztaug(iplon,ikl,iw(iplon)) + ptaua(iplon,ikl,ibm)
               zomcc(iplon,jk) = ztaur(iplon,ikl,iw(iplon)) * 1.0_r8 + ptaua(iplon,ikl,ibm) * pomga(iplon,ikl,ibm)
               zgcc(iplon,jk) = pasya(iplon,ikl,ibm) * pomga(iplon,ikl,ibm) * ptaua(iplon,ikl,ibm) / zomcc(iplon,jk)
               zomcc(iplon,jk) = zomcc(iplon,jk) / ztauc(iplon,jk)
                        ! Pre-delta-scaling clear and cloudy direct beam transmittance (must use 'orig', unscaled cloud OD)
                        !   \/\/\/ This block of code is only needed for unscaled direct beam calculation
               if (idelm .eq. 0) then
                            !
                  zclear(iplon) = 1.0_r8 - pcldfmc(iplon,ikl,iw(iplon))
                  zcloud(iplon) = pcldfmc(iplon,ikl,iw(iplon))
                            ! Clear
                            !                   zdbtmc(iplon) = exp(-ztauc(iplon,jk) / prmu0)
                            ! Use exponential lookup table for transmittance, or expansion of exponential for low tau
                  ze1(iplon) = ztauc(iplon,jk) / prmu0(iplon)
                  if (ze1(iplon) .le. od_lo) then
                     zdbtmc(iplon) = 1._r8 - ze1(iplon) + 0.5_r8 * ze1(iplon) * ze1(iplon)
                  else
                     tblind(iplon) = ze1(iplon) / (bpade + ze1(iplon))
                     itind(iplon) = tblint * tblind(iplon) + 0.5_r8
                     zdbtmc(iplon) = exp_tbl(itind(iplon))
                  endif
                  zdbtc_nodel(iplon,jk) = zdbtmc(iplon)
                  ztdbtc_nodel(iplon,jk+1) = zdbtc_nodel(iplon,jk) * ztdbtc_nodel(iplon,jk)
                            ! Clear + Cloud
                  tauorig(iplon) = ztauc(iplon,jk) + ptaormc(iplon,ikl,iw(iplon))
                            !                   zdbtmo(iplon) = exp(-tauorig(iplon) / prmu0)
                            ! Use exponential lookup table for transmittance, or expansion of exponential for low tau
                  ze1(iplon) = tauorig(iplon) / prmu0(iplon)
                  if (ze1(iplon) .le. od_lo) then
                     zdbtmo(iplon) = 1._r8 - ze1(iplon) + 0.5_r8 * ze1(iplon) * ze1(iplon)
                  else
                     tblind(iplon) = ze1(iplon) / (bpade + ze1(iplon))
                     itind(iplon) = tblint * tblind(iplon) + 0.5_r8
                     zdbtmo(iplon) = exp_tbl(itind(iplon))
                  endif
                  zdbt_nodel(iplon,jk) = zclear(iplon)*zdbtmc(iplon) + zcloud(iplon)*zdbtmo(iplon)
                  ztdbt_nodel(iplon,jk+1) = zdbt_nodel(iplon,jk) * ztdbt_nodel(iplon,jk)
               endif
                        !   /\/\/\ Above code only needed for unscaled direct beam calculation
                        ! Delta scaling - clear
               zf = zgcc(iplon,jk) * zgcc(iplon,jk)
               zwf = zomcc(iplon,jk) * zf
               ztauc(iplon,jk) = (1.0_r8 - zwf) * ztauc(iplon,jk)
               zomcc(iplon,jk) = (zomcc(iplon,jk) - zwf) / (1.0_r8 - zwf)
               zgcc (iplon,jk) = (zgcc(iplon,jk) - zf) / (1.0_r8 - zf)
                        ! Total sky optical parameters (cloud properties already delta-scaled)
                        !   Use this code if cloud properties are derived in rrtmg_sw_cldprop
               if (icpr .ge. 1) then
                  ztauo(iplon,jk) = ztauc(iplon,jk) + ptaucmc(iplon,ikl,iw(iplon))
                  zomco(iplon,jk) = ztauc(iplon,jk) * zomcc(iplon,jk) + ptaucmc(iplon,ikl,iw(iplon)) * pomgcmc(iplon,ikl,iw(iplon))
                  zgco (iplon,jk) = (ptaucmc(iplon,ikl,iw(iplon)) * pomgcmc(iplon,ikl,iw(iplon)) * pasycmc(iplon,ikl,iw(iplon)) + &
                              ztauc(iplon,jk) * zomcc(iplon,jk) * zgcc(iplon,jk)) / zomco(iplon,jk)
                  zomco(iplon,jk) = zomco(iplon,jk) / ztauo(iplon,jk)
                            ! Total sky optical parameters (if cloud properties not delta scaled)
                            !   Use this code if cloud properties are not derived in rrtmg_sw_cldprop
               elseif (icpr .eq. 0) then
                  ztauo(iplon,jk) = ztaur(iplon,ikl,iw(iplon)) + ztaug(iplon,ikl,iw(iplon)) + ptaua(iplon,ikl,ibm) + ptaucmc(iplon,ikl,iw(iplon))
                  zomco(iplon,jk) = ptaua(iplon,ikl,ibm) * pomga(iplon,ikl,ibm) + ptaucmc(iplon,ikl,iw(iplon)) * pomgcmc(iplon,ikl,iw(iplon)) + &
                              ztaur(iplon,ikl,iw(iplon)) * 1.0_r8
                  zgco (iplon,jk) = (ptaucmc(iplon,ikl,iw(iplon)) * pomgcmc(iplon,ikl,iw(iplon)) * pasycmc(iplon,ikl,iw(iplon)) + &
                              ptaua(iplon,ikl,ibm)*pomga(iplon,ikl,ibm)*pasya(iplon,ikl,ibm)) / zomco(iplon,jk)
                  zomco(iplon,jk) = zomco(iplon,jk) / ztauo(iplon,jk)
                            ! Delta scaling - clouds
                            !   Use only if subroutine rrtmg_sw_cldprop is not used to get cloud properties and to apply delta scaling
                  zf = zgco(iplon,jk) * zgco(iplon,jk)
                  zwf = zomco(iplon,jk) * zf
                  ztauo(iplon,jk) = (1._r8 - zwf) * ztauo(iplon,jk)
                  zomco(iplon,jk) = (zomco(iplon,jk) - zwf) / (1.0_r8 - zwf)
                  zgco (iplon,jk) = (zgco(iplon,jk) - zf) / (1.0_r8 - zf)
               endif
                        ! End of layer loop
            enddo
				END DO
				DO iplon=1,ncol

                    ! Clear sky reflectivities
            call reftra_sw (klev,ncol, &
lrtchkclr, zgcc, prmu0, ztauc, zomcc, &
zrefc, zrefdc, ztrac, ztradc)
                    ! Total sky reflectivities
            call reftra_sw (klev, ncol, &
lrtchkcld, zgco, prmu0, ztauo, zomco, &
zrefo, zrefdo, ztrao, ztrado)
			END DO
				DO iplon=1,ncol
            do jk=1,klev
                        ! Combine clear and cloudy contributions for total sky
               ikl = klev+1-jk
               zclear(iplon) = 1.0_r8 - pcldfmc(iplon,ikl,iw(iplon))
               zcloud(iplon) = pcldfmc(iplon,ikl,iw(iplon))
               zref(iplon,jk) = zclear(iplon)*zrefc(iplon,jk) + zcloud(iplon)*zrefo(iplon,jk)
               zrefd(iplon,jk)= zclear(iplon)*zrefdc(iplon,jk) + zcloud(iplon)*zrefdo(iplon,jk)
               ztra(iplon,jk) = zclear(iplon)*ztrac(iplon,jk) + zcloud(iplon)*ztrao(iplon,jk)
               ztrad(iplon,jk)= zclear(iplon)*ztradc(iplon,jk) + zcloud(iplon)*ztrado(iplon,jk)
                        ! Direct beam transmittance
                        ! Clear
                        !                zdbtmc(iplon) = exp(-ztauc(iplon,jk) / prmu0)
                        ! Use exponential lookup table for transmittance, or expansion of
                        ! exponential for low tau
               ze1(iplon) = ztauc(iplon,jk) / prmu0(iplon)
               if (ze1(iplon) .le. od_lo) then
                  zdbtmc(iplon) = 1._r8 - ze1(iplon) + 0.5_r8 * ze1(iplon) * ze1(iplon)
               else
                  tblind(iplon) = ze1(iplon) / (bpade + ze1(iplon))
                  itind(iplon) = tblint * tblind(iplon) + 0.5_r8
                  zdbtmc(iplon) = exp_tbl(itind(iplon))
               endif
               zdbtc(iplon,jk) = zdbtmc(iplon)
               ztdbtc(iplon,jk+1) = zdbtc(iplon,jk)*ztdbtc(iplon,jk)
                        ! Clear + Cloud
                        !                zdbtmo(iplon) = exp(-ztauo(iplon,jk) / prmu0)
                        ! Use exponential lookup table for transmittance, or expansion of
                        ! exponential for low tau
               ze1(iplon) = ztauo(iplon,jk) / prmu0(iplon)
               if (ze1(iplon) .le. od_lo) then
                  zdbtmo(iplon) = 1._r8 - ze1(iplon) + 0.5_r8 * ze1(iplon) * ze1(iplon)
               else
                  tblind(iplon) = ze1(iplon) / (bpade + ze1(iplon))
                  itind(iplon) = tblint * tblind(iplon) + 0.5_r8
                  zdbtmo(iplon) = exp_tbl(itind(iplon))
               endif
               zdbt(iplon,jk) = zclear(iplon)*zdbtmc(iplon) + zcloud(iplon)*zdbtmo(iplon)
               ztdbt(iplon,jk+1) = zdbt(iplon,jk)*ztdbt(iplon,jk)
            enddo
                    ! Vertical quadrature for clear-sky fluxes
					END DO
!					DO iplon=1,ncol
            		call vrtqdr_sw(ncol,klev, iw, &
zrefc, zrefdc, ztrac, ztradc, &
zdbtc, zrdndc, zrupc, zrupdc, ztdbtc, &
zcd, zcu)
                    ! Vertical quadrature for cloudy fluxes
            		call vrtqdr_sw(ncol,klev, iw, &
zref, zrefd, ztra, ztrad, &
zdbt, zrdnd, zrup, zrupd, ztdbt, &
zfd, zfu)
!					END DO
					DO iplon=1,ncol
								! Upwelling and downwelling fluxes at levels
                    !   Two-stream calculations go from top to bottom;
                    !   layer indexing is reversed to go bottom to top for output arrays
            do jk=1,klev+1
               ikl=klev+2-jk
                        ! Accumulate spectral fluxes over bands - inactive
                        !               zbbfu(ikl) = zbbfu(ikl) + zincflx(iplon,iw(iplon))*zfu(iplon,jk,iw(iplon))
                        !               zbbfd(ikl) = zbbfd(ikl) + zincflx(iplon,iw(iplon))*zfd(iplon,jk,iw(iplon))
                        !               zbbcu(ikl) = zbbcu(ikl) + zincflx(iplon,iw(iplon))*zcu(iplon,jk,iw(iplon))
                        !               zbbcd(ikl) = zbbcd(ikl) + zincflx(iplon,iw(iplon))*zcd(iplon,jk,iw(iplon))
                        !               zbbfddir(ikl) = zbbfddir(ikl) + zincflx(iplon,iw(iplon))*ztdbt_nodel(iplon,jk)
                        !               zbbcddir(ikl) = zbbcddir(ikl) + zincflx(iplon,iw(iplon))*ztdbtc_nodel(iplon,jk)
               pbbfsu(iplon,ibm,ikl) = pbbfsu(iplon,ibm,ikl) + zincflx(iplon,iw(iplon))*zfu(iplon,jk,iw(iplon))
               pbbfsd(iplon,ibm,ikl) = pbbfsd(iplon,ibm,ikl) + zincflx(iplon,iw(iplon))*zfd(iplon,jk,iw(iplon))
                        ! Accumulate spectral fluxes over whole spectrum
               pbbfu(iplon,ikl) = pbbfu(iplon,ikl) + zincflx(iplon,iw(iplon))*zfu(iplon,jk,iw(iplon))
               pbbfd(iplon,ikl) = pbbfd(iplon,ikl) +zincflx(iplon,iw(iplon))*zfd(iplon,jk,iw(iplon))
               pbbcu(iplon,ikl) = pbbcu(iplon,ikl) + zincflx(iplon,iw(iplon))*zcu(iplon,jk,iw(iplon))
               pbbcd(iplon,ikl) = pbbcd(iplon,ikl) + zincflx(iplon,iw(iplon))*zcd(iplon,jk,iw(iplon))
               if (idelm .eq. 0) then
                  pbbfddir(iplon,ikl) = pbbfddir(iplon,ikl) + zincflx(iplon,iw(iplon))*ztdbt_nodel(iplon,jk)
                  pbbcddir(iplon,ikl) = pbbcddir(iplon,ikl) + zincflx(iplon,iw(iplon))*ztdbtc_nodel(iplon,jk)
               elseif (idelm .eq. 1) then
                  pbbfddir(iplon,ikl) = pbbfddir(iplon,ikl) + zincflx(iplon,iw(iplon))*ztdbt(iplon,jk)
                  pbbcddir(iplon,ikl) = pbbcddir(iplon,ikl) + zincflx(iplon,iw(iplon))*ztdbtc(iplon,jk)
               endif
                        ! Accumulate direct fluxes for UV/visible bands
               if (ibm >= 10 .and. ibm <= 13) then
                  puvcd(iplon,ikl) = puvcd(iplon,ikl) + zincflx(iplon,iw(iplon))*zcd(iplon,jk,iw(iplon))
                  puvfd(iplon,ikl) = puvfd(iplon,ikl) + zincflx(iplon,iw(iplon))*zfd(iplon,jk,iw(iplon))
                  if (idelm .eq. 0) then
                     puvfddir(iplon,ikl) = puvfddir(iplon,ikl) + zincflx(iplon,iw(iplon))*ztdbt_nodel(iplon,jk)
                     puvcddir(iplon,ikl) = puvcddir(iplon,ikl) + zincflx(iplon,iw(iplon))*ztdbtc_nodel(iplon,jk)
                  elseif (idelm .eq. 1) then
                     puvfddir(iplon,ikl) = puvfddir(iplon,ikl) + zincflx(iplon,iw(iplon))*ztdbt(iplon,jk)
                     puvcddir(iplon,ikl) = puvcddir(iplon,ikl) + zincflx(iplon,iw(iplon))*ztdbtc(iplon,jk)
                  endif
                            ! band 9 is half-NearIR and half-Visible
               else if (ibm == 9) then
                  puvcd(iplon,ikl) = puvcd(iplon,ikl) + 0.5_r8*zincflx(iplon,iw(iplon))*zcd(iplon,jk,iw(iplon))
                  puvfd(iplon,ikl) = puvfd(iplon,ikl) + 0.5_r8*zincflx(iplon,iw(iplon))*zfd(iplon,jk,iw(iplon))
                  pnicd(iplon,ikl) = pnicd(iplon,ikl) + 0.5_r8*zincflx(iplon,iw(iplon))*zcd(iplon,jk,iw(iplon))
                  pnifd(iplon,ikl) = pnifd(iplon,ikl) + 0.5_r8*zincflx(iplon,iw(iplon))*zfd(iplon,jk,iw(iplon))
                  if (idelm .eq. 0) then
                     puvfddir(iplon,ikl) = puvfddir(iplon,ikl) + 0.5_r8*zincflx(iplon,iw(iplon))*ztdbt_nodel(iplon,jk)
                     puvcddir(iplon,ikl) = puvcddir(iplon,ikl) + 0.5_r8*zincflx(iplon,iw(iplon))*ztdbtc_nodel(iplon,jk)
                     pnifddir(iplon,ikl) = pnifddir(iplon,ikl) + 0.5_r8*zincflx(iplon,iw(iplon))*ztdbt_nodel(iplon,jk)
                     pnicddir(iplon,ikl) = pnicddir(iplon,ikl) + 0.5_r8*zincflx(iplon,iw(iplon))*ztdbtc_nodel(iplon,jk)
                  elseif (idelm .eq. 1) then
                     puvfddir(iplon,ikl) = puvfddir(iplon,ikl) + 0.5_r8*zincflx(iplon,iw(iplon))*ztdbt(iplon,jk)
                     puvcddir(iplon,ikl) = puvcddir(iplon,ikl) + 0.5_r8*zincflx(iplon,iw(iplon))*ztdbtc(iplon,jk)
                     pnifddir(iplon,ikl) = pnifddir(iplon,ikl) + 0.5_r8*zincflx(iplon,iw(iplon))*ztdbt(iplon,jk)
                     pnicddir(iplon,ikl) = pnicddir(iplon,ikl) + 0.5_r8*zincflx(iplon,iw(iplon))*ztdbtc(iplon,jk)
                  endif
                  pnicu(iplon,ikl) = pnicu(iplon,ikl) + 0.5_r8*zincflx(iplon,iw(iplon))*zcu(iplon,jk,iw(iplon))
                  pnifu(iplon,ikl) = pnifu(iplon,ikl) + 0.5_r8*zincflx(iplon,iw(iplon))*zfu(iplon,jk,iw(iplon))
                            ! Accumulate direct fluxes for near-IR bands
               else if (ibm == 14 .or. ibm <= 8) then
                  pnicd(iplon,ikl) = pnicd(iplon,ikl) + zincflx(iplon,iw(iplon))*zcd(iplon,jk,iw(iplon))
                  pnifd(iplon,ikl) = pnifd(iplon,ikl) + zincflx(iplon,iw(iplon))*zfd(iplon,jk,iw(iplon))
                  if (idelm .eq. 0) then
                     pnifddir(iplon,ikl) = pnifddir(iplon,ikl) + zincflx(iplon,iw(iplon))*ztdbt_nodel(iplon,jk)
                     pnicddir(iplon,ikl) = pnicddir(iplon,ikl) + zincflx(iplon,iw(iplon))*ztdbtc_nodel(iplon,jk)
                  elseif (idelm .eq. 1) then
                     pnifddir(iplon,ikl) = pnifddir(iplon,ikl) + zincflx(iplon,iw(iplon))*ztdbt(iplon,jk)
                     pnicddir(iplon,ikl) = pnicddir(iplon,ikl) + zincflx(iplon,iw(iplon))*ztdbtc(iplon,jk)
                  endif
                            ! Added for net near-IR flux diagnostic
                  pnicu(iplon,ikl) = pnicu(iplon,ikl) + zincflx(iplon,iw(iplon))*zcu(iplon,jk,iw(iplon))
                  pnifu(iplon,ikl) = pnifu(iplon,ikl) + zincflx(iplon,iw(iplon))*zfu(iplon,jk,iw(iplon))
               endif
            enddo
                    ! End loop on jg, g-point interval
         enddo
                ! End loop on jb, spectral band
      enddo
        end do
        END SUBROUTINE spcvmc_sw
    END MODULE rrtmg_sw_spcvmc

max_stars_repo_name: kammerdienerb/flang

max_issues_repo_name: kammerdienerb/flang

max_forks_repo_name: kammerdienerb/flang