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)is LMD constancy preserving
is LMD constancy preserving
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
Re: is LMD constancy preserving
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: 1367
- Joined: Wed Feb 26, 2003 4:41 pm
- Location: DMCS, Rutgers University
- Contact:
Re: is LMD constancy preserving
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.
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.
Re: is LMD constancy preserving
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.
Re: is LMD constancy preserving
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.
Re: is LMD constancy preserving
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: 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)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
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 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
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.
Re: is LMD constancy preserving
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?) ...
Re: is LMD constancy preserving
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 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
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, 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)))
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 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)
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 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
except short-wave radiation. In other words
Code: Select all
CCSM STF(:,:,1) is the same as ROMS (stflx(:,:,itemp)-srflx(:,:))
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
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
[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 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)
"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).
Re: is LMD constancy preserving
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)
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
....
....
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 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
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.
Re: is LMD constancy preserving
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