-
Notifications
You must be signed in to change notification settings - Fork 26
/
Copy pathmo_ocprod.F90
1658 lines (1480 loc) · 74.6 KB
/
mo_ocprod.F90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
! Copyright (C) 2001 Ernst Maier-Reimer, S. Legutke
! Copyright (C) 2020 K. Assmann, J. Tjiputra, J. Schwinger, I. Kriest,
! A. Moree, C. Heinze
!
! This file is part of BLOM/iHAMOCC.
!
! BLOM is free software: you can redistribute it and/or modify it under the
! terms of the GNU Lesser General Public License as published by the Free
! Software Foundation, either version 3 of the License, or (at your option)
! any later version.
!
! BLOM is distributed in the hope that it will be useful, but WITHOUT ANY
! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
! FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for
! more details.
!
! You should have received a copy of the GNU Lesser General Public License
! along with BLOM. If not, see https://www.gnu.org/licenses/.
module mo_ocprod
implicit none
private
public :: ocprod
contains
subroutine ocprod(kpie,kpje,kpke,kbnd,pdlxp,pdlyp,pddpo,omask,ptho,pi_ph,psao,ppao,prho)
!***********************************************************************************************
! Biological production, remineralization and particle sinking.
!
! Ernst Maier-Reimer, *MPI-Met, HH* 10.04.01
!
! Modified
! S.Legutke, *MPI-MaD, HH* 2010-04-01
! J.Schwinger, *GFI, UiB* 2013-04-22
! - Corrected bug in light penetration formulation
! - Cautious code clean-up
! J.Tjiputra, *UNI-RESEARCH* 2015-11-25
! - Implemented natural DIC/ALK/CALC
! I.Kriest, *GEOMAR* 2016-08-11
! - Modified stoichiometry for denitrification (affects NO3, N2, Alk)
! J.Schwinger, *UNI-RESEARCH* 2017-08-30
! - Removed split of the layer that only partly falls into the euphotic zone. Loops are
! now calculated over
! (1) layers that are completely or partly in the euphotoc zone
! (2) layers that do not lie within the euphotic zone.
! - Moved the accumulation of global fields for output to routine hamocc4bgc.
! The accumulation of local fields has been moved to the end of this routine.
! A.Moree, *GFI, Bergen* 2018-04-12
! - new version of carbon isotope code
! J.Schwinger, *Uni Research, Bergen* 2018-04-12
! - moved accumulation of all output fields to seperate subroutine,
! related code-restructuring
! - added sediment bypass preprocessor option and related code
! J.Schwinger, *NORCE Climate, Bergen* 2020-05-29
! - Cleaned up parameter list
! - Dust deposition field now passed as an argument
!***********************************************************************************************
use mod_xc, only: mnproc
use mo_carbch, only: ocetra,satoxy,hi,co2star
use mo_sedmnt, only: prcaca,produs,prorca,silpro,pror13,pror14,prca13,prca14
use mo_param_bgc, only: drempoc,drempoc_anaerob,bkox_drempoc,dremn2o,dremopal,dremsul, &
dyphy,ecan,epsher,fesoly, &
gammap,gammaz,grami,grazra,pi_alpha,phytomi, &
rcalc,rcar,rdn2o1,rdn2o2,rdnit0,rdnit1,rdnit2, &
relaxfe,remido,riron,rnit,rnoi,ro2ut,ropal, &
spemor,wcal_const,wdust_const,wopal_const,wpoc_const, &
zinges,alar1,alar2,alar3, &
alow1,alow2,alow3,calmax,cellmass, &
cellsink,dustd1,dustd2,dustd3,dustsink,fractdim, &
fse,fsh,nmldmin,plower,pupper,sinkexp,stick,tmfac, &
tsfac,vsmall,zdis,wmin,wmax,wlin,rbro, &
dmsp1,dmsp2,dmsp3,dmsp4,dmsp5,dmsp6,dms_gamma, &
fbro1,fbro2,atten_f,atten_c,atten_uv,atten_w,bkopal,bkphy,bkzoo, &
POM_remin_q10,POM_remin_Tref,opal_remin_q10,opal_remin_Tref, &
bkphyanh4,bkphyano3,bkphosph,bkiron,ro2utammo,max_limiter, &
O2thresh_aerob,O2thresh_hypoxic,NO3thresh_sulf
use mo_biomod, only: bsiflx0100,bsiflx0500,bsiflx1000,bsiflx2000,bsiflx4000,bsiflx_bot, &
calflx0100,calflx0500,calflx1000,calflx2000,calflx4000,calflx_bot, &
carflx0100,carflx0500,carflx1000,carflx2000,carflx4000,carflx_bot, &
dustflx0100,dustflx0500,dustflx1000,dustflx2000,dustflx4000, &
dustflx_bot, &
expoor,exposi,expoca,intdnit,intdms_bac,intdmsprod,intdms_uv, &
intphosy,int_chbr3_prod,int_chbr3_uv, &
phosy3d,abs_oce,strahl,asize3d,wmass,wnumb,eps3d,phosy_NH4, &
phosy_NO3, remin_aerob,remin_sulf
use mo_param1_bgc, only: ialkali,ian2o,iano3,icalc,idet,idms,idoc,ifdust, &
igasnit,iiron,iopal,ioxygen,iphosph,iphy,isco212, &
isilica,izoo,iadust,inos,ibromo, &
icalc13,icalc14,idet13,idet14,idoc13,idoc14, &
iphy13,iphy14,isco213,isco214,izoo13,izoo14,safediv, &
inatalkali,inatcalc,inatsco212,ianh4
use mo_control_bgc, only: dtb,io_stdo_bgc,with_dmsph, &
use_BROMO,use_AGG,use_PBGC_OCNP_TIMESTEP,use_FB_BGC_OCE, &
use_AGG,use_cisonew,use_natDIC, use_WLIN,use_sedbypass,use_M4AGO, &
use_extNcycle,lkwrbioz_off,lTO2depremin
use mo_vgrid, only: dp_min,dp_min_sink,k0100,k0500,k1000,k2000,k4000,kwrbioz,ptiestu
use mo_vgrid, only: kmle
use mo_clim_swa, only: swa_clim
use mo_inventory_bgc, only: inventory_bgc
use mo_ihamocc4m4ago, only: ihamocc_mean_aggregate_sinking_speed,ws_agg
use mo_extNwatercol, only: nitrification,denit_NO3_to_NO2,anammox,denit_dnra,extN_inv_check
! Arguments
integer, intent(in) :: kpie ! 1st dimension of model grid.
integer, intent(in) :: kpje ! 2nd dimension of model grid.
integer, intent(in) :: kpke ! 3rd (vertical) dimension of model grid.
integer, intent(in) :: kbnd ! nb of halo grid points
real, intent(in) :: pdlxp(kpie,kpje) ! size of grid cell (1st dimension) [m].
real, intent(in) :: pdlyp(kpie,kpje) ! size of grid cell (2nd dimension) [m].
real, intent(in) :: pddpo(kpie,kpje,kpke) ! size of grid cell (3rd dimension) [m].
real, intent(in) :: omask(kpie,kpje) ! land/ocean mask (1=ocean)
real, intent(in) :: ptho(1-kbnd:kpie+kbnd,1-kbnd:kpje+kbnd,kpke) ! potential temperature [deg C].
real, intent(in) :: pi_ph(kpie,kpje)
real, intent(in) :: psao(1-kbnd:kpie+kbnd,1-kbnd:kpje+kbnd,kpke) ! salinity [psu].
real, intent(in) :: ppao(1-kbnd:kpie+kbnd,1-kbnd:kpje+kbnd) ! sea level pressure [Pascal].
real, intent(in) :: prho(kpie,kpje,kpke) ! density [g/cm^3].
! Local variables
integer, parameter :: nsinkmax = 12
integer :: i,j,k,l
integer :: is,kdonor
real :: abs_bgc(kpie,kpje,kpke)
real :: tco(nsinkmax),tcn(nsinkmax),q(nsinkmax)
real :: atten,avphy,avanut,avanfe,pho,xa,xn,ya,yn,phosy
real :: avgra,grazing,avsil,avdic,graton
real :: gratpoc,grawa,bacfra,phymor,zoomor,excdoc,exud
real :: export, delsil, delcar, sterph, sterzo, remin
real :: docrem, opalrem, remin2o, aou,refra,pocrem,phyrem
real :: zoothresh,phythresh
real :: temp,temfa,phofa ! temperature and irradiation factor for photosynthesis
real :: absorption,absorption_uv
real :: dmsprod,dms_bac,dms_uv,dms_ph
real :: dtr,dz
real :: wpocd,wcald,wopald,wdustd,dagg
real :: wcal,wdust,wopal,wpoc
real :: o2lim ! O2 limitation of ammonification (POC remin)
! sedbypass
real :: florca,flcaca,flsil
! cisonew
real :: phygrowth
real :: phosy13,phosy14
real :: growth_co2
real :: bifr13,bifr14,bifr13_perm
real :: grazing13,grazing14
real :: graton13,graton14
real :: gratpoc13,gratpoc14
real :: bacfra13,bacfra14
real :: phymor13,phymor14
real :: grawa13,grawa14
real :: zoomor13,zoomor14
real :: excdoc13,excdoc14
real :: exud13,exud14
real :: export13,export14
real :: delcar13,delcar14
real :: dtr13,dtr14
real :: sterph13,sterph14
real :: sterzo13,sterzo14
real :: pocrem13,pocrem14
real :: docrem13,docrem14
real :: phyrem13,phyrem14
real :: rem13,rem14
real :: rco213,rco214,rdoc13,rdoc14,rdet13,rdet14
real :: rphy13,rphy14,rzoo13,rzoo14
! sedbypass
real :: flor13,flor14,flca13,flca14
! AGG
real :: aggregate(kpie,kpje,kpke)
real :: dustagg(kpie,kpje,kpke)
real :: avmass, avnos, anosloss
real :: zmornos, eps, e1,e2,e3,e4,es1,es3
real :: TopM,TopF, snow,fshear,sagg1,sagg2,sagg4
real :: sett_agg,shear_agg,effsti,dfirst,dshagg,dsett
real :: wnos,wnosd
! BROMO
real :: bro_beta,bro_uv
real :: abs_uv(kpie,kpje,kpke)
! extNcycle
character(len=:), allocatable :: inv_message
real :: ano3up_inh,nutlim,anh4lim,nlim,grlim,nh4uptfrac
! set variables for diagnostic output to zero
expoor (:,:) = 0.
expoca (:,:) = 0.
exposi (:,:) = 0.
carflx0100 (:,:) = 0.
carflx0500 (:,:) = 0.
carflx1000 (:,:) = 0.
carflx2000 (:,:) = 0.
carflx4000 (:,:) = 0.
carflx_bot (:,:) = 0.
bsiflx0100 (:,:) = 0.
bsiflx0500 (:,:) = 0.
bsiflx1000 (:,:) = 0.
bsiflx2000 (:,:) = 0.
bsiflx4000 (:,:) = 0.
bsiflx_bot (:,:) = 0.
calflx0100 (:,:) = 0.
calflx0500 (:,:) = 0.
calflx1000 (:,:) = 0.
calflx2000 (:,:) = 0.
calflx4000 (:,:) = 0.
calflx_bot (:,:) = 0.
dustflx0100(:,:) = 0.
dustflx0500(:,:) = 0.
dustflx1000(:,:) = 0.
dustflx2000(:,:) = 0.
dustflx4000(:,:) = 0.
dustflx_bot(:,:) = 0.
intdnit (:,:) = 0.
intphosy (:,:) = 0.
intdmsprod (:,:) = 0.
intdms_bac (:,:) = 0.
intdms_uv (:,:) = 0.
phosy3d (:,:,:) = 0.
if (use_BROMO) then
int_chbr3_uv (:,:) = 0.
int_chbr3_prod(:,:) = 0.
end if
if (use_AGG) then
eps3d(:,:,:) = 0.
asize3d(:,:,:) = 0.
endif
if (use_extNcycle) then
phosy_NH4(:,:,:) = 0.
phosy_NO3(:,:,:) = 0.
remin_aerob(:,:,:) = 0.
remin_sulf(:,:,:) = 0.
endif
if (use_PBGC_OCNP_TIMESTEP) then
if (mnproc == 1) then
write(io_stdo_bgc,*)' '
write(io_stdo_bgc,*)'beginning of OCRPOD '
endif
call inventory_bgc(kpie,kpje,kpke,pdlxp,pdlyp,pddpo,omask,0)
endif
! Calculate swr absorption by water and phytoplankton
abs_bgc(:,:,:) = 0.
if (use_BROMO) then
abs_uv(:,:,:) = 0.
endif
if (use_FB_BGC_OCE) then
abs_oce(:,:,:) = 0.
abs_oce(:,:,1) = 1.
endif
!$OMP PARALLEL DO PRIVATE(i,k,absorption,absorption_uv,atten,dz)
do j = 1,kpje
do i = 1,kpie
if(omask(i,j) > 0.5) then
absorption = 1.
absorption_uv = 1.
vloop: do k = 1,merge(kpke,kwrbioz(i,j),lkwrbioz_off)
if(pddpo(i,j,k) > 0.0) then
dz = pddpo(i,j,k)
! Average light intensity in layer k
atten = atten_w + atten_c * max(0.,ocetra(i,j,k,iphy))
abs_bgc(i,j,k) = ((absorption/atten)* (1.-exp(-atten*dz)))/dz
if (use_BROMO) then
abs_uv(i,j,k) = ((absorption_uv/atten_uv)*(1.-exp(-atten_uv*dz)))/dz
endif
if (use_FB_BGC_OCE) then
abs_oce(i,j,k) = abs_oce(i,j,k) * absorption
if (k == 2) then
abs_oce(i,j,2) = atten_f * absorption
endif
endif
! Radiation intensity I_0 at the top of next layer
absorption = absorption * exp(-atten*dz)
absorption_uv = absorption_uv * exp(-atten_uv*dz)
endif
enddo vloop
endif ! omask > 0.5
enddo
enddo
!$OMP END PARALLEL DO
if (use_M4AGO) then
! even though we loose detritus, etc. we call the calculation for settling velocity by M4AGO here
! to enable further future development... - assuming that the operator splitting decently functions
call ihamocc_mean_aggregate_sinking_speed(kpie, kpje, kpke, kbnd, pddpo, omask, ptho, psao, ppao, prho)
endif
!$OMP PARALLEL DO PRIVATE(avphy,avgra,avsil,avanut,avanfe,pho,xa,xn &
!$OMP ,phosy,ya,yn,grazing,graton,gratpoc,grawa,bacfra,phymor &
!$OMP ,zoomor,excdoc,exud,export,delsil,delcar,dmsprod &
!$OMP ,dms_bac,dms_uv,dtr,phofa,temfa,zoothresh,dms_ph,dz,opalrem &
!$OMP ,avmass,avnos,zmornos &
!$OMP ,rco213,rco214,rphy13,rphy14,rzoo13,rzoo14,grazing13,grazing14 &
!$OMP ,graton13,graton14,gratpoc13,gratpoc14,grawa13,grawa14 &
!$OMP ,phosy13,phosy14,bacfra13,bacfra14,phymor13,phymor14,zoomor13 &
!$OMP ,zoomor14,excdoc13,excdoc14,exud13,exud14,export13,export14 &
!$OMP ,delcar13,delcar14,dtr13,dtr14,bifr13,bifr14,bifr13_perm &
!$OMP ,growth_co2,phygrowth &
!$OMP ,bro_beta,bro_uv &
!$OMP ,ano3up_inh,nutlim,anh4lim,nlim,grlim,nh4uptfrac &
!$OMP ,i,k)
loop1: do j = 1,kpje
do i = 1,kpie
do k = 1,merge(kpke,kwrbioz(i,j),lkwrbioz_off)
if(pddpo(i,j,k) > dp_min .and. omask(i,j) > 0.5) then
if (use_AGG) then
avmass = ocetra(i,j,k,iphy) + ocetra(i,j,k,idet)
endif
temp = min(40.,max(-3.,ptho(i,j,k)))
phofa = pi_alpha * strahl(i,j) * abs_bgc(i,j,k)
temfa = 0.6 * 1.066**temp
!taylor: temfa= 0.6*(1. + 0.0639*ptho(i,j,k) * &
! & (1. + 0.0639*ptho(i,j,k)/2. * (1. + 0.0639*ptho(i,j,k)/3.)))
pho = dtb * phofa * temfa / sqrt(phofa**2 + temfa**2)
avphy = max(phytomi,ocetra(i,j,k,iphy)) ! 'available' phytoplankton
avgra = max(grami,ocetra(i,j,k,izoo)) ! 'available' zooplankton
avsil = max(0.,ocetra(i,j,k,isilica))
avdic = max(0.,ocetra(i,j,k,isco212))
if (use_extNcycle)then
ano3up_inh = bkphyanh4/(bkphyanh4 + ocetra(i,j,k,ianh4)) ! inhibition of NO3 uptake
nutlim = min(ocetra(i,j,k,iphosph)/(ocetra(i,j,k,iphosph)+bkphosph), &
ocetra(i,j,k,iiron)/(ocetra(i,j,k,iiron)+bkiron))
anh4lim = ocetra(i,j,k,ianh4)/(ocetra(i,j,k,ianh4) + bkphyanh4)
nlim = ano3up_inh*ocetra(i,j,k,iano3)/(ocetra(i,j,k,iano3) + bkphyano3) + anh4lim
grlim = min(nutlim,nlim) ! growth limitation
nh4uptfrac = anh4lim/(nlim+epsilon(1.))
! re-check avnut - can sum N avail exceed indiv. contrib?
avanut = max(0.,min(ocetra(i,j,k,iphosph), ocetra(i,j,k,iiron)/riron, &
& rnoi*((1.-nh4uptfrac)*ocetra(i,j,k,iano3) + nh4uptfrac*ocetra(i,j,k,ianh4))))
xn = avphy/(1. - pho*grlim) ! phytoplankton growth
phosy = max(0.,min(xn-avphy,max_limiter*avanut)) ! limit PP growth to available nutr.
else
avanut = max(0.,min(ocetra(i,j,k,iphosph),rnoi*ocetra(i,j,k,iano3)))
avanfe = max(0.,min(avanut,ocetra(i,j,k,iiron)/riron))
xa = avanfe
xn = xa/(1.+pho*avphy/(xa+bkphy))
phosy = max(0.,xa-xn)
endif
phosy = merge(avdic/rcar, phosy, avdic <= rcar*phosy) ! limit phosy by available DIC
ya = avphy+phosy
yn = (ya+grazra*avgra*phytomi/(avphy+bkzoo))/(1.+grazra*avgra/(avphy+bkzoo))
grazing = max(0.,ya-yn)
graton = epsher*(1.-zinges)*grazing
gratpoc = (1.-epsher)*grazing
grawa = epsher*zinges*grazing
phythresh = max(0.,(ocetra(i,j,k,iphy)-2.*phytomi))
phymor = dyphy*phythresh
zoothresh = max(0.,(ocetra(i,j,k,izoo)-2.*grami))
if (lkwrbioz_off) then
bacfra = 0.
else
bacfra = remido*ocetra(i,j,k,idoc)
endif
exud = gammap*phythresh
zoomor = spemor*zoothresh*zoothresh ! *10 compared to linear in tropics (tinka)
excdoc = gammaz*zoothresh ! excretion of doc by zooplankton
export = zoomor*(1.-ecan) + phymor + gratpoc ! ecan=.95, gratpoc= .2*grazing
if (use_cisonew) then
! calculation of isotope fractionation during photosynthesis (Laws 1997)
if(ocetra(i,j,k,iphy) < phytomi) then
bifr13 = 1.
else
phygrowth = ((ocetra(i,j,k,iphy)+phosy)/ocetra(i,j,k,iphy))/dtb ! Growth rate phytoplankton [1/d]
growth_co2 = phygrowth/(co2star(i,j,k)*1.e6+safediv) ! CO2* in [mol/kg]
bifr13_perm = (6.03 + 5.5*growth_co2)/(0.225 + growth_co2) ! Permil (~20)
bifr13_perm = max(5.,min(26.,bifr13_perm)) ! Limit the range to [5,26]
bifr13 = (1000. - bifr13_perm) / 1000. ! Fractionation factor 13c (~0.98)
endif
bifr14 = bifr13**2
! calculation of 13C and 14C equivalent of biology
rco213 = ocetra(i,j,k,isco213)/(ocetra(i,j,k,isco212)+safediv)
rco214 = ocetra(i,j,k,isco214)/(ocetra(i,j,k,isco212)+safediv)
rphy13 = ocetra(i,j,k,iphy13)/(ocetra(i,j,k,iphy)+safediv)
rphy14 = ocetra(i,j,k,iphy14)/(ocetra(i,j,k,iphy)+safediv)
rzoo13 = ocetra(i,j,k,izoo13)/(ocetra(i,j,k,izoo)+safediv)
rzoo14 = ocetra(i,j,k,izoo14)/(ocetra(i,j,k,izoo)+safediv)
phosy13 = phosy*bifr13*rco213
phosy14 = phosy*bifr14*rco214
grazing13 = grazing*rphy13
grazing14 = grazing*rphy14
graton13 = epsher*(1.-zinges)*grazing13
graton14 = epsher*(1.-zinges)*grazing14
gratpoc13 = (1.-epsher)*grazing13
gratpoc14 = (1.-epsher)*grazing14
grawa13 = epsher*zinges*grazing13
grawa14 = epsher*zinges*grazing14
if (lkwrbioz_off) then
bacfra13 = 0.
bacfra14 = 0.
else
bacfra13 = remido*ocetra(i,j,k,idoc13)
bacfra14 = remido*ocetra(i,j,k,idoc14)
endif
phymor13 = phymor*rphy13
phymor14 = phymor*rphy14
zoomor13 = zoomor*rzoo13
zoomor14 = zoomor*rzoo14
excdoc13 = excdoc*rzoo13
excdoc14 = excdoc*rzoo14
exud13 = exud*rphy13
exud14 = exud*rphy14
export13 = zoomor13*(1.-ecan) + phymor13 + gratpoc13
export14 = zoomor14*(1.-ecan) + phymor14 + gratpoc14
endif
if (use_AGG) then
delsil = min(ropal*phosy*avsil/(avsil+bkopal),0.5*avsil)
delcar = rcalc*min(calmax*phosy,(phosy-delsil/ropal))
! definition of delcar13/14 for the AGG scheme currently missing
else
delsil = min(ropal*export*avsil/(avsil+bkopal),0.5*avsil)
delcar = rcalc * export * bkopal/(avsil+bkopal)
if (use_cisonew) then
delcar13 = rcalc * export13 * bkopal/(avsil+bkopal)
delcar14 = rcalc * export14 * bkopal/(avsil+bkopal)
endif
endif
if(with_dmsph) then
dms_ph = 1. + (-log10(hi(i,j,1)) - pi_ph(i,j))*dms_gamma
else
dms_ph = 1.
endif
dmsprod = (dmsp5*delsil+dmsp4*delcar)*(1.+1./(temp+dmsp1)**2)*dms_ph
if (lkwrbioz_off) then
dms_bac = 0.
else
dms_bac = dmsp3*abs(temp+3.)*ocetra(i,j,k,idms) &
& *(ocetra(i,j,k,idms)/(dmsp6+ocetra(i,j,k,idms)))
endif
dms_uv = dmsp2*phofa/pi_alpha*ocetra(i,j,k,idms)
dtr = bacfra-phosy+graton+ecan*zoomor
ocetra(i,j,k,iphosph) = ocetra(i,j,k,iphosph)+dtr
if (.not. use_extNcycle) then
ocetra(i,j,k,iano3) = ocetra(i,j,k,iano3)+dtr*rnit
ocetra(i,j,k,ialkali) = ocetra(i,j,k,ialkali)-2.*delcar-(rnit+1)*dtr
ocetra(i,j,k,ioxygen) = ocetra(i,j,k,ioxygen)-dtr*ro2ut
else
ocetra(i,j,k,iano3) = ocetra(i,j,k,iano3) - (1.-nh4uptfrac)*phosy*rnit
ocetra(i,j,k,ianh4) = ocetra(i,j,k,ianh4) - nh4uptfrac*phosy*rnit + (dtr+phosy)*rnit
ocetra(i,j,k,ialkali) = ocetra(i,j,k,ialkali) - nh4uptfrac*phosy*(rnit-1.) & ! NH4 + PO4 Uptake
& + (1.-nh4uptfrac)*phosy*(rnit+1.) & ! NO3 + PO4 Uptake
& + (dtr+phosy)*(rnit-1.) - 2.*delcar ! Remin to (NH4 + PO4) and CaCO3 formation
ocetra(i,j,k,ioxygen) = ocetra(i,j,k,ioxygen) + nh4uptfrac*phosy*ro2utammo & ! NH4 uptake
& + (1.-nh4uptfrac)*phosy*ro2ut & ! NO3 uptake
& - (dtr+phosy)*ro2utammo ! Remin to NH4
! Output
phosy_NH4(i,j,k) = nh4uptfrac*phosy*rnit ! kmol N/m3/dtb - NH4 uptake during PP growth
phosy_NO3(i,j,k) = (1.-nh4uptfrac)*phosy*rnit ! kmol N/m3/dtb - NO3 uptake during PP growth
remin_aerob(i,j,k) = (dtr+phosy)*rnit ! kmol N/m3/dtb - Aerob remin to ammonium (var. sources)
endif
ocetra(i,j,k,idet) = ocetra(i,j,k,idet)+export
ocetra(i,j,k,idms) = ocetra(i,j,k,idms)+dmsprod-dms_bac-dms_uv
ocetra(i,j,k,isco212) = ocetra(i,j,k,isco212)-delcar+rcar*dtr
ocetra(i,j,k,iphy) = ocetra(i,j,k,iphy)+phosy-grazing-phymor-exud
ocetra(i,j,k,izoo) = ocetra(i,j,k,izoo)+grawa-excdoc-zoomor
ocetra(i,j,k,idoc) = ocetra(i,j,k,idoc)-bacfra+excdoc+exud
ocetra(i,j,k,icalc) = ocetra(i,j,k,icalc)+delcar
if (use_cisonew) then
dtr13 = bacfra13-phosy13+graton13+ecan*zoomor13
dtr14 = bacfra14-phosy14+graton14+ecan*zoomor14
ocetra(i,j,k,idet13) = ocetra(i,j,k,idet13)+export13
ocetra(i,j,k,idet14) = ocetra(i,j,k,idet14)+export14
ocetra(i,j,k,isco213) = ocetra(i,j,k,isco213)-delcar13+rcar*dtr13
ocetra(i,j,k,isco214) = ocetra(i,j,k,isco214)-delcar14+rcar*dtr14
ocetra(i,j,k,iphy13) = ocetra(i,j,k,iphy13)+phosy13-grazing13-phymor13-exud13
ocetra(i,j,k,iphy14) = ocetra(i,j,k,iphy14)+phosy14-grazing14-phymor14-exud14
ocetra(i,j,k,izoo13) = ocetra(i,j,k,izoo13)+grawa13-excdoc13-zoomor13
ocetra(i,j,k,izoo14) = ocetra(i,j,k,izoo14)+grawa14-excdoc14-zoomor14
ocetra(i,j,k,idoc13) = ocetra(i,j,k,idoc13)-bacfra13+excdoc13+exud13
ocetra(i,j,k,idoc14) = ocetra(i,j,k,idoc14)-bacfra14+excdoc14+exud14
ocetra(i,j,k,icalc13) = ocetra(i,j,k,icalc13)+delcar13
ocetra(i,j,k,icalc14) = ocetra(i,j,k,icalc14)+delcar14
endif
if (use_natDIC) then
ocetra(i,j,k,inatsco212) = ocetra(i,j,k,inatsco212)-delcar+rcar*dtr
ocetra(i,j,k,inatalkali) = ocetra(i,j,k,inatalkali)-2.*delcar-(rnit+1)*dtr
ocetra(i,j,k,inatcalc) = ocetra(i,j,k,inatcalc)+delcar
endif
if (lkwrbioz_off) then
opalrem = 0.
else
if (use_M4AGO) then
opalrem = dremopal*opal_remin_q10**((ptho(i,j,k)-opal_remin_Tref)/10.)*ocetra(i,j,k,iopal)
else
opalrem = dremopal*ocetra(i,j,k,iopal)
endif
endif
ocetra(i,j,k,isilica) = ocetra(i,j,k,isilica)-delsil+opalrem
ocetra(i,j,k,iopal) = ocetra(i,j,k,iopal)+delsil-opalrem
ocetra(i,j,k,iiron) = ocetra(i,j,k,iiron)+dtr*riron &
& - relaxfe*max(ocetra(i,j,k,iiron)-fesoly,0.)
if (use_BROMO) then
! Bromo source from phytoplankton production and sink to photolysis
! Hense and Quack (200) Pg537 Decay time scale is 30days =0.0333/day
! sinks owing to degradation by nitrifiers (Pg 538 of Hense and Quack,
! 2009) is omitted because the magnitude is more than 2 order smaller
! than sink through halide substitution & hydrolysis (Fig. 3)
! Assume that only 30% of incoming radiation are UV (i.e. 50% of non-PAR
! radiation; PAR radiationis assume to be 40% of incoming radiation)
bro_beta = rbro*(fbro1*avsil/(avsil+bkopal)+fbro2*bkopal/(avsil+bkopal))
if (swa_clim(i,j,1) > 0.) then
bro_uv = 0.0333*dtb*0.3*(strahl(i,j)/swa_clim(i,j,1))*abs_uv(i,j,k)*ocetra(i,j,k,ibromo)
else
bro_uv = 0.0
endif
ocetra(i,j,k,ibromo) = ocetra(i,j,k,ibromo)+bro_beta*phosy-bro_uv
endif
if (use_AGG) then
!***********************************************************************
! effects of biological processes on number of particles:
! photosynthesis creates POM
! exudation deletes POM
! grazing deletes POM; but only the fraction that is not egested as
! fecal pellets again (grawa remains in zoo, graton goes to po4)
! none of the processes at the current time is assumed to change
! the size distribution (subject to change)
! NOTE that phosy, exud etc. are in kmol/m3!
! Thus divide by avmass (kmol/m3)
!**********************************************************************
if(avmass > 0.) then
avnos = ocetra(i,j,k,inos)
anosloss = (phosy-exud-graton-grawa)*avnos/avmass
ocetra(i,j,k,inos) = ocetra(i,j,k,inos)+anosloss
endif
!***********************************************************************
! dead zooplankton corpses come with their own, flat distribution
! this flow even takes place if there is neither nos nor mass
! NOTE: zoomor is in kmol/m3!! Thus multiply flow by 1.e+6
!***********************************************************************
zmornos = zoomor * (1.-ecan) * zdis * 1.e+6
ocetra(i,j,k,inos) = ocetra(i,j,k,inos)+zmornos
endif
! add up for total inventory and output
dz = pddpo(i,j,k)
expoor(i,j) = expoor(i,j) +export*rcar*dz
expoca(i,j) = expoca(i,j) +delcar*dz
exposi(i,j) = exposi(i,j) +delsil*dz
intdmsprod(i,j) = intdmsprod(i,j)+dmsprod*dz
intdms_bac(i,j) = intdms_bac(i,j)+dms_bac*dz
intdms_uv(i,j) = intdms_uv (i,j)+dms_uv*dz
if (use_BROMO) then
int_chbr3_uv(i,j) = int_chbr3_uv (i,j) + bro_uv*dz
int_chbr3_prod(i,j) = int_chbr3_prod (i,j) + bro_beta*phosy*dz
endif
intphosy(i,j) = intphosy(i,j) +phosy*rcar*dz ! primary production in kmol C m-2
phosy3d(i,j,k) = phosy*rcar ! primary production in kmol C m-3
endif ! pddpo(i,j,k) > dp_min
enddo ! kwrbioz
enddo ! kpie
enddo loop1 ! kpje
!$OMP END PARALLEL DO
if (use_PBGC_OCNP_TIMESTEP) then
if (mnproc == 1) then
write(io_stdo_bgc,*)' '
write(io_stdo_bgc,*)'in OCRPOD after 1st bio prod'
endif
call inventory_bgc(kpie,kpje,kpke,pdlxp,pdlyp,pddpo,omask,0)
endif
!$OMP PARALLEL DO PRIVATE(phythresh,zoothresh,sterph,sterzo,remin &
!$OMP ,opalrem,aou,refra,dms_bac,pocrem,docrem,phyrem,dz,o2lim &
!$OMP ,avmass,avnos,zmornos &
!$OMP ,rphy13,rphy14,rzoo13,rzoo14,rdet13,rdet14,rdoc13,rdoc14 &
!$OMP ,sterph13,sterph14,sterzo13,sterzo14,pocrem13,pocrem14 &
!$OMP ,docrem13,docrem14,phyrem13,phyrem14 &
!$OMP ,i,k)
loop2: do j = 1,kpje
do i = 1,kpie
do k = merge(1,kwrbioz(i,j)+1,lkwrbioz_off),kpke
if(pddpo(i,j,k) > dp_min .and. omask(i,j) > 0.5) then
if (use_AGG) then
avmass = ocetra(i,j,k,iphy)+ocetra(i,j,k,idet)
endif
temp = min(40.,max(-3.,ptho(i,j,k)))
phythresh = max(0.,(ocetra(i,j,k,iphy)-2.*phytomi))
zoothresh = max(0.,(ocetra(i,j,k,izoo)-2.*grami))
sterph = 0.5*dyphy*phythresh ! phytoplankton to detritus
sterzo = spemor*zoothresh*zoothresh ! quadratic mortality
if (use_cisonew) then
rphy13 = ocetra(i,j,k,iphy13)/(ocetra(i,j,k,iphy)+safediv)
rphy14 = ocetra(i,j,k,iphy14)/(ocetra(i,j,k,iphy)+safediv)
rzoo13 = ocetra(i,j,k,izoo13)/(ocetra(i,j,k,izoo)+safediv)
rzoo14 = ocetra(i,j,k,izoo14)/(ocetra(i,j,k,izoo)+safediv)
rdet13 = ocetra(i,j,k,idet13)/(ocetra(i,j,k,idet)+safediv)
rdet14 = ocetra(i,j,k,idet14)/(ocetra(i,j,k,idet)+safediv)
rdoc13 = ocetra(i,j,k,idoc13)/(ocetra(i,j,k,idoc)+safediv)
rdoc14 = ocetra(i,j,k,idoc14)/(ocetra(i,j,k,idoc)+safediv)
sterph13 = sterph*rphy13
sterph14 = sterph*rphy14
sterzo13 = sterzo*rzoo13
sterzo14 = sterzo*rzoo14
endif
if (lkwrbioz_off) then ! dying before in PP loop
sterph = 0.
sterzo = 0.
sterph13 = 0.
sterph14 = 0.
sterzo13 = 0.
sterzo14 = 0.
endif
ocetra(i,j,k,iphy) = ocetra(i,j,k,iphy)-sterph
ocetra(i,j,k,izoo) = ocetra(i,j,k,izoo)-sterzo
if (use_cisonew) then
ocetra(i,j,k,iphy13) = ocetra(i,j,k,iphy13)-sterph13
ocetra(i,j,k,iphy14) = ocetra(i,j,k,iphy14)-sterph14
ocetra(i,j,k,izoo13) = ocetra(i,j,k,izoo13)-sterzo13
ocetra(i,j,k,izoo14) = ocetra(i,j,k,izoo14)-sterzo14
endif
if(ocetra(i,j,k,ioxygen) > O2thresh_aerob) then
if (lTO2depremin) then
! Both, use_M4AGO and use_extNcycle switch lTO2depremin to true!
o2lim = ocetra(i,j,k,ioxygen)/(ocetra(i,j,k,ioxygen) + bkox_drempoc)
pocrem = drempoc*o2lim*POM_remin_q10**((ptho(i,j,k)-POM_remin_Tref)/10.)*ocetra(i,j,k,idet)
else
pocrem = drempoc*ocetra(i,j,k,idet)
endif
if (.not. use_extNcycle) then
pocrem = min(pocrem, 0.33*ocetra(i,j,k,ioxygen)/ro2ut)
docrem = min(remido*ocetra(i,j,k,idoc), 0.33*ocetra(i,j,k,ioxygen)/ro2ut)
phyrem = min(0.5*dyphy*phythresh, 0.33*ocetra(i,j,k,ioxygen)/ro2ut)
else
pocrem = min(pocrem, 0.33*ocetra(i,j,k,ioxygen)/ro2utammo)
docrem = min(remido*ocetra(i,j,k,idoc), 0.33*ocetra(i,j,k,ioxygen)/ro2utammo)
phyrem = min(0.5*dyphy*phythresh, 0.33*ocetra(i,j,k,ioxygen)/ro2utammo)
endif
if (lkwrbioz_off) then ! dying before in PP loop
phyrem = 0.
endif
if (use_cisonew) then
pocrem13 = pocrem*rdet13
pocrem14 = pocrem*rdet14
docrem13 = docrem*rdoc13
docrem14 = docrem*rdoc14
phyrem13 = phyrem*rphy13
phyrem14 = phyrem*rphy14
endif
else
pocrem = 0.
docrem = 0.
phyrem = 0.
if (use_cisonew) then
pocrem13 = 0.
docrem13 = 0.
phyrem13 = 0.
pocrem14 = 0.
docrem14 = 0.
phyrem14 = 0.
endif
endif
ocetra(i,j,k,idet) = ocetra(i,j,k,idet) - pocrem + sterph + sterzo
ocetra(i,j,k,idoc) = ocetra(i,j,k,idoc) - docrem
ocetra(i,j,k,iphy) = ocetra(i,j,k,iphy) - phyrem
remin = pocrem + docrem + phyrem
ocetra(i,j,k,iphosph) = ocetra(i,j,k,iphosph)+remin
if (.not. use_extNcycle) then
ocetra(i,j,k,iano3) = ocetra(i,j,k,iano3)+remin*rnit
ocetra(i,j,k,ialkali) = ocetra(i,j,k,ialkali)-(rnit+1)*remin
ocetra(i,j,k,ioxygen) = ocetra(i,j,k,ioxygen)-ro2ut*remin
else
ocetra(i,j,k,ianh4) = ocetra(i,j,k,ianh4) + remin*rnit
ocetra(i,j,k,ialkali) = ocetra(i,j,k,ialkali) + (rnit-1.)*remin
ocetra(i,j,k,ioxygen) = ocetra(i,j,k,ioxygen) - ro2utammo*remin
remin_aerob(i,j,k) = remin_aerob(i,j,k)+remin*rnit ! kmol/NH4/dtb - remin to NH4 from various sources
endif
ocetra(i,j,k,isco212) = ocetra(i,j,k,isco212)+rcar*remin
ocetra(i,j,k,iiron) = ocetra(i,j,k,iiron)+remin*riron &
& -relaxfe*max(ocetra(i,j,k,iiron)-fesoly,0.)
if (use_natDIC) then
ocetra(i,j,k,inatsco212) = ocetra(i,j,k,inatsco212)+rcar*remin
ocetra(i,j,k,inatalkali) = ocetra(i,j,k,inatalkali)-(rnit+1)*remin
endif
if (use_cisonew) then
ocetra(i,j,k,idet13) = ocetra(i,j,k,idet13)-pocrem13+sterph13+sterzo13
ocetra(i,j,k,idet14) = ocetra(i,j,k,idet14)-pocrem14+sterph14+sterzo14
ocetra(i,j,k,idoc13) = ocetra(i,j,k,idoc13)-docrem13
ocetra(i,j,k,idoc14) = ocetra(i,j,k,idoc14)-docrem14
ocetra(i,j,k,iphy13) = ocetra(i,j,k,iphy13)-phyrem13
ocetra(i,j,k,iphy14) = ocetra(i,j,k,iphy14)-phyrem14
ocetra(i,j,k,isco213) = ocetra(i,j,k,isco213)+rcar*(pocrem13+docrem13+phyrem13)
ocetra(i,j,k,isco214) = ocetra(i,j,k,isco214)+rcar*(pocrem14+docrem14+phyrem14)
endif
!***********************************************************************
! as ragueneau (2000) notes, Si(OH)4sat is about 1000 umol, but
! Si(OH)4 varies only between 0-100 umol
! so the expression dremopal*(Si(OH)4sat-Si(OH)4) would change the
! rate only from 0 to 100%
!***********************************************************************
if (use_M4AGO) then
opalrem = dremopal*opal_remin_q10**((ptho(i,j,k)-opal_remin_Tref)/10.)*ocetra(i,j,k,iopal)
else
opalrem = dremopal*0.1*(temp+3.)*ocetra(i,j,k,iopal)
endif
ocetra(i,j,k,iopal) = ocetra(i,j,k,iopal)-opalrem
ocetra(i,j,k,isilica) = ocetra(i,j,k,isilica)+opalrem
if (.not. use_extNcycle) then
!***********************************************************************
! There is about 1.e4 O2 on 1 N2O molecule (Broeker&Peng)
! refra : Tim Rixton, private communication
!***********************************************************************
aou = satoxy(i,j,k)-ocetra(i,j,k,ioxygen)
refra = 1.+3.*(0.5+sign(0.5,aou-1.97e-4))
ocetra(i,j,k,ian2o) = ocetra(i,j,k,ian2o)+remin*1.e-4*ro2ut*refra
ocetra(i,j,k,igasnit) = ocetra(i,j,k,igasnit)-remin*1.e-4*ro2ut*refra
ocetra(i,j,k,ioxygen) = ocetra(i,j,k,ioxygen)-remin*1.e-4*ro2ut*refra*0.5
endif
dms_bac = dmsp3 * abs(temp+3.) * ocetra(i,j,k,idms) &
& * (ocetra(i,j,k,idms) / (dmsp6+ocetra(i,j,k,idms)))
ocetra(i,j,k,idms) = ocetra(i,j,k,idms)-dms_bac
dz = pddpo(i,j,k)
intdms_bac(i,j) = intdms_bac(i,j)+dms_bac*dz
if (use_AGG) then
!***********************************************************************
! loss of snow numbers due to remineralization of poc
! gain of snow numbers due to zooplankton mortality
! NOTE that remin is in kmol/m3. Thus divide by avmass (kmol/m3)
!***********************************************************************
if(avmass > 0.) then
avnos = ocetra(i,j,k,inos)
ocetra(i,j,k,inos) = ocetra(i,j,k,inos)-remin*avnos/avmass
endif
!***********************************************************************
! dead zooplankton corpses come with their own, flat distribution
! this flow even takes place if there is neither nos nor mass
! NOTE: zoomor is in kmol/m3!! Thus multiply flow by 1.e+6
!***********************************************************************
zmornos = sterzo * zdis * 1.e+6
ocetra(i,j,k,inos) = ocetra(i,j,k,inos) + zmornos
endif/*AGG*/
endif
enddo
enddo
enddo loop2
!$OMP END PARALLEL DO
if (use_PBGC_OCNP_TIMESTEP) then
if (mnproc == 1) then
write(io_stdo_bgc,*)' '
write(io_stdo_bgc,*)'in OCRPOD after poc remin'
endif
call inventory_bgc(kpie,kpje,kpke,pdlxp,pdlyp,pddpo,omask,0)
endif
if (.not. use_extNcycle) then
! =====>>>> Regular CMIP6 iHAMOCC version for denitrification wo extended nitrogen cycle =====>>>>
!$OMP PARALLEL DO PRIVATE(remin,remin2o,dz,avmass,avnos,rem13,rem14,i,k)
loop3: do j = 1,kpje
do i = 1,kpie
do k = merge(1,kwrbioz(i,j)+1,lkwrbioz_off),kpke
if(omask(i,j) > 0.5) then
if(ocetra(i,j,k,ioxygen) < O2thresh_hypoxic .and. pddpo(i,j,k) > dp_min) then
if (use_AGG) then
avmass = ocetra(i,j,k,iphy) + ocetra(i,j,k,idet)
endif
remin = drempoc_anaerob*min(ocetra(i,j,k,idet),0.5 *ocetra(i,j,k,iano3)/rdnit1)
remin2o = dremn2o*min(ocetra(i,j,k,idet),0.003 *ocetra(i,j,k,ian2o)/rdn2o1)
if (use_cisonew) then
rem13 = (remin+remin2o)*ocetra(i,j,k,idet13)/(ocetra(i,j,k,idet)+safediv)
rem14 = (remin+remin2o)*ocetra(i,j,k,idet14)/(ocetra(i,j,k,idet)+safediv)
endif
ocetra(i,j,k,ialkali) = ocetra(i,j,k,ialkali)+(rdnit1-1)*remin-remin2o
ocetra(i,j,k,isco212) = ocetra(i,j,k,isco212)+rcar*(remin+remin2o)
ocetra(i,j,k,idet) = ocetra(i,j,k,idet)-(remin+remin2o)
ocetra(i,j,k,iphosph) = ocetra(i,j,k,iphosph)+(remin+remin2o)
ocetra(i,j,k,iano3) = ocetra(i,j,k,iano3)-rdnit1*remin
ocetra(i,j,k,igasnit) = ocetra(i,j,k,igasnit)+rdnit2*remin+rdn2o2*remin2o
ocetra(i,j,k,iiron) = ocetra(i,j,k,iiron)+riron*(remin+remin2o)
ocetra(i,j,k,ian2o) = ocetra(i,j,k,ian2o)-rdn2o1*remin2o
if (use_natDIC) then
ocetra(i,j,k,inatalkali) = ocetra(i,j,k,inatalkali)+(rdnit1-1)*remin-remin2o
ocetra(i,j,k,inatsco212) = ocetra(i,j,k,inatsco212)+rcar*(remin+remin2o)
endif
if (use_cisonew) then
ocetra(i,j,k,isco213) = ocetra(i,j,k,isco213)+rcar*rem13
ocetra(i,j,k,isco214) = ocetra(i,j,k,isco214)+rcar*rem14
ocetra(i,j,k,idet13) = ocetra(i,j,k,idet13)-rem13
ocetra(i,j,k,idet14) = ocetra(i,j,k,idet14)-rem14
endif
! nitrate loss through denitrification in kmol N m-2
dz = pddpo(i,j,k)
intdnit(i,j) = intdnit(i,j) + rdnit0*remin*dz
if (use_AGG) then
!***********************************************************************
! loss of snow numbers due to remineralization of poc
! NOTE that remin is in kmol/m3. Thus divide by avmass (kmol/m3)
!***********************************************************************
if(avmass > 0.) then
avnos = ocetra(i,j,k,inos)
ocetra(i,j,k,inos) = ocetra(i,j,k,inos)-(remin+remin2o)*avnos/avmass
endif
endif/*AGG*/
endif
endif
enddo
enddo
enddo loop3
!$OMP END PARALLEL DO
if (use_PBGC_OCNP_TIMESTEP) then
if (mnproc == 1) then
write(io_stdo_bgc,*)' '
write(io_stdo_bgc,*)'in OCRPOD after remin n2o'
endif
call inventory_bgc(kpie,kpje,kpke,pdlxp,pdlyp,pddpo,omask,0)
endif
! <<<<===== end of CMIP6 version denitrification processes without extended nitrogen cycle <<<<=====
else
!======>>>> extended nitrogen cycle processes (aerobic and anaerobic) that follow ammonification
inv_message = 'in OCPROD after extNcycle nitrification'
call nitrification(kpie,kpje,kpke,kbnd,pddpo,omask,ptho)
call extN_inv_check(kpie,kpje,kpke,pdlxp,pdlyp,pddpo,omask,inv_message)
inv_message = 'in OCPROD after extNcycle denitrification NO3 -> NO2'
call denit_NO3_to_NO2(kpie,kpje,kpke,kbnd,pddpo,omask,ptho)
call extN_inv_check(kpie,kpje,kpke,pdlxp,pdlyp,pddpo,omask,inv_message)
inv_message = 'in OCPROD after extNcycle anammox'
call anammox(kpie,kpje,kpke,kbnd,pddpo,omask,ptho)
call extN_inv_check(kpie,kpje,kpke,pdlxp,pdlyp,pddpo,omask,inv_message)
inv_message = 'in OCPROD after extNcycle denitrification / DNRA'
call denit_dnra(kpie,kpje,kpke,kbnd,pddpo,omask,ptho)
call extN_inv_check(kpie,kpje,kpke,pdlxp,pdlyp,pddpo,omask,inv_message)
endif
!sulphate reduction ! introduced 11.5.2007 to improve poc-remineralisation in the
! oxygen minimum zone in the subsurface equatorial Pacific
! assumption of endless pool of SO4 (typical concentration are on the order of mmol/l)
! js 02072007: for other runs than current millenium (cosmos-setup) experiments this seems
! to cause trouble as phosphate concentrations are too high at the depth of the oxygen
! minimum in the equatorial pacific/atlantic
! does it make sense to check for oxygen and nitrate deficit?
!$OMP PARALLEL DO PRIVATE(remin,avmass,avnos,rem13,rem14,i,k)
loop4: do j = 1,kpje
do i = 1,kpie
do k = merge(1,kwrbioz(i,j)+1,lkwrbioz_off),kpke
if(omask(i,j) > 0.5 .and. pddpo(i,j,k) > dp_min) then
if(ocetra(i,j,k,ioxygen) < O2thresh_hypoxic .and. ocetra(i,j,k,iano3) < NO3thresh_sulf ) then
if (use_AGG) then
avmass = ocetra(i,j,k,iphy)+ocetra(i,j,k,idet)
endif
remin = dremsul*ocetra(i,j,k,idet)
if (use_cisonew) then
rem13 = remin*ocetra(i,j,k,idet13)/(ocetra(i,j,k,idet)+safediv)
rem14 = remin*ocetra(i,j,k,idet14)/(ocetra(i,j,k,idet)+safediv)
endif
ocetra(i,j,k,idet) = ocetra(i,j,k,idet)-remin
ocetra(i,j,k,ialkali) = ocetra(i,j,k,ialkali)-(rnit+1)*remin
ocetra(i,j,k,isco212) = ocetra(i,j,k,isco212)+rcar*remin
ocetra(i,j,k,iphosph) = ocetra(i,j,k,iphosph)+remin
ocetra(i,j,k,iano3) = ocetra(i,j,k,iano3)+rnit*remin
ocetra(i,j,k,iiron) = ocetra(i,j,k,iiron)+riron*remin
if (use_natDIC) then
ocetra(i,j,k,inatalkali) = ocetra(i,j,k,inatalkali)-(rnit+1)*remin
ocetra(i,j,k,inatsco212) = ocetra(i,j,k,inatsco212)+rcar*remin
endif
if (use_cisonew) then
ocetra(i,j,k,idet13) = ocetra(i,j,k,idet13)-rem13
ocetra(i,j,k,idet14) = ocetra(i,j,k,idet14)-rem14
ocetra(i,j,k,isco213) = ocetra(i,j,k,isco213)+rcar*rem13
ocetra(i,j,k,isco214) = ocetra(i,j,k,isco214)+rcar*rem14
endif
if (use_extNcycle) then
! Output
remin_sulf(i,j,k) = remin ! kmol P/m3/dtb
endif
if (use_AGG) then
!***********************************************************************
! loss of snow numbers due to remineralization of poc
! NOTE that remin is in kmol/m3. Thus divide by avmass (kmol/m3)
!***********************************************************************
if(avmass > 0.) then
avnos = ocetra(i,j,k,inos)
ocetra(i,j,k,inos) = ocetra(i,j,k,inos)-remin*avnos/avmass
endif
endif
endif
endif
enddo
enddo
enddo loop4
!$OMP END PARALLEL DO
! end sulphate reduction
if (use_PBGC_OCNP_TIMESTEP) then
if (mnproc == 1) then
write(io_stdo_bgc,*)' '
write(io_stdo_bgc,*)'in OCRPOD after sulphate reduction '
endif
call inventory_bgc(kpie,kpje,kpke,pdlxp,pdlyp,pddpo,omask,0)
endif
if (use_AGG) then
!**********************AGGREGATION***************************************
! General:
! Sinking speed, size distribution and aggregation are calculated
! as in Kriest and Evans, 2000. I assume that opal and calcium carbonate
! sink at the same speed as P (mass).
!
! Sinking speed and aggregation: I assume that if there is no phosphorous mass,
! the sinking speed is the minimum sinking speed of aggregates. I further
! assume that then there are no particles, and that the rate of aggregation
! is 0. This scheme removes no P in the absence of P, but still opal and/or
! calcium carbonate.
! This could or should be changed, because silica as well as carbonate
! shell will add to the aggregate mass, and should be considered.
! Puh. Does anyone know functional relationships between
! size and Si or CaCO3? Perhaps on a later version, I have to
! take the relationship bewteen weight and size?
!
! Size distribution and resulting loss of marine snow aggregates due to
! aggregation (aggregate(i,j,k)) and sinking speed of mass and numbers
! (wmass(i,j,k) and wnumb(i,j,k) are calculated in a loop over 2-kpke.
!
!************************************************************************
wmass(:,:,:) = 0.0
wnumb(:,:,:) = 0.0