source: branches/dev_001_GM/NEMO/TOP_SRC/PISCES/trcini.pisces.h90 @ 764

   !!======================================================================
   !!                         ***  trcini.pisces.h90  ***
   !! TOP :   Initialisation of PISCES biological model
   !!======================================================================
   !! History :    -   !  1988-07  (E. Maier-Reiner) Original code
   !!              -   !  1999-10  (O. Aumont, C. Le Quere)
   !!              -   !  2002     (O. Aumont)  PISCES
   !!             1.0  !  2005-03  (O. Aumont, A. El Moussaoui) F90
   !!----------------------------------------------------------------------

#  include "domzgr_substitute.h90"
#  include "passivetrc_substitute.h90"
   !!----------------------------------------------------------------------
   !! NEMO/TOP 1.0 , LOCEAN-IPSL (2005) 
   !! $Id$ 
   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
   !!----------------------------------------------------------------------

CONTAINS

   SUBROUTINE trc_ini
      !!----------------------------------------------------------------------
      !!                   ***  ROUTINE trc_ini ***
      !!
      !! ** Purpose :   Initialisation of PISCES biological and chemical variables
      !!----------------------------------------------------------------------
      USE iom
      !!
      INTEGER :: ji,jj,jk
      INTEGER :: ichl,iband,jm
      INTEGER , PARAMETER :: jpmois = 12, jpan   = 1 

      REAL(wp) :: zcoef
      REAL(wp) :: ztoto,expide,denitide,ztra,zmaskt
      REAL(wp) , DIMENSION (jpi,jpj) :: riverdoc,river,ndepo
      REAL(wp) , DIMENSION (jpi,jpj,jpk) :: cmask

      INTEGER :: numriv, numdust, numbath, numdep
      INTEGER :: numlight 

#if defined key_trc_kriest
      INTEGER  ::   jn, kiter
      REAL(wp) ::   znum, zdiv
      REAL(wp) ::   zws,zwr, zwl,wmax, xnummax, &
      REAL(wp) ::   zmin, zmax, zl, zr, xacc
#endif
      !!----------------------------------------------------------------------


      IF(lwp) WRITE(numout,*)
      IF(lwp) WRITE(numout,*) ' trc_ini :   PISCES biological and chemical initialisation'
      IF(lwp) WRITE(numout,*) ' ~~~~~~~'


      !                                            ! Time-step
      rfact   = rdttra(1) * float(ndttrc)          ! ---------
      rfactr  = 1. / rfact
      rfact2  = rfact / float(nrdttrc)
      rfact2r = 1. / rfact2

      IF(lwp) WRITE(numout,*) '    Tracer  time step=', rfact, ' rdt = ', rdt
      IF(lwp) write(numout,*) '    Biology time step=', rfact2


      !                                            ! Dust input from the atmosphere
      IF( bdustfer ) THEN                          ! ------------------------------
         IF(lwp) WRITE(numout,*) '    Initialize dust input from atmosphere '
         CALL iom_open ( 'dust.orca.nc', numdust )
         DO jm = 1, jpmois
            CALL iom_get( numdust, jpdom_data, 'dust', dustmo(:,:,jm), jm )
         END DO
         CALL iom_close( numdust )
      ELSE
         dustmo(:,:,:) = 0.e0
      ENDIF


      !                                            ! Nutrient input from rivers
      IF( briver ) THEN                            ! --------------------------
         IF(lwp) WRITE(numout,*) '    Initialize the nutrient input by rivers from river.orca.nc file'
         CALL iom_open ( 'river.orca.nc', numriv )
         CALL iom_get  ( numriv, jpdom_data, 'riverdic', river   (:,:), jpan )
         CALL iom_get  ( numriv, jpdom_data, 'riverdoc', riverdoc(:,:), jpan )
         CALL iom_close( numriv )
      ELSE
         river   (:,:) = 0.e0
         riverdoc(:,:) = 0.e0
      endif

      !                                            ! Nutrient input from dust
      IF( bndepo ) THEN                            ! ------------------------
         IF(lwp) WRITE(numout,*) '    Initialize the nutrient input by dust from ndeposition.orca.nc'
         CALL iom_open ( 'ndeposition.orca.nc', numdep )
         CALL iom_get  ( numdep, jpdom_data, 'ndep', ndepo(:,:), jpan )
         CALL iom_close( numdep )
      ELSE
         ndepo(:,:) = 0.e0
      ENDIF

      !                                            ! Coastal and island masks
      IF( bsedinput ) THEN                         ! ------------------------
         IF(lwp) WRITE(numout,*) '    Computation of an island mask to enhance coastal supply of iron '
         IF(lwp) WRITE(numout,*) '       from bathy.orca.nc file '
         CALL iom_open ( 'bathy.orca.nc', numbath )
         CALL iom_get  ( numbath, jpdom_data, 'bathy', cmask(:,:,:), jpan )
         CALL iom_close( numbath )
         !
         DO jk = 1, 5
            DO jj = 2, jpjm1
               DO ji = 2, jpim1
                  IF( tmask(ji,jj,jk) /= 0. ) THEN
                     zmaskt = tmask(ji+1,jj,jk) * tmask(ji-1,jj,jk) * tmask(ji,jj+1,jk)    &
                        &                       * tmask(ji,jj-1,jk) * tmask(ji,jj,jk+1)
                     IF( zmaskt == 0. )   cmask(ji,jj,jk ) = 0.1
                  ENDIF
               END DO
            END DO
         END DO
         DO jk = 1, jpk
            DO jj = 1, jpj
               DO ji = 1, jpi
                  expide   = MIN( 8.,( fsdept(ji,jj,jk) / 500. )**(-1.5) )
                  denitide = -0.9543 + 0.7662 * LOG( expide ) - 0.235 * LOG( expide )**2
                  cmask(ji,jj,jk) = cmask(ji,jj,jk) * MIN( 1., EXP( denitide ) / 0.5 )
               END DO
            END DO
         END DO    
      ELSE
         cmask(:,:,:) = 0.e0
      ENDIF

      CALL lbc_lnk( cmask , 'T', 1. )      ! Lateral boundary conditions on cmask   (sign unchanged)


      !                                            ! total atmospheric supply of Si
      !                                            ! ------------------------------
      sumdepsi = 0.e0
      DO jm = 1, jpmois
         DO jj = 2, jpjm1
            DO ji = 2, jpim1
               sumdepsi = sumdepsi + dustmo(ji,jj,jm) / (12.*rmoss) * 8.8        &
                  &     * 0.075/28.1 * e1t(ji,jj) * e2t(ji,jj) * tmask(ji,jj,1) * tmask_i(ji,jj)
            END DO
         END DO
      END DO
      IF( lk_mpp )   CALL mpp_sum( sumdepsi )  ! sum over the global domain

      !                                            ! N/P and Si releases due to coastal rivers
      !                                            ! -----------------------------------------
      DO jj = 1, jpj
         DO ji = 1, jpi
            zcoef = raass * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,1)
            cotdep(ji,jj) =  river(ji,jj)                  *1E9 / ( 12. * zcoef + rtrn ) * tmask(ji,jj,1)
            rivinp(ji,jj) = (river(ji,jj)+riverdoc(ji,jj)) *1E9 / ( 31.6* zcoef + rtrn ) * tmask(ji,jj,1)
            nitdep(ji,jj) = 7.6 * ndepo(ji,jj)                  / ( 14E6*raass*fse3t(ji,jj,1) + rtrn ) * tmask(ji,jj,1)
         END DO
      END DO
      ! Lateral boundary conditions on ( cotdep, rivinp, nitdep )   (sign unchanged)
      CALL lbc_lnk( cotdep , 'T', 1. )  ;  CALL lbc_lnk( rivinp , 'T', 1. )  ;  CALL lbc_lnk( nitdep , 'T', 1. )

      rivpo4input=0.e0
      rivalkinput=0.e0
      nitdepinput=0.e0
      DO jj = 2 , jpjm1
         DO ji = 2, jpim1
            zcoef = e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,1)) * tmask(ji,jj,1) * tmask_i(ji,jj) * raass
            rivpo4input = rivpo4input + rivinp(ji,jj) * zcoef
            rivalkinput = rivalkinput + cotdep(ji,jj) * zcoef
            nitdepinput = nitdepinput + nitdep(ji,jj) * zcoef
         END DO
      END DO
      IF( lk_mpp ) THEN
         CALL mpp_sum( rivpo4input )  ! sum over the global domain
         CALL mpp_sum( rivalkinput )  ! sum over the global domain
         CALL mpp_sum( nitdepinput )  ! sum over the global domain
      ENDIF


      !                                            !    Coastal supply of iron
      !                                            ! -------------------------
      DO jk = 1, jpkm1
         ironsed(:,:,jk) = sedfeinput * cmask(:,:,jk) / ( fse3t(:,:,jk) * rjjss )
      END DO
      CALL lbc_lnk( ironsed , 'T', 1. )      ! Lateral boundary conditions on ( ironsed )   (sign unchanged)



#if defined key_trc_kriest
      !----------------------------------------------------------------------
      !  COMPUTATION OF THE VERTICAL PROFILE OF MAXIMUM SINKING SPEED
      !  Search of the maximum number of particles in aggregates for each k-level.  
      !  Bissection Method
      !--------------------------------------------------------------------
      WRITE(numout,*)
      WRITE(numout,*)'    kriest : Compute maximum number of particles in aggregates'
            
      xacc     =  0.001
      kiter    = 50
      zmin     =  1.10
      zmax     = xkr_mass_max / xkr_mass_min
      xkr_frac = zmax

      DO jk =1,jpk
         zl = zmin
         zr = zmax
         wmax = 0.5 * fse3t(1,1,jk) * rjjss / rfact2
         zdiv = xkr_zeta + xkr_eta - xkr_eta * zl
         znum = zl - 1.
         zwl =  xkr_wsbio_min * xkr_zeta / zdiv &
            & - ( xkr_wsbio_max * xkr_eta * znum * &
            &     xkr_frac**( -xkr_zeta / znum ) / zdiv ) &
            & - wmax
  
         zdiv = xkr_zeta + xkr_eta - xkr_eta * zr
         znum = zr - 1.
         zwr =  xkr_wsbio_min * xkr_zeta / zdiv &
            & - ( xkr_wsbio_max * xkr_eta * znum * &
            &     xkr_frac**( -xkr_zeta / znum ) / zdiv ) &
            & - wmax

iflag:  DO jn = 1, kiter               
           IF( zwl == 0.e0 ) THEN
              xnummax = zl
           ELSE IF ( zwr == 0.e0 ) THEN
              xnummax = zr
           ELSE
              xnummax = ( zr + zl ) / 2.
              zdiv = xkr_zeta + xkr_eta - xkr_eta * xnummax
              znum = xnummax - 1.
              zws =  xkr_wsbio_min * xkr_zeta / zdiv &
                 & - ( xkr_wsbio_max * xkr_eta * znum * &
                 &     xkr_frac**( -xkr_zeta / znum ) / zdiv ) &
                 & - wmax
              IF( zws * zwl < 0. ) THEN
                 zr = xnummax
              ELSE
                 zl = xnummax
              ENDIF
              zdiv = xkr_zeta + xkr_eta - xkr_eta * zl
              znum = zl - 1.
              zwl =  xkr_wsbio_min * xkr_zeta / zdiv &
                 & - ( xkr_wsbio_max * xkr_eta * znum * &
                 &     xkr_frac**( -xkr_zeta / znum ) / zdiv ) &
                 & - wmax
              
              zdiv = xkr_zeta + xkr_eta - xkr_eta * zr
              znum = zr - 1.
              zwr =  xkr_wsbio_min * xkr_zeta / zdiv &
                 & - ( xkr_wsbio_max * xkr_eta * znum * &
                 &     xkr_frac**( -xkr_zeta / znum ) / zdiv ) &
                 & - wmax

              IF ( ABS ( zws )  <= xacc ) EXIT iflag

           ENDIF
                    
        END DO iflag
               
        xnumm(jk) = xnummax
        WRITE(numout,*) '       jk = ', jk, ' wmax = ', wmax,' xnum max = ', xnumm(jk)
        
     END DO

#endif

      !----------------------------------------------------------------------
      ! Initialize biological variables 
      !----------------------------------------------------------------------
      ! Set biological ratios
      ! ---------------------
      rno3   = (16.+2.) / 122.
      po4r   =   1.e0   / 122.
      o2nit  =  32.     / 122.
      rdenit =  97.6    /  16.
      o2ut   = 140.     / 122.

      !----------------------------------------------------------------------
      ! Initialize chemical variables 
      !----------------------------------------------------------------------

      ! set pre-industrial atmospheric [co2] (ppm) and o2/n2 ratio
      ! ----------------------------------------------------------
      atcox = 0.20946

      ! Set lower/upper limits for temperature and salinity
      ! ---------------------------------------------------
      salchl = 1.e0 / 1.80655
      calcon = 1.03e-2

      ! Set coefficients for apparent solubility equilibrium of calcite
      ! Millero et al. 1995 from Mucci 1983
      ! --------------------------------------------------------------
      akcc1 = -171.9065
      akcc2 =   -0.077993
      akcc3 = 2839.319
      akcc4 =   71.595
      akcc5 =   -0.77712
      akcc6 =    0.0028426
      akcc7 =  178.34
      akcc8 =   -0.07711
      akcc9 =    0.0041249

      ! Set coefficients for seawater pressure correction
      ! -------------------------------------------------
      devk1(1) = -25.5
      devk2(1) =   0.1271
      devk3(1) =   0.e0
      devk4(1) =  -3.08E-3
      devk5(1) =   0.0877E-3
      !
      devk1(2) = -15.82
      devk2(2) =  -0.0219
      devk3(2) =   0.e0
      devk4(2) =   1.13E-3
      devk5(2) =  -0.1475E-3
      !
      devk1(3) = -29.48
      devk2(3) =   0.1622
      devk3(3) =   2.608E-3
      devk4(3) =  -2.84E-3
      devk5(3) =   0.e0
      !
      devk1(4) = -14.51
      devk2(4) =   0.1211
      devk3(4) =  -0.321E-3
      devk4(4) =  -2.67E-3
      devk5(4) =   0.0427E-3
      !
      devk1(5) = -23.12
      devk2(5) =   0.1758
      devk3(5) =  -2.647E-3
      devk4(5) =  -5.15E-3
      devk5(5) =   0.09E-3
      !
      devk1(6) = -26.57
      devk2(6) =   0.2020
      devk3(6) =  -3.042E-3
      devk4(6) =  -4.08E-3
      devk5(6) =   0.0714E-3
      !
      devk1(7) = -25.60
      devk2(7) =   0.2324
      devk3(7) =  -3.6246E-3
      devk4(7) =  -5.13E-3
      devk5(7) =   0.0794E-3
      !
      ! For calcite with Edmond and Gieske 1970
      !     devkst = 0.23
      !     devks  = 35.4
      ! Millero 95 takes this depth dependance for calcite
      devk1(8) = -48.76
      devk2(8) =   0.5304
      devk3(8) =   0.e0
      devk4(8) = -11.76E-3
      devk5(8) =   0.3692E-3
      !
      ! Coefficients for sulfate and fluoride
      devk1(9) = -18.03
      devk2(9) =   0.0466
      devk3(9) =   0.316e-3
      devk4(9) =  -4.53e-3
      devk5(9) =   0.09e-3

      devk1(10) = -9.78
      devk2(10) = -0.0090
      devk3(10) = -0.942e-3
      devk4(10) = -3.91e-3
      devk5(10) =  0.054e-3


      ! Set universal gas constants
      ! ---------------------------
      rgas  = 83.143
      oxyco =  1.e0 / 22.4144

      ! Set boron constants
      ! -------------------
      bor1 = 0.00023
      bor2 = 1.e0 / 10.82

      ! Set volumetric solubility constants for co2 in ml/l (Weiss, 1974)
      ! -----------------------------------------------------------------
      c00 = -60.2409
      c01 =  93.4517
      c02 =  23.3585
      c03 =   0.023517
      c04 =  -0.023656
      c05 =   0.0047036
      !
      ca0 = -162.8301
      ca1 =  218.2968
      ca2 =   90.9241
      ca3 =   -1.47696
      ca4 =    0.025695
      ca5 =   -0.025225
      ca6 =    0.0049867

      ! Set coeff. for 1. dissoc. of carbonic acid (Edmond and Gieskes, 1970)
      ! ---------------------------------------------------------------------
      c10 = -3670.7
      c11 =    62.008
      c12 =    -9.7944
      c13 =     0.0118
      c14 =    -0.000116

      ! Set coeff. for 2. dissoc. of carbonic acid (Edmond and Gieskes, 1970)
      ! ---------------------------------------------------------------------
      c20 = -1394.7
      c21 =    -4.777
      c22 =     0.0184
      c23 =    -0.000118

      ! Set constants for calculate concentrations for sulfate and fluoride
      ! sulfates (Morris & Riley 1966)
      !----------------------------------------------------------------------
      st1 = 0.14
      st2 = 1.e0 / 96.062

      ! fluoride
      ! --------
      ft1 = 0.000067
      ft2 = 1.e0 / 18.9984

      ! sulfates (Dickson 1990 change to mol:kg soln, idem OCMIP)
      !----------------------------------------------------------
      ks0  =    141.328
      ks1  =  -4276.1
      ks2  =    -23.093
      ks3  = -13856.
      ks4  =    324.57
      ks5  =    -47.986
      ks6  =  35474.
      ks7  =   -771.54
      ks8  =    114.723
      ks9  =  -2698.
      ks10 =   1776.
      ks11 =      1.
      ks12 =     -0.001005

      ! fluorides (Dickson & Riley 1979 change to mol/kg soln)
      !-------------------------------------------------------
      kf0 =  -12.641
      kf1 = 1590.2
      kf2 =    1.525
      kf3 =    1.0
      kf4 =   -0.001005

      ! Set coeff. for 1. dissoc. of boric acid (Edmond and Gieskes, 1970)
      ! ------------------------------------------------------------------
      cb0  = -8966.90
      cb1  = -2890.53
      cb2  =   -77.942
      cb3  =     1.728
      cb4  =    -0.0996
      cb5  =   148.0248
      cb6  =   137.1942
      cb7  =     1.62142
      cb8  =   -24.4344
      cb9  =   -25.085
      cb10 =    -0.2474
      cb11 =     0.053105

      ! Set coeff. for dissoc. of water (Dickson and Riley, 1979, 
      !   eq. 7, coefficient cw2 corrected from 0.9415 to 0.09415 
      !   after pers. commun. to B. Bacastow, 1988)
      ! ---------------------------------------------------------
      cw0 = -13847.26
      cw1 =    148.9652
      cw2 =    -23.6521
      cw3 =    118.67
      cw4 =     -5.977
      cw5 =      1.0495
      cw6 =     -0.01615

      ! Set coeff. for dissoc. of phosphate (Millero (1974)
      ! ---------------------------------------------------
      cp10 =   115.54
      cp11 = -4576.752
      cp12 =   -18.453
      cp13 =  -106.736
      cp14 =     0.69171
      cp15 =    -0.65643
      cp16 =    -0.01844
      !
      cp20 =   172.1033
      cp21 = -8814.715
      cp22 =   -27.927
      cp23 =  -160.340
      cp24 =     1.3566
      cp25 =     0.37335
      cp26 =    -0.05778
      !
      cp30 =   -18.126
      cp31 = -3070.75
      cp32 =    17.27039
      cp33 =     2.81197
      cp34 =   -44.99486
      cp35 =    -0.09984

      ! Set coeff. for dissoc. of phosphate (Millero (1974)
      ! ---------------------------------------------------
      cs10 =   117.385
      cs11 = -8904.2
      cs12 =   -19.334
      cs13 =  -458.79
      cs14 =     3.5913
      cs15 =   188.74
      cs16 =    -1.5998
      cs17 =   -12.1652
      cs18 =     0.07871
      cs19 =     0.e0
      cs20 =     1.e0
      cs21 =    -0.001005


      ! Set volumetric solubility constants for o2 in ml/l (Weiss, 1970)
      ! ----------------------------------------------------------------
      ox0 = -58.3877
      ox1 =  85.8079
      ox2 =  23.8439
      ox3 =  -0.034892
      ox4 =   0.015568
      ox5 =  -0.0019387

      !  FROM THE NEW BIOOPTIC MODEL PROPOSED JM ANDRE, WE READ HERE
      !  A PRECOMPUTED ARRAY CORRESPONDING TO THE ATTENUATION COEFFICIENT

      CALL ctlopn( numlight, 'kRGB61.txt', 'OLD', 'FORMATTED', 'SEQUENTIAL',   &
         &           1, numout, .TRUE., 1 )
      DO ichl = 1,61
         READ(numlight,*) ztoto, ( xkrgb(iband,ichl), iband = 1,3 )
      END DO
      CLOSE(numlight)


      CALL p4zche        ! initialize the chemical constants


      ndayflxtr = 0      !  Initialize a counter for the computation of chemistry


      IF(lwp) WRITE(numout,*) ' Initialisation of PISCES done'
      !
   END SUBROUTINE trc_ini