00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038 module ice_transport_remap
00039
00040
00041
00042 use ice_kinds_mod
00043 use ice_communicate, only: my_task, master_task, MPI_COMM_ICE
00044 use ice_domain_size
00045 use ice_constants
00046 use ice_fileunits, only: nu_diag
00047 use perf_mod, only: t_startf, t_stopf, t_barrierf
00048
00049
00050
00051 implicit none
00052 save
00053 private
00054 public :: init_remap, horizontal_remap, make_masks
00055
00056
00057
00058
00059
00060
00061 integer (kind=int_kind), parameter ::
00062 ntrace = 2+ntrcr+nilyr+nslyr
00063
00064 integer (kind=int_kind), parameter ::
00065 ngroups = 6 ,
00066
00067 nvert = 3
00068
00069
00070 real (kind=dbl_kind), parameter ::
00071 p5625m = -9._dbl_kind/16._dbl_kind ,
00072 p52083 = 25._dbl_kind/48._dbl_kind
00073
00074 logical (kind=log_kind), parameter :: bugcheck = .false.
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105
00106
00107
00108
00109
00110
00111
00112
00113
00114
00115
00116
00117
00118
00119
00120
00121
00122
00123
00124
00125
00126
00127
00128
00129
00130
00131
00132
00133
00134
00135
00136
00137
00138
00139
00140
00141
00142
00143
00144
00145
00146
00147
00148
00149
00150
00151
00152
00153
00154
00155
00156
00157
00158
00159
00160
00161
00162
00163
00164
00165
00166
00167
00168
00169
00170
00171
00172
00173
00174
00175
00176
00177
00178
00179
00180
00181
00182
00183
00184
00185
00186
00187
00188
00189
00190
00191
00192
00193
00194
00195
00196
00197
00198
00199
00200
00201
00202
00203
00204
00205
00206
00207
00208
00209
00210
00211
00212
00213
00214
00215
00216
00217
00218
00219
00220
00221
00222
00223
00224
00225
00226
00227
00228
00229
00230
00231
00232
00233
00234
00235
00236
00237
00238
00239
00240
00241
00242
00243
00244
00245
00246
00247
00248
00249
00250
00251
00252
00253
00254
00255
00256
00257
00258
00259
00260
00261
00262
00263
00264
00265 contains
00266
00267
00268
00269
00270
00271
00272
00273
00274
00275 subroutine init_remap
00276
00277
00278
00279
00280
00281
00282
00283
00284
00285
00286
00287 use ice_boundary
00288 use ice_domain
00289 use ice_blocks
00290 use ice_grid, only: dxt, dyt, &
00291 xav, yav, xxav, xyav, yyav, &
00292 xxxav, xxyav, xyyav, yyyav
00293 use ice_exit
00294
00295
00296
00297
00298
00299 integer (kind=int_kind) ::
00300 i, j, iblk
00301
00302
00303
00304
00305
00306
00307
00308
00309 do iblk = 1, nblocks
00310 do j = 1, ny_block
00311 do i = 1, nx_block
00312 xav(i,j,iblk) = c0
00313 yav(i,j,iblk) = c0
00314
00315
00316
00317
00318 xxav(i,j,iblk) = c1/c12
00319 yyav(i,j,iblk) = c1/c12
00320 xyav(i,j,iblk) = c0
00321 xxxav(i,j,iblk) = c0
00322 xxyav(i,j,iblk) = c0
00323 xyyav(i,j,iblk) = c0
00324 yyyav(i,j,iblk) = c0
00325 enddo
00326 enddo
00327 enddo
00328
00329
00330 end subroutine init_remap
00331
00332
00333
00334
00335
00336
00337
00338
00339 subroutine horizontal_remap (dt, &
00340 uvel, vvel, &
00341 mm, tm, &
00342 l_fixed_area, &
00343 edgearea_e, edgearea_n, &
00344 tracer_type_in, depend_in, &
00345 has_dependents_in, &
00346 integral_order_in, &
00347 l_dp_midpt_in)
00348
00349
00350
00351
00352
00353
00354
00355
00356
00357
00358
00359
00360
00361
00362
00363
00364
00365
00366
00367
00368
00369
00370
00371 use ice_boundary
00372 use ice_global_reductions
00373 use ice_domain
00374 use ice_blocks
00375 use ice_grid, only: HTE, HTN, dxt, dyt, dxu, dyu, &
00376 tarea, tarear, hm, &
00377 xav, yav, xxav, xyav, yyav, &
00378 xxxav, xxyav, xyyav, yyyav
00379 use ice_exit
00380 use ice_calendar, only: istep1
00381 use ice_timers
00382
00383
00384
00385 real (kind=dbl_kind), intent(in) ::
00386 dt
00387
00388 real (kind=dbl_kind), intent(in),
00389 dimension(nx_block,ny_block,max_blocks) ::
00390 uvel ,
00391 vvel
00392
00393 real (kind=dbl_kind), intent(inout),
00394 dimension (nx_block,ny_block,0:ncat,max_blocks) ::
00395 mm
00396
00397 real (kind=dbl_kind), intent(inout),
00398 dimension (nx_block,ny_block,ntrace,ncat,max_blocks) ::
00399 tm
00400
00401
00402
00403
00404
00405
00406
00407
00408
00409 logical, intent(in) ::
00410 l_fixed_area
00411
00412
00413 real (kind=dbl_kind), dimension(nx_block,ny_block,max_blocks),
00414 intent(inout) ::
00415 edgearea_e ,
00416 edgearea_n
00417
00418 integer (kind=int_kind), dimension (ntrace), intent(in),
00419 optional ::
00420 tracer_type_in ,
00421 depend_in
00422
00423 logical (kind=log_kind), dimension (ntrace), intent(in),
00424 optional ::
00425 has_dependents_in
00426
00427 integer (kind=int_kind), intent(in), optional ::
00428 integral_order_in
00429
00430 logical (kind=log_kind), intent(in), optional ::
00431 l_dp_midpt_in
00432
00433
00434
00435
00436
00437
00438 integer (kind=int_kind), dimension (ntrace) ::
00439 tracer_type ,
00440 depend
00441
00442 logical (kind=log_kind), dimension (ntrace) ::
00443 has_dependents
00444
00445 integer (kind=int_kind) ::
00446 integral_order
00447
00448 logical (kind=log_kind) ::
00449 l_dp_midpt
00450
00451
00452 integer (kind=int_kind) ::
00453 i, j ,
00454 iblk ,
00455 ilo,ihi,jlo,jhi,
00456 n
00457
00458 integer (kind=int_kind), dimension(0:ncat,max_blocks) ::
00459 icellsnc
00460
00461 integer (kind=int_kind),
00462 dimension(nx_block*ny_block,0:ncat,max_blocks) ::
00463 indxinc, indxjnc
00464
00465 type (block) ::
00466 this_block
00467
00468 real (kind=dbl_kind), dimension (nx_block,ny_block,max_blocks) ::
00469 dpx ,
00470 dpy
00471
00472 real (kind=dbl_kind), dimension(nx_block,ny_block,0:ncat,max_blocks) ::
00473 mc ,
00474 mx, my ,
00475 mmask
00476
00477 real (kind=dbl_kind),
00478 dimension (nx_block,ny_block,ntrace,ncat,max_blocks) ::
00479 tc ,
00480 tx, ty ,
00481 tmask
00482
00483 real (kind=dbl_kind), dimension (nx_block,ny_block,0:ncat,max_blocks) ::
00484 mflxe, mflxn
00485
00486 real (kind=dbl_kind), dimension (nx_block,ny_block,ntrace,ncat,max_blocks) ::
00487 mtflxe, mtflxn
00488
00489 real (kind=dbl_kind), dimension (nx_block,ny_block,ngroups,max_blocks) ::
00490 triarea
00491
00492 real (kind=dbl_kind), dimension (nx_block,ny_block,0:nvert,ngroups,max_blocks) ::
00493 xp, yp
00494
00495
00496 integer (kind=int_kind),
00497 dimension (nx_block,ny_block,ngroups,max_blocks) ::
00498 iflux ,
00499 jflux
00500
00501 integer (kind=int_kind), dimension(ngroups,max_blocks) ::
00502 icellsng
00503
00504 integer (kind=int_kind),
00505 dimension(nx_block*ny_block,ngroups,max_blocks) ::
00506 indxing, indxjng
00507
00508 logical (kind=log_kind), dimension(max_blocks) ::
00509 l_stop
00510
00511 integer (kind=int_kind), dimension(max_blocks) ::
00512 istop, jstop
00513
00514 character (len=char_len), dimension(max_blocks) ::
00515 edge
00516
00517 real (kind=dbl_kind), dimension (nx_block,ny_block) ::
00518 worka,
00519 workb,
00520 workc,
00521 workd
00522
00523 real (kind=dbl_kind), dimension (nx_block,ny_block,2,max_blocks) ::
00524 dpwork
00525
00526 real (kind=dbl_kind), dimension(nx_block,ny_block,2,0:ncat,max_blocks) ::
00527 mwork
00528
00529
00530
00531
00532 l_stop = .false.
00533 istop = 0
00534 jstop = 0
00535
00536
00537
00538
00539
00540
00541
00542 if (present(tracer_type_in)) then
00543 tracer_type(:) = tracer_type_in(:)
00544 else
00545 tracer_type(:) = 1
00546 endif
00547
00548 if (present(depend_in)) then
00549 depend(:) = depend_in(:)
00550 else
00551 depend(:) = 0
00552 endif
00553
00554 if (present(has_dependents_in)) then
00555 has_dependents(:) = has_dependents_in(:)
00556 else
00557 has_dependents(:) = .false.
00558 endif
00559
00560 if (present(integral_order_in)) then
00561 integral_order = integral_order_in
00562 else
00563 integral_order = 2
00564 endif
00565
00566 if (present(l_dp_midpt_in)) then
00567 l_dp_midpt = l_dp_midpt_in
00568 else
00569 l_dp_midpt = .false.
00570 endif
00571
00572 worka(:,:) = c1
00573 workb(:,:) = c1
00574 workc(:,:) = c1
00575 workd(:,:) = c1
00576
00577
00578
00579
00580
00581
00582
00583 do iblk = 1, nblocks
00584
00585 this_block = get_block(blocks_ice(iblk),iblk)
00586 ilo = this_block%ilo
00587 ihi = this_block%ihi
00588 jlo = this_block%jlo
00589 jhi = this_block%jhi
00590
00591
00592
00593
00594
00595
00596
00597 call make_masks (nx_block, ny_block, &
00598 ilo, ihi, jlo, jhi, &
00599 nghost, has_dependents, &
00600 icellsnc (:,iblk), &
00601 indxinc(:,:,iblk), indxjnc (:,:,iblk), &
00602 mm (:,:,:,iblk), mmask (:,:,:,iblk), &
00603 tm (:,:,:,:,iblk), tmask(:,:,:,:,iblk))
00604
00605
00606
00607
00608
00609
00610
00611
00612
00613 call construct_fields(nx_block, ny_block, &
00614 ilo, ihi, jlo, jhi, &
00615 nghost, &
00616 tracer_type, depend, &
00617 has_dependents, icellsnc (0,iblk), &
00618 indxinc(:,0,iblk), indxjnc(:,0,iblk), &
00619
00620 worka (:,:), workb (:,:), &
00621 hm (:,:,iblk), xav (:,:,iblk), &
00622 yav (:,:,iblk), xxav (:,:,iblk), &
00623 xyav (:,:,iblk), yyav (:,:,iblk), &
00624 xxxav (:,:,iblk), xxyav (:,:,iblk), &
00625 xyyav (:,:,iblk), yyyav (:,:,iblk), &
00626
00627 workc (:,:), workd (:,:), &
00628 mm (:,:,0,iblk), mc (:,:,0,iblk), &
00629 mx (:,:,0,iblk), my (:,:,0,iblk), &
00630 mmask(:,:,0,iblk) )
00631
00632
00633
00634 do n = 1, ncat
00635
00636 call construct_fields(nx_block, ny_block, &
00637 ilo, ihi, jlo, jhi, &
00638 nghost, &
00639 tracer_type, depend, &
00640 has_dependents, icellsnc (n,iblk), &
00641 indxinc (:,n,iblk), indxjnc(:,n,iblk), &
00642
00643 worka (:,:), workb (:,:), &
00644 hm (:,:,iblk), xav (:,:,iblk), &
00645 yav (:,:,iblk), xxav (:,:,iblk), &
00646 xyav (:,:,iblk), yyav (:,:,iblk), &
00647 xxxav (:,:,iblk), xxyav (:,:,iblk), &
00648 xyyav (:,:,iblk), yyyav (:,:,iblk), &
00649
00650 workc (:,:), workd (:,:), &
00651 mm (:,:,n,iblk), mc (:,:,n,iblk), &
00652 mx (:,:,n,iblk), my (:,:,n,iblk), &
00653 mmask (:,:,n,iblk), &
00654 tm (:,:,:,n,iblk), tc (:,:,:,n,iblk), &
00655 tx (:,:,:,n,iblk), ty (:,:,:,n,iblk), &
00656 tmask(:,:,:,n,iblk) )
00657
00658 enddo
00659
00660
00661
00662
00663
00664
00665 call departure_points(nx_block, ny_block, &
00666 ilo, ihi, jlo, jhi, &
00667 nghost, dt, &
00668 uvel(:,:,iblk), vvel(:,:,iblk), &
00669 dxu (:,:,iblk), dyu (:,:,iblk), &
00670 HTN (:,:,iblk), HTE (:,:,iblk), &
00671 dpx (:,:,iblk), dpy (:,:,iblk), &
00672 l_dp_midpt, l_stop (iblk), &
00673 istop (iblk), jstop (iblk))
00674
00675 if (l_stop(iblk)) then
00676 this_block = get_block(blocks_ice(iblk),iblk)
00677 write(nu_diag,*) 'istep1, my_task, iblk =', &
00678 istep1, my_task, iblk
00679 write (nu_diag,*) 'Global block:', this_block%block_id
00680 if (istop(iblk) > 0 .and. jstop(iblk) > 0) &
00681 write(nu_diag,*) 'Global i and j:', &
00682 this_block%i_glob(istop(iblk)), &
00683 this_block%j_glob(jstop(iblk))
00684 call abort_ice('remap transport: bad departure points')
00685 endif
00686
00687 enddo
00688
00689
00690
00691
00692
00693
00694
00695
00696 if (nghost==1) then
00697
00698
00699
00700 call ice_timer_start(timer_bound)
00701
00702
00703
00704
00705
00706 do iblk = 1, nblocks
00707 do j = 1, ny_block
00708 do i = 1, nx_block
00709 dpwork(i,j,1, iblk) = dpx(i,j,iblk)
00710 dpwork(i,j,2, iblk) = dpy(i,j,iblk)
00711 mwork (i,j,1,:,iblk) = mx (i,j,:,iblk)
00712 mwork (i,j,2,:,iblk) = my (i,j,:,iblk)
00713 enddo
00714 enddo
00715 enddo
00716
00717
00718
00719
00720
00721 call ice_HaloUpdate (dpwork, halo_info, &
00722 field_loc_NEcorner, field_type_vector)
00723
00724
00725 call ice_HaloUpdate (mc, halo_info, &
00726 field_loc_center, field_type_scalar)
00727
00728
00729
00730
00731 call ice_HaloUpdate (mwork, halo_info, &
00732 field_loc_center, field_type_vector)
00733
00734
00735 do iblk = 1, nblocks
00736 do j = 1, ny_block
00737 do i = 1, nx_block
00738 dpx(i,j, iblk) = dpwork(i,j,1, iblk)
00739 dpy(i,j, iblk) = dpwork(i,j,2, iblk)
00740 mx (i,j,:,iblk) = mwork (i,j,1,:,iblk)
00741 my (i,j,:,iblk) = mwork (i,j,2,:,iblk)
00742 enddo
00743 enddo
00744 enddo
00745
00746
00747
00748
00749
00750
00751
00752 call ice_HaloUpdate (tc, halo_info, &
00753 field_loc_center, field_type_scalar)
00754 call ice_HaloUpdate (tx, halo_info, &
00755 field_loc_center, field_type_vector)
00756 call ice_HaloUpdate (ty, halo_info, &
00757 field_loc_center, field_type_vector)
00758 call ice_timer_stop(timer_bound)
00759
00760
00761 endif
00762
00763
00764
00765
00766
00767 do iblk = 1, nblocks
00768
00769 this_block = get_block(blocks_ice(iblk),iblk)
00770 ilo = this_block%ilo
00771 ihi = this_block%ihi
00772 jlo = this_block%jlo
00773 jhi = this_block%jhi
00774
00775
00776
00777
00778
00779
00780
00781
00782
00783 edge(iblk) = 'east'
00784 call locate_triangles(nx_block, ny_block, &
00785 ilo, ihi, jlo, jhi, &
00786 nghost, edge(iblk), &
00787 icellsng (:,iblk), &
00788 indxing (:,:,iblk), indxjng(:,:,iblk), &
00789 dpx (:,:,iblk), dpy (:,:,iblk), &
00790 dxu (:,:,iblk), dyu (:,:,iblk), &
00791 xp (:,:,:,:,iblk), yp (:,:,:,:,iblk), &
00792 iflux (:,:,:,iblk), jflux(:,:,:,iblk), &
00793 triarea (:,:,:,iblk), l_fixed_area, &
00794 edgearea_e(:,:,iblk))
00795
00796
00797
00798
00799
00800
00801 call triangle_coordinates (nx_block, ny_block, &
00802 integral_order, icellsng (:,iblk), &
00803 indxing(:,:,iblk), indxjng(:,:,iblk), &
00804 xp (:,:,:,:,iblk), yp (:,:,:,:,iblk))
00805
00806
00807
00808
00809
00810
00811
00812 call transport_integrals(nx_block, ny_block, &
00813 icellsng (:,iblk), &
00814 indxing(:,:,iblk), indxjng (:,:,iblk), &
00815 tracer_type, depend, &
00816 integral_order, triarea(:,:,:,iblk), &
00817 iflux(:,:,:,iblk), jflux (:,:,:,iblk), &
00818 xp (:,:,:,:,iblk), yp (:,:,:,:,iblk), &
00819 mc (:,:,0,iblk), mx (:,:,0,iblk), &
00820 my (:,:,0,iblk), mflxe (:,:,0,iblk))
00821
00822
00823 do n = 1, ncat
00824 call transport_integrals &
00825 (nx_block, ny_block, &
00826 icellsng (:,iblk), &
00827 indxing(:,:,iblk), indxjng (:,:,iblk), &
00828 tracer_type, depend, &
00829 integral_order, triarea (:,:,:,iblk), &
00830 iflux(:,:,:,iblk), jflux (:,:,:,iblk), &
00831 xp (:,:,:,:,iblk), yp (:,:,:,:,iblk), &
00832 mc (:,:,n,iblk), mx (:,:,n,iblk), &
00833 my (:,:,n,iblk), mflxe (:,:,n,iblk), &
00834 tc (:,:,:,n,iblk), tx (:,:,:,n,iblk), &
00835 ty (:,:,:,n,iblk), mtflxe(:,:,:,n,iblk))
00836
00837 enddo
00838
00839
00840
00841
00842
00843 edge(iblk) = 'north'
00844 call locate_triangles(nx_block, ny_block, &
00845 ilo, ihi, jlo, jhi, &
00846 nghost, edge(iblk), &
00847 icellsng (:,iblk), &
00848 indxing (:,:,iblk), indxjng(:,:,iblk), &
00849 dpx (:,:,iblk), dpy (:,:,iblk), &
00850 dxu (:,:,iblk), dyu (:,:,iblk), &
00851 xp (:,:,:,:,iblk), yp (:,:,:,:,iblk), &
00852 iflux (:,:,:,iblk), jflux(:,:,:,iblk), &
00853 triarea (:,:,:,iblk), l_fixed_area, &
00854 edgearea_n(:,:,iblk))
00855
00856 call triangle_coordinates (nx_block, ny_block, &
00857 integral_order, icellsng (:,iblk), &
00858 indxing(:,:,iblk), indxjng(:,:,iblk), &
00859 xp (:,:,:,:,iblk), yp (:,:,:,:,iblk))
00860
00861
00862 call transport_integrals(nx_block, ny_block, &
00863 icellsng (:,iblk), &
00864 indxing(:,:,iblk), indxjng(:,:,iblk), &
00865 tracer_type, depend, &
00866 integral_order, triarea(:,:,:,iblk), &
00867 iflux(:,:,:,iblk), jflux (:,:,:,iblk), &
00868 xp (:,:,:,:,iblk), yp (:,:,:,:,iblk), &
00869 mc (:,:,0,iblk), mx (:,:,0,iblk), &
00870 my (:,:,0,iblk), mflxn (:,:,0,iblk))
00871
00872
00873 do n = 1, ncat
00874 call transport_integrals &
00875 (nx_block, ny_block, &
00876 icellsng (:,iblk), &
00877 indxing(:,:,iblk), indxjng (:,:,iblk), &
00878 tracer_type, depend, &
00879 integral_order, triarea (:,:,:,iblk), &
00880 iflux(:,:,:,iblk), jflux (:,:,:,iblk), &
00881 xp (:,:,:,:,iblk), yp (:,:,:,:,iblk), &
00882 mc (:,:,n,iblk), mx (:,:,n,iblk), &
00883 my (:,:,n,iblk), mflxn (:,:,n,iblk), &
00884 tc (:,:,:,n,iblk), tx (:,:,:,n,iblk), &
00885 ty (:,:,:,n,iblk), mtflxn(:,:,:,n,iblk))
00886
00887 enddo
00888
00889
00890
00891
00892
00893
00894
00895 call update_fields (nx_block, ny_block, &
00896 ilo, ihi, jlo, jhi, &
00897 tracer_type, depend, &
00898 tarear(:,:,iblk), l_stop(iblk), &
00899 istop(iblk), jstop (iblk), &
00900 mflxe(:,:,0,iblk), mflxn(:,:,0,iblk), &
00901 mm (:,:,0,iblk))
00902
00903 if (l_stop(iblk)) then
00904 this_block = get_block(blocks_ice(iblk),iblk)
00905 write (nu_diag,*) 'istep1, my_task, iblk, cat =', &
00906 istep1, my_task, iblk, '0'
00907 write (nu_diag,*) 'Global block:', this_block%block_id
00908 if (istop(iblk) > 0 .and. jstop(iblk) > 0) &
00909 write(nu_diag,*) 'Global i and j:', &
00910 this_block%i_glob(istop(iblk)), &
00911 this_block%j_glob(jstop(iblk))
00912 call abort_ice ('ice remap_transport: negative area (open water)')
00913 endif
00914
00915
00916
00917 do n = 1, ncat
00918
00919 call update_fields(nx_block, ny_block, &
00920 ilo, ihi, jlo, jhi, &
00921 tracer_type, depend, &
00922 tarear(:,:,iblk), l_stop(iblk), &
00923 istop(iblk), jstop(iblk), &
00924 mflxe (:,:, n,iblk), mflxn (:,:, n,iblk), &
00925 mm (:,:, n,iblk), &
00926 mtflxe(:,:,:,n,iblk), mtflxn(:,:,:,n,iblk), &
00927 tm (:,:,:,n,iblk))
00928
00929 if (l_stop(iblk)) then
00930 this_block = get_block(blocks_ice(iblk),iblk)
00931 write (nu_diag,*) 'istep1, my_task, iblk, cat =', &
00932 istep1, my_task, iblk, n
00933 write (nu_diag,*) 'Global block:', this_block%block_id
00934 if (istop(iblk) > 0 .and. jstop(iblk) > 0) &
00935 write(nu_diag,*) 'Global i and j:', &
00936 this_block%i_glob(istop(iblk)), &
00937 this_block%j_glob(jstop(iblk))
00938 call abort_ice ('ice remap_transport: negative area (ice)')
00939 endif
00940 enddo
00941
00942 enddo
00943
00944
00945
00946
00947 end subroutine horizontal_remap
00948
00949
00950
00951
00952
00953
00954
00955
00956
00957 subroutine make_masks (nx_block, ny_block, &
00958 ilo, ihi, jlo, jhi, &
00959 nghost, has_dependents, &
00960 icells, &
00961 indxi, indxj, &
00962 mm, mmask, &
00963 tm, tmask)
00964
00965
00966
00967
00968
00969
00970
00971
00972
00973
00974
00975
00976
00977
00978
00979
00980
00981
00982
00983
00984
00985
00986
00987 integer (kind=int_kind), intent(in) ::
00988 nx_block, ny_block ,
00989 ilo,ihi,jlo,jhi ,
00990 nghost
00991
00992 logical (kind=log_kind), dimension (ntrace), intent(in) ::
00993 has_dependents
00994
00995 integer (kind=int_kind), dimension(0:ncat), intent(out) ::
00996 icells
00997
00998 integer (kind=int_kind), dimension(nx_block*ny_block,0:ncat),
00999 intent(out) ::
01000 indxi ,
01001 indxj
01002
01003 real (kind=dbl_kind), dimension (nx_block,ny_block,0:ncat),
01004 intent(in) ::
01005 mm
01006
01007 real (kind=dbl_kind), dimension (nx_block,ny_block,0:ncat),
01008 intent(out) ::
01009 mmask
01010
01011 real (kind=dbl_kind), dimension (nx_block, ny_block, ntrace, ncat),
01012 intent(in), optional ::
01013 tm
01014
01015 real (kind=dbl_kind), dimension (nx_block, ny_block, ntrace, ncat),
01016 intent(out), optional ::
01017 tmask
01018
01019
01020
01021 integer (kind=int_kind) ::
01022 i, j, ij ,
01023 n ,
01024 nt
01025
01026 do n = 0, ncat
01027 do ij = 1, nx_block*ny_block
01028 indxi(ij,n) = 0
01029 indxj(ij,n) = 0
01030 enddo
01031 enddo
01032
01033
01034
01035
01036
01037 icells(0) = 0
01038 do j = 1, ny_block
01039 do i = 1, nx_block
01040 if (mm(i,j,0) > puny) then
01041 mmask(i,j,0) = c1
01042 icells(0) = icells(0) + 1
01043 ij = icells(0)
01044 indxi(ij,0) = i
01045 indxj(ij,0) = j
01046 else
01047 mmask(i,j,0) = c0
01048 endif
01049 enddo
01050 enddo
01051
01052 do n = 1, ncat
01053
01054
01055
01056
01057
01058 icells(n) = 0
01059 do j = 1, ny_block
01060 do i = 1, nx_block
01061 if (mm(i,j,n) > puny) then
01062 icells(n) = icells(n) + 1
01063 ij = icells(n)
01064 indxi(ij,n) = i
01065 indxj(ij,n) = j
01066 endif
01067 enddo
01068 enddo
01069
01070
01071
01072
01073
01074 mmask(:,:,n) = c0
01075 do ij = 1, icells(n)
01076 i = indxi(ij,n)
01077 j = indxj(ij,n)
01078 mmask(i,j,n) = c1
01079 enddo
01080
01081
01082
01083
01084
01085 if (present(tm)) then
01086
01087 tmask(:,:,:,n) = c0
01088 do nt = 1, ntrace
01089 if (has_dependents(nt)) then
01090 do ij = 1, icells(n)
01091 i = indxi(ij,n)
01092 j = indxj(ij,n)
01093 if (abs(tm(i,j,nt,n)) > puny) then
01094 tmask(i,j,nt,n) = c1
01095 endif
01096 enddo
01097 endif
01098 enddo
01099
01100 endif
01101
01102
01103
01104
01105
01106
01107
01108 icells(n) = 0
01109 do j = jlo-nghost+1, jhi+nghost-1
01110 do i = ilo-nghost+1, ihi+nghost-1
01111 if (mm(i,j,n) > puny) then
01112 icells(n) = icells(n) + 1
01113 ij = icells(n)
01114 indxi(ij,n) = i
01115 indxj(ij,n) = j
01116 endif
01117 enddo
01118 enddo
01119
01120 enddo
01121
01122 end subroutine make_masks
01123
01124
01125
01126
01127
01128
01129
01130
01131
01132 subroutine construct_fields (nx_block, ny_block, &
01133 ilo, ihi, jlo, jhi, &
01134 nghost, &
01135 tracer_type, depend, &
01136 has_dependents, icells, &
01137 indxi, indxj, &
01138 HTN, HTE, &
01139 hm, xav, &
01140 yav, xxav, &
01141 xyav, yyav, &
01142 xxxav, xxyav, &
01143 xyyav, yyyav, &
01144 dxt, dyt, &
01145 mm, mc, &
01146 mx, my, &
01147 mmask, &
01148 tm, tc, &
01149 tx, ty, &
01150 tmask)
01151
01152
01153
01154
01155
01156
01157
01158
01159
01160
01161
01162
01163
01164
01165 integer (kind=int_kind), intent(in) ::
01166 nx_block, ny_block ,
01167 ilo,ihi,jlo,jhi ,
01168 nghost ,
01169 icells
01170
01171 integer (kind=int_kind), dimension (ntrace), intent(in) ::
01172 tracer_type ,
01173 depend
01174
01175 logical (kind=log_kind), dimension (ntrace), intent(in) ::
01176 has_dependents
01177
01178 integer (kind=int_kind), dimension(nx_block*ny_block), intent(in) ::
01179 indxi ,
01180 indxj
01181
01182 real (kind=dbl_kind), dimension (nx_block,ny_block),
01183 intent(in) ::
01184 hm ,
01185 HTN ,
01186 HTE ,
01187 xav, yav ,
01188 xxav, xyav, yyav ,
01189 xxxav,xxyav,xyyav,yyyav,
01190 dxt ,
01191 dyt
01192
01193 real (kind=dbl_kind), dimension (nx_block,ny_block),
01194 intent(in) ::
01195 mm ,
01196 mmask
01197
01198 real (kind=dbl_kind), dimension (nx_block,ny_block,ntrace),
01199 intent(in), optional ::
01200 tm ,
01201 tmask
01202
01203 real (kind=dbl_kind), dimension (nx_block,ny_block),
01204 intent(out) ::
01205 mc ,
01206 mx, my
01207
01208 real (kind=dbl_kind), dimension (nx_block,ny_block,ntrace),
01209 intent(out), optional ::
01210 tc ,
01211 tx, ty
01212
01213
01214
01215 integer (kind=int_kind) ::
01216 i, j ,
01217 nt, nt1 ,
01218 ij
01219
01220 real (kind=dbl_kind), dimension (nx_block,ny_block) ::
01221 mxav ,
01222 myav
01223
01224 real (kind=dbl_kind), dimension (nx_block,ny_block,ntrace) ::
01225 mtxav ,
01226 mtyav
01227
01228 real (kind=dbl_kind) ::
01229 w1, w2, w3, w4, w5, w6, w7
01230
01231
01232
01233
01234
01235
01236
01237
01238
01239
01240
01241
01242
01243
01244
01245
01246
01247
01248
01249
01250
01251
01252
01253
01254
01255
01256
01257
01258
01259
01260
01261
01262
01263
01264
01265
01266
01267
01268
01269
01270
01271
01272
01273
01274 do j = 1, ny_block
01275 do i = 1, nx_block
01276 mc(i,j) = c0
01277 mx(i,j) = c0
01278 my(i,j) = c0
01279 mxav(i,j) = c0
01280 myav(i,j) = c0
01281 enddo
01282 enddo
01283
01284 if (present(tm)) then
01285 do nt = 1, ntrace
01286 do j = 1, ny_block
01287 do i = 1, nx_block
01288 tc(i,j,nt) = c0
01289 tx(i,j,nt) = c0
01290 ty(i,j,nt) = c0
01291 enddo
01292 enddo
01293 enddo
01294 endif
01295
01296
01297
01298
01299
01300 call limited_gradient (nx_block, ny_block, &
01301 ilo, ihi, jlo, jhi, &
01302 nghost, &
01303 mm, hm, &
01304 xav, yav, &
01305 HTN, HTE, &
01306 dxt, dyt, &
01307 mx, my)
01308
01309 do ij = 1,icells
01310 i = indxi(ij)
01311 j = indxj(ij)
01312
01313
01314 mc(i,j) = mm(i,j) - xav(i,j)*mx(i,j) &
01315 - yav(i,j)*my(i,j)
01316
01317 enddo
01318
01319
01320
01321 if (present(tm)) then
01322
01323 do ij = 1,icells
01324 i = indxi(ij)
01325 j = indxj(ij)
01326
01327
01328 mxav(i,j) = (mx(i,j)*xxav(i,j) &
01329 + my(i,j)*xyav(i,j) &
01330 + mc(i,j)*xav (i,j)) / mm(i,j)
01331 myav(i,j) = (mx(i,j)*xyav(i,j) &
01332 + my(i,j)*yyav(i,j) &
01333 + mc(i,j)*yav(i,j)) / mm(i,j)
01334 enddo
01335
01336 do nt = 1, ntrace
01337
01338 if (tracer_type(nt)==1) then
01339
01340 call limited_gradient(nx_block, ny_block, &
01341 ilo, ihi, jlo, jhi, &
01342 nghost, &
01343 tm(:,:,nt), mmask, &
01344 mxav, myav, &
01345 HTN, HTE, &
01346 dxt, dyt, &
01347 tx(:,:,nt), ty(:,:,nt))
01348
01349 if (has_dependents(nt)) then
01350
01351 do j = 1, ny_block
01352 do i = 1, nx_block
01353 mtxav(i,j,nt) = c0
01354 mtyav(i,j,nt) = c0
01355 enddo
01356 enddo
01357
01358
01359
01360
01361 do ij = 1, icells
01362
01363 i = indxi(ij)
01364 j = indxj(ij)
01365
01366
01367 tc(i,j,nt) = tm(i,j,nt) - tx(i,j,nt)*mxav(i,j) &
01368 - ty(i,j,nt)*myav(i,j)
01369
01370 if (tmask(i,j,nt) > puny) then
01371
01372
01373 w1 = mc(i,j)*tc(i,j,nt)
01374 w2 = mc(i,j)*tx(i,j,nt) &
01375 + mx(i,j)*tc(i,j,nt)
01376 w3 = mc(i,j)*ty(i,j,nt) &
01377 + my(i,j)*tc(i,j,nt)
01378 w4 = mx(i,j)*tx(i,j,nt)
01379 w5 = mx(i,j)*ty(i,j,nt) &
01380 + my(i,j)*tx(i,j,nt)
01381 w6 = my(i,j)*ty(i,j,nt)
01382 w7 = c1 / (mm(i,j)*tm(i,j,nt))
01383 mtxav(i,j,nt) = (w1*xav (i,j) + w2*xxav (i,j) &
01384 + w3*xyav (i,j) + w4*xxxav(i,j) &
01385 + w5*xxyav(i,j) + w6*xyyav(i,j)) &
01386 * w7
01387 mtyav(i,j,nt) = (w1*yav(i,j) + w2*xyav (i,j) &
01388 + w3*yyav(i,j) + w4*xxyav(i,j) &
01389 + w5*xyyav(i,j) + w6*yyyav(i,j)) &
01390 * w7
01391 endif
01392
01393 enddo
01394
01395 else
01396
01397 do ij = 1, icells
01398 i = indxi(ij)
01399 j = indxj(ij)
01400
01401
01402 tc(i,j,nt) = tm(i,j,nt) - tx(i,j,nt)*mxav(i,j) &
01403 - ty(i,j,nt)*myav(i,j)
01404 enddo
01405
01406 endif
01407
01408 elseif (tracer_type(nt)==2) then
01409 nt1 = depend(nt)
01410
01411 call limited_gradient(nx_block, ny_block, &
01412 ilo, ihi, jlo, jhi, &
01413 nghost, &
01414 tm(:,:,nt), tmask(:,:,nt1), &
01415 mtxav(:,:,nt1), mtyav(:,:,nt1), &
01416 HTN, HTE, &
01417 dxt, dyt, &
01418 tx(:,:,nt), ty(:,:,nt))
01419
01420 do ij = 1, icells
01421 i = indxi(ij)
01422 j = indxj(ij)
01423 tc(i,j,nt) = tm(i,j,nt) &
01424 - tx(i,j,nt) * mtxav(i,j,nt1) &
01425 - ty(i,j,nt) * mtyav(i,j,nt1)
01426 enddo
01427
01428 elseif (tracer_type(nt)==3) then
01429
01430 do ij = 1, icells
01431 i = indxi(ij)
01432 j = indxj(ij)
01433
01434 tc(i,j,nt) = tm(i,j,nt)
01435
01436
01437 enddo
01438
01439 endif
01440 enddo
01441
01442 endif
01443
01444 end subroutine construct_fields
01445
01446
01447
01448
01449
01450
01451
01452
01453
01454 subroutine limited_gradient (nx_block, ny_block, &
01455 ilo, ihi, jlo, jhi, &
01456 nghost, &
01457 phi, phimask, &
01458 cnx, cny, &
01459 HTN, HTE, &
01460 dxt, dyt, &
01461 gx, gy)
01462
01463
01464
01465
01466
01467
01468
01469
01470
01471
01472
01473
01474
01475
01476
01477
01478
01479
01480 integer (kind=int_kind), intent(in) ::
01481 nx_block, ny_block,
01482 ilo,ihi,jlo,jhi ,
01483 nghost
01484
01485 real (kind=dbl_kind), dimension (nx_block,ny_block),
01486 intent (in) ::
01487 phi ,
01488 cnx ,
01489 cny ,
01490 dxt ,
01491 dyt ,
01492 phimask ,
01493
01494
01495
01496 HTN ,&
01497 HTE
01498
01499 real (kind=dbl_kind), dimension (nx_block,ny_block),
01500 intent(out) ::
01501 gx ,
01502 gy
01503
01504
01505
01506 integer (kind=int_kind) ::
01507 i, j, ij ,
01508 icells
01509
01510 integer (kind=int_kind), dimension(nx_block*ny_block) ::
01511 indxi, indxj
01512
01513 real (kind=dbl_kind) ::
01514 phi_nw, phi_n, phi_ne ,
01515 phi_w, phi_e ,
01516 phi_sw, phi_s, phi_se ,
01517 qmn, qmx ,
01518 pmn, pmx ,
01519 w1, w2, w3, w4
01520
01521 real (kind=dbl_kind) ::
01522 gxtmp, gytmp
01523
01524 gx(:,:) = c0
01525 gy(:,:) = c0
01526
01527
01528
01529
01530 icells = 0
01531 do j = jlo-nghost+1, jhi+nghost-1
01532 do i = ilo-nghost+1, ihi+nghost-1
01533 if (phimask(i,j) > puny) then
01534
01535
01536
01537
01538
01539
01540
01541
01542
01543
01544
01545
01546
01547
01548
01549
01550 phi_nw = phimask(i-1,j+1) * phi(i-1,j+1) &
01551 + (c1-phimask(i-1,j+1))* phi(i,j)
01552 phi_n = phimask(i,j+1) * phi(i,j+1) &
01553 + (c1-phimask(i,j+1)) * phi(i,j)
01554 phi_ne = phimask(i+1,j+1) * phi(i+1,j+1) &
01555 + (c1-phimask(i+1,j+1))* phi(i,j)
01556 phi_w = phimask(i-1,j) * phi(i-1,j) &
01557 + (c1-phimask(i-1,j)) * phi(i,j)
01558 phi_e = phimask(i+1,j) * phi(i+1,j) &
01559 + (c1-phimask(i+1,j)) * phi(i,j)
01560 phi_sw = phimask(i-1,j-1) * phi(i-1,j-1) &
01561 + (c1-phimask(i-1,j-1))* phi(i,j)
01562 phi_s = phimask(i,j-1) * phi(i,j-1) &
01563 + (c1-phimask(i,j-1)) * phi(i,j)
01564 phi_se = phimask(i+1,j-1) * phi(i+1,j-1) &
01565 + (c1-phimask(i+1,j-1))* phi(i,j)
01566
01567
01568
01569
01570 gxtmp = (phi_e - phi(i,j)) / (dxt(i,j) + dxt(i+1,j)) &
01571 + (phi(i,j) - phi_w) / (dxt(i-1,j) + dxt(i,j) )
01572 gytmp = (phi_n - phi(i,j)) / (dyt(i,j) + dyt(i,j+1)) &
01573 + (phi(i,j) - phi_s) / (dyt(i,j-1) + dyt(i,j) )
01574
01575
01576 pmn = min (phi_nw, phi_n, phi_ne, phi_w, phi(i,j), &
01577 phi_e, phi_sw, phi_s, phi_se)
01578 pmx = max (phi_nw, phi_n, phi_ne, phi_w, phi(i,j), &
01579 phi_e, phi_sw, phi_s, phi_se)
01580
01581 pmn = pmn - phi(i,j)
01582 pmx = pmx - phi(i,j)
01583
01584
01585
01586 w1 = (p5*HTN(i,j) - cnx(i,j)) * gxtmp &
01587 + (p5*HTE(i,j) - cny(i,j)) * gytmp
01588 w2 = (p5*HTN(i,j-1) - cnx(i,j)) * gxtmp &
01589 - (p5*HTE(i,j) + cny(i,j)) * gytmp
01590 w3 = -(p5*HTN(i,j-1) + cnx(i,j)) * gxtmp &
01591 - (p5*HTE(i-1,j) + cny(i,j)) * gytmp
01592 w4 = (p5*HTE(i-1,j) - cny(i,j)) * gytmp &
01593 - (p5*HTN(i,j) + cnx(i,j)) * gxtmp
01594
01595 qmn = min (w1, w2, w3, w4)
01596 qmx = max (w1, w2, w3, w4)
01597
01598
01599
01600
01601 if (abs(qmn) > 1.0e-300_dbl_kind) then
01602 w1 = max(c0, pmn/qmn)
01603 else
01604 w1 = c1
01605 endif
01606
01607 if (abs(qmx) > 1.0e-300_dbl_kind) then
01608 w2 = max(c0, pmx/qmx)
01609 else
01610 w2 = c1
01611 endif
01612
01613 w1 = min(c1, w1, w2)
01614
01615
01616 gx(i,j) = w1 * gxtmp
01617 gy(i,j) = w1 * gytmp
01618
01619
01620 endif
01621 enddo
01622 enddo
01623
01624 end subroutine limited_gradient
01625
01626
01627
01628
01629
01630
01631
01632
01633 subroutine departure_points (nx_block, ny_block, &
01634 ilo, ihi, jlo, jhi, &
01635 nghost, dt, &
01636 uvel, vvel, &
01637 dxu, dyu, &
01638 HTN, HTE, &
01639 dpx, dpy, &
01640 l_dp_midpt, l_stop, &
01641 istop, jstop)
01642
01643
01644
01645
01646
01647
01648
01649
01650
01651
01652
01653
01654
01655
01656 integer (kind=int_kind), intent(in) ::
01657 nx_block, ny_block,
01658 ilo,ihi,jlo,jhi,
01659 nghost
01660
01661 real (kind=dbl_kind), intent(in) ::
01662 dt
01663
01664 real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) ::
01665 uvel ,
01666 vvel ,
01667 dxu ,
01668 dyu ,
01669 HTN ,
01670 HTE
01671
01672 real (kind=dbl_kind), dimension (nx_block,ny_block), intent(out) ::
01673 dpx ,
01674 dpy
01675
01676 logical (kind=log_kind), intent(in) ::
01677 l_dp_midpt
01678
01679
01680 logical (kind=log_kind), intent(inout) ::
01681 l_stop
01682
01683 integer (kind=int_kind), intent(inout) ::
01684 istop, jstop
01685
01686
01687
01688 integer (kind=int_kind) ::
01689 i, j, i2, j2
01690
01691 real (kind=dbl_kind) ::
01692 mpx, mpy ,
01693
01694 mpxt, mpyt ,&
01695 ump, vmp
01696
01697
01698
01699
01700
01701
01702
01703
01704
01705 dpx(:,:) = c0
01706 dpy(:,:) = c0
01707
01708 do j = jlo-nghost+1, jhi+nghost-1
01709 do i = ilo-nghost+1, ihi+nghost-1
01710
01711 dpx(i,j) = -dt*uvel(i,j)
01712 dpy(i,j) = -dt*vvel(i,j)
01713
01714
01715 if (dpx(i,j) < -HTN(i,j) .or. dpx(i,j) > HTN(i+1,j) .or. &
01716 dpy(i,j) < -HTE(i,j) .or. dpy(i,j) > HTE(i,j+1)) then
01717 l_stop = .true.
01718 istop = i
01719 jstop = j
01720 endif
01721
01722 enddo
01723 enddo
01724
01725 if (l_stop) then
01726 i = istop
01727 j = jstop
01728 write (nu_diag,*) ' '
01729 write (nu_diag,*) &
01730 'Warning: Departure points out of bounds in remap'
01731 write (nu_diag,*) 'my_task, i, j =', my_task, i, j
01732 write (nu_diag,*) 'dpx, dpy =', dpx(i,j), dpy(i,j)
01733 write (nu_diag,*) 'HTN(i,j), HTN(i+1,j) =', HTN(i,j), HTN(i+1,j)
01734 write (nu_diag,*) 'HTE(i,j), HTE(i,j+1) =', HTE(i,j), HTE(i,j+1)
01735 return
01736 endif
01737
01738 if (l_dp_midpt) then
01739
01740 do j = jlo-nghost+1, jhi+nghost-1
01741 do i = ilo-nghost+1, ihi+nghost-1
01742 if (uvel(i,j)/=c0 .or. vvel(i,j)/=c0) then
01743
01744
01745
01746
01747
01748
01749 dpx(i,j) = dpx(i,j) / dxu(i,j)
01750 dpy(i,j) = dpy(i,j) / dyu(i,j)
01751
01752
01753
01754
01755
01756 mpx = p5 * dpx(i,j)
01757 mpy = p5 * dpy(i,j)
01758
01759
01760
01761
01762
01763
01764
01765
01766 if (mpx >= c0 .and. mpy >= c0) then
01767 i2 = i+1
01768 j2 = j+1
01769 mpxt = mpx - p5
01770 mpyt = mpy - p5
01771 elseif (mpx < c0 .and. mpy < c0) then
01772 i2 = i
01773 j2 = j
01774 mpxt = mpx + p5
01775 mpyt = mpy + p5
01776 elseif (mpx >= c0 .and. mpy < c0) then
01777 i2 = i+1
01778 j2 = j
01779 mpxt = mpx - p5
01780 mpyt = mpy + p5
01781 elseif (mpx < c0 .and. mpy >= c0) then
01782 i2 = i
01783 j2 = j+1
01784 mpxt = mpx + p5
01785 mpyt = mpy - p5
01786 endif
01787
01788
01789
01790
01791
01792
01793 ump = uvel(i2-1,j2-1)*(mpxt-p5)*(mpyt-p5) &
01794 - uvel(i2, j2-1)*(mpxt+p5)*(mpyt-p5) &
01795 + uvel(i2, j2 )*(mpxt+p5)*(mpyt+p5) &
01796 - uvel(i2-1,j2 )*(mpxt-p5)*(mpyt+p5)
01797
01798 vmp = vvel(i2-1,j2-1)*(mpxt-p5)*(mpyt-p5) &
01799 - vvel(i2, j2-1)*(mpxt+p5)*(mpyt-p5) &
01800 + vvel(i2, j2 )*(mpxt+p5)*(mpyt+p5) &
01801 - vvel(i2-1,j2 )*(mpxt-p5)*(mpyt+p5)
01802
01803
01804
01805
01806
01807
01808 dpx(i,j) = -dt * ump
01809 dpy(i,j) = -dt * vmp
01810
01811 endif
01812
01813 enddo
01814 enddo
01815
01816 endif
01817
01818 end subroutine departure_points
01819
01820
01821
01822
01823
01824
01825
01826
01827
01828 subroutine locate_triangles (nx_block, ny_block, &
01829 ilo, ihi, jlo, jhi, &
01830 nghost, edge, &
01831 icells, &
01832 indxi, indxj, &
01833 dpx, dpy, &
01834 dxu, dyu, &
01835 xp, yp, &
01836 iflux, jflux, &
01837 triarea, &
01838 l_fixed_area, edgearea)
01839
01840
01841
01842
01843
01844
01845
01846
01847
01848
01849
01850
01851
01852
01853
01854
01855
01856 integer (kind=int_kind), intent(in) ::
01857 nx_block, ny_block,
01858 ilo,ihi,jlo,jhi ,
01859 nghost
01860
01861 character (len=char_len), intent(in) ::
01862 edge
01863
01864 real (kind=dbl_kind), dimension(nx_block,ny_block), intent(in) ::
01865 dpx ,
01866 dpy ,
01867 dxu ,
01868 dyu
01869
01870 real (kind=dbl_kind), dimension (nx_block,ny_block,0:nvert,ngroups),
01871 intent(out) ::
01872 xp, yp
01873
01874 real (kind=dbl_kind), dimension (nx_block,ny_block,ngroups),
01875 intent(out) ::
01876 triarea
01877
01878 integer (kind=int_kind), dimension (nx_block,ny_block,ngroups),
01879 intent(out) ::
01880 iflux ,
01881 jflux
01882
01883 integer (kind=int_kind), dimension (ngroups), intent(out) ::
01884 icells
01885
01886 integer (kind=int_kind), dimension (nx_block*ny_block,ngroups),
01887 intent(out) ::
01888 indxi ,
01889 indxj
01890
01891 logical, intent(in) ::
01892 l_fixed_area
01893
01894
01895
01896
01897 real (kind=dbl_kind), dimension(nx_block,ny_block), intent(inout) ::
01898 edgearea
01899
01900
01901
01902
01903 integer (kind=int_kind) ::
01904 i, j, ij, ic ,
01905 ib, ie, jb, je ,
01906 ng, nv ,
01907 ishift, jshift
01908
01909 integer (kind=int_kind) ::
01910 icellsd
01911
01912 integer (kind=int_kind), dimension (nx_block*ny_block) ::
01913 indxid ,
01914 indxjd
01915
01916 real (kind=dbl_kind), dimension(nx_block,ny_block) ::
01917 dx, dy ,
01918 areafac_c ,
01919 areafac_l ,
01920 areafac_r
01921
01922 real (kind=dbl_kind) ::
01923 xcl, ycl ,
01924
01925 xdl, ydl ,&
01926 xil, yil ,&
01927 xcr, ycr ,&
01928 xdr, ydr ,&
01929 xir, yir ,&
01930 xic, yic ,&
01931 xicl, yicl ,&
01932 xicr, yicr ,&
01933 xdm, ydm ,&
01934
01935 dxc ,&
01936 dxd ,&
01937 md ,&
01938 mdl ,&
01939 mdr ,&
01940 ishift_tl, jshift_tl ,&
01941 ishift_bl, jshift_bl ,&
01942 ishift_tr, jshift_tr ,&
01943 ishift_br, jshift_br ,&
01944 ishift_tc, jshift_tc ,&
01945 ishift_bc, jshift_bc ,&
01946 area1, area2 ,&
01947 area3, area4 ,&
01948 area_c ,&
01949 w1, w2
01950
01951 real (kind=dbl_kind), dimension (nx_block,ny_block,ngroups) ::
01952 areafact
01953
01954 real (kind=dbl_kind), dimension(nx_block,ny_block) ::
01955 areasum
01956
01957
01958
01959
01960
01961
01962
01963
01964
01965
01966
01967
01968
01969
01970
01971
01972
01973
01974
01975
01976
01977
01978
01979
01980
01981
01982
01983
01984
01985
01986
01987
01988
01989
01990
01991
01992
01993
01994
01995
01996
01997
01998
01999
02000
02001
02002
02003
02004
02005
02006
02007
02008
02009
02010 areafac_c(:,:) = c0
02011 areafac_l(:,:) = c0
02012 areafac_r(:,:) = c0
02013 do ng = 1, ngroups
02014 do j = 1, ny_block
02015 do i = 1, nx_block
02016 triarea (i,j,ng) = c0
02017 areafact(i,j,ng) = c0
02018 iflux (i,j,ng) = i
02019 jflux (i,j,ng) = j
02020 enddo
02021 enddo
02022 do nv = 0, nvert
02023 do j = 1, ny_block
02024 do i = 1, nx_block
02025 xp(i,j,nv,ng) = c0
02026 yp(i,j,nv,ng) = c0
02027 enddo
02028 enddo
02029 enddo
02030 enddo
02031
02032 if (trim(edge) == 'north') then
02033
02034
02035
02036 ib = ilo
02037 ie = ihi
02038 jb = jlo - nghost
02039 je = jhi
02040
02041
02042
02043 ishift_tl = -1
02044 jshift_tl = 1
02045 ishift_bl = -1
02046 jshift_bl = 0
02047 ishift_tr = 1
02048 jshift_tr = 1
02049 ishift_br = 1
02050 jshift_br = 0
02051 ishift_tc = 0
02052 jshift_tc = 1
02053 ishift_bc = 0
02054 jshift_bc = 0
02055
02056
02057
02058 do j = jb, je
02059 do i = ib, ie
02060 areafac_l(i,j) = dxu(i-1,j)*dyu(i-1,j)
02061 areafac_r(i,j) = dxu(i,j)*dyu(i,j)
02062 areafac_c(i,j) = p5*(areafac_l(i,j) + areafac_r(i,j))
02063 enddo
02064 enddo
02065
02066 else
02067
02068
02069
02070 ib = ilo - nghost
02071 ie = ihi
02072 jb = jlo
02073 je = jhi
02074
02075
02076
02077 ishift_tl = 1
02078 jshift_tl = 1
02079 ishift_bl = 0
02080 jshift_bl = 1
02081 ishift_tr = 1
02082 jshift_tr = -1
02083 ishift_br = 0
02084 jshift_br = -1
02085 ishift_tc = 1
02086 jshift_tc = 0
02087 ishift_bc = 0
02088 jshift_bc = 0
02089
02090
02091
02092 do j = jb, je
02093 do i = ib, ie
02094 areafac_l(i,j) = dxu(i,j)*dyu(i,j)
02095 areafac_r(i,j) = dxu(i,j-1)*dyu(i,j-1)
02096 areafac_c(i,j) = p5 * (areafac_l(i,j) + areafac_r(i,j))
02097 enddo
02098 enddo
02099
02100 endif
02101
02102
02103
02104
02105
02106 if (l_fixed_area) then
02107 icellsd = 0
02108 do j = jb, je
02109 do i = ib, ie
02110 if (edgearea(i,j) /= c0) then
02111 icellsd = icellsd + 1
02112 indxid(icellsd) = i
02113 indxjd(icellsd) = j
02114 endif
02115 enddo
02116 enddo
02117 else
02118 icellsd = 0
02119 if (trim(edge) == 'north') then
02120 do j = jb, je
02121 do i = ib, ie
02122 if (dpx(i-1,j)/=c0 .or. dpy(i-1,j)/=c0 &
02123 .or. &
02124 dpx(i,j)/=c0 .or. dpy(i,j)/=c0) then
02125 icellsd = icellsd + 1
02126 indxid(icellsd) = i
02127 indxjd(icellsd) = j
02128 endif
02129 enddo
02130 enddo
02131 else
02132 do j = jb, je
02133 do i = ib, ie
02134 if (dpx(i,j-1)/=c0 .or. dpy(i,j-1)/=c0 &
02135 .or. &
02136 dpx(i,j)/=c0 .or. dpy(i,j)/=c0) then
02137 icellsd = icellsd + 1
02138 indxid(icellsd) = i
02139 indxjd(icellsd) = j
02140 endif
02141 enddo
02142 enddo
02143 endif
02144 endif
02145
02146
02147
02148
02149
02150 do j = 1, je
02151 do i = 1, ie
02152 dx(i,j) = dpx(i,j) / dxu(i,j)
02153 dy(i,j) = dpy(i,j) / dyu(i,j)
02154 enddo
02155 enddo
02156
02157
02158
02159
02160
02161
02162
02163
02164
02165
02166
02167 do ij = 1, icellsd
02168 i = indxid(ij)
02169 j = indxjd(ij)
02170
02171 xcl = -p5
02172 ycl = c0
02173
02174 xcr = p5
02175 ycr = c0
02176
02177
02178
02179 if (trim(edge) == 'north') then
02180 xdl = xcl + dx(i-1,j)
02181 ydl = ycl + dy(i-1,j)
02182 xdr = xcr + dx(i,j)
02183 ydr = ycr + dy(i,j)
02184 else
02185 xdl = xcl - dy(i,j)
02186 ydl = ycl + dx(i,j)
02187 xdr = xcr - dy(i,j-1)
02188 ydr = ycr + dx(i,j-1)
02189 endif
02190
02191 xdm = p5 * (xdr + xdl)
02192 ydm = p5 * (ydr + ydl)
02193
02194
02195
02196 xil = xcl
02197 yil = (xcl*(ydm-ydl) + xdm*ydl - xdl*ydm) / (xdm - xdl)
02198
02199 xir = xcr
02200 yir = (xcr*(ydr-ydm) - xdm*ydr + xdr*ydm) / (xdr - xdm)
02201
02202 md = (ydr - ydl) / (xdr - xdl)
02203
02204 if (abs(md) > puny) then
02205 xic = xdl - ydl/md
02206 else
02207 xic = c0
02208 endif
02209 yic = c0
02210
02211 xicl = xic
02212 yicl = yic
02213 xicr = xic
02214 yicr = yic
02215
02216
02217
02218
02219
02220
02221 if (yil > c0 .and. xdl < xcl .and. ydl >= c0) then
02222
02223
02224
02225 ng = 1
02226 xp (i,j,1,ng) = xcl
02227 yp (i,j,1,ng) = ycl
02228 xp (i,j,2,ng) = xil
02229 yp (i,j,2,ng) = yil
02230 xp (i,j,3,ng) = xdl
02231 yp (i,j,3,ng) = ydl
02232 iflux (i,j,ng) = i + ishift_tl
02233 jflux (i,j,ng) = j + jshift_tl
02234 areafact(i,j,ng) = -areafac_l(i,j)
02235
02236 elseif (yil < c0 .and. xdl < xcl .and. ydl < c0) then
02237
02238
02239
02240 ng = 1
02241 xp (i,j,1,ng) = xcl
02242 yp (i,j,1,ng) = ycl
02243 xp (i,j,2,ng) = xdl
02244 yp (i,j,2,ng) = ydl
02245 xp (i,j,3,ng) = xil
02246 yp (i,j,3,ng) = yil
02247 iflux (i,j,ng) = i + ishift_bl
02248 jflux (i,j,ng) = j + jshift_bl
02249 areafact(i,j,ng) = areafac_l(i,j)
02250
02251 elseif (yil < c0 .and. xdl < xcl .and. ydl >= c0) then
02252
02253
02254
02255 ng = 1
02256 xp (i,j,1,ng) = xcl
02257 yp (i,j,1,ng) = ycl
02258 xp (i,j,2,ng) = xdl
02259 yp (i,j,2,ng) = ydl
02260 xp (i,j,3,ng) = xic
02261 yp (i,j,3,ng) = yic
02262 iflux (i,j,ng) = i + ishift_tl
02263 jflux (i,j,ng) = j + jshift_tl
02264 areafact(i,j,ng) = areafac_l(i,j)
02265
02266
02267
02268 ng = 3
02269 xp (i,j,1,ng) = xcl
02270 yp (i,j,1,ng) = ycl
02271 xp (i,j,2,ng) = xic
02272 yp (i,j,2,ng) = yic
02273 xp (i,j,3,ng) = xil
02274 yp (i,j,3,ng) = yil
02275 iflux (i,j,ng) = i + ishift_bl
02276 jflux (i,j,ng) = j + jshift_bl
02277 areafact(i,j,ng) = areafac_l(i,j)
02278
02279 elseif (yil > c0 .and. xdl < xcl .and. ydl < c0) then
02280
02281
02282
02283 ng = 3
02284 xp (i,j,1,ng) = xcl
02285 yp (i,j,1,ng) = ycl
02286 xp (i,j,2,ng) = xil
02287 yp (i,j,2,ng) = yil
02288 xp (i,j,3,ng) = xic
02289 yp (i,j,3,ng) = yic
02290 iflux (i,j,ng) = i + ishift_tl
02291 jflux (i,j,ng) = j + jshift_tl
02292 areafact(i,j,ng) = -areafac_l(i,j)
02293
02294
02295
02296 ng = 1
02297 xp (i,j,1,ng) = xcl
02298 yp (i,j,1,ng) = ycl
02299 xp (i,j,2,ng) = xic
02300 yp (i,j,2,ng) = yic
02301 xp (i,j,3,ng) = xdl
02302 yp (i,j,3,ng) = ydl
02303 iflux (i,j,ng) = i + ishift_bl
02304 jflux (i,j,ng) = j + jshift_bl
02305 areafact(i,j,ng) = -areafac_l(i,j)
02306
02307 endif
02308
02309
02310
02311
02312
02313
02314 if (yir > c0 .and. xdr >= xcr .and. ydr >= c0) then
02315
02316
02317
02318 ng = 2
02319 xp (i,j,1,ng) = xcr
02320 yp (i,j,1,ng) = ycr
02321 xp (i,j,2,ng) = xdr
02322 yp (i,j,2,ng) = ydr
02323 xp (i,j,3,ng) = xir
02324 yp (i,j,3,ng) = yir
02325 iflux (i,j,ng) = i + ishift_tr
02326 jflux (i,j,ng) = j + jshift_tr
02327 areafact(i,j,ng) = -areafac_r(i,j)
02328
02329 elseif (yir < c0 .and. xdr >= xcr .and. ydr < c0) then
02330
02331
02332
02333 ng = 2
02334 xp (i,j,1,ng) = xcr
02335 yp (i,j,1,ng) = ycr
02336 xp (i,j,2,ng) = xir
02337 yp (i,j,2,ng) = yir
02338 xp (i,j,3,ng) = xdr
02339 yp (i,j,3,ng) = ydr
02340 iflux (i,j,ng) = i + ishift_br
02341 jflux (i,j,ng) = j + jshift_br
02342 areafact(i,j,ng) = areafac_r(i,j)
02343
02344 elseif (yir < c0 .and. xdr >= xcr .and. ydr >= c0) then
02345
02346
02347
02348 ng = 2
02349 xp (i,j,1,ng) = xcr
02350 yp (i,j,1,ng) = ycr
02351 xp (i,j,2,ng) = xic
02352 yp (i,j,2,ng) = yic
02353 xp (i,j,3,ng) = xdr
02354 yp (i,j,3,ng) = ydr
02355 iflux (i,j,ng) = i + ishift_tr
02356 jflux (i,j,ng) = j + jshift_tr
02357 areafact(i,j,ng) = areafac_r(i,j)
02358
02359
02360
02361 ng = 3
02362 xp (i,j,1,ng) = xcr
02363 yp (i,j,1,ng) = ycr
02364 xp (i,j,2,ng) = xir
02365 yp (i,j,2,ng) = yir
02366 xp (i,j,3,ng) = xic
02367 yp (i,j,3,ng) = yic
02368 iflux (i,j,ng) = i + ishift_br
02369 jflux (i,j,ng) = j + jshift_br
02370 areafact(i,j,ng) = areafac_r(i,j)
02371
02372 elseif (yir > c0 .and. xdr >= xcr .and. ydr < c0) then
02373
02374
02375
02376 ng = 3
02377 xp (i,j,1,ng) = xcr
02378 yp (i,j,1,ng) = ycr
02379 xp (i,j,2,ng) = xic
02380 yp (i,j,2,ng) = yic
02381 xp (i,j,3,ng) = xir
02382 yp (i,j,3,ng) = yir
02383 iflux (i,j,ng) = i + ishift_tr
02384 jflux (i,j,ng) = j + jshift_tr
02385 areafact(i,j,ng) = -areafac_r(i,j)
02386
02387
02388
02389 ng = 2
02390 xp (i,j,1,ng) = xcr
02391 yp (i,j,1,ng) = ycr
02392 xp (i,j,2,ng) = xdr
02393 yp (i,j,2,ng) = ydr
02394 xp (i,j,3,ng) = xic
02395 yp (i,j,3,ng) = yic
02396 iflux (i,j,ng) = i + ishift_br
02397 jflux (i,j,ng) = j + jshift_br
02398 areafact(i,j,ng) = -areafac_r(i,j)
02399
02400 endif
02401
02402
02403
02404
02405
02406 if (xdl < xcl) then
02407 xdl = xil
02408 ydl = yil
02409 endif
02410
02411 if (xdr > xcr) then
02412 xdr = xir
02413 ydr = yir
02414 endif
02415
02416
02417
02418
02419
02420
02421 if (l_fixed_area) then
02422
02423
02424
02425
02426
02427 ng = 1
02428 area1 = p5 * ( (xp(i,j,2,ng)-xp(i,j,1,ng)) * &
02429 yp(i,j,3,ng) &
02430 - yp(i,j,2,ng) * &
02431 (xp(i,j,3,ng)-xp(i,j,1,ng)) ) &
02432 * areafact(i,j,ng)
02433
02434 ng = 2
02435 area2 = p5 * ( (xp(i,j,2,ng)-xp(i,j,1,ng)) * &
02436 yp(i,j,3,ng) &
02437 - yp(i,j,2,ng) * &
02438 (xp(i,j,3,ng)-xp(i,j,1,ng)) ) &
02439 * areafact(i,j,ng)
02440
02441 ng = 3
02442 area3 = p5 * ( (xp(i,j,2,ng)-xp(i,j,1,ng)) * &
02443 yp(i,j,3,ng) &
02444 - yp(i,j,2,ng) * &
02445 (xp(i,j,3,ng)-xp(i,j,1,ng)) ) &
02446 * areafact(i,j,ng)
02447
02448
02449
02450
02451
02452
02453
02454
02455
02456
02457
02458
02459
02460
02461 if (ydl*ydr >= c0) then
02462
02463
02464 area_c = edgearea(i,j) - area1 - area2 - area3
02465
02466
02467 w1 = c2*area_c/areafac_c(i,j) &
02468 + (xdr-xcl)*ydl + (xcr-xdl)*ydr
02469 w2 = (xdr-xdl)**2 + (ydr-ydl)**2
02470 w1 = w1/w2
02471 xdm = xdm + (ydr - ydl) * w1
02472 ydm = ydm - (xdr - xdl) * w1
02473
02474
02475 mdl = (ydm - ydl) / (xdm - xdl)
02476 mdr = (ydr - ydm) / (xdr - xdm)
02477
02478 if (abs(mdl) > puny) then
02479 xicl = xdl - ydl/mdl
02480 else
02481 xicl = c0
02482 endif
02483 yicl = c0
02484
02485 if (abs(mdr) > puny) then
02486 xicr = xdr - ydr/mdr
02487 else
02488 xicr = c0
02489 endif
02490 yicr = c0
02491
02492 elseif (xic < c0) then
02493
02494 xicl = xic
02495 yicl = yic
02496
02497
02498 xdm = p5 * (xdr + xicl)
02499 ydm = p5 * ydr
02500
02501
02502 area4 = p5 * (xcl - xic) * ydl * areafac_l(i,j)
02503 area_c = edgearea(i,j) - area1 - area2 - area3 - area4
02504
02505
02506 w1 = c2*area_c/areafac_c(i,j) + (xcr-xic)*ydr
02507 w2 = (xdr-xic)**2 + ydr**2
02508 w1 = w1/w2
02509 xdm = xdm + ydr*w1
02510 ydm = ydm - (xdr - xic) * w1
02511
02512
02513 mdr = (ydr - ydm) / (xdr - xdm)
02514 if (abs(mdr) > puny) then
02515 xicr = xdr - ydr/mdr
02516 else
02517 xicr = c0
02518 endif
02519 yicr = c0
02520
02521 elseif (xic >= c0) then
02522
02523 xicr = xic
02524 yicr = yic
02525
02526
02527 xdm = p5 * (xicr + xdl)
02528 ydm = p5 * ydl
02529
02530 area4 = p5 * (xic - xcr) * ydr * areafac_r(i,j)
02531 area_c = edgearea(i,j) - area1 - area2 - area3 - area4
02532
02533
02534 w1 = c2*area_c/areafac_c(i,j) + (xic-xcl)*ydl
02535 w2 = (xic-xdl)**2 + ydl**2
02536 w1 = w1/w2
02537 xdm = xdm - ydl*w1
02538 ydm = ydm - (xic - xdl) * w1
02539
02540
02541
02542 mdl = (ydm - ydl) / (xdm - xdl)
02543 if (abs(mdl) > puny) then
02544 xicl = xdl - ydl/mdl
02545 else
02546 xicl = c0
02547 endif
02548 yicl = c0
02549
02550 endif
02551
02552 endif
02553
02554
02555
02556
02557
02558
02559
02560
02561 if (ydl >= c0 .and. ydr >= c0 .and. ydm >= c0) then
02562
02563
02564
02565 ng = 4
02566 xp (i,j,1,ng) = xcl
02567 yp (i,j,1,ng) = ycl
02568 xp (i,j,2,ng) = xcr
02569 yp (i,j,2,ng) = ycr
02570 xp (i,j,3,ng) = xdl
02571 yp (i,j,3,ng) = ydl
02572 iflux (i,j,ng) = i + ishift_tc
02573 jflux (i,j,ng) = j + jshift_tc
02574 areafact(i,j,ng) = -areafac_c(i,j)
02575
02576
02577
02578 ng = 5
02579 xp (i,j,1,ng) = xcr
02580 yp (i,j,1,ng) = ycr
02581 xp (i,j,2,ng) = xdr
02582 yp (i,j,2,ng) = ydr
02583 xp (i,j,3,ng) = xdl
02584 yp (i,j,3,ng) = ydl
02585 iflux (i,j,ng) = i + ishift_tc
02586 jflux (i,j,ng) = j + jshift_tc
02587 areafact(i,j,ng) = -areafac_c(i,j)
02588
02589
02590 ng = 6
02591 xp (i,j,1,ng) = xdl
02592 yp (i,j,1,ng) = ydl
02593 xp (i,j,2,ng) = xdr
02594 yp (i,j,2,ng) = ydr
02595 xp (i,j,3,ng) = xdm
02596 yp (i,j,3,ng) = ydm
02597 iflux (i,j,ng) = i + ishift_tc
02598 jflux (i,j,ng) = j + jshift_tc
02599 areafact(i,j,ng) = -areafac_c(i,j)
02600
02601 elseif (ydl >= c0 .and. ydr >= c0 .and. ydm < c0) then
02602
02603
02604
02605 ng = 4
02606 xp (i,j,1,ng) = xcl
02607 yp (i,j,1,ng) = ycl
02608 xp (i,j,2,ng) = xicl
02609 yp (i,j,2,ng) = yicl
02610 xp (i,j,3,ng) = xdl
02611 yp (i,j,3,ng) = ydl
02612 iflux (i,j,ng) = i + ishift_tc
02613 jflux (i,j,ng) = j + jshift_tc
02614 areafact(i,j,ng) = -areafac_c(i,j)
02615
02616
02617
02618 ng = 5
02619 xp (i,j,1,ng) = xcr
02620 yp (i,j,1,ng) = ycr
02621 xp (i,j,2,ng) = xdr
02622 yp (i,j,2,ng) = ydr
02623 xp (i,j,3,ng) = xicr
02624 yp (i,j,3,ng) = yicr
02625 iflux (i,j,ng) = i + ishift_tc
02626 jflux (i,j,ng) = j + jshift_tc
02627 areafact(i,j,ng) = -areafac_c(i,j)
02628
02629
02630
02631 ng = 6
02632 xp (i,j,1,ng) = xicr
02633 yp (i,j,1,ng) = yicr
02634 xp (i,j,2,ng) = xicl
02635 yp (i,j,2,ng) = yicl
02636 xp (i,j,3,ng) = xdm
02637 yp (i,j,3,ng) = ydm
02638 iflux (i,j,ng) = i + ishift_bc
02639 jflux (i,j,ng) = j + jshift_bc
02640 areafact(i,j,ng) = areafac_c(i,j)
02641
02642 elseif (ydl < c0 .and. ydr < c0 .and. ydm < c0) then
02643
02644
02645
02646 ng = 4
02647 xp (i,j,1,ng) = xcl
02648 yp (i,j,1,ng) = ycl
02649 xp (i,j,2,ng) = xdl
02650 yp (i,j,2,ng) = ydl
02651 xp (i,j,3,ng) = xcr
02652 yp (i,j,3,ng) = ycr
02653 iflux (i,j,ng) = i + ishift_bc
02654 jflux (i,j,ng) = j + jshift_bc
02655 areafact(i,j,ng) = areafac_c(i,j)
02656
02657
02658
02659 ng = 5
02660 xp (i,j,1,ng) = xcr
02661 yp (i,j,1,ng) = ycr
02662 xp (i,j,2,ng) = xdl
02663 yp (i,j,2,ng) = ydl
02664 xp (i,j,3,ng) = xdr
02665 yp (i,j,3,ng) = ydr
02666 iflux (i,j,ng) = i + ishift_bc
02667 jflux (i,j,ng) = j + jshift_bc
02668 areafact(i,j,ng) = areafac_c(i,j)
02669
02670
02671
02672 ng = 6
02673 xp (i,j,1,ng) = xdl
02674 yp (i,j,1,ng) = ydl
02675 xp (i,j,2,ng) = xdm
02676 yp (i,j,2,ng) = ydm
02677 xp (i,j,3,ng) = xdr
02678 yp (i,j,3,ng) = ydr
02679 iflux (i,j,ng) = i + ishift_bc
02680 jflux (i,j,ng) = j + jshift_bc
02681 areafact(i,j,ng) = areafac_c(i,j)
02682
02683 elseif (ydl < c0 .and. ydr < c0 .and. ydm >= c0) then
02684
02685
02686
02687 ng = 4
02688 xp (i,j,1,ng) = xcl
02689 yp (i,j,1,ng) = ycl
02690 xp (i,j,2,ng) = xdl
02691 yp (i,j,2,ng) = ydl
02692 xp (i,j,3,ng) = xicl
02693 yp (i,j,3,ng) = yicl
02694 iflux (i,j,ng) = i + ishift_bc
02695 jflux (i,j,ng) = j + jshift_bc
02696 areafact(i,j,ng) = areafac_c(i,j)
02697
02698
02699
02700 ng = 5
02701 xp (i,j,1,ng) = xcr
02702 yp (i,j,1,ng) = ycr
02703 xp (i,j,2,ng) = xicr
02704 yp (i,j,2,ng) = yicr
02705 xp (i,j,3,ng) = xdr
02706 yp (i,j,3,ng) = ydr
02707 iflux (i,j,ng) = i + ishift_bc
02708 jflux (i,j,ng) = j + jshift_bc
02709 areafact(i,j,ng) = areafac_c(i,j)
02710
02711
02712
02713 ng = 6
02714 xp (i,j,1,ng) = xicl
02715 yp (i,j,1,ng) = yicl
02716 xp (i,j,2,ng) = xicr
02717 yp (i,j,2,ng) = yicr
02718 xp (i,j,3,ng) = xdm
02719 yp (i,j,3,ng) = ydm
02720 iflux (i,j,ng) = i + ishift_tc
02721 jflux (i,j,ng) = j + jshift_tc
02722 areafact(i,j,ng) = -areafac_c(i,j)
02723
02724
02725
02726
02727
02728 elseif (ydl >= c0 .and. ydr < c0 .and. xic >= c0 &
02729 .and. ydm >= c0) then
02730
02731
02732
02733 ng = 4
02734 xp (i,j,1,ng) = xcl
02735 yp (i,j,1,ng) = ycl
02736 xp (i,j,2,ng) = xicr
02737 yp (i,j,2,ng) = yicr
02738 xp (i,j,3,ng) = xdl
02739 yp (i,j,3,ng) = ydl
02740 iflux (i,j,ng) = i + ishift_tc
02741 jflux (i,j,ng) = j + jshift_tc
02742 areafact(i,j,ng) = -areafac_c(i,j)
02743
02744
02745
02746 ng = 5
02747 xp (i,j,1,ng) = xcr
02748 yp (i,j,1,ng) = ycr
02749 xp (i,j,2,ng) = xicr
02750 yp (i,j,2,ng) = yicr
02751 xp (i,j,3,ng) = xdr
02752 yp (i,j,3,ng) = ydr
02753 iflux (i,j,ng) = i + ishift_bc
02754 jflux (i,j,ng) = j + jshift_bc
02755 areafact(i,j,ng) = areafac_r(i,j)
02756
02757
02758
02759 ng = 6
02760 xp (i,j,1,ng) = xdl
02761 yp (i,j,1,ng) = ydl
02762 xp (i,j,2,ng) = xicr
02763 yp (i,j,2,ng) = yicr
02764 xp (i,j,3,ng) = xdm
02765 yp (i,j,3,ng) = ydm
02766 iflux (i,j,ng) = i + ishift_tc
02767 jflux (i,j,ng) = j + jshift_tc
02768 areafact(i,j,ng) = -areafac_c(i,j)
02769
02770 elseif (ydl >= c0 .and. ydr < c0 .and. xic >= c0 &
02771 .and. ydm < c0 ) then
02772
02773
02774
02775 ng = 4
02776 xp (i,j,1,ng) = xcl
02777 yp (i,j,1,ng) = ycl
02778 xp (i,j,2,ng) = xicl
02779 yp (i,j,2,ng) = yicl
02780 xp (i,j,3,ng) = xdl
02781 yp (i,j,3,ng) = ydl
02782 iflux (i,j,ng) = i + ishift_tc
02783 jflux (i,j,ng) = j + jshift_tc
02784 areafact(i,j,ng) = -areafac_c(i,j)
02785
02786
02787
02788 ng = 5
02789 xp (i,j,1,ng) = xcr
02790 yp (i,j,1,ng) = ycr
02791 xp (i,j,2,ng) = xicr
02792 yp (i,j,2,ng) = yicr
02793 xp (i,j,3,ng) = xdr
02794 yp (i,j,3,ng) = ydr
02795 iflux (i,j,ng) = i + ishift_bc
02796 jflux (i,j,ng) = j + jshift_bc
02797 areafact(i,j,ng) = areafac_r(i,j)
02798
02799
02800
02801 ng = 6
02802 xp (i,j,1,ng) = xicr
02803 yp (i,j,1,ng) = yicr
02804 xp (i,j,2,ng) = xicl
02805 yp (i,j,2,ng) = yicl
02806 xp (i,j,3,ng) = xdm
02807 yp (i,j,3,ng) = ydm
02808 iflux (i,j,ng) = i + ishift_bc
02809 jflux (i,j,ng) = j + jshift_bc
02810 areafact(i,j,ng) = areafac_c(i,j)
02811
02812 elseif (ydl >= c0 .and. ydr < c0 .and. xic < c0 &
02813 .and. ydm < c0) then
02814
02815
02816
02817 ng = 4
02818 xp (i,j,1,ng) = xcl
02819 yp (i,j,1,ng) = ycl
02820 xp (i,j,2,ng) = xicl
02821 yp (i,j,2,ng) = yicl
02822 xp (i,j,3,ng) = xdl
02823 yp (i,j,3,ng) = ydl
02824 iflux (i,j,ng) = i + ishift_tc
02825 jflux (i,j,ng) = j + jshift_tc
02826 areafact(i,j,ng) = -areafac_l(i,j)
02827
02828
02829
02830 ng = 5
02831 xp (i,j,1,ng) = xcr
02832 yp (i,j,1,ng) = ycr
02833 xp (i,j,2,ng) = xicl
02834 yp (i,j,2,ng) = yicl
02835 xp (i,j,3,ng) = xdr
02836 yp (i,j,3,ng) = ydr
02837 iflux (i,j,ng) = i + ishift_bc
02838 jflux (i,j,ng) = j + jshift_bc
02839 areafact(i,j,ng) = areafac_c(i,j)
02840
02841
02842
02843 ng = 6
02844 xp (i,j,1,ng) = xdr
02845 yp (i,j,1,ng) = ydr
02846 xp (i,j,2,ng) = xicl
02847 yp (i,j,2,ng) = yicl
02848 xp (i,j,3,ng) = xdm
02849 yp (i,j,3,ng) = ydm
02850 iflux (i,j,ng) = i + ishift_bc
02851 jflux (i,j,ng) = j + jshift_bc
02852 areafact(i,j,ng) = areafac_c(i,j)
02853
02854 elseif (ydl >= c0 .and. ydr < c0 .and. xic < c0 &
02855 .and. ydm >= c0) then
02856
02857
02858
02859 ng = 4
02860 xp (i,j,1,ng) = xcl
02861 yp (i,j,1,ng) = ycl
02862 xp (i,j,2,ng) = xicl
02863 yp (i,j,2,ng) = yicl
02864 xp (i,j,3,ng) = xdl
02865 yp (i,j,3,ng) = ydl
02866 iflux (i,j,ng) = i + ishift_tc
02867 jflux (i,j,ng) = j + jshift_tc
02868 areafact(i,j,ng) = -areafac_l(i,j)
02869
02870
02871
02872 ng = 5
02873 xp (i,j,1,ng) = xcr
02874 yp (i,j,1,ng) = ycr
02875 xp (i,j,2,ng) = xicr
02876 yp (i,j,2,ng) = yicr
02877 xp (i,j,3,ng) = xdr
02878 yp (i,j,3,ng) = ydr
02879 iflux (i,j,ng) = i + ishift_bc
02880 jflux (i,j,ng) = j + jshift_bc
02881 areafact(i,j,ng) = areafac_c(i,j)
02882
02883
02884
02885 ng = 6
02886 xp (i,j,1,ng) = xicl
02887 yp (i,j,1,ng) = yicl
02888 xp (i,j,2,ng) = xicr
02889 yp (i,j,2,ng) = yicr
02890 xp (i,j,3,ng) = xdm
02891 yp (i,j,3,ng) = ydm
02892 iflux (i,j,ng) = i + ishift_tc
02893 jflux (i,j,ng) = j + jshift_tc
02894 areafact(i,j,ng) = -areafac_c(i,j)
02895
02896 elseif (ydl < c0 .and. ydr >= c0 .and. xic < c0 &
02897 .and. ydm >= c0) then
02898
02899
02900
02901 ng = 4
02902 xp (i,j,1,ng) = xcl
02903 yp (i,j,1,ng) = ycl
02904 xp (i,j,2,ng) = xdl
02905 yp (i,j,2,ng) = ydl
02906 xp (i,j,3,ng) = xicl
02907 yp (i,j,3,ng) = yicl
02908 iflux (i,j,ng) = i + ishift_bc
02909 jflux (i,j,ng) = j + jshift_bc
02910 areafact(i,j,ng) = areafac_l(i,j)
02911
02912
02913
02914 ng = 5
02915 xp (i,j,1,ng) = xcr
02916 yp (i,j,1,ng) = ycr
02917 xp (i,j,2,ng) = xdr
02918 yp (i,j,2,ng) = ydr
02919 xp (i,j,3,ng) = xicl
02920 yp (i,j,3,ng) = yicl
02921 iflux (i,j,ng) = i + ishift_tc
02922 jflux (i,j,ng) = j + jshift_tc
02923 areafact(i,j,ng) = -areafac_c(i,j)
02924
02925
02926
02927 ng = 6
02928 xp (i,j,1,ng) = xicl
02929 yp (i,j,1,ng) = yicl
02930 xp (i,j,2,ng) = xdr
02931 yp (i,j,2,ng) = ydr
02932 xp (i,j,3,ng) = xdm
02933 yp (i,j,3,ng) = ydm
02934 iflux (i,j,ng) = i + ishift_tc
02935 jflux (i,j,ng) = j + jshift_tc
02936 areafact(i,j,ng) = -areafac_c(i,j)
02937
02938 elseif (ydl < c0 .and. ydr >= c0 .and. xic < c0 &
02939 .and. ydm < c0) then
02940
02941
02942
02943 ng = 4
02944 xp (i,j,1,ng) = xcl
02945 yp (i,j,1,ng) = ycl
02946 xp (i,j,2,ng) = xdl
02947 yp (i,j,2,ng) = ydl
02948 xp (i,j,3,ng) = xicl
02949 yp (i,j,3,ng) = yicl
02950 iflux (i,j,ng) = i + ishift_bc
02951 jflux (i,j,ng) = j + jshift_bc
02952 areafact(i,j,ng) = areafac_l(i,j)
02953
02954
02955
02956 ng = 5
02957 xp (i,j,1,ng) = xcr
02958 yp (i,j,1,ng) = ycr
02959 xp (i,j,2,ng) = xdr
02960 yp (i,j,2,ng) = ydr
02961 xp (i,j,3,ng) = xicr
02962 yp (i,j,3,ng) = yicr
02963 iflux (i,j,ng) = i + ishift_tc
02964 jflux (i,j,ng) = j + jshift_tc
02965 areafact(i,j,ng) = -areafac_c(i,j)
02966
02967
02968
02969 ng = 6
02970 xp (i,j,1,ng) = xicr
02971 yp (i,j,1,ng) = yicr
02972 xp (i,j,2,ng) = xicl
02973 yp (i,j,2,ng) = yicl
02974 xp (i,j,3,ng) = xdm
02975 yp (i,j,3,ng) = ydm
02976 iflux (i,j,ng) = i + ishift_bc
02977 jflux (i,j,ng) = j + jshift_bc
02978 areafact(i,j,ng) = areafac_c(i,j)
02979
02980 elseif (ydl < c0 .and. ydr >= c0 .and. xic >= c0 &
02981 .and. ydm < c0) then
02982
02983
02984
02985 ng = 4
02986 xp (i,j,1,ng) = xcl
02987 yp (i,j,1,ng) = ycl
02988 xp (i,j,2,ng) = xdl
02989 yp (i,j,2,ng) = ydl
02990 xp (i,j,3,ng) = xicr
02991 yp (i,j,3,ng) = yicr
02992 iflux (i,j,ng) = i + ishift_bc
02993 jflux (i,j,ng) = j + jshift_bc
02994 areafact(i,j,ng) = areafac_c(i,j)
02995
02996
02997
02998 ng = 5
02999 xp (i,j,1,ng) = xcr
03000 yp (i,j,1,ng) = ycr
03001 xp (i,j,2,ng) = xdr
03002 yp (i,j,2,ng) = ydr
03003 xp (i,j,3,ng) = xicr
03004 yp (i,j,3,ng) = yicr
03005 iflux (i,j,ng) = i + ishift_tc
03006 jflux (i,j,ng) = j + jshift_tc
03007 areafact(i,j,ng) = -areafac_r(i,j)
03008
03009
03010
03011 ng = 6
03012 xp (i,j,1,ng) = xicr
03013 yp (i,j,1,ng) = yicr
03014 xp (i,j,2,ng) = xdl
03015 yp (i,j,2,ng) = ydl
03016 xp (i,j,3,ng) = xdm
03017 yp (i,j,3,ng) = ydm
03018 iflux (i,j,ng) = i + ishift_bc
03019 jflux (i,j,ng) = j + jshift_bc
03020 areafact(i,j,ng) = areafac_c(i,j)
03021
03022 elseif (ydl < c0 .and. ydr >= c0 .and. xic >= c0 &
03023 .and. ydm >= c0) then
03024
03025
03026
03027 ng = 4
03028 xp (i,j,1,ng) = xcl
03029 yp (i,j,1,ng) = ycl
03030 xp (i,j,2,ng) = xdl
03031 yp (i,j,2,ng) = ydl
03032 xp (i,j,3,ng) = xicl
03033 yp (i,j,3,ng) = yicl
03034 iflux (i,j,ng) = i + ishift_bc
03035 jflux (i,j,ng) = j + jshift_bc
03036 areafact(i,j,ng) = areafac_c(i,j)
03037
03038
03039
03040 ng = 5
03041 xp (i,j,1,ng) = xcr
03042 yp (i,j,1,ng) = ycr
03043 xp (i,j,2,ng) = xdr
03044 yp (i,j,2,ng) = ydr
03045 xp (i,j,3,ng) = xicr
03046 yp (i,j,3,ng) = yicr
03047 iflux (i,j,ng) = i + ishift_tc
03048 jflux (i,j,ng) = j + jshift_tc
03049 areafact(i,j,ng) = -areafac_r(i,j)
03050
03051
03052
03053 ng = 6
03054 xp (i,j,1,ng) = xicl
03055 yp (i,j,1,ng) = yicl
03056 xp (i,j,2,ng) = xicr
03057 yp (i,j,2,ng) = yicr
03058 xp (i,j,3,ng) = xdm
03059 yp (i,j,3,ng) = ydm
03060 iflux (i,j,ng) = i + ishift_tc
03061 jflux (i,j,ng) = j + jshift_tc
03062 areafact(i,j,ng) = -areafac_c(i,j)
03063
03064 endif
03065
03066 enddo
03067
03068
03069
03070
03071
03072
03073
03074
03075
03076
03077
03078
03079
03080
03081
03082
03083
03084
03085
03086
03087
03088
03089 areasum(:,:) = c0
03090
03091 do ng = 1, ngroups
03092 icells(ng) = 0
03093
03094 do ij = 1, icellsd
03095 i = indxid(ij)
03096 j = indxjd(ij)
03097
03098 triarea(i,j,ng) = p5 * ( (xp(i,j,2,ng)-xp(i,j,1,ng)) * &
03099 (yp(i,j,3,ng)-yp(i,j,1,ng)) &
03100 - (yp(i,j,2,ng)-yp(i,j,1,ng)) * &
03101 (xp(i,j,3,ng)-xp(i,j,1,ng)) ) &
03102 * areafact(i,j,ng)
03103
03104 if (abs(triarea(i,j,ng)) < eps16*areafac_c(i,j)) then
03105 triarea(i,j,ng) = c0
03106 else
03107 icells(ng) = icells(ng) + 1
03108 ic = icells(ng)
03109 indxi(ic,ng) = i
03110 indxj(ic,ng) = j
03111 endif
03112
03113 areasum(i,j) = areasum(i,j) + triarea(i,j,ng)
03114
03115 enddo
03116 enddo
03117
03118 if (l_fixed_area) then
03119 if (bugcheck) then
03120 do ij = 1, icellsd
03121 i = indxid(ij)
03122 j = indxjd(ij)
03123 if (abs(areasum(i,j) - edgearea(i,j)) > eps13*areafac_c(i,j)) then
03124 print*, ''
03125 print*, 'Areas do not add up: m, i, j, edge =', &
03126 my_task, i, j, trim(edge)
03127 print*, 'edgearea =', edgearea(i,j)
03128 print*, 'areasum =', areasum(i,j)
03129 print*, 'areafac_c =', areafac_c(i,j)
03130 print*, ''
03131 print*, 'Triangle areas:'
03132 do ng = 1, ngroups
03133 if (abs(triarea(i,j,ng)) > eps16*abs(areafact(i,j,ng))) then
03134 print*, ng, triarea(i,j,ng)
03135 endif
03136 enddo
03137 endif
03138 enddo
03139 endif
03140
03141 else
03142 do ij = 1, icellsd
03143 i = indxid(ij)
03144 j = indxjd(ij)
03145 edgearea(i,j) = areasum(i,j)
03146 enddo
03147 endif
03148
03149
03150
03151
03152
03153
03154 if (trim(edge) == 'north') then
03155 do ng = 1, ngroups
03156 do nv = 1, nvert
03157 do ij = 1, icells(ng)
03158 i = indxi(ij,ng)
03159 j = indxj(ij,ng)
03160 ishift = iflux(i,j,ng) - i
03161 jshift = jflux(i,j,ng) - j
03162 xp(i,j,nv,ng) = xp(i,j,nv,ng) - c1*ishift
03163 yp(i,j,nv,ng) = yp(i,j,nv,ng) + p5 - c1*jshift
03164 enddo
03165 enddo
03166 enddo
03167 else
03168 do ng = 1, ngroups
03169 do nv = 1, nvert
03170 do ij = 1, icells(ng)
03171 i = indxi(ij,ng)
03172 j = indxj(ij,ng)
03173 ishift = iflux(i,j,ng) - i
03174 jshift = jflux(i,j,ng) - j
03175
03176 w1 = xp(i,j,nv,ng)
03177 xp(i,j,nv,ng) = yp(i,j,nv,ng) + p5 - c1*ishift
03178 yp(i,j,nv,ng) = -w1 - c1*jshift
03179 enddo
03180 enddo
03181 enddo
03182 endif
03183
03184 if (bugcheck) then
03185 do ng = 1, ngroups
03186 do nv = 1, nvert
03187 do j = jb, je
03188 do i = ib, ie
03189 if (abs(triarea(i,j,ng)) > puny) then
03190 if (abs(xp(i,j,nv,ng)) > p5+puny) then
03191 print*, ''
03192 print*, 'WARNING: xp =', xp(i,j,nv,ng)
03193 print*, 'm, i, j, ng, nv =', my_task, i, j, ng, nv
03194
03195
03196
03197
03198 endif
03199 if (abs(yp(i,j,nv,ng)) > p5+puny) then
03200 print*, ''
03201 print*, 'WARNING: yp =', yp(i,j,nv,ng)
03202 print*, 'm, i, j, ng, nv =', my_task, i, j, ng, nv
03203 endif
03204 endif
03205 enddo
03206 enddo
03207 enddo
03208 enddo
03209 endif
03210
03211 end subroutine locate_triangles
03212
03213
03214
03215
03216
03217
03218
03219
03220 subroutine triangle_coordinates (nx_block, ny_block, &
03221 integral_order, icells, &
03222 indxi, indxj, &
03223 xp, yp)
03224
03225
03226
03227
03228
03229
03230
03231
03232
03233
03234
03235
03236
03237
03238
03239
03240
03241
03242
03243
03244
03245
03246
03247
03248
03249
03250
03251
03252
03253
03254
03255
03256
03257
03258
03259
03260 integer (kind=int_kind), intent(in) ::
03261 nx_block, ny_block,
03262 integral_order
03263
03264 integer (kind=int_kind), dimension (ngroups), intent(in) ::
03265 icells
03266
03267 integer (kind=int_kind), dimension (nx_block*ny_block,ngroups),
03268 intent(in) ::
03269 indxi ,
03270 indxj
03271
03272 real (kind=dbl_kind), intent(inout),
03273 dimension (nx_block, ny_block, 0:nvert, ngroups) ::
03274 xp, yp
03275
03276
03277
03278 integer (kind=int_kind) ::
03279 i, j, ij ,
03280 ng
03281
03282
03283 if (integral_order == 1) then
03284
03285 do ng = 1, ngroups
03286 do ij = 1, icells(ng)
03287 i = indxi(ij,ng)
03288 j = indxj(ij,ng)
03289
03290
03291 xp(i,j,0,ng) = p333 &
03292 * (xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng))
03293 yp(i,j,0,ng) = p333 &
03294 * (yp(i,j,1,ng) + yp(i,j,2,ng) + yp(i,j,3,ng))
03295
03296 enddo
03297 enddo
03298
03299 elseif (integral_order == 2) then
03300
03301
03302
03303
03304
03305 do ng = 1, ngroups
03306 do ij = 1, icells(ng)
03307 i = indxi(ij,ng)
03308 j = indxj(ij,ng)
03309
03310
03311 xp(i,j,0,ng) = p333 &
03312 * (xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng))
03313 yp(i,j,0,ng) = p333 &
03314 * (yp(i,j,1,ng) + yp(i,j,2,ng) + yp(i,j,3,ng))
03315
03316
03317
03318 xp(i,j,1,ng) = p5*xp(i,j,1,ng) + p5*xp(i,j,0,ng)
03319 yp(i,j,1,ng) = p5*yp(i,j,1,ng) + p5*yp(i,j,0,ng)
03320
03321 xp(i,j,2,ng) = p5*xp(i,j,2,ng) + p5*xp(i,j,0,ng)
03322 yp(i,j,2,ng) = p5*yp(i,j,2,ng) + p5*yp(i,j,0,ng)
03323
03324 xp(i,j,3,ng) = p5*xp(i,j,3,ng) + p5*xp(i,j,0,ng)
03325 yp(i,j,3,ng) = p5*yp(i,j,3,ng) + p5*yp(i,j,0,ng)
03326
03327 enddo
03328 enddo
03329
03330 else
03331
03332
03333
03334
03335 do ng = 1, ngroups
03336 do ij = 1, icells(ng)
03337 i = indxi(ij,ng)
03338 j = indxj(ij,ng)
03339
03340
03341 xp(i,j,0,ng) = p333 &
03342 * (xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng))
03343 yp(i,j,0,ng) = p333 &
03344 * (yp(i,j,1,ng) + yp(i,j,2,ng) + yp(i,j,3,ng))
03345
03346
03347
03348 xp(i,j,1,ng) = p4*xp(i,j,1,ng) + p6*xp(i,j,0,ng)
03349 yp(i,j,1,ng) = p4*yp(i,j,1,ng) + p6*yp(i,j,0,ng)
03350
03351 xp(i,j,2,ng) = p4*xp(i,j,2,ng) + p6*xp(i,j,0,ng)
03352 yp(i,j,2,ng) = p4*yp(i,j,2,ng) + p6*yp(i,j,0,ng)
03353
03354 xp(i,j,3,ng) = p4*xp(i,j,3,ng) + p6*xp(i,j,0,ng)
03355 yp(i,j,3,ng) = p4*yp(i,j,3,ng) + p6*yp(i,j,0,ng)
03356
03357 enddo
03358 enddo
03359
03360 endif
03361
03362 end subroutine triangle_coordinates
03363
03364
03365
03366
03367
03368
03369
03370
03371
03372 subroutine transport_integrals (nx_block, ny_block, &
03373 icells, &
03374 indxi, indxj, &
03375 tracer_type, depend, &
03376 integral_order, triarea, &
03377 iflux, jflux, &
03378 xp, yp, &
03379 mc, mx, &
03380 my, mflx, &
03381 tc, tx, &
03382 ty, mtflx)
03383
03384
03385
03386
03387
03388
03389
03390
03391
03392
03393
03394
03395
03396
03397
03398
03399
03400 integer (kind=int_kind), intent(in) ::
03401 nx_block, ny_block ,
03402 integral_order
03403
03404 integer (kind=int_kind), dimension (ntrace), intent(in) ::
03405 tracer_type ,
03406 depend
03407
03408 integer (kind=int_kind), dimension (ngroups), intent(in) ::
03409 icells
03410
03411 integer (kind=int_kind), dimension (nx_block*ny_block,ngroups),
03412 intent(in) ::
03413 indxi ,
03414 indxj
03415
03416 real (kind=dbl_kind), intent(in),
03417 dimension (nx_block, ny_block, 0:nvert, ngroups) ::
03418 xp, yp
03419
03420 real (kind=dbl_kind), intent(in),
03421 dimension (nx_block, ny_block, ngroups) ::
03422 triarea
03423
03424 integer (kind=int_kind), intent(in),
03425 dimension (nx_block, ny_block, ngroups) ::
03426 iflux ,
03427 jflux
03428
03429 real (kind=dbl_kind), intent(in),
03430 dimension (nx_block, ny_block) ::
03431 mc, mx, my
03432
03433 real (kind=dbl_kind), intent(out),
03434 dimension (nx_block, ny_block) ::
03435 mflx
03436
03437 real (kind=dbl_kind), intent(in),
03438 dimension (nx_block, ny_block, ntrace), optional ::
03439 tc, tx, ty
03440
03441 real (kind=dbl_kind), intent(out),
03442 dimension (nx_block, ny_block, ntrace), optional ::
03443 mtflx
03444
03445
03446
03447 integer (kind=int_kind) ::
03448 i, j, ij ,
03449 i2, j2 ,
03450 ng ,
03451 nt, nt1 ,
03452 ilo,ihi,jlo,jhi
03453
03454 real (kind=dbl_kind) ::
03455 m0, m1, m2, m3 ,
03456 w0, w1, w2, w3
03457
03458 real (kind=dbl_kind), dimension (nx_block, ny_block) ::
03459 msum, mxsum, mysum ,
03460 mxxsum, mxysum, myysum
03461
03462 real (kind=dbl_kind), dimension (nx_block, ny_block, ntrace) ::
03463 mtsum ,
03464 mtxsum ,
03465 mtysum
03466
03467
03468
03469
03470
03471 mflx(:,:) = c0
03472 if (present(mtflx)) then
03473 do nt = 1, ntrace
03474 mtflx(:,:,nt) = c0
03475 enddo
03476 endif
03477
03478
03479
03480
03481
03482 do ng = 1, ngroups
03483
03484 if (integral_order == 1) then
03485
03486
03487
03488
03489 do ij = 1, icells(ng)
03490 i = indxi(ij,ng)
03491 j = indxj(ij,ng)
03492
03493 i2 = iflux(i,j,ng)
03494 j2 = jflux(i,j,ng)
03495
03496
03497
03498 m0 = mc(i2,j2) + xp(i,j,0,ng)*mx(i2,j2) &
03499 + yp(i,j,0,ng)*my(i2,j2)
03500 msum(i,j) = m0
03501
03502 mflx(i,j) = mflx(i,j) + triarea(i,j,ng)*msum(i,j)
03503
03504
03505 mxsum(i,j) = m0*xp(i,j,0,ng)
03506 mxxsum(i,j) = mxsum(i,j)*xp(i,j,0,ng)
03507 mxysum(i,j) = mxsum(i,j)*yp(i,j,0,ng)
03508 mysum(i,j) = m0*yp(i,j,0,ng)
03509 myysum(i,j) = mysum(i,j)*yp(i,j,0,ng)
03510 enddo
03511
03512 elseif (integral_order == 2) then
03513
03514
03515
03516
03517 do ij = 1, icells(ng)
03518 i = indxi(ij,ng)
03519 j = indxj(ij,ng)
03520
03521 i2 = iflux(i,j,ng)
03522 j2 = jflux(i,j,ng)
03523
03524
03525
03526
03527 m1 = p333 * (mc(i2,j2) + xp(i,j,1,ng)*mx(i2,j2) &
03528 + yp(i,j,1,ng)*my(i2,j2))
03529 m2 = p333 * (mc(i2,j2) + xp(i,j,2,ng)*mx(i2,j2) &
03530 + yp(i,j,2,ng)*my(i2,j2))
03531 m3 = p333 * (mc(i2,j2) + xp(i,j,3,ng)*mx(i2,j2) &
03532 + yp(i,j,3,ng)*my(i2,j2))
03533 msum(i,j) = m1 + m2 + m3
03534 mflx(i,j) = mflx(i,j) + triarea(i,j,ng)*msum(i,j)
03535
03536
03537 w1 = m1 * xp(i,j,1,ng)
03538 w2 = m2 * xp(i,j,2,ng)
03539 w3 = m3 * xp(i,j,3,ng)
03540
03541 mxsum(i,j) = w1 + w2 + w3
03542
03543 mxxsum(i,j) = w1*xp(i,j,1,ng) + w2*xp(i,j,2,ng) &
03544 + w3*xp(i,j,3,ng)
03545
03546 mxysum(i,j) = w1*yp(i,j,1,ng) + w2*yp(i,j,2,ng) &
03547 + w3*yp(i,j,3,ng)
03548
03549 w1 = m1 * yp(i,j,1,ng)
03550 w2 = m2 * yp(i,j,2,ng)
03551 w3 = m3 * yp(i,j,3,ng)
03552
03553 mysum(i,j) = w1 + w2 + w3
03554
03555 myysum(i,j) = w1*yp(i,j,1,ng) + w2*yp(i,j,2,ng) &
03556 + w3*yp(i,j,3,ng)
03557 enddo
03558
03559 else
03560
03561
03562
03563
03564 do ij = 1, icells(ng)
03565 i = indxi(ij,ng)
03566 j = indxj(ij,ng)
03567
03568 i2 = iflux(i,j,ng)
03569 j2 = jflux(i,j,ng)
03570
03571
03572
03573
03574
03575 m0 = p5625m * (mc(i2,j2) + xp(i,j,0,ng)*mx(i2,j2) &
03576 + yp(i,j,0,ng)*my(i2,j2))
03577 m1 = p52083 * (mc(i2,j2) + xp(i,j,1,ng)*mx(i2,j2) &
03578 + yp(i,j,1,ng)*my(i2,j2))
03579 m2 = p52083 * (mc(i2,j2) + xp(i,j,2,ng)*mx(i2,j2) &
03580 + yp(i,j,2,ng)*my(i2,j2))
03581 m3 = p52083 * (mc(i2,j2) + xp(i,j,3,ng)*mx(i2,j2) &
03582 + yp(i,j,3,ng)*my(i2,j2))
03583 msum(i,j) = m0 + m1 + m2 + m3
03584 mflx(i,j) = mflx(i,j) + triarea(i,j,ng)*msum(i,j)
03585
03586
03587 w0 = m0 * xp(i,j,0,ng)
03588 w1 = m1 * xp(i,j,1,ng)
03589 w2 = m2 * xp(i,j,2,ng)
03590 w3 = m3 * xp(i,j,3,ng)
03591
03592 mxsum(i,j) = w0 + w1 + w2 + w3
03593
03594 mxxsum(i,j) = w0*xp(i,j,0,ng) + w1*xp(i,j,1,ng) &
03595 + w2*xp(i,j,2,ng) + w3*xp(i,j,3,ng)
03596
03597 mxysum(i,j) = w0*yp(i,j,0,ng) + w1*yp(i,j,1,ng) &
03598 + w2*yp(i,j,2,ng) + w3*yp(i,j,3,ng)
03599
03600 w0 = m0 * yp(i,j,0,ng)
03601 w1 = m1 * yp(i,j,1,ng)
03602 w2 = m2 * yp(i,j,2,ng)
03603 w3 = m3 * yp(i,j,3,ng)
03604
03605 mysum(i,j) = w0 + w1 + w2 + w3
03606
03607 myysum(i,j) = w0*yp(i,j,0,ng) + w1*yp(i,j,1,ng) &
03608 + w2*yp(i,j,2,ng) + w3*yp(i,j,3,ng)
03609
03610 enddo
03611
03612 endif
03613
03614
03615
03616 if (present(mtflx)) then
03617
03618 do nt = 1, ntrace
03619 if (tracer_type(nt)==1) then
03620
03621
03622
03623
03624 do ij = 1, icells(ng)
03625 i = indxi(ij,ng)
03626 j = indxj(ij,ng)
03627
03628 i2 = iflux(i,j,ng)
03629 j2 = jflux(i,j,ng)
03630
03631 mtsum(i,j,nt) = msum(i,j) * tc(i2,j2,nt) &
03632 + mxsum(i,j) * tx(i2,j2,nt) &
03633 + mysum(i,j) * ty(i2,j2,nt)
03634
03635 mtflx(i,j,nt) = mtflx(i,j,nt) &
03636 + triarea(i,j,ng) * mtsum(i,j,nt)
03637
03638
03639
03640 mtxsum(i,j,nt) = mxsum(i,j) * tc(i2,j2,nt) &
03641 + mxxsum(i,j) * tx(i2,j2,nt) &
03642 + mxysum(i,j) * ty(i2,j2,nt)
03643
03644 mtysum(i,j,nt) = mysum(i,j) * tc(i2,j2,nt) &
03645 + mxysum(i,j) * tx(i2,j2,nt) &
03646 + myysum(i,j) * ty(i2,j2,nt)
03647 enddo
03648
03649 elseif (tracer_type(nt)==2) then
03650 nt1 = depend(nt)
03651
03652
03653
03654
03655 do ij = 1, icells(ng)
03656 i = indxi(ij,ng)
03657 j = indxj(ij,ng)
03658
03659 i2 = iflux(i,j,ng)
03660 j2 = jflux(i,j,ng)
03661
03662 mtsum(i,j,nt) = mtsum(i,j,nt1) * tc(i2,j2,nt) &
03663 + mtxsum(i,j,nt1) * tx(i2,j2,nt) &
03664 + mtysum(i,j,nt1) * ty(i2,j2,nt)
03665
03666 mtflx(i,j,nt) = mtflx(i,j,nt) &
03667 + triarea(i,j,ng) * mtsum(i,j,nt)
03668 enddo
03669
03670
03671 elseif (tracer_type(nt)==3) then
03672 nt1 = depend(nt)
03673
03674
03675
03676
03677 do ij = 1, icells(ng)
03678 i = indxi(ij,ng)
03679 j = indxj(ij,ng)
03680
03681 i2 = iflux(i,j,ng)
03682 j2 = jflux(i,j,ng)
03683
03684
03685 mtsum(i,j,nt) = mtsum(i,j,nt1) * tc(i2,j2,nt)
03686
03687 mtflx(i,j,nt) = mtflx(i,j,nt) &
03688 + triarea(i,j,ng) * mtsum(i,j,nt)
03689 enddo
03690
03691 endif
03692 enddo
03693 endif
03694 enddo
03695
03696 end subroutine transport_integrals
03697
03698
03699
03700
03701
03702
03703
03704
03705
03706 subroutine update_fields (nx_block, ny_block, &
03707 ilo, ihi, jlo, jhi, &
03708 tracer_type, depend, &
03709 tarear, l_stop, &
03710 istop, jstop, &
03711 mflxe, mflxn, &
03712 mm, &
03713 mtflxe, mtflxn, &
03714 tm)
03715
03716
03717
03718
03719
03720
03721
03722
03723
03724
03725
03726
03727
03728 integer (kind=int_kind), intent(in) ::
03729 nx_block, ny_block,
03730 ilo,ihi,jlo,jhi
03731
03732 integer (kind=int_kind), dimension (ntrace), intent(in) ::
03733 tracer_type ,
03734 depend
03735
03736 real (kind=dbl_kind), dimension (nx_block, ny_block),
03737 intent(in) ::
03738 mflxe, mflxn ,
03739 tarear
03740
03741 real (kind=dbl_kind), dimension (nx_block, ny_block),
03742 intent(inout) ::
03743 mm
03744
03745 real (kind=dbl_kind), dimension (nx_block, ny_block, ntrace),
03746 intent(in), optional ::
03747 mtflxe, mtflxn
03748
03749 real (kind=dbl_kind), dimension (nx_block, ny_block, ntrace),
03750 intent(inout), optional ::
03751 tm
03752
03753 logical (kind=log_kind), intent(inout) ::
03754 l_stop
03755
03756 integer (kind=int_kind), intent(inout) ::
03757 istop, jstop
03758
03759
03760
03761
03762 integer (kind=int_kind) ::
03763 i, j ,
03764 nt, nt1, nt2
03765
03766 real (kind=dbl_kind), dimension(nx_block,ny_block,ntrace) ::
03767 mtold
03768
03769 real (kind=dbl_kind) ::
03770 w1, w2
03771
03772 integer (kind=int_kind), dimension(nx_block*ny_block) ::
03773 indxi ,
03774 indxj
03775
03776 integer (kind=int_kind) ::
03777 icells ,
03778 ij
03779
03780
03781
03782
03783
03784 if (present(tm)) then
03785 do nt = 1, ntrace
03786 if (tracer_type(nt)==1) then
03787 do j = jlo, jhi
03788 do i = ilo, ihi
03789 mtold(i,j,nt) = mm(i,j) * tm(i,j,nt)
03790 enddo
03791 enddo
03792 elseif (tracer_type(nt)==2) then
03793 nt1 = depend(nt)
03794 do j = jlo, jhi
03795 do i = ilo, ihi
03796 mtold(i,j,nt) = mm(i,j) * tm(i,j,nt1) * tm(i,j,nt)
03797 enddo
03798 enddo
03799 elseif (tracer_type(nt)==3) then
03800 nt1 = depend(nt)
03801 nt2 = depend(nt1)
03802 do j = jlo, jhi
03803 do i = ilo, ihi
03804 mtold(i,j,nt) = mm(i,j) &
03805 * tm(i,j,nt2) * tm(i,j,nt1) * tm(i,j,nt)
03806 enddo
03807 enddo
03808
03809
03810 endif
03811 enddo
03812 endif
03813
03814
03815
03816
03817
03818 do j = jlo, jhi
03819 do i = ilo, ihi
03820
03821 w1 = mflxe(i,j) - mflxe(i-1,j) &
03822 + mflxn(i,j) - mflxn(i,j-1)
03823 mm(i,j) = mm(i,j) - w1*tarear(i,j)
03824
03825 if (mm(i,j) < -puny) then
03826 l_stop = .true.
03827 istop = i
03828 jstop = j
03829 elseif (mm(i,j) < c0) then
03830 mm(i,j) = c0
03831 endif
03832
03833 enddo
03834 enddo
03835
03836 if (l_stop) then
03837 i = istop
03838 j = jstop
03839 w1 = mflxe(i,j) - mflxe(i-1,j) &
03840 + mflxn(i,j) - mflxn(i,j-1)
03841 write (nu_diag,*) ' '
03842 write (nu_diag,*) 'New mass < 0, i, j =', i, j
03843 write (nu_diag,*) 'Old mass =', mm(i,j) + w1*tarear(i,j)
03844 write (nu_diag,*) 'New mass =', mm(i,j)
03845 write (nu_diag,*) 'Net transport =', -w1*tarear(i,j)
03846 return
03847 endif
03848
03849
03850
03851
03852
03853 if (present(tm)) then
03854
03855 icells = 0
03856 do j = jlo, jhi
03857 do i = ilo, ihi
03858 if (mm(i,j) > c0) then
03859 icells = icells + 1
03860 indxi(icells) = i
03861 indxj(icells) = j
03862 endif
03863 enddo
03864 enddo
03865
03866 do nt = 1, ntrace
03867
03868 do j = jlo, jhi
03869 do i = ilo, ihi
03870 tm(i,j,nt) = c0
03871 enddo
03872 enddo
03873
03874 if (tracer_type(nt)==1) then
03875
03876 do ij = 1, icells
03877 i = indxi(ij)
03878 j = indxj(ij)
03879
03880 w1 = mtflxe(i,j,nt) - mtflxe(i-1,j,nt) &
03881 + mtflxn(i,j,nt) - mtflxn(i,j-1,nt)
03882 tm(i,j,nt) = (mtold(i,j,nt) - w1*tarear(i,j)) &
03883 / mm(i,j)
03884 enddo
03885
03886
03887 elseif (tracer_type(nt)==2) then
03888 nt1 = depend(nt)
03889
03890
03891
03892
03893 do ij = 1, icells
03894 i = indxi(ij)
03895 j = indxj(ij)
03896
03897 if (abs(tm(i,j,nt1)) > c0) then
03898 w1 = mtflxe(i,j,nt) - mtflxe(i-1,j,nt) &
03899 + mtflxn(i,j,nt) - mtflxn(i,j-1,nt)
03900 tm(i,j,nt) = (mtold(i,j,nt) - w1*tarear(i,j)) &
03901 / (mm(i,j) * tm(i,j,nt1))
03902 endif
03903 enddo
03904
03905 elseif (tracer_type(nt)==3) then
03906 nt1 = depend(nt)
03907 nt2 = depend(nt1)
03908
03909
03910
03911
03912 do ij = 1, icells
03913 i = indxi(ij)
03914 j = indxj(ij)
03915
03916 if (abs(tm(i,j,nt1)) > c0 .and. &
03917 abs(tm(i,j,nt2)) > c0) then
03918 w1 = mtflxe(i,j,nt) - mtflxe(i-1,j,nt) &
03919 + mtflxn(i,j,nt) - mtflxn(i,j-1,nt)
03920 tm(i,j,nt) = (mtold(i,j,nt) - w1*tarear(i,j)) &
03921 / (mm(i,j) * tm(i,j,nt2) * tm(i,j,nt1))
03922 endif
03923 enddo
03924 endif
03925 enddo
03926 endif
03927
03928 end subroutine update_fields
03929
03930
03931
03932 end module ice_transport_remap
03933
03934