so4bnd.f90

CCSM Research Tools: Community Atmosphere Model (CAM)

!
! Normal interpolation between consecutive time slices.
!
      do n=1,timsiz-1
         np1 = n + 1
         call bnddyi(date_sulf(n  ), sec_sulf(n  ), cdaysulfm)
         call bnddyi(date_sulf(np1), sec_sulf(np1), cdaysulfp)
         if (.not.sulfcyc) then
            yr = date_sulf(n)/10000
            cdaysulfm = cdaysulfm + yr*365.
            yr = date_sulf(np1)/10000
            cdaysulfp = cdaysulfp + yr*365.
         end if
         if (caldayloc.gt.cdaysulfm .and. caldayloc.le.cdaysulfp) then
            strt4(4) = n
            call wrap_get_vara_realx (ncid_sulf,sulfbio_id,strt4,cnt4,xsulfbioi(1,1,1,nm))
            call wrap_get_vara_realx (ncid_sulf,sulfant_id,strt4,cnt4,xsulfanti(1,1,1,nm))
            strt4(4) = np1
            call wrap_get_vara_realx (ncid_sulf,sulfbio_id,strt4,cnt4,xsulfbioi(1,1,1,np))
            call wrap_get_vara_realx (ncid_sulf,sulfant_id,strt4,cnt4,xsulfanti(1,1,1,np))
            goto 10
         end if
      end do
      write(6,*)'SULFINI: Failed to find dates bracketing ncdate, ncsec=', ncdate, ncsec
      call endrun

10    continue
      write(6,*)'SULFINI: Read sulfate data for dates ',date_sulf(n), &
           sec_sulf(n),' and ',date_sulf(np1),sec_sulf(np1)

#if (defined SPMD )
     call mpibcast( timsiz, 1, mpiint, 0, mpicom )
     call mpibcast( date_sulf, 1000, mpiint, 0, mpicom )
     call mpibcast( sec_sulf, 1000, mpiint, 0, mpicom )
     call mpibcast( cdaysulfm, 1, mpir8, 0, mpicom )
     call mpibcast( cdaysulfp, 1, mpir8, 0, mpicom )
     call mpibcast( np1, 1, mpiint, 0, mpicom )
  else
     call mpibcast( timsiz, 1, mpiint, 0, mpicom )
     call mpibcast( date_sulf, 1000, mpiint, 0, mpicom )
     call mpibcast( sec_sulf, 1000, mpiint, 0, mpicom )
     call mpibcast( cdaysulfm, 1, mpir8, 0, mpicom )
     call mpibcast( cdaysulfp, 1, mpir8, 0, mpicom )
     call mpibcast( np1, 1, mpiint, 0, mpicom )
#endif
   endif

   call scatter_field_to_chunk(1,pver,2,plon,xsulfbioi,sulfbioi)
   call scatter_field_to_chunk(1,pver,2,plon,xsulfanti,sulfanti)

   return
end subroutine sulfini

!###############################################################################

subroutine sulfint
!----------------------------------------------------------------------- 
! 
! Purpose: Interpolate sulfate mixing ratios to current time, reading in new monthly
!          data if necessary, and spatially interpolating it.
! 
!-----------------------------------------------------------------------
   use time_manager, only: get_curr_date, get_perp_date, get_curr_calday, &
                           is_perpetual
#if ( defined SPMD )
   use mpishorthand
#endif

   implicit none

#include <comctl.h>
#include <comlun.h>
!-----------------------------------------------------------------------
   include 'netcdf.inc'
!-----------------------------------------------------------------------
!
! Local workspace
!
   integer cnt4(4)            ! array of counts for each dimension
   integer strt4(4)           ! array of starting indices
   integer i, k, lchnk        ! column, level, chunk indices
   integer ncol               ! number of columns in current chunk
   integer ntmp               ! temporary
   integer  :: yr, mon, day   ! components of a date
   integer  :: ncdate         ! current date in integer format [yyyymmdd]
   integer  :: ncsec          ! current time of day [seconds]
   real(r8) :: calday         ! current calendar day
   real(r8) fact1, fact2      ! time interpolation factors
   real(r8) caldayloc         ! calendar day (includes yr if no cycling)
   real(r8) deltat            ! time (days) between interpolating sulfate data
   real(r8) xsulfbioi(plon,pver,plat)  ! input sulfate bio mixing ratios
   real(r8) xsulfanti(plon,pver,plat)  ! input sulfate ant mixing ratios

!-----------------------------------------------------------------------
!
! Use year information only if a multiyear dataset
!
   calday = get_curr_calday()
   if ( is_perpetual() ) then
      call get_perp_date(yr, mon, day, ncsec)
   else
      call get_curr_date(yr, mon, day, ncsec)
   end if
   ncdate = yr*10000 + mon*100 + day
   if (sulfcyc) then
      caldayloc = calday
   else
      yr = ncdate/10000
      caldayloc = calday + yr*365.
   end if
!
! Initialize hyperslab corners
!
   if (masterproc) then

      strt4(1) = 1
      strt4(2) = 1
      strt4(3) = 1
      cnt4(1)  = lonsiz
      cnt4(2)  = levsiz
      cnt4(3)  = latsiz
      cnt4(4)  = 1

   endif
!
! If model time is past current forward sulfate timeslice, read in the next
! timeslice for time interpolation.  Messy logic is for sulfcyc = .true. 
! interpolation between December and January (np1.eq.1).  Note that 
! np1 is never 1 when sulfcyc is .false.
!
   if (caldayloc .gt. cdaysulfp .and. .not. (np1.eq.1 .and. caldayloc.gt.cdaysulfm)) then
      if (sulfcyc) then
         np1 = mod(np1,12) + 1
      else
         np1 = np1 + 1
      end if
      if (np1.gt.timsiz) then
         if (masterproc) then
            write(6,*)'SULFINT: Attempt to read past end of dataset'
         endif
         call endrun
      end if
      cdaysulfm = cdaysulfp
      call bnddyi(date_sulf(np1), sec_sulf(np1), cdaysulfp)
      if (.not.sulfcyc) then
         yr = date_sulf(np1)/10000
         cdaysulfp = cdaysulfp + yr*365.
      end if
      if (np1.eq.1 .or. caldayloc.le.cdaysulfp) then
         ntmp = nm
         nm = np
         np = ntmp
         if (masterproc) then
            strt4(4) = np1
            call wrap_get_vara_realx (ncid_sulf,sulfbio_id,strt4,cnt4,xsulfbioi(1,1,1))
            call wrap_get_vara_realx (ncid_sulf,sulfant_id,strt4,cnt4,xsulfanti(1,1,1))
            write(6,*)'SULFINT: Read sulfate for date (yyyymmdd) ',date_sulf(np1),' sec ', &
                      sec_sulf(np1)
         endif
         call scatter_field_to_chunk(1,pver,1,plon,xsulfbioi,sulfbioi(1,1,begchunk,np))
         call scatter_field_to_chunk(1,pver,1,plon,xsulfanti,sulfanti(1,1,begchunk,np))
      else
         if (masterproc) then
            write(6,*)'SULFINT: Input sulfate for date',date_sulf(np1), &
                      ' sec ',sec_sulf(np1), 'does not exceed model date', &
                      ncdate,' sec ',ncsec,' Stopping.'
         endif
         call endrun
      end if
   end if
!
! Determine time interpolation factor.  Account for December-January 
! interpolation if cycling sulfate dataset.  Again note that np1 is never 1 
! when sulfcyc is false
!
   if (np1.eq.1) then                    ! Dec-Jan interpolation
      deltat = cdaysulfp + 365. - cdaysulfm
      if (caldayloc.gt.cdaysulfp) then   ! We're in December
         fact1 = (cdaysulfp + 365. - caldayloc)/deltat
         fact2 = (caldayloc - cdaysulfm)/deltat
      else                               ! We're in January
         fact1 = (cdaysulfp - caldayloc)/deltat
         fact2 = (caldayloc + 365. - cdaysulfm)/deltat
      end if
   else
      deltat = cdaysulfp - cdaysulfm
      fact1 = (cdaysulfp - caldayloc)/deltat
      fact2 = (caldayloc - cdaysulfm)/deltat
   end if
!
! Check sanity of time interpolation calculation to within 32-bit roundoff
!
   if (abs(fact1+fact2-1.).gt.1.e-6 .or.  &
        fact1.gt.1.000001 .or. fact1.lt.-1.e-6 .or. &
        fact2.gt.1.000001 .or. fact2.lt.-1.e-6) then
      if (masterproc) then
         write(6,*)'SULFINT: Bad fact1 and/or fact2=',fact1,fact2
      endif
      call endrun
   end if
!
! Time interpolation.
!
   do lchnk=begchunk,endchunk
      ncol = get_ncols_p(lchnk)
      do k=1,pver
         do i=1,ncol
            sulfbio(i,k,lchnk) = sulfbioi(i,k,lchnk,nm)*fact1 +  &
                               sulfbioi(i,k,lchnk,np)*fact2
            sulfant(i,k,lchnk) = sulfanti(i,k,lchnk,nm)*fact1 +  &
                               sulfanti(i,k,lchnk,np)*fact2
         end do
      end do
   end do

   return
end subroutine sulfint

!###############################################################################

subroutine getso4bnd( lchnk, ncol, bio, anth )

!----------------------------------------------------------------------- 
! 
! Purpose: Return slice of time interpolated sulfate aerosol.
! 
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
   implicit none
!-----------------------------------------------------------------------
!
! Arguments
!
  integer , intent(in)   :: lchnk            ! chunk identifier
  integer , intent(in)   :: ncol             ! number of atmospheric columns

   real(r8), intent(out) :: bio(pcols,pver)  ! biogenic sulfate
   real(r8), intent(out) :: anth(pcols,pver) ! anthropogenic sulfate
!
! Local variables.
!
   integer :: i, k                 ! longitude, level indices
!-----------------------------------------------------------------------
   do k = 1, pver
      do i = 1, ncol
         bio(i,k)  = sulfbio(i,k,lchnk)
         anth(i,k) = sulfant(i,k,lchnk)
      end do
   end do
   return
end subroutine getso4bnd

!###############################################################################

subroutine setso4ramp( x )

!----------------------------------------------------------------------- 
! 
! Purpose: Set so4 ramp value.
! 
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
   implicit none
!-----------------------------------------------------------------------
!
! Arguments
!
   real(r8),intent(in) :: x  ! sulfate scale factor computed in ramp subroutine
!-----------------------------------------------------------------------
   sulfscalef = x
   return
end subroutine setso4ramp

!###############################################################################

real*8 function so4ramp()
!----------------------------------------------------------------------- 
! 
! Purpose: Return so4 ramp value.
! 
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
   implicit none
!-----------------------------------------------------------------------
   so4ramp = sulfscalef

   return
end function so4ramp

end module so4bnd