#include <misc.h>#include <params.h>subroutine oznint!----------------------------------------------------------------------- ! ! Purpose: Interpolate ozone mixing ratios to current time, reading in new monthly!          data if necessary, and spatially interpolating it.! ! Method: Find next month of ozone data to interpolate.  Linearly interpolate !         vertically and horizontally! ! Author: CCM Core Group! !-----------------------------------------------------------------------   use precision   use pmgrid   use comozp   use pspect   use rgrid   use commap   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>!! Local workspace!   integer cnt4(4)                ! array of counts for each dimension   integer strt4(4)               ! array of starting indices   integer i, k, lat              ! longitude, level, latitude indices   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) fact1, fact2          ! time interpolation factors   real(r8) :: calday             ! current calendar day   real(r8) caldayloc             ! calendar day (includes yr if no cycling)   real(r8) tmpozmix(levsiz,plat) ! temporary ozmix array   real(r8) deltat                ! time (days) between interpolating ozone data   real(r8), allocatable :: oznbdyp(:,:,:) ! ozone data from next time sample!! SPMD: Master does all the work.  Sends needed info to slaves!   if (masterproc) then!! 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 (ozncyc) then         caldayloc = calday      else         caldayloc = calday + yr*365.      end if      strt4(1) = 1      strt4(2) = 1      strt4(3) = 1      cnt4(1)  = lonsiz      cnt4(2)  = levsiz      cnt4(3)  = latsiz      cnt4(4)  = 1!! If model time is past current forward ozone timeslice, read in the next! timeslice for time interpolation.  Messy logic is for ozncyc = .true. ! interpolation between December and January (np1 == 1).  Note that ! np1 is never 1 when ozncyc is .false.!      if (caldayloc > cdayozp .and. .not. (np1 == 1 .and. caldayloc > cdayozm)) then!! Allocate memory for dynamic local workspace!         allocate (oznbdyp(lonsiz,levsiz,latsiz))         if (ozncyc) then            np1 = mod(np1,12) + 1         else            np1 = np1 + 1         end if         if (np1 > timesiz) then            write(6,*)'OZNINT: Attempt to read past end of O3 dataset'            call endrun         end if         cdayozm = cdayozp         call bnddyi(date_oz(np1), sec_oz(np1), cdayozp)         if (.not.ozncyc) then            yr = date_oz(np1)/10000            cdayozp = cdayozp + yr*365.         end if         if (np1 == 1 .or. caldayloc <= cdayozp) then            ntmp = nm            nm = np            np = ntmp            strt4(4) = np1            call wrap_get_vara_realx (ncid_oz,oznid,strt4,cnt4,oznbdyp)            write(6,*)'OZNINT: Read ozone for date (yyyymmdd) ', date_oz(np1),' sec ',sec_oz(np1)!! Spatial interpolation.  If ozone dataset is only 2-d (i.e. lonsiz = 1) and ! thus only latitude interpolation is necessary, expand to 3-d after ! interpolation.!            if (lonsiz == 1) then               call lininterp (oznbdyp, ozlat, levsiz, latsiz, tmpozmix, &                               latdeg, plat)               do lat=1,plat                  do k=1,levsiz                     do i=1,nlon(lat)                        ozmixm(i,k,lat,np) = tmpozmix(k,lat)                     end do                  end do               end do            else               call bilin (oznbdyp ,ozlon   ,ozlat   ,lonsiz  ,lonsiz  , &                           levsiz  ,levsiz  ,latsiz  ,ozmixm(1,1,1,np) , londeg  , &                           latdeg  ,plond   ,nlon    ,levsiz  ,plat    )            end if         else            write(6,*)'OZNINT: Input ozone for date',date_oz(np1),' sec ',sec_oz(np1), &                 'does not exceed model date',ncdate,' sec ',ncsec,' Stopping.'            call endrun         end if!! Deallocate dynamic memory for local workspace.!         deallocate (oznbdyp)      end if!! Determine time interpolation factor.  Account for December-January ! interpolation if cycling ozone dataset.  Again note that np1 is never 1 ! when ozncyc is false!      if (np1 == 1) then                    ! Dec-Jan interpolation         deltat = cdayozp + 365. - cdayozm         if (caldayloc > cdayozp) then      ! We're in December            fact1 = (cdayozp + 365. - caldayloc)/deltat            fact2 = (caldayloc - cdayozm)/deltat         else                                ! We're in January            fact1 = (cdayozp - caldayloc)/deltat            fact2 = (caldayloc + 365. - cdayozm)/deltat         end if      else         deltat = cdayozp - cdayozm         fact1 = (cdayozp - caldayloc)/deltat         fact2 = (caldayloc - cdayozm)/deltat      end if!! Check sanity of time interpolation calculation to within 32-bit roundoff!      if (abs(fact1+fact2-1.) > 1.e-6 .or. fact1 > 1.000001 .or. fact1 < -1.e-6 .or. &          fact2 > 1.000001 .or. fact2 < -1.e-6) then         write(6,*)'OZNINT: Bad fact1 and/or fact2=',fact1,fact2         call endrun      end if!! Time interpolation.!      do lat=1,plat         do k=1,levsiz            do i=1,nlon(lat)               ozmix(i,k,lat) = ozmixm(i,k,lat,nm)*fact1 + ozmixm(i,k,lat,np)*fact2            end do         end do      end do   end if#if ( defined SPMD )   call mpibcast(ozmix, plond*levsiz*plat, mpir8, 0, mpicom)#endif   returnend subroutine oznint