bug in ana_dqdsst.F

[jpringle]

Dear Hernan et al --

I think I have found a bug in ana_dqdsst.F. The stub code from the
most recent model version on subversion reads:

Code: Select all

!
!-----------------------------------------------------------------------
!  Set surface heat flux sensitivity to SST (m/s/degC).
!-----------------------------------------------------------------------
!
#ifdef MY_APPLICATION
      fac=day2sec/30.0_r8            ! 30 day relaxation scale 1/s/decC
      DO j=JstrT,JendT
        DO i=IstrT,IendT
          dqdt(i,j)=fac*Hz(i,j,N(ng))
        END DO
      END DO
#else

This cannot be right, for the units of dqdt are then
(m*seconds/day**2) since the units of Hz are meters and the units of
day2sec are seconds/day, and the units of the 30 are days.

Also, dqdt needs to be negative, or it drives the model
exponentially away from the SST climatology.

A corrected version, also with a minor tweak to the comment, is

Code: Select all

!
!-----------------------------------------------------------------------
!  Set surface heat flux sensitivity to SST (m/s/degC).
!-----------------------------------------------------------------------
!
#ifdef MY_APPLICATION
      fac=1.0_r8/(day2sec*30.0_r8)            ! 30 day relaxation scale 1/s/decC, in the absence of vertical mixing
      DO j=JstrT,JendT
        DO i=IstrT,IendT
          dqdt(i,j)=-fac*Hz(i,j,N(ng))
        END DO
      END DO
#else

I hope this is helpful.

Jamie

[zhang_1988]

The forcing dQdSST is usually computed in units of (Watts/m2/degC). In roms it needs to be scaled by dividing by rho0*Cp .the units of rho0 is kg/m3 and Cp is joule/kg/degcC

so i think the units of dQdSST should be m/s and not m/s/degC

[jpringle]

Dear Zhang --

You are correct. My code is correct, but the comments are wrong. The
comments should be changed to

Code: Select all

!
!-----------------------------------------------------------------------
!  Set surface heat flux sensitivity to SST (m/s).
!-----------------------------------------------------------------------
!
#ifdef MY_APPLICATION
      fac=1.0_r8/(day2sec*30.0_r8)           ! 30 day relaxation scale 1/s, in the absence of vertical mixing
      DO j=JstrT,JendT
        DO i=IstrT,IendT
          dqdt(i,j)=-fac*Hz(i,j,N(ng))
        END DO
      END DO
#else

for others to see that Zhang and the new comments are correct, see the
code below from set_vbc.F

Code: Select all

#  ifdef QCORRECTION
!
!-----------------------------------------------------------------------
!  Add in flux correction to surface net heat flux (degC m/s).
!-----------------------------------------------------------------------
!
! Add in net heat flux correction.
!
      DO j=JstrR,JendR
        DO i=IstrR,IendR
          stflx(i,j,itemp)=stflx(i,j,itemp)+                            &
     &          dqdt(i,j)*(t(i,j,N(ng),nrhs,itemp)-sst(i,j)) &
        END DO
      END DO
 #  endif

You see that dqdt must have units of m/s for stflx to have units of
degC m/s.

Thanks,
Jamie

[wilkin]

Jamie,

You appear to be correct on all counts. I don't know how long the code has been this way.

Arguably, it could be better to leave dqdt > 0 (sensitivity of heat flux to temperature is usually reported as a positive number) and reverse the factor it multiplies in set_vbc.F. This is the way the algorithm is described by Barnier et al. (1995), who estimate dqdt = 50 W/m2/K in most sub-tropical gyres.

B. Barnier, L. Siefridt, P. Marchesiello, Thermal forcing for a global ocean circulation model using a three-year climatology of ECMWF analyses, J. Marine Sys., 1995, 6, 363-380

So in set_vbc.F have instead stflx = stflx + dqdt*(Tobs - Tmodel):

Code: Select all

          stflx(i,j,itemp)=stflx(i,j,itemp)+                            &
     &                     dqdt(i,j)*(sst(i,j-t(i,j,N(ng),nrhs,itemp))

But really, having dqdt > 0 makes no sense.

Consider, say, the sensible heat flux. If Tair > SST we have a sensible heat flux into the ocean (heating: positive to ROMS) so if SST increases the air-sea difference decreases, and sensible heat flux decreases. d(Qsensible)/dt is negative.

The folk I know that use QCORRECTION regularly are ROMS Agrif users, so I checked their code. In the Roms_tools Matlab function that evaluates Barnier's algorithm for dqdt we see:

Code: Select all

%  Compute the kinematic surface net heat flux sensitivity to the
%  the sea surface temperature: dQdSST.
%  Q_model ~ Q + dQdSST * (T_model - SST)
%  dQdSST = - 4 * eps * stef * T^3  - rho_atm * Cp * CH * U
%           - rho_atm * CE * L * U * 2353 * ln (10) * q_s / T^2

So for this to make sense dQdSST must be negative -- if the model is warmer than the obs we decrease the heat flux to cool off the ocean.

Looking in the Roms_Agrif source code (it is deep in analytical.f) the default values for the equatorial demonstration case are indeed negative:

Code: Select all

# if defined EQUATOR
...
    dqdt(i,j)=-50.0/(rho0*Cp)

and there we see Barnier's typical value of 50 W/m2/K (divided by rho*Cp) to put the flux in kinematic units of m/s.

The only thing a bit odd here is that the definition above is independent of the surface grid cell thickness Hz(i,j). Not sure if that's a good thing.

[jpringle]

John --

Thanks. You are right, the inclusion of the surface layer depth in
the formulation is deeply problamatic. However, it is a science
issue not a programming issue.

In this code we prescribe the surface heat flux as a function of
SST. Thus the time rate of change of temperature is set by the
ratio of flux to mix layer depth, not the depth of the first
layer. Deeper mixed layer gives slower time rate of change of SST
for a given heat flux. So as given in the code, the timescale of
SST relaxation is the _minimum_ timescale, for when there is no
vertical mixing and all surface heat flux remains in the top
level. Hence my comment in the code.

(I know you know this John, but I think others might find it
helpful.)

So why not just relax SST directly, so the time rate of change of
temperature depends directly on the SST of the ocean (instead of
the heat flux depending on SST)? Because since the change in SST is
distributed over the entire mixed layer, the effective heat flux in
this case depends on the thickness of the mixed layer. Is this
physical? Not for most terms in the surface heat flux.

But most of this is just a kludge anyway. The only way to do it
right is to have a full atmospheric model bolted on top of the
ocean model with two way coupling. Since for this project I don't
want to do that, I am exploring the least harmful
parameterization. Hence a parameterization in which a missmatch
between the observed and desired SST drives a heat flux.

(actually, I am doing something more complicated -- I am using the
zonal average density mismatch to drive a zonally uniform heat flux
-- if anyone cares, they can see my code.)

Cheers,
Jamie

[wilkin]

So why not just relax SST directly, so the time rate of change of
temperature depends directly on the SST of the ocean (instead of
the heat flux depending on SST)?

Applying a relaxation/nudging term solely to the difference in model and observed temperature and ignoring the "known" heat flux leads to the paradox that if the model is perfect and the temperatures match, then the heat flux is zero.

In a situation where circulation is steadily fluxing heat away - say by vertical mixing out of the bottom of the mixed layer - the overall heat budget will be wrong if the temperature is right. Not good.

Barnier's approach of specifying:
stflux = stflux(observed) + fac*(T_model - SSTobserved)
means that if the model is perfect, then stflux is as observed. So all is good.

I greatly favor the choice of the relaxation fac following Barnier's approach, which recognizes that the equation looks a lot like a one term Taylor series expansion for the heat flux if the model-observed temperature difference is small. So then fac is calculated by differentiating the bulk formulae with respect to temperature and evaluating regionally in concert with the observed heat flux calculation.

This returns a number in W/m2/K which seems to me a more robust approach than an ad hoc choice of time scale.

You make a good point in highlighting a shortcoming in the approach of choosing a single time scale, which is that it leads to a relaxation fac that depends on vertical layer thickness rather than mixed layer depth. Barnier's approach doesn't have that shortcoming, so to answer my own question from yesterday I think it IS a good thing to define dqdt this way:

Code: Select all

dqdt(i,j)=-50.0/(rho0*Cp) ! sensitivity dQdsst = 50 W/m2/K 

and a bad thing to define it this way

Code: Select all

fac=1.0_r8/(day2sec*30.0_r8) ! 30 day relaxation scale 1/s, in the absence of vertical mixing
      ...
          dqdt(i,j)=-fac*Hz(i,j,N(ng))

[jpringle]

Agree.