pkg/bling/bling_dvm.F

C $Header: /u/gcmpack/MITgcm/pkg/bling/bling_dvm.F,v 1.5 2016/09/12 20:00:28 mmazloff Exp $
C $Name:  $

#include "BLING_OPTIONS.h"

CBOP
      subroutine BLING_DVM(
     I           N_dvm,P_dvm,Fe_dvm,
     I           PTR_O2, mld,
     O           N_remindvm, P_remindvm, Fe_remindvm,
     I           bi, bj, imin, imax, jmin, jmax,
     I           myIter, myTime, myThid )

C     =================================================================
C     | subroutine bling_dvm
C     | o Diel Vertical Migration
C     =================================================================

      implicit none

C     === Global variables ===

#include "SIZE.h"
#include "DYNVARS.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"
#include "BLING_VARS.h"
#include "PTRACERS_SIZE.h"
#include "PTRACERS_PARAMS.h"
#ifdef ALLOW_AUTODIFF
# include "tamc.h"
#endif

C     === Routine arguments ===
C     bi,bj         :: tile indices
C     iMin,iMax     :: computation domain: 1rst index range
C     jMin,jMax     :: computation domain: 2nd  index range
C     myTime        :: current time
C     myIter        :: current timestep
C     myThid        :: thread Id. number
      INTEGER bi, bj, imin, imax, jmin, jmax
      _RL     myTime
      INTEGER myIter
      INTEGER myThid
C     === Input ===
C     N_dvm         :: vertical transport of nitrogen by DVM
C     P_dvm         :: vertical transport of phosphorus by DVM
C     Fe_dvm        :: vertical transport of iron by DVM
C     PTR_O2        :: nitrate concentration
C     mld           :: mixed layer depth
      _RL     N_dvm (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
      _RL     P_dvm (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
      _RL     Fe_dvm(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
      _RL     PTR_O2(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
      _RL     mld   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      
C     === Output ===
C     N_remindvm    :: nitrogen remineralization due to diel vertical migration
C     P_remindvm    :: phosphorus remineralization due to diel vertical migration
C     Fe_remindvm   :: iron remineralization due to diel vertical migration
      _RL     N_remindvm     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
      _RL     P_remindvm     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
      _RL     Fe_remindvm     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)

#ifdef ALLOW_BLING
C     === Local variables ===
C     i,j,k         :: loop indices
      INTEGER i,j,k
      INTEGER tmp
      _RL depth_l
      _RL o2_upper
      _RL o2_lower
      _RL dz_upper
      _RL dz_lower
      _RL temp_upper
      _RL temp_lower
      _RL z_dvm_regr
      _RL frac_migr
      _RL fdvm_migr
      _RL fdvm_stat
      _RL fdvmn_vint
      _RL fdvmp_vint
      _RL fdvmfe_vint
      _RL z_dvm
      _RL dvm(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
      _RL Fe_burial(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL x_erfcc,z_erfcc,t_erfcc,erfcc
CEOP

c ---------------------------------------------------------------------
c  Initialize output and diagnostics
       DO k=1,Nr
        DO j=jmin,jmax
          DO i=imin,imax
              N_remindvm(i,j,k)        = 0. _d 0
              P_remindvm(i,j,k)        = 0. _d 0
              Fe_remindvm(i,j,k)       = 0. _d 0
              dvm(i,j,k)               = 0. _d 0
          ENDDO
          Fe_burial(i,j)               = 0. _d 0
        ENDDO
       ENDDO
       
C  ---------------------------------------------------------------------
c DIEL VERTICAL MIGRATOR EXPORT
c The effect of vertically-migrating animals on the export flux of organic
c matter from the ocean surface is treated similarly to the scheme of
c Bianchi et al., Nature Geoscience 2013.
c This involves calculating the stationary depth of vertical migrators, using
c an empirical multivariate regression, and ensuring that this remains
c above the bottom as well as any suboxic waters.
c The total DVM export flux is partitioned between a swimming migratory
c component and the stationary component, and these are summed.

C$TAF LOOP = parallel
      DO j=jmin,jmax
C$TAF LOOP = parallel
       DO i=imin,imax

c  Initialize 
        o2_upper = 0.
        o2_lower = 0.
        dz_upper = 0.
        dz_lower = 0.
        temp_upper = 0.
        temp_lower = 0.
        z_dvm_regr = 0.
        z_dvm     = 0.
        frac_migr = 0.
        fdvm_migr = 0.
        fdvm_stat = 0.
        fdvmn_vint = 0.
        fdvmp_vint = 0.
        fdvmfe_vint = 0.

        DO k=1,Nr

         IF ( hFacC(i,j,k,bi,bj).gt.0. _d 0 ) THEN

c  Calculate the depth of migration based on linear regression.

        depth_l=-rF(k+1)

c  Average temperature and oxygen over upper 35 m, and 140-515m.
c  Also convert O2 to mmol m-3.

          if ( abs(depth_l) .lt. 35.) then
            dz_upper = dz_upper + drf(k)
            temp_upper = temp_upper + theta(i,j,k,bi,bj)*drf(k)
            o2_upper = o2_upper + PTR_O2(i,j,k) * drf(k)*1.0 _d 3
          endif
          if ( (abs(depth_l) .gt. 140.0 _d 0) .and.
     &          (abs(depth_l) .lt. 515. _d 0)) then
            dz_lower = dz_lower + drf(k)
            temp_lower = temp_lower + theta(i,j,k,bi,bj)*drf(k)
            o2_lower = o2_lower + PTR_O2(i,j,k) * drf(k)*1.0 _d 3
          endif

         ENDIF
        ENDDO

        o2_upper = o2_upper / (epsln + dz_upper)
        temp_upper = temp_upper / (epsln + dz_upper)
        o2_lower = o2_lower / (epsln + dz_lower)
        temp_lower = temp_lower / (epsln + dz_lower)

c  Calculate the regression, using the constants given in Bianchi et al. (2013).
c  The variable values are bounded to lie within reasonable ranges: 
c         O2 gradient   : [-10,300] mmol/m3
c         Log10 Chl     : [-1.8,0.85] log10(mg/m3)
c         mld           : [0,500] m
c         T gradient    : [-3,20] C

c!! replacing hblt_depth(i,j) with mld... 

#ifdef BLING_ADJOINT_SAFE

        z_dvm = 300. _d 0

#else

        z_dvm_regr = 398. _d 0
     &   - 0.56 _d 0*min(300. _d 0,max(-10. _d 0,(o2_upper - o2_lower)))
     &   - 115. _d 0*min(0.85 _d 0,max(-1.80 _d 0,
     &   log10(max(chl(i,j,1,bi,bj),chl_min))))
     &   + 0.36 _d 0*min(500. _d 0,max(epsln,mld(i,j)))
     &   - 2.40 _d 0*min(20. _d 0,max(-3. _d 0,(temp_upper-temp_lower)))

c  Limit the depth of migration in polar winter.
c  Use irr_mem since this is averaged over multiple days, dampening the 
c  diurnal cycle.
c  Tapers Z_DVM to the minimum when surface irradince is below a given  
c  threshold (here 10 W/m2).

        if ( irr_mem(i,j,1,bi,bj) .lt. 10. ) then
          z_dvm_regr = 150. _d 0 + (z_dvm_regr - 150. _d 0) *
     &       irr_mem(i,j,1,bi,bj) / 10. _d 0
        endif

c  Check for suboxic water within the column. If found, set dvm
c  stationary depth to 2 layers above it. This is not meant to
c  represent a cessation of downward migration, but rather the
c  requirement for aerobic DVM respiration to occur above the suboxic
c  water, where O2 is available.

        tmp = 0
        DO k=1,Nr-2

         IF ( (hFacC(i,j,k,bi,bj).gt.0. _d 0) .and. (tmp.eq.0)) THEN

          z_dvm = -rf(k+1)
          if (PTR_O2(i,j,k+2) .lt. (5. _d 0*oxic_min)) tmp = 1

         ENDIF

        enddo

c  The stationary depth is constrained between 150 and 700, above any
c  anoxic waters found, and above the bottom.

        z_dvm = min(700. _d 0,max(150. _d 0,z_dvm_regr),z_dvm,-rf(k+1))
c!!           bling%zbot(i,j,grid_kmt(i,j))) * grid_tmask(i,j,1)
c!! what is grid_kmt?

#endif

c  Calculate the fraction of migratory respiration that occurs
c  during upwards and downwards swimming. The remainder is
c  respired near the stationary depth.
c  Constants for swimming speed and resting time are hard-coded
c  after Bianchi et al, Nature Geoscience 2013.

        frac_migr = max( 0.0 _d 0, min( 1.0 _d 0, (2.0 _d 0 * z_dvm) /
     &        (epsln + 0.05 _d 0 * 0.5 _d 0 * 86400. _d 0)))

c  Calculate the vertical profile shapes of DVM fluxes.
c  These are given as the downward organic flux due to migratory
c  DVM remineralization, defined at the bottom of each layer k.

        tmp = 0
        DO k=1,Nr

         IF ( (hFacC(i,j,k,bi,bj).gt.0. _d 0) .and. (tmp.eq.0)) THEN

          ! First, calculate the part due to active migration above
          ! the stationary depth.
          if (-rf(k+1) .lt. z_dvm) then
            fdvm_migr = frac_migr / (epsln + z_dvm - (-rf(2))) *
     &            (z_dvm - (-rf(k+1)) )
          else
            fdvm_migr  = 0.0
          endif

c  Then, calculate the part at the stationary depth.

c  Approximation of the complementary error function
c  From Numerical Recipes (F90, Ch. 6, p. 216)
c  Returns the complementary error function erfc(x)
c  with fractional error everywhere less than 1.2e-7
           x_erfcc = (-rf(k) - z_dvm) /
     &        ( (epsln + 2. _d 0 * sigma_dvm**2. _d 0)**0.5)

           z_erfcc = abs(x_erfcc)

           t_erfcc = 1. _d 0/(1. _d 0+0.5 _d 0*z_erfcc)

           erfcc = t_erfcc*exp(-z_erfcc*z_erfcc-1.26551223+t_erfcc*
     &       (1.00002368+t_erfcc*(0.37409196+t_erfcc*
     &       (.09678418+t_erfcc*(-.18628806+t_erfcc*(.27886807+
     &       t_erfcc*(-1.13520398+t_erfcc*(1.48851587+
     &       t_erfcc*(-0.82215223+t_erfcc*0.17087277)))))))))

          if (x_erfcc .lt. 0.0) then
            erfcc = 2.0 - erfcc
          endif

          fdvm_stat = (1. _d 0 - frac_migr) / 2. _d 0 * erfcc

c  Add the shapes, resulting in the 3-d DVM flux operator. If the
c  current layer is the bottom layer, or the layer beneath the
c  underlying layer is suboxic, all fluxes at and below the current
c  layer remain at the initialized value of zero. This will cause all
c  remaining DVM remineralization to occur in this layer.
          IF (k.LT.NR-1) THEN
            if (PTR_O2(i,j,k+2) .lt. (5. _d 0*oxic_min)) tmp = 1
          ENDIF
c!!          if (k .eq. grid_kmt(i,j)) exit
          dvm(i,j,k)  = fdvm_migr + fdvm_stat

         ENDIF

        enddo

c  Sum up the total organic flux to be transported by DVM

        do k = 1, nr
          fdvmn_vint  = fdvmn_vint  + N_dvm(i,j,k)  * drf(k)
          fdvmp_vint  = fdvmp_vint  + P_dvm(i,j,k)  * drf(k)
          fdvmfe_vint = fdvmfe_vint + Fe_dvm(i,j,k) * drf(k)
        enddo

c  Calculate the remineralization terms as the divergence of the flux

        N_remindvm(i,j,1)  = fdvmn_vint *  (1 - dvm(i,j,1)) /
     &     (epsln + drf(1))
        P_remindvm(i,j,1)  = fdvmp_vint *  (1 - dvm(i,j,1)) /
     &     (epsln + drf(1))
        Fe_remindvm(i,j,1) = fdvmfe_vint * (1 - dvm(i,j,1)) /
     &     (epsln + drf(1))

        do k = 2, nr
          N_remindvm(i,j,k)  = fdvmn_vint  *
     &       (dvm(i,j,k-1) - dvm(i,j,k)) / (epsln + drf(k))
          P_remindvm(i,j,k)  = fdvmp_vint  *
     &       (dvm(i,j,k-1) - dvm(i,j,k)) / (epsln + drf(k))
          Fe_remindvm(i,j,k) = fdvmfe_vint *
     &       (dvm(i,j,k-1) - dvm(i,j,k)) / (epsln + drf(k))
        enddo

       enddo
      enddo

#endif /* ALLOW_BLING */

      RETURN
      END
1	mmazloff	1.6	C $Header: /u/gcmpack/MITgcm/pkg/bling/bling_dvm.F,v 1.5 2016/09/12 20:00:28 mmazloff Exp $
2	mmazloff	1.1	C $Name: $
3
4			#include "BLING_OPTIONS.h"
5
6			CBOP
7			subroutine BLING_DVM(
8	jmc	1.2	I N_dvm,P_dvm,Fe_dvm,
9	mmazloff	1.1	I PTR_O2, mld,
10			O N_remindvm, P_remindvm, Fe_remindvm,
11			I bi, bj, imin, imax, jmin, jmax,
12			I myIter, myTime, myThid )
13
14			C =================================================================
15	mmazloff	1.5	C \| subroutine bling_dvm
16			C \| o Diel Vertical Migration
17	mmazloff	1.1	C =================================================================
18
19			implicit none
20	jmc	1.2
21	mmazloff	1.1	C === Global variables ===
22
23			#include "SIZE.h"
24			#include "DYNVARS.h"
25			#include "EEPARAMS.h"
26			#include "PARAMS.h"
27			#include "GRID.h"
28			#include "BLING_VARS.h"
29			#include "PTRACERS_SIZE.h"
30			#include "PTRACERS_PARAMS.h"
31			#ifdef ALLOW_AUTODIFF
32			# include "tamc.h"
33			#endif
34
35			C === Routine arguments ===
36			C bi,bj :: tile indices
37			C iMin,iMax :: computation domain: 1rst index range
38			C jMin,jMax :: computation domain: 2nd index range
39			C myTime :: current time
40			C myIter :: current timestep
41			C myThid :: thread Id. number
42			INTEGER bi, bj, imin, imax, jmin, jmax
43			_RL myTime
44			INTEGER myIter
45			INTEGER myThid
46			C === Input ===
47	mmazloff	1.5	C N_dvm :: vertical transport of nitrogen by DVM
48			C P_dvm :: vertical transport of phosphorus by DVM
49			C Fe_dvm :: vertical transport of iron by DVM
50			C PTR_O2 :: nitrate concentration
51			C mld :: mixed layer depth
52			_RL N_dvm (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
53			_RL P_dvm (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
54			_RL Fe_dvm(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
55			_RL PTR_O2(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
56			_RL mld (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
57
58	mmazloff	1.1	C === Output ===
59	mmazloff	1.5	C N_remindvm :: nitrogen remineralization due to diel vertical migration
60			C P_remindvm :: phosphorus remineralization due to diel vertical migration
61			C Fe_remindvm :: iron remineralization due to diel vertical migration
62			_RL N_remindvm (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
63			_RL P_remindvm (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
64			_RL Fe_remindvm (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
65	jmc	1.2
66	mmazloff	1.1	#ifdef ALLOW_BLING
67			C === Local variables ===
68			C i,j,k :: loop indices
69	jmc	1.2	INTEGER i,j,k
70	mmazloff	1.1	INTEGER tmp
71			_RL depth_l
72			_RL o2_upper
73	jmc	1.2	_RL o2_lower
74			_RL dz_upper
75			_RL dz_lower
76			_RL temp_upper
77			_RL temp_lower
78			_RL z_dvm_regr
79			_RL frac_migr
80			_RL fdvm_migr
81			_RL fdvm_stat
82			_RL fdvmn_vint
83			_RL fdvmp_vint
84			_RL fdvmfe_vint
85	mmazloff	1.1	_RL z_dvm
86			_RL dvm(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
87			_RL Fe_burial(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
88			_RL x_erfcc,z_erfcc,t_erfcc,erfcc
89			CEOP
90	jmc	1.2
91	mmazloff	1.1	c ---------------------------------------------------------------------
92			c Initialize output and diagnostics
93			DO k=1,Nr
94			DO j=jmin,jmax
95			DO i=imin,imax
96	mmazloff	1.5	N_remindvm(i,j,k) = 0. _d 0
97			P_remindvm(i,j,k) = 0. _d 0
98			Fe_remindvm(i,j,k) = 0. _d 0
99			dvm(i,j,k) = 0. _d 0
100	mmazloff	1.1	ENDDO
101	mmazloff	1.5	Fe_burial(i,j) = 0. _d 0
102	mmazloff	1.1	ENDDO
103			ENDDO
104	mmazloff	1.5
105	mmazloff	1.1	C ---------------------------------------------------------------------
106			c DIEL VERTICAL MIGRATOR EXPORT
107			c The effect of vertically-migrating animals on the export flux of organic
108			c matter from the ocean surface is treated similarly to the scheme of
109	jmc	1.2	c Bianchi et al., Nature Geoscience 2013.
110	mmazloff	1.1	c This involves calculating the stationary depth of vertical migrators, using
111	jmc	1.2	c an empirical multivariate regression, and ensuring that this remains
112	mmazloff	1.1	c above the bottom as well as any suboxic waters.
113	jmc	1.2	c The total DVM export flux is partitioned between a swimming migratory
114	mmazloff	1.1	c component and the stationary component, and these are summed.
115
116			C$TAF LOOP = parallel
117			DO j=jmin,jmax
118			C$TAF LOOP = parallel
119			DO i=imin,imax
120
121	mmazloff	1.5	c Initialize
122	jmc	1.2	o2_upper = 0.
123			o2_lower = 0.
124			dz_upper = 0.
125			dz_lower = 0.
126			temp_upper = 0.
127			temp_lower = 0.
128			z_dvm_regr = 0.
129	jmc	1.3	z_dvm = 0.
130	mmazloff	1.1	frac_migr = 0.
131	jmc	1.2	fdvm_migr = 0.
132			fdvm_stat = 0.
133			fdvmn_vint = 0.
134			fdvmp_vint = 0.
135	mmazloff	1.1	fdvmfe_vint = 0.
136
137			DO k=1,Nr
138
139			IF ( hFacC(i,j,k,bi,bj).gt.0. _d 0 ) THEN
140
141	mmazloff	1.5	c Calculate the depth of migration based on linear regression.
142	mmazloff	1.1
143			depth_l=-rF(k+1)
144	jmc	1.2
145	mmazloff	1.5	c Average temperature and oxygen over upper 35 m, and 140-515m.
146			c Also convert O2 to mmol m-3.
147	mmazloff	1.1
148			if ( abs(depth_l) .lt. 35.) then
149			dz_upper = dz_upper + drf(k)
150			temp_upper = temp_upper + theta(i,j,k,bi,bj)*drf(k)
151			o2_upper = o2_upper + PTR_O2(i,j,k) * drf(k)*1.0 _d 3
152			endif
153	jmc	1.2	if ( (abs(depth_l) .gt. 140.0 _d 0) .and.
154	mmazloff	1.1	& (abs(depth_l) .lt. 515. _d 0)) then
155			dz_lower = dz_lower + drf(k)
156			temp_lower = temp_lower + theta(i,j,k,bi,bj)*drf(k)
157			o2_lower = o2_lower + PTR_O2(i,j,k) * drf(k)*1.0 _d 3
158			endif
159
160			ENDIF
161			ENDDO
162	jmc	1.2
163	mmazloff	1.1	o2_upper = o2_upper / (epsln + dz_upper)
164			temp_upper = temp_upper / (epsln + dz_upper)
165			o2_lower = o2_lower / (epsln + dz_lower)
166			temp_lower = temp_lower / (epsln + dz_lower)
167
168	mmazloff	1.5	c Calculate the regression, using the constants given in Bianchi et al. (2013).
169			c The variable values are bounded to lie within reasonable ranges:
170			c O2 gradient : [-10,300] mmol/m3
171			c Log10 Chl : [-1.8,0.85] log10(mg/m3)
172			c mld : [0,500] m
173			c T gradient : [-3,20] C
174	mmazloff	1.1
175	mmazloff	1.5	c!! replacing hblt_depth(i,j) with mld...
176	mmazloff	1.1
177			#ifdef BLING_ADJOINT_SAFE
178	mmazloff	1.5
179	mmazloff	1.1	z_dvm = 300. _d 0
180
181			#else
182
183	jmc	1.2	z_dvm_regr = 398. _d 0
184	mmazloff	1.1	& - 0.56 _d 0*min(300. _d 0,max(-10. _d 0,(o2_upper - o2_lower)))
185			& - 115. _d 0*min(0.85 _d 0,max(-1.80 _d 0,
186	mmazloff	1.5	& log10(max(chl(i,j,1,bi,bj),chl_min))))
187	mmazloff	1.1	& + 0.36 _d 0*min(500. _d 0,max(epsln,mld(i,j)))
188			& - 2.40 _d 0*min(20. _d 0,max(-3. _d 0,(temp_upper-temp_lower)))
189
190	mmazloff	1.5	c Limit the depth of migration in polar winter.
191			c Use irr_mem since this is averaged over multiple days, dampening the
192			c diurnal cycle.
193			c Tapers Z_DVM to the minimum when surface irradince is below a given
194			c threshold (here 10 W/m2).
195	mmazloff	1.1
196	jmc	1.2	if ( irr_mem(i,j,1,bi,bj) .lt. 10. ) then
197			z_dvm_regr = 150. _d 0 + (z_dvm_regr - 150. _d 0) *
198	mmazloff	1.1	& irr_mem(i,j,1,bi,bj) / 10. _d 0
199	jmc	1.2	endif
200	mmazloff	1.1
201	mmazloff	1.5	c Check for suboxic water within the column. If found, set dvm
202			c stationary depth to 2 layers above it. This is not meant to
203			c represent a cessation of downward migration, but rather the
204			c requirement for aerobic DVM respiration to occur above the suboxic
205			c water, where O2 is available.
206	jmc	1.2
207	mmazloff	1.1	tmp = 0
208			DO k=1,Nr-2
209
210			IF ( (hFacC(i,j,k,bi,bj).gt.0. _d 0) .and. (tmp.eq.0)) THEN
211
212			z_dvm = -rf(k+1)
213			if (PTR_O2(i,j,k+2) .lt. (5. _d 0*oxic_min)) tmp = 1
214	jmc	1.2
215	mmazloff	1.1	ENDIF
216
217	jmc	1.2	enddo
218
219	mmazloff	1.5	c The stationary depth is constrained between 150 and 700, above any
220			c anoxic waters found, and above the bottom.
221	mmazloff	1.1
222			z_dvm = min(700. _d 0,max(150. _d 0,z_dvm_regr),z_dvm,-rf(k+1))
223			c!! bling%zbot(i,j,grid_kmt(i,j))) * grid_tmask(i,j,1)
224	mmazloff	1.5	c!! what is grid_kmt?
225	mmazloff	1.1
226			#endif
227
228	mmazloff	1.5	c Calculate the fraction of migratory respiration that occurs
229			c during upwards and downwards swimming. The remainder is
230			c respired near the stationary depth.
231			c Constants for swimming speed and resting time are hard-coded
232			c after Bianchi et al, Nature Geoscience 2013.
233	mmazloff	1.1
234			frac_migr = max( 0.0 _d 0, min( 1.0 _d 0, (2.0 _d 0 * z_dvm) /
235			& (epsln + 0.05 _d 0 * 0.5 _d 0 * 86400. _d 0)))
236
237	mmazloff	1.5	c Calculate the vertical profile shapes of DVM fluxes.
238			c These are given as the downward organic flux due to migratory
239			c DVM remineralization, defined at the bottom of each layer k.
240	mmazloff	1.1
241			tmp = 0
242			DO k=1,Nr
243
244			IF ( (hFacC(i,j,k,bi,bj).gt.0. _d 0) .and. (tmp.eq.0)) THEN
245
246	jmc	1.3	! First, calculate the part due to active migration above
247			! the stationary depth.
248	jmc	1.2	if (-rf(k+1) .lt. z_dvm) then
249			fdvm_migr = frac_migr / (epsln + z_dvm - (-rf(2))) *
250	mmazloff	1.1	& (z_dvm - (-rf(k+1)) )
251	jmc	1.2	else
252	mmazloff	1.1	fdvm_migr = 0.0
253	jmc	1.2	endif
254	mmazloff	1.1
255	mmazloff	1.5	c Then, calculate the part at the stationary depth.
256
257			c Approximation of the complementary error function
258			c From Numerical Recipes (F90, Ch. 6, p. 216)
259			c Returns the complementary error function erfc(x)
260			c with fractional error everywhere less than 1.2e-7
261	jmc	1.2	x_erfcc = (-rf(k) - z_dvm) /
262	mmazloff	1.1	& ( (epsln + 2. _d 0 * sigma_dvm2. _d 0)0.5)
263
264			z_erfcc = abs(x_erfcc)
265
266			t_erfcc = 1. _d 0/(1. _d 0+0.5 _d 0*z_erfcc)
267
268			erfcc = t_erfccexp(-z_erfccz_erfcc-1.26551223+t_erfcc*
269			& (1.00002368+t_erfcc(0.37409196+t_erfcc
270			& (.09678418+t_erfcc(-.18628806+t_erfcc(.27886807+
271			& t_erfcc(-1.13520398+t_erfcc(1.48851587+
272			& t_erfcc(-0.82215223+t_erfcc0.17087277)))))))))
273
274			if (x_erfcc .lt. 0.0) then
275			erfcc = 2.0 - erfcc
276			endif
277	jmc	1.2
278			fdvm_stat = (1. _d 0 - frac_migr) / 2. _d 0 * erfcc
279
280	mmazloff	1.5	c Add the shapes, resulting in the 3-d DVM flux operator. If the
281			c current layer is the bottom layer, or the layer beneath the
282			c underlying layer is suboxic, all fluxes at and below the current
283			c layer remain at the initialized value of zero. This will cause all
284			c remaining DVM remineralization to occur in this layer.
285	mmazloff	1.1	IF (k.LT.NR-1) THEN
286			if (PTR_O2(i,j,k+2) .lt. (5. _d 0*oxic_min)) tmp = 1
287			ENDIF
288			c!! if (k .eq. grid_kmt(i,j)) exit
289			dvm(i,j,k) = fdvm_migr + fdvm_stat
290	jmc	1.2
291	mmazloff	1.1	ENDIF
292
293	jmc	1.2	enddo
294
295	mmazloff	1.5	c Sum up the total organic flux to be transported by DVM
296	mmazloff	1.1
297	jmc	1.2	do k = 1, nr
298	mmazloff	1.1	fdvmn_vint = fdvmn_vint + N_dvm(i,j,k) * drf(k)
299			fdvmp_vint = fdvmp_vint + P_dvm(i,j,k) * drf(k)
300			fdvmfe_vint = fdvmfe_vint + Fe_dvm(i,j,k) * drf(k)
301	jmc	1.2	enddo
302	mmazloff	1.1
303	mmazloff	1.5	c Calculate the remineralization terms as the divergence of the flux
304	mmazloff	1.1
305			N_remindvm(i,j,1) = fdvmn_vint * (1 - dvm(i,j,1)) /
306			& (epsln + drf(1))
307			P_remindvm(i,j,1) = fdvmp_vint * (1 - dvm(i,j,1)) /
308			& (epsln + drf(1))
309	jmc	1.2	Fe_remindvm(i,j,1) = fdvmfe_vint * (1 - dvm(i,j,1)) /
310	mmazloff	1.1	& (epsln + drf(1))
311
312			do k = 2, nr
313	jmc	1.2	N_remindvm(i,j,k) = fdvmn_vint *
314	mmazloff	1.1	& (dvm(i,j,k-1) - dvm(i,j,k)) / (epsln + drf(k))
315	jmc	1.2	P_remindvm(i,j,k) = fdvmp_vint *
316	mmazloff	1.1	& (dvm(i,j,k-1) - dvm(i,j,k)) / (epsln + drf(k))
317			Fe_remindvm(i,j,k) = fdvmfe_vint *
318			& (dvm(i,j,k-1) - dvm(i,j,k)) / (epsln + drf(k))
319	jmc	1.2	enddo
320	mmazloff	1.1
321	jmc	1.2	enddo
322	mmazloff	1.1	enddo
323
324			#endif /* ALLOW_BLING */
325
326			RETURN
327			END