Line data Source code
1 : !=======================================================================
2 : !
3 : ! Drivers for remapping and upwind ice transport
4 : !
5 : ! authors: Elizabeth C. Hunke and William H. Lipscomb, LANL
6 : !
7 : ! 2004: Revised by William Lipscomb from ice_transport_mpdata.
8 : ! Stripped out mpdata, retained upwind, and added block structure.
9 : ! 2006: Incorporated remap transport driver and renamed from
10 : ! ice_transport_upwind.
11 : ! 2011: ECH moved edgearea arrays into ice_transport_remap.F90
12 :
13 : module ice_transport_driver
14 :
15 : use ice_kinds_mod
16 : use ice_communicate, only: my_task, master_task
17 : use ice_constants, only: c0, c1, p5, &
18 : field_loc_center, & ! LCOV_EXCL_LINE
19 : field_type_scalar, field_type_vector, & ! LCOV_EXCL_LINE
20 : field_loc_NEcorner, & ! LCOV_EXCL_LINE
21 : field_loc_Nface, field_loc_Eface
22 : use ice_diagnostics, only: diagnostic_abort
23 : use ice_fileunits, only: nu_diag
24 : use ice_exit, only: abort_ice
25 : use icepack_intfc, only: icepack_warnings_flush, icepack_warnings_aborted
26 : use icepack_intfc, only: icepack_compute_tracers
27 : use icepack_intfc, only: icepack_query_tracer_flags, &
28 : icepack_query_tracer_sizes, icepack_query_tracer_indices, & ! LCOV_EXCL_LINE
29 : icepack_query_parameters
30 :
31 : implicit none
32 : private
33 : public :: init_transport, transport_remap, transport_upwind
34 :
35 : character (len=char_len), public :: &
36 : advection ! type of advection scheme used
37 : ! 'upwind' => 1st order donor cell scheme
38 : ! 'remap' => remapping scheme
39 :
40 : ! NOTE: For remapping, hice and hsno are considered tracers.
41 : ! ntrace is not equal to ntrcr!
42 :
43 : integer (kind=int_kind) :: &
44 : ntrace ! number of tracers in use
45 :
46 : integer (kind=int_kind), dimension(:), allocatable, public :: &
47 : tracer_type , & ! = 1, 2, or 3 (depends on 0, 1 or 2 other tracers) ! LCOV_EXCL_LINE
48 : depend ! tracer dependencies (see below)
49 :
50 : logical (kind=log_kind), dimension (:), allocatable, public :: &
51 : has_dependents ! true if a tracer has dependent tracers
52 :
53 : logical (kind=log_kind), public :: &
54 : conserv_check ! if true, check conservation
55 :
56 : integer (kind=int_kind), parameter :: &
57 : integral_order = 3 ! polynomial order of quadrature integrals
58 : ! linear=1, quadratic=2, cubic=3
59 :
60 : logical (kind=log_kind), parameter :: &
61 : l_dp_midpt = .true. ! if true, find departure points using
62 : ! corrected midpoint velocity
63 :
64 : !=======================================================================
65 :
66 : contains
67 :
68 : !=======================================================================
69 : !
70 : ! This subroutine is a wrapper for init_remap, which initializes the
71 : ! remapping transport scheme. If the model is run with upwind
72 : ! transport, no initializations are necessary.
73 : !
74 : ! authors William H. Lipscomb, LANL
75 :
76 37 : subroutine init_transport
77 :
78 : use ice_state, only: trcr_depend
79 : use ice_timers, only: ice_timer_start, ice_timer_stop, timer_advect
80 : use ice_transport_remap, only: init_remap
81 : use ice_grid, only: grid_ice
82 :
83 : integer (kind=int_kind) :: &
84 : k, nt, nt1 ! tracer indices
85 :
86 : integer (kind=int_kind) :: &
87 : ntrcr , nt_Tsfc , nt_qice , nt_qsno , & ! LCOV_EXCL_LINE
88 : nt_sice , nt_fbri , nt_iage , nt_FY , & ! LCOV_EXCL_LINE
89 : nt_alvl , nt_vlvl , & ! LCOV_EXCL_LINE
90 : nt_apnd , nt_hpnd , nt_ipnd , nt_fsd , & ! LCOV_EXCL_LINE
91 : nt_smice , nt_smliq , nt_rhos , nt_rsnw , & ! LCOV_EXCL_LINE
92 : nt_isosno, nt_isoice, nt_bgc_Nit
93 :
94 : character(len=*), parameter :: subname = '(init_transport)'
95 :
96 37 : call ice_timer_start(timer_advect) ! advection
97 :
98 37 : call icepack_query_tracer_sizes(ntrcr_out=ntrcr)
99 : call icepack_query_tracer_indices(nt_Tsfc_out=nt_Tsfc, nt_qice_out=nt_qice, &
100 : nt_qsno_out=nt_qsno, nt_sice_out=nt_sice, nt_fbri_out=nt_fbri, & ! LCOV_EXCL_LINE
101 : nt_iage_out=nt_iage, nt_FY_out=nt_FY, nt_fsd_out=nt_fsd, & ! LCOV_EXCL_LINE
102 : nt_alvl_out=nt_alvl, nt_vlvl_out=nt_vlvl, & ! LCOV_EXCL_LINE
103 : nt_apnd_out=nt_apnd, nt_hpnd_out=nt_hpnd, nt_ipnd_out=nt_ipnd, & ! LCOV_EXCL_LINE
104 : nt_smice_out=nt_smice, nt_smliq_out=nt_smliq, nt_rhos_out=nt_rhos, & ! LCOV_EXCL_LINE
105 : nt_rsnw_out=nt_rsnw, nt_bgc_Nit_out=nt_bgc_Nit, & ! LCOV_EXCL_LINE
106 37 : nt_isosno_out=nt_isosno, nt_isoice_out=nt_isoice)
107 37 : call icepack_warnings_flush(nu_diag)
108 37 : if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
109 0 : file=__FILE__, line=__LINE__)
110 :
111 37 : ntrace = 2 + ntrcr ! hice,hsno,trcr
112 :
113 37 : if (allocated(tracer_type)) deallocate(tracer_type)
114 37 : if (allocated(depend)) deallocate(depend)
115 37 : if (allocated(has_dependents)) deallocate(has_dependents)
116 :
117 0 : allocate (tracer_type (ntrace), &
118 : depend (ntrace), & ! LCOV_EXCL_LINE
119 37 : has_dependents(ntrace))
120 :
121 : ! define tracer dependency arrays
122 : ! see comments in remapping routine
123 :
124 111 : depend(1:2) = 0 ! hice, hsno
125 111 : tracer_type(1:2) = 1 ! no dependency
126 :
127 37 : k = 2
128 :
129 925 : do nt = 1, ntrcr
130 888 : depend(k+nt) = trcr_depend(nt) ! 0 for ice area tracers
131 : ! 1 for ice volume tracers
132 : ! 2 for snow volume tracers
133 888 : tracer_type(k+nt) = 2 ! depends on 1 other tracer
134 925 : if (trcr_depend(nt) == 0) then
135 148 : tracer_type(k+nt) = 1 ! depends on no other tracers
136 740 : elseif (trcr_depend(nt) > 2) then
137 111 : if (trcr_depend(trcr_depend(nt)-2) > 0) then
138 74 : tracer_type(k+nt) = 3 ! depends on 2 other tracers
139 : endif
140 : endif
141 : enddo
142 :
143 999 : has_dependents = .false.
144 999 : do nt = 1, ntrace
145 999 : if (depend(nt) > 0) then
146 740 : nt1 = depend(nt)
147 740 : has_dependents(nt1) = .true.
148 740 : if (nt1 > nt) then
149 : write(nu_diag,*) &
150 0 : 'Tracer nt2 =',nt,' depends on tracer nt1 =',nt1
151 : call abort_ice(subname// &
152 0 : 'ERROR: remap transport: Must have nt2 > nt1')
153 : endif
154 : endif
155 : enddo ! ntrace
156 :
157 : ! diagnostic output
158 37 : if (my_task == master_task) then
159 7 : write (nu_diag, *) 'tracer index depend type has_dependents'
160 7 : nt = 1
161 7 : write(nu_diag,1000) 'hi ',nt,depend(nt),tracer_type(nt),&
162 14 : has_dependents(nt)
163 7 : nt = 2
164 7 : write(nu_diag,1000) 'hs ',nt,depend(nt),tracer_type(nt),&
165 14 : has_dependents(nt)
166 7 : k=2
167 175 : do nt = k+1, k+ntrcr
168 168 : if (nt-k==nt_Tsfc) &
169 : write(nu_diag,1000) 'nt_Tsfc ',nt,depend(nt),tracer_type(nt),& ! LCOV_EXCL_LINE
170 14 : has_dependents(nt)
171 168 : if (nt-k==nt_qice) &
172 : write(nu_diag,1000) 'nt_qice ',nt,depend(nt),tracer_type(nt),& ! LCOV_EXCL_LINE
173 14 : has_dependents(nt)
174 168 : if (nt-k==nt_qsno) &
175 : write(nu_diag,1000) 'nt_qsno ',nt,depend(nt),tracer_type(nt),& ! LCOV_EXCL_LINE
176 14 : has_dependents(nt)
177 168 : if (nt-k==nt_sice) &
178 : write(nu_diag,1000) 'nt_sice ',nt,depend(nt),tracer_type(nt),& ! LCOV_EXCL_LINE
179 14 : has_dependents(nt)
180 168 : if (nt-k==nt_fbri) &
181 : write(nu_diag,1000) 'nt_fbri ',nt,depend(nt),tracer_type(nt),& ! LCOV_EXCL_LINE
182 14 : has_dependents(nt)
183 168 : if (nt-k==nt_iage) &
184 : write(nu_diag,1000) 'nt_iage ',nt,depend(nt),tracer_type(nt),& ! LCOV_EXCL_LINE
185 14 : has_dependents(nt)
186 168 : if (nt-k==nt_FY) &
187 : write(nu_diag,1000) 'nt_FY ',nt,depend(nt),tracer_type(nt),& ! LCOV_EXCL_LINE
188 14 : has_dependents(nt)
189 168 : if (nt-k==nt_alvl) &
190 : write(nu_diag,1000) 'nt_alvl ',nt,depend(nt),tracer_type(nt),& ! LCOV_EXCL_LINE
191 14 : has_dependents(nt)
192 168 : if (nt-k==nt_vlvl) &
193 : write(nu_diag,1000) 'nt_vlvl ',nt,depend(nt),tracer_type(nt),& ! LCOV_EXCL_LINE
194 14 : has_dependents(nt)
195 168 : if (nt-k==nt_apnd) &
196 : write(nu_diag,1000) 'nt_apnd ',nt,depend(nt),tracer_type(nt),& ! LCOV_EXCL_LINE
197 14 : has_dependents(nt)
198 168 : if (nt-k==nt_hpnd) &
199 : write(nu_diag,1000) 'nt_hpnd ',nt,depend(nt),tracer_type(nt),& ! LCOV_EXCL_LINE
200 14 : has_dependents(nt)
201 168 : if (nt-k==nt_ipnd) &
202 : write(nu_diag,1000) 'nt_ipnd ',nt,depend(nt),tracer_type(nt),& ! LCOV_EXCL_LINE
203 14 : has_dependents(nt)
204 168 : if (nt-k==nt_smice) &
205 : write(nu_diag,1000) 'nt_smice ',nt,depend(nt),tracer_type(nt),& ! LCOV_EXCL_LINE
206 14 : has_dependents(nt)
207 168 : if (nt-k==nt_smliq) &
208 : write(nu_diag,1000) 'nt_smliq ',nt,depend(nt),tracer_type(nt),& ! LCOV_EXCL_LINE
209 14 : has_dependents(nt)
210 168 : if (nt-k==nt_rhos) &
211 : write(nu_diag,1000) 'nt_rhos ',nt,depend(nt),tracer_type(nt),& ! LCOV_EXCL_LINE
212 14 : has_dependents(nt)
213 168 : if (nt-k==nt_rsnw) &
214 : write(nu_diag,1000) 'nt_rsnw ',nt,depend(nt),tracer_type(nt),& ! LCOV_EXCL_LINE
215 14 : has_dependents(nt)
216 168 : if (nt-k==nt_fsd) &
217 : write(nu_diag,1000) 'nt_fsd ',nt,depend(nt),tracer_type(nt),& ! LCOV_EXCL_LINE
218 14 : has_dependents(nt)
219 168 : if (nt-k==nt_isosno) &
220 : write(nu_diag,1000) 'nt_isosno ',nt,depend(nt),tracer_type(nt),& ! LCOV_EXCL_LINE
221 14 : has_dependents(nt)
222 168 : if (nt-k==nt_isoice) &
223 : write(nu_diag,1000) 'nt_isoice ',nt,depend(nt),tracer_type(nt),& ! LCOV_EXCL_LINE
224 14 : has_dependents(nt)
225 168 : if (nt-k==nt_bgc_Nit) &
226 : write(nu_diag,1000) 'nt_bgc_Nit ',nt,depend(nt),tracer_type(nt),& ! LCOV_EXCL_LINE
227 7 : has_dependents(nt)
228 : enddo
229 7 : write(nu_diag,*) ' '
230 : endif ! master_task
231 : 1000 format (1x,a,2x,i6,2x,i6,2x,i4,4x,l4)
232 :
233 37 : if (trim(advection)=='remap') call init_remap ! grid quantities
234 :
235 37 : call ice_timer_stop(timer_advect) ! advection
236 :
237 37 : end subroutine init_transport
238 :
239 : !=======================================================================
240 : !
241 : ! This subroutine solves the transport equations for one timestep
242 : ! using the conservative remapping scheme developed by John Dukowicz
243 : ! and John Baumgardner (DB) and modified for sea ice by William
244 : ! Lipscomb and Elizabeth Hunke.
245 : !
246 : ! This scheme preserves monotonicity of ice area and tracers. That is,
247 : ! it does not produce new extrema. It is second-order accurate in space,
248 : ! except where gradients are limited to preserve monotonicity.
249 : !
250 : ! authors William H. Lipscomb, LANL
251 :
252 5784 : subroutine transport_remap (dt)
253 :
254 : use ice_blocks, only: nx_block, ny_block
255 : use ice_boundary, only: ice_HaloUpdate
256 : use ice_global_reductions, only: global_sum, global_sum_prod
257 : use ice_domain, only: nblocks, distrb_info, blocks_ice, halo_info
258 : use ice_domain_size, only: ncat, max_blocks
259 : use ice_blocks, only: nx_block, ny_block, block, get_block, nghost
260 : use ice_state, only: aice0, aicen, vicen, vsnon, trcrn, &
261 : uvel, vvel, bound_state, uvelE, vvelN
262 : use ice_grid, only: tarea, grid_ice
263 : use ice_calendar, only: istep1
264 : use ice_timers, only: ice_timer_start, ice_timer_stop, &
265 : timer_advect, timer_bound
266 : use ice_transport_remap, only: horizontal_remap, make_masks
267 :
268 : real (kind=dbl_kind), intent(in) :: &
269 : dt ! time step
270 :
271 : ! local variables
272 :
273 : integer (kind=int_kind) :: &
274 : iblk , & ! block index ! LCOV_EXCL_LINE
275 : ilo,ihi,jlo,jhi, & ! beginning and end of physical domain ! LCOV_EXCL_LINE
276 : n , & ! ice category index ! LCOV_EXCL_LINE
277 : nt, nt1, nt2 ! tracer indices
278 :
279 : real (kind=dbl_kind), dimension (nx_block,ny_block,0:ncat,max_blocks) :: &
280 : aim , & ! mean ice category areas in each grid cell ! LCOV_EXCL_LINE
281 30048528 : aimask ! = 1. if ice is present, = 0. otherwise
282 :
283 : real (kind=dbl_kind), dimension (nx_block,ny_block,ntrace,ncat,max_blocks) :: &
284 : trm , & ! mean tracer values in each grid cell ! LCOV_EXCL_LINE
285 650751888 : trmask ! = 1. if tracer is present, = 0. otherwise
286 :
287 : logical (kind=log_kind) :: &
288 : ckflag ! if true, abort the model
289 :
290 : integer (kind=int_kind) :: &
291 : istop, jstop ! indices of grid cell where model aborts
292 :
293 : integer (kind=int_kind), dimension(0:ncat,max_blocks) :: &
294 11568 : icellsnc ! number of cells with ice
295 :
296 : integer (kind=int_kind), dimension(nx_block*ny_block,0:ncat,max_blocks) :: &
297 11568 : indxinc, indxjnc ! compressed i/j indices
298 :
299 : integer (kind=int_kind) :: &
300 : ntrcr ! number of tracers
301 :
302 : type (block) :: &
303 : this_block ! block information for current block
304 :
305 : ! variables related to optional bug checks
306 :
307 : logical (kind=log_kind), parameter :: &
308 : l_monotonicity_check = .false. ! if true, check monotonicity
309 :
310 : real (kind=dbl_kind), dimension(0:ncat) :: &
311 : asum_init , & ! initial global ice area ! LCOV_EXCL_LINE
312 18768 : asum_final ! final global ice area
313 :
314 : real (kind=dbl_kind), dimension(ntrace,ncat) :: &
315 : atsum_init , & ! initial global ice area*tracer ! LCOV_EXCL_LINE
316 204528 : atsum_final ! final global ice area*tracer
317 :
318 : real (kind=dbl_kind), dimension (:,:,:,:,:), allocatable :: &
319 : tmin , & ! local min tracer ! LCOV_EXCL_LINE
320 5784 : tmax ! local max tracer
321 :
322 : integer (kind=int_kind) :: &
323 : alloc_error
324 :
325 : real (kind=dbl_kind), dimension (nx_block,ny_block,max_blocks) :: &
326 5015568 : work1
327 :
328 : character(len=char_len_long) :: &
329 : fieldid
330 :
331 : character(len=*), parameter :: subname = '(transport_remap)'
332 :
333 5784 : call ice_timer_start(timer_advect) ! advection
334 5784 : call icepack_query_tracer_sizes(ntrcr_out=ntrcr)
335 5784 : call icepack_warnings_flush(nu_diag)
336 5784 : if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
337 0 : file=__FILE__, line=__LINE__)
338 :
339 : !-------------------------------------------------------------------
340 : ! Prepare for remapping.
341 : ! Initialize, update ghost cells, fill tracer arrays.
342 : !-------------------------------------------------------------------
343 :
344 5784 : ckflag = .false.
345 5784 : istop = 0
346 5784 : jstop = 0
347 :
348 : !-------------------------------------------------------------------
349 : ! Compute open water area in each grid cell.
350 : ! Note: An aggregate_area call is needed only if the open
351 : ! water area has changed since the previous call.
352 : ! Here we assume that aice0 is up to date.
353 : !-------------------------------------------------------------------
354 :
355 : ! !$OMP PARALLEL DO PRIVATE(i,j,iblk) SCHEDULE(runtime)
356 : ! do iblk = 1, nblocks
357 : ! do j = 1, ny_block
358 : ! do i = 1, nx_block
359 : ! call aggregate_area (ncat,
360 : ! aicen(i,j,:,iblk), &
361 : ! aice (i,j, iblk), & ! LCOV_EXCL_LINE
362 : ! aice0(i,j, iblk))
363 : ! enddo
364 : ! enddo
365 : ! enddo
366 : ! !$OMP END PARALLEL DO
367 :
368 : !-------------------------------------------------------------------
369 : ! Ghost cell updates for state variables.
370 : ! Commented out because ghost cells are updated after cleanup_itd.
371 : !-------------------------------------------------------------------
372 : ! call ice_timer_start(timer_bound)
373 :
374 : ! call ice_HaloUpdate (aice0, halo_info, &
375 : ! field_loc_center, field_type_scalar)
376 :
377 : ! call bound_state (aicen, &
378 : ! vicen, vsnon, & ! LCOV_EXCL_LINE
379 : ! ntrcr, trcrn)
380 :
381 : ! call ice_timer_stop(timer_bound)
382 :
383 : !-------------------------------------------------------------------
384 : ! Ghost cell updates for ice velocity.
385 : ! Commented out because ghost cell velocities are computed
386 : ! in ice_dyn_evp.
387 : !-------------------------------------------------------------------
388 :
389 : ! call ice_timer_start(timer_bound)
390 : ! call ice_HaloUpdate (uvel, halo_info, &
391 : ! field_loc_NEcorner, field_type_vector)
392 : ! call ice_HaloUpdate (vvel, halo_info, &
393 : ! field_loc_NEcorner, field_type_vector)
394 : ! call ice_timer_stop(timer_bound)
395 :
396 :
397 2880 : !$OMP PARALLEL DO PRIVATE(iblk) SCHEDULE(runtime)
398 8688 : do iblk = 1, nblocks
399 :
400 : !-------------------------------------------------------------------
401 : ! Fill arrays with fields to be remapped.
402 : !-------------------------------------------------------------------
403 :
404 : call state_to_tracers(nx_block, ny_block, &
405 : ntrcr, ntrace, & ! LCOV_EXCL_LINE
406 : aice0(:,:, iblk), aicen(:,:,:, iblk), & ! LCOV_EXCL_LINE
407 : trcrn(:,:,:,:,iblk), & ! LCOV_EXCL_LINE
408 : vicen(:,:,:, iblk), vsnon(:,:,:, iblk), & ! LCOV_EXCL_LINE
409 11568 : aim (:,:,:, iblk), trm (:,:,:,:,iblk))
410 :
411 : enddo
412 : !$OMP END PARALLEL DO
413 :
414 : !-------------------------------------------------------------------
415 : ! Optional conservation and monotonicity checks.
416 : !-------------------------------------------------------------------
417 :
418 5784 : if (conserv_check) then
419 :
420 : !-------------------------------------------------------------------
421 : ! Compute initial values of globally conserved quantities.
422 : !-------------------------------------------------------------------
423 :
424 0 : do n = 0, ncat
425 0 : asum_init(n) = global_sum(aim(:,:,n,:), distrb_info, &
426 0 : field_loc_center, tarea)
427 : enddo
428 :
429 0 : do n = 1, ncat
430 0 : do nt = 1, ntrace
431 0 : if (tracer_type(nt)==1) then ! does not depend on another tracer
432 0 : atsum_init(nt,n) = &
433 : global_sum_prod(trm(:,:,nt,n,:), aim(:,:,n,:), & ! LCOV_EXCL_LINE
434 : distrb_info, field_loc_center, & ! LCOV_EXCL_LINE
435 0 : tarea)
436 0 : elseif (tracer_type(nt)==2) then ! depends on another tracer
437 0 : nt1 = depend(nt)
438 0 : work1(:,:,:) = trm(:,:,nt,n,:)*trm(:,:,nt1,n,:)
439 0 : atsum_init(nt,n) = &
440 : global_sum_prod(work1(:,:,:), aim(:,:,n,:), & ! LCOV_EXCL_LINE
441 : distrb_info, field_loc_center, & ! LCOV_EXCL_LINE
442 0 : tarea)
443 0 : elseif (tracer_type(nt)==3) then ! depends on two tracers
444 0 : nt1 = depend(nt)
445 0 : nt2 = depend(nt1)
446 0 : work1(:,:,:) = trm(:,:,nt,n,:)*trm(:,:,nt1,n,:) &
447 0 : *trm(:,:,nt2,n,:)
448 0 : atsum_init(nt,n) = &
449 : global_sum_prod(work1(:,:,:), aim(:,:,n,:), & ! LCOV_EXCL_LINE
450 : distrb_info, field_loc_center, & ! LCOV_EXCL_LINE
451 0 : tarea)
452 : endif ! tracer_type
453 : enddo ! nt
454 : enddo ! n
455 :
456 : endif ! conserv_check
457 :
458 : if (l_monotonicity_check) then
459 :
460 : allocate(tmin(nx_block,ny_block,ntrace,ncat,max_blocks), &
461 : tmax(nx_block,ny_block,ntrace,ncat,max_blocks), & ! LCOV_EXCL_LINE
462 : STAT=alloc_error)
463 :
464 : if (alloc_error /= 0) &
465 : call abort_ice (subname//'ERROR: allocation error')
466 :
467 : tmin(:,:,:,:,:) = c0
468 : tmax(:,:,:,:,:) = c0
469 :
470 : !$OMP PARALLEL DO PRIVATE(iblk,ilo,ihi,jlo,jhi,this_block,n) SCHEDULE(runtime)
471 : do iblk = 1, nblocks
472 : this_block = get_block(blocks_ice(iblk),iblk)
473 : ilo = this_block%ilo
474 : ihi = this_block%ihi
475 : jlo = this_block%jlo
476 : jhi = this_block%jhi
477 :
478 : !-------------------------------------------------------------------
479 : ! Compute masks.
480 : ! Masks are used to prevent tracer values in cells without ice
481 : ! from being used in the monotonicity check.
482 : !-------------------------------------------------------------------
483 :
484 : call make_masks (nx_block, ny_block, &
485 : ilo, ihi, jlo, jhi, & ! LCOV_EXCL_LINE
486 : nghost, ntrace, & ! LCOV_EXCL_LINE
487 : has_dependents, & ! LCOV_EXCL_LINE
488 : icellsnc (:,iblk), & ! LCOV_EXCL_LINE
489 : indxinc(:,:,iblk), indxjnc(:,:, iblk), & ! LCOV_EXCL_LINE
490 : aim(:,:,:, iblk), aimask(:,:,:, iblk), & ! LCOV_EXCL_LINE
491 : trm(:,:,:,:,iblk), trmask(:,:,:,:,iblk))
492 :
493 : !-------------------------------------------------------------------
494 : ! Compute local max and min of tracer fields.
495 : !-------------------------------------------------------------------
496 :
497 : do n = 1, ncat
498 : call local_max_min &
499 : (nx_block, ny_block, & ! LCOV_EXCL_LINE
500 : ilo, ihi, jlo, jhi, & ! LCOV_EXCL_LINE
501 : trm (:,:,:,n,iblk), & ! LCOV_EXCL_LINE
502 : tmin(:,:,:,n,iblk), tmax (:,:,:,n,iblk), & ! LCOV_EXCL_LINE
503 : aimask(:,:,n,iblk), trmask(:,:,:,n,iblk))
504 : enddo
505 : enddo
506 : !$OMP END PARALLEL DO
507 :
508 : call ice_timer_start(timer_bound)
509 : call ice_HaloUpdate (tmin, halo_info, &
510 : field_loc_center, field_type_scalar)
511 : call ice_HaloUpdate (tmax, halo_info, &
512 : field_loc_center, field_type_scalar)
513 : call ice_timer_stop(timer_bound)
514 :
515 : !$OMP PARALLEL DO PRIVATE(iblk,ilo,ihi,jlo,jhi,this_block,n) SCHEDULE(runtime)
516 : do iblk = 1, nblocks
517 : this_block = get_block(blocks_ice(iblk),iblk)
518 : ilo = this_block%ilo
519 : ihi = this_block%ihi
520 : jlo = this_block%jlo
521 : jhi = this_block%jhi
522 :
523 : do n = 1, ncat
524 : call quasilocal_max_min (nx_block, ny_block, &
525 : ilo, ihi, jlo, jhi, & ! LCOV_EXCL_LINE
526 : tmin(:,:,:,n,iblk), & ! LCOV_EXCL_LINE
527 : tmax(:,:,:,n,iblk))
528 : enddo
529 : enddo
530 : !$OMP END PARALLEL DO
531 :
532 : endif ! l_monotonicity_check
533 :
534 : !-------------------------------------------------------------------
535 : ! Main remapping routine: Step ice area and tracers forward in time.
536 : !-------------------------------------------------------------------
537 :
538 5784 : if (grid_ice == 'CD' .or. grid_ice == 'C') then
539 : call horizontal_remap (dt, ntrace, &
540 : uvel (:,:,:), vvel (:,:,:), & ! LCOV_EXCL_LINE
541 : aim (:,:,:,:), trm(:,:,:,:,:), & ! LCOV_EXCL_LINE
542 : tracer_type, depend, & ! LCOV_EXCL_LINE
543 : has_dependents, integral_order, & ! LCOV_EXCL_LINE
544 : l_dp_midpt, & ! LCOV_EXCL_LINE
545 0 : uvelE (:,:,:), vvelN (:,:,:))
546 : else
547 : call horizontal_remap (dt, ntrace, &
548 : uvel (:,:,:), vvel (:,:,:), & ! LCOV_EXCL_LINE
549 : aim (:,:,:,:), trm(:,:,:,:,:), & ! LCOV_EXCL_LINE
550 : tracer_type, depend, & ! LCOV_EXCL_LINE
551 : has_dependents, integral_order, & ! LCOV_EXCL_LINE
552 5784 : l_dp_midpt)
553 : endif
554 :
555 : !-------------------------------------------------------------------
556 : ! Given new fields, recompute state variables.
557 : !-------------------------------------------------------------------
558 :
559 2880 : !$OMP PARALLEL DO PRIVATE(iblk) SCHEDULE(runtime)
560 8688 : do iblk = 1, nblocks
561 :
562 : call tracers_to_state (nx_block, ny_block, &
563 : ntrcr, ntrace, & ! LCOV_EXCL_LINE
564 : aim (:,:,:, iblk), trm (:,:,:,:,iblk), & ! LCOV_EXCL_LINE
565 : aice0(:,:, iblk), aicen(:,:,:, iblk), & ! LCOV_EXCL_LINE
566 : trcrn(:,:,:,:,iblk), & ! LCOV_EXCL_LINE
567 11568 : vicen(:,:,:, iblk), vsnon(:,:, :,iblk))
568 :
569 : enddo ! iblk
570 : !$OMP END PARALLEL DO
571 :
572 : !-------------------------------------------------------------------
573 : ! Ghost cell updates for state variables.
574 : !-------------------------------------------------------------------
575 :
576 5784 : call ice_timer_start(timer_bound)
577 :
578 : call bound_state (aicen, &
579 : vicen, vsnon, & ! LCOV_EXCL_LINE
580 5784 : ntrcr, trcrn)
581 :
582 5784 : call ice_timer_stop(timer_bound)
583 :
584 : !-------------------------------------------------------------------
585 : ! Optional conservation and monotonicity checks
586 : !-------------------------------------------------------------------
587 :
588 : !-------------------------------------------------------------------
589 : ! Compute final values of globally conserved quantities.
590 : ! Check global conservation of area and area*tracers. (Optional)
591 : !-------------------------------------------------------------------
592 :
593 5784 : if (conserv_check) then
594 :
595 0 : do n = 0, ncat
596 0 : asum_final(n) = global_sum(aim(:,:,n,:), distrb_info, &
597 0 : field_loc_center, tarea)
598 : enddo
599 :
600 0 : do n = 1, ncat
601 0 : do nt = 1, ntrace
602 0 : if (tracer_type(nt)==1) then ! does not depend on another tracer
603 0 : atsum_final(nt,n) = &
604 : global_sum_prod(trm(:,:,nt,n,:), aim(:,:,n,:), & ! LCOV_EXCL_LINE
605 : distrb_info, field_loc_center, & ! LCOV_EXCL_LINE
606 0 : tarea)
607 0 : elseif (tracer_type(nt)==2) then ! depends on another tracer
608 0 : nt1 = depend(nt)
609 0 : work1(:,:,:) = trm(:,:,nt,n,:)*trm(:,:,nt1,n,:)
610 0 : atsum_final(nt,n) = &
611 : global_sum_prod(work1(:,:,:), aim(:,:,n,:), & ! LCOV_EXCL_LINE
612 : distrb_info, field_loc_center, & ! LCOV_EXCL_LINE
613 0 : tarea)
614 0 : elseif (tracer_type(nt)==3) then ! depends on two tracers
615 0 : nt1 = depend(nt)
616 0 : nt2 = depend(nt1)
617 0 : work1(:,:,:) = trm(:,:,nt,n,:)*trm(:,:,nt1,n,:) &
618 0 : *trm(:,:,nt2,n,:)
619 0 : atsum_final(nt,n) = &
620 : global_sum_prod(work1(:,:,:), aim(:,:,n,:), & ! LCOV_EXCL_LINE
621 : distrb_info, field_loc_center, & ! LCOV_EXCL_LINE
622 0 : tarea)
623 : endif ! tracer_type
624 : enddo ! nt
625 : enddo ! n
626 :
627 0 : if (my_task == master_task) then
628 0 : fieldid = subname//':000'
629 : call global_conservation (ckflag, fieldid, &
630 0 : asum_init(0), asum_final(0))
631 :
632 0 : if (ckflag) then
633 0 : write (nu_diag,*) 'istep1, my_task =', &
634 0 : istep1, my_task
635 0 : write (nu_diag,*) 'transport: conservation error, cat 0'
636 0 : call abort_ice(subname//'ERROR: conservation error1')
637 : endif
638 :
639 0 : do n = 1, ncat
640 0 : write(fieldid,'(a,i3.3)') subname,n
641 : call global_conservation &
642 : (ckflag, fieldid, & ! LCOV_EXCL_LINE
643 : asum_init(n), asum_final(n), & ! LCOV_EXCL_LINE
644 0 : atsum_init(:,n), atsum_final(:,n))
645 :
646 0 : if (ckflag) then
647 0 : write (nu_diag,*) 'istep1, my_task, cat =', &
648 0 : istep1, my_task, n
649 0 : write (nu_diag,*) 'transport: conservation error, cat ',n
650 0 : call abort_ice(subname//'ERROR: conservation error2')
651 : endif
652 : enddo ! n
653 :
654 : endif ! my_task = master_task
655 :
656 : endif ! conserv_check
657 :
658 : !-------------------------------------------------------------------
659 : ! Check tracer monotonicity. (Optional)
660 : !-------------------------------------------------------------------
661 :
662 : if (l_monotonicity_check) then
663 : !$OMP PARALLEL DO PRIVATE(iblk,ilo,ihi,jlo,jhi,this_block,n,ckflag,istop,jstop) SCHEDULE(runtime)
664 : do iblk = 1, nblocks
665 : this_block = get_block(blocks_ice(iblk),iblk)
666 : ilo = this_block%ilo
667 : ihi = this_block%ihi
668 : jlo = this_block%jlo
669 : jhi = this_block%jhi
670 :
671 : ckflag = .false.
672 : istop = 0
673 : jstop = 0
674 :
675 : do n = 1, ncat
676 : call check_monotonicity (nx_block, ny_block, &
677 : ilo, ihi, jlo, jhi, & ! LCOV_EXCL_LINE
678 : tmin(:,:,:,n,iblk), tmax(:,:,:,n,iblk), & ! LCOV_EXCL_LINE
679 : aim (:,:, n,iblk), trm (:,:,:,n,iblk), & ! LCOV_EXCL_LINE
680 : ckflag, & ! LCOV_EXCL_LINE
681 : istop, jstop)
682 :
683 : if (ckflag) then
684 : write (nu_diag,*) 'istep1, my_task, iblk, cat =', &
685 : istep1, my_task, iblk, n
686 : call diagnostic_abort(istop,jstop,iblk,' monotonicity error')
687 : endif
688 : enddo ! n
689 :
690 : enddo ! iblk
691 : !$OMP END PARALLEL DO
692 :
693 : deallocate(tmin, tmax, STAT=alloc_error)
694 : if (alloc_error /= 0) call abort_ice (subname//'ERROR: deallocation error')
695 :
696 : endif ! l_monotonicity_check
697 :
698 5784 : call ice_timer_stop(timer_advect) ! advection
699 :
700 5784 : end subroutine transport_remap
701 :
702 : !=======================================================================
703 : !
704 : ! Computes the transport equations for one timestep using upwind. Sets
705 : ! several fields into a work array and passes it to upwind routine.
706 :
707 0 : subroutine transport_upwind (dt)
708 :
709 : use ice_boundary, only: ice_HaloUpdate
710 : use ice_blocks, only: nx_block, ny_block, block, get_block, nx_block, ny_block
711 : use ice_domain, only: blocks_ice, halo_info, nblocks
712 : use ice_domain_size, only: ncat, max_blocks
713 : use ice_state, only: aice0, aicen, vicen, vsnon, trcrn, &
714 : uvel, vvel, trcr_depend, bound_state, trcr_base, & ! LCOV_EXCL_LINE
715 : n_trcr_strata, nt_strata, uvelE, vvelN
716 : use ice_flux, only: Tf
717 : use ice_grid, only: HTE, HTN, tarea, tmask, grid_ice
718 : use ice_timers, only: ice_timer_start, ice_timer_stop, &
719 : timer_bound, timer_advect
720 :
721 : real (kind=dbl_kind), intent(in) :: &
722 : dt ! time step
723 :
724 : ! local variables
725 :
726 : integer (kind=int_kind) :: &
727 : ntrcr , & ! ! LCOV_EXCL_LINE
728 : narr ! max number of state variable arrays
729 :
730 : integer (kind=int_kind) :: &
731 : i, j, iblk , & ! horizontal indices ! LCOV_EXCL_LINE
732 : ilo,ihi,jlo,jhi ! beginning and end of physical domain
733 :
734 : real (kind=dbl_kind), dimension (nx_block,ny_block,nblocks) :: &
735 0 : uee, vnn ! cell edge velocities
736 :
737 : real (kind=dbl_kind), dimension (:,:,:,:), allocatable :: &
738 0 : works ! work array
739 :
740 : type (block) :: &
741 : this_block ! block information for current block
742 :
743 : character(len=*), parameter :: subname = '(transport_upwind)'
744 :
745 0 : call ice_timer_start(timer_advect) ! advection
746 :
747 0 : call icepack_query_tracer_sizes(ntrcr_out=ntrcr)
748 0 : call icepack_warnings_flush(nu_diag)
749 0 : if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
750 0 : file=__FILE__, line=__LINE__)
751 :
752 0 : narr = 1 + ncat*(3+ntrcr) ! max number of state variable arrays
753 :
754 0 : allocate (works(nx_block,ny_block,narr,max_blocks))
755 :
756 : !-------------------------------------------------------------------
757 : ! Get ghost cell values of state variables.
758 : ! (Assume velocities are already known for ghost cells, also.)
759 : !-------------------------------------------------------------------
760 : ! call bound_state (aicen, &
761 : ! vicen, vsnon, & ! LCOV_EXCL_LINE
762 : ! ntrcr, trcrn)
763 :
764 : ! call ice_timer_start(timer_bound)
765 : ! call ice_HaloUpdate (uvel, halo_info, &
766 : ! field_loc_NEcorner, field_type_vector)
767 : ! call ice_HaloUpdate (vvel, halo_info, &
768 : ! field_loc_NEcorner, field_type_vector)
769 : ! call ice_timer_stop(timer_bound)
770 :
771 : !-------------------------------------------------------------------
772 : ! Average corner velocities to edges.
773 : !-------------------------------------------------------------------
774 0 : if (grid_ice == 'CD' .or. grid_ice == 'C') then
775 0 : uee(:,:,:)=uvelE(:,:,:)
776 0 : vnn(:,:,:)=vvelN(:,:,:)
777 : else
778 0 : !$OMP PARALLEL DO PRIVATE(iblk,i,j,ilo,ihi,jlo,jhi,this_block) SCHEDULE(runtime)
779 0 : do iblk = 1, nblocks
780 0 : this_block = get_block(blocks_ice(iblk),iblk)
781 0 : ilo = this_block%ilo
782 0 : ihi = this_block%ihi
783 0 : jlo = this_block%jlo
784 0 : jhi = this_block%jhi
785 :
786 0 : do j = jlo, jhi
787 0 : do i = ilo, ihi
788 0 : uee(i,j,iblk) = p5*(uvel(i,j,iblk) + uvel(i ,j-1,iblk))
789 0 : vnn(i,j,iblk) = p5*(vvel(i,j,iblk) + vvel(i-1,j ,iblk))
790 : enddo
791 : enddo
792 : enddo
793 : !$OMP END PARALLEL DO
794 :
795 0 : call ice_timer_start(timer_bound)
796 0 : call ice_HaloUpdate (uee, halo_info, &
797 0 : field_loc_Eface, field_type_vector)
798 0 : call ice_HaloUpdate (vnn, halo_info, &
799 0 : field_loc_Nface, field_type_vector)
800 0 : call ice_timer_stop(timer_bound)
801 : endif
802 :
803 0 : !$OMP PARALLEL DO PRIVATE(iblk,ilo,ihi,jlo,jhi,this_block) SCHEDULE(runtime)
804 0 : do iblk = 1, nblocks
805 0 : this_block = get_block(blocks_ice(iblk),iblk)
806 0 : ilo = this_block%ilo
807 0 : ihi = this_block%ihi
808 0 : jlo = this_block%jlo
809 0 : jhi = this_block%jhi
810 :
811 : !-----------------------------------------------------------------
812 : ! fill work arrays with fields to be advected
813 : !-----------------------------------------------------------------
814 :
815 : call state_to_work (nx_block, ny_block, &
816 : ntrcr, & ! LCOV_EXCL_LINE
817 : narr, trcr_depend, & ! LCOV_EXCL_LINE
818 : aicen (:,:, :,iblk), trcrn (:,:,:,:,iblk), & ! LCOV_EXCL_LINE
819 : vicen (:,:, :,iblk), vsnon (:,:, :,iblk), & ! LCOV_EXCL_LINE
820 0 : aice0 (:,:, iblk), works (:,:, :,iblk))
821 :
822 : !-----------------------------------------------------------------
823 : ! advect
824 : !-----------------------------------------------------------------
825 :
826 : call upwind_field (nx_block, ny_block, &
827 : ilo, ihi, jlo, jhi, & ! LCOV_EXCL_LINE
828 : dt, & ! LCOV_EXCL_LINE
829 : narr, works(:,:,:,iblk), & ! LCOV_EXCL_LINE
830 : uee (:,:,iblk), vnn (:,:,iblk), & ! LCOV_EXCL_LINE
831 : HTE (:,:,iblk), HTN (:,:,iblk), & ! LCOV_EXCL_LINE
832 0 : tarea(:,:,iblk))
833 :
834 : !-----------------------------------------------------------------
835 : ! convert work arrays back to state variables
836 : !-----------------------------------------------------------------
837 :
838 : call work_to_state (nx_block, ny_block, &
839 : ntrcr, narr, & ! LCOV_EXCL_LINE
840 : trcr_depend(:), trcr_base(:,:), & ! LCOV_EXCL_LINE
841 : n_trcr_strata(:), nt_strata(:,:), & ! LCOV_EXCL_LINE
842 : tmask(:,:, iblk), Tf (:,:,iblk), & ! LCOV_EXCL_LINE
843 : aicen(:,:, :,iblk), trcrn (:,:,:,:,iblk), & ! LCOV_EXCL_LINE
844 : vicen(:,:, :,iblk), vsnon (:,:, :,iblk), & ! LCOV_EXCL_LINE
845 0 : aice0(:,:, iblk), works (:,:, :,iblk))
846 :
847 : enddo ! iblk
848 : !$OMP END PARALLEL DO
849 :
850 0 : deallocate (works)
851 :
852 : !-------------------------------------------------------------------
853 : ! Ghost cell updates for state variables.
854 : !-------------------------------------------------------------------
855 :
856 0 : call ice_timer_start(timer_bound)
857 :
858 : call bound_state (aicen, &
859 : vicen, vsnon, & ! LCOV_EXCL_LINE
860 0 : ntrcr, trcrn)
861 :
862 0 : call ice_timer_stop(timer_bound)
863 :
864 0 : call ice_timer_stop(timer_advect) ! advection
865 :
866 0 : end subroutine transport_upwind
867 :
868 : !=======================================================================
869 : ! The next few subroutines (through check_monotonicity) are called
870 : ! by transport_remap.
871 : !=======================================================================
872 : !
873 : ! Fill ice area and tracer arrays.
874 : ! Assume that the advected tracers are hicen, hsnon, trcrn,
875 : ! qicen(1:nilyr), and qsnon(1:nslyr).
876 : ! This subroutine must be modified if a different set of tracers
877 : ! is to be transported. The rule for ordering tracers
878 : ! is that a dependent tracer (such as qice) must have a larger
879 : ! tracer index than the tracer it depends on (i.e., hice).
880 : !
881 : ! author William H. Lipscomb, LANL
882 :
883 22944 : subroutine state_to_tracers (nx_block, ny_block, &
884 : ntrcr, ntrace, & ! LCOV_EXCL_LINE
885 : aice0, aicen, & ! LCOV_EXCL_LINE
886 : trcrn, & ! LCOV_EXCL_LINE
887 : vicen, vsnon, & ! LCOV_EXCL_LINE
888 22944 : aim, trm)
889 :
890 : use ice_domain_size, only: ncat, nslyr
891 :
892 : integer (kind=int_kind), intent(in) :: &
893 : nx_block, ny_block, & ! block dimensions ! LCOV_EXCL_LINE
894 : ntrcr , & ! number of tracers in use ! LCOV_EXCL_LINE
895 : ntrace ! number of tracers in use incl. hi, hs
896 :
897 : real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) :: &
898 : aice0 ! fractional open water area
899 :
900 : real (kind=dbl_kind), dimension (nx_block,ny_block,ncat), intent(in) :: &
901 : aicen , & ! fractional ice area ! LCOV_EXCL_LINE
902 : vicen , & ! volume per unit area of ice (m) ! LCOV_EXCL_LINE
903 : vsnon ! volume per unit area of snow (m)
904 :
905 : real (kind=dbl_kind), dimension (nx_block,ny_block,ntrcr,ncat), intent(in) :: &
906 : trcrn ! ice area tracers
907 :
908 : real (kind=dbl_kind), dimension (nx_block,ny_block,0:ncat), intent(out) :: &
909 : aim ! mean ice area in each grid cell
910 :
911 : real (kind=dbl_kind), dimension (nx_block,ny_block,ntrace,ncat), intent(out) :: &
912 : trm ! mean tracer values in each grid cell
913 :
914 : ! local variables
915 :
916 : integer (kind=int_kind) :: &
917 : nt_qsno , & ! ! LCOV_EXCL_LINE
918 : i, j, n , & ! standard indices ! LCOV_EXCL_LINE
919 : it, kt , & ! tracer indices ! LCOV_EXCL_LINE
920 : ij ! combined i/j index
921 :
922 : real (kind=dbl_kind) :: &
923 : puny , & ! ! LCOV_EXCL_LINE
924 : rhos , & ! snow density (km/m^3) ! LCOV_EXCL_LINE
925 : Lfresh , & ! latent heat of melting fresh ice (J/kg) ! LCOV_EXCL_LINE
926 5760 : w1 ! work variable
927 :
928 : integer (kind=int_kind), dimension(nx_block*ny_block,0:ncat) :: &
929 : indxi , & ! compressed i/j indices ! LCOV_EXCL_LINE
930 45888 : indxj
931 :
932 : integer (kind=int_kind), dimension(0:ncat) :: &
933 45888 : icells ! number of cells with ice
934 :
935 : character(len=*), parameter :: subname = '(state_to_tracers)'
936 :
937 : call icepack_query_parameters(puny_out=puny, rhos_out=rhos, &
938 22944 : Lfresh_out=Lfresh)
939 22944 : call icepack_query_tracer_indices(nt_qsno_out=nt_qsno)
940 22944 : call icepack_warnings_flush(nu_diag)
941 22944 : if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
942 0 : file=__FILE__, line=__LINE__)
943 :
944 16186320 : aim(:,:,0) = aice0(:,:)
945 :
946 137664 : do n = 1, ncat
947 :
948 2104336320 : trm(:,:,:,n) = c0
949 :
950 : !-------------------------------------------------------------------
951 : ! Find grid cells where ice is present and fill area array.
952 : !-------------------------------------------------------------------
953 :
954 114720 : icells(n) = 0
955 3767880 : do j = 1, ny_block
956 80931600 : do i = 1, nx_block
957 77163720 : aim(i,j,n) = aicen(i,j,n)
958 80816880 : if (aim(i,j,n) > puny) then
959 37655263 : icells(n) = icells(n) + 1
960 37655263 : ij = icells(n)
961 37655263 : indxi(ij,n) = i
962 37655263 : indxj(ij,n) = j
963 : endif ! aim > puny
964 : enddo
965 : enddo
966 :
967 : !-------------------------------------------------------------------
968 : ! Fill tracer array
969 : ! Note: If aice > 0, then hice > 0, but we can have hsno = 0.
970 : ! Alse note: We transport qice*nilyr rather than qice, so as to
971 : ! avoid extra operations here and in tracers_to_state.
972 : !-------------------------------------------------------------------
973 :
974 37769983 : do ij = 1, icells(n)
975 37655263 : i = indxi(ij,n)
976 37655263 : j = indxj(ij,n)
977 37655263 : w1 = c1 / aim(i,j,n)
978 37655263 : trm(i,j,1,n) = vicen(i,j,n) * w1 ! hice
979 37769983 : trm(i,j,2,n) = vsnon(i,j,n) * w1 ! hsno
980 : enddo
981 114720 : kt = 2
982 :
983 2890944 : do it = 1, ntrcr
984 2868000 : if (it >= nt_qsno .and. it < nt_qsno+nslyr) then
985 37769983 : do ij = 1, icells(n)
986 37655263 : i = indxi(ij,n)
987 37655263 : j = indxj(ij,n)
988 37769983 : trm(i,j,kt+it,n) = trcrn(i,j,it,n) + rhos*Lfresh ! snow enthalpy
989 : enddo
990 : else
991 868709609 : do ij = 1, icells(n)
992 866071049 : i = indxi(ij,n)
993 866071049 : j = indxj(ij,n)
994 868709609 : trm(i,j,kt+it,n) = trcrn(i,j,it,n) ! other tracers
995 : enddo
996 : endif
997 : enddo
998 : enddo ! ncat
999 :
1000 22944 : end subroutine state_to_tracers
1001 :
1002 : !=======================================================================
1003 : !
1004 : ! Convert area and tracer arrays back to state variables.
1005 : !
1006 : ! author William H. Lipscomb, LANL
1007 :
1008 22944 : subroutine tracers_to_state (nx_block, ny_block, &
1009 : ntrcr, ntrace, & ! LCOV_EXCL_LINE
1010 : aim, trm, & ! LCOV_EXCL_LINE
1011 : aice0, aicen, & ! LCOV_EXCL_LINE
1012 : trcrn, & ! LCOV_EXCL_LINE
1013 22944 : vicen, vsnon)
1014 :
1015 : use ice_domain_size, only: ncat, nslyr
1016 :
1017 : integer (kind=int_kind), intent(in) :: &
1018 : nx_block, ny_block, & ! block dimensions ! LCOV_EXCL_LINE
1019 : ntrcr , & ! number of tracers in use ! LCOV_EXCL_LINE
1020 : ntrace ! number of tracers in use incl. hi, hs
1021 :
1022 : real (kind=dbl_kind), dimension (nx_block,ny_block,0:ncat), intent(in) :: &
1023 : aim ! fractional ice area
1024 :
1025 : real (kind=dbl_kind), dimension (nx_block,ny_block,ntrace,ncat), intent(in) :: &
1026 : trm ! mean tracer values in each grid cell
1027 :
1028 : real (kind=dbl_kind), dimension (nx_block,ny_block), intent(inout) :: &
1029 : aice0 ! fractional ice area
1030 :
1031 : real (kind=dbl_kind), dimension (nx_block,ny_block,ncat), intent(inout) :: &
1032 : aicen , & ! fractional ice area ! LCOV_EXCL_LINE
1033 : vicen , & ! volume per unit area of ice (m) ! LCOV_EXCL_LINE
1034 : vsnon ! volume per unit area of snow (m)
1035 :
1036 : real (kind=dbl_kind), dimension (nx_block,ny_block,ntrcr,ncat), intent(inout) :: &
1037 : trcrn ! tracers
1038 :
1039 : ! local variables
1040 :
1041 : integer (kind=int_kind) :: &
1042 : nt_qsno , & ! ! LCOV_EXCL_LINE
1043 : i, j, n , & ! standard indices ! LCOV_EXCL_LINE
1044 : it, kt , & ! tracer indices ! LCOV_EXCL_LINE
1045 : icells , & ! number of cells with ice ! LCOV_EXCL_LINE
1046 : ij
1047 :
1048 : real (kind=dbl_kind) :: &
1049 : rhos , & ! ! LCOV_EXCL_LINE
1050 5760 : Lfresh !
1051 :
1052 : integer (kind=int_kind), dimension (nx_block*ny_block) :: &
1053 45888 : indxi, indxj ! compressed indices
1054 :
1055 : character(len=*), parameter :: subname = '(tracers_to_state)'
1056 :
1057 22944 : call icepack_query_parameters(rhos_out=rhos, Lfresh_out=Lfresh)
1058 22944 : call icepack_query_tracer_indices(nt_qsno_out=nt_qsno)
1059 22944 : call icepack_warnings_flush(nu_diag)
1060 22944 : if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
1061 0 : file=__FILE__, line=__LINE__)
1062 :
1063 16186320 : aice0(:,:) = aim(:,:,0)
1064 :
1065 137664 : do n = 1, ncat
1066 :
1067 114720 : icells = 0
1068 3767880 : do j = 1, ny_block
1069 80931600 : do i = 1, nx_block
1070 80816880 : if (aim(i,j,n) > c0) then
1071 37664329 : icells = icells + 1
1072 37664329 : indxi(icells) = i
1073 37664329 : indxj(icells) = j
1074 : endif
1075 : enddo
1076 : enddo
1077 :
1078 : !-------------------------------------------------------------------
1079 : ! Compute state variables.
1080 : !-------------------------------------------------------------------
1081 :
1082 37779049 : do ij = 1, icells
1083 37664329 : i = indxi(ij)
1084 37664329 : j = indxj(ij)
1085 37664329 : aicen(i,j,n) = aim(i,j,n)
1086 37664329 : vicen(i,j,n) = aim(i,j,n)*trm(i,j,1,n) ! aice*hice
1087 37779049 : vsnon(i,j,n) = aim(i,j,n)*trm(i,j,2,n) ! aice*hsno
1088 : enddo ! ij
1089 114720 : kt = 2
1090 :
1091 2890944 : do it = 1, ntrcr
1092 2868000 : if (it >= nt_qsno .and. it < nt_qsno+nslyr) then
1093 37779049 : do ij = 1, icells
1094 37664329 : i = indxi(ij)
1095 37664329 : j = indxj(ij)
1096 37779049 : trcrn(i,j,it,n) = trm(i,j,kt+it,n) - rhos*Lfresh ! snow enthalpy
1097 : enddo
1098 : else
1099 868918127 : do ij = 1, icells
1100 866279567 : i = indxi(ij)
1101 866279567 : j = indxj(ij)
1102 868918127 : trcrn(i,j,it,n) = trm(i,j,kt+it,n) ! other tracers
1103 : enddo
1104 : endif
1105 : enddo
1106 : enddo ! ncat
1107 :
1108 22944 : end subroutine tracers_to_state
1109 :
1110 : !=======================================================================
1111 : !
1112 : ! Check whether values of conserved quantities have changed.
1113 : ! An error probably means that ghost cells are treated incorrectly.
1114 : !
1115 : ! author William H. Lipscomb, LANL
1116 :
1117 0 : subroutine global_conservation (ckflag, fieldid, &
1118 : asum_init, asum_final, & ! LCOV_EXCL_LINE
1119 0 : atsum_init, atsum_final)
1120 :
1121 : character(len=*), intent(in) :: &
1122 : fieldid ! field information string
1123 :
1124 : real (kind=dbl_kind), intent(in) :: &
1125 : asum_init , & ! initial global ice area ! LCOV_EXCL_LINE
1126 : asum_final ! final global ice area
1127 :
1128 : real (kind=dbl_kind), dimension(ntrace), intent(in), optional :: &
1129 : atsum_init, & ! initial global ice area*tracer ! LCOV_EXCL_LINE
1130 : atsum_final ! final global ice area*tracer
1131 :
1132 : logical (kind=log_kind), intent(inout) :: &
1133 : ckflag ! if true, abort on return
1134 :
1135 : ! local variables
1136 :
1137 : integer (kind=int_kind) :: &
1138 : nt ! tracer index
1139 :
1140 : real (kind=dbl_kind) :: &
1141 : puny , & ! ! LCOV_EXCL_LINE
1142 0 : diff ! difference between initial and final values
1143 :
1144 : character(len=*), parameter :: subname = '(global_conservation)'
1145 :
1146 0 : call icepack_query_parameters(puny_out=puny)
1147 0 : call icepack_warnings_flush(nu_diag)
1148 0 : if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
1149 0 : file=__FILE__, line=__LINE__)
1150 :
1151 0 : if (asum_init > puny) then
1152 0 : diff = asum_final - asum_init
1153 0 : if (abs(diff/asum_init) > puny) then
1154 0 : ckflag = .true.
1155 0 : write (nu_diag,*)
1156 0 : write (nu_diag,*) subname,'Ice area conserv error ', trim(fieldid)
1157 0 : write (nu_diag,*) subname,' Initial global area =', asum_init
1158 0 : write (nu_diag,*) subname,' Final global area =', asum_final
1159 0 : write (nu_diag,*) subname,' Fractional error =', abs(diff)/asum_init
1160 0 : write (nu_diag,*) subname,' asum_final-asum_init =', diff
1161 : endif
1162 : endif
1163 :
1164 0 : if (present(atsum_init)) then
1165 0 : do nt = 1, ntrace
1166 0 : if (abs(atsum_init(nt)) > puny) then
1167 0 : diff = atsum_final(nt) - atsum_init(nt)
1168 0 : if (abs(diff/atsum_init(nt)) > puny) then
1169 0 : ckflag = .true.
1170 0 : write (nu_diag,*)
1171 0 : write (nu_diag,*) subname,'Ice area*tracer conserv error ', trim(fieldid),nt
1172 0 : write (nu_diag,*) subname,' Tracer index =', nt
1173 0 : write (nu_diag,*) subname,' Initial global area*tracer =', atsum_init(nt)
1174 0 : write (nu_diag,*) subname,' Final global area*tracer =', atsum_final(nt)
1175 0 : write (nu_diag,*) subname,' Fractional error =', abs(diff)/atsum_init(nt)
1176 0 : write (nu_diag,*) subname,' atsum_final-atsum_init =', diff
1177 : endif
1178 : endif
1179 : enddo
1180 : endif ! present(atsum_init)
1181 :
1182 0 : end subroutine global_conservation
1183 :
1184 : !=======================================================================
1185 : !
1186 : ! At each grid point, compute the local max and min of a scalar
1187 : ! field phi: i.e., the max and min values in the nine-cell region
1188 : ! consisting of the home cell and its eight neighbors.
1189 : !
1190 : ! To extend to the neighbors of the neighbors (25 cells in all),
1191 : ! follow this call with a call to quasilocal_max_min.
1192 : !
1193 : ! author William H. Lipscomb, LANL
1194 :
1195 0 : subroutine local_max_min (nx_block, ny_block, &
1196 : ilo, ihi, jlo, jhi, & ! LCOV_EXCL_LINE
1197 : trm, & ! LCOV_EXCL_LINE
1198 : tmin, tmax, & ! LCOV_EXCL_LINE
1199 0 : aimask, trmask)
1200 :
1201 : integer (kind=int_kind), intent(in) :: &
1202 : nx_block, ny_block, & ! block dimensions ! LCOV_EXCL_LINE
1203 : ilo,ihi,jlo,jhi ! beginning and end of physical domain
1204 :
1205 : real (kind=dbl_kind), intent(in), dimension(nx_block,ny_block) :: &
1206 : aimask ! ice area mask
1207 :
1208 : real (kind=dbl_kind), intent(in), dimension (nx_block,ny_block,ntrace) :: &
1209 : trm , & ! tracer fields ! LCOV_EXCL_LINE
1210 : trmask ! tracer mask
1211 :
1212 : real (kind=dbl_kind), intent(out), dimension (nx_block,ny_block,ntrace) :: &
1213 : tmin , & ! local min tracer ! LCOV_EXCL_LINE
1214 : tmax ! local max tracer
1215 :
1216 : ! local variables
1217 :
1218 : integer (kind=int_kind) :: &
1219 : i, j , & ! horizontal indices ! LCOV_EXCL_LINE
1220 : nt, nt1 ! tracer indices
1221 :
1222 : real (kind=dbl_kind), dimension(nx_block,ny_block) :: &
1223 0 : phimask ! aimask or trmask, as appropriate
1224 :
1225 : real (kind=dbl_kind) :: &
1226 : phi_nw, phi_n, phi_ne , & ! field values in 8 neighbor cells ! LCOV_EXCL_LINE
1227 : phi_w , phi_e , & ! LCOV_EXCL_LINE
1228 0 : phi_sw, phi_s, phi_se
1229 :
1230 : character(len=*), parameter :: subname = '(local_max_min)'
1231 :
1232 0 : do nt = 1, ntrace
1233 :
1234 0 : if (tracer_type(nt)==1) then ! does not depend on another tracer
1235 :
1236 0 : do j = 1, ny_block
1237 0 : do i = 1, nx_block
1238 0 : phimask(i,j) = aimask(i,j)
1239 : enddo
1240 : enddo
1241 :
1242 : else ! depends on another tracer
1243 :
1244 0 : nt1 = depend(nt)
1245 0 : do j = 1, ny_block
1246 0 : do i = 1, nx_block
1247 0 : phimask(i,j) = trmask(i,j,nt1)
1248 : enddo
1249 : enddo
1250 :
1251 : endif
1252 :
1253 : !-----------------------------------------------------------------------
1254 : ! Store values of trm in the 8 neighbor cells.
1255 : ! If aimask = 1, use the true value; otherwise use the home cell value
1256 : ! so that non-physical values of phi do not contribute to the gradient.
1257 : !-----------------------------------------------------------------------
1258 :
1259 0 : do j = jlo, jhi
1260 0 : do i = ilo, ihi
1261 :
1262 0 : phi_nw = phimask(i-1,j+1) * trm(i-1,j+1,nt) &
1263 0 : + (c1-phimask(i-1,j+1))* trm(i, j, nt)
1264 0 : phi_n = phimask(i, j+1) * trm(i, j+1,nt) &
1265 0 : + (c1-phimask(i, j+1))* trm(i, j, nt)
1266 0 : phi_ne = phimask(i+1,j+1) * trm(i+1,j+1,nt) &
1267 0 : + (c1-phimask(i+1,j+1))* trm(i, j, nt)
1268 0 : phi_w = phimask(i-1,j) * trm(i-1,j, nt) &
1269 0 : + (c1-phimask(i-1,j)) * trm(i, j, nt)
1270 0 : phi_e = phimask(i+1,j) * trm(i+1,j, nt) &
1271 0 : + (c1-phimask(i+1,j)) * trm(i, j, nt)
1272 0 : phi_sw = phimask(i-1,j-1) * trm(i-1,j-1,nt) &
1273 0 : + (c1-phimask(i-1,j-1))* trm(i, j, nt)
1274 0 : phi_s = phimask(i, j-1) * trm(i, j-1,nt) &
1275 0 : + (c1-phimask(i, j-1))* trm(i, j, nt)
1276 0 : phi_se = phimask(i+1,j-1) * trm(i+1,j-1,nt) &
1277 0 : + (c1-phimask(i+1,j-1))* trm(i, j, nt)
1278 :
1279 : !-----------------------------------------------------------------------
1280 : ! Compute the minimum and maximum among the nine local cells.
1281 : !-----------------------------------------------------------------------
1282 :
1283 0 : tmax(i,j,nt) = max (phi_nw, phi_n, phi_ne, phi_w, &
1284 0 : trm(i,j,nt), phi_e, phi_sw, phi_s, phi_se)
1285 :
1286 0 : tmin(i,j,nt) = min (phi_nw, phi_n, phi_ne, phi_w, &
1287 0 : trm(i,j,nt), phi_e, phi_sw, phi_s, phi_se)
1288 :
1289 : enddo ! i
1290 : enddo ! j
1291 :
1292 : enddo ! nt
1293 :
1294 0 : end subroutine local_max_min
1295 :
1296 : !=======================================================================
1297 : !
1298 : ! Extend the local max and min by one grid cell in each direction.
1299 : ! Incremental remapping is monotone for the "quasilocal" max and min,
1300 : ! but in rare cases may violate monotonicity for the local max and min.
1301 : !
1302 : ! author William H. Lipscomb, LANL
1303 :
1304 0 : subroutine quasilocal_max_min (nx_block, ny_block, &
1305 : ilo, ihi, jlo, jhi, & ! LCOV_EXCL_LINE
1306 0 : tmin, tmax)
1307 :
1308 : integer (kind=int_kind), intent(in) :: &
1309 : nx_block, ny_block, & ! block dimensions ! LCOV_EXCL_LINE
1310 : ilo,ihi,jlo,jhi ! beginning and end of physical domain
1311 :
1312 : real (kind=dbl_kind), intent(inout), dimension (nx_block,ny_block,ntrace) :: &
1313 : tmin , & ! local min tracer ! LCOV_EXCL_LINE
1314 : tmax ! local max tracer
1315 :
1316 : ! local variables
1317 :
1318 : integer (kind=int_kind) :: &
1319 : i, j , & ! horizontal indices ! LCOV_EXCL_LINE
1320 : nt ! tracer index
1321 :
1322 : character(len=*), parameter :: subname = '(quasilocal_max_min)'
1323 :
1324 0 : do nt = 1, ntrace
1325 :
1326 0 : do j = jlo, jhi
1327 0 : do i = ilo, ihi
1328 :
1329 0 : tmax(i,j,nt) = &
1330 : max (tmax(i-1,j+1,nt), tmax(i,j+1,nt), tmax(i+1,j+1,nt), & ! LCOV_EXCL_LINE
1331 : tmax(i-1,j, nt), tmax(i,j, nt), tmax(i+1,j, nt), & ! LCOV_EXCL_LINE
1332 0 : tmax(i-1,j-1,nt), tmax(i,j-1,nt), tmax(i+1,j-1,nt))
1333 :
1334 0 : tmin(i,j,nt) = &
1335 : min (tmin(i-1,j+1,nt), tmin(i,j+1,nt), tmin(i+1,j+1,nt), & ! LCOV_EXCL_LINE
1336 : tmin(i-1,j, nt), tmin(i,j, nt), tmin(i+1,j, nt), & ! LCOV_EXCL_LINE
1337 0 : tmin(i-1,j-1,nt), tmin(i,j-1,nt), tmin(i+1,j-1,nt))
1338 :
1339 : enddo ! i
1340 : enddo ! j
1341 :
1342 : enddo
1343 :
1344 0 : end subroutine quasilocal_max_min
1345 :
1346 : !======================================================================
1347 : !
1348 : ! At each grid point, make sure that the new tracer values
1349 : ! fall between the local max and min values before transport.
1350 : !
1351 : ! author William H. Lipscomb, LANL
1352 :
1353 0 : subroutine check_monotonicity (nx_block, ny_block, &
1354 : ilo, ihi, jlo, jhi, & ! LCOV_EXCL_LINE
1355 : tmin, tmax, & ! LCOV_EXCL_LINE
1356 : aim, trm, & ! LCOV_EXCL_LINE
1357 : ckflag, & ! LCOV_EXCL_LINE
1358 : istop, jstop)
1359 :
1360 : integer (kind=int_kind), intent(in) :: &
1361 : nx_block, ny_block, & ! block dimensions ! LCOV_EXCL_LINE
1362 : ilo,ihi,jlo,jhi ! beginning and end of physical domain
1363 :
1364 : real (kind=dbl_kind), intent(in), dimension (nx_block,ny_block) :: &
1365 : aim ! new ice area
1366 :
1367 : real (kind=dbl_kind), intent(in), dimension (nx_block,ny_block,ntrace) :: &
1368 : trm ! new tracers
1369 :
1370 : real (kind=dbl_kind), intent(in), dimension (nx_block,ny_block,ntrace) :: &
1371 : tmin , & ! local min tracer ! LCOV_EXCL_LINE
1372 : tmax ! local max tracer
1373 :
1374 : logical (kind=log_kind), intent(inout) :: &
1375 : ckflag ! if true, abort on return
1376 :
1377 : integer (kind=int_kind), intent(inout) :: &
1378 : istop, jstop ! indices of grid cell where model aborts
1379 :
1380 : ! local variables
1381 :
1382 : integer (kind=int_kind) :: &
1383 : i, j , & ! horizontal indices ! LCOV_EXCL_LINE
1384 : nt, nt1, nt2 ! tracer indices
1385 :
1386 : real (kind=dbl_kind) :: &
1387 : puny , & ! ! LCOV_EXCL_LINE
1388 0 : w1, w2 ! work variables
1389 :
1390 : logical (kind=log_kind), dimension (nx_block, ny_block) :: &
1391 0 : l_check ! if true, check monotonicity
1392 :
1393 : character(len=*), parameter :: subname = '(check_monotonicity)'
1394 :
1395 0 : call icepack_query_parameters(puny_out=puny)
1396 0 : call icepack_warnings_flush(nu_diag)
1397 0 : if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
1398 0 : file=__FILE__, line=__LINE__)
1399 :
1400 0 : do nt = 1, ntrace
1401 :
1402 : !-------------------------------------------------------------------
1403 : ! Load logical array to identify tracers that need checking.
1404 : !-------------------------------------------------------------------
1405 :
1406 0 : if (tracer_type(nt)==1) then ! does not depend on another tracer
1407 :
1408 0 : do j = jlo, jhi
1409 0 : do i = ilo, ihi
1410 0 : if (aim(i,j) > puny) then
1411 0 : l_check(i,j) = .true.
1412 : else
1413 0 : l_check(i,j) = .false.
1414 : endif
1415 : enddo
1416 : enddo
1417 :
1418 0 : elseif (tracer_type(nt)==2) then ! depends on another tracer
1419 :
1420 0 : nt1 = depend(nt)
1421 0 : do j = jlo, jhi
1422 0 : do i = ilo, ihi
1423 0 : if (abs(trm(i,j,nt1)) > puny) then
1424 0 : l_check(i,j) = .true.
1425 : else
1426 0 : l_check(i,j) = .false.
1427 : endif
1428 : enddo
1429 : enddo
1430 :
1431 0 : elseif (tracer_type(nt)==3) then ! depends on two tracers
1432 :
1433 0 : nt1 = depend(nt)
1434 0 : nt2 = depend(nt1)
1435 0 : do j = jlo, jhi
1436 0 : do i = ilo, ihi
1437 0 : if (abs(trm(i,j,nt1)) > puny .and. &
1438 0 : abs(trm(i,j,nt2)) > puny) then
1439 0 : l_check(i,j) = .true.
1440 : else
1441 0 : l_check(i,j) = .false.
1442 : endif
1443 : enddo
1444 : enddo
1445 : endif
1446 :
1447 : !-------------------------------------------------------------------
1448 : ! Make sure new values lie between tmin and tmax
1449 : !-------------------------------------------------------------------
1450 :
1451 0 : do j = jlo, jhi
1452 0 : do i = ilo, ihi
1453 :
1454 0 : if (l_check(i,j)) then
1455 : ! w1 and w2 allow for roundoff error when abs(trm) is big
1456 0 : w1 = max(c1, abs(tmin(i,j,nt)))
1457 0 : w2 = max(c1, abs(tmax(i,j,nt)))
1458 0 : if (trm(i,j,nt) < tmin(i,j,nt)-w1*puny) then
1459 0 : ckflag = .true.
1460 0 : istop = i
1461 0 : jstop = j
1462 0 : write (nu_diag,*) ' '
1463 0 : write (nu_diag,*) 'new tracer < tmin'
1464 0 : write (nu_diag,*) 'i, j, nt =', i, j, nt
1465 0 : write (nu_diag,*) 'new tracer =', trm (i,j,nt)
1466 0 : write (nu_diag,*) 'tmin =' , tmin(i,j,nt)
1467 0 : write (nu_diag,*) 'ice area =' , aim(i,j)
1468 0 : elseif (trm(i,j,nt) > tmax(i,j,nt)+w2*puny) then
1469 0 : ckflag = .true.
1470 0 : istop = i
1471 0 : jstop = j
1472 0 : write (nu_diag,*) ' '
1473 0 : write (nu_diag,*) 'new tracer > tmax'
1474 0 : write (nu_diag,*) 'i, j, nt =', i, j, nt
1475 0 : write (nu_diag,*) 'new tracer =', trm (i,j,nt)
1476 0 : write (nu_diag,*) 'tmax =' , tmax(i,j,nt)
1477 0 : write (nu_diag,*) 'ice area =' , aim(i,j)
1478 : endif
1479 : endif
1480 :
1481 : enddo ! i
1482 : enddo ! j
1483 :
1484 : enddo ! nt
1485 :
1486 0 : end subroutine check_monotonicity
1487 :
1488 : !=======================================================================
1489 : ! The remaining subroutines are called by transport_upwind.
1490 : !=======================================================================
1491 : !
1492 : ! Fill work array with state variables in preparation for upwind transport
1493 :
1494 0 : subroutine state_to_work (nx_block, ny_block, &
1495 : ntrcr, & ! LCOV_EXCL_LINE
1496 : narr, trcr_depend, & ! LCOV_EXCL_LINE
1497 : aicen, trcrn, & ! LCOV_EXCL_LINE
1498 : vicen, vsnon, & ! LCOV_EXCL_LINE
1499 0 : aice0, works)
1500 :
1501 : use ice_domain_size, only: ncat
1502 :
1503 : integer (kind=int_kind), intent(in) :: &
1504 : nx_block, ny_block, & ! block dimensions ! LCOV_EXCL_LINE
1505 : ntrcr , & ! number of tracers in use ! LCOV_EXCL_LINE
1506 : narr ! number of 2D state variable arrays in works array
1507 :
1508 : integer (kind=int_kind), dimension (ntrcr), intent(in) :: &
1509 : trcr_depend ! = 0 for aicen tracers, 1 for vicen, 2 for vsnon
1510 :
1511 : real (kind=dbl_kind), dimension (nx_block,ny_block,ncat), intent(in) :: &
1512 : aicen , & ! concentration of ice ! LCOV_EXCL_LINE
1513 : vicen , & ! volume per unit area of ice (m) ! LCOV_EXCL_LINE
1514 : vsnon ! volume per unit area of snow (m)
1515 :
1516 : real (kind=dbl_kind), dimension (nx_block,ny_block,ntrcr,ncat), intent(in) :: &
1517 : trcrn ! ice tracers
1518 :
1519 : real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) :: &
1520 : aice0 ! concentration of open water
1521 :
1522 : real (kind=dbl_kind), dimension(nx_block,ny_block,narr), intent (out) :: &
1523 : works ! work array
1524 :
1525 : ! local variables
1526 :
1527 : integer (kind=int_kind) :: &
1528 : nt_alvl, nt_apnd, nt_fbri
1529 :
1530 : logical (kind=log_kind) :: &
1531 : tr_pond_lvl, tr_pond_topo
1532 :
1533 : integer (kind=int_kind) :: &
1534 : i, j, n, it, & ! counting indices ! LCOV_EXCL_LINE
1535 : narrays ! counter for number of state variable arrays
1536 :
1537 : character(len=*), parameter :: subname = '(state_to_work)'
1538 :
1539 : call icepack_query_tracer_flags(tr_pond_lvl_out=tr_pond_lvl, &
1540 0 : tr_pond_topo_out=tr_pond_topo)
1541 : call icepack_query_tracer_indices(nt_alvl_out=nt_alvl, nt_apnd_out=nt_apnd, &
1542 0 : nt_fbri_out=nt_fbri)
1543 0 : call icepack_warnings_flush(nu_diag)
1544 0 : if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
1545 0 : file=__FILE__, line=__LINE__)
1546 :
1547 : !-----------------------------------------------------------------
1548 : ! This array is used for performance (balance memory/cache vs
1549 : ! number of bound calls); a different number of arrays may perform
1550 : ! better depending on the machine used, number of processors, etc.
1551 : ! --tested on SGI R2000, using 4 pes for the ice model under MPI
1552 : !-----------------------------------------------------------------
1553 :
1554 0 : do j = 1, ny_block
1555 0 : do i = 1, nx_block
1556 0 : works(i,j,1) = aice0(i,j)
1557 : enddo
1558 : enddo
1559 0 : narrays = 1
1560 :
1561 0 : do n=1, ncat
1562 :
1563 0 : do j = 1, ny_block
1564 0 : do i = 1, nx_block
1565 0 : works(i,j,narrays+1) = aicen(i,j,n)
1566 0 : works(i,j,narrays+2) = vicen(i,j,n)
1567 0 : works(i,j,narrays+3) = vsnon(i,j,n)
1568 : enddo ! i
1569 : enddo ! j
1570 0 : narrays = narrays + 3
1571 :
1572 0 : do it = 1, ntrcr
1573 0 : if (trcr_depend(it) == 0) then
1574 0 : do j = 1, ny_block
1575 0 : do i = 1, nx_block
1576 0 : works(i,j,narrays+it) = aicen(i,j,n)*trcrn(i,j,it,n)
1577 : enddo
1578 : enddo
1579 0 : elseif (trcr_depend(it) == 1) then
1580 0 : do j = 1, ny_block
1581 0 : do i = 1, nx_block
1582 0 : works(i,j,narrays+it) = vicen(i,j,n)*trcrn(i,j,it,n)
1583 : enddo
1584 : enddo
1585 0 : elseif (trcr_depend(it) == 2) then
1586 0 : do j = 1, ny_block
1587 0 : do i = 1, nx_block
1588 0 : works(i,j,narrays+it) = vsnon(i,j,n)*trcrn(i,j,it,n)
1589 : enddo
1590 : enddo
1591 0 : elseif (trcr_depend(it) == 2+nt_alvl) then
1592 0 : do j = 1, ny_block
1593 0 : do i = 1, nx_block
1594 0 : works(i,j,narrays+it) = aicen(i,j ,n) &
1595 : * trcrn(i,j,nt_alvl,n) & ! LCOV_EXCL_LINE
1596 0 : * trcrn(i,j,it ,n)
1597 : enddo
1598 : enddo
1599 0 : elseif (trcr_depend(it) == 2+nt_apnd .and. &
1600 : tr_pond_topo) then
1601 0 : do j = 1, ny_block
1602 0 : do i = 1, nx_block
1603 0 : works(i,j,narrays+it) = aicen(i,j ,n) &
1604 : * trcrn(i,j,nt_apnd,n) & ! LCOV_EXCL_LINE
1605 0 : * trcrn(i,j,it ,n)
1606 : enddo
1607 : enddo
1608 0 : elseif (trcr_depend(it) == 2+nt_apnd .and. &
1609 : tr_pond_lvl) then
1610 0 : do j = 1, ny_block
1611 0 : do i = 1, nx_block
1612 0 : works(i,j,narrays+it) = aicen(i,j ,n) &
1613 : * trcrn(i,j,nt_alvl,n) & ! LCOV_EXCL_LINE
1614 : * trcrn(i,j,nt_apnd,n) & ! LCOV_EXCL_LINE
1615 0 : * trcrn(i,j,it ,n)
1616 : enddo
1617 : enddo
1618 0 : elseif (trcr_depend(it) == 2+nt_fbri) then
1619 0 : do j = 1, ny_block
1620 0 : do i = 1, nx_block
1621 0 : works(i,j,narrays+it) = vicen(i,j ,n) &
1622 : * trcrn(i,j,nt_fbri,n) & ! LCOV_EXCL_LINE
1623 0 : * trcrn(i,j,it ,n)
1624 : enddo
1625 : enddo
1626 : endif
1627 : enddo
1628 0 : narrays = narrays + ntrcr
1629 :
1630 : enddo ! n
1631 :
1632 0 : if (narr /= narrays) write(nu_diag,*) &
1633 0 : "Wrong number of arrays in transport bound call"
1634 :
1635 0 : end subroutine state_to_work
1636 :
1637 : !=======================================================================
1638 : !
1639 : ! Convert work array back to state variables
1640 :
1641 0 : subroutine work_to_state (nx_block, ny_block, &
1642 : ntrcr, narr, & ! LCOV_EXCL_LINE
1643 : trcr_depend, & ! LCOV_EXCL_LINE
1644 : trcr_base, & ! LCOV_EXCL_LINE
1645 : n_trcr_strata, & ! LCOV_EXCL_LINE
1646 : nt_strata, & ! LCOV_EXCL_LINE
1647 : tmask, Tf, & ! LCOV_EXCL_LINE
1648 : aicen, trcrn, & ! LCOV_EXCL_LINE
1649 : vicen, vsnon, & ! LCOV_EXCL_LINE
1650 0 : aice0, works)
1651 :
1652 : use ice_domain_size, only: ncat
1653 :
1654 : integer (kind=int_kind), intent (in) :: &
1655 : nx_block, ny_block, & ! block dimensions ! LCOV_EXCL_LINE
1656 : ntrcr , & ! number of tracers in use ! LCOV_EXCL_LINE
1657 : narr ! number of 2D state variable arrays in works array
1658 :
1659 : integer (kind=int_kind), dimension (ntrcr), intent(in) :: &
1660 : trcr_depend, & ! = 0 for aicen tracers, 1 for vicen, 2 for vsnon ! LCOV_EXCL_LINE
1661 : n_trcr_strata ! number of underlying tracer layers
1662 :
1663 : real (kind=dbl_kind), dimension (ntrcr,3), intent(in) :: &
1664 : trcr_base ! = 0 or 1 depending on tracer dependency
1665 : ! argument 2: (1) aice, (2) vice, (3) vsno
1666 :
1667 : integer (kind=int_kind), dimension (ntrcr,2), intent(in) :: &
1668 : nt_strata ! indices of underlying tracer layers
1669 :
1670 : logical (kind=log_kind), intent (in) :: &
1671 : tmask (nx_block,ny_block)
1672 :
1673 : real (kind=dbl_kind), intent (in) :: &
1674 : Tf (nx_block,ny_block), & ! LCOV_EXCL_LINE
1675 : works (nx_block,ny_block,narr)
1676 :
1677 : real (kind=dbl_kind), dimension (nx_block,ny_block,ncat), intent(out) :: &
1678 : aicen , & ! concentration of ice ! LCOV_EXCL_LINE
1679 : vicen , & ! volume per unit area of ice (m) ! LCOV_EXCL_LINE
1680 : vsnon ! volume per unit area of snow (m)
1681 :
1682 : real (kind=dbl_kind), dimension (nx_block,ny_block,ntrcr,ncat),intent(out) :: &
1683 : trcrn ! ice tracers
1684 :
1685 : real (kind=dbl_kind), dimension (nx_block,ny_block), intent(out) :: &
1686 : aice0 ! concentration of open water
1687 :
1688 : ! local variables
1689 :
1690 : integer (kind=int_kind) :: &
1691 : i, j, ij, n, & ! counting indices ! LCOV_EXCL_LINE
1692 : narrays , & ! counter for number of state variable arrays ! LCOV_EXCL_LINE
1693 : nt_Tsfc , & ! Tsfc tracer number ! LCOV_EXCL_LINE
1694 : icells ! number of ocean/ice cells
1695 :
1696 : integer (kind=int_kind), dimension (nx_block*ny_block) :: &
1697 0 : indxi, indxj
1698 :
1699 : real (kind=dbl_kind), dimension (nx_block*ny_block,narr) :: &
1700 0 : work
1701 :
1702 : character(len=*), parameter :: subname = '(work_to_state)'
1703 :
1704 0 : call icepack_query_tracer_indices(nt_Tsfc_out=nt_Tsfc)
1705 0 : call icepack_warnings_flush(nu_diag)
1706 0 : if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
1707 0 : file=__FILE__, line=__LINE__)
1708 :
1709 : ! for call to compute_tracers
1710 0 : icells = 0
1711 0 : do j = 1, ny_block
1712 0 : do i = 1, nx_block
1713 0 : icells = icells + 1
1714 0 : indxi(icells) = i
1715 0 : indxj(icells) = j
1716 0 : work (icells,:) = works(i,j,:)
1717 : enddo
1718 : enddo
1719 :
1720 0 : do j=1,ny_block
1721 0 : do i=1,nx_block
1722 0 : aice0(i,j) = works(i,j,1)
1723 : enddo
1724 : enddo
1725 0 : narrays = 1 ! aice0 is first array
1726 :
1727 0 : do n=1,ncat
1728 :
1729 0 : do j=1,ny_block
1730 0 : do i=1,nx_block
1731 0 : aicen(i,j,n) = works(i,j,narrays+1)
1732 0 : vicen(i,j,n) = works(i,j,narrays+2)
1733 0 : vsnon(i,j,n) = works(i,j,narrays+3)
1734 : enddo
1735 : enddo
1736 0 : narrays = narrays + 3
1737 :
1738 0 : do ij = 1, icells
1739 0 : i = indxi(ij)
1740 0 : j = indxj(ij)
1741 :
1742 : call icepack_compute_tracers(ntrcr = ntrcr, &
1743 : trcr_depend = trcr_depend(:), & ! LCOV_EXCL_LINE
1744 : atrcrn = work (ij,narrays+1:narrays+ntrcr), & ! LCOV_EXCL_LINE
1745 : aicen = aicen(i,j,n), & ! LCOV_EXCL_LINE
1746 : vicen = vicen(i,j,n), & ! LCOV_EXCL_LINE
1747 : vsnon = vsnon(i,j,n), & ! LCOV_EXCL_LINE
1748 : trcr_base = trcr_base(:,:), & ! LCOV_EXCL_LINE
1749 : n_trcr_strata = n_trcr_strata(:), & ! LCOV_EXCL_LINE
1750 : nt_strata = nt_strata(:,:), & ! LCOV_EXCL_LINE
1751 : trcrn = trcrn(i,j,:,n), & ! LCOV_EXCL_LINE
1752 0 : Tf = Tf(i,j))
1753 :
1754 : ! tcraig, don't let land points get non-zero Tsfc
1755 0 : if (.not.tmask(i,j)) then
1756 0 : trcrn(i,j,nt_Tsfc,n) = c0
1757 : endif
1758 :
1759 : enddo
1760 :
1761 0 : narrays = narrays + ntrcr
1762 :
1763 : enddo ! ncat
1764 :
1765 0 : call icepack_warnings_flush(nu_diag)
1766 0 : if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
1767 0 : file=__FILE__, line=__LINE__)
1768 :
1769 0 : end subroutine work_to_state
1770 :
1771 : !=======================================================================
1772 : !
1773 : ! upwind transport algorithm
1774 :
1775 0 : subroutine upwind_field (nx_block, ny_block, &
1776 : ilo, ihi, jlo, jhi, & ! LCOV_EXCL_LINE
1777 : dt, & ! LCOV_EXCL_LINE
1778 : narrays, phi, & ! LCOV_EXCL_LINE
1779 : uee, vnn, & ! LCOV_EXCL_LINE
1780 : HTE, HTN, & ! LCOV_EXCL_LINE
1781 0 : tarea)
1782 :
1783 : integer (kind=int_kind), intent (in) :: &
1784 : nx_block, ny_block, & ! block dimensions ! LCOV_EXCL_LINE
1785 : ilo,ihi,jlo,jhi , & ! beginning and end of physical domain ! LCOV_EXCL_LINE
1786 : narrays ! number of 2D arrays to be transported
1787 :
1788 : real (kind=dbl_kind), intent(in) :: &
1789 : dt ! time step
1790 :
1791 : real (kind=dbl_kind), dimension(nx_block,ny_block,narrays), intent(inout) :: &
1792 : phi ! scalar field
1793 :
1794 : real (kind=dbl_kind), dimension(nx_block,ny_block), intent(in) :: &
1795 : uee, vnn ! cell edge velocities
1796 :
1797 : real (kind=dbl_kind), dimension(nx_block,ny_block), intent(in) :: &
1798 : HTE , & ! length of east cell edge ! LCOV_EXCL_LINE
1799 : HTN , & ! length of north cell edge ! LCOV_EXCL_LINE
1800 : tarea ! grid cell area
1801 :
1802 : ! local variables
1803 :
1804 : integer (kind=int_kind) :: &
1805 : i, j, n ! standard indices
1806 :
1807 : real (kind=dbl_kind), dimension (nx_block,ny_block) :: &
1808 0 : worka, workb
1809 :
1810 : character(len=*), parameter :: subname = '(upwind_field)'
1811 :
1812 : !-------------------------------------------------------------------
1813 : ! upwind transport
1814 : !-------------------------------------------------------------------
1815 :
1816 0 : do n = 1, narrays
1817 :
1818 0 : do j = 1, jhi
1819 0 : do i = 1, ihi
1820 0 : worka(i,j)= &
1821 0 : upwind(phi(i,j,n),phi(i+1,j ,n),uee(i,j),HTE(i,j),dt)
1822 0 : workb(i,j)= &
1823 0 : upwind(phi(i,j,n),phi(i ,j+1,n),vnn(i,j),HTN(i,j),dt)
1824 : enddo
1825 : enddo
1826 :
1827 0 : do j = jlo, jhi
1828 0 : do i = ilo, ihi
1829 0 : phi(i,j,n) = phi(i,j,n) - ( worka(i,j)-worka(i-1,j ) &
1830 : + workb(i,j)-workb(i ,j-1) ) & ! LCOV_EXCL_LINE
1831 0 : / tarea(i,j)
1832 : enddo
1833 : enddo
1834 :
1835 : enddo ! narrays
1836 :
1837 0 : end subroutine upwind_field
1838 :
1839 : !=======================================================================
1840 : !
1841 : ! Define upwind function
1842 : !
1843 :
1844 0 : real(kind=dbl_kind) function upwind(y1,y2,a,h,dt)
1845 :
1846 : real(kind=dbl_kind), intent(in) :: y1,y2,a,h,dt
1847 :
1848 0 : upwind = p5*dt*h*((a+abs(a))*y1+(a-abs(a))*y2)
1849 :
1850 0 : end function upwind
1851 :
1852 : !=======================================================================
1853 :
1854 : end module ice_transport_driver
1855 :
1856 : !=======================================================================
|
