Processing ECMWF ERA forcing fields for ROMS

[arango]

In ROMS, we have the option to compute the surface forcing fluxes using BULK_FLUXES or providing directly the fluxes (surface wind stress sustr, svstr, surface net heat flux shflux, and surface freshwater swflux) from an atmospheric model dataset. In some applications with open boundary conditions like T/S radiation, the recommendation is to not use BULK_FLUXES.

The explanation is kind of complicated. The bulk flux formulation in ROMS uses the Monin-Obukhov similarity parameters of Liu et al. (1979) to compute stability functions that compute the turbulent fluxes for wind (Wstar), heat (Tstar), and moisture (Qstar). There are stable and unstable regimes for these functions. This computation can be highly nonlinear. Any bias or errors in the open boundary conditions for temperature, due to radiation conditions, may result in a loss or gain of heat at the boundary. This may create bogus upwelling/downwelling at the open boundary edges. The model becomes unstable quite rapidly and it blows up.

Therefore, I avoid to use BULK_FLUXES in applications with dynamically active circulation regimes. I get instead the total net heat flux and wind stress from a dataset or interpolate them from the coarse model. I usually use the fluxes from ECMWF's ERA-Iterim Dataset. This is one of my favorite datasets. I think that the net incoming shortwave radiation at the ocean surface that we get from other datasets is too high resulting in excessive heat flux in the ocean.

I usually get the following variables from the ERA dataset: The @ denotes accumulated quantity that must be divided by the time interval into the cycle 3, 6, 9 or 12 hours:

Code: Select all

 Select time:   00:00:00     12:00:00
 Select step:   0  3  6  9  12

 @  sshf    W m-2 s        surface sensible heat flux
 @  slhf    W m-2 s        surface latent heat flux
 @  ssr     W m-2 s        surface net solar radiation (shortwave)
 @  str     W m-2 s        surface net thermal radiation (longwave)
 @  strd    W m-2 s        surface thermal radiation downwards
 @  ewss    N m-2 s        east-west surface stress
 @  nsss    N m-2 s        north-south surface stress
 @  e       m              evaporation (downward flux is positive)
 @  ro      m              runoff
 @  tcc     nondimensional total cloud cover [0:1]
 @  tp      m              total precipitation
 @  par     W m-2 s        photosynthetically active radiation at surface
    msl     Pa             mean sea level pressure
    v10v    m s-1          10 metre U wind component
    vl0u    m s-1          10 metre V wind component
    v2t     K              2 metre temperature
    v2d     K              2 metre dewpoint temperature
 
  This dataset is written in compact way (short numbers). We need to convert to floating-point data and scale to ROMS units:

    Uwind       (m s-1)         v10u
    Vwind       (m s-1)         v10v
    sustr       (N m-2)         ewss / (3*3600);   3-hour step
    svstr       (N m-2)         nsss / (3*3600);   3-hour step 
    shflux      (W m-2)         (ssr+str+sshf+slhf) / (3*3000);   3-hour step
    swrad       (W m-2)         ssr  / (3*3600);   3-hour step
    lwrad_down  (W m-2)         strd / (3*3600);   3-hour step
    latent      (W m-2)         slhf / (3*3600);   3-hour step
    sensible    (W m-2)         sshf / (3*3600):   3-hour step
    rain        (kg m-2 s-1)    tp * Rho_w / (3*3600)
    evaporation (kg m-2 s-1)    e  * Rho_w / (3*3600)
    swflux      (cm day-1)      (-e - tp) * 100 / (3/24);  0.125 day step
    cloud       (nondimesional) tcc
    Pair        (mb)            msl / 100;   (1 mb = 100 Pa)
    Tair        (Celsius)       t2m - 273.15;   (1 C = 273.15 K)
    Qair        (percentage)    100 * (E/Es)
 
  where

    Rho_w = 1000 kg m-3  (water density)

    E  = 6.11 * 10.0 ^ (7.5 * v2d / (237.7 + v2d))    vapor pressure (mb)
                                                      v2d in Celsius
 
    Es = 6.11 * 10.0 ^ (7.5 * v2t / (237.7 + v2t))    saturation vapor
                                                      pressure (mb)
                                                      v2t in Celsius

If we denote V3, V6, V9, and V12 as the 3, 6, 9, and 12 hour accumulated values, respectively, the forcing fields averages (A3, A6, A9, and A12) for the three-hour interval are:

Code: Select all

   A3  = V3 / (3*3600)
   A6  = (V6 - V3) / (3*3600)
   A9  = (V9 - V6) / (3*3600)
   A12 = (V12 - V9) / (3*3600)

You can check out this by comparing the processed instantaneous and accumulated fields for surface momentum stress (iews vs ewss; inss vs nsss).



Notice that I usually get the variables for BULK_FLUXES in case that I want to compare solutions with or without it. If we are not using BULK_FLUXES, ROMS just needs shflux, swflux, swrad (with some CPP options), sustr, and svstr.

I added the Matlab script template forcing/d_ecmwf2roms.m to the matlab repository to guide you how to process these forcing fields.