/[MITgcm]/MITgcm/pkg/bling/bling_dvm.F
ViewVC logotype

Diff of /MITgcm/pkg/bling/bling_dvm.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph | View Patch Patch

revision 1.1 by mmazloff, Thu May 19 20:29:26 2016 UTC revision 1.3 by jmc, Mon May 23 21:47:56 2016 UTC
# Line 5  C $Name$ Line 5  C $Name$
5    
6  CBOP  CBOP
7        subroutine BLING_DVM(        subroutine BLING_DVM(
8       I           N_dvm,P_dvm,Fe_dvm,       I           N_dvm,P_dvm,Fe_dvm,
9       I           PTR_O2, mld,       I           PTR_O2, mld,
10       O           N_remindvm, P_remindvm, Fe_remindvm,       O           N_remindvm, P_remindvm, Fe_remindvm,
11       I           bi, bj, imin, imax, jmin, jmax,       I           bi, bj, imin, imax, jmin, jmax,
# Line 14  CBOP Line 14  CBOP
14  C     =================================================================  C     =================================================================
15  C     | subroutine bling_prod  C     | subroutine bling_prod
16  C     | o Nutrient uptake and partitioning between organic pools.  C     | o Nutrient uptake and partitioning between organic pools.
17  C     | - Phytoplankton biomass-specific growth rate is calculated  C     | - Phytoplankton biomass-specific growth rate is calculated
18  C     |   as a function of light, nutrient limitation, and    C     |   as a function of light, nutrient limitation, and
19  C     |   temperature.  C     |   temperature.
20  C     | - Biomass growth xxx  C     | - Biomass growth xxx
21  C     | o Organic matter export, remineralization, and recycling.  C     | o Organic matter export, remineralization, and recycling.
22  C     | - Sinking particulate flux and diel migration contribute to  C     | - Sinking particulate flux and diel migration contribute to
23  C     |   export.  C     |   export.
24  C     | - Denitrification xxx  C     | - Denitrification xxx
25  C     =================================================================  C     =================================================================
26    
27        implicit none        implicit none
28          
29  C     === Global variables ===  C     === Global variables ===
30  C     P_sm          :: Small phytoplankton biomass  C     P_sm          :: Small phytoplankton biomass
31  C     P_lg          :: Large phytoplankton biomass  C     P_lg          :: Large phytoplankton biomass
# Line 77  C     G_xxx         :: Tendency of xxx Line 77  C     G_xxx         :: Tendency of xxx
77        _RL     G_DON     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     G_DON     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
78        _RL     G_DOP     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     G_DOP     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
79        _RL     G_CACO3   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     G_CACO3   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
80    
81  #ifdef ALLOW_BLING  #ifdef ALLOW_BLING
82  C     === Local variables ===  C     === Local variables ===
83  C     i,j,k         :: loop indices  C     i,j,k         :: loop indices
84  C     irr_eff       :: effective irradiance  C     irr_eff       :: effective irradiance
85  C     NO3_lim       :: nitrate limitation  C     NO3_lim       :: nitrate limitation
86  C     PO4_lim       :: phosphate limitation            C     PO4_lim       :: phosphate limitation
87  C     Fe_lim        :: iron limitation for phytoplankton      C     Fe_lim        :: iron limitation for phytoplankton
88  C     Fe_lim_diaz   :: iron limitation for diazotrophs  C     Fe_lim_diaz   :: iron limitation for diazotrophs
89  C     alpha_Fe      :: initial slope of the P-I curve  C     alpha_Fe      :: initial slope of the P-I curve
90  C     theta_Fe      :: Chl:C ratio  C     theta_Fe      :: Chl:C ratio
91  C     theta_Fe_max  :: Fe-replete maximum Chl:C ratio    C     theta_Fe_max  :: Fe-replete maximum Chl:C ratio
92  C     irrk          :: nut-limited efficiency of algal photosystems  C     irrk          :: nut-limited efficiency of algal photosystems
93  C     Pc_m          :: light-saturated max photosynthesis rate for phyt  C     Pc_m          :: light-saturated max photosynthesis rate for phyt
94  C     Pc_m_diaz     :: light-saturated max photosynthesis rate for diaz  C     Pc_m_diaz     :: light-saturated max photosynthesis rate for diaz
# Line 97  C     expkT         :: temperature funct Line 97  C     expkT         :: temperature funct
97  C     mu            :: net carbon-specific growth rate for phyt  C     mu            :: net carbon-specific growth rate for phyt
98  C     mu_diaz       :: net carbon-specific growth rate for diaz  C     mu_diaz       :: net carbon-specific growth rate for diaz
99  C     biomass_sm    :: nitrogen concentration in small phyto biomass  C     biomass_sm    :: nitrogen concentration in small phyto biomass
100  C     biomass_lg    :: nitrogen concentration in large phyto biomass  C     biomass_lg    :: nitrogen concentration in large phyto biomass
101  C     N_uptake      :: nitrogen uptake  C     N_uptake      :: nitrogen uptake
102  C     N_fix         :: nitrogen fixation  C     N_fix         :: nitrogen fixation
103  C     P_uptake      :: phosphorus uptake  C     P_uptake      :: phosphorus uptake
104  C     POC_flux      :: carbon export flux 3d field  C     POC_flux      :: carbon export flux 3d field
105  C     PtoN          :: variable ratio of phosphorus to nitrogen  C     PtoN          :: variable ratio of phosphorus to nitrogen
106  C     FetoN         :: variable ratio of iron to nitrogen  C     FetoN         :: variable ratio of iron to nitrogen
107  C     N_spm         :: particulate sinking of nitrogen  C     N_spm         :: particulate sinking of nitrogen
108  C     P_spm         :: particulate sinking of phosphorus  C     P_spm         :: particulate sinking of phosphorus
109  C     Fe_spm        :: particulate sinking of iron    C     Fe_spm        :: particulate sinking of iron
110  C     N_dvm         :: vertical transport of nitrogen by DVM  C     N_dvm         :: vertical transport of nitrogen by DVM
111  C     P_dvm         :: vertical transport of phosphorus by DVM  C     P_dvm         :: vertical transport of phosphorus by DVM
112  C     Fe_dvm        :: vertical transport of iron by DVM  C     Fe_dvm        :: vertical transport of iron by DVM
113  C     N_recycle     :: recycling of newly-produced organic nitrogen  C     N_recycle     :: recycling of newly-produced organic nitrogen
114  C     P_recycle     :: recycling of newly-produced organic phosphorus  C     P_recycle     :: recycling of newly-produced organic phosphorus
115  C     Fe_recycle    :: recycling of newly-produced organic iron  C     Fe_recycle    :: recycling of newly-produced organic iron
116  c xxx to be completed  c xxx to be completed
117        INTEGER i,j,k                INTEGER i,j,k
118        INTEGER tmp        INTEGER tmp
119        _RL irr_eff(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL irr_eff(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
120        _RL NO3_lim                          _RL NO3_lim
121        _RL PO4_lim        _RL PO4_lim
122        _RL Fe_lim        _RL Fe_lim
123        _RL Fe_lim_diaz                          _RL Fe_lim_diaz
124        _RL expkT        _RL expkT
125        _RL Pc_m        _RL Pc_m
126        _RL Pc_m_diaz        _RL Pc_m_diaz
127        _RL theta_Fe_max              _RL theta_Fe_max
128        _RL theta_Fe        _RL theta_Fe
129        _RL irrk          _RL irrk
130        _RL mu(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL mu(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
131        _RL mu_diaz(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL mu_diaz(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
132        _RL PtoN(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL PtoN(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
# Line 180  c xxx to be completed Line 180  c xxx to be completed
180        _RL Fe_reminp(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL Fe_reminp(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
181        _RL CacO3_diss(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL CacO3_diss(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
182        _RL o2_upper        _RL o2_upper
183        _RL o2_lower        _RL o2_lower
184        _RL dz_upper        _RL dz_upper
185        _RL dz_lower        _RL dz_lower
186        _RL temp_upper        _RL temp_upper
187        _RL temp_lower        _RL temp_lower
188        _RL z_dvm_regr        _RL z_dvm_regr
189        _RL frac_migr        _RL frac_migr
190        _RL fdvm_migr        _RL fdvm_migr
191        _RL fdvm_stat        _RL fdvm_stat
192        _RL fdvmn_vint        _RL fdvmn_vint
193        _RL fdvmp_vint        _RL fdvmp_vint
194        _RL fdvmfe_vint        _RL fdvmfe_vint
195        _RL z_dvm        _RL z_dvm
196        _RL N_remindvm(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL N_remindvm(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
197        _RL P_remindvm(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL P_remindvm(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
# Line 202  c xxx to be completed Line 202  c xxx to be completed
202        _RL x_erfcc,z_erfcc,t_erfcc,erfcc        _RL x_erfcc,z_erfcc,t_erfcc,erfcc
203  cxx order  cxx order
204  CEOP  CEOP
205      
206  c ---------------------------------------------------------------------  c ---------------------------------------------------------------------
207  c  Initialize output and diagnostics  c  Initialize output and diagnostics
208         DO k=1,Nr         DO k=1,Nr
# Line 261  c  Initialize output and diagnostics Line 261  c  Initialize output and diagnostics
261          ENDDO          ENDDO
262  cxx order  cxx order
263    
   
264  C  ---------------------------------------------------------------------  C  ---------------------------------------------------------------------
265  c DIEL VERTICAL MIGRATOR EXPORT  c DIEL VERTICAL MIGRATOR EXPORT
266  c The effect of vertically-migrating animals on the export flux of organic  c The effect of vertically-migrating animals on the export flux of organic
267  c matter from the ocean surface is treated similarly to the scheme of  c matter from the ocean surface is treated similarly to the scheme of
268  c Bianchi et al., Nature Geoscience 2013.  c Bianchi et al., Nature Geoscience 2013.
269  c This involves calculating the stationary depth of vertical migrators, using  c This involves calculating the stationary depth of vertical migrators, using
270  c an empirical multivariate regression, and ensuring that this remains  c an empirical multivariate regression, and ensuring that this remains
271  c above the bottom as well as any suboxic waters.  c above the bottom as well as any suboxic waters.
272  c The total DVM export flux is partitioned between a swimming migratory  c The total DVM export flux is partitioned between a swimming migratory
273  c component and the stationary component, and these are summed.  c component and the stationary component, and these are summed.
274    
   
   
275  C$TAF LOOP = parallel  C$TAF LOOP = parallel
276        DO j=jmin,jmax        DO j=jmin,jmax
277  C$TAF LOOP = parallel  C$TAF LOOP = parallel
278         DO i=imin,imax         DO i=imin,imax
279    
   
280  C        ! Initialize the working variables to zero  C        ! Initialize the working variables to zero
281          o2_upper = 0.          o2_upper = 0.
282          o2_lower = 0.          o2_lower = 0.
283          dz_upper = 0.          dz_upper = 0.
284          dz_lower = 0.          dz_lower = 0.
285          temp_upper = 0.          temp_upper = 0.
286          temp_lower = 0.          temp_lower = 0.
287          z_dvm_regr = 0.          z_dvm_regr = 0.
288            z_dvm     = 0.
289          frac_migr = 0.          frac_migr = 0.
290          fdvm_migr = 0.          fdvm_migr = 0.
291          fdvm_stat = 0.          fdvm_stat = 0.
292          fdvmn_vint = 0.          fdvmn_vint = 0.
293          fdvmp_vint = 0.          fdvmp_vint = 0.
294          fdvmfe_vint = 0.          fdvmfe_vint = 0.
295    
296          DO k=1,Nr          DO k=1,Nr
# Line 303  C        ! Initialize the working variab Line 300  C        ! Initialize the working variab
300        ! Calculate the depth of migration based on linear regression.        ! Calculate the depth of migration based on linear regression.
301    
302          depth_l=-rF(k+1)          depth_l=-rF(k+1)
303      
304          ! Average temperature and oxygen over upper 35 m, and 140-515m. Also convert O2 to mmol m-3.          ! Average temperature and oxygen over upper 35 m, and 140-515m.
305            ! Also convert O2 to mmol m-3.
306    
307            if ( abs(depth_l) .lt. 35.) then            if ( abs(depth_l) .lt. 35.) then
308              dz_upper = dz_upper + drf(k)              dz_upper = dz_upper + drf(k)
309              temp_upper = temp_upper + theta(i,j,k,bi,bj)*drf(k)              temp_upper = temp_upper + theta(i,j,k,bi,bj)*drf(k)
310              o2_upper = o2_upper + PTR_O2(i,j,k) * drf(k)*1.0 _d 3              o2_upper = o2_upper + PTR_O2(i,j,k) * drf(k)*1.0 _d 3
311            endif            endif
312            if ( (abs(depth_l) .gt. 140.0 _d 0) .and.            if ( (abs(depth_l) .gt. 140.0 _d 0) .and.
313       &          (abs(depth_l) .lt. 515. _d 0)) then       &          (abs(depth_l) .lt. 515. _d 0)) then
314              dz_lower = dz_lower + drf(k)              dz_lower = dz_lower + drf(k)
315              temp_lower = temp_lower + theta(i,j,k,bi,bj)*drf(k)              temp_lower = temp_lower + theta(i,j,k,bi,bj)*drf(k)
# Line 320  C        ! Initialize the working variab Line 318  C        ! Initialize the working variab
318    
319           ENDIF           ENDIF
320          ENDDO          ENDDO
321              
322          o2_upper = o2_upper / (epsln + dz_upper)          o2_upper = o2_upper / (epsln + dz_upper)
323          temp_upper = temp_upper / (epsln + dz_upper)          temp_upper = temp_upper / (epsln + dz_upper)
324          o2_lower = o2_lower / (epsln + dz_lower)          o2_lower = o2_lower / (epsln + dz_lower)
325          temp_lower = temp_lower / (epsln + dz_lower)          temp_lower = temp_lower / (epsln + dz_lower)
326    
327          ! Calculate the regression, using the constants given in Bianchi et al. (2013).          ! Calculate the regression, using the constants given
328          ! The variable values are bounded to lie within reasonable ranges:          ! in Bianchi et al. (2013).
329          ! O2 gradient   : [-10,300] mmol/m3          ! The variable values are bounded to lie within reasonable
330          ! Log10 Chl     : [-1.8,0.85] log10(mg/m3)          ! ranges: O2 gradient   : [-10,300] mmol/m3
331          ! mld           : [0,500] m          !         Log10 Chl     : [-1.8,0.85] log10(mg/m3)
332          ! T gradient    : [-3,20] C          !         mld           : [0,500] m
333            !         T gradient    : [-3,20] C
334    
335  C!! I'm replacing hblt_depth(i,j) with mld... not sure what hblt_depth is  C!! I am replacing hblt_depth(i,j) with mld... not sure what hblt_depth is
336    
337  #ifdef BLING_ADJOINT_SAFE  #ifdef BLING_ADJOINT_SAFE
338          z_dvm = 300. _d 0          z_dvm = 300. _d 0
339    
340  #else  #else
341    
342          z_dvm_regr = 398. _d 0          z_dvm_regr = 398. _d 0
343       &   - 0.56 _d 0*min(300. _d 0,max(-10. _d 0,(o2_upper - o2_lower)))       &   - 0.56 _d 0*min(300. _d 0,max(-10. _d 0,(o2_upper - o2_lower)))
344       &   - 115. _d 0*min(0.85 _d 0,max(-1.80 _d 0,       &   - 115. _d 0*min(0.85 _d 0,max(-1.80 _d 0,
345       &                            log10(max(chl(i,j,1,bi,bj),chl_min))))       &                            log10(max(chl(i,j,1,bi,bj),chl_min))))
# Line 351  c        ! Limit the depth of migration Line 350  c        ! Limit the depth of migration
350  c        ! Use irr_mem since this is averaged over multiple days, dampening the diurnal cycle.  c        ! Use irr_mem since this is averaged over multiple days, dampening the diurnal cycle.
351  c        ! Tapers Z_DVM to the minimum when surface irradince is below a given threshold (here 10 W/m2).  c        ! Tapers Z_DVM to the minimum when surface irradince is below a given threshold (here 10 W/m2).
352    
353          if ( irr_mem(i,j,1,bi,bj) .lt. 10. ) then                                if ( irr_mem(i,j,1,bi,bj) .lt. 10. ) then
354            z_dvm_regr = 150. _d 0 + (z_dvm_regr - 150. _d 0) *            z_dvm_regr = 150. _d 0 + (z_dvm_regr - 150. _d 0) *
355       &       irr_mem(i,j,1,bi,bj) / 10. _d 0       &       irr_mem(i,j,1,bi,bj) / 10. _d 0
356          endif                                                                    endif
357    
358  C Check for suboxic water within the column. If found, set dvm    C Check for suboxic water within the column. If found, set dvm
359  C stationary depth to 2 layers above it. This is not meant to  C stationary depth to 2 layers above it. This is not meant to
360  C represent a cessation of downward migration, but rather the  C represent a cessation of downward migration, but rather the
361  C requirement for aerobic DVM respiration to occur above the suboxic  C requirement for aerobic DVM respiration to occur above the suboxic
362  C water, where O2 is available.  C water, where O2 is available.
363    
364          tmp = 0          tmp = 0
365          DO k=1,Nr-2          DO k=1,Nr-2
366    
# Line 369  C water, where O2 is available. Line 368  C water, where O2 is available.
368    
369            z_dvm = -rf(k+1)            z_dvm = -rf(k+1)
370            if (PTR_O2(i,j,k+2) .lt. (5. _d 0*oxic_min)) tmp = 1            if (PTR_O2(i,j,k+2) .lt. (5. _d 0*oxic_min)) tmp = 1
371              
372           ENDIF           ENDIF
             
         enddo    
373    
374  C The stationary depth is constrained between 150 and 700, above any            enddo
375    
376    C The stationary depth is constrained between 150 and 700, above any
377  C anoxic waters found, and above the bottom.  C anoxic waters found, and above the bottom.
378    
379          z_dvm = min(700. _d 0,max(150. _d 0,z_dvm_regr),z_dvm,-rf(k+1))          z_dvm = min(700. _d 0,max(150. _d 0,z_dvm_regr),z_dvm,-rf(k+1))
# Line 383  c!! what is grid_mkt?? Line 382  c!! what is grid_mkt??
382    
383  #endif  #endif
384    
385          ! Calculate the fraction of migratory respiration that occurs during upwards          ! Calculate the fraction of migratory respiration that occurs
386          ! and downwards swimming. The remainder is respired near the stationary depth.          ! during upwards and downwards swimming. The remainder is
387          ! Constants for swimming speed and resting time are hard-coded after Bianchi          ! respired near the stationary depth.
388          ! et al, Nature Geoscience 2013.          ! Constants for swimming speed and resting time are hard-coded
389            ! after Bianchi et al, Nature Geoscience 2013.
390    
391          frac_migr = max( 0.0 _d 0, min( 1.0 _d 0, (2.0 _d 0 * z_dvm) /          frac_migr = max( 0.0 _d 0, min( 1.0 _d 0, (2.0 _d 0 * z_dvm) /
392       &        (epsln + 0.05 _d 0 * 0.5 _d 0 * 86400. _d 0)))       &        (epsln + 0.05 _d 0 * 0.5 _d 0 * 86400. _d 0)))
393    
394          ! Calculate the vertical profile shapes of DVM fluxes. These are given as          ! Calculate the vertical profile shapes of DVM fluxes.
395          ! the downward organic flux due to migratory DVM remineralization, defined at            ! These are given as the downward organic flux due to migratory
396          ! the bottom of each layer k.          ! DVM remineralization, defined at the bottom of each layer k.
   
397    
398          tmp = 0          tmp = 0
399          DO k=1,Nr          DO k=1,Nr
400    
401           IF ( (hFacC(i,j,k,bi,bj).gt.0. _d 0) .and. (tmp.eq.0)) THEN           IF ( (hFacC(i,j,k,bi,bj).gt.0. _d 0) .and. (tmp.eq.0)) THEN
402    
403            ! First, calculate the part due to active migration above the stationary depth.            ! First, calculate the part due to active migration above
404            if (-rf(k+1) .lt. z_dvm) then              ! the stationary depth.
405              fdvm_migr = frac_migr / (epsln + z_dvm - (-rf(2))) *                  if (-rf(k+1) .lt. z_dvm) then
406                fdvm_migr = frac_migr / (epsln + z_dvm - (-rf(2))) *
407       &            (z_dvm - (-rf(k+1)) )       &            (z_dvm - (-rf(k+1)) )
408            else              else
409              fdvm_migr  = 0.0              fdvm_migr  = 0.0
410            endif              endif
411    
412            ! Then, calculate the part at the stationary depth.            ! Then, calculate the part at the stationary depth.
413  c          fdvm_stat = (1. - frac_migr) / 2. * erfcc((-rf(k) - z_dvm) /  c          fdvm_stat = (1. - frac_migr) / 2. * erfcc((-rf(k) - z_dvm) /
414  c     &        ( (epsln + 2. * sigma_dvm**2.)**0.5))  c     &        ( (epsln + 2. * sigma_dvm**2.)**0.5))
415    
   
416  c  ! Approximation of the complementary error function  c  ! Approximation of the complementary error function
417  c  ! From Numerical Recipes (F90, Ch. 6, p. 216)  c  ! From Numerical Recipes (F90, Ch. 6, p. 216)
418  c  ! Returns the complementary error function erfc(x)  c  ! Returns the complementary error function erfc(x)
419  c    with fractional error everywhere less than 1.2e-7      c    with fractional error everywhere less than 1.2e-7
420             x_erfcc = (-rf(k) - z_dvm) /             x_erfcc = (-rf(k) - z_dvm) /
421       &        ( (epsln + 2. _d 0 * sigma_dvm**2. _d 0)**0.5)       &        ( (epsln + 2. _d 0 * sigma_dvm**2. _d 0)**0.5)
422    
423             z_erfcc = abs(x_erfcc)             z_erfcc = abs(x_erfcc)
# Line 434  c    with fractional error everywhere le Line 433  c    with fractional error everywhere le
433            if (x_erfcc .lt. 0.0) then            if (x_erfcc .lt. 0.0) then
434              erfcc = 2.0 - erfcc              erfcc = 2.0 - erfcc
435            endif            endif
436        
437                  fdvm_stat = (1. _d 0 - frac_migr) / 2. _d 0 * erfcc
438            fdvm_stat = (1. _d 0 - frac_migr) / 2. _d 0 * erfcc  
439          C Add the shapes, resulting in the 3-d DVM flux operator. If the
440  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
441  C current layer is the bottom layer, or the layer beneath the  C underlying layer is suboxic, all fluxes at and below the current
442  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 layer remain at the initialized value of zero. This will cause all  
443  C remaining DVM remineralization to occur in this layer.  C remaining DVM remineralization to occur in this layer.
444            IF (k.LT.NR-1) THEN            IF (k.LT.NR-1) THEN
445              if (PTR_O2(i,j,k+2) .lt. (5. _d 0*oxic_min)) tmp = 1              if (PTR_O2(i,j,k+2) .lt. (5. _d 0*oxic_min)) tmp = 1
446            ENDIF            ENDIF
447  c!!          if (k .eq. grid_kmt(i,j)) exit  c!!          if (k .eq. grid_kmt(i,j)) exit
448            dvm(i,j,k)  = fdvm_migr + fdvm_stat            dvm(i,j,k)  = fdvm_migr + fdvm_stat
449              
450           ENDIF           ENDIF
             
         enddo    
451    
452  c Sum up the total organic flux to be transported by DVM          enddo
453    
454    c Sum up the total organic flux to be transported by DVM
455    
456          do k = 1, nr          do k = 1, nr
457            fdvmn_vint  = fdvmn_vint  + N_dvm(i,j,k)  * drf(k)            fdvmn_vint  = fdvmn_vint  + N_dvm(i,j,k)  * drf(k)
458            fdvmp_vint  = fdvmp_vint  + P_dvm(i,j,k)  * drf(k)            fdvmp_vint  = fdvmp_vint  + P_dvm(i,j,k)  * drf(k)
459            fdvmfe_vint = fdvmfe_vint + Fe_dvm(i,j,k) * drf(k)            fdvmfe_vint = fdvmfe_vint + Fe_dvm(i,j,k) * drf(k)
460          enddo            enddo
461    
462  c Calculate the remineralization terms as the divergence of the flux  c Calculate the remineralization terms as the divergence of the flux
463    
# Line 467  c Calculate the remineralization terms a Line 465  c Calculate the remineralization terms a
465       &     (epsln + drf(1))       &     (epsln + drf(1))
466          P_remindvm(i,j,1)  = fdvmp_vint *  (1 - dvm(i,j,1)) /          P_remindvm(i,j,1)  = fdvmp_vint *  (1 - dvm(i,j,1)) /
467       &     (epsln + drf(1))       &     (epsln + drf(1))
468          Fe_remindvm(i,j,1) = fdvmfe_vint * (1 - dvm(i,j,1)) /          Fe_remindvm(i,j,1) = fdvmfe_vint * (1 - dvm(i,j,1)) /
469       &     (epsln + drf(1))       &     (epsln + drf(1))
470    
471          do k = 2, nr          do k = 2, nr
472            N_remindvm(i,j,k)  = fdvmn_vint  *            N_remindvm(i,j,k)  = fdvmn_vint  *
473       &       (dvm(i,j,k-1) - dvm(i,j,k)) / (epsln + drf(k))       &       (dvm(i,j,k-1) - dvm(i,j,k)) / (epsln + drf(k))
474            P_remindvm(i,j,k)  = fdvmp_vint  *            P_remindvm(i,j,k)  = fdvmp_vint  *
475       &       (dvm(i,j,k-1) - dvm(i,j,k)) / (epsln + drf(k))       &       (dvm(i,j,k-1) - dvm(i,j,k)) / (epsln + drf(k))
476            Fe_remindvm(i,j,k) = fdvmfe_vint *            Fe_remindvm(i,j,k) = fdvmfe_vint *
477       &       (dvm(i,j,k-1) - dvm(i,j,k)) / (epsln + drf(k))       &       (dvm(i,j,k-1) - dvm(i,j,k)) / (epsln + drf(k))
478          enddo            enddo
479    
480         enddo         enddo
481        enddo        enddo
482    
483  c ---------------------------------------------------------------------  c ---------------------------------------------------------------------
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.3
  ViewVC Help
Powered by ViewVC 1.1.22