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

Correction Units_ Wd4 Variable

parent f5cf019f
......@@ -55,7 +55,8 @@
!-----LOCAL VARIABLES:---------- from a namelist : bio_polynow.nml----------------------
REALTYPE :: dt_bio
REALTYPE :: splitfac_bio
REALTYPE :: splitfac_bio
REALTYPE :: depth_bio
!INITIAL and Minimum concentration for the variable
REALTYPE :: p_init_value=1.0
......@@ -239,13 +240,15 @@ REALTYPE ::coef2
REALTYPE ::RFV
REALTYPE ::coef3
REALTYPE ::coef4
REALTYPE ::coef5
REALTYPE ::coef5
REALTYPE ::cons_max
integer :: out_unit
! integer, parameter :: p=1,z=2,b=3,d1=4,n=5,a=6,l=7,d2=8,d3=9,d4=10,aug_si_d4=11,taille_msn=12
! GG d1= dph , d2= dzo , d3= fp , d4 = msn, aug_si_d4= size_msn
integer, parameter :: p=1,z=2,b=3,d1=4,n=5,a=6,l=7,d2=8,d3=9,d4=10, &
taille_intrm=11,taille_coag=12, taille_frag=13,taille_msn=14
taille_intrm=11,taille_coag=12, taille_frag=13,taille_msn=14,settl_msn=15
!EOP
......@@ -281,7 +284,7 @@ REALTYPE ::coef5
! !LOCAL VARIABLES:
namelist /bio_polynow_nml/ numc, &
dt_bio,splitfac_bio,p_init_value,p_initial,z_p_gauss_init,sigma_p, &
dt_bio,splitfac_bio,depth_bio,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, &
......@@ -310,7 +313,7 @@ rho_p,rho_dph,rho_dzo,rho_fp,rho_msn,
Floc_coef, &
betaBr,betaSh,betaDs, &
eps_const,eps_n,cons_min, &
coef1,coef2,RFV,coef3,coef4,coef5,&
coef1,coef2,RFV,coef3,coef4,coef5,cons_max,&
write_screen
!EOP
......@@ -565,7 +568,8 @@ Floc_coef, &
!cons_min
!Floc_coef
!coef1,coef2,coef3,coef4,coef5
!RFV
!RFV
!cons_max
900 format (a,f8.5)
901 format (f8.5)
......@@ -809,7 +813,10 @@ end if
cc(taille_frag,i)=_ZERO_
cc(taille_msn,i)=_ZERO_
ws(taille_msn,i)=_ZERO_
ws(taille_msn,i)=_ZERO_
cc(settl_msn,i)= _ZERO_
ws(settl_msn,i) =_ZERO_
end do
......@@ -929,6 +936,12 @@ end if
var_names(14) = 'taille_msn'
var_units(14) = 'm'
var_long(14) = 'Final size of marine snow'
var_names(15) = 'settl_msn'
var_units(15) = 'm/s'
var_long(15) = 'Settling Velocity of marine snow'
return
end subroutine var_info_polynow
......@@ -1149,8 +1162,7 @@ end if
REALTYPE ::prob_break_msn,r_msn
REALTYPE :: diam_msn_max,RFV_msn,min_size_msn
! Sedimentation rate msn
REALTYPE :: w_msn_m,Re_msn,Rmsn,w_min_m
! REALTYPE :: sensi_w
REALTYPE :: w_msn_m,Re_msn,Rmsn,w_min_m
!EOP
!-----------------------------------------------------------------------
......@@ -1545,6 +1557,7 @@ endif !#9
w_min_m = max(w_p_m,w_dph_m,w_dzo_m,w_fp_m)
!Now that we defined a settling velocity depending of the size, let's make it dependent of the depth (Martin et al 1987, Paul Nicot 2013 )
......@@ -2204,43 +2217,46 @@ endif !#1
!-----------------------------------------------------------------------------------
! ---Maximum size at dissagregation (Alldredge 1990)
if(dm_msn .eq. 1.0) then
if(dm_msn .eq. 1.0) then
diam_msn_max = coef1*(eps(ci))**(-coef2) ! [m]
write(*,*)'At depth-----',ci
write(*,*)'Maximum size diameter of Msn from eps :', diam_msn_max
else
diam_msn_max = diam_msn_us ! [m] (which means ~ 10mm)
diam_msn_max = diam_msn_us ! [m]
endif
write(*,*)'At depth-----',ci
write(*,*)'Maximum size diameter - Marine snow :', diam_msn_max
! ---Flux from particules to msn [mmol/m3/d]
! ---Flux from particules to msn [mmol/m3/s]
Flux_P(ci) = (agg_pp)+(agg_pdph)+(agg_pdzo)+(agg_pfp)
Flux_D1(ci) = (agg_pdph)+(agg_dphdph)+(agg_dphdzo)+(agg_dphfp)
Flux_D2(ci) = (agg_pdzo)+(agg_dphdzo)+(agg_dzodzo)+(agg_dzofp)
Flux_D3(ci) = (agg_pfp)+(agg_dphfp)+(agg_dzofp)+(agg_fpfp)
! ---Total Flux received by msn [mmol/m3/d]
flux_msn(ci) = Flux_P(ci)+Flux_D1(ci)+Flux_D2(ci)+Flux_D3(ci)
! write(*,*)'flux_msn(ci) ',flux_msn(ci)
! ---Calcul of the ratios [-]
Ratio_P = Flux_P(ci)/flux_msn(ci)
Ratio_D1 = Flux_D1(ci)/flux_msn(ci)
Ratio_D2 = Flux_D2(ci)/flux_msn(ci)
Ratio_D3 = Flux_D3(ci)/flux_msn(ci)
! ---Maximum Flux of Particles to msn [mmol/m3/s]
max_flux= max(Flux_P(ci),Flux_D1(ci),Flux_D2(ci),Flux_D3(ci))
! ---Safety check - Sum of ratio should equal 1
if(Ratio_P+Ratio_D1+Ratio_D2+Ratio_D3 .gt. 1.0) then
! ---Total Flux received by msn [mmol/m3/s]
flux_msn(ci) = Flux_P(ci)+Flux_D1(ci)+Flux_D2(ci)+Flux_D3(ci)
!-----------------------------------------------------------------------------------
!! ---Safety check
!-----------------------------------------------------------------------------------
! ---Sum of ratio should equal 1
if(Flux_P(ci)+Flux_D1(ci)+Flux_D2(ci)+Flux_D3(ci).gt. 1.0) then
write(*,*)'At depth-----',ci
write(*,*)'Warning: sum of Ratio > 1.0', Ratio_P+Ratio_D1+Ratio_D2+Ratio_D3
write(*,*)'Warning: sum of Ratio > 1.0', Flux_P(ci)+Flux_D1(ci)+Flux_D2(ci)+Flux_D3(ci)
endif
if (flux_msn(ci) .ne. Flux_P(ci)+Flux_D1(ci)+Flux_D2(ci)+Flux_D3(ci))then
if (flux_msn(ci) .ne. Flux_P(ci)+Flux_D1(ci)+Flux_D2(ci)+Flux_D3(ci)) then
write(*,*)'Warning :msn incomming flux .ne. to the sum of particles flux at depth ',ci
endif
! ---Size modification due to Coagulation
if (Coag_coef .eq. 1.0) then ! #1 Size modified only if coagulation happenned
! ---Initialize the size radius augmentation variable for the msn and the size for the coagulation
cc(taille_coag,ci) = 0.0
cc(taille_intrm,ci) = 0.0
......@@ -2251,58 +2267,73 @@ if (Coag_coef .eq. 1.0) then ! #1 Size modified only if coagulation happenned
!Calculation of the cubic racin of X --» SIGN(ABS(X)**(1./3.),X) [FORTRAN, le langage normalisé, Michel Dubesset,Jean Vignes]
!-----------------------------------------------------------------------------------
! ---Depend on what composed it
cc(taille_coag,ci)= Ratio_P*diam_phy+Ratio_D1*diam_dph+Ratio_D2*diam_dzo+Ratio_D3*diam_fp ![m]
cc(taille_coag,ci) = cc(taille_coag,ci)/secs_pr_day !Growth of size per day [m/d]
! ---Depends on what composed it
cc(taille_coag,ci)= Flux_P(ci)*diam_phy+Flux_D1(ci)*diam_dph+Flux_D2(ci)*diam_dzo+Flux_D3(ci)*diam_fp ![m/s]
! ---Definition of the diameter size of marine snow
!-----------------------------------------------------------------------------------
! ---Definition of the diameter size of marine snow
!-----------------------------------------------------------------------------------
!At this stage 'taille_coag' is still the radius and diam_msn the equivalent spherical diameter
!We defined L,I,S for Marine snow L=I=size_x*4 and S= size_X*2
!--» With size_X which represents the radius (cc(taille_coag,ci)) of our marine snow.
!-----------------------------------------------------------------------------------
!-----------------------------------------------------------------------------------
if (cc(taille_coag,ci) .ne. 0.0 .and. cc(taille_coag,ci).gt. 0.0) then !#2
diam_msn = sign(abs(((cc(taille_coag,ci)*4.)*(cc(taille_coag,ci)*4.)*(cc(taille_coag,ci)*2.)))**(1./3.),&
((cc(taille_coag,ci)*4.)*(cc(taille_coag,ci)*4.)*(cc(taille_coag,ci)*2.))) ![m/d]
((cc(taille_coag,ci)*4.)*(cc(taille_coag,ci)*4.)*(cc(taille_coag,ci)*2.))) ![m/s]
else !#2
diam_msn = 0.0
endif !#2
! ---The projected area equivalent diameter of our particles of marine snow - Unit should be of size [m]
cc(taille_coag,ci) = (diam_msn**(CSF))*((dt_bio/splitfac_bio)*flux_msn(ci))*secs_pr_day ![m/mmol.m-3.d-1]
! ---The projected area equivalent diameter of our particles of marine snow
cc(taille_coag,ci) = (diam_msn**(CSF))*(dt_bio/splitfac_bio) ![m]
!--» At this points 'taille_coag' is the projected area equivalent diameter in m.
! ---Elaboration of the size variable (diameter), available for the fragmentation eventually (taille_intrm).
if(cc(d4,ci) .gt. cons_min) then
cc(taille_intrm,ci) = cc(taille_coag,ci) + cc(taille_msn,ci)
cc(taille_intrm,ci) = cc(taille_coag,ci) + cc(taille_msn,ci) ![m]
!--» Msn size after coagulation will depend on the increase of size due to coagaluation (taille_coag) &
!--» + the size at the previous time step (taille_msn) of the particles
!-----------------------------------------------------------------------------------
!! ---Safety check
!-----------------------------------------------------------------------------------
if(cc(taille_intrm,ci).lt. 0.0 .OR. &
cc(taille_coag,ci) .lt. 0.0 .OR. &
cc(taille_msn,ci) .lt. 0.0) then
write(*,*)'Taille _intrm - May be equal to 0 :', cc(taille_intrm,ci)
write(*,*)'Taille _coag - May be equal to 0 :', cc(taille_coag,ci)
write(*,*)'Taille _msn - May be equal to 0 :', cc(taille_msn,ci)
endif
else
cc(taille_intrm,ci) = cc(taille_msn,ci)
!--» If there is not need to have a size augmentation due to low concentration of marine snow, the size of previous step remain
cc(taille_intrm,ci) = cc(taille_msn,ci) ![m]
write(*,*)'No size augmentation as [msn] < cons_min - The size stays the same', cc(taille_msn,ci)
endif
else !#1 Size not modified because no coagulation happenned
cc(taille_coag,ci) = 0.0
cc(taille_intrm,ci) = cc(taille_msn,ci)
cc(taille_intrm,ci) = cc(taille_msn,ci) ![m]
write(*,*)'Size not modified because no coagulation happenned'
endif !#1
if(cc(taille_intrm,ci).gt.diam_msn_max) then
cc(taille_intrm,ci) = diam_msn_max
write(*,*)'Maximum size reach - Marine snow size equal :', cc(taille_intrm,ci)
endif
!-----------------------------------------------------------------------------------
!! ---Safety check
!-----------------------------------------------------------------------------------
!--» Msn size can not over pass the maximum diameter specified by the user.
if(cc(taille_intrm,ci).gt.diam_msn_max) then
cc(taille_intrm,ci) = diam_msn_max
write(*,*)'Maximum size reach - Marine snow size equal :', cc(taille_intrm,ci)
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
write(*,*)'Size but no msn'
cc(taille_msn,ci) = 0.0
cc(taille_intrm,ci) = 0.0
endif
if (cc(d4,ci).lt. cons_min)then
cc(taille_msn,ci) = _ZERO_
cc(taille_intrm,ci) = _ZERO_
write(*,*)'Size but no msn at :',ci
endif
!-----------------------------------------------------------------------------------
! FRAGMENTATION
......@@ -2322,13 +2353,13 @@ if (Frag_meth .eq. 1.0) then !Fragmentation happens
kolmog = (kinvis**3/eps(ci))**0.25 ![m]
! ---Minimum size of what composed the marine snow - Be careful to be sure that their is this particles in the field
if (max_flux .eq. flux_P(ci)) then
min_size_msn = size_phy
min_size_msn = Ratio_P*diam_phy !size_phy
elseif (max_flux .eq. flux_D1(ci)) then
min_size_msn = size_dph
min_size_msn = Ratio_D1*diam_dph !size_dph
elseif (max_flux .eq. flux_D2(ci)) then
min_size_msn = size_dzo
min_size_msn = Ratio_D2*diam_dzo !size_dzo
elseif (max_flux .eq. flux_D3(ci)) then
min_size_msn = size_fp
min_size_msn = Ratio_D3* diam_fp ! size_fp
else
write(*,*)'Error for estimation of min_size_msn'
endif
......@@ -2428,12 +2459,7 @@ if (Frag_meth .eq. 1.0) then !Fragmentation happens
!-----------------------------------------------------------------------------------
!! ---Safety check
!-----------------------------------------------------------------------------------
!--» If the size of msn after the cut is lower thant the smallest particle that may composed it,
!--» the minimum size of its constituant are set for the marine snow
if(size_msn_m .lt. min_size_msn)then
write(*,*)'Smaller than the smallest initial particles'
size_msn_m = min_size_msn
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
write(*,*)'Size but no msn'
......@@ -2480,8 +2506,14 @@ endif
!--» Taille_msn --> is the diameter after coag and frag
Rmsn = (rho_msn-densFlu)
write(*,*)'[msn]', cc(d4,ci)
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
if (cc(d4,ci) .lt. 0.2490) then !#1
elseif (cc(d4,ci) .lt. cons_max) then !#1
!--» We set here the concentration limit in the water column from free settling of marine snow vs flocculation settling (Metha 1989)
!! --- Calculated by value from the .nml
......@@ -2530,20 +2562,26 @@ if (cc(d4,ci) .lt. 0.2490) then !#1
write(*,*) 'No settling velocity choice for msn, value set to 0'
w_msn_m = 0.0
endif!#2
else !#1
w_msn_m=-(coef5*cc(d4,ci)**(1.6)) ! Metha 1989
write(*,*)'w_msn_m6',w_msn_m
endif !#1
elseif (cc(d4,ci) .ge. cons_max)then !#1
w_msn_m=-(coef5*cc(d4,ci)**(1.6)) ! Metha 1989
write(*,*)'[msn] > cons_max at :',ci
else
write(*,*)'At Depth : ',ci
write(*,*)'Error with the specification of [msn] at',cc(d4,ci)
w_msn_m = 0.0
endif !#1
!-----------------------------------------------------------------------------------
!! --- Convertion & Comparaison with the Stokes range
!-----------------------------------------------------------------------------------
!--» Now we put everything in /d not matter the scenario.
w_msn_m = w_msn_m/secs_pr_day
! w_msn_m = w_msn_m/secs_pr_day
!--» Comparaison with the Stokes range ! VanRijn 1993
Re_msn = w_msn_m*(cc(taille_msn,ci)/kinvis)
Re_msn = abs(w_msn_m*(cc(taille_msn,ci)/kinvis))
if (Re_msn .lt. 1.0)then !#3
w_msn_m = w_msn_m
......@@ -2551,8 +2589,6 @@ else !#3
w_msn_m = (cc(taille_msn,ci)**0.5)
write(*,*) 'Reynolds number for marine snow > 1'
endif !#3
!!! COMPARER ICI avec ! MESSAGE d'ERREUR SI w < 1/2 Delta z /Delta t
!-----------------------------------------------------------------------------------
......@@ -2560,16 +2596,45 @@ endif !#3
!-----------------------------------------------------------------------------------
!--» let them settling(m/s)
!--» The ncdf output do ws(x,ci) *86400 to have the ouptut in m/j.
!--» The ncdf output do ws(x,ci) *86400 to have the ouptut in m/j.
if (rho_p .gt. densFlu .AND. w_p_m .gt. 0.0) then
w_p_m = -w_p_m
endif
if (rho_dph .gt. densFlu .AND. w_dph_m .gt. 0.0) then
w_dph_m = -w_dph_m
endif
if (rho_dzo .gt. densFlu .AND. w_dzo_m .gt. 0.0) then
w_dzo_m = -w_dzo_m
endif
if (rho_fp .gt. densFlu .AND. w_fp_m .gt. 0.0) then
w_fp_m = -w_fp_m
endif
ws(p,ci) = w_p_m
ws(d1,ci) = w_dph_m
ws(d2,ci) = w_dzo_m
ws(d3,ci) = w_fp_m
!--» Case of marine snow (already /86400)
ws(d4,ci) = w_msn_m
ws(taille_msn,ci) = ws(d4,ci)
!--» Case of marine snow
!--Concentration - Vs CFL limit
ws(d4,ci)= min(w_msn_m,(depth_bio/nlev)/dt_bio)
if (ws(d4,ci) .eq. w_msn_m) then
write(*,*) 'CFL limit is respected'
else
write(*,*) 'CFL Constrain is OVERPASS'
endif
!--Size - Vs CFL limit
ws(taille_msn,ci) = min(w_msn_m,(depth_bio/nlev)/dt_bio)
!--Settling velocity trait - Vs CFL limit
cc(settl_msn,ci)= w_msn_m
ws(settl_msn,ci) = min(w_msn_m,(depth_bio/nlev)/dt_bio)
!-----------------------------------------------------------------------------------
! NET PRIMARY PRODUCTION
......@@ -3106,5 +3171,25 @@ endif !#A
!!cc(taille_msn,ci)= size_msn_m ! Keep in memory for the next step
!size_msn_m = cc(taille_msn,ci)+ cc(aug_si_d4,ci)
!cc(taille_coag,ci)= size_msn_m
!end if !#3
!end if !#3
!--» If the size of msn after the cut is lower thant the smallest particle that may composed it,
!--» the minimum size of its constituant are set for the marine snow
! if(size_msn_m .lt. min_size_msn)then
! write(*,*)'Smaller than the smallest initial particles'
! size_msn_m = min_size_msn
! endif
! write(*,*)'flux_msn(ci) ',flux_msn(ci)
! ---Calcul of the ratios [-]
! Ratio_P = Flux_P(ci)/flux_msn(ci)
! Ratio_D1 = Flux_D1(ci)/flux_msn(ci)
! Ratio_D2 = Flux_D2(ci)/flux_msn(ci)
! Ratio_D3 = Flux_D3(ci)/flux_msn(ci)
Supports Markdown
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