Merge branch 'master' of gitlasso.uqar.ca:dumoda01/gotm_ismer

nml/bio_gsj.nml

@@ -111,6 +111,7 @@
   g2max   = 1.0
   k3      = 1.0
   beta    = 0.625
+  txc4zoo = .true.
   mu21    = 0.3
   mu21max = 1.5
   bhc     = 11.1

nml/bio_polynow.nml

+#$Id$
+!-------------------------------------------------------------------------------
+! Fasham et al. biological model with modifications by Kuehn and Radach
+! 
+!  numc=        number of compartments for geobiochemical model
+!
+!  p_initial= initial phytoplankton concentration                  [mmol n/m3]
+!  z_initial= initial zooplakton concentration                     [mmol n/m3] 
+!  b_initial= initial bacteria concentration                       [mmol n/m3]
+!  dph_initial= initial dph    concentration                       [mmol n/m3]
+!  dzo_initial= initial dzo    concentration                       [mmol n/m3]
+!  fp_initial= initial fp    concentration                         [mmol n/m3]
+!  msn_initial= initial dzo    concentration                       [mmol n/m3]
+!  n_initial= *** see obs.nml ***                                  [mmol n/m3]
+!  a_initial= *** see obs.nml ***                                  [mmol n/m3]
+!  l_initial= initial LDON concentration                           [mmol n/m3]
+!  p0       = minimum phytoplankton concentration                  [mmol n/m3]
+!  z0       = minimum zooplakton concentration                     [mmol n/m3]
+!  b0       = minimum bacteria concentration                       [mmol n/m3]
+!  mu5      = nitrification      rate                              [1/day] !!GG-CHG1
+!  vp       = maximum phytoplankton uptake rate                    [1/day] 
+!  alpha    = slope of the PI-curvea                               [m2/(W day)] 
+!  inib     = inhibition slope of the PI-curve (positive)          [m2/(W day)] 
+!  kc       = attenuation constant for the self shading effect     [m**2/mmol N]
+!  k1       = half saturation constant nitrate uptake              [mmol n/m3]
+!  k2       = half saturation constant ammonium uptake             [mmol n/m3]
+!  mu1      = phytoplankton mortality rate                         [1/day] 
+!  k5       = half saturation constant phytoplankton mortality     [mmol n/m3]
+!  gamma    = exudation fraction                                   [-]
+!  capt     = Percentage of living phyto. captured in msn	   [%/day] - GG-D 
+! sp_mort   = Phytoplankton non grazing mortality		   [1/day] - GG-D 
+! sp_mort2  = Quadratic mortality/aggrgation rate coefficient	   [1/day] - GG-D 	
+!  w_p      = phytoplankton settling velocity                      [m/day]
+! stip      = Stickiness of choosen Phytoplankton species          [-]     - GG-D
+! stidph    = Stickiness of dead Phytoplankton                     [-]     - GG-D
+!  gmax     = maximum ingestion rate                               [1/day]
+!  k3       = half saturation constant ingestion                   [mmol n/m3]
+!  beta     = grazing efficiency                                   [-]
+!  betamsn  = grazing efficiency on marine snow                    [-]     - GG-D 
+!  mu2      = maximum zooplankton loss rate                        [1/day]
+!  k6       = half saturation zooplankton loss                     [mmol n/m3]
+!  delta    = fractional zooplankton loss to LDON                  [-]
+!  epsi     = fractional zooplankton loss to ammonium              [-]
+!  r1       = grazing preference phytoplankton                     [-]
+!  r2       = grazing preference bacteria                          [-] 
+!  r3       = grazing preference dead phytoplankton (dph)          [-]     - GG-D  
+!  r4       = grazing preference dead zooplankton (dzo)            [-]     - GG-D 
+!  r5       = grazing preference fecal pellets (fp)                [-]     - GG-D  
+!  r6       = grazing preference marine snow (msn)                 [-]     - GG-D 
+!zingest    = Proportion of grazed matter added to zoo. biomass    [-]     - GG-D 
+!  eg       = Percentage of eggestion			           [%/day] - GG-D 
+! stidzo    = Stickiness of choosen zooplankton species            [-]     - GG-D
+! stifp     = Stickiness of fecal pellets         		   [-]     - GG-D
+!  pmin     = intensity rariadtion in the water column
+!  w_zmax   = Max. swimming speed for zoo (not exceed 320)          [1/day]
+!  bertha   = amplitud coefficient
+!  parcrit  =min. phyto concentration needed for zooplankton
+!  vb       = maximum bacterial uptake rate                        [1/day]
+!  k4       = half saturation bacterial uptake                     [mmol n/m3]
+!  mu3      = bacteria excretion rate                              [1/day]  
+!  eta      = uptake ratio ammonium:LDON                           [-]  
+!  mbac     = Percentage of bacerias captured in msn               [%/day] - GG-D 
+!  dphlossl = Dph breakdown rate                                   [%/day] - GG-D
+!  dzolossl = Percentage of Dzo loss in LDON                       [%/day] - GG-D
+!  dzolossb = Percentage of Dzo loss in bacteria respiration       [%/day] - GG-D
+!  fplossl  = Percentage of fp loss in LDON                        [%/day] - GG-D
+!  fplossb  = Percentage of fp loss in bacteria respiration        [%/day] - GG-D
+!  w_d1     =        [1/day]
+!  w_d2     =        [1/day]
+!  w_d3     =        [1/day]
+!  w_d4     =        [1/day]
+! mldon     = Percentage of LDON used to constitute matrice of msn  [%/day] - GG-D
+!lmin       = LDON limit concentration to become the matrix of msn  [mmol n/m3]- GG-D
+!Coag_coef  = Specification by user of utilisation of model calculated Coagulation coefficient(false) or given parameter(true) [true.false] - GG-D
+!betap      = Coagulation coefficient of Phytoplankton             [1/day]- GG-D 
+!betadph    = Coagulation coefficient of dead Phytoplankton        [1/day]- GG-D
+!betadzo    = Coagulation coefficient of deadzooplankton           [1/day]- GG-D 
+!betafp     = Coagulation coefficient of fecal pellets             [1/day]- GG-D 
+!-------------------------------------------------------------------------------
+ &bio_polynow_nml
+  numc=      10
+  p_initial= 0.012
+  z_initial= 0.012
+  b_initial= 0.001
+  dph_initial= 0.012
+  dzo_initial= 0.012
+  fp_initial= 0.012
+  msn_initial= 0.012
+  l_initial= 0.1
+  p0=        0.0001
+  z0=        0.0001
+  b0=        0.0001
+  mu5=       0.1
+  vp=        0.3
+  alpha=     0.04
+  inib=      0.06
+  kc=        0.03
+  k1=        1.0
+  k2=        0.8
+  mu1=       0.05
+  k5=        0.2
+  gamma=     0.0
+  capt=      0.0
+  sp_mort=   0.0
+  sp_mort2=  0.0
+  w_p=       0.5
+ stip=       0.0
+ stidph=     0.0
+  gmax=      0.6
+  k3=        1.0
+  beta=      0.625
+  betamsn=   0.0
+  mu2=       0.3
+  k6=        0.2
+  delta=     0.0
+  epsi=      0.70
+  r1=        0.1
+  r2=        0.2
+  r3=        0.3
+  r4=        0.1
+  r5=        0.2
+  r6=        0.1
+  zingest=   0.0
+  eg=        0.0
+  stidzo=    0.0
+  stifp=     0.0
+  pmin=      0.7
+  w_zmax=    320
+  bertha=    0.05
+  parcrit=   10
+  vb=        0.
+  k4=        0.5
+  mu3=       0.03
+  eta=       0.0
+  mbac=      0.0
+  dphlossl=  0.0
+  dzolossl=  0.0
+  dzolossb=  0.0
+  fplossl=   0.0
+  fplossb=   0.0
+  w_d1=      0     
+  w_d2=      0
+  w_d3=      0 
+  w_d4=      0 
+  mldon=     0.0
+  lmin=      0.0
+  coag_coef= .false.
+  betap=     0.0
+  betadph=   0.0
+  betadzo=   0.0
+  betafp=    0.0
+ /

scripts/shell/run_enp.sh

+#!/bin/bash
+
+export EXPDIR=$(pwd)
+echo $EXPDIR
+
+if [ $# -ne 1 ]
+then
+  echo ""
+  echo "Usage: $0 <exp>"
+  echo ""
+  echo " <exp> should be one of the following:"
+  echo "    0910_E1                           "
+  echo "    0910_E2                           "
+  echo "    1011_E1                           "
+  echo "    1011_E2                           "
+  exit
+fi
+
+exp=$1
+
+echo " Linking entry parameters and files for $exp"
+rm -f *.dat
+rm -f gotmrun.mal
+rm -f obs.nml
+rm -f airsea.mal
+ln -sf $exp/enp/*.mal .
+ln -sf $exp/*.nml .
+ln -sf $exp/*.dat .
+ln -sf enp_files/*.dat .
+
+if [ ! airsea.mal ] || [ ! gotmrun.mal ]
+then
+  echo " ERROR : .mal files missing"
+  exit
+else
+
+echo "-----------------------------------------------------"
+echo " run_enp.sh"
+echo "-----------------------------------------------------"
+
+
+# Parameters for the forcings
+
+
+count=1
+
+for file in ecmwf_meteo_*.dat
+do
+
+count=`expr 000$count | tail -4c`
+#count=`expr 000$count | tail -8c`
+echo $count
+ 
+if [ -f airsea.nml ]
+then
+  rm -f airsea.nml
+fi
+
+cat airsea.mal | sed \
+    -e "s/meteofile/ecmwf_meteo_$count.dat/g" \
+    > airsea.nml
+
+if [ -f airsea.nml ]
+then
+  echo "   airsea.nml created"
+fi
+
+if [ -f gotmrun.nml ]
+then
+  rm -f gotmrun.nml
+fi
+
+cat gotmrun.mal | sed \
+    -e "s/experiment_name/$count/g" \
+    > gotmrun.nml 
+
+if [ -f gotmrun.nml ]
+then
+  echo "   gotmrun.nml    created"
+fi
+
+strt=`date +%H:%M:%S`
+echo "   runnning with "$file
+echo "   started at "$strt
+
+./gotm_prod_IFORT >& $count.out
+
+end=`date +%H:%M:%S`
+echo "   ended   at "$end
+
+mv $count.nc $exp/enp
+mv $count.out $exp/enp/
+mv ${count}_fasham.par $exp/enp/
+
+count=`expr $count + 1`
+
+# Uncoment for debugging
+if [ `expr $count` == 4 ]
+then
+  echo "   Exiting"
+  exit
+fi
+
+done
+
+fi
+
+exit

src/extras/bio/Makefile

@@ -16,8 +16,9 @@ bio_sed.F90		\
 bio_fasham.F90		\
 bio_ismer.F90		\
 bio_gsj.F90             \
-bio_save.F90		\
 bio_nocera.F90		\
+bio_polynow.F90		\
+bio_save.F90		\
 bio_npzd4.F90		\
 nitrate.F90		\
 ammonium.F90            \
@@ -34,6 +35,7 @@ ${LIB}(bio_ismer.o)	\
 ${LIB}(bio_gsj.o)	\
 ${LIB}(bio_fasham.o)	\
 ${LIB}(bio_nocera.o)	\
+${LIB}(bio_polynow.o)	\
 ${LIB}(bio_npzd4.o)	\
 ${LIB}(bio_sed.o)	\
 ${LIB}(bio_save.o)	\

src/extras/bio/ammonium.F90

@@ -180,6 +180,11 @@
       do i=1,nlev
          amm(i) = cc(6,i)
       end do
+   else if (bio_model.eq.9) then  ! GG-A
+     do i=1,nlev
+        amm(i) = cc(6,i)
+     end do
+
    end if
 #endif
 

src/extras/bio/bio.F90

@@ -48,9 +48,13 @@
    use bio_mab,     only : init_bio_mab,init_var_mab,var_info_mab
    use bio_mab,     only : light_mab,surface_fluxes_mab,do_bio_mab
 
+   use bio_polynow,  only : init_bio_polynow,init_var_polynow,var_info_polynow  !GG-A
+   use bio_polynow,  only : light_polynow,do_bio_polynow                        !GG-A
+
    use output,      only : out_fmt,write_results,ts
 
    use util
+
 !
 !  default: all is private.
    private
@@ -58,6 +62,9 @@
 ! !PUBLIC MEMBER FUNCTIONS:
    public init_bio, set_env_bio, do_bio, get_bio_updates, clean_bio
    logical, public                     :: bio_calc=.false.
+
+
+!  REALTYPE,public, dimension(:), allocatable :: eps_ !GG-T
 !
 ! !REVISION HISTORY:!
 !  Original author(s): Hans Burchard & Karsten Bolding
@@ -353,6 +360,17 @@
 
          call var_info_npzd4()
 
+      case (10)  ! The model for the north Water polynya - NPZD- Sedimentation GG-A
+      
+         call init_bio_polynow(namlst,'bio_polynow.nml',unit)
+
+         call allocate_memory(nlev)
+      
+         call init_var_polynow(nlev)
+     
+         call var_info_polynow()
+
+
       case default
          stop "bio: no valid biomodel specified in bio.nml !"
       end select
@@ -448,9 +466,9 @@
 ! modules
 !
 ! !INTERFACE: 
-  subroutine set_env_bio(nlev,h_,t_,s_,nit_,amm_,hcb_,rho_,nuh_,rad_,wind_, &
+ subroutine set_env_bio(nlev,h_,t_,s_,nit_,amm_,hcb_,rho_,nuh_,rad_,wind_, &
                           I_0_,w_,w_adv_ctr_)
-!
+
 ! !DESCRIPTION:
 !
 ! !USES:
@@ -471,6 +489,8 @@
    REALTYPE, intent(in)                :: I_0_
    REALTYPE, optional, intent(in)      :: w_(0:nlev)
    integer, optional, intent(in)       :: w_adv_ctr_
+
+
 !
 ! !REVISION HISTORY:
 !  Original author(s): Hans Burchard & Karsten Bolding
@@ -494,6 +514,9 @@
    if (present(w_)) w = w_
    if (present(w_adv_ctr_)) w_adv_ctr = w_adv_ctr_
 
+
+  
+
    return
    end subroutine set_env_bio
 !EOC
@@ -742,6 +765,10 @@
             case (9)
                call light_npzd4(nlev,bioshade_feedback)
                call ode_solver(ode_method,numc,nlev,dt_eff,cc,do_bio_npzd4)
+            case (10)
+               call light_polynow(nlev,bioshade_feedback)                     
+               call ode_solver(ode_method,numc,nlev,dt_eff,cc,do_bio_polynow)
+
          end select
 
       end do
@@ -820,6 +847,17 @@
    if (allocated(var_units))      deallocate(var_units)
    if (allocated(var_long))       deallocate(var_long)
 
+!Stickiness
+   if (allocated(sti_2p))             deallocate(sti_2p)
+   if (allocated(sti_pdph))           deallocate(sti_pdph)
+   if (allocated(sti_pdzo))           deallocate(sti_pdzo)
+   if (allocated(sti_pfp))            deallocate(sti_pfp)
+   if (allocated(sti_2dph))           deallocate(sti_2dph)
+   if (allocated(sti_dphdzo))         deallocate(sti_dphdzo)
+   if (allocated(sti_dphfp))          deallocate(sti_dphfp)
+   if (allocated(sti_2dzo))           deallocate(sti_2dzo)
+   if (allocated(sti_dzofp))          deallocate(sti_dzofp)
+   if (allocated(sti_2fp))            deallocate(sti_2fp)
 !DD
    if (allocated(lumlim1))        deallocate(lumlim1)
    if (allocated(nitlim1))        deallocate(nitlim1)
@@ -829,6 +867,15 @@
    if (allocated(ammlim2))        deallocate(ammlim2)
    if (allocated(ppnet))          deallocate(ppnet)
 
+   if (allocated(flux_msn))        deallocate(flux_msn)
+   if (allocated(Flux_P))          deallocate(Flux_P)
+   if (allocated(Flux_D1))         deallocate(Flux_D1)
+   if (allocated(Flux_D2))         deallocate(Flux_D2)
+   if (allocated(Flux_D3))         deallocate(Flux_D3)
+
+   if (allocated(size_msnow))         deallocate(size_msnow)
+   if (allocated(w_msn_lev))         deallocate(w_msn_lev)
+
 !  The external provide arrays
    if (allocated(h))              deallocate(h)
    if (allocated(nuh))            deallocate(nuh)
@@ -912,6 +959,29 @@
    allocate(var_long(numc),stat=rc)
    if (rc /= 0) stop 'init_bio(): Error allocating var_long)'
 
+!Stickiness
+   allocate(sti_2p(0:nlev),stat=rc)
+   if (rc /= 0) STOP 'init_bio: Error allocating (sti_2p)'
+   allocate(sti_pdph(0:nlev),stat=rc)
+   if (rc /= 0) STOP 'init_bio: Error allocating (sti_pdph)'
+   allocate(sti_pdzo(0:nlev),stat=rc)
+   if (rc /= 0) STOP 'init_bio: Error allocating (sti_pdzo)'
+   allocate(sti_pfp(0:nlev),stat=rc)
+   if (rc /= 0) STOP 'init_bio: Error allocating (sti_pfp)'
+   allocate(sti_2dph(0:nlev),stat=rc)
+   if (rc /= 0) STOP 'init_bio: Error allocating (sti_2dph)'
+   allocate(sti_dphdzo(0:nlev),stat=rc)
+   if (rc /= 0) STOP 'init_bio: Error allocating (sti_dphdzo)'
+   allocate(sti_dphfp(0:nlev),stat=rc)
+   if (rc /= 0) STOP 'init_bio: Error allocating (sti_dphfp)'
+   allocate(sti_2dzo(0:nlev),stat=rc)
+   if (rc /= 0) STOP 'init_bio: Error allocating (sti_2dzo)'
+   allocate(sti_dzofp(0:nlev),stat=rc)
+   if (rc /= 0) STOP 'init_bio: Error allocating (sti_dzofp)'
+   allocate(sti_2fp(0:nlev),stat=rc)
+   if (rc /= 0) STOP 'init_bio: Error allocating (sti_2fp)'
+
+
 !DD
 
    allocate(lumlim1(0:nlev),stat=rc)
@@ -935,6 +1005,28 @@
    allocate(ppnet(0:nlev),stat=rc)
    if (rc /= 0) stop 'init_bio(): Error allocating (ppnet)'
 
+
+
+   allocate(flux_msn(0:nlev),stat=rc)
+   if (rc /= 0) stop 'init_bio(): Error allocating (flux_msn)'
+
+   allocate(Flux_P(0:nlev),stat=rc)
+   if (rc /= 0) stop 'init_bio(): Error allocating (Flux_P)'
+   allocate(Flux_D1(0:nlev),stat=rc)
+   if (rc /= 0) stop 'init_bio(): Error allocating (Flux_D1)'
+ allocate(Flux_D2(0:nlev),stat=rc)
+   if (rc /= 0) stop 'init_bio(): Error allocating (Flux_D2)'
+ allocate(Flux_D3(0:nlev),stat=rc)
+   if (rc /= 0) stop 'init_bio(): Error allocating (Flux_D3)'
+
+ allocate(size_msnow(0:nlev),stat=rc)
+   if (rc /= 0) stop 'init_bio(): Error allocating (size_msnow)'
+
+
+ allocate(w_msn_lev(0:nlev),stat=rc)
+   if (rc /= 0) stop 'init_bio(): Error allocating (w_msn_lev)'
+
+
 !  The external provide arrays
    allocate(h(0:nlev),stat=rc)
    if (rc /= 0) stop 'init_bio(): Error allocating (h)'

src/extras/bio/bio_gsj.F90

@@ -64,7 +64,8 @@
    REALTYPE                  ::  p0      = 0.0
    REALTYPE                  ::  z0      = 0.0
    REALTYPE                  ::  b0      = 0.0
-   LOGICAL                   ::  mte     = .true.
+   LOGICAL                   ::  mte     = .false.
+   LOGICAL                   ::  txc4zoo = .false.
    REALTYPE                  ::  ca1     = 3.61
    REALTYPE                  ::  ca2     = 14.58    
    REALTYPE                  ::  ch1     = 3.265
@@ -177,6 +178,7 @@
 ! mu22    : Maximum large zooplankton loss rate
 ! k6      : Half saturation constant of zooplankton loss (small and large)
 !! k7	  : Half saturation constant of zooplankton loss by hydrocarbon toxicity (small and large)
+! txc4zoo : Logical to make hydrocarbon toxic for zooplankton
 ! bhc     : Slope factor of zooplankton mortality for sigmoïd exponential type curve
 ! mhc     : Factor of zooplankton mortality for Michaelis-Menten exponential type curve
 ! hc0     : Hydrocarbon sill concentration for zooplankton mortality by hydrocarbons
@@ -241,7 +243,7 @@
                         r11,r12,r13,r14,r21,r22,r23,r24,                 &
                         vb1,vb2,k4,k10,w_h,mu3,etaa,etah,mu4,w_d,kc,mu5, &
                         theta,w_p1max,w_p1min,w_p2min,w_p2max,           &
-                        mte,ca1,ca2,ch1,ch2,amratio,hmratio
+                        txc4zoo,mte,ca1,ca2,ch1,ch2,amratio,hmratio
 
 !EOP
 !-----------------------------------------------------------------------
@@ -766,10 +768,15 @@
       nitlim2(ci) =qn2
       ammlim2(ci) =qa2
    
-     ! Hydrocarbon inhibition factors
-     !mu21hc = mu21+(mu21max-mu21)/(1d0 + exp(-1d0*(cc(hc,ci)-k6hc)/bhc))
-     mu21hc = mu21+((mu21max-mu21)*(cc(hc,ci))**mhc)/((cc(hc,ci))**mhc + k6hc**mhc)
-     mu22hc = mu21hc
+      ! Hydrocarbon inhibition factors
+      !mu21hc = mu21+(mu21max-mu21)/(1d0 + exp(-1d0*(cc(hc,ci)-k6hc)/bhc))
+      if (txc4zoo) then
+         mu21hc = mu21+((mu21max-mu21)*(cc(hc,ci))**mhc)/((cc(hc,ci))**mhc + k6hc**mhc)
+         mu22hc = mu21hc
+      else
+         mu21hc = mu21
+         mu22hc = mu22
+      endif
       
       ! Nutrient uptake by pico- and nano-phytoplankton 
       dd(n,p1,ci) =amr1*ff1*qn1*(cc(p1,ci)+p0)

src/extras/bio/bio_save.F90

@@ -14,6 +14,7 @@
 ! !USES:
    use bio_var
    use output,  only: out_fmt,ts
+
 #ifdef NETCDF_FMT
    use ncdfout, only: ncid
    use ncdfout, only: lon_dim,lat_dim,z_dim,time_dim,dims
@@ -28,6 +29,7 @@
 ! !INPUT PARAMETERS:
    integer, intent(in)                 :: nlev
    REALTYPE, intent(in)                :: totn
+
 !
 ! !REVISION HISTORY:
 !  Original author(s): Hans Burchard & Karsten Bolding
@@ -48,6 +50,9 @@
       first=.true.
    end if
    select case (out_fmt)
+!-----------------------------------------------------------------------------------
+! Ascii format for the output
+!-----------------------------------------------------------------------------------
       case (ASCII)
          if(first) then
             open(out_unit,file='bio.out',status='unknown')
@@ -64,6 +69,9 @@
          end do
 115 format(F10.4,100(1x,E10.4E2))
 
+!-----------------------------------------------------------------------------------
+! NETCDF format for the output
+!-----------------------------------------------------------------------------------
       case (NETCDF)
 #ifdef NETCDF_FMT
          if(first) then
@@ -82,24 +90,100 @@
                                      units=var_units(n),    &
                                      long_name=var_long(n))
             end do
+!-----------------------------------------------------------------------------------
+!Prepare the new variable for the NETCDF output file as well as informations
+!-----------------------------------------------------------------------------------
+    
+!Density of the fluid at each level
+            iret = new_nc_variable(ncid,'rho',NF_REAL,4,dims,rho_id)
+            iret = set_attributes(ncid,rho_id,units='kg.m-3',long_name='Density')
 
-            ! Sedimentation or swimming rate of particulate matter
+! Sedimentation or swimming rate of particulate matter
             ! Living phytoplankton
             iret = new_nc_variable(ncid,'wp',NF_REAL,4,dims,wp_id)
             iret = set_attributes(ncid,wp_id,units='m/day', &
-                                  long_name='phytoplancton settling velocity')
+                                  long_name='Phytoplancton settling velocity')
 
-            ! Living zooplankton
+           ! Living zooplankton
+        !Nocera Model
             if (bio_model .eq. 8) then
                iret = new_nc_variable(ncid,'wz',NF_REAL,4,dims,wz_id)
                iret = set_attributes(ncid,wz_id,units='m/day', &
                                   long_name='zooplancton swimming velocity')
             end if
-            !CHG4 Diagnostic du PAR
+         
+         !Polynow model
+            if (bio_model .eq. 10) then
+            ! Living zooplankton
+               iret = new_nc_variable(ncid,'wz',NF_REAL,4,dims,wz_id)
+               iret = set_attributes(ncid,wz_id,units='m/day', &
+                                  long_name='Zooplancton swimming velocity')
+            ! Dead Phytoplankton
+               iret = new_nc_variable(ncid,'wd1',NF_REAL,4,dims,wd1_id)
+               iret = set_attributes(ncid,wd1_id,units='m/day', &
+                                  long_name='Dead Phytoplankton settling velocity')
+            ! Dead zooplankton
+               iret = new_nc_variable(ncid,'wd2',NF_REAL,4,dims,wd2_id)
+               iret = set_attributes(ncid,wd2_id,units='m/day', &
+                                  long_name='Dead Zooplankton settling velocity')
+            ! Fecals pellets
+               iret = new_nc_variable(ncid,'wd3',NF_REAL,4,dims,wd3_id)
+               iret = set_attributes(ncid,wd3_id,units='m/day', &
+                                  long_name='Fecal pellets settling velocity')
+            ! Marine Snow
+               iret = new_nc_variable(ncid,'wd4',NF_REAL,4,dims,wd4_id)
+               iret = set_attributes(ncid,wd4_id,units='m/day', &
+                                  long_name='Marine snow settling velocity (Calc)')
+
+       ! Stickiness
+            ! 2 Living phytoplankton
+               iret = new_nc_variable(ncid,'sti_2p',NF_REAL,4,dims,sti_2p_id)
+               iret = set_attributes(ncid,sti_2p_id, &
+                                 long_name='Stickiness Phy & Phy')
+            ! Living Phytoplankton with dead phytoplankton
+               iret = new_nc_variable(ncid,'sti_pdph',NF_REAL,4,dims,sti_pdph_id)
+               iret = set_attributes(ncid,sti_pdph_id, &
+                                 long_name='Stickiness Phy & Dph')
+            ! Living Phytoplankton with dead zooplankton
+               iret = new_nc_variable(ncid,'sti_pdzo',NF_REAL,4,dims,sti_pdzo_id)