text:
      SUBROUTINE STNDRD(LUNIT,MSGIN,LMSGOT,MSGOT)

C$$$  SUBPROGRAM DOCUMENTATION BLOCK
C
C SUBPROGRAM:    STNDRD
C   PRGMMR: ATOR             ORG: NP12       DATE: 2004-08-18
C
C ABSTRACT: THIS SUBROUTINE READS AN INPUT NCEP BUFR MESSAGE CONTAINED
C   WITHIN ARRAY MSGIN AND, USING THE BUFR TABLES INFORMATION ASSOCIATED
C   WITH LOGICAL UNIT LUNIT, OUTPUTS A "STANDARDIZED" VERSION OF THIS
C   SAME MESSAGE WITHIN ARRAY MSGOT.  THIS "STANDARDIZATION" INVOLVES
C   REMOVING ALL OCCURRENCES OF NCEP BUFRLIB-SPECIFIC BYTE COUNTERS AND
C   BIT PADS IN SECTION 4 AS WELL AS REPLACING THE TOP-LEVEL TABLE A FXY
C   NUMBER IN SECTION 3 WITH AN EQUIVALENT SEQUENCE OF LOWER-LEVEL
C   TABLE B, TABLE C, TABLE D AND/OR REPLICATION FXY NUMBERS WHICH
C   DIRECTLY CONSTITUTE THAT TABLE A FXY NUMBER AND WHICH THEMSELVES ARE
C   ALL WMO-STANDARD.  THE RESULT IS THAT THE OUTPUT MESSAGE IN MSGOT IS
C   NOW ENTIRELY COMPLIANT WITH WMO FM-94 BUFR REGULATIONS (I.E. IT IS
C   NOW "STANDARD"). IT IS IMPORTANT TO NOTE THAT THE SEQUENCE EXPANSION
C   WITHIN SECTION 3 MAY CAUSE THE FINAL "STANDARDIZED" BUFR MESSAGE TO
C   BE LONGER THAN THE ORIGINAL INPUT NCEP BUFR MESSAGE BY AS MANY AS
C   1000 BYTES, SO THE USER MUST ALLOW FOR ENOUGH SPACE TO ACCOMODATE
C   SUCH AN EXPANSION WITHIN THE MSGOT ARRAY.
C
C PROGRAM HISTORY LOG:
C 2004-08-18  J. ATOR    -- ORIGINAL AUTHOR
C                           THIS SUBROUTINE IS MODELED AFTER SUBROUTINE
C                           STANDARD; HOWEVER, IT USES SUBROUTINE RESTD
C                           TO EXPAND SECTION 3 AS MANY LEVELS AS
C                           NECESSARY IN ORDER TO ATTAIN TRUE WMO
C                           STANDARDIZATION (WHEREAS STANDARD ONLY
C                           EXPANDED THE TOP-LEVEL TABLE A FXY NUMBER
C                           ONE LEVEL DEEP), AND IT ALSO CONTAINS AN
C                           EXTRA INPUT ARGUMENT LMSGOT WHICH PREVENTS
C                           OVERFLOW OF THE MSGOT ARRAY
C 2005-11-29  J. ATOR    -- USE GETLENS AND IUPBS01; ENSURE THAT BYTE 4
C                           OF SECTION 4 IS ZEROED OUT IN MSGOT; CHECK
C                           EDITION NUMBER OF BUFR MESSAGE BEFORE
C                           PADDING TO AN EVEN BYTE COUNT
C
C USAGE:    CALL STNDRD (LUNIT, MSGIN, LMSGOT, MSGOT)
C   INPUT ARGUMENT LIST:
C     LUNIT    - INTEGER: FORTRAN LOGICAL UNIT NUMBER FOR BUFR FILE
C     MSGIN    - INTEGER: *-WORD ARRAY CONTAINING BUFR MESSAGE IN NCEP
C                BUFR
C     LMSGOT   - INTEGER: DIMENSIONED SIZE (IN INTEGER WORDS) OF MSGOT;
C                USED BY THE SUBROUTINE TO ENSURE THAT IT DOES NOT
C                OVERFLOW THE MSGOT ARRAY
C
C   OUTPUT ARGUMENT LIST:
C     MSGOT    - INTEGER: *-WORD ARRAY CONTAINING INPUT BUFR MESSAGE
C                NOW IN STANDARDIZED BUFR
C
C REMARKS:
C    MSGIN AND MSGOT MUST BE SEPARATE ARRAYS.
C
C    THIS ROUTINE CALLS:        BORT     GETLENS  ISTDESC  IUPB
C                               IUPBS01  MVB      NUMTAB   PKB
C                               PKC      RESTD    STATUS   UPB
C                               UPC
C    THIS ROUTINE IS CALLED BY: MSGWRT
C                               Also called by application programs.
C
C ATTRIBUTES:
C   LANGUAGE: FORTRAN 77
C   MACHINE:  PORTABLE TO ALL PLATFORMS
C
C$$$

      INCLUDE 'bufrlib.prm'

      DIMENSION ICD(MAXNC)

      COMMON /HRDWRD/ NBYTW,NBITW,NREV,IORD(8)

      DIMENSION MSGIN(*),MSGOT(*)

      CHARACTER*128 BORT_STR
      CHARACTER*8   SUBSET
      CHARACTER*4   SEVN
      CHARACTER*1   TAB

C-----------------------------------------------------------------------
C-----------------------------------------------------------------------

C  LUNIT MUST POINT TO AN OPEN BUFR FILE
C  -------------------------------------

      CALL STATUS(LUNIT,LUN,IL,IM)
      IF(IL.EQ.0) GOTO 900

C  IDENTIFY THE SECTION LENGTHS AND ADDRESSES IN MSGIN
C  ---------------------------------------------------

      CALL GETLENS(MSGIN,5,LEN0,LEN1,LEN2,LEN3,LEN4,LEN5)

      IAD3 = LEN0+LEN1+LEN2
      IAD4 = IAD3+LEN3

      LENN = LEN0+LEN1+LEN2+LEN3+LEN4+LEN5

      LENM = IUPBS01(MSGIN,'LENM')

      IF(LENN.NE.LENM) GOTO 901

      MBIT = (LENN-4)*8
      CALL UPC(SEVN,4,MSGIN,MBIT)
      IF(SEVN.NE.'7777') GOTO 902

C  COPY SECTIONS 0 THROUGH PART OF SECTION 3 INTO MSGOT
C  ----------------------------------------------------

      MXBYTO = (LMSGOT*NBYTW) - 8

      LBYTO = IAD3+7
      IF(LBYTO.GT.MXBYTO) GOTO 905
      CALL MVB(MSGIN,1,MSGOT,1,LBYTO)

C  REWRITE NEW SECTION 3 IN A "STANDARD" FORM
C  ------------------------------------------

      NSUB = IUPB(MSGIN,IAD3+ 5,16)
      ISUB = IUPB(MSGIN,IAD3+10,16)
      IBIT = (IAD3+7)*8

      CALL NUMTAB(LUN,ISUB,SUBSET,TAB,ITAB)
      IF(ITAB.EQ.0) GOTO 903

      IF (ISTDESC(ISUB).EQ.0) THEN

C         ISUB IS A NON-STANDARD TABLE A DESCRIPTOR AND NEEDS
C         TO BE EXPANDED INTO AN EQUIVALENT STANDARD SEQUENCE

          CALL RESTD(LUN,ISUB,NCD,ICD)
      ELSE

C         ISUB IS ALREADY A STANDARD DESCRIPTOR, SO JUST COPY
C         IT "AS IS" INTO THE NEW SECTION 3 (I.E. NO EXPANSION
C         IS NECESSARY!)

          NCD = 1
          ICD(NCD) = ISUB
      ENDIF

C     USE THE EDITION NUMBER TO DETERMINE THE LENGTH OF THE
C     NEW SECTION 3

      LEN3 = 7+(NCD*2)
      IBEN = IUPBS01(MSGIN,'BEN')
      IF(IBEN.LT.4) THEN
          LEN3 = LEN3+1
      ENDIF
      LBYTO = LBYTO + LEN3 - 7
      IF(LBYTO.GT.MXBYTO) GOTO 905

C     STORE THE DESCRIPTORS INTO THE NEW SECTION 3

      DO N=1,NCD
          CALL PKB(ICD(N),16,MSGOT,IBIT)
      ENDDO

C     DEPENDING ON THE EDITION NUMBER, PAD OUT THE NEW SECTION 3 WITH AN
C     ADDITIONAL ZEROED-OUT BYTE IN ORDER TO ENSURE AN EVEN BYTE COUNT

      IF(IBEN.LT.4) THEN
          CALL PKB(0,8,MSGOT,IBIT)
      ENDIF

C     STORE THE LENGTH OF THE NEW SECTION 3

      IBIT = IAD3*8
      CALL PKB(LEN3,24,MSGOT,IBIT)

C  NOW THE TRICKY PART - NEW SECTION 4
C  -----------------------------------

      NAD4 = IAD3+LEN3

      IBIT = (IAD4+4)*8
      JBIT = (NAD4+4)*8

      LBYTO = LBYTO + 4

C     COPY THE SUBSETS, MINUS THE BYTE COUNTERS AND BIT PADS,
C     INTO THE NEW SECTION 4

      DO 10 I=1,NSUB
      CALL UPB(LSUB,16,MSGIN,IBIT)

      DO L=1,LSUB-2
      CALL UPB(NVAL,8,MSGIN,IBIT)
      LBYTO = LBYTO + 1
      IF(LBYTO.GT.MXBYTO) GOTO 905
      CALL PKB(NVAL,8,MSGOT,JBIT)
      ENDDO

      DO K=1,8
      KBIT = IBIT-K-8
      CALL UPB(KVAL,8,MSGIN,KBIT)
      IF(KVAL.EQ.K) THEN
         JBIT = JBIT-K-8
         GOTO 10
      ENDIF
      ENDDO
      GOTO 904

10    ENDDO

C     FROM THIS POINT ON, WE WILL NEED (AT MOST) 6 MORE BYTES OF SPACE
C     WITHIN MSGOT IN ORDER TO BE ABLE TO STORE THE ENTIRE STANDARDIZED
C     MESSAGE (I.E. WE WILL NEED (AT MOST) 2 MORE ZEROED-OUT BYTES IN
C     SECTION 4 PLUS THE 4 BYTES '7777' IN SECTION 5), SO DO A FINAL
C     MSGOT OVERFLOW CHECK NOW.

      IF(LBYTO+6.GT.MXBYTO) GOTO 905

C     PAD THE NEW SECTION 4 WITH ZEROES UP TO THE NEXT WHOLE BYTE
C     BOUNDARY.

      DO WHILE(.NOT.(MOD(JBIT,8).EQ.0))
         CALL PKB(0,1,MSGOT,JBIT)
      ENDDO

C     DEPENDING ON THE EDITION NUMBER, WE MAY NEED TO FURTHER PAD THE
C     NEW SECTION 4 WITH AN ADDITIONAL ZEROED-OUT BYTE IN ORDER TO
C     ENSURE THAT THE PADDING IS UP TO AN EVEN BYTE BOUNDARY.

      IF( (IBEN.LT.4) .AND. (MOD(JBIT/8,2).NE.0) ) THEN
         CALL PKB(0,8,MSGOT,JBIT)
      ENDIF

      IBIT = NAD4*8
      LEN4 = JBIT/8 - NAD4
      CALL PKB(LEN4,24,MSGOT,IBIT)
      CALL PKB(0,8,MSGOT,IBIT)

C  FINISH THE NEW MESSAGE WITH AN UPDATED SECTION 0 BYTE COUNT
C  -----------------------------------------------------------

      IBIT = 32
      LENN = LEN0+LEN1+LEN2+LEN3+LEN4+LEN5
      CALL PKB(LENN,24,MSGOT,IBIT)

      CALL PKC('7777',4,MSGOT,JBIT)

C  EXITS
C  -----

      RETURN
900   CALL BORT('BUFRLIB: STNDRD - BUFR FILE IS CLOSED, IT MUST BE'//
     . ' OPEN')
901   WRITE(BORT_STR,'("BUFRLIB: STNDRD - INPUT MESSAGE LENGTH FROM'//
     . ' SECTION 0",I6," DOES NOT EQUAL SUM OF ALL INDIVIDUAL SECTION'//
     . ' LENGTHS (",I6,")")') LENM,LENN
      CALL BORT(BORT_STR)
902   WRITE(BORT_STR,'("BUFRLIB: STNDRD - INPUT MESSAGE DOES NOT '//
     . 'END WITH ""7777"" (ENDS WITH ",A)') SEVN
      CALL BORT(BORT_STR)
903   CALL BORT('BUFRLIB: STNDRD - INPUT MESSAGE TABLE D SUBSET '//
     . 'DESCRIPTOR NOT FOUND IN INTERNAL TABLE D')
904   CALL BORT('BUFRLIB: STNDRD - BIT MISMATCH COPYING SECTION 4 '//
     . 'FROM INPUT TO OUTPUT (STANDARD) MESSAGE')
905   CALL BORT('BUFRLIB: STNDRD - OVERFLOW OF OUTPUT (STANDARD) '//
     . 'MESSAGE ARRAY; TRY A LARGER DIMENSION FOR THIS ARRAY')
      END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


text:
      program t
      TYPE PERSON
      INTEGER AGE
      END TYPE PERSON
      print *,kind(AGE)
      end program t

ext: f

max_stars_repo_name: maddenp/ppp

max_issues_repo_name: maddenp/ppp

max_forks_repo_name: maddenp/ppp


text: #redirect University Pines Apartments

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
c { dg-do run }
c { dg-options "-std=gnu" }
c
* g77 0.5.23 and previous had bugs involving too little space
* allocated for EQUIVALENCE and COMMON areas needing initial
* padding to meet alignment requirements of the system.

      call subr
      end

      subroutine subr
      implicit none
      save

      character c1(11), c2(11), c3(11)
      real r1, r2, r3
      character c4, c5, c6
      equivalence (c1(2), r1)
      equivalence (c2(2), r2)
      equivalence (c3(2), r3)

      c1(1) = '1'
      r1 = 1.
      c1(11) = '1'
      c4 = '4'
      c2(1) = '2'
      r2 = 2.
      c2(11) = '2'
      c5 = '5'
      c3(1) = '3'
      r3 = 3.
      c3(11) = '3'
      c6 = '6'

      call x (c1, r1, c2, r2, c3, r3, c4, c5, c6)

      end

      subroutine x (c1, r1, c2, r2, c3, r3, c4, c5, c6)
      implicit none

      character c1(11), c2(11), c3(11)
      real r1, r2, r3
      character c4, c5, c6

      if (c1(1) .ne. '1') STOP 1
      if (r1 .ne. 1.) STOP 2
      if (c1(11) .ne. '1') STOP 3
      if (c4 .ne. '4') STOP 4
      if (c2(1) .ne. '2') STOP 5
      if (r2 .ne. 2.) STOP 6
      if (c2(11) .ne. '2') STOP 7
      if (c5 .ne. '5') STOP 8
      if (c3(1) .ne. '3') STOP 9
      if (r3 .ne. 3.) STOP 10
      if (c3(11) .ne. '3') STOP 11
      if (c6 .ne. '6') STOP 12

      end


!> This is the equilibrium master subroutine that
!! uses input to determine which type of equilibrium
!! reduced density matrix calculation to do.
!! Only currently calculates equilibrium site populations.
SUBROUTINE equilibrium_master_mpi
    USE kinds
    USE input_output
    USE random_numbers
    USE harmonic_bath
    USE hamiltonians
    USE mpi_variables
    IMPLICIT NONE

    CALL initialize_rn(mype)
    CALL initialize_bath
    CALL initialize_hel

    CALL calculate_equilibrium_site_mpi

    CALL finalize_hel
    CALL finalize_bath

END SUBROUTINE equilibrium_master_mpi


text:
*colib7.f
*colib part=7 of 9.  miscellaneous utilities.
*version=4.1 last modified: 24-apr-92
c
c  RCS $Revision: 1.3 $  $Date: 2001/09/07 02:25:12 $
c
c  see colib1.f for version history info.
c
c  this file is used by cmdc.
c
*deck bummer
*deck ibummr
      subroutine bummer( text, ierr, errtyp )
c
c  process a program error.
c  input:
c  text  = character string to be printed.
c  ierr  = internal program error to be printed.
c  errtyp = 0 for warning.  traceback may be generated. execution
c             continues.
c         = 1 for nonfatal error.  traceback may be generated.
c             execution is stopped. jcl condition code is set to allow
c             subsequent program steps to continue if possible.
c         = 2 for fatal error.  traceback may be generated.
c             execution is stopped. jcl condition code is set to abort
c             subsequent program steps if possible.
c
c  entry ibummr must be called prior to bummer() to set the output
c  unit and to perform any additional initialization.
c
c  version log:
c  24-apr-92 %val() blocks added for ibm rs6000. -rls
c  11-sep-91 parerr() calls added. -rjh
c  13-mar-91 posix code added. -rls
c  01-mar-89 write(stderr,*) and call exit() for unix machines (rls).
c  05-jul-88 unicos version (rls).
c  03-nov-87 ibm version (dcc).
c  10-sep-87 written by ron shepard.
c
      implicit integer(a-z)
c
      character*(*) text
      integer iunit, ierr, errtyp
c
c     # bummer error types.
      integer   wrnerr,  nfterr,  faterr
      parameter(wrnerr=0,nfterr=1,faterr=2)
c
      integer f77err
      integer nlist,        stderr
      save    nlist,        stderr
      data    nlist / 6 /,  stderr / 0 /
c
c*******************************************************************
      if ( errtyp .eq. wrnerr ) then
c*******************************************************************
c
c        # print a warning message and continue execution.
c
         write(nlist,6010) 'bummer (warning):', text, ierr
*mdc*if posix
*         write(stderr,6010) 'bummer (warning):', text, ierr
*mdc*elseif unicos
*         write(stderr,6010) 'bummer (warning):', text, ierr
*         call tracebk( nlist )
*mdc*elseif unix
         write(stderr,6010) 'bummer (warning):', text, ierr
*mdc*endif
         return
c
c*******************************************************************
      elseif ( errtyp .eq. nfterr ) then
c*******************************************************************
c
c        # print a warning message, stop execution.
c
         write(nlist,6010) 'bummer (nonfatal):', text, ierr
*mdc*if posix
*         write(stderr,6010) 'bummer (nonfatal):', text, ierr
*         call parerr( ierr )
*         call f77exit( 0 )
*mdc*elseif unicos
*         write(stderr,6010) 'bummer (nonfatal):', text, ierr
*         call tracebk( nlist )
*         call parerr( ierr )
*         stop 'program error'
*mdc*elseif rs6000
*c        # exit() requires call-by-value argument.
*         write(stderr,6010) 'bummer (nonfatal):', text, ierr
*         call parerr( ierr )
*         call exit( %val(0) )
*mdc*elseif unix
         write(stderr,6010) 'bummer (nonfatal):', text, ierr
         call parerr( ierr )
         call exit( 0 )
*mdc*elseif vax
*         call sys$exit(%val(42))
*mdc*elseif ibm
*         stop 901
*mdc*else
*         call parerr( ierr )
*         stop 'program error'
*mdc*endif
c
c**********************************************************************
      elseif ( errtyp .eq. faterr ) then
c**********************************************************************
c
c        # print an error message, stop execution, and abort job
c        # sequence.
c
         write(nlist,6010) 'bummer (fatal):', text, ierr
*mdc*if posix
*         write(stderr,6010) 'bummer (fatal):', text, ierr
*         call parerr( ierr )
*         call f77exit( 1 )
*mdc*elseif unicos
*         write(stderr,6010) 'bummer (fatal):', text, ierr
*         call parerr( ierr )
*c        # abort() generates tracebacks automatically.
*         call abort
*mdc*elseif rs6000
*c        # exit() requires call-by-value argument.
*         write(stderr,6010) 'bummer (fatal):', text, ierr
*         call parerr( ierr )
*         call exit( %val(1) )
*mdc*elseif unix
         write(stderr,6010) 'bummer (fatal):', text, ierr
         call parerr( ierr )
         call exit( 1 )
*mdc*elseif vax
*         call sys$exit( %val(44) )
*mdc*elseif ibm
*         stop 902
*mdc*else
*         stop 'program error'
*mdc*endif
c
c*******************************************************************
      else
c*******************************************************************
c
c        # unknown error level.  treat as a fatal error.
c
         write(nlist,6020) 'bummer (unknown): errtyp=',
     &    errtyp, text, ierr
*mdc*if posix
*         write(stderr,6020) 'bummer (unknown): errtyp=',
*     &    errtyp, text, ierr
*         call parerr( ierr )
*         call f77exit( 1 )
*mdc*elseif unicos
*         write(stderr,6020) 'bummer (unknown): errtyp=',
*     &    errtyp, text, ierr
*         call parerr( ierr )
*c        # abort() generates tracebacks automatically.
*         call abort
*mdc*elseif rs6000
*c        # exit() requires call-by-value argument.
*         write(stderr,6020) 'bummer (unknown): errtyp=',
*     &    errtyp, text, ierr
*         call parerr( ierr )
*         call exit( %val(1) )
*mdc*elseif unix
         write(stderr,6020) 'bummer (unknown): errtyp=',
     &    errtyp, text, ierr
         call parerr( ierr )
         call exit( 1 )
*mdc*elseif vax
*         call sys$exit( %val(44) )
*mdc*elseif ibm
*         stop 903
*mdc*else
*         stop 'program error'
*mdc*endif
c*******************************************************************
      endif
c*******************************************************************
c
c     # this statement is not executed, it is included
c     # just to avoid compiler warnings. -rls
      stop 'bummer error'
c
c     # initialization...
c
      entry ibummr( iunit )
c
c     # save the listing unit for use later.
c
      nlist = iunit
c
*mdc*if posix
*c     # set stderr to the correct value.
*      call f77const( 'stderr', stderr, f77err )
*      if ( f77err .ne. 0 ) then
*         write(*,6010) 'ibummr f77const() error=', f77err
*         call f77exit( 1 )
*      endif
*mdc*endif
c
      return
6010  format(1x,a,a,i10)
6020  format(1x,a,i10,a,i10)
      end
*deck pfname
c *** this routine is incremental ***
      subroutine pfname( nnames, fname )
c
c  construct a unique filename by appending the parallel process
c  number after the stub name,
c  i.e. <fname> = <name>_<nodeid>
c
c  input:  nnames = number of filenames to convert.
c          fname(1:nnames)(:) = character filename array.
c                           entries must be nonblank and must be
c                           large enough to hold the converted name.
c  output: fname(1:nnames)(:) = updated filenames.
c
c  *** this routine should only be called when running in parallel ***
c
c  11-sep-91 written by r. j. harrison.
c
      implicit none
      integer nnames
      character*(*) fname(nnames)
c
c     # bummer error types.
      integer   wrnerr,  nfterr,  faterr
      parameter(wrnerr=0,nfterr=1,faterr=2)
c
      integer lenf, k
c
      integer  nodeid, strlen
      external nodeid, strlen
      intrinsic len
c
      do 10 k = 1, nnames
c
         lenf = strlen( fname(k) )
c
         if ( lenf .eq. 0 ) then
            call bummer('pfname: name is all blanks', 0, faterr)
         elseif ( lenf+4 .gt. len(fname(k)) ) then
            call bummer('pfname: fname too short for name_id',
     &       len(fname(k)), faterr )
         endif
c
c        # append stubname with unique integer.
c
         write(fname(k)(lenf+1:lenf+4),'(a1,i3.3)') '_', nodeid()
10    continue
c
      return
      end
*deck pbginf
c *** this routine is incremental ***
      subroutine pbginf
c     # stub routine replacing the calls to tcgmsg
      return
      end
*deck pend
c *** this routine is incremental ***
      subroutine pend
c     # stub routine replacing the calls to tcgmsg
      return
      end
*deck parerr
c *** this routine is incremental ***
      subroutine parerr(ierr)
c     # stub routine replacing the calls to tcgmsg
      return
      end
*deck nodeid
c *** this routine is incremental ***
      integer function nodeid()
c
c  dummy nodeid() function.  This function should never be called.
c  the actual function, in the tcgmsg library, should be called instead.
c  this version is included in colib to ensure that libraries have been
c  loaded in the correct order.  -rls
c
c     # bummer error types.
      integer   wrnerr,  nfterr,  faterr
      parameter(wrnerr=0,nfterr=1,faterr=2)
c
      call bummer( 'colib nodeid() called.  Make sure libraries are'
     & // ' loaded in the correct order.', 0, faterr )
      nodeid = -1
      return
      end
*deck isalnm
      logical function isalnm( char )
c
c  returns true if char is alphabetic or numeric.
c
c  13-mar-91 posix code added. -rls
c  23-nov-90 cmdc block structure used. -rls
c
      character char
c
c  the ascii and ebcdic versions exists, not to make the routine faster,
c  but to make it less case sensitive
c
*mdc*if fortran
*c     # case-dependent standard fortran version.
*c     # *** warning: case dependent code ***
*      character*(*) alnmch
*      parameter ( alnmch = '0123456789'
*     $ //                  'ABCDEFGHIJKLMNOPQRSTUVWXYZ'
*     $ //                  'abcdefghijklmnopqrstuvwxyz' )
*c
*      isalnm = ( index(alnmch,char) .ne. 0 )
*mdc*elseif posix
*      logical  lf77isalnum
*      external lf77isalnum
*      isalnm = lf77isalnum( char )
*mdc*elseif ebcdic
*c     # ebcdic version.
*c     # this code assumes that ichar() returns ebcdic values.
*      logical alnm
*      integer i
*      integer iset
*c
*      integer    sets
*      parameter( sets = 7 )
*      integer   low   ( sets )
*      integer   high  ( sets )
*      data low  / 129, 145, 162, 193, 209, 226, 240 /
*      data high / 137, 153, 169, 201, 217, 233, 249 /
*c
*      i = ichar(char)
*      alnm = .false.
*      do 100 iset = 1, sets
*         alnm = alnm .or. ( low(iset) .le. i .and. i .le. high(iset) )
* 100  continue
*      isalnm = alnm
*mdc*else
c     # ascii version.
c     # this code assumes that ichar() returns ascii values.
      integer i
      integer num0, num9, uppa, uppz, lowa, lowz
      parameter ( num0 =  48 )
      parameter ( num9 =  57 )
      parameter ( uppa =  65 )
      parameter ( uppz =  90 )
      parameter ( lowa =  97 )
      parameter ( lowz = 122 )
c
      i = ichar(char)
      isalnm = ( num0 .le. i  .and.  i .le. num9 )
     & .or.    ( uppa .le. i  .and.  i .le. uppz )
     & .or.    ( lowa .le. i  .and.  i .le. lowz )
*mdc*endif
      return
      end
*deck isalph
      logical function isalph( char )
c
c  returns true if char is alphabetic.
c
c  13-mar-91 posix code added. -rls
c  23-nov-90 cmdc block structure used. -rls
c
      character char
c
c  the ascii and ebcdic versions exists, not to make the routine faster,
c  but to make it less case sensitive
c
*mdc*if fortran
*c     # case-dependent standard fortran version.
*c     # *** warning: case dependent code ***
*      character*(*) alphch
*      parameter ( alphch = 'ABCDEFGHIJKLMNOPQRSTUVWXYZ'
*     $ //                  'abcdefghijklmnopqrstuvwxyz' )
*c
*      isalph = ( index(alphch,char) .ne. 0 )
*mdc*elseif posix
*      logical  lf77isalpha
*      external lf77isalpha
*      isalph = lf77isalpha( char )
*mdc*elseif ebcdic
*c     # ebcdic version.
*c     # this code assumes that ichar() returns ebcdic values.
*      logical alph
*      integer i
*      integer iset
*c
*      integer   sets
*      parameter ( sets = 6 )
*      integer   low   ( sets )
*      integer   high  ( sets )
*      data low  / 129, 145, 162, 193, 209, 226 /
*      data high / 137, 153, 169, 201, 217, 233 /
*c
*      i = ichar(char)
*      alph = .false.
*      do 100 iset = 1, sets
*         alph = alph .or. ( low(iset) .le. i .and. i .le. high(iset) )
* 100  continue
*      isalph = alph
*mdc*else
c     # ascii version.
c     # this code assumes that ichar() returns ascii values.
      integer i
      integer uppa, uppz, lowa, lowz
      parameter ( uppa =  65 )
      parameter ( uppz =  90 )
      parameter ( lowa =  97 )
      parameter ( lowz = 122 )
c
      i = ichar(char)
      isalph = ( uppa .le. i  .and.  i .le. uppz )
     & .or.    ( lowa .le. i  .and.  i .le. lowz )
*mdc*endif
      return
      end
*deck streq
      logical function streq( str1, str2 )
c
c  returns true if two strings are equal using
c  a case-insensitive comparison.
c
      character*(*)     str1
      character*(*)     str2
c
c  --case insensitive string comparision
c  --need scratch arrays because we don't want to upper case constants
c  --hate having fixed length temp strings, but no way to avoid
c  --could ignore leading spaces, but i think better elsewhere
c    (if someone else adds, i won't delete)
c  --the ascii version exists solely for speed
c
      integer           len1
      integer           len2
*mdc*if fortran
*      integer           mstrln
*      parameter         ( mstrln = 512 )
*      character*(mstrln) tstr1
*      character*(mstrln) tstr2
*mdc*elseif ebcdic
*      integer           mstrln
*      parameter         ( mstrln = 500 )
*      character*(mstrln) tstr1
*      character*(mstrln) tstr2
*mdc*else
c     # ascii version.
      integer           i
      integer           nstr1
      integer           nstr2
      integer           uldiff,      uppa,      uppz
      parameter       ( uldiff = 32, uppa = 65, uppz = 90 )
*mdc*endif
c
c     # bummer error types.
      integer   wrnerr,  nfterr,  faterr
      parameter(wrnerr=0,nfterr=1,faterr=2)
c
c
      len1 = len(str1)
      len2 = len(str2)
c
*mdc*if ebcdic
*c     # convert to upper case and use fortran comparison.
*c-      if ( len1 .gt. mstrln  .or.  len2 .gt. mstrln ) stop 909
*      if ( len1 .gt. mstrln )
*     & call bummer('streq: string too long, len1=',len1,faterr)
*      if ( len2 .gt. mstrln )
*     & call bummer('streq: string too long, len2=',len2,faterr)
*c
*      tstr1(1:len1) = str1
*      tstr2(1:len2) = str2
*      call strupp( tstr1(1:len1) )
*      call strupp( tstr2(1:len2) )
*c
*      streq = tstr1(1:len1) .eq. tstr2(1:len2)
*mdc*else
c
c     # ascii version
c     # this code assumes that ichar() returns ascii values.
c
      streq = .false.
c
      do 100 i = 1, min( len1, len2 )
         nstr1 = ichar(str1(i:i))
         nstr2 = ichar(str2(i:i))
c return if the characters don't "match"
c they match if they're the same or upper/lower of same character
c if they are upper/lower of same char, they differ by uldiff
c   and the smallest one is an uppercase character
         if ( .not. ( nstr1 .eq. nstr2  .or.  (
     $        abs( nstr1 - nstr2 ) .eq. uldiff .and.
     $        min(nstr1,nstr2).ge.uppa .and. min(nstr1,nstr2).le.uppz )
     $        ) ) return
 100  continue
c
      if ( len1 .eq. len2 ) then
         streq = .true.
      elseif ( len1 .lt. len2 ) then
         streq = str2(len1+1:) .eq. ' '
      else
         streq = str1(len2+1:) .eq. ' '
      endif
*mdc*endif
c
      return
      end
*deck streqs
      logical function streqs( str1, str2 )
c
c  compare short string one to the initial characters of str two
c
      character*(*)     str1
      character*(*)     str2
c
c  --the strings are "equal" if the leading characters of the (possibly)
c    longer second string equal the shorter first string
c  --originally the two strings were treated equivalently.  while this
c    feels good, i needed the current behavior
c  --trailing blanks are ignored
c
      integer           strl
c
      external          streq
      logical           streq
      external          strlen
      integer           strlen
c
c
      strl = min( strlen(str1), len(str2) )
c
      streqs = streq( str1, str2(1:strl) )
c
      return
      end
*deck strlen
      integer function strlen( string )
c
c  returns "length" of string -- position of last non-space character
c  (at times it would be convienent to return a character variable of
c  the proper length.  unfortunately, that is not possible in fortran.)
c
c  13-mar-91 posix code added. -rls
c
      character*(*)     string
c
      integer           istr
c
*mdc*if posix
*      integer  if77lentrim
*      external if77lentrim
*      strlen = if77lentrim( string )
*mdc*else
c
c     # assume string has no characters.
      strlen = 0
c
      do 100 istr = len(string), 1, -1
         if ( string(istr:istr) .ne. ' ' ) then
            strlen = istr
            return
         endif
 100  continue
*mdc*endif
c
      return
      end
*deck strlow
      subroutine strlow( string )
c
c  change all upper case characters in string to lower
c
      character*(*)     string
c
c  the ascii and ebcdic versions exists, not to make the routine faster,
c  but to make it less case sensitive
c
      integer           istr
c
*mdc*if fortran
*c     # case-dependent standard fortran version.
*c     # warning: case dependent code
*      integer          indx
*      character*26     lower
*      character*26     upperc
*      data lower  / 'abcdefghijklmnopqrstuvwxyz' /
*      data upperc / 'ABCDEFGHIJKLMNOPQRSTUVWXYZ' /
*c
*      do 100 istr = 1, len(string)
*         indx = index(upperc, string(istr:istr) )
*         if ( indx .ne. 0 )
*     $        string(istr:istr) = lower(indx:indx)
* 100  continue
*mdc*elseif ebcdic
*c     # ebcdic version.
*c     # this code assumes that ichar() returns ebcdic values.
*      integer           ich
*      integer           uldiff
*      parameter ( uldiff = -64 )
*      integer   sets
*      parameter ( sets = 3 )
*      integer   low   ( sets )
*      integer   high  ( sets )
*      data low  / 193, 209, 226 /
*      data high / 201, 217, 233 /
*c
*c
*      do 100 istr = 1, len(string)
*         ich = ichar( string(istr:istr) )
*         if ( ( low(1) .le. ich  .and.  ich .le. high(1) ) .or.
*     $        ( low(2) .le. ich  .and.  ich .le. high(2) ) .or.
*     $        ( low(3) .le. ich  .and.  ich .le. high(3) ) )
*     $        string(istr:istr) = char( ich + uldiff )
* 100  continue
*mdc*else
c     # ascii version.
c     # this code assumes that ichar() returns ascii values.
      integer           ich
      integer uppa, uppz, uldiff
      parameter ( uppa = 65 )
      parameter ( uppz = 90 )
      parameter ( uldiff = 32 )
c
      do 100 istr = 1, len(string)
         ich = ichar( string(istr:istr) )
         if ( uppa .le. ich  .and.  ich .le. uppz )
     $        string(istr:istr) = char( ich + uldiff )
 100  continue
*mdc*endif
      return
      end
*deck strnds
      integer function strnds ( number, target, list, inc )
      implicit integer  ( a-z )
c
c     function strnds returns the element number (index) in
c     array 'list' containing the value 'target' using streqs (short
c     equals) for equality check
c
      integer           number, inc
      character*(*)     list( inc, number ), target
c
      integer           index
c
      external          streqs
      logical           streqs
c
c
c     assume target will not be found
c
      strnds = 0
c
      do 100 index = 1, number
         if ( streqs( list( 1, index ), target ) ) then
c
c           target found
c
            strnds = index
            return
         endif
 100  continue
c
c     target not found
c
      return
      end
*deck strndx
      integer function strndx ( number, target, list, inc )
      implicit integer  ( a-z )
c
c     function strndx returns the element number (index) in
c     array 'list' containing the value 'target' using streq (case
c     insenstive comparison) for equality check
c
      integer           number, inc
      character*(*)     list( inc, number ), target
c
      integer           index
c
      external          streq
      logical           streq
c
c
c     # assume target will not be found.
c
      strndx = 0
c
      do 100 index = 1, number
         if ( streq( list( 1, index ), target ) ) then
c
c           # target found.
c
            strndx = index
            return
         endif
 100  continue
c
c     # target not found
c
      return
      end
*deck strskp
      subroutine strskp( line, lptr, chars )
c
c  skips characters in line starting at lptr
c  until a character not in chars is found
c
c  lptr = len(line) + 1 if all characters in line are in chars
c
c
      implicit integer(a-z)
      character*(*)     line
      integer           lptr
      character*(*)     chars
c
c
  100 continue
c
      if ( lptr .le. len(line)  .and.
     $     index( chars, line(lptr:lptr) ) .gt. 0 ) then
c
         lptr = lptr + 1
c
         goto 100
      endif
c
c
      return
      end
*deck strtil
      subroutine strtil( line, lptr, chars )
c
c  skips characters in line starting at lptr
c  until a character in chars is found
c
c  lptr = len(line) + 1 if no more characters in line are in chars
c
c
      implicit integer(a-z)
      character*(*)     line
      integer           lptr
      character*(*)     chars
c
c
  100 continue
c
      if ( lptr .le. len(line)  .and.
     $     index( chars, line(lptr:lptr) ) .eq. 0 ) then
c
         lptr = lptr + 1
c
         goto 100
      endif
c
c
      return
      end
*deck strtok
      subroutine strtok( line, lptr, token, delim, error )
c
c  parse a token delimited by characters in delim
c
      implicit integer(a-z)
      character*(*)     line
      integer           lptr
      character*(*)     token
      character*(*)     delim
      integer           error
c
      integer           sptr
c
      integer     noerr,     longtk,      linend
      parameter ( noerr = 0, longtk = -2, linend =1 )
c
      external         strskp, strtil
c
c
      error = noerr
c
      call strskp( line, lptr, delim )
      if ( lptr .gt. len(line) ) then
         error = linend
         return
      endif
c
      sptr = lptr
      call strtil( line, lptr, delim )
c
      if ( lptr - sptr .gt. len( token ) ) error = longtk
c
      token = line( sptr : min( lptr, sptr + len(token) ) - 1 )
c
      return
      end
*deck strupp
*deck allcap
      subroutine allcap( string )
c
c  change all lower case characters in string to upper
c
      character*(*)     string
c
c  the ascii and ebcdic versions exists, not to make the routine faster,
c  but to make it less case sensitive
c
      integer           istr
c
c     # strupp() and allcap() do the same thing. -rls
c
      entry strupp( string )
c
*mdc*if fortran
*c     # case-dependent standard fortran version.
*c     # warning: case dependent code
*      integer           indx
*c     # warning: case dependent code
*c     # vax quirk; upperc(:) cannot be a parameter constant (rls).
*      character*(*)     lower
*      character*26     upperc
*      parameter ( lower  = 'abcdefghijklmnopqrstuvwxyz' )
*      data        upperc / 'ABCDEFGHIJKLMNOPQRSTUVWXYZ' /
*c
*      do 100 istr = 1, len(string)
*         indx = index(lower, string(istr:istr) )
*         if ( indx .ne. 0 )
*     $        string(istr:istr) = upperc(indx:indx)
* 100  continue
*mdc*elseif ebcdic
*c     # ebcdic version.
*c     # this code assumes that ichar() returns ebcdic values.
*      integer           ich
*      integer           uldiff
*      parameter ( uldiff = 64 )
*      integer   sets
*      parameter ( sets = 3 )
*      integer   low   ( sets )
*      integer   high  ( sets )
*      data low  / 129, 145, 162 /
*      data high / 137, 153, 169 /
*c
*c
*      do 100 istr = 1, len(string)
*         ich = ichar( string(istr:istr) )
*         if ( ( low(1) .le. ich  .and.  ich .le. high(1) ) .or.
*     $        ( low(2) .le. ich  .and.  ich .le. high(2) ) .or.
*     $        ( low(3) .le. ich  .and.  ich .le. high(3) ) )
*     $        string(istr:istr) = char( ich + uldiff )
* 100  continue
*mdc*else
c     # ascii version.
c     # this code assumes that ichar() returns ascii values.
      integer           ich
      integer lowa, lowz, uldiff
      parameter ( lowa =  97 )
      parameter ( lowz = 122 )
      parameter ( uldiff = -32 )
c
      do 100 istr = 1, len(string)
         ich = ichar( string(istr:istr) )
         if ( lowa .le. ich  .and.  ich .le. lowz )
     $        string(istr:istr) = char( ich + uldiff )
 100  continue
*mdc*endif
      return
      end
*deck strwr
      subroutine strwr( unit, string )
c
c  write string to file -- strips trailing blanks.
c  note: fortran carriage control is *not* used by this routine.
c
      integer           unit
      character*(*)     string
c
      integer           istr
c
      external          strlen
      integer           strlen
c
c
      istr = strlen( string )
      if ( istr .gt. 0 ) then
         write ( unit, '(a)' ) string(1:istr)
      else
         write ( unit , '(a)' )
      endif
c
      return
      end
*deck unicas
      character*(*) function unicas ( string )
c
      character*(*) string
c
c  return string converted to site-dependent consistent case
c
c  14-dec-90 temp*150 removed. -rls
c
      unicas = string
*mdc*if osu
*      call strupp( unicas )
*mdc*else
      call strlow( unicas )
*mdc*endif
c
      return
      end
*deck ifchl
      integer function ifchl( line, chlist, match )
c
c  find the position of the first character in line(1:) that is
c  contained in chlist(1:).  match is set to the index of the
c  matching character in chlist(1:).
c
c  *** beware of side-effects with the argument match ***
c
c  written by ron shepard 1-oct-87.
c
      implicit integer(a-z)
c
      character*(*) line, chlist
      integer match
c
      integer i
c
      do 10 i = 1, len(line)
         ifchl = i
         match = index(chlist,line(i:i))
         if ( match .ne. 0 ) return
10    continue
      ifchl = 0
c
      return
      end
*deck ifnch
      integer function ifnch( line, chlist )
c
c  find the position of the first character in line(1:) that is
c  not contained in chlist(1:).
c
c  written by ron shepard 1-oct-87.
c
      implicit integer(a-z)
c
      character*(*) line, chlist
c
      integer i
c
      do 10 i = 1, len(line)
         ifnch = i
         if ( index( chlist, line(i:i) ) .eq. 0 ) return
10    continue
      ifnch = 0
c
      return
      end
*deck who2c
c *** this routine is incremental ***
      subroutine who2c( prognm, nlist )
c
c  print out a short banner containing the name and address of the
c  local programmer responsible for maintaining the calling program.
c
c  input:
c  prognm = left-justified character string containing the
c           name of the calling program.
c  nlist  = output unit number for the listing file.
c
c  14-feb-91 written by rls/eas/rmp.
c
      implicit integer(a-z)
c
      integer nlist
      character*(*) prognm
c
      logical  streq
      external streq
c
*mdc*if osu
*      integer i
*      integer    rmpnum
*      parameter( rmpnum=5 )
*      character*6 rmpprg(rmpnum)
*c     # the following codes are maintained at OSU by R. M. Pitzer:
*      data rmpprg / 'argos', 'cnvrt', 'scfpq', 'tran', 'cidbg' /
*mdc*endif
c
      write (nlist,6000) prognm
6000  format(/' This Version of Program ', a, ' is Maintained by:')
c
*mdc*if argonne
c     # alliant, sun, stardent titan, vax, cray y-mp, cray 2,
c     # intel hypercube, intel delta (in progress), ibm rs/6000
      write (nlist,6010)
6010  format(
     & t5,'Ron Shepard'/
     & t5,'Theoretical Chemistry Group'/
     & t5,'Chemistry Division'/
     & t5,'Argonne National Laboratory'/
     & t5,'Argonne, IL  60439'/
     & t5,'Internet: shepard@tcg.anl.gov'/)
c
*mdc*elseif wien
*c     # sun, ibm, convex, ibm 3090, ibm rs/6000, sgi,
*c     # intel delta (in progress)
*      write (nlist,6010)
*6010  format(
*     & t5,'Hans Lischka'/
*     & t5,'Institute for Theoretical Chemistry'/
*     & t5,'University of Vienna'/
*     & t5,'Waeringerstr 17, A-1090 Wien, Austria'/
*     & t5,'Internet: hans@itc.univie.ac.at'/)
*c
*mdc*elseif osu
*c     # sun, vax, ibm, cray y-mp, stardent gs
*      do 10 i = 1, rmpnum
*         if ( streq( prognm, rmpprg(i) ) ) then
*            write (nlist,6010)
*            return
*         endif
*10    continue
*c
*      write (nlist,6020)
*c
*6010  format(
*     & t5,'Russell M. Pitzer'/
*     & t5,'Department of Chemistry'/
*     & t5,'The Ohio State University'/
*     & t5,'Columbus, OH  43210'/
*     & t5,'Internet: pitzer@neon.mps.ohio-state.edu'/)
*6020  format(
*     & t5,'Gary Kedziora'/
*     & t5,'Department of Chemistry'/
*     & t5,'The Ohio State University'/
*     & t5,'Columbus, OH  43210'/
*     & t5,'Internet: kedziora@mps.ohio-state.edu'/)
*c
*mdc*elseif radlab
*c     # convex c-1
*      write (nlist,6010)
*6010  format(
*     & t5,'John Bentley'/
*     & t5,'Radiation Laboratory',/
*     & t5,'University of Notre Dame'/
*     & t5,'Notre Dame, IN  46556-0768'/
*     & t5,'Bitnet: bentley@ndradlab'/)
*c
*mdc*elseif anu
*c     # fujitsu vp2200
*      write (nlist,6010)
*6010  format(
*     & t5,'Ross Nobes'/
*     & t5,'ANU Supercomputer Facility',/
*     & t5,'GPO Box 4'/
*     & t5,'Canberra  ACT  2601'/
*     & t5,'Australia'/
*     & t5,'Internet: R.Nobes@anu.edu.au'/)
*c
*mdc*else
*c
*c     # this is just to keep everyone honest. (:-)  -rls
*c
*      write (nlist,6010)
*6010  format(
*     & t5,52('?')/
*     & t5,'???  Tell your local COLUMBUS programmer to add  ???',/
*     & t5,'???  his name to routine who2c() in colib*.f     ???'/
*     & t5,52('?')/)
*c
*mdc*endif
c
      return
      end

max_stars_repo_name: MOLFDIR/MOLFDIR

max_issues_repo_name: yingxingcheng/MOLFDIR

max_forks_repo_name: yingxingcheng/MOLFDIR


text:
module wrapping1
implicit none
integer, parameter :: c_int = 4

contains

subroutine sub1(a, b) bind(c)
integer(c_int), intent(in) :: a
integer(c_int), intent(out) :: b
b = a + 1
end subroutine

end module

max_stars_repo_name: Thirumalai-Shaktivel/lfortran

max_issues_repo_name: Thirumalai-Shaktivel/lfortran

max_forks_repo_name: Thirumalai-Shaktivel/lfortran


text:
subroutine sigma_function(latitude,lhour,dayear,altitude,fpsolar,kpindex,stdev)

    ! last version updated : 05/01/2021
    ! Author : Claude Boniface (CNES)

    !----------------------------------!
    ! ONE-SIGMA ERROR EVALUATION MODULE!
    !----------------------------------!

    ! The model is well justified in the range [200 km - 500 km].
    ! However, strong assumptions must be added to extend the altitude range from 100 km to 1000 km due to unavailability of data:
    ! at lower altitudes (< 200 km), a worst-case is chosen by fixing the evolution of 1-sigma to their values at the altitude of 200 km,
    ! at higher altitudes (> 500 km), an optimistic case is assumed by keeping the evolution of 1-sigma for the altitude of 500 km.

    !     latitude= latitude (in deg, -90;90)
    !     lhour   = local time (in hr, 0-24)
    !     dayear  = day of year (1-366)
    !     altitude= altitude (in km) greater than 120 km
    !     fpsolar = mean flux of last 81 days at t   =   fbar(1)
    !     kpindex = kp delayed by 3 hours = akp(1)

    !----------------------!
    ! DECLARATION VARIABLES!
    !----------------------!

    implicit none
    real :: latitude, lhour, dayear, altitude, fpsolar, kpindex, stdev ! argument parameter values
    real :: Copt_L, Copt_M, Copt_H ! projection coefficient in hyperplane equations
    real :: Res_L, Res_M, Res_H    ! direction coefficient in hyperplane equations
    real :: stdev_L, stdev_M, stdev_H  ! sigma values at different range of altitude
    ! sigma at low altitude (sigma < 200 km)
    real :: sigma_FP_L ! 1-sigma (Fp proxy)
    real :: sigma_KP_L ! 1-sigma (Kp index)
    real :: sigma_alat_L  !1-sigma in latitude
    real :: sigma_hl_L    !1-sigma in local hour
    real :: sigma_day_L   ! 1-sigma in day of year
    real :: sigma_alti_L  ! 1-sigma in altitude
    ! sigma at middle altitude (200 km <= sigma =< 500 km)
    real :: sigma_FP_M ! 1-sigma (Fp proxy)
    real :: sigma_KP_M ! 1-sigma (Kp index)
    real :: sigma_alat_M  !1-sigma in latitude
    real :: sigma_hl_M    !1-sigma in local hour
    real :: sigma_day_M   ! 1-sigma in day of year
    real :: sigma_alti_M  ! 1-sigma in altitude
    ! sigma at high altitude (sigma > 500 km)
    real :: sigma_FP_H ! 1-sigma (Fp proxy)
    real :: sigma_KP_H ! 1-sigma (Kp index)
    real :: sigma_alat_H  !1-sigma in latitude
    real :: sigma_hl_H    !1-sigma in local hour
    real :: sigma_day_H   ! 1-sigma in day of year
    real :: sigma_alti_H  ! 1-sigma in altitude

    !-----------------------!
    ! ONE-SIGMA CALCULATIONS!
    !-----------------------!

    altit : If(altitude.LT.200) then

        ! Global 1-sigma calculation by means of statistical binning of data (GOCE, CHAMP, GRACE and SWARM) and Hyperplan equation (6 dimensions)
	! Polynomial extrapolation by fixing altitude value at 200 km

        Copt_L = 0.929
        Res_L = 112.79
        sigma_alat_L = std_alat_M(latitude)
        sigma_hl_L = std_hl_M(lhour)
        sigma_day_L = std_day_M(dayear)
        sigma_alti_L = std_alti_M(200.)
        sigma_FP_L = std_FP_M(fpsolar)
        sigma_KP_L = std_KP_M(kpindex)

        stdev = sigma_alat_L + sigma_hl_L + sigma_day_L + sigma_alti_L + sigma_FP_L + sigma_KP_L + Copt_L - Res_L

    ElseIf(altitude.GE.200.AND.altitude.LE.500) then

        ! Global 1-sigma calculation by means of statistical binning data (GOCE, CHAMP, GRACE and SWARM)) and Hyperplan equation (6 dimensions)
	! Polynomial interpolation

        Copt_M = 0.929
        Res_M = 112.79
        sigma_alat_M = std_alat_M(latitude)
        sigma_hl_M = std_hl_M(lhour)
        sigma_day_M = std_day_M(dayear)
        sigma_alti_M = std_alti_M(altitude)
        sigma_FP_M = std_FP_M(fpsolar)
        sigma_KP_M = std_KP_M(kpindex)

        stdev = sigma_alat_M + sigma_hl_M + sigma_day_M + sigma_alti_M + sigma_FP_M + sigma_KP_M + Copt_M - Res_M

    ElseIf(altitude.GT.500) then

        ! Global 1-sigma calculation by means of statistical binning data (GOCE, CHAMP, GRACE and SWARM) and Hyperplan equation (6 dimensions)
	! Polynomial extrapolation by fixing altitude value at 500 km

        Copt_H = 0.929
        Res_H = 112.79
        sigma_alat_H = std_alat_M(latitude)
        sigma_hl_H = std_hl_M(lhour)
        sigma_day_H = std_day_M(dayear)
        sigma_alti_H = std_alti_M(500.)
        sigma_FP_H = std_FP_M(fpsolar)
        sigma_KP_H = std_KP_M(kpindex)

        stdev = sigma_alat_H + sigma_hl_H + sigma_day_H + sigma_alti_H + sigma_FP_H + sigma_KP_H + Copt_H - Res_H

    Endif altit

contains

    !--------------------!
    ! ANALYTICAL FUNCTONS!
    !--------------------!

    ! General polynomial function

    function std_poly(x,n,a)
        integer :: h
        real :: std_poly
        integer, intent(in) :: n
        real, intent(in) :: x
        real, dimension(0:n) :: a
        h = 0
        std_poly = 0.
        Do h = 0,n
            std_poly = std_poly + a(h)*x**(h)
        enddo
    end function

    ! Solar Flux Fp
    !--------------

    function std_FP_M(x)
        real :: std_FP_M
        real, intent(in) :: x
        integer :: d
        real, dimension(0:1) :: c
        d = 1
        c(0) = 26.04
        c(1) = -3.36e-2
        std_FP_M = std_poly(x,d,c)
    end function std_FP_M

    ! Magnetic index Kp
    !------------------

    function std_KP_M(x)
        real :: std_KP_M
        real, intent(in) :: x
        integer :: d
        real, dimension(0:3) :: c
        d = 3
        c(0) = 22.47
        c(1) = 7.46e-1
        c(2) = -1.99e-1
        c(3) = 6.56e-2
        std_KP_M = std_poly(x,d,c)
    end function std_KP_M

    ! Latitude
    !---------

    function std_alat_M(x)
        real :: std_alat_M
        real, intent(in) :: x
        integer :: d
        real, dimension(0:2) :: c
        d = 2
        c(0) = 22.33
        c(1) = -7.19e-3
        c(2) = 4.72e-4
        std_alat_M = std_poly(x,d,c)
    end function std_alat_M

    ! Local hour
    !-----------

    function std_hl_M(x)
        real :: std_hl_M
        real, intent(in) :: x
        integer :: d
        real, dimension(0:3) :: c
        d = 3
        c(0) = 28.53
        c(1) = -3.95e-1
        c(2) = -5.64e-2
        c(3) = 3.15e-3
        std_hl_M = std_poly(x,d,c)
    end function std_hl_M

    ! Altitude
    !---------

    function std_alti_M(x)
        real :: std_alti_M
        real, intent(in) :: x
        integer :: d
        real, dimension(0:1) :: c
        d = 1
        c(0) = -9.73
        c(1) = 7.92e-2
        std_alti_M = std_poly(x,d,c)
    end function std_alti_M

    ! Day of year
    !------------
    ! Low solar flux regime

    function std_day_M(x)
        real :: std_day_M
        real, intent(in) :: x
        integer :: d
        real, dimension(0:1) :: c
        d = 1
        c(0) = 23.15
        c(1) = -3.91e-3
        std_day_M = std_poly(x,d,c)
    end function std_day_M

end subroutine

max_stars_repo_name: YoucunWang/mcm

max_stars_repo_licenses: DOC Unlicense

max_issues_repo_name: YoucunWang/mcm

max_issues_repo_licenses: DOC Unlicense

max_forks_repo_name: YoucunWang/mcm

max_forks_repo_licenses: DOC Unlicense


text:
      PROGRAM OFFDIF
      USE MOD_ZA  ! HYCOM array I/O interface
      IMPLICIT NONE
C
C     WIND ARRAYS.
C
      INTEGER, ALLOCATABLE :: MSK(:,:)
      REAL*4,  ALLOCATABLE :: TXO(:,:),TYO(:,:),
     +                        TXN(:,:),TYN(:,:),
     +                        TXH(:,:),TYH(:,:)
C
      REAL*4  WDYNXT
C
C     NAMELISTS.
C
      CHARACTER*40     CTITLE
      REAL*4           SCALEN,SCALEO
      NAMELIST/WTITLE/ CTITLE
      NAMELIST/WOFSET/ SCALEN,SCALEO
C
C**********
C*
C  1) CREATE A WIND STRESS OFFSET, FROM A WIND SET'S ANNUAL CLIMATOLOGY
C      (OLD) AND THE DESIRED REPLACEMENT (NEW) ANNUAL CLIMATOLOGY.
C
C  2) REGION SPECIFICATION (FROM THE OCEAN MODEL):
C
C        IDM    = 1ST DIMENSION OF MAJOR MODEL ARRAYS
C        JDM    = 2ND DIMENSION OF MAJOR MODEL ARRAYS
C
C  3) INPUT:
C        ON UNIT  5:  NAMELISTS /WTITLE/ AND /WOFSET/
C        ON UNIT 44:  NEW .a/.b FORMAT MODEL MEAN TAUEWD
C        ON UNIT 45:  NEW .a/.b FORMAT MODEL MEAN TAUNWD
C        ON UNIT 54:  OLD .a/.b FORMAT MODEL MEAN TAUEWD
C        ON UNIT 55:  OLD .a/.b FORMAT MODEL MEAN TAUNWD
C     OUTPUT:
C        ON UNIT 64:  .a/.b FORMAT MODEL OFFSET TAUEWD
C        ON UNIT 65:  .a/.b FORMAT MODEL OFFSET TAUNWD
C
C     NAMELIST /WTITLE/:
C        CTITLE - 40 CHARACTER WIND OFFSET TITLE
C
C     NAMELIST /WOFSET/:
C        SCALEN - SCALE FACTOR FOR THE NEW ANNUAL CLIMATOLOGY, SEE (5)
C        SCALEO - SCALE FACTOR FOR THE OLD ANNUAL CLIMATOLOGY, SEE (5)
C
C  4) THE INPUT AND OUTPUT WIND STRESSES ARE ALWAYS ON THE HYCOM P-GRID
C      AND ARE IN HYCOM ARRAY FILE FORMAT.
C
C  5) 'SCALEN' AND 'SCALEO' ARE SUPPLIED TO SCALE THE DESIRED MEAN
C      SO THAT IT IS COMPATIBLE WITH THE WIND SET.  USUALLY FOR HYCOM,
C      SCALEN=SCALEO=1.0 (THE DEFAULT).  HOWEVER WHEN CONBINING
C      PARTIAL MEANS SCALEN-SCALEO=1.0 WOULD BE APPROPRIATE.
C
C     THE OFFSET PRODUCED IS:  SCALEN*NEW - SCALEO*OLD.
C
C     WHEN THIS IS ADDED TO THE WIND SET (ASSUMING SCALEO=1.0),
C      IT WILL REPLACE THE OLD MEAN WITH SCALEN*NEW.
C
C  6) ALAN J. WALLCRAFT,  NAVAL RESEARCH LABORATORY,  JANUARY 2002.
C*
C**********
C
      INTEGER I,J
      REAL*4  XMIN,XMAX,XAVE,XRMS,YMIN,YMAX,YAVE,YRMS
      REAL*4  HMINA,HMINB,HMAXA,HMAXB
C
      CHARACTER PREAMBL(5)*79,CTITLO*40,CLINE*80
C
      CALL XCSPMD  !define idm,jdm
      ALLOCATE( TXO(IDM,JDM),TYO(IDM,JDM) )
      ALLOCATE( TXN(IDM,JDM),TYN(IDM,JDM) )
      ALLOCATE( TXH(IDM,JDM),TYH(IDM,JDM) )
C
C     READ IN, AND OUTPUT, THE WIND TITLE(S).
C
      CALL ZAIOST
C
      CALL ZHOPEN(6, 'FORMATTED', 'UNKNOWN', 0)
C
      CALL ZAIOPN('OLD', 44)
      CALL ZAIOPN('OLD', 45)
      CALL ZHOPEN(44, 'FORMATTED', 'OLD', 0)
      CALL ZHOPEN(45, 'FORMATTED', 'OLD', 0)
      READ(44,4101) PREAMBL
      READ(45,4101) PREAMBL
      CTITLO = PREAMBL(1)
      WRITE(6,6000) 'NEW MEAN:',CTITLO
      CALL ZHFLSH(6)
C
      CALL ZAIOPN('OLD', 54)
      CALL ZAIOPN('OLD', 55)
      CALL ZHOPEN(54, 'FORMATTED', 'OLD', 0)
      CALL ZHOPEN(55, 'FORMATTED', 'OLD', 0)
      READ(54,4101) PREAMBL
      READ(55,4101) PREAMBL
      CTITLO = PREAMBL(1)
      WRITE(6,6000) 'OLD MEAN:',CTITLO
      CALL ZHFLSH(6)
C
      READ( 5,WTITLE)
      WRITE(6,6000) 'OFFSET:  ',CTITLE
      WRITE(6,*)
      CALL ZHFLSH(6)
C
      CALL ZAIOPN('NEW', 64)
      CALL ZAIOPN('NEW', 65)
      CALL ZHOPEN(64, 'FORMATTED', 'NEW', 0)
      CALL ZHOPEN(65, 'FORMATTED', 'NEW', 0)
      PREAMBL(1) = CTITLE
      WRITE(64,4101) PREAMBL
      WRITE(65,4101) PREAMBL
      WRITE(6,*)
      WRITE(6, 4101) PREAMBL
      WRITE(6,*)
      CALL ZHFLSH(6)
C
C     NAMELIST /WOFSET/.
C
      SCALEO = 1.0
      SCALEN = 1.0
      READ( 5,WOFSET)
      WRITE(6,WOFSET)
      WRITE(6,*)
      CALL ZHFLSH(6)
C
C     READ IN THE NEW MEAN.
C
      READ(44,'(A)') CLINE
      I = INDEX(CLINE,'=')
      J = INDEX(CLINE,'month')
      IF     (J.NE.0) THEN
        READ(CLINE(I+1:),*) J,HMINB,HMAXB
      ELSE
        READ(CLINE(I+1:),*)   HMINB,HMAXB
      ENDIF
      CALL ZAIORD(TXN,MSK,.FALSE., HMINA,HMAXA, 44)
      IF     (ABS(HMINA-HMINB).GT.ABS(HMINB)*1.E-4 .OR.
     &        ABS(HMAXA-HMAXB).GT.ABS(HMAXB)*1.E-4     ) THEN
        WRITE(6,'(/ a / a,1p3e14.6 / a,1p3e14.6 /)')
     &    'error - .a and .b grid files not consistent (txo):',
     &    '.a,.b min = ',HMINA,HMINB,HMINA-HMINB,
     &    '.a,.b max = ',HMAXA,HMAXB,HMAXA-HMAXB
        CALL ZHFLSH(6)
        STOP
      ENDIF
      CLOSE(UNIT=44)
      CALL ZAIOCL(44)
C
      READ(45,'(A)') CLINE
      I = INDEX(CLINE,'=')
      J = INDEX(CLINE,'month')
      IF     (J.NE.0) THEN
        READ(CLINE(I+1:),*) J,HMINB,HMAXB
      ELSE
        READ(CLINE(I+1:),*)   HMINB,HMAXB
      ENDIF
      CALL ZAIORD(TYN,MSK,.FALSE., HMINA,HMAXA, 45)
      IF     (ABS(HMINA-HMINB).GT.ABS(HMINB)*1.E-4 .OR.
     &        ABS(HMAXA-HMAXB).GT.ABS(HMAXB)*1.E-4     ) THEN
        WRITE(6,'(/ a / a,1p3e14.6 / a,1p3e14.6 /)')
     &    'error - .a and .b grid files not consistent (tyo):',
     &    '.a,.b min = ',HMINA,HMINB,HMINA-HMINB,
     &    '.a,.b max = ',HMAXA,HMAXB,HMAXA-HMAXB
        CALL ZHFLSH(6)
        STOP
      ENDIF
      CLOSE(UNIT=45)
      CALL ZAIOCL(45)
C
      CALL MINMAX(TXN,IDM,JDM, XMIN,XMAX)
      CALL AVERMS(TXN,IDM,JDM, XAVE,XRMS)
      CALL MINMAX(TYN,IDM,JDM, YMIN,YMAX)
      CALL AVERMS(TYN,IDM,JDM, YAVE,YRMS)
      WRITE(6,*)
      WRITE(6,*) 'NEW MEAN STATISTICS.'
      WRITE(6,*)
      WRITE(6,8100) SCALEN*XMIN,SCALEN*XMAX,SCALEN*XAVE,SCALEN*XRMS
      WRITE(6,8200) SCALEN*YMIN,SCALEN*YMAX,SCALEN*YAVE,SCALEN*YRMS
      CALL ZHFLSH(6)
C
C     READ IN THE OLD MEAN.
C
      READ(54,'(A)') CLINE
      I = INDEX(CLINE,'=')
      J = INDEX(CLINE,'month')
      IF     (J.NE.0) THEN
        READ(CLINE(I+1:),*) J,HMINB,HMAXB
      ELSE
        READ(CLINE(I+1:),*)   HMINB,HMAXB
      ENDIF
      CALL ZAIORD(TXO,MSK,.FALSE., HMINA,HMAXA, 54)
      IF     (ABS(HMINA-HMINB).GT.ABS(HMINB)*1.E-4 .OR.
     &        ABS(HMAXA-HMAXB).GT.ABS(HMAXB)*1.E-4     ) THEN
        WRITE(6,'(/ a / a,1p3e14.6 / a,1p3e14.6 /)')
     &    'error - .a and .b grid files not consistent (txo):',
     &    '.a,.b min = ',HMINA,HMINB,HMINA-HMINB,
     &    '.a,.b max = ',HMAXA,HMAXB,HMAXA-HMAXB
        CALL ZHFLSH(6)
        STOP
      ENDIF
      CLOSE(UNIT=54)
      CALL ZAIOCL(54)
C
      READ(55,'(A)') CLINE
      I = INDEX(CLINE,'=')
      J = INDEX(CLINE,'month')
      IF     (J.NE.0) THEN
        READ(CLINE(I+1:),*) J,HMINB,HMAXB
      ELSE
        READ(CLINE(I+1:),*)   HMINB,HMAXB
      ENDIF
      CALL ZAIORD(TYO,MSK,.FALSE., HMINA,HMAXA, 55)
      IF     (ABS(HMINA-HMINB).GT.ABS(HMINB)*1.E-4 .OR.
     &        ABS(HMAXA-HMAXB).GT.ABS(HMAXB)*1.E-4     ) THEN
        WRITE(6,'(/ a / a,1p3e14.6 / a,1p3e14.6 /)')
     &    'error - .a and .b grid files not consistent (tyo):',
     &    '.a,.b min = ',HMINA,HMINB,HMINA-HMINB,
     &    '.a,.b max = ',HMAXA,HMAXB,HMAXA-HMAXB
        CALL ZHFLSH(6)
        STOP
      ENDIF
      CLOSE(UNIT=55)
      CALL ZAIOCL(55)
C
      CALL MINMAX(TXO,IDM,JDM, XMIN,XMAX)
      CALL AVERMS(TXO,IDM,JDM, XAVE,XRMS)
      CALL MINMAX(TYO,IDM,JDM, YMIN,YMAX)
      CALL AVERMS(TYO,IDM,JDM, YAVE,YRMS)
      WRITE(6,*)
      WRITE(6,*) 'OLD MEAN STATISTICS.'
      WRITE(6,*)
      WRITE(6,8100) XMIN,XMAX,XAVE,XRMS
      WRITE(6,8200) YMIN,YMAX,YAVE,YRMS
      WRITE(6,*)
      WRITE(6,*)
      CALL ZHFLSH(6)
C
C     OFFSET.
C
      DO 110 J= 1,JDM
        DO 111 I= 1,IDM
          TXH(I,J) = SCALEN*TXN(I,J) - SCALEO*TXO(I,J)
          TYH(I,J) = SCALEN*TYN(I,J) - SCALEO*TYO(I,J)
  111   CONTINUE
  110 CONTINUE
      CALL MINMAX(TXH,IDM,JDM, XMIN,XMAX)
      CALL AVERMS(TXH,IDM,JDM, XAVE,XRMS)
      CALL MINMAX(TYH,IDM,JDM, YMIN,YMAX)
      CALL AVERMS(TYH,IDM,JDM, YAVE,YRMS)
      WRITE(6,*)
      WRITE(6,*) 'OFFSET STATISTICS.'
      WRITE(6,*)
      WRITE(6,8100) XMIN,XMAX,XAVE,XRMS
      WRITE(6,8200) YMIN,YMAX,YAVE,YRMS
      WRITE(6,*)
      WRITE(6,*)
C
C     WRITE OUT OFFSET.
C
      CALL ZAIOWR(TXH,MSK,.FALSE., HMINA,HMAXA, 64, .FALSE.)
      WRITE(64,4102) ' tau_ewd',1,HMINA,HMAXA
      WRITE( 6,4102) ' tau_ewd',1,HMINA,HMAXA
      CALL ZAIOCL(64)
      CLOSE(UNIT=64)
C
      CALL ZAIOWR(TYH,MSK,.FALSE., HMINA,HMAXA, 65, .FALSE.)
      WRITE(65,4102) ' tau_nwd',1,HMINA,HMAXA
      WRITE( 6,4102) ' tau_nwd',1,HMINA,HMAXA
      CALL ZAIOCL(65)
      CLOSE(UNIT=65)
C
      WRITE(6,*)
      CALL ZHFLSH(6)
      STOP
C
 4101 FORMAT(A79)
 4102 FORMAT(A,': month,range = ',I2.2,1P2E16.7)
 6000 FORMAT(1X,A,2X,A40 //)
 6300 FORMAT(10X,'WRITING OFFSET RECORD' /)
 8100 FORMAT(1X,'TAU-X:  MIN =',F10.4,'   MAX =',F10.4,
     +              '    AVE =',F10.4,'   RMS =',F10.4)
 8200 FORMAT(1X,'TAU-Y:  MIN =',F10.4,'   MAX =',F10.4,
     +              '    AVE =',F10.4,'   RMS =',F10.4)
C     END OF PROGRAM OFFDIF.
      END

max_stars_repo_name: TillRasmussen/HYCOM-tools

max_issues_repo_name: TillRasmussen/HYCOM-tools

max_forks_repo_name: TillRasmussen/HYCOM-tools


text:
include 'FUNCTIONS.f90'

Program MARKS

   implicit none

   integer :: EoF
   integer, parameter :: N = 10
   character (len=500) :: line_soluciones
   character (len=8) :: DNI_soluciones
   real :: soluciones_correctas(N)
   character (len=2) :: nota_final
   character (len=500) :: line_grupos
   character (len=8) :: DNI_grupos

   open (UNIT=1, FILE="p1.TXT",STATUS="old")
   open (UNIT=4, FILE="Temp.csv", STATUS="new")
   open (UNIT=5, FILE="NOTASPERDIDAS.csv",STATUS="new")

   call calcular_soluciones('solutions.csv', 3, N, soluciones_correctas)

   do
      read (UNIT=1, IOSTAT=EoF, FMT='(A)') line_soluciones

      if(EoF < 0) then
         exit
      endif

      DNI_soluciones = line_soluciones(1:8)

      call calcular_nota(soluciones_correctas, line_soluciones, N, nota_final)

      open (UNIT=2, FILE='GRUPOS.csv', STATUS='old')
      read (UNIT=2, IOSTAT=EoF, FMT='(A)') line_grupos

      do
         read (UNIT=2, IOSTAT=EoF, FMT='(A)') line_grupos

         ! Notas perdidas
         if(EoF < 0) then
            if(DNI_soluciones /= "") then
                write(UNIT=5,FMT='(A)') DNI_soluciones // ';' // nota_final
            endif
            close(UNIT=2)
            exit
         endif

         call obtener_dni(line_grupos, DNI_grupos)

         if(DNI_soluciones == DNI_grupos) then
            write(UNIT=4, FMT='(A)') line_grupos // ';' // nota_final
            close(UNIT=2)
            exit
         endif

      end do

   end do

   open (UNIT=2, FILE='GRUPOS.csv', STATUS='old')
   close(UNIT=2, status='delete')

   close(UNIT=1)
   close(UNIT=4)
   close(UNIT=5)

   call RENAME('Temp.csv','GRUPOS.csv')
end program

max_stars_repo_name: jmorenov/notas-alumnos

max_issues_repo_name: jmorenov/notas-alumnos

max_forks_repo_name: jmorenov/notas-alumnos


text:
SUBROUTINE DRESSEDBASIS(D,ID,NM,MODES_NUM,FIELDS,U_FD,E_DRESSED,INFO)

!!$THIS SUBROUTINES CALCULATES THE TRANSFORMATION BETWEEN THE BARE BASIS TO THE DRESSED BASIS DEFINDED BY THE FULL SET OF DRIVING FIELDS.
!!$ D                            : DIMENSION OF THE MULTIMODE EXTENDED HILBERT SPACE
!!$ ID (IN)                      : TYPE OF QUATUM SYSTEM
!!$ NM (IN)                      : NUMBER OF MODES == NUMBER OF DRIVING FIELDS
!!$ MODES_NUM                    : VECTOR INDICATING THE NUMBER OF HARMONICS OF EACH DRESSING FIELD
!!$ FIELDS (IN)                  : AMPLITUDE, FREQUENCY AND PHASES OF ALL DRIVING FIELDS
!!$ U_FD (OUT)                   : THIS IS THE TRANSFORMATION WE ARE LOOKING FOR
!!$ E_DRESSED (OUT)              : DRESSED ENERGIES
!!$ INFO (INOUT)                 : INFO = 0 MEANS SUCESS

  USE ATOMIC_PROPERTIES
  USE TYPES
  USE SUBINTERFACE
  USE SUBINTERFACE_LAPACK
!  USE FLOQUETINIT_
  USE ARRAYS

  IMPLICIT NONE
  TYPE(MODE), DIMENSION(NM),     INTENT(IN)    :: FIELDS
  TYPE(ATOM),                    INTENT(IN)    :: ID
  INTEGER,    DIMENSION(NM),     INTENT(IN)    :: MODES_NUM
  COMPLEX*16, DIMENSION(D,D),       INTENT(OUT)   :: U_FD
  DOUBLE PRECISION, DIMENSION(D), INTENT(OUT)   :: E_DRESSED
  INTEGER,                       INTENT(IN)    :: NM,D
  INTEGER,                       INTENT(INOUT) :: INFO

  INTEGER r,m,FIELD_INDEX
  !INTEGER ND_,NM_
  INTEGER TOTAL_FREQUENCIES
  !TYPE(MODE), DIMENSION(DRESSINGFIELDS) :: FIELDS_
  !TYPE(MODE), DIMENSION(:),ALLOCATABLE  :: FIELDS_
  !INTEGER,    DIMENSION(DRESSINGFIELDS) :: MODES_NUM_

  !NM_ = DRESSINGFIELDS
  !ND_ = DRESSINGFIELDS

  TOTAL_FREQUENCIES = SUM(MODES_NUM,1)

  CALL SETHAMILTONIANCOMPONENTS(ID,NM,TOTAL_FREQUENCIES,MODES_NUM,FIELDS,INFO)

  !--- FIND THE MULTIMODE FLOQUET SPECTRUM
  CALL MULTIMODEFLOQUETMATRIX(ID,NM,TOTAL_FREQUENCIES,MODES_NUM,FIELDS,INFO)

  E_DRESSED = 0.0
  CALL LAPACK_FULLEIGENVALUES(H_FLOQUET,SIZE(H_FLOQUET,1),E_DRESSED,INFO)

  U_FD = H_FLOQUET ! FOURIER DECOMPOSITION OF THE STATES DRESSED BY MODE NUMBER
  DEALLOCATE(H_FLOQUET)

END SUBROUTINE DRESSEDBASIS

SUBROUTINE DRESSEDBASIS_SUBSET(ID,DRESSINGFLOQUETDIMENSION,DRESSINGFIELDS,NM,&
        & DRESSINGFIELDS_INDICES,MODES_NUM,FIELDS,U_FD,E_DRESSED,INFO)

!!$THIS SUBROUTINES CALCULATES THE TRANSFORMATION BETWEEN THE BARE BASIS TO THE DRESSED BASIS DEFINDED BY A SUBSET OF THE DRIVING FIELDS.
!!$
!!$ ID (IN)                      : TYPE OF QUATUM SYSTEM
!!$ DRESSINGFLOQUETDIMENSION(IN) : DIMENSION OF THE DRESSED FLOQUET HILBERT SPACE
!!$ DRESSINGFIELDS (IN)          : NUMBER OF DRESSING FIELDS
!!$ NM (IN)                      : NUMBER OF MODES == NUMBER OF DRIVING FIELDS
!!$ DRESSINGFIELDS_INDICES       : OUT OF ALL THE DRIVING FIELDS, THIS VECTOR INDICATES THE INDEX OF THE DRESSING FIELD
!!$ MODES_NUM                    : VECTOR INDICATING THE NUMBER OF HARMONICS OF EACH DRESSING FIELD
!!$ FIELDS (IN)                  : AMPLITUDE, FREQUENCY AND PHASES OF ALL DRIVING FIELDS
!!$ U_FD (OUT)                   : THIS IS THE TRANSFORMATION WE ARE LOOKING FOR
!!$ E_DRESSED (OUT)              : DRESSED ENERGIES
!!$ INFO (INOUT)                 : INFO = 0 MEANS SUCESS

  USE ATOMIC_PROPERTIES
  USE TYPES
  USE SUBINTERFACE
  USE SUBINTERFACE_LAPACK
!  USE FLOQUETINIT_
  USE ARRAYS

  IMPLICIT NONE
  INTEGER,                                               INTENT(IN)    :: NM
  INTEGER,                                               INTENT(IN)    :: DRESSINGFLOQUETDIMENSION,DRESSINGFIELDS
  TYPE(MODE), DIMENSION(NM),                             INTENT(IN)    :: FIELDS
  TYPE(ATOM),                                            INTENT(IN)    :: ID
  INTEGER,    DIMENSION(DRESSINGFIELDS),                 INTENT(IN)    :: DRESSINGFIELDS_INDICES
  INTEGER,    DIMENSION(NM),                             INTENT(IN)    :: MODES_NUM
  COMPLEX*16, DIMENSION(DRESSINGFLOQUETDIMENSION,DRESSINGFLOQUETDIMENSION),       INTENT(OUT)   :: U_FD
  DOUBLE PRECISION, DIMENSION(DRESSINGFLOQUETDIMENSION), INTENT(OUT)   :: E_DRESSED
  INTEGER,                                               INTENT(INOUT) :: INFO

  INTEGER r,m,FIELD_INDEX
  INTEGER ND_,NM_
  INTEGER TOTAL_FREQUENCIES_
  TYPE(MODE), DIMENSION(:),ALLOCATABLE  :: FIELDS_
  INTEGER,    DIMENSION(DRESSINGFIELDS) :: MODES_NUM_

!  write(*,*)"non-SP A:", DRESSINGFIELDS_INDICES

  DO r=1,DRESSINGFIELDS
     MODES_NUM_(r)=modes_num(DRESSINGFIELDS_INDICES(r))
  END DO

  TOTAL_FREQUENCIES_ = SUM(MODES_NUM_,1)
  ALLOCATE(FIELDS_(TOTAL_FREQUENCIES_))
  DO m=1,TOTAL_FREQUENCIES_
     ALLOCATE(FIELDS_(m)%V(ID%D_BARE,ID%D_BARE))
  END DO
 ! write(*,*) "non-SP:",dressingfields,modes_num_,total_frequencies_

  FIELD_INDEX = 1
  DO r=1,DRESSINGFIELDS
     DO m=1,MODES_NUM_(r)
        FIELDS_(FIELD_INDEX) = FIELDS(DRESSINGFIELDS_INDICES(r)+m-1)
        !FIELDS_(FIELD_INDEX)%N_Floquet = FIELDS(DRESSINGFIELDS_INDICES(r)+m-1)%N_FLOQUET
        !write(*,*) r,m,field_index,DRESSINGFIELDS_INDICES
        FIELD_INDEX = FIELD_INDEX+1
     END DO
  END DO

  CALL SETHAMILTONIANCOMPONENTS(ID,size(modes_num_,1),total_frequencies_,MODES_NUM_,FIELDS_,INFO)

!!$  !--- FIND THE MULTIMODE FLOQUET SPECTRUM
!  write(*,*) fields_(2)%N_Floquet, fields(2)%N_floquet
  !fields_(2)%N_Floquet = fields(2)%N_floquet
  CALL MULTIMODEFLOQUETMATRIX(ID,size(modes_num_,1),total_frequencies_,MODES_NUM_,FIELDS_,INFO)
  E_DRESSED = 0.0
  CALL LAPACK_FULLEIGENVALUES(H_FLOQUET,SIZE(H_FLOQUET,1),E_DRESSED,INFO)
  U_FD = H_FLOQUET ! FOURIER DECOMPOSITION OF THE STATES DRESSED BY MODE NUMBER
  DEALLOCATE(H_FLOQUET)
  DEALLOCATE(FIELDS_)

END SUBROUTINE DRESSEDBASIS_SUBSET

max_stars_repo_name: gsinuco/OPENMMF

max_issues_repo_name: gsinuco/OPENMMF

max_forks_repo_name: gsinuco/OPENMMF


text:
! { dg-do compile { target i?86-*-* x86_64-*-* } }
! { dg-require-effective-target vect_double }
! { dg-require-effective-target sse2 }
! { dg-options "-O3 -ffast-math -msse2 -fpredictive-commoning -ftree-vectorize -fdump-tree-optimized" }

******* RESID COMPUTES THE RESIDUAL:  R = V - AU
*
*      THIS SIMPLE IMPLEMENTATION COSTS  27A + 4M PER RESULT, WHERE
*      A AND M DENOTE THE COSTS OF ADDITION (OR SUBTRACTION) AND
*      MULTIPLICATION, RESPECTIVELY.  BY USING SEVERAL TWO-DIMENSIONAL
*      BUFFERS ONE CAN REDUCE THIS COST TO  13A + 4M IN THE GENERAL
*      CASE, OR  10A + 3M WHEN THE COEFFICIENT A(1) IS ZERO.
*
      SUBROUTINE RESID(U,V,R,N,A)
      INTEGER N
      REAL*8 U(N,N,N),V(N,N,N),R(N,N,N),A(0:3)
      INTEGER I3, I2, I1
C
      DO 600 I3=2,N-1
      DO 600 I2=2,N-1
      DO 600 I1=2,N-1
 600  R(I1,I2,I3)=V(I1,I2,I3)
     >      -A(0)*( U(I1,  I2,  I3  ) )
     >      -A(1)*( U(I1-1,I2,  I3  ) + U(I1+1,I2,  I3  )
     >                 +  U(I1,  I2-1,I3  ) + U(I1,  I2+1,I3  )
     >                 +  U(I1,  I2,  I3-1) + U(I1,  I2,  I3+1) )
     >      -A(2)*( U(I1-1,I2-1,I3  ) + U(I1+1,I2-1,I3  )
     >                 +  U(I1-1,I2+1,I3  ) + U(I1+1,I2+1,I3  )
     >                 +  U(I1,  I2-1,I3-1) + U(I1,  I2+1,I3-1)
     >                 +  U(I1,  I2-1,I3+1) + U(I1,  I2+1,I3+1)
     >                 +  U(I1-1,I2,  I3-1) + U(I1-1,I2,  I3+1)
     >                 +  U(I1+1,I2,  I3-1) + U(I1+1,I2,  I3+1) )
     >      -A(3)*( U(I1-1,I2-1,I3-1) + U(I1+1,I2-1,I3-1)
     >                 +  U(I1-1,I2+1,I3-1) + U(I1+1,I2+1,I3-1)
     >                 +  U(I1-1,I2-1,I3+1) + U(I1+1,I2-1,I3+1)
     >                 +  U(I1-1,I2+1,I3+1) + U(I1+1,I2+1,I3+1) )
C
      RETURN
      END
! we want to check that predictive commoning did something on the
! vectorized loop, which means we have to have exactly 13 vector
! additions.
! { dg-final { scan-tree-dump-times "vect_var\[^\\n\]*\\+ " 13 "optimized" } }
! { dg-final { cleanup-tree-dump "vect" } }
! { dg-final { cleanup-tree-dump "optimized" } }

max_stars_repo_name: maurizioabba/rose

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


text:
      SUBROUTINE SIGTST (EXPSIG, FAILED)

C  Sigtst is used to compare an signalled error code with an expected
C  value.  If not equal, failure and informative messages are generated,
C  including the error code and function name.  If equal, pass is
C  invoked.

C  Input parameter:
C    EXPSIG : expected value of error code (usually 0)
C  Output parameter
C    FAILED : true if expected not equal to actual, else false

      COMMON /GLOBNU/ CTLHND, ERRSIG, ERRFIL, IERRCT, UNERR,
     1        TESTCT, IFLERR, PASSSW, ERRSW, MAXLIN,
     2        CONID, MEMUN, WKID, WTYPE, GLBLUN, INDLUN,
     3        DUMINT, DUMRL
      INTEGER         CTLHND, ERRSIG, ERRFIL, IERRCT, UNERR,
     1        TESTCT, IFLERR, PASSSW, ERRSW, MAXLIN,
     2        CONID, MEMUN, WKID, WTYPE, GLBLUN, INDLUN,
     3        DUMINT(20), ERRIND
      REAL    DUMRL(20)

      COMMON /GLOBCH/ PIDENT,    GLBERR,    TSTMSG,     FUNCID,
     1                DUMCH
      CHARACTER       PIDENT*40, GLBERR*60, TSTMSG*900, FUNCID*80,
     1                DUMCH(20)*20

      INTEGER    EXPSIG
      CHARACTER  MSG*300
      LOGICAL    FAILED

      IF (EXPSIG .EQ. ERRSIG) THEN
         FAILED = .FALSE.
         CALL PASS
      ELSE
         FAILED = .TRUE.
         CALL FAIL
         CALL SIGMSG (ERRSIG, FUNCID, MSG)
         CALL INMSG (MSG)
      ENDIF

      END