ROMS_r215: /home/cholod/wd/peps/ROMS/romsagrif_r215/Roms_tools/Roms_Agrif/bio

00001 !
00002 ! $Id$
00003 !
00004 #include "cppdefs.h"
00005 #if defined BIOLOGY && defined BIO_NPZD
00006 
00007       subroutine biology_tile (Istr,Iend,Jstr,Jend)
00008 !
00009 ! Compute biological forcing functions
00010 !
00011 ! In this particular implementation there is 4 compartments:
00012 ! NO3, PHYTOplankton, ZOOplanknton, DETritus.
00013 !
00014       implicit none
00015       integer Istr,Iend,Jstr,Jend
00016 # include "param.h"
00017 # include "grid.h"
00018 # include "ocean3d.h"
00019 # include "ocean2d.h"
00020 # include "diagnostics.h"
00021 # include "scalars.h"
00022 # include "forces.h"
00023 # include "mixing.h"
00024       real kwater, palpha, kChla, CN_Phyt, theta_m, opc,
00025      &     K_NO3, mu_P_D, gmax, beta, K_Phyt,
00026      &     mu_Z_A, mu_Z_D, mu_D_N,
00027      &     wPhyt, wDet
00028 # ifdef OXYGEN
00029      &   , CN_Z
00030 # endif
00031       integer ITERMAX
00032       integer nsink
00033       parameter (
00034      &  ITERMAX = 3,      ! number of small implicit time step
00035      &  nsink   = NumVSinkTerms + 1, ! add'lly: Chlorophyll
00036 !
00037 ! Parameters as in Table 1; Fasham et al. [JMR, 48, 591-639, 1990]
00038 !
00039      &  kwater  = 0.04,    ! light attenuation due to sea water  [m-1]
00040                            ! range:(0.04<==>0.04];    units:[m-1]
00041      &  palpha  = 1.0,     ! initial slope of the P-I curve
00042                            ! range:(1.00<==>1.00);       [(W m-2 d)-1]
00043      &  kChla  = 0.024,    ! light attenuation by Chlorophyl  
00044                            !                         [(m^2 mg Chla)-1]
00045      &  CN_Phyt= 6.625,    ! C:N ratio for phytoplankton
00046                            !                       [mMol C (mMol N)-1]
00047 # ifdef OXYGEN
00048      &  CN_Z   = 6.625,    ! C:N ratio for zoo
00049                            ! range:(4.<==>6.);     [mol-C (mol-N)-1]
00050 # endif
00051      &  theta_m= 0.0535,   ! max Cellular Chlorophyll to Carbon Ratio
00052                            ! range:(0.015<==>0.072);    [mg Chla/mg C]
00053      &  K_NO3   = 1./0.5,  ! inverse half-saturation for Phytoplankton
00054                            ! range:(1./.0 <==> 1./.9);[1/(mmol-N m-3)]
00055      &  mu_P_D  = 0.03,    ! Phyto mortality to Det rate        [d-1]
00056      &  gmax    = 0.9,     ! maximum Zooplankton growth rate     [d-1]
00057      &  beta    = 0.75,    ! Zooplankton assimilation efficiency of 
00058                            !                       Phytoplankton [n.d.]
00059      &  K_Phyt  = 1.0,     ! Zooplankton half-saturation constant
00060                            ! for ingestion of phyto [d-1]
00061      &  mu_Z_A  = 0.10,    ! Zooplankton specific excretion rate [d-1]
00062      &  mu_Z_D  = 0.10,    ! Zooplankton mortality to Detritus   [d-1]
00063      &  mu_D_N  = 0.05,    ! Detrital remineralization to NO3 rate
00064                            ! range:( <==>   );        units:[d-1]
00065      &  wPhyt   = 0.5,     ! sinking velocities for Phytoplankton [m.d-1]
00066      &  wDet    = 5.0  )   !    ''       ''      '' Detritus
00067 !
00068       integer i,j,k, ITER, iB
00069       real    NO3(N), Phyt(N), Zoo(N), Det(N), Chla(N), 
00070      &        aJ(N),FC(0:N),
00071      &        PAR, PARsup, attn, Vp, Epp, cu, aL,aR, dtdays, L_NO3,
00072      &        E_NO3,cff,cff0,cff1,cff2,cff6,
00073      &        SB(N,nsink),dSB(0:N,nsink),wSB(nsink)
00074 # ifdef OXYGEN
00075      &      , O2(N), tem(N), sal(N), den(N)
00076      &      , O2satu_loc, Kv_O2_loc
00077      &      , eos80, u10_loc, Sc
00078 #  ifdef OCMIP_OXYGENSAT
00079      &       , o2sato   ! OCMIP function, calculates O2 saturation
00080 #  else /* OCMIP_OXYGENSAT */
00081      &       , satpc    ! oxygen saturation in % (calculated, but unused)
00082      &       , AOU      ! Apparent oxygen utilization (calc., but unused)
00083 #  endif /* OCMIP_OXYGENSAT */
00084 # endif /* OXYGEN */
00085 # if defined OXYGEN || defined DIAGNOSTICS_BIO
00086      &      , dtsec     ! length of time step in seconds (for gas exchange)
00087 # endif
00088 # ifdef DIAGNOSTICS_BIO
00089       integer l, iflux
00090       real    trend_no3,trend_phy,trend_zoo,trend_det,somme
00091      &      , bilan_no3,bilan_phy,bilan_zoo,bilan_det
00092      &      , ThisVSinkFlux(N, NumVSinkTerms)  ! [mmol m-2 s-1], upward flux is positive
00093      &      , ThisFlux(N, NumFluxTerms)
00094 #  ifdef OXYGEN
00095      &      , ThisGasExcFlux(NumGasExcTerms)
00096 #  endif
00097      &      , LastVSinkFlux,ColumnMassOld(NumVSinkTerms)
00098      &      , ColumnMassNew(NumVSinkTerms)
00099 # endif /* DIAGNOSTICS_BIO */
00100 !
00101 # include "compute_auxiliary_bounds.h"
00102 !
00103 # undef DEBUG
00104 
00105       dtdays=dt/(24.*3600.*float(ITERMAX))  ! time step as fraction of day.
00106 # if defined DIAGNOSTICS_BIO || defined OXYGEN
00107       dtsec = dt / float(ITERMAX)           ! time step in seconds
00108 # endif /* DIAGNOSTICS_BIO || OXYGEN */
00109 !
00110 !
00111 ! Since the following solver is iterative to achieve implicit
00112 ! discretization of the biological interaction, two time slices are
00113 ! required, BIO where BIO is understood as vector of
00114 ! biological state variables: BIO=[NO3,Phyt,Zoo,Det]. Assume
00115 ! that the iterations converge, the newly obtained state variables
00116 ! satisfy equations
00117 !
00118 !           BIO = BIO + dtdays * rhs(BIO)
00119 ! 
00120 ! where rhs(BIO) is the vector of biological r.h.s. computed at
00121 ! the new time step. During the iterative procedure a series of
00122 ! fractional time steps is performed in a chained mode (splitting
00123 ! by different biological conversion processes) in sequence NO3 -- 
00124 ! Phyt -- Zoo -- Det, that is the main food chain. In all 
00125 ! stages the concentration of the component being consumed is
00126 ! treated in fully implicit manner, so that the algorithm guarantees
00127 ! non-negative values, no matter how strong is the concentration of
00128 ! active consuming component (Phyto or Zoo).
00129 !
00130 ! The overall algorithm, as well as any stage of it is formulated
00131 ! in conservative form (except explicit sinking) in sense that the
00132 ! sum of concentration of all five components is conserved.
00133 !
00134 !# ifdef EW_PERIODIC
00135 !#  define I_RANGE Istr,Iend
00136 !# else
00137 !#  define I_RANGE IstrR,IendR
00138 !# endif
00139 !# ifdef NS_PERIODIC
00140 !#  define J_RANGE Jstr,Jend
00141 !# else
00142 !#  define J_RANGE JstrR,JendR
00143 !# endif
00144 # define I_RANGE Istr,Iend
00145 # define J_RANGE Jstr,Jend
00146 
00147       do j=J_RANGE
00148         do i=I_RANGE
00149 # ifdef DIAGNOSTICS_BIO
00150 ! Reset the biogeochemical fluxes. This is necessary because the
00151 ! biological routine uses multiple. time steps for each physical time
00152 ! step.
00153           do k=1,N
00154             do l=1,NumFluxTerms
00155               bioFlux(i,j,k,l) = 0.0
00156             enddo
00157           enddo
00158           do k=0,N
00159             do l=1,NumVSinkTerms
00160               bioVSink(i,j,k,l) = 0.0
00161             enddo
00162           enddo
00163 #  ifdef OXYGEN
00164           do l=1,NumGasExcTerms
00165             GasExcFlux(i,j,l) = 0.0
00166           enddo
00167 #  endif
00168 # endif /* DIAGNOSTICS_BIO */
00169 !
00170 ! Extract biological variables from tracer arrays; place them into
00171 ! scratch variables; restrict their values to be positive definite.
00172 !
00173           do k=1,N
00174             NO3(k) =max(t(i,j,k,nnew,iNO3_)  ,0.)   ! Nitrate
00175             Phyt(k)=max(t(i,j,k,nnew,iPhy1)  ,0.)   ! Phytoplankton
00176             Chla(k)=max(t(i,j,k,nnew,iChla)  ,0.)   ! Chlor a
00177             Zoo(k) =max(t(i,j,k,nnew,iZoo1)  ,0.)   ! Zooplankton
00178             Det(k) =max(t(i,j,k,nnew,iDet1)  ,0.)   ! Detritus
00179 !
00180             if (Phyt(k) .gt. 0.001 .and. Chla(k) .gt. 0.001) then
00181              theta(i,j,k) = Chla(k)/(Phyt(k)*CN_Phyt*12.)      ! Chla/Phyt ratio
00182              if (theta(i,j,k).gt.theta_m) theta(i,j,k)=theta_m ! [mg Chla(mg C)-1
00183             else                                     ! [mg Chla (mg C)-1]
00184              theta(i,j,k) = theta_m
00185             endif
00186 # ifdef OXYGEN
00187             tem(k) =max(t(i,j,k,nnew,itemp),0.)      ! temperature; [deg. C]
00188             sal(k) =max(t(i,j,k,nnew,isalt),0.)      ! salinity; [PSU]
00189 #  ifndef OCMIP_OXYGENSAT
00190             den(k) =1000.+eos80(0.,tem(k),sal(k))    ! potential density; [kg m-3]
00191 #  endif
00192             O2(k)  =max(t(i,j,k,nnew,iO2),0.)        ! Oxygen;  [mmol O2 m-3]
00193 # endif 
00194           enddo
00195 
00196 
00197           DO ITER=1,ITERMAX      !--> Start internal iterations to achieve
00198                              !    nonlinear backward-implicit solution.
00199 
00200             PAR=srflx(i,j)*rho0*Cp*0.43
00201             opc=0.01*PAR
00202 
00203             if (PAR.gt.0.) then
00204 !
00205 !   *** SUN IS UP ***
00206 !
00207 ! Calulate aJ: Set Photosynthetically Available Radiation (PAR) at
00208 ! surface from solar radiation x 0.43. Then, within each grid box
00209 ! compute attenuation coefficient based on the concentration of
00210 ! Phytoplankton inside the grid box, and attenuate PAR from surface
00211 ! down (thus, PAR at certain depth depends on the whole distribution
00212 ! of Phytoplankton above). To compute aJ, one needs PAR somewhat in
00213 ! the middle of the gridbox, so that attenuation "attn" corresponds
00214 ! to half of the grid box height, while PAR is multiplied by it
00215 ! twice: once to get it in the middle of grid-box and once the
00216 ! compute on trhe lower grid-box interface;
00217 !
00218               do k=N,1,-1     !<-- irreversible
00219 
00220                 attn=exp(-0.5*(kwater+kChla*Chla(k))*
00221      &                 (z_w(i,j,k)-z_w(i,j,k-1)))
00222 
00223                 PARsup=PAR*attn
00224                 Vp=0.59*(1.066**t(i,j,k,nnew,itemp))   ! From Eppley
00225 !                Vp=0.90*(1.066**t(i,j,k,nnew,itemp))  ! a.b^cT=µ_max=3.0
00226                 cff0=PARsup*palpha*theta(i,j,k)        ! for diatoms & 
00227                 Epp=Vp/sqrt(Vp*Vp+cff0*cff0)           ! 2.0 for flagelates
00228                 aJ(k)=Epp*cff0
00229 !
00230 ! theta adaptation
00231                 L_NO3=K_NO3*NO3(k)/(1+K_NO3*NO3(k))
00232                 cff=dtdays*Epp*cff0*L_NO3
00233                 theta(i,j,k)=(theta(i,j,k)+theta_m*Epp*cff*L_NO3)
00234      &                                                  /(1.+cff)
00235 ! (1) NO3 uptake by Phyto
00236 !
00237                 E_NO3=K_NO3/(1+K_NO3*NO3(k)) ! Parker 1993 Ecol Mod. 66 113-120
00238                 cff=dtdays*aJ(k)*Phyt(k)*E_NO3
00239                 NO3(k)=NO3(k)/(1.+cff)
00240 # ifdef DIAGNOSTICS_BIO
00241                 ThisFlux(k, NFlux_NewProd) = cff*NO3(k)
00242 !
00243 #  ifdef OXYGEN
00244 ! production of O2 by phyto growth
00245                 ThisFlux(k, OGain_NewProd) =
00246      &             ThisFlux(k, NFlux_NewProd) * (CN_Phyt + 2.)
00247 #  endif /* OXYGEN */
00248 # endif /* DIAGNOSTICS_BIO */
00249 # ifdef OXYGEN
00250                 O2(k) = O2(k) + cff*NO3(k)*(CN_Phyt + 2.)
00251 # endif
00252                 Phyt(k)=Phyt(k)+cff*NO3(k)
00253 
00254                 PAR=PARsup*attn
00255 
00256                 if (PARsup.ge.opc) then  !  Compute euphotic depth
00257                   if (PAR.ge.opc) then
00258                     hel(i,j)=-z_w(i,j,k-1)
00259                   else
00260                     hel(i,j)=-z_r(i,j,k)
00261                   endif
00262                 endif
00263 !
00264               enddo
00265 !
00266             else
00267 !
00268 !   *** SUN IS DOWN ***
00269 !
00270 # ifdef DIAGNOSTICS_BIO
00271               do k=N,1,-1
00272                 ThisFlux(k, NFlux_NewProd) = 0.0
00273 #  ifdef OXYGEN
00274                 ThisFlux(k, OGain_NewProd) = 0.0
00275 #  endif /* OXYGEN */
00276               enddo
00277 # endif /* DIAGNOSTICS_BIO */
00278               hel(i,j)=0.0
00279             endif
00280 
00281 !
00282 ! (1) Phytoplankton grazing by Zooplankton to Zoo and Detr
00283 ! (2) Phytoplankton mortality to Detr (mu_P_D)
00284 !
00285             do k=1,N
00286               cff1=dtdays*gmax*Zoo(k)/(K_Phyt+Phyt(k))
00287               cff2=dtdays*mu_P_D
00288               Phyt(k)=Phyt(k)/(1.+cff1+cff2)
00289               Zoo(k)=Zoo(k)+Phyt(k)*cff1*beta
00290 # ifdef DIAGNOSTICS_BIO
00291               ThisFlux(k, NFlux_Grazing)=Phyt(k)*cff1*beta
00292               ThisFlux(k, NFlux_SlopFeed) = Phyt(k) * cff1 * (1.-beta)
00293               ThisFlux(k, NFlux_Pmort) = Phyt(k) * cff2
00294 # endif /* DIAGNOSTICS_BIO */
00295               Det(k)=Det(k)+Phyt(k)*(cff1*(1.-beta)+cff2)
00296 !
00297 ! (1) Zoo excretion to NO3  (rate mu_Z_A)
00298 ! (2) Zoo mortality to Det (rate mu_Z_D)
00299 !
00300               cff1=dtdays*mu_Z_A
00301               cff2=dtdays*mu_Z_D*Zoo(k)
00302               Zoo(k)=Zoo(k)/(1.+cff1+cff2)
00303 # ifdef DIAGNOSTICS_BIO
00304               ThisFlux(k, NFlux_Zmetab)=cff1*Zoo(k)
00305               ThisFlux(k, NFlux_Zmort)=cff2*Zoo(k)
00306 #  ifdef OXYGEN
00307 !         Zoo uptake of O2 (rate t_Zbmet + R_C)
00308 !         there is no control yet for assuring non-negative Oxygen
00309 !         values!
00310               ThisFlux(k, OLoss_Zmetab) =
00311      &           ThisFlux(k, NFlux_Zmetab) * CN_Z
00312 #  endif /* OXYGEN */
00313 # endif /* DIAGNOSTICS_BIO */
00314 # ifdef OXYGEN 
00315               O2(k)  = O2(k) - cff1 * Zoo(k) * CN_Z
00316 # endif
00317               NO3(k)=NO3(k)+Zoo(k)*cff1
00318               Det(k)=Det(k)+Zoo(k)*cff2
00319 !
00320 ! (1) Det remineralization to N03
00321 !
00322               cff1=dtdays*mu_D_N
00323               Det(k)=Det(k)/(1.+cff1)
00324 # ifdef DIAGNOSTICS_BIO
00325               ThisFlux(k, NFlux_ReminD)=Det(k)*cff1
00326 #  ifdef OXYGEN
00327 !         Loss of O2 in Det/Det remineralization
00328               ThisFlux(k, OLoss_ReminD)=
00329      &             ThisFlux(k,NFlux_ReminD) * (CN_Phyt+2.)
00330 #  endif /* OXYGEN */
00331 # endif /* DIAGNOSTICS_BIO */
00332 # ifdef OXYGEN
00333               O2(k) = O2(k) - Det(k) * cff1 * (CN_Phyt+2.)
00334 # endif
00335               NO3(k)=NO3(k)+Det(k)*cff1
00336             enddo
00337 !
00338 # ifdef OXYGEN
00339 #  ifdef OCMIP_OXYGEN_SC
00340 !*********************************************************************
00341 !  alternative formulation (Sc will be slightly smaller up to about 35
00342 !  C)
00343 !  Computes the Schmidt number of oxygen in seawater using the
00344 !  formulation proposed by Keeling et al. (1998, Global Biogeochem.
00345 !  Cycles, 12, 141-163).  Input is temperature in deg C.
00346 !
00347             Sc = 1638.0 - 81.83*Tem(N) + 1.483*(Tem(N)**2) -
00348      &             0.008004*(Tem(N)**3)
00349 !*********************************************************************
00350 #  else /* OCMIP_OXYGEN_SC */
00351 !       calculate the Schmidt number for O2 in sea water [Wanninkhof,
00352 !       1992]
00353             Sc=1953.4 - 128.0*Tem(N) + 3.9918*(Tem(N)**2) -
00354      &             0.050091*(Tem(N)**3)
00355 #  endif /* OCMIP_OXYGEN_SC */
00356 !
00357 !       calculate the wind speed from the surface stress values
00358             u10_loc = sqrt(sqrt( (0.5*(sustr(i,j)+sustr(i+1,j)))**2
00359      &                          +(0.5*(svstr(i,j)+svstr(i,j+1)))**2)
00360      &       * rho0 * 550.)   ! 550 = 1 / (1.3 * 0.0014) (=rho_air * CD)
00361 !       calculate the gas transfer coef for O2
00362             Kv_O2_loc=0.31*u10_loc*u10_loc*sqrt(660./Sc)/(100.*3600.) 
00363 !  denominator: convert Kv from [cm/h] to [m/s]
00364 !       calculate the saturation oxygen level
00365 #  ifdef OCMIP_OXYGENSAT
00366             O2satu_loc = o2sato(Tem(N), Sal(N))
00367 #  else /* OCMIP_OXYGENSAT */
00368             call O2sato(O2(N),Tem(N),Sal(N),den(N),O2satu_loc,satpc,AOU)
00369 #  endif /* OCMIP_OXYGENSAT */
00370 !       air-sea flux of O2
00371 !       abs(z_w(i,j,N-1))==> volume of upper layer
00372 #  ifdef DIAGNOSTICS_BIO
00373             ThisGasExcFlux(OFlux_GasExc)=Kv_O2_loc*(O2satu_loc-O2(N))
00374 !       ThisGasExcFlux is positive if ocean takes up O2 from the
00375 !       atmosphere
00376 #  endif
00377             O2(N) = O2(N) + Kv_O2_loc * (O2satu_loc - O2(N)) * dtsec /
00378      &                  abs(z_w(i,j,N-1))
00379 # endif /* OXYGEN */
00380 
00381 !
00382 ! Vertical sinking: Vertical advection algorithm based on monotonic,
00383 ! continuous conservative parabolic splines.
00384 !
00385             do k=1,N
00386               SB(k,1)=theta(i,j,k)*Phyt(k)*CN_Phyt*12.
00387               SB(k,2)=Phyt(k)
00388               SB(k,3)=Det(k)
00389             enddo
00390             wSB(1)=wPhyt
00391             wSB(2)=wPhyt
00392             wSB(3)=wDet
00393 
00394             do iB=1,nsink
00395 
00396 ! Part (i): Construct parabolic splines: compute vertical derivatives
00397 ! of the fields SB. The derivatives are located at W-points;
00398 ! Neumann boundary conditions are assumed on top and bottom.
00399 !
00400               dSB(0,iB)=0.
00401               FC(0)=0.
00402               cff6=6.
00403               do k=1,N-1
00404                 cff=1./(2.*Hz(i,j,k+1)+Hz(i,j,k)*(2.-FC(k-1)))
00405                 FC(k)=cff*Hz(i,j,k+1)
00406                 dSB(k,iB)=cff*(cff6*(SB(k+1,iB)-SB(k,iB))
00407      &                             -Hz(i,j,k)*dSB(k-1,iB))
00408               enddo
00409               dSB(N,iB)=0.
00410               do k=N-1,1,-1     !<-- irreversible
00411                 dSB(k,iB)=dSB(k,iB)-FC(k)*dSB(k+1,iB)
00412               enddo
00413 !
00414 ! Part (ii): Convert dSB [which are now vertical derivatives
00415 ! of fields SB at the grid box interfaces] into field values
00416 ! at these interfaces, assuming parabolic profiles within each grid
00417 ! box. Restrict these values to lie between bounds determined from
00418 ! box-averaged values of grid boxes adjscent from above and below.
00419 ! (This restriction is part of PPM-like monotonization procedure.)
00420 !
00421               cff=1./3.
00422               dSB(0,iB)=SB(1,iB) !-cff*Hz(1)*(dSB(0,iB)+0.5*dSB(1,iB))
00423               dSB(N,iB)=SB(N,iB) !+cff*Hz(N)*(dSB(N,iB)+0.5*dSB(N-1,iB))
00424               do k=2,N          !<-- irreversible
00425                 dSB(k-1,iB)=SB(k,iB)
00426      &                   -cff*Hz(i,j,k)*(0.5*dSB(k,iB)+dSB(k-1,iB))
00427                 dSB(k-1,iB)=max(dSB(k-1,iB),min(SB(k-1,iB),SB(k,iB)))
00428                 dSB(k-1,iB)=min(dSB(k-1,iB),max(SB(k-1,iB),SB(k,iB)))
00429               enddo
00430 !  
00431 ! Part (iii): Convert dSB into flux-integrated values,
00432 ! complete PPM flux limiting. This procedure starts from assigning
00433 ! Left and Right (aR,aL) values of the interpolating parabolae, then
00434 ! monotonicity conditions are checked and aL,aR are modified to fit.
00435 ! Overall, from this moment and further on it follows Colella--
00436 ! --Woodward, 1984 bombmaking code almost exactly.
00437 !
00438               do k=1,N           !<-- irreversible
00439                 FC(k)=dtdays/Hz(i,j,k)
00440                 aR=dSB(k,iB)
00441                 aL=dSB(k-1,iB)
00442                 cff1=(aR-aL)*6.*(SB(k,iB)-.5*(aR+aL))
00443                 cff2=(aR-aL)**2
00444                 if ((aR-SB(k,iB))*(SB(k,iB)-aL).lt.0.) then
00445                   aL=SB(k,iB)
00446                   aR=SB(k,iB)
00447                 elseif (cff1.gt.cff2) then
00448                   aL=3.*SB(k,iB)-2.*aR
00449                 elseif (cff1.lt.-cff2) then
00450                   aR=3.*SB(k,iB)-2.*aL
00451                 endif
00452                 cu=wSB(iB)*FC(k)
00453                 dSB(k-1,iB)=SB(k,iB)-(1.-cu)*(.5*(aR-aL)-(.5*(aR+aL)
00454      &                                   -SB(k,iB) )*(1.-2.*cu))
00455               enddo
00456               dSB(N,iB)=0.   ! Set no-flux boundary conditions at top.
00457 !
00458 ! Apply fluxes:
00459 !
00460               do k=1,N
00461                 SB(k,iB)=SB(k,iB)+wSB(iB)*FC(k)*(dSB(k,iB)-dSB(k-1,iB))
00462               enddo
00463             enddo  ! <-- iB
00464 
00465 # ifdef DIAGNOSTICS_BIO
00466             do iflux = 1, NumVSinkTerms
00467               ColumnMassOld(iflux) = 0.0
00468               ColumnMassNew(iflux) = 0.0
00469             end do
00470 # endif /* DIAGNOSTICS_BIO */
00471 
00472             do k=1,N
00473               theta(i,j,k)= SB(k,1)/(SB(k,2)*CN_Phyt*12.+1.E-20)
00474               if (theta(i,j,k).gt.theta_m) theta(i,j,k)=theta_m
00475 # ifdef DIAGNOSTICS_BIO
00476 ! ColumnMassOld and ColumnMassNew are needed to compute the sinking flux
00477 ! into the sediment
00478               ColumnMassOld(1)=ColumnMassOld(1)
00479      &                    +Phyt(k)*Hz(i,j,k)/(pn(i,j)*pm(i,j))
00480               ThisVSinkFlux(k, NFlux_VSinkP1)=Phyt(k)-SB(k,2)
00481 # endif /* DIAGNOSTICS_BIO */
00482               Phyt(k) = SB(k,2)
00483 # ifdef DIAGNOSTICS_BIO
00484               ColumnMassNew(1)=ColumnMassNew(1)
00485      &                    +Phyt(k)*Hz(i,j,k)/(pn(i,j)*pm(i,j))
00486 # endif /* DIAGNOSTICS_BIO */
00487           ! detritus
00488 # ifdef DIAGNOSTICS_BIO
00489               ColumnMassOld(2)=ColumnMassOld(2)
00490      &                    +Det(k)*Hz(i,j,k)/(pn(i,j)*pm(i,j))
00491               ThisVSinkFlux(k, NFlux_VSinkD1)=Det(k)-SB(k,3)
00492 # endif /* DIAGNOSTICS_BIO */
00493               Det(k) = SB(k,3)
00494 # ifdef DIAGNOSTICS_BIO
00495               ColumnMassNew(2)=ColumnMassNew(2)
00496      &                    +Det(k)*Hz(i,j,k)/(pn(i,j)*pm(i,j))
00497 # endif /* DIAGNOSTICS_BIO */
00498             enddo
00499 !
00500 # ifdef DIAGNOSTICS_BIO
00501 ! Transfer fluxes to global arrays at the end of each biological time step
00502 ! for computational efficiency, divide now by dtsec to get the correct units
00503             do iflux = 1, NumFluxTerms
00504               do k = 1, N
00505                 bioFlux(i,j,k,iflux) = ( bioFlux(i,j,k,iflux) +
00506 !     &         ThisFlux(k, iflux)*Hz(i,j,k)/(pn(i,j)*pm(i,j)*dt) )
00507      &           ThisFlux(k, iflux)/dt )
00508 #  ifdef MASKING
00509      &             * rmask(i,j)
00510 #  endif /* MASKING */
00511               enddo
00512             enddo
00513             do iflux = 1, NumVSinkTerms
00514 ! Compute the vertical sinking flux into the sediment by comparing
00515 ! previous and current mass in this (i,j) column
00516 ! The flux is positive if upward, so usually it will be
00517 ! negative, i.e. into the sediment.
00518               LastVSinkFlux = ( ColumnMassNew(iflux) -
00519      &                        ColumnMassOld(iflux) ) / dtsec
00520               bioVSink(i,j,0,iflux) = (  bioVSink(i,j,0,iflux)
00521      &                             + LastVSinkFlux / float(ITERMAX) )
00522 #  ifdef MASKING
00523      &          * rmask(i,j)
00524 #  endif /* MASKING */
00525               do k = 1, N
00526                 LastVSinkFlux = LastVSinkFlux +
00527      &               ( ThisVSinkFlux(k,iflux) 
00528 !     &            *Hz(i,j,k) / (pn(i,j)*pm(i,j)*dtsec) )
00529      &                                 /dtsec )
00530                 bioVSink(i,j,k,iflux) =(  bioVSink(i,j,k,iflux)
00531      &                              + LastVSinkFlux / float(ITERMAX) )
00532 #  ifdef MASKING
00533      &                                   * rmask(i,j)
00534 #  endif /* MASKING */
00535               end do
00536             end do
00537 #  ifdef OXYGEN
00538             do iflux = 1, NumGasExcTerms
00539                GasExcFlux(i,j,nnew,iflux)=( GasExcFlux(i,j,nnew,iflux)
00540      &                      + ThisGasExcFlux(iflux) / float(ITERMAX) )
00541      &                                 *pn(i,j)*pm(i,j)/Hz(i,j,k)  !PM
00542 #   ifdef MASKING
00543      &                                   * rmask(i,j)
00544 #   endif /* MASKING */
00545             end do
00546 #  endif
00547 # endif /* DIAGNOSTICS_BIO */
00548 !
00549           ENDDO  ! <-- ITER
00550 !
00551 # if defined DIAGNOSTICS_BIO && defined DEBUG
00552           do k=1,N
00553             do iflux=1,NumFluxTerms
00554               bioFlux(i,j,k,iflux) = bioFlux(i,j,k,iflux)
00555 !     &                    *pn(i,j)*pm(i,j)*dt / Hz(i,j,k)
00556             enddo
00557             do iflux=1,NumVSinkTerms
00558               bioVSink(i,j,k,iflux) = bioVSink(i,j,k,iflux)
00559 !     &                    *pn(i,j)*pm(i,j)*dtsec / Hz(i,j,k)
00560             enddo
00561           enddo
00562           if ((i.eq.13) .and. (j.eq.15)) then
00563             bilan_no3=bioFlux(i,j,N,5)+bioFlux(i,j,N,7)
00564      &               -bioFlux(i,j,N,1)
00565             bilan_phy=bioFlux(i,j,N,1)-bioFlux(i,j,N,4)
00566      &               -bioFlux(i,j,N,2)-bioFlux(i,j,N,3)
00567             bilan_zoo=bioFlux(i,j,N,2)-bioFlux(i,j,N,5)
00568      &               -bioFlux(i,j,N,6)
00569             bilan_det=bioFlux(i,j,N,3)+bioFlux(i,j,N,6)
00570      &               +bioFlux(i,j,N,4)-bioFlux(i,j,N,7)
00571             somme=bilan_no3+bilan_phy+bilan_zoo+bilan_det
00572             trend_no3= ( (min(t(i,j,N,nnew,iNO3_),0.) +NO3(N))
00573      &                   - t(i,j,N,nnew,iNO3_) )
00574             trend_phy= ( (min(t(i,j,N,nnew,iPhy1),0.) +Phyt(N))
00575      &                   - t(i,j,N,nnew,iPhy1) )
00576             trend_zoo= ( (min(t(i,j,N,nnew,iZoo1),0.) +Zoo(N))
00577      &                   - t(i,j,N,nnew,iZoo1) )
00578             trend_det= ( (min(t(i,j,N,nnew,iDet1),0.) +Det(N))
00579      &                   - t(i,j,N,nnew,iDet1) )
00580             print*, 'balance = ',somme
00581             print*, 'bilan_no3 - trend_no3 = ',bilan_no3-trend_no3
00582             print*, 'bilan_phy - trend_phy = ',bilan_phy-trend_phy
00583             print*, 'bilan_zoo - trend_zoo = ',bilan_zoo-trend_zoo
00584             print*, 'bilan_det - trend_det = ',bilan_det-trend_det
00585             print*, 'bioFlux  1 = ',bioFlux(i,j,N,1)
00586             print*, 'bioFlux  2 = ',bioFlux(i,j,N,2)
00587             print*, 'bioFlux  3 = ',bioFlux(i,j,N,3)
00588             print*, 'bioFlux  4 = ',bioFlux(i,j,N,4)
00589             print*, 'bioFlux  5 = ',bioFlux(i,j,N,5)
00590             print*, 'bioFlux  6 = ',bioFlux(i,j,N,6)
00591             print*, 'bioFlux  7 = ',bioFlux(i,j,N,7)
00592             print*, 'bioVSink 1 = ',bioVSink(i,j,N,1)
00593             print*, 'bioVSink 2 = ',bioVSink(i,j,N,2)
00594           endif
00595 # endif /* DIAGNOSTICS_BIO && DEBUG */
00596 !
00597           do k=1,N
00598             t(i,j,k,nnew,iNO3_)=min(t(i,j,k,nnew,iNO3_),0.) +NO3(k)
00599             t(i,j,k,nnew,iPhy1)=min(t(i,j,k,nnew,iPhy1),0.) +Phyt(k)
00600             t(i,j,k,nnew,iZoo1)=min(t(i,j,k,nnew,iZoo1),0.) +Zoo(k)
00601             t(i,j,k,nnew,iDet1)=min(t(i,j,k,nnew,iDet1),0.) +Det(k)
00602             t(i,j,k,nnew,iChla)=min(t(i,j,k,nnew,iChla),0.) +
00603      &                            CN_Phyt*12.*Phyt(k)*theta(i,j,k)
00604 # ifdef OXYGEN
00605             t(i,j,k,nnew,iO2)  =min(t(i,j,k,nnew,iO2),0.)   +O2(k)
00606 # endif
00607           enddo
00608 # ifdef OXYGEN
00609           O2satu(i,j) = O2satu_loc
00610           Kv_O2(i,j) = Kv_O2_loc
00611           u10(i,j) = u10_loc
00612 # endif /* OXYGEN */
00613 
00614         enddo
00615       enddo       !<-- i,j 
00616 
00617 # undef DEBUG
00618 #else
00619       subroutine biology_empty ()
00620 #endif
00621       return
00622       end
00623 !
00624 !====================================================================
00625 !
00626       function o2sato(T,S)
00627 !
00628 ! ********************************************************************
00629 ! Computes the oxygen saturation concentration at 1 atm total pressure
00630 ! in mol/m^3 given the temperature (t, in deg C) and the salinity (s,
00631 ! in permil). 
00632 !
00633 ! FROM GARCIA AND GORDON (1992), LIMNOLOGY and OCEANOGRAPHY.
00634 ! THE FORMULA USED IS FROM PAGE 1310, EQUATION (8).
00635 !
00636 ! *** NOTE: THE "A3*TS^2" TERM (IN THE PAPER) IS INCORRECT. ***
00637 ! *** IT SHOULDN'T BE THERE.                                ***
00638 !
00639 ! o2sato IS DEFINED BETWEEN T(freezing) <= T <= 40(deg C) AND
00640 ! 0 permil <= S <= 42 permil
00641 ! C
00642 ! CHECK VALUE:  T = 10.0 deg C, S = 35.0 permil, 
00643 ! o2sato = 0.282015 mol/m^3
00644 !
00645 ! The unit for the return value is now [mmol/m^3], as needed by ROMS
00646 !
00647 ! H. Frenzel, UCLA 1999
00648 ! ********************************************************************
00649 !
00650       implicit none
00651       real    o2sato
00652 
00653 ! input variables
00654       real   T   ! Temperature [C]
00655       real   S   ! Salinity [PSU]
00656 
00657       real A0, A1, A2, A3, A4, A5, B0, B1, B2, B3, C0
00658       parameter (A0 = 2.00907, A1 = 3.22014, A2 = 4.05010,
00659      &           A3 = 4.94457, A4 = -2.56847E-1, A5 = 3.88767,
00660      &           B0=-6.24523E-3, B1=-7.37614E-3, B2=-1.03410E-2,
00661      &           B3=-8.17083E-3, C0=-4.88682E-7)
00662 
00663 ! other variables
00664       real    TT      
00665       real    TK        ! Temperature [K]
00666       real    TS, TS2, TS3, TS4, TS5
00667       real    CO
00668 !
00669       TT  = 298.15-T
00670       TK  = 273.15+T
00671       TS  = LOG(TT/TK)
00672       TS2 = TS**2
00673       TS3 = TS**3
00674       TS4 = TS**4
00675       TS5 = TS**5
00676       CO  = A0 + A1*TS + A2*TS2 + A3*TS3 + A4*TS4 + A5*TS5
00677      &     + S*(B0 + B1*TS + B2*TS2 + B3*TS3)
00678      &     + C0*(S*S)
00679       o2sato = EXP(CO)
00680 !
00681 !  Convert from ml/l to mmol/m^-3 for ROMS
00682 !
00683       o2sato = o2sato/22.3916*1000.0
00684       return
00685       end
00686
/home/cholod/wd/peps/ROMS/romsagrif_r215/Roms_tools/Roms_Agrif/bio_NPZD.F
