Fix porosity read (#201)

* Fixing the reading of variable porosity input field in preparation for the NorESM 2.0.6 release

Cherry-picked from private Ncycleprivate branch 0d56930e2fdd62caba964d375b57304942568926

* Provide number of layers (3rd dim) via ks and not hard-coded

* minor clean-up

hamocc/mo_read_sedpor.F90

@@ -51,9 +51,10 @@ module mo_read_sedpor
 
 subroutine read_sedpor(kpie,kpje,ks,omask,sed_por)
   use mod_xc,         only: mnproc,xchalt
-  use mod_dia,        only: iotype
   use mo_control_bgc, only: io_stdo_bgc,l_3Dvarsedpor
-  use mod_nctools,    only: ncfopn,ncread,ncfcls
+  use netcdf,         only: nf90_noerr,nf90_nowrite,nf90_close,nf90_open
+
+
 
   implicit none
 
@@ -62,9 +63,10 @@ subroutine read_sedpor(kpie,kpje,ks,omask,sed_por)
   real,    intent(inout) :: sed_por(kpie,kpje,ks)
 
   !local variables
-  integer :: i,j,k,errstat,dummymask(2)
+  integer :: i,j,k
   real    :: sed_por_in(kpie,kpje,ks)
   logical :: file_exists = .false.
+  integer :: ncid,ncstat
 
   ! Return if l_3Dvarsedpor is turned off
   if (.not. l_3Dvarsedpor) then
@@ -90,15 +92,36 @@ subroutine read_sedpor(kpie,kpje,ks,omask,sed_por)
     write(io_stdo_bgc,*) 'read_sedpor: read sediment porosity from ',       & 
                           trim(sedporfile)
   endif
-  call ncfopn(trim(sedporfile),'r',' ',1,iotype)
-  call ncread('sedpor',sed_por_in,dummymask,0,0.)
-  call ncfcls
+
+  ! Open netCDF data file     
+  IF(mnproc==1) THEN
+    ncstat = NF90_OPEN(trim(sedporfile),NF90_NOWRITE, ncid)
+    IF (ncstat.NE.NF90_NOERR ) THEN
+      CALL xchalt('(read_sedpor: Problem with netCDF1)')
+             stop '(read_sedpor: Problem with netCDF1)'
+    END IF
+  END IF
+
+  ! Read  data
+  call read_netcdf_var(ncid,'sedpor',sed_por_in(1,1,1),ks,0,0)
+
+  ! Close file
+  IF(mnproc==1) THEN
+    ncstat = NF90_CLOSE(ncid)
+    IF ( ncstat .NE. NF90_NOERR ) THEN
+      CALL xchalt('(read_sedpor: Problem with netCDF200)')
+             stop '(read_sedpor: Problem with netCDF200)'
+    END IF
+  END IF
+
 
   do k=1,ks
   do j=1,kpje
   do i=1,kpie
      if(omask(i,j).gt. 0.5)then
        sed_por(i,j,k)=sed_por_in(i,j,k)
+     else
+       sed_por(i,j,k)=0.
      endif
   enddo
   enddo