Page 1 of 1

Subroutine calc_coszen seems to have an offset of 1 day

PostPosted: Thu Sep 28, 2017 8:19 pm
by cvrodrigues
Dear WRF users and developers,

Subroutine calc_coszen seems to have an offset of 1 day.

In phys/module_radiation_driver.F there is subroutine calc_coszen.
This is called only two times in the code (in phys/module_radiation_driver.F):

Code: Select all
         CALL radconst(XTIME,DECLIN,SOLCON,JULIAN,               &
                       DEGRAD,DPD                                )
         call calc_coszen(ims,ime,jms,jme,its,ite,jts,jte, &
                          julian,xtime,gmt,declin,degrad,  &
                          xlong,xlat,coszen_loc,hrang_loc)
...
! added coszen subroutine : jararias, 2013/08/10
!   outputs: coszen, hrang
     call calc_coszen(ims,ime,jms,jme,its,ite,jts,jte,  &
                      julian,xtime+radt*0.5,gmt, &
                      declin,degrad,xlong,xlat,coszen,hrang)


As you may see the variables JULIAN, XTIME and GMT are arguments to the function.

Before proceeding some discussion is necessary. GMT is the time at the simulation start and together with JULDAY (the "day of the year" at simulation start) define the date & time at the beginning of the simulation. As such, both GMT and JULDAY are constants.

XTIME is the time in minutes since the simulation started, being updated as the simulation progresses (I presume at each time step). Now JULIAN is the date & time referred to the beginning of the year (at least when the simulation started). As such one way to compute JULIAN is:

Code: Select all
JULIAN = XTIME(:) / 1440. + GMT / 24. + JULDAY - 1


Mind that JULIAN is defined such that January 1st is 0, conversely to JULDAY that would have a value of 1 instead. As it can be observed, the argument to calc_coszen is JULIAN (not JULDAY). At the start of subroutine calc_coszen we have:

Code: Select all
   SUBROUTINE calc_coszen(ims,ime,jms,jme,its,ite,jts,jte,  &
                          julian,xtime,gmt, &
                          declin,degrad,xlon,xlat,coszen,hrang)
...
       real, intent(in)    :: julian,declin,xtime,gmt,degrad
...
       da=6.2831853071795862*(julian-1)/365.
...


Thus, da = 2 * pi * (JULIAN-1) / 365., meaning the current angular position relative to the Sun of Earth's translation. However (JULIAN-1)/365 will return a value in [-0.0027, 0.9945], instead of [0, 1]. Moreover it does not account for leap years.

If indeed this analysis is correct, a quick fix would be achieved doing:
Code: Select all
da = 2 * pi * JULIAN / (365. - 1)


Best regards,

C V Rodrigues

Re: Subroutine calc_coszen seems to have an offset of 1 day

PostPosted: Fri Sep 29, 2017 3:24 pm
by kwthomas
Hi,..

I looked at the 3.9.1.1 code and I agree with you.

Even the code documentation (second comment) seems strange.

! JULIAN starts from 0.0 at 0Z on 1 Jan.
intJULIAN = JULIAN + 1.0 ! offset by one day
! jan 1st 00z is julian=1.0 here
IJUL=INT(intJULIAN)

Send email to wrfhelp@ucar.edu and report what you observed.