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