SUBROUTINE solard(r1)
  !
  ! $Author: CPTEC $
  ! $Date: 2000/02/24 13:15:43 $
  ! $Revision: 1.0 $
  !
  !******************************************************************
  ! subprogram:    solard      compute the solar-earth distance
  !   prgrmmr: q.zhao           org: w/nmc2     date: 96-7-23
  ! abstract:
  !   solard calculates the solar-earth distance on each day
  !   for use in short-wave radiation.
  ! program history log:
  !   96-07-23  q.zhao      - originator
  ! usage: call solard from subroutine radtn
  !   output argument list:
  !    r1   - the non-dimensional distance between sun and the earth
  !          (less than 1.0 in summer and larger than 1.0 in winter).
  IMPLICIT NONE
  REAL, PARAMETER ::pi=3.1415926,pi2=2.*pi

  INCLUDE "ctlblk.f90"

  REAL, INTENT(inout) :: r1
  INTEGER, DIMENSION(12) :: ndm
  INTEGER                 :: jyr19,jmn,jdor1,jdor2,jd,jyr,jmnth,jday,jhr,jyrm9
  INTEGER                 :: lp, ic, ny, iter
  REAL                    :: ccr,tpp,fjd,dat,t,ec,em,date,cr,ep,e, year
  DATA jyr19/1900/,jmn/0/,ccr/1.3e-6/
  DATA ndm/0,31,59,90,120,151,181,212,243,273,304,334/

  !tpp = days between epoch and perihelion passage of 1900
  !jdor1 = jd of december 30, 1899 at 12 hours ut
  !jdor2 = jd of epoch which is january 0, 1990 at 12 hours ut
  DATA tpp/1.55/,jdor1/2415019/,jdor2/2415020/

  !computes julian day and fraction from year, month, day and time ut
  !accurate only between march 1, 1900 and february 28, 2100
  !based on julian calendar corrected to correspond to gregorian
  !calendar during this period
  jd=jdor1
  jyr=idat(3)
  jmnth=idat(1)
  jday=idat(2)
  jhr=ihrst
  jyrm9=jyr-jyr19
  lp=jyrm9/4
  IF(lp > 0) jd=jd+1461*lp
  ny=jyrm9-4*lp
  ic=0
  IF(ny > 0) THEN 
     jd=jd+365*ny+1
  ELSE
     IF(jmnth > 2) ic=1
     jd=jd+ndm(jmnth)+jday+ic
  END IF
  IF(jhr >= 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

  !calculate the solar-earth distance
  dat=float(jd-jdor2)-tpp+fjd

  !computes time in julian centuries after epoch
  t=float(jd-jdor2)/36525.e0

  !computes length of anomalistic and tropical years (minus 365 days)
  year=.25964134e0+.304e-5*t

  !computes orbit eccentricity  from t
  ec=.01675104e0-(.418e-4+.126e-6*t)*t
  year=year+365.e0

  !date=days since last perihelion passage
  date = MOD(dat,year)

  !solve orbit equations by newton's method
  em=pi2*date/year
  e=1.e0
  iter = 0
  DO
     ep=e-(e-ec*SIN(e)-em)/(1.e0-ec*COS(e))
     cr=ABS(e-ep)
     e=ep
     iter = iter + 1
     IF(iter > 10) EXIT
     IF(cr > ccr) CYCLE
  END DO
  r1=1.e0-ec*COS(e)

  !write(6,1000)jyr,jmnth,jday,jhr,r1
  !1000 format('sun-earth distance calculation finished in solard'/  &
  !    'year=',i5,'  month=',i3,'  day=',i3,' hour=' ,      i3,' r1=',f9.4)
END SUBROUTINE solard
