Stand-alone CICE6 Enhancement -- building greater flexibility for forcing/initial conditions

[dpath2o]

Hello,

This is somewhat a follow-on to the work that I'm doing with CICE6 and my research interest into building a prognostic Antarctic fast ice model. As I've mentioned before I'm in the beginning stages of this work and working part-time on the topic. As such over the past couple of months I've been making some in-roads into CICE6 ingesting (reading-in) two global atmosphere and ocean climatological datasets, ERA5 and BRAN, respectively, and have done a little bit towards documenting this effort. This thread that I'm opening here as I could not find an `issue' or `project' publicly available on CICE GitHub that is aligned with this topic.

Essentially, I would like to build greater flexibility into CICE6 for processing other climatological datasets in stand-alone mode. I know that CICE benefits greatly from running in a coupled mode, but I still think, with the quality climatological datasets that are available, that there is still a scientific application for running CICE in this mode. Especially for testing CICE and adding features, ya?

Cheers,
Dan

 

[dbailey]

Hi Dan,

In terms of reading in other datasets, CICE standalone still requires that the data be interpolated to the CICE grid you are using. We have talked about some type of on the fly regridding of data, but in the end this is more for the coupled models. For example, E3SM and CESM have a data atmosphere component that use mapping files to regrid atmospheric datasets on the fly. This might be an avenue to consider.

Dave

 

[eclare@lanl_gov]

Hi Dan,

We do recommend that CICE be used in coupled mode, or at least that any developments or results generated in standalone mode be tested in a coupled configuration. The forcing data and routines that are currently available in CICE and Icepack are intended primarily for software testing. However, they do provide a "model" for introducing forcing from other sources, which you are welcome to do. And from there, of course you can introduce other forcing variables to see what difference they might make in the stand-alone configuration. It is a legitimate use of the model as a tool. However, once again I will emphasize that care must be taken when interpreting standalone CICE results.

In your initial question above, it's not completely clear what you are referring to with "build greater flexibility into CICE6 for processing other climatological datasets". There is a balance between the complexity of the model and all its infrastructure, and its usefulness. As Dave alludes to above, generalizing the CICE infrastructure to use other datasets more easily by adding a regridding functionality is an option, but one that we would need to weigh in light of increased code complexity. On the other hand, adding another dataset with expanded functionality in terms of variables and accessible physics, but which is already interpolated to our test grids, would be a lighter lift for us. Groups using other meshes have to regrid all of their data to those meshes anyhow, so we don't necessarily need to take on that responsibility. Is this what you have in mind?

Best,
e

 

[dpath2o]

Hi Elizabeth and Dave,

Thank you very much for your comments and thoughts on the matter.

My question was exploratory -- to open a discussion and see the value for effort into to testing sub-routine(s) functionality to the stand-alone forcing of other gridded climatological datasets, namely ERA5 and BRAN. However, I really need to concentrate more on the scientific outputs rather than the technical parts of work. So I think for this first study I will make my re-gridded datasets conform to what the stand-alone HYCOM input forcing routines are wanting (i.e. what CICE is expecting the stand-alone HYCOM to look like, possibly there is another forcing sub-routine that will meet these needs, like the NCAR forcing sub-routine).

Cheers,
Dan

 

[dpath2o]

Hi again,

Hopefully a quick question.

When I compare the data variables in

../CICE_data/forcing/gx1/JRA55/8XDAILY/JRA55_03hr_forcing_2005.nc

with those listed in CICE documentation 4.5 under "Atmosphere" there is a significant difference.

Explicitly, here is the difference:

JRA55_03hr_forcing_YYYY.nc	CICE documentation atmosphere list
airtmp	Tair = air temperature
dlwsfc	flw = incoming longwave radiation
glbrad	sw = incoming shortwave radiation
spchmd	Qa = specific humidity
ttlpcp
wndewd	uatm = model grid i-direction wind velocity component
wndnwd	vatm = model grid j-direction wind velocity component​

From the CICE documentation the remainder of atmospheric variables that should be provided are:

• zlvl = atmosphere level height (m)
• strax = model grid i-direction wind stress (N/m^2)
• stray = model grid j-direction wind stress (N/m^2)
• potT = air potential temperature (K)
• rhoa = air density (kg/m^3)
• swvdr = sw down, visible, direct (W/m^2)
• swvdf = sw down, visible, diffuse (W/m^2)
• swidr = sw down, near IR, direct (W/m^2)
• swidf = sw down, near IR, diffuse (W/m^2)
• frain = rainfall rate (kg/m^2 s)
• fsnow = snowfall rate (kg/m^2 s)

The last two ('frain' and 'fsnow') I suppose are derived/calculated internally in CICE if the forcing field provides a precipitation rate (i.e. 'ttlpcp' in the JRA55 forcing file)? This gets at my main questions here:

• are the remaining fields (those in the dot points above) are derived/calculated internally in CICE from the seven input variables that are in the JRA55 forcing file?
• how would one provide CICE6 with all the atmospheric variables that are in that dot-point list above?
• "zlvl", is this a static parameter (i.e. 10 meters or maybe 2 meters)?

I'm most interested in providing CICE with accurate solar radiation both incident and diffuse as well as parsed by wavelength. However, it might be interesting to see if providing CICE with parsed rainfall rate versus snowfall rate is of benefit to growth rate accuracy. I imagine the other input variables, names "potT", "rhoa", "strax", and "stray" are derived/calculated by CICE6 as well, but I am interested on providing these to CICE instead of CICE computing them.

Lastly, I think there is a minor discrepancy with the documentation. The above dot-point list is provided in CICE documentation 4.5 and all the variable names match those provided in the Index of primary variables and parameters , with the exception of "swidr" and "swidf". In the index those variables indicated as "swv(n)dr(f)". So by this convention the CICE documentation 4.5 should have "swndr" and "swndf". However, and not to confuse this very minor point, but in the code, namely file ice_flux.F90 these variables are named "swidr" and "swidf". Obviously a very minor point, but it was a bit annoying when trying to search for where these variables are getting created in the code to try to track down answer to my main questions above. Unfortunately, I was not successful in coming to a conclusion and hence why I've asked the above questions.

As always, thank you for your time.

Cheers,
Dan

 

[david_hebert@nrlssc_navy_mil]

Hi again,

Hopefully a quick question.

When I compare the data variables in

with those listed in CICE documentation 4.5 under "Atmosphere" there is a significant difference.

Explicitly, here is the difference:

JRA55_03hr_forcing_YYYY.nc	CICE documentation atmosphere list
airtmp	Tair = air temperature
dlwsfc	flw = incoming longwave radiation
glbrad	sw = incoming shortwave radiation
spchmd	Qa = specific humidity
ttlpcp
wndewd	uatm = model grid i-direction wind velocity component
wndnwd	vatm = model grid j-direction wind velocity component​

From the CICE documentation the remainder of atmospheric variables that should be provided are:

• zlvl = atmosphere level height (m)
• strax = model grid i-direction wind stress (N/m^2)
• stray = model grid j-direction wind stress (N/m^2)
• potT = air potential temperature (K)
• rhoa = air density (kg/m^3)
• swvdr = sw down, visible, direct (W/m^2)
• swvdf = sw down, visible, diffuse (W/m^2)
• swidr = sw down, near IR, direct (W/m^2)
• swidf = sw down, near IR, diffuse (W/m^2)
• frain = rainfall rate (kg/m^2 s)
• fsnow = snowfall rate (kg/m^2 s)

The last two ('frain' and 'fsnow') I suppose are derived/calculated internally in CICE if the forcing field provides a precipitation rate (i.e. 'ttlpcp' in the JRA55 forcing file)? This gets at my main questions here:

• are the remaining fields (those in the dot points above) are derived/calculated internally in CICE from the seven input variables that are in the JRA55 forcing file?
• how would one provide CICE6 with all the atmospheric variables that are in that dot-point list above?
• "zlvl", is this a static parameter (i.e. 10 meters or maybe 2 meters)?

I'm most interested in providing CICE with accurate solar radiation both incident and diffuse as well as parsed by wavelength. However, it might be interesting to see if providing CICE with parsed rainfall rate versus snowfall rate is of benefit to growth rate accuracy. I imagine the other input variables, names "potT", "rhoa", "strax", and "stray" are derived/calculated by CICE6 as well, but I am interested on providing these to CICE instead of CICE computing them.

Lastly, I think there is a minor discrepancy with the documentation. The above dot-point list is provided in CICE documentation 4.5 and all the variable names match those provided in the Index of primary variables and parameters , with the exception of "swidr" and "swidf". In the index those variables indicated as "swv(n)dr(f)". So by this convention the CICE documentation 4.5 should have "swndr" and "swndf". However, and not to confuse this very minor point, but in the code, namely file ice_flux.F90 these variables are named "swidr" and "swidf". Obviously a very minor point, but it was a bit annoying when trying to search for where these variables are getting created in the code to try to track down answer to my main questions above. Unfortunately, I was not successful in coming to a conclusion and hence why I've asked the above questions.

As always, thank you for your time.

Cheers,
Dan

Hi Dan,

First, thanks for pointing out the documentation error. We will correct the variable names in the docs to read 'swidr, swidf' shortly.

For the forcing, using JRA55 in an uncoupled CICE, the data read in are further processed in subroutine ice_forcing.F90:prepare_forcing. In prepare_forcing zlvl0 is set, as well as checking air temp to convert precip to rain or snow. It also breaks up the shortwave into visibile and IR, does wind rotations, etc here. In a coupled setting, most of those would be passed into CICE via a coupling interface. Hopefully that will help with some of the forcing questions you have. Let us know if we can help more. Thanks!

David

 

[eclare@lanl_gov]

Dan, let us know if this answers your question or if something about the different forcing configurations is still confusing. In a nutshell, the answer is "users can change or write code to do whatever they want"!

 

[dpath2o]

Hi David and Elizabeth,

Thank you for your responses. I appreciate them very much.

David, I've had a cursory look into the subroutine prepare_forcing.

Elizabeth, I think your comment in quotes, is pretty much where I need to go with my effort.

Where I think I might get a good deal benefit from is if I create a diagram that maps the forcing_nml options to ice_forcing.F90 subroutines and the various input parameters being read-in (forced) at each time step, additionally with those parameters/variables that are static/hard-wired/defaulted if not present in the forcing file (gridded dataset). There is of course orders of magnitude differences that I'm trying to keep in mind here as well -- i.e. solar radiation significance vs deep ocean heat flux.

In fact, with regards to the later (deep ocean heat flux), I have created a function (sub-routine) in my ocean re-gridding script that you may care to comment on or incorporate into CICE ice_forcing.F90 sub-routine. If you find a mistake I'd be very grateful to know.

def compute_ocn_heat_flux( T, S, mld, temp_var_name='temp' , salt_var_name='salt', mld_var_name='mld', lat_name='xt_ocean' ):
    '''
    Input 'T', 'S' and 'mld' are assumed to be xarray datasets.

    Optional inputs are provided to allow user to use different xarray datasets as opposed
    to the dataset that the function is originally written for -- BRAN. So provided the variable
    names to each of the three input datasets, as well as the name of latitude coordinate.

    Latitude/longitude (horizontal spatial coordinates) are assumed to be in decimal degrees.

    Given ocean temperature and salinity, compute the deep ocean heat flux. The
    approach is first to accurately compute heat capacity and then use the relationship
    Q = cp + dT ; where T is the temperature across the layer and cp is the
    heat capacity calculated at a depth. The heat capacity comes from the reference:
    Fofonff, P. and Millard, R.C. Jr
    Unesco 1983. Algorithms for computation of fundamental properties of seawater,
    1983. Unesco Tech. Pap. in Mar. Sci., No. 44, 53 pp.
    '''
    import numpy as np
    sst    = T.sel(st_ocean=0      ,method='nearest')[temp_var_name]
    T_mld  = T.sel(st_ocean=mld[mld_var_name],method='nearest')[temp_var_name]
    S_mld  = S.sel(st_ocean=mld[mld_var_name],method='nearest')[salt_var_name]
    # depth to pressure, convert lats to radians
    latrad = np.sin(np.abs(T[lat_name])*(np.pi/180))
    C1     = 5.92e-3 + latrad**(2*5.25e-3)
    P_mld  = ( (1-C1) - np.sqrt( ( (1-C1)**2 ) - (8.84e-6*mld[mld_var_name]) ) ) / 4.42e-6
    # to convert db to Bar as used in Unesco routines
    P_mld  = P_mld/10
    # HEAT CAPACITY
    c0     =  4.2174e3
    c1     = -3.720283
    c2     =   .1412855
    c3     = -2.654387e-3
    c4     = 2.093236e-5
    a0     = -7.64357
    a1     =   .1072763
    a2     = -1.38385e-3
    b0     =   .1770383
    b1     = -4.07718e-3
    b2     =  5.148e-5
    Cpst0  = c0 + c1*T_mld + c2*T_mld**2 + c3*T_mld**3 + c4*T_mld**4 + (a0 + a1*T_mld + a2*T_mld**2)*S_mld + (b0 + b1*T_mld + b2*T_mld**2)*S_mld*np.sqrt(S_mld)
    a0     = -4.9592e-1
    a1     =  1.45747e-2
    a2     = -3.13885e-4
    a3     =  2.0357e-6
    a4     =  1.7168e-8
    b0     =  2.4931e-4
    b1     = -1.08645e-5
    b2     =  2.87533e-7
    b3     = -4.0027e-9
    b4     =  2.2956e-11
    c0     = -5.422e-8
    c1     =  2.6380e-9
    c2     = -6.5637e-11
    c3     =  6.136e-13
    dCp0t0 = (a0 + a1*T_mld + a2*T_mld**2 + a3*T_mld**3 + a4*T_mld**4)*P_mld + (b0 + b1*T_mld + b2*T_mld**2 + b3*T_mld**3 + b4*T_mld**4)*P_mld**2 + (c0 + c1*T_mld + c2*T_mld**2 + c3*T_mld**3)*P_mld**3
    d0     =  4.9247e-3
    d1     = -1.28315e-4
    d2     =  9.802e-7
    d3     =  2.5941e-8
    d4     = -2.9179e-10
    e0     = -1.2331e-4
    e1     = -1.517e-6
    e2     =  3.122e-8
    f0     = -2.9558e-6
    f1     =  1.17054e-7
    f2     = -2.3905e-9
    f3     =  1.8448e-11
    g0     =  9.971e-8
    h0     =  5.540e-10
    h1     = -1.7682e-11
    h2     =  3.513e-13
    j1     = -1.4300e-12
    S3_2   = S_mld*np.sqrt(S_mld)
    dCpstp = ((d0 + d1*T_mld + d2*T_mld**2 + d3*T_mld**3 + d4*T_mld**4)*S_mld + (e0 + e1*T_mld + e2*T_mld**2)*S3_2)*P_mld + ((f0 + f1*T_mld + f2*T_mld**2 + f3*T_mld**3)*S_mld + g0*S3_2)*P_mld**2 + ((h0 + h1*T_mld + h2*T_mld**2)*S_mld + j1*T_mld*S3_2)*P_mld**3
    cp     = Cpst0 + dCp0t0 + dCpstp
    # HEAT FLUX
    return cp + (sst-T_mld)

 

[eclare@lanl_gov]

@dbailey Is this roughly what's done when computing deep ocean heat flux from CESM output?

 

[dbailey]

The best reference for the Q-flux calculation is here:

FAQ: Data ocean slab mode (DOCN-SOM)

I get numerous questions about the new SOM in the data ocean model of CESM. So, I thought it appropriate to add a sort of FAQ here. 1. Why does the SOM mode no longer run out of the box? Earlier versions were set up with a test SOM forcing file (pop_frc.1x1d.090130.nc) This forcing file was...

bb.cgd.ucar.edu

This points to a white paper on how we compute this along with some scripts. Also, there is a pointer to the Bitz et al. manuscript on the method.

 

[dpath2o]

Thank you David. I read through that white paper and there are some differences in how I was calculating Q_flx to what is shown as a more scientifically robust method in the white paper and associated Bitz et. al. paper. So I truly appreciate that.

If of any interest to you or others, I'll include my work here with results.

def compute_ocn_heat_flux_at_mixed_layer_depth(T, S, M, F_net, dTdt,
                                               T_var_name='temp',
                                               S_var_name='salt',
                                               M_var_name='mld',
                                               lat_name='xt_ocean'):
    '''
    '''
    T_mld  = T.sel(st_ocean=M[M_var_name],method='nearest')[T_var_name]
    S_mld  = S.sel(st_ocean=M[M_var_name],method='nearest')[S_var_name]
    cp     = compute_ocn_heat_capacity_at_depth(T_mld, S_mld, M,
                                                D_var_name=M_var_name,
                                                lat_name=lat_name)
    rho    = compute_ocn_density_at_depth(T_mld, S_mld, M,
                                          D_var_name=M_var_name,
                                          lat_name=lat_name)
    return ( F_net - rho*cp*M[M_var_name] * dTdt[T_var_name] )

############################################################################
def compute_ocn_density_at_depth(T, S, D, D_var_name='mld', lat_name='xt_ocean'):
    '''
    UNESCO (1981) Tenth report of the joint panel on oceanographic
    tables and stan- dards. UNESCO Technical Papers in Marine Science, Paris, 25 p
    '''
    # depth to pressure (bar)
    P = (D.mld * 1027 * 9.8)/1000
    # standar ocean water
    a0       =   999.842594
    a1       =     6.793953e-2
    a2       =    -9.095290e-3
    a3       =     1.001685e-4
    a4       =    -1.120083e-6
    a5       =     6.536332e-9
    rho_SMOW = a0 + a1*T + a2*T**2 + a3*T**3 + a4*T**4 +a5*T**5
    b0       =     8.2449e-1
    b1       =    -4.0899e-3
    b2       =     7.6438e-5
    b3       =    -8.2467e-7
    b4       =     5.3875e-9
    c0       =    -5.7246e-3
    c1       =     1.0227e-4
    c2       =    -1.6546e-6
    d0       =     4.8314e-4
    B1       = b0 + b1*T + b2*T**2 + b3*T**3 + b4*T**4
    C1       = c0 + c1*T + c2*T**2
    rho_p0   = rho_SMOW + B1*S + C1*S**1.5 + d0*S**2
    e0       = 19652.21
    e1       =   148.4206
    e2       =    -2.327105
    e3       =     1.360477e-2
    e4       =    -5.155288e-5
    f0       =    54.674600
    f1       =    -0.603459
    f2       =     1.099870e-2
    f3       =    -6.167e-5
    g0       =     7.9440e-2
    g1       =     1.6483e-2
    g2       =    -5.3009e-4
    G1       = g0 + g1*T + g2*T**2
    F1       = f0 + f1*T + f2*T**2 + f3*T**3
    Kw       = e0 + e1*T + e2*T**2 + e3*T**3 + e4*T**4
    K0       = Kw + F1*S + G1*S**1.5
    h0       =     3.23990
    h1       =     1.43713e-3
    h2       =     1.16092e-4
    h3       =    -5.77905e-7
    i0       =     2.28380e-3
    i1       =    -1.09810e-5
    i2       =    -1.60780e-6
    j0       =     1.91075e-4
    k0       =     8.50935e-5
    k1       =    -6.12293e-6
    k2       =     5.27870e-8
    m0       =    -9.9348e-7
    m1       =     2.0816e-8
    m2       =     9.1697e-10
    Bw       = k0 + k1*T + k2*T**2
    B2       = Bw + (m0 + m1*T + m2*T**2)*S
    Aw       = h0 + h1*T + h2*T**2 + h3*T**3
    A1       = Aw + (i0 + i1*T + i2*T**2)*S + j0*S**1.5
    K        = K0 + A1*P + B2*P**2
    return ( rho_p0 / (1 - ( P / K )) )

############################################################################
def compute_ocn_heat_capacity_at_depth(T, S, D, D_var_name='mld', lat_name='xt_ocean'):
    '''
    Given ocean temperature, salinity, and mixed layer depth, compute the deep
    ocean heat capacity at the mixed layer depth. The heat capacity comes from
    the reference:
    Fofonff, P. and Millard, R.C. Jr
    Unesco 1983. Algorithms for computation of fundamental properties of seawater,
    1983. Unesco Tech. Pap. in Mar. Sci., No. 44, 53 pp.
    '''
    # depth to pressure (bar)
    P = (D.mld * 1027 * 9.8)/1000
    # HEAT CAPACITY
    a0   = -7.64357
    a1   =   .1072763
    a2   = -1.38385e-3
    b0   =   .1770383
    b1   = -4.07718e-3
    b2   =  5.148e-5
    c0   =  4.2174e3
    c1   = -3.720283
    c2   =   .1412855
    c3   = -2.654387e-3
    c4   = 2.093236e-5
    A1   = c0 + c1*T + c2*T**2 + c3*T**3 + c4*T**4
    A2   = (a0 + a1*T + a2*T**2)*S
    A3   = (b0 + b1*T + b2*T**2)*S*np.sqrt(S)
    Cp0  = A1 + A2 + A3
    d0   = -4.9592e-1
    d1   =  1.45747e-2
    d2   = -3.13885e-4
    d3   =  2.0357e-6
    d4   =  1.7168e-8
    e0   =  2.4931e-4
    e1   = -1.08645e-5
    e2   =  2.87533e-7
    e3   = -4.0027e-9
    e4   =  2.2956e-11
    f0   = -5.422e-8
    f1   =  2.6380e-9
    f2   = -6.5637e-11
    f3   =  6.136e-13
    B1   = (d0 + d1*T + d2*T**2 + d3*T**3 + d4*T**4)*P
    B2   = (e0 + e1*T + e2*T**2 + e3*T**3 + e4*T**4)*P**2
    B3   = (f0 + f1*T + f2*T**2 + f3*T**3)*P**3
    dCp0 =  B1 + B2 + B3 
    d0   =  4.9247e-3
    d1   = -1.28315e-4
    d2   =  9.802e-7
    d3   =  2.5941e-8
    d4   = -2.9179e-10
    e0   = -1.2331e-4
    e1   = -1.517e-6
    e2   =  3.122e-8
    f0   = -2.9558e-6
    f1   =  1.17054e-7
    f2   = -2.3905e-9
    f3   =  1.8448e-11
    g0   =  9.971e-8
    h0   =  5.540e-10
    h1   = -1.7682e-11
    h2   =  3.513e-13
    j1   = -1.4300e-12
    S_3  = S*np.sqrt(S)
    C1   = ((d0 + d1*T + d2*T**2 + d3*T**3 + d4*T**4)*S + (e0 + e1*T + e2*T**2)*S_3)*P
    C2   = ((f0 + f1*T + f2*T**2 + f3*T**3)*S + g0*S_3)*P**2
    C3   = ((h0 + h1*T + h2*T**2)*S + j1*T*S_3)*P**3
    dCp  = C1 + C2 + C3
    return (Cp0 + dCp0 + dCp)

def compute_era5_radiation_surface_heat_flux(self,method='bilinear',periodic=True,reuse_weights=True,F_weights=''):
        '''
        '''
        # first re-grid
        G_CICE = self.cice_grid_preperation()
        stdts  = self.time_series_day_month_start()
        yr_str = self.define_year_string(stdts[0])
        #P_ERA5 = self.define_era5_variable_path(var_name=self.era5_base_var_name)
        G_ERA5 = self.era5_grid_prep()
        lw     = xesmf_regrid_dataset(G_ERA5,G_CICE,
                                        os.path.join(self.D_ERA5,'msdwlwrf',yr_str,'*.nc'),
                                        method=method,
                                        periodic=periodic,
                                        reuse_weights=reuse_weights,
                                        F_weights=self.F_ERA5_weights,
                                        multi_file=True)
        sw     = xesmf_regrid_dataset(G_ERA5,G_CICE,
                                        os.path.join(self.D_ERA5,'msdwswrf',yr_str,'*.nc'),
                                        method=method,
                                        periodic=periodic,
                                        reuse_weights=reuse_weights,
                                        F_weights=self.F_ERA5_weights,
                                        multi_file=True)
        msshf  = xesmf_regrid_dataset(G_ERA5,G_CICE,
                                        os.path.join(self.D_ERA5,'msshf',yr_str,'*.nc'),
                                        method=method,
                                        periodic=periodic,
                                        reuse_weights=reuse_weights,
                                        F_weights=self.F_ERA5_weights,
                                        multi_file=True)
        mslhf  = xesmf_regrid_dataset(G_ERA5,G_CICE,
                                        os.path.join(self.D_ERA5,'mslhf',yr_str,'*.nc'),
                                        method=method,
                                        periodic=periodic,
                                        reuse_weights=reuse_weights,
                                        F_weights=self.F_ERA5_weights,
                                        multi_file=True)
        F_net  = sw.msdwswrf - lw.msdwlwrf - mslhf.mslhf - msshf.msshf
        F_net  = F_net.assign_coords(lon=G_CICE.lon[0,:],lat=G_CICE.lat[:,0])
        F_net  = F_net.rename({'x':'lon','y':'lat'})
        F_net  = F_net.set_index(lon='lon',lat='lat')
        return F_net

dTdt   = temp.differentiate("time").sel(st_ocean=mld_win,method='nearest')
print("computing radiation surface heat flux from ERA5")
F_net  = self.compute_era5_radiation_surface_heat_flux()
F_net  = F_net.rename({'lat':'yt_ocean','lon':'xt_ocean'})
F_net  = F_net.sortby('yt_ocean').reindex()
qdp    = np.zeros( ( len(temp.time), len(temp.yt_ocean), len(temp.xt_ocean) ) )
for j in range(len(temp.time)):
    print("computing ocean heat flux at the mixed layer depth for ",temp.time[j].values)
    qdp[j,:,:] = compute_ocn_heat_flux_at_mixed_layer_depth(temp.isel(time=j),
                                                            salt.isel(time=j),
                                                            mld.isel(time=j),
                                                            F_net.isel(time=j),
                                                            dTdt.isel(time=j))