! source file: /sfs/fs1/work-geomar6/smomw113/UVic_ESCM/2.9_dk_0814/updates/07_aouassim/source/mom/respidata.F
      subroutine respidata(respireturn)

!=======================================================================
!     routine to read and interpolate three dimensional respiration data
!     The original version (March 2017) actually read in anomalies of o_detrremi.
!     In Feb 2018, the code has been modified in order to
!
!     wkoeve@geomar.de, 08.03.17, 23.02.2018
!
!=======================================================================

      implicit none

      character(120) :: fname, name, new_file_name, text, taxname
     &,    variname

      integer iou, n, respidata_n, ln,ib(10), ic(10)
      integer c,i,j,k
      integer :: AllocateStatus, DeAllocateStatus

      logical inqvardef, exists
      logical show_this, nfound

      real data_time, wt1, wt3, c1c3

      real, allocatable :: time(:)
!      real, allocatable :: tim(:,:,:,:)
      real tim(3)
      real, allocatable :: dat(:,:,:,:)
      real, allocatable :: datarespi(:,:,:,:)

      save dat, tim, ln, respidata_n, time, iou
      save datarespi

!     All variable values stored via save are available each time sstdata()
!     is being entered, kind of internal global variables.

      include "size.h"
      include "param.h"
      include "pconst.h"
      include "stdunits.h"
      include "calendar.h"
      include "cembm.h"
      include "switch.h"
      include "tmngr.h"
      include "levind.h"

      real respireturn(imt,jmt,km)

! Objective of code: Read in detrremi data which (currently) reflect the
! difference between two model runs, e.g. a fully coupled and a 'w const anthro' run.
! See further handling in tracer.F
!

! Notes: O_detrremi is stored in tavg with the factor c1e3. Hence the averaged
!        values of remi are devided by c1e3 when saving. Storage units are
!        mol N m-3 s-1 (equivalent to umol N cm-3 s-1). Units at runtime are
!        [nmol cm-3 s-1], see common/npzd.h, rremi = rate of remineralization.
!        Hence we mulitply the return data with c1c3

        name = "O_detrremi.nc"
        taxname='MYTIME'
        variname='O_REMIFORC'
        fname = new_file_name (name)

      c1c3=1000.

      if (.not. allocated (time)) then
!       variable time has not yet been alocated, hence this section of the code is only entered
!       once during the model run
        inquire (file=trim(fname), exist=exists)
        if (.not. exists) then
          print*, "==> Fatal Error: ", trim(fname), " does not exist."
          print*, "stopping UVic in respidata.F"
          stop

        else
          call openfile (fname, iou)
          call getdimlen (taxname, iou, ln)
          write(20,*) 'respidata: dim of Mytime: ', ln
          allocate ( time(ln) )

!         note: the file has been written with ferret and has an offset of 1700 yrs
          call getvara (taxname, iou, ln, 1, ln, time, c1, c0)
          text = 'years'
          call getatttext (iou, taxname, 'units', text)
          if (trim(text) .eq. "days since 1-1-1")
     &      time(:) = time(:)/yrlen - 1.
          if (trim(text) .eq. "days since 0-1-1")
     &       time(:) = time(:)/yrlen
          if (trim(text) .eq. "years since 1-1-1")
     &      time(:) = time(:) - 1.
          do c=1,ln
             time(c)=time(c)+1700
          enddo

          allocate ( dat(imtm2,jmtm2,km,3), STAT = AllocateStatus )
          IF (AllocateStatus /= 0) STOP "*** Not enough memory ***"
          allocate (datarespi(imtm2,jmtm2,km,2),STAT = AllocateStatus )
          IF (AllocateStatus /= 0) STOP "*** Not enough memory ***"

          datarespi(:,:,:,:) = c0
          dat(:,:,:,:) = c0
          tim(:) = time(1)
          respidata_n=2.
!         make sure there are some data in datarespi
          write(20,*) 'load datarespi for the first time'
          call respidata_load(iou,variname,show_this,2,datarespi)

          write(20,*)
     &       'Note: O_detrremi time-data have been loaded into memory.'

          show_this=.false.
          if (show_this) then
            do c=1,ln
              write(20,*) c, time(c)
            enddo
          else
            write(20,*) 'initial probing of time(c) is turned off'
          endif
          show_this=.false.

        endif
      endif

      print*, "==> Fatal error: O_respi_data_transient has to be ON"
      stop

      data_time = year0 + accel_yr0 + (relyr - accel_yr0)*accel
      tim(2) = min(time(ln), max(time(1), data_time))

      if (tim(2) .le. time(1)) then
          call openfile (fname, iou)
          write(20,*) 'time prior to respi-forc data, use first value'
          call respidata_load(iou,variname,show_this
     &,            2,datarespi)
          do i=1,imtm2
            do j=1,jmtm2
              do k=1,km
                  dat(i,j,k,2) =datarespi(i,j,k,1)
              enddo
            enddo
          enddo

      elseif (tim(2) .ge. time(ln)) then
             call openfile (fname, iou)
             write(20,*) 'time post to respi-forc data, use last value'
             call respidata_load(iou,variname,show_this
     &,            ln,datarespi)
             do i=1,imtm2
               do j=1,jmtm2
                 do k=1,km
                  dat(i,j,k,2) =datarespi(i,j,k,2)
                 enddo
               enddo
             enddo

      else
          if (tim(2) .gt. tim(3)) then
            do n=respidata_n-1,ln
              nfound=.false.
              if (time(n-1).le.tim(2)
     &           .and. time(n).ge.tim(2)) then
                tim(1) = time(n-1)
                tim(3) = time(n)
                if (n .ne. respidata_n) then
!                  load the two data slices (n-1) and (n)
                   respidata_n=n
                   call openfile (fname, iou)
                   write(20,*) 'load new datarespi, n', respidata_n
                   write(20,*) 'variname: ', trim(variname)
                   write(20,*) 'iou:      ', iou
                   call respidata_load(iou,variname,show_this
     &,               respidata_n,datarespi)
                else
!                  use the field datarespi which is already in memory
                endif
!               do the interpolation
                do i=1,imtm2
                  do j=1,jmtm2
                    do k=1,km
                       dat(i,j,k,1) = datarespi(i,j,k,1)
                       dat(i,j,k,3) = datarespi(i,j,k,2)
                       dat(i,j,k,2) = dat(i,j,k,1)*wt1
     &                    + dat(i,j,k,3)*wt3
                    enddo
                  enddo
                enddo
                nfound=.true.
!             next line, of: (time(n-1).le.tim(2) .and. time(n).ge.tim(2))
              endif
              if (nfound) exit
!           next line, of:  do n=respidata_n-1,ln
            enddo
          else
             write(20,*) ' tim(2) lt tim(3)'
!         next line, of: (tim(2) .gt. tim(3))
          endif

          wt1 = 1.
          if (tim(3) .ne. tim(1))
     &         wt1 = (tim(3)-tim(2))/(tim(3)-tim(1))
          wt1 = max(0., min(1., wt1))
          wt3 = 1. - wt1
          do i=1,imtm2
             do j=1,jmtm2
                do k=1,km
                  dat(i,j,k,2) = dat(i,j,k,1)*wt1
     &                    + dat(i,j,k,3)*wt3
                enddo
             enddo
          enddo
          show_this=.true.
          if (show_this) then
             write(20,*) 'data_t: ', data_time
!             write(20,*) 'time(1):', time(1)
             write(20,*) 'n:      ', respidata_n
!             write(20,*) 'km:     ', km
!             write(20,*) 'tim(1): ', tim(1)
!             write(20,*) 'tim(2): ', tim(2)
!             write(20,*) 'tim(3): ', tim(3)
!             write(20,*) 'dat(1): ', dat(50,50,1,1)
             write(20,*) 'dat(2): ', dat(50,50,1,2)
!             write(20,*) 'dat(2): ', dat(50,50,2,2)
!             write(20,*) 'dat(2): ', dat(50,50,3,2)
!             write(20,*) 'dat(2): ', dat(50,50,4,2)
!             write(20,*) 'dat(2): ', dat(50,50,5,2)
!             write(20,*) 'dat(2): ', dat(50,50,6,2)
!             write(20,*) 'dat(2): ', dat(50,50,7,2)
!             write(20,*) 'dat(3): ', dat(50,50,1,3)
          endif
          show_this=.false.
!     next line, of: if, elseif, else big loop
      endif

!     Prepare return variable, incl conversion to run time units
      do i=2,imtm1
         do j=2,jmtm1
            do k=1,km

                 respireturn(i,j,k) = dat(i-1,j-1,k,2)*c1c3

            enddo
         enddo
      enddo

! set zonal boundary conditions of respireturn
      do k=1,km
        call setbcx (respireturn(1,1,k), imt, jmt)
      enddo

      return
      end

      subroutine respidata_load(iou,variname,show_this,n,datarespi)

      implicit none

      character(120) :: fname, name, new_file_name, text
     &,     variname

      integer iou, n, ln, nfound,ib(10), ic(10)
      integer c,i,j,k , kmx
      integer :: AllocateStatus, DeAllocateStatus

      logical inqvardef, exists
      logical show_this

      include "size.h"
      include "param.h"
      include "pconst.h"
      include "stdunits.h"
      include "calendar.h"
      include "levind.h"

      real tmpij(1:imtm2,1:jmtm2)
      real datarespi(1:imtm2,1:jmtm2,1:km,1:2)

      datarespi(1:imtm2,1:jmtm2,1:km,1:2) = c0
      tmpij(1:imtm2,1:jmtm2) = c0

      exists = inqvardef(trim(variname), iou)
      if (.not. exists) then
         write(20,*) '==>  Fatal Error: ',trim(variname)
     &,      ' data does not exist.'
         stop
      else

!         write(20,*) 'respidata_load value of n:', n
         do k=1,km
!          write(20,*) 'n-1, k-values:', k
           ib(:) = 1
           ic(:) = 1
           ib(1)=1
           ib(2)=1
           ib(3)=k
           ib(4)=n-1
           ic(1)=imtm2
           ic(2)=jmtm2
           ic(3)=1
           ic(4)=1
           call getvara (variname, iou, imtm2*jmtm2, ib, ic
     &,       tmpij, c1, c0)

           do i=1,imtm2
             do j=1,jmtm2
                  datarespi(i,j,k,1) =
     &            tmpij(i,j)
             enddo
           enddo
         enddo

         do k=1,km
!           write(20,*) 'n, k-values:', k
           ib(:) = 1
           ic(:) = 1
           ib(1)=1
           ib(2)=1
           ib(3)=k
           ib(4)=n
           ic(1)=imtm2
           ic(2)=jmtm2
           ic(3)=1
           ic(4)=1
           call getvara (variname, iou, imtm2*jmtm2, ib, ic
     &,       tmpij, c1, c0)

           do i=1,imtm2
             do j=1,jmtm2
                  datarespi(i,j,k,2) =
     &            tmpij(i,j)
             enddo
           enddo
         enddo

            write(20,*)
     &        'respidata_load: data have been loaded into memory.'
      endif

      return
      end