112
113
118
119
120
121 integer, intent(in) :: ng, tile
122 integer, intent(in) :: LBi, UBi, LBj, UBj, UBk, UBt
123 integer, intent(in) :: IminS, ImaxS, JminS, JmaxS
124 integer, intent(in) :: nstp, nnew
125
126#ifdef ASSUMED_SHAPE
127# ifdef MASKING
128 real(r8), intent(in) :: rmask(LBi:,LBj:)
129# endif
130# if defined IRON_LIMIT && defined IRON_RELAX
131 real(r8), intent(in) :: h(LBi:,LBj:)
132# endif
133 real(r8), intent(in) :: Hz(LBi:,LBj:,:)
134 real(r8), intent(in) :: z_r(LBi:,LBj:,:)
135 real(r8), intent(in) :: z_w(LBi:,LBj:,0:)
136 real(r8), intent(in) :: srflx(LBi:,LBj:)
137 real(r8), intent(in) :: t(LBi:,LBj:,:,:,:)
138
139 real(r8), intent(in) :: tl_Hz(LBi:,LBj:,:)
140 real(r8), intent(in) :: tl_z_r(LBi:,LBj:,:)
141 real(r8), intent(in) :: tl_z_w(LBi:,LBj:,0:)
142 real(r8), intent(in) :: tl_srflx(LBi:,LBj:)
143 real(r8), intent(inout) :: tl_t(LBi:,LBj:,:,:,:)
144#else
145# ifdef MASKING
146 real(r8), intent(in) :: rmask(LBi:UBi,LBj:UBj)
147# endif
148# if defined IRON_LIMIT && defined IRON_RELAX
149 real(r8), intent(in) :: h(LBi:UBi,LBj:UBj)
150# endif
151 real(r8), intent(in) :: Hz(LBi:UBi,LBj:UBj,UBk)
152 real(r8), intent(in) :: z_r(LBi:UBi,LBj:UBj,UBk)
153 real(r8), intent(in) :: z_w(LBi:UBi,LBj:UBj,0:UBk)
154 real(r8), intent(in) :: srflx(LBi:UBi,LBj:UBj)
155 real(r8), intent(in) :: t(LBi:UBi,LBj:UBj,UBk,3,UBt)
156
157 real(r8), intent(in) :: tl_Hz(LBi:UBi,LBj:UBj,UBk)
158 real(r8), intent(in) :: tl_z_r(LBi:UBi,LBj:UBj,UBk)
159 real(r8), intent(in) :: tl_z_w(LBi:UBi,LBj:UBj,0:UBk)
160 real(r8), intent(in) :: tl_srflx(LBi:UBi,LBj:UBj)
161 real(r8), intent(inout) :: tl_t(LBi:UBi,LBj:UBj,UBk,3,UBt)
162#endif
163
164
165
166 integer, parameter :: Nsink = 2
167
168 integer :: Iter, i, ibio, isink, itime, itrc, iTrcMax, j, k, ks
169 integer :: Iteradj
170
171 integer, dimension(Nsink) :: idsink
172
173 real(r8), parameter :: MinVal = 1.0e-6_r8
174
175 real(r8) :: Att, ExpAtt, Itop, PAR
176 real(r8) :: tl_Att, tl_ExpAtt, tl_Itop, tl_PAR
177 real(r8) :: cff, cff1, cff2, cff3, cff4, cff5, cff6, dtdays
178 real(r8) :: tl_cff, tl_cff1, tl_cff4, tl_cff5, tl_cff6
179 real(r8) :: cffL, cffR, cu, dltL, dltR
180 real(r8) :: tl_cffL, tl_cffR, tl_cu, tl_dltL, tl_dltR
181 real(r8) :: fac, fac1, fac2
182 real(r8) :: tl_fac, tl_fac1, tl_fac2
183#ifdef IRON_LIMIT
184 real(r8) :: Nlimit, FNlim
185 real(r8) :: tl_Nlimit, tl_FNlim
186 real(r8) :: FNratio, FCratio, FCratioE, Flimit
187 real(r8) :: tl_FNratio, tl_FCratio, tl_FCratioE, tl_Flimit
188 real(r8) :: FeC2FeN, FeN2FeC
189# ifdef IRON_RELAX
190 real(r8) :: FeNudgCoef
191# endif
192#endif
193 real(r8), dimension(Nsink) :: Wbio
194 real(r8), dimension(Nsink) :: tl_Wbio
195
196 integer, dimension(IminS:ImaxS,N(ng)) :: ksource
197
198 real(r8), dimension(IminS:ImaxS) :: PARsur
199 real(r8), dimension(IminS:ImaxS) :: tl_PARsur
200
201 real(r8), dimension(NT(ng),2) :: BioTrc
202 real(r8), dimension(NT(ng),2) :: BioTrc1
203 real(r8), dimension(NT(ng),2) :: tl_BioTrc
204 real(r8), dimension(IminS:ImaxS,N(ng),NT(ng)) :: Bio
205 real(r8), dimension(IminS:ImaxS,N(ng),NT(ng)) :: Bio1
206 real(r8), dimension(IminS:ImaxS,N(ng),NT(ng)) :: Bio2
207 real(r8), dimension(IminS:ImaxS,N(ng),NT(ng)) :: Bio_old
208
209 real(r8), dimension(IminS:ImaxS,N(ng),NT(ng)) :: tl_Bio
210 real(r8), dimension(IminS:ImaxS,N(ng),NT(ng)) :: tl_Bio_old
211
212 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: FC
213 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: tl_FC
214
215 real(r8), dimension(IminS:ImaxS,N(ng)) :: Hz_inv
216 real(r8), dimension(IminS:ImaxS,N(ng)) :: Hz_inv2
217 real(r8), dimension(IminS:ImaxS,N(ng)) :: Hz_inv3
218 real(r8), dimension(IminS:ImaxS,N(ng)) :: Light
219 real(r8), dimension(IminS:ImaxS,N(ng)) :: WL
220 real(r8), dimension(IminS:ImaxS,N(ng)) :: WR
221 real(r8), dimension(IminS:ImaxS,N(ng)) :: bL
222 real(r8), dimension(IminS:ImaxS,N(ng)) :: bL1
223 real(r8), dimension(IminS:ImaxS,N(ng)) :: bR
224 real(r8), dimension(IminS:ImaxS,N(ng)) :: bR1
225 real(r8), dimension(IminS:ImaxS,N(ng)) :: qc
226
227 real(r8), dimension(IminS:ImaxS,N(ng)) :: tl_Hz_inv
228 real(r8), dimension(IminS:ImaxS,N(ng)) :: tl_Hz_inv2
229 real(r8), dimension(IminS:ImaxS,N(ng)) :: tl_Hz_inv3
230 real(r8), dimension(IminS:ImaxS,N(ng)) :: tl_Light
231 real(r8), dimension(IminS:ImaxS,N(ng)) :: tl_WL
232 real(r8), dimension(IminS:ImaxS,N(ng)) :: tl_WR
233 real(r8), dimension(IminS:ImaxS,N(ng)) :: tl_bL
234 real(r8), dimension(IminS:ImaxS,N(ng)) :: tl_bR
235 real(r8), dimension(IminS:ImaxS,N(ng)) :: tl_qc
236
237#include "set_bounds.h"
238
239
240
241
242
243
244
246
247
248
250
251#if defined IRON_LIMIT && defined IRON_RELAX
252
253
254
256#endif
257#ifdef IRON_LIMIT
258
259
260
261 fen2fec=(16.0_r8/106.0_r8)*1.0e3_r8
262 fec2fen=(106.0_r8/16.0_r8)*1.0e-3_r8
263#endif
264
265
266
269
270
271
272
275# ifdef TL_IOMS
278# else
281# endif
282
283 j_loop : DO j=jstr,jend
284
285
286
288 DO i=istr,iend
289 hz_inv(i,k)=1.0_r8/hz(i,j,k)
290 tl_hz_inv(i,k)=-hz_inv(i,k)*hz_inv(i,k)*tl_hz(i,j,k)+ &
291#ifdef TL_IOMS
292 & 2.0_r8*hz_inv(i,k)
293#endif
294 END DO
295 END DO
297 DO i=istr,iend
298 hz_inv2(i,k)=1.0_r8/(hz(i,j,k)+hz(i,j,k+1))
299 tl_hz_inv2(i,k)=-hz_inv2(i,k)*hz_inv2(i,k)* &
300 & (tl_hz(i,j,k)+tl_hz(i,j,k+1))+ &
301#ifdef TL_IOMS
302 & 2.0_r8*hz_inv2(i,k)
303#endif
304 END DO
305 END DO
307 DO i=istr,iend
308 hz_inv3(i,k)=1.0_r8/(hz(i,j,k-1)+hz(i,j,k)+hz(i,j,k+1))
309 tl_hz_inv3(i,k)=-hz_inv3(i,k)*hz_inv3(i,k)* &
310 & (tl_hz(i,j,k-1)+tl_hz(i,j,k)+ &
311 & tl_hz(i,j,k+1))+ &
312#ifdef TL_IOMS
313 & 2.0_r8*hz_inv3(i,k)
314#endif
315 END DO
316 END DO
317
318
319
323 DO i=istr,iend
324 bio(i,k,ibio)=0.0_r8
325 bio1(i,k,ibio)=0.0_r8
326 bio2(i,k,ibio)=0.0_r8
327 tl_bio(i,k,ibio)=0.0_r8
328 END DO
329 END DO
330 END DO
331
332
333
334
335
336
338 DO i=istr,iend
339
340
341
342
343
344
345
346
347
348
349
350
353
354
355 biotrc(ibio,nstp)=t(i,j,k,nstp,ibio)
356 tl_biotrc(ibio,nstp)=tl_t(i,j,k,nstp,ibio)
357
358
359 biotrc(ibio,nnew)=t(i,j,k,nnew,ibio)
360 tl_biotrc(ibio,nnew)=tl_t(i,j,k,nnew,ibio)* &
361 & hz_inv(i,k)+ &
362 & t(i,j,k,nnew,ibio)*hz(i,j,k)* &
363 & tl_hz_inv(i,k)- &
364# ifdef TL_IOMS
365 & biotrc(ibio,nnew)
366# endif
367 END DO
368
369
370
371 cff2=0.0_r8
372 DO itime=1,2
373 cff1=0.0_r8
374 tl_cff1=0.0_r8
376#ifdef IRON_LIMIT
378#else
380#endif
382 cff1=cff1+max(0.0_r8,minval-biotrc(ibio,itime))
383 tl_cff1=tl_cff1- &
384 & (0.5_r8-sign(0.5_r8, &
385 & biotrc(ibio,itime)-minval))* &
386 & tl_biotrc(ibio,itime)+ &
387# ifdef TL_IOMS
388 & (0.5_r8-sign(0.5_r8, &
389 & biotrc(ibio,itime)-minval))* &
390 & minval
391# endif
392 IF (biotrc(ibio,itime).gt.biotrc(itrcmax,itime)) THEN
393 itrcmax=ibio
394 END IF
395 biotrc1(ibio,itime)=biotrc(ibio,itime)
396 biotrc(ibio,itime)=max(minval,biotrc1(ibio,itime))
397 tl_biotrc(ibio,itime)=(0.5_r8- &
398 & sign(0.5_r8, &
399 & minval- &
400 & biotrc1(ibio,itime)))* &
401 & tl_biotrc(ibio,itime)+ &
402# ifdef TL_IOMS
403 & (0.5_r8+ &
404 & sign(0.5_r8, &
405 & minval- &
406 & biotrc1(ibio,itime)))* &
407 & minval
408# endif
409 END DO
410 IF (biotrc(itrcmax,itime).gt.cff1) THEN
411 biotrc(itrcmax,itime)=biotrc(itrcmax,itime)-cff1
412 tl_biotrc(itrcmax,itime)=tl_biotrc(itrcmax,itime)- &
413 & tl_cff1
414 END IF
415#ifdef IRON_LIMIT
418 biotrc1(ibio,itime)=biotrc(ibio,itime)
419 biotrc(ibio,itime)=max(minval,biotrc1(ibio,itime))
420 tl_biotrc(ibio,itime)=(0.5_r8- &
421 & sign(0.5_r8, &
422 & minval- &
423 & biotrc1(ibio,itime)))* &
424 & tl_biotrc(ibio,itime)+ &
425# ifdef TL_IOMS
426 & (0.5_r8+ &
427 & sign(0.5_r8, &
428 & minval- &
429 & biotrc1(ibio,itime)))* &
430 & minval
431# endif
432 END DO
433#endif
434 END DO
435
436
437
440 bio_old(i,k,ibio)=biotrc(ibio,nstp)
441 tl_bio_old(i,k,ibio)=tl_biotrc(ibio,nstp)
442 bio(i,k,ibio)=biotrc(ibio,nstp)
443 tl_bio(i,k,ibio)=tl_biotrc(ibio,nstp)
444 END DO
445
446#if defined IRON_LIMIT && defined IRON_RELAX
447
448
449
450
451
452 IF (h(i,j).le.
fehmin(ng))
THEN
453
454
455
457 & fenudgcoef*tl_bio(i,k,
ifdis)+ &
458# ifdef TL_IOMS
459 & fenudgcoef*
femax(ng)
460# endif
461 END IF
462#endif
463 END DO
464 END DO
465
466
467
468
469
470 DO i=istr,iend
471#ifdef CONST_PAR
472
473
474
475 parsur(i)=158.075_r8
476# ifdef TL_IOMS
477 tl_parsur(i)=0.0_r8
478# else
479 tl_parsur(i)=0.0_r8
480# endif
481#else
485# ifdef TL_IOMS
486 & parsur(i)
487# endif
488#endif
489 END DO
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 iter_loop:
DO iter=1,
bioiter(ng)
546
547
548
550 DO i=istr,iend
551
552
553
554
555
556
557
558
559
560
561
562
565
566
567 biotrc(ibio,nstp)=t(i,j,k,nstp,ibio)
568
569
570 biotrc(ibio,nnew)=t(i,j,k,nnew,ibio)
571 END DO
572
573
574
575 cff2=0.0_r8
576 DO itime=1,2
577 cff1=0.0_r8
579#ifdef IRON_LIMIT
581#else
583#endif
585 cff1=cff1+max(0.0_r8,minval-biotrc(ibio,itime))
586 IF (biotrc(ibio,itime).gt.biotrc(itrcmax,itime)) THEN
587 itrcmax=ibio
588 END IF
589 biotrc(ibio,itime)=max(minval,biotrc(ibio,itime))
590 END DO
591 IF (biotrc(itrcmax,itime).gt.cff1) THEN
592 biotrc(itrcmax,itime)=biotrc(itrcmax,itime)-cff1
593 END IF
594#ifdef IRON_LIMIT
597 biotrc1(ibio,itime)=biotrc(ibio,itime)
598 biotrc(ibio,itime)=max(minval,biotrc(ibio,itime))
599 END DO
600#endif
601 END DO
602
603
604
607 bio_old(i,k,ibio)=biotrc(ibio,nstp)
608 bio(i,k,ibio)=biotrc(ibio,nstp)
609 END DO
610
611#if defined IRON_LIMIT && defined IRON_RELAX
612
613
614
615
616
617 IF (h(i,j).le.
fehmin(ng))
THEN
620 END IF
621#endif
622 END DO
623 END DO
624
625
626
627
628
629 DO i=istr,iend
630#ifdef CONST_PAR
631
632
633
634 parsur(i)=158.075_r8
635#else
637#endif
638 END DO
639
640
641
642
643
644
645 DO iteradj=1,iter
646
647
648
649 DO i=istr,iend
650 par=parsur(i)
651 IF (parsur(i).gt.0.0_r8) THEN
653
654
655
656
657
659 & (z_w(i,j,k)-z_w(i,j,k-1))
660 expatt=exp(-att)
661 itop=par
662 par=itop*(1.0_r8-expatt)/att
663 light(i,k)=par
664
665
666
667
668 par=itop*expatt
669 END DO
670 ELSE
672 light(i,k)=0.0_r8
673 END DO
674 END IF
675 END DO
676
677
678
679
680
681
682#ifdef IRON_LIMIT
683
684
685
686
687
688
689
690
691
692
693#endif
694
699 DO i=istr,iend
700#ifdef IRON_LIMIT
701
702
703
704 fnratio=bio(i,k,
ifphy)/max(minval,bio(i,k,
iphyt))
705 fcratio=fnratio*fen2fec
707 flimit=fcratio*fcratio/ &
709
711 fnlim=min(1.0_r8,flimit/(bio(i,k,
ino3_)*nlimit))
712#endif
713 cff4=1.0_r8/sqrt(cff2+cff3*light(i,k)*light(i,k))
714 cff=bio(i,k,
iphyt)* &
715#ifdef IRON_LIMIT
716 & cff1*cff4*light(i,k)*fnlim*nlimit
717#else
718 & cff1*cff4*light(i,k)/ &
720#endif
726
727#ifdef IRON_LIMIT
728
729
730
731 fac=cff*bio(i,k,
ino3_)*fnratio/ &
732 & max(minval,bio(i,k,
ifdis))
740
741
742
743 cff5=dtdays*(fcratioe-fcratio)/
t_fe(ng)
744 cff6=bio(i,k,
iphyt)*cff5*fec2fen
745 IF (cff6.ge.0.0_r8) THEN
746 cff=cff6/max(minval,bio(i,k,
ifdis))
750 ELSE
751 cff=-cff6/max(minval,bio(i,k,
ifphy))
755 END IF
756#endif
757 END DO
758 END DO
759
760 IF (iteradj.ne.iter) THEN
761
762
763
764
765
766#ifdef IRON_LIMIT
767
768
769#endif
770
771 cff1=dtdays*
zoogr(ng)
774 DO i=istr,iend
775 cff=bio(i,k,
izoop)* &
776 & cff1*(1.0_r8-exp(-
ivlev(ng)*bio(i,k,
iphyt)))/ &
780 & bio(i,k,
iphyt)*cff2*cff
785#ifdef IRON_LIMIT
789#endif
790 END DO
791 END DO
792
793
794
795
798 cff1=1.0_r8/(1.0_r8+cff2+cff3)
800 DO i=istr,iend
803 & bio(i,k,
iphyt)*cff2
805 & bio(i,k,
iphyt)*cff3
806#ifdef IRON_LIMIT
810#endif
811 END DO
812 END DO
813
814
815
816
819 cff1=1.0_r8/(1.0_r8+cff2+cff3)
821 DO i=istr,iend
824 & bio(i,k,
izoop)*cff2
826 & bio(i,k,
izoop)*cff3
827 END DO
828 END DO
829
830
831
832 cff2=dtdays*
detrr(ng)
833 cff1=1.0_r8/(1.0_r8+cff2)
835 DO i=istr,iend
838 & bio(i,k,
isdet)*cff2
839 END DO
840 END DO
841
842
843
844
845
846
847
848
849
850 DO isink=1,nsink
851 ibio=idsink(isink)
852
853
854
855
856
857
859 DO i=istr,iend
860 qc(i,k)=bio(i,k,ibio)
861 END DO
862 END DO
863
865 DO i=istr,iend
866 fc(i,k)=(qc(i,k+1)-qc(i,k))*hz_inv2(i,k)
867 END DO
868 END DO
870 DO i=istr,iend
871 dltr=hz(i,j,k)*fc(i,k)
872 dltl=hz(i,j,k)*fc(i,k-1)
873 cff=hz(i,j,k-1)+2.0_r8*hz(i,j,k)+hz(i,j,k+1)
874 cffr=cff*fc(i,k)
875 cffl=cff*fc(i,k-1)
876
877
878
879
880 IF ((dltr*dltl).le.0.0_r8) THEN
881 dltr=0.0_r8
882 dltl=0.0_r8
883 ELSE IF (abs(dltr).gt.abs(cffl)) THEN
884 dltr=cffl
885 ELSE IF (abs(dltl).gt.abs(cffr)) THEN
886 dltl=cffr
887 END IF
888
889
890
891
892
893
894
895
896
897
898 cff=(dltr-dltl)*hz_inv3(i,k)
899 dltr=dltr-cff*hz(i,j,k+1)
900 dltl=dltl+cff*hz(i,j,k-1)
901 br(i,k)=qc(i,k)+dltr
902 bl(i,k)=qc(i,k)-dltl
903 wr(i,k)=(2.0_r8*dltr-dltl)**2
904 wl(i,k)=(dltr-2.0_r8*dltl)**2
905 END DO
906 END DO
907 cff=1.0e-14_r8
909 DO i=istr,iend
910 dltl=max(cff,wl(i,k ))
911 dltr=max(cff,wr(i,k+1))
912 br(i,k)=(dltr*br(i,k)+dltl*bl(i,k+1))/(dltr+dltl)
913 bl(i,k+1)=br(i,k)
914 END DO
915 END DO
916 DO i=istr,iend
918#if defined LINEAR_CONTINUATION
919 bl(i,
n(ng))=br(i,
n(ng)-1)
920 br(i,
n(ng))=2.0_r8*qc(i,
n(ng))-bl(i,
n(ng))
921#elif defined NEUMANN
922 bl(i,
n(ng))=br(i,
n(ng)-1)
923 br(i,
n(ng))=1.5*qc(i,
n(ng))-0.5_r8*bl(i,
n(ng))
924#else
925 br(i,
n(ng))=qc(i,
n(ng))
926 bl(i,
n(ng))=qc(i,
n(ng))
927 br(i,
n(ng)-1)=qc(i,
n(ng))
928#endif
929#if defined LINEAR_CONTINUATION
930 br(i,1)=bl(i,2)
931 bl(i,1)=2.0_r8*qc(i,1)-br(i,1)
932#elif defined NEUMANN
933 br(i,1)=bl(i,2)
934 bl(i,1)=1.5_r8*qc(i,1)-0.5_r8*br(i,1)
935#else
936 bl(i,2)=qc(i,1)
937 br(i,1)=qc(i,1)
938 bl(i,1)=qc(i,1)
939#endif
940 END DO
941
942
943
944
945
947 DO i=istr,iend
948 dltr=br(i,k)-qc(i,k)
949 dltl=qc(i,k)-bl(i,k)
950 cffr=2.0_r8*dltr
951 cffl=2.0_r8*dltl
952 IF ((dltr*dltl).lt.0.0_r8) THEN
953 dltr=0.0_r8
954 dltl=0.0_r8
955 ELSE IF (abs(dltr).gt.abs(cffl)) THEN
956 dltr=cffl
957 ELSE IF (abs(dltl).gt.abs(cffr)) THEN
958 dltl=cffr
959 END IF
960 br(i,k)=qc(i,k)+dltr
961 bl(i,k)=qc(i,k)-dltl
962 END DO
963 END DO
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979 cff=dtdays*abs(wbio(isink))
981 DO i=istr,iend
982 fc(i,k-1)=0.0_r8
983 wl(i,k)=z_w(i,j,k-1)+cff
984 wr(i,k)=hz(i,j,k)*qc(i,k)
985 ksource(i,k)=k
986 END DO
987 END DO
990 DO i=istr,iend
991 IF (wl(i,k).gt.z_w(i,j,ks)) THEN
992 ksource(i,k)=ks+1
993 fc(i,k-1)=fc(i,k-1)+wr(i,ks)
994 END IF
995 END DO
996 END DO
997 END DO
998
999
1000
1002 DO i=istr,iend
1003 ks=ksource(i,k)
1004 cu=min(1.0_r8,(wl(i,k)-z_w(i,j,ks-1))*hz_inv(i,ks))
1005 fc(i,k-1)=fc(i,k-1)+ &
1006 & hz(i,j,ks)*cu* &
1007 & (bl(i,ks)+ &
1008 & cu*(0.5_r8*(br(i,ks)-bl(i,ks))- &
1009 & (1.5_r8-cu)* &
1010 & (br(i,ks)+bl(i,ks)- &
1011 & 2.0_r8*qc(i,ks))))
1012 END DO
1013 END DO
1015 DO i=istr,iend
1016 bio(i,k,ibio)=qc(i,k)+ &
1017 & (fc(i,k)-fc(i,k-1))*hz_inv(i,k)
1018 END DO
1019 END DO
1020 END DO
1021 END IF
1022 END DO
1023
1024
1025
1026
1027
1028 DO i=istr,iend
1029 par=parsur(i)
1030# ifdef TL_IOMS
1031 tl_par=parsur(i)
1032# else
1033 tl_par=tl_parsur(i)
1034# endif
1035 IF (parsur(i).gt.0.0_r8) THEN
1037
1038
1039
1040
1041
1043 & (z_w(i,j,k)-z_w(i,j,k-1))
1045 & (z_w(i,j,k)-z_w(i,j,k-1))+ &
1047 & (tl_z_w(i,j,k)-tl_z_w(i,j,k-1))- &
1048# ifdef TL_IOMS
1050 & (z_w(i,j,k)-z_w(i,j,k-1))
1051# endif
1052 expatt=exp(-att)
1053 tl_expatt=-expatt*tl_att+ &
1054# ifdef TL_IOMS
1055 & (1.0_r8+att)*expatt
1056# endif
1057 itop=par
1058 tl_itop=tl_par
1059 par=itop*(1.0_r8-expatt)/att
1060 tl_par=(-tl_att*par+tl_itop*(1.0_r8-expatt)- &
1061 & itop*tl_expatt)/att+ &
1062# ifdef TL_IOMS
1063 & itop/att
1064# endif
1065
1066
1067 tl_light(i,k)=tl_par
1068
1069
1070
1071
1072 par=itop*expatt
1073 tl_par=tl_itop*expatt+itop*tl_expatt- &
1074# ifdef TL_IOMS
1075 & par
1076# endif
1077 END DO
1078 ELSE
1080
1081
1082 tl_light(i,k)=0.0_r8
1083 END DO
1084 END IF
1085 END DO
1086
1087
1088
1089
1090
1091
1092#ifdef IRON_LIMIT
1093
1094
1095
1096
1097
1098
1099
1100
1101
1102
1103#endif
1104
1109 DO i=istr,iend
1110#ifdef IRON_LIMIT
1111
1112
1113
1114
1115
1116 fac1=max(minval,bio1(i,k,
iphyt))
1117 tl_fac1=(0.5_r8-sign(0.5_r8,minval-bio1(i,k,
iphyt)))* &
1118 & tl_bio(i,k,
iphyt)+ &
1119# ifdef TL_IOMS
1120 & (0.5_r8+sign(0.5_r8,minval-bio1(i,k,
iphyt)))* &
1121 & minval
1122# endif
1123 fnratio=bio1(i,k,
ifphy)/fac1
1124 tl_fnratio=(tl_bio(i,k,
ifphy)-tl_fac1*fnratio)/fac1+ &
1125# ifdef TL_IOMS
1126 & fnratio
1127# endif
1128 fcratio=fnratio*fen2fec
1129 tl_fcratio=tl_fnratio*fen2fec
1133 & tl_bio(i,k,
ifdis)- &
1134# ifdef TL_IOMS
1135 & (
a_fe(ng)-1.0_r8)*fcratioe
1136# endif
1137 flimit=fcratio*fcratio/ &
1139 tl_flimit=2.0_r8*(tl_fcratio*fcratio- &
1140 & tl_fcratio*fcratio*flimit)/ &
1142# ifdef TL_IOMS
1143 & flimit*(fcratio*fcratio-
k_fec(ng)*
k_fec(ng))/ &
1145# endif
1146
1148 tl_nlimit=-tl_bio(i,k,
ino3_)*nlimit*nlimit+ &
1149# ifdef TL_IOMS
1150 & (
k_no3(ng)+2.0_r8*bio1(i,k,
ino3_))*nlimit*nlimit
1151# endif
1152
1153
1154 fac1=flimit/(bio1(i,k,
ino3_)*nlimit)
1155 tl_fac1=tl_flimit/(bio1(i,k,
ino3_)*nlimit)- &
1156 & (tl_bio(i,k,
ino3_)*nlimit+ &
1157 & bio1(i,k,
ino3_)*tl_nlimit)*fac1/ &
1158 & (bio1(i,k,
ino3_)*nlimit)+ &
1159# ifdef TL_IOMS
1160 & 2.0_r8*fac1
1161# endif
1162 fnlim=min(1.0_r8,fac1)
1163 tl_fnlim=(0.5_r8+sign(0.5_r8,1.0_r8-fac1))*tl_fac1+ &
1164# ifdef TL_IOMS
1165 & (0.5_r8-sign(0.5_r8,1.0_r8-fac1))
1166# endif
1167#endif
1168 cff4=1.0_r8/sqrt(cff2+cff3*light(i,k)*light(i,k))
1169 tl_cff4=-cff3*tl_light(i,k)*light(i,k)*cff4*cff4*cff4+ &
1170#ifdef TL_IOMS
1171 & (cff2+2.0_r8*cff3*light(i,k)*light(i,k))* &
1172 & cff4*cff4*cff4
1173#endif
1174 cff=bio1(i,k,
iphyt)* &
1175#ifdef IRON_LIMIT
1176 & cff1*cff4*light(i,k)*fnlim*nlimit
1177#else
1178 & cff1*cff4*light(i,k)/ &
1180#endif
1181#ifdef IRON_LIMIT
1182 tl_cff=tl_bio(i,k,
iphyt)* &
1183 & cff1*cff4*light(i,k)*fnlim*nlimit+ &
1184 & bio1(i,k,
iphyt)*cff1*cff4* &
1185 & (tl_light(i,k)*fnlim*nlimit+ &
1186 & light(i,k)*tl_fnlim*nlimit+ &
1187 & light(i,k)*fnlim*tl_nlimit)+ &
1188 & bio1(i,k,
iphyt)*cff1*tl_cff4* &
1189 & light(i,k)*fnlim*nlimit- &
1190# ifdef TL_IOMS
1191 & 4.0_r8*cff
1192# endif
1193#else
1194 tl_cff=(tl_bio(i,k,
iphyt)*cff1*cff4*light(i,k)+ &
1195 & bio1(i,k,
iphyt)*cff1* &
1196 & (tl_cff4*light(i,k)+cff4*tl_light(i,k))- &
1197 & tl_bio(i,k,
ino3_)*cff)/ &
1199# ifdef TL_IOMS
1202# endif
1203#endif
1204
1205
1207 & tl_cff*bio(i,k,
ino3_))/ &
1208 & (1.0_r8+cff)+ &
1209#ifdef TL_IOMS
1210 & cff*bio(i,k,
ino3_)/ &
1211 & (1.0_r8+cff)
1212#endif
1213
1214
1215
1217 & tl_bio(i,k,
ino3_)*cff+ &
1218 & bio(i,k,
ino3_)*tl_cff- &
1219#ifdef TL_IOMS
1220 & bio(i,k,
ino3_)*cff
1221#endif
1222#ifdef IRON_LIMIT
1223
1224
1225
1226
1227
1228 fac1=max(minval,bio1(i,k,
ifdis))
1229 tl_fac1=(0.5_r8-sign(0.5_r8,minval-bio1(i,k,
ifdis)))* &
1230 & tl_bio(i,k,
ifdis)+ &
1231# ifdef TL_IOMS
1232 & (0.5_r8+sign(0.5_r8,minval-bio1(i,k,
ifdis)))* &
1233 & minval
1234# endif
1235 fac2=1.0_r8/fac1
1236 tl_fac2=-fac2*fac2*tl_fac1+ &
1237# ifdef TL_IOMS
1238 & 2.0_r8*fac2
1239# endif
1240 fac=cff*bio(i,k,
ino3_)*fnratio*fac2
1241 tl_fac=fnratio*fac2*(tl_cff*bio(i,k,
ino3_)+ &
1242 & cff*tl_bio(i,k,
ino3_))+ &
1243 & cff*bio(i,k,
ino3_)*(tl_fnratio*fac2+ &
1244 & fnratio*tl_fac2)- &
1245# ifdef TL_IOMS
1246 & 3.0_r8*fac
1247# endif
1248
1249
1251 & tl_fac*bio2(i,k,
ifdis))/ &
1252 & (1.0_r8+fac)+ &
1253# ifdef TL_IOMS
1254 & fac*bio2(i,k,
ifdis)/(1.0_r8+fac)
1255# endif
1256
1257
1258
1260 & tl_bio(i,k,
ifdis)*fac+ &
1261 & bio2(i,k,
ifdis)*tl_fac- &
1262# ifdef TL_IOMS
1263 & bio2(i,k,
ifdis)*fac
1264# endif
1265
1266
1267
1268 cff5=dtdays*(fcratioe-fcratio)/
t_fe(ng)
1269 tl_cff5=dtdays*(tl_fcratioe-tl_fcratio)/
t_fe(ng)
1270 cff6=bio(i,k,
iphyt)*cff5*fec2fen
1271 tl_cff6=(tl_bio(i,k,
iphyt)*cff5+ &
1272 & bio(i,k,
iphyt)*tl_cff5)*fec2fen- &
1273# ifdef TL_IOMS
1274 & cff6
1275# endif
1276 IF (cff6.ge.0.0_r8) THEN
1277
1278
1279 fac1=max(minval,bio2(i,k,
ifdis))
1280 tl_fac1=(0.5_r8-sign(0.5_r8,minval-bio2(i,k,
ifdis)))* &
1281 & tl_bio(i,k,
ifdis)+ &
1282# ifdef TL_IOMS
1283 & (0.5_r8+sign(0.5_r8,minval-bio2(i,k,
ifdis)))* &
1284 & minval
1285# endif
1286 cff=cff6/fac1
1287 tl_cff=(tl_cff6-tl_fac1*cff)/fac1+ &
1288# ifdef TL_IOMS
1289 & cff
1290# endif
1291
1292
1294 & tl_cff*bio(i,k,
ifdis))/ &
1295 & (1.0_r8+cff)+ &
1296# ifdef TL_IOMS
1297 & cff*bio(i,k,
ifdis)/(1.0_r8+cff)
1298# endif
1299
1300
1301
1303 & tl_bio(i,k,
ifdis)*cff+ &
1304 & bio(i,k,
ifdis)*tl_cff- &
1305# ifdef TL_IOMS
1306 & bio(i,k,
ifdis)*cff
1307# endif
1308 ELSE
1309
1310
1311 fac1=-max(minval,bio2(i,k,
ifphy))
1312 tl_fac1=-(0.5_r8-sign(0.5_r8,minval-bio2(i,k,
ifphy)))* &
1313 & tl_bio(i,k,
ifphy)- &
1314# ifdef TL_IOMS
1315 & (0.5_r8+sign(0.5_r8,minval-bio2(i,k,
ifphy)))* &
1316 & minval
1317# endif
1318 cff=cff6/fac1
1319 tl_cff=(tl_cff6-tl_fac1*cff)/fac1+ &
1320# ifdef TL_IOMS
1321 & cff
1322# endif
1323
1324
1326 & tl_cff*bio(i,k,
ifphy))/ &
1327 & (1.0_r8+cff)+ &
1328# ifdef TL_IOMS
1329 & cff*bio(i,k,
ifphy)/(1.0_r8+cff)
1330# endif
1331
1332
1333
1335 & tl_bio(i,k,
ifphy)*cff+ &
1336 & bio(i,k,
ifphy)*tl_cff- &
1337# ifdef TL_IOMS
1338 & bio(i,k,
ifphy)*cff
1339# endif
1340 END IF
1341#endif
1342 END DO
1343 END DO
1344
1345
1346
1348 DO i=istr,iend
1349
1350
1351
1352
1353
1354
1355
1356
1357
1358
1359
1360
1363
1364
1365 biotrc(ibio,nstp)=t(i,j,k,nstp,ibio)
1366
1367
1368 biotrc(ibio,nnew)=t(i,j,k,nnew,ibio)
1369 END DO
1370
1371
1372
1373 cff2=0.0_r8
1374 DO itime=1,2
1375 cff1=0.0_r8
1377#ifdef IRON_LIMIT
1379#else
1381#endif
1383 cff1=cff1+max(0.0_r8,minval-biotrc(ibio,itime))
1384 IF (biotrc(ibio,itime).gt.biotrc(itrcmax,itime)) THEN
1385 itrcmax=ibio
1386 END IF
1387 biotrc(ibio,itime)=max(minval,biotrc(ibio,itime))
1388 END DO
1389 IF (biotrc(itrcmax,itime).gt.cff1) THEN
1390 biotrc(itrcmax,itime)=biotrc(itrcmax,itime)-cff1
1391 END IF
1392#ifdef IRON_LIMIT
1395 biotrc1(ibio,itime)=biotrc(ibio,itime)
1396 biotrc(ibio,itime)=max(minval,biotrc(ibio,itime))
1397 END DO
1398#endif
1399 END DO
1400
1401
1402
1405 bio_old(i,k,ibio)=biotrc(ibio,nstp)
1406 bio(i,k,ibio)=biotrc(ibio,nstp)
1407 END DO
1408
1409#if defined IRON_LIMIT && defined IRON_RELAX
1410
1411
1412
1413
1414
1415 IF (h(i,j).le.
fehmin(ng))
THEN
1418 END IF
1419#endif
1420 END DO
1421 END DO
1422
1423
1424
1425
1426
1427 DO i=istr,iend
1428#ifdef CONST_PAR
1429
1430
1431
1432 parsur(i)=158.075_r8
1433#else
1435#endif
1436 END DO
1437
1438
1439
1440
1441
1442
1443 DO iteradj=1,iter
1444
1445
1446
1447 DO i=istr,iend
1448 par=parsur(i)
1449 IF (parsur(i).gt.0.0_r8) THEN
1451
1452
1453
1454
1455
1457 & (z_w(i,j,k)-z_w(i,j,k-1))
1458 expatt=exp(-att)
1459 itop=par
1460 par=itop*(1.0_r8-expatt)/att
1461 light(i,k)=par
1462
1463
1464
1465
1466 par=itop*expatt
1467 END DO
1468 ELSE
1470 light(i,k)=0.0_r8
1471 END DO
1472 END IF
1473 END DO
1474
1475
1476
1477
1478
1479
1480#ifdef IRON_LIMIT
1481
1482
1483
1484
1485
1486
1487
1488
1489
1490
1491#endif
1492
1497 DO i=istr,iend
1498#ifdef IRON_LIMIT
1499
1500
1501
1502 fnratio=bio(i,k,
ifphy)/max(minval,bio(i,k,
iphyt))
1503 fcratio=fnratio*fen2fec
1505 flimit=fcratio*fcratio/ &
1507
1509 fnlim=min(1.0_r8,flimit/(bio(i,k,
ino3_)*nlimit))
1510#endif
1511 cff4=1.0_r8/sqrt(cff2+cff3*light(i,k)*light(i,k))
1512 cff=bio(i,k,
iphyt)* &
1513#ifdef IRON_LIMIT
1514 & cff1*cff4*light(i,k)*fnlim*nlimit
1515#else
1516 & cff1*cff4*light(i,k)/ &
1518#endif
1521 & bio(i,k,
ino3_)*cff
1522
1523#ifdef IRON_LIMIT
1524
1525
1526
1527 fac=cff*bio(i,k,
ino3_)*fnratio/ &
1528 & max(minval,bio(i,k,
ifdis))
1531 & bio(i,k,
ifdis)*fac
1532
1533
1534
1535 cff5=dtdays*(fcratioe-fcratio)/
t_fe(ng)
1536 cff6=bio(i,k,
iphyt)*cff5*fec2fen
1537 IF (cff6.ge.0.0_r8) THEN
1538 cff=cff6/max(minval,bio(i,k,
ifdis))
1541 & bio(i,k,
ifdis)*cff
1542 ELSE
1543 cff=-cff6/max(minval,bio(i,k,
ifphy))
1546 & bio(i,k,
ifphy)*cff
1547 END IF
1548#endif
1549 END DO
1550 END DO
1551
1552
1553
1554
1555
1556#ifdef IRON_LIMIT
1557
1558
1559#endif
1560
1561 cff1=dtdays*
zoogr(ng)
1564 DO i=istr,iend
1565 cff=bio(i,k,
izoop)* &
1566 & cff1*(1.0_r8-exp(-
ivlev(ng)*bio(i,k,
iphyt)))/ &
1572 & bio(i,k,
iphyt)*cff2*cff
1577#ifdef IRON_LIMIT
1583#endif
1584 END DO
1585 END DO
1586
1587
1588
1589
1592 cff1=1.0_r8/(1.0_r8+cff2+cff3)
1594 DO i=istr,iend
1597 & bio(i,k,
iphyt)*cff2
1599 & bio(i,k,
iphyt)*cff3
1600#ifdef IRON_LIMIT
1604#endif
1605 END DO
1606 END DO
1607
1608 IF (iteradj.ne.iter) THEN
1609
1610
1611
1612
1615 cff1=1.0_r8/(1.0_r8+cff2+cff3)
1617 DO i=istr,iend
1620 & bio(i,k,
izoop)*cff2
1622 & bio(i,k,
izoop)*cff3
1623 END DO
1624 END DO
1625
1626
1627
1628 cff2=dtdays*
detrr(ng)
1629 cff1=1.0_r8/(1.0_r8+cff2)
1631 DO i=istr,iend
1634 & bio(i,k,
isdet)*cff2
1635 END DO
1636 END DO
1637
1638
1639
1640
1641
1642
1643
1644
1645
1646 DO isink=1,nsink
1647 ibio=idsink(isink)
1648
1649
1650
1651
1652
1653
1655 DO i=istr,iend
1656 qc(i,k)=bio(i,k,ibio)
1657 END DO
1658 END DO
1659
1661 DO i=istr,iend
1662 fc(i,k)=(qc(i,k+1)-qc(i,k))*hz_inv2(i,k)
1663 END DO
1664 END DO
1666 DO i=istr,iend
1667 dltr=hz(i,j,k)*fc(i,k)
1668 dltl=hz(i,j,k)*fc(i,k-1)
1669 cff=hz(i,j,k-1)+2.0_r8*hz(i,j,k)+hz(i,j,k+1)
1670 cffr=cff*fc(i,k)
1671 cffl=cff*fc(i,k-1)
1672
1673
1674
1675
1676 IF ((dltr*dltl).le.0.0_r8) THEN
1677 dltr=0.0_r8
1678 dltl=0.0_r8
1679 ELSE IF (abs(dltr).gt.abs(cffl)) THEN
1680 dltr=cffl
1681 ELSE IF (abs(dltl).gt.abs(cffr)) THEN
1682 dltl=cffr
1683 END IF
1684
1685
1686
1687
1688
1689
1690
1691
1692
1693
1694 cff=(dltr-dltl)*hz_inv3(i,k)
1695 dltr=dltr-cff*hz(i,j,k+1)
1696 dltl=dltl+cff*hz(i,j,k-1)
1697 br(i,k)=qc(i,k)+dltr
1698 bl(i,k)=qc(i,k)-dltl
1699 wr(i,k)=(2.0_r8*dltr-dltl)**2
1700 wl(i,k)=(dltr-2.0_r8*dltl)**2
1701 END DO
1702 END DO
1703 cff=1.0e-14_r8
1705 DO i=istr,iend
1706 dltl=max(cff,wl(i,k ))
1707 dltr=max(cff,wr(i,k+1))
1708 br(i,k)=(dltr*br(i,k)+dltl*bl(i,k+1))/(dltr+dltl)
1709 bl(i,k+1)=br(i,k)
1710 END DO
1711 END DO
1712 DO i=istr,iend
1714#if defined LINEAR_CONTINUATION
1715 bl(i,
n(ng))=br(i,
n(ng)-1)
1716 br(i,
n(ng))=2.0_r8*qc(i,
n(ng))-bl(i,
n(ng))
1717#elif defined NEUMANN
1718 bl(i,
n(ng))=br(i,
n(ng)-1)
1719 br(i,
n(ng))=1.5*qc(i,
n(ng))-0.5_r8*bl(i,
n(ng))
1720#else
1721 br(i,
n(ng))=qc(i,
n(ng))
1722 bl(i,
n(ng))=qc(i,
n(ng))
1723 br(i,
n(ng)-1)=qc(i,
n(ng))
1724#endif
1725#if defined LINEAR_CONTINUATION
1726 br(i,1)=bl(i,2)
1727 bl(i,1)=2.0_r8*qc(i,1)-br(i,1)
1728#elif defined NEUMANN
1729 br(i,1)=bl(i,2)
1730 bl(i,1)=1.5_r8*qc(i,1)-0.5_r8*br(i,1)
1731#else
1732 bl(i,2)=qc(i,1)
1733 br(i,1)=qc(i,1)
1734 bl(i,1)=qc(i,1)
1735#endif
1736 END DO
1737
1738
1739
1740
1741
1743 DO i=istr,iend
1744 dltr=br(i,k)-qc(i,k)
1745 dltl=qc(i,k)-bl(i,k)
1746 cffr=2.0_r8*dltr
1747 cffl=2.0_r8*dltl
1748 IF ((dltr*dltl).lt.0.0_r8) THEN
1749 dltr=0.0_r8
1750 dltl=0.0_r8
1751 ELSE IF (abs(dltr).gt.abs(cffl)) THEN
1752 dltr=cffl
1753 ELSE IF (abs(dltl).gt.abs(cffr)) THEN
1754 dltl=cffr
1755 END IF
1756 br(i,k)=qc(i,k)+dltr
1757 bl(i,k)=qc(i,k)-dltl
1758 END DO
1759 END DO
1760
1761
1762
1763
1764
1765
1766
1767
1768
1769
1770
1771
1772
1773
1774
1775 cff=dtdays*abs(wbio(isink))
1777 DO i=istr,iend
1778 fc(i,k-1)=0.0_r8
1779 wl(i,k)=z_w(i,j,k-1)+cff
1780 wr(i,k)=hz(i,j,k)*qc(i,k)
1781 ksource(i,k)=k
1782 END DO
1783 END DO
1786 DO i=istr,iend
1787 IF (wl(i,k).gt.z_w(i,j,ks)) THEN
1788 ksource(i,k)=ks+1
1789 fc(i,k-1)=fc(i,k-1)+wr(i,ks)
1790 END IF
1791 END DO
1792 END DO
1793 END DO
1794
1795
1796
1798 DO i=istr,iend
1799 ks=ksource(i,k)
1800 cu=min(1.0_r8,(wl(i,k)-z_w(i,j,ks-1))*hz_inv(i,ks))
1801 fc(i,k-1)=fc(i,k-1)+ &
1802 & hz(i,j,ks)*cu* &
1803 & (bl(i,ks)+ &
1804 & cu*(0.5_r8*(br(i,ks)-bl(i,ks))- &
1805 & (1.5_r8-cu)* &
1806 & (br(i,ks)+bl(i,ks)- &
1807 & 2.0_r8*qc(i,ks))))
1808 END DO
1809 END DO
1811 DO i=istr,iend
1812 bio(i,k,ibio)=qc(i,k)+ &
1813 & (fc(i,k)-fc(i,k-1))*hz_inv(i,k)
1814 END DO
1815 END DO
1816 END DO
1817 END IF
1818 END DO
1819
1820
1821
1822
1823
1824
1825
1826#ifdef IRON_LIMIT
1827
1828
1829#endif
1830
1831 cff1=dtdays*
zoogr(ng)
1834 DO i=istr,iend
1835 cff=bio1(i,k,
izoop)* &
1836 & cff1*(1.0_r8-exp(-
ivlev(ng)*bio1(i,k,
iphyt)))/ &
1838 tl_cff=(tl_bio(i,k,
izoop)* &
1839 & cff1*(1.0_r8-exp(-
ivlev(ng)*bio1(i,k,
iphyt)))+ &
1842 & tl_bio(i,k,
iphyt)*cff)/ &
1843 & bio1(i,k,
iphyt)- &
1844#ifdef TL_IOMS
1845 & bio1(i,k,
izoop)* &
1848 & 1.0_r8)/ &
1850#endif
1851
1852
1854 & tl_cff*bio(i,k,
iphyt))/ &
1855 & (1.0_r8+cff)+ &
1856#ifdef TL_IOMS
1857 & cff*bio(i,k,
iphyt)/ &
1858 & (1.0_r8+cff)
1859#endif
1860
1861
1862
1864 & cff2*(tl_bio(i,k,
iphyt)*cff+ &
1865 & bio(i,k,
iphyt)*tl_cff)- &
1866#ifdef TL_IOMS
1867 & bio(i,k,
iphyt)*cff2*cff
1868#endif
1869
1870
1871
1874 & bio(i,k,
iphyt)*tl_cff)- &
1875#ifdef TL_IOMS
1877#endif
1878
1879
1880
1883 & bio(i,k,
iphyt)*tl_cff)- &
1884#ifdef TL_IOMS
1886#endif
1887
1888#ifdef IRON_LIMIT
1889
1890
1892 & tl_cff*bio2(i,k,
ifphy))/ &
1893 & (1.0_r8+cff)+ &
1894# ifdef TL_IOMS
1895 & cff*bio2(i,k,
ifphy)/(1.0_r8+cff)
1896# endif
1897
1898
1899
1901 & (tl_bio(i,k,
ifphy)*cff+ &
1903# ifdef TL_IOMS
1905# endif
1906#endif
1907 END DO
1908 END DO
1909
1910
1911
1912
1915 cff1=1.0_r8/(1.0_r8+cff2+cff3)
1917 DO i=istr,iend
1918
1919
1921
1922
1923
1925 & tl_bio(i,k,
iphyt)*cff2
1926
1927
1928
1930 & tl_bio(i,k,
iphyt)*cff3
1931
1932#ifdef IRON_LIMIT
1933
1934
1936
1937
1938
1940 & tl_bio(i,k,
ifphy)*(cff2+cff3)*
ferr(ng)
1941#endif
1942 END DO
1943 END DO
1944
1945
1946
1947
1950 cff1=1.0_r8/(1.0_r8+cff2+cff3)
1952 DO i=istr,iend
1953
1954
1956
1957
1958
1960 & tl_bio(i,k,
izoop)*cff2
1961
1962
1963
1965 & tl_bio(i,k,
izoop)*cff3
1966 END DO
1967 END DO
1968
1969
1970
1971 cff2=dtdays*
detrr(ng)
1972 cff1=1.0_r8/(1.0_r8+cff2)
1974 DO i=istr,iend
1975
1976
1978
1979
1980
1982 & tl_bio(i,k,
isdet)*cff2
1983 END DO
1984 END DO
1985
1986
1987
1989 DO i=istr,iend
1990
1991
1992
1993
1994
1995
1996
1997
1998
1999
2000
2001
2004
2005
2006 biotrc(ibio,nstp)=t(i,j,k,nstp,ibio)
2007
2008
2009 biotrc(ibio,nnew)=t(i,j,k,nnew,ibio)
2010 END DO
2011
2012
2013
2014 cff2=0.0_r8
2015 DO itime=1,2
2016 cff1=0.0_r8
2018#ifdef IRON_LIMIT
2020#else
2022#endif
2024 cff1=cff1+max(0.0_r8,minval-biotrc(ibio,itime))
2025 IF (biotrc(ibio,itime).gt.biotrc(itrcmax,itime)) THEN
2026 itrcmax=ibio
2027 END IF
2028 biotrc(ibio,itime)=max(minval,biotrc(ibio,itime))
2029 END DO
2030 IF (biotrc(itrcmax,itime).gt.cff1) THEN
2031 biotrc(itrcmax,itime)=biotrc(itrcmax,itime)-cff1
2032 END IF
2033#ifdef IRON_LIMIT
2036 biotrc1(ibio,itime)=biotrc(ibio,itime)
2037 biotrc(ibio,itime)=max(minval,biotrc(ibio,itime))
2038 END DO
2039#endif
2040 END DO
2041
2042
2043
2046 bio_old(i,k,ibio)=biotrc(ibio,nstp)
2047 bio(i,k,ibio)=biotrc(ibio,nstp)
2048 END DO
2049
2050#if defined IRON_LIMIT && defined IRON_RELAX
2051
2052
2053
2054
2055
2056 IF (h(i,j).le.
fehmin(ng))
THEN
2059 END IF
2060#endif
2061 END DO
2062 END DO
2063
2064
2065
2066
2067
2068 DO i=istr,iend
2069#ifdef CONST_PAR
2070
2071
2072
2073 parsur(i)=158.075_r8
2074#else
2076#endif
2077 END DO
2078
2079
2080
2081
2082
2083
2084 DO iteradj=1,iter
2085
2086
2087
2088 DO i=istr,iend
2089 par=parsur(i)
2090 IF (parsur(i).gt.0.0_r8) THEN
2092
2093
2094
2095
2096
2098 & (z_w(i,j,k)-z_w(i,j,k-1))
2099 expatt=exp(-att)
2100 itop=par
2101 par=itop*(1.0_r8-expatt)/att
2102 light(i,k)=par
2103
2104
2105
2106
2107 par=itop*expatt
2108 END DO
2109 ELSE
2111 light(i,k)=0.0_r8
2112 END DO
2113 END IF
2114 END DO
2115
2116
2117
2118
2119
2120
2121#ifdef IRON_LIMIT
2122
2123
2124
2125
2126
2127
2128
2129
2130
2131
2132#endif
2133
2138 DO i=istr,iend
2139#ifdef IRON_LIMIT
2140 fnratio=bio(i,k,
ifphy)/max(minval,bio(i,k,
iphyt))
2141 fcratio=fnratio*fen2fec
2143 flimit=fcratio*fcratio/ &
2145
2147 fnlim=min(1.0_r8,flimit/(bio(i,k,
ino3_)*nlimit))
2148#endif
2149 cff4=1.0_r8/sqrt(cff2+cff3*light(i,k)*light(i,k))
2150 cff=bio(i,k,
iphyt)* &
2151#ifdef IRON_LIMIT
2152 & cff1*cff4*light(i,k)*fnlim*nlimit
2153#else
2154 & cff1*cff4*light(i,k)/ &
2156#endif
2159 & bio(i,k,
ino3_)*cff
2160
2161#ifdef IRON_LIMIT
2162
2163
2164
2165 fac=cff*bio(i,k,
ino3_)*fnratio/ &
2166 & max(minval,bio(i,k,
ifdis))
2169 & bio(i,k,
ifdis)*fac
2170
2171
2172
2173 cff5=dtdays*(fcratioe-fcratio)/
t_fe(ng)
2174 cff6=bio(i,k,
iphyt)*cff5*fec2fen
2175 IF (cff6.ge.0.0_r8) THEN
2176 cff=cff6/max(minval,bio(i,k,
ifdis))
2179 & bio(i,k,
ifdis)*cff
2180 ELSE
2181 cff=-cff6/max(minval,bio(i,k,
ifphy))
2184 & bio(i,k,
ifphy)*cff
2185 END IF
2186#endif
2187 END DO
2188 END DO
2189
2190
2191
2192
2193
2194#ifdef IRON_LIMIT
2195
2196
2197#endif
2198
2199 cff1=dtdays*
zoogr(ng)
2202 DO i=istr,iend
2203 cff=bio(i,k,
izoop)* &
2204 & cff1*(1.0_r8-exp(-
ivlev(ng)*bio(i,k,
iphyt)))/ &
2208 & bio(i,k,
iphyt)*cff2*cff
2213#ifdef IRON_LIMIT
2217#endif
2218 END DO
2219 END DO
2220
2221
2222
2223
2226 cff1=1.0_r8/(1.0_r8+cff2+cff3)
2228 DO i=istr,iend
2231 & bio(i,k,
iphyt)*cff2
2233 & bio(i,k,
iphyt)*cff3
2234#ifdef IRON_LIMIT
2238#endif
2239 END DO
2240 END DO
2241
2242
2243
2244
2247 cff1=1.0_r8/(1.0_r8+cff2+cff3)
2249 DO i=istr,iend
2252 & bio(i,k,
izoop)*cff2
2254 & bio(i,k,
izoop)*cff3
2255 END DO
2256 END DO
2257
2258
2259
2260 cff2=dtdays*
detrr(ng)
2261 cff1=1.0_r8/(1.0_r8+cff2)
2263 DO i=istr,iend
2266 & bio(i,k,
isdet)*cff2
2267 END DO
2268 END DO
2269
2270 IF (iteradj.ne.iter) THEN
2271
2272
2273
2274
2275
2276
2277
2278
2279
2280 DO isink=1,nsink
2281 ibio=idsink(isink)
2282
2283
2284
2285
2286
2287
2289 DO i=istr,iend
2290 qc(i,k)=bio(i,k,ibio)
2291 END DO
2292 END DO
2293
2295 DO i=istr,iend
2296 fc(i,k)=(qc(i,k+1)-qc(i,k))*hz_inv2(i,k)
2297 END DO
2298 END DO
2300 DO i=istr,iend
2301 dltr=hz(i,j,k)*fc(i,k)
2302 dltl=hz(i,j,k)*fc(i,k-1)
2303 cff=hz(i,j,k-1)+2.0_r8*hz(i,j,k)+hz(i,j,k+1)
2304 cffr=cff*fc(i,k)
2305 cffl=cff*fc(i,k-1)
2306
2307
2308
2309
2310 IF ((dltr*dltl).le.0.0_r8) THEN
2311 dltr=0.0_r8
2312 dltl=0.0_r8
2313 ELSE IF (abs(dltr).gt.abs(cffl)) THEN
2314 dltr=cffl
2315 ELSE IF (abs(dltl).gt.abs(cffr)) THEN
2316 dltl=cffr
2317 END IF
2318
2319
2320
2321
2322
2323
2324
2325
2326
2327
2328 cff=(dltr-dltl)*hz_inv3(i,k)
2329 dltr=dltr-cff*hz(i,j,k+1)
2330 dltl=dltl+cff*hz(i,j,k-1)
2331 br(i,k)=qc(i,k)+dltr
2332 bl(i,k)=qc(i,k)-dltl
2333 wr(i,k)=(2.0_r8*dltr-dltl)**2
2334 wl(i,k)=(dltr-2.0_r8*dltl)**2
2335 END DO
2336 END DO
2337 cff=1.0e-14_r8
2339 DO i=istr,iend
2340 dltl=max(cff,wl(i,k ))
2341 dltr=max(cff,wr(i,k+1))
2342 br(i,k)=(dltr*br(i,k)+dltl*bl(i,k+1))/(dltr+dltl)
2343 bl(i,k+1)=br(i,k)
2344 END DO
2345 END DO
2346 DO i=istr,iend
2348#if defined LINEAR_CONTINUATION
2349 bl(i,
n(ng))=br(i,
n(ng)-1)
2350 br(i,
n(ng))=2.0_r8*qc(i,
n(ng))-bl(i,
n(ng))
2351#elif defined NEUMANN
2352 bl(i,
n(ng))=br(i,
n(ng)-1)
2353 br(i,
n(ng))=1.5*qc(i,
n(ng))-0.5_r8*bl(i,
n(ng))
2354#else
2355 br(i,
n(ng))=qc(i,
n(ng))
2356 bl(i,
n(ng))=qc(i,
n(ng))
2357 br(i,
n(ng)-1)=qc(i,
n(ng))
2358#endif
2359#if defined LINEAR_CONTINUATION
2360 br(i,1)=bl(i,2)
2361 bl(i,1)=2.0_r8*qc(i,1)-br(i,1)
2362#elif defined NEUMANN
2363 br(i,1)=bl(i,2)
2364 bl(i,1)=1.5_r8*qc(i,1)-0.5_r8*br(i,1)
2365#else
2366 bl(i,2)=qc(i,1)
2367 br(i,1)=qc(i,1)
2368 bl(i,1)=qc(i,1)
2369#endif
2370 END DO
2371
2372
2373
2374
2375
2377 DO i=istr,iend
2378 dltr=br(i,k)-qc(i,k)
2379 dltl=qc(i,k)-bl(i,k)
2380 cffr=2.0_r8*dltr
2381 cffl=2.0_r8*dltl
2382 IF ((dltr*dltl).lt.0.0_r8) THEN
2383 dltr=0.0_r8
2384 dltl=0.0_r8
2385 ELSE IF (abs(dltr).gt.abs(cffl)) THEN
2386 dltr=cffl
2387 ELSE IF (abs(dltl).gt.abs(cffr)) THEN
2388 dltl=cffr
2389 END IF
2390 br(i,k)=qc(i,k)+dltr
2391 bl(i,k)=qc(i,k)-dltl
2392 END DO
2393 END DO
2394
2395
2396
2397
2398
2399
2400
2401
2402
2403
2404
2405
2406
2407
2408
2409 cff=dtdays*abs(wbio(isink))
2411 DO i=istr,iend
2412 fc(i,k-1)=0.0_r8
2413 wl(i,k)=z_w(i,j,k-1)+cff
2414 wr(i,k)=hz(i,j,k)*qc(i,k)
2415 ksource(i,k)=k
2416 END DO
2417 END DO
2420 DO i=istr,iend
2421 IF (wl(i,k).gt.z_w(i,j,ks)) THEN
2422 ksource(i,k)=ks+1
2423 fc(i,k-1)=fc(i,k-1)+wr(i,ks)
2424 END IF
2425 END DO
2426 END DO
2427 END DO
2428
2429
2430
2432 DO i=istr,iend
2433 ks=ksource(i,k)
2434 cu=min(1.0_r8,(wl(i,k)-z_w(i,j,ks-1))*hz_inv(i,ks))
2435 fc(i,k-1)=fc(i,k-1)+ &
2436 & hz(i,j,ks)*cu* &
2437 & (bl(i,ks)+ &
2438 & cu*(0.5_r8*(br(i,ks)-bl(i,ks))- &
2439 & (1.5_r8-cu)* &
2440 & (br(i,ks)+bl(i,ks)- &
2441 & 2.0_r8*qc(i,ks))))
2442 END DO
2443 END DO
2445 DO i=istr,iend
2446 bio(i,k,ibio)=qc(i,k)+ &
2447 & (fc(i,k)-fc(i,k-1))*hz_inv(i,k)
2448 END DO
2449 END DO
2450 END DO
2451 END IF
2452 END DO
2453
2454
2455
2456
2457
2458
2459
2460
2461
2462
2463
2464 sink_loop: DO isink=1,nsink
2465 ibio=idsink(isink)
2466
2467
2468
2469
2470
2471
2473 DO i=istr,iend
2474 qc(i,k)=bio(i,k,ibio)
2475 tl_qc(i,k)=tl_bio(i,k,ibio)
2476 END DO
2477 END DO
2478
2480 DO i=istr,iend
2481 fc(i,k)=(qc(i,k+1)-qc(i,k))*hz_inv2(i,k)
2482 tl_fc(i,k)=(tl_qc(i,k+1)-tl_qc(i,k))*hz_inv2(i,k)+ &
2483 & (qc(i,k+1)-qc(i,k))*tl_hz_inv2(i,k)- &
2484#ifdef TL_IOMS
2485 & fc(i,k)
2486#endif
2487 END DO
2488 END DO
2490 DO i=istr,iend
2491 dltr=hz(i,j,k)*fc(i,k)
2492 tl_dltr=tl_hz(i,j,k)*fc(i,k)+hz(i,j,k)*tl_fc(i,k)- &
2493#ifdef TL_IOMS
2494 & dltr
2495#endif
2496 dltl=hz(i,j,k)*fc(i,k-1)
2497 tl_dltl=tl_hz(i,j,k)*fc(i,k-1)+hz(i,j,k)*tl_fc(i,k-1)- &
2498#ifdef TL_IOMS
2499 & dltl
2500#endif
2501 cff=hz(i,j,k-1)+2.0_r8*hz(i,j,k)+hz(i,j,k+1)
2502 tl_cff=tl_hz(i,j,k-1)+2.0_r8*tl_hz(i,j,k)+tl_hz(i,j,k+1)
2503 cffr=cff*fc(i,k)
2504 tl_cffr=tl_cff*fc(i,k)+cff*tl_fc(i,k)- &
2505#ifdef TL_IOMS
2506 & cffr
2507#endif
2508 cffl=cff*fc(i,k-1)
2509 tl_cffl=tl_cff*fc(i,k-1)+cff*tl_fc(i,k-1)- &
2510#ifdef TL_IOMS
2511 & cffl
2512#endif
2513
2514
2515
2516
2517 IF ((dltr*dltl).le.0.0_r8) THEN
2518 dltr=0.0_r8
2519 tl_dltr=0.0_r8
2520 dltl=0.0_r8
2521 tl_dltl=0.0_r8
2522 ELSE IF (abs(dltr).gt.abs(cffl)) THEN
2523 dltr=cffl
2524 tl_dltr=tl_cffl
2525 ELSE IF (abs(dltl).gt.abs(cffr)) THEN
2526 dltl=cffr
2527 tl_dltl=tl_cffr
2528 END IF
2529
2530
2531
2532
2533
2534
2535
2536
2537
2538
2539 cff=(dltr-dltl)*hz_inv3(i,k)
2540 tl_cff=(tl_dltr-tl_dltl)*hz_inv3(i,k)+ &
2541 & (dltr-dltl)*tl_hz_inv3(i,k)- &
2542#ifdef TL_IOMS
2543 & cff
2544#endif
2545 dltr=dltr-cff*hz(i,j,k+1)
2546 tl_dltr=tl_dltr-tl_cff*hz(i,j,k+1)-cff*tl_hz(i,j,k+1)+ &
2547#ifdef TL_IOMS
2548 & cff*hz(i,j,k+1)
2549#endif
2550 dltl=dltl+cff*hz(i,j,k-1)
2551 tl_dltl=tl_dltl+tl_cff*hz(i,j,k-1)+cff*tl_hz(i,j,k-1)- &
2552#ifdef TL_IOMS
2553 & cff*hz(i,j,k-1)
2554#endif
2555 br(i,k)=qc(i,k)+dltr
2556 tl_br(i,k)=tl_qc(i,k)+tl_dltr
2557 bl(i,k)=qc(i,k)-dltl
2558 tl_bl(i,k)=tl_qc(i,k)-tl_dltl
2559 wr(i,k)=(2.0_r8*dltr-dltl)**2
2560 tl_wr(i,k)=2.0_r8*(2.0_r8*dltr-dltl)* &
2561 & (2.0_r8*tl_dltr-tl_dltl)- &
2562#ifdef TL_IOMS
2563 & wr(i,k)
2564#endif
2565 wl(i,k)=(dltr-2.0_r8*dltl)**2
2566 tl_wl(i,k)=2.0_r8*(dltr-2.0_r8*dltl)* &
2567 & (tl_dltr-2.0_r8*tl_dltl)- &
2568#ifdef TL_IOMS
2569 & wl(i,k)
2570#endif
2571 END DO
2572 END DO
2573 cff=1.0e-14_r8
2575 DO i=istr,iend
2576 dltl=max(cff,wl(i,k ))
2577 tl_dltl=(0.5_r8-sign(0.5_r8,cff-wl(i,k )))* &
2578 & tl_wl(i,k )+ &
2579#ifdef TL_IOMS
2580 & cff*(0.5_r8+sign(0.5_r8,cff-wl(i,k )))
2581#endif
2582 dltr=max(cff,wr(i,k+1))
2583 tl_dltr=(0.5_r8-sign(0.5_r8,cff-wr(i,k+1)))* &
2584 & tl_wr(i,k+1)+ &
2585# ifdef TL_IOMS
2586 & cff*(0.5_r8+sign(0.5_r8,cff-wr(i,k+1)))
2587# endif
2588 br1(i,k)=br(i,k)
2589 bl1(i,k+1)=bl(i,k+1)
2590 br(i,k)=(dltr*br(i,k)+dltl*bl(i,k+1))/(dltr+dltl)
2591 tl_br(i,k)=(tl_dltr*br1(i,k )+dltr*tl_br(i,k )+ &
2592 & tl_dltl*bl1(i,k+1)+dltl*tl_bl(i,k+1))/ &
2593 & (dltr+dltl)- &
2594 & (tl_dltr+tl_dltl)*br(i,k)/(dltr+dltl)
2595 bl(i,k+1)=br(i,k)
2596 tl_bl(i,k+1)=tl_br(i,k)
2597 END DO
2598 END DO
2599 DO i=istr,iend
2601 tl_fc(i,
n(ng))=0.0_r8
2602#if defined LINEAR_CONTINUATION
2603 bl(i,
n(ng))=br(i,
n(ng)-1)
2604 tl_bl(i,
n(ng))=tl_br(i,
n(ng)-1)
2605 br(i,
n(ng))=2.0_r8*qc(i,
n(ng))-bl(i,
n(ng))
2606 tl_br(i,
n(ng))=2.0_r8*tl_qc(i,
n(ng))-tl_bl(i,
n(ng))
2607#elif defined NEUMANN
2608 bl(i,
n(ng))=br(i,
n(ng)-1)
2609 tl_bl(i,
n(ng))=tl_br(i,
n(ng)-1)
2610 br(i,
n(ng))=1.5_r8*qc(i,
n(ng))-0.5_r8*bl(i,
n(ng))
2611 tl_br(i,
n(ng))=1.5_r8*tl_qc(i,
n(ng))-0.5_r8*tl_bl(i,
n(ng))
2612#else
2613 br(i,
n(ng))=qc(i,
n(ng))
2614 bl(i,
n(ng))=qc(i,
n(ng))
2615 br(i,
n(ng)-1)=qc(i,
n(ng))
2616 tl_br(i,
n(ng))=tl_qc(i,
n(ng))
2617 tl_bl(i,
n(ng))=tl_qc(i,
n(ng))
2618 tl_br(i,
n(ng)-1)=tl_qc(i,
n(ng))
2619#endif
2620#if defined LINEAR_CONTINUATION
2621 br(i,1)=bl(i,2)
2622 tl_br(i,1)=tl_bl(i,2)
2623 bl(i,1)=2.0_r8*qc(i,1)-br(i,1)
2624 tl_bl(i,1)=2.0_r8*tl_qc(i,1)-tl_br(i,1)
2625#elif defined NEUMANN
2626 br(i,1)=bl(i,2)
2627 tl_br(i,1)=tl_bl(i,2)
2628 bl(i,1)=1.5_r8*qc(i,1)-0.5_r8*br(i,1)
2629 tl_bl(i,1)=1.5_r8*tl_qc(i,1)-0.5_r8*tl_br(i,1)
2630#else
2631 bl(i,2)=qc(i,1)
2632 br(i,1)=qc(i,1)
2633 bl(i,1)=qc(i,1)
2634 tl_bl(i,2)=tl_qc(i,1)
2635 tl_br(i,1)=tl_qc(i,1)
2636 tl_bl(i,1)=tl_qc(i,1)
2637#endif
2638 END DO
2639
2640
2641
2642
2643
2645 DO i=istr,iend
2646 dltr=br(i,k)-qc(i,k)
2647 tl_dltr=tl_br(i,k)-tl_qc(i,k)
2648 dltl=qc(i,k)-bl(i,k)
2649 tl_dltl=tl_qc(i,k)-tl_bl(i,k)
2650 cffr=2.0_r8*dltr
2651 tl_cffr=2.0_r8*tl_dltr
2652 cffl=2.0_r8*dltl
2653 tl_cffl=2.0_r8*tl_dltl
2654 IF ((dltr*dltl).lt.0.0_r8) THEN
2655 dltr=0.0_r8
2656 tl_dltr=0.0_r8
2657 dltl=0.0_r8
2658 tl_dltl=0.0_r8
2659 ELSE IF (abs(dltr).gt.abs(cffl)) THEN
2660 dltr=cffl
2661 tl_dltr=tl_cffl
2662 ELSE IF (abs(dltl).gt.abs(cffr)) THEN
2663 dltl=cffr
2664 tl_dltl=tl_cffr
2665 END IF
2666 br(i,k)=qc(i,k)+dltr
2667 tl_br(i,k)=tl_qc(i,k)+tl_dltr
2668 bl(i,k)=qc(i,k)-dltl
2669 tl_bl(i,k)=tl_qc(i,k)-tl_dltl
2670 END DO
2671 END DO
2672
2673
2674
2675
2676
2677
2678
2679
2680
2681
2682
2683
2684
2685
2686
2687 cff=dtdays*abs(wbio(isink))
2688 tl_cff=dtdays*sign(1.0_r8,wbio(isink))*tl_wbio(isink)
2690 DO i=istr,iend
2691 fc(i,k-1)=0.0_r8
2692 tl_fc(i,k-1)=0.0_r8
2693 wl(i,k)=z_w(i,j,k-1)+cff
2694 tl_wl(i,k)=tl_z_w(i,j,k-1)+tl_cff
2695 wr(i,k)=hz(i,j,k)*qc(i,k)
2696 tl_wr(i,k)=tl_hz(i,j,k)*qc(i,k)+hz(i,j,k)*tl_qc(i,k)- &
2697#ifdef TL_IOMS
2698 & wr(i,k)
2699#endif
2700 ksource(i,k)=k
2701 END DO
2702 END DO
2705 DO i=istr,iend
2706 IF (wl(i,k).gt.z_w(i,j,ks)) THEN
2707 ksource(i,k)=ks+1
2708 fc(i,k-1)=fc(i,k-1)+wr(i,ks)
2709 tl_fc(i,k-1)=tl_fc(i,k-1)+tl_wr(i,ks)
2710 END IF
2711 END DO
2712 END DO
2713 END DO
2714
2715
2716
2718 DO i=istr,iend
2719 ks=ksource(i,k)
2720 cu=min(1.0_r8,(wl(i,k)-z_w(i,j,ks-1))*hz_inv(i,ks))
2721 tl_cu=(0.5_r8+sign(0.5_r8, &
2722 & (1.0_r8-(wl(i,k)-z_w(i,j,ks-1))* &
2723 & hz_inv(i,ks))))* &
2724 & ((tl_wl(i,k)-tl_z_w(i,j,ks-1))*hz_inv(i,ks)+ &
2725 & (wl(i,k)-z_w(i,j,ks-1))*tl_hz_inv(i,ks)- &
2726#ifdef TL_IOMS
2727 & (wl(i,k)-z_w(i,j,ks-1))*hz_inv(i,ks) &
2728#endif
2729 & )+ &
2730#ifdef TL_IOMS
2731 & (0.5_r8-sign(0.5_r8, &
2732 & (1.0_r8-(wl(i,k)-z_w(i,j,ks-1))* &
2733 & hz_inv(i,ks))))
2734#endif
2735 fc(i,k-1)=fc(i,k-1)+ &
2736 & hz(i,j,ks)*cu* &
2737 & (bl(i,ks)+ &
2738 & cu*(0.5_r8*(br(i,ks)-bl(i,ks))- &
2739 & (1.5_r8-cu)* &
2740 & (br(i,ks)+bl(i,ks)- &
2741 & 2.0_r8*qc(i,ks))))
2742 tl_fc(i,k-1)=tl_fc(i,k-1)+ &
2743 & (tl_hz(i,j,ks)*cu+hz(i,j,ks)*tl_cu)* &
2744 & (bl(i,ks)+ &
2745 & cu*(0.5_r8*(br(i,ks)-bl(i,ks))- &
2746 & (1.5_r8-cu)* &
2747 & (br(i,ks)+bl(i,ks)- &
2748 & 2.0_r8*qc(i,ks))))+ &
2749 & hz(i,j,ks)*cu* &
2750 & (tl_bl(i,ks)+ &
2751 & tl_cu*(0.5_r8*(br(i,ks)-bl(i,ks))- &
2752 & (1.5_r8-cu)* &
2753 & (br(i,ks)+bl(i,ks)- &
2754 & 2.0_r8*qc(i,ks)))+ &
2755 & cu*(0.5_r8*(tl_br(i,ks)-tl_bl(i,ks))+ &
2756 & tl_cu* &
2757 & (br(i,ks)+bl(i,ks)-2.0_r8*qc(i,ks))- &
2758 & (1.5_r8-cu)* &
2759 & (tl_br(i,ks)+tl_bl(i,ks)- &
2760 & 2.0_r8*tl_qc(i,ks))))- &
2761#ifdef TL_IOMS
2762 & hz(i,j,ks)*cu* &
2763 & (2.0_r8*bl(i,ks)+ &
2764 & cu*(1.5_r8*(br(i,ks)-bl(i,ks))- &
2765 & (4.5_r8-4.0_r8*cu)* &
2766 & (br(i,ks)+bl(i,ks)- &
2767 & 2.0_r8*qc(i,ks))))
2768#endif
2769 END DO
2770 END DO
2772 DO i=istr,iend
2773 bio(i,k,ibio)=qc(i,k)+(fc(i,k)-fc(i,k-1))*hz_inv(i,k)
2774 tl_bio(i,k,ibio)=tl_qc(i,k)+ &
2775 & (tl_fc(i,k)-tl_fc(i,k-1))*hz_inv(i,k)+ &
2776 & (fc(i,k)-fc(i,k-1))*tl_hz_inv(i,k)- &
2777#ifdef TL_IOMS
2778 & (fc(i,k)-fc(i,k-1))*hz_inv(i,k)
2779#endif
2780 END DO
2781 END DO
2782
2783 END DO sink_loop
2784 END DO iter_loop
2785
2786
2787
2788
2789
2790
2791
2792
2793
2794
2795
2796
2797
2798
2799
2800
2804 DO i=istr,iend
2805 cff=bio(i,k,ibio)-bio_old(i,k,ibio)
2806 tl_cff=tl_bio(i,k,ibio)-tl_bio_old(i,k,ibio)
2807
2808
2809 tl_t(i,j,k,nnew,ibio)=tl_t(i,j,k,nnew,ibio)+ &
2810 & tl_cff*hz(i,j,k)+cff*tl_hz(i,j,k)- &
2811#ifdef TL_IOMS
2812 & cff*hz(i,j,k)
2813#endif
2814 END DO
2815 END DO
2816 END DO
2817
2818 END DO j_loop
2819
2820 RETURN
real(r8), dimension(:), allocatable parfrac
real(r8), dimension(:), allocatable phymrd
real(r8), dimension(:), allocatable zooeed
real(r8), dimension(:), allocatable ferr
real(r8), dimension(:), allocatable phyis
real(r8), dimension(:), allocatable attsw
real(r8), dimension(:), allocatable a_fe
real(r8), dimension(:), allocatable tl_parfrac
real(r8), dimension(:), allocatable tl_wphy
real(r8), dimension(:), allocatable fenudgtime
real(r8), dimension(:), allocatable ivlev
real(r8), dimension(:), allocatable attphy
real(r8), dimension(:), allocatable b_fe
real(r8), dimension(:), allocatable wphy
real(r8), dimension(:), allocatable fehmin
real(r8), dimension(:), allocatable k_fec
real(r8), dimension(:), allocatable t_fe
real(r8), dimension(:), allocatable zoomrn
real(r8), dimension(:), allocatable femax
real(r8), dimension(:), allocatable phymrn
real(r8), dimension(:), allocatable zoomrd
real(r8), dimension(:), allocatable zooeen