Question about using skin temperature as sst input for noskin function in ECMWF algorithm

[hhweiCT]

Hi, thank you so much for providing this useful python package. I used noskin function and it works well. Since I would like to use the ECMWF algorithm, the instruction said that we should use the skin function. When I use skin function, it crashed without any error message. I can make it not crash by reducing the input array dimension size but the output value does not look reasonable.

I saw some discussion saying the skin function is not working yet. Therefore, I am wondering whether I can use the noskin function but use "skin temperature" as the sst input (instead of the bulk sst) for the ECMWF algorithm. I checked the code and it seems that the skin function would update the T_s from sst to skin temperature with the cold-skin+warm-layer scheme. Thanks!

[jbusecke]

Interesting question @hhweiCT. Re the skin function I will raise a detailed issue shortly.

I am concerned about this:

When I use skin function, it crashed without any error message.

This is not what I observed last I checked. Could you confirm your version? Are you running this in a notebook or script? The latter often provides better error messages so that might be worth a try.

I can make it not crash by reducing the input array dimension size

What sizes are we talking here? This might just run out of memory on your machine.

but the output value does not look reasonable.

Can you provide a minimal example which I could reproduce?

Assuming that we figure out what is leading to the crash I suppose you could actually see if your approach is good enough yourself: Compute some representative values for noskin using skin_sst and skin with bulk sst, then compare the results.

ccing @Brodeau who might have more insight into the actual algorithm and whether my naive suggestion is reasonable or not.

[hhweiCT]

Thank you so much! I will provide more information for the crash soon. (Also, try to use script again instead of notebook)

Yes, that is what I wanted to try too (compare "noskin using skin_sst" and "skin with bulk sst"). I will for sure do that after making the skin function work.

[hhweiCT]

This is my code I attach a shared link for the input data (only one time and over North Pacific), which are from the 6-hourly ERA-Interim data. https://drive.google.com/drive/folders/1-LuRfYtKhujrIEseSU3Gdpc_jMOo5Bz7?usp=sharing

from aerobulk import skin
from aerobulk import noskin

import xarray as xr
import numpy as np

ds_NPac_1979010100=xr.open_dataset('ds_NPac_alldata_1979010100.nc')

da_sst_Rechunked_NPac=ds_NPac_1979010100['da_sst_Rechunked_NPac']
da_Tskin_Rechunked_NPac=ds_NPac_1979010100['da_Tskin_Rechunked_NPac']
da_T2m_Rechunked_NPac=ds_NPac_1979010100['da_T2m_Rechunked_NPac']
da_uwnd_Rechunked_NPac=ds_NPac_1979010100['da_uwnd_Rechunked_NPac']
da_vwnd_Rechunked_NPac=ds_NPac_1979010100['da_vwnd_Rechunked_NPac']
da_q2m_Rechunked_NPac=ds_NPac_1979010100['da_q2m_Rechunked_NPac']
da_ps_NPac=ds_NPac_1979010100['da_ps_NPac']

ds_NPac_radiation_1979010100=xr.open_dataset('ds_NPac_radiation_1979010100.nc')

da_ssrd_NPac_Wm2=ds_NPac_radiation_1979010100['da_ssrd_NPac_Wm2']
da_strd_NPac_Wm2=ds_NPac_radiation_1979010100['da_strd_NPac_Wm2']

# using this range would trigger this error
 # *** E R R O R :
 # BULK_FORMULA_VCTR()@mod_phymbl: wind stress too strong!
 #  =>   836.18 N/m^2 ! At ji, jj = 0001, 0001

# llat1=0
# llat2=-1
# llon1=0
# llon2=-1

## using this range (chose randomly) would work but the output is not reasonable 
llat1=25
llat2=45
llon1=60
llon2=90

LH_skin,SH_skin,taux_skin,tauy_skin,T_s_output_skin,Evap_skin=skin(
       da_sst_Rechunked_NPac[llat1:llat2,llon1:llon2],
       da_T2m_Rechunked_NPac[llat1:llat2,llon1:llon2],
       da_q2m_Rechunked_NPac[llat1:llat2,llon1:llon2],
       da_uwnd_Rechunked_NPac[llat1:llat2,llon1:llon2],
       da_vwnd_Rechunked_NPac[llat1:llat2,llon1:llon2],
       da_ssrd_NPac_Wm2[llat1:llat2,llon1:llon2],
       da_strd_NPac_Wm2[llat1:llat2,llon1:llon2],
       slp=da_ps_NPac[llat1:llat2,llon1:llon2],
       algo='ecmwf',
       zt=2,
       zu=10,
       niter=6,
       input_range_check=False)

# LH_noskin,SH_noskin,taux_noskin,tauy_noskin,Evap_noskin=noskin(
#        da_sst_Rechunked_NPac[llat1:llat2,llon1:llon2],
#        da_T2m_Rechunked_NPac[llat1:llat2,llon1:llon2],
#        da_q2m_Rechunked_NPac[llat1:llat2,llon1:llon2],
#        da_uwnd_Rechunked_NPac[llat1:llat2,llon1:llon2],
#        da_vwnd_Rechunked_NPac[llat1:llat2,llon1:llon2],
#        slp=da_ps_NPac[llat1:llat2,llon1:llon2],
#        algo='ecmwf',
#        zt=2,
#        zu=10,
#        niter=6,
#        input_range_check=False)

As I wrote in the comments, if I use the whole North Pacific spatial range, I encounter error message (after using script, the error message did show up! thanks! )

*** E R R O R : BULK_FORMULA_VCTR()@mod_phymbl: wind stress too strong! => 836.18 N/m^2 ! At ji, jj = 0001, 0001

If I use a smaller spatial range, it does not trigger error but the output skin temperature (T_s_output_skin) is around 400K. The flux LH (order of 10^5) and SH (order of -10^3) are also not with reasonable values.

I am not sure where I did wrong.

I also check the code by running the noskin function with the same input, that one does output reasonable LH, SH values.

I tried to go into the Fortran code and it seems that the slp values when using the skin function is not correct (I tried to output the slp in the BULK_FORMULA, it is in the order of 100, instead of 10^5). I am not sure whether it is the reason. I also cannot figure out why the slp is different.

Thanks for your help!!