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

!===============================================================================
! Read in rank-specific hotstart outputs and combine them into hotstart.in.
! Gobal-local mappings are read in from separate files.

! Usage: ./combine_hotstart6 -h for help, inside outputs/
! Inputs:
!        rank-specific hotstart files;
!        local_to_global_*;
!        screen: ntracers; it_char (iteration #)
! Output: hotstart_it=[time step].nc
!
!  ifort -O2 -cpp -CB -mcmodel=medium -assume byterecl -g -traceback -o combine_hotstart7.exe ../UtilLib/argparse.f90 combine_hotstart7.f90 -I$NETCDF/include -I$NETCDF_FORTRAN/include -L$NETCDF_FORTRAN/lib -L$NETCDF/lib -lnetcdf -lnetcdff

! Revisions: v6 with nc
!================================================================================
program combine_hotstart
integer :: istep
integer :: nproc
integer :: ntracer
call combine_hotstart7_input(istep)
call combine_hotstart7(istep)
end program combine_hotstart

subroutine combine_hotstart7_input(istep)
use argparse
implicit none
integer       :: istep   !< global iteration/step
!integer       :: nproc   !< number of processors in original simulation
!integer       :: ntracer !< number of tracers

cmd_name = "combine_hotstart7"
call cla_init(cmd_name)
call cla_register('-i','--iteration', 'global iteration to combine (before _00*_hotstart', cla_int,'')
!call cla_register('-p','--nproc', 'number of procs to combine', cla_int,'1')
!call cla_register('-t','--ntracer','number of tracers (including T,S)', cla_int  ,'2')
call cla_validate

call cla_get("--iteration",istep)
!call cla_get("--nproc",nproc)
!call cla_get("--ntracer",ntracer)
end subroutine

subroutine combine_hotstart7(istep)
!-------------------------------------------------------------------------------
  use netcdf
  implicit real(8)(a-h,o-z),integer(i-n)
  parameter(nbyte=4)
  character(len=12) :: it_char
  character(len=72) :: fgb,fgb2,fdb,variable_nm
  character(len=72),allocatable :: dim_name(:)
  integer :: lfgb,lfdb       ! Length of processor specific global output file name
  integer :: one_dim,xtype,var1d_dims(1),var2d_dims(2),var3d_dims(3), &
 &data_start_1d(1),data_start_2d(2),data_start_3d(3), &
 &data_count_1d(1),data_count_2d(2),data_count_3d(3)
  integer,allocatable :: dimids(:),idims(:),iwork1(:),iwork2(:), &
 &ner(:),npr(:),nsr(:), & !resident (no ghosts)
 &ielg(:,:),iplg(:,:),islg(:,:)
  real(kind=8),allocatable :: work1(:,:,:),work2(:,:,:)
!-------------------------------------------------------------------------------

!-------------------------------------------------------------------------------
! Aquire user inputs
!-------------------------------------------------------------------------------

  !Debug
  print*, 'istep=',istep

! Read mapping info
  open(10,file='local_to_global_0000',status='old')
  read(10,*)ns_global,ne_global,np_global,nvrt,nproc,ntracers
  close(10)
  !For dim purpose
  if(np_global>ns_global.or.ne_global>ns_global) stop 'ns_global not max'
  allocate(ner(0:nproc-1),npr(0:nproc-1),nsr(0:nproc-1))

  fdb='local_to_global_0000'; fdb=adjustl(fdb)
  lfdb=len_trim(fdb)
  mxner=0; mxnpr=0; mxnsr=0
  do irank=0,nproc-1
    write(fdb(lfdb-3:lfdb),'(i4.4)') irank
    open(10,file=fdb,status='old')
    read(10,*)!ns_global,ne_global,np_global,nvrt,nproc,ntracers

    read(10,*)
    read(10,*)ner(irank); do i=1,ner(irank); read(10,*); enddo;
    read(10,*)npr(irank); do i=1,npr(irank); read(10,*); enddo;
    read(10,*)nsr(irank); do i=1,nsr(irank); read(10,*); enddo;
    close(10)
  enddo !irank
  mxner=maxval(ner)
  mxnpr=maxval(npr)
  mxnsr=maxval(nsr)

  allocate(ielg(0:nproc-1,mxner),iplg(0:nproc-1,mxnpr),islg(0:nproc-1,mxnsr), &
     &stat=istat)
  if(istat/=0) stop 'Allocation error (3)'

  do irank=0,nproc-1
    write(fdb(lfdb-3:lfdb),'(i4.4)') irank
    open(10,file=fdb,status='old')
    read(10,*); read(10,*)
    read(10,*)ner(irank)
    do i=1,ner(irank)
      read(10,*)j,ielg(irank,i)
    enddo
    read(10,*)npr(irank)
    do i=1,npr(irank)
      read(10,*)j,iplg(irank,i)
    enddo
    read(10,*)nsr(irank)
    do i=1,nsr(irank)
      read(10,*)j,islg(irank,i)
    enddo
    close(10)
  enddo !irank

  print*, 'Global quantities:',ne_global,np_global,ns_global

!-------------------------------------------------------------------------------
! Read/write hotstart files
! See schism_step.F90 for rules
!-------------------------------------------------------------------------------
  ! Open file
  write(it_char,'(i12)') istep
  it_char=adjustl(it_char)  !place blanks at end
  it_len=len_trim(it_char)  !length without trailing blanks
  fgb='0000_'//it_char(1:it_len);
  fgb=adjustl(fgb); lfgb=len_trim(fgb);
  !print*, 'suffix is:',fgb(1:lfgb)

  !Query # of vars
  iret=nf90_open('hotstart_0000_'//it_char(1:it_len)//'.nc',OR(NF90_NETCDF4,NF90_NOWRITE),ncid2)
  iret=nf90_inquire(ncid2,nDimensions=ndimensions,nVariables=nvars)
  write(*,*)'nvars=',nvars,ndimensions
  allocate(idims(ndimensions),dim_name(ndimensions),dimids(ndimensions))
  do i=1,ndimensions
    iret=nf90_inquire_dimension(ncid2,i,name=dim_name(i),len=idims(i))
    !write(99,*)'dim=',i,idims(i),dim_name(i)
  enddo !i
  !Get type II arrays (including from modules)
  iret=nf90_inq_varid(ncid2, "time",i);
  iret=nf90_get_var(ncid2,i,time);
  iret=nf90_inq_varid(ncid2, "it",i);
  iret=nf90_get_var(ncid2,i,iths);
  iret=nf90_inq_varid(ncid2, "ifile",i);
  iret=nf90_get_var(ncid2,i,ifile);
  !From modules
  ice_free_flag=-99 !flag
  iret=nf90_inq_varid(ncid2, "ice_free_flag",i);
  if(iret==NF90_NOERR) then
    iret=nf90_get_var(ncid2,i,ice_free_flag);
    if(iret/=NF90_NOERR) stop 'cannot get ice_free_flag'
    print*, 'Found ice_free_flag:',ice_free_flag
  endif
  print*, 'static info:',time,iths,ifile
  iret=nf90_close(ncid2)

  !Open output file
  iret=nf90_create('hotstart_it='//it_char(1:it_len)//'.nc',OR(NF90_NETCDF4,NF90_CLOBBER),ncid)
  iret=nf90_def_dim(ncid,'node',np_global,node_dim)
  iret=nf90_def_dim(ncid,'elem',ne_global,nele_dim)
  iret=nf90_def_dim(ncid,'side',ns_global,nedge_dim)
  do i=4,ndimensions !skip local node/elem/side #s
    iret=nf90_def_dim(ncid,trim(adjustl(dim_name(i))),idims(i),itmp)
    !so that we can easily refer to the dim id in writes
    if(itmp/=i) stop 'pre- and post-comb conflict'
  enddo !i
  !Extra dims
  iret=nf90_def_dim(ncid,'one_new',1,one_dim)
  iret=nf90_enddef(ncid)

  loop1: do m=1,nvars
    do irank=0,nproc-1
      fgb2=fgb
      fgb2=adjustl(fgb2)
      write(fgb2(1:4),'(i4.4)') irank
      iret=nf90_open('hotstart_'//fgb2(1:lfgb)//'.nc',OR(NF90_NETCDF4,NF90_NOWRITE),ncid2)
      if(iret.ne.NF90_NOERR) then
        print*, nf90_strerror(iret); stop
      endif
!      write(99,*)'file=',m,irank,'hotstart_'//fgb2(1:lfgb)//'.nc',ncid2

      iret=nf90_inquire_variable(ncid2,m,name=variable_nm,xtype=xtype,ndims=ndims)
      if(ndims<=0.or.ndims>ndimensions) stop '0 dim encountered'
      variable_nm=adjustl(variable_nm); len_var=len_trim(variable_nm)
!      write(99,*)'var dim:',m,irank,variable_nm(1:len_var),ndims

      if(ndims>ndimensions) stop 'dim overflow'
      iret=nf90_inquire_variable(ncid2,m,dimids=dimids(1:ndims))
      if(iret.ne.NF90_NOERR) then
        print*, nf90_strerror(iret); stop
      endif
      !Re-read dims (rank specific)
      do i=1,ndimensions
        iret=nf90_inquire_dimension(ncid2,i,len=idims(i))
      enddo !i
!      write(99,*)'var dim(2):',m,irank,dimids(1:ndims),(idims(dimids(i)),i=1,ndims)

      if(dimids(ndims)>3) then
        iret=nf90_close(ncid2)
        cycle loop1 !skip static arrays
      endif

      !Define output dims
      if(dimids(ndims)==1) then !node array
        npse=np_global
        npse_dim=node_dim
        npse_lcl=npr(irank)
      else if(dimids(ndims)==2) then !elem array
        npse=ne_global
        npse_dim=nele_dim
        npse_lcl=ner(irank)
      else if(dimids(ndims)==3) then !side array
        npse=ns_global
        npse_dim=nedge_dim
        npse_lcl=nsr(irank)
      else
        stop 'last dim is not node/elem/side'
      endif !dimids

      if(idims(dimids(ndims))/=npse_lcl) then
        print*, 'last dim wrong:',m,irank,npse_lcl,idims(dimids(ndims))
        stop
      endif

      !alloc
      if(ndims==1) then !1D
        itmp1=idims(dimids(1))
        allocate(iwork1(npse_lcl),work1(1,1,itmp1),stat=istat)
        if(irank==0) allocate(iwork2(npse),work2(1,1,npse),stat=istat)
      else if(ndims==2) then !2D array
        itmp1=idims(dimids(1)); itmp2=idims(dimids(2))
        allocate(iwork1(npse_lcl),work1(1,itmp1,itmp2),stat=istat)
        if(irank==0) allocate(iwork2(npse),work2(1,itmp1,npse),stat=istat)
      else if(ndims==3) then
        itmp1=idims(dimids(1))
        itmp2=idims(dimids(2)); itmp3=idims(dimids(3))
        allocate(iwork1(npse_lcl),work1(itmp1,itmp2,itmp3),stat=istat)
        if(irank==0) allocate(iwork2(npse),work2(itmp1,itmp2,npse),stat=istat)
      else
       stop '>3D array'
      endif !ndims
      if(istat/=0) stop 'Allocation error (7)'

      !Check if the var is already defined
      iret2=nf90_inq_varid(ncid,variable_nm(1:len_var),i)

      if(xtype==NF90_INT) then !1D arrays only
        data_start_1d(1)=1; data_count_1d(1)=npse_lcl !npr(irank)
        iret=nf90_get_var(ncid2,m,iwork1(1:npse_lcl),data_start_1d,data_count_1d)

        !Define output var
        if(iret2/=NF90_NOERR) then !not defined yet
          iret=nf90_redef(ncid)
          var1d_dims(1)=npse_dim
          iret=nf90_def_var(ncid,variable_nm(1:len_var),NF90_INT,var1d_dims,ivarid)
          iret=nf90_enddef(ncid)
        endif

        !write(99,*)'Im in:',m,iwork1,ivarid
      else !doubles
        if(ndims==1) then
          data_start_1d(1)=1; data_count_1d(1)=npse_lcl !npr(irank)
          iret=nf90_get_var(ncid2,m,work1(1,1,1:npse_lcl),data_start_1d,data_count_1d)

          !Define output var
          if(iret2/=NF90_NOERR) then !not defined yet
            iret=nf90_redef(ncid)
            var1d_dims(1)=npse_dim
            iret=nf90_def_var(ncid,variable_nm(1:len_var),NF90_DOUBLE,var1d_dims,ivarid)
            iret=nf90_enddef(ncid)
          endif
!          write(99,*)'Im in:',m,irank,ivarid,work1(1,1,1:npse_lcl)

        else if(ndims==2) then !2D
          data_start_2d(1:2)=1
          data_count_2d(1)=idims(dimids(1)); data_count_2d(2)=npse_lcl !npr(irank)
          iret=nf90_get_var(ncid2,m,work1(1,1:data_count_2d(1),1:npse_lcl),data_start_2d,data_count_2d)

          !Define output var
          if(iret2/=NF90_NOERR) then !not defined yet
            iret=nf90_redef(ncid)
            !Pre- and post-comb nc share same dims of ID>3
            var2d_dims(1)=dimids(1); var2d_dims(2)=npse_dim
            iret=nf90_def_var(ncid,variable_nm(1:len_var),NF90_DOUBLE,var2d_dims,ivarid)
            iret=nf90_enddef(ncid)
          endif

        else if(ndims==3) then !3D array
          data_start_3d(1:3)=1
          data_count_3d(1:2)=idims(dimids(1:2)); data_count_3d(3)=npse_lcl !npr(irank)
          iret=nf90_get_var(ncid2,m,work1(1:data_count_3d(1),1:data_count_3d(2),1:npse_lcl), &
     &data_start_3d,data_count_3d)

          !Define output var
          if(iret2/=NF90_NOERR) then !not defined yet
            iret=nf90_redef(ncid)
            var3d_dims(1:2)=dimids(1:2); var3d_dims(3)=npse_dim
            iret=nf90_def_var(ncid,variable_nm(1:len_var),NF90_DOUBLE,var3d_dims,ivarid)
            iret=nf90_enddef(ncid)
          endif
        else
          stop 'unknown array(2)'
        endif !ndims
      endif !xtype

      !Put into global array
      do i=1,npse_lcl !npr(irank)
        if(dimids(ndims)==1) then !node
          ip=iplg(irank,i)
        else if(dimids(ndims)==2) then
          ip=ielg(irank,i)
        else !side
          ip=islg(irank,i)
        endif
        iwork2(ip)=iwork1(i)
        work2(:,:,ip)=work1(:,:,i)
      enddo !i

      iret=nf90_close(ncid2)
      deallocate(iwork1,work1)
    enddo !irank

    !write
    if(xtype==NF90_INT) then !1D arrays only
      data_start_1d(1)=1; data_count_1d(1)=npse
      iret=nf90_put_var(ncid,ivarid,iwork2,data_start_1d,data_count_1d)
    else !doubles
      if(ndims==1) then
        data_start_1d(1)=1; data_count_1d(1)=npse
        iret=nf90_put_var(ncid,ivarid,work2(1,1,1:npse),data_start_1d,data_count_1d)
      else if(ndims==2) then
        data_start_2d(1:2)=1
        !dimids(1) same across all ranks
        data_count_2d(1)=idims(dimids(1)); data_count_2d(2)=npse
        iret=nf90_put_var(ncid,ivarid,work2(1,1:data_count_2d(1),1:npse),data_start_2d,data_count_2d)
      else if(ndims==3) then
        data_start_3d(1:3)=1
        data_count_3d(1:2)=idims(dimids(1:2)); data_count_3d(3)=npse
        iret=nf90_put_var(ncid,ivarid,work2(1:data_count_3d(1),1:data_count_3d(2),1:npse), &
     &data_start_3d,data_count_3d)
      endif !ndims
    endif !xtype

    !Debug
!    if(variable_nm(1:len_var).eq.'eta2') then
!      write(99,*)'Check output:',m,variable_nm(1:len_var)
!      do i=1,npse
!        write(99,*)i,work2(1,1,i)
!      enddo !i

      !iret=nf90_close(ncid)
      !stop
!    endif

    deallocate(iwork2,work2)
  enddo loop1 !m: vars

  !Append type II arrays
  iret=nf90_redef(ncid)
  var1d_dims(1)=one_dim
  iret=nf90_def_var(ncid,'time',NF90_DOUBLE,var1d_dims,i1)
  iret=nf90_def_var(ncid,'iths',NF90_INT,var1d_dims,i2)
  iret=nf90_def_var(ncid,'ifile',NF90_INT,var1d_dims,i3)
  if(ice_free_flag>=0) iret=nf90_def_var(ncid,'ice_free_flag',NF90_INT,var1d_dims,i4)
  iret=nf90_enddef(ncid)
  iret=nf90_put_var(ncid,i1,time)
  iret=nf90_put_var(ncid,i2,iths)
  iret=nf90_put_var(ncid,i3,ifile)
  if(ice_free_flag>=0) iret=nf90_put_var(ncid,i4,ice_free_flag)

  iret=nf90_close(ncid)

end subroutine combine_hotstart7


text:
! @@name:	fort_race.1f
! @@type:	F-free
! @@compilable:	yes
! @@linkable:	no
! @@expect:	rt-error
SUBROUTINE SHARED_RACE

  INCLUDE "omp_lib.h"      ! or USE OMP_LIB

  REAL A(20)
  INTEGER MYTHREAD

!$OMP PARALLEL SHARED(A) PRIVATE(MYTHREAD)

  MYTHREAD = OMP_GET_THREAD_NUM()
  IF (MYTHREAD .EQ. 0) THEN
     CALL SUB(A(1:10)) ! compiler may introduce writes to A(6:10)
  ELSE
     A(6:10) = 12
  ENDIF

!$OMP END PARALLEL

END SUBROUTINE SHARED_RACE

SUBROUTINE SUB(X)
  REAL X(*)
  X(1:5) = 4
END SUBROUTINE SUB

ext: f

max_stars_repo_name: kammerdienerb/flang

max_issues_repo_name: kammerdienerb/flang

max_forks_repo_name: kammerdienerb/flang


module tachedata

  !-----------------------------------------------------------------------
  ! Data module for input parameters to TACHE and tensors
  ! PJC 20/05/2008
  ! DHF 01/09/2009
  ! DHF 26/02/2015 to incorporate modern sphNG output
  !-----------------------------------------------------------------------

  implicit none
  save

  !-------------------Single values---------------------------------------
  ! Integers
  integer :: nfiles

  integer :: ielement, nelement, neigen

  ! Reals
  real(kind=8) :: threshold

  ! Characters
  character(12),parameter :: paramfile='tache.params'
  character(100) :: listfile, filetype,fileformat,tensorchoice
  character(1) :: splitdumpchoice, tensorchar

  ! Parameter decides how neighbours are calculated:
  ! 'g' - particle neighbours found using regular grid
  ! 'o' - neighbours found by building an octree (O(N log N))
  ! 'b' - neighbours found by brute force (O(N^2))

  character(1),parameter :: use_octree_grid = 'b'

  ! Parameter decides how gravitational force is calculated
  ! 'p' - if potential exists, take its derivative
  ! 'o' - use an octree
  ! 'b' - use brute force (on neighbour list)

  character(1),parameter :: grav_calc_choice = 'b'

  ! When splitting SPH dump into constituents, recalculate element densities?
  character(1),parameter :: density_recalc = 'n'

  character(100),allocatable,dimension(:) :: filename,gravfile,potfile
  character(100),allocatable,dimension(:) :: eigenfile,vectorfile,memberfile

! Logicals
  logical :: splitdump

!-------------------Array values----------------------------------------
! Integers
  ! Reals

  integer,allocatable,dimension(:) :: class
  integer,dimension(4) :: classnum

  integer,parameter:: nclasses=4
  character(4),dimension(nclasses) :: prefixes

  integer,allocatable,dimension(:) :: eigenelement
  real,allocatable,dimension(:,:) :: eigenvalues
  real,allocatable,dimension(:,:,:) :: eigenvectors
  real, allocatable,dimension(:,:,:) :: tensor

end module tachedata


text:
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C       Written by: Anshuk Attri
C
C       Contact: contact@anshukattri.in
C       Website: www.anshukattri.in/research
C       GITHUB: github.com/aa-tree/
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

        SUBROUTINE START_TWOBODYSYSTEM(VBODYI)
        IMPLICIT NONE

        INCLUDE 'common.dat'

        INTEGER VBODYI
        LOGICAL BFOUND

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C       SUMMARY:
C
C       This is the routine that must be called in case you are working
C       with a two body system.
C
C       VBODYI                  INTEGER INPUT. It is the SpiceLib ID of
C                               the body that is to be used as a primary
C                               body in a two body system. For example,
C                               for integration of a satellite around
C                               Earth, you must pass in 399.
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

        INTEGER I

        BFOUND = .FALSE.

        DO I=1,MMAX

            IF (LIST_I(I,1) .EQ. VBODYI) THEN

                IPB = VBODYI            !Initiate common variable IPB.
                BFOUND=.TRUE.
                SYS_TYPE=2

            END IF
        END DO

        IF(BFOUND .EQV. .FALSE.) THEN
            WRITE (*,*) 'Error in START_TWOBODYSYSTEM:'
            WRITE (*,*) 'Two body system initialisation error!'
            WRITE (*,*) 'Primary body supplied is not in the model.'
            WRITE (*,*) 'Killing program.'

            STOP
         END IF

        RETURN
        END

max_stars_repo_name: aa-tree/SPICEReader

max_issues_repo_name: aa-tree/SPICEReader

max_forks_repo_name: aa-tree/SPICEReader


! RUN: %S/test_errors.sh %s %t %flang_fc1
! Test SELECT TYPE and ASSOCIATE errors: C1103

subroutine s1()
  class(*),allocatable :: calc[:]
  integer,save :: icoa[*]
  !ERROR: Selector must not be a coindexed object
  associate(sel=>icoa[2])
  end associate
  icoa = 2
  allocate(integer::calc[*])
  !ERROR: Selector must not be a coindexed object
  select type(sel=>calc[2])
  end select
end subroutine


text:
C
C $Id: mdglty.f,v 1.2 2008-07-27 00:17:06 haley Exp $
C
C                Copyright (C)  2000
C        University Corporation for Atmospheric Research
C                All Rights Reserved
C
C The use of this Software is governed by a License Agreement.
C
      SUBROUTINE MDGLTY (IARG)
C
        INTEGER IARG
C
C This routine may be called by a user from the line-processing routine
C specified as the final argument in a call to MDLNDM in order to obtain
C the type of the line being processed.
C
        COMMON /MAPCMZ/  NNMS,ILTY,IAIL,IAIR,BLAG,SLAG,BLOG,SLOG,
     +                   PNTS(200)
        INTEGER          NNMS,ILTY,IAIL,IAIR
        REAL             BLAG,SLAG,BLOG,SLOG,PNTS
        SAVE   /MAPCMZ/
C
        IARG=ILTY
C
        RETURN
C
      END

max_stars_repo_name: tenomoto/ncl

max_issues_repo_name: tenomoto/ncl

max_forks_repo_name: tenomoto/ncl


subroutine  read_goesimgr_skycover(nread,ndata,nodata,infile,obstype,lunout,gstime,twind,sis,&
                        prsl_full,nobs)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    read_goesimgr_skycover                    read GOES Imager sky cover product
!   prgmmr: Jacob Carley                               date: 2014-11-07
!
! abstract:  This routine reads GOES Imager sky cover data from bufr_d files.
!            It also has options to thin the data by using conventional
!            thinning programs, though data are only thinned in 2D and not 3D
!            since sky cover is a 2D data set.
!            When running the gsi in regional mode, the code only
!            retains those observations that fall within the regional
!            domain
!
! program history log:
!   2014-11-07 J. Carley - Initial code
!   2015-03-06 C. Thomas - Add thin4d logical for removal of time thinning
!   2016-01-11 D. Keyser - Enable use of efclam dump as the primary observation
!                          source file. However, if ob is missing, will look for it in old
!                          BUFR mnemonic sequence
!   2016-03-11 j. guo    - Fixed {dlat,dlon}_earth_deg in the obs data stream
!   2016-04-22 M. Pondeca  - Replace "if (goescld(3)==bmiss)" condition with "if (goescld(3) > r0_01_bmiss)"
!
!   input argument list:
!     ithin    - flag to thin data
!     rmesh    - thinning mesh size (km)
!     gstime   - analysis time in minutes from reference date
!     infile   - unit from which to read BUFR data
!     lunout   - unit to which to write data for further processing
!     obstype  - observation type to process
!     twind    - input group time window (hours)
!     sis      - satellite/instrument/sensor indicator
!
!   output argument list:
!     nread    - number of obs read
!     ndata    - number of obs retained for further processing
!     nodata   - number of obs retained for further processing
!     nobs     - array of observations on each subdomain for each processor
!
! attributes:
!   language: f95/2003
!   machine:  WCOSS
!
!$$$

  use kinds, only: r_single,r_kind,r_double,i_kind
  use constants, only: zero,one_tenth,one,deg2rad,half,&
      three,four, r60inv,r10,r100,r2000

  use convinfo, only: nconvtype, &
      icuse,ictype,ioctype,&
      ithin_conv,rmesh_conv,pmesh_conv,ctwind
  use convthin, only: make3grids,map3grids,del3grids,use_all
  use gridmod, only: regional,nlon,nlat,nsig,tll2xy,txy2ll,&
      rlats,rlons
  use deter_sfc_mod, only: deter_sfc2
  use obsmod, only: bmiss,ran01dom
  use gsi_4dvar, only: l4dvar,l4densvar,iwinbgn,winlen,time_4dvar,thin4d
  use adjust_cloudobs_mod, only: adjust_goescldobs
  use mpimod, only: npe

  implicit none

! Declare passed variables
  character(len=*)                      ,intent(in   ) :: infile,obstype
  character(len=20)                     ,intent(in   ) :: sis
  integer(i_kind)                       ,intent(in   ) :: lunout
  integer(i_kind)                       ,intent(inout) :: nread,ndata,nodata
  integer(i_kind),dimension(npe)        ,intent(inout) :: nobs
  real(r_kind)                          ,intent(in   ) :: twind,gstime
  real(r_kind),dimension(nlat,nlon,nsig),intent(in   ) :: prsl_full

! Declare local parameters
  real(r_kind),parameter:: r90  = 90.0_r_kind
  real(r_kind),parameter:: r0_01 = 0.01_r_kind
  real(r_kind),parameter:: r0_1_bmiss=one_tenth*bmiss
  real(r_kind),parameter:: r0_01_bmiss=r0_01*bmiss
  real(r_kind),parameter:: r1200= 1200.0_r_kind
  real(r_kind),parameter:: r6= 6.0_r_kind
  real(r_kind),parameter:: r360 = 360.0_r_kind
  character(8),parameter:: cspval= '88888888'

! Declare local variables
  character(len=80) :: hdrstr,goescldstr,goescldstr_new
  character(len=8) ::  subset
  character(len=22) :: myname
  character(len=8) :: c_prvstg,c_sprvstg ,c_station_id

  integer(i_kind) :: nmsub,ireadmg,ireadsb,nreal,nc,i,lunin,nmsg,ntb
  integer(i_kind) :: iret,kx,pflag,nlevp,nmind,levs,idomsfc
  integer(i_kind) :: low_cldamt_qc,mid_cldamt_qc,hig_cldamt_qc,tcamt_qc
  integer(i_kind) :: ithin,klat1,klon1,klonp1,klatp1,kk,k,ilat,ilon,nchanl
  integer(i_kind) :: iout,ntmp,iiout,maxobs,icount,itx,iuse,idate,ierr
  integer(i_kind),dimension(5) :: idate5
  integer(i_kind),allocatable,dimension(:):: isort,iloc
  real(r_kind) :: dlat,dlon,dlat_earth,dlon_earth,toff,t4dv
  real(r_kind) :: dlat_earth_deg,dlon_earth_deg
  real(r_kind) :: dx,dx1,dy,dy1,w00,w10,w01,w11,crit1,timedif,tdiff
  real(r_kind) :: rmesh,pmesh,xmesh,tcamt,tcamt_oe,ff10,tsavg
  real(r_kind) :: rminobs,ppb
  real(r_kind) :: low_cldamt,mid_cldamt,hig_cldamt,usage,zz,sfcr,rstation_id
  real(r_kind),allocatable,dimension(:):: presl_thin
  real(r_kind),dimension(nsig):: presl
  real(r_kind),allocatable,dimension(:,:):: cdata_all,cdata_out
  real(r_double),dimension(9):: hdr
  real(r_double),dimension(3):: goescld

  logical :: outside,ithinp,luse

  real(r_double),dimension(1,1):: r_prvstg,r_sprvstg

!  equivalence to handle character names
  equivalence(rstation_id,c_station_id)

  lunin=11_i_kind
  ithin=-9_i_kind
  rmesh=-99.999_r_kind
  myname='READ_GOESIMGR_SKYCOVER'
  hdrstr='NUL  YEAR MNTH DAYS HOUR MINU SECO CLATH CLONH'
  goescldstr='SAID TOCC TOCC_AVG'
  goescldstr_new='SAID ECAS ECAM'
  nreal=20

  nc=0
  conv: do i=1,nconvtype
     if(trim(obstype) == trim(ioctype(i)) .and. ictype(i)==154_i_kind) then
        nc=i
        exit conv
     end if
  end do conv
  if(nc == 0)then
     write(6,*) myname,' no matching obstype found in convinfo ',obstype
     return
  end if

  ! Try opening bufr file, if unable print error to the screen
  !  and return to read_obs.F90
  open(lunin,file=trim(infile),form='unformatted',iostat=ierr)
  if (ierr/=0) then
     write(6,'(A)')myname,':ERROR: Trouble opening input file: ',trim(infile),' returning to read_obs.F90...'
     return
  end if

  call openbf(lunin,'IN',lunin)
  call datelen(10)

  ! Set up thinning params
  use_all = .true.
  ithin=ithin_conv(nc)
  if (ithin > 0 ) then
     rmesh=rmesh_conv(nc)
     pmesh=pmesh_conv(nc)
     use_all = .false.
     if(pmesh > zero) then
        pflag=1
        nlevp=r1200/pmesh
     else
        pflag=0
        nlevp=nsig
     endif
     pflag=0
     nlevp=nsig
     xmesh=rmesh
     call make3grids(xmesh,nlevp)
     if (.not.use_all) then
        allocate(presl_thin(nlevp))
        if (pflag==1) then
           do k=1,nlevp
              presl_thin(k)=(r1200-(k-1)*pmesh)*one_tenth
           enddo
        endif
     endif
     write(6,'(A,1x,A,1x,A,I4,1x,f8.2,1x,I2,1x,I3,1x,f8.2,1x,I3)')myname,': ioctype(nc),ictype(nc),rmesh,pflag,nlevp,pmesh,nc ',&
                 trim(ioctype(nc)),ictype(nc),rmesh,pflag,nlevp,pmesh,nc
  endif

  nmsg=0
  ntb = 0
  ! Find number of reports and messages so we know how much memory to allocate
  do while (ireadmg(lunin,subset,idate) == 0)
!    Time offset
     if(nmsg == 0) call time_4dvar(idate,toff)
     nmsg=nmsg+1
     ntb = ntb + nmsub(lunin) !nmsub is a bufrlib function which returns the number of subsets in
                              !  a bufr message open for input via a previous call to a bufrlib
                              !  routine readmg or equivalent.  The subsets are not required to be read (saves time).
  end do
  maxobs=ntb

  allocate(cdata_all(nreal,maxobs),isort(maxobs))
  isort = 0
  cdata_all=zero
  nread=0
  nchanl=0
  ilon=2
  ilat=3
  ntb=0

   close(lunin)
   call closbf(lunin)
   open(lunin,file=trim(infile),form='unformatted')
   call openbf(lunin,'IN',lunin)
   call datelen(10)

   loop_msg: do while (ireadmg(lunin,subset,idate) == 0)
         loop_readsb: do while (ireadsb(lunin) == 0)
            ntb=ntb+1
            ! - Extract type, date, and location information
            call ufbint(lunin,hdr,9,1,iret,hdrstr)

            ! - Compare relative obs time with window.  If obs
            ! -  falls outside of window, don't use this obs
            idate5(1) = hdr(2)     ! year
            idate5(2) = hdr(3)     ! month
            idate5(3) = hdr(4)     ! day
            idate5(4) = hdr(5)     ! hours
            idate5(5) = hdr(6)     ! minutes
            call w3fs21(idate5,nmind)
            rminobs=real(nmind,8)+(real(hdr(7),8)*r60inv)!convert the seconds of the ob to minutes and store to rminobs
            t4dv = (rminobs-real(iwinbgn,r_kind))*r60inv
            tdiff=(rminobs-gstime)*r60inv  !GS time is the analysis time in minutes from w3fs21

            if (l4dvar.or.l4densvar) then
               if (t4dv<zero .OR. t4dv>winlen) cycle loop_readsb
            else
               ! - Check to make sure ob is within convinfo time window (ctwind) and
               ! -  is within overwall time window twind (usually +-3)
               if( (abs(tdiff) > ctwind(nc)) .or. (abs(tdiff) > twind) )cycle loop_readsb
            endif

            kx=999_i_kind !hardwire typ to 999
            if(abs(hdr(8))>r90 .or. abs(hdr(9))>r360) cycle loop_readsb
            if(hdr(9)== r360)hdr(9)=hdr(9)-r360
            if(hdr(9) < zero)hdr(9)=hdr(9)+r360
            dlon_earth_deg = hdr(9)
            dlat_earth_deg = hdr(8)
            dlon_earth=hdr(9)*deg2rad
            dlat_earth=hdr(8)*deg2rad
            nread=nread+1
            if(regional)then
               call tll2xy(dlon_earth,dlat_earth,dlon,dlat,outside)    ! convert to rotated coordinate
               if(outside) cycle loop_readsb   ! check to see if outside regional domain
            else
               dlat = dlat_earth
               dlon = dlon_earth
               call grdcrd1(dlat,rlats,nlat,1)
               call grdcrd1(dlon,rlons,nlon,1)
            endif

            ! Read in the obs
            goescld=bmiss
            call ufbint(lunin,goescld,3,1,levs,goescldstr_new)
            if (goescld(3) > r0_01_bmiss) then
! if ob is missing, look for it in old BUFR mnemonic sequence
               goescld=bmiss
               call ufbint(lunin,goescld,3,1,levs,goescldstr)
               if (goescld(3) > r0_01_bmiss) cycle loop_readsb !If obs are missing, cycle
            endif
            c_prvstg=cspval
            c_sprvstg=cspval

            ! -  Set station ID
            rstation_id=goescld(1)

            ithin=ithin_conv(nc)
            ithinp = ithin > 0 .and. pflag /= 0

          ! - Thin in vertical  - note we can only thin in the horizontal
          ! -   since sky cover is a 2D field.  So this branch should never run
          ! -   unless we get info about the vertical location of the clouds in the
          ! -   future.  Leaving here as a 'just-in-case' measure.
            if(ithinp   )then
!           Interpolate guess pressure profile to observation location
               klon1= int(dlon);  klat1= int(dlat)
               dx   = dlon-klon1; dy   = dlat-klat1
               dx1  = one-dx;     dy1  = one-dy
               w00=dx1*dy1; w10=dx1*dy; w01=dx*dy1; w11=dx*dy

               klat1=min(max(1,klat1),nlat); klon1=min(max(0,klon1),nlon)
               if (klon1==0) klon1=nlon
               klatp1=min(nlat,klat1+1); klonp1=klon1+1
               if (klonp1==nlon+1) klonp1=1
               do kk=1,nsig
                  presl(kk)=w00*prsl_full(klat1 ,klon1 ,kk) +  &
                            w10*prsl_full(klatp1,klon1 ,kk) + &
                            w01*prsl_full(klat1 ,klonp1,kk) + &
                            w11*prsl_full(klatp1,klonp1,kk)
               end do
            end if

            iuse=icuse(nc)

            ! General block for data thinning - if requested
            if (ithin > 0 .and. iuse >=0) then
               ntmp=ndata  ! counting moved to map3gridS
            ! - Set data quality index for thinning
               if (thin4d) then
                  timedif = zero
               else
                  timedif=abs(t4dv-toff)
               endif

               crit1 = timedif/r6+half

               ! - simple 1-to-1 mapping of vertical levels when no thinning in the vertical
               if (pflag==0) then
                  do kk=1,nsig
                     presl_thin(kk)=presl(kk)
                  end do
               endif
               ppb=one_tenth*1013.25_r_kind !number is irrelevant for 2D - set to standard SLP -> 1013.25 and convert from mb to cb
               call map3grids(-1,pflag,presl_thin,nlevp,dlat_earth,dlon_earth,&
                                 ppb,crit1,ndata,iout,ntb,iiout,luse,.false.,.false.)

               if (.not. luse) cycle loop_readsb
               if(iiout > 0) isort(iiout)=0
               if (ndata > ntmp) then
                  nodata=nodata+1
               endif
               isort(ntb)=iout
            else  ! - no thinnning
               ndata=ndata+1
               nodata=nodata+1
               iout=ndata
               isort(ntb)=iout
         endif

         !-  Set usage variable
         usage = 0
         if(iuse <= 0)usage=r100

         ! Get information from surface file necessary for conventional data here
         call deter_sfc2(dlat_earth,dlon_earth,t4dv,idomsfc,tsavg,ff10,sfcr,zz)

         ! - Obtain the ob and tune the QC marks for ob error tuning a bit later
         call adjust_goescldobs(goescld(3),tdiff,dlat_earth,dlon_earth, &
                                low_cldamt,low_cldamt_qc,mid_cldamt,mid_cldamt_qc, &
                                hig_cldamt,hig_cldamt_qc,tcamt,tcamt_qc)

         if(tcamt_qc==15 .or. tcamt_qc==12 .or. tcamt_qc==9) usage=r100
         tcamt_oe=20.0_r_kind
         if(tcamt_qc==1) tcamt_oe=tcamt_oe*1.25_r_kind
         if(tcamt_qc==2) tcamt_oe=tcamt_oe*1.50_r_kind
         if(tcamt_qc==3) tcamt_oe=tcamt_oe*1.75_r_kind

         cdata_all( 1,iout)=tcamt_oe               !  obs error
         cdata_all( 2,iout)=dlon                   !  grid relative longitude
         cdata_all( 3,iout)=dlat                   !  grid relative latitude
         cdata_all( 4,iout)=tcamt                  !  total cloud amount (%)
         cdata_all( 5,iout)=rstation_id            !  station ID
         cdata_all( 6,iout)=t4dv                   !  time
         cdata_all( 7,iout)=nc                     !  type
         cdata_all( 8,iout)=tcamt_qc               !  quality mark
         cdata_all( 9,iout)=usage                  !  usage parameter
         cdata_all(10,iout)=idomsfc                !  dominate surface type
         cdata_all(11,iout)=tsavg                  !  skin temperature
         cdata_all(12,iout)=ff10                   !  10 meter wind factor
         cdata_all(13,iout)=sfcr                   !  surface roughness
         cdata_all(14,iout)=dlon_earth_deg         !  earth relative longitude (degrees)
         cdata_all(15,iout)=dlat_earth_deg         !  earth relative latitude (degrees)
         cdata_all(16,iout)=bmiss                  !  station elevation (m)
         cdata_all(17,iout)=bmiss                  !  observation height (m)
         cdata_all(18,iout)=zz                     !  terrain height at ob location
         cdata_all(19,iout)=r_prvstg(1,1)          !  provider name
         cdata_all(20,iout)=r_sprvstg(1,1)         !  subprovider name

        enddo  loop_readsb

     enddo loop_msg

!    Close unit to bufr file
     call closbf(lunin)
!    Deallocate arrays used for thinning data
     if (.not.use_all) then
        deallocate(presl_thin)
        call del3grids
     endif
! Normal exit

! Write header record and data to output file for further processing
  allocate(iloc(ndata))
  icount=0
  do i=1,maxobs
     if(isort(i) > 0)then
       icount=icount+1
       iloc(icount)=isort(i)
     end if
  end do
  if(ndata /= icount)then
     write(6,*) myname,': ndata and icount do not match STOPPING...ndata,icount ',ndata,icount
     call stop2(50)
  end if
  allocate(cdata_out(nreal,ndata))
  do i=1,ndata
     itx=iloc(i)
     do k=1,nreal
        cdata_out(k,i)=cdata_all(k,itx)
     end do
  end do
  deallocate(iloc,isort,cdata_all)

  call count_obs(ndata,nreal,ilat,ilon,cdata_out,nobs)
  write(lunout) obstype,sis,nreal,nchanl,ilat,ilon,ndata
  write(lunout) cdata_out

  deallocate(cdata_out)

  if (ndata == 0) then
     write(6,*)myname,':  closbf(',lunin,')'
  endif

  close(lunin)

end subroutine  read_goesimgr_skycover


text:
      SUBROUTINE SHIRBL(NAME,NNAM)
C
C        SHIFT DATA 16 PLACES RIGHT AND ZERO FILL
C
      PARAMETER(WORDSZ=32)
      DIMENSION NAME(NNAM)
      NUMSHFT=WORDSZ-16
      DO 100 I=1,NNAM
         NAME(I)=ICEDSHFT(NAME(I),-NUMSHFT)
C         NAME(I)=RSHIFT(NAME(I),NUMSHFT)
 100  CONTINUE
      RETURN
      END


text:
      subroutine SVECPR(V, N, TEXT, LWIDTH, LUNIT, NUMDIG)
c Copyright (c) 1996 California Institute of Technology, Pasadena, CA.
c ALL RIGHTS RESERVED.
c Based on Government Sponsored Research NAS7-03001.
c>> 2010-02-22 SVECPR Krogh Initliazed TTEXT to avoid compiler worrirs.
C>> 2000-12-01 SVECPR Krogh Removed unused parameter METXTF.
C>> 1994-11-11 SVECPR Krogh Declared all vars.
c>> 1994-10-20 SVECPR Krogh Added M77CON code.
c>> 1992-05-03 SVECPR Krogh Convert to use MESSFT for Fortran text.
c>> 1992-04-08 SVECPR Krogh Replaced dummy K with MACT in ?MESS call.
c>> 1991-11-22 SVECPR F. Krogh   Initial code
c++I DEFAULT MACT1D=3, TAIL=", LUNIT"
c++  DEFAULT MACT1D=5, TAIL=", LUNIT, NUMDIG"
c++  REPLACE ", LUNIT, NUMDIG" = TAIL
c--S replaces "?": ?VECPR, ?MESS
c
c     *****     Formal Arguments     ***********************************
c
c V      Vector to be output, V = V(I), I = 1, N
c N      Number of vector components.
c TEXT   a variable length character type that gives a message to print.
c LWIDTH Line width in characters.  If this or any of the following
c        parameter are < 0, then current defaults set in MESS are used.
c LUNIT  Logical unit number.  (0 prints to the standard output.)
c NUMDIG Number of significant digits to print (Not used for integer).
c
c  ******************** Parameter for interfacing to MESS  *************
c
      integer MEMLIN, MEMUNI, METDIG, MECONT, MERET, MACT1D
      parameter (MEMLIN =13)
      parameter (MEMUNI =15)
      parameter (METDIG =22)
      parameter (MECONT =50)
      parameter (MERET  =51)
c++ Substitute for MACT1D below
      parameter (MACT1D = 5)
c
      integer           N, LWIDTH, LUNIT, NUMDIG
      character*(*)     TEXT
      character*2       TTEXT(1)
      data TTEXT(1) / '  ' /
      integer           K, MACT(9), MACT1(MACT1D)
c--S Next line special: I
      real              V(N)
c++  Code for {I} is inactive
C      integer MEIVEC
C      parameter (MEIVEC=57)
C      data MACT1 / MEIVEC, 0, MERET /
C      MACT1(2) = max(N, 0)
c++  Code for ~{I} is active
      integer MEFVEC
      parameter (MEFVEC=61)
      data MACT1 / METDIG, 0, MEFVEC, 0, MERET /
      MACT1(2) = max(0, NUMDIG)
      MACT1(4) = max(N, 0)
c++  End
c
      K = 1
      if (LWIDTH .gt. 20) then
         MACT(1) = -MEMLIN
         MACT(3) = MEMLIN
         MACT(4) = LWIDTH
         K = 5
      end if
      if (LUNIT .ge. 0) then
         MACT(K) = -MEMUNI
         MACT(K+2) = MEMUNI
         MACT(K+3) = LUNIT
         K = K + 4
      end if
      MACT(K) = MECONT
      if (K .gt. 1) call MESSFT(MACT, TEXT)
c++  Code for {I} is inactive
C      call MESS(MACT1, TTEXT, V)
c++  Code for ~{I} is active
      call SMESS(MACT1, TTEXT, MACT, V)
c++  End
      if (MACT(1) .lt. 0) then
c                           Restore MESS parameters to original state
         MACT(1) = - MACT(1)
         MACT(3) = MERET
         if (MACT(5) .lt. 0) then
            MACT(3) = MEMUNI
            MACT(4) = MACT(6)
            MACT(5) = MERET
         end if
         call MESS(MACT, TTEXT, MACT)
      end if
      return
      end


!> This module contains a function to perform unit conversions
MODULE unit_conversions
  IMPLICIT NONE
  PRIVATE
  PUBLIC convert

CONTAINS

    !> Function that converts units from "a" to "b"
  FUNCTION convert(a,b) RESULT(res)
      USE parameters
    USE input_output
    IMPLICIT NONE
    CHARACTER(len=*), INTENT(in) :: a,b
    REAL(dp) :: res

    ! Femtoseconds converted to atomic units of time
    IF ( a == 'fs' .AND. b == 'au_time' ) THEN
       res = fs2s * s2autime
    ELSE IF ( a == 'au_time' .AND. b == 'fs' ) THEN
       res = 1.0_dp / ( fs2s * s2autime )

    ! Kelvin converted to atomic units of energy
    ELSE IF ( a == 'kelvin' .AND. b == 'au_energy' ) THEN
       res = kb * joules2au
    ELSE IF ( a == 'au_energy' .AND. b == 'kelvin' ) THEN
       res = 1.0_dp / ( kb * joules2au )

    ! Wavenumbers converted to atomic units of angular frequency
    ELSE IF ( a == 'wvnbr' .AND. b == 'au_ang_freq' ) THEN
       res = wvnbr2Hz/s2autime*2.0_dp*pi
    ELSE IF ( a == 'au_ang_freq' .AND. b == 'wvnbr' ) THEN
       res = 1.0_dp/(wvnbr2Hz/s2autime*2.0_dp*pi)

    ! Wavenumbers converted to atomic units of energy
    ELSE IF ( a == 'wvnbr' .AND. b == 'au_energy' ) THEN
       res = wvnbr2joules*joules2au
    ELSE IF ( a == 'au_energy' .AND. b == 'wvnbr' ) THEN
       res = 1.0_dp/(wvnbr2joules*joules2au)

    ELSE
       OPEN(UNIT=10, FILE=ERRORLOG)
       WRITE(10,*) 'Undefined unit conversion attempted, see file unit_conversions.f90'
       CLOSE(10)
    END IF

  END FUNCTION convert

END MODULE unit_conversions


text:
	SUBROUTINE IM_DROP  ( iret )
C************************************************************************
C* IM_DROP								*
C*									*
C* This subroutine drops the requested image.  It has to be called 	*
C* after IM_SIMG is called. 						*
C*									*
C* IM_DROP  ( IRET )							*
C*									*
C* Output parameters:							*
C*	IRET		INTEGER		Return code			*
C*					  0 = normal return		*
C**									*
C* Log:									*
C* J. Cowie/COMET	 2/95						*
C* C. Lin/EAI	 	 6/95						*
C* C. Lin/EAI	 	 9/95	Add ER_WMSG after GSATIM		*
C* J. Cowie/COMET	10/95	Fixed arg in call to ER_WMSG()		*
C* C. Lin/EAI	         1/96	pass the error back if GSATIM fails	*
C* J. Cowie/COMET	 7/96	Updated color bar to use calibration	*
C* S. Jacobs/NCEP	 8/96	Commented call to IM_CBAR, for now	*
C* J. Cowie/COMET	10/96	Fixed IFINVD typo, minpv fix for brit	*
C* J. Cowie/COMET	 1/97	Changed IMGDEF common variable names	*
C* S. Jacobs/NCEP	 4/97	Removed call to IM_RTBL			*
C* J. Cowie/COMET	12/97	Replaced imsorc check with imradf	*
C* S. Jacobs/NCEP	12/97	Added color bar for radar images only	*
C* D.W.Plummer/NCEP	 3/03	Change calling seq of IM_BTOT & IM_TTOB	*
C* M. Li/SAIC		11/03	Removed IM_CBAR				*
C* M. Li/SAIC		11/03	Color bar value plotted as temperature	*
C*				only for imtype = 8 or 128		*
C* S. Chiswell/Unidata	 2/04	Added imcmn.cmn for calibration bar flag*
C* D.W.Plummer/NCEP	 7/04	Reverse order of IM_ICMN and GSICMN	*
C* S. Chiswell/Unidata	 7/04	Added plotting of additional NIDS parts *
C* T. Piper/SAIC	07/06	Moved color bar legend to im_labl	*
C************************************************************************
	INCLUDE		'GEMPRM.PRM'
	INCLUDE		'IMGDEF.CMN'
	INCLUDE		'imcmn.cmn'
C
	CHARACTER	nidsfile*132
C*
C------------------------------------------------------------------------
	iret = 0
C
C*	Check if IM_SIMG was called for seting up the navigation
C
	IF ( imftyp .eq. IFINVD ) RETURN
C
C*      Set image common information in APPL imgdef.h
C*      Send through first entry in COMMON /IMGDEF/ which is 'imftyp'.
C*      The imgdef.h information is necessary for IM_GVTOTA
C
	CALL IM_ICMN ( imftyp, iret )
	IF ( iret .ne. 0 ) RETURN
C
C*      Set image common information in GEMPLT
C
	CALL GSICMN ( iret )
	IF ( iret .ne. 0 ) RETURN
C
C*	Set non-image start flag if NIDS level III
C
	istart = 1
	IF ( ( imftyp .eq. IFNIDS ) .or. ( imftyp .eq. IFNEXZ ) ) THEN
	    IF ( imrdfl .ge. 2 ) istart = 0
	ENDIF
C
C*	Drop the image
C
	IF ( istart .gt. 0 ) CALL GSATIM ( cmfile, iret )
	IF ( iret .ne. 0 ) THEN
	    CALL ER_WMSG ( 'GEMPLT', iret, ' ', ierr )
	    RETURN
	ENDIF
C
C*	If a NIDS level III product, plot the auxilary parts- if any
C
	IF ( ( imftyp .eq. IFNIDS ) .or. ( imftyp .eq. IFNEXZ ) ) THEN
	    IF ( imnblk .ge. ( istart + 3 ) ) THEN
		CALL ST_NULL  ( cmfile, nidsfile, lens, iret )
		CALL im_nidsdriver ( nidsfile, imftyp, imprsz,
     +			imlenf, imnblk, imboffs, istart )
	    ENDIF
	ENDIF
	CALL IM_LABL()
C
	RETURN
	END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


text:
!***********************************************************************
!  Integrated Water Flow Model (IWFM)
!  Copyright (C) 2005-2021
!  State of California, Department of Water Resources
!
!  This program is free software; you can redistribute it and/or
!  modify it under the terms of the GNU General Public License
!  as published by the Free Software Foundation; either version 2
!  of the License, or (at your option) any later version.
!
!  This program is distributed in the hope that it will be useful,
!  but WITHOUT ANY WARRANTY; without even the implied warranty of
!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!  GNU General Public License for more details.
!  (http://www.gnu.org/copyleft/gpl.html)
!
!  You should have received a copy of the GNU General Public License
!  along with this program; if not, write to the Free Software
!  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
!
!  For tecnical support, e-mail: IWFMtechsupport@water.ca.gov
!***********************************************************************
PROGRAM BUDGET
  USE MessageLogger    , ONLY: SetLogFileName , &
                               LogMessage     , &
                               LogLastMessage , &
                               PrintRunTime   , &
                               FILE           , &
                               MessageArray   , &
                               iMessage
  USE IWFM_Core_Version
  USE BudgetControls
  USE GeneralUtilities
  USE ProgramTimer     , ONLY: StartTimer
  USE IOInterface
  IMPLICIT NONE

  !Local variables
  TYPE(GenericFileType) :: DummyFile !Dummy file to re-set the standard output file
  INTEGER               :: iStat

  !Start CPU timer
  CALL StartTimer()

  !Refresh the standard output file
  CALL DummyFile%New(FileName='BudgetMessages.out',InputFile=.FALSE.,iStat=iStat)
  IF (iStat .EQ. -1) CALL EndExecution(iStat)
  CALL DummyFile%Kill(Status='DELETE')

  !Set the file for messages
  CALL SetLogFileName('BudgetMessages.out',iStat)
  IF (iStat .EQ. -1) CALL EndExecution(iStat)

  !Print-out date and time of the execution
  MessageArray(1) = 'PROGRAM: IWFM Budget Post-Processor ' // TRIM(IWFM_Core%GetVersion())
  MessageArray(2) = 'This run is made on '//TRIM(GetDate())//' at '//TRIM(GetTime())
  MessageArray(3) = ''
  CALL LogMessage(MessageArray(1:3),iMessage,'',Destination=FILE)

  !Print budget tables
  CALL PrintBudgetTables(iStat=iStat)

  !End execution
  CALL EndExecution(iStat)

END

max_stars_repo_name: giorgk/IWFM-SAFE

max_issues_repo_name: giorgk/IWFM-SAFE

max_forks_repo_name: giorgk/IWFM-SAFE


text:
*----------------------------------------------------------------------*
      subroutine strip_contr(contr,mode)
*----------------------------------------------------------------------*
*     mode = 1:  remove (x 0 [C]) and (0 x [C]) proto-arcs
*     mode = .... to come
*----------------------------------------------------------------------*
      implicit none

      include 'opdim.h'
      include 'def_contraction.h'
      include 'ifc_operators.h'

      type(contraction), intent(inout) ::
     &     contr
      integer, intent(in) ::
     &     mode

      integer ::
     &     narc, narc_new, iarc
      type(cntr_arc), pointer ::
     &     arc(:)

      narc = contr%narc
      arc => contr%arc

      select case (mode)
      case(1)
        narc_new = 0
        do iarc = 1, narc
          if (arc(iarc)%link(1).eq.0.or.
     &         arc(iarc)%link(2).eq.0) cycle
          narc_new = narc_new + 1
          if (narc_new.lt.iarc)
     &         arc(narc_new) = arc(iarc)
        end do
        contr%narc = narc_new
      case default
        call quit(1,'strip_contr','undefined mode')
      end select

      return
      end

max_stars_repo_name: ak-ustutt/GeCCo-public

max_issues_repo_name: ak-ustutt/GeCCo-public

max_forks_repo_name: ak-ustutt/GeCCo-public


text:
C23456789012345678901234567890123456789012345678901234567890123456789012
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                      C
C SUBROUTINE NODINP                                                    C
C                                                                      C
C                                                                      C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
      SUBROUTINE NODINP(NUMNP,ID,IBC,MXDOFDIM,MCRD,COORD,NDOFN,NEQ,IIN,
     1                  IOUT)
C
      IMPLICIT REAL*8 (A-H,O-Z)
C
      DIMENSION ID(MXDOFDIM,NUMNP),COORD(MCRD,NUMNP),NDOFN(NUMNP),
     1          IBC(MXDOFDIM,NUMNP)
C
      CALL CLEAR(COORD,MCRD,NUMNP)
      CALL CLEARIM(ID,MXDOFDIM,NUMNP)
      CALL CLEARIM(IBC,MXDOFDIM,NUMNP)
C
C     READ AND GENERATE NODAL POINT DATA
C
      WRITE(IOUT,1000)
      WRITE(IOUT,1005)
C
      DO I=1,NUMNP
        READ(IIN,     *) N,NDOFN(N),(ID(II,I),II=1,NDOFN(N)),
     1  (COORD(JJ,I),JJ=1,MCRD)
cc        WRITE(IOUT,1010)N,NDOFN(N),(IBC(II,I),II=1, MXDOFDIM),
cc     1  (COORD(JJ,I),JJ=1,MCRD)
      END DO
C
      ICOUNT=1
      DO K1=1,NUMNP
        NDOF=NDOFN(K1)
        DO K2=1,NDOFN(K1)
           IF(ID(K2,K1).EQ.0) THEN
              ID(K2,K1)=ICOUNT
              ICOUNT=ICOUNT+1
           ELSE
              ID(K2,K1)=0
           END IF
        END DO
      END DO
      NEQ=ICOUNT-1
C
C     WRITE EQUATION NUMBERS
C
      WRITE(IOUT,1020)
      WRITE(IOUT,1025)
      WRITE(IOUT,1030) (N,(ID(I,N),I=1,MXDOFDIM),N=1,NUMNP)
C
 1000 FORMAT(//,6X,' N O D A L  D A T A',//)
 1005 FORMAT('ID',2X,'NDOF',2X,'BC-X',2X,'BC-Y',2X,'BC-R',
     1        6X,'COORD-X',6X,'COORD-Y',/)
 1010 FORMAT(I3,X,I2,4X,I2,4X,I2,4X,I2,6X,F5.2,6X,F5.2)
 1020 FORMAT(/,'E Q U A T I O N  N U M B E R S',/)
 1025 FORMAT(/,'NODE',6X,'DEGREES OF FREEDOM',/,
     1      'NUMBER',8X,'U',9X,'V',9X,'R'/)
 1030 FORMAT(I5,5X,I5,5X,I5,5X,I5)
C
      RETURN
C
      END
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                      C
C SUBROUTINE MATINP                                                    C
C                                                                      C
C                                                                      C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
      SUBROUTINE MATINP(NUMAT,NMPR,AMATE,NMATP,IIN,IOUT)
C
      IMPLICIT REAL*8 (A-H,O-Z)
C
      DIMENSION AMATE(NMPR,NUMAT),NMATP(NUMAT)
C
      CALL CLEAR(AMATE,NMPR,NUMAT)
C
C     READS AND GENERATES MATERIAL DATA
C
      WRITE (IOUT,1000)
      DO I=1,NUMAT
        READ  (IIN, *) N,NMATP(N),(AMATE(JJ,N),JJ=1,NMATP(N))
        WRITE (IOUT,3020) N,NMATP(N),(AMATE(JJ,N),JJ=1,NMATP(N))
      END DO
C
 1000 FORMAT(//,4X,'M A T E R I A L  D A T A',/)
 3020 FORMAT(I5,3X,I5,12X,6(5X,F13.3))
C
      RETURN
C
      END
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                      C
C     SUBROUTINE ELEINP                                                C
C                                                                      C
C                                                                      C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
      SUBROUTINE ELEINP(NUMNP,NUMEL,NNE,IELT,NDOFEL,NMNE,MATP,IELCON,
     1                  IIN,IOUT)
C
      IMPLICIT REAL*8 (A-H,O-Z)
C
      DIMENSION NNE(NUMEL),IELT(NUMEL),MATP(NUMEL),IELCON(NMNE,NUMEL),
     1          NDOFEL(NUMEL)
C
C     Creates IELCON()
C
      WRITE(IOUT,1000)
      DO I=1,NUMEL
        READ(IIN,*) M,IELT(M),NDOFEL(M),MATP(M),NNE(M),
     1  (IELCON(J,M),J=1,NNE(M))
        WRITE(IOUT,1010) M,IELT(M),NDOFEL(M),MATP(M),NNE(M),
     1  (IELCON(J,M),J=1,NNE(M))
      END DO
C
 1000 FORMAT(///,8X,'E L E M E N T  I N F O R M A T I O N',//)
 1010 FORMAT(5(2X,I2),X,8(2X,I3))
C
      RETURN
C
      END
C

ext: for


text:
c
c     file mud3sa.f
c
c  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
c  .                                                             .
c  .                  copyright (c) 1999 by UCAR                 .
c  .                                                             .
c  .       University Corporation for Atmospheric Research       .
c  .                                                             .
c  .                      all rights reserved                    .
c  .                                                             .
c  .                                                             .
c  .                      MUDPACK  version 5.0                   .
c  .                                                             .
c  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
c
c ... author and specialist
c
c          John C. Adams (National Center for Atmospheric Research)
c          email: johnad@ucar.edu, phone: 303-497-1213

c ... For MUDPACK 5.0 information, visit the website:
c     (http://www.scd.ucar.edu/css/software/mudpack)
c
c ... purpose (see mud3sa.d for details)
c
c     mud3sa attempts to produce a second order finite difference
c     approximation to the three dimensional nonseparable selfadjoint
c     partial differential equation of the form:
c
c
c     (sigx(x,y,z)*px)x + (sigy(x,y,z)*py)y + (sigz(x,y,z)*pz)z
c
c     xlmbda(x,y,z)*p(x,y,z) = r(x,y,z)
c
c
c ... documentation and test files
c
c     see the documentation file "mud3sa.d" for a complete discussion
c     of how to use subroutine mud3sa.  file "tmud3sa.f" is a test/driver
c     sample program illustrating use of mud3sa
c
c ... required MUDPACK files
c
c     mudcom.f, mud3ln.f, mud3pn.f
c
      subroutine mud3sa(iparm,fparm,work,sgx,sgy,sgz,xlmbda,bndyc,rhs,
     +                  phi,mgopt,ierror)
      implicit none
      integer iparm(23),mgopt(4),ierror
      real work(*),phi(*),rhs(*),fparm(8)
      integer intl,nxa,nxb,nyc,nyd,nze,nzf,ixp,jyq,kzr,iex,jey,kez,
     +nfx,nfy,nfz,iguess,maxcy,method,meth2,nwork,lwork,itero,
     +kcycle,iprer,ipost,intpol
      common/imud3/intl,nxa,nxb,nyc,nyd,nze,nzf,ixp,jyq,kzr,iex,jey,kez,
     +nfx,nfy,nfz,iguess,maxcy,method,meth2,nwork,lwork,itero,
     +kcycle,iprer,ipost,intpol
      real xa,xb,yc,yd,ze,zf,tolmax,relmax
      common/fmud3/xa,xb,yc,yd,ze,zf,tolmax,relmax
      integer kpbgn,kcbgn,ktxbgn,ktybgn,ktzbgn,nxk,nyk,nzk,ngrid,
     +        klevel,kcur,kps
      common/mud3c/kpbgn(50),kcbgn(50),ktxbgn(50),ktybgn(50),
     +ktzbgn(50),nxk(50),nyk(50),nzk(50),ngrid,klevel,kcur,kps
      integer sxy,sxz,syz,kpspace,kcspace,int
      integer m,lxy,lxz,lyz,isx,jsy,ksz,iisx,jjsy,kksz
      integer nx,ny,nz,itx,ity,itz,k,kb,kkb,kk,ip,ic
      real sgx,sgy,sgz,xlmbda
      external sgx,sgy,sgz,xlmbda,bndyc
      data int / 0 /
      save int
      ierror = 1
      intl = iparm(1)    ! set and check intl on ALL calls
      if (intl*(intl-1).ne.0) return
      if (int.eq.0) then
	int = 1
	if (intl.ne.0) return  ! very first call is not intl=0
      end if
      ierror = 0
c
c     set input parameters from iparm,fparm internally
c
      intl = iparm(1)
      nxa = iparm(2)
      nxb = iparm(3)
      nyc = iparm(4)
      nyd = iparm(5)
      nze = iparm(6)
      nzf = iparm(7)
c
c     set grid size params
c
      ixp = iparm(8)
      jyq = iparm(9)
      kzr = iparm(10)
      iex = iparm(11)
      jey = iparm(12)
      kez = iparm(13)
c
c     set number of subgrids for mg cycling
c
      ngrid = max0(iex,jey,kez)
      nfx = iparm(14)
      nfy = iparm(15)
      nfz = iparm(16)

      iguess = iparm(17)
      maxcy = iparm(18)
      method = iparm(19)
      meth2 = iparm(20)
      nwork = iparm(21)
c
c     set floating point params
c
      xa = fparm(1)
      xb = fparm(2)
      yc = fparm(3)
      yd = fparm(4)
      ze = fparm(5)
      zf = fparm(6)
      tolmax = fparm(7)
c
c     set multigrid option parameters
c
      kcycle = mgopt(1)
      if (kcycle .eq. 0) then
c
c     use default settings
c
	kcycle = 2
	iprer = 2
	ipost = 1
	intpol = 3
      else
	iprer = mgopt(2)
	ipost = mgopt(3)
	intpol = mgopt(4)
      end if
      if (intl .eq. 0) then  ! intialization call
c
c     check input arguments
c
	ierror = 2   ! check boundary condition flags
	if (max0(nxa,nxb,nyc,nyd,nze,nzf).gt.2) return
	if (min0(nxa,nxb,nyc,nyd,nze,nzf).lt.0) return
	if (nxa.eq.0.and.nxb.ne.0) return
	if (nxa.ne.0.and.nxb.eq.0) return
	if (nyc.eq.0.and.nyd.ne.0) return
	if (nyc.ne.0.and.nyd.eq.0) return
	if (nze.eq.0.and.nzf.ne.0) return
	if (nze.ne.0.and.nzf.eq.0) return
	ierror = 3   ! check grid sizes
	if (ixp.lt.2) return
	if (jyq.lt.2) return
	if (kzr.lt.2) return
	ierror = 4
	ngrid = max0(iex,jey,kez)
	if (iex.lt.1) return
	if (jey.lt.1) return
	if (kez.lt.1) return
	if (ngrid.gt.50) return
	ierror = 5
	if (nfx.ne.ixp*2**(iex-1)+1) return
	if (nfy.ne.jyq*2**(jey-1)+1) return
	if (nfz.ne.kzr*2**(kez-1)+1) return
	ierror = 6
	if (iguess*(iguess-1).ne.0) return
	ierror = 7
	if (maxcy.lt.1) return
	ierror = 8
	if (method.lt.0 .or. method.gt.10) return
	if (meth2.lt.0 .or. meth2.gt.3) return
	ierror = 9
c
c     compute required work space length
c
	m = method
	isx = 0
	if ((m-1)*(m-4)*(m-5)*(m-7).eq.0) then
	  isx = 3
	  if (nxa.eq.0) then
	    isx = 5
	  end if
	end if
	jsy = 0
	if ((m-2)*(m-4)*(m-6)*(m-7).eq.0) then
	  jsy = 3
	  if (nyc.eq.0) then
	    jsy = 5
	  end if
	end if
	ksz = 0
	if ((m-3)*(m-5)*(m-6)*(m-7).eq.0) then
	  ksz = 3
	  if (nze.eq.0) then
	    ksz = 5
	  end if
	end if
c
c     set scales for planar relaxation
c
	iisx = 0
	jjsy = 0
	kksz = 0
	lxy = 0
	lxz = 0
	lyz = 0
	nx = nfx
	ny = nfy
	nz = nfz
	if (method.eq.8) then
	  lxy = 1
	  if (meth2.eq.1.or.meth2.eq.3) then
	    iisx = 3
	    if (nxa.eq.0) iisx = 5
	  end if
	  if (meth2.eq.2.or.meth2.eq.3) then
	    jjsy = 3
	    if (nyc.eq.0) jjsy = 5
	  end if
	end if
	if (method.eq.9) then
	  lxz = 1
	  if (meth2.eq.1.or.meth2.eq.3) then
	    iisx = 3
	    if (nxa.eq.0) iisx = 5
	  end if
	  if (meth2.eq.2.or.meth2.eq.3) then
	    kksz = 3
	    if (nze.eq.0) kksz = 5
	    end if
	end if
	if (method.eq.10) then
	lyz = 1
	  if (meth2.eq.1.or.meth2.eq.3) then
	    jjsy = 3
	    if (nyc.eq.0) jjsy = 5
	  end if
	  if (meth2.eq.2.or.meth2.eq.3) then
	    kksz = 3
	    if (nze.eq.0) kksz = 5
	  end if
	end if
	sxy = 0
	sxz = 0
	syz = 0
c
c     set subgrid sizes
c
	do k=1,ngrid
	  nxk(k) = ixp*2**(max0(k+iex-ngrid,1)-1)+1
	  nyk(k) = jyq*2**(max0(k+jey-ngrid,1)-1)+1
	  nzk(k) = kzr*2**(max0(k+kez-ngrid,1)-1)+1
	end do
	kps = 1
	do kb=1,ngrid
	  k = ngrid-kb+1
	  kpspace = 0
	  kcspace = 0
	  if (method.gt.7) then
c
c     set spacers for planar relaxation
c
	    do kkb=2,k
	      kk = k-kkb+1
	      nx = nxk(kk)
	      ny = nyk(kk)
	      nz = nzk(kk)
	      if (method.eq.8) nz = nzk(k)
	      if (method.eq.9) ny = nyk(k)
	      if (method.eq.10) nx = nxk(k)
	      kpspace = kpspace + (nx+2)*(ny+2)*(nz+2)
	      kcspace = kcspace + 8*nx*ny*nz
	    end do
	  end if
	  nx = nxk(k)
	  ny = nyk(k)
	  nz = nzk(k)
c
c     set pointers
c
	  kpbgn(k) = kps
	  kcbgn(k) = kpbgn(k)+(nx+2)*(ny+2)*(nz+2)+kpspace
	  ktxbgn(k) = kcbgn(k) + 8*nx*ny*nz + kcspace
	  ktybgn(k) = ktxbgn(k) + isx*nx*ny*nz
	  ktzbgn(k) = ktybgn(k) + jsy*nx*ny*nz
	  kps = ktzbgn(k) + ksz*nx*ny*nz
c
c     sum space in case planar relaxation in xy or xz or yz
c
	  sxy = sxy + (6+iisx+jjsy)*nx*ny + (nx+2)*(ny+2)
	  sxz = sxz + (6+iisx+kksz)*nx*nz + (nx+2)*(nz+2)
	  syz = syz + (6+jjsy+kksz)*ny*nz + (ny+2)*(nz+2)
	end do
c
c     set and check minimal work space
c
	nx = nxk(ngrid)
	ny = nyk(ngrid)
	nz = nzk(ngrid)
	iparm(22)=kps+max0((nx+2)*(ny+2)*(nz+2),lxy*sxy,lxz*sxz,lyz*syz)
	lwork = iparm(22)
	if (lwork .gt. nwork) return
	ierror = 10   ! check solution region
	if (xb.le.xa .or. yd.le.yc .or. zf.le.ze) return
	ierror = 11
	if (tolmax .lt. 0.0) return
	ierror = 12   ! multigrid parameters
	if (kcycle.lt.0) return
	if (min0(iprer,ipost).lt.1) return
	if ((intpol-1)*(intpol-3).ne.0) return
	if (max0(kcycle,iprer,ipost).gt.2) then
	  ierror = -5   ! inefficient multigrid cycling
	end if
	if (ierror .gt. 0) ierror = 0   ! no fatal errors
c
c     discretize pde at each grid level
c
	do kb=1,ngrid
	  k = ngrid-kb+1
	  nx = nxk(k)
	  ny = nyk(k)
	  nz = nzk(k)
	  ip = kpbgn(k)
	  ic = kcbgn(k)
	  itx = ktxbgn(k)
	  ity = ktybgn(k)
	  itz = ktzbgn(k)
	  call dismd3sa(nx,ny,nz,work(ic),work(itx),work(ity),
     +    work(itz),bndyc,sgx,sgy,sgz,xlmbda,work,ierror)
	  if (method.gt.7) then
c
c      discretize for planar coarsening
c
	    do kkb=2,k
	      kk = k-kkb+1
	      ip = ip+(nx+2)*(ny+2)*(nz+2)
	      ic = ic + 8*nx*ny*nz
	      nx = nxk(kk)
	      ny = nyk(kk)
	      nz = nzk(kk)
	      if (method.eq.8) nz = nzk(k)
	      if (method.eq.9) ny = nyk(k)
	      if (method.eq.10) nx = nxk(k)
	      call dismd3sa(nx,ny,nz,work(ic),work(itx),work(ity),
     +        work(itz),bndyc,sgx,sgy,sgz,xlmbda,work,ierror)
	    end do
	  end if
	end do
	return
      end if   ! end of intl=0 initialization call block
      nx = nfx
      ny = nfy
      nz = nfz
      call mud3sa1(nx,ny,nz,rhs,phi,sgx,sgy,sgz,xlmbda,bndyc,work)
      iparm(23) = itero
      if (ierror .le. 0) then
	if (tolmax.gt.0.0) then
c
c     set final computed maximum relative difference
c
	  fparm(8) = relmax
	  if (relmax.gt.tolmax .and. ierror.eq.0) then
c
c     flag convergence failure
c
	    ierror = -1
	  end if
	end if
      end if
      return
      end

      subroutine mud3sa1(nx,ny,nz,rhsf,phif,sgx,sgy,sgz,xlmbda,bndyc,wk)
      implicit none
      integer intl,nxa,nxb,nyc,nyd,nze,nzf,ixp,jyq,kzr,iex,jey,kez,
     +nfx,nfy,nfz,iguess,maxcy,method,meth2,nwork,lwork,itero,
     +kcycle,iprer,ipost,intpol
      real xa,xb,yc,yd,ze,zf,tolmax,relmax
      integer kpbgn,kcbgn,ktxbgn,ktybgn,ktzbgn,nxk,nyk,nzk,ngrid,
     +        klevel,kcur,kps
      integer nx,ny,nz,ip,ic,ir,icc,irc,ncx,ncy,ncz,ipc
      integer i,j,k,kb,jk,kk,kkb,ijk,iter
      real phif(nx,ny,nz),rhsf(nx,ny,nz),wk(*),phmax
      common/imud3/intl,nxa,nxb,nyc,nyd,nze,nzf,ixp,jyq,kzr,iex,jey,kez,
     +nfx,nfy,nfz,iguess,maxcy,method,meth2,nwork,lwork,itero,
     +kcycle,iprer,ipost,intpol
      common/fmud3/xa,xb,yc,yd,ze,zf,tolmax,relmax
      common/mud3c/kpbgn(50),kcbgn(50),ktxbgn(50),ktybgn(50),
     +ktzbgn(50),nxk(50),nyk(50),nzk(50),ngrid,klevel,kcur,kps
      real sgx,sgy,sgz,xlmbda
      external sgx,sgy,sgz,xlmbda,bndyc
      nx = nxk(ngrid)
      ny = nyk(ngrid)
      nz = nzk(ngrid)
      ip = kpbgn(ngrid)
      ic = kcbgn(ngrid)
      ir = ic+7*nx*ny*nz
c
c     set phif,rhsf in wk
c
      call swk3(nx,ny,nz,phif,rhsf,wk(ip),wk(ir))
      if (iguess.eq.0) then
c
c     no initial guess at finest grid level!
c
	do kb=2,ngrid
	  k = ngrid-kb+1
	  nx = nxk(k+1)
	  ny = nyk(k+1)
	  nz = nzk(k+1)
	  ip = kpbgn(k+1)
	  ic = kcbgn(k+1)
	  ir = kcbgn(k+1)+7*nx*ny*nz
	  ncx = nxk(k)
	  ncy = nyk(k)
	  ncz = nzk(k)
	  ipc = kpbgn(k)
	  icc = kcbgn(k)
	  irc = icc+7*ncx*ncy*ncz
c
c     transfer down to all grid levels
c
	  call trsfc3(nx,ny,nz,wk(ip),wk(ir),ncx,ncy,ncz,
     +                wk(ipc),wk(irc))
	  if (method.gt.7) then
c
c     transfer down for planar coarsening
c
	    do kkb=1,k-1
	      kk = k-kkb+1
	      ipc = ip + (nx+2)*(ny+2)*(nz+2)
	      icc = ic + 8*nx*ny*nz
	      ncx = nxk(kk)
	      ncy = nyk(kk)
	      ncz = nzk(kk)
	      if (method.eq.8) then
		ncz = nzk(k+1)
	      else if (method.eq.9) then
		ncy = nyk(k+1)
	      else
		ncx = nxk(k+1)
	      end if
	      irc = icc+7*ncx*ncy*ncz
	      call trsfc3(nx,ny,nz,wk(ip),wk(ir),ncx,ncy,ncz,
     +                    wk(ipc),wk(irc))
	      nx = ncx
	      ny = ncy
	      nz = ncz
	      ip = ipc
	      ir = irc
	      ic = icc
	    end do
	  end if
	end do
c
c     adjust right hand side at all grid levels in case
c     rhs or specified b.c. in phi or gbdy changed
c
	do kb=1,ngrid
	  k = ngrid-kb+1
	  nx = nxk(k)
	  ny = nyk(k)
	  nz = nzk(k)
	  ip = kpbgn(k)
	  ic = kcbgn(k)
	  call adjmd3sa(nx,ny,nz,wk(ip),wk(ic),bndyc,sgx,sgy,sgz,xlmbda)
c
c     adjust for planar grid levels if necessary
c
	  if (method.gt.7) then
	    do kkb=2,k
	      kk = k-kkb+1
	      ip = ip+(nx+2)*(ny+2)*(nz+2)
	      ic = ic + 8*nx*ny*nz
	      nx = nxk(kk)
	      ny = nyk(kk)
	      nz = nzk(kk)
	      if (method.eq.8) then
		nz = nzk(k)
	      else if (method.eq.9) then
		ny = nyk(k)
	      else
		nx = nxk(k)
	      end if
	      call adjmd3sa(nx,ny,nz,wk(ip),wk(ic),bndyc,sgx,sgy,sgz,
     +                      xlmbda)
	    end do
	  end if
	end do
c
c     execute one full multigrid cycle
c
	do k=1,ngrid-1
	  kcur = k
	  call kcymd3sa(wk)
	  nx = nxk(k+1)
	  ny = nyk(k+1)
	  nz = nzk(k+1)
	  ip = kpbgn(k+1)
	  ipc = kpbgn(k)
	  ncx = nxk(k)
	  ncy = nyk(k)
	  ncz = nzk(k)
c
c     lift or prolong approximation from k to k+1
c
	  call prolon3(ncx,ncy,ncz,wk(ipc),nx,ny,nz,wk(ip),
     +                 nxa,nxb,nyc,nyd,nze,nzf,intpol)
	end do

      else
c
c     adjust rhs at finest grid level only
c
	nx = nxk(ngrid)
	ny = nyk(ngrid)
	nz = nzk(ngrid)
	ip = kpbgn(ngrid)
	ic = kcbgn(ngrid)
	call adjmd3sa(nx,ny,nz,wk(ip),wk(ic),bndyc,sgx,sgy,sgz,xlmbda)
      end if
c
c     execute maxcy more multigrid k cycles from finest level
c
      kcur = ngrid
      do iter=1,maxcy
	itero = iter
	call kcymd3sa(wk)
	if (tolmax.gt.0.0) then
c
c      error control
c
	  relmax = 0.0
	  phmax = 0.0
	  do k=1,nfz
	    kk = k*(nfx+2)*(nfy+2)
	    do j=1,nfy
	      jk = kk+j*(nfx+2)
	      do i=1,nfx
		ijk = jk+i+1
		phmax = amax1(phmax,abs(wk(ijk)))
		relmax = amax1(relmax,abs(wk(ijk)-phif(i,j,k)))
		phif(i,j,k) = wk(ijk)
	      end do
	    end do
	  end do
c
c     set maximum relative difference and check for convergence
c
	  if (phmax.gt.0.0) relmax = relmax/phmax
	  if (relmax.le.tolmax) return
	end if
      end do
c
c     set final iterate in phif
c
      do k=1,nfz
	kk = k*(nfx+2)*(nfy+2)
	do j=1,nfy
	  jk = kk+j*(nfx+2)
	  do i=1,nfx
	    ijk = jk+i+1
	    phif(i,j,k) = wk(ijk)
	  end do
	end do
      end do
      return
      end

      subroutine kcymd3sa(wk)
c
c     perform multigrid k-cycle at kcur level
c     kcycle = 1 corresponds to v cycles
c     kcycle = 2 corresponds to w cycles
c
      implicit none
      real wk(*)
      integer intl,nxa,nxb,nyc,nyd,nze,nzf,ixp,jyq,kzr,iex,jey,kez,
     +nfx,nfy,nfz,iguess,maxcy,method,meth2,nwork,lwork,itero,
     +kcycle,iprer,ipost,intpol
      common/imud3/intl,nxa,nxb,nyc,nyd,nze,nzf,ixp,jyq,kzr,iex,jey,kez,
     +nfx,nfy,nfz,iguess,maxcy,method,meth2,nwork,lwork,itero,
     +kcycle,iprer,ipost,intpol
      integer kpbgn,kcbgn,ktxbgn,ktybgn,ktzbgn,nxk,nyk,nzk,ngrid,
     +        klevel,kcur,kps
      common/mud3c/kpbgn(50),kcbgn(50),ktxbgn(50),ktybgn(50),
     +ktzbgn(50),nxk(50),nyk(50),nzk(50),ngrid,klevel,kcur,kps
      integer nx,ny,nz,ncx,ncy,ncz
      integer kount(50),ip,ic,ipc,irc,nrel,l
      klevel = kcur
c
c     pre-relax at current finest grid level
c
      do l = 1,iprer
	call relmd3sa(wk)
      end do
c
c     if at coarsest level post-relax
c
      if (kcur .eq. 1) go to 2
c
c     restrict residual to kcur-1
c
      nx = nxk(klevel)
      ny = nyk(klevel)
      nz = nzk(klevel)
      ip = kpbgn(klevel)
      ic = kcbgn(klevel)
      ipc = kpbgn(klevel-1)
      ncx = nxk(klevel-1)
      ncy = nyk(klevel-1)
      ncz = nzk(klevel-1)
      irc = kcbgn(klevel-1)+7*ncx*ncy*ncz
c
c     use full weighting with residual restriction
c
      call resmd3sa(nx,ny,nz,wk(ip),wk(ic),ncx,ncy,ncz,wk(ipc),wk(irc),
     +              wk(kps))
c
c     set counter for grid levels to zero
c
      do l = 1,kcur
	kount(l) = 0
      end do
c
c    set new level and continue k-cycling
c
      klevel = kcur-1
      nrel = iprer
c
c     kcycle control point
c
    1 continue
c
c     post-relax when kcur revisited
c
      if (klevel .eq. kcur) go to 2
c
c     count "hit" at current level
c
      kount(klevel) = kount(klevel)+1
c
c     relax at current level
c
      do l = 1,nrel
	call relmd3sa(wk)
      end do
      if (kount(klevel) .eq. kcycle+1) then
c
c     kcycle(iprer,ipost) complete at klevel
c     inject correction to finer grid
c
	nx = nxk(klevel+1)
	ny = nyk(klevel+1)
	nz = nzk(klevel+1)
	ip = kpbgn(klevel+1)
	ncx = nxk(klevel)
	ncy = nyk(klevel)
	ncz = nzk(klevel)
	ipc = kpbgn(klevel)
	call cor3(nx,ny,nz,wk(ip),ncx,ncy,ncz,wk(ipc),
     +            nxa,nxb,nyc,nyd,nze,nzf,intpol,wk(kps))
c
c     reset counter to zero at klevel
c
	kount(klevel) = 0
c
c     ascend to next higher level and set to post-relax there
c
	klevel = klevel+1
	nrel = ipost
	go to 1
      else
c
c     kcycle not complete so descend unless at coarsest
c
	if (klevel .gt. 1) then
	  nx = nxk(klevel)
	  ny = nyk(klevel)
	  nz = nzk(klevel)
	  ip = kpbgn(klevel)
	  ic = kcbgn(klevel)
	  ncx = nxk(klevel-1)
	  ncy = nyk(klevel-1)
	  ncz = nzk(klevel-1)
	  irc = kcbgn(klevel-1)+7*ncx*ncy*ncz
	  ipc = kpbgn(klevel-1)
	  call resmd3sa(nx,ny,nz,wk(ip),wk(ic),ncx,ncy,ncz,wk(ipc),
     +                  wk(irc),wk(kps))
c
c     pre-relax at next coarser level
c
	  klevel = klevel-1
	  nrel = iprer
	  go to 1
	else
c
c     post-relax at coarsest level (klevel=1)
c
	  do l = 1,ipost
	    call relmd3sa(wk)
	  end do
c
c     inject correction to grid level 2
c
	  ipc = kpbgn(1)
	  ncx = nxk(1)
	  ncy = nyk(1)
	  ncz = nzk(1)
	  ip = kpbgn(2)
	  nx = nxk(2)
	  ny = nyk(2)
	  nz = nzk(2)
	  call cor3(nx,ny,nz,wk(ip),ncx,ncy,ncz,wk(ipc),
     +              nxa,nxb,nyc,nyd,nze,nzf,intpol,wk(kps))
c
c     set to post-relax at level 2
c
	  nrel = ipost
	  klevel = 2
	  go to 1
	end if
      end if
    2 continue
c
c     post-relax at kcur level
c
      do l = 1,ipost
	call relmd3sa(wk)
      end do
      return
      end

      subroutine resmd3sa(nx,ny,nz,phi,cof,ncx,ncy,ncz,phic,rhsc,resf)
c
c     compute fully weighted residual restriction in rhsc
c
      implicit none
      integer nx,ny,nz,ncx,ncy,ncz
      integer intl,nxa,nxb,nyc,nyd,nze,nzf,ixp,jyq,kzr,iex,jey,kez,
     +nfx,nfy,nfz,iguess,maxcy,method,meth2,nwork,lwork,itero,
     +kcycle,iprer,ipost,intpol
      common/imud3/intl,nxa,nxb,nyc,nyd,nze,nzf,ixp,jyq,kzr,iex,jey,kez,
     +nfx,nfy,nfz,iguess,maxcy,method,meth2,nwork,lwork,itero,
     +kcycle,iprer,ipost,intpol
      real phi(0:nx+1,0:ny+1,0:nz+1),phic(0:ncx+1,0:ncy+1,0:ncz+1)
      real rhsc(ncx,ncy,ncz),resf(nx,ny,nz),cof(nx,ny,nz,8)
      integer ic,jc,kc,i,j,k
c
c     initialize phic to zero
c
      do kc=0,ncz+1
	do jc=0,ncy+1
	  do ic=0,ncx+1
	    phic(ic,jc,kc) = 0.0
	  end do
	end do
      end do
c
c     compute fine grid residual
c
!$OMP PARALLEL DO PRIVATE(i,j,k), SHARED(resf,cof,phi,nx,ny,nz)
      do k=1,nz
	do j=1,ny
	  do i=1,nx
	    resf(i,j,k) =  cof(i,j,k,8)-(
     +                     cof(i,j,k,1)*phi(i-1,j,k)+
     +                     cof(i,j,k,2)*phi(i+1,j,k)+
     +                     cof(i,j,k,3)*phi(i,j-1,k)+
     +                     cof(i,j,k,4)*phi(i,j+1,k)+
     +                     cof(i,j,k,5)*phi(i,j,k-1)+
     +                     cof(i,j,k,6)*phi(i,j,k+1)+
     +                     cof(i,j,k,7)*phi(i,j,k))
	  end do
	end do
      end do
c
c     restrict resf to interior coarse mesh in rhsc
c     using fully weighted residual restriction in 3-d
c
      call res3(nx,ny,nz,resf,ncx,ncy,ncz,rhsc,nxa,nxb,nyc,nyd,nze,nzf)
      return
      end

      subroutine dismd3sa(nx,ny,nz,cof,tx,ty,tz,bndyc,sgx,sgy,sgz,
     +                    xlmbda,wk,ier)
c
c     discretize the 3-d elliptic pde
c
      implicit none
      integer nx,ny,nz,ier
      real cof(nx,ny,nz,8)
      real tx(nx,ny,nz,*),ty(ny,nx,nz,*),tz(nz,nx,ny,*),wk(*)
      integer intl,nxa,nxb,nyc,nyd,nze,nzf,ixp,jyq,kzr,iex,jey,kez,
     +nfx,nfy,nfz,iguess,maxcy,method,meth2,nwork,lwork,itero,
     +kcycle,iprer,ipost,intpol
      real xa,xb,yc,yd,ze,zf,tolmax,relmax
      integer kpbgn,kcbgn,ktxbgn,ktybgn,ktzbgn,nxk,nyk,nzk,ngrid,
     +        klevel,kcur,kps
      common/imud3/intl,nxa,nxb,nyc,nyd,nze,nzf,ixp,jyq,kzr,iex,jey,kez,
     +nfx,nfy,nfz,iguess,maxcy,method,meth2,nwork,lwork,itero,
     +kcycle,iprer,ipost,intpol
      common/fmud3/xa,xb,yc,yd,ze,zf,tolmax,relmax
      common/mud3c/kpbgn(50),kcbgn(50),ktxbgn(50),ktybgn(50),
     +ktzbgn(50),nxk(50),nyk(50),nzk(50),ngrid,klevel,kcur,kps
      real dlx,dly,dlz,dlx2,dly2,dlz2,dlxx,dlyy,dlzz,alfmax
      real odlxx,odlyy,odlzz
      real sgxm,sgxp,sgym,sgyp,sgzm,sgzp,sigmin,xlm,xlmin,xlmax
      real alfa,gbdy,x,y,z,c1,c2,c3,c4,c5,c6,c7
      integer i,j,k,l,ist,ifn,jst,jfn,kst,kfn,kbdy
      integer nxny,nxnz,nynz,im1,jm1,km1
      real sgx,sgy,sgz,xlmbda
      external sgx,sgy,sgz,xlmbda,bndyc
c
c     set current grid increments
c
      dlx = (xb-xa)/(nx-1)
      dlx2 = dlx+dlx
      dlxx = dlx*dlx
      dly = (yd-yc)/(ny-1)
      dly2 = dly+dly
      dlyy = dly*dly
      dlz = (zf-ze)/(nz-1)
      dlz2 = dlz+dlz
      dlzz = dlz*dlz
      odlxx = 1.0/dlxx
      odlyy = 1.0/dlyy
      odlzz = 1.0/dlzz
c
c     set x,y,z subscript limits to bypass specified boundaries
c     when calling sgx,sgy,sgz,xlmbda or bndyc
c
      jst = 1
      jfn = ny
      ist = 1
      ifn = nx
      kst = 1
      kfn = nz
      if (nxa.eq.1) ist = 2
      if (nxb.eq.1) ifn = nx-1
      if (nyc.eq.1) jst = 2
      if (nyd.eq.1) jfn = ny-1
      if (nze.eq.1) kst = 2
      if (nzf.eq.1) kfn = nz-1
      sigmin = 1.0
      xlmin = 1.0
      xlmax = 0.0
      do k=kst,kfn
	z = ze+(k-1)*dlz
	do j=jst,jfn
	  y = yc+(j-1)*dly
	  do i=ist,ifn
	    x = xa+(i-1)*dlx
	    sgxm = sgx(x-0.5*dlx,y,z)
	    sgxp = sgx(x+0.5*dlx,y,z)
	    sgym = sgy(x,y-0.5*dly,z)
	    sgyp = sgy(x,y+0.5*dly,z)
	    sgzm = sgz(x,y,z-0.5*dlz)
	    sgzp = sgz(x,y,z+0.5*dlz)
	    xlm = xlmbda(x,y,z)
c
c     monitor possible nonfatal coefficient errors
c
	    sigmin = amin1(sigmin,sgx(x,y,z),sgy(x,y,z),sgz(x,y,z))
	    xlmin = amin1(xlmin,xlm)
	    xlmax = amax1(abs(xlm),xlmax)
	    c1 = sgxm*odlxx
	    c2 = sgxp*odlxx
	    c3 = sgym*odlyy
	    c4 = sgyp*odlyy
	    c5 = sgzm*odlzz
	    c6 = sgzp*odlzz
	    c7 = -xlm-(c1+c2+c3+c4+c5+c6)
	    cof(i,j,k,1) = c1
	    cof(i,j,k,2) = c2
	    cof(i,j,k,3) = c3
	    cof(i,j,k,4) = c4
	    cof(i,j,k,5) = c5
	    cof(i,j,k,6) = c6
	    cof(i,j,k,7) = c7
	  end do
	end do
      end do
c
c     flag non-ellipticity
c
      if (sigmin.le.0.0) then
	ier = -2
      end if
c
c     set warning flag if xlmbda is less then zero at some point
c
      if (xlmin.lt.0.0) then
	if (ier.eq.0) then
	  ier = -4
	end if
      end if
      alfmax = 0.0
c
c     adjust equation at mixed b.c.
c
      if (nxa.eq.2) then
	kbdy = 1
	x = xa
	i = 1
	do k=kst,kfn
	  z = ze+(k-1)*dlz
	  do j=jst,jfn
	    y = yc+(j-1)*dly
	    call bndyc(kbdy,y,z,alfa,gbdy)
	    alfmax = amax1(abs(alfa),alfmax)
	    c1 = cof(i,j,k,1)
	    cof(i,j,k,1) = 0.0
	    cof(i,j,k,2) = cof(i,j,k,2)+c1
	    cof(i,j,k,7) = cof(i,j,k,7)+dlx2*alfa*c1
	  end do
	end do
      end if
      if (nxb.eq.2) then
	kbdy = 2
	x = xb
	i = nx
	do k=kst,kfn
	  z = ze+(k-1)*dlz
	  do j=jst,jfn
	    y = yc+(j-1)*dly
	    call bndyc(kbdy,y,z,alfa,gbdy)
	    alfmax = amax1(abs(alfa),alfmax)
	    c2 = cof(i,j,k,2)
	    cof(i,j,k,1) = cof(i,j,k,1)+c2
	    cof(i,j,k,2) = 0.0
	    cof(i,j,k,7) = cof(i,j,k,7)-dlx2*alfa*c2
	  end do
	end do
      end if
      if (nyc.eq.2) then
	kbdy = 3
	y = yc
	j = 1
	do k=kst,kfn
	  z = ze+(k-1)*dlz
	  do i=ist,ifn
	    x = xa+(i-1)*dlx
	    call bndyc(kbdy,x,z,alfa,gbdy)
	    alfmax = amax1(abs(alfa),alfmax)
	    c3 = cof(i,j,k,3)
	    cof(i,j,k,3) = 0.0
	    cof(i,j,k,4) = cof(i,j,k,4)+c3
	    cof(i,j,k,7) = cof(i,j,k,7)+dly2*alfa*c3
	  end do
	end do
      end if
      if (nyd.eq.2) then
	kbdy = 4
	y = yd
	j = ny
	do k=kst,kfn
	z = ze+(k-1)*dlz
	do i=ist,ifn
	    x = xa+(i-1)*dlx
	    call bndyc(kbdy,x,z,alfa,gbdy)
	    alfmax = amax1(abs(alfa),alfmax)
	    c4 = cof(i,j,k,4)
	    cof(i,j,k,3) = cof(i,j,k,3)+c4
	    cof(i,j,k,4) = 0.0
	    cof(i,j,k,7) = cof(i,j,k,7)-dly2*c4*alfa
	  end do
	end do
      end if
      if (nze.eq.2) then
	kbdy = 5
	z = ze
	k = 1
	do j=jst,jfn
	  y = yc+(j-1)*dly
	  do i=ist,ifn
	    x = xa+(i-1)*dlx
	    call bndyc(kbdy,x,y,alfa,gbdy)
	    alfmax = amax1(abs(alfa),alfmax)
	    c5 = cof(i,j,k,5)
	    cof(i,j,k,5) = 0.0
	    cof(i,j,k,6) = cof(i,j,k,6)+c5
	    cof(i,j,k,7) = cof(i,j,k,7)+dlz2*c5*alfa
	  end do
	end do
      end if
      if (nzf.eq.2) then
	kbdy = 6
	z = zf
	k = nz
	do j=jst,jfn
	  y = yc+(j-1)*dly
	  do i=ist,ifn
	    x = xa+(i-1)*dlx
	    call bndyc(kbdy,x,y,alfa,gbdy)
	    alfmax = amax1(abs(alfa),alfmax)
	    c6 = cof(i,j,k,6)
	    cof(i,j,k,5) = cof(i,j,k,5)+c6
	    cof(i,j,k,6) = 0.0
	    cof(i,j,k,7) = cof(i,j,k,7)-dlz2*c6*alfa
	  end do
	end do
      end if
c
c     flag continuous singular elliptic pde if detected
c
      if (ier .ne. -4) then
	if (xlmax.eq.0.0.and.alfmax.eq.0.0) then
	  if (nxa.eq.0.or.(nxa.eq.2.and.nxb.eq.2)) then
	    if (nyc.eq.0.or.(nyc.eq.2.and.nyd.eq.2)) then
	      if (nze.eq.0.or.(nze.eq.2.and.nzf.eq.2)) then
		ier = -3
	      end if
	    end if
	  end if
	end if
      end if
c
c     reset cof for specified b.c.
c
      if (nxa.eq.1) then
	i = 1
	do j=1,ny
	  do k=1,nz
	    do l=1,7
	      cof(i,j,k,l) = 0.0
	    end do
	    cof(i,j,k,7) = 1.0
	  end do
	end do
      end if
      if (nxb.eq.1) then
	i = nx
	do k=1,nz
	  do j=1,ny
	    do l=1,7
	      cof(i,j,k,l) = 0.0
	    end do
	    cof(i,j,k,7) = 1.0
	  end do
	end do
      end if
      if (nyc.eq.1) then
	j = 1
	do k=1,nz
	  do i=1,nx
	    do l=1,7
	      cof(i,j,k,l) = 0.0
	    end do
	    cof(i,j,k,7) = 1.0
	  end do
	end do
      end if
      if (nyd.eq.1) then
	j = ny
	do i=1,nx
	  do k=1,nz
	    do l=1,7
	      cof(i,j,k,l) = 0.0
	    end do
	    cof(i,j,k,7) = 1.0
	  end do
	end do
      end if
      if (nze.eq.1) then
	k = 1
	do j=1,ny
	  do i=1,nx
	    do l=1,7
	      cof(i,j,k,l) = 0.0
	    end do
	    cof(i,j,k,7) = 1.0
	  end do
	end do
      end if
      if (nzf.eq.1) then
	k = nz
	do j=1,ny
	  do i=1,nx
	    do l=1,7
	      cof(i,j,k,l) = 0.0
	    end do
	    cof(i,j,k,7) = 1.0
	  end do
	end do
      end if
      if (method*(method-8)*(method-9)*(method-10) .eq. 0) return
c
c     set,factor tridiagonal matrices for line relaxation
c
      if ((method-1)*(method-4)*(method-5)*(method-7).eq.0) then
c
c     line relaxation in x used
c
	if (nxa.ne.0) then
c
c     set non-periodic tridiagonal matrices in tx and factor
c
	  do i=1,nx
	    im1 = max0(i-1,1)
	    do k=1,nz
	      do j=1,ny
		tx(im1,j,k,1) = cof(i,j,k,1)
		tx(i,j,k,2) = cof(i,j,k,7)
		tx(i,j,k,3) = cof(i,j,k,2)
	      end do
	    end do
	  end do
	  nynz = ny*nz
	  call factri(nynz,nx,tx(1,1,1,1),tx(1,1,1,2),tx(1,1,1,3))
	else
	  if (nx .gt. 3) then
c
c     set "periodic" tridiagonal matrices in tx and factor when nx > 3
c
	    do k=1,nz
	      do j=1,ny
		do i=1,nx-1
		  tx(i,j,k,1) = cof(i,j,k,1)
		  tx(i,j,k,2) = cof(i,j,k,7)
		  tx(i,j,k,3) = cof(i,j,k,2)
		end do
	      end do
	    end do
	    nynz = ny*nz
	    call factrp(nynz,nx,tx,tx(1,1,1,2),tx(1,1,1,3),tx(1,1,1,4),
     +                  tx(1,1,1,5),wk(kps))
	  end if
	end if
      end if
      if ((method-2)*(method-4)*(method-6)*(method-7).eq.0) then
c
c     line relaxation in y used
c
	if (nyc.ne.0) then
c
c     set non-periodic tridiagonal matrices and factor
c
	  do j=1,ny
	    jm1 = max0(j-1,1)
	    do k=1,nz
	      do i=1,nx
		ty(jm1,i,k,1) = cof(i,j,k,3)
		ty(j,i,k,2) = cof(i,j,k,7)
		ty(j,i,k,3) = cof(i,j,k,4)
	      end do
	    end do
	  end do
	  nxnz = nx*nz
	  call factri(nxnz,ny,ty(1,1,1,1),ty(1,1,1,2),ty(1,1,1,3))
	else
	  if (ny .gt. 3) then
c
c     set and factor periodic "tridiagonal" matrices when ny > 3
c
	    do k=1,nz
	      do i=1,nx
		do j=1,ny-1
		  ty(j,i,k,1) = cof(i,j,k,3)
		  ty(j,i,k,2) = cof(i,j,k,7)
		  ty(j,i,k,3) = cof(i,j,k,4)
		end do
	      end do
	    end do
	    nxnz = nx*nz
	    call factrp(nxnz,ny,ty,ty(1,1,1,2),ty(1,1,1,3),ty(1,1,1,4),
     +                  ty(1,1,1,5),wk(kps))
	  end if
	end if
      end if
      if ((method-3)*(method-5)*(method-6)*(method-7).eq.0) then
c
c     line relaxation in z used
c
	if (nze.ne.0) then
c
c     set and factor non-periodic tridiagonal matrices
c
	  do k=1,nz
	    km1 = max0(k-1,1)
	    do j=1,ny
	      do i=1,nx
		tz(km1,i,j,1) = cof(i,j,k,5)
		tz(k,i,j,2) = cof(i,j,k,7)
		tz(k,i,j,3) = cof(i,j,k,6)
	      end do
	    end do
	  end do
	  nxny = nx*ny
	  call factri(nxny,nz,tz(1,1,1,1),tz(1,1,1,2),tz(1,1,1,3))
	else
	  if (nz .gt. 3) then
c
c     set and factor periodic "tridiagonal matrices when nz > 3
c
	    do j=1,ny
	      do i=1,nx
		do k=1,nz-1
		  tz(k,i,j,1) = cof(i,j,k,5)
		  tz(k,i,j,2) = cof(i,j,k,7)
		  tz(k,i,j,3) = cof(i,j,k,6)
		end do
	      end do
	    end do
	    nxny = nx*ny
	    call factrp(nxny,nz,tz(1,1,1,1),tz(1,1,1,2),tz(1,1,1,3),
     +                  tz(1,1,1,4),tz(1,1,1,5),wk(kps))
	  end if
	end if
      end if
      return
      end

      subroutine adjmd3sa(nx,ny,nz,phi,cof,bndyc,sgx,sgy,sgz,xlmbda)
c
c     adjust rhs for solution in cof(i,j,k,8) on non-initial calls
c     (i.e., values in cof have not changed)
c
      implicit none
      integer nx,ny,nz
      real phi(0:nx+1,0:ny+1,0:nz+1),cof(nx,ny,nz,8)
      integer intl,nxa,nxb,nyc,nyd,nze,nzf,ixp,jyq,kzr,iex,jey,kez,
     +nfx,nfy,nfz,iguess,maxcy,method,meth2,nwork,lwork,itero,
     +kcycle,iprer,ipost,intpol
      real xa,xb,yc,yd,ze,zf,tolmax,relmax
      integer kpbgn,kcbgn,ktxbgn,ktybgn,ktzbgn,nxk,nyk,nzk,ngrid,
     +        klevel,kcur,kps
      common/imud3/intl,nxa,nxb,nyc,nyd,nze,nzf,ixp,jyq,kzr,iex,jey,kez,
     +nfx,nfy,nfz,iguess,maxcy,method,meth2,nwork,lwork,itero,
     +kcycle,iprer,ipost,intpol
      common/fmud3/xa,xb,yc,yd,ze,zf,tolmax,relmax
      common/mud3c/kpbgn(50),kcbgn(50),ktxbgn(50),ktybgn(50),
     +ktzbgn(50),nxk(50),nyk(50),nzk(50),ngrid,klevel,kcur,kps
      real dlx,dly,dlz,dlx2,dly2,dlz2,dlxx,dlyy,dlzz,odlxx,odlyy,odlzz
      real alfa,gbdy,x,y,z,c1,c2,c3,c4,c5,c6
      real sgxm,sgxp,sgym,sgyp,sgzm,sgzp
      integer i,j,k,ist,ifn,jst,jfn,kst,kfn,kbdy
      real sgx,sgy,sgz,xlmbda
      external sgx,sgy,sgz,xlmbda,bndyc
c
c     set current grid increments
c
      dlx = (xb-xa)/(nx-1)
      dlx2 = dlx+dlx
      dlxx = dlx*dlx
      dly = (yd-yc)/(ny-1)
      dly2 = dly+dly
      dlyy = dly*dly
      dlz = (zf-ze)/(nz-1)
      dlz2 = dlz+dlz
      dlzz = dlz*dlz
      odlxx = 1.0/dlxx
      odlyy = 1.0/dlyy
      odlzz = 1.0/dlzz
c
c     set x,y,z subscript limits for calls to sgx,sgy,sgz,xlmbda,bndyc
c
      jst = 1
      jfn = ny
      ist = 1
      ifn = nx
      kst = 1
      kfn = nz
      if (nxa.eq.1) ist = 2
      if (nxb.eq.1) ifn = nx-1
      if (nyc.eq.1) jst = 2
      if (nyd.eq.1) jfn = ny-1
      if (nze.eq.1) kst = 2
      if (nzf.eq.1) kfn = nz-1
c
c     adjust for derivative b.c.
c
      if (nxa.eq.2) then
	kbdy=1
	x=xa
	i = 1
	do k=kst,kfn
	  z=ze+(k-1)*dlz
	  do j=jst,jfn
	    y=yc+(j-1)*dly
	    call bndyc(kbdy,y,z,alfa,gbdy)
	    sgxm = sgx(x-0.5*dlx,y,z)
	    c1 = sgxm*odlxx
	    cof(i,j,k,8) = cof(i,j,k,8)+dlx2*c1*gbdy
	  end do
	end do
      end if
      if (nxb.eq.2) then
	kbdy=2
	x=xb
	i = nx
	do k=kst,kfn
	  z=ze+(k-1)*dlz
	  do j=jst,jfn
	    y=yc+(j-1)*dly
	    call bndyc(kbdy,y,z,alfa,gbdy)
	    sgxp = sgx(x+0.5*dlx,y,z)
	    c2 = sgxp*odlxx
	    cof(i,j,k,8) = cof(i,j,k,8)-dlx2*c2*gbdy
	  end do
	end do
      end if
      if (nyc.eq.2) then
	kbdy = 3
	y=yc
	j = 1
	do k=kst,kfn
	  z=ze+(k-1)*dlz
	  do i=ist,ifn
	    x=xa+(i-1)*dlx
	    call bndyc(kbdy,x,z,alfa,gbdy)
	    sgym = sgy(x,y-0.5*dly,z)
	    c3 = sgym*odlyy
	    cof(i,j,k,8) = cof(i,j,k,8)+dly2*c3*gbdy
	  end do
	end do
      end if
      if (nyd.eq.2) then
	kbdy = 4
	y=yd
	j = ny
	do k=kst,kfn
	  z=ze+(k-1)*dlz
	  do i=ist,ifn
	    x=xa+(i-1)*dlx
	    call bndyc(kbdy,x,z,alfa,gbdy)
	    sgyp = sgy(x,y+0.5*dly,z)
	    c4 = sgyp*odlyy
	    cof(i,j,k,8) = cof(i,j,k,8)-dly2*c4*gbdy
	  end do
	end do
      end if
      if (nze.eq.2) then
	kbdy = 5
	k = 1
	z=ze
	do j=jst,jfn
	  y=yc+(j-1)*dly
	  do i=ist,ifn
	    x=xa+(i-1)*dlx
	    call bndyc(kbdy,x,y,alfa,gbdy)
	    sgzm = sgz(x,y,z-0.5*dlz)
	    c5 = sgzm*odlzz
	    cof(i,j,k,8) = cof(i,j,k,8)+dlz2*c5*gbdy
	  end do
	end do
      end if
      if (nzf.eq.2) then
	kbdy = 6
	z=zf
	k = nz
	do j=jst,jfn
	  y=yc+(j-1)*dly
	  do i=ist,ifn
	    x=xa+(i-1)*dlx
	    call bndyc(kbdy,x,y,alfa,gbdy)
	    sgzp = sgz(x,y,z+0.5*dlz)
	    c6 = sgzp*odlzz
	    cof(i,j,k,8) = cof(i,j,k,8)-dlz2*c6*gbdy
	  end do
	end do
      end if
c
c     set specified b.c.
c
      if (nxa.eq.1) then
	i = 1
	do j=1,ny
	  do k=1,nz
	    cof(i,j,k,8) = phi(i,j,k)
	  end do
	end do
      end if
      if (nxb.eq.1) then
	i = nx
	do j=1,ny
	  do k=1,nz
	    cof(i,j,k,8) = phi(i,j,k)
	  end do
	end do
      end if
      if (nyc.eq.1) then
	j = 1
	do k=1,nz
	  do i=1,nx
	    cof(i,j,k,8) = phi(i,j,k)
	  end do
	end do
      end if
      if (nyd.eq.1) then
	j = ny
	do k=1,nz
	  do i=1,nx
	    cof(i,j,k,8) = phi(i,j,k)
	  end do
	end do
      end if
      if (nze.eq.1) then
	k = 1
	do j=1,ny
	  do i=1,nx
	    cof(i,j,k,8) = phi(i,j,k)
	  end do
	end do
      end if
      if (nzf.eq.1) then
	k = nz
	do j=1,ny
	  do i=1,nx
	    cof(i,j,k,8) = phi(i,j,k)
	  end do
	end do
      end if
      return
      end

      subroutine relmd3sa(wk)
c
c     use point or line relaxation in the x and/or y and/or z
c     or planar relaxation in the x,y or x,z or y,z planes
c
      implicit none
      real wk(*)
      integer intl,nxa,nxb,nyc,nyd,nze,nzf,ixp,jyq,kzr,iex,jey,kez,
     +nfx,nfy,nfz,iguess,maxcy,method,meth2,nwork,lwork,itero,
     +kcycle,iprer,ipost,intpol
      real xa,xb,yc,yd,ze,zf,tolmax,relmax
      integer kpbgn,kcbgn,ktxbgn,ktybgn,ktzbgn,nxk,nyk,nzk,ngrid,
     +        klevel,kcur,kps
      common/imud3/intl,nxa,nxb,nyc,nyd,nze,nzf,ixp,jyq,kzr,iex,jey,kez,
     +nfx,nfy,nfz,iguess,maxcy,method,meth2,nwork,lwork,itero,
     +kcycle,iprer,ipost,intpol
      common/fmud3/xa,xb,yc,yd,ze,zf,tolmax,relmax
      common/mud3c/kpbgn(50),kcbgn(50),ktxbgn(50),ktybgn(50),
     +ktzbgn(50),nxk(50),nyk(50),nzk(50),ngrid,klevel,kcur,kps
      integer nx,ny,nz,ip,ic,m,itx,ity,itz
      nx = nxk(klevel)
      ny = nyk(klevel)
      nz = nzk(klevel)
      ip = kpbgn(klevel)
      ic = kcbgn(klevel)
      if (method.eq.0) then
c
c     gauss-seidel pointwise red/black relaxation
c
	call relmd3sap(nx,ny,nz,wk(ip),wk(ic))
	return
      end if
      itx = ktxbgn(klevel)
      ity = ktybgn(klevel)
      itz = ktzbgn(klevel)
      m = method
c
c     check for line relaxation(s) (in combinations)
c
      if ((m-1)*(m-4)*(m-5)*(m-7) .eq. 0 ) then
c
c     line - x relaxation
c
	if (nxa .ne. 0 .or. nx .gt. 3) then
	 itx = ktxbgn(klevel)
	 call slxmd3(nx,ny,nz,wk(ip),wk(ic),wk(itx),wk(kps),nxa,nyc,nze)
	else
c
c     replace by point if x-periodic and nx=3
c
	  call relmd3sap(nx,ny,nz,wk(ip),wk(ic))
	end if
	if (method .eq. 1) return
      end if

      if ((m-2)*(m-4)*(m-6)*(m-7) .eq. 0 ) then
c
c     line - y relaxation
c
	if (nyc .ne. 0 .or. ny .gt. 3) then
	 ity = ktybgn(klevel)
	 call slymd3(nx,ny,nz,wk(ip),wk(ic),wk(ity),wk(kps),nxa,nyc,nze)
	else
c
c     replace by point if y-periodic and ny=3
c
	  call relmd3sap(nx,ny,nz,wk(ip),wk(ic))
	end if
	if ((m-2)*(m-4) .eq. 0) return
      end if

      if ((m-3)*(m-5)*(m-6)*(m-7) .eq. 0 ) then
c
c     line - z relaxation
c
	if (nze .ne. 0 .or. nz .gt. 3) then
	 itz = ktzbgn(klevel)
	 call slzmd3(nx,ny,nz,wk(ip),wk(ic),wk(itz),wk(kps),nxa,nyc,nze)
	else
c
c     replace by point if z-periodic and nz=3
c
	call relmd3sap(nx,ny,nz,wk(ip),wk(ic))
	end if
	return
      end if
c
c     planar relaxation
c
      if (method.eq.8) then
c
c     planar relaxation in the x,y plane
c
	call planxy(wk)
      else if (method.eq.9) then
c
c     planar relaxation in the x,z plane
c
	call planxz(wk)
      else if (method.eq.10) then
c
c     planar relaxation in the y,z plane
c
	call planyz(wk)
      end if
      return
      end

      subroutine relmd3sap(nx,ny,nz,phi,cof)
c
c     gauss-seidel point relaxation with red/black ordering
c     in three dimensions for nonseparable pde
c
      implicit none
      integer nx,ny,nz
      integer intl,nxa,nxb,nyc,nyd,nze,nzf,ixp,jyq,kzr,iex,jey,kez,
     +nfx,nfy,nfz,iguess,maxcy,method,meth2,nwork,lwork,itero,
     +kcycle,iprer,ipost,intpol
      common/imud3/intl,nxa,nxb,nyc,nyd,nze,nzf,ixp,jyq,kzr,iex,jey,kez,
     +nfx,nfy,nfz,iguess,maxcy,method,meth2,nwork,lwork,itero,
     +kcycle,iprer,ipost,intpol
      real phi(0:nx+1,0:ny+1,0:nz+1),cof(nx,ny,nz,8)
      integer i,j,k,nper
c
c     set periodic b.c. indicator
c
      nper = nxa*nyc*nze
c
c     set periodic boundaries as necessary
c
      if (nper.eq.0) call per3vb(nx,ny,nz,phi,nxa,nyc,nze)
c
c     relax in order:
c     (1) red (x,y) on odd z planes
c     (2) black (x,y) on even z planes
c     (3) black (x,y) on odd z planes
c     (4) red (x,y) on even z planes
c
!$OMP PARALLEL DO PRIVATE(i,j,k), SHARED(phi,cof,nx,ny,nz)
      do k=1,nz,2
c
c     red (x,y) points on odd z planes
c
	do i=1,nx,2
	  do j=1,ny,2
	    phi(i,j,k) = (cof(i,j,k,8) - (
     +                    cof(i,j,k,1)*phi(i-1,j,k)+
     +                    cof(i,j,k,2)*phi(i+1,j,k)+
     +                    cof(i,j,k,3)*phi(i,j-1,k)+
     +                    cof(i,j,k,4)*phi(i,j+1,k)+
     +                    cof(i,j,k,5)*phi(i,j,k-1)+
     +                    cof(i,j,k,6)*phi(i,j,k+1)))
     +                   /cof(i,j,k,7)
	  end do
	end do
	do i=2,nx,2
	  do j=2,ny,2
	    phi(i,j,k) = (cof(i,j,k,8) - (
     +                    cof(i,j,k,1)*phi(i-1,j,k)+
     +                    cof(i,j,k,2)*phi(i+1,j,k)+
     +                    cof(i,j,k,3)*phi(i,j-1,k)+
     +                    cof(i,j,k,4)*phi(i,j+1,k)+
     +                    cof(i,j,k,5)*phi(i,j,k-1)+
     +                    cof(i,j,k,6)*phi(i,j,k+1)))
     +                   /cof(i,j,k,7)
	  end do
	end do
      end do
      if (nper.eq.0) call per3vb(nx,ny,nz,phi,nxa,nyc,nze)
c
c    black (x,y) points on even z planes
c
!$OMP PARALLEL DO PRIVATE(i,j,k), SHARED(phi,cof,nx,ny,nz)
      do k=2,nz,2
	do i=1,nx,2
	  do j=2,ny,2
	    phi(i,j,k) = (cof(i,j,k,8) - (
     +                    cof(i,j,k,1)*phi(i-1,j,k)+
     +                    cof(i,j,k,2)*phi(i+1,j,k)+
     +                    cof(i,j,k,3)*phi(i,j-1,k)+
     +                    cof(i,j,k,4)*phi(i,j+1,k)+
     +                    cof(i,j,k,5)*phi(i,j,k-1)+
     +                    cof(i,j,k,6)*phi(i,j,k+1)))
     +                   /cof(i,j,k,7)
	  end do
	end do
	do i=2,nx,2
	  do j=1,ny,2
	    phi(i,j,k) = (cof(i,j,k,8) - (
     +                    cof(i,j,k,1)*phi(i-1,j,k)+
     +                    cof(i,j,k,2)*phi(i+1,j,k)+
     +                    cof(i,j,k,3)*phi(i,j-1,k)+
     +                    cof(i,j,k,4)*phi(i,j+1,k)+
     +                    cof(i,j,k,5)*phi(i,j,k-1)+
     +                    cof(i,j,k,6)*phi(i,j,k+1)))
     +                   /cof(i,j,k,7)
	  end do
	end do
      end do
      if (nper.eq.0) call per3vb(nx,ny,nz,phi,nxa,nyc,nze)
c
c     black (x,y) points on odd z planes
c
!$OMP PARALLEL DO PRIVATE(i,j,k), SHARED(phi,cof,nx,ny,nz)
      do k=1,nz,2
	do i=1,nx,2
	  do j=2,ny,2
	    phi(i,j,k) = (cof(i,j,k,8) - (
     +                    cof(i,j,k,1)*phi(i-1,j,k)+
     +                    cof(i,j,k,2)*phi(i+1,j,k)+
     +                    cof(i,j,k,3)*phi(i,j-1,k)+
     +                    cof(i,j,k,4)*phi(i,j+1,k)+
     +                    cof(i,j,k,5)*phi(i,j,k-1)+
     +                    cof(i,j,k,6)*phi(i,j,k+1)))
     +                   /cof(i,j,k,7)
	  end do
	end do
	do i=2,nx,2
	  do j=1,ny,2
	    phi(i,j,k) = (cof(i,j,k,8) - (
     +                    cof(i,j,k,1)*phi(i-1,j,k)+
     +                    cof(i,j,k,2)*phi(i+1,j,k)+
     +                    cof(i,j,k,3)*phi(i,j-1,k)+
     +                    cof(i,j,k,4)*phi(i,j+1,k)+
     +                    cof(i,j,k,5)*phi(i,j,k-1)+
     +                    cof(i,j,k,6)*phi(i,j,k+1)))
     +                   /cof(i,j,k,7)
	  end do
	end do
      end do
      if (nper.eq.0) call per3vb(nx,ny,nz,phi,nxa,nyc,nze)
c
c     red (x,y) points on even z planes
c
!$OMP PARALLEL DO PRIVATE(i,j,k), SHARED(phi,cof,nx,ny,nz)
      do k=2,nz,2
	do i=1,nx,2
	  do j=1,ny,2
	    phi(i,j,k) = (cof(i,j,k,8) - (
     +                    cof(i,j,k,1)*phi(i-1,j,k)+
     +                    cof(i,j,k,2)*phi(i+1,j,k)+
     +                    cof(i,j,k,3)*phi(i,j-1,k)+
     +                    cof(i,j,k,4)*phi(i,j+1,k)+
     +                    cof(i,j,k,5)*phi(i,j,k-1)+
     +                    cof(i,j,k,6)*phi(i,j,k+1)))
     +                   /cof(i,j,k,7)
	  end do
	end do
	do i=2,nx,2
	  do j=2,ny,2
	    phi(i,j,k) = (cof(i,j,k,8) - (
     +                    cof(i,j,k,1)*phi(i-1,j,k)+
     +                    cof(i,j,k,2)*phi(i+1,j,k)+
     +                    cof(i,j,k,3)*phi(i,j-1,k)+
     +                    cof(i,j,k,4)*phi(i,j+1,k)+
     +                    cof(i,j,k,5)*phi(i,j,k-1)+
     +                    cof(i,j,k,6)*phi(i,j,k+1)))
     +                   /cof(i,j,k,7)
	  end do
	end do
      end do
c
c     final set of periodic virtual boundaries if necessary
c
      if (nper.eq.0) call per3vb(nx,ny,nz,phi,nxa,nyc,nze)
      return
      end

max_stars_repo_name: utastudents/selalib

max_issues_repo_name: utastudents/selalib

max_forks_repo_name: utastudents/selalib


text:
C$Procedure PLTNRM ( DSK, compute outward normal of plate )

      SUBROUTINE PLTNRM ( V1, V2, V3, NORMAL )

C$ Abstract
C
C     Compute an outward normal vector of a triangular plate.
C     The vector does not necessarily have unit length.
C
C$ Disclaimer
C
C     THIS SOFTWARE AND ANY RELATED MATERIALS WERE CREATED BY THE
C     CALIFORNIA INSTITUTE OF TECHNOLOGY (CALTECH) UNDER A U.S.
C     GOVERNMENT CONTRACT WITH THE NATIONAL AERONAUTICS AND SPACE
C     ADMINISTRATION (NASA). THE SOFTWARE IS TECHNOLOGY AND SOFTWARE
C     PUBLICLY AVAILABLE UNDER U.S. EXPORT LAWS AND IS PROVIDED "AS-IS"
C     TO THE RECIPIENT WITHOUT WARRANTY OF ANY KIND, INCLUDING ANY
C     WARRANTIES OF PERFORMANCE OR MERCHANTABILITY OR FITNESS FOR A
C     PARTICULAR USE OR PURPOSE (AS SET FORTH IN UNITED STATES UCC
C     SECTIONS 2312-2313) OR FOR ANY PURPOSE WHATSOEVER, FOR THE
C     SOFTWARE AND RELATED MATERIALS, HOWEVER USED.
C
C     IN NO EVENT SHALL CALTECH, ITS JET PROPULSION LABORATORY, OR NASA
C     BE LIABLE FOR ANY DAMAGES AND/OR COSTS, INCLUDING, BUT NOT
C     LIMITED TO, INCIDENTAL OR CONSEQUENTIAL DAMAGES OF ANY KIND,
C     INCLUDING ECONOMIC DAMAGE OR INJURY TO PROPERTY AND LOST PROFITS,
C     REGARDLESS OF WHETHER CALTECH, JPL, OR NASA BE ADVISED, HAVE
C     REASON TO KNOW, OR, IN FACT, SHALL KNOW OF THE POSSIBILITY.
C
C     RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF
C     THE SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY
C     CALTECH AND NASA FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE
C     ACTIONS OF RECIPIENT IN THE USE OF THE SOFTWARE.
C
C$ Required_Reading
C
C     DSK
C
C$ Keywords
C
C     DSK
C     FILES
C     TOPOGRAPHY
C
C$ Declarations

      IMPLICIT NONE

      DOUBLE PRECISION      V1     ( 3 )
      DOUBLE PRECISION      V2     ( 3 )
      DOUBLE PRECISION      V3     ( 3 )
      DOUBLE PRECISION      NORMAL ( 3 )

C$ Brief_I/O
C
C     Variable  I/O  Description
C     --------  ---  --------------------------------------------------
C     V1,
C     V2,
C     V3         I   Vertices of a plate.
C     NORMAL     O   Plate's outward normal vector.
C
C$ Detailed_Input
C
C     V1,
C     V2,
C     V3             are vertices of a triangular plate.
C
C$ Detailed_Output
C
C     NORMAL         is an outward normal vector of the plate defined by
C                    the input vertices. The order of the vertices is
C                    used to determine the choice of normal direction:
C                    the normal vector is
C
C                       ( V2 - V1 ) x ( V3 - V2 )
C
C
C$ Parameters
C
C     None.
C
C$ Exceptions
C
C     1) The input plate may be degenerate: it may be a line segment
C        or a point. These are not considered to be erroneous inputs.
C
C$ Files
C
C     None.
C
C$ Particulars
C
C     This routine saves computation time by not scaling the output
C     vector to unit length. The caller can scale the vector using
C     the routine VHAT.
C
C$ Examples
C
C
C     1) Compute an upward normal of an equilateral triangle lying
C        in the X-Y plane and centered at the origin.
C
C        Example code begins here.
C
C
C           PROGRAM EX1
C           IMPLICIT NONE
C
C           DOUBLE PRECISION      NORMAL ( 3 )
C           DOUBLE PRECISION      S
C           DOUBLE PRECISION      V1     ( 3 )
C           DOUBLE PRECISION      V2     ( 3 )
C           DOUBLE PRECISION      V3     ( 3 )
C
C
C           S = SQRT(3.D0)/2
C
C           CALL VPACK (    S,  -0.5D0,  0.D0, V1 )
C           CALL VPACK ( 0.D0,    1.D0,  0.D0, V2 )
C           CALL VPACK (   -S,  -0.5D0,  0.D0, V3 )
C
C
C           CALL PLTNRM ( V1, V2, V3, NORMAL )
C
C           WRITE (*, '(1X,A,3(3E20.12))' ) 'NORMAL = ', NORMAL
C
C           END
C
C        When run on a PC/Linux/gfortran/64-bit platform, the output
C        from this program was:
C
C
C  NORMAL =   0.000000000000E+00  0.000000000000E+00  0.259807621135E+01
C
C
C$ Restrictions
C
C     None.
C
C$ Literature_References
C
C     None.
C
C$ Author_and_Institution
C
C     N.J. Bachman    (JPL)
C
C$ Version
C
C-    SPICELIB Version 1.0.0, 26-JAN-2016 (NJB)
C
C-&

C$ Index_Entries
C
C     compute normal vector of triangular plate from vertices
C
C-&

C
C     SPICELIB functions
C

C
C     Local variables
C
      DOUBLE PRECISION      EDGE1  ( 3 )
      DOUBLE PRECISION      EDGE2  ( 3 )

C
C     This routine is error-free.
C

C
C     Type 2 plate vertices are ordered in the positive
C     (right-handed) sense about the outward normal.
C
      CALL VSUB ( V2, V1, EDGE1 )
      CALL VSUB ( V3, V2, EDGE2 )

      CALL VCRSS ( EDGE1, EDGE2, NORMAL )

      END

max_stars_repo_name: agforero/FTFramework

max_issues_repo_name: agforero/fortran-testing-framework

max_forks_repo_name: agforero/fortran-testing-framework