110
111
116
117
118
119 integer, intent(in) :: ng, tile
120 integer, intent(in) :: LBi, UBi, LBj, UBj, UBk, UBt
121 integer, intent(in) :: IminS, ImaxS, JminS, JmaxS
122 integer, intent(in) :: nstp, nnew
123
124#ifdef ASSUMED_SHAPE
125# ifdef MASKING
126 real(r8), intent(in) :: rmask(LBi:,LBj:)
127# endif
128#if defined IRON_LIMIT && defined IRON_RELAX
129 real(r8), intent(in) :: h(LBi:,LBj:)
130# endif
131 real(r8), intent(in) :: Hz(LBi:,LBj:,:)
132 real(r8), intent(in) :: z_r(LBi:,LBj:,:)
133 real(r8), intent(in) :: z_w(LBi:,LBj:,0:)
134 real(r8), intent(in) :: srflx(LBi:,LBj:)
135 real(r8), intent(inout) :: t(LBi:,LBj:,:,:,:)
136
137 real(r8), intent(inout) :: ad_Hz(LBi:,LBj:,:)
138 real(r8), intent(in) :: ad_z_r(LBi:,LBj:,:)
139 real(r8), intent(inout) :: ad_z_w(LBi:,LBj:,0:)
140 real(r8), intent(inout) :: ad_srflx(LBi:,LBj:)
141 real(r8), intent(inout) :: ad_t(LBi:,LBj:,:,:,:)
142#else
143# ifdef MASKING
144 real(r8), intent(in) :: rmask(LBi:UBi,LBj:UBj)
145# endif
146#if defined IRON_LIMIT && defined IRON_RELAX
147 real(r8), intent(in) :: h(LBi:UBi,LBj:UBj)
148# endif
149 real(r8), intent(in) :: Hz(LBi:UBi,LBj:UBj,UBk)
150 real(r8), intent(in) :: z_r(LBi:UBi,LBj:UBj,UBk)
151 real(r8), intent(in) :: z_w(LBi:UBi,LBj:UBj,0:UBk)
152 real(r8), intent(in) :: srflx(LBi:UBi,LBj:UBj)
153 real(r8), intent(inout) :: t(LBi:UBi,LBj:UBj,UBk,3,UBt)
154
155 real(r8), intent(inout) :: ad_Hz(LBi:UBi,LBj:UBj,UBk)
156 real(r8), intent(in) :: ad_z_r(LBi:UBi,LBj:UBj,UBk)
157 real(r8), intent(inout) :: ad_z_w(LBi:UBi,LBj:UBj,0:UBk)
158 real(r8), intent(inout) :: ad_srflx(LBi:UBi,LBj:UBj)
159 real(r8), intent(inout) :: ad_t(LBi:UBi,LBj:UBj,UBk,3,UBt)
160#endif
161
162
163
164 integer, parameter :: Nsink = 2
165
166 integer :: Iter, i, ibio, isink, itime, itrc, iTrcMax, j, k, ks
167 integer :: Iteradj, kk
168
169 integer, dimension(Nsink) :: idsink
170
171 real(r8), parameter :: MinVal = 1.0e-6_r8
172
173 real(r8) :: Att, ExpAtt, Itop, PAR, PAR1
174 real(r8) :: ad_Att, ad_ExpAtt, ad_Itop, ad_PAR
175 real(r8) :: cff, cff1, cff2, cff3, cff4, cff5, cff6, dtdays
176 real(r8) :: ad_cff, ad_cff1, ad_cff4, ad_cff5, ad_cff6
177 real(r8) :: cffL, cffR, cu, dltL, dltR
178 real(r8) :: ad_cffL, ad_cffR, ad_cu, ad_dltL, ad_dltR
179 real(r8) :: fac, fac1, fac2
180 real(r8) :: ad_fac, ad_fac1, ad_fac2
181 real(r8) :: adfac, adfac1, adfac2, adfac3
182#ifdef IRON_LIMIT
183 real(r8) :: Nlimit, FNlim
184 real(r8) :: ad_Nlimit, ad_FNlim
185 real(r8) :: FNratio, FCratio, FCratioE, Flimit
186 real(r8) :: ad_FNratio, ad_FCratio, ad_FCratioE, ad_Flimit
187 real(r8) :: FeC2FeN, FeN2FeC
188# ifdef IRON_RELAX
189 real(r8) :: FeNudgCoef
190# endif
191#endif
192 real(r8), dimension(Nsink) :: Wbio
193 real(r8), dimension(Nsink) :: ad_Wbio
194
195 integer, dimension(IminS:ImaxS,N(ng)) :: ksource
196
197 real(r8), dimension(IminS:ImaxS) :: PARsur
198 real(r8), dimension(IminS:ImaxS) :: ad_PARsur
199
200 real(r8), dimension(NT(ng),2) :: BioTrc
201 real(r8), dimension(NT(ng),2) :: BioTrc1
202 real(r8), dimension(NT(ng),2) :: ad_BioTrc
203 real(r8), dimension(IminS:ImaxS,N(ng),NT(ng)) :: Bio
204 real(r8), dimension(IminS:ImaxS,N(ng),NT(ng)) :: Bio1
205 real(r8), dimension(IminS:ImaxS,N(ng),NT(ng)) :: Bio2
206 real(r8), dimension(IminS:ImaxS,N(ng),NT(ng)) :: Bio_old
207
208 real(r8), dimension(IminS:ImaxS,N(ng),NT(ng)) :: ad_Bio
209 real(r8), dimension(IminS:ImaxS,N(ng),NT(ng)) :: ad_Bio_old
210
211 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: FC
212 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: ad_FC
213
214 real(r8), dimension(IminS:ImaxS,N(ng)) :: Hz_inv
215 real(r8), dimension(IminS:ImaxS,N(ng)) :: Hz_inv2
216 real(r8), dimension(IminS:ImaxS,N(ng)) :: Hz_inv3
217 real(r8), dimension(IminS:ImaxS,N(ng)) :: Light
218 real(r8), dimension(IminS:ImaxS,N(ng)) :: WL
219 real(r8), dimension(IminS:ImaxS,N(ng)) :: WR
220 real(r8), dimension(IminS:ImaxS,N(ng)) :: bL
221 real(r8), dimension(IminS:ImaxS,N(ng)) :: bL1
222 real(r8), dimension(IminS:ImaxS,N(ng)) :: bR
223 real(r8), dimension(IminS:ImaxS,N(ng)) :: bR1
224 real(r8), dimension(IminS:ImaxS,N(ng)) :: qc
225
226 real(r8), dimension(IminS:ImaxS,N(ng)) :: ad_Hz_inv
227 real(r8), dimension(IminS:ImaxS,N(ng)) :: ad_Hz_inv2
228 real(r8), dimension(IminS:ImaxS,N(ng)) :: ad_Hz_inv3
229 real(r8), dimension(IminS:ImaxS,N(ng)) :: ad_Light
230 real(r8), dimension(IminS:ImaxS,N(ng)) :: ad_WL
231 real(r8), dimension(IminS:ImaxS,N(ng)) :: ad_WR
232 real(r8), dimension(IminS:ImaxS,N(ng)) :: ad_bL
233 real(r8), dimension(IminS:ImaxS,N(ng)) :: ad_bR
234 real(r8), dimension(IminS:ImaxS,N(ng)) :: ad_qc
235
236#include "set_bounds.h"
237
238
239
240
241
242
243
245
246
247
249
250#if defined IRON_LIMIT && defined IRON_RELAX
251
252
253
255#endif
256#ifdef IRON_LIMIT
257
258
259
260 fen2fec=(16.0_r8/106.0_r8)*1.0e3_r8
261 fec2fen=(106.0_r8/16.0_r8)*1.0e-3_r8
262#endif
263
264
265
268
269
270
271
274
275 ad_wbio(1)=0.0_r8
276 ad_wbio(2)=0.0_r8
277
278 j_loop : DO j=jstr,jend
279
280
281
282
283
284 ad_par=0.0_r8
285 ad_att=0.0_r8
286 ad_expatt=0.0_r8
287 ad_itop=0.0_r8
288 ad_dltl=0.0_r8
289 ad_dltr=0.0_r8
290 ad_cu=0.0_r8
291 ad_cff=0.0_r8
292 ad_cff1=0.0_r8
293 ad_cff4=0.0_r8
294 ad_cff5=0.0_r8
295 ad_cff6=0.0_r8
296 ad_fac=0.0_r8
297 ad_fac1=0.0_r8
298 ad_fac2=0.0_r8
299 ad_cffl=0.0_r8
300 ad_cffr=0.0_r8
301 adfac=0.0_r8
302 adfac1=0.0_r8
303 adfac2=0.0_r8
304 adfac3=0.0_r8
305#ifdef IRON_LIMIT
306 ad_fnratio=0.0_r8
307 ad_fcratio=0.0_r8
308 ad_fcratioe=0.0_r8
309 ad_flimit=0.0_r8
310 ad_fnlim=0.0_r8
311 ad_nlimit=0.0_r8
312#endif
313
315 DO i=imins,imaxs
316 ad_hz_inv(i,k)=0.0_r8
317 ad_hz_inv2(i,k)=0.0_r8
318 ad_hz_inv3(i,k)=0.0_r8
319 ad_wl(i,k)=0.0_r8
320 ad_wr(i,k)=0.0_r8
321 ad_bl(i,k)=0.0_r8
322 ad_br(i,k)=0.0_r8
323 ad_qc(i,k)=0.0_r8
324 ad_light(i,k)=0.0_r8
325 END DO
326 END DO
329 ad_biotrc(ibio,1)=0.0_r8
330 ad_biotrc(ibio,2)=0.0_r8
331 END DO
332 DO i=imins,imaxs
333 ad_parsur(i)=0.0_r8
334 END DO
336 DO i=imins,imaxs
337 ad_fc(i,k)=0.0_r8
338 END DO
339 END DO
340
341
342
346 DO i=istr,iend
347 bio(i,k,ibio)=0.0_r8
348 bio1(i,k,ibio)=0.0_r8
349 bio2(i,k,ibio)=0.0_r8
350 ad_bio(i,k,ibio)=0.0_r8
351 ad_bio_old(i,k,ibio)=0.0_r8
352 END DO
353 END DO
354 END DO
355
356
357
359 DO i=istr,iend
360 hz_inv(i,k)=1.0_r8/hz(i,j,k)
361 END DO
362 END DO
364 DO i=istr,iend
365 hz_inv2(i,k)=1.0_r8/(hz(i,j,k)+hz(i,j,k+1))
366 END DO
367 END DO
369 DO i=istr,iend
370 hz_inv3(i,k)=1.0_r8/(hz(i,j,k-1)+hz(i,j,k)+hz(i,j,k+1))
371 END DO
372 END DO
373
374
375
376
377
378
379
380
382 DO i=istr,iend
383
384
385
386
387
388
389
390
391
392
393
394
397
398
399 biotrc(ibio,nstp)=t(i,j,k,nstp,ibio)
400
401
402 biotrc(ibio,nnew)=t(i,j,k,nnew,ibio)
403 END DO
404
405
406
407 cff2=0.0_r8
408 DO itime=1,2
409 cff1=0.0_r8
411#ifdef IRON_LIMIT
413#else
415#endif
417 cff1=cff1+max(0.0_r8,minval-biotrc(ibio,itime))
418 IF (biotrc(ibio,itime).gt.biotrc(itrcmax,itime)) THEN
419 itrcmax=ibio
420 END IF
421 biotrc(ibio,itime)=max(minval,biotrc(ibio,itime))
422 END DO
423 IF (biotrc(itrcmax,itime).gt.cff1) THEN
424 biotrc(itrcmax,itime)=biotrc(itrcmax,itime)-cff1
425 END IF
426#ifdef IRON_LIMIT
429 biotrc(ibio,itime)=max(minval,biotrc(ibio,itime))
430 END DO
431#endif
432 END DO
433
434
435
438 bio_old(i,k,ibio)=biotrc(ibio,nstp)
439 bio(i,k,ibio)=biotrc(ibio,nstp)
440 END DO
441
442#if defined IRON_LIMIT && defined IRON_RELAX
443
444
445
446
447
448 IF (h(i,j).le.
fehmin(ng))
THEN
451 END IF
452#endif
453 END DO
454 END DO
455
456
457
458
459
460 DO i=istr,iend
461#ifdef CONST_PAR
462
463
464
465 parsur(i)=158.075_r8
466#else
468#endif
469 END DO
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 iter_loop:
DO iter=1,
bioiter(ng)
526
527
528
529 DO i=istr,iend
530 par=parsur(i)
531 IF (parsur(i).gt.0.0_r8) THEN
533
534
535
536
537
539 & (z_w(i,j,k)-z_w(i,j,k-1))
540 expatt=exp(-att)
541 itop=par
542 par=itop*(1.0_r8-expatt)/att
543 light(i,k)=par
544
545
546
547
548 par=itop*expatt
549 END DO
550 ELSE
552 light(i,k)=0.0_r8
553 END DO
554 END IF
555 END DO
556
557
558
559
560
561
562#ifdef IRON_LIMIT
563
564
565
566
567
568
569
570
571
572
573#endif
574
579 DO i=istr,iend
580#ifdef IRON_LIMIT
581 fnratio=bio(i,k,
ifphy)/max(minval,bio(i,k,
iphyt))
582 fcratio=fnratio*fen2fec
584 flimit=fcratio*fcratio/ &
586
588 fnlim=min(1.0_r8,flimit/(bio(i,k,
ino3_)*nlimit))
589#endif
590 cff4=1.0_r8/sqrt(cff2+cff3*light(i,k)*light(i,k))
591 cff=bio(i,k,
iphyt)* &
592#ifdef IRON_LIMIT
593 & cff1*cff4*light(i,k)*fnlim*nlimit
594#else
595
596 & cff1*cff4*light(i,k)/ &
598#endif
599
603#ifdef IRON_LIMIT
604
605
606
607 fac=cff*bio(i,k,
ino3_)*fnratio/ &
608 & max(minval,bio(i,k,
ifdis))
612
613
614
615 cff5=dtdays*(fcratioe-fcratio)/
t_fe(ng)
616 cff6=bio(i,k,
iphyt)*cff5*fec2fen
617 IF (cff6.ge.0.0_r8) THEN
618 cff=cff6/max(minval,bio(i,k,
ifdis))
622 ELSE
623 cff=-cff6/max(minval,bio(i,k,
ifphy))
627 END IF
628#endif
629 END DO
630 END DO
631
632
633
634
635
636#ifdef IRON_LIMIT
637
638
639#endif
640
641 cff1=dtdays*
zoogr(ng)
644 DO i=istr,iend
645 cff=bio(i,k,
izoop)* &
646 & cff1*(1.0_r8-exp(-
ivlev(ng)*bio(i,k,
iphyt)))/ &
650 & bio(i,k,
iphyt)*cff2*cff
655#ifdef IRON_LIMIT
659#endif
660 END DO
661 END DO
662
663
664
665
668 cff1=1.0_r8/(1.0_r8+cff2+cff3)
670 DO i=istr,iend
673 & bio(i,k,
iphyt)*cff2
675 & bio(i,k,
iphyt)*cff3
676#ifdef IRON_LIMIT
680#endif
681 END DO
682 END DO
683
684
685
686
689 cff1=1.0_r8/(1.0_r8+cff2+cff3)
691 DO i=istr,iend
694 & bio(i,k,
izoop)*cff2
696 & bio(i,k,
izoop)*cff3
697 END DO
698 END DO
699
700
701
702 cff2=dtdays*
detrr(ng)
703 cff1=1.0_r8/(1.0_r8+cff2)
705 DO i=istr,iend
708 & bio(i,k,
isdet)*cff2
709 END DO
710 END DO
711
712
713
714
715
716
717
718
719
720 sink_loop: DO isink=1,nsink
721 ibio=idsink(isink)
722
723
724
725
726
727
729 DO i=istr,iend
730 qc(i,k)=bio(i,k,ibio)
731 END DO
732 END DO
733
735 DO i=istr,iend
736 fc(i,k)=(qc(i,k+1)-qc(i,k))*hz_inv2(i,k)
737 END DO
738 END DO
740 DO i=istr,iend
741 dltr=hz(i,j,k)*fc(i,k)
742 dltl=hz(i,j,k)*fc(i,k-1)
743 cff=hz(i,j,k-1)+2.0_r8*hz(i,j,k)+hz(i,j,k+1)
744 cffr=cff*fc(i,k)
745 cffl=cff*fc(i,k-1)
746
747
748
749
750 IF ((dltr*dltl).le.0.0_r8) THEN
751 dltr=0.0_r8
752 dltl=0.0_r8
753 ELSE IF (abs(dltr).gt.abs(cffl)) THEN
754 dltr=cffl
755 ELSE IF (abs(dltl).gt.abs(cffr)) THEN
756 dltl=cffr
757 END IF
758
759
760
761
762
763
764
765
766
767
768 cff=(dltr-dltl)*hz_inv3(i,k)
769 dltr=dltr-cff*hz(i,j,k+1)
770 dltl=dltl+cff*hz(i,j,k-1)
771 br(i,k)=qc(i,k)+dltr
772 bl(i,k)=qc(i,k)-dltl
773 wr(i,k)=(2.0_r8*dltr-dltl)**2
774 wl(i,k)=(dltr-2.0_r8*dltl)**2
775 END DO
776 END DO
777 cff=1.0e-14_r8
779 DO i=istr,iend
780 dltl=max(cff,wl(i,k ))
781 dltr=max(cff,wr(i,k+1))
782 br(i,k)=(dltr*br(i,k)+dltl*bl(i,k+1))/(dltr+dltl)
783 bl(i,k+1)=br(i,k)
784 END DO
785 END DO
786 DO i=istr,iend
788#if defined LINEAR_CONTINUATION
789 bl(i,
n(ng))=br(i,
n(ng)-1)
790 br(i,
n(ng))=2.0_r8*qc(i,
n(ng))-bl(i,
n(ng))
791#elif defined NEUMANN
792 bl(i,
n(ng))=br(i,
n(ng)-1)
793 br(i,
n(ng))=1.5*qc(i,
n(ng))-0.5_r8*bl(i,
n(ng))
794#else
795 br(i,
n(ng))=qc(i,
n(ng))
796 bl(i,
n(ng))=qc(i,
n(ng))
797 br(i,
n(ng)-1)=qc(i,
n(ng))
798#endif
799#if defined LINEAR_CONTINUATION
800 br(i,1)=bl(i,2)
801 bl(i,1)=2.0_r8*qc(i,1)-br(i,1)
802#elif defined NEUMANN
803 br(i,1)=bl(i,2)
804 bl(i,1)=1.5_r8*qc(i,1)-0.5_r8*br(i,1)
805#else
806 bl(i,2)=qc(i,1)
807 br(i,1)=qc(i,1)
808 bl(i,1)=qc(i,1)
809#endif
810 END DO
811
812
813
814
815
817 DO i=istr,iend
818 dltr=br(i,k)-qc(i,k)
819 dltl=qc(i,k)-bl(i,k)
820 cffr=2.0_r8*dltr
821 cffl=2.0_r8*dltl
822 IF ((dltr*dltl).lt.0.0_r8) THEN
823 dltr=0.0_r8
824 dltl=0.0_r8
825 ELSE IF (abs(dltr).gt.abs(cffl)) THEN
826 dltr=cffl
827 ELSE IF (abs(dltl).gt.abs(cffr)) THEN
828 dltl=cffr
829 END IF
830 br(i,k)=qc(i,k)+dltr
831 bl(i,k)=qc(i,k)-dltl
832 END DO
833 END DO
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849 cff=dtdays*abs(wbio(isink))
851 DO i=istr,iend
852 fc(i,k-1)=0.0_r8
853 wl(i,k)=z_w(i,j,k-1)+cff
854 wr(i,k)=hz(i,j,k)*qc(i,k)
855 ksource(i,k)=k
856 END DO
857 END DO
860 DO i=istr,iend
861 IF (wl(i,k).gt.z_w(i,j,ks)) THEN
862 ksource(i,k)=ks+1
863 fc(i,k-1)=fc(i,k-1)+wr(i,ks)
864 END IF
865 END DO
866 END DO
867 END DO
868
869
870
872 DO i=istr,iend
873 ks=ksource(i,k)
874 cu=min(1.0_r8,(wl(i,k)-z_w(i,j,ks-1))*hz_inv(i,ks))
875 fc(i,k-1)=fc(i,k-1)+ &
876 & hz(i,j,ks)*cu* &
877 & (bl(i,ks)+ &
878 & cu*(0.5_r8*(br(i,ks)-bl(i,ks))- &
879 & (1.5_r8-cu)* &
880 & (br(i,ks)+bl(i,ks)- &
881 & 2.0_r8*qc(i,ks))))
882 END DO
883 END DO
885 DO i=istr,iend
886 bio(i,k,ibio)=qc(i,k)+(fc(i,k)-fc(i,k-1))*hz_inv(i,k)
887 END DO
888 END DO
889
890 END DO sink_loop
891 END DO iter_loop
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
911 DO i=istr,iend
912 cff=bio(i,k,ibio)-bio_old(i,k,ibio)
913
914
915
916 ad_hz(i,j,k)=ad_hz(i,j,k)+cff*ad_t(i,j,k,nnew,ibio)
917 ad_cff=ad_cff+hz(i,j,k)*ad_t(i,j,k,nnew,ibio)
918
919
920 ad_bio_old(i,k,ibio)=ad_bio_old(i,k,ibio)-ad_cff
921 ad_bio(i,k,ibio)=ad_bio(i,k,ibio)+ad_cff
922 ad_cff=0.0_r8
923 END DO
924 END DO
925 END DO
926
927
928
929
930
931
932 iter_loop1:
DO iter=
bioiter(ng),1,-1
933
934
935
936
937
938
939
940
941
942
943
945 DO i=istr,iend
946
947
948
949
950
951
952
953
954
955
956
957
960
961
962 biotrc(ibio,nstp)=t(i,j,k,nstp,ibio)
963
964
965 biotrc(ibio,nnew)=t(i,j,k,nnew,ibio)
966 END DO
967
968
969
970 cff2=0.0_r8
971 DO itime=1,2
972 cff1=0.0_r8
974#ifdef IRON_LIMIT
976#else
978#endif
980 cff1=cff1+max(0.0_r8,minval-biotrc(ibio,itime))
981 IF (biotrc(ibio,itime).gt.biotrc(itrcmax,itime)) THEN
982 itrcmax=ibio
983 END IF
984 biotrc(ibio,itime)=max(minval,biotrc(ibio,itime))
985 END DO
986 IF (biotrc(itrcmax,itime).gt.cff1) THEN
987 biotrc(itrcmax,itime)=biotrc(itrcmax,itime)-cff1
988 END IF
989#ifdef IRON_LIMIT
992 biotrc(ibio,itime)=max(minval,biotrc(ibio,itime))
993 END DO
994#endif
995 END DO
996
997
998
1001 bio_old(i,k,ibio)=biotrc(ibio,nstp)
1002 bio(i,k,ibio)=biotrc(ibio,nstp)
1003 END DO
1004
1005#if defined IRON_LIMIT && defined IRON_RELAX
1006
1007
1008
1009
1010
1011 IF (h(i,j).le.
fehmin(ng))
THEN
1014 END IF
1015#endif
1016 END DO
1017 END DO
1018
1019
1020
1021
1022
1023 DO i=istr,iend
1024#ifdef CONST_PAR
1025
1026
1027
1028 parsur(i)=158.075_r8
1029#else
1031#endif
1032 END DO
1033
1034
1035
1036
1037
1038
1039 DO iteradj=1,iter
1040
1041
1042
1043 DO i=istr,iend
1044 par=parsur(i)
1045 IF (parsur(i).gt.0.0_r8) THEN
1047
1048
1049
1050
1051
1053 & (z_w(i,j,k)-z_w(i,j,k-1))
1054 expatt=exp(-att)
1055 itop=par
1056 par=itop*(1.0_r8-expatt)/att
1057 light(i,k)=par
1058
1059
1060
1061
1062 par=itop*expatt
1063 END DO
1064 ELSE
1066 light(i,k)=0.0_r8
1067 END DO
1068 END IF
1069 END DO
1070
1071
1072
1073
1074
1075
1076#ifdef IRON_LIMIT
1077
1078
1079
1080
1081
1082
1083
1084
1085
1086
1087#endif
1088
1093 DO i=istr,iend
1094#ifdef IRON_LIMIT
1095
1096
1097
1098 fnratio=bio(i,k,
ifphy)/max(minval,bio(i,k,
iphyt))
1099 fcratio=fnratio*fen2fec
1101 flimit=fcratio*fcratio/ &
1103
1105 fnlim=min(1.0_r8,flimit/(bio(i,k,
ino3_)*nlimit))
1106#endif
1107 cff4=1.0_r8/sqrt(cff2+cff3*light(i,k)*light(i,k))
1108 cff=bio(i,k,
iphyt)* &
1109#ifdef IRON_LIMIT
1110 & cff1*cff4*light(i,k)*fnlim*nlimit
1111#else
1112 & cff1*cff4*light(i,k)/ &
1114#endif
1117 & bio(i,k,
ino3_)*cff
1118#ifdef IRON_LIMIT
1119
1120
1121
1122 fac=cff*bio(i,k,
ino3_)*fnratio/ &
1123 & max(minval,bio(i,k,
ifdis))
1126 & bio(i,k,
ifdis)*fac
1127
1128
1129
1130 cff=cff*bio(i,k,
ino3_)* &
1131 & fnratio/max(minval,bio(i,k,
ifdis))
1134 & bio(i,k,
ifdis)*cff
1135
1136
1137
1138 cff5=dtdays*(fcratioe-fcratio)/
t_fe(ng)
1139 cff6=bio(i,k,
iphyt)*cff5*fec2fen
1140 IF (cff6.ge.0.0_r8) THEN
1141 cff=cff6/max(minval,bio(i,k,
ifdis))
1144 & bio(i,k,
ifdis)*cff
1145 ELSE
1146 cff=-cff6/max(minval,bio(i,k,
ifphy))
1149 & bio(i,k,
ifphy)*cff
1150 END IF
1151#endif
1152 END DO
1153 END DO
1154
1155
1156
1157
1158
1159#ifdef IRON_LIMIT
1160
1161
1162#endif
1163
1164 cff1=dtdays*
zoogr(ng)
1167 DO i=istr,iend
1168 cff=bio(i,k,
izoop)* &
1169 & cff1*(1.0_r8-exp(-
ivlev(ng)*bio(i,k,
iphyt)))/ &
1173 & bio(i,k,
iphyt)*cff2*cff
1178#ifdef IRON_LIMIT
1182#endif
1183 END DO
1184 END DO
1185
1186
1187
1188
1191 cff1=1.0_r8/(1.0_r8+cff2+cff3)
1193 DO i=istr,iend
1196 & bio(i,k,
iphyt)*cff2
1198 & bio(i,k,
iphyt)*cff3
1199#ifdef IRON_LIMIT
1203#endif
1204 END DO
1205 END DO
1206
1207
1208
1209
1212 cff1=1.0_r8/(1.0_r8+cff2+cff3)
1214 DO i=istr,iend
1217 & bio(i,k,
izoop)*cff2
1219 & bio(i,k,
izoop)*cff3
1220 END DO
1221 END DO
1222
1223
1224
1225 cff2=dtdays*
detrr(ng)
1226 cff1=1.0_r8/(1.0_r8+cff2)
1228 DO i=istr,iend
1231 & bio(i,k,
isdet)*cff2
1232 END DO
1233 END DO
1234
1235 IF (iteradj.ne.iter) THEN
1236
1237
1238
1239
1240
1241
1242
1243
1244
1245 DO isink=1,nsink
1246 ibio=idsink(isink)
1247
1248
1249
1250
1251
1252
1254 DO i=istr,iend
1255 qc(i,k)=bio(i,k,ibio)
1256 END DO
1257 END DO
1258
1260 DO i=istr,iend
1261 fc(i,k)=(qc(i,k+1)-qc(i,k))*hz_inv2(i,k)
1262 END DO
1263 END DO
1265 DO i=istr,iend
1266 dltr=hz(i,j,k)*fc(i,k)
1267 dltl=hz(i,j,k)*fc(i,k-1)
1268 cff=hz(i,j,k-1)+2.0_r8*hz(i,j,k)+hz(i,j,k+1)
1269 cffr=cff*fc(i,k)
1270 cffl=cff*fc(i,k-1)
1271
1272
1273
1274
1275 IF ((dltr*dltl).le.0.0_r8) THEN
1276 dltr=0.0_r8
1277 dltl=0.0_r8
1278 ELSE IF (abs(dltr).gt.abs(cffl)) THEN
1279 dltr=cffl
1280 ELSE IF (abs(dltl).gt.abs(cffr)) THEN
1281 dltl=cffr
1282 END IF
1283
1284
1285
1286
1287
1288
1289
1290
1291
1292
1293 cff=(dltr-dltl)*hz_inv3(i,k)
1294 dltr=dltr-cff*hz(i,j,k+1)
1295 dltl=dltl+cff*hz(i,j,k-1)
1296 br(i,k)=qc(i,k)+dltr
1297 bl(i,k)=qc(i,k)-dltl
1298 wr(i,k)=(2.0_r8*dltr-dltl)**2
1299 wl(i,k)=(dltr-2.0_r8*dltl)**2
1300 END DO
1301 END DO
1302 cff=1.0e-14_r8
1304 DO i=istr,iend
1305 dltl=max(cff,wl(i,k ))
1306 dltr=max(cff,wr(i,k+1))
1307 br(i,k)=(dltr*br(i,k)+dltl*bl(i,k+1))/(dltr+dltl)
1308 bl(i,k+1)=br(i,k)
1309 END DO
1310 END DO
1311 DO i=istr,iend
1313#if defined LINEAR_CONTINUATION
1314 bl(i,
n(ng))=br(i,
n(ng)-1)
1315 br(i,
n(ng))=2.0_r8*qc(i,
n(ng))-bl(i,
n(ng))
1316#elif defined NEUMANN
1317 bl(i,
n(ng))=br(i,
n(ng)-1)
1318 br(i,
n(ng))=1.5*qc(i,
n(ng))-0.5_r8*bl(i,
n(ng))
1319#else
1320 br(i,
n(ng))=qc(i,
n(ng))
1321 bl(i,
n(ng))=qc(i,
n(ng))
1322 br(i,
n(ng)-1)=qc(i,
n(ng))
1323#endif
1324#if defined LINEAR_CONTINUATION
1325 br(i,1)=bl(i,2)
1326 bl(i,1)=2.0_r8*qc(i,1)-br(i,1)
1327#elif defined NEUMANN
1328 br(i,1)=bl(i,2)
1329 bl(i,1)=1.5_r8*qc(i,1)-0.5_r8*br(i,1)
1330#else
1331 bl(i,2)=qc(i,1)
1332 br(i,1)=qc(i,1)
1333 bl(i,1)=qc(i,1)
1334#endif
1335 END DO
1336
1337
1338
1339
1340
1342 DO i=istr,iend
1343 dltr=br(i,k)-qc(i,k)
1344 dltl=qc(i,k)-bl(i,k)
1345 cffr=2.0_r8*dltr
1346 cffl=2.0_r8*dltl
1347 IF ((dltr*dltl).lt.0.0_r8) THEN
1348 dltr=0.0_r8
1349 dltl=0.0_r8
1350 ELSE IF (abs(dltr).gt.abs(cffl)) THEN
1351 dltr=cffl
1352 ELSE IF (abs(dltl).gt.abs(cffr)) THEN
1353 dltl=cffr
1354 END IF
1355 br(i,k)=qc(i,k)+dltr
1356 bl(i,k)=qc(i,k)-dltl
1357 END DO
1358 END DO
1359
1360
1361
1362
1363
1364
1365
1366
1367
1368
1369
1370
1371
1372
1373
1374 cff=dtdays*abs(wbio(isink))
1376 DO i=istr,iend
1377 fc(i,k-1)=0.0_r8
1378 wl(i,k)=z_w(i,j,k-1)+cff
1379 wr(i,k)=hz(i,j,k)*qc(i,k)
1380 ksource(i,k)=k
1381 END DO
1382 END DO
1385 DO i=istr,iend
1386 IF (wl(i,k).gt.z_w(i,j,ks)) THEN
1387 ksource(i,k)=ks+1
1388 fc(i,k-1)=fc(i,k-1)+wr(i,ks)
1389 END IF
1390 END DO
1391 END DO
1392 END DO
1393
1394
1395
1397 DO i=istr,iend
1398 ks=ksource(i,k)
1399 cu=min(1.0_r8,(wl(i,k)-z_w(i,j,ks-1))*hz_inv(i,ks))
1400 fc(i,k-1)=fc(i,k-1)+ &
1401 & hz(i,j,ks)*cu* &
1402 & (bl(i,ks)+ &
1403 & cu*(0.5_r8*(br(i,ks)-bl(i,ks))- &
1404 & (1.5_r8-cu)* &
1405 & (br(i,ks)+bl(i,ks)- &
1406 & 2.0_r8*qc(i,ks))))
1407 END DO
1408 END DO
1410 DO i=istr,iend
1411 bio(i,k,ibio)=qc(i,k)+ &
1412 & (fc(i,k)-fc(i,k-1))*hz_inv(i,k)
1413 END DO
1414 END DO
1415 END DO
1416 END IF
1417 END DO
1418
1419
1420
1421 sink_loop1: DO isink=1,nsink
1422 ibio=idsink(isink)
1423
1424
1425
1426
1427
1428
1429
1430
1431
1433 DO i=istr,iend
1434 qc(i,k)=bio(i,k,ibio)
1435 END DO
1436 END DO
1437
1439 DO i=istr,iend
1440 fc(i,k)=(qc(i,k+1)-qc(i,k))*hz_inv2(i,k)
1441 END DO
1442 END DO
1444 DO i=istr,iend
1445 dltr=hz(i,j,k)*fc(i,k)
1446 dltl=hz(i,j,k)*fc(i,k-1)
1447 cff=hz(i,j,k-1)+2.0_r8*hz(i,j,k)+hz(i,j,k+1)
1448 cffr=cff*fc(i,k)
1449 cffl=cff*fc(i,k-1)
1450
1451
1452
1453
1454 IF ((dltr*dltl).le.0.0_r8) THEN
1455 dltr=0.0_r8
1456 dltl=0.0_r8
1457 ELSE IF (abs(dltr).gt.abs(cffl)) THEN
1458 dltr=cffl
1459 ELSE IF (abs(dltl).gt.abs(cffr)) THEN
1460 dltl=cffr
1461 END IF
1462
1463
1464
1465
1466
1467
1468
1469
1470
1471
1472 cff=(dltr-dltl)*hz_inv3(i,k)
1473 dltr=dltr-cff*hz(i,j,k+1)
1474 dltl=dltl+cff*hz(i,j,k-1)
1475 br(i,k)=qc(i,k)+dltr
1476 bl(i,k)=qc(i,k)-dltl
1477 wr(i,k)=(2.0_r8*dltr-dltl)**2
1478 wl(i,k)=(dltr-2.0_r8*dltl)**2
1479 END DO
1480 END DO
1481 cff=1.0e-14_r8
1483 DO i=istr,iend
1484 dltl=max(cff,wl(i,k ))
1485 dltr=max(cff,wr(i,k+1))
1486 br(i,k)=(dltr*br(i,k)+dltl*bl(i,k+1))/(dltr+dltl)
1487 bl(i,k+1)=br(i,k)
1488 END DO
1489 END DO
1490 DO i=istr,iend
1492#if defined LINEAR_CONTINUATION
1493 bl(i,
n(ng))=br(i,
n(ng)-1)
1494 br(i,
n(ng))=2.0_r8*qc(i,
n(ng))-bl(i,
n(ng))
1495#elif defined NEUMANN
1496 bl(i,
n(ng))=br(i,
n(ng)-1)
1497 br(i,
n(ng))=1.5*qc(i,
n(ng))-0.5_r8*bl(i,
n(ng))
1498#else
1499 br(i,
n(ng))=qc(i,
n(ng))
1500 bl(i,
n(ng))=qc(i,
n(ng))
1501 br(i,
n(ng)-1)=qc(i,
n(ng))
1502#endif
1503#if defined LINEAR_CONTINUATION
1504 br(i,1)=bl(i,2)
1505 bl(i,1)=2.0_r8*qc(i,1)-br(i,1)
1506#elif defined NEUMANN
1507 br(i,1)=bl(i,2)
1508 bl(i,1)=1.5_r8*qc(i,1)-0.5_r8*br(i,1)
1509#else
1510 bl(i,2)=qc(i,1)
1511 br(i,1)=qc(i,1)
1512 bl(i,1)=qc(i,1)
1513#endif
1514 END DO
1515
1516
1517
1518
1519
1521 DO i=istr,iend
1522 dltr=br(i,k)-qc(i,k)
1523 dltl=qc(i,k)-bl(i,k)
1524 cffr=2.0_r8*dltr
1525 cffl=2.0_r8*dltl
1526 IF ((dltr*dltl).lt.0.0_r8) THEN
1527 dltr=0.0_r8
1528 dltl=0.0_r8
1529 ELSE IF (abs(dltr).gt.abs(cffl)) THEN
1530 dltr=cffl
1531 ELSE IF (abs(dltl).gt.abs(cffr)) THEN
1532 dltl=cffr
1533 END IF
1534 br(i,k)=qc(i,k)+dltr
1535 bl(i,k)=qc(i,k)-dltl
1536 END DO
1537 END DO
1538
1539
1540
1541
1542
1543
1544
1545
1546
1547
1548
1549
1550
1551
1552
1553 cff=dtdays*abs(wbio(isink))
1555 DO i=istr,iend
1556 fc(i,k-1)=0.0_r8
1557 wl(i,k)=z_w(i,j,k-1)+cff
1558 wr(i,k)=hz(i,j,k)*qc(i,k)
1559 ksource(i,k)=k
1560 END DO
1561 END DO
1564 DO i=istr,iend
1565 IF (wl(i,k).gt.z_w(i,j,ks)) THEN
1566 ksource(i,k)=ks+1
1567 fc(i,k-1)=fc(i,k-1)+wr(i,ks)
1568 END IF
1569 END DO
1570 END DO
1571 END DO
1572
1573
1574
1576 DO i=istr,iend
1577 ks=ksource(i,k)
1578 cu=min(1.0_r8,(wl(i,k)-z_w(i,j,ks-1))*hz_inv(i,ks))
1579 fc(i,k-1)=fc(i,k-1)+ &
1580 & hz(i,j,ks)*cu* &
1581 & (bl(i,ks)+ &
1582 & cu*(0.5_r8*(br(i,ks)-bl(i,ks))- &
1583 & (1.5_r8-cu)* &
1584 & (br(i,ks)+bl(i,ks)- &
1585 & 2.0_r8*qc(i,ks))))
1586 END DO
1587 END DO
1589 DO i=istr,iend
1590
1591
1592
1593
1594 ad_qc(i,k)=ad_qc(i,k)+ad_bio(i,k,ibio)
1595 ad_fc(i,k)=ad_fc(i,k)+hz_inv(i,k)*ad_bio(i,k,ibio)
1596 ad_fc(i,k-1)=ad_fc(i,k-1)-hz_inv(i,k)*ad_bio(i,k,ibio)
1597 ad_hz_inv(i,k)=ad_hz_inv(i,k)+ &
1598 & (fc(i,k)-fc(i,k-1))*ad_bio(i,k,ibio)
1599 ad_bio(i,k,ibio)=0.0_r8
1600 END DO
1601 END DO
1602
1603
1604
1606 DO i=istr,iend
1607 ks=ksource(i,k)
1608 cu=min(1.0_r8,(wl(i,k)-z_w(i,j,ks-1))*hz_inv(i,ks))
1609
1610
1611
1612
1613
1614
1615
1616
1617
1618
1619
1620
1621
1622
1623
1624
1625
1626
1627
1628
1629 adfac=(bl(i,ks)+ &
1630 & cu*(0.5_r8*(br(i,ks)-bl(i,ks))- &
1631 & (1.5_r8-cu)* &
1632 & (br(i,ks)+bl(i,ks)- &
1633 & 2.0_r8*qc(i,ks))))*ad_fc(i,k-1)
1634 adfac1=hz(i,j,ks)*cu*ad_fc(i,k-1)
1635 adfac2=adfac1*cu
1636 adfac3=adfac2*(1.5_r8-cu)
1637 ad_hz(i,j,ks)=ad_hz(i,j,ks)+cu*adfac
1638 ad_cu=ad_cu+hz(i,j,ks)*adfac
1639 ad_bl(i,ks)=ad_bl(i,ks)+adfac1
1640 ad_cu=ad_cu+ &
1641 & adfac1*(0.5_r8*(br(i,ks)-bl(i,ks))- &
1642 & (1.5_r8-cu)* &
1643 & (br(i,ks)+bl(i,ks)- &
1644 & 2.0_r8*qc(i,ks)))+ &
1645 & adfac2*(br(i,ks)+bl(i,ks)-2.0_r8*qc(i,ks))
1646 ad_br(i,ks)=ad_br(i,ks)+0.5_r8*adfac2-adfac3
1647 ad_bl(i,ks)=ad_bl(i,ks)-0.5_r8*adfac2-adfac3
1648 ad_qc(i,ks)=ad_qc(i,ks)+2.0_r8*adfac3
1649
1650
1651
1652
1653
1654
1655 adfac=(0.5_r8+sign(0.5_r8, &
1656 & (1.0_r8-(wl(i,k)-z_w(i,j,ks-1))* &
1657 & hz_inv(i,ks))))*ad_cu
1658 adfac1=adfac*hz_inv(i,ks)
1659 ad_wl(i,k)=ad_wl(i,k)+adfac1
1660 ad_z_w(i,j,ks-1)=ad_z_w(i,j,ks-1)-adfac1
1661 ad_hz_inv(i,ks)=ad_hz_inv(i,ks)+ &
1662 & (wl(i,k)-z_w(i,j,ks-1))*adfac
1663 ad_cu=0.0_r8
1664 END DO
1665 END DO
1666
1667
1668
1669
1670
1671
1672
1673
1674
1675
1676
1677
1678
1679
1680
1681 cff=dtdays*abs(wbio(isink))
1684 DO i=istr,iend
1685 IF (wl(i,k).gt.z_w(i,j,ks)) THEN
1686
1687
1688 ad_wr(i,ks)=ad_wr(i,ks)+ad_fc(i,k-1)
1689 END IF
1690 END DO
1691 END DO
1692 END DO
1694 DO i=istr,iend
1695
1696
1697 ad_hz(i,j,k)=ad_hz(i,j,k)+qc(i,k)*ad_wr(i,k)
1698 ad_qc(i,k)=ad_qc(i,k)+hz(i,j,k)*ad_wr(i,k)
1699 ad_wr(i,k)=0.0_r8
1700
1701
1702 ad_z_w(i,j,k-1)=ad_z_w(i,j,k-1)+ad_wl(i,k)
1703 ad_cff=ad_cff+ad_wl(i,k)
1704 ad_wl(i,k)=0.0_r8
1705
1706
1707 ad_fc(i,k-1)=0.0_r8
1708 END DO
1709 END DO
1710
1711
1712 ad_wbio(isink)=ad_wbio(isink)+ &
1713 & dtdays*sign(1.0_r8,wbio(isink))*ad_cff
1714 ad_cff=0.0_r8
1715
1716
1717
1718
1719
1720
1721
1722
1724 DO i=istr,iend
1725 qc(i,k)=bio(i,k,ibio)
1726 END DO
1727 END DO
1728
1730 DO i=istr,iend
1731 fc(i,k)=(qc(i,k+1)-qc(i,k))*hz_inv2(i,k)
1732 END DO
1733 END DO
1735 DO i=istr,iend
1736 dltr=hz(i,j,k)*fc(i,k)
1737 dltl=hz(i,j,k)*fc(i,k-1)
1738 cff=hz(i,j,k-1)+2.0_r8*hz(i,j,k)+hz(i,j,k+1)
1739 cffr=cff*fc(i,k)
1740 cffl=cff*fc(i,k-1)
1741
1742
1743
1744
1745 IF ((dltr*dltl).le.0.0_r8) THEN
1746 dltr=0.0_r8
1747 dltl=0.0_r8
1748 ELSE IF (abs(dltr).gt.abs(cffl)) THEN
1749 dltr=cffl
1750 ELSE IF (abs(dltl).gt.abs(cffr)) THEN
1751 dltl=cffr
1752 END IF
1753
1754
1755
1756
1757
1758
1759
1760
1761
1762
1763 cff=(dltr-dltl)*hz_inv3(i,k)
1764 dltr=dltr-cff*hz(i,j,k+1)
1765 dltl=dltl+cff*hz(i,j,k-1)
1766 br(i,k)=qc(i,k)+dltr
1767 bl(i,k)=qc(i,k)-dltl
1768 wr(i,k)=(2.0_r8*dltr-dltl)**2
1769 wl(i,k)=(dltr-2.0_r8*dltl)**2
1770 END DO
1771 END DO
1772 cff=1.0e-14_r8
1774 DO i=istr,iend
1775 dltl=max(cff,wl(i,k ))
1776 dltr=max(cff,wr(i,k+1))
1777 br(i,k)=(dltr*br(i,k)+dltl*bl(i,k+1))/(dltr+dltl)
1778 bl(i,k+1)=br(i,k)
1779 END DO
1780 END DO
1781 DO i=istr,iend
1783#if defined LINEAR_CONTINUATION
1784 bl(i,
n(ng))=br(i,
n(ng)-1)
1785 br(i,
n(ng))=2.0_r8*qc(i,
n(ng))-bl(i,
n(ng))
1786#elif defined NEUMANN
1787 bl(i,
n(ng))=br(i,
n(ng)-1)
1788 br(i,
n(ng))=1.5*qc(i,
n(ng))-0.5_r8*bl(i,
n(ng))
1789#else
1790 br(i,
n(ng))=qc(i,
n(ng))
1791 bl(i,
n(ng))=qc(i,
n(ng))
1792 br(i,
n(ng)-1)=qc(i,
n(ng))
1793#endif
1794#if defined LINEAR_CONTINUATION
1795 br(i,1)=bl(i,2)
1796 bl(i,1)=2.0_r8*qc(i,1)-br(i,1)
1797#elif defined NEUMANN
1798 br(i,1)=bl(i,2)
1799 bl(i,1)=1.5_r8*qc(i,1)-0.5_r8*br(i,1)
1800#else
1801 bl(i,2)=qc(i,1)
1802 br(i,1)=qc(i,1)
1803 bl(i,1)=qc(i,1)
1804#endif
1805 END DO
1806
1807
1808
1809
1810
1812 DO i=istr,iend
1813 dltr=br(i,k)-qc(i,k)
1814 dltl=qc(i,k)-bl(i,k)
1815 cffr=2.0_r8*dltr
1816 cffl=2.0_r8*dltl
1817
1818
1819 ad_qc(i,k)=ad_qc(i,k)+ad_bl(i,k)
1820 ad_dltl=ad_dltl-ad_bl(i,k)
1821 ad_bl(i,k)=0.0_r8
1822
1823
1824 ad_qc(i,k)=ad_qc(i,k)+ad_br(i,k)
1825 ad_dltr=ad_dltr+ad_br(i,k)
1826 ad_br(i,k)=0.0_r8
1827 IF ((dltr*dltl).lt.0.0_r8) THEN
1828
1829
1830 ad_dltr=0.0_r8
1831
1832
1833 ad_dltl=0.0_r8
1834 ELSE IF (abs(dltr).gt.abs(cffl)) THEN
1835
1836
1837 ad_cffl=ad_cffl+ad_dltr
1838 ad_dltr=0.0_r8
1839 ELSE IF (abs(dltl).gt.abs(cffr)) THEN
1840
1841
1842 ad_cffr=ad_cffr+ad_dltl
1843 ad_dltl=0.0_r8
1844 END IF
1845
1846
1847 ad_dltl=ad_dltl+2.0_r8*ad_cffl
1848 ad_cffl=0.0_r8
1849
1850
1851 ad_dltr=ad_dltr+2.0_r8*ad_cffr
1852 ad_cffr=0.0_r8
1853
1854
1855 ad_qc(i,k)=ad_qc(i,k)+ad_dltl
1856 ad_bl(i,k)=ad_bl(i,k)-ad_dltl
1857 ad_dltl=0.0_r8
1858
1859
1860 ad_br(i,k)=ad_br(i,k)+ad_dltr
1861 ad_qc(i,k)=ad_qc(i,k)-ad_dltr
1862 ad_dltr=0.0_r8
1863 END DO
1864 END DO
1865 DO i=istr,iend
1866#if defined LINEAR_CONTINUATION
1867
1868
1869 ad_bl(i,2)=ad_bl(i,2)+ad_br(i,1)
1870 ad_br(i,1)=0.0_r8
1871
1872
1873 ad_qc(i,1)=ad_qc(i,1)+2.0_r8*ad_bl(i,1)
1874 ad_br(i,1)=ad_br(i,1)-ad_bl(i,1)
1875 ad_bl(i,1)=0.0_r8
1876#elif defined NEUMANN
1877
1878
1879 ad_bl(i,2)=ad_bl(i,2)+ad_br(i,1)
1880 ad_br(i,1)=0.0_r8
1881
1882
1883 ad_qc(i,1)=ad_qc(i,1)+1.5_r8*ad_bl(i,1)
1884 ad_br(i,1)=ad_br(i,1)-0.5_r8*ad_bl(i,1)
1885 ad_bl(i,1)=0.0_r8
1886#else
1887
1888
1889
1890
1891 ad_qc(i,1)=ad_qc(i,1)+ad_bl(i,1)+ &
1892 & ad_br(i,1)+ &
1893 & ad_bl(i,2)
1894 ad_bl(i,1)=0.0_r8
1895 ad_br(i,1)=0.0_r8
1896 ad_bl(i,2)=0.0_r8
1897#endif
1898#if defined LINEAR_CONTINUATION
1899
1900
1901 ad_br(i,
n(ng)-1)=ad_br(i,
n(ng)-1)+ad_bl(i,
n(ng))
1902 ad_bl(i,
n(ng))=0.0_r8
1903
1904
1905 ad_qc(i,
n(ng))=ad_qc(i,
n(ng))+2.0_r8*ad_br(i,
n(ng))
1906 ad_bl(i,
n(ng))=ad_bl(i,
n(ng))-ad_br(i,
n(ng))
1907 ad_br(i,
n(ng))=0.0_r8
1908#elif defined NEUMANN
1909
1910
1911 ad_br(i,
n(ng)-1)=ad_br(i,
n(ng)-1)+ad_bl(i,
n(ng))
1912 ad_bl(i,
n(ng))=0.0_r8
1913
1914
1915 ad_qc(i,
n(ng))=ad_qc(i,
n(ng))+1.5_r8*ad_br(i,
n(ng))
1916 ad_bl(i,
n(ng))=ad_bl(i,
n(ng))-0.5_r8*ad_br(i,
n(ng))
1917 ad_br(i,
n(ng))=0.0_r8
1918#else
1919
1920
1921
1922
1923 ad_qc(i,
n(ng))=ad_qc(i,
n(ng))+ad_br(i,
n(ng)-1)+ &
1926 ad_br(i,
n(ng)-1)=0.0_r8
1927 ad_bl(i,
n(ng))=0.0_r8
1928 ad_br(i,
n(ng))=0.0_r8
1929#endif
1930 END DO
1931
1932
1933
1934
1935
1936
1937
1938
1940 DO i=istr,iend
1941 qc(i,k)=bio(i,k,ibio)
1942 END DO
1943 END DO
1944
1946 DO i=istr,iend
1947 fc(i,k)=(qc(i,k+1)-qc(i,k))*hz_inv2(i,k)
1948 END DO
1949 END DO
1951 DO i=istr,iend
1952 dltr=hz(i,j,k)*fc(i,k)
1953 dltl=hz(i,j,k)*fc(i,k-1)
1954 cff=hz(i,j,k-1)+2.0_r8*hz(i,j,k)+hz(i,j,k+1)
1955 cffr=cff*fc(i,k)
1956 cffl=cff*fc(i,k-1)
1957
1958
1959
1960
1961 IF ((dltr*dltl).le.0.0_r8) THEN
1962 dltr=0.0_r8
1963 dltl=0.0_r8
1964 ELSE IF (abs(dltr).gt.abs(cffl)) THEN
1965 dltr=cffl
1966 ELSE IF (abs(dltl).gt.abs(cffr)) THEN
1967 dltl=cffr
1968 END IF
1969
1970
1971
1972
1973
1974
1975
1976
1977
1978
1979 cff=(dltr-dltl)*hz_inv3(i,k)
1980 dltr=dltr-cff*hz(i,j,k+1)
1981 dltl=dltl+cff*hz(i,j,k-1)
1982 br(i,k)=qc(i,k)+dltr
1983 bl(i,k)=qc(i,k)-dltl
1984 wr(i,k)=(2.0_r8*dltr-dltl)**2
1985 wl(i,k)=(dltr-2.0_r8*dltl)**2
1986 END DO
1987 END DO
1988
1989 cff=1.0e-14_r8
1991 DO i=istr,iend
1992 dltl=max(cff,wl(i,k ))
1993 dltr=max(cff,wr(i,k+1))
1994 br1(i,k)=br(i,k)
1995 br(i,k)=(dltr*br(i,k)+dltl*bl(i,k+1))/(dltr+dltl)
1996 bl1(i,k+1)=bl(i,k+1)
1997 bl(i,k+1)=br(i,k)
1998
1999
2000 ad_br(i,k)=ad_br(i,k)+ad_bl(i,k+1)
2001 ad_bl(i,k+1)=0.0_r8
2002
2003
2004
2005
2006
2007 adfac=ad_br(i,k)/(dltr+dltl)
2008 adfac1=ad_br(i,k)*br(i,k)/(dltr+dltl)
2009 ad_dltr=ad_dltr+adfac*br1(i,k)
2010 ad_dltl=ad_dltl+adfac*bl1(i,k+1)
2011 ad_bl(i,k+1)=ad_bl(i,k+1)+dltl*adfac
2012 ad_dltr=ad_dltr-adfac1
2013 ad_dltl=ad_dltl-adfac1
2014 ad_br(i,k)=dltr*adfac
2015
2016
2017
2018 ad_wr(i,k+1)=ad_wr(i,k+1)+ &
2019 & (0.5_r8-sign(0.5_r8,cff-wr(i,k+1)))* &
2020 & ad_dltr
2021 ad_dltr=0.0_r8
2022
2023
2024
2025 ad_wl(i,k )=ad_wl(i,k )+ &
2026 & (0.5_r8-sign(0.5_r8,cff-wl(i,k )))* &
2027 & ad_dltl
2028 ad_dltl=0.0_r8
2029 END DO
2030 END DO
2031
2033 DO i=istr,iend
2034
2035
2036
2037 dltr=hz(i,j,k)*fc(i,k)
2038 dltl=hz(i,j,k)*fc(i,k-1)
2039 cff=hz(i,j,k-1)+2.0_r8*hz(i,j,k)+hz(i,j,k+1)
2040 cffr=cff*fc(i,k)
2041 cffl=cff*fc(i,k-1)
2042
2043
2044
2045
2046 IF ((dltr*dltl).le.0.0_r8) THEN
2047 dltr=0.0_r8
2048 dltl=0.0_r8
2049 ELSE IF (abs(dltr).gt.abs(cffl)) THEN
2050 dltr=cffl
2051 ELSE IF (abs(dltl).gt.abs(cffr)) THEN
2052 dltl=cffr
2053 END IF
2054
2055
2056
2057
2058
2059
2060
2061
2062
2063
2064 cff=(dltr-dltl)*hz_inv3(i,k)
2065 dltr=dltr-cff*hz(i,j,k+1)
2066 dltl=dltl+cff*hz(i,j,k-1)
2067
2068
2069
2070 adfac=ad_wl(i,k)*2.0_r8*(dltr-2.0_r8*dltl)
2071 ad_dltr=ad_dltr+adfac
2072 ad_dltl=ad_dltl-2.0_r8*adfac
2073 ad_wl(i,k)=0.0_r8
2074
2075
2076
2077
2078 adfac=ad_wr(i,k)*2.0_r8*(2.0_r8*dltr-dltl)
2079 ad_dltr=ad_dltr+2.0_r8*adfac
2080 ad_dltl=ad_dltl-adfac
2081 ad_wr(i,k)=0.0_r8
2082
2083
2084 ad_qc(i,k)=ad_qc(i,k)+ad_bl(i,k)
2085 ad_dltl=ad_dltl-ad_bl(i,k)
2086 ad_bl(i,k)=0.0_r8
2087
2088
2089 ad_qc(i,k)=ad_qc(i,k)+ad_br(i,k)
2090 ad_dltr=ad_dltr+ad_br(i,k)
2091 ad_br(i,k)=0.0_r8
2092
2093
2094
2095
2096 dltr=hz(i,j,k)*fc(i,k)
2097 dltl=hz(i,j,k)*fc(i,k-1)
2098 cff=hz(i,j,k-1)+2.0_r8*hz(i,j,k)+hz(i,j,k+1)
2099 cffr=cff*fc(i,k)
2100 cffl=cff*fc(i,k-1)
2101
2102
2103
2104
2105 IF ((dltr*dltl).le.0.0_r8) THEN
2106 dltr=0.0_r8
2107 dltl=0.0_r8
2108 ELSE IF (abs(dltr).gt.abs(cffl)) THEN
2109 dltr=cffl
2110 ELSE IF (abs(dltl).gt.abs(cffr)) THEN
2111 dltl=cffr
2112 END IF
2113
2114 cff=(dltr-dltl)*hz_inv3(i,k)
2115
2116
2117 ad_cff=ad_cff+ad_dltl*hz(i,j,k-1)
2118 ad_hz(i,j,k-1)=ad_hz(i,j,k-1)+cff*ad_dltl
2119
2120
2121 ad_cff=ad_cff-ad_dltr*hz(i,j,k+1)
2122 ad_hz(i,j,k+1)=ad_hz(i,j,k+1)-cff*ad_dltr
2123
2124
2125
2126 adfac=ad_cff*hz_inv3(i,k)
2127 ad_dltr=ad_dltr+adfac
2128 ad_dltl=ad_dltl-adfac
2129 ad_hz_inv3(i,k)=ad_hz_inv3(i,k)+(dltr-dltl)*ad_cff
2130 ad_cff=0.0_r8
2131
2132
2133
2134 dltr=hz(i,j,k)*fc(i,k)
2135 dltl=hz(i,j,k)*fc(i,k-1)
2136 cff=hz(i,j,k-1)+2.0_r8*hz(i,j,k)+hz(i,j,k+1)
2137 cffr=cff*fc(i,k)
2138 cffl=cff*fc(i,k-1)
2139
2140 IF ((dltr*dltl).le.0.0_r8) THEN
2141
2142
2143 ad_dltr=0.0_r8
2144
2145
2146 ad_dltl=0.0_r8
2147 ELSE IF (abs(dltr).gt.abs(cffl)) THEN
2148
2149
2150 ad_cffl=ad_cffl+ad_dltr
2151 ad_dltr=0.0_r8
2152 ELSE IF (abs(dltl).gt.abs(cffr)) THEN
2153
2154
2155 ad_cffr=ad_cffr+ad_dltl
2156 ad_dltl=0.0_r8
2157 END IF
2158
2159
2160 ad_cff=ad_cff+ad_cffl*fc(i,k-1)
2161 ad_fc(i,k-1)=ad_fc(i,k-1)+cff*ad_cffl
2162 ad_cffl=0.0_r8
2163
2164
2165 ad_cff=ad_cff+ad_cffr*fc(i,k)
2166 ad_fc(i,k)=ad_fc(i,k)+cff*ad_cffr
2167 ad_cffr=0.0_r8
2168
2169
2170 ad_hz(i,j,k-1)=ad_hz(i,j,k-1)+ad_cff
2171 ad_hz(i,j,k)=ad_hz(i,j,k)+2.0_r8*ad_cff
2172 ad_hz(i,j,k+1)=ad_hz(i,j,k+1)+ad_cff
2173 ad_cff=0.0_r8
2174
2175
2176 ad_hz(i,j,k)=ad_hz(i,j,k)+ad_dltl*fc(i,k-1)
2177 ad_fc(i,k-1)=ad_fc(i,k-1)+ad_dltl*hz(i,j,k)
2178 ad_dltl=0.0_r8
2179
2180
2181 ad_hz(i,j,k)=ad_hz(i,j,k)+ad_dltr*fc(i,k)
2182 ad_fc(i,k)=ad_fc(i,k)+ad_dltr*hz(i,j,k)
2183 ad_dltr=0.0_r8
2184 END DO
2185 END DO
2187 DO i=istr,iend
2188
2189
2190
2191 adfac=ad_fc(i,k)*hz_inv2(i,k)
2192 ad_qc(i,k+1)=ad_qc(i,k+1)+adfac
2193 ad_qc(i,k)=ad_qc(i,k)-adfac
2194 ad_hz_inv2(i,k)=ad_hz_inv2(i,k)+(qc(i,k+1)-qc(i,k))* &
2195 & ad_fc(i,k)
2196 ad_fc(i,k)=0.0_r8
2197 END DO
2198 END DO
2200 DO i=istr,iend
2201
2202
2203 ad_bio(i,k,ibio)=ad_bio(i,k,ibio)+ad_qc(i,k)
2204 ad_qc(i,k)=0.0_r8
2205 END DO
2206 END DO
2207
2208 END DO sink_loop1
2209
2210
2211
2213 DO i=istr,iend
2214
2215
2216
2217
2218
2219
2220
2221
2222
2223
2224
2225
2228
2229
2230 biotrc(ibio,nstp)=t(i,j,k,nstp,ibio)
2231
2232
2233 biotrc(ibio,nnew)=t(i,j,k,nnew,ibio)
2234 END DO
2235
2236
2237
2238 cff2=0.0_r8
2239 DO itime=1,2
2240 cff1=0.0_r8
2242#ifdef IRON_LIMIT
2244#else
2246#endif
2248 cff1=cff1+max(0.0_r8,minval-biotrc(ibio,itime))
2249 IF (biotrc(ibio,itime).gt.biotrc(itrcmax,itime)) THEN
2250 itrcmax=ibio
2251 END IF
2252 biotrc(ibio,itime)=max(minval,biotrc(ibio,itime))
2253 END DO
2254 IF (biotrc(itrcmax,itime).gt.cff1) THEN
2255 biotrc(itrcmax,itime)=biotrc(itrcmax,itime)-cff1
2256 END IF
2257#ifdef IRON_LIMIT
2260 biotrc(ibio,itime)=max(minval,biotrc(ibio,itime))
2261 END DO
2262#endif
2263 END DO
2264
2265
2266
2269 bio_old(i,k,ibio)=biotrc(ibio,nstp)
2270 bio(i,k,ibio)=biotrc(ibio,nstp)
2271 END DO
2272
2273#if defined IRON_LIMIT && defined IRON_RELAX
2274
2275
2276
2277
2278
2279 IF (h(i,j).le.
fehmin(ng))
THEN
2282 END IF
2283#endif
2284 END DO
2285 END DO
2286
2287
2288
2289
2290
2291 DO i=istr,iend
2292#ifdef CONST_PAR
2293
2294
2295
2296 parsur(i)=158.075_r8
2297#else
2299#endif
2300 END DO
2301
2302
2303
2304
2305
2306
2307 DO iteradj=1,iter
2308
2309
2310
2311 DO i=istr,iend
2312 par=parsur(i)
2313 IF (parsur(i).gt.0.0_r8) THEN
2315
2316
2317
2318
2319
2321 & (z_w(i,j,k)-z_w(i,j,k-1))
2322 expatt=exp(-att)
2323 itop=par
2324 par=itop*(1.0_r8-expatt)/att
2325 light(i,k)=par
2326
2327
2328
2329
2330 par=itop*expatt
2331 END DO
2332 ELSE
2334 light(i,k)=0.0_r8
2335 END DO
2336 END IF
2337 END DO
2338
2339
2340
2341
2342
2343
2344#ifdef IRON_LIMIT
2345
2346
2347
2348
2349
2350
2351
2352
2353
2354
2355#endif
2356
2361 DO i=istr,iend
2362#ifdef IRON_LIMIT
2363 fnratio=bio(i,k,
ifphy)/max(minval,bio(i,k,
iphyt))
2364 fcratio=fnratio*fen2fec
2366 flimit=fcratio*fcratio/ &
2368
2370 fnlim=min(1.0_r8,flimit/(bio(i,k,
ino3_)*nlimit))
2371#endif
2372
2373 cff4=1.0_r8/sqrt(cff2+cff3*light(i,k)*light(i,k))
2374 cff=bio(i,k,
iphyt)* &
2375#ifdef IRON_LIMIT
2376 & cff1*cff4*light(i,k)*fnlim*nlimit
2377#else
2378 & cff1*cff4*light(i,k)/ &
2380#endif
2383 & bio(i,k,
ino3_)*cff
2384
2385#ifdef IRON_LIMIT
2386
2387
2388
2389 fac=cff*bio(i,k,
ino3_)*fnratio/ &
2390 & max(minval,bio(i,k,
ifdis))
2393 & bio(i,k,
ifdis)*fac
2394
2395
2396
2397 cff5=dtdays*(fcratioe-fcratio)/
t_fe(ng)
2398 cff6=bio(i,k,
iphyt)*cff5*fec2fen
2399 IF (cff6.ge.0.0_r8) THEN
2400 cff=cff6/max(minval,bio(i,k,
ifdis))
2403 & bio(i,k,
ifdis)*cff
2404 ELSE
2405 cff=-cff6/max(minval,bio(i,k,
ifphy))
2408 & bio(i,k,
ifphy)*cff
2409 END IF
2410#endif
2411 END DO
2412 END DO
2413
2414
2415
2416
2417
2418#ifdef IRON_LIMIT
2419
2420
2421#endif
2422
2423 cff1=dtdays*
zoogr(ng)
2426 DO i=istr,iend
2427 cff=bio(i,k,
izoop)* &
2428 & cff1*(1.0_r8-exp(-
ivlev(ng)*bio(i,k,
iphyt)))/ &
2434 & bio(i,k,
iphyt)*cff2*cff
2439#ifdef IRON_LIMIT
2445#endif
2446 END DO
2447 END DO
2448
2449
2450
2451
2454 cff1=1.0_r8/(1.0_r8+cff2+cff3)
2456 DO i=istr,iend
2459 & bio(i,k,
iphyt)*cff2
2461 & bio(i,k,
iphyt)*cff3
2462#ifdef IRON_LIMIT
2466#endif
2467 END DO
2468 END DO
2469
2470 IF (iteradj.ne.iter) THEN
2471
2472
2473
2474
2477 cff1=1.0_r8/(1.0_r8+cff2+cff3)
2479 DO i=istr,iend
2482 & bio(i,k,
izoop)*cff2
2484 & bio(i,k,
izoop)*cff3
2485 END DO
2486 END DO
2487
2488
2489
2490 cff2=dtdays*
detrr(ng)
2491 cff1=1.0_r8/(1.0_r8+cff2)
2493 DO i=istr,iend
2496 & bio(i,k,
isdet)*cff2
2497 END DO
2498 END DO
2499
2500
2501
2502
2503
2504
2505
2506
2507
2508 DO isink=1,nsink
2509 ibio=idsink(isink)
2510
2511
2512
2513
2514
2515
2517 DO i=istr,iend
2518 qc(i,k)=bio(i,k,ibio)
2519 END DO
2520 END DO
2521
2523 DO i=istr,iend
2524 fc(i,k)=(qc(i,k+1)-qc(i,k))*hz_inv2(i,k)
2525 END DO
2526 END DO
2528 DO i=istr,iend
2529 dltr=hz(i,j,k)*fc(i,k)
2530 dltl=hz(i,j,k)*fc(i,k-1)
2531 cff=hz(i,j,k-1)+2.0_r8*hz(i,j,k)+hz(i,j,k+1)
2532 cffr=cff*fc(i,k)
2533 cffl=cff*fc(i,k-1)
2534
2535
2536
2537
2538 IF ((dltr*dltl).le.0.0_r8) THEN
2539 dltr=0.0_r8
2540 dltl=0.0_r8
2541 ELSE IF (abs(dltr).gt.abs(cffl)) THEN
2542 dltr=cffl
2543 ELSE IF (abs(dltl).gt.abs(cffr)) THEN
2544 dltl=cffr
2545 END IF
2546
2547
2548
2549
2550
2551
2552
2553
2554
2555
2556 cff=(dltr-dltl)*hz_inv3(i,k)
2557 dltr=dltr-cff*hz(i,j,k+1)
2558 dltl=dltl+cff*hz(i,j,k-1)
2559 br(i,k)=qc(i,k)+dltr
2560 bl(i,k)=qc(i,k)-dltl
2561 wr(i,k)=(2.0_r8*dltr-dltl)**2
2562 wl(i,k)=(dltr-2.0_r8*dltl)**2
2563 END DO
2564 END DO
2565 cff=1.0e-14_r8
2567 DO i=istr,iend
2568 dltl=max(cff,wl(i,k ))
2569 dltr=max(cff,wr(i,k+1))
2570 br(i,k)=(dltr*br(i,k)+dltl*bl(i,k+1))/(dltr+dltl)
2571 bl(i,k+1)=br(i,k)
2572 END DO
2573 END DO
2574 DO i=istr,iend
2576#if defined LINEAR_CONTINUATION
2577 bl(i,
n(ng))=br(i,
n(ng)-1)
2578 br(i,
n(ng))=2.0_r8*qc(i,
n(ng))-bl(i,
n(ng))
2579#elif defined NEUMANN
2580 bl(i,
n(ng))=br(i,
n(ng)-1)
2581 br(i,
n(ng))=1.5*qc(i,
n(ng))-0.5_r8*bl(i,
n(ng))
2582#else
2583 br(i,
n(ng))=qc(i,
n(ng))
2584 bl(i,
n(ng))=qc(i,
n(ng))
2585 br(i,
n(ng)-1)=qc(i,
n(ng))
2586#endif
2587#if defined LINEAR_CONTINUATION
2588 br(i,1)=bl(i,2)
2589 bl(i,1)=2.0_r8*qc(i,1)-br(i,1)
2590#elif defined NEUMANN
2591 br(i,1)=bl(i,2)
2592 bl(i,1)=1.5_r8*qc(i,1)-0.5_r8*br(i,1)
2593#else
2594 bl(i,2)=qc(i,1)
2595 br(i,1)=qc(i,1)
2596 bl(i,1)=qc(i,1)
2597#endif
2598 END DO
2599
2600
2601
2602
2603
2605 DO i=istr,iend
2606 dltr=br(i,k)-qc(i,k)
2607 dltl=qc(i,k)-bl(i,k)
2608 cffr=2.0_r8*dltr
2609 cffl=2.0_r8*dltl
2610 IF ((dltr*dltl).lt.0.0_r8) THEN
2611 dltr=0.0_r8
2612 dltl=0.0_r8
2613 ELSE IF (abs(dltr).gt.abs(cffl)) THEN
2614 dltr=cffl
2615 ELSE IF (abs(dltl).gt.abs(cffr)) THEN
2616 dltl=cffr
2617 END IF
2618 br(i,k)=qc(i,k)+dltr
2619 bl(i,k)=qc(i,k)-dltl
2620 END DO
2621 END DO
2622
2623
2624
2625
2626
2627
2628
2629
2630
2631
2632
2633
2634
2635
2636
2637 cff=dtdays*abs(wbio(isink))
2639 DO i=istr,iend
2640 fc(i,k-1)=0.0_r8
2641 wl(i,k)=z_w(i,j,k-1)+cff
2642 wr(i,k)=hz(i,j,k)*qc(i,k)
2643 ksource(i,k)=k
2644 END DO
2645 END DO
2648 DO i=istr,iend
2649 IF (wl(i,k).gt.z_w(i,j,ks)) THEN
2650 ksource(i,k)=ks+1
2651 fc(i,k-1)=fc(i,k-1)+wr(i,ks)
2652 END IF
2653 END DO
2654 END DO
2655 END DO
2656
2657
2658
2660 DO i=istr,iend
2661 ks=ksource(i,k)
2662 cu=min(1.0_r8,(wl(i,k)-z_w(i,j,ks-1))*hz_inv(i,ks))
2663 fc(i,k-1)=fc(i,k-1)+ &
2664 & hz(i,j,ks)*cu* &
2665 & (bl(i,ks)+ &
2666 & cu*(0.5_r8*(br(i,ks)-bl(i,ks))- &
2667 & (1.5_r8-cu)* &
2668 & (br(i,ks)+bl(i,ks)- &
2669 & 2.0_r8*qc(i,ks))))
2670 END DO
2671 END DO
2673 DO i=istr,iend
2674 bio(i,k,ibio)=qc(i,k)+ &
2675 & (fc(i,k)-fc(i,k-1))*hz_inv(i,k)
2676 END DO
2677 END DO
2678 END DO
2679 END IF
2680 END DO
2681
2682
2683
2684
2685
2686
2687 cff2=dtdays*
detrr(ng)
2688 cff1=1.0_r8/(1.0_r8+cff2)
2690 DO i=istr,iend
2691
2692
2693
2695 & cff2*ad_bio(i,k,
ino3_)
2696
2697
2699 END DO
2700 END DO
2701
2702
2703
2704
2707 cff1=1.0_r8/(1.0_r8+cff2+cff3)
2709 DO i=istr,iend
2710
2711
2712
2714 & cff3*ad_bio(i,k,
isdet)
2715
2716
2717
2719 & cff2*ad_bio(i,k,
ino3_)
2720
2721
2723 END DO
2724 END DO
2725
2726
2727
2728
2729
2732 cff1=1.0_r8/(1.0_r8+cff2+cff3)
2734 DO i=istr,iend
2735#ifdef IRON_LIMIT
2736
2737
2738
2740 & (cff2+cff3)*
ferr(ng)*ad_bio(i,k,
ifdis)
2741
2742
2744#endif
2745
2746
2747
2749 & cff3*ad_bio(i,k,
isdet)
2750
2751
2752
2754 & cff2*ad_bio(i,k,
ino3_)
2755
2756
2758 END DO
2759 END DO
2760
2761
2762
2763
2764
2765#ifdef IRON_LIMIT
2766
2767
2768#endif
2769
2770 cff1=dtdays*
zoogr(ng)
2773 DO i=istr,iend
2774 cff=bio1(i,k,
izoop)* &
2775 & cff1*(1.0_r8-exp(-
ivlev(ng)*bio1(i,k,
iphyt)))/ &
2777#ifdef IRON_LIMIT
2778
2779
2780
2781
2785
2786
2787
2788
2789 adfac=ad_bio(i,k,
ifphy)/(1.0_r8+cff)
2790 ad_cff=ad_cff-bio2(i,k,
ifphy)*adfac
2791 ad_bio(i,k,
ifphy)=adfac
2792#endif
2793
2794
2795
2796
2800
2801
2802
2803
2804 ad_cff=ad_cff+ &
2808
2809
2810
2811
2812 ad_cff=ad_cff+ &
2815 & cff2*cff*ad_bio(i,k,
izoop)
2816
2817
2818
2819
2820 adfac=ad_bio(i,k,
iphyt)/(1.0_r8+cff)
2821 ad_cff=ad_cff-bio(i,k,
iphyt)*adfac
2822 ad_bio(i,k,
iphyt)=adfac
2823
2824
2825
2826
2827
2828
2829
2831 adfac=ad_cff/bio1(i,k,
iphyt)
2832 ad_bio(i,k,
iphyt)=ad_bio(i,k,
iphyt)-adfac*cff
2835 & adfac*cff1*fac
2837 & adfac*cff1*(1.0_r8-fac)
2838 ad_cff=0.0_r8
2839 END DO
2840 END DO
2841
2842
2843
2845 DO i=istr,iend
2846
2847
2848
2849
2850
2851
2852
2853
2854
2855
2856
2857
2860
2861
2862 biotrc(ibio,nstp)=t(i,j,k,nstp,ibio)
2863
2864
2865 biotrc(ibio,nnew)=t(i,j,k,nnew,ibio)
2866 END DO
2867
2868
2869
2870 cff2=0.0_r8
2871 DO itime=1,2
2872 cff1=0.0_r8
2874#ifdef IRON_LIMIT
2876#else
2878#endif
2880 cff1=cff1+max(0.0_r8,minval-biotrc(ibio,itime))
2881 IF (biotrc(ibio,itime).gt.biotrc(itrcmax,itime)) THEN
2882 itrcmax=ibio
2883 END IF
2884 biotrc(ibio,itime)=max(minval,biotrc(ibio,itime))
2885 END DO
2886 IF (biotrc(itrcmax,itime).gt.cff1) THEN
2887 biotrc(itrcmax,itime)=biotrc(itrcmax,itime)-cff1
2888 END IF
2889#ifdef IRON_LIMIT
2892 biotrc(ibio,itime)=max(minval,biotrc(ibio,itime))
2893 END DO
2894#endif
2895 END DO
2896
2897
2898
2901 bio_old(i,k,ibio)=biotrc(ibio,nstp)
2902 bio(i,k,ibio)=biotrc(ibio,nstp)
2903 END DO
2904
2905#if defined IRON_LIMIT && defined IRON_RELAX
2906
2907
2908
2909
2910
2911 IF (h(i,j).le.
fehmin(ng))
THEN
2914 END IF
2915#endif
2916 END DO
2917 END DO
2918
2919
2920
2921
2922
2923 DO i=istr,iend
2924#ifdef CONST_PAR
2925
2926
2927
2928 parsur(i)=158.075_r8
2929#else
2931#endif
2932 END DO
2933
2934
2935
2936
2937
2938
2939 DO iteradj=1,iter
2940
2941
2942
2943 DO i=istr,iend
2944 par=parsur(i)
2945 IF (parsur(i).gt.0.0_r8) THEN
2947
2948
2949
2950
2951
2953 & (z_w(i,j,k)-z_w(i,j,k-1))
2954 expatt=exp(-att)
2955 itop=par
2956 par=itop*(1.0_r8-expatt)/att
2957 light(i,k)=par
2958
2959
2960
2961
2962 par=itop*expatt
2963 END DO
2964 ELSE
2966 light(i,k)=0.0_r8
2967 END DO
2968 END IF
2969 END DO
2970
2971
2972
2973
2974
2975
2976#ifdef IRON_LIMIT
2977
2978
2979
2980
2981
2982
2983
2984
2985
2986
2987#endif
2988
2993 DO i=istr,iend
2994#ifdef IRON_LIMIT
2995
2996
2997
2998 fnratio=bio(i,k,
ifphy)/max(minval,bio(i,k,
iphyt))
2999 fcratio=fnratio*fen2fec
3001 flimit=fcratio*fcratio/ &
3003
3005 fnlim=min(1.0_r8,flimit/(bio(i,k,
ino3_)*nlimit))
3006#endif
3007 cff4=1.0_r8/sqrt(cff2+cff3*light(i,k)*light(i,k))
3008 cff=bio(i,k,
iphyt)* &
3009#ifdef IRON_LIMIT
3010 & cff1*cff4*light(i,k)*fnlim*nlimit
3011#else
3012 & cff1*cff4*light(i,k)/ &
3014#endif
3019 & bio(i,k,
ino3_)*cff
3020
3021#ifdef IRON_LIMIT
3022
3023
3024
3025 fac=cff*bio(i,k,
ino3_)*fnratio/ &
3026 & max(minval,bio(i,k,
ifdis))
3032 & bio(i,k,
ifdis)*fac
3034
3035
3036
3037 cff5=dtdays*(fcratioe-fcratio)/
t_fe(ng)
3038 cff6=bio(i,k,
iphyt)*cff5*fec2fen
3039 IF (cff6.ge.0.0_r8) then
3040 cff=cff6/max(minval,bio(i,k,
ifdis))
3043 & bio(i,k,
ifdis)*cff
3044 ELSE
3045 cff=-cff6/max(minval,bio(i,k,
ifphy))
3048 & bio(i,k,
ifphy)*cff
3049 END IF
3050#endif
3051 END DO
3052 END DO
3053
3054 IF (iteradj.ne.iter) THEN
3055
3056
3057
3058
3059
3060#ifdef IRON_LIMIT
3061
3062
3063#endif
3064
3065 cff1=dtdays*
zoogr(ng)
3068 DO i=istr,iend
3069 cff=bio(i,k,
izoop)* &
3070 & cff1*(1.0_r8-exp(-
ivlev(ng)*bio(i,k,
iphyt)))/ &
3074 & bio(i,k,
iphyt)*cff2*cff
3079#ifdef IRON_LIMIT
3083#endif
3084 END DO
3085 END DO
3086
3087
3088
3089
3092 cff1=1.0_r8/(1.0_r8+cff2+cff3)
3094 DO i=istr,iend
3097 & bio(i,k,
iphyt)*cff2
3099 & bio(i,k,
iphyt)*cff3
3100#ifdef IRON_LIMIT
3104#endif
3105 END DO
3106 END DO
3107
3108
3109
3110
3113 cff1=1.0_r8/(1.0_r8+cff2+cff3)
3115 DO i=istr,iend
3118 & bio(i,k,
izoop)*cff2
3120 & bio(i,k,
izoop)*cff3
3121 END DO
3122 END DO
3123
3124
3125
3126 cff2=dtdays*
detrr(ng)
3127 cff1=1.0_r8/(1.0_r8+cff2)
3129 DO i=istr,iend
3132 & bio(i,k,
isdet)*cff2
3133 END DO
3134 END DO
3135
3136
3137
3138
3139
3140
3141
3142
3143
3144 DO isink=1,nsink
3145 ibio=idsink(isink)
3146
3147
3148
3149
3150
3151
3153 DO i=istr,iend
3154 qc(i,k)=bio(i,k,ibio)
3155 END DO
3156 END DO
3157
3159 DO i=istr,iend
3160 fc(i,k)=(qc(i,k+1)-qc(i,k))*hz_inv2(i,k)
3161 END DO
3162 END DO
3164 DO i=istr,iend
3165 dltr=hz(i,j,k)*fc(i,k)
3166 dltl=hz(i,j,k)*fc(i,k-1)
3167 cff=hz(i,j,k-1)+2.0_r8*hz(i,j,k)+hz(i,j,k+1)
3168 cffr=cff*fc(i,k)
3169 cffl=cff*fc(i,k-1)
3170
3171
3172
3173
3174 IF ((dltr*dltl).le.0.0_r8) THEN
3175 dltr=0.0_r8
3176 dltl=0.0_r8
3177 ELSE IF (abs(dltr).gt.abs(cffl)) THEN
3178 dltr=cffl
3179 ELSE IF (abs(dltl).gt.abs(cffr)) THEN
3180 dltl=cffr
3181 END IF
3182
3183
3184
3185
3186
3187
3188
3189
3190
3191
3192 cff=(dltr-dltl)*hz_inv3(i,k)
3193 dltr=dltr-cff*hz(i,j,k+1)
3194 dltl=dltl+cff*hz(i,j,k-1)
3195 br(i,k)=qc(i,k)+dltr
3196 bl(i,k)=qc(i,k)-dltl
3197 wr(i,k)=(2.0_r8*dltr-dltl)**2
3198 wl(i,k)=(dltr-2.0_r8*dltl)**2
3199 END DO
3200 END DO
3201 cff=1.0e-14_r8
3203 DO i=istr,iend
3204 dltl=max(cff,wl(i,k ))
3205 dltr=max(cff,wr(i,k+1))
3206 br(i,k)=(dltr*br(i,k)+dltl*bl(i,k+1))/(dltr+dltl)
3207 bl(i,k+1)=br(i,k)
3208 END DO
3209 END DO
3210 DO i=istr,iend
3212#if defined LINEAR_CONTINUATION
3213 bl(i,
n(ng))=br(i,
n(ng)-1)
3214 br(i,
n(ng))=2.0_r8*qc(i,
n(ng))-bl(i,
n(ng))
3215#elif defined NEUMANN
3216 bl(i,
n(ng))=br(i,
n(ng)-1)
3217 br(i,
n(ng))=1.5*qc(i,
n(ng))-0.5_r8*bl(i,
n(ng))
3218#else
3219 br(i,
n(ng))=qc(i,
n(ng))
3220 bl(i,
n(ng))=qc(i,
n(ng))
3221 br(i,
n(ng)-1)=qc(i,
n(ng))
3222#endif
3223#if defined LINEAR_CONTINUATION
3224 br(i,1)=bl(i,2)
3225 bl(i,1)=2.0_r8*qc(i,1)-br(i,1)
3226#elif defined NEUMANN
3227 br(i,1)=bl(i,2)
3228 bl(i,1)=1.5_r8*qc(i,1)-0.5_r8*br(i,1)
3229#else
3230 bl(i,2)=qc(i,1)
3231 br(i,1)=qc(i,1)
3232 bl(i,1)=qc(i,1)
3233#endif
3234 END DO
3235
3236
3237
3238
3239
3241 DO i=istr,iend
3242 dltr=br(i,k)-qc(i,k)
3243 dltl=qc(i,k)-bl(i,k)
3244 cffr=2.0_r8*dltr
3245 cffl=2.0_r8*dltl
3246 IF ((dltr*dltl).lt.0.0_r8) THEN
3247 dltr=0.0_r8
3248 dltl=0.0_r8
3249 ELSE IF (abs(dltr).gt.abs(cffl)) THEN
3250 dltr=cffl
3251 ELSE IF (abs(dltl).gt.abs(cffr)) THEN
3252 dltl=cffr
3253 END IF
3254 br(i,k)=qc(i,k)+dltr
3255 bl(i,k)=qc(i,k)-dltl
3256 END DO
3257 END DO
3258
3259
3260
3261
3262
3263
3264
3265
3266
3267
3268
3269
3270
3271
3272
3273 cff=dtdays*abs(wbio(isink))
3275 DO i=istr,iend
3276 fc(i,k-1)=0.0_r8
3277 wl(i,k)=z_w(i,j,k-1)+cff
3278 wr(i,k)=hz(i,j,k)*qc(i,k)
3279 ksource(i,k)=k
3280 END DO
3281 END DO
3284 DO i=istr,iend
3285 IF (wl(i,k).gt.z_w(i,j,ks)) THEN
3286 ksource(i,k)=ks+1
3287 fc(i,k-1)=fc(i,k-1)+wr(i,ks)
3288 END IF
3289 END DO
3290 END DO
3291 END DO
3292
3293
3294
3296 DO i=istr,iend
3297 ks=ksource(i,k)
3298 cu=min(1.0_r8,(wl(i,k)-z_w(i,j,ks-1))*hz_inv(i,ks))
3299 fc(i,k-1)=fc(i,k-1)+ &
3300 & hz(i,j,ks)*cu* &
3301 & (bl(i,ks)+ &
3302 & cu*(0.5_r8*(br(i,ks)-bl(i,ks))- &
3303 & (1.5_r8-cu)* &
3304 & (br(i,ks)+bl(i,ks)- &
3305 & 2.0_r8*qc(i,ks))))
3306 END DO
3307 END DO
3309 DO i=istr,iend
3310 bio(i,k,ibio)=qc(i,k)+ &
3311 & (fc(i,k)-fc(i,k-1))*hz_inv(i,k)
3312 END DO
3313 END DO
3314 END DO
3315 END IF
3316 END DO
3317
3318
3319
3320
3321
3322
3323
3324
3325#ifdef IRON_LIMIT
3326
3327
3328
3329
3330
3331
3332
3333
3334
3335
3336#endif
3337
3342 DO i=istr,iend
3343#ifdef IRON_LIMIT
3344
3345
3346
3347 fnratio=bio1(i,k,
ifphy)/max(minval,bio1(i,k,
iphyt))
3348 fcratio=fnratio*fen2fec
3350 flimit=fcratio*fcratio/ &
3352
3354 fnlim=min(1.0_r8,flimit/(bio1(i,k,
ino3_)*nlimit))
3355
3356 cff5=dtdays*(fcratioe-fcratio)/
t_fe(ng)
3357 cff6=bio(i,k,
iphyt)*cff5*fec2fen
3358 IF (cff6.ge.0.0_r8) THEN
3359 fac1=max(minval,bio2(i,k,
ifdis))
3360 cff=cff6/fac1
3361
3362
3363
3364
3366 & cff*ad_bio(i,k,
ifphy)
3367 ad_cff=ad_cff+bio(i,k,
ifdis)*ad_bio(i,k,
ifphy)
3368
3369
3370
3371
3372 adfac=ad_bio(i,k,
ifdis)/(1.0_r8+cff)
3373 ad_cff=ad_cff-bio(i,k,
ifdis)*adfac
3374 ad_bio(i,k,
ifdis)=adfac
3375
3376
3377 adfac=ad_cff/fac1
3378 ad_cff6=ad_cff6+adfac
3379 ad_fac1=ad_fac1-cff*adfac
3380 ad_cff=0.0_r8
3381
3382
3383
3385 & (0.5_r8- &
3386 & sign(0.5_r8, &
3387 & minval-bio2(i,k,
ifdis)))*ad_fac1
3388 ad_fac1=0.0_r8
3389 ELSE
3390 fac1=-max(minval,bio2(i,k,
ifphy))
3391 cff=cff6/fac1
3392
3393
3394
3395
3397 & cff*ad_bio(i,k,
ifdis)
3398 ad_cff=ad_cff+bio(i,k,
ifphy)*ad_bio(i,k,
ifdis)
3399
3400
3401
3402
3403 adfac=ad_bio(i,k,
ifphy)/(1.0_r8+cff)
3404 ad_cff=ad_cff-bio(i,k,
ifphy)*adfac
3405 ad_bio(i,k,
ifphy)=adfac
3406
3407
3408 adfac=ad_cff/fac1
3409 ad_cff6=ad_cff6+adfac
3410 ad_fac1=ad_fac1-cff*adfac
3411 ad_cff=0.0_r8
3412
3413
3414
3416 & (0.5_r8- &
3417 & sign(0.5_r8, &
3418 & minval-bio2(i,k,
ifphy)))*ad_fac1
3419 ad_fac1=0.0_r8
3420 END IF
3421
3422
3423
3424 adfac=ad_cff6*fec2fen
3425 ad_bio(i,k,
iphyt)=ad_bio(i,k,
iphyt)+cff5*adfac
3426 ad_cff5=ad_cff5+bio(i,k,
iphyt)*adfac
3427 ad_cff6=0.0_r8
3428
3429
3430 adfac=dtdays*ad_cff5/
t_fe(ng)
3431 ad_fcratioe=ad_fcratioe+adfac
3432 ad_fcratio=ad_fcratio-adfac
3433 ad_cff5=0.0_r8
3434#endif
3435 cff4=1.0_r8/sqrt(cff2+cff3*light(i,k)*light(i,k))
3436 cff=bio1(i,k,
iphyt)* &
3437#ifdef IRON_LIMIT
3438 & cff1*cff4*light(i,k)*fnlim*nlimit
3439#else
3440 & cff1*cff4*light(i,k)/ &
3442#endif
3443#ifdef IRON_LIMIT
3444
3445
3446
3447 fac1=max(minval,bio1(i,k,
ifdis))
3448 fac2=1.0_r8/fac1
3449 fac=cff*bio(i,k,
ino3_)*fnratio*fac2
3450
3451
3452
3453
3454 ad_fac=ad_fac+bio2(i,k,
ifdis)*ad_bio(i,k,
ifphy)
3456 & fac*ad_bio(i,k,
ifphy)
3457
3458
3459
3460
3461 adfac=ad_bio(i,k,
ifdis)/(1.0_r8+fac)
3462 ad_fac=ad_fac-bio2(i,k,
ifdis)*adfac
3463 ad_bio(i,k,
ifdis)=adfac
3464
3465
3466
3467
3468
3469 adfac1=fnratio*fac2*ad_fac
3470 adfac2=cff*bio(i,k,
ino3_)*ad_fac
3471 ad_cff=ad_cff+bio(i,k,
ino3_)*adfac1
3472 ad_bio(i,k,
ino3_)=ad_bio(i,k,
ino3_)+cff*adfac1
3473 ad_fnratio=ad_fnratio+fac2*adfac2
3474 ad_fac2=ad_fac2+fnratio*adfac2
3475 ad_fac=0.0_r8
3476
3477
3478 ad_fac1=ad_fac1-fac2*fac2*ad_fac2
3479 ad_fac2=0.0_r8
3480
3481
3482
3484 & (0.5_r8- &
3485 & sign(0.5_r8,minval-bio1(i,k,
ifdis)))* &
3486 & ad_fac1
3487 ad_fac1=0.0_r8
3488#endif
3489
3490
3491
3492
3493
3494
3495
3496 ad_cff=ad_cff+bio(i,k,
ino3_)*ad_bio(i,k,
iphyt)
3498 & cff*ad_bio(i,k,
iphyt)
3499
3500
3501
3502
3503 adfac=ad_bio(i,k,
ino3_)/(1.0_r8+cff)
3504 ad_cff=ad_cff-bio(i,k,
ino3_)*adfac
3505 ad_bio(i,k,
ino3_)=adfac
3506
3507#ifdef IRON_LIMIT
3508
3509
3510
3511
3512
3513
3514
3515
3516
3517 adfac1=cff1*cff4*ad_cff
3518 adfac2=adfac1*bio1(i,k,
iphyt)
3520 & light(i,k)*fnlim*nlimit*adfac1
3521 ad_light(i,k)=ad_light(i,k)+ &
3522 & fnlim*nlimit*adfac2
3523 ad_fnlim=ad_fnlim+ &
3524 & light(i,k)*nlimit*adfac2
3525 ad_nlimit=ad_nlimit+ &
3526 & light(i,k)*fnlim*adfac2
3527 ad_cff4=ad_cff4+ &
3528 & bio1(i,k,
iphyt)*cff1*light(i,k)*fnlim*nlimit* &
3529 & ad_cff
3530 ad_cff=0.0_r8
3531#else
3532
3533
3534
3535
3536
3537
3539 adfac1=adfac*bio1(i,k,
iphyt)*cff1
3541 & cff1*cff4*light(i,k)*adfac
3542 ad_cff4=adfac1*light(i,k)
3543 ad_light(i,k)=ad_light(i,k)+ &
3544 & adfac1*cff4
3546 & adfac*cff
3547 ad_cff=0.0_r8
3548#endif
3549
3550
3551 ad_light(i,k)=ad_light(i,k)- &
3552 & cff3*light(i,k)* &
3553 & cff4*cff4*cff4*ad_cff4
3554 ad_cff4=0.0_r8
3555
3556#ifdef IRON_LIMIT
3557
3558
3559
3561 fac1=flimit/(bio1(i,k,
ino3_)*nlimit)
3562
3563
3564 ad_fac1=(0.5_r8+sign(0.5_r8,1.0_r8-fac1))*ad_fnlim
3565 ad_fnlim=0.0_r8
3566
3567
3568
3569
3570
3571 adfac1=ad_fac1/(bio1(i,k,
ino3_)*nlimit)
3572 adfac2=adfac1*fac1
3573 ad_flimit=ad_flimit+adfac1
3574 ad_bio(i,k,
ino3_)=ad_bio(i,k,
ino3_)-nlimit*adfac2
3575 ad_nlimit=ad_nlimit-bio1(i,k,
ino3_)*adfac2
3576 ad_fac1=0.0_r8
3577
3578
3580 & ad_nlimit*nlimit*nlimit
3581 ad_nlimit=0.0_r8
3582
3583
3584
3585
3586 adfac=2.0_r8*ad_flimit/ &
3588 ad_fcratio=ad_fcratio+ &
3589 & fcratio*adfac- &
3590 & fcratio*flimit*adfac
3591 ad_flimit=0.0_r8
3592
3593
3594
3595
3599 & ad_fcratioe
3600 ad_fcratioe=0.0_r8
3601
3602
3603 ad_fnratio=ad_fnratio+fen2fec*ad_fcratio
3604 ad_fcratio=0.0_r8
3605
3606 fac1=max(minval,bio1(i,k,
iphyt))
3607
3608
3609 adfac=ad_fnratio/fac1
3611 ad_fac1=ad_fac1-ad_fnratio*fnratio/fac1
3612 ad_fnratio=0.0_r8
3613
3614
3615
3617 & (0.5_r8- &
3618 & sign(0.5_r8,minval-bio1(i,k,
iphyt)))* &
3619 & ad_fac1
3620 ad_fac1=0.0_r8
3621#endif
3622 END DO
3623 END DO
3624
3625
3626
3627 DO i=istr,iend
3628 par=parsur(i)
3629 IF (parsur(i).gt.0.0_r8) THEN
3631
3632
3633
3634 par=parsur(i)
3636
3637
3638
3639
3640
3642 & (z_w(i,j,kk)-z_w(i,j,kk-1))
3643 expatt=exp(-att)
3644 itop=par
3645 par=itop*(1.0_r8-expatt)/att
3646 par1=par
3647
3648
3649
3650
3651 par=itop*expatt
3652 END DO
3653
3654
3655
3656
3657
3658
3659 ad_expatt=ad_expatt+itop*ad_par
3660 ad_itop=ad_itop+expatt*ad_par
3661 ad_par=0.0_r8
3662
3663
3664
3665
3666
3667
3668
3669 ad_par=ad_par+ad_light(i,k)
3670 ad_light(i,k)=0.0_r8
3671
3672
3673
3674 adfac=ad_par/att
3675 ad_att=ad_att-par1*adfac
3676 ad_expatt=ad_expatt-itop*adfac
3677 ad_itop=ad_itop+(1.0_r8-expatt)*adfac
3678 ad_par=0.0_r8
3679
3680
3681 ad_par=ad_par+ad_itop
3682 ad_itop=0.0_r8
3683
3684
3685 ad_att=ad_att-expatt*ad_expatt
3686 ad_expatt=0.0_r8
3687
3688
3689
3690
3691
3694 &
attphy(ng)*(z_w(i,j,k)-z_w(i,j,k-1))* &
3695 & ad_att
3696 ad_z_w(i,j,k-1)=ad_z_w(i,j,k-1)-adfac
3697 ad_z_w(i,j,k )=ad_z_w(i,j,k )+adfac
3698 ad_att=0.0_r8
3699 END DO
3700 ELSE
3702
3703
3704 ad_light(i,k)=0.0_r8
3705 END DO
3706 END IF
3707
3708
3709 ad_parsur(i)=ad_parsur(i)+ad_par
3710 ad_par=0.0_r8
3711 END DO
3712
3713 END DO iter_loop1
3714
3715
3716
3717
3718
3719
3720 DO i=istr,iend
3721#ifdef CONST_PAR
3722
3723
3724
3725
3726
3727
3728#else
3729
3730
3731
3732 adfac=
rho0*
cp*ad_parsur(i)
3733 ad_srflx(i,j)=ad_srflx(i,j)+
parfrac(ng)*adfac
3735
3736#endif
3737 END DO
3738
3739
3740
3741
3742
3743
3745 DO i=istr,iend
3746
3747#if defined IRON_LIMIT && defined IRON_RELAX
3748
3749
3750
3751
3752
3753 IF (h(i,j).le.
fehmin(ng))
THEN
3754
3755
3756
3758 & fenudgcoef*ad_bio(i,k,
ifdis)
3759 END IF
3760#endif
3761
3762
3763
3766
3767
3768 ad_biotrc(ibio,nstp)=ad_biotrc(ibio,nstp)+ &
3769 & ad_bio(i,k,ibio)
3770 ad_bio(i,k,ibio)=0.0_r8
3771
3772
3773 ad_biotrc(ibio,nstp)=ad_biotrc(ibio,nstp)+ &
3774 & ad_bio_old(i,k,ibio)
3775 ad_bio_old(i,k,ibio)=0.0_r8
3776 END DO
3777
3778
3779
3780 DO itime=1,2
3783
3784
3785
3786
3787
3788 biotrc(ibio,itime)=t(i,j,k,itime,ibio)
3789 biotrc1(ibio,itime)=biotrc(ibio,itime)
3790 END DO
3791 END DO
3792
3793 cff2=0.0_r8
3794 DO itime=1,2
3795 cff1=0.0_r8
3797#ifdef IRON_LIMIT
3799#else
3801#endif
3803 cff1=cff1+max(0.0_r8,minval-biotrc(ibio,itime))
3804 IF (biotrc(ibio,itime).gt.biotrc(itrcmax,itime)) THEN
3805 itrcmax=ibio
3806 END IF
3807 biotrc(ibio,itime)=max(minval,biotrc(ibio,itime))
3808 END DO
3809 IF (biotrc(itrcmax,itime).gt.cff1) THEN
3810 biotrc(itrcmax,itime)=biotrc(itrcmax,itime)-cff1
3811 END IF
3812
3813#ifdef IRON_LIMIT
3816 biotrc(ibio,itime)=max(minval,biotrc1(ibio,itime))
3817
3818
3819
3820
3821
3822
3823 ad_biotrc(ibio,itime)=(0.5_r8- &
3824 & sign(0.5_r8, &
3825 & minval- &
3826 & biotrc1(ibio,itime)))* &
3827 & ad_biotrc(ibio,itime)
3828 END DO
3829#endif
3830
3831
3832
3833 cff1=0.0_r8
3835#ifdef IRON_LIMIT
3837#else
3839#endif
3841
3842
3843
3844
3845
3846 biotrc(ibio,itime)=t(i,j,k,itime,ibio)
3847 cff1=cff1+max(0.0_r8,minval-biotrc(ibio,itime))
3848 IF (biotrc(ibio,itime).gt.biotrc(itrcmax,itime)) THEN
3849 itrcmax=ibio
3850 END IF
3851 biotrc1(ibio,itime)=biotrc(ibio,itime)
3852 biotrc(ibio,itime)=max(minval,biotrc(ibio,itime))
3853 END DO
3854 IF (biotrc(itrcmax,itime).gt.cff1) THEN
3855
3856
3857
3858 ad_cff1=-ad_biotrc(itrcmax,itime)
3859 END IF
3860
3861 cff1=0.0_r8
3862#ifdef IRON_LIMIT
3864#else
3866#endif
3868
3869
3870
3871
3872
3873 biotrc(ibio,itime)=t(i,j,k,itime,ibio)
3874 cff1=cff1+max(0.0_r8,minval-biotrc(ibio,itime))
3875
3876
3877
3878
3879
3880
3881 ad_biotrc(ibio,itime)=(0.5_r8- &
3882 & sign(0.5_r8, &
3883 & minval- &
3884 & biotrc1(ibio,itime)))* &
3885 & ad_biotrc(ibio,itime)
3886
3887
3888
3889
3890
3891 ad_biotrc(ibio,itime)=ad_biotrc(ibio,itime)- &
3892 & (0.5_r8-sign(0.5_r8, &
3893 & biotrc(ibio,itime)- &
3894 & minval))*ad_cff1
3895 END DO
3896 ad_cff1=0.0_r8
3897 END DO
3898
3899
3900
3901
3902
3903
3904
3905
3906
3907
3908
3909
3912
3913
3914
3915
3916
3917 ad_hz_inv(i,k)=ad_hz_inv(i,k)+ &
3918 & t(i,j,k,nnew,ibio)*hz(i,j,k)* &
3919 & ad_biotrc(ibio,nnew)
3920 ad_t(i,j,k,nnew,ibio)=ad_t(i,j,k,nnew,ibio)+ &
3921 & hz_inv(i,k)*ad_biotrc(ibio,nnew)
3922 ad_biotrc(ibio,nnew)=0.0_r8
3923
3924
3925 ad_t(i,j,k,nstp,ibio)=ad_t(i,j,k,nstp,ibio)+ &
3926 & ad_biotrc(ibio,nstp)
3927 ad_biotrc(ibio,nstp)=0.0_r8
3928 END DO
3929 END DO
3930 END DO
3931
3932
3933
3935 DO i=istr,iend
3936
3937
3938
3939
3940 adfac=hz_inv3(i,k)*hz_inv3(i,k)*ad_hz_inv3(i,k)
3941 ad_hz(i,j,k-1)=ad_hz(i,j,k-1)-adfac
3942 ad_hz(i,j,k )=ad_hz(i,j,k )-adfac
3943 ad_hz(i,j,k+1)=ad_hz(i,j,k+1)-adfac
3944 ad_hz_inv3(i,k)=0.0_r8
3945 END DO
3946 END DO
3948 DO i=istr,iend
3949
3950
3951
3952 adfac=hz_inv2(i,k)*hz_inv2(i,k)*ad_hz_inv2(i,k)
3953 ad_hz(i,j,k )=ad_hz(i,j,k )-adfac
3954 ad_hz(i,j,k+1)=ad_hz(i,j,k+1)-adfac
3955 ad_hz_inv2(i,k)=0.0_r8
3956 END DO
3957 END DO
3959 DO i=istr,iend
3960
3961
3962 ad_hz(i,j,k)=ad_hz(i,j,k)- &
3963 & hz_inv(i,k)*hz_inv(i,k)*ad_hz_inv(i,k)
3964 ad_hz_inv(i,k)=0.0_r8
3965 END DO
3966 END DO
3967
3968 END DO j_loop
3969
3970
3971
3972
3973
3974
3976 ad_wbio(2)=0.0_r8
3977
3978
3980 ad_wbio(1)=0.0_r8
3981
3982 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 fenudgtime
real(r8), dimension(:), allocatable ivlev
real(r8), dimension(:), allocatable attphy
real(r8), dimension(:), allocatable b_fe
real(r8), dimension(:), allocatable ad_parfrac
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
real(r8), dimension(:), allocatable ad_wphy