C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
                         SUBROUTINE SOLARD(R1)
C     ******************************************************************
C$$$  SUBPROGRAM DOCUMENTATION BLOCK
C                .      .    .                               .
C SUBPROGRAM:    SOLARD      COMPUTE THE SOLAR-EARTH DISTANCE
C   PRGRMMR: Q.ZHAO           ORG: W/NMC2     DATE: 96-7-23       
C     
C ABSTRACT:
C     SOLARD CALCULATES THE SOLAR-EARTH DISTANCE ON EACH DAY
C     FOR USE IN SHORT-WAVE RADIATION.
C     
C PROGRAM HISTORY LOG:
C   96-07-23  Q.ZHAO      - ORIGINATOR
C   98-10-09  Q.ZHAO      - CHANGED TO USE IW3JDN IN W3LIB TO
C                           CALCULATE JD.
C     
C USAGE: CALL SOLARD FROM SUBROUTINE INIT
C
C   INPUT ARGUMENT LIST:
C       NONE
C  
C   OUTPUT ARGUMENT LIST: 
C       R1   - THE NON-DIMENSIONAL DISTANCE BETWEEN SUN AND THE EARTH
C              (LESS THAN 1.0 IN SUMMER AND LARGER THAN 1.0 IN WINTER).
C     
C   INPUT FILES:
C     NONE
C        
C   OUTPUT FILES:
C     NONE
C     
C   SUBPROGRAMS CALLED:
C  
C     UNIQUE: NONE
C  
C     LIBRARY: IW3JDN
C  
C   COMMON BLOCKS: CTLBLK
C   
C ATTRIBUTES:
C   LANGUAGE: FORTRAN 90
C   MACHINE : IBM SP
C$$$  
#include "sp.h"
C***********************************************************************
                             P A R A M E T E R
     & (PI=3.1415926,PI2=2.*PI)
C-----------------------------------------------------------------------
                             L O G I C A L
     & RUN,FIRST,RESTRT,SIGMA,LEAP
C-----------------------------------------------------------------------
      INCLUDE "CTLBLK.comm"
C-----------------------------------------------------------------------
                             D I M E N S I O N
     1 NDM (12)
C-----------------------------------------------------------------------
                             D A T A
     1   JYR19/1900/,     JMN/0/,  CCR/1.3E-6/
C
                             D A T A
     1   NDM/0,31,59,90,120,151,181,212,243,273,304,334/
C
C.....TPP = DAYS BETWEEN EPOCH AND PERIHELION PASSAGE OF 1900
C.....JDOR1 = JD OF DECEMBER 30, 1899 AT 12 HOURS UT
C.....JDOR2 = JD OF EPOCH WHICH IS JANUARY 0, 1990 AT 12 HOURS UT
C
                                   D A T A
     1   TPP/1.55/,    JDOR1/2415019/,     JDOR2/2415020/
C
C    *******************************************************************
C     COMPUTES JULIAN DAY AND FRACTION FROM YEAR, MONTH, DAY AND TIME UT
C     ACCURATE ONLY BETWEEN MARCH 1, 1900 AND FEBRUARY 28, 2100
C     BASED ON JULIAN CALENDAR CORRECTED TO CORRESPOND TO GREGORIAN
C        CALENDAR DURING THIS PERIOD
C    *******************************************************************
C
      INTEGER iyr,imo,idy,cyr,cmo,cdy,utc,cutc

C------------------------------Arguments--------------------------------
C Input arguments
C
      real calday,     ! Calendar day, including fraction
     $     eccen,      ! Eccentricity
     $     obliqr,     ! Earth's obliquity in radians
     $     lambm0,     ! Mean longitude of perihelion at the 
C                      ! vernal equinox (radians)
     $     mvelpp      ! Earth's moving vernal equinox longitude
C                      ! of perihelion plus pi (radians)

C
C Output arguments
C
      real delta,   ! Solar declination angle in radians
     $     eccf     ! Earth-sun distance factor ( i.e. (1/r)**2 )
C
C---------------------------Local variables-----------------------------
C
      real
     $     pie,                ! pi
     $     dayspy,              ! days per year
     $     DOY
      data dayspy / 365. /
      real ve          ! Calday of vernal equinox
      parameter( ve = 80.5 ) ! correct for Jan 1 = calday 1
      parameter( pie = 3.14159265358979323846 )
C
      real   lambm,   ! Lambda m, earth's mean longitude of perihelion (radians)
     $       lmm,     ! Intermediate argument involving lambm
     $       lamb,    ! Lambda, the earth's longitude of perihelion
     $       invrho,  ! Inverse normalized sun/earth distance
     $       sinl     ! Sine of lmm
C
C Compute eccentricity factor and solar declination using
C day value where a round day (such as 213.0) refers to 0z at
C Greenwich longitude.
C
C Use formulas from Berger, Andre 1978: Long-Term Variations of Daily
C Insolation and Quaternary Climatic Changes. J. of the Atmo. Sci.
C 35:2362-2367.
C
C To get the earth's true longitude (position in orbit; lambda in Berger 1978),
C which is necessary to find the eccentricity factor and declination,
C must first calculate the mean longitude (lambda m in Berger 1978) at
C the present day.  This is done by adding to lambm0 (the mean longitude
C at the vernal equinox, set as March 21 at noon, when lambda = 0; in radians)
C an increment (delta lambda m in Berger 1978) that is the number of
C days past or before (a negative increment) the vernal equinox divided by
C the days in a model year times the 2*pi radians in a complete orbit.
C       
  
C      do i=1,365

CJFP
C     inicio: consegue a data completa atual na integracao

      iyr=IDAT(3)
      imo=IDAT(1)
      idy=IDAT(2)
      
C      IF(mod(NTSD*96/3600,6).EQ.0) THEN
       utc=(NTSD*DT)/3600
       CALL getdate(iyr,imo,idy,utc,cyr,cmo,cdy,cutc)
C       write(*,*) 'data inicial:',iyr,imo,idy,utc,NTSD*DT/3600
C       write(*,*) 'data atual:',cyr,cmo,cdy,cutc
       JYR=cyr
       JMNTH=cmo
       JDAY=cdy
       JHR=cutc
C      final:ja possui a data completa
CJFP
C
      JD=IW3JDN(JYR,JMNTH,JDAY)
      IF(JHR.GE.12) THEN
      JD=JD-1
      FJD=.5E0+.041666667E0*FLOAT(JHR)+.00069444444E0*FLOAT(JMN)
      ELSE
  7   FJD=.041666667E0*FLOAT(JHR-12)+.00069444444E0*FLOAT(JMN)
      END IF
      DAYINC=JHR/24.0
      JD=JD+FJD+DAYINC
      FJD=JD+FJD+DAYINC-JD
C***
C*** CALCULATE THE SOLAR-EARTH DISTANCE
C***
      DAT=FLOAT(JD-JDOR2)-TPP+FJD
C***
C    COMPUTES TIME IN JULIAN CENTURIES AFTER EPOCH
C***
      T=FLOAT(JD-JDOR2)/36525.E0
C***
C    COMPUTES LENGTH OF ANOMALISTIC AND TROPICAL YEARS (MINUS 365 DAYS)
C***
      YEAR=.25964134E0+.304E-5*T
C***
C    COMPUTES ORBIT ECCENTRICITY  FROM T
C***
      EC=.01675104E0-(.418E-4+.126E-6*T)*T
      YEAR=YEAR+365.E0
C***
C    DATE=DAYS SINCE LAST PERIHELION PASSAGE
C***
      DATE = MOD(DAT,YEAR)
C***
C    SOLVE ORBIT EQUATIONS BY NEWTON'S METHOD     '
C***
      EM=PI2*DATE/YEAR
      E=1.E0
      ITER = 0
 31   EP=E-(E-EC*SIN(E)-EM)/(1.E0-EC*COS(E))
      CR=ABS(E-EP)
      E=EP
      ITER = ITER + 1
      IF(ITER.GT.10) GOTO 1031
      IF(CR.GT.CCR) GO TO 31
 1031 CONTINUE
      R1=1.E0-EC*COS(E)
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      DOY = ((275*cmo)/9) - 2*((cmo+9)/12) + cdy - 30
      calday = DOY               ! Calendar day, including fraction
      eccen  =1.8681820E-02      ! Eccentricity
      obliqr =0.4207183          ! Earth's obliquity in radians
      lambm0 =3.156770           ! Mean longitude of perihelion at the vernal equinox (radians)
      mvelpp =-5.7511014E-04     ! Earth's moving vernal equinox longitude of perihelion plus pi (radians)

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

      lambm  = lambm0 + (calday - ve)*2.*pie/dayspy
      lmm    = lambm  - mvelpp
C
C The earth's true longitude, in radians, is then found from
C the formula in Berger 1978:
C
      sinl   = sin(lmm)
      lamb   = lambm  + eccen*(2.*sinl 
     $         + eccen*(1.25*sin(2.*lmm)
     $         + eccen*((13.0/12.0)*sin(3.*lmm) - 0.25*sinl)))
C
C Using the obliquity, eccentricity, moving vernal equinox longitude of
C perihelion (plus), and earth's true longitude, the declination (delta)
C and the normalized earth/sun distance (rho in Berger 1978; actually inverse
C rho will be used), and thus the eccentricity factor (eccf), can be calculated
C from formulas given in Berger 1978.
C
      invrho = (1. + eccen*cos(lamb - mvelpp))
     $         / (1. - eccen*eccen)
C
C Set solar declination and eccentricity factor
C
      delta  = asin(sin(obliqr)*sin(lamb))
      eccf   = invrho*invrho
C
      print *, 'R1 and DAY'
      print *, eccf, DOY

      R1=eccf
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC



C
c      WRITE(6,1000)JYR,JMNTH,JDAY,JHR,R1
c 1000 FORMAT('SUN-EARTH DISTANCE CALCULATION FINISHED IN SOLARD'/
c     &       'YEAR=',I5,'  MONTH=',I3,'  DAY=',I3,' HOUR='
c    &,      I3,' R1=',F9.4)
C***
C    RETURN TO RADTN
C***
      RETURN
      END
