forked from CPRA-MP/ICM_Hydro
-
Notifications
You must be signed in to change notification settings - Fork 0
/
main.f
1291 lines (1159 loc) · 63.7 KB
/
main.f
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
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
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
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
180
181
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
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
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
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
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
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
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
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
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
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
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
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
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
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
1000
! This section contains the Doxygen comments for
!> @mainpage
!> @section Introduction
!> This is the Fortran code for the Coastal Master Plan Ecohydrology routines,
!> originally prepared by Dr. J. Alex McCorquodale at the University of New Orleans.
!> This code has been edited and updated by Eric White at The Water Institute of the Gulf.
!> @section Compiler
!> This source code was originally tested and compiled using Intel's Fortran compiler.
!BUG! EVENTUALLY ADD THIS PORTION TO THE DOXYGEN COMMENTS
! @section Participating Organizations
! @image html CPRAlogo.png
! <a href="http:\\www.coastal.LA.gov">Coastal Protection and Restoration Authority</a>
! @image html WIlogo.png
! @image html UNOlogo.png
!> @file
!> @brief This is the main program/control subroutine of the Ecohydrology model.
!> @details This is the control portion of the Ecohydrology model which controls the model I/O,
!> determines the simulation timesteps, assigns initial conditions for each timestep,
!> steps through each simulation timestep, and eventually post-processes model output to be read into
!> other ICM routines.
! Ecohydro program is set up solely to call the 'main' subroutine. The only reason 'main' is a subroutine the program is to incorporate the
! Doxygen documentation logic structure, which needs to be embedded within a subroutine, not a program.
program Ecohydro
!ZW 3/13/2015 move the log file open statement to here
!open (unit=1, file='hydro_run.log', status='unknown')
call main
write(1,*) '----------------------------------------------------'
write(1,*)
write(1,*) '----------------------------------------------------'
write(1,*) ' HYDRO MODEL RUN IS COMPLETE.'
write(1,*) '----------------------------------------------------'
write(*,*) '----------------------------------------------------'
write(*,*)
write(*,*) '----------------------------------------------------'
write(*,*) ' HYDRO MODEL RUN IS COMPLETE.'
write(*,*) '----------------------------------------------------'
! pause
stop
end
!> @author J. Alex McCorquodale - University of New Orleans
!> @author Eric White - The Water Institute of the Gulf
ccccccccccccccccccLower Breton/ Lower Barataria 10/31/13
c
c Main MASS BALANCE MODEL FOR Lake Pontchartrain WATER QUALITY
c
c NOTE: All distance units input as km are converted to m in code
c
c Elements
c i - M Internal links
c j - N Storage cells
c ib - Mb Boundary Links
c +ve direction into the cell
c
c NOTATION: follows EPA SWMM 3 and 4
c* Link no., u/s Cell no., d/s Cell No., Length, Width, Deptho, n, Km, Ken, Kx, Kb
c* Cell no., As, Eso, Bed, Ahydro,dAdz, CSS, S
c* Boundary link no., type(+1 u/s,-1 d/s), Cell no.,itype, Length, Width,depth, n, Km, Ken, Kx, Kb
c icc(j,10) link numbers connected to cell (- in, + out)
c* Cell no., link in -, link out + !up to 6 numbers fo rthe links connected to a cell
c* !connectivity
c* u/s BC Daily flows --> File=TimeSeries1.dat FORMAT free --> Day,flow m3/s, SS mg/L, S ppt
c* d/s BC Stage --> File=TimeSeries2.dat FORMAT free --> hours,stage m, SS mg/L, S ppt
c***********************************************************************************************
c* Matrix Variables:
c
c Cells;- j
!> @param[out] As(N,3) surface area of cells (m2)
!> @param[out] Eso(N) initial stage of storage cells (m)
!> @param[out] Bed(N) bed elevation of storage cells (m)
!> @param[out] ds(N,3) depth in storage cells (m)
!> @param[out] Es(N,3) stage in storage cells (m)
!> @param[out] Ahydro(N) area of hydrologically connected area for water balance (m2)
!> @param[out] ow % open water in attach drainage basin
!> @param[out] dAdz(N) change in As with depth (m)
!> @param[out] CSS(N,3) sediment concentration (mg/L)
!> @param[out] S(N,3) salinity (TDS) (ppt)
!> @param[out] acss(N) resuspension parameter
!> @param[out] bcss(N) resuspension parameter
!> @param[out] Vss(N) deposition velocity (m/d)
!> @param[out] tau(N,3) bed shear stress (Pa)
!> @param[out] Tauc(N) critical shear stress (Pa)
!> @param[out] Fetch(N,10) Fetch EW (km)
!> @param[out] Fetch(N,10) Fetch NS (km)
c Cellsh - Properties for marsh compontent of cells
!> @param[out] Ahf(N) Marsh area subject to flooding
!> @param[out] flood(N) Fraction of Ahydro subject to flooding
!> @param[out] Eho(N) Bed marsh elevation in flood prone area (m)
!> @param[out] Eh(N,3) stage in flooded area of marsh.
!> @param[out] CSSh(N,3) TSS in flood area of marsh
!> @param[out] Sh(N,3) Salinity in flood area of marsh (ppt)
!> @param[out] acssh(N) resupension in flooded area of marsh
!> @param[out] Vsh(N) settling velocity in flooded area of marsh (m/d)
!> @param[out] Esho(N) dead storage in marshes (m)
c
cc***********************************************************************************************
c Links:- i
!> @param[out] Deptho(M) initial channel depth (m)
!> @param[out] Depth(M,3) channel depth (m)
!> @param[out] Length(M) channel length (m)
!> @param[out] Width(M) channel width (m)
!> @param[out] an(M) channel Mannings n
!> @param[out] Km(M) minor loss coefficient due to structures based on Km(Vch)^2/2g
!> @param[out] Ken(100) entrance loss coefficient
!> @param[out] Kx(M) exit loss coefficient based on Kx(Vch-Vrec)^2/2g
!> @param[out] Exy(M) diffusion coefficient thru link between cells (m2/s)
!> @param[out] Q(M,3) channel flow in m3/s
!> @param[out] Qss(M,3) channel sediment flow (kg/s)
!> @param[out] QSal(M,3) channel TDS transport (kg/s)
!> @param[out] jus(M) u/s junction no.
!> @param[out] jds(M) d/s junction no.
!> @param[out] itype(M) "-1" is u/s BC, "+1" is d/s BC, "2" is two cell link.
c***********************************************************************************************
c Forcing:-
!> @param[out] Qtrib(26,366*25) Mean daily tributary flow (m3/s)
!> @param[out] CssT(26,366*25) Mean daily tributary sediment conc (mg/L)
!> @param[out] Strib(30,366*25) Mean daily tributary salinity (ppt)
!> @param[out] Tcoef(M,25) Tide generation coefficients
!> @param[out] Sal(1,1) d/s BC on Salinity (ppt)
!> @param[out] BCTSS(10) d/s BC on TSS (mg/L) !added TSS=CSS+Algae+Detritus XX!!NOT USED JAM July 09
!> @param[out] PET(366) daily potential evapotranspiration (m/d)
!> @param[out] Precip(366) daily rainfall (m/d)
!> @param[out] QCHEM(N,15,366*25) chemical load in tributaries (g/sec)
!> @param[out] CHEM(N,15,2) water quality concentration
!> @param[out] QChemSUM(15)
!> @param[out] DChem(15)
!> @param[out] CSSo(N)
!> @param[out] QMult(N,40) Multiplier on tributary-compartment lookup matrix
!> @param[out] PET(days,ETgages)
!> @param[out] Rain(days,raingages)
!> @param[out] Chem(N,15,2)
!> @param[out] Qmultdiv(N,40) Multiplier on diversion-compartment lookup matrix
!> @param[out] QChemdiv(26,15,366*25) chemical load in diversions (g/sec)
!> @param[out] cChemdiv(1,15,366*25) WQ concentrations in diversion input data (mg/L)
!> @param[out] QAtm(1,15,366*25)
!> @param[out] TSS(N)
!> @param[out] cChemface(15)
!> @param[out] T7(366*25)
!> @param[out] TKN(M,2)
!> @param[out] uplandNP(M,14)
c
c***********************************************************************************************
c Control parameters:-
!> @param[out] N number of compartments
!> @param[out] M number of links
!> @param[out] Mus number of u/s BCs
!> @param[out] Mds number of d/s BCs
!> @param[out] dt computational time step (s)
!> @param[out] dTprint output time step (hours)
!> @param[out] startrun Julian day of start of run
!> @param[out] endrun Julian day at end of run
!> @param[out] Specg Specific gravity of sediment
!> @param[out] time elapsed model time seconds
!> @param[out] day elapsed model time in decimals days
!> @param[out] thour elapsed model time in decimal hours
!> @param[out] tmon elapsed model time in decimal months
!> @param[out] kthr integer of hourly elapsed model time +1
!> @param[out] kmon integer of monthly elapsed model time +1
!> @param[out] kday integer of daily elapsed model time +1
!> @param[out] dday decimal portion of day, dday=0.0 at 0:00 (midnight)
!> @param[out] hday decimal portion of day normalized to noon, hday = 0.0 at 12:00 (noon)
!> @param ddayhr hour of day (dday/24.)
!> @param windupdate flag used to update wind data
!> @param dtwind timestep of observed wind data file (hr)
!> @param windrow row of wind data array that corresponds to hour of simulation
c
c***********************************************************************************************
c Chemical Inputs
c 1 = NO3 + NO2
c 2 = NH4
c 3 = DIN (1+2)
c 4 = Organic N
c 5 = TIP
c 6 = TOC
c 7 = DO
c 8 = Live Algae (ALG)
c 9 = Dead Algae (DET)
c 10= DON
c 11= DOP
c 12= DIP !PArtition SRP = ParP*TP
c 13= ChLa !PArtition Chla= ParCla*LivA
c 14= POP !-EDW used to say 14=TKN !-EDW
c***********************************************************************************************
! Main subroutine is the main 'control' subroutine for the ecohydro model.
! Embedded Doxygen comments only work within subroutines, not 'programs'.
subroutine main
use params
integer :: i,it,j,jj,jjj,sedclass !iterators used in main.f
Character*100 header
!>@par General Structure of Subroutine Logic:
! determine start time for calculating runtimes
call SYSTEM_CLOCK(runtime_start,count_rate1,count_max1)
!>> Open file with model information regarding file formatting in ICM (written only before first model year - NOT updated yearly).
open (unit=300, file= 'ICM_info_into_EH.txt', status ='unknown')
!>> Read in names of Veg files to write output to.
read(300,*) VegWaveAmpFile
read(300,*) VegMeanSalFile
read(300,*) VegSummerDepthFile
read(300,*) VegSummerSalFile
read(300,*) VegSummerTempFile
read(300,*) VegTreeEstCondFile
read(300,*) VegBIHeightFile
read(300,*) VegPerLandFile
!>> Read in grid dimensions for 500m grid used by LAVegMod.
read(300,*) n_500m_cells
read(300,*) veg_matrix_rows
read(300,*) veg_matrix_cols
read(300,*) veg_xllcorner
read(300,*) veg_yllcorner
!>> Read in grid dimensions for 1000m grid used by EwE.
read(300,*) n_1000m_cells
read(300,*) ewe_matrix_rows
read(300,*) ewe_matrix_cols
read(300,*) ewe_xllcorner
read(300,*) ewe_yllcorner
!>> Read in name of log file to write output to
read(300,*) Hydro_log
! Hydro_log = 'ICM_NoTimeStamp_Hydro.log'
close(300)
open (unit=1, file=TRIM(ADJUSTL(Hydro_log)),position='append')
write(1,*) '----------------------------------------------------'
write(1,*) '------------- RUNNING HYDROLOGY MODEL --------------'
write(1,*) '----------------------------------------------------'
write(1,*)
write(*,*) '----------------------------------------------------'
write(*,*) '------------- RUNNING HYDROLOGY MODEL --------------'
write(*,*) '----------------------------------------------------'
write(*,*)
!>> Read RunControlR.dat and set simulation variables (RunControlR.dat is updated yearly by ICM.py).
allocate(D50(4)) !sediment class parameters that need to be allocated for control data inputs (e.g. before allocate_params is called)
allocate(Tcrit(4)) !sediment class parameters that need to be allocated for control data inputs (e.g. before allocate_params is called)
open (unit=30, file= 'RuncontrolR.dat', status = 'unknown')
READ(30,*) year ! 1 year of simulation (e.g. 2015)
READ(30,*) Nyear ! 2 number of years simulated before exiting back to ICM
READ(30,*) dt ! 3 Computational time step
READ(30,*) M ! 4 # of links
READ(30,*) N ! 5 # of cells
READ(30,*) Mus ! 6 # of u/s BCs
READ(30,*) Mds ! 7 # of d/s BCs
READ(30,*) dTprint ! 8 Print time interval/ output time step (hr)
READ(30,*) startrun ! 9 Julian day start of run
READ(30,*) endrun ! 10 Julian day end of run
READ(30,*) Ntrib ! 11 Number of tributaries
READ(30,*) Ndiv ! 12 Number of diversions.
READ(30,*) Specg ! 13 Specific gravity of sediment
READ(30,*) Cp ! 14 Calibration for wind re-suspension
READ(30,*) dtwind ! 15 Time step for wind time series (hr)
READ(30,*) windgages ! 16 Number of windgages read in
READ(30,*) raingages ! 17 Number of precip gages read in
READ(30,*) etgages ! 18 Number of ET gages read in
READ(30,*) tidegages ! 19 Number of tide gages read in
! READ(30,*) Ndttrib ! Tributary Q & TSS u/s time interval (day)
READ(30,*) dttide ! 20 Time Step for tidal time series (hr)
READ(30,*) flocA ! 21 clay flocculation parameter (McAnally A)
READ(30,*) flocB ! 22 clay flocculation parameter (McAnally B)
READ(30,*) flocN ! 23 clay flocculation exponent (McAnally n)
READ(30,*) flocM ! 24 clay flocculation exponent (McAnally m)
READ(30,*) flocC1 ! 25 clay concentration where floculation settling begins (McAnally C1)(kg/m**3)
READ(30,*) flocC3 ! 26 clay concentration where floculation settling becominbegins to hinder settling (McAnally C3) (kg/m**3)
READ(30,*) Csalmax ! 27 upper limit to salinity impact on floc (salinity concentration - ppt)
READ(30,*) Pflocmax ! 28 upper limit on portion of fines available to floc (ratio - unitless)
READ(30,*) D50(1) ! 29 sand D50 (m) -- USDA particle size definition: 0.05 mm < sand < 2 mm
READ(30,*) D50(2) ! 30 silt D50 (m) -- USDA particle size definition: 0.002 mm < silt < 0.05 mm
READ(30,*) D50(3) ! 31 unflocculated clay D50 (m) -- USDA particle size definition: clay < 0.002 mm
READ(30,*) D50(4) ! 32 flocculated clay D50(m) -- USDA particle size definition: clay < 0.002 mm
READ(30,*) D90x ! 33 representative D90/D50 ratio
READ(30,*) Tcrit(1) ! 34 critical shear stress for sand particles (N/m**2)
READ(30,*) Tcrit(2) ! 35 critical shear stress for silt particlss (N/m**2)
READ(30,*) Tcrit(3) ! 36 critical shear stress for unflocculated clay particles (N/m**2)
READ(30,*) Tcrit(4) ! 37 critical shear stress for flocculated clay particles (N/m**2)
READ(30,*) saltox ! 38 Salinity concentration at which point algal growth is halved
READ(30,*) kphy20 ! 39 temperature-dependent phytoplankton mortality rate @ 20 degC (kphy20) (1/day)
READ(30,*) thetaphy ! 40 temperature-dependent phytoplankton mortality rate constant (thetaphy)
READ(30,*) kdet20 ! 41 temperature-dependent detritus hydrolysis rate @ 20 degC (kdet20) (1/day)
READ(30,*) thetadet ! 42 temperature-dependent detritus hydrolysis rate constant (thetaDET)
READ(30,*) kresp20 ! 43 temperature-dependent phytoplankton respiration rate @ 20 degC (kresp20) (1/day)
READ(30,*) thetaresp ! 44 temperature-dependent phytoplankton respirtion rate constant (thetaresp)
READ(30,*) kpo420 ! 45 temperature-dependent orthophosphate release rate @ 20 degC (kpo4) (1/day)
READ(30,*) thetapo4 ! 46 temperature-dependent orthophosphate rate coefficient (thetapo4)
READ(30,*) kdenit20 ! 47 temperature-dependent denitrification rate @ 20 degC (kdenit20) (1/day)
READ(30,*) thetadenit ! 48 temperature-dependent denitrification rate coefficient (thetadenit) (1/day)
READ(30,*) kphot20 ! 49 temperature-dependent photosynthetic rate @ 20 degC (kphot20)(1/day)
READ(30,*) thetaphot ! 50 temperature-dependent photosynthetic rate coefficient (thetaphot)(1/day)
READ(30,*) kdon20 ! 51 temperature-dependent dissolved organic N hydrolysis rate @ 20 degC (kdon20)(1/day)
READ(30,*) thetadon ! 52 temperature-dependent dissolved organic N hydrolysis rate coefficient (thetadon)(1/day)
READ(30,*) knit20num ! 53 depth-dependent nitrification rate numerator @ 20 degC (knit20num) (knit20 = knit20num/depth) (1/day)
READ(30,*) thetanit ! 54 temperature-dependent nitrification rate coefficient (thetanit) (1/day)
READ(30,*) knitmin ! 55 minimum threshold to apply to calculated nitrification rate @ 20-deg (knitmin) (1/day)
READ(30,*) knitmax ! 56 maximum threshold to apply to calculated nitrification rate @ 20-deg(knitmax) (1/day)
READ(30,*) KKexp ! 57 exponent on Kadlec & Knight depth term for flow into marsh (KKexp)
READ(30,*) defbedelev ! 58 default bed elevation of open water compartments - used to fill any missing bed elevation values (m NAVD88)
READ(30,*) defmarelev ! 59 default marsh elevation of open water compartments - used to fill any missing marsh elevation values (m NAVD88)
READ(30,*) def_n ! 60 default roughness value to use if input roughness for link is less than 0.0, or greater than 1.0
READ(30,*) maxconnect ! 61 maximum number of hydraulic connections per compartment (must be <= 20)
READ(30,*) stagemax ! 62 maximum water surface elevation value allowed to be passed to other ICM routines (m NAVD88) (999999999 if not to be used)
READ(30,*) salmax ! 63 maximum salinity concentration value allowed to be passed to other ICM routines (ppt) (999999999 if not to be used)
READ(30,*) tmpmax ! 64 maximum water temperature value allowed to be passed to other ICM routines (deg C) (999999999 if not to be used)
READ(30,*) sedmaxow ! 65 maximum annual open water sediment accumulation allowed to be passed to other ICM routines (kg/m2/yr) (999999999 if not to be used)
READ(30,*) sedmaxmi ! 66 maximum annual marsh interior sediment accumulation allowed to be passed to other ICM routines (kg/m2/yr) (999999999 if not to be used)
READ(30,*) sedmaxme ! 67 maximum annual marsh edge sediment accumulation allowed to be passed to other ICM routines (kg/m2/yr) (999999999 if not to be used)
READ(30,*) prismmax ! 68 maximum tidal prism volume allowed to be passed to other ICM routines (m3) (999999999 if not to be used)
READ(30,*) rangemax ! 69 maximum range in water surface elevation allowed to be passed to other ICM routines (m) (999999999 if not to be used)
READ(30,*) depthmax ! 70 maximum water depth value allowed to be passed to other ICM routines (m) (-9999 if not to be used)
READ(30,*) algmax ! 71 maximum algae concentration value allowed to be passed to other ICM routines (kg/m3) (999999999 if not to be used)
READ(30,*) tssmax ! 72 maximum TSS concentration value allowed to be passed to other ICM routines (kg/m3) (999999999 if not to be used)
READ(30,*) tknmax ! 73 maximum TKN concentration value allowed to be passed to other ICM routines (kg/m3) (999999999 if not to be used)
READ(30,*) nlinksw ! 74 number of links to write to FLO.out file (must match number of links in 'links_to_write.csv', set to 0 if no link flows are to be written)
READ(30,*) nstghr ! 75 number of compartments to write to STGhr.out file (must match number of compartments in 'hourly_stage_to_write.csv', set to 0 if no hourly stages are to be written)
READ(30,*) modeloverland ! 76 should overland flow be modeled in link types 8 and 9? (1=yes, 0=no)
READ(30,*) minwater ! 77 minimum water area required in each compartment (sq. meter)
READ(30,*) oscilflag ! 78 vertical oscillation in water level between timesteps that will flag 'dz' warning
READ(30,*) maxdz ! 79 maximum change in water level allowed between between timesteps - unrealistic but set to flag instabilities (maxdz) (meters)
READ(30,*) iTemp ! 80 modelling Temperature (No=0, Yes=1)
! READ(30,*) temp0 ! initial system water temp
READ(30,*) iSed ! 81 modelling Sediment (No=0, Yes=1)
READ(30,*) iSal ! 82 modelling Salinity (No=0, Yes=1)
READ(30,*) iWQ ! 83 modelling Water Quality chemical constituents (No=0, Yes=1)
READ(30,*) fpet ! 84 overwater evapotranspiration model (fpet) (1.0 = use PET input data 0.0 = Use average of PET input data)
READ(30,*) fa_def ! 85 upwind factor for weighting upwind/downwind contributions to WQ,salinity,TSS and temp transport calculations (fa_def) (range from 0.5 to 1) (0.5 gives equal weight to upwind and downwind, 1 will only use upwind)
READ(30,*) upwind_vel ! 86 channel velocity at which upwind factor, fa, is set to 1.0 (upwind_vel) (m/sec)
READ(30,*) fe ! 87 dispersion calibration factor (fe)
READ(30,*) lockdays ! 88 number of days in lock control observation timeseries data file LockControlObservedData.csv (lockdays) (days)
READ(30,*) dtlock ! 89 timestep of lock control observation timeseries data (dtlock) (hours)
READ(30,*) CSSmin ! 90 minimum CSS concentration allowed to be reached - calculated CSS below this value will be increased - this is applied to each sediment class separately (cssmin) (mg/L)
READ(30,*) CSSmax ! 91 maximum CSS concentration allowed to be reached - calculated CSS above this value will be throttled - this is applied to each sediment class separately (cssmax) (mg/L)
READ(30,*) TSSresusOff ! 92 TSS concentration threshold - used to determine maximum concentration of suspended sediment for each size class which deactivates bed resuspension (mg/L)
READ(30,*) sal_0_1_error ! 93 error term for freshest locations (0-1 ppt) - to be used to perturb output files (sal_0_1_error)
READ(30,*) sal_1_5_error ! 94 error term for intermediate salinity locations (1-5 ppt) - to be used to perturb output files (sal_1_5_error)
READ(30,*) sal_5_20_error ! 95 error term for saline locations (5-20 ppt) - to be used to perturb output files (sal_5_20_error)
READ(30,*) sal_20_35_error ! 96 error term for most saline locations (>20 ppt) - to be used to perturb output files (sal_20_35_error)
READ(30,*) tss_error ! 97 error term for TSS - to be used to perturb output files -percentage adjustment (between -1 and 1) (tss_error)
READ(30,*) stage_error ! 98 error term for stage - to be used to perturb output files (stage_error)
READ(30,*) stvar_error ! 99 error term for water level variability - to be used to perturb output files (stvar_error)
close(30)
! sal_0_1_error = 0.0
! sal_1_5_error = 0.0
! sal_5_20_error = 0.0
! sal_20_35_error = 0.0
! tss_error = 0.0
! stage_error = 0.0
! stvar_error = 0.0
!>> Determine number of simulation timesteps
simdays = endrun-startrun+1
NTs =int(simdays*24*60*60/dt)
NNN=NTs-int(24*60*60/dt)
lastdaystep = 24*60*60/dt !dt is in seconds - number of timesteps in a day
!>> Determine number of simulation steps needed before updating tide and wind data
lasttidestep = dttide*60*60/dt !dttide is in hours - number of timesteps before updating tide data
lastwindstep = dtwind*60*60/dt !dtwind is in hours - number of timesteps before updating wind data
lastlockstep = dtlock*60*60/dt !dtlock is in hours - number of timesteps before updating lock control data
!> Initialize counters for updating tide and wind data - initial values are set to 1 instead of zerob/c first day starts one timestep AFTER midnight - when midnight is hit at end of day this counter is then reset
daystep = 1
tidestep = 1
windstep = 1
lockstep = 1
!>> Call 'allocate_params' subroutine which allocates memory for almost all arrays used by this program
!>> (some temporary arrays are allocated in ICM_InterpolateToGrid subroutine which postprocesses model results).
call allocate_params
!>> Set initial conditions for constants and model parameters that are not included in input text files
! Initial Conditions
g=9.81 ! Gravity (m/s2)
pi=4.0*atan(1.0)
TemI = 15. ! Initial Water Temperature
KnN= 20. ! (ug/L) DIN Michaelis Constant Thomann & Mueller
KnP= 3. ! (ug/L) P Michaelis Constant
KnSS= 50. ! (mg/L) SS Michaelis Constant chged 30 to 50 JAM March 2011
KnSal=4. ! (ppt) Salinity Michaelis Constant
ParP= 0.4 ! SRP/TP in Tribs J. Day 1994 BCS trial
ParDOP= 0.1 ! DOP/TP in Tribs J. Day 1994 BCS trial
PARPOP=1.-ParDOP-ParP ! POP/TP in Tribs J. Day 1994 BCS trial
ParPMR= 0.2 ! SRP/TP in Diversions J. Day 1994 BCS trial
Pardpmr= 0.05
PPMR=1.-ParPMR-Pardpmr
ParSand=0.05 ! sand/TSS in Tribs and MR typical
ParCLa=0.03 ! Partition LivA --> ChlA
consd=24*3600 ! sec to days
conv=0.001*consd ! (mg/L)*(m3/s) --> kg/d
floodf(:)=0.0
nuo=0.000001 ! DEFAULT Viscosity
!>> Generate array for Day of Year value for the first day of each month
month_DOY(1) = 1
month_DOY(2) = 32
month_DOY(3) = 60
month_DOY(4) = 91
month_DOY(5) = 121
month_DOY(6) = 152
month_DOY(7) = 182
month_DOY(8) = 213
month_DOY(9) = 244
month_DOY(10) = 274
month_DOY(11) = 305
month_DOY(12) = 335
!>> Update first day of month for March through December during a leap year
if(simdays == 366) then
do j = 3,12
month_DOY(j) = month_DOY(j)+1
enddo
endif
C> initially set GrowAlgae array equal to zero
do j=1,N
do ichem=1,14
do me=1,14
GrowAlgae(j,ichem,me) = 0.
enddo
enddo
enddo
!>> determine upper CSS threshold for each sediment class where resuspension of bed material will be deactivated
!>> These values are based on an assumption that at a TSS value equal to the threshold concentration, the particle size distribution is 10% sand, 45% silt, and 45% clay.
CSSresusOff(1) = TSSresusOff*0.1
CSSresusOff(2) = TSSresusOff*0.45
CSSresusOff(3) = TSSresusOff*0.45
CSSresusOff(4) = TSSresusOff*0.45
!>> Open input text files.
open (unit=32, file= 'Cells.csv', status = 'unknown')
open (unit=323, file ='Fetch.csv', status = 'unknown')
open (unit=33, file= 'Links.csv', status = 'unknown') ! JAM Oct 2010
open (unit=34, file= 'LinksClosedHours.dat', status = 'unknown') !-EDW !value of 1 means link is closed for the hour, zero means it is open
open (unit=35,file='LockControlObservedData.csv',status='unknown')
! open (unit=36, file= 'SWR.dat', status = 'unknown')
open (unit=74, file= 'MissRToC.csv', status = 'unknown') ! JAM Oct 2010
open (unit=39, file= 'TribQ.csv', status = 'unknown') !Tributary flow (m3/s)
open (unit=40, file= 'PET.csv', status = 'unknown')
open (unit=42, file= 'Precip.csv', status='unknown')
open (unit=45, file= 'Meteorology.csv', status = 'unknown')
open (unit=44, file= 'AnthL.csv', status = 'unknown') ! Farm and Urban WW Loads (kg/d)
open (unit=43, file= 'WindVectorsX.csv',status= 'unknown')
open (unit=46, file= 'WindVectorsY.csv',status= 'unknown')
open (unit=47, file= 'TideData.csv',status='unknown')
open (unit=48, file= 'TideTranspose.csv',status='unknown')
open (unit=49, file= 'TideWeight.csv',status='unknown')
open (unit=77, file= 'QMult.csv', form= 'formatted',
& status = 'unknown')
open (unit=55, file= 'TribS.csv', status = 'unknown') ! Tributary sand concentration (mg/L)
open (unit=555,file= 'TribF.csv',status='unknown') ! Tributary fines concentration (mg/L)
open (unit=56, file= 'SBC.dat', status = 'unknown')
open (unit=80, file= 'NO2NO3Data.csv', status = 'unknown')
open (unit=81, file= 'NH4Data.csv', status = 'unknown')
open (unit=82, file= 'OrgNData.csv', status = 'unknown')
open (unit=83, file= 'PhosphorusData.csv', status ='unknown')
open (unit=84, file= 'AtmChemData.csv', status ='unknown')
open (unit=85, file= 'Decay.csv', status ='unknown')
open (unit=86, file= 'DivQm.csv', status ='unknown') ! diversion flow multiplier on Miss Riv flow
open (unit=87, file= 'DivWQ.csv', status ='unknown')
open (unit=88, file= 'QMult_div.csv', form= 'formatted',
& status ='unknown')
open (unit=89, file= 'DivSW.csv', status = 'unknown')
open (unit=101, file= 'BCToC2.dat', form = 'formatted')
open (unit=110, file= 'surge.csv', form = 'formatted')
! open (unit=117, file= 'AsedOW.csv',form ='formatted', ! Sediment Accretion !Status='unknown' added by Joao Pereira 5/17/2011
! & status ='unknown')
! open (unit=118, file= 'UplandNP.dat', form ='formatted')
open (unit=125, file= 'KBC.dat', status = 'unknown') !node numbers of open boundary
open (unit=126, file= 'links_to_write.csv',status='unknown') !input csv file with the link ID numbers of links to write flowrate to output file
open (unit=127, file='hourly_stage_to_write.csv',status='unknown') !input csv file with the ID numbers of compartments to write hourly stage to output file
!>> Open output text files (in append mode, if needed).
open (unit=70,file='DIN.out',form ='formatted',position='append')
open (unit=71,file='OrgN.out',form='formatted',position='append')
open (unit=72,file='TPH.out',form='formatted',position='append') ! TP.out
open (unit=73,file='TOC.out',form='formatted',position='append')
open (unit=75,file='SAL.out',form ='formatted',position='append') ! Salinity.out
open (unit=91,file='NO3.out',form='formatted',position='append') ! NO2NO3.out
open (unit=92,file='NH4.out',form = 'formatted',position='append')
open (unit=93,file='O2Sat.out',form='formatted',position='append')
open (unit=94,file='ALG.out',form='formatted',position='append')
open (unit=95,file='DO.out',form='formatted',position='append')
open (unit=96,file='TSS.out',form='formatted',position='append')
open (unit=97,file='DET.out',form='formatted',position='append') ! DeadAlgae.out
open (unit=100,file='TMP.out',form='formatted',position='append')
open (unit=103,file='SedAcc.out',form='formatted',position='append') ! Last row will be used to compute 20yr open water Acc.
open (unit=105,file='fflood.out',form='formatted',position='append')
open (unit=111,file='STG.out',form='formatted',position='append') ! ESAVE.OUT
open (unit=112,file='TRG.out',form='formatted',position='append') ! Range.out
open (unit=113,file='DON.out',form='formatted',position='append')
open (unit=119,file='SPH.out',form='formatted',position='append') ! SRP.out
open (unit=121,file='NRM.out',form='formatted',position='append') ! NRAcc.out -> Denitrification
open (unit=123,file='TKN.out',form='formatted',position='append')
open (unit=124,file='FLOm.out',form='formatted',position='append')
! read in information for grid cells used to pass data to other ICM routines !-EDW
open (unit=200, file='grid_lookup_500m.csv', form='formatted') ! compartment and link lookup table for 500-m grid cells
open (unit=201, file='grid_interp_dist_500m.csv',form='formatted') ! distance from each 500-m grid cell centroid to the compartment and link centroids
open (unit=202, file='grid_data_500m.csv', form='formatted') ! mean elevation for 500 m grid cells
open (unit=203, file='grid_IDs_Veg_matrix.csv', form='formatted') ! 500m grid cell names formatted in the matrix used by Vegetation ICM routine
open (unit=204, file='grid_500m_out.csv', form='formatted') ! output file for 500 m grid cells - in list form
open (unit=205, file='compartment_out.csv',form='formatted') ! output file for hydro compartments - summary values for ICM in list form
open(unit=206,file='sal_monthly_ave_500m.out',form='formatted')
open(unit=207,file='tmp_monthly_ave_500m.out',form='formatted')
open(unit=208,file='tkn_monthly_ave_500m.out',form='formatted')
open(unit=209,file='TSS_monthly_ave_500m.out',form='formatted')
open(unit=210,file='STGhr.out',form='formatted',position='append') !output file for hourly water level in Boundary Condition cells
open(unit=211,file='FLO.out',form='formatted',position='append') !output file for flowrate
open(unit=212,file='STGm.out',form='formatted',position='append')
! output files for use in the Vegetation ICM routine !-EDW
! these are written in append mode. ICM checks when first run as to whethere these files exist.
! If Ecohydro is run outside of the ICM these files may be erroneously appended to if they contain data and the model is re-run.
open (unit=301, file=TRIM(ADJUSTL(VegMeanSalFile)), ! mean salinity formatted for input into Vegetation ICM routine
& form='formatted', position='append')
open (unit=302, file=TRIM(ADJUSTL(VegSummerSalFile)), ! mean summertime salinity for input into Vegetation ICM routine
& form='formatted', position='append')
open (unit=303, file=TRIM(ADJUSTL(VegSummerDepthFile)), ! mean summertime water depth for input into Vegetation ICM routine
& form='formatted', position='append')
open (unit=304, file=TRIM(ADJUSTL(VegWaveAmpFile)), ! variance in daily water depth for input into Vegetation ICM routine
& form='formatted', position='append')
open (unit=305, file=TRIM(ADJUSTL(VegSummerTempFile)), ! mean summertime water temperature for input into Vegetation ICM routine
& form='formatted', position='append')
open (unit=306, file=TRIM(ADJUSTL(VegTreeEstCondFile)), ! tree establishment conditions for input into Vegetation ICM routine
& form='formatted', position='append')
open (unit=307, file=TRIM(ADJUSTL(VegBIHeightFile)), ! tree establishment conditions for input into Vegetation ICM routine
& form='formatted', position='append')
open (unit=308, file=TRIM(ADJUSTL(VegPerLandFile)), ! tree establishment conditions for input into Vegetation ICM routine
& form='formatted', position='append')
! hotstart files
open(unit=400,file='hotstart_in.dat',form='formatted')
open(unit=401,file='hotstart_out.dat',form='formatted')
!>> Read Boundary Conditions file
Read(125,*)(KBC(jj), jj=1,mds) !AMc Oct 8 2013
close(125)
! Initialize some variables and arrays
NR(:)=0.0
stds=0.
do j=1,N
accsed(j)=0.0
sal_ave(j)=0.0
cumul_retreat(j) = 0.0
enddo
do i=1,M
asedout(i)=0.0
enddo
!>> Call 'infile' subroutine, which reads input text files and stores data into arrays
call infile
!>> write header row for hourly stage output file (since not all compartments are printed)
if (nstghr > 0) then
write(210,908) 'Compartment:',(stghrwrite(jjk),jjk=1,nstghr)
else
write(210,*) 'Hourly water level not saved to file.'
endif
908 FORMAT(A,<nstghr-1>(I0,','),I0) ! first column has 'Compartment:##', followed by comma delimited list of boundary compartments
!>> write header row for flowrate output file (since not all links are printed)
! if (nlinksw > 0) then
! write(211,909) 'Link:',(linkswrite(jjk),jjk=1,nlinksw)
! else
! write(211,*) 'No links chosen to have flowrate outputs saved.'
! endif
!909 format(A,<nlinksw-1>(I0,','),I0) ! first column has 'Link:##', followed by comma delimited list of links
!>> Close input files that were imported in infile subroutine
close(32)
close(33)
close(34)
close(39)
close(40)
close(42)
close(44)
close(45)
close(43)
close(46)
close(47)
close(48)
close(49)
close(55)
close(56)
close(74)
close(77)
close(80)
close(81)
close(82)
close(83)
close(84)
close(85)
close(86)
close(87)
close(88)
close(89)
close(90)
close(101)
close(110)
! close(118)
close(200)
close(201)
close(202)
!>> Take first timestep of imported wind data and save into windx and windy arrays.
!>> These arrays will be overwritten at a delta t that matches the wind data timestep (this update occurs immediately prior to calling hydrod)
!>> 'windrow' is a counter that is incrementally updated each time a wind data timestep is reached
windrow = 1
do jjj = 1,N
windx(jjj) = windx_data(windrow,jwind(jjj))
windy(jjj) = windy_data(windrow,jwind(jjj))
enddo
!>> 'tiderow','surgerow', and 'lockrow' are counters that are incrementally updated each time a tide or lock control data timestep is reached
tiderow = 1
surgerow = 1
lockrow = 1
!>> Set initial conditions for links (from input files)
do i=1,M
fa(i) = fa_def*fa_mult(i) !Set array of initial upwind factor to default value
Q(i,1)=0.1
Q(i,2)=Q(i,1)
enddo
!>> Initialize hardcoded salinity and stage control trigger flags
SalLockStatusHNC = 1 !Open = 1 Close = -1
SalLockTriggerHNC = 1
StgTriggerSuperiorCanal = 1
StgTriggerStatusSuperiorCanal = 1
!SalLockStatusCSC = 1
!SalLockTriggerCSC = 1
!Atch_div_onoff = 1
!>> Read in hotstart file and set initial conditions (will overwrite some ICs set previously from input files)
write(1,*)
write(1,*)'-----------------------------------------------'
write(1,*)'Reading in hotstart file and setting
& values as initial conditons.'
write(1,*)'-----------------------------------------------'
write(*,*)
write(*,*)'-----------------------------------------------'
write(*,*)'Reading in hotstart file and setting
& values as initial conditons.'
write(*,*)'-----------------------------------------------'
read(400,*) ! ignore header row
do j=1,N
read(400,*) dump_int, !no need to save compartment number - read in to a dummy integer variable
& Es(j,1),
& S(j,1),
& Css(j,1,1),
& Css(j,1,2),
& Css(j,1,3),
& Css(j,1,4),
& Tempw(j,1),
& Chem(j,1,1),
& Chem(j,2,1),
& Chem(j,3,1),
& Chem(j,4,1),
& Chem(j,5,1),
& Chem(j,6,1),
& Chem(j,7,1),
& Chem(j,8,1),
& Chem(j,9,1),
& Chem(j,10,1),
& Chem(j,11,1),
& Chem(j,12,1),
& Chem(j,13,1),
& Chem(j,14,1) !Chem unit = mg/L
!>> Initialize some variables and arrays
Sandacc(j,1) = 0.0
Siltacc(j,1) = 0.0
Clayacc(j,1) = 0.0
As(j,2)= As(j,1) ! Surface area of cells (m2)
Es(j,2)= Es(j,1) ! Stage in storage cells (m)
ds(j,1)= Es(j,2)-Bed(j) ! Depth in storage cells (m)
Eh(j,2)=Eh(j,1) ! Stage in Marsh storage (m) !JAM Oct 2010
BCnosurge(j,1) = 0.0 ! Initialize no surge BC to zero for all compartments - only BC nodes will be updated - rest of array will be 0.0
BCnosurge(j,2) = BCnosurge(j,1) ! boundary conditions stage(m) before surge is added
Qmarsh(j,2) = Qmarsh(j,1) ! Flow in marsh cell !JAM Oct 2010
S(j,2) = S(j,1)
Tempw(j,2) = Tempw(j,1)
enddo
close(400)
Emax=0.0
Emin=0.0
!*****************************Start of the unsteady model run***********************************
!>> Start time stepping through model. MAIN LOOP THAT IS COMPLETED FOR EACH SIMULATION TIMESTEP.
!>> Main DO loop. Looped over each simulation timestep.
write(1,*)
write(1,*) '----------------------------------------------------'
write(1,*) 'START MAIN HYDRODYNAMIC MODEL'
write(1,*) '----------------------------------------------------'
write(*,*)
write(*,*) '----------------------------------------------------'
write(*,*) 'START MAIN HYDRODYNAMIC MODEL'
write(*,*) '----------------------------------------------------'
do mm= 1, NTs-1 ! ******do loop ends at ~ line 438
t=float(mm)*dt ! lapse time in seconds JAM 5/25/2011
!>> calculate various versions of time to be used as flags throughout program
day = t/3600./24.
dday=day-int(day) ! decimal portion of day, dday=0.0 at 0:00 (midnight), 0.99 at end of day
ddayhr = dday*24. ! decimal portion of day converted to hours dday = 0.00 @ midnight, 23.99 at end of day
!>> update counters for current timestep in day, tide data, and wind data
!>> counters are reset to zero at end of main.f timestepping Do loop when value meets the laststep values calculated at start of main.f
daystep = daystep + 1
tidestep = tidestep + 1
windstep = windstep + 1
lockstep = lockstep + 1
!>> -- Update counter that tracks with row of the tide data to use for the current model timestep
if (tidestep == 1) then
tiderow = tiderow + 1
surgerow = tiderow
endif
!>> -- Update wind arrays if timestep matches wind data timestep
! initial wind data for timestep = 0 was read into windx and windy arrays prior to timestepping began
! all subsequent timesteps that are divisible by dtwind will update the wind data arrays
if (windstep == 1) then
windrow = windrow + 1
do jjj = 1,N
windx(jjj) = windx_data(windrow,jwind(jjj))
windy(jjj) = windy_data(windrow,jwind(jjj))
enddo
endif
!>> -- Update lock control arrays if timesteps matches lock control data timestep
if (lockstep == 1) then
lockrow = lockrow + 1
endif
!>> -- Loop over links and number of diversions and calculate sediment accumulation timeseries for each model link.
kt=1+ifix(t/24./3600./365.25)
! do j=1,N
! do it=1,Ndiv
! ASANDD(j,kt)=ParSandD(kt)*Qdiv(it,kt)*cssTdiv(it,kt)
! & *Qmultdiv(j,it)
! enddo
! enddo
!>> -- Call 'hydrod' subroutine, which is another 'control' routine which calls all hydrodynamic subroutines at each simulation timestep.
call hydrod(mm)
!>> -- Loop over links. Save calculated flowrates, Q as the initial condition for the next simulation timestep.
do i=1,M
Q(i,1)=Q(i,2)
enddo
!>> -- Loop over compartments. Save numerous variables that were just calculated by 'hydrod' as the initial condition for the next simulation timestep.
do j=1,N
Es(j,1)=Es(j,2)
S(j,1)=S(j,2) ! resetting ICs
SL(j,1) = SL(j,2)
BCnosurge(j,1) = BCnosurge(j,2)
do sedclass=1,4
CSS(j,1,sedclass) = CSS(j,2,sedclass)
CSSh(j,1,sedclass) = CSSh(j,2,sedclass)
enddo
Tempw(j,1)=Tempw(j,2)
! Age(j,1)=Age(j,2)
Sacc(j,1)=Sacc(j,2)
Sacch_int(j,1)=Sacch_int(j,2)
Sacch_edge(j,1)=Sacch_edge(j,2)
Sandacc(j,1) = Sandacc(j,2)
Siltacc(j,1) = Siltacc(j,2)
Clayacc(j,1) = Clayacc(j,2)
do ichem = 1,14
Chem(j,ichem,1)=Chem(j,ichem,2)
enddo
CSSvRs(j,1)= CSSvRs(j,2)
Qmarsh(j,1) = Qmarsh(j,2) ! added EDW/ZW -02/16/2015
Eh(j,1) = Eh(j,2) ! added EDW/ZW - 02/16/2015
!>> Check if any water surface elevation values are NaN - if so, pause model run
if(isNAN(Es(j,2))) then
write(1,*)'Compartment',j,
& 'WSEL is NaN @ end of timestep=',mm
write(*,*)'Compartment',j,
& 'WSEL is NaN @ end of timestep=',mm
pause
endif
!!>> -- Check that some calculated water quality values do not exceed default threshold values. If they do, set equal to the threshold
! if(Chem(j,10,2) > 0.00025) then
! Chem(j,10,2) = 0.00025
! endif
!
! if(Chem(j,11,2) > 0.00025) then
! Chem(j,11,2) = 0.00025
! endif
!
! if(Chem(j,12,2) > 0.00025) then
! Chem(j,12,2) = 0.00025
! endif
!!>> -- Repeat water quality threshold check, so that the initial conditions array for the next model timestep matches the 'current' array.
! if(Chem(j,10,1) > 0.00025) then
! Chem(j,10,1)=0.00025
! endif
!
! if(Chem(j,11,1) > 0.00025) then
! Chem(j,11,1)=0.00025
! endif
!
! if(Chem(j,12,1) > 0.00025) then
! Chem(j,12,1)=0.00025
! endif
!zw 4/28/2015 remove low and high pass filters
enddo
!>> Check if any link flowrate values are NaN - if so, pause model run
do i=1,M
Q(i,1)=Q(i,2)
if(isNAN(Q(i,2))) then
write(1,*)'Link',i,'flow is NaN @ end of timestep=',mm
write(*,*)'Link',i,'flow is NaN @ end of timestep=',mm
write(*,*) ' Linkt=',linkt(i)
pause
endif
enddo
!>> Reset tidestep counter because end of observed tidal timestep is met
if (tidestep == lasttidestep) then
tidestep = 0
endif
!>> Reset windstep counter because end of observed tidal timestep is met
if (windstep == lastwindstep) then
windstep = 0
endif
!>> Reset lockstep counter because end of observed lock control timestep i smet
if (lockstep == lastlockstep) then
lockstep = 0
endif
!>> Reset daystep counter because end of day is met
if (daystep == lastdaystep) then
write(1,3333)'Year',year,'Day',int(day),'complete.'
write(*,3333)'Year',year,'Day',int(day),'complete.'
daystep = 0
endif
!>> End main model DO loop that is looped over each simulation timestep
enddo
!>> Add error terms to last timestep value before writing hotstart file
write(1,*) 'Adding error adjustment to last timestep values.'
write(*,*) 'Adding error adjustment to last timestep values.'
do j=1,N
Es(j,2) = Es(j,2) + stage_error
if(S(j,2) < 1.0) then
S(j,2) = max(0.0,S(j,2) + sal_0_1_error)
elseif(S(j,2) < 5.0) then
S(j,2) = max(0.0,S(j,2) + sal_1_5_error)
elseif(S(j,2) < 20.0) then
S(j,2) = max(0.0,S(j,2) + sal_5_20_error)
else
S(j,2) = max(0.0,S(j,2) + sal_20_35_error)
endif
Css(j,2,1) = max(0.0,Css(j,2,1)*(1.0+tss_error))
Css(j,2,2) = max(0.0,Css(j,2,2)*(1.0+tss_error))
Css(j,2,3) = max(0.0,Css(j,2,3)*(1.0+tss_error))
Css(j,2,4) = max(0.0,Css(j,2,4)*(1.0+tss_error))
enddo
!>> Save output values of last simulation timestep in a hotstart file
write(1,*)'Saving values from last timestep to a hotstart file.'
write(*,*)'Saving values from last timestep to a hotstart file.'
write(401,11140)'Compartment',
& 'Stage',
& 'Salinity',
& 'CSS_sand',
& 'CSS_silt',
& 'CSS_clay',
& 'CSS_clayfloc',
& 'WaterTemp',
& 'NO3_NO2',
& 'NH4',
& 'DIN',
& 'OrgN',
& 'TIP',
& 'TOC',
& 'DO',
& 'LiveAlg',
& 'DeadAlg',
& 'DON',
& 'DOP',
& 'DIP',
& 'ChlA',
& 'TKN'
do j=1,N
write(401,11142) j,
& Es(j,2),
& S(j,2),
& Css(j,2,1),
& Css(j,2,2),
& Css(j,2,3),
& Css(j,2,4),
& Tempw(j,2),
& min(Chem(j,1,2),1000.0),
& min(Chem(j,2,2),1000.0),
& min(Chem(j,3,2),1000.0),
& min(Chem(j,4,2),1000.0),
& min(Chem(j,5,2),1000.0),
& min(Chem(j,6,2),1000.0),
& min(Chem(j,7,2),1000.0),
& min(Chem(j,8,2),1000.0),
& min(Chem(j,9,2),1000.0),
& min(Chem(j,10,2),1000.0),
& min(Chem(j,11,2),1000.0),
& min(Chem(j,12,2),1000.0),
& min(Chem(j,13,2),1000.0),
& min(Chem(j,14,2),1000.0) ! chem unit = mg/L
enddo
close(401)
! Write sediment accumulation in open water output file - 1 value per year
! pm1 = (1.-Apctmarsh(j))
! WRITE(117,9229)(((ASandA(j)/As(j,1)+clams
! & +max(0.,Sacc(j,2)))*pm1 !modified June 20, 2011 JAM for testing removed /1000 added ASANDA(
! & +Sacch(j,2)*Apctmarsh(j)+476.), j=1,N) !increase fines captured
! close (117)
!>> Calculate percentage of sand in open water bed sediments
do j = 1,N
pct_sand_bed(j) = max(0.0,min(100.0,Sandacc(j,2)/
& (Sandacc(j,2)+Siltacc(j,2)+Clayacc(j,2))))
enddo
!>> Call 'ICM_Summaries' subroutine which calculates the average output values used by the other ICM routines (Wetland Morph, Vegetation, and Ecosytems)
call ICM_Summaries
!>> Call 'ICM_SalinityThreshold_WM' subroutine which calculates the maximum two-week salinity average for each hydro compartment (used in Marsh Collapse calculations in Wetland Morph)
call ICM_SalinityThreshold_WM
!>> Call 'ICM_MapToGrid' subroutine which maps (without interpolation) the hydro compartment results to the gridded structrure used by the other ICM routines.
call ICM_MapToGrid(1,500) ! Map annual mean salinity to grid
call ICM_MapToGrid(2,500) ! Map summer salinity to grid
call ICM_MapToGrid(3,500) ! Interpolate annual mean temperature to grid
call ICM_MapToGrid(4,500) ! Interpolate summer mean temperature to grid
call ICM_MapToGrid(5,500) ! Map annual mean water stage to grid
call ICM_MapToGrid(6,500) ! Map summer water stage to grid
call ICM_MapToGrid(7,500) ! Map summer water stage variance to grid
call ICM_MapToGrid(8,500) ! Map summer maximum 2-wk mean salinity to grid