Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Bug fixed 1d evp #568

Merged
merged 41 commits into from
Aug 5, 2021
Merged

Bug fixed 1d evp #568

Show file tree
Hide file tree
Changes from 3 commits
Commits
Show all changes
41 commits
Select commit Hold shift + click to select a range
0ba9ddf
Initialization of global HTE and HTN arrays with ghost cells on maste…
srethmeier Mar 1, 2021
f5c1c1f
Reverted configuration/scripts/cice.batch.csh to resolve bad auto mer…
srethmeier Mar 2, 2021
17d48d8
Merge pull request #13 from TillRasmussen/evp-1d
TillRasmussen Mar 2, 2021
b4bf34b
Spell error corrected (https://github.com/CICE-Consortium/CICE/pull/5…
srethmeier Mar 12, 2021
4b702b8
Redundant initializations to default values removed (https://github.c…
srethmeier Mar 12, 2021
a711bdc
Subroutine primary_grid_lengths_global_ext moved from ice_grid module…
srethmeier Mar 13, 2021
d669c8d
Merge pull request #14 from TillRasmussen/evp-1d
TillRasmussen Mar 14, 2021
aabb12f
Added 1d solver to be fully integrated. Removed kevp=102 bypass.
TillRasmussen Mar 16, 2021
1194e0c
General cleanup of ice_dyn_evp_1d module (https://github.com/CICE-Con…
srethmeier Mar 21, 2021
6f11797
Modified documentation in order to comply with modifications of evp 1d
Mar 23, 2021
462e7f0
Uncommenting of numainit parallelization removed and domp_init moved …
srethmeier Mar 23, 2021
8b02825
Merge pull request #15 from TillRasmussen/evp-1d
TillRasmussen Mar 23, 2021
bc79443
renamed set_nml.evp_algorithm to set_nml.evp1d
Mar 23, 2021
d0d3523
Merge branch 'master' of github.com:TillRasmussen/CICE
Mar 23, 2021
5bb4cba
rm set_nml.kevp102
Mar 23, 2021
37f42d9
Merge branch 'master' of https://github.com/cice-consortium/CICE into…
srethmeier Mar 24, 2021
0a5f938
Merge branch 'evp-1d'
srethmeier Mar 24, 2021
f624275
Merge remote-tracking branch 'upstream/master'
TillRasmussen Mar 25, 2021
09be327
Naming of variable se is not logical. Changed to sse which is the same
TillRasmussen Apr 6, 2021
bdc65bd
Alignment of NUMA-aware array splitting between threads
srethmeier Apr 7, 2021
3b8f6c7
Initialization of global HTE and HTN arrays with ghost cells on maste…
srethmeier Mar 1, 2021
ec162b2
Reverted configuration/scripts/cice.batch.csh to resolve bad auto mer…
srethmeier Mar 2, 2021
cbc4832
Spell error corrected (https://github.com/CICE-Consortium/CICE/pull/5…
srethmeier Mar 12, 2021
c3663ec
Redundant initializations to default values removed (https://github.c…
srethmeier Mar 12, 2021
bbe2e13
Subroutine primary_grid_lengths_global_ext moved from ice_grid module…
srethmeier Mar 13, 2021
15a4e9c
Added 1d solver to be fully integrated. Removed kevp=102 bypass.
TillRasmussen Mar 16, 2021
8e7970c
General cleanup of ice_dyn_evp_1d module (https://github.com/CICE-Con…
srethmeier Mar 21, 2021
6cccdaa
Modified documentation in order to comply with modifications of evp 1d
Mar 23, 2021
462ed9e
Uncommenting of numainit parallelization removed and domp_init moved …
srethmeier Mar 23, 2021
6b97cf0
renamed set_nml.evp_algorithm to set_nml.evp1d
Mar 23, 2021
08debdc
rm set_nml.kevp102
Mar 23, 2021
2f849ed
Naming of variable se is not logical. Changed to sse which is the same
TillRasmussen Apr 6, 2021
e3b00e4
Alignment of NUMA-aware array splitting between threads
srethmeier Apr 7, 2021
e2f49f9
Merge branch 'master' of github.com:TillRasmussen/CICE
TillRasmussen Apr 14, 2021
e4bedb4
Fix bad scaling of calculation of tinyarea. Different order of multip…
TillRasmussen Jul 11, 2021
fb2b3ed
Merge branch 'master' into master
TillRasmussen Jul 11, 2021
3eb47b2
Merge branch 'master' into master
srethmeier Jul 16, 2021
1262455
Fix evp 1d masking issues
apcraig Jul 16, 2021
064bc47
Merge pull request #17 from apcraig/evp1da
srethmeier Jul 17, 2021
1869e93
Some minor revisions in the documentation as suggested by @apcraig (h…
srethmeier Jul 17, 2021
92c49bb
Merge branch 'master' of https://github.com/TillRasmussen/CICE
srethmeier Jul 17, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
3 changes: 1 addition & 2 deletions cicecore/cicedynB/dynamics/ice_dyn_evp.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -88,7 +88,7 @@ subroutine evp (dt)
stress12_1, stress12_2, stress12_3, stress12_4
use ice_grid, only: tmask, umask, dxt, dyt, dxhy, dyhx, cxp, cyp, cxm, cym, &
tarear, uarear, tinyarea, to_ugrid, t2ugrid_vector, u2tgrid_vector, &
grid_type, HTE, HTN
grid_type
use ice_state, only: aice, vice, vsno, uvel, vvel, divu, shear, &
aice_init, aice0, aicen, vicen, strength
use ice_timers, only: timer_dynamics, timer_bound, &
Expand Down Expand Up @@ -366,7 +366,6 @@ subroutine evp (dt)
endif
call ice_dyn_evp_1d_copyin( &
nx_block,ny_block,nblocks,nx_global+2*nghost,ny_global+2*nghost, &
HTE,HTN, &
!v1 dxhy,dyhx,cyp,cxp,cym,cxm,tinyarea, &
!v1 waterx,watery, &
icetmask, iceumask, &
Expand Down
93 changes: 45 additions & 48 deletions cicecore/cicedynB/dynamics/ice_dyn_evp_1d.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1135,7 +1135,6 @@ end subroutine dealloc1d
!----------------------------------------------------------------------------

subroutine evp_copyin_v2(nx,ny,nblk,nx_glob,ny_glob, &
I_HTE,I_HTN, &
!v1 I_dxhy,I_dyhx,I_cyp,I_cxp,I_cym,I_cxm,I_tinyarea, &
!v1 I_waterx,I_watery, &
I_icetmask,I_iceumask, &
Expand All @@ -1149,14 +1148,15 @@ subroutine evp_copyin_v2(nx,ny,nblk,nx_glob,ny_glob,
use ice_gather_scatter, only: gather_global_ext
use ice_domain, only: distrb_info
use ice_communicate, only: my_task, master_task
use ice_grid, only: G_HTE, G_HTN
use ice_constants, only: c0

implicit none

integer(int_kind), intent(in) :: nx, ny, nblk, nx_glob, ny_glob
integer (kind=int_kind),dimension (nx,ny,nblk), intent(in) :: I_icetmask
logical (kind=log_kind),dimension (nx,ny,nblk), intent(in) :: I_iceumask
real (kind=dbl_kind), dimension(nx,ny,nblk), intent(in) :: &
I_HTE,I_HTN, &
!v1 I_dxhy,I_dyhx,I_cyp,I_cxp,I_cym,I_cxm,I_tinyarea, &
!v1 I_waterx,I_watery, &
I_cdn_ocn,I_aiu,I_uocn,I_vocn,I_forcex,I_forcey,I_Tbu, &
Expand All @@ -1171,7 +1171,6 @@ subroutine evp_copyin_v2(nx,ny,nblk,nx_glob,ny_glob,
integer (kind=int_kind),dimension (nx_glob,ny_glob) :: G_icetmask
logical (kind=log_kind),dimension (nx_glob,ny_glob) :: G_iceumask
real (kind=dbl_kind), dimension(nx_glob,ny_glob) :: &
G_HTE,G_HTN, &
!v1 G_dxhy,G_dyhx,G_cyp,G_cxp,G_cym,G_cxm,G_tinyarea, &
!v1 G_waterx,G_watery, &
G_cdn_ocn,G_aiu,G_uocn,G_vocn,G_forcex,G_forcey,G_Tbu, &
Expand All @@ -1186,51 +1185,49 @@ subroutine evp_copyin_v2(nx,ny,nblk,nx_glob,ny_glob,
!---------------------------------------
!-- Gather data into one single block --

call gather_global_ext(G_icetmask, I_icetmask, master_task, distrb_info)
call gather_global_ext(G_iceumask, I_iceumask, master_task, distrb_info)
call gather_global_ext(G_HTE, I_HTE, master_task, distrb_info)
call gather_global_ext(G_HTN, I_HTN, master_task, distrb_info)
!v1 call gather_global_ext(G_dxhy, I_dxhy, master_task, distrb_info)
!v1 call gather_global_ext(G_dyhx, I_dyhx, master_task, distrb_info)
!v1 call gather_global_ext(G_cyp, I_cyp, master_task, distrb_info)
!v1 call gather_global_ext(G_cxp, I_cxp, master_task, distrb_info)
!v1 call gather_global_ext(G_cym, I_cym, master_task, distrb_info)
!v1 call gather_global_ext(G_cxm, I_cxm, master_task, distrb_info)
!v1 call gather_global_ext(G_tinyarea, I_tinyarea, master_task, distrb_info)
!v1 call gather_global_ext(G_waterx, I_waterx, master_task, distrb_info)
!v1 call gather_global_ext(G_watery, I_watery, master_task, distrb_info)
call gather_global_ext(G_cdn_ocn, I_cdn_ocn, master_task, distrb_info)
call gather_global_ext(G_aiu, I_aiu, master_task, distrb_info)
call gather_global_ext(G_uocn, I_uocn, master_task, distrb_info)
call gather_global_ext(G_vocn, I_vocn, master_task, distrb_info)
call gather_global_ext(G_forcex, I_forcex, master_task, distrb_info)
call gather_global_ext(G_forcey, I_forcey, master_task, distrb_info)
call gather_global_ext(G_Tbu, I_Tbu, master_task, distrb_info)
call gather_global_ext(G_umassdti, I_umassdti, master_task, distrb_info)
call gather_global_ext(G_fm, I_fm, master_task, distrb_info)
call gather_global_ext(G_uarear, I_uarear, master_task, distrb_info)
call gather_global_ext(G_tarear, I_tarear, master_task, distrb_info)
call gather_global_ext(G_strintx, I_strintx, master_task, distrb_info)
call gather_global_ext(G_strinty, I_strinty, master_task, distrb_info)
call gather_global_ext(G_uvel_init, I_uvel_init, master_task, distrb_info)
call gather_global_ext(G_vvel_init, I_vvel_init, master_task, distrb_info)
call gather_global_ext(G_strength, I_strength, master_task, distrb_info)
call gather_global_ext(G_uvel, I_uvel, master_task, distrb_info)
call gather_global_ext(G_vvel, I_vvel, master_task, distrb_info)
call gather_global_ext(G_dxt, I_dxt, master_task, distrb_info)
call gather_global_ext(G_dyt, I_dyt, master_task, distrb_info)
call gather_global_ext(G_stressp_1, I_stressp_1, master_task, distrb_info)
call gather_global_ext(G_stressp_2, I_stressp_2, master_task, distrb_info)
call gather_global_ext(G_stressp_3, I_stressp_3, master_task, distrb_info)
call gather_global_ext(G_stressp_4, I_stressp_4, master_task, distrb_info)
call gather_global_ext(G_stressm_1, I_stressm_1, master_task, distrb_info)
call gather_global_ext(G_stressm_2, I_stressm_2, master_task, distrb_info)
call gather_global_ext(G_stressm_3, I_stressm_3, master_task, distrb_info)
call gather_global_ext(G_stressm_4, I_stressm_4, master_task, distrb_info)
call gather_global_ext(G_stress12_1, I_stress12_1, master_task, distrb_info)
call gather_global_ext(G_stress12_2, I_stress12_2, master_task, distrb_info)
call gather_global_ext(G_stress12_3, I_stress12_3, master_task, distrb_info)
call gather_global_ext(G_stress12_4, I_stress12_4, master_task, distrb_info)
call gather_global_ext(G_icetmask, I_icetmask, master_task, distrb_info, 0 )
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It's curious that icetmask has 0 at the end of the argument list and iceumask has .false.! I see that is how they are defined. No need to fix it here...

call gather_global_ext(G_iceumask, I_iceumask, master_task, distrb_info, .false.)
!v1 call gather_global_ext(G_dxhy, I_dxhy, master_task, distrb_info )
!v1 call gather_global_ext(G_dyhx, I_dyhx, master_task, distrb_info )
!v1 call gather_global_ext(G_cyp, I_cyp, master_task, distrb_info )
!v1 call gather_global_ext(G_cxp, I_cxp, master_task, distrb_info )
!v1 call gather_global_ext(G_cym, I_cym, master_task, distrb_info )
!v1 call gather_global_ext(G_cxm, I_cxm, master_task, distrb_info )
!v1 call gather_global_ext(G_tinyarea, I_tinyarea, master_task, distrb_info )
!v1 call gather_global_ext(G_waterx, I_waterx, master_task, distrb_info )
!v1 call gather_global_ext(G_watery, I_watery, master_task, distrb_info )
call gather_global_ext(G_cdn_ocn, I_cdn_ocn, master_task, distrb_info )
call gather_global_ext(G_aiu, I_aiu, master_task, distrb_info )
call gather_global_ext(G_uocn, I_uocn, master_task, distrb_info )
call gather_global_ext(G_vocn, I_vocn, master_task, distrb_info )
call gather_global_ext(G_forcex, I_forcex, master_task, distrb_info )
call gather_global_ext(G_forcey, I_forcey, master_task, distrb_info )
call gather_global_ext(G_Tbu, I_Tbu, master_task, distrb_info )
call gather_global_ext(G_umassdti, I_umassdti, master_task, distrb_info )
call gather_global_ext(G_fm, I_fm, master_task, distrb_info )
call gather_global_ext(G_uarear, I_uarear, master_task, distrb_info )
call gather_global_ext(G_tarear, I_tarear, master_task, distrb_info )
call gather_global_ext(G_strintx, I_strintx, master_task, distrb_info )
call gather_global_ext(G_strinty, I_strinty, master_task, distrb_info )
call gather_global_ext(G_uvel_init, I_uvel_init, master_task, distrb_info )
call gather_global_ext(G_vvel_init, I_vvel_init, master_task, distrb_info )
call gather_global_ext(G_strength, I_strength, master_task, distrb_info )
call gather_global_ext(G_uvel, I_uvel, master_task, distrb_info, c0 )
call gather_global_ext(G_vvel, I_vvel, master_task, distrb_info, c0 )
call gather_global_ext(G_dxt, I_dxt, master_task, distrb_info )
call gather_global_ext(G_dyt, I_dyt, master_task, distrb_info )
call gather_global_ext(G_stressp_1, I_stressp_1, master_task, distrb_info )
call gather_global_ext(G_stressp_2, I_stressp_2, master_task, distrb_info )
call gather_global_ext(G_stressp_3, I_stressp_3, master_task, distrb_info )
call gather_global_ext(G_stressp_4, I_stressp_4, master_task, distrb_info )
call gather_global_ext(G_stressm_1, I_stressm_1, master_task, distrb_info )
call gather_global_ext(G_stressm_2, I_stressm_2, master_task, distrb_info )
call gather_global_ext(G_stressm_3, I_stressm_3, master_task, distrb_info )
call gather_global_ext(G_stressm_4, I_stressm_4, master_task, distrb_info )
call gather_global_ext(G_stress12_1, I_stress12_1, master_task, distrb_info )
call gather_global_ext(G_stress12_2, I_stress12_2, master_task, distrb_info )
call gather_global_ext(G_stress12_3, I_stress12_3, master_task, distrb_info )
call gather_global_ext(G_stress12_4, I_stress12_4, master_task, distrb_info )

!-- All calculations has to be done on the master-task --

Expand Down
6 changes: 5 additions & 1 deletion cicecore/cicedynB/general/ice_init.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -94,7 +94,8 @@ subroutine input_data
bathymetry_file, use_bathymetry, &
bathymetry_format, &
grid_type, grid_format, &
dxrect, dyrect
dxrect, dyrect, &
pgl_global_ext
use ice_dyn_shared, only: ndte, kdyn, revised_evp, yield_curve, &
kevp_kernel, &
seabed_stress, seabed_stress_method, &
Expand Down Expand Up @@ -314,6 +315,7 @@ subroutine input_data
ndtd = 1 ! dynamic time steps per thermodynamic time step
ndte = 120 ! subcycles per dynamics timestep: ndte=dt_dyn/dte
kevp_kernel = 0 ! EVP kernel (0 = 2D, >0: 1D. Only ver. 2 is implemented yet)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The comment here is kind of confusing. Does kevp_kernel=102 mean 1D, version 02? Do we have a version 1? Should we simplify the kevp_kernel choices? I remember that we chose 102 because it wasn't really ready, so this was a way to 'hide' it from users, or at least make it less likely they'd set it up.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Originally several versions of 1d evp was created. Version 1 included most of the refactoring. Version 2 moved derived grid parameters of HTE and HTN to the 1d solver as a scalar.
Version 3 changed some of the variables from real4 to real8. The influence of this was limited change especially when compared to the accuracy of the iteration. In the end only v2 was implemented. I think that kevp_kernel could be included in kdyn as option 4. I have not thought this thrue

Copy link
Contributor

@mhrib mhrib Mar 6, 2021 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Right.
102="Our version 2"
Version 1 and version 2 gave identical results (maybe except for really aggressive flags, I do not recall exactly). But v2 only takes about half of memory and is a bit faster. For conservative flags we were also able to produce BFB results.
That's not the case for v3, where many internal variables was calculated as real*4. But it was again faster and also takes up less memory. As I recall, the uncertainties calculated using uvel,vvel was less or comparable to the uncertainties obtained across different computational architecture.
v3 is an interesting exercise I think and worth to considering for special cases to reduce memory load (and gain of speed as well).

pgl_global_ext = .false. ! if true, init primary grid lebgths (global ext.)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

spell lengths

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
Expand Down Expand Up @@ -637,6 +639,7 @@ subroutine input_data
call broadcast_scalar(ndtd, master_task)
call broadcast_scalar(ndte, master_task)
call broadcast_scalar(kevp_kernel, 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)
Expand Down Expand Up @@ -1749,6 +1752,7 @@ subroutine input_data
if (kevp_kernel /= 0) then
if (kevp_kernel == 102) then
kevp_kernel = 2
if (my_task == master_task) pgl_global_ext = .true.
else
if (my_task == master_task) write(nu_diag,*) subname//' ERROR: kevp_kernel = ',kevp_kernel
if (kevp_kernel == 2) then
Expand Down
134 changes: 132 additions & 2 deletions cicecore/cicedynB/infrastructure/ice_grid.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -77,6 +77,10 @@ module ice_grid
ocn_gridcell_frac ! only relevant for lat-lon grids
! gridcell value of [1 - (land fraction)] (T-cell)

real (kind=dbl_kind), dimension (:,:), allocatable, public :: &
G_HTE , & ! length of eastern edge of T-cell (global ext.)
G_HTN ! length of northern edge of T-cell (global ext.)

real (kind=dbl_kind), dimension (:,:,:), allocatable, public :: &
cyp , & ! 1.5*HTE - 0.5*HTE
cxp , & ! 1.5*HTN - 0.5*HTN
Expand Down Expand Up @@ -125,7 +129,8 @@ module ice_grid
kmt ! ocean topography mask for bathymetry (T-cell)

logical (kind=log_kind), public :: &
use_bathymetry ! flag for reading in bathymetry_file
use_bathymetry, & ! flag for reading in bathymetry_file
pgl_global_ext ! flag for init primary grid lengths (global ext.)

logical (kind=log_kind), &
dimension (:,:,:), allocatable, public :: &
Expand Down Expand Up @@ -153,6 +158,8 @@ subroutine alloc_grid

integer (int_kind) :: ierr

character(len=*), parameter :: subname = '(alloc_grid)'

allocate( &
dxt (nx_block,ny_block,max_blocks), & ! width of T-cell through the middle (m)
dyt (nx_block,ny_block,max_blocks), & ! height of T-cell through the middle (m)
Expand Down Expand Up @@ -203,7 +210,15 @@ subroutine alloc_grid
mse (2,2,nx_block,ny_block,max_blocks), &
msw (2,2,nx_block,ny_block,max_blocks), &
stat=ierr)
if (ierr/=0) call abort_ice('(alloc_grid): Out of memory')
if (ierr/=0) call abort_ice(subname//'ERROR: Out of memory')

if (pgl_global_ext) then
allocate( &
G_HTE(nx_global+2*nghost, ny_global+2*nghost), & ! length of eastern edge of T-cell (global ext.)
G_HTN(nx_global+2*nghost, ny_global+2*nghost), & ! length of northern edge of T-cell (global ext.)
stat=ierr)
if (ierr/=0) call abort_ice(subname//'ERROR: Out of memory')
endif

end subroutine alloc_grid

Expand Down Expand Up @@ -1498,6 +1513,9 @@ subroutine primary_grid_lengths_HTN(work_g)
enddo
enddo
endif
if (pgl_global_ext) then
call primary_grid_lengths_global_ext(G_HTN, work_g)
endif
call scatter_global(HTN, work_g, master_task, distrb_info, &
field_loc_Nface, field_type_scalar)
call scatter_global(dxu, work_g2, master_task, distrb_info, &
Expand Down Expand Up @@ -1572,6 +1590,9 @@ subroutine primary_grid_lengths_HTE(work_g)
enddo
endif
endif
if (pgl_global_ext) then
call primary_grid_lengths_global_ext(G_HTE, work_g)
endif
call scatter_global(HTE, work_g, master_task, distrb_info, &
field_loc_Eface, field_type_scalar)
call scatter_global(dyu, work_g2, master_task, distrb_info, &
Expand Down Expand Up @@ -2555,6 +2576,115 @@ subroutine read_seabedstress_bathy

end subroutine read_seabedstress_bathy

!=======================================================================
! Initialize global primary grid lengths array with ghost cells from
! global primary grid lengths array

subroutine primary_grid_lengths_global_ext(ARRAY_O, ARRAY_I)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It makes more sense to me for this routine to be put in the boundary modules, since it fills the ghost cells. This is all done on master task so MPI doesn't matter, and yet it would need to be put in both mpi and serial directories... That's the down side of how the code is currently structured, but I think we should stick with it for this, otherwise things can get really confusing.


use ice_constants, only: c0

real (kind=dbl_kind), dimension(:,:), intent(in) :: &
ARRAY_I

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)'

if ((ns_boundary_type == 'tripole' ) .or. &
(ns_boundary_type == 'tripoleT')) then
call abort_ice(subname // 'ERROR: ' // &
ns_boundary_type // ' bndy type not impl for cfg')
endif

do jo = 1, (ny_global + 2 * nghost)
ji = -nghost + jo

! Southern ghost cells

if (ji < 1) then
select case (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 n-s bndy type')
end select
endif

! Northern ghost cells

if (ji > ny_global) then
select case (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 n-s bndy type')
end select
endif

do io = 1, (nx_global + 2 * nghost)
ii = -nghost + io

! Western ghost cells

if (ii < 1) then
select case (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 e-w bndy type')
end select
endif

! Eastern ghost cells

if (ii > nx_global) then
select case (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 e-w bndy 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 module ice_grid
Expand Down