nitrate.F90 6.03 KB
Newer Older
dumoda01's avatar
dumoda01 committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66
!$Id: nitrate.F90,v 1.11 2007-01-06 11:49:15 kbk Exp $
#include"cppdefs.h"
!-----------------------------------------------------------------------
!BOP
!
! !ROUTINE: The nitrate equation \label{sec:nitrate}
!
! !INTERFACE:
   subroutine nitrate(nlev,dt,cnpar,nus,gams,cc)
!
! !DESCRIPTION:
! This subroutine computes the balance of nitrate in the form
!  \begin{equation}
!   \label{SEq}
!    \dot{S}
!    = {\cal D}_S
!    - \frac{1}{\tau^S_R}(S-S_{obs})
!    \comma
!  \end{equation}
!  where $\dot{S}$ denotes the material derivative of the nitrate $N$, and
!  ${\cal D}_N$ is the sum of the turbulent and viscous transport
!  terms modelled according to
!  \begin{equation}
!   \label{DN}
!    {\cal D}_N
!    = \frstder{z}
!     \left(
!        \left( \nu^N_t + \nu^N \right) \partder{N}{z} - \tilde{\Gamma}_N
!        \right)
!    \point
!  \end{equation}
!  In this equation, $\nu^N_t$ and $\nu^N$ are the turbulent and
!  molecular diffusivities of nitrate, respectively,
!  and $\tilde{\Gamma}_N$
!  denotes the non-local flux of nitrate, see
!  \sect{sec:turbulenceIntro}. In the current version of GOTM,
!  we set $\nu^N_t = \nu^\N$ for simplicity.
!
!  Horizontal advection is optionally
!  included  (see {\tt obs.nml}) by means of prescribed
!  horizontal gradients $\partial_xN$ and $\partial_yN$ and
!  calculated horizontal mean velocities $U$ and $V$.
!  Relaxation with the time scale $\tau^N_R$
!  towards a precribed (changing in time)
!  profile $N_{obs}$ is possible.

!  Inner sources or sinks are not considered.
!  The surface freshwater flux is given by means of the precipitation
!  - evaporation data read in as $P-E$ through the {\tt airsea.nml} namelist:
!  \begin{equation}
!     \label{S_sbc}
!    {\cal D}_S =  S (P-E),
!    \qquad \mbox{at } z=\zeta,
!  \end{equation}
!  with $P-E$ given as a velocity (note that ${\cal D}_S$ is the flux in the
!  direction of $z$, and thus positive for a \emph{loss} of nitrate) .
!  Diffusion is numerically treated implicitly,
!  see equations (\ref{sigmafirst})-(\ref{sigmalast}).
!  The tri-diagonal matrix is solved then by a simplified Gauss elimination.
!  Vertical advection is included, and it must be non-conservative,
!  which is ensured by setting the local variable {\tt adv\_mode=0},
!  see section \ref{sec:advectionMean} on page \pageref{sec:advectionMean}.
!
! !USES:

   use meanflow,     only: avmoln                       !CHG3
67
   use meanflow,     only: h,u,v,w,nit,avh              !CHG3
dumoda01's avatar
dumoda01 committed
68 69 70 71 72 73
   use observations, only: dndx,dndy,n_adv              !CHG3
   use observations, only: w_adv_discr,w_adv_method
   use observations, only: nprof,NRelaxTau              !CHG3
   use airsea,       only: p_e
   use util,         only: Dirichlet,Neumann
   use util,         only: oneSided,zeroDivergence
74
!   use bio_var,      only: bio_model
dumoda01's avatar
dumoda01 committed
75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97

   IMPLICIT NONE
!

! !INPUT PARAMETERS:


!  number of vertical layers
   integer, intent(in)                 :: nlev

!   time step (s)
   REALTYPE, intent(in)                :: dt

!  numerical "implicitness" parameter
   REALTYPE, intent(in)                :: cnpar

!  diffusivity of nitrate (m^2/s) (!CHG3 same as salinity)
   REALTYPE, intent(in)                :: nus(0:nlev)

!  non-local salinity flux (psu m/s)
   REALTYPE, intent(in)                :: gams(0:nlev)

!  nitrate concentration after bio loop
98 99 100
!DD Not independent of the bio model, only works with 7 compartments model (Fasham)
!   REALTYPE, intent(in), optional      :: cc(1:10,0:nlev)   !CHG3
   REALTYPE, intent(in), optional      :: cc(1:9,0:nlev)   !CHG3
dumoda01's avatar
dumoda01 committed
101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159
!
! !REVISION HISTORY:
!  Original author(s): Dany Dumont (dany_dumont@ete.inrs.ca)
!
!  $Log: nitrate.F90,v $
!  Creation 1.0  2008-06-11 14:27:00  dd
!  based on salinity.F90
!
!EOP
!
! !LOCAL VARIABLES:
   integer                   :: adv_mode=0
   integer                   :: posconc=1
   integer                   :: i
   integer                   :: DiffBcup,DiffBcdw
   integer                   :: AdvBcup,AdvBcdw
   REALTYPE                  :: DiffSup,DiffSdw
   REALTYPE                  :: AdvSup,AdvSdw
   REALTYPE                  :: Lsour(0:nlev)
   REALTYPE                  :: Qsour(0:nlev)
!
!-----------------------------------------------------------------------
!BOC
!
!  set boundary conditions
   DiffBcup       = Neumann
   DiffBcdw       = Neumann
   DiffSup        = -1.*nit(nlev)*p_e     !CHG3
   DiffSdw        = _ZERO_

   AdvBcup       = zeroDivergence
   AdvBcdw       = oneSided
   AdvSup        = _ZERO_
   AdvSdw        = _ZERO_

!  compute total diffusivity
   do i=0,nlev
      avh(i)=nus(i)+avmoln            !CHG3
   end do

!  add contributions to source term
   Lsour=_ZERO_
   Qsour=_ZERO_

   do i=1,nlev
!     from non-local turbulence
#ifdef NONLOCAL
      Qsour(i) = Qsour(i) - ( gams(i) - gams(i-1) )/h(i)
#endif
   end do

!  ... and from lateral advection
   if (n_adv) then
      do i=1,nlev
         Qsour(i) = Qsour(i) - u(i)*dndx(i) - v(i)*dndy(i)        !CHG3
      end do
   end if

! redefinir nit apres un cyle bio
160
#ifdef BIO
161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181
do i=1,nlev
   !nit(i) = cc(7,i)    ! bio_model=2,4
   nit(i) =  cc(1,i)   ! bio_model=1,6
end do
!   if (bio_model.eq.1) then
!      do i=1,nlev
!         nit(i) = cc(1,i)
!      end do
!   else if (bio_model.eq.2) then
!      do i=1,nlev
!         nit(i) = cc(7,i)
!      end do
!   else if (bio_model.eq.4) then
!      do i=1,nlev
!         nit(i) = cc(7,i)
!      end do
!  else if (bio_model.eq.6) then
!      do i=1,nlev
!         nit(i) = cc(1,i)
!      end do
!   end if
182
#endif
dumoda01's avatar
dumoda01 committed
183 184 185 186

!  do advection step
   if (w_adv_method .ne. 0) then
      call adv_center(nlev,dt,h,h,w,AdvBcup,AdvBcdw,                    &
187
                          AdvSup,AdvSdw,w_adv_discr,adv_mode,nit) !CHG3
dumoda01's avatar
dumoda01 committed
188 189 190 191
   end if

!  do diffusion step
   call diff_center(nlev,dt,cnpar,posconc,h,DiffBcup,DiffBcdw,          &
192
                    DiffSup,DiffSdw,avh,LSour,Qsour,NRelaxTau,nprof,nit)     !CHG3
dumoda01's avatar
dumoda01 committed
193 194 195 196 197 198 199 200

   return
   end subroutine nitrate
!EOC

!-----------------------------------------------------------------------
! Copyright by the GOTM-team under the GNU Public License - www.gnu.org
!-----------------------------------------------------------------------