! source file: /Users/nmengis/Documents/UVic_ESCM/2.10/updates/02_CE_permafrost_merge_CMIP6forcing/source/mtlm/mtlm_tavg.F

      subroutine mtlm_tavg_def (fname, imt, jmt, NPFT, NTYPE, NGND
     &,                         xt, yt, calendar, expnam, runstamp)

!=======================================================================
!     definition routine for land time averages

!   inputs:
!     fname        = file name
!     imt, jmt ... = global array dimensions
!     xt, yt ...   = global axes
!     calendar     = calendar
!     expnam       = experiment name
!     runstamp     = run stamp
!=======================================================================

      implicit none

      character(*) :: fname, calendar, expnam, runstamp

      integer iou, j, imt, jmt, igs, ige, ig, jgs, jge, jg
      integer id_time, id_xt, id_yt, id_pft, id_type
      integer id_xt_e, id_yt_e, id_pft_e, id_type_e
      integer it(10)

      integer id_wetm, id_wetm_e
      integer id_frozenm, id_frozenm_e
      integer id_gnd, id_gnd_e, NGND

      integer NPFT, NTYPE

      real xt(imt), yt(jmt)

      real c0, c1, c1e6, c1e20

      c0 = 0.
      c1 = 1.
      c1e6 = 1.e6
      c1e20 = 1.e20

!-----------------------------------------------------------------------
!     open file
!-----------------------------------------------------------------------
      call openfile (fname, iou)

!-----------------------------------------------------------------------
!     set global write domain size (may be less than global domain)
!-----------------------------------------------------------------------
      igs = 1
      ige = imt
      if (xt(1) + 360. lt. xt(imt)) then
!       assume cyclic boundary
        igs = 2
        ige = imt-1
      endif
      ig  = ige-igs+1
      jgs = 1
      jge = jmt
      do j=2,jmt
        if (yt(j-1) .lt. -90. .and. yt(j) .gt. -90.) jgs = j
        if (yt(j-1) .lt.  90. .and. yt(j) .gt. 90.) jge = j-1
      enddo
      jg  = jge-jgs+1

!-----------------------------------------------------------------------
!     start definitions
!-----------------------------------------------------------------------
      call redef (iou)

!-----------------------------------------------------------------------
!     write global attributes
!-----------------------------------------------------------------------
      call putatttext (iou, 'global', 'Conventions', 'CF-1.0')
      call putatttext (iou, 'global', 'experiment_name', expnam)
      call putatttext (iou, 'global', 'run_stamp', runstamp)

!-----------------------------------------------------------------------
!     define dimensions
!-----------------------------------------------------------------------
      call defdim ('time', iou, 0, id_time)
      call defdim ('longitude', iou, ig, id_xt)
      call defdim ('latitude', iou, jg, id_yt)
      call defdim ('pft', iou, NPFT, id_pft)
      call defdim ('type', iou, NTYPE, id_type)
      call defdim ('longitude_edges', iou, ig+1, id_xt_e)
      call defdim ('latitude_edges', iou, jg+1, id_yt_e)
      call defdim ('pft_edges', iou, NPFT+1, id_pft_e)
      call defdim ('type_edges', iou, NTYPE+1, id_type_e)

      call defdim ('wetm',iou, 3, id_wetm)
      call defdim ('wetm_edges', iou, 4, id_wetm_e)
      call defdim ('frozenm', iou, 7, id_frozenm)
      call defdim ('frozenm_edges', iou, 8, id_frozenm_e)
      call defdim ('gnd', iou, NGND, id_gnd)
      call defdim ('gnd_edges', iou, NGND+1, id_gnd_e)

!-----------------------------------------------------------------------
!     define 1d data (t)
!-----------------------------------------------------------------------
      it(1) = id_time
      call defvar ('time', iou, 1, it, c0, c0, 'T', 'D'

     &, 'time', 'time', 'years since 0-1-1')

      call putatttext (iou, 'time', 'calendar', calendar)
      call defvar ('T_avgper', iou, 1, it, c0, c0, ' ', 'F'
     &,   'averaging period', ' ','day')

!-----------------------------------------------------------------------
!     define 1d data (x, y or n)
!-----------------------------------------------------------------------
      it(1) = id_xt
      call defvar ('longitude', iou, 1, it, c0, c0, 'X', 'D'
     &, 'longitude', 'longitude', 'degrees_east')
      call defvar ('G_dxt', iou, 1, it, c0, c0, ' ', 'D'
     &, 'width t grid', ' ', 'm')
      it(1) = id_yt
      call defvar ('latitude', iou, 1, it, c0, c0, 'Y', 'D'
     &, 'latitude', 'latitude', 'degrees_north')
      call defvar ('G_dyt', iou, 1, it, c0, c0, ' ', 'D'
     &, 'height t grid', ' ', 'm')
      it(1) = id_pft
      call defvar ('pft', iou, 1, it, c0, c0, 'Z', 'F'
     &,   'plant functional type', ' ', ' ')
      it(1) = id_type
      call defvar ('type', iou, 1, it, c0, c0, 'Z', 'F'
     &,   'land type', ' ', ' ')
      it(1) = id_xt_e
      call defvar ('longitude_edges', iou, 1, it, c0, c0, 'X', 'D'
     &,   'longitude edges', 'longitude', 'degrees_east')
      it(1) = id_yt_e
      call defvar ('latitude_edges', iou, 1, it, c0, c0, 'Y', 'D'
     &,   'latitude edges', 'latitude', 'degrees_north')
      it(1) = id_pft_e
      call defvar ('pft_edges', iou, 1, it, c0, c0, 'Z', 'F'
     &,   'plant functional type', ' ', ' ')
      it(1) = id_type_e
      call defvar ('type_edges', iou, 1, it, c0, c0, 'Z'
     &,   'F', 'land type', ' ', ' ')

      it(1) = id_wetm
      call defvar ('wetm', iou, 1, it, c0, c0, 'Z', 'F'
     &,   'wetland fraction', ' ', ' ')
      it(1) = id_frozenm
      call defvar ('frozenm', iou, 1, it, c0, c0, 'Z', 'F'
     &,   'frozen ground map', ' ', ' ')
      it(1) = id_gnd
      call defvar ('gnd', iou, 1, it, c0, c0, ' ', 'F'
     &,   'gnd hydrological level', ' ', ' ')
      it(1) = id_wetm_e
      call defvar ('wetm_edges', iou, 1, it, c0, c0, 'Z', 'F'
     &,   'wetland fraction', ' ', ' ')
      it(1) = id_frozenm_e
      call defvar ('frozenm_edges', iou, 1, it, c0, c0, 'Z', 'F'
     &,   'frozen ground map', ' ', ' ')
      it(1) = id_gnd_e
      call defvar ('gnd_edges', iou, 1, it, c0, c0, ' ', 'F'
     &,   'gnd hydrological level', ' ', ' ')

!-----------------------------------------------------------------------
!     define 2d data (x,y)
!-----------------------------------------------------------------------
      it(1) = id_xt
      it(2) = id_yt
      call defvar ('G_latT', iou, 2, it, -c1e6, c1e6, ' ', 'F'
     &,   'tracer grid latitude', 'latitude', 'degrees_north')
       call defvar ('G_lonT', iou, 2, it, -c1e6, c1e6, ' ', 'F'
     &,   'tracer grid longitude', 'longitude', 'degrees_east')
      call defvar ('G_areaT', iou, 2, it, -c1e6, c1e6, ' ', 'F'
     &, 'tracer grid area', ' ', 'm2')

!-----------------------------------------------------------------------
!     define 3d data (x,y,t)
!-----------------------------------------------------------------------
      it(1) = id_xt
      it(2) = id_yt
      it(3) = id_time

      call defvar ('L_soilcarb', iou , 3, it, -c1e20, c1e20, ' ', 'F'
     &,   'soil carbon', 'soil_carbon_content', 'kg m-2 ')
      call defvar ('L_soilresp', iou , 3, it, -c1e20, c1e20, ' ', 'F'
     &,   'soil respiration', 'soil_respiration_carbon_flux'
     &,   'kg m-2 s-1')
      call defvar ('L_veglit', iou , 3, it, -c1e20, c1e20, ' ', 'F'
     &,   'total leaf litter', 'leaf_litter_carbon_flux'
     &,   'kg m-2 s-1')
      call defvar ('L_vegburn', iou , 3, it, -c1e20, c1e20, ' ', 'F'
     &,   'total vegetation burning', 'burning_carbon_flux'
     &,   'kg m-2 s-1')

      call defvar ('L_schDmin', iou , 3, it, -c1e20, c1e20, ' ', 'F'
     &,   'Depth of active C', 'Depth of active C', 'layer')
      call defvar ('L_permacarb', iou , 3, it, -c1e20, c1e20, ' ', 'F'
     &,   'permafrost carbon', 'released PF carbon', 'kg m-2 ')
      call defvar ('L_squedcarb', iou , 3, it, -c1e20, c1e20, ' ', 'F'
     &,   'frozen carbon', 'sequestered (frozen) carbon', 'kg m-2 ')
      call defvar ('L_perma', iou , 3, it, -c1e20, c1e20, ' ', 'F'
     &,   'depth to PF', 'depth to PF', 'm')
      call defvar ('L_activelt', iou , 3, it, -c1e20, c1e20, ' ', 'F'
     &,   'active layer thickness', 'active layer thickness', 'm')
      call defvar ('L_talik', iou , 3, it, -c1e20, c1e20, ' ', 'F'
     &,   'talik thickness', 'talik thickness', 'm')
      call defvar ('L_dwet', iou , 3, it, -c1e20, c1e20, ' ', 'F'
     &,   'days wet', 'average number of days wet', 'm')
      call defvar ('L_PFthick', iou , 3, it, -c1e20, c1e20, ' ', 'F'
     &,   'PF thickness', 'permafrost thickness', 'm')

!-----------------------------------------------------------------------
!     define 4d data (x,y,pft,t)
!-----------------------------------------------------------------------
      it(1) = id_xt
      it(2) = id_yt
      it(3) = id_pft
      it(4) = id_time
      call defvar ('L_veggpp', iou , 4, it, -c1e20, c1e20, ' ', 'F'
     &,   'gross primary productivity'
     &,   'gross_primary_productivity_of_carbon', 'kg m-2 s-1')
      call defvar ('L_vegnpp', iou , 4, it, -c1e20, c1e20, ' ', 'F'
     &,   'net primary productivity'
     &,   'net_primary_productivity_of_carbon', 'kg m-2 s-1')
      call defvar ('L_veghgt', iou , 4, it, -c1e20, c1e20, ' ', 'F'
     &,   'canopy height', 'canopy_height', 'm')
      call defvar ('L_veglai', iou , 4, it, -c1e20, c1e20, ' ', 'F'
     &,   'leaf area index ', 'leaf_area_index', '1')
      call defvar ('L_vegcarb', iou , 4, it, -c1e20, c1e20, ' ', 'F'
     &,   'vegetation carbon', 'vegetation_carbon_content', 'kg m-2 ')

       call defvar ('L_vegfsmc', iou, 4, it, -c1e20, c1e20, ' ', 'F'
     &, 'fsmc', 'vegetation_fsmc', '1')

!-----------------------------------------------------------------------
!     define 4d data (x,y,wetland,t)
!-----------------------------------------------------------------------
      it(1) = id_xt
      it(2) = id_yt
      it(3) = id_wetm
      it(4) = id_time
      call defvar ('L_wetmap', iou , 4, it, -c1e20, c1e20, ' ', 'F'
     &,   'fraction of gridcell containing wetlands'
     &,   'fraction of gridcell containing wetlands', '')

!-----------------------------------------------------------------------
!     define 4d data (x,y,frozenmap,t)
!-----------------------------------------------------------------------
      it(1) = id_xt
      it(2) = id_yt
      it(3) = id_frozenm
      it(4) = id_time
      call defvar ('L_frozenmap', iou , 4, it, -c1e20, c1e20, ' ', 'F'
     &,   'map of frozen ground types'
     &,   'map of frozen ground types', '')

!---------------------------------------------------------------------
!     define 4d data (x,y,ngnd,t)
!---------------------------------------------------------------------
      it(1) = id_xt
      it(2) = id_yt
      it(3) = id_gnd
      it(4) = id_time
      call defvar('L_soilcarbD', iou, 4, it, -c1e20, c1e20, ' ', 'F'
     &,   'soil carbon w depth'
     &,  'soil_carbon_density', 'kg m-3')
      call defvar('L_soilrespD', iou, 4, it, -c1e20, c1e20, ' ', 'F'
     &,   'soil respiration w depth'
     &,  'soil_respiraction', 'kg m-3 s-1')

      call defvar('L_PermaCarbD', iou, 4, it, -c1e20, c1e20, ' ', 'F'
     &,   'perma carbon w depth'
     &,  'perma_carbon_density', 'kg m-3')
      call defvar('L_PermaResp', iou, 4, it, -c1e20, c1e20, ' ', 'F'
     &,   'perma respiration w depth'
     &,  'perma_respiraction', 'kg m-3 s-1')
      call defvar('L_Available_PC', iou, 4, it, -c1e20, c1e20, ' ', 'F'
     &,   'fraction perma C available for respiration'
     &,  'Available_Fraction', '%')

      call defvar ('L_gndtemp', iou , 4, it, -c1e20, c1e20, ' ', 'F'

     &,   'gnd temperature', 'gnd_temperature', 'C')

      call defvar('L_soilSU', iou, 4, it, -c1e20, c1e20, ' ', 'F'
     &,   'unfrozen soil moisture fraction'
     &,  'unfrozen soil moisture fraction', ' ')
      call defvar('L_soilSF', iou, 4, it, -c1e20, c1e20, ' ', 'F'
     &,   'frozen soil moisture fraction'
     &,  'frozen soil moisture fraction', ' ')

!-----------------------------------------------------------------------
!     define 4d data (x,y,type,t)
!-----------------------------------------------------------------------
      it(1) = id_xt
      it(2) = id_yt
      it(3) = id_type
      it(4) = id_time
       call defvar ('L_vegfra', iou , 4, it, -c1e20, c1e20, ' ', 'F'
     &, 'areal coverage', 'vegetation_area_fraction', '1')

!-----------------------------------------------------------------------
!     end definitions
!-----------------------------------------------------------------------
      call enddef (iou)

      return
      end

      subroutine mtlm_tavg_out (fname, ids, ide, jds, jde, imt, jmt
     &,                         POINTS, NPFT, NTYPE, xt, yt, xu
     &,                         yu, dxt, dyt, dxu, dyu, avgper, time
     &,                         stamp, land_map, TS1, CS, BURN, RESP_S
     &,                         LIT_C_T, FRAC, GPP, NPP, HT, LAI, C_VEG

     &,                         NGND, TGND, FSMC, SCHDMIN, WETMAP, PERMA
     &,                         FROZENMAP, ANN_DWET, PFTHICK, ACTIVELT
     &,                         TALIK, CS_D, RESP_SD, TRACK_PFC, SQU_PFC

     &,                         CS_P,RESP_SP, AFp

     &,                         SU, SF

     &,                         tlat, tlon, tgarea, ntrec)
!=======================================================================
!     output routine for land time averages

!     data may be sized differently in x and y from the global fields.
!     fields may be written with or without a time dimension. data
!     should be defined with the routine defvar and written with putvar.
!     if no time dimension, then data is only written once per file.
!     make sure the it, iu, ib, and ic arrays and are defining the
!     correct dimensions. ln may also need to be recalculated.

!   inputs:
!     fname        = file name
!     ids, ide ... = start and end index for data domain
!     imt, jmt ... = global array dimensions
!     xt, yt ...   = global axes
!     dxt, dyt ... = grid widths
!     avgper       = length of averaging period
!     time         = time in years
!     stamp        = time stamp
!     land_map     = land map
!     TS1, ...     = data to be written

!   outputs:
!     ntrec        = number of time record in file
!=======================================================================

      implicit none

      character(*) :: fname, stamp

      integer iou, j, ln, n, ntrec, imt, jmt, ids, ide, jds, jde, igs
      integer ige, ig, jgs, jge, jg, ils, ile, jls, jle, ib(10), ic(10)
      integer nyear, nmonth, nday, nhour, nmin, nsec
      integer POINTS, NPFT, NTYPE, land_map(imt,jmt)

      real xt(imt), xu(imt), yt(jmt), yu(jmt)
      real dxt(imt), dxu(imt), dyt(jmt), dyu(jmt)
      real tmpmask(imt,jmt), data(imt,jmt), pft(NPFT), type(NTYPE)
      real xt_e(imt+1), yt_e(jmt+1), pft_e(NPFT+1), type_e(NTYPE+1)
      real avgper, time, tmp, c0, c1, c100, c1e4, C2K
      real TS1(POINTS), CS(POINTS), RESP_S(POINTS), LIT_C_T(POINTS)
      real BURN(POINTS), FRAC(POINTS,NTYPE), GPP(POINTS,NPFT)
      real NPP(POINTS,NPFT), HT(POINTS,NPFT), LAI(POINTS,NPFT)
      real C_VEG(POINTS,NPFT)
      real tlat(ids:ide,jds:jde), tlon(ids:ide,jds:jde)
      real tgarea(ids:ide,jds:jde)

      integer NGND

      real gnd(NGND), gnd_e(NGND+1)
      real wetm(3), wetm_e(4), frozenm(7), frozenm_e(8)

      real TGND(POINTS, NGND), SCHDMIN(POINTS)
      real PERMA(POINTS), FSMC(POINTS, NPFT)
      real ACTIVELT(POINTS), TALIK(POINTS)
      real PFTHICK(POINTS), ANN_DWET(POINTS)
      real WETMAP(POINTS,3), FROZENMAP(POINTS,7)
      real TRACK_PFC(POINTS), CS_D(POINTS, NGND)
      real RESP_SD(POINTS, NGND)
      real SQU_PFC(POINTS)

      real CS_P(POINTS, NGND),RESP_SP(POINTS, NGND)
      real AFp(POINTS, NGND)

      real SU(POINTS, NGND), SF(POINTS, NGND)

      real, allocatable :: tmpij(:,:), tmpijm(:,:)
      real, allocatable :: tmpi(:), tmpj(:)
      real, allocatable :: tmpie(:), tmpje(:)

      c0 = 0.
      c1 = 1.
      c100 = 100.
      c1e4 = 1.e4
      C2K = 273.15

!-----------------------------------------------------------------------
!     open file and get latest record number
!-----------------------------------------------------------------------
      call opennext (fname, time, ntrec, iou)
      if (ntrec .le. 0) ntrec = 1

!-----------------------------------------------------------------------
!     set global write domain size (may be less than global domain)
!-----------------------------------------------------------------------
      igs = 1
      ige = imt
      if (xt(1) + 360. lt. xt(imt)) then
!       assume cyclic boundary
        igs = 2
        ige = imt-1
      endif
      ig  = ige-igs+1
      jgs = 1
      jge = jmt
      do j=2,jmt
        if (yt(j-1) .lt. -90. .and. yt(j) .gt. -90.) jgs = j
        if (yt(j-1) .lt.  90. .and. yt(j) .gt. 90.) jge = j-1
      enddo
      jg  = jge-jgs+1

!-----------------------------------------------------------------------
!     local domain size (minimum of data domain and global write domain)
!-----------------------------------------------------------------------
      ils = max(ids,igs)
      ile = min(ide,ige)
      jls = max(jds,jgs)
      jle = min(jde,jge)

      allocate ( tmpij(ils:ile,jls:jle) )
      allocate ( tmpijm(ils:ile,jls:jle) )

!-----------------------------------------------------------------------
!     write 1d data (t)
!-----------------------------------------------------------------------
      call putvars ('time', iou, ntrec, time, c1, c0)
      call rdstmp (stamp, nyear, nmonth, nday, nhour, nmin, nsec)
      call putvars ('T_avgper', iou, ntrec, avgper, c1, c0)

      if (ntrec .eq. 1) then

!-----------------------------------------------------------------------
!       write 1d data (x, y or z)
!-----------------------------------------------------------------------
        allocate ( tmpi(igs:ige) )
        allocate ( tmpj(jgs:jge) )
        allocate ( tmpie(igs:ige+1) )
        allocate ( tmpje(jgs:jge+1) )

        ib(1) = 1
        ic(1) = ig
        tmpi(igs:ige) = xt(igs:ige)
        call putvara ('longitude', iou, ig, ib, ic, tmpi, c1, c0)
        tmpi(igs:ige) = dxt(igs:ige)
        call putvara ('G_dxt', iou, ig, ib, ic, tmpi, c100, c0)

        ic(1) = jg
        tmpj(jgs:jge) = yt(jgs:jge)
        call putvara ('latitude', iou, jg, ib, ic, tmpj, c1, c0)
        tmpj(jgs:jge) = dyt(jgs:jge)
        call putvara ('G_dyt', iou, jg, ib, ic, tmpj, c100, c0)

        ic(1) = ig + 1
        call edge_maker (1, xt_e, xt, dxt, xu, dxu, imt)
        tmpie(igs:ige+1) = xt_e(igs:ige+1)
        call putvara ('longitude_edges', iou, ig+1, ib, ic, tmpie
     &,   c1, c0)

        ic(1) = jg + 1
        call edge_maker (1, yt_e, yt, dyt, yu, dyu, jmt)
        tmpje(jgs:jge+1) = yt_e(jgs:jge+1)
        call putvara ('latitude_edges', iou, jg+1, ib, ic, tmpje
     &,   c1, c0)

        deallocate ( tmpi )
        deallocate ( tmpj )
        deallocate ( tmpie )
        deallocate ( tmpje )

        do n=1, NPFT
          pft(n) = n
          pft_e(n) = pft(n) - 0.5
        enddo
        pft_e(NPFT+1) = pft(NPFT) + 0.5
        ic(1) = NPFT
        call putvara ('pft', iou, NPFT, ib, ic, pft, c1, c0)
        ic(1) = NPFT + 1
        call putvara ('pft_edges', iou, NPFT+1, ib, ic, pft_e, c1, c0)
        do n=1, NTYPE
          type(n) = float(n)
          type_e(n) = type(n) - 0.5
        enddo
        type_e(NTYPE+1) = type(NTYPE) + 0.5
        ic(1) = NTYPE
        call putvara ('type', iou, NTYPE, ib, ic, type, c1, c0)
        ic(1) = NTYPE + 1
        call putvara ('type_edges', iou, NTYPE+1, ib, ic, type_e
     &,   c1, c0)

        do n = 1, 3
          wetm(n) = n
          wetm_e(n) = wetm(n) - 0.5
        enddo
        wetm_e(4) = wetm(3) + 0.5
        ic(1) = 3
        call putvara ('wetm', iou, 3, ib, ic, wetm, c1, c0)
        ic(1) = 4
        call putvara ('wetm_edges', iou, 4, ib, ic, wetm_e, c1, c0)
        do n = 1, 7
          frozenm(n) = n
          frozenm_e(n) = frozenm(n) - 0.5
        enddo
        frozenm_e(8) = frozenm(7) + 0.5
        ic(1) = 7
        call putvara ('frozenm', iou, 7, ib, ic, frozenm, c1, c0)
        ic(1) = 8
        call putvara ('frozenm_edges', iou, 8, ib, ic, frozenm_e, c1,c0)
        do n=1, NGND
          gnd(n) = n
          gnd_e(n) = gnd(n) - 0.5
        enddo
        gnd_e(NGND+1) = gnd(NGND) + 0.5
        ic(1) = NGND
        call putvara ('gnd', iou, NGND, ib, ic, gnd, c1, c0)
        ic(1) = NGND + 1
        call putvara ('gnd_edges', iou, NGND+1,ib, ic, gnd_e
     &, c1, c0)

!-----------------------------------------------------------------------
!       write 2d data (x,y)
!-----------------------------------------------------------------------
        ib(1) = ils-igs+1
        ic(1) = ile-ils+1
        ib(2) = jls-jgs+1
        ic(2) = jle-jls+1
        ln = ic(1)*ic(2)
        tmpij(ils:ile,jls:jle) = tlat(ils:ile,jls:jle)
        call putvara ('G_latT', iou, ln, ib, ic, tmpij, c1, c0)
        tmpij(ils:ile,jls:jle) = tlon(ils:ile,jls:jle)
        call putvara ('G_lonT', iou, ln, ib, ic, tmpij, c1, c0)
        tmpij(ils:ile,jls:jle) = tgarea(ils:ile,jls:jle)
        call putvara ('G_areaT', iou, ln, ib, ic, tmpij, c1e4, c0)

      endif

!-----------------------------------------------------------------------
!     write 3d data (x,y,t)
!-----------------------------------------------------------------------
      ib(1) = ils-igs+1
      ic(1) = ile-ils+1
      ib(2) = jls-jgs+1
      ic(2) = jle-jls+1
      ib(3) = ntrec
      ic(3) = 1
      ln = ic(1)*ic(2)*ic(3)
      tmpmask(:,:) = 0.
      where (land_map(:,:) .gt. 0) tmpmask(:,:) = 1.
      tmpijm(ils:ile,jls:jle) = tmpmask(ils:ile,jls:jle)

      call unloadland (POINTS, CS, imt, jmt, land_map, data)
      tmpij(ils:ile,jls:jle) = data(ils:ile,jls:jle)
      call putvaramsk ('L_soilcarb', iou, ln, ib, ic, tmpij, tmpijm
     &, c1, c0)
      call unloadland (POINTS, RESP_S, imt, jmt, land_map, data)
      tmpij(ils:ile,jls:jle) = data(ils:ile,jls:jle)
      call putvaramsk ('L_soilresp', iou, ln, ib, ic, tmpij, tmpijm
     &, c1, c0)
      call unloadland (POINTS, LIT_C_T, imt, jmt, land_map, data)
      tmpij(ils:ile,jls:jle) = data(ils:ile,jls:jle)
      call putvaramsk ('L_veglit', iou, ln, ib, ic, tmpij, tmpijm
     &, c1, c0)
      call unloadland (POINTS, BURN, imt, jmt, land_map, data)
      tmpij(ils:ile,jls:jle) = data(ils:ile,jls:jle)
      call putvaramsk ('L_vegburn', iou, ln, ib, ic, tmpij, tmpijm
     &, c1, c0)

      call unloadland (POINTS, PERMA, imt, jmt, land_map, data)
      tmpij(ils:ile,jls:jle) = data(ils:ile,jls:jle)
      call putvaramsk ('L_perma', iou, ln, ib, ic, tmpij, tmpijm
     &, c1, c0)
      call unloadland (POINTS, ACTIVELT, imt, jmt, land_map, data)
      tmpij(ils:ile,jls:jle) = data(ils:ile,jls:jle)
      call putvaramsk ('L_activelt', iou, ln, ib, ic, tmpij, tmpijm
     &, c1, c0)
      call unloadland (POINTS, TALIK, imt, jmt, land_map, data)
      tmpij(ils:ile,jls:jle) = data(ils:ile,jls:jle)
      call putvaramsk ('L_talik', iou, ln, ib, ic, tmpij, tmpijm
     &, c1, c0)
      call unloadland (POINTS, ANN_DWET, imt, jmt, land_map, data)
      tmpij(ils:ile,jls:jle) = data(ils:ile,jls:jle)
      call putvaramsk ('L_dwet', iou, ln, ib, ic, tmpij, tmpijm
     &, c1, c0)
      call unloadland (POINTS, PFTHICK, imt, jmt, land_map, data)
      tmpij(ils:ile,jls:jle) = data(ils:ile,jls:jle)
      call putvaramsk ('L_PFthick', iou, ln, ib, ic, tmpij, tmpijm
     &, c1, c0)
      call unloadland (POINTS, SCHDMIN, imt, jmt, land_map, data)
      tmpij(ils:ile,jls:jle) = data(ils:ile,jls:jle)
      call putvaramsk ('L_schDmin', iou, ln, ib, ic, tmpij, tmpijm
     &, c1, c0)
      call unloadland (POINTS, TRACK_PFC, imt, jmt, land_map, data)
      tmpij(ils:ile,jls:jle) = data(ils:ile,jls:jle)
      call putvaramsk ('L_permacarb', iou, ln, ib, ic, tmpij, tmpijm
     &, c1, c0)
      call unloadland (POINTS, SQU_PFC, imt, jmt, land_map, data)
      tmpij(ils:ile,jls:jle) = data(ils:ile,jls:jle)
      call putvaramsk ('L_squedcarb', iou, ln, ib, ic, tmpij, tmpijm
     &, c1, c0)

!-----------------------------------------------------------------------
!     write 4d data (x,y,pft,t)
!-----------------------------------------------------------------------
      ib(1) = ils-igs+1
      ic(1) = ile-ils+1
      ib(2) = jls-jgs+1
      ic(2) = jle-jls+1
      ic(3) = 1
      ib(4) = ntrec
      ic(4) = 1
      ln = ic(1)*ic(2)*ic(3)*ic(4)
      do n=1,npft
        ib(3) = n
        call unloadland (POINTS, GPP(1,n), imt, jmt, land_map, data)
        tmpij(ils:ile,jls:jle) = data(ils:ile,jls:jle)
        call putvaramsk ('L_veggpp', iou, ln, ib, ic, tmpij, tmpijm
     &,   c1, c0)
        call unloadland (POINTS, NPP(1,n), imt, jmt, land_map, data)
        tmpij(ils:ile,jls:jle) = data(ils:ile,jls:jle)
        call putvaramsk ('L_vegnpp', iou, ln, ib, ic, tmpij, tmpijm
     &,   c1, c0)
        call unloadland (POINTS, HT(1,n), imt, jmt, land_map, data)
        tmpij(ils:ile,jls:jle) = data(ils:ile,jls:jle)
        call putvaramsk ('L_veghgt', iou, ln, ib, ic, tmpij, tmpijm
     &,   c1, c0)
        call unloadland (POINTS, LAI(1,n), imt, jmt, land_map, data)
        tmpij(ils:ile,jls:jle) = data(ils:ile,jls:jle)
        call putvaramsk ('L_veglai', iou, ln, ib, ic, tmpij, tmpijm
     &,   c1, c0)
        call unloadland (POINTS, C_VEG(1,n), imt, jmt, land_map, data)
        tmpij(ils:ile,jls:jle) = data(ils:ile,jls:jle)
        call putvaramsk ('L_vegcarb', iou, ln, ib, ic, tmpij, tmpijm
     &,   c1, c0)

        call unloadland (POINTS, FSMC(1,n), imt, jmt, land_map, data)
        tmpij(ils:ile,jls:jle) = data(ils:ile,jls:jle)
        call putvaramsk ('L_vegfsmc', iou, ln, ib, ic, tmpij, tmpijm
     &,   c1, c0)

      enddo

!-----------------------------------------------------------------------
!     write 4d data (x,y,wetm,t)
!-----------------------------------------------------------------------
      ib(1) = ils-igs+1
      ic(1) = ile-ils+1
      ib(2) = jls-jgs+1
      ic(2) = jle-jls+1
      ic(3) = 1
      ib(4) = ntrec
      ic(4) = 1
      ln = ic(1)*ic(2)*ic(3)*ic(4)
      do n=1,3
        ib(3) = n
        call unloadland (POINTS,WETMAP(1,n), imt, jmt, land_map, data)
        tmpij(ils:ile,jls:jle) = data(ils:ile,jls:jle)
        call putvaramsk ('L_wetmap', iou, ln, ib, ic, tmpij, tmpijm
     &,   c1, c0)
      enddo

!-----------------------------------------------------------------------
!     write 4d data (x,y,frozenm,t)
!-----------------------------------------------------------------------
      ib(1) = ils-igs+1
      ic(1) = ile-ils+1
      ib(2) = jls-jgs+1
      ic(2) = jle-jls+1
      ic(3) = 1
      ib(4) = ntrec
      ic(4) = 1
      ln = ic(1)*ic(2)*ic(3)*ic(4)
      do n=1,7
        ib(3) = n
        call unloadland (POINTS,FROZENMAP(1,n),imt, jmt, land_map, data)
        tmpij(ils:ile,jls:jle) = data(ils:ile,jls:jle)
        call putvaramsk ('L_frozenmap', iou, ln, ib, ic, tmpij, tmpijm
     &,   c1, c0)
      enddo

!-----------------------------------------------------------------------
!     write 4d data (x,y,gnd,t)
!-----------------------------------------------------------------------
      ib(1) = ils-igs+1
      ic(1) = ile-ils+1
      ib(2) = jls-jgs+1
      ic(2) = jle-jls+1
      ic(3) = 1
      ib(4) = ntrec
      ic(4) = 1
      ln = ic(1)*ic(2)*ic(3)*ic(4)
      do n=1,ngnd
        ib(3) = n
        call unloadland (POINTS,TGND(1,n),imt,jmt,land_map, data)
        tmpij(ils:ile,jls:jle) = data(ils:ile,jls:jle)
        call putvaramsk ('L_gndtemp', iou,ln,ib,ic,tmpij,tmpijm

     &, c1, C2K)

        call unloadland (POINTS, CS_D(1,n), imt, jmt, land_map,data)
        tmpij(ils:ile,jls:jle) = data(ils:ile,jls:jle)
        call putvaramsk ('L_soilcarbD', iou, ln, ib, ic, tmpij, tmpijm
     &,   c1, c0)
        call unloadland (POINTS, RESP_SD(1,n), imt, jmt, land_map,data)
        tmpij(ils:ile,jls:jle) = data(ils:ile,jls:jle)
        call putvaramsk ('L_soilrespD', iou, ln, ib, ic, tmpij, tmpijm
     &,   c1, c0)

        call unloadland (POINTS, CS_P(1,n), imt, jmt, land_map,data)
        tmpij(ils:ile,jls:jle) = data(ils:ile,jls:jle)
        call putvaramsk ('L_PermaCarbD', iou, ln, ib, ic, tmpij, tmpijm
     &,   c1, c0)
        call unloadland (POINTS, RESP_SP(1,n), imt, jmt, land_map,data)
        tmpij(ils:ile,jls:jle) = data(ils:ile,jls:jle)
        call putvaramsk ('L_PermaResp', iou, ln, ib, ic, tmpij, tmpijm
     &,   c1, c0)
        call unloadland (POINTS, AFp(1,n), imt, jmt, land_map,data)
        tmpij(ils:ile,jls:jle) = data(ils:ile,jls:jle)
        call putvaramsk ('L_Available_PC', iou, ln, ib, ic, tmpij, tmpijm
     &,   c1, c0)

        call unloadland (POINTS, SU(1,n), imt, jmt, land_map, data)
        tmpij(ils:ile,jls:jle) = data(ils:ile,jls:jle)
        call putvaramsk ('L_soilSU', iou, ln, ib, ic, tmpij, tmpijm
     &,   c1, c0)
        call unloadland (POINTS, SF(1,n), imt, jmt, land_map, data)
        tmpij(ils:ile,jls:jle) = data(ils:ile,jls:jle)
        call putvaramsk ('L_soilSF', iou, ln, ib, ic, tmpij, tmpijm
     &,   c1, c0)

      enddo

!-----------------------------------------------------------------------
!     write 4d data (x,y,type,t)
!-----------------------------------------------------------------------
      ib(1) = ils-igs+1
      ic(1) = ile-ils+1
      ib(2) = jls-jgs+1
      ic(2) = jle-jls+1
      ic(3) = 1
      ib(4) = ntrec
      ic(4) = 1
      do n=1,ntype
        ib(3) = n
        call unloadland (POINTS, FRAC(1,n), imt, jmt, land_map, data)
        tmpij(ils:ile,jls:jle) = data(ils:ile,jls:jle)
        where (tmpij(ils:ile,jls:jle) .lt. 0.)
     &    tmpij(ils:ile,jls:jle) = c0
        call putvaramsk ('L_vegfra', iou, ln, ib, ic, tmpij, tmpijm
     &,   c1, c0)
      enddo

      deallocate ( tmpij )
      deallocate ( tmpijm )

      return
      end

