00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022 module ice_transport_driver
00023
00024
00025 use ice_kinds_mod
00026 use ice_communicate, only: my_task, master_task
00027 use ice_domain_size
00028 use ice_constants
00029 use ice_fileunits
00030
00031
00032
00033 implicit none
00034 save
00035
00036 character (len=char_len) ::
00037 advection
00038
00039
00040
00041 logical, parameter ::
00042 l_fixed_area = .false.
00043
00044
00045
00046
00047 integer (kind=int_kind), parameter ::
00048 ntrace = 2+ntrcr+nilyr+nslyr
00049
00050 integer (kind=int_kind), dimension (ntrace) ::
00051 tracer_type ,
00052 depend
00053
00054 logical (kind=log_kind), dimension (ntrace) ::
00055 has_dependents
00056
00057 integer (kind=int_kind), parameter ::
00058 integral_order = 3
00059
00060
00061 logical (kind=log_kind), parameter ::
00062 l_dp_midpt = .true.
00063
00064
00065
00066
00067 contains
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077 subroutine init_transport
00078
00079
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091 use ice_state, only: trcr_depend
00092 use ice_exit
00093 use ice_timers
00094 use ice_transport_remap, only: init_remap
00095
00096
00097
00098 integer (kind=int_kind) ::
00099 k, nt, nt1
00100
00101 call ice_timer_start(timer_advect)
00102
00103 if (trim(advection)=='remap') then
00104
00105
00106
00107
00108
00109
00110 depend(1:2) = 0
00111 tracer_type(1:2) = 1
00112
00113 k = 2
00114
00115 do nt = 1, ntrcr
00116 depend(k+nt) = trcr_depend(nt)
00117
00118
00119 if (trcr_depend(nt) == 0) then
00120 tracer_type(k+nt) = 1
00121 else
00122 tracer_type(k+nt) = 2
00123 endif
00124 enddo
00125
00126 k = k + ntrcr
00127
00128 depend(k+1:k+nilyr) = 1
00129 tracer_type(k+1:k+nilyr) = 2
00130
00131 k = k + nilyr
00132
00133 depend(k+1:k+nslyr) = 2
00134 tracer_type(k+1:k+nslyr) = 2
00135
00136 has_dependents = .false.
00137 do nt = 1, ntrace
00138 if (depend(nt) > 0) then
00139 nt1 = depend(nt)
00140 has_dependents(nt1) = .true.
00141 if (nt1 > nt) then
00142 write(nu_diag,*) &
00143 'Tracer nt2 =',nt,' depends on tracer nt1 =',nt1
00144 call abort_ice &
00145 ('ice: remap transport: Must have nt2 > nt1')
00146 endif
00147 endif
00148 enddo
00149
00150 call init_remap
00151
00152 else
00153
00154 continue
00155
00156 endif
00157
00158 call ice_timer_stop(timer_advect)
00159
00160 end subroutine init_transport
00161
00162
00163
00164
00165
00166
00167
00168
00169 subroutine transport_remap (dt)
00170
00171
00172
00173
00174
00175
00176
00177
00178
00179
00180
00181
00182
00183
00184
00185
00186
00187
00188 use ice_boundary
00189 use ice_global_reductions
00190 use ice_domain
00191 use ice_blocks
00192 use ice_state
00193 use ice_grid, only: tarea, HTE, HTN
00194 use ice_exit
00195 use ice_calendar, only: istep1, istep, diagfreq
00196 use ice_timers
00197 use ice_transport_remap, only: horizontal_remap, make_masks
00198
00199
00200
00201 real (kind=dbl_kind), intent(in) ::
00202 dt
00203
00204
00205
00206
00207
00208 integer (kind=int_kind) ::
00209 i, j ,
00210 iblk ,
00211 ilo,ihi,jlo,jhi,
00212 n ,
00213 nt, nt1, nt2
00214
00215 real (kind=dbl_kind),
00216 dimension (nx_block,ny_block,0:ncat,max_blocks) ::
00217 aim ,
00218 aimask
00219
00220 real (kind=dbl_kind),
00221 dimension (nx_block,ny_block,ntrace,ncat,max_blocks) ::
00222 trm ,
00223 trmask
00224
00225 logical (kind=log_kind) ::
00226 l_stop
00227
00228 integer (kind=int_kind) ::
00229 istop, jstop
00230
00231 integer (kind=int_kind), dimension(0:ncat,max_blocks) ::
00232 icellsnc
00233
00234 integer (kind=int_kind),
00235 dimension(nx_block*ny_block,0:ncat,max_blocks) ::
00236 indxinc, indxjnc
00237
00238 type (block) ::
00239 this_block
00240
00241
00242
00243
00244
00245
00246
00247
00248
00249 real (kind=dbl_kind), dimension(nx_block,ny_block,max_blocks) ::
00250 edgearea_e ,
00251 edgearea_n
00252
00253
00254
00255 logical (kind=log_kind), parameter ::
00256 l_conservation_check = .true. ,
00257 l_monotonicity_check = .true.
00258
00259 real (kind=dbl_kind), dimension(0:ncat) ::
00260 asum_init ,
00261 asum_final
00262
00263 real (kind=dbl_kind), dimension(ntrace,ncat) ::
00264 atsum_init ,
00265 atsum_final
00266
00267 real (kind=dbl_kind), dimension (:,:,:,:,:), allocatable ::
00268 tmin ,
00269 tmax
00270
00271 integer (kind=int_kind) :: alloc_error
00272
00273 real (kind=dbl_kind), dimension(nx_block,ny_block,max_blocks) ::
00274 work1
00275
00276 call ice_timer_start(timer_advect)
00277
00278
00279
00280
00281
00282
00283 l_stop = .false.
00284 istop = 0
00285 jstop = 0
00286
00287
00288
00289
00290
00291
00292
00293
00294
00295
00296
00297
00298
00299
00300
00301
00302
00303
00304
00305
00306
00307
00308
00309
00310
00311
00312
00313
00314
00315
00316
00317
00318
00319
00320
00321
00322
00323
00324
00325
00326
00327
00328
00329
00330
00331
00332 do iblk = 1, nblocks
00333
00334
00335
00336
00337
00338 call state_to_tracers(nx_block, ny_block, &
00339 aice0(:,:, iblk), &
00340 aicen(:,:,:,iblk), trcrn(:,:,:,:,iblk), &
00341 vicen(:,:,:,iblk), vsnon(:,:, :,iblk), &
00342 eicen(:,:,:,iblk), esnon(:,:, :,iblk), &
00343 aim (:,:,:,iblk), trm (:,:,:,:,iblk))
00344
00345 enddo
00346
00347
00348
00349
00350
00351
00352 if (l_conservation_check .and. mod(istep,diagfreq) == 0) then
00353
00354
00355
00356
00357
00358 do n = 0, ncat
00359 asum_init(n) = global_sum(aim(:,:,n,:), distrb_info, &
00360 field_loc_center, tarea)
00361 enddo
00362
00363 do n = 1, ncat
00364 do nt = 1, ntrace
00365 if (tracer_type(nt)==1) then
00366 atsum_init(nt,n) = &
00367 global_sum_prod(trm(:,:,nt,n,:), aim(:,:,n,:), &
00368 distrb_info, field_loc_center, &
00369 tarea)
00370 elseif (tracer_type(nt)==2) then
00371 nt1 = depend(nt)
00372 do iblk = 1, nblocks
00373 do j= 1,ny_block
00374 do i = 1,nx_block
00375 work1(i,j,iblk) = trm(i,j,nt,n,iblk)*trm(i,j,nt1,n,iblk)
00376 end do
00377 end do
00378 end do
00379 atsum_init(nt,n) = &
00380 global_sum_prod(work1(:,:,:), aim(:,:,n,:), &
00381 distrb_info, field_loc_center, &
00382 tarea)
00383 elseif (tracer_type(nt)==3) then
00384 nt1 = depend(nt)
00385 nt2 = depend(nt1)
00386 do iblk = 1, nblocks
00387 do j= 1,ny_block
00388 do i = 1,nx_block
00389 work1(i,j,iblk) = trm(i,j,nt,n,iblk)*trm(i,j,nt1,n,iblk) &
00390 *trm(i,j,nt2,n,iblk)
00391 end do
00392 end do
00393 end do
00394 atsum_init(nt,n) = &
00395 global_sum_prod(work1(:,:,:), aim(:,:,n,:), &
00396 distrb_info, field_loc_center, &
00397 tarea)
00398 endif
00399 enddo
00400 enddo
00401
00402 endif
00403
00404 if (l_monotonicity_check .and. mod(istep,diagfreq) == 0) then
00405
00406 allocate(tmin(nx_block,ny_block,ntrace,ncat,max_blocks), &
00407 tmax(nx_block,ny_block,ntrace,ncat,max_blocks), &
00408 STAT=alloc_error)
00409
00410 if (alloc_error /= 0) &
00411 call abort_ice ('ice: allocation error')
00412
00413 tmin(:,:,:,:,:) = c0
00414 tmax(:,:,:,:,:) = c0
00415
00416
00417 do iblk = 1, nblocks
00418 this_block = get_block(blocks_ice(iblk),iblk)
00419 ilo = this_block%ilo
00420 ihi = this_block%ihi
00421 jlo = this_block%jlo
00422 jhi = this_block%jhi
00423
00424
00425
00426
00427
00428
00429
00430 call make_masks (nx_block, ny_block, &
00431 ilo, ihi, jlo, jhi, &
00432 nghost, has_dependents, &
00433 icellsnc(:,iblk), &
00434 indxinc(:,:,iblk), indxjnc(:,:,iblk), &
00435 aim(:,:,:,iblk), aimask(:,:,:,iblk), &
00436 trm(:,:,:,:,iblk), trmask(:,:,:,:,iblk))
00437
00438
00439
00440
00441
00442 do n = 1, ncat
00443 call local_max_min &
00444 (nx_block, ny_block, &
00445 ilo, ihi, jlo, jhi, &
00446 trm (:,:,:,n,iblk), &
00447 tmin(:,:,:,n,iblk), tmax (:,:,:,n,iblk), &
00448 aimask(:,:,n,iblk), trmask(:,:,:,n,iblk))
00449 enddo
00450 enddo
00451
00452
00453 call ice_timer_start(timer_bound)
00454 call ice_HaloUpdate (tmin, halo_info, &
00455 field_loc_center, field_type_scalar)
00456 call ice_HaloUpdate (tmax, halo_info, &
00457 field_loc_center, field_type_scalar)
00458 call ice_timer_stop(timer_bound)
00459
00460
00461 do iblk = 1, nblocks
00462 this_block = get_block(blocks_ice(iblk),iblk)
00463 ilo = this_block%ilo
00464 ihi = this_block%ihi
00465 jlo = this_block%jlo
00466 jhi = this_block%jhi
00467
00468 do n = 1, ncat
00469 call quasilocal_max_min (nx_block, ny_block, &
00470 ilo, ihi, jlo, jhi, &
00471 tmin(:,:,:,n,iblk), &
00472 tmax(:,:,:,n,iblk))
00473 enddo
00474 enddo
00475
00476
00477 endif
00478
00479
00480
00481
00482
00483
00484
00485
00486
00487 do iblk = 1, nblocks
00488 do j = 1, ny_block
00489 do i = 1, nx_block
00490 edgearea_e(i,j,iblk) = c0
00491 edgearea_n(i,j,iblk) = c0
00492 enddo
00493 enddo
00494 enddo
00495
00496
00497 if (l_fixed_area) then
00498
00499
00500 do iblk = 1, nblocks
00501 this_block = get_block(blocks_ice(iblk),iblk)
00502 ilo = this_block%ilo
00503 ihi = this_block%ihi
00504 jlo = this_block%jlo
00505 jhi = this_block%jhi
00506
00507 do j = jlo, jhi
00508 do i = ilo-1, ihi
00509 edgearea_e(i,j,iblk) = (uvel(i,j,iblk) + uvel(i,j-1,iblk)) &
00510 * p5 * HTE(i,j,iblk) * dt
00511 enddo
00512 enddo
00513
00514 do j = jlo-1, jhi
00515 do i = ilo, ihi
00516 edgearea_n(i,j,iblk) = (vvel(i,j,iblk) + vvel(i-1,j,iblk)) &
00517 * p5 * HTN(i,j,iblk) * dt
00518 enddo
00519 enddo
00520
00521 enddo
00522
00523
00524 endif
00525
00526 call horizontal_remap (dt, &
00527 uvel (:,:,:), vvel (:,:,:), &
00528 aim (:,:,:,:), trm (:,:,:,:,:), &
00529 l_fixed_area, &
00530 edgearea_e(:,:,:), edgearea_n(:,:,:), &
00531 tracer_type, depend, &
00532 has_dependents, integral_order, &
00533 l_dp_midpt)
00534
00535
00536
00537
00538
00539
00540 do iblk = 1, nblocks
00541
00542 call tracers_to_state (nx_block, ny_block, &
00543 aim (:,:,:,iblk), trm (:,:,:,:,iblk), &
00544 aice0(:,:, iblk), &
00545 aicen(:,:,:,iblk), trcrn(:,:,:,:,iblk), &
00546 vicen(:,:,:,iblk), vsnon(:,:, :,iblk), &
00547 eicen(:,:,:,iblk), esnon(:,:, :,iblk))
00548
00549 enddo
00550
00551
00552
00553
00554
00555
00556 call ice_timer_start(timer_bound)
00557
00558 call bound_state (aicen, trcrn, &
00559 vicen, vsnon, &
00560 eicen, esnon)
00561
00562 call ice_timer_stop(timer_bound)
00563
00564
00565
00566
00567
00568
00569
00570
00571
00572
00573 if (l_conservation_check .and. mod(istep,diagfreq) == 0) then
00574
00575 do n = 0, ncat
00576 asum_final(n) = global_sum(aim(:,:,n,:), distrb_info, &
00577 field_loc_center, tarea)
00578 enddo
00579
00580 do n = 1, ncat
00581 do nt = 1, ntrace
00582 if (tracer_type(nt)==1) then
00583 atsum_final(nt,n) = &
00584 global_sum_prod(trm(:,:,nt,n,:), aim(:,:,n,:), &
00585 distrb_info, field_loc_center, &
00586 tarea)
00587 elseif (tracer_type(nt)==2) then
00588 nt1 = depend(nt)
00589 do iblk = 1, nblocks
00590 do j= 1,ny_block
00591 do i = 1,nx_block
00592 work1(i,j,iblk) = trm(i,j,nt,n,iblk)*trm(i,j,nt1,n,iblk)
00593 end do
00594 end do
00595 end do
00596 atsum_final(nt,n) = &
00597 global_sum_prod(work1(:,:,:), aim(:,:,n,:), &
00598 distrb_info, field_loc_center, &
00599 tarea)
00600 elseif (tracer_type(nt)==3) then
00601 nt1 = depend(nt)
00602 nt2 = depend(nt1)
00603 do iblk = 1, nblocks
00604 do j= 1,ny_block
00605 do i = 1,nx_block
00606 work1(i,j,iblk) = trm(i,j,nt,n,iblk)*trm(i,j,nt1,n,iblk) &
00607 *trm(i,j,nt2,n,iblk)
00608 end do
00609 end do
00610 end do
00611 atsum_final(nt,n) = &
00612 global_sum_prod(work1(:,:,:), aim(:,:,n,:), &
00613 distrb_info, field_loc_center, &
00614 tarea)
00615 endif
00616 enddo
00617 enddo
00618
00619
00620 if (my_task == master_task) then
00621 call global_conservation (l_stop, &
00622 asum_init(0), asum_final(0))
00623
00624 if (l_stop) then
00625 write (nu_diag,*) 'istep1 =', istep1
00626 write (nu_diag,*) 'transport: conservation error, cat 0'
00627 l_stop = .false.
00628 call abort_ice('ice remap transport: conservation error')
00629 endif
00630
00631 do n = 1, ncat
00632 call global_conservation &
00633 (l_stop, &
00634 asum_init(n), asum_final(n), &
00635 atsum_init(:,n), atsum_final(:,n))
00636
00637 if (l_stop) then
00638 write (nu_diag,*) 'istep1, cat =', &
00639 istep1, n
00640 write (nu_diag,*) 'transport: conservation error, cat ',n
00641 l_stop = .false.
00642 call abort_ice &
00643 ('ice remap transport: conservation error')
00644 endif
00645 enddo
00646
00647 endif
00648
00649 endif
00650
00651
00652
00653
00654
00655 if (l_monotonicity_check .and. mod(istep,diagfreq) == 0) then
00656 do iblk = 1, nblocks
00657 this_block = get_block(blocks_ice(iblk),iblk)
00658 ilo = this_block%ilo
00659 ihi = this_block%ihi
00660 jlo = this_block%jlo
00661 jhi = this_block%jhi
00662
00663 do n = 1, ncat
00664 call check_monotonicity &
00665 (nx_block, ny_block, &
00666 ilo, ihi, jlo, jhi, &
00667 iblk, &
00668 tmin(:,:,:,n,iblk), tmax(:,:,:,n,iblk), &
00669 aim (:,:, n,iblk), trm (:,:,:,n,iblk), &
00670 l_stop, &
00671 istop, jstop)
00672
00673 if (l_stop) then
00674 write (nu_diag,*) 'istep1, my_task, iblk, cat =', &
00675 istep1, my_task, iblk, n
00676 write (nu_diag,*) 'i_glob, j_glob',this_block%i_glob(istop), &
00677 this_block%j_glob(jstop)
00678 call abort_ice('ice remap transport: monotonicity error')
00679 endif
00680
00681 enddo
00682
00683 enddo
00684
00685 deallocate(tmin, tmax, STAT=alloc_error)
00686 if (alloc_error /= 0) call abort_ice ('deallocation error')
00687
00688 endif
00689
00690 call ice_timer_stop(timer_advect)
00691
00692 end subroutine transport_remap
00693
00694
00695
00696
00697
00698
00699
00700
00701 subroutine transport_upwind (dt)
00702
00703
00704
00705
00706
00707
00708
00709
00710
00711
00712
00713
00714 use ice_boundary
00715 use ice_blocks
00716 use ice_domain
00717 use ice_state
00718 use ice_grid, only: HTE, HTN, tarea
00719 use ice_timers
00720
00721
00722
00723 real (kind=dbl_kind), intent(in) ::
00724 dt
00725
00726
00727
00728 integer (kind=int_kind), parameter ::
00729 narr = 1 + ncat*(3+ntrcr)
00730
00731
00732 integer (kind=int_kind) ::
00733 i, j, iblk ,
00734 ilo,ihi,jlo,jhi
00735
00736 real (kind=dbl_kind), dimension (nx_block,ny_block,nblocks) ::
00737 uee, vnn
00738
00739 real (kind=dbl_kind),
00740 dimension (nx_block,ny_block,narr,max_blocks) ::
00741 works
00742
00743 type (block) ::
00744 this_block
00745
00746 call ice_timer_start(timer_advect)
00747
00748
00749
00750
00751
00752
00753
00754
00755
00756 uee(:,:,:) = c0
00757 vnn(:,:,:) = c0
00758
00759
00760
00761
00762
00763
00764 do iblk = 1, nblocks
00765 this_block = get_block(blocks_ice(iblk),iblk)
00766 ilo = this_block%ilo
00767 ihi = this_block%ihi
00768 jlo = this_block%jlo
00769 jhi = this_block%jhi
00770
00771 do j = jlo, jhi
00772 do i = ilo, ihi
00773 uee(i,j,iblk) = p5*(uvel(i,j,iblk) + uvel(i,j-1,iblk))
00774 vnn(i,j,iblk) = p5*(vvel(i,j,iblk) + vvel(i-1,j,iblk))
00775 enddo
00776 enddo
00777 enddo
00778
00779
00780 call ice_timer_start(timer_bound)
00781 call ice_HaloUpdate (uee, halo_info, &
00782 field_loc_Eface, field_type_vector)
00783 call ice_HaloUpdate (vnn, halo_info, &
00784 field_loc_Nface, field_type_vector)
00785 call ice_timer_stop(timer_bound)
00786
00787
00788 do iblk = 1, nblocks
00789 this_block = get_block(blocks_ice(iblk),iblk)
00790 ilo = this_block%ilo
00791 ihi = this_block%ihi
00792 jlo = this_block%jlo
00793 jhi = this_block%jhi
00794
00795
00796
00797
00798
00799
00800 call state_to_work (nx_block, ny_block, &
00801 narr, trcr_depend, &
00802 aicen (:,:, :,iblk), trcrn (:,:,:,:,iblk), &
00803 vicen (:,:, :,iblk), vsnon (:,:, :,iblk), &
00804 aice0 (:,:, iblk), works (:,:, :,iblk))
00805
00806
00807
00808
00809
00810 call upwind_field (nx_block, ny_block, &
00811 ilo, ihi, jlo, jhi, &
00812 dt, &
00813 narr, works(:,:,:,iblk), &
00814 uee(:,:,iblk), vnn (:,:,iblk), &
00815 HTE(:,:,iblk), HTN (:,:,iblk), &
00816 tarea(:,:,iblk))
00817
00818 call upwind_field (nx_block, ny_block, &
00819 ilo, ihi, jlo, jhi, &
00820 dt, &
00821 ntilyr, eicen(:,:,:,iblk), &
00822 uee(:,:,iblk), vnn (:,:,iblk), &
00823 HTE(:,:,iblk), HTN (:,:,iblk), &
00824 tarea(:,:,iblk))
00825
00826 call upwind_field (nx_block, ny_block, &
00827 ilo, ihi, jlo, jhi, &
00828 dt, &
00829 ntslyr, esnon(:,:,:,iblk), &
00830 uee(:,:,iblk), vnn (:,:,iblk), &
00831 HTE(:,:,iblk), HTN (:,:,iblk), &
00832 tarea(:,:,iblk))
00833
00834
00835
00836
00837
00838 call work_to_state (nx_block, ny_block, &
00839 narr, trcr_depend, &
00840 aicen(:,:, :,iblk), trcrn (:,:,:,:,iblk), &
00841 vicen(:,:, :,iblk), vsnon (:,:, :,iblk), &
00842 aice0(:,:, iblk), works (:,:, :,iblk))
00843
00844 enddo
00845
00846
00847
00848
00849
00850
00851 call ice_timer_start(timer_bound)
00852
00853 call bound_state (aicen, trcrn, &
00854 vicen, vsnon, &
00855 eicen, esnon)
00856
00857 call ice_timer_stop(timer_bound)
00858
00859 call ice_timer_stop(timer_advect)
00860
00861 end subroutine transport_upwind
00862
00863
00864
00865
00866
00867
00868
00869
00870
00871
00872
00873
00874 subroutine state_to_tracers (nx_block, ny_block, &
00875 aice0, &
00876 aicen, trcrn, &
00877 vicen, vsnon, &
00878 eicen, esnon, &
00879 aim, trm)
00880
00881
00882
00883
00884
00885
00886
00887
00888
00889
00890
00891
00892
00893
00894
00895
00896
00897 use ice_itd, only: ilyr1, slyr1
00898
00899
00900
00901 integer (kind=int_kind), intent(in) ::
00902 nx_block, ny_block
00903
00904 real (kind=dbl_kind), dimension (nx_block,ny_block),
00905 intent(in) ::
00906 aice0
00907
00908 real (kind=dbl_kind), dimension (nx_block,ny_block,ncat),
00909 intent(in) ::
00910 aicen ,
00911 vicen ,
00912 vsnon
00913
00914 real (kind=dbl_kind), dimension (nx_block,ny_block,ntrcr,ncat),
00915 intent(in) ::
00916 trcrn
00917
00918 real (kind=dbl_kind), dimension (nx_block,ny_block,ntilyr),
00919 intent(in) ::
00920 eicen
00921
00922 real (kind=dbl_kind), dimension (nx_block,ny_block,ntslyr),
00923 intent(in) ::
00924 esnon
00925
00926 real (kind=dbl_kind), dimension (nx_block,ny_block,0:ncat),
00927 intent(out)::
00928 aim
00929
00930 real (kind=dbl_kind), dimension (nx_block,ny_block,ntrace,ncat),
00931 intent(out) ::
00932 trm
00933
00934 real (kind=dbl_kind), dimension (nx_block,ny_block) ::
00935 worka,
00936 workb
00937
00938
00939
00940 integer (kind=int_kind) ::
00941 i, j, k, n ,
00942 it, kt ,
00943 ij
00944
00945 real (kind=dbl_kind) ::
00946 w1
00947
00948 integer (kind=int_kind), dimension(nx_block*ny_block,0:ncat) ::
00949 indxi ,
00950 indxj
00951
00952 integer (kind=int_kind), dimension(0:ncat) ::
00953 icells
00954
00955 worka(:,:) = c0
00956 workb(:,:) = c0
00957
00958 aim(:,:,0) = aice0(:,:)
00959
00960 do n = 1, ncat
00961
00962 trm(:,:,:,n) = c0
00963
00964
00965
00966
00967
00968 icells(n) = 0
00969 do j = 1, ny_block
00970 do i = 1, nx_block
00971 aim(i,j,n) = aicen(i,j,n)
00972 if (aim(i,j,n) > puny) then
00973 icells(n) = icells(n) + 1
00974 ij = icells(n)
00975 indxi(ij,n) = i
00976 indxj(ij,n) = j
00977 endif
00978 enddo
00979 enddo
00980
00981
00982
00983
00984
00985
00986
00987
00988 do ij = 1, icells(n)
00989 i = indxi(ij,n)
00990 j = indxj(ij,n)
00991 w1 = c1 / aim(i,j,n)
00992 worka(i,j) = c1 / vicen(i,j,n)
00993 trm(i,j,1,n) = vicen(i,j,n) * w1
00994 trm(i,j,2,n) = vsnon(i,j,n) * w1
00995 if (trm(i,j,2,n) > puny) then
00996 workb(i,j) = c1 / vsnon(i,j,n)
00997 else
00998 workb(i,j) = c0
00999 endif
01000 enddo
01001 kt = 2
01002
01003 do it = 1, ntrcr
01004 do ij = 1, icells(n)
01005 i = indxi(ij,n)
01006 j = indxj(ij,n)
01007 trm(i,j,kt+it,n) = trcrn(i,j,it,n)
01008 enddo
01009 enddo
01010 kt = kt + ntrcr
01011
01012 do k =1, nilyr
01013 do ij = 1, icells(n)
01014 i = indxi(ij,n)
01015 j = indxj(ij,n)
01016 trm(i,j,kt+k,n) = eicen(i,j,ilyr1(n)+k-1)*worka(i,j)
01017 enddo
01018 enddo
01019 kt = kt + nilyr
01020
01021 do k = 1, nslyr
01022 do ij = 1, icells(n)
01023 i = indxi(ij,n)
01024 j = indxj(ij,n)
01025 if (trm(i,j,2,n) > puny) &
01026 trm(i,j,kt+k,n) = esnon(i,j,slyr1(n)+k-1)*workb(i,j) &
01027 + rhos*Lfresh
01028 enddo
01029 enddo
01030
01031 enddo
01032
01033 end subroutine state_to_tracers
01034
01035
01036
01037
01038
01039
01040
01041
01042 subroutine tracers_to_state (nx_block, ny_block, &
01043 aim, trm, &
01044 aice0, &
01045 aicen, trcrn, &
01046 vicen, vsnon, &
01047 eicen, esnon)
01048
01049
01050
01051
01052
01053
01054
01055
01056
01057
01058
01059 use ice_itd, only: ilyr1, slyr1
01060
01061
01062
01063 integer (kind=int_kind), intent(in) ::
01064 nx_block, ny_block
01065
01066 real (kind=dbl_kind), dimension (nx_block,ny_block,0:ncat),
01067 intent(in) ::
01068 aim
01069
01070 real (kind=dbl_kind), dimension (nx_block,ny_block,ntrace,ncat),
01071 intent(in) ::
01072 trm
01073
01074 real (kind=dbl_kind), dimension (nx_block,ny_block),
01075 intent(inout) ::
01076 aice0
01077
01078 real (kind=dbl_kind), dimension (nx_block,ny_block,ncat),
01079 intent(inout) ::
01080 aicen ,
01081 vicen ,
01082 vsnon
01083
01084 real (kind=dbl_kind), dimension (nx_block,ny_block,ntrcr,ncat),
01085 intent(inout) ::
01086 trcrn
01087
01088 real (kind=dbl_kind), dimension (nx_block,ny_block,ntilyr),
01089 intent(inout) ::
01090 eicen
01091
01092 real (kind=dbl_kind), dimension (nx_block,ny_block,ntslyr),
01093 intent(inout) ::
01094 esnon
01095
01096
01097
01098 integer (kind=int_kind) ::
01099 i, j, k, n ,
01100 it, kt ,
01101 icells ,
01102 ij
01103
01104 integer (kind=int_kind), dimension (nx_block*ny_block) ::
01105 indxi, indxj
01106
01107 aice0(:,:) = aim(:,:,0)
01108
01109 do n = 1, ncat
01110
01111 icells = 0
01112 do j = 1, ny_block
01113 do i = 1, nx_block
01114 if (aim(i,j,n) > c0) then
01115 icells = icells + 1
01116 indxi(icells) = i
01117 indxj(icells) = j
01118 endif
01119 enddo
01120 enddo
01121
01122
01123
01124
01125
01126 do ij = 1, icells
01127 i = indxi(ij)
01128 j = indxj(ij)
01129 aicen(i,j,n) = aim(i,j,n)
01130 vicen(i,j,n) = aim(i,j,n)*trm(i,j,1,n)
01131 vsnon(i,j,n) = aim(i,j,n)*trm(i,j,2,n)
01132 enddo
01133 kt = 2
01134
01135 do it = 1, ntrcr
01136 do ij = 1, icells
01137 i = indxi(ij)
01138 j = indxj(ij)
01139 trcrn(i,j,it,n) = trm(i,j,kt+it,n)
01140 enddo
01141 enddo
01142 kt = kt + ntrcr
01143
01144 do k = 1, nilyr
01145 do ij = 1, icells
01146 i = indxi(ij)
01147 j = indxj(ij)
01148 eicen(i,j,ilyr1(n)+k-1) = vicen(i,j,n)*trm(i,j,kt+k,n)
01149 enddo
01150 enddo
01151 kt = kt + nilyr
01152
01153 do k = 1, nslyr
01154 do ij = 1, icells
01155 i = indxi(ij)
01156 j = indxj(ij)
01157 esnon(i,j,slyr1(n)+k-1) = (trm(i,j,kt+k,n) - rhos*Lfresh) &
01158 * vsnon(i,j,n)
01159 enddo
01160 enddo
01161
01162 enddo
01163
01164 end subroutine tracers_to_state
01165
01166
01167
01168
01169
01170
01171
01172
01173
01174 subroutine global_conservation (l_stop, &
01175 asum_init, asum_final, &
01176 atsum_init, atsum_final)
01177
01178
01179
01180
01181
01182
01183
01184
01185
01186
01187
01188
01189
01190
01191 real (kind=dbl_kind), intent(in) ::
01192 asum_init ,
01193 asum_final
01194
01195 real (kind=dbl_kind), dimension(ntrace), intent(in), optional ::
01196 atsum_init ,
01197 atsum_final
01198
01199 logical (kind=log_kind), intent(inout) ::
01200 l_stop
01201
01202
01203
01204 integer (kind=int_kind) ::
01205 nt
01206
01207 real (kind=dbl_kind) ::
01208 diff
01209
01210
01211 if (asum_init > puny) then
01212 diff = asum_final - asum_init
01213 if (abs(diff/asum_init) > puny) then
01214 l_stop = .true.
01215 write (nu_diag,*)
01216 write (nu_diag,*) 'Ice area conserv error'
01217 write (nu_diag,*) 'Initial global area =', asum_init
01218 write (nu_diag,*) 'Final global area =', asum_final
01219 write (nu_diag,*) 'Fractional error =', abs(diff)/asum_init
01220 write (nu_diag,*) 'asum_final-asum_init =', diff
01221 endif
01222 endif
01223
01224 if (present(atsum_init)) then
01225 do nt = 1, ntrace
01226 if (abs(atsum_init(nt)) > puny) then
01227 diff = atsum_final(nt) - atsum_init(nt)
01228 if (abs(diff/atsum_init(nt)) > puny) then
01229 l_stop = .true.
01230 write (nu_diag,*)
01231 write (nu_diag,*) 'area*tracer conserv error'
01232 write (nu_diag,*) 'tracer index =', nt
01233 write (nu_diag,*) 'Initial global area*tracer =', &
01234 atsum_init(nt)
01235 write (nu_diag,*) 'Final global area*tracer =', &
01236 atsum_final(nt)
01237 write (nu_diag,*) 'Fractional error =', &
01238 abs(diff)/atsum_init(nt)
01239 write (nu_diag,*) 'atsum_final-atsum_init =', diff
01240 endif
01241 endif
01242 enddo
01243 endif
01244
01245 end subroutine global_conservation
01246
01247
01248
01249
01250
01251
01252
01253
01254 subroutine local_max_min (nx_block, ny_block, &
01255 ilo, ihi, jlo, jhi, &
01256 trm, &
01257 tmin, tmax, &
01258 aimask, trmask)
01259
01260
01261
01262
01263
01264
01265
01266
01267
01268
01269
01270
01271
01272
01273
01274
01275
01276
01277 integer (kind=int_kind), intent(in) ::
01278 nx_block, ny_block,
01279 ilo,ihi,jlo,jhi
01280
01281 real (kind=dbl_kind), intent(in),
01282 dimension(nx_block,ny_block) ::
01283 aimask
01284
01285 real (kind=dbl_kind), intent(in),
01286 dimension (nx_block,ny_block,ntrace) ::
01287 trm ,
01288 trmask
01289
01290 real (kind=dbl_kind), intent(out),
01291 dimension (nx_block,ny_block,ntrace) ::
01292 tmin ,
01293 tmax
01294
01295
01296
01297 integer (kind=int_kind) ::
01298 i, j ,
01299 nt, nt1
01300
01301 real (kind=dbl_kind), dimension(nx_block,ny_block) ::
01302 phimask
01303
01304 real (kind=dbl_kind) ::
01305 phi_nw, phi_n, phi_ne ,
01306 phi_w, phi_e ,
01307 phi_sw, phi_s, phi_se
01308
01309 do nt = 1, ntrace
01310
01311 if (tracer_type(nt)==1) then
01312
01313 do j = 1, ny_block
01314 do i = 1, nx_block
01315 phimask(i,j) = aimask(i,j)
01316 enddo
01317 enddo
01318
01319 else
01320
01321 nt1 = depend(nt)
01322 do j = 1, ny_block
01323 do i = 1, nx_block
01324 phimask(i,j) = trmask(i,j,nt1)
01325 enddo
01326 enddo
01327
01328 endif
01329
01330
01331
01332
01333
01334
01335
01336 do j = jlo, jhi
01337 do i = ilo, ihi
01338
01339 phi_nw = phimask(i-1,j+1) * trm(i-1,j+1,nt) &
01340 + (c1-phimask(i-1,j+1))* trm(i, j, nt)
01341 phi_n = phimask(i, j+1) * trm(i, j+1,nt) &
01342 + (c1-phimask(i, j+1))* trm(i, j, nt)
01343 phi_ne = phimask(i+1,j+1) * trm(i+1,j+1,nt) &
01344 + (c1-phimask(i+1,j+1))* trm(i, j, nt)
01345 phi_w = phimask(i-1,j) * trm(i-1,j, nt) &
01346 + (c1-phimask(i-1,j)) * trm(i, j, nt)
01347 phi_e = phimask(i+1,j) * trm(i+1,j, nt) &
01348 + (c1-phimask(i+1,j)) * trm(i, j, nt)
01349 phi_sw = phimask(i-1,j-1) * trm(i-1,j-1,nt) &
01350 + (c1-phimask(i-1,j-1))* trm(i, j, nt)
01351 phi_s = phimask(i, j-1) * trm(i, j-1,nt) &
01352 + (c1-phimask(i, j-1))* trm(i, j, nt)
01353 phi_se = phimask(i+1,j-1) * trm(i+1,j-1,nt) &
01354 + (c1-phimask(i+1,j-1))* trm(i, j, nt)
01355
01356
01357
01358
01359
01360 tmax(i,j,nt) = max (phi_nw, phi_n, phi_ne, phi_w, &
01361 trm(i,j,nt), phi_e, phi_sw, phi_s, phi_se)
01362
01363 tmin(i,j,nt) = min (phi_nw, phi_n, phi_ne, phi_w, &
01364 trm(i,j,nt), phi_e, phi_sw, phi_s, phi_se)
01365
01366 enddo
01367 enddo
01368
01369 enddo
01370
01371 end subroutine local_max_min
01372
01373
01374
01375
01376
01377
01378
01379
01380 subroutine quasilocal_max_min (nx_block, ny_block, &
01381 ilo, ihi, jlo, jhi, &
01382 tmin, tmax)
01383
01384
01385
01386
01387
01388
01389
01390
01391
01392
01393
01394
01395
01396
01397
01398 integer (kind=int_kind), intent(in) ::
01399 nx_block, ny_block,
01400 ilo,ihi,jlo,jhi
01401
01402 real (kind=dbl_kind), intent(inout),
01403 dimension (nx_block,ny_block,ntrace) ::
01404 tmin ,
01405 tmax
01406
01407
01408
01409 integer (kind=int_kind) ::
01410 i, j ,
01411 nt
01412
01413 do nt = 1, ntrace
01414
01415 do j = jlo, jhi
01416 do i = ilo, ihi
01417
01418 tmax(i,j,nt) = &
01419 max (tmax(i-1,j+1,nt), tmax(i,j+1,nt), tmax(i+1,j+1,nt), &
01420 tmax(i-1,j, nt), tmax(i,j, nt), tmax(i+1,j, nt), &
01421 tmax(i-1,j-1,nt), tmax(i,j-1,nt), tmax(i+1,j-1,nt))
01422
01423 tmin(i,j,nt) = &
01424 min (tmin(i-1,j+1,nt), tmin(i,j+1,nt), tmin(i+1,j+1,nt), &
01425 tmin(i-1,j, nt), tmin(i,j, nt), tmin(i+1,j, nt), &
01426 tmin(i-1,j-1,nt), tmin(i,j-1,nt), tmin(i+1,j-1,nt))
01427
01428 enddo
01429 enddo
01430
01431 enddo
01432
01433 end subroutine quasilocal_max_min
01434
01435
01436
01437
01438
01439
01440
01441
01442
01443 subroutine check_monotonicity (nx_block, ny_block, &
01444 ilo, ihi, jlo, jhi, &
01445 iblk, &
01446 tmin, tmax, &
01447 aim, trm, &
01448 l_stop, &
01449 istop, jstop)
01450
01451
01452
01453
01454
01455
01456
01457
01458
01459
01460
01461
01462
01463
01464 integer (kind=int_kind), intent(in) ::
01465 nx_block, ny_block,
01466 ilo,ihi,jlo,jhi ,
01467 iblk
01468
01469 real (kind=dbl_kind), intent(in),
01470 dimension (nx_block,ny_block) ::
01471 aim
01472
01473 real (kind=dbl_kind), intent(in),
01474 dimension (nx_block,ny_block,ntrace) ::
01475 trm
01476
01477 real (kind=dbl_kind), intent(in),
01478 dimension (nx_block,ny_block,ntrace) ::
01479 tmin ,
01480 tmax
01481
01482 logical (kind=log_kind), intent(inout) ::
01483 l_stop
01484
01485 integer (kind=int_kind), intent(inout) ::
01486 istop, jstop
01487
01488
01489
01490 integer (kind=int_kind) ::
01491 i, j ,
01492 nt, nt1, nt2
01493
01494 real (kind=dbl_kind) ::
01495 w1, w2
01496
01497 logical (kind=log_kind), dimension (nx_block, ny_block) ::
01498 l_check
01499
01500 do nt = 1, ntrace
01501
01502
01503
01504
01505
01506 if (tracer_type(nt)==1) then
01507
01508 do j = jlo, jhi
01509 do i = ilo, ihi
01510 if (aim(i,j) > puny) then
01511 l_check(i,j) = .true.
01512 else
01513 l_check(i,j) = .false.
01514 endif
01515 enddo
01516 enddo
01517
01518 elseif (tracer_type(nt)==2) then
01519
01520 nt1 = depend(nt)
01521 do j = jlo, jhi
01522 do i = ilo, ihi
01523 if (abs(trm(i,j,nt1)) > puny .and. aim(i,j) > puny) then
01524 l_check(i,j) = .true.
01525 else
01526 l_check(i,j) = .false.
01527 endif
01528 enddo
01529 enddo
01530
01531 elseif (tracer_type(nt)==3) then
01532
01533 nt1 = depend(nt)
01534 nt2 = depend(nt1)
01535 do j = jlo, jhi
01536 do i = ilo, ihi
01537 if (abs(trm(i,j,nt1)) > puny .and. &
01538 abs(trm(i,j,nt2)) > puny .and. aim(i,j) > puny) then
01539 l_check(i,j) = .true.
01540 else
01541 l_check(i,j) = .false.
01542 endif
01543 enddo
01544 enddo
01545 endif
01546
01547
01548
01549
01550
01551 do j = jlo, jhi
01552 do i = ilo, ihi
01553
01554 if (l_check(i,j)) then
01555
01556 w1 = max(c1, abs(tmin(i,j,nt)))
01557 w2 = max(c1, abs(tmax(i,j,nt)))
01558 if (trm(i,j,nt) < tmin(i,j,nt)-w1*puny) then
01559 l_stop = .true.
01560 istop = i
01561 jstop = j
01562 write (nu_diag,*) ' '
01563 write (nu_diag,*) 'new tracer < tmin'
01564 write (nu_diag,*) 'i, j, nt =', i, j, nt
01565 write (nu_diag,*) 'new tracer =', trm (i,j,nt)
01566 write (nu_diag,*) 'tmin =' , tmin(i,j,nt)
01567 write (nu_diag,*) 'ice area =' , aim(i,j)
01568 elseif (trm(i,j,nt) > tmax(i,j,nt)+w2*puny) then
01569 l_stop = .true.
01570 istop = i
01571 jstop = j
01572 write (nu_diag,*) ' '
01573 write (nu_diag,*) 'new tracer > tmax'
01574 write (nu_diag,*) 'i, j, nt =', i, j, nt
01575 write (nu_diag,*) 'new tracer =', trm (i,j,nt)
01576 write (nu_diag,*) 'tmax =' , tmax(i,j,nt)
01577 write (nu_diag,*) 'ice area =' , aim(i,j)
01578 endif
01579 endif
01580
01581 enddo
01582 enddo
01583
01584 enddo
01585
01586 end subroutine check_monotonicity
01587
01588
01589
01590
01591
01592
01593
01594
01595
01596
01597 subroutine state_to_work (nx_block, ny_block, &
01598 narr, trcr_depend, &
01599 aicen, trcrn, &
01600 vicen, vsnon, &
01601 aice0, works)
01602
01603
01604
01605
01606
01607
01608
01609
01610
01611
01612
01613
01614 use ice_itd, only: ilyr1, slyr1
01615
01616
01617
01618 integer (kind=int_kind), intent(in) ::
01619 nx_block, ny_block ,
01620 narr
01621
01622 integer (kind=int_kind), dimension (ntrcr), intent(in) ::
01623 trcr_depend
01624
01625 real (kind=dbl_kind), dimension (nx_block,ny_block,ncat),
01626 intent(in) ::
01627 aicen ,
01628 vicen ,
01629 vsnon
01630
01631 real (kind=dbl_kind), dimension (nx_block,ny_block,ntrcr,ncat),
01632 intent(in) ::
01633 trcrn
01634
01635 real (kind=dbl_kind), dimension (nx_block,ny_block),
01636 intent(in) ::
01637 aice0
01638
01639 real (kind=dbl_kind), dimension(nx_block,ny_block,narr),
01640 intent (out) ::
01641 works
01642
01643
01644
01645 integer (kind=int_kind) ::
01646 i, j, k, n, it ,
01647 narrays
01648
01649
01650
01651
01652
01653
01654
01655
01656 do j = 1, ny_block
01657 do i = 1, nx_block
01658 works(i,j,1) = aice0(i,j)
01659 enddo
01660 enddo
01661 narrays = 1
01662
01663 do n=1, ncat
01664
01665 do j = 1, ny_block
01666 do i = 1, nx_block
01667 works(i,j,narrays+1) = aicen(i,j,n)
01668 works(i,j,narrays+2) = vicen(i,j,n)
01669 works(i,j,narrays+3) = vsnon(i,j,n)
01670 enddo
01671 enddo
01672 narrays = narrays + 3
01673
01674 do it = 1, ntrcr
01675 if (trcr_depend(it) == 0) then
01676 do j = 1, ny_block
01677 do i = 1, nx_block
01678 works(i,j,narrays+it) = aicen(i,j,n)*trcrn(i,j,it,n)
01679 enddo
01680 enddo
01681 elseif (trcr_depend(it) == 1) then
01682 do j = 1, ny_block
01683 do i = 1, nx_block
01684 works(i,j,narrays+it) = vicen(i,j,n)*trcrn(i,j,it,n)
01685 enddo
01686 enddo
01687 elseif (trcr_depend(it) == 2) then
01688 do j = 1, ny_block
01689 do i = 1, nx_block
01690 works(i,j,narrays+it) = vsnon(i,j,n)*trcrn(i,j,it,n)
01691 enddo
01692 enddo
01693 endif
01694 enddo
01695 narrays = narrays + ntrcr
01696
01697 enddo
01698
01699 if (narr /= narrays) write(nu_diag,*) &
01700 "Wrong number of arrays in transport bound call"
01701
01702 end subroutine state_to_work
01703
01704
01705
01706
01707
01708
01709
01710
01711 subroutine work_to_state (nx_block, ny_block, &
01712 narr, trcr_depend, &
01713 aicen, trcrn, &
01714 vicen, vsnon, &
01715 aice0, works)
01716
01717
01718
01719
01720
01721
01722
01723
01724
01725
01726
01727
01728 use ice_itd, only: ilyr1, slyr1, compute_tracers
01729
01730
01731
01732 integer (kind=int_kind), intent (in) ::
01733 nx_block, ny_block ,
01734 narr
01735
01736 integer (kind=int_kind), dimension (ntrcr), intent(in) ::
01737 trcr_depend
01738
01739 real (kind=dbl_kind), intent (in) ::
01740 works (nx_block,ny_block,narr)
01741
01742 real (kind=dbl_kind), dimension (nx_block,ny_block,ncat),
01743 intent(out) ::
01744 aicen ,
01745 vicen ,
01746 vsnon
01747
01748 real (kind=dbl_kind), dimension (nx_block,ny_block,ntrcr,ncat),
01749 intent(out) ::
01750 trcrn
01751
01752 real (kind=dbl_kind), dimension (nx_block,ny_block),
01753 intent(out) ::
01754 aice0
01755
01756
01757
01758 integer (kind=int_kind) ::
01759 i, j, k, n , it,
01760 narrays ,
01761 icells
01762
01763 integer (kind=int_kind), dimension (nx_block*ny_block) ::
01764 indxi, indxj
01765
01766 real (kind=dbl_kind), dimension (nx_block*ny_block,narr) ::
01767 work
01768
01769
01770 icells = 0
01771 do j = 1, ny_block
01772 do i = 1, nx_block
01773 icells = icells + 1
01774 indxi(icells) = i
01775 indxj(icells) = j
01776 work (icells,:) = works(i,j,:)
01777 enddo
01778 enddo
01779
01780 do j=1,ny_block
01781 do i=1,nx_block
01782 aice0(i,j) = works(i,j,1)
01783 enddo
01784 enddo
01785 narrays = 1
01786
01787 do n=1,ncat
01788
01789 do j=1,ny_block
01790 do i=1,nx_block
01791 aicen(i,j,n) = works(i,j,narrays+1)
01792 vicen(i,j,n) = works(i,j,narrays+2)
01793 vsnon(i,j,n) = works(i,j,narrays+3)
01794 enddo
01795 enddo
01796 narrays = narrays + 3
01797
01798 call compute_tracers (nx_block, ny_block, &
01799 icells, indxi, indxj, &
01800 trcr_depend, &
01801 work (:,narrays+1:narrays+ntrcr), &
01802 aicen(:,:,n), &
01803 vicen(:,:,n), vsnon(:,:,n), &
01804 trcrn(:,:,:,n))
01805
01806 narrays = narrays + ntrcr
01807
01808 enddo
01809
01810 end subroutine work_to_state
01811
01812
01813
01814
01815
01816
01817
01818
01819 subroutine upwind_field (nx_block, ny_block, &
01820 ilo, ihi, jlo, jhi, &
01821 dt, &
01822 narrays, phi, &
01823 uee, vnn, &
01824 HTE, HTN, &
01825 tarea)
01826
01827
01828
01829
01830
01831
01832
01833
01834
01835
01836
01837
01838
01839
01840 integer (kind=int_kind), intent (in) ::
01841 nx_block, ny_block ,
01842 ilo,ihi,jlo,jhi ,
01843 narrays
01844
01845 real (kind=dbl_kind), intent(in) ::
01846 dt
01847
01848 real (kind=dbl_kind), dimension(nx_block,ny_block,narrays),
01849 intent(inout) ::
01850 phi
01851
01852 real (kind=dbl_kind), dimension(nx_block,ny_block),
01853 intent(in)::
01854 uee, vnn
01855
01856 real (kind=dbl_kind), dimension(nx_block,ny_block), intent(in) ::
01857 HTE ,
01858 HTN ,
01859 tarea
01860
01861
01862
01863 integer (kind=int_kind) ::
01864 i, j, k, n
01865
01866 real (kind=dbl_kind) ::
01867 upwind, y1, y2, a, h
01868
01869 real (kind=dbl_kind), dimension (nx_block,ny_block) ::
01870 worka,
01871 workb
01872
01873
01874
01875
01876
01877 upwind(y1,y2,a,h) = p5*dt*h*((a+abs(a))*y1+(a-abs(a))*y2)
01878
01879
01880
01881
01882
01883 worka(:,:) = c0
01884 workb(:,:) = c0
01885
01886 do n = 1, narrays
01887
01888 do j = 1, jhi
01889 do i = 1, ihi
01890 worka(i,j)= &
01891 upwind(phi(i,j,n),phi(i+1,j,n),uee(i,j),HTE(i,j))
01892 workb(i,j)= &
01893 upwind(phi(i,j,n),phi(i,j+1,n),vnn(i,j),HTN(i,j))
01894 enddo
01895 enddo
01896
01897 do j = jlo, jhi
01898 do i = ilo, ihi
01899 phi(i,j,n) = phi(i,j,n) - ( worka(i,j)-worka(i-1,j) &
01900 + workb(i,j)-workb(i,j-1) ) &
01901 / tarea(i,j)
01902 enddo
01903 enddo
01904
01905 enddo
01906
01907 end subroutine upwind_field
01908
01909
01910
01911 end module ice_transport_driver
01912
01913