Cosine solar zenith angle forcing is off by a time-step

[ekluzek]

Moving this over from the cime issue:

ESMCI/cime#3380

From some analysis by Sean Swenson, we see that the cosine solar zenith angle forcing is off by a time-step.

I think the problem is in shr_tInterpMod.F90 in shr_tInterp_getAvgCosz, where a time update happens and then the solar-zenith angle is calculated. I think that should be reversed the solar zenith angle is calculated and then the time-update happens.

This affects CLM standalone simulations for all forcing types (I compsets). But, also CORE_IAF and CORE2_IAF forcing (JRA and IAF C and G compsets as well as the DIAF compset

[mvertens]

@ekluzek @swensosc - can you please create a PR with the relevant changes. We can then create new baselines and get this in.

[jedwards4b]

@swensosc @ekluzek is this what you need?

    do while( reday < reday2) ! mid-interval t-steps thru interval [LB,UB]    
       !--- get next cosz value for t-avg ---                                                                                                                                                                                                  
       call shr_tInterp_getCosz(cosz,lon,lat,ymd,tod,eccen,mvelpp,lambm0,obliqr,calendar)                                                                                                                                                      
       n = n + ldt                                                                                                                                                                                                                             
       tavCosz = tavCosz + cosz*real(ldt,r8)  ! add to partial sum                                                                                                                                                                                                     
                                                                                                                                                                                                                                               
       !--- advance to next time in [LB,UB] ---                                                                                                                                                                                                
       ymd0 = ymd                                                                                                                                                                                                                              
       tod0 = tod                                                                                                                                                                                                                              
       reday0 = reday                                                                                                                                                                                                                          
       call shr_cal_advDateInt(ldt,'seconds',ymd0,tod0,ymd,tod,calendar)                                                                                                                                                                       
       call shr_cal_timeSet(reday,ymd,tod,calendar)                                                                                                                                                                                            
                                                                                                                                                                                                                                               
       if (reday > reday2) then                                                                                                                                                                                                                
          ymd = ymd2                                                                                                                                                                                                                           
          tod = tod2                                                                                                                                                                                                                           
          timeint = reday2-reday0                                                                                                                                                                                                              
          call ESMF_TimeIntervalGet(timeint, s_i8=dtsec, rc=rc)                                                                                                                                                                                
          if (ChkErr(rc,__LINE__,u_FILE_u)) return                                                                                                                                                                                             
          ldt = int(dtsec, shr_kind_in)                                                                                                                                                                                                        
       endif                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                           
    end do                                          

[swensosc]

yes, the advance should occur after the getCosz call. In this order, does the 'if (reday > reday2) then' actually do anything? It looks like it would exit the while loop immediately after it would potentially be triggered (?)

…

On Thu, Sep 30, 2021 at 10:30 AM Jim Edwards ***@***.***> wrote: @swensosc <https://github.com/swensosc> @ekluzek <https://github.com/ekluzek> is this what you need? do while( reday < reday2) ! mid-interval t-steps thru interval [LB,UB] !--- get next cosz value for t-avg --- call shr_tInterp_getCosz(cosz,lon,lat,ymd,tod,eccen,mvelpp,lambm0,obliqr,calendar) n = n + ldt tavCosz = tavCosz + cosz*real(ldt,r8) ! add to partial sum !--- advance to next time in [LB,UB] --- ymd0 = ymd tod0 = tod reday0 = reday call shr_cal_advDateInt(ldt,'seconds',ymd0,tod0,ymd,tod,calendar) call shr_cal_timeSet(reday,ymd,tod,calendar) if (reday > reday2) then ymd = ymd2 tod = tod2 timeint = reday2-reday0 call ESMF_TimeIntervalGet(timeint, s_i8=dtsec, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return ldt = int(dtsec, shr_kind_in) endif end do — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#120 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AGRN57GTDOJK4YCISOUCEFTUESGD5ANCNFSM5EM3WSTA> . Triage notifications on the go with GitHub Mobile for iOS <https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675> or Android <https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub>.

[jedwards4b]

I agree, that block of code does nothing:

   do while( reday < reday2) ! mid-interval t-steps thru interval [LB,UB]                                                                                                                                                                     
                                                                                                                                                                                                                                               
       !--- get next cosz value for t-avg ---                                                                                                                                                                                                  
       call shr_tInterp_getCosz(cosz,lon,lat,ymd,tod,eccen,mvelpp,lambm0,obliqr,calendar)                                                                                                                                                      
       n = n + ldt                                                                                                                                                                                                                             
       tavCosz = tavCosz + cosz*real(ldt,r8)  ! add to partial sum                                                                                                                                                                                                                                                                                                                                                                                                                 

       !--- advance to next time in [LB,UB] ---                                                                                                                                                                                                
       ymd0 = ymd                                                                                                                                                                                                                              
       tod0 = tod                                                                                                                                                                                                                                                                                                                                                                                                                                                
       call shr_cal_advDateInt(ldt,'seconds',ymd0,tod0,ymd,tod,calendar)                                                                                                                                                                       
       call shr_cal_timeSet(reday,ymd,tod,calendar)                                                                                                                                                                                            
                                                                                                                                                                                                                                               
    end do                                                                                                                                                                                                                                     
    tavCosz = tavCosz/real(n,r8) ! form t-avg                                                                                                                                                                                                  

[jedwards4b]

Can you provide a test that it's now correct?

[swensosc]

I will have to go back and see how I checked this; it has been a couple of years. I believe that what I did was run the model a few days with time step level output, then check that the interpolated solar was more smoothly varying at the forcing time boundaries. Would something like that be a sufficient test?

…

On Thu, Sep 30, 2021 at 10:59 AM Jim Edwards ***@***.***> wrote: Can you provide a test that it's now correct? — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#120 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AGRN57DJRXU4GRUKDDRTRFLUESJQRANCNFSM5EM3WSTA> . Triage notifications on the go with GitHub Mobile for iOS <https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675> or Android <https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub>.

[jedwards4b]

Sure - here is the cdeps fork https://github.com/jedwards4b/CDEPS and branch jedwards/solar_zenith_angle_correction

[swensosc]

I think that the other change, in cdeps/streams/dshr_strdata_mod.F90, needs to be made at the same time as this. Iirc, the error in dshr_tinterp_mod.F90 was probably made in order to compensate for the error in dshr_strdata_mod.F90. That bug in dshr_strdata_mod.F90 is in this line: if (rDateM < rDateLB .or. rDateM > rDateUB) then which should be: if (rDateM < rDateLB .or. rDateM >= rDateUB) then (the > should be a >=)

…

On Thu, Sep 30, 2021 at 11:18 AM Jim Edwards ***@***.***> wrote: Sure - here is the cdeps fork https://github.com/jedwards4b/CDEPS and branch jedwards/solar_zenith_angle_correction — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#120 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AGRN57FHUBCFWZ2VZSDACALUESLYDANCNFSM5EM3WSTA> . Triage notifications on the go with GitHub Mobile for iOS <https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675> or Android <https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub>.

[jedwards4b]

I've pushed this additional change to the branch.

[jedwards4b]

From @swensosc "I ran two five day simulations, with and without your cdeps changes. The changes correct the bug."

[mvertens]

Fixed by #123