SUBROUTINE SETUP #include "com2d.h" COMMON/SCFLUX/SCF(IL$) COMMON/SCUV/F10P7(427),SCVAR(24,427),SCDATE(427) DIMENSION RLAMB(6) DATA RLAMB/1700.,1900.,2100.,2400.,2600.,3000./ C DECLINATION, ZENITH ANGLE, C CHAPMAN FUNCTION AT LATITUDE LATP AND ELAPSED TIME T. C RESET DECLINATION WITH TIME IF REQUESTED. FORMULA AND C CONSTANTS ARE FROM THE ASTRONOMICAL ALMANAC, WITH C MEANL=MEAN SOLAR LONGITUDE (CORRECTED FOR ABERRATION), C MEANA=MEAN ANOMALY, ECLIPL=ECLIPTIC LONGITUDE. MEANL=DEC1+DEC2*DAY MEANA=DEC3+DEC4*DAY ECLIPL=MEANL+DEC5*SIN(MEANA)+DEC6*SIN(2.*MEANA) SIND=DEC7*SIN(ECLIPL) R=ASIN(SIND) COSD=COS(R) DECD=R*RTD c print *,' dec1,dec2,dec3,dec4,dec5,dec6,meanl,meana, c * eclipl,dec7,dec8,rtd,r',dec1,dec2,dec3,dec4,dec5, c * dec6,meanl,meana,eclipl,dec7,dec8,rtd,r c print *,' day,decd,sind,cosd=',day,decd,sind,cosd C CORRECT SOLAR FLUXES AT THE TOP OF THE ATMOSPHERE FOR C THE CURRENT ORBITAL POSITION OF THE EARTH. CORB=1.0/(1.0-DEC8*COS(MEANA))**2 IYRMIN=5.5 IYRMAX1=0.0 IYRMAX2=11.0 PIYR=IYR PIYRUSE=PIYR + (DAY360/360.) - 1. C THE 1 CALIBRATES IT TO HAVE SOLAR MAX IN 1980 (45), 1969 (34), 1958 (23), C 1947 (12), AND 1936 (1). 3/6/90 PIYRSC=PIYR + (DAY360/360.) c PIYRSC=PIYR + (DAY360/360.) + .5 C C THIS NEXT SECTION REDUCES THE YEAR DOWN TO A SINGLE SOLAR CYCLE C DO 200 JL=1,IL$ SOLCYC=1.0 FMAX=1.0 FMIN=1.0 C DATA RLAMB/1700.,1900.,2100.,2400.,2600.,3000./ IF(JL .EQ. 1)THEN FMAX=1.3 FMIN=0.7 END IF IF(WVL(JL).GE.RLAMB(1) .AND. WVL(JL).LT.RLAMB(2))THEN FMAX=1.05 FMIN=0.95 END IF IF(WVL(JL).GE.RLAMB(2) .AND. WVL(JL).LT.RLAMB(3))THEN FMAX=1.03 FMIN=0.97 END IF IF(WVL(JL).GE.RLAMB(3) .AND. WVL(JL).LT.RLAMB(4))THEN FMAX=1.01 FMIN=0.99 END IF IF(WVL(JL).GE.RLAMB(4) .AND. WVL(JL).LT.RLAMB(5))THEN FMAX=1.01 FMIN=0.99 END IF IF(WVL(JL).GE.RLAMB(5) .AND. WVL(JL).LT.RLAMB(6))THEN FMAX=1.00 FMIN=1.00 END IF FMEAN=1.0 SOLCYC=FMEAN + (FMAX-FMEAN)*COS(PIYRUSE * (2.*PI/11.)) if(piyrsc .gt. scdate(427))go to 2211 if(jl .le. 24)then if(piyrsc .ge. 25.)then if(piyrsc .lt. scdate(1))then percuv=scvar(jl,1) solcyc=1. + (percuv/100.) c if(jl .eq. 19)then c print *,' scvar=',scvar(jl,1) c endif if(1 .gt. 0)go to 2211 endif do 2233 iscd=2,427 itemp=iscd if(piyrsc .le. scdate(iscd))go to 2255 2233 continue 2255 fdelta=scdate(itemp) - scdate(itemp-1) fscdup=(scdate(itemp) - piyrsc)/fdelta fscdlow=(piyrsc-scdate(itemp-1))/fdelta percuv=fscdup*scvar(jl,itemp-1) + * fscdlow*scvar(jl,itemp) solcyc=1. + (percuv/100.) c if(jl .eq. 19)then c print *,' percuv=',percuv,' fscdup=',fscdup, c * ' fscdlow=',fscdlow,' scvar(it=',scvar(jl,itemp), c * ' scvar(it-1=',scvar(jl,itemp-1) c endif endif endif 2211 continue c if(jl.eq.1) c * WRITE(6,8290)SOLCYC,IYR,PIYRUSE,FMAX,FMIN,jl,wvl(jl), c * corb 8290 FORMAT(' SOLCYC=',1PE13.4,' IYR=',I5,' PIYRUSE=',1PE13.4, * /,' FMAX=',1PE13.4,' FMIN=',1PE13.4,' jl=',i3,' wvl=', * 1pe13.4,' corb=',1pe13.4) c if(jl .eq. 19)then c print *,' piyrsc=',piyrsc,' solcyc=',solcyc c endif SOLCYC = 1.0 SCF(JL)=CORB*SOLCYC FLUX(JL)=CORB*SOLCYC*FLUX0(JL) c if(jl.eq.1)print *,' flux=',flux(jl),' flux0=',flux0(jl) 200 continue C DO 2100 IJ=1,L$ DO 2100 IK=1,Z$ GCRMEAN=(GCRMIN(IJ,IK)+GCRMAX(IJ,IK))/2. GCR(IJ,IK)=GCRMEAN 2100 CONTINUE DO 3000 IJ=1,L$ CALL FILL(IJ,Z$) SINL=SIN(LATP*DTR) COSL=COS(LATP*DTR)+1.e-20 c print *,' ij latp cosl sinl',ij,latp,cosl,sinl C SOLAR ZENITH ANGLE CHID. COSCB=SIND*SINL+COSD*COSL*COS(OMEGAD*DTR*(CLOCK+12.0)) R=ACOS(COSCB) SINC=SIN(R) CHID(IJ)=R*RTD C GRAZING ALTITUDE ZGRZ. AT ALTITUDES FOR WHICH IT IS ZERO OR C LESS, THE SUN IS NOT UP; THE J COEFFICIENTS SHOULD BE SET TO C ZERO, AND THE CHAPMAN FUNCTION IS NOT REQUIRED. IF(CHID(IJ).LE.90.0)THEN ZGRZ(IJ)=1000.0 ELSE ZGRZ(IJ)=SINC*(ZP+R0)-R0 END IF c type *,'ij,zgrz(ij),zp,r0,chid(ij)',ij,zgrz(ij),zp,r0,chid(ij) C CHAPMAN FUNCTION. APPROXIMATE IT WITH THE SECANT(CHI) FOR C CHI LESS THAN 70 DEGREES. IF(CHID(IJ).LT.70.0)THEN SFAC(IJ)=1.0/COSCB c print *,' ij=',ij,' coscb=',coscb,' sfac=',sfac(ij) ELSE IF(ZGRZ(IJ).GT.0.0)THEN C CHAPMAN FUNCTION SFAC FOR CHID GREATER THAN OR EQUAL TO 70 C DEGREES WHEN ZGRZ IS POSITIVE. R=SQRT(0.50*R0/HBAR) S=R*ABS(COSCB) IF(S.LE.8.0)THEN S=(D1+D2*S)/(D3+D4*S+S**2) ELSE S=D5/(D6+S) END IF R=R*SQRT(PI) SFAC(IJ)=R*S IF(CHID(IJ).GT.90.0)SFAC(IJ)=2.0*R*EXP((R0+ZBAR)*( . 1.0-SINC)/HBAR) - SFAC(IJ) ELSE SFAC(IJ)=-9999.0 END IF C ADJUST THE LENGTH OF DAY AND LENGTH OF NIGHT FOR THE C CURRENT LATITUDE. TOTAL DAYLIGHT=24*ACOS(-TAN(L)*TAN(D))/PI C FROM RUNDEL (1977), WITH L=LATITUDE, D=DECLINATION. R=SINL*SIND/(COSL*COSD) IF(R.GE.1.e0)THEN TDF=1. ELSE IF(R.LE.-1.)THEN TDF=0. ELSE TDF=ACOS(-R)*RPI END IF TNF=1.-TDF TFD(IJ)=TDF c print *,' ij,sfac, zgrz tfd=',ij,sfac(ij),zgrz(ij), c * tfd(ij) 3000 CONTINUE C print *,' sfac=',sfac C print *,' zgrz=',zgrz SAVE RETURN END