Original FORTRAN Code
Below is the original FORTRAN source code, portions of which are
replicated in the SIGMA TGM data flow graph, partcularly the code
corresponding to the subroutine ATMSETUP.
Related information
PROGRAM TITAN
C23456789012345678901234567890123456789012345678901234567890123456789012
C PROGRAM TO CALCULATE THE SCATTERING GREENHOUSE EFFECT
C ON TITAN CODE BEGUN 18 MARCH 1985
C C.P. McKAY
C.......................................................................
PARAMETER (NLAYER=30,NLEVEL=31)
PARAMETER (NSPECI=46,NSPC1I=47,NSPECV=24,NSPC1V=25)
COMMON /UBARED/ UBARI,UBARV,UBAR0
COMMON /ATM/ Z(NLEVEL),PRESS(NLEVEL),DEN(NLEVEL),TEMP(NLEVEL)
COMMON /LAPSE/ DTDP(NLAYER),CONVEQ
COMMON /GASS/ CH4(NLEVEL),XN2(NLEVEL),H2(NLEVEL),AR(NLEVEL)
& ,XMU(NLEVEL),GAS1(NLAYER),COLDEN(NLAYER)
COMMON /STRATO/ C2H2(NLAYER),C2H6(NLAYER)
COMMON /VISGAS/SOLARF(NSPECV),NTERM(NSPECV),PEXPON(NSPECV),
& ATERM(4,NSPECV),BTERM(4,NSPECV)
COMMON /AERSOL/ RADIUS(NLAYER), XNUMB(NLAYER)
& , REALI(NSPECI), XIMGI(NSPECI), REALV(NSPECV), XIMGV(NSPECV)
COMMON /CLOUD/ RADCLD(NLAYER), XNCLD(NLAYER)
& , RCLDI(NSPECI), XICLDI(NSPECI), RCLDV(NSPECV), XICLDV(NSPECV)
COMMON /TAUS/ TAUHI(NSPECI),TAUCI(NSPECI),TAUGI(NSPECI),
& TAURV(NSPECV),TAUHV(NSPECV),TAUCV(NSPECV),TAUGV(NSPECV)
COMMON /OPTICI/ DTAUI(NLAYER,NSPECI),TAUI(NLEVEL,NSPECI)
& ,WBARI(NLAYER,NSPECI), COSBI(NLAYER,NSPECI)
COMMON /SPECTI/ BWNI(NSPC1I),WNOI(NSPECI),DWNI(NSPECI)
& ,WLNI(NSPECI)
COMMON /FLUXI/ FNETI(NLEVEL), FUPI(NLEVEL,NSPECI),
& FDI(NLEVEL,NSPECI),FMNETI(NLEVEL),FMUPI(NLEVEL),FMDI(NLEVEL)
COMMON /OPTICV/ DTAUV(NLAYER,NSPECV,4),TAUV(NLEVEL,NSPECV,4)
& ,WBARV(NLAYER,NSPECV,4), COSBV(NLAYER,NSPECV,4)
COMMON /SPECTV/ BWNV(NSPC1V),WNOV(NSPECV),DWNV(NSPECV)
& ,WLNV(NSPECV)
COMMON /FLUXV/ FNETV(NLEVEL), FUPV(NLEVEL,NSPECV),
& FDV(NLEVEL,NSPECV), FMNETV(NLEVEL),FMUPV(NLEVEL),FMDV(NLEVEL)
COMMON /FLUX/ FNET(NLEVEL), FMNET(NLEVEL), THEAT(NLAYER)
COMMON /PLANT/ CSUBP,RSFI,RSFV,F0PI
COMMON /ADJUST/ RHCH4,FH2,FHAZE,FHVIS,FHIR,TAUFAC,RCLOUD,FARGON
COMMON /CONST/RGAS,RHOP,PI,SIGMA
COMMON /IO/ TIDAL
C*
C* IPRINT CONTOLS OUTPUT AMOUNT:0=IRREDUCIBLE OUTPUT,LESS THAN 1 PAGE
C* PER RUN, 0=MINIMAL OUTPUT, 1=BACKGROUND ATM AND SPEC; 10=FULL DEBUG
WRITE(6,1111)
1111 FORMAT(' TITAN GREENHOUSE MODEL TUCSON 87 RUNS ',
&' WITH VALUES CONSTRAINED BY GEOMETRIC ALBEDO WITH GAS2 OPACITY
& ' //)
IPRINT=1
C*
C* SET UP BACKGROUNG ATMOSHPERE
C*
C* SET UP BACKGROUNG ATMOSHPERE
C FAC IS THE RATIO BETWEEN GRID SPACING TOP/BOTTOM (1 IS EVEN)
FAC=15.
ZTOP=300.
ZTOP1=250.
C MODIFY ADJUSTABLE NUMBERS HERE -- NOT IN COMMON
taufac=0.0
C fhir=0.6
c ?flag?
C print *,'enter fhaze'
C accept *,fhaze
C print *,'enter N-IR factor'
C accept *,fhir
c print *,'enter cloud radius'
c accept *,rcloud
C print *,'enter cloud opacity'
C accept *,taufac
C
CALL ZGRID(NLEVEL,FAC,ZTOP,ZTOP1,Z)
CALL ATMSETUP(NLEVEL,Z,RHCH4,FH2,FARGON,TEMP,PRESS,DEN,XMU,
& CH4,H2,XN2,AR,IPRINT)
C
C NOW CALCULATE THE LAYER AVERAGE GAS MIXING RATIOS.
CALL GASSES(IPRINT)
C
C* CALL A ROUTINE THAT SETS UP THE IR SPECTRAL INTERVALS
CALL SETSPI(IPRINT)
CALL SETSPV(IPRINT)
C SET UP PIA COEFFICIENTS
CALL SETPIA(IPRINT,1)
C
C CALL A ROTUINE THAT SETS UP THE AEROSOL DIST AND OPTICAL PROP.
IF ( FHAZE .GT. 0.) CALL HAZE2(IPRINT,NLAYER,Z,RADIUS,XNUMB)
IF (TAUFAC .GT. 0.) CALL CLD(IPRINT)
C
C* CALL A SUBROUTINE THAT SETS UP THE OPTICAL PROPERTIES IN THE
C INFRARED. AND THEN IN THE VISIBLE.
CALL OPTCI(IPRINT)
CALL OPTCV(IPRINT)
C*******************************************************************
C
C CALL A ROUTINE THAT GETS THE NET FLUXES IN THE VISIBLE
C INTEGRATING OVER ALL SPECTRAL INTERVALS AND THE IR NEXT
CALL SFLUXV(IPRINT)
CALL SFLUXI(IPRINT)
C
C
C* NOW CALCULATE THE TOTAL NET FLUX AND THE HEATING RATES
C * WE UNLY USE .5 * THE VISIBLE TO MAKE A GLOBAL AVERAGE...SEE NOTES
C ADD TIDAL HEATING IN CGS UNITS NEGATIVE MEANS ADDED HEAT FLOW.
C 1500 IS ONE IO UNIT OF HEAT FLOW. ?FLAG?
TIDAL =-1.5E3
TIDAL = 0.*TIDAL
DO 551 J=1,NLEVEL
FNET(J)=FNETI(J)+FNETV(J)*0.5 + TIDAL
FMNET(J)=FMNETI(J)+FMNETV(J)*0.5 +TIDAL
551 CONTINUE
C
CALL PRTOUT(IPRINT)
C
C BRIT2.F GOES HERE WHEN MAKING UP TGMBRIT.FOR
C.......................................................
C OUTPUT THE TEMP PROFILE
WRITE(6,891)
891 FORMAT(//' INTITAL TEMPERATURE PROFILE')
CALL FOUT(0,0)
C DISABLE CONVERGENCE, PLOT MODEL ATMOSPHERE: GO TO 9876
c go to 9876
C
C CALL A ROUTINE THAT GETS A FIRST GUESS AT A
C CONVERGED RADIATIVE CONVECTIVE TEMPERATURE PROFILE.
NSTRAT=NLEVEL-1
ITS=4
DO 1234 III=1,5
CALL TSTART(IPRINT,NSTRAT,ITS)
CALL FOUT(NSTRAT,III)
CALL ATMTEMP(IPRINT)
CALL SFLUXV(0)
1234 CONTINUE
C CALL PRTOUT(IPRINT)
C
C OUTPUT THE FIRST GUESS TEMPERATURE PROFILE
WRITE(6,899)
899 FORMAT(///' INTITIAL RADIATIVE-CONVECTIVE TEMPERATURE PROFILE')
CALL FOUT(NSTRAT,ITS)
C JUMP TO 9876 TO GET PLOTS
C THIS IS A SUITABLE JUMP POINT IF THE COMPOSITION OR THE TROPOPAUSE
C ARE NOT TO BE DETERMINED INTERATIVELY.
GO TO 9876
C
C************************************END PHASE I***********************
C FIRST DO A COMPOSITION-TEMPERATURE ADJUSTEMENT AT FIXED NSTRAT
C
IP2=-10
ITS=2
DO III=1,5
CALL TSTART(IPRINT,NSTRAT,ITS)
CALL FOUT(NSTRAT,III)
CALL ATMTEMP(IP2)
CALL SFLUXV(IP2)
ENDDO
C
C DO LOOP IS THE CONVERGENCE LOOP
C INCLUDING ADJUSTING NSTRAT
ITMAX=10
ICON=1
DO 529 ITS=1,ITMAX
C LETS CHECK FOR INSTABILITY UP TO THREE LAYERS ABOVE THE TROPOPAUSE
IF( DTDP(NSTRAT-1) .GT. CONVEQ
& .OR. DTDP(NSTRAT-2) .GT. CONVEQ
& .OR. DTDP(NSTRAT-3) .GT. CONVEQ .OR. NSTRAT .EQ. NLEVEL) THEN
C
C WE HAVE AN UNSTABLE CONDITION OR WE ARE AT THE BONDARY
IF (ICON .EQ. 1 .AND. NSTRAT .LT. NLEVEL) THEN
C IF WE REACH HERE WE DO NOT KNOW IF WE'RE HIGH ENOUGH YET
C SINCE WE HAVE YET TO HIT A STABLE LAYER SO WE KEEP RAISING
C THE TROPOPAUSE UNTIL WE HIT AN STABLE LAYER WHICH SETS ICON=0
C WE ALWAYS APPROACH THE STABILITY CRITERIA FROM ABOVE
C WE MAY BE AT THE BOUNDARY ALSO,
C EITHER WAY WE HAVE TO RAISE THE TROPOPAUSE (LOWER NSTRAT)
NSTRAT=NSTRAT-1
C CALL SUB TO SET UP NEW TEMP DEPENDENT
C
CALL ATMTEMP (IP2)
CALL SFLUXV(IP2)
CALL TSTART(IPRINT,NSTRAT,2)
CALL FOUT(NSTRAT,ITS)
CALL ATMTEMP (IP2)
CALL SFLUXV(IP2)
CALL TSTART(IPRINT,NSTRAT,2)
CALL FOUT(NSTRAT,ITS)
ELSE
C IF WE REACH HERE WE HAVE THE FIRST UNSTABLE LOOP AFTER STABLE ONES
C OR THAT THE TROPOPAUSE HAS BEEN LOWERED TO THE GROUND. EITHER WAY
C THAT MUST MEAN THAT THE CORRECT VALUE OF NSTRAT IS THE PREVIOUS ONE
C
NSTRAT=NSTRAT-1
C CALL SUB TO SET UP NEW TEMP DEPENDENT
C
CALL TSTART(IPRINT,NSTRAT,2)
CALL ATMTEMP (IP2)
CALL SFLUXV(IP2)
CALL TSTART(IPRINT,NSTRAT,2)
CALL ATMTEMP (IP2)
CALL SFLUXV(IP2)
WRITE(6,*)' BEGIN FINAL CONVERGENCE'
CALL FOUT(NSTRAT,ITS)
CALL TSTART(IPRINT,NSTRAT,5)
WRITE (6,821)
CALL ATMTEMP (IPRINT)
CALL SFLUXV(IP2)
CALL PRTOUT(IPRINT)
821 FORMAT('1 CONVERGENCED SOLUTION')
CALL FOUT(NSTRAT,ITS)
GO TO 9876
ENDIF
C
ELSE
C IF WE REACH HERE WE HAVE STABLE CONDITIONS BUT WE LOWER THE TROPOPAUSE
C TO SEE IF WE MAINTAIN STABILITY (WE NEVER REACH HERE IF NSTRAT=NLEVEL)
NSTRAT=NSTRAT+1
C CALL SUB TO SET UP NEW TEMP DEPENDENT
C
C AND THERE IS NO POINT IN DOING CALCULATIONS IF IT EQUALS NLEVEL
IF (NSTRAT .LT. NLEVEL) THEN
CALL TSTART(IPRINT,NSTRAT,2)
CALL FOUT(NSTRAT,ITS)
CALL ATMTEMP (IP2)
CALL SFLUXV(IP2)
CALL TSTART(IPRINT,NSTRAT,2)
CALL ATMTEMP (IP2)
CALL SFLUXV(IP2)
CALL TSTART(IPRINT,NSTRAT,2)
CALL FOUT(NSTRAT,ITS)
ENDIF
ICON=0
ENDIF
529 CONTINUE
C PLOT ROUTINE OUTPUTS
C JUMP TO 9876 TO GET PLOTS
9876 CONTINUE
IP1=7
WRITE(IP1,*)Z,PRESS,TEMP,WLNV,FNETV,SOLARF,
& TAUHV,FDV,FUPV,TAUGI,WNOI,FNET,FNETI
& ,(FDI(NLEVEL,K),K=1,NSPECI)
& ,(FUPI(NLEVEL,K),K=1,NSPECI)
& ,(FUPI(1,K),K=1,NSPECI)
& ,(((1.-RSFV)*FDV(NLEVEL,K)*1.E-4/(1./BWNV(K+1)-1./BWNV(K)))
& ,K=1,NSPECV)
& ,TAUHI,TAUCI,TAURV,TAUCV,TAUGV
STOP 'TGM DONE'
END
SUBROUTINE FOUT(NSTRAT,ITS)
PARAMETER (NLAYER=30,NLEVEL=31)
PARAMETER (NSPECI=46,NSPC1I=47,NSPECV=24,NSPC1V=25)
COMMON /ATM/ Z(NLEVEL),PRESS(NLEVEL),DEN(NLEVEL),TEMP(NLEVEL)
COMMON /LAPSE/ DTDP(NLAYER),CONVEQ
COMMON /FLUXI/ FNETI(NLEVEL), FUPI(NLEVEL,NSPECI),
& FDI(NLEVEL,NSPECI),FMNETI(NLEVEL),FMUPI(NLEVEL),FMDI(NLEVEL)
COMMON /FLUX/ FNET(NLEVEL), FMNET(NLEVEL), THEAT(NLAYER)
WRITE (6,814) TEMP(NLEVEL), NSTRAT, ITS
WRITE(6,892) (J,Z(J),TEMP(J),PRESS(J)
&,FMNETI(J),DTDP(J),FMNET(J),J=1,NLAYER),NLEVEL,Z(NLEVEL)
&,TEMP(NLEVEL),PRESS(NLEVEL),FNETI(NLEVEL)
&,DTDP(NLEVEL-1),FNET(NLEVEL)
814 FORMAT(////' SURFACE TEMP=', F8.3,' NSTRAT=',I3,' ITS=',I2)
892 FORMAT(/' LEVEL ALT TEMP PRESS ',
& 'IR-FNET DT/DP FNET-MID'/
& (1X,I2,F7.2,2F8.3,1PE11.3,0P,F11.5,1PE11.3,0P))
RETURN
END
SUBROUTINE MATSOL(A,NLEVEL,NSTRAT,DFLUX)
C THIS SUBROUTINE SOLVES AX=B MATRIX EQUATION
DIMENSION A(1,1),DFLUX(1)
DIMENSION IPVT(101),WKAREA (1152)
C CALL SGEFA(A,NLEVEL,NSTRAT,IPVT,INFO)
C IF (INFO .NE. 0) WRITE (6,1927)
C 1927 FORMAT(/' ****ERROR IN SGEFA IN TITAN MAIN CODE AT ITERATION')
C CALL SGESL(A,NLEVEL,NSTRAT,IPVT,DFLUX,0)
C
C IMSL ROUTINES
M=1
IDGT=3
CALL LEQT2F(A,M,NSTRAT,NLEVEL,DFLUX,IDGT,WKAREA,IER)
IF (IER .NE. 0) WRITE (6,1927)IER
1927 FORMAT(/' *****ERROR IN MATRIX SOLVER IN TITAN MAIN CODE ',I5)
RETURN
END
BLOCK DATA TGMDAT
PARAMETER (NLAYER=30)
COMMON /UBARED/ UBARI,UBARV,UBAR0
COMMON /LAPSE/ DTDP(NLAYER),CONVEQ
COMMON /PLANT/ CSUBP,RSFI,RSFV,F0PI
COMMON /ADJUST/ RHCH4,FH2,FHAZE,FHVIS,FHIR,TAUFAC,RCLOUD,FARGON
COMMON /CONST/RGAS,RHOP,PI,SIGMA
C*
DATA PI/3.14159265358979323846/
C RGAS IS THE UNIVERSAL GAS CONSTANT IN UNITS OF: M SEC-2 AMU K-1 KM
DATA RGAS/8.31432/
C SIGMA IS THE STEFAN-BOLTZMAN CONSTATNT IN CGS UNITS
DATA SIGMA/5.6677E-5/
C*
C RT CONSTANTS
DATA UBARI,UBARV,UBAR0/0.5,0.5,0.5/
C* PLANET SPECIFIC CONSTANTS
C
C CONVEQ IS THE DRY N2 ADIABATE DLNT/DLNP
DATA CONVEQ/0.2994/
C CSUBP IS THE SPECIFIC HEAT AT CONSTANT P OF THE ATMOSPHERE
C IN UNITS OF ERGS K-1 G-1
DATA CSUBP/1.05E7/
C RHOP IS THE UNITS CONVERSION FROM TO GET MASS UNITS (G CM-2)
C FROM PRESSURE (BARS) DEVIDED BY GRAVITY (M SEC-2)
C IS EQUAL TO ONE GM CM-2 BARS-1 M SEC-2
C IF ONE WHATS TO CHANGE UNITS ON PRESSURE THIS
C CONSTANT MUST BE CHANGED
DATA RHOP/1.E4/
C RSF IS THE SURFACE REFLECTANCE FOR VIS AND IR
DATA RSFV,RSFI/0.1,0.0/
C FOPI IS THE ACTUAL SOLAR FLUX IN ERGS/CM2
DATA F0PI/1.5E4/
C RHCH4 IS THE METHANE RH AT THE SURFACE
DATA RHCH4 /0.60/
C FH2 IS THE CONSTANT MIXING RATIO OF H2
DATA FH2/0.003/
C FHAZE IS THE HAZE PRODUCTION SCALING FACTOR
DATA FHAZE/0.35/
C FHIR IS THE HAZE INFRARED ABSORPTION SCALE FACTOR
DATA FHIR/0.5/
C FHVIS IS THE HAZE INFRARED ABSORPTION SCALE FACTOR
DATA FHVIS/1.333333333/
C TAUFAC IS THE 200 CM-1 SCALING FACTOR
DATA TAUFAC/2.00/
C RCLOUD IS THE PARTICLE SIZE IN THE CLOUD IN MICRONS
DATA RCLOUD/100./
END
SUBROUTINE XMIE(RAD,RREAL,XIMG,QEXT,QSCT,QABS,CBAR,WAVEN)
C CHANGED TO WORK WITH THE VAX MIE VERSION
COMPLEX*16 ACAP(7000)
REAL*8 ELTRMX(4,1,2),PIE(3,1),TAU(3,1),CSTHT(1),SI2THT(1)
REAL*8 RO,REALO,XIMGO,QE,QS,CB,RD2,RN,CN,TPILAM
RO=RAD
REALO=RREAL
XIMGO=XIMG
RD2=0.
RN=1.0
CN=0.
WAVEL=1.E4/WAVEN
TPILAM=2.0*3.14159/WAVEL
CALL DMIESS (RO, REALO, XIMGO, 0.0, 1,
* QE, QS,CB,
* ELTRMX,PIE,TAU,CSTHT,SI2THT,ACAP,1,
* 7000, RD2, RN, CN, TPILAM)
CBAR=CB/QS
QEXT=QE*3.14159*RAD*RAD*1.E-8
QSCT=QS*3.14159*RAD*RAD*1.E-8
QABS=QEXT-QSCT
SIZE=3.14159*RAD*RAD*1.E-8
RETURN
END
SUBROUTINE TSTART(IPRINT,NSTRAT,ITMAX)
PARAMETER (NLAYER=30,NLEVEL=31)
PARAMETER (NSPECI=46,NSPC1I=47,NSPECV=24,NSPC1V=25)
COMMON /ATM/ Z(NLEVEL),PRESS(NLEVEL),DEN(NLEVEL),TEMP(NLEVEL)
COMMON /LAPSE/ DTDP(NLAYER),CONVEQ
COMMON /FLUXI/ FNETI(NLEVEL), FUPI(NLEVEL,NSPECI),
& FDI(NLEVEL,NSPECI),FMNETI(NLEVEL),FMUPI(NLEVEL),FMDI(NLEVEL)
COMMON /FLUXV/ FNETV(NLEVEL), FUPV(NLEVEL,NSPECV),
& FDV(NLEVEL,NSPECV), FMNETV(NLEVEL),FMUPV(NLEVEL),FMDV(NLEVEL)
COMMON /FLUX/ FNET(NLEVEL), FMNET(NLEVEL), THEAT(NLAYER)
DIMENSION BETA(NLEVEL), DFLUX(NLEVEL), A(NLEVEL,NLEVEL)
& ,FMIP(NLEVEL), FNETIP(NLEVEL)
COMMON /IO/ TIDAL
DATA DELT /10./
C NSTRAT IS DEFINED TO BE THE LEVEL NUMBER WHICH
C CORRESPONDS TO THE TROPOPAUSE. LAYERS WITH NUMBERS
C LESS THAN NSTRAT ARE IN RADIATIVE EQUILIBRIUM. THERE
C ARE NSTRAT-1 OF THESE LAYERS. LEVELS AND MIDPOINTS
C WITH NUMBERS LESS THAN NSTRAT SHOULD HAVE ZERO FLUX.
C LEVEL NSTRAT ITSELF IS NOT CONTRAINED TO HAVE ZERO FLUX.
C
C IN THIS SUBROUTINE NSTRAT IS DEFINED AND CONSTANT
C
DO 9999 ITS = 1, ITMAX
C
C SET UP THE INDEPENDENT VARIABLES FOR THE ITERATION PROCESESS
C THESE ARE TEMPERATURES AT LEVELS UP TO NSTRAT, THESE ARE PUT
C SEQUENTIALLY IN BETA VECTOR ALSO SET UP THE DFLUX VECTOR
C INTITAL SET UP OF THESE VECTORS
DO 1025 J=1,NSTRAT
BETA(J)= TEMP(J)
1025 CONTINUE
C DFLUX VECTOR CONSISTES OF FLUX AT TOP AND MIDPOINT FLUXES
C FOR THE NSTRAT-1 LAYERS IN THE STRATOSPHERE
C
DFLUX(1) = FNET(1)
DFLUX(NSTRAT)=FNET(NSTRAT)
DO 1026 J=2,NSTRAT
DFLUX(J) = FMNET(J-1)
1026 CONTINUE
C
C
C BEGIN JACOBIAN LOOP****************************************
C AT THIS POINT WE HAVE VALUES FOR THE NET FLUXES AND THE VECTOR
C DFLUX IS SET UP AS THE APPROPRIATE NET FLUXES. WE NOW WISH TO CALULATE
C THE JACOBIAN
C
C WE MUST STORE THE VALUE OF THE FLUX CALCULATED WITHOUT PERTURBATION
C IN FNETP AND FMNETP
DO 219 J=1,NLEVEL
FNETIP(J)=FNETI(J)
FMIP (J)=FMNETI(J)
219 CONTINUE
C
C BEGIN PERTURBATION CALCULATION
C
DO 1020 JM = 1, NSTRAT
C CHOSE A UNIT PERTURBATION IT IS FIVE DEG
BETA(JM)=BETA(JM) + DELT
C NOW RECONSTRUCT THE TEMPERATURE PROFILE
DO 1019 J1 = 1,NSTRAT
TEMP(J1) = BETA(J1)
1019 CONTINUE
DO 1582 J1 = NSTRAT+1, NLEVEL
TEMP(J1)=EXP(LOG(TEMP(J1-1))+CONVEQ*
& (ALOG(PRESS(J1))-ALOG(PRESS(J1-1))))
1582 CONTINUE
C
C NOW RECALULATE THE IR FLUXES GIVEN THE PURTURBATION IN TEMP
C
CALL SFLUXI(IPRINT)
C
C******************************************************************
C
C NOW WE CALCULATE THE JACOBIAN TERMS IN THE SAME MANNER
C AS THE DFLUX TERMS -- THESE ARE LINKED!!! --
C
A(1,JM)=(FNETI(1) - FNETIP(1))/DELT
A(NSTRAT,JM) = (FNETI(NSTRAT) - FNETIP(NSTRAT) )/DELT
DO 1018 IM=1,NSTRAT-1
A(IM+1,JM) = ( FMNETI(IM) - FMIP(IM) )/DELT
1018 CONTINUE
C******************************************************************
C UNDO BETA VECTOR PERTURBATION
BETA(JM) = BETA(JM) - DELT
C***********************************************************************
C234567..............................................................012
1020 CONTINUE
C
C PRINT OUT JACOBIAN MATRIX
IF (ITS .EQ. 1 .AND. IPRINT .GT. 9) THEN
WRITE (6,1400) 1,FNETIP(1),FNETI(1), (J+1,FMIP(J),FMNETI(J)
& ,J=1,NSTRAT)
1400 FORMAT (//' J AFTER BEFORE'/(0PI3,1P2E12.3))
DO 1225 IM = 1, NSTRAT
WRITE (6,1411) IM
1411 FORMAT(' JACOBIAN ROW: ',I3)
WRITE (6,1311) ( A(IM,JM),JM=1,NSTRAT)
1311 FORMAT (1X, 1P8E10.2)
1225 CONTINUE
END IF
C
C NOW CALL MATRIX SOLVER AND GET THE CHANGE IN THE
C BETA VECTOR AND USE TO COMPUTE A NEW BETA VECTOR
CALL MATSOL(A,NLEVEL,NSTRAT,DFLUX)
C
C **********************************************
DO 1239 J=1,NSTRAT
C THIS IS AN ARTIFICAL DAMPER .....
IF ( ABS(DFLUX(J)) .GT. 10.0)
& DFLUX(J)=10.*DFLUX(J)/ABS(DFLUX(J))
C
BETA(J)=BETA(J)-DFLUX(J)*.5
1239 CONTINUE
C
C WE NOW HAVE A NEW BETA VECTOR - AN ITERATION HAS BEEN COMPLETED
C
C NOW RECOMPUTE THE TEMPERATURE PROFILE
DO 1119 J1 = 1,NSTRAT
TEMP(J1) = BETA(J1)
1119 CONTINUE
DO 1120 J1 = NSTRAT+1,NLEVEL
TEMP(J1)=EXP(LOG(TEMP(J1-1))+CONVEQ*
& (ALOG(PRESS(J1))-ALOG(PRESS(J1-1))))
1120 CONTINUE
C THIS IS AN ARTIFICAL DAMPER ....
DO 1121 J1=1,NLEVEL
IF (TEMP(J1) .LT. 10.) TEMP(J1)=50.00
1121 IF (TEMP(J1) .GT. 300.) TEMP(J1)=300.00
C
C RECALCULATE THE NET IR FLUXES
CALL SFLUXI(IPRINT)
C
C* NOW RECALCULATE THE TOTAL NET FLUX
DO 553 J=1,NLEVEL
FNET(J)=FNETI(J)+FNETV(J)*0.5 +TIDAL
FMNET(J)=FMNETI(J)+FMNETV(J)*0.5 + TIDAL
553 CONTINUE
C
9999 CONTINUE
C
C END OF ITERATION LOOP ***********************************
C
C NOW RECALCULATE THE DFLUX VECTOR
DFLUX(1) = FNET(1)
DFLUX(NSTRAT) = FNET(NSTRAT)
DO 4025 J=2,NSTRAT
DFLUX(J) = FMNET(J-1)
4025 CONTINUE
C NOW CALCULATE THE ACTUAL LAPSE RATES
DO 529 J=1,NLAYER
DTDP(J)=(LOG( TEMP(J)) - LOG( TEMP(J+1)))
& /(LOG(PRESS(J)) - LOG(PRESS(J+1)))
529 CONTINUE
RETURN
END
SUBROUTINE SETSPI(IPRINT)
PARAMETER (NSPECI=46,NSPC1I=47)
COMMON /SPECTI/ BWNI(NSPC1I),WNOI(NSPECI),DWNI(NSPECI)
& ,WLNI(NSPECI)
C SET UP SPECRAL INTERVALS
BWNI(1)=8.
BWNI(2)=15.
BWNI(3)=25.
BWNI(4)=37.5
DO 100 K=5,27
BWNI(K)=BWNI(K-1)+25.
100 CONTINUE
BWNI(28)=645.
BWNI(29)=645.+35.6760
BWNI(30)=BWNI(29)+4.324/2.
BWNI(31)=BWNI(30)+4.324/2.
BWNI(32)=BWNI(31)+12.327/4.
BWNI(33)=BWNI(32)+12.327/4.
BWNI(34)=BWNI(33)+12.327/4.
BWNI(35)=BWNI(34)+12.327/4.
BWNI(36)=BWNI(35)+35.6290
BWNI(37)=BWNI(36)+7.044/4.
BWNI(38)=BWNI(37)+7.044/4.
BWNI(39)=BWNI(38)+7.044/4.
BWNI(40)=BWNI(39)+7.044/4.
BWNI(41)=BWNI(40)+16.32/3.
BWNI(42)=BWNI(41)+16.32/3.
BWNI(43)=BWNI(42)+16.32/3.
BWNI(44)=BWNI(43)+156.48
BWNI(45)=BWNI(44)+27.2/3.
BWNI(46)=BWNI(45)+27.2/3.
BWNI(47)=BWNI(46)+27.2/3.
C
C SET UP MEAN WAVENUMBERS AND DELTAS
C UNITS ON WAVENUMBER ARE CM-1
C UNITS ON WAVELN ARE MICRONS
DO 160 K=1,NSPECI
WNOI(K)=( BWNI(K+1)+BWNI(K) )*0.5
DWNI(K)=BWNI(K+1)-BWNI(K)
WLNI(K)=1.E+4/WNOI(K)
160 CONTINUE
C IF THERE IS ONLY ONE INTERVAL THEN TOTAL ENERGY IS USED AND
IF (NSPECI .EQ. 1) DWNI(1) = 1.0
C PRINT OUT SPECTRAL INTERVALS
IF (IPRINT .GT. 1) THEN
WRITE (6,190)
DO 200 K=1,NSPECI
WRITE (6,210)K,WLNI(K),WNOI(K),BWNI(K)
& ,BWNI(K+1),DWNI(K)
200 CONTINUE
END IF
210 FORMAT(1X,I3,F10.3,F10.2,F10.2,'-',F8.2,F10.3)
190 FORMAT(///' SPECTRAL INTERVALS'//
& ' SNUM MICRONS WAVENU INTERVAL',11X,'DELTA-WN')
C****** END SPECTRAL INTERVAL SET UP *************
RETURN
END
SUBROUTINE AEROSL(IPRINT)
C THIS ROUTINE SETS UP THE AEROSOL DISTRIBUTION
C
PARAMETER (NLAYER=30,NLEVEL=31)
PARAMETER (NSPECI=46,NSPECV=24)
COMMON /ATM/ Z(NLEVEL),PRESS(NLEVEL),DEN(NLEVEL),TEMP(NLEVEL)
COMMON /GASS/ CH4(NLEVEL),XN2(NLEVEL),H2(NLEVEL),AR(NLEVEL)
& ,XMU(NLEVEL),GAS1(NLAYER),COLDEN(NLAYER)
COMMON /AERSOL/ RADIUS(NLAYER), XNUMB(NLAYER)
& , REALI(NSPECI), XIMGI(NSPECI), REALV(NSPECV), XIMGV(NSPECV)
COMMON /PLANT/ CSUBP,RSFI,RSFV,F0PI
COMMON /ADJUST/ RHCH4,FH2,FHAZE,FHVIS,FHIR,TAUFAC,RCLOUD,FARGON
COMMON /CONST/RGAS,RHOP,PI,SIGMA
DIMENSION PROD(NLEVEL)
C
C LETS REDEFINE THE AEROSOLS BASED ON THE MCKAY SIMPLE
C AEROSOL MODEL... SEE NOTES
C PRODUCTION PROFILE FROM TOON ET AL. (1980)
XKB=1.380E-16
RHOH=1.
PX=1.E-7
RCRIT=0.09E-4
c
PY=PX*EXP(-1.124)
C
C
DO 500 J=1,NLEVEL
500 PROD(J)=FHAZE*2.E-7*3.5E-14*EXP(-0.5*( (PRESS(J)-PX)/PY)**2)
C
C DUE TO NOLINEAR AVERAGEING WE NEED TO DETERMINE THE TOTAL
C COLUMN PRODUCTION NUMERICALLY, AND ITERATIVELY
DO ITS=1,4
CLEVEL=0.0
DO J=2,NLEVEL
PRODMID=AVERGE( PROD(J),PROD(J-1) )
CLEVEL=CLEVEL+PRODMID*(Z(J-1)-Z(J))*1E5
ENDDO
FACTOR=FHAZE*3.5E-14/CLEVEL
C
C SCALE THE RESULTS TO ADJUST THE PRODUCTION RATE
DO J=1,NLEVEL
PROD(J)=PROD(J)*FACTOR
ENDDO
CC
ENDDO
C
TAU=0.
IA=6
IFLAG=0
IF (IPRINT .GT. 0) WRITE(IA,389)
389 FORMAT('1 LVL ALT PRESS TAU(.635) PROD CUM-PROD'
&,' DEN FALL-V R (CM)')
C.........
DO 501 J=1,NLEVEL
GASDEN=1.E6*PRESS(J)/(XKB*TEMP(J))
EG = 100.*EFFG(Z(J))
XMASS=XMU(J)/6.022E23
B=EG*RHOH*0.5*SQRT(PI/(2.*XMASS*XKB*TEMP(J)))/GASDEN
COAG=4.*SQRT(3.*XKB*TEMP(J)/RHOH)
IF (J .GT. 1 ) THEN
PRODMID=AVERGE( PROD(J),PROD(J-1) )
CMID=CLEVEL+0.5*PRODMID*(Z(J-1)-Z(J))*1E5
CLEVEL=CLEVEL+PRODMID*(Z(J-1)-Z(J))*1E5
C DO THE STOKES LIMIT ADJUSTMENT
C AS DISCRIBED IN WRITE UP.
PATH=1.689E14/GASDEN
XKNUD=PATH/RAD2
ETA=1.78E-4*SQRT(TEMP(J)/300.)
IK=1
IF (XKNUD .LT. 100.) THEN
B= (B*XKNUD**IK + 5.4*2.*RHOH*EG*RAD2/(9.*ETA) )/(XKNUD**IK +1.)
COAG=(COAG*XKNUD**IK + 4.*XKB*TEMP(J)/(3.*ETA*SQRT(RAD2)))
& /(XKNUD**IK +1.)
ENDIF
C END OF STOKES FLOW ADJUSTMENT
C
C PUT IN CHARGE REDUCTION IN COAGULATION ?FLAG?
COAG=COAG*EXP(-RAD2/RCRIT)
C END OF CHARGE ADJUST
R=( R**4.5 + ( 9.*COAG*CMID/(8.*PI*RHOH*B*B)
& + 9.*SQRT(R)*PRODMID/(8.*PI*RHOH*B*DEN1)) *
& (Z(J-1)-Z(J))*1.E5 )**(2./9.)
DEN1=CLEVEL*3./(4.*PI*RHOH*B*R**4)
C USE AVERAGEING PROCESS DENSITY AS DETERMINED IN NOTES
C TO GET VALUES FOR THE LAYERS
XNUMB(J-1)=AVERGE(DEN1,DEN2)*1.E5*( Z(J-1) - Z(J) )
RADIUS(J-1)=1.E4 * AVERGE(R,RAD2)
ELSE
C TOP LAYER CALCULATIONS (J .EQ. 1)
CLEVEL=0.0
R= ( 9.*COAG*CLEVEL/(8.*PI*RHOH*B*B) )**(2./9.)
CC BRITNESS SCALE FACTOR FOR TOP LAYER ?FLAG??
R=2.E-6
DEN1=4.*FHAZE/20.
ENDIF
RAD2=R
DEN2=DEN1
C CALL THE MIE CODE TO GIVE THE AEROSOL PROPERTIES AT A REF WAVELN
WVLN=.635
RMICRON=R*1E4
WNO=1.E4/WVLN
CALL THOLIN(WVLN,REAL,XIMG)
XIMG=XIMG*FHVIS
CALL XMIE(RMICRON,REAL,XIMG,
& QEXT,QSCT,QABS,CBAR,WNO)
DTAU=QEXT*XNUMB(J-1)
TAU=TAU+DTAU
C
C LETS DROP OUT THE HAZE FOR ALL ALTITUDES BELOW THE CONDENSATION POINT
C WITH A ADD-ON HERE. ?FLAG??
IF ( IFLAG .EQ. 1 ) XNUMB(J-1)=0.0
IF ( CH4(J)*PRESS(J)/PCH4(TEMP(J)) .GT. 0.9 ) IFLAG=1
C
IF (IPRINT .GT. 0)
& WRITE(IA,390) J,Z(J),PRESS(J),TAU,PROD(J),CLEVEL,DEN1,B,R
501 CONTINUE
390 FORMAT(I3,F6.0,1P7E10.2)
RETURN
END
SUBROUTINE CLD(IPRINT)
C PUT IN A METHANE CLOUD HERE
C THIS ROUTINE SETS UP THE CLOUD DISTRIBUTION
C
PARAMETER (NLAYER=30,NLEVEL=31)
PARAMETER (NSPECI=46,NSPC1I=47,NSPECV=24,NSPC1V=25)
COMMON /ATM/ Z(NLEVEL),PRESS(NLEVEL),DEN(NLEVEL),TEMP(NLEVEL)
COMMON /GASS/ CH4(NLEVEL),XN2(NLEVEL),H2(NLEVEL),AR(NLEVEL)
& ,XMU(NLEVEL),GAS1(NLAYER),COLDEN(NLAYER)
COMMON /CLOUD/ RADCLD(NLAYER), XNCLD(NLAYER)
& , RCLDI(NSPECI), XICLDI(NSPECI), RCLDV(NSPECV), XICLDV(NSPECV)
COMMON /PLANT/ CSUBP,RSFI,RSFV,F0PI
COMMON /ADJUST/ RHCH4,FH2,FHAZE,FHVIS,FHIR,TAUFAC,RCLOUD,FARGON
COMMON /CONST/RGAS,RHOP,PI,SIGMA
TOTALC=0.0
CCC
XC=.95
DO 190 J=1,NLAYER
XNCLD(J)=0.
RADCLD(J)=0.
IF ( CH4(J)*PRESS(J)/PCH4(TEMP(J)) .GT. XC) THEN
RADCLD(J)=RCLOUD
C TO COLAPSE THE CLOUD INTO ONE LAYER: ?FLAG? XC=9.
C LET 1% OF THE GAS BE CLOUD AS AN INTITIAL GUESS
XNCLD(J)=.01*COLDEN(J)*GAS1(J)/((4.*PI/3.)*RADCLD(J)**3*1.E-12)
IF (IPRINT .GT. 0 ) WRITE(6,95) J,RADCLD(J),XNCLD(J),Z(J)
95 FORMAT(' CLOUD INSERTED: ',I3,F8.2,1P5E10.3)
TOTALC=TOTALC+XNCLD(J)
ENDIF
190 CONTINUE
C CALL THE MIE CODE TO GIVE THE AEROSOL PROPERTIES AT A REF WAVENO
C WHICH IS THE REF WAVENO OF TOON ET AL.
WNOREF=200.
RREF=1.27
XIREF=REFLIQ(WNOREF)
CALL XMIE(RCLOUD,RREF,XIREF,
& QEXT,QSCT,QABS,CBAR,WNOREF)
CTAU=QEXT*TOTALC
IF (IPRINT .GT. 0) WRITE(6,98) WNOREF,RREF,XIREF,TOTALC,CTAU
98 FORMAT(' CLOUD AT REFERENCE WAVENUMBER OF ',F7.2,' REAL, IMG =',
& 1P2E10.2/' COLUMN DENSITY , OPTICAL DEPTH= ' 2E10.2)
C SCALE THE CLOUD DENSITIES TO THE REFERENCE WAVENUMBER
DO 145 J=1,NLAYER
XNCLD(J)=XNCLD(J)*TAUFAC/CTAU
145 CONTINUE
RETURN
END
SUBROUTINE OPTCI(IPRINT)
PARAMETER (NLAYER=30,NLEVEL=31)
PARAMETER (NSPECI=46,NSPC1I=47,NSPECV=24,NSPC1V=25)
COMMON /ATM/ Z(NLEVEL),PRESS(NLEVEL),DEN(NLEVEL),TEMP(NLEVEL)
COMMON /GASS/ CH4(NLEVEL),XN2(NLEVEL),H2(NLEVEL),AR(NLEVEL)
& ,XMU(NLEVEL),GAS1(NLAYER),COLDEN(NLAYER)
COMMON /STRATO/ C2H2(NLAYER),C2H6(NLAYER)
COMMON /AERSOL/ RADIUS(NLAYER), XNUMB(NLAYER)
& , REALI(NSPECI), XIMGI(NSPECI), REALV(NSPECV), XIMGV(NSPECV)
COMMON /CLOUD/ RADCLD(NLAYER), XNCLD(NLAYER)
& , RCLDI(NSPECI), XICLDI(NSPECI), RCLDV(NSPECV), XICLDV(NSPECV)
COMMON /TAUS/ TAUHI(NSPECI),TAUCI(NSPECI),TAUGI(NSPECI),
& TAURV(NSPECV),TAUHV(NSPECV),TAUCV(NSPECV),TAUGV(NSPECV)
COMMON /OPTICI/ DTAUI(NLAYER,NSPECI),TAUI(NLEVEL,NSPECI)
& ,WBARI(NLAYER,NSPECI), COSBI(NLAYER,NSPECI)
COMMON /SPECTI/ BWNI(NSPC1I),WNOI(NSPECI),DWNI(NSPECI)
& ,WLNI(NSPECI)
COMMON /PLANT/ CSUBP,RSFI,RSFV,F0PI
COMMON /ADJUST/ RHCH4,FH2,FHAZE,FHVIS,FHIR,TAUFAC,RCLOUD,FARGON
COMMON /CONST/RGAS,RHOP,PI,SIGMA
DIMENSION PROD(NLEVEL)
C THE PRESSURE INDUCED TRANSITIONS ARE FROM REGIS
C THE LAST SEVENTEEN INTERVALS ARE THE BANDS FROM GNF.
C*
C THIS SUBROUTINE SETS THE OPTICAL CONSTANTS IN THE INFRARED
C IT CALCUALTES FOR EACH LAYER, FOR EACH SPECRAL INTERVAL IN THE IR
C LAYER: WBAR, DTAU, COSBAR
C LEVEL: TAU
C
DO 420 K=1,NSPECI
C LETS USE THE THOLIN OPTICAL CONSTANTS FOR THE HAZE.
CALL THOLIN(WLNI(K),TNR,TNI)
REALI(K)=TNR
XIMGI(K)=TNI*FHIR
C SET UP THE OPTICAL CONSTANTS FOR THE CLOUD
RCLDI(K)=1.27
XICLDI(K)=REFLIQ(WNOI(K))
420 CONTINUE
C
C ZERO ALL OPTICAL DEPTHS.
C ?FLAG? FOR SOME APPLCIATIONS THE TOP OPACITY MAY NOT VANISH
DO 80 K=1,NSPECI
TAUHI(K)=0.
TAUCI(K)=0.
TAUGI(K)=0.
80 CONTINUE
C
DO 100 J=1,NLAYER
C SET UP THE COEFFICIENT TO REDUCE MASS PATH TO STP ...SEE NOTES
C T0 =273.15 PO=1.01325 BAR
TBAR=0.5*(TEMP(J)+TEMP(J+1))
PBAR=SQRT(PRESS(J)*PRESS(J+1))
BMU=0.5*(XMU(J+1)+XMU(J))
COEF1=RGAS*273.15**2*.5E5* (PRESS(J+1)**2 - PRESS(J)**2)
& /(1.01325**2 *EFFG(Z(J))*TBAR*BMU)
IF (IPRINT .GT. 9) WRITE(6,21) J,EFFG(Z(J)),TBAR,BMU,COEF1
21 FORMAT(' J, EFFG, TBAR, BMU, COEF1,: 'I3,1P6E10.3)
DO 100 K=1,NSPECI
C1
C FIRST COMPUTE TAU AEROSOL
CALL XMIE(RADIUS(J),REALI(K),XIMGI(K),
& QEXT,QSCT,QABS,CBAR,WNOI(K))
TAEROS=QEXT*XNUMB(J)
C2
C NEXT COMPUTE TAU CLOUD
TAUCLD=0.0
IF ( XNCLD(J) .GT. 0.) THEN
CALL XMIE(RADCLD(J),RCLDI(K),XICLDI(K),
& QEXTC,QSCTC,QABSC,CBARC,WNOI(K))
TAUCLD=QEXTC*XNCLD(J)
ENDIF
C3
C NOW COMPUTE TAUGAS DUE TO THE PIA TERM ONLY FOR LAMDA LT 940
TAUGAS=0.0
IF (WNOI(K) .LT. 940. ) THEN
CALL PIA(K,TBAR,PNN,PCC,PCN,PHN)
C HERE IS WHERE WE COULD SCALE THE PIA COEFFICEINTS
C BASED ON REGIS' NOTES. ---TGM HAS NO ADJUST IN IT AS DEFAULT
C 1.25 FACTOR SUGGESTED FOR HERE BY TOON
C PCN=PCN*1.25*MIN(1.75 , MAX(1.0,WNOI(K)/200.))
C PCN=PCN*MIN(1.75 , MAX(1.0,WNOI(K)/200.))
TAUGAS=COEF1*
& (XN2(J)*XN2(J)*PNN + CH4(J)*CH4(J)*PCC
& + XN2(J)*CH4(J)*PCN + XN2(J)*H2(J)*PHN)
IF (J .EQ. NLAYER .AND. IPRINT .GT. 9)
& WRITE (6,22) WNOI(K),TAUGAS,XN2(J),CH4(J),H2(J),
& TBAR, PNN,PCC,PCN, PHN,
& XN2(J)*XN2(J)*PNN , CH4(J)*CH4(J)*PCC ,
& XN2(J)*CH4(J)*PCN , XN2(J)*H2(J)*PHN
22 FORMAT(1X,1P8E10.2)
ENDIF
IF (K .GT. 28) THEN
KGAS=K-28
C ?FLAG? HERE MUST BE WATCHED CAREFULLY
U=COLDEN(J)*6.02204E23/BMU
CALL GAS2(J, KGAS,TBAR,PBAR,U,TAU2)
TAUGAS=TAUGAS+TAU2
ENDIF
C
IF (TAEROS + TAUCLD .GT. 0.) THEN
COSBI(J,K)=(CBAR*TAEROS + CBARC*TAUCLD )/(TAEROS + TAUCLD)
ELSE
COSBI(J,K)=0.0
ENDIF
DTAUI(J,K)=TAUGAS+TAEROS+TAUCLD
TAUHI(K)=TAUHI(K) + TAEROS
TAUGI(K)=TAUGI(K) + TAUGAS
TAUCI(K)=TAUCI(K) + TAUCLD
C ?FLAG? SERIOUS PROBLEM WITH THE CODE HERE!
TLIMIT=1.E-6
IF (DTAUI(J,K) .GT. TLIMIT) THEN
WBARI(J,K)=(QSCT*XNUMB(J) + QSCTC*XNCLD(J))
& /DTAUI(J,K)
ELSE
WBARI(J,K)=0.0
WRITE (6,999) J,K,DTAUI(J,K)
999 FORMAT (' WARNING TAU MINIMUM AT J,K,DTAU:',2I3,1PE10.3)
DTAUI(J,K)=TLIMIT
ENDIF
IF (IPRINT .GT. 9)
& WRITE(6,73)J,K,TAUGAS,TAEROS,QEXT,QSCT
73 FORMAT(2I3,1P8E10.3)
100 CONTINUE
C TOTAL EXTINCTION OPTICAL DEPTHS
DO 119 K=1,NSPECI
TAUI(1,K)=0.0
DO 119 J=1,NLAYER
TAUI(J+1,K)=TAUI(J,K)+DTAUI(J,K)
119 CONTINUE
IF (IPRINT .GT. 2) THEN
WRITE (6,120)
120 FORMAT(///' OPTICAL CONSTANTS IN THE INFRARED')
DO 200 K=1,NSPECI
WRITE (6,190)
WRITE (6,210)K,WLNI(K),WNOI(K),BWNI(K)
& ,BWNI(K)+DWNI(K),DWNI(K)
WRITE (6,230)REALI(K),XIMGI(K)
DO 195 J=1,NLAYER
WRITE (6,220)XNUMB(J), WBARI(J,K),COSBI(J,K),DTAUI(J,K)
& ,TAUI(J,K)
195 CONTINUE
WRITE (6,240) TAUI(NLEVEL,K)
200 CONTINUE
END IF
210 FORMAT(1X,I3,F10.3,F10.2,F10.2,'-',F8.2,F10.3)
190 FORMAT(1X//' SNUM MICRONS WAVENU INTERVAL DELTA-WN')
230 FORMAT(1X,'NREAL(LAYER)= ',1PE10.3,' NIMG(LAYER)= ',E10.3/
&' #AEROSOLS WBAR COSBAR DTAU TAU')
220 FORMAT(5(1X,G9.3))
240 FORMAT(41X,G9.3)
RETURN
END
SUBROUTINE SETSPV(IPRINT)
PARAMETER (NSPECV=24,NSPC1V=25)
COMMON /SPECTV/ BWNV(NSPC1V),WNOV(NSPECV),DWNV(NSPECV)
& ,WLNV(NSPECV)
COMMON /VISGAS/SOLARF(NSPECV),NTERM(NSPECV),PEXPON(NSPECV),
& ATERM(4,NSPECV),BTERM(4,NSPECV)
DATA WLNV/
& 0.325, 0.375, 0.425, 0.475, 0.525, 0.575,
& 0.640, 0.715, 0.789, 0.850, 0.891, 0.935,
& 0.998, 1.073, 1.144, 1.213, 1.292, 1.381,
& 1.484, 1.603, 1.742, 1.909, 2.111, 2.361/
data solarf/
& 525.430, 680.440, 1051.900, 1173.300, 1082.000, 1056.500,
& 1495.100, 1123.100, 1053.200, 541.820, 397.630, 484.210,
& 624.540, 528.390, 383.250, 372.260, 360.090, 339.230,
& 315.600, 292.580, 267.600, 239.700, 208.040, 180.550/
data NTERM/1,1,4,4,4,3,4,4,3,4,3,3,3,4,4,3,3,4,3,4,4,3,4,4/
data pexpon/6*0.0,
& 0.000, 0.000, 0.000, 0.149, 0.156, 0.186,
& 0.302, 0.097, 1.150, 1.040, 1.030, 1.040,
& 1.080, 1.070, 1.090, 1.050, 1.050, 0.959/
data ATERM/
& 1.000000, 0.000000, 0.000000, 0.000000, 1.000000, 0.000000,
& 0.000000, 0.000000, 0.700000, 0.093900, 0.015000, 0.191100,
& 0.300030, 0.135014, 0.460546, 0.104410, 0.164416, 0.375038,
& 0.295630, 0.164917, 0.444755, 0.355864, 0.199380, 0.000000,
& 0.198120, 0.331633, 0.335834, 0.134413, 0.327300, 0.335500,
& 0.146800, 0.190400, 0.277055, 0.333367, 0.389578, 0.000000,
& 0.126666, 0.416016, 0.364590, 0.092728, 0.286216, 0.492476,
& 0.221308, 0.000000, 0.445402, 0.453717, 0.100882, 0.000000,
& 0.351017, 0.371854, 0.277129, 0.000000, 0.434400, 0.477200,
& 0.077600, 0.010800, 0.292343, 0.374023, 0.242032, 0.091602,
& 0.501604, 0.385625, 0.112771, 0.000000, 0.597075, 0.308155,
& 0.094771, 0.000000, 0.116916, 0.447207, 0.338013, 0.097864,
& 0.541475, 0.367862, 0.090663, 0.000000, 0.468164, 0.212875,
& 0.213276, 0.105685, 0.245440, 0.416617, 0.242734, 0.095209,
& 0.330423, 0.503512, 0.166065, 0.000000, 0.307538, 0.361097,
& 0.232456, 0.098909, 0.115609, 0.353355, 0.413319, 0.117718/
data BTERM/
&0.0000E+00,0.0000E+00,0.0000E+00,0.0000E+00,0.0000E+00,0.0000E+00,
&0.0000E+00,0.0000E+00,0.0000E+00,7.8260E-04,2.3210E-03,5.5120E-03,
&0.0000E+00,9.4390E-04,4.5750E-03,4.2690E-02,1.9500E-04,2.5340E-03,
&1.5020E-02,8.9550E-02,2.5430E-03,1.5420E-02,3.1900E-02,0.0000E+00,
&7.1740E-03,2.6010E-02,9.7160E-02,4.4030E-01,4.4550E-02,1.8750E-01,
&7.9000E-01,2.8660E+00,5.6920E-02,1.9440E-01,7.9630E-01,0.0000E+00,
&3.1190E-02,4.9600E-01,2.2970E+00,2.8660E+01,1.2740E+00,1.4250E+01,
&1.4070E+02,0.0000E+00,9.4170E-02,6.1850E-01,6.9190E+00,0.0000E+00,
&1.8690E+00,1.2560E+01,1.4730E+02,0.0000E+00,1.3600E-01,9.3830E-01,
&4.2070E+00,1.4080E+02,6.1300E-02,1.3220E+00,2.2710E+01,8.5640E+02,
&0.0000E+00,6.0630E-01,3.9880E+01,0.0000E+00,0.0000E+00,6.0850E-01,
&4.5140E+01,0.0000E+00,2.7190E-01,1.9420E+00,3.1990E+01,1.5370E+03,
&0.0000E+00,5.8440E-01,5.0650E+01,0.0000E+00,6.2920E-03,7.8350E-01,
&3.5630E+01,1.2870E+03,1.7680E+00,3.0850E+01,4.6500E+02,1.3940E+04,
&0.0000E+00,6.6970E-01,4.6030E+01,0.0000E+00,0.0000E+00,7.3460E-01,
&3.9680E+01,3.1920E+03,9.5930E+00,3.7750E+02,5.1130E+03,2.7940E+05/
C*******
C SET UP SPECRAL INTERVALS THESE ARE BASED ON KATHY RAGES PROGRAM
C CONVERT PER KM-AMAGATS TO PER GM CM-3 (SEE NOTES FOR CONSTANT)
DO K=1,NSPECV
DO NT=1,4
BTERM(NT,K)=BTERM(NT,K)/(1.E5 * 16./22.4E3)
ENDDO
ENDDO
C
C SET UP MEAN WAVENUMBERS AND DELTAS
C UNITS ON WAVENUMBER ARE CM-1
C UNITS ON WAVELN ARE MICRONS
BWNV(1)=1.E4/.3
DO 100 K=2,NSPC1V
BWLN=1.E4/BWNV(K-1)
EWLN=2.*WLNV(K-1)-BWLN
BWNV(K)=1.E4/EWLN
100 CONTINUE
DO 160 K=1,NSPECV
WNOV(K)=1.E4/WLNV(K)
DWNV(K)=BWNV(K)-BWNV(K+1)
160 CONTINUE
C IF THERE IS ONLY ONE SPECTRAL INTERVAL THEN TOTAL E IS USED AND
IF (NSPECV .EQ. 1) DWNV(1) =1.0
C PRINT OUT SPECTRAL INTERVALS
IF (IPRINT .GT. 1) THEN
WRITE (6,190)
DO 200 K=1,NSPECV
WRITE (6,210)K,WLNV(K),WNOV(K),BWNV(K)
& ,BWNV(K)+DWNV(K),DWNV(K)
200 CONTINUE
WRITE(6,320)
320 FORMAT (///' J WAVELN SOLAR FLUX N NP ['
& 16X,'A TERMS',14X,'] [',16X,'B TERMS',14X,']'/)
DO 300 J=1,24
K=1
SUM=ATERM(K,J)+ATERM(K+1,J)+ATERM(K+2,J)+ATERM(K+3,J)
WRITE(6,20)J,WLNV(J),SOLARF(J),NTERM(J),PEXPON(J)
&,(ATERM(K,J),K=1,4),(BTERM(K,J),K=1,4)
300 CONTINUE
20 FORMAT(1X,I2,F7.4,1PE14.5,0PI2,F6.3,3X,1P8E10.3)
END IF
210 FORMAT(1X,I3,F10.3,F10.2,F10.2,'-',F8.2,F10.3)
190 FORMAT(///' SPECTRAL INTERVALS'//
& ' SNUM MICRONS WAVENU INTERVAL',11X,'DELTA-WN')
C****** END SPECTRAL INTERVAL SET UP *************
RETURN
END
SUBROUTINE OPTCV(IPRINT)
PARAMETER (NLAYER=30,NLEVEL=31)
PARAMETER (NSPECI=46,NSPC1I=47,NSPECV=24,NSPC1V=25)
COMMON /ATM/ Z(NLEVEL),PRESS(NLEVEL),DEN(NLEVEL),TEMP(NLEVEL)
COMMON /GASS/ CH4(NLEVEL),XN2(NLEVEL),H2(NLEVEL),AR(NLEVEL)
& ,XMU(NLEVEL),GAS1(NLAYER),COLDEN(NLAYER)
COMMON /VISGAS/SOLARF(NSPECV),NTERM(NSPECV),PEXPON(NSPECV),
& ATERM(4,NSPECV),BTERM(4,NSPECV)
COMMON /AERSOL/ RADIUS(NLAYER), XNUMB(NLAYER)
& , REALI(NSPECI), XIMGI(NSPECI), REALV(NSPECV), XIMGV(NSPECV)
COMMON /CLOUD/ RADCLD(NLAYER), XNCLD(NLAYER)
& , RCLDI(NSPECI), XICLDI(NSPECI), RCLDV(NSPECV), XICLDV(NSPECV)
COMMON /TAUS/ TAUHI(NSPECI),TAUCI(NSPECI),TAUGI(NSPECI),
& TAURV(NSPECV),TAUHV(NSPECV),TAUCV(NSPECV),TAUGV(NSPECV)
COMMON /OPTICV/ DTAUV(NLAYER,NSPECV,4),TAUV(NLEVEL,NSPECV,4)
& ,WBARV(NLAYER,NSPECV,4), COSBV(NLAYER,NSPECV,4)
COMMON /SPECTV/ BWNV(NSPC1V),WNOV(NSPECV),DWNV(NSPECV)
& ,WLNV(NSPECV)
COMMON /PLANT/ CSUBP,RSFI,RSFV,F0PI
COMMON /ADJUST/ RHCH4,FH2,FHAZE,FHVIS,FHIR,TAUFAC,RCLOUD,FARGON
COMMON /CONST/RGAS,RHOP,PI,SIGMA
C*
C*
C THIS SUBROUTINE SETS THE OPTICAL CONSTANTS IN THE VISIBLE
C IT CALCUALTES FOR EACH LAYER, FOR EACH SPECRAL INTERVAL IN THE VIS
C LAYER: WBAR, DTAU, COSBAR
C LEVEL: TAU
C
DO 130 K=1,NSPECV
C LETS USE THE OPTICAL CONSTANTS FOR THOLIN
CALL THOLIN(WLNV(K),TNR,TNI)
REALV(K)=TNR
XIMGV(K)=TNI*FHVIS
C BUT WE NOW USE THE GEOMETRIC ALBEDO FITTED RESULTS
C XIMGV(K)=FITEDT(WLNV(K))
C XIMGV(K)=FITEDN(WLNV(K))
C THE CLOUD IS CLEAR IN THE VISIBLE
RCLDV(K)=1.27
XICLDV(K)=1.E-7
130 CONTINUE
C
DO 100 K=1,NSPECV
C ZERO THE COLUMN OPTICAL DEPTHS OF EACH TYPE
C ?FLAG? THE OPTICAL DEPTH OF THE TOP OF THE MODEL
C MAY NOT BE ZERO.
TAURV(K)=0.
TAUHV(K)=0.
TAUCV(K)=0.
TAUGV(K)=0.
DO 100 J=1,NLAYER
C1 HAZE
C CALL THE MIE CODE TO GIVE THE AEROSOL PROPERTIES
CALL XMIE(RADIUS(J),REALV(K),XIMGV(K),
& QEXT,QSCT,QABS,CBAR,WNOV(K))
TAEROS=QEXT*XNUMB(J)
C2 RAYLEIGH
C RAYLEIGH SCATTERING STRAIGHT FROM HANSEN AND TRAVIS...SEE NOTES
C RATIOED BY THE LAYER COLUMN NUMBER TO THE TOTAL
C COLUMN NUMBER ON EARTH. CM-2
TAURAY=(COLDEN(J)*28.9/(XMU(J)*1013.25))*
&(.008569/WLNV(K)**4)*(1.+.0113/WLNV(K)**2+.00013/WLNV(K)**4)
C3 CLOUD
C NEXT COMPUTE TAU CLOUD
TAUCLD=0.0
IF ( XNCLD(J) .GT. 0.) THEN
CALL XMIE(RADCLD(J),RCLDV(K),XICLDV(K),
& QEXTC,QSCTC,QABSC,CBARC,WNOV(K))
TAUCLD=QEXTC*XNCLD(J)
ENDIF
C
TAURV(K)=TAURV(K)+TAURAY
TAUHV(K)=TAUHV(K)+TAEROS
TAUCV(K)=TAUCV(K)+TAUCLD
C4 TAUGAS
C LOOP OVER THE NTERMS
DO 909 NT=1,NTERM(K)
TAUGAS=COLDEN(J)*GAS1(J)*BTERM(NT,K)*
& ( (PRESS(J+1) + PRESS(J))*.5 )**PEXPON(K)
C COMPUTE THE AVERAGE COSBAR AND WBAR
COSBV(J,K,NT)=(CBAR*TAEROS + CBARC*TAUCLD)
& /(TAEROS+TAUCLD+TAURAY)
DTAUV(J,K,NT)=TAUGAS+TAEROS+TAURAY+TAUCLD
TAUGV(K)=TAUGV(K)+TAUGAS*ATERM(NT,K)
C WE LET W RAYLEIGH BE .999 OR W=1 WHEN ONLY RAYLEIGH=PROBLEM FOR TRID
C WE HAVE ASSUMED ABOVE THAT COSBAR FOR RAYLEIGH IS ZERO.
WBARV(J,K,NT)=(QSCT*XNUMB(J)+TAURAY*0.9999999 + QSCTC*XNCLD(J) )
& /(TAUGAS+TAEROS+TAURAY+TAUCLD)
909 CONTINUE
100 CONTINUE
C TOTAL EXTINCTION OPTICAL DEPTHS
DO 119 K=1,NSPECV
C LOOP OVER NTERMS
DO 119 NT=1,NTERM(K)
TAUV(1,K,NT)=0.0
DO 119 J=1,NLAYER
TAUV(J+1,K,NT)=TAUV(J,K,NT)+DTAUV(J,K,NT)
119 CONTINUE
C
IF (IPRINT .GT. 1) THEN
WRITE (6,120)
120 FORMAT(///' OPTICAL CONSTANTS IN THE VISIBLE ')
DO 200 K=1,NSPECV
WRITE (6,190)
WRITE (6,210)K,WLNV(K),WNOV(K),BWNV(K)
& ,BWNV(K)+DWNV(K),DWNV(K)
WRITE (6,230)REALV(K),XIMGV(K)
DO 195 J=1,NLAYER
C RECALCULATE FOR PRINT OUT ONLY, ONLY FIRST NTERM
WRITE (6,220)XNUMB(J), WBARV(J,K,NT),COSBV(J,K,NT)
& ,DTAUV(J,K,NT),TAUV(J,K,NT)
195 CONTINUE
WRITE (6,240) TAUV(NLEVEL,K,NT)
200 CONTINUE
END IF
210 FORMAT(1X,I3,F10.3,F10.2,F10.2,'-',F8.2,F10.3)
190 FORMAT(1X//' SNUM MICRONS WAVENU INTERVAL DELTA-WN')
230 FORMAT(1X,'NREAL(LAYER)= ',1PE10.3,' NIMG(LAYER)= ',E10.3/
&' #AEROSOLS WBAR COSBAR DTAU TAU'
& ,9X,'RAY GAS AEROSOL')
220 FORMAT(8(1X,G9.3))
240 FORMAT(41X,G9.3)
RETURN
END
SUBROUTINE SFLUXV(IPRINT)
PARAMETER (NLAYER=30,NLEVEL=31)
PARAMETER (NSPECV=24,NSPC1V=25)
COMMON /VISGAS/SOLARF(NSPECV),NTERM(NSPECV),PEXPON(NSPECV),
& ATERM(4,NSPECV),BTERM(4,NSPECV)
COMMON /OPTICV/ DTAUV(NLAYER,NSPECV,4),TAUV(NLEVEL,NSPECV,4)
& ,WBARV(NLAYER,NSPECV,4), COSBV(NLAYER,NSPECV,4)
COMMON /SPECTV/ BWNV(NSPC1V),WNOV(NSPECV),DWNV(NSPECV)
& ,WLNV(NSPECV)
COMMON /FLUXV/ FNETV(NLEVEL), FUPV(NLEVEL,NSPECV),
& FDV(NLEVEL,NSPECV), FMNETV(NLEVEL),FMUPV(NLEVEL),FMDV(NLEVEL)
COMMON /PLANT/ CSUBP,RSFI,RSFV,F0PI
COMMON /UBARED/ UBARI,UBARV,UBAR0
DIMENSION FUW(NLEVEL),FDW(NLEVEL)
C ZERO THE NET FLUXES
DO 212 J=1,NLEVEL
FNETV(J)=-0.
FMNETV(J)=-0.
212 CONTINUE
C
C WE NOW ENTER A MAJOR LOOP OVER SPECRAL INTERVALS IN THE VISIBLE
C TO CALCULATE THE NET FLUX IN EACH SPECTRAL INTERVAL
C
C***************************************************************
DO 500 K=1,NSPECV
C ZERO THE SPECTRAL FLUXES IN ANTCIPATION OF SUMMING OVER NTERMS
DO 214 J=1,NLEVEL
FUPV(J,K)=0.
FDV(J,K)=0.
214 CONTINUE
C
C SET UP THE UPPER AND LOWER BOUNDARY CONDITIONS ON THE VISIBLE
F0PI=SOLARF(K)
BTOP=0.0
C
C LOOP OVER THE NTERMS BEGINING HERE
DO 912 NT=1,NTERM(K)
BSURF=0.+ RSFV*UBAR0*F0PI*EXP(-TAUV(NLEVEL,K,NT)/UBAR0)
C
C* WE CAN NOW SOLVE FOR THE COEFFICIENTS OF THE TWO STREAM
C CALL A SUBROUTINE THAT SOLVES FOR THE FLUX TERMS
C WITHIN EACH INTERVAL AT THE MIDPOINT WAVENUMBER
C
C FUW AND FDW ARE WORKING FLUX ARRAYS THAT WILL BE USED TO
C RETURN FLUXES FOR A GIVEN NT
C
C23456789012345678901234567890123456789012345678901234567890123456789012
CALL GFLUXV(NLEVEL,WNOV(K),DTAUV(1,K,NT),TAUV(1,K,NT),
& WBARV(1,K,NT),COSBV(1,K,NT),F0PI,RSFV,BTOP,BSURF,FUW,FDW,FMUPV,
& FMDV,IPRINT)
C
C NOW CALACULTE THE CUMULATIVE VISIBLE NET FLUX
DO 300 J=1,NLEVEL
FMNETV(J)=FMNETV(J)+( FMUPV(J)-FMDV(J) )*ATERM(NT,K)
FNETV(J)=FNETV(J)+( FUW(J)-FDW(J) )*ATERM(NT,K)
C AND THE SPECTRAL FLUXES SUMMED OVER THE NTERMS
FUPV(J,K)=FUPV(J,K)+FUW(J)*ATERM(NT,K)
FDV(J,K)=FDV(J,K)+FDW(J)*ATERM(NT,K)
300 CONTINUE
C
C
912 CONTINUE
500 CONTINUE
C*** END OF MAJOR SPECTRAL INTERVAL LOOP IN THE VISIBLE*****
RETURN
END
SUBROUTINE SFLUXI(IPRINT)
PARAMETER (NLAYER=30,NLEVEL=31)
PARAMETER (NSPECI=46,NSPC1I=47)
COMMON /ATM/ Z(NLEVEL),PRESS(NLEVEL),DEN(NLEVEL),TEMP(NLEVEL)
COMMON /OPTICI/ DTAUI(NLAYER,NSPECI),TAUI(NLEVEL,NSPECI)
& ,WBARI(NLAYER,NSPECI), COSBI(NLAYER,NSPECI)
COMMON /SPECTI/ BWNI(NSPC1I),WNOI(NSPECI),DWNI(NSPECI)
& ,WLNI(NSPECI)
COMMON /FLUXI/ FNETI(NLEVEL), FUPI(NLEVEL,NSPECI),
& FDI(NLEVEL,NSPECI),FMNETI(NLEVEL),FMUPI(NLEVEL),FMDI(NLEVEL)
COMMON /PLANT/ CSUBP,RSFI,RSFV,F0PI
COMMON /UBARED/ UBARI,UBARV,UBAR0
DATA PI/3.14159265358979323846/
CC ZERO THE NET FLUXES
DO 212 J=1,NLEVEL
FNETI(J)=0.
FMNETI(J)=0.
212 CONTINUE
C
C WE NOW ENTER A MAJOR LOOP OVER SPECRAL INTERVALS IN THE INFRARED
C TO CALCULATE THE NET FLUX IN EACH SPECTRAL INTERVAL
C
C***************************************************************
DO 501 K=1,NSPECI
C
C SET UP THE UPPER AND LOWER BOUNDARY CONDITIONS ON THE IR
C CALCULATE THE DOWNWELLING RADIATION AT THE TOP OF THE MODEL
C OR THE TOP LAYER WILL COOL TO SPACE UNPHYSICALLY
TTOP=TEMP(1)
TAUTOP=TAUI(NLEVEL,K)*PTOP/P0
BTOP=(1. - EXP(-TAUTOP/UBARI))*PLANCK(WNOI(K),TTOP,DWNI(K))
btop=0.
C SURFACE EMISSIONS
TSURF=TEMP(NLEVEL)
BSURF=PLANCK( WNOI(K), TSURF, DWNI(K) )
C* WE CAN NOW SOLVE FOR THE COEFFICIENTS OF THE TWO STREAM
C CALL A SUBROUTINE THAT SOLVES FOR THE FLUX TERMS
C WITHIN EACH INTERVAL AT THE MIDPOINT WAVENUMBER
C
C
CALL GFLUXI(NLEVEL,TEMP,WNOI(K),DWNI(K),DTAUI(1,K),WBARI(1,K)
& ,COSBI(1,K),RSFI,BTOP,BSURF,FUPI(1,K),FDI(1,K),
& FMUPI,FMDI,IPRINT)
C
C NOW CALACULTE THE CUMULATIVE IR NET FLUX
DO J=1,NLEVEL
C CORRECT FOR THE WAVENUMBER INTERVALS
FMNETI(J)=FMNETI(J)+( FMUPI(J)-FMDI(J) )*DWNI(K)
FNETI(J)=FNETI(J)+( FUPI(J,K)-FDI(J,K) )*DWNI(K)
ENDDO
C
501 CONTINUE
C*** END OF MAJOR SPECTRAL INTERVAL LOOP IN THE INFRARED****
RETURN
END
SUBROUTINE GFLUXI(NLL,TEMP,WAVEN,DW,DTAU,W0,COSBAR
& ,RSF,BTOP,BSURF,FP,FM,FMIDP,FMIDM,IPRINT)
PARAMETER (NL=35)
C* THIS SUBROUTINE TAKES THE OPTICAL CONSTANTS AND BOUNDARY CONDITONS
C FOR THE INFRARED FLUX AT ONE WAVELENGTH AND SOLVES FOR THE FLUXES AT
C THE LEVELS. THIS VERSION IS SET UP TO WORK WITH LAYER OPTICAL DEPTHS
C MEASURED FROM THE TOP OF EACH LAEYER. THE TOP OF EACH LAYER HAS
C OPTICAL DEPTH ZERO. IN THIS SUB LEVEL N IS ABOVE LAYER N. THAT IS LAYER N
C HAS LEVEL N ON TOP AND LEVEL N+1 ON BOTTOM. OPTICAL DEPTH INCREASES
C FROM TOP TO BOTTOM. SEE C.P. MCKAY, TGM NOTES.
REAL LAMDA
REAL*8 GAMA,CP,CM,CPM1,CMM1,E1,E2,E3,E4,B81,B82,R81,XK1,XK2
REAL*8 EM,EP
COMMON /UBARED/ UBARI,UBARV,UBAR0
DIMENSION TEMP(1),W0(1),COSBAR(1),DTAU(1),FM(1),FP(1)
& ,FMIDP(1),FMIDM(1)
DIMENSION B0(NL),B1(NL),ALPHA(NL),LAMDA(NL),XK1(NL),XK2(NL)
DIMENSION GAMA(NL),CP(NL),CM(NL),CPM1(NL),CMM1(NL)
&,E1(NL),E2(NL),E3(NL),E4(NL)
DATA PI/3.14159265358979323846/
C WE GO WITH THE HEMISPHERIC CONSTANT APPROACH
NLAYER=NLL-1
DO 20 J=1,NLAYER
ALPHA(J)=SQRT( (1.-W0(J))/(1.-W0(J)*COSBAR(J)) )
LAMDA(J)=ALPHA(J)*(1.-W0(J)*COSBAR(J))/UBARI
B1(J)=(PLANCK(WAVEN,TEMP(J+1),DW) - PLANCK(WAVEN,TEMP(J),DW))
& /(DTAU(J))
B0(J)=PLANCK(WAVEN,TEMP(J),DW)
20 CONTINUE
DO 7 J=1,NLAYER
GAMA(J)=(1-ALPHA(J))/(1.+ALPHA(J))
TERM=UBARI/(1.-W0(J)*COSBAR(J))
C* CP AND CM ARE THE CPLUS AND CMINUS TERMS EVALUATED AT THE
C BOTTOM OF THE LAYER. THAT IS AT DTAU OPTICAL DEPTH
CP(J)=B0(J)+B1(J)*DTAU(J) +B1(J)*TERM
CM(J)=B0(J)+B1(J)*DTAU(J) -B1(J)*TERM
C* CPM1 AND CMM1 ARE THE CPLUS AND CMINUS TERMS EVALUATED
C AT THE TOP OF THE LAYER, THAT IS ZERO OPTICAL DEPTH
CPM1(J)=B0(J)+B1(J)*TERM
CMM1(J)=B0(J)-B1(J)*TERM
7 CONTINUE
C*
C* NOW CALCULATE THE EXPONTNEIAL TERMS NEEDED
C* FOR THE TRIDIAGONAL ROTATED LAYERED METHOD
C* WARNING IF DTAU(J) IS GREATER THAN ABOUT 35
C* WE CLIPP IT TO AVOID OVERFLOW.
C* EXP (TAU) - EXP(-TAU) WILL BE NONSENSE THIS IS
C* CORRECTED IN THE DSOLVER ROUTINE. ?FLAG?
DO 8 J=1,NLAYER
C IF (LAMDA(J)*DTAU(J) .GE. 35.) WRITE(6,109) J,LAMDA(J)*DTAU(J)
C 109 FORMAT(' WATCHOUT FROM GFLUXI, DTAU .GT. 35:',I5,F20.4)
EP=EXP(35.)
C* NEED TO CLIPP THE EXPONENTIAL HERE.
IF ( LAMDA(J)*DTAU(J) .LT. 35. ) EP=EXP(LAMDA(J)*DTAU(J))
EM=1.0/EP
E1(J)=EP+GAMA(J)*EM
E2(J)=EP-GAMA(J)*EM
E3(J)=GAMA(J)*EP+EM
E4(J)=GAMA(J)*EP-EM
8 CONTINUE
C
B81=BTOP
B82=BSURF
R81=RSF
CALL DSOLVER(NLAYER,GAMA,CP,CM,CPM1,CMM1
& ,E1,E2,E3,E4,BTOP,BSURF,RSF,XK1,XK2)
C
C EVALUATE THE NLEVEL FLUXES THROUGH THE NLAYER LAYERS
C USE THE TOP (TAU=0) OPTICAL DEPTH EXPRESSIONS TO EVALUATE FP AND FM
C AT THE THE TOP OF EACH LAYER,J = LEVEL J
DO 50 J=1,NLAYER
FP(J)=XK1(J) + GAMA(J)*XK2(J) + CPM1(J)
FM(J)=GAMA(J)*XK1(J) + XK2(J) + CMM1(J)
50 CONTINUE
C USE EXPRESSION FOR BOTTOM FLUX TO GET THE FP AND FM AT NLEVEL
J=NLAYER
EP=EXP(35.)
C* NEED TO CLIPP THE EXPONENTIAL HERE.
IF ( LAMDA(J)*DTAU(J) .LT. 35. ) EP=EXP(LAMDA(J)*DTAU(J))
EM=1.0/EP
FP(J+1)=XK1(J)*EP + GAMA(J)*XK2(J)*EM + CP(J)
FM(J+1)=XK1(J)*EP*GAMA(J) + XK2(J)*EM + CM(J)
C TO GET THE FLUX WE NEED TO INTERGRATE OVER THE HEMISPHERE
C UBARI = .5 IS HEMISPHERIC CONSTANT
C UBARI = SQRT(3) IS GAUSS QUADRATURE
DO 60 J=1,NLL
FP(J)=FP(J)*2.*PI*UBARI
FM(J)=FM(J)*2.*PI*UBARI
60 CONTINUE
C
C NOW WE CALCULATE THE FLUXES AT THE MIDPOINTS OF THE LAYERS.
C
DO 1982 J=1,NLAYER
EP=EXP(35.)
C* NEED TO CLIPP THE EXPONENTIAL HERE.
IF (0.5*LAMDA(J)*DTAU(J) .LT. 35.)EP=EXP(0.5*LAMDA(J)*DTAU(J))
EM=1.0/EP
TERM=UBARI/(1.-W0(J)*COSBAR(J))
C* CP AND CM ARE THE CPLUS AND CMINUS TERMS EVALUATED AT THE
C BOTTOM OF THE LAYER. THAT IS AT DTAU OPTICAL DEPTH
CPMID=B0(J)+B1(J)*0.5*DTAU(J) +B1(J)*TERM
CMMID=B0(J)+B1(J)*0.5*DTAU(J) -B1(J)*TERM
FMIDP(J)=XK1(J)*EP + GAMA(J)*XK2(J)*EM + CPMID
FMIDM(J)=XK1(J)*EP*GAMA(J) + XK2(J)*EM + CMMID
C TO GET THE FLUX WE NEED TO INTERGRATE OVER THE HEMISPHERE
C UBARI = .5 IS HEMISPHERIC CONSTANT
C UBARI = SQRT(3) IS GAUSS QUADRATURE
FMIDP(J)=FMIDP(J)*2.*PI*UBARI
FMIDM(J)=FMIDM(J)*2.*PI*UBARI
1982 CONTINUE
IF (IPRINT .LT. 10) RETURN
WRITE(6,601)RSF,BTOP,BSURF
TAU=0.
DO 120 J=1,NLAYER
WRITE (6,301) TEMP(J),TAU,FP(J),FM(J),DTAU(J),W0(J),COSBAR(J)
& ,ALPHA(J), LAMDA(J),B0(J),B1(J)
TAU=TAU+DTAU(J)
120 CONTINUE
J=NLL
WRITE (6,301) TEMP(J),TAU,FP(J),FM(J)
WRITE (6,602)
DO 130 J=1,NLAYER
WRITE (6,301) CP(J),CM(J),E1(J),E2(J),E3(J),E4(J),CPM1(J)
& ,CMM1(J),GAMA(J)
130 CONTINUE
301 FORMAT(1X,1P13E10.3)
602 FORMAT(' CP(J),CM(J),E1(J),E2(J),EM1(J),EM2(J),CPM1(J),CMM1(J)'
& ,',GAMA(J)')
601 FORMAT(1X,'RSF,BTOP,BSURF= ',1P3E10.3,/
& ' TEMP(J),TAU,FUP,FDOWN,DTAU(J),W0(J),COSBAR(J)'
& ,',ALPHA(J),LAMDA(J),B0(J),B1(J)')
C C.P. McKAY
C234567..............................................................012
C*******************************************************
RETURN
END
SUBROUTINE GFLUXV(NTODO,WAVEN,DTDEL,TDEL,WDEL,CDEL
& , F0PI,RSF,BTOP,BSURF,FP,FM,FMIDP,FMIDM,IPRINT)
PARAMETER (NL=101)
C ?FLAG? THIS VALUE (101) MUST BE .GE. NTODO
C* THIS SUBROUTINE TAKES THE OPTICAL CONSTANTS AND BOUNDARY CONDITONS
C FOR THE VISIBLE FLUX AT ONE WAVELENGTH AND SOLVES FOR THE FLUXES AT
C THE LEVELS. THIS VERSION IS SET UP TO WORK WITH LAYER OPTICAL DEPTHS
C MEASURED FROM THE TOP OF EACH LAEYER. (DTAU) TOP OF EACH LAYER HAS
C OPTICAL DEPTH TAU(N).IN THIS SUB LEVEL N IS ABOVE LAYER N. THAT IS LAYER N
C HAS LEVEL N ON TOP AND LEVEL N+1 ON BOTTOM. OPTICAL DEPTH INCREASES
C FROM TOP TO BOTTOM. SEE C.P. MCKAY, TGM NOTES.
C THIS SUBROUTINE DIFFERS FROM ITS IR CONTERPART IN THAT HERE WE SOLVE FOR
C THE FLUXES DIRECTLY USING THE GENERALIZED NOTATION OF MEADOR AND WEAVOR
C J.A.S., 37, 630-642, 1980.
REAL LAMDA
REAL*8 GAMA,CP,CM,CPM1,CMM1,E1,E2,E3,E4,B81,B82,R81,XK1,XK2
REAL*8 EM,EP
COMMON /UBARED/ UBARI,UBARV,UBAR0
C THIS NEXT ROW OF VARIABLES ARE THOSE ACTUALLY USED IN THE
C ROUTINE
DIMENSION W0(NL), COSBAR(NL),DTAU(NL),TAU(NL)
DIMENSION WDEL(1),CDEL(1), DTDEL(1),TDEL(1),FM(1),FP(1)
& ,FMIDM(1),FMIDP(1)
DIMENSION ALPHA(NL),LAMDA(NL),XK1(NL),XK2(NL)
DIMENSION G1(NL),G2(NL),G3(NL)
DIMENSION GAMA(NL),CP(NL),CM(NL),CPM1(NL),CMM1(NL)
&,E1(NL),E2(NL),E3(NL),E4(NL),EXPTRM(NL)
DATA PI/3.14159265358979323846/
DATA IDELTA/1/
NAYER=NTODO-1
C TURN ON THE DELTA-FUNCTION IF REQUIRED HERE
IF (IDELTA .EQ. 0) THEN
DO J=1,NAYER
W0(J)=WDEL(J)
COSBAR(J)=CDEL(J)
DTAU(J)=DTDEL(J)
TAU(J)=TDEL(J)
ENDDO
TAU(NTODO)=TDEL(NTODO)
ELSE
C FOR THE DELTA FUNCTION HERE...
TAU(1)=TDEL(1)*(1.-WDEL(1)*CDEL(1)**2)
DO J=1,NAYER
W0(J)=WDEL(J)*(1.-CDEL(J)**2)/(1.-WDEL(J)*CDEL(J)**2)
COSBAR(J)=CDEL(J)/(1.+CDEL(J))
DTAU(J)=DTDEL(J)*(1.-WDEL(J)*CDEL(J)**2)
TAU(J+1)=TAU(J)+DTAU(J)
ENDDO
ENDIF
C WE GO WITH THE HEMISPHERIC CONSTANT APPROACH
C AS DEFINED BY M&W - THIS IS THE WAY THE IR IS DONE
DO 20 J=1,NAYER
ALPHA(J)=SQRT( (1.-W0(J))/(1.-W0(J)*COSBAR(J)) )
C THIS SET OF G'S IS FOR THE TWO STREAM HEMI-CONSTANT
C G1(J)=(1.-W0(J)*0.5*(1.+COSBAR(J)))/UBARV
C G2(J)=W0(J)*0.5*(1.-COSBAR(J))/UBARV
C G3(J)=0.5*(1.-COSBAR(J))
C SET OF CONSTANTS DETERMINED BY DOM
G1(J)= (SQRT(3.)*0.5)*(2. - W0(J)*(1.+COSBAR(J)))
G2(J)= (SQRT(3.)*W0(J)*0.5)*(1.-COSBAR(J))
G3(J)=0.5*(1.-SQRT(3.)*COSBAR(J)*UBAR0)
LAMDA(J)=SQRT(G1(J)**2 - G2(J)**2)
GAMA(J)=(G1(J)-LAMDA(J))/G2(J)
20 CONTINUE
DO 7 J=1,NAYER
G4=1.-G3(J)
DENOM=LAMDA(J)**2 - 1./UBAR0**2
C THERE IS A PONTENTIAL PROBLEM HERE IF W0=0 AND UBARV=UBAR0
C THEN DENOM WILL VANISH. THIS ONLY HAPPENS PHYSICALLY WHEN
C THE SCATTERING GOES TO ZERO??
C PREVENT THIS WITH A IF STATEMENT
IF ( DENOM .EQ. 0.) THEN
DENOM=1.E-10
WRITE (6,99)
99 FORMAT (' DENOM ZERO; RESET IN GFLUXV, W0=0?')
END IF
AM=F0PI*W0(J)*(G4 *(G1(J)+1./UBAR0) +G2(J)*G3(J) )/DENOM
AP=F0PI*W0(J)*(G3(J)*(G1(J)-1./UBAR0) +G2(J)*G4 )/DENOM
C* CPM1 AND CMM1 ARE THE CPLUS AND CMINUS TERMS EVALUATED
C AT THE TOP OF THE LAYER, THAT IS LOWER OPTICAL DEPTH TAU(J)
CPM1(J)=AP*EXP(-TAU(J)/UBAR0)
CMM1(J)=AM*EXP(-TAU(J)/UBAR0)
C* CP AND CM ARE THE CPLUS AND CMINUS TERMS EVALUATED AT THE
C BOTTOM OF THE LAYER. THAT IS AT HIGHER OPTICAL DEPTH TAU(J+1)
CP(J)=AP*EXP(-TAU(J+1)/UBAR0)
CM(J)=AM*EXP(-TAU(J+1)/UBAR0)
7 CONTINUE
C*
C* NOW CALCULATE THE EXPONTNEIAL TERMS NEEDED
C* FOR THE TRIDIAGONAL ROTATED LAYERED METHOD
C* WARNING IF DTAU(J) IS GREATER THAN ABOUT 35
C* WE CLIPP IT TO AVOID OVERFLOW.
C* EXP (TAU) - EXP(-TAU) WILL BE NONSENSE THIS IS
C* CORRECTED IN THE DSOLVER ROUTINE. ?FLAG?
DO J=1,NAYER
EXPTRM(J)=35.
IF ( LAMDA(J)*DTAU(J) .LT. 35. ) EXPTRM(J)=LAMDA(J)*DTAU(J)
ENDDO
C* NEED TO CLIPP THE EXPONENTIAL HERE.
DO 8 J=1,NAYER
EP=EXP(EXPTRM(J))
EM=1.0/EP
E1(J)=EP+GAMA(J)*EM
E2(J)=EP-GAMA(J)*EM
E3(J)=GAMA(J)*EP+EM
E4(J)=GAMA(J)*EP-EM
8 CONTINUE
C
B81=BTOP
B82=BSURF
R81=RSF
CALL DSOLVER(NAYER,GAMA,CP,CM,CPM1,CMM1
& ,E1,E2,E3,E4,B81,B82,R81,XK1,XK2)
C
C EVALUATE THE NTODO FLUXES THROUGH THE NAYER LAYERS
C USE THE TOP (TAU=0) OPTICAL DEPTH EXPRESSIONS TO EVALUATE FP AND FM
C AT THE THE TOP OF EACH LAYER,J = LEVEL J
DO 46 J=1,NAYER
FP(J)= XK1(J) + GAMA(J)*XK2(J) + CPM1(J)
FM(J)=GAMA(J)*XK1(J) + XK2(J) + CMM1(J)
46 CONTINUE
C USE EXPRESSION FOR BOTTOM FLUX TO GET THE FP AND FM AT NTODO
J=NAYER
EP=EXP(EXPTRM(J))
EM=1.0/EP
FP(J+1)=XK1(J)*EP + GAMA(J)*XK2(J)*EM + CP(J)
FM(J+1)=XK1(J)*EP*GAMA(J) + XK2(J)*EM + CM(J)
C NOTE THAT WE HAVE SOLVED FOR THE FLUXES DIRECTLY AND NO
C FURTHER INTEGRATION IS NEEDED, THE UBARV TERM IS ABSORBED
C INTO THE DEFINITION OF G1 THU G4
C UBARV = .5 IS HEMISPHERIC CONSTANT
C UBARV = SQRT(3) IS GAUSS QUADRATURE
C AND OTHER CASES AS PER MEADOR AND WEAVOR, JAS, 37, 630-643,1980.
C
C ADD THE DIRECT FLUX TERM TO THE DOWNWELLING RADIATION, LIOU 182
DO 80 J=1,NTODO
FM(J)=FM(J)+UBAR0*F0PI*EXP(-TAU(J)/UBAR0)
80 CONTINUE
C
C NOW WE CALCULATE THE FLUXES AT THE MIDPOINTS OF THE LAYERS.
C EXACLTY ANALOGOUS TO THE ABOVE COMPUTATION
C
DO 1982 J=1,NAYER
EP=EXP(0.5*EXPTRM(J))
EM=1.0/EP
G4=1.-G3(J)
DENOM=LAMDA(J)**2 - 1./UBAR0**2
C THERE IS A PONTENTIAL PROBLEM HERE IF W0=0 AND UBARV=UBAR0
C THEN DENOM WILL VANISH. THIS ONLY HAPPENS PHYSICALLY WHEN
C THE SCATTERING GOES TO ZERO??
C PREVENT THIS WITH A IF STATEMENT
IF ( DENOM .EQ. 0.) THEN
DENOM=1.E-10
WRITE (6,78)
78 FORMAT (' DENOM ZERO; RESET IN GFLUXV, W0=0?')
END IF
AM=F0PI*W0(J)*(G4 *(G1(J)+1./UBAR0) +G2(J)*G3(J) )/DENOM
AP=F0PI*W0(J)*(G3(J)*(G1(J)-1./UBAR0) +G2(J)*G4 )/DENOM
C* CPMID AND CMMID ARE THE CPLUS AND CMINUS TERMS EVALUATED
C AT THE MIDDLE OF THE LAYER, THAT IS LOWER OPTICAL DEPTH TAU(J)+
TAUMID= (TAU(J)+0.5*DTAU(J) )
CPMID=AP*EXP(-TAUMID/UBAR0)
CMMID=AM*EXP(-TAUMID/UBAR0)
FMIDP(J)=XK1(J)*EP + GAMA(J)*XK2(J)*EM + CPMID
FMIDM(J)=XK1(J)*EP*GAMA(J) + XK2(J)*EM + CMMID
C ADD THE DIRECT FLUX TO THE DOWNWELLING TERM
FMIDM(J)= FMIDM(J) +UBAR0*F0PI*EXP(-TAUMID/UBAR0)
1982 CONTINUE
C** NOW PRINTOUT IF NECESSARY
IF (IPRINT .LT. 9) RETURN
WRITE(6,601) F0PI,RSF,BTOP,BSURF
DO 120 J=1,NAYER
WRITE (6,301) TAU(J),FP(J),FM(J),DTAU(J),W0(J),COSBAR(J)
& ,ALPHA(J), LAMDA(J),G1(J),G2(J),G3(J)
120 CONTINUE
J=NTODO
WRITE (6,301) TAU(J),FP(J),FM(J)
WRITE (6,602)
DO 130 J=1,NAYER
WRITE (6,301) CP(J),CM(J),E1(J),E2(J),E3(J),E4(J),CPM1(J)
& ,CMM1(J),GAMA(J)
130 CONTINUE
301 FORMAT(1X,1P13E10.3)
602 FORMAT(' CP(J),CM(J),E1(J),E2(J),EM1(J),EM2(J),CPM1(J),CMM1(J)'
& ,',GAMA(J)')
601 FORMAT(1X,'F0PI,RSF,BTOP,BSURF= ',1P4E10.3,/
& ' TAU,FUP,FDOWN,DTAU(J),W0(J),COSBAR(J)'
& ,',ALPHA(J),LAMDA(J),G1,G2,G3')
C********************************************************************012
RETURN
END
FUNCTION AVERGE(X,Y)
AVERGE = .5*SQRT(X*Y)+0.25*(X+Y)
RETURN
END
SUBROUTINE PIA(K,TBAR,PNN,PCC,PCN,PHN)
PARAMETER (NSPECI=46, NTEMPS=14)
REAL*8 PIANN,PIACC,PIACN,PIAHN
COMMON /PIAC/ PIANN(NSPECI,NTEMPS),PIACC(NSPECI,NTEMPS)
& ,PIACN(NSPECI,NTEMPS),PIAHN(NSPECI,NTEMPS),TMIN,TMAX
TD=(TMAX-TMIN)/(NTEMPS-1)
IF (TBAR .LT. TMIN) TBAR=TMIN
IF (TBAR .GT. TMAX) TBAR=TMAX
DO 100 I=1,NTEMPS
T=TMIN + (I-1)*TD
IF (TBAR .GT. T) GO TO 100
FACTOR= (T-TBAR)/TD
PNN= PIANN(K,I) - FACTOR*(PIANN(K,I)-PIANN(K,I-1))
PCC= PIACC(K,I) - FACTOR*(PIACC(K,I)-PIACC(K,I-1))
PCN= PIACN(K,I) - FACTOR*(PIACN(K,I)-PIACN(K,I-1))
PHN= PIAHN(K,I) - FACTOR*(PIAHN(K,I)-PIAHN(K,I-1))
RETURN
100 CONTINUE
RETURN
END
SUBROUTINE SETPIA(IPRINT,IFLAG)
PARAMETER (NSPECI=46,NSPC1I=47, NTEMPS=14)
REAL*8 PIANN,PIACC,PIACN,PIAHN,T,W
COMMON /PIAC/ PIANN(NSPECI,NTEMPS),PIACC(NSPECI,NTEMPS)
& ,PIACN(NSPECI,NTEMPS),PIAHN(NSPECI,NTEMPS),TMIN,TMAX
COMMON /SPECTI/ BWNI(NSPC1I),WNOI(NSPECI),DWNI(NSPECI)
& ,WLNI(NSPECI)
DIMENSION T(NTEMPS),W(NSPECI)
DATA TMAX,TMIN/190.,60./
DO 11 K=1,NSPECI
DO 11 NT=1,NTEMPS
PIANN(K,NT)=0.0
PIACC(K,NT)=0.0
PIACN(K,NT)=0.0
PIAHN(K,NT)=0.0
11 CONTINUE
NWAVES=NSPECI
DO 99 INU=1,NSPECI
W(INU)=WNOI(INU)
IF (WNOI(INU) .LT. 1000.) GOTO 99
NWAVES=INU-1
GO TO 991
99 CONTINUE
991 CONTINUE
DT=(TMAX-TMIN)/(NTEMPS-1)
DO 12 I=1,NTEMPS
T(I)=TMIN + (I-1)*DT
12 CONTINUE
CALL REGIS (W,NSPECI,NWAVES,T,NTEMPS,NTEMPS,
& PIANN,PIACC,PIACN,PIAHN)
RETURN
CC TO BE DELETED...
1234 NT=NTEMPS
NNU=NSPECI
WRITE(6,101)
101 FORMAT(/55X,'N2-N2 ABSORPTION'//)
WRITE(6,100) (T(IT),IT=1,NT)
100 FORMAT(2X,'NU\T',13F9.0/)
DO 3 INU=1,NNU
3 WRITE(6,110) WNOI(INU),(PIANN(INU,IT),IT=1,NT)
WRITE(6,103)
103 FORMAT(/55X,'CH4-CH4 ABSORPTION'//)
WRITE(6,100) (T(IT),IT=1,NT)
DO 5 INU=1,NNU
5 WRITE(6,110) WNOI(INU),(PIACC(INU,IT),IT=1,NT)
WRITE(6,105)
105 FORMAT(/50X,'N2-CH4 + CH4-N2 ABSORPTION'//)
DO 7 INU=1,NNU
7 WRITE(6,110) WNOI(INU),(PIACN(INU,IT),IT=1,NT)
WRITE(6,102)
102 FORMAT(/55X,'N2-H2 + H2-N2 ABSORPTION'//)
WRITE(6,100) (T(IT),IT=1,NT)
DO 4 INU=1,NNU
4 WRITE(6,110) WNOI(INU),(PIAHN(INU,IT),IT=1,NT)
110 FORMAT(1X,F5.0,2X,13(1PD9.2))
RETURN
END
FUNCTION FITEDN(WAVELN)
C THIS FUNCTION USESS THE OUTPUT FROM THE GEOMETRIC
C ALBEDO FITTING PROGRAM TO DERIVE A SET OF
C INDECIES OF REFRACTION FOR THE HAZE IN THE VISIBLE
DIMENSION W(7),XK(7)
REWIND 021
DO K=1,7
READ (21,*) W(K),XK(K)
ENDDO
IF (WAVELN .LT. W(1)) THEN
CALL THOLIN(W(1),R1,XI1)
CALL THOLIN(WAVELN,R,XI)
FITEDN=XI*XK(1)/XI1
RETURN
ELSE
IF (WAVELN .GE. W(7)) THEN
FITEDN=XK(7)
RETURN
ELSE
C ALL INTERPOLATION IS IN LOG LAMBDA
DO 100 I=2,7
IF (WAVELN .LT. W(I) ) GO TO 101
100 CONTINUE
101 CONTINUE
C ALL INTERPOLATION IS IN LOG LAMBDA
FACTOR= (alog(WAVELN) - alog(W(I)) ) / (alog(W(I-1)) - alog(W(I)))
C IMAGINARY PART IS LOG INTERPOLATED
XNI=alog(XK(I)) + FACTOR*(alog(XK(I-1)) - alog(XK(I)))
FITEDN=exp(XNI)
RETURN
ENDIF
ENDIF
END
SUBROUTINE HAZE2(IPRINT,NINPUT,Z2,R2,XNUMB2)
C THIS ROUTINE SETS UP THE AEROSOL DISTRIBUTION
C
PARAMETER (NLHAZE=50)
COMMON /ADJUST/ RHCH4,FH2,FHAZE,FHVIS,FHIR,TAUFAC,RCLOUD,FARGON
DIMENSION Z(NLHAZE),DENH(NLHAZE),PRESS(NLHAZE),TEMP(NLHAZE)
DIMENSION CH4(NLHAZE),XN2(NLHAZE),H2(NLHAZE),AR(NLHAZE)
& ,XMU(NLHAZE)
DIMENSION RADIUS(NLHAZE), XNUMB(NLHAZE) ,TAU(NLHAZE)
DIMENSION PROD(NLHAZE)
DIMENSION Z2(1),R2(1),XNUMB2(1)
DATA PI/3.1415926535/
NLEVEL=NLHAZE
C
C* SET UP BACKGROUNG ATMOSHPERE
C FAC IS THE RATIO BETWEEN GRID SPACING TOP/BOTTOM (1 IS EVEN)
FAC=5.
ZTOP=670.
ZTOP1=635.
C
IP=IPRINT-3
CALL ZGRID(NLEVEL,FAC,ZTOP,ZTOP1,Z)
CALL ATMSETUP(NLEVEL,Z,RHCH4,FH2,FARGON,TEMP,PRESS,DENH,XMU,
& CH4,H2,XN2,AR,IP)
C
C
C
C LETS REDEFINE THE AEROSOLS BASED ON THE MCKAY SIMPLE
C AEROSOL MODEL... SEE NOTES
C PRODUCTION PROFILE FROM TOON ET AL. (1980)
XKB=1.380E-16
RHOH=1.
PX=1.E-7
RCRIT=0.09E-4
PY=PX*EXP(-1.124)
C
C
DO 500 J=1,NLEVEL
500 PROD(J)=FHAZE*2.E-7*3.5E-14*EXP(-0.5*( (PRESS(J)-PX)/PY)**2)
C
C DUE TO NOLINEAR AVERAGEING WE NEED TO DETERMINE THE TOTAL
C COLUMN PRODUCTION NUMERICALLY, AND ITERATIVELY
DO ITS=1,4
CLEVEL=0.0
DO J=2,NLEVEL
PRODMID=AVERGE( PROD(J),PROD(J-1) )
CLEVEL=CLEVEL+PRODMID*(Z(J-1)-Z(J))*1E5
ENDDO
FACTOR=FHAZE*3.5E-14/CLEVEL
C
C SCALE THE RESULTS TO ADJUST THE PRODUCTION RATE
DO J=1,NLEVEL
PROD(J)=PROD(J)*FACTOR
ENDDO
CC
ENDDO
C
TAU(1)=0.
IFLAG=0.
IA=6
IF (IPRINT .GT. 2) WRITE(IA,389)
389 FORMAT('1 LVL ALT PRESS TAU(.64) PROD CUM-PROD'
&,' DEN FALL-V R (CM)')
C.........
DO 501 J=1,NLEVEL
GASDEN=1.E6*PRESS(J)/(XKB*TEMP(J))
EG = 100.*EFFG(Z(J))
XMASS=XMU(J)/6.022E23
B=EG*RHOH*0.5*SQRT(PI/(2.*XMASS*XKB*TEMP(J)))/GASDEN
COAG=4.*SQRT(3.*XKB*TEMP(J)/RHOH)
IF (J .GT. 1 ) THEN
PRODMID=AVERGE( PROD(J),PROD(J-1) )
CMID=CLEVEL+0.5*PRODMID*(Z(J-1)-Z(J))*1E5
CLEVEL=CLEVEL+PRODMID*(Z(J-1)-Z(J))*1E5
C DO THE STOKES LIMIT ADJUSTMENT
C AS DISCRIBED IN WRITE UP.
PATH=1.689E14/GASDEN
XKNUD=PATH/RAD2
ETA=1.78E-4*SQRT(TEMP(J)/300.)
IK=1
IF (XKNUD .LT. 100.) THEN
B= (B*XKNUD**IK + 5.4*2.*RHOH*EG*RAD2/(9.*ETA) )/(XKNUD**IK +1.)
COAG=(COAG*XKNUD**IK + 4.*XKB*TEMP(J)/(3.*ETA*SQRT(RAD2)))
& /(XKNUD**IK +1.)
ENDIF
C END OF STOKES FLOW ADJUSTMENT
C
C PUT IN CHARGE REDUCTION IN COAGULATION ?FLAG?
COAG=COAG*EXP(-RAD2/RCRIT)
C END OF CHARGE ADJUST
R=( R**4.5 + ( 9.*COAG*CMID/(8.*PI*RHOH*B*B)
& + 9.*SQRT(R)*PRODMID/(8.*PI*RHOH*B*DEN1)) *
& (Z(J-1)-Z(J))*1.E5 )**(2./9.)
DEN1=CLEVEL*3./(4.*PI*RHOH*B*R**4)
C USE AVERAGEING PROCESS DENSITY AS DETERMINED IN NOTES
C TO GET VALUES FOR THE LAYERS
XNUMB(J-1)=AVERGE(DEN1,DEN2)*1.E5*( Z(J-1) - Z(J) )
CC ?FLAG? CHANGE THIS TO BE LEVEL VALUE SO THAT INTERPOLATON WILL
CC WORK
RADIUS(J)=1.E4 * R
C LETS DROP OUT THE HAZE FOR ALL ALTITUDES BELOW THE CONDENSATION POINT
C WITH A ADD-ON HERE. ?FLAG??
IF ( IFLAG .EQ. 1 ) XNUMB(J-1)=0.0
IF ( CH4(J)*PRESS(J)/PCH4(TEMP(J)) .GT. 0.9 ) IFLAG=1
ELSE
C TOP LAYER CALCULATIONS (J .EQ. 1)
CLEVEL=0.0
R= ( 9.*COAG*CLEVEL/(8.*PI*RHOH*B*B) )**(2./9.)
CC BRITNESS SCALE FACTOR FOR TOP LAYER ?FLAG??
R=1.E-7
DEN1=1.*FHAZE
RADIUS(1)=1.E4* R
ENDIF
RAD2=R
DEN2=DEN1
C
C CALL THE MIE CODE TO GIVE THE AEROSOL PROPERTIES AT A REF WAVELN
WVLN=.64
RMICRON=R*1E4
WNO=1.E4/WVLN
CALL THOLIN(WVLN,REAL,XIMG)
XIMG=XIMG*FHVIS
CALL XMIE(RMICRON,REAL,XIMG,
& QEXT,QSCT,QABS,CBAR,WNO)
DTAU=QEXT*XNUMB(J-1)
TAU(J)=TAU(J-1)+DTAU
C
IF (IPRINT .GT. 2)
& WRITE(IA,390) J,Z(J),PRESS(J),TAU(J),PROD(J),CLEVEL,DEN1,B,R
501 CONTINUE
390 FORMAT(I3,F6.0,1P7E10.2)
C
C NOW INTERPOLATE TO GET THE RETURNED VALUES.
DO 194 I=1,NINPUT
IF (Z2(I) .GT. Z(1) ) GO TO 194
J1=2
DO 195 J=J1,NLEVEL-1
IF (Z2(I) .LT. Z(J+1) )GO TO 195
FACTOR = (Z2(I) - Z(J)) / (Z(J+1) - Z(J))
XNUMB2(I)= TAU(J) + ( TAU(J+1) - TAU(J) ) * FACTOR
R2(I) = RADIUS(J) + (RADIUS(J+1) - RADIUS(J) ) * FACTOR
J1=J-1
GO TO 194
195 CONTINUE
194 CONTINUE
C
IF (IPRINT .GT. 0) THEN
WRITE(6,839)
839 FORMAT('1 LVL ALT TAU(.64) R (UM) IN LAYER')
WRITE(6,930) (J,Z2(J),XNUMB2(J),R2(J),J=1,NINPUT)
WRITE(6,930) NINPUT+1,Z2(NINPUT+1),TAU(NLEVEL)
930 FORMAT(0PI3,F6.1,1P2E10.2)
ENDIF
C
C NOW ADJUST THE NUMBER DEN SO THAT THE CORRECT OPACTIY IS
C OBTAINED
DO J=1,NINPUT-1
RMICRON=R2(J)
C CALL THE MIE CODE TO GIVE THE AEROSOL PROPERTIES AT A REF WAVELN
CALL XMIE(RMICRON,REAL,XIMG,
& QEXT,QSCT,QABS,CBAR,WNO)
XNUMB2(J)= (XNUMB2(J+1) - XNUMB2(J))/QEXT
ENDDO
C C?FLAG? QUICK FIX TO NLEVEL PROBLEM
XNUMB2(NINPUT)=(TAU(NLEVEL)-XNUMB2(NINPUT))/QEXT
RETURN
END
SUBROUTINE GAS2(J,K,T,P,U,TAUGAS)
PARAMETER (NLAYER=30)
COMMON /STRATO/ C2H2(NLAYER),C2H6(NLAYER)
DIMENSION A(18),A1(18),A2(18),B(18),B1(18),B2(18)
DIMENSION C(18),C1(18),C2(18)
DATA P0/0.1/
C P0 IS THE REFERENCE PRESSURE USED IN FORMULA
C PRESSURE IS IN BARS!!! ?FLAG?
C U IS IN MOLECULES PER CM2
DATA a /
&-58.07850,-73.66951,-22.00411,-26.61158,-27.04048,-25.84024,
& 0.00000,-18.70291,-16.08370,-23.03309,-31.18288,-19.96538,
&-24.37518,-17.94444, 0.00000,-17.46335,-20.93769,-10.74428/
DATA a1/
& 31.24782, 45.09350, 3.53728, 7.34218, 7.09808, 5.35252,
& 0.00000, 1.58486, -2.62994, 2.88087, 11.38614, 0.33491,
& 3.36078, -3.94288, 0.00000, -1.57041, 0.28770,-10.54978/
DATA a2/
& -6.15763, -9.47812, -0.72611, -1.73805, -1.65804, -1.23760,
& 0.00000, -0.33713, 0.84310, -0.44196, -2.79135, 0.00277,
& -0.60197, 1.27842, 0.00000, 0.37686, 0.08828, 2.77701/
DATA b/
& 5.27820, 2.66675, 1.31741, -5.05287, 1.88008, 4.49298,
& 0.00000, -0.94667, 2.77477, 15.80634, 33.83065, -2.29512,
& 2.90790, 0.77265, 0.00000, -0.95471, 3.30490, 0.43886/
DATA b1/
& -4.00372, -1.49356, -0.80081, 5.86571, -0.95994, -3.38754,
& 0.00000, 1.24800, -1.75869,-13.07462,-29.97486, 2.67427,
& -1.73763, 0.87904, 0.00000, 1.35597, -2.09025, 1.40393/
DATA b2/
& 0.81995, 0.32331, 0.16170, -1.42280, 0.26043, 0.87377,
& 0.00000, -0.31593, 0.36935, 2.86686, 6.83669, -0.75690,
& 0.24220, -0.46079, 0.00000, -0.42681, 0.32905, -0.61443/
DATA c/
& -8.99733, 6.51736, -0.54907, -4.80457, 2.83010, -2.12019,
& 0.00000, -1.53211,-13.48286,-18.91273,-16.89597, -1.02362,
& 0.25710, -6.27612, 0.00000, -0.56746, 0.14651, -8.53711/
DATA c1/
& 8.43865, -6.62977, 0.48722, 4.43910, -3.21534, 1.92581,
& 0.00000, 1.33927, 12.54849, 17.74174, 14.03467, 0.84636,
& -0.30260, 5.98928, 0.00000, 0.40909, -0.17410, 8.13665/
DATA c2/
& -1.99052, 1.64435, -0.13150, -1.02683, 0.80781, -0.49989,
& 0.00000, -0.30569, -2.91842, -4.14473, -2.91521, -0.18050,
& 0.06981, -1.43493, 0.00000, -0.07578, 0.03729, -1.93800/
XKV = A(K) + A1(K)*LOG10(T) + A2(K)*LOG10(T)**2
&+ LOG10(P/P0)*(B(K) + B1(K)*LOG10(T) + B2(K)*LOG10(T)**2)
&+ (C(K) + C1(K)*LOG10(T) + C2(K)*LOG10(T)**2)
& *LOG10(P/P0)**2
GASMIX=C2H2(J)
IF (K .GT. 11) GASMIX=C2H6(J)
XKV=10.**XKV
C CLIP OUT THE CLEAR INTERVALS.
IF (A(K) .EQ. 0.0) XKV=0.0
TAUGAS=U*GASMIX*XKV
RETURN
END
FUNCTION FITEDT(WAVELN)
C THIS FUNCTION USESS THE OUTPUT FROM THE GEOMETRIC
C ALBEDO FITTING PROGRAM TO DERIVE A SET OF
C INDECIES OF REFRACTION FOR THE HAZE IN THE VISIBLE
DIMENSION W(7),XK(7)
REWIND 021
DO K=1,7
READ (21,*) W(K),XK(K)
ENDDO
IF (WAVELN .LT. W(1)) THEN
CALL THOLIN(W(1),R1,XI1)
CALL THOLIN(WAVELN,R,XI)
FITEDT=XI*XK(1)/XI1
RETURN
ELSE
IF (WAVELN .GE. W(7)) THEN
CALL THOLIN(W(7),R1,XI1)
CALL THOLIN(WAVELN,R,XI)
FITEDT=XI*XK(7)/XI1
RETURN
ELSE
C ALL INTERPOLATION IS IN LOG LAMBDA
DO 100 I=2,7
IF (WAVELN .LT. W(I) ) GO TO 101
100 CONTINUE
101 CONTINUE
C ALL INTERPOLATION IS IN LOG LAMBDA
FACTOR= (alog(WAVELN) - alog(W(I)) ) / (alog(W(I-1)) - alog(W(I)))
C IMAGINARY PART IS LOG INTERPOLATED
XNI=alog(XK(I)) + FACTOR*(alog(XK(I-1)) - alog(XK(I)))
FITEDT=exp(XNI)
RETURN
ENDIF
ENDIF
END
SUBROUTINE ATMTEMP(IPRINT)
PARAMETER (NLAYER=30,NLEVEL=31)
PARAMETER (NSPECI=46,NSPC1I=47,NSPECV=24,NSPC1V=25)
COMMON /UBARED/ UBARI,UBARV,UBAR0
COMMON /ATM/ Z(NLEVEL),PRESS(NLEVEL),DEN(NLEVEL),TEMP(NLEVEL)
COMMON /LAPSE/ DTDP(NLAYER),CONVEQ
COMMON /GASS/ CH4(NLEVEL),XN2(NLEVEL),H2(NLEVEL),AR(NLEVEL)
& ,XMU(NLEVEL),GAS1(NLAYER),COLDEN(NLAYER)
COMMON /AERSOL/ RADIUS(NLAYER), XNUMB(NLAYER)
& , REALI(NSPECI), XIMGI(NSPECI), REALV(NSPECV), XIMGV(NSPECV)
COMMON /CLOUD/ RADCLD(NLAYER), XNCLD(NLAYER)
& , RCLDI(NSPECI), XICLDI(NSPECI), RCLDV(NSPECV), XICLDV(NSPECV)
COMMON /PLANT/ CSUBP,RSFI,RSFV,F0PI
COMMON /ADJUST/ RHCH4,FH2,FHAZE,FHVIS,FHIR,TAUFAC,RCLOUD,FARGON
COMMON /CONST/RGAS,RHOP,PI,SIGMA
DIMENSION XLNP(101)
DO J=1,NLEVEL
XLNP(J)=ALOG(PRESS(J))
ENDDO
C INPUTS:
C NLEVEL: NUMBER OF ALTITUDE LEVELS, J=1 IS AT THE TOP
C Z ALTITUDE GRID IN KM
C RHCH4 RELATIVE HUMIDITY OF METHANE AT Z=0.
C
C NOTES: METHANE FOLLOWS CONSTANT MIXING RATIO UNLESS SATUARTION
C VALUE IS LOWER.
C
C
C
C FIRST COMPUTE THE MEAN MOLECULAR WEIGHT AT EACH LEVEL
C NOW SET UP THE MIXING RATIOS OF THE GASSES BASED ON SATURATION
C CURVE OF CH4 AND SPECIFIED MEAN MOLECULAR WEIGHT AND H2 CONSTANT
CH4(NLEVEL)=PCH4(TEMP(NLEVEL))*RHCH4/PRESS(NLEVEL)
DO 134 J=NLEVEL-1,1,-1
CH4SAT=PCH4(TEMP(J))/PRESS(J)
CH4(J)=AMIN1(CH4SAT,CH4(NLEVEL),CH4(J+1))
134 CONTINUE
C ARGON AND H2 ARE NOT CHANGED AND ASSUMED TO BE SET
DO J=1,NLEVEL
XN2(J)=1.0 - H2(J) - CH4(J) -AR(J)
XMU(J)=28.0134*XN2(J)+2.158*H2(J)+16.0426*CH4(J)+39.948*AR(J)
ENDDO
C NOW THAT WE HAVE A NEW SET OF MEAN MOLECULAR WEIGHTS
C WE NEED TO RECOMPUTE THE ALT GRID.
C START AT THE BOTTOM LEVEL AND THEN CALCULATE THE
C ALTITUDE DROP ACROSS EACH LAYER
C GIVEN A T-P PROFILE CALCULATE THE HEIGHT
C PRESSURE IN BARS, T IN KELVIN, XMU IN AMUS, G IN M/SS
C RETURNS Z IN KILOMETERS
Z(NLEVEL)=0.0
DO 130 J=NLAYER,1,-1
C WE USE THE ALTITUDE BELOW (NOT AVG) THE LAYER TO GET G, SMALL ERROR
C IN THE SCALE HEIGHT H
C ASSUMES EQUATION OF STATE GIVEN BY IDEAL GAS LAW. ?FLAG?
C WHICH IS PROBABLY OK SINCE THE Z IS USED ONLY FOR G(Z)
BETA=(TEMP(J) - TEMP(J+1))/(TEMP(J+1)*XLNP(J) - TEMP(J)*XLNP(J+1))
TO=TEMP(J)/(1. + BETA * XLNP(J))
Z(J)=Z(J+1)-(RGAS*TO/(XMU(J)*EFFG(Z(J+1))))*(XLNP(J) - XLNP(J+1)+
& 0.5*BETA*(XLNP(J)**2 - XLNP(J+1)**2) )
130 CONTINUE
C NOW WE MUST RECOMPUTE THE OPTICALLY ACTIVE GASSES MASS MIXING RATIOS
CALL GASSES(IPRINT)
C CALL A ROTUINE THAT SETS UP THE AEROSOL DIST AND OPTICAL PROP.
C?? IF ( FHAZE .GT. 0.) CALL HAZE2(IPRINT,NLAYER,Z,RADIUS,XNUMB)
C?? IF (TAUFAC .GT. 0.) CALL CLD(IPRINT)
C
C* CALL A SUBROUTINE THAT SETS UP THE OPTICAL PROPERTIES IN THE
C INFRARED. AND THEN IN THE VISIBLE.
CALL OPTCI(IPRINT)
CALL OPTCV(IPRINT)
C* PRINT OUT THE AEROSOL AND CLOUD INFORMATION
IF (IPRINT .GT. 0) THEN
WRITE (6,139)
DO 135 J=1,NLAYER
WRITE(6,140)J,Z(J),PRESS(J),TEMP(J),
& CH4(J)*PRESS(J)/PCH4(TEMP(J)),CH4(J)*100.
135 CONTINUE
139 FORMAT(///' AEROSOLS AND CLOUDS AT LEVELS'/
&' LVL ALTITUDE P(BARS) TEMP RH-CH4'
& , ' %CH4 ')
140 FORMAT (1X,I3,F8.3,1PE10.3,0P2F7.2,F6.2,2(6X,0PF7.3,1PE10.3))
ENDIF
RETURN
END
FUNCTION PC2H6(T)
C THIS FUNCTION RETURNS THE VAPOR PRESSURE IN BARS.
C BASED ON THE FORMULA NUMBER (7) PAGE 11 OF
C NBS TECHNICAL PB-151 363
PC2H6=(1./760.)*10.**(8.5812- (965.4/T) )
RETURN
END
FUNCTION PC2H2(T)
C THIS FUNCTION RETURNS THE VAPOR PRESSURE IN BARS.
C BASED ON A LEAST SQUARES FIT TO DATA IN TABLE III
C NBS TECHNICAL PB-151 363
DATA A,B/21.68901,-2814.018/
PC2H2=(1./760.)*EXP(A+B/T)
RETURN
END
SUBROUTINE GASSES(IPRINT)
C THIS SUBROUTINE SETS UP THE MASS MIXING RATIOS OF THE
C OPTICALLY ACTIVE GASES: CH4, C2H2, AND C2H6
PARAMETER (NLAYER=30,NLEVEL=31)
COMMON /ATM/ Z(NLEVEL),PRESS(NLEVEL),DEN(NLEVEL),TEMP(NLEVEL)
COMMON /GASS/ CH4(NLEVEL),XN2(NLEVEL),H2(NLEVEL),AR(NLEVEL)
& ,XMU(NLEVEL),GAS1(NLAYER),COLDEN(NLAYER)
COMMON /STRATO/ C2H2(NLAYER),C2H6(NLAYER)
COMMON /ADJUST/ RHCH4,FH2,FHAZE,FHVIS,FHIR,TAUFAC,RCLOUD,FARGON
COMMON /CONST/RGAS,RHOP,PI,SIGMA
C*
C NOW CALCULATE THE LAYER AVERAGE GAS MIXING RATIOS.
C OF THE ABSORBING GAS IN UNITS OF GRAMS PER GRAM
C AND THE TOTAL LAYER COLUMN MASS GRAMS CM-2.
DO 159 J=1,NLAYER
EMU=(XMU(J+1)+XMU(J))*0.5
COLDEN(J)=RHOP*(PRESS(J+1)-PRESS(J))/EFFG(Z(J))
GAS1(J)=(16./EMU)*AVERGE(CH4(J+1),CH4(J))
159 CONTINUE
C WE NOW ALSO CALCULTE THE MASS MIXING RATIOS OF THE
C STRATOSPHERIC GASES USED IN THE IR WITHIN EACH LAYER.
J=1
FC2H2=2.E-6
FC2H6=2.E-5
C2H2(J) = MIN(FC2H2,PC2H2(TEMP(J))/PRESS(J))
C2H6(J) = MIN(FC2H6,PC2H6(TEMP(J))/PRESS(J))
DO J=2,NLAYER
C2H2(J) = MIN(FC2H2,PC2H2(TEMP(J))/PRESS(J),C2H2(J-1))
C2H6(J) = MIN(FC2H6,PC2H6(TEMP(J))/PRESS(J),C2H6(J-1))
ENDDO
C NOW CONVERT TO MASS MIXING RATIO
DO J=1,NLAYER
EMU=(XMU(J+1)+XMU(J))*0.5
C2H2(J)=C2H2(J)*26.0/EMU
C2H6(J)=C2H6(J)*30.0/EMU
ENDDO
C
IF (IPRINT .LT. 1) RETURN
WRITE (6,9)
9 FORMAT(///' ALT CH4 C2H2 C2H6: MASS MIXING RATIOS')
DO J=1,NLAYER
WRITE (6,10)Z(J),GAS1(J),C2H2(J),C2H6(J)
ENDDO
10 FORMAT(1X,F6.2,1P3E9.1)
RETURN
END
SUBROUTINE PRTOUT(IPRINT)
PARAMETER (NLAYER=30,NLEVEL=31)
PARAMETER (NSPECI=46,NSPC1I=47,NSPECV=24,NSPC1V=25)
COMMON /UBARED/ UBARI,UBARV,UBAR0
COMMON /ATM/ Z(NLEVEL),PRESS(NLEVEL),DEN(NLEVEL),TEMP(NLEVEL)
COMMON /LAPSE/ DTDP(NLAYER),CONVEQ
COMMON /GASS/ CH4(NLEVEL),XN2(NLEVEL),H2(NLEVEL),AR(NLEVEL)
& ,XMU(NLEVEL),GAS1(NLAYER),COLDEN(NLAYER)
COMMON /STRATO/ C2H2(NLAYER),C2H6(NLAYER)
COMMON /VISGAS/SOLARF(NSPECV),NTERM(NSPECV),PEXPON(NSPECV),
& ATERM(4,NSPECV),BTERM(4,NSPECV)
COMMON /AERSOL/ RADIUS(NLAYER), XNUMB(NLAYER)
& , REALI(NSPECI), XIMGI(NSPECI), REALV(NSPECV), XIMGV(NSPECV)
COMMON /CLOUD/ RADCLD(NLAYER), XNCLD(NLAYER)
& , RCLDI(NSPECI), XICLDI(NSPECI), RCLDV(NSPECV), XICLDV(NSPECV)
COMMON /TAUS/ TAUHI(NSPECI),TAUCI(NSPECI),TAUGI(NSPECI),
& TAURV(NSPECV),TAUHV(NSPECV),TAUCV(NSPECV),TAUGV(NSPECV)
COMMON /OPTICI/ DTAUI(NLAYER,NSPECI),TAUI(NLEVEL,NSPECI)
& ,WBARI(NLAYER,NSPECI), COSBI(NLAYER,NSPECI)
COMMON /SPECTI/ BWNI(NSPC1I),WNOI(NSPECI),DWNI(NSPECI)
& ,WLNI(NSPECI)
COMMON /FLUXI/ FNETI(NLEVEL), FUPI(NLEVEL,NSPECI),
& FDI(NLEVEL,NSPECI),FMNETI(NLEVEL),FMUPI(NLEVEL),FMDI(NLEVEL)
COMMON /OPTICV/ DTAUV(NLAYER,NSPECV,4),TAUV(NLEVEL,NSPECV,4)
& ,WBARV(NLAYER,NSPECV,4), COSBV(NLAYER,NSPECV,4)
COMMON /SPECTV/ BWNV(NSPC1V),WNOV(NSPECV),DWNV(NSPECV)
& ,WLNV(NSPECV)
COMMON /FLUXV/ FNETV(NLEVEL), FUPV(NLEVEL,NSPECV),
& FDV(NLEVEL,NSPECV), FMNETV(NLEVEL),FMUPV(NLEVEL),FMDV(NLEVEL)
COMMON /FLUX/ FNET(NLEVEL), FMNET(NLEVEL), THEAT(NLAYER)
COMMON /PLANT/ CSUBP,RSFI,RSFV,F0PI
COMMON /ADJUST/ RHCH4,FH2,FHAZE,FHVIS,FHIR,TAUFAC,RCLOUD,FARGON
COMMON /CONST/RGAS,RHOP,PI,SIGMA
COMMON /IO/ TIDAL
C
C* PRINT OUT THE AEROSOL AND CLOUD INFORMATION
IF (IPRINT .GT. 0) THEN
WRITE (6,139)
DO 135 J=1,NLAYER
WRITE(6,140)J,Z(J),PRESS(J),TEMP(J),
& CH4(J)*PRESS(J)/PCH4(TEMP(J)),CH4(J)*100.,RADIUS(J),XNUMB(J)
& ,RADCLD(J),XNCLD(J)
135 CONTINUE
139 FORMAT(///' AEROSOLS AND CLOUDS AT LEVELS'/
&' LVL ALTITUDE P(BARS) TEMP RH-CH4'
& , ' %CH4 HAZE:SIZE COL NUMB CLOUD:SIZE COL NUM')
140 FORMAT (1X,I3,F8.3,1PE10.3,0P2F7.2,F6.2,2(6X,0PF7.3,1PE10.3))
ENDIF
IF (IPRINT .GT. 0) THEN
C PRINT OUT SPECTRAL INTERVALS IN THE INFRARED
WRITE (6,190) 'INFR'
WRITE (6,210)(K,WLNI(K),WNOI(K),DWNI(K),TAUHI(K),
& TAUCI(K),0.0,TAUGI(K),TAUI(NLEVEL,K),REALI(K)
& ,XIMGI(K),RCLDI(K),XICLDI(K), K=1,NSPECI)
C PRINT OUT SPECTRAL INTERVALS IN THE VISIBLE
WRITE (6,190) 'VIS '
WRITE (6,210)(K,WLNV(K),WNOV(K),DWNV(K),TAUHV(K),
& TAUCV(K),TAURV(K),TAUGV(K),TAUV(NLEVEL,K,1),REALV(K)
& ,XIMGV(K),RCLDV(K),XICLDV(K), K=1,NSPECV)
END IF
210 FORMAT(1X,I2,F8.3,F9.1,F8.1,1P9E10.3,0P)
190 FORMAT(///' SPECTRAL INTERVALS IN ',A4//
& ' MICRONS WAVENU DEL-WN',
& ' TAU HAZE TAU CLOUD TAU RAY TAU GAS TAU-TOTAL ',
& 'HAZE:REAL IMG, CLOUD:REAL IMG')
C
C LETS OUTPUT THE SOLAR CONSTANT
SOLARC=0.0
DO 552 K=1,NSPECV
SOLARC=SOLARC+SOLARF(K)
552 CONTINUE
SOLALB=1.0 + 2.* FNETV(1)/SOLARC
C COMPUTE FRACTION OF TOTAL SUNLIT ABSORBED AT SURFACE, SEE NOTES
C THE 2 IS THE RATIO OF HEMISPHERIC AVERAGE AREA TO CROSS SECTIONAL
SURABS=-100.*2.*FNETV(NLEVEL)/SOLARC
WRITE(6,167) SOLARC, SOLALB,SURABS,100.*4.*FNETI(NLEVEL)/SOLARC
& ,RSFV,RSFI
167 FORMAT(/' SOLAR CONSTANT=', 1PE13.5,' BOND ALBEDO=',0PF7.4
&/' SURFACE ABSORPTION=',F7.4,'% OF TOTAL LIGHT INCIDENT ON TITAN'
&/' SURFACE RADIATION =',F7.4,'% OF TOTAL LIGHT INCIDENT ON TITAN'
& ,' RSFV=',F5.2,' RSFI=',F5.2)
C
C OUTPUT THE VISIBLE FLUXES AND PRESENT NET FLUX.
WRITE(6,892) (J,Z(J),TEMP(J),PRESS(J)
&,FNETV(J),FNETV(J)/(SOLARC*UBAR0),FNET(J),J=1,NLEVEL)
892 FORMAT(/' LEVEL ALT TEMP PRESS ',
& 'VIS FNET FNET/SO FNET-LEVEL'/
& (1X,I2,F7.2,2F8.3,1P3E11.3,0P))
C OUTPUT SELECTED SOLAR FLUX DATA FOR ALL SPECV'S
C NOTE THAT NLEVEL SHOULD NOT BE DIVISIBLE BY 6 OR FORMAT ERROR HERE
WRITE(6,289)' UP','SOLAR',(J,J=1,NLEVEL,NLEVEL/6),(K,WLNV(K),
& (FUPV(J,K),J=1,NLEVEL,NLEVEL/6), K=1,NSPECV)
WRITE(6,289)' DOWN','SOLAR',(J,J=1,NLEVEL,NLEVEL/6),(K,WLNV(K),
& (FDV(J,K),J=1,NLEVEL,NLEVEL/6), K=1,NSPECV)
C OUTPUT SELECTED IR FLUX DATA FOR ALL SPECI'S
WRITE(6,289)' UP','THERM',(J,J=1,NLEVEL,NLEVEL/6),(K,WLNI(K),
& (FUPI(J,K),J=1,NLEVEL,NLEVEL/6), K=1,NSPECI)
WRITE(6,289)' DOWN','THERM',(J,J=1,NLEVEL,NLEVEL/6),(K,WLNI(K),
& (FDI(J,K),J=1,NLEVEL,NLEVEL/6), K=1,NSPECI)
289 FORMAT (//A5,'WARD ',A5,' FLUXES',12X,'LEVELS'/
& ' WAVELENGTH',7I10,/(1X,I2,F8.3,1P7E10.3,0P))
C
RETURN
END
C ----------------------------------
C THIS IS THE START OF TGMSUBS.FOR
C ----------------------------------
SUBROUTINE OTMSETUP(NLEVEL,Z,RHCH4,FH2,FARGON,TEMP,PRESS,DEN,XMU,
& CH4,H2,XN2,AR,IPRINT)
C THIS IS THE OLD ONE WITH THE LINDAL SCALING
PARAMETER (NMAX=201)
C THIS SUBROUTINE SETS UP THE INTITAL ATMOSPHERIC PROFILE FOR TITAN
C BASED ON THE LINDAL ET AL DATA. THE ROUTINE STARTS WITH INPUTS
C INPUTS:
C NLEVEL: NUMBER OF ALTITUDE LEVELS, J=1 IS AT THE TOP
C Z ALTITUDE GRID IN KM
C RHCH4 RELATIVE HUMIDITY OF METHANE AT Z=0.
C FH2 MIXING RATIO BY NUMBER OF H2
C FARGON ARGON FLAG 0 = NO ARGON IN ATMSOPHERE
C -Y = ADJUST AR TO GIVE MEAN WEIGHT=Y
C +X = MIXING RATIO OF ARGON = X
C OUTPUTS: AT EACH LEVEL (NOT LAYER AVERAGES)
C TEMP (K), PRESS(BARS), DEN(CM-3), XMU = MEAN MOLECUALR WEIGHT
C CH4, H2, XN2, AR ARE THE NUMBER MIXING RATIOS OF THE GASES
C INTERNAL VARIABLES
C TLINAL, DLINAL, PLINAL : THE LINDAL INGRESS VALUES ON THE Z GRID
C
C NOTES: METHANE FOLLOWS CONSTANT MIXING RATIO UNLESS SATUARTION
C VALUE IS LOWER.
C
DIMENSION Z(NLEVEL),TEMP(NLEVEL),PRESS(NLEVEL),DEN(1),XMU(1)
DIMENSION CH4(1),H2(1),XN2(1),AR(1)
DIMENSION TLINAL(NMAX),DLINAL(NMAX),PLINAL(NMAX)
C
C FIRST SET UP THE LINDAL PURE N2 VALUES
CALL LINDAL(NLEVEL,Z,TLINAL,DLINAL,PLINAL)
C LOAD THE LINDAL VALUES INTO THE ARRAYS
DO J=1,NLEVEL
TEMP(J)= TLINAL(J)
DEN(J)= DLINAL(J)
PRESS(J)=PLINAL(J)
ENDDO
C
DO 1000 ITS =1,20
C
C NOW COMPUTE THE MEAN MOLECULAR WEIGHT AT EACH LEVEL
C NOW SET UP THE MIXING RATIOS OF THE GASSES BASED ON SATURATION
C CURVE OF CH4 AND SPECIFIED MEAN MOLECULAR WEIGHT AND H2 CONSTANT
CH4(NLEVEL)=PCH4(TEMP(NLEVEL))*RHCH4/PRESS(NLEVEL)
DO 134 J=NLEVEL-1,1,-1
CH4SAT=PCH4(TEMP(J))/PRESS(J)
CH4(J)=AMIN1(CH4SAT,CH4(NLEVEL),CH4(J+1))
134 CONTINUE
DO 20 J=1,NLEVEL
H2(J)=FH2
IF (FARGON .LT. 0.) THEN
C WE DECIDED TO KEEP CONSTANT MIXING RATIO = -FARGON
AR(J)=(-FARGON-28.0134+25.8554*H2(J)+11.9708*CH4(J))/11.9346
ELSE
IF (FARGON .EQ. 0.) THEN
C WE DECIDED TO DROP THE ARGON...
AR(J)=0.0
ELSE
C ARGON GIVEN BY A CONSTANT MIXING RATIO
AR(J)=FARGON
ENDIF
ENDIF
XN2(J)=1.0 - H2(J) - CH4(J) -AR(J)
XMU(J)=28.0134*XN2(J)+2.158*H2(J)+16.0426*CH4(J)+39.948*AR(J)
20 CONTINUE
C AT THIS POINT WE HAVE THE TEMP, PRESS, DEN AND XMU VALUES.
C ADJUST THE DEN BY THE FACTOR DUE TO THE MEAN REFRACTIVITY
C NOW LETS INTERGRATE THE DENSITY WITH ALTITUDE TO GET THE PRESS
SUMT=PLINAL(1)*6.02E23/10.
SUMB=SUMT
TLAST=TEMP(NLEVEL)
DO J=2,NLEVEL
C DENSITY ADJUSTMENT BASED ON REFRACTIVITIES ... SEE NOTES
DENF=294.1/(XN2(J)*294.1 + CH4(J)*410. + H2(J)*136. + AR(J)*277.8)
DEN(J) = DLINAL(J)*DENF
C PERFORM THE INTEGRALS LISTED IN NOTES
C SUMT IS ACTUAL PRESSURE
ADEN=(DEN(J)-DEN(J-1))/ALOG(DEN(J)/DEN(J-1))
SUMT=SUMT+(EFFG(Z(J))*ADEN)*( Z(J-1)-Z(J))*XMU(J)
ADEN=(DLINAL(J)-DLINAL(J-1))/ALOG(DLINAL(J)/DLINAL(J-1))
SUMB=SUMB+(EFFG(Z(J))*ADEN)*( Z(J-1)-Z(J))*28.01340
C
C NON IDEAL GAS CORRECTION IS 3.5% TIMES DT = 0.03% :NEGLECTED
PRESS(J)=PLINAL(J)*SUMT/SUMB
TEMP(J) =TLINAL(J)*(SUMT/SUMB)*(1./DENF)
C
ENDDO
30 CONTINUE
C
C HOW WELL WE DO ON CONVERGENCE
DT= ABS(TEMP(NLEVEL)-TLAST)
IF (DT .LT. 0.001) GO TO 1001
1000 CONTINUE
1001 IF (IPRINT .LT. 0) RETURN
WRITE (6,139)RHCH4,FH2,FARGON,DT
DO 135 J=1,NLEVEL-1
WRITE(6,140)J,Z(J),PRESS(J),DEN(J),TEMP(J),
& CH4(J)*PRESS(J)/PCH4(TEMP(J))
& ,CH4(J)*100.,XN2(J)*100.,H2(J)*100.,AR(J)*100.,XMU(J)
& ,(TEMP(J+1)-TEMP(J))/(Z(J+1)-Z(J))
135 CONTINUE
J=NLEVEL
WRITE(6,140)J,Z(J),PRESS(J),DEN(J),TEMP(J),
& CH4(J)*PRESS(J)/PCH4(TEMP(J))
& ,CH4(J)*100.,XN2(J)*100.,H2(J)*100.,AR(J)*100.,XMU(J)
139 FORMAT(///' BACKGROUNG ATMOSPHERE AT LEVELS'/
& ' SURFACE HUMIDITY OF CH4:',F5.3,' H2 MIXING RATIO:',F6.4,
& ' ARGON SETTING:',F8.4/' FINAL CONVERGENCE ON TEMP:',F10.5
& ' LINDAL ET AL SCALING'/
&' LVL ALTITUDE P(BARS) DEN(CM-3) TEMP RH-CH4'
& , ' %CH4 %N2 %H2 %AR MU DT/DZ' )
140 FORMAT(1X,I3,F8.3,1P2E10.3,0PF7.2,F5.2,2F6.2,2F5.2,4F6.2)
C END INTITAL BACKGROUND ATMOSPHERE SETUP FOR TITAN
RETURN
END
SUBROUTINE ATMSETUP(NLEVEL,Z,RHCH4,FH2,FARGON,TEMP,PRESS,DEN,XMU,
& CH4,H2,XN2,AR,IPRINT)
C THIS IS THE NEW ONE WITH THE BWC EQ OF STATE
PARAMETER (NMAX=201)
C THIS SUBROUTINE SETS UP THE INTITAL ATMOSPHERIC PROFILE FOR TITAN
C BASED ON THE LINDAL ET AL DATA. THE ROUTINE STARTS WITH INPUTS
C INPUTS:
C NLEVEL: NUMBER OF ALTITUDE LEVELS, J=1 IS AT THE TOP
C Z ALTITUDE GRID IN KM
C RHCH4 RELATIVE HUMIDITY OF METHANE AT Z=0.
C FH2 MIXING RATIO BY NUMBER OF H2
C FARGON ARGON FLAG 0 = NO ARGON IN ATMSOPHERE
C -Y = ADJUST AR TO GIVE MEAN WEIGHT=Y
C +X = MIXING RATIO OF ARGON = X
C OUTPUTS: AT EACH LEVEL (NOT LAYER AVERAGES)
C TEMP (K), PRESS(BARS), DEN(CM-3), XMU = MEAN MOLECUALR WEIGHT
C CH4, H2, XN2, AR ARE THE NUMBER MIXING RATIOS OF THE GASES
C INTERNAL VARIABLES
C TLINAL, DLINAL, PLINAL : THE LINDAL INGRESS VALUES ON THE Z GRID
C
C NOTES: METHANE FOLLOWS CONSTANT MIXING RATIO UNLESS SATUARTION
C VALUE IS LOWER.
C
DIMENSION Z(NLEVEL),TEMP(NLEVEL),PRESS(NLEVEL),DEN(1),XMU(1)
DIMENSION CH4(1),H2(1),XN2(1),AR(1)
DIMENSION TLINAL(NMAX),DLINAL(NMAX),PLINAL(NMAX)
C
C FIRST SET UP THE LINDAL PURE N2 VALUES
CALL LINDAL(NLEVEL,Z,TLINAL,DLINAL,PLINAL)
C LOAD THE LINDAL VALUES INTO THE ARRAYS
DO J=1,NLEVEL
TEMP(J)= TLINAL(J)
DEN(J)= DLINAL(J)
PRESS(J)=PLINAL(J)
ENDDO
C
DO 1000 ITS =1,20
C
C NOW COMPUTE THE MEAN MOLECULAR WEIGHT AT EACH LEVEL
C NOW SET UP THE MIXING RATIOS OF THE GASSES BASED ON SATURATION
C CURVE OF CH4 AND SPECIFIED MEAN MOLECULAR WEIGHT AND H2 CONSTANT
CH4(NLEVEL)=PCH4(TEMP(NLEVEL))*RHCH4/PRESS(NLEVEL)
DO 134 J=NLEVEL-1,1,-1
CH4SAT=PCH4(TEMP(J))/PRESS(J)
CH4(J)=AMIN1(CH4SAT,CH4(NLEVEL),CH4(J+1))
134 CONTINUE
DO 20 J=1,NLEVEL
H2(J)=FH2
IF (FARGON .LT. 0.) THEN
C WE DECIDED TO KEEP CONSTANT MIXING RATIO = -FARGON
AR(J)=(-FARGON-28.0134+25.8554*H2(J)+11.9708*CH4(J))/11.9346
ELSE
IF (FARGON .EQ. 0.) THEN
C WE DECIDED TO DROP THE ARGON...
AR(J)=0.0
ELSE
C ARGON GIVEN BY A CONSTANT MIXING RATIO
AR(J)=FARGON
ENDIF
ENDIF
XN2(J)=1.0 - H2(J) - CH4(J) -AR(J)
XMU(J)=28.0134*XN2(J)+2.158*H2(J)+16.0426*CH4(J)+39.948*AR(J)
20 CONTINUE
C AT THIS POINT WE HAVE THE TEMP, PRESS, DEN AND XMU VALUES.
C ADJUST THE DEN BY THE FACTOR DUE TO THE MEAN REFRACTIVITY
C NOW LETS INTERGRATE THE DENSITY WITH ALTITUDE TO GET THE PRESS
SUMT=PLINAL(1)*6.02E23/10.
SUMB=SUMT
TLAST=TEMP(NLEVEL)
DO J=2,NLEVEL
C DENSITY ADJUSTMENT BASED ON REFRACTIVITIES ... SEE NOTES
DENF=294.1/(XN2(J)*294.1 + CH4(J)*410. + H2(J)*136. + AR(J)*277.8)
DEN(J) = DLINAL(J)*DENF
C PERFORM THE INTEGRALS LISTED IN NOTES
C SUMT IS ACTUAL PRESSURE
ADEN=(DEN(J)-DEN(J-1))/ALOG(DEN(J)/DEN(J-1))
SUMT=SUMT+(EFFG(Z(J))*ADEN)*( Z(J-1)-Z(J))*XMU(J)
ADEN=(DLINAL(J)-DLINAL(J-1))/ALOG(DLINAL(J)/DLINAL(J-1))
SUMB=SUMB+(EFFG(Z(J))*ADEN)*( Z(J-1)-Z(J))*28.01340
C
C
C PRESSURE IS IN BARS
C DETERMINE THE TEMPERATURE WITH BWR EQ OF STATE
PRESS(J)=SUMT*1.E7/6.02E23 *1E-6
TEMP(J) =EQSTATE(XN2(J),CH4(J),H2(J),AR(J),PRESS(J),DEN(J))
C
C
ENDDO
30 CONTINUE
C
C HOW WELL WE DO ON CONVERGENCE
DT= ABS(TEMP(NLEVEL)-TLAST)
IF (DT .LT. 0.001) GO TO 1001
1000 CONTINUE
1001 IF (IPRINT .LT. 0) RETURN
WRITE (6,139)RHCH4,FH2,FARGON,DT
DO 135 J=1,NLEVEL-1
WRITE(6,140)J,Z(J),PRESS(J),DEN(J),TEMP(J),
& CH4(J)*PRESS(J)/PCH4(TEMP(J))
& ,CH4(J)*100.,XN2(J)*100.,H2(J)*100.,AR(J)*100.,XMU(J)
& ,(TEMP(J+1)-TEMP(J))/(Z(J+1)-Z(J))
135 CONTINUE
J=NLEVEL
WRITE(6,140)J,Z(J),PRESS(J),DEN(J),TEMP(J),
& CH4(J)*PRESS(J)/PCH4(TEMP(J))
& ,CH4(J)*100.,XN2(J)*100.,H2(J)*100.,AR(J)*100.,XMU(J)
139 FORMAT(///' BACKGROUNG ATMOSPHERE AT LEVELS'/
& ' SURFACE HUMIDITY OF CH4:',F5.3,' H2 MIXING RATIO:',F6.4,
& ' ARGON SETTING:',F8.4/' FINAL CONVERGENCE ON TEMP:',F10.5/
&' LVL ALTITUDE P(BARS) DEN(CM-3) TEMP RH-CH4'
& , ' %CH4 %N2 %H2 %AR MU DT/DZ' )
140 FORMAT(1X,I3,F8.3,1P2E10.3,0PF7.2,F5.2,2F6.2,2F5.2,4F6.2)
C END INTITAL BACKGROUND ATMOSPHERE SETUP FOR TITAN
RETURN
END
FUNCTION EQSTATE(FN2,FCH4,FH2,FAR,PIN,DENIN)
PARAMETER (NGAS=2)
DIMENSION BO(NGAS),AO(NGAS),CO(NGAS),DO(NGAS),B(NGAS),A(NGAS),
& C(NGAS),DELTA(NGAS),GAMMA(NGAS),ALPHA(NGAS),GAS(NGAS)
C
C NITROGEN
DATA BO(1),AO(1),CO(1),DO(1),B(1),
& A(1),C(1),DELTA(1),GAMMA(1),ALPHA(1)/
& 0.0484824, 1.27389, 4273.00, 7.61781E6, 0.00232373,
& 0.0178444, 475.000, 0.8320E6, 0.006500, 0.00015300/
C METHANE
DATA BO(2),AO(2),CO(2),DO(2),B(2),
& A(2),C(2),DELTA(2),GAMMA(2),ALPHA(2)/
& 0.042600, 1.85500, 22570.0, 0.000000, 0.00338004,
& 0.0494000, 2545.00, 0.00000, 0.006000, 0.000124359/
C ABOVE NUMBERS FROM IGT PUBLICATION: ELLINGTION ET AL, 1959
C
C INPUTS: DEN=DENSITY MOLECULES/CM3, PRESS = BARS,
C OUTPUT: EQSTATE IS THE TEMPERATURE OF THE MIXTURE
C
DATA RGAS/0.08206/
C R IS DIFFERENT ON ACCOUNT OF ATMS INSTEAD OF BARS
C CONVERT TO CGS UNITS DEN= MOLES/LITERS, PRESS=ATMS
DEN =1.E3*DENIN/6.0220E23
PRESS = PIN/1.01325
C
C BENEDICT, WEBB AND RUBIN EQUATION OF STATE (ELLINGTON ET AL.,1959)
C TO COMPUTE THE VIRIAL COEFFICEINTS (REIF, P422)
C
C FIRST SET UP THE COEFFICIENTS FOR THE MIXING RATIOS
GAS(1)=FN2+FH2+FAR
GAS(2)=FCH4
C ZERO THE MEANS
BOM = 0.
AOM = 0.
COM = 0.
DOM = 0.
AM = 0.
BM = 0.
CM = 0.
DELTAM = 0.
GAMMAM = 0.
ALPHAM = 0.
C COMPUTE SUMS
DO 100 N=1,NGAS
BOM = BOM + GAS(N)*BO(N)
AOM = AOM + GAS(N)*SQRT(AO(N))
COM = COM + GAS(N)*SQRT(CO(N))
DOM = DOM + GAS(N)*SQRT(DO(N))
AM = AM + GAS(N)*A(N)**(1./3.)
BM = BM + GAS(N)*B(N)**(1./3.)
CM = CM + GAS(N)*C(N)**(1./3.)
DELTAM = DELTAM + GAS(N)*DELTA(N)**(1./3.)
GAMMAM = GAMMAM + GAS(N)*GAMMA(N)**(1./3.)
ALPHAM = ALPHAM + GAS(N)*ALPHA(N)**(1./3.)
100 CONTINUE
C
AOM = AOM**2
COM = COM**2
DOM = DOM**2
AM = AM**3
BM = BM**3
CM = CM**3
DELTAM = DELTAM**3
GAMMAM = GAMMAM**3
ALPHAM = ALPHAM**3
C WE NOW HAVE THE TERMS FOR THE EQUATION OF STATE
C NEXT GET THE IDEAL GAS TEMPERATURE
T1=PRESS/(RGAS*DEN)
T=T1
EQSTATE=T
DO ITS =1,50
C1=BOM*RGAS*T - AOM - COM/T**2
C2=(BM*RGAS*T-AM)+CM*EXP(-GAMMAM*DEN**2)*(1.+GAMMAM*DEN**2)/(T*T)
C3=AM*ALPHAM
PNEW=RGAS*DEN*T+C1*DEN**2 + C2*DEN**3 + C3*DEN**6
T= T*PRESS/PNEW
IF (ABS((T-EQSTATE)/T) .LT. 1.E-5) RETURN
EQSTATE=T
ENDDO
RETURN
END
FUNCTION PCH4(T)
C THIS SUBROUTINE RETURNS THE VAPOR PRESSURE OF METHANE
C OVER THE TEMPERATURE INTERVAL (74-97K) IN UNITS OF BARS
C IT IS BASED ON DATA FROM THE MATHASEN GAS DATA BOOK P 351
C FITTED BY CP MCKAY 10/85
PCH4=3.4543E4 * EXP(-1145.705/T)
RETURN
END
FUNCTION TITANT(P)
DIMENSION PRESS(107),TEMP(107)
C THIS SUBROUTINE RETURNS THE LINDAL ET AL INGRESS
C T,P AND DEN PROFILE ON THE INPUT ZGRID.
DIMENSION ZLIN(106),DLIN(106)
CALL INGRESS(TEMP(2),ZLIN,DLIN,PRESS(2))
DATA PRESS(1)/0.0/
TEMP(1)=TEMP(2)
DO 100 I=2,107
IF (PRESS(I) .GT. P) GO TO 101
100 CONTINUE
101 CONTINUE
FACTOR= (PRESS(I) - P )/(PRESS(I) - PRESS(I-1))
TITANT=TEMP(I) + FACTOR*(TEMP(I-1) - TEMP(I))
IF ( P .GT. PRESS(107)) TITANT=94.0
RETURN
END
SUBROUTINE THOLIN(WAVELN,XNR,XNI)
DIMENSION W(90),XN(90),XK(90)
DATA W/
&920.0000,850.0000,774.9000,688.8000,563.5000,387.4000,229.6000,
&172.2000,140.9000,121.5000, 81.5700, 56.3500, 36.4600, 31.0000,
& 22.1400, 18.2300, 17.7100, 14.4200, 12.9100, 11.7000, 11.0700,
& 10.5100, 8.7310, 7.6530, 7.0440, 6.6660, 6.4570, 6.3260,
& 5.9610, 5.8480, 5.7400, 5.4380, 5.2530, 5.1660, 4.8810,
& 4.6260, 4.5920, 4.4920, 4.4280, 4.2170, 3.9480, 3.6680,
& 3.4630, 3.2460, 3.0090, 2.9380, 2.8180, 2.7430, 2.6950,
& 2.4220, 2.4030, 2.3930, 2.2140, 2.0190, 1.8730, 1.8230,
& 1.8130, 1.8020, 1.3810, 1.3570, 1.1920, 1.1480, 1.0160,
& 0.8731, 0.6888, 0.5635, 0.4428, 0.4133, 0.3874, 0.3542,
& 0.3263, 0.2952, 0.2638, 0.2384, 0.1968, 0.1631, 0.1362,
& 0.1215, 0.1181, 0.1159, 0.1097, 0.1016, 0.0925, 0.0800,
& 0.0785, 0.0588, 0.0449, 0.0415, 0.0312, 0.0207/
DATA XN/
&2.170,2.170,2.160,2.160,2.150,2.120,2.070,2.040,2.030,2.020,1.930,
&1.860,1.810,1.810,1.800,1.760,1.740,1.670,1.670,1.640,1.660,1.670,
&1.710,1.720,1.690,1.690,1.640,1.580,1.430,1.440,1.480,1.550,1.580,
&1.580,1.610,1.620,1.610,1.610,1.610,1.630,1.640,1.650,1.650,1.650,
&1.610,1.590,1.580,1.590,1.600,1.620,1.620,1.620,1.630,1.630,1.630,
&1.640,1.640,1.640,1.640,1.640,1.650,1.650,1.650,1.660,1.680,1.700,
&1.720,1.690,1.660,1.630,1.640,1.660,1.680,1.680,1.660,1.650,1.700,
&1.740,1.750,1.750,1.720,1.670,1.580,1.370,1.330,0.963,0.812,0.802,
&0.850,0.920/
DATA XK/
&3.0E-3,1.0E-2,2.9E-2,4.7E-2,7.0E-2,1.0E-1,1.4E-1,1.6E-1,1.6E-1,
&1.9E-1,2.1E-1,1.9E-1,1.5E-1,1.4E-1,1.8E-1,2.1E-1,2.1E-1,1.7E-1,
&1.4E-1,9.7E-2,7.9E-2,7.5E-2,9.2E-2,1.3E-1,1.7E-1,2.2E-1,2.6E-1,
&2.8E-1,1.5E-1,7.0E-2,2.9E-2,1.1E-2,8.7E-3,7.6E-3,1.0E-2,2.7E-2,
&2.8E-2,1.4E-2,1.1E-2,1.0E-2,1.3E-2,2.1E-2,3.5E-2,5.6E-2,7.5E-2,
&6.0E-2,2.4E-2,1.1E-2,4.1E-3,1.2E-3,8.5E-4,8.0E-4,8.9E-4,7.2E-4,
&5.2E-4,4.4E-4,4.2E-4,4.0E-4,4.1E-4,4.2E-4,5.2E-4,6.4E-4,1.0E-3,
&2.4E-3,8.8E-3,2.3E-2,6.0E-2,7.6E-2,9.1E-2,1.1E-1,1.3E-1,1.5E-1,
&1.8E-1,2.1E-1,2.2E-1,2.4E-1,2.7E-1,3.7E-1,4.0E-1,4.3E-1,5.0E-1,
&5.8E-1,6.7E-1,7.7E-1,7.7E-1,6.2E-1,3.8E-1,3.1E-1,1.4E-1,4.9E-2/
XNR=XN(1)
XNI=XK(1)
IF (WAVELN .GT. W(1)) RETURN
XNR=XN(90)
XNI=XK(90)
IF (WAVELN .LT. W(90)) RETURN
DO 100 I=2,90
IF (WAVELN .GT. W(I) ) GO TO 101
100 CONTINUE
101 CONTINUE
C ALL INTERPOLATION IS IN LOG LAMBDA
FACTOR= (alog(WAVELN) - alog(W(I)) ) / (alog(W(I-1)) - alog(W(I)))
C REAL PART IS LINEARLY INTERPOLATED
XNR=XN(I) + FACTOR*(XN(I-1) - XN(I))
C IMAGINARY PART IS LOG INTERPOLATED
XNI=alog(XK(I)) + FACTOR*(alog(XK(I-1)) - alog(XK(I)))
XNI=exp(XNI)
RETURN
END
FUNCTION REFLIQ(W)
DIMENSION WAVENO(55),XIMG(55)
DATA WAVENO/0., 10., 40., 60., 80., 100., 120., 140.,
&160., 180., 200., 220., 240., 260., 280., 300., 325.,
&340., 350., 360., 380., 400., 450., 480., 500., 520.,
&540., 560., 600., 640., 684., 720., 760., 780., 800.,
&840., 875., 920., 960., 1000., 1040., 1080., 1120., 1160.,
&1200., 1220., 1240., 1260., 1280., 1300., 1320., 1340.,
&1400., 1800., 2000./
DATA XIMG/1.0E-5, 5.0E-5, 1.1E-3, 1.5E-3, 1.9E-3, 2.3E-3,
&2.4E-3, 2.5E-3, 2.4E-3, 2.2E-3, 1.8E-3, 1.3E-3, 1.1E-3,
&8.6E-4, 6.8E-4, 5.0E-4, 4.0E-4, 2.9E-4, 2.0E-4, 1.8E-4,
&1.2E-4, 5.0E-5, 2.0E-5, 2.1E-5, 2.3E-5, 2.7E-5, 3.1E-5,
&3.5E-5, 5.0E-5, 6.2E-5, 9.0E-5, 1.2E-4, 1.7E-4, 2.0E-4,
&2.4E-4, 3.3E-4, 4.5E-4, 6.2E-4, 9.0E-4, 1.3E-3, 2.0E-3,
&2.9E-3, 4.4E-3, 6.6E-3, 1.0E-2, 1.3E-2, 1.9E-2, 3.0E-2,
&7.0E-2, 1.6E-1, 7.0E-2, 6.0E-3, 6.0E-3, 6.0E-3, 6.0E-3/
DO 100 I=2,55
IF (W .GT. WAVENO(I)) GO TO 100
FACTOR= (WAVENO(I) - W )/(WAVENO(I) - WAVENO(I-1))
REFLIQ=XIMG(I) + FACTOR*(XIMG(I-1) - XIMG(I))
RETURN
100 CONTINUE
REFLIQ=XIMG(55)
RETURN
END
SUBROUTINE ZGRID(NLEVEL,FAC,ZTOP,ZTOP1,Z)
DIMENSION Z(NLEVEL)
C CALCULATE THE GRID SPACING AT BOTTOM
Z(1)=ZTOP
Z(2)=ZTOP1
Z(NLEVEL)=0.0
C SUM UP THE TOTAL EQUIVALENT NUMBER OF BOTTOM GRIDS
C THAT FIT IN THE MODEL
SUM=0.
DO 99 J=NLEVEL,3,-1
SUM=SUM + ((1.-FAC)*(J-3.)/(NLEVEL-3.)+FAC)
99 CONTINUE
C NOTE THAT IF FAC=1 THEN THIS IS JUST NLEVEL-2
DZ=( Z(2)-Z(NLEVEL) )/SUM
DO 101 J=NLEVEL,3,-1
Z(J-1)= Z(J)
& +DZ* ((1.-FAC)*(J-3.)/(NLEVEL-3.)+FAC)
101 CONTINUE
RETURN
END
SUBROUTINE LINDAL(NLEVEL,Z,T,DEN,P)
C THIS SUBROUTINE RETURNS THE LINDAL ET AL INGRESS
C T,P AND DEN PROFILE ON THE INPUT ZGRID.
C RETURNS PRESSURE IN BARS
DIMENSION Z(1),T(1),DEN(1),P(1)
DIMENSION ZLIN(106),TLIN(106),DLIN(106),PLIN(106)
CALL INGRESS(TLIN,ZLIN,DLIN,PLIN)
DO J=1,106
C CONVERT TO BARS
PLIN(J)=PLIN(J)*0.001
ENDDO
DO 100 J=1,NLEVEL
C EXTRAPOLATE WITH ISOTHERMAL ATM ABOVE DATA POINTS
ISTART=1
IF (Z(J) .GT. ZLIN(1) ) THEN
T(J)=TLIN(1)
P(J)= PLIN(1)*EXP(-0.024733*(Z(J)-ZLIN(1)))
DEN(J)=DLIN(1)*EXP(-0.024733*(Z(J)-ZLIN(1)))
ELSE
DO 101 I=ISTART,105
C INTERPOLATE LINEAR IN T EXPONENTIAL IN P AND DEN
IF (Z(J) .LT. ZLIN(I+1) ) GO TO 101
FACTOR= (Z(J)-ZLIN(I) )/(ZLIN(I+1) - ZLIN(I))
T(J)=TLIN(I) + FACTOR*(TLIN(I+1) - TLIN(I))
P(J) =EXP(ALOG(PLIN(I))+FACTOR*(ALOG(PLIN(I+1))-ALOG(PLIN(I))))
DEN(J)=EXP(ALOG(DLIN(I))+FACTOR*(ALOG(DLIN(I+1))-ALOG(DLIN(I))))
ISTART=I
GO TO 100
101 CONTINUE
ENDIF
100 CONTINUE
RETURN
END
SUBROUTINE INGRESS(TEMP,ALT,DEN,PRESS)
DIMENSION TEMP(1),ALT(1),DEN(1),PRESS(1)
DIMENSION TL(106),ZL(106),DENL(106),PL(106)
C THIS SUBROUTINE RETURNS THE DATA FROM LINDAL ET AL (1983) ICARUS 53,
C 348-363, FOR THE INGRESS PROFILE: TEMPERATURE (K), ALTITUDE (KM),
C DENSITY (CM-3), PRESS (MB).
C234567
DATA TL/
&169.4,169.1,168.8,168.6,168.3,168.1,167.9,167.6,167.4,167.3,167.2,
&167.2,167.2,167.1,167.0,166.9,166.9,166.4,165.7,165.4,165.3,165.0,
&164.5,163.9,163.1,162.6,162.3,162.0,161.5,160.9,160.1,159.4,159.0,
&158.4,157.8,157.0,156.1,155.2,154.0,153.0,152.1,151.3,150.7,150.0,
&149.2,148.5,147.8,147.0,145.6,144.3,143.1,141.9,140.6,139.0,137.1,
&135.9,135.1,133.2,130.4,127.4,124.4,121.4,118.4,115.4,110.8,104.8,
& 99.1, 92.1, 85.4, 80.5, 77.5, 75.1, 73.6, 72.7, 72.2, 71.9, 71.6,
& 71.5, 71.4, 71.2, 71.2, 71.4, 71.7, 71.9, 72.3, 72.9, 73.5, 74.2,
& 75.1, 76.2, 77.0, 78.2, 79.5, 80.8, 82.2, 83.6, 85.3, 87.1, 88.1,
& 88.9, 89.9, 91.2, 91.9, 92.6, 93.3, 94.0/
DATA ZL/200.0, 198.0, 196.0, 194.0, 192.0, 190.0,
&188.0, 186.0, 184.0, 182.0, 180.0, 178.0, 176.0, 174.0,
&172.0, 170.0, 168.0, 166.0, 164.0, 162.0, 160.0, 158.0,
&156.0, 154.0, 152.0, 150.0, 148.0, 146.0, 144.0, 142.0,
&140.0, 138.0, 136.0, 134.0, 132.0, 130.0, 128.0, 126.0,
&124.0, 122.0, 120.0, 118.0, 116.0, 114.0, 112.0, 110.0,
&108.0, 106.0, 104.0, 102.0, 100.0, 98.0, 96.0, 94.0,
&92.0, 90.0, 88.0, 86.0, 84.0, 82.0, 80.0, 78.0, 76.0,
&74.0, 72.0, 70.0, 68.0, 66.0, 64.0, 62.0, 60.0, 58.0,
&56.0, 54.0, 52.0, 50.0, 48.0, 46.0, 44.0, 42.0, 40.0,
&38.0, 36.0, 34.0, 32.0, 30.0, 28.0, 26.0, 24.0, 22.0,
&20.0, 18.0, 16.0, 14.0, 12.0, 10.0, 8.0, 6.0, 5.0, 4.0,
&3.0, 2.0, 1.5, 1.0, 0.5, 0.0/
DATA DENL/0.32, 0.34, 0.35, 0.37, 0.39, 0.41, 0.43, 0.45,
&0.47, 0.50, 0.52, 0.55, 0.57, 0.60, 0.63, 0.66, 0.70,
&0.73, 0.77, 0.81, 0.85, 0.90, 0.95, 1.00, 1.06, 1.11,
&1.17, 1.24, 1.30, 1.38, 1.46, 1.54, 1.63, 1.72, 1.82,
&1.92, 2.04, 2.17, 2.30, 2.45, 2.60, 2.76, 2.93, 3.11,
&3.31, 3.52, 3.74, 3.98, 4.26, 4.55, 4.87, 5.21, 5.58,
&6.00, 6.47, 6.95, 7.44, 8.05, 8.77, 9.60, 10.51, 11.57,
&12.74, 14.07, 15.83, 18.14, 20.87, 24.61, 29.29, 34.51,
&40.07, 46.34, 53.23, 60.73, 69.07, 78.48, 89.07, 101.09,
&114.65, 130.33, 147.77, 167.00, 188.70, 213.27, 240.49,
&270.56, 304.10, 340.67, 380.79, 423.99, 473.54, 525.79,
&582.13, 643.58, 710.23, 781.93, 856.89, 937.22, 979.00,
&1023.97, 1067.09, 1108.70, 1129.37, 1150.16, 1171.26,
&1192.38/
DATA PL/
& 0.75, 0.79, 0.83, 0.87, 0.91, 0.95, 1.00, 1.04,
& 1.10, 1.15, 1.20, 1.26, 1.33, 1.39, 1.46, 1.53,
& 1.61, 1.69, 1.77, 1.86, 1.95, 2.05, 2.15, 2.26,
& 2.38, 2.50, 2.63, 2.76, 2.91, 3.06, 3.22, 3.39,
& 3.57, 3.76, 3.96, 4.17, 4.40, 4.64, 4.89, 5.17,
& 5.46, 5.76, 6.09, 6.44, 6.81, 7.21, 7.63, 8.08,
& 8.56, 9.07, 9.62, 10.21, 10.84, 11.52, 12.25, 13.03,
& 13.88, 14.80, 15.79, 16.88, 18.06, 19.38, 20.81, 22.41,
& 24.20, 26.21, 28.54, 31.25, 34.48, 38.31, 42.78, 47.97,
& 53.97, 60.84, 68.67, 77.60, 87.75, 99.29, 112.40, 127.30,
& 144.26, 163.48, 185.22, 209.84, 237.68, 269.08, 304.45, 344.16,
& 388.68, 438.44, 494.04, 556.03, 624.87, 701.09, 785.43, 878.54,
& 980.98,1093.37,1153.44,1216.28,1282.01,1350.40,1385.60,1421.48,
&1458.02,1495.26/
DO 10 J=1,106
TEMP(J)=TL(J)
ALT(J)=ZL(J)
PRESS(J)=PL(J)
10 DEN(J)=DENL(J)*1.E17
RETURN
END
FUNCTION EFFG(Z)
DATA G/1.35/
DATA RPLANT/2575./
EFFG = G * (RPLANT/(RPLANT + Z ) )**2
RETURN
END
SUBROUTINE DSOLVER(NL,GAMA,CP,CM,CPM1,CMM1
,,E1,E2,E3,E4,BTOP,BSURF,RSF,XK1,XK2)
C DOUBLE PRECISION VERSION OF SOLVER
PARAMETER (NMAX=201)
IMPLICIT REAL*8 (A-H,O-Z)
DIMENSION GAMA(NL),CP(NL),CM(NL),
,CPM1(NL),CMM1(NL),XK1(NL),XK2(NL)
,,E1(NL),E2(NL),E3(NL),E4(NL)
DIMENSION AF(NMAX),BF(NMAX),CF(NMAX),DF(NMAX),XK(NMAX)
C*********************************************************
C* THIS SUBROUTINE SOLVES FOR THE COEFFICIENTS OF THE *
C* TWO STREAM SOLUTION FOR GENERAL BOUNDARY CONDITIONS *
C* NO ASSUMPTION OF THE DEPENDENCE ON OPTICAL DEPTH OF *
C* C-PLUS OR C-MINUS HAS BEEN MADE. *
C* NL = NUMBER OF LAYERS IN THE MODEL *
C* CP = C-PLUS EVALUATED AT TAO=0 (TOP) *
C* CM = C-MINUS EVALUATED AT TAO=0 (TOP) *
C* CPM1 = C-PLUS EVALUATED AT TAOSTAR (BOTTOM) *
C* CMM1 = C-MINUS EVALUATED AT TAOSTAR (BOTTOM) *
C* EP = EXP(LAMDA*DTAU) *
C* EM = 1/EP *
C* E1 = EP + GAMA *EM *
C* E2 = EP - GAMA *EM *
C* E3 = GAMA*EP + EM *
C* E4 = GAMA*EP - EM *
C* BTOP = THE DIFFUSE RADIATION INTO THE MODEL AT TOP *
C* BSURF = THE DIFFUSE RADIATION INTO THE MODEL AT *
C* THE BOTTOM: INCLUDES EMMISION AND REFLECTION *
C* OF THE UNATTENUATED PORTION OF THE DIRECT *
C* BEAM. BSTAR+RSF*FO*EXP(-TAOSTAR/U0) *
C* RSF = REFLECTIVITY OF THE SURFACE *
C* XK1 = COEFFICIENT OF THE POSITIVE EXP TERM *
C* XK2 = COEFFICIENT OF THE NEGATIVE EXP TERM *
C*********************************************************
L=2*NL
C************MIXED COEFFICENTS**********
C* THIS VERSION AVOIDS SINGULARITIES ASSOC.
C* WIRH W0=0 BY SOLVING FOR XK1+XK2, AND XK1-XK2.
AF(1)=0.0
BF(1)=GAMA(1)+1.
CF(1)=GAMA(1)-1.
DF(1)=BTOP-CMM1(1)
N=0
LM2=L-2
C* EVEN TERMS
DO 10 I=2,LM2,2
N=N+1
AF(I)=(E1(N)+E3(N))*(GAMA(N+1)-1.)
BF(I)=(E2(N)+E4(N))*(GAMA(N+1)-1.)
CF(I)=2.*(1.-GAMA(N+1)**2)
DF(I)=(GAMA(N+1)-1.) * (CPM1(N+1) - CP(N))
& + (1.-GAMA(N+1))* (CM(N)-CMM1(N+1))
10 CONTINUE
N=0
LM1=L-1
DO 20 I=3,LM1,2
N=N+1
AF(I)=2.*(1.-GAMA(N)**2)
BF(I)=(E1(N)-E3(N))*(1.+GAMA(N+1))
CF(I)=(E1(N)+E3(N))*(GAMA(N+1)-1.)
DF(I)=E3(N)*(CPM1(N+1) - CP(N))
& + E1(N)*(CM(N) - CMM1(N+1))
20 CONTINUE
AF(L)=E1(NL)-RSF*E3(NL)
BF(L)=E2(NL)-RSF*E4(NL)
CF(L)=0.0
DF(L)=BSURF-CP(NL)+RSF*CM(NL)
CALL DTRIDGL(L,AF,BF,CF,DF,XK)
C***UNMIX THE COEFFICIENTS****
DO 28 N=1,NL
XK1(N)=XK(2*N-1)+XK(2*N)
XK2(N)=XK(2*N-1)-XK(2*N)
C NOW TEST TO SEE IF XK2 IS REALLY ZERO TO THE LIMIT OF THE
C MACHINE ACCURACY = 1 .E -30
C XK2 IS THE COEFFICEINT OF THE GROWING EXPONENTIAL AND MUST
C BE TREATED CAREFULLY
IF (XK2(N) .EQ. 0.0) GO TO 28
IF (ABS (XK2(N)/XK(2*N-1)) .LT. 1.E-30) XK2(N)=0.0
28 CONTINUE
RETURN
END
SUBROUTINE DTRIDGL(L,AF,BF,CF,DF,XK)
C DOUBLE PRESCISION VERSION OF TRIDGL
PARAMETER (NMAX=201)
IMPLICIT REAL*8 (A-H,O-Z)
DIMENSION AF(L),BF(L),CF(L),DF(L),XK(L)
DIMENSION AS(NMAX),DS(NMAX)
C* THIS SUBROUTINE SOLVES A SYSTEM OF TRIDIAGIONAL MATRIX
C* EQUATIONS. THE FORM OF THE EQUATIONS ARE:
C* A(I)*X(I-1) + B(I)*X(I) + C(I)*X(I+1) = D(I)
C* WHERE I=1,L LESS THAN 103.
C* ..............REVIEWED -CP........
AS(L) = AF(L)/BF(L)
DS(L) = DF(L)/BF(L)
DO 10 I=2,L
X=1./(BF(L+1-I) - CF(L+1-I)*AS(L+2-I))
AS(L+1-I)=AF(L+1-I)*X
DS(L+1-I)=(DF(L+1-I)-CF(L+1-I)*DS(L+2-I))*X
10 CONTINUE
XK(1)=DS(1)
DO 20 I=2,L
XKB=XK(I-1)
XK(I)=DS(I)-AS(I)*XKB
20 CONTINUE
910 FORMAT(/,8X,'AF(I)',7X,'BF(I)',7X,'CF(I)',7X,'DF(I)',7X,
* 'AS(I)',7X,'DS(I)',7X,'XK(I)',/)
915 FORMAT((3X,7(1X,1PE11.4)))
RETURN
END
FUNCTION PLANCK(WAV,T,DW)
C* PLANCK FUNCTION RETURNS B IN CGS UNITS, ERGS CM-2 WAVENUMBER-1
C* WAVNUM IS WAVENUMBER IN CM-1
C* T IS IN KELVIN
DATA NBB/2/
PLANCK=0.
DO 100 I=-NBB,NBB,1
WAVNUM=WAV + I*DW/(2.*NBB)
PLANCK=PLANCK+1.191E-5* WAVNUM**3 * EXP(-1.4388 * WAVNUM/T )/
& ( 1.0- EXP(-1.4388 * WAVNUM/T ) )
100 CONTINUE
PLANCK=PLANCK/(2*NBB+1)
RETURN
END
SUBROUTINE SOLVER(NL,GAMA,CP,CM,CPM1,CMM1
,,E1,E2,E3,E4,BTOP,BSURF,RSF,XK1,XK2)
PARAMETER (NMAX=201)
DIMENSION GAMA(NL),CP(NL),CM(NL),
,CPM1(NL),CMM1(NL),XK1(NL),XK2(NL)
,,E1(NL),E2(NL),E3(NL),E4(NL)
DIMENSION AF(NMAX),BF(NMAX),CF(NMAX),DF(NMAX),XK(NMAX)
C*********************************************************
C* THIS SUBROUTINE SOLVES FOR THE COEFFICIENTS OF THE *
C* TWO STREAM SOLUTION FOR GENERAL BOUNDARY CONDITIONS *
C* NO ASSUMPTION OF THE DEPENDENCE ON OPTICAL DEPTH OF *
C* C-PLUS OR C-MINUS HAS BEEN MADE. *
C* NL = NUMBER OF LAYERS IN THE MODEL *
C* CP = C-PLUS EVALUATED AT TAO=0 (TOP) *
C* CM = C-MINUS EVALUATED AT TAO=0 (TOP) *
C* CPM1 = C-PLUS EVALUATED AT TAOSTAR (BOTTOM) *
C* CMM1 = C-MINUS EVALUATED AT TAOSTAR (BOTTOM) *
C* EP = EXP(LAMDA*DTAU) *
C* EM = 1/EP *
C* E1 = EP + GAMA *EM *
C* E2 = EP - GAMA *EM *
C* E3 = GAMA*EP + EM *
C* E4 = GAMA*EP - EM *
C* BTOP = THE DIFFUSE RADIATION INTO THE MODEL AT TOP *
C* BSURF = THE DIFFUSE RADIATION INTO THE MODEL AT *
C* THE BOTTOM: INCLUDES EMMISION AND REFLECTION *
C* OF THE UNATTENUATED PORTION OF THE DIRECT *
C* BEAM. BSTAR+RSF*FO*EXP(-TAOSTAR/U0) *
C* RSF = REFLECTIVITY OF THE SURFACE *
C* XK1 = COEFFICIENT OF THE POSITIVE EXP TERM *
C* XK2 = COEFFICIENT OF THE NEGATIVE EXP TERM *
C*********************************************************
L=2*NL
C************MIXED COEFFICENTS**********
C* THIS VERSION AVOIDS SINGULARITIES ASSOC.
C* WIRH W0=0 BY SOLVING FOR XK1+XK2, AND XK1-XK2.
AF(1)=0.0
BF(1)=GAMA(1)+1.
CF(1)=GAMA(1)-1.
DF(1)=BTOP-CMM1(1)
N=0
LM2=L-2
C* EVEN TERMS
DO 10 I=2,LM2,2
N=N+1
AF(I)=(E1(N)+E3(N))*(GAMA(N+1)-1.)
BF(I)=(E2(N)+E4(N))*(GAMA(N+1)-1.)
CF(I)=2.*(1.-GAMA(N+1)**2)
DF(I)=(GAMA(N+1)-1.) * (CPM1(N+1) - CP(N))
& + (1.-GAMA(N+1))* (CM(N)-CMM1(N+1))
10 CONTINUE
N=0
LM1=L-1
DO 20 I=3,LM1,2
N=N+1
AF(I)=2.*(1.-GAMA(N)**2)
BF(I)=(E1(N)-E3(N))*(1.+GAMA(N+1))
CF(I)=(E1(N)+E3(N))*(GAMA(N+1)-1.)
DF(I)=E3(N)*(CPM1(N+1) - CP(N))
& + E1(N)*(CM(N) - CMM1(N+1))
20 CONTINUE
AF(L)=E1(NL)-RSF*E3(NL)
BF(L)=E2(NL)-RSF*E4(NL)
CF(L)=0.0
DF(L)=BSURF-CP(NL)+RSF*CM(NL)
CALL TRIDGL(L,AF,BF,CF,DF,XK)
C***UNMIX THE COEFFICIENTS****
DO 28 N=1,NL
XK1(N)=XK(2*N-1)+XK(2*N)
XK2(N)=XK(2*N-1)-XK(2*N)
28 CONTINUE
RETURN
END
SUBROUTINE TRIDGL(L,AF,BF,CF,DF,XK)
PARAMETER (NMAX=201)
DIMENSION AF(L),BF(L),CF(L),DF(L),XK(L)
DIMENSION AS(NMAX),DS(NMAX)
C* THIS SUBROUTINE SOLVES A SYSTEM OF TRIDIAGIONAL MATRIX
C* EQUATIONS. THE FORM OF THE EQUATIONS ARE:
C* A(I)*X(I-1) + B(I)*X(I) + C(I)*X(I+1) = D(I)
C* WHERE I=1,L LESS THAN 103.
C* ..............REVIEWED -CP........
AS(L) = AF(L)/BF(L)
DS(L) = DF(L)/BF(L)
DO 10 I=2,L
X=1./(BF(L+1-I) - CF(L+1-I)*AS(L+2-I))
AS(L+1-I)=AF(L+1-I)*X
DS(L+1-I)=(DF(L+1-I)-CF(L+1-I)*DS(L+2-I))*X
10 CONTINUE
XK(1)=DS(1)
DO 20 I=2,L
XKB=XK(I-1)
XK(I)=DS(I)-AS(I)*XKB
20 CONTINUE
910 FORMAT(/,8X,'AF(I)',7X,'BF(I)',7X,'CF(I)',7X,'DF(I)',7X,
* 'AS(I)',7X,'DS(I)',7X,'XK(I)',/)
915 FORMAT((3X,7(1X,1PE11.4)))
RETURN
END
SUBROUTINE DMIESS( RO, RFR, RFI, THETD, JX, 2
2 QEXT, QSCAT, CTBRQS, ELTRMX, PI, 3
3 TAU, CSTHT, SI2THT, ACAP, IT, 4
4 LL, R, RE2, TMAG2, WVNO ) 5
C 6
C THIS VERSION OF THE FAMOUS DMIESS CODE IS ADAPTED FOR THE VAX
C ITS ACCURACY HAS BEEN CHECKED BY COMPARISON TO CRAY RUNS.
C -C.P. MCKAY
C
IMPLICIT REAL*8 (A-H,O-Z)
C 13
COMPLEX*16 FNAP, FNBP, ACAP(LL), W, 14
2 FNA, FNB, RF, RRF, 15
3 RRFX, WM1, FN1, FN2, 16
4 TC1, TC2, WFN(2), Z(4), 17
5 K1, K2, K3, SIN, 18
6 COS, RC, U(8), DH1, 19
7 DH2, DH4, P24H24, P24H21, 20
8 PSTORE, HSTORE, DUMMY, DUMSQ 21
C 22
DIMENSION W(3,9000) 12
DIMENSION T(5), TA(4), TB(2), TC(2), 23
2 TD(2), TE(2), PI( 3,IT ), TAU( 3,IT ), 24
3 CSTHT(IT), THETD(IT), SI2THT(IT), ELTRMX( 4,IT,2 ) 25
C 26
EQUIVALENCE ( WFN(1),TA(1) ), ( FNA,TB(1) ), 27
2 ( FNB,TC(1) ), ( FNAP,TD(1) ), 28
3 ( FNBP,TE(1) ) 29
C 30
C 31
C THIS SUBROUTINE COMPUTES MIE SCATTERING BY A STRATIFIED SPHERE, 32
C I.E. A PARTICLE CONSISTING OF A SPHERICAL CORE SURROUNDED BY A 33
C SPHERICAL SHELL. THE BASIC CODE USED WAS THAT DESCRIBED IN THE 34
C REPORT: " SUBROUTINES FOR COMPUTING THE PARAMETERS OF THE 35
C ELECTROMAGNETIC RADIATION SCATTERED BY A SPHERE " J.V. DAVE, 36
C I B M SCIENTIFIC CENTER, PALO ALTO , CALIFORNIA. 37
C REPORT NO. 320 - 3236 .. MAY 1968 . 38
C 39
C THE MODIFICATIONS FOR STRATIFIED SPHERES ARE DESCRIBED IN 40
C TOON AND ACKERMAN, APPL. OPTICS, IN PRESS, 1981 41
C 42
C 43
C THE PARAMETERS IN THE CALLING STATEMENT ARE DEFINED AS FOLLOWS : 44
C RO IS THE OUTER (SHELL) RADIUS; 45
C R IS THE CORE RADIUS; 46
C RFR, RFI ARE THE REAL AND IMAGINARY PARTS OF THE SHELL INDEX 47
C OF REFRACTION IN THE FORM (RFR - I* RFI); 48
C RE2, TMAG2 ARE THE INDEX PARTS FOR THE CORE; 49
C ( WE ASSUME SPACE HAS UNIT INDEX. ) 50
C THETD(J): ANGLE IN DEGREES BETWEEN THE DIRECTIONS OF THE INCIDENT 51
C AND THE SCATTERED RADIATION. THETD(J) IS< OR= 90.0 52
C IF THETD(J) SHOULD HAPPEN TO BE GREATER THAN 90.0, ENTER WITH 53
C SUPPLEMENTARY VALUE, SEE COMMENTS BELOW ON ELTRMX; 54
C JX: TOTAL NUMBER OF THETD FOR WHICH THE COMPUTATIONS ARE 55
C REQUIRED. JX SHOULD NOT EXCEED IT UNLESS THE DIMENSIONS 56
C STATEMENTS ARE APPROPRIATEDLY MODIFIED; 57
C 58
C THE DEFINITIONS FOR THE FOLLOWING SYMBOLS CAN BE FOUND IN"LIGHT 59
C SCATTERING BY SMALL PARTICLES,H.C.VAN DE HULST, JOHN WILEY ' 60
C SONS, INC., NEW YORK, 1957" . 61
C QEXT: EFFIECIENCY FACTOR FOR EXTINCTION,VAN DE HULST,P.14 ' 127. 62
C QSCAT: EFFIECINCY FACTOR FOR SCATTERING,V.D. HULST,P.14 ' 127. 63
C CTBRQS: AVERAGE(COSINE THETA) * QSCAT,VAN DE HULST,P.128 64
C ELTRMX(I,J,K): ELEMENTS OF THE TRANSFORMATION MATRIX F,V.D.HULST 65
C ,P.34,45 ' 125. I=1: ELEMENT M SUB 2..I=2: ELEMENT M SUB 1.. 66
C I = 3: ELEMENT S SUB 21.. I = 4: ELEMENT D SUB 21.. 67
C ELTRMX(I,J,1) REPRESENTS THE ITH ELEMENT OF THE MATRIX FOR 68
C THE ANGLE THETD(J).. ELTRMX(I,J,2) REPRESENTS THE ITH ELEMENT 69
C OF THE MATRIX FOR THE ANGLE 180.0 - THETD(J) .. 70
C 71
C IT: IS THE DIMENSION OF THETD, ELTRMX, CSTHT, PI, TAU, SI2THT, 72
C IT MUST CORRESPOND EXACTLY TO THE SECOND DIMENSION OF ELTRMX. 73
C LL: IS THE DIMENSION OF ACAP 74
C IN THE ORIGINAL PROGRAM THE DIMENSION OF ACAP WAS 7000. 75
C FOR CONSERVING SPACE THIS SHOULD BE NOT MUCH HIGHER THAN 76
C THE VALUE, N=1.1*(NREAL**2 + NIMAG**2)**.5 * X + 1 77
C WVNO: 2*PI / WAVELENGTH 78
C 79
C THIS SUBROUTINE COMPUTES THE CAPITAL A FUNCTION BY MAKING USE OF 80
C DOWNWARD RECURRENCE RELATIONSHIP. 81
C 82
C TA(1): REAL PART OF WFN(1). TA(2): IMAGINARY PART OF WFN(1). 83
C TA(3): REAL PART OF WFN(2). TA(4): IMAGINARY PART OF WFN(2). 84
C TB(1): REAL PART OF FNA. TB(2): IMAGINARY PART OF FNA. 85
C TC(1): REAL PART OF FNB. TC(2): IMAGINARY PART OF FNB. 86
C TD(1): REAL PART OF FNAP. TD(2): IMAGINARY PART OF FNAP. 87
C TE(1): REAL PART OF FNBP. TE(2): IMAGINARY PART OF FNBP. 88
C FNAP, FNBP ARE THE PRECEDING VALUES OF FNA, FNB RESPECTIVELY. 89
C 90
C IF THE CORE IS SMALL SCATTERING IS COMPUTED FOR THE SHELL ONLY 91
IFLAG = 1 92
IF ( R/RO .LT. 1.0E-06 ) IFLAG = 2 93
C 94
IF ( JX .LE. IT ) GO TO 20 95
WRITE( 6,7 ) 96
WRITE( 6,6 ) 97
STOP 30 98
20 RF = DCMPLX( RFR, -RFI ) 99
RC = DCMPLX( RE2,-TMAG2 ) 100
X = RO * WVNO 101
K1 = RC * WVNO 102
K2 = RF * WVNO
ZET=0.0 103
K3 = DCMPLX( WVNO, ZET ) 104
Z(1) = K2 * RO 105
Z(2) = K3 * RO 106
Z(3) = K1 * R 107
Z(4) = K2 * R 108
X1 = DREAL( Z(1) ) 109
Y1 = DIMAG( Z(1) ) 110
X4 = DREAL( Z(4) ) 111
Y4 = DIMAG( Z(4) ) 112
RRF = 1.0 / RF 113
RX = 1.0 / X 114
RRFX = RRF * RX 115
T(1) = ( X**2 ) * ( RFR**2 + RFI**2 ) 116
T(1) = SQRT( T(1) ) 117
NMX1 = 1.10 * T(1) 118
C 119
IF ( NMX1 .LE. LL-1 ) GO TO 21 120
WRITE(6,8) 121
STOP 32 122
21 NMX2 = T(1) 123
IF ( NMX1 .GT. 150 ) GO TO 22 124
NMX1 = 150 125
NMX2 = 135 126
C 127
22 ACAP( NMX1+1 ) = ( 0.0,0.0 ) 128
IF ( IFLAG .EQ. 2 ) GO TO 26 129
DO 29 N = 1,3 130
29 W( N,NMX1+1 ) = ( 0.0,0.0 ) 131
26 CONTINUE 132
DO 23 N = 1,NMX1 133
NN = NMX1 - N + 1 134
ACAP(NN) = (NN+1) * RRFX - 1.0 / ( (NN+1) * RRFX + ACAP(NN+1) ) 135
IF ( IFLAG .EQ. 2 ) GO TO 23 136
DO 31 M = 1,3 137
31 W( M,NN ) = (NN+1) / Z(M+1) - 138
1 1.0 / ( (NN+1) / Z(M+1) + W( M,NN+1 ) ) 139
23 CONTINUE 140
DO 30 J = 1,JX 141
IF ( THETD(J) .LT. 0.0 ) THETD(J) = ABS( THETD(J) ) 142
IF ( THETD(J) .GT. 0.0 ) GO TO 24 143
CSTHT(J) = 1.0 144
SI2THT(J) = 0.0 145
GO TO 30 146
24 IF ( THETD(J) .GE. 90.0 ) GO TO 25 147
T(1) = ( 3.14159265359 * THETD(J) ) / 180.0 148
CSTHT(J) = COS( T(1) ) 149
SI2THT(J) = 1.0 - CSTHT(J)**2 150
GO TO 30 151
25 IF ( THETD(J) .GT. 90.0 ) GO TO 28 152
CSTHT(J) = 0.0 153
SI2THT(J) = 1.0 154
GO TO 30 155
28 WRITE( 6,5 ) THETD(J) 156
WRITE( 6,6 ) 157
STOP 34 158
30 CONTINUE 159
C 160
DO 35 J = 1,JX 161
PI(1,J) = 0.0 162
PI(2,J) = 1.0 163
TAU(1,J) = 0.0 164
TAU(2,J) = CSTHT(J) 165
35 CONTINUE 166
C 167
C INITIALIZATION OF HOMOGENEOUS SPHERE 168
T(1) = COS(X) 169
T(2) = SIN(X) 170
WM1 = DCMPLX( T(1),-T(2) ) 171
WFN(1) = DCMPLX( T(2), T(1) ) 172
WFN(2) = RX * WFN(1) - WM1 173
C 174
IF ( IFLAG .EQ. 2 ) GO TO 560 175
N = 1 176
C 177
C INITIALIZATION PROCEDURE FOR STRATIFIED SPHERE BEGINS HERE 178
C 179
SINX1 = SIN( X1 ) 180
SINX4 = SIN( X4 ) 181
COSX1 = COS( X1 ) 182
COSX4 = COS( X4 ) 183
EY1 = EXP( Y1 ) 184
E2Y1 = EY1 * EY1 185
EY4 = EXP( Y4 ) 186
EY1MY4 = EXP( Y1 - Y4 ) 187
EY1PY4 = EY1 * EY4 188
EY1MY4 = EXP( Y1 - Y4 ) 189
AA = SINX4 * ( EY1PY4 + EY1MY4 ) 190
BB = COSX4 * ( EY1PY4 - EY1MY4 ) 191
CC = SINX1 * ( E2Y1 + 1.0 ) 192
DD = COSX1 * ( E2Y1 - 1.0 ) 193
DENOM = 1.0 + E2Y1 * ( 4.0 * SINX1 * SINX1 - 2.0 + E2Y1 ) 194
REALP = ( AA * CC + BB * DD ) / DENOM 195
AMAGP = ( BB * CC - AA * DD ) / DENOM 196
DUMMY = DCMPLX( REALP, AMAGP ) 197
AA = SINX4 * SINX4 - 0.5 198
BB = COSX4 * SINX4 199
P24H24 = 0.5 + DCMPLX( AA,BB ) * EY4 * EY4 200
AA = SINX1 * SINX4 - COSX1 * COSX4 201
BB = SINX1 * COSX4 + COSX1 * SINX4 202
CC = SINX1 * SINX4 + COSX1 * COSX4 203
DD = -SINX1 * COSX4 + COSX1 * SINX4 204
P24H21 = 0.5 * DCMPLX( AA,BB ) * EY1 * EY4 + 205
2 0.5 * DCMPLX( CC,DD ) * EY1MY4 206
DH4 = Z(4) / ( 1.0 + ( 0.0,1.0 ) * Z(4) ) - 1.0 / Z(4) 207
DH1 = Z(1) / ( 1.0 + ( 0.0,1.0 ) * Z(1) ) - 1.0 / Z(1) 208
DH2 = Z(2) / ( 1.0 + ( 0.0,1.0 ) * Z(2) ) - 1.0 / Z(2) 209
PSTORE = ( DH4 + N / Z(4) ) * ( W(3,N) + N / Z(4) ) 210
P24H24 = P24H24 / PSTORE 211
HSTORE = ( DH1 + N / Z(1) ) * ( W(3,N) + N / Z(4) ) 212
P24H21 = P24H21 / HSTORE 213
PSTORE = ( ACAP(N) + N / Z(1) ) / ( W(3,N) + N / Z(4) ) 214
DUMMY = DUMMY * PSTORE 215
DUMSQ = DUMMY * DUMMY 216
C 217
C NOTE: THE DEFINITIONS OF U(I) IN THIS PROGRAM ARE NOT THE SAME AS 218
C THE USUBI DEFINED IN THE ARTICLE BY TOON AND ACKERMAN. THE 219
C CORRESPONDING TERMS ARE: 220
C USUB1 = U(1) USUB2 = U(5) 221
C USUB3 = U(7) USUB4 = DUMSQ 222
C USUB5 = U(2) USUB6 = U(3) 223
C USUB7 = U(6) USUB8 = U(4) 224
C RATIO OF SPHERICAL BESSEL FTN TO SPHERICAL HENKAL FTN = U(8) 225
C 226
U(1) = K3 * ACAP(N) - K2 * W(1,N) 227
U(2) = K3 * ACAP(N) - K2 * DH2 228
U(3) = K2 * ACAP(N) - K3 * W(1,N) 229
U(4) = K2 * ACAP(N) - K3 * DH2 230
U(5) = K1 * W(3,N) - K2 * W(2,N) 231
U(6) = K2 * W(3,N) - K1 * W(2,N) 232
U(7) = ( 0.0,-1.0 ) * ( DUMMY * P24H21 - P24H24 ) 233
U(8) = TA(3) / WFN(2) 234
C 235
FNA = U(8) * ( U(1)*U(5)*U(7) + K1*U(1) - DUMSQ*K3*U(5) ) / 236
2 ( U(2)*U(5)*U(7) + K1*U(2) - DUMSQ*K3*U(5) ) 237
FNB = U(8) * ( U(3)*U(6)*U(7) + K2*U(3) - DUMSQ*K2*U(6) ) / 238
2 ( U(4)*U(6)*U(7) + K2*U(4) - DUMSQ*K2*U(6) ) 239
GO TO 561 240
560 TC1 = ACAP(1) * RRF + RX 241
TC2 = ACAP(1) * RF + RX 242
FNA = ( TC1 * TA(3) - TA(1) ) / ( TC1 * WFN(2) - WFN(1) ) 243
FNB = ( TC2 * TA(3) - TA(1) ) / ( TC2 * WFN(2) - WFN(1) ) 244
561 CONTINUE 245
FNAP = FNA 246
FNBP = FNB 247
T(1) = 1.50 248
C 249
C FROM HERE TO THE STATMENT NUMBER 90, ELTRMX(I,J,K) HAS 250
C FOLLOWING MEANING: 251
C ELTRMX(1,J,K): REAL PART OF THE FIRST COMPLEX AMPLITUDE. 252
C ELTRMX(2,J,K): IMAGINARY PART OF THE FIRST COMPLEX AMPLITUDE. 253
C ELTRMX(3,J,K): REAL PART OF THE SECOND COMPLEX AMPLITUDE. 254
C ELTRMX(4,J,K): IMAGINARY PART OF THE SECOND COMPLEX AMPLITUDE. 255
C K = 1 : FOR THETD(J) AND K = 2 : FOR 180.0 - THETD(J) 256
C DEFINITION OF THE COMPLEX AMPLITUDE: VAN DE HULST,P.125. 257
TB(1) = T(1) * TB(1) 258
TB(2) = T(1) * TB(2) 259
TC(1) = T(1) * TC(1) 260
TC(2) = T(1) * TC(2) 261
DO 60 J = 1,JX 262
ELTRMX(1,J,1) = TB(1) * PI(2,J) + TC(1) * TAU(2,J) 263
ELTRMX(2,J,1) = TB(2) * PI(2,J) + TC(2) * TAU(2,J) 264
ELTRMX(3,J,1) = TC(1) * PI(2,J) + TB(1) * TAU(2,J) 265
ELTRMX(4,J,1) = TC(2) * PI(2,J) + TB(2) * TAU(2,J) 266
ELTRMX(1,J,2) = TB(1) * PI(2,J) - TC(1) * TAU(2,J) 267
ELTRMX(2,J,2) = TB(2) * PI(2,J) - TC(2) * TAU(2,J) 268
ELTRMX(3,J,2) = TC(1) * PI(2,J) - TB(1) * TAU(2,J) 269
ELTRMX(4,J,2) = TC(2) * PI(2,J) - TB(2) * TAU(2,J) 270
60 CONTINUE 271
C 272
QEXT = 2.0 * ( TB(1) + TC(1)) 273
QSCAT = ( TB(1)**2 + TB(2)**2 + TC(1)**2 + TC(2)**2 ) / 0.75 274
CTBRQS = 0.0 275
N = 2 276
65 T(1) = 2*N - 1 277
T(2) = N - 1 278
T(3) = 2*N + 1 279
DO 70 J = 1,JX 280
PI(3,J) = ( T(1) * PI(2,J) * CSTHT(J) - N * PI(1,J) ) / T(2) 281
TAU(3,J) = CSTHT(J) * ( PI(3,J) - PI(1,J) ) - 282
1 T(1) * SI2THT(J) * PI(2,J) + TAU(1,J) 283
70 CONTINUE 284
C 285
C HERE SET UP HOMOGENEOUS SPHERE 286
WM1 = WFN(1) 287
WFN(1) = WFN(2) 288
WFN(2) = T(1) * RX * WFN(1) - WM1 289
IF ( IFLAG .EQ. 2 ) GO TO 1000 290
C 291
C HERE SET UP STRATIFIED SPHERE 292
C 293
DH2 = - N / Z(2) + 1.0 / ( N / Z(2) - DH2 ) 294
DH4 = - N / Z(4) + 1.0 / ( N / Z(4) - DH4 ) 295
DH1 = - N / Z(1) + 1.0 / ( N / Z(1) - DH1 ) 296
PSTORE = ( DH4 + N / Z(4) ) * ( W(3,N) + N / Z(4) ) 297
P24H24 = P24H24 / PSTORE 298
HSTORE = ( DH1 + N / Z(1) ) * ( W(3,N) + N / Z(4) ) 299
P24H21 = P24H21 / HSTORE 300
PSTORE = ( ACAP(N) + N / Z(1) ) / ( W(3,N) + N / Z(4) ) 301
DUMMY = DUMMY * PSTORE 302
DUMSQ = DUMMY * DUMMY 303
C 304
U(1) = K3 * ACAP(N) - K2 * W(1,N) 305
U(2) = K3 * ACAP(N) - K2 * DH2 306
U(3) = K2 * ACAP(N) - K3 * W(1,N) 307
U(4) = K2 * ACAP(N) - K3 * DH2 308
U(5) = K1 * W(3,N) - K2 * W(2,N) 309
U(6) = K2 * W(3,N) - K1 * W(2,N) 310
U(7) = ( 0.0,-1.0 ) * ( DUMMY * P24H21 - P24H24 ) 311
U(8) = TA(3) / WFN(2) 312
C 313
FNA = U(8) * ( U(1)*U(5)*U(7) + K1*U(1) - DUMSQ*K3*U(5) ) / 314
2 ( U(2)*U(5)*U(7) + K1*U(2) - DUMSQ*K3*U(5) ) 315
FNB = U(8) * ( U(3)*U(6)*U(7) + K2*U(3) - DUMSQ*K2*U(6) ) / 316
2 ( U(4)*U(6)*U(7) + K2*U(4) - DUMSQ*K2*U(6) ) 317
C 318
1000 CONTINUE 319
TC1 = ACAP(N) * RRF + N * RX 320
TC2 = ACAP(N) * RF + N * RX 321
FN1 = ( TC1 * TA(3) - TA(1) ) / ( TC1 * WFN(2) - WFN(1) ) 322
FN2 = ( TC2 * TA(3) - TA(1) ) / ( TC2 * WFN(2) - WFN(1) ) 323
M = WVNO * R 324
IF ( N .LT. M ) GO TO 1002 325
IF ( IFLAG .EQ. 2 ) GO TO 1001 326
IF ( ABS( ( FN1-FNA ) / FN1 ) .LT. 1.0E-09 .AND. 327
1 ABS( ( FN2-FNB ) / FN2 ) .LT . 1.0E-09 ) IFLAG = 2 328
IF ( IFLAG .EQ. 1 ) GO TO 1002 329
1001 FNA = FN1 330
FNB = FN2 331
1002 CONTINUE 332
T(5) = N 333
T(4) = T(1) / ( T(5) * T(2) ) 334
T(2) = ( T(2) * ( T(5) + 1.0 ) ) / T(5) 335
C 336
CTBRQS = CTBRQS + T(2) * ( TD(1) * TB(1) + TD(2) * TB(2) + 337
1 TE(1) * TC(1) + TE(2) * TC(2) ) 338
2 + T(4) * ( TD(1) * TE(1) + TD(2) * TE(2) ) 339
QEXT = QEXT + T(3) * ( TB(1) + TC(1) ) 340
T(4) = TB(1)**2 + TB(2)**2 + TC(1)**2 + TC(2)**2 341
QSCAT = QSCAT + T(3) * T(4) 342
T(2) = N * (N+1) 343
T(1) = T(3) / T(2) 344
K = (N/2)*2 345
DO 80 J = 1,JX 346
ELTRMX(1,J,1) = ELTRMX(1,J,1)+T(1)*(TB(1)*PI(3,J)+TC(1)*TAU(3,J)) 347
ELTRMX(2,J,1) = ELTRMX(2,J,1)+T(1)*(TB(2)*PI(3,J)+TC(2)*TAU(3,J)) 348
ELTRMX(3,J,1) = ELTRMX(3,J,1)+T(1)*(TC(1)*PI(3,J)+TB(1)*TAU(3,J)) 349
ELTRMX(4,J,1) = ELTRMX(4,J,1)+T(1)*(TC(2)*PI(3,J)+TB(2)*TAU(3,J)) 350
IF ( K .EQ. N ) GO TO 75 351
ELTRMX(1,J,2) = ELTRMX(1,J,2)+T(1)*(TB(1)*PI(3,J)-TC(1)*TAU(3,J)) 352
ELTRMX(2,J,2) = ELTRMX(2,J,2)+T(1)*(TB(2)*PI(3,J)-TC(2)*TAU(3,J)) 353
ELTRMX(3,J,2) = ELTRMX(3,J,2)+T(1)*(TC(1)*PI(3,J)-TB(1)*TAU(3,J)) 354
ELTRMX(4,J,2) = ELTRMX(4,J,2)+T(1)*(TC(2)*PI(3,J)-TB(2)*TAU(3,J)) 355
GO TO 80 356
75 ELTRMX(1,J,2) =ELTRMX(1,J,2)+T(1)*(-TB(1)*PI(3,J)+TC(1)*TAU(3,J)) 357
ELTRMX(2,J,2) =ELTRMX(2,J,2)+T(1)*(-TB(2)*PI(3,J)+TC(2)*TAU(3,J)) 358
ELTRMX(3,J,2) =ELTRMX(3,J,2)+T(1)*(-TC(1)*PI(3,J)+TB(1)*TAU(3,J)) 359
ELTRMX(4,J,2) =ELTRMX(4,J,2)+T(1)*(-TC(2)*PI(3,J)+TB(2)*TAU(3,J)) 360
80 CONTINUE 361
C 362
IF ( T(4) .LT. 1.0E-14 ) GO TO 100 363
N = N + 1 364
DO 90 J = 1,JX 365
PI(1,J) = PI(2,J) 366
PI(2,J) = PI(3,J) 367
TAU(1,J) = TAU(2,J) 368
TAU(2,J) = TAU(3,J) 369
90 CONTINUE 370
FNAP = FNA 371
FNBP = FNB 372
IF ( N .LE. NMX2 ) GO TO 65 373
WRITE( 6,8 ) 374
STOP 36 375
100 DO 120 J = 1,JX 376
DO 120 K = 1,2 377
DO 115 I= 1,4 378
T(I) = ELTRMX(I,J,K) 379
115 CONTINUE 380
ELTRMX(2,J,K) = T(1)**2 + T(2)**2 381
ELTRMX(1,J,K) = T(3)**2 + T(4)**2 382
ELTRMX(3,J,K) = T(1) * T(3) + T(2) * T(4) 383
ELTRMX(4,J,K) = T(2) * T(3) - T(4) * T(1) 384
120 CONTINUE 385
T(1) = 2.0 * RX**2 386
QEXT = QEXT * T(1) 387
QSCAT = QSCAT * T(1) 388
CTBRQS = 2.0 * CTBRQS * T(1) 389
C 390
RETURN 391
C 392
5 FORMAT( 10X,' THE VALUE OF THE SCATTERING ANGLE IS GREATER THAN 393
1 90.0 DEGREES. IT IS ', E15.4 ) 394
6 FORMAT( // 10X, 'PLEASE READ COMMENTS.' // ) 395
7 FORMAT( // 10X, 'THE VALUE OF THE ARGUMENT JX IS GREATER THAN IT') 396
8 FORMAT( // 10X, 'THE UPPER LIMIT FOR ACAP IS NOT ENOUGH. SUGGEST 397
1 GET DETAILED OUTPUT AND MODIFY SUBROUTINE' // ) 398
C 399
END 400
SUBROUTINE REGIS (FNU,ND1,NNU,T,ND2,NT,XN2N2,XCH4CH4,XN2CH4,XN2H2)
C THIS SUBROUTINE RETURNS THE PRESSURE INDUCED ABSORBTION COEFFICIENTS
C FOR N2-N2, CH4-CH4, N2-CH4 + CH4-N2, H2-N2 + N2-H2.
C FNU IS THE WAVENUMBER ARRAY
C T IS THE TEMPERATURE ARRAY
C NNU IS THE NUMBER OF WAVENUMBERS
C NT IS THE NUMBER OF TEMPERATURES
C ND1 IS THE DIMENSION OF THE WAVENUMERS
C ND2 IS THE DIMENSION OF THE TEMPERATURES
C XN2N2(ND1,ND2) IS THE N2-N2 COEFFICIENT
C XCH4CH4(ND1,ND2) IS THE CH4-CH4 COEFFICIENT
C XN2CH4(ND1,ND2) IS THE N2-CH4 + CH4-N2 COEFFICIENT
C XN2H2(ND1,ND2) IS THE N2-H2 + H2-N2 COEFFICIENT
C FROM REGIS COURTIN ADAPTED BY CPM.
IMPLICIT REAL*8 (A-H,O-Z)
DIMENSION FNU(NNU),T(NT),
& XN2N2(ND1,ND2),XCH4CH4(ND1,ND2),XN2H2(ND1,ND2),XN2CH4(ND1,ND2)
DIMENSION F8MN(100),F8NH(100),F8HN(100),F10(100),
2 XN2R(100),XN2T(100),XCH4R(100),XCH4T(100),XH2R(100),XH2T(100)
DATA EPS1/71.4D0/,EPS2/148.6D0/,EPS3/36.8D0/
DO 1 IT=1,NT
EPS4=DSQRT(EPS1*EPS2)
EPS5=DSQRT(EPS1*EPS3)
Y=4.D0*EPS4/T(IT)
F8MN(IT)=FI8(Y)
F10(IT)=FI10(Y)
Y=4.D0*EPS5/T(IT)
Z=FI8(Y)
F8NH(IT)=Z
F8HN(IT)=Z
Y=4.D0*EPS1/T(IT)
Z=FI8(Y)
F8MN(IT)=F8MN(IT)/Z
F8HN(IT)=F8HN(IT)/Z
Y=4.D0*EPS2/T(IT)
F10(IT)=F10(IT)/FI10(Y)
Y=4.D0*EPS3/T(IT)
F8NH(IT)=F8NH(IT)/FI8(Y)
1 CONTINUE
CALL PIAN2(0.D0,NT,XN2R,XN2T,T)
CALL PIACH4(0.D0,NT,XCH4R,XCH4T,T)
CALL PIAH2(0.D0,NT,XH2R,XH2T,T)
DO 2 INU=1,NNU
CALL OPAN2(FNU(INU),NT,XN2R,XN2T)
CALL OPACH4(FNU(INU),NT,XCH4R,XCH4T)
CALL OPAH2(FNU(INU),NT,XH2R,XH2T)
DO 11 IT=1,NT
XN2N2(INU,IT)=XN2R(IT)+XN2T(IT)
XCH4CH4(INU,IT)=XCH4R(IT)+XCH4T(IT)
C
C LINES 1 AND 3 ARE CORRECTION FACTORS INTRODUCED IN ORDER TO FIT
C THE MEASUREMENTS OF DAGG ET AL (1986) BETWEEN 126 AND 212 K.
C
XN2CH4(INU,IT)=F8MN(IT)*(2.067D0*XN2R(IT)+2.865D0*XN2T(IT))
1 * 2.48D0/(T(IT)**0.184) +
2 F10(IT)*(0.480D0*XCH4R(IT)+0.374D0*XCH4T(IT))
3 * 1.54D0
C
C LINE 3 IS A CORRECTION FACTOR INTRODUCED IN ORDER TO FIT
C THE MEASUREMENTS OF DORE ET AL (1986) BETWEEN 91 AND 298 K.
C OTHER CORRECTION FACTORS ARE IMBEDDED IN THE ROUTINE PIAH2.
C
XN2H2(INU,IT)= F8NH(IT)*(2.543D0*XH2R(IT)+1.445D1*XH2T(IT)) +
2 F8HN(IT)*(0.360D0*XN2R(IT)+0.246D0*XN2T(IT))
3 * 0.30D0
11 CONTINUE
2 CONTINUE
RETURN
END
SUBROUTINE PIAN2(NU,NT,XROT,XTRA,T)
C
C THIS ROUTINE WAS FIRST WRITTEN BY REGIS COURTIN IN 1980 :
C THE PARAMETERS WERE DETERMINED BY THE MEASUREMENTS OF
C BUONTEMPO ET AL (1975), JOURNAL OF CHEMICAL PHYSICS 63, 2570.
C LAST MODIFIED BY REGIS COURTIN (APRIL 23,1986) :
C THE QUADS1 PARAMETER WAS ADJUSTED TO FIT THE MEASUREMENTS OF
C DAGG ET AL (1985), CANADIAN JOURNAL OF PHYSICS 63, 625.
C
IMPLICIT REAL*8 (A-H,O-Z) 192.
REAL*8 NU,I8,NSQR,K 193.
DIMENSION T(100) 194.
DIMENSION XROT(100),XTRA(100) 195.
DIMENSION XJ(32),G(32),E(32),WW(32),CG(32),RHO(100,32), 196.
$I8(100),QUADS1(100),TAU1(100),TAU2(100),A2(100) 197.
DATA JMAX1/30/,JMAX2/32/,PI/3.141592654D0/,H/1.05459D-27/, 198.
$C/2.997925D10/,K/1.38054D-16/,HCK/1.43892D0/,XNLOS/2.687D19/, 199.
$ASQR/3.3124D0/,EPSK/71.4D0/,PW1/.639D0/,PW2/.423D0/,
$DSDT/-.0204D0/
DO 8 IT=1,NT 201.
X=4.D0*EPSK/T(IT) 202.
I8(IT)=FI8(X) 203.
8 CONTINUE 204.
DO 10 J=1,JMAX2 205.
XJ(J)=DFLOAT(J-1) 206.
G(J)=6.D0 207.
IF (MOD(J,2).EQ.0) G(J)=3.D0 208.
E(J)=AZENER(0.D0,XJ(J)) 209.
10 CONTINUE 210.
DO 20 J=1,JMAX1 211.
CG(J)=1.5D0*(XJ(J)+1.D0)*(XJ(J)+2.D0)/(2.D0*XJ(J)+3.D0) 212.
WW(J)=2.D0*PI*C*(E(J+2)-E(J)) 213.
20 CONTINUE 214.
NSQR=(XNLOS*1.0D-24)**2 215.
C1=4.D0*PI*PI*NSQR/H/C 216.
DO 90 IT=1,NT 217.
QUADS1(IT)=27.72D0+DSDT*T(IT) 218.
TAU1(IT)=9.797D-12/(T(IT)**PW1) 219.
TAU2(IT)=1.518D-12/(T(IT)**PW2) 220.
SUM=0.D0 221.
DO 50 J=1,JMAX2 222.
ARG=HCK*E(J)/T(IT) 223.
RHO(IT,J)=DEXP(-ARG) 224.
SUM=SUM+(2.D0*XJ(J)+1.D0)*G(J)*RHO(IT,J) 225.
50 CONTINUE 226.
DO 55 J=1,JMAX2 227.
RHO(IT,J)=G(J)*RHO(IT,J)/SUM 228.
55 CONTINUE 229.
A2(IT)=0.D0 230.
DO 60 J=1,JMAX1 231.
A2(IT)=A2(IT)+XJ(J)*(XJ(J)+1.D0)*(2.D0*XJ(J)+1.D0)*RHO(IT,J)/ 232.
$ ((2.D0*XJ(J)-1.D0)*(2.D0*XJ(J)+3.D0)) 233.
60 CONTINUE 234.
90 CONTINUE 235.
RETURN 236.
ENTRY OPAN2(NU,NT,XROT,XTRA) 237.
W=2.D0*PI*C*NU 238.
DO 200 IT=1,NT 239.
CZ=HCK*NU/T(IT) 240.
CW=W*(1.D0-DEXP(-CZ)) 241.
XTRA(IT)=A2(IT)*GAMMB(W,T(IT),TAU1(IT),TAU2(IT)) 242.
XROT(IT)=0.D0 243.
DO 125 J=1,JMAX1 244.
XROT(IT)=XROT(IT)+CG(J)* 245.
$(RHO(IT,J )*GAMMB(W-WW(J),T(IT),TAU1(IT),TAU2(IT)) 246.
$+RHO(IT,J+2)*GAMMB(W+WW(J),T(IT),TAU1(IT),TAU2(IT))) 247.
125 CONTINUE 248.
COEF=CW*C1*K*ASQR*QUADS1(IT)*I8(IT) 249.
XROT(IT)=XROT(IT)*COEF 250.
XTRA(IT)=XTRA(IT)*COEF 251.
200 CONTINUE 252.
RETURN 253.
END 254.
DOUBLE PRECISION FUNCTION AZENER(V,J) 255.
IMPLICIT REAL*8 (A-H,O-Z) 256.
REAL*8 J,JP 257.
VP=V+0.5D0 258.
JP=J+1.D0 259.
E=2359.61D0*VP-14.456D0*VP**2+7.51D-3*VP**3+ 260.
$(1.9987D0-1.87D-2*VP)*J*JP-(5.8D-6-1.D-9*VP)*J**2*JP**2 261.
AZENER=E 262.
RETURN 263.
END 264.
SUBROUTINE PIACH4(NU,NT,XROT,XTRA,T) 265.
IMPLICIT REAL*8 (A-H,O-Z) 266.
REAL*8 NU,I10,NSQR,K 267.
DIMENSION T(100) 268.
DIMENSION FOX(60),XJ(23),G(23),E(23),WW(3,23),RHO(100,23), 269.
$I10(100),XROT(100),XTRA(100) 270.
DATA FOX /1.00,1.00,4.12,1.00,0.52,0.81,1.00,0.70,0.75,0.80, 271.
$0.69,1.24,0.44,1.29,1.05,0.56,0.79,0.79,0.89,1.16,1.27,0.86, 272.
$0.95,1.00,0.72,1.16,0.92,1.00,0.92,1.18,0.90,1.18,1.09,0.85, 273.
$0.96,0.93,0.82,1.18,1.17,0.95,0.96,1.00,0.85,1.06,0.90,15*1.00/ 274.
DATA JMAX1/20/,JMAX2/23/,PI/3.141592654D0/,H/1.05459D-27/, 275.
$C/2.997925D10/,K/1.38054D-16/,HCK/1.43892D0/,XNLOS/2.687D19/, 276.
$ASQR/6.7081D0/,EPSK/148.6D0/,B/5.25D0/,ANORM/4.797D17/ 277.
DO 8 IT=1,NT 278.
X=4.D0*EPSK/T(IT) 279.
I10(IT)=FI10(X) 280.
8 CONTINUE 281.
DO 10 J=1,JMAX2 282.
XJ(J)=DFLOAT(J-1) 283.
E(J)=B*XJ(J)*(XJ(J)+1.D0) 284.
10 CONTINUE 285.
G(1)=5.D0 286.
G(2)=3.D0 287.
G(3)=5.D0 288.
G(4)=11.D0 289.
G(5)=13.D0 290.
G(6)=11.D0 291.
DO 15 J=7,JMAX2 292.
15 G(J)=G(J-6)+16.D0 293.
DO 20 J=1,JMAX1 294.
DO 20 I=1,3 295.
WW(I,J)=2.D0*PI*C*(E(J+I)-E(J)) 296.
20 CONTINUE 297.
NSQR=(XNLOS*1.0D-24)**2 298.
C1=64.D0/63.D0*PI*PI*NSQR/H/C/ANORM 299.
QUADS1=11.10D0 298.
TAU1=9.00D-14 301.
TAU2=13.60D-14 302.
DO 90 IT=1,NT 303.
SUM=0.D0 304.
DO 50 J=1,JMAX2 305.
ARG=HCK*E(J)/T(IT) 306.
RHO(IT,J)=DEXP(-ARG) 307.
SUM=SUM+(2.D0*XJ(J)+1.D0)*G(J)*RHO(IT,J) 308.
50 CONTINUE 309.
DO 55 J=1,JMAX2 310.
RHO(IT,J)=RHO(IT,J)/SUM 311.
55 CONTINUE 312.
90 CONTINUE 313.
RETURN 314.
ENTRY OPACH4(NU,NT,XROT,XTRA) 315.
W=2.D0*PI*C*NU 316.
DO 200 IT=1,NT 317.
XROT(IT)=0.D0 318.
XTRA(IT)=0.D0 319.
CZ=HCK*NU/T(IT) 320.
CW=W*(1.D0-DEXP(-CZ)) 321.
DO 125 J=1,JMAX1 322.
XTRA(IT)=XTRA(IT)+(2.D0*XJ(J)+1.D0)*(2.D0*XJ(J)+1.D0)* 323.
$ RHO(IT,J)*GAMMB(W,T(IT),TAU1,TAU2) 324.
DO 25 I=1,3 325.
K=3*(J-1)+I 326.
XROT(IT)=XROT(IT)+FOX(K)*(2.D0*XJ(J)+1.D0)*(2.D0*XJ(J+I)+1.D0)* 327.
$(RHO(IT,J )*GAMMB(W-WW(I,J),T(IT),TAU1,TAU2) 328.
$+RHO(IT,J+I)*GAMMB(W+WW(I,J),T(IT),TAU1,TAU2)) 329.
25 CONTINUE 330.
125 CONTINUE 331.
COEF=CW*C1*K*ASQR*QUADS1*I10(IT) 332.
XROT(IT)=XROT(IT)*COEF 333.
XTRA(IT)=XTRA(IT)*COEF 334.
200 CONTINUE 335.
RETURN 336.
END 337.
DOUBLE PRECISION FUNCTION FI8(X) 355.
IMPLICIT REAL*8 (A-H,O-Z) 356.
FI8=.3737D0/X+3.718903D0-.1908D0*X+.1617D0*X*X-3.633D-3*X**3 357.
RETURN 358.
END 359.
DOUBLE PRECISION FUNCTION FI10(X) 360.
IMPLICIT REAL*8 (A-H,O-Z) 361.
FI10=.353D0/X+3.154000D0-.3060D0*X+.1410D0*X*X-2.887D-3*X**3 362.
RETURN 363.
END 364.
DOUBLE PRECISION FUNCTION GAMMB(W,T,TAU1,TAU2) 365.
IMPLICIT REAL*8 (A-H,O-Z) 366.
DATA HK/7.638967D-12/,PI/3.141592654D0/
cfix ,TAU10/0.D0/,TAU20/0.D0/ 367.
C NOTE. HK = 1.05459D-27/1.38054D-16 368.
cfix IF (TAU1.NE.TAU10) GO TO 10 369.
cfix IF (TAU2.EQ.TAU20) GO TO 20 370.
cfix 10 CONTINUE 371.
TAU12=TAU1*TAU1 372.
TAU22=TAU2*TAU2 373.
HBH = 0.5D0*HK/T 374.
Z2=DSQRT ( TAU22 + HBH**2) / TAU1 375.
cfix TAU10 = TAU1 376.
cfix TAU20=TAU2 377.
cfix 20 CONTINUE 378.
WSQR=W*W 379.
Z = DSQRT (1.D0+ WSQR*TAU12) * Z2 380.
C COMPUTE THE MODIFIED BESSEL FUNCTION OF THE SECOND KIND (K1) 381.
C USING AN UNPUBLISHED APPROXIMATION GIVEN BY COHEN. 382.
F=1.5707963D0*(Z+0.5616D0)/(Z+0.4619D0) 383.
BK1=DSQRT(1.D0+Z*F) 384.
BK1=BK1*DEXP(-Z) 385.
GAMMA=TAU1/PI*DEXP(TAU2/TAU1+HBH*W)*BK1 / (1.D0+WSQR*TAU12) 386.
GAMMB=GAMMA 387.
RETURN 388.
END 389.
SUBROUTINE PIAH2(NU,NT,XROT,XTRA,T) 00001000
C 00002000
C XKH2 IS CALLED BY TH2HE TO COMPUT XK1, A FACTOR IN THE 00003000
C H2-H2 COLLISION-INDUCED ABSORPTION. 00004000
C 00005000
C THIS VERSION INCLUDES RESULTS OF RECENT LOW-TEMPERATURE FITS 00006000
C BY G. BIRNBAUM ET AL (NBS), AS IMPLEMENTED BY VIRGIL KUNDE 00007000
C AND GORDON BJORAKER, AUGUST 1982. 00008000
C 00009000
C THIS VERSION ALLOWS FOR A NON-EQUILIBRIUM PARA-HYDROGEN FRACTION. 00010000
C 00011000
C LAST CHANGED BY JOHN HORNSTEIN CSC FEB 28,1983 00012000
C THIS CHANGE DELETED THE DOUBLE-TRANSITIONS, WHICH HAD BEEN 00013000
C INCORRECTLY FORMULATED IN BACHET ET AL. THE DELETION IS BY 00014000
C COMMENTING OUT WITH A "CD". 00015000
C 00016000
C THREE PARAMETERS ARE USED TO FIT THE H2-H2 OPACITY FROM 00017000
C 0 TO 2000 CM-1: STREN, TAU1, & TAU2. THE TEMPERATURE 00018000
C DEPENDENCE OF EACH PARAMETER IS HANDLED DIFFERENTLY. 00019000
C STREN WAS FIT TO EXPERIMENTAL VALUES IN DORE ET AL AND 00020000
C BACHET ET AL (1982) USING A LINEAR TEMPERATURE RELATION. 00021000
C TAU1 & TAU2 FOLLOW A POWER LAW OF THE FORM 00022000
C TAU = TAU0 * (T/T0)**BT AS SUGGESTED BY DORE ET AL. 00023000
C CONSTANTS FOR TAU1 ARE FROM DORE ET AL 00024000
C TAU2 = 1.23 * (TAU2 EVALUATED USING DORE POWER LAW). THIS MATCHES 00025000
C BACHET ET AL VALUES AT 195K & 297K (TAIL WT SINGLE + 00026000
C DOUBLE TRANSITION COEFFICIENTS). 00027000
C 00028000
C PARAMETER DESCRIPTION 00029000
C 00030000
C STRENGTH PARAMETER: STREN (SEE BACHET ET AL EQUATION 4B) 00031000
C STREN IS REALLY S/KB, TO AVOID LARGE POWERS OF 10. 00032000
C 00033000
C STREN = 6 * * OMEGA**2 * I8(R/SIGMA) / K 00034000
C (UNITS: KELVIN*ANGSTROM**6) 00035000
C 00036000
C IS MEAN SQUARE POLARIZABILITY UNIT: ANGSTROM**6 00037000
C OMEGA IS QUADRUPOLE MOMENT UNIT: ESU*CM**2 00038000
C R IS DISTANCE, SIGMA IS MOLECULAR SEPARATION, X = R/SIGMA 00039000
C I8(X) = 4*PI*INT(X**-8*EXP(-4.*E/(K*T)*(X**-12-X**-6))*X**2*DX) 00040000
C E & SIGMA ARE PARAMETERS OF LENNARD-JONES INTERMOLECULAR POTENTIAL 00041000
C INT(...DX) IS INTEGRAL 0 TO INFINITY 00042000
C 00043000
C TIME PARAMETERS: TAU1 & TAU2 00044000
C 00045000
C TAU1 CONTROLS HALF WIDTH NEAR LINE CENTER 00046000
C TAU2 CONTROLS EXPONENTIAL DECAY OF WINGS 00047000
C 00048000
C 00049000
C DOUBLE TRANSITIONS ARE INCLUDED. 00050000
C THE ABSORPTION DUE TO DOUBLE TRANSITIONS IS PROPORTIONAL TO C2 00051000
C C2 = 4./45. * KAPPA**2 00052000
C KAPPA = (APL-APP)/(1./3. * (APL+2.*APP)) 00053000
C APL & APP ARE POLARIZABILITIES PARALLEL 00054000
C AND PERPENDICULAR TO INTERNUCLEAR AXIS 00055000
C KAPPA = .375 FOR GROUND VIB STATE & J=0,1 OR 2 SEE KOLOS & WOLNIEWICZ00056000
C JOURNAL OF CHEMICAL PHYSICS 46, 1426 (1967) TABLE 3
IMPLICIT REAL*8 (A-H,O-Z)
REAL*8 NU,NSQR,K 00058000
LOGICAL ILIST
DIMENSION XROT(251),XTRA(251),XJ(10),G(10),A2(250),RHO(250,10), 0006000
A E(10),WN(10),WW(10),CG(10),BH(250),STREN(250),TAU1R(250),
B TAU2R(250),TAU1T(250),TAU2T(250),COREG0(2),COREG(10,2) 0006100
DIMENSION T(251) 0006200
DATA JMAX1/8/,JMAX2/10/, 00063000
. PI/3.141593/,HBAR/1.05450D-27/,C/2.997925D+10/,K/1.38054D-16/, 00064000
. HCK/1.43879/,XNLOS/2.687E19/,S0/178./,DSDT/.4091/, 00065000
. TAU1R0/7.25/,TAU2R0/3.45/,TAU1T0/7.25/,TAU2T0/3.45/,
. BT1/-.605/,BT2/-.607/,T0/273.15/,C2/0.0125/
C 00067000
DIMENSION FPARA(99) 00068000
C 00070000
C***********************************************************************00071000
C 00072000
C 00073000
C G(J) ARE STATISTICAL WTS: 00075000
C G(J)=1 EVEN ROTATIONAL STATES, 00076000
C G(J)=3 ODD ROTATIONAL STATES. 00077000
ILIST=.FALSE.
C SET THE ORTHO/PARA TO EQUILIBRIUM -CPM
NPARA=-1.
DO 10 J = 1,JMAX2 00078000
COREG(J,1)=1.D0
COREG(J,2)=0.D0
XJ(J) = FLOAT(J-1) 00079000
G(J) = 1.0 00080000
IF (MOD(J,2) .EQ. 0) G(J) = 3.0 00081000
10 E(J) = H2ENER(0.0,XJ(J)) 00082000
C ****************** WARNING ***********************
C CHANGE MADE BY REGIS COURTIN (APRIL 1986).
C THE FOLLOWING COEFFICIENTS ARE INTRODUCED FOR THE CASE
C OF COLLISIONS BETWEEN H2 AND N2 MOLECULES.
C THE CORRECTED ABSORPTION COEFFICIENTS FIT THE DATA OF
C DORE ET AL (1986), PREPRINT.
C ******************************************************
COREG0(1)=22.35D0
COREG0(2)=-0.56D0
COREG(1,1)=0.36D0
COREG(1,2)=0.18D0
COREG(2,1)=10.91D0
COREG(2,2)=-0.433D0
C ******************************************************
C CG(J) = (2*J+1) * (CLEBSCH-GORDAN COEFF )**2 00083000
DO 20 J = 1,JMAX1 00084000
CG(J) = (2.0*XJ(J)+1.0) * (3.0*(XJ(J)+1.0)*(XJ(J)+2.0)) / 00085000
. (2.0*(2.0*XJ(J)+1.0)*(2.0*XJ(J)+3.0)) 00086000
WN(J) = E(J+2) - E(J) 00087000
20 WW(J) = 2.0 * PI * C * WN(J) 00088000
NSQR = (XNLOS*1.0D-24)**2 00089000
C1 = 2. * PI**2 * NSQR/(3. * HBAR * C) 00090000
C LIST HEADING FOR ORTHO AND PARA PROFILES: 00091000
IF ( (NPARA .LE. 0) .OR. (.NOT. ILIST)) GO TO 7000 00092000
WRITE(6,5000) 00093000
5000 FORMAT(//,' EQUILIBRIUM AND ACTUAL FRACTIONS OF PARA-', 00094000
. ' AND ORTHO-HYDROGEN:',/,' LAYER',7X,'FPARA(EQU)', 00095000
. 7X,'FPARA(HERE)',7X,'FORTHO(EQU)',7X,'FORTHO(HERE)') 00096000
7000 DO 90 IT = 1,NT 00097000
BH(IT) = HBAR / (K*T(IT)) 00098000
C EVALUATE STREN, TAU1R, TAU2R AT DESIRED T 00099000
STREN(IT) = (S0 + DSDT*T(IT)) 00100000
TAU1R(IT) = TAU1R0 * (T(IT)/T0)**BT1 * 1.0D-14 00101000
TAU2R(IT) = TAU2R0 * (T(IT)/T0)**BT2 * 1.0D-14 00102000
TAU1T(IT) = TAU1T0 * (T(IT)/T0)**BT1 * 1.0D-14
TAU2T(IT) = TAU2T0 * (T(IT)/T0)**BT2 * 1.0D-14
C COMPUTE THE FULL PARTION FUNCTION Z, USED IN EQUILIBRIUM, 00103000
C WHERE ORTHO AND PARA ARE CONVENIENTLY TREATED AS DIFFERENT 00104000
C STATES OF THE SAME SPECIES. 00105000
Z = 0.0 00106000
DO 50 J = 1,JMAX2 00107000
RHO(IT,J) = EXP(-1.*HCK*E(J)/T(IT)) 00108000
50 Z = Z + (2.0*XJ(J)+1.0)*G(J)*RHO(IT,J) 00109000
IF (NPARA .LE. 0) GO TO 54 00110000
C COMPUTE THE PARTITION FUNCTIONS ZPARA AND ZORTHO USED FOR 00111000
C NON-EQUILIBRIUM RATIOS, WHERE IT IS CONVENIENT TO TREAT 00112000
C ORTHO AND PARA AS DISTINCT SPECIES. THE NUCLEAR SPIN 00113000
C FACTORS G(J) CANCEL OUT IN THIS CASE. 00114000
ZPARA = 0. 00115000
ZORTHO = 0. 00116000
DO 1000 J=1,JMAX2,2 00117000
ZPARA = ZPARA + (2.*XJ(J) + 1.)*RHO(IT,J) 00118000
JJ = J + 1 00119000
1000 ZORTHO = ZORTHO + (2.*XJ(JJ) + 1.)*RHO(IT,JJ) 00120000
C COMPUTE AND LIST THE EQUILIBRIUM AND ACTUAL PROFILES: 00121000
FEPARA = ZPARA/Z 00122000
FEORTH = 3.*ZORTHO/Z 00123000
FORTHO = 1. - FPARA(IT) 00124000
IF (ILIST) WRITE(6,2000) IT, FEPARA, FPARA(IT), FEORTH, FORTHO00125000
2000 FORMAT(' ',I4,2F15.5,10X,2F15.5) 00126000
C FORM A NEW RHO WHICH EQUALS RHO(PARA) WHEN XJ IS EVEN 00127000
C (INDEX J IS ODD) AND EQUALS RHO(ORTHO) WHEN XJ IS ODD. 00128000
DO 3000 J=1,JMAX2,2 00129000
RHO(IT,J) = FPARA(IT)*RHO(IT,J)/ZPARA 00130000
JJ = J + 1 00131000
3000 RHO(IT,JJ) = FORTHO*RHO(IT,JJ)/ZORTHO 00132000
GO TO 57 00133000
C EQUILIBRIUM HYDROGEN STATISTICAL WEIGHTS 00134000
54 DO 55 J = 1,JMAX2 00135000
55 RHO(IT,J) = G(J)*RHO(IT,J) / Z 00136000
57 A2(IT) = 0.0 00137000
DO 60 J = 1,JMAX1 00138000
60 A2(IT) = A2(IT) + XJ(J)*(XJ(J)+1.0)*(2.0*XJ(J)+1.0)*RHO(IT,J)/00139000
. ((2.0*XJ(J)-1.0)*(2.0*XJ(J)+3.0)) 00140000
90 CONTINUE 00141000
RETURN
ENTRY OPAH2(NU,NT,XROT,XTRA)
W = 2.0 * PI * C * NU 00142000
DO 200 IT = 1,NT 00143000
DBL = 0. 00144000
FR = 0. 00145000
C EVALUATE TRANSLATIONAL PART OF SHAPE FACTOR F 00146000
FT = A2(IT) * GAMFCN(W,T(IT),TAU1T(IT),TAU2T(IT)) 00147000
.* COREG0(1)*T(IT)**COREG0(2)
DO 125 J = 1,JMAX1 00148000
C EVALUATE THE ROTATIONAL PART OF F 00149000
FR = FR + CG(J) * (RHO(IT,J)*GAMFCN(W-WW(J),T(IT),TAU1R(IT), 00150000
. TAU2R(IT))+RHO(IT,J+2)*GAMFCN(W+WW(J),T(IT),TAU1R(IT), 00151000
. TAU2R(IT))) 00152000
. * COREG(J,1)*T(IT)**COREG(J,2)
C DBL = DBL + CG(J) * (RHO(IT,J) + RHO(IT,J+2)) 00153000
125 CONTINUE 00154000
C ADD ON PART OF DOUBLE TRANSITION OPACITY (BACHET ET AL EQN 11) 00155000
CD FR = FR * (1.0 + C2 * A2(IT)) 00156000
C MORE DOUBLE TRANS (BACHET ET AL EQN 12) 00157000
CD FT = FT * (1.0 + C2 * A2(IT)) 00158000
C AND STILL MORE DOUBLE TRANS (BACHET ET AL EQN 13) 00159000
CD FT = FT + C2 * DBL*DBL*GAMFCN(W,T(IT),TAU1T(IT),TAU2T(IT)) 00160000
F = FT + FR 00161000
C EVALUATE & ADD ON DOUBLE TRANSITIONS (BACHET ET AL EQN 14) 00162000
CD DO 130 J1=1,4 00163000
CD DO 130 J2=1,4 00164000
CD IF (J1.NE.(J1+2).OR.J2.NE.(J1+2)) 00165000
CD A F = F + C2*RHO(IT,J1)*CG(J1)*RHO(IT,J2)*CG(J2) 00166000
CD B * GAMFCN(W-WW(J1)-WW(J2),T(IT),TAU1R(IT),TAU2R(IT)) 00167000
CD IF ((J1+2).NE.(J2+2).OR.J2.NE.J1) 00168000
CD A F = F + C2*RHO(IT,J1+2)*CG(J1)*RHO(IT,J2)*CG(J2) 00169000
CD B * GAMFCN(W+WW(J1)-WW(J2),T(IT),TAU1R(IT),TAU2R(IT)) 00170000
CD IF ((J1+2).NE.J2.OR.(J2+2).NE.J1) 00171000
CD A F = F + C2*RHO(IT,J1+2)*CG(J1)*RHO(IT,J2+2)*CG(J2) 00172000
CD B * GAMFCN(W+WW(J1)+WW(J2),T(IT),TAU1R(IT),TAU2R(IT)) 00173000
CD IF (J1.NE.J2.OR.(J2+2).NE.(J1+2)) 00174000
CD A F = F + C2*RHO(IT,J1)*CG(J1)*RHO(IT,J2+2)*CG(J2) 00175000
CD B * GAMFCN(W-WW(J1)+WW(J2),T(IT),TAU1R(IT),TAU2R(IT) 00176000
CD130 CONTINUE 00177000
XROT(IT) = C1 * K * W * (1.0-EXP(-BH(IT)*W)) * STREN(IT)*FR 00178000
XTRA(IT) = XROT(IT)*FT/FR
200 CONTINUE 00179000
RETURN 00180000
END 00183000
DOUBLE PRECISION FUNCTION GAMFCN(W,T,TAU1,TAU2)
IMPLICIT REAL*8 (A-H,O-Z)
C
C COMPUTES THE LINE SHAPE FOR PRESSURE-INDUCED H2-H2 AND H2-HE
C TRANSITIONS, FROM THE SEMI-EMIRICAL FORMULAE OF BIRNBAUM ET AL.
C SEE EG., GEORGE BIRNBAUM AND E. RICHARD COHEN, CANADIAN JOURNAL
C OF PHYSICS, VOL. 54, 593 (1976).
C NOTE THAT "GAMFCN" IS A POOR NAME FOR THIS ROUTINE; THE GAMMA
C OF BIRNBAUM AND COHEN IS NOT THE USUAL "GAMMA FUNCTION".
C
DATA HK/7.638315D-12/,PI/3.141593/,TAU10/0.0/,TAU20/0.0/
C NOTE: HK = 1.05450D-27 / 1.38054D-16
C
C***********************************************************************
IF (TAU1 .NE. TAU10) GO TO 10
IF (TAU2 .EQ. TAU20) GO TO 20
10 TAU12 = TAU1 * TAU1
TAU22 = TAU2 * TAU2
HBH = 0.5 * HK / T
Z2 = DSQRT(TAU22 + HBH**2) / TAU1
TAU10 = TAU1
TAU20 = TAU2
20 WSQR = W * W
Z = DSQRT(1.0+WSQR*TAU12) * Z2
IF (Z .LE. 1.0) GO TO 50
C COMPUTE K1 BESSEL FUNCTION USING POLYNOMIAL APPROXIMATION
A = 1.0 / Z
BK1 = 1.253314 + .4699927*A
B = A * A
BK1 = BK1 - .1468583*B
B = B * A
BK1 = BK1 + .1280427*B
B = B * A
BK1 = BK1 - .1736432*B
B = B * A
BK1 = BK1 + .2847618*B
B = B * A
BK1 = BK1 - .4594342*B
B = B * A
BK1 = BK1 + .6283381*B
B = B * A
BK1 = BK1 - .6632295*B
B = B * A
BK1 = BK1 + .5050239*B
B = B * A
BK1 = BK1 - .2581304*B
B = B * A
BK1 = BK1 + .7880001D-01*B
B = B * A
BK1 = BK1 - .1082418D-01*B
BK1 = EXP(-Z) * BK1 * DSQRT(A)
GO TO 100
C COMPUTE K1 BESSEL FUNCTION USING SERIES EXPANSION
50 A = 0.5 * Z
B = .5772157 + DLOG(A)
C = A * A
BK1 = 1.0/Z + A*(B-0.5)
A = A * C
BK1 = BK1 + A*.2500000D+00*(0.5+(B-1.500000)*2.0)
A = A * C
BK1 = BK1 + A*.2777777D-01*(0.5+(B-1.833333)*3.0)
A = A * C
BK1 = BK1 + A*.1736110D-02*(0.5+(B-2.083333)*4.0)
A = A * C
BK1 = BK1 + A*.6944439D-04*(0.5+(B-2.283333)*5.0)
A = A * C
BK1 = BK1 + A*.1929009D-05*(0.5+(B-2.449999)*6.0)
A = A * C
BK1 = BK1 + A*.3936752D-07*(0.5+(B-2.592855)*7.0)
A = A * C
BK1 = BK1 + A*.6151173D-09*(0.5+(B-2.717855)*8.0)
100 CONTINUE
GAMFCN = TAU1/PI * EXP(TAU2/TAU1+HBH*W) * Z*BK1 / (1.0+WSQR*TAU12)
RETURN
END
DOUBLE PRECISION FUNCTION H2ENER(V,J)
C H2ENER SUBROUTINE OF THE INV PROGRAM. COMPUTES THE ENERGY
C (IN CM-1) OF A VIBRATION-ROTATION STATE OF A HYDROGEN MOLECULE.
C THE VIBRATION QUANTUM NUMBER IS V, AND THE QUANTUM NUMBER FOR
C THE RIGID BODY ANGULAR MOMENTUM IS J (A REAL*4 QUANTITY).
C THE FIRST LINE OF THE FORMULA FOR E IS THE CONTRIBUTION FROM
C PURE VIBRATION, INCLUDING ANHARMONIC TERMS. THE OTHER LINES
C ACCOUNT FOR COUPLED VIBRATION AND ROTATION, INCLUDING
C CENTRIFUGAL DISTORTION. (THE (J(J+1))**2 AND (J(J+1))**3
C PROVIDE FOR CENTRIFUGAL DISTORTION; THE VP AND VP**2 IN THE
C ROTATION TERMS PROVIDE FOR COUPLING BETWEEN VIBRATION AND
C ROTATION.)
C THE FORMULA OF COHEN AND BIRNBAUM(1981),
C NU = 59.3392*(J(J+1)) - 0.04599*(J(J+1))**2
C + 0.000052*(J(J+1))**3 CM-1,
C IS A SPECIAL CASE OF THE FORMULA USED HERE, OBTAINED BY
C SETTING V=0 (IE., VP = 1/2); THE VP TERMS IN THE VIBRATION-
C ROTATION CONTRIBUTIONS CORRECT THE INITIAL TERMS TO
C PRODUCE THE COEFFICIENTS IN THE COHEN AND BIRNBAUM FORMULA.
C FOR THE COLD ATMOSPHERES OF THE OUTER PLANETS, THE SIGNIF-
C ICANTLY POPULATED LEVELS HAVE V=0.
C
C***********************************************************************
C
IMPLICIT REAL*8 (A-H,O-Z)
REAL*8 J,JP
C
VP = V + 0.5
JP = J + 1.0
E = 4400.39*VP - 120.815*VP**2 + 0.7242*VP**3 +
A (60.841 - 3.0177*VP + 0.0286*VP**2)*J*JP -
B (0.04684 - 0.00171*VP + 3.1D-05*VP**2)*J**2*JP**2 +
C 5.2D-05*J**3*JP**3 - 2170.08
H2ENER = E
RETURN
END
SUBROUTINE IKSUM(K,EK,NT,T,IK)
IMPLICIT REAL*8 (A-H,O-Z)
REAL*8 IK
DOUBLE PRECISION TSUM
DIMENSION T(250),IK(250),TSUM(250),EKT(250)
DATA PI/3.14159/,X1/0.1/,DX/0.05/,NMAX/199/
C
DO 10 I = 1,NT
EKT(I) = -4.0*EK/T(I)
10 TSUM(I) = 0.D0
C
DO 100 N = 1,NMAX
X = X1 + FLOAT(N-1)*DX
XKR = X**(2-K)
XDIF = X**(-6)
XDIF = XDIF*(XDIF - 1.)
DO 50 I = 1,NT
ARG = EKT(I) * XDIF
IF (ARG .LT. -20.0) GO TO 50
TNOW = XKR * EXP(ARG)
TSUM(I) = TSUM(I) + DBLE(TNOW)
50 CONTINUE
100 CONTINUE
C
COEF = 4.0*PI*DX
DO 125 I = 1,NT
125 IK(I) = COEF*SNGL(TSUM(I))
C
RETURN
END