MODULE OMSAO_solcomp_module

  ! =================================================================
  !
  ! This module defines variables associated with parameterization of
  ! the composite solar spectrum.
  !
  ! =================================================================

  USE OMSAO_precision_module,  ONLY: i4, r8, C_LONG
  USE OMSAO_omidata_module,    ONLY: nxtrack_max
  USE OMSAO_errstat_module

  IMPLICIT NONE

  ! ----------------------------------------------------------------------
  ! Type declaration for variable with composite solar spectrum parameters
  ! ----------------------------------------------------------------------
  TYPE, PUBLIC :: SolarCompositeParameters
     CHARACTER (LEN=maxchlen)                            :: Title
     INTEGER   (KIND=I4)                                 :: nBins, nPol, nXtr
     REAL      (KIND=r8)                                 :: SolarNorm
     REAL      (KIND=r8), DIMENSION (2)                  :: FirstLastWav
     REAL      (KIND=r8), DIMENSION (:),     ALLOCATABLE :: SolComParWavs
     REAL      (KIND=r8), DIMENSION (:,:,:), ALLOCATABLE :: SolComParData
  END TYPE SolarCompositeParameters

  ! -----------------------------------------------------------
  ! Variable that holds the tabulated solar spectrum parameters
  ! -----------------------------------------------------------
  TYPE (SolarCompositeParameters) :: solarcomp_pars

  ! ---------------------------------------------------
  ! The division orbit between OPF versions < and >= 25
  ! ---------------------------------------------------
  INTEGER (KIND=i4), PARAMETER :: opf_ge25 = 25

  ! ----------------------------------------------------------------
  ! Strings to aid in the construction of the HE5 field name to read
  ! ----------------------------------------------------------------
  CHARACTER (LEN= 4), DIMENSION (3), PARAMETER :: ch_typ = (/ 'MEA_', 'MED_', 'PC0_' /)
  CHARACTER (LEN= 8), DIMENSION (2), PARAMETER :: ch_opf = (/ 'OPFlt25_', 'OPFge25_' /)
  CHARACTER (LEN= 5), DIMENSION (2), PARAMETER :: ch_orb = (/ '1st7_',    'comp_'    /)
  CHARACTER (LEN= 3), DIMENSION (3), PARAMETER :: ch_cha = (/ 'UV1',  'UV2',   'VIS' /)

  ! ---------------------------
  ! 32bit/64bit C_LONG integers
  ! ---------------------------
  INTEGER (KIND=C_LONG), PARAMETER, PRIVATE :: zerocl = 0, onecl = 1

CONTAINS

  SUBROUTINE soco_pars_read ( &
       soco_fname, typ_idx, l1b_channel, winwav_min, winwav_max, errstat )

    USE OMSAO_parameters_module, ONLY: maxchlen
    IMPLICIT NONE

    INCLUDE 'hdfeos5.inc'

    ! ---------------------------------------------------------------------
    ! Explanation of subroutine arguments:
    !
    !    soco_fname ........... file name with Solar Composite spectra
    !    typ_idx .............. whether "MEA", "MED" or "PC0" is requested
    !    l1b_channel .......... OMI channel ("UV1", "UV2", or "VIS")
    !    winwav_min ........... starting wavelength to read from file
    !    winwav_max ........... final wavelength to read from file
    !    errstat .............. error status returned from the subroutine
    ! ---------------------------------------------------------------------

    CHARACTER (LEN=14), PARAMETER :: modulename = 'soco_pars_read'

    ! ---------------
    ! Input variables
    ! ---------------
    INTEGER (KIND=i4), INTENT (IN) :: typ_idx
    REAL    (KIND=r8), INTENT (IN) :: winwav_min, winwav_max
    CHARACTER (LEN=*), INTENT (IN) :: soco_fname
    CHARACTER (LEN=3), INTENT (IN) :: l1b_channel

    ! -----------------
    ! Modified variable
    ! -----------------
    INTEGER (KIND=i4), INTENT (INOUT) :: errstat

    ! --------------
    ! Local variable
    ! --------------
    INTEGER (KIND=i4) :: n_xtr, n_wvl, n_pars
    INTEGER (KIND=i4) :: ios, ntmp, j1, j2, swlen
    REAL    (KIND=r8) :: soco_norm
    REAL    (KIND=r8), DIMENSION (:),     ALLOCATABLE :: tmp_wvls
    REAL    (KIND=r8), DIMENSION (:,:,:), ALLOCATABLE :: tmp_pars

    INTEGER   (KIND=i4)                :: he5stat, swath_file_id, swath_id
    INTEGER   (KIND=C_LONG)               :: he5statcl, he5_start_1, he5_stride_1, he5_edge_1
    INTEGER   (KIND=C_LONG), DIMENSION(3) :: he5_start_3, he5_stride_3, he5_edge_3
    CHARACTER (LEN=maxchlen)           :: soco_he5_field, swath_name

    ! ---------------------------------
    ! External OMI and Toolkit routines
    ! ---------------------------------
    INTEGER (KIND=i4),  EXTERNAL :: &
         HE5_SWopen, HE5_SWattach, HE5_SWrdfld, HE5_SWrdlattr, HE5_SWdetach, HE5_SWclose
    INTEGER (KIND=C_LONG), EXTERNAL :: &
         HE5_SWinqswath

    ! -----------------------
    ! Initialize error status
    ! -----------------------
    errstat = pge_errstat_ok

    ! -----------------------------------------------------------
    ! Open HE5 output file and check SWATH_FILE_ID ( -1 if error)
    ! -----------------------------------------------------------
    swath_file_id = HE5_SWopen ( TRIM(ADJUSTL(soco_fname)), he5f_acc_rdonly )
    IF ( swath_file_id == he5_stat_fail ) THEN
       CALL error_check ( &
            0, 1, pge_errstat_error, OMSAO_E_HE5SWOPEN, modulename, vb_lev_default, he5stat )
       errstat = MAX ( errstat, he5stat )
       RETURN
    END IF
    ! ---------------------------------------------
    ! Check for existing HE5 swath and attach to it
    ! ---------------------------------------------
    swath_name = "" !JED fix
    he5statcl  = HE5_SWinqswath  ( TRIM(ADJUSTL(soco_fname)), swath_name, swlen )
    swath_id = HE5_SWattach ( swath_file_id, TRIM(ADJUSTL(swath_name)) )
    IF ( swath_id == he5_stat_fail ) THEN
       CALL error_check ( &
            0, 1, pge_errstat_error, OMSAO_E_HE5SWATTACH, modulename, vb_lev_default, he5stat )
       errstat = MAX ( errstat, he5stat )
       RETURN
    END IF

    ! -----------------------------
    ! Read dimensions from HE5 file
    ! -----------------------------
    CALL soco_get_dims ( swath_id, l1b_channel, n_xtr, n_wvl, n_pars )

    ! ------------------------------------
    ! Allocate memory for temporary arrays
    ! ------------------------------------
    ALLOCATE ( tmp_wvls(1:n_wvl),                  STAT=ios )
    ALLOCATE ( tmp_pars(1:n_xtr,1:n_wvl,1:n_pars), STAT=ios )

    ! --------------------
    ! Read wavelength bins
    ! --------------------
    he5_start_1 = zerocl ; he5_stride_1 = onecl ; he5_edge_1 = INT(n_wvl,KIND=C_LONG)
    he5stat = HE5_SWrdfld ( swath_id, "SoCo_Wvl_"//l1b_channel,   &
         he5_start_1, he5_stride_1, he5_edge_1, tmp_wvls(1:n_wvl) )
    IF ( he5stat /= pge_errstat_ok ) &
         CALL error_check ( he5stat, OMI_S_SUCCESS, pge_errstat_error, OMSAO_E_PREFITCOL, &
         modulename//f_sep//'Solar Composite access failed.', vb_lev_default, errstat )

    ! ----------------------
    ! Compose HE5 field name
    ! ----------------------
    soco_he5_field = ' '
    soco_he5_field = ch_typ(typ_idx)//l1b_channel

    ! -------------------------------
    ! Read Solar Composite parameters
    ! -------------------------------
    he5_start_3  = (/               zerocl,               zerocl,                zerocl /)
    he5_stride_3 = (/               onecl,               onecl,                onecl /)
    he5_edge_3   = (/ INT(n_xtr,KIND=C_LONG), INT(n_wvl,KIND=C_LONG), INT(n_pars,KIND=C_LONG) /)
    he5stat = HE5_SWrdfld ( swath_id, TRIM(ADJUSTL(soco_he5_field)),   &
         he5_start_3, he5_stride_3, he5_edge_3, tmp_pars(1:n_xtr,1:n_wvl,1:n_pars) )
    IF ( he5stat /= pge_errstat_ok ) &
         CALL error_check ( he5stat, OMI_S_SUCCESS, pge_errstat_error, OMSAO_E_PREFITCOL, &
         modulename//f_sep//'Solar Composite access failed.', vb_lev_default, errstat )

    ! -------------------------
    ! Read Solar Norm attribute
    ! -------------------------
    he5stat = HE5_SWrdlattr ( &
         swath_id, TRIM(ADJUSTL(soco_he5_field)), "SolarSpectrumNorm", soco_norm )

    ! -----------------------------------------------
    ! Detach from HE5 swath and close HE5 output file
    ! -----------------------------------------------
    he5stat = HE5_SWdetach ( swath_id )
    he5stat = HE5_SWclose  ( swath_file_id )

    ! --------------------------------------------
    ! Find the slice to store in the global arrays
    ! --------------------------------------------
    CALL array_locate_r8 ( n_wvl, tmp_wvls(1:n_wvl), winwav_min, 'LE', j1 )
    CALL array_locate_r8 ( n_wvl, tmp_wvls(1:n_wvl), winwav_max, 'GE', j2 )

    j1 = MAX ( j1, 1 )
    IF ( j2 < 1 ) j2 = n_wvl

    ntmp = j2 - j1 + 1

    ! ----------------
    ! Assign the title
    ! ----------------
    solarcomp_pars%Title = 'OMI Solar Composite Spectrum --- '//TRIM(ADJUSTL(soco_he5_field))

    ! --------------------------------
    ! Assign variables to global array
    ! --------------------------------
    solarcomp_pars%nXtr  = ntmp
    solarcomp_pars%nBins = ntmp
    solarcomp_pars%nPol  = n_pars

    ALLOCATE ( solarcomp_pars%SolComParWavs(1:ntmp),                  STAT=ios )
    ALLOCATE ( solarcomp_pars%SolComParData(1:n_xtr,1:ntmp,1:n_pars), STAT=ios )

    solarcomp_pars%SolComParWavs(1:ntmp)                  = tmp_wvls(j1:j2)
    solarcomp_pars%SolComParData(1:n_xtr,1:ntmp,1:n_pars) = tmp_pars(1:n_xtr,j1:j2,1:n_pars)

    solarcomp_pars%SolarNorm = soco_norm

    ! ---------------------------------------
    ! De-allocate memory for temporary arrays
    ! ---------------------------------------
    IF ( ALLOCATED(tmp_wvls) ) DEALLOCATE (tmp_wvls)
    IF ( ALLOCATED(tmp_pars) ) DEALLOCATE (tmp_pars)

    RETURN
  END SUBROUTINE soco_pars_read

  SUBROUTINE soco_pars_read_v002 ( &
       soco_fname, typ_idx, sel_idx, l1b_channel, l1br_opf_version, winwav_min, winwav_max, errstat )

    ! -----------------------------------------------------------
    ! Old version (Collection 2 L1b data) of SOCO_PARS_READ above
    ! -----------------------------------------------------------

    USE OMSAO_parameters_module, ONLY: maxchlen
    IMPLICIT NONE

    INCLUDE 'hdfeos5.inc'

    ! ---------------------------------------------------------------------
    ! Explanation of subroutine arguments:
    !
    !    soco_fname ........... file name with Solar Composite spectra
    !    typ_idx .............. whether "MED" or "PC0" is requested
    !    sel_idx .............. whether "1st7" or "comp" is requested
    !    l1b_channel .......... OMI channel ("UV1", "UV2", or "VIS")
    !    l1br_opf_version ..... OPF version of L1B radiance file
    !    winwav_min ........... starting wavelength to read from file
    !    winwav_max ........... final wavelength to read from file
    !    errstat .............. error status returned from the subroutine
    ! ---------------------------------------------------------------------

    CHARACTER (LEN=19), PARAMETER :: modulename = 'soco_pars_read_v003'

    ! ---------------
    ! Input variables
    ! ---------------
    INTEGER (KIND=i4), INTENT (IN) :: typ_idx, sel_idx, l1br_opf_version
    REAL    (KIND=r8), INTENT (IN) :: winwav_min, winwav_max
    CHARACTER (LEN=*), INTENT (IN) :: soco_fname
    CHARACTER (LEN=3), INTENT (IN) :: l1b_channel

    ! -----------------
    ! Modified variable
    ! -----------------
    INTEGER (KIND=i4), INTENT (INOUT) :: errstat

    ! --------------
    ! Local variable
    ! --------------
    INTEGER (KIND=i4) :: n_xtr, n_wvl, n_pars
    INTEGER (KIND=i4) :: ios, ntmp, j1, j2, opf_idx, swlen
    REAL    (KIND=r8) :: soco_norm
    REAL    (KIND=r8), DIMENSION (:),     ALLOCATABLE :: tmp_wvls
    REAL    (KIND=r8), DIMENSION (:,:,:), ALLOCATABLE :: tmp_pars

    INTEGER   (KIND=i4)                :: he5stat, swath_file_id, swath_id
    INTEGER   (KIND=C_LONG)               :: he5statcl, he5_start_1, he5_stride_1, he5_edge_1
    INTEGER   (KIND=C_LONG), DIMENSION(3) :: he5_start_3, he5_stride_3, he5_edge_3
    CHARACTER (LEN=maxchlen)           :: soco_he5_field, swath_name

    ! ---------------------------------
    ! External OMI and Toolkit routines
    ! ---------------------------------
    INTEGER (KIND=i4),  EXTERNAL :: &
         HE5_SWopen, HE5_SWattach, HE5_SWrdfld, HE5_SWrdlattr, HE5_SWdetach, HE5_SWclose
    INTEGER (KIND=C_LONG), EXTERNAL :: &
         HE5_SWinqswath

    ! -----------------------
    ! Initialize error status
    ! -----------------------
    errstat = pge_errstat_ok

    ! --------------------------------------------------------------------------
    ! Select either pre- or post-OPFv25 data. The OPF version has been read from
    ! the L1b radiance file in INIT_METADATA.
    ! --------------------------------------------------------------------------
    IF ( l1br_opf_version >= opf_ge25 ) THEN
       opf_idx = 2
    ELSE
       opf_idx = 1
    END IF

    ! -----------------------------------------------------------
    ! Open HE5 output file and check SWATH_FILE_ID ( -1 if error)
    ! -----------------------------------------------------------
    swath_file_id = HE5_SWopen ( TRIM(ADJUSTL(soco_fname)), he5f_acc_rdonly )
    IF ( swath_file_id == he5_stat_fail ) THEN
       CALL error_check ( &
            0, 1, pge_errstat_error, OMSAO_E_HE5SWOPEN, modulename, vb_lev_default, he5stat )
       errstat = MAX ( errstat, he5stat )
       RETURN
    END IF

    ! ---------------------------------------------
    ! Check for existing HE5 swath and attach to it
    ! ---------------------------------------------
    he5stat  = HE5_SWinqswath  ( TRIM(ADJUSTL(soco_fname)), swath_name, swlen )
    swath_id = HE5_SWattach ( swath_file_id, TRIM(ADJUSTL(swath_name)) )
    IF ( swath_id == he5_stat_fail ) THEN
       CALL error_check ( &
            0, 1, pge_errstat_error, OMSAO_E_HE5SWATTACH, modulename, vb_lev_default, he5stat )
       errstat = MAX ( errstat, he5stat )
       RETURN
    END IF

    ! -----------------------------
    ! Read dimensions from HE5 file
    ! -----------------------------
    CALL soco_get_dims ( swath_id, l1b_channel, n_xtr, n_wvl, n_pars )
    ! ------------------------------------
    ! Allocate memory for temporary arrays
    ! ------------------------------------
    ALLOCATE ( tmp_wvls(1:n_wvl),                  STAT=ios )
    ALLOCATE ( tmp_pars(1:n_xtr,1:n_wvl,1:n_pars), STAT=ios )

    ! --------------------
    ! Read wavelength bins
    ! --------------------
    he5_start_1 = 0 ; he5_stride_1 = 1 ; he5_edge_1 = n_wvl
    he5stat = HE5_SWrdfld ( swath_id, "SoCo_Wvl_"//l1b_channel,   &
         he5_start_1, he5_stride_1, he5_edge_1, tmp_wvls(1:n_wvl) )
    IF ( he5stat /= pge_errstat_ok ) &
         CALL error_check ( he5stat, OMI_S_SUCCESS, pge_errstat_error, OMSAO_E_PREFITCOL, &
         modulename//f_sep//'Solar Composite access failed.', vb_lev_default, errstat )

    ! ----------------------
    ! Compose HE5 field name
    ! ----------------------
    soco_he5_field = ' '
    soco_he5_field = ch_typ(typ_idx)//ch_opf(opf_idx)//ch_orb(sel_idx)//l1b_channel

    ! -------------------------------
    ! Read Solar Composite parameters
    ! -------------------------------
    he5_start_3  = (/     0,     0,      0 /)
    he5_stride_3 = (/     1,     1,      1 /)
    he5_edge_3   = (/ n_xtr, n_wvl, n_pars /)
    he5stat = HE5_SWrdfld ( swath_id, TRIM(ADJUSTL(soco_he5_field)),   &
         he5_start_3, he5_stride_3, he5_edge_3, tmp_pars(1:n_xtr,1:n_wvl,1:n_pars) )
    IF ( he5stat /= pge_errstat_ok ) &
         CALL error_check ( he5stat, OMI_S_SUCCESS, pge_errstat_error, OMSAO_E_PREFITCOL, &
         modulename//f_sep//'Solar Composite access failed.', vb_lev_default, errstat )

    ! -------------------------
    ! Read Solar Norm attribute
    ! -------------------------
    he5stat = HE5_SWrdlattr ( &
         swath_id, TRIM(ADJUSTL(soco_he5_field)), "SolarSpectrumNorm", soco_norm )

    ! -----------------------------------------------
    ! Detach from HE5 swath and close HE5 output file
    ! -----------------------------------------------
    he5stat = HE5_SWdetach ( swath_id )
    he5stat = HE5_SWclose  ( swath_file_id )

    ! --------------------------------------------
    ! Find the slice to store in the global arrays
    ! --------------------------------------------
    CALL array_locate_r8 ( n_wvl, tmp_wvls(1:n_wvl), winwav_min, 'LE', j1 )
    CALL array_locate_r8 ( n_wvl, tmp_wvls(1:n_wvl), winwav_max, 'GE', j2 )

    j1 = MAX ( j1, 1 )
    IF ( j2 < 1 ) j2 = n_wvl

    ntmp = j2 - j1 + 1

    ! ----------------
    ! Assign the title
    ! ----------------
    solarcomp_pars%Title = 'OMI Solar Composite Spectrum --- '//TRIM(ADJUSTL(soco_he5_field))

    ! --------------------------------
    ! Assign variables to global array
    ! --------------------------------
    solarcomp_pars%nXtr  = ntmp
    solarcomp_pars%nBins = ntmp
    solarcomp_pars%nPol  = n_pars

    ALLOCATE ( solarcomp_pars%SolComParWavs(1:ntmp),                  STAT=ios )
    ALLOCATE ( solarcomp_pars%SolComParData(1:n_xtr,1:ntmp,1:n_pars), STAT=ios )

    solarcomp_pars%SolComParWavs(1:ntmp)                  = tmp_wvls(j1:j2)
    solarcomp_pars%SolComParData(1:n_xtr,1:ntmp,1:n_pars) = tmp_pars(1:n_xtr,j1:j2,1:n_pars)

    solarcomp_pars%SolarNorm = soco_norm

    ! ---------------------------------------
    ! De-allocate memory for temporary arrays
    ! ---------------------------------------
    IF ( ALLOCATED(tmp_wvls) ) DEALLOCATE (tmp_wvls)
    IF ( ALLOCATED(tmp_pars) ) DEALLOCATE (tmp_pars)

    RETURN
  END SUBROUTINE soco_pars_read_v002

  SUBROUTINE soco_get_dims ( swath_id, l1b_channel, n_xtr, n_wvl, n_pars )

    IMPLICIT NONE

    ! ---------------
    ! Input variables
    ! ---------------
    INTEGER (KIND=i4), INTENT (IN) :: swath_id
    CHARACTER (LEN=3), INTENT (IN) :: l1b_channel

    ! ----------------
    ! Output variables
    ! ----------------
    INTEGER (KIND=i4), INTENT (OUT) :: n_xtr, n_wvl, n_pars

    ! ---------------
    ! Local variables
    ! ---------------
    INTEGER (KIND=C_LONG)                  :: ndimcl
    INTEGER (KIND=i4)                   :: ndim, swlen, istart, iend, nsep, i
    INTEGER (KIND=i4),  DIMENSION (0:9) :: dim_array, dim_seps
    INTEGER (KIND=C_LONG), DIMENSION (0:9) :: dim_arraycl
    CHARACTER (LEN=maxchlen)            :: dim_chars

    ! ------------------
    ! EXTERNAL functions
    ! ------------------
    INTEGER (KIND=C_LONG), EXTERNAL :: HE5_SWinqdims

    ! ------------------------------
    ! Inquire about swath dimensions
    ! ------------------------------
    ndimcl = HE5_SWinqdims  ( swath_id, dim_chars, dim_arraycl(0:6) )
    ndim   = INT ( ndimcl, KIND=i4 )
    dim_array(0:6) = INT ( dim_arraycl(0:6), KIND=i4 )

    ! ---------------------------
    ! Initialize output variables
    ! ---------------------------
    n_xtr = 0 ; n_wvl = 0 ; n_pars = 0

    ! ----------------------------------------------------------------------
    ! Find the positions of separators (commas, ",") between the dimensions.
    ! Add a "pseudo separator" at the end to fully automate the consecutive
    ! check for nTimes and nXtrack.
    ! ----------------------------------------------------------------------
    dim_chars = TRIM(ADJUSTL(dim_chars))
    swlen = LEN_TRIM(ADJUSTL(dim_chars))
    istart = 1  ;  iend = 1 ; nsep = 0 ; dim_seps = 0

    getseps: DO i = 1, swlen
       IF ( dim_chars(i:i) == ',' ) THEN
          nsep = nsep + 1
          dim_seps(nsep) = i
       END IF
    END DO getseps
    nsep = nsep + 1 ; dim_seps(nsep) = swlen+1

    ! --------------------------------------------------------------------
    ! Hangle along the NSEP indices until we have found the two dimensions
    ! we are interested in.
    ! --------------------------------------------------------------------
    getdims:DO i = 0, nsep-1
       istart = dim_seps(i)+1 ; iend = dim_seps(i+1)-1
       IF  ( dim_chars(istart:iend) == "N_Xtr_"//l1b_channel  ) n_xtr  = dim_array(i)
       IF  ( dim_chars(istart:iend) == "N_Wvl_"//l1b_channel  ) n_wvl  = dim_array(i)
       IF  ( dim_chars(istart:iend) == "N_SoCo_Pars"          ) n_pars = dim_array(i)
       IF ( n_xtr > 0 .AND. n_wvl > 0 .AND. n_pars > 0 ) EXIT getdims
    END DO getdims

    RETURN
  END SUBROUTINE soco_get_dims

  SUBROUTINE soco_compute ( yn_norm, ixt, nwvl, wvl, spc )

    IMPLICIT NONE

    ! ------------------------------------------------------------------------------
    ! Explanation of subroutine arguments:
    !
    !    ixt .................. current cross-track position
    !    nwvl ................. number of wavlengths to be computed
    !    wvl .................. wavlength at which to compute spectrum
    !    spc .................. the solar composite spectrum
    ! ------------------------------------------------------------------------------

    ! ---------------
    ! Input variables
    ! ---------------
    LOGICAL,                             INTENT (IN) :: yn_norm
    INTEGER (KIND=i4),                   INTENT (IN) :: ixt, nwvl
    REAL    (KIND=r8), DIMENSION (nwvl), INTENT (IN) :: wvl

    ! ----------------
    ! Output variables
    ! ----------------
    REAL    (KIND=r8), DIMENSION (nwvl), INTENT (OUT) :: spc

    ! ---------------
    ! Local variables
    ! ---------------
    INTEGER (KIND=i4) :: iwvl, iloc, ipow, nbin, npol
    REAL    (KIND=r8) :: binwvl, binpow, snorm

    ! -------------------------------
    ! Name of the function/subroutine
    ! -------------------------------
    CHARACTER (LEN=12), PARAMETER :: modulename = 'soco_compute'

    nbin  = solarcomp_pars%nBins
    npol  = solarcomp_pars%nPol
    snorm = solarcomp_pars%SolarNorm

    ! ---------------------------------------------------------------------------
    ! Adding the contributions one-by-one and power-by-power is a bit pedestrian,
    ! but anything else involves allocating and deallocating memory, and we don't
    ! really want to go to that trouble.
    ! ---------------------------------------------------------------------------
    DO iwvl = 1, nwvl
       ! ----------------------------------------------------------------------------
       ! CAREFUL: Prevent OOB if "iloc = n_soco_bin" or "wvl = soco_wvl(n_soco_bins)"
       ! ----------------------------------------------------------------------------
       IF ( wvl(iwvl) <= solarcomp_pars%SolComParWavs(1) ) THEN
          iloc = 1 ; binwvl = -1.0_r8
       ELSE IF ( wvl(iwvl) >= solarcomp_pars%SolComParWavs(nbin) ) THEN
          iloc = nbin ; binwvl = 1.0_r8
       ELSE
          CALL array_locate_r8 ( nbin, solarcomp_pars%SolComParWavs(1:nbin), wvl(iwvl), 'LE', iloc )
          binwvl = 2.0_r8*(wvl(iwvl)-solarcomp_pars%SolComParWavs(iloc)) / &
               (solarcomp_pars%SolComParWavs(iloc+1)-solarcomp_pars%SolComParWavs(iloc)) - 1.0_r8
       END IF
       binpow     = 1.0_r8
       spc (iwvl) = solarcomp_pars%SolComParData(ixt,iloc,1)
       DO ipow = 2, npol
          binpow    = binpow * binwvl
          spc(iwvl) = spc(iwvl) + solarcomp_pars%SolComParData(ixt,iloc,ipow) * binpow
       END DO

    END DO

    ! -----------------------------------------------------------------
    ! Multiply by norm of spectrum, unless spectrum is to be normalized
    ! (YN_NORM = .TRUE.)
    ! -----------------------------------------------------------------
    !!! tpk(11/09/2010) IF ( .NOT. yn_norm ) spc(1:nwvl) = spc(1:nwvl) * snorm
    IF ( .NOT. yn_norm ) spc(1:nwvl) = spc(1:nwvl) * snorm

    ! -------------------------------------------------------------------
    ! Multiply by norm of spectrum anyway: It is better to take care of
    ! normalization in the routine that computes the final spectrum norm.
    ! -------------------------------------------------------------------
    !!!spc(1:nwvl) = spc(1:nwvl) * snorm

    RETURN
  END SUBROUTINE soco_compute

  SUBROUTINE soco_pars_deallocate ( errstat )

    IMPLICIT NONE

    ! ---------------------------------------------------------------------
    ! Explanation of subroutine arguments:
    !
    !    errstat .............. error status returned from the subroutine
    ! ---------------------------------------------------------------------

    ! -----------------
    ! Modified variable
    ! -----------------
    INTEGER (KIND=i4), INTENT (INOUT) :: errstat

    errstat = pge_errstat_ok

    ! ------------------------------------
    ! De-allocate memory for global arrays
    ! ------------------------------------
    IF ( ALLOCATED(solarcomp_pars%SolComParWavs) ) DEALLOCATE (solarcomp_pars%SolComParWavs)
    IF ( ALLOCATED(solarcomp_pars%SolComParData) ) DEALLOCATE (solarcomp_pars%SolComParData)

    RETURN
  END SUBROUTINE soco_pars_deallocate

END MODULE OMSAO_solcomp_module


text:
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!
!! Copyright (c) 2018, Neil N. Carlson
!!
!! Permission is hereby granted, free of charge, to any person obtaining a
!! copy of this software and associated documentation files (the "Software"),
!! to deal in the Software without restriction, including without limitation
!! the rights to use, copy, modify, merge, publish, distribute, sublicense,
!! and/or sell copies of the Software, and to permit persons to whom the
!! Software is furnished to do so, subject to the following conditions:
!!
!! The above copyright notice and this permission notice shall be included
!! in all copies or substantial portions of the Software.
!!
!! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
!! IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
!! FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
!! THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
!! LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
!! FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
!! DEALINGS IN THE SOFTWARE.
!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

use json
implicit none

class(json_value), allocatable :: val
type(json_array_iterator) :: iter1, iter3
type(json_object_iterator)  :: iter2
character(:), allocatable :: errmsg
integer :: stat, n1, n3

call json_from_string('[{"a":[42]}]', val, stat, errmsg)
if (stat /= 0) stop 1

select type (val)
type is (json_array)

  n1 = 0
  iter1 = json_array_iterator(val)
  do while (.not.iter1%at_end())
    n1 = n1 + 1
    select type (ival1 => iter1%value())
    type is (json_object)
      if (n1 /= 1) stop 11
      iter2 = json_object_iterator(ival1)
      do while (.not.iter2%at_end())
        select case (iter2%name())
        case ('a')
          select type (ival2 => iter2%value())
          type is (json_array)
            n3 = 0
            iter3 = json_array_iterator(ival2)
            do while (.not.iter3%at_end())
              n3 = n3 + 1
              select type (ival3 => iter3%value())
              type is (json_integer)
                if (n3 /= 1) stop 10
                if (ival3%value /= 42) stop 9
              class default
                stop 8
              end select
              call iter3%next
            end do
            if (n3 /= 1) stop 7
          class default
            stop 6
          end select
        case default
          stop 5
        end select
        call iter2%next
      end do
    class default
      stop 4
    end select
    call iter1%next
  end do
  if (n1 /= 1) stop 3

class default
  stop 2
end select

end

max_stars_repo_name: nncarlson/yajl-fort

max_issues_repo_name: nncarlson/yajl-fort

max_forks_repo_name: nncarlson/yajl-fort


text:
      program detilt_hgt

c****************************************************************
c**
c**   FILE NAME: detilt_hgt.f
c**
c**   DATE WRITTEN:
c**
c**   PROGRAMMER: Scott Hensley
c**
c**   FUNCTIONAL DESCRIPTION: Remove cross track, along track
c**   tilts in height data and adds an offset.
c**
c**   ROUTINES CALLED:
c**
c**   NOTES:
c**
c**   UPDATE LOG:
c**
c**   Date Changed        Reason Changed                  CR # and Version #
c**   ------------       ----------------                 -----------------
c**
c*****************************************************************

      implicit none

c     INCLUDE FILES:

c     PARAMETER STATEMENTS:

      integer SIZE
      parameter(SIZE=30000)

c     INPUT VARIABLES:

c     OUTPUT VARIABLES:

c     LOCAL VARIABLES:

      character*3 a_type
      real a(SIZE),b(SIZE),sum,pha,phd,c(SIZE),d(SIZE),e(SIZE),suma
      character*120 fin,fout,a_string
      real*4 pa,pd,r_h0,r_null

      integer na,la,ld,lineout,line,nd,k,jpix,i,j,i_recmult

c     COMMON BLOCKS:

c     EQUIVALENCE STATEMENTS:

c     DATA STATEMENTS:

c     FUNCTION STATEMENTS:

      integer iargc

c     SAVE STATEMENTS:

c     PROCESSING STEPS:

      write(6,*) ' '
      write(6,*) '   <<   DELTILT HEIGHT DATA   >> '
      write(6,*) ' '

      if(iargc() .lt. 10)then
         write(6,'(a)') 'Usage: detilt_hgt infile outfile type samples looks_xt looks_at delhgt_xt delhgt_at h0 null'
         write(6,'(a)') 'Units: delhgt_at and delhgt_xt are meters/pixel '
         write(6,'(a)') 'type: s = separate file , rmg = RMG file'
         write(6,*) ' '
         stop
      endif

      call getarg(1,fin)
      call getarg(2,fout)
      call getarg(3,a_type)
      call getarg(4,a_string)
      read(a_string,*) na
      call getarg(5,a_string)
      read(a_string,*) la
      call getarg(6,a_string)
      read(a_string,*) ld
      call getarg(7,a_string)
      read(a_string,*) pa
      call getarg(8,a_string)
      read(a_string,*) pd
      call getarg(9,a_string)
      read(a_string,*) r_h0
      call getarg(10,a_string)
      read(a_string,*) r_null

c     compute height ramps ramps

      pha = pa
      phd = pd

c     opne input and output files

      if(index(a_type,'rmg') .ne. 0)then
         open(21,file=fin,form='unformatted',access='direct',recl=na*8)
         open(22,file=fout,form='unformatted',access='direct',recl=na/la*8)
         i_recmult = 2
      elseif(index(a_type,'s') .ne. 0)then
         open(21,file=fin,form='unformatted',access='direct',recl=na*4)
         open(22,file=fout,form='unformatted',access='direct',recl=na/la*4)
         i_recmult = 1
      endif

      lineout=0
      nd = 10000000

      write(6,*) ' '
      write(6,'(a,x,i5)') 'Number of output samples: ',na/la
      write(6,*) ' '

      do line=1,nd,ld

         if(mod(line,16*ld) .eq. 1)then
            write(6,'(a,x,i10)') 'At line: ',line
         endif

         lineout = lineout + 1

         do j=1,2*na
            b(j) = 0.0
            c(j) = 0.0
         end do

c     take looks down

         do i=0,ld-1

            read(21,rec=line+i,err=999)(a(k),k=1,i_recmult*na)
            if(index(a_type,'rmg') .ne. 0)then
               do j=1,na
                  if(a(j+na) .ne. r_null)then
                     c(j) = c(j) + 1
                     b(j) = b(j) + a(j)
                     b(j+na) = b(j+na) + a(j+na) + pha*j + phd*(line+i)
                  endif
               end do
            elseif(index(a_type,'s') .ne. 0)then
               do j=1,na
                  if(a(j) .ne. r_null)then
                     c(j) = c(j) + 1
                     b(j) = b(j) + a(j) + pha*j + phd*(line+i)
                  endif
               end do
            endif

         end do

c  take looks across

         jpix = 0
         do j=1,na,la
            jpix = jpix + 1
            e(jpix) = 0.0
            sum = 0.0
            suma = 0.0
            do k=0,la-1
               if(index(a_type,'rmg') .ne. 0)then
                  e(jpix) = e(jpix) + c(j+k)
                  sum = sum + b(j+k+na)
                  suma = suma + b(j+k)
               elseif(index(a_type,'s') .ne. 0)then
                  e(jpix) = e(jpix) + c(j+k)
                  sum = sum + b(j+k)
               endif
            end do
            if(index(a_type,'rmg') .ne. 0)then
               if(e(jpix) .ne. 0)then
                  d(jpix+na/la) = sum/e(jpix) + r_h0
                  d(jpix) = suma/e(jpix)
               else
                  d(jpix+na/la) = r_null
                  d(jpix) = 0.0
               endif
            elseif(index(a_type,'s') .ne. 0)then
               if(e(jpix) .ne. 0)then
                  d(jpix) = sum/e(jpix) + r_h0
               else
                  d(jpix) = r_null
               endif
            endif
         end do

         write(22,rec=lineout) (d(k),k=1,i_recmult*(na/la))

      end do

 999  write(6,*) ' '
      write(6,'(a)') 'Less resolution and with luck less noise...'

      end

ext: f

max_stars_repo_name: glshort/MapReady

max_issues_repo_name: glshort/MapReady

max_forks_repo_name: glshort/MapReady


text:
*
*     -------------------------------------------------------------
*      S S C
*     -------------------------------------------------------------
*                                                                  *
*                                                                  *
*     THIS PACKAGE DETERMINES THE VALUES OF MATRIX ELEMENTS        *
*     OF ONE PARTICLE OPERATOR IN CASE :           N'1 = N1        *
*                                                                  *
*                                                                  *
*     Written by G. Gaigalas,                                      *
*     Vanderbilt University Nashville, USA           October 1995  *
*
      SUBROUTINE SSC(IG,KL,IA,IB,IC,ID,IIA,IIB,IIC,IID,IREZ,XXX)
      IMPLICIT DOUBLEPRECISION (A-H,O-Z)
      COMMON/CONSTS/ZERO,TENTH,HALF,ONE,TWO,THREE,FOUR,SEVEN,ELEVEN,EPS
      COMMON/MEDEFN/IHSH,NJ(16),LJ(16),NOSH1(16),NOSH2(16),J1QN1(31,3),
     :     J1QN2(31,3),IJFUL(16)
      COMMON/DIAGNL/IDIAG,JA,JB
      EXTERNAL XXX
      KLP=KL+2
C      IF(IG-1.LT.KLP) RETURN
      L1=LJ(IIA)
      L2=LJ(IIB)
      L3=LJ(IIC)
      L4=LJ(IID)
      CALL SSA(L1,L2,L3,L4,KLP,KL,A1)
      CALL SSA(L2,L1,L4,L3,KLP,KL,A2)
      IF((DABS(A1)+DABS(A2)).LT.TWO*EPS) RETURN
      LA=IJFUL(IIA)
      LB=IJFUL(IIB)
      LC=IJFUL(IIC)
      LD=IJFUL(IID)
      CALL XXX(KLP,1,KL,1,2,IA,IB,IC,ID,IREZ,C,CC)
      C=C*A1
      CC=CC*A2
      IF(DABS(C).GT.EPS)CALL SAVENON(8,C,KL,LA,LB,LC,LD,JA,JB,0)
      IF(DABS(CC).GT.EPS)CALL SAVENON(8,CC,KL,LB,LA,LD,LC,JA,JB,0)
      RETURN
      END


! { dg-do run }
! This program tests that use associated common blocks work.
!
! provided by Paul Thomas - pault@gcc.gnu.org
!
module m1
  common /x/ a
end module m1
module m2
  common /x/ a
end module m2

subroutine foo ()
  use m2
  if (a.ne.99.0) call abort ()
end subroutine foo

program collision
  use m1
  use m2, only: b=>a
  b = 99.0
  call foo ()
end program collision


text:
      subroutine pinner(a,b,c,n,m,ndim)
      implicit none
C*DEC added by Y.Tange 10-Jan-1995
      integer*4 n,m,ndim,i,j
      real*8 a(ndim,m),b(m),c(n),s
C*DEC End
C*HP
C     real*8 a(ndim,m),b(m),c(n),s
C     integer*4 n,m,ndim,i,j
C*HP
      do 11 i=1,n
        s=0d0
        do 10 j=1,m
          s=s+a(i,j)*b(j)
   10   continue
        c(i)=s
   11 continue
      return
      end

max_stars_repo_name: noboruatkek/SAD

max_issues_repo_name: noboruatkek/SAD

max_forks_repo_name: noboruatkek/SAD


text:
      Subroutine screen(x,y,w)
      dimension xx(10), yy(10)
*
      w5 = 0.2*w
      call Igset('BORD',1.)
      call Igset('FAIS',0.)
      call Igset('LTYP',1.)
      Call Igbox(x,x+w,y+0.4*w,y+w)
      xx(1) = x
      xx(2) = x + w5
      xx(3) = x + 0.8*w
      xx(4) = x + w
      xx(5) = xx(1)
      yy(1) = y
      yy(2) = y + 0.2*w
      yy(3) = yy(2)
      yy(4) = y
      yy(5) = yy(1)
      call ipl(5,xx,yy)
      r = 0.66*w5
      call igarc(x+w5,y+4*w5,r,r,90.,180.)
      call igarc(x+w5,y+3*w5,r,r,180.,270.)
      call igarc(x+4*w5,y+4*w5,r,r,0.,90.)
      call igarc(x+4*w5,y+3*w5,r,r,270.,0.)
      xx(1) = x + w5
      xx(2) = x + 4*w5
      yy(1) = y +4.66*w5
      yy(2) = yy(1)
      call ipl(2,xx,yy)
      yy(1) = y +2.33*w5
      yy(2) = yy(1)
      call ipl(2,xx,yy)
      xx(1) = x + 0.33*w5
      xx(2) = xx(1)
      yy(1) = y + 3*w5
      yy(2) = y + 4*w5
      call ipl(2,xx,yy)
      xx(1) = x + 4.66*w5
      xx(2) = xx(1)
      call ipl(2,xx,yy)
*
      wb = 0.1*w5
      do j=0,3
         x1 = x + w5 + 0.5*wb - 2*j*wb
         y1 = y + w5 - (2.5+2*j)*wb
         do i=1,15+2*j
            call Igbox(x1,x1+wb,y1,y1+wb)
            x1 = x1 + 2*wb
         enddo
      enddo
      call igbox (x+2*w5,x+3*w5,y+w5,y+2*w5)
*
      End

max_stars_repo_name: berghaus/cernlib-docs

max_issues_repo_name: berghaus/cernlib-docs

max_forks_repo_name: berghaus/cernlib-docs


text:
	SUBROUTINE SNEOPU  ( snfile, timstn, isnfln, newfil, iret )
C************************************************************************
C* SNEOPU								*
C*									*
C* This subroutine opens an unmerged sounding dataset. 			*
C* If SNFILE does not exist, a new sounding file with the given 	*
C* parameters is created.						*
C*									*
C* SNEOPU  ( SNFILE, TIMSTN, ISNFLN, NEWFIL, IRET )			*
C*									*
C* Input parameters:							*
C*	SNFILE		CHAR*	 	Sounding file name		*
C*      TIMSTN		CHAR*           # Times, stations in new file   *
C*									*
C* Output parameters:							*
C*	ISNFLN		INTEGER		Sounding file number		*
C*	NEWFIL		LOGICAL		New file created flag		*
C* 	IRET		INTEGER		Return code			*
C*					  0 = normal return		*
C*					 -7 = error creating file	*
C*					 -8 = different parms		*
C*					 -9 = invalid vertical coord	*
C*					-10 = error opening file	*
C*					-20 = merged dataset		*
C*					-12 = NTIME is non positive	*
C*					-13 = NSTNS is non positive	*
C**									*
C* Log:									*
C* S. Schotz/GSC	12/89						*
C************************************************************************
	INCLUDE		'GEMPRM.PRM'
C*
	CHARACTER*(*)	snfile, timstn
	LOGICAL		newfil
C*
	CHARACTER	dsparm (MMPARM)*4
	LOGICAL		exist, mrgdat
C*
	INTEGER 	iarr (2)
	EQUIVALENCE 	( iarr (1), ntime ) , ( iarr (2), nstns )
C*
C-------------------------------------------------------------------------
	iret   = 0
	newfil = .false.
C
C*	First check to see if the file exists.
C
	INQUIRE  ( FILE = snfile, EXIST = exist )
C
C*	If the file exists, open and check parameters.
C
	IF  ( exist )  THEN
	    CALL SN_OPNF  ( snfile, .true., isnfln, ifl, mparms,
     +			    dsparm, jvert, mrgdat, iret )
            IF ( iret .ne. 0) THEN
		CALL ER_WMSG ( 'SN', iret, snfile, ier )
		iret = -10
		RETURN
	    ELSE IF  ( mrgdat )  THEN
		iret = -20
		RETURN
            END IF
C
C*	    Otherwise, create a new file.
C
        ELSE
	    CALL ST_ILST ( timstn, '/', 0, 2, iarr, n, ier )
C
C*	    Check for errors in ntime and nstns.
C
	    IF ( ntime .le. 0 ) THEN
		iret = -12
	    ELSE IF ( nstns .le. 0 ) THEN
		iret = -13
	    ELSE
                newfil = .true.
	        CALL SN_CRUA ( snfile, MFRAOB, 3, nstns, ntime,
     +                      .false., .true., .true., isnfln, iret )
                IF  ( iret .ne. 0 )  THEN
	           iret = -7
	        END IF
            END IF
	END IF
C*
	RETURN
	END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


text:
Derek Moore won UC Davis UC Davis first Division I DI wrestling championship in 2007.  He claimed his win in the 2007 NCAA Wrestling Finals for 141 lbs., not only beating the 1st seed but controlling the match completely and finishing with a technical fall (scoring more than 15 points over the opponent, considered the most dominating victory possible).

Media

  http://www.youtube.com/watch?vxvc2LZIgsgQ Video of the final match.

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
C Copyright(C) 1999-2020 National Technology & Engineering Solutions
C of Sandia, LLC (NTESS).  Under the terms of Contract DE-NA0003525 with
C NTESS, the U.S. Government retains certain rights in this software.
C
C See packages/seacas/LICENSE for details

C=======================================================================
      SUBROUTINE RWSTEP (NTXT, NDB, istep, idelb,
     &  NVARGL, NVARNP, NUMNP, NVAREL, NELBLK, NUMELB,
     &  ISEVOK, TIME, VAR, *)
C=======================================================================

C   --*** RDSTEP *** (TXTEXO) Read database variables for one time step
C   --   Written by Amy Gilkey - revised 03/02/88
C   --
C   --RDSTEP reads the database history, global, nodal, and element variables
C   --for one time step.
C   --
C   --Parameters:
C   --   NTXT - IN - the text file
C   --   NVARGL - IN - the number of global variables
C   --   NVARNP - IN - the number of nodal variables
C   --   NUMNP - IN - the number of nodes
C   --   NVAREL - IN - the number of element variables
C   --   NELBLK - IN - the number of element blocks
C   --   NUMELB - IN - the number of elements per block
C   --   ISEVOK - IN - the element block variable truth table;
C   --      variable i of block j exists iff ISEVOK(j,i)
C   --   TIME - OUT - the time step time
C   --   VAR - OUT - the global variables for the time step
C   --   * - return statement if error encountered, including end-of-file;
C   --      message is printed
C   --
C   --Database must be positioned in front of time step upon entry;
C   --upon exit positioned after time step.

      INTEGER NUMELB(*)
      integer idelb(*)
      LOGICAL ISEVOK(NVAREL, *)
      REAL VAR(*)

      CHARACTER*32 STRA

C   --Read step time

      READ (NTXT, *, END=180, ERR=130)
      READ (NTXT, *, END=180, ERR=130) TIME
      call exptim(ndb, istep, time, ierr)

C      --Read global variables

         IF (NVARGL .GT. 0) THEN
            READ (NTXT, *, END=150, ERR=150)
            READ (NTXT, *, END=150, ERR=150) (VAR(I), I=1,NVARGL)
            call expgv(ndb, istep, nvargl, var, ierr)
         END IF

C      --Read nodal variables

         IF (NVARNP .GT. 0) THEN
            READ (NTXT, *, END=160, ERR=160)
            DO 100 INP = 1, NUMNP
               READ (NTXT, *, END=160, ERR=160)
     &          (VAR(inp+(i-1)*numnp), I=1,NVARNP)
  100       CONTINUE

            ioff = 1
            do 105 i=1, nvarnp
              call expnv(ndb, istep, i, numnp, var(ioff), ierr)
              ioff = ioff + numnp
 105        continue
         END IF

C      --Read element variables
C ... Values for all element variables are stored in text file
C     Read an element block at a time and then write that element block
         IF (NVAREL .GT. 0) THEN
C ... Skip comment line
            READ (NTXT, *, END=170, ERR=170)
            DO 120 IELB = 1, NELBLK
C ... Skip comment line
              READ (NTXT, *, END=170, ERR=170)
              READ (NTXT, *, END=170, ERR=170)
               DO 110 N = 1, NUMELB(IELB)
                  READ (NTXT, *, END=170, ERR=170)
     &               (VAR(N+(I-1)*numelb(ielb)), I=1,NVAREL)
  110          CONTINUE
C ... Write this blocks variables
               ioff = 1
               do 115 i=1, nvarel
                 if (isevok(i,ielb)) then
                   call expev(ndb, istep, i, idelb(ielb), numelb(ielb),
     &               var(ioff), ierr)
                 end if
                 ioff = ioff + numelb(ielb)
 115           continue
  120       CONTINUE
         END IF

      RETURN

  130 CONTINUE
      CALL PRTERR ('FATAL', 'Reading TIME STEP TIME')
      GOTO 180
  150 CONTINUE
      CALL PRTERR ('FATAL', 'Reading GLOBAL VARIABLES')
      GOTO 180
  160 CONTINUE
      CALL INTSTR (1, 0, INP, STRA, LSTRA)
      CALL PRTERR ('FATAL',
     &   'Reading NODAL VARIABLES for node ' // STRA(:LSTRA))
      GOTO 180
  170 CONTINUE
      CALL INTSTR (1, 0, N, STRA, LSTRA)
      CALL PRTERR ('FATAL',
     &   'Reading ELEMENT VARIABLES for element ' // STRA(:LSTRA))
      GOTO 180
  180 CONTINUE
      RETURN 1
      END

max_stars_repo_name: jschueller/seacas

max_issues_repo_name: jschueller/seacas

max_forks_repo_name: jschueller/seacas


text: include "../constant.fppi"

ext: fpp

max_stars_repo_name: kevin-f-ortega/fpp

max_issues_repo_name: kevin-f-ortega/fpp

max_forks_repo_name: kevin-f-ortega/fpp


text:
      SUBROUTINE MB03TS( ISHAM, WANTU, N, A, LDA, G, LDG, U1, LDU1, U2,
     $                   LDU2, J1, N1, N2, DWORK, INFO )
C
C     SLICOT RELEASE 5.7.
C
C     Copyright (c) 2002-2020 NICONET e.V.
C
C     PURPOSE
C
C     To swap diagonal blocks A11 and A22 of order 1 or 2 in the upper
C     quasi-triangular matrix A contained in a skew-Hamiltonian matrix
C
C                   [  A   G  ]          T
C             X  =  [       T ],   G = -G,
C                   [  0   A  ]
C
C     or in a Hamiltonian matrix
C
C                   [  A   G  ]          T
C             X  =  [       T ],   G =  G.
C                   [  0  -A  ]
C
C     This routine is a modified version of the LAPACK subroutine
C     DLAEX2.
C
C     The matrix A must be in Schur canonical form (as returned by the
C     LAPACK routine DHSEQR), that is, block upper triangular with
C     1-by-1 and 2-by-2 diagonal blocks; each 2-by-2 diagonal block has
C     its diagonal elements equal and its off-diagonal elements of
C     opposite sign.
C
C     ARGUMENTS
C
C     Mode Parameters
C
C     ISHAM   LOGIGAL
C             Specifies the type of X:
C             = .TRUE.:   X is a Hamiltonian matrix;
C             = .FALSE.:  X is a skew-Hamiltonian matrix.
C
C     WANTU   LOGIGAL
C             = .TRUE.:   update the matrices U1 and U2 containing the
C                         Schur vectors;
C             = .FALSE.:  do not update U1 and U2.
C
C     Input/Output Parameters
C
C     N       (input) INTEGER
C             The order of the matrix A. N >= 0.
C
C     A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
C             On entry, the leading N-by-N part of this array must
C             contain the upper quasi-triangular matrix A, in Schur
C             canonical form.
C             On exit, the leading N-by-N part of this array contains
C             the reordered matrix A, again in Schur canonical form.
C
C     LDA     INTEGER
C             The leading dimension of the array A.  LDA >= MAX(1,N).
C
C     G       (input/output) DOUBLE PRECISION array, dimension (LDG,N)
C             On entry, the leading N-by-N part of this array must
C             contain the upper triangular part of the symmetric
C             matrix G, if ISHAM = .TRUE., or the strictly upper
C             triangular part of the skew-symmetric matrix G, otherwise.
C             The rest of this array is not referenced.
C             On exit, the leading N-by-N part of this array contains
C             the upper or strictly upper triangular part of the
C             symmetric or skew-symmetric matrix G, respectively,
C             updated by the orthogonal transformation which reorders A.
C
C     LDG     INTEGER
C             The leading dimension of the array G.  LDG >= MAX(1,N).
C
C     U1      (input/output) DOUBLE PRECISION array, dimension (LDU1,N)
C             On entry, if WANTU = .TRUE., the leading N-by-N part of
C             this array must contain the matrix U1.
C             On exit, if WANTU = .TRUE., the leading N-by-N part of
C             this array contains U1, postmultiplied by the orthogonal
C             transformation which reorders A. See the description in
C             the SLICOT subroutine MB03TD for further details.
C             If WANTU = .FALSE., this array is not referenced.
C
C     LDU1    INTEGER
C             The leading dimension of the array U1.
C             LDU1 >= MAX(1,N),  if WANTU = .TRUE.;
C             LDU1 >= 1,         otherwise.
C
C     U2      (input/output) DOUBLE PRECISION array, dimension (LDU2,N)
C             On entry, if WANTU = .TRUE., the leading N-by-N part of
C             this array must contain the matrix U2.
C             On exit, if WANTU = .TRUE., the leading N-by-N part of
C             this array contains U2, postmultiplied by the orthogonal
C             transformation which reorders A.
C             If WANTU = .FALSE., this array is not referenced.
C
C     LDU2    INTEGER
C             The leading dimension of the array U2.
C             LDU2 >= MAX(1,N),  if WANTU = .TRUE.;
C             LDU2 >= 1,         otherwise.
C
C     J1      (input) INTEGER
C             The index of the first row of the first block A11.
C             If J1+N1 < N, then A11 is swapped with the block starting
C             at (J1+N1+1)-th diagonal element.
C             If J1+N1 > N, then A11 is the last block in A and swapped
C             with -A11', if ISHAM = .TRUE.,
C             or    A11', if ISHAM = .FALSE..
C
C     N1      (input) INTEGER
C             The order of the first block A11. N1 = 0, 1 or 2.
C
C     N2      (input) INTEGER
C             The order of the second block A22. N2 = 0, 1 or 2.
C
C     Workspace
C
C     DWORK   DOUBLE PRECISION array, dimension (N)
C
C     Error Indicator
C
C     INFO    INTEGER
C             = 0:  successful exit;
C             = 1:  the transformed matrix A would be too far from Schur
C                   form; the blocks are not swapped and A, G, U1 and
C                   U2 are unchanged.
C
C     REFERENCES
C
C     [1] Bai, Z., and Demmel, J.W.
C        On swapping diagonal blocks in real Schur form.
C        Linear Algebra Appl., 186, pp. 73-95, 1993.
C
C     [2] Benner, P., Kressner, D., and Mehrmann, V.
C         Skew-Hamiltonian and Hamiltonian Eigenvalue Problems: Theory,
C         Algorithms and Applications. Techn. Report, TU Berlin, 2003.
C
C     CONTRIBUTORS
C
C     D. Kressner, Technical Univ. Berlin, Germany, and
C     P. Benner, Technical Univ. Chemnitz, Germany, December 2003.
C
C     REVISIONS
C
C     V. Sima, May 2008 (SLICOT version of the HAPACK routine DHAEX2).
C
C     ******************************************************************
C
C     .. Parameters ..
      DOUBLE PRECISION   ZERO, HALF, ONE, TWO, THIRTY, FORTY
      PARAMETER          ( ZERO  = 0.0D+0, HALF = 0.5D+0, ONE = 1.0D+0,
     $                     TWO   = 2.0D+0, THIRTY = 3.0D+1,
     $                     FORTY = 4.0D+1 )
      INTEGER            LDD, LDX
      PARAMETER          ( LDD = 4, LDX = 2 )
C     .. Scalar Arguments ..
      LOGICAL            ISHAM, WANTU
      INTEGER            INFO, J1, LDA, LDG, LDU1, LDU2, N, N1, N2
C     .. Array Arguments ..
      DOUBLE PRECISION   A(LDA,*), DWORK(*), G(LDG,*), U1(LDU1,*),
     $                   U2(LDU2,*)
C     .. Local Scalars ..
      LOGICAL            LBLK
      INTEGER            IERR, J2, J3, J4, K, ND
      DOUBLE PRECISION   A11, A22, A33, CS, DNORM, EPS, SCALE, SMLNUM,
     $                   SN, TAU, TAU1, TAU2, TEMP, THRESH, WI1, WI2,
     $                   WR1, WR2, XNORM
C     .. Local Arrays ..
      DOUBLE PRECISION   D(LDD,4), V(3), V1(3), V2(3), X(LDX,2)
C     .. External Functions ..
      DOUBLE PRECISION   DDOT, DLAMCH, DLANGE
      EXTERNAL           DDOT, DLAMCH, DLANGE
C     .. External Subroutines ..
      EXTERNAL           DAXPY, DLACPY, DLANV2, DLARFG, DLARFX, DLARTG,
     $                   DLASET, DLASY2, DROT, DSCAL, DSWAP, DSYMV,
     $                   DSYR2, MB01MD, MB01ND
C     .. Intrinsic Functions ..
      INTRINSIC          ABS, MAX
C
C     .. Executable Statements ..
C
      INFO = 0
C
C     Quick return if possible.
C
      IF ( N.EQ.0 .OR. N1.EQ.0 .OR. N2.EQ.0 )
     $   RETURN
      LBLK = ( J1+N1.GT.N )
C
      J2 = J1 + 1
      J3 = J1 + 2
      J4 = J1 + 3
C
      IF ( LBLK .AND. N1.EQ.1 ) THEN
C
         IF ( ISHAM ) THEN
            A11 = A(N,N)
            CALL DLARTG( G(N,N), -TWO*A11, CS, SN, TEMP )
            CALL DROT( N-1, A(1,N), 1, G(1,N), 1, CS, SN )
            A(N,N) = -A11
            IF ( WANTU )
     $         CALL DROT( N, U1(1,N), 1, U2(1,N), 1, CS, SN )
         ELSE
            CALL DSWAP( N-1, A(1,N), 1, G(1,N), 1 )
            CALL DSCAL( N-1, -ONE, A(1,N), 1 )
            IF ( WANTU ) THEN
               CALL DSWAP( N, U1(1,N), 1, U2(1,N), 1 )
               CALL DSCAL( N, -ONE, U1(1,N), 1 )
            END IF
         END IF
C
      ELSE IF ( LBLK .AND. N1.EQ.2 ) THEN
C
         IF ( ISHAM ) THEN
C
C           Reorder Hamiltonian matrix:
C
C                      [ A11  G11  ]
C                      [         T ].
C                      [  0  -A11  ]
C
            ND = 4
            CALL DLACPY( 'Full', 2, 2, A(N-1,N-1), LDA, D, LDD )
            CALL DLASET( 'All', 2, 2, ZERO, ZERO, D(3,1), LDD )
            CALL DLACPY( 'Upper', 2, 2, G(N-1,N-1), LDG, D(1,3), LDD )
            D(2,3) = D(1,4)
            D(3,3) = -D(1,1)
            D(4,3) = -D(1,2)
            D(3,4) = -D(2,1)
            D(4,4) = -D(2,2)
            DNORM = DLANGE( 'Max', ND, ND, D, LDD, DWORK )
C
C           Compute machine-dependent threshold for test for accepting
C           swap.
C
            EPS = DLAMCH( 'P' )
            SMLNUM = DLAMCH( 'S' ) / EPS
            THRESH = MAX( FORTY*EPS*DNORM, SMLNUM )
C
C           Solve A11*X + X*A11' = scale*G11 for X.
C
            CALL DLASY2( .FALSE., .FALSE., -1, 2, 2, D, LDD, D(3,3),
     $                   LDD, D(1,3), LDD, SCALE, X, LDX, XNORM, IERR )
C
C           Compute symplectic QR decomposition of
C
C                  (  -X11  -X12 )
C                  (  -X21  -X22 ).
C                  ( scale    0  )
C                  (    0  scale )
C
            TEMP = -X(1,1)
            CALL DLARTG( TEMP, SCALE, V1(1), V2(1), X(1,1) )
            CALL DLARTG( X(1,1), -X(2,1), V1(2), V2(2), TEMP )
            X(1,2) = -X(1,2)
            X(2,2) = -X(2,2)
            X(1,1) = ZERO
            X(2,1) = SCALE
            CALL DROT( 1, X(1,2), 1, X(1,1), 1, V1(1), V2(1) )
            CALL DROT( 1, X(1,2), 1, X(2,2), 1, V1(2), V2(2) )
            CALL DROT( 1, X(1,1), 1, X(2,1), 1, V1(2), V2(2) )
            CALL DLARTG( X(2,2), X(2,1), V1(3), V2(3), TEMP )
C
C           Perform swap provisionally on D.
C
            CALL DROT( 4, D(1,1), LDD, D(3,1), LDD, V1(1), V2(1) )
            CALL DROT( 4, D(1,1), LDD, D(2,1), LDD, V1(2), V2(2) )
            CALL DROT( 4, D(3,1), LDD, D(4,1), LDD, V1(2), V2(2) )
            CALL DROT( 4, D(2,1), LDD, D(4,1), LDD, V1(3), V2(3) )
            CALL DROT( 4, D(1,1), 1, D(1,3), 1, V1(1), V2(1) )
            CALL DROT( 4, D(1,1), 1, D(1,2), 1, V1(2), V2(2) )
            CALL DROT( 4, D(1,3), 1, D(1,4), 1, V1(2), V2(2) )
            CALL DROT( 4, D(1,2), 1, D(1,4), 1, V1(3), V2(3) )
C
C           Test whether to reject swap.
C
            IF ( MAX( ABS( D(3,1) ), ABS( D(3,2) ), ABS( D(4,1) ),
     $                ABS( D(4,2) ) ).GT.THRESH ) GO TO 50
C
            CALL DLACPY( 'All', 2, 2, D(1,1), LDD, A(N-1,N-1), LDA )
            CALL DLACPY( 'Upper', 2, 2, D(1,3), LDD, G(N-1,N-1), LDG )
C
            IF ( N.GT.2 ) THEN
               CALL DROT( N-2, A(1,N-1), 1, G(1,N-1), 1, V1(1), V2(1) )
               CALL DROT( N-2, A(1,N-1), 1, A(1,N), 1, V1(2), V2(2) )
               CALL DROT( N-2, G(1,N-1), 1, G(1,N), 1, V1(2), V2(2) )
               CALL DROT( N-2, A(1,N), 1, G(1,N), 1, V1(3), V2(3) )
            END IF
C
            IF ( WANTU ) THEN
               CALL DROT( N, U1(1,N-1), 1, U2(1,N-1), 1, V1(1), V2(1) )
               CALL DROT( N, U1(1,N-1), 1, U1(1,N), 1, V1(2), V2(2) )
               CALL DROT( N, U2(1,N-1), 1, U2(1,N), 1, V1(2), V2(2) )
               CALL DROT( N, U1(1,N), 1, U2(1,N), 1, V1(3), V2(3) )
            END IF
C
         ELSE
C
            IF ( ABS( A(N-1,N) ).GT.ABS( A(N,N-1) ) ) THEN
               TEMP = G(N-1,N)
               CALL DLARTG( TEMP, A(N-1,N), CS, SN, G(N-1,N) )
               SN = -SN
               CALL DROT(N-2, A(1,N), 1, G(1,N), 1, CS, SN )
C
               A(N-1,N) = -SN*A(N,N-1)
               TEMP = -CS*A(N,N-1)
               A(N,N-1) = G(N-1,N)
               G(N-1,N) = TEMP
               IF ( WANTU )
     $            CALL DROT( N, U1(1,N), 1, U2(1,N), 1, CS, SN )
               CALL DSWAP( N-2, A(1,N-1), 1, G(1,N-1), 1 )
               CALL DSCAL( N-2, -ONE, A(1,N-1), 1 )
               IF ( WANTU ) THEN
                  CALL DSWAP( N, U1(1,N-1), 1, U2(1,N-1), 1 )
                  CALL DSCAL( N, -ONE, U1(1,N-1), 1 )
               END IF
            ELSE
               TEMP = G(N-1,N)
               CALL DLARTG( TEMP, A(N,N-1), CS, SN, G(N-1,N) )
               CALL DROT( N-2, A(1,N-1), 1, G(1,N-1), 1, CS, SN )
               A(N,N-1) = -SN*A(N-1,N)
               A(N-1,N) = CS*A(N-1,N)
               IF ( WANTU )
     $            CALL DROT( N, U1(1,N-1), 1, U2(1,N-1), 1, CS, SN )
               CALL DSWAP( N-1, A(1,N), 1, G(1,N), 1 )
               CALL DSCAL( N-1, -ONE, A(1,N), 1 )
               IF ( WANTU ) THEN
                  CALL DSWAP( N, U1(1,N), 1, U2(1,N), 1 )
                  CALL DSCAL( N, -ONE, U1(1,N), 1 )
               END IF
            END IF
         END IF
C
C        Standardize new 2-by-2 block.
C
         CALL DLANV2( A(N-1,N-1), A(N-1,N), A(N,N-1),
     $                A(N,N), WR1, WI1, WR2, WI2, CS, SN )
         CALL DROT( N-2, A(1,N-1), 1, A(1,N), 1, CS, SN )
         IF ( ISHAM ) THEN
            TEMP = G(N-1,N)
            CALL DROT( N-1, G(1,N-1), 1, G(1,N), 1, CS, SN )
            TAU = CS*TEMP + SN*G(N,N)
            G(N,N) = CS*G(N,N) - SN*TEMP
            G(N-1,N-1) = CS*G(N-1,N-1) + SN*TAU
            CALL DROT( 1, G(N-1,N), LDG, G(N,N), LDG, CS, SN )
         ELSE
            CALL DROT( N-2, G(1,N-1), 1, G(1,N), 1, CS, SN )
         END IF
         IF ( WANTU ) THEN
            CALL DROT( N, U1(1,N-1), 1, U1(1,N), 1, CS, SN )
            CALL DROT( N, U2(1,N-1), 1, U2(1,N), 1, CS, SN )
         END IF
C
      ELSE IF ( N1.EQ.1 .AND. N2.EQ.1 ) THEN
C
C        Swap two 1-by-1 blocks.
C
         A11 = A(J1,J1)
         A22 = A(J2,J2)
C
C        Determine the transformation to perform the interchange.
C
         CALL DLARTG( A(J1,J2), A22-A11, CS, SN, TEMP )
C
C        Apply transformation to the matrix A.
C
         IF ( J3.LE.N )
     $      CALL DROT( N-J1-1, A(J1,J3), LDA, A(J2,J3), LDA, CS, SN )
         CALL DROT( J1-1, A(1,J1), 1, A(1,J2), 1, CS, SN )
C
         A(J1,J1) = A22
         A(J2,J2) = A11
C
C        Apply transformation to the matrix G.
C
         IF ( ISHAM ) THEN
            TEMP = G(J1,J2)
            CALL DROT( J1, G(1,J1), 1, G(1,J2), 1, CS, SN )
            TAU = CS*TEMP + SN*G(J2,J2)
            G(J2,J2) = CS*G(J2,J2) - SN*TEMP
            G(J1,J1) = CS*G(J1,J1) + SN*TAU
            CALL DROT( N-J1, G(J1,J2), LDG, G(J2,J2), LDG, CS, SN )
         ELSE
            IF ( N.GT.J1+1 )
     $         CALL DROT( N-J1-1, G(J1,J1+2), LDG, G(J2,J1+2), LDG, CS,
     $                    SN )
            CALL DROT( J1-1, G(1,J1), 1, G(1,J2), 1, CS, SN )
         END IF
         IF ( WANTU ) THEN
C
C           Accumulate transformation in the matrices U1 and U2.
C
            CALL DROT( N, U1(1,J1), 1, U1(1,J2), 1, CS, SN )
            CALL DROT( N, U2(1,J1), 1, U2(1,J2), 1, CS, SN )
         END IF
C
      ELSE
C
C        Swapping involves at least one 2-by-2 block.
C
C        Copy the diagonal block of order N1+N2 to the local array D
C        and compute its norm.
C
         ND = N1 + N2
         CALL DLACPY( 'Full', ND, ND, A(J1,J1), LDA, D, LDD )
         DNORM = DLANGE( 'Max', ND, ND, D, LDD, DWORK )
C
C        Compute machine-dependent threshold for test for accepting
C        swap.
C
         EPS = DLAMCH( 'P' )
         SMLNUM = DLAMCH( 'S' ) / EPS
         THRESH = MAX( THIRTY*EPS*DNORM, SMLNUM )
C
C        Solve A11*X - X*A22 = scale*A12 for X.
C
         CALL DLASY2( .FALSE., .FALSE., -1, N1, N2, D, LDD,
     $                D(N1+1,N1+1), LDD, D(1,N1+1), LDD, SCALE, X, LDX,
     $                XNORM, IERR )
C
C        Swap the adjacent diagonal blocks.
C
         K = N1 + N1 + N2 - 3
         GO TO ( 10, 20, 30 )K
C
   10    CONTINUE
C
C        N1 = 1, N2 = 2: generate elementary reflector H so that:
C
C        ( scale, X11, X12 ) H = ( 0, 0, * ).
C
         V(1) = SCALE
         V(2) = X(1,1)
         V(3) = X(1,2)
         CALL DLARFG( 3, V(3), V, 1, TAU )
         V(3) = ONE
         A11 = A(J1,J1)
C
C        Perform swap provisionally on diagonal block in D.
C
         CALL DLARFX( 'Left', 3, 3, V, TAU, D, LDD, DWORK )
         CALL DLARFX( 'Right', 3, 3, V, TAU, D, LDD, DWORK )
C
C        Test whether to reject swap.
C
         IF ( MAX( ABS( D(3,1) ), ABS( D(3,2) ), ABS( D(3,3)-A11 ) )
     $        .GT.THRESH ) GO TO 50
C
C        Accept swap: apply transformation to the entire matrix A.
C
         CALL DLARFX( 'Left', 3, N-J1+1, V, TAU, A(J1,J1), LDA, DWORK )
         CALL DLARFX( 'Right', J2, 3, V, TAU, A(1,J1), LDA, DWORK )
C
         A(J3,J1) = ZERO
         A(J3,J2) = ZERO
         A(J3,J3) = A11
C
C        Apply transformation to G.
C
         IF ( ISHAM ) THEN
            CALL DLARFX( 'Right', J1-1, 3, V, TAU, G(1,J1), LDG, DWORK )
            CALL DSYMV( 'Upper', 3, TAU, G(J1,J1), LDG, V, 1, ZERO,
     $                  DWORK, 1 )
            TEMP = -HALF*TAU*DDOT( 3, DWORK, 1, V, 1 )
            CALL DAXPY( 3, TEMP, V, 1, DWORK, 1 )
            CALL DSYR2( 'Upper', 3, -ONE, V, 1, DWORK, 1,
     $                  G(J1,J1), LDG )
            IF ( N.GT.J1+2 )
     $         CALL DLARFX( 'Left', 3, N-J1-2, V, TAU, G(J1,J1+3), LDG,
     $                      DWORK )
         ELSE
            CALL DLARFX( 'Right', J1-1, 3, V, TAU, G(1,J1), LDG, DWORK )
            CALL MB01MD( 'Upper', 3, TAU, G(J1,J1), LDG, V, 1, ZERO,
     $                   DWORK, 1 )
            CALL MB01ND( 'Upper', 3, ONE, V, 1, DWORK, 1, G(J1,J1),
     $                   LDG )
            IF ( N.GT.J1+2 )
     $         CALL DLARFX( 'Left', 3, N-J1-2, V, TAU, G(J1,J1+3), LDG,
     $                      DWORK )
         END IF
C
         IF ( WANTU ) THEN
C
C           Accumulate transformation in the matrices U1 and U2.
C
            CALL DLARFX( 'R', N, 3, V, TAU, U1(1,J1), LDU1, DWORK )
            CALL DLARFX( 'R', N, 3, V, TAU, U2(1,J1), LDU2, DWORK )
         END IF
         GO TO 40
C
   20    CONTINUE
C
C        N1 = 2, N2 = 1: generate elementary reflector H so that:
C
C        H (  -X11 ) = ( * )
C          (  -X21 ) = ( 0 ).
C          ( scale ) = ( 0 )
C
         V(1) = -X(1,1)
         V(2) = -X(2,1)
         V(3) = SCALE
         CALL DLARFG( 3, V(1), V(2), 1, TAU )
         V(1) = ONE
         A33 = A(J3,J3)
C
C        Perform swap provisionally on diagonal block in D.
C
         CALL DLARFX( 'L', 3, 3, V, TAU, D, LDD, DWORK )
         CALL DLARFX( 'R', 3, 3, V, TAU, D, LDD, DWORK )
C
C        Test whether to reject swap.
C
         IF ( MAX( ABS( D(2,1) ), ABS( D(3,1) ), ABS( D(1,1)-A33 ) )
     $       .GT. THRESH ) GO TO 50
C
C        Accept swap: apply transformation to the entire matrix A.
C
         CALL DLARFX( 'Right', J3, 3, V, TAU, A(1,J1), LDA, DWORK )
         CALL DLARFX( 'Left', 3, N-J1, V, TAU, A(J1,J2), LDA, DWORK )
C
         A(J1,J1) = A33
         A(J2,J1) = ZERO
         A(J3,J1) = ZERO
C
C        Apply transformation to G.
C
         IF ( ISHAM ) THEN
            CALL DLARFX( 'Right', J1-1, 3, V, TAU, G(1,J1), LDG, DWORK )
            CALL DSYMV( 'Upper', 3, TAU, G(J1,J1), LDG, V, 1, ZERO,
     $                  DWORK, 1 )
            TEMP = -HALF*TAU*DDOT( 3, DWORK, 1, V, 1 )
            CALL DAXPY( 3, TEMP, V, 1, DWORK, 1 )
            CALL DSYR2( 'Upper', 3, -ONE, V, 1, DWORK, 1, G(J1,J1),
     $                  LDG )
            IF ( N.GT.J1+2 )
     $         CALL DLARFX( 'Left', 3, N-J1-2, V, TAU, G(J1,J1+3), LDG,
     $                      DWORK )
         ELSE
            CALL DLARFX( 'Right', J1-1, 3, V, TAU, G(1,J1), LDG, DWORK )
            CALL MB01MD( 'Upper', 3, TAU, G(J1,J1), LDG, V, 1, ZERO,
     $                   DWORK, 1 )
            CALL MB01ND( 'Upper', 3, ONE, V, 1, DWORK, 1, G(J1,J1),
     $                   LDG )
            IF ( N.GT.J1+2 )
     $         CALL DLARFX( 'Left', 3, N-J1-2, V, TAU, G(J1,J1+3), LDG,
     $                      DWORK )
         END IF
C
         IF ( WANTU ) THEN
C
C           Accumulate transformation in the matrices U1 and U2.
C
            CALL DLARFX( 'R', N, 3, V, TAU, U1(1,J1), LDU1, DWORK )
            CALL DLARFX( 'R', N, 3, V, TAU, U2(1,J1), LDU2, DWORK )
         END IF
         GO TO 40
C
   30    CONTINUE
C
C        N1 = 2, N2 = 2: generate elementary reflectors H(1) and H(2) so
C        that:
C
C        H(2) H(1) (  -X11  -X12 ) = (  *  * )
C                  (  -X21  -X22 )   (  0  * ).
C                  ( scale    0  )   (  0  0 )
C                  (    0  scale )   (  0  0 )
C
         V1(1) = -X(1,1)
         V1(2) = -X(2,1)
         V1(3) = SCALE
         CALL DLARFG( 3, V1(1), V1(2), 1, TAU1 )
         V1(1) = ONE
C
         TEMP  = -TAU1*( X(1,2)+V1(2)*X(2,2) )
         V2(1) = -TEMP*V1(2) - X(2,2)
         V2(2) = -TEMP*V1(3)
         V2(3) = SCALE
         CALL DLARFG( 3, V2(1), V2(2), 1, TAU2 )
         V2(1) = ONE
C
C        Perform swap provisionally on diagonal block in D.
C
         CALL DLARFX( 'L', 3, 4, V1, TAU1, D, LDD, DWORK )
         CALL DLARFX( 'R', 4, 3, V1, TAU1, D, LDD, DWORK )
         CALL DLARFX( 'L', 3, 4, V2, TAU2, D(2,1), LDD, DWORK )
         CALL DLARFX( 'R', 4, 3, V2, TAU2, D(1,2), LDD, DWORK )
C
C        Test whether to reject swap.
C
         IF ( MAX( ABS( D(3,1) ), ABS( D(3,2) ), ABS( D(4,1) ),
     $             ABS( D(4,2) ) ).GT.THRESH ) GO TO 50
C
C        Accept swap: apply transformation to the entire matrix A.
C
         CALL DLARFX( 'L', 3, N-J1+1, V1, TAU1, A(J1,J1), LDA, DWORK )
         CALL DLARFX( 'R', J4, 3, V1, TAU1, A(1,J1), LDA, DWORK )
         CALL DLARFX( 'L', 3, N-J1+1, V2, TAU2, A(J2,J1), LDA, DWORK )
         CALL DLARFX( 'R', J4, 3, V2, TAU2, A(1,J2), LDA, DWORK )
C
         A(J3,J1) = ZERO
         A(J3,J2) = ZERO
         A(J4,J1) = ZERO
         A(J4,J2) = ZERO
C
C        Apply transformation to G.
C
         IF ( ISHAM ) THEN
            CALL DLARFX( 'Right', J1-1, 3, V1, TAU1, G(1,J1), LDG,
     $                   DWORK )
            CALL DSYMV( 'Upper', 3, TAU1, G(J1,J1), LDG, V1, 1, ZERO,
     $                  DWORK, 1 )
            TEMP = -HALF*TAU1*DDOT( 3, DWORK, 1, V1, 1 )
            CALL DAXPY( 3, TEMP, V1, 1, DWORK, 1 )
            CALL DSYR2( 'Upper', 3, -ONE, V1, 1, DWORK, 1,
     $                  G(J1,J1), LDG )
            IF ( N.GT.J1+2 )
     $         CALL DLARFX( 'Left', 3, N-J1-2, V1, TAU1, G(J1,J1+3),
     $                      LDG, DWORK )
C
            CALL DLARFX( 'Right', J2-1, 3, V2, TAU2, G(1,J2), LDG,
     $                   DWORK )
            CALL DSYMV( 'Upper', 3, TAU2, G(J2,J2), LDG, V2, 1, ZERO,
     $                  DWORK, 1 )
            TEMP = -HALF*TAU2*DDOT( 3, DWORK, 1, V2, 1 )
            CALL DAXPY( 3, TEMP, V2, 1, DWORK, 1 )
            CALL DSYR2( 'Upper', 3, -ONE, V2, 1, DWORK, 1, G(J2,J2),
     $                  LDG )
            IF ( N.GT.J2+2 )
     $         CALL DLARFX( 'Left', 3, N-J2-2, V2, TAU2, G(J2,J2+3),
     $                      LDG, DWORK )
         ELSE
            CALL DLARFX( 'Right', J1-1, 3, V1, TAU1, G(1,J1), LDG,
     $                   DWORK )
            CALL MB01MD( 'Upper', 3, TAU1, G(J1,J1), LDG, V1, 1, ZERO,
     $                   DWORK, 1 )
            CALL MB01ND( 'Upper', 3, ONE, V1, 1, DWORK, 1, G(J1,J1),
     $                   LDG )
            IF ( N.GT.J1+2 )
     $         CALL DLARFX( 'Left', 3, N-J1-2, V1, TAU1, G(J1,J1+3),
     $                      LDG, DWORK )
            CALL DLARFX( 'Right', J2-1, 3, V2, TAU2, G(1,J2), LDG,
     $                   DWORK )
            CALL MB01MD( 'Upper', 3, TAU2, G(J2,J2), LDG, V2, 1, ZERO,
     $                   DWORK, 1 )
            CALL MB01ND( 'Upper', 3, ONE, V2, 1, DWORK, 1, G(J2,J2),
     $                   LDG )
            IF ( N.GT.J2+2 )
     $         CALL DLARFX( 'Left', 3, N-J2-2, V2, TAU2, G(J2,J2+3),
     $                      LDG, DWORK )
         END IF
C
         IF ( WANTU ) THEN
C
C           Accumulate transformation in the matrices U1 and U2.
C
            CALL DLARFX( 'R', N, 3, V1, TAU1, U1(1,J1), LDU1, DWORK )
            CALL DLARFX( 'R', N, 3, V2, TAU2, U1(1,J2), LDU1, DWORK )
            CALL DLARFX( 'R', N, 3, V1, TAU1, U2(1,J1), LDU2, DWORK )
            CALL DLARFX( 'R', N, 3, V2, TAU2, U2(1,J2), LDU2, DWORK )
         END IF
C
   40    CONTINUE
C
         IF ( N2.EQ.2 ) THEN
C
C           Standardize new 2-by-2 block A11.
C
            CALL DLANV2( A(J1,J1), A(J1,J2), A(J2,J1), A(J2,J2), WR1,
     $                   WI1, WR2, WI2, CS, SN )
            CALL DROT( N-J1-1, A(J1,J1+2), LDA, A(J2,J1+2), LDA, CS,
     $                 SN )
            CALL DROT( J1-1, A(1,J1), 1, A(1,J2), 1, CS, SN )
            IF ( ISHAM ) THEN
               TEMP = G(J1,J2)
               CALL DROT( J1, G(1,J1), 1, G(1,J2), 1, CS, SN )
               TAU = CS*TEMP + SN*G(J2,J2)
               G(J2,J2) = CS*G(J2,J2) - SN*TEMP
               G(J1,J1) = CS*G(J1,J1) + SN*TAU
               CALL DROT( N-J1, G(J1,J2), LDG, G(J2,J2), LDG, CS, SN )
            ELSE
               IF ( N.GT.J1+1 )
     $            CALL DROT( N-J1-1, G(J1,J1+2), LDG, G(J2,J1+2), LDG,
     $                       CS, SN )
               CALL DROT( J1-1, G(1,J1), 1, G(1,J2), 1, CS, SN )
            END IF
            IF ( WANTU ) THEN
               CALL DROT( N, U1(1,J1), 1, U1(1,J2), 1, CS, SN )
               CALL DROT( N, U2(1,J1), 1, U2(1,J2), 1, CS, SN )
            END IF
         END IF
C
         IF ( N1.EQ.2 ) THEN
C
C           Standardize new 2-by-2 block A22.
C
            J3 = J1 + N2
            J4 = J3 + 1
            CALL DLANV2( A(J3,J3), A(J3,J4), A(J4,J3), A(J4,J4), WR1,
     $                   WI1, WR2, WI2, CS, SN )
            IF ( J3+2.LE.N )
     $         CALL DROT( N-J3-1, A(J3,J3+2), LDA, A(J4,J3+2), LDA, CS,
     $                    SN )
            CALL DROT( J3-1, A(1,J3), 1, A(1,J4), 1, CS, SN )
            IF ( ISHAM ) THEN
               TEMP = G(J3,J4)
               CALL DROT( J3, G(1,J3), 1, G(1,J4), 1, CS, SN )
               TAU = CS*TEMP + SN*G(J4,J4)
               G(J4,J4) = CS*G(J4,J4) - SN*TEMP
               G(J3,J3) = CS*G(J3,J3) + SN*TAU
               CALL DROT( N-J3, G(J3,J4), LDG, G(J4,J4), LDG, CS, SN )
            ELSE
               IF ( N.GT.J3+1 )
     $            CALL DROT( N-J3-1, G(J3,J3+2), LDG, G(J4,J3+2), LDG,
     $                       CS, SN )
               CALL DROT( J3-1, G(1,J3), 1, G(1,J4), 1, CS, SN )
            END IF
            IF ( WANTU ) THEN
               CALL DROT( N, U1(1,J3), 1, U1(1,J4), 1, CS, SN )
               CALL DROT( N, U2(1,J3), 1, U2(1,J4), 1, CS, SN )
            END IF
         END IF
C
      END IF
      RETURN
C
C     Exit with INFO = 1 if swap was rejected.
C
   50 CONTINUE
      INFO = 1
      RETURN
C *** Last line of MB03TS ***
      END

max_stars_repo_name: bnavigator/SLICOT-Reference

max_issues_repo_name: RJHKnight/slicotr

max_forks_repo_name: RJHKnight/slicotr


text:
      SUBROUTINE harfun(jacfac, hiota, gpsi, ipsi, js, nznt, xlt,
     1   xlz, xl, wt, wz, w, uboz, vboz, xjac)
      USE stel_kinds
      IMPLICIT NONE
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER :: js, nznt
      REAL(rprec) :: jacfac
      REAL(rprec), DIMENSION(*) :: hiota, gpsi, ipsi
      REAL(rprec), DIMENSION(nznt) ::
     1   xl, xlt, xlz, w, wt, wz, uboz, vboz, xjac
C-----------------------------------------------
C   L o c a l   P a r a m e t e r s
C-----------------------------------------------
      REAL(rprec), PARAMETER :: one = 1, zero = 0
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      REAL(rprec), ALLOCATABLE, DIMENSION(:) ::
     1   bsupu, bsupv, psubu, psubv
      REAL(rprec) :: dem, dem2, gpsi1, hiota1, ipsi1
	INTEGER     :: istat
C-----------------------------------------------
!
!     HERE, W IS THE PART OF THE TRANSFORMATION FUNCTION P
!     WHICH DEPENDS ON THE SURFACE-AVERAGED COVARIANT B COMPONENTS [RIGHT-SIDE IN Eq(10)],
!     COMPUTED IN transpmn AND vcoord_w.
!
!     THE FULL TRANSFORMATION FUNCTION IS THEN [Eq(10 in paper], with D=gpsi+iota*Ipsi:
!
!     P = (w - Ipsi*Lambda) / D
!
!     ALSO [Eq.(3) in paper]:
!
!     uboz = lambda + iota*P   (non-secular piece of boozer theta in VMEC coordinates)
!
!          = (gpsi*lambda + iota*w)/D
!
!     vboz = P                 (non-secular piece of boozer phi in VMEC coordinates)
!
!     FINALLY, XJAC IS THE JACOBIAN BETWEEN BOOZER, VMEC COORDINATES
!
!
!     NOTE THAT LAMBDA = xl, d(lambda)/du = xlt, d(lambda)/dv = xlz
!
      ALLOCATE(bsupu(nznt), bsupv(nznt), psubu(nznt), psubv(nznt),
     1         stat=istat)
	IF (istat .ne. 0) STOP 'Allocation error in harfun!'

      jacfac = gpsi(js) + hiota(js)*ipsi(js)
      IF (jacfac .eq. zero)
     1   STOP 'Boozer coordinate XFORM failed, jacfac = 0!'
      dem = one/jacfac
      gpsi1 = gpsi(js)*dem
      hiota1 = hiota(js)*dem
      ipsi1 = ipsi(js)*dem

      vboz = dem*w - ipsi1*xl        !TOTAL p in Eq.(10)
      uboz = xl + hiota(js)*vboz
      psubv = dem*wz - ipsi1*xlz
      psubu = dem*wt - ipsi1*xlt
      bsupv = 1 + xlt
      bsupu = hiota(js) - xlz

!
!     Eq. (12)
!
      xjac = bsupv*(1+psubv) + bsupu*psubu

      dem = MINVAL(xjac)
      dem2= MAXVAL(xjac)
!     SAL 07/06/16
!      IF (dem*dem2 .le. zero) PRINT *,
!     1   ' Jacobian xjac changed sign in harfun in xbooz_xform'

      DEALLOCATE(bsupu, bsupv, psubu, psubv, stat=istat)

      END SUBROUTINE harfun

      SUBROUTINE modbooz(bmnc, bmns, bmod, xmb, xnb,
     1   u_b, v_b, cosmm, sinmm, cosnn, sinnn,
     2   mnmax, mboz, nboz, nfp, lasym)
      USE stel_kinds
      IMPLICIT NONE
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      INTEGER :: mnmax, mboz, nboz, nfp
      REAL(rprec), DIMENSION(mnmax), INTENT(in) ::
     1   bmnc, bmns
      REAL(rprec), DIMENSION(mnmax), INTENT(in) :: xmb, xnb
      REAL(rprec), DIMENSION(4), INTENT(in)  :: u_b, v_b
      REAL(rprec), DIMENSION(4), INTENT(out) :: bmod
      REAL(rprec), DIMENSION(0:mboz) :: cosmm, sinmm
      REAL(rprec), DIMENSION(0:nboz) :: cosnn, sinnn
      LOGICAL, INTENT(in) :: lasym
!-----------------------------------------------
!   L o c a l   P a r a m e t e r s
!-----------------------------------------------
      REAL(rprec), PARAMETER :: one = 1.0_dp, zero = 0.0_dp
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      INTEGER :: mn, m, n, angles
      REAL(rprec) :: cost, sint
      REAL(rprec) :: sgn
!-----------------------------------------------
      bmod = zero

      ANGLE_LOOP: DO angles=1,4

      cosmm(0) = one
      sinmm(0) = zero
      cosmm(1) = COS(u_b(angles))
      sinmm(1) = SIN(u_b(angles))

      cosnn(0) = one
      sinnn(0) = zero
      IF (nboz .ge. 1) THEN
      cosnn(1) = COS(v_b(angles)*nfp)
      sinnn(1) = SIN(v_b(angles)*nfp)
      END IF

      DO m = 2,mboz
         cosmm(m) = cosmm(m-1)*cosmm(1)
     1            - sinmm(m-1)*sinmm(1)
         sinmm(m) = sinmm(m-1)*cosmm(1)
     1            + cosmm(m-1)*sinmm(1)
      END DO

      DO n = 2,nboz
         cosnn(n) = cosnn(n-1)*cosnn(1)
     1            - sinnn(n-1)*sinnn(1)
         sinnn(n) = sinnn(n-1)*cosnn(1)
     1            + cosnn(n-1)*sinnn(1)
      END DO

      DO mn=1,mnmax
        m = NINT(xmb(mn))
        n = NINT(ABS(xnb(mn)))/nfp
        sgn = SIGN(one,xnb(mn))
        cost = cosmm(m)*cosnn(n)
     1       + sinmm(m)*sinnn(n)*sgn
        bmod(angles) = bmod(angles) + bmnc(mn)*cost
        IF (lasym) THEN
        sint = sinmm(m)*cosnn(n)
     1       - cosmm(m)*sinnn(n)*sgn
        bmod(angles) = bmod(angles) + bmns(mn)*sint
        END IF
      END DO

      END DO ANGLE_LOOP

      END SUBROUTINE modbooz

max_stars_repo_name: joseluisvelasco/STELLOPT

max_issues_repo_name: joseluisvelasco/STELLOPT

max_forks_repo_name: joseluisvelasco/STELLOPT


! { dg-do compile }
!
! PR fortran/43591
!
! Pureness check for TPB/PPC in specification expressions
!
! Based on a test case of Thorsten Ohl
!
!

module m
  implicit none
  type t
     procedure(p1_type), nopass, pointer :: p1 => NULL()
  contains
     procedure, nopass :: tbp => p1_type
  end type t
contains
  subroutine proc (t1, t2)
    type(t), intent(in) :: t1, t2
    integer, dimension(t1%p1(), t2%tbp()) :: table
  end subroutine proc
  pure function p1_type()
   integer :: p1_type
   p1_type = 42
  end function p1_type
  pure subroutine p(t1)
    type(t), intent(inout) :: t1
    integer :: a(t1%p1())
  end subroutine p
end module m

module m2
  implicit none
  type t
     procedure(p1_type), nopass, pointer :: p1 => NULL()
  contains
     procedure, nopass :: tbp => p1_type
  end type t
contains
  subroutine proc (t1, t2)
    type(t), intent(in) :: t1, t2
    integer, dimension(t1%p1()) :: table1 ! { dg-error "must be PURE" }
    integer, dimension(t2%tbp()) :: table2 ! { dg-error "must be PURE" }
  end subroutine proc
  function p1_type()
    integer :: p1_type
    p1_type = 42
  end function p1_type
end module m2


text:
!     Fix bad quality quads based on angle ratio criterion (like
!     xmgredit5).
!     Inputs: screen; hgrid.gr3 (projection or lon/lat)
!     Output: hgrid.gr3.new (node order not changed); split_loc.bp (location of quads that r
!     split)

!     ifort -O2 -CB -g -traceback -o fix_bad_quads fix_bad_quads.f90

      implicit real*8(a-h,o-z)
      integer :: nwild(3),nwild2(4)
      real*8 :: aspect_ratio(2,2)
      integer, allocatable :: i34(:),elnode(:,:),i34_new(:),elnode_new(:,:)
      real*8, allocatable :: xnd(:),ynd(:),dp(:),area(:),quad_loc(:,:), &
     &xctr(:),yctr(:)

      pi=acos(-1.d0)
      !Quality is measured by the ratio of smallest to largest interior
      !angle of quads
      print*, 'Input threshold for angle ratio (e.g. 0.5):'
      read*, rat_angle

      open(14,file='hgrid.gr3',status='old')
      read(14,*)
      read(14,*) ne,np
      allocate(xnd(np),ynd(np),dp(np),area(ne),i34(ne),elnode(4,ne), &
     &i34_new(2*ne),elnode_new(4,2*ne),quad_loc(2,ne),xctr(ne),yctr(ne))

      do i=1,np
        read(14,*) j,xnd(i),ynd(i),dp(i)
      enddo !i

      do i=1,ne
        read(14,*) j,i34(i),elnode(1:i34(i),i)

        n1=elnode(1,i)
        n2=elnode(2,i)
        n3=elnode(3,i)
        area(i)=signa(xnd(n1),xnd(n2),xnd(n3),ynd(n1),ynd(n2),ynd(n3))
        if(area(i)<=0) then
          write(*,*)'Negative area at elem:',i
          stop
        endif

        if(i34(i)==4) then
          n4=elnode(4,i)
          tmp=signa(xnd(n1),xnd(n3),xnd(n4),ynd(n1),ynd(n3),ynd(n4))
          if(tmp<=0) then
            write(*,*)'Negative area at elem:',i
            stop
          endif
          area(i)=area(i)+tmp
        endif

        xctr(i)=sum(xnd(elnode(1:i34(i),i)))/i34(i)
        yctr(i)=sum(ynd(elnode(1:i34(i),i)))/i34(i)
      enddo !i=1,ne
      i34_new(1:ne)=i34
      elnode_new(:,1:ne)=elnode(:,1:ne)

      ne_extra=0
      do i=1,ne
        if(i34(i)/=4) cycle

        !Quads; interior angles
        ang_max=-1; ang_min=pi*1.1
        !local index where splitting occurs; quad is split into tri
        !(isplit,isplit+1,isplit+2)
        isplit=0
        do j=1,4
          n1=elnode(j,i) !angle centered @ n1
          j1=j+1; j3=j+3
          if(j1>4) j1=j1-4
          if(j3>4) j3=j3-4
          n2=elnode(j1,i)
          n3=elnode(j3,i)

          ar2=signa(xnd(n1),xnd(n2),xnd(n3),ynd(n1),ynd(n2),ynd(n3))
          if(ar2<=0) then !concave quad; split from this node
            isplit=j; exit
          endif

          !Convext @ n1; interior angle
          dot=(xnd(n2)-xnd(n1))*(xnd(n3)-xnd(n1))+(ynd(n2)-ynd(n1))*(ynd(n3)-ynd(n1))
          rl1=sqrt((xnd(n2)-xnd(n1))**2+(ynd(n2)-ynd(n1))**2)
          rl3=sqrt((xnd(n3)-xnd(n1))**2+(ynd(n3)-ynd(n1))**2)
          if(rl1==0.or.rl3==0) stop '0 side'
          costmp=dot/rl1/rl3
          if(abs(costmp)>1) then
            print*, 'Impossible cosine:',i,j,costmp
            stop
          endif
          ang=acos(costmp)
          if(ang>ang_max) ang_max=ang
          if(ang<ang_min) ang_min=ang
        enddo !j=1,4

        if(isplit==0) then !calc angle ratio for convex quad
          if(ang_max==0.or.ang_max<ang_min) stop 'failed to find max/min'
!'
          if(ang_min/ang_max<rat_angle) then
            !Decide which node (1 or 2) to split from, based on best aspect ratios
            !calc 2 A.R.'s
            do j=1,2 !2 options of splitting
              if(j==1) then
                nwild2(1:4)=(/1,2,3,4/)
              else
                nwild2(1:4)=(/2,3,4,1/)
              endif

              do k=1,2 !pair of tri's
                if(k==1) then
                  nwild(1:3)=(/1,2,3/)
                else
                  nwild(1:3)=(/1,3,4/)
                endif
                n1=elnode(nwild2(nwild(1)),i)
                n2=elnode(nwild2(nwild(2)),i)
                n3=elnode(nwild2(nwild(3)),i)
                ar3=signa(xnd(n1),xnd(n2),xnd(n3),ynd(n1),ynd(n2),ynd(n3))
                if(ar3<=0) then
                  print*, 'ar3<=0:',i,j,k,n1,n2,n3
                  stop
                endif
                rl1=sqrt((xnd(n1)-xnd(n2))**2+(ynd(n1)-ynd(n2))**2)
                rl2=sqrt((xnd(n1)-xnd(n3))**2+(ynd(n1)-ynd(n3))**2)
                rl3=sqrt((xnd(n2)-xnd(n3))**2+(ynd(n2)-ynd(n3))**2)
                rlmax=max(rl1,rl2,rl3)
                aspect_ratio(k,j)=rlmax/sqrt(ar3/pi)
              enddo !k
            enddo !j=1,2

            !Change here if you want to use different criteria
            arat1=sum(aspect_ratio(:,1))
            arat2=sum(aspect_ratio(:,2))
            if(arat1<arat2) then
              isplit=1 !split at node 1
            else
              isplit=2
            endif
          endif !ang_min/
        endif !isplit==0

        if(isplit/=0) then !add extra tri to the end of conn table
          ne_extra=ne_extra+1
          if(ne_extra>ne) stop 'overflow(1)'
          i34_new(i)=3; i34_new(ne+ne_extra)=3
          i1=isplit; i2=isplit+1; i3=isplit+2; i4=isplit+3
          if(i2>4) i2=i2-4
          if(i3>4) i3=i3-4
          if(i4>4) i4=i4-4

          !Check
          n1=elnode(i1,i); n2=elnode(i2,i); n3=elnode(i3,i); n4=elnode(i4,i)
          ar3=signa(xnd(n1),xnd(n2),xnd(n3),ynd(n1),ynd(n2),ynd(n3))
          ar4=signa(xnd(n1),xnd(n3),xnd(n4),ynd(n1),ynd(n3),ynd(n4))
          if(ar3<=0.or.ar4<=0) then
            print*, 'Final orientation wrong:',i,ar3,ar4,n1,n2,n3,n4
            stop
          endif

          nwild(1:3)=(/i1,i2,i3/)
          elnode_new(1:3,i)=elnode(nwild(1:3),i)
          nwild(1:3)=(/i1,i3,i4/)
          elnode_new(1:3,ne+ne_extra)=elnode(nwild(1:3),i)

          quad_loc(1,ne_extra)=xctr(i)
          quad_loc(2,ne_extra)=yctr(i)
        endif !isplit/=0
      enddo !i=1,ne

      print*, ne_extra,' quads split; see split_loc.bp'
      ne=ne+ne_extra

!     Output
      open(13,file='split_loc.bp',status='replace')
      write(13,*); write(13,*)ne_extra
      do i=1,ne_extra
        write(13,'(i12,2(1x,e22.14),1x,f8.1)')i,quad_loc(1:2,i),0.
      enddo !i
      close(13)

      open(12,file='hgrid.gr3.new',status='replace')
      write(12,*)'Threshold of ang ratio=',rat_angle
      write(12,*)ne,np
      do i=1,np
        write(12,'(i12,2(1x,e22.14),1x,e16.5)')i,xnd(i),ynd(i),dp(i)
      enddo !i
      do i=1,ne
        write(12,*)i,i34_new(i),elnode_new(1:i34_new(i),i)
      enddo !i
      close(12)

      stop
      end

      function signa(x1,x2,x3,y1,y2,y3)
!...  Compute signed area formed by pts 1,2,3
      implicit real*8(a-h,o-z)

      signa=((x1-x3)*(y2-y3)-(x2-x3)*(y1-y3))/2

      return
      end

max_stars_repo_name: rustychris/schism

max_issues_repo_name: rustychris/schism

max_forks_repo_name: rustychris/schism


text:
	SUBROUTINE SHN_BUFR  ( iret )
C************************************************************************
C* SHN_BUFR								*
C*									*
C* This subroutine retrieves data from the interface arrays, converts	*
C* it into BUFR output, and then writes the BUFR output to the BUFR	*
C* output stream.							*
C*									*
C* SHN_BUFR  ( IRET )							*
C*									*
C* Input parameters:							*
C*									*
C* Output parameters:							*
C*	IRET		INTEGER		Return code:			*
C*					  0 = normal return		*
C**									*
C* Log:									*
C* J. Ator/NCEP		04/05						*
C* J. Ator/NCEP		07/06	Set BUFR types/subtypes based upon STYP	*
C************************************************************************
	INCLUDE		'GEMPRM.PRM'
	INCLUDE		'BRIDGE.PRM'
	INCLUDE		'shncmn.cmn'
C*
	PARAMETER	( MAXTD = 10 )
C*
	CHARACTER	ctds(MAXTD)*2, ctd*2, cimn*8, bfrstr*4, iqual,
     +			bfstyp*8, timctd*2
C*
	INTEGER		itds(MAXTD)
C*
	LOGICAL		newtd
C*
        REAL*8          r8ary (4), UT_RIBM, PKFTBV
C*
	INCLUDE		'ERMISS.FNC'
C*-----------------------------------------------------------------------
	iret = 0
C
C*	Divide the interface values into groups based upon their
C*	associated time durations.  Then, encode a separate BUFR
C*	report for each such group of values.
C
	ntd = 0
	DO WHILE ( .true. )
C
C*	    Determine the next time duration for which to encode
C*	    a BUFR report.
C
	    irptr = 9 + ( ntd * 3 )
	    newtd = .false.
	    DO WHILE  ( ( .not. newtd ) .and.
     +		       ( ( irptr + 2 ) .le. nimn ) )
		ctd = rimnem ( irptr )(5:6)
		itd = INT ( rivals ( irptr ) )
                timdur = rivals (irptr)
                timctd = ctd
C
C*		Is this a "new" time duration (for which we have not
C*		yet encoded a BUFR report)?
C
		i = 1
		newtd = .true.
		DO WHILE  ( ( i .le. ntd ) .and. ( newtd ) )
		    IF  ( ( ctd .eq. ctds (i) ) .and.
     +			  ( itd .eq. itds (i) ) )  THEN
			newtd = .false.
		        irptr = irptr + 3
		    ELSE
			i = i + 1
		    END IF
		END DO
	    END DO
	    IF  ( .not. newtd )  THEN
C
C*		We are done!
C
		RETURN
	    END IF
C
C*	    Initialize a BUFR report for this time duration.
C
	    year = rivals (1)
	    rmth = rivals (2)
	    days = rivals (3)
	    hour = rivals (4)
	    rmin = rivals (5)
	    IF  ( ( ERMISS ( year ) ) .or. ( ERMISS ( rmth ) ) .or.
     +		  ( ERMISS ( days ) ) .or. ( ERMISS ( hour ) )  )  THEN
		RETURN
	    END IF
C
	    ibfdt = ( INT ( year ) * 1000000 ) +
     +		    ( INT ( rmth ) * 10000 )  +
     +		    ( INT ( days ) * 100 ) + INT ( hour )
C
	    IF  ( ( civals (2)(1:4) .eq. 'COOP' ) .or.
     +		  ( civals (2)(1:5) .eq. 'UCOOP' ) )  THEN
		bfstyp = 'NC255102'
	    ELSE IF ( civals (2)(1:5) .eq. 'C-MAN' ) THEN
		bfstyp = 'NC001010'
	    ELSE IF ( civals (2)(1:4) .eq. 'BUOY' ) THEN
		bfstyp = 'NC001011'
	    ELSE IF ( ( civals (2)(1:4) .eq. 'CTGN' ) .or.
     +		      ( civals (2)(1:3) .eq. 'HTG' ) )  THEN
		bfstyp = 'NC001012'
	    ELSE
		bfstyp = 'NC000010'
	    END IF
C
	    CALL OPENMB  ( iunbfo, bfstyp, ibfdt )
C
C*	    Report ID.
C
            IF ( (timdur .eq. 24.0) .and. (timctd .eq. 'HR') ) then
	    CALL UT_CIBF  ( iunbfo, 'RPID', civals (1), 8, iercbf )
C
C*	    Report date-time.
C
	    CALL UT_RIBF  ( iunbfo, 'YEAR', year, ierrbf )
	    CALL UT_RIBF  ( iunbfo, 'MNTH', rmth, ierrbf )
	    CALL UT_RIBF  ( iunbfo, 'DAYS', days, ierrbf )
	    CALL UT_RIBF  ( iunbfo, 'HOUR', hour, ierrbf )
	    CALL UT_RIBF  ( iunbfo, 'MINU', rmin, ierrbf )
C
C*	    Time duration.
C
	    IF  ( .not. ERMISS ( rivals ( irptr ) ) )  THEN
		CALL UT_RIBF  ( iunbfo, 'TP' // ctd ,
     +				( rivals ( irptr ) * (-1) ) , ierrbf )
	    END IF
C
C*	    Latitude, longitude and elevation.
C
	    CALL UT_RIBF  ( iunbfo, 'CLAT' , rivals (6), ierrbf )
	    CALL UT_RIBF  ( iunbfo, 'CLON' , rivals (7), ierrbf )
	    CALL UT_RIBF  ( iunbfo, 'HSMSL', rivals (8), ierrbf )
C
C*	    Receipt date-time.
C
	    CALL UT_RIBF ( iunbfo, 'RCYR', FLOAT ( irundt(1) ), ierrbf )
	    CALL UT_RIBF ( iunbfo, 'RCMO', FLOAT ( irundt(2) ), ierrbf )
	    CALL UT_RIBF ( iunbfo, 'RCDY', FLOAT ( irundt(3) ), ierrbf )
	    CALL UT_RIBF ( iunbfo, 'RCHR', FLOAT ( irundt(4) ), ierrbf )
	    CALL UT_RIBF ( iunbfo, 'RCMI', FLOAT ( irundt(5) ), ierrbf )
	    CALL UT_RIBF ( iunbfo, 'RCTS', 0., ierrbf )
C
C*	    Bulletin ID information.
C
	    CALL UT_CIBF  ( iunbfo, 'SEQNUM', seqnum, 8, iercbf )
	    CALL UT_CIBF  ( iunbfo, 'BUHD', buhd, 8, iercbf )
	    CALL UT_CIBF  ( iunbfo, 'BORG', cborg, 8, iercbf )
	    CALL UT_CIBF  ( iunbfo, 'BULTIM', bulldt, 8, iercbf )
	    CALL UT_CIBF  ( iunbfo, 'BBB', bbb, 8, iercbf )
C
C*	    Locate all other interface values which have this same time
C*	    duration and encode them into this BUFR report.
C
	    nst = 0
	    nmxst = 0
	    nmist = 0
	    DO WHILE  ( ( irptr + 2 ) .le. nimn )
		IF  ( ( rimnem ( irptr )(5:6) .eq. ctd ) .and.
     +		      ( INT ( rivals ( irptr ) ) .eq. itd ) )  THEN
C
		    rvalu = rivals ( irptr + 1 )
C
C*		    Encode this interface value, along with its
C*		    associated SHEF revision and qualifier flags,
C*		    into the BUFR report.
C*
		    cimn  = rimnem ( irptr + 1 )
		    IF  ( ( cimn(1:2) .eq. 'MX' ) .or.
     +			  ( cimn(1:2) .eq. 'MN' ) )  THEN
			ipt = 3
		    ELSE
			ipt = 1
		    END IF
C
C*		    SHEF revision flag.
C
		    r8ary (3) = UT_RIBM ( rivals ( irptr + 2 ) )
C
C*		    SHEF qualifier flag.
C
		    iqual = civals ( irptr + 1 )(1:1)
		    r8ary (4) =
     +			UT_RIBM ( FLOAT ( MOVA2I ( iqual ) - 64 ) )
C
		    IF  ( cimn(1:4) .eq. 'TPCI' )  THEN
C
C*			Total precipitation.
C
			r8ary (2) = UT_RIBM ( PR_INMM ( rvalu ) )
			CALL SHN_BFSQ  ( r8ary, 3, 'TOPC', ierfsq )
C
		    ELSE IF  ( cimn .eq. 'SHPT    ' )  THEN
C
C*			Type of precipitation.
C
			ishpt = INT ( rvalu )
			IF  ( ( ishpt .ge. 0 ) .and.
     +			      ( ishpt .le. 9 ) ) THEN
			  IF  ( ishpt .eq. 0 )  THEN
			    r8ary (2) = PKFTBV ( 30, 19 )
			  ELSE IF  ( ishpt .eq. 1 )  THEN
			    r8ary (2) = PKFTBV ( 30,  5 )
			  ELSE IF  ( ishpt .eq. 2 )  THEN
			    r8ary (2) = PKFTBV ( 30,  3 ) +
     +					PKFTBV ( 30,  5 )
			  ELSE IF  ( ishpt .eq. 3 )  THEN
			    r8ary (2) = PKFTBV ( 30,  4 )
			  ELSE IF  ( ishpt .eq. 4 )  THEN
			    r8ary (2) = PKFTBV ( 30,  3 ) +
     +					PKFTBV ( 30,  4 )
			  ELSE IF  ( ishpt .eq. 5 )  THEN
			    r8ary (2) = PKFTBV ( 30,  7 )
			  ELSE IF  ( ishpt .eq. 6 )  THEN
			    r8ary (2) = PKFTBV ( 30,  9 )
			  ELSE IF  ( ishpt .eq. 7 )  THEN
			    r8ary (2) = PKFTBV ( 30,  8 )
			  ELSE IF  ( ishpt .eq. 8 )  THEN
			    r8ary (2) = PKFTBV ( 30, 10 )
			  ELSE
			    r8ary (2) = PKFTBV ( 30, 14 )
			  END IF
			  CALL SHN_BFSQ  ( r8ary, 3, 'PRTP', ierfsq )
			END IF
C
		    ELSE IF  ( cimn .eq. 'RSHF    ' )  THEN
C
C*			River stage height.
C
			r8ary (2) = UT_RIBM ( PR_HGFM ( rvalu ) )
			CALL SHN_BFSQ  ( r8ary, 3, 'RSHM', ierfsq )
C
		    ELSE IF  ( ( cimn(1:2) .eq. 'SN' ) .and.
     +			       ( cimn(4:4) .eq. 'W' ) )  THEN
C
C*			Snow depth (fresh or total).
C
			r8ary (2) =
     +			    UT_RIBM ( PR_HGMK ( PR_INMM ( rvalu ) ) )
			IF  ( cimn(3:3) .eq. 'E' )  THEN
			  bfrstr = 'DOFS'
			ELSE
			  bfrstr = 'TOSD'
			END IF
			CALL SHN_BFSQ  ( r8ary, 3, bfrstr, ierfsq )
C
		    ELSE IF  ( cimn .eq. 'WEQS    ' )  THEN
C
C*			Snow water equivalent.
C
			r8ary (2) = UT_RIBM ( PR_INMM ( rvalu ) )
			CALL SHN_BFSQ  ( r8ary, 3, 'SWEM', ierfsq )
C
		    ELSE IF  ( cimn(ipt:ipt+3) .eq. 'TMPF' )  THEN
C
C*			Temperature.
C
			r8ary (2) = UT_RIBM ( PR_TMFK ( rvalu ) )
			IF  ( cimn(2:2) .eq. 'X' )  THEN
			  bfrstr = 'MXTM'
			ELSE IF  ( cimn(2:2) .eq. 'N' )  THEN
			  bfrstr = 'MITM'
			ELSE
			  bfrstr = 'TMDB'
			END IF
			CALL SHN_BFSQ  ( r8ary, 3, bfrstr, ierfsq )
C
		    ELSE IF  ( cimn .eq. 'TMWF    ' )  THEN
C
C*			Wet bulb temperature.
C
			r8ary (2) = UT_RIBM ( PR_TMFK ( rvalu ) )
			CALL SHN_BFSQ  ( r8ary, 3, 'TMWB', ierfsq )
C
		    ELSE IF  ( cimn .eq. 'DWPF    ' )  THEN
C
C*			Dew point temperature.
C
			r8ary (2) = UT_RIBM ( PR_TMFK ( rvalu ) )
			CALL SHN_BFSQ  ( r8ary, 3, 'TMDP', ierfsq )
C
		    ELSE IF  ( cimn .eq. 'WTMF    ' )  THEN
C
C*			Water temperature.
C
			r8ary (2) = UT_RIBM ( PR_TMFK ( rvalu ) )
			CALL SHN_BFSQ  ( r8ary, 3, 'WATM', ierfsq )
C
		    ELSE IF  ( ( cimn(ipt:ipt+3) .eq. 'CSTD' ) .or.
     +			       ( cimn(ipt:ipt+3) .eq. 'STMF' ) )  THEN
C
C*			Soil temperature (with or without depth).
C
			IF  ( cimn (ipt:ipt) .eq. 'C' )  THEN
			  icstdn = NINT ( rvalu * 1000. )
			  dblsi = FLOAT ( icstdn / 1000 )
			  stmf = ( FLOAT ( MOD ( icstdn, 1000 ) ) )
			ELSE
			  dblsi = RMISSD
			  stmf = rvalu
			END IF
C
			r8ary (1) =
     +			    UT_RIBM ( PR_HGMK ( PR_INMM ( dblsi ) ) )
			r8ary (2) = UT_RIBM ( PR_TMFK ( stmf ) )
C
			IF  ( cimn(2:2) .eq. 'X' )  THEN
			  nmxst = nmxst + 1
			  IF  ( nmxst .eq. 1 )  THEN
			    bfrstr = 'XST1'
			  ELSE
			    bfrstr = 'XST2'
			  END IF
			ELSE IF  ( cimn(2:2) .eq. 'N' )  THEN
			  nmist = nmist + 1
			  IF  ( nmist .eq. 1 )  THEN
			    bfrstr = 'IST1'
			  ELSE
			    bfrstr = 'IST2'
			  END IF
			ELSE
			  nst = nst + 1
			  IF  ( nst .eq. 1 )  THEN
			    bfrstr = 'STM1'
			  ELSE
			    bfrstr = 'STM2'
			  END IF
			END IF
C
			CALL SHN_BFSQ  ( r8ary, 4, bfrstr, ierfsq )
C
		    ELSE IF  ( cimn(ipt:ipt+3) .eq. 'RELH' )  THEN
C
C*			Relative humidity.
C
			r8ary (2) = UT_RIBM ( rvalu )
			IF  ( cimn(2:2) .eq. 'X' )  THEN
			  bfrstr = 'MXRH'
			ELSE IF  ( cimn(2:2) .eq. 'N' )  THEN
			  bfrstr = 'MIRH'
			ELSE
			  bfrstr = 'REHU'
			END IF
			CALL SHN_BFSQ  ( r8ary, 3, bfrstr, ierfsq )
C
		    ELSE IF  ( ( cimn .eq. 'DRCT    ' ) .or.
     +			       ( cimn .eq. 'PWDR    ' ) )  THEN
C
C*			Wind direction (normal or peak).
C
			r8ary (2) = UT_RIBM ( rvalu )
			IF  ( cimn(1:1) .eq. 'D' )  THEN
			  bfrstr = 'WDIR'
			ELSE
			  bfrstr = 'PKWD'
			END IF
			CALL SHN_BFSQ  ( r8ary, 3, bfrstr, ierfsq )
C
		    ELSE IF  ( cimn(2:4) .eq. 'MPH' )  THEN
C
C*			Wind speed (normal, peak or gusts).
C
			r8ary (2) =
     +			    UT_RIBM ( PR_KNMS ( PR_MHKN ( rvalu ) ) )
			IF  ( cimn(1:1) .eq. 'S' )  THEN
			  bfrstr = 'WSPD'
			ELSE IF  ( cimn(1:1) .eq. 'P' )  THEN
			  bfrstr = 'PKWS'
			ELSE
			  bfrstr = 'MXGS'
			END IF
			CALL SHN_BFSQ  ( r8ary, 3, bfrstr, ierfsq )
C
		    ELSE IF  ( cimn .eq. 'CWDS    ' )  THEN
C
C*			Combined wind direction and speed.
C
			icwdsn = NINT ( rvalu * 10000. )
C
			r8ary (2) =
     +			    UT_RIBM ( FLOAT ( MOD ( icwdsn, 1000 ) ) )
			CALL SHN_BFSQ  ( r8ary, 3, 'WDIR', ierfsq )
C
			smph = FLOAT ( icwdsn / 1000 ) / 10.
			r8ary (2) =
     +			    UT_RIBM ( PR_KNMS ( PR_MHKN ( smph ) ) )
			CALL SHN_BFSQ  ( r8ary, 3, 'WSPD', ierfsq )
C
		    ELSE IF  ( ( cimn(1:1) .eq. 'P' ) .and.
     +			       ( cimn(3:4) .eq. 'SI' ) )  THEN
C
C*			Pressure (normal or reduced to MSL).
C
			IF  ( rvalu .lt. 900. )  THEN
C
C*			  The value is in inches, so convert to mb.
C
			  prmb = PR_ALTM ( rvalu )
			ELSE
			  prmb = rvalu
			END IF
C
			r8ary (2) = UT_RIBM ( PR_M100 ( prmb ) )
			IF  ( cimn(2:2) .eq. 'R' )  THEN
			  bfrstr = 'PRES'
			ELSE
			  bfrstr = 'PMSL'
			END IF
			CALL SHN_BFSQ  ( r8ary, 3, bfrstr, ierfsq )
C
		    ELSE IF  ( INDEX ( cimn, 'WWM' ) .ne. 0 )  THEN
C
C*			Past or present weather.
C
			r8ary (2) = UT_RIBM ( rvalu )
			IF  ( cimn(1:1) .eq. 'P' )  THEN
			  bfrstr = 'PSW1'
			ELSE
			  bfrstr = 'PRWE'
			END IF
			CALL SHN_BFSQ  ( r8ary, 3, bfrstr, ierfsq )
C
		    ELSE IF  ( cimn .eq. 'VSBY    ' )  THEN
C
C*			Horizontal visibility.
C
			r8ary (2) =
     +			    UT_RIBM ( PR_HGFM ( PR_HGSF ( rvalu ) ) )
			CALL SHN_BFSQ  ( r8ary, 3, 'HOVI', ierfsq )
C
		    ELSE IF  ( ( cimn .eq. 'SOGR    ' ) .or.
C
C*				State of the ground.
C
     +			       ( cimn .eq. 'CLAM    ' ) .or.
C
C*				Cloud amount.
C
     +			       ( cimn .eq. 'ALBD    ' ) .or.
C
C*				Albedo.
C
     +			       ( cimn .eq. 'TOSH    ' ) )  THEN
C
C*				Total sunshine.
C
			r8ary (2) = UT_RIBM ( rvalu )
			CALL SHN_BFSQ  ( r8ary, 3, cimn(1:4), ierfsq )
C
		    ELSE
C
			logmsg = '  unable to create BUFR output for '
     +				// cimn
			CALL DC_WLOG ( 2, 'DC', 2, logmsg, ierwlg)
C
		    END IF
C
		END IF
C
		irptr = irptr + 3
C
	    END DO
C
C*	    Store the BUFR report for this time duration.
C
	    CALL UT_WBFR  ( iunbfo, 'shef', .false., ierwbf )

C
C* TIME DURATION END IF
C
            END IF
C
C*	    Update the ctds and itds arrays, then go back and determine
C*	    the next time duration for which to encode a BUFR report.
C
	    ntd = ntd + 1
	    IF  ( ntd .gt. MAXTD )  THEN
		logmsg = '  MAXTD overflow when generating BUFR output.'
		CALL DC_WLOG ( 2, 'DC', 2, logmsg, ierwlg)
		RETURN
	    END IF
	    ctds ( ntd ) = ctd
	    itds ( ntd ) = itd
	END DO
C*
	RETURN
	END


text:
only forth definitions

\ test
create s 100 *stack drop
: numbers  locals| c |  c vacate  c capacity 0 do  i c push  loop ;
s numbers

:noname  %node venery:sizeof allocate throw dup /node ; is new-node
:noname  free throw ; is free-node

new-node constant p
new-node constant n1  n1 p push
new-node constant n2  n2 p push
new-node constant n3  n3 p push
new-node constant n4  n4 p push
new-node constant p2
new-node constant n5  n5 p2 push
new-node constant n6  n6 p2 push
new-node constant n7  n7 p2 push
new-node constant n8  n8 p2 push

only forth definitions

max_stars_repo_name: RogerLevy/RamenEngine

max_issues_repo_name: Alim-Oezdemir/RamenEngine

max_forks_repo_name: Alim-Oezdemir/RamenEngine