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)`

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
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 (?).

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)`

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.

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.

- arango
- Site Admin
**Posts:**1108**Joined:**Wed Feb 26, 2003 4:41 pm**Location:**IMCS, Rutgers University-
**Contact:**

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

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.

We are planning to have a better model for the

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

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.

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.

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
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
where
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
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):
And the relevant portion of file "step3d_t3S.F" where it is used:
There are no other instances throughout the entire code where CPP-switch LMD_NONLOCAL is present.

...Hopefully this will help.

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:

flux for a tracer (temperature and salinity) is

tracer provided that the net surface buoyancy forcing flux is classified as

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

flux acts to increase the stability of near surface stratification

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

the nonlocal flux is zero as well.

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

Now back to LMD 1994 paper. What makes it confusing is that the first time

in Eq. (9) inside the brackets together with the usual tracer gradient

multiplied by diffusion coefficient

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)`

What you should do instead when reading the LMD 1994 paper is

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

adjusted by stratification) multiplied by the non-dimensional shape function "G(sigma)", and, at

the same time

Eq. (20) basically defines "ghat" as the ratio of surface forcing flux divided by

the whole thing is multiplied by a non-dimensional constant.

Once you substitute Eqs. (10) and (20) into Eq. (9), the product of

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
```

Code: Select all

```
NonLocalFlux = C * G(sigma) * SurfaceForcingFlux
```

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

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

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
```

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
```

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
```

...Hopefully this will help.

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?) ...

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?) ...

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:
which is later scaled (essentially multiplied by parameter *Cg* divided by *hbl*wscale*)

at about line 890
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
same as
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,
which are later, around line 315, are recombined into Bflux
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
from where it is clear that STF(:,:,1) is the surface forcing heat flux which contains everything

**except short-wave radiation**. In other words
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):
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
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:
And finally the usage of GHAT(:,:,k) at around line 910
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.
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).

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

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
```

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
```

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,

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)))
```

Code: Select all

```
stflx(i,j,itemp) - srflx(i,j)*swdk(i,j)
```

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

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

middle of boundary layer

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

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)
```

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)))
```

and Bosol is always positive, Bflux increases with depth (going downward) and may change its

sign from negative to positive.

Nothing unusual, just another

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
```

Code: Select all

```
CCSM STF(:,:,1) is the same as ROMS (stflx(:,:,itemp)-srflx(:,:))
```

Then BFSFC and stability switch STABLE

within the subroutine bldepth (toward the end of it) between lines 1893 and 1937, where,

depth used for the purpose of computing absorbed short-wave fraction is the depth of

HBL

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
```

STABLE are two-dimensional arrays. Obviously, this part is different than what we have in ROMS.

Then goes the computation of GHAT.

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
```

[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)
```

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
```

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

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)
```

"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

tracers as well, as long as there is non-zero surface forcing flux for them (say precipitation of

nutrients or acid rain or something).

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
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:
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.,

and
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,
for temperature, and
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
in the case where there is no mixing below PBL, or
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
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.

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
```

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
```

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)
```

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
....
....
```

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
....
```

[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

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
```

Code: Select all

```
Akt = wscale * hbl * G(sigma)
```

as G(sigma) ~ sigma, when sigma --> 0, and, on the other hand, ghats is basically

Code: Select all

```
ghats = C_g / (wscale * hbl)
```

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.

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

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