SUBROUTINE RATEKIN INCLUDE 'PARAMETERS.INC' INCLUDE 'PLANET.INC' CHARACTER*8 AREACT(IB,6) COMMON & /SWITCHES/ ZLAT , DELTAT , DELT1 , DAYSTART , & TIMESTART , TIMESTOP , OUTTIME , & MODE , NTIME , NITER , IDIFF , istepdump & /BLOCK1/ NDISSC & /ONED15/ THOUR & /ONED32/ SUNRISE(NALT) , SUNSET(NALT) & /BLOCK2/ A(IB) , B(IB) , C(IB) & /BLOCK7/ RATE(IB,NALT) & /ONED28/ DEN(NALT) & /BLOCK11/ T(NALT) & /BLOCK18/ CONC(MB,NALT) & /BLOCK38/ NREACT & /ONED20/ RCONST(IB,NALT) C###################################################################### DO K = 1 , NALT C CALCULATE ALTITUDE/LATITUDE-DEPENDENT REACTION CONSTANTS C THE FIRST NDISSC REACTIONS ARE PHOTODISSOCIATIVE REACTIONS. C THEIR J VALUES ARE COMPUTED USING THE RESULTS FROM SUBROUTINE PHOTO IF( (MODE .NE. 0) .AND. ( (THOUR .GT. SUNSET(K)) .OR. & (THOUR .LT. SUNRISE(K)) ) ) THEN DO I = 1 , NDISSC RATE(I,K) = 0.0D0 ! IN DARKNESS J VALUES ARE ZERO ENDDO ELSE CALL JVALUE(K) ENDIF C CALCULATE CHEMICAL REACTION RATES C THE REACTION CONSTANT IS A STANDARD ARRHENIUS TYPE EXPRESSION: C K = A*(T**B)*EXP(C/T) N1 = NDISSC + 1 TEMP = T(K) DO I = N1 , NREACT RATE(I,K) = A(I)*TEMP**B(I)*DEXP(C(I)/TEMP) ENDDO C INCORPORATE CORRECTIONS DUE TO LYMAN ALPHA, SCHUMANN-RUNGE BANDS, ETC. CALL EARTH(K) DO I = 1 , NREACT RCONST(I,K) = RATE(I,K) ! ADJUSTED RATE CONSTANTS ENDDO ENDDO RETURN END C------------------------------------------------------------------------------ C------------------------------------------------------------------------------ C------------------------------------------------------------------------------ SUBROUTINE PMINUSL INCLUDE 'PARAMETERS.INC' INCLUDE 'PLANET.INC' CHARACTER*8 AREACT(IB,6) COMMON & /ONED17/ NTRANS & /ONED18/ NFORM(MB) , NDEST(MB) , & IFORM(IB,MB) , IDEST(IB,MB) & /ONED19/ KFORM(IB,MB) , KDEST(IB,MB) & /ONED27/ MOL(6,IB) , KOEFF(6,IB) & /BLOCK3/ NOREAC(IB) & /ONED20/ RCONST(IB,NALT) & /BLOCK7/ RATE(IB,NALT) & /BLOCK18/ CONC(MB,NALT) & /ONED40/ MTRANS(MB) & /BLOCK38/ NREACT & /ONED22/ RNET(MB,NALT) & /DEBUG/ RPROD(MB,NALT) & /ONED59/ RLOSS(MB,NALT) C CALCULATE PRODUCTION/LOSS RATES C FOR A REACTION, A + B + C = D + F + G, UP TO A SIGN, THE RATE OF C PRODUCTION/LOSS FOR ANY SPECIES IN THE REACTION IS GIVEN BY C K*[A]*[B]*[C], WHERE K IS THE REACTION CONSTANT AND [X] REPRESENTS C THE CONCENTRATION OF SPECIES X. C THE EFFECTIVE FORMATION/DESTRUCTION COEFFICIENT IS USED TO ACCOUNT C FOR THE FACT THAT SOME REACTIONS PRODUCE/DESTROY MORE THAN ONE C MOLECULE OF A GIVEN SPECIES OR INVOLVE A CATALYST c$doacross c$& local(k,i,r,nrct,j,c0,c1,kk, c$& m,mvary,rp,nf,ifrm, c$& rl,nd,idst) , c$& share(rconst,noreac,conc,mol,koeff,rate, c$& mtrans,nform,iform,kform,rprod, c$& ndest,kdest,rloss,rnet) DO K = 1 , NALT DO I = 1 , NREACT R = RCONST(I,K) NRCT = NOREAC(I) DO J = 1 , NRCT C0 = CONC(MOL(J,I),K) C1 = C0 DO KK = 1 , KOEFF(J,I)-1 ! MORE EFFICIENT WAY OF DOING CONC**KOEFF C1 = C1*C0 ENDDO R = R*C1 ENDDO RATE(I,K) = R ENDDO DO M = 1 , NTRANS MVARY = MTRANS(M) C CALCULATE TOTAL PRODUCTION RATE FOR THE M-TH TRANSPORTED SPECIES RP = 0.0D0 NF = NFORM(MVARY) DO I = NF , 1 , -1 IFRM = IFORM(I,MVARY) RP = RP + RATE(IFRM,K)*KFORM(IFRM,MVARY) ENDDO RPROD(MVARY,K) = RP C CALCULATE TOTAL DESTRUCTION RATE FOR THE M-TH TRANSPORTED SPECIES RL = 0.0D0 ND = NDEST(MVARY) DO I = ND , 1 , -1 IDST = IDEST(I,MVARY) RL = RL + RATE(IDST,K)*KDEST(IDST,MVARY) ENDDO RLOSS(MVARY,K) = RL C CALCULATE P-L RNET(M,K) = RP - RL ENDDO ENDDO ! END LOOP OVER K RETURN END C------------------------------------------------------------------------------ C------------------------------------------------------------------------------ C------------------------------------------------------------------------------ DOUBLE PRECISION FUNCTION CHAP(COSTOP,SINTOP,SCALE) INCLUDE 'PARAMETERS.INC' INCLUDE 'PLANET.INC' COMMON & /BLOCK10/ Z(NALT) & /BLOCK27/ PI C THIS FUNCTION COMPUTES THE CHAPMAN FUNCTION CORRECTION FOR C SPHERICAL ATMOSPHERES. THE CHAPMAN FUNCTION IS COMPUTED USING THE C APPROXIMATIONS FOUND IN "JOURNAL OF GEOPHYSICAL RESEARCH" VOLUME 77, C NUMBER 19 , PAGES 3592-3597 (1972) DATA A1,A2,A3,A4/1.0606963d0,0.55643831d0,1.0619896d0,1.7245609d0/ DATA B1,B2/0.56498823d0,0.06651874d0/ X = (RADIUS + Z(NALT))/(SCALE*1.0d-5) Y = SQRT( 0.5d0*X*(COSTOP)**2 ) C THE COMPLEMENTARY ERROR FUNCTION IS APPROXIMATED BY A PADE' C APPROXIMATION. THIS GREATLY SIMPLIFIES THE FORMULA. IF( Y .LE. 8.0d0 ) THEN CHAP = SQRT(PI*X/2.0d0)*( (A1+A2*Y)/(A3+A4*Y+Y**2) ) ELSE CHAP = SQRT(PI*X/2.0d0)*( B1/(B2 + Y) ) ENDIF IF( COSTOP .LT. 0.0d0) THEN POWER = X*(1.0d0-SINTOP) CHAP = SQRT(2.0d0*PI*X*SINTOP)*EXP( POWER ) - CHAP ENDIF RETURN END C------------------------------------------------------------------------------ C------------------------------------------------------------------------------ C------------------------------------------------------------------------------ DOUBLE PRECISION FUNCTION E(N,X) DOUBLE PRECISION X C THIS FUNCTION COMPUTES THE EXPONENTIAL INTEGRAL OF ORDER N. THIS ROUTINE HAS C BEEN ADAPTED FROM "Numerical Recipes in FORTRAN" WHERE IT IS CALLED EXPINT. C SMALL MODIFICATIONS HAVE BEEN MADE SINCE WE ONLY NEED IT FOR N = 2 AND 3 PARAMETER (MAXIT = 100 , EPS = 1.0d-1 , FPMIN = 1.0d-30 ) PARAMETER (EULER = 0.5772156649d0 ) NM1 = N - 1 IF ( X .EQ. 0.0d0 ) THEN E = 1.0d0/NM1 ! SPECIAL CASE RETURN ENDIF IF (X .GT. 1.0d0 ) THEN B = X + N C = 1.0d0/FPMIN D = 1.0d0/B H = D DO I = 1 , MAXIT A = -I*(NM1 + I) B = B + 2.0d0 D = 1.0d0/(A*D + B) C = B + A/C DEL = C*D H = H*DEL IF( ABS(DEL-1.0d0) .LT. EPS ) THEN E = H*DEXP(-X) RETURN ENDIF ENDDO ELSE E = 1.0d0/NM1 FACT = 1.0d0 DO I = 1 , MAXIT FACT = -FACT*X/I IF( I .NE. NM1 ) THEN DEL = -FACT/(I - NM1) ELSE PSI = -EULER DO II = 1 , NM1 PSI = PSI + 1.0d0/II ENDDO DEL = FACT*(-DLOG(X) + PSI) ENDIF E = E + DEL IF( ABS(DEL) .LT. ABS(E)*EPS ) RETURN ENDDO ENDIF RETURN END C------------------------------------------------------------------------------ C------------------------------------------------------------------------------ C------------------------------------------------------------------------------ SUBROUTINE INPUT INCLUDE 'PLANET.INC' INCLUDE 'PARAMETERS.INC' CHARACTER*8 AMOL(MB) , AREACT(IB,6) , STOP , ATRANS(MB) , & ATAU(MB) , APHOTO(MB) , ABNDRY , APCE(MB) , H2O CHARACTER SYMBOL(IB,6) , PLUS , EQUAL , ACOEFF(IB,6) , BLANK CHARACTER*2 ATOM(MB),QUIT LOGICAL FINISH DIMENSION EMASS(10) , NOPROD(IB) COMMON & /SWITCHES/ ZLAT , DELTAT , DELT1 , DAYSTART , & TIMESTART , TIMESTOP , OUTTIME , & MODE , NTIME , NITER , IDIFF , istepdump & /BLOCK1/ NDISSC & /ONED12/ ILOWER(MB) , IUPPER(MB) , & XLOWER(MB) , XUPPER(MB) & /ONED17/ NTRANS & /ONED18/ NFORM(MB) , NDEST(MB) , & IFORM(IB,MB) , IDEST(IB,MB) & /BLOCK2/ A(IB) , B(IB) , C(IB) & /ONED40/ MTRANS(MB) & /ONED19/ KFORM(IB,MB) , KDEST(IB,MB) & /ONED27/ MOL(6,IB) , KOEFF(6,IB) & /ONED42/ NMOL & /BLOCK3/ NOREAC(IB) & /ONED43/ AMOL & /BLOCK4/ ATOM , ATRANS & /BLOCK5/ ACOEFF , SYMBOL , APHOTO , APCE , NPHOTO & /BLOCK10/ Z(NALT) & /BLOCK11/ T(NALT) & /BLOCK13/ PRESS(NALT) & /BLOCK14/ MTAU(MB) & /BLOCK15/ NTAU & /BLOCK18/ CONC(MB,NALT) & /BLOCK28/ SCALE(MB) & /BLOCK34/ ATAU & /BLOCK37/ AREACT & /BLOCK38/ NREACT & /BLOCK39/ NX(10,MB) , HF(MB) , DELH(IB) , & NATOM , MAXFRM , MAXDES , MFORM , MDEST & /ONED26/ XMASS(MB) ! USED IN SUBROUTINE COEFF & /ONED29/ EDDY(NALT) ! USED IN SUBROUTINE COEFF & /ONED53/ MPCE(MB) & /ONED54/ NPCE & /ONED55/ MH2O , MH2OTRAN ! H2O INDEX IN SPECIES LISTING AND IN ! TRANSPORTED SPECIES LISTING dimension temp(mb,nalt) DATA STOP/'STOP '/ , BLANK/' '/ , PLUS/'+'/ , EQUAL/'='/ , & QUIT/'Q '/ , H2O/'H2O '/ C PROGRAM IS SET UP TO HANDLE 10 ATOMS, 100 SPECIES, AND 300 REACTIONS C READ IN SWITCHES READ(4,15) MODE READ(4,15) NTIME READ(4,15) NITER READ(4,15) IDIFF READ(4,25) ZLAT READ(4,25) DELTAT READ(4,25) DELT1 READ(4,25) DAYSTART READ(4,25) TIMESTART READ(4,25) TIMESTOP READ(4,25) OUTTIME read(4,15) istepdump 15 FORMAT(I4) 25 FORMAT(F10.1) C READ IN ELEMENT LIST FROM JPL90 FILE N = 0 FINISH = .FALSE. DO WHILE(FINISH .NE. .TRUE.) N = N + 1 READ(23,200) ATOM(N) , EMASS(N) IF(ATOM(N) .EQ. QUIT) FINISH = .TRUE. ENDDO NATOM = N - 1 C READ IN MOLECULE LIST AND INITIALIZE NFORM(M) AND NDEST(M) M = 0 FINISH = .FALSE. DO WHILE (FINISH .NE. .TRUE.) M = M + 1 NFORM(M) = 0 NDEST(M) = 0 READ (23,210) AMOL(M) , (NX(J,M),J=1,10) print*,'input',m,amol(m) IF(AMOL(M) .EQ. STOP) FINISH = .TRUE. IF(AMOL(M) .EQ. H2O ) MH2O = M ENDDO NMOL = M - 1 C IDENTIFY THE TRANSPORTED SPECIES M = 0 MH2OTRAN = 0 FINISH = .FALSE. DO WHILE (FINISH .NE. .TRUE.) M = M + 1 READ(23,220) ATRANS(M) IF(ATRANS(M) .EQ. STOP) FINISH = .TRUE. IF(ATRANS(M) .EQ. H2O ) MH2OTRAN = M ENDDO NTRANS = M - 1 DO N = 1 , NTRANS M = 1 DO WHILE(ATRANS(N) .NE. AMOL(M)) M = M + 1 ENDDO MTRANS(N) = M ENDDO C IDENTIFY THE SPECIES KEPT IN PHOTOCHEMICAL EQUILIBRIUM (PCE) M = 0 FINISH = .FALSE. DO WHILE (FINISH .NE. .TRUE.) M = M + 1 READ(23,220) APCE(M) IF(APCE(M) .EQ. STOP) FINISH = .TRUE. ENDDO NPCE = M - 1 DO N = 1 , NPCE M = 1 DO WHILE(APCE(N) .NE. AMOL(M)) M = M + 1 ENDDO MPCE(N) = M ENDDO C IDENTIFY PHOTODISSOCIATIVE SPECIES WHICH CONTRIBUTE TO THE C ATMOSPHERIC OPACITY. THAT IS, WHICH SPECIES ABSORB A SIGNIFICANT C PORTION OF THE SOLAR FLUX N = 0 FINISH = .FALSE. DO WHILE (FINISH .NE. .TRUE.) N = N + 1 READ(23,220) ATAU(N) print*,'opacity',n,atau(N) IF(ATAU(N) .EQ. STOP) FINISH = .TRUE. ENDDO NTAU = N - 1 DO N = 1 , NTAU M = 1 DO WHILE(ATAU(N) .NE. AMOL(M)) M = M + 1 ENDDO MTAU(N) = M C COMPUTE MOLECULAR MASS FOR USE IN OPACITY CALCULATIONS AMU = 0.0d0 DO J = 1 , 7 AMU = AMU + NX(J,M)*EMASS(J) ENDDO SCALE(M) = BOLTZ/(AMU*GRAV) ENDDO C COMPUTE MASS OF EACH MOLECULE. THIS INFORMATION IS NEEDED IN SUBROUTINE COEFF DO M = 1 , NMOL AMU = 0.0d0 DO J = 1 , 7 AMU = AMU + NX(J,M)*EMASS(J) ENDDO XMASS(M) = AMU ENDDO C READ IN REACTION LIST C THE REACTIONS ARE INPUT IN THE FORM: C KOEFF(1,I)*AREACT(I,1)*SYMBOL(I,1),KOEFF(2,I)*AREACT(I,2)*SYMBOL(I,2) C ... WHERE KOEFF(1,I) IS THE COEFFICIENT OF THE FIRST REACTANT, C AREACT(I,1) IS THE FIRST REACTANT, AND SYMBOL(I,1) IS EITHER C "+" , "=" , OR " ". FOR EXAMPLE, THE REACTION: Cl + OClO = 2ClO HAS C KOEFF(1,I) = 1, KOEFF(2,I) = 1, KOEFF(3,I) = 2 C AREACT(I,1) = Cl, AREACT(I,2) = OClO, AREACT(I,3) = ClO C SYMBOL(I,1) = "+", SYMBOL(I,2) = "=", AND SYMBOL(I,3) = " ". C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C NOTE: ****THE REACTIONS MUST BE INPUT IN THE FORM 2X , NOT X + X.***** C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! I = 0 FINISH = .FALSE. DO WHILE(FINISH .NE. .TRUE.) I = I + 1 READ (23,230) (KOEFF(J,I),AREACT(I,J),SYMBOL(I,J),J=1,6),A(I), & B(I) , C(I) IF(AREACT(I,1) .EQ. STOP) FINISH = .TRUE. ENDDO NREACT = I - 1 C READ IN BOUNDRY CONDITIONS DO M = 1 , NMOL ILOWER(M) = 4 ! INTIALIZE DEFAULT VALUES XLOWER(M) = 0.0d0 IUPPER(M) = 4 XUPPER(M) = 0.0d0 ENDDO DO WHILE(ABNDRY .NE. STOP) READ(23,235) ABNDRY , ILOW , XLOW , IUP , XUP DO M = 1 , NMOL ! MATCH UP B.C. TO APPROPRIATE SPECIES IF(AMOL(M) .EQ. ABNDRY) THEN ILOWER(M) = ILOW XLOWER(M) = XLOW IUPPER(M) = IUP XUPPER(M) = XUP ENDIF ENDDO ENDDO C THE REACTION LIST CONTAINS BLANKS FOR THE COEFFICIENT OF MOST C REACTANTS. FOR EXAMPLE: NO3 = NO2 + O INSTEAD OF 1NO3 = 1NO2 + 1O. C THIS IS DONE IN ORDER TO MAKE THE REACTIONS LESS CONFUSING AND THUS C LESS PRONE TO MISTAKES. C FOR OBVIOUS REASONS, WE WISH TO READ THE COEFFICIENTS IN AS INTEGERS. C AS AN INTEGER, THE BLANK IS READ IN AS A ZERO. THUS WE NEED TO CONVERT C THEM BACK TO ONES. SIMILARLY, FOR OUTPUT PURPOSES WE WISH TO CONVERT C COEFFICIENTS OF ONE BACK INTO BLANKS. FOR EXAMPLE, O + O2 INSTEAD OF C 1O + 1O2. THIS IS DONE BY ASSIGNING A NEW CHARACTER VARIABLE, ACOEFF, C WHICH IS USED FOR OUTPUT PURPOSES ONLY. INTEGER VALUES OF KOEFF ARE C CONVERTED TO CHARACTER VALUES BY WAY OF THE FORTRAN FUNCTION 'CHAR'. C 'CHAR' CONVERTS AN ASCII INTEGER VALUE INTO A CHARACTER VALUE. FOR C SOME STRANGE REASON, THE ASCII CHARACTER FOR 48 IS '0', 49 IS '1', ETC. DO I = 1, NREACT DO J = 1 , 6 KOEFF(J,I) = MAX0(1,KOEFF(J,I)) ACOEFF(I,J) = CHAR( 48 + KOEFF(J,I)) IF(KOEFF(J,I).LE.1) ACOEFF(I,J) = ' ' MOL(J,I) = 0 ! INTIALIZE FOR USE IN A LATER LOOP ENDDO C CALCULATE NUMBER OF REACTANTS IN ITH REACTION J = 1 DO WHILE(SYMBOL(I,J) .NE. EQUAL) J = J + 1 ENDDO NOREAC(I) = J NRCT = NOREAC(I) C CALCULATE NUMBER OF PRODUCTS IN ITH REACTION J = NRCT + 1 DO WHILE(SYMBOL(I,J) .NE. BLANK) J = J + 1 ENDDO NOPROD(I) = J - NRCT NPRD = NOPROD(I) C THE FOLLOWING SETS OF LOOPS COUNT UP THE NUMBER OF REACTIONS IN WHICH C THE MTH SPECIES IS A REACTANT/PRODUCT. IF THERE IS NO MATCH WITH ANY C INPUT SPECIES, THEN AN ERROR MESSAGE IS PRINTED OUT. IF A MATCH IS C FOUND IN THE INNER LOOP OVER THE INPUT SPECIES, THEN CONTROL JUMPS TO C CHECKING A NEW REACTANT/PRODUCT. C CHECK REACTANTS DO J = 1 ,NRCT DO M = 1 , NMOL IF(AREACT(I,J) .EQ. AMOL(M)) THEN NDEST(M) = NDEST(M) + 1 NDE = NDEST(M) IDEST(NDE,M) = I MOL(J,I) = M ! THE JTH TERM IS SPECIES M GOTO 10 ENDIF ENDDO C PRINT OUT ERROR MESSAGE, IF NO MATCH IS FOUND IF(MOL(J,I) .EQ. 0) THEN WRITE(25,250) I, & (ACOEFF(I,L),AREACT(I,L),SYMBOL(I,L), L=1,6) & , AREACT(I,J) MOL(J,I) = NMOL-2 ! SET UNKNOWNS TO "PROD" ( CONC = 0 ) A(I) = 0.0d0 ! UNBALANCED CHEMISTRY - SET RATE TO ZERO ENDIF 10 ENDDO C CHECK PRODUCTS DO J = 1 , NPRD DO M = 1 , NMOL IF(AREACT(I,NRCT+J) .EQ. AMOL(M)) THEN NFORM(M) = NFORM(M) + 1 NFO = NFORM(M) IFORM(NFO,M) = I MOL(NRCT+J,I) = M GOTO 20 ENDIF ENDDO C PRINT OUT ERROR MESSAGE, IF NO MATCH IS FOUND IF(MOL(NRCT+J,I).EQ.0) THEN WRITE(25,250) I, & (ACOEFF(I,L),AREACT(I,L),SYMBOL(I,L), L=1,6) & , AREACT(I,NRCT+J) MOL(NRCT+J,I) = NMOL-2 ! SET UNKNOWNS TO "PROD" ( CONC = 0 ) A(I) = 0.0d0 ! UNBALANCED CHEMISTRY - SET RATE TO ZERO ENDIF 20 ENDDO C CALCULATE REACTION ENTHALPY. C HFR = 0.0d0 C HFP = 0.0d0 C DO J = 1 , NRCT C HFR = HFR + KOEFF(J,I)*HF(MOL(J,I)) C ENDDO C DO J = 1 , NPRD C HFP = HFP + KOEFF(NRCT+J,I)*HF(MOL(NRCT+J,I)) C ENDDO C DELH(I) = HFP - HFR C PERFORM MASS BALANCE CHECK DO N = 1 , NATOM NELEM = 0 DO J = 1 , NPRD NELEM = NELEM + KOEFF(NRCT+J,I)*NX(N,MOL(NRCT+J,I)) ENDDO DO J = 1 , NRCT NELEM = NELEM - KOEFF(J,I)*NX(N,MOL(J,I)) ENDDO IF( NELEM.NE.0 ) WRITE(25,240) I, & (ACOEFF(I,L),AREACT(I,L),SYMBOL(I,L), & L=1,6) , ATOM(N) ENDDO ENDDO ! END LOOP OVER I C DETERMINE WHICH REACTIONS ARE PHOTODISSOCIATIVE REACTIONS C PHOTODISSOCIATIVE REACTIONS HAVE ONLY ONE REACTANT. THUS, WE ONLY C NEED TO CHECK THE NUMBER OF REACTANTS IN A REACTION IN ORDER TO C DETERMINE IF IT IS A PHOTODISSOCIATIVE REACTION. NOTICE THAT WE ALSO C NEED TO CHECK THE COEFFICIENTS OF THE REACTANTS. FOR EXAMPLE: 2O = O2 C REALLY INVOLVES TWO REACTANTS AND IS NOT A PHOTODISSOCIATIVE REACTION. C WE TAKE INTO ACCOUNT THE FACT THAT SOME SPECIES PHOTODISSOCIATE IN C MORE THAN ONE WAY. THE NUMBER OF PHOTODISSOCIATING SPECIES WILL BE C GIVEN BY NPHOTO WHILE THE NUMBER OF PHOTODISSOCIATING REACTIONS WILL C BE GIVEN BY NDISSC. C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C THE FIRST SET OF REACTIONS IN THE INPUT FILE WILL BE PHOTODISSOCIATIVE C REACTIONS C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! I = 1 L = 0 DO WHILE( (NOREAC(I).EQ.1) .AND. (KOEFF(1,I).EQ.1) ) K = 1 DO WHILE (AREACT(K,1) .NE. AREACT(I,1) ) K = K + 1 ENDDO IF( K .EQ. I) THEN L = L + 1 APHOTO(L) = AREACT(I,1) ENDIF I = I + 1 ENDDO NPHOTO = L NDISSC = I - 1 C CALCULATE EFFECTIVE (I.E. NET) FORMATION AND DESTRUCTION COEFFICIENTS C FOR EACH SPECIES IN REACTION. THIS IS NECESSARY BECAUSE SOME REACTIONS C E.G. 2O + O2 = O + O3 INVOLVE THE SAME SPECIES ON BOTH SIDES OF THE C EQUATION. THUS, IN THE ABOVE REACTION, THE EFFECTIVE DESTRUCTION C COFFECIENT FOR OXYGEN IS ONE - NOT TWO - BECAUSE THE NET LOSS IS C ONLY ONE OXYGEN ATOM. IN THE ABOVE REACTION OXYGEN IS NOT REALLY C FORMED AND MUST THEREFORE BE DELETED FROM THE FORMATION LIST. THESE C INSTANCES CAN BE EASILY IDENTIFIED SINCE THEIR EFFECTIVE C FORMATION/DESTRUCTION COEFFICIENTS WILL BE NEGATIVE. DO M = 1 , NMOL I = 1 DO WHILE (I .LE. NFORM(M)) IFRM = IFORM(I,M) NRCT = NOREAC(IFRM) KFORM(IFRM,M) = 0 C DETERMINE IF SPECIES M IS ALSO A REACTANT IN THE IFRM-TH REACTION C (BY DEFINITION OF IFORM IT IS ALREADY KNOWN TO BE A PRODUCT) DO J = 1 , NRCT IF ( MOL(J,IFRM).EQ.M ) THEN KFORM(IFRM,M) = - KOEFF(J,IFRM) GOTO 30 ENDIF ENDDO C DETERMINE COEFFICIENT OF MTH SPECIES IN THE IFRM-TH REACTION. NOTICE C THAT BY DEFINITION OF IFRM, THE MTH SPECIES MUST OCCUR IN THE PRODUCT C SIDE OF THE REACTION 30 J = NRCT + 1 DO WHILE( MOL(J,IFRM) .NE. M) J = J + 1 ENDDO KFORM(IFRM,M) = KFORM(IFRM,M) + KOEFF(J,IFRM) C ADJUST NFORM AND IFORM LIST TO ACCOUNT FOR REACTIONS THAT INVOLVE C THE SAME REACTANT ON BOTH SIDES OF THE EQUATION IF(KFORM(IFRM,M).LE.0) THEN NFORM(M) = NFORM(M) - 1 NF = NFORM(M) DO K = I , NF IFORM(K,M) = IFORM(K+1,M) ENDDO I = I - 1 ENDIF I = I + 1 ENDDO ! END LOOP OVER I C PERFORM A SIMILAR CALCULATION FOR THE EFFECTIVE DESTRUCTION C COEFFICIENTS I = 1 DO WHILE (I .LE. NDEST(M)) IDST = IDEST(I,M) NRCT = NOREAC(IDST) NPRD = NOPROD(IDST) KDEST(IDST,M) = 0 C DETERMINE IF THE MTH SPECIES APPEARS AS A PRODUCT IN THE IDST-TH C REACTION N1 = NRCT + 1 N2 = NRCT + NPRD DO J = N1 , N2 IF( MOL(J,IDST) .EQ. M) THEN KDEST(IDST,M) = -KOEFF(J,IDST) GOTO 40 ENDIF ENDDO C DETERMINE WHERE THE MTH SPECIES APPEARS AS A REACTANT IN THE IDST-TH C REACTION 40 J = 1 DO WHILE ( MOL(J,IDST) .NE. M) J = J + 1 ENDDO KDEST(IDST,M) = KDEST(IDST,M) + KOEFF(J,IDST) C ADJUST NDEST AND IDEST LIST TO ACCOUNT FOR REACTIONS THAT INVOLVE C THE SAME REACTANT ON BOTH SIDES OF THE EQUATION IF(KDEST(IDST,M).LE.0) THEN NDEST(M) = NDEST(M) - 1 ND = NDEST(M) DO K = I , ND IDEST(K,M) = IDEST(K+1,M) ENDDO I = I - 1 ENDIF I = I + 1 ENDDO ! END LOOP OVER I ENDDO ! END LOOP OVER M C DETERMINE LARGEST VALUES FOR NFORM AND NDEST MAXFRM = 0 MAXDES = 0 DO M = 1 , NMOL IF( NFORM(M) .GT. MAXFRM) THEN MFORM = M MAXFRM = NFORM(M) ENDIF IF( NDEST(M) .GT. MAXDES) THEN MDEST = M MAXDES = NDEST(M) ENDIF ENDDO C############################################################################## C READ IN INITIAL TEMPERATURE AND PRESSURE PROFILES C READ(1,100) ! SKIP OVER HEADER INFORMATION C DO K = 1 , NALT C READ(1,110) Z(K) , T(K) , PRESS(K) C ENDDO C READ IN EDDY DIFFUSION COEFFICIENT PRINT*,'GETTING READY TO READ IN Z,EDDY',NALT DO K = 1 , NALT READ(31,120) Z(K) , EDDY(K) print*,k,z(k),eddy(k) CC EDDY(K) = EDDY(K) !what the heck was this for? ENDDO c WRITE(6,*) Z c WRITE(6,*) NALT C READ IN INITIAL CONCENTRATION PROFILES NSTOP = NMOL/10 IF( MOD(NMOL,10) .NE. 0 ) NSTOP = NSTOP + 1 DO N = 1 , NSTOP READ(2,130) ! SKIP OVER HEADER INFORMATION DO K = 1 , NALT c print*,K READ(2,140) ( CONC( (N-1)*10+M ,K) , M = 1 , 10) ENDDO ENDDO IFLAG = 0 C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 100 FORMAT(//) 110 FORMAT( 11X , F5.1 , 6X , F5.1 , 17X , E10.3 ) 120 FORMAT( 1X , F5.1 , 3X , E9.3 ) 130 FORMAT( ) 140 FORMAT( 19X , 10(E9.3 , 2X) ) C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 200 FORMAT (A2,1X,F5.2) C210 FORMAT (A8,2X,10(I2,3X),4X,F6.1) 210 FORMAT (A8,2X,10(I2,3X)) 220 FORMAT(A8) 230 FORMAT (6(I1,A8,A1),E9.2,F5.1,F7.0) 235 FORMAT(A8,1X,I1,1X,E9.2,2X,I1,1X,E9.2) 240 FORMAT (1X,I3,'.',2X, 6(A1,A8,A1),/,3X, & 'MASS BALANCE ERROR FOR',A8/) 250 FORMAT (1X,I3,'.',2X, 6(A1,A8,A1),/,3X,'CANNOT FIND ',A8, & ' IN MOLECULE INPUT LIST'/) C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ RETURN END C------------------------------------------------------------------------------ C------------------------------------------------------------------------------ C------------------------------------------------------------------------------ SUBROUTINE PHOTO INCLUDE 'PARAMETERS.INC' CHARACTER*8 AREACT(IB ,6) , ADISSC(IB,4) , ACROSS , & ATAU(MB) , APHOTO(MB) , APCE(MB) , STOP , BLANK CHARACTER ACOEFF(IB,6) , SYMBOL(IB,6) DIMENSION ICROSS(MB) , ISIGMA(IB) , FRACT(IB,176) DIMENSION WTEMP(176) , CTEMP(MB,176) , STEMP(176) , SFTEMP(176) COMMON & /BLOCK1/ NDISSC & /BLOCK5/ ACOEFF , SYMBOL , APHOTO , APCE , NPHOTO & /BLOCK6/ MCROSS(MB) & /BLOCK8/ IPHOTO(IB) , LAST(IB) , & LOW(IB) , SIGMA(NWAVE,IB) & /BLOCK9/ WAVE(NWAVE) & /BLOCK12/ CROSS(NWAVE,MB) & /BLOCK15/ NTAU & /BLOCK29/ FLUX(NWAVE) & /BLOCK34/ ATAU & /BLOCK37/ AREACT LOGICAL FINISH , DEPTH(MB) DATA STOP/'STOP '/ C THIS SUBROUTINE CALCULATES THE PHOTODISSOCIATION RATES FOR THE C PHOTODISSOCIATIVE REACTIONS WHICH APPEAR IN THE MAIN PROGRAM. C THIS SUBROUTINE FIRST READS IN A PHOTODISSOCIATING/ABSORBING SPECIES, C THE NUMBER OF ASSOCIATED BRANCHES, AND ITS CROSS-SECTION AS A C FUNCTION OF WAVELENGTH. INDIVIDUAL REACTIONS CORRESPONDING TO THE C PHOTODISSOCIATIVE SPECIES ARE THEN READ IN ALONG WITH THEIR C FRACTIONAL CONTRIBUTION TO THE TOTAL CROSS SECTION. C THE FLUX AND CROSXX FILES CONTAIN INFORMATION FROM ZERO TO 8,000 C ANGSTROMS. SINCE ONLY WAVELENGTHS GREATER THAN LYMAN ALPHA C CAN PENETRATE INTO THE MIDDLE ATMOSPHERE, WE WILL PASS THAT PORTION C ONLY BACK TO THE MAIN ROUTINE. WE SET UP TEMPORARY ARRAYS - CTEMP, C STEMP, WTEMP, AND SFTEMP WHICH STORE ALL OF THE DATA FROM THE CROSXX FILE. C THE CROSS SECTIONS AND FLUXES WHICH ARE OF INTEREST TO US ARE THEN C CREATED FROM THESE TEMPORARY ARRAYS. C INITIALIZE OPACITY FLAG DO N = 1 , NTAU DEPTH(N) = .FALSE. ENDDO C INITIALIZE STARTING VALUES FINISH = .FALSE. M = 0 I = 0 NCROSS = 0 NSIGMA = 0 C READ IN VALUES OF WAVELENGTH BINS DO LAMBDA = 1 , 22 READ(24,300) (WTEMP(L+(LAMBDA-1)*8) , L = 1 , 8) ENDDO READ(24,310) C READ IN VALUES OF THE SOLAR FLUX DO LAMBDA = 1 , 22 READ(24,320) (SFTEMP(L+(LAMBDA-1)*8) , L = 1 , 8) ENDDO DO LAMBDA = 78 , 176 WAVE(LAMBDA-77) = WTEMP(LAMBDA) FLUX(LAMBDA-77) = SFTEMP(LAMBDA) ENDDO C READ IN CROSS-SECTION SPECIES (ACROSS) , UPPER AND LOWER WAVELENGTH C LIMITS (L1 , L2) AND THE NUMBER OF REACTIONS ASSOCIATED WITH EACH C ABSORBING SPECIES (NBRNCH) DO WHILE (FINISH .NE. .TRUE.) M = M + 1 READ(27,330) ACROSS , NBRNCH , L1 , L2 IF(ACROSS .EQ. STOP) THEN FINISH = .TRUE. GOTO 20 ENDIF LEND = ((L2 - L1 + 1)/8) LASTROW = MOD(L2 - L1 + 1 ,8) C THE CROSS SECTIONS ARE READ IN ROW BY ROW. C THE DATA IS STORED IN AN ARRAY CONSISTING OF 8 COLUMNS. SINCE C FORTRAN THROWS THE REMAINDER AWAY, (L2-L1)/8 GIVES THE NUMBER C OF COMPLETE ROWS. C LEFT IS THE INDEX CORRESPONDING TO THE LEFTMOST VALUE OF THE C LAMBDA-TH ROW. DO LAMBDA = 1 , 176 CTEMP(M,LAMBDA) = 0.0d0 ! ZERO OUT TEMPORARY ARRAY ENDDO DO LAMBDA = 1 , LEND LEFT = L1 + (LAMBDA-1)*8 - 1 READ(27,320)(CTEMP(M,LEFT + L),L=1,8) ENDDO IF( LASTROW .GE. 1 ) THEN LEFT = L1 + LEND*8 - 1 READ(27,320)(CTEMP(M,LEFT+L),L=1,LASTROW) ! READ IN REMAINING ROW ENDIF DO LAMBDA = 78 , 176 CROSS(LAMBDA-77,M) = CTEMP(M,LAMBDA) ENDDO C THIS SECTION COMPARES THE REACTIONS AGAINST PHOTODISSOCIATION AND C OPACITY LISTS FROM THE JPL INPUT FILE. SEPARATE CHECKS ARE C PERFORMED ON THE OPACITY LIST AND PHOTODISSOCIATION LISTS. THIS IS C NECESSARY BECAUSE SOME SPECIES HAVE A LISTING FOR ABSORPTION VALUES C AND NO CORRESPONDING SET OF PHOTODISSOCIATION BRANCHES C COMPARE WITH OPACITY LISTS. IF A MATCH IS FOUND, THEN SIGNIFY THIS C BY SETTING DEPTH = .TRUE. AFTER ALL SPECIES HAVE BEEN READ IN, C UNMATCHED VALUES WILL BE FLAGGED WITH DEPTH = .FALSE. DO N = 1 , NTAU IF(ACROSS .EQ. ATAU(N)) THEN DEPTH(N) = .TRUE. MCROSS(N) = M GOTO 10 ENDIF ENDDO C THIS SECTION READS IN EVERY REACTION ASSOCIATED WITH A C PHOTODISSOCIATING SPECIES. SOME OF THESE REACTIONS WILL NOT BE C RELEVANT TO THE MODEL BEING STUDIED. THE LIST OF REACTIONS WHICH C PAIR UP WITH THE JPL LIST IS COMMUNICATED TO THE MAIN PROGRAM C BY WAY OF THE FUNCTION IPHOTO(I). 10 DO N = 1 , NPHOTO IF(ACROSS .EQ. APHOTO(N)) THEN DO NBR = 1 , NBRNCH I = I + 1 READ(27,340)(ADISSC(I,J) ,J = 1,4) , LOW(I) , LAST(I) LEND = (( LAST(I) - LOW(I) + 1)/8) LASTROW = MOD(LAST(I) - LOW(I) + 1 ,8) DO LAMBDA = 1 , 176 STEMP(LAMBDA) = 0.0d0 !ZERO OUT TEMPORARY ARRAY ENDDO DO LAMBDA = 1 , LEND LEFT = LOW(I) + (LAMBDA-1)*8 - 1 READ(27,350)(FRACT(I , LEFT + L) , L = 1,8) DO L = 1 , 8 STEMP(LEFT+L) = FRACT(I,LEFT+L)*CTEMP(M,LEFT+L) ENDDO ENDDO IF( LASTROW .GE. 1 ) THEN LEFT = LOW(I) + LEND*8 - 1 READ(27,350)(FRACT(I,LEFT+L),L=1,LASTROW) ! READ IN LAST ROW DO L = 1 , LASTROW STEMP(LEFT+L) = FRACT(I,LEFT+L)*CTEMP(M,LEFT+L) ENDDO ENDIF DO LAMBDA = 78 , 176 SIGMA(LAMBDA-77,I) = STEMP(LAMBDA) ENDDO LOW(I) = MAX0(1,LOW(I)-77) ! RESET LIMITS IN TERMS OF LY-ALPHA LAST(I) = MAX0(0,LAST(I)-77) ENDDO GOTO 20 ENDIF ENDDO C IF THE MTH BLOCK OF REACTIONS DOES NOT CORRESPOND TO A PHOTODISSOCIATING C SPECIES USED IN THE JPL MODEL, THEN SKIP OVER THE REACTIONS (NBRNCH IN C NUMBER) AND CHECK A NEW BLOCK OF REACTIONS. DO N = 1 , NBRNCH READ(27,360) L1 , L2 LEND = ((L2 - L1)/8) + 1 DO L = 1 , LEND READ(27, 370) ENDDO ENDDO 20 ENDDO ! END INPUT C NSIGMA IS THE TOTAL NUMBER OF PHOTODISSOCIATING REACTIONS. NSIGMA = I C THIS SECTION OF SUBROUTINE PHOTO COMPARES INFORMATION READ IN C FROM THE CROSXX.INP FILE WITH CORRESPONDING INFORMATION FROM C THE JPL90.INP FILE. C WRITE ERROR MESSAGE IF THE CROSXX.INP FILE IS MISSING VALUES C FOR AN OPACITY SPECIES IDENTIFIED IN THE JPL90.INP FILE DO N = 1 , NTAU IF(DEPTH(N) .EQ. .FALSE.) WRITE(25,380) ATAU(N) ENDDO C COMPARE PHOTOTDISSOCIATING REACTIONS WITH THOSE FOUND IN THE JPL C INPUT LIST. IDENTIFYING THE CORRESPONDING REACTION NUMBER IS C COMPLICATED BY THE FACT THAT THE PRODUCTS MAY APPEAR IN A C DIFFERENT ORDER IN THE JPL LIST. FOR EXAMPLE, THE JPL LISTING C MIGHT READ A = B + C + D WHEREAS THE CROSS-SECTION LISTING C MIGHT APPEAR AS A = C + B + D C THE FOLLOWING LOOP PAIRS UP THE CROSS-SECTION INDEXING C WITH THE JPL INDEXING ( IPHOTO(JPL) ) DO JPL = 1 , NDISSC DO I = 1 , NSIGMA IF( ADISSC(I,1) .EQ. AREACT(JPL,1)) THEN C CHECK REMAINDER OF REACTION TO SEE IF A MATCH IS MADE. IT IS C NOT SUFFICIENT TO SEE IF A MATCH IS MADE FOR EVERY PRODUCT IN C THE CROSS-SECTION LISTING. FOR EXAMPLE, IN THE TWO REACTIONS: C O2 =2O , AND O2 = O + O(1D) , EVERY PRODUCT IN THE FIRST C REACTION CAN BE MATCHED UP WITH A PRODUCT IN THE SECOND C REACTION. THE EASIEST WAY TO AVOID THIS COMPLICATION , WITHOUT C COMPARING COEFFICIENTS, IS TO ALSO DO THE COMPARISON OF THE C PRODUCTS OF THE SECOND REACTION AGAINST THE FIRST. DO J = 2 , 4 IF( ( ADISSC(I,J) .NE. AREACT(JPL,2) ) .AND. & ( ADISSC(I,J) .NE. AREACT(JPL,3) ) .AND. & ( ADISSC(I,J) .NE. AREACT(JPL,4) ) )GOTO 30 ENDDO DO J = 2 , 4 IF( ( AREACT(JPL,J) .NE. ADISSC(I,2) ) .AND. & ( AREACT(JPL,J) .NE. ADISSC(I,3) ) .AND. & ( AREACT(JPL,J) .NE. ADISSC(I,4) ) )GOTO 30 ENDDO IPHOTO(JPL) = I C WRITE (6,*) IPHOTO(JPL), I GOTO 40 ENDIF 30 ENDDO WRITE (25,390) (ACOEFF(JPL,J) , AREACT(JPL,J) , SYMBOL(JPL,J) & , J = 1 , 4) 40 ENDDO C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 300 FORMAT(8(4X,F6.1)) 310 FORMAT() 320 FORMAT(8E10.3) 330 FORMAT(1X , A8 , 14X , I1 , 36X , I3 , 2X , I3) 340 FORMAT( 4(1X,A8,1X) , 20X , I3 , 2X , I3 ) 350 FORMAT( 8( F7.4 , 3X) ) 360 FORMAT(60X , I3 , 2X , I3) 370 FORMAT() 380 FORMAT('0' , A8, ' CONTRIBUTES TO THE OPTICAL DEPTH ', & 'BUT IS NOT INCLUDED IN THE OPACITY TABLES') 390 FORMAT('0' , 'CAN NOT FIND OPACITY VALUES FOR THE FOLLOWING', & ' REACTION:' , 4(A1,A8,A1)) C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ RETURN END C------------------------------------------------------------------------------ C------------------------------------------------------------------------------ C------------------------------------------------------------------------------ SUBROUTINE OUTPUT INCLUDE 'PARAMETERS.INC' CHARACTER*8 AMOL(MB) , AREACT(IB,6) , ATRANS(MB) , ATAU(MB) & , APHOTO(MB) , APCE(MB) CHARACTER SYMBOL(IB,6) , ACOEFF(IB,6) CHARACTER*2 ATOM(MB) CHARACTER*14 ABOUNDRY(5) COMMON & /BLOCK1/ NDISSC & /ONED12/ ILOWER(MB) , IUPPER(MB) , & XLOWER(MB) , XUPPER(MB) & /ONED17/ NTRANS & /ONED18/ NFORM(MB) , NDEST(MB) , & IFORM(IB,MB) , IDEST(IB,MB) & /BLOCK2/ A(IB) , B(IB) , C(IB) & /ONED40/ MTRANS(MB) & /ONED19/ KFORM(IB,MB) , KDEST(IB,MB) & /ONED27/ MOL(6,IB) , KOEFF(6,IB) & /ONED42/ NMOL & /BLOCK3/ NOREAC(IB) & /ONED43/ AMOL & /ONED54/ NPCE & /BLOCK4/ ATOM , ATRANS & /BLOCK5/ ACOEFF , SYMBOL , APHOTO , APCE , NPHOTO & /BLOCK15/ NTAU & /BLOCK34/ ATAU & /BLOCK37/ AREACT & /BLOCK38/ NREACT & /BLOCK39/ NX(10,MB) , HF(MB) , DELH(IB) , & NATOM , MAXFRM , MAXDES , MFORM , MDEST DATA ABOUNDRY/'NUMBER DENSITY' , 'VELOCITY ' , & 'MAX. VELOCITY ' , 'FLUX ' , & 'MIXING RATIO ' / C THIS SUBROUTINE PRINTS CORRELATION LISTS AND REACTION DATA SETS WRITE (26,400) NATOM , NMOL , NREACT , NTRANS , NPHOTO , NTAU WRITE (26,410) (ATOM(N) , N = 1 , NATOM) WRITE (26,420) DO M = 1 , NTRANS MVARY = MTRANS(M) WRITE(26,430) MVARY , AMOL(MVARY) , ABOUNDRY(ILOWER(MVARY)) , & XLOWER(MVARY) , ABOUNDRY(IUPPER(MVARY)) , & XUPPER(MVARY) ENDDO WRITE (26,440) AMOL(MFORM) , MAXFRM , AMOL(MDEST) , MAXDES C PRINT OUT INFORMATION ABOUT TRANSPORTED SPECIES IF ( NTRANS .NE. 0 ) THEN WRITE (26,450) DO N = 1 , NTRANS WRITE (26,460) ATRANS(N) ENDDO ENDIF C PRINT OUT INFORMATION ABOUT SPECIES WHICH ARE IN PHOTOCHEMICAL EQUILIBRIUM IF ( NPCE .NE. 0 ) THEN WRITE (26,470) DO N = 1 , NPCE WRITE (26,460) APCE(N) ENDDO ENDIF C PRINT OUT INFORMATION ABOUT PHOTODISSOCIATIVE SPECIES IF(NPHOTO.NE.0) THEN WRITE(26,480) DO I = 1 , NPHOTO WRITE(26,460) APHOTO(I) ENDDO ENDIF C PRINT OUT INFORMATION ABOUT SPECIES CONTRIBUTING TO ATMOSPHERIC OPACITY IF(NTAU .NE. 0 ) THEN WRITE(26,490) DO N = 1 , NTAU WRITE(26,460) ATAU(N) ENDDO ENDIF C ECHO INPUT REACTION LIST WRITE (26,500) DO I = 1 , NDISSC WRITE(26,510) I, (ACOEFF(I,J),AREACT(I,J),SYMBOL(I,J), J =1,6) ENDDO N1 = NDISSC + 1 WRITE (26,520) DO I = N1 , NREACT IF ( B(I) .NE. 0.0 ) THEN IF( C(I) .NE. 0.0) THEN WRITE (26,530) I , (ACOEFF(I,J),AREACT(I,J),SYMBOL(I,J), & J = 1,6), A(I) , B(I) ,C(I) ELSE WRITE (26,540) I ,(ACOEFF(I,J),AREACT(I,J),SYMBOL(I,J), & J = 1,6) , A(I) , B(I) ENDIF ELSE IF( C(I) .NE. 0.0) THEN WRITE (26,550) I , (ACOEFF(I,J),AREACT(I,J),SYMBOL(I,J), & J = 1,6), A(I) , C(I) ELSE WRITE (26,560) I ,(ACOEFF(I,J),AREACT(I,J),SYMBOL(I,J), & J = 1,6) , A(I) ENDIF ENDIF ENDDO C SUPPLY SPECIES SPECIFIC INFORMATION DO M = 1 , NMOL NF = NFORM(M) ND = NDEST(M) WRITE (26,570) M , AMOL(M) WRITE (26,580) NF IF(NF.NE.0) THEN DO J = 1 , NF I = IFORM(J,M) IF ( B(I).NE.0.0d0 ) THEN IF(C(I) .NE. 0.0d0) THEN WRITE (26,590) I, (ACOEFF(I,K),AREACT(I,K), & SYMBOL(I,K),K=1,6),A(I),B(I),C(I) ELSE WRITE (26,600) I , (ACOEFF(I,K),AREACT(I,K), & SYMBOL(I,K),K=1,6),A(I), B(I) ENDIF ELSE IF(C(I) .NE. 0.0d0) THEN WRITE (26,610) I, (ACOEFF(I,K),AREACT(I,K), & SYMBOL(I,K),K=1,6),A(I), C(I) ELSE WRITE (26,620) I , (ACOEFF(I,K),AREACT(I,K), & SYMBOL(I,K),K=1,6),A(I) ENDIF ENDIF ENDDO ENDIF WRITE(26, 630) ND IF( ND.NE.0) THEN DO J = 1 , ND I = IDEST(J,M) IF ( B(I).NE.0.0d0 ) THEN IF(C(I) .NE. 0.0d0) THEN WRITE (26,590) I, (ACOEFF(I,K),AREACT(I,K), & SYMBOL(I,K),K=1,6),A(I),B(I), C(I) ELSE WRITE (26,600) I , (ACOEFF(I,K),AREACT(I,K), & SYMBOL(I,K),K=1,6),A(I), B(I) ENDIF ELSE IF(C(I) .NE. 0.0d0) THEN WRITE (26,610) I, (ACOEFF(I,K),AREACT(I,K), & SYMBOL(I,K),K=1,6),A(I), C(I) ELSE WRITE (26,620) I , (ACOEFF(I,K),AREACT(I,K), & SYMBOL(I,K),K=1,6),A(I) ENDIF ENDIF ENDDO ENDIF ENDDO C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 400 FORMAT ( '1-D PHOTOCHEMISTRY MODEL INCLUDES ',I3,' ATOMS, ',I3, & ' SPECIES, AND ',I3,' REACTIONS'//1X,'NUMBER OF SPECIES', & ' WHICH ARE TRANSPORTED = ',I3/1X, & 'NUMBER OF PHOTODISSOCIATIVE SPECIES = ' , I3/1X, & 'NUMBER OF SPECIES CONTRIBUTING TO THE ATMOSPHERIC', & ' OPACITY = ' , I3) 410 FORMAT( /, 'THE ATOMS ARE:',3X,10(A2,1X)) 420 FORMAT(/, 'THE SPECIES ARE:',/) 430 FORMAT(' ', 1X,I3,'.',1X, A8, 1X , 'BOTTOM ', & A14, 1X , ' IS' , 1X , 1PE9.2 , ' TOP ' , & A14 , 1X , ' IS' 1X , 1PE9.2) 440 FORMAT(/,A8,' IS FORMED IN THE LARGEST NUMBER OF REACTIONS - ', & I3,/,1X, A8,' IS DESTROYED IN THE LARGEST NUMBER OF', & ' REACTIONS - ', I3) 450 FORMAT(/,'THE FOLLOWING SPECIES ARE TRANSPORTED: ' ) 460 FORMAT(' ' , 2X , A8) 470 FORMAT(/, 'THE SPECIES WHICH ARE IN PHOTOCHEMICAL ' , & 'EQUILIBRIUM ARE:' ) 480 FORMAT(/, 'THE PHOTODISSOCIATIVE SPECIES ARE:') 490 FORMAT(/, 'THE SPECIES WHICH MAKE A SIGNIFICANT CONTRIBUTION' & ,' TO THE ATMOSPHERIC OPACITY ARE: ') 500 FORMAT (/, 'PHOTODISSOCIATIVE REACTIONS:'///) 510 FORMAT (/,1X,I3,'.',1X,6(A1,A8,A1)) 520 FORMAT (/,'NON-PHOTODISSOCIATIVE REACTIONS:'///) 530 FORMAT (/,1X,I3,'.',1X,6(A1,A8,A1),/,5X,'K = ',1PE8.2, & '*T**',0PF4.1,'*EXP(',0PF6.0,'/T)',5X) 540 FORMAT (/,1X,I3,'.',1X,6(A1,A8,A1),/,5X,'K = ',1PE8.2, & '*T**',0PF4.1) 550 FORMAT (/,1X,I3,'.',1X,6(A1,A8,A1),/,5X,'K = ',1PE8.2, & '*EXP(',0PF6.0,'/T)') 560 FORMAT (/,1X,I3,'.',1X,6(A1,A8,A1),/,5X,'K = ',1PE8.2) 570 FORMAT (/,I3,'.',5X,A8) 580 FORMAT (/,5X,'FORMED IN ',I3,' REACTIONS :'/) 590 FORMAT (6X,'REACTION NO. ',I3,'.',5X,6(A1,A8,A1),/,6X, & 'K = ',1PE8.2,'*T**',0PF4.1,'*EXP(',0PF6.0,'/T)') 600 FORMAT (6X,'REACTION NO. ',I3,'.',5X,6(A1,A8,A1),/,6X, & 'K = ',1PE8.2,'*T**',0PF4.1,1X) 610 FORMAT (6X,'REACTION NO. ',I3,'.',5X,6(A1,A8,A1),/,6X, & 'K = ',1PE8.2,'*EXP(',0PF6.0,'/T)') 620 FORMAT (6X,'REACTION NO. ',I3,'.',5X,6(A1,A8,A1),/,6X, & 'K = ',1PE8.2) 630 FORMAT (/,5X,'DESTROYED IN ',I3, ' REACTIONS:'/) C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ RETURN END C------------------------------------------------------------------------------ C------------------------------------------------------------------------------ C------------------------------------------------------------------------------ SUBROUTINE EARTH(K) INCLUDE 'PARAMETERS.INC' INCLUDE 'PLANET.INC' COMMON & /SWITCHES/ ZLAT , DELTAT , DELT1 , DAYSTART , & TIMESTART , TIMESTOP , OUTTIME , & MODE , NTIME , NITER , IDIFF , istepdump & /BLOCK7/ RATE(IB,NALT) & /ONED28/ DEN(NALT) & /BLOCK8/ IPHOTO(IB) , LAST(IB) , & LOW(IB) , SIGMA(NWAVE,IB) & /BLOCK9/ WAVE(NWAVE) & /BLOCK11/ T(NALT) & /BLOCK12/ CROSS(NWAVE,MB) & /BLOCK13/ PRESS(NALT) & /BLOCK18/ CONC(MB,NALT) & /BLOCK19/ TRANS(NWAVE,NALT) & /BLOCK21/ CRSO2(11) & /BLOCK22/ O2SRB(6,16) & /BLOCK23/ TRANSRB(6,16,NALT) & /BLOCK24/ TRANLYM(11,NALT) & /BLOCK25/ H2OSRB(6,16) , SRBNO(2,6,16) , WNO(2,6,16) & /BLOCK26/ CRSH2O(11) , TO3(NWAVE) , TNO2(NWAVE), & TH2CO(NWAVE) , TCLONO2(NWAVE) , & ITHREE(IB) , IENDO(IB) , T1(IB) , T2(IB) , & T3(IB) , T4(IB) , F3(IB) , IEXO(IB) , & SIGN2O(NWAVE) , POWN2O(NWAVE) , T295(NWAVE) , & T250(NWAVE) , TCH3CL(NWAVE) , NTHREE , NENDOTH & /BLOCK27/ PI !!!mes 10/1/97 & /ONED35/ FLUX(NWAVE) & /BLOCK38/ NREACT & /ONED15/ THOUR & /ONED32/ SUNRISE(NALT) , SUNSET(NALT) & /ONED44/ SIGSRB(16,IB) & /ONED45/ SRBFLUX(16) DIMENSION SIGH2O(11) , RAINOUT(7) , CRSN2O(NWAVE) , & TROPNO2(7) , CRNO(26) , PERLYM(11) , SRBN2O(16) , & SCF2CL2(NWAVE) , SRBFREON(16) , SRB34(16) , & SCLONO2(NWAVE) , SRB40(16) , ADJ(NWAVE) , ADJSRB(16) , & N2ION(NALT) , O2ION(NALT) , AH2O2(8) , BH2O2(5) DIMENSION WEIGHT(6) DATA WEIGHT/0.05d0, 0.20d0, 0.25d0, 0.25d0, 0.20d0, 0.05d0 / DATA RAINOUT/ 2.31d-6 , 2.31d-6 , 2.31d-6 , 1.04d-6 , 4.66d-7 , & 2.10d-7 , 9.42d-8/ DATA PERLYM/ 6.6d-2 , 8.9d-2 , 11.1d-2 , 11.5d-2 , 8.9d-2 , & 5.6d-2 , 9.6d-2 , 11.5d-2 , 11.1d-2 , 8.9d-2 , & 6.3d-2 / DATA TROPNO2/ 4.0d+3 , 4.0d+3 , 4.0d+3 , 2.0d+3 , 9.0d+2 , & 4.0d+2 , 1.0d+2/ DATA CRNO/ 2.50d+0 , 3.50d+0 , 5.50d+0 , 9.50d+0 , 1.60d+1 , & 2.04d+1 , 2.21d+1 , 1.97d+1 , 1.48d+1 , 1.23d+1 , & 8.00d+0 , 4.90d+0 , 4.00d+0 , 2.50d+0 , 1.90d+0 , & 1.10d+0 , 8.50d-1 , 6.50d-1 , 4.50d-1 , 3.50d-1 , & 2.80d-1 , 2.20d-1 , 1.70d-1 , 1.30d-1 , 1.00d-1 , & 8.00d-2/ DATA N2ION/ 36*0.0d0 , 1.30d-19 , 1.00d-17 , 1.10d-15 , 6.00d-14 , & 4.90d-12 , 3.80d-10 , 3.00d-8 , 2.40d-6 , 1.80d-4 , & 1.40d-2 , 0.31d0 , 2.80d0 , 15.4d0 , 62.5d0 , & 183.0d0 , 363.0d0 , 545.0d0 , 632.0d0 , 635.0d0 , & 594.0d0 , 532.0d0 , 465.0d0 , 401.0d0 , 346.0d0 , & 302.0d0 , 270.0d0 , 245.0d0 , 230.0d0 , 220.0d0 , & 211.0d0 / DATA O2ION /36*0.0d0 , 3.80d-20 , 3.00d-18 , 3.40d-16 , 1.80d-14 , & 1.35d-12 , 1.07d-10 , 8.00d-9 , 6.40d-6 , 4.95d-5 , & 3.80d-3 , 0.08d0 , 1.60d0 , 5.80d0 , 40.2d0 , & 183.0d0 , 482.0d0 , 802.0d0 , 951.0d0 , 926.0d0 , & 821.0d0 , 702.0d0 , 588.0d0 , 485.0d0 , 397.0d0 , & 324.0d0 , 266.0d0 , 220.0d0 , 185.0d0 , 157.0d0 , & 134.0d0/ DATA AH2O2/ 6.4761d4 , -9.2170972d2 , 4.535649d0 , & -4.4589016d-3 , -4.035101d-5 , 1.6878206d-7 , & -2.652014d-10 , 1.5534675d-13 / DATA BH2O2/ 6.8123d3 , -5.1351d1 , 1.1522d-1 , & -3.0493d-5 , -1.0924d-7 / C SIGADJ MODIFIES CROSS-SECTIONS IN ORDER TO ACCOUNT FOR TEMPERATURE C DEPENDENCIES SIGADJ(S,TREF,TVAL,F) = S*( 1.0d0 + F*(TREF - TVAL) ) C RATE3 IS AN INTERPOLATION FUNCTION WHICH MODIFIES SELECTED 3-BODY C RATE CONSTANTS IN ORDER TO TAKE TEMPERATURE AND DENSITY EFFECTS INTO C ACCOUNT RATE3(Ra,Ta,Da,F3a) = ( Ra/( 1.0d0 + (Ra*Da/Ta) ) )* & F3a**( 1.0d0/( 1.0d0 + & (DLOG10(Ra*Da/Ta))*(DLOG10(Ra*Da/Ta)) ) ) C FUNCTIONS USED TO COMPUTE QUANTUM YIELD OF O(1D) NEAR 3000 ANGSTROMS A0(Ta) = 0.94932d0 - 1.7039d-4*Ta + 1.4072d-6*Ta*Ta A1(Ta) = -2.4052d-2 + 1.0479d-3*Ta - 1.0655d-5*Ta*Ta A2(Ta) = 1.8771d-2 - 3.6401d-4*Ta - 1.8587d-5*Ta*Ta A3(Ta) = -1.454d-2 - 4.7787d-5*Ta + 8.1277d-6*Ta*Ta A4(Ta) = 2.3287d-3 + 1.9891d-5*Ta - 1.1801d-6*Ta*Ta A5(Ta) = -1.4471d-4 - 1.7188d-6*Ta + 7.2661d-8*Ta*Ta A6(Ta) = 3.183d-6 + 4.6209d-8*Ta - 1.6266d-9*Ta*Ta W(Ta) = 308.2d0 + 4.4871d-2*Ta + 6.938d-5*Ta*Ta - 2.5452d-6*Ta*Ta*Ta TempValue = T(K) C----------------------------------------------------------------->>>>>> C PHOTOLYSIS REACTIONS - SKIP IF IN DARKNESS C-----------------------------------------------------------------<<<<<< IF( (MODE .NE. 0) .AND. ( (THOUR .GT. SUNSET(K)) .OR. & (THOUR .LT. SUNRISE(K)) ) ) GOTO 10 C******************************************* C SCHUMANN-RUNGE BAND CORRECTIONS TO O2 = 2O C******************************************* I = 1 IP = IPHOTO(I) ! REACTION NUMBER ACCORDING TO THE CROSS-SECTION LISTING R = RATE(I,K) C ADJUST SCHUMANN-RUNGE BANDS WITH IMPROVED O2 CROSS-SECTION DO LSRB = 1 , 16 SUM = 0.0D0 DO L = 1 , 6 ! DO SUBINTERVAL INTEGRATION SUM = SUM + WEIGHT(L)*(O2SRB(L,LSRB) - SIGSRB(LSRB,I))* & TRANSRB(L,LSRB,K) ENDDO R = R + SUM*SRBFLUX(LSRB) ENDDO RATE(I,K) = R C****************************** C CORRECTIONS TO O2 = O + O(1D) C****************************** I = 2 IP = IPHOTO(I) ! REACTION NUMBER ACCORDING TO THE CROSS-SECTION LISTING R = RATE(I,K) C REMOVE OLD LYMAN ALPHA PORTION R = R - SIGMA(LYMAN,IP)*FLUX(LYMAN)*TRANS(LYMAN,K) C ADD ON IMPROVED LYMAN ALPHA CONTRIBUTION DO LALPHA = 1 , 11 R = R + PERLYM(LALPHA)*FLUX(LYMAN)*TRANLYM(LALPHA,K)* & CRSO2(LALPHA) ENDDO RATE(I,K) = R C**************************************** C TEMPERATURE DEPENDENT OZONE CORRECTIONS C**************************************** DTEMP = 298.0d0 - TempValue b0 = A0(DTEMP) b1 = A1(DTEMP) b2 = A2(DTEMP) b3 = A3(DTEMP) b4 = A4(DTEMP) b5 = A5(DTEMP) b6 = A6(DTEMP) DO I = 3 , 4 IP = IPHOTO(I) ! REACTION NUMBER ACCORDING TO THE CROSS-SECTION LISTING IP3 = IPHOTO(3) IP4 = IPHOTO(4) R = RATE(I,K) DO LAMBDA = LYMAN + 30 , LYMAN + 45 C TOTAL O3 CROSS-SECTION BETWEEN 3025 AND 3175 ANGSTROMS SIGTOT = SIGMA(LAMBDA,IP3) + SIGMA(LAMBDA,IP4) CRO3 = SIGADJ(SIGTOT,273.0d0,TempValue,TO3(LAMBDA)) IF( (LAMBDA .GE. LYMAN + 30) .AND. & (LAMBDA .LE. LYMAN + 34) ) PHIO1D = 0.9d0 IF( (LAMBDA .GE. LYMAN + 35) .AND. & (LAMBDA .LE. LYMAN + 37) ) PHIO1D = 0.95d0 IF( (LAMBDA .GE. LYMAN + 41) ) PHIO1D = 0.0d0 IF( (LAMBDA .GE. LYMAN + 38) .AND. & (LAMBDA .LE. LYMAN + 40) ) THEN DL = 0.1d0*WAVE(LAMBDA) - 305.0d0 PHIO1D = b0 + b1*DL + b2*DL**2 + b3*DL**3 + & b4*DL**4 + b5*DL**5 + b6*DL**6 PHIO1D = DMIN1( PHIO1D , 0.95D0 ) PHIO1D = DMAX1( PHIO1D , 0.0D0 ) ENDIF IF( I .EQ. 4 )THEN P = PHIO1D ELSE P = 1.0d0 - PHIO1D ENDIF SIG = P*CRO3 R = R + FLUX(LAMBDA)*TRANS(LAMBDA,K)* & (SIG - SIGMA(LAMBDA,IP)) ENDDO RATE(I,K) = R ENDDO C**************************** C CORRECTIONS TO H2O = H + OH C**************************** I = 5 IP = IPHOTO(I) ! REACTION NUMBER ACCORDING TO THE CROSS-SECTION LISTING R = RATE(I,K) C REMOVE OLD LYMAN ALPHA PORTION R = R - SIGMA(LYMAN,IP)*FLUX(LYMAN)*TRANS(LYMAN,K) C ADD ON IMPROVED LYMAN ALPHA CONTRIBUTION DO LALPHA = 1 , 11 C CRSH2O IS THE TOTAL CROSS-SECTION OF H2O AT LYMAN ALPHA. WE NEED C THE CONTRIBUTIONS DUE TO EACH INDIVIDUAL BRANCH SO WE MUST MULTIPLY C CRSH2O BY 0.78 . TRANLYM IS THE TRANSMITANCE AT LYMAN ALPHA AND C IS COMPUTED IN SUBROUTINE DIULYM. SIGH2O(LALPHA) = 0.78d0*CRSH2O(LALPHA) R = R + (FLUX(LYMAN))*TRANLYM(LALPHA,K)* & PERLYM(LALPHA)*SIGH2O(LALPHA) ENDDO C ADJUST SRB CONTRIBUTION WITH IMPROVED H2O CROSS-SECTION DO LSRB = 1 , 16 SUM = 0.0D0 DO L = 1 , 6 ! DO SUBINTERVAL INTEGRATION SUM = SUM + WEIGHT(L)*(H2OSRB(L,LSRB) - SIGSRB(LSRB,I))* & TRANSRB(L,LSRB,K) ENDDO R = R + SUM*SRBFLUX(LSRB) ENDDO RATE(I,K) = R C************************************************* C CORRECTIONS TO H2O = H2 + O(1D) AND H2O = 2H + O C************************************************* DO I = 6 , 7 IP = IPHOTO(I) ! REACTION NUMBER ACCORDING TO THE CROSS-SECTION LISTING R = RATE(I,K) IF( I .EQ. 6) THEN FRACT = 0.10d0 ELSE FRACT = 0.12d0 ENDIF C REMOVE OLD LYMAN ALPHA PORTION R = R - SIGMA(LYMAN,IP)*FLUX(LYMAN)*TRANS(LYMAN,K) C ADD ON IMPROVED LYMAN ALPHA CONTRIBUTION DO LALPHA = 1 , 11 C CRSH2O IS THE TOTAL CROSS-SECTION OF H2O AT LYMAN ALPHA. WE NEED C THE CONTRIBUTIONS DUE TO EACH INDIVIDUAL BRANCH SO WE MUST MULTIPLY C CRSH2O BY FRACT. TRANLYM IS THE C TRANSMITANCE AT LYMAN ALPHA AND IS COMPUTED IN SUBROUTINE DIULYM. SIGH2O(LALPHA) = FRACT*CRSH2O(LALPHA) R = R + (FLUX(LYMAN))*TRANLYM(LALPHA,K)* & PERLYM(LALPHA)*SIGH2O(LALPHA) ENDDO RATE(I,K) = R ENDDO C************************************************************** C TEMPERATURE CORRECTIONS TO H2O2 = 2OH C************************************************************** I = 8 IP = IPHOTO(I) ! REACTION NUMBER ACCORDING TO THE CROSS-SECTION LISTING R = RATE(I,K) CHI = 1.0D0/( 1.0D0 + EXP(-1265D0/TempValue) ) DO LAMBDA = LYMAN+29 , LYMAN+46 IF( TempValue .GE. 200d0 .AND. TempValue .LE. 400d0 ) THEN WAVEnm = 0.1d0*WAVE(LAMBDA) ! wavelength in nanometers TERM1 = AH2O2(1) DO N = 2 , 8 TERM1 = TERM1 + AH2O2(N)*WAVEnm**(N-1) ENDDO TERM2 = BH2O2(1) DO N = 2 , 5 TERM2 = TERM2 + BH2O2(N)*WAVEnm**(N-1) ENDDO SIGH2O2 = 1.0D-21*(CHI*TERM1 + (1.0D0 - CHI)*TERM2) ELSE SIGH2O2 = SIGMA(LAMBDA,IP) ENDIF R = R + FLUX(LAMBDA)*TRANS(LAMBDA,K)* & (SIGH2O2 - SIGMA(LAMBDA,IP)) ENDDO RATE(I,K) = R C************************************************************** C NO PHOTODISSOCIATION IN SRB REGION. C************************************************************** I = 9 IP = IPHOTO(I) ! REACTION NUMBER ACCORDING TO THE CROSS-SECTION LISTING C THE CROSXX.INP FILE DOES NOT CONTAIN CROSS-SECTIONS FOR NO. R = 0.0D0 DO LSRB = 6 , 11 SUM = 0.0D0 DO L = 1 , 6 CNO = 0.0D0 DO N = 1 , 2 CNO = CNO + WNO(N,L,LSRB)*SRBNO(N,L,LSRB) ENDDO SUM = SUM + CNO*TRANSRB(L,LSRB,K) ENDDO R = R + SUM*SRBFLUX(LSRB) ENDDO RATE(I,K) = R RATE(I,K) = 0.0D0 C***************************************************** C TEMPERATURE DEPENDENT PHOTODISSOCIATION RATE FOR NO2 C***************************************************** I = 10 IP = IPHOTO(I) ! REACTION NUMBER ACCORDING TO THE CROSS-SECTION LISTING R = RATE(I,K) DO LAMBDA = LYMAN + 35 , LYMAN + 58 SIGNO2 = SIGADJ(SIGMA(LAMBDA,IP),298.0d0,TempValue,TNO2(LAMBDA)) R = R + FLUX(LAMBDA)*TRANS(LAMBDA,K)* & (SIGNO2 - SIGMA(LAMBDA,IP)) ENDDO RATE(I,K) = R C****************************************** C TEMPERATURE/WAVELENGTH CORRECTIONS TO N2O C****************************************** I = 13 IP = IPHOTO(I) ! REACTION NUMBER ACCORDING TO THE CROSS-SECTION LISTING TN2O = TempValue - 300.0d0 TN2O = DMAX1(TN2O , 194.0D0 - 300.0D0 ) TN2O = DMIN1(TN2O , 320.0D0 - 300.0D0 ) R = RATE(I,K) C N2O CROSS-SECTION IS A FUNCTION OF TEMPERATURE. THIS FACT MUST BE C INCLUDED IN THE SRB CALCULATION C CALCULATE TEMPERATURE DEPENDENT CROSS-SECTIONS DO LAMBDA = LYMAN + 11 , LYMAN + 24 CRSN2O(LAMBDA) = SIGN2O(LAMBDA) + & (TN2O)*EXP(POWN2O(LAMBDA)) CRSN2O(LAMBDA) = EXP(CRSN2O(LAMBDA)) ENDDO C INTERPOLATE ONTO SRB GRID DO LSRB = 1 , 16 SRBN2O(LSRB) = SIGINT(CRSN2O,LSRB) ENDDO C ADJUST SRB PORTION WITH IMPROVED CROSS-SECTION DO LSRB = 1 , 16 SUM = 0.0D0 DO L= 1 , 6 ! DO SUBINTERVAL INTEGRATION SUM = SUM + WEIGHT(L)*TRANSRB(L,LSRB,K) ENDDO R = R + SUM*SRBFLUX(LSRB)*(SRBN2O(LSRB) - SIGSRB(LSRB,I)) ENDDO C ADJUST REMAINING PORTION C THE LYMAN+11 INTERVAL IS A FRACTION INTERVAL, SINCE SRB BEGINS C HALFWAY THROUGH IT. R = R + ((1754.0d0-1725.0d0)/(1775.0d0-1725.0d0))* & FLUX(LYMAN+11)*TRANS(LYMAN+11,K)* & (CRSN2O(LYMAN+11) - SIGMA(LYMAN+11,IP)) DO LAMBDA = LYMAN + 17 , LYMAN + 24 R = R + FLUX(LAMBDA)*TRANS(LAMBDA,K)* & (CRSN2O(LAMBDA) - SIGMA(LAMBDA,IP)) ENDDO RATE(I,K) = R C******************************************* C TEMPERATURE/WAVELENGTH CORRECTIONS TO N2O5 C******************************************* I = 14 IP = IPHOTO(I) ! REACTION NUMBER ACCORDING TO THE CROSS-SECTION LISTING R = RATE(I,K) DO LAMBDA = LYMAN + 33 , LYMAN + 53 SN2O5 = (1.0d-20)*EXP( 2.735d0 + & (4728.5d0 - 1.7127d0*WAVE(LAMBDA))/TempValue ) R = R + FLUX(LAMBDA)*TRANS(LAMBDA,K)* & (SN2O5 - SIGMA(LAMBDA,IP)) ENDDO RATE(I,K) = R C****** C test of sensitivity to J_HCL C******* I = 25 C RATE(I,K) = 0.5*RATE(I,K) C******************************************** C TEMPERATURE/WAVELENGTH CORRECTIONS TO CH3CL C******************************************** I = 27 IP = IPHOTO(I) ! REACTION NUMBER ACCORDING TO THE CROSS-SECTION LISTING R = RATE(I,K) DO LAMBDA = LYMAN + 18 , LYMAN + 19 SCH3CL = SIGADJ(SIGMA(LAMBDA,IP),296.0d0, & TempValue,TCH3CL(LAMBDA)) R = R + FLUX(LAMBDA)*TRANS(LAMBDA,K)* & (SCH3CL - SIGMA(LAMBDA,IP)) ENDDO RATE(I,K) = R C******************************************************** C TEMPERATURE/WAVELENGTH CORRECTIONS TO FREON 12 (CF2CL2) C******************************************************** I = 33 IP = IPHOTO(I) ! REACTION NUMBER ACCORDING TO THE CROSS-SECTION LISTING C CALCULATE TEMPERATURE/WAVE DEPENDENT CROSS-SECTIONS DO LAMBDA = LYMAN + 10 , LYMAN + 24 POWF12 = (4.10d-4*(TempValue-298.0d0))* & (-184.9d0 + 0.1d0*WAVE(LAMBDA)) SCF2CL2(LAMBDA) = SIGMA(LAMBDA,IP)*EXP(POWF12) ENDDO C INTERPOLATE CROSS-SECTIONS ONTO SRB GRID DO LSRB = 1 , 16 SRBFREON(LSRB) = SIGINT(SCF2CL2,LSRB) ENDDO R = 0.0d0 LAMBDA = LYMAN + 10 POWF12 = (4.10d-4*(TempValue-298.0d0))* & (-184.9d0 + 0.1d0*WAVE(LAMBDA)) SCF2CL2(LAMBDA) = SIGMA(LAMBDA,IP)*EXP(POWF12) R = R + FLUX(LAMBDA)*TRANS(LAMBDA,K)*SCF2CL2(LAMBDA) LAMBDA = LYMAN + 11 C THE LYMAN+11 INTERVAL IS A FRACTION INTERVAL, SINCE SRB BEGINS C HALFWAY THROUGH IT. POWF12 = (4.10d-4*(TempValue-298.0d0))* & (-184.9d0 + 0.1d0*WAVE(LAMBDA)) SCF2CL2(LAMBDA) = SIGMA(LAMBDA,IP)*EXP(POWF12) R = R + ((1754.0d0-1725.0d0)/(1775.0d0-1725.0d0))* & FLUX(LAMBDA)*TRANS(LAMBDA,K)*SCF2CL2(LAMBDA) C ADJUST SRB REGION WITH IMPROVED CROSS-SECTIONS DO LSRB = 1 , 16 SUM = 0.0D0 DO L = 1 , 6 ! DO SUBINTERVAL INTEGRATION SUM = SUM + WEIGHT(L)*TRANSRB(L,LSRB,K) ENDDO R = R + SUM*SRBFREON(LSRB)*SRBFLUX(LSRB) ENDDO C ADD ON REMAINING PORTION DO LAMBDA = LYMAN + 17 , LYMAN + 24 R = R + FLUX(LAMBDA)*TRANS(LAMBDA,K)*SCF2CL2(LAMBDA) ENDDO RATE(I,K) = R C********************************************** C TEMPERATURE/WAVELENGTH CORRECTIONS TO C2H3CL3 C********************************************** I = 34 IP = IPHOTO(I) ! REACTION NUMBER ACCORDING TO THE CROSS-SECTION LISTING TEMP = DMAX1(TempValue , 250.0D0) R = RATE(I,K) C THE TEMPERATURE ADJUSTMENTS COVER THE REGION [LYMAN+16,LYMAN+24] C SINCE LYMAN+16 IS WITHIN THE SCHUMANN-RUNGE BANDS WE MUST FIRST C CLACULATE THE TEMPERATURE DEPENDENT CROSS-SECTION AT LYMAN+16 C AND THEN INTERPOLATE IT ONTO THE SRB GRID C COMPUTE TEMPERATURE DEPENDENT CROSS-SECTION AT LYMAN+16 SC2H3CL3 = SIGADJ(SIGMA(LYMAN+16,IP), & 295.0d0,TEMP,T295(LYMAN+16)) IF( TempValue .LT. 250.0d0 ) THEN SC2H3CL3 = SIGADJ(SC2H3CL3,250.0d0,TempValue,T250(LYMAN+16)) ENDIF C INTERPOLATE ONTO SRB GRID SRB34(14) = ((1975.0d0-1972.0d0)/(1985.0d0-1972.0d0))* & SIGSRB(14,I) + & ((1985.0d0-1975.0d0)/(1985.0d0-1972.0d0))* & SC2H3CL3 SRB34(15) = SC2H3CL3 SRB34(16) = SC2H3CL3 C ADJUST SRB PORTION DO LSRB = 14 , 16 ! LYMAN+16 REGION SUM = 0.0D0 DO L = 1 , 6 ! DO SUBINTERVAL INTEGRATION SUM = SUM + WEIGHT(L)*TRANSRB(L,LSRB,K) ENDDO R = R + SUM*(SRB34(LSRB) - SIGSRB(LSRB,I))*SRBFLUX(LSRB) ENDDO C CORRECT CROSS-SECTIONS BEYOND SRB DO LAMBDA = LYMAN + 17 , LYMAN + 24 SC2H3CL3 = SIGADJ(SIGMA(LAMBDA,IP),295.0d0,TEMP,T295(LAMBDA)) IF( TempValue .LT. 250.0d0 ) THEN SC2H3CL3 = SIGADJ(SC2H3CL3,250.0d0,TempValue,T250(LAMBDA)) ENDIF R = R + FLUX(LAMBDA)*TRANS(LAMBDA,K)* & (SC2H3CL3 - SIGMA(LAMBDA,IP)) ENDDO RATE(I,K) = R C************************************************************ C TEMPERATURE CORRECTIONS TO PHOTODISSOCIATION RATE OF CLONO2 C************************************************************ I = 40 IP = IPHOTO(I) ! REACTION NUMBER ACCORDING TO THE CROSS-SECTION LISTING R = RATE(I,K) C THE TEMPERATURE ADJUSTMENTS COVER THE REGION [LYMAN+14,LYMAN+56] C SINCE [LYMAN+14,LYMAN+16] IS WITHIN THE SCHUMANN-RUNGE BANDS WE MUST FIRST C CLACULATE THE TEMPERATURE DEPENDENT CROSS-SECTION AT [LYMAN+14,LYMAN+16] C AND THEN INTERPOLATE IT ONTO THE SRB GRID C COMPUTE TEMPERATURE ADJUSTED CROSS-SECTIONS SCLONO2(LYMAN+13) = SIGMA(LYMAN+13,IP) ! NEEDED FOR INTERPOLATION DO LAMBDA = LYMAN + 14 , LYMAN + 56 SCLONO2(LAMBDA) = SIGADJ(SIGMA(LAMBDA,IP),296.0d0, & TempValue,TCLONO2(LAMBDA)) ENDDO DO LSRB = 9 , 16 ! [LYMAN+14,LYMAN+16] SRB40(LSRB) = SIGINT(SCLONO2,LSRB) ENDDO C ADJUST SRB PORTION DO LSRB = 9 , 16 SUM = 0.0D0 DO L = 1 , 6 ! DO SUBINTERVAL INTEGRATION SUM = SUM + WEIGHT(L)*TRANSRB(L,LSRB,K) ENDDO R = R + SUM*(SRB40(LSRB) - SIGSRB(LSRB,I))*SRBFLUX(LSRB) ENDDO C CORRECT CROSS-SECTIONS BEYOND SRB DO LAMBDA = LYMAN + 17 , LYMAN + 56 R = R + FLUX(LAMBDA)*TRANS(LAMBDA,K)* & (SCLONO2(LAMBDA) - SIGMA(LAMBDA,IP)) ENDDO RATE(I,K) = R C**************************************************** C TEMPERATURE/WAVELENGTH DEPENDENT CORRECTIONS TO CO2 C**************************************************** I = 48 IP = IPHOTO(I) ! REACTION NUMBER ACCORDING TO THE CROSS-SECTION LISTING C CO2 PHOTOLYSIZES IN THE WAVELENGTH RANGE [LYMAN-38,LYMAN+17] C SINCE CHEM1D ONLY CALCULATES FROM SRB ON UP, THE ONLY REGION OF C INTEREST TO US IS [LYMAN,LYMAN+17]. THE CROSS-SECTION HAS A WAVELENGTH C DEPENDENT ADJUSTMENT FACTOR. C COMPUTE ADJUSTMENT FACTOR DO LAMBDA = LYMAN , LYMAN + 17 IF(WAVE(LAMBDA) .GT. 1970.0d0) THEN P = 0.9874d0 - 0.0015d0*((WAVE(LAMBDA) - 1950.0d0)/50.0d0) ENDIF IF( (WAVE(LAMBDA) .GE. 1450.0d0) .AND. & (WAVE(LAMBDA) .LE. 1970.0d0) )THEN DW = WAVE(LAMBDA) - 1450.0d0 P = 6.2548608d-2 - 5.0750811d-4*DW + & 1.1945104d-5*DW**2 - 1.2276127d-8*DW**3 P = 1.0d0 - DMAX1(P,0.0D0)/100.0d0 ENDIF IF( (WAVE(LAMBDA) .GE. 1200.0d0) .AND. & (WAVE(LAMBDA) .LT. 1450.0d0) )THEN DW = WAVE(LAMBDA) - 1200.0d0 P = 0.1590215d0 - 6.1266460d-3*DW + & 1.4502171d-4*DW**2 - 1.3434882d-6*DW**3 + & 5.1814042d-9*DW**4 - 7.0675488d-12*DW**5 P = 1.0d0 - DMAX1(P,0.0D0)/100.0d0 ENDIF ADJ(LAMBDA) = P**(298.0d0 - TempValue) ENDDO R = 0.0d0 C COMPUTE LYMAN ALPHA CONTRIBUTION C ADD ON LYMAN ALPHA CONTRIBUTION DO LALPHA = 1 , 11 R = R + PERLYM(LALPHA)*FLUX(LYMAN)*TRANLYM(LALPHA,K)* & ADJ(LYMAN)*SIGMA(LYMAN,IP) ENDDO C ADD ON SCHUMAN-RUNGE CONTINUUM DO LAMBDA = LYMAN + 1 , LYMAN + 10 R = R + FLUX(LAMBDA)* & ADJ(LAMBDA)*SIGMA(LAMBDA,IP)*TRANS(LAMBDA,K) ENDDO LAMBDA = LYMAN + 11 C THE LYMAN+11 INTERVAL IS A FRACTION INTERVAL, SINCE SRB BEGINS C HALFWAY THROUGH IT. R = R + ((1754.0d0-1725.0d0)/(1775.0d0-1725.0d0))*FLUX(LAMBDA)* & ADJ(LAMBDA)*SIGMA(LAMBDA,IP)*TRANS(LAMBDA,K) C COMPUTE SRB PORTION C INTERPOLATE ADJUSTMENT FACTOR ONTO SRB GRID DO LSRB = 1 , 16 ADJSRB(LSRB) = SIGINT(ADJ,LSRB) ENDDO DO LSRB = 1 , 16 SUM = 0.0D0 DO L = 1 , 6 ! DO SUBINTERVAL INTEGRATION SUM = SUM + WEIGHT(L)*TRANSRB(L,LSRB,K) ENDDO R = R + SUM*ADJSRB(LSRB)*SIGSRB(LSRB,I)* & SRBFLUX(LSRB) ENDDO C ADD ON LYMAN+17 PORTION LAMBDA = LYMAN + 17 R = R + FLUX(LAMBDA)* & ADJ(LAMBDA)*SIGMA(LAMBDA,IP)*TRANS(LAMBDA,K) RATE(I,K) = R C****************************************************** C TEMPERATURE DEPENDENT PHOTODISSOCIATION RATE FOR H2CO C****************************************************** C THERE ARE TWO SUCH REACTIONS: #49 H2CO = HCO + H C AND #50 H2CO = H2 + CO. THERE IS AN ADDITIONAL PRESSURE DEPENDENCE FOR C REACTION NUMBER 50. I = 49 IP = IPHOTO(I) ! REACTION NUMBER ACCORDING TO THE CROSS-SECTION LISTING R = RATE(I,K) DO LAMBDA = LYMAN + 24 , LYMAN + 48 SIGH2CO = SIGADJ(SIGMA(LAMBDA,IP),290.0d0, & TempValue,TH2CO(LAMBDA)) R = R + FLUX(LAMBDA)*TRANS(LAMBDA,K)* & (SIGH2CO - SIGMA(LAMBDA,IP)) ENDDO RATE(I,K) = R I = 50 IP = IPHOTO(I) ! REACTION NUMBER ACCORDING TO THE CROSS-SECTION LISTING R = RATE(I,K) DO LAMBDA = LYMAN + 24 , LYMAN + 48 SIGH2CO = SIGADJ(SIGMA(LAMBDA,IP),290.0d0, & TempValue,TH2CO(LAMBDA)) IF( LAMBDA .GE. (LYMAN + 43) ) THEN RATIO = (WAVE(LAMBDA) - 3290.0d0)/(3640.0d0 - WAVE(LAMBDA)) ADJUST = (1.0d0 + RATIO)/(1.0d0 + (PRESS(K)/1013.0d0)*RATIO ) SIGH2CO = SIGH2CO*ADJUST ENDIF R = R + FLUX(LAMBDA)*TRANS(LAMBDA,K)* & (SIGH2CO - SIGMA(LAMBDA,IP)) ENDDO RATE(I,K) = R 10 CONTINUE C************************************ C ADJUSTMENTS TO THREE-BODY REACTIONS C************************************ DO I3 = 1 , NTHREE I = ITHREE(I3) RATE(I,K) = RATE3( RATE(I,K) , T1(I3)/TempValue**T2(I3) , & DEN(K) , F3(I3) ) ENDDO DO IE = 1 , NENDOTH I = IENDO(IE) RATE(I,K) = RATE(IEXO(IE),K)/( T3(IE)*EXP(T4(IE)/TempValue) ) ENDDO C**************************** C COSMIC RAY PRODUCTION OF NO C**************************** I = 117 IF( K .LE. 26 ) RATE(I,K) = CRNO(K) C*************************** C TROPOSPHERIC SOURCE OF NO2 C*************************** I = 118 IF( K .LE. 7 ) RATE(I,K) = TROPNO2(K) C******************************** C CORRECTION TO CO + OH = CO2 + H C******************************** I = 227 c print*,'in earth',i,k,rate(i,k) c stop RATE(I,K) = RATE(I,K)*(1.0d0 + 0.6d0*PRESS(K)/1013.0d0) C******************************************************************* C RAINOUT RATE CONSTANTS FOR H2O2, HNO2, HNO3, HO2NO2, HF, HBR, HCL, C H2CO, AND CH3OOH C******************************************************************* DO I = 242,254 IF (K .LE. 7) RATE(I,K) = RAINOUT(K) ENDDO C*********************************** C SURFACE CATALYSIS ON METEORIC DUST C*********************************** C --- h2 + o --dust--> h2o c I = 250 c -------------------- mes 10/2/97 --------------------------------- c -- rate for surface catalysis of O and H2 on meteoric dust -- c idmode = 0 c kdustbot = 30 c kdusttop = 35 c zmassh = 1.67d-24 !mass of hydrogen atom c zmasso = 16.0d0*zmassh c zmassh2 = 2.0d0*zmassh c rcc = 1.0d0 !reaction coefficient c surfa = 0.0d0 c -- constant dust surface area c if (idmode.eq.0) then c surfa = RATE(I,K) !actually it's the surface area/unit volume c endif c -- gaussian distribution of dust surface area c if (idmode.eq.1) then c zkdp = 33.0d0 ! dust peak grid point c zkdw = 3.0d0 ! dust "width" c zzk = float(K) c surfa = RATE(I,K)*exp(-abs(((zzk-zkdp)/zkdw))**2.0d0) c endif c -- calculate surface reaction rate c if(K.lt.kdustbot.or.K.gt.kdusttop) rate(i,K) = 0.0d0 c if((K.ge.kdustbot).and.(K.le.kdusttop)) then c cbaro = sqrt(BOLTZ*T(K)/(2.0d0*PI*zmasso)) !O mean speed c cbarh2 = sqrt(BOLTZ*T(K)/(2.0d0*PI*zmassh2)) !H2 mean speed c rateohits = conc(1,K)*cbaro !*surface area c rateh2hits = conc(16,K)*cbarh2 !*surface area c if(rateohits.gt.rateh2hits) rate(I,K)=cbarh2*rcc*surfa/conc(1,K) c if(rateh2hits.gt.rateohits) rate(I,K)=cbaro*rcc*surfa/conc(16,K) c endif c RATE(I,K) = 0.0d0 c c -------------------------------------------------------------------- C************************** C E-REGION IONIZATION RATES C************************** RATE(259,K) = 0.0 !N2 IONIZATION (IN2 + PEIN2) RATE(260,K) = 0.0 !O2 IONIZATION (IO2 + PEIO2) RATE(261,K) = 0.0 !N2 DISSOCIATION W/ BRANCHING RATE(262,K) = 0.0 C RATE(259,K) = N2ION(K) !N2 IONIZATION (IN2 + PEIN2) C RATE(260,K) = O2ION(K) !O2 IONIZATION (IO2 + PEIO2) C RATE(261,K) = N2ION(K)*0.3d0 !N2 DISSOCIATION W/ BRANCHING C RATE(262,K) = N2ION(K)*0.7d0 RETURN END C------------------------------------------------------------------------------ C------------------------------------------------------------------------------ C------------------------------------------------------------------------------ SUBROUTINE ADJUST INCLUDE 'PARAMETERS.INC' DIMENSION O2LYM(5) , H2OLYM(5) DIMENSION FNO2(5,7) , FCLONO2(5,5) , FH2CO(5,3) , FO3(5,2) , & F295(5) , F250(5) , FCH3CL(5) DIMENSION LNO2(8) , LCLONO2(6) , LH2CO(4) , LO3(3) DIMENSION AN2O(5) , BN2O(4) , ODF(6) DIMENSION SIGTEMP(NWAVE) COMMON & /BLOCK1/ NDISSC & /BLOCK6/ MCROSS(MB) & /BLOCK8/ IPHOTO(IB) , LAST(IB) , & LOW(IB) , SIGMA(NWAVE,IB) & /BLOCK9/ WAVE(NWAVE) & /BLOCK12/ CROSS(NWAVE,MB) & /BLOCK21/ CRSO2(11) & /BLOCK22/ O2SRB(6,16) & /BLOCK25/ H2OSRB(6,16) , SRBNO(2,6,16) , WNO(2,6,16) & /BLOCK26/ CRSH2O(11) , TO3(NWAVE) , TNO2(NWAVE) , & TH2CO(NWAVE) , TCLONO2(NWAVE) , & ITHREE(IB) , IENDO(IB) , T1(IB) , T2(IB) , & T3(IB) , T4(IB) , F3(IB) , IEXO(IB) , & SIGN2O(NWAVE) , POWN2O(NWAVE) , T295(NWAVE) , & T250(NWAVE) , TCH3CL(NWAVE) , NTHREE , NENDOTH & /BLOCK40/ O3SRB(16) & /ONED44/ SIGSRB(16,IB) & /ONED45/ SRBFLUX(16) & /ONED48/ RAYSRB(16) & /ONED51/ ALBEDO(NWAVE) DATA O2LYM/8.9470039d-21 , -1.1441759d-20 , 1.4286004d-20 , & -2.1569135d-21 , -2.0319808d-21/ DATA H2OLYM/1.5037800d-17 , 2.8891091d-18 , 1.1857752d-18 , & -8.2043638d-19 , -1.2199025d-18/ DATA FNO2/ 6.617924d+00 , -5.331151d-03 , 8.912312d-07 , & 2.000057d-10 , -4.995888d-14 , & -1.147617d+01 , 8.789331d-03 , -1.373991d-06 , & -3.143766d-10 , 7.362594d-14 , & 2.104963d+01 , -1.481654d-02 , 2.176160d-06 , & 4.233010d-10 , -9.331119d-14 , & -8.030047d+01 , 5.382994d-02 , -7.524939d-06 , & -1.396228d-09 , 2.927822d-13 , & 3.415726d+01 , -2.181919d-02 , 2.929183d-06 , & 5.003484d-10 , -1.006985d-13 , & 1.323501d+00 , -6.782789d-04 , 8.683617d-08 , & 0.000000d+00 , 0.000000d+00 , & 3.983415d-02 , -1.003122d-05 , 0.000000d+00 , & 0.000000d+00 , 0.000000d+00/ DATA FCLONO2/ 7.656614d-01 , -9.061447d-04 , 2.219188d-07 , & 7.314442d-11 , -2.667751d-14 , & 5.892475d-03 , 4.179946d-06 , -1.580935d-09 , & -1.854288d-12 , 5.348531d-16 , & 1.593791d+00 , -1.215318d-03 , 1.936901d-07 , & 4.022887d-11 , -9.637668d-15 , & 5.892661d-02 , -5.395108d-05 , 8.388010d-09 , & 2.044443d-12 , -4.244845d-16 , & 2.215180d+00 , -1.111861d-03 , 1.394046d-07 , & 0.000000d+00 , 0.000000d+00/ DATA FH2CO/ 1.251090d+02 , -1.183672d-01 , 2.387435d-05 , & 5.684730d-09 , -1.721222d-12 , & -2.432510d-01 , 1.918508d-04 , -3.297522d-08 , & -5.543820d-12 , 1.407500d-15 , & 1.469056d+01 , -1.037230d-02 , 1.502076d-06 , & 3.122914d-10 , -6.782839d-14/ DATA FO3/ 4.378123d-02 , -3.859602d-04 , 5.299228d-07 , & 3.326920d-09 , -8.250185d-12 , & 3.836849d-01 , -1.152642d-03 , 0.000000d+00 , & 0.000000d+00 , 0.000000d+00/ DATA F295/ 3.519039d-01 , -3.842405d-03 , 8.032397d-06 , & 3.686217d-08 , -1.247817d-10/ DATA F250/ 4.868972d-01 , -5.398480d-03 , 1.152140d-05 , & 5.104976d-08 , -1.728080d-10/ DATA FCH3CL/ 7.123261d+01 , -8.130570d-02 , 1.923681d-05 , & 6.236362d-09 , -2.215346d-12/ DATA LNO2/ 35 , 38 , 42 , 46 , 49 , 52 , 56 , 58 / DATA LCLONO2/ 14 , 20 , 35 , 42 , 54 , 56 / DATA LH2CO/ 24 , 27 , 43 , 48 / DATA LO3/ 30 , 43 , 45 / DATA AN2O/ 68.21023d0 , -4.071805d0 , 4.301146d-2 , -1.777846d-4 , & 2.520672d-7/ DATA BN2O/ 123.4014d0 , -2.116255d0 , 1.111572d-2 , -1.881058d-5/ C THIS SUBROUTINE SETS UP PARAMETERS USED FOR MAKING ADJUSTMENTS TO THE C REACTION RATES. THESE MODIFICATIONS ARE PRIMARILY DUE TO C TEMPERATURE/PRESSURE EFFECTS, LYMAN ALPHA AND SCHUMANN RUNGE BAND C DESCRIPTIONS. MOST OF THE PARAMETERS COMPUTED IN THIS PORTION ARE FEED C DIRECTLY TO SUBROUTINE EARTH. THESE PARAMETERS ARE COMPUTED IN THIS C SECTION IN ORDER TO AVOID HAVING TO RECOMPUTE THEM ON EACH NEW CALL C TO SUBROUTINE EARTH. C SET UP THE GROUND ALBEDO FOR DIFFUSE FLUX ADJUSTMENTS DO LAMBDA = 1 , NWAVE IF( WAVE(LAMBDA) .LE. 3000.0d0 ) ALBEDO(LAMBDA) = 0.0D0 IF( (WAVE(LAMBDA) .GT. 3000.0d0 ) .AND. & (WAVE(LAMBDA) .LE. 4000.0d0 ) ) THEN ALBEDO(LAMBDA) = (3.0D-4)*WAVE(LAMBDA) - 0.9D0 ENDIF IF( WAVE(LAMBDA) .GT. 4000.0d0 ) ALBEDO(LAMBDA) = 0.3D0 ENDDO C************************************************* C READ IN CORRECTIONS FOR THE THREE-BODY REACTIONS C************************************************* DO I = 1 , IB ITHREE(I) = 0 IENDO(I) = 0 ENDDO I = -1 I3 = 0 DO WHILE( I .NE. 0 ) READ(29,130) I IF ( I .NE. 0 ) THEN I3 = I3 + 1 ITHREE(I3) = I READ(29,140) T1(I3) , T2(I3) , F3(I3) ENDIF ENDDO NTHREE = I3 C READ IN ENDOTHERMIC 3-BODY REACTIONS. THESE REACTION RATES ARE EXPRESSED IN C TERMS OF THE CORRESPONDING FORWARD REACTION. E.G. THE RATE FOR C A + B + M = C + D + M (ENDOTHERMIC) IS EXPRESSED IN TERMS OF THE RATE FOR C C + D + M = A + B + M (EXOTHERMIC). C REACTIONS TO BE MODIFIED ARE REPRESENTED BY IENDO(I) AND THE CORRESPONDING C RELATED REACTION IS GIVEN BY IEXO(I) I = -1 IE = 0 DO WHILE(I .NE. 0 ) READ(29,130) I IF( I .NE. 0 ) THEN IE = IE + 1 IENDO(IE) = I READ(29,130) IEXO(IE) READ(29,150) T3(IE) , T4(IE) ENDIF ENDDO NENDOTH = IE C**************************** C SCHUMANN-RUNGE BAND EFFECTS C**************************** C THERE ARE 6 COURSE INTERVALS COVERING THE SCHUMANN-RUNGE BANDS C (LAMBDA = 90 TO LAMBDA = 95). THIS REGION IS DIVIDED UP INTO 16 FINER C INTERVALS AND THE O2 CROSS-SECTION IS COPUTED WITHIN EACH OF THESE INTERVALS C USING THE MINSCHWANER OPACITY DISTRIBUTION FUNCTION. C THIS PORTION READS IN THE CROSS-SECTIONS AND QUADRATURE RULE USED IN C CALCULATING THE OPACITY DISTRIBUTION FUNCTION READ(30,100) DO LSRB = 1 , 16 READ(30,110) (ODF(L) , L = 1 , 6) DO L = 1 , 6 O2SRB(L,LSRB) = ODF(L) ENDDO ENDDO READ(30,100) DO LSRB = 1 , 16 READ(30,110) (ODF(L) , L = 1 , 6) DO L = 1 , 6 H2OSRB(L,LSRB) = ODF(L) ENDDO ENDDO READ(30,100) DO LSRB = 1 , 16 READ(30,120) WAVE1 , WAVE2 , SUSIM SRBFLUX(LSRB) = (WAVE2-WAVE1)*SUSIM ENDDO C READ IN NO CROSS-SECTIONS AND QUADRATURE WEIGHTS DO LSRB = 1 , 16 DO L = 1 , 6 DO I = 1 ,2 SRBNO(I,L,LSRB) = 0.0D0 WNO(I,L,LSRB) = 0.0D0 ENDDO ENDDO ENDDO LSRB = 6 DO I = 1 , 2 READ(32,160)(SRBNO(I,L,LSRB) , L = 1 , 6) READ(32,170)(WNO(I,L,LSRB) , L = 1 , 6) ENDDO LSRB = 7 DO I = 1 , 2 READ(32,160)(SRBNO(I,L,LSRB) , L = 1 , 6) READ(32,170)(WNO(I,L,LSRB) , L = 1 , 6) ENDDO LSRB = 11 DO I = 1 , 2 READ(32,160)(SRBNO(I,L,LSRB) , L = 1 , 6) READ(32,170)(WNO(I,L,LSRB) , L = 1 , 6) ENDDO C INTERPOLATE O3 CROSS-SECTIONS ONTO MINSCHWANER GRID DO LAMBDA = 1 , NWAVE SIGTEMP(LAMBDA) = CROSS(LAMBDA,MCROSS(2)) ENDDO DO LSRB = 1 , 16 O3SRB(LSRB) = SIGINT(SIGTEMP,LSRB) ENDDO C INTERPOLATE RAYLEIGH SCATTERING CROSS-SECTIONS ONTO SRB GRID DO LAMBDA = 1 , NWAVE SIGTEMP(LAMBDA) = CROSS(LAMBDA,MCROSS(3)) ENDDO DO LSRB = 1 , 16 RAYSRB(LSRB) = SIGINT(SIGTEMP,LSRB) ENDDO C INTERPOLATE CROSS-SECTIONS ONTO SRB GRID DO I = 1 , NDISSC DO LSRB = 1 , 16 SIGSRB(LSRB,I) = 0.0d0 ENDDO ENDDO DO I = 1 , NDISSC DO LAMBDA = 1 , NWAVE IF(IPHOTO(I) .GT. 0 ) THEN SIGTEMP(LAMBDA) = SIGMA(LAMBDA,IPHOTO(I)) ELSE SIGTEMP(LAMBDA) = 0.0d0 ENDIF ENDDO DO LSRB = 1 , 16 SIGSRB(LSRB,I) = SIGINT(SIGTEMP,LSRB) ENDDO ENDDO C********************************************************************** C ADJUST HERZBERG CONTINUUM CROSS-SECTIONS C********************************************************************** C THIS FORMULA COMES FROM: Planet Space Science, 36, 1469-1475, 1988 S0 = 6.884D-24 W0 = 2048.7D0 P0 = -69.7374D0 DO LAMBDA = 1 , NWAVE IF( WAVE(LAMBDA) .GE. W0 ) THEN CROSS(LAMBDA,MCROSS(1)) = S0*(W0/WAVE(LAMBDA))*DEXP( & P0*(DLOG(W0/WAVE(LAMBDA)))**2 ) SIGMA(LAMBDA,IPHOTO(1)) = CROSS(LAMBDA,MCROSS(1)) ENDIF ENDDO C********************************************************************** C ADJUST RAYLEIGH SCATTERING CROSS-SECTIONS C********************************************************************** C THIS FORMULA COMES FROM: Planet. Space Science, 30, 935-983, 1982 DO LAMBDA = 1 , NWAVE W = WAVE(LAMBDA)*1.0d-4 ! CONVERT TO MICRONS X = 0.074d0*W + 0.05d0/W - 0.084d0 CROSS(LAMBDA,MCROSS(3)) = 3.93d-28/W**(4.0d0+X) ENDDO C********************************************************************** C THIS SECTION SETS UP THE QUANTITIES NEEDED TO INCORPORATE TEMPERATURE C EFFECTS INTO THE CROSS-SECTION OF OZONE ( I = 3 AND 4 ) C********************************************************************** DO N = 1 , 2 L1 = LYMAN + LO3(N) L2 = LYMAN + LO3(N+1) DO LAMBDA = L1 , L2 TO3(LAMBDA) = FO3(1,N) DO J = 2 , 5 TO3(LAMBDA) = TO3(LAMBDA) + & FO3(J,N)*(WAVE(LAMBDA)/10.0d0)**(J-1) ENDDO TO3(LAMBDA) = TO3(LAMBDA)/(1.0d0 - 70.0d0*TO3(LAMBDA)) ENDDO ENDDO C********************************************************************** C CALCULATE O2 AND H2O CROSS-SECTIONS NEAR LYMAN ALPHA ( I = 5 , 6 , 7) C REF: LEWIS ET. AL. JGR 88 , 4935 (1983) C********************************************************************** C THIS PORTION PRODUCES A FOURTH-ORDER POLYNOMIAL FIT TO THE H2O C CROSS-SECTION ABOUT THE LYMAN ALPHA LINE (1215.67 A). THE C CROSS-SECTION IS SAMPLED WITH A 0.1 ANGSTROM RESOLUTION IN THE C INTERVAL [1215.2 , 1216.2]. C PERLYM(L) IS THE PORTION OF THE SOLAR FLUX IN EACH OF THE C 0.1 ANGSTROM INTERVALS. CONSEQUENTLY, PERLYM(L)*FLUX(LYMAN) WILL C SPECIFY THE TOTAL SOLAR FLUX CORRESPONDING TO EACH OF THE 0.1 ANGSTROM C INTERVALS. DO L = 1 , 11 CRSO2(L) = O2LYM(1) CRSH2O(L) = H2OLYM(1) DL = 0.1d0*(L-1) - 0.5d0 DO J = 2 , 5 CRSO2(L) = CRSO2(L) + O2LYM(J)*DL**(J-1) CRSH2O(L) = CRSH2O(L) + H2OLYM(J)*DL**(J-1) ENDDO ENDDO C**************************************************** C TEMPERATURE/WAVELENGTH CORRECTIONS TO NO2 ( I = 10) C**************************************************** DO N = 1 , 7 L1 = LYMAN + LNO2(N) L2 = LYMAN + LNO2(N+1) DO LAMBDA = L1 , L2 TNO2(LAMBDA) = FNO2(1,N) DO J = 2 , 5 TNO2(LAMBDA) = TNO2(LAMBDA) + & FNO2(J,N)*WAVE(LAMBDA)**(J-1) ENDDO ENDDO ENDDO C**************************************************** C TEMPERATURE/WAVELENGTH CORRECTIONS TO N2O ( I = 13) C**************************************************** DO LAMBDA = LYMAN + 11 , LYMAN + 24 SIGN2O(LAMBDA) = AN2O(1) DO J = 2 , 5 SIGN2O(LAMBDA) = SIGN2O(LAMBDA) + & AN2O(J)*(WAVE(LAMBDA)/10.0d0)**(J-1) ENDDO POWN2O(LAMBDA) = BN2O(1) DO J = 2 , 4 POWN2O(LAMBDA) = POWN2O(LAMBDA) + & BN2O(J)*(WAVE(LAMBDA)/10.0d0)**(J-1) ENDDO ENDDO C****************************************************** C TEMPERATURE/WAVELENGTH CORRECTIONS TO CH3CL ( I = 27) C****************************************************** DO LAMBDA = LYMAN + 18 , LYMAN + 19 TCH3CL(LAMBDA) = FCH3CL(1) DO J = 2 , 5 TCH3CL(LAMBDA) = TCH3CL(LAMBDA) + & FCH3CL(J)*WAVE(LAMBDA)**(J-1) ENDDO ENDDO C******************************************************** C TEMPERATURE/WAVELENGTH CORRECTIONS TO C2H3CL3 ( I = 34) C******************************************************** DO LAMBDA = LYMAN + 16 , LYMAN + 24 T295(LAMBDA) = F295(1) DO J = 2 , 5 T295(LAMBDA) = T295(LAMBDA) + & F295(J)*(WAVE(LAMBDA)/10.0d0)**(J-1) ENDDO T250(LAMBDA) = F250(1) DO J = 2 , 5 T250(LAMBDA) = T250(LAMBDA) + & F250(J)*(WAVE(LAMBDA)/10.0d0)**(J-1) ENDDO ENDDO C******************************************************* C TEMPERATURE/WAVELENGTH CORRECTIONS TO CLONO2 ( I = 40) C******************************************************* DO N = 1 , 5 L1 = LYMAN + LCLONO2(N) L2 = LYMAN + LCLONO2(N+1) DO LAMBDA = L1 , L2 TCLONO2(LAMBDA) = FCLONO2(1,N) DO J = 2 , 5 TCLONO2(LAMBDA) = TCLONO2(LAMBDA) + & FCLONO2(J,N)*WAVE(LAMBDA)**(J-1) ENDDO ENDDO ENDDO C********************************************************** C TEMPERATURE/WAVELENGTH CORRECTIONS TO H2CO ( I = 49 , 50) C********************************************************** DO N = 1 , 3 L1 = LYMAN + LH2CO(N) L2 = LYMAN + LH2CO(N+1) DO LAMBDA = L1 , L2 TH2CO(LAMBDA) = FH2CO(1,N) DO J = 2 , 5 TH2CO(LAMBDA) = TH2CO(LAMBDA) + & FH2CO(J,N)*WAVE(LAMBDA)**(J-1) ENDDO ENDDO ENDDO C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 100 FORMAT(//) 110 FORMAT(15X,6(E9.3,2X)) 120 FORMAT(F6.1,1X,F6.1,2X,E9.3) 130 FORMAT(I3) 140 FORMAT(E8.2 , 1X , F4.1 , 2X , F3.1) 150 FORMAT(E8.2 , 2X , F7.1) 160 FORMAT(6(E8.2,2X)) 170 FORMAT(6(G8.2,2X)) C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ RETURN END C------------------------------------------------------------------------------ C------------------------------------------------------------------------------ C------------------------------------------------------------------------------ DOUBLE PRECISION FUNCTION PCE(M,K) INCLUDE 'PARAMETERS.INC' INCLUDE 'PLANET.INC' COMMON & /ONED18/ NFORM(MB) , NDEST(MB) , & IFORM(IB,MB) , IDEST(IB,MB) & /BLOCK7/ RATE(IB,NALT) & /ONED19/ KFORM(IB,MB) , KDEST(IB,MB) & /ONED27/ MOL(6,IB) , KOEFF(6,IB) & /BLOCK3/ NOREAC(IB) & /BLOCK18/ CONC(MB,NALT) & /ONED20/ RCONST(IB,NALT) & /DEBUG/ RPROD(MB,NALT) C THIS MODULE COMPUTES THE CONCENTRATION OF A GIVEN SPECIES WHICH C IS IN PHOTOCHEMICAL EQUILIBRIUM. C RPROD AND RLOSS_FREQ ARE BOTH CALCULATED HERE. C THE LOSS FREQUENCY OF "A" IN THE REACTION A + B + C -> D + E + F C IS GIVEN BY K[B][C] , WHERE K IS THE REACTION CONSTANT RP = 0.0D0 NF = NFORM(M) DO I = 1 , NF IFRM = IFORM(I,M) RTEMP = RATE(IFRM,K)*KFORM(IFRM,M) RP = RP + RTEMP ENDDO RPROD(M,K) = RP C CALCULATE LOSS FREQUENCY - RFREQ SUM = 0.0D0 ND = NDEST(M) ! NUMBER OF REACTIONS WHICH DESTROY MOLECULE M DO I = 1 , ND IDST = IDEST(I,M) ! ITH REACTION NUMBER IN WHICH M IS DESTROYED NRCT = NOREAC(IDST) ! NUMBER OF REACTANTS IN REACTION # IDST FREQ = RCONST(IDST,K) DO J = 1 , NRCT IF( M .NE. MOL(J,IDST) ) THEN C0 = CONC(MOL(J,IDST),K) C1 = C0 DO KK = 1 , KOEFF(J,IDST)-1 ! MORE EFFICIENT WAY TO DO CONC**KOEFF C1 = C1*C0 ENDDO FREQ = FREQ*C1 ENDIF ENDDO SUM = SUM + FREQ*KDEST(IDST,M) ENDDO IF ((RP .EQ. 0.0D0) .OR. (SUM .EQ. 0.0D0)) THEN PCE = 0.0D0 !SET TO ZERO IF NO PRODUCTION OR LOSS RETURN ENDIF IF (K .EQ. 10) THEN C PRINT*,'IN PCE',M,RP,SUM ENDIF PCE = RP/SUM ! PHOTOCHEMICAL EQUILIB = P/L RETURN END C------------------------------------------------------------------------------ C------------------------------------------------------------------------------ C------------------------------------------------------------------------------ SUBROUTINE MODATM INCLUDE 'PARAMETERS.INC' INCLUDE 'PLANET.INC' COMMON & /SWITCHES/ ZLAT , DELTAT , DELT1 , DAYSTART , & TIMESTART , TIMESTOP , OUTTIME , & MODE , NTIME , NITER , IDIFF , istepdump & /BLOCK18/ CONC(MB,NALT) & /ONED36/ DAY & /BLOCK10/ Z(NALT) & /BLOCK11/ T(NALT) & /BLOCK13/ PRESS(NALT) & /ONED28/ DEN(NALT) & /ONED42/ NMOL & /ATMRATS/ O2RAT(NALT), ORAT(NALT), ON2RAT(NALT) C THIS ROUTINE CALLS THE MAP/MSIS ATMOSPHERE TO GET DENSITY AND C TEMPERATURE FOR A GIVEN DAY OF THE YEAR. CURRENTLY SET FOR 3PM C AT GREENWICH, WITH F107=150 AND AP=6 (AVE VALUES) C EDDY DIFFUSION IS READ IN FROM A SEPARATE FILE. REAL TEMP(2) , D(9) IYD = DAY + 81 !DAY IS DEFINED REL TO EQUINOX SEC = 5.4d4 !FOR 3 PM XLONG = 0.0d0 !GREENWICH MERIDIAN STL = SEC/3600.0d0 + XLONG/15.0d0 F107A = 150.0d0 F107D = 150.0d0 AP = 6.0d0 DO K = 1 , NALT CALL GTD7(IYD,SNGL(SEC),SNGL(Z(K)),SNGL(ZLAT), & SNGL(XLONG),SNGL(STL),SNGL(F107A),SNGL(F107D), & SNGL(AP),48,D,TEMP) C STORE OLD DENSITY WHICH WAS READ IN. THIS WILL BE USED IN FIXING C THE VALUES OF THE MIXING RATIOS. DenOld = 0.0d0 DO M = 1 , NMOL-3 DenOld = DenOld + CONC(M,K) ENDDO C D(3) = N2 , D(4) = O2 , D(5) = ARGON DEN(K) = D(3) + D(4) + D(5) + D(2) T(K) = TEMP(2) PRESS(K) = DEN(K)*(1.38d-19*T(K)) C FIX MIXING RATIOS INSTEAD OF FIXING THE CONCENTRATION. THIS WILL C ALSO SCALE THE SPECIES WHICH VARY, BUT THEY WILL CHANGE ANYWAY DO M = 1 , NMOL-3 CONC(M,K) = CONC(M,K)*(DEN(K)/DenOld) ENDDO CONC(NMOL,K)=DEN(K) CONC(1,K) = D(2) !LOAD THE ATOMIC O DENS CONC(6,K) = D(3) !LOAD THE NITROGEN DENSITY CONC(3,K) = D(4) !LOAD THE MOLEC OXYGEN DENS O2RAT(K) = D(4) / DEN(K) !AJK- careful, NRLMSIS00 O goes to zero below 75 km ORAT(K) = D(2) / DEN(K) !AJK ON2RAT(K) = D(3) / DEN(K) !AJK ENDDO RETURN END C------------------------------------------------------------------------------ C------------------------------------------------------------------------------ C------------------------------------------------------------------------------