SST big problem use ROMS sea ice define BULK_FLUX

[zhaobiao]

could anyone help me? i am in trouble when i make a climatology run in Arctic with ROMS sea ice.Fist, I make forcing file which contain many variabls include Pair Tair Qair Uwind Vwind lwrf_down swrad rain.The data i used is monthly average COADS data. my netcdf file attributes:
dimensions:
xi_rho = 1159 ;
eta_rho = 999 ;
xi_psi = 1158 ;
eta_psi = 998 ;
xi_u = 1158 ;
eta_u = 999 ;
xi_v = 1159 ;
eta_v = 998 ;
bulk_time = 12 ;
pair_time = 12 ;
tair_time = 12 ;
qair_time = 12 ;
srf_time = 12 ;
albedo_time = 12 ;
lrf_time = 12 ;
rain_time = 12 ;
wind_time = 12 ;
cloud_time = 12 ;
evap_time = 12 ;
sst_time = 12 ;
variables:
double lon_rho(eta_rho, xi_rho) ;
lon_rho:long_name = "lon at rho points" ;
lon_rho:coordinates = "lon_rho lat_rho" ;
double lat_rho(eta_rho, xi_rho) ;
lat_rho:long_name = "lat at rho points" ;
lat_rho:coordinates = "lon_rho lat_rho" ;
double bulk_time(bulk_time) ;
bulk_time:long_name = "bulk formulation execution time" ;
bulk_time:units = "days" ;
bulk_time:cycle_length = 360. ;
double pair_time(pair_time) ;
pair_time:long_name = "surface air pressure time" ;
pair_time:units = "days" ;
pair_time:cycle_length = 360. ;
double tair_time(tair_time) ;
tair_time:long_name = "surface air temerature time" ;
tair_time:units = "days" ;
tair_time:cycle_length = 360. ;
double sst_time(sst_time) ;
sst_time:long_name = "sea surface temerature time" ;
sst_time:units = "days" ;
sst_time:cycle_length = 360. ;
double qair_time(qair_time) ;
qair_time:long_name = "surface air humidity time" ;
qair_time:units = "days" ;
qair_time:cycle_length = 360. ;
double srf_time(srf_time) ;
srf_time:long_name = "solar shortwave radiation time" ;
srf_time:units = "days" ;
srf_time:cycle_length = 360. ;
double lrf_time(lrf_time) ;
lrf_time:long_name = "longwave radiation time" ;
lrf_time:units = "days" ;
lrf_time:cycle_length = 360. ;
double albedo_time(albedo_time) ;
albedo_time:long_name = "albedo time" ;
albedo_time:units = "days" ;
albedo_time:cycle_length = 360. ;
double wind_time(wind_time) ;
wind_time:long_name = "wind time" ;
wind_time:units = "days" ;
wind_time:cycle_length = 360. ;
double rain_time(rain_time) ;
rain_time:long_name = "rain fall rate time" ;
rain_time:units = "days" ;
rain_time:cycle_length = 360. ;
double evap_time(evap_time) ;
evap_time:long_name = "evaporation rate time" ;
evap_time:units = "days" ;
evap_time:cycle_length = 360. ;
double cloud_time(cloud_time) ;
cloud_time:long_name = "cloud time" ;
cloud_time:units = "days" ;
cloud_time:cycle_length = 360. ;
double Tair(tair_time, eta_rho, xi_rho) ;
Tair:long_name = "surface air temperature" ;
Tair:units = "Celsius" ;
Tair:coordinates = "lon_rho lat_rho" ;
Tair:time = "tair_time" ;
double Qair(qair_time, eta_rho, xi_rho) ;
Qair:long_name = "surface relative humidity" ;
Qair:units = "kg/kg" ;
Qair:coordinates = "lon_rho lat_rho" ;
Qair:time = "qair_time" ;
double rain(rain_time, eta_rho, xi_rho) ;
rain:long_name = "rain fall rate" ;
rain:units = "kilogram meter-2 second-1" ;
rain:coordinates = "lon_rho lat_rho" ;
rain:time = "rain_time" ;
double lwrad_down(lrf_time, eta_rho, xi_rho) ;
lwrad_down:long_name = "downwelling longwave radiation flux" ;
lwrad_down:units = "Watts meter-2" ;
lwrad_down:coordinates = "lon_rho lat_rho" ;
lwrad_down:time = "lrf_time" ;
lwrad_down:positive = "downward flux, warming water" ;
double swrad(srf_time, eta_rho, xi_rho) ;
swrad:long_name = "solar shortwave radiation" ;
swrad:units = "Watts meter-2" ;
swrad:coordinates = "lon_rho lat_rho" ;
swrad:time = "srf_time" ;
swrad:positive = "downward flux, heating water" ;
double SST(sst_time, eta_rho, xi_rho) ;
SST:long_name = "sea surface temperature" ;
SST:units = "Celsius" ;
SST:coordinates = "lon_rho lat_rho" ;
SST:time = "sst_time" ;
double Pair(pair_time, eta_rho, xi_rho) ;
Pair:long_name = "surface air pressure" ;
Pair:units = "Pascal" ;
Pair:coordinates = "lon_rho lat_rho" ;
Pair:time = "pair_time" ;
double Uwind(wind_time, eta_rho, xi_rho) ;
Uwind:long_name = "u-wind" ;
Uwind:units = "m/s" ;
Uwind:coordinates = "lon_rho lat_rho" ;
Uwind:time = "wind_time" ;
double Vwind(wind_time, eta_rho, xi_rho) ;
Vwind:long_name = "v-wind" ;
Vwind:units = "m/s" ;
Vwind:coordinates = "lon_rho lat_rho" ;
Vwind:time = "wind_time" ;
my CPP as follows:
#define CURVGRID
#define MASKING
#define NONLIN_EOS
#define SOLVE3D
#define SALINITY
#ifdef SOLVE3D
# define SPLINES
#endif
#undef FLOATS
#undef STATIONS
#undef WET_DRY

#undef T_PASSIVE
#ifdef T_PASSIVE
# define ANA_PASSIVE
# define TRC_PSOURCE
# define ANA_TRC_PSOURCE
# define AGE_PASSIVE
#endif

/* ice */

#ifdef SOLVE3D
# define ICE_MODEL
# ifdef ICE_MODEL
# define ANA_ICE
# define OUTFLOW_MASK
# undef FASTICE_CLIMATOLOGY
# define ICE_THERMO
# define ICE_MK
# undef ICE_SMOOTH
# define ICE_MOMENTUM
# define ICE_MOM_BULK
# define ICE_EVP
# define ICE_ADVECT
# define ICE_SMOLAR
# define ICE_UPWIND
# define ICE_BULK_FLUXES
# define ANA_AIOBC
# define ANA_HIOBC
define ANA_HSNOBC
# endif
#endif

/* output stuff */

#define NO_WRITE_GRID
#define INLINE_2DIO
#define OUT_DOUBLE
#define AVERAGES
#define AVERAGES2
#ifdef SOLVE3D
# undef AVERAGES_DETIDE
# undef DIAGNOSTICS_TS
#endif
#undef DIAGNOSTICS_UV

/* advection, dissipation, pressure grad, etc. */

#ifdef SOLVE3D
# define DJ_GRADPS
#endif

#define UV_ADV
#define UV_COR
#undef UV_SADVECTION

#ifdef SOLVE3D
# define TS_U3HADVECTION
# define TS_C4VADVECTION
# undef TS_MPDATA
#endif

#define UV_VIS2
#undef UV_SMAGORINSKY
#undef VISC_3DCOEF
#define MIX_S_UV
#define VISC_GRID
#undef SPONGE

#ifdef SOLVE3D
# define TS_DIF2
# define MIX_GEO_TS
# define DIFF_GRID
#endif

/* vertical mixing */

#ifdef SOLVE3D
# undef WTYPE_GRID

# undef LMD_MIXING
# ifdef LMD_MIXING
# define LMD_RIMIX
# define LMD_CONVEC
# define LMD_SKPP
# undef LMD_BKPP
# define LMD_NONLOCAL
# define LMD_SHAPIRO
# undef LMD_DDMIX
# endif

# undef GLS_MIXING
# define MY25_MIXING
# if defined GLS_MIXING || defined MY25_MIXING
# define KANTHA_CLAYSON
# define N2S2_HORAVG
# endif
#endif

/* surface forcing */

#ifdef SOLVE3D
# define CORE_FORCING
# define BULK_FLUXES
# define CCSM_FLUXES
# if defined BULK_FLUXES || defined CCSM_FLUXES
# define LONGWAVE_OUT
# define ANA_CLOUD
# define DIURNAL_SRFLUX
# define SOLAR_SOURCE
# define EMINUSP
# undef ANA_SRFLUX
# undef ALBEDO
# define ALBEDO_CURVE /* for water */
# define ICE_ALB_EC92 /* for ice */
# undef ALBEDO_CSIM /* for ice */
# undef ALBEDO_FILE /* for both */
# undef LONGWAVE
# endif
#endif

/* surface and side corrections */

#ifdef SOLVE3D
# define SCORRECTION
# undef QCORRECTION
#endif
#define RADIATION_2D

/* roms quirks */

#ifdef SOLVE3D
# define ANA_BSFLUX
# define ANA_BTFLUX
#else
# define ANA_SMFLUX
#endif


i set ice free at the beginning of model start. after integrated 12 days, SST can upto 96 degree centigrade. I check the latent,sensible and net heat flux which cumputed by bulk formula.It's not normal. the values of net heat flux are over 1000 at some regions which correspond to unnormal SST. I don't know why. it's look like that latent and sensible heat flux are vary sensitive to wind speeds. what cause simulated SST abnormal ? the forcing data i used ? or bulk formula? so I upload some figures for you to diagnose. I would appreciate it if sombody can figure it out. thanks in advance!

[Paul_Budgell]

Looks like some issue with units. Have you converted sea level pressure from Pascals to Hectopascals (mb) in varinfo.dat?

[zhaobiao]

Looks like some issue with units. Have you converted sea level pressure from Pascals to Hectopascals (mb) in varinfo.dat?

Dear Paul.Thank you for your reply. I got this code from Kate. the description of Pair in my varinfo.dat is like this:

'Pair' ! Input
'surface air pressure'
'Pascal' ! [ROMS wants millibar]
'Pair, scalar, series'
'pair_time'
'idPair'
'r2dvar'
0.01d0
did you mean that i should convert it into this?

'Pair' ! Input
'surface air pressure'
'millibar' ! [ROMS wants millibar]
'Pair, scalar, series'
'pair_time'
'idPair'
'r2dvar'
0.01d0
or should i convert sea level pressure from Pascals to Hectopascals (mb) in my forcing netcdf file?

[Paul_Budgell]

No, that should be OK - the important thing is the 0.01 factor.
You need to see what the forcing variables look like as they are being read in (their ranges are written out at input times). That will tell you if your variables are in the right range, i.e., you can catch problems like air temps being K instead of C, or specific humidities being inconsistent with the cppdef options, for example.

BTW, you need a # in front of define ANA_HSNOBC, but I don't think that has anything to do with your problem.

[zhaobiao]

No, that should be OK - the important thing is the 0.01 factor.
You need to see what the forcing variables look like as they are being read in (their ranges are written out at input times). That will tell you if your variables are in the right range, i.e., you can catch problems like air temps being K instead of C, or specific humidities being inconsistent with the cppdef options, for example.

BTW, you need a # in front of define ANA_HSNOBC, but I don't think that has anything to do with your problem.

I confgure a new experiment. I multiply the Pair value in my focing netcdf file by 0.01. the Pair valube will be in the range of 9.72~10.4. The following are ouputs from the model:
GET_2DFLD - surface u-wind component, t = 345 00:00:00
(Rec=0000012, Index=1, File: roms_blk.nc)
(Tmin= 15.0000 Tmax= 345.0000)
(Min = -8.64797974E+00 Max = 1.76889019E+01)
GET_2DFLD - surface v-wind component, t = 345 00:00:00
(Rec=0000012, Index=1, File: roms_blk.nc)
(Tmin= 15.0000 Tmax= 345.0000)
(Min = -2.44668770E+01 Max = 1.03587303E+01)
GET_2DFLD - surface air pressure, t = 345 00:00:00
(Rec=0000012, Index=1, File: roms_blk.nc)
(Tmin= 15.0000 Tmax= 345.0000)
(Min = 9.72266123E+00 Max = 1.04074050E+01)
GET_2DFLD - surface air temperature, t = 345 00:00:00
(Rec=0000012, Index=1, File: roms_blk.nc)
(Tmin= 15.0000 Tmax= 345.0000)
(Min = -4.40002213E+01 Max = 8.36931229E+00)
GET_2DFLD - surface air relative humidity, t = 345 00:00:00
(Rec=0000012, Index=1, File: roms_blk.nc)
(Tmin= 15.0000 Tmax= 345.0000)
(Min = 5.44482613E-01 Max = 1.00000000E+00)
GET_2DFLD - solar shortwave radiation flux, t = 345 00:00:00
(Rec=0000012, Index=1, File: roms_blk.nc)
(Tmin= 15.0000 Tmax= 345.0000)
(Min = 0.00000000E+00 Max = 6.15148763E-06)
GET_2DFLD - downwelling longwave radiation flux, t = 345 00:00:00
(Rec=0000012, Index=1, File: roms_blk.nc)
(Tmin= 15.0000 Tmax= 345.0000)
(Min = 3.02729104E-05 Max = 8.81157885E-05)
GET_2DFLD - rain fall rate, t = 345 00:00:00
(Rec=0000012, Index=1, File: roms_blk.nc)
(Tmin= 15.0000 Tmax= 345.0000)
(Min = -8.48331606E-05 Max = 2.65164552E-04)
GET_2DFLD - sea surface salinity climatology, t = 345 00:00:00
(Rec=0000012, Index=1, File: roms_frc.nc)
(Tmin= 15.0000 Tmax= 345.0000)
(Min = 3.59050012E+00 Max = 3.55107002E+01)

What is really weird is that it seems the simulated SST result is normal now. please see the picture below. But the Pair i made is wrong. is'nt it? that make me confused!
why?

[Paul_Budgell]

I think I see your problem now. Just put your Pair the way you originally had it and #undef CORE_FORCING. You will also need to make sure you have the right setting for relative humidity in your varinfo.dat.

You clearly are using relative, not specific humidity, so the cppdef above will cause ROMS to use the wrong formula for vapour pressure in the computation of latent heat flux. Your reducing Pair by 0.01 partially compensated for the wrong humidity in the vapour pressure calculation, which is why your results from that run look more reasonable.

What atmospheric forcing data set are you using? It doesn't look like it would be CORE, if I am interpreting the humidity values correctly.

[zhaobiao]

I think I see your problem now. Just put your Pair the way you originally had it and #undef CORE_FORCING. You will also need to make sure you have the right setting for relative humidity in your varinfo.dat.

You clearly are using relative, not specific humidity, so the cppdef above will cause ROMS to use the wrong formula for vapour pressure in the computation of latent heat flux. Your reducing Pair by 0.01 partially compensated for the wrong humidity in the vapour pressure calculation, which is why your results from that run look more reasonable.

What atmospheric forcing data set are you using? It doesn't look like it would be CORE, if I am interpreting the humidity values correctly.

you are right!I used ralative humidity. I have checked the bulk_flux.F and ccsm_bulk.F. i found that ccsm_bulk.F can not convert relative humidity into specific humidity. so i have three ways to solve this problem.
Fist, As you said #undef CORE_FORCING #define BULK_FLUX and #undef CCSM_BULK.
Second, add transformation equation to ccsm_bulk.F. #define CCSM_BULK #def CORE_FORCING #define BULK_FLUX and #define CCSM_FLUX.
The last, I can use specific humidity instead of relative humidity in my forcing data.
now I am testing the fist way. The results are normal. Thank you Paul! thanks for you patience!