Commit c946efdb0003e6deaf62241c8317c292d0e79c1a

Authored by Gwenaelle Gremion
1 parent e48ad91f
Exists in master and in 1 other branch snow

Set up Wd4 Max

Showing 1 changed file with 50 additions and 38 deletions   Show diff stats
src/extras/bio/bio_polynow.F90
... ... @@ -605,7 +605,7 @@ leak=leak/secs_pr_day
605 605 w_dph = w_dph /secs_pr_day
606 606 w_dzo = w_dzo /secs_pr_day
607 607 w_fp = w_fp /secs_pr_day
608   -! w_msn = w_msn /secs_pr_day
  608 + w_msn = w_msn /secs_pr_day
609 609  
610 610 litt_msn_w =litt_msn_w /secs_pr_day
611 611 !Phytoplankton
... ... @@ -1164,7 +1164,7 @@ end if
1164 1164 REALTYPE ::prob_break_msn,r_msn
1165 1165 REALTYPE :: diam_msn_max,RFV_msn,min_size_msn
1166 1166 ! Sedimentation rate msn
1167   - REALTYPE :: w_msn_m,Re_msn,Rmsn,w_min_m,CFL
  1167 + REALTYPE :: w_msn_m,Re_msn,Rmsn,w_min_m,CFL,max_w_msn
1168 1168  
1169 1169 !EOP
1170 1170 !-----------------------------------------------------------------------
... ... @@ -2291,7 +2291,7 @@ if (Coag_coef .eq. 1.0) then ! #1 Size modified only if coagulation happenned
2291 2291 !--» At this points 'taille_coag' is the projected area equivalent diameter in m.
2292 2292  
2293 2293 ! ---Elaboration of the size variable (diameter), available for the fragmentation eventually (taille_intrm).
2294   - if(cc(d4,ci) .gt. cons_min) then
  2294 + if(cc(d4,ci) .gt. cons_min) then
2295 2295 cc(taille_intrm,ci) = cc(taille_coag,ci) + cc(taille_msn,ci) ![m]
2296 2296 !--» Msn size after coagulation will depend on the increase of size due to coagaluation (taille_coag) &
2297 2297 !--» + the size at the previous time step (taille_msn) of the particles
... ... @@ -2307,11 +2307,11 @@ if (Coag_coef .eq. 1.0) then ! #1 Size modified only if coagulation happenned
2307 2307 write(*,*)'Taille _msn - May be equal to 0 :', cc(taille_msn,ci)
2308 2308 endif
2309 2309  
2310   - else
  2310 + else
2311 2311 !--» If there is not need to have a size augmentation due to low concentration of marine snow, the size of previous step remain
2312 2312 cc(taille_intrm,ci) = cc(taille_msn,ci) ![m]
2313   - write(*,*)'No size augmentation as [msn] < cons_min - The size stays the same', cc(taille_msn,ci)
2314   - endif
  2313 + ! write(*,*)'No size augmentation as [msn] < cons_min - The size stays the same', cc(taille_msn,ci)
  2314 + endif
2315 2315  
2316 2316 else !#1 Size not modified because no coagulation happenned
2317 2317 cc(taille_coag,ci) = 0.0
... ... @@ -2334,8 +2334,8 @@ endif !#1
2334 2334 if (cc(d4,ci).lt. cons_min)then
2335 2335 cc(taille_msn,ci) = _ZERO_
2336 2336 cc(taille_intrm,ci) = _ZERO_
2337   - write(*,*)'Size but no msn at :',ci
2338   - endif
  2337 + ! write(*,*)'Size but no msn at :',ci
  2338 + endif
2339 2339  
2340 2340 !-----------------------------------------------------------------------------------
2341 2341 ! FRAGMENTATION
... ... @@ -2511,15 +2511,16 @@ endif
2511 2511  
2512 2512 if (cc(d4,ci) .le. cons_min) then
2513 2513 w_msn_m = 0.0
2514   - write(*,*)'[msn] lt Cons_min, settling velocity set to :',w_msn_m
  2514 + ! write(*,*)'[msn] lt Cons_min, settling velocity set to :',w_msn_m
2515 2515  
2516 2516 elseif (cc(d4,ci) .lt. cons_max) then !#1
  2517 +
2517 2518 !--» We set here the concentration limit in the water column from free settling of marine snow vs flocculation settling (Metha 1989)
2518 2519  
2519 2520 !! --- Calculated by value from the .nml
2520 2521 if (w_msnow .eq. 0.0 ) then !#2
2521 2522 !w_msn_m = w_msn*sqrt(w_depth/100)
2522   - w_msn_m = w_msn_m
  2523 + w_msn_m = w_msn
2523 2524  
2524 2525 !! --- Stokes's law from Ghosh et al 2013
2525 2526 !--» msn is considered as sphere
... ... @@ -2564,7 +2565,7 @@ elseif (cc(d4,ci) .lt. cons_max) then !#1
2564 2565 endif!#2
2565 2566 elseif (cc(d4,ci) .ge. cons_max)then !#1
2566 2567 w_msn_m=-(coef5*cc(d4,ci)**(1.6)) ! Metha 1989
2567   - write(*,*)'[msn] > cons_max at :',ci
  2568 +! write(*,*)'[msn] > cons_max at :',ci
2568 2569  
2569 2570 else
2570 2571 write(*,*)'At Depth : ',ci
... ... @@ -2638,6 +2639,27 @@ endif !#3
2638 2639 w_msn_m = -w_msn_m
2639 2640 endif
2640 2641  
  2642 + !--» Definition of the maximum diameter size from our model calculation
  2643 +
  2644 +
  2645 +If (litt_msn_w .lt. max_w_msn) then
  2646 +max_w_msn = litt_msn_w
  2647 +else
  2648 +max_w_msn = max_w_msn
  2649 +endif
  2650 +
  2651 + !--» Does CFL limit respected ?
  2652 + CFL = (depth_bio/nlev)/dt_bio
  2653 +
  2654 +if(abs(w_msn_m) .gt. CFL) then
  2655 +max_w_msn = max_w_msn
  2656 +else
  2657 + if(w_msn_m .lt. max_w_msn) then
  2658 + max_w_msn = w_msn_m
  2659 + else
  2660 + max_w_msn = max_w_msn
  2661 + endif
  2662 +endif
2641 2663  
2642 2664 !--» let them settling(m/s)
2643 2665 !--» The ncdf output(bio_save) do ws(x,ci) *86400 to have the ouptut in m/j.
... ... @@ -2645,8 +2667,8 @@ endif !#3
2645 2667 !-----------------------------------------------------------------------------------
2646 2668 !! ---Safety check
2647 2669 !-----------------------------------------------------------------------------------
  2670 +
2648 2671 !--» Does CFL limit respected ?
2649   - CFL = (depth_bio/nlev)/dt_bio
2650 2672  
2651 2673 !--» Case of Phytoplankton
2652 2674 if (abs(w_p_m) .gt. (depth_bio/nlev)/dt_bio )then
... ... @@ -2700,34 +2722,16 @@ endif
2700 2722 !--Concentration - Vs CFL limit
2701 2723 if (abs(w_msn_m) .gt. (depth_bio/nlev)/dt_bio )then
2702 2724  
2703   - if (w_msnow .eq. 4.0) then !#2
2704   -
2705   - if(CSF .ge. 0.0 .and. CSF.lt. 0.4) then !#2.3
2706   - CSF_msn = 2.18-(2.09*CSF)
2707   - else if (CSF .ge. 0.4 .and. CSF .lt. 0.8) then!#2.3
2708   - CSF_msn = 0.946*(CSF)**(-0.378)
2709   - else if (CSF .ge. 0.8 .and. CSF .le. 1.0) then!#2.3
2710   - CSF_msn = 1.0 ! consider as a sphere
2711   - else !#2.3
2712   - write (*,*) 'Error in CSF values for Marine snow : CSF .NE. [0-1]'
2713   - endif !#2.3
2714   -
2715   - w_msn_m=(1/(18*(kinvis*densFlu)))*(1/CSF_msn)*(Rmsn)*g*(diam_msn_max)**2 ![m/s]
2716   - if (w_msn_m .gt. 0.0) then
2717   - ws(d4,ci) = -w_msn_m
2718   - else
2719   - ws(d4,ci) = w_msn_m
2720   - endif
2721   - else
2722 2725 ! -- As settling velocity overpass the CFL constrain we will attribute to the msn
2723   -! the maximum settling velocity found in the litterature
2724   - ws(d4,ci) = litt_msn_w ![m/s]
  2726 +! the maximum settling velocity found between our simulation and the litterature
  2727 + ! ws(d4,ci) = litt_msn_w ![m/s]
  2728 + ws(d4,ci) = max_w_msn ![m/s]
2725 2729 write(*,*) 'CFL Constrain is OVERPASS, w_msn use from nml',ws(d4,ci)
2726 2730 if (abs(litt_msn_w) .gt. (depth_bio/nlev)/dt_bio )then
2727 2731 ws(d4,ci) = 0.0
2728 2732 write(*,*) 'CFL Constrain is OVERPASS for w_msn,from the nml value', ws(d4,ci)
2729 2733 endif
2730   - endif
  2734 +
2731 2735 else
2732 2736 ws(d4,ci) = w_msn_m
2733 2737 write(*,*) 'CFL limit is respected, ws(D4,ci) = ',ws(d4,ci)
... ... @@ -2737,11 +2741,10 @@ endif
2737 2741 ws(taille_msn,ci) = ws(d4,ci)
2738 2742  
2739 2743 !--Settling velocity trait - Vs CFL limit
  2744 +
2740 2745 cc(settl_msn,ci) = ws(d4,ci)
2741 2746 ws(settl_msn,ci) = ws(d4,ci)
2742 2747  
2743   -
2744   -
2745 2748  
2746 2749 !-----------------------------------------------------------------------------------
2747 2750 ! NET PRIMARY PRODUCTION
... ... @@ -3052,8 +3055,17 @@ endif !#A
3052 3055 do j=1,numc
3053 3056 pp(i,j,ci)=dd(j,i,ci)
3054 3057 end do
3055   - end do
3056   - end do
  3058 + end do
  3059 + end do
  3060 +
  3061 +!-----------------------------------------------------------------------
  3062 +! If you want to do something over this timestep do it here !
  3063 +!-----------------------------------------------------------------------
  3064 +
  3065 +
  3066 +write(*,*)'-------------------------NEW TIMESTEP-----------------------'
  3067 +
  3068 +!-----------------------------------------------------------------------
3057 3069  
3058 3070 return
3059 3071 end subroutine do_bio_polynow
... ...