
c=======================================================================
c
c     Setup for a simple biogeochemical model:
c      (P: phytoplankton, N: Nitrate and D: Detritus)
c
#define enable_zooplankton
c
#ifdef enable_zooplankton
c      dN/dt = r D + d Z - u(N,Temp.,Light) P  
c      dP/dt = u(N,Temp.,Light) P - p_P P       - g(P) P Z 
c      dD/dt =  p_P P -r D + w_D dD/dz          + (1-a) g(P) P Z + p_Z Z^2
c      dZ/dt = a g(P) P Z - d Z - p_Z Z^2
#else
c      dN/dt = r D  - u(N,Temp.,Light) P  
c      dP/dt = 0.2 u(N,Temp.,Light) P - p_P P      
c      dD/dt =  p_P P -r D + w_D dD/dz   +0.8 u(N,Temp.,Light) P
#endif
c
c     linked in the code in ../extra_modules/tracer.F
c=======================================================================

      module bgc_module
      use cpflame_module
      use tracer_module
      implicit none
c ----------------------------------
c    Parameters for the biogeochemical model
c ----------------------------------
#ifdef enable_zooplankton
      integer :: nut=1,phy=2,zoo=3,det=4,oxy=5,dic=6
#else
      integer :: nut=1,phy=2,det=3,oxy=4,dic=5
#endif
      real, allocatable :: solar(:,:) ! solar radiation in W/m^2
      real, allocatable :: sss(:,:)   ! sea surface salinity in psu
      real, allocatable :: u10(:,:)  ! wind speed in 10 m height in m/s
      real, allocatable :: u_nut(:,:,:)  ! phytoplankton growth rate
      real, allocatable :: u_light(:,:,:)  ! phytoplankton growth rate
      real, allocatable :: o2_sflx(:,:), o2_sat(:,:)
      real, allocatable :: dic_sflx(:,:) ,pco2(:,:)
      real, parameter :: phiphy   = 0.03 /86400.  ! phytoplankton death rate
      real, parameter :: remina =   0.05*0.8 /86400. ! detritus remineralisation rate 
      real, parameter :: w_detr   = -5.0  /86400.  ! sinking velocity of detritus at surface
      real, parameter :: w_detr2  = -5.0  /86400.  ! sinking velocity of detritus at depth
      real, parameter :: abio     = 0.6 /86400.  ! saturation growth rate
      real, parameter :: bbio     = 1.066 ! constants for saturation growth rate
      real, parameter :: cbio     = 1.0
#ifdef enable_zooplankton
      real, parameter  :: phizoo   = 0.20 /86400.   ! zooplankton mortality rate
      real, parameter  :: d_npz    = 0.03 /86400.  ! zooplankton excretation rate
      real, parameter  :: a_npz    = 0.75             ! zooplanktion grazing fraction
      real, parameter  :: gbio      = 2.0  /86400.   ! zooplankton grazing rate
      real, parameter  :: epsbio   = 1.0  /86400.   ! zooplankton grazing rate
#endif
      real, parameter :: rk1bio   = 0.5
      real, parameter :: dbio     = 25.0  ! penetration depth of solar radiation in meters
      real, parameter :: alphabio = 0.025 /86400. ! add docu for biological parameters
      real, parameter :: parbio   = 0.43

      real, parameter :: c_ON   = 10.0
      real, parameter :: c_CN   = 6.6
      real, parameter :: max_dt = 900.0 ! maximal time step in biog. model
      integer :: ts_npz                       ! additional time steps in bio. module 
      end module bgc_module


      subroutine tracer_set_number
c ----------------------------------
c     set number of tracers
c ----------------------------------
      use bgc_module
      implicit none
#ifdef enable_zooplankton
      nt=6
#else
      nt=5
#endif
      ts_npz=max(1,nint(dt_in*2. / max_dt))
      if (my_pe==0) then
       print*,' there are ',nt,' biogeochem. variables '
       print*,' biogeochemical model will do ',ts_npz,' time steps'
       print*,' time step in biochem. model will be ',
     &          dt_in/ts_npz,'sec.'

       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 
#ifdef enable_zooplankton
       print*,' phizoo  = ',phizoo
       print*,' a_npz   = ',a_npz
       print*,' d_npz   = ',1./d_npz/86400.,' days'
       print*,' gbio     = ', gbio,' days' 
       print*,' epsbio    = ',epsbio,' days'
#endif
       print*,' parbio   = ',parbio
       print*,' rk1bio   = ',rk1bio
       print*,' dbio     = ', dbio,' m'
       print*,' phiphy   = ',1./phiphy/86400.,' days'
       print*,' remina   = ',1./remina/86400.,' days'
       print*,' w_detr   = ',w_detr *86400. ,' m/day'
       print*,' w_detr2  = ',w_detr2*86400. ,' m/day'
       print*,' c_ON     = ',c_ON
       print*,' c_CN     = ',c_CN
       print*,' '
      endif
      end subroutine tracer_set_number


      subroutine tracer_initial_conditions
c ----------------------------------
c     initial conditions for tracer
c ----------------------------------
      use bgc_module
      implicit none
      integer :: i,j,k
      real :: L_y,y2
      allocate( u_nut(imt,jmt,km) ) ; u_nut = 0.0
      allocate( u_light(imt,jmt,km) ) ; u_light = 0.0
      allocate( solar(imt,jmt) ) ; solar = 0.0
      allocate( sss(imt,jmt) ) ; sss = 0.0
      allocate( u10(imt,jmt) ) ; u10 = 0.0
      allocate( o2_sflx(imt,jmt),dic_sflx(imt,jmt) )  
      o2_sflx = 0.0;  dic_sflx=0.0
      allocate( o2_sat(imt,jmt),pco2(imt,jmt) ) 
      o2_sat = 0.0;   pco2=0.0
c ----------------------------------
c     specify surface boundary conditions
c ----------------------------------
      L_y = yt(jmt)-yt(1)
      y2=L_y*0.75
      do j=1,jmt
       do i=1,imt
        solar(i,j)=200.
        if (yt(j)<L_y/2.0) solar(i,j)=200*yt(j)/(L_y/2.0)
        if (yt(j)>y2)      solar(i,j)=200*(1-(yt(j)-y2)/(yt(jmt)-y2))
        sss(i,j)=35.0
        u10(i,j)=10.0
       enddo
      enddo
c ----------------------------------
c     specify initial conditions
c ----------------------------------
      do k=1,km
       do j=1,jmt
        do i=1,imt
         tr(i,j,k,:,phy)=0.0001*maskT(i,j,k)
#ifdef enable_zooplankton
         tr(i,j,k,:,zoo)=0.0001*maskT(i,j,k)
#endif
         tr(i,j,k,:,det)=0.0000*maskT(i,j,k)
         if (k<km-2) tr(i,j,k,:,nut)=30*maskT(i,j,k)
         tr(i,j,k,:,oxy)=500*maskT(i,j,k)
         tr(i,j,k,:,dic)=2100*maskT(i,j,k)
        enddo
       enddo
      enddo
      if (my_pe==0) call init_bgc_diag
      end subroutine tracer_initial_conditions


      subroutine tracer_sources
c ----------------------------------
c     Sources and sinks for tracers
c ----------------------------------
      use bgc_module
      implicit none
      integer :: i,j,k,js,je,n
      real :: fxa,fxb,dtsb,dnut(imt,jmt,km),tnut,tphy,tzoo,tdet,graz
      real :: sst(imt),sc(imt),kw(imt),dco2(imt),sdic(imt)
      js=max(2,js_pe); je = min(je_pe,jmt-1)
c ----------------------------------
c    interior sources of biogeochem. variables
c ----------------------------------
      dtsb=c2dt/float(ts_npz)  
      call light_limitation
      dnut(:,js_pe:je_pe,:)=0.
      do n = 1,ts_npz
       call phyto_growth_rate
       do k=2,km-1
        do j=js,je
         do i=2,imt-1
          fxa=min(u_nut(i,j,k),u_light(i,j,k))                             ! growth rate
          tnut=tr(i,j,k,taup1,nut)
          tdet=tr(i,j,k,taup1,det)
          tphy=tr(i,j,k,taup1,phy)
#ifdef enable_zooplankton
          tzoo=tr(i,j,k,taup1,zoo)
          graz = gbio*epsbio*tphy/(gbio+epsbio*tphy*tphy)    ! grazing
c        dP/dt = u(N,Temp.,Light) P - p_P P       - g(P) P Z 
          tr(i,j,k,taup1,phy)= tphy + maskT(i,j,k)*dtsb*
     &             tphy*(fxa - phiphy - graz*tzoo) 
c        dZ/dt = a g(P) P Z - d Z - p_Z Z^2
          tr(i,j,k,taup1,zoo) = tzoo + maskT(i,j,k)*dtsb* 
     &             tzoo*(a_npz*graz*tphy - d_npz - phizoo*tzoo )
c        dN/dt = r D + d Z - u(N,Temp.,Light) P  
          fxb =remina*tdet +d_npz*tzoo -fxa*tphy                    ! nutrient change
          dnut(i,j,k)=dnut(i,j,k)+fxb                                          ! accumulate nutrient change
          tr(i,j,k,taup1,nut)= tnut + maskT(i,j,k)*dtsb*fxb
c        dD/dt =  p_P P -r D + w_D dD/dz          + (1-a) g(P) P Z + p_Z Z^2
          tr(i,j,k,taup1,det)= tdet + maskT(i,j,k)*dtsb* 
     &         (phiphy*tphy-remina*tdet
     &         +(1-a_npz)*graz*tphy*tzoo+phizoo*tzoo**2 )
#else
c        dP/dt = 0.2*u(N,Temp.,Light) P - p_P P     
          tr(i,j,k,taup1,phy)= tphy + maskT(i,j,k)*dtsb*
     &             tphy*(0.2*fxa - phiphy) 
c        dN/dt = r D - u(N,Temp.,Light) P  
          fxb =remina*tdet -fxa*tphy                    ! nutrient change
          dnut(i,j,k)=dnut(i,j,k)+fxb                       ! accumulate nutrient change
          tr(i,j,k,taup1,nut)= tnut + maskT(i,j,k)*dtsb*fxb
c        dD/dt =  p_P P -r D + w_D dD/dz   +0.8 u(N,Temp.,Light) P
          tr(i,j,k,taup1,det)= tdet + maskT(i,j,k)*dtsb* 
     &         (phiphy*tphy-remina*tdet+0.8*fxa*tphy)
#endif
         enddo
        enddo
       enddo
      enddo
c ----------------------------------
c    other tracers
c ----------------------------------
      do k=2,km-1
       do j=js,je
        do i=2,imt-1
          tr(i,j,k,taup1,oxy)=max(0.,tr(i,j,k,taup1,oxy)+
     &          maskT(i,j,k)*c2dt*(-c_ON*dnut(i,j,k)) )
          tr(i,j,k,taup1,dic)=max(0.,tr(i,j,k,taup1,dic)+
     &          maskT(i,j,k)*c2dt*(c_CN*dnut(i,j,k)) )
        enddo
       enddo
      enddo
c ----------------------------------
c     sinking of organic materal
c ----------------------------------
      call detritus_sinking
c ----------------------------------
c    surface fluxes of biogeochem. variables
c ----------------------------------
       k=km-1
       do j=js,je
        sst=-b(:,j,k,tau)*1024./9.81/0.2
        sdic=tr(:,j,k,taum1,dic)
        call wanninkhof(sst,u10(:,j),kw,imt)
        call schmidt_o2(sst,sc,imt)
        call eqconc_o2(sst,sss(:,j),o2_sat(:,j),imt)
        o2_sflx(:,j)=kw*sqrt(660.0/sc)*(o2_sat(:,j)-tr(:,j,k,taum1,oxy))
        tr(:,j,k,taup1,oxy)=tr(:,j,k,taup1,oxy)+
     &                             maskT(:,j,k)*c2dt*o2_sflx(:,j)/dz
        call schmidt_co2(sst,sc,imt)
        call eqconc_co2(sst,sss(:,j),sdic,imt,dco2,pco2(:,j))
        dic_sflx(:,j)=  kw*sqrt(660.0/sc)*dco2
        tr(:,j,k,taup1,dic)=tr(:,j,k,taup1,dic)+
     &                           maskT(:,j,k)*c2dt*dic_sflx(:,j)/dz 
      enddo
c ----------------------------------
c    diagnostic output
c ----------------------------------
      fxa=0.0;fxb=0.0
      do k=2,km-1
       do j=js,je
        do i=2,imt-1
          fxa=fxa+tr(i,j,k,tau,nut)*dx**2*dz*maskT(i,j,k)
          fxa=fxa+tr(i,j,k,tau,det)*dx**2*dz*maskT(i,j,k)
          fxa=fxa+tr(i,j,k,tau,phy)*dx**2*dz*maskT(i,j,k)
#ifdef enable_zooplankton
          fxa=fxa+tr(i,j,k,tau,zoo)*dx**2*dz*maskT(i,j,k)
#endif
          fxb=fxb+tr(i,j,k,tau,oxy)*dx**2*dz*maskT(i,j,k)
        enddo
       enddo
      enddo
      call global_sum(fxa)
      call global_sum(fxb)
      fxa = fxa/1000. /1000/1000/1000
      fxb = fxb/1000. /1000/1000/1000
      if (my_pe==0) print*, 'total nutrient content ',fxa,' Giga mol'
      if (my_pe==0) print*, 'total oxygen  content ',fxb,' Giga mol'
      call sub_flush(6)
      if (snapshot_time_step.or.initial_time==current_time) then
        call diag_bgc
      endif
      end subroutine tracer_sources

c ----------------------------------
c    the rest is specific for biogeochem. model
c ----------------------------------

      subroutine light_limitation
c ----------------------------------
c    compute daily averaged light-limited growth rate
c    analytical integration over layer thickness and day
c    after Evans and Parslow (1985)
c ----------------------------------
      use bgc_module
      implicit none
      integer :: i,j,k,js,je
      real :: fxa,fxb,fx1,fx3,fx4,fu1,fu2
      real :: daylen,cobeta,radbio,rayb(km)
      js=max(2,js_pe); je = min(je_pe,jmt-1)
c    day length from astronomic functions
      call get_daylen(180.,30.0,daylen,cobeta)
      daylen=max(1.e-12,daylen)
c    depth penetration of solar insulation
      do k=1,km-1
        rayb(k)= exp(zw(k)/(cobeta*dbio))
      enddo
      rayb(km) = 1.0
      fx1 = 0.5*daylen*daylen/alphabio
      do j=js,je
       do k=2,km-1
        do i=2,imt-1
          radbio = max(1.0,parbio*solar(i,j))
          fxa = abio*bbio**(-cbio*b(i,j,k,tau)*1024./9.81/0.2)
          fx3 = fx1*fxa/(radbio*rayb(k-1))
          fx4 = fx1*fxa/(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)
          u_light(i,j,k) =fxa/dz*cobeta*dbio*(fu1 - fu2 - fx3 + fx4) 
        enddo
       enddo
      enddo
      end subroutine light_limitation


      subroutine phyto_growth_rate
c ----------------------------------
c      nutrient limited growth rate
c ----------------------------------
      use bgc_module
      implicit none
      integer :: i,j,k,js,je
      real :: fxa,fxb
      js=max(2,js_pe); je = min(je_pe,jmt-1)
      do k=2,km-1
       do j=js,je
        do i=2,imt-1
c       saturation growth rate (infinite light, infinite nutrients)
         fxa = abio*bbio**(-cbio*b(i,j,k,tau)*1024./9.81/0.2)
         fxb = tr(i,j,k,taup1,nut)/(rk1bio+tr(i,j,k,taup1,nut))
c       nutrient limited growth rate
         u_nut(i,j,k) = fxa*fxb
        enddo
       enddo
      enddo
      end subroutine phyto_growth_rate


      subroutine detritus_sinking
c ----------------------------------
c     sinking of detritus
c ----------------------------------
      use bgc_module
      implicit none
      integer :: i,j,k,js,je
      real :: adv_ft(imt,jmt,km),betaf(imt,jmt,km),denom,ww
      js=max(2,js_pe); je = min(je_pe,jmt-1)
c    add tendency due to sinking of detritus
      adv_ft(:,js:je,:)=0.0
      do k=1,km-1
        ww = w_detr+min(0.,(w_detr2-w_detr)*(zw(k)+100.)/zw(1))
       do j=js,je
        do i=2,imt-1
          adv_ft(i,j,k)=maskW(i,j,k)*
     &   0.5*ww*(tr(i,j,k,tau,det) + tr(i,j,k+1,tau,det))
         enddo
        enddo
       enddo
c    compute delimiter for outgoing flux
      betaf(:,js:je,:)=2.
      do k=2,km
       do j=js,je
        do i=2,imt-1
         denom =       max(0.,adv_ft(i,j,k)/dz)
         denom = denom-min(0.,adv_ft(i,j,k-1)/dz) 
         denom = denom*c2dt
         denom = sign(1.,denom)*max(1.e-15,abs(denom)) ! should be positive anyway
         betaf(i,j,k)=tr(i,j,k,taup1,det)/denom
        enddo
       enddo
      enddo
c    delimit vertical flux
      do k=1,km-1
       do j=js,je
        do i=2,imt-1
           adv_ft(i,j,k)=
     &         min(1.,betaf(i,j,k   ))*max(0.,adv_ft(i,j,k))
     &       + min(1.,betaf(i,j,k+1))*min(0.,adv_ft(i,j,k))
        enddo
       enddo
      enddo
      do k=2,km-1
       do j=js,je
        do i=2,imt-1
          tr(i,j,k,taup1,det) = tr(i,j,k,taup1,det) - 
     &       c2dt*maskT(i,j,k)*(adv_ft(i,j,k) - adv_ft(i,j,k-1))/dz
        enddo
       enddo
      enddo
      end subroutine detritus_sinking


      subroutine get_daylen(day,lat,daylen,cobeta)
c ----------------------------------
c     length of day as function of latitude and day of year
c ----------------------------------
      implicit none
      real, intent(in) :: day,lat
      real, intent(out)   :: daylen,cobeta
      real :: fx1,fx2,timeyear,declin
      real, parameter :: pi=3.1415
c    solar declination 
      timeyear = day/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    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)
      fx1    = pi/180.*lat
      cobeta = max(0.,sin(fx1)*sin(declin)+cos(fx1)*cos(declin))
c    2. cobeta is cos(incidence angle in water)
c      using Snells law: sin(water)=sin(air)/1.33
      cobeta = sqrt(1.-(1.-cobeta**2.)/1.33**2.)
c    finally,length of day
      fx2    = max(-1., min( 1.,-tan(fx1)*tan(declin)) )
      daylen = acos(fx2)/pi
      end subroutine get_daylen



      subroutine wanninkhof(tc,u,kw,len)
c-----------------------------------------------------------------------
c   Piston velocity in m/s for surface fluxes a la Wanninkhof , 1992
c-----------------------------------------------------------------------
      implicit none
      integer :: len,i
      real :: tc(len),u(len)   ! temperature, wind speed in m/s (input)
      real :: kw(len)          ! piston velocity (output)
      do i=1,len
       kw(i) = 0.39*u(i)*u(i)   
      enddo
      kw=kw/100./3600.   ! convert from cm/hr to m/s
      end subroutine wanninkhof


      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 mmol/m^3 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)
         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)
         osa(i)=osa(i)/22391.6*1000.0 *1.e3
      enddo
      end subroutine eqconc_o2


      subroutine eqconc_co2(t,s_in,dic_in,imt,dco2,pco2_water)
c-------------------------------------------------------------------------
c   Calculate pco2* from total alkalinity and total CO2 at
c   temperature (t), salinity (s) and "atmpres" atmosphere total pressure. 
c   using a fit of pH as a function of DIC, SALT and TEMP
c   valid range is 1940<DIC<2200 mmol/m^3
c                  -5<T<25 degC and 31.5<S<37.5 PSU
c   accuracy within this range is as good as for other routine
c--------------------------------------------------------------------------
       implicit none
c      INPUT
       integer :: imt         ! vector length
       real   :: t(imt)        ! temperature deg C
       real   :: s_in(imt)     ! salinity in PSU
       real   :: dic_in(imt)   ! diss. inorg. carbon  in mmol/m^3
c      OUTPUT
       real   :: dco2(imt)     ! difference in CO2 at sea surface and saturation CO2 in mmol/m^3
       real   :: pco2_water(imt) ! partial pressure of CO2 at sea surface  (ppmv)

      real, parameter :: rho0 = 1024.0    ! unit conversion factor
      real, parameter :: rrho02 = 1.0/rho0**2
      real, parameter :: permeg=1.e-6, permil = 1.0/rho0 ! unit conversion
      real   :: slp=1013.25      ! Sea level pressure in mbar
      real   :: pco2_air=278.0 ! atmospheric mole fraction CO2 in dry air (ppmv) 
      
c      local work variables
       real :: tk,tk100,tk1002,invtk,dlogtk,s2
       real :: htotal,s,htotal2,co2star,pco2
       real :: dic,ff,k1,k2,atmpres,pH_fit,xx,yy,zz
       real :: xx2,yy2,zz2
       integer :: i

       do i=1,imt
        s = s_in(i)
        dic = dic_in(i)*rrho02
        tk = 273.15 + t(i)
        tk100 = tk/100.0
        tk1002=tk100*tk100
        invtk=1.0/tk
        dlogtk=log(tk)
        s2=s*s
C        f = k0(1-pH2O)*correction term for non-ideality
C        Weiss & Price (1980, Mar. Chem., 8, 347-359; Eq 13 with table 6 values)
        ff = exp(-162.8301 + 218.2968/tk100  +
     &		90.9241*log(tk100) - 1.47696*tk1002 +
     &		s*(.025695 - .025225*tk100 + 
     &		0.0049867*tk1002))
C        k1 = [H][HCO3]/[H2CO3]
C        k2 = [H][CO3]/[HCO3]
C        Millero p.664 (1995) using Mehrbach et al. data on seawater scale 
        k1=10**(-1*(3670.7*invtk - 62.008 + 9.7944*dlogtk -
     &		0.0118 * s + 0.000116*s2))
        k2=10**(-1*(1394.7*invtk + 4.777 - 
     &		0.0184*s + 0.000118*s2))
c    calculate pH from estimated fit
      xx=dic_in(i)-2000.0; yy=s-35.0; zz=t(i);
      xx2=xx*xx; yy2=yy*yy; zz2=zz*zz;
      pH_fit=8.456655e+00+(-1.688880e-03)*xx 
     & +(7.390822e-02)*yy+(-1.739867e-02)*zz 
     & +(-1.864020e-06)*xx2+(-2.034953e-03)*yy2 
     & +(6.480540e-05)*zz2 
     & +(8.483447e-07)*xx*zz+(1.945036e-04)*xx*yy 
     & +(-4.480909e-05)*yy*zz 
     & +(-2.946383e-09)*xx**3+(-1.420486e-03)*yy**3 
     & +(-4.430561e-07)*zz**3 
     & +(-1.515261e-05)*xx*yy2+(5.880059e-07)*yy*xx2  
     & +(4.673103e-08)*xx*zz2+(1.772955e-08)*zz*xx2 
     & +(-1.897002e-06)*yy*zz2+(2.194509e-05)*zz*yy2 
     & +(-1.493190e-06)*xx*yy*zz   
	htotal = 10**(-pH_fit)
C       Calculate [CO2*] as defined in DOE Methods Handbook 1994 Ver.2, 
C       ORNL/CDIAC-74, Dickson and Goyet, eds. (Ch 2 p 10, Eq A.49)
	htotal2=htotal*htotal
	co2star=dic*htotal2/(htotal2 + k1*htotal + k1*k2) ! mol/kg
        atmpres=slp/1013.25 ! mbar to atm
        dco2(i)= ( pco2_air*permeg*ff*atmpres - co2star )/permil *1000.0  ! mmol/m^3
        pco2_water(i) = co2star/ff/permeg/atmpres                  ! ppm or uatm
      enddo
      end subroutine eqconc_co2



      subroutine init_bgc_diag
c-----------------------------------------------------------------------
c     initialize NetCDF snapshot file
c-----------------------------------------------------------------------
      use bgc_module
      implicit none
#include "netcdf.inc"
      integer :: ncid,iret,i,j,k,n
      integer :: lon_tdim,lon_udim,z_tdim,z_udim,itimedim
      integer :: lat_tdim,lat_udim, id
      integer :: dims(4), corner(4), edges(4)
      character :: name*24, unit*16

      call def_grid_cdf('bgc.cdf')
      iret=nf_open('bgc.cdf',NF_WRITE,ncid)
      iret=nf_set_fill(ncid, NF_NOFILL, iret)
      call ncredf(ncid, iret)
      iret=nf_inq_dimid(ncid,'xt',lon_tdim)
      iret=nf_inq_dimid(ncid,'yt',lat_tdim)
      iret=nf_inq_dimid(ncid,'zt',z_tdim)
      iret=nf_inq_dimid(ncid,'Time',itimedim)
      dims = (/Lon_tdim,lat_tdim, iTimedim,0/)
      id = ncvdef (ncid,'o2_sflx', NCFLOAT,3,dims,iret)
      name = 'O2 surface flux'; unit = 'mmol/m^2/s'
      call dvcdf(ncid,id,name,24,unit,16,spval)
      id = ncvdef (ncid,'o2_sat', NCFLOAT,3,dims,iret)
      name = 'O2 saturation conc.'; unit = 'mmol/m^3'
      call dvcdf(ncid,id,name,24,unit,16,spval)
      id = ncvdef (ncid,'dic_sflx', NCFLOAT,3,dims,iret)
      name = 'CO2 surface flux'; unit = 'mmol/m^2/s'
      call dvcdf(ncid,id,name,24,unit,16,spval)
      id = ncvdef (ncid,'pco2', NCFLOAT,3,dims,iret)
      name = 'pCO2'; unit = 'ppmv'
      call dvcdf(ncid,id,name,24,unit,16,spval)
      call ncclos (ncid, iret)
      end subroutine init_bgc_diag


      subroutine diag_bgc
c-----------------------------------------------------------------------
c     write to NetCDF snapshot file
c-----------------------------------------------------------------------
      use bgc_module
      implicit none
#include "netcdf.inc"
      integer :: ncid,iret,n,npe, corner(4), edges(4)
      real :: a(imt,js_pe:je_pe),fxa
      integer :: itdimid,ilen,rid,itimeid,fnid,id
      integer :: i,j,is,ie,js,je
      character :: name*24, unit*16
      type(time_type) :: time

      js=max(2,js_pe); je = min(je_pe,jmt-1)
      do npe=0,n_pes
       if (my_pe==npe) then
        iret=nf_open('bgc.cdf',NF_WRITE,ncid)
        iret=nf_set_fill(ncid, NF_NOFILL, iret)
        iret=nf_inq_dimid(ncid,'Time',itdimid)
        iret=nf_inq_dimlen(ncid, itdimid,ilen)
        iret=nf_inq_varid(ncid,'Time',itimeid)
        if (my_pe==0) then
         ilen=ilen+1
         time = current_time-initial_time
         fxa = time%days + time%seconds/86400.
         iret= nf_put_vara_double(ncid,itimeid,ilen,1,fxa)
        endif
        Corner = (/1,js_pe,ilen,1/); edges  = (/imt,j_blk,1,1/)
        iret=nf_inq_varid(ncid,'o2_sflx',id)
        a=o2_sflx(:,js_pe:je_pe)
        where( maskT(:,js_pe:je_pe,km-1) == 0.) a = spval
        iret= nf_put_vara_double(ncid,id,corner,edges,a)

        iret=nf_inq_varid(ncid,'dic_sflx',id)
        a=dic_sflx(:,js_pe:je_pe)
        where( maskT(:,js_pe:je_pe,km-1) == 0.) a = spval
        iret= nf_put_vara_double(ncid,id,corner,edges,a)

        iret=nf_inq_varid(ncid,'o2_sat',id)
        a=o2_sat(:,js_pe:je_pe)
        where( maskT(:,js_pe:je_pe,km-1) == 0.) a = spval
        iret= nf_put_vara_double(ncid,id,corner,edges,a)

        iret=nf_inq_varid(ncid,'pco2',id)
        a=pco2(:,js_pe:je_pe)
        where( maskT(:,js_pe:je_pe,km-1) == 0.) a = spval
        iret= nf_put_vara_double(ncid,id,corner,edges,a)

        call ncclos (ncid, iret)
       endif
       call barrier
      enddo
      end subroutine diag_bgc




