text:
      SUBROUTINE XYZINT(XYZ,NUMAT,NA,NB,NC,DEGREE,GEO)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION XYZ(3,*), NA(*), NB(*), NC(*), GEO(3,*)
***********************************************************************
*
* XYZINT WORKS OUT THE INTERNAL COORDINATES OF A MOLECULE.
*        THE "RULES" FOR THE CONNECTIVITY ARE AS FOLLOWS:
*        ATOM I IS DEFINED AS BEING AT A DISTANCE FROM THE NEAREST
*        ATOM J, ATOM J ALREADY HAVING BEEN DEFINED.
*        ATOM I MAKES AN ANGLE WITH ATOM J AND THE ATOM K, WHICH HAS
*        ALREADY BEEN DEFINED, AND IS THE NEAREST ATOM TO J
*        ATOM I MAKES A DIHEDRAL ANGLE WITH ATOMS J, K, AND L. L HAVING
*        BEEN DEFINED AND IS THE NEAREST ATOM TO K, AND J, K AND L
*        HAVE A CONTAINED ANGLE IN THE RANGE 15 TO 165 DEGREES,
*        IF POSSIBLE.
*
*        IF(NA(2).EQ.-1 OR -2 THEN THE ORIGINAL CONNECTIVITY IS USED.
*
*        NOTE THAT GEO AND XYZ MUST NOT BE THE SAME IN THE CALL.
*
*   ON INPUT XYZ    = CARTESIAN ARRAY OF NUMAT ATOMS
*            DEGREE = 1 IF ANGLES ARE TO BE IN RADIANS
*            DEGREE = 57.29578 IF ANGLES ARE TO BE IN DEGREES
*
***********************************************************************
      COMMON /GEOOK/ IGEOOK
      COMMON /NUMCAL/ NUMCAL
      DATA ICALCN/0/
      IGEOOK=99
      IF(.NOT.(ICALCN.NE.NUMCAL).AND.NA(2).EQ.-1 .OR. NA(2).EQ.-2)THEN
         NA(2)=1
         DO 10 I=2,NUMAT
            J=NA(I)
            IF(I.GT.3)CALL DIHED(XYZ,I,J,NB(I),NC(I),GEO(3,I))
            IF(I.GT.2)CALL BANGLE(XYZ,I,J,NB(I),GEO(2,I))
            GEO(1,I)=SQRT((XYZ(1,I)-XYZ(1,J))**2+
     1          (XYZ(2,I)-XYZ(2,J))**2+
     2          (XYZ(3,I)-XYZ(3,J))**2)
   10    CONTINUE
      ELSE
         IF(NA(2).EQ.-1)ICALCN=NUMCAL
         DO 30 I=1,NUMAT
            NA(I)=2
            NB(I)=3
            NC(I)=4
            IM1=I-1
            IF(IM1.EQ.0)GOTO 30
            SUM=1.D30
            DO 20 J=1,IM1
               R=(XYZ(1,I)-XYZ(1,J))**2+
     1          (XYZ(2,I)-XYZ(2,J))**2+
     2          (XYZ(3,I)-XYZ(3,J))**2
               IF(R.LT.SUM.AND.NA(J).NE.J.AND.NB(J).NE.J) THEN
                  SUM=R
                  K=J
               ENDIF
   20       CONTINUE
C
C   ATOM I IS NEAREST TO ATOM K
C
            NA(I)=K
            IF(I.GT.2)NB(I)=NA(K)
            IF(I.GT.3)NC(I)=NB(K)
C
C   FIND ANY ATOM TO RELATE TO NA(I)
C
   30    CONTINUE
      ENDIF
      NA(1)=0
      NB(1)=0
      NC(1)=0
      NB(2)=0
      NC(2)=0
      NC(3)=0
      CALL XYZGEO(XYZ,NUMAT,NA,NB,NC,DEGREE,GEO)
      RETURN
      END
      SUBROUTINE XYZGEO(XYZ,NUMAT,NA,NB,NC,DEGREE,GEO)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION XYZ(3,*), NA(*), NB(*), NC(*), GEO(3,*)
***********************************************************************
*
*   XYZGEO CONVERTS COORDINATES FROM CARTESIAN TO INTERNAL.
*
*     ON INPUT XYZ  = ARRAY OF CARTESIAN COORDINATES
*              NUMAT= NUMBER OF ATOMS
*              NA   = NUMBERS OF ATOM TO WHICH ATOMS ARE RELATED
*                     BY DISTANCE
*              NB   = NUMBERS OF ATOM TO WHICH ATOMS ARE RELATED
*                     BY ANGLE
*              NC   = NUMBERS OF ATOM TO WHICH ATOMS ARE RELATED
*                     BY DIHEDRAL
*
*    ON OUTPUT GEO  = INTERNAL COORDINATES IN ANGSTROMS, RADIANS,
*                     AND RADIANS
*
***********************************************************************
      DO 30 I=2,NUMAT
         J=NA(I)
         K=NB(I)
         L=NC(I)
         IF(I.LT.3) GOTO 30
         II=I
         CALL BANGLE(XYZ,II,J,K,GEO(2,I))
         GEO(2,I)=GEO(2,I)*DEGREE
         IF(I.LT.4) GOTO 30
C
C   MAKE SURE DIHEDRAL IS MEANINGLFUL
C
         CALL BANGLE(XYZ,J,K,L,ANGL)
         TOL=0.2617994D0
         IF(ANGL.GT.3.1415926D0-TOL.OR.ANGL.LT.TOL)THEN
C
C  ANGLE IS UNSATISFACTORY, LET'S SEARCH FOR ANOTHER ATOM FOR
C  DEFINING THE DIHEDRAL.
   10       SUM=100.D0
            DO 20 I1=1,II-1
               R=(XYZ(1,I1)-XYZ(1,K))**2+
     1          (XYZ(2,I1)-XYZ(2,K))**2+
     2          (XYZ(3,I1)-XYZ(3,K))**2
               IF(R.LT.SUM.AND.I1.NE.J.AND.I1.NE.K) THEN
                  CALL BANGLE(XYZ,J,K,I1,ANGL)
                  IF(ANGL.LT.3.1415926D0-TOL.AND.ANGL.GT.TOL)THEN
                     SUM=R
                     L=I1
                     NC(II)=L
                  ENDIF
               ENDIF
   20       CONTINUE
            IF(SUM.GT.99.D0.AND.TOL.GT.0.1D0)THEN
C
C ANYTHING WITHIN 5 DEGREES?
C
               TOL=0.087266D0
               GOTO 10
            ENDIF
         ENDIF
         CALL DIHED(XYZ,II,J,K,L,GEO(3,I))
         GEO(3,I)=GEO(3,I)*DEGREE
   30 GEO(1,I)= SQRT((XYZ(1,I)-XYZ(1,J))**2+
     1                   (XYZ(2,I)-XYZ(2,J))**2+
     2                   (XYZ(3,I)-XYZ(3,J))**2)
      GEO(1,1)=0.D0
      GEO(2,1)=0.D0
      GEO(3,1)=0.D0
      GEO(2,2)=0.D0
      GEO(3,2)=0.D0
      GEO(3,3)=0.D0
      RETURN
      END
      SUBROUTINE BANGLE(XYZ,I,J,K,ANGLE)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION XYZ(3,*)
*********************************************************************
*
* BANGLE CALCULATES THE ANGLE BETWEEN ATOMS I,J, AND K. THE
*        CARTESIAN COORDINATES ARE IN XYZ.
*
*********************************************************************
      D2IJ = (XYZ(1,I)-XYZ(1,J))**2+
     1       (XYZ(2,I)-XYZ(2,J))**2+
     2       (XYZ(3,I)-XYZ(3,J))**2
      D2JK = (XYZ(1,J)-XYZ(1,K))**2+
     1       (XYZ(2,J)-XYZ(2,K))**2+
     2       (XYZ(3,J)-XYZ(3,K))**2
      D2IK = (XYZ(1,I)-XYZ(1,K))**2+
     1       (XYZ(2,I)-XYZ(2,K))**2+
     2       (XYZ(3,I)-XYZ(3,K))**2
      XY = SQRT(D2IJ*D2JK)
      TEMP = 0.5D0 * (D2IJ+D2JK-D2IK) / XY
      IF (TEMP .GT. 1.0D0) TEMP=1.0D0
      IF (TEMP .LT. -1.0D0) TEMP=-1.0D0
      ANGLE = ACOS( TEMP )
      RETURN
      END
      SUBROUTINE DIHED(XYZ,I,J,K,L,ANGLE)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION XYZ(3,*)
*********************************************************************
*
*      DIHED CALCULATES THE DIHEDRAL ANGLE BETWEEN ATOMS I, J, K,
*            AND L.  THE CARTESIAN COORDINATES OF THESE ATOMS
*            ARE IN ARRAY XYZ.
*
*     DIHED IS A MODIFIED VERSION OF A SUBROUTINE OF THE SAME NAME
*           WHICH WAS WRITTEN BY DR. W. THEIL IN 1973.
*
*********************************************************************
      XI1=XYZ(1,I)-XYZ(1,K)
      XJ1=XYZ(1,J)-XYZ(1,K)
      XL1=XYZ(1,L)-XYZ(1,K)
      YI1=XYZ(2,I)-XYZ(2,K)
      YJ1=XYZ(2,J)-XYZ(2,K)
      YL1=XYZ(2,L)-XYZ(2,K)
      ZI1=XYZ(3,I)-XYZ(3,K)
      ZJ1=XYZ(3,J)-XYZ(3,K)
      ZL1=XYZ(3,L)-XYZ(3,K)
C      ROTATE AROUND Z AXIS TO PUT KJ ALONG Y AXIS
      DIST= SQRT(XJ1**2+YJ1**2+ZJ1**2)
      COSA=ZJ1/DIST
      IF(COSA.GT.1.0D0) COSA=1.0D0
      IF(COSA.LT.-1.0D0) COSA=-1.0D0
      DDD=1.0D0-COSA**2
      IF(DDD.LE.0.0) GO TO 10
      YXDIST=DIST* SQRT(DDD)
      IF(YXDIST.GT.1.0D-6) GO TO 20
   10 CONTINUE
      XI2=XI1
      XL2=XL1
      YI2=YI1
      YL2=YL1
      COSTH=COSA
      SINTH=0.D0
      GO TO 30
   20 COSPH=YJ1/YXDIST
      SINPH=XJ1/YXDIST
      XI2=XI1*COSPH-YI1*SINPH
      XL2=XL1*COSPH-YL1*SINPH
      YI2=XI1*SINPH+YI1*COSPH
      YJ2=XJ1*SINPH+YJ1*COSPH
      YL2=XL1*SINPH+YL1*COSPH
C      ROTATE KJ AROUND THE X AXIS SO KJ LIES ALONG THE Z AXIS
      COSTH=COSA
      SINTH=YJ2/DIST
   30 CONTINUE
      YI3=YI2*COSTH-ZI1*SINTH
      YL3=YL2*COSTH-ZL1*SINTH
      CALL DANG(XL2,YL3,XI2,YI3,ANGLE)
      IF (ANGLE .LT. 0.) ANGLE=4.0D0* ASIN(1.0D00)+ANGLE
      IF (ANGLE .GE. 6.2831853D0 ) ANGLE=0.D0
      RETURN
      END
      SUBROUTINE DANG(A1,A2,B1,B2,RCOS)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
**********************************************************************
*
*    DANG  DETERMINES THE ANGLE BETWEEN THE POINTS (A1,A2), (0,0),
*          AND (B1,B2).  THE RESULT IS PUT IN RCOS.
*
**********************************************************************
      ZERO=1.0D-6
      IF( ABS(A1).LT.ZERO.AND. ABS(A2).LT.ZERO) GO TO 10
      IF( ABS(B1).LT.ZERO.AND. ABS(B2).LT.ZERO) GO TO 10
      ANORM=1.0D0/ SQRT(A1**2+A2**2)
      BNORM=1.0D0/ SQRT(B1**2+B2**2)
      A1=A1*ANORM
      A2=A2*ANORM
      B1=B1*BNORM
      B2=B2*BNORM
      SINTH=(A1*B2)-(A2*B1)
      COSTH=A1*B1+A2*B2
      IF(COSTH.GT.1.0D0) COSTH=1.0D0
      IF(COSTH.LT.-1.0D0) COSTH=-1.0D0
      RCOS= ACOS(COSTH)
      IF( ABS(RCOS).LT.4.0D-4) GO TO 10
      IF(SINTH.GT.0.D0) RCOS=4.0D0* ASIN(1.0D00)-RCOS
      RCOS=-RCOS
      RETURN
   10 RCOS=0.0D0
      RETURN
      END

max_stars_repo_name: openmopac/MOPAC-archive

max_issues_repo_name: openmopac/MOPAC-archive

max_forks_repo_name: openmopac/MOPAC-archive


text:
CHARACTER(3) :: a, my_string
a = my_string(2:3)
END PROGRAM

max_stars_repo_name: OpenFortranProject/ofp-sdf

max_issues_repo_name: OpenFortranProject/ofp-sdf

max_forks_repo_name: OpenFortranProject/ofp-sdf


text:
// Includes
+incdir+../source/

//Source File (PKG, RTL, MACROS)
../FPGA_CPUGarage/sv/cpu_pkg.sv
../FPGA_CPUGarage/sv/cpu.sv
../FPGA_CPUGarage/sv/alu.sv
../FPGA_CPUGarage/sv/div.sv
../FPGA_CPUGarage/sv/ctrl.sv
../FPGA_CPUGarage/sv/acc_machine.sv

//Test Bench - simulation
../source/cpu_garage_tb.sv
../source/cpu_garage.sv
../source/definitions.sv
../source/ram.sv
../source/rom.sv

ext: f

max_stars_repo_name: amichai-bd/CPUGarageChallenge

max_issues_repo_name: amichai-bd/CPUGarageChallenge

max_forks_repo_name: amichai-bd/CPUGarageChallenge


text:
module M {

  enum E { A, B, C }

}

ext: fpp

max_stars_repo_name: kevin-f-ortega/fpp

max_issues_repo_name: kevin-f-ortega/fpp

max_forks_repo_name: kevin-f-ortega/fpp


text:
!
!     CalculiX - A 3-dimensional finite element program
!              Copyright (C) 1998-2018 Guido Dhondt
!
!     This program is free software; you can redistribute it and/or
!     modify it under the terms of the GNU General Public License as
!     published by the Free Software Foundation(version 2);
!
!
!     This program is distributed in the hope that it will be useful,
!     but WITHOUT ANY WARRANTY; without even the implied warranty of
!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
!     GNU General Public License for more details.
!
!     You should have received a copy of the GNU General Public License
!     along with this program; if not, write to the Free Software
!     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
!
      subroutine desiperelem(ndesi,istartdesi,ialdesi,ipoeldi,ieldi,
     &      ne,istartelem,ialelem)
!
      implicit none
!
      integer ndesi,istartdesi(*),ialdesi(*),ipoeldi(*),ieldi(2,*),
     &  ieldifree,i,j,nelem,ne,ifree,index,istartelem(*),ialelem(*)
!
      intent(in) ndesi,istartdesi,ialdesi,ne
!
      intent(inout) ipoeldi,ieldi,istartelem,ialelem
!
!     storing the design variables per element
!
      ieldifree=1
      do i=1,ndesi
         do j=istartdesi(i),istartdesi(i+1)-1
            nelem=ialdesi(j)
            ieldi(1,ieldifree)=i
            ieldi(2,ieldifree)=ipoeldi(nelem)
            ipoeldi(nelem)=ieldifree
            ieldifree=ieldifree+1
         enddo
      enddo
!
!     adding the zero design variable to all elements with
!     a nonzero ipoeldi value
!
      do i=1,ne
         if(ipoeldi(i).eq.0) cycle
         ieldi(1,ieldifree)=0
         ieldi(2,ieldifree)=ipoeldi(i)
         ipoeldi(i)=ieldifree
         ieldifree=ieldifree+1
      enddo
!
!     determining the design variables belonging to a given
!     element i. They are stored in ialelem(istartelem(i))..
!     ...up to..... ialdesi(istartelem(i+1)-1)
!
      ifree=1
      do i=1,ne
         istartelem(i)=ifree
         index=ipoeldi(i)
         do
            if(index.eq.0) exit
            ialelem(ifree)=ieldi(1,index)
c            write(*,*) 'desiperelem ',i,ialelem(ifree)
            ifree=ifree+1
            index=ieldi(2,index)
         enddo
      enddo
      istartelem(ne+1)=ifree
!
      return
      end

max_stars_repo_name: alleindrach/calculix-desktop

max_issues_repo_name: alleindrach/calculix-desktop

max_forks_repo_name: alleindrach/calculix-desktop


text:
C  =====================================================================
C  pgm: SHDECB .. Decode shef .B format messages
C
C  use:     CALL SHDECB(KHAR,KHPOS,IREV,ICONT,QUO)
C
C  i/o: KHAR ...... last buffer char obtained - CHAR*1
C  i/o: KHPOS ..... last char loc: 2=eol,1=err-eol,0=eof,neg=err - INT
C   in: IREV ...... revision code (0 = no, 1 = yes) - INT
C   in: ICONT ..... continuation code (0 = no, 1 = yes) - INT
C   in: QUO ....... string for quotes, only carried thru this rtn in
C   in:             order to define its size in one place - CHAR*(*)
C
C  rqd: SH2BEG,SH2ST2,SHDOTB,SHERR,SHGETK,SHPCOD,SHPDEC,SHPOS
C  rqd: SHTYPA,SHSAVQ,SHSAVD,SHSAVI,SHSAVK,SHSAVL,SHSAVM,SHSAVP
C  rqd: SHSAVJ,SHSAVN,SHSAVA,SHSAVS
C
C  cmt: Hardwired for 200 fields (in stmt comparing "MREC .GE. 200").
C  =====================================================================
      SUBROUTINE SHDECB(KHAR,KHPOS,IREV,ICONT,QUO)

      EXTERNAL       SH2BEG,SH2ST2,SHDOTB,SHERR,SHGETK,SHSAVA,SHSAVS
      EXTERNAL       SHPCOD,SHPDEC,SHPOS,SHTYPA,SHSAVQ,SHSAVJ,SHSAVN
      EXTERNAL       SHSAVD,SHSAVI,SHSAVK,SHSAVL,SHSAVM,SHSAVP

      CHARACTER*1    KHAR,KKK
      CHARACTER*8    KHWRD,PARCOD,JKHID,PPP
      CHARACTER*(*)  QUO
      REAL           CODP
      INTEGER        KHPOS,IREV,ICONT,MREC,NOC,III

      SAVE           JKHID,MREC,CODP
C
C    ================================= RCS keyword statements ==========
      CHARACTER*68     RCSKW1,RCSKW2
      DATA             RCSKW1,RCSKW2 /                                 '
     .$Source: /fs/hseb/ob5/rfc/ofs/src/shefpars_driv/RCS/shdecb.f,v $
     . $',                                                             '
     .$Id: shdecb.f,v 1.6 2000/03/17 19:58:12 dws Exp $
     . $' /
C    ===================================================================
C

C                   If first line messg, set defaults, get pos data

      IF (ICONT .EQ. 0) THEN
        CALL SHSAVK('I',III,III,III,III,III,III,III)
        CALL SHSAVM('I',III,III,III,III,III,III,III,III)
        CALL SHSAVQ('I',III,KKK)
        CALL SHSAVD('I',III,III)
        CALL SHSAVI('I',III,III)
        CALL SHSAVP('I',III,PPP)
         MREC   = 0
        CALL SHPOS(KHAR,KHPOS,JKHID)
      ENDIF

C                   Get the date and data type elements, parameter code

   30 IF (KHPOS .LE. 2) GOTO 34

        CALL SHTYPA(KHAR,KHPOS)
        CALL SHPCOD(KHAR,KHPOS)
        CALL SHPDEC(KHAR,KHPOS,CODP)

C                   Save data for next field using "SHSAV" rtns

        CALL SHSAVP('G',III,PARCOD)
        IF (KHPOS.GT.2 .OR. PARCOD.NE.'        ') THEN
            IF (PARCOD .EQ. '        ') THEN
              CALL SHERR('E',50,KHPOS,KHAR)
            ELSEIF (KHPOS .NE. 1) THEN
              IF (MREC .GE. 200) THEN
                CALL SHERR('E',40,KHPOS,KHAR)
              ELSE
                MREC = MREC+1
                CALL SHSAVP('S',MREC,PPP)
                CALL SHSAVL('S',MREC,III,III,III,III,III,III)
                CALL SHSAVK('S',MREC,III,III,III,III,III,III)
                CALL SHSAVM('S',MREC,III,III,III,III,III,III,III)
                CALL SHSAVA('S',MREC,III)
                CALL SHSAVS('S',MREC,III)
                CALL SHSAVN('S',MREC,III)
                CALL SHSAVJ('S',MREC,III)
                CALL SHSAVD('S',MREC,III)
                CALL SHSAVI('S',MREC,III)
                CALL SHSAVQ('S',MREC,KKK)
              ENDIF

              CALL SHSAVP('I',III,PPP)
            ENDIF

            IF (KHAR .EQ. '/') THEN
              CALL SHGETK(KHAR,KHPOS)
            ELSEIF (KHPOS .GT. 2) THEN
              CALL SHERR('E',33,KHPOS,KHAR)
            ENDIF
        ENDIF

        GOTO 30
   34 CONTINUE

C                   Check next line to see if it contains body data

      IF (KHPOS .NE. 1) THEN
        CALL SH2BEG(KHAR,KHPOS)

        IF (KHPOS .EQ. 0) THEN
          CALL SHERR('W',46,KHPOS,KHAR)
        ELSEIF (MREC.GT.0 .AND. KHPOS.GT.0 .AND. KHAR.NE.'.') THEN
          CALL SHDOTB(KHAR,KHPOS,MREC,IREV,JKHID,QUO)

C                   If KHPOS=0 then eof, else have dot line so see
C                   if ".END" otherwise the ".END" line is missing.
C                   This is a fatal error if we have a dot line
C                   with a continuation code, else just a warning.

          IF (KHPOS .EQ. 0) THEN
            CALL SHERR('W',46,KHPOS,KHAR)
          ELSEIF (KHPOS .GT. 1) THEN
            CALL SHGETK(KHAR,KHPOS)
            CALL SH2ST2(KHAR,KHPOS,NOC,KHWRD)
            IF (KHWRD .EQ. 'END     ') THEN
              KHPOS = 2
              KHAR  = ':'
            ELSEIF (KHWRD(1:1).EQ.'B' .AND. KHWRD(2:2).NE.' ') THEN
              CALL SHERR('E',46,KHPOS,KHAR)
            ELSE
              CALL SHERR('W',46,KHPOS,KHAR)
              KHPOS = 3
              CALL SHGETK(KHAR,KHPOS)
            ENDIF
          ENDIF

        ELSEIF (MREC.EQ.0 .AND. KHAR.NE.'.') THEN
            CALL SHERR('W',82,KHPOS,KHAR)

        ENDIF
      ENDIF

      RETURN
      END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


text:
!!  Copyright (C)  Stichting Deltares, 2012-2016.
!!
!!  This program is free software: you can redistribute it and/or modify
!!  it under the terms of the GNU General Public License version 3,
!!  as published by the Free Software Foundation.
!!
!!  This program is distributed in the hope that it will be useful,
!!  but WITHOUT ANY WARRANTY; without even the implied warranty of
!!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
!!  GNU General Public License for more details.
!!
!!  You should have received a copy of the GNU General Public License
!!  along with this program. If not, see <http://www.gnu.org/licenses/>.
!!
!!  contact: delft3d.support@deltares.nl
!!  Stichting Deltares
!!  P.O. Box 177
!!  2600 MH Delft, The Netherlands
!!
!!  All indications and logos of, and references to registered trademarks
!!  of Stichting Deltares remain the property of Stichting Deltares. All
!!  rights reserved.

!
!  *********************************************************************
!  *        SUBROUTINE TO PRINT MAXIMUM SOLUTION ON TAPE IOU(6)        *
!  *********************************************************************
!
      SUBROUTINE PRINMA(X,BIO2,TOTAL,NI,NIN,INT)
      IMPLICIT REAL*8 (A-H,O-Z)
      INCLUDE 'blmdim.inc'
      INCLUDE 'ioblck.inc'
      INCLUDE 'phyt1.inc'
      INCLUDE 'phyt2.inc'
      DIMENSION X(MX)
!
!  PRINT MAXIMUM SOLUTION ON UNIT IOU(6)
!
      IF (IOFLAG .EQ. 1) CALL MORESC
      WRITE (IOU(6),10)
   10 FORMAT (12X,'******* MAXIMUM SOLUTION *******')
      WRITE(IOU(6),20)
   20 FORMAT (2X,'Species',34X,'Types',/,26X,'1',13X,'2',13X,'3',
     1        13X,'4')
      DO 60 J=1,NUECOG
      L1=IT2(J,1)
      L2=IT2(J,2)
      WRITE (IOU(6),50) GRNAME(J), (X(K+NUROWS),K=L1,L2)
   50 FORMAT (2X,A8,11X,4(F11.4,3X))
   60 CONTINUE
      BIOPOS = BIO2
      IF (BIOPOS .LT. 0.0) BIOPOS = 0.0
      WRITE (IOU(6),70) BIOPOS
   70 FORMAT (2X,'Total biomass',6X,F11.4,3X,'g/m3')
      WRITE (IOU(6),90) TOTAL
   90 FORMAT (2X,'Chlorophyll',8X,F11.4,3X,'mg/m3',/)
!
!  PRINT NUTRIENT CONCENTRATIONS
!
      WRITE (IOU(6),100)
  100 FORMAT (2X,'Nutrient',14X,'Total',9X,'Slacks')
      WRITE (IOU(6),110) (CSTRA(K),CONCEN(K),X(K),K=1,NUNUCO)
  110 FORMAT (6(2X,A8,11X,F11.4,3X,F11.4,/))
      WRITE(IOU(6),120) NI,NIN,INT
  120 FORMAT ('  Number of intervals:',I3,2X,'Infeasible:',I3,
     1        2X,'Maximum interval:',I3,//)
      IF (IOFLAG .EQ. 1) CALL MORESC
      RETURN
      END

max_stars_repo_name: liujiamingustc/phd

max_issues_repo_name: liujiamingustc/phd

max_forks_repo_name: liujiamingustc/phd


text:
      SUBROUTINE CMSFIL
C
C     THIS SUBROUTINE GENERATES THE WORKING SUBSTRUCTURE FILE AND
C     APPLIES ALL TRANSFORMATIONS
C
      EXTERNAL        RSHIFT,ANDF
      LOGICAL         TDAT,XTRAN,XCSTM,IAUTO
CWKBI ALPHA-OSF 9/94
      INTEGER         SCR3, SCSFIL
      INTEGER         AAA(2),ANDF,RSHIFT,BUFEX,OUTT,
     1                SCBDAT,BUF4,SCORE,TRN,SYM,COMBO,NAM(2),Z,SCMCON,
     2                SCCSTM,CSTMID,CGID,BUF2,BUF3,BUF1,ECPT1,TWOJM1
      DIMENSION       RZ(1),ECPT(4),DOFN(6),LIST(32),TG(3,3),TG6(6,6),
     1                TSAVE(6,6),TMAT(6,6),TC(3,3),TT(3,3),XX(3)
      COMMON /CMB001/ SCR1,SCR2,SCBDAT,SCSFIL,SCCONN,SCMCON,SCTOC,
     1                GEOM4,CASECC,SCCSTM,SCR3
      COMMON /CMB002/ BUF1,BUF2,BUF3,BUF4,BUF5,SCORE,LCORE,INTP,OUTT
      COMMON /CMB003/ COMBO(7,5),CONSET,IAUTO,TOLER,NPSUB,CONECT,TRAN,
     1                MCON,RESTCT(7,7),ISORT,ORIGIN(7,3),IPRINT
      COMMON /CMB004/ TDAT(6)
      COMMON /GTMATX/ LOC1,KLTRAN,TRN,TT6(6,6),TC6(6,6)
      COMMON /ZZZZZZ/ Z(1)
      EQUIVALENCE     (ECPT1,ECPT(1)), (RZ(1),Z(1))
      DATA    AAA   / 4H CMS, 4HFIL  /
      DATA    NHBGSS, NHCSTM,NHEQSS  / 4HBGSS,4HCSTM,4HEQSS /
C
      BUFEX  = LCORE - BUF1 + BUF2
      LCORE  = BUFEX - 1
      IF (LCORE .LT. 0) GO TO 620
      LLCO   = LCORE
      IOEFIL = 310
      IFILE  = SCBDAT
      CALL OPEN (*600,SCBDAT,Z(BUF4),0)
C
C     READ GTRAN DATA INTO OPEN CORE
C
      IF (TDAT(3) .OR. TDAT(6)) CALL SKPFIL (SCBDAT,1)
      IF (.NOT.TDAT(6)) GO TO 20
      KSGTRN = SCORE
      CALL READ (*610,*10,SCBDAT,Z(KSGTRN),LLCO,1,KLGTRN)
      GO TO 620
   10 LLCO   = LLCO - KLGTRN
      KFGTRN = KSGTRN + KLGTRN - 1
      GO TO 30
C
C     READ TRANS DATA INTO OPEN CORE
C
   20 KSTRAN = SCORE
      GO TO 40
   30 KSTRAN = KFGTRN + 1
   40 IF (TDAT(3) .OR. TDAT(6)) CALL SKPFIL (SCBDAT,1)
      IF (.NOT. TDAT(3)) GO TO 60
      CALL READ (*610,*50,SCBDAT,Z(KSTRAN),LLCO,1,KLTRAN)
      GO TO 620
   50 LOC1   = KSTRAN
      LLCO   = LLCO - KLTRAN
      XTRAN  = .TRUE.
      KFTRAN = KSTRAN + KLTRAN - 1
      GO TO 70
   60 LOC1   = 0
      XTRAN  = .FALSE.
   70 CALL CLOSE (SCBDAT,1)
      IFILE  = SCSFIL
      CALL OPEN (*600,SCSFIL,Z(BUF4),1)
      IFILE  = SCCSTM
      CALL OPEN (*600,SCCSTM,Z(BUF3),1)
      KKC    = 0
C
C     LOOP ON EACH PSEUDOSTRUCTURE
C
      IFILE  = SCR3
      CALL OPEN (*600,SCR3,Z(BUF1),1)
      IFILE  = IOEFIL
      CALL OPEN (*600,IOEFIL,Z(BUFEX),1)
      LLCOLD = LLCO
C
      DO 550 I = 1,NPSUB
      LLCO   = LLCOLD
      NAM(1) = COMBO(I,1)
      NAM(2) = COMBO(I,2)
      TRN    = COMBO(I,3)
      SYM    = COMBO(I,4)
      NCOMP  = COMBO(I,5)
C
C     READ BGSS FOR I-TH PSEUDOSTRUCTURE
C
      KSBGSS = KSTRAN
      IF (XTRAN) KSBGSS = KSTRAN + KLTRAN
      CALL SFETCH (NAM,NHBGSS,1,ITEST)
      NGRP   = 1
      CALL SJUMP (NGRP)
      CALL SUREAD (Z(KSBGSS),LLCO,KLBGSS,ITEST)
      IF (KLBGSS.EQ.LLCO .AND. ITEST.NE.2) GO TO 620
      LLCO   = LLCO - KLBGSS
      KFBGSS = KSBGSS + KLBGSS - 1
      NIP    = KLBGSS/4
C
C     READ CSTM FOR THIS PSEUDOSTRUCTURE
C
      CALL SFETCH (NAM,NHCSTM,1,ITEST)
      XCSTM  = .FALSE.
      LOC2   = 0
      IF (ITEST .EQ. 3) GO TO 80
      KSCSTM = KFBGSS + 1
      NGRP   = 1
      CALL SJUMP (NGRP)
      CALL SUREAD (Z(KSCSTM),LLCO,KLCSTM,ITEST)
      IF (KLCSTM.EQ.LLCO .AND. ITEST.NE.2) GO TO 620
      LLCO   = LLCO - KLCSTM
      LOC2   = KSCSTM
      XCSTM  = .TRUE.
      KFCSTM = KSCSTM + KLCSTM - 1
   80 CONTINUE
C
C     DEFINE OPEN CORE ARRAYS FOR NEW CID AND HPTR
C
      KSNCID = KFBGSS + 1
      IF (LOC2 .NE. 0) KSNCID = KFCSTM + 1
      KLNCID = NIP
      LLCO   = LLCO - KLNCID
      KFNCID = KSNCID + KLNCID - 1
C
      KSHPTR = KFNCID + 1
      KLHPTR = NIP
      LLCO   = LLCO - KLHPTR
      KFHPTR = KSHPTR + KLHPTR - 1
C
C     SET ARRAYS TO ZERO
C
      DO 90 J = KSNCID,KFNCID
      Z(J) = 0
   90 CONTINUE
      DO 100 J = KSHPTR,KFHPTR
      Z(J) = 0
  100 CONTINUE
C
C     GET THE TRANS AND SYMT MATRIX FOR THIS PSEUDOSTRUCTURE
C
      CALL GTMAT1 (SYM,TT)
C
C     TRANSFORM THE COORDINATES IN THE BGSS, NOTE THAT THE ORIGINS
C     FOR TRANSLATION ARE STORED IN ARRAY ORIGIN.
C
      IF (TRN+SYM .EQ. 0) GO TO 130
      DO 120 J = KSBGSS,KFBGSS,4
      IF (Z(J) .EQ. -1) GO TO 120
      CALL GMMATS (TT,3,3,0 ,RZ(J+1),3,1,0 ,XX)
      DO 110 JJ = 1,3
      RZ(J+JJ) = XX(JJ) + ORIGIN(I,JJ)
  110 CONTINUE
  120 CONTINUE
  130 CONTINUE
C
C     TRANSFORM DEGREES OF FREEDOM FOR EACH EQSS CONTAINED
C     IN THE PSEUDOSTRUCTURE.
C
      CALL WRITE (SCR3,TT6,36,1)
      NHMAT = 1
      CALL SFETCH (NAM,NHEQSS,1,ITEST)
      KNEQSS = KFHPTR + 1
      CALL SUREAD (Z(KNEQSS),4,KLEQSS,ITEST)
      CALL SUREAD (Z(KNEQSS),-1,KLEQSS,ITEST)
      LLCO  = LLCO - 2*NCOMP
      IFILE = SCMCON
      CALL OPEN (*600,SCMCON,Z(BUF2),1)
      DO 350 J = 1,NCOMP
      KSEQSS = KFHPTR + 2*NCOMP
      CALL SUREAD (Z(KSEQSS),LLCO,KLEQSS,ITEST)
      IF (KLEQSS .EQ. 0) GO TO 340
      IF (KLEQSS.EQ.LLCO .AND. ITEST.NE.2) GO TO 620
      KFEQSS = KSEQSS + KLEQSS - 1
C
C     LOOP ON EACH IP IN THE EQSS AND GENERATE TRANSFORMATION MATRIX
C
      DO 330 JJ = KSEQSS,KFEQSS,3
      IP    = Z(JJ+1)
      ICOMP = Z(JJ+2)
C
C     GET CSTM FOR THIS IP
C
      CSTMID = Z(KSBGSS+4*IP-4)
      ECPT1  = CSTMID
      IF (CSTMID .LT. 0) ECPT1 = 0
      DO 140 JDH = 1,3
      ECPT(JDH+1) = RZ(KSBGSS+4*IP-4+JDH)
  140 CONTINUE
      CALL GTMAT2 (LOC2,KLCSTM,ECPT,TC)
C
C     TEST FOR POSSIBLE GTRAN
C
      IGTRAN = 0
      IF (.NOT.TDAT(6)) GO TO 170
      CGID = 1000000*J + Z(JJ)
      DO 150 K = KSGTRN,KFGTRN,5
      IF (Z(K+3).EQ.CGID .AND. Z(K).EQ.I .AND. Z(K+1).EQ.J) GO TO 160
  150 CONTINUE
C
C     NO GTRAN
C
      GO TO 170
  160 CALL GTMAT3 (Z(K+4),TG,TG6,IKIND)
      IGTRAN = 1
      GO TO 180
  170 CALL GTMAT3 (-1,TG,TG6,IKIND)
C
C     ALL TRANSFORMATIONS HAVE BEEN FOUND, COMPUTE THE FINAL MATRIX TMAT
C
  180 CALL GMMATS (TG6  ,6,6,1,TT6,6,6,0,TSAVE)
      CALL GMMATS (TSAVE,6,6,0,TC6,6,6,0,TMAT )
C
C     DECODE DEGREES OF FREEDOM AND FORM VECTOR
C
      CALL DECODE (ICOMP,LIST,NDOF)
C
C     FIND NEW DEGREES OF FREEDOM AND UPDATE EQSS
C
      IF (CSTMID.NE.0 .AND. IGTRAN.EQ.0) GO TO 220
      DO 200 I1 = 1,6
      DOFN(I1) = 0.0
      DO 190 I2 = 1,NDOF
      L = LIST(I2) + 1
      IF (ABS(TMAT(L,I1)) .LT. 1.0E-4) GO TO 190
      DOFN(I1) = 1.0
      GO TO 200
  190 CONTINUE
  200 CONTINUE
      ICODE = 0
      DO 210 I1 = 1,6
      ICODE = ICODE + DOFN(I1)*2**(I1-1)
  210 CONTINUE
      GO TO 230
  220 ICODE   = ICOMP
  230 Z(JJ+2) = ICODE
C
C     WRITE IP,C ON SCRATCH TO COMPUTE NEW SIL,C
C
      CALL WRITE (SCMCON,IP,1,0)
      CALL WRITE (SCMCON,ICODE,1,0)
C
C     UPDATE CID NUMBERS
C
      IADD   = KSBGSS + 4*IP - 4
      IADD1  = KSNCID + IP - 1
      IKKIND = IKIND  + 1
      GO TO (240,240,250,250,240,240,260,260,290,290,290,290,270,270,
     1       280,280,290,290,290,290,290,290,290,290,290,290,290,290,
     2       280,280,280,280,280,280,280), IKKIND
  240 Z(IADD1) = Z(IADD)
      IF (Z(IADD) .EQ. -1) Z(IADD1) = -100000000
      IF (Z(IADD) .EQ. -2) Z(IADD1) = -200000000
      GO TO 290
C
C     COMMENTS FROM G.CHAN/UNISYS  9/92
C     250 AND 240 ARE IDENTICAL HERE. IS IT POSSIBLY AN ERROR HERE?
C
  250 Z(IADD1) = Z(IADD)
      IF (Z(IADD) .EQ. -1) Z(IADD1) = -100000000
      IF (Z(IADD) .EQ. -2) Z(IADD1) = -200000000
      GO TO 290
  260 Z(IADD1) = 0
      GO TO 290
  270 Z(IADD1) = -TRN
      GO TO 290
  280 Z(IADD1) = -Z(K+4)
  290 CONTINUE
C
C     SET POINTERS FOR H MATRIX
C
      ITIS = 0
      IADD2 = KSHPTR + IP - 1
      IF (CSTMID   .LT. 0) GO TO 300
      IF (Z(IADD2) .GT. 2) GO TO 330
      IF (IKKIND.EQ. 3 .OR. IKKIND.EQ.4 .OR. IKKIND.EQ.13 .OR.
     1    IKKIND.EQ.14) ITIS = 1
      IF (IKKIND.EQ. 1 .OR. IKKIND.EQ.2 .OR. IKKIND.EQ. 5 .OR.
     1    IKKIND.EQ. 6) ITIS = 2
      IF (ITIS - 1) 320,300,310
  300 Z(IADD2) = 0
      GO TO 330
  310 Z(IADD2) = 1
      GO TO 330
  320 CONTINUE
      NHMAT = NHMAT + 1
      Z(IADD2) = NHMAT
      CALL WRITE (SCR3,TMAT,36,1)
  330 CONTINUE
C
C     INSERT MULTIPLE IP CODE
C
      IF (NCOMP .NE. 1) CALL EQSCOD (KSEQSS,KLEQSS,Z(1))
C
C     WRITE EQSS ON FILE SCSFIL
C
  340 CALL WRITE (SCSFIL,Z(KSEQSS),KLEQSS,1)
      TWOJM1 = 2*(J-1)
      IF (ANDF(RSHIFT(IPRINT,19),1) .EQ. 1)
     1    CALL CMIWRT (1,NAM,Z(KNEQSS+TWOJM1),KSEQSS,KLEQSS,Z,Z)
  350 CONTINUE
      CALL EOF (SCR3)
      CALL CLOSE (SCMCON,1)
C
C     GENERATE NEW SIL,C LIST
C
      IFILE = SCMCON
      CALL OPEN (*600,SCMCON,Z(BUF2),0)
      CALL READ (*360,*360,SCMCON,Z(KSEQSS),LLCO,1,NNN)
      GO TO 620
  360 CALL SORT (0,0,2,1,Z(KSEQSS),NNN)
      KSEJ = KSEQSS + NNN - 1
      I1 = KSEQSS
      I2 = KSEQSS + 2
  370 IF (I2-KSEQSS .GE. NNN) GO TO 400
      IF (Z(I1)   .EQ. Z(I2)) GO TO 380
      I1 = I1 + 2
      I2 = I2 + 2
      GO TO 370
  380 DO 390 J = I2,KSEJ
  390 Z(J-2) = Z(J)
      KSEJ = KSEJ - 2
      NNN  = NNN  - 2
      GO TO 370
  400 CONTINUE
      Z(KSEQSS) = 1
      DO 410 J = 3,NNN,2
      JJ = J - 1
      ICODE = Z(KSEQSS+JJ-1)
      CALL DECODE (ICODE,LIST,NDOF)
      Z(KSEQSS+JJ) = Z(KSEQSS+JJ-2) + NDOF
  410 CONTINUE
      CALL WRITE  (SCSFIL,Z(KSEQSS),NNN,1)
      CALL SUREAD (Z(KSEQSS),LLCO,KLEQSS,ITEST)
      CALL WRITE  (IOEFIL,Z(KSEQSS),KLEQSS,1)
      CALL CLOSE  (SCMCON,1)
C
C     PRINT EQSS SIL LIST IF REQUESTED
C
      IF (ANDF(RSHIFT(IPRINT,19),1) .EQ. 1)
     1    CALL CMIWRT (8,NAM,0,KSEQSS,KLEQSS,Z,Z)
C
C     UPDATE CSTM NUMBERING SYSTEM
C     KKC IS TRANSFORMED SYSTEM COORD. ID
C
      IP  = 0
      DO 540 I6 = KSNCID,KFNCID
      IP  = IP + 1
      LOC = KSBGSS + 4*(IP-1)
      IF (Z(I6) .EQ. 100000000) GO TO 540
      IF (Z(I6)) 430,520,420
  420 I1  = KSCSTM
      I2  = KFCSTM
      GO TO 440
  430 IF (Z(I6) .EQ. -100000000) GO TO 530
      I1  = KSTRAN
      I2  = KFTRAN
  440 KKC = KKC + 1
      IF (IAUTO) GO TO 480
      IF (KKC .GT. 1) GO TO 460
      CALL PAGE1
      CALL PAGE2 (5)
      WRITE  (OUTT,450)
  450 FORMAT (//45X,'SUMMARY OF OVERALL SYSTEM COORDINATES', //36X,
     1      'PSEUDO STRUCTURE ID.   SYSTEM COORD.ID    USER COORD.ID',/)
  460 CALL PAGE2 (1)
      WRITE  (OUTT,470) I,KKC,Z(LOC)
  470 FORMAT (43X,I6,14X,I6,11X,I6)
  480 LOOK4 = Z(I6)
      DO 500  J6 = I1,I2,14
      IF (IABS(Z(I6)) .NE. Z(J6)) GO TO 500
      IF (Z(I6) .LE. 0) GO TO 490
      CALL GMMATS (TT,3,3,0,Z(J6+5),3,3,0,TC)
      CALL WRITE (SCCSTM,KKC,1,0)
      CALL WRITE (SCCSTM,Z(J6+1),4,0)
      CALL WRITE (SCCSTM,TC,9,0)
      GO TO 500
  490 CONTINUE
      CALL WRITE (SCCSTM,KKC,1,0)
      CALL WRITE (SCCSTM,Z(J6+1),13,0)
  500 CONTINUE
C
C     FIND OTHER CIDS THAT ARE THE SAME
C
      IIP = 0
      DO 510 J6 = KSNCID,KFNCID
      IIP = IIP + 1
      IF (Z(J6).NE. LOOK4) GO TO 510
      LOC = KSBGSS + 4*(IIP-1)
      Z(LOC) = KKC
      Z(J6) = 100000000
  510 CONTINUE
      GO TO 540
  520 Z(LOC) = 0
      GO TO 540
  530 Z(LOC) = -1
      IF (Z(I6) .EQ. -200000000) Z(LOC) = -2
  540 CONTINUE
C
C     WRITE PROCESSED BGSS
C
      CALL WRITE (SCSFIL,Z(KSBGSS),KLBGSS,1)
      IF (ANDF(RSHIFT(IPRINT,18),1) .EQ.  1)
     1    CALL CMIWRT (2,NAM,NAM,KSBGSS,KLBGSS,Z(1),Z(1) )
C
C     WRITE ARRAY OF H POINTERS
C
      CALL WRITE (SCSFIL,Z(KSHPTR),KLHPTR,1)
      CALL EOF (SCSFIL)
  550 CONTINUE
C
      CALL CLOSE (SCR3,1)
      CALL CLOSE (SCSFIL,1)
      CALL WRITE (SCCSTM,TMAT,0,1)
      CALL CLOSE (SCCSTM,1)
      CALL CLOSE (IOEFIL,1)
      LCORE = BUFEX + BUF1 - BUF2
      RETURN
C
  600 IMSG = -1
      GO TO 630
  610 IMSG = -2
      GO TO 630
  620 IMSG = -8
  630 CALL MESAGE (IMSG,IFILE,AAA)
      RETURN
      END


text:
! For verifying parsing of various syntax
MODULE syntax

  IMPLICIT NONE

CONTAINS

  FUNCTION add_ints1(a,b) RESULT(y)
    INTEGER ( KIND = 1 ) , INTENT(in) :: a  ! Comment with & character
    INTEGER(1), INTENT(in) :: b
    INTEGER*1 :: y
    y = a+b
  END FUNCTION add_ints1

  FUNCTION add_ints2(a,b) RESULT(y)
    INTEGER*2, INTENT(in) :: a,& ! Line continuation with comment
      b
    INTEGER(KIND=2  ) :: y
    y = a+b
  END FUNCTION add_ints2

  FUNCTION add_ints4(a,b) RESULT(y)
    INTEGER (KIND  =4), INTENT(in) :: a  ! Comment with colons
    INTEGER, INTENT(in) :: b
    INTEGER :: y
    y = a+b
  END FUNCTION add_ints4

  ! Verifies argument clashes in contained procedures, and also that we
  ! don't wrap contained procedures
  FUNCTION contains_arg_clash(a,b) RESULT(y)
    INTEGER, INTENT(in) :: a, b
    INTEGER :: y

    y = a + b

  CONTAINS

    SUBROUTINE foo(a)
      REAL*8, INTENT(in) :: a
    END SUBROUTINE foo

    SUBROUTINE bar(b)
      REAL*8, INTENT(in) :: b
    END SUBROUTINE bar
  END FUNCTION contains_arg_clash

  FUNCTION argument_case_sensitivity(x) RESULT(Y)
    INTEGER, INTENT(in) :: X
    INTEGER :: y
    y = x
  END FUNCTION argument_case_sensitivity

END MODULE syntax

max_stars_repo_name: mcfarljm/fortwrap

max_issues_repo_name: mcfarljm/fortwrap

max_forks_repo_name: mcfarljm/fortwrap


text:
C  ****   New calculation of average temperature (TAVG)   **********
C  ****   in freezing/thawing layer using UP, DOWN, and MIDDLE   ***
c  ****   layer temperatures (TUP, TDN, TM)               **********
c!!!  Ken's version and early my version was wrong because all IF
c     statements compared to 0. not T0 temperature. However, results were
c     reasonably good due to a reasonable assumption that TAVG in all
c     combinations of TUP, TM & TDN equals (TUP+2*TM+TDN)/4.
c     In this version I replaced all IF statements to compare to T0,
c     and corrected one old statement (marked c!!!) that was wrong.

      FUNCTION ST_AVG1(TUP,TM,TDN,DZ)

      PARAMETER (T0=273.16)

      DZH=DZ*0.5
      IF(TUP .LT. T0) THEN

        IF(TM .LT. T0) THEN

          IF(TDN .LT. T0) THEN
C1  ******   TUP, TM, TDN < 0.  CASE  ************
            TAVG=(TUP+2*TM+TDN)/4.
            GOTO 777
C       C111  ****   END TDN_111  THEN

          ELSE
C2  ******   TUP & TM < 0.,  TDN > 0. CASE   *****
            X0=(T0-TM)*DZH/(TDN-TM)
            TAVG=0.5*(TUP*DZH+TM*(DZH+X0)+T0*(2.*DZH-X0))/DZ
            GOTO 777
C       C112   ****  END TDN_111  ELSE   ****
          ENDIF

C   C11  ***********  END TM_1 THEN   ***********
        ELSE

          IF(TDN .LT. T0) THEN
C3  *******  TUP < 0.  TM > 0.  TDN < 0. CASE   ****
            XUP=(T0-TUP)*DZH/(TM-TUP)
            XDN=DZH-(T0-TM)*DZH/(TDN-TM)
            TAVG=0.5*(TUP*XUP+T0*(2.*DZ-XUP-XDN)+TDN*XDN)/DZ
            GOTO 777
C       C121   ****   END TDN_121  THEN  *****

          ELSE
C4   ******  TUP < 0  TM > 0  TDN > 0  CASE   ******
            XUP=(T0-TUP)*DZH/(TM-TUP)
            TAVG=0.5*(TUP*XUP+T0*(2.*DZ-XUP))/DZ
            GOTO 777

C       C122   ***   END TDN_121  ELSE   ***
          ENDIF

C   C12  ***********  END TM_1 ELSE   **********
        ENDIF

C1    ********************  END TUP THEN   ***********
      ELSE

        IF(TM .LT. T0) THEN

          IF(TDN .LT. T0) THEN
C5   *****  TUP > 0  TM < 0 TDN < 0 CASE   *********
            XUP=DZH-(T0-TUP)*DZH/(TM-TUP)
            TAVG=0.5*(T0*(DZ-XUP)+TM*(DZH+XUP)+TDN*DZH)/DZ
            GOTO 777
C       C211   ****   END TDN_211  THEN  *****

          ELSE
C6   *****  TUP > 0  TM < 0  TDN > 0 CASE  *********
            XUP=DZH-(T0-TUP)*DZH/(TM-TUP)
            XDN=(T0-TM)*DZH/(TDN-TM)
c!!!  It was wrong statement
c!!!            TAVG=0.5*(T0*(2.*DZ-XUP+XDN)+TM*(XUP+XDN))/DZ
            TAVG=0.5*(T0*(2.*DZ-XUP-XDN)+TM*(XUP+XDN))/DZ
            GOTO 777

C       C212   ***   END TDN_211  ELSE   ***
          ENDIF

C   C21   ************   END TM_2  THEN   ************
        ELSE

          IF(TDN .LT. T0) THEN
C7   *****  TUP > 0  TM > 0  TDN < 0 CASE   ********
            XDN=DZH-(T0-TM)*DZH/(TDN-TM)
            TAVG=(T0*(DZ-XDN)+0.5*(T0+TDN)*XDN)/DZ
            GOTO 777
C       C221   ****   END TDN_221  THEN  *****

          ELSE
C8   *****  TUP > 0  TM > 0  TDN > 0 CASE   ********
            TAVG=(TUP+2.*TM+TDN)/4.
            GOTO 777
C       C222   ***   END TDN_221  ELSE   ***
          ENDIF

C   C22   ************   END TM_2  ELSE   ************
        ENDIF

C   *********************  END TUP ELSE  **************
      ENDIF
777   ST_AVG1=TAVG

      RETURN
      END

max_stars_repo_name: KathyNie/LISF

max_issues_repo_name: dmocko/LISF

max_forks_repo_name: dmocko/LISF


text:
      SUBROUTINE DRWINT (NUMBUN, BUNDIS, BUNDIF)

C  DRWINT draws a set of rectangles where the actual interiors
C  are drawn using the specified bundle values, and a set of
C  expected interiors are drawn using the attribute values
C  returned by <inquire predefined interior representation>.
C  It draws the expected interior incorrectly at the specified
C  position.
C
C  Input parameters:
C    NUMBUN : number of entries in BUNDIS
C    BUNDIS : list of interior bundle indices to be displayed
C    BUNDIF : position of interior to be drawn incorrectly

      COMMON /DIALOG/ DOUTYP, DINTYP, DSTDNR, DSTRID, PSTRID, DTCLIM,
     1                SCRMOD, DTXCI, SPECWT,
     2                DSIZE, EFRAC, DYXRAT, SYXRAT, MTRPDC, WCPDC, QVIS
      INTEGER         DOUTYP, DINTYP, DSTDNR, DSTRID, PSTRID, DTCLIM,
     1                SCRMOD, DTXCI, SPECWT
      REAL            DSIZE, EFRAC, DYXRAT, SYXRAT, MTRPDC, WCPDC, QVIS

C aspect source
C                bundled     individual
      INTEGER    PBUNDL,     PINDIV
      PARAMETER (PBUNDL = 0, PINDIV = 1)
C interior style
      INTEGER    PHOLLO,   PSOLID,   PPATTR,   PHATCH,   PISEMP
      PARAMETER (PHOLLO=0, PSOLID=1, PPATTR=2, PHATCH=3, PISEMP=4)

      INTEGER  NUMBUN, BUNDIS(NUMBUN), ERRIND, PDIS, PDISI,
     1         PDCI, BUNDIF, IX, SIZ, NPTS(1)

      REAL     XA(5),YA(5), YINCR,YLOC

C  use bundled attributes, except for edge attributes
      CALL SETASF (PBUNDL)
C  to ensure edges predictable - set ASF individual:
C  set edge flag (pedfg=14) = off
      CALL PSIASF (14, PINDIV)
      CALL PSEDFG (0)

C  draw and label actual
      CALL SETRVS ('0.2, 0.2, 0.45, 0.45', XA, SIZ)
      YINCR = 0.8/NUMBUN
      YLOC  = 0.9
      CALL NUMLAB (NUMBUN, 0.15, YLOC, YINCR)
      NPTS(1) = 4
      DO 100 IX = 1,NUMBUN
         YA(1) = YLOC - 0.3*YINCR
         YA(2) = YLOC + 0.3*YINCR
         YA(3) = YA(2)
         YA(4) = YA(1)
         YA(5) = YA(1)
         CALL PSII (BUNDIS(IX))
C switch between fill area and fill area set
         IF (MOD(IX,2) .EQ. 0) THEN
            CALL PFAS (1, NPTS, XA, YA)
         ELSE
            CALL PFA (4, XA, YA)
         ENDIF
         YLOC  = YLOC-YINCR
100   CONTINUE

C  now draw expected results - with individual attributes from WDT
C  use individual attributes
      CALL SETASF (PINDIV)
      CALL SETRVS ('0.6, 0.6, 0.85, 0.85', XA, SIZ)
      YLOC  = .9
      DO 200 IX = 1, NUMBUN
         CALL PQPIR (SPECWT, BUNDIS(IX), ERRIND, PDIS, PDISI, PDCI)
         CALL CHKINQ ('pqpir', ERRIND)
         IF (IX .EQ. BUNDIF) THEN
            IF (PDIS .EQ. PHOLLO) THEN
               PDIS = PISEMP
            ELSE
               PDIS = MOD(PDIS,4) + 1
            ENDIF
            PDISI = MOD(PDISI,4) +1
            PDCI = MOD(PDCI,5) + 1
         ENDIF
         CALL PSIS   (PDIS)
         CALL PSISI (PDISI)
         CALL PSICI (PDCI)
         YA(1) = YLOC - 0.3*YINCR
         YA(2) = YLOC + 0.3*YINCR
         YA(3) = YA(2)
         YA(4) = YA(1)
C switch between fill area and fill area set
         IF (MOD(IX, 2) .EQ. 0) THEN
            CALL PFAS (1, NPTS, XA, YA)
         ELSE
            CALL PFA (4, XA, YA)
         ENDIF
         YLOC  = YLOC-YINCR
200   CONTINUE

      END


text:
Kevin Massoudi is a current ASUCD ASUCD Senate Senator. Kevin has served on the Internal Affairs Commission since 2008 and has been active in campus affairs since his freshman year. He has also served on the Telecommunications Advisory Board and as a senate intern. He is the Social Chair of the Tau Kappa Epsilon fraternity.

Platform
1) Support The Green Initiative Fund (TGIF), and promote a sustainable campus by having ASUCD act as an example.
  (i.e. work to accomplish a nearpaperless ASUCD)

2) Work with taps to give more incentives for students who carpool.
  (Lower rates, better parking spots)

3)Work to expand Wireless internet through continued work with CCFIT and the Telecommunications Advisory board.

4)Bring back the epic quad battle!

Candidates Statement:

Hello fellow UC Davis students! My name is Kevin Massoudi and I am running for ASUCD Senate as an Independent candidate. Since arriving here at Davis, I have served the campus and students to the best of my ability, working both within and outside ASUCD to better student life as well as give presence to student voices during administrative discussions. A vote for me will be a vote for experience and dedication, with no constraining party line ideals.

As a Senator I will have four main goals that will be accomplished. First, I want to work to reduce ASUCD’s paper waste dramatically. By reversing a ban of electronics at Senate meetings, we can allow Senators and others sitting in at meetings to opt for an electronic copy of documents being circulated, saving precious resources and money for the students.

Second, in line with my commitment to reducing UC Davis’s environmental impact, I want to work with TAPS to eliminate the restrictions placed on carpool permits.  At this time, only students living outside of Davis city limits are allowed to apply for the discounted carpool permits being offered. By lifting this ban students stand to save $300 a year, and a ridesharing community will be encouraged.

Third, it is my goal to have wireless internet available throughout the entirety of the UC Davis central campus. We are a leading national research university and the necessary educational resources we provide our students are essential to the continued success of the university.

Finally, I will work with the administration to bring back one of the unique aspects of student life that went missing just last year. I speak of course, of the EPIC QUAD BATTLE!

Thank you for your consideration. Please vote Independent KEVIN MASSOUDI #1 for ASUCD Senate!

For more information on his platform and experience please visit http://www.kevinmassoudi.com

20090116 00:13:35 nbsp Kevin is a standup guy that everybody thinks is awesome. He has real experience, I know this because hes been serving with me on the Internal Affairs Commission for the past year. Not only does he hold the position, he uses it with the authority of change. Kevin was one of the loudest voices against closeddoor student government policies. He also works hard outside of his commission work to make positive change for students. Hes made important strides in working towards increased wireless coverage for UCD. A lot of candidates are going to make promises and will try to make up experience. Kevins been doing stuff and has experience that is beyond titles and positions.

Theres a reason Im backing Kevin Massoudi  the students. Users/GregWebb

20090218 23:52:25 nbsp Good luck covering the insurance on EQB I remember when the student was injured, and I also remember and agreed with the Administrations decision to prohibit it, on the basis that there was no organization or accountability. Users/MaxMikalonis

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


MODULE dtatsd
   !!======================================================================
   !!                     ***  MODULE  dtatsd  ***
   !! Ocean data  :  read ocean Temperature & Salinity Data from gridded data
   !!======================================================================
   !! History :  OPA  ! 1991-03  ()  Original code
   !!             -   ! 1992-07  (M. Imbard)
   !!            8.0  ! 1999-10  (M.A. Foujols, M. Imbard)  NetCDF FORMAT
   !!   NEMO     1.0  ! 2002-06  (G. Madec)  F90: Free form and module
   !!            3.3  ! 2010-10  (C. Bricaud, S. Masson)  use of fldread
   !!            3.4  ! 2010-11  (G. Madec, C. Ethe) Merge of dtatem and dtasal + suppression of CPP keys
   !!----------------------------------------------------------------------

   !!----------------------------------------------------------------------
   !!   dta_tsd      : read and time interpolated ocean Temperature & Salinity Data
   !!----------------------------------------------------------------------
   USE oce             ! ocean dynamics and tracers
   USE dom_oce         ! ocean space and time domain
   USE fldread         ! read input fields
   USE in_out_manager  ! I/O manager
   USE phycst          ! physical constants
   USE lib_mpp         ! MPP library
   USE wrk_nemo        ! Memory allocation
   USE timing          ! Timing

   IMPLICIT NONE
   PRIVATE

   PUBLIC   dta_tsd_init   ! called by opa.F90
   PUBLIC   dta_tsd        ! called by istate.F90 and tradmp.90

   LOGICAL , PUBLIC ::   ln_tsd_init      !: T & S data flag
   LOGICAL , PUBLIC ::   ln_tsd_tradmp    !: internal damping toward input data flag

   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_tsd   ! structure of input SST (file informations, fields read)

   !!----------------------------------------------------------------------
   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
   !! $Id: dtatsd.F90 6140 2015-12-21 11:35:23Z timgraham $
   !! Software governed by the CeCILL licence     (./LICENSE)
   !!----------------------------------------------------------------------
CONTAINS

   SUBROUTINE dta_tsd_init( ld_tradmp )
      !!----------------------------------------------------------------------
      !!                   ***  ROUTINE dta_tsd_init  ***
      !!
      !! ** Purpose :   initialisation of T & S input data
      !!
      !! ** Method  : - Read namtsd namelist
      !!              - allocates T & S data structure
      !!----------------------------------------------------------------------
      LOGICAL, INTENT(in), OPTIONAL ::   ld_tradmp   ! force the initialization when tradp is used
      !
      INTEGER ::   ierr0, ierr1, ierr2, ierr3   ! temporary integers
      !
      CHARACTER(len=100)            ::   cn_dir          ! Root directory for location of ssr files
      TYPE(FLD_N), DIMENSION( jpts) ::   slf_i           ! array of namelist informations on the fields to read
      TYPE(FLD_N)                   ::   sn_tem, sn_sal
      !!
      NAMELIST/namtsd/   ln_tsd_init, ln_tsd_tradmp, cn_dir, sn_tem, sn_sal
      INTEGER  ::   ios
      !!----------------------------------------------------------------------
      !
      IF( nn_timing == 1 )  CALL timing_start('dta_tsd_init')
      !
      !  Initialisation
      ierr0 = 0  ;  ierr1 = 0  ;  ierr2 = 0  ;  ierr3 = 0
      !
      REWIND( numnam_ref )              ! Namelist namtsd in reference namelist :
      READ  ( numnam_ref, namtsd, IOSTAT = ios, ERR = 901)
901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtsd in reference namelist', lwp )

      REWIND( numnam_cfg )              ! Namelist namtsd in configuration namelist : Parameters of the run
      READ  ( numnam_cfg, namtsd, IOSTAT = ios, ERR = 902 )
902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtsd in configuration namelist', lwp )
      IF(lwm) WRITE ( numond, namtsd )

      IF( PRESENT( ld_tradmp ) )   ln_tsd_tradmp = .TRUE.     ! forces the initialization when tradmp is used

      IF(lwp) THEN                  ! control print
         WRITE(numout,*)
         WRITE(numout,*) 'dta_tsd_init : Temperature & Salinity data '
         WRITE(numout,*) '~~~~~~~~~~~~ '
         WRITE(numout,*) '   Namelist namtsd'
         WRITE(numout,*) '      Initialisation of ocean T & S with T &S input data   ln_tsd_init   = ', ln_tsd_init
         WRITE(numout,*) '      damping of ocean T & S toward T &S input data        ln_tsd_tradmp = ', ln_tsd_tradmp
         WRITE(numout,*)
         IF( .NOT.ln_tsd_init .AND. .NOT.ln_tsd_tradmp ) THEN
            WRITE(numout,*)
            WRITE(numout,*) '   T & S data not used'
         ENDIF
      ENDIF
      !
      IF( ln_rstart .AND. ln_tsd_init ) THEN
         CALL ctl_warn( 'dta_tsd_init: ocean restart and T & S data intialisation, ',   &
            &           'we keep the restart T & S values and set ln_tsd_init to FALSE' )
         ln_tsd_init = .FALSE.
      ENDIF
      !
      !                             ! allocate the arrays (if necessary)
      IF(  ln_tsd_init .OR. ln_tsd_tradmp  ) THEN
         !
         ALLOCATE( sf_tsd(jpts), STAT=ierr0 )
         IF( ierr0 > 0 ) THEN
            CALL ctl_stop( 'dta_tsd_init: unable to allocate sf_tsd structure' )   ;   RETURN
         ENDIF
         !
                                ALLOCATE( sf_tsd(jp_tem)%fnow(jpi,jpj,jpk)   , STAT=ierr0 )
         IF( sn_tem%ln_tint )   ALLOCATE( sf_tsd(jp_tem)%fdta(jpi,jpj,jpk,2) , STAT=ierr1 )
                                ALLOCATE( sf_tsd(jp_sal)%fnow(jpi,jpj,jpk)   , STAT=ierr2 )
         IF( sn_sal%ln_tint )   ALLOCATE( sf_tsd(jp_sal)%fdta(jpi,jpj,jpk,2) , STAT=ierr3 )
         !
         IF( ierr0 + ierr1 + ierr2 + ierr3 > 0 ) THEN
            CALL ctl_stop( 'dta_tsd : unable to allocate T & S data arrays' )   ;   RETURN
         ENDIF
         !                         ! fill sf_tsd with sn_tem & sn_sal and control print
         slf_i(jp_tem) = sn_tem   ;   slf_i(jp_sal) = sn_sal
         CALL fld_fill( sf_tsd, slf_i, cn_dir, 'dta_tsd', 'Temperature & Salinity data', 'namtsd' )
         !
      ENDIF
      !
      IF( nn_timing == 1 )  CALL timing_stop('dta_tsd_init')
      !
   END SUBROUTINE dta_tsd_init

   SUBROUTINE dta_tsd( kt, ptsd )
      !!----------------------------------------------------------------------
      !!                   ***  ROUTINE dta_tsd  ***
      !!
      !! ** Purpose :   provides T and S data at kt
      !!
      !! ** Method  : - call fldread routine
      !!              - ORCA_R2: add some hand made alteration to read data
      !!              - 'key_orca_lev10' interpolates on 10 times more levels
      !!              - s- or mixed z-s coordinate: vertical interpolation on model mesh
      !!              - ln_tsd_tradmp=F: deallocates the T-S data structure
      !!                as T-S data are no are used
      !!
      !! ** Action  :   ptsd   T-S data on medl mesh and interpolated at time-step kt
      !!----------------------------------------------------------------------
      INTEGER                              , INTENT(in   ) ::   kt     ! ocean time-step
      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(  out) ::   ptsd   ! T & S data
      !
      INTEGER ::   ji, jj, jk, jl, jkk   ! dummy loop indicies
      INTEGER ::   ik, il0, il1, ii0, ii1, ij0, ij1   ! local integers
      REAL(wp)::   zl, zi
      REAL(wp), POINTER, DIMENSION(:) ::  ztp, zsp   ! 1D workspace
      !!----------------------------------------------------------------------
      !
      IF( nn_timing == 1 )  CALL timing_start('dta_tsd')
      !
      CALL fld_read( kt, 1, sf_tsd )      !==   read T & S data at kt time step   ==!
      !
      !
      !                                   !==   ORCA_R2 configuration and T & S damping   ==!
      IF( cp_cfg == "orca" .AND. jp_cfg == 2 .AND. ln_tsd_tradmp ) THEN    ! some hand made alterations
         !
         ij0 = 101   ;   ij1 = 109                       ! Reduced T & S in the Alboran Sea
         ii0 = 141   ;   ii1 = 155
         DO jj = mj0(ij0), mj1(ij1)
            DO ji = mi0(ii0), mi1(ii1)
               sf_tsd(jp_tem)%fnow(ji,jj,13:13) = sf_tsd(jp_tem)%fnow(ji,jj,13:13) - 0.20_wp
               sf_tsd(jp_tem)%fnow(ji,jj,14:15) = sf_tsd(jp_tem)%fnow(ji,jj,14:15) - 0.35_wp
               sf_tsd(jp_tem)%fnow(ji,jj,16:25) = sf_tsd(jp_tem)%fnow(ji,jj,16:25) - 0.40_wp
               !
               sf_tsd(jp_sal)%fnow(ji,jj,13:13) = sf_tsd(jp_sal)%fnow(ji,jj,13:13) - 0.15_wp
               sf_tsd(jp_sal)%fnow(ji,jj,14:15) = sf_tsd(jp_sal)%fnow(ji,jj,14:15) - 0.25_wp
               sf_tsd(jp_sal)%fnow(ji,jj,16:17) = sf_tsd(jp_sal)%fnow(ji,jj,16:17) - 0.30_wp
               sf_tsd(jp_sal)%fnow(ji,jj,18:25) = sf_tsd(jp_sal)%fnow(ji,jj,18:25) - 0.35_wp
            END DO
         END DO
         ij0 =  87   ;   ij1 =  96                          ! Reduced temperature in Red Sea
         ii0 = 148   ;   ii1 = 160
         sf_tsd(jp_tem)%fnow( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ,  4:10 ) = 7.0_wp
         sf_tsd(jp_tem)%fnow( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 11:13 ) = 6.5_wp
         sf_tsd(jp_tem)%fnow( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 14:20 ) = 6.0_wp
      ENDIF
      !
      ptsd(:,:,:,jp_tem) = sf_tsd(jp_tem)%fnow(:,:,:)    ! NO mask
      ptsd(:,:,:,jp_sal) = sf_tsd(jp_sal)%fnow(:,:,:)
      !
      IF( ln_sco ) THEN                   !==   s- or mixed s-zps-coordinate   ==!
         !
         CALL wrk_alloc( jpk, ztp, zsp )
         !
         IF( kt == nit000 .AND. lwp )THEN
            WRITE(numout,*)
            WRITE(numout,*) 'dta_tsd: interpolates T & S data onto the s- or mixed s-z-coordinate mesh'
         ENDIF
         !
         DO jj = 1, jpj                         ! vertical interpolation of T & S
            DO ji = 1, jpi
               DO jk = 1, jpk                        ! determines the intepolated T-S profiles at each (i,j) points
                  zl = gdept_0(ji,jj,jk)
                  IF(     zl < gdept_1d(1  ) ) THEN          ! above the first level of data
                     ztp(jk) =  ptsd(ji,jj,1    ,jp_tem)
                     zsp(jk) =  ptsd(ji,jj,1    ,jp_sal)
                  ELSEIF( zl > gdept_1d(jpk) ) THEN          ! below the last level of data
                     ztp(jk) =  ptsd(ji,jj,jpkm1,jp_tem)
                     zsp(jk) =  ptsd(ji,jj,jpkm1,jp_sal)
                  ELSE                                      ! inbetween : vertical interpolation between jkk & jkk+1
                     DO jkk = 1, jpkm1                                  ! when  gdept(jkk) < zl < gdept(jkk+1)
                        IF( (zl-gdept_1d(jkk)) * (zl-gdept_1d(jkk+1)) <= 0._wp ) THEN
                           zi = ( zl - gdept_1d(jkk) ) / (gdept_1d(jkk+1)-gdept_1d(jkk))
                           ztp(jk) = ptsd(ji,jj,jkk,jp_tem) + ( ptsd(ji,jj,jkk+1,jp_tem) - ptsd(ji,jj,jkk,jp_tem) ) * zi
                           zsp(jk) = ptsd(ji,jj,jkk,jp_sal) + ( ptsd(ji,jj,jkk+1,jp_sal) - ptsd(ji,jj,jkk,jp_sal) ) * zi
                        ENDIF
                     END DO
                  ENDIF
               END DO
               DO jk = 1, jpkm1
                  ptsd(ji,jj,jk,jp_tem) = ztp(jk) * tmask(ji,jj,jk)     ! mask required for mixed zps-s-coord
                  ptsd(ji,jj,jk,jp_sal) = zsp(jk) * tmask(ji,jj,jk)
               END DO
               ptsd(ji,jj,jpk,jp_tem) = 0._wp
               ptsd(ji,jj,jpk,jp_sal) = 0._wp
            END DO
         END DO
         !
         CALL wrk_dealloc( jpk, ztp, zsp )
         !
      ELSE                                !==   z- or zps- coordinate   ==!
         !
         ptsd(:,:,:,jp_tem) = ptsd(:,:,:,jp_tem) * tmask(:,:,:)    ! Mask
         ptsd(:,:,:,jp_sal) = ptsd(:,:,:,jp_sal) * tmask(:,:,:)
         !
         IF( ln_zps ) THEN                      ! zps-coordinate (partial steps) interpolation at the last ocean level
            DO jj = 1, jpj
               DO ji = 1, jpi
                  ik = mbkt(ji,jj)
                  IF( ik > 1 ) THEN
                     zl = ( gdept_1d(ik) - gdept_0(ji,jj,ik) ) / ( gdept_1d(ik) - gdept_1d(ik-1) )
                     ptsd(ji,jj,ik,jp_tem) = (1.-zl) * ptsd(ji,jj,ik,jp_tem) + zl * ptsd(ji,jj,ik-1,jp_tem)
                     ptsd(ji,jj,ik,jp_sal) = (1.-zl) * ptsd(ji,jj,ik,jp_sal) + zl * ptsd(ji,jj,ik-1,jp_sal)
                  ENDIF
                  ik = mikt(ji,jj)
                  IF( ik > 1 ) THEN
                     zl = ( gdept_0(ji,jj,ik) - gdept_1d(ik) ) / ( gdept_1d(ik+1) - gdept_1d(ik) )
                     ptsd(ji,jj,ik,jp_tem) = (1.-zl) * ptsd(ji,jj,ik,jp_tem) + zl * ptsd(ji,jj,ik+1,jp_tem)
                     ptsd(ji,jj,ik,jp_sal) = (1.-zl) * ptsd(ji,jj,ik,jp_sal) + zl * ptsd(ji,jj,ik+1,jp_sal)
                  END IF
               END DO
            END DO
         ENDIF
         !
      ENDIF
      !
      IF( .NOT.ln_tsd_tradmp ) THEN                   !==   deallocate T & S structure   ==!
         !                                              (data used only for initialisation)
         IF(lwp) WRITE(numout,*) 'dta_tsd: deallocte T & S arrays as they are only use to initialize the run'
                                        DEALLOCATE( sf_tsd(jp_tem)%fnow )     ! T arrays in the structure
         IF( sf_tsd(jp_tem)%ln_tint )   DEALLOCATE( sf_tsd(jp_tem)%fdta )
                                        DEALLOCATE( sf_tsd(jp_sal)%fnow )     ! S arrays in the structure
         IF( sf_tsd(jp_sal)%ln_tint )   DEALLOCATE( sf_tsd(jp_sal)%fdta )
                                        DEALLOCATE( sf_tsd              )     ! the structure itself
      ENDIF
      !
      IF( nn_timing == 1 )  CALL timing_stop('dta_tsd')
      !
   END SUBROUTINE dta_tsd

   !!======================================================================
END MODULE dtatsd


text:
c
c
c     ###################################################
c     ##  COPYRIGHT (C)  1990  by  Jay William Ponder  ##
c     ##              All Rights Reserved              ##
c     ###################################################
c
c     ##################################################################
c     ##                                                              ##
c     ##  subroutine lattice  --  setup periodic boundary conditions  ##
c     ##                                                              ##
c     ##################################################################
c
c
c     "lattice" stores the periodic box dimensions and sets angle
c     values to be used in computing fractional coordinates
c
c
      subroutine lattice
      implicit none
      include 'sizes.i'
      include 'boxes.i'
      include 'cell.i'
      include 'inform.i'
      include 'iounit.i'
      include 'math.i'
      real*8 alpha_cos
      real*8 ar1,ar2,ar3
      real*8 br1,br2,br3
      real*8 cr1,cr2,cr3
c
c
c     compute and store the half box length values
c
      xbox2 = 0.5d0 * xbox
      ybox2 = 0.5d0 * ybox
      zbox2 = 0.5d0 * zbox
      if (octahedron)  box34 = 0.75d0 * xbox
c
c     set replicated cell dimensions equal to the unitcell
c
      xcell = xbox
      ycell = ybox
      zcell = zbox
      xcell2 = xbox2
      ycell2 = ybox2
      zcell2 = zbox2
c
c     get values needed for fractional coordinate computations
c
      if (orthogonal .or. octahedron) then
         alpha_cos = 0.0d0
         beta_sin = 1.0d0
         beta_cos = 0.0d0
         gamma_sin = 1.0d0
         gamma_cos = 0.0d0
         beta_term = 0.0d0
         gamma_term = 1.0d0
      else if (monoclinic) then
         alpha_cos = 0.0d0
         beta_sin = sin(beta/radian)
         beta_cos = cos(beta/radian)
         gamma_sin = 1.0d0
         gamma_cos = 0.0d0
         beta_term = 0.0d0
         gamma_term = beta_sin
      else if (triclinic) then
         alpha_cos = cos(alpha/radian)
         beta_sin = sin(beta/radian)
         beta_cos = cos(beta/radian)
         gamma_sin = sin(gamma/radian)
         gamma_cos = cos(gamma/radian)
         beta_term = (alpha_cos - beta_cos*gamma_cos) / gamma_sin
         gamma_term = sqrt(beta_sin**2 - beta_term**2)
      end if
c
c     determine the volume of the parent periodic box
c
      volbox = 0.0d0
      if (orthogonal .or. octahedron) then
         volbox = xbox * ybox * zbox
      else if (monoclinic) then
         volbox = beta_sin * xbox * ybox * zbox
      else if (triclinic) then
         volbox = (gamma_sin*gamma_term) * xbox * ybox * zbox
      end if
c
c     compute and store real space lattice vectors as rows
c
      ar1 = xbox
      ar2 = 0.0d0
      ar3 = 0.0d0
      br1 = ybox * gamma_cos
      br2 = ybox * gamma_sin
      br3 = 0.0d0
      cr1 = zbox * beta_cos
      cr2 = zbox * beta_term
      cr3 = zbox * gamma_term
      lvec(1,1) = ar1
      lvec(1,2) = ar2
      lvec(1,3) = ar3
      lvec(2,1) = br1
      lvec(2,2) = br2
      lvec(2,3) = br3
      lvec(3,1) = cr1
      lvec(3,2) = cr2
      lvec(3,3) = cr3
c
c     compute and store reciprocal lattice vectors as columns
c
      if (volbox .ne. 0.0d0) then
         recip(1,1) = (br2*cr3 - cr2*br3) / volbox
         recip(2,1) = (br3*cr1 - cr3*br1) / volbox
         recip(3,1) = (br1*cr2 - cr1*br2) / volbox
         recip(1,2) = (cr2*ar3 - ar2*cr3) / volbox
         recip(2,2) = (cr3*ar1 - ar3*cr1) / volbox
         recip(3,2) = (cr1*ar2 - ar1*cr2) / volbox
         recip(1,3) = (ar2*br3 - br2*ar3) / volbox
         recip(2,3) = (ar3*br1 - br3*ar1) / volbox
         recip(3,3) = (ar1*br2 - br1*ar2) / volbox
      end if
c
c     volume of truncated octahedron is half of cubic parent
c
      if (octahedron)  volbox = 0.5d0 * volbox
      return
      end

max_stars_repo_name: htna/HCsbLib

max_issues_repo_name: htna/HCsbLib

max_forks_repo_name: htna/HCsbLib


text:
      SUBROUTINE DINTF (ANSWER, WORK, IOPT)
c Copyright (c) 1996 California Institute of Technology, Pasadena, CA.
c ALL RIGHTS RESERVED.
c Based on Government Sponsored Research NAS7-03001.
C>> 1996-04-26 DINTF  Krogh  Added ANSWER=0, for getting prototypes.
C>> 1994-10-20 DINTF  Krogh  Changes to use M77CON
C>> 1990-01-23 DINTF  CLL Added type stmts. Avoids msgs from compilers.
C>> 1988-06-07 DINTF  Snyder  Replaced " by '.
C>> 1987-08-14 DINTF  Snyder  Initial code.
c--D replaces "?": ?INTF, ?INT
C
C     Dummy DINTF subprogram that prints
C     'You have neglected either to write the DINTF subprogram'
C     'or to select reverse communication before invoking DINT.'
C     and then STOPs.
C
      double precision ANSWER, WORK(*)
      integer IOPT(*)
      ANSWER = 0.D0
      PRINT *,'You have neglected either to write the DINTF subprogram'
      PRINT *,'or to select reverse communication before invoking DINT*'
      STOP
      END


text:
C    Copyright(C) 1988-2017 National Technology & Engineering Solutions
C    of Sandia, LLC (NTESS).  Under the terms of Contract DE-NA0003525 with
C    NTESS, the U.S. Government retains certain rights in this software.
C
C    Redistribution and use in source and binary forms, with or without
C    modification, are permitted provided that the following conditions are
C    met:
C
C    * Redistributions of source code must retain the above copyright
C       notice, this list of conditions and the following disclaimer.
C
C    * Redistributions in binary form must reproduce the above
C      copyright notice, this list of conditions and the following
C      disclaimer in the documentation and/or other materials provided
C      with the distribution.
C
C    * Neither the name of NTESS nor the names of its
C      contributors may be used to endorse or promote products derived
C      from this software without specific prior written permission.
C
C    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
C    "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
C    LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
C    A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
C    OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
C    SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
C    LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
C    DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
C    THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
C    (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
C    OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

C $Id: elcent.f,v 1.1 1991/02/21 15:43:04 gdsjaar Exp $
C $Log: elcent.f,v $
C Revision 1.1  1991/02/21 15:43:04  gdsjaar
C Initial revision
C
      SUBROUTINE ELCENT (ELCEN, IX, COORD, NDIM, NUMEL, NELNOD, NUMNP)
      DIMENSION ELCEN(NUMEL, *), IX(NELNOD, *), COORD(NUMNP, *)
C
      IF (NDIM .EQ. 2) THEN
          DO 10 I=1, NUMEL
              ELCEN(I,1) = (COORD(IX(1,I),1) + COORD(IX(2,I),1) +
     *                      COORD(IX(3,I),1) + COORD(IX(4,I),1))/4.0
              ELCEN(I,2) = (COORD(IX(1,I),2) + COORD(IX(2,I),2) +
     *                      COORD(IX(3,I),2) + COORD(IX(4,I),2))/4.0
   10     CONTINUE
      ELSE
          DO 20 I=1, NUMEL
              ELCEN(I,1) = (COORD(IX(1,I),1) + COORD(IX(2,I),1) +
     *                      COORD(IX(3,I),1) + COORD(IX(4,I),1) +
     *                      COORD(IX(5,I),1) + COORD(IX(6,I),1) +
     *                      COORD(IX(7,I),1) + COORD(IX(8,I),1))/8.0
C
              ELCEN(I,2) = (COORD(IX(1,I),2) + COORD(IX(2,I),2) +
     *                      COORD(IX(3,I),2) + COORD(IX(4,I),2) +
     *                      COORD(IX(5,I),2) + COORD(IX(6,I),2) +
     *                      COORD(IX(7,I),2) + COORD(IX(8,I),2))/8.0
C
              ELCEN(I,3) = (COORD(IX(1,I),3) + COORD(IX(2,I),3) +
     *                      COORD(IX(3,I),3) + COORD(IX(4,I),3) +
     *                      COORD(IX(5,I),3) + COORD(IX(6,I),3) +
     *                      COORD(IX(7,I),3) + COORD(IX(8,I),3))/8.0
   20     CONTINUE
      END IF
      RETURN
      END

max_stars_repo_name: brian-kelley/seacas

max_issues_repo_name: brian-kelley/seacas

max_forks_repo_name: brian-kelley/seacas


module comm_br_mod
  implicit none

  ! *********************************************************************
  ! *  comm_br_mod -- An F90 module for computing the Blackwell-Rao     *
  ! *                 estimator given signal samples from the posterior *
  ! *                                                                   *
  ! *  H. K. Eriksen, E. Gjerløw, K. Mikkelsen (University of Oslo)     *
  ! *                                                                   *
  ! *                               and                                 *
  ! *                                                                   *
  ! *  I. K. Wehus, D. Pietrobon, K. M. Gorski, L. Colombo, G. Rocha,   *
  ! *  J. Bartlett, C. R. Lawrence, J. B. Jewell, S. Hildebrandt (JPL)  *
  ! *                                                                   *
  ! *                                                                   *
  ! *  Please cite the following papers when using this code:           *
  ! *                                                                   *
  ! *  - H. K. Eriksen et al. 2008, ApJ, 676, 10  (Commander)           *
  ! *  - M. Chu et al. 2005, Phys. Rev. D, 71, 103002 (Blackwell-Rao)   *
  ! *                                                                   *
  ! *********************************************************************

  ! ================================================================================================
  ! User routines:
  !
  !  ** subroutine comm_br_initialize_object(sigmafile, clfile, lmin, lmax, firstchain, lastchain, &
  !           & firstsample, lastsample, thinstep, handle)
  !
  !    Input parameters:
  !      sigmafile   (character=*) :: FITS file containing sigma_l samples in units of l(l+1)/2pi
  !      clfile      (character=*) :: ASCII file containing (l, C_l) in units of l(l+1)/2pi;
  !                                   used only for initialization (ie., offset estimation) to avoid
  !                                   numerical errors during exponentiation
  !      lmin        (i4b)         :: Lowest multipole to include in BR estimator
  !      lmax        (i4b)         :: Highest multipole to include in BR estimator
  !      firstchain  (i4b)         :: First chain to include in BR estimator
  !      lastchain   (i4b)         :: Last chain to include in BR estimator
  !      firstsample (i4b)         :: First sample to include in BR estimator; allow for burn-in if necessary
  !      lastsample  (i4b)         :: Last sample to include in BR estimator
  !      thinstep    (i4b)         :: Sample thinning factor -- typically 1
  !
  !    Output parameters:
  !      handle      (i4b)         :: this code supports multiple BR objects (optional; default=1), and
  !                                   the handle identifies specific BR data sets
  !
  !
  !  ** function comm_br_compute_lnL(cls, handle)
  !
  !    Input parameters:
  !      cls         (dp)          :: array containing at least cls(2:lmax) in units of l(l+1)/2pi
  !      handle      (i4b)         :: handle selecting which BR estimator to use (optional; default = 1)
  !
  !
  !  ** subroutine comm_br_deallocate_object(handle)
  !
  !    Input parameters:
  !      handle      (i4b)         :: handle of object to deallocate (optional; default=remove all)
  !
  ! ================================================================================================

  integer, parameter, private :: i4b = selected_int_kind(9)
  integer, parameter, private :: sp  = selected_real_kind(5,30)
  integer, parameter, private :: dp  = selected_real_kind(12,200)
  integer, parameter, private :: lgt = kind(.true.)

  type comm_br_data
     integer(i4b) :: lmin, lmax, numsamples, numchain
     integer(i4b) :: firstchain, lastchain, firstsample, lastsample, thinstep
     real(dp)     :: offset
     real(dp), allocatable, dimension(:,:,:) :: sigmas
  end type comm_br_data

  integer(i4b),       parameter       :: N_br = 100
  type(comm_br_data), dimension(N_br) :: comm_br

contains

  ! Initialization routines
  subroutine comm_br_initialize_object(sigmafile, clfile, lmin, lmax, firstchain, lastchain, &
       & firstsample, lastsample, thinstep, handle)
    implicit none

    integer(i4b),      intent(in)  :: lmin, lmax, firstchain, lastchain, firstsample, lastsample, thinstep
    character(len=*),  intent(in)  :: sigmafile, clfile
    integer(i4b),      intent(out), optional :: handle

    real(dp)           :: cl
    integer(i4b)       :: i, j, l, q, unit, numsamples, numchains, lmax_chain, id
    character(len=512) :: line
    real(dp), allocatable, dimension(:) :: cls

    id = 1
    if (present(handle)) then
       do while (allocated(comm_br(id)%sigmas) .and. id < N_br)
          id = id+1
       end do
       handle = id
    end if
    if (id < 1 .or. id > N_br) stop 'Error -- planck_br_mod: BR id out of range'

    ! Initialize BR object
    unit = comm_br_getlun()
    if (allocated(comm_br(id)%sigmas)) deallocate(comm_br(id)%sigmas)
    call comm_br_read_chain(sigmafile, unit, lmax_chain, numchains, numsamples, comm_br(id)%sigmas)
    if (lmax_chain < lmax .or. numchains < lastchain .or. numsamples < lastsample) then
       write(*,*) 'Error -- comm_br_mod: Requested parameters inconsistent with sigma file content'
       write(*,*) ''
       write(*,*) '    sigmafile  = ', trim(sigmafile)
       write(*,*) '    lmax_chain = ', lmax_chain, ' vs requested lmax       = ', lmax
       write(*,*) '    numchains  = ', numchains,  ' vs requested lastchain  = ', lastchain
       write(*,*) '    numsamples = ', numsamples, ' vs requested lastsample = ', lastsample
       stop
    end if
    comm_br(id)%lmin        = lmin
    comm_br(id)%lmax        = lmax
    comm_br(id)%firstsample = firstsample
    comm_br(id)%lastsample  = lastsample
    comm_br(id)%firstchain  = firstchain
    comm_br(id)%lastchain   = lastchain
    comm_br(id)%thinstep    = thinstep

    do j = comm_br(id)%firstsample, comm_br(id)%lastsample
       do i = comm_br(id)%firstchain, comm_br(id)%lastchain
          do l = comm_br(id)%lmin, comm_br(id)%lmax
             if (comm_br(id)%sigmas(l,i,j) <= 0.0) then
                write(*,*) 'Error -- comm_br_mod: Negative sigma_l = ', comm_br(id)%sigmas(l,i,j)
                write(*,*) '   l      = ', l
                write(*,*) '   chain  = ', i
                write(*,*) '   sample = ', j
                stop
             endif
          enddo
       enddo
    enddo

    ! Compute offset relative to input C_l spectrum
    allocate(cls(2:lmax))
    open(unit,file=trim(clfile))
    do while (.true.)
       read(unit,'(a)',end=99) line
       line = trim(line)
       if (line(1:1) == '#') cycle
       read(line,*) l, cl
       if (l >= lmin .and. l <= lmax) cls(l) = cl
    end do
99  close(unit)

    comm_br(id)%offset = comm_br_compute_offset(cls(2:lmax), handle)

  end subroutine comm_br_initialize_object

  subroutine comm_br_deallocate_object(handle)
    implicit none

    integer(i4b), optional :: handle

    integer(i4b) :: i, id_min, id_max

    id_min = 1; id_max = N_br;
    if (present(handle)) then
       id_min = handle
       id_max = handle
    end if

    do i = id_min, id_max
       if (allocated(comm_br(i)%sigmas)) deallocate(comm_br(i)%sigmas)
       comm_br(i)%lmin        = -1
       comm_br(i)%lmax        = -1
       comm_br(i)%firstchain  = -1
       comm_br(i)%lastchain   = -1
       comm_br(i)%firstsample = -1
       comm_br(i)%lastsample  = -1
       comm_br(i)%thinstep    =  0
    end do

  end subroutine comm_br_deallocate_object

  ! Base computation routine
  function comm_br_compute_lnL(cls, handle)
    implicit none

    real(dp),     dimension(2:), intent(in)           :: cls
    real(dp)                                          :: comm_br_compute_lnL
    integer(i4b),                intent(in), optional :: handle

    integer(i4b) :: i, j, l, id
    real(dp)     :: subtotal, x, lnL

    id = 1; if (present(handle)) id = handle
    if (.not. allocated(comm_br(id)%sigmas)) then
       write(*,*) 'Error -- comm_br_mod: Requested BR handle ', id, ' is not allocated'
       stop
    end if

    ! Compute the Blackwell-Rao estimator
    lnL = 0.d0
    do j = comm_br(id)%firstsample, comm_br(id)%lastsample, comm_br(id)%thinstep
       do i = comm_br(id)%firstchain, comm_br(id)%lastchain
          subtotal = 0.d0
          do l = comm_br(id)%lmin, comm_br(id)%lmax
             x = comm_br(id)%sigmas(l,i,j)/cls(l)
             subtotal = subtotal + &
                  & 0.5d0 * real(2*l+1,dp) * (-x + log(x)) - log(comm_br(id)%sigmas(l,i,j))
          end do
          lnL = lnL + exp(subtotal-comm_br(id)%offset)
       end do
    end do

    if (lnL > 1e-20) then
       lnL = log(lnL)
    else
       lnL = log(1e-30)
    end if

    comm_br_compute_lnL = lnL

  end function comm_br_compute_lnL

  ! Utility routine for initializing the offset to be subtracted from each term
  ! to avoid overflow errors. Only called with the first power spectrum
  function comm_br_compute_offset(cls, handle)
    implicit none

    real(dp),     dimension(2:), intent(in)           :: cls
    real(dp)                                          :: comm_br_compute_offset
    integer(i4b),                intent(in), optional :: handle

    integer(i4b) :: i, j, l, id
    real(dp)     :: subtotal, x, offset

    id = 1; if (present(handle)) id = handle

    ! Compute the Blackwell-Rao estimator
    offset = -1.6375d30
    do j = comm_br(id)%firstsample, comm_br(id)%lastsample, comm_br(id)%thinstep
       do i = comm_br(id)%firstchain, comm_br(id)%lastchain
          subtotal = 0.d0
          do l = comm_br(id)%lmin, comm_br(id)%lmax
             x = comm_br(id)%sigmas(l,i,j)/cls(l)
             subtotal = subtotal + &
                  & 0.5d0 * real(2*l+1,dp) * (-x + log(x)) - log(comm_br(id)%sigmas(l,i,j))
          end do
          offset = max(offset,subtotal)
       end do
    end do

    if (offset < -1.637e30) then
       write(*,*) 'Error -- comm_br_mod: Offset evaluation failed'
       write(*,*) '    handle      = ', id
       write(*,*) '    lmin        = ', comm_br(id)%lmin
       write(*,*) '    lmax        = ', comm_br(id)%lmax
       write(*,*) '    firstchain  = ', comm_br(id)%firstchain
       write(*,*) '    lastchain   = ', comm_br(id)%lastchain
       write(*,*) '    firstsample = ', comm_br(id)%firstsample
       write(*,*) '    lastsample  = ', comm_br(id)%lastsample
       write(*,*) '    thinstep    = ', comm_br(id)%thinstep
       write(*,*) ''
       write(*,*) '    cls         = ', cls(comm_br(id)%lmin:comm_br(id)%lmax)
       write(*,*) '    offset      = ', offset
       stop
    endif

    comm_br_compute_offset = offset

  end function comm_br_compute_offset

  function comm_br_getlun()
    implicit none
    integer(i4b) :: comm_br_getlun
    logical(lgt) :: exists, isopen
    comm_br_getlun = 9
    do
       comm_br_getlun = comm_br_getlun+1
       inquire(unit=comm_br_getlun,exist=exists)
       if(exists) then
          inquire(unit=comm_br_getlun,opened=isopen)
          if(.not. isopen) return
       end if
    end do
  end function comm_br_getlun

  ! Routine for reading the Gibbs sigma samples
  subroutine comm_br_read_chain(filename, unit, lmax, numchains, numsamples, data)
    implicit none

    character(len=*),                            intent(in)  :: filename
    integer(i4b),                                intent(in)  :: unit
    integer(i4b),                                intent(out) :: lmax, numchains, numsamples
    real(dp),         allocatable, dimension(:,:,:)          :: data

    integer(i4b)         :: l, status, blocksize, readwrite, numspec, i, j, k
    integer(i4b)         :: fpixel, group, numargs
    logical(lgt)         :: anyf
    real(dp)             :: nullval
    character(len=80)    :: comment

    integer(i4b),          dimension(4)     :: naxes
    real(dp),     pointer, dimension(:,:,:,:) :: indata

    status = 0
    readwrite = 0
    nullval = 0.

    ! numargs = 1
    numargs = 0

    ! Open the result file
    call ftopen(unit,trim(filename),readwrite,blocksize,status)

    ! Read keywords
    call ftgkyj(unit,'LMAX',     lmax,       comment,status)
    call ftgkyj(unit,'NUMSAMP',  numsamples, comment,status)
    call ftgkyj(unit,'NUMCHAIN', numchains,  comment,status)
    call ftgkyj(unit,'NUMSPEC',  numspec,    comment,status)

    allocate(data(0:lmax,numchains,numsamples))
    nullify(indata)
    allocate(indata(0:lmax,1:1,1:numchains,1:numargs+numsamples))

    ! Read the binned power spectrum array
    group  = 1
    fpixel = 1
    call ftgpvd(unit,group,fpixel,size(indata),nullval,indata,anyf,status)
    call ftclos(unit,status)

    do k = numargs+1, numargs+numsamples
       do j = 1, numchains
          do i = 0, lmax
             data(i, j, k) = indata(i, 1, j, k)
          enddo
       enddo
    enddo
    deallocate(indata)

  end subroutine comm_br_read_chain

end module comm_br_mod


subroutine calctends_no_tl(st,vp,t,p,mype,u_t,v_t,t_t,p_t,uvflag)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    calctends_tl       tlm of calctends
!   prgmmr: kleist           org: np20                date: 2005-09-29
!
! abstract: TLM of routine that compute tendencies for pressure, Tv, u, v
!
! program history log:
!   2005-09-29  kleist
!   2005-10-17  kleist - changes to improve computational efficiency
!   2005-11-21  kleist - add tracer tendencies, use new module
!   2006-04-12  treadon - replace sigi with bk5
!   2006-04-21  kleist - add divergence tendency parts
!   2006-07-31  kleist - changes to use ps instead of ln(ps)
!   2006-09-21  kleist - add rescaling to divergence tendency formulation
!   2006-10-04  rancic - correct bug in tracer advection terms
!   2007-04-16  kleist - move constraint specific items elsewhere
!   2007-05-08  kleist - add bits for fully generalized vertical coordinate
!   2007-06-21  rancic - add pbl
!   2007-07-26  cucurull - add 3d pressure pri in argument list ;
!                          move getprs_tl outside calctends_tl;
!                          call getprs_horiz_tl;
!                          remove ps from argument list
!   2007-08-08  derber - optimize
!   2008-06-05  safford - rm unused uses
!   2009-04-21  derber - remove call to getuv and modify get_derivatives to include uv
!   2009-08-20  parrish - replace curvfct with curvx, curvy.  this allows tendency computation to
!                          work for any general orthogonal coordinate.
!   2010-11-03  derber - moved threading calculations to gridmod and modified
!   2012-02-08  kleist - add uvflag to argument list
!   2013-10-19  todling - derivatives and guess fields now in bundles
!
! usage:
!   input argument list:
!     u        - zonal wind on subdomain
!     v        - meridional wind on subdomain
!     t        - virtual temperature on subdomain
!     mype     - task id
!     uvflag   - logical, set to true for st,vp wind components, instead of stream/potential function
!
!   output argument list:
!     u_t      - time tendency of u
!     v_t      - time tendency of v
!     t_t      - time tendency of t
!     p_t      - time tendency of 2d surface pressure
!
!   notes:
!     TLM check performed & verified on 2005-09-29 by d. kleist
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$
  use kinds,only: r_kind,i_kind
  use gridmod, only: lat2,lon2,nsig,istart,nlat,idvc5,bk5,&
      wrf_nmm_regional,nems_nmmb_regional,eta2_ll,regional,nthreads,jtstart,jtstop
  use constants, only: zero,half,two,rd,rcp
  use tendsmod, only: what9,prsth9,r_prsum9,prdif9,r_prdif9,pr_xsum9,pr_xdif9,&
      pr_ysum9,pr_ydif9,curvx,curvy,coriolis
  use guess_grids, only: ntguessig,ges_teta,ges_prsi
  use gsi_metguess_mod, only: gsi_metguess_bundle
  use gsi_bundlemod, only: gsi_bundlegetpointer
  use derivsmod, only: gsi_xderivative_bundle
  use derivsmod, only: gsi_yderivative_bundle
  implicit none

! Declare passed variables
  real(r_kind),dimension(lat2,lon2,nsig),intent(in   ) :: st,vp,t
  real(r_kind),dimension(lat2,lon2)     ,intent(in   ) :: p
  integer(i_kind)                       ,intent(in   ) :: mype
  real(r_kind),dimension(lat2,lon2,nsig),intent(  out) :: u_t,v_t,t_t
  real(r_kind),dimension(lat2,lon2)     ,intent(  out) :: p_t
  logical                               ,intent(   in) :: uvflag

! Declare local variables
  character(len=*),parameter::myname='calctends_no_tl'
  real(r_kind),dimension(lat2,lon2,nsig+1):: pri
  real(r_kind),dimension(lat2,lon2,nsig):: u,v
  real(r_kind),dimension(lat2,lon2,nsig+1):: pri_x,pri_y
  real(r_kind),dimension(lat2,lon2,nsig+1):: prsth,what
  real(r_kind),dimension(lat2,lon2,nsig):: prsum,prdif,pr_xsum,pr_xdif,&
       pr_ysum,pr_ydif
  real(r_kind),dimension(lat2,lon2):: sumkm1,sumvkm1,sum2km1,sum2vkm1
  real(r_kind),dimension(lat2,lon2,nsig):: t_thor9
  real(r_kind),dimension(lat2,lon2,nsig):: u_x,u_y,v_x,v_y,t_x,t_y

  real(r_kind) tmp,tmp2,tmp3,sumk,sumvk,sum2k,sum2vk,uduvdv
  integer(i_kind) i,j,k,ix,it,kk,ier,istatus

  real(r_kind),pointer,dimension(:,:,:) :: ges_u=>NULL()
  real(r_kind),pointer,dimension(:,:,:) :: ges_v=>NULL()
  real(r_kind),pointer,dimension(:,:,:) :: ges_tv=>NULL()

  real(r_kind),pointer,dimension(:,:,:) :: ges_u_lon=>NULL()
  real(r_kind),pointer,dimension(:,:,:) :: ges_v_lon=>NULL()
  real(r_kind),pointer,dimension(:,:,:) :: ges_tv_lon=>NULL()

  real(r_kind),pointer,dimension(:,:,:) :: ges_u_lat=>NULL()
  real(r_kind),pointer,dimension(:,:,:) :: ges_v_lat=>NULL()
  real(r_kind),pointer,dimension(:,:,:) :: ges_tv_lat=>NULL()

! linearized about guess solution, so set it flag accordingly
  it=ntguessig

  ier=0
  call gsi_bundlegetpointer(gsi_metguess_bundle(it),'u' ,ges_u, istatus);ier=istatus+ier
  call gsi_bundlegetpointer(gsi_metguess_bundle(it),'v' ,ges_v, istatus);ier=istatus+ier
  call gsi_bundlegetpointer(gsi_metguess_bundle(it),'tv',ges_tv,istatus);ier=istatus+ier
  if(ier/=0) then
     write(6,*) myname, ': pointers not found in met-guess, ier=', ier
     call stop2(999)
  endif

  ier=0
  call gsi_bundlegetpointer(gsi_xderivative_bundle(it),'u' ,ges_u_lon, istatus);ier=istatus+ier
  call gsi_bundlegetpointer(gsi_xderivative_bundle(it),'v' ,ges_v_lon, istatus);ier=istatus+ier
  call gsi_bundlegetpointer(gsi_xderivative_bundle(it),'tv',ges_tv_lon,istatus);ier=istatus+ier
  if(ier/=0) then
     write(6,*) myname, ': pointers not found in lon-derivatives, ier=', ier
     call stop2(999)
  endif

  ier=0
  call gsi_bundlegetpointer(gsi_yderivative_bundle(it),'u' ,ges_u_lat, istatus);ier=istatus+ier
  call gsi_bundlegetpointer(gsi_yderivative_bundle(it),'v' ,ges_v_lat, istatus);ier=istatus+ier
  call gsi_bundlegetpointer(gsi_yderivative_bundle(it),'tv',ges_tv_lat,istatus);ier=istatus+ier
  if(ier/=0) then
     write(6,*) myname, ': pointers not found in lat-derivatives, ier=', ier
     call stop2(999)
  endif

! preliminaries:

! get 3d pressure
  call getprs_tl(p,t,pri)

  if(uvflag)then
    call get_derivatives2uv(st,vp,t,pri,u,v,u_x,v_x,t_x,pri_x, &
                                        u_y,v_y,t_y,pri_y)
  else
    call get_derivatives2(st,vp,t,pri,u,v,u_x,v_x,t_x,pri_x, &
                                        u_y,v_y,t_y,pri_y)
  end if

!$omp parallel do private(i,j,k,kk,tmp,tmp2,uduvdv, &
!$omp                  tmp3,sumk,sumvk,sum2k,sum2vk,ix)
  do kk=1,nthreads

    do k=1,nsig
      do j=jtstart(kk),jtstop(kk)
        do i=1,lat2
           prsum  (i,j,k)=pri  (i,j,k)+pri  (i,j,k+1)
           prdif  (i,j,k)=pri  (i,j,k)-pri  (i,j,k+1)
           pr_xsum(i,j,k)=pri_x(i,j,k)+pri_x(i,j,k+1)
           pr_xdif(i,j,k)=pri_x(i,j,k)-pri_x(i,j,k+1)
           pr_ysum(i,j,k)=pri_y(i,j,k)+pri_y(i,j,k+1)
           pr_ydif(i,j,k)=pri_y(i,j,k)-pri_y(i,j,k+1)
        end do
      end do
    end do

!   Compute horizontal part of tendency for 3d pressure

    do j=jtstart(kk),jtstop(kk)
      do i=1,lat2
        prsth(i,j,nsig+1)=zero
      end do
    end do
    do k=nsig,1,-1
      do j=jtstart(kk),jtstop(kk)
        do i=1,lat2
           prsth(i,j,k)=prsth(i,j,k+1) -  ( u(i,j,k)*pr_xdif9(i,j,k) + &
              ges_u(i,j,k)*pr_xdif(i,j,k) + v(i,j,k)*pr_ydif9(i,j,k) + &
              ges_v(i,j,k)*pr_ydif(i,j,k) + &
              (u_x      (i,j,k) + v_y      (i,j,k))*(prdif9(i,j,k)) + &
              (ges_u_lon(i,j,k) + ges_v_lat(i,j,k))*(prdif (i,j,k)) )
        end do
      end do
    end do

!   Get horizontal part of temperature tendency for vertical velocity term

    do k=1,nsig
      do j=jtstart(kk),jtstop(kk)
        do i=1,lat2
           tmp=rd*ges_tv(i,j,k)*r_prsum9(i,j,k)
           t_t(i,j,k)=-u(i,j,k)*ges_tv_lon(i,j,k) - ges_u(i,j,k)*t_x(i,j,k) - &
              v(i,j,k)*ges_tv_lat(i,j,k) - ges_v(i,j,k)*t_y(i,j,k) +          &
              tmp*rcp*( ges_u(i,j,k)*pr_xsum(i,j,k) + &
              u(i,j,k)*pr_xsum9(i,j,k) + &
              ges_v(i,j,k)*pr_ysum(i,j,k) + &
              v(i,j,k)*pr_ysum9(i,j,k) + &
              prsth(i,j,k) + prsth(i,j,k+1)) + &
              rcp*( ges_u(i,j,k)*pr_xsum9(i,j,k) + &
              ges_v(i,j,k)*pr_ysum9(i,j,k) + &
              prsth9(i,j,k)+prsth9(i,j,k+1) ) * &
              ( rd*t(i,j,k)*r_prsum9(i,j,k) - tmp*prsum(i,j,k)*r_prsum9(i,j,k) )

        end do
      end do
    end do

!   calculate vertical velocity term:  z(dp/dz) (zero at top/bottom interfaces)
!   if running global, and there is a c(k) coefficient, we call the vvel subroutine

    if ( (.not.regional) .AND. (idvc5==3)) then

!     Get horizontal part of temperature tendency for vertical velocity term

      do k=1,nsig
        do j=jtstart(kk),jtstop(kk)
           do i=1,lat2
              tmp=rd*ges_tv(i,j,k)*r_prsum9(i,j,k)

              t_thor9(i,j,k)=-ges_u(i,j,k)*ges_tv_lon(i,j,k) - &
                   ges_v(i,j,k)*ges_tv_lat(i,j,k)
              t_thor9(i,j,k)=t_thor9(i,j,k) -tmp*rcp * ( ges_u(i,j,k)*pr_xsum9(i,j,k) + &
                   ges_v(i,j,k)*pr_ysum9(i,j,k) + &
                   prsth9(i,j,k) + prsth9(i,j,k+1) )
           end do
        end do
      end do
      call getvvel_tl(t,t_t,t_thor9,prsth,prdif,what,jtstart(kk),jtstop(kk))
    else
      if(wrf_nmm_regional.or.nems_nmmb_regional) then
        do k=2,nsig
           do j=jtstart(kk),jtstop(kk)
              do i=1,lat2
                 what(i,j,k)=prsth(i,j,k)-eta2_ll(k)*prsth(i,j,1)
              end do
           end do
        end do
      else
        do k=2,nsig
           do j=jtstart(kk),jtstop(kk)
              do i=1,lat2
                 what(i,j,k)=prsth(i,j,k)-bk5(k)*prsth(i,j,1)
              end do
           end do
        end do
      end if
    end if

!   top/bottom boundary condition:

    do j=jtstart(kk),jtstop(kk)
      do i=1,lat2
        what(i,j,1)=zero
        what(i,j,nsig+1)=zero
      enddo
    enddo

!   load actual dp/dt

    do j=jtstart(kk),jtstop(kk)
      do i=1,lat2
        p_t(i,j)=prsth(i,j,1)
      end do
    end do

!   before big k loop, zero out the km1 summation arrays

    do j=jtstart(kk),jtstop(kk)
      do i=1,lat2
        sumkm1  (i,j)=zero
        sum2km1 (i,j)=zero
        sumvkm1 (i,j)=zero
        sum2vkm1(i,j)=zero
      end do
    end do

!   Compute terms for tendencies of wind components & Temperature

    do k=1,nsig
      do j=jtstart(kk),jtstop(kk)
        do i=1,lat2
           uduvdv=two*(ges_u(i,j,k)*u(i,j,k) + ges_v(i,j,k)*v(i,j,k))
           u_t(i,j,k)=-u(i,j,k)*ges_u_lon(i,j,k) - ges_u(i,j,k)*u_x(i,j,k) - &
              v(i,j,k)*ges_u_lat(i,j,k) - ges_v(i,j,k)*u_y(i,j,k) + &
              coriolis(i,j)*v(i,j,k) + curvx(i,j)*uduvdv
           v_t(i,j,k)=-u(i,j,k)*ges_v_lon(i,j,k) - ges_u(i,j,k)*v_x(i,j,k) - &
              v(i,j,k)*ges_v_lat(i,j,k) - ges_v(i,j,k)*v_y(i,j,k) - &
              coriolis(i,j)*u(i,j,k) + curvy(i,j)*uduvdv

           tmp=rd*ges_tv(i,j,k)*r_prsum9(i,j,k)
           tmp2=rd*t(i,j,k)*r_prsum9(i,j,k)

           u_t(i,j,k) = u_t(i,j,k)-tmp*( pr_xsum(i,j,k) - &
              (prsum(i,j,k)*pr_xsum9(i,j,k)*r_prsum9(i,j,k)) ) - &
              tmp2*pr_xsum9(i,j,k)

           v_t(i,j,k) = v_t(i,j,k)-tmp*( pr_ysum(i,j,k) - &
              (prsum(i,j,k)*pr_ysum9(i,j,k)*r_prsum9(i,j,k)) ) - &
              tmp2*pr_ysum9(i,j,k)

!          vertical flux terms
           if (k > 1) then
              tmp=half*what(i,j,k)*r_prdif9(i,j,k)
              tmp2=half*what9(i,j,k)*r_prdif9(i,j,k)

              u_t(i,j,k) = u_t(i,j,k) - tmp*(ges_u (i,j,k-1)-ges_u (i,j,k)) - &
                 tmp2*( (u(i,j,k-1)-u(i,j,k)) - (ges_u (i,j,k-1)-ges_u (i,j,k))* &
                 prdif(i,j,k)*r_prdif9(i,j,k))
              v_t(i,j,k) = v_t(i,j,k) - tmp*(ges_v (i,j,k-1)-ges_v (i,j,k)) - &
                 tmp2*( (v(i,j,k-1)-v(i,j,k)) - (ges_v (i,j,k-1)-ges_v (i,j,k))* &
                 prdif(i,j,k)*r_prdif9(i,j,k))
              t_t(i,j,k) = t_t(i,j,k) - tmp*(ges_tv(i,j,k-1)-ges_tv(i,j,k)) - &
                 tmp2*( (t(i,j,k-1)-t(i,j,k)) - (ges_tv(i,j,k-1)-ges_tv(i,j,k))* &
                 prdif(i,j,k)*r_prdif9(i,j,k))
           end if
           if (k < nsig) then
              tmp=half*what(i,j,k+1)*r_prdif9(i,j,k)
              tmp2=half*what9(i,j,k+1)*r_prdif9(i,j,k)
              u_t(i,j,k) = u_t(i,j,k) - tmp*(ges_u (i,j,k)-ges_u (i,j,k+1)) - &
                 tmp2*( (u(i,j,k)-u(i,j,k+1)) - (ges_u (i,j,k)-ges_u (i,j,k+1))* &
                 prdif(i,j,k)*r_prdif9(i,j,k))
              v_t(i,j,k) = v_t(i,j,k) - tmp*(ges_v (i,j,k)-ges_v (i,j,k+1)) - &
                 tmp2*( (v(i,j,k)-v(i,j,k+1)) - (ges_v (i,j,k)-ges_v (i,j,k+1))* &
                 prdif(i,j,k)*r_prdif9(i,j,k))
              t_t(i,j,k) = t_t(i,j,k) - tmp*(ges_tv(i,j,k)-ges_tv(i,j,k+1)) - &
                 tmp2*( (t(i,j,k)-t(i,j,k+1)) - (ges_tv(i,j,k)-ges_tv(i,j,k+1))* &
                 prdif(i,j,k)*r_prdif9(i,j,k))
           end if
        end do   !end do i
      end do     !end do j

!     first sum to level k-1

      do j=jtstart(kk),jtstop(kk)
        do i=1,lat2
           tmp=rd*t(i,j,k)*r_prsum9(i,j,k)
           tmp2=rd*ges_tv(i,j,k)*r_prsum9(i,j,k)
           tmp3=prdif9(i,j,k)*r_prsum9(i,j,k)

           sumk = sumkm1(i,j) + ( tmp - tmp2*prsum(i,j,k)*r_prsum9(i,j,k) ) * &
              ( pr_xdif9(i,j,k) - ( pr_xsum9(i,j,k)*tmp3 ) )
           sumk = sumk + tmp2*( pr_xdif(i,j,k) - &
              (tmp3 *pr_xsum(i,j,k) + pr_xsum9(i,j,k)*( ( &
              prdif(i,j,k) - tmp3* prsum(i,j,k) )*r_prsum9(i,j,k) ) ) )

           sumvk = sumvkm1(i,j) + ( tmp - tmp2*prsum(i,j,k)*r_prsum9(i,j,k) ) * &
              ( pr_ydif9(i,j,k) - ( pr_ysum9(i,j,k)*tmp3 ) )
           sumvk = sumvk + tmp2*( pr_ydif(i,j,k) - &
              (tmp3*pr_ysum(i,j,k) + pr_ysum9(i,j,k)*( ( &
              prdif(i,j,k) - tmp3* prsum(i,j,k) )*r_prsum9(i,j,k) ) ) )

           sum2k = sum2km1(i,j) + t_x(i,j,k)*tmp3 + &
              ges_tv_lon(i,j,k)*( (prdif(i,j,k) - &
              tmp3*prsum(i,j,k))*r_prsum9(i,j,k))

           sum2vk = sum2vkm1(i,j) + t_y(i,j,k)*tmp3 + &
              ges_tv_lat(i,j,k)*( (prdif(i,j,k) - &
              tmp3*prsum(i,j,k))*r_prsum9(i,j,k))

           u_t(i,j,k) = u_t(i,j,k) - sumkm1 (i,j) - rd*sum2km1 (i,j) - &
              sumk  - rd*sum2k
           v_t(i,j,k) = v_t(i,j,k) - sumvkm1(i,j) - rd*sum2vkm1(i,j) - &
              sumvk - rd*sum2vk

!          load up the km1 arrays for next k loop
           sumkm1  (i,j)=sumk
           sumvkm1 (i,j)=sumvk
           sum2km1 (i,j)=sum2k
           sum2vkm1(i,j)=sum2vk

        end do !end do i
      end do   !end do j
    end do  !end do k

    call turbl_tl(ges_prsi(1,1,1,it),ges_tv,ges_teta(1,1,1,it),&
               u,v,pri,t,u_t,v_t,t_t,jtstart(kk),jtstop(kk))

    if(.not.wrf_nmm_regional.and..not.nems_nmmb_regional)then
      do k=1,nsig

!       Zero out time derivatives at poles

        do j=jtstart(kk),jtstop(kk)
           do i=1,lat2
              ix=istart(mype+1)+i-2
              if (ix == 1 .or. ix == nlat) then
                 u_t(i,j,k)=zero
                 v_t(i,j,k)=zero
              end if
           end do
        end do

      end do  !end do k
    end if

  end do
!  end threading loop

  return
end subroutine calctends_no_tl


text:
function plm(l,m,x)
!----------------------------------------------------
! Subroutine:
!   plm
!
! Purpose:
!  Calculate assiociated Legendre polynomials
!----------------------------------------------------
  implicit none
  integer,parameter::dp=selected_real_kind(8) !double precision
  integer::l,m
  integer::i,ll
  real(kind=dp)::x,fact,pll,pmm,pmmp1,somx2
  real(kind=dp)::plm
  IF((m.lt.0).or.(m.gt.l).or.(abs(x).gt.1.D0))THEN
      WRITE(*,*) 'bad arguments in plm'
  ENDIF
  pmm=1.D0
  if(m.gt.0)then
     somx2=sqrt((1.D0-x)*(1.D0+x))
     fact=1.D0
     do 11 i=1,m
        pmm=-pmm*fact*somx2
        fact=fact+2.D0
11   continue
  endif
  if(l.eq.m)then
     plm=pmm
  else
     pmmp1=x*dble(2*m+1)*pmm
     if(l.eq.m+1)then
        plm=pmmp1
     else
        do 12 ll=m+2,l
           pll=(x*dble(2*ll-1)*pmmp1-dble(ll+m-1)*pmm)/dble(ll-m)
           pmm=pmmp1
           pmmp1=pll
12      continue
        plm=pll
     endif
  endif
  return
end function plm

max_stars_repo_name: rutgervandeelen/das

max_issues_repo_name: rutgervandeelen/das

max_forks_repo_name: rutgervandeelen/das


SUBROUTINE getname(argv,nlen)
!
! This subroutine reads the base name of data file.
!
 IMPLICIT NONE
 INTEGER::narg
 INTEGER,INTENT(OUT)::nlen
 INTEGER::lnblnk,iargc
 CHARACTER(*),INTENT(OUT)::argv
 LOGICAL found
 narg=IARGC()
 IF(narg<1)THEN
   WRITE(*,*)'Please enter the base name of data file: '
   READ(*,*)argv
  ELSE
   CALL getarg(1,argv)
 ENDIF
 nlen=LNBLNK(argv)
 INQUIRE(FILE='.//subordinate//'+argv(1:nlen)//'.dat',EXIST=found)
 IF(.NOT.found)THEN
  WRITE(*,*)'Data file not found: ',argv(1:nlen)//'.dat'
  WRITE(*,*)'Please create or check spelling.'
  STOP
 ENDIF
RETURN
END SUBROUTINE getname


!
!     -----------------------------------------------------------------
!      R M E W D 2 L S
!     -----------------------------------------------------------------
!
!     Written by G. Gaigalas,                                      *
!     Vilnius,  Lithuania                              June 1995   *
!
      SUBROUTINE RMEWD2LS(J, L, COEF)
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE vast_kind_param, ONLY:  DOUBLE
!...Translated by Pacific-Sierra Research 77to90  4.3E  13:13:23  11/16/01
!...Switches:
      IMPLICIT NONE
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      INTEGER , INTENT(IN) :: J
      INTEGER , INTENT(IN) :: L
      REAL(DOUBLE) , INTENT(OUT) :: COEF
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      INTEGER , DIMENSION(136) :: ID111L, ID030L, ID131L
      INTEGER :: ID, IV
!-----------------------------------------------
      DATA ID111L/ 3*0, -2700, -420, 600, 0, 480, 3360, 240, 2*0, 4800, 1050, &
         3375, 3*0, 2430, -2625, -75, 2*0, -2100, 240, 0, -2400, 3*0, -1680, 30&
         , 0, 300, 6*0, -4200, 3000, 2400, 120, 8400, 4*0, 4200, 120, -960, &
         -1920, 840, -210, 5*0, -1200, 2*0, 2100, -525, 9*0, -5400, -1350, &
         -1215, 1458, 8*0, -1800, 4500, 0, -540, 9*0, -1980, -495, 0, -1485, 13&
         *0, 4752, 3960, 1440, 1188, 14*0, -3510, 0/
      DATA ID030L/ 14*0, -53900, 3*0, 39600, 0, 5940, 6*0, 28160, 7*0, -880, 2*&
         0, -77616, 5*0, 38808, 3*0, -693, 0, 34650, 3*0, -43659, 0, -12320, 5*&
         0, -28875, 2*0, 1617, 3*0, 40095, 0, -990, 3*0, -38115, 0, 16335, 6*0&
         , -39600, 5*0, -43560, 7*0, -27225, 2*0, -3465, 2*0, 23805, 5*0, 43560&
         , 3*0, 38115, 0, 70785, 13*0, 27027, 2*0, 28665, 0, 20384/
      DATA ID131L/ 12*0, 58800, -16800, 110250, 2*0, -42000, -7680, 15750, &
         -36450, 2*0, 33600, 29040, 0, 14400, 3*0, -1680, -5880, 0, 28800, 2*0&
         , -88200, 2520, 117600, 3360, 235200, 0, 8400, 15120, -58800, 0, 10080&
         , 11760, -35490, 14700, 15120, -41160, -945, 94080, 1470, 2*0, -29400&
         , -13125, 0, 9450, 3*0, 33075, 3*0, 75600, -1350, 94500, -10800, &
         -48600, 16875, 0, 103950, 10395, 16038, 2*0, -25200, 18000, 0, 36000, &
         2*0, -138600, 2*0, 95040, 3*0, -27720, 24255, 0, -24750, 2*0, 20160, &
         25515, 0, -135, 6*0, -115500, -29700, -26400, 0, 184800, -11550, 10395&
         , 7722, -102960, 28665, -168168, 8*0, 109200, -6825, 0, -85995, 2*0, &
         76440, 0/
      SELECT CASE (L)
      CASE (0)
         ID = ID111L(J)
         IV = 25
      CASE (1)
         ID = ID030L(J)
         IV = 770
      CASE (2)
         ID = ID131L(J)
         IV = 350
      END SELECT
      IF (ID >= 0) THEN
         COEF = DSQRT(DBLE(ID))
      ELSE
         COEF = -DSQRT((-DBLE(ID)))
      ENDIF
      COEF = COEF/DSQRT(DBLE(IV))
      RETURN
      END SUBROUTINE RMEWD2LS