changeset 56:d6f060519f45

Added diagnostic output features: + separate off-line sediment output configuration (namelist group &DIASEDOFFL); + added support for decadal, centenial and millenial output (for off- and online). Compiles and runs, but didn't check if iteration works well, and there is a runtime variable initialisation issue (if that debug option is enabled).
author Marco van Hulten <marco@hulten.org>
date Fri, 14 Dec 2018 14:06:54 +0100
parents 373476c063cd
children 354c3347b4a3
files common_bgcs.h90 mo_bgcmean.F90 mo_sedmnt_offline.F90 mod_dia.F ncout_hamocc.F ncout_onlysed.F90
diffstat 6 files changed, 431 insertions(+), 266 deletions(-) [+]
line wrap: on
line diff
--- a/common_bgcs.h90	Thu Dec 06 14:29:09 2018 +0100
+++ b/common_bgcs.h90	Fri Dec 14 14:06:54 2018 +0100
@@ -5,9 +5,9 @@
 ! Also, they written/read to and from restart files.
 ! There are probably more efficient/elegant solutions.
 
-real :: sedlay2,powtra2,burial2
-
 common /bgcs_hamocc/                                                &
    &   sedlay2(idm,jdm,2*ks,nsedtra) ,powtra2(idm,jdm,2*ks,npowtra) &
    &  ,burial2(idm,jdm,2,   nsedtra)
 
+real :: sedlay2,powtra2,burial2
+
--- a/mo_bgcmean.F90	Thu Dec 06 14:29:09 2018 +0100
+++ b/mo_bgcmean.F90	Fri Dec 14 14:06:54 2018 +0100
@@ -13,6 +13,8 @@
 !        index arrays
 !     Tjiputra          *UNI-RESEARCH    25.11.15
 !      - added natural DIC/ALK/CALC/OMEGAC variables 
+!     Marco van Hulten  *GFI, UiB        28.11.18
+!      - added support for decadal, centenial and millenial output
 !  
 !     Purpose
 !     -------
@@ -35,8 +37,9 @@
       INTEGER, PARAMETER :: nbgcmax=10
       REAL, DIMENSION(nbgcmax), SAVE ::  diagfq_bgc,filefq_bgc
       INTEGER, DIMENSION(nbgcmax), SAVE :: nacc_bgc
-      LOGICAL, DIMENSION(nbgcmax), SAVE :: diagmon_bgc,diagann_bgc,     &
-     &  filemon_bgc,fileann_bgc,bgcwrt
+      LOGICAL, DIMENSION(nbgcmax), SAVE :: bgcwrt,                      &
+         & diagmon_bgc,diagann_bgc,diagdec_bgc,diagcen_bgc,diagmil_bgc, &
+         & filemon_bgc,fileann_bgc,filedec_bgc,filecen_bgc,filemil_bgc
  
 ! --- Namelist for diagnostic output 
       INTEGER, DIMENSION(nbgcmax), SAVE ::                              &
@@ -403,10 +406,19 @@
         ENDIF
         diagmon_bgc(n)=.false.
         diagann_bgc(n)=.false.
+        diagdec_bgc(n)=.false.
+        diagcen_bgc(n)=.false.
+        diagmil_bgc(n)=.false.
         IF (GLB_AVEPERIO(n).EQ.30) THEN
           diagmon_bgc(n)=.true.
         ELSEIF (GLB_AVEPERIO(n).EQ.365) THEN
           diagann_bgc(n)=.true.
+        ELSEIF (GLB_AVEPERIO(n).EQ.3650) THEN
+          diagdec_bgc(n)=.true.
+        ELSEIF (GLB_AVEPERIO(n).EQ.36500) THEN
+          diagcen_bgc(n)=.true.
+        ELSEIF (GLB_AVEPERIO(n).EQ.365000) THEN
+          diagmil_bgc(n)=.true.
         ENDIF
         IF (GLB_FILEFREQ(n).LT.0) THEN
           filefq_bgc(n)=-real(nstepinday)/GLB_FILEFREQ(n)
@@ -415,10 +427,19 @@
         ENDIF
         filemon_bgc(n)=.false.
         fileann_bgc(n)=.false.
+        filedec_bgc(n)=.false.
+        filecen_bgc(n)=.false.
+        filemil_bgc(n)=.false.
         IF (GLB_FILEFREQ(n).EQ.30) THEN
           filemon_bgc(n)=.true.
         ELSEIF (GLB_FILEFREQ(n).EQ.365) THEN
           fileann_bgc(n)=.true.
+        ELSEIF (GLB_FILEFREQ(n).EQ.3650) THEN
+          filedec_bgc(n)=.true.
+        ELSEIF (GLB_FILEFREQ(n).EQ.36500) THEN
+          filecen_bgc(n)=.true.
+        ELSEIF (GLB_FILEFREQ(n).EQ.365000) THEN
+          filemil_bgc(n)=.true.
         ENDIF
       ENDDO
 
--- a/mo_sedmnt_offline.F90	Thu Dec 06 14:29:09 2018 +0100
+++ b/mo_sedmnt_offline.F90	Fri Dec 14 14:06:54 2018 +0100
@@ -73,13 +73,15 @@
 use mo_control_bgc
 use mo_carbch,      only: ocetra, keqb, co3
 use mo_param1_bgc,  only: nocetra, idet, icalc, iopal, ifdust        &
-   &                    , nsedtra, npowtra
+   &                    , ks, nsedtra, npowtra
 use mo_biomod,      only: kbo, bolay, wpoc, wmin, wmax, wdust        &
    &                    , wlin, wopal, wcal
 use mod_xc
 use mo_common_bgc
-use mo_bgcmean
-
+use mod_dia,        only: ddm,depthslev,depthslev_bnds,nstepinday,pbath &
+   &                    , diafnm,iotype
+use mod_nctools
+use mo_bgcmean,     only: nacc_bgc
 
 implicit none
 
@@ -121,11 +123,12 @@
    &                                 diagmil_sed,                 &
    &                                 filemon_sed, fileann_sed,    &
    &                                 filedec_sed, filecen_sed,    &
-   &                                 filemil_sed, sedwrt
+   &                                 filemil_sed
 
 ! namelist for diagnostic output
 !
 integer, dimension(nsedmax), save ::                                &
+   & CARFLX_BOT    =0    ,BSIFLX_BOT    =0    ,CALFLX_BOT    =0  ,  &
    & SDM_POWAIC    =0    ,SDM_POWAAL    =0    ,SDM_POWAPH    =0  ,  &
    & SDM_POWAOX    =0    ,SDM_POWN2     =0    ,SDM_POWNO3    =0  ,  &
    & SDM_POWASI    =0    ,SDM_SSSO12    =0    ,SDM_SSSSIL    =0  ,  &
@@ -135,7 +138,8 @@
    & GLB_AVEPERIO  =0    ,GLB_FILEFREQ  =0    ,GLB_COMPFLAG  =0  ,  &
    & GLB_NCFORMAT  =0    ,GLB_INVENTORY =0
 character(len=10), dimension(nsedmax), save :: GLB_FNAMETAG
-namelist /DIASED/                                                   &
+namelist /DIASEDOFFL/                                               &
+   & CARFLX_BOT        ,BSIFLX_BOT        ,CALFLX_BOT        ,      &
    & SDM_POWAIC        ,SDM_POWAAL        ,SDM_POWAPH        ,      &
    & SDM_POWAOX        ,SDM_POWN2         ,SDM_POWNO3        ,      &
    & SDM_POWASI        ,SDM_SSSO12        ,SDM_SSSSIL        ,      &
@@ -145,6 +149,49 @@
    & GLB_AVEPERIO      ,GLB_FILEFREQ      ,GLB_COMPFLAG      ,      &
    & GLB_NCFORMAT      ,GLB_FNAMETAG      ,GLB_INVENTORY
 
+! bottom fluxes
+!
+integer, save :: i_bscoffl_m2d
+integer, dimension(nsedmax), save ::                                 &
+   &          jcarflx_bot = 0 ,                                      &
+   &          jbsiflx_bot = 0 ,                                      &
+   &          jcalflx_bot = 0
+
+integer, save :: nsedm2d
+
+! sediment
+!
+integer, save :: i_bscoffl_sed
+integer, dimension(nsedmax), save ::                                 &
+   &          jpowaic = 0 ,                                          &
+   &          jpowaal = 0 ,                                          &
+   &          jpowaph = 0 ,                                          &
+   &          jpowaox = 0 ,                                          &
+   &          jpown2  = 0 ,                                          &
+   &          jpowno3 = 0 ,                                          &
+   &          jpowasi = 0 ,                                          &
+   &          jssso12 = 0 ,                                          &
+   &          jssssil = 0 ,                                          &
+   &          jsssc12 = 0 ,                                          &
+   &          jssster = 0
+
+integer, save :: nsedt_sed
+
+! burial
+!
+integer, save :: i_bscoffl_bur
+integer, dimension(nsedmax), save ::                                 &
+   &          jburssso12 = 0 ,                                       &
+   &          jbursssc12 = 0 ,                                       &
+   &          jburssssil = 0 ,                                       &
+   &          jburssster = 0
+
+integer, save :: nsedt_bur
+
+real, dimension (:,:,:),   allocatable :: sedm2d
+real, dimension (:,:,:,:), allocatable :: sedt_sed
+real, dimension (:,:,:),   allocatable :: sedt_bur
+
 ! private variables and subprograms
 !
 integer, private :: i,j,iocetra,n
@@ -373,6 +420,13 @@
       nacc_bgc_1        = nacc_bgc(1)
       nacc_bgc_2        = nacc_bgc(2)
 
+      ! initialise #ts for every group to zero
+      do n = 1, nsedmax
+         if (GLB_AVEPERIO(n) /= 0) then
+            nacc_sed(n) = 0
+         endif
+      enddo
+
       do nyear = 1, maxyear
          if (mnproc == 1) write(io_stdo_bgc,'(a,i6)')                      &
                &         'sedmnt_offline(): nyear_global = ', nyear_global
@@ -393,39 +447,30 @@
                & silpro_clim(:,:,nmonth), produs_clim(:,:,nmonth),         &
                & co3_kbo_clim(:,:,nmonth))
 
-!TODO: 1. update sedwrt() instead of setting..? (inner/outer loop)
-!      2. inner/outer loop
-!      3. no nstep, nday_... but use nyear and nmonth counters.
-
             ! write monthly output fields
             do n = 1, nsed
-               nacc_sed(n) = 1 ! mo (#timesteps)
-               if diagmon_sed(n) then
-                  if (GLB_INVENTORY(n) /= 0) then                          &
+               nacc_sed(n) = nacc_sed(n) + 1 ! mo (#timesteps)
+               if ( diagmon_sed(n) ) then
+                  if (GLB_INVENTORY(n) /= 0) then
                      call INVENTORY_BGC(kpie,kpje,kpke,pdlxp,pdlyp,pddpo,omask,0)
                   endif
-                  call ncwrt_onlysed(n) ! ncout_onlysed.F90
+                  call ncwrt_onlysed(n)
                   nacc_sed(n)=0
                endif
             enddo
 
          enddo
 
-         !FIXME: Following not very useful because no special code in ncout_onlysed.F90.
-         !       Moreover, now all other cases are excluded.  Better: something with fq.
-         !TODO:  Can also combine below in one loop / if .. or .. or ... -> DONE.
-
          ! write yearly and multi-yearly output fields
          do n = 1, nsed
-            nacc_sed(n) = nacc_sed(n) + 12 ! mo (#timesteps)
-            if ( (diagann_sed(n)                              .or.      &
+            if ( (diagann_sed(n)                            ) .or.      &
                & (diagdec_sed(n) .and. mod(nyear,   10) == 0) .or.      &
                & (diagcen_sed(n) .and. mod(nyear,  100) == 0) .or.      &
                & (diagmil_sed(n) .and. mod(nyear, 1000) == 0) ) then
-               if (GLB_INVENTORY(n) /= 0) then                          &
+               if (GLB_INVENTORY(n) /= 0) then
                   call INVENTORY_BGC(kpie,kpje,kpke,pdlxp,pdlyp,pddpo,omask,0)
                endif
-               call ncwrt_onlysed(n) ! ncout_onlysed.F90
+               call ncwrt_onlysed(n)
                nacc_sed(n) = 0
             endif
          enddo
@@ -463,25 +508,44 @@
 
 subroutine alloc_mem_sedmnt_offline(kpie, kpje)
    integer :: kpie, kpje
-   integer :: errstat
+   integer :: errstat, iounit
+   logical :: isopen, exists
+
+!  Read namelist for diagnostic output
+   GLB_AVEPERIO = 0
+   do iounit = 10, 99
+      inquire(unit=iounit, opened=isopen)
+      if (.not.isopen) EXIT
+   enddo
+   inquire(file='ocn_in', exist=exists)
+   if (exists) then
+      open(iounit, file='ocn_in', status='old', action='read', recl=80)
+   else
+      inquire(file='limits',exist=exists)
+      if (exists) then
+         open(iounit, file='limits', status='old', action='read', recl=80)
+      else
+        stop 'cannot find limits file'
+      endif
+   endif
+   read(iounit, nml=DIASEDOFFL)
+   close(iounit)
 
    ! determine number of output groups
    nsed = 0
    do n = 1, nsedmax
       if (GLB_AVEPERIO(n) /= 0) then
          nsed = nsed + 1
-         nacc_sed(n) = 0
       endif
    enddo
 
    do n = 1, nsed
       GLB_FILEFREQ(n) = max(GLB_AVEPERIO(n), GLB_FILEFREQ(n))
-      !diagfq_sed(n) = GLB_AVEPERIO(n)/30 !FIXME: Do we want to use it?
       diagmon_sed(n) = .false.
       diagann_sed(n) = .false.
-      diagdec_sed(n) = .false. ! isn't actually used now
-      diagcen_sed(n) = .false. ! isn't actually used now
-      diagmil_sed(n) = .false. ! isn't actually used now
+      diagdec_sed(n) = .false.
+      diagcen_sed(n) = .false.
+      diagmil_sed(n) = .false.
       if (GLB_AVEPERIO(n) == 30) then
          diagmon_sed(n) = .true.
       elseif (GLB_AVEPERIO(n) == 365) then
@@ -507,83 +571,106 @@
          filemon_sed(n) = .true.
       elseif (GLB_FILEFREQ(n) == 365) then
          fileann_sed(n) = .true.
-      elseif (GLB_FILEPERIO(n) == 3650) then
+      elseif (GLB_FILEFREQ(n) == 3650) then
          filedec_sed(n) = .true.
-      elseif (GLB_FILEPERIO(n) == 36500) then
+      elseif (GLB_FILEFREQ(n) == 36500) then
          filecen_sed(n) = .true.
-      elseif (GLB_FILEPERIO(n) == 365000) then
+      elseif (GLB_FILEFREQ(n) == 365000) then
          filemil_sed(n) = .true.
       endif
    enddo
 
-   i_bsc_sed = 0
+   ! Re-define index variables according to namelist
+   i_bscoffl_m2d = 0
+   do n = 1,nsed
+     if (CARFLX_BOT(n) > 0) i_bscoffl_m2d = i_bscoffl_m2d+1
+     jcarflx_bot(n) = i_bscoffl_m2d*min(1,CARFLX_BOT(n))
+     if (BSIFLX_BOT(n) > 0) i_bscoffl_m2d = i_bscoffl_m2d+1
+     jbsiflx_bot(n) = i_bscoffl_m2d*min(1,BSIFLX_BOT(n))
+     if (CALFLX_BOT(n) > 0) i_bscoffl_m2d = i_bscoffl_m2d+1
+     jcalflx_bot(n) = i_bscoffl_m2d*min(1,CALFLX_BOT(n))
+   enddo
+
+   i_bscoffl_sed = 0
    do n = 1,nsed
-      if (SDM_POWAIC(n) > 0) i_bsc_sed = i_bsc_sed+1
-      jpowaic(n) = i_bsc_sed*min(1,SDM_POWAIC(n))
-      if (SDM_POWAAL(n) > 0) i_bsc_sed = i_bsc_sed+1
-      jpowaal(n) = i_bsc_sed*min(1,SDM_POWAAL(n))
-      if (SDM_POWAPH(n) > 0) i_bsc_sed = i_bsc_sed+1
-      jpowaph(n) = i_bsc_sed*min(1,SDM_POWAPH(n))
-      if (SDM_POWAOX(n) > 0) i_bsc_sed = i_bsc_sed+1
-      jpowaox(n) = i_bsc_sed*min(1,SDM_POWAOX(n))
-      if (SDM_POWN2(n)  > 0) i_bsc_sed = i_bsc_sed+1
-      jpown2(n)  = i_bsc_sed*min(1,SDM_POWN2(n))
-      if (SDM_POWNO3(n) > 0) i_bsc_sed = i_bsc_sed+1
-      jpowno3(n) = i_bsc_sed*min(1,SDM_POWNO3(n))
-      if (SDM_POWASI(n) > 0) i_bsc_sed = i_bsc_sed+1
-      jpowasi(n) = i_bsc_sed*min(1,SDM_POWASI(n))
-      if (SDM_SSSO12(n) > 0) i_bsc_sed = i_bsc_sed+1
-      jssso12(n) = i_bsc_sed*min(1,SDM_SSSO12(n))
-      if (SDM_SSSSIL(n) > 0) i_bsc_sed = i_bsc_sed+1
-      jssssil(n) = i_bsc_sed*min(1,SDM_SSSSIL(n))
-      if (SDM_SSSC12(n) > 0) i_bsc_sed = i_bsc_sed+1
-      jsssc12(n) = i_bsc_sed*min(1,SDM_SSSC12(n))
-      if (SDM_SSSTER(n) > 0) i_bsc_sed = i_bsc_sed+1
-      jssster(n) = i_bsc_sed*min(1,SDM_SSSTER(n))
+      if (SDM_POWAIC(n) > 0) i_bscoffl_sed = i_bscoffl_sed+1
+      jpowaic(n) = i_bscoffl_sed*min(1,SDM_POWAIC(n))
+      if (SDM_POWAAL(n) > 0) i_bscoffl_sed = i_bscoffl_sed+1
+      jpowaal(n) = i_bscoffl_sed*min(1,SDM_POWAAL(n))
+      if (SDM_POWAPH(n) > 0) i_bscoffl_sed = i_bscoffl_sed+1
+      jpowaph(n) = i_bscoffl_sed*min(1,SDM_POWAPH(n))
+      if (SDM_POWAOX(n) > 0) i_bscoffl_sed = i_bscoffl_sed+1
+      jpowaox(n) = i_bscoffl_sed*min(1,SDM_POWAOX(n))
+      if (SDM_POWN2(n)  > 0) i_bscoffl_sed = i_bscoffl_sed+1
+      jpown2(n)  = i_bscoffl_sed*min(1,SDM_POWN2(n))
+      if (SDM_POWNO3(n) > 0) i_bscoffl_sed = i_bscoffl_sed+1
+      jpowno3(n) = i_bscoffl_sed*min(1,SDM_POWNO3(n))
+      if (SDM_POWASI(n) > 0) i_bscoffl_sed = i_bscoffl_sed+1
+      jpowasi(n) = i_bscoffl_sed*min(1,SDM_POWASI(n))
+      if (SDM_SSSO12(n) > 0) i_bscoffl_sed = i_bscoffl_sed+1
+      jssso12(n) = i_bscoffl_sed*min(1,SDM_SSSO12(n))
+      if (SDM_SSSSIL(n) > 0) i_bscoffl_sed = i_bscoffl_sed+1
+      jssssil(n) = i_bscoffl_sed*min(1,SDM_SSSSIL(n))
+      if (SDM_SSSC12(n) > 0) i_bscoffl_sed = i_bscoffl_sed+1
+      jsssc12(n) = i_bscoffl_sed*min(1,SDM_SSSC12(n))
+      if (SDM_SSSTER(n) > 0) i_bscoffl_sed = i_bscoffl_sed+1
+      jssster(n) = i_bscoffl_sed*min(1,SDM_SSSTER(n))
    enddo
 
-   i_bsc_bur = 0
+   i_bscoffl_bur = 0
    do n = 1,nsed
-      if (BUR_SSSO12(n) > 0) i_bsc_bur = i_bsc_bur+1
-      jburssso12(n) = i_bsc_bur*min(1,BUR_SSSO12(n))
-      if (BUR_SSSC12(n) > 0) i_bsc_bur = i_bsc_bur+1
-      jbursssc12(n) = i_bsc_bur*min(1,BUR_SSSC12(n))
-      if (BUR_SSSSIL(n) > 0) i_bsc_bur = i_bsc_bur+1
-      jburssssil(n) = i_bsc_bur*min(1,BUR_SSSSIL(n))
-      if (BUR_SSSTER(n) > 0) i_bsc_bur = i_bsc_bur+1
-      jburssster(n) = i_bsc_bur*min(1,BUR_SSSTER(n))
+      if (BUR_SSSO12(n) > 0) i_bscoffl_bur = i_bscoffl_bur+1
+      jburssso12(n) = i_bscoffl_bur*min(1,BUR_SSSO12(n))
+      if (BUR_SSSC12(n) > 0) i_bscoffl_bur = i_bscoffl_bur+1
+      jbursssc12(n) = i_bscoffl_bur*min(1,BUR_SSSC12(n))
+      if (BUR_SSSSIL(n) > 0) i_bscoffl_bur = i_bscoffl_bur+1
+      jburssssil(n) = i_bscoffl_bur*min(1,BUR_SSSSIL(n))
+      if (BUR_SSSTER(n) > 0) i_bscoffl_bur = i_bscoffl_bur+1
+      jburssster(n) = i_bscoffl_bur*min(1,BUR_SSSTER(n))
    enddo
 
-   nsedt_sed  = i_bsc_sed
-   nsedt_bur  = i_bsc_bur
+   nsedm2d    = i_bscoffl_m2d
+   nsedt_sed  = i_bscoffl_sed
+   nsedt_bur  = i_bscoffl_bur
 
-   if (mnproc == 1) THEN
-        write(io_stdo_bgc,*)'Memory allocation for variable sedtsed ...'
-        write(io_stdo_bgc,*)'First dimension    : ',kpie
-        write(io_stdo_bgc,*)'Second dimension   : ',kpje
-        write(io_stdo_bgc,*)'Third dimension    : ',ks
-        write(io_stdo_bgc,*)'Forth dimension    : ',nsedt_sed
+   if (mnproc == 1) then
+     write(io_stdo_bgc,*)'Memory allocation for variable sedm2d ...'
+     write(io_stdo_bgc,*)'First dimension    : ',kpie
+     write(io_stdo_bgc,*)'Second dimension   : ',kpje
+     write(io_stdo_bgc,*)'Forth dimension    : ',nsedm2d
    endif
 
-   allocate (sedt_sed(1-nbdy:kpie+nbdy,1-nbdy:kpje+nbdy,ks,          &
-      &      nsedt_sed),stat=errstat)
-   if (errstat /= 0) STOP 'not enough memory sedt_sed'
+   allocate(sedm2d(1-nbdy:kpie+nbdy,1-nbdy:kpje+nbdy,nsedm2d),stat=errstat)
+   if (errstat /= 0) stop 'not enough memory sedm2d'
+   if (nsedt_sed /= 0) sedm2d = 0.
+
+   if (mnproc == 1) then
+     write(io_stdo_bgc,*)'Memory allocation for variable sedt_sed ...'
+     write(io_stdo_bgc,*)'First dimension    : ',kpie
+     write(io_stdo_bgc,*)'Second dimension   : ',kpje
+     write(io_stdo_bgc,*)'Third dimension    : ',ks
+     write(io_stdo_bgc,*)'Forth dimension    : ',nsedt_sed
+   endif
+
+   allocate(sedt_sed(1-nbdy:kpie+nbdy,1-nbdy:kpje+nbdy,ks,nsedt_sed),stat=errstat)
+   if (errstat /= 0) stop 'not enough memory sedt_sed'
    if (nsedt_sed /= 0) sedt_sed = 0.
 
-   if (mnproc == 1) THEN
-        write(io_stdo_bgc,*)'Memory allocation for variable sedtbur ...'
-        write(io_stdo_bgc,*)'First dimension    : ',kpie
-        write(io_stdo_bgc,*)'Second dimension   : ',kpje
-        write(io_stdo_bgc,*)'Third dimension    : ',nsedt_bur
+   if (mnproc == 1) then
+     write(io_stdo_bgc,*)'Memory allocation for variable sedt_bur ...'
+     write(io_stdo_bgc,*)'First dimension    : ',kpie
+     write(io_stdo_bgc,*)'Second dimension   : ',kpje
+     write(io_stdo_bgc,*)'Third dimension    : ',nsedt_bur
    endif
 
-   allocate (sedt_bur(1-nbdy:kpie+nbdy,1-nbdy:kpje+nbdy,             &
-      &      nsedt_bur),stat=errstat)
-   if (errstat /= 0) STOP 'not enough memory sedt_sed'
+   allocate(sedt_bur(1-nbdy:kpie+nbdy,1-nbdy:kpje+nbdy,nsedt_bur),stat=errstat)
+   if (errstat /= 0) stop 'not enough memory sedt_bur'
    if (nsedt_bur /= 0) sedt_bur = 0.
 
+
    !-- In HAMOCC, above is done by ALLOC_MEM_BGCMEAN() below by ALLOC_MEM_SEDMNT()
 
+
    if (mnproc == 1) then
       write(io_stdo_bgc,*)' '
       write(io_stdo_bgc,*)'Memory allocation for spin-up-specific sediment modules:'
@@ -697,5 +784,228 @@
 
 end subroutine alloc_mem_sedmnt_offline
 
+subroutine ncwrt_onlysed(iogrp)
+
+!-----------------------------------------------------------------------
+!
+! Write diagnostic sediment fields during off-line sediment spin-up
+!
+! Copyright (C) 2018 Marco van Hulten <Marco.Hulten@uib.no> et al.
+!                    Geophysical Institute @ University of Bergen
+!
+! This subroutine is free software: you can redistribute it and/or modify
+! it under the terms of the GNU General Public License as published by
+! the Free Software Foundation, either version 3 of the License, or
+! (at your option) any later version.
+!
+!
+! Method:
+!
+! History:
+!
+! This is based on ncout_hamocc.F that contains ncwrt_bgc().
+!
+! Marco van Hulten <Marco.Hulten@uib.no>           2018-11-20
+! - Initial code based on ncout_hamocc.F.
+! Marco van Hulten <Marco.Hulten@uib.no>           2018-11-29
+! - Added a second routine vardef_onlysed() to this file akin to
+!   ncout_hamocc.F.
+!
+!-----------------------------------------------------------------------
+
+use mo_bgcmean    , only: inisdm, inibur, accsdm, accbur, wrtsdm, wrtbur, logsdm
+
+implicit none
+
+!#include "common_clndr.h90"
+#include "common_blocks.h90"
+
+integer, intent(in) :: iogrp
+
+integer i,j,k,l,nt
+integer nhour,ny,nm,nd,dayfrac,irec(nsedmax),cmpflg
+character(len= 2) seqstring
+character(len=80) fname(nsedmax)
+character(len=20) startdate
+character(len=30) timeunits
+real datenum,rnacc
+logical append2file(nsedmax)
+data append2file /nsedmax*.false./
+save fname,irec,append2file
+
+! set time information
+timeunits=' '
+startdate=' '
+fname=' '
+write(timeunits,'(a11,i4.4,a1,i2.2,a1,i2.2,a6)')                           &
+   &            'days since ',min(1800,nyear0),'-',1,'-',1,' 00:00'
+write(startdate,'(i4.4,a1,i2.2,a1,i2.2,a6)')                               &
+   &            nyear0,'-',nmonth0,'-',nday0,' 00:00'
+datenum=time-time0-0.5*diagfq_sed(iogrp)/nstep_in_day
+
+! get file name
+write (seqstring,'(I0.2)') nburst
+if (.not.append2file(iogrp)) then
+   call diafnm(runid,runid_len,expcnf,trim(GLB_FNAMETAG(iogrp))//"."//seqstring,nstep, &
+      &        filefq_sed(iogrp)/real(nstep_in_day),filemon_sed(iogrp),       &
+      &        fileann_sed(iogrp),filedec_sed(iogrp),filecen_sed(iogrp),      &
+      &        filemil_sed(iogrp),fname(iogrp))
+   append2file(iogrp)=.true.
+   irec(iogrp)=1
+else
+   irec(iogrp)=irec(iogrp)+1
+endif
+if ( (filemil_sed(iogrp) .or. filecen_sed(iogrp) .or. filemil_sed(iogrp) .or. &
+   &  fileann_sed(iogrp) .or. filemon_sed(iogrp))                             &
+   & .or. .not.(fileann_sed(iogrp) .or. filemon_sed(iogrp)) .and.             &
+   &  mod(nstep+.5,filefq_sed(iogrp))<2.) then
+   append2file(iogrp) = .false.
+endif
+
+! prepare output fields
+if (mnproc==1) then
+   write (lp,'(a,i6,a)') ' ncwrt_onlysed: fields averaged over ',              &
+      &                    nacc_sed(iogrp),' steps'
+   write(lp,*) 'irec(iogrp)',irec(iogrp)
+endif
+rnacc=1./real(nacc_sed(iogrp))
+cmpflg=GLB_COMPFLAG(iogrp)
+
+! create output file
+if (GLB_NCFORMAT(iogrp)==1) then
+   call ncfopn(path1(1:path1_len)//fname(iogrp),'w','6', irec(iogrp),iotype)
+elseif (GLB_NCFORMAT(iogrp)==2) then
+   call ncfopn(path1(1:path1_len)//fname(iogrp),'w','h', irec(iogrp),iotype)
+else
+   call ncfopn(path1(1:path1_len)//fname(iogrp),'w','c', irec(iogrp),iotype)
+endif
+
+! define spatial and time dimensions
+if (cmpflg/=0) then
+   call ncdimc('pcomp',ip,0)
+else
+   call ncdims('x',itdm)
+   call ncdims('y',jtdm)
+endif
+!call ncdims('sigma',kdm)
+call ncdims('depth',ddm)
+call ncdims('ks',ks)
+call ncdims('bounds',2)
+call ncdims('time',0)
+call vardef_onlysed(iogrp,timeunits,calendar,cmpflg)
+call nctime(datenum,calendar,timeunits,startdate)
+
+! write auxillary dimension information
+!call ncwrt1('sigma','sigma',sigmar1)
+call ncwrt1('depth','depth',depthslev)
+call ncwrt1('depth_bnds','bounds depth',depthslev_bnds)
+
+! Store sediment fields
+call wrtsdm(jpowaic(iogrp), SDM_POWAIC(iogrp), rnacc*1e3,0., cmpflg,       &
+   &        'powdic','Porewater DIC',       ' ', 'mol C m-3')
+call wrtsdm(jpowaal(iogrp), SDM_POWAAL(iogrp), rnacc*1e3,0., cmpflg,       &
+   &        'powalk','Porewater alkalinity',' ', 'eq m-3')
+call wrtsdm(jpowaph(iogrp), SDM_POWAPH(iogrp), rnacc*1e3,0., cmpflg,       &
+   &        'powpho','Porewater phosphorus',' ', 'mol P m-3')
+call wrtsdm(jpowaox(iogrp), SDM_POWAOX(iogrp), rnacc*1e3,0., cmpflg,       &
+   &        'powox', 'Porewater oxygen',    ' ', 'mol O2 m-3')
+call wrtsdm(jpown2(iogrp),  SDM_POWN2(iogrp),  rnacc*1e3,0., cmpflg,       &
+   &        'pown2', 'Porewater N2',        ' ', 'mol N2 m-3')
+call wrtsdm(jpowno3(iogrp), SDM_POWNO3(iogrp), rnacc*1e3,0., cmpflg,       &
+   &        'powno3','Porewater nitrate',   ' ', 'mol N m-3')
+call wrtsdm(jpowasi(iogrp), SDM_POWASI(iogrp), rnacc*1e3,0., cmpflg,       &
+   &        'powsi', 'Porewater silicate',  ' ', 'mol Si m-3')
+call wrtsdm(jssso12(iogrp), SDM_SSSO12(iogrp), rnacc*1e3,0., cmpflg,       &
+   &        'ssso12','Sediment detritus',   ' ', 'mol P m-3')
+call wrtsdm(jssssil(iogrp), SDM_SSSSIL(iogrp), rnacc*1e3,0., cmpflg,       &
+   &        'ssssil','Sediment silicate',   ' ', 'mol Si m-3')
+call wrtsdm(jsssc12(iogrp), SDM_SSSC12(iogrp), rnacc*1e3,0., cmpflg,       &
+   &        'sssc12','Sediment CaCO3',      ' ', 'mol C m-3')
+call wrtsdm(jssster(iogrp), SDM_SSSTER(iogrp), rnacc*1e3,0., cmpflg,       &
+   &        'ssster','Sediment clay',       ' ', 'mol m-3')
+
+! Store sediment burial fields
+call wrtbur(jburssso12(iogrp),BUR_SSSO12(iogrp),rnacc*1e3,0.,cmpflg,       &
+   & 'buro12','Burial org carbon',' ','mol P m-2')
+call wrtbur(jbursssc12(iogrp),BUR_SSSC12(iogrp),rnacc*1e3,0.,cmpflg,       &
+   & 'burc12','Burial calcium ',' ','mol C m-2')
+call wrtbur(jburssssil(iogrp),BUR_SSSSIL(iogrp),rnacc*1e3,0.,cmpflg,       &
+   & 'bursil','Burial silicate',' ','mol Si m-2')
+call wrtbur(jburssster(iogrp),BUR_SSSTER(iogrp),rnacc*1e3,0.,cmpflg,       &
+   & 'burter','Burial clay',' ','mol  m-2')
+
+! close netcdf file
+call ncfcls
+
+call inisdm(jpowaic(iogrp),0.)
+call inisdm(jpowaal(iogrp),0.)
+call inisdm(jpowaph(iogrp),0.)
+call inisdm(jpowaox(iogrp),0.)
+call inisdm(jpown2(iogrp),0.)
+call inisdm(jpowno3(iogrp),0.)
+call inisdm(jpowasi(iogrp),0.)
+call inisdm(jssso12(iogrp),0.)
+call inisdm(jssssil(iogrp),0.)
+call inisdm(jsssc12(iogrp),0.)
+call inisdm(jssster(iogrp),0.)
+
+call inibur(jburssso12(iogrp),0.)
+call inibur(jbursssc12(iogrp),0.)
+call inibur(jburssssil(iogrp),0.)
+call inibur(jburssster(iogrp),0.)
+
+nacc_sed(iogrp)=0
+
+end subroutine ncwrt_onlysed
+
+subroutine vardef_onlysed(iogrp,timeunits,calendar,cmpflg)
+
+implicit none
+
+integer iogrp,cmpflg
+character timeunits*30,calendar*19
+
+call ncdefvar('time','time',ndouble,0)
+call ncattr('long_name','time')
+call ncattr('units',timeunits)
+call ncattr('calendar',calendar)
+
+call ncdefvar('depth','depth',ndouble,8)
+call ncattr('long_name','z level')
+call ncattr('units','m')
+call ncattr('positive','down')
+call ncattr('bounds','depth_bnds')
+
+call ncdefvar('depth_bnds','bounds depth',ndouble,8)
+call ncdefvar3d(CARFLX_BOT(iogrp),cmpflg,'p','carflx_bot',           &
+   &            'C flux to sediment',' ','mol C m-2 s-1',0)
+call ncdefvar3d(BSIFLX_BOT(iogrp),cmpflg,'p','bsiflx_bot',           &
+   &            'Opal flux to sediment',' ','mol Si m-2 s-1',0)
+call ncdefvar3d(CALFLX_BOT(iogrp),cmpflg,'p','calflx_bot',           &
+   &            'CaCO3 flux to sediment',' ','mol Ca m-2 s-1',0)
+
+! define sediment fields
+call ncdefvar3d(SDM_POWAIC(iogrp),cmpflg,'p','powdic','PoWa DIC',' ','mol C m-3',3)
+call ncdefvar3d(SDM_POWAAL(iogrp),cmpflg,'p','powalk','PoWa alkalinity',' ','eq m-3',3)
+call ncdefvar3d(SDM_POWAPH(iogrp),cmpflg,'p','powpho','PoWa phosphorus',' ','mol P m-3',3)
+call ncdefvar3d(SDM_POWAOX(iogrp),cmpflg,'p','powox','PoWa oxygen',' ','mol O2 m-3',3)
+call ncdefvar3d(SDM_POWN2(iogrp), cmpflg,'p','pown2','PoWa N2',' ','mol N2 m-3',3)
+call ncdefvar3d(SDM_POWNO3(iogrp),cmpflg,'p','powno3','PoWa nitrate',' ','mol N m-3',3)
+call ncdefvar3d(SDM_POWASI(iogrp),cmpflg,'p','powsi','PoWa silicate',' ','mol Si m-3',3)
+call ncdefvar3d(SDM_SSSO12(iogrp),cmpflg,'p','ssso12','Sediment detritus',' ','mol P m-3',3)
+call ncdefvar3d(SDM_SSSSIL(iogrp),cmpflg,'p','ssssil','Sediment silicate',' ','mol Si m-3',3)
+call ncdefvar3d(SDM_SSSC12(iogrp),cmpflg,'p','sssc12','Sediment CaCO3',' ','mol C m-3',3)
+call ncdefvar3d(SDM_SSSTER(iogrp),cmpflg,'p','ssster','Sediment clay',' ','mol m-3',3)
+
+! define sediment burial fields
+call ncdefvar3d(BUR_SSSO12(iogrp),cmpflg,'p','buro12','Burial org carbon',' ','mol P m-2',4)
+call ncdefvar3d(BUR_SSSC12(iogrp),cmpflg,'p','burc12','Burial calcium ',' ','mol C m-2',4)
+call ncdefvar3d(BUR_SSSSIL(iogrp),cmpflg,'p','bursil','Burial silicate',' ','mol Si m-2',4)
+call ncdefvar3d(BUR_SSSTER(iogrp),cmpflg,'p','burter','Burial clay',' ','mol  m-2',4)
+
+! enddef netcdf file
+call ncedef
+end subroutine vardef_onlysed
+
 end module mo_sedmnt_offline
 #endif
--- a/mod_dia.F	Thu Dec 06 14:29:09 2018 +0100
+++ b/mod_dia.F	Fri Dec 14 14:06:54 2018 +0100
@@ -198,7 +198,7 @@
 
 
       subroutine diafnm(runid,runid_len,expcnf,ctag,nstep,diagfq,diagmon
-     .  ,diagann,fname)
+     .  ,diagann,diagdec,diagcen,diagmil,fname)
 c
 c --- ------------------------------------------------------------------
 c --- Description: creates file name for the diagnostic output
@@ -212,6 +212,9 @@
 c ---   real diagfq    (in)    : diagnostic frequency     
 c ---   logi diagmon   (in)    : switch to show whether diagfq=month 
 c ---   logi diagann   (in)    : switch to show whether diagfq=year 
+c ---   logi diagdec   (in)    : switch to show whether diagfq=decade
+c ---   logi diagcen   (in)    : switch to show whether diagfq=century
+c ---   logi diagmil   (in)    : switch to show whether diagfq=millenium
 c ---   char fname     (out)   : file name 
 c --- ------------------------------------------
 c
@@ -222,7 +225,7 @@
       character runid*(*),expcnf*(*),ctag*(*),fname*(*)
       integer runid_len,nstep
       real diagfq
-      logical diagmon,diagann
+      logical diagmon,diagann,diagdec,diagcen,diagmil
 c
       real epsil
       parameter (epsil=1.e-11)
@@ -269,7 +272,7 @@
                write(fname,'(4a,i4.4,a,i2.2,a)')
      .           runid(1:runid_len),'.micom.',ctag,'.',ny,'-',nm,'.nc'
             endif
-          elseif (diagann) then
+          elseif (diagann .or. diagdec .or. diagcen .or. diagmil) then
             if ( lspinning_up_sed ) then
                write(fname,'(4a,i6.6,a)')
      .           runid(1:runid_len),'.micom.',ctag,'.',ny,'.nc'
@@ -288,6 +291,7 @@
                write(fname,'(4a,i4.4,a,i2.2,a,i2.2,a)')
      .           runid(1:runid_len),'.micom.',ctag,'.',ny,'-',nm,'-',nd,
      .           '.nc'
+            endif
           endif
         else
           if (mnproc.eq.1) then
@@ -1682,7 +1686,7 @@
       if (.not.append2file(iogrp)) then
         call diafnm(runid,runid_len,expcnf,trim(GLB_FNAMETAG(IOGRP)),
      .  nstep,filefq_phy(iogrp)/real(nstep_in_day),filemon_phy(iogrp),
-     .  fileann_phy(iogrp),fname(iogrp))
+     .  fileann_phy(iogrp),.false.,.false.,.false.,fname(iogrp))
         append2file(iogrp)=.true.
         irec(iogrp)=1
       else
--- a/ncout_hamocc.F	Thu Dec 06 14:29:09 2018 +0100
+++ b/ncout_hamocc.F	Fri Dec 14 14:06:54 2018 +0100
@@ -41,15 +41,20 @@
       if (.not.append2file(iogrp)) then
          call diafnm(runid,runid_len,expcnf,trim(GLB_FNAMETAG(iogrp)),nstep,
      .               filefq_bgc(iogrp)/real(nstep_in_day),filemon_bgc(iogrp),
-     .               fileann_bgc(iogrp),fname(iogrp))
+     .               fileann_bgc(iogrp),filedec_bgc(iogrp),filecen_bgc(iogrp),
+     .               filemil_bgc(iogrp),fname(iogrp))
          append2file(iogrp)=.true.
          irec(iogrp)=1
       else
          irec(iogrp)=irec(iogrp)+1
       endif
       if (((fileann_bgc(iogrp).and.nday_of_year.eq.1.or.
+     .  filedec_bgc(iogrp).and.mod(nyear,10).and.nday_of_year.eq.1.or.
+     .  filecen_bgc(iogrp).and.mod(nyear,100).and.nday_of_year.eq.1.or.
+     .  filemil_bgc(iogrp).and.mod(nyear,1000).and.nday_of_year.eq.1.or.
      .  filemon_bgc(iogrp).and.nday.eq.1).and.mod(nstep,nstep_in_day)
-     .  .le.1).or..not.(fileann_bgc(iogrp).or.filemon_bgc(iogrp)).and.
+     .  .le.1).or..not.(filemil_bgc(iogrp).or.filecen_bgc(iogrp).or.
+     .  filedec_bgc(iogrp).or.fileann_bgc(iogrp).or.filemon_bgc(iogrp)).and.
      .  mod(nstep+.5,filefq_bgc(iogrp)).lt.2.) then
         append2file(iogrp)=.false.
       endif
--- a/ncout_onlysed.F90	Thu Dec 06 14:29:09 2018 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,175 +0,0 @@
-#if defined(SED_OFFLINE)
-subroutine ncwrt_onlysed(iogrp)
-
-!-----------------------------------------------------------------------
-!
-! Write diagnostic sediment fields during off-line sediment spin-up
-!
-! Copyright (C) 2018 Marco van Hulten <Marco.Hulten@uib.no> et al.
-!                    Geophysical Institute @ University of Bergen
-!
-! This subroutine is free software: you can redistribute it and/or modify
-! it under the terms of the GNU General Public License as published by
-! the Free Software Foundation, either version 3 of the License, or
-! (at your option) any later version.
-!
-!
-! Method:
-!
-! History:
-!
-! This is based on ncout_hamocc.F that contains ncwrt_bgc().
-!
-! Marco van Hulten <Marco.Hulten@uib.no>           2018-11-20
-! - Initial code based on ncout_hamocc.F.
-!
-!-----------------------------------------------------------------------
-
-use mod_xc
-use mod_dia       , only: diafnm,sigmar1,iotype
-use mo_control_bgc, only: dtbgc, nburst
-use mo_biomod     , only: k0100,k0500,k1000,k2000,k4000
-!use mo_bgcmean   ! FIXME: otherwise collision of GLB_FNAMETAG, ...?
-use mo_sedmnt_offline
-
-implicit none
-
-#include "common_clndr.h90"
-#include "common_blocks.h90"
-
-integer, intent(in) :: iogrp
-
-integer i,j,k,l,nt
-integer nhour,ny,nm,nd,dayfrac,irec(nsedmax),cmpflg
-character(len= 2) seqstring
-character(len=80) fname(nsedmax)
-character(len=20) startdate
-character(len=30) timeunits
-real datenum,rnacc
-logical append2file(nsedmax)
-data append2file /nsedmax*.false./
-save fname,irec,append2file
-
-! set time information
-timeunits=' '
-startdate=' '
-fname=' '
-write(timeunits,'(a11,i4.4,a1,i2.2,a1,i2.2,a6)')                           &
-   &            'days since ',min(1800,nyear0),'-',1,'-',1,' 00:00'
-write(startdate,'(i4.4,a1,i2.2,a1,i2.2,a6)')                               &
-   &            nyear0,'-',nmonth0,'-',nday0,' 00:00'
-datenum=time-time0-0.5*diagfq_sed(iogrp)/nstep_in_day
-
-! get file name
-write (seqstring,'(I0.2)') nburst
-if (.not.append2file(iogrp)) then
-   call diafnm(runid,runid_len,expcnf,trim(GLB_FNAMETAG(iogrp))//"."//seqstring,nstep, &
-      &         filefq_sed(iogrp)/real(nstep_in_day),filemon_sed(iogrp),       &
-      &         fileann_sed(iogrp),fname(iogrp))
-   append2file(iogrp)=.true.
-   irec(iogrp)=1
-else
-   irec(iogrp)=irec(iogrp)+1
-endif
-if ( (fileann_sed(iogrp) .or. filemon_sed(iogrp))                       &
-   & .or. .not.(fileann_sed(iogrp) .or. filemon_sed(iogrp)) .and.       &
-   &  mod(nstep+.5,filefq_sed(iogrp))<2.) then
-   append2file(iogrp) = .false.
-endif
-
-! prepare output fields
-if (mnproc==1) then
-   write (lp,'(a,i6,a)') ' ncwrt_sed: fields averaged over ',              &
-      &                    nacc_sed(iogrp),' steps'
-   write(lp,*) 'irec(iogrp)',irec(iogrp)
-endif
-rnacc=1./real(nacc_sed(iogrp))
-cmpflg=GLB_COMPFLAG(iogrp)
-
-! create output file
-if (GLB_NCFORMAT(iogrp)==1) then
-   call ncfopn(path1(1:path1_len)//fname(iogrp),'w','6', irec(iogrp),iotype)
-elseif (GLB_NCFORMAT(iogrp)==2) then
-   call ncfopn(path1(1:path1_len)//fname(iogrp),'w','h', irec(iogrp),iotype)
-else
-   call ncfopn(path1(1:path1_len)//fname(iogrp),'w','c', irec(iogrp),iotype)
-endif
-
-! define spatial and time dimensions
-if (cmpflg/=0) then
-   call ncdimc('pcomp',ip,0)
-else
-   call ncdims('x',itdm)
-   call ncdims('y',jtdm)
-endif
-call ncdims('sigma',kdm)
-call ncdims('depth',ddm)
-call ncdims('ks',ks)
-call ncdims('bounds',2)
-call ncdims('time',0)
-call hamoccvardef(iogrp,timeunits,calendar,cmpflg)
-call nctime(datenum,calendar,timeunits,startdate)
-
-! write auxillary dimension information
-call ncwrt1('sigma','sigma',sigmar1)
-call ncwrt1('depth','depth',depthslev)
-call ncwrt1('depth_bnds','bounds depth',depthslev_bnds)
-
-! Store sediment fields
-call wrtsdm(jpowaic(iogrp), SDM_POWAIC(iogrp), rnacc*1e3,0., cmpflg,       &
-   &        'powdic','Porewater DIC',       ' ', 'mol C m-3')
-call wrtsdm(jpowaal(iogrp), SDM_POWAAL(iogrp), rnacc*1e3,0., cmpflg,       &
-   &        'powalk','Porewater alkalinity',' ', 'eq m-3')
-call wrtsdm(jpowaph(iogrp), SDM_POWAPH(iogrp), rnacc*1e3,0., cmpflg,       &
-   &        'powpho','Porewater phosphorus',' ', 'mol P m-3')
-call wrtsdm(jpowaox(iogrp), SDM_POWAOX(iogrp), rnacc*1e3,0., cmpflg,       &
-   &        'powox', 'Porewater oxygen',    ' ', 'mol O2 m-3')
-call wrtsdm(jpown2(iogrp),  SDM_POWN2(iogrp),  rnacc*1e3,0., cmpflg,       &
-   &        'pown2', 'Porewater N2',        ' ', 'mol N2 m-3')
-call wrtsdm(jpowno3(iogrp), SDM_POWNO3(iogrp), rnacc*1e3,0., cmpflg,       &
-   &        'powno3','Porewater nitrate',   ' ', 'mol N m-3')
-call wrtsdm(jpowasi(iogrp), SDM_POWASI(iogrp), rnacc*1e3,0., cmpflg,       &
-   &        'powsi', 'Porewater silicate',  ' ', 'mol Si m-3')
-call wrtsdm(jssso12(iogrp), SDM_SSSO12(iogrp), rnacc*1e3,0., cmpflg,       &
-   &        'ssso12','Sediment detritus',   ' ', 'mol P m-3')
-call wrtsdm(jssssil(iogrp), SDM_SSSSIL(iogrp), rnacc*1e3,0., cmpflg,       &
-   &        'ssssil','Sediment silicate',   ' ', 'mol Si m-3')
-call wrtsdm(jsssc12(iogrp), SDM_SSSC12(iogrp), rnacc*1e3,0., cmpflg,       &
-   &        'sssc12','Sediment CaCO3',      ' ', 'mol C m-3')
-call wrtsdm(jssster(iogrp), SDM_SSSTER(iogrp), rnacc*1e3,0., cmpflg,       &
-   &        'ssster','Sediment clay',       ' ', 'mol m-3')
-
-! Store sediment burial fields
-call wrtbur(jburssso12(iogrp),BUR_SSSO12(iogrp),rnacc*1e3,0.,cmpflg,       &
-   & 'buro12','Burial org carbon',' ','mol P m-2')
-call wrtbur(jbursssc12(iogrp),BUR_SSSC12(iogrp),rnacc*1e3,0.,cmpflg,       &
-   & 'burc12','Burial calcium ',' ','mol C m-2')
-call wrtbur(jburssssil(iogrp),BUR_SSSSIL(iogrp),rnacc*1e3,0.,cmpflg,       &
-   & 'bursil','Burial silicate',' ','mol Si m-2')
-call wrtbur(jburssster(iogrp),BUR_SSSTER(iogrp),rnacc*1e3,0.,cmpflg,       &
-   & 'burter','Burial clay',' ','mol  m-2')
-
-! close netcdf file
-call ncfcls
-
-call inisdm(jpowaic(iogrp),0.)
-call inisdm(jpowaal(iogrp),0.)
-call inisdm(jpowaph(iogrp),0.)
-call inisdm(jpowaox(iogrp),0.)
-call inisdm(jpown2(iogrp),0.)
-call inisdm(jpowno3(iogrp),0.)
-call inisdm(jpowasi(iogrp),0.)
-call inisdm(jssso12(iogrp),0.)
-call inisdm(jssssil(iogrp),0.)
-call inisdm(jsssc12(iogrp),0.)
-call inisdm(jssster(iogrp),0.)
-
-call inibur(jburssso12(iogrp),0.)
-call inibur(jbursssc12(iogrp),0.)
-call inibur(jburssssil(iogrp),0.)
-call inibur(jburssster(iogrp),0.)
-
-nacc_sed(iogrp)=0
-
-end subroutine ncwrt_onlysed
-#endif