Shortwave radiation and net surface heat flux inconsistency

[lcbernardo]

Dear ROMS users,

I've encountered some confusion simulating water temperatures for my site. I have been running both Rutgers ROMS (circa October 2018) and ROMS-only COAWST (version 3.3) and there was a striking difference in the results (more than 1 deg C temperature difference), even when using the similar cpp definitions and the same input files. The Rutgers ROMS results are actually closer to the observations, and I'm trying to understand why there's an inconsistency. I tried outputting the met parameters like incoming shortwave radiation and net surface heat flux for a specific grid cell, and the values I get for the ROMS simulation are twice or more that as the COAWST result (see attached figure).

I think that this issue is related to BULK_FLUXES and ANA_SRFLUX, and for each configuration, the CPP defs related to the surface heat flux calculations are as follows:

(ROMS)
#define SOLAR_SOURCE
#define ATM_PRESS
#define BULK_FLUXES
#ifdef BULK_FLUXES
# define LONGWAVE
# define EMINUSP
# define COOL_SKIN
# define ANA_SRFLUX
# define ALBEDO
#else
# define ANA_SMFLUX
# define ANA_STFLUX
#endif

(COAWST)
#define SOLAR_SOURCE
#define ATM_PRESS
#define BULK_FLUXES
#ifdef BULK_FLUXES
# define LONGWAVE
# define EMINUSP
# define ANA_SRFLUX
# define COOL_SKIN
# define ALBEDO_CLOUD
#else
# define ANA_SMTFLUX
# define ANA_STFLUX
#endif

The only difference is ALBEDO is ALBEDO_CLOUD in COAWST, which I believe to be the correct equivalent after lookin through ana_srflux.f. As for my input met data for the model, they are obtained from an atmospheric model, and the parameters provided are the wind components, air pressure, air temperature, relative humidity, rainfall rate, and cloud fraction. And I am using the same input files for both ROMS and COAWST runs.

Can anyone help me find the possible source or sources of this inconsistency? Thanks!

Lawrence

[wilkin]

COAWST handles the humidity slightly differently (in ana_srflux.h):

Code: Select all

#  ifdef SPECIFIC_HUMIDITY
!  With this directive specific humidity is input as kg/kg
            vap_p=Pair(i,j)*Hair(i,j)/(0.62197_r8+0.378_r8*Hair(i,j))
#  else
            vap_p=e_sat*Hair(i,j) ! water vapor pressure (hPa=mbar)
#  endif

and in bulk_flux.F:

Code: Select all

# if defined LONGWAVE
!
!  Use Berliand (1952) formula to calculate net longwave radiation.
!  The equation for saturation vapor pressure is from Gill (Atmosphere-
!  Ocean Dynamics, pp 606). Here the coefficient in the cloud term
!  is assumed constant, but it is a function of latitude varying from
!  1.0 at poles to 0.5 at the Equator).
!
#  ifdef SPECIFIC_HUMIDITY
!  specific humidity in units of kg/kg
          vap_p=PairM*RH/(1.0_r8+0.378_r8*RH)             !WPB
#  else
      IF(RH.lt.2.0_r8) THEN                               !WPB
          cff=(0.7859_r8+0.03477_r8*TairC(i))/                          &
     &        (1.0_r8+0.00412_r8*TairC(i))
          e_sat=10.0_r8**cff   ! saturation vapor pressure (hPa or mbar)
          vap_p=e_sat*RH       ! water vapor pressure (hPa or mbar)
      ELSE                                                !WPB
          vap_p=0.001_r8*PairM*RH/(1.0_r8+0.000378_r8*RH) !WPB
      ENDIF                                               !WPB
#  endif

I forget the history of this, but the myroms.org code was written to trap inconsistent humidity units with the test

Code: Select all

IF(RH.lt.2.0_r8) THEN 

which is a bit of a kludge. You might look into this as a source of the divergent swrad calculations.

[lcbernardo]

Thank you Dr. Wilkin for directing me to these sections in the code! I think I can see now how the inconsistency arises, though I'm also curious as to the reason/basis for the kludge. As you mentioned, it's something for me to look into further and to tweak.

Best,
Lawrence

[kate]

The code I was first given had the test on units for humidity cited above. A colleague of mine didn't like how that test is done and wanted to use only kg/kg as the specific humidity units. The forcing files we've been using for the past decade + have all been specific humidity in kg/kg, so nothing else has been tested in my code. Other choices may fail.

[lcbernardo]

Hi Kate, thanks for sharing a bit of history on the code. I think indeed other choices may have some problems.

As a follow up to this issue, I think I found the source of the inconsistency and it was in the Nonlinear/bulk_flux.F of the COAWST version I'm using, in particular the lines:

# if defined ALBEDO && defined SHORTWAVE
# if defined ALBEDO_CLOUD && !defined BENCHMARK
cawdif = 1._r8-albw_d
srflx(i,j) = srflx(i,j)*(cawdir(i,j)* &
& (1._r8-cloud(i,j))+cawdif*cloud(i,j))
# else
srflx(i,j) = srflx(i,j)*(1._r8-albedo(i,j))
# endif
# endif

It seems that cawdir wasn't being computed correctly and defaulted to zero, and srflx ended up being a function of cawdif*cloud. This seems to have resulted in the reduced magnitudes compared to ROMS, which does not have similar code.

As to why cawdir all had zero values, it turns out that some years back we tried to turn off all references in the code to cawdir_eval, which computes cawdir, because the code wouldn't compile. I tried restoring these references, but always end up having the error message below upon compiling:

cawdir_eval.f90(184): error #6404: This name does not have a type, and must have an explicit type. [R_DATE]
call caldate (r_date,tdays(ng),year,yday,month,iday,hour)
--------------------^
compilation aborted for cawdir_eval.f90 (code 1)

Perhaps anyone who was used Utility/cawdir_eval.F previously can share some thoughts? Looking at newer versions of the COAWST code however, I can't seem to find this file anymore, and I'm trying still trying to figure out how it is being implemented.

[kate]

Hernan rewrote the caldate routine at some point. All instances of calling it with r_date come from my moldy code (grep -r r_date ROMS). The new style should look like:

Code: Select all

ROMS/Functionals/ana_specir.h:      USE dateclock_mod,   ONLY : caldate
ROMS/Functionals/ana_specir.h:      CALL caldate (tdays(ng), yd_r8=yday, h_r8=hour)

[lcbernardo]

Thank you so much Kate! Indeed, adjusting to the new style was able to solve the issue. Now cawdir_eval.F seems to be working as it was designed to do. I greatly appreciate everyone's kind assistance on this issue!

Best,
Lawrence