ufs-community / ufs-weather-model

UFS Weather Model
Other
137 stars 247 forks source link

SW net to ocean is too large at low sun angles #1768

Closed DeniseWorthen closed 1 year ago

DeniseWorthen commented 1 year ago

Description

The SW net sent to the ocean in CMEPS uses a constant ocean albedo of 0.06 for the ocean fraction. This replicated the existing NEMS mediator functionality at the time we switched to CMEPS, but results in too much SW net to the ocean when the albedo > 0.06 due to a low sun angle.

This is an issue with both the DATM application as well as the fully coupled. It is not an issue for the HAFS model because HAFS is able to use SWnet from the ATM since there is no time-varying ice fraction.

To Reproduce:

Additional context

When running fully coupled, the SW net to the ocean depends on both the time-varying ice-fraction, which determines the SW penetrating through the sea ice into the ocean, as well the SW net at the open ocean surface. Since the ocean runs at a slower coupling frequency than ATM-ICE, this is a mean value over the slow coupling frequency. That is

SWnet to ocean = mean (swnet ocean * ofrac) + mean (ifrac * sw penetration)

Currently, in UFS the swup over the ocean = (1-.06) * sw down. At low sun angles, when albedo>0.06, this over-estimates the open water SWnet.

Note that we cannot use just the SWnet provided by the ATM because we need to weight by the time-varying ice fraction. Otherwise, the value sent to the ocean would be

SWnet to ocean = mean(swnet) * mean(ofrac) + mean(ifrac) * mean(sw penetration)

To correctly account for both sun angle and ice-fraction, we need to calculate the ocean albedo at the end of each fast coupling step and then accumulate the net SW from both ATM and ICE. The ocean then receives the mean net SW to the ocean at its coupling frequency

The code to implement an ocean albedo phase is already present in the CMEPS code. It estimates ocean albedo as

                ocnalb%anidr(n) = (.026_r8/(cosz**1.7_r8 + 0.065_r8)) +   &
                                  (.150_r8*(cosz         - 0.100_r8 ) *   &
                                  (cosz - 0.500_r8 ) * (cosz - 1.000_r8 )  )
                ocnalb%avsdr(n) = ocnalb%anidr(n)
                ocnalb%anidf(n) = albdif
                ocnalb%avsdf(n) = albdif

with albdif= 0.06 and albdir = 0.07.

This appears to be the same formula used in ccpp/physics/physics/GFS_surface_generic_post.F90, which I think is where the albedo is being calculated for the global case.

              ocalnirbm_cpl = max(albdf, 0.026_kind_phys/(xcosz_loc**1.7_kind_phys+0.065_kind_phys)     &
       &                       + 0.15_kind_phys * (xcosz_loc-0.1_kind_phys) * (xcosz_loc-0.5_kind_phys) &
       &                       * (xcosz_loc-one))

with albdf = 0.06 (It isn't clear why the bm albedo is using the diffuse value here as a max, and not the direct value as a max).

The use of the max function is not present in the current CMEPS albedo code, but can easily be set as a configurable option so that the code still works for the CESM case.

The ocean albedo calculation occurs with in a new ocnalb run phase which occurs after the ATM and ICE both run and before the fields are accumulated for the ocean. The SWnet for the open ocean is then

SW net to ocean = Sw from atm (1-ifrac) (1-albedo) + ifrac * SW penetration

This should result in a lower SW net to ocean in the mean.

I will post figures for a coupled and DATM case in the next few days.

Output

Qingfu-Liu commented 1 year ago

I searched online, in the CESM code seq_flux_mct.F90, there is no max function in the calculation of the ocean surface albedo (not sure the code has been changed). May we can remove the max function so that the calculation in ccpp/physics/physics/GFS_surface_generic_post.F90 to be consistent with those in CMEPS

(line 876 - 901 in seq_flux_mct.F90):

! Compute albedos do n=1,nloc_o rlat = const_deg2rad lats(n) rlon = const_deg2rad lons(n) cosz = shr_orb_cosz( calday, rlat, rlon, delta ) if (cosz > 0.0_r8) then !--- sun hit -- anidr = (.026_r8/(cosz*1.7_r8 + 0.065_r8)) + & (.150_r8(cosz - 0.100_r8 ) & (cosz - 0.500_r8 ) & (cosz - 1.000_r8 ) ) avsdr = anidr anidf = seq_flux_mct_albdif avsdf = seq_flux_mct_albdif else !--- dark side of earth --- anidr = 1.0_r8 avsdr = 1.0_r8 anidf = 1.0_r8 avsdf = 1.0_r8 end if xao_o%rAttr(index_xao_So_avsdr,n) = avsdr xao_o%rAttr(index_xao_So_anidr,n) = anidr xao_o%rAttr(index_xao_So_avsdf,n) = avsdf xao_o%rAttr(index_xao_So_anidf,n) = anidf end do ! nloc_o

DeniseWorthen commented 1 year ago

@Qingfu-Liu The MCT code is not used by CMEPS (MCT is the previous mediator). The ocean albedo calculation is in med_phases_ocnalb_mod.F90. In my working branch (PR #1773), I have:

          ! Compute albedos
          do n = 1,lsize
             rlat = const_deg2rad * ocnalb%lats(n)
             rlon = const_deg2rad * ocnalb%lons(n)
             cosz = shr_orb_cosz( nextsw_cday, rlat, rlon, delta )
             if (cosz  >  0.0_r8) then !--- sun hit --
                ocnalb%anidr(n) = (.026_r8/(cosz**1.7_r8 + 0.065_r8)) +   &
                                  (.150_r8*(cosz         - 0.100_r8 ) *   &
                                  (cosz - 0.500_r8 ) * (cosz - 1.000_r8 )  )
                if (use_min_albedo) then
                   !TODO: why does fv3atm use albdif here and not albdir ?
                   ocnalb%anidr(n) = max (ocnalb%anidr(n), albdif)
                end if
                ocnalb%avsdr(n) = ocnalb%anidr(n)
                ocnalb%anidf(n) = albdif
                ocnalb%avsdf(n) = albdif
             else !--- dark side of earth ---
                ocnalb%anidr(n) = 1.0_r8
                ocnalb%avsdr(n) = 1.0_r8
                ocnalb%anidf(n) = 1.0_r8
                ocnalb%avsdf(n) = 1.0_r8
             end if
          end do

One question you could help answer for me is why, in ccpp/physics/physics/GFS_surface_generic_post.F90, the max beam (ie, direct) albedo uses albdf?

              ocalnirbm_cpl = max(albdf, 0.026_kind_phys/(xcosz_loc**1.7_kind_phys+0.065_kind_phys)     &
       &                       + 0.15_kind_phys * (xcosz_loc-0.1_kind_phys) * (xcosz_loc-0.5_kind_phys) &
       &                       * (xcosz_loc-one))

In the ocnalb phase, the direct and diffuse have different values

    real(R8), parameter     :: albdif = 0.06_r8 ! 60 deg reference albedo, diffuse
    real(R8), parameter     :: albdir = 0.07_r8 ! 60 deg reference albedo, direct

Is the use of max(albdf,...) in the physic routine correct?

Qingfu-Liu commented 1 year ago

@Denise Worthen - NOAA Affiliate @.***> Both variables albdif and albdir are constants. I think someone might want to use a smaller value (0.06 instead of 0.07) to limit the ocean surface albedo. So I think the code should be fine.

On Wed, May 31, 2023 at 8:53 AM Denise Worthen @.***> wrote:

@Qingfu-Liu https://github.com/Qingfu-Liu The MCT code is not used by CMEPS (MCT is the previous mediator). The ocean albedo calculation is in med_phases_ocnalb_mod.F90. In my working branch (PR #1773 https://github.com/ufs-community/ufs-weather-model/pull/1773), I have:

      ! Compute albedos
      do n = 1,lsize
         rlat = const_deg2rad * ocnalb%lats(n)
         rlon = const_deg2rad * ocnalb%lons(n)
         cosz = shr_orb_cosz( nextsw_cday, rlat, rlon, delta )
         if (cosz  >  0.0_r8) then !--- sun hit --
            ocnalb%anidr(n) = (.026_r8/(cosz**1.7_r8 + 0.065_r8)) +   &
                              (.150_r8*(cosz         - 0.100_r8 ) *   &
                              (cosz - 0.500_r8 ) * (cosz - 1.000_r8 )  )
            if (use_min_albedo) then
               !TODO: why does fv3atm use albdif here and not albdir ?
               ocnalb%anidr(n) = max (ocnalb%anidr(n), albdif)
            end if
            ocnalb%avsdr(n) = ocnalb%anidr(n)
            ocnalb%anidf(n) = albdif
            ocnalb%avsdf(n) = albdif
         else !--- dark side of earth ---
            ocnalb%anidr(n) = 1.0_r8
            ocnalb%avsdr(n) = 1.0_r8
            ocnalb%anidf(n) = 1.0_r8
            ocnalb%avsdf(n) = 1.0_r8
         end if
      end do

One question you could help answer for me is why, in ccpp/physics/physics/GFS_surface_generic_post.F90, the max beam (ie, direct) albedo uses albdf?

          ocalnirbm_cpl = max(albdf, 0.026_kind_phys/(xcosz_loc**1.7_kind_phys+0.065_kind_phys)     &
   &                       + 0.15_kind_phys * (xcosz_loc-0.1_kind_phys) * (xcosz_loc-0.5_kind_phys) &
   &                       * (xcosz_loc-one))

In the ocnalb phase, the direct and diffuse have different values

real(R8), parameter     :: albdif = 0.06_r8 ! 60 deg reference albedo, diffuse
real(R8), parameter     :: albdir = 0.07_r8 ! 60 deg reference albedo, direct

Is the use of max(albdf,...) in the physic routine correct?

— Reply to this email directly, view it on GitHub https://github.com/ufs-community/ufs-weather-model/issues/1768#issuecomment-1570178625, or unsubscribe https://github.com/notifications/unsubscribe-auth/AGTS6UUJMNP4NPRN3WICAHDXI45LPANCNFSM6AAAAAAYN2ZWUA . You are receiving this because you were mentioned.Message ID: @.***>

DeniseWorthen commented 1 year ago

Thanks, I realize these are constants, I'm just curious as to why the lower albedo was set. It seems very much a way to tune the response, which may or may not be appropriate when we run the coupled model vs standalone.

A second issue is that w/in the ocnalb routine, we have the capability to set the albedos as just averages, similar to what is being done now. But that feature uses the two different values of constant albedo (0.06, 0.07).

Qingfu-Liu commented 1 year ago

@Denise Worthen - NOAA Affiliate @.***> Agree with you. I think the max function should be removed, or give a very small number (much smaller than 0.06).

On Wed, May 31, 2023 at 9:31 AM Denise Worthen @.***> wrote:

Thanks, I realize these are constants, I'm just curious as to why the lower albedo was set. It seems very much a way to tune the response, which may or may not be appropriate when we run the coupled model vs standalone.

— Reply to this email directly, view it on GitHub https://github.com/ufs-community/ufs-weather-model/issues/1768#issuecomment-1570242853, or unsubscribe https://github.com/notifications/unsubscribe-auth/AGTS6UU4YJCO3BEK3Q4HZO3XI5B3VANCNFSM6AAAAAAYN2ZWUA . You are receiving this because you were mentioned.Message ID: @.***>

yangfanglin commented 1 year ago

Is it possible to trace the history of ocalnirbm_cpl and check who coded this at the first place and ask the question why the max was used ?

Qingfu-Liu commented 1 year ago

@Fanglin Yang - NOAA Federal @.> I traced back and found that the section of the code first appeared in Aug 2018 ( grantfirl https://github.com/ufs-community/ccpp-physics/commits?author=grantfirl committed on Aug 8, 2018 ). The code might imported from other place. I will add Grant to this email. @Grant Firl - NOAA Affiliate @.> Hi Grant, do you know why we have the limit in the calculation of Ocean surface albedo? Thanks ocalnirbm_cpl = max(albdf, 0.026/(xcosz_loc*1.7+0.065) & & + 0.15 (xcosz_loc-0.1) (xcosz_loc-0.5) & & (xcosz_loc-1.0))

On Wed, May 31, 2023 at 10:38 AM Fanglin Yang @.***> wrote:

Is it possible to trace the history of ocalnirbm_cpl and check who coded this at the first place and ask the question why the max was used ?

— Reply to this email directly, view it on GitHub https://github.com/ufs-community/ufs-weather-model/issues/1768#issuecomment-1570365852, or unsubscribe https://github.com/notifications/unsubscribe-auth/AGTS6URO25KMXYT3V63RAY3XI5JXRANCNFSM6AAAAAAYN2ZWUA . You are receiving this because you were mentioned.Message ID: @.***>

DeniseWorthen commented 1 year ago

The max function will mean that the ocean absorbs at a maximum 94% of the downward SW radiation, ie (1-.06)*dswrf, right? I'm guessing it was put in to reduce ocean warming. But the fix I'm proposing using the cmeps ocean albedo calculation also has a cooling effect. So removing the max function along with that might result in a cold bias. I don't think there is any way to know, short of doing some test runs.

jiandewang commented 1 year ago

The max function will mean that the ocean absorbs at a maximum 94% of the downward SW radiation, ie (1-.06)*dswrf, right? I'm guessing it was put in to reduce ocean warming. But the fix I'm proposing using the cmeps ocean albedo calculation also has a cooling effect. So removing the max function along with that might result in a cold bias. I don't think there is any way to know, short of doing some test runs.

we can make some realtive long runs (i.e HR1 or even back to P8 using g-w) to assess the impact by the time

DeniseWorthen commented 1 year ago

@LydiaStefanova-NOAA Prepared some figures for me using a 20-day run of the cpld_control_nowave_noaero_p8 RT.

These figures show area mean (30S:30N) comparisons of the net SW in the ocean output (ocnSW) and SST between the case using the ocean albedo calculation (fix) vs the current constant (0.06) ocean albedo (nofix).

Screen Shot 2023-06-02 at 8 52 20 AM

The figures below show the imbalance between the netSW that the ocean receives and the netSW that the ATM has internally (using it's own time-varying ocean albedo). Ideally, the netSW received by the ocean should match the netSW that the ATM has at the ocean surface. Using the current constant albedo in calculating the netSW to the ocean produces about a 1-2 W/m2 excess SW to the ocean.

Screen Shot 2023-06-02 at 9 00 02 AM
DeniseWorthen commented 1 year ago

@denise Worthen - NOAA Affiliate @.***> Agree with you. I think the max function should be removed, or give a very small number (much smaller than 0.06).

@Qingfu-Liu To retain current default behaviour, perhaps there could be a namelist setting for whether to apply this maximum or not and maybe even a way to choose what the minimum value is (0, .06, .07 etc). That way it would be easy to match what the ATM is doing and what the CMEPS ocean albedo routine is doing.

junwang-noaa commented 1 year ago

@DeniseWorthen From the top right plot on SST, the SST is dropping with the fix instead of warming up the ocean, is this what we expected? Also you mentioned the fix you have in CMEPS also has a cooling effect.

The max function will mean that the ocean absorbs at a maximum 94% of the downward SW radiation, ie (1-.06)*dswrf, right? I'm guessing it was put in to reduce ocean warming. But the fix I'm proposing using the cmeps ocean albedo calculation also has a cooling effect. So removing the max function along with that might result in a cold bias. I don't think there is any way to know, short of doing some test runs.

DeniseWorthen commented 1 year ago

@junwang-noaa I did expect a cooling effect from using the CMEPS mediator to calculate the ocean albedo. I have not changed anything in FV3 at this point. In fact, both CMEPS and FV3 at this point are setting the max direct albedo = max(0.06, f(cosz)). If we removed the max(), then when the sun is high, the albedo can be lower than 0.06 and my thinking would be that would have an opposite warming effect.

DeniseWorthen commented 1 year ago

This is the daily mean SST difference (swfix-nofix) at the end of the 20day test case

Screen Shot 2023-06-05 at 11 42 25 AM
Qingfu-Liu commented 1 year ago

@DeniseWorthen if you change the 0.06 to 0.00, do you think we can remove some of the large negative bias in the SST between 30S-30N in the long time forecast?

DeniseWorthen commented 1 year ago

@Qingfu-Liu Certainly the new case is cooler (on average). But this is a single, low-resolution run. Wouldn't you need multiple longer runs which are compared to obs to say whether it is biased?

Or are you just interested in how sensitive this case would be to removing the 0.06? I could do that test, removing the 0.06 from both CMEPS and FV3.

jiandewang commented 1 year ago

@Qingfu-Liu Certainly the new case is cooler (on average). But this is a single, low-resolution run. Wouldn't you need multiple longer runs which are compared to obs to say whether it is biased?

Or are you just interested in how sensitive this case would be to removing the 0.06? I could do that test, removing the 0.06 from both CMEPS and FV3.

I can make some long runs (repeat P8 or HR1 so that we have something to compare) using g-w if needed

Qingfu-Liu commented 1 year ago

This comment moved to issues #1767, and posted to wrong place

DeniseWorthen commented 1 year ago

@Qingfu-Liu The issue of the imbalance between the ATM's internal SWnet over ocean and what the OCN receives as the SWnet is fundamentally about the ocean albedo that CMEPS is using when exporting the netSW to the ocean. It is not directly related to the issue of how the ATM is calculating the SW down at each timestep (Issue #1767).

Qingfu-Liu commented 1 year ago

@DeniseWorthen @jiandewang I think we should run the low resolution (and high resolution) tests to see the SST bias for long time forecast, then try to do research on what should be changed in order to remove the SST bias.

LydiaStefanova-NOAA commented 1 year ago

@Qingfu-Liu I think the main point for this particular issue would be to make sure that the calculations are done so that there is no mismatch between the atm and ocn SW. What impact that ends up having on the SST bias, in my opinion, would be a separate question.

jiandewang commented 1 year ago

@Qingfu-Liu I think you are mixing three separated issue here:

  1. SW inconsistency bwtween ATM and OCN
  2. calculation of DSW in ATM itslef
  3. SST bias in UFS

for (1) @DeniseWorthen has identified that the problem happened in CEMP. The receiving of DSW in CMEP from ATM is correct. It is the CMEP that didn't translate the correct one to OCN. The fixing will require some long run to assess the impact and I will do that by the time. for (2) it is purely in ATM side, and we will re-do the assessment of impact in long runs. for (3) it is due to the whole physical package (including the calculation of DSW in ATM) in UFS

these three are kind of independent. @DeniseWorthen @LydiaStefanova-NOAA and I will try to solve (1) at this moment and we can work with you to repeat the long run when you have new methid in (2)

Qingfu-Liu commented 1 year ago

@jiandewang @LydiaStefanova-NOAA The calculation of DSW in ATM itself is not accurate, but it is not a bug. The improvement will take time (I do not expect the new version will be ready for HR2).

jiandewang commented 1 year ago

@Qingfu-Liu got you, let's try to slove item 1 first.

grantfirl commented 1 year ago

@fanglin Yang - NOAA Federal @.> I traced back and found that the section of the code first appeared in Aug 2018 ( grantfirl https://github.com/ufs-community/ccpp-physics/commits?author=grantfirl committed on Aug 8, 2018 ). The code might imported from other place. I will add Grant to this email. @grant Firl - NOAA Affiliate @.> Hi Grant, do you know why we have the limit in the calculation of Ocean surface albedo? Thanks ocalnirbm_cpl = max(albdf, 0.026/(xcosz_loc1.7+0.065) & & + 0.15 (xcosz_loc-0.1) (xcosz_loc-0.5) & & (xcosz_loc-1.0)) On Wed, May 31, 2023 at 10:38 AM Fanglin Yang @.> wrote: Is it possible to trace the history of ocalnirbm_cpl and check who coded this at the first place and ask the question why the max was used ? — Reply to this email directly, view it on GitHub <#1768 (comment)>, or unsubscribe https://github.com/notifications/unsubscribe-auth/AGTS6URO25KMXYT3V63RAY3XI5JXRANCNFSM6AAAAAAYN2ZWUA . You are receiving this because you were mentioned.Message ID: @.***>

@Qingfu-Liu If my name is on a commit from that time, it would have likely been the initial work of putting this code into a CCPP scheme from the original GFS physics driver. I have not had any input on the actual algorithm. Also, be sure to reference @grantfirl instead of the other GitHub user name (the one that thumbs downed your comment) since that is not me.

Qingfu-Liu commented 1 year ago

@Grant Firl - NOAA Affiliate @.***> Sorry to mess up your user name. I did not pay detailed attention to the user name when I quoted the name. I will double check the user name next time. Thank you very much.

@grantfirl Thanks for remind me your user name. I really did not pay attention to it. You have the same name. I just use the @ and the user name pop up. Qingfu

On Tue, Jun 6, 2023 at 3:04 PM Grant Firl @.***> wrote:

@fanglin https://github.com/fanglin Yang - NOAA Federal @.> I traced back and found that the section of the code first appeared in Aug 2018 ( grantfirl https://github.com/ufs-community/ccpp-physics/commits?author=grantfirl https://github.com/ufs-community/ccpp-physics/commits?author=grantfirl committed on Aug 8, 2018 ). The code might imported from other place. I will add Grant to this email. @grant https://github.com/grant Firl - NOAA Affiliate @.> Hi Grant, do you know why we have the limit in the calculation of Ocean surface albedo? Thanks ocalnirbm_cpl = max(albdf, 0.026/(xcosz_loc

1.7+0.065) & & + 0.15 (xcosz_loc-0.1) (xcosz_loc-0.5) & & (xcoszloc-1.0)) … <#m-2080750619557003295_> On Wed, May 31, 2023 at 10:38 AM Fanglin Yang @.**> wrote: Is it possible to trace the history of ocalnirbm_cpl and check who coded this at the first place and ask the question why the max was used ? — Reply to this email directly, view it on GitHub <#1768 (comment) https://github.com/ufs-community/ufs-weather-model/issues/1768#issuecomment-1570365852>, or unsubscribe https://github.com/notifications/unsubscribe-auth/AGTS6URO25KMXYT3V63RAY3XI5JXRANCNFSM6AAAAAAYN2ZWUA https://github.com/notifications/unsubscribe-auth/AGTS6URO25KMXYT3V63RAY3XI5JXRANCNFSM6AAAAAAYN2ZWUA . You are receiving this because you were mentioned.Message ID: @.***>

@Qingfu-Liu https://github.com/Qingfu-Liu If my name is on a commit from that time, it would have likely been the initial work of putting this code into a CCPP scheme from the original GFS physics driver. I have not had any input on the actual algorithm. Also, be sure to reference @grantfirl https://github.com/grantfirl instead of the other GitHub user name (the one that thumbs downed your comment) since that is not me.

— Reply to this email directly, view it on GitHub https://github.com/ufs-community/ufs-weather-model/issues/1768#issuecomment-1579300992, or unsubscribe https://github.com/notifications/unsubscribe-auth/AGTS6UXKUQTWCZHUALIXMYTXJ55KTANCNFSM6AAAAAAYN2ZWUA . You are receiving this because you were mentioned.Message ID: @.***>

jkbk2004 commented 1 year ago

@jieshunzhu FYI about this SST cold bias discussion

DeniseWorthen commented 1 year ago

The current implementation I have working allows for verification against the existing SWnet calculation for the CDEPS and coupled tests.

1) There is already an existing configuration within the ocnalb calculation to use fixed albedos in the SWnet calculation, triggered by setting flux_albav = .true. This sets albdif = 0.06 and albdir = 0.07. I've added the ability to instead set albdif = albdir = 0.06 via configuration. This essentially replicates (although not B4B) the current calculation of SWnet exported to the ocean.

2) For the coupled case, a new configuration variable which sets the ocean albedo limit ocean_albedo_limit = 0.06 can be used to force the ocean albedo calculation to do what FV3 is currently doing when it applies the max function to the calculated albedo

ocalnirbm_cpl = max(albdf, 0.026_kind_phys/(xcosz_loc**1.7_kind_phys+0.065_kind_phys)     &
       &                       + 0.15_kind_phys * (xcosz_loc-0.1_kind_phys) * (xcosz_loc-0.5_kind_phys) &
       &                       * (xcosz_loc-one))

Unless an optional configuration for ocean albedo is introduced in FV3, the coupled tests will require the use of this albedo limit for consistency between ATM and OCN.

The control GEFS CDEPS test was run (24hr) for the current develop branch and two cases using the ocnalb calculation. The upper panel shows the SST difference when albdir=albdif=0.06 is used. This essentially replicates the current behavior in CMEPS, which does not use the ocean albedo calculation. The lower panel shows the difference when the ocean albedo is calculated as a function of time-of-day. In this case, no limit was placed on the albedo.

Screen Shot 2023-06-12 at 11 39 24 AM
jiandewang commented 1 year ago

@DeniseWorthen looks we got cooling everything (as expected). Let me know when time is mature for me to make some long runs.

DeniseWorthen commented 1 year ago

@jiandewang For this DATM case, the results are pretty unambiguous. With the coupled model, it takes more time for a consistent pattern to develop. I posted an SST difference after 20days for the coupled case last week.

Hopefully PR #1752 will be merged by tomorrow; after that I'll be able to provide a testing branch with the required nems.configure changes for you to test.

DeniseWorthen commented 1 year ago

@jiandewang You should be able to test this issue using my feature/fixsw branch. You'll need to use this branch because it contains some additional share code and configuration variables that are needed for the ocean albedo calculation. You'll be able to change the nems.configure in order to test make a direct comparison between using the average albedos (equivalent to what we are currently doing) versus calculating the ocean albedos using the coszen.

In the MED_attributes of nems.configure use

      ocean_albedo_limit = 0.06

This will calculate the ocean albedos but enforce a minimum value of 0.06, which is equivalent to what the ATM is doing.

To replicate what we are currently doing, use

      flux_albav = .true.
      albdir = 0.06

This will use constant albedos of 0.06 for both the direct and diffuse components.

jiandewang commented 1 year ago

@DeniseWorthen thanks for the information. Do you have a run directory on disk for me to take a look on key files such as input.nml, et. al ?

DeniseWorthen commented 1 year ago

I have a benchmark RT run directory on hera

/scratch1/NCEPDEV/stmp2/Denise.Worthen/FV3_RT/rt_12453/cpld_bmark_p8_intel
jiandewang commented 1 year ago

@DeniseWorthen I have made a set of mini runs (8 cases) for SW fixed vs non-fixed using your feature/fixsw branch in g-w. Can you double check the nems.configure I used. They are on HERA /scratch1/NCEPDEV/stmp2/Jiande.Wang/S2S_RT/SW-problem/CONFIG

@LydiaStefanova-NOAA the results are saved on HPSS at /NCEPDEV/emc-climate/1year/Jiande.Wang/HERA/SW-problem note ocn-ice-grib2-tar have some missing files (~ 4 files each) for all cases. Apparently there are some imperfectness in g-w script. The runs used P8 IC with C384+OCN025+ICE025.

the way I did is checkout the latest g-w develop branch (I did this about 5 days ago, so it's not what you will see today's commit), then replaced the UWM with Denise's branch. I found it's impossible to just simply repeat P8 or HR1 but replace CMEP with Denise's branch as latest CMEPS branch must come with latest WW.

DeniseWorthen commented 1 year ago

@jiandewang Yes, those two configurations are correct. The "nofix" case replicates the current code since it uses average albedos for both direct and diffuse and sets both averages to 0.06. The "fix" case uses the ocean albedo calculation but limits the calculated direct albedo to be at least 0.06 (consistent w/ the ATM).

LydiaStefanova-NOAA commented 1 year ago

Looking at the set of 8 P8 runs with/without the fix, the results show that:

a) The proposed fix noticeably improves the agreement between ATM and OCN net SW radiation b) The impact of the fix on SST is a slight cooling of the domain-average (consistent with expectations), but the magnitude of this cooling is quite small on the time scales of up to 1 month. The impact on SST bias appears very small.

The plots on which these statements are based can be seen here: https://docs.google.com/presentation/d/1bwUgW_HCfg9Tp-kXzzz7J4JI2-k4zd5D4V0o70Enb9k/edit?usp=sharing

junwang-noaa commented 1 year ago

@LydiaStefanova-NOAA Thanks for verifying the results. I am curious from which files you are looking at the atm/ocn radiation fields. Are the ATM net SW radiation fields from the atm sfc files or they are from the mediator files? Also where are the OCN net SW radiation fields from?

LydiaStefanova-NOAA commented 1 year ago

@junwang-noaa I used the 0.25 deg grib2 files for both ATM (from gfs.t00z.pgrb2.0p25. ) and OCN (from ocn_ice_0p25x0p25.grb2)

DeniseWorthen commented 1 year ago

@LydiaStefanova-NOAA @jiandewang Thanks for doing those test runs and analyzing the results. I'm going to make a few more changes to the CMEPS PR after consultation w/ ESCOMP and then plan to open the PR #1773 for commit.

This PR will implement the ocean albedo calculation in the mediator but it will contain via config the ability to revert to the current SWnet calculation (ie, constant, equal ocean albedos).

DeniseWorthen commented 1 year ago

@jiandewang @LydiaStefanova-NOAA @JessicaMeixner-NOAA Because the changes to make this fix will also impact code that CESM uses, my plan is to first make the PR to ESCOMP from my feature branch of CMEPS. I want to make sure this fix is OK w/ the coupled group first before PRing ESCOMP. Note that committing this change will allow UFS to continue using albdir=albdif=0.06 always, if desired. But the answers will not be B4B with current code.

jiandewang commented 1 year ago

@jiandewang @LydiaStefanova-NOAA @JessicaMeixner-NOAA Because the changes to make this fix will also impact code that CESM uses, my plan is to first make the PR to ESCOMP from my feature branch of CMEPS. I want to make sure this fix is OK w/ the coupled group first before PRing ESCOMP. Note that committing this change will allow UFS to continue using albdir=albdif=0.06 always, if desired. But the answers will not be B4B with current code.

@DeniseWorthen sounds good for me

LydiaStefanova-NOAA commented 1 year ago

@DeniseWorthen sounds good to me too

yangfanglin commented 1 year ago

@DeniseWorthen Denise, to confirm, this update has been committed to ESCOMP and/or ufs-weather model. The upcoming HR3 will be using it without taking any extra action, right ?

DeniseWorthen commented 1 year ago

@yangfanglin Yes, this PR is complete. It has been merged to ESCOMP and is currently being used in all the cpld RT tests for UFS. I had created a G-W issue for making the changes there NOAA-EMC/global-workflow/issues/1811

junwang-noaa commented 1 year ago

@yangfanglin The GW PR#1889 will bring the changes to GW.