[MITgcm-support] Friction velocity ustar in S/R exf_wind.F

Mon Apr 23 12:58:37 EDT 2018

Hi Martin!

Yes, this possible issue may not affect the further calculations
because the call to s/r exf_wind (where ustar is calculated in such a
manner) preceeds the call to s/r exf_bulkformulae with correct ustar
calculation. Still, I found that (possible?) bug with ustar by
accident and decided to bring it to light, may be it can be important
in any cases.

I went to git repository, but the bug-fix procedure seems a little bit
too complicated...

Stanislav M.

> Hi Stanislav,

>I think you got me convinced that there is a problem. The code has been around for a long time (the routine was first introduced in 2006), but the code was there even earlier in the ecco branch. I can’t figure out who added the code originally, in order ask if we are overlooking >something, but maybe someone on this list feels responsible and can check …
>In a quick check, it seems that this does not affect any of the forward tests that we run regularily. That means that this code is probably not used very often, probably a reason why this did go unnoticed for such a long time.

>If you want to, you could contribute the bug-fix yourself to the still relatively new git repository, see here how to do that: <https://mitgcm.readthedocs.io/en/latest/contributing/contributing.html>

>Martin

> See further comments below:

> On 22. Apr 2018, at 15:04, Stanislav Martyanov <martyanov.sd at gmail.com> wrote:
>
> Hi Martin!
>
> Thank your for your reply. Yes, I understand that some optimization
> was made for the computational reasons. But I still have some doubts
> about the used formula for ustar in s/r exf_wind. I'll try to explain
> what I mean in detail.
>
> So, let's see the code in exf_wind.f with my brief comments:
>
> Line 146-153: We calculate usSq = (ustress*ustress) +
> (vstress*vstress). Both ustress and vstress have units [N/m^2] because
> they are the wind stress vector's components along X and Y axes. After
> that we obtain variable usSq which has unit [ (N/m^2)^2 ]. It is OK.
>
> Line 156: wStress = SQRT(usSq). So wStress is now in [N/m^2]. That is
> OK too, because wStress is just a module of wind vector.
>
This is already wrong, because usSq is stress squared. It should be
SQRT(SQRT(usSq))/atmrho). The rest are follow up errors.
> Line 157: c            ustar = SQRT(usSq/atmrho). This is a correct
> statement, but it is commented.
>
> Line 176: recip_sqrtRhoA = 1. _d 0 / SQRT(atmrho). That is OK, we are
> optimizing the code for the further computations.
>
> Line 179: ustar = wStress*recip_sqrtRhoA.    Now that is a problem, as
> I think. Because after that, ustar has units (N/m^2) / (kg/m^3)^0.5.
> But it has to have units as friction velocity - [m/s]. The problem is
> that we have not taken another one SQRT.
>
I agree that this should fix it:
> If we rewrite, adding additional SQRT, as:
> ustar = SQRT(wStress)*recip_sqrtRhoA,
>
> then the units of ustar will be [m/s] as it should be. Friction
> velocity is basically determined as SQRT(stress/density).
>
> Best regards,
> Stanislav

2018-04-23 16:30 GMT+03:00 Martin Losch <Martin.Losch at awi.de>:
> Hi Stanislav,
>
> I think you got me convinced that there is a problem. The code has been around for a long time (the routine was first introduced in 2006), but the code was there even earlier in the ecco branch. I can’t figure out who added the code originally, in order ask if we are overlooking something, but maybe someone on this list feels responsible and can check …
> In a quick check, it seems that this does not affect any of the forward tests that we run regularily. That means that this code is probably not used very often, probably a reason why this did go unnoticed for such a long time.
>
> If you want to, you could contribute the bug-fix yourself to the still relatively new git repository, see here how to do that: <https://mitgcm.readthedocs.io/en/latest/contributing/contributing.html>
>
> Martin
>
> See further comments below:
>
>> On 22. Apr 2018, at 15:04, Stanislav Martyanov <martyanov.sd at gmail.com> wrote:
>>
>> Hi Martin!
>>
>> Thank your for your reply. Yes, I understand that some optimization
>> was made for the computational reasons. But I still have some doubts
>> about the used formula for ustar in s/r exf_wind. I'll try to explain
>> what I mean in detail.
>>
>> So, let's see the code in exf_wind.f with my brief comments:
>>
>> Line 146-153: We calculate usSq = (ustress*ustress) +
>> (vstress*vstress). Both ustress and vstress have units [N/m^2] because
>> they are the wind stress vector's components along X and Y axes. After
>> that we obtain variable usSq which has unit [ (N/m^2)^2 ]. It is OK.
>>
>> Line 156: wStress = SQRT(usSq). So wStress is now in [N/m^2]. That is
>> OK too, because wStress is just a module of wind vector.
>>
> This is already wrong, because usSq is stress squared. It should be SQRT(SQRT(usSq))/atmrho). The rest are follow up errors.
>> Line 157: c            ustar = SQRT(usSq/atmrho). This is a correct
>> statement, but it is commented.
>>
>> Line 176: recip_sqrtRhoA = 1. _d 0 / SQRT(atmrho). That is OK, we are
>> optimizing the code for the further computations.
>>
>> Line 179: ustar = wStress*recip_sqrtRhoA.    Now that is a problem, as
>> I think. Because after that, ustar has units (N/m^2) / (kg/m^3)^0.5.
>> But it has to have units as friction velocity - [m/s]. The problem is
>> that we have not taken another one SQRT.
>>
> I agree that this should fix it:
>> If we rewrite, adding additional SQRT, as:
>> ustar = SQRT(wStress)*recip_sqrtRhoA,
>>
>> then the units of ustar will be [m/s] as it should be. Friction
>> velocity is basically determined as SQRT(stress/density).
>>
>> Best regards,
>> Stanislav
>>
>>
>>
>>
>>
>>> Hi Stanislav,
>>>
>>> I am no expert on this, but did you see lines 156/157:
>>>             wStress(i,j,bi,bj) = SQRT(usSq)
>>> c            ustar = SQRT(usSq/atmrho)
>>>
>>> so ustar should be sqrt(usSq/atmrho). for algorithmic reasons wstress=sqrt(usSq) is computed first and then divided by sqrt(atmrho) later in lin 179. To this makes sense, although it’s not very easy to read.
>>>
>>> Martin
>>
>>> On 4. Mar 2018, at 13:58, Stanislav Martyanov <martyanov.sd at gmail.com> wrote:
>>>
>>> Hello!
>>>
>>> I think it should be
>>>
>>> ustar = SQRT (wStress(i,j,bi,bj) / atmrho)
>>>
>>> rather than
>>>
>>> ustar = wStress(i,j,bi,bj)*recip_sqrtRhoA
>>>
>>> in exf package in S/R exf_winf.F (line 179)?
>>>
>>> In the current version wStress is not affected by the function SQRT,
>>> but should be, in my opinion.
>>>
>>> Best regards,
>>> Stanislav
>> _______________________________________________
>> MITgcm-support mailing list
>> MITgcm-support at mitgcm.org
>> http://mailman.mitgcm.org/mailman/listinfo/mitgcm-support
>
> _______________________________________________
> MITgcm-support mailing list
> MITgcm-support at mitgcm.org
> http://mailman.mitgcm.org/mailman/listinfo/mitgcm-support