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
277
278 j_loop : DO j=jstr,jend
279
280
281
283 DO i=istr,iend
284 hz_inv(i,k)=1.0_r8/hz(i,j,k)
285 tl_hz_inv(i,k)=-hz_inv(i,k)*hz_inv(i,k)*tl_hz(i,j,k)
286 END DO
287 END DO
289 DO i=istr,iend
290 hz_inv2(i,k)=1.0_r8/(hz(i,j,k)+hz(i,j,k+1))
291 tl_hz_inv2(i,k)=-hz_inv2(i,k)*hz_inv2(i,k)* &
292 & (tl_hz(i,j,k)+tl_hz(i,j,k+1))
293 END DO
294 END DO
296 DO i=istr,iend
297 hz_inv3(i,k)=1.0_r8/(hz(i,j,k-1)+hz(i,j,k)+hz(i,j,k+1))
298 tl_hz_inv3(i,k)=-hz_inv3(i,k)*hz_inv3(i,k)* &
299 & (tl_hz(i,j,k-1)+tl_hz(i,j,k)+ &
300 & tl_hz(i,j,k+1))
301 END DO
302 END DO
303
304
305
309 DO i=istr,iend
310 bio(i,k,ibio)=0.0_r8
311 bio1(i,k,ibio)=0.0_r8
312 bio2(i,k,ibio)=0.0_r8
313 tl_bio(i,k,ibio)=0.0_r8
314 END DO
315 END DO
316 END DO
317
318
319
320
321
322
324 DO i=istr,iend
325
326
327
328
329
330
331
332
333
334
335
336
339
340
341 biotrc(ibio,nstp)=t(i,j,k,nstp,ibio)
342 tl_biotrc(ibio,nstp)=tl_t(i,j,k,nstp,ibio)
343
344
345 biotrc(ibio,nnew)=t(i,j,k,nnew,ibio)
346 tl_biotrc(ibio,nnew)=tl_t(i,j,k,nnew,ibio)* &
347 & hz_inv(i,k)+ &
348 & t(i,j,k,nnew,ibio)*hz(i,j,k)* &
349 & tl_hz_inv(i,k)
350 END DO
351
352
353
354 cff2=0.0_r8
355 DO itime=1,2
356 cff1=0.0_r8
357 tl_cff1=0.0_r8
359#ifdef IRON_LIMIT
361#else
363#endif
365 cff1=cff1+max(0.0_r8,minval-biotrc(ibio,itime))
366 tl_cff1=tl_cff1- &
367 & (0.5_r8-sign(0.5_r8, &
368 & biotrc(ibio,itime)-minval))* &
369 & tl_biotrc(ibio,itime)
370 IF (biotrc(ibio,itime).gt.biotrc(itrcmax,itime)) THEN
371 itrcmax=ibio
372 END IF
373 biotrc1(ibio,itime)=biotrc(ibio,itime)
374 biotrc(ibio,itime)=max(minval,biotrc1(ibio,itime))
375 tl_biotrc(ibio,itime)=(0.5_r8- &
376 & sign(0.5_r8, &
377 & minval- &
378 & biotrc1(ibio,itime)))* &
379 & tl_biotrc(ibio,itime)
380 END DO
381 IF (biotrc(itrcmax,itime).gt.cff1) THEN
382 biotrc(itrcmax,itime)=biotrc(itrcmax,itime)-cff1
383 tl_biotrc(itrcmax,itime)=tl_biotrc(itrcmax,itime)- &
384 & tl_cff1
385 END IF
386#ifdef IRON_LIMIT
389 biotrc1(ibio,itime)=biotrc(ibio,itime)
390 biotrc(ibio,itime)=max(minval,biotrc1(ibio,itime))
391 tl_biotrc(ibio,itime)=(0.5_r8- &
392 & sign(0.5_r8, &
393 & minval- &
394 & biotrc1(ibio,itime)))* &
395 & tl_biotrc(ibio,itime)
396 END DO
397#endif
398 END DO
399
400
401
404 bio_old(i,k,ibio)=biotrc(ibio,nstp)
405 tl_bio_old(i,k,ibio)=tl_biotrc(ibio,nstp)
406 bio(i,k,ibio)=biotrc(ibio,nstp)
407 tl_bio(i,k,ibio)=tl_biotrc(ibio,nstp)
408 END DO
409
410#if defined IRON_LIMIT && defined IRON_RELAX
411
412
413
414
415
416 IF (h(i,j).le.
fehmin(ng))
THEN
417
418
419
421 & fenudgcoef*tl_bio(i,k,
ifdis)
422 END IF
423#endif
424 END DO
425 END DO
426
427
428
429
430
431 DO i=istr,iend
432#ifdef CONST_PAR
433
434
435
436 parsur(i)=158.075_r8
437 tl_parsur(i)=0.0_r8
438#else
442#endif
443 END DO
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 iter_loop:
DO iter=1,
bioiter(ng)
500
501
502
504 DO i=istr,iend
505
506
507
508
509
510
511
512
513
514
515
516
519
520
521 biotrc(ibio,nstp)=t(i,j,k,nstp,ibio)
522
523
524 biotrc(ibio,nnew)=t(i,j,k,nnew,ibio)
525 END DO
526
527
528
529 cff2=0.0_r8
530 DO itime=1,2
531 cff1=0.0_r8
533#ifdef IRON_LIMIT
535#else
537#endif
539 cff1=cff1+max(0.0_r8,minval-biotrc(ibio,itime))
540 IF (biotrc(ibio,itime).gt.biotrc(itrcmax,itime)) THEN
541 itrcmax=ibio
542 END IF
543 biotrc(ibio,itime)=max(minval,biotrc(ibio,itime))
544 END DO
545 IF (biotrc(itrcmax,itime).gt.cff1) THEN
546 biotrc(itrcmax,itime)=biotrc(itrcmax,itime)-cff1
547 END IF
548#ifdef IRON_LIMIT
551 biotrc1(ibio,itime)=biotrc(ibio,itime)
552 biotrc(ibio,itime)=max(minval,biotrc(ibio,itime))
553 END DO
554#endif
555 END DO
556
557
558
561 bio_old(i,k,ibio)=biotrc(ibio,nstp)
562 bio(i,k,ibio)=biotrc(ibio,nstp)
563 END DO
564
565#if defined IRON_LIMIT && defined IRON_RELAX
566
567
568
569
570
571 IF (h(i,j).le.
fehmin(ng))
THEN
574 END IF
575#endif
576 END DO
577 END DO
578
579
580
581
582
583 DO i=istr,iend
584#ifdef CONST_PAR
585
586
587
588 parsur(i)=158.075_r8
589#else
591#endif
592 END DO
593
594
595
596
597
598
599 DO iteradj=1,iter
600
601
602
603 DO i=istr,iend
604 par=parsur(i)
605 IF (parsur(i).gt.0.0_r8) THEN
607
608
609
610
611
613 & (z_w(i,j,k)-z_w(i,j,k-1))
614 expatt=exp(-att)
615 itop=par
616 par=itop*(1.0_r8-expatt)/att
617 light(i,k)=par
618
619
620
621
622 par=itop*expatt
623 END DO
624 ELSE
626 light(i,k)=0.0_r8
627 END DO
628 END IF
629 END DO
630
631
632
633
634
635
636#ifdef IRON_LIMIT
637
638
639
640
641
642
643
644
645
646
647#endif
648
653 DO i=istr,iend
654#ifdef IRON_LIMIT
655
656
657
658 fnratio=bio(i,k,
ifphy)/max(minval,bio(i,k,
iphyt))
659 fcratio=fnratio*fen2fec
661 flimit=fcratio*fcratio/ &
663
665 fnlim=min(1.0_r8,flimit/(bio(i,k,
ino3_)*nlimit))
666#endif
667 cff4=1.0_r8/sqrt(cff2+cff3*light(i,k)*light(i,k))
668 cff=bio(i,k,
iphyt)* &
669#ifdef IRON_LIMIT
670 & cff1*cff4*light(i,k)*fnlim*nlimit
671#else
672 & cff1*cff4*light(i,k)/ &
674#endif
680
681#ifdef IRON_LIMIT
682
683
684
685 fac=cff*bio(i,k,
ino3_)*fnratio/ &
686 & max(minval,bio(i,k,
ifdis))
694
695
696
697 cff5=dtdays*(fcratioe-fcratio)/
t_fe(ng)
698 cff6=bio(i,k,
iphyt)*cff5*fec2fen
699 IF (cff6.ge.0.0_r8) THEN
700 cff=cff6/max(minval,bio(i,k,
ifdis))
704 ELSE
705 cff=-cff6/max(minval,bio(i,k,
ifphy))
709 END IF
710#endif
711 END DO
712 END DO
713
714 IF (iteradj.ne.iter) THEN
715
716
717
718
719
720#ifdef IRON_LIMIT
721
722
723#endif
724
725 cff1=dtdays*
zoogr(ng)
728 DO i=istr,iend
729 cff=bio(i,k,
izoop)* &
730 & cff1*(1.0_r8-exp(-
ivlev(ng)*bio(i,k,
iphyt)))/ &
734 & bio(i,k,
iphyt)*cff2*cff
739#ifdef IRON_LIMIT
743#endif
744 END DO
745 END DO
746
747
748
749
752 cff1=1.0_r8/(1.0_r8+cff2+cff3)
754 DO i=istr,iend
757 & bio(i,k,
iphyt)*cff2
759 & bio(i,k,
iphyt)*cff3
760#ifdef IRON_LIMIT
764#endif
765 END DO
766 END DO
767
768
769
770
773 cff1=1.0_r8/(1.0_r8+cff2+cff3)
775 DO i=istr,iend
778 & bio(i,k,
izoop)*cff2
780 & bio(i,k,
izoop)*cff3
781 END DO
782 END DO
783
784
785
786 cff2=dtdays*
detrr(ng)
787 cff1=1.0_r8/(1.0_r8+cff2)
789 DO i=istr,iend
792 & bio(i,k,
isdet)*cff2
793 END DO
794 END DO
795
796
797
798
799
800
801
802
803
804 DO isink=1,nsink
805 ibio=idsink(isink)
806
807
808
809
810
811
813 DO i=istr,iend
814 qc(i,k)=bio(i,k,ibio)
815 END DO
816 END DO
817
819 DO i=istr,iend
820 fc(i,k)=(qc(i,k+1)-qc(i,k))*hz_inv2(i,k)
821 END DO
822 END DO
824 DO i=istr,iend
825 dltr=hz(i,j,k)*fc(i,k)
826 dltl=hz(i,j,k)*fc(i,k-1)
827 cff=hz(i,j,k-1)+2.0_r8*hz(i,j,k)+hz(i,j,k+1)
828 cffr=cff*fc(i,k)
829 cffl=cff*fc(i,k-1)
830
831
832
833
834 IF ((dltr*dltl).le.0.0_r8) THEN
835 dltr=0.0_r8
836 dltl=0.0_r8
837 ELSE IF (abs(dltr).gt.abs(cffl)) THEN
838 dltr=cffl
839 ELSE IF (abs(dltl).gt.abs(cffr)) THEN
840 dltl=cffr
841 END IF
842
843
844
845
846
847
848
849
850
851
852 cff=(dltr-dltl)*hz_inv3(i,k)
853 dltr=dltr-cff*hz(i,j,k+1)
854 dltl=dltl+cff*hz(i,j,k-1)
855 br(i,k)=qc(i,k)+dltr
856 bl(i,k)=qc(i,k)-dltl
857 wr(i,k)=(2.0_r8*dltr-dltl)**2
858 wl(i,k)=(dltr-2.0_r8*dltl)**2
859 END DO
860 END DO
861 cff=1.0e-14_r8
863 DO i=istr,iend
864 dltl=max(cff,wl(i,k ))
865 dltr=max(cff,wr(i,k+1))
866 br(i,k)=(dltr*br(i,k)+dltl*bl(i,k+1))/(dltr+dltl)
867 bl(i,k+1)=br(i,k)
868 END DO
869 END DO
870 DO i=istr,iend
872#if defined LINEAR_CONTINUATION
873 bl(i,
n(ng))=br(i,
n(ng)-1)
874 br(i,
n(ng))=2.0_r8*qc(i,
n(ng))-bl(i,
n(ng))
875#elif defined NEUMANN
876 bl(i,
n(ng))=br(i,
n(ng)-1)
877 br(i,
n(ng))=1.5*qc(i,
n(ng))-0.5_r8*bl(i,
n(ng))
878#else
879 br(i,
n(ng))=qc(i,
n(ng))
880 bl(i,
n(ng))=qc(i,
n(ng))
881 br(i,
n(ng)-1)=qc(i,
n(ng))
882#endif
883#if defined LINEAR_CONTINUATION
884 br(i,1)=bl(i,2)
885 bl(i,1)=2.0_r8*qc(i,1)-br(i,1)
886#elif defined NEUMANN
887 br(i,1)=bl(i,2)
888 bl(i,1)=1.5_r8*qc(i,1)-0.5_r8*br(i,1)
889#else
890 bl(i,2)=qc(i,1)
891 br(i,1)=qc(i,1)
892 bl(i,1)=qc(i,1)
893#endif
894 END DO
895
896
897
898
899
901 DO i=istr,iend
902 dltr=br(i,k)-qc(i,k)
903 dltl=qc(i,k)-bl(i,k)
904 cffr=2.0_r8*dltr
905 cffl=2.0_r8*dltl
906 IF ((dltr*dltl).lt.0.0_r8) THEN
907 dltr=0.0_r8
908 dltl=0.0_r8
909 ELSE IF (abs(dltr).gt.abs(cffl)) THEN
910 dltr=cffl
911 ELSE IF (abs(dltl).gt.abs(cffr)) THEN
912 dltl=cffr
913 END IF
914 br(i,k)=qc(i,k)+dltr
915 bl(i,k)=qc(i,k)-dltl
916 END DO
917 END DO
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933 cff=dtdays*abs(wbio(isink))
935 DO i=istr,iend
936 fc(i,k-1)=0.0_r8
937 wl(i,k)=z_w(i,j,k-1)+cff
938 wr(i,k)=hz(i,j,k)*qc(i,k)
939 ksource(i,k)=k
940 END DO
941 END DO
944 DO i=istr,iend
945 IF (wl(i,k).gt.z_w(i,j,ks)) THEN
946 ksource(i,k)=ks+1
947 fc(i,k-1)=fc(i,k-1)+wr(i,ks)
948 END IF
949 END DO
950 END DO
951 END DO
952
953
954
956 DO i=istr,iend
957 ks=ksource(i,k)
958 cu=min(1.0_r8,(wl(i,k)-z_w(i,j,ks-1))*hz_inv(i,ks))
959 fc(i,k-1)=fc(i,k-1)+ &
960 & hz(i,j,ks)*cu* &
961 & (bl(i,ks)+ &
962 & cu*(0.5_r8*(br(i,ks)-bl(i,ks))- &
963 & (1.5_r8-cu)* &
964 & (br(i,ks)+bl(i,ks)- &
965 & 2.0_r8*qc(i,ks))))
966 END DO
967 END DO
969 DO i=istr,iend
970 bio(i,k,ibio)=qc(i,k)+ &
971 & (fc(i,k)-fc(i,k-1))*hz_inv(i,k)
972 END DO
973 END DO
974 END DO
975 END IF
976 END DO
977
978
979
980
981
982 DO i=istr,iend
983 par=parsur(i)
984 tl_par=tl_parsur(i)
985 IF (parsur(i).gt.0.0_r8) THEN
987
988
989
990
991
993 & (z_w(i,j,k)-z_w(i,j,k-1))
995 & (z_w(i,j,k)-z_w(i,j,k-1))+ &
997 & (tl_z_w(i,j,k)-tl_z_w(i,j,k-1))
998 expatt=exp(-att)
999 tl_expatt=-expatt*tl_att
1000 itop=par
1001 tl_itop=tl_par
1002 par=itop*(1.0_r8-expatt)/att
1003 tl_par=(-tl_att*par+tl_itop*(1.0_r8-expatt)- &
1004 & itop*tl_expatt)/att
1005
1006
1007 tl_light(i,k)=tl_par
1008
1009
1010
1011
1012 par=itop*expatt
1013 tl_par=tl_itop*expatt+itop*tl_expatt
1014 END DO
1015 ELSE
1017
1018
1019 tl_light(i,k)=0.0_r8
1020 END DO
1021 END IF
1022 END DO
1023
1024
1025
1026
1027
1028
1029#ifdef IRON_LIMIT
1030
1031
1032
1033
1034
1035
1036
1037
1038
1039
1040#endif
1041
1046 DO i=istr,iend
1047#ifdef IRON_LIMIT
1048
1049
1050
1051
1052
1053 fac1=max(minval,bio1(i,k,
iphyt))
1054 tl_fac1=(0.5_r8-sign(0.5_r8,minval-bio1(i,k,
iphyt)))* &
1056 fnratio=bio1(i,k,
ifphy)/fac1
1057 tl_fnratio=(tl_bio(i,k,
ifphy)-tl_fac1*fnratio)/fac1
1058 fcratio=fnratio*fen2fec
1059 tl_fcratio=tl_fnratio*fen2fec
1064 flimit=fcratio*fcratio/ &
1066 tl_flimit=2.0_r8*(tl_fcratio*fcratio- &
1067 & tl_fcratio*fcratio*flimit)/ &
1070 tl_nlimit=-tl_bio(i,k,
ino3_)*nlimit*nlimit
1071
1072
1073 fac1=flimit/(bio1(i,k,
ino3_)*nlimit)
1074 tl_fac1=tl_flimit/(bio1(i,k,
ino3_)*nlimit)- &
1075 & (tl_bio(i,k,
ino3_)*nlimit+ &
1076 & bio1(i,k,
ino3_)*tl_nlimit)*fac1/ &
1077 & (bio1(i,k,
ino3_)*nlimit)
1078 fnlim=min(1.0_r8,fac1)
1079 tl_fnlim=(0.5_r8+sign(0.5_r8,1.0_r8-fac1))*tl_fac1
1080#endif
1081 cff4=1.0_r8/sqrt(cff2+cff3*light(i,k)*light(i,k))
1082 tl_cff4=-cff3*tl_light(i,k)*light(i,k)*cff4*cff4*cff4
1083 cff=bio1(i,k,
iphyt)* &
1084#ifdef IRON_LIMIT
1085 & cff1*cff4*light(i,k)*fnlim*nlimit
1086#else
1087 & cff1*cff4*light(i,k)/ &
1089#endif
1090#ifdef IRON_LIMIT
1091 tl_cff=tl_bio(i,k,
iphyt)* &
1092 & cff1*cff4*light(i,k)*fnlim*nlimit+ &
1093 & bio1(i,k,
iphyt)*cff1*cff4* &
1094 & (tl_light(i,k)*fnlim*nlimit+ &
1095 & light(i,k)*tl_fnlim*nlimit+ &
1096 & light(i,k)*fnlim*tl_nlimit)+ &
1097 & bio1(i,k,
iphyt)*cff1*tl_cff4* &
1098 & light(i,k)*fnlim*nlimit
1099#else
1100 tl_cff=(tl_bio(i,k,
iphyt)*cff1*cff4*light(i,k)+ &
1101 & bio1(i,k,
iphyt)*cff1* &
1102 & (tl_cff4*light(i,k)+cff4*tl_light(i,k))- &
1103 & tl_bio(i,k,
ino3_)*cff)/ &
1105#endif
1106
1107
1109 & tl_cff*bio(i,k,
ino3_))/ &
1110 & (1.0_r8+cff)
1111
1112
1113
1115 & tl_bio(i,k,
ino3_)*cff+ &
1116 & bio(i,k,
ino3_)*tl_cff
1117
1118#ifdef IRON_LIMIT
1119
1120
1121
1122
1123
1124 fac1=max(minval,bio1(i,k,
ifdis))
1125 tl_fac1=(0.5_r8-sign(0.5_r8,minval-bio1(i,k,
ifdis)))* &
1127 fac2=1.0_r8/fac1
1128 tl_fac2=-fac2*fac2*tl_fac1
1129 fac=cff*bio(i,k,
ino3_)*fnratio*fac2
1130 tl_fac=fnratio*fac2*(tl_cff*bio(i,k,
ino3_)+ &
1131 & cff*tl_bio(i,k,
ino3_))+ &
1132 & cff*bio(i,k,
ino3_)*(tl_fnratio*fac2+ &
1133 & fnratio*tl_fac2)
1134
1135
1137 & tl_fac*bio2(i,k,
ifdis))/ &
1138 & (1.0_r8+fac)
1139
1140
1141
1143 & tl_bio(i,k,
ifdis)*fac+ &
1144 & bio2(i,k,
ifdis)*tl_fac
1145
1146
1147
1148 cff5=dtdays*(fcratioe-fcratio)/
t_fe(ng)
1149 tl_cff5=dtdays*(tl_fcratioe-tl_fcratio)/
t_fe(ng)
1150 cff6=bio(i,k,
iphyt)*cff5*fec2fen
1151 tl_cff6=(tl_bio(i,k,
iphyt)*cff5+ &
1152 & bio(i,k,
iphyt)*tl_cff5)*fec2fen
1153 IF (cff6.ge.0.0_r8) THEN
1154
1155
1156 fac1=max(minval,bio2(i,k,
ifdis))
1157 tl_fac1=(0.5_r8-sign(0.5_r8,minval-bio2(i,k,
ifdis)))* &
1159 cff=cff6/fac1
1160 tl_cff=(tl_cff6-tl_fac1*cff)/fac1
1161
1162
1164 & tl_cff*bio(i,k,
ifdis))/ &
1165 & (1.0_r8+cff)
1166
1167
1168
1170 & tl_bio(i,k,
ifdis)*cff+ &
1171 & bio(i,k,
ifdis)*tl_cff
1172 ELSE
1173
1174
1175 fac1=-max(minval,bio2(i,k,
ifphy))
1176 tl_fac1=-(0.5_r8-sign(0.5_r8,minval-bio2(i,k,
ifphy)))* &
1178 cff=cff6/fac1
1179 tl_cff=(tl_cff6-tl_fac1*cff)/fac1
1180
1181
1183 & tl_cff*bio(i,k,
ifphy))/ &
1184 & (1.0_r8+cff)
1185
1186
1187
1189 & tl_bio(i,k,
ifphy)*cff+ &
1190 & bio(i,k,
ifphy)*tl_cff
1191 END IF
1192#endif
1193 END DO
1194 END DO
1195
1196
1197
1199 DO i=istr,iend
1200
1201
1202
1203
1204
1205
1206
1207
1208
1209
1210
1211
1214
1215
1216 biotrc(ibio,nstp)=t(i,j,k,nstp,ibio)
1217
1218
1219 biotrc(ibio,nnew)=t(i,j,k,nnew,ibio)
1220 END DO
1221
1222
1223
1224 cff2=0.0_r8
1225 DO itime=1,2
1226 cff1=0.0_r8
1228#ifdef IRON_LIMIT
1230#else
1232#endif
1234 cff1=cff1+max(0.0_r8,minval-biotrc(ibio,itime))
1235 IF (biotrc(ibio,itime).gt.biotrc(itrcmax,itime)) THEN
1236 itrcmax=ibio
1237 END IF
1238 biotrc(ibio,itime)=max(minval,biotrc(ibio,itime))
1239 END DO
1240 IF (biotrc(itrcmax,itime).gt.cff1) THEN
1241 biotrc(itrcmax,itime)=biotrc(itrcmax,itime)-cff1
1242 END IF
1243#ifdef IRON_LIMIT
1246 biotrc1(ibio,itime)=biotrc(ibio,itime)
1247 biotrc(ibio,itime)=max(minval,biotrc(ibio,itime))
1248 END DO
1249#endif
1250 END DO
1251
1252
1253
1256 bio_old(i,k,ibio)=biotrc(ibio,nstp)
1257 bio(i,k,ibio)=biotrc(ibio,nstp)
1258 END DO
1259
1260#if defined IRON_LIMIT && defined IRON_RELAX
1261
1262
1263
1264
1265
1266 IF (h(i,j).le.
fehmin(ng))
THEN
1269 END IF
1270#endif
1271 END DO
1272 END DO
1273
1274
1275
1276
1277
1278 DO i=istr,iend
1279#ifdef CONST_PAR
1280
1281
1282
1283 parsur(i)=158.075_r8
1284#else
1286#endif
1287 END DO
1288
1289
1290
1291
1292
1293
1294 DO iteradj=1,iter
1295
1296
1297
1298 DO i=istr,iend
1299 par=parsur(i)
1300 IF (parsur(i).gt.0.0_r8) THEN
1302
1303
1304
1305
1306
1308 & (z_w(i,j,k)-z_w(i,j,k-1))
1309 expatt=exp(-att)
1310 itop=par
1311 par=itop*(1.0_r8-expatt)/att
1312 light(i,k)=par
1313
1314
1315
1316
1317 par=itop*expatt
1318 END DO
1319 ELSE
1321 light(i,k)=0.0_r8
1322 END DO
1323 END IF
1324 END DO
1325
1326
1327
1328
1329
1330
1331#ifdef IRON_LIMIT
1332
1333
1334
1335
1336
1337
1338
1339
1340
1341
1342#endif
1343
1348 DO i=istr,iend
1349#ifdef IRON_LIMIT
1350
1351
1352
1353 fnratio=bio(i,k,
ifphy)/max(minval,bio(i,k,
iphyt))
1354 fcratio=fnratio*fen2fec
1356 flimit=fcratio*fcratio/ &
1358
1360 fnlim=min(1.0_r8,flimit/(bio(i,k,
ino3_)*nlimit))
1361#endif
1362 cff4=1.0_r8/sqrt(cff2+cff3*light(i,k)*light(i,k))
1363 cff=bio(i,k,
iphyt)* &
1364#ifdef IRON_LIMIT
1365 & cff1*cff4*light(i,k)*fnlim*nlimit
1366#else
1367 & cff1*cff4*light(i,k)/ &
1369#endif
1372 & bio(i,k,
ino3_)*cff
1373
1374#ifdef IRON_LIMIT
1375
1376
1377
1378 fac=cff*bio(i,k,
ino3_)*fnratio/ &
1379 & max(minval,bio(i,k,
ifdis))
1382 & bio(i,k,
ifdis)*fac
1383
1384
1385
1386 cff5=dtdays*(fcratioe-fcratio)/
t_fe(ng)
1387 cff6=bio(i,k,
iphyt)*cff5*fec2fen
1388 IF (cff6.ge.0.0_r8) THEN
1389 cff=cff6/max(minval,bio(i,k,
ifdis))
1392 & bio(i,k,
ifdis)*cff
1393 ELSE
1394 cff=-cff6/max(minval,bio(i,k,
ifphy))
1397 & bio(i,k,
ifphy)*cff
1398 END IF
1399#endif
1400 END DO
1401 END DO
1402
1403
1404
1405
1406
1407#ifdef IRON_LIMIT
1408
1409
1410#endif
1411
1412 cff1=dtdays*
zoogr(ng)
1415 DO i=istr,iend
1416 cff=bio(i,k,
izoop)* &
1417 & cff1*(1.0_r8-exp(-
ivlev(ng)*bio(i,k,
iphyt)))/ &
1423 & bio(i,k,
iphyt)*cff2*cff
1428#ifdef IRON_LIMIT
1434#endif
1435 END DO
1436 END DO
1437
1438
1439
1440
1443 cff1=1.0_r8/(1.0_r8+cff2+cff3)
1445 DO i=istr,iend
1448 & bio(i,k,
iphyt)*cff2
1450 & bio(i,k,
iphyt)*cff3
1451#ifdef IRON_LIMIT
1455#endif
1456 END DO
1457 END DO
1458
1459 IF (iteradj.ne.iter) THEN
1460
1461
1462
1463
1466 cff1=1.0_r8/(1.0_r8+cff2+cff3)
1468 DO i=istr,iend
1471 & bio(i,k,
izoop)*cff2
1473 & bio(i,k,
izoop)*cff3
1474 END DO
1475 END DO
1476
1477
1478
1479 cff2=dtdays*
detrr(ng)
1480 cff1=1.0_r8/(1.0_r8+cff2)
1482 DO i=istr,iend
1485 & bio(i,k,
isdet)*cff2
1486 END DO
1487 END DO
1488
1489
1490
1491
1492
1493
1494
1495
1496
1497 DO isink=1,nsink
1498 ibio=idsink(isink)
1499
1500
1501
1502
1503
1504
1506 DO i=istr,iend
1507 qc(i,k)=bio(i,k,ibio)
1508 END DO
1509 END DO
1510
1512 DO i=istr,iend
1513 fc(i,k)=(qc(i,k+1)-qc(i,k))*hz_inv2(i,k)
1514 END DO
1515 END DO
1517 DO i=istr,iend
1518 dltr=hz(i,j,k)*fc(i,k)
1519 dltl=hz(i,j,k)*fc(i,k-1)
1520 cff=hz(i,j,k-1)+2.0_r8*hz(i,j,k)+hz(i,j,k+1)
1521 cffr=cff*fc(i,k)
1522 cffl=cff*fc(i,k-1)
1523
1524
1525
1526
1527 IF ((dltr*dltl).le.0.0_r8) THEN
1528 dltr=0.0_r8
1529 dltl=0.0_r8
1530 ELSE IF (abs(dltr).gt.abs(cffl)) THEN
1531 dltr=cffl
1532 ELSE IF (abs(dltl).gt.abs(cffr)) THEN
1533 dltl=cffr
1534 END IF
1535
1536
1537
1538
1539
1540
1541
1542
1543
1544
1545 cff=(dltr-dltl)*hz_inv3(i,k)
1546 dltr=dltr-cff*hz(i,j,k+1)
1547 dltl=dltl+cff*hz(i,j,k-1)
1548 br(i,k)=qc(i,k)+dltr
1549 bl(i,k)=qc(i,k)-dltl
1550 wr(i,k)=(2.0_r8*dltr-dltl)**2
1551 wl(i,k)=(dltr-2.0_r8*dltl)**2
1552 END DO
1553 END DO
1554 cff=1.0e-14_r8
1556 DO i=istr,iend
1557 dltl=max(cff,wl(i,k ))
1558 dltr=max(cff,wr(i,k+1))
1559 br(i,k)=(dltr*br(i,k)+dltl*bl(i,k+1))/(dltr+dltl)
1560 bl(i,k+1)=br(i,k)
1561 END DO
1562 END DO
1563 DO i=istr,iend
1565#if defined LINEAR_CONTINUATION
1566 bl(i,
n(ng))=br(i,
n(ng)-1)
1567 br(i,
n(ng))=2.0_r8*qc(i,
n(ng))-bl(i,
n(ng))
1568#elif defined NEUMANN
1569 bl(i,
n(ng))=br(i,
n(ng)-1)
1570 br(i,
n(ng))=1.5*qc(i,
n(ng))-0.5_r8*bl(i,
n(ng))
1571#else
1572 br(i,
n(ng))=qc(i,
n(ng))
1573 bl(i,
n(ng))=qc(i,
n(ng))
1574 br(i,
n(ng)-1)=qc(i,
n(ng))
1575#endif
1576#if defined LINEAR_CONTINUATION
1577 br(i,1)=bl(i,2)
1578 bl(i,1)=2.0_r8*qc(i,1)-br(i,1)
1579#elif defined NEUMANN
1580 br(i,1)=bl(i,2)
1581 bl(i,1)=1.5_r8*qc(i,1)-0.5_r8*br(i,1)
1582#else
1583 bl(i,2)=qc(i,1)
1584 br(i,1)=qc(i,1)
1585 bl(i,1)=qc(i,1)
1586#endif
1587 END DO
1588
1589
1590
1591
1592
1594 DO i=istr,iend
1595 dltr=br(i,k)-qc(i,k)
1596 dltl=qc(i,k)-bl(i,k)
1597 cffr=2.0_r8*dltr
1598 cffl=2.0_r8*dltl
1599 IF ((dltr*dltl).lt.0.0_r8) THEN
1600 dltr=0.0_r8
1601 dltl=0.0_r8
1602 ELSE IF (abs(dltr).gt.abs(cffl)) THEN
1603 dltr=cffl
1604 ELSE IF (abs(dltl).gt.abs(cffr)) THEN
1605 dltl=cffr
1606 END IF
1607 br(i,k)=qc(i,k)+dltr
1608 bl(i,k)=qc(i,k)-dltl
1609 END DO
1610 END DO
1611
1612
1613
1614
1615
1616
1617
1618
1619
1620
1621
1622
1623
1624
1625
1626 cff=dtdays*abs(wbio(isink))
1628 DO i=istr,iend
1629 fc(i,k-1)=0.0_r8
1630 wl(i,k)=z_w(i,j,k-1)+cff
1631 wr(i,k)=hz(i,j,k)*qc(i,k)
1632 ksource(i,k)=k
1633 END DO
1634 END DO
1637 DO i=istr,iend
1638 IF (wl(i,k).gt.z_w(i,j,ks)) THEN
1639 ksource(i,k)=ks+1
1640 fc(i,k-1)=fc(i,k-1)+wr(i,ks)
1641 END IF
1642 END DO
1643 END DO
1644 END DO
1645
1646
1647
1649 DO i=istr,iend
1650 ks=ksource(i,k)
1651 cu=min(1.0_r8,(wl(i,k)-z_w(i,j,ks-1))*hz_inv(i,ks))
1652 fc(i,k-1)=fc(i,k-1)+ &
1653 & hz(i,j,ks)*cu* &
1654 & (bl(i,ks)+ &
1655 & cu*(0.5_r8*(br(i,ks)-bl(i,ks))- &
1656 & (1.5_r8-cu)* &
1657 & (br(i,ks)+bl(i,ks)- &
1658 & 2.0_r8*qc(i,ks))))
1659 END DO
1660 END DO
1662 DO i=istr,iend
1663 bio(i,k,ibio)=qc(i,k)+ &
1664 & (fc(i,k)-fc(i,k-1))*hz_inv(i,k)
1665 END DO
1666 END DO
1667 END DO
1668 END IF
1669 END DO
1670
1671
1672
1673
1674
1675
1676
1677#ifdef IRON_LIMIT
1678
1679
1680#endif
1681
1682 cff1=dtdays*
zoogr(ng)
1685 DO i=istr,iend
1686 cff=bio1(i,k,
izoop)* &
1687 & cff1*(1.0_r8-exp(-
ivlev(ng)*bio1(i,k,
iphyt)))/ &
1689 tl_cff=(tl_bio(i,k,
izoop)* &
1690 & cff1*(1.0_r8-exp(-
ivlev(ng)*bio1(i,k,
iphyt)))+ &
1693 & tl_bio(i,k,
iphyt)*cff)/ &
1695
1696
1698 & tl_cff*bio(i,k,
iphyt))/ &
1699 & (1.0_r8+cff)
1700
1701
1702
1704 & cff2*(tl_bio(i,k,
iphyt)*cff+ &
1705 & bio(i,k,
iphyt)*tl_cff)
1706
1707
1708
1711 & bio(i,k,
iphyt)*tl_cff)
1712
1713
1714
1717 & bio(i,k,
iphyt)*tl_cff)
1718
1719#ifdef IRON_LIMIT
1720
1721
1723 & tl_cff*bio2(i,k,
ifphy))/ &
1724 & (1.0_r8+cff)
1725
1726
1727
1729 & (tl_bio(i,k,
ifphy)*cff+ &
1731#endif
1732 END DO
1733 END DO
1734
1735
1736
1737
1740 cff1=1.0_r8/(1.0_r8+cff2+cff3)
1742 DO i=istr,iend
1743
1744
1746
1747
1748
1750 & tl_bio(i,k,
iphyt)*cff2
1751
1752
1753
1755 & tl_bio(i,k,
iphyt)*cff3
1756
1757#ifdef IRON_LIMIT
1758
1759
1761
1762
1763
1765 & tl_bio(i,k,
ifphy)*(cff2+cff3)*
ferr(ng)
1766#endif
1767 END DO
1768 END DO
1769
1770
1771
1772
1775 cff1=1.0_r8/(1.0_r8+cff2+cff3)
1777 DO i=istr,iend
1778
1779
1781
1782
1783
1785 & tl_bio(i,k,
izoop)*cff2
1786
1787
1788
1790 & tl_bio(i,k,
izoop)*cff3
1791 END DO
1792 END DO
1793
1794
1795
1796 cff2=dtdays*
detrr(ng)
1797 cff1=1.0_r8/(1.0_r8+cff2)
1799 DO i=istr,iend
1800
1801
1803
1804
1805
1807 & tl_bio(i,k,
isdet)*cff2
1808 END DO
1809 END DO
1810
1811
1812
1814 DO i=istr,iend
1815
1816
1817
1818
1819
1820
1821
1822
1823
1824
1825
1826
1829
1830
1831 biotrc(ibio,nstp)=t(i,j,k,nstp,ibio)
1832
1833
1834 biotrc(ibio,nnew)=t(i,j,k,nnew,ibio)
1835 END DO
1836
1837
1838
1839 cff2=0.0_r8
1840 DO itime=1,2
1841 cff1=0.0_r8
1843#ifdef IRON_LIMIT
1845#else
1847#endif
1849 cff1=cff1+max(0.0_r8,minval-biotrc(ibio,itime))
1850 IF (biotrc(ibio,itime).gt.biotrc(itrcmax,itime)) THEN
1851 itrcmax=ibio
1852 END IF
1853 biotrc(ibio,itime)=max(minval,biotrc(ibio,itime))
1854 END DO
1855 IF (biotrc(itrcmax,itime).gt.cff1) THEN
1856 biotrc(itrcmax,itime)=biotrc(itrcmax,itime)-cff1
1857 END IF
1858#ifdef IRON_LIMIT
1861 biotrc1(ibio,itime)=biotrc(ibio,itime)
1862 biotrc(ibio,itime)=max(minval,biotrc(ibio,itime))
1863 END DO
1864#endif
1865 END DO
1866
1867
1868
1871 bio_old(i,k,ibio)=biotrc(ibio,nstp)
1872 bio(i,k,ibio)=biotrc(ibio,nstp)
1873 END DO
1874
1875#if defined IRON_LIMIT && defined IRON_RELAX
1876
1877
1878
1879
1880
1881 IF (h(i,j).le.
fehmin(ng))
THEN
1884 END IF
1885#endif
1886 END DO
1887 END DO
1888
1889
1890
1891
1892
1893 DO i=istr,iend
1894#ifdef CONST_PAR
1895
1896
1897
1898 parsur(i)=158.075_r8
1899#else
1901#endif
1902 END DO
1903
1904
1905
1906
1907
1908
1909 DO iteradj=1,iter
1910
1911
1912
1913 DO i=istr,iend
1914 par=parsur(i)
1915 IF (parsur(i).gt.0.0_r8) THEN
1917
1918
1919
1920
1921
1923 & (z_w(i,j,k)-z_w(i,j,k-1))
1924 expatt=exp(-att)
1925 itop=par
1926 par=itop*(1.0_r8-expatt)/att
1927 light(i,k)=par
1928
1929
1930
1931
1932 par=itop*expatt
1933 END DO
1934 ELSE
1936 light(i,k)=0.0_r8
1937 END DO
1938 END IF
1939 END DO
1940
1941
1942
1943
1944
1945
1946#ifdef IRON_LIMIT
1947
1948
1949
1950
1951
1952
1953
1954
1955
1956
1957#endif
1958
1963 DO i=istr,iend
1964#ifdef IRON_LIMIT
1965 fnratio=bio(i,k,
ifphy)/max(minval,bio(i,k,
iphyt))
1966 fcratio=fnratio*fen2fec
1968 flimit=fcratio*fcratio/ &
1970
1972 fnlim=min(1.0_r8,flimit/(bio(i,k,
ino3_)*nlimit))
1973#endif
1974 cff4=1.0_r8/sqrt(cff2+cff3*light(i,k)*light(i,k))
1975 cff=bio(i,k,
iphyt)* &
1976#ifdef IRON_LIMIT
1977 & cff1*cff4*light(i,k)*fnlim*nlimit
1978#else
1979 & cff1*cff4*light(i,k)/ &
1981#endif
1984 & bio(i,k,
ino3_)*cff
1985
1986#ifdef IRON_LIMIT
1987
1988
1989
1990 fac=cff*bio(i,k,
ino3_)*fnratio/ &
1991 & max(minval,bio(i,k,
ifdis))
1994 & bio(i,k,
ifdis)*fac
1995
1996
1997
1998 cff5=dtdays*(fcratioe-fcratio)/
t_fe(ng)
1999 cff6=bio(i,k,
iphyt)*cff5*fec2fen
2000 IF (cff6.ge.0.0_r8) THEN
2001 cff=cff6/max(minval,bio(i,k,
ifdis))
2004 & bio(i,k,
ifdis)*cff
2005 ELSE
2006 cff=-cff6/max(minval,bio(i,k,
ifphy))
2009 & bio(i,k,
ifphy)*cff
2010 END IF
2011#endif
2012 END DO
2013 END DO
2014
2015
2016
2017
2018
2019#ifdef IRON_LIMIT
2020
2021
2022#endif
2023
2024 cff1=dtdays*
zoogr(ng)
2027 DO i=istr,iend
2028 cff=bio(i,k,
izoop)* &
2029 & cff1*(1.0_r8-exp(-
ivlev(ng)*bio(i,k,
iphyt)))/ &
2033 & bio(i,k,
iphyt)*cff2*cff
2038#ifdef IRON_LIMIT
2042#endif
2043 END DO
2044 END DO
2045
2046
2047
2048
2051 cff1=1.0_r8/(1.0_r8+cff2+cff3)
2053 DO i=istr,iend
2056 & bio(i,k,
iphyt)*cff2
2058 & bio(i,k,
iphyt)*cff3
2059#ifdef IRON_LIMIT
2063#endif
2064 END DO
2065 END DO
2066
2067
2068
2069
2072 cff1=1.0_r8/(1.0_r8+cff2+cff3)
2074 DO i=istr,iend
2077 & bio(i,k,
izoop)*cff2
2079 & bio(i,k,
izoop)*cff3
2080 END DO
2081 END DO
2082
2083
2084
2085 cff2=dtdays*
detrr(ng)
2086 cff1=1.0_r8/(1.0_r8+cff2)
2088 DO i=istr,iend
2091 & bio(i,k,
isdet)*cff2
2092 END DO
2093 END DO
2094
2095 IF (iteradj.ne.iter) THEN
2096
2097
2098
2099
2100
2101
2102
2103
2104
2105 DO isink=1,nsink
2106 ibio=idsink(isink)
2107
2108
2109
2110
2111
2112
2114 DO i=istr,iend
2115 qc(i,k)=bio(i,k,ibio)
2116 END DO
2117 END DO
2118
2120 DO i=istr,iend
2121 fc(i,k)=(qc(i,k+1)-qc(i,k))*hz_inv2(i,k)
2122 END DO
2123 END DO
2125 DO i=istr,iend
2126 dltr=hz(i,j,k)*fc(i,k)
2127 dltl=hz(i,j,k)*fc(i,k-1)
2128 cff=hz(i,j,k-1)+2.0_r8*hz(i,j,k)+hz(i,j,k+1)
2129 cffr=cff*fc(i,k)
2130 cffl=cff*fc(i,k-1)
2131
2132
2133
2134
2135 IF ((dltr*dltl).le.0.0_r8) THEN
2136 dltr=0.0_r8
2137 dltl=0.0_r8
2138 ELSE IF (abs(dltr).gt.abs(cffl)) THEN
2139 dltr=cffl
2140 ELSE IF (abs(dltl).gt.abs(cffr)) THEN
2141 dltl=cffr
2142 END IF
2143
2144
2145
2146
2147
2148
2149
2150
2151
2152
2153 cff=(dltr-dltl)*hz_inv3(i,k)
2154 dltr=dltr-cff*hz(i,j,k+1)
2155 dltl=dltl+cff*hz(i,j,k-1)
2156 br(i,k)=qc(i,k)+dltr
2157 bl(i,k)=qc(i,k)-dltl
2158 wr(i,k)=(2.0_r8*dltr-dltl)**2
2159 wl(i,k)=(dltr-2.0_r8*dltl)**2
2160 END DO
2161 END DO
2162 cff=1.0e-14_r8
2164 DO i=istr,iend
2165 dltl=max(cff,wl(i,k ))
2166 dltr=max(cff,wr(i,k+1))
2167 br(i,k)=(dltr*br(i,k)+dltl*bl(i,k+1))/(dltr+dltl)
2168 bl(i,k+1)=br(i,k)
2169 END DO
2170 END DO
2171 DO i=istr,iend
2173#if defined LINEAR_CONTINUATION
2174 bl(i,
n(ng))=br(i,
n(ng)-1)
2175 br(i,
n(ng))=2.0_r8*qc(i,
n(ng))-bl(i,
n(ng))
2176#elif defined NEUMANN
2177 bl(i,
n(ng))=br(i,
n(ng)-1)
2178 br(i,
n(ng))=1.5*qc(i,
n(ng))-0.5_r8*bl(i,
n(ng))
2179#else
2180 br(i,
n(ng))=qc(i,
n(ng))
2181 bl(i,
n(ng))=qc(i,
n(ng))
2182 br(i,
n(ng)-1)=qc(i,
n(ng))
2183#endif
2184#if defined LINEAR_CONTINUATION
2185 br(i,1)=bl(i,2)
2186 bl(i,1)=2.0_r8*qc(i,1)-br(i,1)
2187#elif defined NEUMANN
2188 br(i,1)=bl(i,2)
2189 bl(i,1)=1.5_r8*qc(i,1)-0.5_r8*br(i,1)
2190#else
2191 bl(i,2)=qc(i,1)
2192 br(i,1)=qc(i,1)
2193 bl(i,1)=qc(i,1)
2194#endif
2195 END DO
2196
2197
2198
2199
2200
2202 DO i=istr,iend
2203 dltr=br(i,k)-qc(i,k)
2204 dltl=qc(i,k)-bl(i,k)
2205 cffr=2.0_r8*dltr
2206 cffl=2.0_r8*dltl
2207 IF ((dltr*dltl).lt.0.0_r8) THEN
2208 dltr=0.0_r8
2209 dltl=0.0_r8
2210 ELSE IF (abs(dltr).gt.abs(cffl)) THEN
2211 dltr=cffl
2212 ELSE IF (abs(dltl).gt.abs(cffr)) THEN
2213 dltl=cffr
2214 END IF
2215 br(i,k)=qc(i,k)+dltr
2216 bl(i,k)=qc(i,k)-dltl
2217 END DO
2218 END DO
2219
2220
2221
2222
2223
2224
2225
2226
2227
2228
2229
2230
2231
2232
2233
2234 cff=dtdays*abs(wbio(isink))
2236 DO i=istr,iend
2237 fc(i,k-1)=0.0_r8
2238 wl(i,k)=z_w(i,j,k-1)+cff
2239 wr(i,k)=hz(i,j,k)*qc(i,k)
2240 ksource(i,k)=k
2241 END DO
2242 END DO
2245 DO i=istr,iend
2246 IF (wl(i,k).gt.z_w(i,j,ks)) THEN
2247 ksource(i,k)=ks+1
2248 fc(i,k-1)=fc(i,k-1)+wr(i,ks)
2249 END IF
2250 END DO
2251 END DO
2252 END DO
2253
2254
2255
2257 DO i=istr,iend
2258 ks=ksource(i,k)
2259 cu=min(1.0_r8,(wl(i,k)-z_w(i,j,ks-1))*hz_inv(i,ks))
2260 fc(i,k-1)=fc(i,k-1)+ &
2261 & hz(i,j,ks)*cu* &
2262 & (bl(i,ks)+ &
2263 & cu*(0.5_r8*(br(i,ks)-bl(i,ks))- &
2264 & (1.5_r8-cu)* &
2265 & (br(i,ks)+bl(i,ks)- &
2266 & 2.0_r8*qc(i,ks))))
2267 END DO
2268 END DO
2270 DO i=istr,iend
2271 bio(i,k,ibio)=qc(i,k)+ &
2272 & (fc(i,k)-fc(i,k-1))*hz_inv(i,k)
2273 END DO
2274 END DO
2275 END DO
2276 END IF
2277 END DO
2278
2279
2280
2281
2282
2283
2284
2285
2286
2287
2288
2289 sink_loop: DO isink=1,nsink
2290 ibio=idsink(isink)
2291
2292
2293
2294
2295
2296
2298 DO i=istr,iend
2299 qc(i,k)=bio(i,k,ibio)
2300 tl_qc(i,k)=tl_bio(i,k,ibio)
2301 END DO
2302 END DO
2303
2305 DO i=istr,iend
2306 fc(i,k)=(qc(i,k+1)-qc(i,k))*hz_inv2(i,k)
2307 tl_fc(i,k)=(tl_qc(i,k+1)-tl_qc(i,k))*hz_inv2(i,k)+ &
2308 & (qc(i,k+1)-qc(i,k))*tl_hz_inv2(i,k)
2309 END DO
2310 END DO
2312 DO i=istr,iend
2313 dltr=hz(i,j,k)*fc(i,k)
2314 tl_dltr=tl_hz(i,j,k)*fc(i,k)+hz(i,j,k)*tl_fc(i,k)
2315 dltl=hz(i,j,k)*fc(i,k-1)
2316 tl_dltl=tl_hz(i,j,k)*fc(i,k-1)+hz(i,j,k)*tl_fc(i,k-1)
2317 cff=hz(i,j,k-1)+2.0_r8*hz(i,j,k)+hz(i,j,k+1)
2318 tl_cff=tl_hz(i,j,k-1)+2.0_r8*tl_hz(i,j,k)+tl_hz(i,j,k+1)
2319 cffr=cff*fc(i,k)
2320 tl_cffr=tl_cff*fc(i,k)+cff*tl_fc(i,k)
2321 cffl=cff*fc(i,k-1)
2322 tl_cffl=tl_cff*fc(i,k-1)+cff*tl_fc(i,k-1)
2323
2324
2325
2326
2327 IF ((dltr*dltl).le.0.0_r8) THEN
2328 dltr=0.0_r8
2329 tl_dltr=0.0_r8
2330 dltl=0.0_r8
2331 tl_dltl=0.0_r8
2332 ELSE IF (abs(dltr).gt.abs(cffl)) THEN
2333 dltr=cffl
2334 tl_dltr=tl_cffl
2335 ELSE IF (abs(dltl).gt.abs(cffr)) THEN
2336 dltl=cffr
2337 tl_dltl=tl_cffr
2338 END IF
2339
2340
2341
2342
2343
2344
2345
2346
2347
2348
2349 cff=(dltr-dltl)*hz_inv3(i,k)
2350 tl_cff=(tl_dltr-tl_dltl)*hz_inv3(i,k)+ &
2351 & (dltr-dltl)*tl_hz_inv3(i,k)
2352 dltr=dltr-cff*hz(i,j,k+1)
2353 tl_dltr=tl_dltr-tl_cff*hz(i,j,k+1)-cff*tl_hz(i,j,k+1)
2354 dltl=dltl+cff*hz(i,j,k-1)
2355 tl_dltl=tl_dltl+tl_cff*hz(i,j,k-1)+cff*tl_hz(i,j,k-1)
2356 br(i,k)=qc(i,k)+dltr
2357 tl_br(i,k)=tl_qc(i,k)+tl_dltr
2358 bl(i,k)=qc(i,k)-dltl
2359 tl_bl(i,k)=tl_qc(i,k)-tl_dltl
2360 wr(i,k)=(2.0_r8*dltr-dltl)**2
2361 tl_wr(i,k)=2.0_r8*(2.0_r8*dltr-dltl)* &
2362 & (2.0_r8*tl_dltr-tl_dltl)
2363 wl(i,k)=(dltr-2.0_r8*dltl)**2
2364 tl_wl(i,k)=2.0_r8*(dltr-2.0_r8*dltl)* &
2365 & (tl_dltr-2.0_r8*tl_dltl)
2366 END DO
2367 END DO
2368 cff=1.0e-14_r8
2370 DO i=istr,iend
2371 dltl=max(cff,wl(i,k ))
2372 tl_dltl=(0.5_r8-sign(0.5_r8,cff-wl(i,k )))* &
2373 & tl_wl(i,k )
2374 dltr=max(cff,wr(i,k+1))
2375 tl_dltr=(0.5_r8-sign(0.5_r8,cff-wr(i,k+1)))* &
2376 & tl_wr(i,k+1)
2377 br1(i,k)=br(i,k)
2378 bl1(i,k+1)=bl(i,k+1)
2379 br(i,k)=(dltr*br(i,k)+dltl*bl(i,k+1))/(dltr+dltl)
2380 tl_br(i,k)=(tl_dltr*br1(i,k )+dltr*tl_br(i,k )+ &
2381 & tl_dltl*bl1(i,k+1)+dltl*tl_bl(i,k+1))/ &
2382 & (dltr+dltl)- &
2383 & (tl_dltr+tl_dltl)*br(i,k)/(dltr+dltl)
2384 bl(i,k+1)=br(i,k)
2385 tl_bl(i,k+1)=tl_br(i,k)
2386 END DO
2387 END DO
2388 DO i=istr,iend
2390 tl_fc(i,
n(ng))=0.0_r8
2391#if defined LINEAR_CONTINUATION
2392 bl(i,
n(ng))=br(i,
n(ng)-1)
2393 tl_bl(i,
n(ng))=tl_br(i,
n(ng)-1)
2394 br(i,
n(ng))=2.0_r8*qc(i,
n(ng))-bl(i,
n(ng))
2395 tl_br(i,
n(ng))=2.0_r8*tl_qc(i,
n(ng))-tl_bl(i,
n(ng))
2396#elif defined NEUMANN
2397 bl(i,
n(ng))=br(i,
n(ng)-1)
2398 tl_bl(i,
n(ng))=tl_br(i,
n(ng)-1)
2399 br(i,
n(ng))=1.5_r8*qc(i,
n(ng))-0.5_r8*bl(i,
n(ng))
2400 tl_br(i,
n(ng))=1.5_r8*tl_qc(i,
n(ng))-0.5_r8*tl_bl(i,
n(ng))
2401#else
2402 br(i,
n(ng))=qc(i,
n(ng))
2403 bl(i,
n(ng))=qc(i,
n(ng))
2404 br(i,
n(ng)-1)=qc(i,
n(ng))
2405 tl_br(i,
n(ng))=tl_qc(i,
n(ng))
2406 tl_bl(i,
n(ng))=tl_qc(i,
n(ng))
2407 tl_br(i,
n(ng)-1)=tl_qc(i,
n(ng))
2408#endif
2409#if defined LINEAR_CONTINUATION
2410 br(i,1)=bl(i,2)
2411 tl_br(i,1)=tl_bl(i,2)
2412 bl(i,1)=2.0_r8*qc(i,1)-br(i,1)
2413 tl_bl(i,1)=2.0_r8*tl_qc(i,1)-tl_br(i,1)
2414#elif defined NEUMANN
2415 br(i,1)=bl(i,2)
2416 tl_br(i,1)=tl_bl(i,2)
2417 bl(i,1)=1.5_r8*qc(i,1)-0.5_r8*br(i,1)
2418 tl_bl(i,1)=1.5_r8*tl_qc(i,1)-0.5_r8*tl_br(i,1)
2419#else
2420 bl(i,2)=qc(i,1)
2421 br(i,1)=qc(i,1)
2422 bl(i,1)=qc(i,1)
2423 tl_bl(i,2)=tl_qc(i,1)
2424 tl_br(i,1)=tl_qc(i,1)
2425 tl_bl(i,1)=tl_qc(i,1)
2426#endif
2427 END DO
2428
2429
2430
2431
2432
2434 DO i=istr,iend
2435 dltr=br(i,k)-qc(i,k)
2436 tl_dltr=tl_br(i,k)-tl_qc(i,k)
2437 dltl=qc(i,k)-bl(i,k)
2438 tl_dltl=tl_qc(i,k)-tl_bl(i,k)
2439 cffr=2.0_r8*dltr
2440 tl_cffr=2.0_r8*tl_dltr
2441 cffl=2.0_r8*dltl
2442 tl_cffl=2.0_r8*tl_dltl
2443 IF ((dltr*dltl).lt.0.0_r8) THEN
2444 dltr=0.0_r8
2445 tl_dltr=0.0_r8
2446 dltl=0.0_r8
2447 tl_dltl=0.0_r8
2448 ELSE IF (abs(dltr).gt.abs(cffl)) THEN
2449 dltr=cffl
2450 tl_dltr=tl_cffl
2451 ELSE IF (abs(dltl).gt.abs(cffr)) THEN
2452 dltl=cffr
2453 tl_dltl=tl_cffr
2454 END IF
2455 br(i,k)=qc(i,k)+dltr
2456 tl_br(i,k)=tl_qc(i,k)+tl_dltr
2457 bl(i,k)=qc(i,k)-dltl
2458 tl_bl(i,k)=tl_qc(i,k)-tl_dltl
2459 END DO
2460 END DO
2461
2462
2463
2464
2465
2466
2467
2468
2469
2470
2471
2472
2473
2474
2475
2476 cff=dtdays*abs(wbio(isink))
2477 tl_cff=dtdays*sign(1.0_r8,wbio(isink))*tl_wbio(isink)
2479 DO i=istr,iend
2480 fc(i,k-1)=0.0_r8
2481 tl_fc(i,k-1)=0.0_r8
2482 wl(i,k)=z_w(i,j,k-1)+cff
2483 tl_wl(i,k)=tl_z_w(i,j,k-1)+tl_cff
2484 wr(i,k)=hz(i,j,k)*qc(i,k)
2485 tl_wr(i,k)=tl_hz(i,j,k)*qc(i,k)+hz(i,j,k)*tl_qc(i,k)
2486 ksource(i,k)=k
2487 END DO
2488 END DO
2491 DO i=istr,iend
2492 IF (wl(i,k).gt.z_w(i,j,ks)) THEN
2493 ksource(i,k)=ks+1
2494 fc(i,k-1)=fc(i,k-1)+wr(i,ks)
2495 tl_fc(i,k-1)=tl_fc(i,k-1)+tl_wr(i,ks)
2496 END IF
2497 END DO
2498 END DO
2499 END DO
2500
2501
2502
2504 DO i=istr,iend
2505 ks=ksource(i,k)
2506 cu=min(1.0_r8,(wl(i,k)-z_w(i,j,ks-1))*hz_inv(i,ks))
2507 tl_cu=(0.5_r8+sign(0.5_r8, &
2508 & (1.0_r8-(wl(i,k)-z_w(i,j,ks-1))* &
2509 & hz_inv(i,ks))))* &
2510 & ((tl_wl(i,k)-tl_z_w(i,j,ks-1))*hz_inv(i,ks)+ &
2511 & (wl(i,k)-z_w(i,j,ks-1))*tl_hz_inv(i,ks))
2512 fc(i,k-1)=fc(i,k-1)+ &
2513 & hz(i,j,ks)*cu* &
2514 & (bl(i,ks)+ &
2515 & cu*(0.5_r8*(br(i,ks)-bl(i,ks))- &
2516 & (1.5_r8-cu)* &
2517 & (br(i,ks)+bl(i,ks)- &
2518 & 2.0_r8*qc(i,ks))))
2519 tl_fc(i,k-1)=tl_fc(i,k-1)+ &
2520 & (tl_hz(i,j,ks)*cu+hz(i,j,ks)*tl_cu)* &
2521 & (bl(i,ks)+ &
2522 & cu*(0.5_r8*(br(i,ks)-bl(i,ks))- &
2523 & (1.5_r8-cu)* &
2524 & (br(i,ks)+bl(i,ks)- &
2525 & 2.0_r8*qc(i,ks))))+ &
2526 & hz(i,j,ks)*cu* &
2527 & (tl_bl(i,ks)+ &
2528 & tl_cu*(0.5_r8*(br(i,ks)-bl(i,ks))- &
2529 & (1.5_r8-cu)* &
2530 & (br(i,ks)+bl(i,ks)- &
2531 & 2.0_r8*qc(i,ks)))+ &
2532 & cu*(0.5_r8*(tl_br(i,ks)-tl_bl(i,ks))+ &
2533 & tl_cu* &
2534 & (br(i,ks)+bl(i,ks)-2.0_r8*qc(i,ks))- &
2535 & (1.5_r8-cu)* &
2536 & (tl_br(i,ks)+tl_bl(i,ks)- &
2537 & 2.0_r8*tl_qc(i,ks))))
2538 END DO
2539 END DO
2541 DO i=istr,iend
2542 bio(i,k,ibio)=qc(i,k)+(fc(i,k)-fc(i,k-1))*hz_inv(i,k)
2543 tl_bio(i,k,ibio)=tl_qc(i,k)+ &
2544 & (tl_fc(i,k)-tl_fc(i,k-1))*hz_inv(i,k)+ &
2545 & (fc(i,k)-fc(i,k-1))*tl_hz_inv(i,k)
2546 END DO
2547 END DO
2548
2549 END DO sink_loop
2550 END DO iter_loop
2551
2552
2553
2554
2555
2556
2557
2558
2559
2560
2561
2562
2563
2564
2565
2566
2570 DO i=istr,iend
2571 cff=bio(i,k,ibio)-bio_old(i,k,ibio)
2572 tl_cff=tl_bio(i,k,ibio)-tl_bio_old(i,k,ibio)
2573
2574
2575 tl_t(i,j,k,nnew,ibio)=tl_t(i,j,k,nnew,ibio)+ &
2576 & tl_cff*hz(i,j,k)+cff*tl_hz(i,j,k)
2577 END DO
2578 END DO
2579 END DO
2580
2581 END DO j_loop
2582
2583 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