       program insola
*******************************************************************************
*
*              Computation of the various insolation functions
*
*              this program uses
*                   insolsub.f   : subroutines for the computation
*                                  of insolation
*                   insola.par   : parameters for insola.f
*                   nomfich      : binary file, created by
*                                  prepinsol
*                                  and which contains the necessary
*                                  orbital and precession quantities
*
*
*
*   the following parameters are taken from the NAMELIST insola.par
*
*           nomfich     : binary file for the elements t,k,h,eps,psi          *
*           pas         : stepsize (years)                                    *
*           datefin     : final time (million of years)                       *
*           datedebut   : starting time  (million of years)                   *
*           statut      : status for the opening of results file              *
*                         'new' ou 'unknown'                                  *
*           so          : solar constant   W*m^-2                             *
*                                                                             *
*   other parameters:                                                         *
*                                                                             *
*           nbf         : 4 (number of variables)                             *
*           nbel1       : 4 + 1                                               *
*           nacd        : numbers of n-uplets in a record                     *
*                                                                             *
*  (c) Astronomie et Systemes Dynamiques/ Bureau des Longitudes (1993)        *
*  J. Laskar, F. Joutel                                                       *
*  version 0.82 (7/7/1993)                                                    *
*  version 0.83 (27/1/93)                                                     *
*******************************************************************************

      implicit double precision(a-h,o-z)
      parameter(pi=3.1415926535897932d0)
      parameter(nacd=100,nbel=5,nbf=4)
      dimension Tel(nbf)
      character *50 nomfich,fichout
      character *20 statut
      common/cdo/so
      common/date/pas,itdebut
c
      NAMELIST/NAMSTD/nomfich,pas,datedebut,datefin,statut,so
      open(18,file='insola.par',form='formatted',status='old')
      read(18,NAMSTD)
c      write (*,NAMSTD)
      close(18)
c
c Les fichiers ont une longueur d'enregistrement correspondant
c a 100 000 ans.
      nrecl=8*nbel*nacd
      open(10,file=nomfich,status='old',access='direct',recl=nrecl)
c
      read(10,rec=1) debut
      itdebut=int(debut)

c Nombre de lignes lues:
      datedebut=datedebut*1.D6
      datefin=datefin*1.D6
      aux=dabs(datefin-datedebut)
      npt=int((aux/dabs(pas))+1)
*-----------------------------------------------------------------------
*   principal menu
*-----------------------------------------------------------------------
      write (*,*)
      write (*,*)
      write (*,*) '***************************************************'
      write (*,*) '*                                                 *'
      write (*,*) '*                                                 *'
      write (*,*) '*                                                 *'
      write (*,*) '*                  INSOLA                         *'
      write (*,*) '*                                                 *'
      write (*,*) '*   computation of various insolation quantities  *'
      write (*,*) '*                                                 *'
      write (*,*) '*              (c) ASD/BdL (1993)                 *'
      write (*,*) '*                                                 *'
      write (*,*) '*                                                 *'
      write (*,*) '*                                                 *'
      write (*,*) '*                                                 *'
      write (*,*) '***************************************************'
999   continue
***** V 0.83 suppression des lignes suivantes
***      insojour =.FALSE.
***      insojcal =.FALSE.
***      insomcal =.FALSE.
***      insoann  =.FALSE.
      write (*,*)
      write (*,*)
      write (*,*) '    1 : mean daily    insolation/true longitude '
      write (*,*) '    2 : mean daily    insolation/mean longitude '
      write (*,*) '    3 : mean monthly  insolation '
      write (*,*) '    4 : mean annual   insolation '
      write (*,*) '    5 : HELP'
      write (*,*) '    6 : QUIT'
      write (*,*)
      write (*,*) ' your choice ? '
      read  (*,*) ICHOI
      GOTO (1000,2000,3000,4000,5000,7000), ICHOI
      write (*,*) ' incorrect choice '
      GOTO 999

*---------------------------------- mean daily/ true longitude
1000  CONTINUE
      write (*,*)' Latitude on the Earth (in degrees) ? '
      read (*,*)  phi
      write (*,*) ' latitude = ',phi
      write (*,*)' True longitude (in degrees) ?'
      read (*,*) datejour
      write (*,*) 'true longitude = ',datejour
      phi=phi*pi/180.d0
      GOTO 6000


*---------------------------------- mean daily/ mean longitude
2000  CONTINUE
      write (*,*)' Latitude on the Earth (in degrees) ? '
      read (*,*)  phi
      write (*,*) ' latitude = ',phi
      write (*,*)
     $'approximate conventional dates for a given mean longitude'
      write (*,*)
      write (*,*)'  0: 21 march    '
      write (*,*)' 30: 21 april    '
      write (*,*)' 60: 21 may      '
      write (*,*)' 90: 21 june     '
      write (*,*)'120: 21 july     '
      write (*,*)'150: 21 august   '
      write (*,*)'180: 21 september'
      write (*,*)'210: 21 october  '
      write (*,*)'240: 21 november '
      write (*,*)'270: 21 december  '
      write (*,*)'300: 21 january   '
      write (*,*)'330: 21 february  '
      write (*,*)
      write (*,*)' MEAN longitude (in degrees)  (0-360)?'
      read (*,*) datecal
      write (*,*) 'MEANlongitude = ',datecal
      phi=phi*pi/180.d0
      GOTO 6000

*---------------------------------- mean monthly
3000  CONTINUE
      write (*,*)' Latitude on the Earth (in degrees) ? '
      read (*,*)  phi
      write (*,*) ' latitude = ',phi
      write (*,*)
      write (*,*)'approximate conventional dates of the months'
      write (*,*)
      write (*,*)' 1: 21 december  - 20 january'
      write (*,*)' 2: 21 january   - 20 february'
      write (*,*)' 3: 21 february  - 20 march'
      write (*,*)' 4: 21 march     - 20 april'
      write (*,*)' 5: 21 april     - 20 may'
      write (*,*)' 6: 21 may       - 20 june'
      write (*,*)' 7: 21 june      - 20 july'
      write (*,*)' 8: 21 july      - 20 august'
      write (*,*)' 9: 21 august    - 20 september'
      write (*,*)'10: 21 september - 20 october'
      write (*,*)'11: 21 october   - 20 november'
      write (*,*)'12: 21 november  - 20 december'
      write (*,*)
      write (*,*) 'Your choice ? (1-12)'
      read (*,*) moiscal
      write (*,*) 'month = ',moiscal
      phi=phi*pi/180.d0
      GOTO 6000
*---------------------------------- mean annual
4000  CONTINUE
      GOTO 6000


*--------------------------- Computation of insolation ------------
6000  continue
      write (*,*)
      write (*,*) 'Name of the output file'
      read(*,1010) fichout
1010  format (A)
      open (22,file = fichout,form='formatted')
      tps=datedebut
      do i=1,npt
         call Telinsol(tps,Tel)
         ak=Tel(1)
         ah=Tel(2)
         eps=Tel(3)
         psi=Tel(4)
         goto (6100,6200,6300,6400) ICHOI

6100     continue
         call wjour(datejour,ak,ah,eps,psi,phi,w)
         goto 6666

6200     continue
         call wjcal(datecal,ak,ah,eps,psi,phi,w)
         goto 6666

6300     continue
         call wmcal(moiscal,ak,ah,eps,psi,phi,w)
         goto 6666

6400     continue
         call wam(ak,ah,w)
         goto 6666

6666     continue
         write(22,1020) tps*1D-3,w 
1020  format(1x,f10.0,d20.12)
         tps=tps+pas
      end do
      close (22)
      GOTO 999
*--------------------------- HELP ------------------------------------
5000  write (*,*)
      write (*,*)' 1,2 :'
      write (*,*)'The daily insolation is computed for a given latitude'
      write (*,*)'of the Earth, and for a given position of the Earth '
      write (*,*)'on its orbit. This position is sometime given by a '
      write (*,*)'conventional date where the origine is at the moving'
      write (*,*)'equinox with conventional date 21 march'
      write (*,*)'As this terminology is confusing, in the present '
      write (*,*)'program, the position of the Earth on its orbit is'
      write (*,*)'given either by its true longitude (1), or by the'
      write (*,*)'mean longitude (2), both expressed in degrees'
      write (*,*)'with origine at the moving equinox'
      write (*,*)
      write (*,*)'The true longitude is the true position angle of the'
      write (*,*)'Sun-Earth direction with respect to the equinox. '
      write (*,*)'The mean longitude is a fictif angle which is '
      write (*,*)'described by the Earth at constant velocity.'
      write (*,*)'The mean longitude is thus proportional to the time,'
      write (*,*)'and  should be the normal choice in many cases.'
      write (*,*)
      write (*,*)
      write (*,*) ' press RETURN for MORE '
      PAUSE
      write (*,*)' 3 :'
      write (*,*)' The monthly insolation correspond to the mean daily'
      write (*,*)'insolation over 1/12 of a year, that is 30 degrees'
      write (*,*)'of mean longitude, with origine at the equinox.'
      write (*,*)
      write (*,*)' 4 :'
      write (*,*)'The mean annual insolation is computed over a full'
      write (*,*)'revolution of the Earth around the Sun.'
      write (*,*)
      write (*,*)
      GOTO 999
7000  CONTINUE
      end

