From 2872cda347ef3a02bb4a61404b68246d8fa02d5e Mon Sep 17 00:00:00 2001 From: Philippe Blain Date: Mon, 10 Jun 2019 13:34:12 -0400 Subject: [PATCH] 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 | 59 ++++++++++++++++++- cicecore/cicedynB/infrastructure/ice_grid.F90 | 9 ++- .../drivers/direct/hadgem3/CICE_InitMod.F90 | 3 + cicecore/drivers/mct/cesm1/CICE_InitMod.F90 | 3 + .../drivers/standalone/cice/CICE_InitMod.F90 | 17 +----- cicecore/shared/ice_constants.F90 | 15 ++--- 6 files changed, 75 insertions(+), 31 deletions(-) diff --git a/cicecore/cicedynB/dynamics/ice_dyn_vp.F90 b/cicecore/cicedynB/dynamics/ice_dyn_vp.F90 index f07fa1ecc..3bab1290e 100644 --- a/cicecore/cicedynB/dynamics/ice_dyn_vp.F90 +++ b/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 diff --git a/cicecore/cicedynB/infrastructure/ice_grid.F90 b/cicecore/cicedynB/infrastructure/ice_grid.F90 index 21230b8b8..28ed56492 100644 --- a/cicecore/cicedynB/infrastructure/ice_grid.F90 +++ b/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)) diff --git a/cicecore/drivers/direct/hadgem3/CICE_InitMod.F90 b/cicecore/drivers/direct/hadgem3/CICE_InitMod.F90 index b208bcbef..3494e0f2d 100644 --- a/cicecore/drivers/direct/hadgem3/CICE_InitMod.F90 +++ b/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 diff --git a/cicecore/drivers/mct/cesm1/CICE_InitMod.F90 b/cicecore/drivers/mct/cesm1/CICE_InitMod.F90 index b72745e30..f14354af3 100644 --- a/cicecore/drivers/mct/cesm1/CICE_InitMod.F90 +++ b/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 diff --git a/cicecore/drivers/standalone/cice/CICE_InitMod.F90 b/cicecore/drivers/standalone/cice/CICE_InitMod.F90 index 1dc98eea5..3dd426765 100644 --- a/cicecore/drivers/standalone/cice/CICE_InitMod.F90 +++ b/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 diff --git a/cicecore/shared/ice_constants.F90 b/cicecore/shared/ice_constants.F90 index b0370b0b1..c49732e35 100644 --- a/cicecore/shared/ice_constants.F90 +++ b/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,7 +132,7 @@ 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) @@ -143,8 +140,7 @@ subroutine ice_init_constants( & 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,7 +158,7 @@ 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) @@ -171,8 +166,7 @@ subroutine ice_query_constants( & 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