Commit 8cdf10f3 authored by Dany Dumont's avatar Dany Dumont
Browse files

Premier depot

parents
This source diff could not be displayed because it is too large. You can view the blob instead.
1 24.9997616 24.9997616 1.25000000 13.7087755 1.00000000 3.17407346 1.24999997E-03 0.00000000
1 24.9997616 24.9997616 3.75000000 90.7216568 1.00000000 3.13634229 3.74999992E-03 2.50000000
1 24.9997616 24.9997616 6.25000000 288.376648 1.00000000 3.09815192 6.24999963E-03 5.00000000
1 24.9997616 24.9997616 8.75000000 644.005127 1.00000000 3.05948496 8.75000004E-03 7.50000000
1 0.00000000 0.00000000 15.0000000 2863.12500 1.00000000 2.96060848 1.12500004E-02 10.0000000
0 1 24.9997616 1.25000000 13.7087755 1.00000000 3.17407346
0 1 24.9997616 3.75000000 90.7216568 1.00000000 3.13634229
0 1 24.9997616 6.25000000 288.376648 1.00000000 3.09815192
0 1 24.9997616 8.75000000 644.005127 1.00000000 3.05948496
0 1 0.00000000 15.0000000 2863.12500 1.00000000 2.96060848
1 1 23.4998455 1.25000000 13.7087755 1.00000000 3.17407346
1 1 23.4998169 3.75000000 90.7216568 1.00000000 3.13634229
1 1 23.4997883 6.25000000 288.376648 1.00000000 3.09815192
1 1 17.4998188 8.75000000 644.005127 1.00000000 3.05948496
1 1 11.9997730 15.0000000 2863.12500 1.00000000 2.96060848
!'********************************************************'
! COAGULATION - FRAGMENTATION - 0D
! Co-build by Gremion Gwenaëlle and Louis-Philippe Nadeau
! Institut des Sciences de la Mer - Université du Québec à Rimouski
!
! Contact: gwenaelle.gremion@gmail.com , Louis-Philippe_Nadeau@uqar.ca
! Last update: 4 DEC - 2020
!'********************************************************'
!%%%% SET UP %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!-------------------------------------------------------------
! DECLARATIONS of Model Parameters
!-------------------------------------------------------------
MODULE data_initial
IMPLICIT NONE
include 'parameters.f90' ! List of parameters to declare for the model
END MODULE data_initial
!%%%% PROGRAM %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
program coagfrag
USE data_initial
IMPLICIT NONE
!-------------------------------------------------------------
! DECLARATIONS OF VARIABLES
!-------------------------------------------------------------
include 'variables.f90' ! List of Variables for the model
!-------------------------------------------------------------
! INITIALISATION OF VARIABLES
!------------------------------------------------------------
print*,' INITIALISATION '
print*,'********************************************************'
include 'initialisation.f90' ! Initialisation of the Simulation + Initial profile of the Highest Resolution (†)
include 'outputs_profile_init.f90' ! Creation of the output files containing the initial profile
!----------------------------------------------------------
!! REACTIONS [ADVECTION - DIFFUSION]
!----------------------------------------------------------
print*,' LOOP '
print*,'********************************************************'
include 'outputs_profile_results.f90' ! Creation of the output files containing the initial profile + results
!-------------------------------------------------------------------------------------
do t = 1,Nt ! Timestepping
do k= 1,nlev ! Vertical loop
include 'reactions.f90' ! Reactions of Coagulation and Fragmentation (related to functions.f90)
! include 'adv_flux.f90' ! Advective fluxes if needed (file not provided)
end do !--> nlev
! include 'diff_adv.f90' ! Diffusion & Advection if need (file not provided)
!-- When only reactions apply
C1_sp_mid_lowRes(:,:) = C2_sp_mid_lowRes(:,:)! [mmolN m-3(w)]
!-- When reactions - advection - diffusion apply
! C1_sp_mid_lowRes(:,:) = C2_sp_mid_lowRes(:,:)+ dt * (-Adv_sp(:,:) + DIFF_sp(:,:)
include 'outputs_profile_results.f90' ! Creation of output files containing the initial profile + results at the end of each timestep
end do !--> End of Timestepping
!-------------------------------------------------------------------------------------
!----------------------------------------------------------
!! SUM UP ON SCREEN
!----------------------------------------------------------
! print*,'********************************************************'
print*,' SUM UP '
print*,'********************************************************'
!print*, ' Total concentration Spectra Middle bin-LowRes '
print*,sum(C1_sp_mid_lowRes(:,1))
!print*, ' Concentration Spectra Middle bin -LowRes'
print*,C1_sp_mid_lowRes(:,1)
end program coagfrag
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!-------------------------------------------------------------------------------------
!! SUB ROUTINES
!-------------------------------------------------------------------------------------
include 'functions.f90'
!#####################################################################################
!! SUB ROUTINES
!#####################################################################################
!-------------------------------------------------------------
! Diameter and Volume
!-------------------------------------------------------------
!Call Line : Series_DandV(nc_sp_lowRes,D_sp_lowRes,V_sp_lowRes,D_sp_mid_lowRes,V_sp_mid_lowRes)
subroutine Series_DandV(loop,D_lim,V_lim,D_mid,V_mid)
USE data_initial
IMPLICIT NONE
!--IN
Integer loop
!--OUT
real D_lim(1:loop+1)
real V_lim(1:loop+1)
real D_mid(1:loop)
real V_mid(1:loop)
!--INSIDE
Integer i,k
real Step
real x(1:loop)
! Diameter at the Bounds
!-------------------------------------------------------------------------------------
!--Logarithmic Scale - E2 & E2_res_dep_func
! Desire that Initial and Final size included in the simulation
!Step = (Final - Initial)/(loop-1)
!x = (/(Initial + (i-1)*Step, i=1, loop)/)
! D_lim = 10**x ! [m]
!--Linear Scale - E1
do i = 1,loop
D_lim(i) = 10.*(i-1)/(loop-1) ! [m]
enddo
!--For N*, upper limit equal 2 times the max size consider (loop)
D_lim(loop+1) = 2*D_lim(loop)
! Volume at the Bounds
!-------------------------------------------------------------------------------------
V_lim = 2.8*D_lim**2.49 ! [m3]
! Diameter AND Volume at the Middle of the size bins
!-------------------------------------------------------------------------------------
k=1
do i = 2,loop+1
D_mid(k) = (D_lim(i)+D_lim(i-1))/2.
V_mid(k) = (V_lim(i)+V_lim(i-1))/2.
k=k+1
enddo
end subroutine Series_DandV
!-------------------------------------------------------------
! MATRICE VOLUME REACTIONS
!-------------------------------------------------------------
!Call Line : Volume_reactions(nc_sp_lowRes-1,V_sp_mid_lowRes,Vr_sp_mid_lowRes,Vr_ind_sp_mid_lowRes)
subroutine Volume_reactions(loop,V,Vr,Vr_ind)
USE data_initial
IMPLICIT NONE
!--IN
Integer loop
!--OUT
real Vr(1:loop+1,1:loop+1)
real V(1:loop+1)
INTEGER Vr_ind(1:loop,1:loop)
!--INSIDE
Integer i,j,k
do j = 1,loop+1
do i =1,loop+1
Vr(i,j) = V(i) + V(j)
if (Vr(i,j) > V(loop) ) then
Vr(i,j) = V(loop+1) ! [m3]
end if
end do
end do
do j = 1,loop
do i =1,loop
k =1
do while (Vr(i,j) > V(k))
k=k+1
end do
Vr_ind(i,j) = k
end do
end do
end subroutine Volume_reactions
!-------------------------------------------------------------
! ELEMENT PARTICLES CONTENT & PARTICLES DENSITY
!-------------------------------------------------------------
!Call Line : Contentrho(nc_sp-1,V_sp_mid,Nit_sp_mid,rho_sp_mid)
subroutine Contentrho(loop,V,Nit,rho)
USE data_initial
IMPLICIT NONE
!--IN
Integer loop
!Note : We do not consider the element particles content dependant of V for size bin 3/2N .
real V(1:loop)
!--OUT
real Nit(1:loop)
real rho(1:loop+1)
!--INSIDE
Integer i
! Element particles content the Middle of the size bins[From Alldredge 1998]
!-------------------------------------------------------------------------------------
do i = 1,loop
!--Linear Scale - E1
Nit(i) = 1. ![mmol element]
!--Logarithmic Scale - E2 & E2_res_dep_func
!Nit(i) = 1.09*V(i)**0.52 ! [mmol element] Parameters are set regarding the element of interest (e.g., Nitrogen, Carbon)
enddo
! Density of particle for the Middle of the size bins
!-------------------------------------------------------------------------------------
do i = 1,loop+1
rho(i) = 1050. ! Density [kg.m-3(w)]
enddo
end subroutine Contentrho
!-------------------------------------------------------------
! PARTICLES SETTLING VELOCITY
!-------------------------------------------------------------
!Call Line : Speed(nc_sp_lowRes, D_sp_mid_lowRes, rho_sp_mid_lowRes, w_sp_mid_lowRes )
subroutine Speed( loop, D, rho, w )
USE data_initial
IMPLICIT NONE
!--In
Integer loop
real rho(1:loop)
!--OUT
real w(1:loop)
real D(1:loop)
!--INSIDE
Integer i
do i = 1,loop
if (D(i) .le. 0.0001) then
w(i) = ((g*D(i)**2)/(18*kinvis*1))*((rho(i)-rho_w)/(rho_w))
elseif (D(i) .gt. 0.0001 .and. D(i) .lt. 0.001) then
w(i) = ((10.0*kinvis)/D(i))*(((1+((0.01*g*D(i)**3.0) / &
(kinvis**2.0))*((rho(i)-rho_w)/rho_w))**0.5)-1) ![m/s]
! write(*,*)'D(i) .gt. 0.0001 .and. D(i) .lt. 0.001 for i =' ,i
elseif (D(i) .ge. 0.001 )then
w(i)= ((((((10.0*kinvis)/0.001)* &
(((1+((0.01*(rho(i)-rho_w/rho_w)*g*(0.001**3.0))/(kinvis**2.0)))**0.5)-1)) &
/sqrt(rho(i)-rho_w*g*0.001/rho_w)))*sqrt(rho(i)-rho_w*g*D(i)/rho_w))
! write(*,*)'D(i) .ge. 0.001 for i =' ,i
endif
! Verification that Settling velocities are positive (to go downward)
! (So particles' density need to be higher than the one of seawater
if (w(i) .lt. 0.0 ) then
! w(i) = abs(w(i))
print*, 'Settling velocity going upward'
endif
! Check reagrding the CFL limit
if(w(i) .gt. 1.0 / dt ) then
! w(i) = 1. / dt
! write(*,*)'CFL overpass for w(i) with i = ' ,i, w(i)
! D(i) = 0.035
! write(*,*)'CFL overpass for D(i) with i = ' ,i, D(i)
endif
end do
end subroutine Speed
!-------------------------------------------------------------
! COLLISION KERNELS FOR COAGULATION REACTION
!-------------------------------------------------------------
!Call Line : ColKernels(nc_sp_lowRes-1,w_sp_mid_lowRes,D_sp_mid_lowRes,kernel_coef_sp_mid_lowRes,&
! betaBr_sp_mid_lowRes,betaSh_sp_mid_lowRes,betaDs_sp_mid_lowRes,q_sp_mid_lowRes,abs_w_sp_mid_lowRes)
subroutine ColKernels(loop,w,D,kernel_coef,betaBr,betaSh,betaDs,q,abs_w)
USE data_initial
IMPLICIT NONE
!--IN
Integer loop,k
real w(1:loop),D(1:loop)
!--OUT
REAL kernel_coef(1:loop, 1:loop)!Collision kernel (sum of Br,Sh,Ds)[m3.s-1]
REAL betaBr(1:loop, 1:loop) !Collision kernel for brownian motion[m3.s-1]
REAL betaSh(1:loop, 1:loop) !collision kernel for shear[m3.s-1]
REAL betaDs(1:loop, 1:loop) !collision kernel for differential settling [m3.s-1]
REAL q(1:loop, 1:loop) !Parameter of collision kernel [m]
REAL abs_w(1:loop, 1:loop) !Parameter of collision kernel [m.s-1]
!--INSIDE
Integer i,j
do j = 1,loop
do i = j,loop
betaBr(i,j)=(physic_Br)*(4.*(0.5*(D(i)+D(j))**2)/(D(i)*D(j)))
q(i,j)= (min(D(i),D(j)))/(max(D(i),D(j)))
betaSh(i,j)=9.8*((q(i,j)**2)/(1+2* q(i,j)**2))*(sqrt(bio_eps/kinvis))*(0.5*(D(i)+D(j)))**3
abs_w(i,j)= ABS(w(i)-w(j))
! betaDs(i,j)=(1/2)*2.*asin(1.)*(0.5*min(D(i),D(j))**2)*abs_w(i,j)
! betaDs(i,j)=(pi/2)*(min(D(i),D(j))**2)*abs_w(i,j)
betaDs(i,j)= 0.
kernel_coef(i,j)= (betaBr(i,j)+betaSh(i,j)+betaDs(i,j) ) ![m3.s-1]
end do
end do
end subroutine ColKernels
!-------------------------------------------------------------
! REACTIONS
!-------------------------------------------------------------
!Call Line : Koalaful(nc_sp_lowRes-1, &kernel_coef_sp_mid_lowRes(:,:),&
! C2_sp_mid_lowRes(:,k), D_sp_mid_lowRes(:),Vr_ind_sp_mid_lowRes(:,:),&
! Coag1, Frag1)
subroutine Koalaful(loop,kernel_coef,C2,D,Vr_ind,coag,frag)
USE data_initial
IMPLICIT NONE
!--IN
Integer loop,n
REAL kernel_coef(1:loop, 1:loop)
REAL C2(1:loop)
Real D(1:loop)
integer Vr_ind(1:loop,1:loop)
!--OUT
!--INSIDE
Integer i,j
REAL coag(1:loop, 1:loop) ,frag(1:loop, 1:loop)
do j = 1,loop
do i =j,loop
! ----COAGULATION = Number of Collision = Number of particles formed [Jackson 2001]
!Constant Value over the range
! ----E1
coag(i,j) = ((Tx_coag / dt)* C2(i) * C2(j) )*(loop/(loop+1.))
! ----E2
! coag(i,j) = ((Tx_coag ) * C2(i) * C2(j) )*(loop/(loop+1.))
! ----E2_res_dep_func
! coag(i,j) = ((Tx_coag ) * C2(i) * C2(j) )*(loop/(loop+1.))*&
! (fac1*(1-exp(fac2*loop))+(1-fac1))
!Depending on Kernel Coefficient
! coag(i,j) = (kernel_coef(i,j) * C2(i) * C2(j) )*(loop/(loop+1.))
!** UNITS **
! coag(i,j) = alpha_ij* kernel_coef(i,j)* C2(i) * C2(j)
! [(mmolN)² m-3 s-1] = [ - * m3.s-1 * mmolN.m-3 * mmolN.m-3]
! ----FRAGMENTATION
if (Vr_ind(i,j) .gt. loop) then
! Constant Value over the range
! ----E1
frag(i,j) = ((Tx_frag/dt) *( C2(Vr_ind(i,j)) /( loop*((loop+1.)/2.)) ) ) * (1./(loop+1.))! E1
! ----E2
! frag(i,j) = ((Tx_frag) *( C2(Vr_ind(i,j)) /( loop*((loop+1.)/2.)) ) ) * (1./(loop+1.))
! ----E2_res_dep_func
! frag(i,j) = ((Tx_frag) *( C2(Vr_ind(i,j)) /( loop*((loop+1.)/2.)) ) ) * (1./(loop+1.))*&
! (fac1*(1-exp(fac2*loop))+(1-fac1))
else
!Constant Value over the range
! ----E1
frag(i,j) = (Tx_frag/dt) * C2(Vr_ind(i,j)) * (1./(loop+1.))
! ----E2
! frag(i,j) = (Tx_frag) * C2(Vr_ind(i,j)) * (1./(loop+1.))
! ----E2_res_dep_func
!frag(i,j) = (Tx_frag) * C2(Vr_ind(i,j)) * (1./(loop+1.))*&
! (fac1*(1-exp(fac2*loop))+(1-fac1))
endif
!** UNITS **
! frag(i,j) = Tx_Frag*C2(Vr_ind(i,j))
! [mmolN m-3 s-1] = [s-1 * mmolN.m-3 ]
enddo
enddo
end subroutine Koalaful
!-------------------------------------------------------------
! Concentrations
!-------------------------------------------------------------
!Call Line : Concentrations(nc_sp_lowRes-1,Coag1,Frag1,Vr_ind_sp_mid_lowRes(:,:),&
! Nit_sp_mid_lowRes,C2_sp_mid_lowRes(:,k))
subroutine Concentrations(loop,coag,frag,Vr_ind,Nit,C2)
USE data_initial
IMPLICIT NONE
!--IN
integer loop
integer Vr_ind(1:loop,1:loop)
REAL Nit(1:loop)
!--OUT
REAL C2(1:loop+1)
!--INSIDE
Integer i,j,k
REAL tmpi, tmpj , tmpk
REAL coag(1:loop, 1:loop) ,frag(1:loop, 1:loop)
select case (run_type)
!-------------------------------------------------------------------------------------
case (1) ! Coagulation ONLY
frag(:,:)= 0.
case (2) ! Fragmentation ONLY
coag(:,:)= 0.
case (3) ! BOTH (Coagulation & fragmentation)
!-------------------------------------------------------------------------------------
end select
!--- Attribution to concentration to the size classes---------------------------
! To be equal we will make half of the matrice reacting from littles ones to big ones (Even)
! and the other half from big ones to little one (Odd) to assure to not favorise any of the reactions
!-------------------------------------------------------------------------------------
!--- i=j !!
do j = 1, loop ,1
i = j
k = Vr_ind(i,j)
! -- Calculation of what have reacted for each class --
tmpi = dt* ( coag(i,j)/Nit(j) - frag(i,j) )
! [mmolN m-3] = [ s * [(mmolN)² m-3 s-1 / mmolN - mmolN m-3 s-1]
if( C2(i)-tmpi <= tolerance ) then
C2(k) = C2(k) + C2(i)
C2(i) = 0.
elseif( C2(k)+tmpi <= tolerance ) then
C2(i) = C2(i) + C2(k)
C2(k) = 0.
else
C2(i) = C2(i) - tmpi
C2(k) = C2(k) + tmpi
end if
end do
!--- i~=j !!
do j = 1, loop ,1
do i = j+1,loop,1
k = Vr_ind(i,j)
! -- Calculation of what have reacted for each class --
tmpi = dt* ( 0.5*coag(i,j)/Nit(j) - (Nit(i)/(Nit(i)+Nit(j)))*frag(i,j) )
tmpj = dt* ( 0.5*coag(i,j)/Nit(i) - (Nit(j)/(Nit(i)+Nit(j)))*frag(i,j) )
tmpk = tmpi+tmpj
! [mmolN m-3] = [ s * [(mmolN)² m-3 s-1 / mmolN - mmolN m-3 s-1]
if( C2(i)-tmpi <= tolerance .AND. C2(j)-tmpj <= tolerance ) then
C2(k) = C2(k) + C2(i) + C2(j)
C2(i) = 0.
C2(j) = 0.
!print*,'4','sum',sum(C2(:)),C2(:)
elseif( C2(i)-tmpi <= tolerance ) then
if( i .ne. k) then
C2(j) = C2(j) - tmpj
C2(k) = C2(k) + tmpj+C2(i)
C2(i) = 0.
else
C2(j) = C2(j) - tmpj
C2(k) = C2(k) + tmpj
endif
!print*,'5','sum',sum(C2(:)),C2(:)
elseif( C2(j)-tmpj <= tolerance ) then
if( j .ne. k) then
C2(i) = C2(i) - tmpi
C2(k) = C2(k) + tmpi+C2(j)
C2(j) = 0.
else
C2(i) = C2(i) - tmpi
C2(k) = C2(k) + tmpi
endif
!print*,'6','sum',sum(C2(:)),C2(:)
elseif( C2(k)+tmpk <= tolerance ) then
C2(i) = C2(i) + (Nit(i)/(Nit(i)+Nit(j)))*C2(k)
C2(j) = C2(j) + (Nit(j)/(Nit(i)+Nit(j)))*C2(k)
C2(k) = 0.
!print*,'7' ,'sum',sum(C2(:)),C2(:)
else
C2(i) = C2(i) - tmpi
C2(j) = C2(j) - tmpj
C2(k) = C2(k) + tmpk
!print*,'8','sum',sum(C2(:)),C2(:)
end if
end do
end do
end subroutine Concentrations
!**************************************************************************************
!! INITIALISATION OF VARIABLES WITH 0
!**************************************************************************************
! %%%%%% HIGHEST resolution (†,nc_sp) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
! Diameter to the Bounds of our size bins
!-------------------------------------------------------------------------------------
D_sp(:) = 0.0
V_sp(:) = 0.0
! Diameter to the middle of our size bins
!-------------------------------------------------------------------------------------
D_sp_mid(:) = 0.0
V_sp_mid(:) = 0.
Nit_sp_mid(:) = 0.
rho_sp_mid(:) = 0.
C1_sp_mid(:,:) = 0.
C1_in_Interval(:) = 0.
C_Count_Interval(:) = 0.
! %%%%%% LOW resolutions (HR & LR) (nc_sp_lowRes) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
! Variables that require the definition of the 3/2N size bin
!-------------------------------------------------------------------------------------
D_sp_lowRes(:) = 0.
D_sp_mid_lowRes(:) = 0.
V_sp_lowRes(:) = 0.
V_sp_mid_lowRes(:) = 0.
Vr_sp_mid_lowRes(:, :) = 0.
Vr_ind_sp_mid_lowRes(:, :) = 0.