Commit dae7253b authored by Gwenaelle Gremion's avatar Gwenaelle Gremion
Browse files

Nettoyage de adv_center

parents 20015aff 545c01fb
......@@ -1356,11 +1356,12 @@ 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
......@@ -1477,7 +1478,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
......@@ -1512,7 +1513,7 @@ else if (Phys_w .eq. 4.0) then !#5
!--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
......@@ -1558,7 +1559,9 @@ 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
......
......@@ -439,6 +439,41 @@ 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