Commit 8f1bf0691bdffff00b8e566407bf1fc49c1b4ad8

Authored by Gwenaelle Gremion
1 parent 3d7dc6b7
Exists in snow

Test avec les equations de VanRijn-beug pour grande tailles

nml/airsea.nml
... ... @@ -97,11 +97,11 @@
97 97 ice_file= 'ice.dat'
98 98 iceatt= 0.90
99 99 wet_mode= 1
100   - heat_method= 2
  100 + heat_method= 0
101 101 const_swr= 100.0
102 102 const_heat= -100.0
103 103 heatflux_file= 'narr_daily_heatflux_ice.dat'
104   - momentum_method= 2
  104 + momentum_method= 0
105 105 const_tx= 0.1
106 106 const_ty= 0.0
107 107 momentumflux_file='narr_hourly_momentumflux.dat'
... ...
nml/gotmrun.nml
... ... @@ -16,11 +16,11 @@
16 16 !
17 17 !-------------------------------------------------------------------------------
18 18 &model_setup
19   - title= "Arctic SCM"
  19 + title= "exp1"
20 20 nlev= 80
21   - dt= 300.
  21 + dt= 120.
22 22 cnpar= 1.0
23   - buoy_method= 2
  23 + buoy_method= 1
24 24 /
25 25  
26 26 !-------------------------------------------------------------------------------
... ... @@ -57,9 +57,9 @@
57 57 !-------------------------------------------------------------------------------
58 58 &time
59 59 timefmt= 2
60   - MaxN= 1200
61   - start= '2004-01-01 00:00:00'
62   - stop= '2004-12-31 00:00:00'
  60 + MaxN= 10
  61 + start= '2008-01-01 00:00:00'
  62 + stop= '2009-05-15 00:00:00'
63 63 /
64 64  
65 65 !-------------------------------------------------------------------------------
... ... @@ -87,9 +87,9 @@
87 87 !-------------------------------------------------------------------------------
88 88 &output
89 89 out_fmt= 2
90   - out_dir= "."
91   - out_fn= "amdgulf"
92   - nsave= 36
  90 + out_dir= "./"
  91 + out_fn= "exptt"
  92 + nsave= 10
93 93 diagnostics= .false.
94 94 mld_method= 2
95 95 diff_k= 1.e-5
... ...
nml/obs.nml
... ... @@ -48,7 +48,7 @@
48 48 !
49 49 !-------------------------------------------------------------------------------
50 50 &sprofile
51   - s_prof_method= 2
  51 + s_prof_method= 1
52 52 s_analyt_method= 2
53 53 z_s1= 25.
54 54 s_1= 31.
... ... @@ -56,9 +56,9 @@
56 56 s_2= 32.
57 57 s_obs_NN= 2.56e-4
58 58 s_prof_file= 'franklin_sprof_ctd.dat'
59   - SRelaxTauM= 1209600.
60   - SRelaxTauB= 1209600.
61   - SRelaxTauS= 1209600.
  59 + SRelaxTauM= 1.e15
  60 + SRelaxTauB= 1.e15
  61 + SRelaxTauS= 1.e15
62 62 SRelaxBott= 0.
63 63 SRelaxSurf= 0.
64 64 /
... ... @@ -110,7 +110,7 @@
110 110 !
111 111 !-------------------------------------------------------------------------------
112 112 &tprofile
113   - t_prof_method= 2
  113 + t_prof_method= 1
114 114 t_analyt_method= 2
115 115 z_t1= 25.
116 116 t_1= -1.0
... ... @@ -118,9 +118,9 @@
118 118 t_2= 0.0
119 119 t_obs_NN= 2.56e-4
120 120 t_prof_file= 'franklin_tprof_ctd.dat'
121   - TRelaxTauM= 1209600.
122   - TRelaxTauB= 1209600.
123   - TRelaxTauS= 1209600.
  121 + TRelaxTauM= 1.e15
  122 + TRelaxTauB= 1.e15
  123 + TRelaxTauS= 1.e15
124 124 TRelaxBott= 0.
125 125 TRelaxSurf= 0.
126 126 /
... ...
src/extras/bio/bio_polynow.F90
... ... @@ -1184,6 +1184,8 @@ end if
1184 1184 REALTYPE :: diam_msn_max,min_size_msn
1185 1185 ! Sedimentation rate msn
1186 1186 REALTYPE :: w_msn_m,Re_msn,Rmsn,CFL
  1187 +
  1188 +REALTYPE :: ceci
1187 1189  
1188 1190 !EOP
1189 1191 !-----------------------------------------------------------------------
... ... @@ -1360,6 +1362,18 @@ diam_dph = sign(abs(((size_dph*2.)*(size_dph*2.)*(size_dph*2.)))**(1./3.),((size
1360 1362 diam_dzo = sign(abs(((size_dzo*4.)*(size_dzo*4.)*(size_dzo*2.)))**(1./3.),((size_dzo*4.)*(size_dzo*4.)*(size_dzo*2.)))
1361 1363  
1362 1364 diam_fp = sign(abs(((size_fp*4.)*(size_fp*4.)*(size_fp*2.)))**(1./3.),((size_fp*4.)*(size_fp*4.)*(size_fp*2.)))
  1365 +
  1366 +
  1367 +!----------Calculation of the CSF factor
  1368 +!--Phytoplankton
  1369 + CSF_phy = (size_phy*2)/sqrt((size_phy*2)*(size_phy*2))
  1370 +!--Dead Phytoplankton
  1371 + CSF_dph = (size_dph*2)/sqrt((size_dph*2)*(size_dph*2))
  1372 +!--Dead Zooplankton
  1373 + CSF_dzo = (size_dzo*2)/sqrt((size_dzo*4)*(size_dzo*4))
  1374 +!--Fecla Pellets
  1375 + CSF_fp = (size_fp*2)/sqrt((size_fp*4)*(size_fp*4))
  1376 +
1363 1377  
1364 1378 !----------Calcul of densities
1365 1379 !Fluid density
... ... @@ -1458,8 +1472,6 @@ else if (Phys_w .eq. 3.0) then !#5
1458 1472 else if (Phys_w .eq. 4.0) then !#5
1459 1473  
1460 1474 !--Phytoplankton
1461   - !-- Calculation of the CSF factor
1462   - CSF_phy = (size_phy*2)/sqrt((size_phy*2)*(size_phy*2))
1463 1475  
1464 1476 if(CSF_phy .ge. 0.0 .and. CSF_phy.lt. 0.4) then !#5.5
1465 1477 CSF_phy = 2.18-(2.09*CSF_phy)
... ... @@ -1474,8 +1486,6 @@ else if (Phys_w .eq. 4.0) then !#5
1474 1486 w_p_m=(1/(18*(kinvis*densFlu)))*(1/CSF_phy)*(Rp)*g*(diam_phy)**2 ![m/s]
1475 1487  
1476 1488 !--Dead phytoplankton
1477   - !-- Calculation of the CSF factor
1478   - CSF_dph = (size_dph*2)/sqrt((size_dph*2)*(size_dph*2))
1479 1489  
1480 1490 if(CSF_dph .ge. 0.0 .and. CSF_dph.lt. 0.4) then !#5.7
1481 1491 CSF_dph = 2.18-(2.09*CSF_dph)
... ... @@ -1489,11 +1499,7 @@ else if (Phys_w .eq. 4.0) then !#5
1489 1499  
1490 1500 w_dph_m=(1/(18*(kinvis*densFlu)))*(1/CSF_dph)*(Rdph)*g*(diam_dph)**2 ![m/s]
1491 1501  
1492   -
1493 1502 !--Dead Zooplankton
1494   - !-- Calculation of the CSF factor
1495   -
1496   - CSF_dzo = (size_dzo*2)/sqrt((size_dzo*4)*(size_dzo*4))
1497 1503  
1498 1504 if(CSF_dzo .ge. 0.0 .and. CSF_dzo.lt. 0.4) then !#5.8
1499 1505 CSF_dzo = 2.18-(2.09*CSF_dzo)
... ... @@ -1507,10 +1513,7 @@ else if (Phys_w .eq. 4.0) then !#5
1507 1513  
1508 1514 w_dzo_m=(1/(18*(kinvis*densFlu)))*(1/CSF_dzo)*(Rdzo)*g*(diam_dzo)**2 ![m/s]
1509 1515  
1510   -
1511 1516 !--Fecal pellets
1512   - !-- Calculation of the CSF factor
1513   - CSF_fp = (size_fp*2)/sqrt((size_fp*4)*(size_fp*4))
1514 1517  
1515 1518 if(CSF_fp .ge. 0.0 .and. CSF_fp.lt. 0.4) then !#5.9
1516 1519 CSF_fp = 2.18-(2.09*CSF_fp)
... ... @@ -1523,7 +1526,120 @@ else if (Phys_w .eq. 4.0) then !#5
1523 1526 endif !#5.9
1524 1527  
1525 1528 w_fp_m=(1/(18*(kinvis*densFlu)))*(1/CSF_fp)*(Rfp)*g*(diam_fp)**2 ![m/s]
1526   -
  1529 +
  1530 +!!--------Depending on CSF + Re
  1531 +else if (Phys_w .eq. 5.0) then !#5
  1532 +
  1533 + !--Phytoplankton
  1534 + if (CSF_phy .eq. 1.0) then ! Considered as a sphere
  1535 + w_p_m = -(((Rp/densFlu)*g*(diam_phy**2))/(18*kinvis)) ![m/s]
  1536 + !! write(*,*) 'CSF_phy = 0 --> Sphere, then Stokes Law ',w_p_m
  1537 + ! Comparaison with the Stokes range ! VanRijn 1993
  1538 + Re_phy = abs(w_p_m*(diam_phy/kinvis))
  1539 + if (Re_phy .gt. 1.0) then
  1540 + w_p_m = -(diam_phy**0.5/secs_pr_day) ![m/s]
  1541 + !! write(*,*) 'Reynolds number for phytoplankton > 1',w_p_m
  1542 + endif
  1543 + else ! No spherical
  1544 + if(diam_phy .le. 0.0001 )then
  1545 + ! If the diameter of phy (d) --> 1 < d <= 100 micrometers (VanRijn 1993)
  1546 + w_p_m =-(((Rp/densFlu)*g*(diam_phy**2))/(18*kinvis)) ![m/s]
  1547 + !! write(*,*) 'CSF .ne. 0.0 --> Diameter of phy : 1< d <= 100 micrometers, settling =' ,w_p_m
  1548 + elseif(diam_phy .gt. 0.0001 .and. diam_phy .le. 0.001)then
  1549 + ! If the diameter of Phy (d) --> 100< d micrometers (VanRijn 1993)
  1550 + w_p_m= ((10.0*kinvis)/diam_phy)*(((1+((0.01*(Rp/densFlu)*g*((diam_phy)**3.0))/(kinvis**2.0)))**0.5)-1) ![m/s]
  1551 + ! w_p_m = -(((10*kinvis)/diam_phy)*&
  1552 + ! ((1+((0.01*(Rp/densFlu)*g*(diam_phy)**3.0) /kinvis**2.0)**0.5) -1)) ![m/s]
  1553 + !! write(*,*) 'CSF .ne. 0.0 --> Diameter of phy : 100< d micrometers, settling =' ,w_p_m
  1554 + elseif(diam_phy .gt. 0.001) then
  1555 + ! If the diameter of Phy (d) --> d > 1000 micrometers (VanRijn 1993)
  1556 + w_p_m = -(1.1*((Rp/densFlu)*g*diam_phy)**0.5) ![m/s]
  1557 + !! write(*,*) 'Diameter of Phy : d > 1000 micrometers, settling =' ,w_p_m
  1558 + endif
  1559 + endif
  1560 +!--Dead Phytoplankton
  1561 + if (CSF_dph .eq. 1.0) then ! Considered as a sphere
  1562 + w_dph_m = -(((Rdph/densFlu)*g*(diam_dph**2))/(18*kinvis)) ![m/s]
  1563 + !! write(*,*) 'CSF_dph = 0 --> Sphere, then Stokes Law ',w_dph_m
  1564 + ! Comparaison with the Stokes range ! VanRijn 1993
  1565 + Re_dph = abs(w_dph_m*(diam_dph/kinvis))
  1566 + if (Re_dph .gt. 1.0) then
  1567 + w_dph_m = -diam_dph**0.5/secs_pr_day ![m/s]
  1568 + !! write(*,*) 'Reynolds number for Dead-phytoplankton > 1',w_dph_m
  1569 + endif
  1570 + else ! No spherical
  1571 + if(diam_dph .le. 0.0001 )then
  1572 + ! If the diameter of dph (d) --> 1 < d <= 100 micrometers (VanRijn 1993)
  1573 + w_dph_m =-(((Rdph/densFlu)*g*(diam_dph**2))/(18*kinvis)) ![m/s]
  1574 + !! write(*,*) 'CSF .ne. 0.0 --> Diameter of dph : 1< d <= 100 micrometers, settling =' ,w_dph_m
  1575 + elseif(diam_dph .gt. 0.0001 .and. diam_dph .le. 0.001 )then
  1576 + ! If the diameter of dph (d) --> 100< d micrometers (VanRijn 1993)
  1577 + w_dph_m= ((10.0*kinvis)/diam_dph)*(((1+((0.01*(Rdph/densFlu)*g*((diam_dph)**3.0))/(kinvis**2.0)))**0.5)-1) ![m/s]
  1578 + ! w_dph_m = -(((10*kinvis)/diam_dph)*&
  1579 + ! ((1+((0.01*(Rdph/densFlu)*g*(diam_dph)**3.0) /kinvis**2.0)**0.5) -1)) ![m/s]
  1580 + !! write(*,*) 'CSF .ne. 0.0 --> Diameter of dph : 100< d micrometers, settling =' ,w_dph_m
  1581 + elseif(diam_dph .gt. 0.001) then
  1582 + ! If the diameter of dph (d) --> d > 1000 micrometers (VanRijn 1993)
  1583 + w_dph_m = -(1.1*((Rdph/densFlu)*g*diam_dph)**0.5) ![m/s]
  1584 + !! write(*,*) 'Diameter of dph : d > 1000 micrometers, settling =' ,w_dph_m
  1585 + endif
  1586 + endif
  1587 +!--Dead Zooplankton
  1588 + if (CSF_dzo .eq. 1.0) then ! Considered as a sphere
  1589 + w_dzo_m = -(((Rdzo/densFlu)*g*(diam_dzo**2))/(18*kinvis)) ![m/s]
  1590 + !! write(*,*) 'CSF_dzo = 0 --> Sphere, then Stokes Law ',w_dzo_m
  1591 + ! Comparaison with the Stokes range ! VanRijn 1993
  1592 + Re_dzo = abs(w_dzo_m*(diam_dzo/kinvis))
  1593 + if (Re_dzo .gt. 1.0) then
  1594 + w_dzo_m = -diam_dzo**0.5/secs_pr_day ![m/s]
  1595 + !! write(*,*) 'Reynolds number for Dead-zooplankton > 1',w_dzo_m
  1596 + endif
  1597 + else ! No spherical
  1598 + if(diam_dzo .le. 0.0001 )then
  1599 + ! If the diameter of dzo (d) --> 1 < d <= 100 micrometers (VanRijn 1993)
  1600 + w_dzo_m =-(((Rdzo/densFlu)*g*(diam_dzo**2))/(18*kinvis)) ![m/s]
  1601 + !! write(*,*) 'CSF .ne. 0.0 --> Diameter of dzo: 1< d <= 100 micrometers, settling =' ,w_dzo_m
  1602 + elseif(diam_dzo .gt. 0.0001 .and. diam_dzo .le. 0.001)then
  1603 + ! If the diameter of dzo (d) --> 100< d micrometers (VanRijn 1993)
  1604 + w_dzo_m= ((10.0*kinvis)/diam_dzo)*(((1+((0.01*(Rdzo/densFlu)*g*((diam_dzo)**3.0))/(kinvis**2.0)))**0.5)-1) ![m/s]
  1605 + ! w_dzo_m = -(((10*kinvis)/diam_dzo)*&
  1606 + ! ((1+((0.01*(Rdzo/densFlu)*g*(diam_dzo)**3.0) /kinvis**2.0)**0.5) -1)) ![m/s]
  1607 + !! write(*,*) 'CSF .ne. 0.0 --> Diameter of dzo : 100< d micrometers, settling =' ,w_dzo_m
  1608 +
  1609 + elseif(diam_dzo .gt. 0.001) then
  1610 + ! If the diameter of Dzo (d) --> d > 1000 micrometers (VanRijn 1993)
  1611 + w_dzo_m = -(1.1*((Rdzo/densFlu)*g*diam_dzo)**0.5) ![m/s]
  1612 + !! write(*,*) 'Diameter of Dzo : d > 1000 micrometers, settling =' ,w_dzo_m
  1613 + endif
  1614 + endif
  1615 + !--Fecal pellets
  1616 + if (CSF_fp .eq. 1.0) then ! Considered as a sphere
  1617 + w_fp_m = -(((Rfp/densFlu)*g*(diam_fp**2))/(18*kinvis)) ![m/s]
  1618 + !! write(*,*) 'CSF_fp = 0 --> Sphere, then Stokes Law ',w_fp_m
  1619 + ! Comparaison with the Stokes range ! VanRijn 1993
  1620 + Re_fp = abs(w_fp_m*(diam_fp/kinvis))
  1621 + if (Re_fp .gt. 1.0) then
  1622 + w_fp_m = -diam_fp**0.5/secs_pr_day ![m/s]
  1623 + !! write(*,*) 'Reynolds number for Fecal pellets > 1',w_fp_m
  1624 + endif
  1625 + else ! No spherical
  1626 + if( diam_fp .le. 0.0001 )then
  1627 + ! If the diameter of fp (d) --> 1 < d <= 100 micrometers (VanRijn 1993)
  1628 + w_fp_m =-(((Rfp/densFlu)*g*(diam_fp**2))/(18*kinvis)) ![m/s]
  1629 + !! write(*,*) 'CSF .ne. 0.0 --> Diameter of fp: 1< d <= 100 micrometers, settling =' ,w_fp_m
  1630 + elseif(diam_fp .gt. 0.0001 .and. diam_fp .le. 0.001)then
  1631 + ! If the diameter of fp (d) --> 100< d micrometers (VanRijn 1993)
  1632 + w_fp_m= ((10.0*kinvis)/diam_fp)*(((1+((0.01*(Rfp/densFlu)*g*((diam_fp)**3.0))/(kinvis**2.0)))**0.5)-1) ![m/s]
  1633 + ! w_fp_m = -(((10*kinvis)/diam_fp)*&
  1634 + ! ((1+((0.01*(Rfp/densFlu)*g*(diam_fp)**3.0) /kinvis**2.0)**0.5) -1)) ![m/s]
  1635 + !! write(*,*) 'CSF .ne. 0.0 --> Diameter of fp : 100< d micrometers, settling =' ,w_fp_m
  1636 + elseif(diam_fp .gt. 0.001) then
  1637 + ! If the diameter of Fp (d) --> d > 1000 micrometers (VanRijn 1993)
  1638 + ! w_fp_m = -(1.1*((Rfp/densFlu)*g*diam_fp)**0.5) ![m/s]
  1639 + w_fp_m = -(1.1*((Rfp/densFlu)*g*diam_fp)**0.5) ![m/s]
  1640 + !! write(*,*) 'Diameter of Fp : d > 1000 micrometers, settling =' ,w_fp_m
  1641 + endif
  1642 + endif
1527 1643 else !#5
1528 1644 write(*,*) 'No specification for settling velocity of particles,set to nml values'
1529 1645 w_p_m = w_p ![m/s]
... ... @@ -1531,8 +1647,9 @@ else !#5
1531 1647 w_dzo_m = w_dzo![m/s]
1532 1648 w_fp_m = w_fp ![m/s]
1533 1649 endif !#5
1534   -
1535   - ! Comparaison with the Stokes range ! VanRijn 1993
  1650 +
  1651 +If (Phys_w .ne. 5.0) then
  1652 +! ! Comparaison with the Stokes range ! VanRijn 1993
1536 1653 Re_phy = abs(w_p_m*(diam_phy/kinvis))
1537 1654 Re_dph = abs(w_dph_m*(diam_dph/kinvis))
1538 1655 Re_dzo = abs(w_dzo_m*(diam_dzo/kinvis))
... ... @@ -1557,6 +1674,9 @@ if (Re_fp .gt. 1.0)then !#9
1557 1674 w_fp_m = -diam_fp**0.5/secs_pr_day ![m/s]
1558 1675 write(*,*) 'Reynolds number for fecal pellets > 1'
1559 1676 endif !#9
  1677 +
  1678 +endif
  1679 +
1560 1680  
1561 1681 !Everyone will settle at the end of the loop (see below)
1562 1682  
... ... @@ -2200,7 +2320,7 @@ endif !#1
2200 2320 !-----------------------------------------------------------------------------------
2201 2321 ! SIZE Marine Snow with aggregation
2202 2322 !-----------------------------------------------------------------------------------
2203   -
  2323 +
2204 2324 ! ---Maximum size at dissagregation (Alldredge 1990)
2205 2325 if(dm_msn .eq. 1.0) then
2206 2326 diam_msn_max = coef1*(eps(ci))**(-coef2) ! [m]
... ... @@ -2322,7 +2442,6 @@ endif !#1
2322 2442 ! cc(size_intrm,ci) = _ZERO_
2323 2443 ! write(*,*)'Size but no msn at :',ci
2324 2444 ! endif
2325   -
2326 2445 !-----------------------------------------------------------------------------------
2327 2446 ! FRAGMENTATION
2328 2447 ! -----------
... ... @@ -2485,7 +2604,9 @@ endif
2485 2604 ! We want to use a calculated settling velocity for the msn
2486 2605 ! Settling velocity send to the advection scheme as well as the size of the msn
2487 2606 !-------------------------------------------------------------------------------------
2488   -
  2607 +write(*,*)'Depth',ci
  2608 +write(*,*)'[msn]',cc(d4,ci)
  2609 +write(*,*) 'Diameter of msn :',cc(size_msn,ci)
2489 2610 !--» cc(size_msn,ci) --> is the diameter after coag and frag
2490 2611 Rmsn = (rho_msn-densFlu)
2491 2612  
... ... @@ -2538,52 +2659,96 @@ endif
2538 2659 w_msn_m=(1/(18*(kinvis*densFlu)))*(1/CSF_msn)*(Rmsn)*g*(cc(size_msn,ci))**2 ![m/s]
2539 2660  
2540 2661 elseif (w_msnow .eq. 5.0) then !#2!
2541   -!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
2542   - if( diam_msn_us .gt. 0.000001 .and. diam_msn_us .le. 0.0001 )then
  2662 + ! if( diam_msn_us .gt. 0.000001 .and. diam_msn_us .le. 0.0001 )then
  2663 + write(*,*) 'Diameter of msn :',cc(size_msn,ci)
  2664 +
  2665 + if(cc(size_msn,ci) .le. 0.0001 )then
2543 2666 ! If the diameter of Msn (d) --> 1< d <= 100 micrometers (VanRijn 1993)
2544   - w_msn_m = -(((2.65-1)*g*(cc(size_msn,ci)**2))/18*kinvis) ![m/s]
2545   - write(*,*) 'Diameter of msn : 1< d <= 100 micrometers, settling =' ,w_msn_m
  2667 + w_msn_m = -(((Rmsn/densFlu)*g*(cc(size_msn,ci)**2))/(18*kinvis)) ![m/s]
  2668 + write(*,*) 'Diameter of msn : d <= 100 micrometers, settling (Stokes Law is used !) =' ,w_msn_m
2546 2669  
2547   - elseif( diam_msn_us .gt. 0.0001 .and. diam_msn_us .lt. 0.001 ) then
  2670 +
  2671 + elseif(cc(size_msn,ci) .gt. 0.0001 .and. cc(size_msn,ci) .lt. 0.001 ) then
2548 2672 ! If the diameter of Msn (d) --> 100< d <= 1000 micrometers (VanRijn 1993)
2549 2673 w_msn_m = -(((10*kinvis)/cc(size_msn,ci))*&
2550   - ((1+((0.01*(2.65-1)*g*(cc(size_msn,ci))**3.0) /kinvis**2.0)**0.5) -1)) ![m/s]
  2674 + ((1+((0.01*(Rmsn/densFlu)*g*cc(size_msn,ci)**3.0) /kinvis**2.0)**0.5) -1)) ![m/s]
2551 2675 write(*,*) 'Diameter of msn : 100< d <= 1000 micrometers, settling =' ,w_msn_m
2552 2676  
2553   - elseif( diam_msn_us .gt. 0.001) then
  2677 + else !(cc(size_msn,ci).gt. 0.001) then
2554 2678 ! If the diameter of Msn (d) --> d > 1000 micrometers (VanRijn 1993)
2555   - w_msn_m = -(1.1*((2.65-1)*g*cc(size_msn,ci))**0.5)/secs_pr_day ![m/s]
  2679 + w_msn_m = -(1.1*((Rmsn/densFlu)*g*cc(size_msn,ci))**0.5) ![m/s]
2556 2680 write(*,*) 'Diameter of msn : d > 1000 micrometers, settling =' ,w_msn_m
2557   -
2558   - else
2559   - w_msn_m =- (g*((cc(size_msn,ci))**2)*Rmsn)/(18*dynvis) ![m/s]
2560   - write(*,*) 'Diam.of msn < than 1 micrometers: Stokes Law is used ! ',w_msn_m
  2681 +
2561 2682 endif
2562 2683  
2563 2684  
2564   -
2565 2685 elseif (w_msnow .eq. 6.0) then !#2!
  2686 + if (CSF .eq. 1.0) then ! Considered as a sphere
  2687 + w_msn_m = -(((Rmsn/densFlu)*g*(cc(size_msn,ci)**2))/(18*kinvis)) ![m/s]
  2688 + write(*,*) 'CSF_msn = 0 --> Sphere, then Stokes Law ',w_msn_m
  2689 + ! Comparaison with the Stokes range ! VanRijn 1993
  2690 + Re_msn = abs(w_msn_m*(cc(size_msn,ci)/kinvis))
  2691 + if (Re_msn .gt. 1.0) then
  2692 + w_msn_m = -(cc(size_msn,ci)**0.5/secs_pr_day) ![m/s]
  2693 + write(*,*) 'Reynolds number for msn > 1',w_msn_m
  2694 + endif
  2695 + else ! No spherical
  2696 + if(cc(size_msn,ci).le. 0.0001 )then
  2697 + ! If the diameter of msn (d) --> 1 < d <= 100 micrometers (VanRijn 1993)
  2698 + w_msn_m =-(((Rmsn/densFlu)*g*(cc(size_msn,ci)**2))/(18*kinvis)) ![m/s]
  2699 + write(*,*) 'CSF .ne. 0.0 --> Diameter of msn: 1< d <= 100 micrometers, settling =' ,w_msn_m
  2700 + elseif(cc(size_msn,ci) .gt. 0.0001 )then
  2701 + ! If the diameter of dzo (d) --> 100< d micrometers (VanRijn 1993)
  2702 + w_msn_m = -(((10*kinvis)/cc(size_msn,ci))*&
  2703 + ((1+((0.01*(Rmsn/densFlu)*g*(cc(size_msn,ci))**3.0) /kinvis**2.0)**0.5) -1)) ![m/s]
  2704 + write(*,*) 'CSF .ne. 0.0 --> Diameter of msn : 100< d micrometers, settling =' ,w_msn_m
  2705 + endif
  2706 + endif
2566 2707  
2567   - if(cc(size_msn,ci) .gt. 0.000001 .and. cc(size_msn,ci) .le. 0.0001 )then
2568   - ! If the diameter of Msn (d) --> 1< d <= 100 micrometers (VanRijn 1993)
2569   - w_msn_m = -(((2.65-1)*g*(cc(size_msn,ci)**2))/18*kinvis) ![m/s]
2570   - write(*,*) 'Diameter of msn : 1< d <= 100 micrometers, settling =' ,w_msn_m
2571 2708  
2572   - elseif(cc(size_msn,ci) .gt. 0.0001 .and. cc(size_msn,ci) .lt. 0.001 ) then
2573   - ! If the diameter of Msn (d) --> 100< d <= 1000 micrometers (VanRijn 1993)
2574   - w_msn_m = -(((10*kinvis)/cc(size_msn,ci))*&
2575   - ((1+((0.01*(2.65-1)*g*(cc(size_msn,ci))**3.0) /kinvis**2.0)**0.5) -1)) ![m/s]
2576   - write(*,*) 'Diameter of msn : 100< d <= 1000 micrometers, settling =' ,w_msn_m
  2709 + elseif (w_msnow .eq. 7.0) then !#2!
  2710 +
  2711 + ! If the diameter of Msn (d) --> d <= 100 micrometers (Komar 1978)
  2712 + if(cc(size_msn,ci) .le. 0.0001 )then
  2713 +
  2714 + if(CSF .ge. 0.0 .and. CSF.lt. 0.4) then !#2.3
  2715 + CSF_msn = 2.18-(2.09*CSF)
  2716 + else if (CSF .ge. 0.4 .and. CSF .lt. 0.8) then!#2.3
  2717 + CSF_msn = 0.946*(CSF)**(-0.378)
  2718 + else if (CSF .ge. 0.8 .and. CSF .le. 1.0) then!#2.3
  2719 + CSF_msn = 1.0 ! consider as a sphere
  2720 + else !#2.3
  2721 + write (*,*) 'Error in CSF values for Marine snow : CSF .NE. [0-1]'
  2722 + endif !#2.3
  2723 + w_msn_m=(1/(18*(kinvis*densFlu)))*(1/CSF_msn)*(Rmsn/densFlu)*g*(cc(size_msn,ci))**2 ![m/s]
  2724 + write(*,*) 'Diameter of msn : d <= 100 micrometers, settling (Stokes Law is used !) =' ,w_msn_m
  2725 +
  2726 +! If the diameter of Msn (d) --> 100< d <= 1000 micrometers (VanRijn 1993)
  2727 +
  2728 + elseif(cc(size_msn,ci) .gt. 0.0001 .and. cc(size_msn,ci) .lt. 0.001 ) then
  2729 + w_msn_m= ((10.0*kinvis)/cc(size_msn,ci))*(((1+((0.01*(Rmsn/densFlu)*g*(cc(size_msn,ci)**3.0))/(kinvis**2.0)))**0.5)-1) ![m/s]
  2730 + write(*,*) 'Diameter of msn : 100< d <= 1000 micrometers, settling =' ,w_msn_m
  2731 + ! w_msn_m = -(((10.0*kinvis)/cc(size_msn,ci))* &
  2732 + ! ((1+((0.01*(Rmsn/densFlu)*g*(cc(size_msn,ci)**3.0))/kinvis**2.0)**0.5) -1))
  2733 +
  2734 +! If the diameter of Msn (d) --> d > 1000 micrometers (VanRijn 1993)
  2735 + elseif (cc(size_msn,ci) .ge. 0.001 .and. cc(size_msn,ci) .le. diam_msn_max) then
  2736 + ! w_msn_m = -(1.1*(((Rmsn/densFlu)*g*(cc(size_msn,ci)))**0.5)) ![m/s]
  2737 + ! w_msn_m = -(((cc(size_msn,ci))**(1/2))/secs_pr_day)
  2738 + w_msn_m = -(1.1*(((Rmsn/densFlu)*g*cc(size_msn,ci))**(1/2)))
  2739 +
  2740 + ! write(*,*) 'Diameter of msn : d > 1000 micrometers, settling =' ,w_msn_m
  2741 + else
  2742 + w_msn_m = w_msn
  2743 + write(*,*) ' Diam_msn_max is Overpass, value form nml is used',w_msn_m
  2744 + endif
  2745 +
  2746 +if (abs(w_msn_m) .gt. (depth_bio/nlev)/dt_bio )then
  2747 + write(*,*) 'CFL Constrain is OVERPASS for W-msn'
  2748 +w_msn_m = w_msn
  2749 + write(*,*) 'settling',w_msn_m
  2750 +endif
2577 2751  
2578   - elseif(cc(size_msn,ci) .gt. 0.001) then
2579   - ! If the diameter of Msn (d) --> d > 1000 micrometers (VanRijn 1993)
2580   - w_msn_m = -(1.1*((2.65-1)*g*cc(size_msn,ci))**0.5) ![m/s]
2581   - write(*,*) 'Diameter of msn : d > 1000 micrometers, settling =' ,w_msn_m
2582   -
2583   - else
2584   - w_msn_m =- (g*((cc(size_msn,ci))**2)*Rmsn)/(18*dynvis) ![m/s]
2585   - write(*,*) 'Diam.of msn < than 1 micrometers: Stokes Law is used ! ',w_msn_m
2586   - endif
2587 2752  
2588 2753 else !#2!
2589 2754 !if (w_msnow .ne. 0.0 .and. w_msnow .ne. 1.0 .and. w_msnow .ne. 2.0 .and. w_msnow .ne. 3.0 .and. w_msnow .ne. 4.0 .and.w_msnow .ne. 5.0 .and. w_msnow .ne. 6.0 ) then
... ... @@ -2593,31 +2758,33 @@ endif
2593 2758  
2594 2759 elseif (cc(d4,ci) .ge. cons_max)then !#1
2595 2760 w_msn_m=-(coef3*cc(d4,ci)**(1.6))/secs_pr_day ![m/s] ! Metha 1989
2596   -! write(*,*)'[msn] > cons_max at :',ci
  2761 + write(*,*)'[msn] > cons_max at :',ci
2597 2762  
2598 2763 else !#1
2599 2764 write(*,*)'At Depth : ',ci
2600 2765 write(*,*)'Error with the specification of [msn] that equal : ',cc(d4,ci)
2601 2766 w_msn_m = 0.0
  2767 + stop "Nan for cc(d4,ci)"
2602 2768  
2603 2769 endif !#1
2604 2770  
2605 2771 !-----------------------------------------------------------------------------------
2606 2772 !! --- Convertion & Comparaison with the Stokes range
2607 2773 !-----------------------------------------------------------------------------------
2608   -
  2774 +!if (w_msnow .ne. 6.0) then
2609 2775 !--» Comparaison with the Stokes range ! VanRijn 1993
2610   - Re_msn = abs(w_msn_m*(cc(size_msn,ci)/kinvis))
  2776 +! Re_msn = abs(w_msn_m*(cc(size_msn,ci)/kinvis))
2611 2777  
2612   -if (Re_msn .lt. 1.0)then !#3
2613   - ! write(*,*) 'Reynolds number for marine snow allow to respect Stockes range < 1'
2614   -else !#3
2615   - write(*,*) 'Reynolds number for marine snow > 1'
2616   -w_msn_m = -(cc(size_msn,ci)**0.5/secs_pr_day) ![m/s]
2617   - if (w_msnow .eq. 1.0 .or. w_msnow .eq. 2.0 .or. w_msnow .eq. 3.0 )then
2618   - write(*,*) 'Reynolds number for marine snow > 1, another settling velocity scheme should be chose for msn'
2619   - endif
2620   -endif !#3
  2778 +! if (Re_msn .lt. 1.0)then !#3
  2779 + ! write(*,*) 'Reynolds number for marine snow allow to respect Stockes range < 1'
  2780 +! else !#3
  2781 +! write(*,*) 'Reynolds number for marine snow > 1'
  2782 +! w_msn_m = -(cc(size_msn,ci)**0.5/secs_pr_day) ![m/s]
  2783 +! if (w_msnow .eq. 1.0 .or. w_msnow .eq. 2.0 .or. w_msnow .eq. 3.0 )then
  2784 +! write(*,*) 'Reynolds number for marine snow > 1, another settling velocity scheme should be chose for msn'
  2785 +! endif
  2786 +! endif !#3
  2787 +!endif
2621 2788  
2622 2789 !-----------------------------------------------------------------------------------
2623 2790 ! ADVECTION (WS) OF THE PARTICLES
... ... @@ -2645,7 +2812,7 @@ endif !#3
2645 2812 endif
2646 2813  
2647 2814 if (rho_msn .gt. densFlu .AND. w_msn_m .gt. 0.0) then
2648   - write(*,*)' Msn density > dens_fluid + W_msn_m > 0 --> Settling go up, but is constrain to do down !'
  2815 + write(*,*)' Msn density > dens_fluid + W_msn_m > 0 --> Settling go up, but is constrain to go down !'
2649 2816 w_msn_m = -w_msn_m
2650 2817 endif
2651 2818  
... ... @@ -2662,32 +2829,37 @@ endif !#3
2662 2829 !--» Case of Phytoplankton
2663 2830 if (abs(w_p_m) .gt. (depth_bio/nlev)/dt_bio )then
2664 2831 write(*,*) 'CFL Constrain is OVERPASS for W-Phy'
  2832 +stop
2665 2833 endif
2666 2834 ws(p,ci) = w_p_m
2667 2835  
2668 2836 !--» Case of Dead_Phytoplankton
2669 2837 if (abs(w_dph_m) .gt. (depth_bio/nlev)/dt_bio )then
2670 2838 write(*,*) 'CFL Constrain is OVERPASS for W-Dph'
  2839 +stop
2671 2840 endif
2672 2841 ws(d1,ci) = w_dph_m
2673 2842 !--» Case of Dead Zooplankton
2674 2843 if (abs(w_dzo_m) .gt. (depth_bio/nlev)/dt_bio )then
2675 2844 write(*,*) 'CFL Constrain is OVERPASS for W-dzo'
  2845 +stop
2676 2846 endif
2677 2847 ws(d2,ci) = w_dzo_m
2678 2848 !--» Case of Fecal pellets
2679 2849 if (abs(w_fp_m) .gt. (depth_bio/nlev)/dt_bio )then
2680 2850 write(*,*) 'CFL Constrain is OVERPASS for W-fp'
  2851 +stop
2681 2852 endif
2682 2853 ws(d3,ci) = w_fp_m
2683 2854 !--» Case of marine snow
2684 2855 if (abs(w_msn_m) .gt. (depth_bio/nlev)/dt_bio )then
2685 2856 write(*,*) 'CFL Constrain is OVERPASS for W-msn'
  2857 +!stop
2686 2858 endif
2687 2859 ws(d4,ci) = w_msn_m
2688 2860  
2689 2861 !--Size of marine snow
2690   -ws(size_msn,ci) = ws(d4,ci)
  2862 +ws(size_msn,ci) = ws(d4,ci)
2691 2863  
2692 2864  
2693 2865 !-----------------------------------------------------------------------------------
... ... @@ -2745,5 +2917,8 @@ ws(size_msn,ci) = ws(d4,ci)
2745 2917  
2746 2918 !-----------------------------------------------------------------------
2747 2919 ! Copyright by the GOTM-team under the GNU Public License - www.gnu.org
2748   -!-----------------------------------------------------------------------
  2920 +!-----------------------------------------------------------------------
  2921 +
  2922 +
  2923 +
2749 2924  
... ...