is LMD constancy preserving

Bug reports, work arounds and fixes

Moderators: arango, robertson

Post Reply
Message
Author
jcwarner
Posts: 1172
Joined: Wed Dec 31, 2003 6:16 pm
Location: USGS, USA

is LMD constancy preserving

#1 Unread post by jcwarner »

I do not typically use LMD mixing, and am not that familiar with the intended operations. Hopefully others will have a more understanding.

In pre_step3d, there is a computation of

Code: Select all

              FC(i,k)=FC(i,k)-dt(ng)*Akt(i,j,k,ltrc)*ghats(i,j,k,ltrc)
However, ghats is computed for salt and temp. The ltrc is limited to min(NAT,itrc). Here is the problem: If you have a zero passive tracer, and use LMD, the pre-step line shown above will add a non-conservative surface flux and this violates the constancy preserving attribute of the model. Perhaps the correct solution is to compute Akt mixing or ghats for each tracer (?).

User avatar
kate
Posts: 4088
Joined: Wed Jul 02, 2003 5:29 pm
Location: CFOS/UAF, USA

Re: is LMD constancy preserving

#2 Unread post by kate »

Hmm, this might explain a problem I was having where non-NAT tracers got outrageous values from the vertical mixing. I "solved" it by limiting Akt for non-NAT tracers, but there might be a better fix.

OK, so this is a real problem and I've submitted a bug report. Now I'm wondering if the portion of the surface flux that's added via LMD_NONLOCAL is added twice - once at the surface and once in the mixed layer.

User avatar
arango
Site Admin
Posts: 1347
Joined: Wed Feb 26, 2003 4:41 pm
Location: DMCS, Rutgers University
Contact:

Re: is LMD constancy preserving

#3 Unread post by arango »

This is a good one :!: The nonlocal transport formulation to vertical mixing is only available and can be computed for the temperature and salinity equations. We need the turbulent buoyancy fluxes to compute this convective enhancement term. We do not know how to do this in other passive tracers like biology and sediment. However, the surface radiative buoyancy flux is affected by the turbidity of water column that can be associated to sediment of biological activity.

We are planning to have a better model for the swdk term which affects the incoming solar radiation and water clarity. John Wilkin is working on this. We will update the code in the future 8)

I think that the safest thing is to include the nonlocal transport in the temperature and salinity equations. Thank you for reporting this problem. I only use the LMD_NONLOCAL option in applications that just have two tracers: temperature and salinity. I will correct the code in the repository. I need to think what it the better way to do this since I need to correct the tangent linear, representer, and adjoint version of this routine.

zxue
Posts: 5
Joined: Tue May 20, 2008 2:59 pm
Location: Louisiana State University

Re: is LMD constancy preserving

#4 Unread post by zxue »

Hi Kate, I am very interested in your way to "limiting Akt for non-NAT tracers" as I am also experiencing some over vertical mixing in biological applications.


kate wrote:Hmm, this might explain a problem I was having where non-NAT tracers got outrageous values from the vertical mixing. I "solved" it by limiting Akt for non-NAT tracers, but there might be a better fix.

OK, so this is a real problem and I've submitted a bug report. Now I'm wondering if the portion of the surface flux that's added via LMD_NONLOCAL is added twice - once at the surface and once in the mixed layer.

User avatar
kate
Posts: 4088
Joined: Wed Jul 02, 2003 5:29 pm
Location: CFOS/UAF, USA

Re: is LMD constancy preserving

#5 Unread post by kate »

You will note that the prior messages were from over a year ago. Since then, pre_step3d has been fixed so that only active tracers get the nonlocal transport.

User avatar
shchepet
Posts: 188
Joined: Fri Nov 14, 2003 4:57 pm

Re: is LMD constancy preserving

#6 Unread post by shchepet »

Just to reply to the original question of John Warner.

Yes, KPP is constancy preserving in its NONLOCAL transport part, as well as everywhere else,
if you do it right.

Unfortunately, the original 1994 LMD paper buries this aspect into a rather sloppy notation making
it very hard to understand what is intended. [I often say that the rationale of KPP is as elegant
as the the United States Federal Tax Code with the part related to mixing within surface boundary
layer playing the role of Alternative Minimum Tax. I truly mean it.]

Sifting through the original codes of Bill Large and Gokhan Danabasoglu (in today's version it is
NCAR CCSM POP) does not make it easier either: these codes look like elaborate money laundering
schemes designed specifically to make a sane Tax Auditor get tired and eventually abandon his quest
to find the truth. The logic of fragmenting everything into many subroutines along with the habit
of multiplying a term by something, then passing it to a subroutine where it is divided by the same
serves neither easier understanding, nor computational efficiency.

Back to LMD_NONLOCAL: The key point is that the nonlocal (also called "countergradient")
flux for a tracer (temperature and salinity) is proportional to the surface forcing flux for that
tracer provided that the net surface buoyancy forcing flux is classified as unstable.


Unstable means that the forcing acts to make stratification in near-surface layer be less stable
(like surface cooling). Nonlocal flux is always zero for all tracers if the net surface buoyancy forcing
flux acts to increase the stability of near surface stratification even if the individual surface forcing
flux for that tracer acts to decrease stability: say, if you have surface cooling decreasing stability,
but at the same time precipitation (rain) which dilutes surface water making it more stratified, and
the net effect is stabilization, then nonlocal flux is zero for both tracers. This makes sense
because the nonlocal flux is meant to simulate extra penetration of surface forcing flux caused
by destabilization (say enhancement of convection near the surface), and once this effect is
absent, the nonlocal flux should disappear for all tracers.


The proportionality to surface forcing flux above means that once the forcing flux is zero,
the nonlocal flux is zero as well. So constancy is preserved. Once there is non-zero surface
forcing flux and the net buoyancy flux acts to destabilize surface layer, the nonlocal flux kicks
in to propagate surface forcing at the rate higher than diffusion (say parcels of heavier water
sink below and mix with ambient water), the constancy is not preserved because it should not.

Now back to LMD 1994 paper. What makes it confusing is that the first time ghat term appears
in Eq. (9) inside the brackets together with the usual tracer gradient and the whole bracket is
multiplied by diffusion coefficient
. To my experience this particular equation is the source of
confusion and code bugs. For example the line from Rutgers ROMS code

Code: Select all

FC(i,k)=FC(i,k)-dt(ng)*Akt(i,j,k,ltrc)*ghats(i,j,k,ltrc)
quoted by John is an example of falling into this patrick-trap.

What you should do instead when reading the LMD 1994 paper is immediately substitute the
expression for mixing coefficient K_x (in LMD notation) in Eq. (9) with its expression from Eq. (10),
where it is defined as the product of hbl*wscale (thickness of boundary layer times friction velocity
adjusted by stratification) multiplied by the non-dimensional shape function "G(sigma)", and, at
the same time immediately substitute the expression for ghat from Eq. (20) [or (19)].
Eq. (20) basically defines "ghat" as the ratio of surface forcing flux divided by hbl*wscale and
the whole thing is multiplied by a non-dimensional constant.

Once you substitute Eqs. (10) and (20) into Eq. (9), the product of hbl*wscale in the what used
to be the "ghat" part of Eq. (9) cancels out, and the outcome boils down to

Code: Select all

   TracerFlux = - hbl * wscale(sigma) * G(sigma) * d_z TracerConcentration + NonLocalFlux
where

Code: Select all

    NonLocalFlux = C * G(sigma) * SurfaceForcingFlux
At this moment I strongly advise to immediately forget the original Eq. (9) and use the above as
the primary definition: reading LMD 1994 between the lines: by Eq. (9) Bill Large made the nonlocal
flux be proportional to turbulent diffusion coefficient. What he actually meant is that the nonlocal
flux within surface boundary layer should be merely proportional to the same shape function G(sigma),
while the hbl*wscale has nothing to do with it.

Does it have to be the same shape function G(sigma) for both the diffusion coefficient and
the nonlocal flux?


The LMD 1994 answer to that is "tentatively yes" in sense that it follows from the formulae [e.g., Eq. (9)]
and nowhere in the paper it is stated otherwise or where possible alternatives are discussed.


My answer to this question is "No", because the shape function for the diffusion coefficient within "hbl"
is designed to make the mixing coefficient be continuous when going from inside to outside of "hbl",
i.e., to satisfy the boundary conditions Eq. (19) at the lower extent of "hbl". This is physically
justifiable for the mixing coefficient, but this also means that G(sigma) has non-zero limit at the lower
extent when sigma --> 1, since the diffusion coefficient below "hbl" is, generally speaking, nonzero.
However, this also means that the NonlocalFlux abruptly changes from non-zero to zero at the lower
extent, which is obviously physically wrong and leads to creating abrupt temperature and salinity
anomalies in the lowest gridbox within "hbl". [Here it should be noted that the NonlocalFlux is technically
a flux, so temperature and salinity tendencies depend on its vertical derivative, so any vertical
discontinuity of NonlocalFlux results in infinitely large temporal increments in temperature and salinity
values when vertical grid resolution becomes very fine.]

So the NonlocalFlux should smoothly vanish at the lower extent of "hbl" no matter what.

The simplest way to achieve this is to set the shape function for NonLocalFlux as

Code: Select all

 G(sigma) = sigma*(1-sigma)^2
which is similar to the shape function for mixing coefficient in the case when there is no mixing
below "hbl".

Here are the relevant pieces from my code:

from file "lmd_kpp.F" where the shape function for the nonlocal flux is computed (note that despite
using the same name,"ghat" computed below is merely a non-dimensional shape function and has
a different meaning than in all other KPP codes):

Code: Select all

      do j=jstr,jend              ! <-- a very long j-loop starts at the beginning "lmd_kpp.F"
....
....
        do i=istr,iend
          do k=N-1,kbl(i),-1
            sigma=(z_w(i,j,N)-z_w(i,j,k))/max(hbl(i,j),eps)
....
....                                  !<-- many other things are computed here as well.
# ifdef LMD_NONLOCAL
            if (Bfsfc .lt. 0.) then
              ghat(i,j,k)=Cg * sigma*(1.-sigma)**2
            else
              ghat(i,j,k)=0.
            endif
# endif
          enddo
          do k=kbl(i)-1,1,-1
# ifdef LMD_NONLOCAL
            ghat(i,j,k)=0.
# endif
....
....
          enddo
        enddo
....
....
      enddo   !<-- j
And the relevant portion of file "step3d_t3S.F" where it is used:

Code: Select all

      do j=jstr,jend              ! <-- a very long j-loop starts at the beginning
...
...
# ifdef LMD_KPP
! Add the solar radiation flux in temperature equation. Also compute
! the nonlocal transport flux for unstable (convective) forcing
! conditions into matrix DC when using the Large et al. 1994 KPP
! scheme.

          if (itrc.eq.itemp) then
            do k=N-1,1,-1
              do i=istr,iend
                cff=srflx(i,j)*swr_frac(i,j,k)
#  ifdef LMD_NONLOCAL
     &                 -ghat(i,j,k)*(stflx(i,j,itemp)-srflx(i,j))
#  endif
                t(i,j,k+1,nnew,itemp)=t(i,j,k+1,nnew,itemp) -dt*cff
                t(i,j,k  ,nnew,itemp)=t(i,j,k  ,nnew,itemp) +dt*cff
              enddo
            enddo
#  if defined LMD_NONLOCAL && defined SALINITY
          elseif (itrc.eq.isalt) then
            do k=N-1,1,-1
              do i=istr,iend
                cff=-dt*ghat(i,j,k)*stflx(i,j,isalt)
                t(i,j,k+1,nnew,isalt)=t(i,j,k+1,nnew,isalt) -cff
                t(i,j,k  ,nnew,isalt)=t(i,j,k  ,nnew,isalt) +cff
              enddo
            enddo
#  endif
          endif
# endif
...
...
      enddo    !<-- j
There are no other instances throughout the entire code where CPP-switch LMD_NONLOCAL is present.

...Hopefully this will help.

User avatar
patrickm
Posts: 26
Joined: Fri Apr 30, 2004 6:41 pm
Location: IRD, FRANCE
Contact:

Re: is LMD constancy preserving

#7 Unread post by patrickm »

Hi Sasha,

In your formulation you use :
NonLocalFlux = C * G(sigma) * SurfaceForcingFlux

taking:
SurfaceForcingFlux = stflx(i,j,itemp) - srflx(i,j)

But in Large et al. (1994) SurfaceForcingFlux also contains the radiative heat absorbed in the boundary layer (at level sigma) as it was believed to contribute to the nonlocal transport (a difference from atmospheric implementations; Troen and Mahrt, 1986). This radiative part of the SurfaceForcingFlux is still in the Rutgers and Agrif codes. Did you remove it for a particular reason? (or is it another patrick trap?) ...

User avatar
shchepet
Posts: 188
Joined: Fri Nov 14, 2003 4:57 pm

Re: is LMD constancy preserving

#8 Unread post by shchepet »

Patrick,

First, the purpose of the above post was to fix what is obviously wrong -- mainly the discontinuity
of nonlocal flux at the base of PBL -- not to advance the physical accuracy of parameterization
(that would be a separate and a much longer subject), and yes, KPP is full of p-traps to fall into.

In my implementation of KPP I generally conform to what Gokhan and Bill Large are doing unless
I feel a compelling reason to deviate -- yes, there are instances where I do, as well as they
themselves keep updating things relative to 1994 paper.

There was a historic bug in ROMS KPP related to GHAT term -- it was wrong even dimensionally.
This was fixed long time ago by Hernan and this fix differs from what I did, an apparently it
differs from what is going in Bill Large and Gokhan' codes. For some period of time (1999 to
~2004) I was simply ignoring this issue by having LMD_NONLOCAL undefined -- it was doing more
harm than good.

Lets examine the issue of having/not having the radiative flux into GHAT term:

Rutgers ROMS (and I believe AGRIF as well) splits GHAT into two, computing it separately for
temperature and salinity. Thus, at about line 320 of lmd_skpp.F from Rutgers ROMS we find:

Code: Select all

# ifdef LMD_NONLOCAL
            cff=1.0_r8-(0.5_r8+SIGN(0.5_r8,Bflux(i,j,k)))
            ghats(i,j,k,itemp)=-cff*(stflx(i,j,itemp)-srflx(i,j)+       &
     &                               srflx(i,j)*(1.0_r8-swdk(i,j)))
#  ifdef SALINITY
            ghats(i,j,k,isalt)=cff*stflx(i,j,isalt)
#  endif
# endif
which is later scaled (essentially multiplied by parameter Cg divided by hbl*wscale)
at about line 890

Code: Select all

# ifdef LMD_NONLOCAL
!  Compute boundary layer nonlocal transport (m/s2).

              cff=lmd_Cg*(1.0_r8-(0.5_r8+SIGN(0.5_r8,Bflux(i,j,k))))/   &
     &            (zbl*ws(i,j)+eps)
              ghats(i,j,k,itemp)=cff*ghats(i,j,k,itemp)
#  ifdef SALINITY
              ghats(i,j,k,isalt)=cff*ghats(i,j,k,isalt)
#  endif
# endif
where the total buoyancy Bflux(i,j,k) flux acts a a simple switch: Bflux(i,j,k)>0 means warming
up the top (stabilizing) causes both ghats(i,j,k,itemp:isalt) --> 0 as it should. The repeated
application of this logic is obviously redundant -- one time, either upper or lower instance is
enough.

Just to avoid any possible ambiguity: in the code above stflx(i,j,itemp) is TOTAL heat flux coming
to the ocean surface from above. This includes everything, latent, evaporation/rain, and radiation
-- incoming short-wave and outgoing infrared. This is historical ROMS convention -- to put the
total heat flux into input files which is sufficient to drive model with a simple surface layer
parameterization (a parameterization which does not care about partitioning the net heat flux
into different physical fluxes). The incoming short-wave radiation srflx(i,j) is available separately
and used when needed. This convention is, however, different from what other models do as
we will see later.

The meaning of

Code: Select all

                               ....(stflx(i,j,itemp)-srflx(i,j)+       &
     &                               srflx(i,j)*(1.0_r8-swdk(i,j)))
same as

Code: Select all

                              stflx(i,j,itemp) - srflx(i,j)*swdk(i,j)
is the total heat absorbed within the portion of water column from the surface to, inclusive,
grid box Hz(i,j,k+1). As usual, in ROMS ghats(i,j,k,itemp) is placed at W-points between Hz(i,j,k)
and Hz(i,j,k+1). The remaining term srflx(i,j)*swdk(i,j) is the radiative flux of the surviving
light which is still not absorbed yet in the water column above. So once the lower extent of PBL
is reached this net absorbed heat becomes the integral net heat absorbed within the boundary
layer at the whole and it is used to compute ghats(:,:,:,itemp). So far all this makes sense.

However, there is a delicate issue associated with using Bflux(i,j,k) as the switch in computation
of ghats: because Bflux(i,j,k) has natural tendency to go from destabilizing to stabilizing as one
tracks it down to through the water column starting from surface and going downward, the logic
in the code above leaves the possibility that ghats is abruptly cut off somewhere in the
middle of boundary layer
. This is physically wrong and leads to noticeable artifacts due to the
discontinuity of the nonlocal flux. In fact, having/not having the non-local flux depends on whether
the boundary layer buoyancy forcing is instable/stable as applied to the boundary layer as whole,
not grid point by grid point.

Lets examine how Bflux(i,j,k) is computed:

At about line 290 of lmd_skpp.F the incoming heat and freshwater fluxes are combined into two
two parts, Bo and Bosol,

Code: Select all

# ifdef SALINITY
          Bo(i,j)=g*(alpha(i,j)*(stflx(i,j,itemp)-srflx(i,j))-          &
     &               beta (i,j)*stflx(i,j,isalt))
# else
          Bo(i,j)=g*alpha(i,j)*(stflx(i,j,itemp)-srflx(i,j))
# endif
          Bosol(i,j)=g*alpha(i,j)*srflx(i,j)
which are later, around line 315, are recombined into Bflux

Code: Select all

        CALL lmd_swfrac_tile (ng, tile,                                 &
     &                        LBi, UBi, LBj, UBj,                       &
     &                        IminS, ImaxS, JminS, JmaxS,               &
     &                        -1.0_r8, zgrid, swdk)
        DO j=Jstr,Jend
          DO i=Istr,Iend
            Bflux(i,j,k)=(Bo(i,j)+Bosol(i,j)*(1.0_r8-swdk(i,j)))
Because swdk(i,j) starts from 1 at surface and naturally decreases as it progresses downward,
and Bosol is always positive, Bflux increases with depth (going downward) and may change its
sign from negative to positive.

Nothing unusual, just another p-trap to fall into.



Now let's examine the "official" KPP:

The current KPP from CCSM 4.0 which is closest to Bill Large and Gokhan Danabasogly is
available at
https://svn-ccsm-release.cgd.ucar.edu/m ... ix_kpp.F90
where you have to register (takes about 5 minutes).

At first, lets see how BO and BOSOL are computed: at about line 1500 we find

Code: Select all

   do j=1,ny_block
   do i=1,nx_block
      if (RHO1(i,j) /= c0) then
         BO   (i,j) = grav*(-TALPHA(i,j)*STF(i,j,1) - &
                             SBETA (i,j)*STF(i,j,2))/RHO1(i,j)

         BOSOL(i,j) = -grav*TALPHA(i,j)*SHF_QSW(i,j)/RHO1(i,j)
      else
         BO   (i,j) = c0
         BOSOL(i,j) = c0
      endif
   end do
   end do
from where it is clear that STF(:,:,1) is the surface forcing heat flux which contains everything
except short-wave radiation. In other words

Code: Select all

       CCSM   STF(:,:,1)    is the same as  ROMS    (stflx(:,:,itemp)-srflx(:,:))
there is nothing wrong with this, just a different convention.

Then BFSFC and stability switch STABLE relevant to the computation of GHATS are computed
within the subroutine bldepth (toward the end of it) between lines 1893 and 1937, where, the
depth used for the purpose of computing absorbed short-wave fraction is the depth of
HBL
(pay attention to the arguments for the corresponding light attenuation routines):

Code: Select all

   if (lshort_wave) then
      select case (sw_absorption_type)

      case ('top-layer')

        BFSFC   = BO + BOSOL
!       QSW_HBL = SHF_QSW
        if (tavg_requested(tavg_QSW_HBL)) then
           WORK = SHF_QSW/hflux_factor
           call accumulate_tavg_field(WORK,tavg_QSW_HBL,bid,1)
        endif


      case ('jerlov')

         do j = 1,ny_block
         do i = 1,nx_block
            call sw_absorb_frac(HBLT(i,j),absorb_frac)
            BFSFC(i,j)  = BO(i,j) + BOSOL(i,j)*(c1 - absorb_frac)
         enddo
         enddo

         if (tavg_requested(tavg_QSW_HBL)) then
           !QSW_HBL(i,j) = SHF_QSW(i,j)*(c1-absorb_frac)  ! boundary layer sw
            WORK = SHF_QSW*(c1-absorb_frac)/hflux_factor
            call accumulate_tavg_field(WORK,tavg_QSW_HBL,bid,1)
         endif

      case ('chlorophyll')

         ZTRANS(:,:,bid) = HBLT(:,:)
         call sw_trans_chl(0,this_block)
         BFSFC   = BO + BOSOL*(c1-TRANS(:,:,bid))

         if (tavg_requested(tavg_QSW_HBL)) then
           !QSW_HBL = SHF_QSW   *(c1-TRANS(:,:,bid)) ! boundary layer sw heating
            WORK = SHF_QSW*(c1-TRANS(:,:,bid))/hflux_factor
            call accumulate_tavg_field(WORK,tavg_QSW_HBL,bid,1)
         endif

      end select
   endif

   STABLE = MERGE(c1, c0, BFSFC >= c0)

   BFSFC  = BFSFC + STABLE * eps ! ensures bfsfc never=0
where, of course, this is an array-syntax code, so one should remember that both BFSFC and
STABLE are two-dimensional arrays. Obviously, this part is different than what we have in ROMS.

Then goes the computation of GHAT. Note that there is only single GHAT, i.e., there is no
such thing like separate GHATs for temperature and salinity. At about line 2224 we find the
preliminary computation

Code: Select all

   do k = 1,km
      if (partial_bottom_cells) then
         if (k > 1) then
            SIGMA = (-zgrid(k-1) + p5*DZT(:,:,k-1,bid) +  &
                     DZT(:,:,k,bid)) / HBLT
         else
            SIGMA = (-zgrid(k) + p5*hwide(k)) / HBLT
         end if
      else
         SIGMA = (-zgrid(k) + p5*hwide(k)) / HBLT
      endif
      F1 = min(SIGMA,epssfc)

      call wscale(F1, HBLT, USTAR, BFSFC, 3, WM, WS)

      !DIR$ COLLAPSE
      do j=1,ny_block
      do i=1,nx_block
         BLMC(i,j,k,1) = HBLT(i,j)*WM(i,j)*SIGMA(i,j)*       &
                         (c1 + SIGMA(i,j)*((SIGMA(i,j)-c2) + &
                         (c3-c2*SIGMA(i,j))*GAT1(i,j,1) +    &
                         (SIGMA(i,j)-c1)*DAT1(i,j,1)))
         BLMC(i,j,k,2) = HBLT(i,j)*WS(i,j)*SIGMA(i,j)*       &
                         (c1 + SIGMA(i,j)*((SIGMA(i,j)-c2) + &
                         (c3-c2*SIGMA(i,j))*GAT1(i,j,2) +    &
                         (SIGMA(i,j)-c1)*DAT1(i,j,2)))
         BLMC(i,j,k,3) = HBLT(i,j)*WS(i,j)*SIGMA(i,j)*       &
                         (c1 + SIGMA(i,j)*((SIGMA(i,j)-c2) + &
                         (c3-c2*SIGMA(i,j))*GAT1(i,j,3) +    &
                         (SIGMA(i,j)-c1)*DAT1(i,j,3)))

         GHAT(i,j,k) = (c1-STABLE(i,j))* cg/(WS(i,j)*HBLT(i,j) +eps)
      end do
      end do
   end do
then apply a fix to GHAT inside the last grid box which is at the transition from PBL to below it
[This is known an "enhance" routine designed to smooth out the vertical movement of boundary
layer boundary as it progresses along the discrete grid; CASEA is a switch which is 0 or 1
depending on whether the lower extent of PBL fall into upper or lower half of the grid box; not
principle for the present discussion; there is no counterpart to it in ROMS codes].
At about line 2332 we find:

Code: Select all

         GHAT(:,:,k) = (c1-CASEA) * GHAT(:,:,k)
And finally the usage of GHAT(:,:,k) at around line 910

Code: Select all

   do n=1,nt
      mt2=min(n,2)
      KPP_SRC(:,:,1,n,bid) = STF(:,:,n)/dz(1)           &
                             *(-VDC(:,:,1,mt2)*GHAT(:,:,1))
      if (partial_bottom_cells) then
         do k=2,km
            KPP_SRC(:,:,k,n,bid) = STF(:,:,n)/DZT(:,:,k,bid)         &
                                 *( VDC(:,:,k-1,mt2)*GHAT(:,:,k-1)   &
                                   -VDC(:,:,k  ,mt2)*GHAT(:,:,k  ))
         enddo
      else
         do k=2,km
            KPP_SRC(:,:,k,n,bid) = STF(:,:,n)/dz(k)                  &
                                 *( VDC(:,:,k-1,mt2)*GHAT(:,:,k-1)   &
                                   -VDC(:,:,k  ,mt2)*GHAT(:,:,k  ))
         enddo
      endif
   enddo
where VDC(:,:,k-1,mt2) is vertical diffusivity for tracer "mt2" saved from
the previous time step.


From the examination of KPP code from CCSM 4.0 above we can draw the following four conclusions:

1. In CCSM 4.0 the nonlocal flux is tied to the surface forcing flux STF, which, in the case of
temperature EXCLUDES the incoming short-wave radiation, since it is the same STF as a used
in computation of buoyancy forcing.

Code: Select all

      BO    = grav*(-TALPHA*STF(:,:,1) - SBETA*STF(:,:,2))/RHO1

      BOSOL = -grav*TALPHA*SHF_QSW/RHO1

      BFSFC(i,j) = BO(i,j) + BOSOL(i,j)*(c1 - absorb_frac)
Hopefully, this answers Patrick's question without leaving any ambiguity: and yes, what is in the
"official" KPP is different from what is implemented in Rutgers and AGRIF ROMS with respect to
this feature.

I am not saying that one of them is correct and the other is wrong. All I am saying is that they are
different and I do not like them both in sense that there are some merits there, but both of them
are incomplete. Obviously, if we imagine very transparent water so light goes straight through the
boundary layer both versions give the same result (in the case of Rutgers code above swdk(:,:)=1
in this case which effectively results in straight subtraction of srflx from stflx), and intuitively
this is result is correct. Conversely, infinitely fast absorption of light results in no distinction
between latent and short-wave heat flux, so it seems that swrad should be fully applied upper grid
box, then Rutgers ROMS version has more merit in this case. Somewhere between the two limits it
is hard to physical justify either of these approaches.

2. In CCSM 4.0 when computing GHATS the stability switch STABLE(:,:) (a two-dimensional array)
is computed using depth of HBL. This is correct and this is how it is intended to be.

In contrast Rutgers ROMS when computing ghats, the stability switch is effectively set individually
for each vertical grid point. This is wrong and needs to be corrected.

3. In both CCSM 4.0 and Rutgers ROMS nonlocal flux is discontinuous at the base of HBL (unless
vertical mixing coefficient, VDC and Akt in the respective codes, smoothly vanish, which means
that the mixing below HBL is identically zero). This is physically wrong and needs to be corrected.

4. Peculiarly, and this perhaps answers to the questions of John and Kate, CCSM 4.0 applies
GHAT not just to temperature and salinity, but to all tracers, which possible include biological
tracers as well, as long as there is non-zero surface forcing flux for them (say precipitation of
nutrients or acid rain or something).

User avatar
shchepet
Posts: 188
Joined: Fri Nov 14, 2003 4:57 pm

Re: is LMD constancy preserving

#9 Unread post by shchepet »

I had a chance to examine the relevant KPP routines from a recent ROMS AGRIF,
files lmd_skpp.F and step3d_t.F of the following SVN versions
! $Id: lmd_skpp.F 726 2011-08-12 13:48:44Z gcambon $
! $Id: step3d_t.F 697 2011-04-11 12:35:17Z gcambon $
and found the following:


ghats(i,j,k), the same for all tracers is computed inside lmd_skpp.F at about line 570

Code: Select all

# ifdef LMD_NONLOCAL
!
!  Compute boundary layer nonlocal transport [m/s^2]
!
              if (Bflux(i,j).le.0.) then
                ghats(i,j,k)=lmd_Cg/(ws(i,j)*hbl(i,j)+eps)
              else
                ghats(i,j,k)=0.
              endif
# endif
where Bflux(i,j) is computed at the lower extent of PBL, and ws(i,j) is computed at whatever
depth of grid point z_w(i,j,k).

This conforms to what is going on in CCSM 4.0 code and all other earlier versions from Bill Large
and Gokhan Danabasoglu I am aware of.

Then ghats(i,j,k) appears inside step3d_t.F where it is used as follows:

Code: Select all

! Add top and bottom fluxes
!
          do i=Istr,Iend
            FC(i,N)=dt*stflx(i,j,itrc)
            FC(i,0)=-dt*btflx(i,j,itrc)
          enddo
!
! Add solar radiation flux in temperature equation
! Also compute the nonlocal transport flux for unstable
! (convective) forcing conditions into matrix DC when using
! the Large et al. 1994 KPP scheme.
!
          if (itrc.eq.itemp) then
            do k=1,N-1
              do i=Istr,Iend
                FC(i,k)=0.
# if defined LMD_SKPP || defined LMD_BKPP
     &                 +dt*srflx(i,j)*swdk(i,j,k)
#  ifdef LMD_NONLOCAL
     &                 -dt*Akt(i,j,k,itemp)*ghats(i,j,k)
#  endif
# endif
              enddo
            enddo
# ifdef SALINITY
          elseif (itrc.eq.isalt) then
            do k=1,N-1
              do i=Istr,Iend
                FC(i,k)=0.
# if defined LMD_SKPP || defined LMD_BKPP
#  ifdef LMD_NONLOCAL
     &                 -dt*Akt(i,j,k,isalt)*ghats(i,j,k)
#  endif
# endif
              enddo
            enddo
# endif
          endif
# ifdef SALINITY
          if (itrc.eq.itemp .or. itrc.eq.isalt) then
# else
          if (itrc.eq.itemp) then
# endif
            do k=1,N
              do i=Istr,Iend
                t(i,j,k,nnew,itrc)=t(i,j,k,nnew,itrc)+FC(i,k )
     &                                               -FC(i,k-1)
              enddo
            enddo
          endif
This is wrong and, obviously, even dimensionally wrong, because ghats contribution
to FC for temperature and FC for salinity should have different dimensions, but they clearly
have the same dimension, c.f.,

Code: Select all

          FC(i,k) = ...           -dt*Akt(i,j,k,itemp)*ghats(i,j,k)

and

Code: Select all

           FC(i,k) = ...          -dt*Akt(i,j,k,isalt)*ghats(i,j,k)

respectively.

To fix the above, and in fact to make it conform to Gokhan's and Bill Large code, the above
ghats contributions to FC fluxes should be multiplied by their respective surface forcing fluxes,

Code: Select all

....
....
                FC(i,k)=0.
# if defined LMD_SKPP || defined LMD_BKPP
     &                 +dt*srflx(i,j)*swdk(i,j,k)
#  ifdef LMD_NONLOCAL
     &                 -dt*Akt(i,j,k,itemp)*ghats(i,j,k)*(stflx(i,j,itemp)-srflx(i,j))
#  endif
# endif
....
....
for temperature, and

Code: Select all

....
                FC(i,k)=0.
# if defined LMD_SKPP || defined LMD_BKPP
#  ifdef LMD_NONLOCAL
     &                 -dt*Akt(i,j,k,isalt)*ghats(i,j,k)*stflx(i,j,isalt)
#  endif
# endif
....
for salinity.

[Whether in the case of temperature it should be (stflx(i,j,itemp)-srflx(i,j)) or some other
combination of stflx and srflx as suggested by Patrick remains open for the discussion,
but the fact that these FC's are missing dimensional multipliers leaves no ambiguity.]

But then again, the remark 3 at the end of my previous post (the possibility of
discontinuity of the non-local flax at the local extent of PBL) stands.

The products Akt(i,j,k,itemp)*ghats(i,j,k) and Akt(i,j,k,isalt)*ghats(i,j,k) are merely
non-dimensional shape functions because Akt is basically

Code: Select all

   Akt = wscale * hbl * sigma * (1-sigma)^2
in the case where there is no mixing below PBL, or

Code: Select all

   Akt = wscale * hbl * G(sigma)
where G(sigma) is a more complicated shape function in a more general case, but still asymptotes
as G(sigma) ~ sigma, when sigma --> 0, and, on the other hand, ghats is basically

Code: Select all

 ghats = C_g / (wscale * hbl)
so if you cancel out wscale * hbl you will wind up with Akt*ghat = G_g*G(sigma) in basically naked
form, and if you replace G(sigma) from the ones used to compute Kt and Ks with the one which
smoothly vanishes sigma=1, you will end up with what is already in my code.
Last edited by shchepet on Wed Jan 25, 2012 4:52 pm, edited 3 times in total.

User avatar
patrickm
Posts: 26
Joined: Fri Apr 30, 2004 6:41 pm
Location: IRD, FRANCE
Contact:

Re: is LMD constancy preserving

#10 Unread post by patrickm »

Hi Sasha,

Thanks for the details on the radiative heating choices. I understand the issue of Bflux changing sign with depth and how it may affect the physical consistency of nonlocal mixing. This is a good point that leaves us with interesting physical issues.

Note that we are well aware of the "bug" in ROMS_AGRIF concerning the NONLOCAL option. As you did yourself before rewritting KPP, we used LMD_KPP with the NONLOCAL option turned off, waiting for a better understanding of how it should be implemented. This choice dates back to the beginning of ROMS (1998). It is certainly time to update the code at present. Actually, we are heading towards a full upgrade to your new version of KPP...

Patrick

Post Reply