bio_polynow.F90 128 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
!$id: bio_polynow.F90,v 1.11 2016-05-05 11:49:15 dd Exp $
#include"cppdefs.h"
!-----------------------------------------------------------------------
!BOP
!
! !MODULE: bio_fasham --- Fasham et al. biological model \label{sec:bio-fasham}
!
! !INTERFACE:
   module bio_polynow
!
!----DESCRIPTION:--------------------------------------------------------------
!  The model developed by Fashametal1990 
!  uses nitrogen as 'currency' according to the evidence that in
!  most cases nitrogen is the limiting macronutrient. It consists of
!  seven state variables: phytoplankton, zooplankton, bacteria,
!  particulate organic matter (detritus), dissolved organic matter
!  and the nutrients nitrate and ammonium.
Gwenaelle Gremion's avatar
Gwenaelle Gremion committed
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
!  The structure of the biogeochemical model in 14 states variables
!1 - Phytoplankton (phy) 
!2 - Zooplankton (zoo)
!3 - Bacteria (bac),
! NUTRIENTS ---
!5- Nitrate (nit)
!6- Ammonium (amm)
!7- Labile dissolved organic nitrogen (ldn)
! DETRITUS ----
!4 -Dead phytoplankton(Dph)
!8 -Dead zooplankton (Dzo)
!9 -Fecal pellets (Fp)
!10 - Marine snow (Msn)
! TRAITS of Marine Snow ---
!11 - Size of marine snow after coagulation
!12 - Size of marine snow when coagulation happens
!13 - Size of marine snow when fragmentation occurs
!14 - Final size of marine snow

37
38
39
40
41
42
43
44
45
46

! The concentrations are in mmol N\m^-3,and all fluxes are conservative.
!---------------------------------------------------------------------------------

! !USES:
!  default: all is private.
   use bio_var       ! S,T,zlev,rho
   use output
   use observations, only : aa,g2
   use turbulence,   only : eps  
Gwenaelle Gremion's avatar
Gwenaelle Gremion committed
47
   use meanflow,     only : Rho_0
48
49
50
51
52
53
54
 
   private
!
! !PUBLIC MEMBER FUNCTIONS:
   public init_bio_polynow, init_var_polynow, var_info_polynow, &
          light_polynow, do_bio_polynow, end_bio_polynow
   REALTYPE, public,parameter   :: pi= 3.141592654
Gwenaelle Gremion's avatar
Gwenaelle Gremion committed
55
56
   REALTYPE, public,parameter   :: g= 9.80665     !m/s2 
   REALTYPE                     :: pres = 10.1325 !  gauge pressure (absolute pressure - 10.1325 bar)
57
58

!-----LOCAL VARIABLES:---------- from a namelist : bio_polynow.nml----------------------
59

Gwenaelle Gremion's avatar
Gwenaelle Gremion committed
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
!-----These parameters need to be similar to those in gotmrun.nml (dt, depth) or those in bio.nml (split_factor)----------------------
  REALTYPE                  :: dt_bio 
  REALTYPE                  :: splitfac_bio
  REALTYPE                  :: depth_bio

!-----INITIAL and Minimum concentration for the variables 

  REALTYPE                  ::  p_init_value  = 1.0
  REALTYPE                  ::  p_initial     = 0.01
  REALTYPE                  ::  z_p_gauss_init= 2.0
  REALTYPE                  ::  sigma_p       = 2.0

  REALTYPE                  ::  zoo_init_value = 1.0
  REALTYPE                  ::  z_initial      = 0.01
  REALTYPE                  ::  z_zoo_gauss_init= 2.0
  REALTYPE                  ::  sigma_zoo      = 2.0

  REALTYPE                  ::  b_initial      = 0.001

  REALTYPE                  ::  dph_init_value = 1.0
  REALTYPE                  ::  dph_initial    = 0.01
  REALTYPE                  ::  z_dph_gauss_init= 2.0
  REALTYPE                  ::  sigma_dph      = 2.0

  REALTYPE                  ::  dzo_init_value = 1.0
  REALTYPE                  ::  dzo_initial    = 0.01
  REALTYPE                  ::  z_dzo_gauss_init= 2.0
  REALTYPE                  ::  sigma_dzo      = 2.0

  REALTYPE                  ::  fp_init_value  = 1.0
  REALTYPE                  ::  fp_initial     = 0.01
  REALTYPE                  ::  z_fp_gauss_init= 2.0
  REALTYPE                  ::  sigma_fp       = 2.0

  REALTYPE                  ::  msn_init_value = 1.0
  REALTYPE                  ::  msn_initial    = 0.01
  REALTYPE                  ::  z_msn_gauss_init= 2.0
  REALTYPE                  ::  sigma_msn      = 2.0

  REALTYPE                  ::  l_initial      = 0.14
  REALTYPE                  ::  p0             = 0.0
  REALTYPE                  ::  z0             = 0.0
  REALTYPE                  ::  b0             = 0.0
  REALTYPE                  ::  mu5            = 0.02  
!-----Phytoplankton and dead-phytoplankton
  REALTYPE                  ::  vp             = 1.5
  REALTYPE                  ::  alpha          = 0.065
  REALTYPE                  ::  inib           = 0.05
  REALTYPE, public          ::  kc             = 0.03
  REALTYPE                  ::  k1             = 0.2
  REALTYPE                  ::  k2             = 0.8
  REALTYPE                  ::  mu1            = 0.05
  REALTYPE                  ::  k5             = 0.2
  REALTYPE                  ::  gamma          = 0.05
  REALTYPE                  ::  txloss_p       = 0.7  
  REALTYPE                  ::  txloss_dph     = 0.05  
!-----Zooplankton
  REALTYPE                  ::  gmax           = 1.0
  REALTYPE                  ::  k3             = 1.0
  REALTYPE                  ::  beta           = 0.625
  REALTYPE                  ::  mu2            = 0.3
  REALTYPE                  ::  k6             = 0.2
  REALTYPE                  ::  delta          = 0.1
  REALTYPE                  ::  epsi           = 0.70
  REALTYPE                  ::  eg             = 0.05   
  REALTYPE                  ::  r1             = 0.55   ! phy
  REALTYPE                  ::  r2             = 0.4    ! zoo
  REALTYPE                  ::  r3             = 0.05   ! dph  
  REALTYPE                  ::  r4             = 0.05   ! dzo 
  REALTYPE                  ::  r5             = 0.05   ! fp   
  REALTYPE                  ::  r6             = 0.05   ! msn 
  REALTYPE                  ::  txloss_dzo     = 0.05  
  REALTYPE                  ::  txloss_fp      = 0.05 
!-----Zooplankton Vertical migration(From Ariadna Nocera)
  REALTYPE                  ::  Migra_zoo      = 1.0   
  REALTYPE                  ::  pmin           = 0.05
  REALTYPE                  ::  w_zmax         = 100.0  
  REALTYPE                  ::  bertha         = 0.05   
  REALTYPE                  ::  parcrit        = 0.02   
!-----Bacterias and LDON
  REALTYPE                  ::  vb             = 1.2
  REALTYPE                  ::  remi           = 0.1
  REALTYPE                  ::  k4             = 0.5
  REALTYPE                  ::  mu3            = 0.15
  REALTYPE                  ::  eta            = 0.0
  REALTYPE                  ::  mbac           = 0.0    
  REALTYPE                  ::  dphlossb       = 0.4
  REALTYPE                  ::  dphlossl       = 10.0  
  REALTYPE                  ::  dzolossl       = 10.0  
  REALTYPE                  ::  dzolossb       = 10.0  
  REALTYPE                  ::  fplossl        = 10.0 
  REALTYPE                  ::  fplossb        = 10.0 
  REALTYPE                  ::  leak           = 0.1
  REALTYPE                  ::  mldon          = 0.02    
  REALTYPE                  ::  lmin           = 0.02    
!----- Detritus/Settling
  REALTYPE                  ::  Phys_w         = 0.0 
  REALTYPE                  ::  w_msnow        = 0.0 
  REALTYPE                  ::  w_p            = -0.5
  REALTYPE                  ::  w_dph          = -1.0   
  REALTYPE                  ::  w_dzo          = -10.0   
  REALTYPE                  ::  w_fp           = -100.0  
  REALTYPE                  ::  w_msn          = -2.0    
  REALTYPE                  ::  rho_p          = 0.02 
  REALTYPE                  ::  rho_dph        = 0.02 
  REALTYPE                  ::  rho_dzo        = 0.02 
  REALTYPE                  ::  rho_fp         = 0.02 
  REALTYPE                  ::  rho_msn        = 0.02 
!----- Coagulation  
  REALTYPE                  ::  Coag_coef      = 0.0     
  REALTYPE                  ::  betap_p        = 0.2
  REALTYPE                  ::  betap_dph      =  1
  REALTYPE                  ::  betap_dzo      = 1
  REALTYPE                  ::  betap_fp       = 1 
  REALTYPE                  ::  betadph_dph    = 0.2
  REALTYPE                  ::  betadph_dzo    = 1
  REALTYPE                  ::  betadph_fp     = 1
  REALTYPE                  ::  betadzo_dzo    = 0.02    
  REALTYPE                  ::  betadzo_fp     = 0.02  
  REALTYPE                  ::  betafp_fp      = 0.02   
  REALTYPE                  ::  betaBr         = 1.0
  REALTYPE                  ::  betaSh         = 1.0 
  REALTYPE                  ::  betaDs         = 1.0 
  REALTYPE                  ::  eps_const      = 1.0 
  REALTYPE                  ::  eps_n          = 1.0 
  REALTYPE                  ::  cons_min       = 0.00001 
  REALTYPE                  ::  sti_cst        = 0.0
  REALTYPE                  ::  stip_p         = -0.5  
  REALTYPE                  ::  stip_dph       = -0.5  
  REALTYPE                  ::  stip_dzo       = -0.5   
  REALTYPE                  ::  stip_fp        = -0.5   
  REALTYPE                  ::  stidph_dph     = -0.5  
  REALTYPE                  ::  stidph_dzo     = -0.5  
  REALTYPE                  ::  stidph_fp      = -0.5   
  REALTYPE                  ::  stidzo_dzo     = -0.5   
  REALTYPE                  ::  stidzo_fp      = -0.5   
  REALTYPE                  ::  stifp_fp       = -0.5   
!-----  Settling
  REALTYPE                  ::   CSF           =1.0
  REALTYPE                  ::   size_rand     = 1.0
  REALTYPE                  ::   size_phy_us   = -0.5  
  REALTYPE                  ::   size_phy_up   = -0.5 
  REALTYPE                  ::   size_dph_us   = 0.00001
  REALTYPE                  ::   size_dph_up   = 0.00001
  REALTYPE                  ::   size_dzo_us   = 0.000100
  REALTYPE                  ::   size_dzo_up   = 0.000100
  REALTYPE                  ::   size_fp_us    = 0.000050
  REALTYPE                  ::   size_fp_up    = 0.000050 
  REALTYPE                  ::   dm_msn        = 0.0
  REALTYPE                  ::   coef1         = 0.0
  REALTYPE                  ::   coef2         = 0.0
  REALTYPE                  ::   diam_msn_us   = 0.01  
  REALTYPE                  ::   coef3         = -0.5  
  REALTYPE                  ::   cons_max      = -0.5  
  REALTYPE                  ::   dynvis        = -0.5  
  REALTYPE                  ::   kinvis        = -0.5  
  REALTYPE                  ::   kB            = -0.5  
!-----Fragmentation
  REALTYPE                   :: Frag_meth      = 0.0 
  REALTYPE                   :: swim_brk       = 0.8
  REALTYPE                   :: Floc_coef      = 0.1
  REALTYPE                   :: Daugther_part  = 1.0

   integer                   ::  out_unit
!Variables declarations
225
 integer, parameter        ::  p=1,z=2,b=3,d1=4,n=5,a=6,l=7,d2=8,d3=9,d4=10, &
Gwenaelle Gremion's avatar
Gwenaelle Gremion committed
226
                               size_intrm=11,size_coag=12,size_frag=13,size_msn=14
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260

!EOP
!-----------------------------------------------------------------------

   contains

!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: Initialise the bio module
!
! !INTERFACE:
   subroutine init_bio_polynow(namlst,fname,unit)
!
! !DESCRIPTION:
!  Here, the bio namelist bio_polynow.nml is read and
!  various variables (rates and settling velocities)
!  are transformed into SI units.
!
! !USES:
   IMPLICIT NONE
!
! !INPUT PARAMETERS:
   integer,          intent(in)   :: namlst
   character(len=*), intent(in)   :: fname
   character(len=20)              :: pfile
   integer,          intent(in)   :: unit
!
! !REVISION HISTORY:
!  Original author(s): Hans Burchard & Karsten Bolding
!
! !LOCAL VARIABLES:
   namelist /bio_polynow_nml/ numc, &
                                                        
Gwenaelle Gremion's avatar
Gwenaelle Gremion committed
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
   dt_bio,splitfac_bio,depth_bio,&
!-----INITIAL and Minimum concentration for the variables 
   p_init_value,p_initial,z_p_gauss_init,sigma_p,          &
   zoo_init_value,z_initial,z_zoo_gauss_init,sigma_zoo,    &
   b_initial,                                              &
   dph_init_value,dph_initial,z_dph_gauss_init,sigma_dph,  &
   dzo_init_value,dzo_initial,z_dzo_gauss_init,sigma_dzo,  &
   fp_init_value,fp_initial,z_fp_gauss_init,sigma_fp,      &
   msn_init_value,msn_initial,z_msn_gauss_init,sigma_msn,  & 
   l_initial,p0,z0,b0,mu5,                                 &
!-----Phytoplankton and dead-phytoplankton
   vp,alpha,inib,kc,k1,k2,mu1,k5,gamma,txloss_p,txloss_dph,& 
!-----Zooplankton
   gmax,k3,beta,mu2,k6,delta,epsi,eg,r1,r2,r3,r4,r5,r6,    & 
   txloss_dzo,txloss_fp, &
   Migra_zoo,pmin,w_zmax,bertha,parcrit,                   & 
!-----Bacterias and LDON
   vb,remi,k4,mu3,eta,mbac,dphlossb,dphlossl,dzolossl,dzolossb,fplossl,fplossb, & 
   leak,mldon,lmin,                                        & 
!----- Detritus/Settling
   Phys_w,w_msnow,w_p,w_dph,w_dzo,w_fp,w_msn,              & 
   rho_p,rho_dph,rho_dzo,rho_fp,rho_msn,                   &  
!----- Coagulation  
   Coag_coef,                                              &
   betap_p,betap_dph,betap_dzo,betap_fp,                   &
   betadph_dph,betadph_dzo,betadph_fp,                     &
   betadzo_dzo,betadzo_fp,betafp_fp,                       & 
   betaBr,betaSh,betaDs,                                   &
   eps_const,eps_n,cons_min,                               &
   sti_cst,                                                &
   stip_p,stip_dph,stip_dzo,stip_fp,stidph_dph,stidph_dzo,stidph_fp,stidzo_dzo,stidzo_fp,stifp_fp, &
!-----  Settling
   CSF,size_rand,size_phy_us,size_phy_up,size_dph_us,size_dph_up,size_dzo_us,size_dzo_up,&
   size_fp_us,size_fp_up,dm_msn,coef1,coef2,diam_msn_us,   &
   coef3,cons_max,                                         &
   dynvis,kinvis,kB,                                       &
!-----Fragmentation
   Frag_meth,swim_brk,Floc_coef, Daugther_part                
             
300
301
302
303
304
305
306
307
308
309
310
311
312
313
!EOP
!-----------------------------------------------------------------------
!BOC
   LEVEL2 'init_bio_polynow'

   open(namlst,file=fname,action='read',status='old',err=98)
   read(namlst,nml=bio_polynow_nml,err=99)
   close(namlst)

   numcc=numc
!----Print some parameter values in standard output and save them in a separate file [out_fn]_polynow.par
   pfile = trim(out_fn) // '_polynow.par'
   open(10,status='unknown',action='write',file=pfile)
   LEVEL3 'Biogeochemical parameters saved in ', pfile
Gwenaelle Gremion's avatar
Gwenaelle Gremion committed
314
315
316
317
318
319
320
321
322
323
  
   write(*,900) '                dt_bio     = ',dt_bio
   write(10,901) dt_bio
   write(*,900) '                splitfac_bio  = ',splitfac_bio
   write(10,901) splitfac_bio
   write(*,900) '                depth_bio     = ',depth_bio
   write(10,901) depth_bio
!-----INITIAL and Minimum concentration for the variables 
   write(*,900) '                p_init_value      = ',p_init_value 
   write(10,901) p_init_value 
324
   write(*,900) '                p_initial     = ',p_initial
Gwenaelle Gremion's avatar
Gwenaelle Gremion committed
325
326
327
328
329
330
331
   write(10,901) p_initial
   write(*,900) '                z_p_gauss_init     = ',z_p_gauss_init
   write(10,901) z_p_gauss_init
   write(*,900) '                sigma_p     = ',sigma_p
   write(10,901) sigma_p
   write(*,900) '                zoo_init_value      = ',zoo_init_value 
   write(10,901) zoo_init_value 
332
   write(*,900) '                z_initial     = ',z_initial
Gwenaelle Gremion's avatar
Gwenaelle Gremion committed
333
334
335
336
337
   write(10,901) z_initial
   write(*,900) '                z_zoo_gauss_init     = ',z_zoo_gauss_init
   write(10,901) z_zoo_gauss_init
   write(*,900) '                sigma_zoo     = ',sigma_zoo
   write(10,901) sigma_zoo
338
   write(*,900) '                b_initial     = ',b_initial
Gwenaelle Gremion's avatar
Gwenaelle Gremion committed
339
340
341
   write(10,901) b_initial
   write(*,900) '                dph_init_value      = ',dph_init_value 
   write(10,901) dph_init_value 
342
   write(*,900) '                dph_initial     = ',dph_initial
Gwenaelle Gremion's avatar
Gwenaelle Gremion committed
343
344
345
346
347
348
349
   write(10,901) dph_initial
   write(*,900) '                z_dph_gauss_init     = ',z_dph_gauss_init
   write(10,901) z_dph_gauss_init
   write(*,900) '                sigma_dph     = ',sigma_dph
   write(10,901) sigma_dph
   write(*,900) '                dzo_init_value      = ',dzo_init_value 
   write(10,901) dzo_init_value 
350
   write(*,900) '                dzo_initial     = ',dzo_initial
Gwenaelle Gremion's avatar
Gwenaelle Gremion committed
351
352
353
354
355
356
357
   write(10,901) dzo_initial
   write(*,900) '                z_dzo_gauss_init     = ',z_dzo_gauss_init
   write(10,901) z_dzo_gauss_init
   write(*,900) '                sigma_dzo     = ',sigma_dzo
   write(10,901) sigma_dzo
   write(*,900) '                fp_init_value      = ',fp_init_value 
   write(10,901) fp_init_value 
358
   write(*,900) '                fp_initial     = ',fp_initial
Gwenaelle Gremion's avatar
Gwenaelle Gremion committed
359
360
361
362
363
364
365
   write(10,901) fp_initial
   write(*,900) '                z_fp_gauss_init     = ',z_fp_gauss_init
   write(10,901) z_fp_gauss_init
   write(*,900) '                sigma_fp     = ',sigma_fp
   write(10,901) sigma_fp
   write(*,900) '                msn_init_value      = ',msn_init_value 
   write(10,901) msn_init_value 
366
   write(*,900) '                msn_initial     = ',msn_initial
Gwenaelle Gremion's avatar
Gwenaelle Gremion committed
367
368
369
370
371
   write(10,901) msn_initial
   write(*,900) '                z_msn_gauss_init     = ',z_msn_gauss_init
   write(10,901) z_msn_gauss_init
   write(*,900) '                sigma_msn     = ',sigma_msn
   write(10,901) sigma_msn
372
373
374
375
376
377
378
379
380
   write(*,900) '                l_initial     = ',l_initial
   write(10,901) l_initial
   write(*,900) '                p0     = ',p0
   write(10,901) p0
   write(*,900) '                z0     = ',z0
   write(10,901) z0
   write(*,900) '                b0     = ',b0
   write(10,901) b0
   write(*,900) '                mu5     = ',mu5
Gwenaelle Gremion's avatar
Gwenaelle Gremion committed
381
382
   write(10,901) mu5
!-----Phytoplankton and dead-phytoplankton 
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
   write(*,900) '                vp     = ',vp
   write(10,901) vp
   write(*,900) '                alpha  = ',alpha
   write(10,901) alpha
   write(*,900) '                inib   = ',inib
   write(10,901) inib
   write(*,900) '                kc     = ',kc
   write(10,901) kc
   write(*,900) '                k1     = ',k1
   write(10,901) k1
   write(*,900) '                k2     = ',k2
   write(10,901) k2
   write(*,900) '                mu1     = ',mu1
   write(10,901) mu1
   write(*,900) '                k5     = ',k5
   write(10,901) k5
   write(*,900) '               gamma     = ',gamma
   write(10,901) gamma
   write(*,900) '               txloss_p = ',txloss_p
   write(10,901) txloss_p
   write(*,900) '               txloss_dph = ',txloss_dph
Gwenaelle Gremion's avatar
Gwenaelle Gremion committed
404
405
   write(10,901) txloss_dph
!-----Zooplankton
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
   write(*,900) '                gmax  = ',gmax
   write(10,901) gmax
   write(*,900) '                k3     = ',k3
   write(10,901) k3
   write(*,900) '                beta     = ',beta
   write(10,901) beta
   write(*,900) '                mu2      = ',mu2
   write(10,901) mu2
   write(*,900) '                k6     = ',k6
   write(10,901) k6
   write(*,900) '                delta    = ',delta
   write(10,901) delta
   write(*,900) '                epsi    = ',epsi
   write(10,901) epsi
   write(*,900) '                eg     = ',eg
   write(10,901) eg
   write(*,900) '                r1     = ',r1
   write(10,901) r1
   write(*,900) '                r2   = ',r2
   write(10,901) r2
   write(*,900) '                r3    = ',r3
   write(10,901) r3
   write(*,900) '                r4    = ',r4
   write(10,901) r4
   write(*,900) '                r5     = ',r5
   write(10,901) r5
   write(*,900) '                r6     = ',r6
   write(10,901) r6
   write(*,900) '               txloss_dzo = ',txloss_dzo
   write(10,901) txloss_dzo
   write(*,900) '               txloss_fp = ',txloss_fp
   write(10,901) txloss_fp
   write(*,900) '               Migra_zoo = ',Migra_zoo
   write(10,901) Migra_zoo
   write(*,900) '                pmin     = ',pmin
   write(10,901) pmin
   write(*,900) '                w_zmax    = ',w_zmax
   write(10,901) w_zmax
   write(*,900) '                bertha    = ',bertha
   write(10,901) bertha
   write(*,900) '                parcrit     = ',parcrit
Gwenaelle Gremion's avatar
Gwenaelle Gremion committed
447
448
   write(10,901) parcrit
!-----Bacterias and LDON
449
450
   write(*,900) '                vb    = ',vb
   write(10,901) vb
Gwenaelle Gremion's avatar
Gwenaelle Gremion committed
451
452
   write(*,900) '                remi    = ',remi
   write(10,901) remi
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
   write(*,900) '                k4    = ',k4
   write(10,901) k4
   write(*,900) '                mu3    = ',mu3
   write(10,901) mu3
   write(*,900) '                eta    = ',eta
   write(10,901) eta
   write(*,900) '                mbac    = ',mbac
   write(10,901) mbac
   write(*,900) '                dphlossb    = ',dphlossb
   write(10,901) dphlossb
   write(*,900) '               dphlossl    = ',dphlossl
   write(10,901) dphlossl
   write(*,900) '               dzolossl    = ',dzolossl
   write(10,901) dzolossl
   write(*,900) '               dzolossb   = ',dzolossb
   write(10,901) dzolossb
   write(*,900) '               fplossl    = ',fplossl
   write(10,901) fplossl
   write(*,900) '               fplossb   = ',fplossb
   write(10,901) fplossb
Gwenaelle Gremion's avatar
Gwenaelle Gremion committed
473
474
   write(*,900) '                leak    = ',leak
   write(10,901) leak
475
476
477
   write(*,900) '               mldon   = ',mldon
   write(10,901) mldon
   write(*,900) '               lmin   = ',lmin
Gwenaelle Gremion's avatar
Gwenaelle Gremion committed
478
479
   write(10,901) lmin
!----- Detritus/Settling
480
   write(*,900) '               Phys_w   = ',Phys_w
Gwenaelle Gremion's avatar
Gwenaelle Gremion committed
481
482
483
   write(10,901) Phys_w
   write(*,900) '               w_msnow   = ',w_msnow
   write(10,901) w_msnow
484
485
486
487
488
489
490
491
492
493
   write(*,900) '               w_p   = ',w_p
   write(10,901) w_p
   write(*,900) '               w_dph    = ',w_dph
   write(10,901) w_dph
   write(*,900) '               w_dzo    = ',w_dzo
   write(10,901) w_dzo
   write(*,900) '               w_fp  = ',w_fp
   write(10,901) w_fp
   write(*,900) '               w_msn    = ',w_msn
   write(10,901) w_msn
Gwenaelle Gremion's avatar
Gwenaelle Gremion committed
494
495
496
497
   write(*,900) '               rho_p    = ',rho_p
   write(10,901) rho_p
   write(*,900) '               rho_dph    = ',rho_dph
   write(10,901) rho_dph
498
499
500
501
502
   write(*,900) '               rho_dzo    = ',rho_dzo
   write(10,901) rho_dzo
   write(*,900) '               rho_fp    = ',rho_fp
   write(10,901) rho_fp
   write(*,900) '               rho_msn    = ',rho_msn
Gwenaelle Gremion's avatar
Gwenaelle Gremion committed
503
504
505
506
   write(10,901) rho_msn
!----- Coagulation 
   write(*,900) '               coag_coef    = ',Coag_coef
   write(10,901) Coag_coef
507
508
   write(*,900) '               betap_p   = ',betap_p
   write(10,901) betap_p
Gwenaelle Gremion's avatar
Gwenaelle Gremion committed
509
510
511
512
513
514
   write(*,900) '               betap_dph   = ',betap_dph
   write(10,901) betap_dph
   write(*,900) '               betap_dzo   = ',betap_dzo
   write(10,901) betap_dzo
   write(*,900) '               betap_fp   = ',betap_fp
   write(10,901) betap_fp
515
516
   write(*,900) '               betadph_dph   = ',betadph_dph
   write(10,901) betadph_dph
Gwenaelle Gremion's avatar
Gwenaelle Gremion committed
517
518
519
520
   write(*,900) '               betadph_dzo   = ',betadph_dzo
   write(10,901) betadph_dzo
   write(*,900) '               betadph_fp   = ',betadph_fp
   write(10,901) betadph_fp
521
   write(*,900) '               betadzo_dzo   = ',betadzo_dzo
Gwenaelle Gremion's avatar
Gwenaelle Gremion committed
522
523
524
   write(10,901) betadzo_dzo
   write(*,900) '               betadzo_fp   = ',betadzo_fp
   write(10,901) betadzo_fp
525
   write(*,900) '               betafp_fp   = ',betafp_fp
Gwenaelle Gremion's avatar
Gwenaelle Gremion committed
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
   write(10,901) betafp_fp
   write(*,900) '               betaBr   = ',betaBr
   write(10,901) betaBr
   write(*,900) '               betaSh   = ',betaSh
   write(10,901) betaSh
   write(*,900) '               betaDs   = ',betaDs
   write(10,901) betaDs
   write(*,900) '               eps_const   = ',eps_const
   write(10,901)eps_const
   write(*,900) '               eps_n   = ',eps_n
   write(10,901)eps_n
   write(*,900) '               cons_min   = ',cons_min
   write(10,901)cons_min
   write(*,900) '              sti_cst   = ',sti_cst
   write(10,901)sti_cst
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
   write(*,900) '               stip_p   = ',stip_p
   write(10,901) stip_p
   write(*,900) '               stip_dph   = ',stip_dph
   write(10,901) stip_dph
   write(*,900) '               stip_dzo   = ',stip_dzo
   write(10,901) stip_dzo
   write(*,900) '               stip_fp   = ',stip_fp
   write(10,901) stip_fp
   write(*,900) '               stidph_dph   = ',stidph_dph
   write(10,901) stidph_dph
   write(*,900) '               stidph_dzo   = ',stidph_dzo
   write(10,901) stidph_dzo
   write(*,900) '               stidph_fp   = ',stidph_fp
   write(10,901) stidph_fp
   write(*,900) '               stidzo_dzo   = ',stidzo_dzo
   write(10,901) stidzo_dzo
   write(*,900) '               stidzo_fp   = ',stidzo_fp
   write(10,901) stidzo_fp
   write(*,900) '               stifp_fp   = ',stifp_fp
Gwenaelle Gremion's avatar
Gwenaelle Gremion committed
560
561
562
563
564
565
   write(10,901) stifp_fp
!-----  Settling
   write(*,900) '               CSF   = ',CSF
   write(10,901) CSF
   write(*,900) '               size_rand   = ',size_rand
   write(10,901) size_rand
566
567
   write(*,900) '               size_phy_us   = ',size_phy_us
   write(10,901) size_phy_us
Gwenaelle Gremion's avatar
Gwenaelle Gremion committed
568
569
   write(*,900) '               size_phy_up   = ',size_phy_up
   write(10,901) size_phy_up
570
   write(*,900) '               size_dph_us   = ',size_dph_us
Gwenaelle Gremion's avatar
Gwenaelle Gremion committed
571
572
573
   write(10,901) size_dph_us
   write(*,900) '               size_dph_up   = ',size_dph_up
   write(10,901) size_dph_up
574
   write(*,900) '               size_dzo_us   = ',size_dzo_us
Gwenaelle Gremion's avatar
Gwenaelle Gremion committed
575
576
577
   write(10,901) size_dzo_us
   write(*,900) '               size_dzo_up   = ',size_dzo_up
   write(10,901) size_dzo_up
578
   write(*,900) '               size_fp_us  = ',size_fp_us
Gwenaelle Gremion's avatar
Gwenaelle Gremion committed
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
   write(10,901) size_fp_us
   write(*,900) '               size_fp_up  = ',size_fp_up
   write(10,901) size_fp_up
   write(*,900) '               dm_msn  = ',dm_msn
   write(10,901) dm_msn
   write(*,900) '               coef1  = ',coef1 
   write(10,901) coef1 
   write(*,900) '               coef2  = ',coef2
   write(10,901) coef2
   write(*,900) '               diam_msn_us  = ',diam_msn_us
   write(10,901) diam_msn_us
   write(*,900) '               coef3  = ',coef3
   write(10,901) coef3
   write(*,900) '               cons_max  = ',cons_max
   write(10,901) cons_max
594
595
596
597
598
   write(*,900) '               dynvis   = ',dynvis
   write(10,901) dynvis
   write(*,900) '               kinvis   = ',kinvis
   write(10,901) kinvis
   write(*,900) '               kB   = ',kB
Gwenaelle Gremion's avatar
Gwenaelle Gremion committed
599
600
   write(10,901) kB
!-----Fragmentation
601
602
603
604
   write(*,900) '               Frag_meth   = ',Frag_meth
   write(10,901) Frag_meth
   write(*,900) '                swim_brk   = ', swim_brk
   write(10,901)  swim_brk
Gwenaelle Gremion's avatar
Gwenaelle Gremion committed
605
606
607
608
   write(*,900) '               Floc_coef   = ', Floc_coef
   write(10,901)  Floc_coef
   write(*,900) '                Daugther_part     = ', Daugther_part  
   write(10,901)  Daugther_part  
609
610
611
612

900 format (a,f8.5)
901 format (f8.5)

Gwenaelle Gremion's avatar
Gwenaelle Gremion committed
613
614
615
616
617
!  Conversion from day to second (Days in the nml but the model go in seconds !)
 
!-----INITIAL and Minimum concentration for the variables
   mu5      = mu5 /secs_pr_day  
!-----Phytoplankton and dead-phytoplankton    
618
619
620
   vp       = vp /secs_pr_day
   alpha    = alpha /secs_pr_day                
   inib     = inib /secs_pr_day         
Gwenaelle Gremion's avatar
Gwenaelle Gremion committed
621
622
623
624
   mu1      = mu1 /secs_pr_day
   txloss_p = txloss_p /secs_pr_day 
  txloss_dph= txloss_dph /secs_pr_day 
!-----Zooplankton
625
   gmax     = gmax /secs_pr_day
Gwenaelle Gremion's avatar
Gwenaelle Gremion committed
626
627
628
629
630
   mu2      = mu2 /secs_pr_day
txloss_dzo  = txloss_dph /secs_pr_day  
txloss_fp   = txloss_fp /secs_pr_day 
  w_zmax    = w_zmax /secs_pr_day      
!-----Bacterias and LDON
631
   vb       = vb /secs_pr_day
Gwenaelle Gremion's avatar
Gwenaelle Gremion committed
632
  remi      = remi/secs_pr_day
633
634
   mu3      = mu3 /secs_pr_day
   mbac     = mbac /secs_pr_day       
Gwenaelle Gremion's avatar
Gwenaelle Gremion committed
635
  dphlossb  = dphlossb/secs_pr_day 
636
637
638
  dphlossl  = dphlossl/secs_pr_day 
  dzolossl  = dzolossl/secs_pr_day 
  dzolossb  = dzolossb/secs_pr_day 
Gwenaelle Gremion's avatar
Gwenaelle Gremion committed
639
640
641
642
643
  fplossl   = fplossl/secs_pr_day
  fplossb   = fplossb/secs_pr_day
  leak      = leak/secs_pr_day 
  mldon     = mldon /secs_pr_day       
!----- Detritus/Settling
644
645
646
647
   w_p      = w_p /secs_pr_day
   w_dph    = w_dph /secs_pr_day    
   w_dzo    = w_dzo /secs_pr_day    
   w_fp     = w_fp /secs_pr_day    
Gwenaelle Gremion's avatar
Gwenaelle Gremion committed
648
649
650
651
652
653
654
655
656
657
   w_msn    = w_msn /secs_pr_day    
!----- Coagulation  
betap_p     = betap_p /secs_pr_day    
betap_dph   = betap_dph /secs_pr_day   
betap_dzo   = betap_dzo /secs_pr_day   
betap_fp    = betap_fp /secs_pr_day    
betadph_dph = betadph_dph/secs_pr_day 
betadph_dzo = betadph_dzo/secs_pr_day 
betadph_fp  = betadph_fp/secs_pr_day 
betadzo_dzo = betadzo_dzo /secs_pr_day    
658
betadzo_fp  = betadzo_fp /secs_pr_day    
Gwenaelle Gremion's avatar
Gwenaelle Gremion committed
659
660
661
betafp_fp   = betafp_fp /secs_pr_day 
!-----Fragmentation
swim_brk   =  swim_brk /secs_pr_day
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693


   out_unit=unit

   LEVEL3 'polynow bio module initialised ...'

   return

98 LEVEL2 'I could not open bio_polynow.nml'
   LEVEL2 'If thats not what you want you have to supply bio_polynow.nml'
   LEVEL2 'See the bio example on www.gotm.net for a working bio_polynow.nml'
   return
99 FATAL 'I could not read bio_polynow.nml'
   stop 'init_bio_polynow'
   end subroutine init_bio_polynow
!EOC

!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: Initialise the concentration variables
!
! !INTERFACE:
   subroutine init_var_polynow(nlev)
!
! !DESCRIPTION:
!  Here, the the initial conditions are set and the settling velocities are
!  transferred to all vertical levels. All concentrations are declared
!  as non-negative variables, and it is defined which variables would be
!  taken up by benthic filter feeders.
!
! !USES:
Gwenaelle Gremion's avatar
Gwenaelle Gremion committed
694
695
   use observations,    only: nprof,aprof              
   use meanflow,        only: nit,amm,T,S               
696
697
698
699
700
701
702
703
704
705
706
 
   IMPLICIT NONE
!
! !INPUT PARAMETERS:
   integer, intent(in)                 :: nlev
!
! !REVISION HISTORY:
!  Original author(s): Hans Burchard & Karsten Bolding

! !LOCAL VARIABLES:
  integer                    :: i,j
Gwenaelle Gremion's avatar
Gwenaelle Gremion committed
707
  REALTYPE                   :: gauss_p,gauss_dph,gauss_dzo,gauss_fp,gauss_msn,gauss_zoo              
708
709
710
711
712
713
714
715
716
717
718

!EOP
!-----------------------------------------------------------------------
!BOC

!-----------------------------------------------------------------------------------
!-----------------------------------------------------------------------------------
!   Variable INITIAL PROFILES
!-----------------------------------------------------------------------------------
!-----------------------------------------------------------------------------------

Gwenaelle Gremion's avatar
Gwenaelle Gremion committed
719
   do i=1,nlev
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754

!-----------------------------------------------------------------------------------
!   Phytoplankton initial profile
!-----------------------------------------------------------------------------------

!!       ---If p_init_value .eq.1 --» Cst initial value for all the depth as mentionned in the nml.
if (p_init_value .eq. 1.0) then
      cc(p,i)=p_initial
!!       ---If p_init_value .eq.0 --» Initial profile describe according to gaussian curve (Mean + standard error)
else
  gauss_p=p_initial*exp(-((i-z_p_gauss_init)**2/(2*(sigma_p**2))))
     cc(p,i)= gauss_p
endif

!-----------------------------------------------------------------------------------
!   Zooplankton initial profiles
!-----------------------------------------------------------------------------------

!       ---If zoo_init_value .eq.1 --» Cst initial value for all the depth as mentionned in the nml.
if (zoo_init_value .eq. 1.0) then
      cc(z,i)=z_initial
!       ---If zoo_init_value.eq.0 --» Initial profile describe according to gaussian curve (Mean + standard error)
else
  gauss_zoo=z_initial*exp(-((i-z_zoo_gauss_init)**2/(2*(sigma_zoo**2))))
     cc(z,i)= gauss_zoo
endif
!-----------------------------------------------------------------------------------
!   Dead- Phytoplankton initial profile
!-----------------------------------------------------------------------------------

!!       ---If dph_init_value .eq.1 --» Cst initial value for all the depth as mentionned in the nml.
if (dph_init_value .eq. 1.0) then
      cc(d1,i)=dph_initial
!!       ---If dph_init_value .eq.0 --» Initial profile describe according to gaussian curve (Mean + standard error)
else
Gwenaelle Gremion's avatar
Gwenaelle Gremion committed
755
756
  gauss_dph=dph_initial*exp(-((i-z_dph_gauss_init)**2/(2*(sigma_dph**2))))
   cc(d1,i)= gauss_dph
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
endif

!-----------------------------------------------------------------------------------
!   Dead- Zooplankton initial profile
!-----------------------------------------------------------------------------------

!!       ---If dzo_init_value .eq.1 --» Cst initial value for all the depth as mentionned in the nml.
if (dzo_init_value .eq. 1.0) then
      cc(d2,i)=dzo_initial
!!       ---If dzo_init_value .eq.0 --» Initial profile describe according to gaussian curve (Mean + standard error)
else
  gauss_dzo=dzo_initial*exp(-((i-z_dzo_gauss_init)**2/(2*(sigma_dzo**2))))
     cc(d2,i)= gauss_dzo
endif

!-----------------------------------------------------------------------------------
!   Fecal pellets initial profile
!-----------------------------------------------------------------------------------

!!       ---If fp_init_value .eq.1 --» Cst initial value for all the depth as mentionned in the nml.
if (fp_init_value .eq. 1.0) then
      cc(d3,i)=fp_initial
!!       ---If fp_init_value .eq.0 --» Initial profile describe according to gaussian curve (Mean + standard error)
else
  gauss_fp=fp_initial*exp(-((i-z_fp_gauss_init)**2/(2*(sigma_fp**2))))
     cc(d3,i)= gauss_fp
endif

!-----------------------------------------------------------------------------------
!   Marine snow initial profile
!-----------------------------------------------------------------------------------

!!       ---If msn_init_value .eq.1 --» Cst initial value for all the depth as mentionned in the nml.
if (msn_init_value .eq. 1.0) then
      cc(d4,i)=msn_initial
!!       ---If msn_init_value.eq.0 --» Initial profile describe according to gaussian curve (Mean + standard error)
else
  gauss_msn=msn_initial*exp(-((i-z_msn_gauss_init)**2/(2*(sigma_msn**2))))
     cc(d4,i)= gauss_msn
endif


!-----------------------------------------------------------------------------------
!  Bacteria- Nitrate/ammonium / LDON initial profiles
!-----------------------------------------------------------------------------------
      cc(b,i)=b_initial
      cc(n,i)=nprof(i)       
      cc(a,i)=aprof(i)       
      cc(l,i)=l_initial

   end do

!-----------------------------------------------------------------------------------
!   Settling/swimming velocity
!-----------------------------------------------------------------------------------

   do i=0,nlev
      ws(z,i) = _ZERO_
      ws(b,i) = _ZERO_
      ws(n,i) = _ZERO_
      ws(a,i) = _ZERO_
      ws(l,i) = _ZERO_

if(Phys_w .eq. 0.0) then    
      ws(p,i) = w_p
      ws(d1,i)= w_dph   
      ws(d2,i)= w_dzo   
      ws(d3,i)= w_fp    
      ws(d4,i)= w_msn 
else
     ws(p,i) =_ZERO_           
     ws(d1,i)=_ZERO_
     ws(d2,i)=_ZERO_ 
     ws(d3,i)=_ZERO_ 
     ws(d4,i)=_ZERO_ 

end if

!-----------------------------------------------------------------------------------
!   Marine Snow Size and settling velocity
!-----------------------------------------------------------------------------------
Gwenaelle Gremion's avatar
Gwenaelle Gremion committed
838
839
840
  cc(size_intrm,i) =_ZERO_
  cc(size_coag,i)  =_ZERO_
  cc(size_frag,i)  =_ZERO_
841

Gwenaelle Gremion's avatar
Gwenaelle Gremion committed
842
843
  cc(size_msn,i)=_ZERO_
  ws(size_msn,i)=_ZERO_
844

845
846
847
848
849
850
851
852
853
854

   end do

   posconc(p) = 1
   posconc(z) = 1
   posconc(b) = 1
   posconc(n) = 1
   posconc(a) = 1
   posconc(l) = 1
  
Gwenaelle Gremion's avatar
Gwenaelle Gremion committed
855
856
857
858
  posconc(d1) = 1 
  posconc(d2) = 1 
  posconc(d3) = 1 
  posconc(d4) = 1 
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945

#if 0
   mussels_inhale(p) = .true.
   mussels_inhale(z) = .true.
   mussels_inhale(b) = .false.
   mussels_inhale(n) = .false.
   mussels_inhale(a) = .false.
   mussels_inhale(l) = .true.
 
   mussels_inhale(d1) = .true.  
   mussels_inhale(d2) = .true.  
   mussels_inhale(d3) = .true.  
   mussels_inhale(d4) = .true.  

#endif

   LEVEL3 'polynow variables initialised ...'

   return

   end subroutine init_var_polynow
!EOC

!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: Providing info on variables
!
! !INTERFACE:
   subroutine var_info_polynow()
!
! !DESCRIPTION:
!  This subroutine provides information about the variable names as they
!  will be used when storing data in NetCDF files.
!
! !USES:
   IMPLICIT NONE
!
! !REVISION HISTORY:
!  Original author(s): Hans Burchard & Karsten Bolding
!
! !LOCAL VARIABLES:
!EOP
!-----------------------------------------------------------------------
!BOC
   var_names(1) = 'phy'
   var_units(1) = 'mmol/m**3'
   var_long(1)  = 'phytoplankton'

   var_names(2) = 'zoo'
   var_units(2) = 'mmol/m**3'
   var_long(2)  = 'zooplankton'

   var_names(3) = 'bac'
   var_units(3) = 'mmol/m**3'
   var_long(3)  = 'bacteria'

   var_names(4) = 'dph'                
   var_units(4) = 'mmol/m**3'
   var_long(4)  = 'dead phytoplankton'

   var_names(5) = 'nit'
   var_units(5) = 'mmol/m**3'
   var_long(5)  = 'nitrate'

   var_names(6) = 'amm'
   var_units(6) = 'mmol/m**3'
   var_long(6)  = 'ammonium'

   var_names(7) = 'ldn'
   var_units(7) = 'mmol/m**3'
   var_long(7)  = 'labile_dissolved_organic_nitrogen'

   var_names(8) = 'dzo'                            
   var_units(8) = 'mmol/m**3'                      
   var_long(8)  = 'dead zooplankton'           

   var_names(9) = 'fp'                              
   var_units(9) = 'mmol/m**3'
   var_long(9)  = 'fecal pellets'

   var_names(10) = 'msn'                          
   var_units(10) = 'mmol/m**3'
   var_long(10)  = 'marine snow'

!--------------------------------------------------

Gwenaelle Gremion's avatar
Gwenaelle Gremion committed
946
   var_names(11) = 'Size_intrm'                          
947
   var_units(11) = 'm'
948
   var_long(11)  = 'Intermediate size of marine snow'
949
950


Gwenaelle Gremion's avatar
Gwenaelle Gremion committed
951
   var_names(12) = 'size_coag'                          
952
953
954
955
   var_units(12) = 'm'
   var_long(12)  = 'Size augmentation of msn due to Coag.'


Gwenaelle Gremion's avatar
Gwenaelle Gremion committed
956
   var_names(13) = 'size_frag'                          
957
958
959
960
   var_units(13) = 'm'
   var_long(13)  = 'Size of msn due to Frag.'


Gwenaelle Gremion's avatar
Gwenaelle Gremion committed
961
   var_names(14) = 'size_msn'                          
962
963
   var_units(14) = 'm'
   var_long(14)  = 'Final size of marine snow'
964

965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
1001
1002
1003
1004
1005
1006
1007
1008
1009
1010
1011
1012
1013
1014
1015
1016
1017
1018
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034
1035

   return
   end subroutine var_info_polynow
!EOC

!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: Light properties for the Fasham model
!
! !INTERFACE:
   subroutine light_polynow(nlev,bioshade_feedback)
!
! !DESCRIPTION:
! Here, the photosynthetically available radiation is calculated
! by simply assuming that the short wave part of the total
! radiation is available for photosynthesis. 
! The photosynthetically
! available radiation, $I_{PAR}$, follows from (\ref{light}).
! The user should make
! sure that this is consistent with the light class given in the
! {\tt extinct} namelist of the {\tt obs.nml} file.
! The self-shading effect is also calculated in this subroutine,
! which may be used to consider the effect of bio-turbidity also
! in the temperature equation (if {\tt bioshade\_feedback} is set
! to true in {\tt bio.nml}).
! For details, see section \ref{sec:do-bio}.
!
! !USES:
   IMPLICIT NONE
!
! !INPUT PARAMETERS:
  integer, intent(in)                  :: nlev
  logical, intent(in)                  :: bioshade_feedback
!
! !REVISION HISTORY:
!  Original author(s): Hans Burchard, Karsten Bolding
!
! !LOCAL VARIABLES:
   integer                   :: i
   REALTYPE                  :: zz,add
!EOP
!-----------------------------------------------------------------------
!BOC
   zz = _ZERO_
   add = _ZERO_
   do i=nlev,1,-1
      add=add+0.5*h(i)*(cc(p,i)+p0)
      zz=zz+0.5*h(i)
      par(i)=rad(nlev)*(1.-aa)*exp(-zz/g2)*exp(-kc*add)
      add=add+0.5*h(i)*(cc(p,i)+p0)
      zz=zz+0.5*h(i)
      if (bioshade_feedback) bioshade_(i)=exp(-kc*add)
   end do


   return
   end subroutine light_polynow
!EOC

!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: Right hand sides of geobiochemical model \label{sec:bio-fasham-rhs}
!
! !INTERFACE:
   subroutine do_bio_polynow(first,numc,nlev,cc,pp,dd)
!
! !DESCRIPTION:
! The model consisting of the following state variable:
! phytoplankton, bacteria, detritus (4 types), zooplankton, 
Gwenaelle Gremion's avatar
Gwenaelle Gremion committed
1036
! nitrate, ammonium and dissolved organic nitrogen as well as traits link to the size evolution of marine snow (4)
1037
!
Gwenaelle Gremion's avatar
Gwenaelle Gremion committed
1038
1039
1040
1041
1042
1043
1044
1045
1046
!-----BIBLIOGRAPHY:
! Platt et al., 1980
! Moore et al., 2002
! Alldredge and Gotschalk 1988
! Kajihara 1971
! Komar 1981 & 1978 
! Ghosh 2013 
! Mc Donnell 2010 
! Alldredge 1990
1047
1048
1049
! Biddanda 1998
! Ploug and Grossart 2000
! Jokulsdottir 2011
Gwenaelle Gremion's avatar
Gwenaelle Gremion committed
1050
1051
1052
1053
1054
! Jackson 1990
! Bagheri 2015
! VanRijn 1993
! Spicer et al 1996
! Dilling 2000
1055
!
Gwenaelle Gremion's avatar
Gwenaelle Gremion committed
1056
1057
1058
! !USES: Force the user to declare all the variables 
!The variables who do not have declaration will follow the rules : If the variables start with a letters from I to N will be integer, otherwise they will be reals. 
  IMPLICIT NONE
1059
1060

! !INPUT PARAMETERS:
Gwenaelle Gremion's avatar
Gwenaelle Gremion committed
1061
1062
1063
1064
1065
   logical, intent(in)        :: first
   integer, intent(in)        :: numc,nlev

! -- Concentrations
   REALTYPE, intent(out)      :: cc(1:numc,0:nlev)
1066
1067
1068
!With INTENT(IN), you are promising to the compiler that argument cc is not modified nor its definition status changed. 

! !OUTPUT PARAMETERS:
Gwenaelle Gremion's avatar
Gwenaelle Gremion committed
1069
1070
1071
   REALTYPE, intent(out)      :: pp(1:numc,1:numc,0:nlev)
! -- Fluxes
   REALTYPE, intent(out)      :: dd(1:numc,1:numc,0:nlev)
1072
1073
1074
1075
1076
1077
1078
1079
1080
!Similarly, INTENT(OUT) says that variables pp and dd are always assigned a value in the routine.

! !LOCAL VARIABLES:
   integer                    :: i,j,ci
!   Light and nutrient limitation factors & Zooplankton and bacterias needed factors
   REALTYPE                   :: Ps,ff
   REALTYPE                   :: q1,q2
   REALTYPE                   :: fac,min67
!!   Size variation for each particle type
Gwenaelle Gremion's avatar
Gwenaelle Gremion committed
1081
1082
   REALTYPE                   ::size_phy,size_dph,size_dzo,size_fp
   REALTYPE                   ::size_phy_max,size_dph_max,size_dzo_max,size_fp_max
1083
1084
1085
1086
   REALTYPE                   ::diam_phy,diam_dph,diam_dzo,diam_fp,diam_msn 
!Collision & Aggregation
!!   Sedimentation rate for each particle type
     REALTYPE                    ::w_p_m,w_dph_m,w_dzo_m,w_fp_m
Gwenaelle Gremion's avatar
Gwenaelle Gremion committed
1087
     REALTYPE                    ::densFlu
1088
     REALTYPE                    ::Rp,Rdph,Rdzo,Rfp
1089
      REALTYPE                   ::Re_phy,Re_dph,Re_dzo,Re_fp
Gwenaelle Gremion's avatar
Gwenaelle Gremion committed
1090
1091
     REALTYPE                    ::CSF_phy,CSF_dph,CSF_dzo,CSF_fp,CSF_msn
   REALTYPE                    ::CSF_phy_calc,CSF_dph_calc,CSF_dzo_calc,CSF_fp_calc,CSF_msn_calc
1092
1093
1094
1095
1096
1097
1098
!!   Maximum loss for aggregation for each particle type
      REALTYPE                    ::pprime   
      REALTYPE                    ::d1rime  
      REALTYPE                    ::d2rime
      REALTYPE                    ::d3rime
!!   Stickiness for each particle type interaction
      REAL                        :: r_phy, r_dph, r_dzo, r_fp
Gwenaelle Gremion's avatar
Gwenaelle Gremion committed
1099
      REALTYPE                    :: stip ,stidph,stidzo,stifp
1100
1101
1102
1103
1104
1105
1106
1107
1108
1109
1110
1111
1112
1113
1114
1115
1116
1117
1118
1119
1120
1121
1122
1123
1124
1125
1126
1127
1128
1129
1130
1131
1132
1133
1134
1135
1136
1137
1138
1139
1140
1141
1142
1143
1144
1145
1146
1147
1148
1149
1150
1151
1152
1153
1154
1155
1156
1157
1158
1159
1160
1161
1162
1163
1164
1165
1166
1167
1168
1169
1170
1171
1172
1173
1174
1175
1176
1177
1178
1179
1180
1181
1182
!'No physical coagulation happens here [Constant Beta are set]'
! Phytoplankton loss by aggregation
      REALTYPE                    ::coli_pp_NoPhys,agg_pp
      REALTYPE                    ::coli_pdph_NoPhys,agg_pdph
      REALTYPE                    ::coli_pdzo_NoPhys,agg_pdzo
      REALTYPE                    ::coli_pfp_NoPhys,agg_pfp
! Dead Phytoplankton loss by aggregation
      REALTYPE                    ::coli_dphdph_NoPhys,agg_dphdph
      REALTYPE                    ::coli_dphdzo_NoPhys,agg_dphdzo
      REALTYPE                    :: coli_dphfp_NoPhys,agg_dphfp
! Dead Zooplankton loss by aggregation
      REALTYPE                    ::coli_dzodzo_NoPhys,agg_dzodzo
      REALTYPE                    ::coli_dzofp_NoPhys,agg_dzofp
! Fecal pellets loss by aggregation
      REALTYPE                    ::coli_fpfp_NoPhys,agg_fpfp
!Coagulate only with : BetaBR'
      REALTYPE                    ::T_degK
      REALTYPE                    ::physic_Br
  !Phytoplankton
      REALTYPE                    :: betaBr_pp,coli_pp
      REALTYPE                    :: betaBr_pdph,coli_pdph
      REALTYPE                    :: betaBr_pdzo,coli_pdzo
      REALTYPE                    :: betaBr_pfp,coli_pfp
  !Dead Phytoplankton
      REALTYPE                    :: betaBr_dphdph,coli_dphdph
      REALTYPE                    :: betaBr_dphdzo,coli_dphdzo
      REALTYPE                    :: betaBr_dphfp,coli_dphfp
  !Dead Zooplankton
      REALTYPE                    :: betaBr_dzodzo,coli_dzodzo
      REALTYPE                    :: betaBr_dzofp,coli_dzofp
  !Fecal pellets
      REALTYPE                    :: betaBr_fpfp,coli_fpfp
!Coagulate only with : BetaSh'
      REALTYPE                    ::bio_eps
  !Phytoplankton
      REALTYPE                    ::q_pp,betaSh_pp
      REALTYPE                    ::q_pdph,betaSh_pdph
      REALTYPE                    ::q_pdzo,betaSh_pdzo
      REALTYPE                    ::q_pfp,betaSh_pfp
  !Dead Phytoplankton
      REALTYPE                    ::q_dphdph,betaSh_dphdph
      REALTYPE                    ::q_dphdzo,betaSh_dphdzo
      REALTYPE                    ::q_dphfp,betaSh_dphfp
  !Dead Zooplankton
      REALTYPE                    ::q_dzodzo,betaSh_dzodzo
      REALTYPE                    ::q_dzofp,betaSh_dzofp
  !Fecal pellets
      REALTYPE                    ::q_fpfp,betaSh_fpfp
!Coagulate only with : BetaDs'
  !Phytoplankton
      REALTYPE                    ::abs_pp, betaDs_pp
      REALTYPE                    ::abs_pdph,betaDs_pdph
      REALTYPE                    ::abs_pdzo,betaDs_pdzo
      REALTYPE                    ::abs_pfp,betaDs_pfp
  !Dead Phytoplankton
      REALTYPE                    :: abs_dphdph,betaDs_dphdph
      REALTYPE                    :: abs_dphdzo,betaDs_dphdzo
      REALTYPE                    :: abs_dphfp,betaDs_dphfp
  !Dead Zooplankton
      REALTYPE                    :: abs_dzodzo,betaDs_dzodzo
      REALTYPE                    :: abs_dzofp,betaDs_dzofp
  !Fecal pellets
      REALTYPE                    :: abs_fpfp,betaDs_fpfp
!'Physical Coagulation happens here ![Beta are calculated via physic]'
  !Phytoplankton
      REALTYPE                    ::kernel_coef_pp
      REALTYPE                    ::kernel_coef_pdph
      REALTYPE                    ::kernel_coef_pdzo
      REALTYPE                    ::kernel_coef_pfp
  !Dead Phytoplankton
      REALTYPE                    :: kernel_coef_dphdph
      REALTYPE                    :: kernel_coef_dphdzo
      REALTYPE                    :: kernel_coef_dphfp
  !Dead Zooplankton
      REALTYPE                    :: kernel_coef_dzodzo
      REALTYPE                    :: kernel_coef_dzofp
  !Fecal pellets
      REALTYPE                    :: kernel_coef_fpfp
!SIZE Marine Snow with aggregation
      REALTYPE                    :: Ratio_P,Ratio_D1,Ratio_D2,Ratio_D3
      REALTYPE                    :: size_msn_m
      REALTYPE                    :: max_flux
!FRAGMENTATION
Gwenaelle Gremion's avatar
Gwenaelle Gremion committed
1183
1184
1185
1186
      REALTYPE                    :: Frag_bio,Frag_phys,leackage,remineralization
      REALTYPE                    :: kolmog,stable_size
      REALTYPE                    :: prob_break_msn,r_msn
      REALTYPE                    :: diam_msn_max,min_size_msn
1187
!   Sedimentation rate msn
Gwenaelle Gremion's avatar
Gwenaelle Gremion committed
1188
      REALTYPE                    :: w_msn_m,Re_msn,Rmsn,CFL
1189
1190

REALTYPE :: ceci
1191
1192
1193
1194
1195
1196
1197
1198
1199

!EOP
!-----------------------------------------------------------------------
!BOC
!KBK - is it necessary to initialise every time - expensive in a 3D model
   pp = _ZERO_
   dd = _ZERO_

   do ci=1,nlev   !#A 
Gwenaelle Gremion's avatar
Gwenaelle Gremion committed
1200

1201
1202
1203
1204
1205
1206
1207
1208
1209
1210
1211
1212
1213
1214
1215
1216
1217
1218
1219
1220
1221
1222
1223
1224
1225
1226
1227
!-----------------------------------------------------------------------------------
!   Light and nutrient limitation factors & Zooplankton and bacterias needed factors
!-----------------------------------------------------------------------------------

!!       ---Light : Platt et al. (1980) - inhibition
      Ps= vp/((alpha/(alpha+inib))*(alpha/(alpha+inib))**(inib/alpha))
      ff= Ps*(1.-exp(-1.*alpha*par(ci)/Ps))*exp(-1.*inib*par(ci)/Ps)
      lumlim1(ci)=ff
!!        ---Nitrate :
      q1=(cc(n,ci)/k1)/(1.+cc(n,ci)/k1+cc(a,ci)/k2)
      nitlim1(ci)=q1
!!        ---Ammonium :
      q2=(cc(a,ci)/k2)/(1.+cc(n,ci)/k1+cc(a,ci)/k2)
      ammlim1(ci)=q2
!!        ---Zooplankton grazing preference normalization factor  
      fac=(cc(z,ci)+z0)/(k3*(r1*cc(p,ci)+r2*cc(b,ci)+  &
             r3*cc(d1,ci)+r4*cc(d2,ci)+r5*cc(d3,ci)+r6*cc(d4,ci))+  &          
                      r1*cc(p,ci)**2+r2*cc(b,ci)**2+   &
             r3*cc(d1,ci)**2+r4*cc(d2,ci)**2+r5*cc(d3,ci)**2+r6*cc(d4,ci)**2)  
!!        --- ???
      min67=min(cc(a,ci),eta*cc(l,ci))  

!-----------------------------------------------------------------------------------
!   Nutrients 
!-----------------------------------------------------------------------------------

      dd(n,p,ci)=ff*q1*(cc(p,ci)+p0)
Gwenaelle Gremion's avatar
Gwenaelle Gremion committed
1228
1229
      dd(a,p,ci)=ff*q2*(cc(p,ci)+p0)

1230
1231
1232
1233
1234
1235
1236
1237
1238
1239
1240
1241
1242
1243
1244
1245
1246
1247
1248
1249
1250
1251
1252
1253
1254
1255
1256
1257
1258
1259
1260
1261
1262
1263
1264
1265
1266
1267
1268
1269
1270
1271
1272
1273
1274
1275
1276
1277
1278
1279
1280
1281
1282
1283
1284
1285
1286
1287
1288
1289
1290
1291
1292
1293
1294
1295
1296
1297
1298
1299
1300
1301
1302
1303
1304
1305
1306
1307
1308
1309
1310
1311
1312
1313
!!    ---Bacteria respiration [Bacteria biomass growth]
      dd(a,b,ci)=vb*min67/(k4+min67+cc(l,ci))*(cc(b,ci)+b0)
      dd(l,b,ci)=vb*cc(l,ci)/(k4+min67+cc(l,ci))*(cc(b,ci)+b0)

!!    ---Condition of [ldon] to allow ldon to become a matrix for the aggregate 
if (cc(l,ci) .ge. lmin) then        !#1                           
     dd(l,d4,ci)=mldon*cc(l,ci)                            
else                                !#1                       
     dd(l,d4,ci)= 0.0                                         
end if                              !#1                 
!!    ---Nitrification rate 
     dd(a,n,ci)=mu5*cc(a,ci)                                 
      
!-----------------------------------------------------------------------------------
!   Phytoplankton [LOSSES]
!-----------------------------------------------------------------------------------

!!    ---Grazing losses (if beta =1 zooplankton eat everything)
      dd(p,z,ci)=beta*gmax*r1*cc(p,ci)**2*fac

!!    ---Non grazing losses : Excretion
      dd(p,l,ci)=gamma*ff*(q1+q2)*cc(p,ci)

!!    ---Non grazing losses : 
!!            -viral lysis (not taken into acount here)
!!            -internal respiration degradation
!In Moore et al,.2002, all Non grazing losses are estimated at 10% for diatom and small phytoplankton.
     dd(p,d1,ci)=mu1*(cc(p,ci)+p0)/(k5+cc(p,ci)+p0)*cc(p,ci)  &  
                 +(1.-beta)*gmax*r1*cc(p,ci)**2*fac

!!    ---Aggregation losses (See section COLLISION & AGGREGATION PROCESSES)

!-----------------------------------------------------------------------------------
!   Zooplankton [LOSSES]
!-----------------------------------------------------------------------------------

!!    ---Excretion
       dd(z,a,ci)=epsi*mu2*(cc(z,ci)+z0)/(k6+cc(z,ci)+z0)*cc(z,ci)
       dd(z,l,ci)=delta*mu2*(cc(z,ci)+z0)/(k6+cc(z,ci)+z0)*cc(z,ci)
!!    ---Sloppy feeding (z to dph or d1) :Taken into account in dd(p,d1,ci) see above.
!!    ---Fecal pelets
      dd(z,d3,ci)=eg*mu2*(cc(z,ci)+z0)/(k6+cc(z,ci)+z0)*cc(z,ci)
!!    ---Dead Zoo
      dd(z,d2,ci)=(1.-epsi-delta-eg)  &                       
        *mu2*(cc(z,ci)+z0)/(k6+cc(z,ci)+z0)*cc(z,ci) 

!-----------------------------------------------------------------------------------
!   Zooplankton [Diurnal migration]
!               -----------
! As a function of light and [phytoplancton](food)
!-----------------------------------------------------------------------------------

  if (Migra_zoo.eq. 1.0) then  !#2
      if (cc(p,ci) .lt. pmin .or. cc(d1,ci) .lt. dph_initial )  then        !#3    
          ws(z,ci) = -1.0*w_zmax*tanh(bertha*(par(ci)-parcrit)) 
      else !#3
          ws(z,ci) = 0.0
      end if !#3
 else            !#2
  ws(z,ci) = 0.0
 end if          !#2

!-----------------------------------------------------------------------------------
!   Bacterias [LOSSES]
!-----------------------------------------------------------------------------------
!!    ---Bacteria remineralisation    
       dd(b,a,ci)=mu3*cc(b,ci)   
!!    ---Bacteria grazed by zooplankton     
       dd(b,z,ci)=beta*gmax*r2*cc(b,ci)**2*fac
!!    ---Bacteria in aggregate (D4)
       dd(b,d4,ci)=mbac*cc(b,ci)                        
!!    ---Bacteria to detritus-non grazing mortality (D1)
       dd(b,d1,ci)=(1.-beta)*gmax*r2*cc(b,ci)**2*fac   

!-----------------------------------------------------------------------------------
!   Detritus [LOSSES]
!-----------------------------------------------------------------------------------
 
!!    ---Exsudation (see Moore et al., 2002)
      dd(d1,l,ci)=dphlossl*cc(d1,ci)            
      dd(d2,l,ci)=dzolossl*cc(d2,ci)            
      dd(d3,l,ci)=fplossl*cc(d3,ci)             

!!    ---Grazing
Gwenaelle Gremion's avatar
Gwenaelle Gremion committed
1314
1315
1316
1317
      dd(d1,z,ci)=beta*gmax*r3*cc(d1,ci)**2*fac       !Grazing-Saprophagy
      dd(d2,z,ci)=beta*gmax*r4*cc(d2,ci)**2*fac       !Carnivory/cannibalism/necrophagy
      dd(d3,z,ci)=beta*gmax*r5*cc(d3,ci)**2*fac       !Copprophagy/scatophagy
      dd(d4,z,ci)=beta*gmax*r6*cc(d4,ci)**2*fac       !Grazing-Macrophagie
1318
1319
1320
1321
1322
1323
1324
1325
1326
1327
1328
1329
1330
1331
1332
1333
1334
1335
1336
1337
1338
1339
1340
1341
1342
1343
1344
1345
1346
1347
1348
1349


!-------------------------------------------------------------------------------------
!!   Size (radius in m) variation for each particle type 
!-------------------------------------------------------------------------------------

if (size_rand .eq. 1.0) then!#4

  call random_number(r_phy)
size_phy=size_phy_us+((size_phy_up-size_phy_us)*r_phy)
   call random_number(r_dph)
size_dph=size_dph_us+((size_dph_up-size_dph_us)*r_dph)
   call random_number(r_dzo)
size_dzo=size_dzo_us+((size_dzo_up-size_dzo_us)*r_dzo)
  call random_number(r_fp)
size_fp=size_fp_us+((size_fp_up-size_fp_us)*r_fp)
   
else  !#4
size_phy=size_phy_us
size_dph=size_dph_us
size_dzo=size_dzo_us
size_fp=size_fp_us

endif !#4


!-----------------------------------------------------------------------------------
!   COLLISION & AGGREGATION 
!-----------------------------------------------------------------------------------
 

!-------------------------------------------------------------------------------------
Gwenaelle Gremion's avatar
Gwenaelle Gremion committed
1350
!!   Usuefull parameters of Sedimentation rate for each particle type 
1351
1352
1353
1354
1355
1356
1357
1358
1359
1360
1361
1362
!------------------------------------------------------------------------------------- 

!----------Definition of the equivalent diameter of our particles [Bagheri 2015]-Unit of size [m]
!dS,i=racine3(LIS)   
!With L the longest dimension of the particle, I the longest dimension perpendicular to L, and S the longest dimension perpendicular to L and I (Bagheri 2015).
!With for phytoplankton and dead phytoplankton : L=I=S= size_X*2, and for dead zoo plankton and fecal pellets : L=I=size_x*4 and S= size_X*2 --» With size_X which represents the radius of our particle.
!We consider size_X the value of size that we can multiply to have the dimension we are looking for.
!Calculation of the cubic racin of X --» SIGN(ABS(X)**(1./3.),X) [FORTRAN, le langage normalisé, Michel Dubesset,Jean Vignes]

diam_phy = sign(abs(((size_phy*2.)*(size_phy*2.)*(size_phy*2.)))**(1./3.),((size_phy*2.)*(size_phy*2.)*(size_phy*2.)))

diam_dph = sign(abs(((size_dph*2.)*(size_dph*2.)*(size_dph*2.)))**(1./3.),((size_dph*2.)*(size_dph*2.)*(size_dph*2.))) 
1363

1364
1365
1366
diam_dzo = sign(abs(((size_dzo*4.)*(size_dzo*4.)*(size_dzo*2.)))**(1./3.),((size_dzo*4.)*(size_dzo*4.)*(size_dzo*2.)))

diam_fp  = sign(abs(((size_fp*4.)*(size_fp*4.)*(size_fp*2.)))**(1./3.),((size_fp*4.)*(size_fp*4.)*(size_fp*2.)))
1367

Gwenaelle Gremion's avatar
Gwenaelle Gremion committed
1368
1369
!--» Does CFL limit 
  CFL = (depth_bio/nlev)/dt_bio
1370
1371
1372
1373

!----------Calculation of the CSF factor
!--Phytoplankton
  CSF_phy = (size_phy*2)/sqrt((size_phy*2)*(size_phy*2))
Gwenaelle Gremion's avatar
Gwenaelle Gremion committed
1374
1375
1376
1377
1378
1379
1380
1381
1382
1383
       if(CSF_phy .ge. 0.0 .and. CSF_phy.lt. 0.4) then !#5.5
          CSF_phy_calc = 2.18-(2.09*CSF_phy)
       else if (CSF_phy .ge. 0.4 .and. CSF_phy .lt. 0.8) then !#5.6
          CSF_phy_calc = 0.946*(CSF_phy)**(-0.378)
       else if (CSF_phy .ge. 0.8 .and. CSF_phy .le. 1.0) then !#5.6
          CSF_phy_calc = 1.0 ! consider as a sphere
        else  !#5.6
          write (*,*) 'Error in CSF values for phytoplankton : CSF .NE. [0-1]', CSF_phy
       endif !#5.6

1384
!--Dead Phytoplankton
Gwenaelle Gremion's avatar
Gwenaelle Gremion committed
1385
1386
1387
1388
1389
1390
1391
1392
1393
1394
1395
1396
  CSF_dph = (size_dph*2)/sqrt((size_dph*2)*(size_dph*2))
write(*,*) 'CSF_dph',CSF_dph
 if(CSF_dph .ge. 0.0 .and. CSF_dph.lt. 0.4) then !#5.7
          CSF_dph_calc = 2.18-(2.09*CSF_dph)
      else if (CSF_dph .ge. 0.4 .and. CSF_dph .lt. 0.8) then !#5.7
           CSF_dph_calc = 0.946*(CSF_dph)**(-0.378)
      else if (CSF_dph .ge. 0.8 .and. CSF_dph .le. 1.0) then !#5.7
           CSF_dph_calc = 1.0 ! consider as a sphere
      else !#5.7
           write (*,*) 'Error in CSF values for Dead phytoplankton : CSF .NE. [0-1]',CSF_dph
      endif !#5.7

1397
!--Dead Zooplankton
Gwenaelle Gremion's avatar
Gwenaelle Gremion committed
1398
1399
1400
1401
1402
1403
1404
1405
1406
1407
1408
1409
1410
1411
1412
1413
1414
1415
1416
1417
1418
1419
1420
  CSF_dzo = (size_dzo*2)/sqrt((size_dzo*4)*(size_dzo*4))
   if(CSF_dzo .ge. 0.0 .and. CSF_dzo.lt. 0.4) then !#5.8
          CSF_dzo_calc = 2.18-(2.09*CSF_dzo)
      else if (CSF_dzo .ge. 0.4 .and. CSF_dzo .lt. 0.8) then !#5.8
           CSF_dzo_calc = 0.946*(CSF_dzo)**(-0.378)
      else if (CSF_dzo .ge. 0.8 .and. CSF_dzo .le. 1.0) then !#5.8
           CSF_dzo_calc = 1.0 ! consider as a sphere
      else  !#5.8
           write (*,*) 'Error in CSF values for dead zooplankton : CSF .NE. [0-1]',CSF_dzo
      endif !#5.8

!--Fecal Pellets
 CSF_fp = (size_fp*2)/sqrt((size_fp*4)*(size_fp*4))
write(*,*) 'CSF_fp',CSF_fp
      if(CSF_fp .ge. 0.0 .and. CSF_fp.lt. 0.4) then !#5.9
          CSF_fp_calc = 2.18-(2.09*CSF_fp)
      else if (CSF_fp .ge. 0.4 .and. CSF_fp .lt. 0.8) then !#5.9
           CSF_fp_calc = 0.946*(CSF_fp)**(-0.378)
      else if (CSF_fp .ge. 0.8 .and. CSF_fp .le. 1.0) then !#5.9
           CSF_fp_calc = 1.0 ! consider as a sphere
      else !#5.9
            write (*,*) 'Error in CSF values for fecal pellets : CSF .NE. [0-1]',CSF_fp
      endif !#5.9
1421

1422
1423
1424
1425
1426
1427
1428
1429

!----------Calcul of densities
!Fluid density
densFlu = rho_0
!Difference of densities
Rp = (rho_p-densFlu)
Rdph= (rho_dph-densFlu)
Rdzo =(rho_dzo-densFlu)
Gwenaelle Gremion's avatar
Gwenaelle Gremion committed
1430
1431
1432
1433
1434
1435
1436
1437
1438
1439
1440
1441
1442
1443
1444
1445
1446
1447
1448
1449
1450
1451
1452
1453
1454
1455
1456
1457
1458
1459
1460
1461
1462
1463
1464
Rfp =(rho_fp-densFlu)
Rmsn = (rho_msn-densFlu)

        !-----------------------------------------------------------------------------------
        !!    ---Safety check
        !-----------------------------------------------------------------------------------

!----------Estimation of the maximum size allowed for particles
 !--» Comparaison with CFL limit 
    size_phy_max = ((CFL/0.93)**2)*(1/(g*(Rp/densFlu)))
    size_dph_max = ((CFL/0.93)**2)*(1/(g*(Rdph/densFlu)))
    size_dzo_max = ((CFL/0.93)**2)*(1/(g*(Rdzo/densFlu)))
    size_fp_max  = ((CFL/0.93)**2)*(1/(g*(Rfp/densFlu)))

!---------- Attribution of the good size
if (size_phy .gt. size_phy_max) then
  size_phy = size_phy_max
   write(*,*) 'Size of Phy OVERPASSES Max.Size of Phy to respect CFL limit'
endif
if (size_dph .gt. size_dph_max) then
  size_dph = size_dph_max
   write(*,*) 'Size of Dph OVERPASSES Max.Size of Dph to respect CFL limit'
endif
if (size_dzo .gt. size_dzo_max) then
  size_dzo = size_dzo_max
   write(*,*) 'Size of Dzo OVERPASSES Max.Size of Dzo to respect CFL limit'
endif
if (size_fp .gt. size_fp_max) then
  size_fp = size_fp_max
   write(*,*) 'Size of Fp OVERPASSES Max.Size of FP to respect CFL limit'
endif

!-------------------------------------------------------------------------------------
!!   Calculation of Sedimentation rates 
!------------------------------------------------------------------------------------- 
1465
1466
1467
1468

!!--------!Calculated by value from the .nml
if (Phys_w .eq. 0.0 ) then  !#5

Gwenaelle Gremion's avatar
Gwenaelle Gremion committed
1469
1470
1471
1472
      w_p_m   = w_p   ! [m/s]
      w_dph_m = w_dph ! [m/s]
      w_dzo_m = w_dzo ! [m/s]
      w_fp_m  = w_fp  ! [m/s]
1473
1474
1475
1476
1477
1478

!!--------Stokes's law from Ghosh et al 2013
!     all particles are considered as sphere,so the value of diameter needed in the equation = radius*2                     
else if (Phys_w .eq. 1.0) then  !#5

      !--Phytoplankton
Gwenaelle Gremion's avatar
Gwenaelle Gremion committed
1479
      w_p_m=(g*((2*size_phy)**2)*Rp)/(18*dynvis)  ! [m/s]
1480
1481

      !--Dead phytoplankton
Gwenaelle Gremion's avatar
Gwenaelle Gremion committed
1482
1483
      w_dph_m=(g*((2*size_dph)**2)*Rdph)/(18*dynvis)  ! [m/s]
 
1484
      !--Dead Zooplankton
Gwenaelle Gremion's avatar
Gwenaelle Gremion committed
1485
1486
      w_dzo_m=(g*((2*size_dzo)**2)*Rdzo)/(18*dynvis) ! [m/s]
    
1487
      !--Fecal pellets
Gwenaelle Gremion's avatar
Gwenaelle Gremion committed
1488
      w_fp_m=(g*((2*size_fp)**2)*Rfp)/(18*dynvis) ! [m/s]
1489
1490
1491
1492
1493
1494
1495
1496
1497


!!--------Modified Stoke's law (For cylindrical particles) from McDonnell et al 2010
!   Utilisation of spherical diameter as D of particles, and respective values of L,I,S to represent particle shapes
!Even if the values of L,I,S of phy and dph represents a spehrical particles, it will be used in the same way following this special stoke's law for cylindrical
!    [w=((0.079*(r_X-r_F)*L^2*g)/dynvis)*(L/D)^-1.664].

else if (Phys_w .eq. 2.0) then !#5

Gwenaelle Gremion's avatar
Gwenaelle Gremion committed
1498
1499
      !--Phytoplankton 
      w_p_m=((0.079*Rp*((2*size_phy)**2)*g)/dynvis)*((2*size_phy)/(diam_phy))**(-1.664) ! [m/s]
1500

Gwenaelle Gremion's avatar
Gwenaelle Gremion committed
1501
1502
      !--Dead phytoplankton
      w_dph_m=((0.079*Rdph*((2*size_dph)**2)*g)/dynvis)*((2*size_dph)/(diam_dph))**(-1.664)! [m/s]
1503

Gwenaelle Gremion's avatar
Gwenaelle Gremion committed
1504
1505
1506
1507
1508
      !--Dead Zooplankton
      w_dzo_m = ((0.079*Rdzo*((4*size_dzo)**2)*g)/dynvis)*((4*size_dzo)/(diam_dzo))**(-1.664)  ! [m/s]
  
      !--Fecal pellets
      w_fp_m = ((0.079*Rfp*((4*size_fp)**2)*g)/dynvis)*((4*size_fp)/(diam_fp))**(-1.664)  
1509
1510
1511
1512
1513

!!--------Power law evolution from Alldredge and Gotschalk, 1988 based on Kajihara 1971 and from Komar 1981
!                    Calculation of the settling velocity (m/j)/secs_pr_day = m/s - utilisation here of the diameters 
else if (Phys_w .eq. 3.0) then !#5

Gwenaelle Gremion's avatar
Gwenaelle Gremion committed
1514
      !--Phytoplankton ! [m/s]
1515
1516
1517
1518
1519
1520
1521
1522
      if (diam_phy .gt. 0.001) then !#5.1
       w_p_m =-(50*((diam_phy)**0.26))/secs_pr_day 
      else if (diam_phy .lt. 0.001 .and. diam_phy .gt. 0.0) then !#5.1
       w_p_m =-(160*((diam_phy)**0.57))/secs_pr_day  
      else !#5.1
       w_p_m =0.0
      endif!#5.1

Gwenaelle Gremion's avatar
Gwenaelle Gremion committed
1523
   !--Dead Phytoplankton ! [m/s]
1524
1525
1526
1527
1528
1529
1530
1531
      if (diam_dph .gt. 0.001) then !#5.2
      w_dph_m =-(50*((diam_dph)**0.26))/secs_pr_day 
      else if (diam_dph.lt. 0.001 .and. diam_dph.gt. 0.0) then !#5.2
      w_dph_m =-(160*((diam_dph)**0.57))/secs_pr_day  
      else !#5.2
      w_dph_m =0.0
      endif !#5.2

Gwenaelle Gremion's avatar
Gwenaelle Gremion committed
1532
   !--Dead Zooplankton ! [m/s]
1533
1534
1535
1536
1537
1538
1539
1540
      if (diam_dzo .gt. 0.001) then !#5.3
      w_dzo_m =-(50*((diam_dzo)**0.26))/secs_pr_day 
      else if (diam_dzo.lt. 0.001 .and. diam_dzo.gt. 0.0) then !#5.3
      w_dzo_m =-(160*((diam_dzo)**0.57))/secs_pr_day  
      else!#5.3
      w_dzo_m =0.0
      endif !#5.3

Gwenaelle Gremion's avatar
Gwenaelle Gremion committed
1541
   !-- Fecal pellets ! [m/s]
1542
1543
1544
1545
1546
1547
1548
1549
1550
1551
1552
1553
      if (diam_fp .gt. 0.001) then !#5.4
      w_fp_m =-(50*((diam_fp)**0.26))/secs_pr_day 
      else if (diam_fp.lt. 0.001 .and. diam_fp.gt. 0.0) then !#5.4
      w_fp_m =-(160*((diam_fp)**0.57))/secs_pr_day  
      else!#5.4
      w_fp_m =0.0
      endif!#5.4

!!--------Stokes law for high reynolds number + Corey Shape Factor [Komar 1978]
else if (Phys_w .eq. 4.0) then !#5

      !--Phytoplankton 
Gwenaelle Gremion's avatar
Gwenaelle Gremion committed
1554
1555
      ! w_p_m=(1/(18*(kinvis*densFlu)))*(1/CSF_phy_calc)*(Rp)*g*(diam_phy)**2 ![m/s]
       w_p_m=(1/(18*(kinvis*densFlu)))*(1/CSF_phy_calc)*(Rp/densFlu)*g*(diam_phy)**2 ![m/s]
Gwenaelle Gremion's avatar
Gwenaelle Gremion committed
1556
  
1557
      !--Dead phytoplankton
Gwenaelle Gremion's avatar
Gwenaelle Gremion committed
1558
1559
     !  w_dph_m=(1/(18*(kinvis*densFlu)))*(1/CSF_dph_calc)*(Rdph)*g*(diam_dph)**2  ![m/s]
       w_dph_m=(1/(18*(kinvis*densFlu)))*(1/CSF_dph_calc)*(Rdph/densFlu)*g*(diam_dph)**2 ![m/s]
1560
1561

      !--Dead Zooplankton
Gwenaelle Gremion's avatar
Gwenaelle Gremion committed
1562
1563
   !   w_dzo_m=(1/(18*(kinvis*densFlu)))*(1/CSF_dzo_calc)*(Rdzo)*g*(diam_dzo)**2 ![m/s]
      w_dzo_m=(1/(18*(kinvis*densFlu)))*(1/CSF_dzo_calc)*(Rdzo/densFlu)*g*(diam_dzo)**2 ![m/s]
Gwenaelle Gremion's avatar
Gwenaelle Gremion committed
1564

1565
      !--Fecal pellets
Gwenaelle Gremion's avatar
Gwenaelle Gremion committed
1566
1567
   !   w_fp_m=(1/(18*(kinvis*densFlu)))*(1/CSF_fp_calc)*(Rfp)*g*(diam_fp)**2 ![m/s]
      w_fp_m=(1/(18*(kinvis*densFlu)))*(1/CSF_fp_calc)*(Rfp/densFlu)*g*(diam_fp)**2 ![m/s]
1568
1569
1570
1571
1572
1573

!!--------Depending on CSF + Re
else if (Phys_w .eq. 5.0) then !#5

  !--Phytoplankton
  if (CSF_phy .eq. 1.0) then ! Considered as a sphere
Gwenaelle Gremion's avatar
Gwenaelle Gremion committed
1574
1575
1576
        !   w_p_m = -(((Rp/densFlu)*g*(diam_phy**2))/(18*kinvis)) ![m/s]
         w_p_m=(1/(18*(kinvis*densFlu)))*(1/CSF_phy_calc)*(Rp/densFlu)*g*(diam_phy)**2 ![m/s]
    !!        write(*,*) 'CSF_phy = 1 --> Sphere, then Stokes Law ',w_p_m
1577
1578
1579
1580
1581
1582
1583
1584
1585
        ! Comparaison with the Stokes range ! VanRijn 1993
            Re_phy = abs(w_p_m*(diam_phy/kinvis))            
         if (Re_phy .gt. 1.0) then
            w_p_m = -(diam_phy**0.5/secs_pr_day