[MITgcm-support] Friction velocity ustar in S/R exf_wind.F
Martin Losch
Martin.Losch at awi.de
Tue Apr 24 05:21:24 EDT 2018
In fact, I will (o:
Just to practice (or rather develop) my git skills ...
Martin
> On 23. Apr 2018, at 19:06, Jody Klymak <jklymak at uvic.ca> wrote:
>
>
>
>> On 23 Apr 2018, at 09:58, Stanislav Martyanov <martyanov.sd at gmail.com> wrote:
>> I went to git repository, but the bug-fix procedure seems a little bit
>> too complicated…
>
> If you don’t feel like making a full Pull Request, then simply open an issue, with as much specifics as you can, including line numbers and file names if possible. Then someone with a development environment can fix it for you.
>
> Cheers, Jody
>
>> 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
>> _______________________________________________
>> MITgcm-support mailing list
>> MITgcm-support at mitgcm.org
>> http://mailman.mitgcm.org/mailman/listinfo/mitgcm-support
>
> --
> Jody Klymak
> http://web.uvic.ca/~jklymak/
>
>
>
>
>
> _______________________________________________
> MITgcm-support mailing list
> MITgcm-support at mitgcm.org
> http://mailman.mitgcm.org/mailman/listinfo/mitgcm-support
More information about the MITgcm-support
mailing list