C PROGRAM oceanspec C John T.O. Kirk, Kirk Marine Optics, 2004. C MONTE CARLO SIMULATION OF ATTENUATION C OF SOLAR RADIATION (300-4000 NM) WITHIN NATURAL WATERS. C For the surface-incident radiant flux we use the spectral C distribution of solar radiation, direct and diffuse, calculated C for the specified meteorological conditions with Fortran C program solirrad.for. C The concentrations of the optically-significant components of the C water - CDOM, IPM, phytoplankton - are specified, as is also the C scattering coefficient due to particles. C The scattering phase function to be used is in the program under C the name, UPSPEC.DAT, (see under line 6, below). For each run C the phase function appropriate to the water under study, C e.g. San_Diego.dat or baham8.dat, should C be copied under the name, upspec.dat. C THE SPECTRAL DISTRIBUTION OF DOWNWARD AND OF UPWARD C IRRADIANCE IS CALCULATED FOR EVERY TENTH DEPTH INTERVAL. C THE RATE OF ENERGY ABSORPTION PER UNIT AREA IS CALCULATED FOR C EACH OF THE 100 LAYERS INTO WHICH THE TOTAL DEPTH IS DIVIDED. C The spectral distribution of the radiant flux emergent through the C surface is calculated. C IT IS ASSUMED THAT THE BOTTOM CAN REFLECT PHOTONS AND C ITS SPECTRAL REFLECTIVITY IS SPECIFIED. C THE ANGULAR DISTRIBUTION OF SCATTERING USES C EQUAL (5 DEG.) ANGULAR INTERVALS. INTEGER ZSTART,ZFINSH,BOTTOM,S,BTM real lamda,irradprop,irrad,mixlayer,ipm DOUBLE PRECISION PATH,DEPTH,BPATH,SPATH,BDEPTH DIMENSION DFLUX(165,101),UPFLUX(165,101),Z(102),SINTH(36), * SCATT(72),COSANG(36),SEAWATER(36),BENTHIC(165), * CSINTH(36),WATER(20),propheat(102),btotal(165), * ANGLE(36),SINANG(36),A(165),P(165),cdom(165),ipm(165), * TDFLUX(102),TUPFLX(102),lamda(167),HEAT(102),Enet(102), * DSPEC(165,101),UPSPEC(165,101),W(165),refl(165),deltalamda(165), * raman(165),c(165),fcdom(165),fphyto(165),exraman(165), * emraman(165),nlamda(10),Pwater(165),bw(165),Wdspec(165,101), * Wupspec(165,101),Qdspec(165,101),Qupspec(165,101), * emergent(165),emreflect(165),Qemergent(165) dimension Pquantdd(165),Pquantds(165),Pdirect(165),Psky(165), * Ptotal(165),Eu(165),sunwatts(165),escapewatts(165), * sediment(165),sedimwatts(165),wattsdd(165),reflect(165), * Wperphoton(165),reflectwatts(165),wattstotal(165),wattspec(165) dimension phyto(165),phytosp(165),phactiv(165),phinact(165), * cfcdom(165),cfphyto(165),cdomfluor(165),chlfluor(165) dimension suntheta(18),skytheta(9),azimuth(18),cospsi(18,9,18), * psi(18,9,18),skyrad(18,9,18),avskyrad(9,18),irrad(9,18), * irradprop(9,18),skyprop(9) screen=6 open(unit=70,file='Pquanta.dat',access='sequential', * form='formatted',recl=80) open(unit=80,file='ratio.dat',access='sequential', * form='formatted',recl=7) open(unit=90,file='sky.dat',access='sequential', * form='formatted',recl=72) open(unit=100,file='emflux.dat',access='sequential', * form='formatted',recl=50) open(unit=110,file='fluor.dat',access='sequential', * form='formatted',recl=50) open(unit=120,file='raman.dat',access='sequential', * form='formatted',recl=50) open(unit=130,file='spectra.dat',access='sequential', * form='formatted',recl=80) open(unit=140,file='phyto.dat',access='sequential', * form='formatted',recl=50) open(unit=150,file='wattspectra.dat',access='sequential', * form='formatted',recl=80) open(unit=160,file='Qspectra.dat',access='sequential', * form='formatted',recl=80) open(unit=170,file='water.dat',access='sequential', * form='formatted',recl=120) lamda(0)=295 lamda(166)=4100 do 4 L=1,165 read(70,2)lamda(L),Pquantdd(L),Pquantds(L),wattsdd(L), * wattstotal(L) 2 format(f6.1,X,f10.9,X,f10.9,X,f8.3,X,f8.3) 4 continue read(80,6)ratio,sumwatts 6 format(f6.3,/,f8.1) write(6,*)'ratio= ',ratio,'watts= ',sumwatts OPEN(UNIT=10,FILE='UPSPEC.DAT',ACCESS='sequential', * FORM='FORMATTED',RECL=7) OPEN(UNIT=20,FILE='OUTPUT.DAT',ACCESS='sequential', * FORM='FORMATTED',RECL=150) OPEN(UNIT=30,FILE='phytoplankton.dat',ACCESS='sequential', * FORM='FORMATTED',RECL=120) OPEN(UNIT=40,FILE='BENTHIC.DAT',ACCESS='sequential', * FORM='FORMATTED',RECL=12) OPEN(UNIT=50,FILE='ABSCO.DAT',ACCESS='sequential', * FORM='FORMATTED',RECL=80) OPEN(UNIT=60,FILE='input.dat',ACCESS='sequential', * FORM='FORMATTED',RECL=12) DO 15 I=1,72 READ(10,10)SCATT(I) 10 FORMAT(F7.4) 15 CONTINUE DO 20 I=1,36 SINTH(I)=SCATT(I) 20 CSINTH(I)=SCATT(I+36) C The values of the sine of the scattering angle, and the cumulative C distribution for scattering angle have been read in. do 30 L=21,87 READ(40,25)BENTHIC(L) 25 FORMAT(F8.5) 30 CONTINUE do 35 L=1,20 35 benthic(L)=benthic(21) do 40 L=88,165 40 benthic(L)=benthic(87) C We have entered benthic reflectivity values over the spectral range C 400-725.5 nm, based on a measured or assumed spectral reflectivity C curve in the visible/far-red. In the absence of data we have, for C simplicity assumed that the spectral reflectivity curve is flat in C in the UV and the infrared. C do 45 L=1,165 C 45 benthic(L)=0.2 C We can, if we prefer, specify a constant reflectivity at all C wavelengths. At the moment this step is commented out. Read(60,50)Bp,sunalt,DI,SOL,mixlayer,II,wind,cdom(29),ipm(29), * chl 50 format(f5.2,/,f4.1,/,f5.2,/,f10.0,/,f4.1,/,i6,/,f4.1, * /,f8.4,/,f8.4,/,f6.2) write(6,60)Bp,sunalt,DI,SOL,mixlayer,II,wind,cdom(29),ipm(29), * chl 60 format(f5.2,/,f4.1,/,f5.2,/,f10.0,/,f4.1,/,i6,/,f4.1, * /,f8.4,/,f8.4,/,f6.2) slopesquare=0.003+0.00512*wind C We have calculated the mean square slope of the water surface C from the wind speed. MM=100000000 M1=10000 IB=31415821 IR=1234567 IR=IRAN(MM,M1,IB,IR) N=0 80 N=N+1 IF(N.GT.II)GO TO 90 IR=IRAN(MM,M1,IB,IR) GO TO 80 90 CONTINUE w(1)=0.008 w(2)=0.008 w(3)=0.008 w(4)=0.008 w(5)=0.008 w(6)=0.008 w(7)=0.008 w(8)=0.008 w(9)=0.008 w(10)=0.008 w(11)=0.008 w(12)=0.008 w(13)=0.008 w(14)=0.008 w(15)=0.008 w(16)=0.008 w(17)=0.011 w(18)=0.009 w(19)=0.009 w(20)=0.008 W(21)=0.007 w(22)=0.005 W(23)=0.005 w(24)=0.004 W(25)=0.005 w(26)=0.005 W(27)=0.005 w(28)=0.005 W(29)=0.006 w(30)=0.008 W(31)=0.009 w(32)=0.010 W(33)=0.010 w(34)=0.010 W(35)=0.011 w(36)=0.011 W(37)=0.013 w(38)=0.014 W(39)=0.015 w(40)=0.017 W(41)=0.020 w(42)=0.026 W(43)=0.033 w(44)=0.040 W(45)=0.041 w(46)=0.042 W(47)=0.043 w(48)=0.045 W(49)=0.047 w(50)=0.051 W(51)=0.057 w(52)=0.060 W(53)=0.062 w(54)=0.064 W(55)=0.070 w(56)=0.077 W(57)=0.090 w(58)=0.110 W(59)=0.135 w(60)=0.167 W(61)=0.222 w(62)=0.258 W(63)=0.264 w(64)=0.268 W(65)=0.276 w(66)=0.283 W(67)=0.292 w(68)=0.301 W(69)=0.311 w(70)=0.325 W(71)=0.340 w(72)=0.371 W(73)=0.410 w(74)=0.429 W(75)=0.439 w(76)=0.448 W(77)=0.465 w(78)=0.486 W(79)=0.516 w(80)=0.559 W(81)=0.624 w(82)=0.827 W(83)=1.119 w(84)=1.489 W(85)=2.380 w(86)=2.490 W(87)=2.530 w(88)=2.540 W(89)=2.520 w(90)=2.360 W(91)=2.070 w(92)=1.95 w(93)=2.236 w(94)=2.92 w(95)=3.54 w(96)=3.97 w(97)=5.25 w(98)=7.15 w(99)=9.05 w(100)=11.9 w(101)=13.8 w(102)=17.2 w(103)=27.4 w(104)=48.3 w(105)=50.2 w(106)=45.0 w(107)=19.7 w(108)=14.8 w(109)=19.2 w(110)=29.5 w(111)=43.6 w(112)=84.0 do 105 i=113,165 105 w(i)=100.0 C THE ABSORPTION COEFFICIENTS OF PURE WATER HAVE NOW BEEN ENTERED C 300-375 nm, based on Quickenden & Irving (1980). C 380-725 nm, from Pope & Fry (1997). C 740-800 nm, from Smith & Baker (1981). C 816-1145 nm, from Palmer & Williams (1974). C We assume that at all wavelengths > 1145 nm, a(water)=100 m^-1, C i.e. water is effectively black in this region of the infrared. do 107 L=1,37 write(170,106)lamda(L),w(L),lamda(L+37),w(L+37), * lamda(L+74),w(L+74) 106 format(f6.1,10X,f7.4,10X,f6.1,10X,f7.4,10X,f6.1,10X,f7.4) 107 continue DO 110 L=1,165 cdom(L)=cdom(29)*exp(-0.014*(lamda(L)-440.0)) ipm(L)=ipm(29)*exp(-0.01*(lamda(L)-440.0)) 110 CONTINUE C We have calculated the absorption coefficients due to CDOM and to C inanimate particulate matter at all wavelengths. rad=0.0174533 BOTTOM=101 BDEPTH=100*DI mixed=INT(mixlayer/DI)+2 do 120 L=1,165 Pdirect(L)=Pquantdd(L)*sol Psky(L)=Pquantds(L)*ratio*sol Ptotal(L)=Pdirect(L)+Psky(L) deltalamda(L)=(lamda(L+1)-lamda(L-1))/2000.0 Wperphoton(L)=wattsdd(L)/Pdirect(L) 120 continue do 130 L=151,153 130 Wperphoton(L)=Wperphoton(150) C To avoid numerical problems later we ensure here that in the C 2.6, 2.7 and 2.8 micron wavebands, where water vapour absorbs very C strongly, so that there are no watts or photons with which to C calculate 'watts per photon', there is a finite value of C Wperphoton(L) for the computer to work with. C For the irradiance due to both the direct solar beam and the sky C the photons have been distributed amongst the 165 wavebands, in C accordance with the distributions calculated for the specified C meteorological conditions, with program solirrad. C **** We now calculate the angular distribution of irradiance **** C at sea level, originating from the sky, using C the equation of Harrison & Coombes (1988). do 140 k=1,18 suntheta(k)=(k-1)*5.0*rad 140 continue do 160 j=1,9 160 skytheta(j)=(10*j-5.0)*rad do 180 i=1,18 180 azimuth(i)=(10*i-5.0)*rad do 260 k=1,18 do 240 j=1,9 sumskyrad=0.0 do 220 i=1,18 cospsi(i,j,k)=cos(skytheta(j))*cos(suntheta(k)) * +sin(skytheta(j))*sin(suntheta(k))*cos(azimuth(i)) psi(i,j,k)=acos(cospsi(i,j,k)) skyrad(i,j,k)=(1.63+53.7*exp(-5.94*psi(i,j,k)) * +2.04*(cospsi(i,j,k)**2)*cos(suntheta(k))) * *(1.0-exp(-1.9/cos(skytheta(j)))) * *(1.0-exp(-0.53/cos(suntheta(k)))) 220 sumskyrad=sumskyrad+skyrad(i,j,k) avskyrad(j,k)=sumskyrad/18.0 240 continue 260 continue do 320 k=1,18 sumirrad=0.0 do 280 j=1,9 irrad(j,k)=avskyrad(j,k)*cos(skytheta(j))*sin(skytheta(j)) sumirrad=sumirrad+irrad(j,k) 280 continue do 300 j=1,9 irradprop(j,k)=irrad(j,k)/sumirrad 300 continue n=int((90.0-sunalt)/5.0) 320 continue do 350 j=1,9 skyprop(j)=irradprop(j,n) write(6,*)skyprop(j),'n= ',n write(90,340)(skytheta(j)/c),skyprop(j) 340 format(f5.1,'sky degrees ',f6.4) 350 continue C We have calculated a set of numbers which are proportional to C the irradiance values in 9 angular intervals from 0-10 deg to C 80-90 deg, and then calculated the proportion of the total irradiance C which is in each of these angular intervals. C This calculation has been carried out for 18 different solar zenith C angles from 0 deg to 85 deg. C From amongst these we have separated out, as the array skyprop(j), C the angular distribution of sky irradiance corresponding to the C particular solar altitude we have chosen. C ** End of sky irradiance calculation ** C **** Spectral data **** do 360 i=1,35 360 cfcdom(i)=0.0 cfcdom(36)=0.0055 CFCDOM(37)=0.0111 cfcdom(38)=0.0256 CFCDOM(39)=0.0400 cfcdom(40)=0.0657 CFCDOM(41)=0.0913 cfcdom(42)=0.1306 CFCDOM(43)=0.1699 cfcdom(44)=0.2197 CFCDOM(45)=0.2694 cfcdom(46)=0.3224 CFCDOM(47)=0.3753 cfcdom(48)=0.4271 CFCDOM(49)=0.4788 cfcdom(50)=0.5250 CFCDOM(51)=0.5711 cfcdom(52)=0.6116 CFCDOM(53)=0.6521 cfcdom(54)=0.6862 CFCDOM(55)=0.7203 cfcdom(56)=0.7488 CFCDOM(57)=0.7773 cfcdom(58)=0.8010 CFCDOM(59)=0.8246 cfcdom(60)=0.8439 CFCDOM(61)=0.8631 cfcdom(62)=0.8792 CFCDOM(63)=0.8952 cfcdom(64)=0.9077 CFCDOM(65)=0.9201 cfcdom(66)=0.9300 CFCDOM(67)=0.9400 cfcdom(68)=0.9476 CFCDOM(69)=0.9552 cfcdom(70)=0.9614 CFCDOM(71)=0.9676 cfcdom(72)=0.9728 CFCDOM(73)=0.9780 cfcdom(74)=0.9818 CFCDOM(75)=0.9856 cfcdom(76)=0.9883 CFCDOM(77)=0.9909 cfcdom(78)=0.9925 CFCDOM(79)=0.9941 cfcdom(80)=0.9953 CFCDOM(81)=0.9965 cfcdom(82)=0.9981 CFCDOM(83)=0.9992 cfcdom(84)=0.9996 CFCDOM(85)=1.0000 do 370 i=86,165 370 cfcdom(i)=0.0 do 380 i=1,70 380 cfphyto(i)=0.0 cfphyto(71)=0.0 cfphyto(72)=0.0100 CFPHYTO(73)=0.0201 cfphyto(74)=0.0604 CFPHYTO(75)=0.1007 cfphyto(76)=0.2417 CFPHYTO(77)=0.3826 cfphyto(78)=0.5570 CFPHYTO(79)=0.7315 cfphyto(80)=0.7836 CFPHYTO(81)=0.8356 cfphyto(82)=0.9128 CFPHYTO(83)=0.9644 cfphyto(84)=0.9776 CFPHYTO(85)=1.0000 do 390 i=86,165 390 cfphyto(i)=0.0 C We have now entered the cumulative frequency distributions C of CDOM and chlorophyll fluorescence as a function of wavelength. do 400 i=1,8 400 phytosp(i)=0.0150 phytosp(1)=0.02394 phytosp(2)=0.02465 phytosp(3)=0.02535 phytosp(4)=0.02606 phytosp(5)=0.02676 phytosp(6)=0.02711 phytosp(7)=0.02746 phytosp(8)=0.02782 phytosp(9)=0.02817 phytosp(10)=0.02641 phytosp(11)=0.02465 phytosp(12)=0.02324 phytosp(13)=0.02183 phytosp(14)=0.02200 phytosp(15)=0.02218 phytosp(16)=0.02306 phytosp(17)=0.02394 phytosp(18)=0.02500 phytosp(19)=0.02606 phytosp(20)=0.02712 PHYTOSP(21)=0.02817 phytosp(22)=0.03046 PHYTOSP(23)=0.03275 phytosp(24)=0.03451 PHYTOSP(25)=0.03627 phytosp(26)=0.03750 PHYTOSP(27)=0.03873 phytosp(28)=0.04049 PHYTOSP(29)=0.04225 phytosp(30)=0.04172 PHYTOSP(31)=0.04120 phytosp(32)=0.04138 PHYTOSP(33)=0.04155 phytosp(34)=0.04120 PHYTOSP(35)=0.04085 phytosp(36)=0.03908 PHYTOSP(37)=0.03732 phytosp(38)=0.03521 PHYTOSP(39)=0.03310 phytosp(40)=0.03010 PHYTOSP(41)=0.02711 phytosp(42)=0.02464 PHYTOSP(43)=0.02218 phytosp(44)=0.02042 PHYTOSP(45)=0.01866 phytosp(46)=0.01690 PHYTOSP(47)=0.01514 phytosp(48)=0.01345 PHYTOSP(49)=0.01176 phytosp(50)=0.01046 PHYTOSP(51)=0.00915 phytosp(52)=0.00782 PHYTOSP(53)=0.00648 phytosp(54)=0.00612 PHYTOSP(55)=0.00577 phytosp(56)=0.00588 PHYTOSP(57)=0.00599 phytosp(58)=0.00634 PHYTOSP(59)=0.00669 phytosp(60)=0.00634 PHYTOSP(61)=0.00599 phytosp(62)=0.00596 PHYTOSP(63)=0.00592 phytosp(64)=0.00613 PHYTOSP(65)=0.00634 phytosp(66)=0.00669 PHYTOSP(67)=0.00704 phytosp(68)=0.00740 PHYTOSP(69)=0.00775 phytosp(70)=0.00757 PHYTOSP(71)=0.00739 phytosp(72)=0.00862 PHYTOSP(73)=0.00986 phytosp(74)=0.01398 PHYTOSP(75)=0.01810 phytosp(76)=0.02000 PHYTOSP(77)=0.01761 phytosp(78)=0.01169 PHYTOSP(79)=0.00577 phytosp(80)=0.00366 PHYTOSP(81)=0.00155 phytosp(82)=0.00056 PHYTOSP(83)=0.00035 phytosp(84)=0.00017 PHYTOSP(85)=0.00000 do 410 i=86,165 410 PHYTOSP(i)=0.00000 C The specific absorption coefficients of phytoplankton C (m^2 mg chl a^-1) have now been entered do 420 i=1,19 420 phactiv(i)=0.0 phactiv(20)=0.00106 PHACTIV(21)=0.00216 phactiv(22)=0.00289 PHACTIV(23)=0.00366 phactiv(24)=0.00578 PHACTIV(25)=0.00789 phactiv(26)=0.01077 PHACTIV(27)=0.01366 phactiv(28)=0.01543 PHACTIV(29)=0.01718 phactiv(30)=0.01986 PHACTIV(31)=0.02254 phactiv(32)=0.02394 PHACTIV(33)=0.02535 phactiv(34)=0.02422 PHACTIV(35)=0.02310 phactiv(36)=0.02000 PHACTIV(37)=0.01690 phactiv(38)=0.01578 PHACTIV(39)=0.01465 phactiv(40)=0.01508 PHACTIV(41)=0.01550 phactiv(42)=0.01557 PHACTIV(43)=0.01564 phactiv(44)=0.01571 PHACTIV(45)=0.01578 phactiv(46)=0.01529 PHACTIV(47)=0.01479 phactiv(48)=0.01345 PHACTIV(49)=0.01211 phactiv(50)=0.00944 PHACTIV(51)=0.00676 phactiv(52)=0.00408 PHACTIV(53)=0.00141 do 430 i=54,165 430 phactiv(i)=0.0 C The specific absorption coefficients (400-560 nm)of the pigment C system which excites phytoplankton fluorescence have now C been entered. The shape of the spectrum is based on the C the excitation spectrum for diatoms/dinoflagellates/cocco C -lithophores chlorophyll fluorescence given by Yentsch, 1994. do 440 i=1,165 440 phinact(i)=phytosp(i)-phactiv(i) C The specific absorption coefficients of the inactive part of C phytoplankton absorption in the range 400-560 nm have now C been calculated. sumphyto=0.0 sumphactiv=0.0 do 450 L=1,85 sumphyto=sumphyto+phytosp(L) sumphactiv=sumphactiv+phactiv(L) write(140,445)lamda(L),phytosp(L),phactiv(L) 445 format(f6.1,2x,f8.7,2x,f8.7) 450 continue DO 460 L=1,165 PHYTO(L)=CHL*PHYTOSP(L) PHACTIV(L)=CHL*PHACTIV(L) 460 A(L)=W(L)+CDOM(L)+ipm(L)+PHYTO(L) C We have calculated the absorption coefficients C due to phytoplankton at 165 wavelengths and then calculated C the total absorption coefficients at these wavelengths. do 470 L=1,165 bw(L)=0.0045*((450.0/lamda(L))**4.32) btotal(L)=bp*660.0/lamda(L)+bw(L) raman(L)=0.00025*((490.0/lamda(L))**4.77) 470 c(L)=a(L)+btotal(L)+raman(L) C We are assuming that the input (particle) scattering coefficient C is the value at 660 nm, and that bp varies inversely with wavelength. C At each wavelength we add the appropriate value of the scattering C coefficient of water. C We have also calculated the Raman scattering cofficient and beam C attenuation coefficient as a function of wavelength. do 490 L=1,165 WRITE(50,480)A(L),btotal(L),c(L),raman(L) 480 FORMAT(F7.4,2x,f7.4,2x,f7.4,2x,f8.6) 490 continue do 530 L=1,29 write(30,520)lamda(L),phytosp(L),lamda(L+29),phytosp(L+29), * lamda(L+58),phytosp(L+58) 520 format(f6.1,10X,f7.5,10X,f6.1,10X,f7.5,10X,f6.1,10X,f7.5) 530 continue C *** end of spectral data section *** do 1000 L=1,165 1000 Eu(L)=0.0 DO 1020 I=1,101 TDFLUX(I)=0.0 1020 TUPFLX(I)=0.0 1030 ESCAPE=0.0 PN=0.0 photons=0.0 reflect=0.0 sedphotons=0.0 n=0 m=0 1040 m=m+1 if(m.eq.3)go to 5200 L=0 1050 L=L+1 IF(L.GT.165)GO TO 1040 if(m.eq.1)P(L)=Pdirect(L) if(m.eq.2)P(L)=Psky(L) PN=0.0 1100 PN=PN+1.0 photons=photons+1.0 C 'photons' counts the total number of photons in the direct solar beam C and from the sky, incident on the surface. lam=L C 'lam' is the working wavelength of each photon. In the course of the C life history of the photon this can change as the result of C fluorescence or Raman scattering, and the subsequent fate of the C photon will be influenced by this change. When that photon is finally C extinguished, the next photon will once again start off with lam=L. C A photon encounters the surface. atheta=(90.0-sunalt)*rad 1120 IF(PN.GT.P(L))GO TO 1050 if(m.eq.1)go to 1200 C This photon is in the direct solar beam. go to 1160 C This photon is part of the sky irradiance. 1140 n=n+1 skyphoton=0.0 1160 skyphoton=skyphoton+1.0 if(skyphoton.gt.skyprop(n)*P(L))go to 1140 atheta=skytheta(n)*rad 1200 IR=IRAN(MM,M1,IB,IR) R=REAL(IR)/100000000.0 slopeangle=atan(sqrt(-2.0*slopesquare*log(1.0-R))) IR=IRAN(MM,M1,IB,IR) R=REAL(IR)/100000000.0 if(R.gt.0.5)slopeangle=-slopeangle atheta=atheta-slopeangle wtheta=asin(sin(atheta/1.333)) SUM=ATHETA+WTHETA DIFF=ATHETA-WTHETA X=SIN(DIFF)/SIN(SUM) Y=TAN(DIFF)/TAN(SUM) RSURF=0.5*(X*X+Y*Y) IR=IRAN(MM,M1,IB,IR) R=REAL(IR)/100000000.0 if(R.gt.RSURF)go to 1300 reflect(L)=reflect(L)+1.0 go to 1100 1300 wtheta=wtheta+slopeangle C We have selected the slope angle of the surface, and have determined C whether it is towards (slopeangle +ve) or C away from (slopeangle -ve) the photon. By subtracting the slope angle C from the zenith angle of the photon we obtain the angle of the C photon to the normal to the surface. C Application of Snell's Law gives us the angle (wtheta) of the photon C within the water to the normal to the surface. C We use the Fresnel equation to calculate the probability of reflection C and with a random number determine whether reflection takes place: if C it does we note the reflection event and start a new photon. C If the photon passes through the surface we determine its true nadir C angle by correcting for the surface slope angle. Pwater(L)=Pwater(L)+1.0 C We count the number of photons that enter the water in C each waveband. 1400 SINALF=SIN(1.5707963-wtheta) C alpha is the angle of the photon trajectory to the horizontal. 1500 DEPTH=0.0 S=0 BTM=0 NSCAT=0 1600 NSCAT=NSCAT+1 IR=IRAN(MM,M1,IB,IR) R=REAL(IR)/100000000.0 IF((1.0-R).LT.0.00001)R=0.99999 PATH=-ALOG(1.0-R)/(A(lam)+btotal(lam)) C HAVE DETERMINED THE PATHLENGTH OF THE PHOTON C BEFORE ITS NEXT INTERACTION IF(NSCAT.EQ.1) GO TO 2100 C IF PHOTON HAS BEEN SCATTERED DETERMINE NEW TRAJECTORY 1620 IR=IRAN(MM,M1,IB,IR) R=REAL(IR)/100000000.0 if(R.gt.bw(lam)/btotal(lam))go to 1700 C We have determined whether the photon was scattered by water or by C particles. In the case of water we use a calculated cumulative C distribution. In the case of particles we use the distribution C obtained from the measured volume scattering function for C San Diego Harbour. waterscat=waterscat+1.0 IR=IRAN(MM,M1,IB,IR) R=REAL(IR)/100000000.0 S=-4.594*(R-0.5) T=sqrt(21.105*(R-0.5)**2+1.71794) STsum=S+T STdiff=S-T if(STsum.ge.0.0)go to 1640 STsum=-STsum cubicA=-(STsum**0.33333) go to 1660 1640 cubicA=STsum**0.33333 1660 continue if(STdiff.ge.0.0)go to 1670 STdiff=-STdiff cubicB=-(STdiff**0.33333) go to 1680 1670 cubicB=STdiff**0.33333 1680 continue costheta=cubicA+cubicB go to 2020 1700 IR=IRAN(MM,M1,IB,IR) R=REAL(IR)/100000000.0 KA=0 1800 KA=KA+1 IF(R.GT.CSINTH(KA))GO TO 1800 SINTHETA=SINTH(KA) COSTHETA=SQRT(1.0-SINTHETA*SINTHETA) IF(KA.GT.18)COSTHETA=-COSTHETA 2020 COSALF=SQRT(1.0-SINALF*SINALF) IR=IRAN(MM,M1,IB,IR) R=REAL(IR)/100000000.0 PHI=R*6.2832 SINALF=COSTHETA*SINALF+SINTHETA*COSALF*COS(PHI) IF(ABS(SINALF).LT.0.0001)SINALF=SINALF+0.001 C We must make sure that the computer never rounds SINALF C off to zero since we sometimes divide by SINALF. C SIN ALPHA HAS BEEN CALCULATED FOR THE NEW TRAJECTORY IF(SINALF.LE.0.0)GO TO 2700 C POSITIVE ALPHA MEANS PHOTON TRAVELLING DOWNWARDS 2100 ZSTART=DEPTH/DI+2 ZFINSH=(DEPTH+PATH*SINALF)/DI+1 BPATH=(BDEPTH-DEPTH)/SINALF C BPATH is the length of the path to the bottom from C the point at which the current photon trajectory C started (regardless of whether this particular C trajectory actually extends as far as the bottom). DEPTH=DEPTH+PATH*SINALF 2150 IF(ZFINSH.LT.ZSTART)GO TO 4100 C This trajectory starts and ends within the same C depth interval. IF(ZFINSH.LT.BOTTOM)GO TO 2400 C This trajectory ends at some point in another depth C interval above the bottom. 2200 DO 2280 I=ZSTART,BOTTOM 2280 DFLUX(lam,I)=DFLUX(lam,I)+1.0 C PHOTON HAS REACHED BOTTOM 2300 IR=IRAN(MM,M1,IB,IR) R=REAL(IR)/100000000.0 if(R.gt.benthic(lam))sediment(lam)=sediment(lam)+1.0 if(R.gt.benthic(lam))sedphotons=sedphotons+1.0 IF(R.GT.BENTHIC(lam))GO TO 1100 C IF THAT PHOTON WAS ABSORBED BY THE BOTTOM A NEW ONE IS REQUIRED. C IF NOT ABSORBED THE PHOTON IS REFLECTED, AND ITS SUBSEQUENT C FATE IS FOLLOWED. IF(S.EQ.1)BPATH=BDEPTH/SINALF C THIS IS FOR PHOTONS REFLECTED ALL THE WAY FROM THE SURFACE C TO THE BOTTOM S=0 C THIS CANCELS THE FLAG INDICATING REFLECTION FROM THE SURFACE PATH=PATH-BPATH C THIS CALCULATES THE REMAINDER OF THE PATH THAT THE PHOTON HAS C TO COMPLETE AFTER REFLECTION FROM THE BOTTOM 2350 IR=IRAN(MM,M1,IB,IR) R=REAL(IR)/100000000.0 SINALF=-SQRT(R) IF(ABS(SINALF).LT.0.0001)SINALF=SINALF-0.001 C FOR PHOTONS REFLECTED FROM A LAMBERTIAN SURFACE, THE PROPORTION OF C TOTAL FLUX WHICH OCCURS UP TO ALPHA IS (SIN ALPHA) SQUARED. C SIN ALPHA IS NEGATIVE SINCE PHOTON IS TRAVELLING UPWARDS. ZSTART=BOTTOM ZFINSH=(BDEPTH+PATH*SINALF)/DI+2 DEPTH=ABS(BDEPTH+PATH*SINALF) IF(ZFINSH.EQ.BOTTOM)GO TO 4100 C THE REFLECTED PHOTON'S TRAJECTORY ENDS WITHIN THE LOWEST DEPTH INTERVAL IF(1.LT.ZFINSH)GO TO 3900 C THE REFLECTED PHOTON'S TRAJECTORY ENDS AT SOME POINT ABOVE THE C LOWEST DEPTH INTERVAL BUT BELOW THE SURFACE. C IF NEITHER OF THE ABOVE CONDITIONS HOLDS (I.E. IF ZFINSH IS 1,0 OR NEGATIVE) C THE REFLECTED PHOTON MUST REACH THE SURFACE. BTM=1 C THIS IS A FLAG MARKING THE FACT THAT THIS SECTION OF THE PHOTON C TRAJECTORY ORIGINATES AT THE BOTTOM. GO TO 2900 2400 DO 2480 I=ZSTART,ZFINSH 2480 DFLUX(lam,I)=DFLUX(lam,I)+1.0 2600 GO TO 4100 C DETERMINE WHETHER THIS PHOTON IS TO BE SCATTERED AGAIN OR C TERMINATED C NEGATIVE ALPHA MEANS PHOTON TRAVELLING UPWARDS 2700 ZSTART=DEPTH/DI+1 ZFINSH=(DEPTH+PATH*SINALF)/DI+2 SPATH=ABS(DEPTH/SINALF) C SPATH is the length of the path to the surface from C the point at which the current photon trajectory c started (regardless of whether this particular C trajectory actually extends as far as the surface). DEPTH=ABS(DEPTH+PATH*SINALF) IF(ZFINSH.GT.ZSTART)GO TO 4100 C This trajectory starts and ends within the same C depth interval. IF(1.LT.ZFINSH)GO TO 3900 C This trajectory finishes at some point below the surface C If neither of the above conditions holds, the upward C scattered photon must reach the surface. 2900 DO 3000 I=1,ZSTART 3000 UPFLUX(lam,I)=UPFLUX(lam,I)+1.0 WSINTH=SQRT(1.0-SINALF*SINALF) wtheta=asin(wsinth) wtheta=wtheta+slopeangle C For simplicity we are assuming that the surface slope angle C encountered from below by this upward-travelling photon is the same C as that encountered by the photon when it first entered the water. C While for individual photons this will generally not be the case, C the statistical distribution of slope angles will be correct. if(wtheta.GT.(rad*48.6))go to 3200 C ALL PHOTONS INCIDENT ON SURFACE FROM BELOW AT C THETA > 48.6 DEG. ARE REFLECTED DOWNWARDS. ATHETA=ASIN(1.333*sin(wtheta)) SUM=ATHETA+WTHETA IF(SIN(SUM).EQ.0.0)GO TO 3100 C This is the special case of a photon being incident C on the surface at theta = 0.0 degrees. DIFF=WTHETA-ATHETA X=SIN(DIFF)/SIN(SUM) Y=TAN(DIFF)/TAN(SUM) RRSURF=0.5*(X*X+Y*Y) GO TO 3150 3100 RRSURF=0.02 C HAVE CALCULATED THE PROBABILITY OF DOWNWARD REFLECTION C FOR THE PHOTON WHICH HAS JUST REACHED THE SURFACE C FROM BELOW AT THETA < 48.6 DEG. 3150 IR=IRAN(MM,M1,IB,IR) R=REAL(IR)/100000000.0 IF(R.GT.RRSURF)GO TO 3700 C HAVE DETERMINED WHETHER THIS PHOTON WAS REFLECTED C DOWNWARDS OR ESCAPED THROUGH THE SURFACE. IF THE C FORMER, ITS FURTHER FATE IS FOLLOWED. 3200 IF(BTM.EQ.1)SPATH=ABS(BDEPTH/SINALF) C THIS IS FOR PHOTONS WHICH HAVE BEEN REFLECTED ALL THE WAY FROM C THE BOTTOM TO THE SURFACE PATH=PATH-SPATH C HAVE CALCULATED THE REMAINING PATHLENGTH THE PHOTON MUST TRAVERSE C AFTER INTERNAL REFLECTION AT THE SURFACE. BTM=0 C THIS CANCELS THE FLAG INDICATING REFLECTION FROM THE BOTTOM. SINALF=-SINALF C THE SIGN OF ALPHA MUST BE REVERSED SINCE THE PHOTON IS NOW C TRAVELLING DOWNWARDS AGAIN. ZFINSH=DEPTH/DI+1 IF(ZFINSH.LT.BOTTOM)GO TO 3500 C IF THE ABOVE CONDITION IS NOT FULFILLED, THE PHOTON MUST BE C REFLECTED ALL THE WAY TO THE BOTTOM. S=1 C THIS IS A FLAG MARKING THE FACT THAT THIS SECTION OF THE PHOTON C TRAJECTORY ORIGINATES IN REFLECTION AT THE SURFACE. 3300 DO 3380 I=1,BOTTOM 3380 DFLUX(lam,I)=DFLUX(lam,I)+1.0 GO TO 2300 C THIS PHOTON HAS REACHED THE BOTTOM AND SO WHETHER OR NOT C IT IS REFLECTED MUST NOW BE DETERMINED. 3500 DO 3580 I=1,ZFINSH 3580 DFLUX(lam,I)=DFLUX(lam,I)+1.0 GO TO 4100 C IS PHOTON TO BE SCATTERED AGAIN OR TERMINATED? C PHOTON NOT REFLECTED AT SURFACE 3700 Eu(lam)=Eu(lam)+1.0 ESCAPE=ESCAPE+1.0 GO TO 1100 C NEW PHOTON REQUIRED 3900 DO 4000 I=ZFINSH,ZSTART 4000 UPFLUX(lam,I)=UPFLUX(lam,I)+1.0 4100 IR=IRAN(MM,M1,IB,IR) R=REAL(IR)/100000000.0 IF(R.LE.((A(lam)+RAMAN(lam))/C(lam)))GO TO 4900 GO TO 1600 C HAVE DETERMINED WHETHER THAT PHOTON WAS ABSORBED OR IS TO C BE SCATTERED C **** Fluorescence and Raman emission **** C If the photon has undergone elastic scattering, control is C transferred to that part of the program where its new pathlength C and direction are calculated. 4900 IF(R.GT.(A(lam)/C(lam)))TOTRAMAN=TOTRAMAN+1.0 IF(R.GT.(A(lam)/C(lam)))GO TO 4940 C If the photon undergoes Raman scattering we transfer control C to that part of the program where its new wavelength C and direction are calculated. C Also we count the number of Raman scattering events. C Having dealt with scattering we now consider absorption. C We determine whether absorption is by water, CDOM or phytoplankton. IR=IRAN(MM,M1,IB,IR) R=REAL(IR)/100000000.0 IF(R.LE.(W(lam)/A(lam)))GO TO 4904 4904 IF(R.LE.(W(lam)/A(lam)))GO TO 1100 C The photon was absorbed by water, so a new photon is required. IF(R.GT.(W(lam)+CDOM(lam))/A(lam))PHYTEST=PHYTEST+1.0 IF(R.GT.(W(lam)+CDOM(lam))/A(lam))GO TO 4920 C Alternatively the photon was absorbed by phytoplankton C and so control is transferred to the appropriate part of C the program. C The photon we are presently left with has C been absorbed by CDOM, and before extinguishing it we must C determine if it is re-emitted as fluorescence. We assume a C fluorescence yield of 0.01 (Vodacek et al., L & O 42(1997)674). C For simplicity we shall assume that every 100th photon C absorbed by CDOM is re-emitted as fluorescence. NCDOM=NCDOM+1 IF(NCDOM.LT.100)GO TO 1100 FCDOM(lam)=FCDOM(lam)+1.0 TOTFCDOM=TOTFCDOM+1.0 NCDOM=0 IFLAG=1 C All photons absorbed by CDOM, prior to the 100th, we assume to C be extinguished, and so a new photon is required. In the case C of the 100th photon we now determine the wavelength of re-emission. C We set the variable IFLAG=1 as a flag to tell us that the photon C currently being followed was emitted as fluorescence. IR=IRAN(MM,M1,IB,IR) R=REAL(IR)/100000000.0 K=0 4910 K=K+1 IF(R.GT.CFCDOM(K))GO TO 4910 lam=K cdomfluor(lam)=cdomfluor(lam)+1.0 GO TO 4998 C Using R and the cumulative frequency distribution for CDOM C fluorescence versus wavelength, we have selected the C wavelength (lam) of emission. C We now deal with a photon which has been absorbed by C phytoplankton. We must determine if it is extinguished or C re-emitted as fluorescence. We assume that phytoplankton C have a 1% fluorescence yield overall, but that their active C pigment system contributes 27% to total phyto. absorption C and has a fluorescence yield of 3.7%. For simplicity we C shall assume that every 27th photon absorbed by the active C pigment system is re-emitted as fluorescence. First we C determine whether absorption is by active or inactive C phytoplankton pigments. 4920 IR=IRAN(MM,M1,IB,IR) R=REAL(IR)/100000000.0 IF(R.GT.(PHACTIV(lam)/PHYTO(lam)))GO TO 1100 C The photon was absorbed by inactive phytoplankton pigment C and so a new one is required. NPHYTO=NPHYTO+1 IF(NPHYTO.LT.27)GO TO 1100 FPHYTO(lam)=FPHYTO(lam)+1.0 TOTFPHYTO=TOTFPHYTO+1.0 NPHYTO=0 IFLAG=1 C All photons absorbed by phyto. active pigment, prior to C the 27th, we assume to be extinguished, and so a new C photon is required. In the case of the 27th photon we C now determine the wavelength of re-emission. IR=IRAN(MM,M1,IB,IR) R=REAL(IR)/100000000.0 K=0 4930 K=K+1 IF(R.GT.CFPHYTO(K))GO TO 4930 lam=K chlfluor(lam)=chlfluor(lam)+1.0 GO TO 4998 C Using R and the cumulative frequency distribution for C phytoplankton fluorescence versus wavelength, we have C selected the wavelength (lam) of emission. C **** FOR THE PHOTON WHICH UNDERGOES RAMAN SCATTERING WE **** C **** CALCULATE THE WAVELENGTH OF EMISSION **** 4940 IF(lam.GT.61)GO TO 1100 C For photons with wavelength > 600 nm, the wavelength after C Raman scattering is > 750 nm, and will be quickly absorbed. C A new photon is therefore required. EXRAMAN(lam)=EXRAMAN(lam)+1.0 LAMDA2=lam IF(NLAMDA(1).EQ.0)GO TO 4942 GO TO 4944 4942 IF(LAMDA2.EQ.1)lam=33 NLAMDA(1)=1 GO TO 4946 4944 IF(LAMDA2.EQ.1)LAM=35 NLAMDA(1)=0 4946 IF(LAMDA2.EQ.2)LAM=37 IF(LAMDA2.EQ.3)LAM=39 IF(NLAMDA(2).EQ.0)GO TO 4948 GO TO 4950 4948 IF(LAMDA2.EQ.4)LAM=41 NLAMDA(2)=1 GO TO 4952 4950 IF(LAMDA2.EQ.4)LAM=43 NLAMDA(2)=0 4952 IF(LAMDA2.EQ.5)LAM=45 IF(NLAMDA(3).EQ.0)GO TO 4954 GO TO 4956 4954 IF(LAMDA2.EQ.6)LAM=47 NLAMDA(3)=1 GO TO 4958 4956 IF(LAMDA2.EQ.6)LAM=49 NLAMDA(3)=0 4958 IF(LAMDA2.EQ.7)LAM=51 IF(LAMDA2.EQ.8)LAM=53 IF(NLAMDA(4).EQ.0)GO TO 4960 GO TO 4962 4960 IF(LAMDA2.EQ.9)LAM=55 NLAMDA(4)=1 GO TO 4964 4962 IF(LAMDA2.EQ.9)LAM=57 NLAMDA(4)=0 4964 IF(LAMDA2.EQ.10)LAM=59 IF(NLAMDA(5).EQ.0)GO TO 4966 GO TO 4968 4966 IF(LAMDA2.EQ.11)LAM=61 NLAMDA(5)=1 GO TO 4970 4968 IF(LAMDA2.EQ.11)LAM=63 NLAMDA(5)=0 4970 IF(LAMDA2.EQ.12)LAM=65 IF(NLAMDA(6).EQ.0)GO TO 4972 GO TO 4974 4972 IF(LAMDA2.EQ.13)LAM=67 NLAMDA(6)=1 GO TO 4976 4974 IF(LAMDA2.EQ.13)LAM=69 NLAMDA(6)=0 4976 IF(LAMDA2.EQ.14)LAM=71 IF(NLAMDA(7).EQ.0)GO TO 4978 GO TO 4980 4978 IF(LAMDA2.EQ.15)LAM=73 NLAMDA(7)=1 GO TO 4982 4980 IF(LAMDA2.EQ.15)LAM=75 NLAMDA(7)=0 4982 IF(LAMDA2.EQ.16)LAM=77 IF(NLAMDA(8).EQ.0)GO TO 4984 GO TO 4986 4984 IF(LAMDA2.EQ.17)LAM=79 NLAMDA(8)=1 GO TO 4988 4986 IF(LAMDA2.EQ.17)LAM=81 NLAMDA(8)=0 4988 IF(LAMDA2.EQ.18)LAM=83 IF(NLAMDA(9).EQ.0)GO TO 4990 GO TO 4992 4990 IF(LAMDA2.EQ.19)LAM=85 NLAMDA(9)=1 GO TO 4994 4992 IF(LAMDA2.EQ.19)LAM=87 NLAMDA(9)=0 4994 IF(LAMDA2.EQ.20)LAM=89 IF(LAMDA2.EQ.21)LAM=91 EMRAMAN(lam)=EMRAMAN(lam)+1.0 C We have determined the wavelength of the Raman-emitted photon, C assuming a frequency diminution of 3400 cm^-1 (Smith & Marshall). C Lamda(R)=Lamda(initial)/(1-0.00034*Lamda(initial)). IFLAG=1 C We set the variable IFLAG=1 as a flag to tell us that the photon C currently being followed was emitted by Raman scattering (or C fluorescence). GO TO 4998 C We make the approximate assumption that Raman scattering has C (like fluorescence) an C isotropic angular distribution, and so transfer control to that C part of the program where the new direction and pathlength of C fluorescence-derived photons is calculated. C **** FOR THE PHOTON EMITTED IN FLUORESCENCE WE CALCULATE **** C **** THE NEW DIRECTION RELATIVE TO THE ORIGINAL DIRECTION **** C **** AND THE NEW PATHLENGTH **** 4998 IR=IRAN(MM,M1,IB,IR) R=REAL(IR)/100000000.0 PHI=R*6.2832 SINPHI=SIN(PHI) COSPHI=COS(PHI) C We assume that for the new trajectory all angles of rotation C around the original direction are equally likely. IR=IRAN(MM,M1,IB,IR) R=REAL(IR)/100000000.0 costheta=1.0-2*R sintheta=sqrt(1.0-costheta**2) C We also assume that for a photon emitted in fluorescence, all C elements of solid angle are equally likely. It can be shown that C for a point source the proportion of the total surrounding spherical C surface between theta=0 and theta=theta' which is irradiated is C (1-cos theta')/2, and so we set the random number for the C cumulative distribution for cos theta equal to this. IR=IRAN(MM,M1,IB,IR) R=REAL(IR)/100000000.0 PATH=-ALOG(1.0-R)/C(lam) GO TO 2020 C Have determined the pathlength of the photon before its next C interaction, taking account - through c(lam) - of the C change in wavelength. C For further calculation of the path followed by the photon C within the cavity we transfer control to that part of the C program which deals with photons which have been scattered. C ************************* C *** End of fluorescence and Raman section *** 5200 WRITE(20,5210)SOL,ESCAPE 5210 FORMAT(10X,'TOTAL PHOTONS IN',25X,F10.0,/, * 10X,'TOTAL PHOTONS OUT',24X,F10.0,/,5X, * 'Spectral distribution of downward and upward irradiance') c write(6,*)'totraman',totraman,'phytest',phytest, c * 'totfcdom',totfcdom,'totfphyhto',totfphyto C THE SPECTRAL DISTRIBUTION OF DOWNWARD AND OF UPWARD C IRRADIANCE IS NOW CALCULATED AT A SERIES OF DEPTHS. DO 5220 J=1,11 5220 Z(J)=(J-1)*DI*10 DO 5250 L=1,165 5250 DFLUX(L,1)=DFLUX(L,1)+Pwater(L) c upspec(L,2)=upflux(L,2)/(1000.0*deltalamda(L)) c 5250 dspec(L,2)=dflux(L,2)/(1000.0*deltalamda(L)) do 5260 i=1,11 k=10*(i-1)+1 do 5260 L=1,165 dspec(L,k)=dflux(L,k)/(1000.0*deltalamda(L)) 5260 upspec(L,k)=upflux(L,k)/(1000.0*deltalamda(L)) C We have divided the downward and upward irradiance values in each C waveband by the width of the waveband (in nanometres) to give the C spectral irradiance, i.e. the spectral concentration of radiation C ([Monte Carlo] quanta m^-2 s^-1 nm^-1) do 5300 L=1,165 wattspec(L)=wattstotal(L)/deltalamda(L) C We have divided the total surface-incident irradiance in each C waveband by the width of the waveband to give the spectral C irradiance distribution of the incident flux. emergent(L)=Eu(L)*Wperphoton(L)/deltalamda(L) Qemergent(L)=5.03*emergent(L)*lamda(L)/1000.0 emreflect(L)=emergent(L)/wattspec(L) C We have converted the emergent flux value in each waveband from Monte C Carlo quanta to watts m^-2, and divided by the width of the waveband C to give the spectral irradiance distribution of the emergent flux. C We have also calculated the spectral quantum irradiance distribution. write(100,5265)lamda(L),emergent(L),lamda(L),emreflect(L), * Qemergent(L) 5265 format(f6.1,2X,f8.2,2X,f6.1,2X,f8.6,2X,f8.2) write(110,5270)lamda(L),fphyto(L) 5270 format(f6.1,X,f8.1) C We have divided the emergent flux values in each waveband by the C width of the waveband to give the spectral irradiance. In this case C to reduce the size of the numbers, the waveband width is expressed C in nanometres. write(20,5275)lamda(L) 5275 format(/,5X,f7.1,' nm downward upward') do 5290 i=1,11 k=10*(i-1)+1 WRITE(20,5280)z(i),DSPEC(L,k),UPSPEC(L,k) 5280 FORMAT(5X,f6.1,' m',2X,F10.0,2x,f10.0) 5290 continue 5300 CONTINUE do 5320 i=1,101 do 5310 L=1,165 dspec(L,i)=dflux(L,i)/(1000.0*deltalamda(L)) Wdspec(L,i)=Wperphoton(L)*dspec(L,i)*1000.0 Qdspec(L,i)=5.03*Wdspec(L,i)*lamda(L)/1000.0 upspec(L,i)=upflux(L,i)/(1000.0*deltalamda(L)) Wupspec(L,i)=Wperphoton(L)*upspec(L,i)*1000.0 5310 Qupspec(L,i)=5.03*Wupspec(L,i)*lamda(L)/1000.0 5320 continue C We have calculated the spectral distribution of downward and upward C irradiance, in both energy and quantum terms, at every depth. do 5390 L=1,165 write(120,5350)lamda(L),(exraman(L)/(1000.0*deltalamda(L))), * (emraman(L)/(1000.0*deltalamda(L))),lamda(L), * (cdomfluor(L)/(1000.0*deltalamda(L))), * (chlfluor(L)/(1000.0*deltalamda(L))) 5350 format(f6.1,X,f8.1,X,f8.1,X,f6.1,X,f8.1,X,f8.1) write(130,5360)lamda(L),dspec(L,1),upspec(L,1),dspec(L,2), * upspec(L,2),dspec(L,11),upspec(L,11), * (upspec(L,11)/dspec(L,11)) 5360 format(f6.1,X,f8.1,X,f8.1,X,f8.1,X,f8.1,X,f8.1,X,f8.1,X,f5.4) write(150,5370)lamda(L),Wdspec(L,1),Wupspec(L,1),Wdspec(L,2), * Wupspec(L,2),Wdspec(L,11),Wupspec(L,11), * (Wupspec(L,11)/Wdspec(L,11)) 5370 format(f6.1,X,f8.1,X,f8.1,X,f8.1,X,f8.1,X,f8.1,X,f8.1,X,f5.4) write(160,5380)lamda(L),Qdspec(L,1),Qdspec(L,2), * Qdspec(L,11),Qupspec(L,1),Qupspec(L,2),Qupspec(L,11), * (Qupspec(L,11)/Qdspec(L,11)) 5380 format(f6.1,X,f8.1,X,f8.1,X,f8.1,X,f8.1,X,f8.1,X,f8.1,X,f5.4) 5390 continue total=0.0 DO 5420 I=1,101 DO 5400 L=1,150 TDFLUX(I)=TDFLUX(I)+DFLUX(L,I)*Wperphoton(L) TUPFLX(I)=TUPFLX(I)+UPFLUX(L,I)*Wperphoton(L) 5400 continue C By adding the irradiances in all 165 wavebands, after first dividing C by wavelength to make them proportional to energy, we obtain C numbers proportional to the total downward and upward irradiance C at each depth. Enet(i)=TDFLUX(I)-TUPFLX(I) 5420 continue do 5720 i=2,102 5500 HEAT(I)=Enet(i-1)-Enet(i) C The amount by which the net downward irradiance diminishes across a C given layer of water is equal to the rate of energy absorption by C that layer of water. TOTAL=TOTAL+HEAT(I) C THE RATE OF HEATING WITHIN EACH DEPTH INTERVAL (I.E. THE INTERVAL C ABOVE THE SPECIFIED DEPTH) HAS NOW BEEN CALCULATED. C THE TOTAL RATE OF HEAT ABSORPTION BY THE WATER HAS C ALSO BEEN CALCULATED. 5720 CONTINUE do 5740 i=1,101 5740 propheat(i)=heat(i)/total summixed=0.0 do 5750 i=1,mixed 5750 summixed=summixed+propheat(i) c write(6,*)'summixed= ',summixed write(20,5800) 5800 format(5x,'Depth m heating (proportion)') do 6100 i=0,101 z(i)=i*di write(20,6000)z(i),propheat(i+1) 6000 format(5x,f6.2,5x,f6.5) 6100 continue write(20,6400) 6400 format(/,2x,'wavelength nm ', * 'absorption coefficients m^-1',/) DO 6600 L=1,165 WRITE(20,6500)lamda(L),A(L) 6500 FORMAT(/,5X,f6.1,F8.3,) 6600 CONTINUE c write(20,6700)photons,(100.0*reflect/photons), c * (100.0*escape/photons) c 6700 format(1x,'incident photons= ',f10.0,/, c * '% surface reflection= ', c * f5.2,/,'% sub-surface reflection= ',f5.2) WRITE(20,9600)sunalt 9600 FORMAT(1X,'SOLAR ELEVATION = ',F6.2) WRITE(20,9700)II 9700 FORMAT(5X,'II= ',I5) sumsun=0.0 sumescape=0.0 sumsediment=0.0 do 9800 L=1,165 reflectwatts(L)=reflect(L)*Wperphoton(L) escapewatts(L)=Eu(L)*Wperphoton(L) sedimwatts(L)=sediment(L)*Wperphoton(L) sumreflect=sumreflect+reflectwatts(L) sumescape=sumescape+escapewatts(L) 9800 sumsediment=sumsediment+sedimwatts(L) write(6,*)'watts reflected ',sumreflect,'watts escaping ', * sumescape,'watts absorbed in sediment ', sumsediment c write(6,*)'waterscat ',waterscat wattsin=sumwatts-sumreflect c write(6,*)'wattsin ',wattsin wattsretained=wattsin-sumescape c write(6,*)'wattsretained ',wattsretained wattswater=wattsretained-sumsediment c write(6,*)'wattswater ',wattswater WRITE(20,9900)sumWATTS,wattsin,wattsretained,wattswater 9900 format(1x,'Solar irradiance on surface ',f7.1,' W m^-2 ',/ * 'Solar irradiance penetrating surface ',f7.1,' W m^-2 ',/ * 'Solar irradiance retained by system ',f7.1,' W m^-2 ',/ * 'Solar irradiance absorbed by water ',f7.1,' W m^-2 ') WRITE(20,9820)mixlayer 9820 FORMAT(5X,'mixed layer depth = ',F7.4,' M') PRCENT=sumescape*100.0/sumWATTS WRITE(20,10000)sumescape,PRCENT 10000 FORMAT(5X,'EMERGENT IRRADIANCE FROM BELOW SURFACE = ', * F6.0,' WATTS PER SQUARE METRE',/,F6.2,' % OF INCIDENT') PRCENT=sumsediment*100.0/sumWATTS WRITE(20,10200)sumsediment,PRCENT 10200 FORMAT(5X,'SOLAR IRRADIANCE ABSORBED BY SEDIMENTS = ', * F6.0,' WATTS PER SQUARE METRE',/,F6.2,' % OF INCIDENT') wattsmixed=wattswater*summixed PRCENT=wattsmixed*100.0/sumWATTS WRITE(20,10400)wattsmixed,PRCENT 10400 FORMAT(5X,'SOLAR IRRADIANCE ABSORBED BY mixed layer = ', * F6.0,' WATTS PER SQUARE METRE',/,F6.2,' % OF INCIDENT') 10500 CONTINUE END FUNCTION IRAN(MM,M1,IB,IR) IA=IR IA1=IA/M1 IA0=MOD(IA,M1) IB1=IB/M1 IB0=MOD(IB,M1) MX=IA0*IB1+IA1*IB0 MY=MOD(MX,M1)*M1+IA0*IB0 MULT=MOD(MY,MM) IRAN=MOD((MULT+1),MM) END C Solirrad C John T.O. Kirk, Kirk Marine Optics, 2004 C Fortran program solirrad.for calculates the spectral distribution of C the downward solar radiation stream at sea-level C in 165 wavebands over the range 300-4000 nm. C The separate contributions of the direct solar beam and diffuse sky C radiation, and the angular distribution of the latter, are calculated. C The prevailing meteorological conditions used in the C calculation are specified in the data file, solar.dat. integer day,AM real lamda dimension lamda(167),sol(165),aw(165),aoz(165),au(165), * smonth(13),Tr(165),Ta(165),Taa(165),Tas(165),Tw(165), * Toz(165),Tu(165),H(165),Tdirect(165),Edd(165),Eds(165),EdsR(165), * EdsA(165),deltalamda(165),wattsdd(165),quantdd(165),Pquantdd(165) dimension wattsds(165),quantds(165),Pquantds(165),Qsol(195), * Qdd(165),Qds(165),wattstotal(165) dimension skytheta(9),azimuth(18),cospsi(18,9), * psi(18,9),skyrad(18,9),avskyrad(9),irrad(9), * irradprop(9) real irradprop,irrad open(unit=60,file='sky2.dat',access='sequential', * form='formatted',recl=72) open(3,file='prn') screen=6 OPEN(UNIT=15,FILE='solar.dat',ACCESS='SEQUENTIAL', * FORM='FORMATTED',RECL=12) READ(15,50)month,day,sunalt,P,visrange,Angstrom,AM,RH,W,tozone 50 FORMAT(I2,/,I2,/,F3.1,/,F6.1,/,F4.1,/, * F2.1,/,I2,/,F4.1,/,F2.1,/,F5.3) x=0.0174533 suntheta=90.0-sunalt sunthet=suntheta*x C We have calculated the solar zenith angle in both degrees and radians OPEN(UNIT=10,FILE='METER.DAT',ACCESS='SEQUENTIAL', * FORM='FORMATTED',RECL=7) OPEN(UNIT=20,FILE='OUTPUT1.DAT',ACCESS='SEQUENTIAL', * FORM='FORMATTED',RECL=80) OPEN(UNIT=30,FILE='OUTPUT2.DAT',ACCESS='SEQUENTIAL', * FORM='FORMATTED',RECL=80) OPEN(UNIT=40,FILE='OUTPUT3.DAT',ACCESS='SEQUENTIAL', * FORM='FORMATTED',RECL=80) OPEN(UNIT=50,FILE='OUTPUT4.DAT',ACCESS='SEQUENTIAL', * FORM='FORMATTED',RECL=120) OPEN(UNIT=60,FILE='OUTPUT5.DAT',ACCESS='SEQUENTIAL', * FORM='FORMATTED',RECL=120) OPEN(UNIT=70,FILE='solquanta.dat',ACCESS='SEQUENTIAL', * FORM='FORMATTED',RECL=80) OPEN(UNIT=80,FILE='Pquanta.dat',ACCESS='SEQUENTIAL', * FORM='FORMATTED',RECL=80) OPEN(UNIT=90,FILE='ratio.dat',ACCESS='SEQUENTIAL', * FORM='FORMATTED',RECL=50) OPEN(UNIT=100,FILE='Table1.dat',ACCESS='SEQUENTIAL', * FORM='FORMATTED',RECL=120) C ***** WE NOW CALCULATE THE ANGULAR DISTRIBUTION OF SKY RADIATION *** C We shall not be keeping track of azimuth angle. The zenith angle for C sky radiation is divided into 9 equal 10-degree intervals, and the C proportion of total sky irradiance, averaged over all azimuth angles, C in each of these intervals is calculated, using the empirically based C equation for angular distribution of clear sky radiance of C Harrison & Coombes. do 200 j=1,9 200 skytheta(j)=(10*j-5.0)*x do 300 i=1,18 300 azimuth(i)=(10*i-5.0)*x do 1100 j=1,9 sumskyrad=0.0 write(60,400)(sunthet/x),(skytheta(j)/x) 400 format(5X,'sun degrees ',F5.1,'sky degrees ',F5.1) do 500 i=1,18 cospsi(i,j)=cos(skytheta(j))*cos(sunthet) * +sin(skytheta(j))*sin(sunthet)*cos(azimuth(i)) psi(i,j)=acos(cospsi(i,j)) C Psi is the scattering angle - the angle between the solar beam and C any particular sky radiance direction - calculated with Harrison & C Coombes eq. (3). skyrad(i,j)=(1.63+53.7*exp(-5.94*psi(i,j)) * +2.04*(cospsi(i,j)**2)*cos(sunthet)) * *(1.0-exp(-1.9/cos(skytheta(j)))) * *(1.0-exp(-0.53/cos(sunthet))) C We have calculated the angular distribution of sky radiance for a C set of azimuth and zenith angles, using Harrison & Coombes eq. (6). 500 sumskyrad=sumskyrad+skyrad(i,j) avskyrad(j)=sumskyrad/18.0 write(60,600)avskyrad(j) 600 format(5x,F8.4) C We have averaged sky radiance over all azimuth angles, at each C zenith angle. 1100 continue sumirrad=0.0 write(60,1300)sunthet/x 1300 format(5X,'solar zenith angle ',F4.1) do 1400 j=1,9 irrad(j)=avskyrad(j)*cos(skytheta(j))*sin(skytheta(j)) C From the (azimuthally-averaged) radiance at each zenith angle we have C calculated the corresponding irradiance for that particular C angular interval. sumirrad=sumirrad+irrad(j) 1400 continue do 1500 j=1,9 irradprop(j)=irrad(j)/sumirrad 1500 continue do 1700 j=1,9 write(60,1600)skytheta(j)/x,irradprop(j) 1600 format(5X,'sky zenith angle ',F4.1,'irradiance * proportion ',F5.4) 1700 continue C We have calculated a set of numbers which are proportional to C the irradiance values in 9 angular intervals from 0-10 deg to C 80-90 deg, and then calculated the proportion of the total irradiance C which is in each of these angular intervals. smonth(0)=0 smonth(1)=31 smonth(2)=59 smonth(3)=90 smonth(4)=120 smonth(5)=151 smonth(6)=181 smonth(7)=212 smonth(8)=243 smonth(9)=273 smonth(10)=304 smonth(11)=334 smonth(12)=365 daynumber=smonth(month-1)+day lamda(0)=295 lamda(1)=300 sol(1)=535.9 lamda(2)=305 sol(2)=558.3 lamda(3)=310 sol(3)=622.0 lamda(4)=315 sol(4)=692.7 lamda(5)=320 sol(5)=715.1 lamda(6)=325 sol(6)=832.9 lamda(7)=330 sol(7)=961.9 lamda(8)=335 sol(8)=931.9 lamda(9)=340 sol(9)=900.6 lamda(10)=345 sol(10)=911.3 lamda(11)=350 sol(11)=975.5 lamda(12)=355 sol(12)=1078 lamda(13)=360 sol(13)=1062 lamda(14)=365 sol(14)=1108 lamda(15)=370 sol(15)=1197 lamda(16)=375 sol(16)=1004 lamda(17)=380 sol(17)=1139 lamda(18)=385 sol(18)=1029 lamda(19)=390 sol(19)=1152 lamda(20)=395 sol(20)=1250 lamda(21)=400 sol(21)=1674 lamda(22)=405 sol(22)=1651 lamda(23)=410 sol(23)=1621 lamda(24)=415 sol(24)=1781 lamda(25)=420 sol(25)=1724 lamda(26)=425 sol(26)=1748 lamda(27)=430 sol(27)=1407 lamda(28)=435 sol(28)=1767 lamda(29)=440 sol(29)=1711 lamda(30)=445 sol(30)=1963 lamda(31)=450 sol(31)=2135 lamda(32)=455 sol(32)=2030 lamda(33)=460 sol(33)=2029 lamda(34)=465 sol(34)=1987 lamda(35)=470 sol(35)=1946 lamda(36)=475 sol(36)=2029 lamda(37)=480 sol(37)=2067 lamda(38)=485 sol(38)=1950 lamda(39)=490 sol(39)=1938 lamda(40)=495 sol(40)=2003 lamda(41)=500 sol(41)=1871 lamda(42)=505 sol(42)=1924 lamda(43)=510 sol(43)=1939 lamda(44)=515 sol(44)=1858 lamda(45)=520 sol(45)=1828 lamda(46)=525 sol(46)=1873 lamda(47)=530 sol(47)=1959 lamda(48)=535 sol(48)=1951 lamda(49)=540 sol(49)=1840 lamda(50)=545 sol(50)=1911 lamda(51)=550 sol(51)=1878 lamda(52)=555 sol(52)=1896 lamda(53)=560 sol(53)=1823 lamda(54)=565 sol(54)=1844 lamda(55)=570 sol(55)=1805 lamda(56)=575 sol(56)=1857 lamda(57)=580 sol(57)=1852 lamda(58)=585 sol(58)=1859 lamda(59)=590 sol(59)=1761 lamda(60)=595 sol(60)=1798 lamda(61)=600 sol(61)=1752 lamda(62)=605 sol(62)=1768 lamda(63)=610 sol(63)=1726 lamda(64)=615 sol(64)=1693 lamda(65)=620 sol(65)=1731 lamda(66)=625 sol(66)=1669 lamda(67)=630 sol(67)=1656 lamda(68)=635 sol(68)=1658 lamda(69)=640 sol(69)=1630 lamda(70)=645 sol(70)=1619 lamda(71)=650 sol(71)=1564 lamda(72)=655 sol(72)=1462 lamda(73)=660 sol(73)=1505 lamda(74)=665 sol(74)=1566 lamda(75)=670 sol(75)=1531 lamda(76)=675 sol(76)=1500 lamda(77)=680 sol(77)=1472 lamda(78)=685 sol(78)=1457 lamda(79)=690 sol(79)=1441 lamda(80)=695 sol(80)=1426 lamda(81)=700 sol(81)=1411 lamda(82)=710 sol(82)=1399 lamda(83)=718 sol(83)=1374 lamda(84)=724.4 sol(84)=1373 lamda(85)=740 sol(85)=1298 lamda(86)=752.5 sol(86)=1269 lamda(87)=757.5 sol(87)=1245 lamda(88)=762.5 sol(88)=1223 lamda(89)=767.5 sol(89)=1205 lamda(90)=780 sol(90)=1183 lamda(91)=800 sol(91)=1148 lamda(92)=816 sol(92)=1091 lamda(93)=823.7 sol(93)=1062 lamda(94)=831.5 sol(94)=1038 lamda(95)=840 sol(95)=1022 lamda(96)=860 sol(96)=998.7 lamda(97)=880 sol(97)=947.2 lamda(98)=905 sol(98)=893.2 lamda(99)=915 sol(99)=868.2 lamda(100)=925 sol(100)=829.7 lamda(101)=930 sol(101)=830.3 lamda(102)=937 sol(102)=814 lamda(103)=948 sol(103)=786.9 lamda(104)=965 sol(104)=768.3 lamda(105)=980 sol(105)=767 lamda(106)=993.5 sol(106)=757.6 lamda(107)=1040 sol(107)=688.1 lamda(108)=1070 sol(108)=640.7 lamda(109)=1100 sol(109)=606.2 lamda(110)=1120 sol(110)=585.9 lamda(111)=1130 sol(111)=570.2 lamda(112)=1145 sol(112)=564.1 lamda(113)=1161 sol(113)=544.2 lamda(114)=1170 sol(114)=533.4 lamda(115)=1200 sol(115)=501.6 lamda(116)=1240 sol(116)=477.5 lamda(117)=1270 sol(117)=442.7 lamda(118)=1290 sol(118)=440 lamda(119)=1320 sol(119)=416.8 lamda(120)=1350 sol(120)=391.4 lamda(121)=1395 sol(121)=358.9 lamda(122)=1442.5 sol(122)=327.5 lamda(123)=1462.5 sol(123)=317.5 lamda(124)=1477 sol(124)=307.3 lamda(125)=1497 sol(125)=300.4 lamda(126)=1520 sol(126)=292.8 lamda(127)=1539 sol(127)=275.5 lamda(128)=1558 sol(128)=272.1 lamda(129)=1578 sol(129)=259.3 lamda(130)=1592 sol(130)=246.9 lamda(131)=1610 sol(131)=244 lamda(132)=1630 sol(132)=243.5 lamda(133)=1646 sol(133)=234.8 lamda(134)=1678 sol(134)=220.5 lamda(135)=1740 sol(135)=190.8 lamda(136)=1800 sol(136)=171.1 lamda(137)=1860 sol(137)=144.5 lamda(138)=1920 sol(138)=135.7 lamda(139)=1960 sol(139)=123 lamda(140)=1985 sol(140)=123.8 lamda(141)=2005 sol(141)=113 lamda(142)=2035 sol(142)=108.5 lamda(143)=2065 sol(143)=97.5 lamda(144)=2100 sol(144)=92.4 lamda(145)=2148 sol(145)=82.4 lamda(146)=2198 sol(146)=74.6 lamda(147)=2270 sol(147)=68.3 lamda(148)=2360 sol(148)=63.8 lamda(149)=2450 sol(149)=49.5 lamda(150)=2500 sol(150)=48.5 lamda(151)=2600 sol(151)=38.6 lamda(152)=2700 sol(152)=36.6 lamda(153)=2800 sol(153)=32.0 lamda(154)=2900 sol(154)=28.1 lamda(155)=3000 sol(155)=24.8 lamda(156)=3100 sol(156)=22.1 lamda(157)=3200 sol(157)=19.6 lamda(158)=3300 sol(158)=17.5 lamda(159)=3400 sol(159)=15.7 lamda(160)=3500 sol(160)=14.1 lamda(161)=3600 sol(161)=12.7 lamda(162)=3700 sol(162)=11.5 lamda(163)=3800 sol(163)=10.4 lamda(164)=3900 sol(164)=9.5 lamda(165)=4000 sol(165)=8.6 lamda(166)=4100 C We have entered the values of extraterrestrial solar irradiance C in 165 wavebands from the data of Frohlich & Wehrli (1981) C and Neckel & Labs (1981). do 1800 L=1,165 1800 deltalamda(L)=(lamda(L+1)-lamda(L-1))/2000.0 C We have calculated the width of each waveband, in microns. aw(55)=0.034 aw(56)=0.012 aw(57)=0.0 aw(58)=0.003 aw(59)=0.7 aw(60)=0.361 aw(61)=0.023 aw(62)=0.0 aw(63)=0.0 aw(64)=0.0 aw(65)=0.0 aw(66)=0.0 aw(67)=0.005 aw(68)=0.0 aw(69)=0.001 aw(70)=0.011 aw(71)=0.173 aw(72)=0.142 aw(73)=0.008 aw(74)=0.001 aw(75)=0.0 aw(76)=0.0 aw(77)=0.0 aw(78)=0.002 aw(79)=0.001 aw(80)=0.53 aw(81)=0.602 aw(82)=0.0125 aw(83)=1.8 aw(84)=2.5 aw(85)=0.061 aw(86)=0.0008 aw(87)=0.0001 aw(88)=0.00001 aw(89)=0.00001 aw(90)=0.0006 aw(91)=0.036 aw(92)=1.6 aw(93)=2.5 aw(94)=0.5 aw(95)=0.155 aw(96)=0.00001 aw(97)=0.0026 aw(98)=7.0 aw(99)=5.0 aw(100)=5.0 aw(101)=27.0 aw(102)=55.0 aw(103)=45.0 aw(104)=4.0 aw(105)=1.48 aw(106)=0.1 aw(107)=0.00001 aw(108)=0.001 aw(109)=3.2 aw(110)=115.0 aw(111)=70.0 aw(112)=75.0 aw(113)=10.0 aw(114)=5.0 aw(115)=2.0 aw(116)=0.002 aw(117)=0.002 aw(118)=0.1 aw(119)=4.0 aw(120)=200.0 aw(121)=1000.0 aw(122)=185.0 aw(123)=80.0 aw(124)=80.0 aw(125)=12.0 aw(126)=0.16 aw(127)=0.002 aw(128)=0.0005 aw(129)=0.0001 aw(130)=0.00001 aw(131)=0.0001 aw(132)=0.001 aw(133)=0.01 aw(134)=0.036 aw(135)=1.1 aw(136)=130.0 aw(137)=1000.0 aw(138)=500.0 aw(139)=100.0 aw(140)=4.0 aw(141)=2.9 aw(142)=1.0 aw(143)=0.4 aw(144)=0.22 aw(145)=0.25 aw(146)=0.33 aw(147)=0.5 aw(148)=4.0 aw(149)=80.0 aw(150)=310.0 aw(151)=15000.0 aw(152)=22000.0 aw(153)=8000.0 aw(154)=650.0 aw(155)=240.0 aw(156)=230.0 aw(157)=100.0 aw(158)=120.0 aw(159)=19.5 aw(160)=3.6 aw(161)=3.1 aw(162)=2.5 aw(163)=1.4 aw(164)=0.17 aw(165)=0.0045 aoz(1)=10.0 aoz(2)=4.8 aoz(3)=2.7 aoz(4)=1.35 aoz(5)=0.8 aoz(6)=0.38 aoz(7)=0.16 aoz(8)=0.075 aoz(9)=0.04 aoz(10)=0.019 aoz(11)=0.007 aoz(12)=0.002 aoz(13)=0.0 aoz(14)=0.0 aoz(15)=0.0 aoz(16)=0.0 aoz(17)=0.0 aoz(18)=0.0 aoz(19)=0.0 aoz(20)=0.0 aoz(21)=0.0 aoz(22)=0.0 aoz(23)=0.0 aoz(24)=0.0 aoz(25)=0.0 aoz(26)=0.0 aoz(27)=0.001 aoz(28)=0.002 aoz(29)=0.003 aoz(30)=0.003 aoz(31)=0.003 aoz(32)=0.005 aoz(33)=0.008 aoz(34)=0.008 aoz(35)=0.007 aoz(36)=0.012 aoz(37)=0.018 aoz(38)=0.02 aoz(39)=0.018 aoz(40)=0.022 aoz(41)=0.028 aoz(42)=0.038 aoz(43)=0.04 aoz(44)=0.038 aoz(45)=0.045 aoz(46)=0.054 aoz(47)=0.063 aoz(48)=0.07 aoz(49)=0.072 aoz(50)=0.077 aoz(51)=0.084 aoz(52)=0.089 aoz(53)=0.096 aoz(54)=0.107 aoz(55)=0.116 aoz(56)=0.119 aoz(57)=0.117 aoz(58)=0.111 aoz(59)=0.108 aoz(60)=0.115 aoz(61)=0.124 aoz(62)=0.124 aoz(63)=0.118 aoz(64)=0.111 aoz(65)=0.104 aoz(66)=0.097 aoz(67)=0.09 aoz(68)=0.083 aoz(69)=0.077 aoz(70)=0.07 aoz(71)=0.064 aoz(72)=0.059 aoz(73)=0.055 aoz(74)=0.051 aoz(75)=0.046 aoz(76)=0.042 aoz(77)=0.038 aoz(78)=0.034 aoz(79)=0.03 aoz(80)=0.025 aoz(81)=0.022 aoz(82)=0.018 aoz(83)=0.015 aoz(84)=0.012 aoz(85)=0.01 aoz(86)=0.008 aoz(87)=0.007 aoz(88)=0.006 aoz(89)=0.005 aoz(90)=0.0 do 2000 L=91,165 aoz(L)=0.0 2000 continue do 3000 L=1,78 au(L)=0.0 3000 continue au(79)=0.15 do 4000 L=80,87 au(L)=0.0 4000 continue au(88)=4.0 au(89)=0.35 do 5000 L=90,115 au(L)=0.0 5000 continue au(116)=0.05 au(117)=0.3 au(118)=0.02 au(119)=0.0002 au(120)=0.00011 au(121)=0.00001 au(122)=0.05 au(123)=0.011 au(124)=0.005 au(125)=0.0006 au(126)=0.00 au(127)=0.005 au(128)=0.13 au(129)=0.04 au(130)=0.06 au(131)=0.13 au(132)=0.001 au(133)=0.0014 au(134)=0.0001 au(135)=0.00001 au(136)=0.00001 au(137)=0.0001 au(138)=0.001 au(139)=4.3 au(140)=0.2 au(141)=21.0 au(142)=0.13 au(143)=1.0 au(144)=0.08 au(145)=0.001 au(146)=0.00038 au(147)=0.001 au(148)=0.0005 au(149)=0.00015 au(150)=0.00014 au(151)=0.00066 au(152)=100.0 au(153)=150.0 au(154)=0.13 au(155)=0.0095 au(156)=0.001 au(157)=0.8 au(158)=1.9 au(159)=1.3 au(160)=0.075 au(161)=0.01 au(162)=0.00195 au(163)=0.004 au(164)=0.29 au(165)=0.025 C We have entered the values of the absorption coefficients of water C vapour, ozone, and uniformly mixed gases. do 5020 L=1,165 write(100,5010)lamda(L),sol(L),aw(L),aoz(L),au(L) 5010 format(f6.1,20X,f6.1,20X,f10.4,20X,f6.3,20X,f10.5) 5020 continue orbit=(1+0.0167*cos(6.2832*(daynumber-3)/365))**2 do 10010 L=1,165 H(L)=sol(L)*orbit 10010 continue C We have corrected the solar irradiance for variation C in Earth-Sun distance using eq. 2.13. C **** EFFECTIVE AIR MASS **** solpath=1/(cos(sunthet)+0.15*(93.885-sunthet)**(-1.253)) solpathP=solpath*P/1013 c write(6,*)'solpath',solpath,' ','solpathP',solpathP C We have calculated the effective path of the solar beam through C the atmosphere using eq. 2.12 and the pressure-corrected quantity C from eq. 2.18. C **** RAYLEIGH SCATTERING **** do 10030 L=1,165 lamda(L)=lamda(L)/1000.0 10030 Tr(L)=exp(-solpathP/((115.6406*(lamda(L)**4) * -(1.335*lamda(L)**2)))) C We have re-expressed wavelengths in microns, and calculated the C transmittance functions due to Rayleigh scattering, using eq. 2.22. C **** AEROSOL SCATTERING AND ABSORPTION **** beta=3.91*(0.55**Angstrom)/visrange omega=(0.972-0.0032*AM)*exp(0.000306*RH) c write(6,*)'beta',beta,' ','omega',omega C We have calculated the turbidity coefficient using eq. 2.25, and C the single-scattering albedo using eq. 2.31. do 10060 L=1,165 Ta(L)=exp(-beta*(lamda(L)**(-Angstrom))*solpath) Taa(L)=exp(-(1.0-omega)*beta*lamda(L)**(-Angstrom)*solpath) 10060 Tas(L)=exp(-omega*beta*lamda(L)**(-Angstrom)*solpath) C We have calculated the transmittance functions due to aerosol C scattering and absorption combined, and to absorption and C scattering separately, using eqs. 2.26, 2.29 and 2.30. C **** ABSORPTION BY WATER VAPOUR **** do 10080 L=1,165 10080 Tw(L)=exp((-0.2385*aw(L)*W*solpath)/ * ((1+20.07*aw(L)*W*solpath)**(0.45))) C We have calculated the transmittance functions due to absorption by C water vapour using eq. 2.32. C ***** ABSORPTION BY OZONE **** do 10100 L=1,165 10100 Toz(L)=exp(-aoz(L)*tozone*solpath) C We have calculated the transmittance functions due to absorption C by ozone using eq. 2.33. C **** ABSORPTION BY UNIFORMLY MIXED GASES **** do 10120 L=1,165 Tu(L)=exp((-1.41*au(L)*solpathP)/ * ((1.0+118.3*au(L)*solpathP)**(0.45))) if(au(L).eq.0.0)Tu(L)=1.0 10120 continue C We have calculated the transmittance functions due to the uniformly C mixed gases using eq. 2.35. C **** THE DIRECT SOLAR FLUX **** sumquantdd=0.0 do 10150 L=1,165 Tdirect(L)=Tr(L)*Ta(L)*Tw(L)*Toz(L)*Tu(L) C This is the transmittance function for the direct solar beam due C to all components of the atmosphere. Edd(L)=H(L)*cos(sunthet)*Tdirect(L) wattsdd(L)=Edd(L)*deltalamda(L) C We multiply the spectral irradiance (in watts m^-2 microns^-1), C i.e. the spectral concentration of energy, C at each wavelength by the width of the waveband for which that C wavelength is the centre, to give the irradiance (watts m^-2) in C that waveband. write(40,10140)Edd(L) 10140 format(F6.1) quantdd(L)=5.03*wattsdd(L)*lamda(L) C To calculate the number of quanta (m^-2 s^-1 X 10^18) C in each waveband we multiply the energy irradiance C in that waveband by 5.03 X wavelength in microns (Kirk, 1994; p. 5). C We then add all these numbers up and calculate the proportion of the C sum which is in each waveband, and in this way obtain the proportion C of the total quanta in the direct solar beam which is in C each of the 165 wavebands. sumquantdd=sumquantdd+quantdd(L) 10150 continue do 10160 L=1,165 10160 Pquantdd(L)=quantdd(L)/sumquantdd C From the total quanta in each waveband we have now calculated the C quantum spectral irradiance, i.e. the spectral concentration of C radiant energy (per unit wavelength) C expressed in quantum terms, as a function of wavelength. C We have calculated the irradiance due to the direct solar beam C on a horizontal surface at sea level, and also the proportion C of the total quanta in the direct beam irradiance in each of the C 165 wavebands. C **** SKY RADIATION **** asymmfac=-0.1417*Angstrom+0.82 if(Angstrom.gt.1.2)asymmfac=0.65 Fa=(1+asymmfac)/2+0.018 C We have calculated the forward to total scattering ratio for C the aerosol using eqs. 2.48 and 2.45. sumquantds=0.0 do 10300 L=1,165 EdsR(L)=H(L)*cos(sunthet)*Taa(L)*Tw(L)*Toz(L)*Tu(L) * *(1.0-(Tr(L)**0.95))*0.5 EdsA(L)=H(L)*cos(sunthet)*(Tr(L)**1.5)*Taa(L)*Tw(L)* * Toz(L)*Tu(L)*(1.0-Tas(L))*Fa Eds(L)=EdsR(L)+EdsA(L) C We have calculated the irradiance on a horizontal surface at C sea level due to Rayleigh scattering and aerosol scattering C from the sky. wattsds(L)=Eds(L)*deltalamda(L) C We multiply the spectral irradiance (in watts m^-2 microns^-1), C i.e. the spectral concentration of energy, at C each wavelength by the width of the waveband for which that C wavelength is the centre, to give the irradiance (watts m^-2) in C that waveband. write(60,10200)Eds(L) 10200 format(f6.1) quantds(L)=5.03*wattsds(L)*lamda(L) C To calculate the number of quanta (m^-2 s^-1 X 10^18) C in each waveband we multiply the energy irradiance C in that waveband by 5.03 X wavelength (microns). C We then add all these numbers up and calculate the proportion of the C sum which is in each waveband, and in this way obtain the proportion C of the total quanta in the diffuse solar flux which is in C each of the 165 wavebands. sumquantds=sumquantds+quantds(L) 10300 continue do 10400 L=1,165 10400 Pquantds(L)=quantds(L)/sumquantds sumwatts=0.0 do 10700 L=1,165 wattstotal(L)=wattsdd(L)+wattsds(L) sumwatts=sumwatts+wattsdd(L)+wattsds(L) Qsol(L)=5.03*H(L)*lamda(L) Qdd(L)=5.03*Edd(L)*lamda(L) Qds(L)=5.03*Eds(L)*lamda(L) write(70,10500)1000.0*lamda(L),Qsol(L),Qdd(L),Qds(L) 10500 format(f8.1,X,f6.1,X,f6.1,X,f6.1) C From the energy spectral irradiance in each waveband we have now C calculated the quantum spectral irradiance, i.e. the spectral C concentration of radiant energy (per unit wavelength) C expressed in quantum terms, as a function of wavelength. write(80,10600)1000.0*lamda(L),Pquantdd(L),Pquantds(L), * wattsdd(L),wattstotal(L) 10600 format(f6.1,X,f10.9,X,f10.9,X,f8.3,X,f8.3) 10700 continue ratio=sumquantds/sumquantdd write(90,10800)ratio,sumwatts 10800 format(f6.3,/,f8.1) do 30020 L=1,165 c write(20,30010)1000.0*lamda(L) 30010 format(f5.0) 30020 continue do 30040 L=1,165 write(30,30030)L,sol(L),H(L) 30030 format(i3,X,f6.1,X,f6.1) 30040 continue do 30060 L=1,165 write(50,30050)lamda(L),Tr(L),Ta(L),Tw(L),Toz(L),Tu(L) 30050 format(5X,F6.3,6X,F8.6,6X,F8.6,6X,F8.6,6X,F8.6,6X,F8.6) 30060 continue end C PROGRAM Upscat C John T.O. Kirk, Kirk Marine Optics, 1988 & 2004 C This is a Fortran program for the Monte Carlo simulation of the C monochromatic light field established within any natural water of C specified optical properties. The water surface is assumed to be flat C and to be illuminated by a parallel stream of photons at some C specified angle. C The values of absorption and scattering coefficients are specified, C and a scattering phase function with equal (5 degree) angular C intervals is supplied. C A completely absorbing bottom is assumed to be present at C 20 depth intervals. C The data for each run are first entered into the data file C upscatdata.dat. This has 6 lines: C 1. the absorption coefficient of the water (m^-1) C 2. the scattering coefficient of the water (m^-1) C 3. the solar altitude (degrees). N.B. for zenith sun, make C the altitude just less than 90.0 deg., e.g. 89.99 deg. C 4. a large integer to initialize the random number generator C 5. the depth interval (m) C 6. the number of photons to be used in the simulation C A copy of whichever phase function is to be used should be made C under the name, 'phasefn.dat' INTEGER ZSTART,ZFINSH,BOTTOM DOUBLE PRECISION DFLUX,DSCAL,rad DIMENSION DFLUX(21),UPFLUX(21),Z(21),SINTH(36),COSANG(36), * CSINTH(36),REFLEC(21),DSCAL(21),UPSCAL(21),NAME(20), * SCAL(21),AVCOS(21),DAVCOS(21),UAVCOS(21),RAD(21,36), * ANGLE(36),SINANG(36),UPRAD(21,36),UPCOS(21),UP(21),SCUP(21), * scatt(72) OPEN(UNIT=10,FILE='PHASEFN.DAT',ACCESS='sequential', * FORM='FORMATTED',RECL=7) OPEN(UNIT=30,FILE='upscat.DAT',ACCESS='sequential', * FORM='FORMATTED',RECL=150) OPEN(UNIT=20,FILE='upscatdata.DAT',ACCESS='sequential', * FORM='FORMATTED',RECL=12) OPEN(UNIT=40,FILE='test.DAT',ACCESS='sequential', * FORM='FORMATTED',RECL=50) DO 2 I=1,72 READ(10,1)SCATT(I) 1 FORMAT(F7.5) 2 CONTINUE DO 4 I=1,36 SINTH(I)=SCATT(I) CSINTH(I)=SCATT(I+36) 4 CONTINUE C THE VALUES OF THE SCATTERING ANGLE AND THE CUMULATIVE C DISTRIBUTION HAVE BEEN READ IN read(20,10)a,b,sun,ii,di,pm 10 format(f6.3,/,f6.3,/,f5.2,/,i8,/,f7.4,/,f14.0) IJK=II SUN=SUN*1.5708/90.0 SUNINW=ATAN(SQRT(1.0-0.5625*COS(SUN)*COS(SUN))/(0.75*COS(SUN))) SINSNW=SIN(SUNINW) BOTTOM=21 m=100000000 m1=10000 ib=31415821 ir=iran(m,m1,ib,ir) n=0 40 n=n+1 if(n.gt.ii)go to 50 ir=iran(m,m1,ib,ir) go to 40 50 continue do 102 i=1,21 DO 101 J=1,36 UPRAD(I,J)=0.0 101 RAD(I,J)=0.0 SCAL(I)=0.0 DSCAL(I)=0.0 UPSCAL(I)=0.0 REFLEC(I)=0.0 DAVCOS(I)=0.0 UAVCOS(I)=0.0 UPCOS(I)=0.0 UP(I)=0.0 SCUP(I)=0.0 DFLUX(I)=0.0 102 UPFLUX(I)=0.0 ESCAPE=0 PN=0 C A photon enters the water 110 PN=PN+1.0 IF(PN.GT.PM)GO TO 420 140 SINALF=SINSNW IALFA=INT(SUNINW/0.0872665)+1 150 DEPTH=0.0 NSCAT=0 160 NSCAT=NSCAT+1 162 ir=iran(m,m1,ib,ir) r=real(ir)/100000000.0 IF(R.EQ.1.0)GO TO 162 PATH=-LOG(1.0-R)/(A+B) C HAVE DETERMINED THE PATHLENGTH OF THE PHOTON C BEFORE ITS NEXT INTERACTION IF(NSCAT.EQ.1) GO TO 210 C IF PHOTON HAS BEEN SCATTERED DETERMINE NEW TRAJECTORY 170 ir=iran(m,m1,ib,ir) r=real(ir)/100000000.0 KA=0 180 KA=KA+1 IF(R.GT.CSINTH(KA))GO TO 180 SINTHE=SINTH(KA) ir=iran(m,m1,ib,ir) r=real(ir)/100000000.0 PHI=R*6.2832 COSTHE=SQRT(1.0-SINTHE*SINTHE) IF(KA.GT.18)COSTHE=-COSTHE COSALF=SQRT(1.0-SINALF*SINALF) SINALF=COSTHE*SINALF+SINTHE*COSALF*COS(PHI) NIALFA=IALFA C We are noting the previous IALFA ofthe photon so C that we know whether it was travelling up or C down before the present scattering event. IALFA=INT(ABS(ASIN(SINALF)/0.0872665))+1 C SIN ALPHA HAS BEEN CALCULATED FOR THE NEW TRAJECTORY IF(SINALF.LE.0.0)GO TO 270 C POSITIVE ALPHA MEANS PHOTON TRAVELLING DOWNWARDS 210 ZSTART=DEPTH/DI+2 ZFINSH=(DEPTH+PATH*SINALF)/DI+1 DEPTH=DEPTH+PATH*SINALF IF(ZFINSH.LT.ZSTART)GO TO 410 IF(ZFINSH.LT.BOTTOM)GO TO 240 220 DO 228 I=ZSTART,BOTTOM DFLUX(I)=DFLUX(I)+1.0 228 RAD(I,IALFA)=RAD(I,IALFA)+1.0 GO TO 110 C NEW PHOTON REQUIRED(PREVIOUS ONE REACHED BOTTOM) 240 DO 248 I=ZSTART,ZFINSH DFLUX(I)=DFLUX(I)+1.0 248 RAD(I,IALFA)=RAD(I,IALFA)+1.0 260 GO TO 410 C DETERMINE WHETHER THIS PHOTON IS TO BE SCATTERED AGAIN OR C TERMINATED C NEGATIVE ALPHA MEANS PHOTON TRAVELLING UPWARDS 270 IALFA=IALFA+18 IF(NIALFA.GT.18)GO TO 280 C We only wish to study the upward scattering of photons C which were previously travelling downwards. ZSTART=DEPTH/DI+1 UPRAD(ZSTART,IALFA)=UPRAD(ZSTART,IALFA)+1.0 C In this step, for each layer of water we keep a record C of how many photons are scattered into each angular C interval. This gives us what may be regarded as essentially C equivalent to the angular distribution of an upward C irradiance which is due to photons scattered upwards in this C water layer. UP(ZSTART)=UP(ZSTART)+1.0 C In this step we accumulate the number of photons scattered C upwards in that water layer at all angles, ending up with a number C proportional to the notional irradiance due to photons C scattered upwards within that layer. 280 ZSTART=DEPTH/DI+1 ZFINSH=(DEPTH+PATH*SINALF)/DI+2 DEPTH=ABS(DEPTH+PATH*SINALF) IF(ZFINSH.GT.ZSTART)GO TO 410 IF(1.LT.ZFINSH)GO TO 390 C PHOTON REACHES SURFACE 290 DO 300 I=1,ZSTART RAD(I,IALFA)=RAD(I,IALFA)+1.0 300 UPFLUX(I)=UPFLUX(I)+1.0 IF(ABS(SINALF).GE.0.6613)GO TO 370 C INTERNAL REFLECTION AT SURFACE IALFA=IALFA-18 ZFINSH=DEPTH/DI+1 IF(ZFINSH.LT.BOTTOM)GO TO 350 330 DO 338 I=1,BOTTOM DFLUX(I)=DFLUX(I)+1.0 338 RAD(I,IALFA)=RAD(I,IALFA)+1.0 GO TO 110 C NEW PHOTON REQUIRED (PREVIOUS ONE REACHED BOTTOM) 350 DO 358 I=1,ZFINSH DFLUX(I)=DFLUX(I)+1.0 358 RAD(I,IALFA)=RAD(I,IALFA)+1.0 SINALF=-SINALF C SIGN OF ALPHA MUST BE REVERSED SINCE PHOTON IS NOW C TRAVELLING DOWNWARDS AGAIN GO TO 410 C IS PHOTON TO BE SCATTERED AGAIN OR TERMINATED? C PHOTON NOT REFLECTED AT SURFACE 370 ESCAPE=ESCAPE+1.0 GO TO 110 C NEW PHOTON REQUIRED 390 DO 400 I=ZFINSH,ZSTART RAD(I,IALFA)=RAD(I,IALFA)+1.0 400 UPFLUX(I)=UPFLUX(I)+1.0 410 ir=iran(m,m1,ib,ir) r=real(ir)/100000000.0 IF(R.LE.(A/(A+B)))GO TO 110 415 GO TO 160 C HAVE DETERMINED WHETHER THAT PHOTON WAS ABSORBED OR IS TO C BE SCATTERED AGAIN 420 WRITE(30,430)PM,ESCAPE 430 FORMAT(10X,'TOTAL PHOTONS IN',25X,F10.0,/, *10X,'TOTAL PHOTONS OUT',24X,F10.0) DFLUX(1)=DFLUX(1)+PM C THE ANGULAR DISTRIBUTION OF IRRADIANCE IS NOW USED TO CALCULATE C THE DOWNWARDS AND UPWARDS SCALAR IRRADIANCE AT EACH DEPTH DO 436 I=1,21 DO 432 J=1,18 432 DSCAL(I)=DSCAL(I)+RAD(I,J)/SIN((J*5.0-2.5)*3.14159/180.0) DO 434 J=19,36 434 UPSCAL(I)=UPSCAL(I)+RAD(I,J)/SIN(((J-18)*5.0-2.5)*3.14159/180.0) 436 CONTINUE DSCAL(1)=DSCAL(1)+PM/SINSNW C THE AVERAGE COSINES FOR THE DOWNWELLING AND UPWELLING STREAMS C ARE NOW CALCULATED FROM THE RATIOS OF IRRADIANCE TO SCALAR C IRRADIANCE IN THE UPPER AND LOWER HEMISPHERES dscal$2=0.0 i=2 do 439 j=1,18 dscal$2=dscal$2+RAD(I,J)/SIN((J*5.0-2.5)*3.14159/180.0) write(40,438)rad(i,j),(RAD(I,J)/SIN((J*5.0-2.5)*3.14159/180.0)) 438 format(f12.1,X,f12.1) 439 continue I=0 440 I=I+1 IF(I.EQ.22)GO TO 450 IF(DFLUX(I).LT.1.0)GO TO 450 DAVCOS(I)=DFLUX(I)/DSCAL(I) GO TO 440 450 CONTINUE I=0 460 I=I+1 IF (I.EQ.22)GO TO 470 IF(UPFLUX(I).LT.1.0)GO TO 470 UAVCOS(I)=UPFLUX(I)/UPSCAL(I) GO TO 460 470 CONTINUE L=0 482 L=L+1 IF(L.EQ.21)GO TO 483 IF(UPFLUX(L).LT.100.0) GO TO 483 REFLEC(L)=UPFLUX(L)/DFLUX(L) GO TO 482 483 CONTINUE C THE AVERAGE VALUES IN THE EUPHOTIC ZONE OF VERTICAL ATTENUATION C COEFFICIENTS FOR IRRADIANCE AND SCALAR IRRADIANCE ARE NOW C CALCULATED NE=0 485 NE=NE+1 IF(NE.EQ.21) GO TO 487 IF(100.0*DFLUX(NE).GT.DFLUX(1))GO TO 485 IF(DFLUX(NE).LT.1.0)NE=NE-1 IF(NE.LT.2)GO TO 535 487 SUMQ=0.0 SUMZ=0.0 SUMQZ=0.0 SUMZSQ=0.0 DO 520 I=1,NE X=LOG(DFLUX(I)) SUMQ=SUMQ+X Z(I)=DI*(I-1) SUMZ=SUMZ+Z(I) SUMQZ=SUMQZ+X*Z(I) 520 SUMZSQ=SUMZSQ+Z(I)*Z(I) ATTENK=-(SUMQZ-SUMQ*SUMZ/NE)/(SUMZSQ-SUMZ*SUMZ/NE) WRITE(30,530)ATTENK,B 530 FORMAT(5X,'AVERAGE VERTICAL ATTENUATION COEFFICIENT IS',4X,F7.3,/, * 2X,'WHEN SCATTERING COEFFICIENT IS',4X,F6.3) 535 CONTINUE DO 550 I=1,21 550 SCAL(I)=DSCAL(I)+UPSCAL(I) 570 NE=0 580 NE=NE+1 IF(NE.EQ.21)GO TO 590 IF(100.0*SCAL(NE).GT.SCAL(1))GO TO 580 IF(SCAL(NE).LT.1.0)NE=NE-1 IF(NE.LT.2)GO TO 608 590 SUMQ=0.0 SUMZ=0.0 SUMQZ=0.0 SUMZSQ=0.0 DO 600 I=1,NE X=LOG(SCAL(I)) SUMQ=SUMQ+X Z(I)=DI*(I-1) SUMZ=SUMZ+Z(I) SUMQZ=SUMQZ+X*Z(I) 600 SUMZSQ=SUMZSQ+Z(I)*Z(I) ATTENK=-(SUMQZ-SUMQ*SUMZ/NE)/(SUMZSQ-SUMZ*SUMZ/NE) WRITE(30,605)ATTENK,B 605 FORMAT(X,'AVERAGE VERT. ATTENUATION COEFF. FOR SCALAR IRRADIANCE * IS',4X,F7.3,/,2X,'WHEN SCATTERING COEFFICIENT IS',4X,F6.3) C THE AVERAGE COSINE FOR ALL PHOTONS IS NOW OBTAINED AT EACH C DEPTH FROM THE RATIO OF NET DOWNWARD IRRADIANCE TO SCALAR C IRRADIANCE 608 DO 610,I=1,21 610 AVCOS(I)=0.0 I=0 612 I=I+1 IF(I.EQ.22)GO TO 615 IF(DFLUX(I).LT.1.0)GO TO 615 AVCOS(I)=(DFLUX(I)-UPFLUX(I))/SCAL(I) GO TO 612 615 CONTINUE 629 DO 630 I=1,21 630 Z(I)=DI*(I-1) WRITE(30,632) 632 FORMAT(X,'DEPTH',8X,'IRRAD. DOWN',2X,'IRRAD. UP',2X, * 'REFLECTANCE',X,'SCALAR IRRAD.',2X,'AV. COSINE',3X, * 'AV. COSINE DOWN',X,'AV.COSINE UP',4X,/) DO 636 I=1,21 WRITE(30,634)(Z(I)),(DFLUX(I)),(UPFLUX(I)),(REFLEC(I)),(SCAL(I)) * ,(AVCOS(I)),(DAVCOS(I)),(UAVCOS(I)) 634 FORMAT(X,F6.2,X,'M',F16.0,F11.0,F13.4,F14.0,F12.3,2F13.3,F10.3) 636 CONTINUE IALFA=INT(SUNINW/0.0872665)+1 RAD(1,IALFA)=RAD(1,IALFA)+PM WRITE(30,637) 637 FORMAT(/,X,' IRRADIANCE OCCURRING IN EACH ANGULAR INTERVAL') DO 638 I=1,18 638 ANGLE(I)=I*5.0-2.5 DO 639 I=19,36 639 ANGLE(I)=-((I-18)*5.0-2.5) WRITE(30,642) 642 FORMAT(X,'ANGLE TO HORIZONTAL') WRITE(30,643)(ANGLE(I),I=1,9) 643 FORMAT(/,X,'DEGREES',2X,9F12.1,/) DO 645 I=1,21 WRITE(30,644)(Z(I)),(RAD(I,J),J=1,9) 644 FORMAT(X,F6.2,X,'M',X,9F12.0) 645 CONTINUE WRITE(30,646)(ANGLE(I),I=10,18) 646 FORMAT(/,X,'DEGREES',2X,9F12.1,/) DO 648 I=1,21 WRITE(30,647)(Z(I)),(RAD(I,J),J=10,18) 647 FORMAT(X,F6.2,X,'M',X,9F12.0) 648 CONTINUE WRITE(30,655)(ANGLE(I),I=19,27) 655 FORMAT(/,X,'DEGREES',2X,9F12.1,/) DO 670 I=1,21 WRITE(30,660)(Z(I)),(RAD(I,J),J=19,27) 660 FORMAT(X,F6.2,X,'M',X,9F12.0) 670 CONTINUE WRITE(30,675)(ANGLE(I),I=28,36) 675 FORMAT(/,X,'DEGREES',2X,9F12.1,/) DO 690 I=1,21 WRITE(30,680)(Z(I)),(RAD(I,J),J=28,36) 680 FORMAT(X,F6.2,X,'M',9F12.0) 690 CONTINUE C THE (RELATIVE) RADIANCE DISTRIBUTION IS NOW CALCULATED DO 749 I=1,36 SINANG(I)=SIN(ABS(ANGLE(I)*1.5708/90.0)) 749 COSANG(I)=COS(ABS(ANGLE(I)*1.5708/90.0)) WRITE(30,750) 750 FORMAT(/,X,'(RELATIVE) RADIANCE DISTRIBUTION ') DO 751 I=1,21 DO 751 J=1,36 751 RAD(I,J)=0.1*RAD(I,J)/(COSANG(J)*SINANG(J) * *6.2832*(5.0*3.14159/180.0)) WRITE(30,752)(ANGLE(I),I=1,9) 752 FORMAT(/,X,'DEGREES',2X,9F12.1,/) DO 754 I=1,21 WRITE(30,753)(Z(I)),(RAD(I,J),J=1,9) 753 FORMAT(X,F6.2,X,'M',X,9F12.1) 754 CONTINUE WRITE(30,755)(ANGLE(I),I=10,18) 755 FORMAT(/,X,'DEGREES',2X,9F12.1,/) DO 758 I=1,21 WRITE(30,756)(Z(I)),(RAD(I,J),J=10,18) 756 FORMAT(X,F6.2,X,'M',X,9F12.1) 758 CONTINUE WRITE(30,855)(ANGLE(I),I=19,27) 855 FORMAT(/,X,'DEGREES',2X,9F12.1,/) DO 870 I=1,21 WRITE(30,860)(Z(I)),(RAD(I,J),J=19,27) 860 FORMAT(X,F6.2,X,'M',X,9F12.1) 870 CONTINUE WRITE(30,875)(ANGLE(I),I=28,36) 875 FORMAT(/,X,'DEGREES',2X,9F12.1,/) DO 890 I=1,21 WRITE(30,880)(Z(I)),(RAD(I,J),J=28,36) 880 FORMAT(X,F6.2,X,'M',9F12.1) 890 CONTINUE WRITE(30,959)A 959 FORMAT(/,5X,'ABSORPTION COEFFICIENT',F6.3) SUN=SUN*90.0/1.5708 WRITE(30,960)SUN 960 FORMAT(X,'SOLAR ELEVATION = ',F6.2) WRITE(30,970)IJK 970 FORMAT(5X,'II= ',I5) C THE ANGULAR DISTRIBUTION OF UPWARD SCATTERING EVENTS C WITHIN EACH DEPTH INTERVAL IS NOW CALCULATED WRITE(30,1000) 1000 FORMAT(/,X,'ANGULAR DISTRIBUTION OF UPWARD SCATTERING * EVENTS IN EACH DEPTH INTERVAL') WRITE(30,1655)(ANGLE(I),I=19,27) 1655 FORMAT(/,X,'DEGREES',2X,9F12.1,/) DO 1670 I=1,21 WRITE(30,1660)(Z(I)),(UPRAD(I,J),J=19,27) 1660 FORMAT(X,F6.2,X,'M',X,9F12.0) 1670 CONTINUE WRITE(30,1675)(ANGLE(I),I=28,36) 1675 FORMAT(/,X,'DEGREES',2X,9F12.1,/) DO 1690 I=1,21 WRITE(30,1680)(Z(I)),(UPRAD(I,J),J=28,36) 1680 FORMAT(X,F6.2,X,'M',9F12.0) 1690 CONTINUE C THE AVERAGE (UPWARD) COSINE IN EACH DEPTH INTERVAL C OF THE PHOTONS WHICH HAVE JUST BEEN SCATTERED C UPWARDS IS NOW CALCULATED WRITE(30,1800) 1800 FORMAT(/,X,'AVERAGE UPWARD COSINE OF SCATTERING * EVENTS',/,X,'DEPTH',3X,'AV. COSINE',3X,'NUMBER OF EVENTS') DO 1810 I=1,21 DO 1810 J=19,36 1810 SCUP(I)=SCUP(I)+UPRAD(I,J)/SINANG(J) C The equivalent scalar irradiance due to photons scattered C upwards within each layer has been calculated. DO 1840 I=1,21 UPCOS(I)=UP(I)/SCUP(I) C The average cosine of the upward scattering events is C obtained by dividing the (equivalent) irradiance by C the (equivalent) scalar irradiance. WRITE(30,1830)(Z(I)),(UPCOS(I)),(UP(I)) 1830 FORMAT(X,F6.2,X,'M',F10.3,5X,F10.0) 1840 CONTINUE END FUNCTION IRAN(MM,M1,IB,IR) IA=IR IA1=IA/M1 IA0=MOD(IA,M1) IB1=IB/M1 IB0=MOD(IB,M1) MX=IA0*IB1+IA1*IB0 MY=MOD(MX,M1)*M1+IA0*IB0 MULT=MOD(MY,MM) IRAN=MOD((MULT+1),MM) END 0.044 0.131 0.216 0.301 0.383 0.462 0.537 0.609 0.676 0.737 0.793 0.843 0.887 0.924 0.954 0.976 0.991 0.999 0.999 0.991 0.976 0.954 0.924 0.887 0.843 0.793 0.737 0.676 0.609 0.537 0.462 0.383 0.301 0.216 0.131 0.044 0.45798 0.59525 0.68686 0.75292 0.79888 0.83190 0.85741 0.87756 0.89360 0.90638 0.91672 0.92511 0.93220 0.93822 0.94350 0.94809 0.95226 0.95599 0.95952 0.96280 0.96603 0.96913 0.97221 0.97518 0.97822 0.98118 0.98407 0.98672 0.98928 0.99160 0.99380 0.99570 0.99743 0.99877 0.99972 1.0000 The sin(theta)/f(theta) data are for equal (5 degree) intervals. This scattering phase function is derived from the Petzold (1972) data set for Bahamas Islands, Station 8. 0.044 0.131 0.216 0.301 0.383 0.462 0.537 0.609 0.676 0.737 0.793 0.843 0.887 0.924 0.954 0.976 0.991 0.999 0.999 0.991 0.976 0.954 0.924 0.887 0.843 0.793 0.737 0.676 0.609 0.537 0.462 0.383 0.301 0.216 0.131 0.044 0.517 0.666 0.751 0.811 0.851 0.879 0.901 0.918 0.931 0.942 0.950 0.957 0.963 0.968 0.972 0.975 0.978 0.981 0.9832 0.9853 0.9870 0.9888 0.9902 0.9917 0.9928 0.9941 0.9950 0.9960 0.9968 0.9976 0.9982 0.9988 0.9992 0.9997 0.9998 1.0000 Data file Sandiego.DAT contains the sin(theta) and F(theta) data for San Diego harbour (Petzold). The sin(theta)/f(theta) data are for equal (5 degree) intervals. 0.3 60.0 1.0 10000000.0 25.0 20000 5.0 0.2 0.02 0.2 C Data file input.dat. Provides input data for program Oceanspec C 1. particle scattering coefficient, m^-1 C 2. solar altitude C 3. depth interval, m C 4. Number of photons C 5. Depth of mixed layer C 6. Integer to start random number generator C 7. Windspeed (m/s) at 41 ft (12.5 m) above water C 8. Absorption coefficient (m^-1) of CDOM at 440 nm C 9. Absorption coefficient (m^-1) of IPM at 440 nm. C 10. Phytoplankton chl conc., mg chl a m^-3. 1 20 60.0 1010.0 10.0 1.3 5 20.0 2.5 0.250 C Data file solar.dat provides input data for Fortran C program solirrad.for C 1. month C 2. day of the month C 3. solar altitude, degrees C 4. atmospheric pressure, mb C 5. visual range, km C 6. Angstrom exponent C 7. Air mass type C 8. relative humidity, % C 9. precipitable water, cm C 10. ozone, atm-cm (Dobson Units/1000) 0.044,0.131,0.216,0.301,0.383,0.462,0.537,0.609,0.676, 0.737,0.793,0.843,0.887,0.924,0.954,0.976,0.991,0.999, 0.999,0.991,0.976,0.954,0.924,0.887,0.843,0.793,0.737, 0.676,0.609,0.537,0.462,0.383,0.301,0.216,0.131,0.044, 0.517,0.666,0.751,0.811,0.851,0.879,0.901,0.918,0.931, 0.942,0.950,0.957,0.963,0.968,0.972,0.975,0.978,0.981, 0.9832,0.9853,0.9870,0.9888,0.9902,0.9917,0.9928, 0.9941,0.9950,0.9960,0.9968,0.9976,0.9982,0.9988, 0.9992,0.9997,0.9998,1.0000, Data file Sandiego.DAT contains the sin(theta) and F(theta) data for San Diego harbour (Petzold). The sin(theta)/f(theta) data are for equal (5 degree) intervals.