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?

[DeniseWorthen]

C96.

[Qingfu-Liu]

I am not sure following test can fix the problem, but I think you can try it.
This problem you showed might only impact this low resolution run. I have the following reasoning:
In the code dcyc2t3.f line 285-286:
tem1 = fhswr / deltim
nstp = max(6, nint(tem1))
In this low resolution run: tem1 = fhswr / deltim =3600/720=5
If nint(tem1) is less than 6, the xcosz(i) may not calculate correctly.
You can change: nstp = max(5, nint(tem1))
Try to run it again. Not sure if this can fix this problem.

[DeniseWorthen]

@Qingfu-Liu Would it also work to switch dt_atmos=600 to test your dcyc2t3 change? I can do at run time w/o recompile.

[Qingfu-Liu]

Yes. It should work

[yangfanglin]

@Qingfu-Liu what happens if dt_atmos is 1800s ?

[Qingfu-Liu]

@yangfanglin I looked the code again, I am not sure if the calculations of nstp and nstl cause this problem. But the problem related to the calculation of xcosz(i) and coszen(i)

[junwang-noaa]

I was looking at the same code as Qingfu. I think the problem is coszen, which uses at least the 6 substeps within a SW call interval.

[DeniseWorthen]

Setting dt_atmos=600 (and also dt for ice and updating the coupling freq) results in the same odd pattern of both diffuse and direct. Diffuse N-IR shown here:

[Qingfu-Liu]

In our case, nstl=1, so only the first part of the code is used:
if (nstl == 1) then
cns = pid12 * (solhr + deltimf7200 - hour12) + slag
do i = 1, IM
xcosz(i) = sdecsinlat(i) + cdec*coslat(i)*cos(cns+xlon(i))
enddo
elseif (nstl == nstp) then

[Qingfu-Liu]

I am curious why the visible direct SW component and visible diffuse SW component show two different patterns? They are supposed to be multipled by the same weight

[DeniseWorthen]

@Qingfu-Liu That's a good question. I've shown just the visible, but the n-ir show the same difference between dir and dif.

[Qingfu-Liu]

So I think the problem may not come from the weighting xmu(i). I searched two variables sfcnirdfd and sfcvisbmd, there are no other physics programs using them (except in program dcyc2t3.f). Those two variables are carried around from other part of the code, and might be changed between the dynamic steps.

[junwang-noaa]

Do we see the same issue on higher resolution test? Does this also show up in model history tile file field such as sfc down sw?

[Qingfu-Liu]

GFS_surface_generic_post.F90

@DeniseWorthen You are running a couple test (the code you changed need to set the parameter cplflx=.true.)? Can you run an uncoupled test (set the parameter cplflx=.false.) to see if we still have the same problem? (I am still learning how to output those results you showed here)

[junwang-noaa]

Yes, at this level dt_phys is actually the atmosphere time step and is applied to both dycore and phys. Inside physics, it can have its inner time step defined for physics schemes (e.g. dt_inner for Thompson MP).

[Qingfu-Liu]

@junwang-noaa I looked the code fv3atm_history_io.F90, still not sure where the DSWRF_AVE is calculated. So If I want to print out the DSWRF_AVE, which part of the code I should add the printout? Thanks

[junwang-noaa]

You can check the code here:

https://github.com/NOAA-EMC/fv3atm/blob/develop/io/fv3atm_history_io.F90#L450

[Qingfu-Liu]

@junwang-noaa Thank you very much. Some of the code is hard to understand (for example, why the albedo line 361 is multiplied by the lcnvfac). I will add printout to look of the output data

[DeniseWorthen]

@Qingfu-Liu @junwang-noaa I spent some time trying to understand what is happening w/ xmu. I understand now why you say that it is inaccurate, but not a bug.

The part I'm still having trouble understanding is why using the estimated SW for time-steps inside of fswhr is better than just allowing the SW to be a constant, or maybe just adjusting at the mid-point of fswhr. What I see is just a very noisy SW signal. I don't understand the advantage of that.

[Qingfu-Liu]

@DeniseWorthen Allowing the SW to be a constant is OK if the radiation time step is short (say less than 1 hour). However, in earlier numerical modeling, radiation time step is 3 hours (you can found that ECMWF model still use 3 hours several years ago), so interpolation to the dynamic time step using xmu(i) is a better way comparing to use as a constant, especially in the early morning and late afternoon when the radiation flux changes with time is very large. The accuracy of xmu(i) has been improvement, but the improvement has not been added into GFS model yet. See the paper:
Manners, J., Thelen, J.-C., Petch, J., Hill, P. and Edwards, J. M. (2009). Two fast radiative transfer methods
to improve the temporal sampling of clouds in numerical weather prediction and climate models. Q. J. R.
Meteorol. Soc., 135, 457–468.

[Qingfu-Liu]

@DeniseWorthen If you have time, could you please run a test and plot the variables DSWRF_AVE and ATMIMP_FAXA_SWNDF, that will be great help. The new branch here:
https://github.com/Qingfu-Liu/fv3atm/tree/HR2_fluxr_fix
If you are too busy, just neglect this message. Thanks

[DeniseWorthen]

I had done a run but forgot about it. If I can't find it, I'll rerun.

[DeniseWorthen]

@Qingfu-Liu Interesting--the paper you point to references on page 4 a 2016 paper:

The cosine of the solar zenith angle is also passed to the radiation scheme, having been computed as the average value over the sunlit part of the radiation time step as described by Hogan and Hirahara (2016)

And that 2016 paper actually references the paper I posted about above, by Zhang (2015) for calculating the mean zenith angle

I was also looking at this paper which seems somewhat related (?) since it is talking about how to estimate cosz between radiation timesteps.

The code implementing Zhang's 2015 method I found in the CDEPS share code. It's an interesting problem!

[Qingfu-Liu]

@junwang-noaa Considering the large variations of DSWRF_AVE values, I feel that the DSWRF_AVE output does not make much sense, not sure where it is used. I give an example here: assuming we have 6 hours accumulation bucket, DSWRF instantaneous values has a parabolic shape with time. Assuming we have dynamic time step 10 minutes, when we calculating the DSWRF_AVE from 6AM to 11:50AM, we will have DSWRF_AVE values gradually increase with time. However, at 12PM, we will suddenly have a peak value equals the 12PM DSWRF instant value (which is very large compared to the previous 6h average value), then gradually decrease with time. Where we can use those kind of weird shape of DSWRF_AVE ? By the way, any variables with large diurnal cycles will have the same strange results. If those average variables just used for diagnoses, I think we should remove them.

[junwang-noaa]

@Qingfu-Liu This is the bucket average field that GFS have been outputting for a long time. The feature you mentioned is shown in these fields. I am including @WenMeng-NOAA @yangfanglin for their comments of downstream users with these fields from GFS output.

[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-NOAA]

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-NOAA]

@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-NOAA]

@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-NOAA]

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-NOAA]

@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-NOAA]

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!