Commit e48ad91f1409b7c5f24db9716f6bd4eb5426b552

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

Clarification de la contrainte CFL

Showing 2 changed files with 175 additions and 80 deletions   Show diff stats
src/extras/bio/bio_polynow.F90
... ... @@ -213,7 +213,8 @@
213 213 REALTYPE :: size_fp_up= 0.000050
214 214 REALTYPE :: size_msn = -0.5
215 215 REALTYPE :: dm_msn = 0.0
216   - REALTYPE :: diam_msn_us = 0.01
  216 + REALTYPE :: diam_msn_us = 0.01
  217 + REALTYPE :: litt_msn_w = 0.01
217 218 !Parametres pour l'aggregation physique
218 219 REALTYPE :: dynvis = -0.5
219 220 REALTYPE :: kinvis = -0.5
... ... @@ -307,7 +308,7 @@ rho_p,rho_dph,rho_dzo,rho_fp,rho_msn,
307 308 sti_cst, &
308 309 stip_p,stip_dph,stip_dzo,stip_fp,stidph_dph,stidph_dzo,stidph_fp,stidzo_dzo,stidzo_fp,stifp_fp, &
309 310 CSF,size_rand,size_phy_us,size_phy_up,size_dph_us,size_dph_up,size_dzo_us,size_dzo_up,&
310   - size_fp_us,size_fp_up,size_msn,dm_msn,diam_msn_us,&
  311 + size_fp_us,size_fp_up,size_msn,dm_msn,diam_msn_us,litt_msn_w,&
311 312 dynvis,kinvis,kB, &
312 313 Frag_meth,swim_brk, &
313 314 Floc_coef, &
... ... @@ -549,7 +550,8 @@ Floc_coef, &
549 550 write(*,900) ' size_msn = ',size_msn
550 551 write(10,901) size_msn
551 552 !dm_msn
552   -!diam_msn_us
  553 +!diam_msn_us
  554 +!litt_msn_w
553 555 write(*,900) ' dynvis = ',dynvis
554 556 write(10,901) dynvis
555 557 write(*,900) ' kinvis = ',kinvis
... ... @@ -605,7 +607,7 @@ leak=leak/secs_pr_day
605 607 w_fp = w_fp /secs_pr_day
606 608 ! w_msn = w_msn /secs_pr_day
607 609  
608   -
  610 +litt_msn_w =litt_msn_w /secs_pr_day
609 611 !Phytoplankton
610 612 txloss_p = txloss_p /secs_pr_day
611 613 betap_p = betap_p /secs_pr_day ! GG-D
... ... @@ -1162,7 +1164,7 @@ end if
1162 1164 REALTYPE ::prob_break_msn,r_msn
1163 1165 REALTYPE :: diam_msn_max,RFV_msn,min_size_msn
1164 1166 ! Sedimentation rate msn
1165   - REALTYPE :: w_msn_m,Re_msn,Rmsn,w_min_m
  1167 + REALTYPE :: w_msn_m,Re_msn,Rmsn,w_min_m,CFL
1166 1168  
1167 1169 !EOP
1168 1170 !-----------------------------------------------------------------------
... ... @@ -2507,8 +2509,6 @@ endif
2507 2509 !--» Taille_msn --> is the diameter after coag and frag
2508 2510 Rmsn = (rho_msn-densFlu)
2509 2511  
2510   -write(*,*)'[msn]', cc(d4,ci)
2511   -
2512 2512 if (cc(d4,ci) .le. cons_min) then
2513 2513 w_msn_m = 0.0
2514 2514 write(*,*)'[msn] lt Cons_min, settling velocity set to :',w_msn_m
... ... @@ -2577,17 +2577,35 @@ endif !#1
2577 2577 !! --- Convertion & Comparaison with the Stokes range
2578 2578 !-----------------------------------------------------------------------------------
2579 2579  
2580   - !--» Now we put everything in /d not matter the scenario.
2581   - ! w_msn_m = w_msn_m/secs_pr_day
2582   -
2583 2580 !--» Comparaison with the Stokes range ! VanRijn 1993
2584 2581 Re_msn = abs(w_msn_m*(cc(taille_msn,ci)/kinvis))
2585 2582  
2586 2583 if (Re_msn .lt. 1.0)then !#3
2587 2584 w_msn_m = w_msn_m
2588 2585 else !#3
2589   - w_msn_m = (cc(taille_msn,ci)**0.5)
2590   - write(*,*) 'Reynolds number for marine snow > 1'
  2586 + write(*,*) 'Reynolds number for marine snow > 1'
  2587 +
  2588 + if(cc(taille_msn,ci) .gt. 0.000001 .and. cc(taille_msn,ci) .le. 0.0001 )then
  2589 +! If the diameter of Msn (d) --> 1< d <= 100 micrometers (VanRijn 1993)
  2590 + w_msn_m = ((2.65-1)*g*(cc(taille_msn,ci)**2))/18*kinvis
  2591 + write(*,*) 'Diameter of msn : 1< d <= 100 micrometers, settling =' ,w_msn_m
  2592 +
  2593 + elseif(cc(taille_msn,ci) .gt. 0.0001 .and. cc(taille_msn,ci) .lt. 0.001 ) then
  2594 +! If the diameter of Msn (d) --> 100< d <= 1000 micrometers (VanRijn 1993)
  2595 + w_msn_m = ((10*kinvis)/cc(taille_msn,ci))*&
  2596 + ((1+((0.01*(2.65-1)*g*(cc(taille_msn,ci))**3.0) /kinvis**2.0)**0.5) -1)
  2597 + write(*,*) 'Diameter of msn : 100< d <= 1000 micrometers, settling =' ,w_msn_m
  2598 +
  2599 + elseif(cc(taille_msn,ci) .gt. 0.001) then
  2600 +! If the diameter of Msn (d) --> d > 1000 micrometers (VanRijn 1993)
  2601 + w_msn_m = 1.1*((2.65-1)*g*cc(taille_msn,ci))**0.5
  2602 + write(*,*) 'Diameter of msn : d > 1000 micrometers, settling =' ,w_msn_m
  2603 +
  2604 + else
  2605 + write(*,*) 'Diameter of msn lower than 1 micrometers, and RE > 1 : Houston we have a problem ! '
  2606 + w_msn_m = 0.0
  2607 + endif
  2608 +
2591 2609 endif !#3
2592 2610  
2593 2611  
... ... @@ -2595,47 +2613,136 @@ endif !#3
2595 2613 ! ADVECTION (WS) OF THE PARTICLES
2596 2614 !-----------------------------------------------------------------------------------
2597 2615  
  2616 + !-----------------------------------------------------------------------------------
  2617 + !! ---Safety check
  2618 + !-----------------------------------------------------------------------------------
  2619 + !--» Does density difference with the fluid allow to go up ?
  2620 +
  2621 + if (rho_p .gt. densFlu .AND. w_p_m .gt. 0.0) then
  2622 + w_p_m = -w_p_m
  2623 + endif
  2624 +
  2625 + if (rho_dph .gt. densFlu .AND. w_dph_m .gt. 0.0) then
  2626 + w_dph_m = -w_dph_m
  2627 + endif
  2628 +
  2629 + if (rho_dzo .gt. densFlu .AND. w_dzo_m .gt. 0.0) then
  2630 + w_dzo_m = -w_dzo_m
  2631 + endif
  2632 +
  2633 + if (rho_fp .gt. densFlu .AND. w_fp_m .gt. 0.0) then
  2634 + w_fp_m = -w_fp_m
  2635 + endif
  2636 +
  2637 + if (rho_msn .gt. densFlu .AND. w_msn_m .gt. 0.0) then
  2638 + w_msn_m = -w_msn_m
  2639 + endif
  2640 +
  2641 +
2598 2642 !--» let them settling(m/s)
2599   - !--» The ncdf output do ws(x,ci) *86400 to have the ouptut in m/j.
  2643 + !--» The ncdf output(bio_save) do ws(x,ci) *86400 to have the ouptut in m/j.
2600 2644  
2601   -if (rho_p .gt. densFlu .AND. w_p_m .gt. 0.0) then
2602   -w_p_m = -w_p_m
2603   -endif
  2645 + !-----------------------------------------------------------------------------------
  2646 + !! ---Safety check
  2647 + !-----------------------------------------------------------------------------------
  2648 + !--» Does CFL limit respected ?
  2649 + CFL = (depth_bio/nlev)/dt_bio
2604 2650  
2605   -if (rho_dph .gt. densFlu .AND. w_dph_m .gt. 0.0) then
2606   -w_dph_m = -w_dph_m
  2651 + !--» Case of Phytoplankton
  2652 +if (abs(w_p_m) .gt. (depth_bio/nlev)/dt_bio )then
  2653 + ws(p,ci) = w_p
  2654 + write(*,*) 'CFL Constrain is OVERPASS for W-Phy, ws(P,ci) = nml value', ws(p,ci)
  2655 + if (abs(w_p) .gt. (depth_bio/nlev)/dt_bio )then
  2656 + ws(p,ci) = 0.0
  2657 + write(*,*) 'CFL Constrain is OVERPASS for W-Phy,from the nml value', ws(p,ci)
  2658 + endif
  2659 +else
  2660 + ws(p,ci) = w_p_m
2607 2661 endif
2608 2662  
2609   -if (rho_dzo .gt. densFlu .AND. w_dzo_m .gt. 0.0) then
2610   -w_dzo_m = -w_dzo_m
  2663 + !--» Case of Dead_Phytoplankton
  2664 +if (abs(w_dph_m) .gt. (depth_bio/nlev)/dt_bio )then
  2665 + ws(d1,ci) = w_dph
  2666 + write(*,*) 'CFL Constrain is OVERPASS for W-Dph, ws(D1,ci) = nml value', ws(d1,ci)
  2667 + if (abs(w_dph) .gt. (depth_bio/nlev)/dt_bio )then
  2668 + ws(d1,ci) = 0.0
  2669 + write(*,*) 'CFL Constrain is OVERPASS for W-Dph,from the nml value', ws(d1,ci)
  2670 + endif
  2671 +else
  2672 + ws(d1,ci) = w_dph_m
2611 2673 endif
2612 2674  
2613   -if (rho_fp .gt. densFlu .AND. w_fp_m .gt. 0.0) then
2614   -w_fp_m = -w_fp_m
  2675 + !--» Case of Dead Zooplankton
  2676 +if (abs(w_dzo_m) .gt. (depth_bio/nlev)/dt_bio )then
  2677 + ws(d2,ci) = w_dzo
  2678 + write(*,*) 'CFL Constrain is OVERPASS for W-dzo, ws(d2,ci) = nml value', ws(d2,ci)
  2679 + if (abs(w_dzo) .gt. (depth_bio/nlev)/dt_bio )then
  2680 + ws(d2,ci) = 0.0
  2681 + write(*,*) 'CFL Constrain is OVERPASS for W-dzo,from the nml value', ws(d2,ci)
  2682 + endif
  2683 +else
  2684 + ws(d2,ci) = w_dzo_m
2615 2685 endif
2616   -
2617   - ws(p,ci) = w_p_m
2618   - ws(d1,ci) = w_dph_m
2619   - ws(d2,ci) = w_dzo_m
2620   - ws(d3,ci) = w_fp_m
  2686 +
  2687 + !--» Case of Fecal pellets
  2688 +if (abs(w_fp_m) .gt. (depth_bio/nlev)/dt_bio )then
  2689 + ws(d3,ci) = w_fp
  2690 + write(*,*) 'CFL Constrain is OVERPASS for W-fp, ws(d3,ci) = nml value', ws(d3,ci)
  2691 + if (abs(w_fp) .gt. (depth_bio/nlev)/dt_bio )then
  2692 + ws(d3,ci) = 0.0
  2693 + write(*,*) 'CFL Constrain is OVERPASS for W-fp,from the nml value', ws(d3,ci)
  2694 + endif
  2695 +else
  2696 + ws(d3,ci) = w_fp_m
  2697 +endif
2621 2698  
2622 2699 !--» Case of marine snow
  2700 + !--Concentration - Vs CFL limit
  2701 +if (abs(w_msn_m) .gt. (depth_bio/nlev)/dt_bio )then
2623 2702  
2624   - !--Concentration - Vs CFL limit
2625   - ws(d4,ci)= min(w_msn_m,(depth_bio/nlev)/dt_bio)
2626   - if (ws(d4,ci) .eq. w_msn_m) then
2627   - write(*,*) 'CFL limit is respected'
2628   - else
2629   - write(*,*) 'CFL Constrain is OVERPASS'
  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 +! -- 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]
  2725 + write(*,*) 'CFL Constrain is OVERPASS, w_msn use from nml',ws(d4,ci)
  2726 + if (abs(litt_msn_w) .gt. (depth_bio/nlev)/dt_bio )then
  2727 + ws(d4,ci) = 0.0
  2728 + write(*,*) 'CFL Constrain is OVERPASS for w_msn,from the nml value', ws(d4,ci)
2630 2729 endif
2631   - !--Size - Vs CFL limit
2632   - ws(taille_msn,ci) = min(w_msn_m,(depth_bio/nlev)/dt_bio)
2633   -
2634   - !--Settling velocity trait - Vs CFL limit
2635   - cc(settl_msn,ci)= w_msn_m
2636   - ws(settl_msn,ci) = min(w_msn_m,(depth_bio/nlev)/dt_bio)
  2730 + endif
  2731 +else
  2732 +ws(d4,ci) = w_msn_m
  2733 + write(*,*) 'CFL limit is respected, ws(D4,ci) = ',ws(d4,ci)
  2734 +endif
2637 2735  
2638   -
  2736 +!--Size - Vs CFL limit
  2737 + ws(taille_msn,ci) = ws(d4,ci)
  2738 +
  2739 +!--Settling velocity trait - Vs CFL limit
  2740 + cc(settl_msn,ci) = ws(d4,ci)
  2741 + ws(settl_msn,ci) = ws(d4,ci)
  2742 +
  2743 +
  2744 +
  2745 +
2639 2746 !-----------------------------------------------------------------------------------
2640 2747 ! NET PRIMARY PRODUCTION
2641 2748 !-----------------------------------------------------------------------------------
... ... @@ -3183,7 +3290,8 @@ endif !#A
3183 3290 ! endif
3184 3291  
3185 3292  
3186   -
  3293 +
  3294 + ! ws(d4,ci)= min(w_msn_m,(depth_bio/nlev)/dt_bio)
3187 3295  
3188 3296  
3189 3297 ! write(*,*)'flux_msn(ci) ',flux_msn(ci)
... ...
src/extras/bio/bio_save.F90
... ... @@ -94,20 +94,17 @@
94 94 !Prepare the new variable for the NETCDF output file as well as informations
95 95 !-----------------------------------------------------------------------------------
96 96  
97   -
98 97 !Density of the fluid at each level
99 98 iret = new_nc_variable(ncid,'rho',NF_REAL,4,dims,rho_id)
100 99 iret = set_attributes(ncid,rho_id,units='kg.m-3',long_name='Density')
101 100  
102   -
103 101 ! Sedimentation or swimming rate of particulate matter
104 102 ! Living phytoplankton
105 103 iret = new_nc_variable(ncid,'wp',NF_REAL,4,dims,wp_id)
106 104 iret = set_attributes(ncid,wp_id,units='m/day', &
107   - long_name='phytoplancton settling velocity')
108   -
  105 + long_name='Phytoplancton settling velocity')
109 106  
110   - ! Living zooplankton
  107 + ! Living zooplankton
111 108 !Nocera Model
112 109 if (bio_model .eq. 8) then
113 110 iret = new_nc_variable(ncid,'wz',NF_REAL,4,dims,wz_id)
... ... @@ -120,7 +117,7 @@
120 117 ! Living zooplankton
121 118 iret = new_nc_variable(ncid,'wz',NF_REAL,4,dims,wz_id)
122 119 iret = set_attributes(ncid,wz_id,units='m/day', &
123   - long_name='zooplancton swimming velocity')
  120 + long_name='Zooplancton swimming velocity')
124 121 ! Dead Phytoplankton
125 122 iret = new_nc_variable(ncid,'wd1',NF_REAL,4,dims,wd1_id)
126 123 iret = set_attributes(ncid,wd1_id,units='m/day', &
... ... @@ -136,8 +133,7 @@
136 133 ! Marine Snow
137 134 iret = new_nc_variable(ncid,'wd4',NF_REAL,4,dims,wd4_id)
138 135 iret = set_attributes(ncid,wd4_id,units='m/day', &
139   - long_name='Marine snow settling velocity')
140   -
  136 + long_name='Marine snow settling velocity (Calc)')
141 137  
142 138 ! Stickiness
143 139 ! 2 Living phytoplankton
... ... @@ -183,8 +179,6 @@
183 179  
184 180 end if
185 181  
186   -
187   -
188 182 !CHG4 Diagnostic du PAR
189 183 iret = new_nc_variable(ncid,'par',NF_REAL,4,dims,par_id)
190 184 iret = set_attributes(ncid,par_id,units='W/m2',long_name='PAR')
... ... @@ -207,25 +201,24 @@
207 201 iret = new_nc_variable(ncid,'ppnet',NF_REAL,4,dims,ppnet_id)
208 202 iret = set_attributes(ncid,ppnet_id,units='1/day',long_name='net primary production rate')
209 203  
210   -
  204 + !Polynow model
  205 + if (bio_model .eq. 10) then
  206 + ! Total flux of particle reaching msn
211 207 iret = new_nc_variable(ncid,'flux_msn',NF_REAL,4,dims,flux_msn_id)
212 208 iret = set_attributes(ncid,flux_msn_id,units='1/day',long_name='flux_msn')
213   -
  209 + ! Flux of Phyto going to msn
214 210 iret = new_nc_variable(ncid,'Flux_P',NF_REAL,4,dims,Flux_P_id)
215 211 iret = set_attributes(ncid,Flux_P_id,units='1/day',long_name='Flux_P')
  212 + ! Flux of DPH going to msn
216 213 iret = new_nc_variable(ncid,'Flux_D1',NF_REAL,4,dims,Flux_D1_id)
217 214 iret = set_attributes(ncid,Flux_D1_id,units='1/day',long_name='Flux_D1')
218   - iret = new_nc_variable(ncid,'Flux_D2',NF_REAL,4,dims,Flux_D2_id)
  215 + ! Flux of DZO going to msn
  216 + iret = new_nc_variable(ncid,'Flux_D2',NF_REAL,4,dims,Flux_D2_id)
219 217 iret = set_attributes(ncid,Flux_D2_id,units='1/day',long_name='Flux_D2')
220   - iret = new_nc_variable(ncid,'Flux_D3',NF_REAL,4,dims,Flux_D3_id)
  218 + ! Flux of FP going to msn
  219 + iret = new_nc_variable(ncid,'Flux_D3',NF_REAL,4,dims,Flux_D3_id)
221 220 iret = set_attributes(ncid,Flux_D3_id,units='1/day',long_name='Flux_D3')
222   -
223   - iret = new_nc_variable(ncid,'size_msnow',NF_REAL,4,dims,size_msnow_id)
224   - iret = set_attributes(ncid,size_msnow_id,units='m',long_name='size_msnow')
225   -
226   -
227   - iret = new_nc_variable(ncid,'w_msn_lev',NF_REAL,4,dims,w_msn_lev_id)
228   - iret = set_attributes(ncid,w_msn_lev_id,units='m/j',long_name='w_msn_lev')
  221 + endif
229 222  
230 223 !DD Diagnostic de npar (nb de particules lagrangiennes) pour bebogage
231 224 !iret = new_nc_variable(ncid,'npar',NF_REAL,4,dims,npar_id)
... ... @@ -248,12 +241,12 @@
248 241 end do
249 242  
250 243  
251   -!Density of the fluid at each level
252   - iret = store_data(ncid,rho_id,XYZT_SHAPE,nlev,array=rho)
  244 + !Density of the fluid at each level
  245 + iret = store_data(ncid,rho_id,XYZT_SHAPE,nlev,array=rho)
253 246  
254 247  
255 248 ! Sedimentation rate of phytoplankton
256   - iret = store_data(ncid,wp_id,XYZT_SHAPE,nlev,array=secs_pr_day*ws(1,:))
  249 + iret = store_data(ncid,wp_id,XYZT_SHAPE,nlev,array=secs_pr_day*ws(1,:))
257 250  
258 251 ! Swimming velocity of zootoplankton
259 252 !Nocera model
... ... @@ -265,6 +258,8 @@
265 258 !As example : : p=1,z=2,b=3,d1(dph)=4,n=5,a=6,l=7,d2(dzo)=8,d3(fp)=9,d4(msn)=10
266 259  
267 260 if (bio_model .eq. 10) then
  261 +
  262 +! Here it is needed to multiply by 86400(secs_pr_day) in order to have the model ouptu (in s) convert in per day
268 263 ! Living zooplankton
269 264 iret = store_data(ncid,wz_id,XYZT_SHAPE,nlev,array=secs_pr_day*ws(2,:))
270 265 ! Dead Phytoplankton
... ... @@ -307,15 +302,6 @@
307 302 ! 2 fecal pellets
308 303 iret = store_data(ncid,sti_2fp_id,XYZT_SHAPE,nlev,array=sti_2fp(:))
309 304  
310   -
311   - !
312   - iret = store_data(ncid,size_msnow_id,XYZT_SHAPE,nlev,array=size_msnow(:))
313   -
314   -
315   - !
316   - iret = store_data(ncid,w_msn_lev_id,XYZT_SHAPE,nlev,array=w_msn_lev(:))
317   -
318   -
319 305 ! Rho_F
320 306 ! iret = store_data(ncid,rho_F_id,XYZT_SHAPE,nlev,array=)
321 307  
... ... @@ -333,12 +319,13 @@
333 319 iret = store_data(ncid,ammlim2_id,XYZT_SHAPE,nlev,array=ammlim2(:))
334 320 iret = store_data(ncid,ppnet_id,XYZT_SHAPE,nlev,array=ppnet(:))
335 321  
336   - iret = store_data(ncid,flux_msn_id,XYZT_SHAPE,nlev,array=flux_msn(:))
337   - iret = store_data(ncid,Flux_P_id,XYZT_SHAPE,nlev,array=Flux_P(:))
338   - iret = store_data(ncid,Flux_D1_id,XYZT_SHAPE,nlev,array=Flux_D1(:))
339   - iret = store_data(ncid,Flux_D2_id,XYZT_SHAPE,nlev,array=Flux_D2(:))
340   - iret = store_data(ncid,Flux_D3_id,XYZT_SHAPE,nlev,array=Flux_D3(:))
341   -
  322 + if (bio_model .eq. 10) then
  323 + iret = store_data(ncid,flux_msn_id,XYZT_SHAPE,nlev,array=secs_pr_day*flux_msn(:))
  324 + iret = store_data(ncid,Flux_P_id,XYZT_SHAPE,nlev,array=secs_pr_day*Flux_P(:))
  325 + iret = store_data(ncid,Flux_D1_id,XYZT_SHAPE,nlev,array=secs_pr_day*Flux_D1(:))
  326 + iret = store_data(ncid,Flux_D2_id,XYZT_SHAPE,nlev,array=secs_pr_day*Flux_D2(:))
  327 + iret = store_data(ncid,Flux_D3_id,XYZT_SHAPE,nlev,array=secs_pr_day*Flux_D3(:))
  328 +endif
342 329  
343 330  
344 331 !DD
... ...