ice_dyn_vp: refactor 'puny' treatment for implicit solver

The way 'puny_dyn' is set in the 'cice' driver, and in 'ice_constants'
and 'ice_grid' is kind of hacky. Revert those changes.

Instead, add an 'init_vp' subroutine that is called from the drivers and
that recomputes 'tinyarea' with the new 'puny_vp' value.

Call 'init_evp' from 'init_vp', mimicking what is done ine 'init_eap'.

cicecore/cicedynB/dynamics/ice_dyn_vp.F90

@@ -52,7 +52,8 @@ module ice_dyn_vp
 
       implicit none
       private
-      public :: imp_solver, matvec, arrays_to_vec, vec_to_arrays, precond_diag
+      public :: imp_solver, matvec, arrays_to_vec, vec_to_arrays, precond_diag, &
+                init_vp
 
       ! namelist parameters
 
@@ -87,6 +88,62 @@ module ice_dyn_vp
 
 !=======================================================================
 
+! Initialize parameters and variables needed for the vp dynamics
+! author: Philippe Blain, ECCC
+
+      subroutine init_vp (dt)
+
+      use ice_blocks, only: get_block, block
+      use ice_boundary, only: ice_HaloUpdate
+      use ice_constants, only: c1, &
+          field_loc_center, field_type_scalar
+      use ice_domain, only: nblocks, blocks_ice, halo_info
+      use ice_dyn_shared, only: init_evp
+      use ice_grid, only: tarea, tinyarea
+
+      real (kind=dbl_kind), intent(in) :: &
+         dt      ! time step
+
+      ! local variables
+
+      integer (kind=int_kind) :: &
+         i, j, iblk, &
+         ilo,ihi,jlo,jhi      ! beginning and end of physical domain
+
+      type (block) :: &
+         this_block           ! block information for current block
+
+      real (kind=dbl_kind) :: &
+         puny_vp = 2e-09_dbl_kind      ! special puny value for computing tinyarea
+
+      ! Initialize variables shared with evp
+      call init_evp(dt)
+
+      ! Redefine tinyarea using a different puny value
+
+      !$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
+            tinyarea(i,j,iblk) = puny_vp*tarea(i,j,iblk)
+         enddo
+         enddo
+      enddo                     ! iblk
+      !$OMP END PARALLEL DO
+
+      call ice_HaloUpdate (tinyarea,           halo_info, &
+                           field_loc_center,   field_type_scalar, &
+                           fillValue=c1)
+
+      end subroutine init_vp
+!=======================================================================
+
 ! Viscous-plastic dynamics driver
 !
 #ifdef CICE_IN_NEMO

cicecore/cicedynB/infrastructure/ice_grid.F90

@@ -326,8 +326,7 @@ subroutine init_grid2
       use ice_blocks, only: get_block, block, nx_block, ny_block
       use ice_constants, only: c0, c1, c2, p5, p25, c1p5, &
           field_loc_center, field_loc_NEcorner, &
-          field_type_scalar, field_type_vector, field_type_angle, &
-          puny_dyn
+          field_type_scalar, field_type_vector, field_type_angle
       use ice_domain_size, only: max_blocks
 
       integer (kind=int_kind) :: &
@@ -336,7 +335,7 @@ subroutine init_grid2
 
       real (kind=dbl_kind) :: &
          angle_0, angle_w, angle_s, angle_sw, &
-         pi, pi2
+         pi, pi2, puny
 
       logical (kind=log_kind), dimension(nx_block,ny_block,max_blocks):: &
          out_of_range
@@ -353,7 +352,7 @@ subroutine init_grid2
       ! lat, lon, cell widths, angle, land mask
       !-----------------------------------------------------------------
 
-      call icepack_query_parameters(pi_out=pi, pi2_out=pi2)
+      call icepack_query_parameters(pi_out=pi, pi2_out=pi2, puny_out=puny)
       call icepack_warnings_flush(nu_diag)
       if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
          file=__FILE__, line=__LINE__)
@@ -405,7 +404,7 @@ subroutine init_grid2
             else
                uarear(i,j,iblk) = c0 ! possible on boundaries
             endif
-            tinyarea(i,j,iblk) = puny_dyn*tarea(i,j,iblk)
+            tinyarea(i,j,iblk) = puny*tarea(i,j,iblk)
 
             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))

cicecore/drivers/direct/hadgem3/CICE_InitMod.F90

@@ -72,6 +72,7 @@ subroutine cice_init
       use ice_domain_size, only: ncat, nfsd
       use ice_dyn_eap, only: init_eap, alloc_dyn_eap
       use ice_dyn_shared, only: kdyn, init_evp, basalstress, alloc_dyn_shared
+      use ice_dyn_vp, only: init_vp
       use ice_flux, only: init_coupler_flux, init_history_therm, &
           init_history_dyn, init_flux_atm, init_flux_ocn, alloc_flux
       use ice_forcing, only: init_forcing_ocn, init_forcing_atmo, &
@@ -126,6 +127,8 @@ subroutine cice_init
       if (kdyn == 2) then
          call alloc_dyn_eap     ! allocate dyn_eap arrays
          call init_eap (dt_dyn) ! define eap dynamics parameters, variables
+      else if (kdyn == 3) then
+         call init_vp (dt_dyn)  ! define vp dynamics parameters, variables
       else                      ! for both kdyn = 0 or 1
          call init_evp (dt_dyn) ! define evp dynamics parameters, variables
       endif

cicecore/drivers/mct/cesm1/CICE_InitMod.F90

@@ -72,6 +72,7 @@ subroutine cice_init(mpicom_ice)
       use ice_domain_size, only: ncat, nfsd
       use ice_dyn_eap, only: init_eap, alloc_dyn_eap
       use ice_dyn_shared, only: kdyn, init_evp, alloc_dyn_shared
+      use ice_dyn_vp, only: init_vp
       use ice_flux, only: init_coupler_flux, init_history_therm, &
           init_history_dyn, init_flux_atm, init_flux_ocn, alloc_flux
       use ice_forcing, only: init_forcing_ocn, init_forcing_atmo, &
@@ -128,6 +129,8 @@ subroutine cice_init(mpicom_ice)
       if (kdyn == 2) then
          call alloc_dyn_eap     ! allocate dyn_eap arrays
          call init_eap (dt_dyn) ! define eap dynamics parameters, variables
+      else if (kdyn == 3) then
+         call init_vp (dt_dyn)  ! define vp dynamics parameters, variables
       else                      ! for both kdyn = 0 or 1
          call init_evp (dt_dyn) ! define evp dynamics parameters, variables
       endif

cicecore/drivers/standalone/cice/CICE_InitMod.F90

@@ -67,12 +67,12 @@ subroutine cice_init
       use ice_calendar, only: dt, dt_dyn, time, istep, istep1, write_ic, &
           init_calendar, calendar
       use ice_communicate, only: init_communicate, my_task, master_task
-      use ice_constants, only: ice_init_constants
       use ice_diagnostics, only: init_diags
       use ice_domain, only: init_domain_blocks
       use ice_domain_size, only: ncat, nfsd
       use ice_dyn_eap, only: init_eap, alloc_dyn_eap
       use ice_dyn_shared, only: kdyn, init_evp, alloc_dyn_shared
+      use ice_dyn_vp, only: init_vp
       use ice_flux, only: init_coupler_flux, init_history_therm, &
           init_history_dyn, init_flux_atm, init_flux_ocn, alloc_flux
       use ice_forcing, only: init_forcing_ocn, init_forcing_atmo, &
@@ -94,9 +94,6 @@ subroutine cice_init
 
       logical(kind=log_kind) :: tr_aero, tr_zaero, skl_bgc, z_tracers, &
          tr_iso, tr_fsd, wave_spec
-
-      real (kind=dbl_kind) :: puny
-
       character(len=*), parameter :: subname = '(cice_init)'
 
       call init_communicate     ! initial setup for message passing
@@ -124,16 +121,6 @@ subroutine cice_init
       call alloc_flux           ! allocate flux arrays
       call init_ice_timers      ! initialize all timers
       call ice_timer_start(timer_total)   ! start timing entire run
-      ! By default, the puny value used for computing tinyarea in init_grid2 (puny_dyn)
-      ! is set to a special value for use with the implicit solver (kdyn = 3).
-      ! Thus we reset it back to puny if kdyn .ne. 3
-      if (kdyn /= 3) then
-         call icepack_query_parameters(puny_out=puny)
-         call icepack_warnings_flush(nu_diag)
-         if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
-            file=__FILE__, line=__LINE__)
-         call ice_init_constants(puny_dyn_in=puny)
-      endif
       call init_grid2           ! grid variables
       call init_zbgc            ! vertical biogeochemistry initialization
       call init_calendar        ! initialize some calendar stuff
@@ -142,6 +129,8 @@ subroutine cice_init
       if (kdyn == 2) then
          call alloc_dyn_eap     ! allocate dyn_eap arrays
          call init_eap (dt_dyn) ! define eap dynamics parameters, variables
+      else if (kdyn == 3) then
+         call init_vp (dt_dyn)  ! define vp dynamics parameters, variables
       else                      ! for both kdyn = 0 or 1
          call init_evp (dt_dyn) ! define evp dynamics parameters, variables
       endif

cicecore/shared/ice_constants.F90

@@ -33,9 +33,6 @@ module ice_constants
       real (kind=dbl_kind), public :: &
          shlat  =  30.0_dbl_kind   ,&! artificial masking edge (deg)
          nhlat  = -30.0_dbl_kind     ! artificial masking edge (deg)
-
-      real (kind=dbl_kind), public :: &
-         puny_dyn = 2e-09_dbl_kind      ! special puny value for computing tinyarea for implicit solver
 
       !-----------------------------------------------------------------
       ! numbers used outside the column package
@@ -135,16 +132,15 @@ module ice_constants
 ! subroutine to set the cice constants
 
       subroutine ice_init_constants(   &
-         omega_in, radius_in, spval_dbl_in, spval_in, shlat_in, nhlat_in, puny_dyn_in)
+         omega_in, radius_in, spval_dbl_in, spval_in, shlat_in, nhlat_in)
 
       real (kind=dbl_kind), intent(in), optional :: &
          omega_in     , &   ! angular velocity of earth (rad/sec)
          radius_in    , &   ! earth radius (m)
          spval_dbl_in , &   ! special value (double precision)
          spval_in     , &   ! special value for netCDF output
          shlat_in     , &   ! artificial masking edge (deg)
-         nhlat_in     , &   ! artificial masking edge (deg)
-         puny_dyn_in        ! special puny value for computing tinyarea
+         nhlat_in           ! artificial masking edge (deg)
 
       character(len=*),parameter :: subname='(ice_init_constants)'
 
@@ -154,7 +150,6 @@ subroutine ice_init_constants(   &
       if (present(spval_in)) spval = spval_in
       if (present(shlat_in)) shlat = shlat_in
       if (present(nhlat_in)) nhlat = nhlat_in
-      if (present(puny_dyn_in)) puny_dyn = puny_dyn_in
 
       end subroutine ice_init_constants
 
@@ -163,16 +158,15 @@ end subroutine ice_init_constants
 ! subroutine to set the cice constants
 
       subroutine ice_query_constants(   &
-         omega_out, radius_out, spval_dbl_out, spval_out, shlat_out, nhlat_out, puny_dyn_out)
+         omega_out, radius_out, spval_dbl_out, spval_out, shlat_out, nhlat_out)
 
       real (kind=dbl_kind), intent(out), optional :: &
          omega_out     , &   ! angular velocity of earth (rad/sec)
          radius_out    , &   ! earth radius (m)
          spval_dbl_out , &   ! special value (double precision)
          spval_out     , &   ! special value for netCDF output
          shlat_out     , &   ! artificial masking edge (deg)
-         nhlat_out     , &   ! artificial masking edge (deg)
-         puny_dyn_out        ! special puny value for computing tinyarea
+         nhlat_out           ! artificial masking edge (deg)
 
       character(len=*),parameter :: subname='(ice_query_constants)'
 
@@ -182,7 +176,6 @@ subroutine ice_query_constants(   &
       if (present(spval_out)) spval_out = spval
       if (present(shlat_out)) shlat_out = shlat
       if (present(nhlat_out)) nhlat_out = nhlat
-      if (present(puny_dyn_out)) puny_dyn_out = puny_dyn
 
       end subroutine ice_query_constants