CCC-----------------------------------------------------------------
CCC
CCC                     ROUTINE trcini.pisces.h
CCC                     ************************
CCC
CCC  PURPOSE :
CCC  ---------
CCC     Initialisation of PISCES biological and chemical variables
CCC
CC   INPUT :
CC   -----
CC      common
CC              all the common defined in opa 
CC
CC
CC   OUTPUT :                   : no
CC   ------
CC
CC   EXTERNAL :
CC   ----------
CC         p4zche
CC
CC   MODIFICATIONS:
CC   --------------
CC      original  : 1988-07  E. MAIER-REIMER      MPI HAMBURG
CC      additions : 1999-10  O. Aumont and C. Le Quere
CC      additions : 2002     O. Aumont (PISCES)
CC
CCC---------------------------------------------------------------------
CC local declarations
CC ==================
      INTEGER ichl,iband,mo
      INTEGER jpmois,jpan
      PARAMETER (jpan=1, jpmois=12)

      REAL xtoto,expide,denitide
      REAL ndepo(jpi,jpj),river(jpi,jpj)
      REAL riverdoc(jpi,jpj)
      CHARACTER*34 clname

      INTEGER ipi,ipj,ipk,istep(jpmois),itime,istep0(jpan)
      REAL zsecond, zlon(jpi,jpj),zlat(jpi,jpj),zlev(jpk),zdate0
      INTEGER numriv,numdust,numbath,numdep
C
C 1. initialization
C -----------------
C
C computation of the record length for direct access FILE
C this length depend of 512 for the t3d machine
C
      rfact = rdttra(1) * float(ndttrc)
      rfactr = 1./rfact
      WRITE(numout,*) ' Tracer time step=',rfact,' rdt=',rdt
      rfact2= rfact / float(nrdttrc)
      rfact2r = 1./rfact2
      write(numout,*) ' Biology time step=',rfact2
C
C    INITIALISE DUST INPUT FROM ATMOSPHERE
C    -------------------------------------
C
       IF (bdustfer) THEN
         clname='dust.orca.nc'
       CALL flinopen(clname,mig(1),nlci,mjg(1),nlcj,.false.,ipi,ipj,0
     $      ,zlon,zlat,zlev,itime,istep,zdate0,zsecond,numdust)
       CALL flinget(numdust,'dust',jpidta,jpjdta,0,jpmois,1,
     $        12,mig(1),nlci,mjg(1),nlcj,dustmo(1:nlci,1:nlcj,12) )
       CALL flinclo(numdust)
       ELSE
       dustmo(:,:,:)=0.
       ENDIF
C
C    INITIALISE THE NUTRIENT INPUT BY RIVERS
C    ---------------------------------------
C
       IF (briver) THEN
        clname='river.orca.nc'
       CALL flinopen(clname,mig(1),nlci,mjg(1),nlcj,.false.,ipi,ipj,0
     $      ,zlon,zlat,zlev,itime,istep0,zdate0,zsecond,numriv)
       CALL flinget(numriv,'riverdic',jpidta,jpjdta,0,jpan,1,
     $        1,mig(1),nlci,mjg(1),nlcj,river(1:nlci,1:nlcj) )
       CALL flinget(numriv,'riverdoc',jpidta,jpjdta,0,jpan,1,
     $        1,mig(1),nlci,mjg(1),nlcj,riverdoc(1:nlci,1:nlcj) )
       CALL flinclo(numriv)
       ELSE
       river(:,:)=0.
       riverdoc(:,:)=0.
       endif
C
C    INITIALISE THE N INPUT BY DUST
C    ---------------------------------------
C
       IF (bndepo) THEN
        clname='ndeposition.orca.nc'
       CALL flinopen(clname,mig(1),nlci,mjg(1),nlcj,.false.,ipi,ipj,0
     $      ,zlon,zlat,zlev,itime,istep0,zdate0,zsecond,numdep)
       CALL flinget(numdep,'ndep',jpidta,jpjdta,0,jpan,1,
     $        1,mig(1),nlci,mjg(1),nlcj,ndepo(1:nlci,1:nlcj) )
       CALL flinclo(numdep)
       ELSE
       ndepo(:,:)=0.
       ENDIF
C
C    Computation of the coastal mask.
C    Computation of an island mask to enhance coastal supply
C    of iron
C    -------------------------------------------------------
C
       IF (bsedinput) THEN
        clname='bathy.orca.nc'
       CALL flinopen(clname,mig(1),nlci,mjg(1),nlcj,.false.,ipi,ipj,ipk
     $      ,zlon,zlat,zlev,itime,istep0,zdate0,zsecond,numbath)
       CALL flinget(numbath,'bathy',jpidta,jpjdta,jpk,jpan,1,
     $        1,mig(1),nlci,mjg(1),nlcj,cmask(1:nlci,1:nlcj,1:jpk) )
       CALL flinclo(numbath)
C
       DO jk=1,jpk
         DO ji=1,jpi
           DO jj=1,jpj
           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)*exp(denitide)/0.6858
           END DO
         END DO
       END DO
       ELSE
       cmask(:,:,:)=0.
       ENDIF
C
C     Computation of the total atmospheric supply of Si
C     -------------------------------------------------
C
       sumdepsi=0.
       DO mo=1,12
         DO jj=2,jpjm1
           DO ji=2,jpim1
           sumdepsi=sumdepsi+dustmo(ji,jj,mo)/(12.*rmoss)*8.8
     &       *0.075/28.1*e1t(ji,jj)*e2t(ji,jj)*tmask(ji,jj,1)
           END DO
         END DO
       END DO
C
C    COMPUTATION OF THE N/P RELEASE DUE TO COASTAL RIVERS
C    COMPUTATION OF THE Si RELEASE DUE TO COASTAL RIVERS 
C    ---------------------------------------------------
C
       DO jj=1,jpj
         DO ji=1,jpi
       cotdep(ji,jj,1)=river(ji,jj)*1E9/(12.*raass
     &   *e1t(ji,jj)*e2t(ji,jj)*fse3t(ji,jj,1))*tmask(ji,jj,1)
       po4dep(ji,jj,1)=(river(ji,jj)+riverdoc(ji,jj))*1E9
     &   /(31.6*raass*e1t(ji,jj)*e2t(ji,jj)*fse3t(ji,jj,1))
     &   *tmask(ji,jj,1)
       nitdep(ji,jj,1)=7.6*ndepo(ji,jj)*tmask(ji,jj,1)/(14E6*raass
     &   *fse3t(ji,jj,1))
         END DO
       END DO

      rivpo4input=0.
      rivalkinput=0.
      rivnitinput=0.
       DO jj=2,jpjm1
         DO ji=2,jpim1
         rivpo4input=rivpo4input+po4dep(ji,jj,1)*(e1t(ji,jj)*e2t(ji,jj)
     &     *fse3t(ji,jj,1))*tmask(ji,jj,1)*raass
         rivalkinput=rivalkinput+cotdep(ji,jj,1)*(e1t(ji,jj)*e2t(ji,jj)
     &     *fse3t(ji,jj,1))*tmask(ji,jj,1)*raass
         rivnitinput=rivnitinput+nitdep(ji,jj,1)*(e1t(ji,jj)*e2t(ji,jj)
     &     *fse3t(ji,jj,1))*tmask(ji,jj,1)*raass
         END DO
       END DO
C
C    Coastal supply of iron
C    ----------------------
C
      DO jk=1,jpkm1
        DO jj=1,jpj
          DO ji=1,jpi
          ironsed(ji,jj,jk)=sedfeinput*cmask(ji,jj,jk)
     &      /(fse3t(ji,jj,jk)*rjjss)
          END DO
        END DO
      END DO
C
CC----------------------------------------------------------------------
CC
CC Initialize biological variables 
CC
CC----------------------------------------------------------------------
C
      spocri = 0.003
      jkopt = 14
C
C Set biological ratios
C ---------------------
C
      rno3   = (16.+2.)/122.
      po4r   = 1./122.
      o2ut   = 172./122.
      o2nit  = 32./122.
      rdenit = 97.6/16.
      o2ut   = 140./122.
C
CC----------------------------------------------------------------------
CC
CC Initialize chemical variables 
CC
CC----------------------------------------------------------------------
C
C set pre-industrial atmospheric [co2] (ppm) and o2/n2 ratio
C ----------------------------------------------------------
C
      atcox = 0.20946
C
C Set lower/upper limits for temperature and salinity
C ---------------------------------------------------
C
      salchl = 1./1.80655
      calcon = 1.03E-2
C
C Set coefficients for apparent solubility equilibrium
C   of calcite (Ingle, 1800, eq. 6)
C ----------------------------------------------------
C
      akcc1 = -34.452
      akcc2 = -39.866
      akcc3 = 110.21
      akcc4 = -7.5752E-6
C
C
C Set coefficients for seawater pressure correction
C -------------------------------------------------
C
      devk1  = 24.2
      devk2  = 16.4
      devkb  = 27.5
      devk1t = 0.085
      devk2t = 0.04
      devkbt = 0.095
C
      devkst = 0.23
      devks  = 35.4
C
C Set universal gas constants
C ---------------------------
C
      rgas = 83.143
      oxyco = 1./22.4144
C
C Set boron constants
C -------------------
C
      bor1 = 0.00023
      bor2 = 1./10.82
C
C Set volumetric solubility constants for co2 in ml/l (Weiss, 1974)
C -----------------------------------------------------------------
C
      c00 = -58.0931
      c01 = 90.5069
      c02 = 22.2940
      c03 = 0.027766
      c04 = -0.025888
      c05 = 0.0050578
C
C Set coeff. for 1. dissoc. of carbonic acid (Edmond and Gieskes, 1970)
C ---------------------------------------------------------------------
C
      c10 = -2307.1266
      c11 = 2.83655
      c12 = -1.5529413
      c13 = -4.0484
      c14 = -0.20760841
      c15 = 0.08468345
      c16 = -0.00654208
      c17 = -0.001005
C
C Set coeff. for 2. dissoc. of carbonic acid (Edmond and Gieskes, 1970)
C ---------------------------------------------------------------------
C
      c20 = -3351.6106
      c21 = -9.226508
      c22 = -0.2005743
      c23 = -23.9722
      c24 = -0.106901773
      c25 = 0.1130822
      c26 = -0.00846934
      c27 = -0.001005
C
C Set coeff. for 1. dissoc. of boric acid (Edmond and Gieskes, 1970)
C ------------------------------------------------------------------
C
      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
C
C Set coeff. for dissoc. of water (Dickson and Riley, 1979, 
C   eq. 7, coefficient cw2 corrected from 0.9415 to 0.09415 
C   after pers. commun. to B. Bacastow, 1988)
C ---------------------------------------------------------
C
      cw0 = -13847.26
      cw1 = 148.9652
      cw2 = -23.6521
      cw3 = 118.67
      cw4 = -5.977
      cw5 = 1.0495
      cw6 = -0.01615
C
C Set volumetric solubility constants for o2 in ml/l (Weiss, 1970)
C ----------------------------------------------------------------
C
      ox0 = -58.3877
      ox1 = 85.8079
      ox2 = 23.8439
      ox3 = -0.034892
      ox4 = 0.015568
      ox5 = -0.0019387
C
C  FROM THE NEW BIOOPTIC MODEL PROPOSED JM ANDRE, WE READ HERE
C  A PRECOMPUTED ARRAY CORRESPONDING TO THE ATTENUATION COEFFICIENT
C
         open(49,file='kRGB61.txt',form='formatted')
         do ichl=1,61
           READ(49,*) xtoto,(xkrgb(iband,ichl),iband = 1,3)
         end do
         close(49)
C      
#if defined key_off_degrad
C
C Read volume for degraded regions (DEGINIT)
C ------------------------------------------
C
#    if defined key_vpp
      CALL READ3S(902,facvol,jpi,jpj,jpk)
#    else
      READ (902) facvol
#    endif
#endif
C
C
C  Call p4zche to initialize the chemical constants
C  ------------------------------------------------
C
      CALL p4zche
C
      WRITE(numout,*) ' Initialisation of PISCES done'