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
