Skip to content
Snippets Groups Projects
icbrst.F90 22.5 KiB
Newer Older
Guillaume Samson's avatar
Guillaume Samson committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428
MODULE icbrst
   !!======================================================================
   !!                       ***  MODULE  icbrst  ***
   !! Ocean physics:  read and write iceberg restart files
   !!======================================================================
   !! History : 3.3.1 !  2010-01  (Martin&Adcroft) Original code
   !!            -    !  2011-03  (Madec)          Part conversion to NEMO form
   !!            -    !                            Removal of mapping from another grid
   !!            -    !  2011-04  (Alderson)       Split into separate modules
   !!            -    !  2011-04  (Alderson)       Restore restart routine
   !!            -    !                            Currently needs a fixed processor
   !!            -    !                            layout between restarts
   !!            -    !  2015-11  Dave Storkey     Convert icb_rst_read to use IOM so can
   !!                                              read single restart files
   !!----------------------------------------------------------------------

   !!----------------------------------------------------------------------
   !!   icb_rst_read    : read restart file
   !!   icb_rst_write   : write restart file
   !!----------------------------------------------------------------------
   USE par_oce        ! NEMO parameters
   USE dom_oce        ! NEMO domain
   USE in_out_manager ! NEMO IO routines
   USE lib_mpp        ! NEMO MPI library, lk_mpp in particular
   USE netcdf         ! netcdf routines for IO
   USE iom
   USE icb_oce        ! define iceberg arrays
   USE icbutl         ! iceberg utility routines

   IMPLICIT NONE
   PRIVATE

   PUBLIC   icb_rst_read    ! routine called in icbini.F90 module
   PUBLIC   icb_rst_write   ! routine called in icbstp.F90 module
   
   INTEGER ::   nlonid, nlatid, nxid, nyid, nuvelid, nvvelid
   INTEGER ::   nmassid, nthicknessid, nwidthid, nlengthid
   INTEGER ::   nyearid, ndayid
   INTEGER ::   nscaling_id, nmass_of_bits_id, nheat_density_id, numberid
   INTEGER ::   nsiceid, nsheatid, ncalvid, ncalvhid, nkountid
   INTEGER ::   nret, ncid, nc_dim
   
   INTEGER,  DIMENSION(3)                  :: nstrt3, nlngth3

   !!----------------------------------------------------------------------
   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
   !! $Id: icbrst.F90 15088 2021-07-06 13:03:34Z acc $
   !! Software governed by the CeCILL license (see ./LICENSE)
   !!----------------------------------------------------------------------
CONTAINS

   SUBROUTINE icb_rst_read()
      !!----------------------------------------------------------------------
      !!                 ***  SUBROUTINE icb_rst_read  ***
      !!
      !! ** Purpose :   read a iceberg restart file
      !!      NB: for this version, we just read back in the restart for this processor
      !!      so we cannot change the processor layout currently with iceberg code
      !!----------------------------------------------------------------------
      INTEGER                      ::   idim, ivar, iatt
      INTEGER                      ::   jn, iunlim_dim, ibergs_in_file
      INTEGER                      ::   ii, ij, iclass, ibase_err, imax_icb
      REAL(wp), DIMENSION(nkounts) ::   zdata      
      LOGICAL                      ::   ll_found_restart
      CHARACTER(len=256)           ::   cl_path
      CHARACTER(len=256)           ::   cl_filename
      CHARACTER(len=NF90_MAX_NAME) ::   cl_dname
      TYPE(iceberg)                ::   localberg ! NOT a pointer but an actual local variable
      TYPE(point)                  ::   localpt   ! NOT a pointer but an actual local variable
      !!----------------------------------------------------------------------
      ! Find a restart file. Assume iceberg restarts in same directory as ocean restarts
      ! and are called TRIM(cn_ocerst)//'_icebergs'
      cl_path = TRIM(cn_icbrst_indir)
      IF( cl_path(LEN_TRIM(cl_path):) /= '/' ) cl_path = TRIM(cl_path) // '/'
      cl_filename = TRIM(cn_icbrst_in)
      CALL iom_open( TRIM(cl_path)//cl_filename, ncid )

      imax_icb = 0
      IF( iom_file(ncid)%iduld .GE. 0) THEN

         ibergs_in_file = iom_file(ncid)%lenuld
         DO jn = 1,ibergs_in_file

            ! iom_get treats the unlimited dimension as time. Here the unlimited dimension 
            ! is the iceberg index, but we can still use the ktime keyword to get the iceberg we want. 

            CALL iom_get( ncid, 'xi'     ,localpt%xi  , ktime=jn )
            CALL iom_get( ncid, 'yj'     ,localpt%yj  , ktime=jn )

            ii = INT( localpt%xi + 0.5 ) + ( nn_hls-1 )
            ij = INT( localpt%yj + 0.5 ) + ( nn_hls-1 )
            ! Only proceed if this iceberg is on the local processor (excluding halos).
            IF ( ii >= mig(Nis0) .AND. ii <= mig(Nie0) .AND.   &
           &     ij >= mjg(Njs0) .AND. ij <= mjg(Nje0) ) THEN           

               CALL iom_get( ncid, jpdom_unknown, 'number', zdata(:) , ktime=jn, kstart=(/1/), kcount=(/nkounts/) )
               localberg%number(:) = INT(zdata(:))
               imax_icb = MAX( imax_icb, INT(zdata(1)) )
               CALL iom_get( ncid, 'mass_scaling' , localberg%mass_scaling, ktime=jn )
               CALL iom_get( ncid, 'lon'          , localpt%lon           , ktime=jn )
               CALL iom_get( ncid, 'lat'          , localpt%lat           , ktime=jn )
               CALL iom_get( ncid, 'uvel'         , localpt%uvel          , ktime=jn )
               CALL iom_get( ncid, 'vvel'         , localpt%vvel          , ktime=jn )
               CALL iom_get( ncid, 'mass'         , localpt%mass          , ktime=jn )
               CALL iom_get( ncid, 'thickness'    , localpt%thickness     , ktime=jn )
               CALL iom_get( ncid, 'width'        , localpt%width         , ktime=jn )
               CALL iom_get( ncid, 'length'       , localpt%length        , ktime=jn )
               CALL iom_get( ncid, 'year'         , zdata(1)              , ktime=jn )
               localpt%year = INT(zdata(1))
               CALL iom_get( ncid, 'day'          , localpt%day           , ktime=jn )
               CALL iom_get( ncid, 'mass_of_bits' , localpt%mass_of_bits  , ktime=jn )
               CALL iom_get( ncid, 'heat_density' , localpt%heat_density  , ktime=jn )
               !
               CALL icb_utl_add( localberg, localpt )
               !
            ENDIF
            !
         END DO
         !
      ELSE
         ibergs_in_file = 0
      ENDIF 

      ! Gridded variables
      CALL iom_get( ncid, jpdom_auto,    'calving'     , src_calving  )
      CALL iom_get( ncid, jpdom_auto,    'calving_hflx', src_calving_hflx  )
      CALL iom_get( ncid, jpdom_auto,    'stored_heat' , berg_grid%stored_heat  )
      ! with jpdom_auto_xy, ue use only the third element of kstart and kcount.
      CALL iom_get( ncid, jpdom_auto_xy, 'stored_ice'  , berg_grid%stored_ice, kstart=(/-99,-99,1/), kcount=(/-99,-99,nclasses/) )
      
      CALL iom_get( ncid, jpdom_unknown, 'kount' , zdata(:) )
      num_bergs(:) = INT(zdata(:))
      !

      ! Sanity checks
      jn = icb_utl_count()
      IF ( lwp .AND. nn_verbose_level >= 0 )   &
         WRITE(numout,'(2(a,i5))') 'icebergs, read_restart_bergs: # bergs =',jn,' on PE',narea-1
      IF( lk_mpp ) THEN
         ! Only mpp_sum ibergs_in_file if we are reading from multiple restart files. 
         IF( INDEX(iom_file(ncid)%name,'icebergs.nc' ) .EQ. 0 ) CALL mpp_sum('icbrst', ibergs_in_file)
         CALL mpp_sum('icbrst', jn)
      ENDIF
      IF( lwp )   WRITE(numout,'(a,i5,a,i5,a)') 'icebergs, icb_rst_read: there were',ibergs_in_file,   &
         &                                    ' bergs in the restart file and', jn,' bergs have been read'
      ! Close file
      CALL iom_close( ncid )
      !
      ! Confirm that all areas have a suitable base for assigning new iceberg
      ! numbers. This will not be the case if restarting from a collated dataset
      ! (even if using the same processor decomposition)
      !
      ibase_err = 0
      IF( num_bergs(1) < 0 .AND. num_bergs(1) /= narea - jpnij ) THEN
         ! If this area has never calved a new berg then the base should be
         ! set to narea - jpnij. If it is negative but something else then
         ! a new base will be needed to guarantee unique, future iceberg numbers
         ibase_err = 1
      ELSEIF( MOD( num_bergs(1) - narea , jpnij ) /= 0 ) THEN
         ! If this area has a base which is not in the set {narea + N*jpnij}
         ! for positive integers N then a new base will be needed to guarantee 
         ! unique, future iceberg numbers
         ibase_err = 1
      ENDIF
      IF( lk_mpp ) THEN
         CALL mpp_sum('icbrst', ibase_err)
      ENDIF
      IF( ibase_err > 0 ) THEN
         ! 
         ! A new base is needed. The only secure solution is to set bases such that
         ! all future icebergs numbers will be greater than the current global maximum
         IF( lk_mpp ) THEN
            CALL mpp_max('icbrst', imax_icb)
         ENDIF
         num_bergs(1) = imax_icb - jpnij + narea
      ENDIF
      !
      IF( lwp .AND. nn_verbose_level >= 0 )  WRITE(numout,'(a)') 'icebergs, icb_rst_read: completed'
      !
   END SUBROUTINE icb_rst_read


   SUBROUTINE icb_rst_write( kt )
      !!----------------------------------------------------------------------
      !!                 ***  SUBROUTINE icb_rst_write  ***
      !!
      !!----------------------------------------------------------------------
      INTEGER, INTENT( in ) :: kt
      !
      INTEGER ::   jn   ! dummy loop index
      INTEGER ::   idg  ! number of digits
      INTEGER ::   ix_dim, iy_dim, ik_dim, in_dim
      CHARACTER(len=256)     :: cl_path
      CHARACTER(len=256)     :: cl_filename
      CHARACTER(len=8  )     :: cl_kt
      CHARACTER(LEN=12 )     :: clfmt            ! writing format
      TYPE(iceberg), POINTER :: this
      TYPE(point)  , POINTER :: pt
      !!----------------------------------------------------------------------

      ! Following the normal restart procedure, this routine will be called
      ! the timestep before a restart stage as well as the restart timestep.
      ! This is a performance step enabling the file to be opened and contents
      ! defined in advance of the write. This is not possible with icebergs
      ! since the number of bergs to be written could change between timesteps
      IF( kt == nitrst ) THEN
         ! Only operate on the restart timestep itself.
         ! Assume we write iceberg restarts to same directory as ocean restarts.
         !
         ! directory name
         cl_path = TRIM(cn_icbrst_outdir)
         IF( cl_path(LEN_TRIM(cl_path):) /= '/' ) cl_path = TRIM(cl_path) // '/'
         !
         ! file name
         WRITE(cl_kt, '(i8.8)') kt
         cl_filename = TRIM(cexper)//"_"//cl_kt//"_"//TRIM(cn_icbrst_out)
         IF( lk_mpp ) THEN
            idg = MAX( INT(LOG10(REAL(MAX(1,jpnij-1),wp))) + 1, 4 )          ! how many digits to we need to write? min=4, max=9
            WRITE(clfmt, "('(a,a,i', i1, '.', i1, ',a)')") idg, idg          ! '(a,a,ix.x,a)'
            WRITE(cl_filename,  clfmt) TRIM(cl_filename), '_', narea-1, '.nc'
         ELSE
            WRITE(cl_filename,'(a,a)') TRIM(cl_filename),               '.nc'
         ENDIF

         IF ( lwp .AND. nn_verbose_level >= 0) WRITE(numout,'(2a)') 'icebergs, write_restart: creating ',  &
           &                                                         TRIM(cl_path)//TRIM(cl_filename)
   
         nret = NF90_CREATE(TRIM(cl_path)//TRIM(cl_filename), NF90_CLOBBER, ncid)
         IF (nret .ne. NF90_NOERR) CALL ctl_stop('icebergs, write_restart: nf_create failed')
   
         ! Dimensions
         nret = NF90_DEF_DIM(ncid, 'x', Ni_0, ix_dim)
         IF (nret .ne. NF90_NOERR) CALL ctl_stop('icebergs, write_restart: nf_def_dim x failed')
   
         nret = NF90_DEF_DIM(ncid, 'y', Nj_0, iy_dim)
         IF (nret .ne. NF90_NOERR) CALL ctl_stop('icebergs, write_restart: nf_def_dim y failed')
   
         nret = NF90_DEF_DIM(ncid, 'c', nclasses, nc_dim)
         IF (nret .ne. NF90_NOERR) CALL ctl_stop('icebergs, write_restart: nf_def_dim c failed')
   
         nret = NF90_DEF_DIM(ncid, 'k', nkounts, ik_dim)
         IF (nret .ne. NF90_NOERR) CALL ctl_stop('icebergs, write_restart: nf_def_dim k failed')
   
         ! global attributes
         IF( lk_mpp ) THEN
            ! Set domain parameters (assume jpdom_local_full)
            nret = NF90_PUT_ATT( ncid, NF90_GLOBAL, 'DOMAIN_number_total'   , jpnij                        )
            nret = NF90_PUT_ATT( ncid, NF90_GLOBAL, 'DOMAIN_number'         , narea-1                      )
            nret = NF90_PUT_ATT( ncid, NF90_GLOBAL, 'DOMAIN_dimensions_ids' , (/ 1         , 2          /) )
            nret = NF90_PUT_ATT( ncid, NF90_GLOBAL, 'DOMAIN_size_global'    , (/ Ni0glo    , Nj0glo     /) )
            nret = NF90_PUT_ATT( ncid, NF90_GLOBAL, 'DOMAIN_size_local'     , (/ Ni_0      , Nj_0       /) )
            nret = NF90_PUT_ATT( ncid, NF90_GLOBAL, 'DOMAIN_position_first' , (/ mig0(Nis0), mjg0(Njs0) /) )
            nret = NF90_PUT_ATT( ncid, NF90_GLOBAL, 'DOMAIN_position_last'  , (/ mig0(Nie0), mjg0(Nje0) /) )
            nret = NF90_PUT_ATT( ncid, NF90_GLOBAL, 'DOMAIN_halo_size_start', (/ 0         , 0          /) )
            nret = NF90_PUT_ATT( ncid, NF90_GLOBAL, 'DOMAIN_halo_size_end'  , (/ 0         , 0          /) )
            nret = NF90_PUT_ATT( ncid, NF90_GLOBAL, 'DOMAIN_type'           , 'BOX'                        )
         ENDIF
         
         IF (associated(first_berg)) then
            nret = NF90_DEF_DIM(ncid, 'n', NF90_UNLIMITED, in_dim)
            IF (nret .ne. NF90_NOERR) CALL ctl_stop('icebergs, write_restart: nf_def_dim n failed')
         ENDIF
   
         ! Variables
         nret = NF90_DEF_VAR(ncid, 'kount'       , NF90_INT   , (/ ik_dim /), nkountid)
         nret = NF90_DEF_VAR(ncid, 'calving'     , NF90_DOUBLE, (/ ix_dim, iy_dim /), ncalvid)
         nret = NF90_DEF_VAR(ncid, 'calving_hflx', NF90_DOUBLE, (/ ix_dim, iy_dim /), ncalvhid)
         nret = NF90_DEF_VAR(ncid, 'stored_ice'  , NF90_DOUBLE, (/ ix_dim, iy_dim, nc_dim /), nsiceid)
         nret = NF90_DEF_VAR(ncid, 'stored_heat' , NF90_DOUBLE, (/ ix_dim, iy_dim /), nsheatid)
   
         ! Attributes
         nret = NF90_PUT_ATT(ncid, ncalvid , 'long_name', 'iceberg calving')
         nret = NF90_PUT_ATT(ncid, ncalvid , 'units', 'some')
         nret = NF90_PUT_ATT(ncid, ncalvhid, 'long_name', 'heat flux associated with iceberg calving')
         nret = NF90_PUT_ATT(ncid, ncalvhid, 'units', 'some')
         nret = NF90_PUT_ATT(ncid, nsiceid , 'long_name', 'stored ice used to calve icebergs')
         nret = NF90_PUT_ATT(ncid, nsiceid , 'units', 'kg/s')
         nret = NF90_PUT_ATT(ncid, nsheatid, 'long_name', 'heat in stored ice used to calve icebergs')
         nret = NF90_PUT_ATT(ncid, nsheatid, 'units', 'J/kg/s')
   
         IF ( ASSOCIATED(first_berg) ) THEN
   
            ! Only add berg variables for this PE if we have anything to say
   
            ! Variables
            nret = NF90_DEF_VAR(ncid, 'lon', NF90_DOUBLE, in_dim, nlonid)
            nret = NF90_DEF_VAR(ncid, 'lat', NF90_DOUBLE, in_dim, nlatid)
            nret = NF90_DEF_VAR(ncid, 'xi', NF90_DOUBLE, in_dim, nxid)
            nret = NF90_DEF_VAR(ncid, 'yj', NF90_DOUBLE, in_dim, nyid)
            nret = NF90_DEF_VAR(ncid, 'uvel', NF90_DOUBLE, in_dim, nuvelid)
            nret = NF90_DEF_VAR(ncid, 'vvel', NF90_DOUBLE, in_dim, nvvelid)
            nret = NF90_DEF_VAR(ncid, 'mass', NF90_DOUBLE, in_dim, nmassid)
            nret = NF90_DEF_VAR(ncid, 'thickness', NF90_DOUBLE, in_dim, nthicknessid)
            nret = NF90_DEF_VAR(ncid, 'width', NF90_DOUBLE, in_dim, nwidthid)
            nret = NF90_DEF_VAR(ncid, 'length', NF90_DOUBLE, in_dim, nlengthid)
            nret = NF90_DEF_VAR(ncid, 'number', NF90_INT, (/ik_dim,in_dim/), numberid)
            nret = NF90_DEF_VAR(ncid, 'year', NF90_INT, in_dim, nyearid)
            nret = NF90_DEF_VAR(ncid, 'day', NF90_DOUBLE, in_dim, ndayid)
            nret = NF90_DEF_VAR(ncid, 'mass_scaling', NF90_DOUBLE, in_dim, nscaling_id)
            nret = NF90_DEF_VAR(ncid, 'mass_of_bits', NF90_DOUBLE, in_dim, nmass_of_bits_id)
            nret = NF90_DEF_VAR(ncid, 'heat_density', NF90_DOUBLE, in_dim, nheat_density_id)
   
            ! Attributes
            nret = NF90_PUT_ATT(ncid, nlonid, 'long_name', 'longitude')
            nret = NF90_PUT_ATT(ncid, nlonid, 'units', 'degrees_E')
            nret = NF90_PUT_ATT(ncid, nlatid, 'long_name', 'latitude')
            nret = NF90_PUT_ATT(ncid, nlatid, 'units', 'degrees_N')
            nret = NF90_PUT_ATT(ncid, nxid, 'long_name', 'x grid box position')
            nret = NF90_PUT_ATT(ncid, nxid, 'units', 'fractional')
            nret = NF90_PUT_ATT(ncid, nyid, 'long_name', 'y grid box position')
            nret = NF90_PUT_ATT(ncid, nyid, 'units', 'fractional')
            nret = NF90_PUT_ATT(ncid, nuvelid, 'long_name', 'zonal velocity')
            nret = NF90_PUT_ATT(ncid, nuvelid, 'units', 'm/s')
            nret = NF90_PUT_ATT(ncid, nvvelid, 'long_name', 'meridional velocity')
            nret = NF90_PUT_ATT(ncid, nvvelid, 'units', 'm/s')
            nret = NF90_PUT_ATT(ncid, nmassid, 'long_name', 'mass')
            nret = NF90_PUT_ATT(ncid, nmassid, 'units', 'kg')
            nret = NF90_PUT_ATT(ncid, nthicknessid, 'long_name', 'thickness')
            nret = NF90_PUT_ATT(ncid, nthicknessid, 'units', 'm')
            nret = NF90_PUT_ATT(ncid, nwidthid, 'long_name', 'width')
            nret = NF90_PUT_ATT(ncid, nwidthid, 'units', 'm')
            nret = NF90_PUT_ATT(ncid, nlengthid, 'long_name', 'length')
            nret = NF90_PUT_ATT(ncid, nlengthid, 'units', 'm')
            nret = NF90_PUT_ATT(ncid, numberid, 'long_name', 'iceberg number on this processor')
            nret = NF90_PUT_ATT(ncid, numberid, 'units', 'count')
            nret = NF90_PUT_ATT(ncid, nyearid, 'long_name', 'calendar year of calving event')
            nret = NF90_PUT_ATT(ncid, nyearid, 'units', 'years')
            nret = NF90_PUT_ATT(ncid, ndayid, 'long_name', 'year day of calving event')
            nret = NF90_PUT_ATT(ncid, ndayid, 'units', 'days')
            nret = NF90_PUT_ATT(ncid, nscaling_id, 'long_name', 'scaling factor for mass of calving berg')
            nret = NF90_PUT_ATT(ncid, nscaling_id, 'units', 'none')
            nret = NF90_PUT_ATT(ncid, nmass_of_bits_id, 'long_name', 'mass of bergy bits')
            nret = NF90_PUT_ATT(ncid, nmass_of_bits_id, 'units', 'kg')
            nret = NF90_PUT_ATT(ncid, nheat_density_id, 'long_name', 'heat density')
            nret = NF90_PUT_ATT(ncid, nheat_density_id, 'units', 'J/kg')
   
         ENDIF ! associated(first_berg)
   
         ! End define mode
         nret = NF90_ENDDEF(ncid)
   
         ! --------------------------------
         ! now write some data
   
         nstrt3(1) = 1
         nstrt3(2) = 1
         nlngth3(1) = Ni_0
         nlngth3(2) = Nj_0
         nlngth3(3) = 1
   
         DO jn=1,nclasses
            nstrt3(3) = jn
            nret = NF90_PUT_VAR( ncid, nsiceid, berg_grid%stored_ice(Nis0:Nie0,Njs0:Nje0,jn), nstrt3, nlngth3 )
            IF (nret .ne. NF90_NOERR) THEN
               IF( lwp ) WRITE(numout,*) TRIM(NF90_STRERROR( nret ))
               CALL ctl_stop('icebergs, write_restart: nf_put_var stored_ice failed')
            ENDIF
         ENDDO
         IF( lwp ) WRITE(numout,*) 'file: ',TRIM(cl_path)//TRIM(cl_filename),' var: stored_ice  written'
   
         nret = NF90_PUT_VAR( ncid, nkountid, num_bergs(:) )
         IF (nret .ne. NF90_NOERR) CALL ctl_stop('icebergs, write_restart: nf_put_var kount failed')
   
         nret = NF90_PUT_VAR( ncid, nsheatid, berg_grid%stored_heat(Nis0:Nie0,Njs0:Nje0) )
         IF (nret .ne. NF90_NOERR) CALL ctl_stop('icebergs, write_restart: nf_put_var stored_heat failed')
         IF( lwp ) WRITE(numout,*) 'file: ',TRIM(cl_path)//TRIM(cl_filename),' var: stored_heat written'
   
         nret = NF90_PUT_VAR( ncid, ncalvid , src_calving(Nis0:Nie0,Njs0:Nje0) )
         IF (nret .ne. NF90_NOERR) CALL ctl_stop('icebergs, write_restart: nf_put_var calving failed')
         nret = NF90_PUT_VAR( ncid, ncalvhid, src_calving_hflx(Nis0:Nie0,Njs0:Nje0) )
         IF (nret .ne. NF90_NOERR) CALL ctl_stop('icebergs, write_restart: nf_put_var calving_hflx failed')
         IF( lwp ) WRITE(numout,*) 'file: ',TRIM(cl_path)//TRIM(cl_filename),' var: calving written'
   
         IF ( ASSOCIATED(first_berg) ) THEN
   
            ! Write variables
            ! just write out the current point of the trajectory
   
            this => first_berg
            jn = 0
            DO WHILE (ASSOCIATED(this))
               pt => this%current_point
               jn=jn+1
   
               nret = NF90_PUT_VAR(ncid, numberid, this%number, (/1,jn/), (/nkounts,1/) )
               nret = NF90_PUT_VAR(ncid, nscaling_id, this%mass_scaling, (/ jn /) )
   
               nret = NF90_PUT_VAR(ncid, nlonid, pt%lon, (/ jn /) )
               nret = NF90_PUT_VAR(ncid, nlatid, pt%lat, (/ jn /) )
               nret = NF90_PUT_VAR(ncid, nxid, pt%xi, (/ jn /) )
               nret = NF90_PUT_VAR(ncid, nyid, pt%yj, (/ jn /) )
               nret = NF90_PUT_VAR(ncid, nuvelid, pt%uvel, (/ jn /) )
               nret = NF90_PUT_VAR(ncid, nvvelid, pt%vvel, (/ jn /) )
               nret = NF90_PUT_VAR(ncid, nmassid, pt%mass, (/ jn /) )
               nret = NF90_PUT_VAR(ncid, nthicknessid, pt%thickness, (/ jn /) )
               nret = NF90_PUT_VAR(ncid, nwidthid, pt%width, (/ jn /) )
               nret = NF90_PUT_VAR(ncid, nlengthid, pt%length, (/ jn /) )
               nret = NF90_PUT_VAR(ncid, nyearid, pt%year, (/ jn /) )
               nret = NF90_PUT_VAR(ncid, ndayid, pt%day, (/ jn /) )
               nret = NF90_PUT_VAR(ncid, nmass_of_bits_id, pt%mass_of_bits, (/ jn /) )
               nret = NF90_PUT_VAR(ncid, nheat_density_id, pt%heat_density, (/ jn /) )
   
               this=>this%next
            END DO
            !
         ENDIF ! associated(first_berg)
   
         ! Finish up
         nret = NF90_CLOSE(ncid)
         IF (nret /= NF90_NOERR) CALL ctl_stop('icebergs, write_restart: nf_close failed')
   
         ! Sanity check
         jn = icb_utl_count()
         IF ( lwp .AND. nn_verbose_level >= 0)   &
            WRITE(numout,'(2(a,i5))') 'icebergs, icb_rst_write: # bergs =',jn,' on PE',narea-1
         IF( lk_mpp ) THEN
            CALL mpp_sum('icbrst', jn)
         ENDIF
         IF(lwp)   WRITE(numout,'(a,i5,a,i5,a)') 'icebergs, icb_rst_write: ', jn,   &
            &                                    ' bergs in total have been written at timestep ', kt
         !
         ! Finish up
         !
      ENDIF
   END SUBROUTINE icb_rst_write
   !
   !!======================================================================
END MODULE icbrst