module m_ECOSM_initialize_ecosys

      contains

      subroutine initialize_ecosys()
      use mod_xc
      use mod_za
      use mod_necessary_ecovars
      use m_ECOSM_bioini
      use m_ECOSM_riverloads

      implicit none
      include 'common_blocks.h'
      include 'biovar.h'
      include 'ECOSMparam1.h'
      integer :: o
      integer lslot,nk
      character preambl(5)*79
      integer lgth
      real tiny
      call ECOSM_bioini
      !KAL - Should use existing margin, points up to ii+margin should
      !KAL - still be valid at this point 
      !margin=0
      tiny=2.23e-16

!    initialize nitrate from climatology files
      lgth=len_trim(flnmfor)
      call zaiopf(flnmforw(1:lgth)//'relax.nitr.a','old',931)
      if (mnproc.eq.1) then   ! .b file from 1st tile only
        open (unit=2931,file=flnmforw(1:lgth)//'relax.nitr.b',
     &      status='old',action='read')
        read(2931,'(a79)') preambl
      endif  !1st tile
      call preambl_print(preambl)

      do k=1,kdm
        call rdmonthck(bio(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy,k,1,init),
     &                 931,1)
      enddo   
      close(931)
!     from mmolN to mgrC using redfield stoichiometrie:
      bio(:,:,:,1,init)=bio(:,:,:,1,init)*REDF(1)*REDF(6)
      bio(:,:,:,2,init)=bio(:,:,:,1,init)

!
!    initialize phosphate from climatology files
      lgth=len_trim(flnmfor)
      call zaiopf(flnmforw(1:lgth)//'relax.phos.a','old',932)
      if (mnproc.eq.1) then   ! .b file from 1st tile only
        open (unit=2932,file=flnmforw(1:lgth)//'relax.phos.b',
     &      status='old',action='read')
        read(2932,'(a79)') preambl
      endif  !1st tile
      call preambl_print(preambl)

      do k=1,kdm
        call rdmonthck(bio(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy,k,1,ipho),
     &                 932,1)
      enddo   
      close(932)
!     from mmolP to mgrC:      
      bio(:,:,:,1,ipho)=bio(:,:,:,1,ipho)*REDF(2)*REDF(6)
      bio(:,:,:,2,ipho)=bio(:,:,:,1,ipho)
!      

!    initialize silicate from climatology files
      lgth=len_trim(flnmfor)
      call zaiopf(flnmforw(1:lgth)//'relax.sili.a','old',933)
      if (mnproc.eq.1) then   ! .b file from 1st tile only
        open (unit=2933,file=flnmforw(1:lgth)//'relax.sili.b',
     &      status='old',action='read')
        read(2933,'(a79)') preambl
      endif  !1st tile
      call preambl_print(preambl)

      do k=1,kdm
        call rdmonthck(bio(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy,k,1,isil),
     &                 933,1)
      enddo   
      close(933)
!     from mmolSI to mgrC:
      bio(:,:,:,1,isil)=bio(:,:,:,1,isil)*REDF(3)*REDF(6)
      bio(:,:,:,2,isil)=bio(:,:,:,1,isil)

!    initialize oxygen from climatology files
      lgth=len_trim(flnmfor)
      call zaiopf(flnmforw(1:lgth)//'relax.oxyg.a','old',934)
      if (mnproc.eq.1) then   ! .b file from 1st tile only
        open (unit=2934,file=flnmforw(1:lgth)//'relax.oxyg.b',
     &      status='old',action='read')
        read(2934,'(a79)') preambl
      endif  !1st tile
      call preambl_print(preambl)

      do k=1,kdm
        call rdmonthck(bio(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy,k,1,ioxy),
     &                 934,1)
      enddo   
      close(934)
!     from input in ml/l to mmol/m3:
      bio(:,:,:,1,ioxy)=bio(:,:,:,1,ioxy)*REDF(7)
      bio(:,:,:,2,ioxy)=bio(:,:,:,1,ioxy)


      ! KAL - initialization, so do the whole array
      do i=1-nbdy,idm+nbdy
      do j=1-nbdy,jdm+nbdy
      do k=1,kdm
      do o=1,nbio
         if(bio(i,j,k,1,o) .lt. 0.) then
            bio(i,j,k,1,o)=tiny
         endif   
         if(bio(i,j,k,2,o) .lt. 0.) then
            bio(i,j,k,2,o)=tiny
         endif
!Oxygene becomes negative when it is small   
         if (o==ioxy) then
            if (bio(i,j,k,1,o).le.44.0) then
               bio(i,j,k,1,o)=44.0
               bio(i,j,k,2,o)=44.0
            endif
         endif
      enddo   
      enddo   
      enddo   
      enddo  

      bio(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy,1:kdm,1:2,idet)=0.1
      bio(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy,1:kdm,1:2,ifla)=0.1
      bio(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy,1:kdm,1:2,idia)=0.1
      bio(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy,1:kdm,1:2,imeso)=0.1 
      bio(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy,1:kdm,1:2,imicro)=0.1 
      bio(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy,1:kdm,1:2,inh4)=0.1 
      bio(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy,1:kdm,1:2,idom)=0.1 
      bio(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy,1:kdm,1:2,iopa)=0.1 
      bio(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy,1:kdm,1:2,ibg)=0.1 
      bio(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy,1:kdm,1:2,ifish)=0.1 
      bio(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy,1:kdm,1:2,ichlf)=0.1/50.
      bio(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy,1:kdm,1:2,ichld)=0.1/50.
      bio(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy,1:kdm,1:2,ichlb)=0.1/50.
     
#ifdef ECO2
!    initialize alkalinity from climatology files
      lgth=len_trim(flnmfor)
      call zaiopf(flnmforw(1:lgth)//'relax.alkl.a','old',935)
      if (mnproc.eq.1) then   ! .b file from 1st tile only
        open (unit=2935,file=flnmforw(1:lgth)//'relax.alkl.b',
     &      status='old',action='read')
        read(2935,'(a79)') preambl
      endif  !1st tile
      call preambl_print(preambl)

      do k=1,kdm
        call rdmonthck(bio(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy,k,1,ialk),
     &                 935,1)
      enddo   
      close(935)
      bio(:,:,:,2,ialk)=bio(:,:,:,1,ialk)
!   
!    initialize DIC from climatology files
      lgth=len_trim(flnmfor)
      call zaiopf(flnmforw(1:lgth)//'relax.dico.a','old',936)
      if (mnproc.eq.1) then   ! .b file from 1st tile only
        open (unit=2936,file=flnmforw(1:lgth)//'relax.dico.b',
     &      status='old',action='read')
        read(2936,'(a79)') preambl
      endif  !1st tile
      call preambl_print(preambl)

      do k=1,kdm
        call rdmonthck(bio(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy,k,1,idic),
     &                 936,1)
      enddo   
      close(936)
      bio(:,:,:,1,idic)=bio(:,:,:,1,idic)*REDF(6)
      bio(:,:,:,2,idic)=bio(:,:,:,1,idic)
      where(abs(bio)>10e15) bio=0.0
      
!c----initialize alkalinity and DIC based on salinity relationships
!      do j=1-margin,jj+margin
!      do i=1-margin,ii+margin
!      do k=1,kdm
!      do nk=1,2!
!
!      bio(i,j,k,nk,idic)=(46.5*saln(i,j,k,nk)+546.4)*12.01 !DIC from Shadwick et al. 2011 Marine Chemistry
!c      if (saln(i,j,k,nk).gt.34.65)then
! 
!c          TA = 520.1 + 51.24*saln(i,j,k,nk) ! Atlantic (Millero 1998)
!      bio(i,j,k,nk,ialk)=66.96*saln(i,j,k,nk)-36.803    ! North Atlantic values from Bellerby et al . 2005
!c      else
!c      bio(i,j,k,nk,ialk)=25.3*saln(i,j,k,nk)-1470.    ! North Sea values from Hjalmarson et al. 2010
!c      endif      
!      enddo
!      enddo
!      enddo
!      enddo
!
#ifdef ECOCCO
      bio(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy,1:kdm,1:2,icocc)=0.1
      bio(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy,1:kdm,1:2,ichlc)=0.1/50.
      bio(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy,1:kdm,1:2,icaco3)=0.01
#endif ECOCCO
#endif ECO2
      biodt_sou=baclin!AS*nphys_sou


! initialize river loads:
      call ECOSM_readnutrloads()

      if (mod(86400,int(biodt_sou)) .ne. 0) then
         print *, 'mod_necessary_ecovars:'
         print *, 'baclin*nphys_sou should add up to 60*60*24=86400 s'
         print *, 'please redefine nphys_sou to ensure this...'
         stop
      else
         if (mnproc==1) print *, 'biosteps ok...'
      endif

      end subroutine initialize_ecosys
      end module m_ECOSM_initialize_ecosys