fabm: added river and hotstart support - Hofmeister; changed dimension order...

fabm: added river and hotstart support - Hofmeister; changed dimension order to start with i,j,k; use missing_value,minimum,maximum metadata

src/3d/getm_fabm.F90

@@ -18,6 +18,7 @@
 #else
    use domain, only: dx,dy,ard1
 #endif
+   use variables_3d,only: fabm_pel,fabm_ben,fabm_diag,fabm_diag_hz
    use variables_3d, only: uu,vv,ww,hun,hvn,ho,hn
    use variables_3d, only: nuh,T,S,rho,a,g1,g2,taub
    use advection_3d, only: do_advection_3d
@@ -51,8 +52,6 @@
    REALTYPE, dimension(:,:,:), allocatable :: ff
 #endif
 
-   REALTYPE, allocatable, dimension(:,:,:,:) :: cc_pel,cc_diag
-   REALTYPE, allocatable, dimension(:,:,:)   :: cc_ben,cc_diag_hz
 !
 ! !REVISION HISTORY:
 !  Original author(s): Hans Burchard & Karsten Bolding
@@ -97,38 +96,38 @@
    LEVEL2 'init_getm_fabm()'
 
 !  Initialize FABM.
-   call init_gotm_fabm(kmax,NAMLST2,'fabm.nml')
+   call init_gotm_fabm(kmax,NAMLST2,'gotm_fabm.nml')
 
    if (fabm_calc) then
 !     Temporary: make sure diagnostic variables store the last value,
 !     not their time integral. This will be redundant when time-integrating/averaging
 !     is moved from FABM to the physical host.
-      do n=1,ubound(model%info%diagnostic_variables,1)
+      do n=1,size(model%info%diagnostic_variables)
          model%info%diagnostic_variables(n)%time_treatment = time_treatment_last
       end do
-      do n=1,ubound(model%info%diagnostic_variables_hz,1)
+      do n=1,size(model%info%diagnostic_variables_hz)
          model%info%diagnostic_variables_hz(n)%time_treatment = time_treatment_last
       end do
 
 !     Allocate memory for pelagic state variables.
-      allocate(cc_pel(ubound(model%info%state_variables,1),I3DFIELD),stat=rc)
+      allocate(fabm_pel(I3DFIELD,size(model%info%state_variables)),stat=rc)
       if (rc /= 0) stop 'init_getm_fabm: Error allocating memory (cc_pel)'
-      cc_pel = _ZERO_
+      fabm_pel = _ZERO_
 
 !     Allocate memory for benthic state variables.
-      allocate(cc_ben(ubound(model%info%state_variables_ben,1),I2DFIELD),stat=rc)
+      allocate(fabm_ben(I2DFIELD,size(model%info%state_variables_ben)),stat=rc)
       if (rc /= 0) stop 'init_getm_fabm: Error allocating memory (cc_ben)'
-      cc_ben = _ZERO_
+      fabm_ben = _ZERO_
 
 !     Allocate memory for 3D diagnostic variables.
-      allocate(cc_diag(ubound(model%info%diagnostic_variables,1),I3DFIELD),stat=rc)
+      allocate(fabm_diag(I3DFIELD,size(model%info%diagnostic_variables)),stat=rc)
       if (rc /= 0) stop 'init_getm_fabm: Error allocating memory (cc_diag)'
-      cc_diag = _ZERO_
+      fabm_diag = _ZERO_
 
 !     Allocate memory for 2D [horizontal-only] diagnostic variables.
-      allocate(cc_diag_hz(ubound(model%info%diagnostic_variables_hz,1),I2DFIELD),stat=rc)
+      allocate(fabm_diag_hz(I2DFIELD,size(model%info%diagnostic_variables_hz)),stat=rc)
       if (rc /= 0) stop 'init_getm_fabm: Error allocating memory (cc_diag_hz)'
-      cc_diag_hz = _ZERO_
+      fabm_diag_hz = _ZERO_
 
 !     Read settings specific to GETM-FABM interaction.
       open(NAMLST2,status='unknown',file=trim(nml_file))
@@ -151,17 +150,21 @@
             do j=jmin,jmax
                do i=imin,imax
                   if (az(i,j) .ge. 1 ) then
-                     cc_pel(:,i,j,:) = cc_col(1:ubound(model%info%state_variables,1) ,:)
-                     cc_ben(:,i,j)   = cc_col(ubound(model%info%state_variables,1)+1:,1)
+                     do n=1,size(model%info%state_variables)
+                        fabm_pel(i,j,:,n) = cc_col(n,:)
+                     end do
+                     do n=1,size(model%info%state_variables_ben)
+                        fabm_ben(i,j,  n) = cc_col(size(model%info%state_variables)+n,1)
+                     end do
                   end if
                end do
             end do
          case(2)
             LEVEL3 'reading initial biogeochemical fields from ',trim(fabm_init_file)
-            do n=1,ubound(model%info%state_variables,1)
+            do n=1,size(model%info%state_variables)
                LEVEL4 'inquiring ',trim(model%info%state_variables(n)%name)
                call get_field(fabm_init_file,trim(model%info%state_variables(n)%name),fabm_field_no, &
-                              cc_pel(n,:,:,:))
+                              fabm_pel(:,:,:,n))
             end do
          case default
             FATAL 'Not valid fabm_init_method specified'
@@ -169,8 +172,8 @@
       end select
 
 !     Update halos with biogeochemical variable values (distribute initial values).
-      do n=1,ubound(model%info%state_variables,1)
-         call update_3d_halo(cc_pel(n,:,:,:),cc_pel(n,:,:,:),az, &
+      do n=1,size(model%info%state_variables)
+         call update_3d_halo(fabm_pel(:,:,:,n),fabm_pel(:,:,:,n),az, &
                              imin,jmin,imax,jmax,kmax,D_TAG)
          call wait_halo(D_TAG)
       end do
@@ -224,7 +227,7 @@
 ! !DESCRIPTION:
 !
 ! !USES:
-   use getm_timers, only: tic, toc, TIM_GETM_BIO
+   use getm_timers, only: tic, toc, TIM_GETM_FABM, TIM_ADVECTFABM
    IMPLICIT NONE
 !
 ! !INPUT PARAMETERS:
@@ -243,10 +246,14 @@
 !-----------------------------------------------------------------------
 !BOC
 
-   call tic(TIM_GETM_BIO)
+   call tic(TIM_GETM_FABM)
 
 !  First we do all the vertical processes
+#ifdef SLICE_MODEL
+   do j=2,2
+#else
    do j=jmin,jmax
+#endif
       do i=imin,imax
          if (az(i,j) .ge. 1 ) then
 
@@ -274,10 +281,18 @@
             taub_nonnorm = taub(i,j)*rho_0
 
 !           Copy current values of biogeochemical variables from full 3D field to columns.
-            cc_col(1:ubound(model%info%state_variables,1) ,:) = cc_pel(:,i,j,:)
-            cc_col(ubound(model%info%state_variables,1)+1:,1) = cc_ben(:,i,j)
-            cc_diag_col    = cc_diag(:,i,j,:)
-            cc_diag_hz_col = cc_diag_hz(:,i,j)
+            do n=1,size(model%info%state_variables)
+               cc_col(n,:) = fabm_pel(i,j,:,n)
+            end do
+            do n=1,size(model%info%state_variables_ben)
+               cc_col(size(model%info%state_variables)+n,1) = fabm_ben(i,j,n)
+            end do
+            do n=1,size(model%info%diagnostic_variables)
+               cc_diag_col(n,:) = fabm_diag(i,j,:,n)
+            end do
+            do n=1,size(model%info%diagnostic_variables_hz)
+               cc_diag_hz_col(n) = fabm_diag_hz(i,j,n)
+            end do
 
 !           Transfer pointers to physical environment variables to FABM.
             call set_env_gotm_fabm(dt,0,0,T(i,j,1:),S(i,j,1:), &
@@ -289,20 +304,38 @@
             call do_gotm_fabm(kmax)
 
 !           Copy updated column values of biogeochemical variables to full 3D field.
-            cc_pel    (:,i,j,:) = cc_col(1:ubound(model%info%state_variables,1) ,:)
-            cc_ben    (:,i,j)   = cc_col(ubound(model%info%state_variables,1)+1:,1)
-            cc_diag   (:,i,j,:) = cc_diag_col
-            cc_diag_hz(:,i,j)   = cc_diag_hz_col
+            do n=1,size(model%info%state_variables)
+               fabm_pel(i,j,:,n) = cc_col(n,:)
+            end do
+            do n=1,size(model%info%state_variables_ben)
+               fabm_ben(i,j,n) = cc_col(size(model%info%state_variables)+n,1)
+            end do
+            do n=1,size(model%info%diagnostic_variables)
+               fabm_diag(i,j,:,n) = cc_diag_col(n,:)
+            end do
+            do n=1,size(model%info%diagnostic_variables_hz)
+               fabm_diag_hz(i,j,n) = cc_diag_hz_col(n)
+            end do
 
          end if
       end do
    end do
 
+#ifdef SLICE_MODEL
+      do i=imin,imax
+         cc_pel(:,i,3,:)=cc_pel(:,i,2,:)
+         cc_ben(:,i,3)  =cc_ben(:,i,2)
+      end do
+#endif
+
+
+
 !  Advect pelagic biogeochemical variables.
-   do n=1,ubound(model%info%state_variables,1)
+   call tic(TIM_ADVECTFABM)
+   do n=1,size(model%info%state_variables)
 
 #if 1
-      ff = cc_pel(n,:,:,:)
+      ff = fabm_pel(:,:,:,n)
 
       call update_3d_halo(ff,ff,az, &
                           imin,jmin,imax,jmax,kmax,D_TAG)
@@ -312,19 +345,20 @@
               delxu,delxv,delyu,delyv,area_inv,az,au,av, &
               fabm_hor_adv,fabm_ver_adv,fabm_adv_split,fabm_AH)
 
-      cc_pel(n,:,:,:) = ff
+      fabm_pel(:,:,:,n) = ff
 #else
-      call update_3d_halo(cc3d(n,:,:,:),cc3d(n,:,:,:),az, &
+      call update_3d_halo(fabm_pel(:,:,:,n),fabm_pel(:,:,:,n),az, &
                           imin,jmin,imax,jmax,kmax,D_TAG)
       call wait_halo(D_TAG)
 
-      call do_advection_3d(dt,cc3d(n,:,:,:),uu,vv,ww,hun,hvn,ho,hn, &
+      call do_advection_3d(dt,fabm_pel(:,:,:,n),uu,vv,ww,hun,hvn,ho,hn, &
               delxu,delxv,delyu,delyv,area_inv,az,au,av, &
               fabm_hor_adv,fabm_ver_adv,fabm_adv_split,fabm_AH)
 #endif
    end do
+   call toc(TIM_ADVECTFABM)
 
-   call toc(TIM_GETM_BIO)
+   call toc(TIM_GETM_FABM)
 
    return
    end subroutine do_getm_fabm

src/3d/rivers.F90

@@ -42,6 +42,10 @@
    use bio, only: bio_calc
    use bio_var, only: numc
    use variables_3d, only: cc3d
+#endif
+#ifdef _FABM_
+   use gotm_fabm, only: model
+   use getm_fabm, only: fabm_pel
 #endif
    IMPLICIT NONE
 !
@@ -51,6 +55,9 @@
    public init_rivers, do_rivers, clean_rivers
 #ifdef GETM_BIO
    public init_rivers_bio
+#endif
+#ifdef _FABM_
+   public init_rivers_fabm
 #endif
    integer, public                     :: river_method=0,nriver=0,rriver=0
    logical,public                      :: use_river_temp = .false.
@@ -71,6 +78,9 @@
    REALTYPE, public, allocatable       :: river_bio(:,:)
    REALTYPE, public, parameter         :: bio_missing=-9999.0
 #endif
+#ifdef _FABM_
+   REALTYPE, public, allocatable       :: river_fabm(:,:)
+#endif
 !
 ! !PRIVATE DATA MEMBERS:
    integer                   :: river_format=2
@@ -327,6 +337,55 @@
 !EOC
 #endif
 
+#ifdef _FABM_
+!-----------------------------------------------------------------------
+!BOP
+!
+! !IROUTINE: init_rivers_fabm
+!
+! !INTERFACE:
+   subroutine init_rivers_fabm()
+!
+! !DESCRIPTION:
+! First, memory for storing the biological loads from rivers is
+! allocated.
+! The variable - {\tt river\_fabm} - is initialised to  - variable-
+! specific missing values obtained provided by FABM.
+!
+! !USES:
+   IMPLICIT NONE
+!
+! !LOCAL VARIABLES:
+   integer                   :: rc,m
+!EOP
+!-------------------------------------------------------------------------
+!BOC
+#ifdef DEBUG
+   integer, save :: Ncall = 0
+   Ncall = Ncall+1
+   write(debug,*) 'init_rivers_fabm() # ',Ncall
+#endif
+
+   if (allocated(fabm_pel)) then
+      LEVEL1 'init_rivers_fabm()'
+
+      allocate(river_fabm(nriver,size(model%info%state_variables)),stat=rc)
+      if (rc /= 0) stop 'rivers: Error allocating memory (river_fabm)'
+
+      do m=1,size(model%info%state_variables)
+         river_fabm(:,m) = model%info%state_variables(m)%missing_value
+      end do
+   end if
+
+#ifdef DEBUG
+   write(debug,*) 'Leaving init_rivers_fabm()'
+   write(debug,*)
+#endif
+   return
+   end subroutine init_rivers_fabm
+!EOC
+#endif
+
 !-----------------------------------------------------------------------
 !BOP
 !
@@ -415,6 +474,18 @@
                      end do
                   end if
 #endif
+#ifdef _FABM_
+                  if (allocated(fabm_pel)) then
+                     do m=1,size(model%info%state_variables)
+                        if ( river_fabm(n,m) .ne. model%info%state_variables(m)%missing_value ) then
+                           fabm_pel(i,j,1:kmax,m) = &
+                                 (fabm_pel(i,j,1:kmax,m)*(H(i,j)+ssen(i,j)) &
+                                 + river_fabm(n,m)*macro_height(n))      &
+                                 / (H(i,j)+ssen(i,j)+macro_height(n))
+                        end if
+                     end do
+                  end if
+#endif
 !                 Changes of total and layer height due to river inflow:
                   hn(i,j,1:kmax) = hn(i,j,1:kmax)/(H(i,j)+ssen(i,j)) &
                                   *(H(i,j)+ssen(i,j)+macro_height(n))
@@ -475,6 +546,9 @@
                tot = tot+irr(n)
             end if
          end do
+#ifdef _FABM_
+         if (allocated(river_fabm)) deallocate(river_fabm)
+#endif
       case default
          FATAL 'Not valid rivers_method specified'
          stop 'clean_rivers'

src/3d/variables_3d.F90

@@ -129,6 +129,10 @@
 #ifdef GETM_BIO
    REALTYPE, allocatable               :: cc3d(:,:,:,:)
    REALTYPE, allocatable               :: ws3d(:,:,:,:)
+#endif
+#ifdef _FABM_
+   REALTYPE, allocatable, dimension(:,:,:,:) :: fabm_pel,fabm_diag
+   REALTYPE, allocatable, dimension(:,:,:)   :: fabm_ben,fabm_diag_hz
 #endif
    integer                             :: size3d_field
    integer                             :: mem3d

src/futils/getm_timers.F90

@@ -67,6 +67,7 @@
    integer, parameter :: TIM_CALCMEANF   = 68   ! 3d calc_mean_fields
    integer, parameter :: TIM_METEO       = 70   ! do_meteo (could use + halo)
    integer, parameter :: TIM_GETM_BIO    = 72   ! do_getm_bio
+   integer, parameter :: TIM_GETM_FABM   = 73   ! do_getm_fabm
    ! These catch compuations in integrate_[23]d, which are not in other timers:
    integer, parameter :: TIM_INTEGR2D    = 80   ! 2d integrate_2d - remaining stuff
    integer, parameter :: TIM_INTEGR3D    = 81   ! 3d integrate_3d - remaining stuff
@@ -81,6 +82,9 @@
    integer, parameter :: TIM_HALO2D      = 110  ! do halo 2d (initialize comm)
    integer, parameter :: TIM_HALO3D      = 111  ! do halo 3d (initialize comm)
    integer, parameter :: TIM_HALOWAIT    = 112  ! wait_halo (2d+3d both)
+   integer, parameter :: TIM_ADVECTBIO   = 113  ! advection_3d bio
+   integer, parameter :: TIM_ADVECTFABM  = 114  ! advection_3d fabm
+
    ! This is test timers for temporary coding purposes:
    !  Note: All timers with index 170+ (test_timer_first) are
    !  considered test timers, so dont implement your timers here
@@ -223,6 +227,15 @@
    timernames(TIM_HALO3D)      = ' sum do_halo_3d'
 #endif
 
+#ifdef GETM_BIO
+   timernames(TIM_GETM_BIO)    = 'getm_bio'
+   timernames(TIM_ADVECTBIO)   = ' advection getm_bio'
+#endif
+#ifdef _FABM_
+   timernames(TIM_GETM_FABM)   = 'getm_fabm'
+   timernames(TIM_ADVECTFABM)  = ' advection getm_fabm'
+#endif
+
 #ifdef STRUCTURE_FRICTION
    timernames(TIM_STRCTFRICT)  = 'structure_friction_3d'
 #endif

src/getm/initialise.F90

@@ -62,6 +62,7 @@
 #endif
 #ifdef _FABM_
    use getm_fabm, only: init_getm_fabm
+   use rivers, only: init_rivers_fabm 
 #endif
 #ifdef GETM_BIO
    use bio, only: bio_calc
@@ -236,7 +237,7 @@
 #endif
 #ifdef _FABM_
       call init_getm_fabm(trim(input_dir) // 'getm_fabm.inp')
-!KB      call init_rivers_bio
+      call init_rivers_fabm
 #endif
 #ifdef GETM_BIO
       call init_getm_bio(trim(input_dir) // 'getm_bio.inp')

src/ncdf/create_restart_ncdf.F90

@@ -25,6 +25,9 @@
 #ifdef GETM_BIO
    use bio, only: bio_calc
    use bio_var, only: numc
+#endif
+#ifdef _FABM_
+   use gotm_fabm, only: fabm_calc,model
 #endif
    IMPLICIT NONE
 !
@@ -72,6 +75,19 @@
       if (status .NE. NF90_NOERR) go to 10
    end if
 #endif
+#ifdef _FABM_
+   if (fabm_calc) then
+      status = nf90_def_dim(ncid, "fabm_pel_dim", size(model%info%state_variables), fabmpeldim_id)
+      if (status .NE. NF90_NOERR) go to 10
+
+      if (size(model%info%state_variables_ben).gt.0) then
+         status = nf90_def_dim(ncid, "fabm_ben_dim", size(model%info%state_variables_ben), fabmbendim_id)
+         if (status .NE. NF90_NOERR) go to 10
+      else
+         fabmbendim_id = 0
+      end if
+   end if
+#endif
 
    status = nf90_def_var(ncid, "loop", nf90_int, loop_id)
    if (status .NE. NF90_NOERR) go to 10
@@ -233,6 +249,21 @@
                                 bio_id)
          if (status .NE. NF90_NOERR) go to 10
       end if
+#endif
+#ifdef _FABM_
+      if (fabm_calc) then
+         status = nf90_def_var(ncid, "fabm_pel", nf90_double, &
+                               (/ fabmpeldim_id, xdim_id, ydim_id, zdim_id /), &
+                                fabm_pel_id)
+         if (status .NE. NF90_NOERR) go to 10
+
+         if (fabmbendim_id.gt.0) then
+            status = nf90_def_var(ncid, "fabm_ben", nf90_double, &
+                               (/ fabmbendim_id, xdim_id, ydim_id /), &
+                                fabm_ben_id)
+            if (status .NE. NF90_NOERR) go to 10
+         endif
+      end if
 #endif
    end if
 #endif

src/ncdf/init_3d_ncdf.F90

@@ -26,7 +26,7 @@
    use bio_var, only: numc,var_names,var_units,var_long
 #endif
 #ifdef _FABM_
-   use gotm_fabm, only: model
+   use gotm_fabm, only: model,fabm_calc
 #endif
 
    IMPLICIT NONE
@@ -393,54 +393,63 @@
 #endif
 
 #ifdef _FABM_
-   fv = bio_missing
-   mv = bio_missing
-   vr(1) = _ZERO_
-   vr(2) = 9999.
-
-   allocate(fabm_ids(ubound(model%info%state_variables,1)),stat=rc)
-   if (rc /= 0) stop 'init_3d_ncdf(): Error allocating memory (fabm_ids)'
-   do n=1,ubound(model%info%state_variables,1)
-      err = nf90_def_var(ncid,model%info%state_variables(n)%name,NCDF_FLOAT_PRECISION,f4_dims,fabm_ids(n))
-      if (err .NE.  NF90_NOERR) go to 10
-      call set_attributes(ncid,fabm_ids(n), &
-                          long_name=trim(model%info%state_variables(n)%longname), &
-                          units=trim(model%info%state_variables(n)%units), &
-                          FillValue=fv,missing_value=mv,valid_range=vr)
-   end do
-
-   allocate(fabm_ids_ben(ubound(model%info%state_variables_ben,1)),stat=rc)
-   if (rc /= 0) stop 'init_3d_ncdf(): Error allocating memory (fabm_ids_ben)'
-   do n=1,ubound(model%info%state_variables_ben,1)
-      err = nf90_def_var(ncid,model%info%state_variables_ben(n)%name,NCDF_FLOAT_PRECISION,f4_dims,fabm_ids_ben(n))
-      if (err .NE.  NF90_NOERR) go to 10
-      call set_attributes(ncid,fabm_ids_ben(n), &
-                          long_name=trim(model%info%state_variables_ben(n)%longname), &
-                          units=trim(model%info%state_variables_ben(n)%units), &
-                          FillValue=fv,missing_value=mv,valid_range=vr)
-   end do
-
-   allocate(fabm_ids_diag(ubound(model%info%diagnostic_variables,1)),stat=rc)
-   if (rc /= 0) stop 'init_3d_ncdf(): Error allocating memory (fabm_ids_diag)'
-   do n=1,ubound(model%info%diagnostic_variables,1)
-      err = nf90_def_var(ncid,model%info%diagnostic_variables(n)%name,NCDF_FLOAT_PRECISION,f4_dims,fabm_ids_diag(n))
-      if (err .NE.  NF90_NOERR) go to 10
-      call set_attributes(ncid,fabm_ids_diag(n), &
-                          long_name=trim(model%info%diagnostic_variables(n)%longname), &
-                          units=trim(model%info%diagnostic_variables(n)%units), &
-                          FillValue=fv,missing_value=mv,valid_range=vr)
-   end do
-
-   allocate(fabm_ids_diag_hz(ubound(model%info%diagnostic_variables_hz,1)),stat=rc)
-   if (rc /= 0) stop 'init_3d_ncdf(): Error allocating memory (fabm_ids_diag_hz)'
-   do n=1,ubound(model%info%diagnostic_variables_hz,1)
-      err = nf90_def_var(ncid,model%info%diagnostic_variables_hz(n)%name,NCDF_FLOAT_PRECISION,f3_dims,fabm_ids_diag_hz(n))
-      if (err .NE.  NF90_NOERR) go to 10
-      call set_attributes(ncid,fabm_ids_diag_hz(n), &
-                          long_name=trim(model%info%diagnostic_variables_hz(n)%longname), &
-                          units=trim(model%info%diagnostic_variables_hz(n)%units), &
-                          FillValue=fv,missing_value=mv,valid_range=vr)
-   end do
+   if (fabm_calc) then
+      allocate(fabm_ids(size(model%info%state_variables)),stat=rc)
+      if (rc /= 0) stop 'init_3d_ncdf(): Error allocating memory (fabm_ids)'
+      do n=1,size(model%info%state_variables)
+         err = nf90_def_var(ncid,model%info%state_variables(n)%name,NCDF_FLOAT_PRECISION,f4_dims,fabm_ids(n))
+         if (err .NE.  NF90_NOERR) go to 10
+         call set_attributes(ncid,fabm_ids(n), &
+                          long_name    =trim(model%info%state_variables(n)%longname), &
+                          units        =trim(model%info%state_variables(n)%units),    &
+                          FillValue    =model%info%state_variables(n)%missing_value,  &
+                          missing_value=model%info%state_variables(n)%missing_value,  &
+                          valid_min    =model%info%state_variables(n)%minimum,        &
+                          valid_max    =model%info%state_variables(n)%maximum)
+      end do
+
+      allocate(fabm_ids_ben(size(model%info%state_variables_ben)),stat=rc)
+      if (rc /= 0) stop 'init_3d_ncdf(): Error allocating memory (fabm_ids_ben)'
+      do n=1,size(model%info%state_variables_ben)
+         err = nf90_def_var(ncid,model%info%state_variables_ben(n)%name,NCDF_FLOAT_PRECISION,f3_dims,fabm_ids_ben(n))
+         if (err .NE.  NF90_NOERR) go to 10
+         call set_attributes(ncid,fabm_ids_ben(n), &
+                          long_name    =trim(model%info%state_variables_ben(n)%longname), &
+                          units        =trim(model%info%state_variables_ben(n)%units),    &
+                          FillValue    =model%info%state_variables_ben(n)%missing_value,  &
+                          missing_value=model%info%state_variables_ben(n)%missing_value,  &
+                          valid_min    =model%info%state_variables_ben(n)%minimum,        &
+                          valid_max    =model%info%state_variables_ben(n)%maximum)
+      end do
+
+      allocate(fabm_ids_diag(size(model%info%diagnostic_variables)),stat=rc)
+      if (rc /= 0) stop 'init_3d_ncdf(): Error allocating memory (fabm_ids_diag)'
+      do n=1,size(model%info%diagnostic_variables)
+         err = nf90_def_var(ncid,model%info%diagnostic_variables(n)%name,NCDF_FLOAT_PRECISION,f4_dims,fabm_ids_diag(n))
+         if (err .NE.  NF90_NOERR) go to 10
+         call set_attributes(ncid,fabm_ids_diag(n), &
+                          long_name    =trim(model%info%diagnostic_variables(n)%longname), &
+                          units        =trim(model%info%diagnostic_variables(n)%units),    &
+                          FillValue    =model%info%diagnostic_variables(n)%missing_value,  &
+                          missing_value=model%info%diagnostic_variables(n)%missing_value,  &
+                          valid_min    =model%info%diagnostic_variables(n)%minimum,        &
+                          valid_max    =model%info%diagnostic_variables(n)%maximum)
+      end do
+
+      allocate(fabm_ids_diag_hz(size(model%info%diagnostic_variables_hz)),stat=rc)
+      if (rc /= 0) stop 'init_3d_ncdf(): Error allocating memory (fabm_ids_diag_hz)'
+      do n=1,size(model%info%diagnostic_variables_hz)
+         err = nf90_def_var(ncid,model%info%diagnostic_variables_hz(n)%name,NCDF_FLOAT_PRECISION,f3_dims,fabm_ids_diag_hz(n))
+         if (err .NE.  NF90_NOERR) go to 10
+         call set_attributes(ncid,fabm_ids_diag_hz(n), &
+                          long_name    =trim(model%info%diagnostic_variables_hz(n)%longname), &
+                          units        =trim(model%info%diagnostic_variables_hz(n)%units),    &
+                          FillValue    =model%info%diagnostic_variables_hz(n)%missing_value,  &
+                          missing_value=model%info%diagnostic_variables_hz(n)%missing_value,  &
+                          valid_min    =model%info%diagnostic_variables_hz(n)%minimum,        &
+                          valid_max    =model%info%diagnostic_variables_hz(n)%maximum)
+      end do
+   end if
 #endif
 
 !  globals

src/ncdf/ncdf_restart.F90

@@ -58,6 +58,12 @@
    integer                             :: biodim_id
    integer                             :: bio_id
 #endif
+#ifdef _FABM_
+   integer                             :: fabmpeldim_id
+   integer                             :: fabmbendim_id
+   integer                             :: fabm_pel_id
+   integer                             :: fabm_ben_id
+#endif
 #endif
 
    integer                             :: xlen,ylen,zlen

src/ncdf/ncdf_rivers.F90

@@ -22,6 +22,11 @@
    use bio, only: bio_calc
    use bio_var, only: numc,var_names
    use rivers, only: river_bio
+#endif
+#ifdef _FABM_
+   !use gotm_fabm, only: fabm_calc
+   use gotm_fabm, only: model
+   use rivers, only: river_fabm
 #endif
    IMPLICIT NONE
 !
@@ -45,6 +50,10 @@
    integer, allocatable                :: bio_id(:,:)
    integer, allocatable                :: r_bio(:,:)
 #endif
+#ifdef _FABM_