In coupled configs, SW exported by ATM has odd time-dependence

[DeniseWorthen]

Description

The SW direct components exported by FV3 have odd "blips" over the course of a day. This is shown by setting mediator history files to write the instantaneous fields every time through the coupling loop.

Note that "instantaneous" in this case refers to the type of mediator history file that CMEPS is writing out and not the type of field (instantaneous or mean) field that ATM is providing.

For the C96 case on Tile 1, at i=50,j=70 (6E, 0.4N) the visible direct SW component is shown below (near ir shows similar). Every instantaneous value is shown, with values on the hours called out.

The visible diffuse component shows

I do not believe this is due to any cloud interaction, since the blips occur on the hour and are nearly symmetric during the diurnal cycle.

To Reproduce:

Mediator history files for the atmosphere can be added using the following in MED_attributes in nems.configure:

      history_tile_atm = 96
      history_n_atm_inst = 1
      history_option_atm_inst = nsteps

Files are produced for every coupling timestep and can be cat'd together to form a timeseries using nco. The SW fields imported by the mediator are named

near-ir diffuse: atmImp_Faxa_swndf
visible diffuse: atmImp_Faxa_swvdf
near-ir direct: atmImp_Faxa_swndr
visible direct: atmImp_Faxa_swvdr

Additional context

Output

[yangfanglin]

@DeniseWorthen Denise, are these visible direct and diffuse downward SW fluxes at the surface ? Would be interesting to see if the current operational model (GFSv16) has the same issue. @Qingfu-Liu @dustinswales

[DeniseWorthen]

Yes, these are down, not net sw components (or they should be). They are in terms of the export fields

  ! MEAN sfc downward nir direct flux (W/m**2)
  case ('mean_down_sw_ir_dir_flx')
   call block_data_copy(datar82d, GFS_data(nb)%coupling%dnirbm_cpl, Atm_block, nb, scale_factor=rtime, rc=localrc)

I did a test swapping these with the "instantaneous" values GFS_data(nb)%coupling%dnirbmi_cpl and removing rtime and got an identical result.

[yangfanglin]

@DeniseWorthen I understand solar zenith is not a standard output, would it be possible to write out solar zenith angle at each time step to a text file and inspect if there is any odd features. A short run for 24 hours should suffice. I understand you are busy. Let me know if this will take too much of your time. I can ask physics developers to check. Thanks for uncovering this "bug".

[DeniseWorthen]

@yangfanglin If I knew how to translate the i,j location on a tile to the loop in the code where it is calculating cosz then I could do that. I think even just the single location above would show if that is the issue. Is that easy to do?

[yangfanglin]

@DeniseWorthen I think you can pick a random single point (i,j) in the cosz loop and print out the values, assuming you will have the luck of not picking a polar night point.

[junwang-noaa]

l think we have looked the zenith angle issue with the inst sw radiation and fixed it in GFSv16. I checked the code updates on zenith angle, this one comes out:

ufs-community/ccpp-physics@7644b0c#diff-18710705427e7b03b02adfa6bc3969fa8a5484bd18d6c08192fe7001e489f885R890

@grantfirl Could this be related to the zenith angle change Denise showed?

[yangfanglin]

@junwang-noaa See slide 5 of this ppt for the bug fix we made to IPD physics in 2019. The pattern looks to be slightly different from what Denise showed.

The change made to ufs-community/ccpp-physics@7644b0c#diff-18710705427e7b03b02adfa6bc3969fa8a5484bd18d6c08192fe7001e489f885R890 does beg for questions. coszen=0 is physical.

[DeniseWorthen]

That figure looks to be the combined diffuse + direct, so to me it looks like that is the basic same pattern.

[junwang-noaa]

Thanks, Fanglin. What I mean is that the bug fix in GFSv16 resolved the zigzag pattern shown in sw flux. The fixes in GFSv16 are in the current develop branch. I am trying to find changes coming to the model since GFSv16 that may cause the issue in the develop branch now.

[DeniseWorthen]

I was able to produce a cosz plot using a trick to substitute the value in the dnirbm_cpl array

diff --git a/physics/GFS_surface_generic_post.F90 b/physics/GFS_surface_generic_post.F90
index 76d3f570..c7a732cd 100644
--- a/physics/GFS_surface_generic_post.F90
+++ b/physics/GFS_surface_generic_post.F90
@@ -143,7 +143,12 @@
             dnirdfi_cpl (i) = adjnirdfd(i)
             dvisbmi_cpl (i) = adjvisbmd(i)
             dvisdfi_cpl (i) = adjvisdfd(i)
-            dnirbm_cpl  (i) = dnirbm_cpl(i) + adjnirbmd(i)*dtf
+            !dnirbm_cpl  (i) = dnirbm_cpl(i) + adjnirbmd(i)*dtf
+!            if (wet(i)) then                    ! some open water
+!  ---  compute open water albedo
+              xcosz_loc = max( zero, min( one, xcosz(i) ))
+              dnirbm_cpl  (i) = xcosz_loc
+!            endif

I'm not sure of the units, but a plot of the values makes sense in terms of shape for it to be cosz in the export field (I mulitiplied by 1e3 in this plot)

Plotting this field as a time series at the same i=50,j=70 on tile 1 gives

So I don't think it is a cosz issue.

Also, the original point I plotted on tile1 turns out to be over land. But an adjacent point on the ocean (i=48,j=70) shows the same zig-zag effect

[junwang-noaa]

@DeniseWorthen Is the zenith angle plot from every time step? The sw is called every hour, and the radiation heating is then adjusted at every time step with the zenith angle. From the first two plots in the issue, it looks to me the sw radiation flux values coming out of sw radiation are correct, but the value interpolated to each time step is not. Can you print out each time step value and the hourly values? Also can you make a run without the change in:

ufs-community/ccpp-physics@7644b0c#diff-18710705427e7b03b02adfa6bc3969fa8a5484bd18d6c08192fe7001e489f885R890

to see if there is any impact on xcosz?

[DeniseWorthen]

@junwang-noaa This is the cosz every time the ATM passes fields to the mediator. In the C96 configuration, that is every 720s.

I wasn't sure if radiation was calculated every hour or every 30mins, but I also had the same thought. The actual calculation of the radiation on the hour appears correct, but the interpolation to each time step is not.

I'll try to make the test you suggest.

[junwang-noaa]

I just checked the code, the change will impact coszen (average of cosine of zenith angle over daytime shortwave call time interval), not the xcosz, but it will impact the interpolation factor xmu. So you may want to plot a SW field such as the visible direct SW to see the impact.

[DeniseWorthen]

@junwang-noaa I made the following change to radiation_astronomy.f

!  --- ...  compute time averages

      do i = 1, IM
        coszdg(i) = coszen(i) * rstp
        if (istsun(i) > 0) coszen(i) = coszen(i) / istsun(i)
        !if (istsun(i) > 0 .and. coszen(i) /= 0.0_kind_phys) then
        !  coszen(i) = coszen(i) / istsun(i)
        !endif
      enddo

But I still see the same effect as on the original plots (both direct and diffuse)

[junwang-noaa]

The plots you show are from C96 or C384?

[Qingfu-Liu]

@junwang-noaa
continue our yesterday's discussion here regarding to the hourly output of the variable DSWRF_AVE.
After the code changes, the values of the hourly output for DSWRF_AVE will be different compared to those from our operational GFS. This difference can be explained here:
suppose we have 6 hour bucket, and 1 hour radiation time step and 10 minutes dynamic time step. The bucket start accumulation from 6AM. At 6AM, we have the same DSWRF_AVE output after the code change (comparing to the control), however, at 7AM, we will have different values of DSWRF_AVE (so different output). The calculations of the DSWRF_AVE at 7 AM are different
(in the following calculation, 60 is 60 minutes, 10 is 10 minutes, DSWRF_6h and DSWRF_7h are instantaneous values at 6AM and 7AM, respectively):

1. In the control (or current GFS), DSWRF_AVE calculation will be (60 is 60 minutes, 10 is 10 minutes):
   DSWRF_AVE = (DSWRF_6h * 60+DSWRF_7h * 60)/(60+60) = 0.5 * (DSWRF_6h + DSWRF_7h)
2. After the code changes,
   DSWRF_AVE = (DSWRF_6h * 60+DSWRF_7h * 10)/(60+10)
   So the hourly output of DSWRF_AVE will be different after the code changes.
   Since we only have hourly outputs for those radiation related averages, we can postpone the fix until we all agree with a solution.

[junwang-noaa]

@Qingfu-Liu Thanks for the explanation. I want to clarify that the sw is called after the hour, not at the end of the hour. In other words, the DSWRF_7h is output at fh=7h in sfcf007, at which time the sw is not called yet, it is called in the first time step integration after 7. The DSWRF is accumulated from 6 steps at 10mins after 7, 20mins after 7 until fh=8. In the two cases, we should have:
DSWRF_AVE = (DSWRF_6h * 60)/60
At fh=8.
DSWRF_AVE = (DSWRF_6h * 60+DSWRF_7h * 60)/(60+60)
if the model ends at fh=12, then we have:
DSWRF_AVE = (DSWRF_6h * 60+DSWRF_7h * 60+...DSWRF_11f*60)/(360)

[Qingfu-Liu]

@junwang-noaa Thanks for the clarification

[LarissaReames]

I wanted to check if this issue disappeared when I used a non-coupled ATM-only case by comparing output in the sfcf*.nc files produced by ATM between tests cpld_control_nowave_noaero_p8_intel and control_c96 ( I just swapped out C48 files with C96). The results, shown below (dashed lines are column-avg cloud cover to show clear-sky points) from a point off the NE coast of he US (329,50) show that this issue is only present in the coupled simulation.

Note that this issue is also present for average DSWRF at the top of the atmosphere

But not for instantaneous DSWRF at the surface

Again, these are ATM files, not mediator files, since CMEPS isn't used in control_c96. I also, just as a double check test, ran control_C96 with the executable from pld_control_nowave_noaero_p8_intel and got identical results. I also ran a coupled c192 regression test to check resolution dependence and saw similar behavior to cpld_control_nowave_noaero_p8_intel

[junwang-noaa]

@LarissaReames-NOAA Thank you for looking into this issue! It looks to me the DSWRF_ave fields are correct in control_C96. Something happened in C96 coupled test.

[DeniseWorthen]

To see what is exported by the ATM in the coupled configs, you can also use the DumpFields setting in the ATM attributes. I think that feature was broken at the time this issue was first created, but Dusan fixed that it #2355.

The export state shortwave values are what the ATM is actually exporting at each step through the coupling sequence. Whether that is the same field that gets put into the sfc files is something to be determined.

[LarissaReames]

@DeniseWorthen Thanks for pointing out that feature. I didn't know about it.

I ran another run of the cpld_control_nowave_noaero_p8_intel test with DumpFields = true in ATM and do indeed see similar patterns in diffuse (inst_down_sw_vis_dif_flx) and direct (inst_down_sw_vis_dir_flx) SW flux in the fv3_cap_export* files as you reported in your initial post.

However, it looks like all of the fields in the fv3_cap_import* files that come from the OCN component are all 0s. The fields that come from CICE look fine. Is this expected?

[yangfanglin]

So to confirm the issue reported in this thread has already been resolved, right ?

[LarissaReames]

@yangfanglin No, it's still very much an issue in coupled simulations. It just isn't (and was probably never) present in ATM-only forecasts.

[LarissaReames]

I ran another global test, this time with ATM+WAV, which means radiation variables aren't even exchanged and aren't in either ATM fv3_cap_import or export files. The undesirable "peaking" behavior is still present in dswrf_ave and dswrf_avetoa.

[dustinswales]

@LarissaReames-NOAA It's possible that the problem originates somewhere in here.
All of the instantaneous fluxes that are used for coupling use the (rad-2-phys) time-adjusted surface fluxes that come from dcyc2t3.f
The average fluxes are being integrated over the interval and multiplied by the dynamics timestep, which I'm not sure is correct?

[LarissaReames]

@dustinswales That would make sense if it were only the surface fluxes that display the issue and if the issue only appeared when cplflx=true or cpllnd=true. But, we also see the hourly fluctuations in the _avgtoa variables, and these variations are also present in the ATM+WAV simulation when no fluxes are coupled.

I mentioned this to @junwang-noaa yesterday evening and she suggested I try a ATM+AERO (GOCART) simulation as it doesn't use CMEPS. The hourly fluctuations in the averaged variables of the sfcf*.nc files are still present in these simulations.

[DeniseWorthen]

@LarissaReames-NOAA From what I remember (it has been a while), this is fundamentally an issue of the coszen approximation. I think what I did is to fill one of the diffuse arrays w/ the coszen angle --- this will fill the export field with the cosz value instead of the diffuse SW and you'll be able to plot the time-dependence.

diff --git a/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_surface_generic_post.F90 b/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_surface_generic_post.F90
index 648f6bf8..3091efc6 100644
--- a/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_surface_generic_post.F90
+++ b/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_surface_generic_post.F90
@@ -188,7 +188,8 @@
               nnirbmi_cpl(i) = adjnirbmd(i) * (one-ocalnirbm_cpl)
               nnirdfi_cpl(i) = adjnirdfd(i) * (one-ocalnirdf_cpl)
               nvisbmi_cpl(i) = adjvisbmd(i) * (one-ocalvisbm_cpl)
-              nvisdfi_cpl(i) = adjvisdfd(i) * (one-ocalvisdf_cpl)
+              !nvisdfi_cpl(i) = adjvisdfd(i) * (one-ocalvisdf_cpl)
+              nvisdfi_cpl(i) = xcosz_loc
             else
               nnirbmi_cpl(i) = adjnirbmd(i) - adjnirbmu(i)
               nnirdfi_cpl(i) = adjnirdfd(i) - adjnirdfu(i)

EDIT: aplogies..looking back at earlier messages, it is the variable xmu(i) which is the issue.

[LarissaReames]

So, after investigating, it appears the issue I was seeing in the sfcf* files is caused by setting fhzero=6 in the coupled simulations but fhzero=1 in the ATM-only simulations. This changes how often the diag buckets are emptied and causes a mis-calculation of the radiation variables here if you're writing more frequently than fhzero. This explains the fact that dswrf_ave is identical to that produced by ATM-only simulations in the first hour because rtime_radsw is less than rtime_int (which gets reset at every fhzero) in that first hour for both simulations, but for the next 5 hours is greater when fhzero = 6 so you get a slowly decreasing lcnvfrac. It also explains why every 6th hour of the simulation once again looks identical to the ATM only runs with fhzero=1.

You can see in this plot what happens to both dswrf_avg and dswrf_avgtoa when I switch fhzero from 6 to 1 in the ATM+AERO simulations.

Apologies for going down this rabbit hole but the cause of the issue in sfcf files is unrelated to the problem in the mediator files* and just an out-of-design use of high-frequency output with low-frequency diagnostic buckets resets since if you only looked at output every 6 hours then the values would be correct. I think this also means that we won't be able to use ATM only forecasts to help diagnose this issue since you can't get mediator files with these simulations and the instantaneous radiation fields in sfcf* files look fine regardless of settings.

I'm going to look more closely at what @DeniseWorthen pointed to in her latest comment. Thanks Denise for bring that back to the forefront of the discussion.

[dustinswales]

Good detective work!