#include "options.inc"

c---------------------------------------------------------------
c
c     A simple pelagic ecosystem model  (NPZD)
c     taken from Oschlies and Garcon, 1999, 
c                 Glob. Biogeochem. Cycles, 13, 135-160
c     acting as sources for dissolved inorganic carbon and oxygen
c
c     linked in the code in passive_tracer.F and in driver.F
c
c---------------------------------------------------------------



c local cpp options
c#define test_minmax
#define prevent_tracer_from_zero
c
c put any detritus at bottom to zero (into sediment)
c to prevent spurious accumulation at shelfs
c
c#define detritus_sedimentation
c another, more realistic version
c#define detritus_slow_sedimentation
c
c account for denitrification in case of low oxygen concentration
c
c#define denitrification
c
c
c delta CO_2 concentration depends also on sea level pressure
c account for that
c#define account_for_slp_in_co2_flux
c
c atmospheric pCO2 depends on latitude and season
c account for that
c#define account_for_seasonality_in_pco2
c
c
c default is wanninkhof 92 for piston velocity
c#define use_wanninkhof99
c
#define use_fast_co2_subr
c
      module npzd_module
      implicit none

      real :: alphabio = 0.025 /86400. ! add docu for biological parameters
      real :: abio     = 0.6   /86400.
      real :: bbio     = 1.066
      real :: cbio     = 1.0
      real :: dbio     = 2500.0
      real :: gbio     = 2.0  /86400.
      real :: epsbio   = 1.0  /86400.
      real :: a_npz    = 0.75
      real :: d_npz    = 0.03 /86400.
      real :: parbio   = 0.43
      real :: rk1bio   = 0.5
      real :: phiphy   = 0.03 /86400.
      real :: phizoo   = 0.20 /86400.
      real :: remina   = 0.05 /86400.
      real :: remina_slow   = 1./365. /86400.
      real :: w_detr   = 5.0  /86400.* 100.0
      real :: dkcbio   = 0.03

      real :: s_doc = 1./50.0/86400.0 ! production rate of DOC
      real :: r_doc = 1./800.0 /86400.0 ! reminarilisation rate for DOC

      integer :: ts_npzd ! number of time steps within NPZD model
      integer :: kmeuph  ! deepest level of  euphotic zone
c
c   Atmpospheric carbon dioxide partial pressure
c   in yearly intervals, from 1765 extending to the year 2200. Data are
c   from the spline fit to ico-core measurements of Friedli
c   et al., and the Mauna Loa observations by Keeling et al.,
c   combined with the IPCC scenario S650 (This scenario ends
c   at a steady concentration of 650 ppm in the year 2200).
c   PCO2(I) is the partial pressure in the middle of the year I
c   (i.e. year I + 0.5).

      real :: pco2_air_fixed = 278.0 ! fixed value instead of time series
      real :: pco2_air               ! actual value for both time series and fixed value
      integer :: pco2series_index   ! pointer for the interpolation method
      integer,parameter :: pco2series_dim=2200-1765+1
#ifdef notdef
      real :: pco2series(pco2series_dim) = (/
     &    277.965, 277.995, 278.027, 278.063, 278.101, 278.142, 278.186, ! - 1771
     &    278.233, 278.284, 278.337, 278.393, 278.452, 278.514, 278.579, ! - 1778
     &    278.648, 278.719, 278.794, 278.872, 278.953, 279.038, 279.125, ! - 1785
     &    279.216, 279.311, 279.409, 279.510, 279.614, 279.722, 279.833, ! - 1892
     &    279.948, 280.066, 280.188, 280.313, 280.441, 280.573, 280.708, ! - 1799
     &    280.845, 280.984, 281.124, 281.266, 281.408, 281.551, 281.693, ! - 1806
     &    281.835, 281.976, 282.116, 282.254, 282.390, 282.523, 282.653, ! - 1764+7*7
     &    282.779, 282.902, 283.021, 283.134, 283.244, 283.349, 283.452, ! - 1764+7*8
     &    283.551, 283.648, 283.743, 283.836, 283.929, 284.020, 284.112, ! - 1764+7*9
     &    284.203, 284.296, 284.390, 284.486, 284.584, 284.686, 284.791, ! - 1764+7*10 = 1834
     &    284.900, 285.013, 285.132, 285.256, 285.387, 285.523, 285.665, ! - 1764+7*11
     &    285.809, 285.953, 286.097, 286.239, 286.379, 286.516, 286.650, ! - 1764+7*12
     &    286.781, 286.909, 287.033, 287.154, 287.272, 287.387, 287.497, ! - 1764+7*13
     &    287.605, 287.711, 287.816, 287.920, 288.025, 288.131, 288.239, ! - 1764+7*14 = 1862
     &    288.350, 288.465, 288.584, 288.708, 288.838, 288.974, 289.115, ! - 1764+7*15
     &    289.262, 289.416, 289.576, 289.744, 289.919, 290.101, 290.292, ! - 1764+7*16
     &    290.490, 290.695, 290.908, 291.128, 291.354, 291.586, 291.824, ! - 1764+7*17
     &    292.066, 292.312, 292.562, 292.815, 293.071, 293.328, 293.586, ! - 1764+7*18
     &    293.843, 294.098, 294.350, 294.599, 294.842, 295.082, 295.320, ! - 1764+7*19
     &    295.558, 295.796, 296.038, 296.283, 296.535, 296.793, 297.061, ! - 1764+7*20 = 1904
     &    297.336, 297.620, 297.909, 298.204, 298.504, 298.806, 299.111, ! - 1764+7*21
     &    299.419, 299.728, 300.040, 300.352, 300.666, 300.980, 301.294, ! - 1764+7*22 = 1918
     &    301.609, 301.923, 302.237, 302.551, 302.863, 303.172, 303.478, ! - 1764+7*23
     &    303.780, 304.076, 304.367, 304.651, 304.931, 305.207, 305.478, ! - 1764+7*24
     &    305.747, 306.013, 306.280, 306.546, 306.815, 307.087, 307.365, ! - 1764+7*25 = 1939
     &    307.649, 307.942, 308.244, 308.559, 308.886, 309.227, 309.582, ! - 1764+7*26 = 1946
     &    309.954, 310.341, 310.746, 311.170, 311.612, 312.075, 312.559, ! - 1764+7*27 = 1953
     &    313.065, 313.596, 314.151, 314.733, 315.343, 315.981, 316.644, ! - 1764+7*28 = 1960
     &    317.327, 318.024, 318.738, 319.483, 320.275, 321.125, 322.038, ! - 1764+7*29 = 1967
     &    323.013, 324.052, 325.149, 326.294, 327.481, 328.697, 329.930, ! - 1764+7*30 = 1974
     &    331.189, 332.492, 333.847, 335.250, 336.687, 338.148, 339.625, ! - 1764+7*31 = 1981
     &    341.122, 342.646, 344.202, 345.792, 347.419, 349.084, 350.770, ! - 1764+7*32 = 1988
     &    352.470, 354.170, 355.874, 357.587, 359.309, 361.040, 362.779, ! - 1764+7*33 = 1995
     &    364.526, 366.282, 368.046, 369.818, 371.598, 373.386, 375.183, ! - 1764+7*34 = 2002
     &    376.987, 378.799, 380.618, 382.445, 384.280, 386.122, 387.971, ! - 1764+7*35
     &    389.827, 391.691, 393.561, 395.438, 397.321, 399.211, 401.108, ! - 1764+7*36
     &    403.010, 404.919, 406.833, 408.754, 410.680, 412.612, 414.549, ! - 1764+7*37
     &    416.491, 418.438, 420.390, 422.347, 424.308, 426.274, 428.244, ! - 1764+7*38
     &    430.217, 432.195, 434.176, 436.161, 438.149, 440.140, 442.134, ! - 1764+7*39
     &    444.131, 446.130, 448.132, 450.135, 452.141, 454.148, 456.156, ! - 1764+7*40 = 2044
     &    458.166, 460.177, 462.188, 464.200, 466.213, 468.225, 470.238, ! - 1764+7*41
     &    472.250, 474.261, 476.272, 478.281, 480.289, 482.296, 484.301, ! - 1764+7*42
     &    486.304, 488.304, 490.302, 492.298, 494.290, 496.279, 498.264, ! - 1764+7*43
     &    500.246, 502.224, 504.197, 506.166, 508.130, 510.089, 512.043, ! - 1764+7*44
     &    513.991, 515.933, 517.869, 519.799, 521.722, 523.638, 525.548,
     &    527.450, 529.344, 531.230, 533.108, 534.978, 536.839, 538.692,
     &    540.535, 542.369, 544.193, 546.007, 547.811, 549.605, 551.388,
     &    553.160, 554.921, 556.671, 558.409, 560.136, 561.850, 563.552,
     &    565.242, 566.919, 568.583, 570.234, 571.872, 573.496, 575.106,
     &    576.703, 578.285, 579.853, 581.407, 582.945, 584.469, 585.978,
     &    587.471, 588.949, 590.411, 591.858, 593.289, 594.703, 596.101,
     &    597.483, 598.849, 600.197, 601.529, 602.844, 604.142, 605.422,
     &    606.686, 607.932, 609.160, 610.371, 611.564, 612.739, 613.896,
     &    615.035, 616.156, 617.259, 618.344, 619.411, 620.459, 621.488,
     &    622.499, 623.492, 624.466, 625.422, 626.358, 627.276, 628.176,
     &    629.057, 629.919, 630.762, 631.587, 632.392, 633.179, 633.948,
     &    634.697, 635.428, 636.140, 636.834, 637.509, 638.165, 638.803,
     &    639.422, 640.023, 640.605, 641.169, 641.714, 642.242, 642.751,
     &    643.242, 643.714, 644.169, 644.606, 645.025, 645.427, 645.811,
     &    646.177, 646.525, 646.857, 647.170, 647.467, 647.747, 648.010,
     &    648.255, 648.485, 648.697, 648.893, 649.072, 649.236, 649.383,
     &    649.514, 649.629, 649.728, 649.812, 649.880, 649.933, 649.970,
     &    649.993, 650.000 /)
#else
      real :: pco2series(pco2series_dim) 
#endif
c
c   local array for pCO2 air-sea difference (all calulated in npzd_sflx)
c
      real, allocatable :: diff_co2(:,:) ! in mol/m^3
      real, allocatable :: pco2_water(:,:) ! partial pressure of surface CO2  (ppmv)
      real, allocatable :: o2_sat(:,:)  ! O2 saturation concentration at sea surface
c
c    frequency (in time steps) of updates for diff_co2/pco2_water
c
      integer  :: pco2_update_int = 10
c
c   Redfield ratio of C/N
c
      real :: redfield_CN = 6.6
      real :: redfield_CN2 = 15
c
c   Redfield ratio of O/N
chd   oxygen concentrations are in ml/l; donstant
chd   redfieldratio of N/O2=1/10 (ocmip) is assumed;
chd   it is: mlO2/l = 10*22.3916/1000 mmolN/m^3 
c
      real :: redfield_ON = 0.223916

#ifdef denitrification
c
c   treshold of oxygen for denitrification
c
      real :: o2_threshold = 0.08
#endif
c
c    diagnostic array
c
      real, allocatable  :: dnut_diag(:,:,:)

      end module npzd_module


      subroutine npzd_init
c---------------------------------------------------------------
c   Initialization of the NPZD module
c---------------------------------------------------------------
      use spflame_module
      use passive_tracer_module
      use npzd_module
      implicit none

      type( time_type) :: start
      type( time_type), allocatable :: interval(:)
      integer :: indp,k,n

      if (my_pe == 0) then
       print*,''
       print*,' Initialization of NPZD model'
       print*,''
      endif
c
c   display mapping of tracer number to biological variable
c
      if (my_pe==0) then
       print*,' Tracer # ',n_npzd(1),' will be nutrient'
       print*,' Tracer # ',n_npzd(2),' will be phytoplancton'
       print*,' Tracer # ',n_npzd(3),' will be zooplancton'
       print*,' Tracer # ',n_npzd(4),' will be detritus'
       print*,'  (',n_nutrient,n_phyto,n_zoo,n_detritus,')'
       if (enable_co2) then
        print*,' Tracer # ',n_co2,' will be dissolved inorganic carbon'
       endif
       if (enable_o2) then
        print*,' Tracer # ',n_o2,' will be total dissolved oxygen'
       endif
       if (enable_doc) then
        print*,' we have DOC as part of the NPZD model'
        print*,' Tracer # ',n_doc,' will be dissolved organic carbon'
       endif
c    print out values of parameters for NPZD
       print*,' the parameter for the ecosystem model:'
       print*,' alphabio = ',alphabio*86400,' 1/days' 
       print*,' abio     = ', abio* 86400., ' 1/days'
       print*,' bbio     = ', bbio 
       print*,' cbio     = ', cbio 
       print*,' dbio     = ', dbio
       print*,' gbio     = ', gbio*86400.,' 1/days'
       print*,' epsbio   = ',epsbio*86400.,' 1/days'
       print*,' a_npz    = ',a_npz
       print*,' d_npz    = ',d_npz*86400., ' 1/days'
       print*,' parbio   = ',parbio
       print*,' rk1bio   = ',rk1bio
       print*,' dkcbio   = ',dkcbio
       print*,' phiphy   = ',1./phiphy/86400.,' days'
       print*,' phizoo   = ',1./phizoo/86400.,' days'
       print*,' remina   = ',1./remina/86400.,' days'
       print*,' w_detr   = ',w_detr *86400./100.0 ,' m/day'
      endif
c
c   time step for NPZD model
c
      ts_npzd=max(1,nint(time_step*2. / 900.))
      if (my_pe==0) then
       print*,' NPZD model will do ',ts_npzd,' time steps'
       print*,' time step in NPZD model will be ',
     &          time_step/ts_npzd,'sec.'
      endif
c
c   Depth of euphotic zone
c
      kmeuph = indp(200.0e2,zt,km)
      if (my_pe==0) 
     &    print*,' depth of euphotic zone ',zt(kmeuph)/100.0,'m'

      if (enable_co2 .and. enable_pco2_series) then
c
c    interpolation setup for atmospheric partial pressure of CO2
c
       if (my_pe==0) then
         print*,' using realistic time series of atmospheric CO2 '
         print*, 'for suface fluxes of dissolved inorganic carbon'
         print*, 'time series starts in (model) year ',ref_year_co2
       endif

       start = set_date(ref_year_co2,1, 1, 0 , 0 ,0)
       allocate(interval(pco2series_dim))

       if (get_calendar_type() == THIRTY_DAY_MONTHS ) then
        do k=1,pco2series_dim
         call set_time(interval(k), 0, 360 )
        enddo
       elseif (get_calendar_type() == NO_LEAP ) then
        do k=1,pco2series_dim
         call set_time(interval(k), 0, 365 )
        enddo
       else 
        print*,' wrong calendar type for CO2'
        call halt_stop(' in npzd_init')
       endif

       call init_forcing_interp(pco2series_dim,.false.,
     &                          start,interval,pco2series_index)
       deallocate(interval)

       if (my_pe==0) then
        print*,' reading time series of atmospheric CO2 '
        print*,' from file pco2_series.dta'
       endif

       call getunit(n,'pco2_series.dta','fsr')
       do k=1,pco2series_dim
        read(n, * ) pco2series(k)
       enddo
       close(n)

      endif ! enable_co2 and enable pco2_series

      if (enable_co2 .and. .not. enable_pco2_series .and. my_pe==0) then
         print*,' using a fixed value for atmospheric CO2 '
         print*, ' pco2_air_fixed =  ',pco2_air_fixed
      endif

      if (enable_co2) then
       allocate( diff_co2( is_pe:ie_pe, js_pe:je_pe ) )
       allocate( pco2_water( is_pe:ie_pe, js_pe:je_pe ) )
       diff_co2=0.0; pco2_water=0.0
       if (my_pe==0) then
         print*,' updating air-sea difference of CO2 concentration',
     &            ' every ',pco2_update_int,' model time step'
         print*,' Redfield ratio of C/N = ',redfield_CN
       endif
      endif

      if (enable_o2) then
       allocate( o2_sat( is_pe:ie_pe, js_pe:je_pe ) )
       o2_sat = 0.0 ! diagnostic array
       if (my_pe==0) print*,' Redfield ratio of O/N = ',redfield_ON
      endif

      if (enable_doc) then
        if (my_pe == 0) then
          print*,' Okay, DOC is modeled as well ...'
          print*,' decay rate of DOC is ',1./r_doc/86400.,' days'
          print*,' production rate of DOC ',1./s_doc/86400.,'days'
          print*,' CN of DOC is ',redfield_CN2
        endif
      endif

      if (enable_diag_npzd) then
       if (my_pe==0) then
        print*,' diagnosing monthly means of npzd model'
       endif
       call init_npzd_averages()
       allocate( dnut_diag(is_pe:ie_pe,km,js_pe:je_pe) )
       dnut_diag = 0.
      endif

      if (my_pe==0) then
#ifdef detritus_sedimentation
        print*,' put any detritus at bottom to zero (into sediment)'
        print*,' to prevent spurious accumulation at shelfs'
#endif
#ifdef detritus_slow_sedimentation
        print*,' put detritus at bottom to zero (into sediment)'
        print*,' with a e-folding time scale of ',
     &                  1./remina_slow/86400.,' days'
#endif
#ifdef denitrification
        print*,' account for denitrification in case of '
        print*,' low oxygen concentration of ',o2_threshold
#endif
#ifdef prevent_tracer_from_zero
        print*,' put negative values of biological variables to zero'
#endif
        print*,' account for the dilutive effect of salt flux'
        print*,' account for bug in CO2 fluxes'
#ifdef account_for_slp_in_co2_flux
        print*,' account for slp calculating CO2 fluxes'
#endif
#ifdef account_for_seasonality_in_pco2
        print*,' account for seasonality in atmos. pCO2'
#endif
#ifdef use_wanninkhof99
        print*,' using Wanninkhof+McGillis 99 for piston velocity'
#else
        print*,' using Wanninkhof 92 for piston velocity'
#endif
#ifdef use_fast_co2_subr
       print*,' using fast routine for surface DIC-Eq. calculation'
#endif

      endif
      if (my_pe == 0) print*,  'done'

      end subroutine npzd_init


      subroutine npzd_sflx
c---------------------------------------------------------------
c   apply surface fluxes for biological variables
c---------------------------------------------------------------
      use spflame_module
      use passive_tracer_module
      use npzd_module
      implicit none
      integer :: n,i,j,is,ie,js,je,p1,p2
      real, dimension(is_pe:ie_pe) :: kw,sc,slp
      real :: f1,f2,pco2_atmos
      logical, save :: init = .true.

      is=max(is_pe,2); ie=min(ie_pe,imt-1)
      js=max(2,js_pe); je=min(je_pe,jmt-1)

      do n=1,nt_npzd
c    no surface flux for biological tracers
       do j=js,je
         do i=is,ie
           stf(i,j,n_npzd(n)) = 0.0
         enddo
       enddo
      enddo

      if (enable_co2) then
c    calculate CO2 surface flux
       if (enable_pco2_series) then
        call forcing_interp(pco2series_index,p1,p2,f1,f2)
        pco2_air=pco2series(p1)*f1+pco2series(p2)*f2
       else
        pco2_air = pco2_air_fixed
       endif

       do j=js,je
#ifdef use_wanninkhof99
        call wanninkhof99(
#else
        call wanninkhof(
#endif
     &    t(is,1,j,1,tau),stf_rest(is,j,n_co2),kw(is),ie-is+1)

        if (init.or. mod(itt,pco2_update_int)==0)  then
#ifdef account_for_slp_in_co2_flux
         slp(is:ie)=stf_clim(is:ie,j,n_co2)
         where( tmask(is:ie,1,j) == 0.0 ) slp(is:ie)=1.0
#else
         slp(is:ie)=1013.25
#endif

#ifdef account_for_seasonality_in_pco2

       if (get_calendar_type() == THIRTY_DAY_MONTHS ) then
         f1 = get_day_of_year(get_current_time())*1.0/360.
       elseif (get_calendar_type() == NO_LEAP ) then
         f1 = get_day_of_year(get_current_time())*1.0/365.
       else 
         print*,'ERROR: wrong calendar type '
         call halt_stop(' in npzd_flux (file npzd_model.F)')
       endif
c
c    Seasonality in atmos. pCO2 depends on latitude
c    and there is also a constant latitudinal gradient
c    (after Conway et al 1994, JGR)
c
        pco2_atmos =  pco2_air
     &   +(7.6558*(tanh( (yt(j)-24.7554)/ 26.5702 ) +1) +1.3971)
     &      *(0.9626*sin(f1*2*pi-0.1349)
     &       -0.2708*sin(f1*4*pi-0.5242)   )
     &   -yt(j)**3*2.6829e-06 + yt(j)**2*3.3458e-05
     &   +yt(j)*0.0353 + 1.2472
#else
         pco2_atmos = pco2_air
#endif

#ifdef use_fast_co2_subr
         call co2calc_fast(t(is,1,j,1,tau),t(is,1,j,2,tau),
#else
         call co2calc(t(is,1,j,1,tau),t(is,1,j,2,tau),
#endif
     &                t(is,1,j,n_co2,taum1),t(is,1,j,n_nutrient,tau),
     &                slp(is),ie-is+1,pco2_atmos,
     &                diff_co2(is,j),pco2_water(is,j))
        endif
        call schmidt_co2(t(is,1,j,1,tau),sc(is),ie-is+1)
        do i=is,ie
         kw(i) = kw(i)*sqrt(660.0/sc(i))/3600.0
         stf(i,j,n_co2)= kw(i) * diff_co2(i,j) *1000 ! dco2 comes in mol/m^3 instead of mmol/m^3
        enddo
       enddo ! j

c
c    For all rigid lid models we must take into account the concentration
c    dilution effect of the salt flux
c    but add this term as a source term in the first model
c    level instead to the surface fluxes.
c    This is done in routine npzd_model.
c

      endif ! enable_co2

      if (enable_o2) then
c    calculate O2 surface flux
       do j=js,je
#ifdef use_wanninkhof99
        call wanninkhof99(
#else
        call wanninkhof(
#endif
     &     t(is,1,j,1,tau),stf_rest(is,j,n_o2),kw(is),ie-is+1)
        call schmidt_o2(t(is,1,j,1,tau),sc(is),ie-is+1)
        call eqconc_o2(t(is,1,j,1,tau),t(is,1,j,2,tau), 
     &                     o2_sat(is,j),ie-is+1)
        do i=is,ie
         kw(i) = kw(i)*sqrt(660.0/sc(i))/3600.0
         stf(i,j,n_o2)= kw(i) * (o2_sat(i,j) -t(i,1,j,n_o2,taum1) )
        enddo
       enddo ! j
      endif ! enable_o2

      if (init) init = .false.
      end subroutine npzd_sflx


      subroutine npzd_model
c---------------------------------------------------------------
c   time step the NPZD model
c---------------------------------------------------------------
      use spflame_module
      use passive_tracer_module
      use npzd_module
      implicit none
      integer :: i,k,j,n,is,ie,js,je
      real :: daylen,dtsb,cobeta
      real :: fx1,fx2,fx3,fx4,fu1,fu2,vpbio,radbio
      real :: adv_fb(is_pe:ie_pe,0:km)
      real :: betaf(is_pe:ie_pe,0:km+1)
      real :: rayb(0:km)
      real :: avej(is_pe:ie_pe,km)
      real :: biotr(is_pe:ie_pe,km,nt_npzd),dnut(is_pe:ie_pe,km)
      real :: ddet(is_pe:ie_pe,km),biodoc,drem(is_pe:ie_pe,km)
      real :: o2mask(is_pe:ie_pe,km)
      real :: bion,biop,bioz,biod,u_npz,g_npz
      real :: denom
      real :: s_ref,dic_ref,exu,resp,cn2
      real :: maxtr(nt_npzd)
      real :: mintr(nt_npzd)
      integer :: maxtr_ijk(nt_npzd,3)
      integer :: mintr_ijk(nt_npzd,3)

      integer :: nn_o2,nn_co2,nn_doc  ! pointer to CO2, O2, DOC in nt_npzd

      logical, save :: init = .true.

      k=4 ! construct pointer to CO2,O2 in nt_npzd
      if (enable_co2) then
          nn_co2=k+1; k=k+1  ! always # 5
      endif
      if (enable_o2) then
          nn_o2=k+1; k=k+1  ! sometimes 5, else 6
      endif
      if (enable_doc) then
          nn_doc=k+1; k=k+1  ! it depends
      endif

      is=max(is_pe,2); ie=min(ie_pe,imt-1)
      js=max(2,js_pe); je=min(je_pe,jmt-1)
c
c   treatment of open boundaries: 
c   deactivate  biol. model there, cut off negative values
c
      if (my_blk_j == 1 .and. enable_obc_south) then
       j=2 ; js=3
       do n=1,nt_npzd
        do k=1,km
         do i=is,ie
           t(i,k,j,n_npzd(n),taup1) = max(0.0,t(i,k,j,n_npzd(n),taup1))
         enddo
        enddo
       enddo
      endif

      if (my_blk_j == n_pes_j .and. enable_obc_north) then
       j=jmt-1 ; je=jmt-2
       do n=1,nt_npzd
        do k=1,km
         do i=is,ie
           t(i,k,j,n_npzd(n),taup1) = max(0.0,t(i,k,j,n_npzd(n),taup1))
         enddo
        enddo
       enddo
      endif

      if ( my_blk_i == 1 .and. enable_obc_west) then
       i=2; is=3
       do n=1,nt_npzd
        do j=js,je
         do k=1,km
           t(i,k,j,n_npzd(n),taup1) = max(0.0,t(i,k,j,n_npzd(n),taup1))
         enddo
        enddo
       enddo
      endif

      if ( my_blk_i == n_pes_i .and. enable_obc_east) then
       i=imt-1; ie=imt-2
       do n=1,nt_npzd
        do j=js,je
         do k=1,km
           t(i,k,j,n_npzd(n),taup1) = max(0.0,t(i,k,j,n_npzd(n),taup1))
         enddo
        enddo
       enddo
      endif

#ifdef test_minmax
c
c   test global min/max before biol. model
c
      maxtr=-1e15; mintr= 1e15; maxtr_ijk = -1; mintr_ijk = -1
      do n=1,nt_npzd
       do j=js,je
        do k=1,km
          do i=is,ie
            if (maxtr(n) < t(i,k,j,n_npzd(n),taup1) ) then
              maxtr(n) = t(i,k,j,n_npzd(n),taup1)
              maxtr_ijk(n,:) = (/i,j,k/)
            endif
            if (mintr(n) > t(i,k,j,n_npzd(n),taup1) ) then
              mintr(n) = t(i,k,j,n_npzd(n),taup1)
              mintr_ijk(n,:) = (/i,j,k/)
            endif
          enddo
        enddo
       enddo
       print'(a,i3,a,f12.8,a,i3,i3,i3)',
     &   ' min tr(phys) ',n,' : ' ,mintr(n),' (i,j,k)=',mintr_ijk(n,:)
      enddo 
#endif

      if (enable_diag_npzd) dnut_diag=0.0
c
c   time step within NPZD model
c
      dtsb=c2dt/float(ts_npzd)  
      do j=js,je

       call get_daylen(j,daylen,cobeta)
       daylen=max(1.e-12,daylen)
c
c    compute daily averaged light-limited growth rate
c    analytical integration over layer thickness and day
c    after Evans and Parslow (1985)
c
       avej(:,:)= 0.0
       fx1 = 0.5*daylen*daylen/alphabio
       rayb(0) = 1.0
       do k=1,km 
        rayb(k)= exp(-zw(k)/(cobeta*dbio))
       enddo

       do k=1,kmeuph
        do i=is,ie
          radbio = max(1.0,parbio*stf_clim(i,j,n_npzd(1)))
          vpbio = abio*bbio**(cbio*t(i,k,j,1,tau))
          fx3 = fx1*vpbio/(radbio*rayb(k-1))
          fx4 = fx1*vpbio/(radbio*rayb(k))
          fu1 = sqrt(fx3**2.+daylen**2.)
          fu1 = fu1-daylen*log((daylen+fu1)/fx3)
          fu2 = sqrt(fx4**2.+daylen**2.)
          fu2 = fu2-daylen*log((daylen+fu2)/fx4)
          avej(i,k) = -1.*vpbio/dzt(k) * cobeta*dbio
     &                    *(fu1 - fu2 - fx3 + fx4)
        enddo
       enddo
c
c    transfer tracer at tau-1 to temporary arrays
c
       do n=1,nt_npzd
        do k=1,km
          do i=is,ie
            biotr(i,k,n) = t(i,k,j,n_npzd(n),taum1) 
          enddo
        enddo
       enddo
c
c    begin time stepping over 2*dt in NPZD model with Euler forward
c
       do n = 1,ts_npzd
c
c     This is the NPZD model:
c      (P: phytoplankton, Z: Zooplankton, N: Nitrate and D: Detritus)
c
c      dP/dt = u(N,Temp.,Light) P - p_P P - g(P) P Z 
c
c      dZ/dt = a g(P) P Z - d Z - p_Z Z^2
c
c      dN/dt = r D + d Z - u(N,Temp.,Light) P  [ + r_d DOC ]
c
c      dD/dt = (1-s)[ (1-a) g(P) P Z + p_P P + p_Z Z^2] -r D + w_D dD/dz
c
c     optional (for s>0) is dissolved organic carbon (DOC):
c      dDOC/dt =   s[ (1-a) g(P) P Z + p_P P + p_Z Z^2] -r_d DOC  
c
#ifdef denitrification
        o2mask=1.0
        if (enable_o2) then
c
c     If O2 falls below a small threshold, remineralisation
c     is done using NO3 instead of O2. Therefore account for that
c     process by subtracting 6.67 r D from N budget. The grazing function
c     for activity of zoo plankton is also set to zero in that case.
c     Contruct a mask to handle denitrification
c
         do k=1,km
          do i=is,ie
            if (biotr(i,k,nn_o2) <= o2_threshold) 
     &        o2mask(i,k)=1.0-tmask(i,k,j)
          enddo
         enddo
        endif
#endif

        do k=1,km
         do i=is,ie
c
#ifdef prevent_tracer_from_zero
          bion = max(0.0,biotr(i,k,1))
          biop = max(0.0,biotr(i,k,2))
          bioz = max(0.0,biotr(i,k,3))
          biod = max(0.0,biotr(i,k,4))
#else
          bion = biotr(i,k,1)
          biop = biotr(i,k,2)
          bioz = biotr(i,k,3)
          biod = biotr(i,k,4)
#endif
c
c   -- phytoplankton equation
c
c      use Liebigs Law of the Minimum (Liebig, 1845) for growth rate
c      (minimum of light-limited and nutrient limited growth rates; 
c      although chlorophyll is not explicitly considered, this will 
c      later allow for a diagnostic determination of a Chl:N ratio 
c      depending on light- or nutrient-limited growth. 
c      --> Hurtt and Armstrong, 1996)
c      saturation growth rate (infinite light, infinite nutrients)
          vpbio = abio*bbio**(cbio*t(i,k,j,1,tau))
c      growth rate
          u_npz = min(avej(i,k),vpbio*bion/(rk1bio+bion))
c      grazing function
          g_npz = gbio*epsbio*biop/(gbio+epsbio*biop*biop)
#ifdef denitrification
     &         *o2mask(i,k)
c      In case of low oxygen, activity of zoo plankton is set to zero 
#endif
c
          biotr(i,k,2) = biotr(i,k,2) + dtsb * biop * (
     &         u_npz - phiphy - g_npz*bioz )*tmask(i,k,j)
c
c   -- zooplankton equation
c
          biotr(i,k,3) = biotr(i,k,3) + dtsb * bioz * (
     &         a_npz*g_npz*biop - d_npz - phizoo*bioz )*tmask(i,k,j)
c
c   -- nutrient equation
c
          dnut(i,k)    =   (remina*biod + d_npz*bioz - u_npz*biop)
     &                         *dtsb*tmask(i,k,j)
          biotr(i,k,1) = biotr(i,k,1) + dnut(i,k)
c
c   -- detritus equation
c
          ddet(i,k) = (1.-a_npz)*g_npz*biop*bioz 
     &                + phiphy*biop + phizoo*bioz*bioz
          drem(i,k) = remina*biod
          biotr(i,k,4) = biotr(i,k,4) 
     &       + dtsb * (ddet(i,k) - drem(i,k) )*tmask(i,k,j)
         enddo
        enddo

#ifdef denitrification
c
c     If O2 falls below a small threshold, remineralisation
c     is done using NO3 instead of O2. Therefore account for that
c     process by subtracting 6.67 r D from N budget. If NO3 falls
c     below zero, other processes are taking over. Set NO3 to zero
c     in that case.
c
        do k=1,km
         do i=is,ie
          biotr(i,k,1) = max(0.0,   biotr(i,k,1)-
     &         (1.0-o2mask(i,k))*drem(i,k)*6.67*dtsb*tmask(i,k,j))
         enddo
        enddo
c
c     note that the array dnut is not affected by this process
c
#endif



        if (enable_co2) then
c
c     change CO2 according to a constant Redfield ratio 
c
         do k=1,km
          do i=is,ie
           biotr(i,k,nn_co2)=biotr(i,k,nn_co2)+redfield_CN*dnut(i,k)
          enddo
         enddo
        endif

        if (enable_o2) then
c
c     Change oxygen concentration due to biological activity
chd   oxygen concentrations are in ml/l; donstant
chd   redfieldratio of N/O2=1/10 (ocmip) is assumed;
chd   it is: mlO2/l = 10*22.3916/1000 mmolN/m^3 
c
         do k=1,km
          do i=is,ie
           biotr(i,k,nn_o2) = biotr(i,k,nn_o2) - redfield_ON*dnut(i,k)
#ifdef denitrification
     &     + drem(i,k)*dtsb*tmask(i,k,j)*(1.0-o2mask(i,k))*redfield_ON
c      do not change oxygen by remineralisation anymore if o2<o2_threshold
#endif
          enddo
         enddo

        endif


        if (enable_doc) then
c
c   -- DOC equation
c 
         do k=1,km
          do i=is,ie
c          (0.05+50*tanh(N*50))
           bion  = biotr(i,k,1)
           biop  = biotr(i,k,2)
c           exu  = s_doc*biotr(i,k,2)
           exu = 5+99*tanh(bion*50)
           exu = 1./exu/86400.*biop
           resp = r_doc*biotr(i,k,nn_doc)
           biotr(i,k,nn_doc)=biotr(i,k,nn_doc)
     &             +(exu-resp)*dtsb*tmask(i,k,j)
c        correct other tracers:
c            DIN
           biotr(i,k,1)=biotr(i,k,1)+dtsb*tmask(i,k,j)*resp
c            P_N
           biotr(i,k,2)=biotr(i,k,2)-dtsb*tmask(i,k,j)*exu
c             O2
           biotr(i,k,nn_o2)=biotr(i,k,nn_o2)-dtsb*tmask(i,k,j)*
     &              (resp)*redfield_ON
c             CO2
c           cn2=redfield_CN2-(redfield_CN2-redfield_CN)*tanh(bion)
           biotr(i,k,nn_co2)=biotr(i,k,nn_co2)+dtsb*tmask(i,k,j)*
     &      ( redfield_CN2*resp - (redfield_CN2-redfield_CN)*exu )
          enddo
         enddo
        endif


        if (enable_diag_npzd) then
c
c     diagnose sinks/sources of nutrients to be averaged
c
         do k=1,km
          do i=is,ie
           dnut_diag(i,k,j)=dnut_diag(i,k,j)+dnut(i,k)
          enddo
         enddo
        endif

       enddo ! time steps
c
c    transfer from temporary arrays and add tendency to tracer at tau +1
c
       do n=1,nt_npzd
        do k=1,km
          do i=is,ie
#ifdef prevent_tracer_from_zero
           t(i,k,j,n_npzd(n),taup1) = max(0.0,
     &         t(i,k,j,n_npzd(n),taup1)  +
     &       ( biotr(i,k,n) - t(i,k,j,n_npzd(n),taum1 ) ) )
#else
           t(i,k,j,n_npzd(n),taup1) = 
     &         t(i,k,j,n_npzd(n),taup1)  +
     &       ( biotr(i,k,n) - t(i,k,j,n_npzd(n),taum1 ) ) 
#endif
          enddo
        enddo
       enddo
c
c    add tendency due to sinking of detritus
c
       n=n_npzd(4)
       do k=1,km-1
        do i=is,ie
	 adv_fb(i,k)=-w_detr*(t(i,k,j,n,tau) + t(i,k+1,j,n,tau))
        enddo
       enddo

c    no flux boundary conditions
       do i=is,ie
        adv_fb(i,0)       =0.0
        adv_fb(i,km)      =0.0
        adv_fb(i,kmt(i,j))=0.0
       enddo

c    compute delimiter for outgoing flux
       betaf=2.
       do k=1,km
        do i=is,ie
         denom =      -min(0.,adv_fb(i,k  )*dzt2r(k))
         denom = denom+max(0.,adv_fb(i,k-1)*dzt2r(k))

         denom = denom*c2dt
         denom = sign(1.,denom)*max(1.e-15,abs(denom)) ! should be positive anyway

         betaf(i,k)=t(i,k,j,n,taum1)/denom
        enddo
       enddo

c    delimit vertical flux
       do k=0,km
        do i=is,ie
           adv_fb(i,k)=
     &         min(1.,betaf(i,k  ))*min(0.,adv_fb(i,k))
     &       + min(1.,betaf(i,k+1))*max(0.,adv_fb(i,k))
        enddo
       enddo

       do k=1,km
        do i=is,ie
	  t(i,k,j,n,taup1) = t(i,k,j,n,taup1) - 
     &           c2dt*tmask(i,k,j)*(adv_fb(i,k-1) - adv_fb(i,k))
#ifdef partial_cell
     &               *.5/dht(i,k,j)
#else
     &               *dzt2r(k)
#endif
        enddo
       enddo

#ifdef detritus_sedimentation
c    put any detritus at bottom into sediment
       n=n_npzd(4)
       do i=is,ie
        t(i,km,j,n,taup1) = 0.0
       enddo
       do k=1,km-1
        do i=is,ie
         t(i,k,j,n,taup1) = tmask(i,k+1,j)*t(i,k,j,n,taup1)
        enddo
       enddo
#endif
#ifdef detritus_slow_sedimentation
c    different, sligthly more realistic formulation 
       n=n_npzd(4)
       do k=1,km-1
        do i=is,ie
         t(i,k,j,n,taup1) = t(i,k,j,n,taup1) - c2dt*
     &           t(i,k,j,n,taum1)*remina_slow*
     &             (1-tmask(i,k+1,j))*tmask(i,k,j)
        enddo
       enddo
       k=km
       do i=is,ie
         t(i,k,j,n,taup1) = t(i,k,j,n,taup1) - c2dt*
     &           t(i,k,j,n,taum1)*remina_slow*tmask(i,k,j)
       enddo
#endif

       if (enable_co2 .and. .not. enable_freesurf) then
c
c    For all rigid lid models we must take into account the concentration
c    dilution effect of the salt flux
c

c    stf(:,:,2)=dzt(1)/tau ( S* - SSS)
c    after OMIP abiotic-Howto:  E-P =  (S*-SSS)/S_ref *dz/tau
c                         and:  F_dic = DIC_ref*(E-P)
c  
c     (stf(,,2)*1000)/S_ref = E-P
c
c    for S*>SSS :  stf(,,2)>0
c    that means put salinity into ocean, take freshwater out of ocean
c    DIC concentration should increase ->  F_dic >0
c     -> should work 
c
c       do j=js,je
c        do i=is,ie
c         stf(i,j,n_co2)= stf(i,j,n_co2)
c     &                 + dic_ref*(stf(i,j,2)*1000)/s_ref
c        enddo
c       enddo ! j
         s_ref=35; dic_ref=2000.
         do i=is,ie
           s_ref=t(i,1,j,2,tau)*1000.0+35.0
           dic_ref = t(i,1,j,n_co2,tau)
           t(i,1,j,n_co2,taup1) = t(i,1,j,n_co2,taup1) + c2dt*
     &           dic_ref*(stf(i,j,2)*1000)/s_ref/dzt(1)*tmask(i,1,j)
         enddo 
       endif

      enddo ! j loop

#ifdef test_minmax
c
c   global min/max
c
      maxtr=-1e15; mintr= 1e15; maxtr_ijk = -1; mintr_ijk = -1
      do n=1,nt_npzd
       do j=js,je
        do k=1,km
          do i=is,ie
            if (maxtr(n) < t(i,k,j,n_npzd(n),taup1) ) then
              maxtr(n) = t(i,k,j,n_npzd(n),taup1)
              maxtr_ijk(n,:) = (/i,j,k/)
            endif
            if (mintr(n) > t(i,k,j,n_npzd(n),taup1) ) then
              mintr(n) = t(i,k,j,n_npzd(n),taup1)
              mintr_ijk(n,:) = (/i,j,k/)
            endif
          enddo
        enddo
       enddo
      enddo 

      print*,' max tracer : ',maxtr
      do n=1,nt_npzd
       print'(a,i3,a,f12.8,a,i3,i3,i3)',
     &   ' min tracer   ',n,' : ' ,mintr(n),' (i,j,k)=',mintr_ijk(n,:)
      enddo
#endif


       if (enable_diag_npzd) then
         dnut_diag= dnut_diag/c2dt  
         call npzd_averages()
       endif

      end subroutine npzd_model



      subroutine get_daylen(j,daylen,cobeta)
      use spflame_module
      implicit none
      integer, intent(in) :: j
      real, intent(out)   :: daylen,cobeta
      real :: fx1,fx2,timeyear,declin
c
c   solar declination 
c
      timeyear = get_day_of_year(get_current_time()) /365.
      fx1    = 2.0*pi*(timeyear-int(timeyear))
      declin = 0.006918
     &      -0.399912*cos(   fx1)+0.070257*sin(   fx1)
     &      -0.006758*cos(2.*fx1)+0.000907*sin(2.*fx1)
     &      -0.002697*cos(3.*fx1)+0.001480*sin(3.*fx1)
c
c   compute solar angle at noon (and assume that this is the equivalent
c   daily averaged incidence angle for direct+diffuse radiation)
c   1. cobeta is cos(incidence angle of solar radiation at noon)
c
      fx1    = pi/180.*yt(j)
      cobeta = max(0.,sin(fx1)*sin(declin)+cos(fx1)*cos(declin))
c
c   2. cobeta is cos(incidence angle in water)
c      using Snells law: sin(water)=sin(air)/1.33
c
      cobeta = sqrt(1.-(1.-cobeta**2.)/1.33**2.)
c
c   finally,length of day
c
      fx2    = max(-1., min( 1.,-tan(fx1)*tan(declin)) )
      daylen = acos(fx2)/pi
      end subroutine get_daylen




      subroutine schmidt_co2(t,sc,len)
c-----------------------------------------------------------------------
c    Computes the Schmidt number of CO2 in seawater using the 
c    formulation presented by Wanninkhof (1992, J. Geophys. Res., 97,
c    7373-7382).  Input is temperature in deg C.
c-----------------------------------------------------------------------
      implicit none
      integer :: len
      real    :: t(len),sc(len)
      integer :: i
      do i=1,len
       sc(i) = 2073.1 - 125.62*t(i) + 3.6276*t(i)*t(i) 
     &       - 0.043219*t(i)*t(i)*t(i)
      enddo
      end subroutine schmidt_co2




      subroutine schmidt_o2(t,sc,len)
c-----------------------------------------------------------------------
c   Computes the Schmidt number of oxygen in seawater using the
c   formulation proposed by Keeling et al. (1998, Global Biogeochem.
c   Cycles, 12, 141-163).  Input is temperature in deg C.
c-----------------------------------------------------------------------
      implicit none
      integer :: len
      real    :: t(len),sc(len)
      integer :: i
      do i=1,len
       sc(i) = 1638.0-81.83*t(i)+1.483*t(i)*t(i) 
     &             - 0.008004*t(i)*t(i)*t(i)
      enddo
      end subroutine schmidt_o2



      subroutine eqconc_o2(t,s,osa,len)
c-----------------------------------------------------------------------
C   Computes the oxygen saturation concentration at 1 atm 
c   total pressure in ml/l given the temperature 
c   (t, in deg C) and the salinity (s, in permil, or psu). 
C   FROM GARCIA AND GORDON (1992), LIMNOLOGY and OCEANOGRAPHY.
C   THE FORMULA USED IS FROM PAGE 1310, EQUATION (8).
C   *** NOTE: THE "A3*TS^2" TERM (IN THE PAPER) IS INCORRECT. ***
C   *** IT SHOULD not BE THERE.                                ***
C   osa IS DEFINED BETWEEN T(freezing) <= T <= 40(deg C) AND
c   0 permil <= S <= 42 permil
C   CHECK VALUE:  T = 10.0 deg C, S = 35.0 permil,osa = 0.282015 mol/m^3
c     Convert from ml/l to mol/m^3  osa = o2sato/22391.6*1000.0
c-----------------------------------------------------------------------
      implicit none
      integer :: len
      real :: t(len), s(len), osa(len)
      real :: TT, TK, TS, TS2, TS3, TS4, TS5, CO ,SS
      integer :: i
      real :: A0= 2.00907, A1= 3.22014, A2= 4.05010 
      real :: A3= 4.94457, A4=-2.56847E-1, A5= 3.88767 
      real :: B0=-6.24523E-3, B1=-7.37614E-3
      real :: B2=-1.03410E-2, B3=-8.17083E-3, Cc0=-4.88682E-7
      do i=1,len
         TT  = 298.15-t(i)
         TK  = 273.15+t(i)
         TS  = LOG(TT/TK)
         TS2 = TS*TS
         TS3 = TS2*TS
         TS4 = TS3*TS
         TS5 = TS4*TS
         ss  = s(i)*1000.0+35.0
         CO  = A0 + A1*TS + A2*TS2 + A3*TS3 + A4*TS4 + A5*TS5
     $        + ss*(B0 + B1*TS + B2*TS2 + B3*TS3)
     $        + Cc0*(ss*ss)
         osa(i) = EXP(CO)
      enddo
      end subroutine eqconc_o2
