Commit 3d7dc6b7 authored by Gwenaelle Gremion's avatar Gwenaelle Gremion
Browse files

Avant mise a jour Settl. Vel continue

parent 4fa08a35
......@@ -1356,12 +1356,10 @@ endif !#4
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.)))
write(*,*)'diam_dph',diam_dph
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.)))
write(*,*)'diam_fp',diam_fp
!----------Calcul of densities
!Fluid density
......@@ -1478,7 +1476,7 @@ else if (Phys_w .eq. 4.0) then !#5
!--Dead phytoplankton
!-- Calculation of the CSF factor
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 = 2.18-(2.09*CSF_dph)
else if (CSF_dph .ge. 0.4 .and. CSF_dph .lt. 0.8) then !#5.7
......@@ -1513,7 +1511,7 @@ write(*,*)'CSF_dph',CSF_dph
!--Fecal pellets
!-- Calculation of the CSF factor
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 = 2.18-(2.09*CSF_fp)
else if (CSF_fp .ge. 0.4 .and. CSF_fp .lt. 0.8) then !#5.9
......@@ -1559,9 +1557,7 @@ if (Re_fp .gt. 1.0)then !#9
w_fp_m = -diam_fp**0.5/secs_pr_day ![m/s]
write(*,*) 'Reynolds number for fecal pellets > 1'
endif !#9
write(*,*)'w_dph',w_dph
write(*,*)'w_fp',w_fp
!Everyone will settle at the end of the loop (see below)
!! ---Condition of the environment to allow collision and aggregation
......@@ -2240,6 +2236,13 @@ endif !#1
write(*,*)'Warning :msn incomming flux .ne. to the sum of particles flux at depth ',ci
endif
! ---Flux out from D1, d2, d3, or P can't be superior to the concentration of D1,D2,D3 and P
! if(Flux_P(ci) .gt. cc(p,ci) .Or. Flux_D1(ci) .gt. cc(d1,ci) &
! .or. Flux_D2(ci) .gt. cc(d2,ci) .or. Flux_D3(ci) .gt. cc(d3,ci)) then
! 99 FATAL 'Coagulation rate overpass [Variables]'
! stop 'do_bio'
! endif
! ---Size modification due to Coagulation
if (Coag_coef .eq. 1.0) then ! #1 Size modified only if coagulation happenned
......@@ -2314,11 +2317,11 @@ endif !#1
endif
!--» If there is no marine snow, there is no reason to set up a size, then the size will equal 0.
if (cc(d4,ci).lt. cons_min)then
cc(size_msn,ci) = _ZERO_
cc(size_intrm,ci) = _ZERO_
! if (cc(d4,ci).lt.cons_min)then
! cc(size_msn,ci) = _ZERO_
! cc(size_intrm,ci) = _ZERO_
! write(*,*)'Size but no msn at :',ci
endif
! endif
!-----------------------------------------------------------------------------------
! FRAGMENTATION
......@@ -2486,13 +2489,13 @@ endif
!--» cc(size_msn,ci) --> is the diameter after coag and frag
Rmsn = (rho_msn-densFlu)
if (cc(d4,ci) .le. cons_min) then
w_msn_m = 0.0
!if (cc(d4,ci) .le. cons_min) then
! w_msn_m = 0.0
! write(*,*)'[msn] lt Cons_min, settling velocity set to :',w_msn_m
!--» We set here the concentration limit in the water column from free settling of marine snow vs flocculation settling (Metha 1989)
elseif (cc(d4,ci) .lt. cons_max) then !#1
!elseif (cc(d4,ci) .lt. cons_max) then !#1
if (cc(d4,ci) .lt. cons_max) then !#1
!! --- Calculated by value from the .nml
if (w_msnow .eq. 0.0 ) then !#2
w_msn_m = w_msn ![m/s]
......@@ -2535,9 +2538,29 @@ elseif (cc(d4,ci) .lt. cons_max) then !#1
w_msn_m=(1/(18*(kinvis*densFlu)))*(1/CSF_msn)*(Rmsn)*g*(cc(size_msn,ci))**2 ![m/s]
elseif (w_msnow .eq. 5.0) then !#2!
!The choice of settling velocity will depend on the choice for diam_msn_us by the user. This will apply no matter the size of the particle of msn
if( diam_msn_us .gt. 0.000001 .and. diam_msn_us .le. 0.0001 )then
! If the diameter of Msn (d) --> 1< d <= 100 micrometers (VanRijn 1993)
w_msn_m = -(((2.65-1)*g*(cc(size_msn,ci)**2))/18*kinvis) ![m/s]
write(*,*) 'Diameter of msn : 1< d <= 100 micrometers, settling =' ,w_msn_m
elseif( diam_msn_us .gt. 0.0001 .and. diam_msn_us .lt. 0.001 ) then
! If the diameter of Msn (d) --> 100< d <= 1000 micrometers (VanRijn 1993)
w_msn_m = -(((10*kinvis)/cc(size_msn,ci))*&
((1+((0.01*(2.65-1)*g*(cc(size_msn,ci))**3.0) /kinvis**2.0)**0.5) -1)) ![m/s]
write(*,*) 'Diameter of msn : 100< d <= 1000 micrometers, settling =' ,w_msn_m
elseif( diam_msn_us .gt. 0.001) then
! If the diameter of Msn (d) --> d > 1000 micrometers (VanRijn 1993)
w_msn_m = -(1.1*((2.65-1)*g*cc(size_msn,ci))**0.5)/secs_pr_day ![m/s]
write(*,*) 'Diameter of msn : d > 1000 micrometers, settling =' ,w_msn_m
else
w_msn_m =- (g*((cc(size_msn,ci))**2)*Rmsn)/(18*dynvis) ![m/s]
write(*,*) 'Diam.of msn < than 1 micrometers: Stokes Law is used ! ',w_msn_m
endif
! In theory should be use only it if the diameter of Msn (d) --> 1< d <= 100 micrometers (VanRijn 1993)
w_msn_m = -(((2.65-1)*g*(cc(size_msn,ci)**2))/18*kinvis)
elseif (w_msnow .eq. 6.0) then !#2!
......@@ -2590,6 +2613,7 @@ if (Re_msn .lt. 1.0)then !#3
! write(*,*) 'Reynolds number for marine snow allow to respect Stockes range < 1'
else !#3
write(*,*) 'Reynolds number for marine snow > 1'
w_msn_m = -(cc(size_msn,ci)**0.5/secs_pr_day) ![m/s]
if (w_msnow .eq. 1.0 .or. w_msnow .eq. 2.0 .or. w_msnow .eq. 3.0 )then
write(*,*) 'Reynolds number for marine snow > 1, another settling velocity scheme should be chose for msn'
endif
......
......@@ -249,7 +249,7 @@
do k=1,N-1
c=abs(ww(k))*dt/(0.5*(h(k)+h(k+1)))
if (c.gt.cmax) cmax=c
enddo
enddo
it=min(itmax,int(cmax)+1)
......@@ -418,7 +418,7 @@ endif
end if
case (zeroDivergence)
cu(0) = cu(1)
case (zeroDivergence3)
case (zeroDivergence3) ! For Msn_size in Bio_model=10
cu(0) = ww(1)*Y(1)
case default
FATAL 'unkown lower boundary condition type in adv_center()'
......@@ -441,41 +441,6 @@ endif
Y(k)=Y(k)-1./float(it)*dt*((cu(k)-cu(k-1))/h(k))
enddo
! For the other variables of bio_polynow, or other biological models
elseif (mode.eq.3) then
do k=1,N
! compute courant number
c=-ww(k)/float(it)*dt/(0.5*(h(k)+h(k+1)))
if (k .lt. N-1) then
Yu=Y(k+2) ! upstream value
else
Yu=Y(k+1)
endif
Yc=Y(k+1)
Yd=Y(k )
! compute slope ratio
if (abs(Yc-Yd) .gt. 1e-10) then
r=(Yu-Yc)/(Yc-Yd)
else
r=(Yu-Yc)*1.e10
endif
! compute flux limiter
x = one6th*(1.-2.0*c)
Phi = (0.5+x)+(0.5-x)*r
! Superbee limiter
limit=max(_ZERO_, min(_ONE_, 2.0*r), min(r,2.*_ONE_) )
! compute limited flux
cu(k) =ww(k)*(Yc+0.5*limit*(1-c)*(Yd-Yc))
cu(k-1)=ww(k)*(Yd+0.5*limit*(1-c)*(Yd-Y(k-1)))
! advection step
Y(k)=Y(k)-1./float(it)*dt*((cu(k)-cu(k-1))/h(k))
enddo
else ! conservative
do k=1,N
Y(k)=Y(k)-1./float(it)*dt*((cu(k)-cu(k-1))/h(k))
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment