text:
      SUBROUTINE DIRNA
      IMPLICIT NONE
C
C  Add to TOC variables of new equation of equinox procedure
C
C       PROGRAMMER:
C         93.09.01 Norbert Zacharias - Equation of the equinoxes GAST update
C         94.01.04 David Gordon - Cleaned up TOC entries.
C         98.01.20 David Gordon - New TOC entries for New definition of the
C                                 equation of the equinoxes. Removed Lcodes
C                                 for previous EQE (calc 8.x) contributions.
C
       INCLUDE 'ccon.i'
C            VARIABLES 'FROM':
C              1. KDIUC - THE DIURNAL SPIN UTILITY ROUTINE FLOW CONTROL FLAG.
C                         = 0 (default) Use new definition of the equation of
C                             the equinoxes.
C                         = 1 Use old definition of the equation of the
C                             equinoxes.
C
C  Lcode for difference in GAST (UT1)
C    [Even though we change the default to use the New definition, we keep this
C     Lcode the same as before so we don't have to change the name.]
       CALL ADDR (2,'EQE DIFF','GAST diff. (New-Old) from EQE   ',
     *           1,1,1)
C
C   Add contribution Lcode for the delay and rate contribution due to the
C   change in the definition of the equation of the equinox. The default is
C   now to use the New definition. Also remove previous EQE contribution
C   Lcodes that might be in the database.
        CALL DELR(2,'EQE CONT')
C    For the default (New definition), put in Lcode to convert back to the old
C    definition.
      If(KDIUC .eq. 0) Then
       CALL ADDR (2,'OLDEQCON','Old Eqn. of Equinoxes Contrib.  ',
     *           2,1,1)
       CALL DELR(2,'NEWEQCON')
      Endif
C
C    For the Old definition, put in Lcode to convert to the new definition.
      If(KDIUC .eq. 1)  Then
       CALL ADDR (2,'NEWEQCON','New Eqn. of Equinoxes Contrib.  ',
     *           2,1,1)
        CALL DELR(2,'OLDEQCON')
      Endif
C
      Return
      END
C
C***************************************************************************
      SUBROUTINE DIRNL ( DATDCT, DPSI, DUT1AT, EPS, DEPS, FA, UT1,
     1                   XJD, DIURNV, GAST, GMST, GASTD, RS )
      IMPLICIT None
C
C 1.    DIRNL
C
C 1.1   DIRNL PROGRAM SPECIFICATION
C
C 1.1.1 DIRNL is the utility routine which computes the diurnal spin portion of
C       the complete crust fixed to J2000.0 rotation matrix and its first two CT
C       time derivatives. DIRNL also computes the diurnal angular velocity of
C       the Earth, the Greenwich Mean Sidereal Time, and the Greenwich Apparent
C       Siderial Time (GAST) and its CT time derivative.
C
C 1.1.2 RESTRICTIONS - NONE
C
C 1.1.3 REFERENCES - 1) 'THE EXPLANATORY SUPPLEMENT TO THE AMERICAN
C                       EPHEMERIS AND NAUTICAL ALMANAC", P.72-76,
C                    2) AOKI, S. ET AL., "THE NEW DEFINITION OF UNIVERSAL
C                       TIME", ASTRON. ASTROPHYS., ????,1980.
C                    3) McCarthy, D., IERS Technical Note 13, Paris 1992
C
C 1.2   DIRNL PROGRAM INTERFACE
C
C 1.2.1 CALLING SEQUENCE -
C
C           INPUT VARIABLES:
C             1. DATDCT   -  THE PARTIAL DERIVATIVE OF ATOMIC TIME WITH
C                            RESPECT TO COORDINATE TIME. (SEC/SEC)
C             2. DPSI(2)  -  THE NUTATION IN LONGITUDE AND ITS CT
C                            TIME DERIVATIVE. (RAD, RAD/SEC)
C             3. DUT1AT   -  THE PARTIAL DERIVATIVE OF UT1 TIME WITH
C                            RESPECT TO ATOMIC TIME. (SEC/SEC)
C             4. EPS(2)   -  THE TRUE OBLIQUITY OF THE ECLIPTIC AND ITS
C                            CT TIME DERIVATIVE. (RAD, RAD/SEC)
C             5. DEPS(2)  -  THE NUTATION IN OBLIQUITY AND ITS CT TIME
C                            DERIVATIVE COMPUTED FROM WAHR OR TAKEN FROM
C                            THE DATA BASE. (RAD, RAD/SEC)
C             6. FA(5)    -  The fundamental nutation arguments (see NUTFA)
C             6. OMEGA    -  The mean longitude of the ascending node of the
C                            Moon. (arcsec!)
C             7. UT1      -  THE UT1 TIME OF THE DAY. (SEC)
C             8. XJD      -  THE JULIAN DATE AT ZERO HOURS UTC OF THE
C                            DATE IN QUESTION. (DAYS)
C
C           OUTPUT VARIABLES:
C             1. DIURNV    -  THE DIURNAL ANGULAR VELOCITY OF THE EARTH.
C                             (RAD/SEC)
C             2. GAST(2)   -  THE GREENWICH APPARENT SIDEREAL TIME AND
C                             ITS CT TIME DERIVATIVE. (RAD, RAD/SEC)
C             3. GMST      -  THE GREENWICH MEAN SIDEREAL TIME. (RAD)
C             4. GASTD     -  The difference in GAST (new-old version)
C             5. RS(3,3,3) -  THE DIURNAL SPIN PORTION OF THE COMPLETE
C                             CRUST FIXED TO 2000.0 ROTATION MATRIX
C                             AND THE FIRST TWO CT TIME DERIVATIVES OF
C                             THAT MATRIX. (UNITLESS, 1/SEC, 1/SEC**2)
C
      Real*8 DATDCT, DPSI(2), DUT1AT, EPS(2), DEPS(2), FA(5), UT1,
     .       XJD, DIURNV, GAST(2), GMST, GASTD, RS(3,3,3)
C
C 1.2.2 COMMON BLOCKS USED -
C
      Real*8           PI, TWOPI, HALFPI, CONVD, CONVDS, CONVHS, SECDAY
      COMMON / CMATH / PI, TWOPI, HALFPI, CONVD, CONVDS, CONVHS, SECDAY
      EXTERNAL CMATHB
C
C            VARIABLES 'FROM':
C              1.  CONVHS  -  THE CONVERSION FACTOR OF RADIANS PER TIME-SECOND.
C                             (RAD/TIME-SEC)
C
           INCLUDE 'ccon.i'
C            VARIABLES 'FROM':
C              1.  KDIUC  -  THE DIURNAL SPIN UTILITY ROUTINE FLOW CONTROL FLAG.
C              2.  KDIUD  -  THE DIURNAL SPIN UTILITY ROUTINE DEBUG OUTPUT FLAG.
C              3.  KNUTC  -  THE NUTATION MODULE FLOW CONTROL FLAG. USED TO
C                            DETERMINE WHETHER OR NOT DPSI SHOULD BE ZEROED.
C
      Real*8 GMSTC(4), SIDVEL(3), DPSIS(2), epsmd, omega, CENTJ, DJ2000,
     .       GMSTMN, DAYSJ, CENT
      Integer*4 I
C
      DATA  CENTJ  / 36525.D0 /,
     1      DJ2000 /  2451545.D0 /,
     2      GMSTC  / + 24110.54841D0, + 8640184.812866D0,
     3               + 0.093104D0, - 6.2D-6 /,
     4      SIDVEL / + 1.002737909350795D0, + 5.9006D-11, - 5.9D-15 /
C
C 1.2.4 DATA BASE ACCESS - NONE
C
C 1.2.5 EXTERNAL INPUT/OUTPUT -
C            INPUT VARIABLES - NONE
C            OUTPUT VARIABLES - POSSIBLE DEBUG OUTPUT
C
C 1.2.6 SUBROUTINE INTERFACE -
C             CALLER SUBROUTINES: DRIVR
C             CALLED SUBROUTINES: DDROT, DROTT, ROTAT, EQEG
C
C 1.2.7 CONSTANTS USED - CENTJ, CONVHS, DJ2000, GMSTC(4), SIDVEL(3)
C
C 1.2.7.1 CONSTANTS INITIALIZED IN THIS UTILITY ROUTINE -
C           1. CENTJ     -  THE NUMBER OF COORDINATE TIME DAYS PER JULIAN
C                           CENTURY (CENTJ = 36525.D0 DAYS/CENTURY)
C           2. DJ2000    -  THE JULIAN DATE OF THE EPOCH JANUARY 1.5, 2000.
C                           (DJ2000 = 2451545.D0 DAYS)
C           4. GMSTC(4)  -  THE CONSTANTS APPEARING IN TERMS 1-4 IN THE
C                           CALCULATION OF THE GREENWICH MEAN SIDEREAL TIME AT
C                           ZERO HOURS UT1 OF THE DATE IN QUESTION. (TIME-SEC)
C                           (SEE REFERENCES)
C                           ( GMSTC(1) = + 24110.54841D0,
C                             GMSTC(2) = + 8640184.812866D0,
C                             GMSTC(3) = + 0.093104D0,
C                             GMSTC(4) = - 6.2D-6 )
C           4. SIDVEL(3) -  THE CONSTANTS APPEARING IN TERMS 1-3 IN THE
C                           CALCULATION OF THE EARTH'S DIURNAL ANGULAR VELOCITY.
C                           (TIME SEC/SEC) (SEE REFERENCES)
C                           ( SIDVEL(1) = + 1.002737909350795D0,
C                             SIDVEL(2) = + 5.9006D-11,
C                             SIDVEL(3) = - 5.9D-15 )
C
C 1.2.8 PROGRAM VARIABLES -
C           1.  CENT    -  THE NUMBER OF JULIAN CENTURIES ELAPSED SINCE THE
C                          EPOCH JANUARY 1.5, 2000. (CENTURIES)
C           2.  DAYSJ   -  THE NUMBER OF JULIAN DAYS ELAPSED SINCE THE EPOCH
C                          JANUARY 1.5, 2000. (DAYS)
C           3.  GMSTMN  -  THE GREENWICH MEAN SIDEREAL TIME AT ZERO HOURS UT1
C                          OF THE DATE IN QUESTION. (RAD)
C           4. DPSIS(2) -  THE VALUE OF NUTATION USED IN THE CALCULATION OF
C                          SIDEREAL TIME.
C           5. EPSMD    -  Mean obliquity of the ecliptic at date. (rad)
C           6. OMEGA    -  Mean longitude of the ascending node of the Moon
C                          (arcsec)
C
C 1.2.9 PROGRAMMER
C         770207  DALE MARKHAM
C         770718  PETER DENATALE
C         771128  BRUCE SCHUPLER
C         790201  BRUCE SCHUPLER
C         810804  CHOPO MA
C         870604  SAVITA GOEL
C         890726  Jim Ryan:          Documentation simplified.
C         930901  Norbert Zacharias: Equation of equinoxes GAST update
C
C  DIRNL Program Structure
C
C   Compute mean obliquity of date, req. for EQE update in DIRNL
      epsmd = eps (1) - deps (1)
C   Mean longitude of ascending node of the Moon
      omega = fa (5)

C   Compute the number of Julian days elapsed since the epoch
C   JANUARY 1.5, 2000 to 0 hours UT1 of the date in question.
      DAYSJ = XJD  -  DJ2000
C
C   Compute the number of Julian centuries elapsed since Epoch January 1.5, 2000
      CENT = DAYSJ / CENTJ
C
C   Compute the diurnal angular velocity of the Earth.
      DIURNV = ( SIDVEL(1) + SIDVEL(2) * CENT + SIDVEL(3) * CENT**2 )
     1       * CONVHS
C
C   Compute the Greenwich Mean Sidereal Time at 0 hours UT1 of the date
C   in question. (radians)
      GMSTMN = ( GMSTC(1) + GMSTC(2) * CENT + GMSTC(3) * CENT**2 +
     .           GMSTC(4) * CENT**3 ) * CONVHS
C
C   Compute the Greenwich Mean Sidereal Time.
      GMST = GMSTMN  +  DIURNV * UT1
      GMST = DMOD(GMST,TWOPI)
C
C   Check to see if the nutation in longitude is to be turned off.
      DO I=1,2
        IF (KNUTC .EQ. 1) DPSIS(I) = 0.0D0
        IF (KNUTC .NE. 1) DPSIS(I) = DPSI(I)
      ENDDO
C
C  Compute the Greenwich Apparent Sidereal Time by calling the geometry part of
C   the equation of equinox procedure.
      CALL EQEG (gmst, dpsis(1), eps(1), epsmd, omega,
     .           gast, gastd)
C
C  Compute the time derivative of Greenwich Apparent Siderial Time
      GAST(2) = DIURNV * DUT1AT * DATDCT
     1        + DCOS ( EPS(1) ) * DPSIS(2)
     2        - DSIN ( EPS(1) ) * DPSIS(1) * EPS(2)
C
C   Construct the diurnal spin matrix and its first two CT time derivatives.
C
C   Construct the diurnal spin matrix.
      CALL ROTAT ( -GAST(1), 3, RS(1,1,1) )
C
C   Construct the first CT time derivative of the diurnal spin matrix
      CALL DROTT ( -GAST(1), -GAST(2), 3, RS(1,1,2) )
C
C   Construct the second CT time derivative of the diurnal spin matrix.
      CALL DDROT ( -GAST(1), GAST(2)**2, 3, RS(1,1,3) )
C
C****************************************************************************
C   Check KDIUD for debug output.
      IF ( KDIUD .NE. 0 )  THEN
        WRITE ( 6, 9)
    9   FORMAT (1X, "Debug output for utillity DIRNL." )
        WRITE(6,8)' CENT    ',CENT
    8   FORMAT(A,4D25.16/(7X,5D25.16))
        WRITE(6,8)' CENTJ   ',CENTJ
        WRITE(6,8)' DJ2000  ',DJ2000
        WRITE(6,8)' DPSIS   ',DPSIS
        WRITE(6,8)' CONVHS  ',CONVHS
        WRITE(6,8)' DAYSJ   ',DAYSJ
        WRITE(6,8)' GMSTC   ',GMSTC
        WRITE(6,8)' GMSTMN  ',GMSTMN
        WRITE(6,8)' SIDVEL  ',SIDVEL
C
        WRITE ( 6, 9200 )  DATDCT, DPSI, DUT1AT, EPS, UT1, XJD,
     1                     DIURNV, GAST, GMST, RS
 9200   FORMAT (1X, "DATDCT = ", D30.16, /, 1X,
     1            "DPSI   = ", 2 ( D30.16, 10X ), /, 1X,
     2            "DUT1AT = ", D30.16, /, 1X,
     3            "EPS    = ", 2 ( D30.16, 10X ), /, 1X,
     4            "UT1    = ", D30.16, /, 1X,
     5            "XJD    = ", D30.16, /, 1X,
     6            "DIURNV = ", D30.16, /, 1X,
     7            "GAST   = ", 2 ( D30.16, 10X ), /, 1X,
     8            "GMST   = ", D30.16, /, 1X,
     9            "RS     = ", 9 ( 3 ( D30.16, 10X ), /, 1X ) )
C
      ENDIF   ! debug output flag
C
      Return
      END
C
C******************************************************************************
      SUBROUTINE EQEG (gmst,dpsi,eps,epsmd,omega, gast,gastd)
      IMPLICIT NONE
C
C  This subroutine computes the geometry section of the equation of equinox
C  procedure.
C
C  References:  D.McCarthy, IERS Technical Note 13, Paris 1992
C
C  Calling sequence:
C    input variables:
C      gmst  = Greenwich mean siderial time           (radian)
C      dpsi  = nutation in longitude                  (radian)
C      eps   = true obliquity of the ecliptic         (radian)
C      epsmd = obliquity of the mean ecliptic of date (radian)
C      omega = mean longitude of the asc. node of the moon (arcsec)
C
C    output variables:
C      GAST(1)=Greenwich apparent siderial time (New or Old depending on KDIUC)
C               (radians)
C      GASTD = Difference in GAST of two versions (New-Old) (radians)
C
      REAL*8   gmst, dpsi, eps, epsmd, omega, gast(2), gastd
      REAL*8   old_gast, aswas_gast, new_gast,
     .         eqe_const(2), omega_rad
      REAL*8   pi, twopi, halfpi, convd, convds, convhs, secday
      COMMON /CMATH/ pi, twopi, halfpi, convd, convds, convhs, secday
      EXTERNAL CMATHB
C
C     Variables from:  convds = convert arcsec to radian
C
       INCLUDE 'ccon.i'
C        Variables from:
C           KDIUC - Flow control flag for EQE routines
C                   = 0, GAST defined by new equation of equinoxes definition
C                   = 1, GAST defined by old equation of equinoxes definition
C           KDIUD - Debug flag for EQE routines
C
C  Data base access: None
C
C  Subroutine interface:
C     caller subr.: DIRNL
C     called subr.: None
C
C  Program variables:
C     eqe_const   = Constants in new version of GAST formula
C                   depending on sine of omega  (arcsec)
C     aswas_gast  = GAST as calculated incorrectly before (CALC 7.6 and
C                   earlier). (radian)
C     old_gast    = GAST as calculated before (CALC 7.6 and earlier) but
C                   corrected to use mean obliquity of date, to be used until
C                   1997.02.26. (radians)
C     omega_rad   = Omega (longit. of asc. node of moon) (radians)
C     new_gast    = GAST for use after 1997.02.26 as recommended by the IERS.
C
C  Programmer:
C         930901 Norbert Zacharias: first implementation, equat.of eq. GAST
C                                   update.
C         94.01.04 David Gordon - Clean up documentation, change sign of GAST
C                                 correction.
C         98.01.20 David Gordon - Updated to use new definition EQE
C
      DATA eqe_const / 0.00264D0, 0.000063D0 /  ! in arc sec
C
C  compute old GAST as before, use true obliquity eps as before
      aswas_gast = gmst + dpsi * DCOS (eps)
C
C  compute old GAST (formula in effect until 19970226)
C  but with correct epsilon (mean obliquity of date)
      old_gast = gmst + dpsi * DCOS (epsmd)
C
C  compute new GAST as recomm. in IERS techn.note
      omega_rad = omega * convds
      new_gast  = gmst + dpsi * DCOS (epsmd)
     .            + eqe_const(1) * convds * DSIN (omega_rad)
     .            + eqe_const(2) * convds * DSIN (2.0d0*omega_rad)
C
C  Calc 9, use new definition
      If (KDIUC.eq.0) gast(1)  = new_gast
C  For diehards, allow using old definition
      If (KDIUC.eq.1) gast(1)  = old_gast
C  Define difference, Gastd, as (New - Old)
      gastd = new_gast - old_gast
C
C  check debug output flag
      IF (KDIUD.NE.0) THEN
	WRITE (6,'(1x,a)') 'debug output for EQEG'
	WRITE (6,'(1x,a,d25.16)') 'gmst       = ', gmst
	WRITE (6,'(1x,a,d25.16)') 'dpsi       = ', dpsi
	WRITE (6,'(1x,a,d25.16)') 'eps        = ', eps
	WRITE (6,'(1x,a,d25.16)') 'epsmd      = ', epsmd
	WRITE (6,'(1x,a,d25.16)') 'omega      = ', omega
	WRITE (6,'(1x,a,d25.16)') 'aswas_gast = ', aswas_gast
	WRITE (6,'(1x,a,d25.16)') 'old_gast   = ', old_gast
	WRITE (6,'(1x,a,d25.16)') 'new_gast   = ', new_gast
	WRITE (6,'(1x,a,d25.16)') 'gast(1)    = ', gast(1)
	WRITE (6,'(1x,a,d25.16)') 'aswas-old  = ', aswas_gast-old_gast
	WRITE (6,'(1x,a,d25.16)') 'gastd      = ', gastd
      ENDIF
C
      Return
      END
C ****************************************************************************
      SUBROUTINE DIRNC (gastd, dut1p)
      IMPLICIT NONE
C
C  This subr. computes the contributions to the delay and delay rate due to the
C  different versions of the equation of equinox procedure. The New version
C  (to be used after 26 Feb 1997) becomes the default in Calc 9.0. To get back
C  to the old definition, simply add the contribution Lcode OLDEQCON during
C  SOLVE analysis.
C
C  Users can still use the old EQE definition by setting KDIUC=1. In that case,
C  a diferent contribution goes into the database, 'NEWEQCON'. When NEWEQCON
C  is added to the delays and rates in this case, they will be converted to
C  the new definition. Don't count on this being supported in SOLVE though.
C
C  References: D.McCarthy, IERS Technical Note 13, Paris 1992
C              E-Mail Memo, 'UT1 Errors', D. Gordon, 10 April, 1997
C
C  Calling sequence:
C    Input Variables:
C       GASTD      = Differ. in GAST of the two versions (new-old) (radian)
C       DUT1P(2,2) = Partial derivatives of delay and rate wrt (A1-UT1)
C                    (First index - delay and rate, second index - first and
C                     second derivatives) (s/s, s/s**2)
C         [Note: The partial we want here is actually the delay and rate
C                derivatives wrt GAST. We approximate this by dividing
C                dut1p by -1.0027379 sidereal days/solar day.]
C
      REAL*8   gastd, dut1p(2,2), deqec(2)
      REAL*8   pi,twopi, halfpi, convd, convds, convhs, secday
      COMMON /CMATH/ pi,twopi,halfpi, convd, convds,
     .               convhs, secday
      EXtERNAL CMATHB
C      Variables From:
C           Convhs = conversion factor (radian/sec.of time)
C
       INCLUDE 'ccon.i'
C      Variables From:
C           KDIUC - Flow control flag for EQE routines
C                   = 0, GAST defined by new equation of equinoxes definition
C                   = 1, GAST defined by old equation of equinoxes definition
C           KDIUD - debug flag for EQE routines
C
C  Data base access:
C     PUT4 GASTD    = differ. in GAST (new - old) (radian)
C     PUT4 DEQEC(2) = contrib. of GASTD to delay and delay rate
C                     (sec, sec/sec)
C  Subroutine interface:
C     Caller Subr.: DRIVR
C     Called Subr.: None
C
C  Programmer:
C    930901 Norbert Zacharias: first implementation
C    94.01.04 David Gordon - Fixed up documentation, sign of contributions
C                            flipped in EQEG.
C    94.06.02 David Gordon - Changed PUTR to PUT4.
C    95.02.09 David Gordon - Corrected PUT4's of 'EQE DIFF' and 'EQE CONT'.
C                            Only 8 charcters allocated, changed to 14. Bug
C                            found by Bill Petrachenko using a Sun workstation
C                            at the Dominion Radio Astrophysical Obs. (DRAO).
C    98.01.20 David Gordon - Corrected old (Calc 8.x) delay and rate
C                            contribution (Deqec) by dividing by -1.0027379
C                            siderial days/solar day. Changed database Lcode to
C                            'NEWEQCON' (from 'EQE CONT') and made it the
C                            non-default. Added 'OLDEQCON' Lcode contribution
C                            and made it the default.
C
C  Calculate contributions to delay and rate
C
C   Default case, contribution to convert from New to Old EQE definition
      If (KDIUC .eq. 0) then
       deqec (1) = -gastd * dut1p(1,1) / (-1.0027379 * convhs)
       deqec (2) = -gastd * dut1p(2,1) / (-1.0027379 * convhs)
       CALL PUT4 ('OLDEQCON      ', deqec, 2,1,1)
      Endif
C
C   Non-default case, contribution to convert from Old to New EQE definition
      If (KDIUC .eq. 1) then
C    Correct Calc 8.2 error by dividing by -1.0027379
       deqec (1) = gastd * dut1p(1,1) / (-1.0027379 * convhs)
       deqec (2) = gastd * dut1p(2,1) / (-1.0027379 * convhs)
       CALL PUT4 ('NEWEQCON      ', deqec, 2,1,1)
      Endif
C
C  Output EQE difference to data base
      CALL PUT4 ('EQE DIFF      ', gastd, 1,1,1)
C
C  check flag for debug output:
      IF (KDIUD.NE.0) THEN
	WRITE (6,'(1x,a)') 'debug output for EQEC'
	WRITE (6,'(1x,a,2d25.16)')
     .    'dut1p         = ', dut1p(1,1), dut1p(2,1)
	WRITE (6,'(1x,a,2d25.16)')
     .    'gastd, convhs = ', gastd, convhs
	WRITE (6,'(1x,a,2d25.16)') 'deqec         = ', deqec
      ENDIF
C
      Return
      END

max_stars_repo_name: liulei/VOLKS

max_issues_repo_name: liulei/VOLKS

max_forks_repo_name: liulei/VOLKS


text:
 Mission Statement
The Interfaith Campus Council (ICC) is an association of campus religious organizations, represented by their professional staff and/or student leaders who are committed to working together collegially and cooperatively.  Although we come from different religious and denominational traditions, all of our organizations consider interfaith engagement to be an essential part of our own programs as well as a vital shared service that we offer to UC Davis and especially its students. We seek to promote mutual respect, understanding, and appreciation among people of different faiths, and we work together to create shared programs and activities that give witness to our common values and ideals.

If your campus religious organization would like to become a member of the Interfaith Campus Council, or if you represent a campus unit or student organization that would like to work with the ICC, please email our Convener mailto:campusministry@uudavis.org Aidan McCormack for our next meeting date.

 What Our Member Groups Offer
  Worship and Holiday Services
  Retreat Opportunities
  Sacred Text Studies
  Meals and Discussions
  Fellowship and Social Groups
  Social Concerns Activities
  Peer Ministry
  Conversation and Counseling about faith, ethics, personal problems, and career questions

 Annual Interfaith Activities
  Interfaith Progressive Dinner
  Martin Luther King Peace and Justice Week
  Sex and the Texts: Religious Perspectives on Sexual Orientations Panel (Cosponsored by the UCD LGBT Resource Center and Counseling and Psychological Services)
  Humanitarian of the Year Award

 Our BIG Project!
The Multifaith Living Community at UC Davis
The Multifaith Living Community is an innovative project, sponsored by CA House and supported by the Interfaith Campus Council at UC Davis, to gather 34 students per year of different faiths to live in a community together.

Through their own faith journeys, shared meals and programs, and community service and justice work, these students will provide a model of peace, understanding and cooperation among people of different faiths and cultures.

To house this new community, CA House is building six new townhouses on their property at 433 Russell Boulevard, directly across from the UC Davis campus, and organizing student and community supporters around this project.

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
hello people!
this is my first time in a wiki website. so excuse all of the random editing mistakes. i promise to get better in time.
i happen to love making food, eating food, and drinking, so you will see many of my comments on those pages...
i am student at ucdavis, majoring in biology; and i will be a student forever.

Welcome to the Wiki. You might find some of the How to Edit editing tips helpful.  Users/ElleWeber

20080612 06:31:32 nbsp Ahhh, to be a student forever...and a student of food and wine, no less! Users/ScottLay


text:
c
c --------------------------------------------------------------
c
      subroutine bdcrse( xleft,xright,ybot,ytop,hx,hy,time,intime,
     1                   level,nvar,lw2,valdub,midub,mjdub,valbig,
     2                   mi2tot,mj2tot,irr,mitot,mjtot)
c
c  bdcrse: first set double the stencil size worth of boundary values,
c          then coarsen them for the giant step integration.
c
      implicit double precision (a-h,o-z)
      include "cirr.i"
      dimension valdub(midub,mjdub,nvar), valbig(mi2tot,mj2tot,nvar),
     1          irr(mitot,mjtot)

c     left boundary

      xl = xleft - lw2*hx
      yb = ybot - lw2*hy
      xr = xleft
      yt = ytop

      call filpc(xl,xr,yb,yt,level,nvar,valdub,time,intime,midub,
     1           mjdub,1,lw2,1,mjdub-lw2,irr,mitot,mjtot)

c     right boundary

      xl = xright
      yb = ybot
      xr = xright + lw2*hx
      yt = ytop + lw2*hy

      call filpc(xl,xr,yb,yt,level,nvar,valdub,time,intime,midub,
     1           mjdub,midub-lw2+1,midub,lw2+1,mjdub,irr,mitot,mjtot)

c     bottom boundary

      xl = xleft
      yb = ybot - lw2*hy
      xr = xright + lw2*hx
      yt = ybot

      call filpc(xl,xr,yb,yt,level,nvar,valdub,time,intime,midub,
     1           mjdub,lw2+1,midub,1,lw2,irr,mitot,mjtot)

c     top boundary

      xl = xleft - lw2*hx
      yb = ytop
      xr = xright
      yt = ytop + lw2*hy

      call filpc(xl,xr,yb,yt,level,nvar,valdub,time,intime,midub,
     1           mjdub,1,midub-lw2,mjdub-lw2+1,mjdub,irr,mitot,mjtot)

c
c now coarsen by 2 in every direction - conservatively, so that
c irregular and exterior cells count appropriately.
c irr array only exists for lwidth, not 2*lwidth sized grid.
c boundary cells not done correctly here.
c
      lwidth = lw2 / 2
      ar(-1) = 0.d0
      do 100 ivar = 1,nvar
      do 100 j = 1, mj2tot
         jfine = 2*(j - 1) + 1
         do 100 i = 1, mi2tot
            ifine = 2*(i - 1) + 1
	    tarea = 0.d0
	    if ((ifine .le. lwidth) .or. (jfine .le. lwidth) .or.
     .          (ifine+1 .gt. midub-lwidth) .or.
     .          (jfine+1 .gt. mjdub-lwidth)) then
	       a11 = 1.d0
	       a21 = 1.d0
	       a12 = 1.d0
	       a22 = 1.d0
	    else
	       ieff = ifine - lwidth
	       jeff = jfine - lwidth
	       a11   = ar(irr(ieff,jeff))
	       a21   = ar(irr(ieff+1,jeff))
	       a12   = ar(irr(ieff,jeff+1))
	       a22   = ar(irr(ieff+1,jeff+1))
	    endif
c
	    tarea = a11+a12+a21+a22
c
            if (tarea .eq. 0.d0) then
	      valbig(i,j,ivar) = valdub(ifine,jfine,ivar)
	    else
              valbig(i,j,ivar)=(valdub(ifine,jfine,ivar)*a11 +
     &                    valdub(ifine+1,jfine,ivar)*a21+
     &                    valdub(ifine,jfine+1,ivar)*a12 +
     &                    valdub(ifine+1,jfine+1,ivar)*a22)/tarea
	    endif
100   continue
c
c  external boundary conditions
c
       lwidth = lw2 / 2
       hx2    = 2.d0*hx
       hy2    = 2.d0*hy
       xl = xleft  - lwidth*hx2
       yb = ybot   - lwidth*hy2
       xr = xright + lwidth*hx2
       yt = ytop   + lwidth*hy2
       call physbd(xl,xr,yb,yt,level,mi2tot,mj2tot,nvar,valbig,
     1             time,hx2,hy2)
c
      return
      end


text:
      subroutine  dualfa
!***********************************************************************
!  Copyright, 2004,  The  Regents  of the  University of California.
!  This program was prepared by the Regents of the University of
!  California at Los Alamos National Laboratory (the University) under
!  contract No. W-7405-ENG-36 with the U.S. Department of Energy (DOE).
!  All rights in the program are reserved by the DOE and the University.
!  Permission is granted to the public to copy and use this software
!  without charge, provided that this Notice and any statement of
!  authorship are reproduced on all copies. Neither the U.S. Government
!  nor the University makes any warranty, express or implied, or
!  assumes any liability or responsibility for the use of this software.
C**********************************************************************
CD1
CD1 PURPOSE
CD1
CD1 To compute the Jacobian and residual terms of the mass balance
CD1 equations at each node for a dual porosity isothermal air-water
CD1 solution.
CD1
C**********************************************************************
CD2
CD2 REVISION HISTORY
CD2
CD2 Revision                    ECD
CD2 Date         Programmer     Number  Comments
CD2
CD2 05-17-93     G. Zyvoloski   00022   Initial implementation.
CD2
CD2 $Log:   /pvcs.config/fehm90/src/dualfa.f_a  $
!D2
!D2    Rev 2.5   06 Jan 2004 10:42:54   pvcs
!D2 FEHM Version 2.21, STN 10086-2.21-00, Qualified October 2003
!D2
!D2    Rev 2.4   29 Jan 2003 09:03:08   pvcs
!D2 FEHM Version 2.20, STN 10086-2.20-00
!D2
!D2    Rev 2.3   14 Nov 2001 13:08:40   pvcs
!D2 FEHM Version 2.12, STN 10086-2.12-00
!D2
!D2    Rev 2.2   06 Jun 2001 13:23:28   pvcs
!D2 FEHM Version 2.11, STN 10086-2.11-00
!D2
!D2    Rev 2.1   30 Nov 2000 11:59:58   pvcs
!D2 FEHM Version 2.10, STN 10086-2.10-00
!D2
!D2    Rev 2.0   Fri May 07 14:41:00 1999   pvcs
!D2 FEHM Version 2.0, SC-194 (Fortran 90)
CD2
CD2    Rev 1.4   Mon Jan 29 15:34:46 1996   hend
CD2 Updated Requirements Traceability
CD2
CD2    Rev 1.3   08/18/95 10:18:54   llt
CD2 neigh was already defined, removed for cray
CD2
CD2    Rev 1.2   12/09/94 16:02:16   llt
CD2 Changes to fix  radi (made by gaz).
CD2
CD2    Rev 1.1   03/18/94 15:50:06   gaz
CD2 Added solve_new and cleaned up memory management.
CD2
CD2    Rev 1.0   01/20/94 10:23:24   pvcs
CD2 original version in process of being certified
CD2
C**********************************************************************
CD3
CD3 INTERFACES
CD3
CD3 Formal Calling Parameters
CD3
CD3 None
CD3
CD3 Interface Tables
CD3
CD3 None
CD3
CD3 Files
CD3
CD3 None
CD3
C**********************************************************************
CD4
CD4 GLOBAL OBJECTS
CD4
CD4 Global Constants
CD4
CD4 Identifier  Type     Description
CD4
CD4 Global Types
CD4
CD4 None
CD4
CD4 Global Variables
CD4
CD4                      COMMON
CD4 Identifier   Type    Block   Description
CD4
CD4 lenreal, nn, neq, nelmdg, sx1, apuv1, tmch, phi, pcp,
CD4 pnx, pny, pnz, dil, div, s,
CD4 dnwgt, upwgt, aw, ay, a, bp,
CD4 wb11, tb11, a21mpf, a32mpf, rb2mf, r3mf,
CD4  t1, t2, t3, t4, t8, t9, zero_t
CD4
CD4
CD4 ipb, ipt1, ipt2, ipt3, ipt4, ipt8, ipt9, ipt10,
CD4 ipit8, ipit9, ipit10, ipit11, ipit12, istrw, dpcef, enlf,
CD4 delf, delef, envf, devf, devef, dilp, dile,
CD4 divp, dive, thx, thy,
CD4 thz, t, s, perml, permv,
CD4 sk, deni, denj, wb12, wb21, wb22, tb12,
CD4 tb21, tb22, a21mef, a21epf,
CD4 a21eef, a32mef, a32epf, a32eef, rb2ef, r3ef, nbd
CD4
CD4 Global Subprograms
CD4
CD4 Name       Type     Description
CD4
CD4 mmgetblk   N/A      Allocates memory for an array
CD4 mmrelblk   N/A      Frees array space
CD4 thermc     N/A      Computes solute storage and reaction terms
CD4
C**********************************************************************
CD5
CD5 LOCAL IDENTIFIERS
CD5
CD5 Local Constants
CD5
CD5 Identifier   Type        Description
CD5
CD5 Local Types
CD5
CD5 None
CD5
CD5 Local Variables
CD5
CD5 Identifier   Type        Description
CD5
CD5 sx1d         real*8      Volume associated with the current node
CD5 pvii         real*8      Pressure of vapor
CD5 phii         real*8      Pressure of liquid
CD5 swi          real*8      Saturation
CD5 dili         real*8      Liquid transmissibility
CD5 dilkb        real*8      Liquid transmissibility, connected node
CD5 divi         real*8      Vapor transmissibility
CD5 divkb        real*8      Vapor transmissibility, connected node
CD5 id           int         Integer index parameter
CD5 idg          int         Integer index parameter
CD5 jm           int        Integer index parameter
CD5 neqp1        int         Number of equations plus 1
CD5 neigh        int         Integer index parameter
CD5 neighc       int         Integer index parameter
CD5 axkb         real*8      Permeability in x-direction, connected node
CD5 sx2c         real*8      Real parameter used in calculation
CD5 pvikb        real*8      Vapor pressure at connected node
CD5 phikb        real*8      Liquid pressure at connected node
CD5 radi         real*8      Parameter used in calculation
CD5 radkb        real*8      Parameter used in calculation
CD5 fid          real*8      Parameter used in calculation
CD5 fid1         real*8      Parameter used in calculation
CD5 axy          real*8      Parameter used in calculation
CD5 axyd         real*8      Parameter used in calculation
CD5 axyf         real*8      Parameter used in calculation
CD5 vxy          real*8      Parameter used in calculation
CD5 vxyf         real*8      Parameter used in calculation
CD5 pxyh         real*8      Derivative term used in calculation
CD5 vxyd         real*8      Derivative term used in calculation
CD5 icode        int         Returned flag from call to mmgetblk
CD5 neq2         int         Twice the number of equations
CD5 idum         int         Index for arrays
CD5 tot          real*8      Total volume at a node
CD5 idl          int         Unknown number of first matrix node
CD5 idll         int         Unknown number of second matrix node
CD5 frac0        real*8      Volume fraction of primary node
CD5 frac1        real*8      Volume fraction of first matrix node
CD5 frac2        real*8      Volume fraction of second matrix node
CD5 alen         real*8      Volume to area ratio
CD5 area         real*8      Area term used in calculation
CD5 al0          real*8      Length scale for primary node
CD5 al1          real*8      Length scale for first matrix node
CD5 al2          real*8      Length scale for second matrix node
CD5 dist01       real*8      Average length scale used in calculation
CD5 dist12       real*8      Average length scale of matrix nodes
CD5 sx1d0        real*8      Volume of primary node
CD5 sx1d1        real*8      Volume of first matrix node
CD5 sx1d2        real*8      Volume of second matrix node
CD5 coef1        real*8      Length scale used in calculation
CD5 coef2        real*8      Length scale used in calculation
CD5 i            int         Index for first matrix node number
CD5 kb           int         Index for second matrix node number
CD5 iz           int         Index for primary node number
CD5 kz           int         Index for primary node number
CD5 sx2t         real*8      Parameter used in calculation
CD5 pxy          real*8      Parameter used in calculation
CD5 pvxy         real*8      Parameter used in calculation
CD5 pxyi         real*8      Parameter used in calculation
CD5 a1           real*8      Argument to statement function definition
CD5 a2           real*8      Argument to statement function definition
CD5 a3           real*8      Argument to statement function definition
CD5 a4           real*8      Argument to statement function definition
CD5 b1           real*8      Argument to statement function definition
CD5 b2           real*8      Argument to statement function definition
CD5 iwa          int         Integer index used in calculation
CD5 r1m          real*8      Term used in correction of residual of
CD5                             mass balance in matrix nodes
CD5 r2m          real*8      Term used in correction of residual of
CD5                             mass balance in matrix nodes
CD5 r3m          real*8      Term used in correction of residual of
CD5                             mass balance in matrix nodes
CD5 r1e          real*8      Term used in correction of residual of
CD5                             energy balance in matrix nodes
CD5 r2e          real*8      Term used in correction of residual of
CD5                             energy balance in matrix nodes
CD5 r3e          real*8      Term used in correction of residual of
CD5                             energy balance in matrix nodes
CD5 a11mp        real*8      Derivative term used in mass balance
CD5                             (mass with respect to pressure)
CD5 a12mp        real*8      Derivative term used in mass balance
CD5                             (mass with respect to pressure)
CD5 a21mp        real*8      Derivative term used in mass balance
CD5                             (mass with respect to pressure)
CD5 a22mp        real*8      Derivative term used in mass balance
CD5                             (mass with respect to pressure)
CD5 a23mp        real*8      Derivative term used in mass balance
CD5                             (mass with respect to pressure)
CD5 a32mp        real*8      Derivative term used in mass balance
CD5                             (mass with respect to pressure)
CD5 a33mp        real*8      Derivative term used in mass balance
CD5                             (mass with respect to pressure)
CD5 a11me        real*8      Derivative term used in mass balance
CD5                             (mass with respect to energy)
CD5 a12me        real*8      Derivative term used in mass balance
CD5                             (mass with respect to energy)
CD5 a21me        real*8      Derivative term used in mass balance
CD5                             (mass with respect to energy)
CD5 a22me        real*8      Derivative term used in mass balance
CD5                             (mass with respect to energy)
CD5 a23me        real*8      Derivative term used in mass balance
CD5                             (mass with respect to energy)
CD5 a32me        real*8      Derivative term used in mass balance
CD5                             (mass with respect to energy)
CD5 a33me        real*8      Derivative term used in mass balance
CD5                             (mass with respect to energy)
CD5 a11ep        real*8      Derivative term used in energy balance
CD5                             (energy with respect to pressure)
CD5 a12ep        real*8      Derivative term used in energy balance
CD5                             (energy with respect to pressure)
CD5 a21ep        real*8      Derivative term used in energy balance
CD5                             (energy with respect to pressure)
CD5 a22ep        real*8      Derivative term used in energy balance
CD5                             (energy with respect to pressure)
CD5 a23ep        real*8      Derivative term used in energy balance
CD5                             (energy with respect to pressure)
CD5 a32ep        real*8      Derivative term used in energy balance
CD5                             (energy with respect to pressure)
CD5 a33ep        real*8      Derivative term used in energy balance
CD5                             (energy with respect to pressure)
CD5 a11ee        real*8      Derivative term used in energy balance
CD5                             (energy with respect to energy)
CD5 a12ee        real*8      Derivative term used in energy balance
CD5                             (energy with respect to energy)
CD5 a21ee        real*8      Derivative term used in energy balance
CD5                             (energy with respect to energy)
CD5 a22ee        real*8      Derivative term used in energy balance
CD5                             (energy with respect to energy)
CD5 a23ee        real*8      Derivative term used in energy balance
CD5                             (energy with respect to energy)
CD5 a32ee        real*8      Derivative term used in energy balance
CD5                             (energy with respect to energy)
CD5 a33ee        real*8      Derivative term used in energy balance
CD5                             (energy with respect to energy)
CD5 dpvti        real*8      Cap. pressure derivative term
CD5 enli         real*8      Liquid enthalpy
CD5 enlkb        real*8      Liquid enthalpy, connected node
CD5 envkb        real*8      Vapor enthalpy, connected node
CD5 deli         real*8      Enthalpy derivative with pressure
CD5 delkb        real*8      Enthalpy derivative with pressure,
CD5                             connected node
CD5 devkb        real*8      Enthalpy derivative with pressure (vapor),
CD5                             connected node
CD5 delei        real*8      Enthalpy derivative with energy
CD5 delekb       real*8      Enthalpy derivative with energy,
CD5                             connected node
CD5 devekb       real*8      Enthalpy derivative with energy (vapor),
CD5                             connected node
CD5 envi         real*8      Vapor enthalpy
CD5 devi         real*8      Enthalpy derivative with pressure (vapor)
CD5 devei        real*8      Enthalpy derivative with energy (vapor)
CD5 dilpi        real*8      Trans. derivative with pressure
CD5 dilpkb       real*8      Trans. derivative with pressure,
CD5                             connected node
CD5 dilei       real*8      Trans. derivative with enthalpy
CD5 dilekb      real*8      Trans. derivative with enthalpy, connected
CD5                             node
CD5 divpi        real*8      Trans. derivative with pressure (vapor)
CD5 divpkb       real*8      Trans. derivative with pressure (vapor),
CD5                             connected node
CD5 divei        real*8      Trans. derivative with enthalpy (vapor)
CD5 divekb       real*8      Trans. derivative with enthalpy (vapor),
CD5                             connected node
CD5 thxi         real*8      Thermal conductivity, x-direction
CD5 thyi         real*8      Thermal conductivity, y-direction
CD5 thzi         real*8      Thermal conductivity, z-direction
CD5 ti           real*8      Temperature
CD5 athkb        real*8      Thermal conductivity, matrix nodes
CD5 thxkb        real*8      Thermal conductivity, matrix nodes
CD5 isl          int         Flag denoting if liquid or vapor exists
CD5 dlpi         real*8      Derivative term used in calculation
CD5 dvpi         real*8      Derivative term used in calculation
CD5 dvpkb        real*8      Derivative term used in calculation
CD5 dvaekb       real*8      Derivative term used in calculation
CD5 dlpkb        real*8      Derivative term used in calculation
CD5 dlei         real*8      Derivative term used in calculation
CD5 dvei         real*8      Derivative term used in calculation
CD5 dlekb        real*8      Derivative term used in calculation
CD5 dvekb        real*8      Derivative term used in calculation
CD5 dlapi        real*8      Derivative term used in calculation
CD5 dlapkb       real*8      Derivative term used in calculation
CD5 dlaei        real*8      Derivative term used in calculation
CD5 dvaei        real*8      Derivative term used in calculation
CD5 dlaekb       real*8      Derivative term used in calculation
CD5 dvapi        real*8      Derivative term used in calculation
CD5 dvapkb       real*8      Derivative term used in calculation
CD5 t11          real*8      Intermediate term used in calculation
CD5 t12          real*8      Intermediate term used in calculation
CD5 t21          real*8      Intermediate term used in calculation
CD5 t22          real*8      Intermediate term used in calculation
CD5 u11          real*8      Intermediate term used in calculation
CD5 u12          real*8      Intermediate term used in calculation
CD5 u21          real*8      Intermediate term used in calculation
CD5 u22          real*8      Intermediate term used in calculation
CD5 v11          real*8      Intermediate term used in calculation
CD5 v12          real*8      Intermediate term used in calculation
CD5 v21          real*8      Intermediate term used in calculation
CD5 v22          real*8      Intermediate term used in calculation
CD5 ab22mp       real*8      Intermediate term used in calculation
CD5 ab22me       real*8      Intermediate term used in calculation
CD5 ab22ep       real*8      Intermediate term used in calculation
CD5 ab22ee       real*8      Intermediate term used in calculation
CD5 w11          real*8      Intermediate term used in calculation
CD5 w12          real*8      Intermediate term used in calculation
CD5 w21          real*8      Intermediate term used in calculation
CD5 w22          real*8      Intermediate term used in calculation
CD5 x11          real*8      Intermediate term used in calculation
CD5 x12          real*8      Intermediate term used in calculation
CD5 x21          real*8      Intermediate term used in calculation
CD5 x22          real*8      Intermediate term used in calculation
CD5 y11          real*8      Intermediate term used in calculation
CD5 y12          real*8      Intermediate term used in calculation
CD5 y21          real*8      Intermediate term used in calculation
CD5 y22          real*8      Intermediate term used in calculation
CD5 rb2m         real*8      Intermediate term used in calculation
CD5 rb2e         real*8      Intermediate term used in calculation
CD5
CD5
CD5
CD5 Local Subprograms
CD5
CD5 Identifier    Type      Description
CD5
CD5 ali           real*8    Function used in matrix calculation
CD5 alm           real*8    Function used in matrix calculation
CD5
C**********************************************************************
CD6
CD6 FUNCTIONAL DESCRIPTION
CD6
CD6
C**********************************************************************
CD7
CD7 ASSUMPTIONS AND LIMITATIONS
CD7
CD7 None
CD7
C**********************************************************************
CD8
CD8 SPECIAL COMMENTS
CD8
CD8  Requirements from SDN: 10086-RD-2.20-00
CD8    SOFTWARE REQUIREMENTS DOCUMENT (RD) for the
CD8    FEHM Application Version 2.20
CD8
C**********************************************************************
CD9
CD9 REQUIREMENTS TRACEABILITY
CD9
CD9 2.3.2 Heat- and mass-transfer equations
CD9 2.3.3 Noncondensible gas flow equations
CD9 2.4.7 Dual-porosity formulation
CD9
C**********************************************************************
CDA
CDA REFERENCES
CDA
CDA See FEHMN SRS, MMS, and SDD
CDA
C**********************************************************************
CPS
CPS PSEUDOCODE
CPS
CPS BEGIN dualfa
CPS
CPS mmrelblk - deallocate array space not needed here
CPS mmgetblk - allocate space in temporary arrays
CPS
CPS varchk - determine equation of state parameters for first set of...
CPS ... matrix nodes
CPS varchk - determine equation of state parameters for second set of...
CPS ... matrix nodes
CPS
CPS Compute integer parameters needed later
CPS
CPS FOR each node
CPS
CPS   Initialize parameters for computation of mass balance...
CPS   ... calculations at the matrix nodes
CPS
CPS   IF there is no volume in the second set of matrix nodes
CPS     Set parameters accordingly
CPS   ENDIF
CPS
CPS   Compute terms used later in the heat and mass balance...
CPS   ... calculations for the matrix nodes
CPS
CPS   IF there is liquid at this node
CPS     Compute liquid mass balance contributions for...
CPS     ... liquid in the matrix nodes
CPS   ENDIF
CPS
CPS   IF there is vapor at this node
CPS     Compute vapor mass balance contributions for...
CPS     ... vapor in the matrix nodes
CPS   ENDIF
CPS
CPS   Add accumulation term for first and second matrix nodes
CPS
CPS   Compute terms used later in the mass balance...
CPS   ... calculations for the primary nodes
CPS
CPS   IF there is liquid at this node
CPS     Compute liquid mass balance contributions for...
CPS     ... liquid in the primary nodes
CPS   ENDIF
CPS
CPS   IF there is vapor at this node
CPS     Compute vapor mass balance contributions for...
CPS     ... vapor in the primary nodes
CPS   ENDIF
CPS
CPS   Compute terms needed in matrix for computing heat and mass in...
CPS   ... the matrix nodes
CPS
CPS   Save information needed to back-substitute to obtain values in...
CPS   ... matrix nodes
CPS   Add contributions of matrix nodes to derivative terms
CPS   Add contributions of matrix nodes to residual values
CPS
CPS ENDFOR
CPS
CPS mmrelblk - deallocate array space not needed here
CPS mmgetblk - allocate space in temporary arrays
CPS
CPS END dualfa
CPS
C**********************************************************************

      use davidi
      use comji
      use comhi
      use comgi
      use comfi
      use comei
      use comdi
      use comci
      use combi
      use comdti
      use comai
      implicit none

      real*8 sx1d,pvii,phii,swi,dili,dilkb,divi,divkb,axkb,sx2c,pvikb
      real*8 phikb,radi,radkb,fid,fid1,axy,axyd,axyf,vxy,vxyf,pxyh
      real*8 vxyd,tot,frac0,frac1,frac2,alen,area,al0,al1,al2,dist01
      real*8 dist12,sx1d0,sx1d1,sx1d2,coef1,coef2,sx2t,pxy,pvxy,pxyi
      real*8 a1,a2,a3,a4,b1,b2,r1m,r2m,r3m,r1e,r2e,r3e,a11mp,a12mp
      real*8 a21mp,a22mp,a23mp,a32mp,a33mp,a11me,a12me,a21me,a22me
      real*8 a23me,a32me,a33me,a11ep,a12ep,a21ep,a22ep,a23ep,a32ep
      real*8 a33ep,a11ee,a12ee,a21ee,a22ee,a23ee,a32ee,a33ee,dpvti
      real*8 enli,enlkb,envkb,deli,delkb,devkb,delei,delekb,devekb
      real*8 envi,devi,devei,dilpi,dilpkb,dilei,dilekb,divpi,divpkb
      real*8 divei,divekb,thxi,thyi,thzi,ti,athkb,thxkb,dlpi,dvpi
      real*8 dvpkb,dvaekb,dvaei,dlpkb,dlei,dvei,dlekb,dvekb,dlapi
      real*8 dlapkb,dlaei,dlaekb,dvapi,dvapkb,t11,t12,t21,t22,u11
      real*8 u12,u21,u22,v11,v12,v21,v22,ab22mp,ab22me,ab22ep,ab22ee
      real*8 w11,w12,w21,w22,x11,x12,x21,x22,y11,y12,y21,y22
      real*8 rb2m,rb2e,ali,alm
      integer id,idg,jm,neqp1,neighc,icode,neq2,idum,idl,idll,i
      integer kb,iz,kz,iwa,isl

c**** linear algebra ****
      alm(a1,a2,b1,b2)    =  a1*b1+a2*b2
      ali(a1,a2,a3,a4,b1) =  b1/( a1*a4-a2*a3 )
      neq2   =  neq+neq
      neqp1  =  neq+1
c
c**** call thermo properties of matrix material ****
c
      call  varchk  ( 0,neq   )
      call  varchk  ( 0,neq*2 )

      iwa    =  istrw(neq)+neigh
      neq2   =  neq+neq
      neqp1  =  neq+1
c
c     modify geometric terms
c     save nelm terms for later use
c
c     loop on nodes
      do id=1,neq
c     identify diagonal member
         idg =  nelmdg(id)
         idum =  idg-neqp1
c
c     form equations at id+neq
c
         idl =  id+neq
         idll =  id+neq2
         r1m=0.0
         r2m=0.0
         r3m=0.0
         r1e=0.0
         r2e=0.0
         r3e=0.0
         a11mp=0.0
         a12mp=0.0
         a21mp=0.0
         a22mp=0.0
         a23mp=0.0
         a32mp=0.0
         a33mp=0.0
         a11me=0.0
         a12me=0.0
         a21me=0.0
         a22me=0.0
         a23me=0.0
         a32me=0.0
         a33me=0.0
         a11ep=0.0
         a12ep=0.0
         a21ep=0.0
         a22ep=0.0
         a23ep=0.0
         a32ep=0.0
         a33ep=0.0
         a11ee=0.0
         a12ee=0.0
         a21ee=0.0
         a22ee=0.0
         a23ee=0.0
         a32ee=0.0
         a33ee=0.0
c
c     generate transport terms
c
         tot =  sx1(id)+sx1(idl)+sx1(idll)
         frac0 =  sx1(id)/tot
         frac1 =  sx1(idl)/tot
         frac2 =  sx1(idll)/tot
         alen =  apuv1(id)
         area=tot/alen
         al0 =  frac0*alen
         al1 =  frac1*alen
         al2 =  frac2*alen
         dist12=0.5*(al1+al2)
         coef2=-area/dist12
         sx1d0=sx1(id)
         sx1d1=sx1(idl)
         sx1d2=sx1(idll)
         if(frac2.le.tmch) then
            coef2=0.0
            sx1d2=1.0
         endif

         i=idl
         kb=idll
         iz=id
         kz=id

         pvii=phi(i)
         phii=pvii-pcp(i)
         dpvti=dpcef(i)
         enli=enlf(i)
         deli=delf(i)
         delei=delef(i)
         envi=envf(i)
         devi=devf(i)
         devei=devef(i)
         dili=dil(i)
         divi=div(i)
         dilpi=dilp(i)
         dilei=dile(i)
         divpi=divp(i)
         divei=dive(i)
         thxi=thx(i)
         thyi=thy(i)
         thzi=thz(i)
         ti=t(i)
         swi=s(i)
c
c     form constants for i>neq
c
         axkb =  max( pnx(idll  ),pny(idll  ),
     *        pnz(idll  ),zero_t )
         athkb =  max( thx(idll  ),thy(idll  ),
     *        thz(idll  ),zero_t )
c
c     2-d geometry
c
         if(icnl.ne.0) then
            radi=cord(iz,3)
         else
            radi = 1.
         endif
         jm=1
         neighc=1
         perml(1)=axkb
         permv(1)=perml(1)
         radkb=radi
         sx2c=radkb*coef2
         thxkb=athkb
         sx2t=sx2c*thxkb
         pvikb=phi(kb)
         phikb=pvikb-pcp(kb)
         pxy=sx2c*perml(1)
         pvxy=sx2c*permv(1)
         pxyi=pxy*(phikb-phii)
         pxyh=pvxy*(pvikb-pvii)
         t1(neighc)=pxyi
         t2(neighc)=pxyh
         t3(neighc)=pxy
         t4(neighc)=pvxy

c     liquid phase calculations
c
         jm=1
         neighc=1
         pxyi=t1(neighc)
         axyd=pxyi
         t8(neighc)=axyd
c
c     determine upwind nodes and if liquid phase exists
c
         isl=0
         jm=1
         neighc=1
         fid=.5
         axyd=t8(neighc)
         if(axyd.lt.0.0) fid=dnwgt
         if(axyd.gt.0.0) fid=upwgt
         if(t3(neighc).lt.0.0) t9(neighc)=fid
         if(t3(neighc).gt.0.0) t9(neighc)=fid
         if(swi+s(kb).ne.0.0) isl=1
c
c     form equations
c
         if(isl.ne.0) then
            jm=1
            neighc=1
            axyd=t8(neighc)
            fid=t9(neighc)
            fid1=1.0-fid
            pxyi=t1(neighc)
            pxy=t3(neighc)
            dilkb=dil(kb)
            dilpkb=dilp(kb)
            dilekb=dile(kb)
            enlkb=enlf(kb)
            delkb=delf(kb)
            delekb=delef(kb)
            dilekb=dile(kb)
            dlpi=-pxy
            dlpkb=pxy
            dlei=pxy*dpvti
            dlekb=-pxy*dpcef(kb)

            axyf=(fid*dilkb+fid1*dili)
            axy=axyd*axyf
            dlapi=dlpi*axyf+axyd*fid1*dilpi
            dlapkb=dlpkb*axyf+axyd*fid*dilpkb
            dlaei=dlei*axyf+axyd*fid1*dilei
            dlaekb=dlekb*axyf+axyd*fid*dilekb

            r2m=r2m+axy
            r3m=r3m-axy
            a22mp=a22mp+dlapi
            a22me=a22me+dlaei
            a32mp=a32mp-dlapi
            a32me=a32me-dlaei
            a23mp=a23mp+dlapkb
            a23me=a23me+dlaekb
            a33mp=a33mp-dlapkb
            a33me=a33me-dlaekb
         endif
c
c     vapour phase calculations
c
         jm=1
         neighc=1
         pxyh=t2(neighc)
         vxyd=pxyh
         t8(neighc)=vxyd
c
c     determine upwind nodes and if vapour phase exists
c
c     isl=0
         jm=1
         neighc=1
         fid=.5
         vxyd=t8(neighc)
         if(vxyd.lt.0.0) fid=dnwgt
         if(vxyd.gt.0.0) fid=upwgt
         if(t4(neighc).lt.0.0) t9(neighc)=fid
         if(t4(neighc).gt.0.0) t9(neighc)=fid
c     if(swi+s(kb).ne.2.0) isl=1
c
c     form equations
c
         if(isl.ne.0) then
            jm=1
            neighc=1
            fid=t9(neighc)
            fid1=1.0-fid
            pxyh=t2(neighc)
            pvxy=t4(neighc)
            vxyd=t8(neighc)
            divkb=div(kb)
            divpkb=divp(kb)
            divekb=dive(kb)
            envkb=envf(kb)
            devkb=devf(kb)
            devekb=devef(kb)
            divekb=dive(kb)
            dvpi=-pvxy
            dvpkb=pvxy
            dvei=0.
            dvekb=0.
            vxyf=(fid*divkb+fid1*divi)
            vxy=vxyd*vxyf
            dvapi=dvpi*vxyf+vxyd*fid1*divpi
            dvapkb=dvpkb*vxyf+vxyd*fid*divpkb
            dvaei=dvei*vxyf+vxyd*fid1*divei
            dvaekb=dvekb*vxyf+vxyd*fid*divekb
            r2e=r2e+vxy
            r3e=r3e-vxy
            a22ep=a22ep+dvapi
            a22ee=a22ee+dvaei
            a32ep=a32ep-dvapi
            a32ee=a32ee-dvaei
            a23ep=a23ep+dvapkb
            a23ee=a23ee+dvaekb
            a33ep=a33ep-dvapkb
            a33ee=a33ee-dvaekb
         endif
c
c     add heat conduction

c
c     add accumulation terms for second and first matrix levels
c
         sx1d=sx1d1
c     r2m=r2m+sx1d*(aw*deni(i)+ay*denj(i))+sk(i)
c     r2e=r2e+sx1d*(aw*denei(i)+ay*denej(i))+qh(i)
      r2m=r2m+sx1d*deni(i)+sk(i)
      r2e=r2e+sx1d*denei(i)+qh(i)
c     a22mp=a22mp+sx1d*(aw*dmpf(i))+dq(i)
c     a22me=a22me+sx1d*(aw*dmef(i))+dqh(i)
c     a22ep=a22ep+sx1d*(depf(i)*aw)+drc(i)
c     a22ee=a22ee+sx1d*(aw*deef(i))+deqh(i)
      a22mp=a22mp+sx1d*dmpf(i)+dq(i)
      a22me=a22me+sx1d*dmef(i)+dqh(i)
      a22ep=a22ep+sx1d*depf(i)+drc(i)
      a22ee=a22ee+sx1d*deef(i)+deqh(i)

         sx1d=sx1d2
c     r3m=r3m+sx1d*(aw*deni(kb)+ay*denj(kb))+sk(kb)
c     r3e=r3e+sx1d*(aw*denei(kb)+ay*denej(kb))+qh(kb)
      r3m=r3m+sx1d*deni(kb)+sk(kb)
      r3e=r3e+sx1d*denei(kb)+qh(kb)
c     a33mp=a33mp+sx1d*(aw*dmpf(kb))+dq(kb)
c     a33me=a33me+sx1d*(aw*dmef(kb))+dqh(kb)
c     a33ep=a33ep+sx1d*(depf(kb)*aw)+drc(kb)
c     a33ee=a33ee+sx1d*(aw*deef(kb))+deqh(kb)
      a33mp=a33mp+sx1d*dmpf(kb)+dq(kb)
      a33me=a33me+sx1d*dmef(kb)+dqh(kb)
      a33ep=a33ep+sx1d*depf(kb)+drc(kb)
      a33ee=a33ee+sx1d*deef(kb)+deqh(kb)
c
c     form equations at node id
c
c     geometry
c
         dist01=0.5*(al0+al1)
         coef1=-area/dist01
         if(frac1.le.0.0) then
            coef1=0.0
         endif

         i=id
         kb=idl
         iz=id
         kz=id

         pvii=phi(i)
         phii=pvii-pcp(i)
         dpvti=dpcef(i)
         enli=enlf(i)
         deli=delf(i)
         delei=delef(i)
         envi=envf(i)
         devi=devf(i)
         devei=devef(i)
         dili=dil(i)
         divi=div(i)
         dilpi=dilp(i)
         dilei=dile(i)
         divpi=divp(i)
         divei=dive(i)
         thxi=thx(i)
         thyi=thy(i)
         thzi=thz(i)
         ti=t(i)
         swi=s(i)
c
c     form constants for i>neq
c
         axkb =  max( pnx(idl  ),pny(idl  ),
     *        pnz(idl  ),zero_t )
         athkb =  max( thx(idl  ),thy(idl  ),
     *        thz(idl  ),zero_t )
c
c     2-d geometry
c
         if(icnl.ne.0) then
            radi=cord(iz,3)
         else
            radi = 1.
         endif
         jm=1
         neighc=1
         perml(1)=axkb
         permv(1)=perml(1)
         radkb=radi
         sx2c=radkb*coef1
         thxkb=athkb
         sx2t=sx2c*thxkb
         pvikb=phi(kb)
         phikb=pvikb-pcp(kb)
         pxy=sx2c*perml(1)
         pvxy=sx2c*permv(1)
         pxyi=pxy*(phikb-phii)
         pxyh=pvxy*(pvikb-pvii)
         t1(neighc)=pxyi
         t2(neighc)=pxyh
         t3(neighc)=pxy
         t4(neighc)=pvxy

c     liquid phase calculations
         jm=1
         neighc=1
         pxyi=t1(neighc)
         axyd=pxyi
         t8(neighc)=axyd

c     determine upwind nodes and if liquid phase exists
         isl=0
         jm=1
         neighc=1
         fid=.5
         axyd=t8(neighc)
         if(axyd.lt.0.0) fid=dnwgt
         if(axyd.gt.0.0) fid=upwgt
         if(t3(neighc).lt.0.0) t9(neighc)=fid
         if(t3(neighc).gt.0.0) t9(neighc)=fid
         if(swi+s(kb).ne.0.0) isl=1

c     form equations
         if(isl.ne.0) then
            jm=1
            neighc=1
            axyd=t8(neighc)
            fid=t9(neighc)
            fid1=1.0-fid
            pxyi=t1(neighc)
            pxy=t3(neighc)
            dilkb=dil(kb)
            dilpkb=dilp(kb)
            dilekb=dile(kb)
            enlkb=enlf(kb)
            delkb=delf(kb)
            delekb=delef(kb)
            dilekb=dile(kb)
            dlpi=-pxy
            dlpkb=pxy
            dlei=pxy*dpvti
            dlekb=-pxy*dpcef(kb)
            axyf=(fid*dilkb+fid1*dili)
            axy=axyd*axyf
            dlapi=dlpi*axyf+axyd*fid1*dilpi
            dlapkb=dlpkb*axyf+axyd*fid*dilpkb
            dlaei=dlei*axyf+axyd*fid1*dilei
            dlaekb=dlekb*axyf+axyd*fid*dilekb

            r1m=r1m+axy
            r2m=r2m-axy
            a11mp=a11mp+dlapi
            a11me=a11me+dlaei
            a21mp=a21mp-dlapi
            a21me=a21me-dlaei
            a12mp=a12mp+dlapkb
            a12me=a12me+dlaekb
            a22mp=a22mp-dlapkb
            a22me=a22me-dlaekb
         endif
c
c     vapour phase calculations
c
         jm=1
         neighc=1
         pxyh=t2(neighc)
         vxyd=pxyh
         t8(neighc)=vxyd
c
c     determine upwind nodes and if vapour phase exists
c
c     isl=0
         jm=1

         neighc=1
         fid=.5
         vxyd=t8(neighc)
         if(vxyd.lt.0.0) fid=dnwgt
         if(vxyd.gt.0.0) fid=upwgt
         if(t4(neighc).lt.0.0) t9(neighc)=fid
         if(t4(neighc).gt.0.0) t9(neighc)=fid
c
c     form equations
c
         if(isl.ne.0) then
            jm=1
            neighc=1
            fid=t9(neighc)
            fid1=1.0-fid
            pxyh=t2(neighc)
            pvxy=t4(neighc)
            vxyd=t8(neighc)
            divkb=div(kb)
            divpkb=divp(kb)
            divekb=dive(kb)
            envkb=envf(kb)
            devkb=devf(kb)
            devekb=devef(kb)
            divekb=dive(kb)
            dvpi=-pvxy
            dvpkb=pvxy
            dvei=0.
            dvekb=0.
            vxyf=(fid*divkb+fid1*divi)
            vxy=vxyd*vxyf
            dvapi=dvpi*vxyf+vxyd*fid1*divpi
            dvapkb=dvpkb*vxyf+vxyd*fid*divpkb
            dvaei=dvei*vxyf+vxyd*fid1*divei
            dvaekb=dvekb*vxyf+vxyd*fid*divekb
            r1e=r1e+vxy
            r2e=r2e-vxy
            a11ep=a11ep+dvapi
            a11ee=a11ee+dvaei
            a21ep=a21ep-dvapi
            a21ee=a21ee-dvaei
            a12ep=a12ep+dvapkb
            a12ee=a12ee+dvaekb
            a22ep=a22ep-dvapkb
            a22ee=a22ee-dvaekb
         endif
c
c     add heat conduction
c

c     prefactor dual porosity solution in
c
c     form matrix ab22=a22-a23*a33i*a32
c
c     form a33i=t
         t11 =  ali(a33mp,a33me,a33ep,a33ee,a33ee)
         t12 =  -ali(a33mp,a33me,a33ep,a33ee,a33me)
         t21 =  -ali(a33mp,a33me,a33ep,a33ee,a33ep)
         t22 =  ali(a33mp,a33me,a33ep,a33ee,a33mp)
c     form a23*a33i =  u
         u11 =  alm(a23mp,a23me,t11,t21)
         u12 =  alm(a23mp,a23me,t12,t22)
         u21 =  alm(a23ep,a23ee,t11,t21)
         u22 =  alm(a23ep,a23ee,t12,t22)
c     form a23*a33i*a32=v
         v11 =  alm(u11,u12,a32mp,a32ep)
         v12 =  alm(u11,u12,a32me,a32ee)
         v21 =  alm(u21,u22,a32mp,a32ep)
         v22 =  alm(u21,u22,a32me,a32ee)
c     form a22-a23*a33i*a32=ab22
         ab22mp =  a22mp-v11
         ab22me =  a22me-v12
         ab22ep =  a22ep-v21
         ab22ee =  a22ee-v22
c     form rb2=r2-a23*a33i*r3
         rb2m =  r2m-alm(u11,u12,r3m,r3e)
         rb2e =  r2e-alm(u21,u22,r3m,r3e)
c
c     form matrix ab11=a11-a12*ab22i*a21
c
c     form ab22i =  w
         w11 =   ali(ab22mp,ab22me,ab22ep,ab22ee,ab22ee)
         w12 =  -ali(ab22mp,ab22me,ab22ep,ab22ee,ab22me)
         w21 =  -ali(ab22mp,ab22me,ab22ep,ab22ee,ab22ep)
         w22 =   ali(ab22mp,ab22me,ab22ep,ab22ee,ab22mp)
c     form a12*ab22i=x
         x11 =  alm(a12mp,a12me,w11,w21)
         x12 =  alm(a12mp,a12me,w12,w22)
         x21 =  alm(a12ep,a12ee,w11,w21)
         x22 =  alm(a12ep,a12ee,w12,w22)
c     form a12*ab22i*a21=y
         y11 =  alm(x11,x12,a21mp,a21ep)
         y12 =  alm(x11,x12,a21me,a21ee)
         y21 =  alm(x21,x22,a21mp,a21ep)
         y22 =  alm(x21,x22,a21me,a21ee)
c     form a11-a12*a22i*a12
         a(idum+nmat(1)) =  a(idum+nmat(1))-y11+a11mp
         a(idum+nmat(2)) =  a(idum+nmat(2))-y12+a11me
         a(idum+nmat(3)) =  a(idum+nmat(3))-y21+a11ep
         a(idum+nmat(4)) =  a(idum+nmat(4))-y22+a11ee
c     form rb1=r1-a12*ab22i*rb2
         bp(id+nrhs(1)) =  bp(id+nrhs(1))+r1m-alm(x11,x12,rb2m,rb2e)
         bp(id+nrhs(2)) =  bp(id+nrhs(2))+r1e-alm(x21,x22,rb2m,rb2e)
c
c     save matrices to back out solution
c     ab22i=wb,a33i=tb,a21,a32
c
         wb11(id) =  w11
         wb12(id) =  w12
         wb21(id) =  w21
         wb22(id) =  w22
         tb11(id) =  t11
         tb12(id) =  t12
         tb21(id) =  t21
         tb22(id) =  t22
         a21mpf(id) =  a21mp
         a21mef(id) =  a21me
         a21epf(id) =  a21ep
         a21eef(id) =  a21ee
         a32mpf(id) =  a32mp
         a32mef(id) =  a32me
         a32epf(id) =  a32ep
         a32eef(id) =  a32ee
c
c     save vectors neccessary to get solution
c
         rb2mf(id) =  rb2m
         rb2ef(id) =  rb2e
         r3mf(id) =  r3m
         r3ef(id) =  r3e
      enddo

      return
      end

max_stars_repo_name: satkarra/FEHM

max_issues_repo_name: satkarra/FEHM

max_forks_repo_name: satkarra/FEHM


text:
subroutine uv2zeta(mmax        ,nmaxus     ,kcs        ,divByCellWidth   ,ucompu   , &
                 & vcompv      ,ucompz     ,vcompz     ,guu              ,gvv      ,gdp    )
!----- GPL ---------------------------------------------------------------------
!
!  Copyright (C)  Stichting Deltares, 2011-2016.
!
!  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 3.
!
!  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, "Delft3D" and "Deltares"
!  are registered trademarks of Stichting Deltares, and remain the property of
!  Stichting Deltares. All rights reserved.
!
!-------------------------------------------------------------------------------
!  $Id: uv2zeta.f90 5717 2016-01-12 11:35:24Z mourits $
!  $HeadURL: https://svn.oss.deltares.nl/repos/delft3d/tags/6686/src/engines_gpl/flow2d3d/packages/io/src/output/uv2zeta.f90 $
!!--description-----------------------------------------------------------------
!
!    Function: - Averaging of the parameters on which the Fourier Analysis is performed
!                from u- and v-points to zeta-points.
!                In addition, for discharges, the value is scaled to a unit discharge
!                [m3/m]. This is necessary for determining the vector magnitude.
!
!!--pseudo code and references--------------------------------------------------
! NONE
!!--declarations----------------------------------------------------------------
    use precision
    use globaldata
    !
    implicit none
    !
    type(globdat),target :: gdp
    !
    ! The following list of pointer parameters is used to point inside the gdp structure
    !
!
! Global variables
!
    integer  , dimension(gdp%d%nlb:gdp%d%nub, gdp%d%mlb:gdp%d%mub) , intent(in)  :: kcs             !!  Description and declaration in esm_alloc_int.f90
    integer                                                        , intent(in)  :: mmax            !!  Description and declaration in esm_alloc_int.f90
    integer                                                        , intent(in)  :: nmaxus          !!  Description and declaration in esm_alloc_int.f90
    real(fp) , dimension(gdp%d%nlb:gdp%d%nub, gdp%d%mlb:gdp%d%mub) , intent(in)  :: ucompu          !  First component of vector (u, qxk or taubpu) in u-points
    real(fp) , dimension(gdp%d%nlb:gdp%d%nub, gdp%d%mlb:gdp%d%mub) , intent(in)  :: vcompv          !  Second component of vector (v, qyk or taubpv) in v-points
    real(fp) , dimension(gdp%d%nlb:gdp%d%nub, gdp%d%mlb:gdp%d%mub) , intent(out) :: ucompz          !  First component of vector (u, qxk or taubpu) in zeta-points
    real(fp) , dimension(gdp%d%nlb:gdp%d%nub, gdp%d%mlb:gdp%d%mub) , intent(out) :: vcompz          !  Second component of vector (v, qyk or taubpv) in zeta-points
    real(fp) , dimension(gdp%d%nlb:gdp%d%nub, gdp%d%mlb:gdp%d%mub) , intent(in)  :: guu             !  Distance between depth points over u points
    real(fp) , dimension(gdp%d%nlb:gdp%d%nub, gdp%d%mlb:gdp%d%mub) , intent(in)  :: gvv             !  Distance between depth points over v points
    logical                                                        , intent(in)  :: divByCellWidth  !  is TRUE when parameter has be scaled with cell width to obtain the correct dimensions

!
! Local variables
!
    integer         :: m       ! Loop counter over MMAX
    integer         :: md      ! m-1
    integer         :: n       ! Loop counter over NMAXUS
    integer         :: nd      ! n-1
!
!! executable statements -------------------------------------------------------
!
    !
    ! Average the fourier component from u and v points to the centre of the cell (zeta-point)
    !
    do n = 1, nmaxus
       do m = 1, mmax
          if (kcs(n, m) == 1) then
             md = max(1, m - 1)
             nd = max(1, n - 1)
             if (divByCellWidth) then
                !
                ! Scale the parameter with the cell width, to obtain correct dimensions
                !
                ucompz(n,m) = 0.5_fp * ( ucompu(n,md)/guu(n, md) + ucompu(n,m)/guu(n, m) )
                vcompz(n,m) = 0.5_fp * ( vcompv(nd,m)/gvv(nd, m) + vcompv(n,m)/gvv(n, m) )
             else
                ucompz(n,m) = 0.5_fp * ( ucompu(n,md) + ucompu(n,m) )
                vcompz(n,m) = 0.5_fp * ( vcompv(nd,m) + vcompv(n,m) )
             endif
          endif
       enddo
    enddo
    !
end subroutine uv2zeta

max_stars_repo_name: liujiamingustc/phd

max_issues_repo_name: liujiamingustc/phd

max_forks_repo_name: liujiamingustc/phd


text:
*DECK DRJ
      DOUBLE PRECISION FUNCTION DRJ (X, Y, Z, P, IER)
C***BEGIN PROLOGUE  DRJ
C***PURPOSE  Compute the incomplete or complete (X or Y or Z is zero)
C            elliptic integral of the 3rd kind.  For X, Y, and Z non-
C            negative, at most one of them zero, and P positive,
C             RJ(X,Y,Z,P) = Integral from zero to infinity of
C                              -1/2     -1/2     -1/2     -1
C                    (3/2)(t+X)    (t+Y)    (t+Z)    (t+P)  dt.
C***LIBRARY   SLATEC
C***CATEGORY  C14
C***TYPE      DOUBLE PRECISION (RJ-S, DRJ-D)
C***KEYWORDS  COMPLETE ELLIPTIC INTEGRAL, DUPLICATION THEOREM,
C             INCOMPLETE ELLIPTIC INTEGRAL, INTEGRAL OF THE THIRD KIND,
C             TAYLOR SERIES
C***AUTHOR  Carlson, B. C.
C             Ames Laboratory-DOE
C             Iowa State University
C             Ames, IA  50011
C           Notis, E. M.
C             Ames Laboratory-DOE
C             Iowa State University
C             Ames, IA  50011
C           Pexton, R. L.
C             Lawrence Livermore National Laboratory
C             Livermore, CA  94550
C***DESCRIPTION
C
C   1.     DRJ
C          Standard FORTRAN function routine
C          Double precision version
C          The routine calculates an approximation result to
C          DRJ(X,Y,Z,P) = Integral from zero to infinity of
C
C                                -1/2     -1/2     -1/2     -1
C                      (3/2)(t+X)    (t+Y)    (t+Z)    (t+P)  dt,
C
C          where X, Y, and Z are nonnegative, at most one of them is
C          zero, and P is positive.  If X or Y or Z is zero, the
C          integral is COMPLETE.  The duplication theorem is iterated
C          until the variables are nearly equal, and the function is
C          then expanded in Taylor series to fifth order.
C
C
C   2.     Calling Sequence
C          DRJ( X, Y, Z, P, IER )
C
C          Parameters on Entry
C          Values assigned by the calling routine
C
C          X      - Double precision, nonnegative variable
C
C          Y      - Double precision, nonnegative variable
C
C          Z      - Double precision, nonnegative variable
C
C          P      - Double precision, positive variable
C
C
C          On  Return    (values assigned by the DRJ routine)
C
C          DRJ     - Double precision approximation to the integral
C
C          IER    - Integer
C
C                   IER = 0 Normal and reliable termination of the
C                           routine. It is assumed that the requested
C                           accuracy has been achieved.
C
C                   IER >  0 Abnormal termination of the routine
C
C
C          X, Y, Z, P are unaltered.
C
C
C   3.    Error Messages
C
C         Value of IER assigned by the DRJ routine
C
C              Value assigned         Error Message printed
C              IER = 1                MIN(X,Y,Z) .LT. 0.0D0
C                  = 2                MIN(X+Y,X+Z,Y+Z,P) .LT. LOLIM
C                  = 3                MAX(X,Y,Z,P) .GT. UPLIM
C
C
C
C   4.     Control Parameters
C
C                  Values of LOLIM, UPLIM, and ERRTOL are set by the
C                  routine.
C
C
C          LOLIM and UPLIM determine the valid range of X, Y, Z, and P
C
C          LOLIM is not less than the cube root of the value
C          of LOLIM used in the routine for DRC.
C
C          UPLIM is not greater than 0.3 times the cube root of
C          the value of UPLIM used in the routine for DRC.
C
C
C                     Acceptable values for:   LOLIM      UPLIM
C                     IBM 360/370 SERIES   :   2.0D-26     3.0D+24
C                     CDC 6000/7000 SERIES :   5.0D-98     3.0D+106
C                     UNIVAC 1100 SERIES   :   5.0D-103    6.0D+101
C                     CRAY                 :   1.32D-822   1.4D+821
C                     VAX 11 SERIES        :   2.5D-13     9.0D+11
C
C
C
C          ERRTOL determines the accuracy of the answer
C
C                 the value assigned by the routine will result
C                 in solution precision within 1-2 decimals of
C                 "machine precision".
C
C
C
C
C          Relative error due to truncation of the series for DRJ
C          is less than 3 * ERRTOL ** 6 / (1 - ERRTOL) ** 3/2.
C
C
C
C        The accuracy of the computed approximation to the integral
C        can be controlled by choosing the value of ERRTOL.
C        Truncation of a Taylor series after terms of fifth order
C        introduces an error less than the amount shown in the
C        second column of the following table for each value of
C        ERRTOL in the first column.  In addition to the truncation
C        error there will be round-off error, but in practice the
C        total error from both sources is usually less than the
C        amount given in the table.
C
C
C
C          Sample choices:  ERRTOL   Relative truncation
C                                    error less than
C                           1.0D-3    4.0D-18
C                           3.0D-3    3.0D-15
C                           1.0D-2    4.0D-12
C                           3.0D-2    3.0D-9
C                           1.0D-1    4.0D-6
C
C                    Decreasing ERRTOL by a factor of 10 yields six more
C                    decimal digits of accuracy at the expense of one or
C                    two more iterations of the duplication theorem.
C
C *Long Description:
C
C   DRJ Special Comments
C
C
C     Check by addition theorem: DRJ(X,X+Z,X+W,X+P)
C     + DRJ(Y,Y+Z,Y+W,Y+P) + (A-B) * DRJ(A,B,B,A) + 3.0D0 / SQRT(A)
C     = DRJ(0,Z,W,P), where X,Y,Z,W,P are positive and X * Y
C     = Z * W,  A = P * P * (X+Y+Z+W),  B = P * (P+X) * (P+Y),
C     and B - A = P * (P-Z) * (P-W).  The sum of the third and
C     fourth terms on the left side is 3.0D0 * DRC(A,B).
C
C
C          On Input:
C
C     X, Y, Z, and P are the variables in the integral DRJ(X,Y,Z,P).
C
C
C          On Output:
C
C
C          X, Y, Z, P are unaltered.
C
C          ********************************************************
C
C          WARNING: Changes in the program may improve speed at the
C                   expense of robustness.
C
C    -------------------------------------------------------------------
C
C
C   Special double precision functions via DRJ and DRF
C
C
C                  Legendre form of ELLIPTIC INTEGRAL of 3rd kind
C                  -----------------------------------------
C
C
C                          PHI         2         -1
C             P(PHI,K,N) = INT (1+N SIN (THETA) )   *
C                           0
C
C
C                                  2    2         -1/2
C                             *(1-K  SIN (THETA) )     D THETA
C
C
C                                           2          2   2
C                        = SIN (PHI) DRF(COS (PHI), 1-K SIN (PHI),1)
C
C                                   3             2         2   2
C                         -(N/3) SIN (PHI) DRJ(COS (PHI),1-K SIN (PHI),
C
C                                  2
C                         1,1+N SIN (PHI))
C
C
C
C                  Bulirsch form of ELLIPTIC INTEGRAL of 3rd kind
C                  -----------------------------------------
C
C
C                                            2 2    2
C                  EL3(X,KC,P) = X DRF(1,1+KC X ,1+X ) +
C
C                                            3           2 2    2     2
C                               +(1/3)(1-P) X  DRJ(1,1+KC X ,1+X ,1+PX )
C
C
C                                           2
C                  CEL(KC,P,A,B) = A RF(0,KC ,1) +
C
C
C                                                      2
C                                 +(1/3)(B-PA) DRJ(0,KC ,1,P)
C
C
C                  Heuman's LAMBDA function
C                  -----------------------------------------
C
C
C                                2                      2      2    1/2
C                  L(A,B,P) =(COS (A)SIN(B)COS(B)/(1-COS (A)SIN (B))   )
C
C                                            2         2       2
C                            *(SIN(P) DRF(COS (P),1-SIN (A) SIN (P),1)
C
C                                 2       3            2       2
C                            +(SIN (A) SIN (P)/(3(1-COS (A) SIN (B))))
C
C                                    2         2       2
C                            *DRJ(COS (P),1-SIN (A) SIN (P),1,1-
C
C                                2       2          2       2
C                            -SIN (A) SIN (P)/(1-COS (A) SIN (B))))
C
C
C
C                  (PI/2) LAMBDA0(A,B) =L(A,B,PI/2) =
C
C                   2                         2       2    -1/2
C              = COS (A)  SIN(B) COS(B) (1-COS (A) SIN (B))
C
C                           2                  2       2
C                 *DRF(0,COS (A),1) + (1/3) SIN (A) COS (A)
C
C                                      2       2    -3/2
C                 *SIN(B) COS(B) (1-COS (A) SIN (B))
C
C                           2         2       2          2       2
C                 *DRJ(0,COS (A),1,COS (A) COS (B)/(1-COS (A) SIN (B)))
C
C
C                  Jacobi ZETA function
C                  -----------------------------------------
C
C                        2                     2   2    1/2
C             Z(B,K) = (K/3) SIN(B) COS(B) (1-K SIN (B))
C
C
C                                  2      2   2                 2
C                        *DRJ(0,1-K ,1,1-K SIN (B)) / DRF (0,1-K ,1)
C
C
C  ---------------------------------------------------------------------
C
C***REFERENCES  B. C. Carlson and E. M. Notis, Algorithms for incomplete
C                 elliptic integrals, ACM Transactions on Mathematical
C                 Software 7, 3 (September 1981), pp. 398-403.
C               B. C. Carlson, Computing elliptic integrals by
C                 duplication, Numerische Mathematik 33, (1979),
C                 pp. 1-16.
C               B. C. Carlson, Elliptic integrals of the first kind,
C                 SIAM Journal of Mathematical Analysis 8, (1977),
C                 pp. 231-242.
C***ROUTINES CALLED  D1MACH, DRC, XERMSG
C***REVISION HISTORY  (YYMMDD)
C   790801  DATE WRITTEN
C   890531  Changed all specific intrinsics to generic.  (WRB)
C   891009  Removed unreferenced statement labels.  (WRB)
C   891009  REVISION DATE from Version 3.2
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
C   900326  Removed duplicate information from DESCRIPTION section.
C           (WRB)
C   900510  Changed calls to XERMSG to standard form, and some
C           editorial changes.  (RWC)).
C   920501  Reformatted the REFERENCES section.  (WRB)
C***END PROLOGUE  DRJ
      INTEGER IER
      CHARACTER*16 XERN3, XERN4, XERN5, XERN6, XERN7
      DOUBLE PRECISION ALFA, BETA, C1, C2, C3, C4, EA, EB, EC, E2, E3
      DOUBLE PRECISION LOLIM, UPLIM, EPSLON, ERRTOL, D1MACH
      DOUBLE PRECISION LAMDA, MU, P, PN, PNDEV
      DOUBLE PRECISION POWER4, DRC, SIGMA, S1, S2, S3, X, XN, XNDEV
      DOUBLE PRECISION XNROOT, Y, YN, YNDEV, YNROOT, Z, ZN, ZNDEV,
     * ZNROOT
      LOGICAL FIRST
      SAVE ERRTOL,LOLIM,UPLIM,C1,C2,C3,C4,FIRST
      DATA FIRST /.TRUE./
C
C***FIRST EXECUTABLE STATEMENT  DRJ
      IF (FIRST) THEN
         ERRTOL = (D1MACH(3)/3.0D0)**(1.0D0/6.0D0)
         LOLIM  = (5.0D0 * D1MACH(1))**(1.0D0/3.0D0)
         UPLIM  = 0.30D0*( D1MACH(2) / 5.0D0)**(1.0D0/3.0D0)
C
         C1 = 3.0D0/14.0D0
         C2 = 1.0D0/3.0D0
         C3 = 3.0D0/22.0D0
         C4 = 3.0D0/26.0D0
      ENDIF
      FIRST = .FALSE.
C
C         CALL ERROR HANDLER IF NECESSARY.
C
      DRJ = 0.0D0
      IF (MIN(X,Y,Z).LT.0.0D0) THEN
         IER = 1
         WRITE (XERN3, '(1PE15.6)') X
         WRITE (XERN4, '(1PE15.6)') Y
         WRITE (XERN5, '(1PE15.6)') Z
         CALL XERMSG ('SLATEC', 'DRJ',
     *      'MIN(X,Y,Z).LT.0 WHERE X = ' // XERN3 // ' Y = ' // XERN4 //
     *      ' AND Z = ' // XERN5, 1, 1)
         RETURN
      ENDIF
C
      IF (MAX(X,Y,Z,P).GT.UPLIM) THEN
         IER = 3
         WRITE (XERN3, '(1PE15.6)') X
         WRITE (XERN4, '(1PE15.6)') Y
         WRITE (XERN5, '(1PE15.6)') Z
         WRITE (XERN6, '(1PE15.6)') P
         WRITE (XERN7, '(1PE15.6)') UPLIM
         CALL XERMSG ('SLATEC', 'DRJ',
     *      'MAX(X,Y,Z,P).GT.UPLIM WHERE X = ' // XERN3 // ' Y = ' //
     *      XERN4 // ' Z = ' // XERN5 // ' P = ' // XERN6 //
     *      ' AND UPLIM = ' // XERN7, 3, 1)
         RETURN
      ENDIF
C
      IF (MIN(X+Y,X+Z,Y+Z,P).LT.LOLIM) THEN
         IER = 2
         WRITE (XERN3, '(1PE15.6)') X
         WRITE (XERN4, '(1PE15.6)') Y
         WRITE (XERN5, '(1PE15.6)') Z
         WRITE (XERN6, '(1PE15.6)') P
         WRITE (XERN7, '(1PE15.6)') LOLIM
         CALL XERMSG ('SLATEC', 'RJ',
     *      'MIN(X+Y,X+Z,Y+Z,P).LT.LOLIM WHERE X = ' // XERN3 //
     *      ' Y = ' // XERN4 // ' Z = '  // XERN5 // ' P = ' // XERN6 //
     *      ' AND LOLIM = ', 2, 1)
         RETURN
      ENDIF
C
      IER = 0
      XN = X
      YN = Y
      ZN = Z
      PN = P
      SIGMA = 0.0D0
      POWER4 = 1.0D0
C
   30 MU = (XN+YN+ZN+PN+PN)*0.20D0
      XNDEV = (MU-XN)/MU
      YNDEV = (MU-YN)/MU
      ZNDEV = (MU-ZN)/MU
      PNDEV = (MU-PN)/MU
      EPSLON = MAX(ABS(XNDEV), ABS(YNDEV), ABS(ZNDEV), ABS(PNDEV))
      IF (EPSLON.LT.ERRTOL) GO TO 40
      XNROOT =  SQRT(XN)
      YNROOT =  SQRT(YN)
      ZNROOT =  SQRT(ZN)
      LAMDA = XNROOT*(YNROOT+ZNROOT) + YNROOT*ZNROOT
      ALFA = PN*(XNROOT+YNROOT+ZNROOT) + XNROOT*YNROOT*ZNROOT
      ALFA = ALFA*ALFA
      BETA = PN*(PN+LAMDA)*(PN+LAMDA)
      SIGMA = SIGMA + POWER4*DRC(ALFA,BETA,IER)
      POWER4 = POWER4*0.250D0
      XN = (XN+LAMDA)*0.250D0
      YN = (YN+LAMDA)*0.250D0
      ZN = (ZN+LAMDA)*0.250D0
      PN = (PN+LAMDA)*0.250D0
      GO TO 30
C
   40 EA = XNDEV*(YNDEV+ZNDEV) + YNDEV*ZNDEV
      EB = XNDEV*YNDEV*ZNDEV
      EC = PNDEV*PNDEV
      E2 = EA - 3.0D0*EC
      E3 = EB + 2.0D0*PNDEV*(EA-EC)
      S1 = 1.0D0 + E2*(-C1+0.750D0*C3*E2-1.50D0*C4*E3)
      S2 = EB*(0.50D0*C2+PNDEV*(-C3-C3+PNDEV*C4))
      S3 = PNDEV*EA*(C2-PNDEV*C3) - C2*PNDEV*EC
      DRJ = 3.0D0*SIGMA + POWER4*(S1+S2+S3)/(MU* SQRT(MU))
      RETURN
      END


text:
      subroutine genfmt_writefeffdat(ipath, ntext, text,
     $     nleg, deg, reff, rnrmav, edge, rat, iz, ipot,
     $     potlbl, l0, il0, ne, xk, ck, ph, cchi)
c
c  write feffNNNN.dat file for NNNN=ipath
      implicit none

      include 'const.h'
      include 'dim.h'
      include 'vers.h'
      double precision eps
      parameter (eps = 1.0d-16)

      integer ipath, ntext, nleg, l0, il0, ne, it, istrln
      integer ipot(0:legtot), iz(0:npotx), jpot, ileg, ie, ios

      double precision deg, reff, rnrmav, edge, getxk
      double precision rat(3, 0:legtot+1)
      double precision xlam, redfac, phff, phffo, cdelt, cdelto
      double precision xk(nex)

      complex*16 ph(nex, ltot+1, 0:npotx)
      complex*16 cchi(nex), ck(nex), cfms

      character fname*16
      character*80 text(40)
      character*6 potlbl(0:npotx)

      external istrln, getxk

 60   format (1x, a)
 70   format (1x, 79('-'))

      write(fname, 241)  ipath
 241  format ('feff', i4.4, '.dat')
      open (unit=3, file=fname, status='unknown', iostat=ios)
      call chopen (ios, fname, 'genfmt')
c     put header on feff.dat
      do 245  it = 1, ntext
         write(3,60)  text(it)(1:istrln(text(it)))
 245  continue
      write(3,250) ipath, 2, vfeff, vgenfm
 250  format (' Path', i5, '      icalc ', i7, t57, 2a12)
      write(3,70)
      write(3,290)  nleg, deg, reff*bohr, rnrmav, edge*ryd
 290  format (1x, i3, f8.3, f9.4, f10.4, f11.5,
     1     ' nleg, deg, reff, rnrmav(bohr), edge')

c
      write(3,300)
 300  format ('        x         y         z   pot at#')
      jpot = ipot(nleg)
      write(3,310)  rat(1,nleg)*bohr, rat(2,nleg)*bohr,
     $     rat(3,nleg)*bohr, jpot, iz(jpot), potlbl(jpot)
 310  format (1x, 3f10.4, i3, i4, 1x, a6, '   absorbing atom')
 311  format (1x, 3f10.4, i3, i4, 1x, a6)
      do 330  ileg = 1, nleg-1
         jpot = ipot(ileg)
         write(3,311)  rat(1,ileg)*bohr, rat(2,ileg)*bohr,
     $        rat(3,ileg)*bohr, jpot, iz(jpot), potlbl(jpot)
 330  continue

      write(3,340)
 340  format    ('    k   real[2*phc]   mag[feff]  phase[feff]',
     1     ' red factor   lambda      real[p]@#')

c     Make the feff.dat stuff and write it to feff.dat
      do 900  ie = 1, ne
c     Consider chi in the standard XAFS form.  Use R = rtot/2.
         xlam = 1.0e10
         if (abs(dimag(ck(ie))) .gt. eps) xlam = 1/dimag(ck(ie))
         redfac = exp (-2 * dimag (ph(ie,il0,ipot(nleg))))
         cdelt = 2*dble(ph(ie,il0,ipot(nleg)))
         cfms = cchi(ie) * xk(ie) * reff**2 *
     1        exp(2*reff/xlam) / redfac
         if (abs(cchi(ie)) .lt. eps)  then
            phff = 0
         else
            phff = atan2 (dimag(cchi(ie)), dble(cchi(ie)))
         endif
c     remove 2 pi jumps in phases
         if (ie .gt. 1)  then
            call pijump (phff, phffo)
            call pijump (cdelt, cdelto)
         endif
         phffo = phff
         cdelto = cdelt

c     write 1 k, momentum wrt fermi level k=sqrt(p**2-kf**2)
c           2 central atom phase shift (real part),
c           3 magnitude of feff,
c           4 phase of feff,
c           5 absorbing atom reduction factor,
c           6 mean free path = 1/(Im (p))
c           7 real part of local momentum p
         write(3,640)
     1        xk(ie)/bohr,
     2        cdelt + l0*pi,
     3        abs(cfms) * bohr,
     4        phff - cdelt - l0*pi,
     5        redfac,
     6        xlam * bohr,
     7        dble(ck(ie))/bohr
 640     format (1x, f6.3, 1x, 3(1pe11.4,1x),0pe11.4,1x,
     1        2(1pe11.4,1x))
 900  continue

      close(3)
      return
      end

      subroutine calc_zabinsky(ne, ik0, deg, ck, cchi, xportx, crit)
c
c calculate Zabinsky Curved Wave Importance Factor for a Path
c
c  input:  ne    number of energy points
c          ik0
c          deg   degeneracy
c          ck    complex k
c          cchi  complex chi(k)
c          xport max importance factor so far
c  output: xport max importance factor so far
c          crit  importance / max importance

      implicit none
      include 'dim.h'
      complex*16 cchi(nex), ck(nex)
      double precision ckmag(nex), ffmag(nex)
      double precision xport, xportx, deg, crit
      integer ie, ne, ik0, nemax

      do 10  ie = 1, ne
         ckmag(ie) = abs(ck(ie))
         ffmag(ie) = abs(cchi(ie))
 10   continue

c     integrate from edge (ik0) to ne
      nemax = ne - ik0 + 1
      call trap(ckmag(ik0), ffmag(ik0), nemax, xport)
      xport = abs(deg*xport)
      if (xport .gt. xportx)  xportx = xport
      crit = 100 * xport / xportx

      return
      end


program tarefaa
    implicit real(8) (a-h,o-z)
    parameter (pi = 4d0*atan(1d0))

    g = 9.8d0
    a_l = 9.8d0
    a_m = 1d0
    dt = 1d-3
    tmax = 120d0
    omega_i = 0d0
    theta_i = pi/15
    i = 0

    open(10, file='saida-a1-10407962.dat')

    do while ( i*dt.le.tmax )

        theta_i = mod(theta_i, 2*pi)

        omega_j = omega_i - g*theta_i*dt/a_l
        theta_j = theta_i + omega_i*dt

        omega_i = omega_j
        theta_i = theta_j

        Ec = 0.5*a_m*(omega_i*a_l)**2
        Ep = a_m*g*a_l*(1-cos(theta_i))

        write(10, '(F0.16,3(" ",F0.16))') i*dt, theta_i, Ec+Ep
        i = i+1
    end do

    close(10)

    omega_i = 0d0
    theta_i = pi/15
    i = 0

    open(10, file='saida-a2-10407962.dat')

    do while ( i*dt.le.tmax )

        theta_i = mod(theta_i, 2*pi)

        omega_i = omega_i - g*theta_i*dt/a_l
        theta_i = theta_i + omega_i*dt

        Ec = 0.5*a_m*(omega_i*a_l)**2
        Ep = a_m*g*a_l*(1-cos(theta_i))

        write(10, '(F0.16,2(" ",F0.16))') i*dt, theta_i, Ec+Ep
        i = i+1
    end do

    close(10)
end program tarefaa


text: The Jewish Studies program is a minors minor offered at UC Davis with about 25 courses being taught every year.  Most classes in the program focus on the lifestyles and modern social concerns of the Jewish people, but there are also some that discuss the core religious beliefs and historical foundations of Judaism.  As part of the program, students have the option to study the Hebrew language for up to two years.


text:
	SUBROUTINE SNPXAX  ( pname1, pname2, xaxis, iyaxis, xmin, xmax,
     +			     xstrt, xstop, xlbl, nxlbl, ixlbfr, ixglfr,
     +                       ixtmfr, iret )
C************************************************************************
C* SNPXAX								*
C*									*
C* This subroutine defines the x axis to be drawn in SNPROF.		*
C*									*
C* SNPXAX  ( PNAME1, PNAME2, XAXIS, IYAXIS, XMIN, XMAX, XSTRT, XSTOP,	*
C*           XLBL, NXLBL, IXLBFR, IXGLFR, IXTMFR, IRET )		*
C*									*
C* Input parameters:							*
C*	PNAME1		CHAR*		Parameter 1			*
C*	PNAME2		CHAR*		Parameter 2			*
C*	XAXIS		CHAR*		X axis input			*
C*	IYAXIS		INTEGER		Y axis coordinate		*
C*      XMIN            REAL		Minimum x value			*
C*      XMAX        	REAL		Maximum x value			*
C*									*
C* Output parameters:							*
C*	XSTRT		REAL		Left x				*
C*	XSTOP		REAL		Right x				*
C*	XLBL  (NXLBL)	REAL		X labels			*
C*	NXLBL		INTEGER		Number of x labels		*
C*      IXLBFR          INTEGER         Label frequency			*
C*      IXGLFR          INTEGER         Grid line frequency		*
C*      IXTMFR          INTEGER         Tick mark frequency             *
C*	IRET		INTEGER		Return code			*
C*					  0 = normal return		*
C*                                       -5 = need range for axis	*
C*					 -6 = x coords not identical	*
C**									*
C* Log:									*
C* M. desJardins/GSFC	11/88	GEMPAK 4.1				*
C* M. desJardins/GSFC	 4/89	Check for XAXIS separated with ;	*
C* K. Brill/GSC          5/90   Changes for IN_AXIS			*
C* S. Schotz/GSC	 7/90	Update for new IN_AXIS calling sequence	*
C* S. Schotz/GSC	 8/90	Added call to GEPLOT for error write	*
C* K. Brill/NMC		 1/91   Set default grid line freq to zero	*
C************************************************************************
	INCLUDE		'GEMPRM.PRM'
C*
	CHARACTER*(*)	pname1, pname2, xaxis
	REAL		xlbl (*)
C*
	LOGICAL         skew
C----------------------------------------------------------------------
	iret = 0
	ivc = 0
	skew = .false.
	IF  ( iyaxis .eq. 4 )  skew = .true.
C
C*	Check for unit consistency if two temperatures are requested.
C*	Type  0 = none,  1 = degrees C,  2 = degrees F,  3 = K.
C
	IF  ( pname1 .eq. ' ' )  THEN
	    it1 = 0
	  ELSE IF  ( ( pname1 .eq. 'TEMP' ) .or. ( pname1 .eq. 'TMPC' )
     +		     .or. ( pname1 .eq. 'DWPT' )
     +		     .or. ( pname1 .eq. 'DWPC' ) )  THEN
	    it1 = 1
	  ELSE IF ( ( pname1 .eq. 'TMPF' ) .or. ( pname1 .eq. 'DWPF' ) )
     +						    THEN
	    it1 = 2
	  ELSE IF ( ( pname1 .eq. 'TMPK' ) .or. ( pname1 .eq. 'DWPK' ) )
     +						    THEN
	    it1 = 3
	  ELSE IF ( ( pname1 .eq. 'THTA' ) .or. ( pname1 .eq. 'THTE' ) )
     +						    THEN
	    it1 = 4
	  ELSE
	    it1 = 0
	END IF
C*
	IF  ( pname2 .eq. ' ' )  THEN
	    it2 = 0
	  ELSE IF  ( ( pname2 .eq. 'TEMP' ) .or. ( pname2 .eq. 'TMPC' )
     +		     .or. ( pname2 .eq. 'DWPT' ) .or.
     +		     ( pname2 .eq. 'DWPC' ) )  THEN
	    it2 = 1
	  ELSE IF ( ( pname2 .eq. 'TMPF' ) .or. ( pname2 .eq. 'DWPF' ) )
     +						    THEN
	    it2 = 2
	  ELSE IF ( ( pname2 .eq. 'TMPK' ) .or. ( pname2 .eq. 'DWPK' ) )
     +						    THEN
	    it2 = 3
	  ELSE IF ( ( pname2 .eq. 'THTA' ) .or. ( pname2 .eq. 'THTE' ) )
     +						    THEN
	    it2 = 4
	  ELSE
	    it2 = 0
	END IF
C
C*	Check for default parameters.
C
	IF  ( ( it1 .eq. 0 ) .and. ( it2 .ne. 0 ) )  it1 = it2
	IF  ( ( it2 .eq. 0 ) .and. ( it1 .ne. 0 ) )  it2 = it1
C
C*	Check that the types are the same.
C
	IF  ( it1 .ne. it2 )  THEN
	    iret = -6
            CALL GEPLOT ( ier )
	    CALL ER_WMSG  ( 'SNPROF', iret, ' ', ier )
	    RETURN
	END IF
C
C*	Interpret user input and setup the axis.
C
	ilfdef = 1
	igfdef = 0
	itfdef = 1
	CALL IN_AXIS ( xaxis, ivc, skew, pname1, xmin, xmax,
     +                 ilfdef, igfdef, itfdef, xstrt, xstop, xlbl,
     +                 nxlbl, ixlbfr, ixglfr, ixtmfr, iret )
C*
	IF ( iret .ne. 0 ) THEN
	  iret = -5
	  CALL ER_WMSG ( 'SNPROF', iret, ' ', ier )
	END IF
	RETURN
	END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


text:
On 2nd Street in front of the Train Station train station there is a tall piece of Town Art art titled Solar Intersections.  It was created in 1989 by Robert Behrens.

The statue is covered in holographic material that shimmers with rainbows when you look at it, particularly under the bright summer sun.  It is similar to the spinning discs found in nature stores in the mall.

The holographic Mylar was replaced a couple of years ago, but it degrades fairly quickly in sunlight. Many outdoor sculptures in the 1980s used this material, before the product lifespan was known.

20101122 01:13:38 nbsp Its gorgeous and eyecatching when the sunlight hits it just right, but this is just for a few minutes on the right morning. The rest of the time theyre just a bunch of halfhidden white poles nobody notices. Users/SteveDavison


text:
        SUBROUTINE  MA_LSCD  ( mszrpt, marrpt, ipt, iret )
C************************************************************************
C* MA_LSCD					 	                *
C*							                *
C* This subroutine gets the length and starting position of each of the *
C* sections 1-4 in the ZZYY drifting buoy report.  If a section is      *
C* missing, then its length and starting position are set to zero.      *
C*								        *
C* MA_LSCD  ( MSZRPT, MARRPT, IPT, IRET )                               *
C*								        *
C* Input parameters:						        *
C*      MSZRPT          INTEGER         Length of report in bytes       *
C*      MARRPT          CHAR*           Report array                    *
C*      IPT             INTEGER         Pointer to Q(t) in 6Q(i)Qt)//   *
C*								        *
C* Output parameters:						        *
C*	LSEC1           INTEGER         Length of section 1 in report   *
C*      ISEC1           INTEGER         Pointer to start of section 1   *
C*	LSEC2           INTEGER         Length of section 2 in report   *
C*      ISEC2           INTEGER         Pointer to start of section 2   *
C*	LSEC3           INTEGER         Length of section 3 in report   *
C*      ISEC3           INTEGER         Pointer to start of section 3   *
C*	LSEC4           INTEGER         Length of section 4 in report   *
C*      ISEC4           INTEGER         Pointer to start of section 4   *
C*	IRET		INTEGER         Return code                     *
C*				          0 = normal return             *
C*                                        1 = problems                  *
C**								        *
C* Log:								        *
C* R. Hollern/NCEP       6/96                                           *
C* R. Hollern/NCEP       8/96 	Added check on section lengths          *
C* R. Hollern/NCEP      11/96  	Modified logic to compute length of     *
C*                             	sections                                *
C* D. Kidwell/NCEP       4/97	Cleaned up code and documentation       *
C************************************************************************
        INCLUDE 	'macmn.cmn'
C*
        CHARACTER*(*)   marrpt
C------------------------------------------------------------------------
        iret  = 0
        lsec1 = 0
        lsec2 = 0
        lsec3 = 0
        lsec4 = 0
        isec1 = 0
        isec2 = 0
        isec3 = 0
        isec4 = 0
C
C*      If there is a section 444, determine its length.
C
        j1   = ipt
        lrpt = mszrpt
C
        i444 = INDEX ( marrpt (j1:lrpt), ' 444 ' )
C
        IF ( i444 .gt. 0 ) THEN
            ka    = j1 + i444
            isec4 = ka + 2
            lsec4 = lrpt - isec4 + 1
            lrpt  = ka - 1
        END IF
C
C*      If there is a section 333, determine its length.
C
        i333 = INDEX ( marrpt (j1:lrpt), ' 333' )
C
        IF ( i333 .gt. 0 ) THEN
            ka    = j1 + i333
            isec3 = ka + 3
            lsec3 = lrpt - isec3 + 1
            lrpt  = ka - 1
        END IF
C
C*      If there is a section 222, determine its length.
C
        i222 = INDEX ( marrpt (j1:lrpt), ' 222' )
C
        IF ( i222 .gt. 0 ) THEN
            ka    = j1 + i222
            isec2 = ka + 3
            lsec2 = lrpt - isec2 + 1
            lrpt  = ka - 1
        END IF
C
C*      If there is a section 111, determine its length.
C
        i111 = INDEX ( marrpt (j1:lrpt), ' 111' )
C
        IF ( i111 .gt. 0 ) THEN
            ka    = j1 + i111
            isec1 = ka + 3
            lsec1 = lrpt - isec1 + 1
            lrpt  = ka - 1
        END IF
C
C*      If section lengths are too long, reject report.
C*      This can happen if reports are not separated from
C*      each other by a report separator.
C
        IF ( lsec1 .gt. 50 .or. lsec2 .gt. 35 ) iret = 1
C*
        RETURN
        END