C ============================================================================= C simple driver for batch execution: hab_rt IMPLICIT REAL*4 (A-H,O-Z) REAL*8 THERMO_MSIS(7,61) REAL*8 D_EXO,FLUX,D_MAX,ALT_JNT,ALT_MAX,T_EXO INTEGER LINE_LABEL, IAPH, IGEO, MOD_MSIS DIMENSION AP(7) COMMON/PARA_MSIS/ GLAT,GLONG,IDAY,IYEAR,UT,STL,IAPH,AP,F107,F107A COMMON/PARA_EXOS/ IGEO,SATT,SATD,SCALN,ADJT COMMON/H_PROFILE/ D_EXO,FLUX,D_MAX,ALT_JNT,ALT_MAX,T_EXO,MOD_MSIS COMMON/FINE_GRID/ THERMO_MSIS,BASE,TOP,ITHERM C GATHERING INFORMATION C basic flags READ(05,*) LINE_LABEL IF ((LINE_LABEL.NE.2).AND.(LINE_LABEL.NE.3)) THEN WRITE(06,'(A)') 'INVALID BALMER SERIES TRANSITION SPECIFIED.' WRITE(06,'(A)') 'MUST BE EITHER 2 (Ly-beta) OR 3 (Ly-gamma).' GOTO 999 END IF C MSIS parameters READ(05,*) GLAT, GLONG READ(05,*) IDAY, IYEAR READ(05,*) UT READ(05,*) IAPH, AP READ(05,*) F107, F107A C thermospheric atomic hydrogen parameters READ(05,*) MOD_MSIS, D_EXO, FLUX, D_MAX C exospheric extension parameters READ(05,*) IGEO, SATT, SATD CALL ADJUST_LL(IDAY,IYEAR,UT,GLAT,GLONG,SZA,STL) CALL GET_MSIS(LINE_LABEL) CALL H_DENSITY CALL HAB_RT(LINE_LABEL) 999 CONTINUE STOP END C ============================================================================= C xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx C ============================================================================= C | | | C H A B _ R T | ATMOSPHERE & SOURCE ARRAYS | VERSION: 6.0 | C | DISTRIBUTION VERSION | JANUARY 2002 | C | | | C ============================================================================= C | C SUBROUTINES NEEDED: | C CORONA . . . ANALYTIC EXOSPHERE MODELS | C ZSUN_CELL . . POINT-TO-SUN SPECTRAL OPACITIES AND RESULTANT RESONANCE | C LINE TRANSMISSION FUNCTION | C ZONE . . . ZONE-TO-ZONE LOS PROPAGATOR | C FAR_SIDE. . . LOCATION OF CELL PIERCE-POINTS | C LOS_PARAM . . LINE-OF-SIGHT PROPAGATION PARAMETERS | C RT_GRID . . . DENSITY & TEMPERATURE ARRAYS FOR RT-SCATTERING MATRIX | C COLM_DEN. . . VERTICAL COLUMN DENSITIES FOR INPUT MODEL ATMOSPHERE | C MATRIX_INV. . MATRIX INVERSION PROCEDURES FROM NUMERICAL RECIPES | C GAUSS . . . GAUSS-LEGENDRE POINTS & WEIGHTS | C SPEED . . . SPEED POINTS & WEIGHTS | C | C ============================================================================= C C CODE TO CALCULATE SOURCE FUNCTIONS FOR RESONANTLY SCATTERED SOLAR LYMAN BETA C & GAMMA IN A SPHERICALLY SYMMETRIC NONISOTHERMAL THERMOSPHERE & EXOSPHERE. C BASED ON ALGORITHM OF ANDERSON AND HORD [1977]. C C VERSION FOR USE IN H-alpha OR H-beta AIRGLOW MODELING: C -- NONISOTHERMAL LYMAN BETA OR LYMAN GAMMA TRANSPORT ONLY C -- "FINE SCALE" THERMOSPHERE GRID EVALUATED IN SUBROUTINE GET_MSIS C -- DEFAULT OPTION FOR MSIS [H](z): CORRECTED FOR DIFFUSIVE FLOW C WITH MESOSPHERIC SOURCE REGION C -- CHOICE OF ORIGINAL CHAMBERLAIN OR MODIFIED ANALYTIC FORMULATION C FOR EXOSPHERIC EXTENSION C -- N2 & O2 OPACITY COMBINED AS A SINGLE ABSORBER SPECIES ("O2") C -- OUTPUT FILE CONTAINING SOURCE FUNCTION ARRAYS RENAMED H_alpha.source C & H_beta.source C C ADDITIONAL UPDATES FOR THIS TEST VERSION: C -- ZSUN, ZONE & HAB_ ROUTINES REBUILT TO REMOVE REDUNDANT CODE C -- EXOSPHERIC TEMPERATURE PROFILE NO LONGER REPLACED BY CONSTANT C ----------------------------------------------------------------------------- C SECTIONS: C 1 PRELIMINARY SETUP C 2 BASIC PARAMETERS & INFORMATION C 3 PARTITION ZONES & ESTABLISH GRIDS C 4 EXTERNAL SOURCE ARRAYS C 5 CONSTRUCT HOMOGENEOUS MATRIX C 6 MATRIX INVERSION C ----------------------------------------------------------------------------- SUBROUTINE HAB_RT(LINE_LABEL) C 11111111111111111111111111111111111111111111111111111111111111111111111111111 C 111 PRELIMINARY SETUP 111111111111111111111111111111111111111111111111111 C 11111111111111111111111111111111111111111111111111111111111111111111111111111 IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (IKNT = 25) ! FIXED PARAMETER (JKNT = 18) ! MUST BE EVEN PARAMETER (MDIM = 450) ! IKNT * JKNT PARAMETER (LMU = 64) ! OPTIONS: 8 16 32 64 PARAMETER (LPHI = 8) ! MUST BE EVEN PARAMETER (LPROF = 8) ! OPTIONS: 4 8 16 PARAMETER (IEXO = 14) ! EXOBASE INDEX IN RADIAL BINNING PARAMETER (SPEEDC = 2.9979D10) ! SPEED OF LIGHT PARAMETER (BOLTZ = 1.3806D-16) ! BOLTZMANN CONSTANT PARAMETER (ABSIC = 2.654D-2) ! INTEGRATED ABSORPTION CONSTANT PARAMETER (AMASS = 1.6738D-24) ! ATOMIC MASS LOGICAL TRIP(IKNT,JKNT) CHARACTER*13 CHOOZ CHARACTER*15 HISTORY CHARACTER*13 ITORNIT REAL*4 WAVELN, FNUMBER, BRATIO, ABS_N2, ABS_O2 REAL*4 GLAT, GLONG, UT, STL, F107, F107A, AP(7) REAL*4 SATT, SATD, SCALN, ADJT, BASE, TOP DIMENSION IFIBO(IKNT) DIMENSION BRAD(IKNT+1),RADN(IKNT),I_STEP(IKNT) DIMENSION BCHI(JKNT+1),CHIN(JKNT) DIMENSION THERMO(5,61),HDNS(IKNT),O2DNS(IKNT),TATM(IKNT) DIMENSION TSOL(IKNT,JKNT),HOMO(MDIM,MDIM),SOURCE(MDIM) DIMENSION LUI(MDIM),ALUD(MDIM,MDIM),ASRC(MDIM),SRC(IKNT,JKNT) DIMENSION ODLC(2,0 : LPROF,IKNT,JKNT,2),ODO2(2,IKNT,JKNT,2), & STORE(LMU,LPHI,IKNT,JKNT),OSUN(3) DIMENSION SPDPT(LPROF),WSPD(LPROF) DIMENSION THERMO_MSIS(7,61) COMMON/PARA_MSIS/ GLAT,GLONG,IDAY,IYEAR,UT,STL,IAPH,AP,F107,F107A COMMON/PARA_EXOS/ IGEOM,SATT,SATD,SCALN,ADJT COMMON/H_PROFILE/ D_EXO,FLUX,D_MAX,ALT_JNT,ALT_MAX,T_EXO,MOD_MSIS COMMON/FINE_GRID/ THERMO_MSIS,BASE,TOP,ITHERM COMMON /EXOS/ GM,CONL,PLANETR,RBASE,RC,RUPR,RP,TEXO,DEXO,VELT, & RADPF,CENTER,ABSCSX,FTSAT,FDSAT,IGEO COMMON /NMBR/ PI,RTPI,PID2,OFFSET COMMON /GAUS04/WW04(04),XX04(04) COMMON /GAUS08/WW08(08),XX08(08) COMMON /GAUS16/WW16(16),XX16(16) COMMON /GAUS32/WW32(32),XX32(32) COMMON /GAUS64/WW64(64),XX64(64) COMMON /SPW04/ SWGT04(04),SPNT04(04) COMMON /SPW08/ SWGT08(08),SPNT08(08) COMMON /SPW16/ SWGT16(16),SPNT16(16) DATA GM / 3.9898D20/ ! G-CONSTANT * PLANET MASS DATA PLANETR / 6371.0D5 / ! MEAN PLANET RADIUS DATA OFFSET / 1.0D-6 / ! GENERIC OFFSET VALUE C GENERAL/MISC STUFF PID2 = ASIN(1.0D0) PI = 2.0D0 * PID2 RTPI = SQRT(PI) CONL = GM * AMASS / BOLTZ TSAT = DBLE(SATT) DSAT = DBLE(SATD) IGEO = IGEOM C EFFECTIVE ZERO FOR TRANSMISSION FUNCTION EFZERO = EXP(-100.0D0) C FIBONACCI SEQUENCE (USED IN SETTING UP RADIAL ZONES) II1 = 1 II2 = 2 DO II = 1,IKNT IFIBO(II) = II1 II3 = II1 + II2 II1 = II2 II2 = II3 END DO C ALTITUDE STEPS USED IN SETTING UP RADIAL ZONES I_STEP(1) = 1 DO II = 2,IKNT IF (II.LE.7) THEN I_STEP(II) = I_STEP(II-1) + II ELSE I_STEP(II) = I_STEP(II-1) + IFIBO(II-3) END IF END DO C SPEED POINT ARRAYS C (DOPPLER DISPLACEMENTS FOR NONISOTHERMAL CALCULATIONS) DO LP = 1,LPROF IF (LPROF.EQ.4) THEN SPDPT(LP) = SPNT04(LP) WSPD(LP) = SWGT04(LP) ELSE IF (LPROF.EQ.8) THEN SPDPT(LP) = SPNT08(LP) WSPD(LP) = SWGT08(LP) ELSE IF (LPROF.EQ.16) THEN SPDPT(LP) = SPNT16(LP) WSPD(LP) = SWGT16(LP) END IF END DO C 22222222222222222222222222222222222222222222222222222222222222222222222222222 C 222 BASIC PARAMETERS & INFORMATION 22222222222222222222222222222222222222 C 22222222222222222222222222222222222222222222222222222222222222222222222222222 C ----------------------------------------------------------------------------- C PARAMETERS DEFINING CASE | C | C NONIT . . . LOGICAL FLAG FOR NONISOTHERMAL EVALUATION OF SOURCE | C FUNCTIONS (STRICTLY .TRUE. & NOT LONGER A VARIABLE) | C (LPROF = NUMBER OF SPEED POINTS USED IN THE EVALUATION) | C LINE_LABEL INDEX SPECIFYING BALMER SERIES LINE OF INTEREST: | C 2 - LYMAN BETA : BALMER ALPHA | C 3 - LYMAN GAMMA : BALMER BETA | C | C BASE . . . ALTITUDE OF BASE OF RADIATING REGION (FIXED) | C TOP . . . TOP OF THERMOSPHERE OR EXOBASE (FIXED: BASE + 392 KM) | C ITHERM. . . NUMBER OF LAYERS IN REFERENCE THERMOSPHERE MODEL | C (QFCTR,QPIVOT CONTROL THE NONUNIFORM SPACING OF THE | C MSIS REFERENCE LEVELS, MUST BE SPECIFIED FOR EACH | C EMISSION OF INTEREST) | C THERMO. . . ARRAY CONTAINING "FINE GRID" THERMOSPHERIC QUANTITIES | C (ALL PURE ABSORBERS - O2 & N2 - MERGED INTO A COMPOSITE | C "O2" DENSITY) | C | C GLAT . . . GEODETIC LATITUDE OF EVALUATION (DEGREES) | C GLONG . . . GEOGRAPHIC LONGITUDE OF EVALUATION (DEGREES) | C IDAY . . . OBSERVATION UT DAY NUMBER | C IYEAR . . . OBSERVATION YEAR | C UT . . . UNIVERSAL TIME OF EVALUATION | C STL . . . LOCAL APPARENT SOLAR TIME OF EVALUATION | C IAPH . . . OPTION SWITCH FOR AP-HISTORY MODE | C IF 1, NEED TO SPECIFY AP ARRAY ELEMENTS 2-7 | C (SET AP ARRAY ELEMENTS 2-7 TO 0.0 OTHERWISE) | C AP(1) . . . DAILY (AVERAGE) AP | C F107 . . . SOLAR F(10.7) FOR PREVIOUS 24 HOURS | C F107A . . . SOLAR F(10.7) 81-DAY AVERAGE | C | C IGEO . . . FLAG FOR GEOCORONAL MODEL: | C 0 - ORIGINAL CHAMBERLAIN [1963] MODEL | C 1 - MODIFIED GEOCORONAL [BISHOP 1991] MODEL | C TSAT . . . SATELLITE COMPONENT "EXOBASE" T-PARAMETER FOR IGEO=1 | C SATELLITE CRITICAL RADIUS (RE) FOR IGEO=0 | C DSAT . . . SATELLITE COMPONENT "EXOBASE" n-PARAMETER FOR IGEO=1 | C IF VALUE < 0 IS ENTERED, EVAPORATIVE CASE IS EVALUATED | C IF VALUE = 0, BALLISTIC CASE IS EVALUATED | C NOT USED FOR IGEO=0 | C | C ----------------------------------------------------------------------------- C ----------------------------------------------------------------------------- C SELECTED EMISSION C ----------------------------------------------------------------------------- C LYMAN BETA TRANSITION PARAMETERS IF (LINE_LABEL.EQ.2) THEN CHOOZ = 'LYMAN BETA ' WAVELN = 1025.72 FNUMBER = 0.0791 BRATIO = 0.882 ABS_N2 = 0.0 ABS_O2 = 5.46E-18 OPEN (10,FILE='H_alpha.source',STATUS='UNKNOWN') ELSE IF (LINE_LABEL.EQ.3) THEN C LYMAN GAMMA TRANSITION PARAMETERS CHOOZ = 'LYMAN GAMMA ' WAVELN = 972.537 FNUMBER = 0.0290 BRATIO = 0.839 ABS_N2 = 1.04E-16 ABS_O2 = 3.90E-17 OPEN (10,FILE='H_beta.source',STATUS='UNKNOWN') END IF C NOTE: in the Lyman-gamma case, ABS_N2 is temperature dependent: C solmax 1.0E-16 cm^2 C solmin 1.5E-16 cm^2 C ----------------------------------------------------------------------------- C "FINE GRID" THERMOSPHERE MODEL & RELATED QUANTITIES C ----------------------------------------------------------------------------- C WRITE OUT RELEVANT ATMOSPHERIC INFO IF (IAPH.EQ.0) HISTORY = 'NO AP HISTORY ' IF (IAPH.EQ.1) HISTORY = 'AP HISTORY MODE' ITORNIT = 'NONISOTHERMAL' WRITE(10,10)CHOOZ,ITORNIT, & GLAT,GLONG,IDAY,IYEAR,UT,STL,AP(1),F107,F107A,HISTORY 10 FORMAT(1X,'SOURCE EVALUATION CONDITIONS: ',A,2X,A/ & 1X,F6.2,2X,F7.2,2X,I3,2X,I2,2X,5(F6.2,2X),2X,A) C THERMOSPHERE REFERENCE TABLE WRITE(10,11)BASE,TOP,ITHERM,MOD_MSIS,SCALN,ADJT,D_EXO,FLUX,D_MAX 11 FORMAT(1X,'MSIS-90 THERMOSPHERE:'/1X,0P2F9.3,2I4,0P2F9.3,1P3E12.3) 12 FORMAT(1X,5(1PE10.4,1X)) DO IBIN = 1,ITHERM+1 THERMO(1,IBIN) = THERMO_MSIS(1,IBIN) ! z (km) THERMO(2,IBIN) = THERMO_MSIS(2,IBIN) ! r (cm) THERMO(3,IBIN) = THERMO_MSIS(3,IBIN) ! T (K) THERMO(4,IBIN) = THERMO_MSIS(4,IBIN) ! [H] (cm^-3) THERMO(5,IBIN) = THERMO_MSIS(5,IBIN) ! [O2] (cm^-3) & + (ABS_N2/ABS_O2)*THERMO_MSIS(6,IBIN) ! [N2] (cm^-3) WRITE(10,12)(THERMO(JJ,IBIN),JJ=1,5) END DO C ----------------------------------------------------------------------------- C EXOBASE & EXOSPHERE INFORMATION C ----------------------------------------------------------------------------- C TEXO . . . EXOBASE TEMPERATURE C DEXO . . . EXOBASE ATOMIC HYDROGEN DENSITY C FTSAT . . . IGEO=1: SCALING FOR SATELLITE "EXOBASE" TEMPERATURE C IGEO=0: SATELLITE CRITICAL RADIUS C FDSAT . . . IGEO=1: SCALING FOR SATELLITE "EXOBASE" DENSITY C IGEO=0: NOT USED C RBASE . . . INNERMOST RADIUS FOR RT BINNING C RC . . . EXOBASE RADIUS C RUPR . . . OUTERMOST RADIUS FOR RT BINNING C RP . . . RADIUS OF THE EXOPAUSE: USES BOSSY [1983] RELATION C GIVING INTEGRATED SOLAR LYMAN ALPHA FLUX FOR SPECIFIED C F10.7 AND RELYING ON ASSESSMENT OF AJELLO ET AL 1987 C THAT LINE CENTER FLUX ALMOST NUMERICALLY THE SAME C RADPF . . . LYMAN ALPHA LINE CENTER FLUX USED IN CALCULATING RP C ----------------------------------------------------------------------------- TEXO = THERMO(3,ITHERM+1) DEXO = THERMO(4,ITHERM+1) IF (IGEO.EQ.1) THEN IF (DSAT.LT.0.0D0) THEN TSAT = TEXO DSAT = DEXO END IF FTSAT = TSAT/TEXO FDSAT = DSAT/DEXO ELSE IF (IGEO.EQ.0) THEN FTSAT = TSAT * PLANETR FDSAT = 0.0D0 DSAT = ABS(DSAT) END IF RBASE = PLANETR + DBLE(BASE)*1.0D5 RC = PLANETR + DBLE(TOP) *1.0D5 RUPR = RBASE + DBLE(I_STEP(IKNT))*1.0D5 RADPF = DBLE(2.91 * (1.0 + 0.002*(F107-65.0))) RP = SQRT(GM/0.1774D0/RADPF) C ----------------------------------------------------------------------------- C PHOTO CROSS SECTIONS & BRANCHING RATIOS C ----------------------------------------------------------------------------- C (DOPPLER PROFILE AT EXOBASE TEMPERATURE) C VELT . . . MOST PROBABLE SPEED C CENTER . . . LINE CENTER SCATTERING CROSS SECTION C ABSCSX . . . ABSORBER PHOTOABSORPTION CROSS SECTION C BRANCH . . . REEMISSION BRANCHING RATIO C ----------------------------------------------------------------------------- VELT = SQRT(2.0d0 * BOLTZ / AMASS) * SQRT(TEXO) CENTER = ABSIC * DBLE(FNUMBER*WAVELN) * 1.0d-8 / RTPI / VELT ABSCSX = DBLE(ABS_O2) BRANCH = DBLE(BRATIO) C ----------------------------------------------------------------------------- C WRITE OUT EXOSPHERE & OPTICAL QUANTITIES C ----------------------------------------------------------------------------- WRITE(10,13)IGEO,TEXO,DEXO,TSAT,DSAT,RBASE,RC,RUPR,RP, & RADPF,VELT,CENTER,ABSCSX,BRANCH 13 FORMAT(1X,'EXOSPHERE & OPTICAL QUANTITIES:'/ & I3,4(1X,1PE10.4)/10(1X,1PE10.4)) C ----------------------------------------------------------------------------- C EVALUATE ZENITH COLUMN DENSITIES (FOR REFERENCE ONLY) C ----------------------------------------------------------------------------- CALL COLM_DEN(ITHERM,THERMO, & TCOLM_H,TCOLM_O2,TCOLM_TOT,COLM_EXO) C WRITE OUT WRITE(10,14) TCOLM_H, TCOLM_O2, TCOLM_TOT, COLM_EXO 14 FORMAT(1X,'COLUMN DENSITIES:'/1X,4(1PE10.4,2X)) C 33333333333333333333333333333333333333333333333333333333333333333333333333333 C 333 PARTITION ZONES & ESTABLISH GRIDS 33333333333333333333333333333333333 C 33333333333333333333333333333333333333333333333333333333333333333333333333333 C ----------------------------------------------------------------------------- C CELL GRIDS USED IN RT-CALCS C ----------------------------------------------------------------------------- C RADIAL GRID POINTS (BASED ON FIBONACCI SEQUENCE) C II . . . INDEX OVER RADIUS BRAD(1) = RBASE DO II = 1,IKNT BRAD(II+1) = BRAD(1) + DBLE(I_STEP(II))*1.0D5 END DO C _____ DEFINE "NOMINAL" CENTROID RADII DO II = 1,IKNT RADN(II) = (BRAD(II) + BRAD(II+1)) / 2.0D0 END DO C SOLAR ANGLE GRID POINTS C JJ . . . INDEX OVER SOLAR ANGLE CHI C _____ SUBSOLAR HEMISPHERE BCHI(1) = 0.0D0 BCHI(2) = ASIN(RBASE/RUPR) XLOC = COS(BCHI(2)) STEP = XLOC / DBLE(JKNT/2 - 1) DO JJ = 3,JKNT/2 + 1 XLOC = XLOC - STEP BCHI(JJ) = ACOS(XLOC) END DO C _____ FOLDING OVER FOR ANTISOLAR HEMISPHERE DO JJ = JKNT/2 + 1,JKNT JR = JKNT + 1 - JJ BCHI(JJ+1) = PI - BCHI(JR) END DO C _____ DEFINE "NOMINAL" CENTROID SOLAR ANGLES DO JJ = 1,JKNT CHIN(JJ) = (BCHI(JJ) + BCHI(JJ+1)) / 2.0D0 END DO C WRITE OUT RADIAL AND SOLAR ANGLE GRIDS WRITE(10,15)BRAD 15 FORMAT(1X,'RADIAL GRID INFO:'/(10(1X,1PE10.4))) WRITE(10,16)BCHI 16 FORMAT(1X,'SOLAR ANGLE GRID INFO:'/(10(1X,1PE10.4))) C WRITE OUT "NOMINAL" CENTROID RADII AND SOLAR ANGLES WRITE(10,17)RADN 17 FORMAT(1X,'NOMINAL CENTROID RADII:'/(10(1X,1PE10.4))) WRITE(10,18)CHIN 18 FORMAT(1X,'NOMINAL CENTROID SOLAR ANGLES:'/(10(1X,1PE10.4))) C ----------------------------------------------------------------------------- C ZONE DENSITIES & TEMPERATURES: SPHERICALLY SYMMETRIC C ----------------------------------------------------------------------------- CALL RT_GRID(ITHERM,THERMO,IEXO,IKNT,BRAD,HDNS,O2DNS,TATM) C WRITE OUT DENSITY & TEMPERATURE ARRAYS WRITE(10,21)HDNS 21 FORMAT(1X,'HYDROGEN DENSITIES:'/(10(1X,1PE10.4))) WRITE(10,22)O2DNS 22 FORMAT(1X,'MOLECULAR OXYGEN DENSITIES:'/(10(1X,1PE10.4))) WRITE(10,23)TATM 23 FORMAT(1X,'TEMPERATURES:'/(10(1X,1PE10.4))) C 44444444444444444444444444444444444444444444444444444444444444444444444444444 C 444 EXTERNAL SOURCE ARRAYS 4444444444444444444444444444444444444444444444 C 44444444444444444444444444444444444444444444444444444444444444444444444444444 C ----------------------------------------------------------------------------- C INITIALIZE ARRAYS C ----------------------------------------------------------------------------- C CONVERT THERMOSPHERIC DENSITIES TO DENSITY LOGARITHMS FOR USE IN C LINEAR-LOGARITHMIC INTERPOLATIONS DO IBIN = 1,ITHERM+1 THERMO(4,IBIN) = LOG(THERMO(4,IBIN)) THERMO(5,IBIN) = LOG(THERMO(5,IBIN)) END DO DO II = 1,IKNT DO JJ = 1,JKNT MROW = (JJ-1)*IKNT + II TSOL(II,JJ) = 0.0D0 SOURCE(MROW) = 0.0D0 END DO END DO C ----------------------------------------------------------------------------- C CALCULATE SOLAR TRANSMISSION FUNCTIONS AT CELL EFFECTIVE CENTROIDS C ----------------------------------------------------------------------------- C NOTES: C - RBASE=BRAD(1) DEFINES THE PLANETARY SHADOW TERMINATOR RADIUS C LOS PARALLEL TO PLANET-SUN AXIS C - FOR SOLAR IRRADIANCE IN CURRENT GEOMETRY THETA=CHIS AND PHI IS SET C TO A SMALL POSITIVE VALUE C C RADS . . . RADIAL DISTANCE FOR EFFECTIVE CELL CENTROID C CHIS . . . SOLAR ANGLE FOR EFFECTIVE CELL CENTROID C VOLTOT . . MEASURE OF TOTAL CELL VOLUME C VOLILL . . MEASURE OF ILLUMINATED CELL VOLUME C FRAC . . . ILLUMINATED FRACTION OF VOLUME C ----------------------------------------------------------------------------- SOLPHI = 0.001D0 DO II = 1,IKNT DO JJ = 1,JKNT MROW = (JJ-1)*IKNT + II BR1 = BRAD(II) BR2 = BRAD(II+1) BCX1 = COS(BCHI(JJ)) BCX2 = COS(BCHI(JJ+1)) BSX1 = SIN(BCHI(JJ)) BSX2 = SIN(BCHI(JJ+1)) VOLTOT = (BR2**3 - BR1**3) * (BCX1 - BCX2) C ILLUMINATED ZONES C _____ FULLY ILLUMINATED IF ((BCHI(JJ+1).LE.PID2).OR.(RBASE.LE.BR1*BSX2)) THEN RADS = (BRAD(II) + BRAD(II+1)) / 2.0D0 CHIS = (BCHI(JJ) + BCHI(JJ+1)) / 2.0D0 VOLILL = VOLTOT C _____ SHADOW INTERSECTS BRAD(II) BOUNDARY ELSE IF ((BR1*BSX2.LT.RBASE).AND.(RBASE.LT.BR1*BSX1)) THEN RADS = (BRAD(II) + BRAD(II+1)) / 2.0D0 XSM = PI - ASIN(RBASE/RADS) IF (BCHI(JJ+1).LT.XSM) XSM = BCHI(JJ+1) CHIS = (BCHI(JJ) + XSM) / 2.0D0 XS1 = PI - ASIN(RBASE/BR1) XS2 = PI - ASIN(RBASE/BR2) IF (BCHI(JJ+1).LT.XS2) XS2 = BCHI(JJ+1) VOLILL = (BR2**3 - BR1**3) * (BCX1 - COS(XS1)) & + (BR2**3) * (COS(XS1) - COS(XS2)) & + (RBASE**3) * (1.0D0/TAN(XS2) - 1.0D0/TAN(XS1)) C _____ SHADOW INTERSECTS BCHI(JJ) BOUNDARY ELSE IF ((BR1*BSX1.LE.RBASE).AND.(RBASE.LT.BR2*BSX1)) THEN RSM = RBASE / BSX1 RADS = (RSM + BRAD(II+1)) / 2.0D0 XSM = PI - ASIN(RBASE/RADS) IF (BCHI(JJ+1).LT.XSM) XSM = BCHI(JJ+1) CHIS = (BCHI(JJ) + XSM) / 2.0D0 XS2 = PI - ASIN(RBASE/BR2) IF (BCHI(JJ+1).LT.XS2) XS2 = BCHI(JJ+1) VOLILL = (BR2**3) * (BCX1 - COS(XS2)) & + (RBASE**3) & * (1.0D0/TAN(XS2) - 1.0D0/TAN(BCHI(JJ))) C COMPLETELY IN SHADOW ELSE RADS = (BRAD(II) + BRAD(II+1)) / 2.0D0 CHIS = (BCHI(JJ) + BCHI(JJ+1)) / 2.0D0 TSOL(II,JJ) = 0.0D0 GO TO 444 END IF FRAC = VOLILL / VOLTOT C LINE CENTER OPTICAL DEPTHS & NET TRANSMISSION FUNCTION CALL ZSUN_CELL(IKNT,JKNT,BRAD,BCHI,HDNS,O2DNS,TATM, & II,JJ,RADS,CHIS,CHIS,SOLPHI, & LPROF,SPDPT,WSPD,OSUN) TSOL(II,JJ) = FRAC * EXP(-OSUN(2)) * OSUN(3) IF (TSOL(II,JJ).LT.EFZERO) TSOL(II,JJ) = EFZERO 444 CONTINUE C STUFF INTO 1-D ARRAY SOURCE(MROW) = TSOL(II,JJ) END DO END DO C ----------------------------------------------------------------------------- C WRITE OUT SOLAR TRANSMISSION FUNCTIONS C ----------------------------------------------------------------------------- WRITE(10,31) 31 FORMAT(1X,'SOLAR TRANSMISSION FUNCTIONS:') 32 FORMAT(10(1X,1PE10.4)) DO II = 1,IKNT WRITE(10,32)(TSOL(II,JJ),JJ=1,JKNT) END DO C 55555555555555555555555555555555555555555555555555555555555555555555555555555 C 555 CONSTRUCT HOMOGENEOUS MATRIX (USING NOMINAL CENTROID POINTS) 55555555 C 55555555555555555555555555555555555555555555555555555555555555555555555555555 C ----------------------------------------------------------------------------- C INITIALIZE HOMOGENOUS MATRIX & AZIMUTH ANGLE FACTOR C ----------------------------------------------------------------------------- DO MROW = 1,MDIM DO MCOL = 1,MDIM HOMO(MROW,MCOL) = 0.0D0 END DO END DO C AZIMUTHAL WEIGHT (CONSTANT) WPHI = 2.0D0 * PI / DBLE(LPHI) C ----------------------------------------------------------------------------- C LAY OUT INTEGRALS OVER ORIENTATION, FOR EACH BIN C NOTE: INVOKING SYMMETRY IN CHI=PI/2 PLANE. JKNT MUST BE EVEN. C ----------------------------------------------------------------------------- DO 1001 II = 1,IKNT DO 1001 JJ = 1,JKNT/2 MROW = (JJ-1)*IKNT + II R1 = RADN(II) X1 = CHIN(JJ) FCTREF = TEXO/TATM(II) C _____ SUMMATION SWEEPS: THETA IS LOCAL ZENITH ANGLE AND C _____ PHI IS AZIMUTH ABOUT LOCAL ZENITH, ZERO IN SOLAR DIRECTION C _____ (INVOKING SYMMETRY BETWEEN EAST & WEST IN PHI INTEGRAL) DO 1101 KMU = 1,LMU IF (LMU.EQ.8) THEN THETA = ACOS(XX08(KMU)) ELSE IF (LMU.EQ.16) THEN THETA = ACOS(XX16(KMU)) ELSE IF (LMU.EQ.32) THEN THETA = ACOS(XX32(KMU)) ELSE IF (LMU.EQ.64) THEN THETA = ACOS(XX64(KMU)) END IF DO 1102 KPHI = 1,LPHI PHI = PI * DBLE(2*KPHI-1) / DBLE(2*LPHI) C WRITE(06,*) II, JJ, KMU, KPHI, 'START ZONE' CALL ZONE(IKNT,JKNT,BRAD,BCHI,HDNS,O2DNS,TATM, & II,JJ,R1,X1,THETA,PHI,LPROF,SPDPT, & TRIP,ODLC,ODO2,IEXIT,JEXIT) C WRITE(06,*) II, JJ, KMU, KPHI, 'END ZONE' C _____ NOW CONVERTING THE LINE-CENTER OPTICAL DEPTHS TO C _____ LINE-INTEGRATED TRANSMISSION FUNCTIONS, USING PREVIOUSLY C _____ COMPUTED INTERPOLATION TABLE. ALSO COMBINING TRANSMISSION C _____ FUNCTIONS FOR CASE WHERE ZONE IS CROSSED TWICE. C _____ SUBTRACTION IS THEN PERFORMED. C _____ HANDLING OF O2 ABSORPTION CORRECTED: 01 MAY 1991. C _____ CONVERTED TO FULLY NONISOTHERMAL EVALUATION: 26 SEP 1994. DO 1103 III = 1,IKNT DO 1103 JJJ = 1,JKNT IF (TRIP(III,JJJ)) THEN RHO = O2DNS(III)*ABSCSX/(HDNS(III)*CENTER) ABS11 = EXP(-ODO2(1,III,JJJ,1)) ABS21 = EXP(-ODO2(2,III,JJJ,1)) TAU1 = ODLC(1,0,III,JJJ,1) TAU2 = ODLC(2,0,III,JJJ,1) TLC11 = 0.0D0 TLC21 = 0.0D0 DO LP = 1,LPROF XPNT = (FCTREF - 1.0D0) * (SPDPT(LP)**2) TLC11 = TLC11 + WSPD(LP) * EXP(-XPNT) * & EXP(-ODLC(1,LP,III,JJJ,1)) TLC21 = TLC21 + WSPD(LP) * EXP(-XPNT) * & EXP(-ODLC(2,LP,III,JJJ,1)) END DO TLC11 = TLC11 * SQRT(FCTREF) * 2.0D0 / RTPI TLC21 = TLC21 * SQRT(FCTREF) * 2.0D0 / RTPI HOLDIT = ABS11*TLC11 - ABS21*TLC21 IF (RHO.GT.OFFSET) THEN SUM = 0.0D0 DO KK = 1,4 TAU = ((TAU2-TAU1)*XX04(KK)+TAU2+TAU1)/2.0D0 WTAU = (TAU2-TAU1)*WW04(KK)/2.0D0 TABS = EXP(-(ODO2(1,III,JJJ,1) + RHO*(TAU-TAU1))) FCTR = SQRT(TEXO/TATM(III)) TSCAT = 0.0D0 DO LP = 1,LPROF XPNT = (FCTR * SPDPT(LP))**2 TAULP = (TAU-TAU1) * FCTR * EXP(-XPNT) XPNT = (FCTREF - 1.0D0) * (SPDPT(LP)**2) TSCAT = TSCAT + WSPD(LP) * EXP(-XPNT) * & EXP(-(TAULP + ODLC(1,LP,III,JJJ,1))) END DO TSCAT = TSCAT * SQRT(FCTREF) * 2.0D0 / RTPI SUM = SUM + RHO * WTAU * TSCAT * TABS END DO HOLDIT = HOLDIT - SUM END IF IF (ODLC(2,0,III,JJJ,2).GT.0.0D0) THEN ABS12 = EXP(-ODO2(1,III,JJJ,2)) ABS22 = EXP(-ODO2(2,III,JJJ,2)) TAU1 = ODLC(1,0,III,JJJ,2) TAU2 = ODLC(2,0,III,JJJ,2) TLC12 = 0.0D0 TLC22 = 0.0D0 DO LP = 1,LPROF XPNT = (FCTREF - 1.0D0) * (SPDPT(LP)**2) TLC12 = TLC12 + WSPD(LP) * EXP(-XPNT) * & EXP(-ODLC(1,LP,III,JJJ,2)) TLC22 = TLC22 + WSPD(LP) * EXP(-XPNT) * & EXP(-ODLC(2,LP,III,JJJ,2)) END DO TLC12 = TLC12 * SQRT(FCTREF) * 2.0D0 / RTPI TLC22 = TLC22 * SQRT(FCTREF) * 2.0D0 / RTPI HOLDIT = HOLDIT + ABS12*TLC12 - ABS22*TLC22 IF (RHO.GT.OFFSET) THEN SUM = 0.0D0 DO KK = 1,4 TAU = ((TAU2-TAU1)*XX04(KK)+TAU2+TAU1)/2.0D0 WTAU = (TAU2-TAU1)*WW04(KK)/2.0D0 TABS = EXP(-(ODO2(1,III,JJJ,2) + RHO*(TAU-TAU1))) FCTR = SQRT(TEXO/TATM(III)) TSCAT = 0.0D0 DO LP = 1,LPROF XPNT = (FCTR * SPDPT(LP))**2 TAULP = (TAU-TAU1) * FCTR * EXP(-XPNT) XPNT = (FCTREF - 1.0D0) * (SPDPT(LP)**2) TSCAT = TSCAT + WSPD(LP) * EXP(-XPNT) * & EXP(-(TAULP + ODLC(1,LP,III,JJJ,2))) END DO TSCAT = TSCAT * SQRT(FCTREF) * 2.0D0 / RTPI SUM = SUM + RHO * WTAU * TSCAT * TABS END DO HOLDIT = HOLDIT - SUM END IF END IF IF (HOLDIT.LE.0.0D0) HOLDIT = 0.0D0 STORE(KMU,KPHI,III,JJJ) = HOLDIT ELSE STORE(KMU,KPHI,III,JJJ) = 0.0D0 END IF 1103 CONTINUE 1102 CONTINUE 1101 CONTINUE C _____ NOW TO PERFORM SUMMATIONS OVER UMU AND PHI, STORING RESULTS IN C _____ HOMOGENEOUS MATRIX DO III = 1,IKNT DO JJJ = 1,JKNT MCOL = (JJJ-1)*IKNT + III SUM1 = 0.0D0 DO KMU = 1,LMU IF (LMU.EQ.8) THEN WUMU = WW08(KMU) ELSE IF (LMU.EQ.16) THEN WUMU = WW16(KMU) ELSE IF (LMU.EQ.32) THEN WUMU = WW32(KMU) ELSE IF (LMU.EQ.64) THEN WUMU = WW64(KMU) END IF SUM2 = 0.0D0 DO KPHI = 1,LPHI SUM2 = SUM2 + WPHI * STORE(KMU,KPHI,III,JJJ) END DO SUM1 = SUM1 + WUMU * SUM2 END DO HOMO(MROW,MCOL) = SUM1 END DO END DO 1001 CONTINUE C ----------------------------------------------------------------------------- C FOLD MATRIX, TO GET CHI > PI/2 ENTRIES (SPHERICAL SYMMETRY) C ----------------------------------------------------------------------------- DO II = 1,IKNT DO JJ = 1,JKNT/2 JQ = (JKNT+1) - JJ MM = (JJ-1)*IKNT + II MQ = (JQ-1)*IKNT + II DO III = 1,IKNT DO JJJ = 1,JKNT JJQ = (JKNT+1) - JJJ NN = (JJJ-1)*IKNT + III NQ = (JJQ-1)*IKNT + III HOMO(MQ,NQ) = HOMO(MM,NN) END DO END DO END DO END DO C 66666666666666666666666666666666666666666666666666666666666666666666666666666 C 666 MATRIX INVERSION 6666666666666666666666666666666666666666666666666666 C 66666666666666666666666666666666666666666666666666666666666666666666666666666 C ----------------------------------------------------------------------------- C PREPARATIONS C ----------------------------------------------------------------------------- C RECASTING HOMOGENEOUS MATRIX INTO COEFFICIENT MATRIX DO MROW = 1,MDIM DO MCOL = 1,MDIM HOMO(MROW,MCOL) = -BRANCH*HOMO(MROW,MCOL)/(4.0D0*PI) IF (MROW.EQ.MCOL) HOMO(MROW,MCOL) = 1.0D0 + HOMO(MROW,MCOL) END DO END DO C PREPARING TO CALL INVERSION ROUTINES (COPY ARRAYS TO BE DESTROYED) DO MROW = 1,MDIM ASRC(MROW) = SOURCE(MROW) DO MCOL = 1,MDIM ALUD(MROW,MCOL) = HOMO(MROW,MCOL) END DO END DO C ----------------------------------------------------------------------------- C INVERT COEFFICIENT MATRIX C ----------------------------------------------------------------------------- CALL LUDCMP(ALUD,MDIM,LUI,D) CALL LUBKSB(ALUD,MDIM,LUI,ASRC) CALL MPROVE(HOMO,ALUD,MDIM,LUI,SOURCE,ASRC) C ----------------------------------------------------------------------------- C WRAP-UP C ----------------------------------------------------------------------------- C STUFF INTO 2-D ARRAY DO II = 1,IKNT DO JJ = 1,JKNT MROW = (JJ-1)*IKNT + II SRC(II,JJ) = ASRC(MROW) IF (SRC(II,JJ).LT.EFZERO) SRC(II,JJ) = EFZERO END DO END DO C WRITE OUT TOTAL SOURCE FUNCTIONS WRITE(10,33) 33 FORMAT(1X,'TOTAL SOURCE FUNCTIONS:') 34 FORMAT(10(1X,1PE10.4)) DO II = 1,IKNT WRITE(10,34)(SRC(II,JJ),JJ=1,JKNT) END DO CLOSE(10) C ----------------------------------------------------------------------------- WRITE(06,'(A)') 'did we get this far?' GO TO 9009 !! test jump 9009 CONTINUE RETURN END C =============================================================================