Merge pull request CICE-Consortium#66 from apcraig/cgtestI

Fix C/CD restart and decomp bit-for-bit issues

cicecore/cicedynB/dynamics/ice_dyn_evp.F90

@@ -971,7 +971,9 @@ subroutine evp (dt)
 
             call ice_timer_start(timer_bound)
             call stack_velocity_field(uvel, vvel, fld2)
-            if (maskhalo_dyn) then
+            ! maskhalo_dyn causes non bit-for-bit results on different decomps 
+            ! with C/CD in some cases
+            if (grid_ice == 'B' .and. maskhalo_dyn) then
                call ice_HaloUpdate (fld2,               halo_info_mask, &
                                     field_loc_NEcorner, field_type_vector)
             else

cicecore/cicedynB/infrastructure/ice_domain.F90

@@ -287,7 +287,7 @@ end subroutine init_domain_blocks
 
 !***********************************************************************
 
- subroutine init_domain_distribution(KMTG,ULATG)
+ subroutine init_domain_distribution(KMTG,ULATG,grid_ice)
 
 !  This routine calls appropriate setup routines to distribute blocks
 !  across processors and defines arrays with block ids for any local
@@ -302,6 +302,9 @@ subroutine init_domain_distribution(KMTG,ULATG)
       KMTG           ,&! global topography
       ULATG            ! global latitude field (radians)
 
+   character(len=*), intent(in) :: &
+      grid_ice         ! grid_ice, B, C, CD, etc
+
 !----------------------------------------------------------------------
 !
 !  local variables
@@ -319,6 +322,7 @@ subroutine init_domain_distribution(KMTG,ULATG)
    integer (int_kind) :: &
       i,j,n              ,&! dummy loop indices
       ig,jg              ,&! global indices
+      igm1,igp1,jgm1,jgp1,&! global indices
       ninfo              ,&! ice_distributionGet check
       work_unit          ,&! size of quantized work unit
 #ifdef USE_NETCDF
@@ -509,10 +513,25 @@ subroutine init_domain_distribution(KMTG,ULATG)
                   if (this_block%i_glob(i) > 0) then
                      ig = this_block%i_glob(i)
                      jg = this_block%j_glob(j)
-                     if (KMTG(ig,jg) > puny .and.                      &
-                        (ULATG(ig,jg) < shlat/rad_to_deg .or.          &
-                         ULATG(ig,jg) > nhlat/rad_to_deg) )            & 
-                          nocn(n) = nocn(n) + flat(ig,jg)
+                     if (grid_ice == 'C' .or. grid_ice == 'CD') then
+                        ! Have to be careful about block elimination with C/CD
+                        ! Use a bigger stencil
+                        igm1 = mod(ig-2+nx_global,nx_global)+1
+                        igp1 = mod(ig,nx_global)+1
+                        jgm1 = max(jg-1,1)
+                        jgp1 = min(jg+1,ny_global)
+                        if ((KMTG(ig  ,jg  ) > puny .or.                             &
+                             KMTG(igm1,jg  ) > puny .or. KMTG(igp1,jg  ) > puny .or. &
+                             KMTG(ig  ,jgp1) > puny .or. KMTG(ig  ,jgm1) > puny) .and. &
+                            (ULATG(ig,jg) < shlat/rad_to_deg .or.          &
+                             ULATG(ig,jg) > nhlat/rad_to_deg) )            & 
+                             nocn(n) = nocn(n) + flat(ig,jg)
+                     else
+                        if (KMTG(ig,jg) > puny .and.                      &
+                           (ULATG(ig,jg) < shlat/rad_to_deg .or.          &
+                            ULATG(ig,jg) > nhlat/rad_to_deg) )            & 
+                             nocn(n) = nocn(n) + flat(ig,jg)
+                     endif
                   endif
                end do
             endif
@@ -529,10 +548,8 @@ subroutine init_domain_distribution(KMTG,ULATG)
          ! Keep all blocks even the ones only containing land points
          if (distribution_wght == 'block') nocn(n) = nx_block*ny_block
 #else
-         if (distribution_wght == 'block' .and. &   ! POP style
-             nocn(n) > 0) nocn(n) = nx_block*ny_block
-         if (distribution_wght == 'blockall') &
-             nocn(n) = nx_block*ny_block
+         if (distribution_wght == 'block' .and. nocn(n) > 0) nocn(n) = nx_block*ny_block
+         if (.not. landblockelim) nocn(n) = max(nocn(n),1)
 #endif
       end do
    endif  ! distribution_wght = file

cicecore/cicedynB/infrastructure/ice_grid.F90

@@ -400,7 +400,7 @@ subroutine init_grid1
       ! distribute blocks among processors
       !-----------------------------------------------------------------
 
-      call init_domain_distribution(work_g2, work_g1)  ! KMT, ULAT
+      call init_domain_distribution(work_g2, work_g1, grid_ice)  ! KMT, ULAT
 
       deallocate(work_g1)
       deallocate(work_g2)

cicecore/cicedynB/infrastructure/ice_read_write.F90

@@ -1670,7 +1670,7 @@ subroutine ice_read_nc_xyf(fid,  nrec,  varname, work,  diag, &
             amin = minval(work_g1(:,:,n))
             amax = maxval(work_g1(:,:,n), mask = work_g1(:,:,n) /= missingvalue)
             asum = sum   (work_g1(:,:,n), mask = work_g1(:,:,n) /= missingvalue)
-            write(nu_diag,*) subname,' min, max, sum =', amin, amax, asum
+            write(nu_diag,*) subname,' min, max, sum =', amin, amax, asum, trim(varname)
          enddo
       endif
 

cicecore/cicedynB/infrastructure/ice_restart_driver.F90

@@ -19,6 +19,7 @@ module ice_restart_driver
 
       use ice_kinds_mod
       use ice_arrays_column, only: oceanmixed_ice
+      use ice_communicate, only: my_task, master_task
       use ice_constants, only: c0, c1, p5, &
           field_loc_center, field_loc_NEcorner, &
           field_loc_Eface, field_loc_Nface, &
@@ -222,29 +223,29 @@ subroutine dumpfile(filename_spec)
 
       if (grid_ice == 'CD' .or. grid_ice == 'C') then
 
-      !$OMP PARALLEL DO PRIVATE(iblk,i,j)
-      do iblk = 1, nblocks
-         do j = 1, ny_block
-         do i = 1, nx_block
-            work1(i,j,iblk) = c0
-            if (icenmask(i,j,iblk)) work1(i,j,iblk) = c1
-         enddo
-         enddo
-      enddo
-      !$OMP END PARALLEL DO
-      call write_restart_field(nu_dump,0,work1,'ruf8','icenmask',1,diag)
-
-      !$OMP PARALLEL DO PRIVATE(iblk,i,j)
-      do iblk = 1, nblocks
-         do j = 1, ny_block
-         do i = 1, nx_block
-            work1(i,j,iblk) = c0
-            if (iceemask(i,j,iblk)) work1(i,j,iblk) = c1
+         !$OMP PARALLEL DO PRIVATE(iblk,i,j)
+         do iblk = 1, nblocks
+            do j = 1, ny_block
+            do i = 1, nx_block
+               work1(i,j,iblk) = c0
+               if (icenmask(i,j,iblk)) work1(i,j,iblk) = c1
+            enddo
+            enddo
          enddo
+         !$OMP END PARALLEL DO
+         call write_restart_field(nu_dump,0,work1,'ruf8','icenmask',1,diag)
+
+         !$OMP PARALLEL DO PRIVATE(iblk,i,j)
+         do iblk = 1, nblocks
+            do j = 1, ny_block
+            do i = 1, nx_block
+               work1(i,j,iblk) = c0
+               if (iceemask(i,j,iblk)) work1(i,j,iblk) = c1
+            enddo
+            enddo
          enddo
-      enddo
-      !$OMP END PARALLEL DO
-      call write_restart_field(nu_dump,0,work1,'ruf8','iceemask',1,diag)
+         !$OMP END PARALLEL DO
+         call write_restart_field(nu_dump,0,work1,'ruf8','iceemask',1,diag)
 
       endif
 
@@ -266,7 +267,6 @@ subroutine restartfile (ice_ic)
       use ice_boundary, only: ice_HaloUpdate_stress
       use ice_blocks, only: nghost, nx_block, ny_block
       use ice_calendar, only: istep0, npt, calendar
-      use ice_communicate, only: my_task, master_task
       use ice_domain, only: nblocks, halo_info
       use ice_domain_size, only: nilyr, nslyr, ncat, &
           max_blocks
@@ -697,7 +697,6 @@ subroutine restartfile_v4 (ice_ic)
       use ice_blocks, only: nghost, nx_block, ny_block
       use ice_calendar, only: istep0, istep1, timesecs, calendar, npt, &
           set_date_from_timesecs
-      use ice_communicate, only: my_task, master_task
       use ice_domain, only: nblocks, distrb_info
       use ice_domain_size, only: nilyr, nslyr, ncat, nx_global, ny_global, &
           max_blocks

cicecore/cicedynB/infrastructure/io/io_netcdf/ice_restart.F90

@@ -10,6 +10,7 @@
       module ice_restart
 
       use ice_broadcast
+      use ice_communicate, only: my_task, master_task
       use ice_kinds_mod
 #ifdef USE_NETCDF
       use netcdf
@@ -45,7 +46,6 @@ subroutine init_restart_read(ice_ic)
 
       use ice_calendar, only: msec, mmonth, mday, myear, &
                              istep0, istep1, npt
-      use ice_communicate, only: my_task, master_task
 
       character(len=char_len_long), intent(in), optional :: ice_ic
 
@@ -132,7 +132,6 @@ subroutine init_restart_write(filename_spec)
 
       use ice_blocks, only: nghost
       use ice_calendar, only: msec, mmonth, mday, myear, istep1
-      use ice_communicate, only: my_task, master_task
       use ice_domain_size, only: nx_global, ny_global, ncat, nilyr, nslyr, &
                                  n_iso, n_aero, nblyr, n_zaero, n_algae, n_doc,   &
                                  n_dic, n_don, n_fed, n_fep, nfsd
@@ -841,7 +840,6 @@ end subroutine write_restart_field
       subroutine final_restart()
 
       use ice_calendar, only: istep1, idate
-      use ice_communicate, only: my_task, master_task
 
       integer (kind=int_kind) :: status
 
@@ -904,8 +902,11 @@ logical function query_field(nu,vname)
 
       query_field = .false.
 #ifdef USE_NETCDF
-      status = nf90_inq_varid(ncid,trim(vname),varid)
-      if (status == nf90_noerr) query_field = .true.
+      if (my_task == master_task) then
+         status = nf90_inq_varid(ncid,trim(vname),varid)
+         if (status == nf90_noerr) query_field = .true.
+      endif
+      call broadcast_scalar(query_field,master_task)
 #else
       call abort_ice(subname//'ERROR: USE_NETCDF cpp not defined for '//trim(ice_ic), &
           file=__FILE__, line=__LINE__)

cicecore/cicedynB/infrastructure/io/io_pio2/ice_restart.F90

@@ -6,6 +6,7 @@
       module ice_restart
 
       use ice_broadcast
+      use ice_communicate, only: my_task, master_task
       use ice_exit, only: abort_ice
       use ice_fileunits, only: nu_diag, nu_restart, nu_rst_pointer
       use ice_kinds_mod
@@ -44,7 +45,6 @@ subroutine init_restart_read(ice_ic)
 
       use ice_calendar, only: istep0, istep1, myear, mmonth, &
                               mday, msec, npt
-      use ice_communicate, only: my_task, master_task
       use ice_domain_size, only: ncat
       use ice_read_write, only: ice_open
 
@@ -140,7 +140,6 @@ end subroutine init_restart_read
       subroutine init_restart_write(filename_spec)
 
       use ice_calendar, only: msec, mmonth, mday, myear, istep1
-      use ice_communicate, only: my_task, master_task
       use ice_domain_size, only: nx_global, ny_global, ncat, nilyr, nslyr, &
                                  n_iso, n_aero, nblyr, n_zaero, n_algae, n_doc,   &
                                  n_dic, n_don, n_fed, n_fep, nfsd
@@ -696,7 +695,6 @@ subroutine read_restart_field(nu,nrec,work,atype,vname,ndim3,diag, &
                                     field_loc, field_type)
 
       use ice_blocks, only: nx_block, ny_block
-      use ice_communicate, only: my_task, master_task
       use ice_constants, only: c0, field_loc_center
       use ice_boundary, only: ice_HaloUpdate
       use ice_domain, only: halo_info, distrb_info, nblocks
@@ -815,7 +813,6 @@ end subroutine read_restart_field
       subroutine write_restart_field(nu,nrec,work,atype,vname,ndim3,diag)
 
       use ice_blocks, only: nx_block, ny_block
-      use ice_communicate, only: my_task, master_task
       use ice_constants, only: c0, field_loc_center
       use ice_domain, only: distrb_info, nblocks
       use ice_domain_size, only: max_blocks, ncat
@@ -901,7 +898,6 @@ end subroutine write_restart_field
       subroutine final_restart()
 
       use ice_calendar, only: istep1, idate, msec
-      use ice_communicate, only: my_task, master_task
 
       character(len=*), parameter :: subname = '(final_restart)'
 
@@ -951,8 +947,11 @@ logical function query_field(nu,vname)
 
       query_field = .false.
 #ifdef USE_NETCDF
-      status = pio_inq_varid(File,trim(vname),vardesc)
-      if (status == PIO_noerr) query_field = .true.
+      if (my_task == master_task) then
+         status = pio_inq_varid(File,trim(vname),vardesc)
+         if (status == PIO_noerr) query_field = .true.
+      endif
+      call broadcast_scalar(query_field,master_task)
 #else
       call abort_ice(subname//'ERROR: USE_NETCDF cpp not defined for '//trim(ice_ic), &
           file=__FILE__, line=__LINE__)