ROMS_r215: /home/cholod/wd/peps/ROMS/romsagrif_r215/Roms_tools/Roms_Agrif/JUNK/lmd_bkpp

00001 !
00002 ! $Id: lmd_bkpp.F,v 1.6 2005/10/06 08:46:29 pmarches Exp $
00003 !
00004 #include "cppdefs.h"
00005 #ifdef LMD_BKPP
00006       subroutine lmd_bkpp_tile (Istr,Iend,Jstr,Jend, Kv,Kt,Ks,
00007      &                                          bl_dpth,Bflux,
00008      &                     Gm1,dGm1dS, Gt1,dGt1dS, Gs1,dGs1dS,
00009      &                                    wrk1,wrk2,wrk3, Rib)
00010       implicit none
00011 # include "param.h"
00012 # include "grid.h"
00013 # include "ocean3d.h"
00014 # include "forces.h"
00015 # include "mixing.h"
00016 # include "scalars.h"
00017       integer Istr,Iend,Jstr,Jend, i,j,k, ka,ku,ksave
00018       real    Kv(PRIVATE_2D_SCRATCH_ARRAY,0:N),
00019      &        Kt(PRIVATE_2D_SCRATCH_ARRAY,0:N),
00020      &        Ks(PRIVATE_2D_SCRATCH_ARRAY,0:N),
00021      &       bl_dpth(PRIVATE_2D_SCRATCH_ARRAY), 
00022      &         Bflux(PRIVATE_2D_SCRATCH_ARRAY),
00023      &           Gm1(PRIVATE_2D_SCRATCH_ARRAY),
00024      &        dGm1dS(PRIVATE_2D_SCRATCH_ARRAY),
00025      &           Gt1(PRIVATE_2D_SCRATCH_ARRAY),
00026      &        dGt1dS(PRIVATE_2D_SCRATCH_ARRAY),
00027      &           Gs1(PRIVATE_2D_SCRATCH_ARRAY),
00028      &        dGs1dS(PRIVATE_2D_SCRATCH_ARRAY),
00029      &          wrk1(PRIVATE_2D_SCRATCH_ARRAY),
00030      &          wrk2(PRIVATE_2D_SCRATCH_ARRAY),
00031      &          wrk3(PRIVATE_2D_SCRATCH_ARRAY),
00032      &         Rib(PRIVATE_2D_SCRATCH_ARRAY,2)
00033 # undef BKPP_BFLUX
00034 # define tind nstp
00035       real Vtc,    hekman,  hmonob, dVsq,   Vtsq,
00036      &             bl_dnew, sig,    Kv_bl,  Kt_bl,  Ks_bl,
00037      &     cff,    lmd_a1, dKv_bl, dKt_bl, dKs_bl,
00038      &     cff_up, lmd_a2, Gm,     Gt,     Gs,
00039      &     cff_dn, Ritop,   lmd_a3
00040       real zbl, zsbl
00041       real eps
00042       parameter (eps=1.E-20)
00043       real lmd_cs, lmd_Cv, Ric, lmd_betaT, lmd_epsilon,
00044      &     lmd_cekman, lmd_nu0c
00045       parameter (
00046      &   lmd_cs=98.96,     ! see parameter associated with turbulent
00047                            ! velocity scales in lmd_wscale.F
00048                            !
00049      &   lmd_Cv=1.8,       ! Ratio of interior Brunt-Vaisala
00050                            ! frequency "N" to that at the
00051                            ! entrainment depth "he".
00052                            !
00053      &   Ric=0.3,          ! Critical bulk Richardson number.
00054                            !
00055      &   lmd_betaT=-0.2,   ! Ratio of entrainment flux to 
00056                            ! to surface buoyancy flux.
00057                            !
00058      &   lmd_epsilon=0.1,  ! Nondimensional extent of the
00059                            ! bottom layer.
00060                            !
00061      &   lmd_cekman=0.7,   ! Constant used in the computation
00062                            ! of Ekman depth.
00063                            !
00064      &   lmd_nu0c=0.1      ! Maximum interior convective
00065                            ! viscosity and  diffusivity due
00066                            ! to shear instability, [m^2/s];
00067      &                                                           )
00068 
00069 !
00070 ! Initialize relevant parameters
00071 !
00072       Vtc=lmd_Cv*sqrt(-lmd_betaT)/( sqrt(lmd_cs*lmd_epsilon)
00073      &                                        *Ric*vonKar*vonKar )
00074 !
00075 !  Get approximation of surface layer depth using "lmd_eps" and
00076 !  boundary layer depth from previous time step.
00077 !
00078       do j=Jstr,Jend
00079         do i=Istr,Iend
00080           bl_dpth(i,j)=lmd_epsilon*hbbl(i,j)
00081         enddo
00082       enddo
00083 !
00084 !  Compute turbulent friction velocity [m/s] "ustar" from bottom stress
00085 !  at RHO-points. (Bottom stress includes wave effects if activated in BBL)
00086 !
00087       do j=Jstr,Jend
00088         do i=Istr,Iend
00089           ustar(i,j)=sqrt(sqrt( (0.5*(bustr(i,j)+bustr(i+1,j)))**2
00090      &                         +(0.5*(bvstr(i,j)+bvstr(i,j+1)))**2))
00091         enddo
00092       enddo
00093 !
00094 !  Compute bottom turbulent buoyancy forcing "Bflux" [m^2/s^3]
00095 !  using thermal expansion coefficient "alpha" [kg/m^3/decC] and
00096 !  saline contraction coefficient "beta" [kg/m^3/PSU] at the bottom.
00097 !  Caution: alfabeta is currently only valid for surface coefficients
00098 !
00099 # ifdef BKPP_BFLUX
00100 #  define alpha wrk1
00101 #  define  beta wrk2
00102       call alfabeta_tile (Istr,Iend,Jstr,Jend, alpha,beta)
00103       do j=Jstr,Jend
00104         do i=Istr,Iend
00105           Bflux(i,j)=g*alpha(i,j)*btflx(i,j,itemp)
00106 #  ifdef SALINITY
00107           Bflux(i,j)=Bflux(i,j)-g*beta(i,j)*btflx(i,j,isalt)
00108 #  endif
00109         enddo
00110       enddo  !--> discard alpha,beta; keep Bflux
00111 #  undef beta
00112 #  undef alpha
00113 # endif /* BKPP_BFLUX */
00114 # define sigma wrk1
00115 # define wm wrk3
00116 # define ws wrk2
00117 !
00118 # if !defined BKPP_BFLUX
00119       do j=Jstr,Jend
00120         do i=Istr,Iend
00121           Bflux(i,j)=0.
00122           wm(i,j)=vonKar*ustar(i,j)
00123           ws(i,j)=wm(i,j)          
00124         enddo
00125       enddo
00126 # endif /* BKPP_BFLUX */
00127 !
00128 !----------------------------------------------------------------
00129 !  Compute bulk Richardson number "Rib" and then find depth of the
00130 !  oceanic planetary boundary layer "hbbl", such that Rib(hbbl)=Ric.
00131 !----------------------------------------------------------------
00132 !
00133 !  Set indices for array "Rib", the bulk Richardson number.
00134 !
00135       ka=1 
00136       ku=2
00137 !
00138 !  Intialize boundary layer depth "hbbl" and index "kbl" of first grid
00139 !  level above "hbbl" to maximum values.
00140 !
00141       do j=Jstr,Jend
00142         do i=Istr,Iend
00143           hbbl(i,j)=z_r(i,j,N)-z_w(i,j,0)
00144           kbl(i,j)=N
00145           Rib(i,j,ku)=0.
00146         enddo
00147       enddo
00148 !
00149 !  Find bulk Richardson number at every grid level until > Ric.
00150 !
00151       do k=2,N 
00152 !
00153 !  Compute nondimensional vertical coordinate "sigma".
00154 !  Compute turbulent velocity scales (wm,ws) at "sigma".
00155 !
00156 # ifdef BKPP_BFLUX
00157         do j=Jstr,Jend
00158           do i=Istr,Iend
00159             sigma(i,j)=min(z_r(i,j,k)-z_w(i,j,0),bl_dpth(i,j))
00160           enddo
00161         enddo
00162 
00163         call lmd_wscale_tile (Istr,Iend,Jstr,Jend, Bflux,
00164      &                                       sigma,wm,ws)
00165 # endif /* BKPP_BFLUX */
00166 !
00167 !  Compute bulk Richardson number "Rib" 
00168 !---------------------------------------
00169 !
00170 !                  [Br - B(d)] * d
00171 !     Rib(d) = ----------------------- ;       Rib(hbbl)=Ric
00172 !              |Vr - V(d)|^2 + Vt(d)^2
00173 !
00174 !  To do so, first compute numerator of bluk Richardson number,
00175 !  Ritop=(Br-B)*d, where Br is the near-surface reference buoyancy,
00176 !  B is the mean buoyancy as function of d, and d is the distance
00177 !  coordinate from the boundary.
00178 !
00179 !  Then compute the square of velocity shear relative to reference
00180 !  velocities, dVsq=|Vr-V|^2, at horizontal and vertical RHO-points.
00181 !
00182 !  Then compute Vtsq
00183 !
00184         cff=g/rho0
00185         do j=Jstr,Jend
00186           do i=Istr,Iend
00187             Ritop=-cff*(rho1(i,j,k)-rho1(i,j,1))
00188      &                           *(z_r(i,j,k)-z_r(i,j,1))
00189 
00190             dVsq=0.25*( (u(i  ,j,k,tind)-u(i  ,j,1,tind)+
00191      &                   u(i+1,j,k,tind)-u(i+1,j,1,tind))**2
00192      &                 +(v(i,j  ,k,tind)-v(i,j  ,1,tind)+
00193      &                   v(i,j+1,k,tind)-v(i,j+1,1,tind))**2)
00194 
00195             Vtsq=Vtc*(z_r(i,j,k)-z_r(i,j,1))*ws(i,j)
00196      &          *sqrt(max(0.,0.5*(bvf(i,j,k)+bvf(i,j,k-1))))
00197 
00198             Rib(i,j,ka)=Ritop/(dVsq+Vtsq+eps)
00199           enddo
00200         enddo   !--> discard ws
00201 !
00202 !  Linearly interpolate to find "hbbl" where Rib=Ric.
00203 !
00204         do j=Jstr,Jend
00205           do i=Istr,Iend
00206             if (kbl(i,j).eq.N .and. Rib(i,j,ka).gt.Ric) then
00207               zbl=z_r(i,j,k)-(z_r(i,j,k)-z_r(i,j,k-1))*
00208      &           (Ric-Rib(i,j,ka))/(Rib(i,j,ku)-Rib(i,j,ka))
00209               hbbl(i,j)=zbl-z_w(i,j,0)
00210               kbl(i,j)=k
00211             endif
00212           enddo
00213         enddo
00214         ksave=ka
00215         ka=ku
00216         ku=ksave
00217       enddo    !<-- k    !--> discard Rr,Zr,Ur,Vr,Rib
00218 !
00219 !  Correct "hbbl" with physically limiting case (Ekman depth)
00220 !
00221       do j=Jstr,Jend
00222         do i=Istr,Iend
00223           hekman=lmd_cekman*ustar(i,j)/max(abs(f(i,j)),eps)
00224           hbbl(i,j)=min(hbbl(i,j),hekman)
00225 # ifdef MASKING
00226           hbbl(i,j)=hbbl(i,j)*rmask(i,j)
00227 # endif
00228           kbl(i,j)=N
00229         enddo
00230       enddo
00231 !
00232 !  Find new boundary layer index "kbl".
00233 !
00234       do k=N,1,-1
00235         do j=Jstr,Jend
00236           do i=Istr,Iend
00237             if (z_r(i,j,k)-z_w(i,j,0).gt.hbbl(i,j)) then
00238               kbl(i,j)=k
00239             endif
00240           enddo
00241         enddo
00242       enddo
00243 !
00244 !  Compute tubulent velocity scales (wm,ws) at "hbbl".
00245 !
00246 # ifdef BKPP_BFLUX
00247       do j=Jstr,Jend
00248         do i=Istr,Iend
00249           wm(i,j)=vonKar*ustar(i,j)
00250           ws(i,j)=wm(i,j)
00251         enddo
00252       enddo
00253 # endif
00254 !
00255 !-----------------------------------------------------------------
00256 !  Compute nondimensional shape function Gx(sigma) at "hbbl"
00257 !  (sigma=1) in terms of interior diffusivities (Gx1) and
00258 !  its vertical derivative (dGx1dS) via interpolation.
00259 !-----------------------------------------------------------------
00260 !
00261       do j=Jstr,Jend
00262         do i=Istr,Iend
00263           zbl=z_w(i,j,0)+hbbl(i,j)
00264           k=kbl(i,j)
00265           if (zbl.lt.z_w(i,j,k-1)) k=k-1
00266           cff=1./(z_w(i,j,k)-z_w(i,j,k-1))
00267           cff_up=cff*(zbl-z_w(i,j,k-1))
00268           cff_dn=cff*(z_w(i,j,k)-zbl)
00269 
00270           Kv_bl=cff_up*Kv(i,j,k)+cff_dn*Kv(i,j,k-1)
00271           dKv_bl=-cff*(Kv(i,j,k)-Kv(i,j,k-1))
00272           Gm1(i,j)=Kv_bl/(hbbl(i,j)*wm(i,j)+eps)
00273           dGm1dS(i,j)=min(0.,-dKv_bl/(wm(i,j)+eps))
00274 
00275           Kt_bl=cff_up*Kt(i,j,k)+cff_dn*Kt(i,j,k-1)
00276           dKt_bl=-cff*(Kt(i,j,k)-Kt(i,j,k-1))
00277           Gt1(i,j)=Kt_bl/(hbbl(i,j)*ws(i,j)+eps)
00278           dGt1dS(i,j)=min(0.,-dKt_bl/(ws(i,j)+eps))
00279 
00280 # ifdef SALINITY
00281           Ks_bl=cff_up*Ks(i,j,k)+cff_dn*Ks(i,j,k-1)
00282           dKs_bl=-cff*(Ks(i,j,k)-Ks(i,j,k-1))
00283           Gs1(i,j)=Ks_bl/(hbbl(i,j)*ws(i,j)+eps)
00284           dGs1dS(i,j)=min(0.,-dKs_bl/(ws(i,j)+eps))
00285 # endif /* SALINITY */
00286         enddo
00287       enddo
00288 !
00289 !-----------------------------------------------------------------
00290 !  Compute boundary layer mixing coefficients.
00291 !-----------------------------------------------------------------
00292 !
00293       do k=1,N-1
00294 !
00295 !  Compute turbulent velocity scales at vertical W-points.
00296 !
00297 # ifdef BKPP_BFLUX
00298         do j=Jstr,Jend
00299           do i=Istr,Iend
00300             bl_dnew=hbbl(i,j)*lmd_epsilon
00301             sigma(i,j)=min(z_w(i,j,k)-z_w(i,j,0),bl_dnew)
00302           enddo
00303         enddo
00304         call lmd_wscale_tile (Istr,Iend,Jstr,Jend,
00305      &                                 Bflux,sigma,wm,ws)
00306 # endif
00307 !
00308         do j=Jstr,Jend
00309           do i=Istr,Iend
00310             if (k.lt.kbl(i,j)) then
00311 !
00312 !  Set polynomial coefficients for shape function.
00313 !
00314               sig=min((z_w(i,j,k)-z_w(i,j,0))/(hbbl(i,j)+eps),1.)
00315 # ifdef MASKING
00316               sig=sig*rmask(i,j)
00317 # endif
00318               lmd_a1=sig-2.
00319               lmd_a2=3.-2.*sig
00320               lmd_a3=sig-1.
00321 !
00322 !  Compute nondimensional shape functions.
00323 !
00324               Gm=lmd_a1+lmd_a2*Gm1(i,j)+lmd_a3*dGm1dS(i,j)
00325               Gt=lmd_a1+lmd_a2*Gt1(i,j)+lmd_a3*dGt1dS(i,j)
00326 # ifdef SALINITY
00327               Gs=lmd_a1+lmd_a2*Gs1(i,j)+lmd_a3*dGs1dS(i,j)
00328 # endif
00329 !
00330 !  Compute boundary layer mixing coefficients, combine them
00331 !  with interior mixing coefficients.
00332 !
00333               Kv(i,j,k)=hbbl(i,j)*wm(i,j)*sig*(1.+sig*Gm)
00334               Kt(i,j,k)=hbbl(i,j)*ws(i,j)*sig*(1.+sig*Gt)
00335 # ifdef SALINITY
00336               Ks(i,j,k)=hbbl(i,j)*ws(i,j)*sig*(1.+sig*Gs)
00337 # endif
00338 # ifdef LMD_SKPP
00339 !
00340 ! If reach into SBL, take the max of surface and bottom values 
00341 !
00342               zsbl=z_w(i,j,N)-hbl(i,j)
00343               if (z_w(i,j,k).gt.zsbl) then
00344                 Kv(i,j,k)=max(Kv(i,j,k),Akv(i,j,k))
00345                 Kt(i,j,k)=max(Kt(i,j,k),Akt(i,j,k,itemp))
00346 #  ifdef SALINITY
00347                 Ks(i,j,k)=max(Ks(i,j,k),Akt(i,j,k,isalt))
00348 #  endif
00349               endif
00350 # endif /* LMD_SKPP */
00351 
00352             else  !<-- k > kbl(i,j)
00353 
00354 # ifdef  LMD_CONVEC
00355 !
00356 !  Add convective adjustment
00357 !
00358               if (bvf(i,j,k).lt.0.) then
00359 #  ifdef LMD_SKPP
00360                 zsbl=z_w(i,j,N)-hbl(i,j)
00361                 if (z_w(i,j,k).lt.zsbl) then
00362 #  endif
00363                   Kv(i,j,k)=Kv(i,j,k)+lmd_nu0c
00364                   Kt(i,j,k)=Kt(i,j,k)+lmd_nu0c
00365                   Ks(i,j,k)=Ks(i,j,k)+lmd_nu0c
00366 #  ifdef LMD_SKPP
00367                 endif
00368 #  endif
00369               endif
00370 # endif
00371 
00372             endif  !<-- k > kbl(i,j)
00373           enddo
00374         enddo
00375 # undef ws
00376 # undef wm
00377 # undef sigma
00378 !
00379 !-----------------------------------------------------------------
00380 !  Finalize: Set vertical mixing coefficients "Akv" and "Akt" 
00381 !-----------------------------------------------------------------
00382 !
00383         do j=Jstr,Jend
00384           do i=Istr,Iend
00385             Akv(i,j,k)       = Kv(i,j,k) 
00386             Akt(i,j,k,itemp) = Kt(i,j,k)
00387 # ifdef SALINITY
00388             Akt(i,j,k,isalt) = Ks(i,j,k) 
00389 # endif
00390           enddo
00391         enddo
00392 
00393       enddo ! <-- k
00394 
00395 # if defined EW_PERIODIC || defined NS_PERIODIC || defined MPI
00396       call exchange_w3d_tile (Istr,Iend,Jstr,Jend, Akv)
00397       call exchange_w3d_tile (Istr,Iend,Jstr,Jend,
00398      &                     Akt(START_2D_ARRAY,0,itemp))
00399       call exchange_w3d_tile (Istr,Iend,Jstr,Jend,
00400      &                     Akt(START_2D_ARRAY,0,isalt))
00401 # endif
00402 #else
00403       subroutine lmd_bkpp_empty
00404 #endif /* LMD_BKPP */
00405       return
00406       end
/home/cholod/wd/peps/ROMS/romsagrif_r215/Roms_tools/Roms_Agrif/JUNK/lmd_bkpp_BFLUX.F
