text:
*DECK DIPREP
      SUBROUTINE DIPREP (NEQ, Y, RWORK, IA, JA, IPFLAG, F, JAC)
      EXTERNAL F, JAC
      INTEGER NEQ, IA, JA, IPFLAG
      DOUBLE PRECISION Y, RWORK
      DIMENSION NEQ(*), Y(*), RWORK(*), IA(*), JA(*)
      INTEGER IOWND, IOWNS,
     1   ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L,
     2   LYH, LEWT, LACOR, LSAVF, LWM, LIWM, METH, MITER,
     3   MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU
      INTEGER IPLOST, IESP, ISTATC, IYS, IBA, IBIAN, IBJAN, IBJGP,
     1   IPIAN, IPJAN, IPJGP, IPIGP, IPR, IPC, IPIC, IPISP, IPRSP, IPA,
     2   LENYH, LENYHM, LENWK, LREQ, LRAT, LREST, LWMIN, MOSS, MSBJ,
     3   NSLJ, NGP, NLU, NNZ, NSP, NZL, NZU
      DOUBLE PRECISION ROWNS,
     1   CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND
      DOUBLE PRECISION RLSS
      COMMON /DLS001/ ROWNS(209),
     1   CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND,
     2   IOWND(6), IOWNS(6),
     3   ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L,
     4   LYH, LEWT, LACOR, LSAVF, LWM, LIWM, METH, MITER,
     5   MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU
      COMMON /DLSS01/ RLSS(6),
     1   IPLOST, IESP, ISTATC, IYS, IBA, IBIAN, IBJAN, IBJGP,
     2   IPIAN, IPJAN, IPJGP, IPIGP, IPR, IPC, IPIC, IPISP, IPRSP, IPA,
     3   LENYH, LENYHM, LENWK, LREQ, LRAT, LREST, LWMIN, MOSS, MSBJ,
     4   NSLJ, NGP, NLU, NNZ, NSP, NZL, NZU
      INTEGER I, IMAX, LEWTN, LYHD, LYHN
C-----------------------------------------------------------------------
C This routine serves as an interface between the driver and
C Subroutine DPREP.  It is called only if MITER is 1 or 2.
C Tasks performed here are:
C  * call DPREP,
C  * reset the required WM segment length LENWK,
C  * move YH back to its final location (following WM in RWORK),
C  * reset pointers for YH, SAVF, EWT, and ACOR, and
C  * move EWT to its new position if ISTATE = 1.
C IPFLAG is an output error indication flag.  IPFLAG = 0 if there was
C no trouble, and IPFLAG is the value of the DPREP error flag IPPER
C if there was trouble in Subroutine DPREP.
C-----------------------------------------------------------------------
      IPFLAG = 0
C Call DPREP to do matrix preprocessing operations. --------------------
      CALL DPREP (NEQ, Y, RWORK(LYH), RWORK(LSAVF), RWORK(LEWT),
     1   RWORK(LACOR), IA, JA, RWORK(LWM), RWORK(LWM), IPFLAG, F, JAC)
      LENWK = MAX(LREQ,LWMIN)
      IF (IPFLAG .LT. 0) RETURN
C If DPREP was successful, move YH to end of required space for WM. ----
      LYHN = LWM + LENWK
      IF (LYHN .GT. LYH) RETURN
      LYHD = LYH - LYHN
      IF (LYHD .EQ. 0) GO TO 20
      IMAX = LYHN - 1 + LENYHM
      DO 10 I = LYHN,IMAX
 10     RWORK(I) = RWORK(I+LYHD)
      LYH = LYHN
C Reset pointers for SAVF, EWT, and ACOR. ------------------------------
 20   LSAVF = LYH + LENYH
      LEWTN = LSAVF + N
      LACOR = LEWTN + N
      IF (ISTATC .EQ. 3) GO TO 40
C If ISTATE = 1, move EWT (left) to its new position. ------------------
      IF (LEWTN .GT. LEWT) RETURN
      DO 30 I = 1,N
 30     RWORK(I+LEWTN-1) = RWORK(I+LEWT-1)
 40   LEWT = LEWTN
      RETURN
C----------------------- End of Subroutine DIPREP ----------------------
      END

max_stars_repo_name: sbwhitecap/clocc-hg

max_issues_repo_name: sbwhitecap/clocc-hg

max_forks_repo_name: sbwhitecap/clocc-hg


text:
Jump to Timeline #Navigation Navigation

  William Vere Cruess born

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


module h5mod
! https://support.hdfgroup.org/HDF5/doc/HL/RM_H5LT.html
use, intrinsic:: iso_fortran_env, only: dp=>real64
use h5lt, only : hid_t, hsize_t, &
   h5f_acc_trunc_f, &  ! overwrite existing file
   h5open_f, h5fcreate_f, h5ltmake_dataset_double_f, h5fclose_f

implicit none
private

public:: h5write

contains

subroutine h5write(filename,dsname,array)
! convenience function I made, for easy resource.
! it overwrites existing "filename"

  character(*), intent(in) :: filename,dsname
  real(dp), intent(in) :: array(:)

  integer(hid_t) :: fid
  integer :: ioerr
  !---------- initialize fortran interface
  call h5open_f(ioerr)
  if (ioerr /= 0) error stop 'could not open hdf5 library'
  !--- create file
  call h5fcreate_f(filename, h5f_acc_trunc_f,fid,ioerr)
  if (ioerr /= 0) error stop 'could not create file'
  !---------- write array to hdf5
  print *,shape(array)
  call h5ltmake_dataset_double_f(fid,dsname,rank(array),shape(array,hsize_t),array,ioerr)
  if (ioerr /= 0) error stop 'could not write data'
  !----------- close file
  call h5fclose_f(fid,ioerr)
  if (ioerr /= 0) error stop 'could not close file'

  print *,'created ',filename

end subroutine h5write

end module


text:
      subroutine box5(k1,k2,k3,k4,k5,k6,za,zb,app,apm,amp,amm)
      implicit none
      include 'constants.f'
      include 'masses.f'
      include 'sprods_com.f'
      include 'zprods_decl.f'
      double precision mtsq,s134,rat,mtsqps134
      integer k1,k2,k3,k4,k5,k6,i1,i2
      double complex app,apm,amp,amm,zab2,izab2,iza,izb
c--- statement functions
      iza(i1,i2)=cone/za(i1,i2)
      izb(i1,i2)=cone/zb(i1,i2)
      zab2(k1,k2,k3,k4)=za(k1,k2)*zb(k2,k4)+za(k1,k3)*zb(k3,k4)
      izab2(k1,k2,k3,k4)=cone/(za(k1,k2)*zb(k2,k4)+za(k1,k3)*zb(k3,k4))
c--- end statement functions

      mtsq=mt**2
      s134=s(k1,k3)+s(k1,k4)+s(k3,k4)
      rat=mtsq/s(k1,k2)
      mtsqps134=mtsq+s134
      app =  + iza(k1,k2)*mtsq * (  - 2.D0*za(k1,k3)*za(k1,k5)*zb(k1,k2
     &    )*zb(k1,k4)*zb(k1,k6) + 2.D0*za(k1,k3)*za(k2,k5)*zb(k1,k2)*
     &    zb(k1,k4)*zb(k2,k6)*rat + 2.D0*za(k1,k5)*za(k2,k3)*zb(k1,k2)*
     &    zb(k1,k6)*zb(k2,k4)*rat + 2.D0*za(k1,k5)*za(k3,k4)*zb(k1,k2)*
     &    zb(k1,k4)*zb(k4,k6) )
      app = app + iza(k1,k2)*izab2(k1,k3,k4,k2)*mtsq*mtsqps134 * ( za(
     &    k1,k3)*za(k1,k5)*zb(k1,k2)*zb(k1,k4)*zb(k2,k6) + za(k1,k3)*
     &    za(k1,k5)*zb(k1,k2)*zb(k1,k6)*zb(k2,k4) - za(k1,k5)*za(k3,k4)
     &    *zb(k1,k2)*zb(k2,k4)*zb(k4,k6) )
      app = app + iza(k1,k2)*izab2(k1,k3,k4,k2)**2*mtsq*mtsqps134**2
     &  * (  - za(k1,k3)*za(k1,k5)*zb(k1,k2)*zb(k2,k4)*zb(k2,k6) )
      app = app + iza(k1,k2)*izab2(k1,k3,k4,k2)*zab2(k2,k3,k4,k1)*mtsq
     &  * (  - 2.D0*za(k1,k3)*za(k1,k5)*zb(k1,k2)*zb(k2,k4)*zb(k2,k6)*
     &    rat )
      app = app + iza(k1,k2)*izab2(k2,k3,k4,k1)*mtsq*mtsqps134 * ( za(
     &    k1,k3)*za(k2,k5)*zb(k1,k2)*zb(k1,k4)*zb(k1,k6) + za(k1,k5)*
     &    za(k2,k3)*zb(k1,k2)*zb(k1,k4)*zb(k1,k6) - za(k2,k5)*za(k3,k4)
     &    *zb(k1,k2)*zb(k1,k4)*zb(k4,k6) )
      app = app + iza(k1,k2)*izab2(k2,k3,k4,k1)**2*mtsq*mtsqps134**2
     &  * (  - za(k2,k3)*za(k2,k5)*zb(k1,k2)*zb(k1,k4)*zb(k1,k6) )
      app = app + iza(k1,k2)*izab2(k2,k3,k4,k1)*zab2(k1,k3,k4,k2)*mtsq
     &  * (  - 2.D0*za(k2,k3)*za(k2,k5)*zb(k1,k2)*zb(k1,k4)*zb(k1,k6)*
     &    rat )
      amm =  + izb(k1,k2)*mtsq * ( 2.D0*za(k1,k2)*za(k1,k3)*za(k2,k5)*
     &    zb(k1,k4)*zb(k2,k6)*rat + 2.D0*za(k1,k2)*za(k1,k3)*za(k2,k5)*
     &    zb(k1,k6)*zb(k2,k4) + 2.D0*za(k1,k2)*za(k1,k5)*za(k2,k3)*zb(
     &    k1,k6)*zb(k2,k4)*rat - 2.D0*za(k1,k2)*za(k2,k5)*za(k3,k4)*zb(
     &    k2,k4)*zb(k4,k6) )
      amm = amm + izb(k1,k2)*izab2(k1,k3,k4,k2)*mtsq*mtsqps134 * ( za(
     &    k1,k2)*za(k1,k3)*za(k1,k5)*zb(k1,k6)*zb(k2,k4) - za(k1,k2)*
     &    za(k1,k3)*za(k2,k5)*zb(k2,k4)*zb(k2,k6) - za(k1,k2)*za(k1,k5)
     &    *za(k3,k4)*zb(k2,k4)*zb(k4,k6) )
      amm = amm + izb(k1,k2)*izab2(k1,k3,k4,k2)**2*mtsq*mtsqps134**2
     &  * (  - za(k1,k2)*za(k1,k3)*za(k1,k5)*zb(k2,k4)*zb(k2,k6) )
      amm = amm + izb(k1,k2)*izab2(k1,k3,k4,k2)*zab2(k2,k3,k4,k1)*mtsq
     &  * (  - 2.D0*za(k1,k2)*za(k1,k3)*za(k1,k5)*zb(k2,k4)*zb(k2,k6)*
     &    rat )
      amm = amm + izb(k1,k2)*izab2(k2,k3,k4,k1)*mtsq*mtsqps134 * ( za(
     &    k1,k2)*za(k1,k3)*za(k2,k5)*zb(k1,k4)*zb(k1,k6) - za(k1,k2)*
     &    za(k2,k3)*za(k2,k5)*zb(k1,k6)*zb(k2,k4) - za(k1,k2)*za(k2,k5)
     &    *za(k3,k4)*zb(k1,k4)*zb(k4,k6) )
      amm = amm + izb(k1,k2)*izab2(k2,k3,k4,k1)**2*mtsq*mtsqps134**2
     &  * (  - za(k1,k2)*za(k2,k3)*za(k2,k5)*zb(k1,k4)*zb(k1,k6) )
      amm = amm + izb(k1,k2)*izab2(k2,k3,k4,k1)*zab2(k1,k3,k4,k2)*mtsq
     &  * (  - 2.D0*za(k1,k2)*za(k2,k3)*za(k2,k5)*zb(k1,k4)*zb(k1,k6)*
     &    rat )
      apm =  + mtsq * (  - 2.D0*za(k2,k3)*za(k2,k5)*zb(k1,k4)*zb(k1,k6)
     & - 2.D0*za(k2,k3)*za(k2,k5)*zb(k1,k4)*zb(k1,k6)*rat )
      apm = apm + izab2(k1,k3,k4,k2)*mtsqps134 * ( za(k1,k2)*za(k1,k3)*
     &    za(k2,k5)*zb(k1,k2)*zb(k1,k4)*zb(k1,k6) - za(k1,k2)*za(k2,k5)
     &    *za(k3,k4)*zb(k1,k2)*zb(k1,k4)*zb(k4,k6) )
      apm = apm + izab2(k1,k3,k4,k2)*mtsq*mtsqps134 * (  - za(k1,k3)*
     &    za(k2,k5)*zb(k1,k4)*zb(k1,k6) - za(k1,k5)*za(k2,k3)*zb(k1,k4)
     &    *zb(k1,k6) + za(k2,k3)*za(k2,k5)*zb(k1,k6)*zb(k2,k4) + za(k2,
     &    k5)*za(k3,k4)*zb(k1,k4)*zb(k4,k6) )
      apm = apm + izab2(k1,k3,k4,k2)**2*mtsqps134**2 * ( za(k1,k2)*za(
     &    k1,k3)*za(k1,k5)*zb(k1,k2)*zb(k1,k4)*zb(k1,k6) - za(k1,k2)*
     &    za(k1,k3)*za(k2,k5)*zb(k1,k2)*zb(k1,k4)*zb(k2,k6) - za(k1,k2)
     &    *za(k1,k3)*za(k2,k5)*zb(k1,k2)*zb(k1,k6)*zb(k2,k4) - za(k1,k2
     &    )*za(k1,k5)*za(k3,k4)*zb(k1,k2)*zb(k1,k4)*zb(k4,k6) + za(k1,
     &    k2)*za(k2,k5)*za(k3,k4)*zb(k1,k2)*zb(k2,k4)*zb(k4,k6) )
      apm = apm + izab2(k1,k3,k4,k2)**2*mtsq*mtsqps134**2 * ( za(k1,k3)
     &    *za(k2,k5)*zb(k1,k4)*zb(k2,k6) + za(k1,k5)*za(k2,k3)*zb(k1,k6
     &    )*zb(k2,k4) )
      apm = apm + izab2(k1,k3,k4,k2)**3*mtsqps134**3 * (  - za(k1,k2)*
     &    za(k1,k3)*za(k1,k5)*zb(k1,k2)*zb(k1,k4)*zb(k2,k6) - za(k1,k2)
     &    *za(k1,k3)*za(k1,k5)*zb(k1,k2)*zb(k1,k6)*zb(k2,k4) + za(k1,k2
     &    )*za(k1,k3)*za(k2,k5)*zb(k1,k2)*zb(k2,k4)*zb(k2,k6) + za(k1,
     &    k2)*za(k1,k5)*za(k3,k4)*zb(k1,k2)*zb(k2,k4)*zb(k4,k6) )
      apm = apm + izab2(k1,k3,k4,k2)**4*mtsqps134**4 * ( za(k1,k2)*za(
     &    k1,k3)*za(k1,k5)*zb(k1,k2)*zb(k2,k4)*zb(k2,k6) )
      apm = apm + izab2(k1,k3,k4,k2)**3*zab2(k2,k3,k4,k1)*mtsq*
     & mtsqps134**2 * (  - 4.D0*za(k1,k3)*za(k1,k5)*zb(k2,k4)*zb(k2,k6)
     & )
      apm = apm + izab2(k1,k3,k4,k2)**2*zab2(k2,k3,k4,k1)*mtsq*
     & mtsqps134 * ( 3.D0*za(k1,k3)*za(k1,k5)*zb(k1,k4)*zb(k2,k6) + 3.D0
     &    *za(k1,k3)*za(k1,k5)*zb(k1,k6)*zb(k2,k4) - 3.D0*za(k1,k3)*za(
     &    k2,k5)*zb(k2,k4)*zb(k2,k6) - 3.D0*za(k1,k5)*za(k3,k4)*zb(k2,
     &    k4)*zb(k4,k6) )
      apm = apm + izab2(k1,k3,k4,k2)**2*zab2(k2,k3,k4,k1)**2*mtsq * (
     & - 2.D0*za(k1,k3)*za(k1,k5)*zb(k2,k4)*zb(k2,k6)*rat )
      apm = apm + izab2(k1,k3,k4,k2)*zab2(k2,k3,k4,k1)*mtsq * (  - 2.D0
     &    *za(k1,k3)*za(k1,k5)*zb(k1,k4)*zb(k1,k6) + 2.D0*za(k1,k3)*za(
     &    k2,k5)*zb(k1,k4)*zb(k2,k6) + 2.D0*za(k1,k3)*za(k2,k5)*zb(k1,
     &    k4)*zb(k2,k6)*rat + 2.D0*za(k1,k3)*za(k2,k5)*zb(k1,k6)*zb(k2,
     &    k4) + 2.D0*za(k1,k5)*za(k2,k3)*zb(k1,k6)*zb(k2,k4)*rat + 2.D0
     &    *za(k1,k5)*za(k3,k4)*zb(k1,k4)*zb(k4,k6) - 2.D0*za(k2,k5)*za(
     &    k3,k4)*zb(k2,k4)*zb(k4,k6) )
      amp =  + mtsq * (  - 2.D0*za(k1,k3)*za(k1,k5)*zb(k2,k4)*zb(k2,k6)
     &    *rat )
      amp = amp + izab2(k2,k3,k4,k1)*mtsq*mtsqps134 * (  - za(k1,k3)*
     &    za(k1,k5)*zb(k1,k6)*zb(k2,k4) + za(k1,k5)*za(k3,k4)*zb(k2,k4)
     &    *zb(k4,k6) )
      amp = amp + izab2(k2,k3,k4,k1)**2*mtsq*mtsqps134**2 * ( za(k1,k3)
     &    *za(k2,k5)*zb(k1,k4)*zb(k2,k6) + za(k1,k5)*za(k2,k3)*zb(k1,k6
     &    )*zb(k2,k4) )
      amp = amp + izab2(k2,k3,k4,k1)**3*mtsqps134**3 * (  - za(k1,k2)*
     &    za(k1,k3)*za(k2,k5)*zb(k1,k2)*zb(k1,k4)*zb(k1,k6) + za(k1,k2)
     &    *za(k2,k5)*za(k3,k4)*zb(k1,k2)*zb(k1,k4)*zb(k4,k6) )
      amp = amp + izab2(k2,k3,k4,k1)**4*mtsqps134**4 * ( za(k1,k2)*za(
     &    k2,k3)*za(k2,k5)*zb(k1,k2)*zb(k1,k4)*zb(k1,k6) )
      amp = amp + izab2(k2,k3,k4,k1)**3*zab2(k1,k3,k4,k2)*mtsq*
     & mtsqps134**2 * (  - 4.D0*za(k2,k3)*za(k2,k5)*zb(k1,k4)*zb(k1,k6))
      amp = amp + izab2(k2,k3,k4,k1)**2*zab2(k1,k3,k4,k2)*mtsq*
     & mtsqps134 * ( 3.D0*za(k1,k3)*za(k2,k5)*zb(k1,k4)*zb(k1,k6) - 3.D0
     &    *za(k2,k5)*za(k3,k4)*zb(k1,k4)*zb(k4,k6) )
      amp = amp + izab2(k2,k3,k4,k1)**2*zab2(k1,k3,k4,k2)**2*mtsq * (
     &     - 2.D0*za(k2,k3)*za(k2,k5)*zb(k1,k4)*zb(k1,k6)*rat )
      amp = amp + izab2(k2,k3,k4,k1)*zab2(k1,k3,k4,k2)*mtsq * ( 2.D0*
     &    za(k1,k3)*za(k2,k5)*zb(k1,k4)*zb(k2,k6)*rat + 2.D0*za(k1,k5)*
     &    za(k2,k3)*zb(k1,k6)*zb(k2,k4)*rat )
      app=app*(-half*im)/(s(k3,k4)*s(k5,k6))
      amm=amm*(-half*im)/(s(k3,k4)*s(k5,k6))
      apm=apm*(-half*im)/(s(k3,k4)*s(k5,k6))
      amp=amp*(-half*im)/(s(k3,k4)*s(k5,k6))
      return
      end

max_stars_repo_name: tmartini/JHUGen

max_issues_repo_name: tmartini/JHUGen

max_forks_repo_name: tmartini/JHUGen


text:
      subroutine get_arg(i,s)
      integer       i
      character*(*) s
      call getarg(i,s)
      return
      end

      subroutine n_arg(i)
      integer       i
      i=iarcg()
      return
      end

max_stars_repo_name: RWTH-OS/MP-MPICH

max_issues_repo_name: RWTH-OS/MP-MPICH

max_forks_repo_name: RWTH-OS/MP-MPICH


text:
module topography_mod

! <CONTACT EMAIL="Bruce.Wyman@noaa.gov">
!   Bruce Wyman
! </CONTACT>

! <HISTORY SRC="http://www.gfdl.noaa.gov/fms-cgi-bin/cvsweb.cgi/FMS/"/>

! <OVERVIEW>
!   Routines for creating land surface topography fields and land-water masks
!   for latitude-longitude grids.
! </OVERVIEW>

! <DESCRIPTION>
!   This module generates realistic mountains and land-water masks
!   on a specified latitude-longitude grid by interpolating from the
!   1/6 degree Navy mean topography and percent water data sets.
!   The fields that can be generated are mean and standard deviation
!   of topography within the specified grid boxes; and land-ocean (or
!   water) mask and land-ocean (or water) fractional area.
!
!   The interpolation scheme conserves the area-weighted average
!   of the input data by using module horiz_interp.
!
!   The interfaces get_gaussian_topog and gaussian_topog_init are documented in <LINK SRC="gaussian_topog.html">gaussian_topog_mod</LINK>.
! </DESCRIPTION>

use gaussian_topog_mod, only: gaussian_topog_init, get_gaussian_topog
use   horiz_interp_mod, only: horiz_interp

use            fms_mod, only: file_exist, check_nml_error,               &
                              open_namelist_file, close_file, stdlog,    &
                              mpp_pe, mpp_root_pe, write_version_number, &
                              open_ieee32_file, error_mesg, FATAL, NOTE, &
                              mpp_error
use         fms_io_mod, only: read_data
implicit none
private

public :: topography_init,                 &
          get_topog_mean, get_topog_stdev, &
          get_ocean_frac, get_ocean_mask,  &
          get_water_frac, get_water_mask,  &
          gaussian_topog_init, get_gaussian_topog

!-----------------------------------------------------------------------
! <NAMELIST NAME="topography_nml">
!   <DATA NAME="topog_file" TYPE="character" DEFAULT="DATA/navy_topography.data">
!       Name of topography file.
!   </DATA>
!   <DATA NAME="water_file" TYPE="character" DEFAULT="DATA/navy_pctwater.data">
!       Name of percent water file.
!   </DATA>

   character(len=128) :: topog_file = 'DATA/navy_topography.data', &
                         water_file = 'DATA/navy_pctwater.data'
   namelist /topography_nml/ topog_file, water_file
! </NAMELIST>

!-----------------------------------------------------------------------
! --- resolution of the topography data set ---
! <DATASET NAME="">
!   This module uses the 1/6 degree U.S. Navy mean topography
!   and percent water data sets.
!
!   These data sets have been re-formatted to separate 32-bit IEEE files.
!   The names of these files is specified by the <LINK SRC="#NAMELIST">namelist</LINK> input.
!
!The format for both files is as follows:
!

!     record = 1    nlon, nlat
!     record = 2    blon, blat
!     record = 3    data
!

!where:
!

!     nlon, nlat = The number of longitude and latitude points
!                  in the horizontal grid.  For the 1/6 degree
!                  data sets this is 2160 x 1080. [integer]
!     blon, blat = The longitude and latitude grid box boundaries in degrees.
!                     [real :: blon(nlon+1), blat(nlat+1)]
!
!     data       = The topography or percent water data.
!                    [real :: data(nlon,nlat)]
!

! </DATASET>
  integer :: unit
  integer :: ipts, jpts
  integer, parameter :: COMPUTE_STDEV = 123  ! use this flag to
                                             !   compute st dev

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

 character(len=128) :: version = '$Id: topography.f90,v 13.0 2006/03/28 21:43:31 fms Exp $'
 character(len=128) :: tagname = '$Name: latest $'

 logical :: module_is_initialized = .FALSE.

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

 contains

!#######################################################################

   subroutine topography_init ()

     if ( module_is_initialized ) return

     call write_version_number (version,tagname)
     call read_namelist
     module_is_initialized = .TRUE.

   end subroutine topography_init

!#######################################################################

! <FUNCTION NAME="get_topog_mean">

!   <OVERVIEW>
!     Returns a "realistic" mean surface height field.
!   </OVERVIEW>
!   <DESCRIPTION>
!     Returns realistic mountains on a latitude-longtude grid.
!     The returned field is the mean topography for the given grid boxes.
!     Computed using a conserving area-weighted interpolation.
!     The current input data set is the 1/6 degree Navy mean topography.
!   </DESCRIPTION>
!   <TEMPLATE>
!     flag = get_topog_mean ( blon, blat, zmean )
!   </TEMPLATE>

!   <IN NAME="blon" TYPE="real" DIM="(:)">
!     The longitude (in radians) at grid box boundaries.
!   </IN>
!   <IN NAME="blat" TYPE="real" DIM="(:)">
!     The latitude (in radians) at grid box boundaries.
!   </IN>
!   <OUT NAME="zmean" UNIT=" meter" TYPE="real" DIM="(:,:)">
!     The mean surface height (meters).
!     The size of this field must be size(blon)-1 by size(blat)-1.
!   </OUT>
!   <OUT NAME="get_topog_mean" TYPE="logical">
!     A logical value of TRUE is returned if the surface height field
!     was successfully created. A value of FALSE may be returned if the
!     input topography data set was not readable.
!   </OUT>

!   <ERROR MSG="shape(zmean) is not equal to (/size(blon)-1,size(blat)-1/))" STATUS="FATAL">
!     Check the input grid size and output field size.
!   </ERROR>

 function get_topog_mean (blon, blat, zmean)

   real, intent(in),  dimension(:)   :: blon, blat
   real, intent(out), dimension(:,:) :: zmean
   logical :: get_topog_mean

!-----------------------------------------------------------------------
   if (.not. module_is_initialized) call topography_init()

   if ( any(shape(zmean(:,:)) /= (/size(blon(:))-1,size(blat(:))-1/)) ) &
        call error_mesg('get_topog_mean','shape(zmean) is not&
            & equal to (/size(blon)-1,size(blat)-1/))', FATAL)

   if ( open_topog_file(topog_file) ) then
       call interp_topog ( blon, blat, zmean )
       get_topog_mean = .true.
   else
       get_topog_mean = .false.
   endif

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

 end function get_topog_mean
! </FUNCTION>

!#######################################################################

! <FUNCTION NAME="get_topog_stdev">

!   <OVERVIEW>
!     Returns a standard deviation of higher resolution topography with
!     the given model grid boxes.
!   </OVERVIEW>
!   <DESCRIPTION>
!     Returns the standard deviation of the "finer" input topography data set,
!     currently the Navy 1/6 degree mean topography data, within the
!     boundaries of the given input grid.
!   </DESCRIPTION>
!   <TEMPLATE>
!     flag = get_topog_stdev ( blon, blat, stdev )
!   </TEMPLATE>
!   <IN NAME="blon" TYPE="real" DIM="(:)">
!     The longitude (in radians) at grid box boundaries.
!   </IN>
!   <IN NAME="blat" TYPE="real" DIM="(:)">
!     The latitude (in radians) at grid box boundaries.
!   </IN>
!   <OUT NAME="stdev" UNITS="meter" TYPE="real" DIM="(:,:)">
!     The standard deviation of surface height (in meters) within
!     given input model grid boxes.
!     The size of this field must be size(blon)-1 by size(blat)-1.
!   </OUT>
!   <OUT NAME="get_topog_stdev" TYPE="logical">
!     A logical value of TRUE is returned if the output field was
!     successfully created. A value of FALSE may be returned if the
!     input topography data set was not readable.
!   </OUT>

 function get_topog_stdev (blon, blat, stdev)

   real, intent(in),  dimension(:)   :: blon, blat
   real, intent(out), dimension(:,:) :: stdev
   logical :: get_topog_stdev

   real :: zsurf (size(stdev,1),size(stdev,2))
!-----------------------------------------------------------------------
   if (.not. module_is_initialized) call topography_init()

   if ( any(shape(stdev(:,:)) /= (/size(blon(:))-1,size(blat(:))-1/)) ) &
       call error_mesg('get_topog_stdev','shape(stdev) is not&
            & equal to (/size(blon)-1,size(blat)-1/))', FATAL)

   if ( open_topog_file(topog_file) ) then
       call interp_topog ( blon, blat, stdev, flag=COMPUTE_STDEV )
       get_topog_stdev = .true.
   else
       get_topog_stdev = .false.
   endif

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

 end function get_topog_stdev
! </FUNCTION>

!#######################################################################

! <FUNCTION NAME="get_ocean_frac">

!   <OVERVIEW>
!      Returns fractional area covered by ocean in a grid box.
!   </OVERVIEW>
!   <DESCRIPTION>
!     Returns fractional area covered by ocean in the given model grid boxes.
!   </DESCRIPTION>
!   <TEMPLATE>
!     flag = get_ocean_frac ( blon, blat, ocean_frac )
!   </TEMPLATE>

!   <IN NAME="blon" UNITS="radians" TYPE="real" DIM="(:)">
!     The longitude (in radians) at grid box boundaries.
!   </IN>
!   <IN NAME="blat" UNITS="radians" TYPE="real" DIM="(:)">
!     The latitude (in radians) at grid box boundaries.
!   </IN>
!   <OUT NAME="ocean_frac" TYPE="real" DIM="(:,:)">
!     The fractional amount (0 to 1) of ocean in a grid box.
!     The size of this field must be size(blon)-1 by size(blat)-1.
!   </OUT>
!   <OUT NAME="get_ocean_frac" TYPE="logical">
!     A logical value of TRUE is returned if the output field
!     was successfully created. A value of FALSE may be returned
!     if the Navy 1/6 degree percent water data set was not readable.
!   </OUT>

 function get_ocean_frac (blon, blat, ocean_frac)

 real, intent(in),  dimension(:)   :: blon, blat
 real, intent(out), dimension(:,:) :: ocean_frac
 logical :: get_ocean_frac

!-----------------------------------------------------------------------
   if (.not. module_is_initialized) call topography_init()

   if ( any(shape(ocean_frac(:,:)) /= (/size(blon(:))-1,size(blat(:))-1/)) ) &
        call error_mesg('get_ocean_frac','shape(ocean_frac) is not&
                 & equal to (/size(blon)-1,size(blat)-1/))', FATAL)

   if ( open_topog_file(water_file) ) then
       call interp_water ( blon, blat, ocean_frac, do_ocean=.true. )
       get_ocean_frac = .true.
   else
       get_ocean_frac = .false.
   endif

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

 end function get_ocean_frac
! </FUNCTION>

!#######################################################################
! <FUNCTION NAME="get_ocean_mask">

!   <OVERVIEW>
!     Returns a land-ocean mask in a grid box.
!   </OVERVIEW>
!   <DESCRIPTION>
!     Returns a land-ocean mask in the given model grid boxes.
!   </DESCRIPTION>
!   <TEMPLATE>
!     flag = get_ocean_mask ( blon, blat, ocean_mask )
!   </TEMPLATE>

!   <IN NAME="blon" UNITS="radians" TYPE="real" DIM="(:)">
!     The longitude (in radians) at grid box boundaries.
!   </IN>
!   <IN NAME="blat" UNITS="radians" TYPE="real" DIM="(:)">
!     The latitude (in radians) at grid box boundaries.
!   </IN>
!   <OUT NAME="ocean_frac" TYPE="real" DIM="(:,:)">
!     The fractional amount (0 to 1) of ocean in a grid box.
!     The size of this field must be size(blon)-1 by size(blat)-1.
!   </OUT>
!   <OUT NAME="get_ocean_mask" TYPE="logical">
!     A logical value of TRUE is returned if the output field
!     was successfully created. A value of FALSE may be returned
!     if the Navy 1/6 degree percent water data set was not readable.
!   </OUT>

 function get_ocean_mask (blon, blat, ocean_mask)

 real   , intent(in),  dimension(:)   :: blon, blat
 logical, intent(out), dimension(:,:) :: ocean_mask
 logical :: get_ocean_mask

 real, dimension(size(ocean_mask,1),size(ocean_mask,2)) :: ocean_frac
!-----------------------------------------------------------------------
   if (.not. module_is_initialized) call topography_init()

 if ( get_ocean_frac(blon, blat, ocean_frac) ) then
   where (ocean_frac > 0.50)
     ocean_mask = .true.
   elsewhere
     ocean_mask = .false.
   end where
   get_ocean_mask = .true.
 else
   get_ocean_mask = .false.
 endif

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

 end function get_ocean_mask
! </FUNCTION>

!#######################################################################
! <FUNCTION NAME="get_water_frac">

!   <OVERVIEW>
!     Returns fractional area covered by water.
!   </OVERVIEW>
!   <DESCRIPTION>
!     Returns the percent of water in a grid box.
!   </DESCRIPTION>
!   <TEMPLATE>
!     flag = get_water_frac ( blon, blat, water_frac )
!   </TEMPLATE>

!   <IN NAME="blon" UNITS="radians" TYPE="real" DIM="(:)">
!     The longitude (in radians) at grid box boundaries.
!   </IN>
!   <IN NAME="blat" UNITS="radians" TYPE="real" DIM="(:)">
!     The latitude (in radians) at grid box boundaries.
!   </IN>
!   <OUT NAME="water_frac" TYPE="real" DIM="(:,:)">
!     The fractional amount (0 to 1) of water in a grid box.
!     The size of this field must be size(blon)-1 by size(blat)-1.
!   </OUT>
!   <OUT NAME="get_water_frac" TYPE="logical">
!     A logical value of TRUE is returned if the output field
!     was successfully created. A value of FALSE may be returned
!     if the Navy 1/6 degree percent water data set was not readable.
!   </OUT>
!   <ERROR MSG="shape(water_frac) is not equal to (/size(blon)-1,size(blat)-1/))" STATUS="FATAL">
!      Check the input grid size and output field size.
!   </ERROR>

 function get_water_frac (blon, blat, water_frac)

 real, intent(in),  dimension(:)   :: blon, blat
 real, intent(out), dimension(:,:) :: water_frac
 logical :: get_water_frac

!-----------------------------------------------------------------------
   if (.not. module_is_initialized) call topography_init()

   if ( any(shape(water_frac(:,:)) /= (/size(blon(:))-1,size(blat(:))-1/)) ) &
        call error_mesg('get_water_frac','shape(water_frac) is not&
                 & equal to (/size(blon)-1,size(blat)-1/))', FATAL)

   if ( open_topog_file(water_file) ) then
       call interp_water ( blon, blat, water_frac )
       get_water_frac = .true.
   else
       get_water_frac = .false.
   endif

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

 end function get_water_frac
! </FUNCTION>

!#######################################################################
! <FUNCTION NAME="get_water_mask">

!   <OVERVIEW>
!     Returns a land-water mask in a grid box.
!   </OVERVIEW>
!   <DESCRIPTION>
!     Returns a land-water mask in the given model grid boxes.
!   </DESCRIPTION>
!   <TEMPLATE>
!     flag = get_water_mask ( blon, blat, water_mask )
!   </TEMPLATE>

!   <IN NAME="blon" UNITS="radians" TYPE="real" DIM="(:)">
!     The longitude (in radians) at grid box boundaries.
!   </IN>
!   <IN NAME="blat" UNITS="radians" TYPE="real" DIM="(:)">
!     The latitude (in radians) at grid box boundaries.
!   </IN>
!   <OUT NAME="water_mask" TYPE="real" DIM="(:,:)">
!     A binary mask for water (true) or land (false).
!     The size of this field must be size(blon)-1 by size(blat)-1.
!   </OUT>
!   <OUT NAME="get_water_mask" TYPE="logical">
!     A logical value of TRUE is returned if the output field
!     was successfully created. A value of FALSE may be returned
!     if the Navy 1/6 degree percent water data set was not readable.
!   </OUT>

 function get_water_mask (blon, blat, water_mask)

 real   , intent(in),  dimension(:)   :: blon, blat
 logical, intent(out), dimension(:,:) :: water_mask
 logical :: get_water_mask

 real, dimension(size(water_mask,1),size(water_mask,2)) :: water_frac
!-----------------------------------------------------------------------
   if (.not. module_is_initialized) call topography_init()

 if ( get_water_frac(blon, blat, water_frac) ) then
   where (water_frac > 0.50)
     water_mask = .true.
   elsewhere
     water_mask = .false.
   end where
   get_water_mask = .true.
 else
   get_water_mask = .false.
 endif

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

 end function get_water_mask
! </FUNCTION>

!#######################################################################
!##################   private interfaces below here   ##################
!#######################################################################

 function open_topog_file ( filename )
 character(len=*), intent(in) :: filename
 logical :: open_topog_file
 real    :: r_ipts, r_jpts
 integer :: namelen

 namelen = len(trim(filename))
  if ( file_exist(filename) .AND. filename(namelen-2:namelen) == '.nc') then
     if (mpp_pe() == mpp_root_pe()) call mpp_error ('topography_mod', &
            'Reading NetCDF formatted input data file: '//filename, NOTE)
     call read_data(filename, 'ipts', r_ipts, no_domain=.true.)
     call read_data(filename, 'jpts', r_jpts, no_domain=.true.)
     ipts = nint(r_ipts)
     jpts = nint(r_jpts)
     open_topog_file = .true.
  else
     if ( file_exist(filename) ) then
        if (mpp_pe() == mpp_root_pe()) call mpp_error ('topography_mod', &
             'Reading native formatted input data file: '//filename, NOTE)
        unit = open_ieee32_file (trim(filename), 'read')
        read (unit) ipts, jpts
        open_topog_file = .true.
     else
        open_topog_file = .false.
     endif
  endif

 end function open_topog_file

!#######################################################################

 subroutine interp_topog ( blon, blat, zout, flag )
 real   , intent(in)  :: blon(:), blat(:)
 real   , intent(out) :: zout(:,:)
 integer, intent(in), optional :: flag

 real :: xdat(ipts+1), ydat(jpts+1)
 real :: zdat(ipts,jpts)
 real :: zout2(size(zout,1),size(zout,2))
 integer :: namelen

 namelen = len(trim(topog_file))

! note: ipts,jpts,unit are global

  if ( file_exist(topog_file) .AND. topog_file(namelen-2:namelen)=='.nc') then
     call read_data(topog_file, 'xdat', xdat, no_domain=.true.)
     call read_data(topog_file, 'ydat', ydat, no_domain=.true.)
     call read_data(topog_file, 'zdat', zdat, no_domain=.true.)
  else
    read (unit) xdat, ydat    ! read lon/lat edges in radians
    read (unit) zdat          ! read land surface height in meters
    call close_file (unit)
 endif

    call horiz_interp ( zdat, xdat, ydat, blon, blat, zout )

! compute standard deviation if necessary
    if (present(flag)) then
       if (flag == COMPUTE_STDEV) then
           zdat = zdat*zdat
           call horiz_interp ( zdat, xdat, ydat, blon, blat, zout2 )
           zout = zout2 - zout*zout
           where (zout > 0.0)
             zout = sqrt ( zout )
           elsewhere
             zout = 0.0
           endwhere
       endif
    endif

 end subroutine interp_topog

!#######################################################################

 subroutine interp_water ( blon, blat, zout, do_ocean )
 real   , intent(in)  :: blon(:), blat(:)
 real   , intent(out) :: zout(:,:)
 logical, intent(in), optional :: do_ocean

 real :: xdat(ipts+1), ydat(jpts+1), zdat(ipts,jpts)
 integer :: namelen

 namelen = len(trim(water_file))

! note: ipts,jpts,unit are global
   if( file_exist(trim(water_file)) .AND. water_file(namelen-2:namelen)=='.nc') then
      call read_data(water_file, 'xdat', xdat, no_domain=.true.)
      call read_data(water_file, 'ydat', ydat, no_domain=.true.)
      call read_data(water_file, 'zdat', zdat, no_domain=.true.)
   else
    read (unit) xdat, ydat    ! read lon/lat edges in radians
    read (unit) zdat          ! read fractional water
    call close_file (unit)
 endif
! only use designated ocean points
    if (present(do_ocean)) then
        if (do_ocean) call determine_ocean_points (zdat)
    endif

! interpolate onto output grid
    call horiz_interp ( zdat, xdat, ydat, blon, blat, zout )

 end subroutine interp_water

!#######################################################################

 subroutine determine_ocean_points ( pctwater )
 real, intent(inout) :: pctwater(:,:)
 logical :: ocean(size(pctwater,1),size(pctwater,2))
 integer :: i, j, m, n, im, ip, jm, jp, new

 real :: ocean_pct_crit = .500

  ! resolution of the grid
    m = size(pctwater,1)
    n = size(pctwater,2)

  ! the 1/6 degree navy percent water data set
  ! designates ocean grid boxes as 100 percent water
  ! all other grid boxes have <= 99 percent water

  ! set a mask for ocean grid boxes
    ocean = (pctwater > .999)
    new = count(ocean)

  ! set land grid boxes that have sufficient amount of water
  ! to ocean grid boxes when they are adjacent to ocean points
  ! iterate until there are no new ocean points
    do
    if (new == 0) exit
    new = 0

       do j = 1, n
       do i = 1, m
          if (.not.ocean(i,j) .and. pctwater(i,j) > ocean_pct_crit) then
             im = i-1; ip = i+1; jm = j-1; jp = j+1
             if (im == 0)   im = m
             if (ip == m+1) ip = 1
             if (jm == 0)   jm = 1
             if (jp == n+1) jp = n
           ! check the 8 grid boxes that surround this grid box
             if (ocean(im,j ) .or. ocean(ip,j ) .or. ocean(i ,jm) .or. ocean(i ,jp) .or. &
                 ocean(im,jm) .or. ocean(ip,jm) .or. ocean(ip,jp) .or. ocean(im,jp)) then
                 ocean(i,j) = .true.
                 new = new + 1
             endif
          endif
       enddo
       enddo
      !print *, 'new=',new

    enddo

  ! final step is to elimate water percentage if land
    where (.not.ocean) pctwater = 0.

 end subroutine determine_ocean_points

!#######################################################################
! reads the namelist file, write namelist to log file,
! and initializes constants

subroutine read_namelist

   integer :: unit, ierr, io
   real    :: dtr, ght

!  read namelist

   if ( file_exist('input.nml')) then
      unit = open_namelist_file ( )
      ierr=1; do while (ierr /= 0)
         read  (unit, nml=topography_nml, iostat=io, end=10)
         ierr = check_nml_error(io,'topography_nml')
      enddo
 10   call close_file (unit)
   endif

!  write version and namelist to log file

   if (mpp_pe() == mpp_root_pe()) then
     unit = stdlog()
     write (unit, nml=topography_nml)
   endif

end subroutine read_namelist

!#######################################################################

end module topography_mod

! <INFO>

!   <TESTPROGRAM NAME="">
!
!  To run this program you will need the topography and percent water
!  data sets and use the following namelist (in file input.nml).
!
!   &gaussian_topog_nml
!     height = 5000., 3000., 3000., 3000.,
!     olon   =   90.,  255.,  285.,    0.,
!     olat   =   45.,   45.,  -15.,  -90.,
!     wlon   =   15.,   10.,    5.,  180.,
!     wlat   =   15.,   25.,   25.,   20., /
!
!  program test
!
!  ! test program for topography and gaussian_topog modules
!

!  use topography_mod
!  implicit none
!
!  integer, parameter :: nlon=24, nlat=18
!  real :: x(nlon), y(nlat), xb(nlon+1), yb(nlat+1), z(nlon,nlat)
!  real :: hpi, rtd
!  integer :: i,j
!  logical :: a
!
!  ! gaussian mountain parameters
!  real, parameter :: ht=4000.
!  real, parameter :: x0=90., y0=45. ! origin in degrees
!  real, parameter :: xw=15., yw=15. ! half-width in degees
!  real, parameter :: xr=30., yr= 0. ! ridge-width in degrees
!
!  ! create lat/lon grid in radians
!    hpi = acos(0.0)
!    rtd = 90./hpi ! rad to deg
!    do i=1,nlon
!      xb(i) = 4.*hpi*real(i-1)/real(nlon)
!    enddo
!      xb(nlon+1) = xb(1)+4.*hpi
!      yb(1) = -hpi
!    do j=2,nlat
!      yb(j) = yb(j-1) + 2.*hpi/real(nlat)
!    enddo
!      yb(nlat+1) = hpi
!  ! mid-point of grid boxes
!    x(1:nlon) = 0.5*(xb(1:nlon)+xb(2:nlon+1))
!    y(1:nlat) = 0.5*(yb(1:nlat)+yb(2:nlat+1))
!  ! test topography_mod routines
!    a = get_topog_mean(xb,yb,z)
!    call printz ('get_topog_mean')
!
!    a = get_water_frac(xb,yb,z)
!    z = z*100. ! in percent
!    call printz ('get_water_frac')
!
!    a = get_ocean_frac(xb,yb,z)
!    z = z*100. ! in percent
!    call printz ('get_ocean_frac')
!
!  ! test gaussian_topog_mod routines
!    a = .true.
!    z = get_gaussian_topog(x,y,ht,x0,y0,xw,yw,xr,yr)
!    call printz ('get_gaussian_topog')
!
!    call gaussian_topog_init (x,y,z)
!    call printz ('gaussian_topog_init')
!
!  contains
!
!  ! simple printout of topog/water array
!    subroutine printz (lab)
!    character(len=*), intent(in) :: lab
!     if (a) then
!        print '(/a)', trim(lab)
!     else
!        print '(/a)', 'no data available: '//trim(lab)
!        return
!     endif
!      ! print full grid
!        print '(3x,25i5)', (nint(x(i)*rtd),i=1,nlon)
!      do j=nlat,1,-1
!        print '(i3,25i5)', nint(y(j)*rtd), (nint(z(i,j)),i=1,nlon)
!      enddo
!    end subroutine printz
!
!  end program test
!

!   </TESTPROGRAM>

!   <BUG>
!      Water mask produces some possible erroneous water points along
!      the coast of Antarctic (at about 90W).
!   </BUG>

!   <FUTURE>Use of netcdf data sets. </FUTURE>
!   <FUTURE>Incorporate other topography and ocean data sets. </FUTURE>
!
! </INFO>

max_stars_repo_name: JoyMonteiro/FMSDynamicalCore


        !COMPILER-GENERATED INTERFACE MODULE: Fri May 22 16:58:36 2020
        ! This source file is for reference only and may not completely
        ! represent the generated interface used by the compiler.
        MODULE COMPUTE_COUPLING_PORO_VISCOELASTIC__genmod
          INTERFACE
            SUBROUTINE COMPUTE_COUPLING_PORO_VISCOELASTIC
            END SUBROUTINE COMPUTE_COUPLING_PORO_VISCOELASTIC
          END INTERFACE
        END MODULE COMPUTE_COUPLING_PORO_VISCOELASTIC__genmod


module coarsen_coeffs_module

  use bl_constants_module
  use bl_types
  use multifab_module

  implicit none

  private
  public :: coarsen_cell_coeffs, coarsen_edge_coeffs

contains

  subroutine coarsen_cell_coeffs(cf,cc)

    type(multifab), intent(in   ) :: cf
    type(multifab), intent(inout) :: cc

    real(kind=dp_t), pointer :: cfp(:,:,:,:)
    real(kind=dp_t), pointer :: ccp(:,:,:,:)
    integer :: i,ng,dm
    integer :: lof(get_dim(cf))
    integer :: loc(get_dim(cc)),hic(get_dim(cc))

    dm = get_dim(cf)

    if ( ncomp(cf) .ne. ncomp(cc) ) then
       print *,'ncomp_fine not equal to ncomp_crse in coarsen_cell_coeffs'
       print *,'ncomp_fine = ',ncomp(cf)
       print *,'ncomp_crse = ',ncomp(cc)
       call bl_error("coarsen_coeffs.f90 :: coarsen_cell_coeffs")
    end if

    ng = nghost(cc)

    !$OMP PARALLEL DO PRIVATE(i,cfp,ccp,loc,hic,lof)
    do i = 1, nfabs(cf)
       cfp => dataptr(cf, i)
       ccp => dataptr(cc, i)

       loc =  lwb(get_box(cc, i))
       hic =  upb(get_box(cc, i))

       lof =  lwb(get_box(cf, i))

       select case (dm)
       case (1)
          call crse_cell_coeffs_1d(ccp(:,1,1,:), cfp(:,1,1,:), ng, loc, hic, lof)
       case (2)
          call crse_cell_coeffs_2d(ccp(:,:,1,:), cfp(:,:,1,:), ng, loc, hic, lof)
       case (3)
          call crse_cell_coeffs_3d(ccp(:,:,:,:), cfp(:,:,:,:), ng, loc, hic, lof)
       end select
    end do
    !$OMP END PARALLEL DO

  end subroutine coarsen_cell_coeffs

  subroutine coarsen_edge_coeffs(cf,cc)

    type(multifab), intent(in   ) :: cf(:)
    type(multifab), intent(inout) :: cc(:)

    real(kind=dp_t), pointer :: cxp(:,:,:,:)
    real(kind=dp_t), pointer :: cyp(:,:,:,:)
    real(kind=dp_t), pointer :: czp(:,:,:,:)
    real(kind=dp_t), pointer :: fxp(:,:,:,:)
    real(kind=dp_t), pointer :: fyp(:,:,:,:)
    real(kind=dp_t), pointer :: fzp(:,:,:,:)

    integer :: i,ng,dm
    integer :: lof(get_dim(cf(1)))
    integer :: loc(get_dim(cc(1))), hic(get_dim(cc(1)))

    dm = get_dim(cf(1))

    if ( ncomp(cf(1)) .ne. ncomp(cc(1)) ) then
       print *,'ncomp_fine not equal to ncomp_crse in coarsen_edge_coeffs'
       print *,'ncomp_fine = ',ncomp(cf(1))
       print *,'ncomp_crse = ',ncomp(cc(1))
       call bl_error("coarsen_coeffs.f90 :: coarsen_edge_coeffs")
    end if

    ng = nghost(cc(1))

    !$OMP PARALLEL DO PRIVATE(i,loc,hic,lof,cxp,cyp,czp,fxp,fyp,fzp)
    do i = 1, nfabs(cf(1))

       select case (dm)
       case (1)
          loc =  lwb(get_box(cc(1), i))
          hic =  upb(get_box(cc(1), i))
          lof =  lwb(get_box(cf(1), i))
          cxp => dataptr(cc(1), i)
          fxp => dataptr(cf(1), i)
          call crse_xedge_coeffs_1d(cxp(:,1,1,:), fxp(:,1,1,:), ng, loc, hic, lof)
       case (2)
          loc =  lwb(get_box(cc(1), i))
          hic =  upb(get_box(cc(1), i))
          lof =  lwb(get_box(cf(1), i))
          cxp => dataptr(cc(1), i)
          fxp => dataptr(cf(1), i)
          call crse_xedge_coeffs_2d(cxp(:,:,1,:), fxp(:,:,1,:), ng, loc, hic, lof)
          loc =  lwb(get_box(cc(2), i))
          hic =  upb(get_box(cc(2), i))
          lof =  lwb(get_box(cf(2), i))
          cyp => dataptr(cc(2), i)
          fyp => dataptr(cf(2), i)
          call crse_yedge_coeffs_2d(cyp(:,:,1,:), fyp(:,:,1,:), ng, loc, hic, lof)
       case (3)
          loc =  lwb(get_box(cc(1), i))
          hic =  upb(get_box(cc(1), i))
          lof =  lwb(get_box(cf(1), i))
          cxp => dataptr(cc(1), i)
          fxp => dataptr(cf(1), i)
          call crse_xedge_coeffs_3d(cxp, fxp, ng, loc, hic, lof)
          loc =  lwb(get_box(cc(2), i))
          hic =  upb(get_box(cc(2), i))
          lof =  lwb(get_box(cf(2), i))
          cyp => dataptr(cc(2), i)
          fyp => dataptr(cf(2), i)
          call crse_yedge_coeffs_3d(cyp, fyp, ng, loc, hic, lof)
          loc =  lwb(get_box(cc(3), i))
          hic =  upb(get_box(cc(3), i))
          lof =  lwb(get_box(cf(3), i))
          czp => dataptr(cc(3), i)
          fzp => dataptr(cf(3), i)
          call crse_zedge_coeffs_3d(czp, fzp, ng, loc, hic, lof)
       end select
    end do
    !$OMP END PARALLEL DO

  end subroutine coarsen_edge_coeffs

  subroutine crse_cell_coeffs_1d(cc, cf, ng, loc, hic, lof)

    integer        , intent(in   ) :: ng,lof(:),loc(:),hic(:)
    real(kind=dp_t), intent(inout) :: cc(loc(1)-ng:,:)
    real(kind=dp_t), intent(in)    :: cf(lof(1)-ng:,:)
    integer :: i, i2, n

    do n = 1, size(cc,2)
       do i = loc(1),hic(1)
          i2 = 2*i
          cc(i,n) = HALF * (cf(i2,n) + cf(i2+1,n))
       end do
    end do

  end subroutine crse_cell_coeffs_1d

  subroutine crse_xedge_coeffs_1d(cc, cf, ng, loc, hic, lof)

    integer        , intent(in   ) :: ng,lof(:),loc(:),hic(:)
    real(kind=dp_t), intent(inout) :: cc(loc(1)-ng:,:)
    real(kind=dp_t), intent(in)    :: cf(lof(1)-ng:,:)
    integer :: i, i2, n

    do n = 1, size(cc,2)
       do i = loc(1),hic(1)+1
          i2 = 2*i
          cc(i,n) = cf(i2,n)
       end do
    end do

  end subroutine crse_xedge_coeffs_1d

  subroutine crse_cell_coeffs_2d(cc, cf, ng, loc, hic, lof)

    integer        , intent(in   ) :: ng,lof(:),loc(:),hic(:)
    real(kind=dp_t), intent(inout) :: cc(loc(1)-ng:,loc(2)-ng:,:)
    real(kind=dp_t), intent(in)    :: cf(lof(1)-ng:,lof(2)-ng:,:)
    integer :: i, i2, j, j2, n

    do n = 1, size(cc,3)
       do j = loc(2),hic(2)
          do i = loc(1),hic(1)
             i2 = 2*i
             j2 = 2*j
             cc(i,j,n) = FOURTH * (  &
                  cf(i2,j2  ,n) + cf(i2+1,j2  ,n) + &
                  cf(i2,j2+1,n) + cf(i2+1,j2+1,n) )
          end do
       end do
    end do

  end subroutine crse_cell_coeffs_2d

  subroutine crse_xedge_coeffs_2d(cc, cf, ng, loc, hic, lof)

    integer        , intent(in   ) :: ng,lof(:),loc(:),hic(:)
    real(kind=dp_t), intent(inout) :: cc(loc(1)-ng:,loc(2)-ng:,:)
    real(kind=dp_t), intent(in)    :: cf(lof(1)-ng:,lof(2)-ng:,:)
    integer :: i, i2, j, j2, n

    do n = 1, size(cc,3)
       do j = loc(2),hic(2)
          do i = loc(1),hic(1)+1
             i2 = 2*i
             j2 = 2*j
             cc(i,j,n) = HALF * (cf(i2,j2,n) + cf(i2,j2+1,n))
          end do
       end do
    end do

  end subroutine crse_xedge_coeffs_2d

  subroutine crse_yedge_coeffs_2d(cc, cf, ng, loc, hic, lof)

    integer        , intent(in   ) :: ng,lof(:),loc(:),hic(:)
    real(kind=dp_t), intent(inout) :: cc(loc(1)-ng:,loc(2)-ng:,:)
    real(kind=dp_t), intent(in)    :: cf(lof(1)-ng:,lof(2)-ng:,:)
    integer :: i, i2, j, j2, n

    do n = 1, size(cc,3)
       do j = loc(2),hic(2)+1
          do i = loc(1),hic(1)
             i2 = 2*i
             j2 = 2*j
             cc(i,j,n) = HALF * (cf(i2,j2,n) + cf(i2+1,j2,n))
          end do
       end do
    end do

  end subroutine crse_yedge_coeffs_2d

  subroutine crse_cell_coeffs_3d(cc, cf, ng, loc, hic, lof)

    integer        , intent(in   ) :: ng,lof(:),loc(:),hic(:)
    real(kind=dp_t), intent(inout) :: cc(loc(1)-ng:,loc(2)-ng:,loc(3)-ng:,1:)
    real(kind=dp_t), intent(in)    :: cf(lof(1)-ng:,lof(2)-ng:,lof(3)-ng:,1:)
    integer :: i, i2, j, j2, k, k2, n

    do n = 1, size(cc,4)
       do k = loc(3),hic(3)
          k2 = 2*k
          do j = loc(2),hic(2)
             j2 = 2*j
             do i = loc(1),hic(1)
                i2 = 2*i
                cc(i,j,k,n) = EIGHTH * ( &
                     + cf(i2,j2  ,k2  ,n) + cf(i2+1,j2  ,k2  ,n) &
                     + cf(i2,j2+1,k2  ,n) + cf(i2+1,j2+1,k2  ,n) &
                     + cf(i2,j2  ,k2+1,n) + cf(i2+1,j2  ,k2+1,n) &
                     + cf(i2,j2+1,k2+1,n) + cf(i2+1,j2+1,k2+1,n) &
                     )
             end do
          end do
       end do
    end do

  end subroutine crse_cell_coeffs_3d

  subroutine crse_xedge_coeffs_3d(cc, cf, ng, loc, hic, lof)

    integer        , intent(in   ) :: ng,lof(:),loc(:),hic(:)
    real(kind=dp_t), intent(inout) :: cc(loc(1)-ng:,loc(2)-ng:,loc(3)-ng:,1:)
    real(kind=dp_t), intent(in)    :: cf(lof(1)-ng:,lof(2)-ng:,lof(3)-ng:,1:)
    integer :: i, i2, j, j2, k, k2, n

    do n = 1, size(cc,4)
       do k = loc(3),hic(3)
          k2 = 2*k
          do j = loc(2),hic(2)
             j2 = 2*j
             do i = loc(1),hic(1)+1
                i2 = 2*i
                cc(i,j,k,n) = FOURTH * ( &
                     + cf(i2,j2,k2  ,n) + cf(i2,j2+1,k2  ,n) &
                     + cf(i2,j2,k2+1,n) + cf(i2,j2+1,k2+1,n) &
                     )
             end do
          end do
       end do
    end do

  end subroutine crse_xedge_coeffs_3d

  subroutine crse_yedge_coeffs_3d(cc, cf, ng, loc, hic, lof)

    integer        , intent(in   ) :: ng,lof(:),loc(:),hic(:)
    real(kind=dp_t), intent(inout) :: cc(loc(1)-ng:,loc(2)-ng:,loc(3)-ng:,1:)
    real(kind=dp_t), intent(in)    :: cf(lof(1)-ng:,lof(2)-ng:,lof(3)-ng:,1:)
    integer :: i, i2, j, j2, k, k2, n

    do n = 1, size(cc,4)
       do k = loc(3),hic(3)
          k2 = 2*k
          do j = loc(2),hic(2)+1
             j2 = 2*j
             do i = loc(1),hic(1)
                i2 = 2*i
                cc(i,j,k,n) = FOURTH * ( &
                     + cf(i2,j2,k2  ,n) + cf(i2+1,j2,k2  ,n) &
                     + cf(i2,j2,k2+1,n) + cf(i2+1,j2,k2+1,n) &
                     )
             end do
          end do
       end do
    end do

  end subroutine crse_yedge_coeffs_3d

  subroutine crse_zedge_coeffs_3d(cc, cf, ng, loc, hic, lof)

    integer        , intent(in   ) :: ng,lof(:),loc(:),hic(:)
    real(kind=dp_t), intent(inout) :: cc(loc(1)-ng:,loc(2)-ng:,loc(3)-ng:,1:)
    real(kind=dp_t), intent(in)    :: cf(lof(1)-ng:,lof(2)-ng:,lof(3)-ng:,1:)
    integer :: i, i2, j, j2, k, k2, n

    do n = 1, size(cc,4)
       do k = loc(3),hic(3)+1
          k2 = 2*k
          do j = loc(2),hic(2)
             j2 = 2*j
             do i = loc(1),hic(1)
                i2 = 2*i
                cc(i,j,k,n) = FOURTH * ( &
                     + cf(i2,j2  ,k2,n) + cf(i2+1,j2  ,k2,n) &
                     + cf(i2,j2+1,k2,n) + cf(i2+1,j2+1,k2,n) &
                     )
             end do
          end do
       end do
    end do

  end subroutine crse_zedge_coeffs_3d

end module coarsen_coeffs_module


text:
!----------------------------------------------------------------------------- best with 100 columns

function readVTK1() result(ierr)
  use modFileIO
  use modPolyMesh
  integer ierr
  type(polyMesh)::mesh

  ierr=0
  open(10,file='data/mesh_simple.vtk',action='read')
  call readVTK(10,mesh)
  close(10)
  if(any(mesh%iNE/=reshape([1,3,2,1,4,3,1,2,4,2,3,4],shape(mesh%iNE))))then
    ierr=ierr+1
  end if
  if(any(mesh%sE/=TRI))then
    ierr=ierr+10
  end if
  if(any(mesh%nNE/=TRI_N))then
    ierr=ierr+100
  end if
  call mesh%clear()
end function

max_stars_repo_name: mianzhi/fosolvers

max_issues_repo_name: mianzhi/fosolvers

max_forks_repo_name: mianzhi/fosolvers


text:
!---------------------------------------------------------------------!
!                                PHAML                                !
!                                                                     !
! The Parallel Hierarchical Adaptive MultiLevel code for solving      !
! linear elliptic partial differential equations of the form          !
! (PUx)x + (QUy)y + RU = F on 2D polygonal domains with mixed         !
! boundary conditions, and eigenvalue problems where F is lambda*U.   !
!                                                                     !
! PHAML is public domain software.  It was produced as part of work   !
! done by the U.S. Government, and is not subject to copyright in     !
! the United States.                                                  !
!                                                                     !
!     William F. Mitchell                                             !
!     Applied and Computational Mathematics Division                  !
!     National Institute of Standards and Technology                  !
!     william.mitchell@nist.gov                                       !
!     http://math.nist.gov/phaml                                      !
!                                                                     !
!---------------------------------------------------------------------!

!----------------------------------------------------
! This file contains the main program supplied by the user.
! This is a simple example that just solves one linear elliptic pde.
!----------------------------------------------------

!       ------------
program phaml_master
!       ------------

!----------------------------------------------------
! This is the main program.
!----------------------------------------------------

!----------------------------------------------------
! Modules used are:

use phaml
!----------------------------------------------------

implicit none

!----------------------------------------------------
! Local variables

type(phaml_solution_type) :: soln
character(len=32) :: filename
!----------------------------------------------------
! Begin executable code

! set mesh file name or read it from a prompt

filename = "superior.1"
print *,"Enter the root of the mesh file names (include .1 if they contain it)"
print *,"for example, superior.1"
read *, filename

call phaml_create(soln,nproc=1, &
                  triangle_files=filename, &
                  draw_grid_who=MASTER  )

call phaml_solve_pde(soln,                   &
                     max_vert=10000,          &
                     draw_grid_when = PHASES, &
                     pause_after_phases = .true., &
                     pause_at_start = .true., &
                     print_error_who = MASTER, &
                     print_error_when = PHASES, &
                     print_error_what=ENERGY_LINF_ERR, &
                     print_errest_what=ENERGY_LINF_ERREST, &
                     print_grid_when=PHASES    , &
                     print_grid_who=MASTER)

call phaml_destroy(soln)

end program phaml_master

max_stars_repo_name: qsnake/phaml

max_issues_repo_name: qsnake/phaml

max_forks_repo_name: qsnake/phaml


text:
      program prog
      c=a+b
      end


! Using functions
! Bug: the function return type replaced by "OVERLOADED"

program main

    implicit none
    integer :: a,b,c,add1,add2
    a = 10
    b = 20

    c = add1(a,b)
    c = add2(a,b)

end program main

function add1(a, b)
    integer a, b, add1
    add1 = a + b
    return
end function add1

integer function add2(a, b)
    integer a, b
    add2 = a + b
    return
end function add2


text:
-v viterbidec.v
+ncaccess+rwc
viterbidectb.v

max_stars_repo_name: eda-ricercatore/PicenoDecoders

max_issues_repo_name: eda-globetrotter/PicenoDecoders

max_forks_repo_name: eda-globetrotter/PicenoDecoders


text:
       function sub_11(T1, T2)
        implicit none

        type(Tensor1), intent(in) :: T1, T2
        type(Tensor1) :: sub_11

        sub_11%a = T1%a - T2%a

       end function sub_11

       function sub_22(T1, T2)
        implicit none

        type(Tensor2), intent(in) :: T1, T2
        type(Tensor2) :: sub_22

        sub_22%ab = T1%ab - T2%ab

       end function sub_22

       function sub_2s2s(T1, T2)
        implicit none

        type(Tensor2s), intent(in) :: T1, T2
        type(Tensor2s) :: sub_2s2s

        sub_2s2s%a6 = T1%a6 - T2%a6

       end function sub_2s2s

       function sub_22s(T1, T2)
        implicit none

        type(Tensor2),  intent(in) :: T1
        type(Tensor2s), intent(in) :: T2
        type(Tensor2) :: sub_22s

        sub_22s = T1 - astensor(T2)

       end function sub_22s

       function sub_2s2(T1, T2)
        implicit none

        type(Tensor2s),  intent(in) :: T1
        type(Tensor2), intent(in) :: T2
        type(Tensor2) :: sub_2s2

        sub_2s2 = astensor(T1) - T2

       end function sub_2s2

       function sub_44(T1, T2)
        implicit none

        type(Tensor4), intent(in) :: T1, T2
        type(Tensor4) :: sub_44

        sub_44%abcd = T1%abcd - T2%abcd

       end function sub_44

       function sub_4s4s(T1, T2)
        implicit none

        type(Tensor4s), intent(in) :: T1, T2
        type(Tensor4s) :: sub_4s4s

        sub_4s4s%a6b6 = T1%a6b6 - T2%a6b6

       end function sub_4s4s

       function sub_44s(T1, T2)
        implicit none

        type(Tensor4),  intent(in) :: T1
        type(Tensor4s), intent(in) :: T2
        type(Tensor4) :: sub_44s

        sub_44s = T1 - astensor(T2)

       end function sub_44s

       function sub_4s4(T1, T2)
        implicit none

        type(Tensor4s),  intent(in) :: T1
        type(Tensor4), intent(in) :: T2
        type(Tensor4) :: sub_4s4

        sub_4s4 = astensor(T1) - T2

       end function sub_4s4

max_stars_repo_name: kengwit/ttb

max_issues_repo_name: kengwit/ttb

max_forks_repo_name: kengwit/ttb


text:
        SUBROUTINE YPFLAG(GRNDMD,PFLAGS,ISTAT)
*
*  Module number:
*
*  Module name: YPFLAG
*
*  Keyphrase:
*  ----------
*       Get FOS processing flags
*
*  Description:
*  ------------
*       This routine reads the FOS processing flags for the input .d0h
*       header and places them into a vector as boolean values
*
*  FORTRAN name: ypflag.for
*
*  Keywords of accessed files and tables:
*  --------------------------------------
*  Name                         I/O     Description / Comments
*       <rootname>.d0h          I       FOS science data file
*
*  Subroutines Called:
*  -------------------
*  CDBS:
*       ymsput
*  SDAS:
*       uhdgst
*
*  History:
*  --------
*  Version      Date        Author          Description
*       1       Jul 89  D. Lindler      Designed and coded
*	1.1	May 90	S. Hulbert	Added GIMP correction
*	1.2	Aug 90	S. Hulbert	Added scaling of reference background
*	1.3	Apr 93	H. Bushouse	Changed declaration of PFLAG(1) to (*)
*	1.4	Mar 94	H. Bushouse	Added SCT_CORR to PFLAGS and NAMES
*	2	Mar 94	H. Bushouse	Mods to handle PEDIGREE keywords
*					(change PFLAGS to CHAR datatype).
*	3	Oct 94	H. Bushouse	Mods for new aper, focus, flux, and
*					time corrections.
*       4       Feb98   M. De La Pena   PFLAGS(10) - PFLAGS(12) cannot be set
*                                       to PERFORM for spectropolarimetry data.
*       5       Dec01   M. De La Pena   Additional checks to ensure if WAV_CORR
*                                       is OMIT, then no component of the flux
*                                       calibration will be done.  Changes done
*                                       as defined by the STPOA group.  PFLAGS
*                                       and output keywords set to SKIPPED,
*                                       as appropriate.
*-------------------------------------------------------------------------------
*
* Inputs:
*       grndmd - ground mode
*
* Outputs:
*       pflags - boolean vector of flags
*            position           step
*               1       conversion to count rates
*		2	GIMP correction
*               3       paired pulse
*               4       background subtraction
*               5       scattered light subtraction
*               6       flat fielding
*               7       sky subtraction
*               8       wavelengths
*               9       conversion to flux (old method)
*		10	aperture throughput and focus corrections
*		11	conversion to flux (new method)
*		12	sensitivity degradation (time) correction
*               13      error propagation
*               14      mode dependent reductions
*               15      scale reference background
*
*       istat - error status
*
        CHARACTER*18 GRNDMD
        CHARACTER*8 PFLAGS(*)
        INTEGER ISTAT
C
C     UMSPUT DESTINATIONS -- CB, DAO, 4-SEP-87
C
      INTEGER STDOUT
      PARAMETER (STDOUT = 1)
      INTEGER STDERR
      PARAMETER (STDERR = 2)
C
C Common block containing input/output file descriptions
C
C       IDS - file id numbers
C       GCOUNT - group counts
C       NAXIS - naxis
C       NAXIS1 - first dimensions
C       NAXIS2 - second dimensions
C       FILL - Fill values
C
        INTEGER IDS(30),NAXIS(30),NAXIS1(30),NAXIS2(30),GCOUNT(30)
        REAL FILL(30)
        COMMON /IOCOM/IDS,GCOUNT,NAXIS,NAXIS1,NAXIS2,FILL
C
C Local variables
C
        INTEGER I
        CHARACTER*7 VAL
        CHARACTER*80 CONTXT
        CHARACTER*8 NAMES(15)
        DATA NAMES/'CNT_CORR','OFF_CORR','PPC_CORR','BAC_CORR',
     *		   'SCT_CORR','FLT_CORR','SKY_CORR','WAV_CORR',
     *		   'FLX_CORR','APR_CORR','AIS_CORR','TIM_CORR',
     *		   'ERR_CORR','MOD_CORR','GMF_CORR'/
C
C--------------------------------------------------------------------------
        DO 100 I=1,15
           PFLAGS(I)='OMIT'
           CALL UHDGST(IDS(1),NAMES(I),VAL,ISTAT)
C
           IF(ISTAT.NE.0)THEN
		IF(I.EQ.5)THEN
		   CONTXT='WARNING: SCT_CORR keyword not '//
     *			  'found in d0h file;'
		   CALL YMSPUT(CONTXT,STDOUT+STDERR,0,ISTAT)
                   CONTXT='         Scattered light correction will '//
     *                    'be skipped'
		   CALL YMSPUT(CONTXT,STDOUT+STDERR,0,ISTAT)
C These PFLAGS should be set to OMIT as the keywords do not exist in the
C input/output files, and therefore, they cannot be updated.
		   PFLAGS(5)='OMIT'
		ELSE IF(I.GE.10 .AND. I.LE.12)THEN
		   PFLAGS(I)='OMIT'
		ELSE
		   CONTXT='ERROR: getting d0h processing flag '//
     *			   NAMES(I)
		   GO TO 999
		ENDIF
	   ELSE
		IF(VAL.EQ.'PERFORM')PFLAGS(I)='PERFORM'
           ENDIF
100     CONTINUE
C
C Improved consistency check for the background-related switches
C
	IF (PFLAGS(4).EQ.'OMIT'.AND.PFLAGS(15).NE.'OMIT') PFLAGS(15)='SKIPPED'
C
C	check that the various flux calibration switches are consistent:
C	If any of the new (10, 11, 12) switches is on, turn off the
C	old abscal (switch 9).
C       If all three (10, 11, 12) are not all on, issue warning that
C       flux calibration may be incomplete.
C
C       Modified the IF blocks.  If any of PFLAGS(10) - PFLAGS(12) are PERFORM,
C       but the data is spectropolarimetry, these flags are set to SKIPPED, and
C       a warning is issued.
C
	IF (PFLAGS(10).EQ.'PERFORM' .OR. PFLAGS(11).EQ.'PERFORM' .OR.
     *       PFLAGS(12).EQ.'PERFORM') THEN
           IF ((PFLAGS(9).EQ.'PERFORM') .AND.
     *         (GRNDMD.NE.'SPECTROPOLARIMETRY'))PFLAGS(9)='SKIPPED'
           IF ((PFLAGS(9).EQ.'PERFORM') .AND.
     *         (GRNDMD.EQ.'SPECTROPOLARIMETRY')) THEN
              WRITE(CONTXT,110)
 110          FORMAT ('WARNING: APR_CORR, AIS_CORR, and',
     *             ' TIM_CORR are not applicable to')
              CALL YMSPUT(CONTXT,STDOUT+STDERR,0,ISTAT)
              WRITE(CONTXT,113)
 113          FORMAT (' spectropolarimetry data.')
              CALL YMSPUT(CONTXT,STDOUT+STDERR,0,ISTAT)
              PFLAGS(10) = 'SKIPPED'
              PFLAGS(11) = 'SKIPPED'
              PFLAGS(12) = 'SKIPPED'
           ENDIF
        ENDIF
C
	IF (PFLAGS(10).EQ.'PERFORM' .OR. PFLAGS(11).EQ.'PERFORM' .OR.
     *       PFLAGS(12).EQ.'PERFORM') THEN
           IF (PFLAGS(10).NE.'PERFORM' .OR.
     *          PFLAGS(11).NE.'PERFORM' .OR.
     *          PFLAGS(12).NE.'PERFORM')THEN
              WRITE(CONTXT,111)
 111          FORMAT ('WARNING: APR_CORR, AIS_CORR, and',
     *             ' TIM_CORR are not all on together.')
              CALL YMSPUT(CONTXT,STDOUT+STDERR,0,ISTAT)
              WRITE(CONTXT,112)
 112          FORMAT ('         Flux calibration may not',
     *             ' be fully applied.')
              CALL YMSPUT(CONTXT,STDOUT+STDERR,0,ISTAT)
           END IF
        END IF
C
C Begin additions defined by the ST Post-Operational Archive group at ST-ECF.
C MDD NOTE: The check for conversion to absolute flux units (PFLAGS(11)) is
C done in YCLPRC; this check has been instituted since December 1995.
C MDD: Set internal PFLAGS to SKIPPED (not OMIT) so output keywords are updated.
C
C set the aperture throughput to SKIPPED if the WAV_CORR=OMIT
        IF(PFLAGS(8).NE.'PERFORM'.AND.PFLAGS(10).EQ.'PERFORM')THEN
           PFLAGS(10)='SKIPPED'
           WRITE(CONTXT,116)
 116       FORMAT ('WARNING: APR_CORR requires that WAV_CORR',
     *             ' be executed.')
           CALL YMSPUT(CONTXT,STDERR+STDOUT,0,ISTAT)
           WRITE(CONTXT,117)
 117       FORMAT ('         APR_CORR will not be performed.')
           CALL YMSPUT(CONTXT,STDERR+STDOUT,0,ISTAT)
        ENDIF

C set the sensitivity degradation (time) correction to SKIPPED if
C the WAV_CORR=OMIT
        IF(PFLAGS(8).NE.'PERFORM'.AND.PFLAGS(12).EQ.'PERFORM')THEN
           PFLAGS(12)='SKIPPED'
           WRITE(CONTXT,118)
 118       FORMAT ('WARNING: TIM_CORR requires that WAV_CORR',
     *             ' be executed.')
           CALL YMSPUT(CONTXT,STDERR+STDOUT,0,ISTAT)
           WRITE(CONTXT,119)
 119       FORMAT ('         TIM_CORR will not be performed.')
           CALL YMSPUT(CONTXT,STDERR+STDOUT,0,ISTAT)
        ENDIF
C
C End additions defined by the ST Post-Operational Archive group at ST-ECF.
C
        ISTAT=0
        GO TO 1000
999     CALL YMSPUT(CONTXT,STDOUT+STDERR,0,ISTAT)
        ISTAT=1
1000    RETURN
        END

max_stars_repo_name: iraf-community/stsdas


text:
      MODULE dopen_I
      INTERFACE
!...Generated by Pacific-Sierra Research 77to90  4.4G  10:47:01  03/09/06
      SUBROUTINE dopen (C, MDIM, NORBS, NDUBL, NSINGL, FRACT, SDM)
      USE vast_kind_param,ONLY: DOUBLE
      INTEGER, INTENT(IN) :: MDIM
      INTEGER, INTENT(IN) :: NORBS
      INTEGER, INTENT(IN) :: NDUBL
      INTEGER, INTENT(IN) :: NSINGL
      REAL(DOUBLE), INTENT(IN) :: FRACT
      REAL(DOUBLE), DIMENSION(MDIM,MDIM), INTENT(IN) :: C
      REAL(DOUBLE), DIMENSION((MDIM*(MDIM + 1))/2), INTENT(OUT) :: SDM
      END SUBROUTINE
      END INTERFACE
      END MODULE

max_stars_repo_name: openmopac/MOPAC-archive

max_issues_repo_name: openmopac/MOPAC-archive

max_forks_repo_name: openmopac/MOPAC-archive


text:
        subroutine e3wmlt (shp,    shg,    WdetJ,
     &                     ri,     rmi,    rl,
     &                     rml,    stiff,  EGmass )
c
c----------------------------------------------------------------------
c
c Up to now most of the terms have not been multiplied by the
c shape function from the weight space.  Now that we have collected
c all the terms that the shape function (and its derivatives for the
c terms that were integrated by parts) multiplies, in this routine
c we carry out that multiplication. After these operations we will
c have this quadrature points contribution to the integral for the
c residual (i.e.  g^e_b(xi^l) D(xi^l) QW(xi^l)  where e is the element,
c b is the local node number on the element, l is the quadrature point
c D is the determinate of the Jacobian of the mapping, and QW is this
c quadrature points weight....WHEW).  When we add up all of the
c integration points we get G^e_b which we will assemble in the
c subroutine LOCAL to form G_B.
c
c This routine also forms the LHS contribution from the LS term
c and the diffusion term which has been partially constructed and
c placed in stiff.  Those familiar with elasticity might recognize
c this naming convention since this is like a stiffness matrix that
c if you had a linear problem would be calculated once and saved for
c all time.
c
c    ri:  LS, Diffusion, Convective and mass,
c   rmi:  LS, Diffusion and mass,
c stiff:  LS, Diffusion.
c
c input:
c  shp    (nshl)                 : element shape-functions
c  shg    (npro,nshl,nsd)        : element grad-shape-functions
c  WdetJ  (npro)                   : weighted Jacobian
c  ri     (npro,nflow*(nsd+1))      : partial residual
c  rmi    (npro,nflow*(nsd+1))      : partial modified residual
c  stiff  (npro,nsd*nflow,nsd*nflow) : stiffness matrix
c                                    ( K_ij + A_i tau A_j )
c
c output:
c  rl     (npro,nshl,nflow)       : residual
c  rml    (npro,nshl,nflow)       : modified residual
c  EGmass (npro,nedof,nedof)       : element LHS tangent mass matrix
c
c
c Zdenek Johan, Summer 1990.  (Modified from e2assm.f)
c Zdenek Johan, Winter 1991.  (Fortran 90)
c Kenneth Jansen, Winter 1997  Prim. variables
c----------------------------------------------------------------------
c
        include "common.h"
c
        dimension shp(npro,nshl),
     &            shg(npro,nshl,nsd),
     &            WdetJ(npro),             ri(npro,nflow*(nsd+1)),
     &            rmi(npro,nflow*(nsd+1)), stiff(npro,3*nflow,3*nflow),
     &            rl(npro,nshl,nflow),     rml(npro,nshl,nflow),
     &            EGmass(npro,nedof,nedof)
c
        dimension shg1(npro),              shg2(npro),
     &            shg3(npro),              stif1(npro,nflow,nflow),
     &            stif2(npro,nflow,nflow), stif3(npro,nflow,nflow)

        ttim(29) = ttim(29) - secs(0.0)
c
c.... ---------------------------->  RHS  <----------------------------
c
c.... add spatial contribution to rl and rml
c
c.... ires = 1 or 3
c
        if ((ires .eq. 1) .or. (ires .eq. 3)) then
c
          do i = 1, nshl
            rl(:,i,1) = rl(:,i,1) + WdetJ * (
     &                  shg(:,i,1) * ri(:, 1) + shg(:,i,2) * ri(:, 6)
     &                                        + shg(:,i,3) * ri(:,11) )
            rl(:,i,2) = rl(:,i,2) + WdetJ * (
     &                  shg(:,i,1) * ri(:, 2) + shg(:,i,2) * ri(:, 7)
     &                                        + shg(:,i,3) * ri(:,12) )
            rl(:,i,3) = rl(:,i,3) + WdetJ * (
     &                  shg(:,i,1) * ri(:, 3) + shg(:,i,2) * ri(:, 8)
     &                                        + shg(:,i,3) * ri(:,13) )
            rl(:,i,4) = rl(:,i,4) + WdetJ * (
     &                  shg(:,i,1) * ri(:, 4) + shg(:,i,2) * ri(:, 9)
     &                                        + shg(:,i,3) * ri(:,14) )
            rl(:,i,5) = rl(:,i,5) + WdetJ * (
     &                  shg(:,i,1) * ri(:, 5) + shg(:,i,2) * ri(:,10)
     &                                        + shg(:,i,3) * ri(:,15) )
          enddo
c
          flops = flops + 36*nshl*npro
        endif
c
c.... ires = 2 or 3
c
        if ((ires .eq. 2) .or. (ires .eq. 3)) then
          do i = 1, nshl
            rml(:,i,1) = rml(:,i,1) + WdetJ * (
     &                   shg(:,i,1) * rmi(:, 1) + shg(:,i,2) * rmi(:, 6)
     &                 + shg(:,i,3) * rmi(:,11)
     &                   + shp(:,i) * rmi(:,16)  )
            rml(:,i,2) = rml(:,i,2) + WdetJ * (
     &                   shg(:,i,1) * rmi(:, 2) + shg(:,i,2) * rmi(:, 7)
     &                 + shg(:,i,3) * rmi(:,12)
     &                   + shp(:,i) * rmi(:,17)    )
            rml(:,i,3) = rml(:,i,3) + WdetJ * (
     &                   shg(:,i,1) * rmi(:, 3) + shg(:,i,2) * rmi(:, 8)
     &                 + shg(:,i,3) * rmi(:,13)
     &                   + shp(:,i) * rmi(:,18)    )
            rml(:,i,4) = rml(:,i,4) + WdetJ * (
     &                   shg(:,i,1) * rmi(:, 4) + shg(:,i,2) * rmi(:, 9)
     &                 + shg(:,i,3) * rmi(:,14)
     &                   + shp(:,i) * rmi(:,19)    )
            rml(:,i,5) = rml(:,i,5) + WdetJ * (
     &                   shg(:,i,1) * rmi(:, 5) + shg(:,i,2) * rmi(:,10)
     &                 + shg(:,i,3) * rmi(:,15)
     &                   + shp(:,i) * rmi(:,20)    )
          enddo
c
          flops = flops + (15+45*nshl)*npro
        endif
c
c.... add temporal contribution to rl
c
        if (ngauss .eq. 1 .and. nshl.eq.4) then !already ex. integ mass
         if(matflg(5,1).ge.1) then
          do i = 1, nshl
            rl(:,i,2) = rl(:,i,2) + shp(:,i) * WdetJ * ri(:,17)
            rl(:,i,3) = rl(:,i,3) + shp(:,i) * WdetJ * ri(:,18)
            rl(:,i,4) = rl(:,i,4) + shp(:,i) * WdetJ * ri(:,19)
            rl(:,i,5) = rl(:,i,5) + shp(:,i) * WdetJ * ri(:,20)
          enddo
         endif
       else   !check for a body force
          do i = 1, nshl
            rl(:,i,1) = rl(:,i,1) + shp(:,i) * WdetJ * ri(:,16)
            rl(:,i,2) = rl(:,i,2) + shp(:,i) * WdetJ * ri(:,17)
            rl(:,i,3) = rl(:,i,3) + shp(:,i) * WdetJ * ri(:,18)
            rl(:,i,4) = rl(:,i,4) + shp(:,i) * WdetJ * ri(:,19)
            rl(:,i,5) = rl(:,i,5) + shp(:,i) * WdetJ * ri(:,20)
          enddo

          flops = flops + 11*nshl*npro
       endif

c
c.... ---------------------------->  LHS  <----------------------------
c
       if (lhs .eq. 1) then
c
c.... loop through columns (nodes j)
c
          do j = 1, nshl
             j0 = nflow * (j - 1)
c
c.... set up factors
c
             shg1 = WdetJ * shg(:,j,1)
             shg2 = WdetJ * shg(:,j,2)
             shg3 = WdetJ * shg(:,j,3)
c
c.... loop through d.o.f.'s
c
             do jdof = 1, nflow
                do idof = 1, nflow
                   idof2 = idof  + nflow
                   jdof2 = jdof  + nflow

                   idof3 = idof2 + nflow
                   jdof3 = jdof2 + nflow
c
c.... calculate weighted stiffness matrix (first part)
c
                   stif1(:,idof,jdof) = shg1 * stiff(:,idof,jdof)
     &                                + shg2 * stiff(:,idof,jdof2)
     &                                + shg3 * stiff(:,idof,jdof3)
                   stif2(:,idof,jdof) = shg1 * stiff(:,idof2,jdof)
     &                                + shg2 * stiff(:,idof2,jdof2)
     &                                + shg3 * stiff(:,idof2,jdof3)
                   stif3(:,idof,jdof) = shg1 * stiff(:,idof3,jdof)
     &                                + shg2 * stiff(:,idof3,jdof2)
     &                                + shg3 * stiff(:,idof3,jdof3)
                enddo
             enddo
c
c.... loop through rows (nodes i)
c
             do i = 1, nshl
                i0 = nflow * (i - 1)
c
c.... add contribution of stiffness to EGmass
c
                do jdof = 1, nflow
                   EGmass(:,i0+1,j0+jdof) = EGmass(:,i0+1,j0+jdof)
     &                  + shg(:,i,1) * stif1(:,1,jdof)
     &                  + shg(:,i,2) * stif2(:,1,jdof)
     &                  + shg(:,i,3) * stif3(:,1,jdof)
                   EGmass(:,i0+2,j0+jdof) = EGmass(:,i0+2,j0+jdof)
     &                  + shg(:,i,1) * stif1(:,2,jdof)
     &                  + shg(:,i,2) * stif2(:,2,jdof)
     &                  + shg(:,i,3) * stif3(:,2,jdof)
                   EGmass(:,i0+3,j0+jdof) = EGmass(:,i0+3,j0+jdof)
     &                  + shg(:,i,1) * stif1(:,3,jdof)
     &                  + shg(:,i,2) * stif2(:,3,jdof)
     &                  + shg(:,i,3) * stif3(:,3,jdof)
                   EGmass(:,i0+4,j0+jdof) = EGmass(:,i0+4,j0+jdof)
     &                  + shg(:,i,1) * stif1(:,4,jdof)
     &                  + shg(:,i,2) * stif2(:,4,jdof)
     &                  + shg(:,i,3) * stif3(:,4,jdof)
                   EGmass(:,i0+5,j0+jdof) = EGmass(:,i0+5,j0+jdof)
     &                  + shg(:,i,1) * stif1(:,5,jdof)
     &                  + shg(:,i,2) * stif2(:,5,jdof)
     &                  + shg(:,i,3) * stif3(:,5,jdof)
                enddo
c
c.... end loop on rows
c
             enddo
c
c.... end loop on columns
c
          enddo
c
c.... end left hand side assembly
c
       endif

       ttim(29) = ttim(29) + secs(0.0)
c
c.... return
c
        return
        end
c
c
c
        subroutine e3wmltSclr(shp,    shg,    WdetJ,
     &                        rti,    rtl,
     &                        stifft, EGmasst )
c
c----------------------------------------------------------------------
c
c Up to now most of the terms have not been multiplied by the
c shape function from the weight space.  Now that we have collected
c all the terms that the shape function (and its derivatives for the
c terms that were integrated by parts) multiplies, in this routine
c we carry out that multiplication. After these operations we will
c have this quadrature points contribution to the integral for the
c residual (i.e.  g^e_b(xi^l) D(xi^l) QW(xi^l)  where e is the element,
c b is the local node number on the element, l is the quadrature point
c D is the determinate of the Jacobian of the mapping, and QW is this
c quadrature points weight....WHEW).  When we add up all of the
c integration points we get G^e_b which we will assemble in the
c subroutine LOCAL to form G_B.
c
c This routine also forms the LHS contribution from the LS term
c and the diffusion term which has been partially constructed and
c placed in stiff.  Those familiar with elasticity might recognize
c this naming convention since this is like a stiffness matrix that
c if you had a linear problem would be calculated once and saved for
c all time.
c
c    ri:  LS, Diffusion, Convective and mass,
c stiff:  LS, Diffusion.
c
c input:
c  shp     (npro,nshl)            : element shape-functions
c  shg     (npro,nshl,nsd)        : element grad-shape-functions
c  WdetJ   (npro)                 : weighted Jacobian
c  rti     (npro,nsd+1)           : partial residual
c  stifft  (npro,nsd,nsd)         : stiffness matrix
c                                    ( K_ij + A_i tau A_j )
c
c output:
c  rtl     (npro,nshl)            : residual  (will end up being G^e_b)
c  EGmasst (npro,nshape,nshape)   : element LHS tangent mass matrix
c                                  (will end up being dG^e_b/dY_a)
c
c
c Zdenek Johan, Summer 1990.  (Modified from e2assm.f)
c Zdenek Johan, Winter 1991.  (Fortran 90)
c Kenneth Jansen, Winter 1997  Prim. variables
c----------------------------------------------------------------------
c
        include "common.h"
c
        dimension shp(npro,nshl),            shg(npro,nshl,nsd),
     &            WdetJ(npro),               rti(npro,nsd+1),
     &            rmti(npro,nsd+1),          stifft(npro,nsd,nsd),
     &            rtl(npro,nshl),            rmtl(npro,nshl),
     &            EGmasst(npro,nshape,nshape)
c
        dimension shg1(npro),                shg2(npro),
     &            shg3(npro),                stift1(npro),
     &            stift2(npro),              stift3(npro)

        ttim(29) = ttim(29) - tmr(0.0)
c
c.... ---------------------------->  RHS  <----------------------------
c
c.... add spatial contribution to rl and rml
c
c.... ires = 1 or 3
c
        if ((ires .eq. 1) .or. (ires .eq. 3)) then
c
          do i = 1, nshl
            rtl(:,i) = rtl(:,i) + WdetJ * (
     &                  shg(:,i,1) * rti(:,1) + shg(:,i,2) * rti(:,2)
     &                                        + shg(:,i,3) * rti(:,3) )

          enddo
          flops = flops + 36*nshl*npro
        endif
c
c.... ires = 2 or 3
c
c        if ((ires .eq. 2) .or. (ires .eq. 3)) then
c          do i = 1, nshl
c            rml(:,i,1) = rml(:,i,1) + WdetJ * (
c     &                   shg(:,i,1) * rmi(:, 1) + shg(:,i,2) * rmi(:, 6)
c     &                 + shg(:,i,3) * rmi(:,11) + shp(:,i) * rmi(:,16) )
c            rml(:,i,2) = rml(:,i,2) + WdetJ * (
c     &                   shg(:,i,1) * rmi(:, 2) + shg(:,i,2) * rmi(:, 7)
c     &                 + shg(:,i,3) * rmi(:,12) + shp(:,i) * rmi(:,17) )
c            rml(:,i,3) = rml(:,i,3) + WdetJ * (
c     &                   shg(:,i,1) * rmi(:, 3) + shg(:,i,2) * rmi(:, 8)
c     &                 + shg(:,i,3) * rmi(:,13) + shp(:,i) * rmi(:,18) )
c            rml(:,i,4) = rml(:,i,4) + WdetJ * (
c     &                   shg(:,i,1) * rmi(:, 4) + shg(:,i,2) * rmi(:, 9)
c     &                 + shg(:,i,3) * rmi(:,14) + shp(:,i) * rmi(:,19) )
c            rml(:,i,5) = rml(:,i,5) + WdetJ * (
c     &                   shg(:,i,1) * rmi(:, 5) + shg(:,i,2) * rmi(:,10)
c     &                 + shg(:,i,3) * rmi(:,15) + shp(:,i) * rmi(:,20) )
c          enddo
c
c          flops = flops + (15+45*nshl)*npro
c        endif

c
c.... add temporal contribution to rl
c

          do i = 1, nshl
            rtl(:,i) = rtl(:,i) + shp(:,i) * WdetJ * rti(:,4)
          enddo

          flops = flops + 11*nshl*npro

c
c.... ---------------------------->  LHS  <----------------------------
c
       if (lhs .eq. 1) then
c
c.... loop through columns (nodes j)
c
          do j = 1, nshl
c
c.... set up factors
c
             shg1 = WdetJ * shg(:,j,1)
             shg2 = WdetJ * shg(:,j,2)
             shg3 = WdetJ * shg(:,j,3)
c
c.... calculate weighted stiffness matrix (first part)
c
                   stift1(:) = shg1 * stifft(:,1,1)
     &                                + shg2 * stifft(:,1,2)
     &                                + shg3 * stifft(:,1,3)
                   stift2(:) = shg1 * stifft(:,2,1)
     &                                + shg2 * stifft(:,2,2)
     &                                + shg3 * stifft(:,2,3)
                   stift3(:) = shg1 * stifft(:,3,1)
     &                                + shg2 * stifft(:,3,2)
     &                                + shg3 * stifft(:,3,3)
c
c.... loop through rows (nodes i)
c
             do i = 1, nshl
c
c.... add contribution of stiffness to EGmass
c
                   EGmasst(:,i,j) = EGmasst(:,i,j)
     &                  + shg(:,i,1) * stift1(:)
     &                  + shg(:,i,2) * stift2(:)
     &                  + shg(:,i,3) * stift3(:)

c
c.... end loop on rows
c
             enddo
c
c.... end loop on columns
c
          enddo
c
c.... end left hand side assembly
c
       endif

       ttim(29) = ttim(29) + tmr(0.0)
c
c.... return
c
        return
        end


text:
chapter xvii, example 3. bivariate spline interpolation
c  followed by conversion to pp-repr and evaluation
c  from  * a practical guide to splines *  by c. de boor
calls  spli2d(bsplvb,banfac/slv),interv,bspp2d(bsplvb*),ppvalu(interv*)
c     parameter nx=7,kx=3,ny=6,ky=4
c     integer i,iflag,j,jj,lefty,lx,ly,mflag
c     real bcoef(nx,ny),breakx(6),breaky(4),coef(kx,5,ky,3),taux(nx)
c    *	  ,tauy(ny),tx(10),ty(10)
c    *	       ,work1(nx,ny),work2(nx),work3(42)
c    *	       ,work4(kx,kx,ny),work5(ny,kx,5),work6(ky,ky,21)
      integer i,iflag,j,jj,kp1,kx,ky,lefty,lx,ly,mflag,nx,ny
      real bcoef(7,6),breakx(6),breaky(4),coef(3,5,4,3),taux(7)
     *	  ,tauy(6),tx(10),ty(10)
     *	       ,work1(7,6),work2(7),work3(42)
     *	       ,work4(3,3,6),work5(6,3,5),work6(4,4,21)
      real ppvalu,g,x,y
      data nx,kx,ny,ky / 7,3,6,4 /
c     note that , with the above parameters, lx=5, ly = 3
      equivalence (work4,work6)
      g(x,y) = amax1(x-3.5,0.)**2 + amax1(y-3.,0.)**3
c *** set up data points and knots
c     in x, interpolate between knots by parabolic splines, using
c     not-a-knot end condition
      do 10 i=1,nx
   10	 taux(i) = float(i)
      do 11 i=1,kx
	 tx(i) = taux(1)
   11	 tx(nx+i) = taux(nx)
      kp1 = kx+1
      do 12 i=kp1,nx
   12	 tx(i) = (taux(i-kx+1) + taux(i-kx+2))/2.
c     in y, interpolate at knots by cubic splines, using not-a-knot
c     end condition
      do 20 i=1,ny
   20	 tauy(i) = float(i)
      do 21 i=1,ky
	 ty(i) = tauy(1)
   21	 ty(ny+i) = tauy(ny)
      kp1 = ky+1
      do 22 i=kp1,ny
   22	 ty(i) = tauy(i-ky+2)
c  *** generate and print out function values
      print 620,(tauy(i),i=1,ny)
  620 format(11h given data//6f12.1)
      do 32 i=1,nx
	 do 31 j=1,ny
   31	    bcoef(i,j) = g(taux(i),tauy(j))
   32	 print 632,taux(i),(bcoef(i,j),j=1,ny)
  632 format(f5.1,6e12.5)
c
c  *** construct b-coefficients of interpolant
c
      call spli2d(taux,bcoef,tx,nx,kx,ny,work2,work3,work1,iflag)
      call spli2d(tauy,work1,ty,ny,ky,nx,work2,work3,bcoef,iflag)
c
c  *** convert to pp-representation
c
      call bspp2d(tx,bcoef,nx,kx,ny,work4,breakx,work5,lx)
      call bspp2d(ty,work5,ny,ky,lx*kx,work6,breaky,coef,ly)
c
c  *** evaluate interpolation error at mesh points and print out
      print 640,(tauy(j),j=1,ny)
  640 format(//20h interpolation error//6f12.1)
      do 45 j=1,ny
	 call interv(breaky,ly,tauy(j),lefty,mflag)
	 do 45 i=1,nx
	    do 41 jj=1,ky
   41	     work2(jj)=ppvalu(breakx,coef(1,1,jj,lefty),lx,kx,taux(i),0)
   45	    work1(i,j) = g(taux(i),tauy(j)) -
     *	    ppvalu(breaky(lefty),work2,1,ky,tauy(j),0)
      do 46 i=1,nx
   46	 print 632,taux(i),(work1(i,j),j=1,ny)
					stop
      end

max_stars_repo_name: ysBach/irafdocgen

max_issues_repo_name: kirxkirx/iraf

max_forks_repo_name: kirxkirx/iraf


text:
	SUBROUTINE HQCOMP ( icolr, color, ired, igreen, iblue, xname,
     +			    iret )
C************************************************************************
C* HQCOMP - NC								*
C*									*
C* This subroutine returns the red, green, and blue components of a	*
C* color.  The color components are defined in the range 0 - 255.  If	*
C* the color was defined by name, the GEMPAK color name is also		*
C* returned.  The X Window System color name is returned if it is	*
C* available.								*
C*									*
C* HQCOMP  ( ICOLR, COLOR, IRED, IGREEN, IBLUE, XNAME, IRET )		*
C*									*
C* Input parameters:							*
C* 	ICOLR		INTEGER		Color number			*
C*									*
C* Output parameters:							*
C*	COLOR		CHAR*		GEMPAK color name		*
C*	IRED		INTEGER		Red color component		*
C*	IGREEN		INTEGER		Green color component		*
C*	IBLUE		INTEGER		Blue color component		*
C*	XNAME		CHAR*		X Window System color name	*
C*	IRET		INTEGER 	Return code			*
C**									*
C* Log:									*
C* L. Williams/EAI	 3/96						*
C* S. Jacobs/NCEP	 4/96	Copied from the XW to the NC driver	*
C* S. Maxwell/GSC        6/97   Documentation changes                   *
C* T. Lee/GSC		 9/97	Included ERROR.PRM			*
C************************************************************************
	INCLUDE		'ERROR.PRM'
C*
	CHARACTER*(*) 	color, xname
C-----------------------------------------------------------------------
	iret = NORMAL
C
C*	Retrieve the values from structure.
C
	color(1:1) = ' '
	xname(1:1) = ' '

 	CALL MQCOMP( icolr, color, ired, igreen, iblue, xname, ilen,
     +		     ilxn, iret )

C
C*	Clear the unused space of Color name and X color name
C
 	color(ilen+1: ) = ' '
 	xname (ilxn+1: ) = ' '

C*
	RETURN
	END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


text:
!==============================================================================!
  subroutine Compute_F22(grid, ini, phi)
!------------------------------------------------------------------------------!
!   Discretizes and solves eliptic relaxation equations for f22.               !
!------------------------------------------------------------------------------!
!---------------------------------[Modules]------------------------------------!
  use Flow_Mod
  use Les_Mod
  use Rans_Mod
  use Comm_Mod
  use Var_Mod
  use Grid_Mod
  use Grad_Mod
  use Info_Mod
  use Numerics_Mod
  use Solvers_Mod, only: Bicg, Cg, Cgs
  use Control_Mod
  use Work_Mod,    only: phi_x => r_cell_01,  &
                         phi_y => r_cell_02,  &
                         phi_z => r_cell_03
!------------------------------------------------------------------------------!
  implicit none
!--------------------------------[Arguments]-----------------------------------!
  type(Grid_Type) :: grid
  integer         :: ini
  type(Var_Type)  :: phi
!----------------------------------[Locals]------------------------------------!
  integer           :: s, c, c1, c2, niter
  real              :: f_ex, f_im
  real              :: a0, a12, a21
  real              :: ini_res, tol
  real              :: phi_x_f, phi_y_f, phi_z_f
  character(len=80) :: precond
  real              :: urf           ! under-relaxation factor
!==============================================================================!
!                                                                              !
!   The form of equations which are solved:                                    !
!                                                                              !
!        /             /              /                                        !
!       |   df22      |   f22        |  f22hg                                  !
!       | - ---- dS + | ------ dV  = | ------- dV                              !
!       |    dy       |  Lsc^2       |  Lsc^2                                  !
!      /             /              /                                          !
!                                                                              !
!   Dimension of the system under consideration                                !
!                                                                              !
!     [A]{f22} = {b}   [kg K/s]                                                !
!                                                                              !
!   Dimensions of certain variables:                                           !
!                                                                              !
!     f22            [1/s]                                                     !
!     Lsc            [m]                                                       !
!------------------------------------------------------------------------------!

  a % val = 0.0

  b = 0.0

  ! Old values (o) and older than old (oo)
  if(ini .eq. 1) then
    do c = 1, grid % n_cells
      phi % oo(c)   = phi % o(c)
      phi % o (c)   = phi % n(c)
    end do
  end if

  ! New values
  do c = 1, grid % n_cells
    phi % c(c) = 0.0
  end do

  ! Gradients
  call Grad_Mod_For_Phi(grid, phi % n, 1, phi_x, .true.)
  call Grad_Mod_For_Phi(grid, phi % n, 2, phi_y, .true.)
  call Grad_Mod_For_Phi(grid, phi % n, 3, phi_z, .true.)

  !------------------!
  !                  !
  !     Difusion     !
  !                  !
  !------------------!

  !----------------------------!
  !   Spatial discretization   !
  !----------------------------!
  do s = 1, grid % n_faces

    c1 = grid % faces_c(1,s)
    c2 = grid % faces_c(2,s)

    phi_x_f = fw(s) * phi_x(c1) + (1.0 - fw(s)) * phi_x(c2)
    phi_y_f = fw(s) * phi_y(c1) + (1.0 - fw(s)) * phi_y(c2)
    phi_z_f = fw(s) * phi_z(c1) + (1.0 - fw(s)) * phi_z(c2)

    ! Total (exact) diffusive flux
    f_ex = (  phi_x_f * grid % sx(s)   &
            + phi_y_f * grid % sy(s)   &
            + phi_z_f * grid % sz(s) )

    a0 = f_coef(s)

    ! Implicit diffusive flux
    f_im=(   phi_x_f * grid % dx(s)        &
           + phi_y_f * grid % dy(s)        &
           + phi_z_f * grid % dz(s)) * a0

    ! Cross diffusion part
    phi % c(c1) = phi % c(c1) + f_ex - f_im
    if(c2  > 0) then
      phi % c(c2) = phi % c(c2) - f_ex + f_im
    end if

    ! Calculate the coefficients for the sysytem matrix
    a12 = a0
    a21 = a0

    ! Fill the system matrix
    if(c2  > 0) then
      a % val(a % pos(1,s)) = a % val(a % pos(1,s)) - a12
      a % val(a % dia(c1))  = a % val(a % dia(c1))  + a12
      a % val(a % pos(2,s)) = a % val(a % pos(2,s)) - a21
      a % val(a % dia(c2))  = a % val(a % dia(c2))  + a21
    else if(c2  < 0) then

      ! Inflow
      if( (Grid_Mod_Bnd_Cond_Type(grid,c2) .eq. INFLOW)) then
        a % val(a % dia(c1)) = a % val(a % dia(c1)) + a12
        b(c1) = b(c1) + a12 * phi % n(c2)
      end if

      ! Wall and wall flux; solid walls in any case
      if( (Grid_Mod_Bnd_Cond_Type(grid,c2) .eq. WALL).or.       &
          (Grid_Mod_Bnd_Cond_Type(grid,c2) .eq. WALLFL) ) then
        a % val(a % dia(c1)) = a % val(a % dia(c1)) + a12
        !---------------------------------------------------------------!
        !   Source coefficient is filled in SourceF22.f90 in order to   !
        !   get updated values of f22 on the wall.  Otherwise f22       !
        !   equation does not converge very well                        !
        !   b(c1) = b(c1) + a12 * phi % n(c2)                           !
        !---------------------------------------------------------------!
      end if
    end if

  end do  ! through faces

  ! Cross diffusion terms are treated explicity
  do c = 1, grid % n_cells
    b(c) = b(c) + phi % c(c)
  end do

  !-------------------------------------!
  !                                     !
  !   Source terms and wall function    !
  !   (Check if it is good to call it   !
  !    before the under relaxation ?)   !
  !                                     !
  !-------------------------------------!
  if(turbulence_model .eq. RSM_MANCEAU_HANJALIC) then
    call Source_F22_Rsm_Manceau_Hanjalic(grid)
  else
    call Source_F22_K_Eps_Zeta_F(grid)
  end if

  !---------------------------------!
  !                                 !
  !   Solve the equations for phi   !
  !                                 !
  !---------------------------------!

  ! Set under-relaxation factor then overwrite with control file if specified
  urf = 1.0
  call Control_Mod_Simple_Underrelaxation_For_Turbulence(urf)

  do c = 1, grid % n_cells
    b(c) = b(c) + a % val(a % dia(c)) * (1.0 - urf)*phi % n(c) / urf
    a % val(a % dia(c)) = a % val(a % dia(c)) / urf
  end do

  ! Get tolerance for linear solver
  call Control_Mod_Tolerance_For_Turbulence_Solver(tol)

  ! Get matrix precondioner
  call Control_Mod_Preconditioner_For_System_Matrix(precond)

  ! Set the number of iterations then overwrite with control file if specified
  niter = 6
  call Control_Mod_Max_Iterations_For_Turbulence_Solver(niter)

  call Cg(a, phi % n, b, precond, niter, tol, ini_res, phi % res)

  call Info_Mod_Iter_Fill_At(3, 4, phi % name, niter, phi % res)

  call Comm_Mod_Exchange_Real(grid, phi % n)

  end subroutine

max_stars_repo_name: MassimoBorrelliPhysicist/T-Flows

max_issues_repo_name: MassimoBorrelliPhysicist/T-Flows

max_forks_repo_name: MassimoBorrelliPhysicist/T-Flows


text:
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C   WAIT - Wait for a Given Number of Seconds Subroutine
C
C   VERSION:  1.0
C   AUTHOR:  A. Roberts, Tesseract Systems, Inc.
C   CREATION DATE:  05/15/83
C   LAST REVISION DATE:  05/20/83
C   LAST REVISED BY:  A. Roberts
C
C FUNCTION:  WAIT is called to delay for the specified number of
C		seconds, then return to the caller.  It makes use
C		of VAX/VMS timer system services to do this.
C
C USEAGE:  CALL WAIT(SECS)
C
C ARGUMENTS:	SECS - Input Number of Seconds to Wait	(REAL*4)
C
C ERRORS:  None known
C
C NOTES:  Written in VAX/VMS FORTRAN.
C
C PROCESS:  The only reason for the existance of WAIT is the
C		relatively short length of time available with
C		a single system service call, as well as the
C		reasonably difficult time coding.  WAIT will
C		break a long time request into shorter requests
C		if need be, and call timer services for each
C		short request.
C
C REVISION HISTORY:  0.0 - Created on 05/15/83 by A. Roberts,
C			Tesseract Systems, Inc.
C					------
C		     1.0 - Initial release version.  05/20/83,
C			A. Roberts
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

	subroutine WAIT(SECS)

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C Variable Declaration, Documentation, Constants, etc.
C
C Outside of the argument, there is little to WAIT, except a working
C copy of the time and the argument to the system service.
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

	implicit none

	parameter BIG = 2147482752.0		! Max time on a call
	parameter WAIT_EFN = 50			! Event flag used

	real*4 secs				! Input Seconds

	real*4 timecop				! Working Variables
	integer*4 itime(2)

	external sys$setimr			! System services used
	external sys$waitfr

C
C Initialize the ITIME array to be the largest possible
C call for the system timer routine.
C
C Also, convert the SECS argument to the required units for
C the system call.
C

	itime(2) = -1
	itime(1) = int(-1.0 * BIG)

	timecop = secs * 1.0E+07

C
C Now decide what to do.  If the required time is small enough,
C make one call to the system service and then return.  Otherwise,
C loop, making the largest possible calls and removing that
C value from the user specified time until the time has elasped.
C

C
C Time is small enough
C
	if(timecop .lt. BIG) then

		itime(1) = int(-1.0 * timecop)

		call sys$setimr(%val(WAIT_EFN),itime,,)
		call sys$waitfr(%val(WAIT_EFN))

C
C Time is not small enough
C
	else

		do while(timecop .ge. BIG)

			call sys$setimr(%val(WAIT_EFN),itime,,)
			call sys$waitfr(%val(WAIT_EFN))

			timecop = timecop - BIG

		end do

			itime(1) = int(-1.0 * timecop)

			call sys$setimr(%val(WAIT_EFN),itime,,)
			call sys$waitfr(%val(WAIT_EFN))

	end if

C
C Thats all.  Return to the caller.
C

	return
	end

ext: for