New 1d evp solver (CICE-Consortium#895)

* New 1d evp solver

* Small changes incl timer names and inclued private/publice in ice_dyn_core1d

* fixed bug on gnu debug

* moved halo update to evp1d, added deallocation, fixed bug

* fixed deallocation dyn_evp1d

* bugfix deallocate

* Remove gather strintx and strinty

* removed 4 test with evp1d and c/cd grid

* Update of evp1d implementation

- Rename halo_HTE_HTN to global_ext_halo and move into ice_grid.F90
- Generalize global_ext_halo to work with any nghost size (was hardcoded for nghost=1)
- Remove argument from dyn_evp1d_init, change to "use" of global grid variables
- rename pgl_global_ext to save_ghte_ghtn
- Update allocation of G_HTE, G_HTN
- Add dealloc_grid to deallocate G_HTE and G_HTN at end of initialization
- Add calls to dealloc_grid to all CICE_InitMod.F90 subroutines
- Make dimension of evp1d arguments implicit size more consistently
- Clean up indentation and formatting a bit

* Clean up trailing blanks

* resolved name conflicts

* 1d grid var name change

---------

Co-authored-by: apcraig <anthony.p.craig@gmail.com>

cicecore/cicedyn/dynamics/ice_dyn_shared.F90

@@ -1451,10 +1451,10 @@ subroutine seabed_stress_factor_prob (nx_block, ny_block,          &
 
       character(len=*), parameter :: subname = '(seabed_stress_factor_prob)'
 
-      call icepack_query_parameters(rhow_out=rhow, rhoi_out=rhoi)
-      call icepack_query_parameters(gravit_out=gravit)
-      call icepack_query_parameters(pi_out=pi)
-      call icepack_query_parameters(puny_out=puny)
+      call icepack_query_parameters(rhow_out=rhow, rhoi_out=rhoi,gravit_out=gravit,pi_out=pi,puny_out=puny)
+      call icepack_warnings_flush(nu_diag)
+      if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
+         file=__FILE__, line=__LINE__)
 
       Tbt=c0
 
@@ -2302,16 +2302,15 @@ end subroutine strain_rates_U
 ! by combining tensile strength and a parameterization for grounded ridges.
 ! J. Geophys. Res. Oceans, 121, 7354-7368.
 
-      subroutine visc_replpress(strength, DminArea, Delta, &
-                                zetax2, etax2, rep_prs, capping)
+    subroutine visc_replpress(strength, DminArea, Delta, &
+                                zetax2, etax2, rep_prs)
 
       real (kind=dbl_kind), intent(in)::  &
          strength, & !
          DminArea    !
 
       real (kind=dbl_kind), intent(in)::  &
-         Delta   , & !
-         capping     !
+         Delta
 
       real (kind=dbl_kind), intent(out):: &
          zetax2  , & ! bulk viscosity

cicecore/cicedyn/dynamics/ice_dyn_vp.F90

@@ -1221,20 +1221,16 @@ subroutine calc_zeta_dPr (nx_block, ny_block, &
 
          call visc_replpress (strength(i,j)  , DminTarea(i,j)  , &
                               Deltane        , zetax2   (i,j,1), &
-                              etax2   (i,j,1), rep_prs  (i,j,1), &
-                              capping)
+                              etax2   (i,j,1), rep_prs  (i,j,1))
          call visc_replpress (strength(i,j)  , DminTarea(i,j)  , &
                               Deltanw        , zetax2   (i,j,2), &
-                              etax2   (i,j,2), rep_prs  (i,j,2), &
-                              capping)
+                              etax2   (i,j,2), rep_prs  (i,j,2))
          call visc_replpress (strength(i,j)  , DminTarea(i,j)  , &
                               Deltasw        , zetax2   (i,j,3), &
-                              etax2   (i,j,3), rep_prs  (i,j,3), &
-                              capping)
+                              etax2   (i,j,3), rep_prs  (i,j,3))
          call visc_replpress (strength(i,j)  , DminTarea(i,j)  , &
                               Deltase        , zetax2   (i,j,4), &
-                              etax2   (i,j,4), rep_prs  (i,j,4), &
-                              capping)
+                              etax2   (i,j,4), rep_prs  (i,j,4))
 
       !-----------------------------------------------------------------
       ! the stresses                            ! kg/s^2

cicecore/cicedyn/dynamics/ice_transport_remap.F90

@@ -647,11 +647,12 @@ subroutine horizontal_remap (dt,             ntrace,   &
       endif  ! nghost
 
       ! tcraig, this OMP loop sometimes fails with cce/14.0.3, compiler bug??
-      !$OMP PARALLEL DO PRIVATE(iblk,i,j,ilo,ihi,jlo,jhi,this_block,n,  &
-      !$OMP                     edgearea_e,edgearea_n,edge,iflux,jflux, &
-      !$OMP                     xp,yp,indxing,indxjng,mflxe,mflxn, &
-      !$OMP                     mtflxe,mtflxn,triarea,istop,jstop,l_stop) &
-      !$OMP SCHEDULE(runtime)
+      ! TILL I can trigger the same with ifort (IFORT) 18.0.0 20170811
+!TILL      !$OMP PARALLEL DO PRIVATE(iblk,i,j,ilo,ihi,jlo,jhi,this_block,n,  &
+!TILL      !$OMP                     edgearea_e,edgearea_n,edge,iflux,jflux, &
+!TILL      !$OMP                     xp,yp,indxing,indxjng,mflxe,mflxn, &
+!TILL      !$OMP                     mtflxe,mtflxn,triarea,istop,jstop,l_stop) &
+!TILL      !$OMP SCHEDULE(runtime)
       do iblk = 1, nblocks
 
          l_stop = .false.
@@ -865,7 +866,7 @@ subroutine horizontal_remap (dt,             ntrace,   &
          enddo                  ! n
 
       enddo                     ! iblk
-      !$OMP END PARALLEL DO
+!TILL      !$OMP END PARALLEL DO
 
       end subroutine horizontal_remap
 

cicecore/cicedyn/general/ice_init.F90

@@ -105,7 +105,7 @@ subroutine input_data
                           grid_ocn, grid_ocn_thrm, grid_ocn_dynu, grid_ocn_dynv, &
                           grid_atm, grid_atm_thrm, grid_atm_dynu, grid_atm_dynv, &
                           dxrect, dyrect, dxscale, dyscale, scale_dxdy, &
-                          lonrefrect, latrefrect, pgl_global_ext
+                          lonrefrect, latrefrect, save_ghte_ghtn
       use ice_dyn_shared, only: ndte, kdyn, revised_evp, yield_curve, &
                                 evp_algorithm, visc_method,     &
                                 seabed_stress, seabed_stress_method,  &
@@ -375,7 +375,7 @@ subroutine input_data
       ndte = 120         ! subcycles per dynamics timestep:  ndte=dt_dyn/dte
       evp_algorithm = 'standard_2d'  ! EVP kernel (standard_2d=standard cice evp; shared_mem_1d=1d shared memory and no mpi
       elasticDamp = 0.36_dbl_kind    ! coefficient for calculating the parameter E
-      pgl_global_ext = .false.       ! if true, init primary grid lengths (global ext.)
+      save_ghte_ghtn = .false.       ! if true, save global hte and htn (global ext.)
       brlx   = 300.0_dbl_kind ! revised_evp values. Otherwise overwritten in ice_dyn_shared
       arlx   = 300.0_dbl_kind ! revised_evp values. Otherwise overwritten in ice_dyn_shared
       revised_evp = .false.   ! if true, use revised procedure for evp dynamics
@@ -963,7 +963,6 @@ subroutine input_data
       call broadcast_scalar(ndte,                 master_task)
       call broadcast_scalar(evp_algorithm,        master_task)
       call broadcast_scalar(elasticDamp,          master_task)
-      call broadcast_scalar(pgl_global_ext,       master_task)
       call broadcast_scalar(brlx,                 master_task)
       call broadcast_scalar(arlx,                 master_task)
       call broadcast_scalar(revised_evp,          master_task)
@@ -1258,6 +1257,10 @@ subroutine input_data
          abort_list = trim(abort_list)//":5"
       endif
 
+      if (kdyn == 1 .and. evp_algorithm == 'shared_mem_1d') then
+         save_ghte_ghtn = .true.
+      endif
+
       if (kdyn == 2 .and. revised_evp) then
          if (my_task == master_task) then
             write(nu_diag,*) subname//' WARNING: revised_evp = T with EAP dynamics'
@@ -1296,10 +1299,10 @@ subroutine input_data
       endif
 
       if (grid_ice == 'C' .or. grid_ice == 'CD') then
-         if (kdyn > 1) then
+         if (kdyn > 1 .or. (kdyn == 1 .and. evp_algorithm /= 'standard_2d')) then
             if (my_task == master_task) then
-               write(nu_diag,*) subname//' ERROR: grid_ice = C | CD only supported with kdyn<=1 (evp or off)'
-               write(nu_diag,*) subname//' ERROR: kdyn and grid_ice inconsistency'
+              write(nu_diag,*) subname//' ERROR: grid_ice = C | CD only supported with kdyn=1 and evp_algorithm=standard_2d'
+              write(nu_diag,*) subname//' ERROR: kdyn and/or evp_algorithm and grid_ice inconsistency'
             endif
             abort_list = trim(abort_list)//":46"
          endif
@@ -1312,6 +1315,15 @@ subroutine input_data
          endif
       endif
 
+      if (evp_algorithm == 'shared_mem_1d' .and. &
+          grid_type     == 'tripole') then
+          if (my_task == master_task) then
+              write(nu_diag,*) subname//' ERROR: evp_algorithm=shared_mem_1d is not tested for gridtype=tripole'
+              write(nu_diag,*) subname//' ERROR: change evp_algorithm to standard_2d'
+          endif
+          abort_list = trim(abort_list)//":49"
+      endif
+
       capping = -9.99e30
       if (kdyn == 1 .or. kdyn == 3) then
          if (capping_method == 'max') then
@@ -1833,7 +1845,6 @@ subroutine input_data
                   tmpstr2 = ' : standard 2d EVP solver'
                elseif (evp_algorithm == 'shared_mem_1d') then
                   tmpstr2 = ' : vectorized 1d EVP solver'
-                  pgl_global_ext = .true.
                else
                   tmpstr2 = ' : unknown value'
                endif

cicecore/cicedyn/infrastructure/comm/mpi/ice_boundary.F90

@@ -113,8 +113,7 @@ module ice_boundary
              ice_HaloUpdate, &
              ice_HaloUpdate_stress, &
              ice_HaloExtrapolate, &
-             ice_HaloDestroy, &
-             primary_grid_lengths_global_ext
+             ice_HaloDestroy
 
    interface ice_HaloUpdate  ! generic interface
       module procedure ice_HaloUpdate2DR8, &
@@ -7164,134 +7163,8 @@ subroutine ice_HaloDestroy(halo)
       call abort_ice(subname,' ERROR: deallocating')
       return
    endif
-end subroutine ice_HaloDestroy
-
-!***********************************************************************
-
- subroutine primary_grid_lengths_global_ext( &
-   ARRAY_O, ARRAY_I, ew_boundary_type, ns_boundary_type)
-
-!  This subroutine adds ghost cells to global primary grid lengths array
-!  ARRAY_I and outputs result to array ARRAY_O
-
-   use ice_constants, only: c0
-   use ice_domain_size, only: nx_global, ny_global
-
-   real (kind=dbl_kind), dimension(:,:), intent(in) :: &
-      ARRAY_I
-
-   character (*), intent(in) :: &
-      ew_boundary_type, ns_boundary_type
-
-   real (kind=dbl_kind), dimension(:,:), intent(out) :: &
-      ARRAY_O
-
-!-----------------------------------------------------------------------
-!
-!  local variables
-!
-!-----------------------------------------------------------------------
-
-   integer (kind=int_kind) :: &
-      ii, io, ji, jo
-
-   character(len=*), parameter :: &
-      subname = '(primary_grid_lengths_global_ext)'
-
-!-----------------------------------------------------------------------
-!
-!  add ghost cells to global primary grid lengths array
-!
-!-----------------------------------------------------------------------
-
-   if (trim(ns_boundary_type) == 'tripole' .or. &
-       trim(ns_boundary_type) == 'tripoleT') then
-      call abort_ice(subname//' ERROR: '//ns_boundary_type &
-         //' boundary type not implemented for configuration')
-   endif
-
-   do jo = 1,ny_global+2*nghost
-      ji = -nghost + jo
-
-      !*** Southern ghost cells
-
-      if (ji < 1) then
-         select case (trim(ns_boundary_type))
-         case ('cyclic')
-            ji = ji + ny_global
-         case ('open')
-            ji = nghost - jo + 1
-         case ('closed')
-            ji = 0
-         case default
-            call abort_ice( &
-               subname//' ERROR: unknown north-south boundary type')
-         end select
-      endif
 
-      !*** Northern ghost cells
-
-      if (ji > ny_global) then
-         select case (trim(ns_boundary_type))
-         case ('cyclic')
-            ji = ji - ny_global
-         case ('open')
-            ji = 2 * ny_global - ji + 1
-         case ('closed')
-            ji = 0
-         case default
-            call abort_ice( &
-               subname//' ERROR: unknown north-south boundary type')
-         end select
-      endif
-
-      do io = 1,nx_global+2*nghost
-         ii = -nghost + io
-
-         !*** Western ghost cells
-
-         if (ii < 1) then
-            select case (trim(ew_boundary_type))
-            case ('cyclic')
-               ii = ii + nx_global
-            case ('open')
-               ii = nghost - io + 1
-            case ('closed')
-               ii = 0
-            case default
-               call abort_ice( &
-                  subname//' ERROR: unknown east-west boundary type')
-            end select
-         endif
-
-         !*** Eastern ghost cells
-
-         if (ii > nx_global) then
-            select case (trim(ew_boundary_type))
-            case ('cyclic')
-               ii = ii - nx_global
-            case ('open')
-               ii = 2 * nx_global - ii + 1
-            case ('closed')
-               ii = 0
-            case default
-               call abort_ice( &
-                  subname//' ERROR: unknown east-west boundary type')
-            end select
-         endif
-
-         if (ii == 0 .or. ji == 0) then
-            ARRAY_O(io, jo) = c0
-         else
-            ARRAY_O(io, jo) = ARRAY_I(ii, ji)
-         endif
-
-      enddo
-   enddo
-
-!-----------------------------------------------------------------------
-
- end subroutine primary_grid_lengths_global_ext
+end subroutine ice_HaloDestroy
 
 !***********************************************************************
 

cicecore/cicedyn/infrastructure/comm/mpi/ice_timers.F90

@@ -65,8 +65,8 @@ module ice_timers
       timer_bundbound,        &! boundary updates bundling
       timer_bgc,              &! biogeochemistry
       timer_forcing,          &! forcing
-      timer_evp_1d,           &! timer only loop
-      timer_evp_2d,           &! timer including conversion 1d/2d
+      timer_evp1dcore,        &! timer only loop
+      timer_evp,              &! timer including conversion 1d/2d
       timer_updstate           ! update state
 !      timer_updstate,         &! update state
 !      timer_tmp1,             &! for temporary timings
@@ -177,34 +177,34 @@ subroutine init_ice_timers
       nullify(all_timers(n)%block_accum_time)
    end do
 
-   call get_ice_timer(timer_total,    'Total',    nblocks,distrb_info%nprocs)
-   call get_ice_timer(timer_step,     'TimeLoop', nblocks,distrb_info%nprocs)
-   call get_ice_timer(timer_dynamics, 'Dynamics', nblocks,distrb_info%nprocs)
-   call get_ice_timer(timer_advect,   'Advection',nblocks,distrb_info%nprocs)
-   call get_ice_timer(timer_column,   'Column',   nblocks,distrb_info%nprocs)
-   call get_ice_timer(timer_thermo,   'Thermo',   nblocks,distrb_info%nprocs)
-   call get_ice_timer(timer_sw,       'Shortwave',nblocks,distrb_info%nprocs)
+   call get_ice_timer(timer_total     , 'Total'     ,nblocks,distrb_info%nprocs)
+   call get_ice_timer(timer_step      , 'TimeLoop'  ,nblocks,distrb_info%nprocs)
+   call get_ice_timer(timer_dynamics  , 'Dynamics'  ,nblocks,distrb_info%nprocs)
+   call get_ice_timer(timer_advect    , 'Advection' ,nblocks,distrb_info%nprocs)
+   call get_ice_timer(timer_column    , 'Column'    ,nblocks,distrb_info%nprocs)
+   call get_ice_timer(timer_thermo    , 'Thermo'    ,nblocks,distrb_info%nprocs)
+   call get_ice_timer(timer_sw        , 'Shortwave' ,nblocks,distrb_info%nprocs)
 !   call get_ice_timer(timer_ponds,    'Meltponds',nblocks,distrb_info%nprocs)
-   call get_ice_timer(timer_ridge,    'Ridging',  nblocks,distrb_info%nprocs)
+   call get_ice_timer(timer_ridge     , 'Ridging'   ,nblocks,distrb_info%nprocs)
 !   call get_ice_timer(timer_catconv,  'Cat Conv', nblocks,distrb_info%nprocs)
-   call get_ice_timer(timer_fsd,      'FloeSize', nblocks,distrb_info%nprocs)
-   call get_ice_timer(timer_couple,   'Coupling', nblocks,distrb_info%nprocs)
-   call get_ice_timer(timer_readwrite,'ReadWrite',nblocks,distrb_info%nprocs)
-   call get_ice_timer(timer_diags,    'Diags    ',nblocks,distrb_info%nprocs)
-   call get_ice_timer(timer_hist,     'History  ',nblocks,distrb_info%nprocs)
-   call get_ice_timer(timer_bound,    'Bound',    nblocks,distrb_info%nprocs)
-   call get_ice_timer(timer_bundbound,'Bundbound',nblocks,distrb_info%nprocs)
-   call get_ice_timer(timer_bgc,      'BGC',      nblocks,distrb_info%nprocs)
-   call get_ice_timer(timer_forcing,  'Forcing',  nblocks,distrb_info%nprocs)
+   call get_ice_timer(timer_fsd       , 'FloeSize'  ,nblocks,distrb_info%nprocs)
+   call get_ice_timer(timer_couple    , 'Coupling'  ,nblocks,distrb_info%nprocs)
+   call get_ice_timer(timer_readwrite , 'ReadWrite' ,nblocks,distrb_info%nprocs)
+   call get_ice_timer(timer_diags     , 'Diags    ' ,nblocks,distrb_info%nprocs)
+   call get_ice_timer(timer_hist      , 'History  ' ,nblocks,distrb_info%nprocs)
+   call get_ice_timer(timer_bound     , 'Bound'     ,nblocks,distrb_info%nprocs)
+   call get_ice_timer(timer_bundbound , 'Bundbound' ,nblocks,distrb_info%nprocs)
+   call get_ice_timer(timer_bgc       , 'BGC'       ,nblocks,distrb_info%nprocs)
+   call get_ice_timer(timer_forcing   , 'Forcing'   ,nblocks,distrb_info%nprocs)
 #if (defined CESMCOUPLED)
-   call get_ice_timer(timer_cplrecv,  'Cpl-recv', nblocks,distrb_info%nprocs)
-   call get_ice_timer(timer_rcvsnd,   'Rcv->Snd', nblocks,distrb_info%nprocs)
-   call get_ice_timer(timer_cplsend,  'Cpl-Send', nblocks,distrb_info%nprocs)
-   call get_ice_timer(timer_sndrcv,   'Snd->Rcv', nblocks,distrb_info%nprocs)
+   call get_ice_timer(timer_cplrecv   , 'Cpl-recv'  ,nblocks,distrb_info%nprocs)
+   call get_ice_timer(timer_rcvsnd    , 'Rcv->Snd'  ,nblocks,distrb_info%nprocs)
+   call get_ice_timer(timer_cplsend   , 'Cpl-Send'  ,nblocks,distrb_info%nprocs)
+   call get_ice_timer(timer_sndrcv    , 'Snd->Rcv'  ,nblocks,distrb_info%nprocs)
 #endif
-   call get_ice_timer(timer_evp_1d,   '1d-evp',   nblocks,distrb_info%nprocs)
-   call get_ice_timer(timer_evp_2d,   '2d-evp',   nblocks,distrb_info%nprocs)
-   call get_ice_timer(timer_updstate, 'UpdState', nblocks,distrb_info%nprocs)
+   call get_ice_timer(timer_evp1dcore , 'evp1dcore' ,nblocks,distrb_info%nprocs)
+   call get_ice_timer(timer_evp       , 'evp'       ,nblocks,distrb_info%nprocs)
+   call get_ice_timer(timer_updstate  , 'UpdState'  ,nblocks,distrb_info%nprocs)
 !   call get_ice_timer(timer_tmp1,     'tmp1',     nblocks,distrb_info%nprocs)
 !   call get_ice_timer(timer_tmp2,     'tmp2',     nblocks,distrb_info%nprocs)
 !   call get_ice_timer(timer_tmp3,     'tmp3',     nblocks,distrb_info%nprocs)