SUBROUTINE solarlong(pday,psollong,pdist_sol) USE ioipsl IMPLICIT NONE c======================================================================= c c Objet: c ------ c c Calcul de la distance soleil-planete et de la declinaison c en fonction du jour de l'annee. c c c Methode: c -------- c c Calcul complet de l'elipse c c Interface: c ---------- c c Uncommon comprenant les parametres orbitaux. c c Arguments: c ---------- c c Input: c ------ c pday jour de l'annee (le jour 0 correspondant a l'equinoxe) c lwrite clef logique pour sorties de controle c c Output: c ------- c pdist_sol distance entre le soleil et la planete c ( en unite astronomique pour utiliser la constante c solaire terrestre 1370 Wm-2 ) c pdecli declinaison ( en radians ) c c======================================================================= c----------------------------------------------------------------------- c Declarations: c ------------- #include "planete.h" #include "YOMCST.h" include 'iniprint.h' c arguments: c ---------- REAL pday,pdist_sol,pdecli,psollong LOGICAL lwrite c Local: c ------ REAL zanom,xref,zx0,zdx,zteta,zz,pi INTEGER iter REAL :: pyear_day,pperi_day REAL :: jD_eq, jD_peri LOGICAL,SAVE :: first=.TRUE. c$OMP THREADPRIVATE(first) c----------------------------------------------------------------------- c calcul de l'angle polaire et de la distance au soleil : c ------------------------------------------------------- c Initialisation eventuelle: if(first) then call ioget_calendar(pyear_day) call ymds2ju(2000, 3, 21, 0., jD_eq) call ymds2ju(2001, 1, 4, 0., jD_peri) pperi_day = jD_peri - jD_eq pperi_day = R_peri + 180. write(lunout,*)' Number of days in a year = ',pyear_day c call iniorbit(249.22,206.66,669.,485.,25.2) call iniorbit(152.59,146.61,pyear_day,pperi_day,R_incl) first=.FALSE. endif c calcul de l'zanomalie moyenne zz=(pday-peri_day)/year_day pi=2.*asin(1.) zanom=2.*pi*(zz-nint(zz)) xref=abs(zanom) c resolution de l'equation horaire zx0 - e * sin (zx0) = xref c methode de Newton ! zx0=xref+e_elips*sin(xref) zx0=xref+R_ecc*sin(xref) DO 110 iter=1,10 ! zdx=-(zx0-e_elips*sin(zx0)-xref)/(1.-e_elips*cos(zx0)) zdx=-(zx0-R_ecc*sin(zx0)-xref)/(1.-R_ecc*cos(zx0)) if(abs(zdx).le.(1.e-7)) goto 120 zx0=zx0+zdx 110 continue 120 continue zx0=zx0+zdx if(zanom.lt.0.) zx0=-zx0 c zteta est la longitude solaire ! zteta=2.*atan(sqrt((1.+e_elips)/(1.-e_elips))*tan(zx0/2.)) zteta=2.*atan(sqrt((1.+R_ecc)/(1.-R_ecc))*tan(zx0/2.)) psollong=zteta-timeperi IF(psollong.LT.0.) psollong=psollong+2.*pi IF(psollong.GT.2.*pi) psollong=psollong-2.*pi psollong = psollong * 180. / pi c distance soleil pdist_sol = (1-R_ecc*R_ecc) & /(1+R_ecc*COS(pi/180.*(psollong-(R_peri+180.0)))) ! pdist_sol = (1-e_elips*e_elips) ! & /(1+e_elips*COS(pi/180.*(psollong-(R_peri+180.0)))) c----------------------------------------------------------------------- c sorties eventuelles: c --------------------- c IF (lwrite) THEN c PRINT*,'jour de l"annee :',pday c PRINT*,'distance au soleil (en unite astronomique) :',pdist_sol c PRINT*,'declinaison (en degres) :',pdecli*180./pi c ENDIF RETURN END