Fix box2001 forcing, restart fields on land, bathymetry default value…

…, omp_suite (#65)

* Fix several problems noted in PR #64.  This should
address the problems with the boxrestore test errors.
- Update box2001 so it's bit-for-bit with different blocks/decomps/pe counts.
- Initialize bathymetry values at all gridcells when bathymetry_format='default' and use_bathymetry=.false.

* update comparelog to exclude extraneous icepack output, help prod comparisons

* - Zero out certain fields on land on restart files.  Some fields have non-zero values
  over land by default and this causes problems with land block elimination and
  comparisons of different decompositions.  Does not affect science.
- Update omp_suite to use different block sizes in comparisons

* use c0 instead of 0.

* fix OpenMP private variables on new code

* update pio restart diagnostics format

cicecore/cicedynB/general/ice_forcing.F90

@@ -682,7 +682,7 @@ subroutine get_forcing_atmo
       !$OMP PARALLEL DO PRIVATE(iblk,ilo,ihi,jlo,jhi,this_block)
       do iblk = 1, nblocks
 
-         this_block = get_block(blocks_ice(iblk),iblk)         
+         this_block = get_block(blocks_ice(iblk),iblk)
          ilo = this_block%ilo
          ihi = this_block%ihi
          jlo = this_block%jlo
@@ -2438,7 +2438,7 @@ subroutine LY_data
         enddo
 
       ! AOMIP
-        this_block = get_block(blocks_ice(iblk),iblk)         
+        this_block = get_block(blocks_ice(iblk),iblk)
         ilo = this_block%ilo
         ihi = this_block%ihi
         jlo = this_block%jlo
@@ -3480,7 +3480,7 @@ subroutine monthly_data
         enddo
 
       ! AOMIP
-      this_block = get_block(blocks_ice(iblk),iblk)         
+      this_block = get_block(blocks_ice(iblk),iblk)
       ilo = this_block%ilo
       ihi = this_block%ihi
       jlo = this_block%jlo
@@ -5313,10 +5313,10 @@ subroutine box2001_data_atm
 ! these are defined at the u point
 ! authors: Elizabeth Hunke, LANL
 
-      use ice_domain, only: nblocks
+      use ice_domain, only: nblocks, blocks_ice
       use ice_domain_size, only: max_blocks
       use ice_calendar, only: timesecs
-      use ice_blocks, only: nx_block, ny_block, nghost
+      use ice_blocks, only: block, get_block, nx_block, ny_block, nghost
       use ice_flux, only: uatm, vatm, wind, rhoa, strax, stray
       use ice_state, only: aice
 
@@ -5325,8 +5325,15 @@ subroutine box2001_data_atm
       integer (kind=int_kind) :: &
          iblk, i,j           ! loop indices
 
+      integer (kind=int_kind) :: &
+         iglob(nx_block), & ! global indices
+         jglob(ny_block)    ! global indices
+
+      type (block) :: &
+         this_block           ! block information for current block
+
       real (kind=dbl_kind) :: &
-          secday, pi , puny, period, pi2, tau
+         secday, pi , puny, period, pi2, tau
 
       character(len=*), parameter :: subname = '(box2001_data_atm)'
 
@@ -5341,6 +5348,10 @@ subroutine box2001_data_atm
          do j = 1, ny_block   
          do i = 1, nx_block   
 
+         this_block = get_block(blocks_ice(iblk),iblk)
+         iglob = this_block%i_glob
+         jglob = this_block%j_glob
+
 !tcraig, move to box2001_data_ocn
 !         ! ocean current
 !         ! constant in time, could be initialized in ice_flux.F90
@@ -5354,14 +5365,14 @@ subroutine box2001_data_atm
 
          ! wind components
          uatm(i,j,iblk) = c5 + (sin(pi2*timesecs/period)-c3) &
-                              * sin(pi2*real(i-nghost, kind=dbl_kind)  &
+                              * sin(pi2*real(iglob(i), kind=dbl_kind)  &
                                        /real(nx_global,kind=dbl_kind)) &
-                              * sin(pi *real(j-nghost, kind=dbl_kind)  &
+                              * sin(pi *real(jglob(j), kind=dbl_kind)  &
                                        /real(ny_global,kind=dbl_kind))
          vatm(i,j,iblk) = c5 + (sin(pi2*timesecs/period)-c3) &
-                              * sin(pi *real(i-nghost, kind=dbl_kind)  &
+                              * sin(pi *real(iglob(i), kind=dbl_kind)  &
                                        /real(nx_global,kind=dbl_kind)) &
-                              * sin(pi2*real(j-nghost, kind=dbl_kind)  &
+                              * sin(pi2*real(jglob(j), kind=dbl_kind)  &
                                        /real(ny_global,kind=dbl_kind))
          ! wind stress
          wind(i,j,iblk) = sqrt(uatm(i,j,iblk)**2 + vatm(i,j,iblk)**2)
@@ -5408,10 +5419,10 @@ subroutine box2001_data_ocn
 ! these are defined at the u point
 ! authors: Elizabeth Hunke, LANL
 
-      use ice_domain, only: nblocks
+      use ice_domain, only: nblocks, blocks_ice
       use ice_domain_size, only: max_blocks
       use ice_calendar, only: timesecs
-      use ice_blocks, only: nx_block, ny_block, nghost
+      use ice_blocks, only: block, get_block, nx_block, ny_block, nghost
       use ice_flux, only: uocn, vocn
       use ice_grid, only: uvm
 
@@ -5420,8 +5431,15 @@ subroutine box2001_data_ocn
       integer (kind=int_kind) :: &
          iblk, i,j           ! loop indices
 
+      integer (kind=int_kind) :: &
+         iglob(nx_block), & ! global indices
+         jglob(ny_block)    ! global indices
+
+      type (block) :: &
+         this_block           ! block information for current block
+
       real (kind=dbl_kind) :: &
-          secday, pi , puny, period, pi2, tau
+         secday, pi , puny, period, pi2, tau
 
       character(len=*), parameter :: subname = '(box2001_data_ocn)'
 
@@ -5431,12 +5449,16 @@ subroutine box2001_data_ocn
          do j = 1, ny_block   
          do i = 1, nx_block   
 
+         this_block = get_block(blocks_ice(iblk),iblk)
+         iglob = this_block%i_glob
+         jglob = this_block%j_glob
+
          ! ocean current
          ! constant in time, could be initialized in ice_flux.F90
-         uocn(i,j,iblk) =  p2*real(j-nghost, kind=dbl_kind) &
-                            / real(nx_global,kind=dbl_kind) - p1
-         vocn(i,j,iblk) = -p2*real(i-nghost, kind=dbl_kind) &
-                            / real(ny_global,kind=dbl_kind) + p1
+         uocn(i,j,iblk) =  p2*real(jglob(j), kind=dbl_kind) &
+                            / real(ny_global,kind=dbl_kind) - p1
+         vocn(i,j,iblk) = -p2*real(iglob(i), kind=dbl_kind) &
+                            / real(nx_global,kind=dbl_kind) + p1
 
          uocn(i,j,iblk) = uocn(i,j,iblk) * uvm(i,j,iblk)
          vocn(i,j,iblk) = vocn(i,j,iblk) * uvm(i,j,iblk)

cicecore/cicedynB/infrastructure/ice_grid.F90

@@ -4177,6 +4177,8 @@ end subroutine gridbox_verts
 
       subroutine get_bathymetry
 
+      use ice_constants, only: c0
+
       integer (kind=int_kind) :: &
          i, j, k, iblk      ! loop indices
 
@@ -4228,6 +4230,7 @@ subroutine get_bathymetry
             depth(k) = depth(k-1) + thick(k)
          enddo
 
+         bathymetry = c0
          do iblk = 1, nblocks
             do j = 1, ny_block
             do i = 1, nx_block

cicecore/cicedynB/infrastructure/ice_restart_driver.F90

@@ -62,7 +62,7 @@ subroutine dumpfile(filename_spec)
           stresspT, stressmT, stress12T, &
           stresspU, stressmU, stress12U 
       use ice_flux, only: coszen
-      use ice_grid, only: grid_ice
+      use ice_grid, only: grid_ice, tmask
       use ice_state, only: aicen, vicen, vsnon, trcrn, uvel, vvel, &
                            uvelE, vvelE, uvelN, vvelN
 
@@ -97,6 +97,20 @@ subroutine dumpfile(filename_spec)
 
       diag = .true.
 
+      !-----------------------------------------------------------------
+      ! Zero out tracers over land
+      !-----------------------------------------------------------------
+
+      !$OMP PARALLEL DO PRIVATE(iblk,i,j)
+      do iblk = 1, nblocks
+         do j = 1, ny_block
+         do i = 1, nx_block
+            if (.not. tmask(i,j,iblk)) trcrn(i,j,:,:,iblk) = c0
+         enddo
+         enddo
+      enddo
+      !$OMP END PARALLEL DO
+
       !-----------------------------------------------------------------
       ! state variables
       ! Tsfc is the only tracer written to binary files.  All other

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

@@ -788,18 +788,15 @@ subroutine read_restart_field(nu,nrec,work,atype,vname,ndim3,diag, &
                   amax = global_maxval(work(:,:,n,:),distrb_info)
                   asum = global_sum(work(:,:,n,:), distrb_info, field_loc_center)
                   if (my_task == master_task) then
-                     write(nu_diag,*) ' min and max =', amin, amax
-                     write(nu_diag,*) ' sum =',asum
+                     write(nu_diag,*) ' min, max, sum =', amin, amax, asum, trim(vname)
                   endif
                enddo
             else
                amin = global_minval(work(:,:,1,:),distrb_info)
                amax = global_maxval(work(:,:,1,:),distrb_info)
                asum = global_sum(work(:,:,1,:), distrb_info, field_loc_center)
                if (my_task == master_task) then
-                  write(nu_diag,*) ' min and max =', amin, amax
-                  write(nu_diag,*) ' sum =',asum
-                  write(nu_diag,*) ''
+                  write(nu_diag,*) ' min, max, sum =', amin, amax, asum, trim(vname)
                endif
             endif
 
@@ -878,17 +875,15 @@ subroutine write_restart_field(nu,nrec,work,atype,vname,ndim3,diag)
                   amax = global_maxval(work(:,:,n,:),distrb_info)
                   asum = global_sum(work(:,:,n,:), distrb_info, field_loc_center)
                   if (my_task == master_task) then
-                     write(nu_diag,*) ' min and max =', amin, amax
-                     write(nu_diag,*) ' sum =',asum
+                     write(nu_diag,*) ' min, max, sum =', amin, amax, asum, trim(vname)
                   endif
                enddo
             else
                amin = global_minval(work(:,:,1,:),distrb_info)
                amax = global_maxval(work(:,:,1,:),distrb_info)
                asum = global_sum(work(:,:,1,:), distrb_info, field_loc_center)
                if (my_task == master_task) then
-                  write(nu_diag,*) ' min and max =', amin, amax
-                  write(nu_diag,*) ' sum =',asum
+                  write(nu_diag,*) ' min, max, sum =', amin, amax, asum, trim(vname)
                endif
             endif
          endif

cicecore/shared/ice_restart_column.F90

@@ -913,6 +913,7 @@ subroutine write_restart_hbrine()
       use ice_blocks, only: block, get_block
       use ice_domain, only: nblocks, blocks_ice
       use ice_fileunits, only: nu_dump_hbrine
+      use ice_grid, only: tmask
       use ice_state, only: trcrn
       use ice_restart,only: write_restart_field
 
@@ -949,7 +950,8 @@ subroutine write_restart_hbrine()
         do j = jlo, jhi
         do i = ilo, ihi  
            do n = 1, ncat
-              if (first_ice     (i,j,n,iblk)) then
+              ! zero out first_ice over land
+              if (tmask(i,j,iblk) .and. first_ice (i,j,n,iblk)) then
                   first_ice_real(i,j,n,iblk) = c1
               else
                   first_ice_real(i,j,n,iblk) = c0
@@ -983,6 +985,7 @@ subroutine write_restart_bgc()
       use ice_fileunits, only: nu_dump_bgc
       use ice_flux_bgc, only: nit, amm, sil, dmsp, dms, algalN, &
           doc, don, dic, fed, fep, zaeros, hum
+      use ice_grid, only: tmask
       use ice_state, only: trcrn
       use ice_flux, only: sss  
       use ice_restart, only:  write_restart_field
@@ -1058,6 +1061,39 @@ subroutine write_restart_bgc()
 
       diag = .true.
 
+      !-----------------------------------------------------------------
+      ! Zero out tracers over land
+      !-----------------------------------------------------------------
+
+      !$OMP PARALLEL DO PRIVATE(iblk,i,j,ilo,ihi,jlo,jhi,this_block)
+      do iblk = 1, nblocks
+         this_block = get_block(blocks_ice(iblk),iblk)
+         ilo = this_block%ilo
+         ihi = this_block%ihi
+         jlo = this_block%jlo
+         jhi = this_block%jhi
+         do j = jlo, jhi
+         do i = ilo, ihi
+            if (.not. tmask(i,j,iblk)) then
+               if (tr_bgc_N  ) algalN(i,j,:,iblk) = c0
+               if (tr_bgc_C  ) doc   (i,j,:,iblk) = c0
+               if (tr_bgc_C  ) dic   (i,j,:,iblk) = c0
+               if (tr_bgc_Nit) nit   (i,j  ,iblk) = c0
+               if (tr_bgc_Am ) amm   (i,j  ,iblk) = c0
+               if (tr_bgc_Sil) sil   (i,j  ,iblk) = c0
+               if (tr_bgc_hum) hum   (i,j  ,iblk) = c0
+               if (tr_bgc_DMS) dms   (i,j  ,iblk) = c0
+               if (tr_bgc_DMS) dmsp  (i,j  ,iblk) = c0
+               if (tr_bgc_DON) don   (i,j,:,iblk) = c0
+               if (tr_bgc_Fe ) fed   (i,j,:,iblk) = c0
+               if (tr_bgc_Fe ) fep   (i,j,:,iblk) = c0
+               if (solve_zsal) sss   (i,j  ,iblk) = c0
+            endif
+         enddo
+         enddo
+      enddo
+      !$OMP END PARALLEL DO
+
       !-----------------------------------------------------------------
       ! Salinity and extras
       !-----------------------------------------------------------------

configuration/scripts/tests/comparelog.csh

@@ -55,8 +55,8 @@ if (${filearg} == 1) then
       touch ${test_out}
 
       if (${cicefile} == 1) then
-        cat ${base_data} | grep -A 99999999 "total ice area  (km^2)" | grep -e istep1: -e = | grep -iv "min, max, sum" >&! ${base_out}
-        cat ${test_data} | grep -A 99999999 "total ice area  (km^2)" | grep -e istep1: -e = | grep -iv "min, max, sum" >&! ${test_out}
+        cat ${base_data} | grep -A 99999999 "total ice area  (km^2)" | grep -e istep1: -e = | grep -iv "min, max, sum" | grep -iv "init_vert" | grep -iv "ridge_ice"  >&! ${base_out}
+        cat ${test_data} | grep -A 99999999 "total ice area  (km^2)" | grep -e istep1: -e = | grep -iv "min, max, sum" | grep -iv "init_vert" | grep -iv "ridge_ice"  >&! ${test_out}
       else
         sed -n '/RunningUnitTest/,$p' ${base_data} >! ${base_out}
         sed -n '/RunningUnitTest/,$p' ${test_data} >! ${test_out}