! ##################################################################################################################################
! Begin MIT license text.
! _______________________________________________________________________________________________________

! Copyright 2019 Dr William R Case, Jr (dbcase29@gmail.com)

! 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 and documentation.

! 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.
! _______________________________________________________________________________________________________

! End MIT license text.

      SUBROUTINE BD_RBE1 ( CARD, LARGE_FLD_INP )

! Processes RBE1 Bulk Data Cards. Writes RBE1 element records to file L1F for later processing.
! Two records are written for each dependent Grid/DOF pair:

!       1) Record 1 has the rigid element type: 'RBE1    '
!       2) Record 2 has:
!             RELID           : Rigid element ID
!             IGID(i),IDOF(i) : Pairs of indep Grid/DOF (up to 6)
!             DGID, DDOF      : One pair of dependent Grid/DOF

      USE PENTIUM_II_KIND, ONLY       :  BYTE, LONG, DOUBLE
      USE IOUNT1, ONLY                :  WRT_ERR, WRT_LOG, ERR, F04, F06, L1F
      USE SCONTR, ONLY                :  BLNK_SUB_NAM, FATAL_ERR, IERRFL, JCARD_LEN, JF, LRIGEL, NRBE1, NRIGEL, NRECARD
      USE TIMDAT, ONLY                :  TSEC
      USE SUBR_BEGEND_LEVELS, ONLY    :  BD_RBE1_BEGEND
      USE MODEL_STUF, ONLY            :  RIGID_ELEM_IDS

      USE BD_RBE1_USE_IFs

      IMPLICIT NONE

      CHARACTER(LEN=LEN(BLNK_SUB_NAM)):: SUBR_NAME = 'BD_RBE1'
      CHARACTER(LEN=*), INTENT(INOUT) :: CARD              ! A Bulk Data card
      CHARACTER(LEN=*), INTENT(IN)    :: LARGE_FLD_INP     ! If 'Y', CARD is large field format
      CHARACTER(LEN(CARD))            :: CHILD             ! "Child" card read in subr NEXTC, called herein
      CHARACTER(LEN=JCARD_LEN)        :: JCARD(10)         ! The 10 fields of characters making up CARD
      CHARACTER(LEN(JCARD))           :: CHR_IDOF(6)
      CHARACTER( 8*BYTE)              :: IP6TYP            ! An output from subr IP6CHK called herein
      CHARACTER(LEN(JCARD))           :: JCARD1_PARENT     ! Field 1 of parent card
      CHARACTER(LEN(JCARD))           :: JCARD2_PARENT     ! Field 2 of parent card
      CHARACTER(LEN(JCARD))           :: JCARDO            ! An output from subr IP6CHK called herein
      CHARACTER( 1*BYTE)              :: MORE_IDOF = 'N'   ! = 'Y' if a cont card with indep DOF's is found
      CHARACTER( 1*BYTE)              :: UMFND             ! = 'Y' when we find "UM" in field 2 of cont card no. 1 or 2
      CHARACTER( 8*BYTE), PARAMETER   :: RTYPE = 'RBE1    '! Rigid element type

      INTEGER(LONG)                   :: CONT_NO   = 0     ! Count of the continuation cards for this parent
      INTEGER(LONG)                   :: DGID      = 0     ! A dependent grid
      INTEGER(LONG)                   :: DDOF      = 0     ! Dependent DOF's at DGID
      INTEGER(LONG)                   :: DOF_NOS(6)= 0     ! Fields 5 & 6 should fill this in to be: 1 2 3 4 5 6
      INTEGER(LONG)                   :: IDOF_ERR  = 0     ! Count of the no. of DOF component errors in fields 5,6
      INTEGER(LONG)                   :: I,J               ! DO loop indices
      INTEGER(LONG)                   :: ICONT     = 0     ! Indicator of whether a cont card exists. Output from subr NEXTC
      INTEGER(LONG)                   :: IDOF(6)   = 0     ! Independent DOF's at IGID(i)
      INTEGER(LONG)                   :: IDUM              ! Dummy arg in subr IP^CHK not used herein
      INTEGER(LONG)                   :: IERR      = 0     ! Error indicator returned from subr NEXTC called herein
      INTEGER(LONG)                   :: IGID(6)   = 0     ! Independent grid ID's
      INTEGER(LONG)                   :: INT1      = 0     ! An integer 1-6 read from internal file
      INTEGER(LONG)                   :: JERR      = 0     ! A local error count
      INTEGER(LONG)                   :: JFLD1             ! A computed field number on the card
      INTEGER(LONG)                   :: JFLD2             ! A computed field number on the card
      INTEGER(LONG)                   :: NUM_IDOF_FLDS = 0 ! Number of fields that have independent DOF's specified
      INTEGER(LONG)                   :: RELID     = 0     ! This rigid elements' ID
      INTEGER(LONG), PARAMETER        :: SUBR_BEGEND = BD_RBE1_BEGEND

! **********************************************************************************************************************************
      IF (WRT_LOG >= SUBR_BEGEND) THEN
         CALL OURTIM
         WRITE(F04,9001) SUBR_NAME,TSEC
 9001    FORMAT(1X,A,' BEGN ',F10.3)
      ENDIF

! **********************************************************************************************************************************
! RBE1 Bulk Data Card routine

!   FIELD   ITEM
!   -----   ------------
!    2      RELID  , Rigid Elem ID
!    3      IGID(1), Grid ID for 1st independent grid
!    4      IDOF(1), DOF's at 1st independent grid
!    5      IGID(2), Grid ID for 2nd independent grid, if it exists
!    6      IDOF(2), DOF's at 2nd independent grid, if IGID(2) exists
!    7      IGID(3), Grid ID for 3rd independent grid, if it exists
!    8      IDOF(3), DOF's at 3rd independent grid, if IGID(3) exists

! Possible continuation card (if there are 4, 5 or 6 independent grids)
!   FIELD   ITEM
!   -----   ------------
!    3      IGID(4), Grid ID for 4th independent grid, if it exists
!    4      IDOF(4), DOF's at 4th independent grid, if IGID(4) exists
!    5      IGID(5), Grid ID for 5th independent grid, if it exists
!    6      IDOF(5), DOF's at 5th independent grid, if IGID(5) exists
!    7      IGID(6), Grid ID for 6th independent grid, if it exists
!    8      IDOF(6), DOF's at 6th independent grid, if IGID(6) exists
!    The collection of IDOF(i) must yield 6 DOF's that completely  describe a general rigid body motion of the
!    rigid element

! Mandatory continuation card
!   FIELD   ITEM
!   -----   ------------
!    2      "UM"
!    3      DGID, Grid ID for 1st dependent grid
!    4      DDOF, DOF's at 1st dependent grid
!    5-8    Up to 2 more pairs of DGID/DDOF, if they exist, followed by more continuation cards with up to 3 pairs
!           of DGID/DDOF in fields 3-8, if needed

! Data is written to file L1F for later processing after checks on
! format of data.

! Make JCARD from CARD

      CALL MKJCARD ( SUBR_NAME, CARD, JCARD )
      JCARD1_PARENT = JCARD(1)
      JCARD2_PARENT = JCARD(2)

! Check for overflow

      NRBE1  = NRBE1+1
      NRIGEL = NRIGEL+1
!xx   IF (NRIGEL > LRIGEL) THEN
!xx      FATAL_ERR = FATAL_ERR + 1
!xx      WRITE(ERR,1163) SUBR_NAME,JCARD(1),LRIGEL
!xx      WRITE(F06,1163) SUBR_NAME,JCARD(1),LRIGEL
!xx      CALL OUTA_HERE ( 'Y' )                            ! Coding error, so quit
!xx   ENDIF

! Initialize CHR_IDOF

      DO I=1,6
         CHR_IDOF(I)(1:) = ' '
      ENDDO

! Read Elem ID

      CALL I4FLD ( JCARD(2), JF(2), RELID )
      IF (IERRFL(2) /= 'N') THEN
         JERR = JERR + 1
      ELSE
         RIGID_ELEM_IDS(NRIGEL) = RELID
      ENDIF

! Read up to 3 pairs of independent grids and DOF's in fields 3-8 of the parent card.
! For any of the 3 pairs, if either of the 2 fields is blank, the other one must also be blank.

      NUM_IDOF_FLDS = 0
      DO J=1,3
         JFLD1 = 2*J+1
         JFLD2 = 2*J+2
         IF      ((JCARD(JFLD1)(1:) == ' ') .AND. (JCARD(JFLD2)(1:) == ' ')) THEN
            CYCLE
         ELSE IF ((JCARD(JFLD1)(1:) /= ' ') .AND. (JCARD(JFLD2)(1:) == ' ')) THEN
            FATAL_ERR = FATAL_ERR + 1
            WRITE(ERR,1181) JCARD2_PARENT,JFLD1,JFLD2
            WRITE(F06,1181) JCARD2_PARENT,JFLD1,JFLD2
         ELSE IF ((JCARD(JFLD1)(1:) == ' ') .AND. (JCARD(JFLD2)(1:) /= ' ')) THEN
            FATAL_ERR = FATAL_ERR + 1
            WRITE(ERR,1182) JCARD2_PARENT,JFLD1,JFLD2
            WRITE(F06,1182) JCARD2_PARENT,JFLD1,JFLD2
         ELSE IF ((JCARD(JFLD1)(1:) /= ' ') .AND. (JCARD(JFLD2)(1:) /= ' ')) THEN
            NUM_IDOF_FLDS = NUM_IDOF_FLDS + 1
            CHR_IDOF(NUM_IDOF_FLDS) = JCARD(JFLD2)
            CALL I4FLD ( JCARD(JFLD1), JF(JFLD1), IGID(NUM_IDOF_FLDS) )
            IF (IERRFL(JFLD1) /= 'N') THEN
               JERR = JERR + 1
            ENDIF
            CALL I4FLD ( JCARD(JFLD2), JF(JFLD2), IDOF(NUM_IDOF_FLDS) )
            IF (IERRFL(JFLD2) == 'N') THEN
               CALL IP6CHK ( JCARD(JFLD2), JCARDO, IP6TYP, IDUM )
               IF ((IP6TYP == 'COMP NOS') .OR. (IP6TYP == 'BLANK   '))  THEN
                  CONTINUE
               ELSE
                  FATAL_ERR  = FATAL_ERR + 1
                  IDOF_ERR = IDOF_ERR + 1
                  WRITE(ERR,1124) JFLD2,JCARD1_PARENT,JCARD2_PARENT,JFLD2,JCARD(JFLD2)
                  WRITE(F06,1124) JFLD2,JCARD1_PARENT,JCARD2_PARENT,JFLD2,JCARD(JFLD2)
               ENDIF
            ELSE
               JERR = JERR + 1
            ENDIF
         ENDIF
      ENDDO

      CALL BD_IMBEDDED_BLANK ( JCARD,2,3,0,5,0,7,0,0 )     ! Make sure no imbedded blanks in fields 2,3,5,7. Flds 4,6,8 ae DOF's
      CALL CARD_FLDS_NOT_BLANK ( JCARD,0,0,0,0,0,0,0,9 )     ! Issue warning if field 9 not blank
      CALL CRDERR ( CARD )                                 ! CRDERR prints errors found when reading fields

! Now look for a continuation card. It may have up to 3 additional pairs of independent grid/DOF or it may begin the
! dependent grid/DOF list. If it is the later, then field 2 must have "UM"

      CONT_NO = 0
      UMFND   = 'N'
      DO
         IF (LARGE_FLD_INP == 'N') THEN
            CALL NEXTC  ( CARD, ICONT, IERR )
         ELSE
            CALL NEXTC2 ( CARD, ICONT, IERR, CHILD )
            CARD = CHILD
         ENDIF
         CALL MKJCARD ( SUBR_NAME, CARD, JCARD )
         IF (ICONT == 1) THEN
            CONT_NO = CONT_NO + 1
            IF (JCARD(2)(1:2) == 'UM') THEN
               UMFND = 'Y'
            ENDIF
            IF (UMFND =='N') THEN                          ! Read up to 3 additional pairs of indep grid/DOF's
               IF (CONT_NO == 1) THEN
                  MORE_IDOF = 'Y'
                  DO J=1,3
                     JFLD1 = 2*J+1
                     JFLD2 = 2*J+2
                     IF      ((JCARD(JFLD1)(1:) == ' ') .AND. (JCARD(JFLD2)(1:) == ' ')) THEN
                        CYCLE
                     ELSE IF ((JCARD(JFLD1)(1:) /= ' ') .AND. (JCARD(JFLD2)(1:) == ' ')) THEN
                        FATAL_ERR = FATAL_ERR + 1
                        WRITE(ERR,1181) JCARD2_PARENT,JFLD1,JFLD2
                        WRITE(F06,1181) JCARD2_PARENT,JFLD1,JFLD2
                     ELSE IF ((JCARD(JFLD1)(1:) == ' ') .AND. (JCARD(JFLD2)(1:) /= ' ')) THEN
                        FATAL_ERR = FATAL_ERR + 1
                        WRITE(ERR,1182) JCARD2_PARENT,JFLD1,JFLD2
                        WRITE(F06,1182) JCARD2_PARENT,JFLD1,JFLD2
                     ELSE IF ((JCARD(JFLD1)(1:) /= ' ') .AND. (JCARD(JFLD2)(1:) /= ' ')) THEN
                        NUM_IDOF_FLDS = NUM_IDOF_FLDS + 1
                        CHR_IDOF(NUM_IDOF_FLDS) = JCARD(JFLD2)
                        CALL I4FLD ( JCARD(JFLD1), JF(JFLD1), IGID(NUM_IDOF_FLDS) )
                        IF (IERRFL(JFLD1) /= 'N') THEN
                           JERR = JERR + 1
                        ENDIF
                        CALL I4FLD ( JCARD(JFLD2), JF(JFLD2), IDOF(NUM_IDOF_FLDS) )
                        IF (IERRFL(JFLD2) == 'N') THEN
                           CALL IP6CHK ( JCARD(JFLD2), JCARDO, IP6TYP, IDUM )
                           IF ((IP6TYP == 'COMP NOS') .OR. (IP6TYP == 'BLANK   '))  THEN
                              CONTINUE
                           ELSE
                              FATAL_ERR  = FATAL_ERR + 1
                              IDOF_ERR = IDOF_ERR + 1
                              WRITE(ERR,1124) JFLD2,JCARD1_PARENT,JCARD2_PARENT,JFLD2,JCARD(JFLD2)
                              WRITE(F06,1124) JFLD2,JCARD1_PARENT,JCARD2_PARENT,JFLD2,JCARD(JFLD2)
                           ENDIF
                        ELSE
                           JERR = JERR + 1
                        ENDIF
                     ENDIF
                  ENDDO
               ELSE
                  FATAL_ERR = FATAL_ERR + 1
                  WRITE(ERR,1135) JCARD1_PARENT, JCARD2_PARENT, CONT_NO
                  WRITE(F06,1135) JCARD1_PARENT, JCARD2_PARENT, CONT_NO
                  RETURN                                   ! Cannot continue, NUM_IDOF_FLDS cannot exceed 6 (for array CHR_IDOF)
               ENDIF

               CALL BD_IMBEDDED_BLANK ( JCARD,0,3,0,5,0,7,0,0 ) ! Make sure no imbedded blanks in fields 3,5,7 (4,6,8 can be IDOF's)
               CALL CARD_FLDS_NOT_BLANK ( JCARD,2,0,0,0,0,0,0,9 )
               CALL CRDERR ( CARD )                        ! CRDERR prints errors found when reading fields

            ELSE                                           ! Found "UM", so read up to 3 pairs of dep. grid/DOF on this con't card
               DO                                          ! There are any number of continuation cards w/dep. grid/DOF
                  DO J=1,3
                     JFLD1 = 2*J+1
                     JFLD2 = 2*J+2
                     IF      ((JCARD(JFLD1)(1:) == ' ') .AND. (JCARD(JFLD2)(1:) == ' ')) THEN
                        CYCLE
                     ELSE IF ((JCARD(JFLD1)(1:) /= ' ') .AND. (JCARD(JFLD2)(1:) == ' ')) THEN
                        FATAL_ERR = FATAL_ERR + 1
                        WRITE(ERR,1181) JCARD2_PARENT,JFLD1,JFLD2
                        WRITE(F06,1181) JCARD2_PARENT,JFLD1,JFLD2
                     ELSE IF ((JCARD(JFLD1)(1:) == ' ') .AND. (JCARD(JFLD2)(1:) /= ' ')) THEN
                        FATAL_ERR = FATAL_ERR + 1
                        WRITE(ERR,1182) JCARD2_PARENT,JFLD1,JFLD2
                        WRITE(F06,1182) JCARD2_PARENT,JFLD1,JFLD2
                     ELSE IF ((JCARD(JFLD1)(1:) /= ' ') .AND. (JCARD(JFLD2)(1:) /= ' ')) THEN
                        CALL I4FLD ( JCARD(JFLD1), JF(JFLD1), DGID )
                        IF (IERRFL(JFLD1) /= 'N') THEN
                           JERR = JERR + 1
                        ENDIF
                        CALL I4FLD ( JCARD(JFLD2), JF(JFLD2), DDOF )
                        IF (IERRFL(JFLD2) == 'N') THEN
                           CALL IP6CHK ( JCARD(JFLD2), JCARDO, IP6TYP, IDUM )
                           IF ((IP6TYP == 'COMP NOS') .OR. (IP6TYP == 'BLANK   '))  THEN
                              CONTINUE
                           ELSE
                              FATAL_ERR  = FATAL_ERR + 1
                              IDOF_ERR = IDOF_ERR + 1
                              WRITE(ERR,1124) JFLD2,JCARD1_PARENT,JCARD2_PARENT,JFLD2,JCARD(JFLD2)
                              WRITE(F06,1124) JFLD2,JCARD1_PARENT,JCARD2_PARENT,JFLD2,JCARD(JFLD2)
                           ENDIF
                        ELSE
                           JERR = JERR + 1
                        ENDIF
                        IF (JERR == 0) THEN
                           WRITE(L1F) RTYPE
                           WRITE(L1F) RELID,DGID,DDOF,NUM_IDOF_FLDS,(IGID(I),IDOF(I),I=1,NUM_IDOF_FLDS)
                           NRECARD = NRECARD + 1
                        ENDIF
                     ENDIF
                  ENDDO

                  CALL BD_IMBEDDED_BLANK ( JCARD,0,3,0,5,0,5,0,0 ) ! Make sure no imbed blanks in fields 3,5,7 (2,4,6 can be IDOF's)
                  CALL CARD_FLDS_NOT_BLANK ( JCARD,0,0,0,0,0,0,0,9 )
                  CALL CRDERR ( CARD )

                  IF (LARGE_FLD_INP == 'N') THEN
                     CALL NEXTC  ( CARD, ICONT, IERR )
                  ELSE
                     CALL NEXTC2 ( CARD, ICONT, IERR, CHILD )
                     CARD = CHILD
                  ENDIF
                  CALL MKJCARD ( SUBR_NAME, CARD, JCARD )
                  IF (ICONT == 1) THEN
                     CYCLE
                  ELSE
                     EXIT
                  ENDIF
               ENDDO
            ENDIF
         ELSE
            EXIT
         ENDIF
      ENDDO

      IF (UMFND == 'N') THEN
         FATAL_ERR = FATAL_ERR + 1
         WRITE(ERR,1136) JCARD1_PARENT,JCARD2_PARENT
         WRITE(F06,1136) JCARD1_PARENT,JCARD2_PARENT
      ENDIF

! If IDOF fields all have valid DOF no's (or blank), check that the DOF's in those fields specify all 6 DOF no's

      IF (IDOF_ERR == 0) THEN
         DO I=1,NUM_IDOF_FLDS
            DO J=1,JCARD_LEN
               IF (CHR_IDOF(I)(J:J) /= ' ') THEN
                  READ (CHR_IDOF(I)(J:J),'(I8)') INT1
                  IF (DOF_NOS(INT1) == 0) THEN
                     DOF_NOS(INT1) = INT1
                  ELSE
                     FATAL_ERR = FATAL_ERR + 1
                     IDOF_ERR  = IDOF_ERR + 1
                  ENDIF
               ENDIF
            ENDDO
         ENDDO
         DO J=1,6
            IF (DOF_NOS(J) == 0) THEN
               FATAL_ERR = FATAL_ERR + 1
               IDOF_ERR  = IDOF_ERR + 1
            ENDIF
         ENDDO
         IF (IDOF_ERR > 0) THEN
            WRITE(ERR,1143) RELID,(CHR_IDOF(I),I=1,3)
            WRITE(F06,1143) RELID,(CHR_IDOF(I),I=1,3)
            IF (MORE_IDOF == 'N') THEN
               WRITE(ERR,11431)
               WRITE(F06,11431)
            ELSE
               WRITE(ERR,11432) (CHR_IDOF(I),I=4,6)
               WRITE(F06,11432) (CHR_IDOF(I),I=4,6)
            ENDIF
         ENDIF
      ENDIF

      CALL CRDERR ( CARD )

! **********************************************************************
      IF (WRT_LOG >= SUBR_BEGEND) THEN
         CALL OURTIM
         WRITE(F04,9002) SUBR_NAME,TSEC
 9002    FORMAT(1X,A,' END  ',F10.3)
      ENDIF

      RETURN

! **********************************************************************************************************************************
 1124 FORMAT(' *ERROR  1124: INVALID DOF NUMBER IN FIELD ',I3,' ON ',A,' ENTRY WITH ID = ',A                                       &
                    ,/,14X,' MUST BE A COMBINATION OF DIGITS 1-6. HOWEVER, FIELD ',I3, ' HAS: "',A,'"')

 1135 FORMAT(' *ERROR  1135: ',A,' ENTRY WITH ID = ',A,' SHOULD HAVE A MAX OF 1 CONTINUATION ENTRY WITHOUT "UM" IN FIELD 2.'       &
                    ,/,14X,' HOWEVER, CONTINUATION ENTRY NUMBER ',I8,' HAS BEEN FOUND')

 1136 FORMAT(' *ERROR  1136: REQUIRED CONTINUATION FOR ',A,' ID = ',A,' MISSING')

 1143 FORMAT(' *ERROR  1143: RBE1 ELEM NUMBER ',I8,' HAS INCORRECT INDEPENDENT DOFs IDENTIFIED.'                                   &
                    ,/,14X,' THESE (UP TO) 6 FIELDS MUST COMBINE TO SPECIFY DOFs 1,2,3,4,5,6 (EACH DOF ONCE, AND ONLY ONCE).'      &
                    ,/,14X,' HOWEVER, FIELDS 4,6,8 OF THE PARENT     ENTRY HAD: "',A,'", "',A,'", "',A,'"')

11431 FORMAT(          14X,' AND THERE WAS NO CONTINUATION ENTRY WITH ADDITIONAL INDEPENDENT DOFs IDENTIFIED')
11432 FORMAT(          14X,' AND      FIELDS 4,6,8 OF A CONTINUATION ENTRY HAD: "',A,'", "',A,'", "',A,'"')

 1163 FORMAT(' *ERROR  1163: PROGRAMMING ERROR IN SUBROUTINE ',A                                                                   &
                    ,/,14X,' TOO MANY ',A,' ENTRIES; LIMIT = ',I12)

 1181 FORMAT(' *ERROR  1181: RBE1 ELEMENT NUMBER ',A,' HAS NON-BLANK GRID FIELD ',I2,' AND BLANK DOF FIELD ',I2                    &
                    ,/,14X,' THIS PAIR OF GRID/DOF MUST BE EITHER BOTH BLANK OR BOTH NON-BLANK')

 1182 FORMAT(' *ERROR  1182: RBE1 ELEMENT NUMBER ',A,' HAS BLANK GRID FIELD ',I2,' AND NON-BLANK DOF FIELD ',I2                    &
                    ,/,14X,' THIS PAIR OF GRID/DOF MUST BE EITHER BOTH BLANK OR BOTH NON-BLANK')

 1183 FORMAT(' *ERROR  1183: RBE1 ELEMENT NUMBER ',A,' HAS NO CONTINUATION ENTRY. AT LEAST 1 IS REQUIRED')

 1184 FORMAT(' *ERROR  1184: RBE1 ELEMENT NUMBER ',A,' HAS ',I8,' INDEPENDENT DOFs DEFINED. REQUIREMENT IS 6')

! **********************************************************************************************************************************

      END SUBROUTINE BD_RBE1


text:
c     Subroutine to export results from test2 (no source letter)
c     AJ_Kettle, 14Nov2019

      SUBROUTINE export_test2(s_date_st,s_time_st,s_zone_st,
     +    i_bad_sourceletter,
     +    s_directory_output_diagnostics,s_filename_test2,
     +    i_index,s_stnname_single,
     +    i_flag_linenumber,s_flag_varname,
     +    l_rgh_lines,l_lines,
     +    s_vec_date,
     +    s_vec_tmax_c,s_vec_tmin_c,s_vec_tavg_c,
     +    s_vec_prcp_mm,s_vec_snow_mm,s_vec_awnd_ms,
     +    s_vec_tmax_attrib,s_vec_tmin_attrib,
     +    s_vec_tavg_attrib,s_vec_prcp_attrib,
     +    s_vec_snow_attrib,s_vec_awnd_attrib,
     +    s_vec_tmax_sourcelett,s_vec_tmin_sourcelett,
     +    s_vec_tavg_sourcelett,s_vec_prcp_sourcelett,
     +    s_vec_snow_sourcelett,s_vec_awnd_sourcelett)

      IMPLICIT NONE
c************************************************************************
c     Variables passed into subroutine

      CHARACTER(LEN=8)    :: s_date_st
      CHARACTER(LEN=10)   :: s_time_st
      CHARACTER(LEN=5)    :: s_zone_st

      INTEGER             :: i_bad_sourceletter

      CHARACTER(LEN=300)  :: s_directory_output_diagnostics
      CHARACTER(LEN=300)  :: s_filename_test2

      INTEGER             :: i_index
      CHARACTER(LEN=*)    :: s_stnname_single

      INTEGER             :: i_flag_linenumber
      CHARACTER(LEN=4)    :: s_flag_varname

      INTEGER             :: l_rgh_lines
      INTEGER             :: l_lines

      CHARACTER(LEN=*)    :: s_vec_date(l_rgh_lines)

      CHARACTER(LEN=*)    :: s_vec_tmax_c(l_rgh_lines)  !len=10
      CHARACTER(LEN=*)    :: s_vec_tmin_c(l_rgh_lines)
      CHARACTER(LEN=*)    :: s_vec_tavg_c(l_rgh_lines)
      CHARACTER(LEN=*)    :: s_vec_prcp_mm(l_rgh_lines)
      CHARACTER(LEN=*)    :: s_vec_snow_mm(l_rgh_lines)
      CHARACTER(LEN=*)    :: s_vec_awnd_ms(l_rgh_lines)

      CHARACTER(LEN=*)    :: s_vec_tmax_attrib(l_rgh_lines)  !len=10
      CHARACTER(LEN=*)    :: s_vec_tmin_attrib(l_rgh_lines)
      CHARACTER(LEN=*)    :: s_vec_tavg_attrib(l_rgh_lines)
      CHARACTER(LEN=*)    :: s_vec_prcp_attrib(l_rgh_lines)
      CHARACTER(LEN=*)    :: s_vec_snow_attrib(l_rgh_lines)
      CHARACTER(LEN=*)    :: s_vec_awnd_attrib(l_rgh_lines)

      CHARACTER(LEN=*)    :: s_vec_tmax_sourcelett(l_rgh_lines)
      CHARACTER(LEN=*)    :: s_vec_tmin_sourcelett(l_rgh_lines)
      CHARACTER(LEN=*)    :: s_vec_tavg_sourcelett(l_rgh_lines)
      CHARACTER(LEN=*)    :: s_vec_prcp_sourcelett(l_rgh_lines)
      CHARACTER(LEN=*)    :: s_vec_snow_sourcelett(l_rgh_lines)
      CHARACTER(LEN=*)    :: s_vec_awnd_sourcelett(l_rgh_lines)
c*****
c     Variables used within subroutine

      INTEGER             :: i,j,k,ii,jj,kk

      CHARACTER(LEN=300)  :: s_pathandname
      CHARACTER(LEN=3)    :: s_filestatus1
      LOGICAL             :: there

c************************************************************************
      print*,'just entered export_test2'

c      print*,'s_directory_output_diagnostics=',
c     +   s_directory_output_diagnostics
c      print*,'s_filename_test2=',s_filename_test2

c      print*,'i_index=',i_index
c      print*,'s_stnname_single=',TRIM(s_stnname_single)
c      print*,'i_flag_linenumber=',i_flag_linenumber
c      print*,'s_flag_varname=',s_flag_varname

      s_pathandname=
     +  TRIM(s_directory_output_diagnostics)//TRIM(s_filename_test2)
c************************************************************************
      INQUIRE(FILE=TRIM(s_pathandname),
     +   EXIST=there)
      IF (there) THEN
        s_filestatus1='OLD'
      ENDIF
      IF (.NOT.(there)) THEN
        s_filestatus1='NEW'
      ENDIF
c************************************************************************
c     Case of new file
      IF (s_filestatus1.EQ.'NEW') THEN

        OPEN(UNIT=5,
     +   FILE=TRIM(s_pathandname),
     +   FORM='formatted',
     +   STATUS='NEW',ACTION='WRITE')

         WRITE(5,FMT=1002) 'Test 2: no NOAA source letter           '
         WRITE(5,FMT=1002) '                                        '
         WRITE(5,FMT=1006) 'AJ_Kettle',s_date_st,s_time_st,s_zone_st
         WRITE(5,FMT=1002) '                                        '
         WRITE(5,FMT=1002) 'subroutine: export_test2.f              '
         WRITE(5,FMT=1002) '                                        '

 1002    FORMAT(a40)
 1006    FORMAT(t1,a9,t11,a8,t20,a10,t31,a5)

         WRITE(5,FMT=1008) 'Error index=   ',i_bad_sourceletter
         WRITE(5,FMT=1008) 'Station index= ',i_index
 1008    FORMAT(t1,a15,t17,i6)
         WRITE(5,FMT=1009) 'Station name=  ',TRIM(s_stnname_single)
 1009    FORMAT(t1,a15,t17,a11)
         WRITE(5,FMT=1010) 'Variable name= ',s_flag_varname
 1010    FORMAT(t1,a15,t17,a4)
         WRITE(5,FMT=1012) 'Line number=   ',i_flag_linenumber
 1012    FORMAT(t1,a15,t17,i6)
         WRITE(5,FMT=1014) 'Date=          ',
     +     s_vec_date(i_flag_linenumber)
 1014    FORMAT(t1,a15,t17,a7)

c        TMAX
         IF (TRIM(s_flag_varname).EQ.'TMAX') THEN
          WRITE(5,FMT=1016) 'Variable value=',
     +     '<'//TRIM(s_vec_tmax_c(i_flag_linenumber))//'>'
          WRITE(5,FMT=1016) 'Attrib string= ',
     +     '<'//TRIM(s_vec_tmax_attrib(i_flag_linenumber))//'>'
          WRITE(5,FMT=1016) 'Source letter= ',
     +     '<'//TRIM(s_vec_tmax_sourcelett(i_flag_linenumber))//'>'
         ENDIF
c        TMIN
         IF (TRIM(s_flag_varname).EQ.'TMIN') THEN
          WRITE(5,FMT=1016) 'Variable value=',
     +     '<'//TRIM(s_vec_tmin_c(i_flag_linenumber))//'>'
          WRITE(5,FMT=1016) 'Attrib string= ',
     +     '<'//TRIM(s_vec_tmin_attrib(i_flag_linenumber))//'>'
          WRITE(5,FMT=1016) 'Source letter= ',
     +     '<'//TRIM(s_vec_tmin_sourcelett(i_flag_linenumber))//'>'
         ENDIF
c        TAVG
         IF (TRIM(s_flag_varname).EQ.'TAVG') THEN
          WRITE(5,FMT=1016) 'Variable value=',
     +     '<'//TRIM(s_vec_tavg_c(i_flag_linenumber))//'>'
          WRITE(5,FMT=1016) 'Attrib string= ',
     +     '<'//TRIM(s_vec_tavg_attrib(i_flag_linenumber))//'>'
          WRITE(5,FMT=1016) 'Source letter= ',
     +     '<'//TRIM(s_vec_tavg_sourcelett(i_flag_linenumber))//'>'
         ENDIF
c        PRCP
         IF (TRIM(s_flag_varname).EQ.'PRCP') THEN
          WRITE(5,FMT=1016) 'Variable value=',
     +     '<'//TRIM(s_vec_prcp_mm(i_flag_linenumber))//'>'
          WRITE(5,FMT=1016) 'Attrib string= ',
     +     '<'//TRIM(s_vec_prcp_attrib(i_flag_linenumber))//'>'
          WRITE(5,FMT=1016) 'Source letter= ',
     +     '<'//TRIM(s_vec_prcp_sourcelett(i_flag_linenumber))//'>'
         ENDIF
c        SNOW
         IF (TRIM(s_flag_varname).EQ.'SNOW') THEN
          WRITE(5,FMT=1016) 'Variable value=',
     +     '<'//TRIM(s_vec_snow_mm(i_flag_linenumber))//'>'
          WRITE(5,FMT=1016) 'Attrib string= ',
     +     '<'//TRIM(s_vec_snow_attrib(i_flag_linenumber))//'>'
          WRITE(5,FMT=1016) 'Source letter= ',
     +     '<'//TRIM(s_vec_snow_sourcelett(i_flag_linenumber))//'>'
         ENDIF
c        AWND
         IF (TRIM(s_flag_varname).EQ.'AWND') THEN
          WRITE(5,FMT=1016) 'Variable value=',
     +     '<'//TRIM(s_vec_awnd_ms(i_flag_linenumber))//'>'
          WRITE(5,FMT=1016) 'Attrib string= ',
     +     '<'//TRIM(s_vec_awnd_attrib(i_flag_linenumber))//'>'
          WRITE(5,FMT=1016) 'Source letter= ',
     +     '<'//TRIM(s_vec_awnd_ms(i_flag_linenumber))//'>'
         ENDIF

1016     FORMAT(t1,a15,t17,a10)
         WRITE(5,FMT=1002) '                                        '

       CLOSE(UNIT=5)

      ENDIF

c     Case of old file
      IF (s_filestatus1.EQ.'OLD') THEN

        OPEN(UNIT=5,
     +   FILE=TRIM(s_pathandname),
     +   FORM='formatted',
     +   STATUS='OLD',ACTION='WRITE',ACCESS='append')

c 1002    FORMAT(a40)
c 1006    FORMAT(t1,a9,t11,a8,t20,a10,t31,a5)

         WRITE(5,FMT=1008) 'Error index=   ',i_bad_sourceletter
         WRITE(5,FMT=1008) 'Station index= ',i_index
c 1008    FORMAT(t1,a15,t17,i6)
         WRITE(5,FMT=1009) 'Station name=  ',TRIM(s_stnname_single)
c 1009    FORMAT(t1,a15,t17,a11)
         WRITE(5,FMT=1010) 'Variable name= ',s_flag_varname
c 1010    FORMAT(t1,a15,t17,a4)
         WRITE(5,FMT=1012) 'Line number=   ',i_flag_linenumber
c 1012    FORMAT(t1,a15,t17,i6)
         WRITE(5,FMT=1014) 'Date=          ',
     +     s_vec_date(i_flag_linenumber)
c 1014    FORMAT(t1,a15,t17,a7)

c        TMAX
         IF (TRIM(s_flag_varname).EQ.'TMAX') THEN
          WRITE(5,FMT=1016) 'Variable value=',
     +     '<'//TRIM(s_vec_tmax_c(i_flag_linenumber))//'>'
          WRITE(5,FMT=1016) 'Attrib string= ',
     +     '<'//TRIM(s_vec_tmax_attrib(i_flag_linenumber))//'>'
          WRITE(5,FMT=1016) 'Source letter= ',
     +     '<'//TRIM(s_vec_tmax_sourcelett(i_flag_linenumber))//'>'
         ENDIF
c        TMIN
         IF (TRIM(s_flag_varname).EQ.'TMIN') THEN
          WRITE(5,FMT=1016) 'Variable value=',
     +     '<'//TRIM(s_vec_tmin_c(i_flag_linenumber))//'>'
          WRITE(5,FMT=1016) 'Attrib string= ',
     +     '<'//TRIM(s_vec_tmin_attrib(i_flag_linenumber))//'>'
          WRITE(5,FMT=1016) 'Source letter= ',
     +     '<'//TRIM(s_vec_tmin_sourcelett(i_flag_linenumber))//'>'
         ENDIF
c        TAVG
         IF (TRIM(s_flag_varname).EQ.'TAVG') THEN
          WRITE(5,FMT=1016) 'Variable value=',
     +     '<'//TRIM(s_vec_tavg_c(i_flag_linenumber))//'>'
          WRITE(5,FMT=1016) 'Attrib string= ',
     +     '<'//TRIM(s_vec_tavg_attrib(i_flag_linenumber))//'>'
          WRITE(5,FMT=1016) 'Source letter= ',
     +     '<'//TRIM(s_vec_tavg_sourcelett(i_flag_linenumber))//'>'
         ENDIF
c        PRCP
         IF (TRIM(s_flag_varname).EQ.'PRCP') THEN
          WRITE(5,FMT=1016) 'Variable value=',
     +     '<'//TRIM(s_vec_prcp_mm(i_flag_linenumber))//'>'
          WRITE(5,FMT=1016) 'Attrib string= ',
     +     '<'//TRIM(s_vec_prcp_attrib(i_flag_linenumber))//'>'
          WRITE(5,FMT=1016) 'Source letter= ',
     +     '<'//TRIM(s_vec_prcp_sourcelett(i_flag_linenumber))//'>'
         ENDIF
c        SNOW
         IF (TRIM(s_flag_varname).EQ.'SNOW') THEN
          WRITE(5,FMT=1016) 'Variable value=',
     +     '<'//TRIM(s_vec_snow_mm(i_flag_linenumber))//'>'
          WRITE(5,FMT=1016) 'Attrib string= ',
     +     '<'//TRIM(s_vec_snow_attrib(i_flag_linenumber))//'>'
          WRITE(5,FMT=1016) 'Source letter= ',
     +     '<'//TRIM(s_vec_snow_sourcelett(i_flag_linenumber))//'>'
         ENDIF
c        AWND
         IF (TRIM(s_flag_varname).EQ.'AWND') THEN
          WRITE(5,FMT=1016) 'Variable value=',
     +     '<'//TRIM(s_vec_awnd_ms(i_flag_linenumber))//'>'
          WRITE(5,FMT=1016) 'Attrib string= ',
     +     '<'//TRIM(s_vec_awnd_attrib(i_flag_linenumber))//'>'
          WRITE(5,FMT=1016) 'Source letter= ',
     +     '<'//TRIM(s_vec_awnd_ms(i_flag_linenumber))//'>'
         ENDIF

c1016     FORMAT(t1,a15,t17,a10)
         WRITE(5,FMT=1002) '                                        '

       CLOSE(UNIT=5)

      ENDIF

c************************************************************************

      print*,'just leaving export_test2'

      RETURN
      END

ext: f

max_stars_repo_name: glamod/glamod-nuim

max_issues_repo_name: glamod/glamod-nuim

max_forks_repo_name: glamod/glamod-nuim


text:
!
! Copyright 2020-2022 Fluid Numerics LLC
! Author : Joseph Schoonover (joe@fluidnumerics.com)
! Support : support@fluidnumerics.com
!
! //////////////////////////////////////////////////////////////////////////////////////////////// !
MODULE SELF_LinearShallowWater

  USE SELF_Metadata
  USE SELF_Mesh
  USE SELF_MappedData
  USE SELF_Model

  TYPE,EXTENDS(Model2D) :: LinearShallowWater
    !! iVar = 1 ~> u velocity component
    !! iVar = 2 ~> v velocity component
    !! iVar = 3 ~> free surface height
    REAL(prec) :: fCori ! coriolis parameter ( 1/s )
    REAL(prec) :: g     ! gravity ( m/s^2)
    REAL(prec) :: H     ! fluid thickness ( m )

    CONTAINS

    ! Overridden Methods
    PROCEDURE :: Init => Init_LinearShallowWater
    PROCEDURE :: CalculateEntropy => CalculateEntropy_LinearShallowWater

    ! Concretized Methods
    PROCEDURE :: SourceMethod => Source_LinearShallowWater
    PROCEDURE :: FluxMethod => Flux_LinearShallowWater
    PROCEDURE :: RiemannSolver => RiemannSolver_LinearShallowWater
    PROCEDURE :: SetBoundaryCondition => SetBoundaryCondition_LinearShallowWater

  END TYPE LinearShallowWater

  INTERFACE
    SUBROUTINE SetBoundaryCondition_LinearShallowWater_gpu_wrapper(solution, extBoundary, nHat, sideInfo, N, nVar, nEl) &
      bind(c,name="SetBoundaryCondition_LinearShallowWater_gpu_wrapper")
      USE iso_c_binding
      USE SELF_Constants
      IMPLICIT NONE
      TYPE(c_ptr) :: solution, extBoundary, nHat, sideInfo
      INTEGER(C_INT),VALUE :: N,nVar,nEl
    END SUBROUTINE SetBoundaryCondition_LinearShallowWater_gpu_wrapper
  END INTERFACE

  INTERFACE
    SUBROUTINE Source_LinearShallowWater_gpu_wrapper(source, solution, f, N, nVar, nEl) &
      bind(c,name="Source_LinearShallowWater_gpu_wrapper")
      USE iso_c_binding
      USE SELF_Constants
      IMPLICIT NONE
      TYPE(c_ptr) :: source, solution
      INTEGER(C_INT),VALUE :: N,nVar,nEl
      REAL(c_prec),VALUE :: f
    END SUBROUTINE Source_LinearShallowWater_gpu_wrapper
  END INTERFACE

  INTERFACE
    SUBROUTINE Flux_LinearShallowWater_gpu_wrapper(flux, solution, g, H, N, nVar, nEl) &
      bind(c,name="Flux_LinearShallowWater_gpu_wrapper")
      USE iso_c_binding
      USE SELF_Constants
      IMPLICIT NONE
      TYPE(c_ptr) :: flux, solution
      INTEGER(C_INT),VALUE :: N,nVar,nEl
      REAL(c_prec),VALUE :: g, H
    END SUBROUTINE Flux_LinearShallowWater_gpu_wrapper
  END INTERFACE

  INTERFACE
    SUBROUTINE RiemannSolver_LinearShallowWater_gpu_wrapper(flux, solution, extBoundary, nHat, nScale, g, H, N, nVar, nEl) &
      bind(c,name="RiemannSolver_LinearShallowWater_gpu_wrapper")
      USE iso_c_binding
      USE SELF_Constants
      IMPLICIT NONE
      TYPE(c_ptr) :: flux, solution, extBoundary, nHat, nScale
      INTEGER(C_INT),VALUE :: N,nVar,nEl
      REAL(c_prec),VALUE :: g, H
    END SUBROUTINE RiemannSolver_LinearShallowWater_gpu_wrapper
  END INTERFACE

CONTAINS

  SUBROUTINE Init_LinearShallowWater(this,nvar,mesh,geometry,decomp)
    IMPLICIT NONE
    CLASS(LinearShallowWater),INTENT(out) :: this
    INTEGER,INTENT(in) :: nvar
    TYPE(Mesh2D),INTENT(in),TARGET :: mesh
    TYPE(SEMQuad),INTENT(in),TARGET :: geometry
    TYPE(MPILayer),INTENT(in),TARGET :: decomp
    ! Local
    INTEGER :: ivar
    CHARACTER(LEN=3) :: ivarChar
    CHARACTER(LEN=25) :: varname
    INTEGER :: nvarloc

    ! Ensure that the number of variables is 3
    ! nvar is unused in this class extension
    nvarloc = 3

    this % decomp => decomp
    this % mesh => mesh
    this % geometry => geometry
    this % gpuAccel = .FALSE.
    this % fluxDivMethod = SELF_CONSERVATIVE_FLUX
    this % g = 1.0_prec
    this % H = 1.0_prec
    this % fCori = 0.0_prec

    CALL this % solution % Init(geometry % x % interp,nvarloc,this % mesh % nElem)
    CALL this % workSol % Init(geometry % x % interp,nVar,this % mesh % nElem)
    CALL this % velocity % Init(geometry % x % interp,1,this % mesh % nElem)
    CALL this % compVelocity % Init(geometry % x % interp,1,this % mesh % nElem)
    CALL this % dSdt % Init(geometry % x % interp,nvarloc,this % mesh % nElem)
    CALL this % solutionGradient % Init(geometry % x % interp,nvarloc,this % mesh % nElem)
    CALL this % flux % Init(geometry % x % interp,nvarloc,this % mesh % nElem)
    CALL this % source % Init(geometry % x % interp,nvarloc,this % mesh % nElem)
    CALL this % fluxDivergence % Init(geometry % x % interp,nvarloc,this % mesh % nElem)

    ! First three variables are treated as u, v, eta
    ! Any additional are treated as passive tracers
    CALL this % solution % SetName(1,"u")
    CALL this % solution % SetUnits(1,"m/s")
    CALL this % solution % SetDescription(1,"x-component of the barotropic velocity field")

    CALL this % solution % SetName(2,"v")
    CALL this % solution % SetUnits(2,"m/s")
    CALL this % solution % SetDescription(2,"y-component of the barotropic velocity field")

    CALL this % solution % SetName(3,"eta")
    CALL this % solution % SetUnits(3,"m")
    CALL this % solution % SetDescription(3,"Free surface height anomaly")

  END SUBROUTINE Init_LinearShallowWater

  SUBROUTINE CalculateEntropy_LinearShallowWater(this)
  !! Base method for calculating entropy of a model
  !! Calculates the entropy as the integration of the
  !! squared tracer over the domain
    IMPLICIT NONE
    CLASS(LinearShallowWater), INTENT(inout) :: this
    ! Local
    INTEGER :: i, j, iVar, iEl
    REAL(prec) :: Jacobian, u, v, eta
    REAL(prec) :: wi,wj

    ! TO DO : GPU reduction

    this % entropy = 0.0_prec

    DO iEl = 1, this % geometry % x % nElem
        DO j = 0, this % geometry % x % interp % N
          DO i = 0, this % geometry % x % interp % N

            ! Coordinate mapping Jacobian
            Jacobian = this % geometry % J % interior % hostData(i,j,1,iEl)

            ! Quadrature weights
            wi = this % geometry % x % interp % qWeights % hostData(i)
            wj = this % geometry % x % interp % qWeights % hostData(j)

            ! Solution
            u = this % solution % interior % hostData(i,j,1,iEl)
            v = this % solution % interior % hostData(i,j,2,iEl)
            eta = this % solution % interior % hostData(i,j,3,iEl)

            this % entropy = this % entropy + 0.5_prec*( this % H*(u*u + v*v) + this % g*eta*eta )*Jacobian*wi*wj

          ENDDO
        ENDDO
    ENDDO

  END SUBROUTINE CalculateEntropy_LinearShallowWater

  SUBROUTINE SetBoundaryCondition_LinearShallowWater(this)
    IMPLICIT NONE
    CLASS(LinearShallowWater),INTENT(inout) :: this
    ! Local
    INTEGER :: iEl, iSide, i
    INTEGER :: bcid, e2
    REAL(prec) :: u, v, nhat(1:2)

    IF( this % gpuAccel )THEN

      CALL SetBoundaryCondition_LinearShallowWater_gpu_wrapper( this % solution % boundary % deviceData,&
            this % solution % extBoundary % deviceData, &
            this % geometry % nHat % boundary % deviceData, &
            this % mesh % sideInfo % deviceData, &
            this % solution % interp % N, &
            this % solution % nVar, &
            this % solution % nElem)

    ELSE

      DO iEl = 1, this % solution % nElem
        DO iSide = 1, 4
            DO i = 0, this % solution % interp % N

              bcid = this % mesh % sideInfo % hostData(5,iSide,iEl) ! Boundary Condition ID
              e2 = this % mesh % sideInfo % hostData(3,iSide,iEl) ! Neighboring Element ID
              IF( e2 == 0 )THEN
                IF( bcid == SELF_BC_RADIATION )THEN

                  this % solution % extBoundary % hostData(i,1,iSide,iEl) = 0.0_prec
                  this % solution % extBoundary % hostData(i,2,iSide,iEl) = 0.0_prec
                  this % solution % extBoundary % hostData(i,3,iSide,iEl) = 0.0_prec

                ELSEIF( bcid == SELF_BC_NONORMALFLOW )THEN

                  nhat(1:2) = this % geometry % nHat % boundary % hostData(1:2,i,1,iSide,iEl)
                  u = this % solution % boundary % hostData(i,1,iSide,iEl)
                  v = this % solution % boundary % hostData(i,2,iSide,iEl)
                  this % solution % extBoundary % hostData(i,1,iSide,iEl) = (nhat(2)**2 - nhat(1)**2)*u - 2.0_prec*nhat(1)*nhat(2)*v
                  this % solution % extBoundary % hostData(i,2,iSide,iEl) = (nhat(1)**2 - nhat(2)**2)*v - 2.0_prec*nhat(1)*nhat(2)*u
                  this % solution % extBoundary % hostData(i,3,iSide,iEl) = this % solution % boundary % hostData(i,3,iSide,iEl)

                ELSE ! Default boundary condition is radiation

                  this % solution % extBoundary % hostData(i,1,iSide,iEl) = 0.0_prec
                  this % solution % extBoundary % hostData(i,2,iSide,iEl) = 0.0_prec
                  this % solution % extBoundary % hostData(i,3,iSide,iEl) = 0.0_prec

                ENDIF
              ENDIF

            ENDDO
        ENDDO
      ENDDO

    ENDIF

  END SUBROUTINE SetBoundaryCondition_LinearShallowWater

  SUBROUTINE Source_LinearShallowWater(this)
    IMPLICIT NONE
    CLASS(LinearShallowWater),INTENT(inout) :: this
    ! Local
    INTEGER :: i,j,iEl,iVar

    IF( this % gpuAccel )THEN

      CALL Source_LinearShallowWater_gpu_wrapper(this % source % interior % deviceData,&
             this % solution % interior % deviceData, &
             this % fCori,&
             this % source % interp % N, &
             this % solution % nVar, &
             this % solution % nElem )

    ELSE

      DO iEl = 1, this % source % nElem
        DO j = 0, this % source % interp % N
          DO i = 0, this % source % interp % N

            this % source % interior % hostData(i,j,1,iEl) = -this % fCori*this % solution % interior % hostData(i,j,2,iEl)
            this % source % interior % hostData(i,j,2,iEl) = this % fCori*this % solution % interior % hostData(i,j,1,iEl)
            this % source % interior % hostData(i,j,3,iEl) = 0.0_prec

          ENDDO
        ENDDO
      ENDDO

    ENDIF

  END SUBROUTINE Source_LinearShallowWater

  SUBROUTINE Flux_LinearShallowWater(this)
    IMPLICIT NONE
    CLASS(LinearShallowWater),INTENT(inout) :: this
    ! Local
    INTEGER :: i,j,iEl,iVar

    IF( this % gpuAccel )THEN

      CALL Flux_LinearShallowWater_gpu_wrapper(this % flux % interior % deviceData,&
                                               this % solution % interior % deviceData, &
                                               this % g, this % H, this % solution % interp % N, &
                                               this % solution % nVar, this % solution % nElem)

    ELSE
      DO iEl = 1, this % solution % nElem
        DO iVar = 1, this % solution % nVar
          DO j = 0, this % solution % interp % N
            DO i = 0, this % solution % interp % N

              IF ( iVar == 1 )THEN ! u-velocity
                this % flux % interior % hostData(1,i,j,iVar,iEl) = &
                      this % g*this % solution % interior % hostData(i,j,3,iEl)

                this % flux % interior % hostData(2,i,j,iVar,iEl) = 0.0_prec

              ELSEIF ( iVar == 2 )THEN ! v-velocity

                this % flux % interior % hostData(1,i,j,iVar,iEl) = 0.0_prec

                this % flux % interior % hostData(2,i,j,iVar,iEl) = &
                      this % g*this % solution % interior % hostData(i,j,3,iEl)

              ELSEIF ( iVar == 3 )THEN ! free surface height
                this % flux % interior % hostData(1,i,j,iVar,iEl) = &
                      this % H*this % solution % interior % hostData(i,j,1,iEl)

                this % flux % interior % hostData(2,i,j,iVar,iEl) = &
                      this % H*this % solution % interior % hostData(i,j,2,iEl)
              ENDIF

            ENDDO
          ENDDO
        ENDDO
      ENDDO
    ENDIF

  END SUBROUTINE Flux_LinearShallowWater

  SUBROUTINE RiemannSolver_LinearShallowWater(this)
    IMPLICIT NONE
    CLASS(LinearShallowWater),INTENT(inout) :: this
    ! Local
    INTEGER :: i,iSide,iEl
    REAL(prec) :: nhat(1:2), nmag
    REAL(prec) :: c, unL, unR, etaL, etaR, wL, wR

    IF( this % gpuAccel )THEN

      CALL RiemannSolver_LinearShallowWater_gpu_wrapper(this % flux % boundaryNormal % deviceData, &
             this % solution % boundary % deviceData, &
             this % solution % extBoundary % deviceData, &
             this % geometry % nHat % boundary % deviceData, &
             this % geometry % nScale % boundary % deviceData, &
             this % g, this % H, &
             this % solution % interp % N, &
             this % solution % nVar, &
             this % solution % nElem)

    ELSE

      DO iEl = 1, this % solution % nElem
        DO iSide = 1, 4
          DO i = 0, this % solution % interp % N

             ! Get the boundary normals on cell edges from the mesh geometry
             nhat(1:2) = this % geometry % nHat % boundary % hostData(1:2,i,1,iSide,iEl)
             nmag = this % geometry % nScale % boundary % hostData(i,1,iSide,iEl)
             c = sqrt( this % g * this % H )

             ! Calculate the normal velocity at the cell edges
             unL = this % solution % boundary % hostData(i,1,iSide,iEl)*nHat(1)+&
                   this % solution % boundary % hostData(i,2,iSide,iEl)*nHat(2)

             unR = this % solution % extBoundary % hostData(i,1,iSide,iEl)*nHat(1)+&
                   this % solution % extBoundary % hostData(i,2,iSide,iEl)*nHat(2)

             etaL = this % solution % boundary % hostData(i,3,iSide,iEl)
             etaR = this % solution % extBoundary % hostData(i,3,iSide,iEl)

             ! Pull external and internal state for the Riemann Solver (Lax-Friedrichs)
             wL = 0.5_prec*(unL/this % g + etaL/c)
             wR = 0.5_prec*(unR/this % g - etaR/c)

             this % flux % boundaryNormal % hostData(i,1,iSide,iEl) = this % g*c*( wL - wR )*nHat(1)*nmag
             this % flux % boundaryNormal % hostData(i,2,iSide,iEl) = this % g*c*( wL - wR )*nHat(2)*nmag
             this % flux % boundaryNormal % hostData(i,3,iSide,iEl) = c*c*( wL + wR )*nmag

          ENDDO
        ENDDO
      ENDDO

    ENDIF

  END SUBROUTINE RiemannSolver_LinearShallowWater

END MODULE SELF_LinearShallowWater

max_stars_repo_name: FluidNumerics/self-fluids

max_issues_repo_name: FluidNumerics/self-fluids

max_forks_repo_name: FluidNumerics/SELF-Fluids


text:
!##############################################################################
!# ****************************************************************************
!# <name> physics </name>
!# ****************************************************************************
!#
!# <purpose>
!# Encapsules the physics of the problem.
!# </purpose>
!##############################################################################

module physics

  use fsystem
  use storage

  implicit none

  ! Encapsules the physics of the problem.
  type t_physics

    ! Type of the underlying operator.
    ! optimal distributed control of the ...
    ! =0: 2D heat equation
    ! =1: 2D Stokes
    ! =2: 1D heat equation
    integer :: cequation = 0

    ! Type of reference problem. Depending on cequation.
    ! creferenceproblem = 0: no chosen reference problem.
    !
    ! cequation = 0/2: (1D/2D heat equation)
    !
    !   creferenceproblem = 1:
    !      y      = t^2 x
    !      lambda = -2 t x
    !
    !   creferenceproblem = 2:
    !      y      = t^2 (1-t)^2 x
    !      lambda = t (1-t) x
    !
    !   creferenceproblem = 3:
    !      y      = t^2 (1-t)^2 x
    !      lambda = t (1-t)^2 x
    !
    !   creferenceproblem = 4:
    !      y      = t^2 (1-t)^2 x
    !      lambda = passend, so dass f=0; incl. alpha.
    !
    ! cequation = 1: (2D Stokes equation)
    !
    !   creferenceproblem = 1:
    !      y      = t^2 (x1,-x2)
    !      lambda = -2 t (x1,-x2)
    !
    !   creferenceproblem = 2:
    !      y      = t^2 (1-t)^2 x2(1-x2)
    !      lambda = passend incl. alpha
    !
    !   creferenceproblem = 3:
    !      y      = t^2 (x1^2 x2,-x1 x2^2)
    !      lambda = passend incl. alpha
    !
    integer :: creferenceproblem = 0

    ! Equation-specific parameter
    real(DP) :: dpar

    ! Viscosity parameter
    real(DP) :: dviscosity = 0.0_DP

    ! Relaxation parameter ALPHA of the optimal control
    real(DP) :: doptControlAlpha = 0.0_DP

    ! Relaxation parameter GAMMA for the terminal condition
    real(DP) :: doptControlGamma = 0.0_DP

    ! Couple the dual equation to the primal
    real(DP) :: dcoupleDualToPrimal = 1.0_DP

    ! Couple the primal equation to the dual
    real(DP) :: dcouplePrimalToDual = 1.0_DP

    ! Couple the terminal condition to the dual solution.
    real(DP) :: dcoupleTermCond = 1.0_DP

    ! Minimum time
    real(DP) :: dtimemin = 0.0_DP

    ! Maximum time
    real(DP) :: dtimemax = 1.0_DP

  end type

end module

max_stars_repo_licenses: Intel Unlicense

max_issues_repo_licenses: Intel Unlicense

max_forks_repo_licenses: Intel Unlicense


text:
program tester

use mod_global
use mod_kinds
use mod_QCsymmetry
use mod_QCcrystallography
use mod_io

IMPLICIT NONE

type(QCSpaceGroup_T)        :: QCSG
type(IO_T)                  :: Message
integer(kind=irg)           :: i

do i=6,11
  if (i.ne.7) then
    QCSG = QCspacegroup_T( nD = 3, QCtype = 'Ico')
    call QCSG%setSGnum(i)
! call QCcell_icosahedral%setMetricParametersQC()
    call QCSG%GenerateQCSymmetry(dopg=.FALSE.)
    write (*,*) ' order of icosahedral space group ', i, ' equals ',QCSG%getnsym()
  end if
end do

end program

max_stars_repo_name: EMsoft-org/EMsoftOO

max_issues_repo_name: EMsoft-org/EMsoftOO

max_forks_repo_name: EMsoft-org/EMsoftOO


!*==dlarrd.f90  processed by SPAG 7.51RB at 20:08 on  3 Mar 2022
!> \brief \b DLARRD computes the eigenvalues of a symmetric tridiagonal matrix to suitable accuracy.
!
!  =========== DOCUMENTATION ===========
!
! Online html documentation available at
!            http://www.netlib.org/lapack/explore-html/
!
!> \htmlonly
!> Download DLARRD + dependencies
!>
!> [TGZ]
!>
!> [ZIP]
!>
!> [TXT]
!> \endhtmlonly
!
!  Definition:
!  ===========
!
!       SUBROUTINE DLARRD( RANGE, ORDER, N, VL, VU, IL, IU, GERS,
!                           RELTOL, D, E, E2, PIVMIN, NSPLIT, ISPLIT,
!                           M, W, WERR, WL, WU, IBLOCK, INDEXW,
!                           WORK, IWORK, INFO )
!
!       .. Scalar Arguments ..
!       CHARACTER          ORDER, RANGE
!       INTEGER            IL, INFO, IU, M, N, NSPLIT
!       DOUBLE PRECISION    PIVMIN, RELTOL, VL, VU, WL, WU
!       ..
!       .. Array Arguments ..
!       INTEGER            IBLOCK( * ), INDEXW( * ),
!      $                   ISPLIT( * ), IWORK( * )
!       DOUBLE PRECISION   D( * ), E( * ), E2( * ),
!      $                   GERS( * ), W( * ), WERR( * ), WORK( * )
!       ..
!
!
!> \par Purpose:
!  =============
!>
!> \verbatim
!>
!> DLARRD computes the eigenvalues of a symmetric tridiagonal
!> matrix T to suitable accuracy. This is an auxiliary code to be
!> called from DSTEMR.
!> The user may ask for all eigenvalues, all eigenvalues
!> in the half-open interval (VL, VU], or the IL-th through IU-th
!> eigenvalues.
!>
!> To avoid overflow, the matrix must be scaled so that its
!> largest element is no greater than overflow**(1/2) * underflow**(1/4) in absolute value, and for greatest
!> accuracy, it should not be much smaller than that.
!>
!> See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal
!> Matrix", Report CS41, Computer Science Dept., Stanford
!> University, July 21, 1966.
!> \endverbatim
!
!  Arguments:
!  ==========
!
!> \param[in] RANGE
!> \verbatim
!>          RANGE is CHARACTER*1
!>          = 'A': ("All")   all eigenvalues will be found.
!>          = 'V': ("Value") all eigenvalues in the half-open interval
!>                           (VL, VU] will be found.
!>          = 'I': ("Index") the IL-th through IU-th eigenvalues (of the
!>                           entire matrix) will be found.
!> \endverbatim
!>
!> \param[in] ORDER
!> \verbatim
!>          ORDER is CHARACTER*1
!>          = 'B': ("By Block") the eigenvalues will be grouped by
!>                              split-off block (see IBLOCK, ISPLIT) and
!>                              ordered from smallest to largest within
!>                              the block.
!>          = 'E': ("Entire matrix")
!>                              the eigenvalues for the entire matrix
!>                              will be ordered from smallest to
!>                              largest.
!> \endverbatim
!>
!> \param[in] N
!> \verbatim
!>          N is INTEGER
!>          The order of the tridiagonal matrix T.  N >= 0.
!> \endverbatim
!>
!> \param[in] VL
!> \verbatim
!>          VL is DOUBLE PRECISION
!>          If RANGE='V', the lower bound of the interval to
!>          be searched for eigenvalues.  Eigenvalues less than or equal
!>          to VL, or greater than VU, will not be returned.  VL < VU.
!>          Not referenced if RANGE = 'A' or 'I'.
!> \endverbatim
!>
!> \param[in] VU
!> \verbatim
!>          VU is DOUBLE PRECISION
!>          If RANGE='V', the upper bound of the interval to
!>          be searched for eigenvalues.  Eigenvalues less than or equal
!>          to VL, or greater than VU, will not be returned.  VL < VU.
!>          Not referenced if RANGE = 'A' or 'I'.
!> \endverbatim
!>
!> \param[in] IL
!> \verbatim
!>          IL is INTEGER
!>          If RANGE='I', the index of the
!>          smallest eigenvalue to be returned.
!>          1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
!>          Not referenced if RANGE = 'A' or 'V'.
!> \endverbatim
!>
!> \param[in] IU
!> \verbatim
!>          IU is INTEGER
!>          If RANGE='I', the index of the
!>          largest eigenvalue to be returned.
!>          1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
!>          Not referenced if RANGE = 'A' or 'V'.
!> \endverbatim
!>
!> \param[in] GERS
!> \verbatim
!>          GERS is DOUBLE PRECISION array, dimension (2*N)
!>          The N Gerschgorin intervals (the i-th Gerschgorin interval
!>          is (GERS(2*i-1), GERS(2*i)).
!> \endverbatim
!>
!> \param[in] RELTOL
!> \verbatim
!>          RELTOL is DOUBLE PRECISION
!>          The minimum relative width of an interval.  When an interval
!>          is narrower than RELTOL times the larger (in
!>          magnitude) endpoint, then it is considered to be
!>          sufficiently small, i.e., converged.  Note: this should
!>          always be at least radix*machine epsilon.
!> \endverbatim
!>
!> \param[in] D
!> \verbatim
!>          D is DOUBLE PRECISION array, dimension (N)
!>          The n diagonal elements of the tridiagonal matrix T.
!> \endverbatim
!>
!> \param[in] E
!> \verbatim
!>          E is DOUBLE PRECISION array, dimension (N-1)
!>          The (n-1) off-diagonal elements of the tridiagonal matrix T.
!> \endverbatim
!>
!> \param[in] E2
!> \verbatim
!>          E2 is DOUBLE PRECISION array, dimension (N-1)
!>          The (n-1) squared off-diagonal elements of the tridiagonal matrix T.
!> \endverbatim
!>
!> \param[in] PIVMIN
!> \verbatim
!>          PIVMIN is DOUBLE PRECISION
!>          The minimum pivot allowed in the Sturm sequence for T.
!> \endverbatim
!>
!> \param[in] NSPLIT
!> \verbatim
!>          NSPLIT is INTEGER
!>          The number of diagonal blocks in the matrix T.
!>          1 <= NSPLIT <= N.
!> \endverbatim
!>
!> \param[in] ISPLIT
!> \verbatim
!>          ISPLIT is INTEGER array, dimension (N)
!>          The splitting points, at which T breaks up into submatrices.
!>          The first submatrix consists of rows/columns 1 to ISPLIT(1),
!>          the second of rows/columns ISPLIT(1)+1 through ISPLIT(2),
!>          etc., and the NSPLIT-th consists of rows/columns
!>          ISPLIT(NSPLIT-1)+1 through ISPLIT(NSPLIT)=N.
!>          (Only the first NSPLIT elements will actually be used, but
!>          since the user cannot know a priori what value NSPLIT will
!>          have, N words must be reserved for ISPLIT.)
!> \endverbatim
!>
!> \param[out] M
!> \verbatim
!>          M is INTEGER
!>          The actual number of eigenvalues found. 0 <= M <= N.
!>          (See also the description of INFO=2,3.)
!> \endverbatim
!>
!> \param[out] W
!> \verbatim
!>          W is DOUBLE PRECISION array, dimension (N)
!>          On exit, the first M elements of W will contain the
!>          eigenvalue approximations. DLARRD computes an interval
!>          I_j = (a_j, b_j] that includes eigenvalue j. The eigenvalue
!>          approximation is given as the interval midpoint
!>          W(j)= ( a_j + b_j)/2. The corresponding error is bounded by
!>          WERR(j) = abs( a_j - b_j)/2
!> \endverbatim
!>
!> \param[out] WERR
!> \verbatim
!>          WERR is DOUBLE PRECISION array, dimension (N)
!>          The error bound on the corresponding eigenvalue approximation
!>          in W.
!> \endverbatim
!>
!> \param[out] WL
!> \verbatim
!>          WL is DOUBLE PRECISION
!> \endverbatim
!>
!> \param[out] WU
!> \verbatim
!>          WU is DOUBLE PRECISION
!>          The interval (WL, WU] contains all the wanted eigenvalues.
!>          If RANGE='V', then WL=VL and WU=VU.
!>          If RANGE='A', then WL and WU are the global Gerschgorin bounds
!>                        on the spectrum.
!>          If RANGE='I', then WL and WU are computed by DLAEBZ from the
!>                        index range specified.
!> \endverbatim
!>
!> \param[out] IBLOCK
!> \verbatim
!>          IBLOCK is INTEGER array, dimension (N)
!>          At each row/column j where E(j) is zero or small, the
!>          matrix T is considered to split into a block diagonal
!>          matrix.  On exit, if INFO = 0, IBLOCK(i) specifies to which
!>          block (from 1 to the number of blocks) the eigenvalue W(i)
!>          belongs.  (DLARRD may use the remaining N-M elements as
!>          workspace.)
!> \endverbatim
!>
!> \param[out] INDEXW
!> \verbatim
!>          INDEXW is INTEGER array, dimension (N)
!>          The indices of the eigenvalues within each block (submatrix);
!>          for example, INDEXW(i)= j and IBLOCK(i)=k imply that the
!>          i-th eigenvalue W(i) is the j-th eigenvalue in block k.
!> \endverbatim
!>
!> \param[out] WORK
!> \verbatim
!>          WORK is DOUBLE PRECISION array, dimension (4*N)
!> \endverbatim
!>
!> \param[out] IWORK
!> \verbatim
!>          IWORK is INTEGER array, dimension (3*N)
!> \endverbatim
!>
!> \param[out] INFO
!> \verbatim
!>          INFO is INTEGER
!>          = 0:  successful exit
!>          < 0:  if INFO = -i, the i-th argument had an illegal value
!>          > 0:  some or all of the eigenvalues failed to converge or
!>                were not computed:
!>                =1 or 3: Bisection failed to converge for some
!>                        eigenvalues; these eigenvalues are flagged by a
!>                        negative block number.  The effect is that the
!>                        eigenvalues may not be as accurate as the
!>                        absolute and relative tolerances.  This is
!>                        generally caused by unexpectedly inaccurate
!>                        arithmetic.
!>                =2 or 3: RANGE='I' only: Not all of the eigenvalues
!>                        IL:IU were found.
!>                        Effect: M < IU+1-IL
!>                        Cause:  non-monotonic arithmetic, causing the
!>                                Sturm sequence to be non-monotonic.
!>                        Cure:   recalculate, using RANGE='A', and pick
!>                                out eigenvalues IL:IU.  In some cases,
!>                                increasing the PARAMETER "FUDGE" may
!>                                make things work.
!>                = 4:    RANGE='I', and the Gershgorin interval
!>                        initially used was too small.  No eigenvalues
!>                        were computed.
!>                        Probable cause: your machine has sloppy
!>                                        floating-point arithmetic.
!>                        Cure: Increase the PARAMETER "FUDGE",
!>                              recompile, and try again.
!> \endverbatim
!
!> \par Internal Parameters:
!  =========================
!>
!> \verbatim
!>  FUDGE   DOUBLE PRECISION, default = 2
!>          A "fudge factor" to widen the Gershgorin intervals.  Ideally,
!>          a value of 1 should work, but on machines with sloppy
!>          arithmetic, this needs to be larger.  The default for
!>          publicly released versions should be large enough to handle
!>          the worst machine around.  Note that this has no effect
!>          on accuracy of the solution.
!> \endverbatim
!>
!> \par Contributors:
!  ==================
!>
!>     W. Kahan, University of California, Berkeley, USA \n
!>     Beresford Parlett, University of California, Berkeley, USA \n
!>     Jim Demmel, University of California, Berkeley, USA \n
!>     Inderjit Dhillon, University of Texas, Austin, USA \n
!>     Osni Marques, LBNL/NERSC, USA \n
!>     Christof Voemel, University of California, Berkeley, USA \n
!
!  Authors:
!  ========
!
!> \author Univ. of Tennessee
!> \author Univ. of California Berkeley
!> \author Univ. of Colorado Denver
!> \author NAG Ltd.
!
!> \date June 2016
!
!> \ingroup OTHERauxiliary
!
!  =====================================================================
      SUBROUTINE DLARRD(Range,Order,N,Vl,Vu,Il,Iu,Gers,Reltol,D,E,E2,   &
     &                  Pivmin,Nsplit,Isplit,M,W,Werr,Wl,Wu,Iblock,     &
     &                  Indexw,Work,Iwork,Info)
      IMPLICIT NONE
!*--DLARRD332
!
!  -- LAPACK auxiliary routine (version 3.7.1) --
!  -- LAPACK is a software package provided by Univ. of Tennessee,    --
!  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
!     June 2016
!
!     .. Scalar Arguments ..
      CHARACTER Order , Range
      INTEGER Il , Info , Iu , M , N , Nsplit
      DOUBLE PRECISION Pivmin , Reltol , Vl , Vu , Wl , Wu
!     ..
!     .. Array Arguments ..
      INTEGER Iblock(*) , Indexw(*) , Isplit(*) , Iwork(*)
      DOUBLE PRECISION D(*) , E(*) , E2(*) , Gers(*) , W(*) , Werr(*) , &
     &                 Work(*)
!     ..
!
!  =====================================================================
!
!     .. Parameters ..
      DOUBLE PRECISION ZERO , ONE , TWO , HALF , FUDGE
      PARAMETER (ZERO=0.0D0,ONE=1.0D0,TWO=2.0D0,HALF=ONE/TWO,FUDGE=TWO)
      INTEGER ALLRNG , VALRNG , INDRNG
      PARAMETER (ALLRNG=1,VALRNG=2,INDRNG=3)
!     ..
!     .. Local Scalars ..
      LOGICAL ncnvrg , toofew
      INTEGER i , ib , ibegin , idiscl , idiscu , ie , iend , iinfo ,   &
     &        im , in , ioff , iout , irange , itmax , itmp1 , itmp2 ,  &
     &        iw , iwoff , j , jblk , jdisc , je , jee , nb , nwl , nwu
      DOUBLE PRECISION atoli , eps , gl , gu , rtoli , tmp1 , tmp2 ,    &
     &                 tnorm , uflow , wkill , wlu , wul

!     ..
!     .. Local Arrays ..
      INTEGER idumma(1)
!     ..
!     .. External Functions ..
      LOGICAL LSAME
      INTEGER ILAENV
      DOUBLE PRECISION DLAMCH
      EXTERNAL LSAME , ILAENV , DLAMCH
!     ..
!     .. External Subroutines ..
      EXTERNAL DLAEBZ
!     ..
!     .. Intrinsic Functions ..
      INTRINSIC ABS , INT , LOG , MAX , MIN
!     ..
!     .. Executable Statements ..
!
      Info = 0
!
!     Quick return if possible
!
      IF ( N<=0 ) RETURN
!
!     Decode RANGE
!
      IF ( LSAME(Range,'A') ) THEN
         irange = ALLRNG
      ELSEIF ( LSAME(Range,'V') ) THEN
         irange = VALRNG
      ELSEIF ( LSAME(Range,'I') ) THEN
         irange = INDRNG
      ELSE
         irange = 0
      ENDIF
!
!     Check for Errors
!
      IF ( irange<=0 ) THEN
         Info = -1
      ELSEIF ( .NOT.(LSAME(Order,'B') .OR. LSAME(Order,'E')) ) THEN
         Info = -2
      ELSEIF ( N<0 ) THEN
         Info = -3
      ELSEIF ( irange==VALRNG ) THEN
         IF ( Vl>=Vu ) Info = -5
      ELSEIF ( irange==INDRNG .AND. (Il<1 .OR. Il>MAX(1,N)) ) THEN
         Info = -6
      ELSEIF ( irange==INDRNG .AND. (Iu<MIN(N,Il) .OR. Iu>N) ) THEN
         Info = -7
      ENDIF
!
      IF ( Info/=0 ) RETURN

!     Initialize error flags
      Info = 0
      ncnvrg = .FALSE.
      toofew = .FALSE.

!     Quick return if possible
      M = 0
      IF ( N==0 ) RETURN

!     Simplification:
      IF ( irange==INDRNG .AND. Il==1 .AND. Iu==N ) irange = 1

!     Get machine constants
      eps = DLAMCH('P')
      uflow = DLAMCH('U')

!     Special Case when N=1
!     Treat case of 1x1 matrix for quick return
      IF ( N==1 ) THEN
         IF ( (irange==ALLRNG) .OR.                                     &
     &        ((irange==VALRNG) .AND. (D(1)>Vl) .AND. (D(1)<=Vu)) .OR.  &
     &        ((irange==INDRNG) .AND. (Il==1) .AND. (Iu==1)) ) THEN
            M = 1
            W(1) = D(1)
!           The computation error of the eigenvalue is zero
            Werr(1) = ZERO
            Iblock(1) = 1
            Indexw(1) = 1
         ENDIF
         RETURN
      ENDIF

!     NB is the minimum vector length for vector bisection, or 0
!     if only scalar is to be done.
      nb = ILAENV(1,'DSTEBZ',' ',N,-1,-1,-1)
      IF ( nb<=1 ) nb = 0

!     Find global spectral radius
      gl = D(1)
      gu = D(1)
      DO i = 1 , N
         gl = MIN(gl,Gers(2*i-1))
         gu = MAX(gu,Gers(2*i))
      ENDDO
!     Compute global Gerschgorin bounds and spectral diameter
      tnorm = MAX(ABS(gl),ABS(gu))
      gl = gl - FUDGE*tnorm*eps*N - FUDGE*TWO*Pivmin
      gu = gu + FUDGE*tnorm*eps*N + FUDGE*TWO*Pivmin
!     [JAN/28/2009] remove the line below since SPDIAM variable not use
!     SPDIAM = GU - GL
!     Input arguments for DLAEBZ:
!     The relative tolerance.  An interval (a,b] lies within
!     "relative tolerance" if  b-a < RELTOL*max(|a|,|b|),
      rtoli = Reltol
!     Set the absolute tolerance for interval convergence to zero to force
!     interval convergence based on relative size of the interval.
!     This is dangerous because intervals might not converge when RELTOL is
!     small. But at least a very small number should be selected so that for
!     strongly graded matrices, the code can get relatively accurate
!     eigenvalues.
      atoli = FUDGE*TWO*uflow + FUDGE*TWO*Pivmin

      IF ( irange==INDRNG ) THEN

!        RANGE='I': Compute an interval containing eigenvalues
!        IL through IU. The initial interval [GL,GU] from the global
!        Gerschgorin bounds GL and GU is refined by DLAEBZ.
         itmax = INT((LOG(tnorm+Pivmin)-LOG(Pivmin))/LOG(TWO)) + 2
         Work(N+1) = gl
         Work(N+2) = gl
         Work(N+3) = gu
         Work(N+4) = gu
         Work(N+5) = gl
         Work(N+6) = gu
         Iwork(1) = -1
         Iwork(2) = -1
         Iwork(3) = N + 1
         Iwork(4) = N + 1
         Iwork(5) = Il - 1
         Iwork(6) = Iu
!
         CALL DLAEBZ(3,itmax,N,2,2,nb,atoli,rtoli,Pivmin,D,E,E2,Iwork(5)&
     &               ,Work(N+1),Work(N+5),iout,Iwork,W,Iblock,iinfo)
         IF ( iinfo/=0 ) THEN
            Info = iinfo
            RETURN
         ENDIF
!        On exit, output intervals may not be ordered by ascending negcount
         IF ( Iwork(6)==Iu ) THEN
            Wl = Work(N+1)
            wlu = Work(N+3)
            nwl = Iwork(1)
            Wu = Work(N+4)
            wul = Work(N+2)
            nwu = Iwork(4)
         ELSE
            Wl = Work(N+2)
            wlu = Work(N+4)
            nwl = Iwork(2)
            Wu = Work(N+3)
            wul = Work(N+1)
            nwu = Iwork(3)
         ENDIF
!        On exit, the interval [WL, WLU] contains a value with negcount NWL,
!        and [WUL, WU] contains a value with negcount NWU.
         IF ( nwl<0 .OR. nwl>=N .OR. nwu<1 .OR. nwu>N ) THEN
            Info = 4
            RETURN
         ENDIF

      ELSEIF ( irange==VALRNG ) THEN
         Wl = Vl
         Wu = Vu

      ELSEIF ( irange==ALLRNG ) THEN
         Wl = gl
         Wu = gu
      ENDIF

!     Find Eigenvalues -- Loop Over blocks and recompute NWL and NWU.
!     NWL accumulates the number of eigenvalues .le. WL,
!     NWU accumulates the number of eigenvalues .le. WU
      M = 0
      iend = 0
      Info = 0
      nwl = 0
      nwu = 0
!
      DO jblk = 1 , Nsplit
         ioff = iend
         ibegin = ioff + 1
         iend = Isplit(jblk)
         in = iend - ioff
!
         IF ( in==1 ) THEN
!           1x1 block
            IF ( Wl>=D(ibegin)-Pivmin ) nwl = nwl + 1
            IF ( Wu>=D(ibegin)-Pivmin ) nwu = nwu + 1
            IF ( irange==ALLRNG .OR.                                    &
     &           (Wl<D(ibegin)-Pivmin .AND. Wu>=D(ibegin)-Pivmin) ) THEN
               M = M + 1
               W(M) = D(ibegin)
               Werr(M) = ZERO
!              The gap for a single block doesn't matter for the later
!              algorithm and is assigned an arbitrary large value
               Iblock(M) = jblk
               Indexw(M) = 1
            ENDIF

!        Disabled 2x2 case because of a failure on the following matrix
!        RANGE = 'I', IL = IU = 4
!          Original Tridiagonal, d = [
!           -0.150102010615740E+00
!           -0.849897989384260E+00
!           -0.128208148052635E-15
!            0.128257718286320E-15
!          ];
!          e = [
!           -0.357171383266986E+00
!           -0.180411241501588E-15
!           -0.175152352710251E-15
!          ];
!
!         ELSE IF( IN.EQ.2 ) THEN
!*           2x2 block
!            DISC = SQRT( (HALF*(D(IBEGIN)-D(IEND)))**2 + E(IBEGIN)**2 )
!            TMP1 = HALF*(D(IBEGIN)+D(IEND))
!            L1 = TMP1 - DISC
!            IF( WL.GE. L1-PIVMIN )
!     $         NWL = NWL + 1
!            IF( WU.GE. L1-PIVMIN )
!     $         NWU = NWU + 1
!            IF( IRANGE.EQ.ALLRNG .OR. ( WL.LT.L1-PIVMIN .AND. WU.GE.
!     $          L1-PIVMIN ) ) THEN
!               M = M + 1
!               W( M ) = L1
!*              The uncertainty of eigenvalues of a 2x2 matrix is very small
!               WERR( M ) = EPS * ABS( W( M ) ) * TWO
!               IBLOCK( M ) = JBLK
!               INDEXW( M ) = 1
!            ENDIF
!            L2 = TMP1 + DISC
!            IF( WL.GE. L2-PIVMIN )
!     $         NWL = NWL + 1
!            IF( WU.GE. L2-PIVMIN )
!     $         NWU = NWU + 1
!            IF( IRANGE.EQ.ALLRNG .OR. ( WL.LT.L2-PIVMIN .AND. WU.GE.
!     $          L2-PIVMIN ) ) THEN
!               M = M + 1
!               W( M ) = L2
!*              The uncertainty of eigenvalues of a 2x2 matrix is very small
!               WERR( M ) = EPS * ABS( W( M ) ) * TWO
!               IBLOCK( M ) = JBLK
!               INDEXW( M ) = 2
!            ENDIF
         ELSE
!           General Case - block of size IN >= 2
!           Compute local Gerschgorin interval and use it as the initial
!           interval for DLAEBZ
            gu = D(ibegin)
            gl = D(ibegin)
            tmp1 = ZERO

            DO j = ibegin , iend
               gl = MIN(gl,Gers(2*j-1))
               gu = MAX(gu,Gers(2*j))
            ENDDO
!           [JAN/28/2009]
!           change SPDIAM by TNORM in lines 2 and 3 thereafter
!           line 1: remove computation of SPDIAM (not useful anymore)
!           SPDIAM = GU - GL
!           GL = GL - FUDGE*SPDIAM*EPS*IN - FUDGE*PIVMIN
!           GU = GU + FUDGE*SPDIAM*EPS*IN + FUDGE*PIVMIN
            gl = gl - FUDGE*tnorm*eps*in - FUDGE*Pivmin
            gu = gu + FUDGE*tnorm*eps*in + FUDGE*Pivmin
!
            IF ( irange>1 ) THEN
               IF ( gu<Wl ) THEN
!                 the local block contains none of the wanted eigenvalues
                  nwl = nwl + in
                  nwu = nwu + in
                  CYCLE
               ENDIF
!              refine search interval if possible, only range (WL,WU] matters
               gl = MAX(gl,Wl)
               gu = MIN(gu,Wu)
               IF ( gl>=gu ) CYCLE
            ENDIF

!           Find negcount of initial interval boundaries GL and GU
            Work(N+1) = gl
            Work(N+in+1) = gu
            CALL DLAEBZ(1,0,in,in,1,nb,atoli,rtoli,Pivmin,D(ibegin),    &
     &                  E(ibegin),E2(ibegin),idumma,Work(N+1),          &
     &                  Work(N+2*in+1),im,Iwork,W(M+1),Iblock(M+1),     &
     &                  iinfo)
            IF ( iinfo/=0 ) THEN
               Info = iinfo
               RETURN
            ENDIF
!
            nwl = nwl + Iwork(1)
            nwu = nwu + Iwork(in+1)
            iwoff = M - Iwork(1)

!           Compute Eigenvalues
            itmax = INT((LOG(gu-gl+Pivmin)-LOG(Pivmin))/LOG(TWO)) + 2
            CALL DLAEBZ(2,itmax,in,in,1,nb,atoli,rtoli,Pivmin,D(ibegin),&
     &                  E(ibegin),E2(ibegin),idumma,Work(N+1),          &
     &                  Work(N+2*in+1),iout,Iwork,W(M+1),Iblock(M+1),   &
     &                  iinfo)
            IF ( iinfo/=0 ) THEN
               Info = iinfo
               RETURN
            ENDIF
!
!           Copy eigenvalues into W and IBLOCK
!           Use -JBLK for block number for unconverged eigenvalues.
!           Loop over the number of output intervals from DLAEBZ
            DO j = 1 , iout
!              eigenvalue approximation is middle point of interval
               tmp1 = HALF*(Work(j+N)+Work(j+in+N))
!              semi length of error interval
               tmp2 = HALF*ABS(Work(j+N)-Work(j+in+N))
               IF ( j>iout-iinfo ) THEN
!                 Flag non-convergence.
                  ncnvrg = .TRUE.
                  ib = -jblk
               ELSE
                  ib = jblk
               ENDIF
               DO je = Iwork(j) + 1 + iwoff , Iwork(j+in) + iwoff
                  W(je) = tmp1
                  Werr(je) = tmp2
                  Indexw(je) = je - iwoff
                  Iblock(je) = ib
               ENDDO
            ENDDO
!
            M = M + im
         ENDIF
      ENDDO

!     If RANGE='I', then (WL,WU) contains eigenvalues NWL+1,...,NWU
!     If NWL+1 < IL or NWU > IU, discard extra eigenvalues.
      IF ( irange==INDRNG ) THEN
         idiscl = Il - 1 - nwl
         idiscu = nwu - Iu
!
         IF ( idiscl>0 ) THEN
            im = 0
            DO je = 1 , M
!              Remove some of the smallest eigenvalues from the left so that
!              at the end IDISCL =0. Move all eigenvalues up to the left.
               IF ( W(je)<=wlu .AND. idiscl>0 ) THEN
                  idiscl = idiscl - 1
               ELSE
                  im = im + 1
                  W(im) = W(je)
                  Werr(im) = Werr(je)
                  Indexw(im) = Indexw(je)
                  Iblock(im) = Iblock(je)
               ENDIF
            ENDDO
            M = im
         ENDIF
         IF ( idiscu>0 ) THEN
!           Remove some of the largest eigenvalues from the right so that
!           at the end IDISCU =0. Move all eigenvalues up to the left.
            im = M + 1
            DO je = M , 1 , -1
               IF ( W(je)>=wul .AND. idiscu>0 ) THEN
                  idiscu = idiscu - 1
               ELSE
                  im = im - 1
                  W(im) = W(je)
                  Werr(im) = Werr(je)
                  Indexw(im) = Indexw(je)
                  Iblock(im) = Iblock(je)
               ENDIF
            ENDDO
            jee = 0
            DO je = im , M
               jee = jee + 1
               W(jee) = W(je)
               Werr(jee) = Werr(je)
               Indexw(jee) = Indexw(je)
               Iblock(jee) = Iblock(je)
            ENDDO
            M = M - im + 1
         ENDIF

         IF ( idiscl>0 .OR. idiscu>0 ) THEN
!           Code to deal with effects of bad arithmetic. (If N(w) is
!           monotone non-decreasing, this should never happen.)
!           Some low eigenvalues to be discarded are not in (WL,WLU],
!           or high eigenvalues to be discarded are not in (WUL,WU]
!           so just kill off the smallest IDISCL/largest IDISCU
!           eigenvalues, by marking the corresponding IBLOCK = 0
            IF ( idiscl>0 ) THEN
               wkill = Wu
               DO jdisc = 1 , idiscl
                  iw = 0
                  DO je = 1 , M
                     IF ( Iblock(je)/=0 .AND. (W(je)<wkill .OR. iw==0) )&
     &                    THEN
                        iw = je
                        wkill = W(je)
                     ENDIF
                  ENDDO
                  Iblock(iw) = 0
               ENDDO
            ENDIF
            IF ( idiscu>0 ) THEN
               wkill = Wl
               DO jdisc = 1 , idiscu
                  iw = 0
                  DO je = 1 , M
                     IF ( Iblock(je)/=0 .AND. (W(je)>=wkill .OR. iw==0) &
     &                    ) THEN
                        iw = je
                        wkill = W(je)
                     ENDIF
                  ENDDO
                  Iblock(iw) = 0
               ENDDO
            ENDIF
!           Now erase all eigenvalues with IBLOCK set to zero
            im = 0
            DO je = 1 , M
               IF ( Iblock(je)/=0 ) THEN
                  im = im + 1
                  W(im) = W(je)
                  Werr(im) = Werr(je)
                  Indexw(im) = Indexw(je)
                  Iblock(im) = Iblock(je)
               ENDIF
            ENDDO
            M = im
         ENDIF
         IF ( idiscl<0 .OR. idiscu<0 ) toofew = .TRUE.
      ENDIF
!
      IF ( (irange==ALLRNG .AND. M/=N) .OR.                             &
     &     (irange==INDRNG .AND. M/=Iu-Il+1) ) toofew = .TRUE.

!     If ORDER='B', do nothing the eigenvalues are already sorted by
!        block.
!     If ORDER='E', sort the eigenvalues from smallest to largest

      IF ( LSAME(Order,'E') .AND. Nsplit>1 ) THEN
         DO je = 1 , M - 1
            ie = 0
            tmp1 = W(je)
            DO j = je + 1 , M
               IF ( W(j)<tmp1 ) THEN
                  ie = j
                  tmp1 = W(j)
               ENDIF
            ENDDO
            IF ( ie/=0 ) THEN
               tmp2 = Werr(ie)
               itmp1 = Iblock(ie)
               itmp2 = Indexw(ie)
               W(ie) = W(je)
               Werr(ie) = Werr(je)
               Iblock(ie) = Iblock(je)
               Indexw(ie) = Indexw(je)
               W(je) = tmp1
               Werr(je) = tmp2
               Iblock(je) = itmp1
               Indexw(je) = itmp2
            ENDIF
         ENDDO
      ENDIF
!
      Info = 0
      IF ( ncnvrg ) Info = Info + 1
      IF ( toofew ) Info = Info + 2
!
!     End of DLARRD
!
      END SUBROUTINE DLARRD


text:
c     compute zeeman components and strengths from the g-values
c     of the lower and upper levels provided directly by the
c     user. This routine is useful for calculating the Zeeman
c     pattern of the Fe I line at 15652 A.

c     The total angular momentum (tam) of the two levels must
c     be provided in the atomic parameter file. The other values
c     specified in the transition will not be used here.

      SUBROUTINE ZEEMAN_jk(MC,MULT,DESIGN,TAM,JI,JF,DL0,NP,NL,NR,DLP,DLL,
     *DLR,SP,SL,SR,xg1,xg2)
      CHARACTER DESIGN*1
      character*1 CODE(13)/'S','P','D','F','G','H','I','K','L','M','N',
     *'O','Q'/
      dimension G(2),MULT(*),DESIGN(*),TAM(*),JI(*),JF(*),DLP(*),DLL(*),
     *DLR(*),SP(*),SL(*),SR(*)
        character*1 code1(21)/'p','1','f','2','h','3','k','4','m','5','o'
     & ,'6','r','7','t','8','u','9','v','0','w'/

      OAM=-1
      JI1=JI(1)
      JI2=JI(2)
      JF1=JF(1)
      JF2=JF(2)
      DO 58 I=1,MC
         DLL(I)=0.
         SL(I)=0.
         DLR(I)=0.
58       SR(I)=0.

      G(1)=xg1
      G(2)=xg2
      LEVEL=0

      IF((JI1+JF1).EQ.0) LEVEL=1
      IF((JI2+JF2).EQ.0) LEVEL=-1
      IF(LEVEL) 3,4,5
3     I=1                    ! TRIPLETS WITH G=0
      GO TO 6
5     I=2
6     TAM(I)=1.              ! TOTAL ANGULAR MOMENTUM
      SPIN=0.5*FLOAT(MULT(I)-1)

      DO 7 J=1,13
         IF(DESIGN(I).NE.CODE(J)) GO TO 7
         OAM=J-1                ! ORBITAL ANGULAR MOMENTUM
         GO TO 8
7     CONTINUE
      do 70 j=1,21
      if(design(I).ne.code1(j)) goto 70
      oam=float(j)/2.                ! ORBITAL ANGULAR MOMENTUM (semientero)
      go to 8
70    continue

      print*,' '
      print*,'STOP: Check the transitions in the file containing the atomic parameters.'
      print*,'Some letter used to code the orbital angular momentum is erroneous or missing.'
      print*,' '
      print*,'__________________________________________________________________________________'
      stop

8     G(I)=g(i)                  !1.5+(SPIN*(1.+SPIN)-OAM*(1.+OAM))/4.
13    NP=1                   ! TRIPLETS WITH G=0 OR G1=G2
      NL=NP
      NR=NP
      SP(1)=1.
      SL(1)=1.
      SR(1)=1.
      DLP(1)=0.
      DLL(1)=G(I)
      DLR(1)=-DLL(1)
      GO TO 9
4     DO 10 I=1,2
       SPIN=0.5*FLOAT(MULT(I)-1)

      DO 11 J=1,13
         IF(DESIGN(I).NE.CODE(J)) GO TO 11
         OAM=J-1
         GO TO 10
11    CONTINUE
      do 110 j=1,21
      if(design(I).ne.code1(j)) goto 110
      oam=float(j)/2.                ! ORBITAL ANGULAR MOMENTUM (semientero)
      go to 10
110   continue

C	LAMDE FACTOR FOR EACH LABEL
10    G(I)=g(i)   !1.5+(SPIN*(1.+SPIN)-OAM*(1.+OAM))/(2.*TAM(I)*(1.+TAM(I)))
      IF(ABS(G(2)-G(1)).GT.5.E-6) GO TO 12
      I=2
      GO TO 13
12    LEVEL=JI2-JI1
      IF(JF1.EQ.5) GO TO 14
      IF(LEVEL) 15,16,17     ! INTEGRAL J'S
15    NP=2*JI2+1
19    NL=NP
      NR=NP
      IF(NP.LE.MC) GO TO 18
      STOP 'EXIT ZEEMAN'
16    NP=2*JI2
      GO TO 19
17    NP=2*JI2-1
      GO TO 19
18    MUMIN=-JI2
      IF(JI1.LT.JI2) MUMIN=-JI1
      MUMAX=-MUMIN
      I=0
      DO 20 MU=MUMIN,MUMAX   ! PI COMPONENTS
      IF(MU.EQ.0.AND.LEVEL.EQ.0) GO TO 20
      I=I+1
      DLP(I)=FLOAT(MU)*(G(1)-G(2))
      J=MU**2
      IF(LEVEL) 21,22,23
21    SP(I)=2*(JI1**2-J)
      GO TO 20
22    SP(I)=2*J
      GO TO 20
23    SP(I)=2*(JI2**2-J)
20    CONTINUE
      MUMIN=1-JI1
      MUMAX=JI2
      I=0
      DO 24 MU=MUMIN,MUMAX   ! R-SIGMA COMPONENTS
      I=I+1
      DLR(I)=FLOAT(MU)*(G(1)-G(2))-G(1)
      IF(LEVEL) 25,26,27
25    SR(I)=(JI1-MU)*(JI2-MU+2)
      GO TO 24
26    SR(I)=(JI2+MU)*(JI2-MU+1)
      GO TO 24
27    SR(I)=(JI2+MU)*(JI1+MU)
24    CONTINUE
      MUMIN=-JI2
      MUMAX=JI1-1
      I=0
      DO 28 MU=MUMIN,MUMAX   ! L-SIGMA COMPONENTS
      I=I+1
      DLL(I)=FLOAT(MU)*(G(1)-G(2))+G(1)
      IF(LEVEL) 29,30,31
29    SL(I)=(JI1+MU)*(JI2+MU+2)
      GO TO 28
30    SL(I)=(JI2-MU)*(JI2+MU+1)
      GO TO 28
31    SL(I)=(JI2-MU)*(JI1-MU)
28    CONTINUE
      GO TO 57
14    I=2*JI1+1              ! HALF-INTEGRAL J'S
      IF(LEVEL) 32,33,34
32    NP=I-1
      NL=NP
      NR=NP
36    IF(NP.LE.MC) GO TO 35
      STOP 'EXIT ZEEMAN'
33    NP=I+1
      NL=I
      NR=I
      GO TO 36
34    NP=I+1
      NL=NP
      NR=NP
      GO TO 36
35    MUMIN=-JI2
      IF(JI1.LT.JI2) MUMIN=-JI1
      MUMAX=1-MUMIN
      I=0
      DO 37 MU=MUMIN,MUMAX   ! PI COMPONENTS
      I=I+1
      SPIN=FLOAT(MU)-0.5
      DLP(I)=(G(1)-G(2))*SPIN
      SPIN=SPIN**2
      IF(LEVEL) 38,39,40
38    SP(I)=2.*((FLOAT(JI1)+0.5)**2-SPIN)
      GO TO 37
39    SP(I)=2.*SPIN
      GO TO 37
40    SP(I)=2.*((FLOAT(JI2)+0.5)**2-SPIN)
37    CONTINUE
      MUMIN=-JI1
      MUMAX=JI2
      I=0
      DO 41 MU=MUMIN,MUMAX   ! R-SIGMA COMPONENTS
      I=I+1
      DLR(I)=(FLOAT(MU)+0.5)*(G(1)-G(2))-G(1)
      IF(LEVEL) 42,43,44
42    SR(I)=(JI1-MU)*(JI2-MU+2)
      GO TO 41
43    SR(I)=(JI2+MU+1)*(JI2-MU+1)
      GO TO 41
44    SR(I)=(JI2+MU+1)*(JI2+MU)
41    CONTINUE
      MUMIN=-MUMAX
      MUMAX=JI1
      I=0
      DO 45 MU=MUMIN,MUMAX   ! L-SIGMA COMPONENTS
      I=I+1
      DLL(I)=(FLOAT(MU)-0.5)*(G(1)-G(2))+G(1)
      IF(LEVEL) 46,47,48
46    SL(I)=(JI1+MU)*(JI2+MU+2)
      GO TO 45
47    SL(I)=(JI2-MU+1)*(JI2+MU+1)
      GO TO 45
48    SL(I)=(JI2-MU+1)*(JI1-MU+1)
45    CONTINUE
57    SUM=0.
      DO 49 I=1,NP
49    SUM=SUM+SP(I)
      DO 50 I=1,NP
50    SP(I)=SP(I)/SUM
      SPIN=0.
      SUM=0.
      DO 51 I=1,NL           ! NL=NR ASSUMED.
      SPIN=SPIN+SL(I)
51    SUM=SUM+SR(I)
      DO 52 I=1,NL
      SL(I)=SL(I)/SPIN
52    SR(I)=SR(I)/SUM
   9  CONTINUE
c      IF(MC.LE.18) GO TO 59 ! SKIP PRINTING IF SYNTHESIS
c      WRITE(6,53) (MULT(I),DESIGN(I),JI(I),JF(I),I=1,2)
c  53  FORMAT(3X,'Transition:',I3,A1,I2,'.',I1,' ----',I3,A1,I2,'.',I1/
c     *3X,'Zeeman Pattern (Lorentz units) and normalized Intensities:'/
c     *11X,'Pi',17X,'Sigma+',15X,'Sigma-')
 59   DO 54 I=1,NP
c      IF(MC.LE.18) GO TO 60
c      WRITE(6,55) DLP(I),SP(I),DLL(I),SL(I),DLR(I),SR(I)
  55  FORMAT(2X,3(F10.6,F11.8))
  60  DLP(I)=DLP(I)*DL0
      DLL(I)=DLL(I)*DL0
54    DLR(I)=DLR(I)*DL0
      IF(MC.LE.18) RETURN
c      WRITE(6,56) G
  56  FORMAT(3X,'Lande Factors: g(lower)=',F11.8,', g(upper)=',F11.8)
      RETURN
      END