text:
! DART software - Copyright UCAR. This open source software is provided
! by UCAR, "as is", without charge, subject to all terms of use at
! http://www.image.ucar.edu/DAReS/DART/DART_download
!
! $Id$

program CONAGUA_convert_streamflow

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! James McCreight jamesmcc at ucar dot edu
! 9/24/2018
! Modified by TJH 12/2019 - use error_handler
!
! Bring in the CONAGUA daily or hourly streamflow (gasto) observations
! 1. Read the meta data file specified in
!       input.nml:CONAGUA_convert_streamflow_nml:meta_data_files
! 2. Process the files listed in
!       intput.nml:CONAGUA_convert_streamflow_nwml:data_file_list
!    All files must either be daily dd*.csv or hourly hd*.csv
!    If a given gage is not in the meta data file, fatal error.
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

use         types_mod, only : r8

use      location_mod, only : VERTISSURFACE

use     utilities_mod, only : initialize_utilities, finalize_utilities,      &
                              open_file, close_file, find_namelist_in_file,  &
                              check_namelist_read, nmlfileunit, do_nml_file, &
                              do_nml_term, error_handler, E_ERR, E_MSG

use  time_manager_mod, only : time_type, set_calendar_type, set_date, set_time, &
                              operator(>=), increment_time, get_time,           &
                              operator(-), GREGORIAN, operator(+), print_date

use  obs_sequence_mod, only : obs_sequence_type, obs_type, read_obs_seq,     &
                              static_init_obs_sequence, init_obs,            &
                              write_obs_seq, init_obs_sequence, get_num_obs, &
                              set_copy_meta_data, set_qc_meta_data

use      obs_kind_mod, only : STREAM_FLOW

use obs_def_streamflow_mod, only : set_streamflow_metadata, &
                                   missing_gage_ID, missing_link_ID

use obs_utilities_mod, only : create_3d_obs, add_obs_to_seq

implicit none

! version controlled file description for error handling, do not edit
character(len=*), parameter :: source   = 'CONAGUA_convert_streamflow.f90'
character(len=*), parameter :: revision = '$Revision$'
character(len=*), parameter :: revdate  = '$Date$'

!------------------------------------------------------------------------
!  Declare namelist parameters
!------------------------------------------------------------------------

character(len=256) :: meta_data_file         = 'unset'
character(len=256) :: data_file_list         = 'unset'
character(len=256) :: obs_out_file           = 'obs_seq.out'
real(r8)           :: obs_fraction_for_error = 0.01_r8
real(r8)           :: obs_min_err_std        = 0.5_r8
logical            :: debug                  = .false.  ! .true. print info

namelist /CONAGUA_convert_streamflow_nml/  &
     meta_data_file, data_file_list, obs_out_file, &
     obs_fraction_for_error, obs_min_err_std, debug

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

integer :: oday, osec, rcio, iunit, io
integer :: num_copies, num_qc, max_obs

logical  :: first_obs

real(r8) :: obs_err, qc
real(r8) :: lat, lon, vert

type(obs_sequence_type) :: obs_seq
type(obs_type)          :: obs, prev_obs
type(time_type)         :: time_obs, prev_time
integer :: meta_unit, list_unit, stn_unit
integer :: key

integer, parameter :: MAX_STATIONS = 5000
integer, parameter :: MAX_RECORDS = 10000
integer :: i_rec, i_stn
character(len=512) :: line_in, stn_name, station_station
integer :: count_stations, count_meta_stations, ii, stn_name_len, wh_meta_stn

logical :: hourly_resolution

integer :: year, month, day
real(r8) :: obs_val
real(r8), parameter :: MISSING_FLOAT = -88888.0_r8
integer,  parameter :: MISSING_INT = -88888

! Meta data derived type/structure.
!FID,LON,LAT,STATION,Name,HOST
!1,-92.75,17.38333,Oxolotan,30111,BANDAS
type :: meta_data_record
   integer :: fid
   real :: lon
   real :: lat
   character(len=64) :: station
   character(len=64) :: name
   character(len=64) :: host
end type meta_data_record

type(meta_data_record), dimension(MAX_STATIONS) :: meta_data
integer, parameter :: N_META_FIELDS = 6
character(len=32) :: meta_header(N_META_FIELDS)

! data_file_list
character(len=256) :: data_files(MAX_STATIONS)
character(len=256) :: stn_file_name, stn_resolution

! Daily data derived type/structure.
type :: station_daily_data_record
   integer :: year
   integer :: month
   real :: day_of_month_streamflow(31)
end type station_daily_data_record

! Will process one in and one out, so the dimension is only len 1.
type(station_daily_data_record) :: stn_daily_record
type(station_daily_data_record) :: stn_daily_record_null
integer, parameter :: N_STN_DAILY_FIELDS = 33
character(len=32) :: stn_daily_header(N_STN_DAILY_FIELDS)

character(len=512) :: string1, string2, string3

! -------------------------------------------------------
! Start work
! Initialize DART modules.
call initialize_utilities('CONAGUA_convert_streamflow')

call find_namelist_in_file('input.nml', 'CONAGUA_convert_streamflow_nml', iunit)
read(iunit, nml = CONAGUA_convert_streamflow_nml, iostat = io)
call check_namelist_read(iunit, io, 'CONAGUA_convert_streamflow_nml')

! Record the namelist values used for the run
if (do_nml_file()) write(nmlfileunit, nml=CONAGUA_convert_streamflow_nml)
if (do_nml_term()) write(     *     , nml=CONAGUA_convert_streamflow_nml)

! time setup
call set_calendar_type(GREGORIAN)

! each observation in this series will have a single observation value
! and a quality control flag.  the max possible number of obs needs to
! be specified but it will only write out the actual number created.
max_obs    = 100000
num_copies = 1
num_qc     = 1

! call the initialization code, and initialize two empty observation types
call static_init_obs_sequence()
call init_obs(obs,      num_copies, num_qc)
call init_obs(prev_obs, num_copies, num_qc)
first_obs = .true.

! create a new, empty obs_seq file.  you must give a max limit
! on number of obs.  increase the size if too small.
call init_obs_sequence(obs_seq, num_copies, num_qc, max_obs)

! the first one needs to contain the string 'observation' and the
! second needs the string 'QC'.
call set_copy_meta_data(obs_seq, 1, 'observation')
call set_qc_meta_data(obs_seq, 1, 'Data QC')

! Several obs_seq quantities which are static for our observations:

! surface obs.  normally we set vert to the actual surface elevation,
! we do not have it here, so set to 0 for now.
vert = 0.0_r8

! Set the DART data quality control.   0 is good data.
! increasingly larger QC values are more questionable quality data.
qc = 0.0_r8

! -------------------------------------------------------
! Read meta data file
meta_data_file = trim(meta_data_file)
meta_unit = open_file(meta_data_file, 'formatted', 'read')

if (debug)  then
   write(string1,*) 'opened meta data file "'//trim(meta_data_file)//'"'
   call error_handler(E_MSG, source, string1)
endif

read(meta_unit, *) meta_header

count_meta_stations = 0
meta_loop: do i_rec = 1, MAX_STATIONS
   read(meta_unit, *, iostat=rcio) meta_data(i_rec)
   if (rcio /= 0) then
      if (debug) then
         write(string1,*)'meta_loop: got bad read code getting line, rcio = ', rcio
         call error_handler(E_MSG, source, string1)
      endif
      exit meta_loop
   endif
   count_meta_stations = count_meta_stations + 1
end do meta_loop
call close_file(meta_unit)

print *, ''
print *, '----------------------------------'
print *, 'reading meta data file: ', trim(meta_data_file)
print *, ''
print *, 'The following header fields were found: '
do ii=1, N_META_FIELDS
    write(*, fmt="(1x,a)", advance="no") trim(meta_header(ii)) // ', '
end do
print *,''
print *,''
print *, "count_meta_stations: ", count_meta_stations
do ii = 1, count_meta_stations
   print*, trim(meta_data(ii)%name), ': ', trim(meta_data(ii)%station)
end do
print *,''

! -------------------------------------------------------
! Pre-Process the data_file_list and check it's consistent with
! 1) meta_data
! 2) itself (all files at the same time resolution)

list_unit = open_file(data_file_list, 'formatted', 'read')

if (debug)  then
   write(string1,*) 'opened data_file_list "'//trim(data_file_list)//'"'
   call error_handler(E_MSG, source, string1)
endif

count_stations = 0
list_loop: do i_rec = 1, MAX_STATIONS
   read(list_unit, '(a)', iostat=rcio) data_files(i_rec)
   if (rcio /= 0) then
      if (debug) then
         write(string1,*)'list_loop: got bad read code getting line, rcio = ', rcio
         call error_handler(E_MSG, source, string1)
      endif
      exit list_loop
   endif
   count_stations = count_stations + 1
end do list_loop
call close_file(list_unit)

! Validate that all the stations in the list are in the meta data.

valid_loop: do ii = 1, count_stations
   stn_file_name = data_files(ii)
   stn_name_len  = len_trim(stn_file_name)

   ! As far as we can tell, CONAGUA stations are always 5 digits wide.
   stn_resolution = stn_file_name((stn_name_len-10):(stn_name_len-9))
   stn_name       = stn_file_name((stn_name_len-8):(stn_name_len-4))

   if (debug) then
      write(string1,*)'stn_file_name substring is ', stn_file_name((stn_name_len-10):(stn_name_len-4))
      write(string2,*)'station resolution is <',trim(stn_resolution),'>'
      write(string3,*)'station name       is <',trim(stn_name),'>'
      call error_handler(E_MSG, source, string1, text2=string2, text3=string3)
   endif

   if (.not. any(meta_data%name == stn_name)) then
      write(string1,*)'Metadata not found for the requested station file: "'//trim(stn_file_name)//'"'
      call error_handler(E_ERR, source, string1)
   end if

   if (ii == 1) then
      hourly_resolution = (stn_resolution == 'hd')
   else
      if (hourly_resolution .and. stn_resolution .ne. 'hd') then
         write(string1,*)"hourly resolution does not appear consistent across the data_file_list.txt"
         call error_handler(E_ERR, source, string1)
      end if
      if (.not. hourly_resolution .and. stn_resolution .ne. 'dd') then
         write(string1,*) "daily resolution does not appear consistent across the data_file_list.txt"
         call error_handler(E_ERR, source, string1)
      end if
   end if

end do valid_loop

! -------------------------------------------------------
! Daily observation processing.

stn_daily_record_null%year  = MISSING_INT
stn_daily_record_null%month = MISSING_INT
stn_daily_record_null%day_of_month_streamflow = MISSING_FLOAT

station_loop: do i_stn = 1, count_stations

   stn_file_name = data_files(i_stn)
   stn_name_len = len(trim(stn_file_name))
   stn_name = stn_file_name((stn_name_len-8):(stn_name_len-4))

   wh_meta_stn = 0
   wh_loop: do ii = 1, count_meta_stations
      if (trim(meta_data(ii)%name) .eq. trim(stn_name)) then
         wh_meta_stn = ii
         exit wh_loop
      end if
   end do wh_loop

   station_station = meta_data(wh_meta_stn)%station
   lon = meta_data(wh_meta_stn)%lon
   if (lon < 0.0_r8) lon = lon + 360.0_r8
   lat = meta_data(wh_meta_stn)%lat

   write(string1,*) ''
   write(string2,*) '!-----------------------------------'
   write(string3,*) 'Processing station ', i_stn, ': "'//trim(station_station)//'"'
   call error_handler(E_MSG,source,string1,text2=string2,text3=string3)

   stn_unit = open_file(trim(stn_file_name), 'formatted', 'read')

   write(string1,*) 'opened stn record file: "'//trim(stn_file_name)//'"'
   call error_handler(E_MSG,source,string1)

   read(stn_unit, *, iostat=rcio) stn_daily_header

   record_loop: do i_rec = 1, MAX_RECORDS
      stn_daily_record = stn_daily_record_null
      read(stn_unit, '(a)', iostat=rcio) line_in
      if (rcio /= 0) then
         if (debug) then
            write(string1,*) 'recordloop: unable to read record for "'//trim(stn_name)//'"'
            write(string2,*) ' rcio = ', rcio, '; record number = ', i_rec
            call error_handler(E_MSG,source,string1,text2=string2)
         endif
         exit record_loop
      endif

      ! This is a band-aid for records/lines which end with a missing value.
      ! I'd prefer to keep it to a single comma so any "badly mangled" records might cause an error.
      write(line_in, '(a)') trim(line_in) // ','

      read(line_in, *, iostat=rcio) stn_daily_record
      if (rcio /= 0) then
         if (debug) then
            write(string1,*) 'recordloop: unable to parse stn_daily_record for "'//trim(stn_name)//'"'
            write(string2,*) 'rcio = ', rcio, '; record number = ', i_rec
            call error_handler(E_MSG,source,string1,text2=string2)
         endif
         exit record_loop
      endif

      year  = stn_daily_record%year
      month = stn_daily_record%month

      day_loop: do day = 1, 31

         obs_val = stn_daily_record%day_of_month_streamflow(day)

         if (obs_val == MISSING_FLOAT .or. obs_val <= 0.000001_r8) then
            cycle day_loop ! skip this recrod.
         end if

         ! Set the observation error variance here.
         obs_err = max(obs_val*obs_fraction_for_error, obs_min_err_std)

         time_obs = set_date(year, month, day, 0, 0, 0)
         ! extract time of observation into gregorian day, sec.
         call get_time(time_obs, osec, oday)

         ! This routine can help tie gages to model spatial elements which can then
         ! be used for more advanced localization.
         ! Currently it is not used for gridded, though a layer in the Fulldom.nc file
         ! could easily facilitate this. The bigger issue with gridded is lack of
         ! connectivity approach currently (only calculated inside the model).
         ! Just supplying dummy values now...

         call set_streamflow_metadata(key, missing_gage_ID, missing_link_ID)

         call create_3d_obs(                 &
              lat, lon, vert, VERTISSURFACE, &
              obs_val,                       &
              STREAM_FLOW,                   &
              obs_err, oday, osec, qc, obs,  &
              key                            )

         call add_obs_to_seq(obs_seq, obs, time_obs, prev_obs, prev_time, first_obs)

      end do day_loop

   end do record_loop

   call close_file(stn_unit)

end do station_loop

write(string1,*)'writing "'//trim(obs_out_file)//'"'
write(string2,*)'obs_count = ', get_num_obs(obs_seq)

! if we added any obs to the sequence, write it out to a file now.
if ( get_num_obs(obs_seq) > 0 ) then
   call error_handler(E_MSG,source,string1,text2=string2)
   call write_obs_seq(obs_seq, obs_out_file)
else
   write(string3,*)'no observations in output. Something is probably wrong.'
   call error_handler(E_ERR,source,string1,text2=string2,text3=string3)
endif

! end of main program
call finalize_utilities()

end program CONAGUA_convert_streamflow

max_stars_repo_name: hkershaw-brown/feature-preprocess

max_issues_repo_name: hkershaw-brown/feature-preprocess

max_forks_repo_name: hkershaw-brown/feature-preprocess


text:
C     Last change:  ERB  22 Oct 2002    2:04 pm
      SUBROUTINE GWF1DRN6ALP(ISUM,LCDRAI,MXDRN,NDRAIN,IN,IOUT,IDRNCB,
     1        NDRNVL,IDRNAL,IFREFM,NPDRN,IDRNPB,NNPDRN,NOPRDR)
C
C-----VERSION 11JAN2000 GWF1DRN6ALP
C     ******************************************************************
C     ALLOCATE ARRAY STORAGE FOR DRAINS
C     ******************************************************************
C
C     SPECIFICATIONS:
C     ------------------------------------------------------------------
      COMMON /DRNCOM/DRNAUX(5)
      CHARACTER*16 DRNAUX
      CHARACTER*200 LINE
C     ------------------------------------------------------------------
C
C1------IDENTIFY PACKAGE AND INITIALIZE NDRAIN.
      WRITE(IOUT,1)IN
    1 FORMAT(1X,/1X,'DRN6 -- DRAIN PACKAGE, VERSION 6, 1/11/2000',
     1' INPUT READ FROM UNIT ',I4)
      NDRAIN=0
      NNPDRN=0
      NOPRDR=0
C
C2------READ MAXIMUM NUMBER OF DRAINS AND UNIT OR FLAG FOR
C2------CELL-BY-CELL FLOW TERMS.
      CALL URDCOM(IN,IOUT,LINE)
      CALL UPARLSTAL(IN,IOUT,LINE,NPDRN,MXPD)
      IF(IFREFM.EQ.0) THEN
         READ(LINE,'(2I10)') MXACTD,IDRNCB
         LLOC=21
      ELSE
         LLOC=1
         CALL URWORD(LINE,LLOC,ISTART,ISTOP,2,MXACTD,R,IOUT,IN)
         CALL URWORD(LINE,LLOC,ISTART,ISTOP,2,IDRNCB,R,IOUT,IN)
      END IF
      WRITE(IOUT,3) MXACTD
    3 FORMAT(1X,'MAXIMUM OF ',I6,' ACTIVE DRAINS AT ONE TIME')
      IF(IDRNCB.LT.0) WRITE(IOUT,7)
    7 FORMAT(1X,'CELL-BY-CELL FLOWS WILL BE PRINTED WHEN ICBCFL NOT 0')
         IF(IDRNCB.GT.0) WRITE(IOUT,8) IDRNCB
    8 FORMAT(1X,'CELL-BY-CELL FLOWS WILL BE SAVED ON UNIT ',I4)
C
C3------READ AUXILIARY VARIABLES AND CBC ALLOCATION OPTION.
      IDRNAL=0
      NAUX=0
   10 CALL URWORD(LINE,LLOC,ISTART,ISTOP,1,N,R,IOUT,IN)
      IF(LINE(ISTART:ISTOP).EQ.'CBCALLOCATE' .OR.
     1   LINE(ISTART:ISTOP).EQ.'CBC') THEN
         IDRNAL=1
         WRITE(IOUT,11)
   11    FORMAT(1X,'MEMORY IS ALLOCATED FOR CELL-BY-CELL BUDGET TERMS')
         GO TO 10
      ELSE IF(LINE(ISTART:ISTOP).EQ.'AUXILIARY' .OR.
     1        LINE(ISTART:ISTOP).EQ.'AUX') THEN
         CALL URWORD(LINE,LLOC,ISTART,ISTOP,1,N,R,IOUT,IN)
         IF(NAUX.LT.5) THEN
            NAUX=NAUX+1
            DRNAUX(NAUX)=LINE(ISTART:ISTOP)
            WRITE(IOUT,12) DRNAUX(NAUX)
   12       FORMAT(1X,'AUXILIARY DRAIN VARIABLE: ',A)
         END IF
         GO TO 10
      ELSE IF(LINE(ISTART:ISTOP).EQ.'NOPRINT') THEN
         WRITE(IOUT,13)
   13    FORMAT(1X,'LISTS OF DRAIN CELLS WILL NOT BE PRINTED')
         NOPRDR = 1
         GO TO 10
      END IF
      NDRNVL=5+NAUX+IDRNAL
C
C4------ALLOCATE SPACE IN THE RX ARRAY FOR THE DRAI ARRAY.
      IDRNPB=MXACTD+1
      MXDRN=MXACTD+MXPD
      ISOLD=ISUM
      LCDRAI=ISUM
      ISUM=ISUM+NDRNVL*MXDRN
C
C5------PRINT AMOUNT OF SPACE USED BY DRAIN PACKAGE.
      ISP=ISUM-ISOLD
      WRITE (IOUT,14)ISP
   14 FORMAT(1X,I10,' ELEMENTS IN RX ARRAY ARE USED BY DRN')
C
C6------RETURN.
      RETURN
      END
      SUBROUTINE GWF1DRN6RPPD(IN,IOUT,NDRNVL,IDRNAL,NCOL,NROW,NLAY,
     1            NPDRN,DRAI,IDRNPB,MXDRN,IFREFM,ITERP,INAMLOC,NOPRDR)
C
C-----VERSION 20011108 GWF1DRN6RPPD
C     ******************************************************************
C     READ DRAIN PARAMETERS
C     ******************************************************************
C     Modified 11/8/2001 to support parameter instances - ERB
C
C     SPECIFICATIONS:
C     ------------------------------------------------------------------
      DIMENSION DRAI(NDRNVL,MXDRN)
      COMMON /DRNCOM/DRNAUX(5)
      CHARACTER*16 DRNAUX
C     ------------------------------------------------------------------
C
C-------READ NAMED PARAMETERS.
      IF (ITERP.EQ.1) WRITE(IOUT,1000) NPDRN
 1000 FORMAT(1X,//1X,I5,' Drain parameters')
      IF(NPDRN.GT.0) THEN
        NAUX=NDRNVL-5-IDRNAL
        LSTSUM=IDRNPB
        ITERPU = ITERP
        IF (NOPRDR .EQ.1) ITERPU = 99
        DO 20 K=1,NPDRN
          LSTBEG=LSTSUM
          CALL UPARLSTRP(LSTSUM,MXDRN,IN,IOUT,IP,'DRN','DRN',ITERP,
     &                   NUMINST,INAMLOC)
          NLST=LSTSUM-LSTBEG
          IF (NUMINST.GT.1) NLST = NLST/NUMINST
C         ASSIGN STARTING INDEX FOR READING INSTANCES
          IF (NUMINST.EQ.0) THEN
            IB=0
          ELSE
            IB=1
          ENDIF
C         READ LIST(S) OF CELLS, PRECEDED BY INSTANCE NAME IF NUMINST>0
          LB=LSTBEG
          DO 10 I=IB,NUMINST
            IF (I.GT.0) THEN
              CALL UINSRP(I,IN,IOUT,IP,ITERP)
            ENDIF
            CALL ULSTRD(NLST,DRAI,LB,NDRNVL,MXDRN,IDRNAL,IN,IOUT,
     &      'DRAIN NO.  LAYER   ROW   COL     DRAIN EL.  STRESS FACTOR',
     &        DRNAUX,5,NAUX,IFREFM,NCOL,NROW,NLAY,5,5,ITERPU)
            LB=LB+NLST
   10     CONTINUE
   20   CONTINUE
      END IF
C
C6------RETURN
      RETURN
      END
      SUBROUTINE GWF1DRN6RPSS(DRAI,NDRAIN,MXDRN,IN,IOUT,NDRNVL,IDRNAL,
     1       IFREFM,NCOL,NROW,NLAY,NNPDRN,NPDRN,IDRNPB,NOPRDR)
C
C-----VERSION 11JAN2000 GWF1DRN6RPSS
C     ******************************************************************
C     READ DRAIN HEAD, CONDUCTANCE AND BOTTOM ELEVATION
C     ******************************************************************
C
C     SPECIFICATIONS:
C     ------------------------------------------------------------------
      DIMENSION DRAI(NDRNVL,MXDRN)
      COMMON /DRNCOM/DRNAUX(5)
      CHARACTER*16 DRNAUX
C     ------------------------------------------------------------------
C
C1------READ ITMP (NUMBER OF DRAINS OR FLAG TO REUSE DATA) AND
C1------NUMBER OF PARAMETERS.
      IF(NPDRN.GT.0) THEN
         IF(IFREFM.EQ.0) THEN
            READ(IN,'(2I10)') ITMP,NP
         ELSE
            READ(IN,*) ITMP,NP
         END IF
      ELSE
         NP=0
         IF(IFREFM.EQ.0) THEN
            READ(IN,'(I10)') ITMP
         ELSE
            READ(IN,*) ITMP
         END IF
      END IF
C
C------CALCULATE SOME CONSTANTS
      NAUX=NDRNVL-5-IDRNAL
      ITERPU = 1
      IOUTU = IOUT
      IF (NOPRDR.EQ.1) THEN
        ITERPU = 99
        IOUTU = -IOUT
      ENDIF
C
C2------DETERMINE THE NUMBER OF NON-PARAMETER DRAINS.
      IF(ITMP.LT.0) THEN
         WRITE(IOUT,7)
    7    FORMAT(1X,/1X,
     1        'REUSING NON-PARAMETER DRAINS FROM LAST STRESS PERIOD')
      ELSE
         NNPDRN=ITMP
      END IF
C
C3------IF THERE ARE NEW NON-PARAMETER DRAINS, READ THEM.
      MXACTD=IDRNPB-1
      IF(ITMP.GT.0) THEN
         IF(NNPDRN.GT.MXACTD) THEN
            WRITE(IOUT,99) NNPDRN,MXACTD
   99       FORMAT(1X,/1X,'THE NUMBER OF ACTIVE DRAINS (',I6,
     1                     ') IS GREATER THAN MXACTD(',I6,')')
            CALL USTOP(' ')
         END IF
         CALL ULSTRD(NNPDRN,DRAI,1,NDRNVL,MXDRN,IDRNAL,IN,IOUT,
     1     'DRAIN NO.  LAYER   ROW   COL     DRAIN EL.  CONDUCTANCE',
     2     DRNAUX,5,NAUX,IFREFM,NCOL,NROW,NLAY,5,5,ITERPU)
      END IF
      NDRAIN=NNPDRN
C
C1C-----IF THERE ARE ACTIVE DRN PARAMETERS, READ THEM AND SUBSTITUTE
      CALL PRESET('DRN')
      IF(NP.GT.0) THEN
         NREAD=NDRNVL-IDRNAL
         DO 30 N=1,NP
         CALL UPARLSTSUB(IN,'DRN',IOUTU,'DRN',DRAI,NDRNVL,MXDRN,NREAD,
     1                MXACTD,NDRAIN,5,5,
     2     'DRAIN NO.  LAYER   ROW   COL     DRAIN EL.  CONDUCTANCE',
     3            DRNAUX,5,NAUX)
   30    CONTINUE
      END IF
C
C3------PRINT NUMBER OF DRAINS IN CURRENT STRESS PERIOD.
      WRITE (IOUT,101) NDRAIN
  101 FORMAT(1X,/1X,I6,' DRAINS')
C
C8------RETURN.
  260 RETURN
      END
      SUBROUTINE GWF1DRN6FM(NDRAIN,MXDRN,DRAI,HNEW,HCOF,RHS,IBOUND,
     1              NCOL,NROW,NLAY,NDRNVL)
C
C-----VERSION 11JAN2000 DRN5FM
C     ******************************************************************
C     ADD DRAIN FLOW TO SOURCE TERM
C     ******************************************************************
C
C        SPECIFICATIONS:
C     ------------------------------------------------------------------
      DOUBLE PRECISION HNEW,EEL
C
      DIMENSION DRAI(NDRNVL,MXDRN),HNEW(NCOL,NROW,NLAY),
     1         RHS(NCOL,NROW,NLAY),IBOUND(NCOL,NROW,NLAY),
     1         HCOF(NCOL,NROW,NLAY)
C     ------------------------------------------------------------------
C
C1------IF NDRAIN<=0 THERE ARE NO DRAINS. RETURN.
      IF(NDRAIN.LE.0) RETURN
C
C2------PROCESS EACH CELL IN THE DRAIN LIST.
      DO 100 L=1,NDRAIN
C
C3------GET COLUMN, ROW AND LAYER OF CELL CONTAINING DRAIN.
      IL=DRAI(1,L)
      IR=DRAI(2,L)
      IC=DRAI(3,L)
C
C4-------IF THE CELL IS EXTERNAL SKIP IT.
      IF(IBOUND(IC,IR,IL).LE.0) GO TO 100
C
C5-------IF THE CELL IS INTERNAL GET THE DRAIN DATA.
      EL=DRAI(4,L)
      EEL=EL
C
C6------IF HEAD IS LOWER THAN DRAIN THEN SKIP THIS CELL.
      IF(HNEW(IC,IR,IL).LE.EEL) GO TO 100
C
C7------HEAD IS HIGHER THAN DRAIN. ADD TERMS TO RHS AND HCOF.
      C=DRAI(5,L)
      HCOF(IC,IR,IL)=HCOF(IC,IR,IL)-C
      RHS(IC,IR,IL)=RHS(IC,IR,IL)-C*EL
  100 CONTINUE
C
C8------RETURN.
      RETURN
      END
      SUBROUTINE GWF1DRN6BD(NDRAIN,MXDRN,VBNM,VBVL,MSUM,DRAI,DELT,HNEW,
     1        NCOL,NROW,NLAY,IBOUND,KSTP,KPER,IDRNCB,ICBCFL,BUFF,IOUT,
     2        PERTIM,TOTIM,NDRNVL,IDRNAL,IAUXSV)
C-----VERSION 05JUNE2000 GWF1DRN6BD
C     ******************************************************************
C     CALCULATE VOLUMETRIC BUDGET FOR DRAINS
C     ******************************************************************
C
C        SPECIFICATIONS:
C     ------------------------------------------------------------------
      COMMON /DRNCOM/DRNAUX(5)
      CHARACTER*16 DRNAUX
      CHARACTER*16 VBNM(MSUM),TEXT
      DOUBLE PRECISION HNEW,HHNEW,EEL,CC,CEL,RATOUT,QQ
C
      DIMENSION VBVL(4,MSUM),DRAI(NDRNVL,MXDRN),HNEW(NCOL,NROW,NLAY),
     1          IBOUND(NCOL,NROW,NLAY),BUFF(NCOL,NROW,NLAY)
C
      DATA TEXT /'          DRAINS'/
C     ------------------------------------------------------------------
C
C1------INITIALIZE CELL-BY-CELL FLOW TERM FLAG (IBD) AND
C1------ACCUMULATOR (RATOUT).
      ZERO=0.
      RATOUT=ZERO
      IBD=0
      IF(IDRNCB.LT.0 .AND. ICBCFL.NE.0) IBD=-1
      IF(IDRNCB.GT.0) IBD=ICBCFL
      IBDLBL=0
C
C2------IF CELL-BY-CELL FLOWS WILL BE SAVED AS A LIST, WRITE HEADER.
      IF(IBD.EQ.2) THEN
         NAUX=NDRNVL-5-IDRNAL
         IF(IAUXSV.EQ.0) NAUX=0
         CALL UBDSV4(KSTP,KPER,TEXT,NAUX,DRNAUX,IDRNCB,NCOL,NROW,NLAY,
     1          NDRAIN,IOUT,DELT,PERTIM,TOTIM,IBOUND)
      END IF
C
C3------CLEAR THE BUFFER.
      DO 50 IL=1,NLAY
      DO 50 IR=1,NROW
      DO 50 IC=1,NCOL
      BUFF(IC,IR,IL)=ZERO
50    CONTINUE
C
C4------IF THERE ARE NO DRAINS THEN DO NOT ACCUMULATE DRAIN FLOW.
      IF(NDRAIN.LE.0) GO TO 200
C
C5------LOOP THROUGH EACH DRAIN CALCULATING FLOW.
      DO 100 L=1,NDRAIN
C
C5A-----GET LAYER, ROW & COLUMN OF CELL CONTAINING REACH.
      IL=DRAI(1,L)
      IR=DRAI(2,L)
      IC=DRAI(3,L)
      Q=ZERO
C
C5B-----IF CELL IS NO-FLOW OR CONSTANT-HEAD, IGNORE IT.
      IF(IBOUND(IC,IR,IL).LE.0) GO TO 99
C
C5C-----GET DRAIN PARAMETERS FROM DRAIN LIST.
      EL=DRAI(4,L)
      EEL=EL
      C=DRAI(5,L)
      HHNEW=HNEW(IC,IR,IL)
C
C5D-----IF HEAD HIGHER THAN DRAIN, CALCULATE Q=C*(EL-HHNEW).
C5D-----SUBTRACT Q FROM RATOUT.
      IF(HHNEW.GT.EEL) THEN
         CC=C
         CEL=C*EL
         QQ=CEL - CC*HHNEW
         Q=QQ
         RATOUT=RATOUT-QQ
      END IF
C
C5E-----PRINT THE INDIVIDUAL RATES IF REQUESTED(IDRNCB<0).
      IF(IBD.LT.0) THEN
         IF(IBDLBL.EQ.0) WRITE(IOUT,61) TEXT,KPER,KSTP
   61    FORMAT(1X,/1X,A,'   PERIOD ',I4,'   STEP ',I3)
         WRITE(IOUT,62) L,IL,IR,IC,Q
   62    FORMAT(1X,'DRAIN ',I6,'   LAYER ',I3,'   ROW ',I5,'   COL ',I5,
     1       '   RATE ',1PG15.6)
         IBDLBL=1
      END IF
C
C5F-----ADD Q TO BUFFER.
      BUFF(IC,IR,IL)=BUFF(IC,IR,IL)+Q
C
C5G-----IF SAVING CELL-BY-CELL FLOWS IN A LIST, WRITE FLOW.  OR IF
C5G-----RETURNING THE FLOW IN THE DRAI ARRAY, COPY FLOW TO DRAI.
   99 IF(IBD.EQ.2) CALL UBDSVB(IDRNCB,NCOL,NROW,IC,IR,IL,Q,
     1                  DRAI(1,L),NDRNVL,NAUX,6,IBOUND,NLAY)
      IF(IDRNAL.NE.0) DRAI(NDRNVL,L)=Q
  100 CONTINUE
C
C6------IF CELL-BY-CELL FLOW WILL BE SAVED AS A 3-D ARRAY,
C6------CALL UBUDSV TO SAVE THEM.
      IF(IBD.EQ.1) CALL UBUDSV(KSTP,KPER,TEXT,IDRNCB,BUFF,NCOL,NROW,
     1                          NLAY,IOUT)
C
C7------MOVE RATES,VOLUMES & LABELS INTO ARRAYS FOR PRINTING.
  200 ROUT=RATOUT
      VBVL(3,MSUM)=ZERO
      VBVL(4,MSUM)=ROUT
      VBVL(2,MSUM)=VBVL(2,MSUM)+ROUT*DELT
      VBNM(MSUM)=TEXT
C
C8------INCREMENT BUDGET TERM COUNTER.
      MSUM=MSUM+1
C
C9------RETURN.
      RETURN
      END

ext: f

max_stars_repo_name: quanpands/wflow

max_issues_repo_name: quanpands/wflow

max_forks_repo_name: quanpands/wflow


text:
Address(Pamplona Avenue) is a residential double ended Culdesacs culdesac in North Davis.

 Intersecting Streets

  Almeria Street
  Cadiz Street

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
C
C
*+ mapcat_enqsr

       subroutine mapcat_enqsr(imap,source,program,status)
C      ---------------------------------------------------
C
C Enquire the source name for catalogue entry IMAP
C
C Input:
C    Map catalogue entry
       integer             imap
C Returned:
C    Source name for entry IMAP
       character*(*)       source
C    Program name
       character*(*)       program
C    Status
       integer             status
*-
       include 'mapcat_cat.inc'
       if (status.ne.0) return
       call mapcat_chk(imap,'NONE',status)
       call mapcat_open(status)
       call mapcat_read(imap,status)
       source  = current_source
       program = current_program
       call mapcat_err(status,'mapcat_enqsr',' ')
       end

max_stars_repo_name: CavendishAstrophysics/anmap

max_issues_repo_name: CavendishAstrophysics/anmap

max_forks_repo_name: CavendishAstrophysics/anmap


text:
      SUBROUTINE DBAR
C
C     THIS ROUTINE COMPUTES THE 2 6X6 DIFFERENTIAL STIFFNESS MATRICES
C     K(NPVT,NPVT) AND K(NPVT,J) FOR A BEAM HAVING END POINTS OF SIL
C     NOS. NPVT AND J.
C
C     ECPT FOR THE BEAM
C
C     ECPT( 1)  -  IELID          ELEMENT ID. NO.
C     ECPT( 2)  -  ISILNO(2)      SCALAR INDEX NOS.
C     ECPT( 3)  -   ...
C     ECPT( 4)  -  SMALLV(3)      REFERENCE VECTOR
C     ECPT( 5)  -   ...
C     ECPT( 6)  -   ...
C     ECPT( 7)  -  IGSUB0         OPTION FOR DEFINING REFERENCE NUMBER.
C     ECPT( 8)  -  IPINFL(2)      PIN FLAGS
C     ECPT( 9)  -   ...
C     ECPT(10)  -  ZA(3)          OFFSET VECTOR AT POINT A
C     ECPT(11)  -   ...
C     ECPT(12)  -   ...
C     ECPT(13)  -  ZB(3)          OFFSET VECTOR AT POINT B
C     ECPT(14)  -   ...
C     ECPT(15)  -   ...
C     ECPT(16)  -  GEF(4)         ECCENTRICITIES FOR FORCE
C     ECPT(17)  -   ...
C     ECPT(18)  -   ...
C     ECPT(19)  -   ...
C     ECPT(20)  -  IMATID         MATERIAL ID.
C     ECPT(21)  -  A              CROSS-SECTIONAL AREA
C     ECPT(22)  -  C1             STRESS COEFFICIENTS
C     ECPT(23)  -  C2                    ...
C     ECPT(24)  -  I1             AREA MOMENTS OF INERTIA
C     ECPT(25)  -  I2                    ...
C     ECPT(26)  -  I3                    ...
C     ECPT(27)  -  FJ             TORSIONAL CONSTANT
C     ECPT(28)  -  FMU            NON-STRUCTURAL MASS
C     ECPT(29)  -  K1             AREA FACTORS FOR SHEAR
C     ECPT(30)  -  K2                    ...
C     ECPT(31)  -  C3 (D1)        STRESS COEFFICIENTS
C     ECPT(32)  -  C4 (D2)               ...
C     ECPT(33)  -  B1             WIDTHS FOR FORCE
C     ECPT(34)  -  B2                    ...
C     ECPT(35)  -  HS1            DEPTHS FOR FORCE
C     ECPT(36)  -  HS2                   ...
C     ECPT(37)  -  HT1                   ...
C     ECPT(38)  -  HT2                   ...
C     ECPT(39)  -  MCSIDA         COOR. SYS. ID. FOR GRID PT. A
C     ECPT(40)  -  GPA(3)         BASIC COORDINATES FOR PT. A
C     ECPT(41)  -   ...                  ...
C     ECPT(42)  -   ...                  ...
C     ECPT(43)  -  MCSIDB         COOR. SYS. ID. FOR GRID PT. B
C     ECPT(44)  -  GPB(3)         BASIC COORDINATES FOR PT. B
C     ECPT(45)  -   ...                  ...
C     ECPT(46)  -   ...                  ...
C     ECPT(47)  -  ELTEMP         ELEMENT TEMPERATURE
C     ECPT(48)  -  ELDEF          ELEMENT DEFORMATION
C     ECPT(49)  -  TEMPER         ELEMENT LOADING TEMPERATURE
C     ECPT(50)  -  UAS(1)                ...
C     ECPT(51)  -  UAS(2)                ...
C     ECPT(52)  -  UAS(3)         SINGLE PRECISION DISPLACEMENTS
C     ECPT(53)  -  UAS(4)               FOR GRID POINT A
C     ECPT(54)  -  UAS(5)                ...
C     ECPT(55)  -  UAS(6)                ...
C     ECPT(56)  -  UBS(1)                ...
C     ECPT(57)  -  UBS(2)                ...
C     ECPT(58)  -  UBS(3)         SINGLE PRECISION DISPLACEMENTS
C     ECPT(59)  -  UBS(4)               FOR GRID POINT B
C     ECPT(60)  -  UBS(5)                ...
C     ECPT(61)  -  UBS(6)                ...
C
      LOGICAL          ABASIC,BBASIC,BASIC,AOFSET,BOFSET,OFFSET
      REAL             K1,K2,I1,I2
      DOUBLE PRECISION TA(18),TB(9),SMALV0(6),DELA,DELB,KE,KEP,VECI,
     1                 VECJ,VECK,FL,FLL,EI1,EI2,GAK1,GAK2,RRV1,RRV2,
     2                 SK1,SK2,SK3,SK4,TERM1,TERM2,TERM3,TERM4,L,LSQ,
     3                 LCUBE,DP(8)
      DOUBLE PRECISION E,DA,ALPHA,T SUB 0,SA(72),SB(36),UA(6),UB(6),
     1                 DPVECA(6),DPVECB(6),FX,VY,VZ,MAY,MAZ,MBY,MBZ,
     2                 KD(144),KC(12,12),TERM5,TERM6,TERM7,TERM8,TERM9,
     3                 TERM10,TERM11,KES(144),KDP(144),DFJ
      DIMENSION        VECI(3),VECJ(3),VECK(3),ECPT(100),IECPT(100),
     1                 IPIN(10),IZ(1)
      COMMON /ZZZZZZ/  Z(1)
      COMMON /DS1AAA/  NPVT,ICSTM,NCSTM,DUMCL(32),NOGO
      COMMON /DS1AET/  IELID,ISILNO(2),SMALLV(3),IGSUB0,IPINFL(2),ZA(3),
     1                 ZB(3),GEF(4),IMATID,A,DUMMY1,DUMMY2,I1,I2,DUMMY3,
     2                 FJ,FMU,K1,K2,DUM2(8),MCSIDA,GPA(3),MCSIDB,GPB(3),
     3                 TEMPEL,ELDEF,TEMPER,UAS(6),UBS(6),DUM3(38)
      COMMON /DS1ADP/  KE(144),KEP(144),DELA(6),DELB(6)
      COMMON /MATIN /  MATIDC,MATFLG,ELTEMP,STRESS,SINTH,COSTH
      COMMON /MATOUT/  E S,G S,NU,RHO,ALPHA S,T SUB 0 S,GSUBE,SIGT,
     1                 SIGC,SIGS
      EQUIVALENCE      (IELID,ECPT(1),IECPT(1)),(IZ(1),Z(1)),
     1                 (TA(10),TB(1)),(ECPT(71),DP(1)),(KC(1,1),KD(1)),
     2                 (SA(37),SB(1))
C
C     DETERMINE WHICH SIL IS THE PIVOT POINT.
C
C     IPVT = 0
      IPVT = 1
      IF (ISILNO(1) .EQ. NPVT) GO TO 20
      IPVT = 2
      IF (ISILNO(2) .NE. NPVT) CALL MESAGE (-30,34,IECPT(1))
C
C     JCSIDA IS AN INDEX WHICH POINTS TO THE COOR. SYS. ID. OF POINT A.
C     JOFSTA IS AN INDEX WHICH POINTS TO THE OFFSET VECTOR FOR POINT A.
C     SIMILARY FOR JCSIDB AND JOFSTB AND POINT B.
C
   20 JCSIDA = 39
      JCSIDB = 43
      JOFSTA = 10
      JOFSTB = 13
      JPINA  = 8
      JPINB  = 9
      ICSIDA = IECPT(JCSIDA)
      ICSIDB = IECPT(JCSIDB)
C
C     NORMALIZE THE REFERENCE VECTOR WHICH LIES IN THE FIRST PRINCIPAL
C     AXIS PLANE  (FMMS - 36 P. 4)
C     WE STORE SMALLV IN SMALV0 SO THAT ARITHMETIC WILL BE DOUBLE
C     PRECISION
C
      DO 50 I = 1,3
   50 SMALV0(I) = SMALLV(I)
      FL = DSQRT(SMALV0(1)**2 + SMALV0(2)**2 + SMALV0(3)**2)
      IF (FL .LE. 0.0D0) GO TO 700
      DO 60 I = 1,3
   60 SMALV0(I) = SMALV0(I)/FL
C
C     DETERMINE IF POINT A AND B ARE IN BASIC COORDINATES OR NOT.
C
      ABASIC = .TRUE.
      BBASIC = .TRUE.
      IF (ICSIDA .NE. 0) ABASIC = .FALSE.
      IF (ICSIDB .NE. 0) BBASIC = .FALSE.
C
C     COMPUTE THE TRANSFORMATION MATRICES TA AND TB IF NECESSARY
C
      IF (.NOT.ABASIC) CALL TRANSD (ECPT(JCSIDA),TA)
      IF (.NOT.BBASIC) CALL TRANSD (ECPT(JCSIDB),TB)
C
C     DETERMINE IF WE HAVE NON-ZERO OFFSET VECTORS.
C
      AOFSET = .TRUE.
      J = JOFSTA - 1
      DO 70 I = 1,3
      J = J + 1
      IF (ECPT(J) .NE. 0.0) GO TO 80
   70 CONTINUE
      AOFSET = .FALSE.
   80 BOFSET = .TRUE.
      J = JOFSTB - 1
      DO 90 I = 1,3
      J = J + 1
      IF (ECPT(J) .NE. 0.0) GO TO 100
   90 CONTINUE
      BOFSET = .FALSE.
C
C     FORM THE CENTER AXIS OF THE BEAM WITHOUT OFFSETS.
C     FIRST WE STORE THE COORDINATES IN THE ARRAY DP SO THAT ALL
C     ARITHMETIC WILL BE DOUBLE PRECISION.
C
  100 DP(1)   = ECPT(JCSIDA+1)
      DP(2)   = ECPT(JCSIDA+2)
      DP(3)   = ECPT(JCSIDA+3)
      DP(4)   = ECPT(JCSIDB+1)
      DP(5)   = ECPT(JCSIDB+2)
      DP(6)   = ECPT(JCSIDB+3)
      VECI(1) = DP(1) - DP(4)
      VECI(2) = DP(2) - DP(5)
      VECI(3) = DP(3) - DP(6)
C
C     TRANSFORM THE OFFSET VECTORS IF NECESSARY
C
      IF (.NOT.AOFSET .AND. .NOT.BOFSET) GO TO 150
C
C     TRANSFORM THE OFFSET VECTOR FOR POINT A IF NECESSARY.
C
      IDELA = 1
      J = JOFSTA - 1
      DO 110 I = 1,3
      J = J + 1
  110 DELA(I) = ECPT(J)
      IF (ABASIC) GO TO 120
      IDELA = 4
      CALL GMMATD (TA,3,3,0, DELA(1),3,1,0, DELA(4))
C
C     TRANSFORM THE OFFSET VECTOR FOR POINT B IF NECESSARY
C
  120 IDELB = 1
      J = JOFSTB - 1
      DO 130 I = 1,3
      J = J + 1
  130 DELB(I) = ECPT(J)
      IF (BBASIC) GO TO 140
      IDELB = 4
      CALL GMMATD (TB,3,3,0, DELB(1),3,1,0, DELB(4))
C
C     SINCE THERE WAS AT LEAST ONE NON-ZERO OFFSET VECTOR RECOMPUTE VECI
C
  140 VECI(1) = VECI(1) + DELA(IDELA  ) - DELB(IDELB  )
      VECI(2) = VECI(2) + DELA(IDELA+1) - DELB(IDELB+1)
      VECI(3) = VECI(3) + DELA(IDELA+2) - DELB(IDELB+2)
C
C     COMPUTE THE LENGTH OF THE BIG V (VECI) VECTOR AND NORMALIZE
C
  150 FL = DSQRT(VECI(1)**2 + VECI(2)**2 + VECI(3)**2)
      IF (FL .EQ. 0.0D0) GO TO 700
      DO 160 I = 1,3
  160 VECI(I) = VECI(I)/FL
C
C     COMPUTE THE SMALL V SUB 0 VECTOR, SMALV0.  ***CHECK THIS LOGIC***
C
      ITA = 1
      ISV = 1
      IF (MCSIDA.EQ.0 .OR. IGSUB0.EQ.0) GO TO 180
      IF (JCSIDA .NE. 39) ITA = 10
      ISV = 4
      CALL GMMATD (TA(ITA),3,3,0, SMALV0(1),3,1,0, SMALV0(4))
C
C     COMPUTE THE K VECTOR, VECK = VECI X SMALV0, AND NORMALIZE
C
  180 VECK(1) =  VECI(2)*SMALV0(ISV+2) - VECI(3)*SMALV0(ISV+1)
      VECK(2) =  VECI(3)*SMALV0(ISV  ) - VECI(1)*SMALV0(ISV+2)
      VECK(3) =  VECI(1)*SMALV0(ISV+1) - VECI(2)*SMALV0(ISV  )
      FLL = DSQRT(VECK(1)**2 + VECK(2)**2 + VECK(3)**2)
      IF (FLL .EQ. 0.0D0) GO TO 700
      VECK(1) =  VECK(1)/FLL
      VECK(2) =  VECK(2)/FLL
      VECK(3) =  VECK(3)/FLL
C
C     COMPUTE THE J VECTOR, VECJ = VECK X VECI, AND NORMALIZE
C
      VECJ(1) =  VECK(2)*VECI(3) - VECK(3)*VECI(2)
      VECJ(2) =  VECK(3)*VECI(1) - VECK(1)*VECI(3)
      VECJ(3) =  VECK(1)*VECI(2) - VECK(2)*VECI(1)
      FLL = DSQRT(VECJ(1)**2 + VECJ(2)**2 + VECJ(3)**2)
      IF (FLL .EQ. 0.0D0) GO TO 700
      VECJ(1) =  VECJ(1)/FLL
      VECJ(2) =  VECJ(2)/FLL
      VECJ(3) =  VECJ(3)/FLL
C
C     SEARCH THE MATERIAL PROPERTIES TABLE FOR E,G AND THE DAMPING
C     CONSTANT.
C
      MATIDC = IMATID
      MATFLG = 1
      ELTEMP = TEMPEL
      CALL MAT (IECPT(1))
C
C     COMPUTE THE RECIPROCALS OF RV1 AND RV2 (CALLING THEM RRV1 AND
C     RRV2)
C
      L = FL
      LSQ = L**2
      LCUBE = LSQ*L
C
C     STORE ECPT AND MPT VARIABLES IN DOUBLE PRECISION LOCATIONS.
C
      DP(1) = E S
      DP(2) = G S
      DP(3) = I1
      DP(4) = I2
      DP(5) = A
      EI1   = DP(1)*DP(3)
      EI2   = DP(1)*DP(4)
      IF (K1 .EQ. 0.0) GO TO 210
      DP(6) = K1
      GAK1  = DP(2)*DP(5)*DP(6)
      RRV1  = (12.0D0*EI1*GAK1)/(GAK1*LCUBE+12.0D0*L*EI1)
      GO TO 220
  210 RRV1  = 12.0D0*EI1/LCUBE
  220 IF (K2 .EQ. 0.0) GO TO 230
      DP(7) = K2
      GAK2  = DP(2)*DP(5)*DP(7)
      RRV2  = (12.0D0*EI2*GAK2)/(GAK2*LCUBE+12.0D0*L*EI2)
      GO TO 240
  230 RRV2  = 12.0D0*EI2/LCUBE
C
C     COMPUTE THE -SMALL- K-S, SK1, SK2, SK3 AND SK4
C
  240 SK1 = 0.25D0*RRV1*LSQ + EI1/L
      SK2 = 0.25D0*RRV2*LSQ + EI2/L
      SK3 = 0.25D0*RRV1*LSQ - EI1/L
      SK4 = 0.25D0*RRV2*LSQ - EI2/L
C
C     COMPUTE THE TERMS THAT WILL BE NEEDED FOR THE 12 X 12 MATRIX KE
C
      TERM1 = DP(5)*DP(1)/L
      TERM2 = 0.5D0*L*RRV1
      TERM3 = 0.5D0*L*RRV2
      DP(8) = FJ
      TERM4 = DP(2)*DP(8)/L
C
C     CONSTRUCT THE 12 X 12 MATRIX KE
C
      DO 250 I = 1,144
  250 KE(  I) =  0.0D0
      KE(  1) =  TERM1
      KE(  7) = -TERM1
      KE( 14) =  RRV1
      KE( 18) = -TERM2
      KE( 20) = -RRV1
      KE( 24) = -TERM2
      KE( 27) =  RRV2
      KE( 29) =  TERM3
      KE( 33) = -RRV2
      KE( 35) =  TERM3
      KE( 40) =  TERM4
      KE( 46) = -TERM4
      KE( 51) =  TERM3
      KE( 53) =  SK2
      KE( 57) = -TERM3
      KE( 59) =  SK4
      KE( 62) = -TERM2
      KE( 66) =  SK1
      KE( 68) =  TERM2
      KE( 72) =  SK3
      KE( 73) = -TERM1
      KE( 79) =  TERM1
      KE( 86) = -RRV1
      KE( 90) =  TERM2
      KE( 92) =  RRV1
      KE( 96) =  TERM2
      KE( 99) = -RRV2
      KE(101) = -TERM3
      KE(105) =  RRV2
      KE(107) = -TERM3
      KE(112) = -TERM4
      KE(118) =  TERM4
      KE(123) =  TERM3
      KE(125) =  SK4
      KE(129) = -TERM3
      KE(131) =  SK2
      KE(134) = -TERM2
      KE(138) =  SK3
      KE(140) =  TERM2
      KE(144) =  SK1
C
C     DETERMINE IF THERE ARE NON-ZERO PIN FLAGS.
C
      KA = IECPT(JPINA)
      KB = IECPT(JPINB)
      IF (KA.EQ.0 .AND. KB.EQ.0) GO TO 325
C
C     SAVE THE KE (UNPINNED) MATRIX IN KES.
C
      DO 255 I = 1,144
  255 KES(I) = KE(I)
C
C     SET UP THE IPIN ARRAY
C
      DO 260 I = 1,5
      IPIN(I  ) = MOD(KA,10)
      IPIN(I+5) = MOD(KB,10) + 6
      IF (IPIN(I+5) .EQ. 6) IPIN(I+5) = 0
      KA = KA/10
  260 KB = KB/10
C
C     ALTER KE MATRIX DUE TO PIN FLAGS.
C
      DO 320 I = 1,10
      IF (IPIN(I) .EQ. 0) GO TO 320
      II = 13*IPIN(I) - 12
      IF (KE(II) .NE. 0.0D0) GO TO 280
      IL = IPIN(I)
      II = II - IL
      DO 270 J = 1,12
      II = II + 1
      KE(II) = 0.0D0
      KE(IL) = 0.0D0
      IL = IL + 12
  270 CONTINUE
      GO TO 320
  280 DO 300 J = 1,12
      JI  = 12*(J-1) + IPIN(I)
      IJ  = 12*(IPIN(I)-1) + J
      DO 290 LL = 1,12
      JLL = 12*(J-1) + LL
      ILL = 12*(IPIN(I)-1) + LL
      KEP(JLL) = KE(JLL) - (KE(ILL)/KE(II))*KE(JI)
  290 CONTINUE
      KEP(IJ) = 0.0D0
      KEP(JI) = 0.0D0
  300 CONTINUE
      DO 310 K = 1,144
  310 KE(K) = KEP(K)
  320 CONTINUE
C
C            E
C     STORE K   IN KEP(1),...,KEP(36) AND
C            AA
C
C            E
C     STORE K   IN KEP(37),...,KEP(72)
C            AB
C
  325 J = 0
      DO 340 I = 1,72,12
      LOW = I
      LIM = LOW + 5
      DO 330 K = LOW,LIM
      J = J + 1
      KEP(J   ) = KE(K  )
  330 KEP(J+36) = KE(K+6)
  340 CONTINUE
C
C                                                            T
C     STORE VECI, VECJ, VECK IN KE(1),...,KE(9) FORMING THE A  MATRIX.
C
      KE(1) = VECI(1)
      KE(2) = VECI(2)
      KE(3) = VECI(3)
      KE(4) = VECJ(1)
      KE(5) = VECJ(2)
      KE(6) = VECJ(3)
      KE(7) = VECK(1)
      KE(8) = VECK(2)
      KE(9) = VECK(3)
C
C     SET POINTERS SO THAT WE WILL BE WORKING WITH POINT A.
C
      BASIC  = ABASIC
      JCSID  = JCSIDA
      OFFSET = AOFSET
      JOFSET = JOFSTA
      IWBEG  = 0
      IKEL   = 1
      IAB    = 1
C
C     ZERO OUT THE ARRAY WHERE THE 3 X 3 MATRIX H AND THE W  AND W
C     6 X 6 MATRICES WILL RESIDE.                          A      B
C
      DO 360 I = 28,108
  360 KE(I) = 0.0D0
C
C     SET UP THE -G- MATRIX.  IG POINTS TO THE BEGINNING OF THE G MATRIX
C     G = AT X TI
C
  365 IG = 1
      IF (BASIC) GO TO 370
      CALL TRANSD (ECPT(JCSID),KE(10))
      CALL GMMATD (KE(1),3,3,0, KE(10),3,3,0, KE(19))
      IG = 19
C
C     IF THERE IS A NON-ZERO OFFSET FOR THE POINT, SET UP THE D 3 X 3
C     MATRIX.
C
  370 IF (.NOT.OFFSET) GO TO 380
      KE(10) =  0.0D0
      KE(11) =  ECPT(JOFSET+2)
      KE(12) = -ECPT(JOFSET+1)
      KE(13) = -KE(11)
      KE(14) =  0.0D0
      KE(15) =  ECPT(JOFSET)
      KE(16) = -KE(12)
      KE(17) = -KE(15)
      KE(18) =  0.0D0
C
C     FORM THE 3 X 3 PRODUCT H = G X D, I.E., KE(28) = KE(IG) X KE(10)
C
      CALL GMMATD (KE(IG),3,3,0, KE(10),3,3,0, KE(28))
C
C
C     FORM THE W  MATRIX OR THE W  MATRIX IN KE(37) OR KE(73) DEPENDING
C               A                B
C     UPON WHICH POINT - A OR B - IS UNDER CONSIDERATION.  G WILL BE
C     STORED IN THE UPPER LEFT AND LOWER RIGHT CORNERS.  H, IF NON-ZERO,
C     WILL BE STORED IN THE UPPER RIGHT CORNER.
C
  380 KE(IWBEG+37) = KE(IG  )
      KE(IWBEG+38) = KE(IG+1)
      KE(IWBEG+39) = KE(IG+2)
      KE(IWBEG+43) = KE(IG+3)
      KE(IWBEG+44) = KE(IG+4)
      KE(IWBEG+45) = KE(IG+5)
      KE(IWBEG+49) = KE(IG+6)
      KE(IWBEG+50) = KE(IG+7)
      KE(IWBEG+51) = KE(IG+8)
      KE(IWBEG+58) = KE(IG  )
      KE(IWBEG+59) = KE(IG+1)
      KE(IWBEG+60) = KE(IG+2)
      KE(IWBEG+64) = KE(IG+3)
      KE(IWBEG+65) = KE(IG+4)
      KE(IWBEG+66) = KE(IG+5)
      KE(IWBEG+70) = KE(IG+6)
      KE(IWBEG+71) = KE(IG+7)
      KE(IWBEG+72) = KE(IG+8)
      IF (.NOT.OFFSET) GO TO 390
      KE(IWBEG+40) = KE(28)
      KE(IWBEG+41) = KE(29)
      KE(IWBEG+42) = KE(30)
      KE(IWBEG+46) = KE(31)
      KE(IWBEG+47) = KE(32)
      KE(IWBEG+48) = KE(33)
      KE(IWBEG+52) = KE(34)
      KE(IWBEG+53) = KE(35)
      KE(IWBEG+54) = KE(36)
C
C                                 E                     E
C     COMPUTE THE PRODUCT S   =  K   X  W   OR  S   =  K    X  W
C                          A      AA     A       B      AB      B
C     WHERE
C                  T                           T
C           W  =  T   X  C  X  E   AND  W  =  T   X  C   X  E
C            A     EB     A     A        B     EB     B      B
C
C     W AT KE(37) AND W AT KE(73) WILL BE USED AGAIN BEFORE FINAL STEPS.
C      A               B
C
  390 CALL GMMATD (KEP(IKEL),6,6,0, KE(IWBEG+37),6,6,0, SA(IAB))
C
C     IF THE POINT UNDER CONSIDERATION IS POINT B WE ARE FINISHED. IF
C     NOT, SET UP POINTS AND INDICATORS FOR WORKING WITH POINT B.
C
      IF (IWBEG .EQ. 36) GO TO 500
      BASIC  = BBASIC
      JCSID  = JCSIDB
      OFFSET = BOFSET
      JOFSET = JOFSTB
      IWBEG  = 36
      IKEL   = 37
      IAB    = 37
      DO 400 I = 28,36
  400 KE(I)  = 0.0D0
      GO TO 365
C
C     BEGIN DIFFERENTIAL STIFFNESS PORTION OF THIS ROUTINE.
C
C     STORE DISPLACEMENT VECTORS IN DOUBLE PRECISION LOCATIONS
C
  500 DO 510 I = 1,6
      UA(I) = ECPT(I+49)
  510 UB(I) = ECPT(I+55)
C
C     COMPUTE  S  X  U   AND  S  X  U
C               A     A        B     B
C
      CALL GMMATD (SA(1),6,6,0, UA,6,1,0, DPVECA)
      CALL GMMATD (SB(1),6,6,0, UB,6,1,0, DPVECB)
C
C     COMPUTE THE NEEDED COMPONENTS OF THE FORCE VECTOR.
C
      FX  = DPVECA(1) + DPVECB(1)
      VY  = DPVECA(2) + DPVECB(2)
      VZ  = DPVECA(3) + DPVECB(3)
      MAY = DPVECA(5) + DPVECB(5)
      MAZ = DPVECA(6) + DPVECB(6)
      MBZ = -MAZ - VY*L
      MBY = -MAY + VZ*L
      E   = E S
      FX  = FX - E*ELDEF/L
      IF (IECPT(49) .EQ. -1) GO TO 520
      ALPHA = ALPHAS
      TSUB0 = TSUB0S
      DP(1) = TEMPER
      FX  = FX - A*ALPHA*E*(DP(1)-TSUB0)
C
C     ZERO OUT THE KD (KC) MATRIX
C
  520 DO 530 I = 1,144
  530 KD(I) = 0.0D0
C
C     FORM THE ELEMENT DIFFERENTIAL STIFFNESS MATRIX (UPPER HALF)
C
      TERM1  = 6.0D0*FX/(5.0D0*L)
      TERM2  = -MAY/L
      TERM3  = FX/10.0D0
      TERM4  = -MBY/L
      TERM5  = -MAZ/L
      TERM6  = -MBZ/L
      DFJ    = I1 + I2
      DA     = A
      TERM7  = DFJ*FX/(L*DA)
      TERM8  = L*VY/6.0D0
      TERM9  = L*VZ/6.0D0
      TERM10 = 2.0D0*L*FX/15.0D0
      TERM11 = L*FX/30.0D0
      KC( 2, 2) =  TERM1
      KC( 2, 4) =  TERM2
      KC( 2, 6) = -TERM3
      KC( 2, 8) = -TERM1
      KC( 2,10) =  TERM4
      KC( 2,12) = -TERM3
      KC( 3, 3) =  TERM1
      KC( 3, 4) =  TERM5
      KC( 3, 5) =  TERM3
      KC( 3, 9) = -TERM1
      KC( 3,10) =  TERM6
      KC( 3,11) =  TERM3
      KC( 4, 4) =  TERM7
      KC( 4, 5) = -TERM8
      KC( 4, 6) = -TERM9
      KC( 4, 8) = -TERM2
      KC( 4, 9) = -TERM5
      KC( 4,10) = -TERM7
      KC( 4,11) =  TERM8
      KC( 4,12) =  TERM9
      KC( 5, 5) =  TERM10
      KC( 5, 9) = -TERM3
      KC( 5,10) =  TERM8
      KC( 5,11) = -TERM11
      KC( 6, 6) =  TERM10
      KC( 6, 8) =  TERM3
      KC( 6,10) =  TERM9
      KC( 6,12) = -TERM11
      KC( 8, 8) =  TERM1
      KC( 8,10) = -TERM4
      KC( 8,12) =  TERM3
      KC( 9, 9) =  TERM1
      KC( 9,10) = -TERM6
      KC( 9,11) = -TERM3
      KC(10,10) =  TERM7
      KC(10,11) = -TERM8
      KC(10,12) = -TERM9
      KC(11,11) =  TERM10
      KC(12,12) =  TERM10
C
C     STORE THE UPPER HALF IN THE LOWER HALF.
C
      DO 550 I = 2,10
      LOW = I + 1
      DO 540 J = LOW,12
      KC(J,I) = KC(I,J)
  540 CONTINUE
  550 CONTINUE
C
C     IF THERE PIN FLAGS, ALTER THE KD MATRIX
C
      IF (KA.EQ.0 .AND. KB.EQ.0) GO TO 620
C
C     ALTER KD DUE TO PIN FLAGS.
C
      DO 610 J = 1,10
      IF (IPIN(J) .EQ. 0) GO TO 610
      JJ = 12*(IPIN(J)-1) + IPIN(J)
      IF (KES(JJ) .EQ. 0.0D0) GO TO 605
      DO 590 I = 1,12
      JI = 12*(IPIN(J)-1) + I
      IJ = 12*(I-1) + IPIN(J)
      DO 580 L1 = 1,12
      IL = 12*(I-1) + L1
      LJ = 12*(L1-1) + IPIN(J)
      KDP(IL) = KD(IL) - KES(LJ)*KD(JI)/KES(JJ) - KES(JI)*KD(LJ)/KES(JJ)
     2        + KES(LJ)* KES(JI)*KD(JJ)/KES(JJ)**2
  580 CONTINUE
  590 CONTINUE
      DO 600 KK = 1,144
  600 KD(KK) = KDP(KK)
C
C     ZERO OUT THE IPIN(J) TH ROW AND COLUMN OF KD.
C
  605 J1 = JJ - IPIN(J)
      J2 = IPIN(J)
      DO 608 KK = 1,12
      J1 = J1 + 1
      KD(J1) = 0.0D0
      KD(J2) = 0.0D0
  608 J2 = J2 + 12
  610 CONTINUE
C
C            D
C     STORE K        AT KEP(1),...,KEP(36)  AND
C            NPVT,A
C
C            D
C           K        AT KEP(37),...,KEP(72)
C            NPVT,B
C
C
  620 J = 0
      IF (IPVT .EQ. 2) GO TO 625
      ILOW = 1
      ILIM = 72
      GO TO 628
  625 ILOW = 73
      ILIM = 144
  628 DO 640 I = ILOW,ILIM,12
      LOW  = I
      LIM  = LOW +5
      DO 630 K = LOW,LIM
      J = J + 1
      KEP(J   ) = KD(K  )
  630 KEP(J+36) = KD(K+6)
  640 CONTINUE
C
C     COMPUTE THE FINAL 2 6X6 DIFFERENTIAL STIFFNESS MATRICES FOR THIS
C     BEAM.
C
      IWLEFT = 37
      IF (IPVT .EQ. 2) IWLEFT = 73
      I = 1
      IKDE = 1
      IWRGHT = 37
  650 CALL GMMATD (KE(IWLEFT),6,6,1, KEP(IKDE),6,6,0, KEP(73))
      CALL GMMATD (KEP(73),6,6,0, KE(IWRGHT),6,6,0, KEP(109))
      CALL DS1B (KEP(109),ISILNO(I))
      IF (I .EQ. 2) RETURN
      I = 2
      IKDE = 37
      IWRGHT = 73
      GO TO 650
C
C     FATAL ERROR
C
  700 CALL MESAGE (30,26,IECPT(1))
C
C     SET FLAG FOR FATAL ERROR WHILE ALLOWING ERROR MESSAGES TO
C     ACCUMULATE
C
      NOGO = 1
      RETURN
      END


text:
C-----------------------------------------------------------------------
C
C     COPYRIGHT (C) 1992  MASSACHUSETTS INSTITUTE OF TECHNOLOGY
C
C-----------------------------------------------------------------------
C
C     Main program 'sphere'
C
C-----------------------------------------------------------------------
C
C This program generates a geometric data file for a sphere for
C input to the FastLap driver: driver.c.  To modify the output
C for the FastLap driver: driver.f, use c2fdata.
C
C
C Boundary conditions are for a sphere translating in an
C infinite fluid for which the solution is known in closed form.
C
C     Written by T. Korsmeyer, March 1992
C
C-------------------------------------------------------------------------
      PROGRAM SPHERE
      IMPLICIT NONE
      INTEGER MAX_NODE, MAX_PANEL, NODE_TOTAL, PANEL_TOTAL
      PARAMETER (MAX_NODE = 16384)
      PARAMETER (MAX_PANEL = 16384)
      CHARACTER*72 HEAD
      INTEGER NUM_SPHERES
      DOUBLE PRECISION X(3,MAX_NODE)
      DOUBLE PRECISION XC(3,4)
      DOUBLE PRECISION XVER(4,3), XVP(4),YVP(4),XCT(3),DIR(3,3),
     $     SS(15),SID(4),VEL(3),VROT(4,3),DIAG(2)
      INTEGER IDIAG
      DOUBLE PRECISION POTENTIAL, DBYDNPOTEN
      INTEGER FSCON(4,MAX_PANEL), MSCON(4,MAX_PANEL), FSPAN, MSPAN
      INTEGER NFS, NMS, NWL
      INTEGER LFS(MAX_NODE), LMS(MAX_NODE), LWL(MAX_NODE)
      INTEGER INDEX_FS, INDEX_MS, M_FIRST, COUNT_FS, COUNT_MS,
     $     COUNT_WL, NA, NM, OFFSET, I, J, K, N
      DOUBLE PRECISION THETA, PHI, D_THETA, D_PHI, R, D_R, PI, RAD,
     $     CLOSENESS
      DOUBLE PRECISION  PHI_NORM, PHIN_NORM, MAX_PHI, MAX_PHIN
      INTEGER TYPE
      DOUBLE PRECISION DUM
      DATA DUM / 1.0D00 /
      PARAMETER(RAD = 1.0)
      DATA PI / 3.1415927 /
C
C Define output device number and name.
C
      CHARACTER*20 FASTIN
      DATA FASTIN /'fastlap.in'/
C
C Output heading, and query on output files.
C
      WRITE (6,935)
 935  FORMAT(1X,/,
     &' A new panelization will be generated for a sphere.'///)
 1    CONTINUE
C
C Query on panelization.
C
 2    CONTINUE
      WRITE (6,941)
 941  FORMAT(1X,//5X,71(1H-),//
     &  8X,' The number of panels is specified by two inputs:'//
     &  8X,' Number from pole to pole:   ',$)
      READ (5,*) NA
      NA = NA / 2
      WRITE (6,942)
 942  FORMAT(1X,/
     &  8X,' Number around any parallel of latitude:   ',$)
      READ (5,*) NM
      OPEN (UNIT=14,FILE=FASTIN,STATUS='UNKNOWN')
C
C Generate nodes on the upper and lower surfaces.  We do them both
C at once because they have common x,y coordinates, only their elevation
C is different. Bottom dead center and top dead center are common to
C every meridian,  so don't repeat them.
C
      NFS = (NA-1) * NM + 1
      NMS = NFS
      NWL = NM
      FSPAN = NM * NA
      MSPAN = NM * NA
      NODE_TOTAL = NFS + NWL + NMS
      PANEL_TOTAL = FSPAN + MSPAN
      IF (NODE_TOTAL .GT. MAX_NODE) THEN
         WRITE(6,'('' Node dimension of '',I5,
     $        '' is exceeded, try again.'')') MAX_NODE
         GO TO 2
      END IF
      IF (PANEL_TOTAL .GT. MAX_PANEL) THEN
         WRITE(6,'('' Panel dimension of '',I5,
     $        '' is exceeded, try again.'')') MAX_PANEL
         GO TO 2
      END IF
      D_R = RAD / NA
      D_THETA = PI / (2. * NA)
      D_PHI = (2.0 * PI) / NM
      M_FIRST = NFS + NWL + 1
      INDEX_FS = 1
      INDEX_MS = M_FIRST
      COUNT_FS = 1
      COUNT_MS = 1
      COUNT_WL = 1
      X(1,INDEX_FS) = 0.0
      X(2,INDEX_FS) = 0.0
      X(3,INDEX_FS) = RAD
      X(1,INDEX_MS) =  0.0
      X(2,INDEX_MS) =  0.0
      X(3,INDEX_MS) = -RAD
      LFS(COUNT_FS) = INDEX_FS
      LMS(COUNT_MS) = INDEX_MS
      COUNT_FS = COUNT_FS + 1
      COUNT_MS = COUNT_MS + 1
C
C Add these nodes to the node lists
C
C
C Loop along the parallels.
C
      DO J = 0, NM-1
         PHI = J * D_PHI
C
C Loop along the meridians, but don't do the equator.
C
         DO I = 1, NA-1
            R = I * D_R
            THETA = PI - (I * D_THETA)
            INDEX_FS = INDEX_FS + 1
            INDEX_MS = INDEX_MS + 1
            X(1,INDEX_MS) = RAD * SIN(THETA) * COS(PHI)
            X(2,INDEX_MS) = RAD * SIN(THETA) * SIN(PHI)
            X(3,INDEX_MS) = RAD * COS(THETA)
            X(1,INDEX_FS) =  X(1,INDEX_MS)
            X(2,INDEX_FS) =  X(2,INDEX_MS)
            X(3,INDEX_FS) = -X(3,INDEX_MS)
            LFS(COUNT_FS) = INDEX_FS
            LMS(COUNT_MS) = INDEX_MS
            COUNT_FS = COUNT_FS + 1
            COUNT_MS = COUNT_MS + 1
         END DO
C
C Now do the waterline.
C
         THETA = PI - (I * D_THETA)
         INDEX_FS = INDEX_FS + 1
         X(1,INDEX_FS) = RAD * SIN(THETA) * COS(PHI)
         X(2,INDEX_FS) = RAD * SIN(THETA) * SIN(PHI)
         X(3,INDEX_FS) = 0.0
         LWL(COUNT_WL) = INDEX_FS
         COUNT_WL = COUNT_WL + 1
      END DO
C
C Generate the connectivity of the panels.  Vertex numbering is
C clockwise from inside the fluid domain.
C
      INDEX_FS = 0
      INDEX_MS = 0
      DO J = 1, NM
C
C First do the top and bottom dead center panels which are special.
C
         INDEX_FS = INDEX_FS + 1
         INDEX_MS = INDEX_MS + 1
         FSCON(4,INDEX_FS) = 1
         FSCON(3,INDEX_FS) = 2 + NA * (J-1)
         FSCON(2,INDEX_FS) = 2 + NA * J
         FSCON(1,INDEX_FS) = 1
         MSCON(1,INDEX_MS) = M_FIRST
         MSCON(2,INDEX_MS) = M_FIRST + 1 + (NA-1) * (J-1)
         MSCON(3,INDEX_MS) = M_FIRST + 1 + (NA-1) * J
         MSCON(4,INDEX_MS) = M_FIRST
C
C Now do the rest EXCEPT for the panels sharing the equator.
C
         DO I = 2, NA-1
            INDEX_FS = INDEX_FS + 1
            INDEX_MS = INDEX_MS + 1
            FSCON(4,INDEX_FS) = NA * (J-1) + I
            FSCON(1,INDEX_FS) = NA * J + I
            FSCON(3,INDEX_FS) = FSCON(4,INDEX_MS) + 1
            FSCON(2,INDEX_FS) = FSCON(1,INDEX_MS) + 1
            MSCON(1,INDEX_MS) = (NA-1) * (J-1) + (M_FIRST-1) + I
            MSCON(4,INDEX_MS) = (NA-1) * J + (M_FIRST-1) + I
            MSCON(3,INDEX_MS) = MSCON(4,INDEX_MS) + 1
            MSCON(2,INDEX_MS) = MSCON(1,INDEX_MS) + 1
         END DO
C
C Now do the equator panels.
C
         INDEX_FS = INDEX_FS + 1
         INDEX_MS = INDEX_MS + 1
         FSCON(4,INDEX_FS) = NA * (J-1) + NA
         FSCON(1,INDEX_FS) = NA * J + NA
         FSCON(3,INDEX_FS) = FSCON(4,INDEX_MS) + 1
         FSCON(2,INDEX_FS) = FSCON(1,INDEX_MS) + 1
         MSCON(1,INDEX_MS) = (NA-1) * (J-1) + (M_FIRST-1) + NA
         MSCON(4,INDEX_MS) = (NA-1) * J + (M_FIRST-1) + NA
         MSCON(2,INDEX_MS) = 1 + J * NA
         MSCON(3,INDEX_MS) = 1 + (J+1) * NA
      END DO
C
C Now correct the last strip of panels on which the node numbers will exceed
C the max node number unless they are wrapped into the first strip of nodes.
C
      OFFSET = NM*NA-NA
      DO I = 1,NA
         FSCON(2,OFFSET+I) = FSCON(3,I)
         FSCON(1,OFFSET+I) = FSCON(4,I)
         MSCON(4,OFFSET+I) = MSCON(1,I)
         MSCON(3,OFFSET+I) = MSCON(2,I)
      END DO
      WRITE (6,903)
 903  FORMAT (1X,/,' Enter header line for gdf file: ',$)
      READ (5,902) HEAD
 902  FORMAT (A72)
      WRITE (14,905) HEAD
 905  FORMAT ('0  ',A72)
      WRITE(*,'('' What type on trailing half? 0=Dirichlet''
     $     '' 1=Neumann.'')')
      READ(5,*) type
      MAX_PHI = 0.0
      MAX_PHIN = 0.0
      PHI_NORM = 0.0
      PHIN_NORM = 0.0
      DO I = 1, FSPAN
         DO J = 1,4
            XC(1,J) = X(1,FSCON(J,I))
            XC(2,J) = X(2,FSCON(J,I))
            XC(3,J) = X(3,FSCON(J,I))
         END DO
c
c Find the centroid of the panel.
c
         CALL PANEL(XC,XCT,DIR,XVP,YVP,SS,SID)
         POTENTIAL = -0.5 * XCT(3)
         DBYDNPOTEN = -XCT(3)
         IF (type .eq. 1) THEN
C
C This is a panel on which PHI is the unknown.
C
            PHI_NORM = PHI_NORM + ABS(POTENTIAL)
            MAX_PHI = MAX(ABS(POTENTIAL),MAX_PHI)
         ELSE
C
C This is a panel on which PHIN is the unknown.
C
            PHIN_NORM = PHIN_NORM + ABS(DBYDNPOTEN)
            MAX_PHIN = MAX(ABS(DBYDNPOTEN),MAX_PHIN)
         END IF
         IF (FSCON(1,I) .EQ. FSCON(4,I)) THEN
            WRITE(14,951) ((XC(K,J),K=1,3),J=1,3),
     $           POTENTIAL,DBYDNPOTEN,TYPE
         ELSE
            WRITE(14,950) ((XC(K,J),K=1,3),J=1,4),
     $           POTENTIAL,DBYDNPOTEN,TYPE
         END IF
      END DO
      write(*,'('' What type on leading half? 0=Dirichlet''
     $     '' 1=Neumann.'')')
      read(5,*) type
      DO I = 1, MSPAN
         DO J = 1,4
            XC(1,J) = X(1,MSCON(J,I))
            XC(2,J) = X(2,MSCON(J,I))
            XC(3,J) = X(3,MSCON(J,I))
         END DO
c
c Find the centroid of the panel.
c
         CALL PANEL(XC,XCT,DIR,XVP,YVP,SS,SID)
         POTENTIAL = -0.5 * XCT(3)
         DBYDNPOTEN = -XCT(3)
         IF (MSCON(1,I) .EQ. MSCON(4,I)) THEN
            WRITE(14,951) ((XC(K,J),K=1,3),J=1,3),
     $           POTENTIAL,DBYDNPOTEN,TYPE
         ELSE
            WRITE(14,950) ((XC(K,J),K=1,3),J=1,4),
     $           POTENTIAL,DBYDNPOTEN,TYPE
         END IF
      END DO
      WRITE(*,'(''Number of panels in fastlap.in:'',I5)')
     $     FSPAN + MSPAN
 950  FORMAT('Q  ',14f18.14,i3)
 951  FORMAT('T  ',11f18.14,i3)
 952  FORMAT(12E12.4)
      STOP
      END
C
C
C
      SUBROUTINE PANEL(V,XCTP,DIRP,XVP,YVP,S,SIDE)
C-----------------------------------------------------------------------
C
C     COPYRIGHT (C) 1988,1991 MASSACHUSETTS INSTITUTE OF TECHNOLOGY
C
C-----------------------------------------------------------------------
C
C     DESCRIPTION : Subroutine PANEL evaluates geometrical data for a
C     plane quadrilateral approximating in a mean sense the curvilinear
C     surface defined by four input vertices (not lying on a plane).
C     Adjacent vertices are connected by straight segments. The plane of
C     the quadrilateral is defined by the segment midpoints and its
C     vertices by the projection of the original vetrices on that plane.
C
C-----------------------------------------------------------------------
C
C     Input arguments (in alphabetical order) :
C
C     Name              Description
C     ------------------------------------------------------------------
C     V(I,J)            Body-system coordinates (x,y,z)<=>J=1,2,3
C                       of I-th vertex, I=1,..,4 of current panel
C
C     Output Arguments (in alphabetical order) :
C
C     Symbolic       Description
C     name
C     ------------------------------------------------------------------
C     DIRP(I,J)      Direction cosines of axes of panel coordinate
C                      system for current panel, relative to axes of
C                      body coordinate system
C                      I=1,2,3 <=> (x,y,z) axes of panel system,
C                      J=1,2,3 <=> (x,y,z) axes of body system
C     S(K)           Moments of inertia of current quadrilateral
C                      surface w.r.t. panel coordinate system, K=1,..,15
C                      Defined as in subroutine GEOM
C     XVP(I),YVP(I)    Panel coordinates of I-th vertex of current
C                      quadrilateral, I=1,..,4
C     SIDE(I)        Length of four sides (I=1,..,4) of current
C                      quadrilateral numbered couter-clockwise and
C                      defined so that SIDE(1) connects vertex 1 to
C                      vertex 2
C     XCTP(I)        Body-system coordinates of centroid of N-th
C                      quadrilateral, (x,y,z) <=> I=1,2,3
C
C-----------------------------------------------------------------------
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION V(3,4),XCTP(*),XVP(*),YVP(*),DIRP(3,*),S(*),SIDE(*)
C-----------------------------------------------------------------------
C     Allocate space for local arrays
C-----------------------------------------------------------------------
      DIMENSION XORP(3),UL(3),VL(3),N1(4),DX(4),DY(4),
     &          IX(4),IY(4),XP(6,4),YP(6,4),CS(15),K1(5)
      DATA  ONE /1.E+00/, HALF / 0.5E+00/, QUART / 0.25E+00/,
     &      ZERO/ 0.0E+0/, GND/ 0.5773503E+00 /,
     &      TOLS / 1.E-10/, TOLT / 1.E-12 /, THIRD / 0.3333333E+00/,
     &      CS / 1., 1., 1.5, 1.5, 3.75, 1., 3.0, 1.5, 7.5, 1.5, 1.5,
     &           3.75, 1.5, 7.5, 3.75 /,
     &      IX / 1, 1, -1, -1/, IY / 1, -1, -1, 1/,
     &      K1 / 1, 6, 10, 13, 15 /,  N1 / 2, 3, 4, 1 /
      DO 10 J=1,3
        XORP(J)=((V(J,1)+V(J,2))+(V(J,3)+V(J,4)))*QUART
        UL(J)=(V(J,1)+V(J,2))-(V(J,3)+V(J,4))
        VL(J)=(V(J,1)+V(J,4))-(V(J,2)+V(J,3))
10    CONTINUE
C-----------------------------------------------------------------------
C     Evaluation of direction cosines DIRP(3,3) of local coordinate
C     system
C-----------------------------------------------------------------------
      ULSQ=SQRT(UL(1)*UL(1)+UL(2)*UL(2)+UL(3)*UL(3))
      DIRP(1,1)=UL(1)/ULSQ
      DIRP(1,2)=UL(2)/ULSQ
      DIRP(1,3)=UL(3)/ULSQ
      UL(1)=( DIRP(1,2)*VL(3)-DIRP(1,3)*VL(2))
      UL(2)=(-DIRP(1,1)*VL(3)+DIRP(1,3)*VL(1))
      UL(3)=( DIRP(1,1)*VL(2)-DIRP(1,2)*VL(1))
      ULSQ=SQRT(UL(1)*UL(1)+UL(2)*UL(2)+UL(3)*UL(3))
      DIRP(3,1)=UL(1)/ULSQ
      DIRP(3,2)=UL(2)/ULSQ
      DIRP(3,3)=UL(3)/ULSQ
      DIRP(2,1)= DIRP(3,2)*DIRP(1,3)-DIRP(3,3)*DIRP(1,2)
      DIRP(2,2)=-DIRP(3,1)*DIRP(1,3)+DIRP(3,3)*DIRP(1,1)
      DIRP(2,3)= DIRP(3,1)*DIRP(1,2)-DIRP(3,2)*DIRP(1,1)
C-----------------------------------------------------------------------
C     Evaluation of vertex coordinates relative to local coordinate
C     system.
C-----------------------------------------------------------------------
      DO 30 I=1,4
        XVP(I)=DIRP(1,1)*(V(1,I)-XORP(1))+
     &          DIRP(1,2)*(V(2,I)-XORP(2)) + DIRP(1,3)*(V(3,I)-XORP(3))
        YVP(I)=DIRP(2,1)*(V(1,I)-XORP(1))+
     &          DIRP(2,2)*(V(2,I)-XORP(2)) + DIRP(2,3)*(V(3,I)-XORP(3))
30    CONTINUE
C-----------------------------------------------------------------------
C     Calculation of local coordinates of Gaussian integration
C     nodes and corresponding weights. Evaluation of quadrilateral
C     centroid. Transfer of local system origin at centroid
C-----------------------------------------------------------------------
      A1 = (XVP(1)+XVP(2))*HALF
      A2 = (XVP(1)+XVP(4))*HALF
      A3 = (XVP(1)+XVP(3))*HALF
      B1 = (YVP(1)+YVP(2))*HALF
      B2 = (YVP(1)+YVP(4))*HALF
      B3 = (YVP(1)+YVP(3))*HALF
      D0 = A1*B2-A2*B1
      D1 = A1*B3-A3*B1
      D2 = A3*B2-A2*B3
      D0I = THIRD / D0
      UL(1) = (D1*A1+D2*A2)*D0I
      UL(2) = (D1*B1+D2*B2)*D0I
      DO 40 J = 1,4
        XVP(J) = XVP(J) - UL(1)
        YVP(J) = YVP(J) - UL(2)
40    CONTINUE
C-----------------------------------------------------------------------
C    Evaluation of coordinates of centroid and Gaussian nodes relative
C    to reference system
C-----------------------------------------------------------------------
      DO 60 J=1,3
        DSUM = DIRP(1,J) * UL(1) + DIRP(2,J) * UL(2)
        XCTP(J) = DSUM + XORP(J)
60    CONTINUE
C-----------------------------------------------------------------------
C     Evaluation of moments of quadrilateral surface relative to
C     local system, array S(15).  First initialize array
C     Note that S(2)=S(6)=0 due to transfer above
C-----------------------------------------------------------------------
      DO 70 N=1,15
        S(N)=ZERO
70    CONTINUE
C-----------------------------------------------------------------------
C     Evaluation of higher powers of vertex coordinates.
C-----------------------------------------------------------------------
      DO 80 N=1,4
        XP(1,N)=XVP(N)
        YP(1,N)=YVP(N)
        DO 75 M=2,6
          XP(M,N)=XP(1,N)*XP(M-1,N)
          YP(M,N)=YP(1,N)*YP(M-1,N)
75      CONTINUE
80    CONTINUE
      DO 120 N=1,4
        NXT=N1(N)
        DX(N)=XVP(N1(N))-XVP(N)
        DY(N)=YVP(N1(N))-YVP(N)
C-----------------------------------------------------------------------
C   Evaluate length of each side and accumulate area of panel
C-----------------------------------------------------------------------
        S(1)=S(1)+HALF*DX(N)*(YVP(N)+YVP(N1(N)))
        DDX=DX(N)*DX(N)
        DDY=DY(N)*DY(N)
        DDS=DDX+DDY
        IF (DDS.LT.TOLT) THEN
          SIDE(N)=ZERO
          GO TO 120
        ENDIF
        SIDE(N)=SQRT(DDS)
        IF (DDY.GT.DDX) GOTO 100
        DYDX=DY(N)/DX(N)
        DO 90 I=2,4
          I1=I+1
          I2=I+2
          SI=(XP(I1,NXT)*YP(1,NXT)-XP(I1,N)*YP(1,N))/(I1)
     &          +DYDX*(XP(I2,N)-XP(I2,NXT))/(I1*I2)
          K=K1(I1)
          S(K)=S(K)+SI
          DO 85 J=1,I
            J1=J+1
            IJ1=I-J+1
            SI=(XP(IJ1,NXT)*YP(J1,NXT)-XP(IJ1,N)*YP(J1,N))/(IJ1*J1)-
     &                                          DYDX*J*SI/IJ1
            K=J+K1(IJ1)
            S(K)=S(K)+SI
85        CONTINUE
90      CONTINUE
        GO TO 120
100     DXDY=DX(N)/DY(N)
        DO 110 I=2,4
          I1=I+1
          I2=I+2
          SI=DXDY/(I1*I2)*(YP(I2,NXT)-YP(I2,N))
          S(I1)=S(I1)+SI
          DO 105 J=1,I
            IJ=I-J
            IJ1=IJ+1
            IJ2=IJ+2
            SI=DXDY*((XP(J,NXT)*YP(IJ2,NXT)-XP(J,N)*YP(IJ2,N))/
     &                                  (IJ1*IJ2)-J*SI/IJ1)
            K=IJ+K1(J+1)
            S(K)=S(K)+SI
105       CONTINUE
110     CONTINUE
120   CONTINUE
      DO 130 N=3,15
        S(N)=S(N)*CS(N)
130   CONTINUE
      RETURN
      END

max_stars_repo_name: OregonIons/bem

max_issues_repo_name: OregonIons/bem

max_forks_repo_name: OregonIons/bem


! { dg-do run }
program p
 implicit none

 real   , parameter :: arr(3) = [ real    :: 2, 2.5, (1.5, 2.5) ]
 real   , parameter :: ari(3) = [ integer :: 2, 2.5, (1.5, 2.5) ]
 real   , parameter :: arc(3) = [ complex :: 2, 2.5, (1.5, 2.5) ]
 integer, parameter :: air(3) = [ real    :: 2, 2.5, (1.5, 2.5) ]
 integer, parameter :: aii(3) = [ integer :: 2, 2.5, (1.5, 2.5) ]
 integer, parameter :: aic(3) = [ complex :: 2, 2.5, (1.5, 2.5) ]
 complex, parameter :: acr(3) = [ real    :: 2, 2.5, (1.5, 2.5) ]
 complex, parameter :: aci(3) = [ integer :: 2, 2.5, (1.5, 2.5) ]
 complex, parameter :: acc(3) = [ complex :: 2, 2.5, (1.5, 2.5) ]

 real   , parameter :: mrr(3) =  4.5       * [ real    :: 2, 2.5, (3.5, 4.0) ]
 real   , parameter :: mri(3) =  4.5       * [ integer :: 2, 2.5, (3.5, 4.0) ]
 real   , parameter :: mrc(3) =  4.5       * [ complex :: 2, 2.5, (3.5, 4.0) ]
 integer, parameter :: mir(3) =  4         * [ real    :: 2, 2.5, (3.5, 4.0) ]
 integer, parameter :: mii(3) =  4         * [ integer :: 2, 2.5, (3.5, 4.0) ]
 integer, parameter :: mic(3) =  4         * [ complex :: 2, 2.5, (3.5, 4.0) ]
 complex, parameter :: mcr(3) = (4.5, 5.5) * [ real    :: 2, 2.5, (3.5, 4.0) ]
 complex, parameter :: mci(3) = (4.5, 5.5) * [ integer :: 2, 2.5, (3.5, 4.0) ]
 complex, parameter :: mcc(3) = (4.5, 5.5) * [ complex :: 2, 2.5, (3.5, 4.0) ]

 if (any(arr /= [2.00, 2.50, 1.50])) STOP 1
 if (any(ari /= [2.00, 2.00, 1.00])) STOP 2
 if (any(arc /= [2.00, 2.50, 1.50])) STOP 3

 if (any(air /= [2, 2, 1])) STOP 4
 if (any(aii /= [2, 2, 1])) STOP 5
 if (any(aic /= [2, 2, 1])) STOP 6

 if (any(acr /= [(2.00, 0.00), (2.50, 0.00), (1.50, 0.00)])) STOP 7
 if (any(aci /= [(2.00, 0.00), (2.00, 0.00), (1.00, 0.00)])) STOP 8
 if (any(acc /= [(2.00, 0.00), (2.50, 0.00), (1.50, 2.50)])) STOP 9

 if (any(mrr /= [9.00, 11.25, 15.75])) STOP 10
 if (any(mri /= [9.00,  9.00, 13.50])) STOP 11
 if (any(mrc /= [9.00, 11.25, 15.75])) STOP 12

 if (any(mir /= [8, 10, 14])) STOP 13
 if (any(mii /= [8,  8, 12])) STOP 14
 if (any(mic /= [8, 10, 14])) STOP 15

 if (any(mcr /= [(9.00, 11.00), (11.25, 13.75), (15.75, 19.25)])) STOP 16
 if (any(mci /= [(9.00, 11.00), ( 9.00, 11.00), (13.50, 16.50)])) STOP 17
 if (any(mcc /= [(9.00, 11.00), (11.25, 13.75), (-6.25, 37.25)])) STOP 18

end program p