PROGRAM ONED INCLUDE 'PARAMETERS.INC' LOGICAL ITERATE COMMON & /SWITCHES/ ZLAT , DELTAT , DELT1 , DAYSTART , & TIMESTART , TIMESTOP , OUTTIME , & MODE , NTIME , NITER , IDIFF , istepdump & /ONED1/ CONCO(MB,NALT) & /ONED2/ ITERATE & /ONED3/ XAU & /ONED5/ SRTOP , SSTOP & /ONED8/ TIME & /ONED10/ ICONVRG & /ONED14/ DELTNEW & /ONED15/ THOUR & /ONED17/ NTRANS & /BLOCK18/ CONC(MB,NALT) & /ONED25/ JUPDATE & /ONED28/ DEN(NALT) & /ONED35/ FLUX(NWAVE) & /ONED36/ DAY & /ONED37/ ITIME & /ONED38/ DELTOLD , NHALFS & /ONED40/ MTRANS(MB) & /ONED41/ SZA & /ONED42/ NMOL & /ONED53/ MPCE(MB) & /ONED54/ NPCE & /ONED56/ DAYNEW & /BLOCK10/ Z(NALT) & /BLOCK11/ T(NALT) & /BLOCK13/ PRESS(NALT) & /BLOCK27/ PI & /BLOCK29/ FLUX1AU(NWAVE) ! FROM SUBROUTINE PHOTO & /ONED62/ kprint(3) ! altitudes at which we output rates & /BLOCK7/ RATE(IB,NALT) & /ATMRATS/ O2RAT(NALT), ORAT(NALT), ON2RAT(NALT) C THIS IS THE MAIN DRIVER FOR THE CHEM1D CALCULATIONS. THIS SUBROUTINE UPDATES C THE TIME AND CALLS THE SUBROUTINE WHICH SOLVES A 1-D DIFFUSION EQUATION. C THIS EQUATION IS NON-LINEAR, SO IT IS ITERATED A NUMBER OF TIMES. IF THE C SOLUTION HAS NOT CONVERGED AFTER THESE ITERATIONS THEN THE TIME STEP IS C HALVED AND THE ITERATION PROCEDURE IS TRYED AGAIN. C C MODE = 0 FOR DIURNAL AVERAGE CALCULATION C 1 FOR DIURNAL CALCULATION C 2 FOR STEADY STATE CALCULATION C C NHALF = NUMBER OF TIMES THE TIMESTEP IS HALVED BEFORE GIVING UP C C NTIME = NUMBER OF TIME STEPS C C DELTNEW = TIMESTEP AS SEEN BY STEP C C DELTAT = TIMESTEP FOR DIURNAL CALCULATION C C ICONVG = CONVGERGENCE FLAG C 0 CONVERGED C 1 NO CONVERGENCE C C IDIFF = 0 NO DIFFUSE FLUX CONTRIBUTION C = 1 ADD DIFFUSE FLUX/GROUND ALBEDO INTO TRANSMISION OF FLUX C C IUPPER = TYPE OF UPPER BOUNDARY CONDITION FOR EACH SPECIES C C ILOWER = TYPE OF LOWER BOUNDARY CONDITION FOR EACH SPECIES C C============================================================================== C------------------------ INTIALIZE VALUES ------------------------------------ CALL INPUT TIME = TIMESTART ! SUBROUTINE STEP UPDATES THIS VALUE THOUR = TIME - 0.5d0*DELTAT/3600.0d0 DAY = DAYSTART DAYNEW = 0.0d0 ! FIRST DAY CALL MODATM ! MSIS ATMOSPHERE CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Swap MSIS P,D,T with MLS do i=1,NALT read(15,101) z(i), T(i), den(i), press(i) CONC(1,i) = ORAT(i) * den(i) ! AJK CONC(3,i) = O2RAT(i) * den(i) ! AJK CONC(6,i) = ON2RAT(i) * den(i) ! AJK enddo CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC c write(6,*) CONC(1,25), CONC(3,25), CONC(6,25), 'AJK' c write(6,*) orat(25), o2rat(25), on2rat(25) C============ play with model atmosphere ============ c write(10,111) c111 format(/,2x,'Model Atmosphere',/) c zz = 0.0d0 c dz = 2.0d0 do i=1,NALT write(10,101) z(i),T(i),den(i),press(i) c zz = zz + dz enddo 101 format (2x,f10.4,2X,f10.2,2x,2e12.4) c do i=1,NALT c read(21,102) zzz,T(i) c enddo c102 format (f10.3,2x,f12.3) c write(10,112) c112 format(/,2x,'Revised Model Atmosphere',/) c zz = 0.0d0 c dz = 2.0d0 c do i=1,66 c write(10,101) z(i),T(i),den(i),press(i) c zz = zz + dz c enddo C====================================================== CALL PHOTO CALL ADJUST CALL OUTPUT CALL ORBIT C Altitude indices to ouput detailed rates in a diurnal run kprint(1) = 24 ! 46 km kprint(2) = 26 ! 50 km kprint(3) = 28 ! 54 km C SCALE FLUX FOR DISTANCE FROM SUN. FAC = 1.0d0/(XAU**2) DO LAMBDA = 1 , NWAVE FLUX(LAMBDA) = FLUX1AU(LAMBDA)*FAC ENDDO CALL DAYTIME ! CALCULATE LENGTH OF DAY. IF( MODE .EQ. 0 ) CALL GEOMAVG ! DIURNAL AVERAGED SOLUTION IF( (MODE .EQ. 2) .AND. (THOUR .GE. SRTOP) .AND. & (THOUR .LE. SSTOP) ) CALL GEOM ! STEADY-STATE c WRITE (6,222) MODE , THOUR , SRTOP , SSTOP 222 FORMAT('MODE = ',G10.3,/, & 'THOUR = ',G10.3,/, & 'SRTOP = ',G10.3,/, & 'SSTOP = ',G10.3) IF( IDIFF .EQ. 1 ) CALL DIFFLUX ! CALCULATE MULTIPLE SCATTERING MATRIX ! ONLY NEEDS TO BE CALLED ONCE A DAY DELTNEW = DELTAT C------------------------------------------------------------------------------ iprint = 1 ! index for # calls to output C MAIN LOOP - INTEGRATE CONCENTRATIONS OVER NTIME STEPS. DO ITIME = 1 , NTIME ! LOOP OVER TIME C RESET ITERATION FLAGS FOR NEW TIME VALUE NHALFS = NHALF ! RESET NHALF FLAG JUPDATE = 0 ! ONLY USED IF REPEATED ITERATIONS OCCUR ITERATE = .TRUE. ! ITERATION FLAG IF( MODE .EQ. 1 ) THEN DELTNEW = DELTAT ELSE DELTNEW = DELTAT DELTOLD = DELTNEW CCC DELTNEW = 2*DELTNEW CALL SLANTAVG ! COMPUTE UPDATED COLUMN DENSITIES ENDIF C UPDATE TIME VALUES IF( (ITIME .GT. 1) .AND. (MODE .EQ. 1) ) CALL STEP C CALCULATE NEW GEOMETRIC FACTORS. IF( (MODE .EQ. 1) .AND. (THOUR .GE. SRTOP) .AND. & (THOUR .LE. SSTOP) ) CALL GEOM C UPDATE DIFFUSE FLUX MATRIX IF( (MODE .NE. 1 ) .AND. ( IDIFF .EQ. 1 ) .AND. & (ITIME .GT. 1 ) ) CALL DIFFLUX C SAVE ORIGINAL CONCENTRATIONS. THESE ARE REUSED IF THE TIME STEP IS HALVED. DO K = 1 , NALT DO M = 1 , NTRANS MVARY = MTRANS(M) CONCO(MVARY,K) = CONC(MVARY,K) ENDDO DO N = 1 , NPCE M = MPCE(N) CONCO(M,K) = CONC(M,K) ! SPECIES IN PHOTOCHEMICAL EQUILIBRIUM ENDDO CONCO(NMOL-3,K) = CONC(NMOL-3,K) ! ELECTRON DENSITY ENDDO C ITERATE OVER A FIXED TIME STEP. IF THE CALCULATIONS CONVERGE, THEN ITERATE C IS SET TO FALSE. OTHERWISE, THE TIME STEP IS HALVED, ITERATE IS SET TO TRUE C AND ANOTHER PATH THROUGH THE ITERATION LOOP IS DONE. DO WHILE( ITERATE .EQ. .TRUE. ) C CALL ADVECTION-DIFFUSION SOLVER AND CHECK FOR CONVERGENCE. ACTUAL ITERATION C IS PERFORMED WITHIN THIS ROUTINE. CALL NEWTON write (6,120) DELTNEW, DELTOLD, ICONVRG , TIME 120 format (2x,'Called NEWTON: DELTNEW = ',G10.3,', DELTOLD = ', & G10.3,', ICONVRG = ',I5,', TIME = ',G10.3) IF( ICONVRG .EQ. 0 ) THEN WRITE (6,100) THOUR , SZA*180.0d0/PI , time ITERATE = .FALSE. ! CONVERGED - NO MORE ITERATIONS ARE NEEDED ENDIF IF( ICONVRG .EQ. 1 ) THEN CALL TIMEHALF ! REDUCE TIMESTEP AND TRY AGAIN write (6,130) DELTNEW, DELTOLD, ICONVRG , TIME,ITIME 130 format (2x,'Called TIMEHALF: DELTNEW = ',G10.3, & ', DELTOLD = ',G10.3,', ICONVRG = ',I3,', TIME = ',G9.3,I3) IF( MODE .EQ. 1 ) CALL GEOM ! RECALCULATE GEOMETRIC FACTORS ! APPROPRIATE TO NEW SZA C RESET CONCENTRATION TO ORIGINAL VALUE DO K = 1 , NALT DO M = 1 , NTRANS MVARY = MTRANS(M) CONC(MVARY,K) = CONCO(MVARY,K) ENDDO DO N = 1 , NPCE M = MPCE(N) CONC(M,K) = CONCO(M,K) ! SPECIES IN PHOTOCHEMICAL EQUILIBRIUM ENDDO CONC(NMOL-3,K) = CONCO(NMOL-3,K) ! ELECTRON DENSITY ENDDO ENDIF ENDDO ! *END LOOP OVER ITERATE FLAG* IF( ICONVRG .EQ. 1 ) THEN ! NO CONVERGENCE, GIVING UP IF( MODE .EQ. 0 ) THEN CALL CHEMOUT( NTIME,1 ) ! OUTPUT LAST CONVERGED RESULT ENDIF WRITE(6,110) GOTO 10 ENDIF IF( TIME .GE. TIMESTOP) THEN CALL CHEMOUT( NTIME ,1) PRINT * , 'TIME > TIMESTOP, MODEL COMPLETED',time,timestop GOTO 10 ENDIF C WRITE RESULTS OF A CONVERGED CALCULATION CCC if( SZA*180.0d0/PI .lt. 90.0d0 ) then IF( MODE .NE. 0 .and. time .ge. outtime ) then write (6,140) ITIME 140 format (2x,'Calling CHEMOUT: ITIME = 'I5) ccccc & 2x,'atox, 78 km=',e10.3) CALL CHEMOUT( ITIME , iprint ) iprint = iprint + 1 endif CCC endif ENDDO ! END LOOP OVER ITIME IF( MODE .EQ. 0 ) CALL CHEMOUT( NTIME,1) print * , 'O.K. I am done' 100 FORMAT( '***>>>CONVERGED<<<***' , 2X , 'THOUR = ', G10.3 , & 2X , 'SZA = ' , F6.2 , 2X , ' DEGREES,' , & ' time = ' , G10.3 ) 110 FORMAT( '-->***NO CONVERGENCE - GIVING UP!***<--' ) 10 STOP END C------------------------------------------------------------------------------ C------------------------------------------------------------------------------ C------------------------------------------------------------------------------ SUBROUTINE NEWTON INCLUDE 'PARAMETERS.INC' COMMON & /SWITCHES/ ZLAT , DELTAT , DELT1 , DAYSTART , & TIMESTART , TIMESTOP , OUTTIME , & MODE , NTIME , NITER , IDIFF , istepdump & /ONED1/ CONCO(MB,NALT) & /ONED5/ SRTOP , SSTOP & /ONED7/ A(MB,NALT) , C(MB,NALT) & /ONED9/ DCONC(MB,NALT) & /ONED10/ ICONVRG & /ONED14/ DELTNEW & /ONED15/ THOUR & /ONED17/ NTRANS & /ONED22/ RNET(MB,NALT) ! FROM SUBROUTINE CHEM & /BLOCK18/ CONC(MB,NALT) & /ONED25/ JUPDATE & /ONED30/ B(MB,NALT) & /ONED39/ RESID(MB,NALT) & /ONED40/ MTRANS(MB) & /ONED42/ NMOL & /ONED53/ MPCE(MB) & /ONED54/ NPCE C THIS SUBROUTINE SOLVES FOR THE NEW LEVELS OF CONCENTRATION. IT DOES THIS BY C SOLVING THE ASSOCIATED DIFFUSION-TYPE EQUATION. SINCE THE SOURCE TERM C ( PRODUCTION - LOSS ) IS A NONLINEAR FUNCTION OF THE CONCENTRATION, THE C SOLUTION MUST BE OBTAINED BY AN ITERATIVE PROCEEDURE. THE INITAL GUESS FOR C THE CONCENTRATION IS TAKEN TO BE THE VALUE FROM THE PREVIOUS TIME STEP AND IS C THEN ITERATED TO CONVERGENCE. THE RESULTING VALUE IS THE CONCENTRATION AT THE C NEW VALUE OF TIME. THIS ROUTINE USES THE NEWTON-RAPHSON ALGORITHM TO PERFORM C THIS ITERATION. DO ITER = 1 , NITER C PHOTOCHEMICAL RATE CONSTANTS IF( (JUPDATE .EQ. 1 ) .OR. (ITER .EQ. 1) .OR. & ( MODE .EQ. 0 )) THEN IF( MODE .NE. 0) THEN ! DIURNAL RATE CONSTANTS IF( (THOUR .GE. SRTOP) .AND. (THOUR .LE. SSTOP) ) THEN CALL RJDIURN CALL SRB ! SCHUMANN-RUNGE BAND TRANSMITANCE CALL DIULYM ! LYMAN ALPHA TRANSMITANCE ENDIF ELSE CALL RJAVG ! DIURNALLY AVERAGED RATE CONSTANTS CALL SRBAVG ! SCHUMANN-RUNGE BAND TRANSMITANCE CALL LYMAVG ! LYMAN ALPHA TRANSMITANCE ENDIF CALL RATEKIN ! KINETIC RATE CONSTANTS ENDIF C COMPUTE RHS OF CONTINUITY EQUATION USING UPDATED RATE CONSTANTS CALL PMINUSL ! RETURNS RNET(M,K) IF( IMARCH .EQ. 0 ) CALL COEFF C CALCULATE RESIDUE IF( IMARCH .EQ. 0) THEN c$doacross c$& local(k,m,mvary) , c$& share(mtrans,resid,conc,conco,rnet,deltnew,a,b,c) DO K = 1 , NALT-2 DO M = 1 , NTRANS MVARY = MTRANS(M) RESID(M,K+1) = -0.5d0*(A(M,K+1)*CONC(MVARY,K) + & B(M,K+1)*CONC(MVARY,K+1) + & C(M,K+1)*CONC(MVARY,K+2) ) - & 0.5d0*(A(M,K+1)*CONCO(MVARY,K) + & B(M,K+1)*CONCO(MVARY,K+1) + & C(M,K+1)*CONCO(MVARY,K+2) ) + & (CONCO(MVARY,K+1) - & CONC(MVARY,K+1))/DELTNEW + & RNET(M,K+1) ENDDO ENDDO ELSE DO K = 1 , NALT DO M = 1 , NTRANS MVARY = MTRANS(M) RESID(M,K) = (CONCO(MVARY,K) - & CONC(MVARY,K))/DELTNEW + & RNET(M,K) ENDDO ENDDO ENDIF C COMPUTE JACOBIAN FOR COUPLED CHEMISTRY CALL BMATRIX C SOLVE SYSTEM OF EQUATIONS FOR DELTA_CONC IF( IMARCH .EQ. 0 ) THEN CALL BNDRY ! SET UP BOUNDARY CONDITIONS CALL THOMAS ! RETURNS DCONC ELSE CALL DIAG ! RETURNS DCONC ENDIF C SOLUTION UPDATED FOR THIS ITERATION DO K = 1 , NALT DO M = 1 , NTRANS MVARY = MTRANS(M) CONC(MVARY,K) = CONC(MVARY,K) + DCONC(MVARY,K) IF( CONC(MVARY,K) .LT. -TINY ) THEN DUM = PCE(MVARY,K) print* , "MVARY = " ,MVARY , "K = " , k, & " CONC = " ,CONC(MVARY,K),CONC(MVARY,K+1) C IF (K .EQ. 1) THEN !MANY NEGATIVE VALUES ARE AT THE GROUND C CONC(MVARY,1) = CONC(MVARY,2) !TRY AND KLUGE AROUND THIS C print*,'trying to kluge this in newton' C print*,dum,conc(mvary,1),conc(mvary,2),conc(mvary,3) C CONC(MVARY,1) = DUM C DCONC(MVARY,1) = 0.0 C ELSE ICONVRG = 1 ! NEGATIVE VALUES -- CUT TIME STEP AND TRY AGAIN C ENDIF RETURN ENDIF CONC(MVARY,K) = DMAX1(CONC(MVARY,K),0.0D0) ! SETS TINY VALUES TO 0 ENDDO DO N = 1 , NPCE M = MPCE(N) CONC(M,K) = PCE(M,K) ! SPECIES IN PHOTOCHEMICAL EQUILIBRIUM IF ((K.EQ. 10) .and. (M .EQ. 2)) PRINT*,'AFTER PCE',N,M,K,CONC(M,K) ENDDO CONC(NMOL-3,K) = CONC(NMOL-4,K) + CONC(NMOL-5,K) + & CONC(NMOL-6,K) ! ELECTRON DENSITY = SUM OVER IONS ENDDO C ASCERTAIN CONVERGENCE OF DIFFUS CALCULATION VIA C FRACTIONAL-CHANGE-IN-CONCENTRATION TEST. CALL CONVRG(ITER) IF (ICONVRG .EQ. 0) THEN RETURN ! CALCULATIONS CONVERGED. RETURN TO MAIN PROGRAM AND ! ADVANCE TIME STEP ENDIF ENDDO ! END LOOP OVER ITER C CALCULATIONS DID NOT CONVERGE! (ICONVRG = 1) RETURN TO MAIN PROGRAM AND C REDUCE THE TIME STEP RETURN END C------------------------------------------------------------------------------ C------------------------------------------------------------------------------ C------------------------------------------------------------------------------ SUBROUTINE CONVRG(ITER) INCLUDE 'PARAMETERS.INC' CHARACTER*8 AMOL(MB) COMMON & /SWITCHES/ ZLAT , DELTAT , DELT1 , DAYSTART , & TIMESTART , TIMESTOP , OUTTIME , & MODE , NTIME , NITER , IDIFF , istepdump & /ONED9/ DCONC(MB,NALT) & /ONED10/ ICONVRG & /ONED14/ DELTNEW & /ONED15/ THOUR & /ONED17/ NTRANS & /BLOCK10/ Z(NALT) & /BLOCK18/ CONC(MB,NALT) & /ONED40/ MTRANS(MB) & /ONED43/ AMOL C CHECK FOR CONVERGENCE OF CALCULATION BY DETERMINING FRACTIONAL C CHANGE IN VALUES BETWEEN TWO SUCCESSIVE ITERATIONS. ICONVRG = 0 DO K = 1 , NALT DO M = 1 , NTRANS MVARY = MTRANS(M) IF( CONC(MVARY,K) .GE. TINY ) THEN IF (ABS( DCONC(MVARY,K)/CONC(MVARY,K) ) .GT. ERRMAX) THEN ICONVRG = 1 IF(ITER .EQ. NITER) THEN WRITE(6,100) AMOL(MVARY) , Z(K) , THOUR , DELTNEW WRITE(6,110) CONC(MVARY,K) , DCONC(MVARY,K) ENDIF ENDIF ENDIF ENDDO ENDDO 100 FORMAT(A8,' DID NOT CONVERGE AT ALTITUDE ' , F6.2 , ' AND TIME ' , & F6.3 , ' WITH TIME STEP = ' , 1PE10.3 ) 110 FORMAT('CONC = ' , 1PE10.3 , 2X , 'DELTA_CONC = ' , 1PE10.3 , /) RETURN END