Skip to content

Commit

Permalink
updating Interp function
Browse files Browse the repository at this point in the history
  • Loading branch information
adrifoster committed Apr 24, 2024
1 parent 673b90a commit f2d58ca
Show file tree
Hide file tree
Showing 2 changed files with 52 additions and 43 deletions.
83 changes: 46 additions & 37 deletions src/biogeophys/DiurnalOzoneType.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,33 +16,33 @@ module DiurnalOzoneType
implicit none
save
private

! !PUBLIC TYPES:
type, public :: diurnal_ozone_anom_type
private
! Private data members
integer :: ntimes ! size of time dimension
real(r8), public, allocatable :: o3_anomaly_grc(:,:) ! o3 anomaly data [grc, ntimes]
real(r8), public, allocatable :: time_arr(:) ! time dimension (units = seconds of day)
real(r8), public, allocatable :: time_arr(:) ! time dimension (units = seconds of day)

contains
! Public routines
procedure, public :: Init
procedure, private :: InitAllocate
procedure, public :: Interp

end type diurnal_ozone_anom_type

character(len=*), parameter, private :: sourcefile = &
__FILE__


contains

! ========================================================================
! Infrastructure routines (initialization, etc.)
! ========================================================================

!-----------------------------------------------------------------------
subroutine Init(this, bounds, n)
!
Expand All @@ -63,8 +63,8 @@ subroutine Init(this, bounds, n)
call this%InitAllocate(bounds)

end subroutine Init


!-----------------------------------------------------------------------
subroutine InitAllocate(this, bounds)
!
Expand Down Expand Up @@ -99,53 +99,62 @@ subroutine Interp(this, forc_o3, forc_o3_down)
!
! !USES:
use clm_time_manager , only : get_curr_time
use clm_varcon , only : secspday
!
! !ARGUMENTS:
class(diurnal_ozone_anom_type), intent(in) :: this
real(r8), intent(in) :: forc_o3( bounds%begg: ) ! ozone partial pressure (mol/mol)
real(r8), intent(out) :: forc_o3_down( bounds%begg: ) ! ozone partial pressure, downscaled (mol/mol)
class(diurnal_ozone_anom_type), intent(in) :: this
real(r8), intent(in) :: forc_o3(bounds%begg:bounds%endg) ! ozone partial pressure (mol/mol)
real(r8), intent(out) :: forc_o3_down(bounds%begg:bounds%endg) ! ozone partial pressure, downscaled (mol/mol)
!
! LOCAL VARIABLES:
integer :: i ! looping index
integer :: yr ! year
integer :: mon ! month
integer :: day ! day of month
integer :: tod ! time of day (seconds past 0Z)
integer :: t ! looping index
integer :: yr ! year
integer :: mon ! month
integer :: day ! day of month
integer :: tod ! time of day (seconds past 0Z)
integer :: begg, endg ! bounds
integer :: t_prev ! previous time index
real(r8) ::
real(r8) :: anomaly_val_start
real(r8) :: anomaly_val_end
real(r8) :: anomaly_scalar
real(r8) :: tdiff_end, tdiff_start
!-----------------------------------------------------------------------

begg = bounds%begg; endg = bounds%endg

! Get current date/time - we really only need seconds
call get_curr_date(yr, mon, day, tod)

! find the time interval we are in
do i = 1, this%ntimes
if (real(tod) <= this%time_arr(i)) then
do t = 1, this%ntimes
if (real(tod) <= this%time_arr(t)) then
exit
end if
end do

! interpolate, checking for edge cases
!TODO: Add bounds here
!TODO: make seconds in day a parameter
if (i == 1) then
! wrap around to back
forc_o3_down(:) = forc_o3(:)*((this%o3_anomaly_grc(:,this%ntimes)* &
(this%time_arr(i) - real(tod)) + &
this%o3_anomaly_grc(:,i)*((86400.0 - this%time_arr(this%ntimes)) + real(tod)))/ &
(this%time_arr(i) + (86400.0 - this%time_arr(this%ntimes))))
else
! interpolate normally
forc_o3_down(:) = forc_o3(:)*((this%o3_anomaly_grc(:,i-1)* &
(this%time_arr(i) - real(tod)) + &
this%o3_anomaly_grc(:,i)*(real(tod) - this%time_arr(i-1)))/ &
(this%time_arr(i) - this%time_arr(i-1)))
end if



if (t == 1) then
t_prev = this%ntimes
tdiff_end = secspday - this%time_arr(t_prev) + real(tod)
tdiff = this%time_arr(t) + secspday - this%time_arr(t_prev)
else
t_prev = t - 1
tdiff_end = real(tod) - this%time_arr(t_prev)
tdiff = this%time_arr(t) - this%time_arr(t_prev)
end if

anomaly_val_start = this%o3_anomaly_grc(begg:endg, t_prev)
anomaly_val_end = this%o3_anomaly_grc(begg:endg, t)
tdiff_start = this%time_arr(t) - real(tod)

! interpolate
anomaly_scalar = (anomaly_val_start*tdiff_start + anomaly_val_end*tdiff_end)/tdiff

forc_o3_down(begg:endg) = forc_o3(begg:endg)*anomaly_scalar

end subroutine Interp

!-----------------------------------------------------------------------

end module DiurnalOzoneType

12 changes: 6 additions & 6 deletions src/biogeophys/OzoneMod.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -154,14 +154,14 @@ subroutine Init(this, bounds, o3_veg_stress_method)
! Local variables
integer :: atm_ozone_frequency_val
!-----------------------------------------------------------------------
! read what atm ozone frequency we have

! read what atm ozone frequency we have
call shr_ozone_coupling_readnl("drv_flds_in", atm_ozone_frequency_val)
this%atm_ozone_freq = atm_ozone_frequency_val

! if we have multi-day average input ozone, we need to convert to sub-daily using
! an input anomaly file
if (this%atm_ozone_freq == atm_ozone_frequency_multiday_average) then
if (this%atm_ozone_freq == atm_ozone_frequency_multiday_average) then
! initialize and read in data for diurnal O3 anomaly stream
call read_O3_stream(this%diurnalOzoneAnomInst, bounds)
end if
Expand Down Expand Up @@ -431,14 +431,14 @@ subroutine CalcOzoneUptake(this, bounds, num_exposedvegp, filter_exposedvegp, &

! Ozone uptake for shaded leaves
call CalcOzoneUptakeOnePoint( &
forc_ozone=forc_o3(g), forc_pbot=forc_pbot(c), forc_th=forc_th(c), &
forc_ozone=forc_o3_down(g), forc_pbot=forc_pbot(c), forc_th=forc_th(c), &
rs=rssha(p), rb=rb(p), ram=ram(p), &
tlai=tlai(p), tlai_old=tlai_old(p), pft_type=patch%itype(p), &
o3uptake=o3uptakesha(p))

! Ozone uptake for sunlit leaves
call CalcOzoneUptakeOnePoint( &
forc_ozone=forc_o3(g), forc_pbot=forc_pbot(c), forc_th=forc_th(c), &
forc_ozone=forc_o3_down(g), forc_pbot=forc_pbot(c), forc_th=forc_th(c), &
rs=rssun(p), rb=rb(p), ram=ram(p), &
tlai=tlai(p), tlai_old=tlai_old(p), pft_type=patch%itype(p), &
o3uptake=o3uptakesun(p))
Expand Down

0 comments on commit f2d58ca

Please sign in to comment.