Problem with surface heat flux

[cloudy4568]

Hi all,

Thank you for your help in advance. There are some problems with my result of water temperature.

My roms result of water temperature is much higher than data from ECMWF or HYCOM from attached figure 1. And I think there is something wrong with heat flux so I compared each component in heat flux in attached figure 2. The problem is, apparently, the sensible heat flux is almost 0 all the time and latent heat flux is positive (latent heat should be negative).

The matlab code used to generate input netcdf is:

Code: Select all

% D_ECMWF2ROMS:  Driver script to create ROMS forcing NetCDF from ECMWF
%                ERA-Interim Dataset
%
% This is a user modifiable script showing how to create ROMS forcing
% NetCDF files. The data source is the ECMWF's ERA-Interim Dataset

Atmosphere data is from matlab code's recommendation:

Code: Select all

% http://apps.ecmwf.int/datasets/data/interim_full_daily/
%
% This global data is available from Jan 1, 1979 to the present. It is
% usually extracted to a regional grid from the ECMWF data server.  No
% attempt is made to interpolate the forcing fields to a particular
% ROMS grid.  ROMS has the capability to perform the spatial and
% temporal interpolation of 2D forcing fields.
%
% This is a template script showing how to create several NetCDF files
% (one per each atmospheric forcing field) using ROMS metadata structure,
% which follows the schema of "nc_inq" or native 'ncinfo" function.
%
% The following parameters are used to extract ERA-Interim fields:
%
% Select time:   00:00:00     12:00:00
%
% Select step:   0  3  6  9  12

Data includes:

Code: Select all

%   Uwind       (m s-1)         v10u
%   Vwind       (m s-1)         v10v
%   swrad       (W m-2)         ssr  / (3*3600);   3-hour step
%   lwrad_down  (W m-2)         strd / (3*3600);   3-hour step
%   rain        (kg m-2 s-1)    tp * Rho_w / (3*3600)
%   Pair        (mb)            msl / 100;   (1 mb = 100 Pa)
%   Tair        (Celsius)       t2m - 273.15;   (1 C = 273.15 K)
%   Qair        (percentage)    100 * (E/Es)

CCSM_FLUXES is used.
Activated C-preprocessing Options:

Code: Select all

 ADD_FSOBC               Adding tidal elevation to processed OBC data
 ADD_M2OBC               Adding tidal currents to processed OBC data
 ALBEDO_CURVE            Ocean shortwave radiation adjusted by albedo curve
 ANA_BSFLUX              Analytical kinematic bottom salinity flux
 ANA_BTFLUX              Analytical kinematic bottom temperature flux
 ANA_DRAG_GRID           Analytical spatially varying quadratic drag coefficient
 ANA_ICE                 Analytical ice initial conditions
 ASSUMED_SHAPE           Using assumed-shape arrays
 AVERAGES                Writing out time-averaged nonlinear model fields
 BULK_FLUXES             Surface bulk fluxes parameterization
 CCSM_FLUXES             Surface bulk fluxes parameterization
 CHARNOK                 Charnok surface roughness from wind stress
 CRAIG_BANNER            Craig and Banner wave breaking surface flux
 CORE_FORCING            Specific humidity from CORE forcing files
 DJ_GRADPS               Parabolic Splines density Jacobian (Shchepetkin, 2002)
 DOUBLE_PRECISION        Double precision arithmetic numerical kernel
 EMINUSP                 Compute Salt Flux using E-P
 FILTERED                Writing out time-filtered output
 FILTERED_RST            Reading and writing time-filtered restart files
 GLS_MIXING              Generic Length-Scale turbulence closure
 ICE_ADVECT              Advection of ice tracers
 ICE_ALB_EC92            Ebert and Curry albedo formula
 ICE_BULK_FLUXES         Ice bulk fluxes from the atmosphere
 ICE_EVP                 Elastic-viscous-plastic ice rheology
 ICE_I_O                 Put shortwave heat into ice interior
 ICE_MK                  Mellor-Kantha ice thermodynamics
 ICE_MODEL               Include sea ice model
 ICE_MOMENTUM            Compute ice momentum equations
 ICE_SMOLAR              Advect ice tracers with MPDATA scheme
 ICE_STRENGTH_QUAD       Ice strength a quadratic function of thickness
 ICE_THERMO              Include ice thermodynamics
 KANTHA_CLAYSON          Kantha and Clayson stability function formulation
 LIMIT_BSTRESS           Limit bottom stress to maintain bottom velocity direction
 LONGWAVE_OUT            Compute outgoing longwave radiation internally
 MASKING                 Land/Sea masking
 MIX_S_UV                Mixing of momentum along constant S-surfaces
 NONLINEAR               Nonlinear Model
 NONLIN_EOS              Nonlinear Equation of State for seawater
 NO_SCORRECTION_ICE      Not including SCORRECTION under the ice
 NO_WRITE_GRID           Not Writing grid arrays into NetCDF ouput files
 N2S2_HORAVG             Horizontal smoothing of buoyancy and shear
 ONE_TRACER_SOURCE       One value per tracer for all sources
 PERFECT_RESTART         Processing perfect restart variables
 POWER_LAW               Power-law shape time-averaging barotropic filter
 PROFILE                 Time profiling activated
 K_GSCHEME               Third-order upstream advection of TKE fields
 RADIATION_2D            Use tangential phase speed in radiation conditions
 RI_SPLINES              Parabolic Spline Reconstruction for Richardson Number
 !RST_SINGLE             Double precision fields in restart NetCDF file
 SALINITY                Using salinity
 SCORRECTION             Surface salinity flux correction
 SOLAR_SOURCE            Solar Radiation Source Term
 SOLVE3D                 Solving 3D Primitive Equations
 SPLINES_VDIFF           Parabolic Spline Reconstruction for Vertical Diffusion
 SPLINES_VVISC           Parabolic Spline Reconstruction for Vertical Viscosity
 SPHERICAL               Spherical grid configuration
 SSH_TIDES               Add tidal elevation to SSH climatology
 TIDES_ASTRO             Compute astronomical (long-period tide) corrections
 TS_U3HADVECTION         Third-order upstream horizontal advection of tracers
 TS_C4VADVECTION         Fourth-order centered vertical advection of tracers
 UV_ADV                  Advection of momentum
 UV_COR                  Coriolis term
 UV_U3HADVECTION         Third-order upstream horizontal advection of 3D momentum
 UV_C4VADVECTION         Fourth-order centered vertical advection of momentum
 UV_DRAG_GRID            Spatially varying quadratic drag coefficient
 UV_QDRAG                Quadratic bottom stress
 UV_TIDES                Add tidal currents to 2D momentum climatologies
 UV_VIS2                 Harmonic mixing of momentum
 VAR_RHO_2D              Variable density barotropic mode
 VISC_GRID               Horizontal viscosity coefficient scaled by grid size

Do you have any idea about this problem? I will appreciate any of your help.

[arango]

I don't know if this behaviour will be related, but I would turn off SPLINES_VDIFF and SPLINES_VVISC and see what happens. I also noticed that you are using a different version of ROMS that has sea ice modeling. This will complicate your analysis.

[cloudy4568]

Thank you, Hernan.
I turned off SPLINES_VDIFF and SPLINES_VVISC but the result is almost unchanged. I think there are something wrong with latent heat and sensible heat calculation. The CCSM_FLUXES calculates heat flux over open water with iteration which I can't understand

[arango]

Hmmm, I noticed that you activated both BULK_FLUXES and CCSM_FLUXES. I am not familiar with that version of ROMS. However, my understanding is that you should use a single formulation. You need to choose one and not both. If the sea ice model is activated, we need BOTH freshwater (E-P) and saltwater (under ice) fluxes as the upper boundary condition for salinity. You may need to ask Kate about it.

[cloudy4568]

Ok, thank you. In this version, both of BULK_FLUXES and CCSM_FLUXES need to be activated to activate CCSM_FLUXES.

Code: Select all

# ifdef BULK_FLUXES
#  ifdef CCSM_FLUXES
      USE ccsm_flux_mod,        ONLY : ccsm_flux
#  else
      USE bulk_flux_mod,        ONLY : bulk_flux
#  endif
# endif
-----------------
# ifdef BULK_FLUXES
#  if defined FOUR_DVAR && defined NL_BULK_FLUXES
                IF (Nrun.eq.1) CALL bulk_flux (ng, tile)
#  elif defined CCSM_FLUXES
                CALL ccsm_flux (ng, tile)
#  else
                CALL bulk_flux (ng, tile)
#  endif
# endif

It seems that different versions are very different.

[kate]

You could always try turning off CCSM_FLUXES to get the other flux formulation to see if it's any better. Beyond that, I'd worry about units of all the inputs.

[kate]

What are you doing for the ice albedo? I have:

# define ALBEDO_CURVE /* for water */
# undef ICE_ALB_EC92 /* for ice */
# define ALBEDO_CSIM /* for ice */
# undef ALBEDO_FILE /* for both */

Someone pointed out that ICE_ALB_EC92 can be unstable, so maybe I should delete it.

ETA: Oh, I see you are using that one.

[cloudy4568]

Thank you, Kate.

I have turned off ICE_ALB_EC92 and tried BULK_FLUXES. But the latent heat flux calculation of BULK_FLUXES is unstable and the model blew up. And I have checked the units of all inputs, they are reasonable.

I think I'd better try directly input surface total heat flux instead of calculating it with flux formulation. Do you know how to set CPP options to make ROMS directly need shflux, swflux, swrad, sustr, and svstr?

Thank you for your help!

[kate]

Simply turn off BULK_FLUXES. However, if you have ice where the driving model does not, things could get weird.

[cloudy4568]

Yes, the compiling of ice models was failed when both BULK_FLUXES and CCSM_FLUXES are turned off.

May I confirm the units of inputs? Should all the units of input the same with varinfo.dat? But I see two units under variable.

Code: Select all

'Pair'                                             ! Input
  'surface air pressure'
  'Pascal'                                       ! [ROMS wants millibar]
  'Pair, scalar, series'
  'pair_time'
  'idPair'
  'r2dvar'
  0.01d0

'Tair'                                             ! Input
  'surface air temperature'
  'Celsius'                                        ! [Kelvin]
  'Tair, scalar, series'
  'tair_time'
  'idTair'
  'r2dvar'
  1.0d0

'Qair'                                             ! Input
  'surface air relative humidity'
  'kg/kg'                                          ! [kg/kg]
  'Qair, scalar, series'
  'qair_time'
  'idQair'
  'r2dvar'
  1.00d0   

Percentage or kg/kg for Qair? Celsius or Kelvin for Tair? Pascal or millibar for Pair?

[kate]

I am using surface forcing from JRA55-do, which has units of Pascals for pressure, kg/kg for specific humidity, and Kelvin for air temperature. That's how my setup works.

[cloudy4568]

Thank you Kate, I mean from where we can know what unit of each input that Roms wants?

[kate]

CCSM_FLUXES want kg/kg for specific humidity while BULK_FLUXES has some weird check to see if you are providing relative humidity vs. specific humidity.

ROMS wants Celsius for air temp. I hacked the reading of it to accept Kelvin if that's what you have.

ROMS wants millibars, hence the 0.01 factor in varinfo.dat.

[cloudy4568]

Problem solved . I changed unit of specific humidity from percentage to kg/kg, and unit of air pressure from millibar to Pascals, then the result of water temperature is reasonable. Seems that different versions of ROMS wants different units for the same variable.
Thank you very much, Kate, I really appreciate your help!