Line data Source code
1 : #ifdef ncdf
2 : #define USE_NETCDF
3 : #endif
4 : !=======================================================================
5 :
6 : ! Spatial grids, masks, and boundary conditions
7 : !
8 : ! authors: Elizabeth C. Hunke and William H. Lipscomb, LANL
9 : ! Tony Craig, NCAR
10 : !
11 : ! 2004: Block structure added by William Lipscomb
12 : ! init_grid split into two parts as in POP 2.0
13 : ! Boundary update routines replaced by POP versions
14 : ! 2006: Converted to free source form (F90) by Elizabeth Hunke
15 : ! 2007: Option to read from netcdf files (A. Keen, Met Office)
16 : ! Grid reading routines reworked by E. Hunke for boundary values
17 : ! 2021: Add N (center of north face) and E (center of east face) grids
18 : ! to support C and CD solvers. Defining T at center of cells, U at
19 : ! NE corner, N at center of top face, E at center of right face.
20 : ! All cells are quadrilaterals with NE, E, and N associated with
21 : ! directions relative to logical grid.
22 :
23 : module ice_grid
24 :
25 : use ice_kinds_mod
26 : use ice_broadcast, only: broadcast_scalar, broadcast_array
27 : use ice_boundary, only: ice_HaloUpdate, ice_HaloExtrapolate, &
28 : primary_grid_lengths_global_ext
29 : use ice_communicate, only: my_task, master_task
30 : use ice_blocks, only: block, get_block, nx_block, ny_block, nghost
31 : use ice_domain_size, only: nx_global, ny_global, max_blocks
32 : use ice_domain, only: blocks_ice, nblocks, halo_info, distrb_info, &
33 : ew_boundary_type, ns_boundary_type, init_domain_distribution
34 : use ice_fileunits, only: nu_diag, nu_grid, nu_kmt, &
35 : get_fileunit, release_fileunit, flush_fileunit
36 : use ice_gather_scatter, only: gather_global, scatter_global
37 : use ice_read_write, only: ice_read, ice_read_nc, ice_read_global, &
38 : ice_read_global_nc, ice_open, ice_open_nc, ice_close_nc
39 : use ice_timers, only: timer_bound, ice_timer_start, ice_timer_stop
40 : use ice_exit, only: abort_ice
41 : use ice_global_reductions, only: global_minval, global_maxval
42 : use icepack_intfc, only: icepack_warnings_flush, icepack_warnings_aborted
43 : use icepack_intfc, only: icepack_query_parameters, icepack_init_parameters
44 :
45 : implicit none
46 : private
47 : public :: init_grid1, init_grid2, grid_average_X2Y, &
48 : alloc_grid, makemask, grid_neighbor_min, grid_neighbor_max
49 :
50 : character (len=char_len_long), public :: &
51 : grid_format , & ! file format ('bin'=binary or 'nc'=netcdf) ! LCOV_EXCL_LINE
52 : gridcpl_file , & ! input file for POP coupling grid info ! LCOV_EXCL_LINE
53 : grid_file , & ! input file for POP grid info ! LCOV_EXCL_LINE
54 : kmt_file , & ! input file for POP grid info ! LCOV_EXCL_LINE
55 : kmt_type , & ! options are file, default, boxislands ! LCOV_EXCL_LINE
56 : bathymetry_file, & ! input bathymetry for seabed stress ! LCOV_EXCL_LINE
57 : bathymetry_format, & ! bathymetry file format (default or pop) ! LCOV_EXCL_LINE
58 : grid_spacing , & ! default of 30.e3m or set by user in namelist ! LCOV_EXCL_LINE
59 : grid_ice , & ! Underlying model grid structure (A, B, C, CD) ! LCOV_EXCL_LINE
60 : grid_ice_thrm, & ! ocean forcing grid for thermo fields (T, U, N, E) ! LCOV_EXCL_LINE
61 : grid_ice_dynu, & ! ocean forcing grid for dyn U fields (T, U, N, E) ! LCOV_EXCL_LINE
62 : grid_ice_dynv, & ! ocean forcing grid for dyn V fields (T, U, N, E) ! LCOV_EXCL_LINE
63 : grid_atm , & ! atmos forcing grid structure (A, B, C, CD) ! LCOV_EXCL_LINE
64 : grid_atm_thrm, & ! atmos forcing grid for thermo fields (T, U, N, E) ! LCOV_EXCL_LINE
65 : grid_atm_dynu, & ! atmos forcing grid for dyn U fields (T, U, N, E) ! LCOV_EXCL_LINE
66 : grid_atm_dynv, & ! atmos forcing grid for dyn V fields (T, U, N, E) ! LCOV_EXCL_LINE
67 : grid_ocn , & ! ocean forcing grid structure (A B, C, CD) ! LCOV_EXCL_LINE
68 : grid_ocn_thrm, & ! ocean forcing grid for thermo fields (T, U, N, E) ! LCOV_EXCL_LINE
69 : grid_ocn_dynu, & ! ocean forcing grid for dyn U fields (T, U, N, E) ! LCOV_EXCL_LINE
70 : grid_ocn_dynv, & ! ocean forcing grid for dyn V fields (T, U, N, E) ! LCOV_EXCL_LINE
71 : grid_type ! current options are rectangular (default),
72 : ! displaced_pole, tripole, regional
73 :
74 : real (kind=dbl_kind), dimension (:,:,:), allocatable, public :: &
75 : dxT , & ! width of T-cell through the middle (m) ! LCOV_EXCL_LINE
76 : dyT , & ! height of T-cell through the middle (m) ! LCOV_EXCL_LINE
77 : dxU , & ! width of U-cell through the middle (m) ! LCOV_EXCL_LINE
78 : dyU , & ! height of U-cell through the middle (m) ! LCOV_EXCL_LINE
79 : dxN , & ! width of N-cell through the middle (m) ! LCOV_EXCL_LINE
80 : dyN , & ! height of N-cell through the middle (m) ! LCOV_EXCL_LINE
81 : dxE , & ! width of E-cell through the middle (m) ! LCOV_EXCL_LINE
82 : dyE , & ! height of E-cell through the middle (m) ! LCOV_EXCL_LINE
83 : HTE , & ! length of eastern edge of T-cell (m) ! LCOV_EXCL_LINE
84 : HTN , & ! length of northern edge of T-cell (m) ! LCOV_EXCL_LINE
85 : tarea , & ! area of T-cell (m^2), valid in halo ! LCOV_EXCL_LINE
86 : uarea , & ! area of U-cell (m^2), valid in halo ! LCOV_EXCL_LINE
87 : narea , & ! area of N-cell (m^2), valid in halo ! LCOV_EXCL_LINE
88 : earea , & ! area of E-cell (m^2), valid in halo ! LCOV_EXCL_LINE
89 : tarear , & ! 1/tarea, valid in halo ! LCOV_EXCL_LINE
90 : uarear , & ! 1/uarea, valid in halo ! LCOV_EXCL_LINE
91 : narear , & ! 1/narea, valid in halo ! LCOV_EXCL_LINE
92 : earear , & ! 1/earea, valid in halo ! LCOV_EXCL_LINE
93 : tarean , & ! area of NH T-cells ! LCOV_EXCL_LINE
94 : tareas , & ! area of SH T-cells ! LCOV_EXCL_LINE
95 : ULON , & ! longitude of velocity pts, NE corner of T pts (radians) ! LCOV_EXCL_LINE
96 : ULAT , & ! latitude of velocity pts, NE corner of T pts (radians) ! LCOV_EXCL_LINE
97 : TLON , & ! longitude of temp (T) pts (radians) ! LCOV_EXCL_LINE
98 : TLAT , & ! latitude of temp (T) pts (radians) ! LCOV_EXCL_LINE
99 : NLON , & ! longitude of center of north face of T pts (radians) ! LCOV_EXCL_LINE
100 : NLAT , & ! latitude of center of north face of T pts (radians) ! LCOV_EXCL_LINE
101 : ELON , & ! longitude of center of east face of T pts (radians) ! LCOV_EXCL_LINE
102 : ELAT , & ! latitude of center of east face of T pts (radians) ! LCOV_EXCL_LINE
103 : ANGLE , & ! for conversions between POP grid and lat/lon ! LCOV_EXCL_LINE
104 : ANGLET , & ! ANGLE converted to T-cells, valid in halo ! LCOV_EXCL_LINE
105 : bathymetry , & ! ocean depth, for grounding keels and bergs (m) ! LCOV_EXCL_LINE
106 : ocn_gridcell_frac ! only relevant for lat-lon grids
107 : ! gridcell value of [1 - (land fraction)] (T-cell)
108 :
109 : real (kind=dbl_kind), dimension (:,:), allocatable, public :: &
110 : G_HTE , & ! length of eastern edge of T-cell (global ext.) ! LCOV_EXCL_LINE
111 : G_HTN ! length of northern edge of T-cell (global ext.)
112 :
113 : real (kind=dbl_kind), dimension (:,:,:), allocatable, public :: &
114 : cyp , & ! 1.5*HTE(i,j)-0.5*HTW(i,j) = 1.5*HTE(i,j)-0.5*HTE(i-1,j) ! LCOV_EXCL_LINE
115 : cxp , & ! 1.5*HTN(i,j)-0.5*HTS(i,j) = 1.5*HTN(i,j)-0.5*HTN(i,j-1) ! LCOV_EXCL_LINE
116 : cym , & ! 0.5*HTE(i,j)-1.5*HTW(i,j) = 0.5*HTE(i,j)-1.5*HTE(i-1,j) ! LCOV_EXCL_LINE
117 : cxm , & ! 0.5*HTN(i,j)-1.5*HTS(i,j) = 0.5*HTN(i,j)-1.5*HTN(i,j-1) ! LCOV_EXCL_LINE
118 : dxhy , & ! 0.5*(HTE(i,j) - HTW(i,j)) = 0.5*(HTE(i,j) - HTE(i-1,j)) ! LCOV_EXCL_LINE
119 : dyhx ! 0.5*(HTN(i,j) - HTS(i,j)) = 0.5*(HTN(i,j) - HTN(i,j-1))
120 :
121 : real (kind=dbl_kind), dimension (:,:,:), allocatable, public :: &
122 : ratiodxN , & ! - dxN(i+1,j) / dxN(i,j) ! LCOV_EXCL_LINE
123 : ratiodyE , & ! - dyE(i ,j+1) / dyE(i,j) ! LCOV_EXCL_LINE
124 : ratiodxNr , & ! 1 / ratiodxN ! LCOV_EXCL_LINE
125 : ratiodyEr ! 1 / ratiodyE
126 :
127 : ! grid dimensions for rectangular grid
128 : real (kind=dbl_kind), public :: &
129 : dxrect, & ! user_specified spacing (cm) in x-direction (uniform HTN) ! LCOV_EXCL_LINE
130 : dyrect ! user_specified spacing (cm) in y-direction (uniform HTE)
131 :
132 : ! growth factor for variable spaced grid
133 : real (kind=dbl_kind), public :: &
134 : dxscale, & ! scale factor for grid spacing in x direction (e.g., 1.02) ! LCOV_EXCL_LINE
135 : dyscale ! scale factor for gird spacing in y direction (e.g., 1.02)
136 :
137 : real (kind=dbl_kind), public :: &
138 : lonrefrect, & ! lower left lon for rectgrid ! LCOV_EXCL_LINE
139 : latrefrect ! lower left lat for rectgrid
140 :
141 : ! Corners of grid boxes for history output
142 : real (kind=dbl_kind), dimension (:,:,:,:), allocatable, public :: &
143 : lont_bounds, & ! longitude of gridbox corners for T point ! LCOV_EXCL_LINE
144 : latt_bounds, & ! latitude of gridbox corners for T point ! LCOV_EXCL_LINE
145 : lonu_bounds, & ! longitude of gridbox corners for U point ! LCOV_EXCL_LINE
146 : latu_bounds, & ! latitude of gridbox corners for U point ! LCOV_EXCL_LINE
147 : lonn_bounds, & ! longitude of gridbox corners for N point ! LCOV_EXCL_LINE
148 : latn_bounds, & ! latitude of gridbox corners for N point ! LCOV_EXCL_LINE
149 : lone_bounds, & ! longitude of gridbox corners for E point ! LCOV_EXCL_LINE
150 : late_bounds ! latitude of gridbox corners for E point
151 :
152 : ! geometric quantities used for remapping transport
153 : real (kind=dbl_kind), dimension (:,:,:), allocatable, public :: &
154 : xav , & ! mean T-cell value of x ! LCOV_EXCL_LINE
155 : yav , & ! mean T-cell value of y ! LCOV_EXCL_LINE
156 : xxav , & ! mean T-cell value of xx ! LCOV_EXCL_LINE
157 : ! xyav , & ! mean T-cell value of xy ! LCOV_EXCL_LINE
158 : ! yyav , & ! mean T-cell value of yy ! LCOV_EXCL_LINE
159 : yyav ! mean T-cell value of yy
160 : ! xxxav, & ! mean T-cell value of xxx
161 : ! xxyav, & ! mean T-cell value of xxy ! LCOV_EXCL_LINE
162 : ! xyyav, & ! mean T-cell value of xyy ! LCOV_EXCL_LINE
163 : ! yyyav ! mean T-cell value of yyy
164 :
165 : real (kind=dbl_kind), &
166 : dimension (:,:,:,:,:), allocatable, public :: & ! LCOV_EXCL_LINE
167 : mne, & ! matrices used for coordinate transformations in remapping ! LCOV_EXCL_LINE
168 : mnw, & ! ne = northeast corner, nw = northwest, etc. ! LCOV_EXCL_LINE
169 : mse, & ! LCOV_EXCL_LINE
170 : msw
171 :
172 : ! masks
173 : real (kind=dbl_kind), dimension (:,:,:), allocatable, public :: &
174 : hm , & ! land/boundary mask, thickness (T-cell) ! LCOV_EXCL_LINE
175 : bm , & ! task/block id ! LCOV_EXCL_LINE
176 : uvm , & ! land/boundary mask (U-cell) ! LCOV_EXCL_LINE
177 : npm , & ! land/boundary mask (N-cell) ! LCOV_EXCL_LINE
178 : epm , & ! land/boundary mask (E-cell) ! LCOV_EXCL_LINE
179 : kmt ! ocean topography mask for bathymetry (T-cell)
180 :
181 : logical (kind=log_kind), public :: &
182 : use_bathymetry, & ! flag for reading in bathymetry_file ! LCOV_EXCL_LINE
183 : pgl_global_ext, & ! flag for init primary grid lengths (global ext.) ! LCOV_EXCL_LINE
184 : scale_dxdy ! flag to apply scale factor to vary dx/dy in rectgrid
185 :
186 : logical (kind=log_kind), dimension (:,:,:), allocatable, public :: &
187 : tmask , & ! land/boundary mask, thickness (T-cell) ! LCOV_EXCL_LINE
188 : umask , & ! land/boundary mask (U-cell) (1 if all surrounding T cells are ocean) ! LCOV_EXCL_LINE
189 : umaskCD, & ! land/boundary mask (U-cell) (1 if at least two surrounding T cells are ocean) ! LCOV_EXCL_LINE
190 : nmask , & ! land/boundary mask, (N-cell) ! LCOV_EXCL_LINE
191 : emask , & ! land/boundary mask, (E-cell) ! LCOV_EXCL_LINE
192 : lmask_n, & ! northern hemisphere mask ! LCOV_EXCL_LINE
193 : lmask_s ! southern hemisphere mask
194 :
195 : real (kind=dbl_kind), dimension (:,:,:), allocatable, public :: &
196 : rndex_global ! global index for local subdomain (dbl)
197 :
198 : logical (kind=log_kind), private :: &
199 : l_readCenter ! If anglet exist in grid file read it otherwise calculate it
200 :
201 : interface grid_average_X2Y
202 : module procedure grid_average_X2Y_base , &
203 : grid_average_X2Y_userwghts, & ! LCOV_EXCL_LINE
204 : grid_average_X2Y_NEversion
205 : end interface
206 :
207 : !=======================================================================
208 :
209 : contains
210 :
211 : !=======================================================================
212 : !
213 : ! Allocate space for all variables
214 : !
215 37 : subroutine alloc_grid
216 :
217 : integer (int_kind) :: ierr
218 :
219 : character(len=*), parameter :: subname = '(alloc_grid)'
220 :
221 : allocate( &
222 : dxT (nx_block,ny_block,max_blocks), & ! width of T-cell through the middle (m) ! LCOV_EXCL_LINE
223 : dyT (nx_block,ny_block,max_blocks), & ! height of T-cell through the middle (m) ! LCOV_EXCL_LINE
224 : dxU (nx_block,ny_block,max_blocks), & ! width of U-cell through the middle (m) ! LCOV_EXCL_LINE
225 : dyU (nx_block,ny_block,max_blocks), & ! height of U-cell through the middle (m) ! LCOV_EXCL_LINE
226 : dxN (nx_block,ny_block,max_blocks), & ! width of N-cell through the middle (m) ! LCOV_EXCL_LINE
227 : dyN (nx_block,ny_block,max_blocks), & ! height of N-cell through the middle (m) ! LCOV_EXCL_LINE
228 : dxE (nx_block,ny_block,max_blocks), & ! width of E-cell through the middle (m) ! LCOV_EXCL_LINE
229 : dyE (nx_block,ny_block,max_blocks), & ! height of E-cell through the middle (m) ! LCOV_EXCL_LINE
230 : HTE (nx_block,ny_block,max_blocks), & ! length of eastern edge of T-cell (m) ! LCOV_EXCL_LINE
231 : HTN (nx_block,ny_block,max_blocks), & ! length of northern edge of T-cell (m) ! LCOV_EXCL_LINE
232 : tarea (nx_block,ny_block,max_blocks), & ! area of T-cell (m^2) ! LCOV_EXCL_LINE
233 : uarea (nx_block,ny_block,max_blocks), & ! area of U-cell (m^2) ! LCOV_EXCL_LINE
234 : narea (nx_block,ny_block,max_blocks), & ! area of N-cell (m^2) ! LCOV_EXCL_LINE
235 : earea (nx_block,ny_block,max_blocks), & ! area of E-cell (m^2) ! LCOV_EXCL_LINE
236 : tarear (nx_block,ny_block,max_blocks), & ! 1/tarea ! LCOV_EXCL_LINE
237 : uarear (nx_block,ny_block,max_blocks), & ! 1/uarea ! LCOV_EXCL_LINE
238 : narear (nx_block,ny_block,max_blocks), & ! 1/narea ! LCOV_EXCL_LINE
239 : earear (nx_block,ny_block,max_blocks), & ! 1/earea ! LCOV_EXCL_LINE
240 : tarean (nx_block,ny_block,max_blocks), & ! area of NH T-cells ! LCOV_EXCL_LINE
241 : tareas (nx_block,ny_block,max_blocks), & ! area of SH T-cells ! LCOV_EXCL_LINE
242 : ULON (nx_block,ny_block,max_blocks), & ! longitude of U pts, NE corner (radians) ! LCOV_EXCL_LINE
243 : ULAT (nx_block,ny_block,max_blocks), & ! latitude of U pts, NE corner (radians) ! LCOV_EXCL_LINE
244 : TLON (nx_block,ny_block,max_blocks), & ! longitude of T pts (radians) ! LCOV_EXCL_LINE
245 : TLAT (nx_block,ny_block,max_blocks), & ! latitude of T pts (radians) ! LCOV_EXCL_LINE
246 : NLON (nx_block,ny_block,max_blocks), & ! longitude of N pts, N face (radians) ! LCOV_EXCL_LINE
247 : NLAT (nx_block,ny_block,max_blocks), & ! latitude of N pts, N face (radians) ! LCOV_EXCL_LINE
248 : ELON (nx_block,ny_block,max_blocks), & ! longitude of E pts, E face (radians) ! LCOV_EXCL_LINE
249 : ELAT (nx_block,ny_block,max_blocks), & ! latitude of E pts, E face (radians) ! LCOV_EXCL_LINE
250 : ANGLE (nx_block,ny_block,max_blocks), & ! for conversions between POP grid and lat/lon ! LCOV_EXCL_LINE
251 : ANGLET (nx_block,ny_block,max_blocks), & ! ANGLE converted to T-cells ! LCOV_EXCL_LINE
252 : bathymetry(nx_block,ny_block,max_blocks),& ! ocean depth, for grounding keels and bergs (m) ! LCOV_EXCL_LINE
253 : ocn_gridcell_frac(nx_block,ny_block,max_blocks),& ! only relevant for lat-lon grids ! LCOV_EXCL_LINE
254 : cyp (nx_block,ny_block,max_blocks), & ! 1.5*HTE - 0.5*HTW ! LCOV_EXCL_LINE
255 : cxp (nx_block,ny_block,max_blocks), & ! 1.5*HTN - 0.5*HTS ! LCOV_EXCL_LINE
256 : cym (nx_block,ny_block,max_blocks), & ! 0.5*HTE - 1.5*HTW ! LCOV_EXCL_LINE
257 : cxm (nx_block,ny_block,max_blocks), & ! 0.5*HTN - 1.5*HTS ! LCOV_EXCL_LINE
258 : dxhy (nx_block,ny_block,max_blocks), & ! 0.5*(HTE - HTW) ! LCOV_EXCL_LINE
259 : dyhx (nx_block,ny_block,max_blocks), & ! 0.5*(HTN - HTS) ! LCOV_EXCL_LINE
260 : xav (nx_block,ny_block,max_blocks), & ! mean T-cell value of x ! LCOV_EXCL_LINE
261 : yav (nx_block,ny_block,max_blocks), & ! mean T-cell value of y ! LCOV_EXCL_LINE
262 : xxav (nx_block,ny_block,max_blocks), & ! mean T-cell value of xx ! LCOV_EXCL_LINE
263 : yyav (nx_block,ny_block,max_blocks), & ! mean T-cell value of yy ! LCOV_EXCL_LINE
264 : hm (nx_block,ny_block,max_blocks), & ! land/boundary mask, thickness (T-cell) ! LCOV_EXCL_LINE
265 : bm (nx_block,ny_block,max_blocks), & ! task/block id ! LCOV_EXCL_LINE
266 : uvm (nx_block,ny_block,max_blocks), & ! land/boundary mask, velocity (U-cell) - water in case of all water point ! LCOV_EXCL_LINE
267 : npm (nx_block,ny_block,max_blocks), & ! land/boundary mask (N-cell) ! LCOV_EXCL_LINE
268 : epm (nx_block,ny_block,max_blocks), & ! land/boundary mask (E-cell) ! LCOV_EXCL_LINE
269 : kmt (nx_block,ny_block,max_blocks), & ! ocean topography mask for bathymetry (T-cell) ! LCOV_EXCL_LINE
270 : tmask (nx_block,ny_block,max_blocks), & ! land/boundary mask, thickness (T-cell) ! LCOV_EXCL_LINE
271 : umask (nx_block,ny_block,max_blocks), & ! land/boundary mask, velocity (U-cell) ! LCOV_EXCL_LINE
272 : umaskCD (nx_block,ny_block,max_blocks), & ! land/boundary mask, velocity (U-cell) ! LCOV_EXCL_LINE
273 : nmask (nx_block,ny_block,max_blocks), & ! land/boundary mask (N-cell) ! LCOV_EXCL_LINE
274 : emask (nx_block,ny_block,max_blocks), & ! land/boundary mask (E-cell) ! LCOV_EXCL_LINE
275 : lmask_n (nx_block,ny_block,max_blocks), & ! northern hemisphere mask ! LCOV_EXCL_LINE
276 : lmask_s (nx_block,ny_block,max_blocks), & ! southern hemisphere mask ! LCOV_EXCL_LINE
277 : rndex_global(nx_block,ny_block,max_blocks), & ! global index for local subdomain (dbl) ! LCOV_EXCL_LINE
278 : lont_bounds(4,nx_block,ny_block,max_blocks), & ! longitude of gridbox corners for T point ! LCOV_EXCL_LINE
279 : latt_bounds(4,nx_block,ny_block,max_blocks), & ! latitude of gridbox corners for T point ! LCOV_EXCL_LINE
280 : lonu_bounds(4,nx_block,ny_block,max_blocks), & ! longitude of gridbox corners for U point ! LCOV_EXCL_LINE
281 : latu_bounds(4,nx_block,ny_block,max_blocks), & ! latitude of gridbox corners for U point ! LCOV_EXCL_LINE
282 : lonn_bounds(4,nx_block,ny_block,max_blocks), & ! longitude of gridbox corners for N point ! LCOV_EXCL_LINE
283 : latn_bounds(4,nx_block,ny_block,max_blocks), & ! latitude of gridbox corners for N point ! LCOV_EXCL_LINE
284 : lone_bounds(4,nx_block,ny_block,max_blocks), & ! longitude of gridbox corners for E point ! LCOV_EXCL_LINE
285 : late_bounds(4,nx_block,ny_block,max_blocks), & ! latitude of gridbox corners for E point ! LCOV_EXCL_LINE
286 : mne (2,2,nx_block,ny_block,max_blocks), & ! matrices used for coordinate transformations in remapping ! LCOV_EXCL_LINE
287 : mnw (2,2,nx_block,ny_block,max_blocks), & ! ne = northeast corner, nw = northwest, etc. ! LCOV_EXCL_LINE
288 : mse (2,2,nx_block,ny_block,max_blocks), & ! LCOV_EXCL_LINE
289 : msw (2,2,nx_block,ny_block,max_blocks), & ! LCOV_EXCL_LINE
290 37 : stat=ierr)
291 37 : if (ierr/=0) call abort_ice(subname//'ERROR: Out of memory')
292 :
293 37 : if (grid_ice == 'CD' .or. grid_ice == 'C') then
294 : allocate( &
295 : ratiodxN (nx_block,ny_block,max_blocks), & ! LCOV_EXCL_LINE
296 : ratiodyE (nx_block,ny_block,max_blocks), & ! LCOV_EXCL_LINE
297 : ratiodxNr(nx_block,ny_block,max_blocks), & ! LCOV_EXCL_LINE
298 : ratiodyEr(nx_block,ny_block,max_blocks), & ! LCOV_EXCL_LINE
299 0 : stat=ierr)
300 0 : if (ierr/=0) call abort_ice(subname//'ERROR: Out of memory')
301 : endif
302 :
303 37 : if (pgl_global_ext) then
304 : allocate( &
305 : G_HTE(nx_global+2*nghost, ny_global+2*nghost), & ! length of eastern edge of T-cell (global ext.) ! LCOV_EXCL_LINE
306 : G_HTN(nx_global+2*nghost, ny_global+2*nghost), & ! length of northern edge of T-cell (global ext.) ! LCOV_EXCL_LINE
307 0 : stat=ierr)
308 0 : if (ierr/=0) call abort_ice(subname//'ERROR: Out of memory')
309 : endif
310 :
311 37 : end subroutine alloc_grid
312 :
313 : !=======================================================================
314 :
315 : ! Distribute blocks across processors. The distribution is optimized
316 : ! based on latitude and topography, contained in the ULAT and KMT arrays.
317 : !
318 : ! authors: William Lipscomb and Phil Jones, LANL
319 :
320 37 : subroutine init_grid1
321 :
322 : use ice_blocks, only: nx_block, ny_block
323 : use ice_broadcast, only: broadcast_array
324 : use ice_constants, only: c1
325 :
326 : integer (kind=int_kind) :: &
327 : fid_grid, & ! file id for netCDF grid file ! LCOV_EXCL_LINE
328 : fid_kmt ! file id for netCDF kmt file
329 :
330 : character (char_len) :: &
331 : fieldname ! field name in netCDF file
332 :
333 : real (kind=dbl_kind), dimension(:,:), allocatable :: &
334 37 : work_g1, work_g2
335 :
336 : real (kind=dbl_kind) :: &
337 8 : rad_to_deg
338 :
339 : character(len=*), parameter :: subname = '(init_grid1)'
340 :
341 : !-----------------------------------------------------------------
342 : ! Get global ULAT and KMT arrays used for block decomposition.
343 : !-----------------------------------------------------------------
344 :
345 37 : call icepack_query_parameters(rad_to_deg_out=rad_to_deg)
346 37 : call icepack_warnings_flush(nu_diag)
347 37 : if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
348 0 : file=__FILE__, line=__LINE__)
349 :
350 37 : allocate(work_g1(nx_global,ny_global))
351 37 : allocate(work_g2(nx_global,ny_global))
352 :
353 : ! check tripole flags here
354 : ! can't check in init_data because ns_boundary_type is not yet read
355 : ! can't check in init_domain_blocks because grid_type is not accessible due to circular logic
356 :
357 37 : if (grid_type == 'tripole' .and. ns_boundary_type /= 'tripole' .and. &
358 : ns_boundary_type /= 'tripoleT') then
359 : call abort_ice(subname//'ERROR: grid_type tripole needs tripole ns_boundary_type', &
360 0 : file=__FILE__, line=__LINE__)
361 : endif
362 :
363 37 : if (grid_type == 'tripole' .and. (mod(nx_global,2)/=0)) then
364 : call abort_ice(subname//'ERROR: grid_type tripole requires even nx_global number', &
365 0 : file=__FILE__, line=__LINE__)
366 : endif
367 :
368 : if (trim(grid_type) == 'displaced_pole' .or. &
369 : trim(grid_type) == 'tripole' .or. & ! LCOV_EXCL_LINE
370 : trim(grid_type) == 'regional' ) then
371 :
372 21 : if (trim(grid_format) == 'nc') then
373 :
374 0 : call ice_open_nc(grid_file,fid_grid)
375 0 : call ice_open_nc(kmt_file,fid_kmt)
376 :
377 0 : fieldname='ulat'
378 0 : call ice_read_global_nc(fid_grid,1,fieldname,work_g1,.true.)
379 0 : fieldname='kmt'
380 0 : call ice_read_global_nc(fid_kmt,1,fieldname,work_g2,.true.)
381 :
382 0 : if (my_task == master_task) then
383 0 : call ice_close_nc(fid_grid)
384 0 : call ice_close_nc(fid_kmt)
385 : endif
386 :
387 : else
388 :
389 21 : call ice_open(nu_grid,grid_file,64) ! ULAT
390 21 : call ice_open(nu_kmt, kmt_file, 32) ! KMT
391 :
392 21 : call ice_read_global(nu_grid,1,work_g1,'rda8',.true.) ! ULAT
393 21 : call ice_read_global(nu_kmt, 1,work_g2,'ida4',.true.) ! KMT
394 :
395 21 : if (my_task == master_task) then
396 5 : close (nu_grid)
397 5 : close (nu_kmt)
398 : endif
399 :
400 : endif
401 :
402 : else ! rectangular grid
403 :
404 264208 : work_g1(:,:) = 75._dbl_kind/rad_to_deg ! arbitrary polar latitude
405 264208 : work_g2(:,:) = c1
406 :
407 : endif
408 :
409 37 : call broadcast_array(work_g1, master_task) ! ULAT
410 37 : call broadcast_array(work_g2, master_task) ! KMT
411 :
412 : !-----------------------------------------------------------------
413 : ! distribute blocks among processors
414 : !-----------------------------------------------------------------
415 :
416 37 : call init_domain_distribution(work_g2, work_g1, grid_ice) ! KMT, ULAT
417 :
418 37 : deallocate(work_g1)
419 37 : deallocate(work_g2)
420 :
421 : !-----------------------------------------------------------------
422 : ! write additional domain information
423 : !-----------------------------------------------------------------
424 :
425 37 : if (my_task == master_task) then
426 7 : write(nu_diag,'(a26,i6)') ' Block size: nx_block = ',nx_block
427 7 : write(nu_diag,'(a26,i6)') ' ny_block = ',ny_block
428 : endif
429 :
430 74 : end subroutine init_grid1
431 :
432 : !=======================================================================
433 :
434 : ! Horizontal grid initialization:
435 : !
436 : ! U{LAT,LONG} = true {latitude,longitude} of U points
437 : ! HT{N,E} = cell widths on {N,E} sides of T cell
438 : ! ANGLE = angle between local x direction and true east
439 : ! hm = land mask (c1 for ocean points, c0 for land points)
440 : ! D{X,Y}{T,U} = {x,y} spacing centered at {T,U} points
441 : ! T-grid and ghost cell values
442 : ! Various grid quantities needed for dynamics and transport
443 : !
444 : ! author: Elizabeth C. Hunke, LANL
445 :
446 37 : subroutine init_grid2
447 :
448 : use ice_blocks, only: get_block, block, nx_block, ny_block
449 : use ice_constants, only: c0, c1, c2, p5, p25, c1p5, &
450 : field_loc_center, field_loc_NEcorner, field_loc_Nface, field_loc_Eface, & ! LCOV_EXCL_LINE
451 : field_type_scalar, field_type_vector, field_type_angle
452 : use ice_domain_size, only: max_blocks
453 : #if defined (_OPENMP)
454 : use OMP_LIB
455 : #endif
456 :
457 : integer (kind=int_kind) :: &
458 : i, j, iblk, & ! LCOV_EXCL_LINE
459 : ilo,ihi,jlo,jhi ! beginning and end of physical domain
460 :
461 : real (kind=dbl_kind) :: &
462 : angle_0, angle_w, angle_s, angle_sw, & ! LCOV_EXCL_LINE
463 8 : pi, pi2, puny
464 :
465 : logical (kind=log_kind), dimension(nx_block,ny_block,max_blocks):: &
466 74 : out_of_range
467 :
468 : real (kind=dbl_kind), dimension(:,:), allocatable :: &
469 37 : work_g1
470 :
471 : type (block) :: &
472 : this_block ! block information for current block
473 :
474 : #if defined (_OPENMP)
475 : integer(kind=omp_sched_kind) :: ompsk ! openmp schedule
476 : integer(kind=int_kind) :: ompcs ! openmp schedule count
477 : #endif
478 :
479 : character(len=*), parameter :: subname = '(init_grid2)'
480 :
481 : !-----------------------------------------------------------------
482 : ! lat, lon, cell widths, angle, land mask
483 : !-----------------------------------------------------------------
484 :
485 37 : l_readCenter = .false.
486 :
487 37 : call icepack_query_parameters(pi_out=pi, pi2_out=pi2, puny_out=puny)
488 37 : call icepack_warnings_flush(nu_diag)
489 37 : if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
490 0 : file=__FILE__, line=__LINE__)
491 :
492 : if (trim(grid_type) == 'displaced_pole' .or. &
493 : trim(grid_type) == 'tripole' .or. & ! LCOV_EXCL_LINE
494 : trim(grid_type) == 'regional' ) then
495 21 : if (trim(grid_format) == 'nc') then
496 0 : call popgrid_nc ! read POP grid lengths from nc file
497 : else
498 21 : call popgrid ! read POP grid lengths directly
499 : endif
500 : #ifdef CESMCOUPLED
501 : elseif (trim(grid_type) == 'latlon') then
502 : call latlongrid ! lat lon grid for sequential CESM (CAM mode)
503 : return
504 : #endif
505 16 : elseif (trim(grid_type) == 'cpom_grid') then
506 0 : call cpomgrid ! cpom model orca1 type grid
507 : else
508 16 : call rectgrid ! regular rectangular grid
509 : endif
510 :
511 : !-----------------------------------------------------------------
512 : ! Diagnose OpenMP thread schedule, force order in output
513 : !-----------------------------------------------------------------
514 :
515 : #if defined (_OPENMP)
516 20 : !$OMP PARALLEL DO ORDERED PRIVATE(iblk) SCHEDULE(runtime)
517 : do iblk = 1, nblocks
518 : if (my_task == master_task) then
519 : !$OMP ORDERED
520 : if (iblk == 1) then
521 : call omp_get_schedule(ompsk,ompcs)
522 : ! write(nu_diag,*) ''
523 : write(nu_diag,*) subname,' OpenMP runtime thread schedule:'
524 : write(nu_diag,*) subname,' omp schedule = ',ompsk,ompcs
525 : endif
526 : write(nu_diag,*) subname,' block, thread = ',iblk,OMP_GET_THREAD_NUM()
527 : !$OMP END ORDERED
528 : endif
529 : enddo
530 : !$OMP END PARALLEL DO
531 20 : call flush_fileunit(nu_diag)
532 : #endif
533 :
534 : !-----------------------------------------------------------------
535 : ! T-grid cell and U-grid cell quantities
536 : ! Fill halo data locally where possible to avoid missing
537 : ! data associated with land block elimination
538 : ! Note: HTN, HTE, dx*, dy* are all defined from global arrays
539 : ! at halos.
540 : !-----------------------------------------------------------------
541 :
542 20 : !$OMP PARALLEL DO PRIVATE(iblk,i,j,ilo,ihi,jlo,jhi,this_block)
543 50 : do iblk = 1, nblocks
544 33 : this_block = get_block(blocks_ice(iblk),iblk)
545 33 : ilo = this_block%ilo
546 33 : ihi = this_block%ihi
547 33 : jlo = this_block%jlo
548 33 : jhi = this_block%jhi
549 :
550 1239 : do j = 1,ny_block
551 50267 : do i = 1,nx_block
552 49028 : tarea(i,j,iblk) = dxT(i,j,iblk)*dyT(i,j,iblk)
553 49028 : uarea(i,j,iblk) = dxU(i,j,iblk)*dyU(i,j,iblk)
554 49028 : narea(i,j,iblk) = dxN(i,j,iblk)*dyN(i,j,iblk)
555 49028 : earea(i,j,iblk) = dxE(i,j,iblk)*dyE(i,j,iblk)
556 :
557 49028 : if (tarea(i,j,iblk) > c0) then
558 49028 : tarear(i,j,iblk) = c1/tarea(i,j,iblk)
559 : else
560 0 : tarear(i,j,iblk) = c0 ! possible on boundaries
561 : endif
562 49028 : if (uarea(i,j,iblk) > c0) then
563 49028 : uarear(i,j,iblk) = c1/uarea(i,j,iblk)
564 : else
565 0 : uarear(i,j,iblk) = c0 ! possible on boundaries
566 : endif
567 49028 : if (narea(i,j,iblk) > c0) then
568 49028 : narear(i,j,iblk) = c1/narea(i,j,iblk)
569 : else
570 0 : narear(i,j,iblk) = c0 ! possible on boundaries
571 : endif
572 50234 : if (earea(i,j,iblk) > c0) then
573 49028 : earear(i,j,iblk) = c1/earea(i,j,iblk)
574 : else
575 0 : earear(i,j,iblk) = c0 ! possible on boundaries
576 : endif
577 :
578 : enddo
579 : enddo
580 :
581 1173 : do j = jlo, jhi
582 45541 : do i = ilo, ihi
583 44368 : dxhy(i,j,iblk) = p5*(HTE(i,j,iblk) - HTE(i-1,j,iblk))
584 45508 : dyhx(i,j,iblk) = p5*(HTN(i,j,iblk) - HTN(i,j-1,iblk))
585 : enddo
586 : enddo
587 :
588 1206 : do j = jlo, jhi+1
589 47871 : do i = ilo, ihi+1
590 46665 : cyp(i,j,iblk) = (c1p5*HTE(i,j,iblk) - p5*HTE(i-1,j,iblk))
591 46665 : cxp(i,j,iblk) = (c1p5*HTN(i,j,iblk) - p5*HTN(i,j-1,iblk))
592 : ! match order of operations in cyp, cxp for tripole grids
593 46665 : cym(i,j,iblk) = -(c1p5*HTE(i-1,j,iblk) - p5*HTE(i,j,iblk))
594 47838 : cxm(i,j,iblk) = -(c1p5*HTN(i,j-1,iblk) - p5*HTN(i,j,iblk))
595 : enddo
596 : enddo
597 :
598 70 : if (grid_ice == 'CD' .or. grid_ice == 'C') then
599 0 : do j = jlo, jhi
600 0 : do i = ilo, ihi
601 0 : ratiodxN (i,j,iblk) = - dxN(i+1,j ,iblk) / dxN(i,j,iblk)
602 0 : ratiodyE (i,j,iblk) = - dyE(i ,j+1,iblk) / dyE(i,j,iblk)
603 0 : ratiodxNr(i,j,iblk) = c1 / ratiodxN(i,j,iblk)
604 0 : ratiodyEr(i,j,iblk) = c1 / ratiodyE(i,j,iblk)
605 : enddo
606 : enddo
607 : endif
608 :
609 : enddo ! iblk
610 : !$OMP END PARALLEL DO
611 :
612 : !-----------------------------------------------------------------
613 : ! Ghost cell updates
614 : ! On the tripole grid, one must be careful with updates of
615 : ! quantities that involve a difference of cell lengths.
616 : ! For example, dyhx and dxhy are cell-centered vector components.
617 : ! Also note that on the tripole grid, cxp and cxm would swap places,
618 : ! as would cyp and cym. These quantities are computed only
619 : ! in north and east ghost cells (above), not south and west.
620 : !-----------------------------------------------------------------
621 :
622 37 : call ice_timer_start(timer_bound)
623 :
624 : call ice_HaloUpdate (dxhy, halo_info, &
625 : field_loc_center, field_type_vector, & ! LCOV_EXCL_LINE
626 37 : fillValue=c1)
627 : call ice_HaloUpdate (dyhx, halo_info, &
628 : field_loc_center, field_type_vector, & ! LCOV_EXCL_LINE
629 37 : fillValue=c1)
630 :
631 : ! Update just on the tripole seam to ensure bit-for-bit symmetry across seam
632 : call ice_HaloUpdate (tarea, halo_info, &
633 : field_loc_center, field_type_scalar, & ! LCOV_EXCL_LINE
634 37 : fillValue=c1, tripoleOnly=.true.)
635 : call ice_HaloUpdate (uarea, halo_info, &
636 : field_loc_NEcorner, field_type_scalar, & ! LCOV_EXCL_LINE
637 37 : fillValue=c1, tripoleOnly=.true.)
638 : call ice_HaloUpdate (narea, halo_info, &
639 : field_loc_Nface, field_type_scalar, & ! LCOV_EXCL_LINE
640 37 : fillValue=c1, tripoleOnly=.true.)
641 : call ice_HaloUpdate (earea, halo_info, &
642 : field_loc_Eface, field_type_scalar, & ! LCOV_EXCL_LINE
643 37 : fillValue=c1, tripoleOnly=.true.)
644 : call ice_HaloUpdate (tarear, halo_info, &
645 : field_loc_center, field_type_scalar, & ! LCOV_EXCL_LINE
646 37 : fillValue=c1, tripoleOnly=.true.)
647 : call ice_HaloUpdate (uarear, halo_info, &
648 : field_loc_NEcorner, field_type_scalar, & ! LCOV_EXCL_LINE
649 37 : fillValue=c1, tripoleOnly=.true.)
650 : call ice_HaloUpdate (narear, halo_info, &
651 : field_loc_Nface, field_type_scalar, & ! LCOV_EXCL_LINE
652 37 : fillValue=c1, tripoleOnly=.true.)
653 : call ice_HaloUpdate (earear, halo_info, &
654 : field_loc_Eface, field_type_scalar, & ! LCOV_EXCL_LINE
655 37 : fillValue=c1, tripoleOnly=.true.)
656 :
657 37 : call ice_timer_stop(timer_bound)
658 :
659 : !-----------------------------------------------------------------
660 : ! Calculate ANGLET to be compatible with POP ocean model
661 : ! First, ensure that -pi <= ANGLE <= pi
662 : !-----------------------------------------------------------------
663 :
664 111936 : out_of_range = .false.
665 111936 : where (ANGLE < -pi .or. ANGLE > pi) out_of_range = .true.
666 111936 : if (count(out_of_range) > 0) then
667 0 : write(nu_diag,*) subname,' angle = ',minval(ANGLE),maxval(ANGLE),count(out_of_range)
668 : call abort_ice (subname//' ANGLE out of expected range', &
669 0 : file=__FILE__, line=__LINE__)
670 : endif
671 :
672 : !-----------------------------------------------------------------
673 : ! Compute ANGLE on T-grid
674 : !-----------------------------------------------------------------
675 37 : if (trim(grid_type) == 'cpom_grid') then
676 0 : ANGLET(:,:,:) = ANGLE(:,:,:)
677 37 : else if (.not. (l_readCenter)) then
678 111936 : ANGLET = c0
679 :
680 : !$OMP PARALLEL DO PRIVATE(iblk,i,j,ilo,ihi,jlo,jhi,this_block, &
681 20 : !$OMP angle_0,angle_w,angle_s,angle_sw)
682 50 : do iblk = 1, nblocks
683 33 : this_block = get_block(blocks_ice(iblk),iblk)
684 33 : ilo = this_block%ilo
685 33 : ihi = this_block%ihi
686 33 : jlo = this_block%jlo
687 33 : jhi = this_block%jhi
688 :
689 1210 : do j = jlo, jhi
690 45541 : do i = ilo, ihi
691 44368 : angle_0 = ANGLE(i ,j ,iblk) ! w----0
692 44368 : angle_w = ANGLE(i-1,j ,iblk) ! | |
693 44368 : angle_s = ANGLE(i, j-1,iblk) ! | |
694 44368 : angle_sw = ANGLE(i-1,j-1,iblk) ! sw---s
695 : ANGLET(i,j,iblk) = atan2(p25*(sin(angle_0)+ &
696 : sin(angle_w)+ & ! LCOV_EXCL_LINE
697 : sin(angle_s)+ & ! LCOV_EXCL_LINE
698 : sin(angle_sw)),& ! LCOV_EXCL_LINE
699 : p25*(cos(angle_0)+ & ! LCOV_EXCL_LINE
700 : cos(angle_w)+ & ! LCOV_EXCL_LINE
701 : cos(angle_s)+ & ! LCOV_EXCL_LINE
702 45508 : cos(angle_sw)))
703 : enddo
704 : enddo
705 : enddo
706 : !$OMP END PARALLEL DO
707 : endif ! cpom_grid
708 37 : if (trim(grid_type) == 'regional' .and. &
709 : (.not. (l_readCenter))) then
710 : ! for W boundary extrapolate from interior
711 0 : !$OMP PARALLEL DO PRIVATE(iblk,i,j,ilo,ihi,jlo,jhi,this_block)
712 0 : do iblk = 1, nblocks
713 0 : this_block = get_block(blocks_ice(iblk),iblk)
714 0 : ilo = this_block%ilo
715 0 : ihi = this_block%ihi
716 0 : jlo = this_block%jlo
717 0 : jhi = this_block%jhi
718 :
719 0 : i = ilo
720 0 : if (this_block%i_glob(i) == 1) then
721 0 : do j = jlo, jhi
722 0 : ANGLET(i,j,iblk) = c2*ANGLET(i+1,j,iblk)-ANGLET(i+2,j,iblk)
723 : enddo
724 : endif
725 : enddo
726 : !$OMP END PARALLEL DO
727 : endif ! regional
728 :
729 37 : call ice_timer_start(timer_bound)
730 : call ice_HaloUpdate (ANGLET, halo_info, &
731 : field_loc_center, field_type_angle, & ! LCOV_EXCL_LINE
732 37 : fillValue=c1)
733 37 : call ice_timer_stop(timer_bound)
734 :
735 37 : call makemask ! velocity mask, hemisphere masks
736 37 : if (.not. (l_readCenter)) then
737 37 : call Tlatlon ! get lat, lon on the T grid
738 : endif
739 : !-----------------------------------------------------------------
740 : ! bathymetry
741 : !-----------------------------------------------------------------
742 :
743 37 : if (trim(bathymetry_format) == 'default') then
744 37 : call get_bathymetry
745 0 : elseif (trim(bathymetry_format) == 'pop') then
746 0 : call get_bathymetry_popfile
747 : else
748 : call abort_ice(subname//'ERROR: bathymetry_format value must be default or pop', &
749 0 : file=__FILE__, line=__LINE__)
750 : endif
751 :
752 : !----------------------------------------------------------------
753 : ! Corner coordinates for CF compliant history files
754 : !----------------------------------------------------------------
755 :
756 37 : call gridbox_corners
757 37 : call gridbox_edges
758 :
759 : !-----------------------------------------------------------------
760 : ! Compute global index (used for unpacking messages from coupler)
761 : !-----------------------------------------------------------------
762 :
763 37 : if (my_task==master_task) then
764 7 : allocate(work_g1(nx_global,ny_global))
765 843 : do j=1,ny_global
766 91611 : do i=1,nx_global
767 91604 : work_g1(i,j) = real((j-1)*nx_global + i,kind=dbl_kind)
768 : enddo
769 : enddo
770 : else
771 30 : allocate(work_g1(1,1)) ! to save memory
772 : endif
773 :
774 : call scatter_global(rndex_global, work_g1, &
775 : master_task, distrb_info, & ! LCOV_EXCL_LINE
776 37 : field_loc_center, field_type_scalar)
777 :
778 37 : deallocate(work_g1)
779 :
780 37 : end subroutine init_grid2
781 :
782 : !=======================================================================
783 :
784 : ! POP displaced pole grid and land mask (or tripole).
785 : ! Grid record number, field and units are: \\
786 : ! (1) ULAT (radians) \\
787 : ! (2) ULON (radians) \\
788 : ! (3) HTN (cm) \\
789 : ! (4) HTE (cm) \\
790 : ! (5) HUS (cm) \\
791 : ! (6) HUW (cm) \\
792 : ! (7) ANGLE (radians)
793 : !
794 : ! Land mask record number and field is (1) KMT.
795 : !
796 : ! author: Elizabeth C. Hunke, LANL
797 :
798 21 : subroutine popgrid
799 :
800 : use ice_blocks, only: nx_block, ny_block
801 : use ice_constants, only: c0, c1, p5, &
802 : field_loc_center, field_loc_NEcorner, & ! LCOV_EXCL_LINE
803 : field_type_scalar, field_type_angle
804 : use ice_domain_size, only: max_blocks
805 :
806 : integer (kind=int_kind) :: &
807 : i, j, iblk, & ! LCOV_EXCL_LINE
808 : ilo,ihi,jlo,jhi ! beginning and end of physical domain
809 :
810 : logical (kind=log_kind) :: diag
811 :
812 : real (kind=dbl_kind), dimension(:,:), allocatable :: &
813 21 : work_g1
814 :
815 : real (kind=dbl_kind), dimension (nx_block,ny_block,max_blocks) :: &
816 27842 : work1
817 :
818 : type (block) :: &
819 : this_block ! block information for current block
820 :
821 : character(len=*), parameter :: subname = '(popgrid)'
822 :
823 21 : call ice_open(nu_grid,grid_file,64)
824 21 : call ice_open(nu_kmt,kmt_file,32)
825 :
826 21 : diag = .true. ! write diagnostic info
827 :
828 : !-----------------------------------------------------------------
829 : ! topography
830 : !-----------------------------------------------------------------
831 :
832 : call ice_read(nu_kmt,1,work1,'ida4',diag, &
833 : field_loc=field_loc_center, & ! LCOV_EXCL_LINE
834 21 : field_type=field_type_scalar)
835 :
836 73808 : hm (:,:,:) = c0
837 73808 : kmt(:,:,:) = c0
838 20 : !$OMP PARALLEL DO PRIVATE(iblk,i,j,ilo,ihi,jlo,jhi,this_block)
839 2 : do iblk = 1, nblocks
840 1 : this_block = get_block(blocks_ice(iblk),iblk)
841 1 : ilo = this_block%ilo
842 1 : ihi = this_block%ihi
843 1 : jlo = this_block%jlo
844 1 : jhi = this_block%jhi
845 :
846 138 : do j = jlo, jhi
847 11717 : do i = ilo, ihi
848 11600 : kmt(i,j,iblk) = work1(i,j,iblk)
849 11716 : if (kmt(i,j,iblk) >= p5) hm(i,j,iblk) = c1
850 : enddo
851 : enddo
852 : enddo
853 : !$OMP END PARALLEL DO
854 :
855 : !-----------------------------------------------------------------
856 : ! lat, lon, angle
857 : !-----------------------------------------------------------------
858 :
859 21 : allocate(work_g1(nx_global,ny_global))
860 :
861 21 : call ice_read_global(nu_grid,1,work_g1,'rda8',.true.) ! ULAT
862 21 : call gridbox_verts(work_g1,latt_bounds)
863 : call scatter_global(ULAT, work_g1, master_task, distrb_info, &
864 21 : field_loc_NEcorner, field_type_scalar)
865 : call ice_HaloExtrapolate(ULAT, distrb_info, &
866 21 : ew_boundary_type, ns_boundary_type)
867 :
868 21 : call ice_read_global(nu_grid,2,work_g1,'rda8',.true.) ! ULON
869 21 : call gridbox_verts(work_g1,lont_bounds)
870 : call scatter_global(ULON, work_g1, master_task, distrb_info, &
871 21 : field_loc_NEcorner, field_type_scalar)
872 : call ice_HaloExtrapolate(ULON, distrb_info, &
873 21 : ew_boundary_type, ns_boundary_type)
874 :
875 21 : call ice_read_global(nu_grid,7,work_g1,'rda8',.true.) ! ANGLE
876 : call scatter_global(ANGLE, work_g1, master_task, distrb_info, &
877 21 : field_loc_NEcorner, field_type_angle)
878 :
879 : !-----------------------------------------------------------------
880 : ! cell dimensions
881 : ! calculate derived quantities from global arrays to preserve
882 : ! information on boundaries
883 : !-----------------------------------------------------------------
884 :
885 21 : call ice_read_global(nu_grid,3,work_g1,'rda8',.true.) ! HTN
886 21 : call primary_grid_lengths_HTN(work_g1) ! dxU, dxT, dxN, dxE
887 :
888 21 : call ice_read_global(nu_grid,4,work_g1,'rda8',.true.) ! HTE
889 21 : call primary_grid_lengths_HTE(work_g1) ! dyU, dyT, dyN, dyE
890 :
891 21 : deallocate(work_g1)
892 :
893 21 : if (my_task == master_task) then
894 5 : close (nu_grid)
895 5 : close (nu_kmt)
896 : endif
897 :
898 21 : end subroutine popgrid
899 :
900 : !=======================================================================
901 :
902 : ! POP displaced pole grid and land mask.
903 : ! Grid record number, field and units are: \\
904 : ! (1) ULAT (radians) \\
905 : ! (2) ULON (radians) \\
906 : ! (3) HTN (cm) \\
907 : ! (4) HTE (cm) \\
908 : ! (5) HUS (cm) \\
909 : ! (6) HUW (cm) \\
910 : ! (7) ANGLE (radians)
911 : !
912 : ! Land mask record number and field is (1) KMT.
913 : !
914 : ! author: Elizabeth C. Hunke, LANL
915 : ! Revised for netcdf input: Ann Keen, Met Office, May 2007
916 :
917 0 : subroutine popgrid_nc
918 :
919 : use ice_blocks, only: nx_block, ny_block
920 : use ice_constants, only: c0, c1, &
921 : field_loc_center, field_loc_NEcorner, & ! LCOV_EXCL_LINE
922 : field_type_scalar, field_type_angle
923 : use ice_domain_size, only: max_blocks
924 : #ifdef USE_NETCDF
925 : use netcdf
926 : #endif
927 :
928 : integer (kind=int_kind) :: &
929 : i, j, iblk, & ! LCOV_EXCL_LINE
930 : ilo,ihi,jlo,jhi, & ! beginning and end of physical domain ! LCOV_EXCL_LINE
931 : fid_grid, & ! file id for netCDF grid file ! LCOV_EXCL_LINE
932 : fid_kmt ! file id for netCDF kmt file
933 :
934 : logical (kind=log_kind) :: diag
935 :
936 : character (char_len) :: &
937 : fieldname ! field name in netCDF file
938 :
939 : real (kind=dbl_kind) :: &
940 0 : pi
941 :
942 : real (kind=dbl_kind), dimension(:,:), allocatable :: &
943 0 : work_g1
944 :
945 : real (kind=dbl_kind), dimension (nx_block,ny_block,max_blocks) :: &
946 0 : work1
947 :
948 : type (block) :: &
949 : this_block ! block information for current block
950 :
951 : integer(kind=int_kind) :: &
952 : varid
953 : integer (kind=int_kind) :: &
954 : status ! status flag
955 :
956 :
957 : character(len=*), parameter :: subname = '(popgrid_nc)'
958 :
959 : #ifdef USE_NETCDF
960 0 : call icepack_query_parameters(pi_out=pi)
961 0 : call icepack_warnings_flush(nu_diag)
962 0 : if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
963 0 : file=__FILE__, line=__LINE__)
964 :
965 0 : call ice_open_nc(grid_file,fid_grid)
966 0 : call ice_open_nc(kmt_file,fid_kmt)
967 :
968 0 : diag = .true. ! write diagnostic info
969 : !-----------------------------------------------------------------
970 : ! topography
971 : !-----------------------------------------------------------------
972 :
973 0 : fieldname='kmt'
974 : call ice_read_nc(fid_kmt,1,fieldname,work1,diag, &
975 : field_loc=field_loc_center, & ! LCOV_EXCL_LINE
976 0 : field_type=field_type_scalar)
977 :
978 0 : hm (:,:,:) = c0
979 0 : kmt(:,:,:) = c0
980 0 : !$OMP PARALLEL DO PRIVATE(iblk,i,j,ilo,ihi,jlo,jhi,this_block)
981 0 : do iblk = 1, nblocks
982 0 : this_block = get_block(blocks_ice(iblk),iblk)
983 0 : ilo = this_block%ilo
984 0 : ihi = this_block%ihi
985 0 : jlo = this_block%jlo
986 0 : jhi = this_block%jhi
987 :
988 0 : do j = jlo, jhi
989 0 : do i = ilo, ihi
990 0 : kmt(i,j,iblk) = work1(i,j,iblk)
991 0 : if (kmt(i,j,iblk) >= c1) hm(i,j,iblk) = c1
992 : enddo
993 : enddo
994 : enddo
995 : !$OMP END PARALLEL DO
996 :
997 : !-----------------------------------------------------------------
998 : ! lat, lon, angle
999 : !-----------------------------------------------------------------
1000 :
1001 0 : allocate(work_g1(nx_global,ny_global))
1002 :
1003 0 : fieldname='ulat'
1004 0 : call ice_read_global_nc(fid_grid,1,fieldname,work_g1,diag) ! ULAT
1005 0 : call gridbox_verts(work_g1,latt_bounds)
1006 : call scatter_global(ULAT, work_g1, master_task, distrb_info, &
1007 0 : field_loc_NEcorner, field_type_scalar)
1008 : call ice_HaloExtrapolate(ULAT, distrb_info, &
1009 0 : ew_boundary_type, ns_boundary_type)
1010 :
1011 0 : fieldname='ulon'
1012 0 : call ice_read_global_nc(fid_grid,1,fieldname,work_g1,diag) ! ULON
1013 0 : call gridbox_verts(work_g1,lont_bounds)
1014 : call scatter_global(ULON, work_g1, master_task, distrb_info, &
1015 0 : field_loc_NEcorner, field_type_scalar)
1016 : call ice_HaloExtrapolate(ULON, distrb_info, &
1017 0 : ew_boundary_type, ns_boundary_type)
1018 :
1019 0 : fieldname='angle'
1020 0 : call ice_read_global_nc(fid_grid,1,fieldname,work_g1,diag) ! ANGLE
1021 : call scatter_global(ANGLE, work_g1, master_task, distrb_info, &
1022 0 : field_loc_NEcorner, field_type_angle)
1023 : ! fix ANGLE: roundoff error due to single precision
1024 0 : where (ANGLE > pi) ANGLE = pi
1025 0 : where (ANGLE < -pi) ANGLE = -pi
1026 :
1027 : ! if grid file includes anglet then read instead
1028 0 : fieldname='anglet'
1029 0 : if (my_task == master_task) then
1030 0 : status = nf90_inq_varid(fid_grid, trim(fieldname) , varid)
1031 0 : if (status /= nf90_noerr) then
1032 0 : write(nu_diag,*) subname//' CICE will calculate angleT, TLON and TLAT'
1033 : else
1034 0 : write(nu_diag,*) subname//' angleT, TLON and TLAT is read from grid file'
1035 0 : l_readCenter = .true.
1036 : endif
1037 : endif
1038 0 : call broadcast_scalar(l_readCenter,master_task)
1039 0 : if (l_readCenter) then
1040 0 : call ice_read_global_nc(fid_grid,1,fieldname,work_g1,diag)
1041 : call scatter_global(ANGLET, work_g1, master_task, distrb_info, &
1042 0 : field_loc_center, field_type_angle)
1043 0 : where (ANGLET > pi) ANGLET = pi
1044 0 : where (ANGLET < -pi) ANGLET = -pi
1045 0 : fieldname="tlon"
1046 0 : call ice_read_global_nc(fid_grid,1,fieldname,work_g1,diag)
1047 : call scatter_global(TLON, work_g1, master_task, distrb_info, &
1048 0 : field_loc_center, field_type_scalar)
1049 0 : fieldname="tlat"
1050 0 : call ice_read_global_nc(fid_grid,1,fieldname,work_g1,diag)
1051 : call scatter_global(TLAT, work_g1, master_task, distrb_info, &
1052 0 : field_loc_center, field_type_scalar)
1053 : endif
1054 : !-----------------------------------------------------------------
1055 : ! cell dimensions
1056 : ! calculate derived quantities from global arrays to preserve
1057 : ! information on boundaries
1058 : !-----------------------------------------------------------------
1059 :
1060 0 : fieldname='htn'
1061 0 : call ice_read_global_nc(fid_grid,1,fieldname,work_g1,diag) ! HTN
1062 0 : call primary_grid_lengths_HTN(work_g1) ! dxU, dxT, dxN, dxE
1063 0 : fieldname='hte'
1064 0 : call ice_read_global_nc(fid_grid,1,fieldname,work_g1,diag) ! HTE
1065 0 : call primary_grid_lengths_HTE(work_g1) ! dyU, dyT, dyN, dyE
1066 :
1067 0 : deallocate(work_g1)
1068 :
1069 0 : if (my_task == master_task) then
1070 0 : call ice_close_nc(fid_grid)
1071 0 : call ice_close_nc(fid_kmt)
1072 : endif
1073 : #else
1074 : call abort_ice(subname//'ERROR: USE_NETCDF cpp not defined', &
1075 : file=__FILE__, line=__LINE__)
1076 : #endif
1077 :
1078 0 : end subroutine popgrid_nc
1079 :
1080 : #ifdef CESMCOUPLED
1081 : !=======================================================================
1082 :
1083 : ! Read in kmt file that matches CAM lat-lon grid and has single column
1084 : ! functionality
1085 : ! author: Mariana Vertenstein
1086 : ! 2007: Elizabeth Hunke upgraded to netcdf90 and cice ncdf calls
1087 :
1088 : subroutine latlongrid
1089 :
1090 : ! use ice_boundary
1091 : use ice_domain_size
1092 : use ice_scam, only : scmlat, scmlon, single_column
1093 : use ice_constants, only: c0, c1, p5, p25, &
1094 : field_loc_center, field_type_scalar, radius
1095 : #ifdef USE_NETCDF
1096 : use netcdf
1097 : #endif
1098 :
1099 : integer (kind=int_kind) :: &
1100 : i, j, iblk
1101 :
1102 : integer (kind=int_kind) :: &
1103 : ni, nj, ncid, dimid, varid, ier
1104 :
1105 : type (block) :: &
1106 : this_block ! block information for current block
1107 :
1108 : integer (kind=int_kind) :: &
1109 : ilo,ihi,jlo,jhi ! beginning and end of physical domain
1110 :
1111 : real (kind=dbl_kind) :: &
1112 : closelat, & ! Single-column latitude value ! LCOV_EXCL_LINE
1113 : closelon, & ! Single-column longitude value ! LCOV_EXCL_LINE
1114 : closelatidx, & ! Single-column latitude index to retrieve ! LCOV_EXCL_LINE
1115 : closelonidx ! Single-column longitude index to retrieve
1116 :
1117 : integer (kind=int_kind) :: &
1118 : start(2), & ! Start index to read in ! LCOV_EXCL_LINE
1119 : count(2) ! Number of points to read in
1120 :
1121 : integer (kind=int_kind) :: &
1122 : start3(3), & ! Start index to read in ! LCOV_EXCL_LINE
1123 : count3(3) ! Number of points to read in
1124 :
1125 : integer (kind=int_kind) :: &
1126 : status ! status flag
1127 :
1128 : real (kind=dbl_kind), allocatable :: &
1129 : lats(:),lons(:),pos_lons(:), glob_grid(:,:) ! temporaries
1130 :
1131 : real (kind=dbl_kind) :: &
1132 : pos_scmlon,& ! temporary ! LCOV_EXCL_LINE
1133 : pi, & ! LCOV_EXCL_LINE
1134 : puny, & ! LCOV_EXCL_LINE
1135 : scamdata ! temporary
1136 :
1137 : character(len=*), parameter :: subname = '(lonlatgrid)'
1138 :
1139 : #ifdef USE_NETCDF
1140 : !-----------------------------------------------------------------
1141 : ! - kmt file is actually clm fractional land file
1142 : ! - Determine consistency of dimensions
1143 : ! - Read in lon/lat centers in degrees from kmt file
1144 : ! - Read in ocean from "kmt" file (1 for ocean, 0 for land)
1145 : !-----------------------------------------------------------------
1146 :
1147 : call icepack_query_parameters(pi_out=pi, puny_out=puny)
1148 : call icepack_warnings_flush(nu_diag)
1149 : if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
1150 : file=__FILE__, line=__LINE__)
1151 :
1152 : ! Determine dimension of domain file and check for consistency
1153 :
1154 : if (my_task == master_task) then
1155 : call ice_open_nc(kmt_file, ncid)
1156 :
1157 : status = nf90_inq_dimid (ncid, 'ni', dimid)
1158 : status = nf90_inquire_dimension(ncid, dimid, len=ni)
1159 : status = nf90_inq_dimid (ncid, 'nj', dimid)
1160 : status = nf90_inquire_dimension(ncid, dimid, len=nj)
1161 : end if
1162 :
1163 : ! Determine start/count to read in for either single column or global lat-lon grid
1164 : ! If single_column, then assume that only master_task is used since there is only one task
1165 :
1166 : if (single_column) then
1167 : ! Check for consistency
1168 : if (my_task == master_task) then
1169 : if ((nx_global /= 1).or. (ny_global /= 1)) then
1170 : write(nu_diag,*) 'Because you have selected the column model flag'
1171 : write(nu_diag,*) 'Please set nx_global=ny_global=1 in file'
1172 : write(nu_diag,*) 'ice_domain_size.F and recompile'
1173 : call abort_ice (subname//'ERROR: check nx_global, ny_global')
1174 : endif
1175 : end if
1176 :
1177 : ! Read in domain file for single column
1178 : allocate(lats(nj))
1179 : allocate(lons(ni))
1180 : allocate(pos_lons(ni))
1181 : allocate(glob_grid(ni,nj))
1182 :
1183 : start3=(/1,1,1/)
1184 : count3=(/ni,nj,1/)
1185 : status = nf90_inq_varid(ncid, 'xc' , varid)
1186 : if (status /= nf90_noerr) call abort_ice (subname//' inq_varid xc')
1187 : status = nf90_get_var(ncid, varid, glob_grid, start3, count3)
1188 : if (status /= nf90_noerr) call abort_ice (subname//' get_var xc')
1189 : do i = 1,ni
1190 : lons(i) = glob_grid(i,1)
1191 : end do
1192 :
1193 : status = nf90_inq_varid(ncid, 'yc' , varid)
1194 : if (status /= nf90_noerr) call abort_ice (subname//' inq_varid yc')
1195 : status = nf90_get_var(ncid, varid, glob_grid, start3, count3)
1196 : if (status /= nf90_noerr) call abort_ice (subname//' get_var yc')
1197 : do j = 1,nj
1198 : lats(j) = glob_grid(1,j)
1199 : end do
1200 :
1201 : ! convert lons array and scmlon to 0,360 and find index of value closest to 0
1202 : ! and obtain single-column longitude/latitude indices to retrieve
1203 :
1204 : pos_lons(:)= mod(lons(:) + 360._dbl_kind,360._dbl_kind)
1205 : pos_scmlon = mod(scmlon + 360._dbl_kind,360._dbl_kind)
1206 : start(1) = (MINLOC(abs(pos_lons-pos_scmlon),dim=1))
1207 : start(2) = (MINLOC(abs(lats -scmlat ),dim=1))
1208 :
1209 : deallocate(lats)
1210 : deallocate(lons)
1211 : deallocate(pos_lons)
1212 : deallocate(glob_grid)
1213 :
1214 : status = nf90_inq_varid(ncid, 'xc' , varid)
1215 : if (status /= nf90_noerr) call abort_ice (subname//' inq_varid xc')
1216 : status = nf90_get_var(ncid, varid, scamdata, start)
1217 : if (status /= nf90_noerr) call abort_ice (subname//' get_var xc')
1218 : TLON = scamdata
1219 : status = nf90_inq_varid(ncid, 'yc' , varid)
1220 : if (status /= nf90_noerr) call abort_ice (subname//' inq_varid yc')
1221 : status = nf90_get_var(ncid, varid, scamdata, start)
1222 : if (status /= nf90_noerr) call abort_ice (subname//' get_var yc')
1223 : TLAT = scamdata
1224 : status = nf90_inq_varid(ncid, 'area' , varid)
1225 : if (status /= nf90_noerr) call abort_ice (subname//' inq_varid area')
1226 : status = nf90_get_var(ncid, varid, scamdata, start)
1227 : if (status /= nf90_noerr) call abort_ice (subname//' get_var are')
1228 : tarea = scamdata
1229 : status = nf90_inq_varid(ncid, 'mask' , varid)
1230 : if (status /= nf90_noerr) call abort_ice (subname//' inq_varid mask')
1231 : status = nf90_get_var(ncid, varid, scamdata, start)
1232 : if (status /= nf90_noerr) call abort_ice (subname//' get_var mask')
1233 : hm = scamdata
1234 : status = nf90_inq_varid(ncid, 'frac' , varid)
1235 : if (status /= nf90_noerr) call abort_ice (subname//' inq_varid frac')
1236 : status = nf90_get_var(ncid, varid, scamdata, start)
1237 : if (status /= nf90_noerr) call abort_ice (subname//' get_var frac')
1238 : ocn_gridcell_frac = scamdata
1239 : else
1240 : ! Check for consistency
1241 : if (my_task == master_task) then
1242 : if (nx_global /= ni .and. ny_global /= nj) then
1243 : write(nu_diag,*) 'latlongrid: ni,nj = ',ni,nj
1244 : write(nu_diag,*) 'latlongrid: nx_g,ny_g = ',nx_global, ny_global
1245 : call abort_ice (subname//'ERROR: ni,nj not equal to nx_global,ny_global')
1246 : end if
1247 : end if
1248 :
1249 : ! Read in domain file for global lat-lon grid
1250 : call ice_read_nc(ncid, 1, 'xc' , TLON , diag=.true.)
1251 : call ice_read_nc(ncid, 1, 'yc' , TLAT , diag=.true.)
1252 : call ice_read_nc(ncid, 1, 'area', tarea , diag=.true., &
1253 : field_loc=field_loc_center,field_type=field_type_scalar)
1254 : call ice_read_nc(ncid, 1, 'mask', hm , diag=.true.)
1255 : call ice_read_nc(ncid, 1, 'frac', ocn_gridcell_frac, diag=.true.)
1256 : end if
1257 :
1258 : if (my_task == master_task) then
1259 : call ice_close_nc(ncid)
1260 : end if
1261 :
1262 : !$OMP PARALLEL DO PRIVATE(iblk,this_block,ilo,ihi,jlo,jhi,i,j)
1263 : do iblk = 1, nblocks
1264 : this_block = get_block(blocks_ice(iblk),iblk)
1265 : ilo = this_block%ilo
1266 : ihi = this_block%ihi
1267 : jlo = this_block%jlo
1268 : jhi = this_block%jhi
1269 :
1270 : do j = jlo, jhi
1271 : do i = ilo, ihi
1272 : ! Convert from degrees to radians
1273 : TLON(i,j,iblk) = pi*TLON(i,j,iblk)/180._dbl_kind
1274 :
1275 : ! Convert from degrees to radians
1276 : TLAT(i,j,iblk) = pi*TLAT(i,j,iblk)/180._dbl_kind
1277 :
1278 : ! Convert from radians^2 to m^2
1279 : ! (area in domain file is in radians^2 and tarea is in m^2)
1280 : tarea(i,j,iblk) = tarea(i,j,iblk) * (radius*radius)
1281 : end do
1282 : end do
1283 : end do
1284 : !$OMP END PARALLEL DO
1285 :
1286 : !-----------------------------------------------------------------
1287 : ! Calculate various geometric 2d arrays
1288 : ! The U grid (velocity) is not used when run with sequential CAM
1289 : ! because we only use thermodynamic sea ice. However, ULAT is used
1290 : ! in the default initialization of CICE so we calculate it here as
1291 : ! a "dummy" so that CICE will initialize with ice. If a no ice
1292 : ! initialization is OK (or desired) this can be commented out and
1293 : ! ULAT will remain 0 as specified above. ULAT is located at the
1294 : ! NE corner of the grid cell, TLAT at the center, so here ULAT is
1295 : ! hacked by adding half the latitudinal spacing (in radians) to
1296 : ! TLAT.
1297 : !-----------------------------------------------------------------
1298 :
1299 : !$OMP PARALLEL DO PRIVATE(iblk,this_block,ilo,ihi,jlo,jhi,i,j)
1300 : do iblk = 1, nblocks
1301 : this_block = get_block(blocks_ice(iblk),iblk)
1302 : ilo = this_block%ilo
1303 : ihi = this_block%ihi
1304 : jlo = this_block%jlo
1305 : jhi = this_block%jhi
1306 :
1307 : do j = jlo, jhi
1308 : do i = ilo, ihi
1309 :
1310 : if (ny_global == 1) then
1311 : uarea(i,j,iblk) = tarea(i,j, iblk)
1312 : else
1313 : uarea(i,j,iblk) = p25* &
1314 : (tarea(i,j, iblk) + tarea(i+1,j, iblk) & ! LCOV_EXCL_LINE
1315 : + tarea(i,j+1,iblk) + tarea(i+1,j+1,iblk))
1316 : endif
1317 : tarear(i,j,iblk) = c1/tarea(i,j,iblk)
1318 : uarear(i,j,iblk) = c1/uarea(i,j,iblk)
1319 :
1320 : if (single_column) then
1321 : ULAT (i,j,iblk) = TLAT(i,j,iblk)+(pi/nj)
1322 : else
1323 : if (ny_global == 1) then
1324 : ULAT (i,j,iblk) = TLAT(i,j,iblk)
1325 : else
1326 : ULAT (i,j,iblk) = TLAT(i,j,iblk)+(pi/ny_global)
1327 : endif
1328 : endif
1329 : ULON (i,j,iblk) = c0
1330 : NLON (i,j,iblk) = c0
1331 : NLAT (i,j,iblk) = c0
1332 : ELON (i,j,iblk) = c0
1333 : ELAT (i,j,iblk) = c0
1334 : ANGLE (i,j,iblk) = c0
1335 :
1336 : ANGLET(i,j,iblk) = c0
1337 : HTN (i,j,iblk) = 1.e36_dbl_kind
1338 : HTE (i,j,iblk) = 1.e36_dbl_kind
1339 : dxT (i,j,iblk) = 1.e36_dbl_kind
1340 : dyT (i,j,iblk) = 1.e36_dbl_kind
1341 : dxU (i,j,iblk) = 1.e36_dbl_kind
1342 : dyU (i,j,iblk) = 1.e36_dbl_kind
1343 : dxN (i,j,iblk) = 1.e36_dbl_kind
1344 : dyN (i,j,iblk) = 1.e36_dbl_kind
1345 : dxE (i,j,iblk) = 1.e36_dbl_kind
1346 : dyE (i,j,iblk) = 1.e36_dbl_kind
1347 : dxhy (i,j,iblk) = 1.e36_dbl_kind
1348 : dyhx (i,j,iblk) = 1.e36_dbl_kind
1349 : cyp (i,j,iblk) = 1.e36_dbl_kind
1350 : cxp (i,j,iblk) = 1.e36_dbl_kind
1351 : cym (i,j,iblk) = 1.e36_dbl_kind
1352 : cxm (i,j,iblk) = 1.e36_dbl_kind
1353 : enddo
1354 : enddo
1355 : enddo
1356 : !$OMP END PARALLEL DO
1357 :
1358 : call makemask
1359 : #else
1360 : call abort_ice(subname//'ERROR: USE_NETCDF cpp not defined', &
1361 : file=__FILE__, line=__LINE__)
1362 : #endif
1363 :
1364 : end subroutine latlongrid
1365 : #endif
1366 : !=======================================================================
1367 :
1368 : ! Regular rectangular grid and mask
1369 : !
1370 : ! author: Elizabeth C. Hunke, LANL
1371 :
1372 16 : subroutine rectgrid
1373 :
1374 : use ice_constants, only: c0, c1, c2, radius, cm_to_m, &
1375 : field_loc_center, field_loc_NEcorner, field_type_scalar
1376 : use ice_domain, only: close_boundaries
1377 :
1378 : integer (kind=int_kind) :: &
1379 : i, j, & ! LCOV_EXCL_LINE
1380 : imid, jmid
1381 :
1382 : real (kind=dbl_kind) :: &
1383 : length, & ! LCOV_EXCL_LINE
1384 0 : rad_to_deg
1385 :
1386 : real (kind=dbl_kind), dimension(:,:), allocatable :: &
1387 16 : work_g1
1388 :
1389 : character(len=*), parameter :: subname = '(rectgrid)'
1390 :
1391 : !-----------------------------------------------------------------
1392 : ! Calculate various geometric 2d arrays
1393 : !-----------------------------------------------------------------
1394 :
1395 16 : call icepack_query_parameters(rad_to_deg_out=rad_to_deg)
1396 16 : call icepack_warnings_flush(nu_diag)
1397 16 : if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
1398 0 : file=__FILE__, line=__LINE__)
1399 :
1400 38128 : hm (:,:,:) = c0
1401 38128 : kmt(:,:,:) = c0
1402 38128 : angle(:,:,:) = c0 ! "square with the world"
1403 :
1404 16 : allocate(work_g1(nx_global,ny_global))
1405 :
1406 16 : if (scale_dxdy) then
1407 : ! scale grid spacing from center outward.
1408 : ! this different than original method in it
1409 : ! needs to define grid spacing before lat/lon.
1410 : ! original rectgrid defines latlon first
1411 0 : call rectgrid_scale_dxdy
1412 : else
1413 : ! rectgrid no grid spacing.
1414 : ! original method with addition to use namelist lat/lon reference
1415 :
1416 16 : if (my_task == master_task) then
1417 33026 : work_g1 = c0
1418 2 : length = dxrect*cm_to_m/radius*rad_to_deg
1419 :
1420 258 : work_g1(1,:) = lonrefrect ! reference lon from namelist
1421 :
1422 258 : do j = 1, ny_global
1423 32770 : do i = 2, nx_global
1424 32768 : work_g1(i,j) = work_g1(i-1,j) + length ! ULON
1425 : enddo
1426 : enddo
1427 33026 : work_g1(:,:) = work_g1(:,:) / rad_to_deg
1428 : endif
1429 : call scatter_global(ULON, work_g1, master_task, distrb_info, &
1430 16 : field_loc_NEcorner, field_type_scalar)
1431 : call ice_HaloExtrapolate(ULON, distrb_info, &
1432 16 : ew_boundary_type, ns_boundary_type)
1433 :
1434 16 : if (my_task == master_task) then
1435 33026 : work_g1 = c0
1436 2 : length = dyrect*cm_to_m/radius*rad_to_deg
1437 :
1438 258 : work_g1(:,1) = latrefrect ! reference latitude from namelist
1439 :
1440 258 : do i = 1, nx_global
1441 32770 : do j = 2, ny_global
1442 32768 : work_g1(i,j) = work_g1(i,j-1) + length ! ULAT
1443 : enddo
1444 : enddo
1445 33026 : work_g1(:,:) = work_g1(:,:) / rad_to_deg
1446 : endif
1447 : call scatter_global(ULAT, work_g1, master_task, distrb_info, &
1448 16 : field_loc_NEcorner, field_type_scalar)
1449 : call ice_HaloExtrapolate(ULAT, distrb_info, &
1450 16 : ew_boundary_type, ns_boundary_type)
1451 :
1452 16 : if (my_task == master_task) then
1453 258 : do j = 1, ny_global
1454 33026 : do i = 1, nx_global
1455 33024 : work_g1(i,j) = dxrect ! HTN
1456 : enddo
1457 : enddo
1458 : endif
1459 16 : call primary_grid_lengths_HTN(work_g1) ! dxU, dxT, dxN, dxE
1460 :
1461 16 : if (my_task == master_task) then
1462 258 : do j = 1, ny_global
1463 33026 : do i = 1, nx_global
1464 33024 : work_g1(i,j) = dyrect ! HTE
1465 : enddo
1466 : enddo
1467 : endif
1468 16 : call primary_grid_lengths_HTE(work_g1) ! dyU, dyT, dyN, dyE
1469 :
1470 : endif ! scale_dxdy
1471 :
1472 : !-----------------------------------------------------------------
1473 : ! Construct T-cell land mask
1474 : ! Keyed on ew_boundary_type; ns_boundary_type should be 'open'.
1475 : !-----------------------------------------------------------------
1476 :
1477 16 : if (my_task == master_task) then
1478 33026 : work_g1(:,:) = c0 ! initialize hm as land
1479 :
1480 2 : if (trim(kmt_type) == 'boxislands') then
1481 :
1482 0 : call grid_boxislands_kmt(work_g1)
1483 :
1484 2 : elseif (trim(kmt_type) == 'channel') then
1485 :
1486 0 : do j = 3,ny_global-2 ! closed top and bottom
1487 0 : do i = 1,nx_global ! open sides
1488 0 : work_g1(i,j) = c1 ! NOTE nx_global > 5
1489 : enddo
1490 : enddo
1491 :
1492 2 : elseif (trim(kmt_type) == 'channel_oneeast') then
1493 :
1494 0 : do j = ny_global/2,ny_global/2 ! one channel wide
1495 0 : do i = 1,nx_global ! open sides
1496 0 : work_g1(i,j) = c1 ! NOTE nx_global > 5
1497 : enddo
1498 : enddo
1499 :
1500 2 : elseif (trim(kmt_type) == 'channel_onenorth') then
1501 :
1502 0 : do j = 1,ny_global ! open sides
1503 0 : do i = nx_global/2,nx_global/2 ! one channel wide
1504 0 : work_g1(i,j) = c1 ! NOTE nx_global > 5
1505 : enddo
1506 : enddo
1507 :
1508 2 : elseif (trim(kmt_type) == 'wall') then
1509 :
1510 0 : do j = 1,ny_global ! open except
1511 0 : do i = 1,nx_global-2 ! closed east edge
1512 0 : work_g1(i,j) = c1
1513 : enddo
1514 : enddo
1515 :
1516 2 : elseif (trim(kmt_type) == 'default') then
1517 :
1518 : ! land in the upper left and lower right corners,
1519 : ! otherwise open boundaries
1520 2 : imid = nint(aint(real(nx_global)/c2))
1521 2 : jmid = nint(aint(real(ny_global)/c2))
1522 :
1523 250 : do j = 3,ny_global-2
1524 31002 : do i = 3,nx_global-2
1525 31000 : work_g1(i,j) = c1 ! open central domain
1526 : enddo
1527 : enddo
1528 :
1529 2 : if (nx_global > 5 .and. ny_global > 5) then
1530 :
1531 134 : do j = 1, jmid+2
1532 8846 : do i = 1, imid+2
1533 8844 : work_g1(i,j) = c1 ! open lower left corner
1534 : enddo
1535 : enddo
1536 :
1537 136 : do j = max(jmid-2,1), ny_global
1538 9114 : do i = max(imid-2,1), nx_global
1539 9112 : work_g1(i,j) = c1 ! open upper right corner
1540 : enddo
1541 : enddo
1542 :
1543 : endif ! > 5x5 grid
1544 :
1545 : else
1546 :
1547 0 : call abort_ice(subname//'ERROR: unknown kmt_type '//trim(kmt_type))
1548 :
1549 : endif ! kmt_type
1550 :
1551 2 : if (close_boundaries) then
1552 0 : work_g1(:, 1:2) = c0
1553 0 : work_g1(:, ny_global-1:ny_global) = c0
1554 0 : work_g1(1:2, :) = c0
1555 0 : work_g1(nx_global-1:nx_global, :) = c0
1556 : endif
1557 :
1558 : endif
1559 :
1560 : call scatter_global(hm, work_g1, master_task, distrb_info, &
1561 16 : field_loc_center, field_type_scalar)
1562 :
1563 16 : deallocate(work_g1)
1564 :
1565 32 : end subroutine rectgrid
1566 :
1567 : !=======================================================================
1568 :
1569 0 : subroutine rectgrid_scale_dxdy
1570 :
1571 : ! generate a variable spaced rectangluar grid.
1572 : ! extend spacing from center of grid outward.
1573 : use ice_constants, only: c0, c1, c2, radius, cm_to_m, &
1574 : field_loc_center, field_loc_NEcorner, field_type_scalar
1575 :
1576 : integer (kind=int_kind) :: &
1577 : i, j, iblk, & ! LCOV_EXCL_LINE
1578 : imid, jmid, & ! LCOV_EXCL_LINE
1579 : center1, center2 ! array centers for expanding dx, dy
1580 :
1581 : real (kind=dbl_kind) :: &
1582 : length, & ! LCOV_EXCL_LINE
1583 0 : rad_to_deg
1584 :
1585 : real (kind=dbl_kind), dimension(:,:), allocatable :: &
1586 0 : work_g1
1587 :
1588 : character(len=*), parameter :: subname = '(rectgrid_scale_dxdy)'
1589 :
1590 0 : call icepack_query_parameters(rad_to_deg_out=rad_to_deg)
1591 :
1592 0 : allocate(work_g1(nx_global,ny_global))
1593 :
1594 : ! determine dx spacing
1595 : ! strategy: initialize with dxrect.
1596 : ! if want to scale the grid, work from center outwards,
1597 : ! multplying neighbor cell by scale factor.
1598 : ! this assumes dx varies in x direction only.
1599 : ! (i.e, dx is the same across same y location)
1600 0 : if (my_task == master_task) then
1601 :
1602 : ! initialize with initial dxrect
1603 0 : work_g1(:,:) = dxrect
1604 :
1605 : ! check if nx is even or odd
1606 : ! if even, middle 2 columns are center
1607 : ! of odd, middle 1 column is center
1608 0 : if (mod(nx_global,2) == 0) then ! nx_global is even
1609 :
1610 : ! with even number of x locatons,
1611 : ! the center two y columns are center
1612 0 : center1 = nx_global/2 ! integer math
1613 0 : center2 = center1 + 1 ! integer math
1614 :
1615 : else ! nx_global = odd
1616 : ! only one center index. set center2=center1
1617 0 : center1 = ceiling(real(nx_global/2),int_kind)
1618 0 : center2 = center1
1619 : endif
1620 :
1621 : ! note loop over only half the x grid points (center1)-1
1622 : ! working from the center outward.
1623 0 : do j = 1, ny_global
1624 0 : do i = 1, center1-1
1625 : ! work from center1 to left
1626 0 : work_g1(center1-i,j) = dxscale*work_g1(center1-i+1,j)
1627 :
1628 : ! work from center2 to right
1629 0 : work_g1(center2+i,j) = dxscale*work_g1(center2+i-1,j)
1630 : enddo ! i
1631 : enddo ! j
1632 :
1633 : endif ! my_task == master_task
1634 :
1635 :
1636 : ! note work_g1 is converted to meters in primary_grid_lengths_HTN
1637 0 : call primary_grid_lengths_HTN(work_g1) ! dxU, dxT, dxN, dxE
1638 :
1639 : ! make ULON array
1640 0 : if (my_task == master_task) then
1641 :
1642 : ! make first column reference lon in radians.
1643 : ! the remaining work_g1 is still dx in meters
1644 0 : work_g1(1,:) = lonrefrect/rad_to_deg ! radians
1645 :
1646 : ! loop over remaining points and add spacing to successive
1647 : ! x locations
1648 0 : do j = 1, ny_global
1649 0 : do i = 2, nx_global ! start from i=2. i=1 is lonrefrect
1650 0 : length = work_g1(i,j)/radius ! grid spacing in radians
1651 0 : work_g1(i,j) = work_g1(i-1,j) + length ! ULON
1652 : enddo ! i
1653 : enddo ! j
1654 : endif ! mytask == master_task
1655 : call scatter_global(ULON, work_g1, master_task, distrb_info, &
1656 0 : field_loc_NEcorner, field_type_scalar)
1657 : call ice_HaloExtrapolate(ULON, distrb_info, &
1658 0 : ew_boundary_type, ns_boundary_type)
1659 :
1660 : ! determine dy spacing
1661 : ! strategy: initialize with dyrect.
1662 : ! if want to scale the grid, work from center outwards,
1663 : ! multplying neighbor cell by scale factor.
1664 : ! this assumes dy varies in y direction only.
1665 : ! (i.e, dy is the same across same x location)
1666 0 : if (my_task == master_task) then
1667 :
1668 : ! initialize with initial dxrect
1669 0 : work_g1(:,:) = dyrect
1670 :
1671 : ! check if ny is even or odd
1672 : ! if even, middle 2 rows are center
1673 : ! of odd, middle 1 row is center
1674 0 : if (mod(ny_global,2) == 0) then ! ny_global is even
1675 :
1676 : ! with even number of x locatons,
1677 : ! the center two y columns are center
1678 0 : center1 = ny_global/2 ! integer math
1679 0 : center2 = center1 + 1 ! integer math
1680 :
1681 : else ! ny_global = odd
1682 : ! only one center index. set center2=center1
1683 0 : center1 = ceiling(real(ny_global/2),int_kind)
1684 0 : center2 = center1
1685 : endif
1686 :
1687 : ! note loop over only half the y grid points (center1)-1
1688 : ! working from the center outward.
1689 0 : do i = 1, nx_global
1690 0 : do j = 1, center1-1
1691 : ! work from center1 to bottom
1692 0 : work_g1(i,center1-j) = dyscale*work_g1(i,center1-j+1)
1693 :
1694 : ! work from center2 to top
1695 0 : work_g1(i,center2+j) = dyscale*work_g1(i,center2+j-1)
1696 : enddo ! i
1697 : enddo ! j
1698 : endif ! mytask == master_task
1699 : ! note work_g1 is converted to meters primary_grid_lengths_HTE
1700 0 : call primary_grid_lengths_HTE(work_g1) ! dyU, dyT, dyN, dyE
1701 :
1702 : ! make ULAT array
1703 0 : if (my_task == master_task) then
1704 :
1705 : ! make first row reference lat in radians.
1706 : ! the remaining work_g1 is still dy in meters
1707 0 : work_g1(:,1) = latrefrect/rad_to_deg ! radians
1708 :
1709 :
1710 : ! loop over remaining points and add spacing to successive
1711 : ! x locations
1712 0 : do j = 2, ny_global ! start from j=2. j=1 is latrefrect
1713 0 : do i = 1, nx_global
1714 0 : length = work_g1(i,j)/radius ! grid spacing in radians
1715 0 : work_g1(i,j) = work_g1(i,j-1) + length ! ULAT
1716 : enddo ! i
1717 : enddo ! j
1718 : endif ! mytask == master_task
1719 : call scatter_global(ULAT, work_g1, master_task, distrb_info, &
1720 0 : field_loc_NEcorner, field_type_scalar)
1721 : call ice_HaloExtrapolate(ULAT, distrb_info, &
1722 0 : ew_boundary_type, ns_boundary_type)
1723 :
1724 :
1725 0 : deallocate(work_g1)
1726 :
1727 0 : end subroutine rectgrid_scale_dxdy
1728 :
1729 : !=======================================================================
1730 :
1731 : ! Complex land mask for testing box cases
1732 : ! Requires nx_global, ny_global > 20
1733 : ! Assumes work array has been initialized to 1 (ocean) and north and
1734 : ! south land boundaries have been applied (ew_boundary_type='cyclic')
1735 :
1736 0 : subroutine grid_boxislands_kmt (work)
1737 :
1738 : use ice_constants, only: c0, c1, c20
1739 :
1740 : real (kind=dbl_kind), dimension(:,:), intent(inout) :: work
1741 :
1742 : integer (kind=int_kind) :: &
1743 : i, j, k, & ! indices ! LCOV_EXCL_LINE
1744 : nxb, nyb ! convenient cell-block sizes for building the mask
1745 :
1746 : character(len=*), parameter :: subname = '(grid_boxislands_kmt)'
1747 :
1748 : ! number of cells in 5% of global grid x and y lengths
1749 0 : nxb = int(real(nx_global, dbl_kind) / c20, int_kind)
1750 0 : nyb = int(real(ny_global, dbl_kind) / c20, int_kind)
1751 :
1752 0 : if (nxb < 1 .or. nyb < 1) &
1753 0 : call abort_ice(subname//'ERROR: requires larger grid size')
1754 :
1755 : ! initialize work area as all ocean (c1).
1756 0 : work(:,:) = c1
1757 :
1758 : ! now add land points (c0)
1759 : ! northeast triangle
1760 0 : k = 0
1761 0 : do j = ny_global, ny_global-3*nyb, -1
1762 0 : k = k+1
1763 0 : do i = nx_global-3*nxb+k, nx_global
1764 0 : work(i,j) = c0
1765 : enddo
1766 : enddo
1767 :
1768 : ! northwest docks
1769 0 : do j = ny_global-3*nyb, ny_global
1770 0 : do i = 1, 1
1771 0 : work(i,j) = c0
1772 : enddo
1773 : enddo
1774 0 : do i = 1, 2*nxb
1775 0 : do j = ny_global-3*nyb, ny_global-nyb-2
1776 0 : work(i,j) = c0
1777 : enddo
1778 0 : do j = ny_global-nyb, ny_global-nyb+1
1779 0 : work(i,j) = c0
1780 : enddo
1781 : enddo
1782 :
1783 : ! southwest docks
1784 0 : do j = 2*nyb, 3*nyb
1785 0 : do i = 1, 1
1786 0 : work(i,j) = c0
1787 : enddo
1788 : enddo
1789 0 : do j = 1, 2*nyb
1790 0 : do i = 2, nxb
1791 0 : work(i,j) = c0
1792 : enddo
1793 0 : do i = 2*nxb-1, 2*nxb
1794 0 : work(i,j) = c0
1795 : enddo
1796 0 : do i = 2*nxb+2,4*nxb
1797 0 : work(i,j) = c0
1798 : enddo
1799 : enddo
1800 :
1801 : ! tiny island
1802 0 : do j = 14*nyb, 14*nyb+1
1803 0 : do i = 14*nxb, 14*nxb+1
1804 0 : work(i,j) = c0
1805 : enddo
1806 : enddo
1807 :
1808 : ! X islands
1809 : ! left triangle
1810 0 : k = 0
1811 0 : do i = 2*nxb, 4*nxb
1812 0 : k=k+1
1813 0 : do j = 10*nyb+k, 14*nyb-k
1814 0 : work(i,j) = c0
1815 : enddo
1816 : enddo
1817 : ! upper triangle
1818 0 : k = 0
1819 0 : do j = 14*nyb, 12*nyb, -1
1820 0 : k=k+1
1821 0 : do i = 2*nxb+2+k, 6*nxb-2-k
1822 0 : work(i,j) = c0
1823 : enddo
1824 : enddo
1825 : ! diagonal
1826 0 : k = 0
1827 0 : do j = 10*nyb, 14*nyb
1828 0 : k=k+1
1829 0 : do i = 2*nxb+4+k, 2*nxb+6+k
1830 0 : work(i,j) = c0
1831 : enddo
1832 : enddo
1833 : ! lower right triangle
1834 0 : k = 0
1835 0 : do j = 12*nyb, 10*nyb, -1
1836 0 : k=k+1
1837 0 : do i = 5*nxb+k, 8*nxb
1838 0 : work(i,j) = c0
1839 : enddo
1840 : enddo
1841 :
1842 : ! bar islands
1843 0 : do i = 10*nxb, 16*nxb
1844 0 : do j = 4*nyb, 5*nyb
1845 0 : work(i,j) = c0
1846 : enddo
1847 0 : do j = 6*nyb+2, 8*nyb
1848 0 : work(i,j) = c0
1849 : enddo
1850 0 : do j = 8*nyb+2, 8*nyb+3
1851 0 : work(i,j) = c0
1852 : enddo
1853 : enddo
1854 :
1855 0 : end subroutine grid_boxislands_kmt
1856 :
1857 : !=======================================================================
1858 :
1859 : ! CPOM displaced pole grid and land mask. \\
1860 : ! Grid record number, field and units are: \\
1861 : ! (1) ULAT (degrees) \\
1862 : ! (2) ULON (degrees) \\
1863 : ! (3) HTN (m) \\
1864 : ! (4) HTE (m) \\
1865 : ! (7) ANGLE (radians) \\
1866 : !
1867 : ! Land mask record number and field is (1) KMT.
1868 : !
1869 : ! author: Adrian K. Turner, CPOM, UCL, 09/08/06
1870 :
1871 0 : subroutine cpomgrid
1872 :
1873 : use ice_blocks, only: nx_block, ny_block
1874 : use ice_constants, only: c0, c1, m_to_cm, &
1875 : field_loc_NEcorner, field_type_scalar
1876 : use ice_domain_size, only: max_blocks
1877 :
1878 : integer (kind=int_kind) :: &
1879 : i, j, iblk, & ! LCOV_EXCL_LINE
1880 : ilo,ihi,jlo,jhi ! beginning and end of physical domain
1881 :
1882 : logical (kind=log_kind) :: diag
1883 :
1884 : real (kind=dbl_kind), dimension(:,:), allocatable :: &
1885 0 : work_g1
1886 :
1887 : real (kind=dbl_kind), dimension (nx_block,ny_block,max_blocks) :: &
1888 0 : work1
1889 :
1890 : real (kind=dbl_kind) :: &
1891 0 : rad_to_deg
1892 :
1893 : type (block) :: &
1894 : this_block ! block information for current block
1895 :
1896 : character(len=*), parameter :: subname = '(cpomgrid)'
1897 :
1898 0 : call icepack_query_parameters(rad_to_deg_out=rad_to_deg)
1899 0 : call icepack_warnings_flush(nu_diag)
1900 0 : if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
1901 0 : file=__FILE__, line=__LINE__)
1902 :
1903 0 : call ice_open(nu_grid,grid_file,64)
1904 0 : call ice_open(nu_kmt,kmt_file,32)
1905 :
1906 0 : diag = .true. ! write diagnostic info
1907 :
1908 : ! topography
1909 0 : call ice_read(nu_kmt,1,work1,'ida4',diag)
1910 :
1911 0 : hm (:,:,:) = c0
1912 0 : kmt(:,:,:) = c0
1913 0 : !$OMP PARALLEL DO PRIVATE(iblk,i,j,ilo,ihi,jlo,jhi,this_block)
1914 0 : do iblk = 1, nblocks
1915 0 : this_block = get_block(blocks_ice(iblk),iblk)
1916 0 : ilo = this_block%ilo
1917 0 : ihi = this_block%ihi
1918 0 : jlo = this_block%jlo
1919 0 : jhi = this_block%jhi
1920 :
1921 0 : do j = jlo, jhi
1922 0 : do i = ilo, ihi
1923 0 : kmt(i,j,iblk) = work1(i,j,iblk)
1924 0 : if (kmt(i,j,iblk) >= c1) hm(i,j,iblk) = c1
1925 : enddo
1926 : enddo
1927 : enddo
1928 : !$OMP END PARALLEL DO
1929 :
1930 0 : allocate(work_g1(nx_global,ny_global))
1931 :
1932 : ! lat, lon, cell dimensions, angles
1933 0 : call ice_read_global(nu_grid,1,work_g1, 'rda8',diag)
1934 : call scatter_global(ULAT, work_g1, master_task, distrb_info, &
1935 0 : field_loc_NEcorner, field_type_scalar)
1936 :
1937 0 : call ice_read_global(nu_grid,2,work_g1, 'rda8',diag)
1938 : call scatter_global(ULON, work_g1, master_task, distrb_info, &
1939 0 : field_loc_NEcorner, field_type_scalar)
1940 :
1941 0 : call ice_read_global(nu_grid,3,work_g1, 'rda8',diag)
1942 0 : work_g1 = work_g1 * m_to_cm
1943 0 : call primary_grid_lengths_HTN(work_g1) ! dxU, dxT, dxN, dxE
1944 :
1945 0 : call ice_read_global(nu_grid,4,work_g1, 'rda8',diag)
1946 0 : work_g1 = work_g1 * m_to_cm
1947 0 : call primary_grid_lengths_HTE(work_g1) ! dyU, dyT, dyN, dyE
1948 :
1949 0 : call ice_read_global(nu_grid,7,work_g1,'rda8',diag)
1950 : call scatter_global(ANGLE, work_g1, master_task, distrb_info, &
1951 0 : field_loc_NEcorner, field_type_scalar)
1952 :
1953 : ! fix units
1954 0 : ULAT = ULAT / rad_to_deg
1955 0 : ULON = ULON / rad_to_deg
1956 :
1957 0 : deallocate(work_g1)
1958 :
1959 0 : if (my_task == master_task) then
1960 0 : close (nu_grid)
1961 0 : close (nu_kmt)
1962 : endif
1963 :
1964 0 : write(nu_diag,*) "min/max HTN: ", minval(HTN), maxval(HTN)
1965 0 : write(nu_diag,*) "min/max HTE: ", minval(HTE), maxval(HTE)
1966 :
1967 0 : end subroutine cpomgrid
1968 :
1969 : !=======================================================================
1970 :
1971 : ! Calculate dxU and dxT from HTN on the global grid, to preserve
1972 : ! ghost cell and/or land values that might otherwise be lost. Scatter
1973 : ! dxU, dxT and HTN to all processors.
1974 : !
1975 : ! author: Elizabeth C. Hunke, LANL
1976 :
1977 37 : subroutine primary_grid_lengths_HTN(work_g)
1978 :
1979 : use ice_constants, only: p25, p5, c2, cm_to_m, &
1980 : field_loc_center, field_loc_NEcorner, & ! LCOV_EXCL_LINE
1981 : field_loc_Nface, field_type_scalar
1982 :
1983 : real (kind=dbl_kind), dimension(:,:) :: work_g ! global array holding HTN
1984 :
1985 : ! local variables
1986 :
1987 : integer (kind=int_kind) :: &
1988 : i, j, & ! LCOV_EXCL_LINE
1989 : ip1 ! i+1
1990 :
1991 : real (kind=dbl_kind), dimension(:,:), allocatable :: &
1992 37 : work_g2
1993 :
1994 : character(len=*), parameter :: subname = '(primary_grid_lengths_HTN)'
1995 :
1996 37 : if (my_task == master_task) then
1997 7 : allocate(work_g2(nx_global,ny_global))
1998 : else
1999 30 : allocate(work_g2(1,1))
2000 : endif
2001 :
2002 : ! HTN, dxU = average of 2 neighbor HTNs in i
2003 :
2004 37 : if (my_task == master_task) then
2005 843 : do j = 1, ny_global
2006 91611 : do i = 1, nx_global
2007 91604 : work_g(i,j) = work_g(i,j) * cm_to_m ! HTN
2008 : enddo
2009 : enddo
2010 843 : do j = 1, ny_global
2011 91611 : do i = 1, nx_global
2012 : ! assume cyclic; noncyclic will be handled during scatter
2013 90768 : ip1 = i+1
2014 90768 : if (i == nx_global) ip1 = 1
2015 91604 : work_g2(i,j) = p5*(work_g(i,j) + work_g(ip1,j)) ! dxU
2016 : enddo
2017 : enddo
2018 : endif
2019 37 : if (pgl_global_ext) then
2020 : call primary_grid_lengths_global_ext( &
2021 0 : G_HTN, work_g, ew_boundary_type, ns_boundary_type)
2022 : endif
2023 0 : call scatter_global(HTN, work_g, master_task, distrb_info, &
2024 37 : field_loc_Nface, field_type_scalar)
2025 : call scatter_global(dxU, work_g2, master_task, distrb_info, &
2026 37 : field_loc_NEcorner, field_type_scalar)
2027 :
2028 : ! dxT = average of 2 neighbor HTNs in j
2029 :
2030 37 : if (my_task == master_task) then
2031 836 : do j = 2, ny_global
2032 90848 : do i = 1, nx_global
2033 90841 : work_g2(i,j) = p5*(work_g(i,j) + work_g(i,j-1)) ! dxT
2034 : enddo
2035 : enddo
2036 : ! extrapolate to obtain dxT along j=1
2037 763 : do i = 1, nx_global
2038 763 : work_g2(i,1) = c2*work_g(i,2) - work_g(i,3) ! dxT
2039 : enddo
2040 : endif
2041 : call scatter_global(dxT, work_g2, master_task, distrb_info, &
2042 37 : field_loc_center, field_type_scalar)
2043 :
2044 : ! dxN = HTN
2045 :
2046 111936 : dxN(:,:,:) = HTN(:,:,:) ! dxN
2047 :
2048 : ! dxE = average of 4 surrounding HTNs
2049 :
2050 37 : if (my_task == master_task) then
2051 836 : do j = 2, ny_global
2052 90848 : do i = 1, nx_global
2053 : ! assume cyclic; noncyclic will be handled during scatter
2054 90012 : ip1 = i+1
2055 90012 : if (i == nx_global) ip1 = 1
2056 90841 : work_g2(i,j) = p25*(work_g(i,j)+work_g(ip1,j)+work_g(i,j-1)+work_g(ip1,j-1)) ! dxE
2057 : enddo
2058 : enddo
2059 : ! extrapolate to obtain dxT along j=1
2060 763 : do i = 1, nx_global
2061 : ! assume cyclic; noncyclic will be handled during scatter
2062 756 : ip1 = i+1
2063 756 : if (i == nx_global) ip1 = 1
2064 200 : work_g2(i,1) = p5*(c2*work_g(i ,2) - work_g(i ,3) + &
2065 963 : c2*work_g(ip1,2) - work_g(ip1,3)) ! dxE
2066 : enddo
2067 : endif
2068 : call scatter_global(dxE, work_g2, master_task, distrb_info, &
2069 37 : field_loc_center, field_type_scalar)
2070 :
2071 37 : deallocate(work_g2)
2072 :
2073 37 : end subroutine primary_grid_lengths_HTN
2074 :
2075 : !=======================================================================
2076 : ! Calculate dyU and dyT from HTE on the global grid, to preserve
2077 : ! ghost cell and/or land values that might otherwise be lost. Scatter
2078 : ! dyU, dyT and HTE to all processors.
2079 : !
2080 : ! author: Elizabeth C. Hunke, LANL
2081 :
2082 37 : subroutine primary_grid_lengths_HTE(work_g)
2083 :
2084 : use ice_constants, only: p25, p5, c2, cm_to_m, &
2085 : field_loc_center, field_loc_NEcorner, & ! LCOV_EXCL_LINE
2086 : field_loc_Eface, field_type_scalar
2087 :
2088 : real (kind=dbl_kind), dimension(:,:) :: work_g ! global array holding HTE
2089 :
2090 : ! local variables
2091 :
2092 : integer (kind=int_kind) :: &
2093 : i, j, & ! LCOV_EXCL_LINE
2094 : im1 ! i-1
2095 :
2096 : real (kind=dbl_kind), dimension(:,:), allocatable :: &
2097 37 : work_g2
2098 :
2099 : character(len=*), parameter :: subname = '(primary_grid_lengths_HTE)'
2100 :
2101 37 : if (my_task == master_task) then
2102 7 : allocate(work_g2(nx_global,ny_global))
2103 : else
2104 30 : allocate(work_g2(1,1))
2105 : endif
2106 :
2107 : ! HTE, dyU = average of 2 neighbor HTE in j
2108 :
2109 37 : if (my_task == master_task) then
2110 843 : do j = 1, ny_global
2111 91611 : do i = 1, nx_global
2112 91604 : work_g(i,j) = work_g(i,j) * cm_to_m ! HTE
2113 : enddo
2114 : enddo
2115 836 : do j = 1, ny_global-1
2116 90848 : do i = 1, nx_global
2117 90841 : work_g2(i,j) = p5*(work_g(i,j) + work_g(i,j+1)) ! dyU
2118 : enddo
2119 : enddo
2120 : ! extrapolate to obtain dyU along j=ny_global
2121 7 : if (ny_global > 1) then
2122 763 : do i = 1, nx_global
2123 763 : work_g2(i,ny_global) = c2*work_g(i,ny_global-1) - work_g(i,ny_global-2) ! dyU
2124 : enddo
2125 : endif
2126 : endif
2127 37 : if (pgl_global_ext) then
2128 : call primary_grid_lengths_global_ext( &
2129 0 : G_HTE, work_g, ew_boundary_type, ns_boundary_type)
2130 : endif
2131 0 : call scatter_global(HTE, work_g, master_task, distrb_info, &
2132 37 : field_loc_Eface, field_type_scalar)
2133 : call scatter_global(dyU, work_g2, master_task, distrb_info, &
2134 37 : field_loc_NEcorner, field_type_scalar)
2135 :
2136 : ! dyT = average of 2 neighbor HTE in i
2137 :
2138 37 : if (my_task == master_task) then
2139 843 : do j = 1, ny_global
2140 91611 : do i = 1, nx_global
2141 : ! assume cyclic; noncyclic will be handled during scatter
2142 90768 : im1 = i-1
2143 90768 : if (i == 1) im1 = nx_global
2144 91604 : work_g2(i,j) = p5*(work_g(i,j) + work_g(im1,j)) ! dyT
2145 : enddo
2146 : enddo
2147 : endif
2148 : call scatter_global(dyT, work_g2, master_task, distrb_info, &
2149 37 : field_loc_center, field_type_scalar)
2150 :
2151 : ! dyN = average of 4 neighbor HTEs
2152 :
2153 37 : if (my_task == master_task) then
2154 836 : do j = 1, ny_global-1
2155 90848 : do i = 1, nx_global
2156 : ! assume cyclic; noncyclic will be handled during scatter
2157 90012 : im1 = i-1
2158 90012 : if (i == 1) im1 = nx_global
2159 90841 : work_g2(i,j) = p25*(work_g(i,j) + work_g(im1,j) + work_g(i,j+1) + work_g(im1,j+1)) ! dyN
2160 : enddo
2161 : enddo
2162 : ! extrapolate to obtain dyN along j=ny_global
2163 7 : if (ny_global > 1) then
2164 763 : do i = 1, nx_global
2165 : ! assume cyclic; noncyclic will be handled during scatter
2166 756 : im1 = i-1
2167 756 : if (i == 1) im1 = nx_global
2168 400 : work_g2(i,ny_global) = p5*(c2*work_g(i ,ny_global-1) - work_g(i ,ny_global-2) + &
2169 1163 : c2*work_g(im1,ny_global-1) - work_g(im1,ny_global-2)) ! dyN
2170 : enddo
2171 : endif
2172 : endif
2173 : call scatter_global(dyN, work_g2, master_task, distrb_info, &
2174 37 : field_loc_center, field_type_scalar)
2175 :
2176 : ! dyE = HTE
2177 :
2178 111936 : dyE(:,:,:) = HTE(:,:,:)
2179 :
2180 37 : deallocate(work_g2)
2181 :
2182 37 : end subroutine primary_grid_lengths_HTE
2183 :
2184 : !=======================================================================
2185 :
2186 : ! Sets the boundary values for the T cell land mask (hm) and
2187 : ! makes the logical land masks for T and U cells (tmask, umask)
2188 : ! and N and E cells (nmask, emask).
2189 : ! Also creates hemisphere masks (mask-n northern, mask-s southern)
2190 : !
2191 : ! author: Elizabeth C. Hunke, LANL
2192 :
2193 37 : subroutine makemask
2194 :
2195 : use ice_constants, only: c0, p5, c1p5, &
2196 : field_loc_center, field_loc_NEcorner, field_type_scalar, & ! LCOV_EXCL_LINE
2197 : field_loc_Nface, field_loc_Eface
2198 :
2199 : integer (kind=int_kind) :: &
2200 : i, j, iblk, & ! LCOV_EXCL_LINE
2201 : ilo,ihi,jlo,jhi ! beginning and end of physical domain
2202 :
2203 : real (kind=dbl_kind) :: &
2204 8 : puny
2205 :
2206 : real (kind=dbl_kind), dimension(:,:,:), allocatable :: &
2207 37 : uvmCD
2208 :
2209 : type (block) :: &
2210 : this_block ! block information for current block
2211 :
2212 : character(len=*), parameter :: subname = '(makemask)'
2213 :
2214 37 : call icepack_query_parameters(puny_out=puny)
2215 37 : call icepack_warnings_flush(nu_diag)
2216 37 : if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
2217 0 : file=__FILE__, line=__LINE__)
2218 :
2219 37 : call ice_timer_start(timer_bound)
2220 : call ice_HaloUpdate (kmt, halo_info, &
2221 37 : field_loc_center, field_type_scalar)
2222 : call ice_HaloUpdate (hm, halo_info, &
2223 37 : field_loc_center, field_type_scalar)
2224 37 : call ice_timer_stop(timer_bound)
2225 :
2226 : !-----------------------------------------------------------------
2227 : ! construct T-cell and U-cell masks
2228 : !-----------------------------------------------------------------
2229 :
2230 111936 : bm = c0
2231 37 : allocate(uvmCD(nx_block,ny_block,max_blocks))
2232 :
2233 20 : !$OMP PARALLEL DO PRIVATE(iblk,i,j,ilo,ihi,jlo,jhi,this_block)
2234 50 : do iblk = 1, nblocks
2235 33 : this_block = get_block(blocks_ice(iblk),iblk)
2236 33 : ilo = this_block%ilo
2237 33 : ihi = this_block%ihi
2238 33 : jlo = this_block%jlo
2239 33 : jhi = this_block%jhi
2240 :
2241 1210 : do j = jlo, jhi
2242 45541 : do i = ilo, ihi
2243 : uvm(i,j,iblk) = min (hm(i,j, iblk), hm(i+1,j, iblk), &
2244 44368 : hm(i,j+1,iblk), hm(i+1,j+1,iblk))
2245 44368 : npm(i,j,iblk) = min (hm(i,j, iblk), hm(i,j+1,iblk))
2246 44368 : epm(i,j,iblk) = min (hm(i,j, iblk), hm(i+1,j,iblk))
2247 44368 : bm(i,j,iblk) = my_task + iblk/100.0_dbl_kind
2248 : uvmCD(i,j,iblk) = (hm(i,j, iblk)+hm(i+1,j, iblk) &
2249 45508 : + hm(i,j+1,iblk)+hm(i+1,j+1,iblk))
2250 : enddo
2251 : enddo
2252 : enddo
2253 : !$OMP END PARALLEL DO
2254 :
2255 37 : call ice_timer_start(timer_bound)
2256 : call ice_HaloUpdate (uvm, halo_info, &
2257 37 : field_loc_NEcorner, field_type_scalar)
2258 : call ice_HaloUpdate (uvmCD, halo_info, &
2259 37 : field_loc_NEcorner, field_type_scalar)
2260 : call ice_HaloUpdate (npm, halo_info, &
2261 37 : field_loc_Nface, field_type_scalar)
2262 : call ice_HaloUpdate (epm, halo_info, &
2263 37 : field_loc_Eface, field_type_scalar)
2264 : call ice_HaloUpdate (bm, halo_info, &
2265 37 : field_loc_center, field_type_scalar)
2266 37 : call ice_timer_stop(timer_bound)
2267 :
2268 20 : !$OMP PARALLEL DO PRIVATE(iblk,i,j,ilo,ihi,jlo,jhi,this_block)
2269 50 : do iblk = 1, nblocks
2270 33 : this_block = get_block(blocks_ice(iblk),iblk)
2271 33 : ilo = this_block%ilo
2272 33 : ihi = this_block%ihi
2273 33 : jlo = this_block%jlo
2274 33 : jhi = this_block%jhi
2275 :
2276 : ! needs to cover halo (no halo update for logicals)
2277 50267 : tmask(:,:,iblk) = .false.
2278 50267 : umask(:,:,iblk) = .false.
2279 50267 : umaskCD(:,:,iblk) = .false.
2280 50267 : nmask(:,:,iblk) = .false.
2281 50267 : emask(:,:,iblk) = .false.
2282 1276 : do j = jlo-nghost, jhi+nghost
2283 50267 : do i = ilo-nghost, ihi+nghost
2284 49028 : if ( hm(i,j,iblk) > p5 ) tmask (i,j,iblk) = .true.
2285 49028 : if (uvm(i,j,iblk) > p5 ) umask (i,j,iblk) = .true.
2286 49028 : if (uvmCD(i,j,iblk) > c1p5) umaskCD(i,j,iblk) = .true.
2287 49028 : if (npm(i,j,iblk) > p5 ) nmask (i,j,iblk) = .true.
2288 50234 : if (epm(i,j,iblk) > p5 ) emask (i,j,iblk) = .true.
2289 : enddo
2290 : enddo
2291 :
2292 : enddo
2293 : !$OMP END PARALLEL DO
2294 :
2295 20 : !$OMP PARALLEL DO PRIVATE(iblk,i,j,ilo,ihi,jlo,jhi,this_block)
2296 50 : do iblk = 1, nblocks
2297 33 : this_block = get_block(blocks_ice(iblk),iblk)
2298 33 : ilo = this_block%ilo
2299 33 : ihi = this_block%ihi
2300 33 : jlo = this_block%jlo
2301 33 : jhi = this_block%jhi
2302 :
2303 : !-----------------------------------------------------------------
2304 : ! create hemisphere masks
2305 : !-----------------------------------------------------------------
2306 :
2307 50267 : lmask_n(:,:,iblk) = .false.
2308 50267 : lmask_s(:,:,iblk) = .false.
2309 :
2310 50267 : tarean(:,:,iblk) = c0
2311 50267 : tareas(:,:,iblk) = c0
2312 :
2313 1210 : do j = jlo,jhi
2314 45541 : do i = ilo,ihi
2315 :
2316 44368 : if (ULAT(i,j,iblk) >= -puny) then
2317 38968 : lmask_n(i,j,iblk) = .true. ! N. Hem.
2318 : else
2319 5400 : lmask_s(i,j,iblk) = .true. ! S. Hem.
2320 : endif
2321 :
2322 : ! N hemisphere area mask (m^2)
2323 44368 : if (lmask_n(i,j,iblk)) tarean(i,j,iblk) = tarea(i,j,iblk) &
2324 38968 : * hm(i,j,iblk)
2325 :
2326 : ! S hemisphere area mask (m^2)
2327 44368 : if (lmask_s(i,j,iblk)) tareas(i,j,iblk) = tarea(i,j,iblk) &
2328 6540 : * hm(i,j,iblk)
2329 :
2330 : enddo
2331 : enddo
2332 :
2333 : enddo ! iblk
2334 : !$OMP END PARALLEL DO
2335 :
2336 37 : deallocate(uvmCD)
2337 :
2338 74 : end subroutine makemask
2339 :
2340 : !=======================================================================
2341 :
2342 : ! Initializes latitude and longitude on T grid
2343 : !
2344 : ! author: Elizabeth C. Hunke, LANL; code originally based on POP grid
2345 : ! generation routine
2346 :
2347 37 : subroutine Tlatlon
2348 :
2349 : use ice_constants, only: c0, c1, c1p5, c2, c4, p5, &
2350 : field_loc_center, field_loc_Nface, field_loc_Eface, & ! LCOV_EXCL_LINE
2351 : field_type_scalar
2352 :
2353 : integer (kind=int_kind) :: &
2354 : i, j, iblk , & ! horizontal indices ! LCOV_EXCL_LINE
2355 : ilo,ihi,jlo,jhi ! beginning and end of physical domain
2356 :
2357 : real (kind=dbl_kind) :: &
2358 : z1,x1,y1,z2,x2,y2,z3,x3,y3,z4,x4,y4,tx,ty,tz,da, & ! LCOV_EXCL_LINE
2359 8 : rad_to_deg
2360 :
2361 : type (block) :: &
2362 : this_block ! block information for current block
2363 :
2364 : character(len=*), parameter :: subname = '(Tlatlon)'
2365 :
2366 37 : call icepack_query_parameters(rad_to_deg_out=rad_to_deg)
2367 37 : call icepack_warnings_flush(nu_diag)
2368 37 : if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
2369 0 : file=__FILE__, line=__LINE__)
2370 :
2371 111936 : TLAT(:,:,:) = c0
2372 111936 : TLON(:,:,:) = c0
2373 111936 : NLAT(:,:,:) = c0
2374 111936 : NLON(:,:,:) = c0
2375 111936 : ELAT(:,:,:) = c0
2376 111936 : ELON(:,:,:) = c0
2377 :
2378 : !$OMP PARALLEL DO PRIVATE(iblk,i,j,ilo,ihi,jlo,jhi,this_block, &
2379 : !$OMP x1,y1,z1,x2,y2,z2,x3,y3,z3,x4,y4,z4, & ! LCOV_EXCL_LINE
2380 20 : !$OMP tx,ty,tz,da)
2381 50 : do iblk = 1, nblocks
2382 33 : this_block = get_block(blocks_ice(iblk),iblk)
2383 33 : ilo = this_block%ilo
2384 33 : ihi = this_block%ihi
2385 33 : jlo = this_block%jlo
2386 33 : jhi = this_block%jhi
2387 :
2388 1210 : do j = jlo, jhi
2389 45541 : do i = ilo, ihi
2390 :
2391 44368 : z1 = cos(ULAT(i-1,j-1,iblk))
2392 44368 : x1 = cos(ULON(i-1,j-1,iblk))*z1
2393 44368 : y1 = sin(ULON(i-1,j-1,iblk))*z1
2394 44368 : z1 = sin(ULAT(i-1,j-1,iblk))
2395 :
2396 44368 : z2 = cos(ULAT(i,j-1,iblk))
2397 44368 : x2 = cos(ULON(i,j-1,iblk))*z2
2398 44368 : y2 = sin(ULON(i,j-1,iblk))*z2
2399 44368 : z2 = sin(ULAT(i,j-1,iblk))
2400 :
2401 44368 : z3 = cos(ULAT(i-1,j,iblk))
2402 44368 : x3 = cos(ULON(i-1,j,iblk))*z3
2403 44368 : y3 = sin(ULON(i-1,j,iblk))*z3
2404 44368 : z3 = sin(ULAT(i-1,j,iblk))
2405 :
2406 44368 : z4 = cos(ULAT(i,j,iblk))
2407 44368 : x4 = cos(ULON(i,j,iblk))*z4
2408 44368 : y4 = sin(ULON(i,j,iblk))*z4
2409 44368 : z4 = sin(ULAT(i,j,iblk))
2410 :
2411 : ! ---------
2412 : ! TLON/TLAT 4 pt computation (pts 1, 2, 3, 4)
2413 : ! ---------
2414 :
2415 44368 : tx = (x1+x2+x3+x4)/c4
2416 44368 : ty = (y1+y2+y3+y4)/c4
2417 44368 : tz = (z1+z2+z3+z4)/c4
2418 44368 : da = sqrt(tx**2+ty**2+tz**2)
2419 :
2420 44368 : tz = tz/da
2421 :
2422 : ! TLON in radians East
2423 44368 : TLON(i,j,iblk) = c0
2424 44368 : if (tx /= c0 .or. ty /= c0) TLON(i,j,iblk) = atan2(ty,tx)
2425 :
2426 : ! TLAT in radians North
2427 45508 : TLAT(i,j,iblk) = asin(tz)
2428 :
2429 : ! these two loops should be merged to save cos/sin calculations,
2430 : ! but atan2 is not bit-for-bit. This suggests the result for atan2 depends on
2431 : ! the prior atan2 call ??? not sure what's going on.
2432 : #if (1 == 1)
2433 : enddo ! i
2434 : enddo ! j
2435 : enddo ! iblk
2436 : !$OMP END PARALLEL DO
2437 :
2438 : !$OMP PARALLEL DO PRIVATE(iblk,i,j,ilo,ihi,jlo,jhi,this_block, &
2439 : !$OMP x1,y1,z1,x2,y2,z2,x3,y3,z3,x4,y4,z4, & ! LCOV_EXCL_LINE
2440 20 : !$OMP tx,ty,tz,da)
2441 50 : do iblk = 1, nblocks
2442 33 : this_block = get_block(blocks_ice(iblk),iblk)
2443 33 : ilo = this_block%ilo
2444 33 : ihi = this_block%ihi
2445 33 : jlo = this_block%jlo
2446 33 : jhi = this_block%jhi
2447 :
2448 1210 : do j = jlo, jhi
2449 45541 : do i = ilo, ihi
2450 :
2451 44368 : z1 = cos(ULAT(i-1,j-1,iblk))
2452 44368 : x1 = cos(ULON(i-1,j-1,iblk))*z1
2453 44368 : y1 = sin(ULON(i-1,j-1,iblk))*z1
2454 44368 : z1 = sin(ULAT(i-1,j-1,iblk))
2455 :
2456 44368 : z2 = cos(ULAT(i,j-1,iblk))
2457 44368 : x2 = cos(ULON(i,j-1,iblk))*z2
2458 44368 : y2 = sin(ULON(i,j-1,iblk))*z2
2459 44368 : z2 = sin(ULAT(i,j-1,iblk))
2460 :
2461 44368 : z3 = cos(ULAT(i-1,j,iblk))
2462 44368 : x3 = cos(ULON(i-1,j,iblk))*z3
2463 44368 : y3 = sin(ULON(i-1,j,iblk))*z3
2464 44368 : z3 = sin(ULAT(i-1,j,iblk))
2465 :
2466 44368 : z4 = cos(ULAT(i,j,iblk))
2467 44368 : x4 = cos(ULON(i,j,iblk))*z4
2468 44368 : y4 = sin(ULON(i,j,iblk))*z4
2469 44368 : z4 = sin(ULAT(i,j,iblk))
2470 : #endif
2471 : ! ---------
2472 : ! NLON/NLAT 2 pt computation (pts 3, 4)
2473 : ! ---------
2474 :
2475 44368 : tx = (x3+x4)/c2
2476 44368 : ty = (y3+y4)/c2
2477 44368 : tz = (z3+z4)/c2
2478 44368 : da = sqrt(tx**2+ty**2+tz**2)
2479 :
2480 44368 : tz = tz/da
2481 :
2482 : ! NLON in radians East
2483 44368 : NLON(i,j,iblk) = c0
2484 44368 : if (tx /= c0 .or. ty /= c0) NLON(i,j,iblk) = atan2(ty,tx)
2485 :
2486 : ! NLAT in radians North
2487 44368 : NLAT(i,j,iblk) = asin(tz)
2488 :
2489 : ! ---------
2490 : ! ELON/ELAT 2 pt computation (pts 2, 4)
2491 : ! ---------
2492 :
2493 44368 : tx = (x2+x4)/c2
2494 44368 : ty = (y2+y4)/c2
2495 44368 : tz = (z2+z4)/c2
2496 44368 : da = sqrt(tx**2+ty**2+tz**2)
2497 :
2498 44368 : tz = tz/da
2499 :
2500 : ! ELON in radians East
2501 44368 : ELON(i,j,iblk) = c0
2502 44368 : if (tx /= c0 .or. ty /= c0) ELON(i,j,iblk) = atan2(ty,tx)
2503 :
2504 : ! ELAT in radians North
2505 45508 : ELAT(i,j,iblk) = asin(tz)
2506 :
2507 : enddo ! i
2508 : enddo ! j
2509 : enddo ! iblk
2510 : !$OMP END PARALLEL DO
2511 :
2512 37 : if (trim(grid_type) == 'regional') then
2513 : ! for W boundary extrapolate from interior
2514 0 : !$OMP PARALLEL DO PRIVATE(iblk,i,j,ilo,ihi,jlo,jhi,this_block)
2515 0 : do iblk = 1, nblocks
2516 0 : this_block = get_block(blocks_ice(iblk),iblk)
2517 0 : ilo = this_block%ilo
2518 0 : ihi = this_block%ihi
2519 0 : jlo = this_block%jlo
2520 0 : jhi = this_block%jhi
2521 :
2522 0 : i = ilo
2523 0 : if (this_block%i_glob(i) == 1) then
2524 0 : do j = jlo, jhi
2525 : TLON(i,j,iblk) = c2*TLON(i+1,j,iblk) - &
2526 0 : TLON(i+2,j,iblk)
2527 : TLAT(i,j,iblk) = c2*TLAT(i+1,j,iblk) - &
2528 0 : TLAT(i+2,j,iblk)
2529 : NLON(i,j,iblk) = c1p5*TLON(i+1,j,iblk) - &
2530 0 : p5*TLON(i+2,j,iblk)
2531 : NLAT(i,j,iblk) = c1p5*TLAT(i+1,j,iblk) - &
2532 0 : p5*TLAT(i+2,j,iblk)
2533 : enddo
2534 : endif
2535 : enddo
2536 : !$OMP END PARALLEL DO
2537 : endif ! regional
2538 :
2539 37 : call ice_timer_start(timer_bound)
2540 : call ice_HaloUpdate (TLON, halo_info, &
2541 : field_loc_center, field_type_scalar, & ! LCOV_EXCL_LINE
2542 37 : fillValue=c1)
2543 : call ice_HaloUpdate (TLAT, halo_info, &
2544 : field_loc_center, field_type_scalar, & ! LCOV_EXCL_LINE
2545 37 : fillValue=c1)
2546 : call ice_HaloUpdate (NLON, halo_info, &
2547 : field_loc_Nface, field_type_scalar, & ! LCOV_EXCL_LINE
2548 37 : fillValue=c1)
2549 : call ice_HaloUpdate (NLAT, halo_info, &
2550 : field_loc_Nface, field_type_scalar, & ! LCOV_EXCL_LINE
2551 37 : fillValue=c1)
2552 : call ice_HaloUpdate (ELON, halo_info, &
2553 : field_loc_Eface, field_type_scalar, & ! LCOV_EXCL_LINE
2554 37 : fillValue=c1)
2555 : call ice_HaloUpdate (ELAT, halo_info, &
2556 : field_loc_Eface, field_type_scalar, & ! LCOV_EXCL_LINE
2557 37 : fillValue=c1)
2558 : call ice_HaloExtrapolate(TLON, distrb_info, &
2559 37 : ew_boundary_type, ns_boundary_type)
2560 : call ice_HaloExtrapolate(TLAT, distrb_info, &
2561 37 : ew_boundary_type, ns_boundary_type)
2562 : call ice_HaloExtrapolate(NLON, distrb_info, &
2563 37 : ew_boundary_type, ns_boundary_type)
2564 : call ice_HaloExtrapolate(NLAT, distrb_info, &
2565 37 : ew_boundary_type, ns_boundary_type)
2566 : call ice_HaloExtrapolate(ELON, distrb_info, &
2567 37 : ew_boundary_type, ns_boundary_type)
2568 : call ice_HaloExtrapolate(ELAT, distrb_info, &
2569 37 : ew_boundary_type, ns_boundary_type)
2570 37 : call ice_timer_stop(timer_bound)
2571 :
2572 37 : x1 = global_minval(TLON, distrb_info, tmask)
2573 37 : x2 = global_maxval(TLON, distrb_info, tmask)
2574 37 : x3 = global_minval(TLAT, distrb_info, tmask)
2575 37 : x4 = global_maxval(TLAT, distrb_info, tmask)
2576 :
2577 37 : y1 = global_minval(ULON, distrb_info, umask)
2578 37 : y2 = global_maxval(ULON, distrb_info, umask)
2579 37 : y3 = global_minval(ULAT, distrb_info, umask)
2580 37 : y4 = global_maxval(ULAT, distrb_info, umask)
2581 :
2582 37 : if (my_task==master_task) then
2583 7 : write(nu_diag,*) ' '
2584 : ! if (nx_block > 5+2*nghost .and. ny_block > 5+2*nghost) then
2585 7 : write(nu_diag,*) 'min/max ULON:', y1*rad_to_deg, y2*rad_to_deg
2586 7 : write(nu_diag,*) 'min/max ULAT:', y3*rad_to_deg, y4*rad_to_deg
2587 : ! endif
2588 7 : write(nu_diag,*) 'min/max TLON:', x1*rad_to_deg, x2*rad_to_deg
2589 7 : write(nu_diag,*) 'min/max TLAT:', x3*rad_to_deg, x4*rad_to_deg
2590 : endif ! my_task
2591 :
2592 37 : x1 = global_minval(NLON, distrb_info, nmask)
2593 37 : x2 = global_maxval(NLON, distrb_info, nmask)
2594 37 : x3 = global_minval(NLAT, distrb_info, nmask)
2595 37 : x4 = global_maxval(NLAT, distrb_info, nmask)
2596 :
2597 37 : y1 = global_minval(ELON, distrb_info, emask)
2598 37 : y2 = global_maxval(ELON, distrb_info, emask)
2599 37 : y3 = global_minval(ELAT, distrb_info, emask)
2600 37 : y4 = global_maxval(ELAT, distrb_info, emask)
2601 :
2602 37 : if (my_task==master_task) then
2603 7 : write(nu_diag,*) ' '
2604 : ! if (nx_block > 5+2*nghost .and. ny_block > 5+2*nghost) then
2605 7 : write(nu_diag,*) 'min/max NLON:', x1*rad_to_deg, x2*rad_to_deg
2606 7 : write(nu_diag,*) 'min/max NLAT:', x3*rad_to_deg, x4*rad_to_deg
2607 7 : write(nu_diag,*) 'min/max ELON:', y1*rad_to_deg, y2*rad_to_deg
2608 7 : write(nu_diag,*) 'min/max ELAT:', y3*rad_to_deg, y4*rad_to_deg
2609 : ! endif
2610 : endif ! my_task
2611 :
2612 37 : end subroutine Tlatlon
2613 :
2614 : !=======================================================================
2615 :
2616 : ! Shifts quantities from one grid to another
2617 : ! Constructs the shift based on the grid
2618 : ! NOTE: Input array includes ghost cells that must be updated before
2619 : ! calling this routine.
2620 : !
2621 : ! author: T. Craig
2622 :
2623 86760 : subroutine grid_average_X2Y_base(type,work1,grid1,work2,grid2)
2624 :
2625 : character(len=*) , intent(in) :: &
2626 : type, grid1, grid2
2627 :
2628 : real (kind=dbl_kind), intent(in) :: &
2629 : work1(:,:,:)
2630 :
2631 : real (kind=dbl_kind), intent(out) :: &
2632 : work2(:,:,:)
2633 :
2634 : ! local variables
2635 :
2636 : character(len=16) :: X2Y
2637 :
2638 : character(len=*), parameter :: subname = '(grid_average_X2Y_base)'
2639 :
2640 86760 : if (trim(grid1) == trim(grid2)) then
2641 59896128 : work2 = work1
2642 : else
2643 63672 : X2Y = trim(grid1)//'2'//trim(grid2)//trim(type)
2644 63672 : call grid_average_X2Y_1(X2Y,work1,work2)
2645 : endif
2646 :
2647 86760 : end subroutine grid_average_X2Y_base
2648 :
2649 : !=======================================================================
2650 :
2651 : ! Shifts quantities from one grid to another
2652 : ! NOTE: Input array includes ghost cells that must be updated before
2653 : ! calling this routine.
2654 : !
2655 : ! author: T. Craig
2656 :
2657 0 : subroutine grid_average_X2Y_userwghts(type,work1,grid1,wght1,mask1,work2,grid2)
2658 :
2659 : character(len=*) , intent(in) :: &
2660 : type, grid1, grid2
2661 :
2662 : real (kind=dbl_kind), intent(in) :: &
2663 : work1(:,:,:), & ! LCOV_EXCL_LINE
2664 : wght1(:,:,:), & ! LCOV_EXCL_LINE
2665 : mask1(:,:,:)
2666 :
2667 : real (kind=dbl_kind), intent(out) :: &
2668 : work2(:,:,:)
2669 :
2670 : ! local variables
2671 :
2672 : character(len=16) :: X2Y
2673 :
2674 : character(len=*), parameter :: subname = '(grid_average_X2Y_userwghts)'
2675 :
2676 0 : if (trim(grid1) == trim(grid2)) then
2677 0 : work2 = work1
2678 : else
2679 0 : X2Y = trim(grid1)//'2'//trim(grid2)//trim(type)
2680 0 : call grid_average_X2Y_1f(X2Y,work1,wght1,mask1,work2)
2681 : endif
2682 :
2683 0 : end subroutine grid_average_X2Y_userwghts
2684 :
2685 : !=======================================================================
2686 :
2687 : ! Shifts quantities from one grid to another
2688 : ! NOTE: Input array includes ghost cells that must be updated before
2689 : ! calling this routine.
2690 : !
2691 : ! author: T. Craig
2692 :
2693 0 : subroutine grid_average_X2Y_NEversion(type,work1a,grid1a,work1b,grid1b,work2,grid2)
2694 :
2695 : character(len=*) , intent(in) :: &
2696 : type, grid1a, grid1b, grid2
2697 :
2698 : real (kind=dbl_kind), intent(in) :: &
2699 : work1a(:,:,:), work1b(:,:,:)
2700 :
2701 : real (kind=dbl_kind), intent(out) :: &
2702 : work2(:,:,:)
2703 :
2704 : ! local variables
2705 :
2706 : character(len=16) :: X2Y
2707 :
2708 : character(len=*), parameter :: subname = '(grid_average_X2Y_NEversion)'
2709 :
2710 0 : X2Y = trim(grid1a)//trim(grid1b)//'2'//trim(grid2)//trim(type)
2711 :
2712 0 : select case (trim(X2Y))
2713 :
2714 : ! state masked
2715 : case('NE2US')
2716 0 : call grid_average_X2Y_2('NE2US',work1a,narea,npm,work1b,earea,epm,work2)
2717 : case('EN2US')
2718 0 : call grid_average_X2Y_2('NE2US',work1b,narea,npm,work1a,earea,epm,work2)
2719 : case('NE2TS')
2720 0 : call grid_average_X2Y_2('NE2TS',work1a,narea,npm,work1b,earea,epm,work2)
2721 : case('EN2TS')
2722 0 : call grid_average_X2Y_2('NE2TS',work1b,narea,npm,work1a,earea,epm,work2)
2723 :
2724 : ! state unmasked
2725 : case('NE2UA')
2726 0 : call grid_average_X2Y_2('NE2UA',work1a,narea,npm,work1b,earea,epm,work2)
2727 : case('EN2UA')
2728 0 : call grid_average_X2Y_2('NE2UA',work1b,narea,npm,work1a,earea,epm,work2)
2729 : case('NE2TA')
2730 0 : call grid_average_X2Y_2('NE2TA',work1a,narea,npm,work1b,earea,epm,work2)
2731 : case('EN2TA')
2732 0 : call grid_average_X2Y_2('NE2TA',work1b,narea,npm,work1a,earea,epm,work2)
2733 :
2734 : case default
2735 0 : call abort_ice(subname//'ERROR: unknown X2Y '//trim(X2Y))
2736 : end select
2737 :
2738 0 : end subroutine grid_average_X2Y_NEversion
2739 :
2740 : !=======================================================================
2741 :
2742 : ! Shifts quantities from one grid to another
2743 : ! NOTE: Input array includes ghost cells that must be updated before
2744 : ! calling this routine.
2745 : !
2746 : ! author: T. Craig
2747 :
2748 63672 : subroutine grid_average_X2Y_1(X2Y,work1,work2)
2749 :
2750 : character(len=*) , intent(in) :: &
2751 : X2Y
2752 :
2753 : real (kind=dbl_kind), intent(in) :: &
2754 : work1(:,:,:)
2755 :
2756 : real (kind=dbl_kind), intent(out) :: &
2757 : work2(:,:,:)
2758 :
2759 : ! local variables
2760 :
2761 : character(len=*), parameter :: subname = '(grid_average_X2Y_1)'
2762 :
2763 127344 : select case (trim(X2Y))
2764 :
2765 : ! flux unmasked
2766 : case('T2UF')
2767 11568 : call grid_average_X2YF('NE',work1,tarea,work2,uarea)
2768 : case('T2EF')
2769 0 : call grid_average_X2YF('E' ,work1,tarea,work2,earea)
2770 : case('T2NF')
2771 0 : call grid_average_X2YF('N' ,work1,tarea,work2,narea)
2772 : case('U2TF')
2773 11568 : call grid_average_X2YF('SW',work1,uarea,work2,tarea)
2774 : case('U2EF')
2775 0 : call grid_average_X2YF('S' ,work1,uarea,work2,earea)
2776 : case('U2NF')
2777 0 : call grid_average_X2YF('W' ,work1,uarea,work2,narea)
2778 : case('E2TF')
2779 0 : call grid_average_X2YF('W' ,work1,earea,work2,tarea)
2780 : case('E2UF')
2781 0 : call grid_average_X2YF('N' ,work1,earea,work2,uarea)
2782 : case('E2NF')
2783 0 : call grid_average_X2YF('NW',work1,earea,work2,narea)
2784 : case('N2TF')
2785 0 : call grid_average_X2YF('S' ,work1,narea,work2,tarea)
2786 : case('N2UF')
2787 0 : call grid_average_X2YF('E' ,work1,narea,work2,uarea)
2788 : case('N2EF')
2789 0 : call grid_average_X2YF('SE',work1,narea,work2,earea)
2790 :
2791 : ! state masked
2792 : case('T2US')
2793 28968 : call grid_average_X2YS('NE',work1,tarea,hm ,work2)
2794 : case('T2ES')
2795 0 : call grid_average_X2YS('E' ,work1,tarea,hm ,work2)
2796 : case('T2NS')
2797 0 : call grid_average_X2YS('N' ,work1,tarea,hm ,work2)
2798 : case('U2TS')
2799 0 : call grid_average_X2YS('SW',work1,uarea,uvm,work2)
2800 : case('U2ES')
2801 0 : call grid_average_X2YS('S' ,work1,uarea,uvm,work2)
2802 : case('U2NS')
2803 0 : call grid_average_X2YS('W' ,work1,uarea,uvm,work2)
2804 : case('E2TS')
2805 0 : call grid_average_X2YS('W' ,work1,earea,epm,work2)
2806 : case('E2US')
2807 0 : call grid_average_X2YS('N' ,work1,earea,epm,work2)
2808 : case('E2NS')
2809 0 : call grid_average_X2YS('NW',work1,earea,epm,work2)
2810 : case('N2TS')
2811 0 : call grid_average_X2YS('S' ,work1,narea,npm,work2)
2812 : case('N2US')
2813 0 : call grid_average_X2YS('E' ,work1,narea,npm,work2)
2814 : case('N2ES')
2815 0 : call grid_average_X2YS('SE',work1,narea,npm,work2)
2816 :
2817 : ! state unmasked
2818 : case('T2UA')
2819 0 : call grid_average_X2YA('NE',work1,tarea,work2)
2820 : case('T2EA')
2821 0 : call grid_average_X2YA('E' ,work1,tarea,work2)
2822 : case('T2NA')
2823 0 : call grid_average_X2YA('N' ,work1,tarea,work2)
2824 : case('U2TA')
2825 11568 : call grid_average_X2YA('SW',work1,uarea,work2)
2826 : case('U2EA')
2827 0 : call grid_average_X2YA('S' ,work1,uarea,work2)
2828 : case('U2NA')
2829 0 : call grid_average_X2YA('W' ,work1,uarea,work2)
2830 : case('E2TA')
2831 0 : call grid_average_X2YA('W' ,work1,earea,work2)
2832 : case('E2UA')
2833 0 : call grid_average_X2YA('N' ,work1,earea,work2)
2834 : case('E2NA')
2835 0 : call grid_average_X2YA('NW',work1,earea,work2)
2836 : case('N2TA')
2837 0 : call grid_average_X2YA('S' ,work1,narea,work2)
2838 : case('N2UA')
2839 0 : call grid_average_X2YA('E' ,work1,narea,work2)
2840 : case('N2EA')
2841 0 : call grid_average_X2YA('SE',work1,narea,work2)
2842 :
2843 : case default
2844 127344 : call abort_ice(subname//'ERROR: unknown X2Y '//trim(X2Y))
2845 : end select
2846 :
2847 63672 : end subroutine grid_average_X2Y_1
2848 :
2849 : !=======================================================================
2850 :
2851 : ! Shifts quantities from one grid to another
2852 : ! NOTE: Input array includes ghost cells that must be updated before
2853 : ! calling this routine.
2854 : !
2855 : ! author: T. Craig
2856 :
2857 0 : subroutine grid_average_X2Y_1f(X2Y,work1,wght1,mask1,work2)
2858 :
2859 : character(len=*) , intent(in) :: &
2860 : X2Y
2861 :
2862 : real (kind=dbl_kind), intent(in) :: &
2863 : work1(:,:,:), & ! LCOV_EXCL_LINE
2864 : wght1(:,:,:), & ! LCOV_EXCL_LINE
2865 : mask1(:,:,:)
2866 :
2867 : real (kind=dbl_kind), intent(out) :: &
2868 : work2(:,:,:)
2869 :
2870 : ! local variables
2871 :
2872 : character(len=*), parameter :: subname = '(grid_average_X2Y_1f)'
2873 :
2874 0 : select case (trim(X2Y))
2875 :
2876 : ! don't support these for now, requires extra destination wght
2877 : ! ! flux unmasked
2878 : ! case('T2UF')
2879 : ! call grid_average_X2YF('NE',work1,tarea,work2,uarea)
2880 : ! case('T2EF')
2881 : ! call grid_average_X2YF('E' ,work1,tarea,work2,earea)
2882 : ! case('T2NF')
2883 : ! call grid_average_X2YF('N' ,work1,tarea,work2,narea)
2884 : ! case('U2TF')
2885 : ! call grid_average_X2YF('SW',work1,uarea,work2,tarea)
2886 : ! case('U2EF')
2887 : ! call grid_average_X2YF('S' ,work1,uarea,work2,earea)
2888 : ! case('U2NF')
2889 : ! call grid_average_X2YF('W' ,work1,uarea,work2,narea)
2890 : ! case('E2TF')
2891 : ! call grid_average_X2YF('W' ,work1,earea,work2,tarea)
2892 : ! case('E2UF')
2893 : ! call grid_average_X2YF('N' ,work1,earea,work2,uarea)
2894 : ! case('E2NF')
2895 : ! call grid_average_X2YF('NW',work1,earea,work2,narea)
2896 : ! case('N2TF')
2897 : ! call grid_average_X2YF('S' ,work1,narea,work2,tarea)
2898 : ! case('N2UF')
2899 : ! call grid_average_X2YF('E' ,work1,narea,work2,uarea)
2900 : ! case('N2EF')
2901 : ! call grid_average_X2YF('SE',work1,narea,work2,earea)
2902 :
2903 : ! state masked
2904 : case('T2US')
2905 0 : call grid_average_X2YS('NE',work1,wght1,mask1,work2)
2906 : case('T2ES')
2907 0 : call grid_average_X2YS('E' ,work1,wght1,mask1,work2)
2908 : case('T2NS')
2909 0 : call grid_average_X2YS('N' ,work1,wght1,mask1,work2)
2910 : case('U2TS')
2911 0 : call grid_average_X2YS('SW',work1,wght1,mask1,work2)
2912 : case('U2ES')
2913 0 : call grid_average_X2YS('S' ,work1,wght1,mask1,work2)
2914 : case('U2NS')
2915 0 : call grid_average_X2YS('W' ,work1,wght1,mask1,work2)
2916 : case('E2TS')
2917 0 : call grid_average_X2YS('W' ,work1,wght1,mask1,work2)
2918 : case('E2US')
2919 0 : call grid_average_X2YS('N' ,work1,wght1,mask1,work2)
2920 : case('E2NS')
2921 0 : call grid_average_X2YS('NW',work1,wght1,mask1,work2)
2922 : case('N2TS')
2923 0 : call grid_average_X2YS('S' ,work1,wght1,mask1,work2)
2924 : case('N2US')
2925 0 : call grid_average_X2YS('E' ,work1,wght1,mask1,work2)
2926 : case('N2ES')
2927 0 : call grid_average_X2YS('SE',work1,wght1,mask1,work2)
2928 :
2929 : ! state unmasked
2930 : case('T2UA')
2931 0 : call grid_average_X2YA('NE',work1,wght1,work2)
2932 : case('T2EA')
2933 0 : call grid_average_X2YA('E' ,work1,wght1,work2)
2934 : case('T2NA')
2935 0 : call grid_average_X2YA('N' ,work1,wght1,work2)
2936 : case('U2TA')
2937 0 : call grid_average_X2YA('SW',work1,wght1,work2)
2938 : case('U2EA')
2939 0 : call grid_average_X2YA('S' ,work1,wght1,work2)
2940 : case('U2NA')
2941 0 : call grid_average_X2YA('W' ,work1,wght1,work2)
2942 : case('E2TA')
2943 0 : call grid_average_X2YA('W' ,work1,wght1,work2)
2944 : case('E2UA')
2945 0 : call grid_average_X2YA('N' ,work1,wght1,work2)
2946 : case('E2NA')
2947 0 : call grid_average_X2YA('NW',work1,wght1,work2)
2948 : case('N2TA')
2949 0 : call grid_average_X2YA('S' ,work1,wght1,work2)
2950 : case('N2UA')
2951 0 : call grid_average_X2YA('E' ,work1,wght1,work2)
2952 : case('N2EA')
2953 0 : call grid_average_X2YA('SE',work1,wght1,work2)
2954 :
2955 : case default
2956 0 : call abort_ice(subname//'ERROR: unknown X2Y '//trim(X2Y))
2957 : end select
2958 :
2959 0 : end subroutine grid_average_X2Y_1f
2960 :
2961 : !=======================================================================
2962 : ! Shifts quantities from one grid to another
2963 : ! State masked version, simple area weighted averager
2964 : ! NOTE: Input array includes ghost cells that must be updated before
2965 : ! calling this routine.
2966 : !
2967 : ! author: T. Craig
2968 :
2969 28968 : subroutine grid_average_X2YS(dir,work1,wght1,mask1,work2)
2970 :
2971 : use ice_constants, only: c0
2972 :
2973 : character(len=*) , intent(in) :: &
2974 : dir
2975 :
2976 : real (kind=dbl_kind), intent(in) :: &
2977 : work1(:,:,:), & ! LCOV_EXCL_LINE
2978 : wght1(:,:,:), & ! LCOV_EXCL_LINE
2979 : mask1(:,:,:)
2980 :
2981 : real (kind=dbl_kind), intent(out) :: &
2982 : work2(:,:,:)
2983 :
2984 : ! local variables
2985 :
2986 : integer (kind=int_kind) :: &
2987 : i, j, iblk, & ! LCOV_EXCL_LINE
2988 : ilo,ihi,jlo,jhi ! beginning and end of physical domain
2989 :
2990 : real (kind=dbl_kind) :: &
2991 10080 : wtmp
2992 :
2993 : type (block) :: &
2994 : this_block ! block information for current block
2995 :
2996 : character(len=*), parameter :: subname = '(grid_average_X2YS)'
2997 :
2998 86101728 : work2(:,:,:) = c0
2999 :
3000 57936 : select case (trim(dir))
3001 :
3002 : case('NE')
3003 20160 : !$OMP PARALLEL DO PRIVATE(iblk,i,j,ilo,ihi,jlo,jhi,this_block,wtmp)
3004 35064 : do iblk = 1, nblocks
3005 17448 : this_block = get_block(blocks_ice(iblk),iblk)
3006 17448 : ilo = this_block%ilo
3007 17448 : ihi = this_block%ihi
3008 17448 : jlo = this_block%jlo
3009 17448 : jhi = this_block%jhi
3010 618864 : do j = jlo, jhi
3011 20233416 : do i = ilo, ihi
3012 : wtmp = (mask1(i ,j ,iblk)*wght1(i ,j ,iblk) &
3013 : + mask1(i+1,j ,iblk)*wght1(i+1,j ,iblk) & ! LCOV_EXCL_LINE
3014 : + mask1(i ,j+1,iblk)*wght1(i ,j+1,iblk) & ! LCOV_EXCL_LINE
3015 19643520 : + mask1(i+1,j+1,iblk)*wght1(i+1,j+1,iblk))
3016 19643520 : if (wtmp /= c0) &
3017 : work2(i,j,iblk) = (mask1(i ,j ,iblk)*work1(i ,j ,iblk)*wght1(i ,j ,iblk) & ! LCOV_EXCL_LINE
3018 : + mask1(i+1,j ,iblk)*work1(i+1,j ,iblk)*wght1(i+1,j ,iblk) & ! LCOV_EXCL_LINE
3019 : + mask1(i ,j+1,iblk)*work1(i ,j+1,iblk)*wght1(i ,j+1,iblk) & ! LCOV_EXCL_LINE
3020 : + mask1(i+1,j+1,iblk)*work1(i+1,j+1,iblk)*wght1(i+1,j+1,iblk)) & ! LCOV_EXCL_LINE
3021 19403592 : / wtmp
3022 : enddo
3023 : enddo
3024 : enddo
3025 : !$OMP END PARALLEL DO
3026 :
3027 : case('SW')
3028 0 : !$OMP PARALLEL DO PRIVATE(iblk,i,j,ilo,ihi,jlo,jhi,this_block,wtmp)
3029 0 : do iblk = 1, nblocks
3030 0 : this_block = get_block(blocks_ice(iblk),iblk)
3031 0 : ilo = this_block%ilo
3032 0 : ihi = this_block%ihi
3033 0 : jlo = this_block%jlo
3034 0 : jhi = this_block%jhi
3035 0 : do j = jlo, jhi
3036 0 : do i = ilo, ihi
3037 : wtmp = (mask1(i ,j ,iblk)*wght1(i ,j ,iblk) &
3038 : + mask1(i-1,j ,iblk)*wght1(i-1,j ,iblk) & ! LCOV_EXCL_LINE
3039 : + mask1(i ,j-1,iblk)*wght1(i ,j-1,iblk) & ! LCOV_EXCL_LINE
3040 0 : + mask1(i-1,j-1,iblk)*wght1(i-1,j-1,iblk))
3041 0 : if (wtmp /= c0) &
3042 : work2(i,j,iblk) = (mask1(i ,j ,iblk)*work1(i ,j ,iblk)*wght1(i ,j ,iblk) & ! LCOV_EXCL_LINE
3043 : + mask1(i-1,j ,iblk)*work1(i-1,j ,iblk)*wght1(i-1,j ,iblk) & ! LCOV_EXCL_LINE
3044 : + mask1(i ,j-1,iblk)*work1(i ,j-1,iblk)*wght1(i ,j-1,iblk) & ! LCOV_EXCL_LINE
3045 : + mask1(i-1,j-1,iblk)*work1(i-1,j-1,iblk)*wght1(i-1,j-1,iblk)) & ! LCOV_EXCL_LINE
3046 0 : / wtmp
3047 : enddo
3048 : enddo
3049 : enddo
3050 : !$OMP END PARALLEL DO
3051 :
3052 : case('NW')
3053 0 : !$OMP PARALLEL DO PRIVATE(iblk,i,j,ilo,ihi,jlo,jhi,this_block,wtmp)
3054 0 : do iblk = 1, nblocks
3055 0 : this_block = get_block(blocks_ice(iblk),iblk)
3056 0 : ilo = this_block%ilo
3057 0 : ihi = this_block%ihi
3058 0 : jlo = this_block%jlo
3059 0 : jhi = this_block%jhi
3060 0 : do j = jlo, jhi
3061 0 : do i = ilo, ihi
3062 : wtmp = (mask1(i-1,j ,iblk)*wght1(i-1,j ,iblk) &
3063 : + mask1(i ,j ,iblk)*wght1(i ,j ,iblk) & ! LCOV_EXCL_LINE
3064 : + mask1(i-1,j+1,iblk)*wght1(i-1,j+1,iblk) & ! LCOV_EXCL_LINE
3065 0 : + mask1(i ,j+1,iblk)*wght1(i ,j+1,iblk))
3066 0 : if (wtmp /= c0) &
3067 : work2(i,j,iblk) = (mask1(i-1,j ,iblk)*work1(i-1,j ,iblk)*wght1(i-1,j ,iblk) & ! LCOV_EXCL_LINE
3068 : + mask1(i ,j ,iblk)*work1(i ,j ,iblk)*wght1(i ,j ,iblk) & ! LCOV_EXCL_LINE
3069 : + mask1(i-1,j+1,iblk)*work1(i-1,j+1,iblk)*wght1(i-1,j+1,iblk) & ! LCOV_EXCL_LINE
3070 : + mask1(i ,j+1,iblk)*work1(i ,j+1,iblk)*wght1(i ,j+1,iblk)) & ! LCOV_EXCL_LINE
3071 0 : / wtmp
3072 : enddo
3073 : enddo
3074 : enddo
3075 : !$OMP END PARALLEL DO
3076 :
3077 : case('SE')
3078 0 : !$OMP PARALLEL DO PRIVATE(iblk,i,j,ilo,ihi,jlo,jhi,this_block,wtmp)
3079 0 : do iblk = 1, nblocks
3080 0 : this_block = get_block(blocks_ice(iblk),iblk)
3081 0 : ilo = this_block%ilo
3082 0 : ihi = this_block%ihi
3083 0 : jlo = this_block%jlo
3084 0 : jhi = this_block%jhi
3085 0 : do j = jlo, jhi
3086 0 : do i = ilo, ihi
3087 : wtmp = (mask1(i ,j-1,iblk)*wght1(i ,j-1,iblk) &
3088 : + mask1(i+1,j-1,iblk)*wght1(i+1,j-1,iblk) & ! LCOV_EXCL_LINE
3089 : + mask1(i ,j ,iblk)*wght1(i ,j ,iblk) & ! LCOV_EXCL_LINE
3090 0 : + mask1(i+1,j ,iblk)*wght1(i+1,j ,iblk))
3091 0 : if (wtmp /= c0) &
3092 : work2(i,j,iblk) = (mask1(i ,j-1,iblk)*work1(i ,j-1,iblk)*wght1(i ,j-1,iblk) & ! LCOV_EXCL_LINE
3093 : + mask1(i+1,j-1,iblk)*work1(i+1,j-1,iblk)*wght1(i+1,j-1,iblk) & ! LCOV_EXCL_LINE
3094 : + mask1(i ,j ,iblk)*work1(i ,j ,iblk)*wght1(i ,j ,iblk) & ! LCOV_EXCL_LINE
3095 : + mask1(i+1,j ,iblk)*work1(i+1,j ,iblk)*wght1(i+1,j ,iblk)) & ! LCOV_EXCL_LINE
3096 0 : / wtmp
3097 : enddo
3098 : enddo
3099 : enddo
3100 : !$OMP END PARALLEL DO
3101 :
3102 : case('E')
3103 0 : !$OMP PARALLEL DO PRIVATE(iblk,i,j,ilo,ihi,jlo,jhi,this_block,wtmp)
3104 0 : do iblk = 1, nblocks
3105 0 : this_block = get_block(blocks_ice(iblk),iblk)
3106 0 : ilo = this_block%ilo
3107 0 : ihi = this_block%ihi
3108 0 : jlo = this_block%jlo
3109 0 : jhi = this_block%jhi
3110 0 : do j = jlo, jhi
3111 0 : do i = ilo, ihi
3112 : wtmp = (mask1(i ,j,iblk)*wght1(i ,j,iblk) &
3113 0 : + mask1(i+1,j,iblk)*wght1(i+1,j,iblk))
3114 0 : if (wtmp /= c0) &
3115 : work2(i,j,iblk) = (mask1(i ,j,iblk)*work1(i ,j,iblk)*wght1(i ,j,iblk) & ! LCOV_EXCL_LINE
3116 : + mask1(i+1,j,iblk)*work1(i+1,j,iblk)*wght1(i+1,j,iblk)) & ! LCOV_EXCL_LINE
3117 0 : / wtmp
3118 : enddo
3119 : enddo
3120 : enddo
3121 : !$OMP END PARALLEL DO
3122 :
3123 : case('W')
3124 0 : !$OMP PARALLEL DO PRIVATE(iblk,i,j,ilo,ihi,jlo,jhi,this_block,wtmp)
3125 0 : do iblk = 1, nblocks
3126 0 : this_block = get_block(blocks_ice(iblk),iblk)
3127 0 : ilo = this_block%ilo
3128 0 : ihi = this_block%ihi
3129 0 : jlo = this_block%jlo
3130 0 : jhi = this_block%jhi
3131 0 : do j = jlo, jhi
3132 0 : do i = ilo, ihi
3133 : wtmp = (mask1(i-1,j,iblk)*wght1(i-1,j,iblk) &
3134 0 : + mask1(i ,j,iblk)*wght1(i ,j,iblk))
3135 0 : if (wtmp /= c0) &
3136 : work2(i,j,iblk) = (mask1(i-1,j,iblk)*work1(i-1,j,iblk)*wght1(i-1,j,iblk) & ! LCOV_EXCL_LINE
3137 : + mask1(i ,j,iblk)*work1(i ,j,iblk)*wght1(i ,j,iblk)) & ! LCOV_EXCL_LINE
3138 0 : / wtmp
3139 : enddo
3140 : enddo
3141 : enddo
3142 : !$OMP END PARALLEL DO
3143 :
3144 : case('N')
3145 0 : !$OMP PARALLEL DO PRIVATE(iblk,i,j,ilo,ihi,jlo,jhi,this_block,wtmp)
3146 0 : do iblk = 1, nblocks
3147 0 : this_block = get_block(blocks_ice(iblk),iblk)
3148 0 : ilo = this_block%ilo
3149 0 : ihi = this_block%ihi
3150 0 : jlo = this_block%jlo
3151 0 : jhi = this_block%jhi
3152 0 : do j = jlo, jhi
3153 0 : do i = ilo, ihi
3154 : wtmp = (mask1(i,j ,iblk)*wght1(i,j ,iblk) &
3155 0 : + mask1(i,j+1,iblk)*wght1(i,j+1,iblk))
3156 0 : if (wtmp /= c0) &
3157 : work2(i,j,iblk) = (mask1(i,j ,iblk)*work1(i,j ,iblk)*wght1(i,j ,iblk) & ! LCOV_EXCL_LINE
3158 : + mask1(i,j+1,iblk)*work1(i,j+1,iblk)*wght1(i,j+1,iblk)) & ! LCOV_EXCL_LINE
3159 0 : / wtmp
3160 : enddo
3161 : enddo
3162 : enddo
3163 : !$OMP END PARALLEL DO
3164 :
3165 : case('S')
3166 0 : !$OMP PARALLEL DO PRIVATE(iblk,i,j,ilo,ihi,jlo,jhi,this_block,wtmp)
3167 0 : do iblk = 1, nblocks
3168 0 : this_block = get_block(blocks_ice(iblk),iblk)
3169 0 : ilo = this_block%ilo
3170 0 : ihi = this_block%ihi
3171 0 : jlo = this_block%jlo
3172 0 : jhi = this_block%jhi
3173 0 : do j = jlo, jhi
3174 0 : do i = ilo, ihi
3175 : wtmp = (mask1(i,j-1,iblk)*wght1(i,j-1,iblk) &
3176 0 : + mask1(i,j ,iblk)*wght1(i,j ,iblk))
3177 0 : if (wtmp /= c0) &
3178 : work2(i,j,iblk) = (mask1(i,j-1,iblk)*work1(i,j-1,iblk)*wght1(i,j-1,iblk) & ! LCOV_EXCL_LINE
3179 : + mask1(i,j ,iblk)*work1(i,j ,iblk)*wght1(i,j ,iblk)) & ! LCOV_EXCL_LINE
3180 0 : / wtmp
3181 : enddo
3182 : enddo
3183 : enddo
3184 : !$OMP END PARALLEL DO
3185 :
3186 : case default
3187 57936 : call abort_ice(subname//'ERROR: unknown option '//trim(dir))
3188 : end select
3189 :
3190 28968 : end subroutine grid_average_X2YS
3191 :
3192 : !=======================================================================
3193 : ! Shifts quantities from one grid to another
3194 : ! State unmasked version, simple weighted averager
3195 : ! NOTE: Input array includes ghost cells that must be updated before
3196 : ! calling this routine.
3197 : !
3198 : ! author: T. Craig
3199 :
3200 11568 : subroutine grid_average_X2YA(dir,work1,wght1,work2)
3201 :
3202 : use ice_constants, only: c0
3203 :
3204 : character(len=*) , intent(in) :: &
3205 : dir
3206 :
3207 : real (kind=dbl_kind), intent(in) :: &
3208 : work1(:,:,:), & ! LCOV_EXCL_LINE
3209 : wght1(:,:,:)
3210 :
3211 : real (kind=dbl_kind), intent(out) :: &
3212 : work2(:,:,:)
3213 :
3214 : ! local variables
3215 :
3216 : integer (kind=int_kind) :: &
3217 : i, j, iblk, & ! LCOV_EXCL_LINE
3218 : ilo,ihi,jlo,jhi ! beginning and end of physical domain
3219 :
3220 : real (kind=dbl_kind) :: &
3221 2880 : wtmp
3222 :
3223 : type (block) :: &
3224 : this_block ! block information for current block
3225 :
3226 : character(len=*), parameter :: subname = '(grid_average_X2YA)'
3227 :
3228 32443968 : work2(:,:,:) = c0
3229 :
3230 23136 : select case (trim(dir))
3231 :
3232 : case('NE')
3233 5760 : !$OMP PARALLEL DO PRIVATE(iblk,i,j,ilo,ihi,jlo,jhi,this_block,wtmp)
3234 0 : do iblk = 1, nblocks
3235 0 : this_block = get_block(blocks_ice(iblk),iblk)
3236 0 : ilo = this_block%ilo
3237 0 : ihi = this_block%ihi
3238 0 : jlo = this_block%jlo
3239 0 : jhi = this_block%jhi
3240 0 : do j = jlo, jhi
3241 0 : do i = ilo, ihi
3242 : wtmp = (wght1(i ,j ,iblk) &
3243 : + wght1(i+1,j ,iblk) & ! LCOV_EXCL_LINE
3244 : + wght1(i ,j+1,iblk) & ! LCOV_EXCL_LINE
3245 0 : + wght1(i+1,j+1,iblk))
3246 0 : if (wtmp /= c0) &
3247 : work2(i,j,iblk) = (work1(i ,j ,iblk)*wght1(i ,j ,iblk) & ! LCOV_EXCL_LINE
3248 : + work1(i+1,j ,iblk)*wght1(i+1,j ,iblk) & ! LCOV_EXCL_LINE
3249 : + work1(i ,j+1,iblk)*wght1(i ,j+1,iblk) & ! LCOV_EXCL_LINE
3250 : + work1(i+1,j+1,iblk)*wght1(i+1,j+1,iblk)) & ! LCOV_EXCL_LINE
3251 0 : / wtmp
3252 : enddo
3253 : enddo
3254 : enddo
3255 : !$OMP END PARALLEL DO
3256 :
3257 : case('SW')
3258 5760 : !$OMP PARALLEL DO PRIVATE(iblk,i,j,ilo,ihi,jlo,jhi,this_block,wtmp)
3259 23184 : do iblk = 1, nblocks
3260 11568 : this_block = get_block(blocks_ice(iblk),iblk)
3261 11568 : ilo = this_block%ilo
3262 11568 : ihi = this_block%ihi
3263 11568 : jlo = this_block%jlo
3264 11568 : jhi = this_block%jhi
3265 397344 : do j = jlo, jhi
3266 12739056 : do i = ilo, ihi
3267 : wtmp = (wght1(i ,j ,iblk) &
3268 : + wght1(i-1,j ,iblk) & ! LCOV_EXCL_LINE
3269 : + wght1(i ,j-1,iblk) & ! LCOV_EXCL_LINE
3270 12353280 : + wght1(i-1,j-1,iblk))
3271 12353280 : if (wtmp /= c0) &
3272 : work2(i,j,iblk) = (work1(i ,j ,iblk)*wght1(i ,j ,iblk) & ! LCOV_EXCL_LINE
3273 : + work1(i-1,j ,iblk)*wght1(i-1,j ,iblk) & ! LCOV_EXCL_LINE
3274 : + work1(i ,j-1,iblk)*wght1(i ,j-1,iblk) & ! LCOV_EXCL_LINE
3275 : + work1(i-1,j-1,iblk)*wght1(i-1,j-1,iblk)) & ! LCOV_EXCL_LINE
3276 12727488 : / wtmp
3277 : enddo
3278 : enddo
3279 : enddo
3280 : !$OMP END PARALLEL DO
3281 :
3282 : case('NW')
3283 0 : !$OMP PARALLEL DO PRIVATE(iblk,i,j,ilo,ihi,jlo,jhi,this_block,wtmp)
3284 0 : do iblk = 1, nblocks
3285 0 : this_block = get_block(blocks_ice(iblk),iblk)
3286 0 : ilo = this_block%ilo
3287 0 : ihi = this_block%ihi
3288 0 : jlo = this_block%jlo
3289 0 : jhi = this_block%jhi
3290 0 : do j = jlo, jhi
3291 0 : do i = ilo, ihi
3292 : wtmp = (wght1(i-1,j ,iblk) &
3293 : + wght1(i ,j ,iblk) & ! LCOV_EXCL_LINE
3294 : + wght1(i-1,j+1,iblk) & ! LCOV_EXCL_LINE
3295 0 : + wght1(i ,j+1,iblk))
3296 0 : if (wtmp /= c0) &
3297 : work2(i,j,iblk) = (work1(i-1,j ,iblk)*wght1(i-1,j ,iblk) & ! LCOV_EXCL_LINE
3298 : + work1(i ,j ,iblk)*wght1(i ,j ,iblk) & ! LCOV_EXCL_LINE
3299 : + work1(i-1,j+1,iblk)*wght1(i-1,j+1,iblk) & ! LCOV_EXCL_LINE
3300 : + work1(i ,j+1,iblk)*wght1(i ,j+1,iblk)) & ! LCOV_EXCL_LINE
3301 0 : / wtmp
3302 : enddo
3303 : enddo
3304 : enddo
3305 : !$OMP END PARALLEL DO
3306 :
3307 : case('SE')
3308 0 : !$OMP PARALLEL DO PRIVATE(iblk,i,j,ilo,ihi,jlo,jhi,this_block,wtmp)
3309 0 : do iblk = 1, nblocks
3310 0 : this_block = get_block(blocks_ice(iblk),iblk)
3311 0 : ilo = this_block%ilo
3312 0 : ihi = this_block%ihi
3313 0 : jlo = this_block%jlo
3314 0 : jhi = this_block%jhi
3315 0 : do j = jlo, jhi
3316 0 : do i = ilo, ihi
3317 : wtmp = (wght1(i ,j-1,iblk) &
3318 : + wght1(i+1,j-1,iblk) & ! LCOV_EXCL_LINE
3319 : + wght1(i ,j ,iblk) & ! LCOV_EXCL_LINE
3320 0 : + wght1(i+1,j ,iblk))
3321 0 : if (wtmp /= c0) &
3322 : work2(i,j,iblk) = (work1(i ,j-1,iblk)*wght1(i ,j-1,iblk) & ! LCOV_EXCL_LINE
3323 : + work1(i+1,j-1,iblk)*wght1(i+1,j-1,iblk) & ! LCOV_EXCL_LINE
3324 : + work1(i ,j ,iblk)*wght1(i ,j ,iblk) & ! LCOV_EXCL_LINE
3325 : + work1(i+1,j ,iblk)*wght1(i+1,j ,iblk)) & ! LCOV_EXCL_LINE
3326 0 : / wtmp
3327 : enddo
3328 : enddo
3329 : enddo
3330 : !$OMP END PARALLEL DO
3331 :
3332 : case('E')
3333 0 : !$OMP PARALLEL DO PRIVATE(iblk,i,j,ilo,ihi,jlo,jhi,this_block,wtmp)
3334 0 : do iblk = 1, nblocks
3335 0 : this_block = get_block(blocks_ice(iblk),iblk)
3336 0 : ilo = this_block%ilo
3337 0 : ihi = this_block%ihi
3338 0 : jlo = this_block%jlo
3339 0 : jhi = this_block%jhi
3340 0 : do j = jlo, jhi
3341 0 : do i = ilo, ihi
3342 : wtmp = (wght1(i ,j,iblk) &
3343 0 : + wght1(i+1,j,iblk))
3344 0 : if (wtmp /= c0) &
3345 : work2(i,j,iblk) = (work1(i ,j,iblk)*wght1(i ,j,iblk) & ! LCOV_EXCL_LINE
3346 : + work1(i+1,j,iblk)*wght1(i+1,j,iblk)) & ! LCOV_EXCL_LINE
3347 0 : / wtmp
3348 : enddo
3349 : enddo
3350 : enddo
3351 : !$OMP END PARALLEL DO
3352 :
3353 : case('W')
3354 0 : !$OMP PARALLEL DO PRIVATE(iblk,i,j,ilo,ihi,jlo,jhi,this_block,wtmp)
3355 0 : do iblk = 1, nblocks
3356 0 : this_block = get_block(blocks_ice(iblk),iblk)
3357 0 : ilo = this_block%ilo
3358 0 : ihi = this_block%ihi
3359 0 : jlo = this_block%jlo
3360 0 : jhi = this_block%jhi
3361 0 : do j = jlo, jhi
3362 0 : do i = ilo, ihi
3363 : wtmp = (wght1(i-1,j,iblk) &
3364 0 : + wght1(i ,j,iblk))
3365 0 : if (wtmp /= c0) &
3366 : work2(i,j,iblk) = (work1(i-1,j,iblk)*wght1(i-1,j,iblk) & ! LCOV_EXCL_LINE
3367 : + work1(i ,j,iblk)*wght1(i ,j,iblk)) & ! LCOV_EXCL_LINE
3368 0 : / wtmp
3369 : enddo
3370 : enddo
3371 : enddo
3372 : !$OMP END PARALLEL DO
3373 :
3374 : case('N')
3375 0 : !$OMP PARALLEL DO PRIVATE(iblk,i,j,ilo,ihi,jlo,jhi,this_block,wtmp)
3376 0 : do iblk = 1, nblocks
3377 0 : this_block = get_block(blocks_ice(iblk),iblk)
3378 0 : ilo = this_block%ilo
3379 0 : ihi = this_block%ihi
3380 0 : jlo = this_block%jlo
3381 0 : jhi = this_block%jhi
3382 0 : do j = jlo, jhi
3383 0 : do i = ilo, ihi
3384 : wtmp = (wght1(i,j ,iblk) &
3385 0 : + wght1(i,j+1,iblk))
3386 0 : if (wtmp /= c0) &
3387 : work2(i,j,iblk) = (work1(i,j ,iblk)*wght1(i,j ,iblk) & ! LCOV_EXCL_LINE
3388 : + work1(i,j+1,iblk)*wght1(i,j+1,iblk)) & ! LCOV_EXCL_LINE
3389 0 : / wtmp
3390 : enddo
3391 : enddo
3392 : enddo
3393 : !$OMP END PARALLEL DO
3394 :
3395 : case('S')
3396 0 : !$OMP PARALLEL DO PRIVATE(iblk,i,j,ilo,ihi,jlo,jhi,this_block,wtmp)
3397 0 : do iblk = 1, nblocks
3398 0 : this_block = get_block(blocks_ice(iblk),iblk)
3399 0 : ilo = this_block%ilo
3400 0 : ihi = this_block%ihi
3401 0 : jlo = this_block%jlo
3402 0 : jhi = this_block%jhi
3403 0 : do j = jlo, jhi
3404 0 : do i = ilo, ihi
3405 : wtmp = (wght1(i,j-1,iblk) &
3406 0 : + wght1(i,j ,iblk))
3407 0 : if (wtmp /= c0) &
3408 : work2(i,j,iblk) = (work1(i,j-1,iblk)*wght1(i,j-1,iblk) & ! LCOV_EXCL_LINE
3409 : + work1(i,j ,iblk)*wght1(i,j ,iblk)) & ! LCOV_EXCL_LINE
3410 0 : / wtmp
3411 : enddo
3412 : enddo
3413 : enddo
3414 : !$OMP END PARALLEL DO
3415 :
3416 : case default
3417 23136 : call abort_ice(subname//'ERROR: unknown option '//trim(dir))
3418 : end select
3419 :
3420 11568 : end subroutine grid_average_X2YA
3421 :
3422 : !=======================================================================
3423 : ! Shifts quantities from one grid to another
3424 : ! Flux masked, original implementation based on earlier t2u and u2t versions
3425 : ! NOTE: Input array includes ghost cells that must be updated before
3426 : ! calling this routine.
3427 : !
3428 : ! author: T. Craig
3429 :
3430 23136 : subroutine grid_average_X2YF(dir,work1,wght1,work2,wght2)
3431 :
3432 : use ice_constants, only: c0, p25, p5
3433 :
3434 : character(len=*) , intent(in) :: &
3435 : dir
3436 :
3437 : real (kind=dbl_kind), intent(in) :: &
3438 : work1(:,:,:), & ! LCOV_EXCL_LINE
3439 : wght1(:,:,:), & ! LCOV_EXCL_LINE
3440 : wght2(:,:,:)
3441 :
3442 : real (kind=dbl_kind), intent(out) :: &
3443 : work2(:,:,:)
3444 :
3445 : ! local variables
3446 :
3447 : integer (kind=int_kind) :: &
3448 : i, j, iblk, & ! LCOV_EXCL_LINE
3449 : ilo,ihi,jlo,jhi ! beginning and end of physical domain
3450 :
3451 : type (block) :: &
3452 : this_block ! block information for current block
3453 :
3454 : character(len=*), parameter :: subname = '(grid_average_X2YF)'
3455 :
3456 64887936 : work2(:,:,:) = c0
3457 :
3458 46272 : select case (trim(dir))
3459 :
3460 : case('NE')
3461 11520 : !$OMP PARALLEL DO PRIVATE(iblk,i,j,ilo,ihi,jlo,jhi,this_block)
3462 23184 : do iblk = 1, nblocks
3463 11568 : this_block = get_block(blocks_ice(iblk),iblk)
3464 11568 : ilo = this_block%ilo
3465 11568 : ihi = this_block%ihi
3466 11568 : jlo = this_block%jlo
3467 11568 : jhi = this_block%jhi
3468 397344 : do j = jlo, jhi
3469 12739056 : do i = ilo, ihi
3470 : work2(i,j,iblk) = p25 * &
3471 : (work1(i ,j ,iblk)*wght1(i ,j ,iblk) & ! LCOV_EXCL_LINE
3472 : + work1(i+1,j ,iblk)*wght1(i+1,j ,iblk) & ! LCOV_EXCL_LINE
3473 : + work1(i ,j+1,iblk)*wght1(i ,j+1,iblk) & ! LCOV_EXCL_LINE
3474 : + work1(i+1,j+1,iblk)*wght1(i+1,j+1,iblk)) & ! LCOV_EXCL_LINE
3475 12727488 : / wght2(i ,j ,iblk)
3476 : enddo
3477 : enddo
3478 : enddo
3479 : !$OMP END PARALLEL DO
3480 :
3481 : case('SW')
3482 5760 : !$OMP PARALLEL DO PRIVATE(iblk,i,j,ilo,ihi,jlo,jhi,this_block)
3483 23184 : do iblk = 1, nblocks
3484 11568 : this_block = get_block(blocks_ice(iblk),iblk)
3485 11568 : ilo = this_block%ilo
3486 11568 : ihi = this_block%ihi
3487 11568 : jlo = this_block%jlo
3488 11568 : jhi = this_block%jhi
3489 397344 : do j = jlo, jhi
3490 12739056 : do i = ilo, ihi
3491 : work2(i,j,iblk) = p25 * &
3492 : (work1(i ,j ,iblk)*wght1(i ,j ,iblk) & ! LCOV_EXCL_LINE
3493 : + work1(i-1,j ,iblk)*wght1(i-1,j ,iblk) & ! LCOV_EXCL_LINE
3494 : + work1(i ,j-1,iblk)*wght1(i ,j-1,iblk) & ! LCOV_EXCL_LINE
3495 : + work1(i-1,j-1,iblk)*wght1(i-1,j-1,iblk)) & ! LCOV_EXCL_LINE
3496 12727488 : / wght2(i ,j ,iblk)
3497 : enddo
3498 : enddo
3499 : enddo
3500 : !$OMP END PARALLEL DO
3501 :
3502 : case('NW')
3503 0 : !$OMP PARALLEL DO PRIVATE(iblk,i,j,ilo,ihi,jlo,jhi,this_block)
3504 0 : do iblk = 1, nblocks
3505 0 : this_block = get_block(blocks_ice(iblk),iblk)
3506 0 : ilo = this_block%ilo
3507 0 : ihi = this_block%ihi
3508 0 : jlo = this_block%jlo
3509 0 : jhi = this_block%jhi
3510 0 : do j = jlo, jhi
3511 0 : do i = ilo, ihi
3512 : work2(i,j,iblk) = p25 * &
3513 : (work1(i-1,j ,iblk)*wght1(i-1,j ,iblk) & ! LCOV_EXCL_LINE
3514 : + work1(i ,j ,iblk)*wght1(i ,j ,iblk) & ! LCOV_EXCL_LINE
3515 : + work1(i-1,j+1,iblk)*wght1(i-1,j+1,iblk) & ! LCOV_EXCL_LINE
3516 : + work1(i ,j+1,iblk)*wght1(i ,j+1,iblk)) & ! LCOV_EXCL_LINE
3517 0 : / wght2(i ,j ,iblk)
3518 : enddo
3519 : enddo
3520 : enddo
3521 : !$OMP END PARALLEL DO
3522 :
3523 : case('SE')
3524 0 : !$OMP PARALLEL DO PRIVATE(iblk,i,j,ilo,ihi,jlo,jhi,this_block)
3525 0 : do iblk = 1, nblocks
3526 0 : this_block = get_block(blocks_ice(iblk),iblk)
3527 0 : ilo = this_block%ilo
3528 0 : ihi = this_block%ihi
3529 0 : jlo = this_block%jlo
3530 0 : jhi = this_block%jhi
3531 0 : do j = jlo, jhi
3532 0 : do i = ilo, ihi
3533 : work2(i,j,iblk) = p25 * &
3534 : (work1(i ,j-1,iblk)*wght1(i ,j-1,iblk) & ! LCOV_EXCL_LINE
3535 : + work1(i+1,j-1,iblk)*wght1(i+1,j-1,iblk) & ! LCOV_EXCL_LINE
3536 : + work1(i ,j ,iblk)*wght1(i ,j ,iblk) & ! LCOV_EXCL_LINE
3537 : + work1(i+1,j ,iblk)*wght1(i+1,j ,iblk)) & ! LCOV_EXCL_LINE
3538 0 : / wght2(i ,j ,iblk)
3539 : enddo
3540 : enddo
3541 : enddo
3542 : !$OMP END PARALLEL DO
3543 :
3544 : case('E')
3545 0 : !$OMP PARALLEL DO PRIVATE(iblk,i,j,ilo,ihi,jlo,jhi,this_block)
3546 0 : do iblk = 1, nblocks
3547 0 : this_block = get_block(blocks_ice(iblk),iblk)
3548 0 : ilo = this_block%ilo
3549 0 : ihi = this_block%ihi
3550 0 : jlo = this_block%jlo
3551 0 : jhi = this_block%jhi
3552 0 : do j = jlo, jhi
3553 0 : do i = ilo, ihi
3554 : work2(i,j,iblk) = p5 * &
3555 : (work1(i ,j,iblk)*wght1(i ,j,iblk) & ! LCOV_EXCL_LINE
3556 : + work1(i+1,j,iblk)*wght1(i+1,j,iblk)) & ! LCOV_EXCL_LINE
3557 0 : / wght2(i ,j,iblk)
3558 : enddo
3559 : enddo
3560 : enddo
3561 : !$OMP END PARALLEL DO
3562 :
3563 : case('W')
3564 0 : !$OMP PARALLEL DO PRIVATE(iblk,i,j,ilo,ihi,jlo,jhi,this_block)
3565 0 : do iblk = 1, nblocks
3566 0 : this_block = get_block(blocks_ice(iblk),iblk)
3567 0 : ilo = this_block%ilo
3568 0 : ihi = this_block%ihi
3569 0 : jlo = this_block%jlo
3570 0 : jhi = this_block%jhi
3571 0 : do j = jlo, jhi
3572 0 : do i = ilo, ihi
3573 : work2(i,j,iblk) = p5 * &
3574 : (work1(i-1,j,iblk)*wght1(i-1,j,iblk) & ! LCOV_EXCL_LINE
3575 : + work1(i ,j,iblk)*wght1(i ,j,iblk)) & ! LCOV_EXCL_LINE
3576 0 : / wght2(i ,j,iblk)
3577 : enddo
3578 : enddo
3579 : enddo
3580 : !$OMP END PARALLEL DO
3581 :
3582 : case('N')
3583 0 : !$OMP PARALLEL DO PRIVATE(iblk,i,j,ilo,ihi,jlo,jhi,this_block)
3584 0 : do iblk = 1, nblocks
3585 0 : this_block = get_block(blocks_ice(iblk),iblk)
3586 0 : ilo = this_block%ilo
3587 0 : ihi = this_block%ihi
3588 0 : jlo = this_block%jlo
3589 0 : jhi = this_block%jhi
3590 0 : do j = jlo, jhi
3591 0 : do i = ilo, ihi
3592 : work2(i,j,iblk) = p5 * &
3593 : (work1(i,j ,iblk)*wght1(i,j ,iblk) & ! LCOV_EXCL_LINE
3594 : + work1(i,j+1,iblk)*wght1(i,j+1,iblk)) & ! LCOV_EXCL_LINE
3595 0 : / wght2(i ,j,iblk)
3596 : enddo
3597 : enddo
3598 : enddo
3599 : !$OMP END PARALLEL DO
3600 :
3601 : case('S')
3602 0 : !$OMP PARALLEL DO PRIVATE(iblk,i,j,ilo,ihi,jlo,jhi,this_block)
3603 0 : do iblk = 1, nblocks
3604 0 : this_block = get_block(blocks_ice(iblk),iblk)
3605 0 : ilo = this_block%ilo
3606 0 : ihi = this_block%ihi
3607 0 : jlo = this_block%jlo
3608 0 : jhi = this_block%jhi
3609 0 : do j = jlo, jhi
3610 0 : do i = ilo, ihi
3611 : work2(i,j,iblk) = p5 * &
3612 : (work1(i,j-1,iblk)*wght1(i,j-1,iblk) & ! LCOV_EXCL_LINE
3613 : + work1(i,j ,iblk)*wght1(i,j ,iblk)) & ! LCOV_EXCL_LINE
3614 0 : / wght2(i ,j,iblk)
3615 : enddo
3616 : enddo
3617 : enddo
3618 : !$OMP END PARALLEL DO
3619 :
3620 : case default
3621 46272 : call abort_ice(subname//'ERROR: unknown option '//trim(dir))
3622 : end select
3623 :
3624 23136 : end subroutine grid_average_X2YF
3625 :
3626 : !=======================================================================
3627 : ! Shifts quantities from one grid to another
3628 : ! State masked version, simple weighted averager
3629 : ! NOTE: Input array includes ghost cells that must be updated before
3630 : ! calling this routine.
3631 : !
3632 : ! author: T. Craig
3633 :
3634 0 : subroutine grid_average_X2Y_2(dir,work1a,wght1a,mask1a,work1b,wght1b,mask1b,work2)
3635 :
3636 : use ice_constants, only: c0
3637 :
3638 : character(len=*) , intent(in) :: &
3639 : dir
3640 :
3641 : real (kind=dbl_kind), intent(in) :: &
3642 : work1a(:,:,:), work1b(:,:,:), & ! LCOV_EXCL_LINE
3643 : wght1a(:,:,:), wght1b(:,:,:), & ! LCOV_EXCL_LINE
3644 : mask1a(:,:,:), mask1b(:,:,:)
3645 :
3646 : real (kind=dbl_kind), intent(out) :: &
3647 : work2(:,:,:)
3648 :
3649 : ! local variables
3650 :
3651 : integer (kind=int_kind) :: &
3652 : i, j, iblk, & ! LCOV_EXCL_LINE
3653 : ilo,ihi,jlo,jhi ! beginning and end of physical domain
3654 :
3655 : real (kind=dbl_kind) :: &
3656 0 : wtmp
3657 :
3658 : type (block) :: &
3659 : this_block ! block information for current block
3660 :
3661 : character(len=*), parameter :: subname = '(grid_average_X2Y_2)'
3662 :
3663 0 : work2(:,:,:) = c0
3664 :
3665 0 : select case (trim(dir))
3666 :
3667 : case('NE2US')
3668 0 : !$OMP PARALLEL DO PRIVATE(iblk,i,j,ilo,ihi,jlo,jhi,this_block,wtmp)
3669 0 : do iblk = 1, nblocks
3670 0 : this_block = get_block(blocks_ice(iblk),iblk)
3671 0 : ilo = this_block%ilo
3672 0 : ihi = this_block%ihi
3673 0 : jlo = this_block%jlo
3674 0 : jhi = this_block%jhi
3675 0 : do j = jlo, jhi
3676 0 : do i = ilo, ihi
3677 : wtmp = (mask1a(i ,j ,iblk)*wght1a(i ,j ,iblk) &
3678 : + mask1a(i+1,j ,iblk)*wght1a(i+1,j ,iblk) & ! LCOV_EXCL_LINE
3679 : + mask1b(i ,j ,iblk)*wght1b(i ,j ,iblk) & ! LCOV_EXCL_LINE
3680 0 : + mask1b(i ,j+1,iblk)*wght1b(i ,j+1,iblk))
3681 0 : if (wtmp /= c0) &
3682 : work2(i,j,iblk) = (mask1a(i ,j ,iblk)*work1a(i ,j ,iblk)*wght1a(i ,j ,iblk) & ! LCOV_EXCL_LINE
3683 : + mask1a(i+1,j ,iblk)*work1a(i+1,j ,iblk)*wght1a(i+1,j ,iblk) & ! LCOV_EXCL_LINE
3684 : + mask1b(i ,j ,iblk)*work1b(i ,j ,iblk)*wght1b(i ,j ,iblk) & ! LCOV_EXCL_LINE
3685 : + mask1b(i ,j+1,iblk)*work1b(i ,j+1,iblk)*wght1b(i ,j+1,iblk)) & ! LCOV_EXCL_LINE
3686 0 : / wtmp
3687 : enddo
3688 : enddo
3689 : enddo
3690 : !$OMP END PARALLEL DO
3691 :
3692 : case('NE2TS')
3693 0 : !$OMP PARALLEL DO PRIVATE(iblk,i,j,ilo,ihi,jlo,jhi,this_block,wtmp)
3694 0 : do iblk = 1, nblocks
3695 0 : this_block = get_block(blocks_ice(iblk),iblk)
3696 0 : ilo = this_block%ilo
3697 0 : ihi = this_block%ihi
3698 0 : jlo = this_block%jlo
3699 0 : jhi = this_block%jhi
3700 0 : do j = jlo, jhi
3701 0 : do i = ilo, ihi
3702 : wtmp = (mask1a(i ,j-1,iblk)*wght1a(i ,j-1,iblk) &
3703 : + mask1a(i ,j ,iblk)*wght1a(i ,j ,iblk) & ! LCOV_EXCL_LINE
3704 : + mask1b(i-1,j ,iblk)*wght1b(i-1,j ,iblk) & ! LCOV_EXCL_LINE
3705 0 : + mask1b(i ,j ,iblk)*wght1b(i ,j ,iblk))
3706 0 : if (wtmp /= c0) &
3707 : work2(i,j,iblk) = (mask1a(i ,j-1,iblk)*work1a(i ,j-1,iblk)*wght1a(i ,j-1,iblk) & ! LCOV_EXCL_LINE
3708 : + mask1a(i ,j ,iblk)*work1a(i ,j ,iblk)*wght1a(i ,j ,iblk) & ! LCOV_EXCL_LINE
3709 : + mask1b(i-1,j ,iblk)*work1b(i-1,j ,iblk)*wght1b(i-1,j ,iblk) & ! LCOV_EXCL_LINE
3710 : + mask1b(i ,j ,iblk)*work1b(i ,j ,iblk)*wght1b(i ,j ,iblk)) & ! LCOV_EXCL_LINE
3711 0 : / wtmp
3712 : enddo
3713 : enddo
3714 : enddo
3715 : !$OMP END PARALLEL DO
3716 :
3717 : case('NE2UA')
3718 0 : !$OMP PARALLEL DO PRIVATE(iblk,i,j,ilo,ihi,jlo,jhi,this_block,wtmp)
3719 0 : do iblk = 1, nblocks
3720 0 : this_block = get_block(blocks_ice(iblk),iblk)
3721 0 : ilo = this_block%ilo
3722 0 : ihi = this_block%ihi
3723 0 : jlo = this_block%jlo
3724 0 : jhi = this_block%jhi
3725 0 : do j = jlo, jhi
3726 0 : do i = ilo, ihi
3727 : wtmp = (wght1a(i ,j ,iblk) &
3728 : + wght1a(i+1,j ,iblk) & ! LCOV_EXCL_LINE
3729 : + wght1b(i ,j ,iblk) & ! LCOV_EXCL_LINE
3730 0 : + wght1b(i ,j+1,iblk))
3731 0 : if (wtmp /= c0) &
3732 : work2(i,j,iblk) = (work1a(i ,j ,iblk)*wght1a(i ,j ,iblk) & ! LCOV_EXCL_LINE
3733 : + work1a(i+1,j ,iblk)*wght1a(i+1,j ,iblk) & ! LCOV_EXCL_LINE
3734 : + work1b(i ,j ,iblk)*wght1b(i ,j ,iblk) & ! LCOV_EXCL_LINE
3735 : + work1b(i ,j+1,iblk)*wght1b(i ,j+1,iblk)) & ! LCOV_EXCL_LINE
3736 0 : / wtmp
3737 : enddo
3738 : enddo
3739 : enddo
3740 : !$OMP END PARALLEL DO
3741 :
3742 : case('NE2TA')
3743 0 : !$OMP PARALLEL DO PRIVATE(iblk,i,j,ilo,ihi,jlo,jhi,this_block,wtmp)
3744 0 : do iblk = 1, nblocks
3745 0 : this_block = get_block(blocks_ice(iblk),iblk)
3746 0 : ilo = this_block%ilo
3747 0 : ihi = this_block%ihi
3748 0 : jlo = this_block%jlo
3749 0 : jhi = this_block%jhi
3750 0 : do j = jlo, jhi
3751 0 : do i = ilo, ihi
3752 : wtmp = (wght1a(i ,j-1,iblk) &
3753 : + wght1a(i ,j ,iblk) & ! LCOV_EXCL_LINE
3754 : + wght1b(i-1,j ,iblk) & ! LCOV_EXCL_LINE
3755 0 : + wght1b(i ,j ,iblk))
3756 0 : if (wtmp /= c0) &
3757 : work2(i,j,iblk) = (work1a(i ,j-1,iblk)*wght1a(i ,j-1,iblk) & ! LCOV_EXCL_LINE
3758 : + work1a(i ,j ,iblk)*wght1a(i ,j ,iblk) & ! LCOV_EXCL_LINE
3759 : + work1b(i-1,j ,iblk)*wght1b(i-1,j ,iblk) & ! LCOV_EXCL_LINE
3760 : + work1b(i ,j ,iblk)*wght1b(i ,j ,iblk)) & ! LCOV_EXCL_LINE
3761 0 : / wtmp
3762 : enddo
3763 : enddo
3764 : enddo
3765 : !$OMP END PARALLEL DO
3766 :
3767 : case default
3768 0 : call abort_ice(subname//'ERROR: unknown option '//trim(dir))
3769 : end select
3770 :
3771 0 : end subroutine grid_average_X2Y_2
3772 :
3773 : !=======================================================================
3774 : ! Compute the minimum of adjacent values of a field at specific indices,
3775 : ! depending on the grid location (U, E, N)
3776 : !
3777 0 : real(kind=dbl_kind) function grid_neighbor_min(field, i, j, grid_location) result(mini)
3778 :
3779 : real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) :: &
3780 : field ! field defined at T point
3781 :
3782 : integer (kind=int_kind), intent(in) :: &
3783 : i, j
3784 :
3785 : character(len=*), intent(in) :: &
3786 : grid_location ! grid location at which to compute the minumum (U, E, N)
3787 :
3788 : character(len=*), parameter :: subname = '(grid_neighbor_min)'
3789 :
3790 0 : select case (trim(grid_location))
3791 : case('U')
3792 0 : mini = min(field(i,j), field(i+1,j), field(i,j+1), field(i+1,j+1))
3793 : case('E')
3794 0 : mini = min(field(i,j), field(i+1,j))
3795 : case('N')
3796 0 : mini = min(field(i,j), field(i,j+1))
3797 : case default
3798 0 : call abort_ice(subname // ' unknown grid_location: ' // grid_location)
3799 : end select
3800 :
3801 0 : end function grid_neighbor_min
3802 :
3803 : !=======================================================================
3804 : ! Compute the maximum of adjacent values of a field at specific indices,
3805 : ! depending on the grid location (U, E, N)
3806 : !
3807 0 : real(kind=dbl_kind) function grid_neighbor_max(field, i, j, grid_location) result(maxi)
3808 :
3809 : real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) :: &
3810 : field ! field defined at T point
3811 :
3812 : integer (kind=int_kind), intent(in) :: &
3813 : i, j
3814 :
3815 : character(len=*), intent(in) :: &
3816 : grid_location ! grid location at which to compute the maximum (U, E, N)
3817 :
3818 :
3819 : character(len=*), parameter :: subname = '(grid_neighbor_max)'
3820 :
3821 0 : select case (trim(grid_location))
3822 : case('U')
3823 0 : maxi = max(field(i,j), field(i+1,j), field(i,j+1), field(i+1,j+1))
3824 : case('E')
3825 0 : maxi = max(field(i,j), field(i+1,j))
3826 : case('N')
3827 0 : maxi = max(field(i,j), field(i,j+1))
3828 : case default
3829 0 : call abort_ice(subname // ' unknown grid_location: ' // grid_location)
3830 : end select
3831 :
3832 0 : end function grid_neighbor_max
3833 :
3834 : !=======================================================================
3835 : ! The following code is used for obtaining the coordinates of the grid
3836 : ! vertices for CF-compliant netCDF history output. Approximate!
3837 : !=======================================================================
3838 :
3839 : ! These fields are only used for netcdf history output, and the
3840 : ! ghost cell values are not needed.
3841 : ! NOTE: Extrapolations were used: these fields are approximate!
3842 : !
3843 : ! authors: A. McLaren, Met Office
3844 : ! E. Hunke, LANL
3845 :
3846 37 : subroutine gridbox_corners
3847 :
3848 : use ice_blocks, only: nx_block, ny_block
3849 : use ice_constants, only: c0, c2, c360, &
3850 : field_loc_NEcorner, field_type_scalar
3851 : use ice_domain_size, only: max_blocks
3852 :
3853 : integer (kind=int_kind) :: &
3854 : i,j,iblk,icorner,& ! index counters ! LCOV_EXCL_LINE
3855 : ilo,ihi,jlo,jhi ! beginning and end of physical domain
3856 :
3857 : real (kind=dbl_kind), dimension(:,:), allocatable :: &
3858 37 : work_g2
3859 :
3860 : real (kind=dbl_kind), dimension (nx_block,ny_block,max_blocks) :: &
3861 27874 : work1
3862 :
3863 : real (kind=dbl_kind) :: &
3864 8 : rad_to_deg
3865 :
3866 : type (block) :: &
3867 : this_block ! block information for current block
3868 :
3869 : character(len=*), parameter :: subname = '(gridbox_corners)'
3870 :
3871 37 : call icepack_query_parameters(rad_to_deg_out=rad_to_deg)
3872 37 : call icepack_warnings_flush(nu_diag)
3873 37 : if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
3874 0 : file=__FILE__, line=__LINE__)
3875 :
3876 : !-------------------------------------------------------------
3877 : ! Get coordinates of grid boxes for each block as follows:
3878 : ! (1) SW corner, (2) SE corner, (3) NE corner, (4) NW corner
3879 : !-------------------------------------------------------------
3880 :
3881 538192 : latu_bounds(:,:,:,:) = c0
3882 538192 : lonu_bounds(:,:,:,:) = c0
3883 :
3884 20 : !$OMP PARALLEL DO PRIVATE(iblk,i,j,ilo,ihi,jlo,jhi,this_block)
3885 50 : do iblk = 1, nblocks
3886 33 : this_block = get_block(blocks_ice(iblk),iblk)
3887 33 : ilo = this_block%ilo
3888 33 : ihi = this_block%ihi
3889 33 : jlo = this_block%jlo
3890 33 : jhi = this_block%jhi
3891 :
3892 1210 : do j = jlo, jhi
3893 45541 : do i = ilo, ihi
3894 :
3895 44368 : latu_bounds(1,i,j,iblk)=TLAT(i ,j ,iblk)*rad_to_deg
3896 44368 : latu_bounds(2,i,j,iblk)=TLAT(i+1,j ,iblk)*rad_to_deg
3897 44368 : latu_bounds(3,i,j,iblk)=TLAT(i+1,j+1,iblk)*rad_to_deg
3898 44368 : latu_bounds(4,i,j,iblk)=TLAT(i ,j+1,iblk)*rad_to_deg
3899 :
3900 44368 : lonu_bounds(1,i,j,iblk)=TLON(i ,j ,iblk)*rad_to_deg
3901 44368 : lonu_bounds(2,i,j,iblk)=TLON(i+1,j ,iblk)*rad_to_deg
3902 44368 : lonu_bounds(3,i,j,iblk)=TLON(i+1,j+1,iblk)*rad_to_deg
3903 45508 : lonu_bounds(4,i,j,iblk)=TLON(i ,j+1,iblk)*rad_to_deg
3904 :
3905 : enddo
3906 : enddo
3907 : enddo
3908 : !$OMP END PARALLEL DO
3909 :
3910 : !----------------------------------------------------------------
3911 : ! extrapolate on global grid to get edge values
3912 : !----------------------------------------------------------------
3913 :
3914 37 : if (my_task == master_task) then
3915 7 : allocate(work_g2(nx_global,ny_global))
3916 : else
3917 30 : allocate(work_g2(1,1))
3918 : endif
3919 :
3920 111936 : work1(:,:,:) = latu_bounds(2,:,:,:)
3921 : ! work_g2 = c0
3922 :
3923 37 : call gather_global(work_g2, work1, master_task, distrb_info)
3924 37 : if (my_task == master_task) then
3925 843 : do j = 1, ny_global
3926 0 : work_g2(nx_global,j) = c2*work_g2(nx_global-1,j) &
3927 843 : - work_g2(nx_global-2,j)
3928 : enddo
3929 : endif
3930 0 : call scatter_global(work1, work_g2, &
3931 : master_task, distrb_info, & ! LCOV_EXCL_LINE
3932 37 : field_loc_NEcorner, field_type_scalar)
3933 111936 : latu_bounds(2,:,:,:) = work1(:,:,:)
3934 :
3935 111936 : work1(:,:,:) = latu_bounds(3,:,:,:)
3936 37 : call gather_global(work_g2, work1, master_task, distrb_info)
3937 37 : if (my_task == master_task) then
3938 763 : do i = 1, nx_global
3939 0 : work_g2(i,ny_global) = c2*work_g2(i,ny_global-1) &
3940 763 : - work_g2(i,ny_global-2)
3941 : enddo
3942 843 : do j = 1, ny_global
3943 0 : work_g2(nx_global,j) = c2*work_g2(nx_global-1,j) &
3944 843 : - work_g2(nx_global-2,j)
3945 : enddo
3946 : endif
3947 0 : call scatter_global(work1, work_g2, &
3948 : master_task, distrb_info, & ! LCOV_EXCL_LINE
3949 37 : field_loc_NEcorner, field_type_scalar)
3950 111936 : latu_bounds(3,:,:,:) = work1(:,:,:)
3951 :
3952 111936 : work1(:,:,:) = latu_bounds(4,:,:,:)
3953 37 : call gather_global(work_g2, work1, master_task, distrb_info)
3954 37 : if (my_task == master_task) then
3955 763 : do i = 1, nx_global
3956 0 : work_g2(i,ny_global) = c2*work_g2(i,ny_global-1) &
3957 763 : - work_g2(i,ny_global-2)
3958 : enddo
3959 : endif
3960 0 : call scatter_global(work1, work_g2, &
3961 : master_task, distrb_info, & ! LCOV_EXCL_LINE
3962 37 : field_loc_NEcorner, field_type_scalar)
3963 111936 : latu_bounds(4,:,:,:) = work1(:,:,:)
3964 :
3965 111936 : work1(:,:,:) = lonu_bounds(2,:,:,:)
3966 37 : call gather_global(work_g2, work1, master_task, distrb_info)
3967 37 : if (my_task == master_task) then
3968 843 : do j = 1, ny_global
3969 0 : work_g2(nx_global,j) = c2*work_g2(nx_global-1,j) &
3970 843 : - work_g2(nx_global-2,j)
3971 : enddo
3972 : endif
3973 0 : call scatter_global(work1, work_g2, &
3974 : master_task, distrb_info, & ! LCOV_EXCL_LINE
3975 37 : field_loc_NEcorner, field_type_scalar)
3976 111936 : lonu_bounds(2,:,:,:) = work1(:,:,:)
3977 :
3978 111936 : work1(:,:,:) = lonu_bounds(3,:,:,:)
3979 37 : call gather_global(work_g2, work1, master_task, distrb_info)
3980 37 : if (my_task == master_task) then
3981 763 : do i = 1, nx_global
3982 0 : work_g2(i,ny_global) = c2*work_g2(i,ny_global-1) &
3983 763 : - work_g2(i,ny_global-2)
3984 : enddo
3985 843 : do j = 1, ny_global
3986 0 : work_g2(nx_global,j) = c2*work_g2(nx_global-1,j) &
3987 843 : - work_g2(nx_global-2,j)
3988 : enddo
3989 : endif
3990 0 : call scatter_global(work1, work_g2, &
3991 : master_task, distrb_info, & ! LCOV_EXCL_LINE
3992 37 : field_loc_NEcorner, field_type_scalar)
3993 111936 : lonu_bounds(3,:,:,:) = work1(:,:,:)
3994 :
3995 111936 : work1(:,:,:) = lonu_bounds(4,:,:,:)
3996 37 : call gather_global(work_g2, work1, master_task, distrb_info)
3997 37 : if (my_task == master_task) then
3998 763 : do i = 1, nx_global
3999 0 : work_g2(i,ny_global) = c2*work_g2(i,ny_global-1) &
4000 763 : - work_g2(i,ny_global-2)
4001 : enddo
4002 : endif
4003 0 : call scatter_global(work1, work_g2, &
4004 : master_task, distrb_info, & ! LCOV_EXCL_LINE
4005 37 : field_loc_NEcorner, field_type_scalar)
4006 111936 : lonu_bounds(4,:,:,:) = work1(:,:,:)
4007 :
4008 37 : deallocate(work_g2)
4009 :
4010 : !----------------------------------------------------------------
4011 : ! Convert longitude to Degrees East >0 for history output
4012 : !----------------------------------------------------------------
4013 :
4014 37 : allocate(work_g2(nx_block,ny_block)) ! not used as global here
4015 : !OMP fails in this loop
4016 197 : do iblk = 1, nblocks
4017 837 : do icorner = 1, 4
4018 446600 : work_g2(:,:) = lont_bounds(icorner,:,:,iblk) + c360
4019 446600 : where (work_g2 > c360) work_g2 = work_g2 - c360
4020 446600 : where (work_g2 < c0 ) work_g2 = work_g2 + c360
4021 446600 : lont_bounds(icorner,:,:,iblk) = work_g2(:,:)
4022 446600 : work_g2(:,:) = lonu_bounds(icorner,:,:,iblk) + c360
4023 446600 : where (work_g2 > c360) work_g2 = work_g2 - c360
4024 446600 : where (work_g2 < c0 ) work_g2 = work_g2 + c360
4025 446760 : lonu_bounds(icorner,:,:,iblk) = work_g2(:,:)
4026 : enddo
4027 : enddo
4028 37 : deallocate(work_g2)
4029 :
4030 74 : end subroutine gridbox_corners
4031 :
4032 : !=======================================================================
4033 : ! The following code is used for obtaining the coordinates of the grid
4034 : ! vertices for CF-compliant netCDF history output. Approximate!
4035 : !=======================================================================
4036 :
4037 : ! These fields are only used for netcdf history output, and the
4038 : ! ghost cell values are not needed.
4039 : ! NOTE: Extrapolations were used: these fields are approximate!
4040 : !
4041 :
4042 37 : subroutine gridbox_edges
4043 :
4044 : use ice_blocks, only: nx_block, ny_block
4045 : use ice_constants, only: c0, c2, c360, &
4046 : field_loc_NEcorner, field_type_scalar
4047 : use ice_domain_size, only: max_blocks
4048 :
4049 : integer (kind=int_kind) :: &
4050 : i,j,iblk,icorner,& ! index counters ! LCOV_EXCL_LINE
4051 : ilo,ihi,jlo,jhi ! beginning and end of physical domain
4052 :
4053 : real (kind=dbl_kind), dimension(:,:), allocatable :: &
4054 37 : work_g2
4055 :
4056 : real (kind=dbl_kind), dimension (nx_block,ny_block,max_blocks) :: &
4057 27874 : work1
4058 :
4059 : real (kind=dbl_kind) :: &
4060 8 : rad_to_deg
4061 :
4062 : type (block) :: &
4063 : this_block ! block information for current block
4064 :
4065 : character(len=*), parameter :: subname = '(gridbox_edges)'
4066 :
4067 37 : call icepack_query_parameters(rad_to_deg_out=rad_to_deg)
4068 37 : call icepack_warnings_flush(nu_diag)
4069 37 : if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
4070 0 : file=__FILE__, line=__LINE__)
4071 :
4072 : !-------------------------------------------------------------
4073 : ! Get coordinates of grid boxes for each block as follows:
4074 : ! for N pt: (1) W edge, (2) E edge, (3) E edge j+1, (4) W edge j+1
4075 : ! for E pt: (1) S edge, (2) S edge i+1, (3) N edge, i+1 (4) N edge
4076 : !-------------------------------------------------------------
4077 :
4078 538192 : latn_bounds(:,:,:,:) = c0
4079 538192 : lonn_bounds(:,:,:,:) = c0
4080 538192 : late_bounds(:,:,:,:) = c0
4081 538192 : lone_bounds(:,:,:,:) = c0
4082 :
4083 20 : !$OMP PARALLEL DO PRIVATE(iblk,i,j,ilo,ihi,jlo,jhi,this_block)
4084 50 : do iblk = 1, nblocks
4085 33 : this_block = get_block(blocks_ice(iblk),iblk)
4086 33 : ilo = this_block%ilo
4087 33 : ihi = this_block%ihi
4088 33 : jlo = this_block%jlo
4089 33 : jhi = this_block%jhi
4090 :
4091 1210 : do j = jlo, jhi
4092 45541 : do i = ilo, ihi
4093 :
4094 44368 : latn_bounds(1,i,j,iblk)=ELAT(i-1,j ,iblk)*rad_to_deg
4095 44368 : latn_bounds(2,i,j,iblk)=ELAT(i ,j ,iblk)*rad_to_deg
4096 44368 : latn_bounds(3,i,j,iblk)=ELAT(i ,j+1,iblk)*rad_to_deg
4097 44368 : latn_bounds(4,i,j,iblk)=ELAT(i-1,j+1,iblk)*rad_to_deg
4098 :
4099 44368 : lonn_bounds(1,i,j,iblk)=ELON(i-1,j ,iblk)*rad_to_deg
4100 44368 : lonn_bounds(2,i,j,iblk)=ELON(i ,j ,iblk)*rad_to_deg
4101 44368 : lonn_bounds(3,i,j,iblk)=ELON(i ,j+1,iblk)*rad_to_deg
4102 44368 : lonn_bounds(4,i,j,iblk)=ELON(i-1,j+1,iblk)*rad_to_deg
4103 :
4104 44368 : late_bounds(1,i,j,iblk)=NLAT(i ,j-1,iblk)*rad_to_deg
4105 44368 : late_bounds(2,i,j,iblk)=NLAT(i+1,j-1,iblk)*rad_to_deg
4106 44368 : late_bounds(3,i,j,iblk)=NLAT(i+1,j ,iblk)*rad_to_deg
4107 44368 : late_bounds(4,i,j,iblk)=NLAT(i ,j ,iblk)*rad_to_deg
4108 :
4109 44368 : lone_bounds(1,i,j,iblk)=NLON(i ,j-1,iblk)*rad_to_deg
4110 44368 : lone_bounds(2,i,j,iblk)=NLON(i+1,j-1,iblk)*rad_to_deg
4111 44368 : lone_bounds(3,i,j,iblk)=NLON(i+1,j ,iblk)*rad_to_deg
4112 45508 : lone_bounds(4,i,j,iblk)=NLON(i ,j ,iblk)*rad_to_deg
4113 :
4114 : enddo
4115 : enddo
4116 : enddo
4117 : !$OMP END PARALLEL DO
4118 :
4119 : !----------------------------------------------------------------
4120 : ! extrapolate on global grid to get edge values
4121 : !----------------------------------------------------------------
4122 :
4123 37 : if (my_task == master_task) then
4124 7 : allocate(work_g2(nx_global,ny_global))
4125 : else
4126 30 : allocate(work_g2(1,1))
4127 : endif
4128 :
4129 : ! latn_bounds
4130 :
4131 111936 : work1(:,:,:) = latn_bounds(1,:,:,:)
4132 37 : call gather_global(work_g2, work1, master_task, distrb_info)
4133 37 : if (my_task == master_task) then
4134 843 : do j = 1, ny_global
4135 0 : work_g2(1,j) = c2*work_g2(2,j) &
4136 843 : - work_g2(3,j)
4137 : enddo
4138 : endif
4139 0 : call scatter_global(work1, work_g2, &
4140 : master_task, distrb_info, & ! LCOV_EXCL_LINE
4141 37 : field_loc_NEcorner, field_type_scalar)
4142 111936 : latn_bounds(1,:,:,:) = work1(:,:,:)
4143 :
4144 111936 : work1(:,:,:) = latn_bounds(3,:,:,:)
4145 37 : call gather_global(work_g2, work1, master_task, distrb_info)
4146 37 : if (my_task == master_task) then
4147 763 : do i = 1, nx_global
4148 0 : work_g2(i,ny_global) = c2*work_g2(i,ny_global-1) &
4149 763 : - work_g2(i,ny_global-2)
4150 : enddo
4151 : endif
4152 0 : call scatter_global(work1, work_g2, &
4153 : master_task, distrb_info, & ! LCOV_EXCL_LINE
4154 37 : field_loc_NEcorner, field_type_scalar)
4155 111936 : latn_bounds(3,:,:,:) = work1(:,:,:)
4156 :
4157 111936 : work1(:,:,:) = latn_bounds(4,:,:,:)
4158 37 : call gather_global(work_g2, work1, master_task, distrb_info)
4159 37 : if (my_task == master_task) then
4160 763 : do i = 1, nx_global
4161 0 : work_g2(i,ny_global) = c2*work_g2(i,ny_global-1) &
4162 763 : - work_g2(i,ny_global-2)
4163 : enddo
4164 843 : do j = 1, ny_global
4165 0 : work_g2(1,j) = c2*work_g2(2,j) &
4166 843 : - work_g2(3,j)
4167 : enddo
4168 : endif
4169 0 : call scatter_global(work1, work_g2, &
4170 : master_task, distrb_info, & ! LCOV_EXCL_LINE
4171 37 : field_loc_NEcorner, field_type_scalar)
4172 111936 : latn_bounds(4,:,:,:) = work1(:,:,:)
4173 :
4174 : ! lonn_bounds
4175 :
4176 111936 : work1(:,:,:) = lonn_bounds(1,:,:,:)
4177 37 : call gather_global(work_g2, work1, master_task, distrb_info)
4178 37 : if (my_task == master_task) then
4179 843 : do j = 1, ny_global
4180 0 : work_g2(1,j) = c2*work_g2(2,j) &
4181 843 : - work_g2(3,j)
4182 : enddo
4183 : endif
4184 0 : call scatter_global(work1, work_g2, &
4185 : master_task, distrb_info, & ! LCOV_EXCL_LINE
4186 37 : field_loc_NEcorner, field_type_scalar)
4187 111936 : lonn_bounds(1,:,:,:) = work1(:,:,:)
4188 :
4189 111936 : work1(:,:,:) = lonn_bounds(3,:,:,:)
4190 37 : call gather_global(work_g2, work1, master_task, distrb_info)
4191 37 : if (my_task == master_task) then
4192 763 : do i = 1, nx_global
4193 0 : work_g2(i,ny_global) = c2*work_g2(i,ny_global-1) &
4194 763 : - work_g2(i,ny_global-2)
4195 : enddo
4196 : endif
4197 0 : call scatter_global(work1, work_g2, &
4198 : master_task, distrb_info, & ! LCOV_EXCL_LINE
4199 37 : field_loc_NEcorner, field_type_scalar)
4200 111936 : lonn_bounds(3,:,:,:) = work1(:,:,:)
4201 :
4202 111936 : work1(:,:,:) = lonn_bounds(4,:,:,:)
4203 37 : call gather_global(work_g2, work1, master_task, distrb_info)
4204 37 : if (my_task == master_task) then
4205 763 : do i = 1, nx_global
4206 0 : work_g2(i,ny_global) = c2*work_g2(i,ny_global-1) &
4207 763 : - work_g2(i,ny_global-2)
4208 : enddo
4209 843 : do j = 1, ny_global
4210 0 : work_g2(1,j) = c2*work_g2(2,j) &
4211 843 : - work_g2(3,j)
4212 : enddo
4213 : endif
4214 0 : call scatter_global(work1, work_g2, &
4215 : master_task, distrb_info, & ! LCOV_EXCL_LINE
4216 37 : field_loc_NEcorner, field_type_scalar)
4217 111936 : lonn_bounds(4,:,:,:) = work1(:,:,:)
4218 :
4219 : ! late_bounds
4220 :
4221 111936 : work1(:,:,:) = late_bounds(1,:,:,:)
4222 37 : call gather_global(work_g2, work1, master_task, distrb_info)
4223 37 : if (my_task == master_task) then
4224 763 : do i = 1, nx_global
4225 0 : work_g2(i,1) = c2*work_g2(i,2) &
4226 763 : - work_g2(i,3)
4227 : enddo
4228 : endif
4229 0 : call scatter_global(work1, work_g2, &
4230 : master_task, distrb_info, & ! LCOV_EXCL_LINE
4231 37 : field_loc_NEcorner, field_type_scalar)
4232 111936 : late_bounds(1,:,:,:) = work1(:,:,:)
4233 :
4234 111936 : work1(:,:,:) = late_bounds(2,:,:,:)
4235 37 : call gather_global(work_g2, work1, master_task, distrb_info)
4236 37 : if (my_task == master_task) then
4237 763 : do i = 1, nx_global
4238 0 : work_g2(i,1) = c2*work_g2(i,2) &
4239 763 : - work_g2(i,3)
4240 : enddo
4241 843 : do j = 1, ny_global
4242 0 : work_g2(nx_global,j) = c2*work_g2(nx_global-1,j) &
4243 843 : - work_g2(nx_global-2,j)
4244 : enddo
4245 : endif
4246 0 : call scatter_global(work1, work_g2, &
4247 : master_task, distrb_info, & ! LCOV_EXCL_LINE
4248 37 : field_loc_NEcorner, field_type_scalar)
4249 111936 : late_bounds(2,:,:,:) = work1(:,:,:)
4250 :
4251 111936 : work1(:,:,:) = late_bounds(3,:,:,:)
4252 37 : call gather_global(work_g2, work1, master_task, distrb_info)
4253 37 : if (my_task == master_task) then
4254 843 : do j = 1, ny_global
4255 0 : work_g2(nx_global,j) = c2*work_g2(nx_global-1,j) &
4256 843 : - work_g2(nx_global-2,j)
4257 : enddo
4258 : endif
4259 0 : call scatter_global(work1, work_g2, &
4260 : master_task, distrb_info, & ! LCOV_EXCL_LINE
4261 37 : field_loc_NEcorner, field_type_scalar)
4262 111936 : late_bounds(3,:,:,:) = work1(:,:,:)
4263 :
4264 : ! lone_bounds
4265 :
4266 111936 : work1(:,:,:) = lone_bounds(1,:,:,:)
4267 37 : call gather_global(work_g2, work1, master_task, distrb_info)
4268 37 : if (my_task == master_task) then
4269 763 : do i = 1, nx_global
4270 0 : work_g2(i,1) = c2*work_g2(i,2) &
4271 763 : - work_g2(i,3)
4272 : enddo
4273 : endif
4274 0 : call scatter_global(work1, work_g2, &
4275 : master_task, distrb_info, & ! LCOV_EXCL_LINE
4276 37 : field_loc_NEcorner, field_type_scalar)
4277 111936 : lone_bounds(1,:,:,:) = work1(:,:,:)
4278 :
4279 111936 : work1(:,:,:) = lone_bounds(2,:,:,:)
4280 37 : call gather_global(work_g2, work1, master_task, distrb_info)
4281 37 : if (my_task == master_task) then
4282 763 : do i = 1, nx_global
4283 0 : work_g2(i,1) = c2*work_g2(i,2) &
4284 763 : - work_g2(i,3)
4285 : enddo
4286 843 : do j = 1, ny_global
4287 0 : work_g2(nx_global,j) = c2*work_g2(nx_global-1,j) &
4288 843 : - work_g2(nx_global-2,j)
4289 : enddo
4290 : endif
4291 0 : call scatter_global(work1, work_g2, &
4292 : master_task, distrb_info, & ! LCOV_EXCL_LINE
4293 37 : field_loc_NEcorner, field_type_scalar)
4294 111936 : lone_bounds(2,:,:,:) = work1(:,:,:)
4295 :
4296 111936 : work1(:,:,:) = lone_bounds(3,:,:,:)
4297 37 : call gather_global(work_g2, work1, master_task, distrb_info)
4298 37 : if (my_task == master_task) then
4299 843 : do j = 1, ny_global
4300 0 : work_g2(nx_global,j) = c2*work_g2(nx_global-1,j) &
4301 843 : - work_g2(nx_global-2,j)
4302 : enddo
4303 : endif
4304 0 : call scatter_global(work1, work_g2, &
4305 : master_task, distrb_info, & ! LCOV_EXCL_LINE
4306 37 : field_loc_NEcorner, field_type_scalar)
4307 111936 : lone_bounds(3,:,:,:) = work1(:,:,:)
4308 :
4309 37 : deallocate(work_g2)
4310 :
4311 : !----------------------------------------------------------------
4312 : ! Convert longitude to Degrees East >0 for history output
4313 : !----------------------------------------------------------------
4314 :
4315 37 : allocate(work_g2(nx_block,ny_block)) ! not used as global here
4316 : !OMP fails in this loop
4317 197 : do iblk = 1, nblocks
4318 837 : do icorner = 1, 4
4319 446600 : work_g2(:,:) = lonn_bounds(icorner,:,:,iblk) + c360
4320 446600 : where (work_g2 > c360) work_g2 = work_g2 - c360
4321 446600 : where (work_g2 < c0 ) work_g2 = work_g2 + c360
4322 446600 : lonn_bounds(icorner,:,:,iblk) = work_g2(:,:)
4323 446600 : work_g2(:,:) = lone_bounds(icorner,:,:,iblk) + c360
4324 446600 : where (work_g2 > c360) work_g2 = work_g2 - c360
4325 446600 : where (work_g2 < c0 ) work_g2 = work_g2 + c360
4326 446760 : lone_bounds(icorner,:,:,iblk) = work_g2(:,:)
4327 : enddo
4328 : enddo
4329 37 : deallocate(work_g2)
4330 :
4331 74 : end subroutine gridbox_edges
4332 :
4333 : !=======================================================================
4334 :
4335 : ! NOTE: Boundary conditions for fields on NW, SW, SE corners
4336 : ! have not been implemented; using NE corner location for all.
4337 : ! Extrapolations are also used: these fields are approximate!
4338 : !
4339 : ! authors: A. McLaren, Met Office
4340 : ! E. Hunke, LANL
4341 :
4342 42 : subroutine gridbox_verts(work_g,vbounds)
4343 :
4344 : use ice_blocks, only: nx_block, ny_block
4345 : use ice_constants, only: c0, c2, &
4346 : field_loc_NEcorner, field_type_scalar
4347 : use ice_domain_size, only: max_blocks
4348 :
4349 : real (kind=dbl_kind), dimension(:,:), intent(in) :: &
4350 : work_g
4351 :
4352 : real (kind=dbl_kind), dimension(4,nx_block,ny_block,max_blocks), intent(out) :: &
4353 : vbounds
4354 :
4355 : integer (kind=int_kind) :: &
4356 : i,j ! index counters
4357 :
4358 : real (kind=dbl_kind) :: &
4359 16 : rad_to_deg
4360 :
4361 : real (kind=dbl_kind), dimension(:,:), allocatable :: &
4362 42 : work_g2
4363 :
4364 : real (kind=dbl_kind), dimension (nx_block,ny_block,max_blocks) :: &
4365 55684 : work1
4366 :
4367 : character(len=*), parameter :: subname = '(gridbox_verts)'
4368 :
4369 42 : call icepack_query_parameters(rad_to_deg_out=rad_to_deg)
4370 42 : call icepack_warnings_flush(nu_diag)
4371 42 : if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
4372 0 : file=__FILE__, line=__LINE__)
4373 :
4374 42 : if (my_task == master_task) then
4375 10 : allocate(work_g2(nx_global,ny_global))
4376 : else
4377 32 : allocate(work_g2(1,1))
4378 : endif
4379 :
4380 : !-------------------------------------------------------------
4381 : ! Get coordinates of grid boxes for each block as follows:
4382 : ! (1) SW corner, (2) SE corner, (3) NE corner, (4) NW corner
4383 : !-------------------------------------------------------------
4384 :
4385 117266 : work_g2(:,:) = c0
4386 42 : if (my_task == master_task) then
4387 1160 : do j = 2, ny_global
4388 115010 : do i = 2, nx_global
4389 115000 : work_g2(i,j) = work_g(i-1,j-1) * rad_to_deg
4390 : enddo
4391 : enddo
4392 : ! extrapolate
4393 1170 : do j = 1, ny_global
4394 1170 : work_g2(1,j) = c2*work_g2(2,j) - work_g2(3,j)
4395 : enddo
4396 1010 : do i = 1, nx_global
4397 1010 : work_g2(i,1) = c2*work_g2(i,2) - work_g2(i,3)
4398 : enddo
4399 : endif
4400 0 : call scatter_global(work1, work_g2, &
4401 : master_task, distrb_info, & ! LCOV_EXCL_LINE
4402 42 : field_loc_NEcorner, field_type_scalar)
4403 147616 : vbounds(1,:,:,:) = work1(:,:,:)
4404 :
4405 117266 : work_g2(:,:) = c0
4406 42 : if (my_task == master_task) then
4407 1160 : do j = 2, ny_global
4408 116160 : do i = 1, nx_global
4409 116150 : work_g2(i,j) = work_g(i,j-1) * rad_to_deg
4410 : enddo
4411 : enddo
4412 : ! extrapolate
4413 1010 : do i = 1, nx_global
4414 1010 : work_g2(i,1) = (c2*work_g2(i,2) - work_g2(i,3))
4415 : enddo
4416 : endif
4417 0 : call scatter_global(work1, work_g2, &
4418 : master_task, distrb_info, & ! LCOV_EXCL_LINE
4419 42 : field_loc_NEcorner, field_type_scalar)
4420 147616 : vbounds(2,:,:,:) = work1(:,:,:)
4421 :
4422 117266 : work_g2(:,:) = c0
4423 42 : if (my_task == master_task) then
4424 1170 : do j = 1, ny_global
4425 117170 : do i = 1, nx_global
4426 117160 : work_g2(i,j) = work_g(i,j) * rad_to_deg
4427 : enddo
4428 : enddo
4429 : endif
4430 0 : call scatter_global(work1, work_g2, &
4431 : master_task, distrb_info, & ! LCOV_EXCL_LINE
4432 42 : field_loc_NEcorner, field_type_scalar)
4433 147616 : vbounds(3,:,:,:) = work1(:,:,:)
4434 :
4435 117266 : work_g2(:,:) = c0
4436 42 : if (my_task == master_task) then
4437 1170 : do j = 1, ny_global
4438 116010 : do i = 2, nx_global
4439 116000 : work_g2(i,j) = work_g(i-1,j ) * rad_to_deg
4440 : enddo
4441 : enddo
4442 : ! extrapolate
4443 1170 : do j = 1, ny_global
4444 1170 : work_g2(1,j) = c2*work_g2(2,j) - work_g2(3,j)
4445 : enddo
4446 : endif
4447 0 : call scatter_global(work1, work_g2, &
4448 : master_task, distrb_info, & ! LCOV_EXCL_LINE
4449 42 : field_loc_NEcorner, field_type_scalar)
4450 147616 : vbounds(4,:,:,:) = work1(:,:,:)
4451 :
4452 42 : deallocate (work_g2)
4453 :
4454 42 : end subroutine gridbox_verts
4455 :
4456 : !=======================================================================
4457 : ! ocean bathymetry for grounded sea ice (seabed stress) or icebergs
4458 : ! currently hardwired for 40 levels (gx3, gx1 grids)
4459 : ! should be read from a file instead (see subroutine read_seabedstress_bathy)
4460 :
4461 37 : subroutine get_bathymetry
4462 :
4463 : use ice_constants, only: c0
4464 :
4465 : integer (kind=int_kind) :: &
4466 : i, j, k, iblk ! loop indices
4467 :
4468 : integer (kind=int_kind), parameter :: &
4469 : nlevel = 40 ! number of layers (gx3 grid)
4470 :
4471 : real (kind=dbl_kind), dimension(nlevel) :: &
4472 328 : depth ! total depth, m
4473 :
4474 : real (kind=dbl_kind) :: &
4475 8 : puny
4476 :
4477 : logical (kind=log_kind) :: &
4478 : calc_dragio
4479 :
4480 : real (kind=dbl_kind), dimension(nlevel), parameter :: &
4481 : thick = (/ & ! ocean layer thickness, m ! LCOV_EXCL_LINE
4482 : 10.01244_dbl_kind, 10.11258_dbl_kind, 10.31682_dbl_kind, & ! LCOV_EXCL_LINE
4483 : 10.63330_dbl_kind, 11.07512_dbl_kind, 11.66145_dbl_kind, & ! LCOV_EXCL_LINE
4484 : 12.41928_dbl_kind, 13.38612_dbl_kind, 14.61401_dbl_kind, & ! LCOV_EXCL_LINE
4485 : 16.17561_dbl_kind, 18.17368_dbl_kind, 20.75558_dbl_kind, & ! LCOV_EXCL_LINE
4486 : 24.13680_dbl_kind, 28.63821_dbl_kind, 34.74644_dbl_kind, & ! LCOV_EXCL_LINE
4487 : 43.20857_dbl_kind, 55.16812_dbl_kind, 72.30458_dbl_kind, & ! LCOV_EXCL_LINE
4488 : 96.74901_dbl_kind, 130.0392_dbl_kind, 170.0489_dbl_kind, & ! LCOV_EXCL_LINE
4489 : 207.9933_dbl_kind, 233.5694_dbl_kind, 245.2719_dbl_kind, & ! LCOV_EXCL_LINE
4490 : 248.9804_dbl_kind, 249.8322_dbl_kind, 249.9787_dbl_kind, & ! LCOV_EXCL_LINE
4491 : 249.9979_dbl_kind, 249.9998_dbl_kind, 250.0000_dbl_kind, & ! LCOV_EXCL_LINE
4492 : 250.0000_dbl_kind, 250.0000_dbl_kind, 250.0000_dbl_kind, & ! LCOV_EXCL_LINE
4493 : 250.0000_dbl_kind, 250.0000_dbl_kind, 250.0000_dbl_kind, & ! LCOV_EXCL_LINE
4494 : 250.0000_dbl_kind, 250.0000_dbl_kind, 250.0000_dbl_kind, & ! LCOV_EXCL_LINE
4495 : 250.0000_dbl_kind /)
4496 :
4497 : character(len=*), parameter :: subname = '(get_bathymetry)'
4498 :
4499 37 : call icepack_query_parameters(puny_out=puny, calc_dragio_out=calc_dragio)
4500 37 : call icepack_warnings_flush(nu_diag)
4501 37 : if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
4502 0 : file=__FILE__, line=__LINE__)
4503 :
4504 37 : if (use_bathymetry) then
4505 :
4506 0 : call read_seabedstress_bathy
4507 :
4508 : else
4509 :
4510 : ! convert to total depth
4511 37 : depth(1) = thick(1)
4512 1480 : do k = 2, nlevel
4513 1480 : depth(k) = depth(k-1) + thick(k)
4514 : enddo
4515 :
4516 111936 : bathymetry = c0
4517 197 : do iblk = 1, nblocks
4518 5340 : do j = 1, ny_block
4519 111650 : do i = 1, nx_block
4520 106347 : k = min(nint(kmt(i,j,iblk)),nlevel)
4521 106347 : if (k > nlevel) call abort_ice(subname//' kmt gt nlevel error')
4522 111490 : if (k > 0) bathymetry(i,j,iblk) = depth(k)
4523 : enddo
4524 : enddo
4525 : enddo
4526 :
4527 : ! For consistency, set thickness_ocn_layer1 in Icepack if 'calc_dragio' is active
4528 37 : if (calc_dragio) then
4529 0 : call icepack_init_parameters(thickness_ocn_layer1_in=thick(1))
4530 0 : call icepack_warnings_flush(nu_diag)
4531 0 : if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
4532 0 : file=__FILE__, line=__LINE__)
4533 : endif
4534 :
4535 : endif ! bathymetry_file
4536 :
4537 37 : end subroutine get_bathymetry
4538 :
4539 : !=======================================================================
4540 : ! with use_bathymetry = false, vertical depth profile generated for max KMT
4541 : ! with use_bathymetry = true, expects to read in pop vert_grid file
4542 :
4543 0 : subroutine get_bathymetry_popfile
4544 :
4545 : integer (kind=int_kind) :: &
4546 : i, j, k, iblk ! loop indices
4547 :
4548 : integer (kind=int_kind) :: &
4549 : ntmp, nlevel , & ! number of levels (max KMT) ! LCOV_EXCL_LINE
4550 : k1 , & ! levels ! LCOV_EXCL_LINE
4551 : ierr , & ! error tag ! LCOV_EXCL_LINE
4552 : fid ! fid unit number
4553 :
4554 : real (kind=dbl_kind), dimension(:),allocatable :: &
4555 : depth , & ! total depth, m ! LCOV_EXCL_LINE
4556 0 : thick ! layer thickness, cm -> m
4557 :
4558 : logical (kind=log_kind) :: &
4559 : calc_dragio
4560 :
4561 : character(len=*), parameter :: subname = '(get_bathymetry_popfile)'
4562 :
4563 0 : ntmp = maxval(nint(KMT))
4564 0 : nlevel = global_maxval(ntmp,distrb_info)
4565 :
4566 0 : if (my_task==master_task) then
4567 0 : write(nu_diag,*) subname,' KMT max = ',nlevel
4568 : endif
4569 :
4570 0 : allocate(depth(nlevel),thick(nlevel))
4571 0 : thick = -999999.
4572 0 : depth = -999999.
4573 :
4574 0 : if (use_bathymetry) then
4575 :
4576 0 : write (nu_diag,*) subname,' Bathymetry file = ', trim(bathymetry_file)
4577 0 : if (my_task == master_task) then
4578 0 : call get_fileunit(fid)
4579 0 : open(fid,file=bathymetry_file,form='formatted',iostat=ierr)
4580 0 : if (ierr/=0) call abort_ice(subname//' open error')
4581 0 : do k = 1,nlevel
4582 0 : read(fid,*,iostat=ierr) thick(k)
4583 0 : if (ierr/=0) call abort_ice(subname//' read error')
4584 : enddo
4585 0 : call release_fileunit(fid)
4586 : endif
4587 :
4588 0 : call broadcast_array(thick,master_task)
4589 :
4590 : else
4591 :
4592 : ! create thickness profile
4593 0 : k1 = min(5,nlevel)
4594 0 : do k = 1,k1
4595 0 : thick(k) = max(10000._dbl_kind/float(nlevel),500._dbl_kind)
4596 : enddo
4597 0 : do k = k1+1,nlevel
4598 0 : thick(k) = min(thick(k-1)*1.2_dbl_kind,20000._dbl_kind)
4599 : enddo
4600 :
4601 : endif
4602 :
4603 : ! convert thick from cm to m
4604 0 : thick = thick / 100._dbl_kind
4605 :
4606 : ! convert to total depth
4607 0 : depth(1) = thick(1)
4608 0 : do k = 2, nlevel
4609 0 : depth(k) = depth(k-1) + thick(k)
4610 0 : if (depth(k) < 0.) call abort_ice(subname//' negative depth error')
4611 : enddo
4612 :
4613 0 : if (my_task==master_task) then
4614 0 : do k = 1,nlevel
4615 0 : write(nu_diag,'(2a,i6,2f13.7)') subname,' k, thick(m), depth(m) = ',k,thick(k),depth(k)
4616 : enddo
4617 : endif
4618 :
4619 0 : bathymetry = 0._dbl_kind
4620 0 : do iblk = 1, nblocks
4621 0 : do j = 1, ny_block
4622 0 : do i = 1, nx_block
4623 0 : k = nint(kmt(i,j,iblk))
4624 0 : if (k > nlevel) call abort_ice(subname//' kmt gt nlevel error')
4625 0 : if (k > 0) bathymetry(i,j,iblk) = depth(k)
4626 : enddo
4627 : enddo
4628 : enddo
4629 :
4630 : ! For consistency, set thickness_ocn_layer1 in Icepack if 'calc_dragio' is active
4631 0 : call icepack_query_parameters(calc_dragio_out=calc_dragio)
4632 0 : if (calc_dragio) then
4633 0 : call icepack_init_parameters(thickness_ocn_layer1_in=thick(1))
4634 : endif
4635 0 : call icepack_warnings_flush(nu_diag)
4636 0 : if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
4637 0 : file=__FILE__, line=__LINE__)
4638 :
4639 0 : deallocate(depth,thick)
4640 :
4641 0 : end subroutine get_bathymetry_popfile
4642 :
4643 : !=======================================================================
4644 :
4645 : ! Read bathymetry data for seabed stress calculation (grounding scheme for
4646 : ! landfast ice) in CICE stand-alone mode. When CICE is in coupled mode
4647 : ! (e.g. CICE-NEMO), hwater should be uptated at each time level so that
4648 : ! it varies with ocean dynamics.
4649 : !
4650 : ! author: Fred Dupont, CMC
4651 :
4652 0 : subroutine read_seabedstress_bathy
4653 :
4654 : ! use module
4655 : use ice_read_write
4656 : use ice_constants, only: field_loc_center, field_type_scalar
4657 :
4658 : ! local variables
4659 : integer (kind=int_kind) :: &
4660 : fid_init ! file id for netCDF init file
4661 :
4662 : character (char_len_long) :: & ! input data file names
4663 : fieldname
4664 :
4665 : logical (kind=log_kind) :: diag=.true.
4666 :
4667 : character(len=*), parameter :: subname = '(read_seabedstress_bathy)'
4668 :
4669 0 : if (my_task == master_task) then
4670 0 : write (nu_diag,*) ' '
4671 0 : write (nu_diag,*) 'Bathymetry file: ', trim(bathymetry_file)
4672 0 : call icepack_warnings_flush(nu_diag)
4673 : endif
4674 :
4675 0 : call ice_open_nc(bathymetry_file,fid_init)
4676 :
4677 0 : fieldname='Bathymetry'
4678 :
4679 0 : if (my_task == master_task) then
4680 0 : write(nu_diag,*) 'reading ',TRIM(fieldname)
4681 0 : call icepack_warnings_flush(nu_diag)
4682 : endif
4683 : call ice_read_nc(fid_init,1,fieldname,bathymetry,diag, &
4684 : field_loc=field_loc_center, & ! LCOV_EXCL_LINE
4685 0 : field_type=field_type_scalar)
4686 :
4687 0 : call ice_close_nc(fid_init)
4688 :
4689 0 : if (my_task == master_task) then
4690 0 : write(nu_diag,*) 'closing file ',TRIM(bathymetry_file)
4691 0 : call icepack_warnings_flush(nu_diag)
4692 : endif
4693 :
4694 0 : end subroutine read_seabedstress_bathy
4695 :
4696 : !=======================================================================
4697 :
4698 : end module ice_grid
4699 :
4700 : !=======================================================================
|
