C PROGRAM MAIN C MAIN COMMON ASSEMBLY FOR SPEC TRANS 2-D MODEL -- MAIN_BASE8OA.F- C has new PPM transport (TPCORE2D) call C also writes out final 360 days of mix rat at top level #include "com2d.h" COMMON/SPSH/ROCKALT(9),ROCK25(9),ROCK35(9),SSEM(9),T4EM(9), * RKHCL25(Z$),RKHCL35(Z$) C INSERT SUBSONIC/SUPERSONIC 12/24/92 COMMON/EMISSION/H2OEM(L$,Z$), XNOXEM(L$,Z$), HCEM(L$,Z$), 1 COEM(L$,Z$) common/volc/volcano C ------------------------- SPE RUN ------------------------- COMMON/SPE/AIONPAIR,AIONHOX,AIONNOX,IYRSPE,FRCSPE(18), * ALTSPE(56),ALTHOX(56),PAIRION(40,365,60), * PAIRION360(40,360,60) COMMON/SCUV/F10P7(427),SCVAR(24,427),SCDATE(427) REAL H2OJET,CH4JET,NOXJET,COJET c common/bcnox/ibc(l$,z$),bcclono2(l$,z$),bcbrono2(l$,z$) c common/noxc/bcclo24(l$,z$),bcn2o5(l$,z$),cno2(l$,z$), c c bcbro24(l$,z$),gamloc(l$,z$),zzz(20,l$,z$),yy(4,l$,z$), c c xxz(4,l$,z$),alpha(l$,z$),betaz(l$,z$),isw2(l$,z$), c c isw(l$,z$),icount(l$,z$) real extime,cputime(2) dimension conv(l$,z$,t$),xsave(l$,z$),OUTO3(L$,Z$) DIMENSION NEWSP(30,18,Z$),OUTSP(18,Z$),JNEW(30), * PHOTJOUT(20,L$,Z$),WAVELOW(39),WAVEHIGH(39),COLO3(18,360) REAL*8 xpi, qp1, zk, dy0, dz0, cosc(18), cose(19), XMR(L$,Z$58,T$) LOGICAL LDYNOUT c get inputs CALL CONTROL c galactic cosmic rays CALL GCRS c aerosol data - NOTE read both distributions (appropriate for 6 mos) CALL AEROSOLS c reaction rates CALL RCREAD c fixed ozone? YSO3FIX=.false. if(YSO3FIX)CALL O3SBUV c aircraft fleets? c YSHSCT=.true. c if(YSHSCT)CALL EMISSREAD c space shuttle and rocket launches CALL SPSHUTIN c get boundary conditions CALL BCIN C set up some contants by calling input, NOW also SETS up P and PRESS arrays (in COMMON) -- BASE8A CALL INPUT c CALL SPEREAD C READ IN SPEs - 6/14/95 CALL SCUVRD C READ IN SOLAR CYCLE UV VARIATION (FROM S. CHANDRA) - 1/5/96 C READ INITIAL PROFILES BY CALLING UNDUMP. CALL UNDUMP C write(33)c C C Initialize ozone column density array DO 650 ID=1,360 DO 650 IJ=1,L$ COLO3(IJ,ID)=0.0E0 650 CONTINUE C SET IYRLIFE FOR LIFETIME CALCULATION IYRLIFE=IYR c Note: DAY360 now goes 1.5-360.5, so IDAY360 => 1-360 DAY=(365./360.)*DAY360 DAYCHEK=DAY360 IDAY360=INT(DAY360) DT=DT0*365./360. C GET PHOTOCHEMICAL INPUT DATA CALL SOLFLIN CALL CROSSIN c print *,' after call crossin' C READ IN TEMPERATURES, dynamics look-up tables (10-day means) for new circulation 8/95 CALL TEMPIN C READ IN TABLE OF REDUCED SOLAR FLUXES GENERAGED BY RANDY KAWA - 1/24/95 CALL PHOTIN C SET UP DYNAMICS ARRAYS for all times (the 36 10-day means), stored in COMMON C TROPKZZ sets up tropospheric Kzz's, and computes Kyzall from tempall and kyyall arrays do 445 im10 = 1,36 CALL STREAMF(im10) CALL TROPKZZ(im10) 445 continue C CALL DYNOUT to write out all dynamics arrays for the full year if requested (to unit 37) LDYNOUT = .TRUE. c LDYNOUT = .FALSE. IF (LDYNOUT) CALL DYNOUT C Call SETDAILY -- interpolates dynamics to daily values, and sets up aerosol and major species, ZKM, etc. C SETDAILY keys on IDAY360 CALL SETDAILY C READ OTHER INPUTS (???**) C CALL RDATA CALL REACTION CALL COLDEN C Initialize column ozone density into array on first day and store DO 1550 IJ=1,L$ ICOL=IDAY360 IF(ICOL.LT.1)GO TO 1550 IF(ICOL.GT.360)THEN WRITE(6,2100)ICOL,IYR,IDAY360 GO TO 1550 ENDIF COLO3(IJ,ICOL)=NCOLGD(2,IJ,1) 1550 CONTINUE C C Initialize new 58 level arrays of transported species only with initial profiles, C keep constant mixing ratio above level 46 to start C do 2150 ik=1,Z$ do 2150 ij=1,L$ do 2150 ist=1,ITRANS c58(ist,ij,ik) = c(inttra(ist),ij,ik) 2150 cn58(ist,ij,ik) = c(inttra(ist),ij,ik) do 2151 ik=Z$+1,Z$58 do 2151 ij=1,L$ do 2151 ist=1,ITRANS c58(ist,ij,ik) = c(inttra(ist),ij,46)/m(ij,46)*m(ij,ik) 2151 cn58(ist,ij,ik) = c(inttra(ist),ij,46)/m(ij,46)*m(ij,ik) c time loop ccj ict = 0 DO 100 I=1,INDAYS c DO 100 I=1,1 c we are assuming we are using a 1 day timestep for chemistry, also in itday parameter ITDAY = I DAY360=DAY360+1. IF(DAY360 .GT. (YRL+1.))THEN DAY360=DAY360-YRL WRITE(14)COLO3 IYR=IYR+1 ENDIF DAY=(365./360.)*DAY360 DYT360=DAY360-0.5+IYR*360. IDAY360=INT(DAY360) c print *,' timestep, day360, day, iyr, ict=', c * i,day360,day,iyr,ict c set up boundary conditions if this is a time dependent run IF(.NOT.LSTSTATE)CALL BOUNDC c type *,' iday360, iyr = ',IDAY360, IYR if (iday360 .eq. 1) type *,' iday360, iyr = ',IDAY360, IYR C Call SETDAILY -- interpolates dynamics to daily values, and sets up aerosol and major species, ZKM, etc. volcano=1.0 reali=i if(i.gt.8220 .and. i.le.11520) * volcano=(1.0 + 9.0*exp(-(reali-8220.)/360.)) c Put in El Chichon eruption (April 4, 1982) if(i.gt.11520 .and. i.lt.15120) * volcano=(1.0 + 29.0*exp(-(reali-11520.)/360.)) c Put in Mt. Pinatubo eruption (June 1, 1991) CALL SETDAILY C Now CALL REACTION every time step to define reaction rates and het. gammas every C day with changes in temperature and water vapor CALL REACTION CALL COLDEN C Put column ozone density into array each day and store DO 550 IJ=1,L$ ICOL=IDAY360+1 ICOL=IDAY360 IF(ICOL.LT.1)THEN ICOL=360 ENDIF IF(ICOL.GT.360)THEN WRITE(6,2100)ICOL,IYR,IDAY360 2100 FORMAT(' ICOL=',I5,' IYR=',I3,' IDAY360=',I5) GO TO 550 ENDIF c if(ij.eq.1)print *,' iday360(550)=',iday360, c * ' icol=',icol,' day360=',day360 COLO3(IJ,ICOL)=NCOLGD(2,IJ,1) 550 CONTINUE C with frequency IPF, get photolysis rates, IPF=10 days c **** Now we'll update photoloysis rates everyday ***** c ccj if(ict .eq. 0)then ccj ict=ict+1 CALL SETUP c type *,'setup' CALL MULTSC c type *,'multsc' CALL RADT c type *,'radt' DO 1000 IJ=1,L$ DO 1000 IK=1,Z$ CALL FILL(IJ,IK) CALL CROSEC(IJ,IK) CALL SPDR(IJ,IK) 1000 CONTINUE c print *,' photolysis calculation completed at day=',day360 DO 4010 II=1,14 C WRITE(31,4040)II,DAY360 4040 FORMAT(' II=',I5,' DAY360=',F7.2) DO 4010 IJ=1,L$ C WRITE(31,4050)(PHOTJOUT(II,IJ,IK),IK=1,Z$) 4050 FORMAT(1P6E13.6) 4010 CONTINUE ccj ccj ELSE ccj ict=ict+1 ccj if(ict .eq. ipf)ict=0 ccj ccj ENDIF ccj ccj print *,' i, ict, =',i,ict c C ADVECTION - PPM TRANSPORT ROUTINE ************************ c c TPCORE2D uses mixing ratio as input xmr(18,58,T$) is in COMMON, and converts to density in routine c Note, cn58 and c58 are just transported species, cn58, c58 (T$,18,58) DO 5700 IT1=1,itrans DO 5700 IK=1,Z$ DO 5700 IJ=1,L$ xmr(IJ,IK,IT1) = c(inttra(IT1),IJ,IK)/M(IJ,IK) cdbc 01/10/95 correct ClONO2 and N2O5 from daytime average to 24 hour cdbc average value for proper transport: if(it1.eq.26.or.it1.eq.27.and. c rnday(inttra(it1),ij,ik).gt.0.) xmr(ij,ik,it1) = c c(inttra(it1),ij,ik)/M(IJ,IK) c *(tfd(ij)+(1.-tfd(ij))*rnday(inttra(it1),ij,ik)) 5700 CONTINUE c Now extend transported species to 58 levels, use new c58 array c do 5701 IT1=1,itrans do 5701 IK=z$+1,z$58 do 5701 IJ=1,L$ xmr(IJ,IK,IT1) = c58(IT1,IJ,IK)/M(IJ,IK) 5701 CONTINUE c c with a 1/2 time step for advection, loop through 2x DT=DT/2.E0 do 8846 ie=1,2 CALL TPCORE2D 8846 continue DT=DT*2.E0 c convert back to density do 5900 IK=1,z$ do 5900 IJ=1,l$ do 5900 IT1=1,itrans cn(inttra(IT1),IJ,IK) = xmr(IJ,IK,IT1)*m(IJ,IK) 5900 CONTINUE do 5901 IK=1,z$58 do 5901 IJ=1,l$ do 5901 IT1=1,itrans cn58(IT1,IJ,IK) = xmr(IJ,IK,IT1)*m(IJ,IK) 5901 CONTINUE C c diffusion should operate on the intermediary!! DO 6100 IK=1,Z$ DO 6100 IJ=1,L$ DO 6100 IT1=1,ITRANS ST(IT1,IJ,IK)=cn(inttra(it1),IJ,IK)/m(IJ,IK) 6100 CONTINUE C Extend diffusion array up to 58 levels DO 6101 IK=Z$+1,Z$58 DO 6101 IJ=1,L$ DO 6101 IT1=1,ITRANS 6101 ST(IT1,IJ,IK) = cn58(it1,IJ,IK)/m(IJ,IK) c c calculate diffusion terms for each transported constituent c from initial grid (CALL EACH DAY), NOW DONE FOR 58 LEVELS C call newdif1 C write(33)difn c add space shuttle and rocket contribution to hcl CALL SPSHUT do 900 ik=1,Z$ do 900 ij=1,L$ c if(tfd(IJ) .ge. 0.1 .and. tfd(IJ) .le. .999) then c call ncdn(IJ,IK) c else c if(tfd(IJ) .lt. 0.1)call ncnight(IJ,IK) c c if(tfd(IJ) .gt. .999)call ncday(IJ,IK) c endif CALL SPESET(IJ,IK) C Compute NOx and HOx from SPEs - 6/14/95 call solver(IJ,IK) 900 continue c Now call separate SOLVER for levels 47-58 using skeleton chemistry (no meat) do 920 ik=Z$+1,Z$58 do 920 ij=1,L$ call solver58(IJ,IK) 920 continue C Note: ISESN is not used in HETCHEM at present, but will give it a value here if it's needed in the future isesn = int((iday360+29)/30) call alumaer call hetchem(isesn,i,iday360) C write(33)cn CALL RAINOUT IF(IYR.NE.IYRLIFE)YSLFPRNT=.TRUE. CALL LIFETIME c type *,' after call lifetime' IF(IYR.NE.IYRLIFE)THEN IYRLIFE=IYR YSLFPRNT=.FALSE. ENDIF c update constituents, check for NaN's do 950 ik=1,Z$ do 950 ij=1,L$ do 950 is=1,S$ c(is,ij,ik) = cn(is,ij,ik) ccdbc 01/10/95 correct ClONO2 and N2O5 from daytime average to 24 hour ccdbc average value for proper transport: c if(is.eq.30.or.is.eq.8)c(is,ij,ik) c c =cn(is,ij,ik)*(tfd(ij)+(1.-tfd(ij))*rnday(is,ij,ik)) if(c(is,ij,ik) .ne. c(is,ij,ik)) then write(6,866) is, ij, ik, c(is,ij,ik), i 866 format(1x,' is ij ik =',3i6,4x,'c = ',1PE12.3,6x,'i(day) =',I5) endif if (ij .ge. 2) then if ( c(is,ij,ik) .ne. c(is,ij,ik) .and. * c(is,ij-1,ik) .ne. c(is,ij-1,ik)) goto 8888 endif 950 continue c also update constituents for new c58 above and below 90 km and check for NaN's, load in c c58 array below level 46 just for cosmetic puposes and for output do 970 ik=1,Z$58 do 970 ij=1,L$ do 970 is=1,ITRANS if (ik .le. 46) then c58(is,ij,ik) = cn(inttra(is),ij,ik) cn58(is,ij,ik) = cn(inttra(is),ij,ik) endif if (ik .ge. 47) c58(is,ij,ik) = cn58(is,ij,ik) if(c58(is,ij,ik) .ne. c58(is,ij,ik)) then write(6,866) is, ij, ik, c58(is,ij,ik), i endif if (ij .ge. 2) then if ( c58(is,ij,ik) .ne. c58(is,ij,ik) .and. * c58(is,ij-1,ik) .ne. c58(is,ij-1,ik)) goto 8888 endif c load in mixing ratio array for output mr58(is,ij,ik) = cn58(is,ij,ik)/m(ij,ik) 970 continue C C write out to fort.33 and fort.14 on final year of run every 30 days at day 15,45,75,... C and write out extended array (cn58 - up to level 58) of species to fort.34 C if (i .ge. indays-359) then if (mod(i-15,IOUTP) .eq. 0.0) then write(33) cn write(34) mr58 c write(35) ctloss end if end if C 100 CONTINUE C print *,' before call dump' CALL DUMP C print *,' after call dump' c put in Tom's output stuff WRITE(14)COLO3 IF(DAY360 .GT. 90.)THEN C Get printout of lifetime at the end here only if the model has run C for at least 90 days. YSLFPRNT=.TRUE. CALL LIFETIME ENDIF 8888 extime = etime(cputime) write(6,6741) cputime(1)/3600., cputime(2)/3600., extime/3600. 6741 format(1x,'USER CPU =', F9.2,4X,'SYS CPU =',F9.2,10X, * 'TOTAL CPU (hrs) =',F9.2) STOP END