!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! ! THIS IS DragoRT1.0, THE IMPLEMENTATION OF THE PBMC ALGORITHM DESCRIBED IN: ! ! PRE-CONDITIONED BACKWARD MONTE CARLO SOLUTIONS TO RADIATIVE ! TRANSPORT IN PLANETARY ATMOSPHERES. FUNDAMENTALS: ! SAMPLING OF PROPAGATION DIRECTIONS IN POLARISING MEDIA ! BY A. GARCIA MUÑOZ & F.P. MILLS ! ASTRONOMY & ASTROPHYSICS, XXXXXXXX, 2014 ! ! THE MANUSCRIPT IS ALSO AVAILABLE AT ASTRO.PH, ID: XXXXXXXXX ! ! DragoRT1.0 M CALCULATES THE STOKES VECTOR FOR OUTGOING RADIATION ! REFLECTED FROM A PLANE-PARALLEL ATMOSPHERE ! ! THE REFERENCE INCLUDES NUMEROUS TEST CASES THAT EXPLORE THE ! ALGORITHM'S PERFORMANCE FOR BOTH RAYLEIGH AND MIE SCATTERING MEDIA ! ! LAST MODIFIED ON 01/AUGUST/2014 ! ! FOR A SKETCH OF THE RELEVANT GEOMETRICAL ANGLES, SEE FIGURE 5 IN THE REFERENCE. ! IN BRIEF: ! COSPOLOBS & COSPOLSUN ARE THE COSINES OF THE OBSERVER AND SOLAR POLAR ANGLES ! AZIMUTH IS THE AZIMUTH ANGLE BETWEEN THE OBSERVER AND SOLAR PLANES. ! ! THIS PACKAGE IS ACCOMPANIED BY: ! * INPUT FILES: INPUT_hazeL.txt, INPUT_L13.txt, INPUT_L60.txt; ! THEY CONTAIN EXPLANATIONS TO THE INPUT PARAMETERS. ! COPY INPUT_XXXX.txt INTO INPUT.dat TO EXECUTE SOME OF THE EXAMPLES DESCRIBED ! IN THE REFERENCE ! * FILES WITH SCATTERING MATRIX PROPERTIES: phF_hazeL.txt, phF_L13.txt, phF_L60.txt ! * SCRIPT FOR COMPILATION IN GFORTRAN (myscript) ! ! FOR RAYLEIGH ATMOSPHERES, ONE CAN SET UP A SIMILAR phF FILE. OTHERWISE, ONE CAN ! INCORPORATE THE RELEVANT SCATTERING MATRIX BY UNCOMMENTING THE 9 ROWS FOLLOWING: ! "FOR RAYLEIGH SCATTERING ATMOSPHERES, UNCOMMENT THE FOLLOWING 9 ROWS" ! IN THE FORTRAN CODE ! ! THE PROGRAM WILL OUTPUT THE FILE OUT.DAT WITH A SUMMARY OF THE INPUT PARAMETERS AND ! THE OUTPUT STOKES VECTORS. NOTE THAT THE DEFINITION HERE OF POSITIVE V STOKES ELEMENT ! DIFFERS FROM THE DEFINITION IN THE SOURCES FOR THE TEST CASES. ! ! THE PROGRAM IS FREELY AVAILABLE FOR ANY NON-FOR-PROFIT APPLICATION. ! IF YOU USE IT, PLEASE GIVE CREDIT TO THE ABOVE REFERENCE. ! ! FOR COMMENTS, QUESTIONS, ETC, PLEASE CONTACT A. GARCIA MUÑOZ ON: TONHINGM@GMAIL.COM ! ALSO, IF YOU SPOT A BUG OR INCONSISTENCY, PLEASE DROP ME AN EMAIL. ! THE INTENT IS TO UPDATE DragoRT TO IMPROVE ITS PERFORMANCE AND INCLUDE NEW FEATURES, ! SO FEEL WELCOME TO DROP ME AN EMAIL TO FIND OUT ABOUT THE LATEST VERSION. ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! PROGRAM PBMC implicit none COMMON /PARAMETERS1/ NPH, NSUN, NTHETAD, NTHETAI, & ex, ey, ez, WEIGHTMIN, TAU, GAMMA, SSA, ALBGR, DXSLAB, & X0, COSPOL0, SS0, COSPOLSUN, AZIMUTH, SSSUN, FLUXSUN, & MM0, THETAD, DTHETAD, THETAI, FFI, DFFI INTEGER, PARAMETER:: NSUNMAX=10, NTHETADMAX= 1800, NTHETAIMAX=4000 REAL, PARAMETER:: PI=4.*atan(1.) REAL, PARAMETER:: MACHPR=1.D-8 INTEGER NPH, NSUN, NTHETAD, NTHETAI REAL ex(1:3), ey(1:3), ez(1:3) REAL WEIGHTMIN, TAU, GAMMA, SSA, ALBGR, DXSLAB REAL X0(1:3), COSPOL0, SS0(1:3) REAL COSPOLSUN(1:NSUNMAX), AZIMUTH(1:NSUNMAX), SSSUN(1:3,1:NSUNMAX), FLUXSUN(1:4) REAL MM0(1:4,1:4,0:NTHETADMAX), THETAD(0:NTHETADMAX), DTHETAD REAL THETAI(0:NTHETAIMAX), FFI(0:NTHETAIMAX), DFFI REAL IVECT0(1:4,1:NSUNMAX), HH0(1:4,1:4), weight0 REAL IVECT0IPH(1:4,1:NSUNMAX) REAL X(1:3), XA(1:3), XB(1:3), SS(1:3), SP(1:3), PP(1:4,1:4) REAL IVECTK(1:4,1:NSUNMAX), HH(1:4,1:4), weight REAL IVECTCHECK INTEGER ISUN, JJ, IPH, IMOVE REAL AXXB, GEOBIAS, LL CHARACTER GOFOR INTEGER IDUM REAL MYRAN, EPSA, RHO CHARACTER(LEN=100) FILENAME !!!!!!!!!!!!! IDUM=0 CALL SEED(IDUM) !!!!!!!!!!!!! !! CALL SETUP !! !!!!!!!!!!!!! FILENAME='./OUT.dat' OPEN(95,FILE=TRIM(FILENAME), ACCESS='APPEND', STATUS='UNKNOWN') WRITE(95,'(A15,20A15)') 'IPH', 'TAU', 'SSA', 'ALBGR', 'COSPOL0BS', 'COSPOLSUN', & 'AZIMUTH', 'STOKES I', 'STOKES Q', 'STOKES U', 'STOKES V' CLOSE(95) !!!!!!!!!!!!! !! INITIALIZATION IVECT0(1:4,1:NSUNMAX)=0. HH0(1:4,1:4)=0. DO JJ=1,4 HH0(JJ,JJ)=1. ENDDO weight0=1. !! !!!!!!!!!!!!! !!!!!!!!!!!!! !! IPH LOOP IPH=0 DO WHILE (IPH .LT. NPH) X=X0 SS=SS0 HH=HH0 weight=weight0 !!!!!!!!!!!!! !! IMOVE LOOP IMOVE=0 IVECT0IPH=0. DO WHILE (weight .GE. weightmin) CALL FINDXB( X, SS, XB, AXXB ) IF (XB(1) .NE. 0. .OR. ALBGR .EQ. 0.) THEN GEOBIAS=AXXB GOFOR = 'A' ELSE GEOBIAS=1. RHO=MYRAN(IDUM) IF (AXXB .GE. RHO .AND. RHO .GT. 0.) GOFOR='A' IF (AXXB .LT. RHO .AND. RHO .LT. 1.) GOFOR='B' ENDIF !!!! IF ( GOFOR .EQ. 'A' ) THEN EPSA=MYRAN(IDUM) LL= -1. / GAMMA * ALOG( 1. + AXXB * (-1. + EPSA) ) XA= X - ss * LL IF (XA(1) .LT. MACHPR) XA(1)= 0. IF (XA(1) .GT. DXSLAB*(1. - MACHPR)) XA(1)= DXSLAB CALL LLMMLLA( IPH, IDUM, HH, SS, SP, PP ) CALL CONTRIBSUNA( XA, SS, GEOBIAS, weight, HH, IVECTk ) !!! HH=MATMUL(HH,PP) HH=HH/HH(1,1) weight=weight * GEOBIAS * SSA X=XA SS=SP ELSEIF ( GOFOR .EQ. 'B' ) THEN LL= -1. / GAMMA * ALOG( 1. - AXXB ) XB= X - ss * LL IF (XB(1) .LT. MACHPR) XB(1)= 0. CALL LLMMLLB( IDUM, SP, PP ) CALL CONTRIBSUNB( XB, weight, HH, IVECTk ) !!! HH=MATMUL(HH,PP) HH=HH/HH(1,1) weight=weight * ALBGR X=XB SS=SP ENDIF IVECT0IPH= IVECT0IPH + IVECTk IMOVE=IMOVE+1 END DO !! IMOVE LOOP !!!!!!!!!!!!! IVECTCHECK=SUM(IVECT0IPH) IF(ISNAN(IVECTCHECK)) THEN IPH=IPH write(87,*) IPH, 'ISNAN(IVECTkCHECK)' ELSE IPH=IPH+1 IVECT0=IVECT0+IVECT0IPH ENDIF IF ( MOD(ALOG10(FLOAT(IPH)),1.) .EQ. 0. .OR. IPH .EQ. NPH) THEN OPEN(95,FILE=TRIM(FILENAME), ACCESS='APPEND', STATUS='UNKNOWN') DO ISUN=1,NSUN WRITE(95,'(I15,20ES15.6)') IPH, TAU, SSA, ALBGR, COSPOL0, COSPOLSUN(ISUN), & AZIMUTH(ISUN), IVECT0(1:4,ISUN) / FLOAT(IPH) ENDDO CLOSE(95) ENDIF ENDDO !! IPH LOOP !!!!!!!!!!!!! END PROGRAM PBMC SUBROUTINE SETUP implicit none COMMON /PARAMETERS1/ NPH, NSUN, NTHETAD, NTHETAI, & ex, ey, ez, WEIGHTMIN, TAU, GAMMA, SSA, ALBGR, DXSLAB, & X0, COSPOL0, SS0, COSPOLSUN, AZIMUTH, SSSUN, FLUXSUN, & MM0, THETAD, DTHETAD, THETAI, FFI, DFFI INTEGER, PARAMETER:: NSUNMAX=10, NTHETADMAX= 1800, NTHETAIMAX=4000 REAL, PARAMETER:: PI=4.*atan(1.) INTEGER NSUN, NTHETAD, NTHETAI REAL ex(1:3), ey(1:3), ez(1:3) REAL COSPOL0, AZI0 INTEGER ISUN REAL COSPOLSUN(1:NSUNMAX), AZIMUTH(1:NSUNMAX), AZIpSUN REAL SS0(1:3), SSSUN(1:3,NSUNMAX) INTEGER NPH REAL WEIGHTMIN, DXSLAB, TAU, SSA, ALBGR REAL GAMMA REAL x0(3) REAL FLUXSUN(4) INTEGER ITHETAD, ITHETAI REAL THETAD(0:NTHETADMAX), DTHETAD, MUD(0:NTHETADMAX) REAL A1(0:NTHETADMAX), B1(0:NTHETADMAX), A3(0:NTHETADMAX), B2(0:NTHETADMAX) REAL MM0(4,4,0:NTHETADMAX) REAL FFD(0:NTHETADMAX),FFI(0:NTHETAIMAX),DFFI,THETAI(0:NTHETAIMAX) CHARACTER(LEN=20) FILESCATMATRIX DATA ex / 1.,0.,0. /, ey / 0.,1.,0. /, ez / 0.,0.,1. / DATA FLUXSUN / pi,0.,0.,0./ OPEN(99,FILE='./INPUT.dat') READ(99,*) COSPOL0 AZI0=00. ! TO HAVE OBSERVER IN THE XZ PLANE IF(COSPOL0 .EQ. 1.) COSPOL0=0.99999999D0 SS0= COSPOL0 * ex + & SQRT(1.-COSPOL0**2.) * ( COS(AZI0*pi/180D0) * ez - SIN(AZI0*pi/180D0) * ey ) READ(99,*) NSUN ! IF (NSUN .NE. NSUNMAX) THEN ! WRITE(*,*) 'NSUN .NE. NSUNMAX' ! STOP ! ENDIF DO ISUN=1,NSUN READ(99,*) COSPOLSUN(ISUN), AZIMUTH(ISUN) AZIpSUN= 180D0 + AZI0 - AZIMUTH(ISUN) SSSUN(1:3,ISUN)= COSPOLSUN(ISUN) *ex + & SQRT(1.-COSPOLSUN(ISUN)**2.) * ( COS(AZIpSUN*pi/180D0)*ez - SIN(AZIpSUN*pi/180D0)*ey ) ENDDO SSSUN=-SSSUN READ(99,*) NPH READ(99,*) WEIGHTMIN READ(99,*) DXSLAB READ(99,*) TAU READ(99,*) SSA READ(99,*) ALBGR READ(99,*) FILESCATMATRIX x0(1) = DXSLAB x0(2:3)= 0. GAMMA=TAU / DXSLAB READ(99,*) NTHETAD IF (NTHETAD .NE. NTHETADMAX) THEN WRITE(*,*) 'NTHETAD .NE. NTHETADMAX' STOP ENDIF DTHETAD= pi / FLOAT(NTHETAD) OPEN(98,FILE=FILESCATMATRIX) DO ITHETAD=NTHETAD,0,-1 READ(98,*) THETAD(ITHETAD),& A1(ITHETAD),A3(ITHETAD),B1(ITHETAD),B2(ITHETAD) ENDDO CLOSE(98) THETAD=THETAD * pi /180. MUD=COS(THETAD) ! FOR RAYLEIGH SCATTERING ATMOSPHERES, UNCOMMENT THE FOLLOWING 9 ROWS ! DTHETAD= pi / FLOAT(NTHETAD) ! DO ITHETAD=0,NTHETAD ! THETAD(ITHETAD)=pi - FLOAT(ITHETAD) * DTHETAD ! ENDDO ! MUD=COS(THETAD) ! A1(:)=3./4. * ( 1. + MUD(:)**2.) ! B1(:)=3./4. * (-1. + MUD(:)**2.) ! A3(:)=3./2. * MUD(:) ! B2(:)=0. MM0=0. MM0(1,1,:)= A1 MM0(1,2,:)= B1 MM0(2,1,:)= B1 MM0(2,2,:)= A1 MM0(3,3,:)= A3 MM0(3,4,:)= B2 MM0(4,3,:)=-B2 MM0(4,4,:)= A3 MM0=MM0 / 4. / PI READ(99,*) NTHETAI IF (NTHETAI .NE. NTHETAIMAX) THEN WRITE(*,*) 'NTHETAI .NE. NTHETAIMAX' STOP ENDIF CLOSE(99) CALL FFDGENERATOR( A1, NTHETAD, THETAD, DTHETAD, FFD) ! PRODUCES FFD CALL THETAIGENERATOR( NTHETAD, THETAD, FFD, NTHETAI, FFI, DFFI, THETAI ) ! PRODUCES THETAI END SUBROUTINE SETUP SUBROUTINE FFDGENERATOR( A1, NTHETAD, THETAD, DTHETAD, FFD) ! EXTENDED TRAPEZOIDAL RULE implicit none INTEGER ITHETAD, NTHETAD REAL A1(0:NTHETAD), THETAD(0:NTHETAD), FFD(0:NTHETAD), DTHETAD FFD(0)=0. DO ITHETAD=1, NTHETAD FFD(ITHETAD)= FFD(ITHETAD-1) + & ( A1(ITHETAD-1) * SIN(THETAD(ITHETAD-1)) + A1(ITHETAD) * SIN(THETAD(ITHETAD)) ) ENDDO FFD= FFD * DTHETAD / 2. / 2. FFD(NTHETAD)=1. END SUBROUTINE FFDGENERATOR SUBROUTINE THETAIGENERATOR( NTHETAD, THETAD, FFD, NTHETAI, FFI, DFFI, THETAI ) implicit none INTEGER NTHETAD, NTHETAI, ITHETAI REAL FFD(0:NTHETAD), THETAD(0:NTHETAD), DER2(0:NTHETAD), FFI(0:NTHETAI), DFFI, THETAI(0:NTHETAI), DUMMY DFFI=1./FLOAT(NTHETAI) DO ITHETAI=0,NTHETAI FFI(ITHETAI)= FLOAT(ITHETAI) * DFFI ENDDO CALL SPLINE(FFD(0:NTHETAD),THETAD(0:NTHETAD),NTHETAD+1,2.E30,2.E30,DER2(0:NTHETAD)) DO ITHETAI=0,NTHETAI CALL SPLINT(FFD(0:NTHETAD),THETAD(0:NTHETAD),DER2(0:NTHETAD),NTHETAD+1, & FFI(ITHETAI),THETAI(ITHETAI),DUMMY) ENDDO END SUBROUTINE THETAIGENERATOR SUBROUTINE FINDXB( X, SS, XB, AXXB ) implicit none COMMON /PARAMETERS1/ NPH, NSUN, NTHETAD, NTHETAI, & ex, ey, ez, WEIGHTMIN, TAU, GAMMA, SSA, ALBGR, DXSLAB, & X0, COSPOL0, SS0, COSPOLSUN, AZIMUTH, SSSUN, FLUXSUN, & MM0, THETAD, DTHETAD, THETAI, FFI, DFFI INTEGER, PARAMETER:: NSUNMAX=10, NTHETADMAX= 1800, NTHETAIMAX=4000 INTEGER NPH, NSUN, NTHETAD, NTHETAI REAL ex(1:3), ey(1:3), ez(1:3) REAL WEIGHTMIN, TAU, GAMMA, SSA, ALBGR, DXSLAB REAL X0(1:3), COSPOL0, SS0(1:3) REAL COSPOLSUN(1:NSUNMAX), AZIMUTH(1:NSUNMAX), SSSUN(1:3,1:NSUNMAX), FLUXSUN(1:4) REAL MM0(1:4,1:4,0:NTHETADMAX), THETAD(0:NTHETADMAX), DTHETAD REAL THETAI(0:NTHETAIMAX), FFI(0:NTHETAIMAX), DFFI REAL X(1:3), SS(1:3), XB(1:3), AXXB REAL MUX, LL MUX=SS(1) IF (MUX .GT. 0.) THEN ! -SS intersects the surface LL=(X(1)- 0.)/MUX XB= X-SS*LL XB(1)=0. ELSEIF (MUX .LT. 0.) THEN ! -SS intersects the top of atmosphere LL=(X(1)-DXSLAB)/MUX XB= X-SS*LL XB(1)=DXSLAB ELSEIF (MUX .EQ. 0.) THEN ! -SS is parallel to yz plane WRITE(*,*) 'MUX = 0 @ FINDXB' STOP ENDIF AXXB=1. - exp(-LL*GAMMA) END SUBROUTINE FINDXB SUBROUTINE LLMMLLA( IPH, IDUM, HH, SS, SP, PP ) implicit none COMMON /PARAMETERS1/ NPH, NSUN, NTHETAD, NTHETAI, & ex, ey, ez, WEIGHTMIN, TAU, GAMMA, SSA, ALBGR, DXSLAB, & X0, COSPOL0, SS0, COSPOLSUN, AZIMUTH, SSSUN, FLUXSUN, & MM0, THETAD, DTHETAD, THETAI, FFI, DFFI INTEGER, PARAMETER:: NSUNMAX=10, NTHETADMAX= 1800, NTHETAIMAX=4000 INTEGER NPH, NSUN, NTHETAD, NTHETAI, IPH REAL ex(1:3), ey(1:3), ez(1:3) REAL WEIGHTMIN, TAU, GAMMA, SSA, ALBGR, DXSLAB REAL X0(1:3), COSPOL0, SS0(1:3) REAL COSPOLSUN(1:NSUNMAX), AZIMUTH(1:NSUNMAX), SSSUN(1:3,1:NSUNMAX), FLUXSUN(1:4) REAL MM0(1:4,1:4,0:NTHETADMAX), THETAD(0:NTHETADMAX), DTHETAD REAL THETAI(0:NTHETAIMAX), FFI(0:NTHETAIMAX), DFFI INTEGER IDUM REAL HH(1:4,1:4), PP(1:4,1:4), SS(1:3), SP(1:3) REAL MM(1:4,1:4) REAL MYRAN REAL ETAA, THETAA, COSTHETAA REAL ZETAA, YZETAA, PHIA, gg REAL m11, m12, qq, uu INTEGER REJ REAL, PARAMETER:: PI=4.*atan(1.) REAL, PARAMETER:: MACHPR=1.D-8 REAL e1(1:3), e2(1:3), e3(1:3) REAL ii, ip REAL NORM, DOTP REAL COSAA, COSAP, SINAP, THETACORR, COSIP REAL LLkk(1:4,1:4), LLkp(1:4,1:4) INTEGER INDL, INDR ETAA=MYRAN(IDUM) IF(ETAA .EQ. 0.) ETAA=MYRAN(IDUM) IF(ETAA .EQ. 1.) ETAA=MYRAN(IDUM) CALL FINDTHETAA(ETAA, NTHETAI, THETAI, FFI, DFFI, THETAA) !! INVERSE: GIVEN ETAA --> FIND THETAA=THETAI(FFI) CALL MMMATRIX(THETAA, MM) !!! m11=MM(1,1) m12=MM(1,2) qq=HH(1,2)/HH(1,1) uu=HH(1,3)/HH(1,1) REJ=1 DO WHILE (REJ .EQ. 1) ZETAA=MYRAN(IDUM) PHIA =2. * PI * ZETAA YZETAA=2. * MYRAN(IDUM) gg= 1. + m12/m11 * ( qq * cos(2.*PHIA) - uu * sin(2.*PHIA) ) IF (YZETAA .LE. gg) REJ=0 END DO IF (PHIA .LT. PI) THEN ii= PI - PHIA ELSEIF (PHIA .GE. PI) THEN ii= 2. * PI - PHIA ENDIF e3=SS/NORM(SS) CALL CROSSP(e3,ez,e2) e2=e2/NORM(e2) CALL CROSSP(e2,e3,e1) e1=e1/NORM(e1) !!! SP=COS(THETAA)*e3 + SIN(THETAA) * (COS(PHIA)*e1 + SIN(PHIA)*e2) SP=SP/NORM(SP) !!! COSAA=DOTP(SS,ez) COSAP=DOTP(SP,ez) SINAP=SQRT(1. - COSAP**2.) IF (PHIA .LT. PI) THEN THETACORR= 2.* PI - THETAA ELSEIF (PHIA .GE. PI) THEN THETACORR=THETAA ENDIF COSIP=(COSAA - COSAP * COS(THETACORR)) / (SINAP * SIN(THETACORR)) IF(COSIP .GE. 1.) THEN write(87,*) 'COSIP GE 1', COSIP COSIP=1. ! pause ELSEIF(COSIP .LE. -1.) THEN write(87,*) 'COSIP LE -1', COSIP COSIP=-1. ! pause ENDIF IP=ACOS(COSIP) CALL ROTATION (LLkk, PI-ii) CALL ROTATION (LLkp, -ip) PP=MATMUL(LLkk,MM) PP=MATMUL(PP,LLkp) END SUBROUTINE LLMMLLA SUBROUTINE LLMMLLB( IDUM, SP, PP ) implicit none COMMON /PARAMETERS1/ NPH, NSUN, NTHETAD, NTHETAI, & ex, ey, ez, WEIGHTMIN, TAU, GAMMA, SSA, ALBGR, DXSLAB, & X0, COSPOL0, SS0, COSPOLSUN, AZIMUTH, SSSUN, FLUXSUN, & MM0, THETAD, DTHETAD, THETAI, FFI, DFFI INTEGER, PARAMETER:: NSUNMAX=10, NTHETADMAX= 1800, NTHETAIMAX=4000 INTEGER NPH, NSUN, NTHETAD, NTHETAI REAL ex(1:3), ey(1:3), ez(1:3) REAL WEIGHTMIN, TAU, GAMMA, SSA, ALBGR, DXSLAB REAL X0(1:3), COSPOL0, SS0(1:3) REAL COSPOLSUN(1:NSUNMAX), AZIMUTH(1:NSUNMAX), SSSUN(1:3,1:NSUNMAX), FLUXSUN(1:4) REAL MM0(1:4,1:4,0:NTHETADMAX), THETAD(0:NTHETADMAX), DTHETAD REAL THETAI(0:NTHETAIMAX), FFI(0:NTHETAIMAX), DFFI INTEGER IDUM REAL PP(1:4,1:4), SP(1:3) REAL MYRAN REAL ETAB, THETAB, COSTHETAB REAL ZETAB, PHIB REAL NORM REAL, PARAMETER:: PI=4.*atan(1.) ETAB=MYRAN(IDUM) COSTHETAB=SQRT(ETAB) THETAB=ACOS(COSTHETAB) !!! ZETAB=MYRAN(IDUM) PHIB =2. * PI * ZETAB !!! SP=COS(THETAB)*ex + SIN(THETAB) * (COS(PHIB)*ez + SIN(PHIB)*ey) SP=-SP/NORM(SP) !!! PP=0. PP(1,1)=1. END SUBROUTINE LLMMLLB SUBROUTINE FINDTHETAA(ETAA, NTHETAI, THETAI, FFI, DFFI, THETAA) implicit none INTEGER NTHETAI, INDL, INDR REAL ETAA, THETAI(0:NTHETAI), FFI(0:NTHETAI), DFFI, THETAA, COSTHETAA REAL, PARAMETER:: MACHPR=1.D-8 REAL, PARAMETER:: PI=4.*atan(1.) INDL=INT(ETAA/DFFI) INDR=INDL+1 COSTHETAA=COS(THETAI(INDL)) + & (COS(THETAI(INDR))-COS(THETAI(INDL))) / (FFI(INDR)-FFI(INDL)) * & (ETAA-FFI(INDL)) THETAA=ACOS(COSTHETAA) IF(THETAA .LE. 0.) THETAA=PI* MACHPR IF(THETAA .GE. PI) THETAA=PI*(1.-MACHPR) END SUBROUTINE FINDTHETAA SUBROUTINE MMMATRIX(THETA, MM) implicit none COMMON /PARAMETERS1/ NPH, NSUN, NTHETAD, NTHETAI, & ex, ey, ez, WEIGHTMIN, TAU, GAMMA, SSA, ALBGR, DXSLAB, & X0, COSPOL0, SS0, COSPOLSUN, AZIMUTH, SSSUN, FLUXSUN, & MM0, THETAD, DTHETAD, THETAI, FFI, DFFI INTEGER, PARAMETER:: NSUNMAX=10, NTHETADMAX= 1800, NTHETAIMAX=4000 INTEGER NPH, NSUN, NTHETAD, NTHETAI REAL ex(1:3), ey(1:3), ez(1:3) REAL WEIGHTMIN, TAU, GAMMA, SSA, ALBGR, DXSLAB REAL X0(1:3), COSPOL0, SS0(1:3) REAL COSPOLSUN(1:NSUNMAX), AZIMUTH(1:NSUNMAX), SSSUN(1:3,1:NSUNMAX), FLUXSUN(1:4) REAL MM0(1:4,1:4,0:NTHETADMAX), THETAD(0:NTHETADMAX), DTHETAD REAL THETAI(0:NTHETAIMAX), FFI(0:NTHETAIMAX), DFFI REAL, PARAMETER:: PI=4.*atan(1.) INTEGER INDL, INDR REAL MM(1:4,1:4), THETA INDL=NTHETAD-INT(THETA/DTHETAD) INDR=INDL-1 MM(1:4,1:4)= MM0(1:4,1:4,INDL) + & (MM0(1:4,1:4,INDR)-MM0(1:4,1:4,INDL))/(COS(THETAD(INDR))-COS(THETAD(INDL))) * & (COS(THETA)-COS(THETAD(INDL))) END SUBROUTINE MMMATRIX SUBROUTINE CONTRIBSUNA( X, SS, GEOBIAS, weight, HH, IVECTk ) implicit none COMMON /PARAMETERS1/ NPH, NSUN, NTHETAD, NTHETAI, & ex, ey, ez, WEIGHTMIN, TAU, GAMMA, SSA, ALBGR, DXSLAB, & X0, COSPOL0, SS0, COSPOLSUN, AZIMUTH, SSSUN, FLUXSUN, & MM0, THETAD, DTHETAD, THETAI, FFI, DFFI INTEGER, PARAMETER:: NSUNMAX=10, NTHETADMAX= 1800, NTHETAIMAX=4000 REAL, PARAMETER:: PI=4.*atan(1.) REAL, PARAMETER:: MACHPR=1.D-8 INTEGER NPH, NSUN, NTHETAD, NTHETAI REAL ex(1:3), ey(1:3), ez(1:3) REAL WEIGHTMIN, TAU, GAMMA, SSA, ALBGR, DXSLAB REAL X0(1:3), COSPOL0, SS0(1:3) REAL COSPOLSUN(1:NSUNMAX), AZIMUTH(1:NSUNMAX), SSSUN(1:3,1:NSUNMAX), FLUXSUN(1:4) REAL MM0(1:4,1:4,0:NTHETADMAX), THETAD(0:NTHETADMAX), DTHETAD REAL THETAI(0:NTHETAIMAX), FFI(0:NTHETAIMAX), DFFI INTEGER ISUN REAL X(1:3), SS(1:3), GEOBIAS, weight, HH(1:4,1:4), IVECTk(1:4,1:NSUNMAX) REAL XSUN(1:3), AXXSUN, TXXSUN, iiSUN, ipSUN REAL LLkkSUN(1:4,1:4), LLkpSUN(1:4,1:4), MMSUN(1:4,1:4), PPSUN(1:4,1:4) REAL DOTP DO ISUN=1,NSUN CALL FINDXB( X, SSSUN(1:3,ISUN), XSUN, AXXSUN ) TXXSUN=1. - AXXSUN CALL METHOD2( SS, SSSUN(1:3,ISUN), iiSUN, ipSUN ) CALL ROTATION(LLkkSUN, PI-iiSUN) CALL ROTATION(LLkpSUN, -ipSUN) CALL MMMATRIX(ACOS(DOTP(SS,SSSUN(1:3,ISUN))),MMSUN) PPSUN=MATMUL(LLkkSUN,MMSUN) PPSUN=MATMUL(PPSUN,LLkpSUN) IVECTk(1:4,ISUN)= weight * MATMUL(MATMUL(HH,PPSUN), FLUXSUN ) * TXXSUN * SSA * GEOBIAS ENDDO END SUBROUTINE CONTRIBSUNA SUBROUTINE CONTRIBSUNB( X, weight, HH, IVECTk ) implicit none COMMON /PARAMETERS1/ NPH, NSUN, NTHETAD, NTHETAI, & ex, ey, ez, WEIGHTMIN, TAU, GAMMA, SSA, ALBGR, DXSLAB, & X0, COSPOL0, SS0, COSPOLSUN, AZIMUTH, SSSUN, FLUXSUN, & MM0, THETAD, DTHETAD, THETAI, FFI, DFFI INTEGER, PARAMETER:: NSUNMAX=10, NTHETADMAX= 1800, NTHETAIMAX=4000 REAL, PARAMETER:: PI=4.*atan(1.) REAL, PARAMETER:: MACHPR=1.D-8 INTEGER NPH, NSUN, NTHETAD, NTHETAI REAL ex(1:3), ey(1:3), ez(1:3) REAL WEIGHTMIN, TAU, GAMMA, SSA, ALBGR, DXSLAB REAL X0(1:3), COSPOL0, SS0(1:3) REAL COSPOLSUN(1:NSUNMAX), AZIMUTH(1:NSUNMAX), SSSUN(1:3,1:NSUNMAX), FLUXSUN(1:4) REAL MM0(1:4,1:4,0:NTHETADMAX), THETAD(0:NTHETADMAX), DTHETAD REAL THETAI(0:NTHETAIMAX), FFI(0:NTHETAIMAX), DFFI INTEGER ISUN REAL X(1:3), weight, HH(1:4,1:4), IVECTk(1:4,1:NSUNMAX) REAL XSUN(1:3), AXXSUN, TXXSUN REAL DOTP DO ISUN=1,NSUN CALL FINDXB( X, SSSUN(1:3,ISUN), XSUN, AXXSUN ) TXXSUN=1. - AXXSUN IVECTk(1:4,ISUN)= weight * MATMUL(HH,FLUXSUN) * TXXSUN * ALBGR / PI * DOTP(-ex,SSSUN(1:3,ISUN)) ENDDO END SUBROUTINE CONTRIBSUNB SUBROUTINE METHOD2( SS, SP, ii, ip ) implicit none COMMON /PARAMETERS1/ NPH, NSUN, NTHETAD, NTHETAI, & ex, ey, ez, WEIGHTMIN, TAU, GAMMA, SSA, ALBGR, DXSLAB, & X0, COSPOL0, SS0, COSPOLSUN, AZIMUTH, SSSUN, FLUXSUN, & MM0, THETAD, DTHETAD, THETAI, FFI, DFFI INTEGER, PARAMETER:: NSUNMAX=10, NTHETADMAX= 1800, NTHETAIMAX=4000 INTEGER NPH, NSUN, NTHETAD, NTHETAI REAL ex(1:3), ey(1:3), ez(1:3) REAL WEIGHTMIN, TAU, GAMMA, SSA, ALBGR, DXSLAB REAL X0(1:3), COSPOL0, SS0(1:3) REAL COSPOLSUN(1:NSUNMAX), AZIMUTH(1:NSUNMAX), SSSUN(1:3,1:NSUNMAX), FLUXSUN(1:4) REAL MM0(1:4,1:4,0:NTHETADMAX), THETAD(0:NTHETADMAX), DTHETAD REAL THETAI(0:NTHETAIMAX), FFI(0:NTHETAIMAX), DFFI REAL SS(1:3), SP(1:3), ii, ip REAL VMERID(1:3), VSCATT(1:3) REAL DOTP, NORM CALL CROSSP(SS,ez,VMERID) VMERID=VMERID/NORM(VMERID) CALL CROSSP(SS,SP,VSCATT) VSCATT=VSCATT/NORM(VSCATT) IF (DOTP(VSCATT,ez) .LT. 0.) VSCATT=-VSCATT ii=ACOS(DOTP(VMERID,VSCATT)) !!! CALL CROSSP(SP,ez,VMERID) VMERID=VMERID/NORM(VMERID) CALL CROSSP(SP,SS,VSCATT) VSCATT=VSCATT/NORM(VSCATT) IF (DOTP(VSCATT,ez) .GT. 0.) VSCATT=-VSCATT ip=ACOS(DOTP(VMERID,VSCATT)) END SUBROUTINE METHOD2 SUBROUTINE seed(idum) implicit none integer idum ,i,da_ti_zo(8) character(8) date character(10) time character(5) zone idum=0 call date_and_time(date,time,zone,da_ti_zo) do i=1,8 idum=idum+da_ti_zo(i) enddo idum=-idum ! idum=-1 ! THIS LINE MUST BE COMMENTED TO ALLOW THE SEEED TO EFFECTIVELY CHANGE OVER TIME open(99,file='seed.out') write(99,*) 'SEED EFFECTIVELY USED IN THE CALCULATIONS=', idum write(*,*) '' write(*,*) 'SEED EFFECTIVELY USED IN THE CALCULATIONS=', idum ! write(*,*) 'LINE ABOVE WITH idum=-1 SHOULD EVENTUALLY BE COMMENTED' write(*,*) '' close(99) END SUBROUTINE seed SUBROUTINE CROSSP (A, B, C) implicit none REAL A(1:3), B(1:3), C(1:3) C(1) = A(2)*B(3) - A(3)*B(2) C(2) = A(3)*B(1) - A(1)*B(3) C(3) = A(1)*B(2) - A(2)*B(1) RETURN END SUBROUTINE CROSSP SUBROUTINE ROTATION (LL, ii) implicit none REAL LL(1:4,1:4), ii LL(1:4,1:4)=0. LL(1,1)=1. LL(2,2)=cos(2.*ii) LL(3,3)=LL(2,2) LL(2,3)=sin(2.*ii) LL(3,2)=-LL(2,3) LL(4,4)=1. RETURN END SUBROUTINE ROTATION REAL FUNCTION NORM(a) implicit none REAL a(1:3) NORM= SQRT( a(1)**2 + a(2)**2 + a(3)**2 ) RETURN END FUNCTION NORM REAL FUNCTION DOTP(a,b) implicit none REAL a(1:3), b(1:3) DOTP= a(1)*b(1) + a(2)*b(2) + a(3)*b(3) RETURN END FUNCTION DOTP SUBROUTINE spline(x,y,n,yp1,ypn,y2) ! FROM NUMERICAL RECIPES IN FORTRAN 77. PRESS ET AL. ISBN 0 521 43064 X INTEGER n,NMAX REAL yp1,ypn,x(n),y(n),y2(n) PARAMETER (NMAX=500000) !! Given arrays x(1:n) and y(1:n) containing a tabulated function, i.e., y_i=f(x_i)m with !! x1