Dany Dumont / gotm_ismer

Commit 764531e62b643048c17faad187a81afe9676f08d

Authored by Gwenaelle Gremion
1 parent 8f1bf069
Exists in snow

Maybe the last one...

Showing 2 changed files with 200 additions and 213 deletions   Show diff stats
Inline Side-by-side
src/extras/bio/bio_polynow.F90
  Diff comments   View file @764531e
... ... @@ -1078,7 +1078,8 @@ end if
1078 1078 REALTYPE :: q1,q2
1079 1079 REALTYPE :: fac,min67
1080 1080 !! Size variation for each particle type
1081   - REALTYPE ::size_phy,size_dph,size_dzo,size_fp
  1081 + REALTYPE ::size_phy,size_dph,size_dzo,size_fp
  1082 + REALTYPE ::size_phy_max,size_dph_max,size_dzo_max,size_fp_max
1082 1083 REALTYPE ::diam_phy,diam_dph,diam_dzo,diam_fp,diam_msn
1083 1084 !Collision & Aggregation
1084 1085 !! Sedimentation rate for each particle type
... ... @@ -1086,7 +1087,8 @@ end if
1086 1087 REALTYPE ::densFlu
1087 1088 REALTYPE ::Rp,Rdph,Rdzo,Rfp
1088 1089 REALTYPE ::Re_phy,Re_dph,Re_dzo,Re_fp
1089   - REALTYPE ::CSF_phy,CSF_dph,CSF_dzo,CSF_fp,CSF_msn
  1090 + REALTYPE ::CSF_phy,CSF_dph,CSF_dzo,CSF_fp,CSF_msn
  1091 + REALTYPE ::CSF_phy_calc,CSF_dph_calc,CSF_dzo_calc,CSF_fp_calc,CSF_msn_calc
1090 1092 !! Maximum loss for aggregation for each particle type
1091 1093 REALTYPE ::pprime
1092 1094 REALTYPE ::d1rime
... ... @@ -1345,7 +1347,7 @@ endif !#4
1345 1347  
1346 1348  
1347 1349 !-------------------------------------------------------------------------------------
1348   -!! Sedimentation rate for each particle type
  1350 +!! Usuefull parameters of Sedimentation rate for each particle type
1349 1351 !-------------------------------------------------------------------------------------
1350 1352  
1351 1353 !----------Definition of the equivalent diameter of our particles [Bagheri 2015]-Unit of size [m]
... ... @@ -1363,16 +1365,59 @@ diam_dzo = sign(abs(((size_dzo*4.)*(size_dzo*4.)*(size_dzo*2.)))**(1./3.),((size
1363 1365  
1364 1366 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 1367  
  1368 +!--» Does CFL limit
  1369 + CFL = (depth_bio/nlev)/dt_bio
1366 1370  
1367 1371 !----------Calculation of the CSF factor
1368 1372 !--Phytoplankton
1369 1373 CSF_phy = (size_phy*2)/sqrt((size_phy*2)*(size_phy*2))
  1374 + if(CSF_phy .ge. 0.0 .and. CSF_phy.lt. 0.4) then !#5.5
  1375 + CSF_phy_calc = 2.18-(2.09*CSF_phy)
  1376 + else if (CSF_phy .ge. 0.4 .and. CSF_phy .lt. 0.8) then !#5.6
  1377 + CSF_phy_calc = 0.946*(CSF_phy)**(-0.378)
  1378 + else if (CSF_phy .ge. 0.8 .and. CSF_phy .le. 1.0) then !#5.6
  1379 + CSF_phy_calc = 1.0 ! consider as a sphere
  1380 + else !#5.6
  1381 + write (*,*) 'Error in CSF values for phytoplankton : CSF .NE. [0-1]', CSF_phy
  1382 + endif !#5.6
  1383 +
1370 1384 !--Dead Phytoplankton
1371   - CSF_dph = (size_dph*2)/sqrt((size_dph*2)*(size_dph*2))
  1385 + CSF_dph = (size_dph*2)/sqrt((size_dph*2)*(size_dph*2))
  1386 +write(*,*) 'CSF_dph',CSF_dph
  1387 + if(CSF_dph .ge. 0.0 .and. CSF_dph.lt. 0.4) then !#5.7
  1388 + CSF_dph_calc = 2.18-(2.09*CSF_dph)
  1389 + else if (CSF_dph .ge. 0.4 .and. CSF_dph .lt. 0.8) then !#5.7
  1390 + CSF_dph_calc = 0.946*(CSF_dph)**(-0.378)
  1391 + else if (CSF_dph .ge. 0.8 .and. CSF_dph .le. 1.0) then !#5.7
  1392 + CSF_dph_calc = 1.0 ! consider as a sphere
  1393 + else !#5.7
  1394 + write (*,*) 'Error in CSF values for Dead phytoplankton : CSF .NE. [0-1]',CSF_dph
  1395 + endif !#5.7
  1396 +
1372 1397 !--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))
  1398 + CSF_dzo = (size_dzo*2)/sqrt((size_dzo*4)*(size_dzo*4))
  1399 + if(CSF_dzo .ge. 0.0 .and. CSF_dzo.lt. 0.4) then !#5.8
  1400 + CSF_dzo_calc = 2.18-(2.09*CSF_dzo)
  1401 + else if (CSF_dzo .ge. 0.4 .and. CSF_dzo .lt. 0.8) then !#5.8
  1402 + CSF_dzo_calc = 0.946*(CSF_dzo)**(-0.378)
  1403 + else if (CSF_dzo .ge. 0.8 .and. CSF_dzo .le. 1.0) then !#5.8
  1404 + CSF_dzo_calc = 1.0 ! consider as a sphere
  1405 + else !#5.8
  1406 + write (*,*) 'Error in CSF values for dead zooplankton : CSF .NE. [0-1]',CSF_dzo
  1407 + endif !#5.8
  1408 +
  1409 +!--Fecal Pellets
  1410 + CSF_fp = (size_fp*2)/sqrt((size_fp*4)*(size_fp*4))
  1411 +write(*,*) 'CSF_fp',CSF_fp
  1412 + if(CSF_fp .ge. 0.0 .and. CSF_fp.lt. 0.4) then !#5.9
  1413 + CSF_fp_calc = 2.18-(2.09*CSF_fp)
  1414 + else if (CSF_fp .ge. 0.4 .and. CSF_fp .lt. 0.8) then !#5.9
  1415 + CSF_fp_calc = 0.946*(CSF_fp)**(-0.378)
  1416 + else if (CSF_fp .ge. 0.8 .and. CSF_fp .le. 1.0) then !#5.9
  1417 + CSF_fp_calc = 1.0 ! consider as a sphere
  1418 + else !#5.9
  1419 + write (*,*) 'Error in CSF values for fecal pellets : CSF .NE. [0-1]',CSF_fp
  1420 + endif !#5.9
1376 1421  
1377 1422  
1378 1423 !----------Calcul of densities
... ... @@ -1382,7 +1427,41 @@ densFlu = rho_0
1382 1427 Rp = (rho_p-densFlu)
1383 1428 Rdph= (rho_dph-densFlu)
1384 1429 Rdzo =(rho_dzo-densFlu)
1385   -Rfp =(rho_fp-densFlu)
  1430 +Rfp =(rho_fp-densFlu)
  1431 +Rmsn = (rho_msn-densFlu)
  1432 +
  1433 + !-----------------------------------------------------------------------------------
  1434 + !! ---Safety check
  1435 + !-----------------------------------------------------------------------------------
  1436 +
  1437 +!----------Estimation of the maximum size allowed for particles
  1438 + !--» Comparaison with CFL limit
  1439 + size_phy_max = ((CFL/0.93)**2)*(1/(g*(Rp/densFlu)))
  1440 + size_dph_max = ((CFL/0.93)**2)*(1/(g*(Rdph/densFlu)))
  1441 + size_dzo_max = ((CFL/0.93)**2)*(1/(g*(Rdzo/densFlu)))
  1442 + size_fp_max = ((CFL/0.93)**2)*(1/(g*(Rfp/densFlu)))
  1443 +
  1444 +!---------- Attribution of the good size
  1445 +if (size_phy .gt. size_phy_max) then
  1446 + size_phy = size_phy_max
  1447 + write(*,*) 'Size of Phy OVERPASSES Max.Size of Phy to respect CFL limit'
  1448 +endif
  1449 +if (size_dph .gt. size_dph_max) then
  1450 + size_dph = size_dph_max
  1451 + write(*,*) 'Size of Dph OVERPASSES Max.Size of Dph to respect CFL limit'
  1452 +endif
  1453 +if (size_dzo .gt. size_dzo_max) then
  1454 + size_dzo = size_dzo_max
  1455 + write(*,*) 'Size of Dzo OVERPASSES Max.Size of Dzo to respect CFL limit'
  1456 +endif
  1457 +if (size_fp .gt. size_fp_max) then
  1458 + size_fp = size_fp_max
  1459 + write(*,*) 'Size of Fp OVERPASSES Max.Size of FP to respect CFL limit'
  1460 +endif
  1461 +
  1462 +!-------------------------------------------------------------------------------------
  1463 +!! Calculation of Sedimentation rates
  1464 +!-------------------------------------------------------------------------------------
1386 1465  
1387 1466 !!--------!Calculated by value from the .nml
1388 1467 if (Phys_w .eq. 0.0 ) then !#5
... ... @@ -1472,68 +1551,29 @@ else if (Phys_w .eq. 3.0) then !#5
1472 1551 else if (Phys_w .eq. 4.0) then !#5
1473 1552  
1474 1553 !--Phytoplankton
1475   -
1476   - if(CSF_phy .ge. 0.0 .and. CSF_phy.lt. 0.4) then !#5.5
1477   - CSF_phy = 2.18-(2.09*CSF_phy)
1478   - else if (CSF_phy .ge. 0.4 .and. CSF_phy .lt. 0.8) then !#5.6
1479   - CSF_phy = 0.946*(CSF_phy)**(-0.378)
1480   - else if (CSF_phy .ge. 0.8 .and. CSF_phy .le. 1.0) then !#5.6
1481   - CSF_phy = 1.0 ! consider as a sphere
1482   - else !#5.6
1483   - write (*,*) 'Error in CSF values for phytoplankton : CSF .NE. [0-1]', CSF_phy
1484   - endif !#5.6
1485   -
1486   - w_p_m=(1/(18*(kinvis*densFlu)))*(1/CSF_phy)*(Rp)*g*(diam_phy)**2 ![m/s]
  1554 + ! w_p_m=(1/(18*(kinvis*densFlu)))*(1/CSF_phy_calc)*(Rp)*g*(diam_phy)**2 ![m/s]
  1555 + w_p_m=(1/(18*(kinvis*densFlu)))*(1/CSF_phy_calc)*(Rp/densFlu)*g*(diam_phy)**2 ![m/s]
1487 1556  
1488 1557 !--Dead phytoplankton
1489   -
1490   - if(CSF_dph .ge. 0.0 .and. CSF_dph.lt. 0.4) then !#5.7
1491   - CSF_dph = 2.18-(2.09*CSF_dph)
1492   - else if (CSF_dph .ge. 0.4 .and. CSF_dph .lt. 0.8) then !#5.7
1493   - CSF_dph = 0.946*(CSF_dph)**(-0.378)
1494   - else if (CSF_dph .ge. 0.8 .and. CSF_dph .le. 1.0) then !#5.7
1495   - CSF_dph = 1.0 ! consider as a sphere
1496   - else !#5.7
1497   - write (*,*) 'Error in CSF values for Dead phytoplankton : CSF .NE. [0-1]',CSF_dph
1498   - endif !#5.7
1499   -
1500   - w_dph_m=(1/(18*(kinvis*densFlu)))*(1/CSF_dph)*(Rdph)*g*(diam_dph)**2 ![m/s]
  1558 + ! w_dph_m=(1/(18*(kinvis*densFlu)))*(1/CSF_dph_calc)*(Rdph)*g*(diam_dph)**2 ![m/s]
  1559 + w_dph_m=(1/(18*(kinvis*densFlu)))*(1/CSF_dph_calc)*(Rdph/densFlu)*g*(diam_dph)**2 ![m/s]
1501 1560  
1502 1561 !--Dead Zooplankton
1503   -
1504   - if(CSF_dzo .ge. 0.0 .and. CSF_dzo.lt. 0.4) then !#5.8
1505   - CSF_dzo = 2.18-(2.09*CSF_dzo)
1506   - else if (CSF_dzo .ge. 0.4 .and. CSF_dzo .lt. 0.8) then !#5.8
1507   - CSF_dzo = 0.946*(CSF_dzo)**(-0.378)
1508   - else if (CSF_dzo .ge. 0.8 .and. CSF_dzo .le. 1.0) then !#5.8
1509   - CSF_dzo = 1.0 ! consider as a sphere
1510   - else !#5.8
1511   - write (*,*) 'Error in CSF values for dead zooplankton : CSF .NE. [0-1]',CSF_dzo
1512   - endif !#5.8
1513   -
1514   - w_dzo_m=(1/(18*(kinvis*densFlu)))*(1/CSF_dzo)*(Rdzo)*g*(diam_dzo)**2 ![m/s]
  1562 + ! w_dzo_m=(1/(18*(kinvis*densFlu)))*(1/CSF_dzo_calc)*(Rdzo)*g*(diam_dzo)**2 ![m/s]
  1563 + w_dzo_m=(1/(18*(kinvis*densFlu)))*(1/CSF_dzo_calc)*(Rdzo/densFlu)*g*(diam_dzo)**2 ![m/s]
1515 1564  
1516 1565 !--Fecal pellets
1517   -
1518   - if(CSF_fp .ge. 0.0 .and. CSF_fp.lt. 0.4) then !#5.9
1519   - CSF_fp = 2.18-(2.09*CSF_fp)
1520   - else if (CSF_fp .ge. 0.4 .and. CSF_fp .lt. 0.8) then !#5.9
1521   - CSF_fp = 0.946*(CSF_fp)**(-0.378)
1522   - else if (CSF_fp .ge. 0.8 .and. CSF_fp .le. 1.0) then !#5.9
1523   - CSF_fp = 1.0 ! consider as a sphere
1524   - else !#5.9
1525   - write (*,*) 'Error in CSF values for fecal pellets : CSF .NE. [0-1]',CSF_fp
1526   - endif !#5.9
1527   -
1528   - w_fp_m=(1/(18*(kinvis*densFlu)))*(1/CSF_fp)*(Rfp)*g*(diam_fp)**2 ![m/s]
  1566 + ! w_fp_m=(1/(18*(kinvis*densFlu)))*(1/CSF_fp_calc)*(Rfp)*g*(diam_fp)**2 ![m/s]
  1567 + w_fp_m=(1/(18*(kinvis*densFlu)))*(1/CSF_fp_calc)*(Rfp/densFlu)*g*(diam_fp)**2 ![m/s]
1529 1568  
1530 1569 !!--------Depending on CSF + Re
1531 1570 else if (Phys_w .eq. 5.0) then !#5
1532 1571  
1533 1572 !--Phytoplankton
1534 1573 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
  1574 + ! w_p_m = -(((Rp/densFlu)*g*(diam_phy**2))/(18*kinvis)) ![m/s]
  1575 + w_p_m=(1/(18*(kinvis*densFlu)))*(1/CSF_phy_calc)*(Rp/densFlu)*g*(diam_phy)**2 ![m/s]
  1576 + !! write(*,*) 'CSF_phy = 1 --> Sphere, then Stokes Law ',w_p_m
1537 1577 ! Comparaison with the Stokes range ! VanRijn 1993
1538 1578 Re_phy = abs(w_p_m*(diam_phy/kinvis))
1539 1579 if (Re_phy .gt. 1.0) then
... ... @@ -1543,51 +1583,55 @@ else if (Phys_w .eq. 5.0) then !#5
1543 1583 else ! No spherical
1544 1584 if(diam_phy .le. 0.0001 )then
1545 1585 ! 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
  1586 + ! w_p_m =-(((Rp/densFlu)*g*(diam_phy**2))/(18*kinvis)) ![m/s]
  1587 + w_p_m=(1/(18*(kinvis*densFlu)))*(1/CSF_phy_calc)*(Rp/densFlu)*g*(diam_phy)**2 ![m/s]
  1588 + !! write(*,*) 'CSF .ne. 1.0 --> Diameter of phy : 1< d <= 100 micrometers, settling =' ,w_p_m
  1589 + elseif(diam_phy .gt. 0.0001 .and. diam_phy .lt. 0.001)then
1549 1590 ! If the diameter of Phy (d) --> 100< d micrometers (VanRijn 1993)
1550 1591 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
  1592 + !! write(*,*) 'CSF .ne. 1.0 --> Diameter of phy : 100< d micrometers, settling =' ,w_p_m
  1593 + elseif(diam_phy .ge. 0.001) then
1555 1594 ! 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]
  1595 + ! w_p_m = (1.1*((Rp/densFlu)*g*diam_phy)**0.5) ![m/s]
  1596 + w_p_m = (0.93*sqrt(Rp*g*diam_phy/densFlu))![m/s]
1557 1597 !! write(*,*) 'Diameter of Phy : d > 1000 micrometers, settling =' ,w_p_m
1558 1598 endif
1559 1599 endif
1560 1600 !--Dead Phytoplankton
  1601 +write(*,*) 'CSF_dph_Calculated',CSF_dph_calc
  1602 +write(*,*) 'diam_dph',diam_dph
1561 1603 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
  1604 + ! w_dph_m = -(((Rdph/densFlu)*g*(diam_dph**2))/(18*kinvis)) ![m/s]
  1605 + w_dph_m=(1/(18*(kinvis*densFlu)))*(1/CSF_dph_calc)*(Rdph/densFlu)*g*(diam_dph)**2 ![m/s]
  1606 + write(*,*) 'CSF_dph = 1 --> Sphere, then Stokes Law ',w_dph_m
1564 1607 ! Comparaison with the Stokes range ! VanRijn 1993
1565 1608 Re_dph = abs(w_dph_m*(diam_dph/kinvis))
1566 1609 if (Re_dph .gt. 1.0) then
1567 1610 w_dph_m = -diam_dph**0.5/secs_pr_day ![m/s]
1568   - !! write(*,*) 'Reynolds number for Dead-phytoplankton > 1',w_dph_m
  1611 + write(*,*) 'Reynolds number for Dead-phytoplankton > 1',w_dph_m
1569 1612 endif
1570 1613 else ! No spherical
1571 1614 if(diam_dph .le. 0.0001 )then
1572 1615 ! 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
  1616 + ! w_dph_m =-(((Rdph/densFlu)*g*(diam_dph**2))/(18*kinvis)) ![m/s]
  1617 + w_dph_m=(1/(18*(kinvis*densFlu)))*(1/CSF_dph_calc)*(Rdph/densFlu)*g*(diam_dph)**2 ![m/s]
  1618 + write(*,*) 'CSF .ne. 1.0 --> Diameter of dph : 1< d <= 100 micrometers, settling =' ,w_dph_m
  1619 + elseif(diam_dph .gt. 0.0001 .and. diam_dph .lt. 0.001 )then
1576 1620 ! If the diameter of dph (d) --> 100< d micrometers (VanRijn 1993)
1577 1621 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
  1622 + write(*,*) 'CSF .ne. 1.0 --> Diameter of dph : 100< d micrometers, settling =' ,w_dph_m
  1623 + elseif(diam_dph .ge. 0.001) then
1582 1624 ! 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
  1625 + ! w_dph_m = (1.1*((Rdph/densFlu)*g*diam_dph)**0.5) ![m/s]
  1626 + w_dph_m = (0.93*sqrt(Rdph*g*diam_dph/densFlu)) ![m/s]
  1627 + write(*,*) 'Diameter of dph : d > 1000 micrometers, settling =' ,w_dph_m
1585 1628 endif
1586 1629 endif
1587 1630 !--Dead Zooplankton
1588 1631 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
  1632 + ! w_dzo_m = -(((Rdzo/densFlu)*g*(diam_dzo**2))/(18*kinvis)) ![m/s]
  1633 + w_dzo_m=(1/(18*(kinvis*densFlu)))*(1/CSF_dzo_calc)*(Rdzo/densFlu)*g*(diam_dzo)**2 ![m/s]
  1634 + !! write(*,*) 'CSF_dzo = 1 --> Sphere, then Stokes Law ',w_dzo_m
1591 1635 ! Comparaison with the Stokes range ! VanRijn 1993
1592 1636 Re_dzo = abs(w_dzo_m*(diam_dzo/kinvis))
1593 1637 if (Re_dzo .gt. 1.0) then
... ... @@ -1597,47 +1641,48 @@ else if (Phys_w .eq. 5.0) then !#5
1597 1641 else ! No spherical
1598 1642 if(diam_dzo .le. 0.0001 )then
1599 1643 ! 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
  1644 + ! w_dzo_m =-(((Rdzo/densFlu)*g*(diam_dzo**2))/(18*kinvis)) ![m/s]
  1645 + w_dzo_m=(1/(18*(kinvis*densFlu)))*(1/CSF_dzo_calc)*(Rdzo/densFlu)*g*(diam_dzo)**2 ![m/s]
  1646 + !! write(*,*) 'CSF .ne.1.0 --> Diameter of dzo: 1< d <= 100 micrometers, settling =' ,w_dzo_m
  1647 + elseif(diam_dzo .gt. 0.0001 .and. diam_dzo .lt. 0.001)then
1603 1648 ! If the diameter of dzo (d) --> 100< d micrometers (VanRijn 1993)
1604 1649 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
  1650 + !! write(*,*) 'CSF .ne. 1.0 --> Diameter of dzo : 100< d micrometers, settling =' ,w_dzo_m
  1651 + elseif(diam_dzo .ge. 0.001) then
1610 1652 ! 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]
  1653 + ! w_dzo_m = (1.1*((Rdzo/densFlu)*g*diam_dzo)**0.5) ![m/s]
  1654 + w_dzo_m = (0.93*sqrt(Rdzo*g*diam_dzo/densFlu)) ![m/s]
1612 1655 !! write(*,*) 'Diameter of Dzo : d > 1000 micrometers, settling =' ,w_dzo_m
1613 1656 endif
1614 1657 endif
1615 1658 !--Fecal pellets
  1659 +write(*,*) 'CSF_fp_Calculated',CSF_fp_calc
  1660 +write(*,*) 'diam_fp',diam_fp
1616 1661 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
  1662 + ! w_fp_m = -(((Rfp/densFlu)*g*(diam_fp**2))/(18*kinvis)) ![m/s]
  1663 + w_fp_m=(1/(18*(kinvis*densFlu)))*(1/CSF_fp_calc)*(Rfp/densFlu)*g*(diam_fp)**2 ![m/s]
  1664 + write(*,*) 'CSF_fp = 1 --> Sphere, then Stokes Law ',w_fp_m
1619 1665 ! Comparaison with the Stokes range ! VanRijn 1993
1620 1666 Re_fp = abs(w_fp_m*(diam_fp/kinvis))
1621 1667 if (Re_fp .gt. 1.0) then
1622 1668 w_fp_m = -diam_fp**0.5/secs_pr_day ![m/s]
1623   - !! write(*,*) 'Reynolds number for Fecal pellets > 1',w_fp_m
  1669 + write(*,*) 'Reynolds number for Fecal pellets > 1',w_fp_m
1624 1670 endif
1625 1671 else ! No spherical
1626 1672 if( diam_fp .le. 0.0001 )then
1627 1673 ! 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
  1674 + ! w_fp_m =-(((Rfp/densFlu)*g*(diam_fp**2))/(18*kinvis)) ![m/s]
  1675 + w_fp_m=(1/(18*(kinvis*densFlu)))*(1/CSF_fp_calc)*(Rfp/densFlu)*g*(diam_fp)**2 ![m/s]
  1676 + write(*,*) 'CSF .ne. 1.0 --> Diameter of fp: 1< d <= 100 micrometers, settling =' ,w_fp_m
  1677 + elseif(diam_fp .gt. 0.0001 .and. diam_fp .lt. 0.001)then
1631 1678 ! If the diameter of fp (d) --> 100< d micrometers (VanRijn 1993)
1632 1679 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
  1680 + write(*,*) 'CSF .ne. 1.0 --> Diameter of fp : 100< d micrometers, settling =' ,w_fp_m
  1681 + elseif(diam_fp .ge. 0.001) then
1637 1682 ! 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
  1683 + ! w_fp_m = (1.1*((Rfp/densFlu)*g*diam_fp)**0.5) ![m/s]
  1684 + w_fp_m = (0.93*sqrt(Rdzo*g*diam_dzo/densFlu)) ![m/s]
  1685 + write(*,*) 'Diameter of Fp : d > 1000 micrometers, settling =' ,w_fp_m
1641 1686 endif
1642 1687 endif
1643 1688 else !#5
... ... @@ -1676,7 +1721,6 @@ write(*,*) &#39;Reynolds number for fecal pellets &gt; 1&#39;
1676 1721 endif !#9
1677 1722  
1678 1723 endif
1679   -
1680 1724  
1681 1725 !Everyone will settle at the end of the loop (see below)
1682 1726  
... ... @@ -2155,6 +2199,7 @@ else !#2.1
2155 2199 endif !#2.1
2156 2200 !! ---Determination of physic_Br
2157 2201 physic_Br = (2*kB*T_degK)/(3*dynvis)
  2202 +
2158 2203  
2159 2204 !-------------------------------------------------------------------------------------
2160 2205 !-------------------------------------------------------------------------------------
... ... @@ -2320,14 +2365,16 @@ endif !#1
2320 2365 !-----------------------------------------------------------------------------------
2321 2366 ! SIZE Marine Snow with aggregation
2322 2367 !-----------------------------------------------------------------------------------
2323   -
2324   - ! ---Maximum size at dissagregation (Alldredge 1990)
2325   - if(dm_msn .eq. 1.0) then
2326   - diam_msn_max = coef1*(eps(ci))**(-coef2) ! [m]
2327   - write(*,*)'At depth-----',ci
2328   - write(*,*)'Maximum size diameter of Msn from eps :', diam_msn_max
2329   - else
  2368 +
  2369 + ! ---Maximum size at dissagregation
  2370 + if(dm_msn .eq. 1.0) then !--» Comparaison with CFL limit
  2371 + diam_msn_max = ((CFL/0.93)**2)*(1/(g*(Rmsn/densFlu)))
  2372 + else !--» From the nml
2330 2373 diam_msn_max = diam_msn_us ! [m]
  2374 + if (diam_msn_us .gt. ((CFL/0.93)**2)*(1/(g*(Rmsn/densFlu))) ) then
  2375 + write(*,*)'The Max. Diam mentionned in the nml do not allow to respect CFL limit',ci
  2376 + stop
  2377 + endif
2331 2378 endif
2332 2379  
2333 2380 ! ---Flux from particules to msn [mmol/m3/s]
... ... @@ -2604,19 +2651,17 @@ endif
2604 2651 ! We want to use a calculated settling velocity for the msn
2605 2652 ! Settling velocity send to the advection scheme as well as the size of the msn
2606 2653 !-------------------------------------------------------------------------------------
2607   -write(*,*)'Depth',ci
2608   -write(*,*)'[msn]',cc(d4,ci)
2609   -write(*,*) 'Diameter of msn :',cc(size_msn,ci)
2610   - !--» cc(size_msn,ci) --> is the diameter after coag and frag
2611   - Rmsn = (rho_msn-densFlu)
  2654 +
  2655 +!--» cc(size_msn,ci) --> is the diameter after coag and frag
  2656 +
2612 2657  
2613   -!if (cc(d4,ci) .le. cons_min) then
2614   -! w_msn_m = 0.0
2615   - ! write(*,*)'[msn] lt Cons_min, settling velocity set to :',w_msn_m
  2658 +if (cc(d4,ci) .le. cons_min) then
  2659 + w_msn_m = 0.0
  2660 + write(*,*)'[msn] lt Cons_min, settling velocity set to :',w_msn_m
2616 2661  
2617 2662 !--» We set here the concentration limit in the water column from free settling of marine snow vs flocculation settling (Metha 1989)
2618   -!elseif (cc(d4,ci) .lt. cons_max) then !#1
2619   - if (cc(d4,ci) .lt. cons_max) then !#1
  2663 +elseif (cc(d4,ci) .lt. cons_max) then !#1
  2664 + ! if (cc(d4,ci) .lt. cons_max) then !#1
2620 2665 !! --- Calculated by value from the .nml
2621 2666 if (w_msnow .eq. 0.0 ) then !#2
2622 2667 w_msn_m = w_msn ![m/s]
... ... @@ -2647,108 +2692,47 @@ write(*,*) &#39;Diameter of msn :&#39;,cc(size_msn,ci)
2647 2692 else if (w_msnow .eq. 4.0) then !#2
2648 2693  
2649 2694 if(CSF .ge. 0.0 .and. CSF.lt. 0.4) then !#2.3
2650   - CSF_msn = 2.18-(2.09*CSF)
  2695 + CSF_msn_calc = 2.18-(2.09*CSF)
2651 2696 else if (CSF .ge. 0.4 .and. CSF .lt. 0.8) then!#2.3
2652   - CSF_msn = 0.946*(CSF)**(-0.378)
  2697 + CSF_msn_calc = 0.946*(CSF)**(-0.378)
2653 2698 else if (CSF .ge. 0.8 .and. CSF .le. 1.0) then!#2.3
2654   - CSF_msn = 1.0 ! consider as a sphere
  2699 + CSF_msn_calc = 1.0 ! consider as a sphere
2655 2700 else !#2.3
2656 2701 write (*,*) 'Error in CSF values for Marine snow : CSF .NE. [0-1]'
2657 2702 endif !#2.3
2658 2703  
2659   - w_msn_m=(1/(18*(kinvis*densFlu)))*(1/CSF_msn)*(Rmsn)*g*(cc(size_msn,ci))**2 ![m/s]
2660   -
2661   - elseif (w_msnow .eq. 5.0) then !#2!
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
2666   - ! If the diameter of Msn (d) --> 1< d <= 100 micrometers (VanRijn 1993)
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
2669   -
2670   -
2671   - elseif(cc(size_msn,ci) .gt. 0.0001 .and. cc(size_msn,ci) .lt. 0.001 ) then
2672   - ! If the diameter of Msn (d) --> 100< d <= 1000 micrometers (VanRijn 1993)
2673   - w_msn_m = -(((10*kinvis)/cc(size_msn,ci))*&
2674   - ((1+((0.01*(Rmsn/densFlu)*g*cc(size_msn,ci)**3.0) /kinvis**2.0)**0.5) -1)) ![m/s]
2675   - write(*,*) 'Diameter of msn : 100< d <= 1000 micrometers, settling =' ,w_msn_m
2676   -
2677   - else !(cc(size_msn,ci).gt. 0.001) then
2678   - ! If the diameter of Msn (d) --> d > 1000 micrometers (VanRijn 1993)
2679   - w_msn_m = -(1.1*((Rmsn/densFlu)*g*cc(size_msn,ci))**0.5) ![m/s]
2680   - write(*,*) 'Diameter of msn : d > 1000 micrometers, settling =' ,w_msn_m
2681   -
2682   - endif
2683   -
2684   -
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
2707   -
2708   -
2709   - elseif (w_msnow .eq. 7.0) then !#2!
  2704 + w_msn_m=(1/(18*(kinvis*densFlu)))*(1/CSF_msn_calc)*(Rmsn)*g*(cc(size_msn,ci))**2 ![m/s]
  2705 +
  2706 +
  2707 + !! ---[Komar 1978] and VanRijn 1993
  2708 + elseif (w_msnow .eq. 5.0) then !#2!
2710 2709  
2711 2710 ! If the diameter of Msn (d) --> d <= 100 micrometers (Komar 1978)
2712   - if(cc(size_msn,ci) .le. 0.0001 )then
2713   -
  2711 + if(cc(size_msn,ci) .le. 0.0001 )then
2714 2712 if(CSF .ge. 0.0 .and. CSF.lt. 0.4) then !#2.3
2715   - CSF_msn = 2.18-(2.09*CSF)
  2713 + CSF_msn_calc = 2.18-(2.09*CSF)
2716 2714 else if (CSF .ge. 0.4 .and. CSF .lt. 0.8) then!#2.3
2717   - CSF_msn = 0.946*(CSF)**(-0.378)
  2715 + CSF_msn_calc = 0.946*(CSF)**(-0.378)
2718 2716 else if (CSF .ge. 0.8 .and. CSF .le. 1.0) then!#2.3
2719   - CSF_msn = 1.0 ! consider as a sphere
  2717 + CSF_msn_calc = 1.0 ! consider as a sphere
2720 2718 else !#2.3
2721 2719 write (*,*) 'Error in CSF values for Marine snow : CSF .NE. [0-1]'
2722 2720 endif !#2.3
2723   - w_msn_m=(1/(18*(kinvis*densFlu)))*(1/CSF_msn)*(Rmsn/densFlu)*g*(cc(size_msn,ci))**2 ![m/s]
  2721 + w_msn_m=(1/(18*(kinvis*densFlu)))*(1/CSF_msn_calc)*(Rmsn/densFlu)*g*(cc(size_msn,ci))**2 ![m/s]
2724 2722 write(*,*) 'Diameter of msn : d <= 100 micrometers, settling (Stokes Law is used !) =' ,w_msn_m
2725 2723  
2726 2724 ! 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
  2725 + elseif(cc(size_msn,ci) .gt. 0.0001 .and. cc(size_msn,ci) .lt. 0.001 ) then
2729 2726 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
2751   -
  2727 + write(*,*) 'Diameter of msn : 100< d <= 1000 micrometers, settling =' ,w_msn_m
  2728 +
  2729 + elseif( cc(size_msn,ci) .ge. 0.001 ) then
  2730 + w_msn_m = (0.93*sqrt(Rmsn*g*cc(size_msn,ci)/densFlu)) ![m/s]
  2731 + write(*,*) 'Diameter of msn : d > 1000 micrometers, settling =' ,w_msn_m
  2732 + else
  2733 + write(*,*) ' No size mentioned for Msn',cc(size_msn,ci)
  2734 + stop
  2735 + endif
2752 2736  
2753 2737 else !#2!
2754 2738 !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
... ... @@ -2812,7 +2796,7 @@ endif !#1
2812 2796 endif
2813 2797  
2814 2798 if (rho_msn .gt. densFlu .AND. w_msn_m .gt. 0.0) then
2815   - write(*,*)' Msn density > dens_fluid + W_msn_m > 0 --> Settling go up, but is constrain to go down !'
  2799 + write(*,*)' Msn density > dens_fluid AND W_msn_m > 0 --> Settling go up, but is constrain to go down !'
2816 2800 w_msn_m = -w_msn_m
2817 2801 endif
2818 2802  
... ... @@ -2823,38 +2807,40 @@ endif !#1
2823 2807 !! ---Safety check
2824 2808 !-----------------------------------------------------------------------------------
2825 2809  
2826   - !--» Does CFL limit respected ?
2827   - CFL = (depth_bio/nlev)/dt_bio
  2810 +!--» Does CFL limit respected ?
2828 2811  
2829 2812 !--» Case of Phytoplankton
2830 2813 if (abs(w_p_m) .gt. (depth_bio/nlev)/dt_bio )then
2831 2814 write(*,*) 'CFL Constrain is OVERPASS for W-Phy'
2832   -stop
  2815 +w_p_m = -((depth_bio/nlev)/dt_bio)
2833 2816 endif
2834 2817 ws(p,ci) = w_p_m
2835 2818  
2836 2819 !--» Case of Dead_Phytoplankton
2837 2820 if (abs(w_dph_m) .gt. (depth_bio/nlev)/dt_bio )then
2838 2821 write(*,*) 'CFL Constrain is OVERPASS for W-Dph'
2839   -stop
  2822 +w_dph_m = -((depth_bio/nlev)/dt_bio)
2840 2823 endif
2841 2824 ws(d1,ci) = w_dph_m
  2825 +
2842 2826 !--» Case of Dead Zooplankton
2843 2827 if (abs(w_dzo_m) .gt. (depth_bio/nlev)/dt_bio )then
2844 2828 write(*,*) 'CFL Constrain is OVERPASS for W-dzo'
2845   -stop
  2829 +w_dzo_m = -((depth_bio/nlev)/dt_bio)
2846 2830 endif
2847 2831 ws(d2,ci) = w_dzo_m
  2832 +
2848 2833 !--» Case of Fecal pellets
2849 2834 if (abs(w_fp_m) .gt. (depth_bio/nlev)/dt_bio )then
2850 2835 write(*,*) 'CFL Constrain is OVERPASS for W-fp'
2851   -stop
  2836 +w_fp_m = -((depth_bio/nlev)/dt_bio)
2852 2837 endif
2853   -ws(d3,ci) = w_fp_m
  2838 +ws(d3,ci) = w_fp_m
  2839 +
2854 2840 !--» Case of marine snow
2855 2841 if (abs(w_msn_m) .gt. (depth_bio/nlev)/dt_bio )then
2856 2842 write(*,*) 'CFL Constrain is OVERPASS for W-msn'
2857   -!stop
  2843 +w_msn_m = -((depth_bio/nlev)/dt_bio)
2858 2844 endif
2859 2845 ws(d4,ci) = w_msn_m
2860 2846  
... ...
src/util/adv_center.F90
  Diff comments   View file @764531e
... ... @@ -377,6 +377,7 @@ if (mode.eq.3) then
377 377 ! compute limited flux
378 378 cu(k) =ww(k)*(Yc+0.5*limit*(1-c)*(Yd-Yc))
379 379 cu(k-1)=ww(k)*(Yd+0.5*limit*(1-c)*(Yd-Y(k-1)))
  380 +
380 381  
381 382 else ! For the other variables of bio_polynow, or other biological models
382 383 ! compute the limited flux
... ...