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$) COMMON/ALOX/ALO25(9),ALO35(9),RKALO25(Z$),RKALO35(Z$), * AERAL(L$,Z$) COMMON/ALOXA/AERAL1(L$,Z$),AERAL2(L$,Z$),AERAL3(L$,Z$) C INSERT SUBSONIC/SUPERSONIC 12/24/92 COMMON/EMISSION/H2OEM(L$,Z$), XNOXEM(L$,Z$), HCEM(L$,Z$), 1 COEM(L$,Z$) 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 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 c 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. 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 - NEW: Compute Kzz diffusion separately in NEWDIFZ (array DIFNZ), with a 1/IKDT of C a day time step, to accomodate up to 100 m2/sec Kzz in mesosphere, IKDT = number of time steps per C day; calculate diffusion terms for each transported constituent, NOW DONE FOR 58 LEVELS ******* C IKDT = 8 DT = DT/real(IKDT) DO 1551 ied=1,IKDT DO 6150 IK=1,Z$ DO 6150 IJ=1,L$ DO 6150 IT1=1,ITRANS 6150 ST(IT1,IJ,IK) = cn(inttra(it1),IJ,IK)/m(IJ,IK) C Extend diffusion array up to 58 levels DO 6151 IK=Z$+1,Z$58 DO 6151 IJ=1,L$ DO 6151 IT1=1,ITRANS 6151 ST(IT1,IJ,IK) = cn58(it1,IJ,IK)/m(IJ,IK) CALL NEWDIFZ do 7800 IT1=1,itrans do 7810 IK=1,z$ do 7810 IJ=1,l$ 7810 cn(inttra(IT1),IJ,IK)= cn(inttra(IT1),IJ,IK)-DIFNZ(IT1,IJ,IK)*DT do 7850 IK=1,z$58 do 7850 IJ=1,l$ 7850 cn58(IT1,IJ,IK) = cn58(IT1,IJ,IK) - DIFNZ(IT1,IJ,IK)*DT 7800 CONTINUE 1551 continue DT = DT*real(IKDT) C NOW, compute Kyy term with the 1 day time step as before, in NEWDIFY, DO 6100 IK=1,Z$ DO 6100 IJ=1,L$ DO 6100 IT1=1,ITRANS 6100 ST(IT1,IJ,IK) = cn(inttra(it1),IJ,IK)/m(IJ,IK) 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) CALL NEWDIFY do 7700 IT1=1,itrans do 7710 IK=1,z$ do 7710 IJ=1,l$ 7710 cn(inttra(IT1),IJ,IK)= cn(inttra(IT1),IJ,IK)- DIFNY(IT1,IJ,IK)*DT do 7750 IK=1,z$58 do 7750 IJ=1,l$ 7750 cn58(IT1,IJ,IK) = cn58(IT1,IJ,IK) - DIFNY(IT1,IJ,IK)*DT 7700 CONTINUE C C Kyz NOT USED, if needed in future, its done separately in NEWDIFYZ C C CALL NEWDIFYZ C C do 4800 IT1=1,itrans C do 4810 IK=1,z$ C do 4810 IJ=1,l$ C 4810 cn(inttra(IT1),IJ,IK)=cn(inttra(IT1),IJ,IK)- DIFNYZ(IT1,IJ,IK)*DT C C do 4850 IK=1,z$58 C do 4850 IJ=1,l$ C 4850 cn58(IT1,IJ,IK) = cn58(IT1,IJ,IK) - DIFNYZ(IT1,IJ,IK)*DT C 4800 CONTINUE C 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 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 alumaera call alumaerb call alumaerc do 922 ik=1,Z$ do 922 ij=1,L$ aeral(ij,ik)=aeral1(ij,ik)+aeral2(ij,ik)+aeral3(ij,ik) 922 continue 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, also, Set limit to 1.E-12 number density c as was done in SOLVER, since HETCHEM may give zero's for Solid HNO3 and ice C (m at 115 km ~1.e12, and ~3.e19 at the ground, so the minimum mixing ratio will be ~3.E-31) do 950 ik=1,Z$ do 950 ij=1,L$ do 950 is=1,S$ IF (cn(is,ij,ik) .LT. 1.E-12) cn(is,ij,ik) = 1.E-12 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 (cn58(is,ij,ik) .LT. 1.E-12) cn58(is,ij,ik) = 1.E-12 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