, ''11 University of Minnesota St. Anthony Falls Hydraulic Laboratory Project Report No. 263 DYNAMIC LAKE WATER QUALITY SIMULATION MODEL "MINLAKE" by Michael J. Riley and Heinz G. Stefan Prepared for LEGISLATIVE COMMISSION ON MINNESOTA RESOURCES State of Minnesota St. Paul, Minnesota August, 1987 The University of Minnesota is committed to the policy that all persons shall have equal access to its programs, facilities, and employment without regard to race, creed, color, sex, national origin, or handicap. " ,.' II. III. IV. V. TABLE OF CONTENTS Page No. Li at of Figures .•.•......•. "." ............................ " 11 PURPOSE ................ " .......... , .................... . 1 MODEL OVERVIEW ••••••••••••••••••••••••••••••••••••••••• A. Components •• " " ••• " • , •• , , , , •••••• " •••••••••.•• " • , •• • 3 B. Modelling Option~ •••••••••••••••••••••••••••••••••• 6 GOVERNING EQUATIONS . ...... , . , , ..... , ..................... ., . METHODS OF SOLUTION , , .................................. . THE STATE VARIABLES . ............................... , ...... . 1. Temperature ............ , ••• ., •••••••• \I •••••••••••••••• 2. Phytoplankton ................ , •• ~ ••••••.••••••.. • • f • .. :3 • Nutrients •••••••••••••••..•.• , ...................... . 4. Dissolved Oxygen and Detritus (BOD) •••••••••••••••• 5. Zooplankton ••••••••.••••••••••••••••••.•••••••••••• 8 11 12 12 20 35 41 43 VI. MODEL OUTPUT ••••••.•••••••• , •.•••• , •••••••••••.•••••••. 48 VII. c, MODEL OPERATION ., •••• \11.' •• • •• 111 ••••••••• ,., ••••••••••••• REFERENCES ............ , .............. ., ...... , ... ., ..... , APPENDIX A Input Data Requirements ••••••••••••••••••• APPENDIX B - Discretization and Solution Procedures for Diffusion and Advection-Diffusion 49 51 55 Equation ••••••••••• ~ , ••. ., •• '. • • • • • • • • • • • • • • 63 APPENDIX C Tabular Output Data Format •••••••••••••••• 67 APPENDIX D Graphical Output Pormata •••••••••••••••••• APPENDIX E Annotated Program Listing ••••••••••••••••• 77 APPENDIX F Definitions of Variables in Common Blocks • 127 APPENDIX G List of Model Subroutines ••••••••••••••••• 13 ] Lake Specific Subroutine Description •••••• APPENDIX I - Example Lake Specific Subroutine •••••••••• 137 APPENDIX H 133 i Figure No. II-l II-2 V-l V-2 V-3 V-4 V-5 V-6 V-7 V-8 V-9 V-lO LIST OF FIGURES A schematic representation of some physical processes occurring within a lake. Schematic diagram of program. Schematic diagram of the source and sink terms treated in the heat budget submodel. Schematic representation of the differential equation for phytoplankton dynamics. Graphical representation of the Michaelis-Mention growth function. The half-saturation constant (k) corresponds to the nutrient concentration for which the growth rate is half the maximum growth rate. The maximum concentration is approached asymptotically with increasing nutrient concentration. Light inhibition function indicating regions of light limited growth and light inhibited growth. Graphical representation of non-linear growth and linear uptake as a function of the internal nutrient concentration. Graphical representation of the temperature effect on growth for a species with an optimum growth rate at 28°C (Topt) and a 90 percent inhibition at 33°C(Tmax). The low temperature ninety percent inhibition is taken to be at O°C (Tmin). Schematic illustration of the processes affecting the predicted chlorophyll-a concentration in each level of the biological submodel. Differential equations for chlorophyll-a concentration in the different levels of chlorophyll-a simulation. Differential equation for processes affecting the concentration of biological available phosphorus. Illustration of the processes (arrows) and components (circles) comprising the nitrogen cycle as treated in the MINLAKE Water Quality Program .• w Figure No. V-ll Differential equation for processes affecting the ammonium concentration. V-12 V-13 Differential equation for processes affecting the concentration of nitrate-nitrite. Coupled Differential equations for BOD and Dissolved O~ygen. iii I. PURPOSE Lakes are in a continuous state of change. Some changes occur over very short periods of time on the order of hours and minutes. Other changes mark long-term trends in the biological and physical condition of a lake. Among the long-term trends in many lakes are a decrease in depth and an increasing productivity in a process called eutrophication. Typically, eutrophication occurs over a time scale of centuries. However, man-made changes to the watershed of a lake may result in a rapid accelera- tion of eutrophication such that significant changes in the water quality of the lake are noticed in a time span of a few years. Anthropogenic acce- leration of eutrophication, known as cultural eutrophication, is due to agricultural, urban, and recreational development in the watershed of a lake which causes an increase in the nutrient loading to a lake. In response to the concern for the water quality of an eutrophic lake, many lake treatment practices have been developed. Treatment processes can be divided into watershed practice methods, inflow-outflow methods, and in- lake treatment methods. A partial list of these practices is provided in Table I. In some instances, treatment methods have been implemented without significantly improving the quality of a lake. The failure of a treatment method is often a result of a poor understanding of the short- term dynamics of a lake which affect the rate and extent of cycling of nutrients within the water column, sediment, and biological components of a lake. To evaluate the effectiveness of lake treatment methods on Minnesota lakes, the Minnesota Lake Model (MINLAKE) has been proposed to model the principal dynamic relationships in a lake on a daily time scale. Treatment methods can then be modeled by modifying those components of the model which will be directly affected by the treatment process. The model can then predict the indirect and direct effects over a period of several months. TABLE I. PARTIAL LIST OF LAKE TREATMENT ALTERNATIVES In-Flow Controls Iloflltrltlon Ind Uptlke of Hatrlent. 1. Dete"tion poad. uUUzla .... eroph,te. to rellOve Dutrlenu 2. Yetl.nd filtratloa of Ilutrleat. 1. Seep •• e treaehe. 4. .u ... lnua oxide colWln. 5. .u ... lnua .1ud.e 6. 'errlc .1uai" ... block. 7. 'errlc lroa treltaent pilat Sedl .. nt IellOvel 8. Sedl.ent retentloa be.la. 9. 'erl.eter road .ediment trIp. Storwvater 10. .uualaU8 elud.. tr •• t .. nt 11. Deteatloa aad .tora.e lla. Coll.ctloa _,.te. eO"trol 12. Seva •• ov.rflov ellataatioa 13. Under.round flltratioa 14. Wetland filtration Wa.tavat.r 15. CeDtral .ever ••• 16. ZVapotr.n.pir.tloa bad 17. Gr., v.ur _a.le .. nt 11. Mouad .,.t ... 19. S.nd fUt.r ZOo S.ptic t.nk .,.t ... 21. t.rtiar, tre.teeat 22. Wa.tevater modiflc.tloa Waterahed Aarieulture1 23. AD ... 1 v •• t. coataineent 24. AD ... 1 •• t. _nal ..... t 25. Wa.te ItoU.. poad 26. ".at. tre.tee.t lapona 27. 'iltar Itrip. .l!!!!!!. 21. ProdUct .ad u •• IIOdlfic.tion to reduce autrienta 29. Strllt sv.eplal 10. Ule of poraul COBCrete tak. W.terahed 11. Di .. rllo. lZ. Gre.abelca 11. 'ari .. ter road .e.ieeat trIp. 14. ..11' •• 1 lboI'. baakl 15. SlIor. Iro.io. protecUoa Ind lubUhatioa Re .. l.tory Coatrol 16. laatina rlsul.tlona 17. Con.lrvltlo ....... nt. 11. Co.trol of rlsul.tloa reaoval 19. 11'01100 coatrol ordl •• acl 40. 1.,1'0" I.forc ... nc of ettltlna r ... latloa 41. Populltio. d.a.lt, control 42. aanit.r, codel 43. 5II0r.llad prot.ctioa ordlnancl 44. Subdlvilioa relu1.tlon. 45. Wlt1and con •• rv.nc, zonlnl 46. Zonlnl cantrall 2 In-Lake Controls take Deep.nlftl S.dl.eat lesoval 47. Dred.lna 48. lac .. etloa la Situ conloUd.Uon - ph,etcal 49. DravdoVD and de.lceatlon 50. aquifer dev.terlns 1. Situ coaaoUdltloa - dl.uUoa 51. Aeration 5Z. IIltro.e •• ddltioa Sl. H,dro.e. peroxUe .ddlt Ion 54. 0 zoal ••• 1 Uo. 55. Proprlet.", or.aal .. addltio. Se41 .. nt lIutrleat In.cUvltioa ClInical 56. Hypollan.tlc aerltloa 57. H,pollanetlc OllJ •• netioa 51. &1U11laillll lajlction lata .ed taeat 59. &1 ... lal... floc 11,erinl 60. '1, I.h l.,erlftS 61. W.ter tre.t •• at pla.t .1udlle Ph,lic.l .e.U.1 62. Plaeelc. via,l .a. I,athetic 11nerl Il. Cl.rMontllOrllUt. la,erilll 64. S.nd la,erial W.t.r Co1WID Mutrl.at Iaactlvatloa 65. AtUIIlnl .. tre.teeat 66. Pot ••• I ... al ... lnlua aulfate 67. 'errie chlorld. 61. Zireonl ... 69. Lanthen ... 70. Air 11. Water jec .utrlent Oucflow Aceeteracion 7Z. Salecclve discharlle 73. ,1uahial/.ilution 74. H,poUaneUc vithdraval II0101ic.l Probl.. tre.ceent Phytoplanktoa 75. Copper .uU.te 76. Cupro •• 77. Cutr1ne Macrophyte. 78. Pl.at hec" .. Ual 79. Macroph,t. herbicide. 80. Bottos barrlar. 81. Shad .. 12. 1If •• 83. 'Iah Predatlam 16. Sh.Ufi.h 85. In •• ct. 86. Dl..... orsaal ... 87. CoIlp.UU" pleat. 88. .,.raollc dredSlal 89. Dlv.r dre •• ln. 90. Manu.l harve.tlnl 91. U.. of a roto".tor 92. Shalla. veter tiliale 93. Dravdova 94. 'hh ". PhdcUu I" " II. MODEL OVERVIEW Figure 11-1 shows some of the components and dynamic relationships in an actual lake that are simulated in the model. To incorporate all processes occurring in a lake would result in a model that would be too difficult to formulate and probably too expensive to operate on a large number of lakes. Consequently, the various components have been evaluated and only the most significant included in the model. The following sections describe the components included in the model. A. Components The structure of the MNLAKE program is built on five main sections: input, heat budget, biological-nutrient, inflow-outflow, and a lake specific subroutine. A schematic diagram of the program showing the main components of each section is given in Fig. 11-2. The input routine provides the information necessary to begin the computation sequence. The first of the input routines (START) reads the initial conditions which are the depth profiles of each state variable for the first day of computation. The second input sequence provides the rate coefficients and yield coefficients used to calibrate the model to a par- ticular lake. A list of calibration variables is given in Appendix A. Finally, the subroutine MTHDATA provides the meteorological data for each day of the month. The meteorological variables used in the program are listed in Appendix A. The initial conditions and the calibration variables are read only once while the MTHDATA subroutine is called at the 'beginning of each month of simulation. The heat budget section is the first of the modelling sequences. In the heat budget section, meteorological data is used to compute the influx and efflux of heat from the water body and the extent of wind mixing. The heat balance and the wind mixing determine the temperature profile and the depth of the mixed layer. The heat budget routines also simulate convec- tive mixing which occurs when a colder, more dense water layer forms over a warmer water layer. Convective mixing is important during cool weather periods when the surface water cools more rapidly than water below the surface. The inflow-outflow routines simulate the variations in stage through the season. Inflow and outflow can be either input parameters or computed internally. The inflow must be in terms of water volume, the temperature of the inflow, and the concentration of all other state variables. Consequently, if the inflow is read from another file, the input data becomes considerable especially if there is more than one source of inflow. Ideally, a separate hydrologic model would be available to generate the inflow quantity and quality from each source based on meteorological conditions. A known outflow is more easily handled as only the quantity is required. The outflow removes a volume of water without changing the 3 ~ [WINO- MIXING I I .. ENTRAINMENT INTERFLOW SURFACE HEAT EXCHANGE MIXED LAYER THERMOCLINE EUPHOTIC ZONE SEDIMENTS I • ,-I INTER LAYER DIFFUSION I • ,-T I • fr Fig. II-I. A schematic representation of some physical processes occurring within a lake. " \ .. " INPUT HYDROLOGIC MODEL , . ~ INFLOW QUALITY- INITIAL' QUANTITY CONDITIONS DATA HYDRODYNAMIC PROCESSES HEAT WATER BALANCE TEMPERATURE MIXING CONSERVATION OF MASS MODEL WATER DENSITY CURRENT BALANCE MODEL BIOCHEMICAL COMPONENTS NON-LIMITING NUTRIENTS LIMITING· NUTRIENTS INTERNAL SOURCES AND SINKS OUTPUT SUSPENDED SEDIMENT I PREDICTED LAKE QUALITY I METEOROLOGICAL MORPHOMETRY DATA WIND MIXING CONVECTIVE MODEL MIXING MODEL OUTFLOW GROUNDWATER MODEL MODEL I PREDICTED OUTFLOW QUALITY I Figure 1I-2. Schematic diagram of program. 5 concentration of the state variables. It is simply read in as a negative inflow with zero for the concentration of each of the state variables. The inflow-outflow sequence also involves computations to determine that the mass of water and each state variable is conserved. The inflow can be to any layer depending on the density of the inflowing water. Consequently, the volume and concentration of the state variables in the receiving layer must be adjusted based on the ambient conditions and the volume and concentrations of the water added. The new volume in the layer is simply the sum of the initial layer volume and the volume of inflow, including water entrained from previous layers. The new concentration of each state variable in the layer is the volume weighted average of the ini- tial concentration in the layer and the concentration in the inflow. Included in the inflow-outflow sequence is a layer thickness check. Outflow from a layer could reduce a layer to a very thin thickness or low volume. Likewise, inflow volume could result in a layer becoming very thick. Very thin layers are not realistic since these would be easily entrained by adjacent layers. Very thick layers near the surface can cause problems in the prediction of light attenuation and photosynthesis. In addition, the development of many thin layers could result in the number of layers exceeding the maximum allowed by the program. Consequently, the final operation of the inflow-outflow sequence is to merge thin layers with the next deeper layer and to divide thick layers into two or more layers. The maximum and minimum thickness are input variables specified by the user. The biological-nutrient routines are the most complex part of the program. The complexity is due to the interrelationships among the many state variables. The biology is characterized by one to three phytoplank- ton groups and one zooplankton group. The growth of the phytoplankton in the model is potentially limited by phosphorus, nitrogen as either nitrate-nitrite or ammonium, or light. Phosphorus, nitrate-nitrite, and ammonium are state variables. Light penetration is computed based on incoming solar radiation and attentuation of light with depth. In addi- tion, detritus is modeled as a pool of oxygen consuming, decaying organic matter (BOD) which is linked to the dissolved oxygen concentration and releases nutrients during the decay process. The lake specific subroutine allows the model to be configured to one particular lake by inserting equations for depth-area and depth-volume computations, fetch, and special features such as non-point inflows or groundwater flow. The subroutine also allows the modeler to modify or add source and sink terms to the biological-nutrient equations and to add submode1s to simulate different treatment options. B. Modeling Options The MNLAKE model offers two modes of operation in modeling a lake. The least complex model is for temperature only. The purpose of a tem- perature only model is to allow the temperature dynamics to be calibrated first at a lower cost than would be possible for computations involving all the state variables. While the biological and chemical reactions of the biological-nutrient section are functions of temperature, the temperature 6 \If " is only slightly affected by the biological state variable. Therefore the heat budget routines can be applied independently of the rest of the model with little loss in accuracy. The second level of complexity is in the treatment of algae growth and chlorophyll concentration. Chlorophyll can be modeled by either a very simple relationship for chlorophyll in the euphotic zone or .by a more detailed layer by layer treatment using either a Michaelis-Mention (Monod) model or an internal nutrient (cell quota) model. The more detailed and costly models are only justified when adequate in-lake data is available for calibration of the model. For rough estimates of the dynamics of a lake with little supporting data, the less costly chlorophyll model can be used. All of the chlorophyll models use both phosphorus and light as fac- tors contributing to the limitation of algal growth. In the more complex internal nutrient model, nitrogen can be included as a.limiting agent. The nitrogen cycle is simulated with algae utilizing nitrate-nitrite and ammo- nium nitrogen for growth. The nitrogen cycle adds two state variables to the model. Yet, many lakes are considered to be phosphorus limited and not nitrogen limited. Therefore, nitrogen was included as an option to reduce the cost of the model on lakes that are strongly phosphorus limited. In this section a general discussion of the various state variables. with respect to the structure of the program and the sequence of com- putations has been given. A more detailed description of the state variables, governing equations, and methods of solution will be presented in the next section. 7 III. GOVERNING EQUATIONS The dynamics of all substances (state variables) characterizing a water volume whether dissolved or suspended can be expressed by the same one-dimensional advection-diffusion equation: where A ac + v aCA = ~ KA ac + sources/sinks at az az az- -1 C = concentration (mg t ) or intrinsic property of the fluid (such as heat content) v = vertical velocity component of the substance (m/day) = 0 for dissolved substances and temperature z = depth (m) K = turbulent diffusion coefficient (assumed to be identical for each state variable) (m2 /day) A = surface area of the control Iolume (not constant with depth) (m ) In words, the advection-diffusion equation states that the time rate of change of a variable (3c/at) plus the rate of inflow of the variable V(aAc/az) equals the rate of diffusion of the variable out of the control volume [(a/3z)kA(3c/3z)] plus the addition or loss due to sources or sinks within the control volume. In the one-dimensional advection-diffusion equation, the velocity component and the turbulent diffusion coefficient must be ,either known, estimated, or computed. The diffusion coefficient is computed separately in the epilimnion and hypolimnion and assumed to be the same for all substances whether dissolved or suspended. The diffusion coefficient for control volumes in the epilimnion is given as a function of wind speed only. The relationship is given by K = 28 w1.3 where w = wind speed (mph). The formulation yields a high diffusion rate that results in,a nearly uniform distribution of dissolved and suspended substances within the mixed layer. The hypolimnetic diffusion coefficient formulation follows earlier work by Mortimer (1942), Welander (1967), and Jasseby and Powell (1975) by 8 .. II relating the diffusion coefficient to the Brunt-Vasa1a frequency (N): .. K ;;. N-1 in which N ;:: where g gravitational acceleration (m/s2) p ;=: density z ;:: depth (m) The Brunt-Vasala frequency for each control volume is computed using a finite difference approximation for the differential term. The K-N relationship is characterized by a maximum hypo1imnetic diffusion rate which is assumed to be related to the surface area of the lake [Mortimer, 1942]. The magnitude of the density differential between control volumes reduces the diffusion coefficient. where K=K max C N-1 k . K = maximum hypolimnetic diffusion coefficient estimated from max 2 surface area of the lake (m /day) Ck cO,efficient. Since the value of K must not be greater than Kmax' the computation for . the hypolimnetic diffusion coefficient is then given by -1 K =' min [K t K . Ck N ] max max in which min denotes minimum and ooefficient Ck is considered as the minimum value of N at which the maximum hypolimnetic diffusion .rate occurs (taken to be 8.66xlO-3 from Jassby and Powell, 1~n5). The vertical velocity component is either computed or estimated and is given as positive in the downward direction. The. suspended sediment settling rate is computed from the formulation given by Gibbs et al. [1971]. -3d + ... /9d2 + gi Pf{p - Pf )(0.015476 + 0.1984lr ) v = __ ~v~s ____ l__ ~v~s __ ~~s~~~s=-__ ~~~ _____ ~ ________ '___ s~ Pf(0.011607 + 0.14881 rs) * 864 9 where d dynamic viscosity (poises) vs -3 1. 79(10 )p f 1 + 0.03368T + 0.000221T2 T water temperature (OC) r s sediment particle radius (cm) density of fluid (g/cm3 ) 3 density of sediment particle (g/cm ) gravitational acceleration (cm/s2) Most of the parameters in the Gibbs equation are nearly constant. Assuming a water density of 1000 kg/m3 , relative specific weight of sediment of 2.65, and a particle diameter of .5 ~m, the equation is simplified to: and v d vs (-3d + .. /9 d2 + .1565xlO-7 vs l' vs 1. 79 (1 + 0.03368T + .000221~) (74390) The settling velocity of the other suspended state variables (chlorophyll-a and detritus) are estimated and must be specified by the model user. The advection-diffusion equation is made particular to one state variable by the choice of source and sink terms. The sources and sinks determine the complexity of the model and the interrelationship among the state variables. Two state variables that are related will have at least one complimentary source and sink where the same expression is a source or sink in the equation of the first state variable and a source or sink in the equation of the second state variable. Directly related variables include detritus as BOD and dissolved oxygen, phosphorus and biomass as chlorophyll. Most of the state variables are also indirectly related to temperature. The temperature in a layer affects the rate of biological reaction such as algal growth rate and respiration rate. However, the effect is weaker than a direct relationship and the biological reactions do not become a source or sink in the temperature equation. 10 .. IV. METHODS OF SOLUTION Two methods of solution are used in the model depending on whether or not the state variable has a component of velocity in the vertical direc- tion. Both solutions use a fully implicit solution for a control volume. The dissolved substances and temperature use a central difference scheme to discretize the partial differential equations. The discretization of the equation is given in Appendix B. The suspended state variables (chlorophyll, suspended sediment, and BOD) cannot be modeled by a simple central difference scheme due to the dual effect of advection and diffusion. The discretization is taken from the suspended sediment model developed by Dhamotharan et al. (1981) and uses the implicit power-law scheme developed by Patankar (1979). The discretization of the full advection-diffusion equation is presented in Appendix B. The discretization of each equation yields a set of linear equations, one for each layer, which must be solved simultaneously. Because the con- centration of a variable within a layer is only dependent on the con- centration in the two adjacent layers, the set of equations can be assembled into a tri-diagona1 matrix. A modified Gauss elimination proce- dure is used to solve the tri-diagona1 matrix operating only on the three non-zero diagonal terms. The solution algorithm is presented in Appendix B. For each day of simulation, the model solves each state variable in sequence beginning with temperature and wind mixing, followed by suspended and dissolved solids, nutrients, chlorophyll-a, detritus, and dissolved oxygen. In actuality, all the processes are occurring simultaneously and are approximated here in sequential order. All state variables could be handled simultaneously in the solution of a single large matrix. However, the matrix would be sparse, and the non-linear relationships among variables would require an iterative solution. The iterative solution of a sparse matrix would be excessively expensive considering that the errors introduced by the sequential formulation would be small and the errors do not propogate. For example, if a nutrient is slightly over-estimated one day, it will contribute to increased growth, which will in turn cause a decrease in nutrients the next day. Temperature is solved first because it affects most other state variables by affecting reaction rates, but temperature is only weakly affected by other state variables. The strongest interrelationship among variables is in the nutrient and chlorophyll-a equations. If a problem does develop, the nutrient and chlorophyll-a solutions may be done through iterations, but this will also cause some processes such as sediment nutrient exchange to be simulated twice in a time step, and the reaction rate coefficients will have to be adjusted accordingly. 11 V. THE STATE VARIABLES The "state variables" are quantities that vary with time and are used to characterize the state of the lake at any time. In the current model, the state variables are temperature, suspended solids, biomass as chlorophyll-a, available phosphorus, dissolved oxygen, and detritus as BOD. In addition, nitrate-nitrite, ammonium, and internal cellular nutrient levels are optional state variables. Following is a brief discussion of each state variable and the sources and sinks treated in the model. 1. Temperature The depth profile of temperature is computed from a balance between incoming heat from solar and long-wave radiation and the outflow of heat through convection, evaporation, and back radiation. The net increase in heat results in an increase in water temperature. The heat balance equation is given by AH .. H + H - H - H - ll. sac e -0 where the subscripts denote solar radiation (s), long-wave radiation from the atmostphere (a), convection (c), evaporation (e), and back radiation (b). A schematic diagram of the heat flux terms is presented in Figure V-1. Each term must be replaced by a function relating the heat source or sink to known parameters such as the meteorological conditions. The functions presented here are taken from earlier work by Stefan and Ford [1975], Ford [1976], and Dhamotharan [1979]. Solar radiation is a source of heat in the upper layers of the lake. Radiation reaching the surface is either reflected or penetrates the lake surface. The reflected fraction (r), is computed as a function of the solar radiation and the concentration of suspended sediment in the surface layer [Dhamotharan, 1979; Stefan et a1., 1982]. -5 r = 0.087 - 6.76(10 )*RAD + 0.11[1 - exp(-0.01*SS)] where r .. reflected fraction RAD -2 -1 - total daily solar radiation (cal cm day ) S8 .. suspended sediment concentration in the first layer (mg/t) The relationship was developed for a highly turbid southern lake. In most Minnesota lakes, the suspended sediments will be low and the sediment ]2 '< \.I) Atmospheric Radiation Incoming Solar Radiation Reflected /Fraction -, Evaporation ~ I .. Surface Absorption I I • Absorption with Depth I I Fraction Transmitted , to Next Layer I I • Absorption with Depth I I Fraction Transmitted , to Next Layer I I I Convection ~ Back Radiation Fig. V-I. Schematic diagram of the source and sink terms treated in the heat budget submodel. ~ • effect on reflectance can be dropped. The reflectance equation is not general and more work needs to be done to verify that the equation is appropriate for Minnesota lakes. A sizeable portion of the penetrating solar radiation is absorbed within the first few centimeters of the lake surface. This is due to the high attenuation rate for infared radiation which accounts for approxi- mately 40 percent of the energy in solar radiation. Consequently, the sur- face absorption fraction is taken to be a constant of 0.40 [Dake and Harleman, 1969]. The net solar radiation below the surface of the lake is then given by Hsn = (1 - r)(1 - S)Hs where H sn net solar radiation (kcal m-2 day-I) r = reflected fraction S = surface absorption fraction = 0.40 -2 -1 H = incoming solar radiation (kcal m day ) s In each layer of the lake additional solar radiation is absorbed until all is absorbed. Consequently, solar radiation is a heat source in more than just the top layer of the lake. The amount of heat absorbed in a layer is the difference between the solar radiation at the top and the bot- tom of the layer. The attenuation of solar radiation with depth follows the well~known Beers law for the attenuation of light. HS(i) = Hs(i-l) * exp(-kz) where -2 -1 Hs (i-l) - solar radiation at the top of the layer (kcal m day ) -2 -1 HS(i) = solar radiation at the bottom of the layer (kcal m day ) k = extinction coefficient (m-1) z = thickness of the layer (m) The total absorbed solar radiation in the surface layer (kcal) becomes Ql - (l-r)(l-B)*H * [AI - A2 exp(-k zl)] + SH A sn sn ] ]4 i r where 2 lake surface area (m ) surface area of the top of the second layer (m2 ) thickness of the first layer (m) The solar radiation reaching the top of the second layer is and the absorbed solar radiation in the second layer becomes The attenuation-adsorption sequence continues until essentially all the solar radiation is absorbed. The extinction coefficient is made unique to each layer as a linear function of three terms; the extinction coefficient of water, the extinc- tion coefficient of the suspended sediment concentration in the layer, and the extinction coefficient of the chlorophyll-a concentration in the layer. The total extinction coefficient is given as k = kw + 0.043*88 + k2*Chla where k = total attenuation coefficient (m -1) k w k2 88 Chla extinction coefficient of the water (m-1) extinction coefficient due to chlorophyll -1 suspended solids concentration (mg ~ ) -1 chlorophyll-a concentration (mg ~ ) -1 (t(m mg chla) ) The chlorophyll coefficient (k ) of 0.016 is taken from Bannister 11974a] but may range from 0.10 to 0.25. The suspended sediment coef- ficient is taken from Dhamotharan [1979J and actually must represent the combined effect of suspended sediment, detritus, and non-chlorophyll con- taining portions of algal cells. k is not the same as the extinction coefficient for pure water, but includes the combined effect of all dissolved substances especially color caused by dissolved organics. Consequently, k can be expected to be specific to a particular lake. w Long-wave radiation is energy emitted by all objects. The heat budget must account for both atmospheric long wave radiation received by the lake and long wave radiation emitted by the lake. Long wave radiation is given by 15 where a = Stefan-Boltzman constant £ = emissivity T temperature of the object (OK) For an ideal black body the emissivity is unity. Neither the water nor the atmosphere are ideal black boxes, but water does have an approxima- tely constant emissivity of 0.97. The emissivity of the atmosphere is more variable being a function of cloud cover and the temperature of the atmosphere. The model uses a functional relationship between emissivity and temperature developed by Idso and Jackson [1969] with a cloud cover correction where £ a T = atmospheric temperature (OC) a C = cloud cover ratio c K = cloud height coefficient The coefficient K varies between 0.04 and 0.25. A TVA [1968] study recommended using a value of 0.17. The atmospheric radiation uses the average daily temperature from the meteorological record. Back radiation is computed using the temperature of the mixed layer. Evaporative heat transfer is due to the change in state of water from a liquid to a vapor in the evaporation process. The evaporative heat loss can be computed directly if the mass of water evaporated is known. where H = k*E e E = total evaporation (cm) k = latent heat of vaporization 5900 kcal cm-l m-2 day-l In most cases the actual mass lost to evaporation is not known and the eva- porative heat flux must be approximated. The standard empirical for- mulation for evaporative heat loss is given by H kW(e - e ) e s a 16 .. where k - latent heat of vaporization per unit mass W • a wind function e - saturated vapor pressure at the water surface s e - vapor pressure of the air a The saturated vapor pressure at the water surface is computed using the surface temperature of the lake in the Magnus-Tetons formula [see Murray, 1967]. e - 6.1078*exp[17.27*T I(T k - 35.86») s s s where T - surface water temperature (OC) s Tsk = surface water temperature (OK) The actual vapor pressure of the air is computed from the dew point temperature. e -a The wind function used was developed by the U. S. Geological Survey from data taken on Lake Hefner and Lake Mead [Dake, 1972]. The U. S. Geological Survey equation was adapted by Ford [1976] following comparisons of results of several equations in a lake stratification model. The equation is a simple linear relationship W • WCoef u where u - the average daily wind speed in mls WCoef = wind function coefficient (taken as 25.62) Convective heat transfer occurs across the lake surface due to the difference in temperature between the water and the air. The convective heat loss is given by H = 0.618 W (T - T ) c s a 17 where T water surface temperature (OC) s T = air temperature (OC) a W = wind function which is the same as previously described The model does not simulate any heat flux from the water to the lake sediment. Although some heat transfer occurs, the quantity of heat is small and should result in slightly higher predicted temperatures in the hypolimnion during the summer. An essential feature of the present model is the wind induced mixing of the surface layers into one constant temperature layer. The extent of wind mixing is determined from a stability criteria [Ford, 1976] given by (J = v Apg(z - z ) m m g where At = time interval (one day) A = surface area Pw = density of water u* = wind induced water shear velocity V = volume of the mixed layer ID Ap = density difference between the mixed layer and the layer immediately below the mixed layer z = depth of the mixed layer m z = center of gravity of the mixed layer g The mixing computation proceeds from layer to layer entraining one layer at a time and computing the corresponding V ,Ap, z ,and z m m g until the stability criteria is less than unity. The wind shear velocity is computed from the wind velocity over the water surface. The fetch of the lake in the direction of the wind is used to evaluate the wind speed over the lake which is greater than the wind speed measured over the land. The wind speed is computed using the simi- larity between logrithmic velocity profiles as [Ford, 1976]: 18 i " where ~ - wind speed over land zl - height at which the wind speed is measured (assumed to be 10 m) zOl • surface roughness of the land (taken as 0.01) z02 = surface roughness of the water (taken as 0.0001) zb - equivalent boundary layer over the water The boundary layer over the water is given as as function of the fetch and land surface roughness [Elliot, 1958] Zb ;000 0.86 FetchO•8 34 0.8 ;000. Fetch The shear velocity of the water due to the wind is given by P 0.5 0 5 ( ~) (C)· u = 0.0343 CO•5 u Pw Z w Z w C is an empirical coefficient that is dependent on the wind speed. C isZdetermined from z C •• 0005 (u )0.5 z w C - 0.0026 z Solution Procedure for u < 15 m/s w for u > 15 m/s w The combination of all the heat flux terms results in a single source term with depth from the adsorption of solar radiation. The surface layer, however, is subject to two source terms, two sink terms, and the convective term which can be either a source or a sink. The actual surface temperature is used to compute the evaporative heat loss (when actual eva- poration is not known), the convective heat flux, and the long wave emitted radiation. Initially, it would appear that an implicit solution with linear approximations for the surface temperature in the convection and long wave radiation equations would be possible. However, a straight implicit solution will over estimate the surface temperature since the 19 effect of wind mixing cannot be incorporated into the implicit solution technique. Consequently, the evaporative heat flux and long-wave radiation flux will be over estimated and the convective heat flux will be either under estimated or overestimated depending on whether the air temperature is higher or lower than the water temperature. To handle the dual problems of wind mixing and surface temperature, the heat flux is treated in a sequential manner. The sources of heat (solar radiation and atmospheric radiation) are applied and the temperature of each layer is computed directly for the mixed layer and by the fully implicit central difference scheme below the mixed layer (see Appendix B). The surface temperature is then the same as the temperature in the mixed layer. The loss of heat from the surface through long wave radiation, convection, and evaporation is computed using the mixed layer temperature. The temperature of the entire mixed layer is reduced by the loss of heat from the surface. If the mixed layer temperature is lower than the next deeper layer, a density instability exists. Then, the mixed layer entrains deeper layers until the mixed layer has the lowest density of all the layers in the lake. The effect of wind mixing is applied following the computations for heat flux. The effect of wind action at the surface is to deepen the mixed layer. The extent of mixed layer deepening is determined by computing the kinetic energy imparted to the lake surface by wind. If the energy exceeds the potential energy required to lift the layer below the mixed layer to the center of mass of the mixed layer, then that layer is entrained and the mixed layer depth increases by the thickness of the entrained layer. The process is repeated until the kinetic energy is insufficient to vertically mix the layer below the mixed layer. 2. Phytoplankton The prediction of phytoplankton biomass is the main purpose of a eutrophication model. Biomass can be expressed in a number of different ways including milligrams dry weight of algae per liter, carbon content in milligrams per liter, and chlorophyll-a concentration in milligrams (or micrograms) per liter. There are problems and assumptions required to use any of the expressions for algal concentration. Use of dry weight or carbon content assumes that productivity is directly related to biomass. However, productivity is more closely related to the chlorophyll content of an algal cell. Use of chlorophyll-a as a measure of biomass assumes that chlorophyll is directly related to biomass. In fact, chlorophyll content can fluctuate independently of biomass to meet the phytosynthetic needs of the cell [Schlesinger and Shuter, 1981; Bannister, 1974b]. Because the purpose of the present model is more concerned with productivity than actual biomass, the concentration of chlorophyll-a in a water volume is used to represent an algal population. Fluctuations in phytoplankton populations in a volume of water are due to the combined effects of growth, diffusion, settling, respiration, and grazing. The differential equation for the time rate of change of the phy- toplankton population is given in Figure V-2 as a schematic representation of each term. The descretization of the equation follows the form for a state variable with a settling velocity as given in Appendix B. 20 : r" TIME RATE OF CHANGE OF CHLOROPHYLL CONCENTRATION LOSS TO MORTALITY = NET FLUX INTO LAYER FROM ALGAE SETTLING AND DIFFUSSION LOSS DUE TO RESPIRATION + LOSS DUE TO ZOOPLANKTON GRAZING GROWTH Fig. V-2. Schematic rep~esentation of the differential equation for phytoplankton dynamics. 21 I Three levels of complexity for phytoplankton modelling are provided. The simplest form is a single algal group model with productivity distri- buted throughout the mixed layer. Forsberg and Shapiro (1980) presented a suitable relationship given by R [ In(I /1 ,)P ] ~ k ,Chlj 1 o z max 1 _ q [ ] _ D E Chla + E S ze cwo m where R == -1 specific growth rate (day ) Chla I 0 I z' z' P max concentration of chlorophyll-a -1 :: (mg ChI a t ) :: intensity of photosynthetically active radiation (PAR) just below the surface. = intensity of PAR at depth z, :: empirically defined water depth (m) :: maximum specific daily rate of photosynthesis under nutrient saturated conditions (mg C(mgChla day)-l) z = mixed layer depth (m) m D = specific loss rate due to settling, grazing and -1 metabolism (day ) e = carbon to ch1orophy11-a mass ratio e: c e: w S o -1 attenuation coefficient of PAR by ch1orophyll-a (t(m mg Ch1a) ) attenuation coefficient of PAR by water and substances other than ch10rophyll-a suspended and dissolved in water (m-l ) concentration of limiting nutrient:total phosphorus (mg 1-1) k , = minimum S /Chla ratio for photosynthesis to occur. q 0 Although the list of variables appears to be prodigious most of the variables are known or easily determined. The term In(I II .) is approxi- o z mat~lY constant at 1.9 if the incoming solar radiation exceeds 320 day which it does for most of the summer. Because In(IoIIz~) is known, z' does not have to be determined. So, e: w , e:c ~ a, and zm 'are either input parameters or computed by the model. D, kq ., and P are left as , max calibration parameters. The level-one model attempts to combine all the growth and loss terms into one equation. The form of the equation is not as a differential equation, but an explicit prediction of chlorophyll-a. Consequently, the solution procedure does not follow the form described in Appendix B. The chlorophyll-a concentration is computed directly and applied to all the layers within the mixed layer. 22 .. The level-one model equation is subject to two severe assumptions [Forsberg and Shapiro, 1980]. The phytoplankton population is assumed to be uniformly distributed throughout the mixed layer. Secondly, it is assumed that the light intensity at some point in the mixed layer reaches the level for light saturated photosynthesis to occur. Although there is some doubt as to the validity of both assumptions, the equation is adequate for the purpose of quickly and simply calculating the chlorophyll-a concentration. The second level of complexity is also a single algal group model, but the growth term is unique to each layer. The second level model represents algal growth by a Michaelis-Menton equation [Monod, 1949] subject to both light and phosphorus limitation. The Michaelis-Menton equation for growth limitation is given by where 11 max S -1 = maximum (nutrient saturated) growth rate (day ) concentration of a limiting factor (light or phosphorus) k half-saturation constant s The effect of Michaelis-Menton equation is depicted graphically in Fig. V-3. The meaning of k becomes clear in the diagram as the nutrient s concentration at which the growth rate is half the maximum growth rate. Phosphorus and light limitation are applied on the basis of a single limiting factor. The effect of multiple nutrient limitation is often pre- sented as a multiplicative effect [Lehman et al., 1975; Jorgensen, 1981]. In fact, the effect of multiple limiting nutrients is probably not either singular or multiplicative [Monod, 1949]. However, the multiplicative method results in an extremely severe limitation of growth. In a direct c.omparison of limiting nutrient and multiplicative approac.hes, the single limiting nutrient more closely simulates the fluctuations of algal popu- lations [De Groot, 1983]. Light limitation is simulated by a function to describe light limited growth and inhibition [Megard et al., 19841 usinf the average quantum of photosynthetically active radiation (llE mr2 sec- ) within a layer. The equation for the light limitation and inhibition function, after adjustment to yield values between 0 and 1, is given by where I( 1 + 2/Kl/K2 ) I + Kl + i! /K2 LI = light limited growth ratio (between 0 and 1) 23 -I >-ClS '0 ....... N Q) .I:"- +J ClS I-< ..c:: +J ~ I-< Co!) ].lmax ].lmax 2 0, ~ o ks Nutrient concentration (mg ~-I) Fig. V-3. Graphical representation of the Michaelis-Menton growth function. The half-saturation constant (k )corresponds to the nutrient concentration for which the growth rate is ha~f the maximum growth rate. The maximum concen- tration is approached asymptotically with increasing nutrient concentration. .. ,. -2 -1 I - photosynthetically active irradiance (J.1E m s ) K1,K2 = light limitation and inhibition coefficients, -2 -1 respectively (pE m s ) A typical form of the equation is given in Fig. V-4. The curve yields light limitation for light levels below the optimum and light inhibition for light levels above the optimum. The incoming solar radiation is attenuated in the same manner as described in the heat budget section. The conversion from energy units to quantum units is given by Combs [1977] as where I = 27.25 RAn s TD I - average photosynthetically active radiation (PAR) s -2 -1 over the daylight period (pE m s ) TD • photo period (hr) . -2 -1 RAn .. solar radiation (cal cm day' ) The average quantum in a layer is found by integrating the light attenuation function over the thickness of a layer divided by the thickness [Cardoni and Stefan, 1982] and is expressed as where I i [l - exp(-k*DZi )] k*Dz i I ave ( -2-1 • average PAR in layer i )J E 111 s ) Ii .. k .. PAR reaching the top of the layer -1 extinction coefficient (m . ) = thickness of layer i (m) ( ... 2-1 pE m s ) In the model, the light availability is performed on eight sub- divisions of the photoperiod and using 0.2 meter sublayers of each layer [Cardoni and Stefan, 1982]. The second level of phytoplankton modelling does not incorporate the zooplankton grazing term in Fig. V-2. Instead, the grazing is simulated by use of a high mortality rate if zooplankton grazing is considered to be significant. 25 1.0 Q ~ 0:: Z 0 i=; N ~ 0.5 CJ\ -~ ...J .- :I: C) ::i o LIGHT LIMITED ~ GROWTH .1. LIGHT INHIBITED GROWTH ~ o I opt 5 k2 10 IRRADIANCE (E·m-2 hr- I ) 15 kl Fig. V-4. Light inhibition function indicating regions of light limited growth and light inhibited growth. p, The third level of phytoplankton modelling is considerably more complex than the previous models. Instead of one algal group, three classes of algae are utilized. The groups correspond to diatoms, green algae, and blue-green algae. By modelling the three very rough classes, the level three model can be used to simulate the shift from diatoms early in the season to green and then blue-green algae during the summer. The distinction between algal groups is modelled by different rates of growth, respiration, settling, zooplankton grazing, and different nutrient require- ments. Nutrient limitation and growth is presented as a two-step process in the level three model. Growth is treated as iridependent of the external nutrient level, but is dependent on the available light within a layer or on the internal nutrient content of the algal group. The growth rate is determined by the single limiting nutrient concept with limitation from internal phosphorus or light. Diatoms are also potentially limited by silica. An option is available to treat diatoms and green algae as limited by internal nitrogen. The growth limitation by internal nutrients is given by Lehman et a1. [1975] and Jorgensen [1981]. where N - Chla (N ) min II ... llmax N -1 llmax = nutrient saturated growth rate (day ) -1 N = internal nutrient concentration (mg 1 ) -1 Ch1a - concentration of ch10rophy1l-a (mg 1 ) N - minimum ratio of phosphorus ,to chlorophyll-a min The minimum phosphorus/ch10rophy1l-a ratio is estimated stoichiometri- cally using phosphorus/carbon ratios [Strumm and Morgan, 1981] and assuming a constant carbon/ch10rophy11-a ratio of 50. The actual carbon chlorophyll ratio may vary from 25 to 100 [Schlesinger and Schuter, 1981], depending on the light adaptation and the physical condition of the algae with lower ratios in low light adapted algae and higher ratios in highly stresses algae. By setting a constant minimum phosphorus/ch10rophy11-a ratio, some inaccuracy is introduced at low light levels where growth is also low. At low light levels, the phosphorus/ch10rophyll-a ratio should not be· affecting growth since growth will be light limited. Consequently, the error associated with the assumption of a constant minimum phosphorus/ ch10rophyll-a ratio should be small. The internal nutrient pool is modelled in much the same way as chlorophyll. Respiration, mortality, settling, and grazing reduce the con- centration of nutrients in a volume of water. The internal nutrient con- centration is increased through nutrient uptake from the water. The rate of uptake is affected by both the initial internal nutrient concentration and by the concentration of the nutrient in the water. The uptake rela- tionship is given by Lehman et a1. [1975] and Jorgensen [1981] as 27 where (Chla) (N ) - N S max -1 u = maximum uptake rate (day ) max N = maximum ratio of nutrient to chlorophyll-a max Nmin = minimum ratio Chla = chlorophyll-a of nutrient to chlorophyll-a -1 concentration (mg 1 ) S = nutrient concentration in the water (mg 1-1) -1 K = half-saturation constant for uptake (mg 1 ) s N = internal nutrient concentration (mg 1-1) The first ratio in the equation represents a decrease in the maximum uptake rate due to the nutrient deficit within the algal population. The second ratio is a decrease in the rate of nutrient uptake due to the nutrient concentration in the water. In most phytoplankton models, the growth rate is related only to the external nutrient concentration as in the level two model [DiToro et a1., 1980; Scavia, 1980; Walters, 1980]. The assumption is made that internal nutrient content of the phytoplankton is always in equilibrium with the nutrient concentration in the water column. The assumption is based on the high rate of uptake relative to the growth rate [DiToro, 1980]. The assumption is valid for the case of an increasing nutrient concentration where the rapid uptake causes the increased growth rate to be well corre- lated with the increased nutrient concentration. However, equilibrium models moderate increasing growth. The equilibrium approach most notably breaks down under nutrient depletion where growth continues despite the low nutrient content of the water due to the internal abundance of nutrients. The internal nutrient concept has been used in a number of models [Lehman et al., 1975; Jorgensen, 1981; Bierman, 1981]. The previous models have been either zero dimensional models treating the entire lake as a well mixed body or complex three-dimensional models. Use of the internal nutrient concept in a single dimension model in depth causes some problems in the mixing of phytoplankton of different internal nutrient content from different layers due to settling and diffusion. Ideally, the model would track individual cells in diffusion, setting, uptake, and growth (cell division). This of course would be impossible. A second alternative would be to model phytoplankton in a layer by ranges of internal nutrient corttent thus preserving conservation of the mass of the nutrient and approximately preserving the potential for nutrient uptake and growth. This approach requires multiple phytoplankton nutrient subclasses within a layer and an apportioning of the diffusion and settling from each class. The method would require a substantial increase in computations and storage and was deemed to be excessive. 28 i c I A third alternative used in the present model treats all the phy- toplankton of one algal group within a layer as having the same internal nutrient content. Consequently, the solution of the equations conserves mass, but moderates the growth rate and uptake rate by essentially averaging the internal nutrient content of algae transferred into a layer with the algae already in a layer. The effect on growth rate is opposite to the effect on uptake rate as indicated in Figure V-5. The non-linear rela- tionship between internal nutrients and growth is marked by two approxima- tely linear ranges and a non-linear range. The effect of averaging the internal nutrient content will only introduce an error in the non-linear portion of the curve. The magnitude of the error will depend on the dif- ference between the internal nutrient content in adjacent layers. In the mixed layer where most of the growth will occur, the internal nutrient con- tent should not vary much between layers and the effect on growth should be minimal. The growth rate of all phytoplankton populations are also affected by the water temperature. Many functions developed to represent the temperature effect have, been used in other models [Lehman et al., 1975, Scavia, 1980; Walters, 1980; Bierman, 1981]. The temperature function should yield a maximum at the optimal temperature for nutrient saturated growth and decrease both above and below this level. The only function of acceptable form without an excessive number of parameters is given by Lehman [1975] as where x T opt T min T max T - T 2 ( Opt x ~ exp - 2.3 (T _ T opt min ) ) T < T opt x T - T 2 exp (- 2.3 ( --:-T--_..::..o¥,::-,t=-- ) max opt ) T < T opt fraction of reduced growth due to temperature (0 ( x ( 1) optimum temperature for nutrient saturated growth low temperature at which growth is reduced 90 percent .. O°C· high temperature at which growth is reduced 90 percent A typical temperature curve is shown in Fig. V-6. The function requires three reference temperatures: optimum, minimum, and maximum. The minimum was chosen at O°C which results in a 90 percent reduction in growth from the optimum temperature to O°C. The optimum temperature is well docu- mented for a number of algal species. The maximum temperature used is cho- sen such that the growth rate is reduced by 90 percent from the optimum temperature to the maximum. Therefore, the maximum temperature used will be slightly lower than a maximum temperature determined eXperimentally. 29 ,...... ~ cO S ;:L ....... ;:L '-' 0 -.-I +J ],0 0.5 0 0 Chla Internal Nutrient Concentration (N)(mg ~-I) 1.0 0.5 o ~ ______ L-______________ ~ ____________ ~ ______ _ o Chla (Nmin) Chla (Nmax) -I Internal Nutrient Concentration (N)(mg ~ ) Fig. V-5 .• Graphical representation of non-linear growth and linear uptake as a function of the internal nutrient concentration. 30 1.0 -- --- - - - -- --- -- -- -- - -- - - -- -- -- I- Z W -S:2 IJ.. IJ.. W 0 U w w 0.5 a:: ::::> ~ 0:: W 0- ~ W t- O 0 5 10 15 20 25 Topt 30 T max TEMPERATURE (OC) Fig. V-6. Graphical representation of the temperature effect on growth for a species . with an optimum growth rate at 2SoC(T ) and a 90 percent inhibition at 330 C(T ). • • O]?t.. • 0 ( ) max The low temperature n~nety percent ~nhlb~t1on ~s taken to be at 0 CT .• m1.n Phytoplankton populations decrease due to three processes: respira- tion, mortality, and settling. Settling is handled in the solution scheme for state variables with a fall velocity as discussed in Section III and Appendix B. The phytoplankton population is assigned a fixed settling velocity or in the case of the level three model, each algal group is assigned an independent value. The use of a single value is an approxima- tion since the settling rate for an algal species may vary considerably with the physiological condition of the algae [Thomann et al., 1975]. The stressed physiological condition is, however, most severe in the hypolim- nion where little or no growth occurs. Consequently, a single settling rate should not cause a significant error in the epilimnion. Respiration is the reverse process of photosynthesis in which organic matter is transformed into inorganic matter with a utilization of oxygen and a release of carbon dioxide. Respiration is simulated in the model by where k = r e T == Chla == R k 0 (T-20) r (Chla) respiration rate constant -1 (day ) temperature adjustment factor temperature in the layer (OC) concentration of chlorophyll-a -1 (mg R. ) The assumption is made that respiration is related only to the con- centration of chlorophyll=a and temperature. The effect is to treat respiration as indirectly related to productivity. The mortality rate in a phytoplankton population is modeled in the same manner as respiration. M = k e T-20 (Chla) m In model three, the mortality is computed separately from respiration as a non-grazing mortality rate. A second source of mortality is modelled by zooplankton grazing on each of the three algal populations. A schematic of the processes represented in each of the model levels is presented in Figure V-7. Differential equations for each level of modeling are presented in Figure V-So In summary, the biological component can be represented by three models of varying complexity. The first and most simple model simulates growth of one general phytoplankton population in a well mixed euphotic zone. The second simulates the change in phytoplankton population for each layer. The third and most complex model treats three populations of phy- toplankton in each layer in a two stage growth process of nutrient uptake 32 Level 1: Analysis for the entire mixed layer LIGHT AVERAGE TOTAL PHOSPHORUS Level 2: Analysis for individual layers LIGHT AVAILABLE PHOSPHORUS Level 3: Analysis for individual layers LIGHT ~HOSPHORUS AMMONIUM NITRATE PHOSPHORUS UPTAKE . NITROGEN 'UPTAKE PREDICTED CHLA GROWTH t-------... PREDICTED CHLA RESPIRATION MORTALITY SETTLING GRAZING' PREDICTED CHLA Fig. V-~ Schematic illustration of the processes affecting the predicted chlorophyll-a concentration in each level of the biological submodel. 33 w .j::o. Fig. V-So OIFFERENTIAL EQUATIONS FOR CHLOROPHYLL-A CONCENTRATION . IN THE DIFFERENT LEVELS OF CHLOROPHYLL-A SIMULATION Levell: Direct prediction of Chlorophyll concentration (not a differential equation). Level 2: Michaelis-Mention growth function. a Chla _ l L (AK aChla) + ~ aA Chla + K e T-20 Chla + K e T-20 Chla at A az z az A az m m r r Gmax 1(1 + 2I~7K2) P f(T)( r?-)( if+P) I + ~ + /K2 p Chla = 0 Level 3: Internal nutrient growth function aChla _ l L (AK aChla) + ~ aA Chla + K eT- 20 Ch1a + K eT- 20 Ch1a at A az z az A 3z m m r r - Gmax feT) [ 1(1 + 2I~7K2 • NC - Ch1a(Nmin) • NC I + ~+r?-/K2 min -1 Chla .. chlorophyll-a concentration (mg L ) A .. Area em2) 2 -1 Kz .. eddy diffusivity (m day ) v • settling velocity (m day-I) -1 Kmt ~ .. mortality and respiration rate (day ) e .. temperature adjustment coefficient T .. temperature (OC) -1 Gmax .. maximum growth. rate (day ) feT) • temperature function for growth I .. intensity of photosynthetically active radiation (pE.m-2 hr-1) PC - 9t1a(Pmin) ] Ch1a o kltk2 = light 1imiation and inhibition coefficients, respectively (pE.m-2's-1) K .. half saturdation constant for phosphorus p -1 limited growth (mg L ) -1 P .. available phosphorus concentration (mg L ) -1 NC .. internal nitrogen concentration (mg t ) Nmin .. minimum ratio of internal nitrogen to chlorophyll-a -1 PC .. internal phosphorus concentration (mg L ) .Pmin .. minimum internal ratio of internal phosphorus to chlorophy11-a , " and growth. The different levels of complexity can be used to meet d~f­ ferent budget constraints. In addition, the first level can be used as a highly simplified model for a lake where very little calibration data is available and where more detailed modelling is not warranted. 3. Nutrients Central to the modelling of phytoplankton growth is the modelling of the nutrient cycles which contribute to growth. The components actually incorporated into a model depend on the nutrient being modeled and the relative importance of the component in the overall nutrient budget. Typically a eutrophication model treats phosphorus as the most significant nutrient, but many models also include nitrate-nitrite and ammonium for nitrogen limitation. In the present model, phosphorus is the only limiting nutrient used in the level one and two models. Level three incorporates nitrogen limitation as a user requested option in the program. Ph08phorus - The first question to be answered in modeling phosphorus is which forms of phosphorus to model. Most phosphorus monitoring programs measure phosphorus as total phosphorus concentration. Only about 5-50 per- cent of the total phosphorus is in the form that is readily used by phy- toplankton. The readily accessible fraction is composed of orthophosphate and polyphosphate ions and is referred to as soluble reactive phosphorus (SRP) or dissolved reactive phosphorus (DRP). Another fraction of phosphorus exists in dissolved and particulate organic matter and is released as the organic matter decays. Phosphorus may also be adsorbed onto sediment particles or complexed with metal ions in solution and in the sediment. The present model treats only SRP as available phosphorus and, indirectly models organic phosphorus as detritus (BOD). The differen- tial equation representing the phosphorus cycle is present in Fig. V-9. In the level one and level two models, phosphorus consumption is directly related to growth. The ratio of chlorophyll-a to phospho'rus is treated as a constant yield ratio. Therefore, phosphorus consumption is given by where A P = 1l YpChla Chla -1 1l = growth rate (day ) Ypchla - mass yield ratio of phosphorus to chlorophyll-a -1 Chla = concentration of chlorophyll-a (mg t ) The yield ratio is developed from a ratio of carbon to phosphorus [Stumm and Morgan, 1981] and the fixed chlorophyll-a/carbon ratio discussed earlier. In a similar manner, respiration directly releases phosphorus to the water. Mortality does not directly release phosphorus, but contributes to the detrital mass and the phosphorus is released through biological decay. 35 I"J C'I Fig. V-9. DIFFERENTIAL EQUATION FOR PROCESSES AFFECTING THE CONCENTRATION OF BIOLOGICAL AVAILABLE PHOSPHORUS S 3 ap _ !. L (Al{ ap) _ ..R. ( aA) - I K a T-20 at A az z az A az n=1 rn rn 3 P - I K a T-20(pc - P Chia ) - k a T-20(BOD) YPBOD mn mn n min n 2 2 n=1 3 T-20 ChIanPmaxn - PCn P + I UPmax aup (Chia (Pmax - Pmin ) ) (HSCP + p) Chian = 0 n=l n n n n n n -1 P - available phosphorus concentration (mg t ) 2 A - area (m ) 2 -1 K2 = eddy diffusion coefficient (m day ) -2 -1 Sp = sediment release rate (gm m day ) n = algal group n K , K , k- = reaction rate coefficients for respiration, mortality, and organic decay, r m -~ -1 respectively (day ) a = temperature adjustment coefficient -1 PC = internal phosphorus concentration (mg t ) Pmax = maximum internal nutrient ratio of mg phosphorus to mg chlorophyll Pmi = minimum internal nutrient ratio of mg phosphorus to mg chlorophyll n -1 Chla = chlorophyll concentration (mg t ) -1 BOD = concentration of BOD (mg t ) YPBOD = ratio of mg phosphorus released to mg oxygen utilized in organic decay -1 UP = maximum uptake rate for phosphorus (day ) max . -1 HSCP = half saturation constant for uptake of phosphorus (mgt ) In the level three model, uptake of phosphorus by algae has a one-to- one correspondence with phosphorus depletion in the water. The internal phosphorus is independent of the carbon content and, therefore, no yield coefficient is required. Respiration directly releases excess phosphorus into the water. Mortality directly releases excess phosphorus, but the minimum chlorophyll-a/ phosphorus ratio is used to apportion a constant fraction of the internal phosphorus to the detrital mass. This assumes that all excess phosphorus is stored as polyphosphates which go into solu- tion as available phosphorus when the cells die [Jensen et a1., 1976]. Similarly, zooplankton grazing releases excess phosphorus to the water and retains the minimum level of phosphorus. Consequently, the phosphorus/mass ratio is assumed to be constant in both zooplankton and detritus. This assumption ignores any storage of phosphorus within zooplankton, but since the zooplankton biomass is small and the fraction of phosphorus in zooplankton is small, the total stored phosphorus within zooplankton should not introduce a serious error, but only result in slightly faster cycling of phosphorus. The sediment of a lake contains a large store of phosphorus in organic matter, absorbed to sediment, and comp1exed with metal ions. The release of phosphorus from the sediments is generally a rather slow, constant rate. However, under anoxic conditions, the comp1exed portion can be rapidly mobilized causing a significant increase in the release rate. Presently, the release rate is a constant, but there is a need to treat the higher anoxic release rate to accurately model the phosphorus cycle. Bitrogen - Nitrogen limitation is modelled only in the case of the first and second algal classes (diatoms and green algae) in model level three. Blue-green algae are assumed to be capable of fixing dissolved nitrogen (N2) under nitrogen limitation and therefore are considered to be phosphorus and light limited only. Nirogen, like phosphorus, exists in many forms. However, uptake occurs almost exclusively from the nitrate- nitrite fractions and the ammonium fraction. There is also a preference in the uptake of ammonium over the uptake of nitrate-nitrite. In the pre- sent model, both nitrate-nitrite and ammonium are modeled separately. Schematic diagrams of each cycle are presented in Fig. V-lO. The differen- tial equations for ammonium and nitrate-nitrite are present in Figures V-11 and V-12, respectively. Ammonium is modeled in a manner very similar to phosphorus. However, Figure V-lO indicates all respiration, mortality, and grazing contribute to ammonium and only sediment exchange and nitrification contribute to the concentration of nitrate-nitrite. The preference for ammonium over nitrate-nitrite in nitrogen uptake is modelled after Scavia [1980] as u' = ( Chla(N ) - N max c )( NH+NO ) ~+NH+NO 37 g .... u III ... .... ~ iii ., ... aortallt res iration grazing Fig. V-IO. Illustration of the processes (arrows) and components (circles) comprising the nitrogen cycle as treated in the MNLAKE water quality program. 38 W \0 ~ r ~ Fig. V-H. DIFFERENTIAL EQUATION FOR PROCESSES AFFECTING THE AMMONIUM CONCENTRATION aNH 1 a (Al{ aNH) I3 at - A az z az- - S NC K eT- 20 + K...B T- 20NH - .J! ( aA) - k eT- 20 (BOD)YNBOD n rn rn -Y N A a z 2 2 3 Chla Nmax -NC 3 + rUNe T-20 ( n n n .:I ( NH :\ NH+NO T 20 "=1 max UN Chla (Nmax -Nmin )1 HSCNH + NHJ (HSCN +NH+N03 Chlan - 1. Kmn8 - ~NC -ChI a Nmin) = 0 u n n n n n n n n=1 mn n n n -1 NH = ammonium concentration (mg 1. ) 2 A - Area (m ) 2 -1 Kz - eddy diffusion coefficient (m day ) n = algal group -1 NC = internal nitrogen concentration (mg 1. ) k ,~-, k_, K = rate coefficients for respiration, nitrification rn -~ -"2 mn. -1 organic decay and mortality, respectively (day ) e = temperature adjustment coefficient T = temperature (OC) S '" N -2 -1 sediment release coefficient (gm m day ) BOD -1 = concentration of BOD (mg 1. ) YNBOD = yield coefficient for mg N per mg BOD N = maximum internal nutrient ratio of mg nitrogen per mg chlorophyll max N i '" minimum internal nutrient ratio of mg nitrogen per mg chlorophyll mn HSCNH '" half saturation constant for preferential uptake of ammonium -1 over nitrate (mg 1. ) -1 HSCN = half saturation constant for total nitrogen uptake (mg 1. ) -1 UN '" maximum uptake rate for nitrogen (day ) max Fig. V-12. DIFFERENTIAL EQUATION FOR PROCESSES AFFECTING THE CONCENTRATION OF NITRATE-NITRITE S 3 aNO _ !. L (AI{ aNO) _ IL-A T-20NH _ NO ( aA) + I UN 9 T-20(1 _ NH ) at A az z az -~- N A a z n=1 max UNn HSCNHn ~ 1I.'IU Chla Nmin - NC ( Chla (Nmaxn_ N n ) )( NH+NO ) Ch1a = 0 ----- +NH+NO n n n min n -1 NO = nitrate-nitrite concentration (mg t ) 2 A = area (m ) 2 -1 Kz - eddy diffusion coefficient (m day ) z - depth (m) -1 ~ = nitrification rate coefficient (day ) n ~ 9N,SUN = temperature adjustment coeffcients for nitrification and uptake, respectively. -1 NH = ammonium concentration (mg t ) -2 -1 SNO = sediment exchange coefficient (gm m day ) n - algal group -1 UN = maximum uptake rate for nitrogen (day ) max HSCNH - half saturation coefficient for preferential uptake of -1 ammonium over nitrate (mg t ) -1 ChI a - chlorophyll concentration (mg t ) Nmax = maximum internal nutrient ratio of mg nitrogen per mg chlorophyll N i = minimum internal nutrient ratio of mg nitrogen per mg chlorophyll m n -1 HSCN = half saturation constant for nitrogen uptake (mg t ) .. .' and where ~O = u' (1 - l)qH w: NH) nitrate-nitrite ~H .. u' ( l)qH m:. NH) ammonium u' = the internal and external nitrogen limited uptake rate (day-I) Chla a concentration of chlorophyll-a (mg 1-1) N = maximum nitrogen/chlorophyll-a ratio max Nmin = minimum nitrogen/chlorophyll-a ratio N = internal nitrogen concentration (mg 1-1 ) c -1 NH - ammonium concentration (mg 1 ) -1 NO = nitrate-nitrite concentration (mg 1 ) ~ = half saturation constant for nitrogen . -1 ~ - ammonium preference constant (mg 1 ) -1 ~O = uptake rate of nitrate-nitrite (day ) -1 UNH = uptake rate of ammonium (day ) -1 uptake (mg 1 ) The ammonium preference ratio has the form of a Michaelis-Menton ratio and varies from 0 to 1. The preference constant (~~) is then defined as the ammonium concentration at which half of the nit~3gen uptake is from ammonium and half from nitrate-nitrite. The only sources of nitrate-nitrite are sediment exchange and nitrifi- cation. Nitrification is the microbial mediated conversion of ammonium to nitrite and nitrite to nitrate. It is treated here as a one-step COnver- sion given by (Scavia, 1980]. where A NO .. ~O e T-20 (NH) ANO - increase in nitrite-nitrate ~O .. conversion rate coefficient ... 1 -1 (mg 1 day) -1 (day ) e .. temperature adjustment coefficient T .. temperature in a layer (OC) NH - ammonium concentration (mg 1-1, 41 4. Dissolved Oxygen and Detritus (BOD) The concentration of dissolved oxygen is seldom a state variable in eutrophication models. Oxygen does not directly affect phytoplankton. However, dissolved oxygen directly affects biological decay processes and phosphorus release under anoxic conditions. In addition, the vertical migration of zooplankton is inhibited to layers of low oxygen con- centration. It is included as a state variable in the present model to aid in simulating the movement of large zooplankton, to account for the micro- bial decay of organic matter (treated as BOD) which is a major component in the nutrient cycles, and to simulate aerators used in some lake restoration projects. In the model, dissolved oxygen is increased from two sources. First, at the lake surface, diffusion of oxygen occurs to attain an equilibrium between the dissolved oxygen concentration in the air and water. Secondly, the photosynthetic process produces oxygen directly from the breakdown of carbon dioxide and water [Wetzel, 1983]. The surface diffusion is modelled using the saturation deficit by where by AO k -1 -1 = change in dissolved oxygen (mg 1 day) -1 = exchange coefficient (day ) -1 = saturated oxygen concentration (mg 1 ) o sat 01 = oxygen concentration in the first layer (mg 1-1) The saturated oxygen concentration is temperature dependent and given o = 14.652 - 0.41022T + 7.99(10-3)~ - 7.7774(10-5)T3 sat where T = the surface temperature (OC). The oxygen exchange coefficient is not a constant and varies widely [Mattingly, 1977]. Many relationships have been developed to estimate the oxygen exchange coefficient [Yu et al., 1977]. In the present model, the oxygen exchange coefficient is predicted from where k = 0.641 + 0.128(.447 w)2 dz f ~ exchange coefficient w - wind velocity (mph) dz = thickness of first layer (m) 42 .. Oxygen is depleted from a layer by the combined actions of respira- tion, sediment uptake, nitrification, and organic decay of detritus (BOD). Respiration is treated in the same manner as stated previously, but with a different yield coefficient. Likewise, sediment uptake is treated in the same manner as sediment exchange of phosphorus, but with different coef- ficients. Oxygen depletion from nitrification occurs from the utilization of three moles of oxygen for every mole of ammonium oxidized. Therefore, the yield coefficient is simply given by the ratio of molecular weight of oxygen and nitrogen multiplied by three. The microbial decay of organic matter follows the same pattern as nitrification, but is a function of the detrital mass expressed in oxygen equivalents. The mortality of cells and a fraction of the grazed phy- toplankton are converted from chlorophyll-a concentrations to oxygen equivalents from the constant carbon/chlorophyll-a ratio and stoichiometric relationships [Stumm and Morgan, 1981]. The result is a one-to-one correspondence between the decay of detritus and the utilization of oxygen. AO ... ABOD - ~OD e T-20 (BOD) where -1 -1 AO = decrease in oxygen (mg t day) -1 -1 ABOD = decrease in detritus as oxygen equivalents (mg t day) -1 kBOD - first order decay coefficient (day ) e = temperature adjustment coefficient T ... temperature (OC) -1 BOD ... detritus as oxygen equivalents (mg 1 ) The linked BOD and dissolved oxygen equations are presented in Figure V-13. 5. Zooplankton The treatment of zooplankton in the model differs from the treatment of all other state variables. Zooplankton movement is not simply a cOm- bination of settling and diffusion, but the result of the individual mobi- lity of the zooplankton. In addition, the zooplankton population, effect of grazing, and nutrient recycling are related to the number of zooplankton rather than to the mass concentration. Therefore a simple. advection-diffusion model does not apply. The movement of zooplankton (especially large zooplankton) follow a pattern of vertical migration from deep water to the surface on a daily time scale. The assumption is made here that the retreat to deeper waters during the day provides zooplankton with a refuge from visual predators. At dusk, the zooplankton begin a vertical rise to the surface while grazing and return to deeper layers by dawn. The day depth of the zooplankton is fixed in the model at the deepest !tyer in which the dissolved oxygen concentration is greater than 0.5 mg 1 • Predation on zooplankton is taken to be a function of light availability at the day depth [Wright et al., 1980). 43 .j:-o -'=" Fig. V-13. COUPLED DIFFERENTIAL EQUATIONS FOR BOD AND DISSOLVED OXYGEN BOD: 3 Chla aBOD _ !. L (AK BOD) + V aA BOD _ I K a T-20 __ n + k_e T-20(BOD) = 0 at A az z z A az n=l mn mn YCBOD --Z BOD DISSOLVED OXYGEN: aDO 1 a (AK DO) + 1 ~ at - A az z Z YCH02 L . n=l (k eT- 20_ G f(T) min[L,P,N])Chla + k_ eTBO-D20(BOD) rn rn max n --z k m S + ...!. aA + 1 ~_~ e T-20(NH) - k (DO - DO) = 0 A a z YNH02~"NH NH sat surface -1 BOD = concentration of BOD (mg! ) A = Area (m2 ) 2 -1 K = eddy diffusivity (m day ) z -1 V = settling velocity (m day ) n = algal group k = mortality and respiration rate coefficients, r -1 respectively (day ) f(T) = temperature function for growth L = light limited growth ratio P = phosphorus limited growth ratio N = nitrogen limited growth ratio Sb = rate of oxygen utilized at the water -2 -1 sediment interface (gm m day ) YNH02 = ratio of mg ammonium consumed per YCBOD - ratio of mg chlorophyll to mg oxygen mg oxygen utilized in nitrification utilized in organic decay e = temperature adjustment coefficients T = temperature .(OC) -1 k2 = organic decomposition rate (day ) -1 DO = dissolved oxygen concentration (mg! ) G = maximum growth (day -1) max ~ = nitrification rate coefficient DO t = saturated oxygen concentration sa -1 k - oxygen exchange rate (day ) YCH02 ratio of mg chlorophyll per mg -1 (day ) -1 (mg! ) oxygen utilized and released in photosynthesis and respiration . ' The two periods of zooplankton activity are treated separately in the model. First, the day depth, light level at the day depth, and predation on zooplankton are computed. The respiratory and mortality losses contri- bute to the oxygen and nutrient budgets of the day depth layer in propor- tion to the length of the daylight period. Grazing at the day depth is not treated in the model. A second sequence of computations begins with the vertical rise in the evening. Grazing is assumed to occur in propor- tion to the ch1orophy11-a concentration in a layer compared to the total ch1orophy11-a concentration from the day depth to the surface. The effect of the proportion of chlorophyll in each layer is to simulate the time spent in each layer where the time (At) is given by where 3 !: Ch1aik k-l ) "";;';"""::-d ---:3.----- !: !: Ch1aik i=1 k-l At .. time increment for layer i (day) TD .. photoperiod (hr) Ch1aik - ch1orophy11-a concentration of phytoplankton group k in layer i • d .. day depth layer Grazing within a layer is also limited by temperature and ch1orophy11-a concentration. The effect of ch1orophy11-a concentration on grazing is expressed by a Michaelis-Menton ratio with a refugium effect for no grazing below a threshold ch10rophy11-a concentration. The grazing term for each algae class is given by where Chla-Ch1a G(l) .. G min e T(I)-20 ZP At V(lZ) (CF) max Ch1a-Ch1amin + HSCG G Vel) G(l) G max Ch1a i mn -1 .. grazing in layer I (mg Ch1a 1 ) . 1 .. maximum grazing rate (mg Chla (ind day)- ) .. minimum ch1orophy11-a concentration for grazing to occur (mg R. -1) -1 .. half saturation constant for grazing (mg Chla 1 ) .. temperature adjustment coefficient .. temperature in layer I (OC) 45 ZP At V(IZ) V(I) CF -3 = zooplankton concentration at the day depth layer (11m ) ~ time spent in layer I (day) (see above) '" volume of day depth layer (m3) :II volume of layer I (m3 ) '" conversion factor for m3 to !I. The zooplankton population is simulated by a constant reproduction rate and a time varying mortality/predation rate. The constant reproduction rate assumes that grazing will occur on the algae classes as well as on other organic material and~ therefore~ the chlorophyll-a concentration alone will not significantly affect the reproduction rate. The zooplankton population is largely treated as responding to preda- tion pressure. Predation is treated as occurring during the nocturnal migration and at the day depth. The nocturnal predation rate is computed for each layer between the day depth and the surface and is given by P(I) '" P (ZP V(IZ) - ZP i )At n V(I) m n where P(I) - nocturnal predation rate in Layer I during migration (1m-3) P '" nocturnal predation rate (day-I) n ZP = zooplankton population in the day depth layer (#m-3) V(IZ) = volume of the day depth layer (m3 ) V(I) .. volume of layer I (m3) ZP i = minimum concentration of zooplankton for predation to occur mn (fhn-3) At = time spent in layer I (day) (see above) The daytime predation is more complicated than the nocturnal predation and is composed of two components: a time varying maximum predation rate and a light limitation on predation. The time varying predation rate is a linear function describing change in predation rate from a minimum to a maximum between two specified days. The form of the function is DY - D ( min ) Pd .. P i + (P - Pmi ) D - D m n max n max min and 46 DY - Dmin o < < 1 Dmax - Dmin where daytime -1 Pd = predation rate (day ) -1 P = minimum seasonal predation rate (day ) min -1 P = maximum seasonal predation rate (day ) max DY .. Julian day D min .. Julian day of last day of minimum predation rate. D .. Julian day of beginning of maximum predation rate. max The time varying predation rate allows the model to simulate a low predation period and a high predation period with a linear interpolation in between. The time vary rate can also represent a measured population without actually simulating the zooplankton population. The light limitation on predation is based on work by Wright et al. [1980] and assumes a linear variation of predation between two light levels and where XI - XI i m n Pd = Pd XI - XI i max m n XI - XI i m n O UIMAX ::: maximum phosphorus uptake rate (day-!) WCHANL = width of inlet channel . WS'lR ::: wind sheltering ooefficient XIMAX ::: light intensity for maxiDn.ln zooplankton predation XIMIN ::: minimlln light intensity for zooplankton predation XKl ::: light attenuation coefioient for water (m- 1 ) (uE (m' s) -1) (uE (ml a)-1) XK2 ::: light attenuation coeficient for chloropbyll-a (1 (mg m) -1) XKM ::: algal. mortality coeffioient (day-1 ) XKM ::: algal mortality coefficient (day-I) XKNNH ::: nitrification rate (day-1) XKRl ::: algal respiration rate in the euJDotic zone (day-1 ) XKR2 = algal respiration rate below the euphotic zone (day-1 ) XKRZP ::: zooplankton respiration rate ( individual day-1) XNMAX = maximun nitrogen to chlorophyll mass ratio XNMIN ::: miniDn.In nitrogen to chlorophyll mass ratio Y2CH02 = mass ratio of chloro}il,y1l to oxygen in algal respiration 59 YCA YCBOD YCH02 YNHBOD YNHZP YPBOD YPZP yzoo Y'lW Z ZP ZR1IN = carbon to chlorophyll mass ratio = mass ratio of chlorophyll to oxygen deIIand in detritus = mass ratio of chlorophyll to oxygen in photosynthesis = mass ratio of 8IIIIlOIlium to oxygen in detritus = mass of WIII10Ilium per individual zooplankter = mass ratio of phosphorus to oxygen in detritus = mass of phosphorus per individual zoop1ankter = oxygen demand per individual zooplankter in zooplankton decay = mass of individual zooplankter = depths in the initial condi tiona (m) = initial concentration of zooplankton (large Clad.ocera) (# m- a) = minimum zooplankton for predation (# m- 3 ) 60 .. Meteorological Data File ~ IIDAYS MYBAR TAIR (1 TO IIDAYS) TDEW (1 TO KDAYS) PRECIP (1 TO KDAYS) WIND (1 TO KDAYS) DRCT (1 TO IIDAYS) OR (1 TO IIDAYS) RAD (1 TO IIDAYS) [Month, Days in the month, year] [Average Air Temperature (degrees F») [Dew Point Temperature (degrees F)] [Precipi tation (inches)] [Wind Speed (mph)] [Wind direction] [Percent Sunshine] [Solar radiation (Langleys)] [Note: The Meteorological Data File does not have to begin or end on the actual dates simulated, but must cover the entire period of simulation] Example MeteoroloSical Data File 5 31 1 Month, Days in the month, Year 43 49 45 45 55 62 62 64 46 47 42 49 65 54 63 63 75 ----- Average Air Temperature 63 58 51 49 59 55 58 55 59 60 66 71 65 57 ------------ 34 39 24 26 37 42 48 43 36 40 32 34 41 38 28 40 54 ----- Dew Point Temperature 49 43 38 34 45 38 37 33 39 39 42 56 58 44 ------------ .4 .29 O. O. O. O. O •• 25 .52 .4 O. O •• 02 .06 O. O •• 93 --- Precipitation O •• 22 O. O •• 56 O. O. O. O. O. O •• 41 .5 O. -------------- 8.2 11.1 16.0 6.9 12.5 16.5 12.1 8.9 14.2 9.5 14.2 7.3 ---------- 11.2 6.2 6.6 15.8 17.6 8.3 8.9 12.8 4.9 14.4 14.1 6.6 5.5 9.1 Maximum Wind Speed 8.1 8.9 10.5 12.7 10.9 -------------------------------------- 13 36 32 32 36 18 18 4 4 27 27 22 27 27 32 18 22 32 32 36 ----- Wind Direction 18 18 36 4 13 4 32 13 32 32 32 ----------------------------- 21 14 97 61 49 58 60 10 6 6 14 95 84 88 100 89 60 47 80 47 ---- Percent Sunshine 91 15 65 99 86 50 89 96 42 13 30 ----------------------------- 226 151 629 618 512 499 363 183 71 132 196 535 530 642 697 629 --- Solar Radiation 521 436 539 280 650 357 546 713 709 502 695 725 528 318 183 I I (file continues for all DIOIlths in the siDlulation period.) 6J Inflo~ D~t~ File NFLOW KONTH DAY INFLOW (l) TINll) PAIN(1) BODINIl) DOIN(1) 1'88(1) TD8(1) 1Il3IN(1) N84IN(1) CULAIN(1,l ) CHLAIN(1,2) CHLAIN(1,3) I I KONTH DAY INFLOlI(N) TIN(N) PAIN(N) BODIN(N) DOIN(N) 1'88(N) TD8(N)1Il3IN(N) NH4IH(I) aJLAIH(N,l) CULAIH(H,Z) CHLAIHIH,3) [Note: up to 6 inflow sources allowed (N<6). 'Jhe Inflow Data file does not have to begin or end on the actual beginning or ending dates of the simulation, but must cover the period of simulation. Outflow data is input as a negative n\.Ullber and is not followed by a line of inflow concentrations.] Ex~ple Inflo~ D~t~ File 3 5 1 .2 o .006 11.1 0 0.00 0 .073 .054 0 0 0 5 1 -21.7 [Negative n\.Ullber indicates an outflow rate] 5 1 2.7 o .024 13.7 0 0.000 0 .29 .22 0 0 0 5 2 4.6 o .006 11.1 0 0.00 0 .073 .064 0 0 0 5 2 -22.8 5 2 2.7 o .024 13.7 0 0.000 0 .29 .22 0 0 0 I I (file continues to cover all days in the simulation period) BODIN CIU.AIN DAY DOIN INFWW K>NTH NFLOW NH4IN N03IN PAIN TAIN TDS TSS Lake Input D~ta V~ri~bles = inflow concentration of detritus as BOD (mg 1-1) = inflow concentration of chloroJilyll in algal classes 1 thru 3 (mg 1-1) = day on which inflow occurred = inflow concentration of dissolved. oxygen (mg 1-1) = average inflow/outflow rate (cfs) for the day = month in which inflow occurred = number of inflow/outflows per day in present data file = inflow concentration of 8IIIIIOI1il. (mg 1 -1 ) = inflow concentration of nitrate-nitrite (mg 1-1) = inflow concentration of available P1osP'wrus (mg 1-1) = tem.pei-ature of inflow (degrees C) = inflow concentration of total dissolved. solids (. 1-1) = inflow concentration of total suspended solids (Jug 1-1) 62 " APPENDIX B DISCRETIZATION AND SOLUTION PROCEDURES FOR DIFFUSION AND ADVECTION-DIFFUSION EQUATIONS 63 Diffusion Equation (Dissolved Substances) or Initial Form: a Cia (AK ac) at"-Xaz zaz ± Discritized form: S = 0 2 c - c - A K (i,j+1 i-1,j+1 )~ ± S (C ) ± S = 0 i-1/2 i-1/2 ltA + A ) ~ 1 i,j+1 2 2'LlZi Llz i _1 Solution form: A K At ( H1/2 i+1/2 ) C A A 1 i+1 "+1 = Ci , j 4= 82 i LlZi 'f Z(MBOT) USE c************ IF REQUESTED C DUMMY VALUE OF -1 ****************** WRITE(1,REC-IREC) PLT(KK) 405 406 IREC=IREC+1 DO 405 I-1,NSET WRITE(1,REC-IREC) XNUL IREC=IREC+1 IF(MODEL.EQ.3) THEN DO 406 K-1,NCLASS WRITE(1,REC-IREC) XNUL IREC=IREC+1 ENDIF 410 CONTINUE ENDIF WRITE(1,REC-IREC) CHLMEAN IREC=IREC+1 X=ZP/TV(MBOT)*ZMAX WRITE(1,REC=IREC) X IREC-IREC+1 WRITE(1,REC-IREC) DMIX IREC=IREC+1 WRITE(1,REC=IREC) SECCHI IREC"IREC+1 IF(MDAV+MONTH*100.NE.NPRNT(NDAVS» THEN NF'"O WRITE(1,REC=IREC) REAL(NF) IREC"IREC+1 ENDIF RETURN END NAME TYPE CHLMEA DMIX MY EAR NF NSET SECCHI X XNUL 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 REAL REAL INTEGER*2 INTEGER*2 INTEGER*2 REAL TWO METER MEAN CHLOROPHYLL-A FROM SIMULATION (MG/L) MIXED LAVER DEPTH (M) CURRENT VEAR OF SIMULATION NUMBER OF FIELD DATA POINTS ON A PROFILE NUMBER OF STATE VARIABLES SIMULATED SECCHI DEPTH COMPUTED IN MODEL (IN LAKE ROUTINE) REAL REAL DEPTH VARIABLE TO MATCH LAVER DEPTH WITH DEPTHS IN PLT DUMMV OUTPUT VALUE OF -1 FOR DATA THAT DOES NOT EXIST SUBROUTINE PTABLE(XK,GRTOT) C* •• •• C •• *** SEND SIMULATION RESULTS TO THE TABULAR C**.** OUTPUT FILE (TAPE8.DAT) C·**"'* REAL*B A,V,TV,ATOP . COMMON/ZOOPL/IZ,MINDAV,MAXDAV,ZP,ZPMIN,PRMIN,PRMAX,PREDMIN,XIMIN, +XIMAX,XKRZP,GRAZMAX(3),THGRAZ(3),ASM,THRZP,HSCGRAZ(3),CHLAMIN(3), +REPRO,XI,XKMZ,GRAZE(3,40) COMMON/RESULT/ T2(40),CHLATOT(40),PA2(40),PTSUM(40),BOD2(40), +DS02(40)~C2(40),CD2(40),XN02(40),XNH2(40),CHLA2(3,40), +PC2(3,40/,XNC2(3,40),T20(40),SI2(40) COMMON/VOL/ZMAX,DZ(40),Z(40),A(40),V(40),TV(40),ATOP(41),DBL COMMON/STEPS/DZLL,DZUL,MBOT,NM,NPRINT,MDAV,MONTH,ILAY,DV COMMON/CHDICE/MODEL,NITRO, IPRNT(6) ,NDAVS,NPRNT(30) ,NCL ASS,PLT(SO) IF(MDAV.NE.O) THEN 85 60S 609 610 611 612 613 614 615 616 617 61S 619 620 621 622 623 624 625 626 627 62S 629 630 631 632 633 634 635 636 637 63S 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660 661 662 663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 200 199 WRITE(8,3000) ELSE WRITE(8,2999) ENDIF IF(IPRNT(4).LT.1.) THEN WRITE(8, 3008) DO 200 I-1,MBOT WRITE(8,3009) Z(I),T2(I),A(I),V(I) RETURN ENDIF IF(NCLASS.EQ.1)THEN WRITE(S,3001) DO 199 I-1,MBOT WRITE(8,3002) Z(I),T2(I),C2(I),CD2(I),CHLA2(1,I),PC2(1,I), + PA2(I),PTSUM(I),BOD2(I),DS02(I),DZ(I),A(I),V(I) ELSE 101 + 102 + 103 + 104 + WRITE(8,3003) IF(NITRO.GT.O) THEN DO 101 I = 1 , MBOT WRITE(8,3004) Z(I),T2(I),C2(I)tCD2(I),PA2(I),PTSUM(I), BOD2(I),DS02(I),XN02(I),XNH2(IJ,DZ(I),A(I),V(I) ELSE DO 102 I=1,MBOT WRITE(S,3005) Z(I),T2(I),C2(I),CD2(I),PA2(I),PTSUM(I), BOD2(I),DS02(I),DZ(I),A(I),V(I) ENDIF WRITE(8, 3006) IF(NITRO.GT.O) THEN DO 103 I=1,MBOT WRITE(8,3007)Z(I),(CHLA2(K,I),K=1,3),CHLATOT(I), (PC2(K,I),K=1,3), (XNC2(K,I),K a 1,3) ELSE DO 104 I=1,MBOT WRITE(8,3007)Z(I),(CHLA2(K,I),K=1,3),CHLATOT(I), (PC2(K, I) ,K=1,3) ENDIF ENDIF IF(MDAY.GT.O .AND. MODEL.GT.1) THEN X"ZP/TV(MBOT) XK=XK/TV(MBOT) GRTOT=GRTOT/TV(IZ) WRITE(S,1200)X,XK,GRTOT,Z(IZ) ENDIF 1200 FORMAT(///5X,'ZOOPLANKTON PARAMETERS',/,5X,22('-'),/,6X, + 'CONC. (H/M3)',2X,'PREDATION(H/M3)',3X,'GRAZING(MG/M3)',3X, + ' DAYDEPTH(M)',/,9X,F7.0,9X,F7.1,9X,F6.4,11X,F4.1,//) 1201 FORMAT(F11.4) 2999 FORMAT(//5X,'INITIAL CONDITIONS',/5X,1S('-'),/) 3000 FORMAT(//5X,'INFORMATION ON LAKE QUALITY',/5X,27('-'),/) 3001 FORMAT(5X,'Z(M) T(C) SS(PPM) TDS(PPM) CHLA(PPM) PC(PPM)', +' P(PPM) PT(PPM) BOD(MG/L) DO(MG/L) OZ(M)', +' AREA(M2) VOL(M3)') 3002 FORMAT(5X,3(F5.2,2X),1X,F5.1,4X,4(2X,F6.4),6X,F6.2,6X +F5.2,4X,F4.2,2(2X,E10.4» 3003 FORMAT(5X,'Z(M) T(C) SS(PPM) TDS(PPM) PA(PPM) PT(PPM)', +' BOD(MG/L) DO(MG/L) N03(MG/L) NH4(MG/L) DZ(M)', +' AREA(M2) VOL(M3)') 3004 FORMAT(5X,3(F5.2,2X),1X,F5.1,2(4X,F6.3),4X,F6.2,5X,F5.2,4X, +2(F6.3,6X),F4.2,2(2X,E10.4» 3005 FORMAT(5X,3(F5.2,2X),1X,F5.1,2(4X,F6.3),4X,F6.2,SX,FS.2,28X,F4.2, +2(2X,E10.4» 3006 FORMAT(//5X,'BIOLOGICAL PARAMETERS'/5X,21('-'),/,SX,'Z(M)',9X, +'CHLOROPHYLL',6X,'TOTAL',6X,'INERNAL P',11X,'INTERNAL N',/16X, +'1',6X,'2',6X,'3',5X,'CHLA',4X,'1',6X,'2',6X,'3 f ,6X,'1',6X,'2', +6X,'3',/,5X,10(' ',2X» 3007 FORMAT(5X,F5.2,2X,10(F5.3,2X» 3008 FORMAT(//,5X,' TEMPERATURE PROFILES''/,5X,'Z(M)',4X,'T(C)', +6X,'AREA (M2)',5X,'VOL (M3)') 3009 FORMAT(5X,2(F5.2,2X),2(2X,E10.4» RETURN END NAME GRTOT X TYPE REAL REAL REAL TOTAL ZOOPLANKTON GRAZING ON ALL ALGAL CLASSES (MG/M3) ZOOPLANKTON POPULATION (H/M3) XK NUMBER OF ZOOPLANKTON LOST TO PREDATION (H/M3) 679 SUBROUTINE SUBPLOT(NITRO,NF,MYEAR) 680 C***** 681 C***** PRODUCE PROFILE OF STATE VARIABLES AND 682 c***** FIELD DATA (IF AVAILABLE) 683 C***** "684 REAL*8 A,V,TV,ATOP 86 1 1 1 1 1 1 1 1. 1 1 685 CHARACTER*1 XS,CH1,CH2,CH3,CH4 686 CHARACTER*32 TITLE(10) 687 CHARACTER*4 MNTH(12) 688 COMMON/SOURCE/RM(3,40) ,PROD(40) ,XMR(3,40),PRODSUM(40) 689 COMMON/STEPS/DZLL,DZUL,MBOT,NM,NPRINT,MDAY,MONTH,ILAY,DY 690 CDMMON/RESULT/ T2(40),CHLATOT(40) ,PA2(40),PTSUM(40) ,BOD2(40), 691 +DS02(40),C2(40),CD2(40),XN02(40),XNH2(40),CHLA2(3,40), 692 +PC2(3,40),XNC2(3,40),T20(40),SI2(40) 693 COMMON/FIELD/IFLAG(10),FLDATA(10,50),DEPTH(50),NFLD(10) SD 694 COMMON/VOL/ZMAX,DZ(40),Z(40),A(40),V(40),TV(40),ATOP(41),DBL 695 DIMENSION ZD(23),FD(23),VTM(43)IZV(43),VAR(40,10) 696 DIMENSION FCT(10) LEN(10),ISCAL,10),ICX(4) 697 EQUIVALENCE (T2(1),VAR(1,1»,(ZD(1),PROD(1»,(FD(1),PRODSUM(1», 698 +(VTM(1),RM(1,1»,(ZV(1) XMR(1,1» 699 EQUIVALENCE (CH1,ICX(1»,(CH2,ICX(2»,(CH3,ICX(3»,(CH4,ICX(4» 700 DATA ICX/16H1B,16H5B,16H32,16H4A/ 701 DATA ISCAL/1,3*-1,3*1,3*-1/ 702 DATA FCT/1,3*1000,4*1,2*1000/ 703 DATA MNTH/'JAN ','FEB ','MAR ','APR ','MAY ','JUN ','JUL ','AUG' 704 +'SEP ','OCT ','NOV ','DEC ,/ 705 DATA TITLE /'TEMPERATURE (C)' , 'TOTAL CHLOROPHYLL-A (UG/L)', 706 +'AVAILABLE PHOSPHORUS (UG/L)','SUM OF PHOSPHORUS (UG/L)', 707 +'DETRITUS AS BOD (MG/L)','DISSOLVED OXYGEN (MG/L)', 708 +'SUSPENDED SOLIDS (MG/L)','DISSOLVED SDLIDS (MG/L)', 709 +'NITRATE-NITROGEN (UG/L)','AMMONIUM-NITROGEN (UG/L)'/ 710 DATA LEN /16,27,28,25,23,4*24,25/ 711 1000 WRITE(*, 109) CH1,CH2,CH3,CH4 712 109 FORMAT(1X,4A1,/) 713 C ... LIST AND SELECT DESIRED STATE VARiABLE FDR PLOTTING 714 DO 100 1=1,8 715 100 WRITE(*,101) I,TITLE(I) 716 101 FORMAT(1X,I2,' s ',A32) 717 IF(NITRO.GT.O) THEN 718 DO 102 1=9,10 719 102 WRITE(*,101) I,TITLE(I) 720 ENDIF 721 WRITE(*,99) 722 99 FORMAT(/1X, 'CHOOSE DESIRED PLOT (ENTER Q TO QUIT): ',\) 723 READ(*,*,ERR=1001) IC1 724 X-O. 725 C ... CHANGE DEPTH TO NEGATIVE VALUES FOR PLOTTING WITH DEPTH 726 DO 110 I=1 t MBOT 727 ZV(I)=-Z(I) 728 VTM(I)=VAR(I,IC1) 729 IF(X.LT.VAR(I t IC1» THEN 730 X=VAR(I,IC1) 731 IND=IC1 732 ENDIF 733 110 CONTINUE 734 12=0 735 C ... IF FIELD DATA IS AVAILABLE, LOCATE FIELD DATA CORRESPONDING 736 C ... TO SELECTED STATE VARIABLES 737 IF(NF.GT.O) THEN 738 DO 111 I=1,NF 739 IF(FLDATA(IC1,I).GT.O) THEN 740 12=12+1 741 FD(I2)-FLDATA(IC1,I) 742 ZD(I2)=-DEPTH(I) 743 IF(X.LT.FD(I2» THEN 744 X-FD(I2) 745 IND=-1 746 ENDIF 747 ENDIF 748 111. CONTINUE 749 ENDIF 750 NPEN=1 751 NPEN1=1 752 IPORT=99 753 MODL=99 754 ZV(MBOT+1)-O. 755 ZD(I2+1)=0. 756 VTM(MBOT+1)-O. 757 FD(I2+1)=0. 758 FCTOR=0.72 759 C ... BEGIN PLOTTING SEQUENCE 760 990 CALL PLOTS(O,IPORT,MODL) 761 CALL SIMPLX 762 CALL FACTOR(FCTOR) 763 CALL NEWPEN(NPEN) 764 C ... DETERMINE MAXIMUM X & Y VALUES EITHER IN SIMULATED 765 C ... STATE VARIABLES OR IN FIELD DATA 766 IF(IND.GT.O) THEN 767 I3=MBOT+1 768 CALL SCALE(VTM,10.,I3,1) 769 DX-VTM(MBOT+3) 770 ELSE 771 13=12+1 772 CALL SCALE(FD,10.,I3,1) 87 773 774 775 776 777 778 779 780 781 782 783 784 785 786 787 788 789 790 791 792 793 794 795 796 797 798 799 800 801 802 803 804 805 806 807 808 809 810 811 812 813 814 815 816 817 818 819 820 821 822 823 824 825 826 827 NAME CHi CH2 CH3 CH4 DAXIS OX FCT FCTOR 12 IC1 ICX IND IPORT ISCAL LEN MNTH MODL MYEAR NF NITRO NPEN NPEN1 TITLE X XMDAY XMYEAR XS YD YSC DX"FD(I2+3) ENDIF IF(ZV(MBOT).LT.ZD(I2» THEN I3-MBOT+1 CALL SCALE(ZV,6.,I3,-1) YSC=-ZV(MBOT+3) ELSE 13=12+1 CALL SCALE(ZD,6.,I3,-1) YSC=-ZD(I2+3) ENDIF ZV(MBOT+2)=YSC ZD(I2+2)=YSC VTM(MBOT+2)=DX FD(I2+2)=DX DAXIS=DX*FCT(IC1) YD=ANINT(6.*YSC) CALL STAXIS(.2,.25,.2,.15,ISCAL(IC1» C ... DRAW X-AXIS CALL AXIS(1.,7.,TITLE(IC1),LEN(IC1),-10.,O.,O.,DAXIS) CALL STAXIS(.2,.25,.2,.15,1) C ... DRAW Y-AXIS CALL AXIS(1.,1.,'DEPTH (M)',10,-6.,90. ,YD,-YSC) XMDAY=MDAY XMYEAR=MYEAR C ... PRINT TITLE TO DIAGRAM CALL SYMBOL(4.5,.5,.25,MNTH(MONTH),O.t4) CALL NUMBER(999.,999. ,.25,XMDAY,.O,-1} CALL SYMBOL(999.,999.,.25,', ',0.,2) CALL NUMBER(999.,999.,.25,XMYEAR,.O,-1) CALL NEWPEN(NPEN1) CALL PLOT (1.,7.,-3) C ... PLOT SIMULATED PROFILES AS A LINE CALL LINE(VTM,ZV,MBOT,1,O,1) C ... PLOT FIELD DATA WITH A SYMBOL IF(I2.GT.0) CAlL LINE(FD,ZD,I2,1,-1,4) CALL PLOT(0.,O.,999) WRITE(*,130) 130 FORMAT(1X,'SEND TO HARDCOPY DEVICE? (Y/N) ',\) READ(*,'(A),///') XS IF(XS.EQ.'Y' .OR. XS.EQ.'Y') THEN WRITE(*,140) READ(*,*) IPORT,MODL 140 FORMAT(/1X,'ENTER PLOT88 IOPORT AND MODEL ',\) WRITE(*, 143) READ(*,*) NPEN,NPEN1 143 FORMAT(/1X,'ENTER LINE WEIGHT (AXIS,DATA) ',\) WRITE(*,141 ) READ(*,*) FCTOR 141 FORMAT(/1X,'ENTER REDUCTION FACTOR (>1.0 : ',\) GOTO 990 ENDIF GOTO 1000 1001 RETURN END TYPE CHAR*1 SCREEN CONTROL CHARACTER CHAR*1 SCREEN CONTROL CHARACTER CHAR*1 SCREEN CONTROL CHARACTER CHAR*1 SCREEN CONTROL CHARACTER REAL AXIS INCREMENTS AD~USTED FOR SAME UNITS AS IN AXIS LABEL REAL AXIS INCREMENT FROM SCALE REAL FACTORS TO CHANGE DATA TO SAME UNITS AS IN AXIS LABELS REAL REDUCTION FACTOR USED IN HARD COPY OUTPUT INTEGER*2 TOTAL NUMBER OF FIELD DATA PROFILE POINTS INTEGER*2 INDEX OF STATE VARIABLE TO BE PLOTTED INTEGER*2 SCREEN CONTROL VARIABLE INTEGER*2 FLAG TO DETERMINE IF MAXIMUM VALUE IS IN MODEL OF FIELD RESULTS INTEGER*2 OUTPUT CODE FOR PLOT88 INTEGER*2 FLAG FOR NUMBER OF DECIMAL PLACES ON AXIS NUMBERS INTEGER*2 LENGTH OF AXIS LABELS CHAR*4 MONTH INTEGER*2 OUTPUT DEVICE CODE FOR PLOT8S INTEGER*2 YEAR INTEGER*2 NUMBER OF POINTS IN FIELD DATA PROFILE INTEGER*2 FLAG FOR NITROGEN SIMULATION INTEGER*2 PEN NUMBER FOR LINE THICKNESS IN PLOT8S INTEGER*2 PEN NUMBER FOR LINE THICKNESS IN PLOT8S CHAR*32 AXIS LABEL TITLES REAL VARIABLE FOR MAXIMUM VALUE OF FIELD DATA OR SIMULATED RESULTS REAL DAY REAL YEAR CHAR*1 VARIABLE FOR RESPONSES TO INTERACTIVE QUESTIONS REAL MAXIMUM DEPTH ON Y-AXIS REAL INCREMENTS OF DEPTH ON Y-AXIS 88 828 SUBROUTINE TIMEPLT(~REC) 829 C***** 830 C***** PROOUCES TIME SERIES PLOTS ON DAVS SPECIFIED IN ARRAV NPRNT. 831 C***** DATA IS TAKEN FROM THE PLOT DATA FILE (TAPE8.PLT) 832 c***** 833 CHARACTER*3 MNTH(12) 834 CHARACTER*32 TITLE(17) 835 CHARACTER*64 LAKE 836 CHARACTER*1 LABEL(64),C1,C2,C3,C4 837 COMMON/TITL/ LAKE 838 COMMON/SUB/SDZ(60),SZ(60),LAV(40),AVG1(4,60),SVOL(60) 839 COMMON/SOURCE/RM(3,40),PROD(40) tXMR(3,40) ,PRODSUM(40) 840 COMMON/TEMP/PARIO(4) ,PCDUM(3,40},XNHP(40) ,XNOD(40), 841 + CHLADUM(3,40),XNCD(3,40) PADUM(40),SID(40) 842 COMMON/FLOW/HMK(41),QE(40),FVCHLA(5),PE(5,41) 843 COMMON/FIELD/ IFLAG(10),FLDATA(10,50),DEPTH(50),NFLD(10),SD 844 DIMENSION FCT(17),LEN(17),ISCAL(17),KDAVS(12) PLT(10) 845 DIMENSION SVAR(420)'DAV(368)'FVAR(52)'FDAV~52~'ICX(4) 846 EQUIVALENCE(C1,ICX(1»,(C2,ICX(2»,(C3,ICX 3) ,(C4,ICX(4» 847 EQUIVALENCE (DAV(1),PARIO(1»,(FVAR(1),HMK 1) , 848 +(FDAV(1),PE(1,1»,(SVAR(1)tSDZ(1» 849 EQUIVALENCE (LABEL(1),LAKE} 850 DATA ICX/16#1B,16#5B,16#32,16#4A/ 851 DATA KDAVS/31,28,31,30,31,30,31,31,30,31,30,31/ 852 DATA ISCAL/1,3*-1,3*1,7*-1,3*1/ 853 DATA FCT/1,3*1000,4*1,2*1000,4*1000,.001,-1,-1/ 854 DATA MNTH/'~AN','FEB','MAR','APR','MAV','~UN','~UL','AUG', 855 +'SEP','OCT','NOV','DEC'/ 856 DATA TITLE /'TEMPERATURE (C)','TOTAL CHLOROPHVLL-A (UG/L)', 857 +'AVAILABLE PHOSPHORUS (UG/L)','SUM OF PHOSPHORUS (UG/L)', 858 +'DETRITUS AS BOO (MG/L)','DISSOLVED OXYGEN (MG/L)', 859 +'SUSPENDED SOLIDS (MG/L)','DISSOLVED SOLIDS (MG/L)', 860 +'NITRATE-NITROGEN (UG/L)','AMMONIUM-NITROGEN (UG/L)', 861 +'CHLOROPHVLL (1) (UG/L)','CHLOROPHVLL (2) (UG/L)', 862 +'CHLOROPHYLL (3) (UG/L)','CHLOR-A 2M MEAN (UG/L)', 863 +'ZOOPLANKTON (1000/M**2)','MIXED LAYER DEPTH (M)', 864 +'SECCHI DEPTH (M)'/ 865 DATA LEN /16,27,28,25,23,4*24,25,3*23,22,23,22,16/ 866 DO 798 1=0,64 867 KLEN=64-I 868 IF(LABEL(KLEN).NE.' ,) GOTO 800 869 798 CONTINUE 870 800 WRITE(*,109) C1,C2,C3,C4 871 109 FORMAT(1X,4A1,//) 872 IPRNT=O 873 IREC=1 874 READ(1,REC=IREC) X 875 IPRNT6-INT(X) 876 DO 95 I=1,IPRNT6 877 95 READ(1) PLT(I) 878 READ(1) X 879 MONTH-INT(X) 880 READ(1) X 881 NM=INT(X) 882 READ(1) YEAR 883 MVEAR=INT(YEAR) 884 IF(AMOO(YEAR,4.).EQ.0.) KDAYS(2)=29 885 READ(1) X 886 MODEL-INT(X) 887 READ(1) X 888 NCLASS=INT(X) 889 READ(1) X 890 NITRO=INT(X) 891 READ(1) X 892 IPRNT4=INT(X) 893 DO 100 1=1,8 894 100 WRITE(*,101) I,TITLE(I) 895 101 FORMAT(1X,I2,' = ',A32) 896 IF(NITRO.GT.O) THEN B97 DO 102 1=9,10 898 102 WRITE(*,101) I,TITLE(I) 899 ENDIF 900 IF(MODEL.EQ.3) THEN 901 DO 103 1=11 10+NCLASS 902 103 WRITE(*,101) I,TITLE(I) 903 ENDIF 904 C ... LIST AND SELECT DESIRED DEPTHS AND STATE VARIABLES FOR PLOTTING 905 00 104 .1=14,17 906 104 WRITE(*,101) I,TITLE(I) 907 WRITE(*,99) 90B 99 FORMAT(//5X,'CHOOSE DESIRED PLOT . ',\) 909 READ(*,*) lC1 910 IF(IC1.LT.14) THEN 911 WRITE(*,801) 89 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 i 1 1 1 1 1 1 1 1 1 1 2 1 2 2 1 2 3 3 3 3 3 3 1 1 1 1 1 1 1 1 1 1 912 913 914 915 916 917 918 919 920 921 922 923 924 925 926 927 928 929 930 931 932 933 934 935 936 937 938 939 940 941 942 943 944 945 946 947 948 949 950 951 952 953 954 955 956 957 958 959 960 961 962 963 964 965 966 967 968 969 970 971 972 973 974 975 976 977 978 979 980 981 982 983 984 985 986 987 988 989 990 991 992 993 994 995 996 997 998 999 801 FORMAT(III.1X. / INTEGER',3X.'DEPTH'.I.1X,7(/-').3X.5(/- '» DO 803 1=1.IPRNT6 803 WRITE(*,804) I.PLT(I) 804 FORMAT(4X.I2.5X,F4.1) WRITE(*.802) READ(*.*) K2 802 FORMAT(II.1X./ENTER INTEGER FOR DESIRED DEPTH '.\) ENDIF IP=1 IC=IC1 C ... SET INCREMENTS FOR READING THE DIRECT ACCESS PLDT DATA FILE IF(IPRNT4.GT.0) THEN IP=8 IF(NITRO.GT.O) THEN IP=10 ELSE IF(IC1.GT.10 .AND. IC1.LT.14) IC=IC1-2 ENDIF IF(MODEL.EQ.3) IP=IP+NCLASS ENDIF IF(IC1.LT.14) THEN IR1=1+(IP+1)*(K2-1)+IC IR2=IP-IC+(IP+1)*(IPRNT6-K2)+5 ELSE IR1=IPRNT6*(IP+1)+IC1-13 IR2"18-IC1 ENDIF 98 IREC"9+IPRNT6 IC2=0 IC3=0 ISTEP=1 C ... READ SPECIFIED DATA FROM THE PLOT DATA FILE C ... SIMULATED RESULTS STORED IN SVAR C ... FIELD DATA RESULTS STORED IN FVAR DO 105 ..1=1.400 IF(IREC.GE.uREC) GOTO 120 READ(1.REC=IREC) DAVu IREC=IREC+IR1 READ(1.REC=IREC) VARu IF(VARJ.GE.O.O) THEN IC2=IC2+1 DAY(J)=DAYJ SVAR(J)=VARJ*FCT(IC1) ENDIF IREC=IREC+IR2 READ(1.REC=IREC) XNF IF(XNF.GT.O) THEN READ(1) XPRFLE IF(IC1.GT.10 .AND. IC1.LT.16) THEN IREC=IREC+XNF+XPRFLE+XPRFLE*XNF+4 GOTO 105 ENDIF FLD=O.O IF(IC1.GE.16) THEN IREC=IREC+XNF+XPRFLE+XPRFLE*XNF+IC1-14 READ(1.REC=IREC) FLD IREC=IREC+18-IC1 GOTO 106 ENDIF NF=INT(XNF) NPRFLE=INT(XPRFLE) DO 304 I=1.NF 304 READ(1) DEPTH(I) DD 305 I=1.NPRFLE READ(1) XNFLD 305 NFLD(I)=INT(XNFLD) DO 306 L=1.NPRFLE DO 306 K=1.NF IF(DEPTH(K).NE.PLT(K2) .OR. NFLD(L).NE.IC1) GOTO 306 IREC=IREC+NPRFLE+NF+NF*(L-1)+K+1 READ(1.REC=IREC) FLD IREC=IREC+NF*(XPRFLE-L+1)-K+3 GDTO 106 306 CONTINUE IREC=IREC+XNF+XPRFLE+XPRFLE*XNF+4 106 IF(FLD.GE.O) THEN IC3=IC3+1 FDAV(IC3)=DAY(J) FVAR(IC3)=FLD*FCT(IC1) ENDIF ELSE IREC=IREC+1 ENDIF 105 CONTINUE 120 IC2=J-1 NPEN=1 NPEN1=1 MODLE=99 90 " 1000 1001 1002 1003 1004 1005 1006 1007 1008 1009 1010 1011 1012 1013 1014 1015 1016 1017 1018 1019 1020 1021 1022 1023 1024 1025 1026 1027 1028 1029 1030 1031 1032 1033 1034 1035 1036 1037 1038 1039 1040 1041 1042 1043 1044 1045 1046 1047 1048 1049 1050 1051 1052 1053 1054 1055 1056 1057 1058 1059 1060 1061 1062 1063 1064 1065 1066 1067 1068 1069 1070 1071 1072 1073 1074 1075 1076 1077 1078 1079 1080 1081 1082 1083 1084 1085 1086 1087 IPORT=99 FCTOR=l. VO=l.25 IX"'IC2+IC3 C ... SELECT POSITIVE OR NEGATIVE V-AXIS DEPENDING ON C ... STATE VARIABLE SELECTED IF(IC3.GT.0) THEN DO 107 1=1 IC3 107 SVAR(I+IC2).FVAR(I) ENDIF IF(IC1.GT.15) THEN VO"5.25 IX"IX+l SVAR(IX)=O.O ISTEP=-l ENDIF C ..• BEGIN PLOT SEQUENCE 899 CALL PLDTS(O.IPORT,MODLE) CALL SIMPLX CALL FACTOR(FCTOR) CALL NEWPEN(NPEN) C ... SCALE V-VARIABLE TO DETERMINE MAXIMUM AND AXIS INCREMENT IF(IPRNT.LT.l) THEN CALL SCALE(SVAR,4.0,IX,ISTEP) FIRST=SVAR(IX+1) DV=SVAR(IX+2 ) IF(IC1.GT.15) THEN FIRST=ANINT(-4*SVAR(IX+2» DV=-FIRST*0.25 SVAR(IX+1)-O.O SVAR(IX+2)=DV ENDIF SVAR(IC2+1)=SVAR(IX+l! FVAR(IC3+1)·SVAR(IX+l SVAR(IC2+2)=SVAR(IX+2 *ISTEP FVAR(IC3+2)=SVAR(IX+2 *ISTEP ENDIF CALL SVMBOL(1.0,.65,O.15,TITLE(IC1),O.O,LEN(IC1» IF (IC1.LT.14) THEN CALL SVMBDL(999.,999.,O.15,' AT ',0.0,4) CALL NUMBER(999.,999.,O.15,PLT(K2),O.O,l) CALL SVMBOL(999.,999.,O.15,' METERS',O.O,7) ENDIF CALL SVMBOL(1.0,O.40,O.16,LAKE,O.O,KLEN) CALL SVMBOL(999.0,999.0,O.15,' ',O.Ot2) CALL STAXIS(.12,.15,.12,.12,ISCAL(IC1}) C ... DRAW V-AXIS CALL AXIS(1.0,l.25,TITLE(IC1),LEN(IC1),-4.0,90.,FIRST,DV) CALL STAXIS(0.O,l.0,l.0,. 12,0.0) . CALL AXIS(7.5,5.25,' ',1,-4. ,270.,FIRST,DV) MX·NM+MONTH-l XSCALE=O.O C ... DETERMINE LENGTH OF THE X-AXIS I2-MONTH DO 200 I=MONTH,MX IF(I2.EQ.13) THEN 12=1 VEAR=VEAR+l KDAVS(2)=28 IF(AMOD(VEAR,4.).EQ.0.) KDAVS(2)=29 ENDIF XSCALE=XSCALE+KDAVS(I2) 200 12=12+1 XSCALE=XSCALE/6.5 BEGIN-O.O DO 205 I=l,MONTH-l 205 BEGIN=BEGIN+KDAYS(I) DAV(IC2+1)=BEGIN OAV(IC2+2)"XSCALE FDAV(IC3+1)=BEGIN FOAV(IC3+2)-XSCALE CALL PLOT(1.0,l.25,S) CALL PLOT(1.0,l.10,2) KD"'O X=1.0 C ••• DRAW X-AXIS, RIGHT SIDE OF V-AXIS ANO TOP I2-MONTH DO 201 I-MONTH,MX IF(I2.EQ.13) 12=1 XX=X+KDAVS(I2)/XSCALE*.5-.15 CALL SVMBDL(XX,O.90,.12,MNTH(12),O.O,3) CALL PLOT(X,l.10,3) KD"KD+KDA VS (I 2) X=REAL(KD)/XSCALE+l.0 CALL PLOT(X,l.10,2) CALL PLOT(X,l.25,2) 12=12+1 201 CONTINUE 91 1088 1089 1090 1091 1092 1093 1094 1095 1096 1097 1098 1099 1100 1101 1102 1103 1104 1105 1106 1107 1108 1109 1110 1 111 1112 11 13 1114 1115 1116 1117 1118 1119 1120 NAME BEGIN C1 C2 C3 C4 DAYJ DY FCT FCTDR FIRST FLO IC IC1 IC2 IC3 ICX IP IPORT IPRNT IPRNT4 IPRNT6 IR1 IR2 IREC ISCAL ISTEP IX JREC K K2 KD KDAYS KLEN LEN MNTH MODEL MODLE MONTH MX MY EAR NCLASS NF NITRO NM NPEN NPEN1 NPRFLE PLT TITLE VARJ XS CALL PLOT(X,5.4,2) CALL PLOT(1.0,5.4,2) CALL PLOT(1.0,5.23,2) CALL PLOT(1.0,YO,-3) CALL NEWPEN(NPEN1) C ... PLOT LINE FOR SIMULATED VALUES CALL LINE(DAY,SVAR,IC2,1,O,1) CALL STLINE(1,.1,O.) C ... PLOT SYMBOL FOR FIELD DATA VALUES IF(IC3.GT.O) CALL LINE(FDAY,FVAR,IC3,1,-1,O) CALL PLOT(O.,O. ,999) WRITE(*,130) 130 FORMAT(1X,'SEND TO HARDCOPY DEVICE? (Y/N) ',\) READ(*,'(A),///') XS IF (XS. EQ. 'Y , . OR. XS. EQ. ' V') THEN WRITE(*, 140) READ(*,*) IPORT,MODLE 140 FORMAT(/1X,'ENTER PLOT88 10PORT AND MODEL ',\) WRITE(*,143) READ(*,*) NPEN,NPEN1 143 FORMAT(/1X, 'ENTER LINE WEIGHT (AXIS,DATA) ',\) WR I TE ( * , 141 ) REAO(*,*) FCTOR 141 FORMAT(/1X,'ENTER REDUCTION FACTOR ',\) IPRNT=1 GOTO 899 ENDIF WRITE(*,142) READ(*,'(A)') XS IF(XS.EQ.'Y' .OR. XS.EQ.'Y') GOTO 800 142 FORMAT(///,' REQUEST MORE GRAPHS (Y/N) ? ',\) RETURN END TYPE REAL FIRST DAY OF THE TIME SERIES PLOT CHAR*1 SCREEN CONTROL CHARACTER CHAR*1 SCREEN CONTROL CHARACTER CHAR*1 SCREEN CONTROL CHARACTER CHAR*1 SCREEN CONTROL CHARACTER REAL INPUT VARIABLE FOR JULIAN DAY REAL INCREMENT ON Y-AXIS REAL FACTOR TO CHANGE VARIABLE VALUES TO SAME UNITS AS AXIS TITLE REAL REDUCTION FACTOR TO BE USED IN HARD COPY OUTPUT REAL LOWEST VALUE ON Y-AXIS REAL FIELD DATA INPUT VARIABLE INTEGER*2 INDEX USED IN SET INCREMENTS FOR DATA IN DIRECT ACCESS FILE INTEGER*2 INDEX FOR STATE VARIABLE TO BE PLOTTED INTEGER*2 NUMBER OF SIMULATION POINTS ON THE PLOT INTEGER*2 NUMBER OF FIELD DATA POINTS ON THE PLOT INTEGER*2 SCREEN CONTROL VARIABLE INTEGER*2 VARIABLE USEO IN SETTING DATA INCREMENT IN DIRECT ACCESS FILES INTEGER*2 OUTPUT CODE IN PLOT88 INTEGER*2 FLAG FOR FIRST OR SECOND PLOT OF THE SAME DATA INTEGER*2 FLAG FOR TEMPERATURE ONLY SIMULATION INTEGER*2 NUMBER OF DEPTHS STORED ON THE PLOT DATA FILE INTEGER*2 DATA INCREMENT FROM DAYJ TO DESIRED SVAR ON PLOT DATA FILE INTEGER*2 DATA INCREMENT SVAR TO NF ON PLOT DATA FILE INTEGER*2 RECORD NUMBER ON PLOT DATA FILE INTEGER*2 FLAG FOR NUMBER OF DECIMAL PLACES ON AXIS NUMBERS INTEGER*2 FLAG FOR CONCENTRATION PLOTS (1) OR DEPTH PLOTS (-1) INTEGER*2 TOTAL NUMBER OF SIMULATED AND FIELD DATA POINTS ON THE PLOT INTEGER*2 TOTAL NUMBER OF RECORDS ON THE DIRECT ACCESS FILE INTEGER*2 INDEX FOR FIELD DATA INTEGER*2 INDEX FOR SELECTED DEPTH TO BE PLOTTED INTEGER*2 LAST DAY ON THE TIME SERIES INTEGER*2 NUMBER OF DAYS IN EACH MONTH INTEGER*2 LENGTH OF LAKE TITLE INTEGER*2 LENGTH OF AXIS LABELS CHAR*3 NAME OF MONTH INTEGER*2 BIOLOGICAL SIMULATION OPTION IN USE INTEGER*2 OUTPUT DEVICE CODE FOR PLOT88 INTEGER*2 MONTH INTEGER*2 LAST MONTH IN TIME SERIES INTEGER*2 YEAR INTEGER*2 NUMBER OF ALGAL CLASSES SIMULATED INTEGER*2 NUMBER OF FIELD DATA POINTS ON A PROFILE INTEGER*2 FLAG FOR NITROGEN SIMULATION INTEGER*2 NUMBER OF MONTHS TO BE SIMULATED INTEGER*2 PEN NUMBER FOR LINE THICKNESS IN PLOT8S INTEGER*2 PEN NUMBER FOR LINE THICKNESS IN PLOT88 INTEGER*2 NUMBER OF PROFILES IN FIELD DATA REAL DEPTHS AVIALABLE ON THE PLOT DATA FILE CHAR*32 TITLE FROM THE INPUT DATA FILE REAL INPUT VARIABLE FOR SIMULATED DATA TO BE PLOTTED CHAR*1 VARIABLE FOR RESPONSES TO INTERACTIVE QUESTIONS 92 XSCALE REAL X-AXIS INCREMENT VO REAL V-DIRECTION LOCATION OF ORIGIN PASS ONE NO ERRORS DETECTED 1122 SOURCE LINES 93 1 $NOFLOATCALLS 2 $STORAGE:2 3 SUBROUTINE ABOUND 4 C****'" 5 C"''''**''' 6 C"'*"''''''' 7 C"'''''''*''' 8 C"'**"'''' 9 C***"'''' COMPUTES THE SURFACE AREA OF EACH LAYER (ATOP) USING THE DEPTH AREA RELATIONSHIP IN LAKE. ATOP(1) = SURFACE AREA OF THE LAKE ATOP(MBOT+1) = 0.0 10 REAL*8 A,V,TV,ATOP,ADUM 11 COMMON/STEPS/DZLL,DZUL,MBOT,NM,NPRINT,MDAY,MONTH,ILAY,DY 12 COMMON/VOL/ZMAX,DZ(40),Z(40),A(40),V(40),TV(40),ATOP(41),DBL 13 DUM=O. 14 DO 100 I=1,MBOT 15 ZDUM=ZMAX-DUM 16 DUM=DUM+DZ(I) 17 CALL LAKE(ZDUM,ADUM,O,1) 18 100 ATOP(I)=ADUM 19 ATOP(MBOT+1)=O. 20 RETURN 21 END NAME TYPE ADUM DUM ZDUM 1 2 1 1 1 1 1 1 1 1 2 1 1 1 1 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 REAL"'8 REAL REAL AREA FROM LAKE ROUTINE DUMMY VARIABLE FOR DEPTH FROM THE SURFACE DEPTH FROM THE BOTTOM USED IN LAKE SUBROUTINE ADVECT(P,HED,NFLOW,S,FT,WCHANL,ST,MYEAR) C"'''''''** C"'*"'** INFLOW/OUTFLOW ROUTINES USING PLUNGING DENSITY C"'*"'** CURRENT ROUTINES ON INFLOW AND SURFACE OUTFLOW. C"''''*** ADDITIONAL INFLOWS AND OUTFLOWS CAN BE INCLUDED C"'*"'** WITH CALLS TO LAKE (SEE BELOW). C"'''''''** REAL*8 A,V,TV,ATOP,ADUM COMMON/FLOW/HMK(41), QE(40), FVCHLA(5), PE(5 ,41) COMMON/VOL/ZMAX,DZ(40),Z(40),A(40),V(40),TV(40),ATOP(41),DBL COMMON/STEPS/DZLL,DZUL,MBOT,NM,NPRINT,MDAY,MONTH,ILAY,DY COMMON/INFLOW/QIN(5),TIN(5),PAIN(5),BOOIN(5),DOIN(5),CIN(5), +CDIN(5),XNHIN(5),XNOIN(5),CHLAIN(3,5) COMMON/CHOICE/MODEL,NITRO,IPRNT(6),NDAYS,NPRNT(30),NCLASS,PLOT(30) COMMON/RESULT/ VAR(40,21) . CHARACTER*1 FF DIMENSION VOUT(11),IFL(11) DATA IFL/1,3,5,6,7,8,9,10,11,12,13/ IF(IPRNT(1).EQ.1) THEN FF=CHAR(12) WRITE(8,999) FF 999 FORMAT(1X,A1) WRITE(8,1000) MONTH,MDAY,MYEAR 1000 FORMAT(///,3X,I3,1X,I2,1X,I4) ENDIF C ... TAKE CARE OF INFLOW AND OUTFLOW DELZ=P-HED ZMAX=ZMAX+DELZ DZ(1)=DZ(1)+DELZ CALL SETZ(MBOT) V(1)=V(1)+DELZ*ATOP(1) C ... CALL TO LAKE FOR ADDITIONAL SURFACE INFLOWS AND OUTFLOWS CALL LAKE(O. ,0. ,NFLOW,10) IF(NFLOW.GT.O) THEN IF(IPRNT(1).EQ.1) WRITE(B,2235) 2235 FORMAT(//,5X,'INFORMATION ON INFLOW',/,5X,21('-'),/ 6X, + 'ISOPYCNIC',3X,'INITIAL VOLUME',7X,'INFLOW (M3/DAyj',BX,'P',7X, + 'N03', 7X, 'NH4' ,6X, 'DO' ,6X, 'BOD',/ ,6X, 'DEPTH (M)', 7X, + 'OF LAYER (M3)',4X,'INITIAL',8X,'FINAL',5X,5('(MG/L)',3X),/) C ... QIN > 0.0 SIGNIFIES INFLOW. QIN < 0.0 SIGNIFIES OUTFLOW. C ... THIS SECTION SKIPPED IF QIN = 0.0 DO 14 IW a 1,NFLOW DO 15 I=1,MBOT 15 QE(I)=O.O IF(QIN(IW).GT.O.O) THEN CALL DCFLOW(IPL,S,FT,WCHANL,IW) CALL CONSMAS(IW,IPL,QIN(IW» ENDIF IF(QIN(IW).LT.O.O) THEN CALL WDEPTH(ST,QIN(IW),LW) ENDIF DO 17 I=1,MBOT 17 V(I)-V(I)-QE(I) 1-1 C ... MERGE LOW VOLUME LAYERS WITH NEXT HIGHER LAYER 4 IF(V(I).LT.500.) THEN CALL MERGE(I,MBOT,ILAY) 94 79 IF(I.GT.MBOT) GOTO 300 BO GOTO 4 B1 ENDIF B2 1·1+1 B3 IF(I.GT.MBOT) GOTO 300 B4 GOTO 4 B5 300 CALL THICKNS(MBOT) B6 CALL SETZ(MBOT) B7 14 CONTINUE BB ENDIF B9 1= 1 90 C ... DETERMINE THAT ALL LAYERS ARE WITHIN THE LIMITS FOR 91 C ... MAXIMUM AND MINIMUM LAYER THICKNESS 92 5 IF(DZ(I).LT.DZLL) THEN 93 CALL MERGE(ItMBOT,ILAY) 94 IF(I.GT.MBOT, GO TO 301 95 GD TD 5 96 END IF 97 1"1+1 9B IF(I.GT.MBOT) GO TO 301 99 GO TO 5 100 301 IF(MBOT.GE.40) GO TO 2 101 1=1 102 3 IF(DZ(I).GT.DZUL .AND. MBOT.LT.40) THEN 103 CALL SPLIT(I,ILAY) 104 GO TO 3 105 END IF 106 1"1+1 107 IF(I.GT.MBOT .OR. I.GT.40) GO TO 2 10B GO TO 3 109 2 CALL SETZ(MBOT) 110 C ... DETERMINE LAKE STAGE FROM WATER BUDGET 111 ZMAX·Z(MBOT)+DZ(MBOT)*0.5 112 ST-DBL+ZMAX 113 CALL VOLUME(MBOT) 114 CALL AREA 115 CALL ABOUND 116 CALL TVOL(MBOT) 117 RETURN 11B END NAME TYPE DDZ DELZ FF FT HED IPL IW LW REAL REAL CHAR*1 REAL MAXIMUM DEPTH CONTRIBUTING TO SURFACE OUTFLOW CHANGE IN THICKNESS OF TOP LAYER DUE TO PRCIPITATION AND EVAPORATION SCREEN CONTROL VARIABLE MANNINGS ROUGHNESS COEFFICIENT FOR DENSITY CURRENTS LOSS FROM TOP LAYER DUE TO EVAPORATION (M/DAY) INDEX FDR THE LAYER RECEIVING THE DENSITY CURRENT INTERFLOW INDEX FOR INFLOw/oUTFLOW MY EAR NFLOW P REAL INTEGER*2 INTEGER*2 INTEGER*2 INTEGER*2 INTEGER*2 REAL INDEX FOR THE LOWEST LAYER CONTRIBUTING TO SURFACE OUTFLOW CURRENT YEAR OF SIMULATION NUMBER OF INFLOWS AND OUTFLDWS PRECIPITATION (M/DAY) S REAL SLOPE OF INFLOW CHANNEL ST WCHANL REAL REAL STAGE (M) WIDTH OF INFLOW CHANNEL (M) 119 SUBROUTINE AREA 120 C***** 121 C***** COMPUTE THE AREA THROUGH THE MIDDLE OF EACH LAYER 122 C***** USING THE DEPTH-AREA RELATIONSHIP IN LAKE 123 C***** 124 REAL*B A,V,TV,ATOP.AOUM 125 COMMON/sTEPs/DZLL,DZUL,MBOT,NM,NPRINT,MDAY,MONTH,ILAY.DY 126 COMMON/VOL/zMAX,DZ(40),Z(40),A(40),V(40),TV(40),ATOP(41),DBl 127 DUM=O. 12B DO 100 I-1,MBOT 129 ZDUM=ZMAX-DUM-DZ(I)/2. 130 DUM-DUM+DZ(I)· 131 CALL LAKE(ZDUM,ADUM,O,1) 132 100 A(I).ADUM 133 RETURN 134 END NAME TYPE ADUM DUM ZDUM REAL*B REAL REAL AREA FROM lAKE ROUTINE DUMMY VARIABLE FOR DEPTH FROM THE SURFACE DEPTH FROM THE BOTTOM USED IN LAKE 95 1 1 2 1 1 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 C***"'* C***** C***** C***** C***** C***** C***** C***** C***** SUBROUTINE CHLORO(K,TD,IEUPH) COMPUTE THE CHLOROPHVLL-A AND BOD CONCENTRATION PROFILES USING THE CONSTANT VOLUME FINITE DIFFERENCE METHOD (FROM PATANKAR) FOR MODELS 2 AND 3. MODEL 1 COMPUTES THE MIXED LAVER CHLOROPHVLL-A CONCENTRATION FOLLOWING WORK BV FORSBERG AND SHAPIRO. MODEL 2 USES MICHAELIS-MENTON GROWTH KINECTICS AND MODEL 3 USES A CELLULAR NUTRIENT GROWTH FORMULATION (FROM LEHMAN, ET. AL. ) REAL*8 A,V,TV,ATOP,AK,BK,CK,OK,BB COMMON/ZOOPL/IZ,MINDAV,MAXDAV,ZP,ZPMIN,PRMIN,PRMAX,PREDMIN,XIMIN, +XIMAX,XKRZP,GRAZMAX(3) ,THGRAZ(3) ,ASM,THRZP,HSCGRAZ(3) ,CHLAMIN(3), +REPRO,XI,XKMZ,GRAZE(3,40) COMMON/VIELD/VCA,VCH02,V2CH02,VCBOD,VPBOD,VZW,VPZP,VNZP,VZDO, +VSCHL,VNHBOD,BRNO,BRNH,XKNNH,THNNH,VPCHLA,BODK20,SB20,BRR COMMON/PHVTO/PDEL(3),PMAX(3),PMIN(3),THR(3),THM(3),XKR1(3), +XKR2(3),XKM(3),HSCPA(3),HSC1(3),HSC2(3),UPMAX(3),THUP(3), +GROMAX(3),TMAX(3),TOPT(3),XNMAX(3),XNMIN(3),UNMAX(3), THUN(3), +HSCN(3) ,HSCNH(3) ,XNDEL(3) ,IDIATOM,CHLMEAN,CHLMAX,SECC HI COMMON/TEMP/PARIO(4),PCDUM(3,40),XNHD(40),XNOD(40), + CHLAOUM(3,40),XNCD(3,40),PADUM(40),SID(40) COMMON/RESULT/ T2(40) ,CHLATOT(40),PA2(40) ,PTSUM(40) ,BOD2(40), +DS02(40) ,C2(40),CD2(40),XN02(40),XNH2(40) ,CHLA2(3,40) , +PC2(3,40),XNC2(3,40),T20(40),SI2(40) COMMON/VOL/ZMAX,DZ(40),Z(40),A(40),V(40),TV(40),ATOP(41),DBL COMMON/STEPS/DZLL,DZUL,MBOT,NM,NPRINT,MDAV,MONTH,ILAV,DV COMMON/CHOICE/MODEL,NITRO, IPRNT(6) ,NDAYS,NPRNT(30) ,NCL ASS,PLOT(30) COMMON/FLOW/HMK(41),QE(40),FVCHLA(6),PE(6,41) COMMON/SOURCE/RM(3,40) ,PROD(40) ,XMR(3,40) ,PRODSUM(40) COMMON/WATER/BETA,EMISS,XK1,XK2,HKMAX,WCOEF,WSTR,RHOW,HCAP COMMON/SO LV/ AK(40),BK(40),CK(40),DK(40) DIMENSION BB(40) ,PDUM(40) ,UN(40) + IF(MODEL.EQ.1) THEN DMIX=Z(ILAV)+DZ(ILAV)*0.5 IDEPTH=ILAV IF(ILAV.LT.IEUPH) IDEPTH=IEUPH DO 23 I=1,IDEPTH X=1.-PMIN(1)*CHLA2(1,I)/PTSUM(I) IF(X.LT.O) THEN PROD(I)·O.O ELSE PROD(I)=1.9*PMAX(1)*X/ «XK1+XK2*CHLA2(1,I»*DMIX*VCA) ENDIF 180 23 181 CONTINUE IF(IDEPTH.LT.MBOT) THEN 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 24 25 DO 24 I=IDEPTH+1,MBOT PROD (I )=0.0 ENDIF DO 25 I=1,MBOT CHLA2(1,I)=CHLA2(1,I)+CHLA2(1,I)*(PROD(I)-XKM(1» IF(CHLA2(1,I).LT.0.0001) CHLA2(1,I)=0.OO01 CHLATOT(I)=CHLA2(1,I) CHLADUM(1,I)=CHLA2(1,I) RETURN ENDIF DO 30 I=1,MBOT-1 UN (I )"1. DUM1=1./(A(I)*DZ(I» DUM2-ATOP(I)*FVCHLA(K)*DUM1 DUM3=ATOP(I+1)*FVCHLA(K)*DUM1 DUM4=PE(K, I) DUM5"PE(K, 1+1) IF(DUM4.LT.0.0) DUM4"0.0 IF(DUM5.LT.0.0) DUM5=0.0 AK(I)=-DUM2*(1.0+DUM4) BB(I)=1.0+DUM3*(1.0+DUM5)+DUM2*DUM4-FVCHLA(K)*DUM1*(ATOP(I+1) +-ATOP(I» CK(I)=-DUM3*DUM5 30 CONTINUE C-------------------------LOWER BOUNDRV CONOITION----------------- OUM1"1./(A(MBOT)*OZ(MBOT» DUM2"ATOP(MBOT)*FVCHLA(K)*DUM1 DUM4-PE(K,MBOT) IF(DUM4.LT.0.0) DUM4-0.0 AK(MBOT)=-DUM2*(1.0+DUM4) BB(MBOT)=1.0+DUM2*DUM4+FVCHLA(K)*DUM1*ATOP(MBOT) C------------------ COMPUTATION OF BOD -------------------- 102 101 IF(K.EQ.NCLASS+1) THEN ASStM=1.0-ASM IF(MODEL.NE.2) ASSIM=1.0 DO 101 I-1,MBOT CKRM=O. DO 102 K=1,NCLASS CKRM=CKRM+XMR(K,I)*CHLA2(K,I)+GRAZE(K,I)*ASSIM DK(I)-BOD2(I)+CKRM/YCBOD BK(I)=BB(I)+BODK20*1.047**T20(I) 96 223 IF(MOOEL.NE.l) DK(IZ)=DK(IZ)+XKMZ!VZDO*YZW*.OOl 224 CALL SOLVE(BOD2,MBOT) 225 RETURN 226 ENDIF 227 C************ MICHAELIS-MENTON MODEL **************** 228 IF(MODEL.EQ.2) THEN 229 DO 35 I=l,MBOT 230 IF(T2(I).LT.TOPT(1» THEN 231 XT=EXP(-2.3*«T2(I)-TOPT(1»!TOPT(1»**2) 232 ELSE 233 XT-EXP(-2.3*«T2(I)-TOPT(1»!(TMAX(1)-TOPT(1»**2» 234 ENDIF 235 X=PADUM(I)!(HSCPA(l)+PADUM(I» 236 PRDD(I)=AMIN1(PROD(I),X) 237 35 PROD(I)=GROMAX(1)*XT*PROD(I)*TD 238 ENDIF 239 C***************** INTERNAL NUTRIENT MODEL **************** 240 C-------------- COMPUTATION OF NITROGEN UPTAKE------------- 241 IF(NITRO.EQ.l) THEN 242 DO 110 IE1 MBOT 243 X=(XNHD(I)+XNOD(I»!(XNDEL(K)*(HSCN(K)+XNHD(I)+XNOD(I») 244 BK(I)=BB(I)+(RM(K,I)+XMR(K,I)+GRAZE(K,I)!CHLADUM(K,I) 245 + +UNMAX(K)*THUN(K)**T20(I)*X) 246 DK(I)=XNCD(K,I)+UNMAX(K)*THUN(K)**T20(I) 247 + *CHLADUM(K,I)*XNMAX(K)*X 248 110 CONTINUE 249 CALL SOLVE (PDUM,MBOT) 250 DO 200 I=1,IEUPH 251 XNC2(K,I)=PDUM(I) 252 IF(K.LT.3) UN(I) =1-XNMIN(K)*CHLADUM(K,I)!XNC2(K,I) 253 200 CONTINUE 254 DO 201 I=IEUPH+1,MBOT 255 201 XNC2(K,I)=PDUM(I) 256 ENDIF 257 IF(MODEL.EQ.3) THEN 258 C------------- COMPUTATION OF SILICA LIMITATION ON DIATOMS ------ 259 IF(K.EQ.1.AND.IDIATOM.GT.0) THEN 260 DO 140 I=1,IEUPH 261 UN(I)=SID(I)!(SID(I)+HSCSI) 262 140 CONTINUE 263 ENDIF 264 C---------------COMPUTATION OF PHOSPHORUS UPTAKE--------------- 265 DO 120 I=1,MBOT 266 X=PADUM(I)!(PDEL(K)*(HSCPA(K)+PADUM(I») 267 BK(I)=BB(I)+RM(K,I)+XMR(K,I)+GRAZE(K,I)!CHLADUM(K;I) 268 + +UPMAX(K)*THUP(K)**T20(I)*X 269 DK(I)=PCDUM(K.I)+UPMAX(K)*THUP(K)**T20(I)*X 270 + *CHLADUM(K,I)*PMAX(K) 271 120 CONTINUE 272 CALL SOLVE (PDUM,MBOT) 273 C-------------- COMPUTATION OF LIMITING NUTRIENT --------------- 274 DO 210 I=l,IEUPH 275 PC2(K,I)=PDUM(I) 276 UP=1.-PMIN(K)*CHLADUM(K,I)!PC2(K,I) 277 IF(T2(I);LT.TOPT(K» THEN 278 XT=EXP(-2.3*«T2(I)-TOPT(K»!TOPT(K»**2) 279 ELSE 280 XT=EXP(-2.3*«T2(I)-TOPT(K»!(TMAX(K)-TOPT(K»**2» 281 ENDIF 282 PROD(I)=AMIN1(PRDD(I),UN(I),UP) 283 IF(PROD(I).LT. 0.0) PROD(I)=O.O 284 PROD(I)=GROMAX(K)*XT*PROD(I)*TD 285 210 CONTINUE 286 DO 220 I=IEUPH+l,MBOT 287 220 PC2(K,I)-PDUM(I) 288 ENDIF 289 C-------------- COMPUTATION OF CHLOROPHYLL-A CONCENTRATION---- 290 DO 40 I=l,MBOT 291 BK(I)=BB(I)+XMR(K,I)+GRAZE(K,I)!CHLADUM(K,I) 292 DK(I)=CHLADUM(K,I)+PROD(I)*CHLADUM(K,I) 293 40 CONTINUE 294 CALL SOLVE(PDUM,MBDT) 295 DO 600 I=l,MBOT 296 CHLA2(K,I)=PDUM(I) 297 IF(CHLA2(K, I). LT .0.0001) CHLA2(K, I )=.0001 298 PRODSUM(I)=PRODSUM(I)+CHLA2(K,I)*PROD(I) 299 600 CONTINUE 300 RETURN 301 END NAME TYPE ASSIM BB CKRM DMIX DUM1 REAL REAL*9 REAL REAL REAL ZOOPLANKTON GRAZING ASSIMULATION COEFFICIENT FOR NUTRIENT RECYCLING DIAGONAL MATRIX COEFFICIENT WITHOUT SOURCE OR SINK TERMS VARIABLE FOR SUMMING NUTRIENT RECYCLING OVER ALL ALGAL CLASSES MIXED LAYER OEPTH (M) COEFFICIENTS USED TO DESCRITIZE THE ADVECTION-DIFFUSION EQUATION 97 DUM2 DUM3 DUM4 DUM5 HSCSI IDEPTH IEUPH K PDUM TO UN UP X XT REAL COEFFICIENTS USED TO DESCRITIZE THE ADVECTION-DIFFUSION REAL COEFFICIENTS USED TO DESCRITIZE THE ADVECTION-DIFFUSION REAL COEFFICIENTS USED TO DESCRITIZE THE ADVECTION-DIFFUSION REAL COEFFICIENTS USED TO DESCRITIZE THE ADVECTION-DIFFUSION REAL HALF-SATURATION COEFFICIENT FOR SILICA (MG/L) INTEGER*2 INDEX OF LAYER AT THE MIXED LAYER DEPTH INTEGER*2 INDEX OF THE LAYER AT THE EUPHOTIC DEPTH INTEGER*2 INDEX OF ALGAL CLASS REAL DUMMY VARIABLE USED BY TRIDIAGONAL MATRIX SOLVER REAL HOURS OF DAYLIGHT IN A GIVEN DAY REAL NITROGEN GROWTH LIMITATION COEFFICIENT REAL PHOSPHORUS GROWTH LIMITATION COEFFICIENT REAL TEMPORARY VARIABLE FOR GROWTH LIMITATION COEFFICIENTS REAL RATIO OF TEMPERATURE EFFECT ON GROWTH 302 SUBROUTINE COEF(MODEL,MBOT,NCLASS) 303 C***** 304 C***** COMPUTE SOME COEFFICIENTS USED IN THE CONSTANT 305 C***** VOLUME ANO FINITE DIFFERENCE SOLUTIONS 306 C***** 307 REAL*8 A,V,TV,ATOP 308 COMMON/COEFF/ DUM2(40),DUM3(40) EQUATION EQUATION EQUATION EQUATION 309 COMMON/VOL/ ZMAX,DZ(40),Z(40),A(40),V(40),TV(40),ATOP(41),DBL 310 COMMON/FLOW/ HMK(41),QE(40),FVCHLA(5),PE(5,41) 311 DO 100 I=2,MBOT-1 312 DUM1 0 2./(A(I)*DZ(I» 313 DUM2(I)=DUM1*ATOP(I)*HMK(I)/(DZ(I)+DZ(I-1» 314 100 DUM3(I)=DUM1*ATOP(I+1)*HMK(I+1)/(DZ(I)+DZ(I+1» 315 KK=2 316 IF(MODEL.GT.1) KK a 1 317 DO 200 K=KK,NCLASS+1 1 318 DO 200 I=2,MBOT 2 319 X=FVCHLA(K)*(DZ(I-1)+DZ(I»*.5/HMK(I) 2 320 C---------- PE=(1.0-.1*ABS(X»**5/X ----------------------------2 321 AO=1.0-.1*ABS(X) 2 322 A1=AO*AO 2 323 PE(K,I)=A1*A1*AO/X 2 324 PE(K,1)=0.0 2 325 200 PE(K,MBOT+1)=0.0 326 RETURN 327 END NAME TYPE AO A1 I K REAL DUMMY VARIABLE FOR POWER LAW COMPUTATIONS DUMMY VARIABLE FOR AO**2 LAYER INDEX ALGAL CLASS INDEX KK MBOT MODEL NCLASS X REAL INTEGER*2 INTEGER*2 INTEGER*2 INTEGER*2 INTEGER*2 INTEGER*2 REAL DUMMY VARIABLE TO COMPUTE POWER LAW FOR BOD ONLY IF MODEL=1 NUMBER OF LAYERS IN MODEL BIOLOGICAL SIMULATION OPTION NUMBER OF ALGAL CLASSES DUMMY VARIABLE FOR PECLET NUMBER 328 SUBROUTINE CONMIX(ILAY,TMIX,MBOT) 329 C***** 330 C***** REMOVE DENSITY INSTABILITIES BY MIXING UNSTABLE 331 c***** LAYERS DOWNWARD ANO MERGING WITH LOWER LAYERS. 332 C***** 333 REAL*8 A,V,TV,ATDP 334 COMMON/RESULT/ T2(40) ,CHLATOT(40), PA2(40) ,PTSUM(40),BOD2(40), 335 +DS02(40) tC2(40) ,CD2(40) ,XN02(40) ,XNH2(40) ,CHLA2(3,40) , 336 +PC2(3,40} ,XNC2(3,40) ,T20(40) ,SI2(40) 337 COMMON/VOL/ZMAX,DZ(40) ,Z(40) ,A(40),V(40),TV(40),ATOP(4 1),DBL 338 DIMENSION RHOT(40) 339 DO 100 I=1,MaOT 340 100 RHOT(I)=RHO(T2(I),0.,0.) 341 6 IFLAG=O 342 1=0 343 M=MBOT-1 344 I-I+1 345 IF(I.EQ.M) GO TO 5 346 IF(RHOT(I).LE.RHOT(I+1» GO TO 1 347 IFLAG=1 348 IB-I 349 TVDUM=T2(I)*V(I) 350 VDUM-V(I) 351 TDUM=TVDUM/VDUM 352 RHODUM=RHO(TDUM,O.,O.) 353 ~"I-1 354 3 ~=I.I+1 355 IF(RHODUM.LE.RHOT(I.I+1» GO TO 2 356 IFLAG-1 98 ., 357 TVDUM~TVDUM+T2(~+1)*V(~+1) 358 VDUM·VDUM+V(~+1) 359 TDUM"TVDUM/VDUM 360 RHODUM=RHO(TDUM,O.,O.) 361 IF(J.EQ.M) GO TO 2 362 GO TO 3 363 2 IE"~ 364 IF~J.EQ.M) IE=IE+1 365 IF IB.EQ.IE) GO TO 4 366 DO 200 K-IB, IE 367 T2(K)=TDUM 368 200 RHOT(K)=RHODUM 369 4 IF(I.NE.M) GO TO 1 370 5 IF(IFLAG.EQ.1) GO TO 6 371 C ... OETERMINE MIXEO LAYER OEPTH ... 372 DO 700 I-1,MBOT-1 373 IF«T2(I)-T2(I+1».LE .. OO1) GO TO 700 374 ILAY=I 375 GO TO 10 376 700 CONTINUE . 377 10 TMIX"T2( 1) 378 RETURN 379 END NAME TYPE I LAYER INDEX IB IE IFLAG ILAY INTEGER*2 INTEGER*2 INTEGER*2 INTEGER*2 INTEGER*2 INTEGER*2 INTEGER*2 INTEGER*2 REAL INDEX OF UNSTABLE LAYER INDEX FOR STABLY MIXED LAYER FLAG FOR PASS THROUGH ALL LAYERS WITHOUT FINDING AN INSTABILITY INDEX OF LAYER AT THE MIXED LAYER DEPTH J M INDEX OF LOWER LAYER MIXED WITH UNSTABLE LAYERS INDEX FOR LAYER MBOT-1 NUMBER OF LAYERS IN THE MODEL MBOT RHODUM RHOT TDUM TMIX TVDUM VDUM REAL DUMMY VARIABLE FOR DENSITY OF TWO MIXED LAYERS ARRAY OF DENSITY IN EACH LAYER REAL REAL REAL REAL DUMMY TEMPERATURE IN TWO MIXED LAYERS TEMPERATURE IN THE MIXED LAYER SUM OF PRODUCTS OF T2*V IN LAYERS COMBINEO TO REMOVE INSTABILITIES DUMMY VARIABLE FOR VOLUME OF LAYERS COMBINED TO REMOVE INSTABILITIES SUBROUTINE CONSMAS(IW,IPL,DCF) 380 381 382 383 38-4 385 386 387 388 c***** c***** C***"'* C*"'*"'* C*"'''''''* C"'**** COMPUTE THE FINAL CONCENTRATION IN A LAVER USING THE CONCENTRATION IN THE DENSITY CURRENT AND THE CONCENTRATION IN THE LAYER WEIGHTED BY THE VOLUME IN THE LAYER AND THE VOLUME IN THE DENSITY CURRENT 389 390 391 392 393 394 395 396 397 398 399 400 401 402 50 403 404 405 406 407 408 409 410 411 412 REAL*8 A,V,TV,ATOP COMMON/CHOICE/MODEL,NITRO, IPRNT(6) ,NDAYS,NPRNT(30) tNCL ASS,PLOT(30) COMMON/RESULT/ T2(40),CHLATOT(40),PA2(40),PTSUM(40J,BOD2(40), +DS02(40) ,C2(40) ,CD2(40) ,XN02(40) ,XNH2(40),CHLA2(3,40) , +PC2(3,40),XNC2(3,40),T20(40),SI2(40) COMMON/STEPS/DZLL,DZUL,MBOT,NM,NPRINT,MDAY,MONTH,ILAY,DY COMMON!FLOW/HMK(41),QE(40),FVCHLA(5),PE(5,41) COMMON/VOL/ZMAX,DZ(40),Z(40),A(40),V(40),TV(40),ATOP(41),DBL COMMON/INFLOW/QIN(5),TIN(5),PAIN(5),BODIN(5),DOIN(5),CIN(5), +CDIN(5),XNHIN(5),XNOIN(5),CHLAIN(3,5) VRCP-1./(V(IPL)+DCF) T2(IPL)=(T2(IPL)*V(IPL)+TIN(IW)*DCF)"'VRCP C2(IPL)=(C2(IPL)*V(IPL)+CIN(IW)*DCF)"'VRCP CD2(IPL)-(C02(IPL)*V(IPL)+CDIN(IW)*DCF)*VRCP DO 50 K= 1, NCLASS CHLA2(K,IPL)=(CHLA2(K,IPL)*V(IPL)+CHLAIN(K,IW)*DCF)*VRCP PA2(IPL)a(PA2(IPL)"'V(IPL)+PAIN(IW)*DCF)*VRCP BOD2(IPL)=(80D2(IPL)*V(IPL)+BODIN(IW)*DCF)*VRCP DS02(IPL)-(DS02(IPL)*V(IPL)+DOIN(IW)*DCF)*VRCP IF(NITRO.EQ.1) THEN XNH2(IPL)-(XNH2(IPL)*V(IPL)+XNHIN(IW)*DCF)"'VRCP XN02(IPL)=(XN02(IPL)"'V(IPL)+XNOIN(IW)"'OCF)"'VRCP ENDIF V(IPL)-V(IPL)+DCF RETURN END NAME TYPE DCF IPL IW K VRCP REAL INTEGER*2 INTEGER*2 INTEGER*2 REAL VOLUME OF WATER IN THE DENSITY CURRENT (M3) INDEX FOR LAYER AT THE INTERFLOW DEPTH INDEX FOR INFLOW SOURCE FROM THE INFLOW DATA FILE INDEX FOR ALGAL CLASS DUMMY VARIABLE FOR RECIPROCAL OF THE SUM V(IPL) AND DCF 99 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1 1 1 1 1 1 1 413 SUBROUTINE DCFLOW(IPL,S,FT,WIDTH,IW) 414 C***** 415 C***** DETERMINE THE DENSITY CURRENT VOLUME AND CONCENTRATION 416 C***** WITH ENTRAINMENT FROM EACH LAYER PASSED AND THE 417 C***** ISOPYCNIC LAYER THAT RECEIVES THE DENSITY CURRENT. 418 C***** 419 REAL*8 A,V,TV,ATOP 420 COMMON/FILE/ DIN,MET,FLO,TAPE8,TAPE1,IREC 421 COMMON/VOL/ZMAX,DZ(40),Z(40),A(40),V(40),TV(40),ATOP(41),DBL 422 COMMON/RESULT/ T2(40),CHLATOT(40),PA2(40),PTSUM(40),BOD2(40). 423 +DS02(40).C2(40).CD2(40).XN02(40).XNH2(40),CHLA2(3.40), 424 +PC2(3,40).XNC2(3,40).T20(40).SI2(40) 425 COMMON/FLOW/HMK(41),QE(40),FVCHLA(5).PE(5,41) 426 COMMON/CHOICE/MODEL,NITRO,IPRNT(6).NDAYS,NPRNT(30).NCLASS.PLOT(30) 427 COMMON/STEPS/DZLL.DZUL,MBOT,NM.NPRINT.MDAY,MONTH.ILAY,DY 428 COMMON/INFLOW/QIN(5),TIN(5),PAIN(5),BODIN(5),DOIN(5).CIN(5), 429 +CDIN(5),XNHIN(5),XNOIN(5),CHLAIN(3,5) 430 CHARACTER*16 DIN.MET,FLO.TAPE8,TAPE1 431 X-QIN(IW) 432 RHOMIX=RHO(T2(1),C2(1).CD2(1» 433 RHOIN=RHO(TIN(IW).CIN(IW).CDIN(IW» 434 CALL PDEPTH(RHOMIX,RHOIN,QIN(IW),S.IHP,QENIN.MBOT,SUMZ, 435 +FT,WIDTH.HP) 436 DO 100 I=1,MBOT-1 437 C ... COMPUTE THE DENSITY OF INFLOW AS IT ENTERS EACH LAYER 438 RDC=RHO(TIN(IW),CIN(IW),CDIN(IW» 439 RHOAMB=RHO(T2(I),C2(I).CD2(I» 440 C ... LOCATE THE ISOPYCNIC LAYER 441 IF(RDC.GT.RHOAMB) GO TO 5 442 IPL=I-1 443 IF(IPL.LT.1) IPL=1 444 GO TO 3 445 C ... CONTINUE THE SEARCH 446 C ... COMPUTE THE ENTRAINMENT FROM EACH LAYER 447 5 QE(I)=ENTRAIN(I.QIN(IW),RDC,RHOAMB,DZ(I),S,SUMZ,QENIN.IHP. 448 + WIDTH.FT) 449 IF(QE(I).GT.V(I» QE(I)-V(I) 450 RDCF D 1./(QIN(IW)+QE(I» 451 TIN(IW)=(TIN(IW)*QIN(IW)+T2(I)*QE(I»*RDCF 452 CIN(IW)=(CIN(IW)*QIN(IW)+C2(I)*QE(I»*RDCF 453 CDIN(IW)=(CDIN(IW)*QIN(IW)+CD2(I)*QE(I»*RDCF 454 DO 50 K=1,NCLASS 455 50 CHLAIN(K,IW)=(CHLAIN(K,IW)*QIN(IW)+CHLA2(K,I)*QE(I»*RDCF 456 PAIN(IW)=(PAIN(IW)*QIN(IW)+PA2(I)*QE(I»*RDCF 457 BODIN(IW)=(BODIN(IW)*QIN(IW)+BOD2(I)*QE(I»*RDCF 458 DOIN(IW)-(DOIN(IW)*OIN(IW)+DS02(I)*QE(I»*RDCF 459 IF(NITRO.EQ.1) THEN 460 XNHIN(IW)=(XNHIN(IW)*OIN(IW)+XNH2(I)*QE(I»*RDCF 461 XNOIN(IW)=(XNOIN(IW)*QIN(IW)+XN02(I)*QE(I»*RDCF 462 ENDIF 453 QIN(IW)QQIN(IW)+QE(I) 464 100 CONTINUE 465 RDC=RHO(TIN(IW),CIN(IW),CDIN(IW» 466 IPL"MBDT . 467 3 IF(IPRNT(1).EQ.1) WRITE(8,2oo0)Z(IPL)t V(IPL),X,QIN(IW),PAIN(IW), 468 +XNOIN(IW),XNHIN(IW),DOIN(IW),BODIN(IW} 469 2000 FORNAT(7X,F5.2,7X,E11.5,4X,E11.5,3X,E11.5,2X,3(F6.3,3X), 470 +2(F6.2,3X» 471 RETURN 472 END NAME TYPE ENTR FT HP I IHP IPL IW K QENIN ROC RDeF RHOAMB RHOIN RHOMIX S SUMZ WIDTH X REAL REAL REAL INTEGER*2 INTEGER*2 INTEGER*2 INTEGER*2 INTEGER*2 REAL REAL REAL REAL REAL REAL REAL REAL REAL REAL ENTRAINED VOLUME OF WATER FROM A LAYER PASSED BY THE DENSITY CURRENT MANNINGS ROUGHNESS FOR THE INFLOW CHANNEL DEPTH AT PLUNGE POINT (M) INDEX FOR LAYERS INDEX FOR DEEPEST LAYER AT PLUNGE POINT INDEX FOR LAYER RECEIVING THE INTERFLOW INDEX FOR' INFLOW SOURCE FROM INFLOW DATA FILE INDEX FOR ALGAL CLASS INITIAL ENTRAINMENT AT PLUNGING DENSITY OF DENSITY CURRENT RECIPROCAL OF QIN(IW) AND ENTR AMBIENT DENSITY IN A LAYER INITIAL DENSITY OF INFLOW DENSITY IN THE MIXED LAYER SLOPE OF INFLOW CHANNEL DEPTH AT PLUNGE POINT WIDTH OF INFLOW CHANNEL DUMMY VARIABLE FOR VOLUME OF INFLOW BEFORE PLUNGING 473 SUBROUTINE DISOLID (VAR2,WIND,IM,TD) 100 1 1 2 2 2 2 2 2 2 1 1 1 1 1 . 2 1 1 1 1 474 C***** 475 C**"'** 476 C***** 477 C***** 478 C* ... * ...... COMPUTE THE CONCENTRATION OF DISSOLVED SUBSTANCES BY SOLVING THE ONE DIMENSIONAL DIFFUSION EQUATION BY THE FULLV IMPLICIT CENTRAL DIFFERENCE SCHEME. 479 REAL*8 A,V,TV,ATOP,AK,BK,CK,DK 480 DIMENSION VAR2(40) 481 COMMON/SOLV/ AK( 40) , BK( 40) • CK( 40) , DK( 40) 482 COMMON/ZOOPL/IZ,MINDAV,MAXDAV,ZP,ZPMIN,PRMIN,PRMAX,PREDMIN,XIMIN, 483 +XIMAX,XKRZP,GRAZMAX(3),THGRAZ(3),ASM,THRZP,HSCGRAZ(S),CHLAMIN(3), 484 +REPRO,XI,XKMZ,GRAZE(S,40) 485 COMMON/VIELD/VCA,VCH02,V2CH02,VCBOD,VPBOD,VZW,VPZP,YNZP,YZOO, 486 +VSCHL,VNHBOD,BRNO,BRNH,XKNNH,THNNH,VPCHLA,BODK20,SB20,BRR 487 COMMON/TEMP/PARIO(4),PCOUM(S,40),XNHD(40),XNOO(40), 488 + CHLADUM(S,40),XNCD(S,40),PADUM(40),SID(40) 489 COMMON/VOL/ZMAX,OZ(40),Z(40),A(40),V(40),TV(40),ATOP(41),DBL 490 COMMON/RESULT/ T2(40),CHLATOT(40),PA2(40),PTSUM(40),BOD2(40), 491 +DS02(40),C2(40),CD2(40),XN02(40),XNH2(40),CHLA2(S,40), 492 +PC2(S,40),XNC2(3,40),T20(40),SI2(40) 49S COMMON/SOURCE/RM(3,40),PROD(40),XMR(S,40),PRODSUM(40) 494 COMMON/PHYTO/PDEL(3),PMAX(S),PMIN(S),THR(S),THM(3).XKR1(S). 495 +XKR2(3).XKM(S).HSCPA(S).HSC1(S),HSC2(S),UPMAX(S).THUP(S). 496 +GROMAX(3).TMAX(3).TOPT(S).XNMAX(S),XNMIN(S),UNMAX(S).THUN(S). 497 +HSCN(3).HSCNH(S).XNDEL(3).IDIATOM.CHLMEAN.CHLMAX,SECCHI 498 COMMON/STEPS/DZLL.DZUL.MBOT.NM.NPRINT.MOAV.MONTH,ILAV.DY 499 COMMON/FLOW/HMK(41).QE(40).FVCHLA(5),PE(5.41) 500 COMMON/COEFF/DUM2(40).DUMS(40) 501 COMMON/CHOICE/MODEL,NITRO. IPRNT(6) ,NDAVS.NPRNT(SO).NCL ASS.PLOT(SO) 502 C ... THIS SUBROUTINE COMPUTES CONCENTRATION OF OISSOLVED SOLIDS IN THE 503 C ... RESERVOIR BV SOLVING THE ONE DIMENSIONAL DIFFUSION EQUATION BV 504 C ... FULLY-IMPLICIT-CENTRAL DIFFERENCE SCHEME 505 DO 100 I=2.MBOT-1 506 AK(I)=-DUM2(I) 507 BK(I)=1.0+DUM2(I)+DUM3(I) 508 CK(I)=-DUMS(I) 509 DK(I)=VAR2(I) 510 100 CONTINUE 511 DK(1)=VAR2(1) 512 DK(MBOT)=VAR2(MBOT) 51S AK(1)=O. 514 CK(1)=-2.*ATOP(2)*HMK(2)/(A(1) ... DZ(1)*(DZ(1)+DZ(2») 515 BK(1)=1.-CK(1) 516 I=MBOT 517 AK(I)=-2.*ATOP(I)*HMK(I)/(A(I)*DZ(I)*(DZ(I)+DZ(I-1») 518 BK(I)=1.-AK(I) 519 CK(I)=O.O 520 GDTO (550,2000. SOOO. 4000.5000.6000,7000) 1M 521 C------------------COMPUTATION OF SILICA CONC.-------------- 522 7000 DO 115 I=1,MBOT 52S DK(I)=DK(I)+CHLADUM(1.I)*(RM(1,I)+XMR(1,I)-PROD(I»*YSCHL 524 115 CONTINUE 525 GOTO 550 526 C----------------- COMPUTATION OF AVAILABLE PHOSPHORUS--------- 527 2000 IF(MODEL.EQ.S) THEN 528 DO 150 1=1,MBOT 529 CKRM=O. 5S0 DO 151 K=1,NCLASS 5S1 CKRM=CKRM+RM( K. I ) *PC2 (K. I) 5S2 + +XMR(K,I)*(PC2(K.I)-PMIN(K)*CHLA2(K,I» 5SS + +GRAZE(K.I)*(PC2(K.I)/CHLADUM(K,I)-PMIN(K»*(1.-ASM) 5S4 BK(I)=BK(I)+UPMAX(K)*THUP(K)**T20(I)* 5S5 + (CHLADUM(K,I)"'PMAX(K)-PC2(K.I»/ 5S6 + (PDEL(K)"'(HSCPA(K)+PADUM(I») 5S7 151 CONTINUE 5S8 150 DK(I)=DK(I)+BRR"'(ATOP(I)-ATOP(I+1»/V(I)+CKRM+ 5S9 + BODK20"'1.047**T20(I)"'BOD2(I)"'YPBOD 540 ELSE 54 1 DO 160 I = 1 • MBOT 542 160 DK(I)=DK(I)+RM(1.I)*CHLADUM(1.I) 54S + +BRR*(ATOP(I)-ATOP(I+1»/V(I) 544 + +BODK20"'1.047**T20(I)*BOD2(I)*YPBOD 545 + -PROD(I)"'CHLADUM(1.I)*YPCHLA 546 ENDIF 547 CALL LAKE(0.,O.,O,6) 548 GOTO 550 549 SOOO RETURN 550 C-------------COMPUTATION OF DISSOLVED OXYGEN -------------- 551 4000 DOKO=(0.641+0.128*(WIND*.447)**2)/DZ(1) 552 DO 10S I=1,MBOT 55S X=O.O 554 BDDK2=BODK20*1.047*"'T20(I) 555 DO 104 K=1.NCLASS 556 104 X=X+RM(K.I)"'CHLADUM(K,I) 557 DK(I)=DK(I)-BODK2"'BOD2(I)-X/Y2CH02+PRODSUM(I)/VCH02 558 + -SB20"".065**T20(I)*(ATOP(I)-ATOP(I+1»/V(I) 559 IF(NITRO.EQ.1) BK(I)-BK(I)+XKNNH*THNNH**T20(I)*XNH2(I)*S.556 560 103 CONTINUE 561 DK(1)-DK(1)+DOKO*(14.652-T2(1)* 101 562 +(.41022-T2(1)*«7.991E-3)-7.7774E-5*T2(1»» 563 BK(1)=BK(1)+DOKO 564 IF(MODEL.GT.1) 565 + DK(IZ)=DK(IZ)-TD*.04167*XKRZP*THRZP**T20(IZ)*ZP*.001/V(IZ) 566 CALL LAKE(0.,0.,0,9) 567 GOTO 550 568 C-------------- COMPUTATION OF NITRATE-NITRITE ------------ 569 5000 DO 170 I=1,MBOT 1 570 CKRM~O. 1 571 X=XNHD(I)+XNOD(I) 1 572 DO 171 K=1,NCLASS 2 573 CKRM=CKRM-UNMAX(K)*THUN(K)**T20(I)*(1.-XNH2(I)/(HSCNH(K) 2 574 + +XNHD(I»)*(CHLADUM(K,I)*XNMAX(K)-XNCD(K,I»*X/(XNDEL(K) 2 575 + *(HSCN(K)+X» 2 576 171 CONTINUE 1 577 DK(I)=DK(I)+(CKRM*4.4286+XKNNH*THNNH**T20(I)*XNHD(I) 1 578 + +BRNO*(ATOP(I)-ATOP(I+1»/V(I» 1 579 170 CONTINUE 580 CALL LAKE(O.,O. ,0,7) 581 GOTO 550 582 C-------------- COMPUTATION OF AMMONIUM ------------------- 583 6000 DO 180 I=1,MBOT 1 584 CKRM=O. 1 585 X=XNHD(I)+XNOD(I) 1 586 DO 181 K=1,NCLASS 2 587 CKRM=CKRM+XNC2(K,I)*RM(K,I) 2 588 + +XMR(K,I)*(XNC2(K,I)-XNMIN(K)*CHLADUM(K,I» 2 589 + +GRAZE(K,I)*(XNC2(K,I)/CHLA2(K,I)-XNMIN(K»*(1.-ASM) 2 590 BK(I)=BK(I)+UNMAX(K)*THUN(K)**T20(I)/«HSCNH(K)+XNHD(I»* 2 591 + (XNDEL(K)*(HSCN(K)+X»)*(CHLADUM(K,I)*XNMAX(K)-XNC2(K,I» 2 592 + *X*1.286 2 593 181 CONTINUE 1 594 DK(I)=DK(I)+(CKRM*1.286+BRNH*(ATOP(I)-ATOP(I+1»/V(I) 1 595 + +BODK20*1.047**T20(I)*BOD2(I)*YNHBOD) 1 596 BK(I)=BK(I)+XKNNH*THNNH**T20(I) 1 597 180 CONTINUE NAME 598 DK(IZ)=DK(IZ)+TD*.04167*XKRZP*THRZP**T20(IZ)*ZP*YNH~P*.001*YZW 599 CALL LAKE(O. ,0.,0,8) 600 GOTO 550 601 C******************CALL SOLVE SUBROUTINE**************************** 602 550 CALL SOLVE(VAR2,MBOT) 603 DO 500 I=1,MBOT 604 IF(VAR2(I).LT.0.) VAR2(I)=0. 605 500 CONTINUE 606 RETURN 607 END TYPE BODK2 REAL TEMPERATURE CORRECTED DETRITAL DECAY RATE (1/DAY) SUM OF SOURCES AND SINKS OVER SEVERAL ALGAL CLASSES OXYGEN TRANSFER RATE THROUGH THE LAKE SURFACE CKRM DOKO I 1M K TD VAR2 WIND X REAL REAL INTEGER*2 INTEGER*2 INTEGER*2 REAL REAL REAL REAL INDEX OF LAYERS INDEX OF WHICH DISSOLVED SUBSTANCE IS BEING COMPUTED INDEX OF ALGAL CLASSES LENGTH OF DAYLIGHT PERIOD (HOURS) DUMMY VARIABLE FOR DISSOLVED SUBSTANCE BEING COMPUTED WIND SPEED FROM METEOROLOGICAL DATA FILE (MPH) DUMMY VARIABLE YNHZP REAL MASS YIELD RATIO OF AMMONIUM TO ZOOPLANKTON 608 FUNCTION ENTRAIN(I,DCF,RDC,RHOAMB,DELZ,S,SUMZ,Q,IHP,WIDTH,FT) 609 C***** 610 C***** COMPUTE THE ENTRAINMENT FROM A LAYER INTO THE 611 C***** DENSITY CURRENT (FROM AKIYAMA) 612 C***** 613 IF(I.GT.IHP) THEN 614 EPSI=(RDC-RHOAMB)/RHOAMB 615 IF(I.NE.IHP+1) GO TO 3 616 FD=1.875E-4+FT 617 F43=( (FD+SQRT(FD*FD+0.OO45*S»/( 1.5*S) )**( 1.3333) 618 3 X=DCF/(86400*WIDTH) 619 GAMAI=0.0015*DELZ*(9.81*EPSI/(X*X»**(.3333)/(F43*S) 620 ENTRAIN=GAMAI*DCF 621 ELSE 622 ENTRAIN=Q*DELZ/SUMZ 623 ENDIF 624 RETURN 625 . END 626 SUBROUTINE EUPHZ(ZEUPH,MBOT,IEUPH,RAD,ALBEDO,NCLASS) 102 · 627 C***** 628 C***** DETERMINE THE LIMIT OF THE EUPHOTIC ZONE GIVEN AS 629 C***** THE LAYER IN WHICH THE LIGHT INTENSITY IS LESS THAN 630 C***** ONE PERCNET OF THE SURFACE VALUE. ALSO COMPUTES THE 631 C***** LIGHT INTENSITY AT THE ZOOPLANKTON DAY DEPTH. 632 C***** 633 REAL*8 A,V,TV,ATDP 634 COMMDN/TEMP/PARIO(4),PCDUM(S,40) ,XNHD(40),XNOD(40) , 635 + CHLADUM(S,40),XNCD(3,40),PADUM(40) SID(40) 636 COMMDN/PHYTO/PDEL(S),PMAX(S),PMIN(S),THR(S),THM(S),XKR1(S), 637 +XKR2(3),XKM(S),HSCPA(S),HSC1(3),HSC2(3),UPMAX(3),THUP(S), 638 +GROMAX(3),TMAX(3) ,TOPT(3) ,XNMAX(3) ,XNMIN(S) ,UNMAX(S) , THUN(S), 639 +HSCN(3),HSCNH(3),XNDEL(S),IDIATOM,CHLMEAN,CHLMAX,SECCHI 640 COMMON/ZODPL/IZ,MINDAY,MAXDAY,ZP,ZPMIN,PRMIN,PRMAX PREDMIN,XIMIN, 641 +XIMAX,XKRZP,GRAZMAX(S),THGRAZ(S),ASM,THRZP,HSCGRAZlS),CHLAMIN(S), 642 +REPRO,XI,XKMZ,GRAZE(3,40) 643 COMMON!VOL/ZMAX,DZ(40)'Z(40)'A(40)'V(40)~TV(40)'ATOP(41),DBL 644 COMMON/RESULT/ T2(40),CHLATOT(40),PA2(40 ,PTSUM(40),BOD2(40), 645 +DS02(40),C2(40),CD2(40),XN02(40),XNH2(40 ,CHLA2(S,40), 646 +PC2(3,40),XNC2(3,40),T20(40),SI2(40) 647 COMMON!SOURCE/RM(3,40) ,PROD(40) ,XMR(3,40),PRDDSUM(40) 648 COMMON/WATER/BETA,EMISS,XK1,XK2 t HKMAX,WCOEF,WSTR,RHOW,HCAP 649 EK=XK1+.04S*C2(1)+XK2*CHLATOT(1, 650 EX=(1.-ALBEOD)*EXP(-EK*(Z(1)+DZ(1)/2.» 651 IF(IZ.EQ.1) XI=RAD 652 IF(EX.GT.O.01) GO TO 2 653 ZEUPH=Z(1)+OZ(1)*0.5 654 IEUPH-1 655 GOTO 101 656 2 00 100 I=2,MBOT 657 C ..... ATTENUATION COEFF. (EK) ..... 658 EK=XK1+.04S*C2(I)+XK2*CHLATOT(I) 659 C ..... DETERMINATION OF EUPHOTIC DEPTH (ZEUPH) 660 IF(I.EQ.IZ) XI=RAD*EX 661 EX-EX*EXP(-EK*DZ(I» 662 IF(EX.LE.O.01) THEN 663 ZEUPH=Z(I)+DZ(I)*0.5 664 IEUPH=I 665 GOTO 101 666 END IF 667 100 CONTINUE 668 ZEUPH=ZMAX 669 I EUPH=MBOT 670 101 IF(IZ.GT.IEUPH) THEN 671 00 102 I=IEUPH+1,IZ-1 672 EK-XK1+.4S*C2(I)+XK2*CHLATOT(I) 673 EX=EX*EXP(-EK*DZ(I» 674 102 CONTINUE 675 XI=RAD*EX 676 ENDIF 677 00 110 I=1,IEUPH 1 678 PRODSUM(I)-O.O 1 679 00 110 K=1,NCLASS 2 680 .110 RM(K,I)=XKR1(K)*THR(K)**T20(I) 681 IF(IEUPH.LT.MBDT) THEN 682 00 120 I=IEUPH+1,MBOT 1 683 PRODSUM(I)=O.O 1 684 00 120 K=1,NCLASS 2 685 120 RM(K,I)-XKR2(K)*THR(K)**T20(I) 686 ENDIF 687 RETURN 688 END NAME TYPE ALBEDO REAL REFLECTED RADIATION AT THE WATER SURFACE EK REAL EXTINCTION COEFFICIENT (1/M) EX REAL FRACTION OF LIGHT REMAINING AT A GIVEN LAYER I INTEGER*2 INDEX FOR LAYERS IEUPH INTEGER*2 INDEX FOR LAYER AT THE EUPHOTIC DEPTH K INTEGER*2 INDEX FOR ALGAL CLASS MBOT INTEGER*2 NUMBER OF LAYERS NCLASS INTEGER*2 NUMBER OF ALGAL CLASS~S RAD REAL RADIATION AT THE LAKE SURFACE (LANGLEY) ZEUPH REAL EUPHOTIC DEPTH (M) 689 FUNCTION FV(T) 690 C ... CALCULATES FALL VELOCITY OF SEDIMENT USING AN EMPIRICAL FORMULA 691 C ... DUE TO GIBBS ETAL FOR PARTICLE SIZES 0.1 TO 6000 MICRONS 692 C ... FALL VELOCITY IS IN M/DAY 69S 694 FV-.0179/(1.0+0.0S368*T+O.000221*T*T) 695 FV=(-S.*FV+SQRT(9.*FV*FV+.62666E-7»*74390. 696 103 1 1 1 1 1 2 2 2 2 1 1 1 2 2 2 1 1 2 1 1 1 1 2 697 RETURN 698 END 699 SUBROUTINE INTRPOL(IR,MBOT) 700 C***** 701 C***** APPLY RELAXATION FACTOR TO SUCCESSIVE ITERATIONS THROUGH 702 C***** THE BIOLOGICAL-NUTRIENT ROUTINES, IN THE FIRST PASS (IR=1), 703 C***** THE PREVIOUS TIME STEP VALUES FOR ALL VARIABLES ARE USED. 704 C***** IN SUCCESSIVE PASSES (IR>1), THE PREVIOUS AND PRESENT 705 C***** VALUES ARE AVERAGED. 706 C***** 707 COMMON/CHOICE/MODEL'NITRO'IPRNT(6)'NDAYS~NPRNT(30)tNCLASS,PLOr(30) 708 COMMON/RESULT/ T2(40),CHLATOT(40),PA2(40 ,PTSUM(40},BOD2(40), 709 +DS02(40}tC2(40},CD2(40),XN02(40),XNH2(40 ,CHLA2(3,40}, 710 +PC2(3,40} ,XNC2(3,40},T20(40) ,SI2(40) 711 COMMON/TEMP/PARIO(4) ,PCDUM(3,40) ,XNHD(40},XNOD(40) , 712 + CHLADUM(3,40),XNCD(3,40},PADUM(40),SID(40) 713 XX-0.5 714 IF(IR.EQ.1) THEN 715 DO 120 I=1 t MBOT 716 PADUM(I}=O.O 717 XNHD(I)=O.O 718 XNOD(I)=O.O 719 SID(I)=O.O 720 DO 120 K=1,NCLASS 721 PCDUM(K,I)=O.O 722 CHLADUM(K,I)=O.O 723 XNCD(K,I)=O.O 724 120 CONTINUE 725 XX=1.0 726 ENDIF 727 DO 100 I=1,MBOT 728 PADUM(I)=(PA2(I}+PADUM(I})*XX 729 CHLATOT(I)=O.O 730 DO 105 K=1,NCLASS 731 CHLADUM(K t I)=(CHLA2(K,I)+CHLADUM(K,I»*XX 732 CHLATOT(I}=CHLATOT(I)+CHLADUM(K,I) 733 105 CONTINUE 734 100 CONTINUE 735 IF(MODEL.NE.3) RETURN 736 IF(NITRO.EQ.1) THEN 737 DO 106 I=1,MBOT . 738 DO 107 K=1,NCLASS 739 107 XNCD(K,I)-(XNC2(K,I)+XNCD(K,I»*XX 740 XNHD(I)=(XNHD(I)+XNHD(I»*XX 741 106 XNOD(I)=(XN02(I)+XNOD(I»*XX 742 ENDIF 743 DO 108 1-1,MBOT 744 SID(I)=(SI2(I)+SID(I»*XX 745 DO 108 K=i,NCLASS 746 108 PCOUM(K,I)=(PC2(K,I)+PCDUM(K,I)}*XX 747 RETURN 748 END NAME TYPE I IR K MBOT XX INTEGER*2 INTEGER*2 INTEGER*2 INTEGER*2 REAL INDEX FOR LAYERS FLAG FOR NUMBER OF ITERATION INDEX FOR ALGAL CLASS NUMBER OF LAYERS RELAXATION FACTOR 749 SUBROUTINE LIGHT(RAD,TD,C2,ALBEDO) 750 C***** 751 C***** COMPUTE LIGHT DISTRIBUTION OVER A DAY IN 8 SUBSECTIONS. 752 C***** THE DISTRIBUTION IS ASSUMED TO BE SYMMETRICAL ABOUT A 753 C***** MIDDAY MAXIMUM. THE LIGHT SUBSECTIONS ARE USED IN AVELITE. 754 C***** 755 COMMON/TEMP/PARIO(4) ,PCDUM(3,40) ,XNHD(40) ,XNOD(40) , 756 + CHLADUM(3,40),XNCD(3,40),PADUM(40),SID(40) 757 C ..... ESTABLISH ARRAY OF LIGHT INTENSITY IN UNITS OF PHOTOSYNTHETICALLY 758 C ..... ACTIVE RADIATION (PAR) IN THE TOP OF THE FIRST LAyER .... 759 ALBEDO=O.087-.0000676*RAD+0.11*(1.-EXP(-.01*C2» 760 PARIO(1}=RAD*27.25*O.19/~0.125*TD)*(1.-ALBEDOI 761 PARIO(2l-RAD*27.25*O.16/ 0.125*TD)*(1.-ALB. EDO 762 PARIO(3 =RAD*27.25*O.11/ O.125*TD)*(1.-ALBEDO 763 PARIO(4 -RAD*27.25*O.04/(O.125*TD)*(1.-ALBEDO) 764 RETURN 765 END NAME TYPE 104 · . ALBEDO REAL REFLECTED RADIATION AT THE WATER SURFACE C2 REAL SUSPENDED SOLIDS CONCENTRATION AT THE WATER SURFACE RADIATION AT THE LAKE SURFACE (LANGLEY) RAD REAL TO REAL PHOTOPERIOD (HOURS) 1 1 1 1 2 1 1 1 1 1 1 2 1 1 2 1 1 1 1 1 766 SUBROUTINE MERGE(I,MBOT,LW) 767 C***** 768 C***** MERGE LAYERS THAT ARE EITHER LOW VOLUME (V < 500 M3) OR 769 C***** TOO THIN (DZ < DZLL). NEGATIVE LAYERS ARE ALSO HANDLED 770 C***** BY REDUCING THE VOLUME OF THE NEXT LOWER LAYER BY THE 771 C***** NEGATIVE VOLUME. 772 C***** 773 774 775 776 777 778 779 780 781 782 2 783 784 785 786 787 788 789 790 791 792 793 794 795 796 797 798 799 800 801 802 803 804 55 805 806 807 57 808 809 810 811 812 813 814 815 816 56 817 818 819 820 821 822 823 824 825 826 827 828 829 830 831 150 832 833 834 835 836 837 838 151 839 840 841 152 842 843 844 845 846 REAL*8 A,V,TV,ATOP COMMON/VOL/ZMAX,DZ(40),Z(40),A(40),V(40),TV(40),ATOP(41),DBL COMMON/RESULT/ T2(40),CHLATOT(40),PA2(40),PTSUM(40),B002(40), +OS02(40) t C2( 40), CD2(40) ,XN02(40) ,XNH2(40), CHLA2(3, 40) .. +PC2(3,40J,XNC2(3,40),T20(40),SI2(40) COMMON/CHOICE/MODEL,NITRO,IPRNT(6),NDAYS,NPRNT(30),NCLASS,PLOT(aO) IF(V(I).LE.O.) THEN IF(I.EQ.MBOT) THEN II=MBOT V(II-1)=V(II-1)+V(II) MBOT"MBOT-1 IF(V(lI-1).LE.O.O) THEN II"II-1 GOTO :2 ENDIF RETURN ENDIF V(I+1)-V(I)+V(I+1) DZZ=DZ(I) KK=I MBOT-MBOT-1 ELSE II"I IF(I.EQ.MBOT) 11-1-1 KK=II+ 1 VC-V( II )+V(KK) VCOMB=1./VC T2(II)=(T2(II)*V(II)+T2(KK)*V(KK»*VCOMB C2(II)=(C2(II)*V(II)+C:2(KK)*V(KK»*VCOMB CD2(II)=(CD2(II)*V(II)+CD2(KK)*V(KK»*VCOMB DO 55 K=1,NCLASS CHLA2(K,II)=(CHLA:2(K,II)*V(II)+CHLA2(K,KK)*V(KK»*VCOMB IF(MODEL.EQ.3) THEN . DO 57 K-1,NCLASS PC2(K,II)-(PC2(K,II)*V(II)+PC2(K,KK)*V(KK»*VCOMB ENDIF PA2(II)-(PA2(II)*V(II)+PA2(KK)*V(KK»*VCOMB BOD2(II)=(B002(II)*V(II)+BOD2(KK)*V(KK»*VCOMB DS02(II)-(DS02(II)*V(II)+DS02(KK)*V(KK»*VCOMB IF(NITRO.EQ.1) THEN XNH2(II)-(XNH2(II)*V(II)+XNH2(KK)*V(KK»*VCOMB XN02(II)-(XN02(II)*V(II)+XN02(KK)*V(KK»*VCOMB DO 56 K=1,NCLASS XNC2(K,II)=(XNC:2(K,II)*V(II)+XNC2(K,KK)*V(KK»*VCOMB ENDIF V( II )=VC DZ(II)=DZ(II)+DZ(KK) Z(II)=Z(II)+DZ(KK)*0.5 DZZ=O.O MBOT=MBOT-1 IF(LW.GT.I) LW=LW-1 IF(I.EQ.MBOT) GO TO 3 END IF DO 100 K-KK,MBOT T2(K)=T2(K+1) C2(K)"C2(K+1) CD2(K)-CD2(K+1) DO 150 KI=1,NCLASS CHLA2(KI,K)=CHLA2(KI,K+1) PA2(K)-PA2(K+1) BDD2(K)=BOD2(K+1) DS02(K)=DS02(K+1) IF(MDDEL.EQ.3) THEN SI2(K)"SI2(K+1) DO 151 KI=1,3 PC2(KI,K)=PC2(KI,K+1) IF(NITRO.EQ.1) THEN DO 152 KI"1,3 XNC2(KI,K)-XNC2(KI,K+1) XNH2(K)-XNH2(K+1) XND2(K)-XN02(K+1) ENDIF ENDIF V(K)-V(K+1) 105 847 DZ(K)=DZ(K+1) 848 Z(K)=Z(K+1)-DZZ 849 100 CONTINUE 850 3 ZMAX= Z(MBOT) + 0.5*DZ(MBOT) 851 RETURN 852 END NAME TVPE DZZ I II K KI KK LW REAL INTEGER*2 INTEGER*2 INTEGER*2 INTEGER*2 INTEGER*2 INTEGER*2 DUMMY VARIABLE FOR LAYER THICKNESS INDEX FOR LAYER INDEX FOR LAYER TO BE MERGED INDEX FOR LAYERS BELOW THE MERGED LAYERS INDEX FOR ALGAL CLASS MBOT INTEGER*2 INDEX OF LAYER RECEIVING THE LAYER TO BE MERGED INDEX OF LAYER AT THE MIXED LAYER DEPTH NUMBER OF LAYERS VC REAL VCOMB REAL COMBINED VOLUME OF THE TWO LAYERS BEING MERGED RECIPROCAL OF THE COMBINED VOLUMES 1 1 1 1· 1 853 SUBROUTINE HEBUG(IL,TS,QN,HS,HA,HBR,HE,HC, 854 +TAIR,TDEW.CR,RAD.WIND,IPAN,PCOEFF) 855 C***** 856 C***** COMPUTE THE TEMPERATURE PROFILE USING ROUTINES FLXOUT AND 857 C***** FLXIN FOR THE SURFACE HEAT EXCHANGE. SOLUTION IS BY THE 858 C***** IMPLICIT CENTRAL DIFFERENCE FORMULATION. CONMIX CALLED TO 859 C***** CHECK FOR AND RESOLVE DENSITY INSTABLITIES BETWEEN LAYERS. 860 C***** 861 REAL*8 A,V.TV,ATOP,AK BK,CK,DK 862 COMMON!VOL/ZMAX'DZ(40)'Z(40)'A(40)'V(40)~TV(40)'ATOP(41),DBL 863 COMMON!RESULT! T2(40),CHLATOT(40),PA2(40 ,PTSUM(40).BOD2(40). 864 +DS02(40),C2(40),CD2(40),XN02(40),XNH2(40 ,CHLA2(3,40), 865 +PC2(3,40),XNC2(3,40),T20(40),SI2(40) 866 COMMON/TEMP!PARIO(4),PCDUM(3,40),XNHD(40),XNOD(40), 867 + CHLADUM(3.40),XNCD(3,40),PADUM(40),SID(40) 868 COMMDN!STEPS/DZLL,DZUL,MBOT,NM,NPRINT,MDAY,MONTH,ILAY,DY 869 COMMON!FLOW!HMK(41),QE(40),FVCHLA(5),PE(5,41) 870 COMMON!WATER!BETA,EMISS,XK1,XK2,HKMAX,WCOEF,WSTR,RHOW,HCAP 871 COMMON!SOLV! AK(40),BK(40),CK(40),DK(40) 872 DIMENSION Q(40) 873 C ... CALCULATION OF THE HEAT ABSORPTION FROM METEOROLOGICAL PARAMETERS 874 C ..• IN A COLUMN OF WATER 875 C ... CALCULATION OF HEAT FLUXES INTO THE WATER BODY 876 CALL FLXIN(HS,HA,TAIR,RAD,CR,C2) 877 CALL FLXOUT(TS,HBR,HE,HC,TAIR,TOEW,WINO,IPAN,PCOEFF,WCOEF) 878 HQOUT=HBR+HE+HC 879 C ... CALCULATION OF EXTINCTION COEFF. (ETA) AS A FUNCTION OF SUSPENDED 880 C ... SEDIMENT CONCENTRATION 881 EiA-XK1+0.04317*C2(1)+XK2*CHLATOT(1) 882 C ... CALCULATION OF HEAT ABSORBED IN EACH LAYER 883 HQ=(1.-BETA)*HS 884 EX=EXP(-ETA*DZ(1» 885 Q(1)=«BETA*HS+HA-HQOUT)*ATOP(1)+HQ*(ATOP(1)-EX*ATOP(2») 886 + /(1OOO.*V(1» 887 C ... CONVERSION FACTOR OF 1000 USED FOR DENSITY*HEAT CAPACITY OF WATER 888 HQ=HQ*EX 889 C ... CALCULATE THE SOURCE TERM Q FOR EACH LAYER 890 DO 10 I=2,MBOT 891 ETA=XK1+0.04317*C2(I)+XK2*CHLATOT(I) 892 EX-EXP(-ETA*OZ(I» 893 Q(I)=HQ*(ATOP(I)-ATOP(I+1)*EX)!(1000.*V(I» 894 HQ=HQ*EX 895 10 CONTINUE 896 CALL SETAMK(IL,MBOT,WIND,HKMAX) 897 C ... SET-UP COEFFICIENTS FOR TRI-OIAGONAL MATRIX 898 DO 100 I=2,MBOT-1 899 01= 2./(A(I)*DZ(I» 900 D2=D1*ATOP(I)*HMK(I)/(DZ(I)+DZ(I-1» 90103=.D1*ATOP(I+1)*HMK(I+1)/(OZ(I)+OZ(I+1» 902 AK!I!--02 903 BK I =1.0+02+03 904 CK I =-03 905 OK I -T2(I)+Q(I) 906 100 CONTINUE 907 OK(1)-T2(1)+Q(1) 908 DK(MBOT)-T2(MBOT)+Q(MBOT) 909 AK~1).O. 910 CK 1)=-2.*ATOP(2)*HMK(2)!(A(1)*DZ(1)*(DZ(1)+DZ(2») 911 BK 1)-1.-CK(1) 912 I-MBOT 914 BK(I =1.-AK(I) 913 AK(I}--2.*ATOP(I)*HMK(I)!(A(I)*OZ(I)*(DZ(I)+OZ(I-1») 915 CK(I =0.0 916 CALL SOLVE(T2.MBOT) 917 TS -T2(1) 106 918 919 CALL CONMIX(IL,TS,MBOT) DO 90 1=1,IL 920 90 921 T2(I)=TS ON=HS+HA-HQOUT 922 RETURN 923 END NAME TYPE CR REAL D1 REAL D2 REAL D3 REAL ETA REAL EX REAL HA REAL HBR REAL HC REAL HE REAL CLOUDINESS RATIO FROM METEOROLOGICAL DATA FILE FINITE DIFFERENCE COEFFICIENT FINITE DIFFERENCE COEFFICIENT FINITE DIFFERENCE COEFFICIENT SOLAR RADIATION EXTINCTION COEFFICIENT (1/M) FRACTION OF INCOMING SOLAR RADIATION AT A GIVEN LAYER ATMOSPHERIC RADIATION (KCAL/M2/DAY) BACK RADIATION (KCAL/M2/DAY) CONVECTIVE HEAT LOSS (KCAL/M2/DAY) EVAPORATIVE HEAT LOSS (KCAL/M2/DAY) HQ HQOUT HS REAL REAL REAL INTEGER*2 INTEGER*2 REAL SOLAR RADIATION AT A GIVEN LAYER SURFACE HEAT LOSS AS THE SUM OF HC, HE, AND HBR INCOMING SOLAR RADIATION (KCAL/M2/DAY) IL IPAN PCOEFF Q REAL INDEX OF LAYER AT THE MIXED LAYER DEPTH FLAG FOR PAN EVAPORATION COMPUTATIONS IN FLXOUT PAN EVAPORATION COEFFICIENT ENERGY ADDEO TO A LAYER FROM ATTENUATION OF SOLAR QL ON RAD TAIR TDEW TS WIND NAME C CC HAN HSN HSR RAD TC TCA REAL REAL REAL REAL REAL REAL REAL RADIATION (KCAL/M2/DAY) TOTAL SURFACE HEAT LOSS (KCAL/DAY) NET HEAT FLUX (KCAL/M2) INCOMING SOLAR RADIATION (LANGLEY) AVERAGE AIR TEMPERTURE (DEGREES C) DEW POINT TEMPERATURE (DEGREES C) MIXED LAYER TEMPERATURE AFTER HEAT INPUT (DEGREES C) WINO SPEED (MPH) 924 SUBROUTINE FLXIN(HSN,HAN,TC,RAD,CC,C) 925 C ... CALCULATION OF THE TOTAL RADIATION FLUX INTO A BODY OF WATER IN 926 C ... FROM NET SOLAR RADIATION (HSR) AND NET ATMOSPHERIC RADIATION (HSN) 927 C ... IN KCAL/M*M 928 C ... IDSO JACKSON FORMULA USED FOR ATM. RADIATION 929 C ... CONVERT AIR TEMPERATURE IN DEGREE C TO DEGREE ABSOLUTE 930 TCA=TC+273. 931 TCA=TCA*TCA 932 HAN=1.171E-6*(1.-0.261*EXP(-6.77E-4*TC*TC» 933 +*(TCA*TCA)*(1.+0.17*CC*CC) 934 C ... CALCULATION OF NET SOLAR RADIATION AND CONVERSION TO KCAL/M*M/DAY 935 C ... CALCULATION OF REFLECTED SOLAR RADIATION HSR 936 C ... HSRW---- DUE TO CLEAR WATER USING KOBERGS CURVES 937 C ... HSRS---- DUE TO SUSPENDED SEDIMENTS AT THE WATER SURFACE 938 HSR=(0.087-6.76E-5*RAD+0.11*(1.-EXP(-0.01*C»)*RAD 939 HSN=(RAD-HSR)*10. 940 RETURN 941 END TYPE REAL REAL REAL REAL SURFACE LAYER SUSPENDED SOLIDS CONCENTRATION (MG/L) SURFACE LAYER DISSOLVED SOLIDS CONCENTRATION (MG/L) SOLAR RADIATION AT THE WATER SURFACES (KCAL/M2/DAY) REFLECTED SOLAR RADIATION (KCAL/M2/DAY) REAL REAL SOLAR RADIATION ENTERING THE SURFACE LAYER (LANGLEY) INCOMING SOLAR RADIATION (LANGLEY) REAL REAL AVERAGE AIR TEMPERATURE (DEGREES C) AVERAGE AIR TEMPERATURE (DEGREES K) SQUARED 942 SUBROUTINE FLXOUT (TT,HBR,HE,HC,TAIR,TD,WIND,IPAN,PCOEFF,WCOEF) 943 C ... CALCULATES THE ENERGY FLUX OUT OF A BODY OF WATER FROM 944 C ... EVAPORATIVE HEAT LOSS (HE), CONVECTIVE HEAT LOSS (HC), AND 945 C ... BACK RADIATION (HBR) IN KCAL/M*M/DAY. 946 C ... CONVERSION OF TEMP. VALUES FROM OEG. C TO OEG. ABSOLUTE 947 TSK=TT+273.15 948 TAK=TAIR+273.15 949 FCN=WCOEF *WIND 950 IF(IPAN.GT.O) GO TO 20 951 C ... EVALUATES SATURATION VAPOR PRESSURE VALUES IN MB USING MAGNUS TETONS 952 C ... FORMULA 953 ESA=6.1078*EXP(17.2693882*TT/(TSK-35.86» 954 C ... EVALUATES SATURATION VAPOR PRESSURE VALUE 955 EA c 6.0353*10.**(7.45*TO/(235.+TD» 956 C ... CALCULATES EVAPORATIVE HEAT LOSS 957 HE=FCN*(ESA-EA) 958 GO TO 30 107 959 20 HE=2.54*TD*PCOEFF*590.*10. 960 C ... CALCULATES BACK RADIATION 961 30 HBR~(1.171E-6*0.97*TSK*TSK*TSK*TSK) 962 C ... CALCULATES CONDUCTIVE LOSS USING BOWENS RATIO 963 HC=FCN*0.61793*(TSK-TAK) 964 RETURN 965 END NAME TVPE EA ESA FCN HBR HC REAL REAL REAL REAL REAL VAPOR PRESSURE (MB) SATURATED VAPOR PRESSURE (MB) WIND COEFFICIENT ON CONVECTIVE AND EVAPORATIVE HEAT TRANSFER BACK RADIATION (KCAL/M2/DAV) CONVECTIVE HEAT FLUX (KCAL/M2/DAV) HE IPAN PCOEFF TAIR TAK REAL INTEGER*2 REAL REAL REAL EVAPORATIVE HEAT FLUX (KCAL/M2/0AV) FLAG FOR PAN EVAPORATION COMPUTATIONS PAN COEFFICIENT AVERAGE AIR TEMPERATURE (DEGREES KC) AVERAGE AIR TEMPERATURE (DEGREES ) PHOTOPERIOD (HOURS) TO REAL TSK REAL TT WCOEF WIND REAL REAL REAL SURFACE TEMPERATURE (DEGREES K) SURFACE TEMPERATURE (DEGREES C) WIND FUNCTION COEFFICIENT WIND SPEED (MPH) 966 SUBROUTINE SETAMK(ILAV,MBOT,W,HKMAX) 967 C***** 968 C***** COMPUTE VERTICAL DIFFUSION COEFFICIENT IN EACH LAVER. 969 C***** DIFFUSION COEFFICIENT BETWEEN LAVERS AS THE HARMONIC 970 C***** MEAN OF THE DIFFUSION COEFFICIENTS IN AD~ACENT LAVERS. 971 C***** 972 REAL*8 A,V,TV,ATOP 973 DIMENSION AMK(41) 974 COMMON/FLOW/HMK(41),QE(40),FVCHLA(5),PE(5,41) 976 COMMON/VOL/ZMAX'DZ(40)'Z(40)'A(40)'V(40)~TV(40)'ATOP(41),DBL 976 COMMON/RESULT/ T2(40),CHLATOT(40),PA2(40 ,PTSUM(40),B002(40), 977 +DS02(40)tC2(40),CD2(40),XN02(40),XNH2(40 ,CHLA2(3,40), 978 +PC2(3,40J,XNC2(3,40),T20(40),SI2(40) 979 C ... VERTICAl DIFFUSION COEFFICIENT IN THE MIXED LAVER 980 C ... COMPUTED AS A FUNCTION OF THE WIND SPEED 981 HKC=HKMAX*8.66E-3 982 DKM=28.*W**1.3 983 IF(DKM.LT.HKMAX) DKM=HKMAX 984 DO 100 I=1,ILAV 986 100 AMK(I)=DKM 986 C ... HMK TEMPORARILY USED TO STORE DENSITIES. 987 C ... DIFFUSION COEFFICIENT BELOW THE MIXED LAYER COMPUTED AS 988 C .•• A FUNCTION OF THE SQUARE OF THE BRUNT-VASALA FREQUENCV (SQN) 989 C .•. SQN= (G/RHO) * D(RHO)/DZ 990 DO 200 I=ILAY,MBOT 991 200 HMK(I)=RHO(T2(I),C2(I),CD2(I» 992 DO 300 I=ILAY+1,MBOT-1 993 AVRHO=(HMK(I-1)+HMK(I+1»*0.6 994 SQN=ABS(HMK(I-1)-HMK(I+1»/«Z(I+1)-Z(I-1»*AVRHO)*9.81 996 IF(SQN.LT. 0.000075) SQN=0.000076 996 300 AMK(I)=HKC/SQRT(SQN) 997 C--------ASSUME AMK(MBOT-1) = MAX AMK FOR SQN-0.OOOO76 998 AMK(MBOT)=AMK(MBOT-1) 999 DO 400 I=2,MBOT 1000 400 HMK(I)=2.*AMK(I)*AMK(I-1)/(AMK(I)+AMK(I-1» 1001 HMK(1)aO.O 1002 RETURN 1003 END NAME TVPE AMK REAL AVRHO REAL DKM REAL HKC REAL HKMAX ILAY MBOT SQN W REAL INTEGER*2 INTEGER*2 REAL REAL DIFFUSION COEFFICIENT IN A LAYER (M2/DAV) AVERAGE DENSITY IN AD~ACENT LAYERS DIFFUSION COEFFICIENT IN THE MIXED LAYER (M2/0AY) DUMMY VARIABLE FOR MAXIMUM DIFFUSION COEFFICIENT IN THE HYPOLIMNION MAXIMUM DIFFUSION COEFFICIENT IN THE MIXED LAYER (M2/DAV) INDEX OF LAVER AT THE MIXED LAYER DEPTH NUMBER OF LAYERS SQUARE OF THE BRUNT-VASALA FREQUENCY WIND SPEED (MPH) PASS ONE NO ERRORS DETECTED 1003 SOURCE LINES 108 1 $DEBUG 2 $STORAGE:2 3 $NOFLOATCALLS 4 SUBROUTINE MTHDATA(MONTH.KDAYS.MYEAR) 5 C ............ ... 6 C ........... READ MONTHLY METEOROLOGICAL DATA 7 C ........ . 8 COMMON/FILE/ DIN.MET.FLO.TAPE8.TAPE1.IREC 9 COMMON/MTHD/TAIR(31).TDEW(31).RAD(31).CR(31).WIND(31). 10 +PR(31).DRCT(31) 11 CHARACTER.16 DIN.MET.FLO.TAPE8,TAPE1 12 C ... FIND MET DATA FOR FIRST MONTH OF SIMULATION 13 IF(KDAYS.EQ.O) THEN 14 MTH=MONTH 15 MYR~MYEAR 16 READ(9,"')MONTH.KDAYS.MYEAR 17 30 IF(MONTH.EQ.MTH.AND.MYR.EQ.MYEAR) GO TO 20 18 READ(9.*.END~10.ERR~30)MONTH,KDAYS.MYEAR 19 GO TO 30 20 10 WRITE(*.1oo1) 21 WRITE(8.1001) 22 STOP 23 ELSE 24 READ(9.*) MONTH.KDAYS.MVEAR 25 ENDIF 26 C ... READ IN MONTH. NO. OF DAYS IN THE MONTH AND YEAR 27 C ... READ IN AIR TEMP. AND DEW PT. TEMP. IN CELSIUS. WIND VELOCITY IN 28 C ... M.P.H .• SOLAR RADIATION IN LANGE LV/DAY AND CLOUD COVER IN TENTHS 29 20 READ(9 .... ) (TAIR(K).K~1.KDAYS) 30 READ(9 .... ) (TDEW(K).KK1.KDAYS) 31 READ(9.*) (PR(K).K=1.KDAYS) 32 READ(9.*) (WIND(K).K~1.KDAYS) 33 READ(9.*) (DRCT(K).K=1.KDAVS) 34 READ(9 .... ) (CR(K).K=1.KDAYS) 35 READ(9.*) (RAD(K).K=1.KDAYS) 36 DO 100 K~1.KDAYS 37 TAIR(K)~(TAIR(K)-32,)"'0.5556 38 TDEW(K)=(TDEW(K)-32.)*0.5556 39 100 CR(K)=(100.-CR(K» •. 01 40 1001 FORMAT(//.5X.70('.').//.20X.'PROGRAM ABORTED.'./.10X. 41 + 'METEOROLOGICAL DATA FILE DOES NOT'./.15X, 42 + 'MATCH YEAR OF SIMULATION'.//.5X.70('''''» 43 RETURN 44 END NAME TYPE K KDAYS MONTH MTH MYEAR MYR 45 46 47 48 49 50 51 52 53 54 55 56 57 INTEGER"'2 INTEGER"'2 INTEGER"'2 INTEGER"'2 INTEGER"'2 INTEGER.2 INDEX OF DAY OF MONTH NUMBER OF DAYS IN THE MONTH MONTH FIRST MONTH OF SIMULATION YEAR YEAR SUBROUTINE PDEPTH( RHOMIX.RHOIN.Q.S.IHP.QENIN. +M80T.SUMZ.FT.WIDTH,HP) C .......... ... C ......... COMPUTE DEPTH AT PLUNGE POINT AND INITIAL ENTRAINMENT AT C ........... THE PLUNGE POINT (GAMAIN). (BY AKYAMA) C· ......... ... REAL"'8 A.V,TV,ATOP COMMON/VOL/ZMAX,DZ(40),Z(40),A(40),V(40).TV(40),ATOP(41),DBl C ... COMPUTE DEPTH AT PLUNGE POINT EPSIN=(RHOIN-RHOMIX)/RHOMIX IF(EPSIN.GT.O.) GO TO 7 IHP=1 HP=O. 58 r RETURN X~Q/(86400"'WIDTH) IF(S.LE.6.66667E-3) THEN 59 7 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 200 75 HP·1.1 ... (FT ... X ... X/(S ... EPSIN.9.81») ...... (.3333) GAMAIN"0.15 ELSE HP=1.6*(X.X/(EPSIN.9.81»"'*(.3333) GAMAIN= 1. 8 ENDIF QENIN=GAMAIN.Q SUMZ~O. IHP=O DO 200 I~1,MBOT IHP"'IHP+1 SUMZ~SUMZ+DZ(I) IF(SUMZ.GE.HP) GO TO 5 CONTINUE HP=SUMZ 109 76 5 77 RETURN END 78 SUBROUTINE PMETE(HS,RAD,HA,WIND,HBR,P,HE,TAIR,HC,TDEW,HED,HEV, 79 +QNET,DMIX,ZEUPH,SECCHI) 80 c***** 81 C***** OUTPUT TABLE OF METEOROLOGICAL AND HEAT FLUX VALUES 82 C***** 83 COMMON/FILE/DIN,MET,FLO,TAPE8,TAPE1,IREC 84 CHARACTER*16DIN,MET,FLO,TAPE8,TAPE1 85 WRITE(a,2000) HS,RAD,HA,WIND,HBR,P,HE,TAIR,HC,TDEW,HED,HEV, 86 +QNET,DMIX,ZEUPH,SECCHI 87 2000 FORMAT(/,5X,26HMETEOROLOGICAL INFORMATION,/,5X.26(1H-). 88 +/,7X.17HNET SOLAR RAD. = .4X.F9.2.13H KCAL/M*M/DAY. 89 +25X, 13HSOLAR RAD. = • 9X.F6.1.8H LANGLEY. 90 +/.7X.16HNET ATM. RAD. = ,5X,F9.2.13H KCAL/M*M/DAY. 91 +25X,16HWIND VELOCITY = , 7X,F5.1.7H M.P.H .. 92 +/.7X, 12HBACK RAD. = • 9X.F9.2.13H KCAL/M*M/DAY. 93 +25X,16HPRECIPITATION = , 8X.F6.3.13H METER(S)/DAY. 94 +/,7X.19HEVAPORATIVE LOSS = .2X,F9.2.13H KCAL/M*M/DAY, 95 +25X.12HAIR TEMP. = .12X.F5.2.10H DEGREES C. 96 +/,7X.20HCONDENSATION LOSS = , 1X.F9.2.13H KCAL/M*M/DAY. 97 +25X.16HDEW PT. TEMP. = • 8X.F5.2.10H DEGREES C. 98 +/,7X.20HEV. HEAT TRANSFER = .5X,F7.4.10H M/DAY OR ,E10.4.7H M**3/D 99 +2HAY./.7X, 19HNET HEAT FLUX IN = • 2X.F9.2.9H KCAL/M*M. 100 +/,7X.' MIXED LAYER DEPTH (M)=',2X.F5.2. 101 +38X,'EUPHOTIC DEPTH (M) ='.2X.F5.2. 102 +/8X,'SECCHI DEPTH (M) c'.7X.F5.2) 103 RETURN 104 END 105 SUBROUTINE AVGLITE(IEUPH) 106 C***** 107 C***** COMPUTE AVERAGE LIGHT IN A LAYER FOR FOUR DIFFERENT 108 C***** TIME SEGMENTS IN A DAY. (FOLLOWS WORK BY J. CARDONI) 109 C***** 110 REAL*8 A.V.TV,ATOP.SVOL 111 COMMON/TEMP/PARIO(4).PCDUM(3.40).XNHD(40),XNOD(40). 112 + CHLADUM(3,40).XNCD(3.40),PADUM(40).SID(40) 113 COMMON/RESULT/ T2(40),CHLATOT(40).PA2(40).PTSUM(40),BOD2(40), 114 +DS02(40).C2(40).CD2(40).XN02(40).XNH2(40),CHLA2(3.40). 115 +PC2(3.40).XNC2(3,40).T20(40),SI2(40) 116 COMMON/SOURCE/RM(3,40).PROD(40),XMR(3,40).PRODSUM(40) 117 COMMON/SUB/SDZ(60).SZ(60),LAY(40).AVGI(4,60).SVOL(60) 118 COMMON/VOL/ZMAX.DZ(40),Z(40).A(40).V(40),TV(40),ATOP(41).DBL 119 COMMON/WATER/BETA.EMISS,XK1,XK2,HKMAX.WCOEF,WSTR 120 DIMENSION XIZ1(4) 121 UU~O 122 DO 50 K=1.4 123 50 XIZ1(K)=PARIO(K) 124 DO 100 I=1.IEUPH 1 125 EK=XK1+.043*C2(I)+XK2*CHLATOT(I) 1 126 DO 100 J=1.LAY(I) 2 127 JJ=JJ+1 2 128 DO 100 K=1.4 3 129 X=EK*SDZ(JJ) 3 130 X1=EXP(-X) 3 131 AVGI(K.JJ)=XIZ1(K)*(1.-X1)/X 3 132 100 XIZ1(K)=XIZ1(K)*X1 NAME EK I 133 RETURN 134 END TYPE REAL INTEGER*2 LIGHT EXTINCTION COEFFICIENT (1/M) INDEX OF LAYERS IEUPH INTEGER*2 INDEX OF LAYER AT THE EUPHOTIC DEPTH INDEX OF SUBLAYERS FORM SUBLAY U UU K X X1 XIZ1 INTEGER*2 INTEGER*2 INTEGER*2 REAL REAL REAL INDEX OF SUBLAYERS WITHIN A LAYER INDEX OF TIME SEGMENTS IN THE PHOTOPERIOD LIGHT ATTENUATION COEFFICIENT THROUGH A SUBLAYER FRACTION OF LIGHT ATTENUATED THROUGH A SUBLAYER LIGHT ENTERING THE TOP OF THE SUBLAYER (KCAL/M2/DAY) 135 SUBROUTINE PRODAVG(IEUPH.MBOT.KK) 136 C***** 137 C***** COMPUTE LIGHT LIMITATION COEFFICIENT ON ALGAL GROWTH 138 C***** IN EACH SUBLAYER FOR FOUR LIGHT PERIODS DETERMINED IN 139 C***** AVELITE. (FOLLOWS WORK BY U. CARDONI WITH 140 C***** LIMITATION/INIHIBITION FUNCTION FROM MEGARD.) 141 C***** 110 1 1 2 2 2 3 3 3 2 1 NAME I 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 REAL*8 A,V,TV,ATOP,SVOL COMMON/SUB/ SDZ(60),SZ(60),LAY(40),AVGI(4,60),SVOL(60) COMMON/SOURCE/ RM(3,40),PROD(40),XMR(3,40),PRODSUM(40) COMMON/PHYTO/PDEL(3),PMAX(3),PMIN(3),THR(3),THM(3),XKR1(3), +XKR2(3),XKM(3),HSCPA(3),HSC1(3),HSC2(3),UPMAX(3),THUP(3), +GROMAX(3),TMAX(3),TOPT(3),XNMAX(3),XNMIN(3),UNMAX(3), THUN(3), +HSCN(3),HSCNH(3),XNDEL(3),IDIATOM,CHLMEAN,CHLMAX,SECCHI COMMON/VOL/ZMAX,DZ(40),Z(40),A(40),V(40),TV(40),ATOP(41),DBL DIMENSION PSUB(60) 1.11.1-0 X-1.+2.*SQRT(HSC1(KK)/HSC2(KK» DO 200 1-1,IEUPH SUMPV-O. 00 300 1.1-1,LAY(I) 1.11.1"1.11.1+1 PSUB(I.II.I)"O DO 400 K"1,4 P-AVGI(K,I.II.I)*X/(AVGI(K,I.II.I)+HSC1(KK)+AVGI(K,I.II.I) + *AVGI(K,I.II.I)/HSC2(KK»*0.25 400 PSUB(I.II.I)-PSUB(I.II.I)+P 300 SUMP V .. SUMPV+PSUB(I.II.I)*SVOL(I.II.I) 200 PROD(I)- SUMPV/V(I) IF(IEUPH.LT.MBOT) THEN DO 500 I-IEUPH+1,MBOT 500 PRODO )"0. TYPE ENDIF RETURN END INTEGER*2 INDEX OF LAYERS IEUPH INTEGER*2 INDEX OF LAYER AT THE EUPHOTIC DEPTH 1.1 INTEGER*2 INDEX OF SUBLAYERS 1.11.1 INTEGER*2 INDEX OF SUB LAYERS WITHIN A LAYER K INTEGER*2 INDEX OF SEGMENTS OF THE PHOTOPERIOD KK INTEGER*2 INDEX OF ALGAL CLASS MBOT INTEGER*2 NUMBER OF LAYERS P REAL LIGHT LINITATION COEFFICIENT IN ONE PHOTOPERIOD SEGMENT PSUB REAL SUM OF LIGHT LIMITATION COEFFICIENTS IN A SUBLAYER SUMPV REAL SUM OF LIGHT LIMITATION COEFFICIENTS IN A LAVER X 1 1 1 2 2 1 1 1 1 2 2 1 REAL DUMMY VARIABLE USED IN LIGHT LIMITATION RATIO 170 SUBROUTINE PTOTALS(MAXMTH,MXDAV,CHLMAX,CHLMEAN) 171 c***** 172 C***** COMPUTE TOTAL PHOSPHORUS AS A SUM OF ORTHO, INTERNAL, AND 173 C***** DETRITAL PHOSPHORUS. TOTAL CHLOROPHYLL-A AS THE SUM OF 174 C***** CHLOROPHYLL-A IN ALL ALGAL CLASSES. 175 c***** 176 REAL*8 A,V,TV,ATOP 177 COMMON/RESULT/ T2(40),CHLATOT(40),PA2(40),PTSUM(40),BOD2(40), 178 +DS02(40)tC2(40),CD2(40) ,XN02(40),XNH2(40) ,CHLA2(3.40) , 179 +PC2(3,40},XNC2(3,40),T20(40),SI2(40) 180 . COMMON/YIELD/VCA,YCH02,V2CH02,VCBOD,YPBOD,YZW,YPZP,YNZP,YZDO, 181 +YSCHL,YNHBOD,BRNO,BRNH,XKNNH,THNNH,YPCHLA,BODK20,SB20,BRR 182 COMMON/CHOICE/MODEL,NITRO,IPRNT(S),NDAYS,NPRNT(30),NCLASS,PLOT(30) 183 COMMON/VDL/ZMAX, DZ(40), Z(40), A(40), V(40), TV(40), ATDP(41), DBL 184 COMMDN/STEPS/DZLL,DZUL,MBOT,NM,NPRINT,MDAY,MONTH,ILAY,DY 185 IF (MODEL.EQ.3) THEN 186 DO 120 I=1,MBOT 187 CHLATOT(I)-O.O 188 PTSUM(I )-0.0 . 189 DO 350 K3=1,NCLASS 190 CHLATOT(I)-CHLATOT(I)+CHLA2(K3,I) 191 350 PTSUM(I)-PTSUM(I)+PC2(K3,I) 192 120 PTSUM(I)=PA2(I)+PTSUM(I)+BOD2(I)*YPBOD 193 ELSE 194 DO 121 1-1,MBOT 195 CHLATOT(I)"O.O 196 PTSUM(I)=O.O 197 DO 351 K3-1,NCLASS 198 CHLATOT(I)=CHLATDT(I)+CHLA2(K3,I) 199 351 PTSUM(I)=PTSUM(I)+CHLA2(K3,I)*YPCHLA 200 121 PTSUM(I)-PA2(I)+PTSUM(I)+BOD2(I)*VPBOD 201 ENDIF 202 CHL-O.O 203 DO 130 1-1,MBOT 204 CHL-CHL+CHLATOT(I)*DZ(I) 205 IF(Z(I).GT.2.0) GOTO 131 206 130 CONTINUE 207 131 CHLMEAN-CHL/(Z(I)+.5*DZ(I» 208 IF(CHLMEAN.GT.CHLMAX) THEN 209 CHLMAX-CHLMEAN 210 MXDAY-MDAY 211 MAXMTH-MONTH 111 212 ENDIF 213 RETURN 214 END NAME TYPE CHL CHLMAX CHLMEA I REAL DUMMY VARIABLE FOR DEPTH WEIGHTED CHlORO-A IN TOP 2 METERS MAXIMUM CHlORO-A CONCENTRATION IN THE SIMULATION REAL TWO-METER MEAN CHlOROPHYll-A CONCENTRATION (MG/L) INDEX OF lAYERS K3 MAXMTH MXDAY REAL INTEGER*2 INTEGER*2 INTEGER*2 INTEGER*2 INDEX OF ALGAL CLASS MONTH OF MAXIMUM CHlORO-A CONCENTRATION IN SIMULATION DAY OF MAXIMUM CHlORO-A CONCENTRATION IN SIMULATION 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 NAME DUM1 DUM2 DUM3 DUM4 DUM5 I MBOT PV X Xi SUBROUTINE RESSETL (MBOT) C ••• C ... THIS SUBROUTINE COMPUTES CONC. OF SUSPENDED SEDIMENTS IN THE lAKE C ... ADVECTIVE, DIFFUSIVE EQUATION IS SOLVED USING FUllY IMPLICIT HYBRID C ... FINITE DIFFERENCE SCHEME C ••• REAl*8 A,V,TV,ATOP,AK,BK,CK,DK DIMENSION PV(42) COMMON/SO lV/ AK(40) ,BK(40),CK(40) ,DK(40) COMMON/FlOW/HMK(41),QE(40),FVCHlA(5),PE(5,41) COMMON/VOl/ZMAX,DZ(40),Z(40),A(40),V(40),TV(40),ATOP(41),DBl COMMON/RESUlT/ T2(40),CHlATOT(40)tPA2(40),PTSUM(40),BOD2(40), +DS02(40)tVAR2(40),CD2(40),XN02(40J,XNH2(40),CHlA2(3,40), +PC2(3,40J,XNC2(3,40),T20(40),SI2(40) C ... OBTAIN FAll VELOCITY OF SEDIMENT FOR EACH lAyER ..... DO 10 I=1,MBOT 10 PV(I+1)=FV(T2(I» C ... CAlCUlATION OF GRID PEClET NUMBER FOR EACH lAyER ..... DO 20 I=2,MBOT 20 PE(5,I)=«PV(I)+PV(I+1»*.5)*«OZ(I-1)+DZ(I»*.5)/HMK(I) PE(5,1)=10.0 PE(5 t MBOT+1)=10. PV(1J=-PV(2) PV(MBOT+2)=-PV(MBOT+1) DO 30 I=1,MBOT DUM1=1./(A(I)*DZ(I)*2.0) DUM2=ATOP(I)*(PV(I+1)+PV(I»*DUM1 DUM3=ATOP(I+1)*(PV(I+1)+PV(I+2»*DUM1 C------- DUM4=MAX(0.O, «1.0-0.1*ABS(PE(5,I»**5)/PE(5,I) X=1.0-0.1*ABS(PE(5,I» X1=X*X DUM4=X1*X1*X/PE(5,I) X=1.0-0.1*ABS(PE(5,I+1» X1=X*X DUM5~X1*X1*X/PE(5,I+1) IF(DUM5.lT.0.0) DUM5=0.0 IF(DUM4.li.0.0) DUM4=0.0 AK(I)=-DUM2*(1.0+DUM4) BK(I)=1.0+DUM3*(1.0+DUM5)+DUM2*DUM4-PV(I)*2.0*DUM1*(ATOP(I+1) +-ATOP{I » CK(I)=-DUM3*DUM5 OK (I )=VAR2(I) 30 CONTINUE CAll SOlVE(VAR2,MBOT) DO 80 I=1,MBOT IF(VAR2(I).lT.0.) VAR2(I)-0. 80 CONTINUE RETURN END TYPE REAL DUMMY VARIABLE USED FOR REAL DUMMY VARIABLE USED FOR REAL DUMMY VARIABLE USED FOR REAL DUMMY VARIABLE USED FOR REAL DUMMY VARIABLE USED FOR INTEGER*2 INDEX OF lAYERS INTEGER*2 NUMBER OF lAYERS REAL FAll VELOCITY FINITE DIFFERENCE FORMULATION FINITE DIFFERENCE FORMULATION FINITE DIFFERENCE FORMULATION FINITE DIFFERENCE FORMULATION FIN.ITE DIFFERENCE FORMULATION REAL DUMMY VARIABLE USED IN POWER lAW COMPUTATIONS REAL DUMMY VARIABLE USED IN POWER lAW COMPUTATIONS 264 FUNCTION RHO (T,C,CD) 265 C ... CAlCUlATES THE DENSITY OF WATER AS A FUNCTION OF TEMPERATURE PLUS 266 C ... DENSITY DUE TO TOTAL SOLIDS (SUSPENDED AND DISSOLVED) 267 RHO=( .999878+T*(6. i6608E-5+T*(-S. 14577E-6+T*4. 76102E-8»)*1000. 268 ++(C+CD)*O.OO1 112 .' 269 RETURN 270 END NAME TVPE C CD T NAME AZ I II VOUM NAME AZ I MBOT 271 272 273 274 275 276 277 27B 279 280 281 282 283 2B4 285 2B6 2B7 2BB 2B9 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 REAL REAL REAL SUSPENDED SOLIDS CONCENTRATION IN A LAVER (MG/L) DISSOLVED SOLIDS CONCENTRATION IN A LAVER (MG/L) TEMPERATURE IN A LAVER (DEGREES C) SUBROUTINE SETUP C***** C***** DETERMINE THE INITIAL THICKNESS, VOLUME, AND AREA OF LAVERS C***** AND THE TOTAL VOLUME OF ABOVE EACH LAVER FROM THE DEPTHS GIVEN C***** IN THE INPUT DATA FILE. C***** REAL*B A,V,TV,ATOP,VDUM COMMON/VOL/ZMAX,DZ(40),Z(40),A(40),V(40),TV(40),ATOP(41),DBL COMMON/STEPS/DZLL,DZUL,MBOT,NM,NPRINT,MOAV,MONTH,ILAV,DV DZ(MBOT)=ZMAX-(Z(MBOT)+Z(MBOT-1»*.5 Z(MBOT)-ZMAX-DZ(MBOT)*.5 CALL LAKE(DZ(MBOT),VDUM,O,3) V ( MBOT) = VDUM AZ"OZ(MBOT) TV(MBOT)-V(MBOT) DO 10 I=1,MBOT-2 II"MBOT-I DZ(II)"Z(II+1)-(DZ(II+1)+Z(II)+Z(II-1»*.5 Z(II)=(DZ(II)+Z(II)+Z(II-1»*.51 290 AZ-AZ+DZ(II) CALL LAKE(AZ,VDUM,O,3) TV( II )=VDUM V(II)aTV(II)-TV(II+1) 10 CONTINUE TVPE DZ(1)=Z(2)-DZ(2)*.5 AZ=AZ+DZ(1) CALL LAKE(AZ,VDUM,O,3) TV(1)"VDUM V(1)=TV(1)-TV(2) RETURN END REAL DUMMV VARIABLE FOR DEPTH AS A SUM OF' DZ VALUES (USED BV LAKE) INTEGER*2 INDEX FOR LAVERS INTEGER*2 DUMMV INDEX FOR INDEXING LAVERS FROM THE BOTTOM REAL*8 DUMMV VARIABLE FOR VOLUME (USED BV LAKE) SUBROUTINE SETZ(MBOT) C***** C***** COMPUTE Z FROM DZ FOR EACH LAVER C***** REAL*B A,V,TV,ATOP COMMON/VOL/ZMAX,DZ(40),Z(40),A(40),V(40),TV(40),ATOP(41),DBL AZ-O. DO 100 1-1, MBOT Z(I)=AZ+DZ(I)*.5 AZ=AZ+DZ(I) 100 CONTINUE RETURN END TVPE REAL DUMMV VARIABLE FOR DEPTH AS ~ SUM OF DZ VALUES (USED BV LAKE) INTEGER*2 INDEX FOR LAVERS INTEGER*2 NUMBER OF LAVERS 315 SUBROUTINE SOLVE(VAR2,MBOT) 316 C***** 317 C***** TRI-DIAGONAL MATRIX SOLVING ROUTINE 318 C***** 319 REAL*8 AK,BK,CK,DK TT,TX(40) 320 COMMON/SOLV/ AK(40),BK(40),CK(40),DK(40) 321 DIMENSION VAR2(40) 322 DO 60 I-2,MBOT 323 TT=AK(I)/BK(I-1) 324 BK(I)-BK(I)-CK(I-1)*TT 325 60 DK(I)-DK(I)-DK(I-1)*TT 326 C**********BACK SUBSTITUTION************** 113 1 1 1 1 1 2 1 1 1 1 1 1 2 1 1 2 1 1 1 1 1 1 NAME I II K KI LW 327 328 329 TX(MBOT)=DK(MBOT)/BK(MBOT) DO 70 I=1,MBDT-1 330 70 331 .J"MBOT-I TX(.J)=(DK(.J)-CK(u)*TX(.J+1»/BK(.J) DO 80 I=1,MBOT 332 80 333 VAR2(I)"SNGL(TX(I» RETURN 334 END 335 SUBROUTINE SPLIT(I,LW) 336 C***** 337 C***** ROUTINE TO SPLIT THICK LAVERS (DZ > DZUL) INTO TWO OR MORE 338 C***** LAVERS OF EQUAL THICKNESS. ALL STATE VARIABLES ARE THE SAME IN 339 C***** EACH NEW LAYER. VOLUME IS AD.JUSTED LATER. 340 C***** 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 REAL*8 A,V,TV,ATOP COMMON/VOL/ZMAX,DZ(40),Z(40),A(40),V(40),TV(40),ATOP(41),DBL COMMON/CHOICE/MODEL ,NITRO, IPRNT(6) ,NDAVS,NPRNT(30) ,NCL ASS,PLOT(30) COMMON/RESULT/ T2(40),CHLATOT(40),PA2(40),PTSUM(40),B002(40), +OS02(40),C2(40),C02(40),XN02(40),XNH2(40),CHLA2(3,40), +PC2(3,40),XNC2(3,40),T20(40),SI2(40) COMMON/STEPS/OZLL,OZUL,MBOT,NM,NPRINT,MDAY,MONTH,ILAV,OV DO 100 KaI,MBOT II "MBOT+ I-K T2(II+1 )"T2(1I) C2( II+1)=C2( II) CD2(II+1)aC02(II) DO 50 Kl=1,NCLASS 50 CHLA2(KI t II+1)=CHLA2(KI,II) PA2(II+1J"PA2(II) B002(Il+1)=B002(Il) DS02(11+1)=DS02(11) IF(MODEL.EQ.3) THEN SI2(Il+1)=SI2(II) DO 51 KI=1,3 51 PC2(KI,Il+1)=PC2(KI,II) IF(NITRO.EQ.1) THEN DO 52 KI=1,3 52 XNC2(KI,II+1)aXNC2(KI,II) XN02(II+1)=XN02(Il) XNH2(II+1)=XNH2(11) END IF ENDIF DZ(II+1)=DZ(II) 100 CONTINUE MBOTaMBOT+1 DZ(I+1)-DZ(I)*0.5 DZ(I )-D2(I+1) IF(LW.GE.I) LW=LW+1 RETURN END TYPE INTEGER*2 INTEGER*2 INTEGER*2 INTEGER*2 INTEGER*2 INDEX OF LAVER TO BE SPLIT DUMMV INDEX TO INDEX LAYERS FROM THE BOTTOM INDEX OF LAVERS AT AND BELOW THE LAYER TO BE SPLIT INDEX OF ALGAL CLASS INDEX OF LAVER AT THE MIXED LAVER DEPTH 377 SUBROUTINE START(ST,S,FT,ISTART,INFLOW,MYEAR,ITER,IRUN,LEN8) 378 C***** 379 C***** ROUTINE TO READ THE INPUT DATA FILE FOR INITIAL 380 C***** CONDITIONS AND INPUT COEFFICIENTS 381 C***** 382 REAL*8 A,V,TV,ATOP 383 COMMON/RESULT/ T2(40),CHLATOT(40),PA2(40),PTSUM(40),BOD2(40), 384 +DS02(40) ,C2(40) ,CD2(40),XN02(40) ,XNH2(40),CHLA2(3,40) , 385 +PC2(3,40) ,XNC2(3,40),T20(40) ,SI2(40) 386 COMMON/TEMP/PARIO(4),PCDUM(3,40) ,XNHD(40),XNOD(40) , 387 + CHLADUM(3,40),XNCD(3,40),PADUM(40),SID(40) 388 COMMON/CHOICE/MODEL,NITRO,IPRNT(6),NDAVS,NPRNT(30),NCLASS,PLOT(30) 389 COMMON/CHANEL/WCHANL,ELCB,ALPHA,BW,WLAKE 390 COMMON/WATER/BETA,EMISS,XK1,XK2,HKMAX,WCOEF,WSTR 391 COMMON/STEPS/DZLL,DZUL,MBOT NM,NPRINT,MDAV,MONTH,ILAV,DV 392 COMMON/VOL/ZMAX,DZ(.40) ,Z(40) ,A(40), V(40), TV(40) ,ATOP(41) ,DBL 393 COMMON/FLOW/HMK(41),QE(40),FVCHLA(5),PE(5,41) 394 COMMON/YIELD/YCA,VCH02,Y2CH02,YCBOD,YPBOD,V2W,VPZP,YNZP,VZDO, 395 +VSCHL,VNHBOD,BRNO ,BRNH,XKNNH,THNNH,YPCHLA,BODK20, SB20 ,BRR 396 COMMON/PHYTO/PDEL(3),PMAX(3),PMIN(3),THR(3),THM(3),XKR1(3), 397 +XKR2(3) ,XKM(3) ,HSCPA(3) ,HSC1(3) ,HSC2(3) ,UPMAX(3), THUP(3), 398 +GROMAX(3), TMAX(3), TOPT(3) .XNMAX(3) .XNMIN(3),UNMAX(3), THUN(3), 399 +HSCN(3) ,HSCNH(3) ,XNDEL(3), IDIATOM',CHLMEAN, CHLMAX. SECCHI 114 400 COMMON/ZOOPL/IZ.MINOAY.MAXDAY.ZP.ZPMIN.PRMIN.PRMAX.PREDMIN.XIMIN. 401 +XIMAX.XKRZP.GRAZMAX(3).THGRAZ(3).ASM.THRZP.HSCGRAZ(3).CHLAMIN(3). 402 +REPRO.XI.XKMZ.GRAZE(3.40) 403 COMMON/FILE/ DIN.MET.FLO.TAPEB.TAPE1.IREC 404 CHARACTER*16 DIN.MET.FLO.TAPEB.TAPE1 405 CHARACTER*1 T1(16) 406 DIMENSION PSAT(3).XNSAT(3) 407 EQUIVALENCE (T1(1).TAPE1) 40B NITRO-O 409 DO 200 I~1.30 410 200 NPRNT(I)-O 411 C********* INPUT MODEL OPTIONS AND INITIAL CONDITIONS ****** 412 READ(7.*) MODEL.NCLASS.IDIATOM 413 READ(7.*) NITRO. ITER 414 WRITE(*.1006) 415 READ(*./(A)/) X 416 IF(X.EQ./V' .OR. X.EQ.'V/) THEN 417 IPRNT(2)=1 41B TAPE1=TAPE8 419 T1(LEN8+2)./O' 420 T1(LEN8+3)='U' 421 OPEN(2.FILE=TAPE1.STATUS./NEW') 422 ENDIF 423 WRITE(*.1000) 424 READ( *.1001) X 425 IF(X.EQ.'V' .OR. X.EQ./V/) THEN 426 TAPE1=TAPEB 427 T1(LEN8+2)./P I 428 T1(LEN8+3)=/L ' 429 OPEN (1.FILE=TAPE1.ACCESS='DIRECT'.STATUS='NEW'.RECL=4) 430 IREC=1 431 IPRNT(5)-1 432 WRITE(*.1004) 433 READ(*.*.ERR=202) (PLOT(I)'.I=1.10) 434 202 IPRNT(6)=I-1 435 WRITE(1.REC=IREC) REAL(IPRNT(6» 436 IREC=IREC+1 437 DO 201 I-1.IPRNT(6) 438 WRITE(1.REC=IREC) PLOT(I) 439 201 IREC=IREC+1 ENDIF WRITE(*.1002) READ( *.1001) X IPRNT(4)=1 IF(X.EQ.'Y ' .OR. X.EQ./Y/) IPRNT(4)-0 1000 FORMAT!' PLOT FILE TO BE CREATED (V/N) 7 1001 FORMAT Ai) . 1002 FORMAT I TEMPERATURE SIMULATION ONLY (Y/N) 7 1004 FORMAT I ENTER UP TO 10 DEPTHS TO BE SAVED './. +, END WITH A CHARACTER (N.N.N •...• X): '.\) 1006 FORMAT(' OUTFLOW FILE TO BE CREATED (V/N) 7 READ(7.*) NDAYS I, \) I, \) I. \) 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 45B 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 IF(NDAYS.GT.O) READ(7.*)(NPRNT(I).I=1.NDAYS) READ(7.*) DZLL,DZUL,BETA,EMISS.XK1.XK2,HKMAX.WCOEF.WSTR READ(7,*) WCHANL,WLAKE.DBL,ST.S.FT READ(7,*) ELCB,ALPHA.BW READ(7,*) MBOT,NM.NPRINT,INFLDW,DY,MONTH.ISTART,MYEAR IF(IPRNT(5).GT.0) THEN WRITE(i.REC=IREC) REAL(MONTH) IREC"IREC+i WRITE(i,REC-IREC) REAL(NM) IREC-IREC+1 WRITE(i.REC=IREC) REAL(MVEAR) IREC-IREC+i WRITE(i.REC=IREC) REAL(MODEL) IREC"IREC+1 WRITE(1.REC=IREC) REAL(NCLASS) IREC~IREC+1 WRITE(i.REC=IREC) IREC=IREC+1 WRITE(1,REC-IREC) IREC=IREC+1 ENDIF REAL(NITRO) REAL(IPRNT(4» READ(7.*! !Z(I),I=1,MBOT) READ(7.* T2(I),I-1,MBOT) READ(7.* C2(I) I-i,MBDT) READ(7.* CD2(I).I c i,MBOT) DO 100 K=1,NCLASS 100 READ(7,*) (CHLA2(K,I),I.1 t MBOT) READ(7,*) (PA2(I),I=1,MBOTJ READ(7,.) (BOD2(I),I-1,MBOT) READ(7.*) (DS02(I).I-i,MBDT) IF(IDIATOM.EQ.3) READ(7,*) (SI2(I),I-1,MBOT) IF(NITRO.EQ. 1) THEN READ(7,*) (XNH2(I),I-1,MBOT) READ(7.*) (XN02(I),I-1,MBOT) ENDIF 1500 FORMAT( /. 5X, 115 1 1 1 1 2 2 1 1 2 2 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 +/,5X,41HMINIMUM THICKNESS OF EACH LAYER (DZLL) = ,F5.2,9H METER(S) +/,5X,41HMAXIMUM THICKNESS OF EACH LAYER (DZUL) = ,F5.2,9H METER(S) +/,5X,40HSURFACE ABSORPTION COEFFICIENT (BETA) = ,F5.2, +/,5X,30HEMISSIVITY OF WATER (EMISS) = ,F5.2, +/,5X, 'EXTINCTION COEFF. OF WATER (XK1) = ',F5.2,3X,'M**-l', +/,5X,'EXTINCTION COEFF. OF CHLA (XK2) = ',F5.2,3X,'L/MG/M', +/,5X,'MAX. HVPOLIMNETIC DIFFUSIVITY (HKMAX) = ',F7.4,' M**2/D', +/,5X, 'WIND FUNCTION COEFFICIENT (WCOEF) = ',F6.3, +/,5X,'WIND SHELTERING COEFFICIENT (WSTR)= ',F6.3,/, +/,5X,34HWIDTH OF INLET CHANNEL (WCHANL) = ,F6.2,9H METER(S» 1501 FORMAT(5X,'LONGITUDINAL LENGTH OF LAKE (WLAKE) =',Fl0.2,7H METERS, +/,5X,30HDEEPEST BED ELEVATION (DBL) = ,F8.2,17H METERS ABOVE MSL, +/,5X,26HINITIAL LAKE STAGE (ST) = ,F8.2,17H METERS ABOVE MSL, +/,5X,16HBED SLOPE (S) • ,Fl0.S, +/,5X,29HROUGHNESS COEFFICIENT (FT) = ,F6.4, +//,5X,48HELEVATION OF BOTTOM OF OUTFLOW CHANNEL (ELCB) = ,F6.2, +17H METERS ABOVE MSL, +/,5X,40HSIDE-SLOPE OF OUTFLOW CHANNEL (ALPHA) = ,F6.2,8H DEGREES, +/,5X,31HBOTTOM WIDTH OF CHANNEL (BW) = ,F6.2,7H METERS, +//,5X,34HINITIAL NUMBER OF LAVERS (MBOT) = ,12, +/,5X,37HNUMBER OF MONTHS OF SIMULATION (NM) = ,12, +/,5X,54HDAV OF MONTH OF THE FIRST DAY OF SIMULATION (ISTART) = ,12 +/,5X,53HINTERVAL AT WHICH RESULTS WILL BE PRINTED (NPRINT) = ,13, +7H DAV(S» C********** INPUT PARAMETERS FOR BIOLOGICAL ROUTINES ***** READ(7,*) BODK20,SB20,BRR,FVCHLA(NCLASS+l) IF(MODEL.EQ.l) THEN REAO(7,*)VCA,PMAX(1),PMIN(1),XKM(1) READ(7,*)VCH02,V2CH02,YCBOO,VPBOO,YNHBOD,YPZP,YZDO, + YZW,YNHZP IF(IRUN.GT.O) THEN CALL RESTART 4 WRITE(*,30oo)VCA,PMAX(1),PMIN(1),XKM(1) WRITE(*,2) READ(*,*,ERR=3) ID,VAR IF(ID.EQ.l) YCA=VAR IF(ID.EQ.2) PMAX(l)=VAR IF(ID.EQ.3) PMIN(1)=VAR IF(ID.EQ.4) XKM(1)=VAR GOTO 4 ENDIF 3 WRITE(S, 1500) DZLL,DZUL,BETA,EMISS,XK1,XK2,HKMAX,WCOEF,WSTR, + WCHANL WRITE(S, 1501)WLAKE,DBL,ST,S,FT,ELCB,ALPHA,BW,MBOT,NM,I START, + NPRINT ' WRITE(8,2000)BOOK20,SB20,BRR,FVCHLA(NCLASS+l) 2 FORMAT(/ /1X, 'ENTER CHANGES: INTEGER,NEW VALUE".I, + 5X,'(ENTER ANV CHARACTER FOR NO CHANGES): ',\) 3000 FORMAT(10(/),1X, 'MODEL 0 PARAMETERS',//5X,'l : YCA = ',F7.4,/5X, + '2 ; PMAX + ',F7.4,/5X,'3 : PMIN = ',F7.4,/5X,'4 : XKM =' + F7.4) WRITE(S, 1950)VCA,PMAX(1),PMIN(i),XKM(1) 1950 FORMAT(/,2X,'BIOLOGICAL COEFFICIENTS',/,2X,23('-'),//,10X, + 'CARBON-CHLOROPHVLL RATIO',5X,F5.0,/,10X,'MAX. NUTRIENT I + 'SATURATED PHOTOSYNTHETIC RATE',3X,F5.3,' /DAV',/,10X, + 'MINIMUM CELL QUOTA FOR PHOSPHORUS',3X,F5.3,/,10X, + 'MORTALITV RATE',3X,F6.3,' /DAY',/) RETURN ENDIF IF(MODEL.EQ.2) THEN READ(7,*) HSCPA(1),HSC1(1),HSC2(1),GROMAX(1),TMAX(1),TOPT(1) READ(7,*) XKR1(1),XKR2(1),THR(1),XKM(1),THM(1),FVCHLA(1) ELSE DO 501 I=1,NCLASS 501 READ(7,*) UPMAX(I),THUP(I),GROMAX(I),TMAX(I),TOPT(I), + XKR1(I),XKR2(I),THR(I),XKM(I),THM(I) DO 502 I=1,NCLASS 502 READ(7,*)HSCPA(I),HSC1(I),HSC2(I),PMIN(I),PMAX(I),PSAT(I), + , FVCHLA(I) DO 500 K=1,NCLASS PDEL(K)=PMAX(K)-PMIN(K) DO 500 I=1,MBOT PC2(K,I)=CHLA2(K,I)*PSAT(K) 500 CONTINUE IF(NITRO.EQ.1) THEN DO 53 K=1,NCLASS 53 READ(7,*)XNMAX(K) ,XNMIN(K) ,XNSAT(K),UNMAX(K),THUN(K), + HSCN(K),HSCNH(K) READ(7,*) BRNO,BRNH,XKNNH,THNNH DO 51 K=1,NCLASS XNOEL(K)=XNMAX(K)-XNMIN(K) DO 51 I=1,MBOT XNC2(K,I)=CHLA2(K,I)*XNSAT(K) 51 CONTINUE ENDIF ENOIF READ(7,*)ZP,ZPMIN,PRMIN,PRMAX,PREDMIN,XIMIN,XIMAX,MINDAV,MAXDAY, + XKRZP,THRZP,REPRO,ASM 116 ." NAME FT 10 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 INFLOW IRUN ISTART ITER LEN8 MYEAR PSAT S ST VAR X XNSAT YNHZP DO 56 K=1,NCLASS 56 READ(7,*)GRAZMAX(K),THGRAZ(K),HSCGRAZ(K),CHLAMIN(K) REAO(7,*)YCHD2,Y2CH02,YCBOD,YPBOD,YNHBOD,YPZP,YZDO, +YZW,YNHZP IF(IRUN.GT.O) CALL RESTART WRITE(8,1500) DZLL,OZUL,BETA,EMISS,XK1,XK2,HKMAX,WCOEF,WSTR, +WCHANL WRITE(8,1501)WLAKE,DBL,ST,S,FT,ELCB,ALPHA,BW,MBOT,NM,ISTART, +NPRINT WRITE(8,2000)BODK20,SB20,BRR,FVCHLA(NCLASS+1) IF(MODEL.EQ.2) THEN WRITE (8,2202) WRITE(8,2203) HSCPA(1),HSC1(1),HSC2(1),GROMAX(1),TMAX(I), + TOPT(1),XKR1(1),XKR2(1),THR(1),XKM(1),THM(1) ELSE WRITE (8,2009) DO 601 I=1,NCLASS WRITE(B,2015) UPMAX(I),THUP(I),GROMAX(I),TMAX(I),TOPT(I), XKR1(I),XKR2(I),THR(I),XKM(I),THM(I) . 601 + WRITE (8,2008) 602 + DO 602 I-1,NCLASS WRITE(8,2015)HSCPA(I),HSC1(I),HSC2(I),PMIN(I),PMAX(I),PSAT(I), FVCHLA(I) IF(NITRO.EQ.1) THEN WRITE (8,2204) DO 653 I=1,NCLASS 653 + WRITE(8,2015)XNMAX(K),XNMIN(K) ,XNSAT(K),UNMAX(K) ,THUN( K), HSCN(K),HSCNH(K) 656 WRITE(8,2205) WRITE(8,2015)BRNH,BRNO,XKNNH,THNNH ENDIF ENOIF WRITE (8,2220) WRITE(8,2225)ZP,ZPMIN,PRMIN,PRMAX,PREOMIN,XIMIN,XIMAX,MINDAY, +MAXOAY,XKRZP,THRZP,REPRO,ASM WRITE(8,2222) 00 656 ~·1,NCLASS WRITE(8,2221)GRAZMAX(K) ,THGRAZ(K),HSCGRAZ(K) ,CHLAMIN(K ) YPCHLA=YPBOD/YCBOD WRITE(8,2230)YCH02,Y2CH02,YCBOO,YPBOD,YNHBOD,YPZP,YZDO, +YZW,YNHZP XKMZ=O. RETURN 2000 FORMAT(//,4X,'BODK20',4X,'SB20',5X,'BRR',5X, 'FVBOD',/,3X, +'(1/DAY)',2(2X,'(GM/M2)'),3X,'(M/S)',/,2X,5(F7.3,2X» 2008 FORMAT(//,5X,'HSCPA',4X,'HSC1',6X,'HSC2',6X,'PMIN', +6X,'PMAX',6X,'PSAT',6X,'FVCHLA') . 2009 FORMAT(//,' BIOLOGICAL COEFFIC. (DIATOMS,GREENS,BLUE-GREENS)', +/,4X,'UPMAX',6X,'THUP',5X,'GROMAX',5X,'TMAX',6X, +, ' TOPT' ,6X, ' XKR l' , 6X, ' XKR2' , 7X, ' THR' ,8X, ' XKM' ,6X, ' THM' ) 2015 FORMAT(1X,12(F8.3,2X» 2202 FORMAT(//,3X,'BIOLOGICAL COEFFICIENTS',/,4X, +'HSCPA',5X,'HSC1',5X,'HSC2',3X,'GROMAX',3X,'TMAX',5X,'TOPT',5X, +' XKR 1 ' , 5X, ' XKR2' ,6X, ' THR' , 5X, ' XKM' , 6X, ' THM' , / ) 2203 FORMAT(1X,11(F7.3,2X» 2204 FORMAT(//,3X,'NITROGEN BIOLOGICAL COEFFICIENTS',//,3X,'NMAX', +7X, 'NMIN', 7X, 'NSAT' ,5X, 'UNMAX' ,4X, 'THUN' ,6X, 'HSCN', 7X, 'HSCNH') 2205 FORMAT(/,6X,'BRNH',5X,'BRNO',5X,'XKNNH',3X,'THNNH',/, +4X,2('(MG/M2)',2X),'(1/DAY)') 2220 FORMAT(///,1X,'ZOOPLANKTON PARAMETERS',/,3X,'ZP',7X,'ZPMIN', +4X,'PRMIN',4X,'PRMAX',2X,'PREOMIN',4X,'XIMIN',3X,'XIMAX',2X, +'MINDAY',2X, 'MAXDAY',2X,'XKRZP',4X,'THRZP',6X,'REPRO',5X,'ASM') 2221 FORMAT(1X,12(F8.4,2X» 2222 FORMAT(//,3X,'GRAZMAX',3X,'THGRAZ',3X,'HSCGRAZ',3X,'CHLAMIN') 2225 FORMAT(1X,F7.0,2X,6(F7.3,2X),1X,I3,5X,I3,3X,4(F7.4,2X» 2230 FORMAT(/,1X,'YIELD COEFFICIENTS',/,3X,'YCH02',3X,'Y2CH02',4X, +'YCBOD',4X,'YPBOD',4X,'YNBOD',5X,'YPZP',5X,'YZDO',5X,'YZW', +4X,'YNHZP',/,1X,9(F7.4,2X» END TYPE REAL MANNINGS FRICTION FACTOR OF INFLOW CHANNEL INTEGER*2 INDEX FOR CHANGE OF VARIABLE IN MODEL 1 RESTART OPTION INTEGER*2 NUMBER OF INFLOWS/OUTFLOWS ON THE INFLOW DATA FILE INTEGER*2 FLAG FOR INITIAL OR RESTARTED SIMULATION INTEGER*2 DAY DF THE FIRST MONTH WHEN SIMULATION BEGINS INTEGER*2 NUMBER OF ITERATIONS IN THE BIOLOGICAL-NUTRIENT RROUTINES INTEGER*2 . LENGTH NAME FOR OUTPUT FILE TAPE8 INTEGER*2 YEAR OF SIMULATION REAL LEVEL OF INTERNAL PHOSPHORUS SATURATION REAL SLOPE OF INFLOW CHANNEL REAL INITIAL STAGE (M) REAL DUMMY VARIABLE FOR CHANGE OF VARIABLE IN MODEL 1 RESTART REAL DUMMY VARIABLE FOR RESPONSES TO INTERACTIVE QUESTIONS REAL LEVEL OF INTERNAL NITROGEN SATURATION REAL MASS YIELD RATIO OF AMMONIUM TO ZOOPLANKTON 117 OPTION 645 SUBROUTINE RESTART 646 C***** 647 C***** ROUTINE TO RESET VARIABLES AS REQUESTED BY THE USER WHEN 648 C***** OPERATING THE RESTART OPTION. 649 C***** 650 COMMON/WATER/BETA,EMISS,XK1,XK2,HKMAX,WCOEF,WSTR 651 COMMON/PHYTO/PDEL(3),PMAX(3),PMIN(3),THR(3) ,THM(3) ,XKR 1(3), 652 +XKR2(3),XKM(3),HSCPA(3),HSC1(3),HSC2(3).UPMAX(3).THUP(3), 653 +GROMAX(3),TMAX(3),TOPT(3),XNMAX(3) ,XNMIN(3) ,UNMAX(3) , THUN(3), 654 +HSCN(3) ,HSCNH(3) ,XNDEL(3),IDIATOM,CHLMEAN,CHLMAX,SECC HI 655 COMMON/YIELD/YCA,YCH02,Y2CH02,YCBOD,YPBOD,YZW,YPZP,YNZP,YZDO, 656 +YSCHL,YNHBOD,BRNO,BRNH,XKNNH,THNNH,YPCHLA,BODK20,SB20,BRR 657 COMMON/ZOOPL/IZ,MINDAY,MAXDAY,ZP,ZPMIN,PRMIN,PRMAX,PREDMIN,XIMIN, 658 +XIMAX,XKRZP,GRAZMAX(3),THGRAZ(3),ASM,THRZP,HSCGRAZ(3),CHLAMIN(3), 659 +REPRO,XI,XKMZ,GRAZE(3,40) 660 COMMON/FLOW/HMK(41),QE(40),FVCHLA(5),PE(5,41) 661 COMMON/CHOICE/MODEL ,NITRO, IPRNT(6) ,NDAYS,NPRNT(30) ,NCL ASS,PLOT(30) 662 DIMENSION WATR(6),PHTO(63),ZOPL(23),YLD(14),IC(4) 663 CHARACTER*1 C1,C2,C3,C4 664 EQUIVALENCE (BETA,WATR(1»,(PMAX(1),PHTO(1»,(ZP,ZOPL(1», 665 +(YCH02,YLD(1» 666 EQUIVALENCE(C1,IC(1»,(C2,IC(2».(C3,IC(3»,(C4,IC(4» 667 DATA IC/1661B,1665B, 16632,1664A/ 668 100 WRITE(*, 109) C1,C2,C3,C4 669 WRITE(*,105) 670 IF(MODEL.NE.1) WRITE(*,106) 671 IF(NITRO.GT.O) WRITE(*,107) 672 WRITE(*,108) 673 READ(*,*,ERR=110) ISET 674 GOTO (201,202,203,204,205,206) ISET 675 105 FORMAT(///,' SELECT CATEGORY FOR MODIFICATION',/10X, 676 +'1 = PHYSICAL AND MIXING COEFFICIENTS',/10X, 677 +'2 • SETTLING & OXYGEN DEPLETION COEFFICIENTS',/10X, 678 +'3 a YIELD COEFFICIENTS') 679 106 FORMAT(10X,'4 D PHYTOPLANKTON COEFFICIENTS',/10X, 680 +'5 = ZOOPLANKTON COEFFICIENTS') 681 107 FORMAT(10X,'6 = NITROGEN LIMITATION COEFFICIENTS') 682 108 FORMAT(10X,'X = END MODIFICATIONS AND BEGIN NEW SIMULATION') 683 109 FORMAT(1X,4A1) 684 201 WRITE~*'109) C1,C2.C3,C4 685 WRITE *,1)(WATR(I),Ia1,7) 686 WRITE *,2) 687 READ(*,*,ERR=100) ID,VAR 688 IF(IO.LT.1 .OR. ID.GT.7) THEN 689 PAUSE ' INTEGER IS OUTSIDE OF RANGE' 690 GOTO 201 691 ENDIF 692 WATR(ID)=VAR 693 GOTO 201 694 202 WRITE(*, 109) C1,C2,C3,C4 695 WRITE(*,4) BODK20,SB20.BRR,FVCHLA(NCLASS+1) 696 DO 212 1-1,NCLASS 697 212 WRITE(*,9) I+4,I.FVCHLA(I) 698 WRITE(*,2) 699 READ(*,*,ERR-100)ID,VAR . 700 IF(ID.LT.1 .OR. ID.GT.4+NCLASS) THEN 701 PAUSE' INTEGER IS OUTSIDE OF RANGE' 702 GOTO 202 703 ENDIF 704 IF(ID.EQ.1) 705 IF!ID.EQ.2l 706 IF ID.EQ.3 707 IF ID.EQ.4 708 IF ID.GT.4) 709 GOTO 202 BODK20-VAR SB20=VAR BRR=VAR FVCHLA(NCLASS+1)-VAR FVCHLA (10-4) =VAR 710 204 WRITE(*, 109) C1,C2.C3,C4 711 WRITE(*,3)(PHTO(I),PHTO(I+15),PHTO(I+30),I=1,9) 712 WRITE(*,8)(PHTO(I),PHTO(I+15),PHTO(I+30).I=10,15) 713 WRITE(*,2) 714 READ(*,*.ERR=100)ID,VAR 715 IF(ID.LT.1 .OR. ID.GT.4S) THEN. 716 PAUSE ' INTEGER IS OUTSIDE OF RANGE' 717 GO TO 204 718 ENDIF 719 PHTO(ID)-VAR 720 GOTO 204 721 206 WRITE(*, 109) C1,C2.C3,C4 722 WRITE(*.5)(PHTO(I).PHTO(I+9),I-46.54) 723 WRITE(*,2) 724 READ(*, *, ERR=100)ID, VAR 725 IF(ID.LT.46 .OR. ID.GT.5~) THEN 726 PAUSE ' INTEGER IS OUTSIDE OF RANGE' 727 GOTO 206 728 ENDIF 118 729 730 731 732 733 734 735 736 737 738 739 740 741 742 743 744 745 746 747 748 749 750 751 752 753 754 755 756 757 758 759 760 761 762 763 764 765 766 767 768 769 770 771 772 773 774 77S 776 777 778 779 780 781 782 783 784 785 786 787 788 789 790 791 792 793 794 795 796 797 798 799 800 801 802 803 804 80S 806 807 808 809 810 811 812 813 814 815 816 PHTO(ID) =VAR GOTO 206 205 WRITE(*,109) C1,C2,C3,C4 WRITE(*,6)(ZOPL(I),ZDPL(I+8),ZOPL(I+16),I=1,7),XKRZP,THRZP,MINDAY, +MAXDAY WRITE(*,2) READ(*,*,ERR=100)ID,VAR IF(ID.LT.1 .OR. ID.GT.25) THEN PAUSE ' INTEGER IS OUTSIDE OF RANGE' GOTO 205 ENDIF IF(ID.LT.24) THEN ZOPL (ID) =VAR ELSE IF(ID.EQ.24) MINDAY=ANINT(VAR) IF(ID.EQ.25) MAXDAY=ANINT(VAR) ENDIF GOTO 20S 203 WRITE(*,109) C1,C2,C3,C4 WRITE(*,7)(YLD(I),YLD(I+4),YLD(I+9),I=1,4),YLD(9),YLD(14) WRITE(*,2) READ(*,*,ERR=100)ID,VAR IF(ID.LT.1 .OR. ID.GT.14) THEN PAUSE' INTEGER IS OUTSIDE DF RANGE' GOTD 203 ENDIF YLD(ID)=VAR GDTO 203 110 RETURN 1 FDRMAT(10(/),1X,'1 : BETA = ',F6.3,/1X,'2 : EMISS = ',F6.S,/1X, +'3 : XK1 = ',F6.3,/1X,'4 : XK2 = ',F6.l,/lX,'5 : HKMAX = " +F6.3,/1X,'6 : WCOEFF= ',F6.1,/lX,'7 : WSTR = ',F6.3) 2 FORMAT(//1X,'ENTER CHANGES (ANY CHARACTER TO END):', +' INTEGER,NEW VALUE ',\) S FORMAT(2(/),1X,'PHYTOPLANKTON PARAMETERS',//2X,'1 : PMAX(1) " , +FS.2,6X,'16 : XKR2(1) " ',FS.3,6X,'31 : UPMAX(1) = ',FS.2,/2X, +'2 :'5X,'(2) " ',FS.2,6X,'17 :',5X,'(2) " ',F5.2,6X,'32 :',6X, +'(2) = , ,FS.2,/2X,'3 : ',SX,'(3) = ',FS.2,6X, '18 : ',SX, '(3) " " +FS.3,6X,'33 :' ,6X,'(3) " ',FS.2,/2X,'4 : PMIN(1) = ',F5.3,6X, +'19 : XKM(1) "',FS.3,6X,'34: THUP(1) = ',F5.3,/2X,'S :',5X, +'(2) • ',FS.3,6X,'20 :',4X,'(2) • ',FS.3,6X,'3S :',5X,'(2) =' +FS.3,/2X,'6 :',SX,'(S) = ',FS.3,6X,'21 :',4X,'(3) • ',FS.3,6X, +'36 :',SX,'(3) = ',FS.3,/2X,'7 : THR(1) "',FS.3,6X, +'22 : HSCPA(1)= ',FS.3,6X,'37 :GROMAX(1) " ',FS.3,/2X,'8 :',4X, +'(2) = ',FS.3,6X,'23 :', +6X, '(2)= I ,F5.3,6X, '38 :' ,6X,'(2) I: '/ ,F5.3,/2X,'9 ! I ,4X,'(S) =' +FS.3,6X,'24 :',6X,'(3)= ',F5.3,6X,'39 :';6X,'(S) = ',F5.3) 8 FORMAT(1X,'10: THM(1) = ',FS.3,6X,'2S : HSC1(1) • ',FS.O,6X, +'40: TMAX(1)" ',FS.1,/1X,'11 :',4X,'(2) = ',F5.3,6X, +'26 :',5X,'(2) = ',FS.O,6X,'41 :',SX,'(2) • ',FS.1,/1X,'12 :',4X, + ' (3 ) =' , F S . 3 , 6X, ' 27 :', SX, ' ( 3) = " F 5 . 0, 6X, ' 42 :', SX, ' ( 3) = " +FS.1,/1X,'13 : XKR1(1) = ',FS.3,6X,'28 : HSC2(1) • ',FS.O,6X, +'43 : TOPT(1) .. ',F5.1,/1X,'14 :',5X,'(2) = ',F5.3,6X,'29 :',SX, +'(2) = ',FS.O,6X,'44 :',SX,'(2) = ',FS.1,/1X,'15 :',SX,'(3) • " +FS.3,6X,'30 :',SX,'(3) = ',FS.O,6X,'4S :',SX,'(3) = ',FS.1) 4 FORMAT(11(/),1X,'REACTION ANO SETTLING RATE COEFFICIENTS', +//5X,'1 : BODK20 = ',F6.3,/SX,'2 : SB20 = ',F6.3, +/5X,'3 : BRR = ',F6.3,/SX,'4 : FVBOD ',F6.3) S FORMAT(8(/),1X,'NITROGEN PARAMETERS', +//5X,'46 : XNMAX(1)= ',FS.1,10X, +'SS : THUN(1) = ',FS.3,/5X,'47 :',6X,'(2)= ',FS.1,10X,'S6 :',SX, +' (2) = " FS. 3, /5X, '48 :', 6X, ' (3)" " FS. 1 , 1 OX, 'S7 :', SX, ' (3) = " +FS.3,/SX,'49 : XNMIN(1)" ',F5.2,10X,'S8 : HSCN(1) " ',F5.3,/SX, +'SO :',SX,'(2)" ',FS.2,10X,'S9 :',SX,'(2) = ',FS.3,/5X,'S1 :',5X, +'(3) .. ',FS.2,10X,'60 :',SX,'(3) = ',FS.3'/SX,'52 : UNMAX(1) = , +FS.3,9X,'61 : HSCNH(1)= ',FS.3,/5X,'S3 :',6X,'(2) = ',F5.3, +9X,'62 :',6X,'(2)= ',FS.3,/SX,'S4 :',6X,'(3) = ',FS.3, +9X, '63 :' ,6X,' (3)= " FS.3) 6 FORMAT(8(/),1X,'ZOOPLANKTON PARAMETERS',//5X,'1 : ZP',4X,'. " +F6.0,6X,'9 : GRAZMAX(1)" ',FS.4,4X,'17 :HSCGRAZ(1)= '.F5.4,/SX, +'2 : ZPMIN = ',F5.3.6X,'10 :'.8X.'(2)= ',FS.4,4X,'18 :'.7X, +'(2)= '.FS.4,/SX,'3 : PRMIN " ',FS.3,6X,'11 :'.8X,'(3)= ',F5.4, +4X,'19 :',7X,'(3)= ',FS.4,/SX,'4 : PRMAX" ',FS.3,6X, +'12 : THGRAZ(1) " ',FS.3.4X,'20 :CHLAMIN(1)= ',FS.3,/SX, +'S : PREDMIN=',FS.3,6X,'13 :',1X,'(2) = ',FS.3,4X,'21 :',7X, +'(2)= ',FS.3,/5X,'6 : XIMIN = ',FS.1,6X,'14 :',7X,'(3) = '.FS.3, +4X,'22 :',7X,'(3)= ',FS.3,/5X,'7 : XIMAX .. ',FS.1,6X,'15 :. ASM', +7X,'c ',FS.2,4X,'2S : REPRO',5X,'= ',FS.2,/SX,'8 : XKRZP = ',F5.3, +6X,'16 : THRZP',SX,'= ',F5.3,4X,'24 : MINDAY',4X,'= ',13, +/54X,'25 : MAXDAY'.4X,'· ',13) 7 FORMAT(11(/),2X,'NUTRIENT AND MASS RATIOS',//2X, +'1 : YCH02· ',F5.4.5X,'5 : VZW = ',F5.4,5X,'10 : VNHBOD = I +FS.4,/2X,'2 : Y2CH02= ',F5.4,SX,'6 : YPZP "',FS.4,SX, +'11: BRNO .. ',FS.4,/2X,'3: YCBOD = ',FS.4,SX,'7 : YNHZP = " +FS.4,5X,'12 : BRNH • ',FS.4,/2X,'4 : YPBOD = ',FS.4,SX, +'8 : YZOO .. ',FS.2,SX,'13 : XKNNH = ',FS.4,/24X,'9 : YSCHL = ' +FS.2,SX,'14 : THNNH .. ',F5.4) 9 FDRMAT(5X,I1,' : FVCHLA(',I1,') = ',F6.3) 119 817 END NAME TYPE C1 C2 C3 C4 IC ID ISET VAR CHAR*1 CHAR*1 CHAR*1 CHAR*1 INTEGER*2 INTEGER*2 INTEGER*2 REAL SCREEN CONTROL CHARACTER SCREEN CONTROL CHARACTER SCREEN CONTROL CHARACTER SCREEN CONTROL CHARACTER SCREEN CONTROL VARIABLE DUMMY INDEX FOR VARIABLE TO BE CHANGED INDEX FOR CATEGORY OF VARIABLE TO BE CHANGED DUMMY VARIABLE FOR NEW VALUE OF THE VARIABLE TO BE CHANGED 818 SUBROUTINE STATS(FLDATA,XX,IFLAG,DEPTH,I) 819 C***** 820 C***** COMPUTE STATISTICS AND STATISITICAL QUANTITIES WITH 821 C***** Y DESIGNATING FIELD DATA AND X DESIGNATING MODEL RESULTS. 822 C***** 823 CoMMON/STAT/SUMXY(10),SUMX(10),SUMY(10),XSQ(10), 824 +YSQ(10),RSQ(10),RMS(10),RELM(10),MTHREL(10),MDAYREL(1O),ZREL(10), 825 +ZRELM(10),RS(10),REL(10),MTHRMS(10),MDAYRMS(10),ZRS(1O),ZRMS(10) 826 CoMMoN/STEPS/DZLL,DZUL,MBOT,NM,NPRINT,MDAY,MONTH,ILAY,DY 827 SUMXY(I)=SUMXY(I)+FLDATA*XX 828 SUMX(I)=SUMX(I)+XX 829 SUMY(I)=SUMY(I)+FLDATA 830 XSQ(I)=XSQ(I)+XX*XX 831 YSQ(I)=YSQ(I)+FLDATA*FLDATA 832 XX=XX-FLDATA 833 X2 a XX/FLDATA*100. 834 X3 a ABS(XX) 835 IF(X2.GT.ABS(REL(I») THEN 836 REL(I)=X2 837 ZREL(I)=DEPTH 838 ENDIF 839 IF(X2.GT.ABS(RELM(I») THEN 840 RELM(I)=X2 841 MTHREL(I)-MoNTH 842 MDAYREL(I)=MDAY 843 ZREL(I)=DEPTH 844 ENDIF 845 IF(X3.GT.ABS(RS(I») THEN 846 RS(I)=XX 847 ZRS(I)=DEPTH 848 ENDIF 849 IF(X3.GT.ABS(RMS(I») THEN 850 RMS(I)=XX 851 MTHRMS(I)-MONTH 852 MDAYRMS(I)=MDAY 853 ZRMS(I)=DEPTH 854 ENDIF 855 IFLAG=IFLAG+1 856 RETURN 857 END NAME TYPE DEPTH REAL DEPTH OF FIELD DATA MEASUREMENT (M) FLDATA REAL FILED DATA MEASUREMENT USED IN COMPUTATIONS I INTEGER*2 .IFLAG INTEGER*2 INDEX OF STATE VARIABLE FOR WHICH STATISTICS ARE BEING COMPUTED NUMBER OF DATA POINTS USED IN THE STATISTICS ON EACH X2 X3 XX REAL REAL REAL STATE VARIABLE RELATIVE DEVIATION BETWEEN FIELD DATA AND MODEL ABSOLUTE VALUE OF DEVIATION BETWEEN FIELD DATA AND MODEL DEVIATION BETWEEN FIELD DATA AND MODEL 858 SUBROUTINE SUBLAY(IEUPH) 859 C .••.. 860 C ...•• SUBDIVIDE EACH LAYER IN EUPHOTIC ZONE INTO LAYERS OF 0.2 M OR LESS 861 C ....• FoR COMPUTING LIGHT LIMITATION ·oN GROWTH. 862 C ...•• 863 REAL*8 A,V,TV,AToP,SVoL,VDUM,TOV 864 COMMON/VOL/ZMAX,DZ(40),Z(40),A(40) ,V(40) ,TV(40),AToP(4 1),DBL 865 COMMON/SUB/SDZ(60) , SZ(60)·, LAY(40) , AVGI(4,60) ,SVOL(60) 866 NSLAY.O 867 DO 100 I-1,IEUPH 1 868 LAY(I)-INT(DZ(I)*5.+.5) 1 869 DO 200 J-1,LAY(I) 2 870 K-NSLAY+J 2 871 SDZ(K)=0.2 2 872 IF(J.EQ.1) SDZ(K)=DZ(I)-FLOAT(LAY(I)-1)*0.2 2 873 200 CONTINUE 1 874 NSLAY=NSLAV+LAY(I) 120 .. (.1 875 100 CONTINUE 876 C ..... ESTABLISH ARRAY OF DISTANCE OF EACH SUBLAYER BELOW WATER SURFACE 877 SUMSZ=O. 878 DO 300 ~-1,NSLAY 879 SZ(~)=SUMSZ+SDZ(~)/2. 880 SUMSZ=SUMSZ+SDZ(~) 881 300 CONTINUE 882 C ..... ESTABLISH ARRAY OF VOLUME OF EACH SUBLAYER .•... 883 CALL LAKE(ZMAX,VDUM,O,3) 884 TOV=VDUM 885 ZDUM-ZMAX 886 DO 400 I=1,NSLAY 887 ZDUM-ZDUM-SDZ(I) 888 CALL LAKE(ZDUM,VDUM,O,3) 889 SVOL(I)=TOV-VDUM 890 TOV-TOV-SVOL(I) 891 400 CONTINUE 892 RETURN 893 ENO NAME TYPE I IEUPH ~ INTEGER*2 INTEGER*2 INTEGER*2 INTEGER*2 INTEGER*2 REAL REAL*8 REAL*8 REAL INDEX OF LAYERS INDEX OF LAYER AT THE EUPHOTIC DEPTH INDEX OF SUBLAYERS WITHIN REGULAR LAYERS INDEX OF SUBLAYERS K NSLAY SUMSZ TOV VDUM ZDUM COUNTER FOR TOTAL NUMBER OF SUBLAYERS DUMMY VARIABLE OF COMPUTING DEPTH OF EACH SUBLAYER DUMMY VARIABLE FOR COMPUTING THE VOLUME OF EACH SUBLAYER VOLUME COMPUTED BY LAKE (M3) DEPTH USED BY LAKE TO COMPUTE VOLUME (M) 894 SUBROUTINE THICKNS(MBOT) 895 C***** " 896 C***** COMPUTE THICKNESS OF EACH LAYER FROM THE DEPTH 897 C** ••• AREA CURVE IN LAKE. 898 C ••••• 899 REAL.8 A,V,TV,ATOP,AVOL 900 COMMON/VOL/ZMAX,DZ(40),Z(40),A(40),V(40),TV(40),ATOP(41),DBL 901 AZ=O. 902 AVOL-O. 903 "DO 100 I-1,MBOT 904 II=MBOT+1-I 905 AVOL=AVOL+V(II) 906 CALL LAKE(ZDUM,AVOL,O,4) 907 DZ(II)=ZDUM-AZ 908 AZ=AZ+DZ(II) 909 100 CONTINUE 910 RETURN 911 END NAME TYPE AVOL AZ I II MBOT ZOUM VOLUME BELOW A LAYER (M3) REAL.8 REAL INTEGER.2 INTEGER.2 INTEGER.2 REAL DISTANCE OF A LAYER FROM THE BOTTOM (M) INDEX OF LAYERS DUMMY INDEX FOR DISTANCE FROM BOTTOM NUMBER OF LAYERS DEPTH COMPUTED IN LAKE eM) 912 SUBROUTINE TVOL(MBOT) 913 C ..... 914 C •••• * DETERMINE THE V~LUME OF WATER ABOVE A LAYER 915 C ..... 916 REAL"'8 A,V,TV,ATOP,SUM 917 COMMON/VOL/ZMAX,DZ(40),Z(40),A(40),V(40),TV(40),ATOP(41),DBL 918 SUMo:O. 919 DO 100 I=1,MBOT 920 SUM-SUM+V(I) 921 TV(I)-SUM 922 100 CONTINUE 923 RETURN 924 END 925 SUBROUTINE WDEPTH(ST,Q,LW) 926 C.* ..... 927 C."'''''''. DETERMINE LAYERS CONTIBUTING TO OUTFLOW. OUTFLOW FROM 928 C ........... A LAYER IS PROPORTIONAL TO THE THICKNESS OF THE LAYER. 929 C****. COMPUTE AND SEND OUTFLOW CONCENTRATIONS TO THE OUTF.LDW 930 C*"'*"'* DATA FILE (TAPE8.0UT). 121 REAL*8 A,V,TV,ATOP COMMON/CHANEL/WCHANL,ELCB,ALPHA,BW,WLAKE COMMON/RESULT/ VAR(40,21) COMMON/VOL/ZMAX,DZ(40),Z(40),A(40),V(40),TV(40),ATOP(41),DBL COMMON/CHOICE/MODEL ,NITRO , IPRNT(6) ,NDAYS,NPRNT(30) ,NCL ASS,PLOT(30) COMMON/STEPS/DZLL,DZUL,MBOT,NM,NPRINT,MDAY,MONTH,ILAY,DY COMMON/FLOW/HMK(41),QE(40),FVCHLA(5),PE(5,41) DIMENSION VOUT(11),IFL(11) DATA IFL/1,3,5,6,7.8,9,10,11.12,13/ Qp.-Q ERR"QP*0.001 DH=ST-ELCB AREA=(BW+DH/TAN(ALPHA*0.017453293»*DH X=Q/(86400*AREA) VH=X*X* .05097 931 932 933 934 935 936 937 938 939 940 941 942 943 944 945 946 947 948 949 950 951 952 953 954 955 956 957 958 959 960 961 962 963 964 965 966 967 968 969 970 971 972 973 974 975 976 977 978 979 980 981 982 983 984 985 986 987 988 989 990 991 111 QPRIME=O.O SUMZ=O.O SUMV=O. SUMVRHO=O. 100 DO 100 I=1,MBOT IF(V(I).LE.QE(I» GOTO 100 RHOWL=RHO(VAR(I.1),VAR(I.7).VAR(I.8» SUMZ-SUMZ+DZ(I) SUMV=SUMV+V(I) SUMVRHO=SUMVRHO+V(I)*RHoWL RHOM=SUMVRHO/SUMV IF(RHOWL.LE.RHOM)GO TO 100 ZW=DH+RHOM*VH/(RHOWL-RHOM) IF(SUMZ.LE.ZW) GO TO 100 LW"I GOTO 112 CONTINUE LW-MBOT 112 DO 110 I=1,LW IF~QE(I).GE.V(I» GOTO 110 QE I)=QP*DZ(I)/SUMZ+QE(I) IF QE(I).GT.V(I» THEN QPRIME=QE(I)-V(I)+QPRIME QE(I )-V(I) ENDIF 110 CONTINUE IF(QPRIME.GT.ERR) THEN QP=QPRIME GOTO 111 ENDIF C*********** COMPUTATIONS FOR OUTFLOW FILE *********** IF(IPRNT(2).GT.0) THEN . C**~********* CONVERT OUTFLOW VOLUME TO CFS ********** QOUT=-Q*4.0873E-4 1 1 2 120 1002 1001 RQIN=-1/Q DO 120 \J=1,1 VOUT(\J)=O.O DO 120 I=1,MBDT VOUT(u)=VAR(I.IFL(\J»*QE(I)*RQIN+VOUT(\J) WRITE(2,1002) MONTH.MDAY.QOUT WRITE(2,1001) (VOUT(\J),\J=1.11) ENDIF RETURN END FORMAT(1X.2I3.F8.1) FORMAT(1X.F4.1,F7.4,F6.1.F5.1,F6.1.F7.4,5F7.4) NAME TYPE AREA REAL DH REAL G REAL I INTEGER*2 IFL INTEGER*2 \J INTEGER*2 LW INTEGER*2 NUM INTEGER*2 Q REAL QOUT REAL RHDM REAL RHDWL REAL ST REAL SUMV REAL SUMVRH REAL SUMZ REAL VH REAL VOUT REAL X REAL ZW REAL CROSS SECTIONAL AREA OF INFLOW DEPTH OF INFLOW GRAVITATIONAL ACCELERATION (9.81 M/S2) INDEX OF LAYERS INDEX ARRAY OF STATE VARIABLE POSITIONS IN COMMON/RESULT/ INQEX ON IFL INDEX OF LOWEST LAYER CONTRIBUTING TO OUTFLOW NUMBER OF STATE VARIABLES IN MODEL OUTFLOW (M3/DAY) OUTFLOW (CFS) . MEAN DENSITY IN THE LAYERS CONTRIBUTING TO OUTFLOW DENSITY OF NEXT LAYER TO BE CONSIDERED FOR CONTIBUTING TO OUTFLOW STAGE (M) . TOTAL VOLUME IN THE LAYERS CONTRIEUTING TO OUTFLOW SUM OF VOLUME AND DENSITY IN THE LAYERS CONTRIBUTING TO OUTFLOW TOTAL DEPTH IN THE LAYERS CONTRIBUTING TO OUTFLOW VELOCITY HEAD OF INF.LOW STATE VARIABLES LISTED TO THE OUTFLOW DATA FILE INFLOW VELOCITY DEPTH CONTRIBUTING TO OUTFLOW 122 f·' 992 SUBROUTINE WINEN(T,V,WIND) 993 C .. . 994 C ... CALCULATION OF THE SHEAR VELOCITY AND THE WIND SHEAR STRESS 995 C ... CONVERSION OF WIND SPEED FROM M.P.H. TO M/S 996 C ..• DENSITY OF WATER ANO AIR ASSUMED TO BE 1000 AND 1.177 KG/M3 997 C ... 998 W=WIND*0.447 999 CALL LAKE(FTCH,O.O,O,2) 1000 ZB-ALOG(FTCH)*0.8-1.0718 1001 WKW*1.666667*(ZB+4.6052)/(ZB+9.2103) 1002 C ..• CALCULATION OF WIND SHEAR STRESS 1003 CZ-SQRT(W)*.0005 1004 IF(W.GE.15.) CZ-.0026 1005 T=1.177*CZ*W*W 1006 C .•. ASSIGNMENT OF CHARACTERISTIC SURFACE VELOCITY 1007 C ... USING CALCULATION OF SHEAR VELOCITIES 1008 V=.0343*SQRT(CZ)*W 1009 RETURN 1010 END NAME TYPE CZ REAL WINO SHEAR STRESS COEFFICIENT FTCH REAL FETCH (M) RHOA REAL DENSITY OF AIR (1.177 KG/M3) RHOW REAL DENSITY OF WATER (1000 KG/M3) SVEL REAL SHEAR VELOCITY (M/S) T REAL WIND SHEAR STRESS (N/M2) V REAL CHARACTERISTIC SURFACE VELOCITY (M/S) W REAL WINO REAL WIND SPEED ~M/S~ WIND SPEED MPH ZB REAL WIND FETCH FUNCTION 1011 SUBROUTINE WINMIX( E,TS,IL,MBOT,) 1012 C .. . 1013 C ... CALCULATE THE AMOUNT OF ENTRAINMENT RESULTING FROM WINO MIXING. 1014 C ... USE THE DEPTH OF CENTER OF MASS OF MIXED LAYER TO DETERMINE THE 1015 C ... POTENTIAL ENERGY THAT MUST. BE OVERCOME BY THE KINECTIC ENERGY 1016 C .•. OF THE WIND FOR ENTRAINMENT TO OCCUR. 1017 C ... 1018 REAL*8 A,V,TV,ATOP 1019 COMMON/VOL/ZMAX,DZ(40),~(40),A(40),V(40),TV(40),ATOP(41),DBL 1020 COMMON/RESULT/ T2(40),CHLATOT(40),PA2(40),PTSUM(40),BOD2(40), 1021 +DS02(40),C2(40),CD2(40),XN02(40),XNH2(40),CHLA2(3,40), 1022 +PC2(3,40),XNC2(3,40),T20(40),SI2(40) 1023 IF(IL.GE.MBOT) GO TO 35 1024 SUM1=O. 1025 SUM2-0. 1026 DO 10 I=1,IL 1027 ILAY=I 1028 RV=RHO(T2(I),C2(I),CD2(I»*V(I) 1029 SUM1-SUM1+RV 1030 10 SUM2-SUM2+RV*Z(I) 1031 I=ILAY 1032 20 DCM=SUM2/SUM1 1033 TSTEP-T2(I+1) 1034 C ... CALCULATION OF POTENTIAL ENERGY OF MIXED LAYER 1035 PE=9.81*TV(I)*(Z(I)+(DZ(I)/2.)-DCM)*(RHO(TSTEP,C2(I+1),CD2(I+1» 1036 +-RHO(TS,C2(I),CD2(I») 1037 C ... CRITERIA FOR ENTRAINMENT 1038 IF( E.LT.PE) GO TO 40 1039 C.,.ENTRAINMENT OF LAYER 1+1 1040 1=1+1 1041 TS=(TS*TV(I-1)+TSTEP*V(I»/TV(I) 1042 IF(I.GE.MBOT) GO TO 40 1043 RV-RHO(T2(I),C2(I),CD2(I»*V(I) 1044 SUM1-SUM1+RV 1045 SUM2=SUM2+RV*Z(I) 1046 GO TO 20 1047 35 I"IL 1048 40 IL-I 1049 DO 50 K"1,IL 1050 50 T2(K)-TS 1051 RETURN 1052 END NAME TYPE DCM E GRA I IL ILAY REAL REAL REAL INTEGER*2 INTEGER*2 INTEGER*2 DEPTH OF CENTER OF MASS OF THE MIXED LAYER (M) KINECTIC ENERGY OF WINO GRAVITATIONAL ACELERATION (9.81 M/S) INDEX OF LAYERS INDEX OF LAYER AT THE INITIAL MIXED LAYER DEPTH DUMMY INDEX FOR MIXED LAYER DEPTH 123 Il.D K MBOT PE INTEGER*2 INTEGER*2 INTEGER*2 REAl. INDEX OF l.AYER AT THE FINAl. MIXED l.AYER DEPTH INDEX OF l.AYERS NUMBER OF l.AYERS POTENTIAl. ENERGY REQUIRED TO l.IFT THE NEXT l.AYER BEl.OW RV SUM1 SUM2 TS TSD TSTEP REAl. REAl. REAl. REAl. REAl. REAl. THE MIXED l.AYER TO THE CENTER OF MASS OF THE MIXED l.AYER DENSITY OF NEXT l.AYER TO BE ENTRAINED SUM OF DENSITY OF l.AYERS IN THE MIXED l.AYER SUM OF PRODUCT OF DENSITY AND DEPTH OF l.AYERS IN THE MIXED l.AYER INITIAl. MIXED l.AYER TEMPERATURE (DEGREES C) FINAl. MIXED l.AYER TEMPERATURE (DEGREES C) TEMPERATURE OF NEXT l.AYER BEl.OW THE MIXED l.AYER 1053 SUBROUTINE VOl.UME(MBOT) 1054 C***** 1055 C***** COMPUTE THE VOl.UME OF EACH l.AYER BASED ON THE DEPTH-VOl.UME 1056 C***** REl.ATIONSHIP FOUND IN l.AKE. 1057 C***** 1058 REAl.*8 A,V,TV,ATOP,VDUM,VZ 1059 COMMON/VOl./ ZMAX,DZ(40),Z(40),A(40),V(40),TV(40),ATOP(41),DBl. 1060 AZZ=O. 1061 CAl.l. l.AKE(ZMAX,VDUM,O,3) 1062 VZ=VDUM 1063 DO 100 I=1,MBOT-1 1064 AZZ=DZ(I)+AZZ 1065 Z2=ZMAX-AZZ 1066 CAl.l. l.AKE(Z2,VDUM,O,3) 1067 V(I)= VZ-VDUM 1068 100 VZ=VDUM 1069 V(MBOT)=VZ 1070 RETURN 1071 END NAME TYPE AZZ I REAl. INTEGER*2 DUMMY VARIABl.E FOR DEPTH TO BOTTOM OF l.AYERS INDEX OF l.AYERS MBOT INTEGER*2 NUMBER OF l.AYERS VDUM REAl.*8 DUMMY VARIABl.E FOR VOl.UME FROM l.AKE VZ Z2 REAl.*8 DUMMY VARIABl.E FOR VOl.UME TO TOP OF PREVIOUS l.AYER DISTANCE FROM THE BOTTOM TO BOTTOM OF THE CURRENT l.AYER REAl. 1072 1073 1074 1075 1076 1077 1078 1079 1080 1081 1082 1083 1084 1085 1086 1087 1088 1089 1090 1091 1092 1093 1094 1095 1096 1097 1098 1099 1100 1101 1102 1103 1104 1105 1106 1107 1108 1109 1110 1111 1112 1113 SUBROUTINE ZPl.KTN(DY,TD,MBOT,XK,GRTOT) C***** C***** COMPUTE ZOOPANKTON GRAZING AS A FUNCTION OF CHl.ORO-A C***** CONCENTRATION AND MIGRATION TIME SPENT IN EACH l.AYER. C***** TIME IN EACH LAYER PROPORTIONAl. TO CHLORO-A CONCENTRATION C***** AND l.ENGTH OF NON-DAYl.IGHT PERIOD. C***** COMPUTE ZOOPl.ANKTON DAYDEPTH FROM DISSOl.VED OXYGEN PROFIl.E. C***** COMPUTE ZOOPl.ANKTON PREDATION AS A FUNCTION OF l.IGHT C***** INTENSITY AT THE DAYDEPTH AND A l.INEAR TIME VARYING C***** PREDATION RATE. C***** REAl.*8 A,V,TV,ATOP COMMON/CHOICE/MODEL ,NITRO , IPRNT(6) ,NDAYS,NPRNT(30) ,NCl. ASS,PLOT(30) COMMON/ZOOPl./IZ, MINDAY ,MAXDAY , ZP, ZPMIN,PRMIN,PRMAX,PRE DMIN,XIMIN, +XIMAX,XKRZP,GRAZMAX(3),THGRAZ(3),ASM,THRZP,HSCGRAZ(3),CHl.AMIN(3), +REPRO,XI,XKMZ,GRAZE(3,40) . COMMON/TEMP/PARIO(4) ,PCDUM(3,40),XNHD(40) ,XNOD(40) , + CHl.ADUM(3,40),XNCD(3,40),PADUM(40),SID(40) COMMON/RESUl.T/ T2(40),CHl.ATOT(40),PA2(40),PTSUN(40),BOD2(40), +DS02(40) tC2(40) ,C02(40) ,XN02(40) ,XNH2(40) ,CHl.A2(3,40) , +PC2(3,40J,XNC2(3,40),T20(40),SI2(40) COMMON /VOl./ ZMAX,DZ(40),Z(40),A(40),V(40),TV(40),ATOP(41),DBl. COMMON/YIEl.D/YCA,YCH02,Y2CH02,YCBOD,YPBOD,YZW,YPZP,YNZP,VZDO, +YSCHL ,YNHBOD,BRNO,BRNH,XKNNH,THNNH,YPCHLA,BODK20, SB20, BRR COMMON/PHYTO/PDEL(3),PMAX(3),PMIN(3),THR(3),THN(3),XKR1(3). +XKR2(3),XKM(3),HSCPA(3),HSC1(3),HSC2(3),UPMAX(3),THUP(3), +GROMAX(3),TMAX(3),TOPT(3),XNMAX(3),XNMIN(3),UNMAX(3),THUN(3), +HSCN(3).HSCNH(3) ,XNDEL(3) ,IDIATOM,CHLMEAN,CHLNAX,SECCH I C********** DETERMINATIONOF ZOOPLANKTON DAY DEPTH********** DO 201 I=1,MBOT IF(DS02(I).LT.0.5) THEN IZ=I-1 GOTO 200 ENDIF 201 CONTINUE IZ-MBOT 200 IF(IZ.EQ.O) IZ=1 ZP-ZP/V(IZ) CHLASUM=O.O DO 99 I-1.IZ 99 CHLASUM-CHLASUM+CHLATOT(I)*V(I) GRTOT·O. 124 .. 1114 1115 1116 1117 1118 1119 1120 1121 1122 1123 1124 1125 1126 1127 1128 1129 1130 1131 1132 1133 1134 1135 1136 1137 1138 1139 1140 1141 1 1142 1 1143 1 1144 2 1145 2 1146 2 1147 2 1148 2 1149 2 1150 2 1151 2 1152 1 1153 1 1154 1 1155 1 1156 1 1157 1 1158 1 1159 1 1160 1161 1162 1163 1164 1165 NAME IF(IZ.GT.MBOT) THEN DO 97 I"'IZ+1,MBOT 97 GRAZE(K,I)=O.O ENDIF C*********DETERMINATIDN OF DAYTIME PREDATION*************** C************ZP IN N/M**3 AS N/V(IZ)*********************** DDAY-(DV-MINDAV)/(MAXDAV-MINDAY) IF(DDAY.LT.O.) DDAY~O.O IF(DDAY.GT.1.) DDAY-1.0 PREDMAX-AMIN1(PRMAX,(PRMIN+(PRMAX-PRMIN)*DDAY» DI-(XI-XIMIN)/(XIMAX-XIMIN) IF(DI.GT.1.) 01-1. IF(DI.LT.O.O) 01-0.0 XKMZ-PREDMAX*DI*(ZP-ZPMIN)*TD*.04167 IF(XKMZ.LT.O.O) XKMZ-O.O XK=XKMZ*V( IZ) X"1.0 IF(NCLASS.EQ.3) THEN x~o.o DO 98 1-1,IZ 98 X-X+CHLA2(3,I)*V(I) X-1-X/CHLASUM ENOIF ZP-ZP-XKMZ+ZP*REPRD*X IF(ZP.LT.ZPMIN) ZP"ZPMIN C*********DETERMINATION OF GRAZING ON CHLA *********** ZPTDEL-ZP*(1-TD*.04167)*V(IZ)*.001 DO 100 n"1, IZ I-IZ+1-II DELZ=CHLATOT(I)*V(I)/CHLASUM DO 101 K-1,NCLASS DCHLA-CHLA2(K,I)-CHLAMIN(K) IF(DCHLA.LT.O.O) DCHLA-O.O GRAZ=GRAZMAX(K)*DCHLA/(HSCGRAZ(K)+DCHLA) + *DELZ*ZPTDEL/V(I)*THGRAZ(K)**T20(I) GRAZ-AMIN1(GRAZ,DCHLA) IF(GRAZ.LT.O.O) GRAZ=O.O GRTOT-GRTOT+GRAZ*V(I) 101 GRAZE(K,I)-GRAZ C**********COMPUTE CHANGE IN ZOOPL BIOMASS************ C**********AS DECREASE IN CONC. IN LAYER IZ*********** X-PREDMIN*(ZP*V(IZ)/V(I)-ZPMIN)*(1.-TD*0.04167)*DELZ IF(X.LT.O.O) X-O.O XKMZ=XKMZ+X XK"'XK+X*V(I) ZP-ZP-X*V(I)/V(IZ) 100 CONTINUE C********** CONVERT TO N OF ZOOPL. *************************** C********** CONC. COULD CAUSE PROBLEMS IN MERGE AND SPLIT ***** ZP-ZP*V(IZ) TYPE REAL RETURN END TOTAL CHLORO-A ABOVE THE DAY DEPTH (KG) PAGE 40 REAL REAL EXCESS CHLORO-A OVER MINIMUM REQUIRED FOR GRAZING TO OCCUR TIME DEPENDANT GRAZING COEFFICIENT CHLASU DCHLA DDAY DELZ 01 REAL RATIO OF CHLORO-A IN A LAYER (KG) TO TOTAL ABOVE THE DAYDEPTH LIGHT LIMITED PREDATION COEFFICIENT DY GRf4Z GRTOT I II K M60T PREDMA TO X XK ZPTDEL REAL REAL REAL REAL INTEGER*2 INTEGER*2 INTEGER*2 INTEGER*2 REAL REAL REAL REAL REAL \JULIAN DAY COMPUTED GRAZING RATE (MG CHLORO-A/L) TOTAL GRAZING (KG CHLORO-A) INDEX OF LAYERS DUMMY INDEX OF LAYERS FOR LAYERS FROM THE DAYDEPTH INDEX FOR ALGAL CLASS NUMBER OF LAYERS COMPUTED MAXIMUM PREDATION RATE (1/DAY) PHOTOPERIOD (HOURS) COMPUTED PREDATION DURING MIGRATION (N/M3/LAYER) TOTAL PREDATION IN A DAY (N/DAY) DUMMY VARIABLE USED IN GRAZING COMPUTATIONS PASS ONE NO ERRORS DETECTED 1165 SOURCE LINES 125 APPENDIX F DEFINITIONS OF VARIABLES IN COMMON BLOCKS 127 Name COMMON A /VOL ! ALPHA !CHANEL! ASM !ZOOPL ! ATOP /VOL ! AVGI !SUB ! BETA !WATER ! BOD2 !RESULT! BODIN !INFLOW! BODK20 !YIELD ! BRNH !YIELD ! BRNO !YIELD ! BRR !YIELD ! BW !CHANEL! C2 !RESULT! CD2 !RESULT! CDIN !INFLOW! CHLA2 !RESULT! CHLADU !TEMP ! CHLAIN !INFLOW! CHLAMI !ZOOPL ! CHLATO !RESULT! CHLMAX !PHYTO ! CHLMEA !PHYTO ! CIN !INFLOW! CR !MTHD ! DBL /VOL ! DEPTH !FIELD ! DIN !FILE ! DOIN !INFLOW! DRCT !MTHD ! DS02 !RESULT! DY !STEPS ! DZ !VOL ! DZLL !STEPS ! DZUL !STEPS ! ELCB !CHANEL! EMISS !WATER ! FDATA !FIELD ! FLO !FILE ! FVCHLA !FLOW ! GRAZE !ZOOPL ! GROMAX !PHYTO ! HKMAX !WATER ! HMK !FLOW ! HSC1 !PHYTO / HSC2 !PHYlO ! HSCGRA /ZOOPL ! HSCN !PHYTO / HSCNH !PHYTO ! HSCPA !PHYTO ! IDIATO /PHVTO / IFLAG /FIELD / ILAY !STEPS / IPRNT !CHOICE/ IREC /FILE / IZ /ZOOPL / LAY /SUB / MAXDAV /ZOOPL / MBOT /STEPS / MDAY /STEPS / MDAVRE /STAT / MDAVRM /STAT ! MET /FILE / MINDAV /ZOOPL / MODEL /CHDICE/ MTHREL /STAT / MTHRMS /STAT / NCLASS /CHOICE/ NDAVS /CHOICE/ COMMON BLOCK VARIABLES Discr1ption Area through the midp01nt of a layer (m2) Side slope of outflow channel (degrees from horzontal) Coefficient for assimilation of nutrients by zooplankton Surface layer of each layer (m2) Average PhAR in each layer (uE!m2!s) Surface radiation absorption coefficient Detritus as BOD (mg!l) Inflow BOD concentration (mg!l) Detrital decay rate (1!day) Benthic ammonium exchange rate (gm!m2/day) Benthic nitrate-nitrite exchange rate (Qm!m2!day) Benthic phosphorus exchange rate (gm!m2/day) Bottom width of outflow channel (m) Suspended solids concentration (mg!l) Dissolved solids concentration (mg!l) Inflow dissolved solids concentration (mg/l) Chlorophyll-a concentration (mg!l) Chloro-a (mg!l) from previous time step or itertion Inflow chlorophyll-a concentration (mg!l) Minimum chloro-a concentration for grazing (mg/l) Sum of chloro-a in all algal groups (mg!l) Year to date maximum chI oro-a concentration (mg!l) Two meter mean chloro-a concentration (mg!l) Inflow suspended solids concentration (mg!l) Cloudiness ratio from meteorological data file Elevation of lake bottom above some datum (m) Depth of field data sample (m) Name of input data file Inflow dissolved oxygen concentration (mg!l) Wind direction from meteorological data file Dissolved oxygen concentration (mg!l) Counter in Julian days from first day of simulation Layer thickness (m) Minimum allowable layer thickness (m) Maximum allowable layer thickness (m) Elevation of outflow channel bottom (m) Emmissivity of water Values from field data in input data file Name of 1nflow data file Fall velocity of algal cells (m/day) Graz1ng rate of zooplankton on algae (mg chloro-a!1nd.!day) Maximum rate of chI oro-a production by algae (1!day) Maximum hypo11mnetic diffus10n rate (m2/day) Harmonic mean diffusion rate between layers (m2!day) Optimum light level for growth (uE/m2/s) Light 1nhibited growth level (uE/m2!s) Half-saturation chloro-a concentration for grazing (mg/l) Half-saturation concentrat10n for total n1trogen uptake (mg/l) Half-saturatiorl concentration for ammon1um uptake (mg/l) Half-saturation concentration for phosphorus uptake (mg/l) Flag for number of silica 11mited algal groups Counter for number of data p01nts used 1n the regression analysis Index for layer at the mixed layer depth . Flags for d·iffere.nt model operation and output options Record number for direct access plot file Index for layer at the zooplankton day depth Number of subdivisions in each layer from SUBLAV Julian day for maximum seasonal . zooplankton predation Number of layers Current day of the month in the simulat10n Day of maximum relative deviation between model and field data Oay of maximum seasonal deviation between field data and model Meteorological data file name Julian day for minimum seasonal zooplankton predation Selected biological-nutrient modeling option (1,2,or 3) Month of maximum relative deviation between model and field data Month of maximum seasonal deviat10n between field data and model Number of algal classes treated in the model Number of days selected for profile or field data output ]28 .. j. NFLD NITRO NM NPRINT NPRNT PA2 PADUM PAIN PARID PC2 PCDUM PDEL PE PLT PMAX PMIN PR PRED~III PRMAX PRMIN PROD PRODSU PTSUM QE QIN RAD REL RELM REPRO RM RMS RS RSQ SB20 SO SDZ SECCHI 512 SID SUMX SUMlEF , WSTR,RHOW , HCAP, HF , TOW (XHoI)N~/PDEL(3),PHAX(3),PMIN(3),THR(3),THM(3),XKRI(3), +XKR2(3),XKM(3),HSCPA(3),HSCI(3),HSC2(3),UPMAX(3),THUP(3), +GRa1AX(3) ,'lMAX(3) ,TOP1'(3) ,XNMAX(3) ,XNMIN(3) ,UNMAX(3) ,TRUN(3), +HSCN ( 3 ) , HSCNH ( 3 ) ,XNDEL ( 3) , USCI , IDIA'ffM, CHUmAN, CHl.MAX, SECCHI (XHoI)N/yIELD/YCA, YCH02, Y2CH02, YCOOD, YPBOD, vzw , YPlP, YNZP, YZOO, + YSCIffi, YNHBOD, BRNO, BRNH ,XKNNH, '1lINNH, YPCIJLA, BODK20, SB20 , BRR CXM-IlN/ZOOPL/IZ, MINDAY ,MAXDAY , ZP, ZlMIN , PBMIN , PBMAX, PREJXt1IN ,XIMIN, +XIMAX,XKRZP, GRAZMAX ( 3 ) , THGRAZ ( 3 ) ,ASH, THRZP,:HSCXlRAZ ( 3) , CHLAMIN ( 3 ) , +REPRO, XI, XKMZ c:t:J.+Dl/SOLV/ AK(40) ,BK(40) ,CK(40) ,DK(40) cx:MI)N/MmD/TAIR(31) ,TDEW(31) ,RAD(31) ,CR(31) ,WIND(31), + PR(31),DRCT(31) OOMMON/RESULT/ T2(40),CHLATOT(40),PA2(40),PTSUM(40),BOD2(40), +DS02(40),C2(40),CD2(40),XN02(40),XNH2(40),CHLA2(3,40), +PC2(3,40),XNC2(3,40),T20(40),SI2(40) OOMHON/SOURCE{RM(3,40),PROD(40),XMR(3,40),P.RODSUM(40) OOMMON/FLOW/HMK(41),QE(40),FVCHLA(6),PB(5,41) OOMHON/TEMPIPARIO(4),PCDUM(3,40),XNHD(40),XNOD(40), + CHLADUM(3,40),XNCD(3,40),PADUM(40),SID(40) cnHlN/VOL/ZMAX,DZ(40) ,Z(40) ,A(40) ,V(40) ,TV(40) ,ATOP(41) ,DBL OOMMON/SUB/SDZ(60),SZ(60),LAY(40),AVGI(4,60),SVOL(60) OOMMON/CHOlCE/MODEL,NITRO,IPRNT(6),NDAYS,NPRNT(30),NCLASS,PLOT(10) OOMMON/CHANEL/WCHANL,EI£B ,ALPHA,BW, WLAKE OOMMON/STEPS/DZLL, DZUL,MBOT , NM,NPRINT , MDAY , r«N.l'H, ILAY, DY OOMMON/CONVAL/GRA,RIDA,OONST OOHMON/INFLOW/QIN(6),TIN(5),PAIN(5),BODIN(5),DOIN(6),CIN(5), +CDIN(5),XNHIN(5),XNOIN(6),CHLAIN(3,5) OOMMON/STAT/SUMXY(10),SUHX(10),SUHY(10),XSQ(10), +YSQ( 10) ,RSQ( 10) ,RBL( 10) ,RBUf( 10) ,Ml'HREL( 10) ,MDAYREL( 10) , ZREL( 10) , +ZREll1( 10) ,RB(10) ,Nt-5(10) ,Kl'HRMS(10) ,MDAYBMS( 10) , ZRS (10) ,ZRMS(10) GOTO (100,200,300,400,500,600,700,800,900,1000,1100,1200,1300) ID C********** AREA OOHPUTATION SECTION *********** 100 DUM=38481.*ZD+185946 IF(ZD.LT.l.524) ~160493.*ZD RBTURN C********** FETCH OOHPUTATION SECTION ****** 200 X=SIN«DRCT(MDAY)*10-67.6)/57.3) ZD=661.2*X*X+382.2 RBTURN C********** VOLUME OOHPUTATION SECTION ***** 300 DUM=119943.*ZD**1.4819 RBTURN CU*** CXHUl'E DEP1'H FlOt VOLUME ***** 400 ZD=(DUM*8. 33729B-6) **.67481 RE'IURN c***** TREATMENT SECTION ***** 500 RBTURN C***** PHOSPHORUS SOURCES/SINKS ***** 600 RBTURN 139 c***** N02-N03 SOURCES/SINKS ***** 700 RETURN C***** NH4 SOURCES/SINKS ***** 800 RETURN C***** 02 SOURCES/SINKS ***** 900 RETURN c***** 0UTFI..aV CG1PUTATION ***** 1000 DO 1001 I=l,NFLOW IF(QIN(I).GT.O.O) THEN TIN ( I ) =T2 (1 ) DOIN(I)=14.652-TIN(I)*(.41022-TIN(I)* + «7.991E-3)-7.7774E-6*TIN(I») ELSE C***** SET STAGE IN FEET AOOVE GAGE DATUM *** STAGE=ZMAX+DBL-263.41 IF(STAGE.GT.O.O) THEN C***** cx::MPUTE Q our USING STAGE-DISCHARGE BY ERG *** c***** ASSUME .667 OF 0l.1l'FIDl FRQt1 S. BABIN AND C***** OONVERT :li'Ra1 CFS TO M**3/DAY *** QIN(I)=-.18714E7 * STAGE**2.661 ELSE QIN(I)=O.O ENDIF ENDIF 1001 OONTINUE C************ COMPUTE OVERLAND INFLOW *********************** IF(PR(MDAY) .GT.O.O) THBN Q=ABS(PR(MDAY)-0.34)*(PR(MDAY)-0.34)/(PR(MDAY)+1.36) IF(Q.GT.O.O) THEN NFLOW=l +NFLOW C************OONVERT INCHES OF RUNOFF TO VOLUHB (M**3)******* QIN(NFLCW)=Q*38341 TIN(NFLCW)=TAIR(MDAY) CIN(NFLOW)=157. COIN (NFLCW) =0.0 DO 1002 I=l,NCLASS 1002 CHLAIN(I,NFLOW)=O.O PAIN (NFLOW) =0.1 BODIN(NFLOW)=29.8 DOIN(NFLOW)=14.652-TAIR(MOAY)*(.4102-TAIR(MDAY) + *(.OO7991-7.7774B-5*TAIR(MDAY») XNHIN(NFLOW)=.28 XNOIN(NFLOW)=.28 ENDIF ENDlF RETURN C***** GROUNDWATER FIm'