From 8abe7d8e77259b4c71384ed57495043e8beab19c Mon Sep 17 00:00:00 2001 From: Jim Edwards Date: Wed, 15 Apr 2015 08:00:17 -0600 Subject: [PATCH 1/2] new optional solar insolation calculation (shr_orb_cosz) - From Minghua Zhang's email: A new routine (shr_orb_avg_cosz) has been introduced "that eliminates the spurious zonal oscillation of daily insolation caused by discrete time sampling. The algorithm uses the mean of cosz in a time step. The module shr_orb_mod.F90 calls this subroutine" - Added a shr_orb_init routine which is used to set the dt to use in the above calculation. If shr_orb_init is not called, the original shr_orb_cosz calculation is used. - It is important to note that ALL CESM components will use this modified cosz calculation with the specified dt value if it is turned on. --- share/csm_share/shr/shr_orb_mod.F90 | 162 ++++++++++++++++++++++++++-- 1 file changed, 156 insertions(+), 6 deletions(-) diff --git a/share/csm_share/shr/shr_orb_mod.F90 b/share/csm_share/shr/shr_orb_mod.F90 index 7fda841942de..ec9c99672ec5 100644 --- a/share/csm_share/shr/shr_orb_mod.F90 +++ b/share/csm_share/shr/shr_orb_mod.F90 @@ -1,11 +1,11 @@ !=============================================================================== -! SVN $Id: shr_orb_mod.F90 60325 2014-05-16 20:33:31Z santos@ucar.edu $ -! SVN $URL: https://svn-ccsm-models.cgd.ucar.edu/csm_share/trunk_tags/share3_150116/shr/shr_orb_mod.F90 $ +! SVN $Id: shr_orb_mod.F90 25434 2010-11-04 22:46:24Z tcraig $ +! SVN $URL: https://svn-ccsm-models.cgd.ucar.edu/csm_share/release_tags/cesm1_2_x_n02_share3_130715/shr/shr_orb_mod.F90 $ !=============================================================================== MODULE shr_orb_mod - use shr_kind_mod + use shr_kind_mod, only: SHR_KIND_R8, SHR_KIND_IN use shr_sys_mod use shr_const_mod use shr_log_mod, only: s_loglev => shr_log_Level @@ -20,6 +20,7 @@ MODULE shr_orb_mod public :: shr_orb_params public :: shr_orb_decl public :: shr_orb_print + public :: shr_orb_init real (SHR_KIND_R8),public,parameter :: SHR_ORB_UNDEF_REAL = 1.e36_SHR_KIND_R8 ! undefined real integer(SHR_KIND_IN),public,parameter :: SHR_ORB_UNDEF_INT = 2000000000 ! undefined int @@ -37,11 +38,24 @@ MODULE shr_orb_mod real (SHR_KIND_R8),parameter :: SHR_ORB_MVELP_MIN = 0.0_SHR_KIND_R8 ! min value for mvelp real (SHR_KIND_R8),parameter :: SHR_ORB_MVELP_MAX = 360.0_SHR_KIND_R8 ! max value for mvelp + ! dt for averaged cosz calculation + real (SHR_KIND_R8) :: dt_avg=-1 + !=============================================================================== CONTAINS !=============================================================================== -real(SHR_KIND_R8) FUNCTION shr_orb_cosz(jday,lat,lon,declin) +subroutine shr_orb_init(dt_in) + + ! Set the dt for the average cosz calculation if it is supplied + real(SHR_KIND_R8), intent(in) :: dt_in + + dt_avg = dt_in + +end subroutine shr_orb_init + + +real(SHR_KIND_R8) pure FUNCTION shr_orb_cosz(jday,lat,lon,declin) !---------------------------------------------------------------------------- ! @@ -63,11 +77,147 @@ real(SHR_KIND_R8) FUNCTION shr_orb_cosz(jday,lat,lon,declin) !---------------------------------------------------------------------------- - shr_orb_cosz = sin(lat)*sin(declin) - & - & cos(lat)*cos(declin)*cos(jday*2.0_SHR_KIND_R8*pi + lon) + ! If dt for the average cosz is specified, then call the shr_orb_avg_cosz + if (dt_avg == -1) then + shr_orb_cosz = sin(lat)*sin(declin) - & + & cos(lat)*cos(declin)*cos(jday*2.0_SHR_KIND_R8*pi + lon) + else + shr_orb_cosz = shr_orb_avg_cosz(jday, lat, lon, declin) + end if END FUNCTION shr_orb_cosz +!======================================================================= +! A New Algorithm for Calculation of Cosine Solar Zenith Angle +! Author: Linjiong Zhou +! E-mail: linjiongzhou@hotmail.com +! Date : 2015.02.22 +! Ref. : Zhou et al., GRL, 2015 +!======================================================================= + +real (SHR_KIND_R8) pure function shr_orb_avg_cosz(jday, lat, lon, declin) + + use shr_const_mod, only : pi => shr_const_pi + + implicit none + +!----------------------------------------------------------------------- +! In/Out Arguements + + real(SHR_KIND_R8), intent(in) :: jday ! Julian calendar day (1.xx to 365.xx) + real(SHR_KIND_R8), intent(in) :: lat ! latitude (radian) + real(SHR_KIND_R8), intent(in) :: lon ! longitude (radian) + real(SHR_KIND_R8), intent(in) :: declin ! solar declination (radian) + +!----------------------------------------------------------------------- +! Local Arguments + + real(SHR_KIND_R8),parameter :: piover2 = pi/2.0_SHR_KIND_R8 + real(SHR_KIND_R8),parameter :: twopi = pi*2.0_SHR_KIND_R8 + + real(SHR_KIND_R8) :: aa, bb + real(SHR_KIND_R8) :: del, phi + real(SHR_KIND_R8) :: cos_h, h + real(SHR_KIND_R8) :: t1, t2, dt + real(SHR_KIND_R8) :: tt1, tt2, tt3, tt4 + +!----------------------------------------------------------------------- +! Compute Half-day Length + + ! adjust latitude so that its tangent will be defined + if (lat == piover2) then + del = lat - 1.0e-05_SHR_KIND_R8 + else if (lat == -piover2) then + del = lat + 1.0e-05_SHR_KIND_R8 + else + del = lat + end if + + ! adjust declination so that its tangent will be defined + if (declin == piover2) then + phi = declin - 1.0e-05_SHR_KIND_R8 + else if (declin == -piover2) then + phi = declin + 1.0e-05_SHR_KIND_R8 + else + phi = declin + end if + + ! define the cosine of the half-day length + ! adjust for cases of all daylight or all night + cos_h = - tan(del) * tan(phi) + if (cos_h <= -1.0_SHR_KIND_R8) then + h = pi + else if (cos_h >= 1.0_SHR_KIND_R8) then + h = 0.0_SHR_KIND_R8 + else + h = acos(cos_h) + end if + +!----------------------------------------------------------------------- +! Define Local Time t and t + dt + + ! adjust t to be between -pi and pi + t1 = (jday - int(jday)) * twopi + lon - pi + + if (t1 >= pi) then + t1 = t1 - twopi + else if (t1 < -pi) then + t1 = t1 + twopi + end if + + dt = dt_avg / 86400.0_SHR_KIND_R8 * twopi + t2 = t1 + dt + +!----------------------------------------------------------------------- +! Compute Cosine Solar Zenith angle + + ! define terms needed in the cosine zenith angle equation + aa = sin(lat) * sin(declin) + bb = cos(lat) * cos(declin) + + ! define the hour angle + ! force it to be between -h and h + ! consider the situation when the night period is too short + if (t2 >= pi .and. t1 <= pi .and. pi - h <= dt) then + tt2 = h + tt1 = min(max(t1, -h) , h) + tt4 = min(max(t2, twopi - h), twopi + h) + tt3 = twopi - h + else if (t2 >= -pi .and. t1 <= -pi .and. pi - h <= dt) then + tt2 = - twopi + h + tt1 = min(max(t1, -twopi - h), -twopi + h) + tt4 = min(max(t2, -h) , h) + tt3 = -h + else + if (t2 > pi) then + tt2 = min(max(t2 - twopi, -h), h) + else if (t2 < - pi) then + tt2 = min(max(t2 + twopi, -h), h) + else + tt2 = min(max(t2 , -h), h) + end if + if (t1 > pi) then + tt1 = min(max(t1 - twopi, -h), h) + else if (t1 < - pi) then + tt1 = min(max(t1 + twopi, -h), h) + else + tt1 = min(max(t1 , -h), h) + end if + tt4 = 0.0_SHR_KIND_R8 + tt3 = 0.0_SHR_KIND_R8 + end if + + ! perform a time integration to obtain cosz if desired + ! output is valid over the period from t to t + dt + if (tt2 > tt1 .or. tt4 > tt3) then + shr_orb_avg_cosz = (aa * (tt2 - tt1) + bb * (sin(tt2) - sin(tt1))) / dt + & + (aa * (tt4 - tt3) + bb * (sin(tt4) - sin(tt3))) / dt + else + shr_orb_avg_cosz = 0.0_SHR_KIND_R8 + end if + +end function shr_orb_avg_cosz + !=============================================================================== SUBROUTINE shr_orb_params( iyear_AD , eccen , obliq , mvelp , & From c18f254bc38bd85a59df1a8645cb61524eccd1d2 Mon Sep 17 00:00:00 2001 From: Jim Edwards Date: Tue, 21 Apr 2015 07:59:00 -0600 Subject: [PATCH 2/2] modify when updated solar insolation is used - Added dt_avg as an optional parameter to shr_orb_cosz. If it is not present or if the dt_avg == 0, then the old calculation is used. Otherwise, the dt_avg is supplied to the new calculation. - Removed shr_orb_init - Each CESM component will have control over whether the modified cosz routine is used and will supply the time step to be used in the new calculation --- share/csm_share/shr/shr_orb_mod.F90 | 33 +++++++++++++---------------- 1 file changed, 15 insertions(+), 18 deletions(-) diff --git a/share/csm_share/shr/shr_orb_mod.F90 b/share/csm_share/shr/shr_orb_mod.F90 index ec9c99672ec5..1e3da9c993ab 100644 --- a/share/csm_share/shr/shr_orb_mod.F90 +++ b/share/csm_share/shr/shr_orb_mod.F90 @@ -20,7 +20,6 @@ MODULE shr_orb_mod public :: shr_orb_params public :: shr_orb_decl public :: shr_orb_print - public :: shr_orb_init real (SHR_KIND_R8),public,parameter :: SHR_ORB_UNDEF_REAL = 1.e36_SHR_KIND_R8 ! undefined real integer(SHR_KIND_IN),public,parameter :: SHR_ORB_UNDEF_INT = 2000000000 ! undefined int @@ -38,24 +37,12 @@ MODULE shr_orb_mod real (SHR_KIND_R8),parameter :: SHR_ORB_MVELP_MIN = 0.0_SHR_KIND_R8 ! min value for mvelp real (SHR_KIND_R8),parameter :: SHR_ORB_MVELP_MAX = 360.0_SHR_KIND_R8 ! max value for mvelp - ! dt for averaged cosz calculation - real (SHR_KIND_R8) :: dt_avg=-1 !=============================================================================== CONTAINS !=============================================================================== -subroutine shr_orb_init(dt_in) - - ! Set the dt for the average cosz calculation if it is supplied - real(SHR_KIND_R8), intent(in) :: dt_in - - dt_avg = dt_in - -end subroutine shr_orb_init - - -real(SHR_KIND_R8) pure FUNCTION shr_orb_cosz(jday,lat,lon,declin) +real(SHR_KIND_R8) pure FUNCTION shr_orb_cosz(jday,lat,lon,declin,dt_avg) !---------------------------------------------------------------------------- ! @@ -74,15 +61,24 @@ real(SHR_KIND_R8) pure FUNCTION shr_orb_cosz(jday,lat,lon,declin) real (SHR_KIND_R8),intent(in) :: lat ! Centered latitude (radians) real (SHR_KIND_R8),intent(in) :: lon ! Centered longitude (radians) real (SHR_KIND_R8),intent(in) :: declin ! Solar declination (radians) + real (SHR_KIND_R8),intent(in), optional :: dt_avg ! if present and set non-zero, then use in the + ! average cosz calculation + logical :: use_dt_avg !---------------------------------------------------------------------------- + use_dt_avg = .false. + if (present(dt_avg)) then + if (dt_avg /= 0.0_shr_kind_r8) use_dt_avg = .true. + end if + + ! If dt for the average cosz is specified, then call the shr_orb_avg_cosz - if (dt_avg == -1) then + if (use_dt_avg) then + shr_orb_cosz = shr_orb_avg_cosz(jday, lat, lon, declin, dt_avg) + else shr_orb_cosz = sin(lat)*sin(declin) - & & cos(lat)*cos(declin)*cos(jday*2.0_SHR_KIND_R8*pi + lon) - else - shr_orb_cosz = shr_orb_avg_cosz(jday, lat, lon, declin) end if END FUNCTION shr_orb_cosz @@ -95,7 +91,7 @@ END FUNCTION shr_orb_cosz ! Ref. : Zhou et al., GRL, 2015 !======================================================================= -real (SHR_KIND_R8) pure function shr_orb_avg_cosz(jday, lat, lon, declin) +real (SHR_KIND_R8) pure function shr_orb_avg_cosz(jday, lat, lon, declin, dt_avg) use shr_const_mod, only : pi => shr_const_pi @@ -108,6 +104,7 @@ real (SHR_KIND_R8) pure function shr_orb_avg_cosz(jday, lat, lon, declin) real(SHR_KIND_R8), intent(in) :: lat ! latitude (radian) real(SHR_KIND_R8), intent(in) :: lon ! longitude (radian) real(SHR_KIND_R8), intent(in) :: declin ! solar declination (radian) + real(SHR_KIND_R8), intent(in) :: dt_avg ! dt for averaged cosz calculation !----------------------------------------------------------------------- ! Local Arguments