MODULE debmodel
#if defined key_ibm_species
   !!======================================================================
   !!                   ***  MODULE DEB  ***
   !!             DEB anchois/sardine (larves et adultes)
   !!             Martin 09/2011
   !!======================================================================

#include "toolcpp.h"

   !! * Modules used
   USE parameters,  ONLY : kmax,ljmax,limaxp2,ljmaxp2,liminm2,ljminm2,valmanq,rsh,rlg,lchain
   USE_MPI toolmpi, ONLY : iproc, nprocs

   IMPLICIT NONE
   PRIVATE

   !! * Accessibility
   PUBLIC :: deb_init, deb_cycle, deb_init_particle, readfood3d, readtemp3d   ! routines called by ibm

   !! * Shared module variables
   LOGICAL,                PUBLIC :: debuse, F_Fix
   REAL(kind=rsh),         PUBLIC :: ffix
   CHARACTER(LEN=lchain),  PUBLIC :: file_food
   CHARACTER(LEN=lchain),  PUBLIC :: file_NBSS
   REAL(KIND=rsh), ALLOCATABLE, DIMENSION(:,:,:),PUBLIC :: climatemp

   !! * Private variables
   REAL(kind=rsh), ALLOCATABLE, DIMENSION(:,:)   :: slope, intercept
   REAL(kind=rsh), ALLOCATABLE, DIMENSION(:,:,:) :: biomassezoo
   REAL(KIND=rsh), ALLOCATABLE, DIMENSION(:,:,:) :: food1,food2
   REAL(KIND=rsh), ALLOCATABLE, DIMENSION(:,:,:) :: temp1,temp2
   REAL(kind=rsh), PARAMETER                     :: pi=3.14159265358979
   INTEGER                                       :: ilecmemfood,ilecmemtemp

   !! Set species parameters
   REAL(kind=rsh), parameter :: pAm = 780.0_rsh        ! J/d.cm2, spec max assim rate {J_EAm}
   REAL(kind=rsh), parameter :: pMi= 160.0_rsh         ! J/cm3/d, volume-spec maintenance costs
   REAL(kind=rsh), parameter :: EG=2700               ! J/cm3, volume-spec costs of structure
   REAL(kind=rsh), parameter :: vc = 0.25_rsh !      ! cm/j, energy conductance
   REAL(kind=rsh), parameter :: kap = 0.77_rsh         ! -, Fraction to growth + maintenance
   !#  REAL(kind=rsh), parameter :: kM = 0.015_rsh    ! 1/d, somatic maintenance rate coeff
   REAL(kind=rsh), parameter :: kJ = 0.048             ! 1/d, maturity maintenance rate coeff
   !  REAL(kind=rsh), parameter :: Hb=0.004_rsh       ! J, maturity at birth (first feeding)
   !  REAL(kind=rsh), parameter :: Hj = 1000.0_rsh    ! J, maturity at metamorphosis
   !  REAL(kind=rsh), parameter :: Hp =8000.0_rsh     ! J, maturity at puberty
   REAL(kind=rsh), parameter :: Hb=15.0_rsh           ! J, maturity at birth (first feeding)(Nisbet 2012)
   REAL(kind=rsh), public,parameter :: Hj = 30.0_rsh         ! J, maturity at metamorphosis
   REAL(kind=rsh), public,parameter :: Hp =9000.0_rsh        ! J, maturity at puberty (Nisbet 2012)
   REAL(kind=rsh), parameter :: Kx = 0.8_rsh          ! Assimilation efficiency

   !Auxillary and compound parameters
   REAL(kind=rsh), parameter :: TA = 9800.0_rsh       ! K, Arrhenius temperature ;  eggs devlpt -> 9800 = after Regner 1996; larval growth rate : 7500
   REAL(kind=rsh), parameter :: T1 = 293.0_rsh        ! K, Reference temperature ;
  ! REAL(kind=rsh), parameter :: TAL = 30000.0_rsh     ! K, Arrhenius temperature at low temp {Kooy2000}
  ! REAL(kind=rsh), parameter :: TL = 281.0_rsh        ! K, lower boundary temp range
   REAL(kind=rsh), parameter :: shape = 0.188_rsh     ! -, shape coefficient
   ! REAL(kind=rsh), parameter :: d_V  = 1.0_rsh      ! g/cm^3, specific density of structure (Ww)
   REAL(kind=rsh), parameter :: d_V  = 0.22_rsh       ! g/cm^3, specific density of structure (Dw)
   REAL(kind=rsh), parameter :: rho_V = 19000_rsh     !  J/g, = mu_V / w_V; Energy density of structure (DW)
   REAL(kind=rsh), parameter :: rho_E = 37300_rsh     ! J/g, = mu_E / w_E; Energy density of reserve (DW)
   ! REAL(kind=rsh), parameter :: rho_V = 2400.0_rsh  !  J/g, = mu_V / w_V; Energy density of structure (WW)
   ! REAL(kind=rsh), parameter :: rho_E = 11260.0_rsh ! J/g, = mu_E / w_E; Energy density of reserve (WW)
   !#  REAL(kind=rsh), parameter :: Lm=kap*pAm/pMi    ! cm, Maximum volumetric length
   !#  REAL(kind=rsh), parameter :: g=vc*EG/(kap*pAm) ! Energy investment ratio
   REAL(kind=rsh), parameter :: K = 3.0_rsh           ! mgC.m3, half saturation constant
   !#  REAL(kind=rsh), parameter :: Em=pAm/vc         ! J/cm3, max Reserve Energy (max. storage density);
   REAL(kind=rsh), parameter :: Kr=0.95_rsh           ! Fraction of the Gamete buffer fixed into eggs

   ! Ajout Martin
   !  REAL(kind=rsh), parameter :: sh_max =0.75_rsh   !  fraction, max shrinking rate from L_max before dying (Augustine)
   REAL(kind=rsh), parameter :: shapeb = 0.076_rsh    ! -, shape coefficient at birth (early larvae after correction for shrinking)
   REAL(kind=rsh), parameter :: Sizeb = 0.4_rsh       ! -cm, Length at first feeding (birth in DEB)
   REAL(kind=rsh), parameter :: lfactor = 0.23_rsh     ! fraction capacity for larvae, weight of f


   ! Reproduction
   REAL(kind=rsh), parameter :: TR = 286.0_rsh        ! 6, K, temperature at start spawning
   REAL(kind=rsh), parameter :: Eb = 82.5_rsh         ! , J/cm3, specific batch energy (Pethybridge, 0.15(E0)*550 (egg/cm3)
   REAL(kind=rsh), parameter :: E0 = 0.15_rsh           ! J, Energy for one egg, after Heidi (for initialisation only, otherwise dynamic)

CONTAINS

   !!======================================================================

   SUBROUTINE deb_init(restart)
   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE deb_init  ***
   !&E
   !&E ** Purpose : initialize deb values for state variables
   !&E
   !&E ** Description :
   !&E
   !&E ** Called by : ibm
   !&E
   !&E ** External calls :
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E       !  M. Huret (02-2011)
   !&E
   !&E---------------------------------------------------------------------
   !! * Modules used
   USE comvars2d, ONLY : t
   USE ionc4,     ONLY : ionc4_openr,ionc4_read_trajt,ionc4_read_traj,ionc4_close
   USE comtraj,   ONLY : type_patch, patches

   !! * Arguments
   LOGICAL,intent(IN)  :: restart

   !! * Local declarations
   CHARACTER(LEN=lchain)            :: file_inp
   INTEGER                          :: m, n, num, nb_part_nc, is, ie
   TYPE(type_patch), POINTER        :: patch

   INTEGER,        ALLOCATABLE, DIMENSION(:)    :: daysp_nc, yearsp_nc, season_nc
   REAL(KIND=rsh), ALLOCATABLE, DIMENSION(:)    :: edeb_nc, hdeb_nc, rdeb_nc, gam_nc, neggs_nc

   !!----------------------------------------------------------------------
   !! * Executable part


   ! Initialisation of embryo stage
   !!----------------------------------------------------------------------

   IF (.not.F_Fix) THEN
      CALL read_NBSS
   ENDIF

   patch => patches % first
   DO n = 1,patches%nb

      ! Define default values
      CALL deb_init_particle( patch % init_particle )

      IF (restart) THEN

         ! Initialize array with content of NetCDF input file
         file_inp = trim(patch % file_inp)
         CALL ionc4_openr(trim(file_inp),.false.)

         nb_part_nc = patch % nb_part_total
         ALLOCATE( edeb_nc(nb_part_nc), hdeb_nc(nb_part_nc), rdeb_nc(nb_part_nc), gam_nc(nb_part_nc) )
         ALLOCATE( neggs_nc(nb_part_nc), daysp_nc(nb_part_nc), yearsp_nc(nb_part_nc), season_nc(nb_part_nc) )

         CALL ionc4_openr(trim(file_inp), .false.)
         CALL ionc4_read_trajt(file_inp, "EDEB",  edeb_nc,  1, nb_part_nc, 1)
         CALL ionc4_read_trajt(file_inp, "HDEB",  hdeb_nc,  1, nb_part_nc, 1)
         CALL ionc4_read_trajt(file_inp, "RDEB",  rdeb_nc,  1, nb_part_nc, 1)
         CALL ionc4_read_trajt(file_inp, "GAM",   gam_nc,   1, nb_part_nc, 1)
         CALL ionc4_read_trajt(file_inp, "NEGGS", neggs_nc, 1, nb_part_nc, 1)

         CALL ionc4_read_traj (file_inp, "DAYSPAWN",  daysp_nc,  1, nb_part_nc)
         CALL ionc4_read_traj (file_inp, "YEARSPAWN", yearsp_nc, 1, nb_part_nc)
         CALL ionc4_read_traj (file_inp, "SEASON",    season_nc, 1, nb_part_nc)
         CALL ionc4_close(file_inp)

         DO m = 1,patch % nb_part_alloc
            IF ( .NOT. patch % particles(m) % active ) CYCLE
            
            num = patch % particles(m) % num

            patch % particles(m) % E         = edeb_nc(num)
            patch % particles(m) % H         = hdeb_nc(num)
            patch % particles(m) % R         = rdeb_nc(num)
            patch % particles(m) % L         = shape*patch % particles(m) % size
            patch % particles(m) % Gam       = gam_nc(num)
            patch % particles(m) % Neggs     = neggs_nc(num)
            patch % particles(m) % dayspawn  = daysp_nc(num)
            patch % particles(m) % yearspawn = yearsp_nc(num)
            patch % particles(m) % season    = season_nc(num)
            patch % particles(m) % Wdeb      = 0.0_rsh ! a priori pas besoin de recuperer
         END DO

         DEALLOCATE( edeb_nc, hdeb_nc, rdeb_nc, gam_nc, neggs_nc, daysp_nc, yearsp_nc, season_nc )

      ELSE

         ! Initialize array with default values
         patch % particles(:) % R = patch % init_particle % R
         patch % particles(:) % L = patch % init_particle % L
         patch % particles(:) % H = patch % init_particle % H
         patch % particles(:) % E = patch % init_particle % E
         patch % particles(:) % dayspawn  = patch % init_particle % dayspawn
         patch % particles(:) % yearspawn = patch % init_particle % yearspawn

      ENDIF

      patch => patch % next

   END DO  ! loop on patch (n)

   END SUBROUTINE deb_init

!!======================================================================

   SUBROUTINE deb_init_particle(particle)
   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE deb_init_particle  ***
   !&E
   !&E ** Purpose : Initialize DEB part of a given particle with
   !&E              default values.
   !&E
   !&E ** Description :
   !&E
   !&E ** Called by : deb_init, ibm_3d
   !&E
   !&E ** External calls :
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E       !  M. Honnorat (2014-01) Adapt for IBM upgrade
   !&E
   !&E---------------------------------------------------------------------
   !! * Modules used
   USE comvars2d, ONLY : t
   USE comtraj,   ONLY : type_particle

   !! * Arguments
   TYPE(type_particle), INTENT(inout)   :: particle

   !! * Local declarations
   INTEGER              :: jj,mm,aaaa,hh,minu,sec
   CHARACTER(len=19)    :: tool_sectodat

   !!----------------------------------------------------------------------
   !! * Executable part
   
   CALL tool_decompdate(tool_sectodat(t),jj,mm,aaaa,hh,minu,sec)

   ! to start for larvae
   particle % R = 0.0_rsh
   particle % L = Sizeb*shapeb
   particle % H = Hb
!   particle % E = E0
   particle % E = 0.02_rsh
   particle % yearspawn = aaaa+1

   ! to start for adults
!   particle % R = 80000.0_rsh
!   particle % L = particle % size*shape
!   particle % H = Hp
!   particle % E = 100000.0_rsh
!   particle % yearspawn = aaaa

!   particle % L_max = particle % L  ! sauvegarde de L_max atteint au cours du temps

   ! Reproduction
!   particle % Gam = 0.0_rsh
!   particle % Neggs = 0.0_rsh
   particle % dayspawn = 500 ! pour etre sur d etre superieur a jjulien
   particle % dayjuv   = 0

   !Weight
!   particle % Wdeb = 0.0_rsh

   END SUBROUTINE deb_init_particle

!!======================================================================

   SUBROUTINE deb_cycle(particle,dt,year,month,day)
   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE deb_cycle  ***
   !&E
   !&E ** Purpose :  Modele D.E.B. generique adapte anchois
   !&E ** Calculate variation in DEB compartments with ODE
   !&E
   !&E ** Description :
   !&E
   !&E ** Called by : ibm
   !&E
   !&E ** External calls :
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E       !  M. Huret (02-2011) from DEB-Larve L. Pecquerie
   !&E
   !&E---------------------------------------------------------------------
   !! * Modules used
   USE comvars2d, ONLY : jjulien
   USE comtraj,   ONLY : type_particle

   !! * Arguments
   TYPE(type_particle), INTENT(inout)   :: particle
   REAL(KIND=rlg),      INTENT(in)      :: dt
   INTEGER,             INTENT(in)      :: year, month, day

   !! * Local declarations
   REAL(kind=rsh)   :: cor_T, K_food,cor_l,shape_fun
   REAL(kind=rsh)   :: pA,pC,pG,pJ,pR,pMT,pM,pM2,pM3,pR2,pGam
   REAL(kind=rsh)   :: kJT,pAmT,vT,f!,kMT
   REAL(kind=rsh)   :: E,L,R,H,WV,WE,WR,WG,NRJ_V,NRJ_g,Ebatch,Ebatchlast,Gam
   REAL(kind=rsh)   :: T,X,WWoDW
   REAL(kind=rsh)   :: dR,dH,dE,dL,dGam

   !!----------------------------------------------------------------------
   !! * Executable part

   ! Remise a 0 pour flux non ecrases a chaque pas de temps
   pR2 = 0.0_rsh
   pM2 = 0.0_rsh
   pM3 = 0.0_rsh
   pGam = 0.0_rsh
   Ebatch = 0.0_rsh
   Ebatchlast = 0.0_rsh

   ! Affectation ------------------------------------------------------------------
   E = particle % E
   L = particle % L
   R = particle % R
   H = particle % H
   Gam = particle % Gam

   ! Temperature ------------------------------------------------------------------
   T = particle % temp + 273.0_rsh
   cor_T = exp(TA/ T1 - TA/ T) !* (1 + exp(TAL/ T1 - TAL/ TL))/ (1.0_rsh + exp(TAL/ T - TAL/ TL))
        ! more elaborated correction for lower temp limit

   ! Food -------------------------------------------------------------------------
   X = 0.0_rsh
   K_food = K !  + L * 10; % account for changes in preference
   f = ffix

   IF ( .NOT. F_Fix ) THEN
      IF (particle % stage < 5) THEN
         X = get_Xdeb(particle%xpos, particle%ypos, particle%spos, particle%size, 0.0_rsh,0.0_rsh)
      ELSE
         X = get_Xdeb(particle%xpos, particle%ypos, particle%spos, particle%size, 0.0_rsh,40.0_rsh)
      ENDIF
      f = X / (X + K_food) !    functional response
   ENDIF

   !-------------------------------------------------------------------------------
   ! shape correction function
   !cor_shape = (MAX( Hb, MIN(H, Hj))/ Hb)**(1.0_rsh/3.0_rsh)  ! cf p161 Laure ? pour juveniles et adultes Hj=Hb ?
   cor_l = MIN(1.0_rsh,MAX(((H-Hb)+(Hj-H)*lfactor) / (Hj-Hb),lfactor))
   f = f*cor_l

   ! CORRECTED PARAMETERS --------------------------------------
   pAmT = pAm*cor_T !*cor_shape , sous hypothese V1-morph
   vT  = vc*cor_T !*cor_shape , sous hypothese V1-morph
   pMT = pMi*cor_T
   kJT = kJ*cor_T
   !kMT = pMT/EG

   ! FLUXES------------------------------------------------------------------------
   pA = pAmT*f*L**2.0_rsh               ! Assimilation
   pM = pMT*L**3.0_rsh                    ! J/d somatic maintenance;

   pC = (E/L**3.0_rsh)*((vT*EG*L**2.0_rsh)+pM) &
        / (EG+(kap*E/L**3.0_rsh))        ! J/d, Reserve Mobilisation, ref=Kooijman (Book:p37); (ok avec eq Laure 3.12)
   pG = (kap*pC)-pM                     ! J/d, Energy allocated to somatic growth;
   pJ = kJT*H                           ! J/d, allocation to maturity maintenance;
   pR = (1-kap)*pC-pJ                   ! J/d, allocation to maturity or reproduction;

   ! Reproduction and  Emergency maintenance
   IF ( H >= Hp ) THEN
      IF (particle % season==0 .and. year>=particle % yearspawn .and. T>TR .and. month>3 .and. month<9) THEN
         particle % season = 1
      ENDIF

      pR2 = R*(vT/L + pMT/EG)*(1.0_rsh-kap*E/(EG*L**3.0_rsh+kap*E))  ! J/d, potential allocation to gametogenesis (Pethybridge13 or Bernard11)

      IF ( pG < 0.0_rsh ) THEN    ! Bad condition, need to remobilise energy
         pM2 = min(-pG,pR2)       ! J/d, first mobilization from R (same speed as for gametogenesis)
         IF ( particle % season==1 .and. Gam>0.0_rsh )  &
            pM3 = min(Kr*Gam*86400.0_rsh/dt,-pG-pM2)
      ENDIF
      IF ( particle % season==1 )  pGam=Kr*(pR2-pM2)
   ENDIF

   !ODE---------------------------------------------------------------------
   dE = pA-pC                            ! J/d, Dynamics of Reserve Energy;
   dL = 0.0_rsh
   IF (pG>=0.0_rsh) THEN
      dL = pG/(3.0_rsh*EG*L**2.0_rsh)    ! cm/d, growth increment;
   ELSE                                  ! Test for mortality
      IF ((pG+pM2+pM3)<0.0_rsh) THEN
         particle % super = 0            ! maintenance pas assuree then mortality
         particle % flag  = -valmanq     ! to stop transport of particles when dead
         L = 0.0_rsh ! pour optimisation sinon on garde une valeur 'normale' et fcout pas affectee
         R = 0.0_rsh
         E = 0.0_rsh
      ENDIF
   ENDIF

   IF (H<Hp) THEN
      dH = pR                   ! J/d, variation in level of maturity;
      dR = 0.0_rsh              ! J/d, variation in reproduction buffer energy;
   ELSE
      dH = 0.0_rsh
      IF ( particle % season == 1 ) THEN
         IF (Gam>(2.0_rsh*Eb*L**3.0_rsh)) THEN        ! enough energy to spawn 1 batch and same amount remains for vitellogenesis (Somarakis Reproduce Del.), explicit
            Ebatch = Eb*L**3.0_rsh
            particle % Neggs = particle % Neggs+Ebatch/E0
            !test if first batch...
            IF (year>=particle % yearspawn) THEN
               particle % yearspawn=year+1
               particle % dayspawn=jjulien
            ENDIF
         ENDIF
         dGam = pGam-(2-Kr)*pM3
         Gam  = Gam+dGam*dt/86400.0_rsh-Ebatch
         IF ((jjulien>particle % dayspawn+100 .or. month>8)) THEN     ! last batch, implicit
            Ebatchlast = Gam
            particle % dayspawn = 500
            particle % season = 0
         ENDIF
         dR = pR-pR2
         particle % Gam = Gam-Ebatchlast ! on vide Gonades
      ELSE
         dR=pR-pM2 ! only necessary fraction is mobilized for emergency
      ENDIF
   ENDIF

   particle % E = E + dE*dt/86400.0_rsh
   particle % L = L + dL*dt/86400.0_rsh
   particle % H = H + dH*dt/86400.0_rsh
   particle % R = R + dR*dt/86400.0_rsh+Ebatchlast ! on recupere energie dernier batch dans R

   !---------------------------------------------------------
   ! Calcul du shape si variable en fonction du stade de developement
   ! transition entre shapeb et shape entre 2 et 4 cm
   IF ( particle % L <  2.0_rsh*shapeb ) shape_fun = shapeb
   IF ( particle % L >= 2.0_rsh*shapeb .and. particle % L <= 4.0_rsh*shape)                                &
                                         shape_fun = ((L-2.0_rsh*shapeb)*shape+(4.0_rsh*shape-L)*shapeb)   &
                                                      / (4.0_rsh*shape-2.0_rsh*shapeb)
   IF ( particle % L >  4.0_rsh*shape  ) shape_fun = shape

 ! Physical length
   particle % size = particle % L/shape_fun !  in cm

 ! Weight
   WV = d_V*particle % L**3.0_rsh
   WE = particle % E/rho_E ! W of E
   WR = particle % R/rho_E !  W of E_R
   WG = particle % Gam/rho_E !  W of Gam
   particle % Wdeb = WV+WE+WR+WG
!   particle % WV = WV
!   particle % WE = WE
!   particle % WR = WR
!   particle % WG = WG

 ! Energy
   NRJ_V = rho_V*WV
   NRJ_g = (NRJ_V + particle % E + particle % R+particle % Gam)
   particle % NRJ = NRJ_g/(particle % Wdeb)/1000.0_rsh ! Energie densite (dry)

   WWoDW = 7.84-0.18*particle % NRJ
   particle % Wdeb = particle % Wdeb*WWoDW ! conversion Dry to Wet
   particle % NRJ = NRJ_g/(particle % Wdeb)/1000.0_rsh  ! Energie densite (wet)

! Auxillary save
   particle % X = X

   END SUBROUTINE deb_cycle

!!======================================================================

   SUBROUTINE read_NBSS

   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE read_NBSS  ***
   !&E
   !&E ** Purpose : Read NBSS (spatially explicit, but not temporally) and
   !&E              provide 2D slope and intercept
   !&E
   !&E ** Description :
   !&E
   !&E ** Called by :
   !&E
   !&E ** External calls :ionc4_openr, ionc4_read_time, ionc4_read
   !&E                     ionc4_read_subxyt, ionc4_read_subzxyt
   !&E                     ionc4_read_subxyt, ionc4_read_subzxyt
   !&E                     ionc4_close
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E       !  2012-03 (M. Huret) Original code
   !&E---------------------------------------------------------------------
   !! * Modules used
   USE ionc4, ONLY : ionc4_openr,ionc4_read_subxy,ionc4_close

   !! * Arguments

   !! * Local declarations
   CHARACTER(LEN=lchain)    :: name_in_slope
   CHARACTER(LEN=lchain)    :: name_in_inter

   !!----------------------------------------------------------------------
   !! * Executable part

   name_in_slope = 'slope'
   name_in_inter = 'intercept'

   CALL ionc4_openr(file_NBSS,l_out_nc4par=.true.)

   ALLOCATE(slope(liminm2:limaxp2,ljminm2:ljmaxp2),intercept(liminm2:limaxp2,ljminm2:ljmaxp2))

   CALL ionc4_read_subxy(file_NBSS,TRIM(name_in_slope),slope,liminm2,limaxp2,ljminm2,ljmaxp2,0,0)
   CALL ionc4_read_subxy(file_NBSS,TRIM(name_in_inter),intercept,liminm2,limaxp2,ljminm2,ljmaxp2,0,0)

   CALL ionc4_close(file_NBSS)

   END SUBROUTINE read_NBSS

!!======================================================================

   SUBROUTINE readfood3d(jul)

   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE readfood3d  ***
   !&E
   !&E ** Purpose : read model variables from a NetCDF file
   !&E
   !&E ** Description : read micro and zooplankton from a MARS 3D file, and sum those
   !&E
   !&E ** Called by :
   !&E
   !&E ** External calls : ionc4_openr, ionc4_read_time, ionc4_read
   !&E                     ionc4_read_subxyt, ionc4_read_subzxyt
   !&E                     ionc4_read_subxyt, ionc4_read_subzxyt
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E       !           M. Huret
   !&E---------------------------------------------------------------------

   !! * Modules used
   USE comvars2d, ONLY : nbouc
   USE ionc4,     ONLY : ionc4_openr,ionc4_read_subzxyt

   !! * Arguments
   INTEGER, INTENT(inout)                               :: jul

   !! * Local declarations
   INTEGER                                              :: ilec,ilecp1
   CHARACTER(LEN=lchain)                                :: name1_in_food,name2_in_food
   REAL(KIND=rsh), ALLOCATABLE, DIMENSION(:,:,:)        :: food1_1,food2_1,food1_2,food2_2
   REAL(kind=rlg)                                       :: dt1,dt2

   !!----------------------------------------------------------------------
   !! * Executable part

   jul = jul+30   ! patch
   IF ( jul > 365 ) jul = jul-365

   name1_in_food = 'microzooplankton_nitrogen'
   name2_in_food = 'mesozooplankton_nitrogen'

   CALL ionc4_openr(file_food,l_out_nc4par=.true.)

   IF( nbouc == 1 ) THEN
      ALLOCATE(biomassezoo(1:kmax,liminm2:limaxp2,ljminm2:ljmaxp2))
      ALLOCATE(food1(1:kmax,liminm2:limaxp2,ljminm2:ljmaxp2),food2(1:kmax,liminm2:limaxp2,ljminm2:ljmaxp2))
   ENDIF

   ilec = INT(jul/7)   !52 dates fournies ~ milieu de chaque semaine, on suppose premiere date le 7
   IF ( ilec == 0 ) ilec = 52
   ilecp1 = ilec+1
   IF ( ilecp1 > 52 ) ilecp1 = 1

   dt1=(real(ilecp1,rlg)*7.0_rlg-real(jul,rlg))/7.0_rlg
   IF (dt1 < 0) dt1=0.0_rlg
   dt2=1.0_rlg-dt1

   ! lecture des tableaux deux indices qui encadrent le pas de temps courant
   IF( nbouc == 1 ) THEN
      ALLOCATE(food1_1(1:kmax,liminm2:limaxp2,ljminm2:ljmaxp2),food1_2(1:kmax,liminm2:limaxp2,ljminm2:ljmaxp2))
      ALLOCATE(food2_1(1:kmax,liminm2:limaxp2,ljminm2:ljmaxp2),food2_2(1:kmax,liminm2:limaxp2,ljminm2:ljmaxp2))

      CALL ionc4_read_subzxyt(file_food,TRIM(name1_in_food),food1_1,liminm2,limaxp2,ljminm2,ljmaxp2,1,kmax,ilec,0,0,1)
      CALL ionc4_read_subzxyt(file_food,TRIM(name2_in_food),food1_2,liminm2,limaxp2,ljminm2,ljmaxp2,1,kmax,ilec,0,0,1)
      CALL ionc4_read_subzxyt(file_food,TRIM(name1_in_food),food2_1,liminm2,limaxp2,ljminm2,ljmaxp2,1,kmax,ilecp1,0,0,1)
      CALL ionc4_read_subzxyt(file_food,TRIM(name2_in_food),food2_2,liminm2,limaxp2,ljminm2,ljmaxp2,1,kmax,ilecp1,0,0,1)

      ! sum micro and mesozoo from model
      food1(:,:,:) = food1_1(:,:,:) + food1_2(:,:,:)
      food2(:,:,:) = food2_1(:,:,:) + food2_2(:,:,:)
      DEALLOCATE(food1_1,food1_2,food2_1,food2_2)
      ilecmemfood=ilec

   ELSE

      IF(ilecmemfood/=ilec) THEN
         ALLOCATE(food2_1(1:kmax,liminm2:limaxp2,ljminm2:ljmaxp2),food2_2(1:kmax,liminm2:limaxp2,ljminm2:ljmaxp2))

         CALL ionc4_read_subzxyt(file_food,TRIM(name1_in_food),food2_1,liminm2,limaxp2,ljminm2,ljmaxp2,1,kmax,ilecp1,0,0,1)
         CALL ionc4_read_subzxyt(file_food,TRIM(name2_in_food),food2_2,liminm2,limaxp2,ljminm2,ljmaxp2,1,kmax,ilecp1,0,0,1)

         food1(:,:,:) = food2(:,:,:)
         food2(:,:,:) = food2_1(:,:,:) + food2_2(:,:,:)
         DEALLOCATE(food2_1,food2_2)
         ilecmemfood=ilec
      ENDIF

   ENDIF

   biomassezoo(:,:,:) = food1(:,:,:)*dt1 + food2(:,:,:)*dt2

   ! convert in mgC/m3
   biomassezoo(:,:,:) = biomassezoo(:,:,:)*5.45_rsh*12.0_rsh

   END SUBROUTINE readfood3d

!!======================================================================

   SUBROUTINE readtemp3d(jul)

   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE readtemp3d  ***
   !&E
   !&E ** Purpose : read model variables from a NetCDF file
   !&E
   !&E ** Description : read temp from a climatology file
   !&E
   !&E ** Called by :
   !&E
   !&E ** External calls : ionc4_openr, ionc4_read_time, ionc4_read
   !&E                     ionc4_read_subxyt, ionc4_read_subzxyt
   !&E                     ionc4_read_subxyt, ionc4_read_subzxyt
   !&E
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E       !           M. Huret
   !&E
   !&E---------------------------------------------------------------------

   !! * Modules used
   USE comvars2d, ONLY : nbouc
   USE ionc4,     ONLY : ionc4_openr,ionc4_read_subzxyt

   !! * Arguments
   INTEGER, INTENT(inout)                               :: jul

   !! * Local declarations
   INTEGER                                              :: ilec,ilecp1
   CHARACTER(LEN=lchain)                                :: name_in_temp
   REAL(KIND=rsh), ALLOCATABLE, DIMENSION(:,:,:)        :: temp_1,temp_2
   REAL(kind=rlg)                                       :: dt1,dt2
   CHARACTER(LEN=lchain)                                :: file_temp

   !!----------------------------------------------------------------------
   !! * Executable part

   file_temp = '../../inputs/Climatotemp_06_13.nc'
   name_in_temp = 'TEMP'

   CALL ionc4_openr(file_temp,l_out_nc4par=.true.)

   IF( nbouc == 1 ) THEN
      ALLOCATE(climatemp(1:kmax,liminm2:limaxp2,ljminm2:ljmaxp2),   &
                   temp1(1:kmax,liminm2:limaxp2,ljminm2:ljmaxp2),   &
                   temp2(1:kmax,liminm2:limaxp2,ljminm2:ljmaxp2))
   ENDIF

   ilec=INT(jul)   !365 dates fournies
   IF (ilec==0) ilec=365
   ilecp1=ilec+1
   IF (ilecp1 > 365) ilecp1=1

   dt1=ilecp1-jul
   IF (dt1 < 0) dt1=0.0_rlg
   dt2=1.0_rlg-dt1
!   write(*,*) 'dt', dt1,dt2,ilec,ilecp1,jul

     ! lecture des tableaux deux indices qui encadrent le pas de temps courant
   IF( nbouc == 1 ) THEN
      ALLOCATE(temp_1(1:kmax,liminm2:limaxp2,ljminm2:ljmaxp2),temp_2(1:kmax,liminm2:limaxp2,ljminm2:ljmaxp2))
      CALL ionc4_read_subzxyt(file_temp,TRIM(name_in_temp),temp_1,liminm2,limaxp2,ljminm2,ljmaxp2,1,kmax,ilec,0,0,1)
      CALL ionc4_read_subzxyt(file_temp,TRIM(name_in_temp),temp_2,liminm2,limaxp2,ljminm2,ljmaxp2,1,kmax,ilecp1,0,0,1)
      temp1(:,:,:) = temp_1(:,:,:)
      temp2(:,:,:) = temp_2(:,:,:)
      DEALLOCATE(temp_1,temp_2)
      ilecmemtemp=ilec

   ELSE

      IF(ilecmemtemp/=ilec) THEN
         ALLOCATE(temp_2(1:kmax,liminm2:limaxp2,ljminm2:ljmaxp2))

         CALL ionc4_read_subzxyt(file_temp,TRIM(name_in_temp),temp_2,liminm2,limaxp2,ljminm2,ljmaxp2,1,kmax,ilecp1,0,0,1)

         temp1(:,:,:) = temp2(:,:,:)
         temp2(:,:,:) = temp_2(:,:,:)

         DEALLOCATE(temp_2)
         ilecmemtemp=ilec
      ENDIF

   ENDIF

   climatemp(:,:,:) = temp1(:,:,:)*dt1 + temp2(:,:,:)*dt2

   END SUBROUTINE readtemp3d

!!======================================================================

   FUNCTION get_Xdeb(xpos,ypos,spos,size,minprof,maxprof)

   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE  get_Xdeb ***
   !&E
   !&E ** Purpose : calculate f (standardized food response) for the DEB
   !&E
   !&E ** Description :
   !&E
   !&E ** Called by :
   !&E
   !&E ** History :
   !&E       !  M. Huret (07-2012)
   !&E---------------------------------------------------------------------
   !! * Modules used
   USE ibmtools, ONLY : ibm_profmean
   USE trajectools, ONLY : ksupkinf
   USE comvars3d, ONLY : sigw
   USE comvarp2d,  ONLY : ssh
   USE comvars2d,  ONLY : h0_g,jjulien
#if defined key_siggen
   USE comsiggen, ONLY : hc_sig
#endif

   !! * Arguments
   REAL(kind=rsh), INTENT (in)  :: xpos,ypos,spos,size,minprof,maxprof
   REAL(kind=rsh) :: get_Xdeb

   !! * Local declarations
   REAL(kind=rsh) ::  slope_max_preywidth,slope_min_preywidth,maxwidth,esdmax,esdmin,minwidth
   REAL(kind=rsh) ::  biom,sl,inter
   INTEGER        ::  kwm,kw,i,j

   !!----------------------------------------------------------------------
   !! * Executable part

   slope_max_preywidth = 0.035_rsh    ! Delivrable FACTS
   slope_min_preywidth = 0.005        ! Test
   maxwidth = slope_max_preywidth*size*10.0_rsh ! cm -> mm
   maxwidth = min(maxwidth,3.0_rsh) ! arbitraire ! pour coherence a maxzoowidth issu du modele ci-dessous
   esdmax = size2ESD(maxwidth)
   minwidth = 0.04_rsh+slope_min_preywidth*(size*10.0_rsh-4.0_rsh)
   minwidth = min(minwidth,0.2_rsh) !toute gamme mesozoo
   esdmin = size2ESD(minwidth)

   sl = slope(NINT(xpos),NINT(ypos)) ! slope climato LOPC
   inter = intercept(NINT(xpos),NINT(ypos))

   IF ( maxprof == 0.0_rsh )  THEN  ! larva
      kwm = 0
      kw  = kmax
      CALL ksupkinf(spos,sigw,kmax+1,kw,kwm,12)
      biom = biomassezoo(kw,NINT(xpos),NINT(ypos))
   ELSE
      i = NINT(xpos)
      j = NINT(ypos)
      IF ( h0_g(i,j) > 0.0_rsh ) THEN
         biom = ibm_profmean(biomassezoo(:,i,j),minprof,maxprof,h0_g(i,j),ssh(i,j) &
#if defined key_siggen
              ,hc_sig(i,j) &
#endif
           )
      ELSE
         biom = 2.0_rsh
      ENDIF
   ENDIF

   IF( biom <= 0.0_rsh .or. biom >= 999.0_rsh) THEN ! not exactly same grid, so lacking cells along the coast
      sl   = -1.0_rsh   ! slope=-1, default
      biom = 2.0_rsh  ! 2mg/m3 = mini Poulet1996
   ENDIF

   biom = max(2.0_rsh,biom) ! pour eviter valeurs trop basses
   IF ( jjulien > 90 .and. jjulien < 300 ) THEN
      biom = max(10.0_rsh,biom) ! pour eviter valeurs trop basses, mini observe ete
   ENDIF

   inter = biom_sl2inter(sl,biom) ! on recalcule l'intercept du NBSS a partir de biomasse modele
   biom  = NBSS_2biom(sl,inter,esdmin,esdmax) ! on recalcule l'abondance dispo en fonction de la taille

   get_Xdeb=biom

   END FUNCTION get_Xdeb

!!======================================================================

   FUNCTION biom_sl2inter(sl,biom)
   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE  biom_sl2inter ***
   !&E
   !&E ** Purpose : calculate the intercept from NBSS slope and model biomass
   !&E
   !&E ** Description :
   !&E
   !&E ** Called by :
   !&E
   !&E ** History :
   !&E       !  M. Huret (07-2012)
   !&E---------------------------------------------------------------------
   !! * Modules used

   !! * Arguments
   REAL(kind=rsh), INTENT (in)  :: sl, biom
   REAL(kind=rsh)   ::  biom_sl2inter

   !! * Local declarations
   REAL(kind=rsh)   :: minesd,maxesd,minwidth,maxwidth,slin,xx1,xx2

   !!----------------------------------------------------------------------
   !! * Executable part

   minwidth = 0.02_rsh  !  min width zoo dans modele en mm
   maxwidth = 3.0_rsh
   minesd = size2ESD(minwidth) ! transfo en ESD
   maxesd = size2ESD(maxwidth)
   slin = sl
   IF (slin==-1) slin= -1.01

   xx1 = ESD2weight(minesd)
   xx2 = ESD2weight(maxesd)

   biom_sl2inter = log((biom*(slin+1))/(xx2**(slin+1)-xx1**(slin+1)))

   END FUNCTION biom_sl2inter

!!======================================================================

   FUNCTION NBSS_2biom(sl,int,size1,size2)
   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE  NBSS_2biom ***
   !&E
   !&E ** Purpose : calculate the biomasse within a size range given a NBSS (slope + intercept)
   !&E
   !&E ** Description :
   !&E
   !&E ** Called by :
   !&E
   !&E ** History :
   !&E       !  M. Huret (07-2012)
   !&E
   !&E---------------------------------------------------------------------
   !! * Modules used

   !! * Arguments
   REAL(kind=rsh), INTENT(in)   :: sl,int,size1,size2
   REAL(kind=rsh)   :: NBSS_2biom

   !! * Local declarations
   REAL(kind=rsh)   :: xmid,xx1,xx2
   REAL(kind=rsh)   :: slin, interin,biom

   !!----------------------------------------------------------------------
   !! * Executable part

   slin=sl
   interin=int
   IF (slin==-1) slin= -1.01

   interin = exp(int)  ! log(y) = sl * log(x) + int     eq.    y = exp(int) * x^sl

   ! convert size1 and size2 from mm esd to mgC:
   xx1  =ESD2weight(size1)
   xx2 = ESD2weight(size2)

   ! compute xmid in mugC (mean so as to have biom[x1,xmid]=biom[xmid,x2]
!   xmid = ((xx1_2[2]^(sl+1) + xx1_2[1]^(sl+1))/2)^(1/(sl+1));

   ! compute biom in mgC.m^{-3} and ab in #.m^{-3}
   biom = (interin*xx2**(slin+1))/(slin+1) - (interin*xx1**(slin+1))/(slin+1)
!   ab = biom / (xmid/1000)

   NBSS_2biom = biom

   END FUNCTION NBSS_2biom

!=========================================================================

   FUNCTION size2ESD(width)
   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE size2esd ***
   !&E
   !&E ** Purpose : Calculate esd from width
   !&E
   !&E ** Description :
   !&E!      ESD=2*sqrt(ab) ou ESD=sqrt(WL), Herman92 avec w et l largeur et longueur. Pour a=3b, ESD=sqrt(3*w^2)=w*sqrt(3)
   !&E ** Called by :
   !&E
   !&E ** History :
   !&E       !  M. Huret (07-2012)
   !&E---------------------------------------------------------------------
   !! * Modules used

   !! * Arguments
   REAL(kind=rsh), INTENT(in)   :: width
   REAL(kind=rsh)  :: size2ESD

   !! * Local declarations

   !!----------------------------------------------------------------------
   !! * Executable part

   size2ESD = width*3.0_rsh**0.5_rsh

   END FUNCTION size2ESD

!=========================================================================

   FUNCTION ESD2weight(esd)
   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE size2esd ***
   !&E
   !&E ** Purpose : calculate Carbon weight from ESD
   !&E
   !&E ** Description :
   !&E!   ! from Lehette & Hernandez-Leon 2009 (DW general mesozooplankton)
   !&E    ! from Mauchline 1998 (average carbon content of copepods)
   !&E ** Called by :
   !&E
   !&E ** History :
   !&E       !  M. Huret (07-2012)
   !&E---------------------------------------------------------------------
   !! * Modules used

   !! * Arguments
   REAL(kind=rsh), INTENT(in)   :: esd
   REAL(kind=rsh)  :: ESD2weight

   !! * Local declarations
   REAL(kind=rsh)  :: S,DW

   !!----------------------------------------------------------------------
   !! * Executable part
 
   S  = pi*(esd/2.0_rsh)**2.0_rsh
   DW = 43.38_rsh*S**1.5_rsh   !from Lehette & Hernandez-Leon 2009 (general mesozooplankton)
   ESD2weight = 0.447_rsh*DW   !from Mauchline 1998 (average carbon content of copepods)

   END FUNCTION ESD2weight

!=========================================================================

#endif 
!end key_ibm_species

END MODULE debmodel
