bugfixes for cxp, cyp...

CICE-Consortium · Dec 19, 2023 · 18e7161 · 18e7161
1 parent 6eac1a6
commit 18e7161
Show file tree

Hide file tree

Showing 2 changed files with 47 additions and 34 deletions.
diff --git a/cicecore/cicedyn/dynamics/ice_dyn_shared.F90 b/cicecore/cicedyn/dynamics/ice_dyn_shared.F90
@@ -200,12 +200,17 @@ subroutine alloc_dyn_shared
          stat=ierr)
       if (ierr/=0) call abort_ice(subname//': Out of memory')
 
-      if (grid_ice == 'B' .and. evp_algorithm == "standard_2d") then
+!      if (evp_algorithm == "standard_2d") then
          allocate( &
             cyp      (nx_block,ny_block,max_blocks), & ! 1.5*HTE - 0.5*HTW
             cxp      (nx_block,ny_block,max_blocks), & ! 1.5*HTN - 0.5*HTS
             cym      (nx_block,ny_block,max_blocks), & ! 0.5*HTE - 1.5*HTW
             cxm      (nx_block,ny_block,max_blocks), & ! 0.5*HTN - 1.5*HTS
+            stat=ierr)
+         if (ierr/=0) call abort_ice(subname//': Out of memory')
+
+      if (grid_ice == 'B' .and. evp_algorithm == "standard_2d") then
+         allocate( &
             dxhy     (nx_block,ny_block,max_blocks), & ! 0.5*(HTE - HTW)
             dyhx     (nx_block,ny_block,max_blocks), & ! 0.5*(HTN - HTS)
             stat=ierr)
@@ -360,39 +365,46 @@ subroutine init_dyn_shared (dt)
       !$OMP END PARALLEL DO
 
      if (grid_ice == 'B' .and. evp_algorithm == "standard_2d") then
-      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
-            dxhy(i,j,iblk) = p5*(HTE(i,j,iblk) - HTE(i-1,j,iblk))
-            dyhx(i,j,iblk) = p5*(HTN(i,j,iblk) - HTN(i,j-1,iblk))
-         enddo
-         enddo
-
-         do j = jlo, jhi+1
-         do i = ilo, ihi+1
-            cyp(i,j,iblk) = (c1p5*HTE(i,j,iblk) - p5*HTE(i-1,j,iblk))
-            cxp(i,j,iblk) = (c1p5*HTN(i,j,iblk) - p5*HTN(i,j-1,iblk))
-            ! match order of operations in cyp, cxp for tripole grids
-            cym(i,j,iblk) = -(c1p5*HTE(i-1,j,iblk) - p5*HTE(i,j,iblk))
-            cxm(i,j,iblk) = -(c1p5*HTN(i,j-1,iblk) - p5*HTN(i,j,iblk))
-         enddo
-         enddo
-      enddo
-         call ice_HaloUpdate (dxhy,               halo_info, &
-                           field_loc_center,   field_type_vector, &
-                           fillValue=c1)
-         call ice_HaloUpdate (dyhx,               halo_info, &
-                           field_loc_center,   field_type_vector, &
-                           fillValue=c1)
+       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
+             dxhy(i,j,iblk) = p5*(HTE(i,j,iblk) - HTE(i-1,j,iblk))
+             dyhx(i,j,iblk) = p5*(HTN(i,j,iblk) - HTN(i,j-1,iblk))
+          enddo
+          enddo
+       enddo
+       call ice_HaloUpdate (dxhy,               halo_info, &
+                            field_loc_center,   field_type_vector, &
+                            fillValue=c1)
+       call ice_HaloUpdate (dyhx,               halo_info, &
+                            field_loc_center,   field_type_vector, &
+                            fillValue=c1)
 
      endif
 
+     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+1
+        do i = ilo, ihi+1
+           cyp(i,j,iblk) = (c1p5*HTE(i,j,iblk) - p5*HTE(i-1,j,iblk))
+           cxp(i,j,iblk) = (c1p5*HTN(i,j,iblk) - p5*HTN(i,j-1,iblk))
+            ! match order of operations in cyp, cxp for tripole grids
+           cym(i,j,iblk) = -(c1p5*HTE(i-1,j,iblk) - p5*HTE(i,j,iblk))
+           cxm(i,j,iblk) = -(c1p5*HTN(i,j-1,iblk) - p5*HTN(i,j,iblk))
+       enddo
+       enddo
+     enddo
+
   end subroutine init_dyn_shared
 
 !=======================================================================

diff --git a/cicecore/cicedyn/dynamics/ice_dyn_vp.F90 b/cicecore/cicedyn/dynamics/ice_dyn_vp.F90
@@ -43,13 +43,12 @@ module ice_dyn_vp
           field_type_scalar, field_type_vector
       use ice_constants, only: c0, p027, p055, p111, p166, &
           p222, p25, p333, p5, c1
-      use ice_domain, only: nblocks, distrb_info, halo_info
+      use ice_domain, only: nblocks, distrb_info
       use ice_domain_size, only: max_blocks
       use ice_dyn_shared, only: dyn_prep1, dyn_prep2, dyn_finish, &
           cosw, sinw, fcor_blk, uvel_init, vvel_init, &
           seabed_stress_factor_LKD, seabed_stress_factor_prob, seabed_stress_method, &
-          seabed_stress, Ktens, stack_fields,  unstack_fields, fld2, fld3, fld4, &
-          dxhy, dyhx, cxp, cyp, cxm, cym
+          seabed_stress, Ktens, stack_fields,  unstack_fields, fld2, fld3, fld4
       use ice_fileunits, only: nu_diag
       use ice_flux, only: fmU
       use ice_global_reductions, only: global_sum
@@ -2744,6 +2743,7 @@ subroutine fgmres (zetax2   , etax2   , &
       use ice_boundary, only: ice_HaloUpdate
       use ice_domain, only: maskhalo_dyn, halo_info
       use ice_timers, only: ice_timer_start, ice_timer_stop, timer_bound
+      use ice_dyn_shared, only: dxhy, dyhx, cxp, cyp, cxm, cym
 
       real (kind=dbl_kind), dimension(nx_block,ny_block,max_blocks,4), intent(in) :: &
          zetax2   , & ! zetax2 = 2*zeta (bulk viscosity)
@@ -3145,6 +3145,7 @@ subroutine pgmres (zetax2   , etax2   , &
       use ice_boundary, only: ice_HaloUpdate
       use ice_domain, only: maskhalo_dyn, halo_info
       use ice_timers, only: ice_timer_start, ice_timer_stop, timer_bound
+      use ice_dyn_shared, only: dyhx, dxhy, cxp, cyp, cxm, cym
 
       real (kind=dbl_kind), dimension(nx_block,ny_block,max_blocks,4), intent(in) :: &
          zetax2   , & ! zetax2 = 2*zeta (bulk viscosity)