ROMS_r215: /home/cholod/wd/peps/ROMS/romsagrif_r215/Roms_tools/Roms

00001 !
00002 ! $Id: bbl.F,v 1.5 2005/01/19 15:01:24 pmarches Exp $
00003 !
00004 #include "cppdefs.h"
00005 #ifdef BBL
00006       subroutine bblm (tile)
00007 !====================================================================!
00008 !  Copyright (c) 2002 Rutgers/UCLA                                   !
00009 !=============================================== Meinte Blaas 2002 ==!
00010 !                                                                    !
00011 !  Compute bottom stresses for combined waves & currents             !
00012 !  using the parametric approximation by Soulsby 1997:               !
00013 !  t_cw = t_c[1+1.2(t_w/(t_c+t_w))^3.2]                              !
00014 !  in which                                                          !
00015 !     t_cw = the combined wave-averaged stress (in current dir)      !
00016 !     t_c  = stress due to currents if waves would be absent         ! 
00017 !     t_w  = amplitude of stress due to waves without currents       !
00018 !  and                                                               !
00019 !     t_cw_max = SQRT([t_cw+t_w cos(phi_cw)]^2 + [t_w sin(phi_cw)]^2)!
00020 !     in which                                                       !
00021 !       t_cw_max ~ the maximum combined wave-averaged stress         !
00022 !       phi_cw   = the angle between current and waves               !
00023 !--------------------------------------------------------------------!
00024 ! References:
00025 !
00026 ! Dyer 1986, Coastal & Estuarine Sediment Dynamics, Wiley, 342 pp.
00027 ! Harris & Wiberg 2001, Comp. & Geosci. 27, 675-690
00028 ! Li & Amos 2001, Comp. & Geosci. 27, 619-645
00029 ! Soulsby 1997, Dynamics of Marine Sands, Telford  Publ., 249 pp.
00030 ! Soulsby 1995, Bed shear-stresses due to combined waves and currents,
00031 !  in: Stive et al: Advances in Coastal Morphodynamics, Wiley, 4.20-4.23
00032 ! Wiberg & Harris 1994, J. Geophys. Res. 99(C4), 775-789
00033 !====================================================================!
00034       implicit none
00035       integer  tile, trd, omp_get_thread_num
00036 # include "param.h"
00037 # include "private_scratch.h"
00038 !
00039 # include "compute_tile_bounds.h"
00040 !
00041       trd=omp_get_thread_num()
00042       call bblm_tile (Istr,Iend,Jstr,Jend,
00043      &                  A2d(1, 1,trd),A2d(1, 2,trd),A2d(1, 3,trd),
00044      &                  A2d(1, 4,trd),A2d(1, 5,trd),A2d(1, 6,trd),
00045      &                  A2d(1, 7,trd))
00046       return
00047       end
00048 
00049 !********************************************************************!
00050       SUBROUTINE bblm_tile (Istr,Iend,Jstr,Jend,Ub,Zr,Ur,Vr,Umag,
00051      &                      tauc,tauw)
00052 !********************************************************************!
00053 !  output: bu/vstr = tauc = t_cw 
00054 !     effective wave-averaged bottom stress
00055 !     (after applying apparent roughness)  
00056 !     bu/vstrw = tauw = t_cw_max
00057 !     maximum skin-frictional bottom stress
00058 !     rheight = ripple height (m)
00059 !     rlength = ripple length (m)
00060 !                                
00061 !  input: Awave, Pwave, Dwave, h, Zob, d50, u,v
00062 !         NB! Dwave is assumed to be direction waves are coming from
00063 !             positive clockwise wrt North      
00064 !********************************************************************!
00065       implicit none
00066 # include "param.h"
00067 # include "bbl.h"
00068 # include "forces.h"
00069 # include "grid.h"
00070 # include "ocean3d.h"
00071 # include "scalars.h"
00072 # include "sediment.h"
00073 
00074       logical
00075      &        SAND, SILT
00076 
00077       integer
00078      &        Iend, Istr, Jend, Jstr, i, ised, j
00079 
00080       real
00081      &  Ab, anglec, anglew, cff, cff1, cff2, d50,
00082      &  Fwave, Kbh, Kbh2, Kdh, phic, phicw, 
00083      &  rhbio, rhbiomax, rhgt, rhmax, rhmin, rhfac, 
00084      &  rlbio, rlen, rlmin, rhosed, rhow,
00085      &  tau_bf, tau_c, tau_cb, tau_ex, tau_up,tau_w,
00086      &  tau_wb, tau_cs, tau_cw, tau_cwb, tau_cws, tau_en,
00087      &  thetw, twopi, Ucur, Vcur, Ubc, visk, wset, 
00088      &  znot, znotc, znot_bl, znot_rip
00089       parameter (Ubc = 0.01)               
00090 
00091       real
00092      &            tauc(PRIVATE_2D_SCRATCH_ARRAY),
00093      &            tauw(PRIVATE_2D_SCRATCH_ARRAY),
00094      &              Ub(PRIVATE_2D_SCRATCH_ARRAY),
00095      &            Umag(PRIVATE_2D_SCRATCH_ARRAY),
00096      &              Ur(PRIVATE_2D_SCRATCH_ARRAY),
00097      &              Vr(PRIVATE_2D_SCRATCH_ARRAY),
00098      &              Zr(PRIVATE_2D_SCRATCH_ARRAY)
00099      
00100       real
00101      &        K1, K2, K3, K4, K5, K6
00102       parameter (K1=0.6666666666, K2=0.3555555555,
00103      &           K3=0.1608465608, K4=0.0632098765,
00104      &           K5=0.0217540484, K6=0.0065407983)
00105 
00106       real
00107      &           scf1, scf2, scf3, scf4, scf5
00108       parameter (scf1 = 0.5 * 1.39, scf2 = 0.52, 
00109      &           scf3 = 2.0 - scf2, scf4 = 1.2,
00110      &           scf5 = 3.2)
00111 
00112       
00113 # include "compute_auxiliary_bounds.h"
00114 !
00115 !---------------------------------------------------------------------
00116 !  Initalize to default values + set b(u,v)str to old time level value.
00117 !---------------------------------------------------------------------
00118 !
00119       twopi=2.0*pi
00120       rhbiomax = 6.0d-3 ! maximum biogenic ripple height [m]  
00121       rhmin = 1.0d-3    ! user-defined minimum ripple height [m]  
00122       rlmin = 1.0d-2    ! user-defined minimum ripple length [m]  
00123       rhfac = 1.0d0/EXP(4.11d0)
00124 
00125       do j=JstrV-1,Jend
00126         do i=IstrU-1,Iend
00127           tauc(i,j)=0.0       ! [m^2/s^2]
00128           tauw(i,j)=0.0
00129         enddo
00130       enddo
00131 !---------------------------------------------------------------------
00132 !  Set currents above bed.
00133 !---------------------------------------------------------------------
00134       do j=JstrV-1,Jend+1
00135         do i=IstrU-1,Iend+1
00136           Zr(i,j)=z_r(i,j,1)-z_w(i,j,0)
00137           Ur(i,j)=u(i,j,1,nrhs)
00138           Vr(i,j)=v(i,j,1,nrhs)
00139         enddo
00140       enddo
00141     
00142 ! start main jV, iU loop
00143 
00144       do j=JstrV-1,Jend
00145         do i=IstrU-1,Iend
00146           rhbio = 0.0d0
00147           rlbio = 0.0d0
00148           rlen = Lripple(i,j)
00149           rhgt = Hripple(i,j)
00150           rhow = rho(i,j,1)+1000.0
00151           visk = 1.3e-3/rhow         ! kinem. viscosity          
00152 !---------------------------------------------------------------------
00153 !  Compute bed wave orbital velocity (m/s) and excursion amplitude
00154 !  (m) from wind-induced waves.  Use Dean & Dalrymple 1991 6th degree
00155 !  polynomial expr. to approximate wave number on shoaling water
00156 !---------------------------------------------------------------------
00157           Fwave=twopi/Pwave(i,j)
00158           Kdh=h(i,j)*Fwave*Fwave/g
00159           Kbh2=Kdh*Kdh+
00160      &         Kdh/(1.0+Kdh*(K1+Kdh*(K2+Kdh*(K3+Kdh*(K4+
00161      &              Kdh*(K5+K6*Kdh))))))
00162           Kbh = SQRT(Kbh2)
00163 
00164 !  Compute bed wave orbital velocity and excursion amplitude.
00165 !
00166           Ab=Awave(i,j)/SINH(Kbh)
00167           Ub(i,j)=Fwave*Ab
00168 # ifdef MASKING
00169      &                                   *rmask(i,j)
00170 # endif
00171 !
00172 !  Compute bottom current magnitude at RHO-points.
00173 !
00174           Ucur=0.5*(Ur(i,j)+Ur(i+1,j))
00175           Vcur=0.5*(Vr(i,j)+Vr(i,j+1))
00176           Umag(i,j)=SQRT(Ucur*Ucur+Vcur*Vcur)
00177 # ifdef MASKING
00178      &                                   *rmask(i,j)
00179 # endif
00180 !
00181 !  Compute angle between currents and waves (radians)
00182 !
00183           if (Ucur.ne.0.0) then
00184             phic=ATAN2(Vcur,Ucur)
00185           else
00186             phic=(pi/2.0)*SIGN(1.0,Vcur)
00187           endif
00188           phicw=(3.0*pi)/2.0 -Dwave(i,j)-phic
00189 # ifdef CURVGRID
00190      &                                  -angler(i,j)
00191 #endif
00192 ! if Dwave is dir waves travel TO:
00193 !         phicw=pi/2.0 +Dwave(i,j)-angler(i,j)-phic
00194 !---------------------------------------------------------------------
00195 !  Determine skin roughness from sediment size and
00196 !  set default logarithmic profile consistent with current-only case
00197 !  
00198 ! Establish local median grain size for all calculations in this subroutine 
00199 ! and determine local values of critical stresses.
00200 ! Since most parameterizations have been derived ignoring multiple grain sizes, 
00201 ! we apply this single d50 also in the case of mixed beds. 
00202 !---------------------------------------------------------------------
00203 # ifdef ANA_BSEDIM
00204 
00205           d50 = Ssize(i,j)             ! [m]
00206           tau_cb = taucb(i,j)/rhow     ! [m^2/s^2]
00207                                        ! Critical stress for bedload transport
00208 
00209           wset = w_set(i,j)            ! [m/s]
00210           rhosed = Sdens(i,j)/rhow     ! [ ] relative density
00211   
00212  
00213 
00214 # elif defined SEDIMENT
00215 ! Compute representative median d50, rho_grain, tau_cb, w_settle using bed_frac 
00216 ! For local use only, assume log-normal distribution of size classes
00217 
00218           d50   = 1.0
00219           tau_cb= 1.0
00220           wset  = 1.0
00221           rhosed= 1.0
00222        
00223           do ised=1,NST
00224             d50    = d50   *(Sd  (ised))**bed_frac(i,j,1,ised)
00225             tau_cb = tau_cb*(tau_ce(ised))**bed_frac(i,j,1,ised) ! [m2/s2] CHECK tau(d50)
00226             wset   = wset  *(wsed  (ised))**bed_frac(i,j,1,ised) ! [m/s]    
00227             rhosed = rhosed*(Srho  (ised))**bed_frac(i,j,1,ised)
00228           enddo 
00229 
00230           rhosed= rhosed/rhow
00231        
00232 # endif /* ANA_BSEDIM || SEDIMENT */ 
00233 
00234 ! critical stress for transition to sheet flow  (Li&Amos'01, eq.11)
00235           tau_up = 0.172*(rhosed-1.)*g*(d50**0.624) ! [m^2/s^2]
00236 
00237 ! critical stress for initiation to suspended load (Bagnold '66|Li&Amos'01, eq.9)
00238 !         tau_crs = max(0.64*wset*wset,ustr_shields)              ! [m^2/s^2]
00239 
00240 ! threshold stress for break off (Grant & Madsen,1982)   
00241           tau_bf = 0.79*(visk**(-0.6))*(((rhosed-1)*g)**0.3)
00242      &                                  *(d50**0.9)*tau_cb        ! [m^2/s^2]
00243 
00244             znotc= d50/12.0
00245             znot = MAX(Zob,znotc) 
00246              
00247 !---------------------------------------------------------------------
00248 !  set tauc acc. to current-only case (max znot & skin friction) [m/s]
00249 !---------------------------------------------------------------------
00250           cff1 = vonKar/LOG(Zr(i,j)/znot)
00251           cff2 = MIN(Cdb_max,MAX(Cdb_min,cff1*cff1))
00252           tauc(i,j)= cff2*Umag(i,j)*Umag(i,j)        ! [m^2/s^2] 
00253       
00254           cff1 = vonKar/LOG(Zr(i,j)/znotc)
00255           tau_cs = cff1*cff1*Umag(i,j)*Umag(i,j)     ! [m^2/s^2] 
00256  
00257 !---------------------------------------------------------------------
00258 ! If significant waves (Ub > Ubc= 0.01 m/s):
00259 !  Wave-current interaction case according to Soulsby 1995.  
00260 !  Otherwise: tauw = tauc for sediment purposes
00261 !---------------------------------------------------------------------
00262 
00263           if(Ub(i,j).gt.Ubc) then 
00264 
00265 ! Determine skin stresses for pure waves and combined flow (m^2/s^2)
00266 ! using Soulsby approximation of the wave friction factor
00267 !  fw = 2*scf1*(znotc/Ab)**scf2;  tauw = fw*Ub^2/2
00268  
00269            tau_w = scf1*((znotc*Fwave)**scf2)*(Ub(i,j)**scf3)
00270  
00271 !  combined wave-current skin stress 
00272             tau_cw= tau_cs*
00273      &          (1.+scf4*((tau_w/(tau_w+tau_cs))**scf5))
00274  
00275 !  Maximum of combined wave-current skin stress component for sediment
00276             tau_cws = SQRT((tau_cw+tau_w*COS(phicw))**2
00277      &                     + (tau_w*SIN(phicw))**2) 
00278             
00279             tauw(i,j)=tau_cws             ! [m^2/s^2]
00280                        
00281 !  combined wave-current stress using znot
00282             tau_w = scf1*((znot*Fwave)**scf2)*(Ub(i,j)**scf3)
00283     
00284 !  combined wave-current stress 
00285             tau_cw= tauc(i,j)*
00286      &          (1.+scf4*((tau_w/(tau_w+tauc(i,j)))**scf5))
00287 
00288 # ifdef Z0_BL
00289 !---------------------------------------------------------------------
00290 !  Compute bedload roughness for ripple predictor & sediment purposes.
00291 !  At high transport stages, friction depends on thickness of bedload layer.
00292 !  Shear stress due to combined grain & bedload roughness is to be used 
00293 !  to predict ripples and onset of suspension (Li & Amos 2001)
00294 !---------------------------------------------------------------------
00295 
00296             tau_ex=max((tau_cws-tau_cb),0.0)
00297             cff=(1./((rhosed-1.)*g*d50))
00298             znot_bl=17.4*d50*(cff*tau_ex)**0.75
00299             znotc = znotc + znot_bl
00300 
00301 !---------------------------------------------------------------------
00302 !  Compute stresses for sediment purposes [m^2/s^2], 
00303 !  using grain + bedload roughness.
00304 !---------------------------------------------------------------------
00305             cff1=vonKar/LOG(Zr(i,j)/znotc)
00306             tau_c = cff1*cff1*Umag(i,j)*Umag(i,j)    
00307             tau_wb = scf1*((znotc*Fwave)**scf2)*(Ub(i,j)**scf3)
00308             tau_cw= tau_c*(1.+scf4*((tau_wb/(tau_wb+tau_c))**scf5))
00309         
00310 !  Maximum of combined wave-current stress component for sediment purposes
00311 
00312             tau_cwb = SQRT((tau_cw+tau_wb*COS(phicw))**2 +
00313      &                   (tau_wb*SIN(phicw))**2)
00314             tauw(i,j) = tau_cwb                 ! [m^2/s^2]
00315     
00316 
00317 #  ifdef Z0_RIP        
00318 !---------------------------------------------------------------------
00319 !  Determine bedform roughness  ripple height (m)
00320 !  and ripple length (m) for sandy beds.
00321 !
00322 !  Use structure according to Li & Amos 2001
00323 !  See Li & Amos 2001 for further references
00324 !---------------------------------------------------------------------
00325 
00326 ! check median grain diameter
00327 
00328              if (d50.ge.0.063d-3) then
00329     
00330 !enhanced skin stress if pre-exisiting ripples (Nielsen, 1986)
00331               rhmax=0.8*rlen/pi
00332               rhgt=max(min(rhmax,rhgt),rhmin)
00333           
00334               tau_en=max(tau_cws,tau_cws*(rlen/(rlen-pi*rhgt))**2) 
00335           
00336               if ((tau_cws.lt.tau_cb).and.(tau_en.ge.tau_cb)) then 
00337 ! local transp.
00338                 rhgt = (19.6*(SQRT(tau_cws/tau_cb))+20.9)*d50
00339                 rlen = rhgt/0.12
00340               else 
00341      &          if ((tau_cws.ge.tau_cb).and.(tau_cwb.lt.tau_bf)) then 
00342 ! bed load in eq. range
00343                   rhgt = (22.15*(SQRT(tau_cwb/tau_cb))+6.38)*d50
00344                   rlen = rhgt/0.12
00345               else 
00346      &          if ((tau_cwb.ge.tau_bf).and.(tau_cwb.lt.tau_up)) then 
00347 ! break off regime
00348                   rlen = 535.*d50
00349                   rhgt = 0.15*rlen*(SQRT(tau_up)-SQRT(tau_cwb))/
00350      &                          (SQRT(tau_up)-SQRT(tau_bf )) 
00351               else if  (tau_cwb.ge.tau_up) then 
00352 ! sheet flow, plane bed
00353                 rlen = 0.0
00354                 rhgt = 0.0
00355               else 
00356 ! tau_en < tau_cb: no transport, pre-existing conditions
00357                 rhgt=Hripple(i,j) 
00358                 rlen=Lripple(i,j) 
00359               endif !tau_cws
00360             endif !d50
00361        
00362 #  endif /* Z0_BL */
00363 # endif /* Z0_RIP */       
00364 
00365 # ifdef Z0_BIO
00366 !---------------------------------------------------------------------
00367 !  Determine (biogenic) bedform roughness  ripple height (m)
00368 !  and ripple length (m) for silty beds 
00369 !  using Harris & Wiberg 2001.
00370 !---------------------------------------------------------------------
00371 !   rlbio -> user defined or input field...
00372             if (d50.lt.0.063d-3) then
00373               rlbio = 0.1d0      ! 10 cm default biogenic r.length (Wheatcroft 1994) 
00374               thetw = tau_cws*(1./((rhosed-1.)*g*d50))
00375               rhbio = (thetw**(-1.67d0))*rlbio*rhfac
00376               rhgt = min(rhbio,rhbiomax)
00377               rlen = rlbio
00378             endif
00379 
00380 # endif /* Z0_BIO */
00381 
00382 # if defined Z0_RIP || defined Z0_BIO
00383 !ripple roughness using Grant & Madsen (1982) roughness length
00384 
00385             znot_rip = 0.92*rhgt*rhgt/(max(rlen,rlmin))
00386             znotc = znotc+znot_rip
00387 !---------------------------------------------------------------------
00388 ! Compute bottom stress components based on total roughnes [m/s]^2
00389 !---------------------------------------------------------------------
00390             cff1  = vonKar/LOG(Zr(i,j)/znotc)
00391             tau_c = cff1*cff1*Umag(i,j)*Umag(i,j)  
00392             tau_w = scf1*((znotc*Fwave)**scf2)*(Ub(i,j)**scf3)
00393             tau_cw= tau_c*(1.+scf4*((tau_w/(tau_w+tau_c))**scf5))
00394 
00395 # endif /* Z0_RIP || Z0_BIO */
00396 !---------------------------------------------------------------------
00397 !  Compute effective bottom shear velocity [m/s]
00398 !  relevant for flow and eddy-diffusivities/viscosity
00399 !---------------------------------------------------------------------
00400 
00401             tauc(i,j) = tau_cw
00402           
00403           else ! (Ub < Ubc)   
00404 
00405 ! If current-only: tauw = tauc(skin) for use in sediment.F (ifdef BBL)
00406 ! tauc for current still depending on roughness due to current ripples (ifdef Z0_RIP)
00407             
00408             tauw(i,j) = tau_cs
00409 
00410 # ifdef Z0_RIP    
00411             if(tau_cs.gt.tau_up) then
00412               rhgt=0.0
00413               rlen=0.0
00414             else if(tau_cs.lt.tau_cb) then
00415               rhgt=Hripple(i,j) 
00416               rlen=Lripple(i,j) 
00417             else
00418               rlen=1000.*d50               ! Yalin (1964)
00419               rhgt=7.4*(rlen/100.)**1.19 ! Allen (1970)
00420             endif
00421  
00422             znotc = znotc+0.92*rhgt*rhgt/(max(rlen,rlmin))
00423         
00424             cff1=vonKar/LOG(Zr(i,j)/znotc)
00425             cff2=MIN(Cdb_max,MAX(Cdb_min,cff1*cff1))
00426             tauc(i,j) = cff2*Umag(i,j)*Umag(i,j)
00427 # endif 
00428           endif 
00429 !---------------------------------------------------------------------
00430 !  Load variables for output purposes.
00431 !---------------------------------------------------------------------
00432 
00433           Abed(i,j)   = Ab
00434 # ifdef MASKING
00435      &                                   *rmask(i,j)
00436 # endif
00437           Hripple(i,j)= rhgt
00438 # ifdef MASKING
00439      &                                   *rmask(i,j)
00440 # endif
00441           Lripple(i,j)= rlen
00442 # ifdef MASKING
00443      &                                   *rmask(i,j)
00444 # endif
00445          Zbnot(i,j)  = znot
00446 # ifdef MASKING
00447      &                                   *rmask(i,j)
00448 # endif
00449          Zbapp(i,j) = znotc
00450 # ifdef MASKING
00451      &                                   *rmask(i,j)
00452 # endif
00453        enddo
00454       enddo  ! end main loop
00455 
00456 !---------------------------------------------------------------------
00457 !  Compute kinematic bottom stress components for flow
00458 !  due to combined current and wind-induced waves.
00459 !---------------------------------------------------------------------
00460 !
00461       do j=Jstr,Jend
00462         do i=IstrU,Iend
00463           cff1=0.5*(tauc(i-1,j)+tauc(i,j))
00464           cff2=0.5*(tauw(i-1,j)+tauw(i,j))
00465           anglec=Ur(i,j)/(0.5*(Umag(i-1,j)+Umag(i,j)+1.e-10))
00466           anglew=COS(0.5*(Dwave(i-1,j)+Dwave(i,j)))
00467           bustr(i,j)=cff1*anglec
00468 # ifdef MASKING
00469      &                                   *umask(i,j)
00470 # endif
00471           bustrw(i,j)=cff2*anglew
00472 # ifdef MASKING
00473      &                                   *umask(i,j)
00474 # endif
00475           Ubed(i,j)=Ub(i,j)*anglew
00476 # ifdef MASKING
00477      &                                   *umask(i,j)
00478 # endif
00479 
00480         enddo
00481       enddo
00482       do j=JstrV,Jend
00483         do i=Istr,Iend
00484           cff1=0.5*(tauc(i,j-1)+tauc(i,j))
00485           cff2=0.5*(tauw(i,j-1)+tauw(i,j))
00486           anglec=Vr(i,j)/(0.5*(Umag(i,j-1)+Umag(i,j)+1.e-10))
00487           anglew=SIN(0.5*(Dwave(i,j-1)+Dwave(i,j)))
00488           bvstr(i,j)=cff1*anglec
00489 # ifdef MASKING
00490      &                                   *vmask(i,j)
00491 # endif
00492           bvstrw(i,j)=cff2*anglew
00493 # ifdef MASKING
00494      &                                   *vmask(i,j)
00495 # endif
00496           Vbed(i,j)=Ub(i,j)*anglew
00497 # ifdef MASKING
00498      &                                   *vmask(i,j)
00499 # endif
00500         enddo
00501       enddo
00502 !
00503 !
00504       return
00505       end
00506 !====================================================================!
00507 #else
00508       subroutine bblm_empty
00509       end
00510 #endif
00511
/home/cholod/wd/peps/ROMS/romsagrif_r215/Roms_tools/Roms_Agrif/bbl.F
