Changeset 7853 for branches/UKMO/r6232_HZG_WAVE-coupling/NEMOGCM/NEMO
- Timestamp:
- 2017-03-30T16:22:25+02:00 (7 years ago)
- Location:
- branches/UKMO/r6232_HZG_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/UKMO/r6232_HZG_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave.F90
r7809 r7853 226 226 IF( ln_phioc .AND. .NOT. cpl_phioc ) THEN !== Wave to ocean energy ==! 227 227 CALL fld_read( kt, nn_fsbc, sf_phioc ) ! read wave to ocean energy from external forcing 228 rn_crban(:,:) = 29.0 * sf_phioc(1)%fnow(:,:,1) ! ! Alfa is phioc*sqrt(r w/ra)sbc_wa228 rn_crban(:,:) = 29.0 * sf_phioc(1)%fnow(:,:,1) ! ! Alfa is phioc*sqrt(rau0/zrhoa) : rau0=water density, zhroa= air density 229 229 WHERE( rn_crban < 10.0 ) rn_crban = 10.0 230 230 WHERE( rn_crban > 300.0 ) rn_crban = 300.0 … … 235 235 IF( jpfld > 0 ) THEN ! Read from file only if the field is not coupled 236 236 CALL fld_read( kt, nn_fsbc, sf_sd ) ! read wave parameters from external forcing 237 IF( jp_hsw > 0 ) hsw (:,:) = sf_sd(jp_hsw)%fnow(:,:,1) ! significant wave height 238 IF( jp_wmp > 0 ) wmp (:,:) = sf_sd(jp_wmp)%fnow(:,:,1) ! wave mean period 239 IF( jp_usd > 0 ) ut0sd(:,:) = sf_sd(jp_usd)%fnow(:,:,1) ! 2D zonal Stokes Drift at T point 240 IF( jp_vsd > 0 ) vt0sd(:,:) = sf_sd(jp_vsd)%fnow(:,:,1) ! 2D meridional Stokes Drift at T point 237 IF( jp_hsw > 0 ) THEN 238 hsw (:,:) = sf_sd(jp_hsw)%fnow(:,:,1) ! significant wave height 239 WHERE( hsw > 100.0 ) hsw = 0.0 240 WHERE( hsw < 0.0 ) hsw = 0.0 241 ENDIF 242 IF( jp_wmp > 0 ) THEN 243 wmp (:,:) = sf_sd(jp_wmp)%fnow(:,:,1) ! wave mean period 244 WHERE( wmp > 100.0 ) wmp = 0.0 245 WHERE( wmp < 0.0 ) wmp = 0.0 246 ENDIF 247 IF( jp_usd > 0 ) THEN 248 ut0sd(:,:) = sf_sd(jp_usd)%fnow(:,:,1) ! 2D zonal Stokes Drift at T point 249 WHERE( ut0sd < -100.0 ) ut0sd = 1.0 250 WHERE( ut0sd > 100.0 ) ut0sd = 1.0 251 ENDIF 252 IF( jp_vsd > 0 ) THEN 253 vt0sd(:,:) = sf_sd(jp_vsd)%fnow(:,:,1) ! 2D meridional Stokes Drift at T point 254 WHERE( vt0sd < -100.0 ) vt0sd = 1.0 255 WHERE( vt0sd > 100.0 ) vt0sd = 1.0 256 ENDIF 241 257 ENDIF 242 258 ! -
branches/UKMO/r6232_HZG_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90
r7809 r7853 52 52 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: rsbc_tke3 53 53 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: rsbc_psi1 54 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: rtrans 54 55 55 56 ! !! ** Namelist namzdf_gls ** … … 61 62 INTEGER :: nn_stab_func ! stability functions G88, KC or Canuto (=0/1/2) 62 63 INTEGER :: nn_clos ! closure 0/1/2/3 MY82/k-eps/k-w/gen 64 INTEGER :: nn_wmix ! type of wave breaking mixing 65 INTEGER, PUBLIC, PARAMETER :: jp_craigbanner = 0 ! Craig and Banner formulation (original NEMO formulation - 66 ! direct conversion of mechanical to turbulent energy) 67 INTEGER, PUBLIC, PARAMETER :: jp_janssen = 1 ! Janssen formulation - no assumption on direct energy conversion 63 68 REAL(wp) :: rn_clim_galp ! Holt 2008 value for k-eps: 0.267 64 69 REAL(wp) :: rn_epsmin ! minimum value of dissipation (m2/s3) … … 96 101 REAL(wp) :: rm7 = 0.0_wp 97 102 REAL(wp) :: rm8 = 0.318_wp 98 REAL(wp) :: rtrans 103 REAL(wp) :: rtrans_default = 0.1_wp 99 104 REAL(wp) :: rc02, rc02r, rc03, rc04 ! coefficients deduced from above parameters 100 105 REAL(wp) :: rsbc_tke1_default, rsbc_tke2, rfact_tke ! - - - - … … 124 129 & ustars2(jpi,jpj) , ustarb2(jpi,jpj) , & 125 130 & rsbc_tke1(jpi,jpj), rsbc_tke3(jpi,jpj) , & 126 & rsbc_psi1(jpi,jpj), STAT= zdf_gls_alloc )131 & rsbc_psi1(jpi,jpj), rtrans(jpi,jpj), STAT= zdf_gls_alloc ) 127 132 ! 128 133 IF( lk_mpp ) CALL mpp_sum ( zdf_gls_alloc ) … … 172 177 ! variable initialization 173 178 IF( ln_phioc ) THEN 174 rsbc_tke1(:,:) = (-rsc_tke*rn_crban(:,:)/(rcm_sf*ra_sf*rl_sf))**(2._wp/3._wp) ! k_eps = 53.Dirichlet + Wave breaking 175 rsbc_tke3(:,:) = rdt * rn_crban(:,:) ! Neumann + Wave breaking 176 rsbc_psi1(:,:) = rc0**rpp * rsbc_tke1(:,:)**rmm * rl_sf**rnn ! Dirichlet + Wave breaking 179 IF( nn_wmix==jp_janssen ) THEN 180 rsbc_tke1(:,:) = (-rsc_tke*rn_crban(:,:)/(rcm_sf*ra_sf*rl_sf))**(2._wp/3._wp) ! k_eps = 53.Dirichlet + Wave breaking 181 rsbc_tke3(:,:) = rdt * rn_crban(:,:) ! Neumann + Wave breaking 182 rsbc_psi1(:,:) = rc0**rpp * rsbc_tke1(:,:)**rmm * rl_sf**rnn ! Dirichlet + Wave breaking 183 ELSE 184 rsbc_tke1(:,:) = -3._wp/2._wp*rn_crban(:,:)*ra_sf*rl_sf 185 rsbc_tke3(:,:) = rdt * rn_crban(:,:) / rl_sf 186 rtrans(:,:) = 0.2_wp / MAX( rsmall, rsbc_tke1(:,:)**(1./(-ra_sf*3._wp/2._wp))-1._wp ) 187 ENDIF 177 188 ENDIF 178 189 … … 336 347 CASE ( 0 ) ! Dirichlet case 337 348 IF( ln_phioc ) THEN ! wave induced mixing case with forced/coupled fields 338 ! First level 339 en(:,:,1) = MAX( rsbc_tke1(:,:) * ustars2(:,:), rn_emin ) 340 z_elem_a(:,:,1) = en(:,:,1) 341 z_elem_c(:,:,1) = 0._wp 342 z_elem_b(:,:,1) = 1._wp 343 344 ! One level below 345 en(:,:,2) = MAX( rsbc_tke1(:,:) * ustars2(:,:) * ((zhsro(:,:)+fsdepw(:,:,2))/zhsro(:,:) )**ra_sf, rn_emin ) 346 z_elem_a(:,:,2) = 0._wp 347 z_elem_c(:,:,2) = 0._wp 348 z_elem_b(:,:,2) = 1._wp 349 IF( nn_wmix==jp_janssen ) THEN 350 ! First level 351 en(:,:,1) = MAX( rsbc_tke1(:,:) * ustars2(:,:), rn_emin ) 352 z_elem_a(:,:,1) = en(:,:,1) 353 z_elem_c(:,:,1) = 0._wp 354 z_elem_b(:,:,1) = 1._wp 355 356 ! One level below 357 en(:,:,2) = MAX( rsbc_tke1(:,:) * ustars2(:,:) * ((zhsro(:,:)+fsdepw(:,:,2))/zhsro(:,:) )**ra_sf, rn_emin ) 358 z_elem_a(:,:,2) = 0._wp 359 z_elem_c(:,:,2) = 0._wp 360 z_elem_b(:,:,2) = 1._wp 361 ELSE 362 en(:,:,1) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1(:,:))**(2._wp/3._wp) 363 en(:,:,1) = MAX(en(:,:,1), rn_emin) 364 z_elem_a(:,:,1) = en(:,:,1) 365 z_elem_c(:,:,1) = 0._wp 366 z_elem_b(:,:,1) = 1._wp 367 ! 368 ! One level below 369 en(:,:,2) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1(:,:) * ((zhsro(:,:)+fsdepw(:,:,2)) & 370 & / zhsro(:,:) )**(1.5_wp*ra_sf))**(2._wp/3._wp) 371 en(:,:,2) = MAX(en(:,:,2), rn_emin ) 372 z_elem_a(:,:,2) = 0._wp 373 z_elem_c(:,:,2) = 0._wp 374 z_elem_b(:,:,2) = 1._wp 375 ! 376 ENDIF 349 377 ELSE ! wave induced mixing case with default values 350 378 en(:,:,1) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1_default)**(2._wp/3._wp) … … 366 394 CASE ( 1 ) ! Neumann boundary condition on d(e)/dz 367 395 IF( ln_phioc ) THEN ! Shear free case: d(e)/dz=Fw with forced/coupled fields 368 ! Dirichlet conditions at k=1 369 en(:,:,1) = MAX( rsbc_tke1(:,:) * ustars2(:,:), rn_emin ) 370 z_elem_a(:,:,1) = en(:,:,1) 371 z_elem_c(:,:,1) = 0._wp 372 z_elem_b(:,:,1) = 1._wp 373 ! at k=2, set de/dz=Fw 374 z_elem_b(:,:,2) = z_elem_b(:,:,2) + z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b 375 z_elem_a(:,:,2) = 0._wp 376 zflxs(:,:) = rsbc_tke3(:,:) * ustars2(:,:)**1.5_wp * ((zhsro(:,:)+fsdept(:,:,1) ) / zhsro(:,:) )**(1.5*ra_sf) 377 en(:,:,2) = en(:,:,2) + zflxs(:,:) / fse3w(:,:,2) 396 IF( nn_wmix==jp_janssen ) THEN 397 ! Dirichlet conditions at k=1 398 en(:,:,1) = MAX( rsbc_tke1(:,:) * ustars2(:,:), rn_emin ) 399 z_elem_a(:,:,1) = en(:,:,1) 400 z_elem_c(:,:,1) = 0._wp 401 z_elem_b(:,:,1) = 1._wp 402 ! at k=2, set de/dz=Fw 403 z_elem_b(:,:,2) = z_elem_b(:,:,2) + z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b 404 z_elem_a(:,:,2) = 0._wp 405 zflxs(:,:) = rsbc_tke3(:,:) * ustars2(:,:)**1.5_wp * ((zhsro(:,:)+fsdept(:,:,1) ) / zhsro(:,:) )**(1.5*ra_sf) 406 en(:,:,2) = en(:,:,2) + zflxs(:,:) / fse3w(:,:,2) 407 ELSE 408 ! Dirichlet conditions at k=1 409 en(:,:,1) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1(:,:))**(2._wp/3._wp) 410 en(:,:,1) = MAX(en(:,:,1), rn_emin) 411 z_elem_a(:,:,1) = en(:,:,1) 412 z_elem_c(:,:,1) = 0._wp 413 z_elem_b(:,:,1) = 1._wp 414 ! 415 ! at k=2, set de/dz=Fw 416 !cbr 417 z_elem_b(:,:,2) = z_elem_b(:,:,2) + z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b 418 z_elem_a(:,:,2) = 0._wp 419 zkar(:,:) = (rl_sf + (vkarmn-rl_sf)*(1.-exp(-rtrans(:,:)*fsdept(:,:,1)/zhsro(:,:)) )) 420 zflxs(:,:) = rsbc_tke3(:,:) * ustars2(:,:)**1.5_wp * zkar(:,:) & 421 & * ((zhsro(:,:)+fsdept(:,:,1))/zhsro(:,:) )**(1.5_wp*ra_sf) 422 423 en(:,:,2) = en(:,:,2) + zflxs(:,:)/fse3w(:,:,2) 424 ! 425 ENDIF 378 426 ELSE ! Shear free case: d(e)/dz=Fw with default values 379 427 ! Dirichlet conditions at k=1 … … 388 436 z_elem_b(:,:,2) = z_elem_b(:,:,2) + z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b 389 437 z_elem_a(:,:,2) = 0._wp 390 zkar(:,:) = (rl_sf + (vkarmn-rl_sf)*(1.-exp(-rtrans *fsdept(:,:,1)/zhsro(:,:)) ))438 zkar(:,:) = (rl_sf + (vkarmn-rl_sf)*(1.-exp(-rtrans_default*fsdept(:,:,1)/zhsro(:,:)) )) 391 439 zflxs(:,:) = rsbc_tke2 * ustars2(:,:)**1.5_wp * zkar(:,:) & 392 440 & * ((zhsro(:,:)+fsdept(:,:,1))/zhsro(:,:) )**(1.5_wp*ra_sf) … … 586 634 ! balance between the production and the 587 635 ! dissipation terms including the wave effect 588 ! Surface value 589 zdep(:,:) = zhsro(:,:) * rl_sf ! Cosmetic 590 psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 591 z_elem_a(:,:,1) = psi(:,:,1) 592 z_elem_c(:,:,1) = 0._wp 593 z_elem_b(:,:,1) = 1._wp 594 ! 595 ! One level below 596 zex1 = (rmm*ra_sf+rnn) 597 zex2 = (rmm*ra_sf) 598 zdep(:,:) = ( (zhsro(:,:) + fsdepw(:,:,2))**zex1 ) / zhsro(:,:)**zex2 599 psi (:,:,2) = rsbc_psi1(:,:) * ustars2(:,:)**rmm * zdep(:,:) * tmask(:,:,1) 600 z_elem_a(:,:,2) = 0._wp 601 z_elem_c(:,:,2) = 0._wp 602 z_elem_b(:,:,2) = 1._wp 603 ! 604 ! 636 IF( nn_wmix==jp_janssen ) THEN 637 ! Surface value 638 zdep(:,:) = zhsro(:,:) * rl_sf ! Cosmetic 639 psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 640 z_elem_a(:,:,1) = psi(:,:,1) 641 z_elem_c(:,:,1) = 0._wp 642 z_elem_b(:,:,1) = 1._wp 643 ! 644 ! One level below 645 zex1 = (rmm*ra_sf+rnn) 646 zex2 = (rmm*ra_sf) 647 zdep(:,:) = ( (zhsro(:,:) + fsdepw(:,:,2))**zex1 ) / zhsro(:,:)**zex2 648 psi (:,:,2) = rsbc_psi1(:,:) * ustars2(:,:)**rmm * zdep(:,:) * tmask(:,:,1) 649 z_elem_a(:,:,2) = 0._wp 650 z_elem_c(:,:,2) = 0._wp 651 z_elem_b(:,:,2) = 1._wp 652 ! 653 ! 654 ELSE 655 ! Surface value 656 zdep(:,:) = zhsro(:,:) * rl_sf ! Cosmetic 657 psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 658 z_elem_a(:,:,1) = psi(:,:,1) 659 z_elem_c(:,:,1) = 0._wp 660 z_elem_b(:,:,1) = 1._wp 661 ! 662 ! One level below 663 zkar(:,:) = (rl_sf + (vkarmn-rl_sf)*(1._wp-exp(-rtrans(:,:)*fsdepw(:,:,2)/zhsro(:,:) ))) 664 zdep(:,:) = (zhsro(:,:) + fsdepw(:,:,2)) * zkar(:,:) 665 psi (:,:,2) = rc0**rpp * en(:,:,2)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 666 z_elem_a(:,:,2) = 0._wp 667 z_elem_c(:,:,2) = 0._wp 668 z_elem_b(:,:,2) = 1._wp 669 ! 670 ! 671 ENDIF 605 672 ELSE ! Wave induced mixing case with default values 606 673 ! Surface value … … 612 679 ! 613 680 ! One level below 614 zkar(:,:) = (rl_sf + (vkarmn-rl_sf)*(1._wp-exp(-rtrans *fsdepw(:,:,2)/zhsro(:,:) )))681 zkar(:,:) = (rl_sf + (vkarmn-rl_sf)*(1._wp-exp(-rtrans_default*fsdepw(:,:,2)/zhsro(:,:) ))) 615 682 zdep(:,:) = (zhsro(:,:) + fsdepw(:,:,2)) * zkar(:,:) 616 683 psi (:,:,2) = rc0**rpp * en(:,:,2)**rmm * zdep(:,:)**rnn * tmask(:,:,1) … … 623 690 CASE ( 1 ) ! Neumann boundary condition on d(psi)/dz 624 691 IF( ln_phioc ) THEN ! Wave induced mixing case with forced/coupled fields 625 ! 626 zdep(:,:) = rl_sf * zhsro(:,:) 627 psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 628 z_elem_a(:,:,1) = psi(:,:,1) 629 z_elem_c(:,:,1) = 0._wp 630 z_elem_b(:,:,1) = 1._wp 631 ! 632 ! Neumann condition at k=2 633 z_elem_b(:,:,2) = z_elem_b(:,:,2) + z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b 634 z_elem_a(:,:,2) = 0._wp 635 ! 636 ! Set psi vertical flux at the surface: 637 zdep(:,:) = (zhsro(:,:) + fsdept(:,:,1))**(rmm*ra_sf+rnn-1._wp) / zhsro(:,:)**(rmm*ra_sf) 638 zflxs(:,:) = rsbc_psi3 * ( zwall_psi(:,:,1)*avm(:,:,1) + zwall_psi(:,:,2)*avm(:,:,2) ) & 639 & * en(:,:,1)**rmm * zdep 640 psi(:,:,2) = psi(:,:,2) + zflxs(:,:) / fse3w(:,:,2) 641 ! 692 IF( nn_wmix==jp_janssen ) THEN 693 ! 694 zdep(:,:) = rl_sf * zhsro(:,:) 695 psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 696 z_elem_a(:,:,1) = psi(:,:,1) 697 z_elem_c(:,:,1) = 0._wp 698 z_elem_b(:,:,1) = 1._wp 699 ! 700 ! Neumann condition at k=2 701 z_elem_b(:,:,2) = z_elem_b(:,:,2) + z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b 702 z_elem_a(:,:,2) = 0._wp 703 ! 704 ! Set psi vertical flux at the surface: 705 zdep(:,:) = (zhsro(:,:) + fsdept(:,:,1))**(rmm*ra_sf+rnn-1._wp) / zhsro(:,:)**(rmm*ra_sf) 706 zflxs(:,:) = rsbc_psi3 * ( zwall_psi(:,:,1)*avm(:,:,1) + zwall_psi(:,:,2)*avm(:,:,2) ) & 707 & * en(:,:,1)**rmm * zdep 708 psi(:,:,2) = psi(:,:,2) + zflxs(:,:) / fse3w(:,:,2) 709 ! 710 ELSE 711 ! Surface value: Dirichlet 712 zdep(:,:) = zhsro(:,:) * rl_sf 713 psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 714 z_elem_a(:,:,1) = psi(:,:,1) 715 z_elem_c(:,:,1) = 0._wp 716 z_elem_b(:,:,1) = 1._wp 717 ! 718 ! Neumann condition at k=2 719 z_elem_b(:,:,2) = z_elem_b(:,:,2) + z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b 720 z_elem_a(:,:,2) = 0._wp 721 ! 722 ! Set psi vertical flux at the surface: 723 zkar(:,:) = rl_sf + (vkarmn-rl_sf)*(1._wp-exp(-rtrans(:,:)*fsdept(:,:,1)/zhsro(:,:) )) ! Lengh scale slope 724 zdep(:,:) = ((zhsro(:,:) + fsdept(:,:,1)) / zhsro(:,:))**(rmm*ra_sf) 725 zflxs(:,:) = (rnn + rsbc_tke1(:,:) * (rnn + rmm*ra_sf) * zdep(:,:)) * & 726 (1._wp + rsbc_tke1(:,:) * zdep(:,:))**(2._wp*rmm/3._wp-1_wp) 727 zdep(:,:) = rsbc_psi1(:,:) * (zwall_psi(:,:,1)*avm(:,:,1)+zwall_psi(:,:,2)*avm(:,:,2)) * & 728 & ustars2(:,:)**rmm * zkar(:,:)**rnn * (zhsro(:,:) + fsdept(:,:,1))**(rnn-1.) 729 zflxs(:,:) = zdep(:,:) * zflxs(:,:) 730 psi(:,:,2) = psi(:,:,2) + zflxs(:,:) / fse3w(:,:,2) 731 ! 732 ENDIF 642 733 ELSE ! Wave induced mixing case with default values 643 734 ! Surface value: Dirichlet … … 653 744 ! 654 745 ! Set psi vertical flux at the surface: 655 zkar(:,:) = rl_sf + (vkarmn-rl_sf)*(1._wp-exp(-rtrans *fsdept(:,:,1)/zhsro(:,:) )) ! Lengh scale slope746 zkar(:,:) = rl_sf + (vkarmn-rl_sf)*(1._wp-exp(-rtrans_default*fsdept(:,:,1)/zhsro(:,:) )) ! Lengh scale slope 656 747 zdep(:,:) = ((zhsro(:,:) + fsdept(:,:,1)) / zhsro(:,:))**(rmm*ra_sf) 657 748 zflxs(:,:) = (rnn + rsbc_tke1_default * (rnn + rmm*ra_sf) * zdep(:,:)) * & … … 956 1047 & rn_crban_default, rn_charn, rn_frac_hs,& 957 1048 & nn_bc_surf, nn_bc_bot, nn_z0_met, & 958 & nn_stab_func, nn_clos 1049 & nn_stab_func, nn_clos, nn_wmix 959 1050 !!---------------------------------------------------------- 960 1051 ! … … 1215 1306 rsbc_tke2 = rdt * rn_crban_default / rl_sf 1216 1307 zcr = MAX( rsmall, rsbc_tke1_default**(1./(-ra_sf*3._wp/2._wp))-1._wp ) 1217 rtrans 1308 rtrans_default = 0.2_wp / zcr 1218 1309 rsbc_zs1 = rn_charn/grav ! Charnock formula for surface roughness 1219 1310 rsbc_zs2 = rn_frac_hs / 0.85_wp / grav * 665._wp ! Rascle formula for surface roughness
Note: See TracChangeset
for help on using the changeset viewer.