00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011 subroutine psmile_info_trs_loc_gauss2_real (comp_info, &
00012 coords, shape, control, len_cpl, &
00013 var_id, grid_valid_shape, &
00014 search, method_id, &
00015 send_index, ierror)
00016
00017
00018
00019 use PRISM_constants
00020
00021 use PSMILe, dummy_interface => PSMILe_info_trs_loc_gauss2_real
00022 #ifdef DEBUG_TRACE
00023 use psmile_debug_trace
00024 #endif
00025
00026 implicit none
00027
00028
00029
00030 Type (Enddef_comp), Intent (In) :: comp_info
00031
00032
00033
00034
00035 Type (Enddef_search), Intent (InOut) :: search
00036
00037
00038
00039 Type (real_vector), Intent (In) :: coords (ndim_3d, search%npart)
00040
00041
00042
00043 Integer, Intent (In) :: shape (2, ndim_3d, search%npart)
00044
00045
00046
00047
00048
00049 Integer, Intent (In) :: control (2, ndim_3d, search%npart)
00050
00051
00052
00053 Integer, Intent (InOut) :: len_cpl (search%npart)
00054
00055
00056
00057 Integer, Intent (In) :: var_id
00058
00059
00060
00061 Integer, Intent (In) :: grid_valid_shape (2, ndim_3d)
00062
00063
00064
00065 Integer, Intent (In) :: method_id
00066
00067
00068
00069 Integer, Intent (InOut) :: send_index
00070
00071
00072
00073
00074
00075 integer, Intent (Out) :: ierror
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086 Integer, parameter :: lon = 1
00087 Integer, parameter :: lat = 2
00088
00089
00090
00091
00092 Integer, Parameter :: indl = 1
00093 Integer, Parameter :: indz = 2
00094
00095
00096
00097
00098
00099 Type (Gridfunction), Pointer :: field
00100
00101
00102
00103 Type (Corner_Block), Pointer :: cp
00104
00105
00106
00107 Type (Method), Pointer :: mp
00108
00109
00110
00111 Type (Grid), Pointer :: grid_info
00112
00113 Integer :: grid_id
00114 Integer :: ncpl, new_ncpl
00115 Integer :: jj, totsiz, offset
00116
00117
00118
00119 Integer :: i, j, k, n, nn
00120 Integer :: jpart, n1, n2, n3
00121
00122
00123
00124 Type (Send_information), Pointer :: send_info
00125 Type (Extra_search_info) :: extra_search
00126
00127 Integer :: nb_neighbors, nb_extra
00128 Logical :: virtual_cell_available
00129
00130 Integer, Pointer :: virtual_cell (:)
00131 Integer, Pointer :: dstijk (:, :)
00132 Integer, Pointer :: new_dstijk (:, :)
00133
00134 Integer, Target :: dummy_virtual_cell (1)
00135
00136 Integer, Allocatable :: neighbors_3d (:, :, :)
00137 Integer :: interpolation_mode, search_mode
00138 Integer :: interpolation_type
00139 Integer :: interpolation_methods (ndim_3d)
00140 Integer :: n_intmethods
00141 Logical :: interpolation_search (ndim_3d)
00142 Logical :: global_search
00143 Logical :: global_cell_search
00144
00145
00146
00147
00148
00149
00150
00151 Integer :: id_src_mask, id_tgt_mask
00152 Integer :: mask_id
00153 Integer :: use_how (3)
00154 Logical :: src_mask_available, use_mask
00155
00156 Logical, Target :: dummy_mask_array (1)
00157
00158 Integer :: dummy_mask_shape (2, ndim_3d)
00159 Integer, Target :: dummy_mask (1)
00160
00161
00162 Integer :: mask_shape (2, ndim_3d)
00163 Logical, Pointer :: mask_array (:)
00164 Integer, Pointer :: src_mask (:)
00165
00166
00167
00168 Integer :: ibeg, iend, ipart, nprev
00169 Integer :: new_len_cpl(search%npart)
00170
00171
00172
00173 Integer :: coords_shape (2, ndim_3d)
00174 Integer :: grid_shape (2, ndim_3d)
00175
00176 Integer :: cpl_id
00177 Integer :: epio_id, trs_rank
00178 Integer :: local_src_size, shift, src_size
00179 Integer :: src_cell_size
00180 Integer :: id_trans_out
00181 Integer :: ip, len, leni, lenij
00182
00183 Integer, Pointer :: list_entries (:, :)
00184 Type (real_vector) :: src_grid (ndim_3d)
00185 Type (real_vector) :: src_cell (ndim_3d)
00186 Integer, Allocatable :: neigh_bascule (:)
00187 Integer, Allocatable :: neighbors (:, :)
00188 Integer, Allocatable :: neighcells (:, :)
00189 Integer, Allocatable :: source_cell_index (:)
00190 Integer, Pointer :: nbr_cells (:)
00191 Integer, Pointer :: new_nbr_cells (:)
00192 Integer :: nbr_cells_tot, nb_corners
00193 Integer :: nbr_corners (ndim_3d)
00194 Integer :: corner_shape (2, ndim_3d)
00195 Integer :: corner_shape_3d (2, ndim_3d, ndim_3d)
00196
00197
00198
00199 Integer :: id_trans_in, dest_comp_id
00200 Logical :: use_target_grid
00201 Logical :: use_target_cell
00202 Logical :: use_source_cell
00203 Integer :: tgt_corners(ndim_3d)
00204
00205 Type (real_vector) :: tgt_grid (ndim_3d)
00206 Type (real_vector) :: tgt_cell (ndim_3d)
00207 Type (real_vector) :: new_tgt_cell (ndim_3d)
00208 Integer, Pointer :: tgt_mask(:)
00209
00210
00211
00212 logical :: transformed, cubic
00213 Type (real_vector) :: sinvec (ndim_2d), cosvec (ndim_2d)
00214
00215
00216
00217 Integer, parameter :: nerrp = 3
00218 Integer :: ierrp (nerrp)
00219
00220
00221
00222
00223
00224
00225
00226
00227
00228
00229
00230
00231
00232
00233
00234
00235
00236
00237
00238
00239
00240
00241
00242 Character(len=len_cvs_string), save :: mycvs =
00243 '$Id: psmile_info_trs_loc_gauss2_real.F90 2994 2011-02-25 11:26:52Z hanke $'
00244
00245
00246
00247
00248
00249 data dummy_mask_shape /0, 1, 0, 1, 0, 1/
00250
00251
00252
00253
00254
00255 #ifdef VERBOSE
00256 print 9990, trim(ch_id), method_id, search%msg_intersections%first_tgt_method_id, &
00257 search%sender
00258
00259 call psmile_flushstd
00260 #endif /* VERBOSE */
00261
00262 ierror = 0
00263
00264 field => Fields(var_id)
00265 mp => Methods(method_id)
00266 grid_id = mp%grid_id
00267 send_info => mp%send_infos_coupler(send_index)
00268 grid_info => Grids(grid_id)
00269 cp => grid_info%corner_pointer
00270
00271 transformed = .false.
00272 cubic = .false.
00273
00274 global_search = .false.
00275 global_cell_search = .false.
00276 interpolation_search (:) = .false.
00277
00278 use_how(1) = PSMILe_undef
00279 use_how(2) = PRISM_Gaussreduced_Regvrt
00280 use_how(3) = grid_id
00281
00282 #ifdef PRISM_ASSERTION
00283
00284
00285
00286 if ( .not. associated(field%Taskout) ) then
00287 print *, 'var_id', var_id
00288 call psmile_assert (__FILE__, &
00289 __LINE__, "Taskout is not available for Field")
00290 end if
00291
00292 if (var_id < 1 .or. var_id > Number_of_Fields_allocated) then
00293 print *, 'var_id', var_id
00294 call psmile_assert (__FILE__, &
00295 __LINE__, "Incorrect var_id")
00296 endif
00297
00298 if ( mp%status == PSMILe_Status_free .or. &
00299 mp%status == PSMILe_Status_undefined ) then
00300 call psmile_assert (__FILE__, &
00301 __LINE__, "Incorrect method")
00302 endif
00303
00304 do ipart = 1, search%npart
00305 if (control (1,1,ipart) < search%range (1,1,ipart) .or. &
00306 search%range (2,1,ipart) < control (2,1,ipart) .or. &
00307 control (1,2,ipart) < search%range (1,2,ipart) .or. &
00308 search%range (2,2,ipart) < control (2,2,ipart) .or. &
00309 control (1,3,ipart) < search%range (1,3,ipart) .or. &
00310 search%range (2,3,ipart) < control (2,3,ipart)) then
00311 print *, ipart, control (:, :, ipart), search%range (:, :, ipart)
00312 call psmile_assert (__FILE__, &
00313 __LINE__, "control is out of range")
00314 endif
00315 end do
00316
00317 if (send_index < 1) then
00318 print *, send_index
00319 call psmile_assert (__FILE__, &
00320 __LINE__, "send_index is out of range")
00321 endif
00322
00323 if (send_info%nloc < 1) then
00324 print *, "ncpl", send_info%nloc
00325 call psmile_assert (__FILE__, &
00326 __LINE__, "ncpl > 0 expected")
00327 endif
00328
00329 if (send_info%nvec == 2) then
00330 if (search%grid_type /= PRISM_Irrlonlat_Regvrt .and. &
00331 search%grid_type /= PRISM_Gaussreduced_Regvrt) then
00332 print *, "nvec", send_info%nvec, search%grid_type
00333 call psmile_assert (__FILE__, __LINE__, &
00334 "nvec == 2 expected only for PRISM_Irrlonlat_Regvrt")
00335 endif
00336 else if (send_info%nvec == ndim_3d) then
00337 if (search%grid_type /= PRISM_Reglonlatvrt) then
00338 print *, "nvec", send_info%nvec, search%grid_type
00339 call psmile_assert (__FILE__, __LINE__, &
00340 "nvec == ndim_3d expected ony for PRISM_Reglonlatvrt")
00341 endif
00342 else if (send_info%nvec /= 1) then
00343 print *, "nvec", send_info%nvec
00344 call psmile_assert (__FILE__, __LINE__, &
00345 "nvec must be in range [1:3]")
00346 endif
00347
00348 if (send_info%nloc /= Sum(len_cpl)) then
00349 print *, "ncpl", send_info%nloc, Sum (len_cpl)
00350 call psmile_assert (__FILE__, &
00351 __LINE__, "ncpl == SUM(len_cpl) expected")
00352 endif
00353 #endif
00354
00355
00356
00357 ncpl = send_info%nloc
00358 dstijk => send_info%dstijk
00359
00360
00361
00362
00363
00364
00365
00366 mask_id = field%mask_id
00367 src_mask_available = mask_id /= PRISM_UNDEFINED
00368 use_mask = src_mask_available
00369
00370 if (src_mask_available) then
00371 mask_array => Masks(mask_id)%mask_array
00372
00373 mask_shape = Masks(mask_id)%mask_shape
00374 else
00375 mask_array => dummy_mask_array
00376
00377 mask_shape = dummy_mask_shape
00378 endif
00379
00380
00381
00382
00383
00384
00385 id_trans_in = search%msg_intersections%transient_in_id
00386 id_trans_out = search%msg_intersections%transient_out_id
00387 dest_comp_id = all_comp_infos (search%msg_intersections%all_comp_infos_comp_idx)%global_comp_id
00388
00389 do i = 1, size (field%Taskout)
00390 if (field%Taskout(i)%origin_type == PSMILe_comp .and. &
00391 field%Taskout(i)%remote_comp_id == dest_comp_id .and. &
00392 field%Taskout(i)%remote_transi_id == id_trans_in .and. &
00393 field%Taskout(i)%global_transi_id == id_trans_out) exit
00394 end do
00395
00396 if (i > size (field%Taskout)) then
00397 #ifdef DEBUG
00398 print *, 'psmile_info_trs_loc_gauss2_real: method_id ', method_id
00399 print *, 'search for Taskout id ', id_trans_out, ' and Taskin Id ', id_trans_in, &
00400 ' for dest component', dest_comp_id
00401
00402 do i = 1, size (field%Taskout)
00403 print 9960, i, &
00404 field%Taskout(i)%origin_type, &
00405 field%Taskout(i)%remote_comp_id, &
00406 field%Taskout(i)%remote_transi_id, &
00407 field%Taskout(i)%global_transi_id
00408 end do
00409 #endif /* DEBUG */
00410
00411 ierrp (1) = id_trans_in
00412 ierrp (2) = id_trans_out
00413 ierror = PRISM_Error_Internal
00414
00415 call psmile_error ( ierror, 'cannot find Taskout index for transient ids', &
00416 ierrp, 2, __FILE__, __LINE__ )
00417 endif
00418
00419 interpolation_type = field%Taskout(i)%interp%interp_type
00420 interpolation_methods = field%Taskout(i)%interp%interp_meth
00421
00422
00423
00424
00425
00426
00427 #ifdef DEBUG
00428 print *, "interpolation type ", interpolation_type
00429 print *, "interpolation method", interpolation_methods
00430 print *, "number of targets (ncpl) is now ", ncpl
00431 #endif
00432
00433 select case (interpolation_type)
00434 case (PSMILe_3D)
00435 n_intmethods = 1
00436
00437 case (PSMILe_2D1D)
00438 n_intmethods = 2
00439
00440 case (PSMILe_1D1D1d)
00441 n_intmethods = 3
00442
00443
00444
00445 case default
00446 ierrp (1) = interpolation_type
00447 ierror = PRISM_Error_Interp_type
00448
00449 call psmile_error ( ierror, 'unsupported 3d interpolation type', &
00450 ierrp, 1, __FILE__, __LINE__ )
00451 return
00452
00453 end select
00454
00455 do j = n_intmethods, 1, -1
00456 if (interpolation_methods(j) /= PSMILe_none) exit
00457 end do
00458
00459 n_intmethods = j
00460
00461 if (n_intmethods == 0) then
00462 ierror = PRISM_Error_Interp_type
00463 call psmile_error ( ierror, 'All interpolations are PSMILE_none', &
00464 ierrp, 0, __FILE__, __LINE__ )
00465 return
00466 endif
00467
00468
00469
00470 nb_neighbors = 1
00471 nb_extra = 1
00472
00473 do j = 1, n_intmethods
00474 select case (interpolation_methods(j))
00475
00476
00477
00478
00479
00480 case (PSMILe_conserv2D)
00481 nb_neighbors = nb_neighbors * 1
00482 use_how(1) = field%Taskout(i)%interp%arg4(j)
00483 use_how(2) = PSMILe_conserv2D
00484 interpolation_search (j) = &
00485 field%Taskout(i)%interp%arg3 (j) == PSMILE_global
00486 case (PSMILe_conserv3D)
00487 print *, "WARNING: conservative 3D remapping not yet implemented", nb_neighbors
00488 nb_neighbors = nb_neighbors * 2
00489 use_how(1) = field%Taskout(i)%interp%arg4(j)
00490 use_how(2) = PSMILe_conserv3D
00491 interpolation_search (j) = &
00492 field%Taskout(i)%interp%arg3 (j) == PSMILE_global
00493
00494
00495
00496
00497
00498 case (PSMILe_nnghbr3D,PSMILe_nnghbr2D)
00499 #ifdef DEBUG
00500 print *, "Number of nneigbours", field%Taskout(i)%interp%arg2(j)
00501 #endif
00502 use_how(1) = field%Taskout(i)%interp%arg4(j)
00503 nb_neighbors = nb_neighbors * field%Taskout(i)%interp%arg2(j)
00504 nb_extra = nb_extra * field%Taskout(i)%interp%arg2(j)
00505 use_mask = use_mask .and. &
00506 (field%Taskout(i)%interp%arg5(j) /= 1)
00507
00508
00509
00510
00511
00512 case (PSMILe_trilinear)
00513 nb_neighbors = nb_neighbors * 8
00514 use_how(1) = field%Taskout(i)%interp%arg4(j)
00515 interpolation_search (j) = &
00516 field%Taskout(i)%interp%arg3 (j) == PSMILE_global
00517
00518 case (PSMILe_bilinear)
00519 nb_neighbors = nb_neighbors * 4
00520 use_how(1) = field%Taskout(i)%interp%arg4(j)
00521 interpolation_search (j) = &
00522 field%Taskout(i)%interp%arg3 (j) == PSMILE_global
00523
00524 case (PSMILe_linear)
00525
00526
00527 nb_neighbors = nb_neighbors * 2
00528 use_how(1) = field%Taskout(i)%interp%arg4(j)
00529 interpolation_search (j) = &
00530 field%Taskout(i)%interp%arg3 (j) == PSMILE_global
00531
00532
00533
00534
00535 case (PSMILe_bicubic)
00536 nb_neighbors = nb_neighbors * 4 * 4
00537 use_how(1) = field%Taskout(i)%interp%arg4(j)
00538 cubic = .true.
00539 interpolation_search (j) = &
00540 field%Taskout(i)%interp%arg3 (j) == PSMILE_global
00541
00542 case (PSMILe_none)
00543
00544 case (PSMILe_undef)
00545 ierrp (1) = j
00546 ierror = PRISM_Error_Interp_type
00547
00548 call psmile_error ( ierror, 'undefined interpolation method (PSMILE_UNDEF) found', &
00549 ierrp, 1, __FILE__, __LINE__ )
00550 return
00551
00552 case default
00553 ierrp (1) = interpolation_methods(j)
00554 ierror = PRISM_Error_Internal
00555
00556 call psmile_error ( ierror, 'unsupported interpolation method', &
00557 ierrp, 1, __FILE__, __LINE__ )
00558 return
00559
00560 end select
00561 end do
00562
00563 if (nb_neighbors < 1) then
00564 ierror = PRISM_Error_Interp_type
00565 ierrp (1) = nb_neighbors
00566
00567 call psmile_error ( ierror, &
00568 'Number of interpolation bases is less than 1', &
00569 ierrp, 1, __FILE__, __LINE__ )
00570 return
00571 endif
00572
00573 if (comp_info%size > 1) then
00574 if ( interpolation_methods(1) == PSMILe_conserv2D ) then
00575 global_cell_search = ANY (interpolation_search)
00576 else
00577 global_search = ANY (interpolation_search)
00578 endif
00579
00580 #ifdef DEBUG
00581 print *, "interpolation search", interpolation_search (:)
00582 print *, "global search", global_search
00583 print *, "global cell search", global_cell_search
00584 #endif
00585 else
00586 interpolation_search = .false.
00587 endif
00588
00589 #ifdef PRISM_ASSERTION
00590 do j = 1, n_intmethods
00591 select case (interpolation_methods(j))
00592
00593
00594
00595
00596
00597 case (PSMILe_nnghbr3D)
00598 if (n_intmethods > 1) then
00599 call psmile_assert (__FILE__, __LINE__, &
00600 "Don't know how to combine PSMILe_nnghbr3D with other interpolations")
00601 endif
00602
00603
00604
00605
00606
00607 case (PSMILe_trilinear)
00608 if (n_intmethods > 1) then
00609 call psmile_assert (__FILE__, __LINE__, &
00610 "Don't know how to combine PSMILe_trilinear with other interpolations")
00611 endif
00612
00613
00614
00615
00616 case (PSMILe_bilinear)
00617 if (n_intmethods > 2 .or. &
00618 (n_intmethods == 2 .and. &
00619 interpolation_methods(2) /= PSMILe_linear)) then
00620 call psmile_assert (__FILE__, __LINE__, &
00621 "Don't know how to combine PSMILe_bilinear with that interpolation")
00622 endif
00623
00624
00625
00626
00627 case (PSMILe_linear)
00628
00629 do k = 1, n_intmethods
00630 if (interpolation_methods(k) /= PSMILe_linear) exit
00631 end do
00632
00633 if (k <= n_intmethods) then
00634 if (n_intmethods < 2 .or. &
00635 (n_intmethods == 2 .and. &
00636 (interpolation_methods(1) /= PSMILe_bilinear .and. &
00637 interpolation_methods(1) /= PSMILe_bicubic) )) then
00638 call psmile_assert (__FILE__, __LINE__, &
00639 "Don't know how to combine PSMILe_linear with that interpolation")
00640 endif
00641 endif
00642
00643
00644
00645
00646
00647
00648 case (PSMILe_bicubic)
00649 if (n_intmethods > 2 .or. &
00650 (n_intmethods == 2 .and. &
00651 interpolation_methods(2) /= PSMILe_linear)) then
00652 call psmile_assert (__FILE__, __LINE__, &
00653 "Don't know how to combine PSMILe_bicubic with that interpolation")
00654 endif
00655
00656 end select
00657 end do
00658 #endif
00659
00660
00661
00662 Allocate (neighbors_3d(ndim_3d, ncpl, nb_neighbors), STAT = ierror)
00663 if ( ierror > 0 ) then
00664 ierrp (1) = ierror
00665 ierrp (2) = ndim_3d * ncpl * nb_neighbors
00666
00667 ierror = PRISM_Error_Alloc
00668 call psmile_error ( ierror, 'neighbors_3d', &
00669 ierrp, 2, __FILE__, __LINE__ )
00670 return
00671 endif
00672
00673
00674
00675
00676
00677
00678
00679
00680
00681 use_target_cell = ( interpolation_methods(1) == PSMILe_Conserv2D .or. &
00682 interpolation_methods(1) == PSMILe_Conserv3D )
00683
00684 use_source_cell = use_target_cell
00685
00686 if (use_target_cell) then
00687
00688
00689
00690 use_target_grid = .false.
00691
00692 if ( search%grid_type == PRISM_Irrlonlat_regvrt ) then
00693 tgt_corners(1:2) = 4
00694 tgt_corners(3) = 4
00695
00696 do i = 1, ndim_3d
00697 Allocate (tgt_cell(i)%vector(ncpl*tgt_corners(i)), stat = ierror)
00698 if ( ierror > 0 ) then
00699 ierrp (1) = ierror
00700 ierrp (2) = ncpl
00701
00702 ierror = PRISM_Error_Alloc
00703 call psmile_error ( ierror, 'tgt_cell(i)%vector', &
00704 ierrp, 2, __FILE__, __LINE__ )
00705 return
00706 endif
00707 end do
00708
00709 else if ( search%grid_type == PRISM_Reglonlatvrt ) then
00710
00711 tgt_corners(:) = 4
00712
00713 do i = 1, ndim_3d
00714 Allocate (tgt_cell(i)%vector(ncpl*tgt_corners(i)), stat = ierror)
00715 if ( ierror > 0 ) then
00716 ierrp (1) = ierror
00717 ierrp (2) = ncpl
00718
00719 ierror = PRISM_Error_Alloc
00720 call psmile_error ( ierror, 'tgt_cell(i)%vector', &
00721 ierrp, 2, __FILE__, __LINE__ )
00722 return
00723 endif
00724 end do
00725
00726 else if ( search%grid_type == PRISM_Gaussreduced_regvrt ) then
00727
00728 tgt_corners(:) = 4
00729
00730 nprev = 0
00731 ncpl = 0
00732
00733 do ipart = 1, search%npart
00734
00735 totsiz = size(coords(1,ipart)%vector)
00736 offset = totsiz/2
00737
00738 if (len_cpl (ipart) > 0) then
00739 ibeg = nprev + 1
00740 iend = nprev + len_cpl (ipart)
00741 do j = ibeg, iend
00742 if ( dstijk(1,j) /= PSMILe_Undef .and. &
00743 dstijk(2,j) /= PSMILe_Undef .and. &
00744 dstijk(1,j)-shape(1,1,ipart)+1 + offset <= totsiz .and. &
00745 dstijk(3,j) /= PSMILe_Undef ) ncpl = ncpl + 1
00746 enddo
00747
00748 nprev = nprev + len_cpl (ipart)
00749
00750 endif
00751
00752 enddo
00753
00754 send_info%nloc = ncpl
00755
00756 do i = 1, ndim_3d
00757 Allocate (tgt_cell(i)%vector(ncpl*tgt_corners(i)), stat = ierror)
00758 if ( ierror > 0 ) then
00759 ierrp (1) = ierror
00760 ierrp (2) = ncpl
00761
00762 ierror = PRISM_Error_Alloc
00763 call psmile_error ( ierror, 'tgt_cell(i)%vector', &
00764 ierrp, 2, __FILE__, __LINE__ )
00765 return
00766 endif
00767 end do
00768
00769 else
00770 ierrp (1) = search%grid_type
00771 ierror = PRISM_Error_Internal
00772 call psmile_error ( ierror, 'unsupported grid type', &
00773 ierrp, 1, __FILE__, __LINE__ )
00774 return
00775 endif
00776
00777 nprev = 0
00778 jj = 0
00779
00780 do ipart = 1, search%npart
00781
00782 if (len_cpl (ipart) > 0) then
00783
00784 ibeg = nprev + 1
00785 iend = nprev + len_cpl (ipart)
00786
00787 if ( search%grid_type == PRISM_Irrlonlat_regvrt ) then
00788
00789 len = shape(2,1,ipart)-shape(1,1,ipart)+1
00790
00791 do i = 1, ndim_2d
00792 do j = ibeg, iend
00793 if ( dstijk(1,j) /= PSMILe_Undef .and. dstijk(2,j) /= PSMILe_Undef ) then
00794 n = (dstijk(2,j)-shape(1,2,ipart))*len + &
00795 dstijk(1,j)-shape(1,1,ipart)+1
00796 tgt_cell(i)%vector(0*ncpl+j) = coords(i,ipart)%vector(n)
00797 tgt_cell(i)%vector(1*ncpl+j) = coords(i,ipart)%vector(n+1)
00798 tgt_cell(i)%vector(2*ncpl+j) = coords(i,ipart)%vector(n+len+1)
00799 tgt_cell(i)%vector(3*ncpl+j) = coords(i,ipart)%vector(n+len)
00800 else
00801 do nn = 0, nb_neighbors-1
00802 tgt_cell(i)%vector(nn*ncpl+j) = 0.0
00803 enddo
00804 endif
00805 enddo
00806 end do
00807
00808 do j = ibeg, iend
00809 if ( dstijk(3,j) /= PSMILe_Undef ) then
00810 n = dstijk(3,j)-shape(1,3,ipart)+1
00811 do nn = 0, nb_neighbors-1
00812 tgt_cell(3)%vector(nn*ncpl+j) = coords(3,ipart)%vector(n)
00813 enddo
00814 else
00815 do nn = 0, nb_neighbors-1
00816 tgt_cell(3)%vector(nn*ncpl+j) = 0.0
00817 enddo
00818 endif
00819 enddo
00820
00821 if ( search%msg_intersections%requires_conserv_remap == PSMILe_Conserv3D ) then
00822 ierrp (1) = search%grid_type
00823 ierror = PRISM_Error_Internal
00824 call psmile_error ( ierror, 'unsupported interpolation', &
00825 ierrp, 1, __FILE__, __LINE__ )
00826 return
00827 endif
00828
00829 else if ( search%grid_type == PRISM_Reglonlatvrt ) then
00830
00831 do j = ibeg, iend
00832 if ( dstijk(1,j) /= PSMILe_Undef .and. dstijk(2,j) /= PSMILe_Undef ) then
00833 tgt_cell(1)%vector(0*ncpl+j) = coords(1,ipart)%vector(dstijk(1,j)-shape(1,1,ipart)+1)
00834 tgt_cell(1)%vector(1*ncpl+j) = coords(1,ipart)%vector(dstijk(1,j)-shape(1,1,ipart)+2)
00835 tgt_cell(1)%vector(2*ncpl+j) = coords(1,ipart)%vector(dstijk(1,j)-shape(1,1,ipart)+2)
00836 tgt_cell(1)%vector(3*ncpl+j) = coords(1,ipart)%vector(dstijk(1,j)-shape(1,1,ipart)+1)
00837 tgt_cell(2)%vector(0*ncpl+j) = coords(2,ipart)%vector(dstijk(2,j)-shape(1,2,ipart)+1)
00838 tgt_cell(2)%vector(1*ncpl+j) = coords(2,ipart)%vector(dstijk(2,j)-shape(1,2,ipart)+1)
00839 tgt_cell(2)%vector(2*ncpl+j) = coords(2,ipart)%vector(dstijk(2,j)-shape(1,2,ipart)+2)
00840 tgt_cell(2)%vector(3*ncpl+j) = coords(2,ipart)%vector(dstijk(2,j)-shape(1,2,ipart)+2)
00841 else
00842 do nn = 0, nb_neighbors-1
00843 tgt_cell(1)%vector(nn*ncpl+j) = 0.0
00844 tgt_cell(2)%vector(nn*ncpl+j) = 0.0
00845 enddo
00846 endif
00847 end do
00848
00849 if ( search%msg_intersections%requires_conserv_remap == PSMILe_Conserv3D ) then
00850 do j = ibeg, iend
00851 if ( dstijk(3,j) /= PSMILe_Undef ) then
00852 tgt_cell(3)%vector(0*ncpl+j) = coords(3,ipart)%vector(dstijk(3,j)-shape(1,3,ipart)+1)
00853 tgt_cell(3)%vector(1*ncpl+j) = tgt_cell(3)%vector(0*ncpl+j)
00854 tgt_cell(3)%vector(2*ncpl+j) = tgt_cell(3)%vector(0*ncpl+j)
00855 tgt_cell(3)%vector(3*ncpl+j) = tgt_cell(3)%vector(0*ncpl+j)
00856 tgt_cell(3)%vector(4*ncpl+j) = coords(3,ipart)%vector(dstijk(3,j)-shape(1,3,ipart)+2)
00857 tgt_cell(3)%vector(5*ncpl+j) = tgt_cell(3)%vector(4*ncpl+j)
00858 tgt_cell(3)%vector(6*ncpl+j) = tgt_cell(3)%vector(4*ncpl+j)
00859 tgt_cell(3)%vector(7*ncpl+j) = tgt_cell(3)%vector(4*ncpl+j)
00860 else
00861 do nn = 0, nb_neighbors-1
00862 tgt_cell(3)%vector(nn*ncpl+j) = 0.0
00863 enddo
00864 endif
00865 end do
00866
00867 else
00868
00869 do j = ibeg, iend
00870 if ( dstijk(3,j) /= PSMILe_Undef ) then
00871 tgt_cell(3)%vector(0*ncpl+j) = coords(3,ipart)%vector(dstijk(3,j)-shape(1,3,ipart)+1)
00872 tgt_cell(3)%vector(1*ncpl+j) = tgt_cell(3)%vector(0*ncpl+j)
00873 tgt_cell(3)%vector(2*ncpl+j) = tgt_cell(3)%vector(0*ncpl+j)
00874 tgt_cell(3)%vector(3*ncpl+j) = tgt_cell(3)%vector(0*ncpl+j)
00875 else
00876 do nn = 0, nb_neighbors-1
00877 tgt_cell(3)%vector(nn*ncpl+j) = 0.0
00878 enddo
00879 endif
00880 end do
00881 endif
00882
00883 else if (search%grid_type == PRISM_Gaussreduced_regvrt) then
00884
00885 totsiz = size(coords(1,ipart)%vector)
00886 offset = totsiz/2
00887
00888 len_cpl(ipart) = 0
00889
00890 do j = ibeg, iend
00891
00892 if ( dstijk(1,j) /= PSMILe_Undef .and. &
00893 dstijk(2,j) /= PSMILe_Undef .and. &
00894 dstijk(1,j)-shape(1,1,ipart)+1 + offset <= totsiz .and. &
00895 dstijk(3,j) /= PSMILe_Undef ) then
00896 n = dstijk(1,j)-shape(1,1,ipart)+1
00897
00898 jj = jj + 1
00899 len_cpl(ipart) = len_cpl(ipart) + 1
00900
00901 tgt_cell(1)%vector(0*ncpl+jj) = coords(1,ipart)%vector(n)
00902 tgt_cell(1)%vector(1*ncpl+jj) = coords(1,ipart)%vector(n+offset)
00903 tgt_cell(1)%vector(2*ncpl+jj) = coords(1,ipart)%vector(n+offset)
00904 tgt_cell(1)%vector(3*ncpl+jj) = coords(1,ipart)%vector(n)
00905 tgt_cell(2)%vector(0*ncpl+jj) = coords(2,ipart)%vector(n)
00906 tgt_cell(2)%vector(1*ncpl+jj) = coords(2,ipart)%vector(n)
00907 tgt_cell(2)%vector(2*ncpl+jj) = coords(2,ipart)%vector(n+offset)
00908 tgt_cell(2)%vector(3*ncpl+jj) = coords(2,ipart)%vector(n+offset)
00909
00910 n = dstijk(3,j)-shape(1,3,ipart)+1
00911 do nn = 0, nb_neighbors-1
00912 tgt_cell(3)%vector(nn*ncpl+jj) = coords(3,ipart)%vector(n)
00913 enddo
00914 else
00915 send_info%msklocs(1,ipart)%vector(j-ibeg+1) = .false.
00916 endif
00917
00918 enddo
00919
00920 if ( search%msg_intersections%requires_conserv_remap == PSMILe_Conserv3D ) then
00921 ierrp (1) = search%grid_type
00922 ierror = PRISM_Error_Internal
00923 call psmile_error ( ierror, 'unsupported interpolation', &
00924 ierrp, 1, __FILE__, __LINE__ )
00925 return
00926 endif
00927
00928 else
00929
00930 ierrp (1) = search%grid_type
00931 ierror = PRISM_Error_Internal
00932 call psmile_error ( ierror, 'unsupported grid type', &
00933 ierrp, 1, __FILE__, __LINE__ )
00934 return
00935
00936 endif
00937
00938 nprev = nprev + len_cpl (ipart)
00939
00940 endif
00941
00942 end do
00943
00944 else
00945
00946
00947
00948 use_target_grid = .not. ( &
00949 search%grid_type == PRISM_Irrlonlatvrt .and. &
00950 search%npart == 1 .and. &
00951 ncpl == (shape(2,1,1)-shape(1,1,1)+1)* &
00952 (shape(2,2,1)-shape(1,2,1)+1)* &
00953 (shape(2,3,1)-shape(1,3,1)+1) )
00954
00955 if (use_target_grid) then
00956
00957
00958
00959 Allocate (neigh_bascule(ncpl), STAT = ierror)
00960 if ( ierror > 0 ) then
00961 ierrp (1) = ierror
00962 ierrp (2) = ncpl
00963
00964 ierror = PRISM_Error_Alloc
00965 call psmile_error ( ierror, 'neigh_bascule', &
00966 ierrp, 2, __FILE__, __LINE__ )
00967 return
00968 endif
00969
00970 neigh_bascule = 0
00971
00972
00973
00974 do i = 1, ndim_3d
00975 Allocate (tgt_grid(i)%vector(ncpl), stat = ierror)
00976 if ( ierror > 0 ) then
00977 ierrp (1) = ierror
00978 ierrp (2) = ncpl
00979
00980 ierror = PRISM_Error_Alloc
00981 call psmile_error ( ierror, 'tgt_grid(i)%vector', &
00982 ierrp, 2, __FILE__, __LINE__ )
00983 return
00984 endif
00985 end do
00986
00987
00988
00989 if (search%grid_type == PRISM_Reglonlatvrt) then
00990
00991
00992
00993 nprev = 0
00994 do ipart = 1, search%npart
00995 if (len_cpl (ipart) > 0) then
00996 ibeg = nprev + 1
00997 iend = nprev + len_cpl (ipart)
00998
00999 do i = 1, ndim_3d
01000 tgt_grid(i)%vector(ibeg:iend) = &
01001 coords(i,ipart)%vector(dstijk(i,ibeg:iend)- &
01002 shape(1,i,ipart)+1)
01003 end do
01004
01005 #ifdef DEBUGX
01006 print *, 'ibeg, iend', ibeg, iend
01007 do i = ibeg, iend
01008 print *, 'i target coord', i, tgt_grid(1)%vector(i), &
01009 tgt_grid(2)%vector(i), &
01010 tgt_grid(3)%vector(i)
01011 enddo
01012 #endif
01013
01014 nprev = nprev + len_cpl (ipart)
01015 endif
01016 end do
01017
01018 else if (search%grid_type == PRISM_Irrlonlat_Regvrt) then
01019
01020
01021
01022 nprev = 0
01023 do ipart = 1, search%npart
01024 if (len_cpl (ipart) > 0) then
01025 ibeg = nprev + 1
01026 iend = nprev + len_cpl (ipart)
01027
01028 call psmile_extract_indices_2d_real ( &
01029 coords(1,ipart)%vector, shape (:,1:ndim_2d, ipart), &
01030 dstijk (:, ibeg:iend), len_cpl(ipart), &
01031 tgt_grid(1)%vector(ibeg:iend), ierror)
01032 if (ierror /= 0) return
01033
01034 call psmile_extract_indices_2d_real ( &
01035 coords(2,ipart)%vector, shape (:,1:ndim_2d, ipart), &
01036 dstijk (:, ibeg:iend), len_cpl(ipart), &
01037 tgt_grid(2)%vector(ibeg:iend), ierror)
01038 if (ierror /= 0) return
01039
01040 tgt_grid(3)%vector(ibeg:iend) = &
01041 coords(3,ipart)%vector(dstijk(3,ibeg:iend)-shape(1,3,ipart)+1)
01042 #ifdef DEBUGX
01043 print *, 'ibeg, iend', ibeg, iend
01044 do i = ibeg, iend
01045 print *, 'i target coord', i, tgt_grid(1)%vector(i), &
01046 tgt_grid(2)%vector(i), &
01047 tgt_grid(3)%vector(i)
01048 end do
01049 #endif
01050 nprev = nprev + len_cpl (ipart)
01051 endif
01052 end do
01053
01054 else if (search%grid_type == PRISM_Gaussreduced_Regvrt) then
01055
01056
01057
01058 nprev = 0
01059 do ipart = 1, search%npart
01060 if (len_cpl (ipart) > 0) then
01061 ibeg = nprev + 1
01062 iend = nprev + len_cpl (ipart)
01063
01064 call psmile_extract_indices_2d_real ( &
01065 coords(1,ipart)%vector, shape (:,1:ndim_2d, ipart), &
01066 dstijk (:, ibeg:iend), len_cpl(ipart), &
01067 tgt_grid(1)%vector(ibeg:iend), ierror)
01068 if (ierror /= 0) return
01069
01070 call psmile_extract_indices_2d_real ( &
01071 coords(2,ipart)%vector, shape (:,1:ndim_2d, ipart), &
01072 dstijk (:, ibeg:iend), len_cpl(ipart), &
01073 tgt_grid(2)%vector(ibeg:iend), ierror)
01074 if (ierror /= 0) return
01075
01076 tgt_grid(3)%vector(ibeg:iend) = &
01077 coords(3,ipart)%vector(dstijk(3,ibeg:iend)-shape(1,3,ipart)+1)
01078
01079
01080 nprev = nprev + len_cpl (ipart)
01081 endif
01082 end do
01083
01084 else if (search%grid_type == PRISM_Irrlonlatvrt) then
01085
01086
01087
01088 nprev = 0
01089 do ipart = 1, search%npart
01090 if (len_cpl (ipart) > 0) then
01091 ibeg = nprev + 1
01092 iend = nprev + len_cpl (ipart)
01093
01094 do i = 1, ndim_3d
01095 call psmile_extract_indices_3d_real ( &
01096 coords(i,ipart)%vector, shape(:,:, ipart), &
01097 dstijk (:, ibeg:iend), len_cpl(ipart), &
01098 tgt_grid(i)%vector(ibeg:iend), &
01099 ierror)
01100 if (ierror > 0) return
01101 enddo
01102
01103 nprev = nprev + len_cpl (ipart)
01104 endif
01105 end do
01106
01107 else
01108
01109
01110
01111 ierrp (1) = search%grid_type
01112 ierror = PRISM_Error_Internal
01113 call psmile_error ( ierror, 'unsupported grid type', &
01114 ierrp, 1, __FILE__, __LINE__ )
01115 return
01116
01117 endif
01118
01119 else
01120
01121
01122
01123 do i = 1, ndim_3d
01124 tgt_grid(i)%vector => coords(i,1)%vector
01125 end do
01126
01127 endif
01128
01129 endif
01130
01131 #ifdef DEBUG_TRACE
01132 if ( use_target_grid ) then
01133 ictl = 1
01134
01135 do ictl = 1, ncpl
01136 if ( dstijk (1, ictl) == ictl_ind (1) .and. &
01137 dstijk (2, ictl) == ictl_ind (2) .and. &
01138 dstijk (3, ictl) == ictl_ind (3)) exit
01139 end do
01140
01141 if (ictl <= ncpl) then
01142 print *, 'ictl target coord', ictl, tgt_grid(1)%vector(ictl), &
01143 tgt_grid(2)%vector(ictl), &
01144 tgt_grid(3)%vector(ictl)
01145 else
01146 print *, 'psmile_info_trs_gauss2_real: ictl_ind not found', ictl_ind (1:3)
01147 do ictl = 1, ncpl
01148 print *, 'dstijk', ictl, dstijk (:, ictl)
01149 end do
01150 endif
01151
01152 endif
01153 #endif
01154
01155
01156
01157
01158
01159
01160
01161 virtual_cell_available = Associated (send_info%virtual)
01162
01163 extra_search%global_search = global_search
01164
01165 call psmile_neigh_extra_search_init (search, grid_id, extra_search, ierror)
01166 if (ierror > 0) return
01167
01168
01169
01170 search_mode = 3
01171
01172 j = 1
01173 do while (j <= n_intmethods)
01174
01175 select case (interpolation_methods(j))
01176
01177
01178
01179
01180
01181 case (PSMILe_nnghbr3D,PSMILe_nnghbr2D)
01182
01183
01184
01185
01186 coords_shape(1:2,1) = mp%coords_pointer%coords_shape(1:2,1)
01187 coords_shape(1:2,2) = mp%coords_pointer%coords_shape(1:2,1)
01188 coords_shape(1:2,3) = mp%coords_pointer%coords_shape(1:2,3)
01189
01190 call psmile_info_coords_3d_reg_real ( &
01191 mp%coords_pointer%coords_real(1)%vector, &
01192 mp%coords_pointer%coords_real(2)%vector, &
01193 mp%coords_pointer%coords_real(3)%vector, &
01194 coords_shape, &
01195 grid_valid_shape, &
01196 sinvec, cosvec, ierror)
01197 if (ierror > 0) return
01198
01199 transformed = .true.
01200
01201
01202
01203 nprev = 0
01204
01205 if (send_info%nvec == 1) then
01206
01207 do ipart = 1, search%npart
01208 if (len_cpl (ipart) > 0) then
01209
01210 call psmile_neigh_near_3d_irr3_real (grid_id, &
01211 tgt_grid(1)%vector, &
01212 tgt_grid(2)%vector, &
01213 tgt_grid(3)%vector, &
01214 mp%coords_pointer%coords_real(1)%vector, &
01215 mp%coords_pointer%coords_real(2)%vector, &
01216 mp%coords_pointer%coords_real(3)%vector, &
01217 coords_shape, &
01218 mask_array, mask_shape, use_mask, &
01219 sinvec(lon)%vector, cosvec(lon)%vector, &
01220 sinvec(lat)%vector, cosvec(lat)%vector, &
01221 grid_valid_shape, search_mode, &
01222 send_info%srclocs(1,1)%vector, &
01223 len_cpl (ipart), ncpl, nprev, &
01224 neighbors_3d, nb_neighbors, &
01225 extra_search, ierror)
01226 if (ierror > 0) return
01227
01228 nprev = nprev + len_cpl (ipart)
01229 end if
01230 end do
01231
01232 else if (send_info%nvec == 2) then
01233
01234 do ipart = 1, search%npart
01235 if (len_cpl (ipart) > 0) then
01236
01237 call psmile_neigh_near_3d_irr2_real (grid_id, &
01238 tgt_grid(1)%vector(nprev+1:), &
01239 tgt_grid(2)%vector(nprev+1:), &
01240 tgt_grid(3)%vector(nprev+1:), &
01241 mp%coords_pointer%coords_real(1)%vector, &
01242 mp%coords_pointer%coords_real(2)%vector, &
01243 mp%coords_pointer%coords_real(3)%vector, &
01244 coords_shape, &
01245 mask_array, mask_shape, use_mask, &
01246 sinvec(lon)%vector, cosvec(lon)%vector, &
01247 sinvec(lat)%vector, cosvec(lat)%vector, &
01248 grid_valid_shape, search_mode, &
01249 send_info%srclocs(indl, ipart)%vector, &
01250 send_info%srclocs(indz, ipart)%vector, &
01251 send_info%npoints(1:indz, ipart), &
01252 ncpl, nprev, &
01253 neighbors_3d, nb_neighbors, &
01254 extra_search, ierror)
01255
01256 if (ierror > 0) return
01257
01258 nprev = nprev + len_cpl (ipart)
01259 end if
01260 end do
01261
01262 else if (send_info%nvec == ndim_3d) then
01263
01264
01265 do ipart = 1, search%npart
01266 if (len_cpl (ipart) > 0) then
01267
01268 call psmile_neigh_near_3d_reg_real (grid_id, &
01269 tgt_grid(1)%vector(nprev+1:), &
01270 tgt_grid(2)%vector(nprev+1:), &
01271 tgt_grid(3)%vector(nprev+1:), &
01272 mp%coords_pointer%coords_real(1)%vector, &
01273 mp%coords_pointer%coords_real(2)%vector, &
01274 mp%coords_pointer%coords_real(3)%vector, &
01275 coords_shape, &
01276 mask_array, mask_shape, use_mask, &
01277 sinvec(lon)%vector, cosvec(lon)%vector, &
01278 sinvec(lat)%vector, cosvec(lat)%vector, &
01279 grid_valid_shape, search_mode, &
01280 send_info%srclocs(1:ndim_3d, ipart), &
01281 send_info%npoints(1:ndim_3d, ipart), &
01282 ncpl, nprev, &
01283 neighbors_3d, nb_neighbors, &
01284 extra_search, ierror)
01285
01286 if (ierror > 0) return
01287
01288 nprev = nprev + len_cpl (ipart)
01289 end if
01290 end do
01291
01292 else
01293
01294
01295
01296 ierrp (1) = search%grid_type
01297 ierror = PRISM_Error_Internal
01298 call psmile_error ( ierror, 'unsupported search grid type', &
01299 ierrp, 1, __FILE__, __LINE__ )
01300 return
01301
01302 endif
01303
01304
01305
01306
01307
01308 case (PSMILe_trilinear, PSMILe_bilinear, PSMILe_bicubic, PSMILe_linear)
01309
01310
01311 interpolation_mode = 3
01312
01313 if (interpolation_methods(j) == PSMILe_bilinear .or. &
01314 interpolation_methods(j) == PSMILe_bicubic ) then
01315 if (interpolation_methods(j+1) == PSMILe_linear) then
01316 j = j + 1
01317 else
01318
01319 interpolation_mode = 2
01320
01321 search_mode = 2
01322 endif
01323 else if (interpolation_methods(j) == PSMILe_linear) then
01324
01325 do k = j+1, n_intmethods
01326 if (interpolation_methods(k) /= PSMILe_linear) exit
01327 end do
01328
01329 interpolation_mode = k - j
01330 j = k - 1
01331 endif
01332
01333 if (send_info%nvec == 1) then
01334
01335
01336 if (virtual_cell_available) then
01337 virtual_cell => send_info%virtual(1,1)%vector
01338 else
01339 virtual_cell => dummy_virtual_cell
01340 endif
01341
01342 if ( cubic ) then
01343 call psmile_neigh_tricu_gauss2 ( &
01344 mp%grid_id, grid_valid_shape, &
01345 interpolation_mode, &
01346 send_info%srclocs(1,1)%vector, &
01347 virtual_cell, ncpl, &
01348 virtual_cell_available, &
01349 neighbors_3d, nb_neighbors, &
01350 neigh_bascule, ierror)
01351 else
01352 call psmile_neigh_trili_gauss2 ( &
01353 mp%grid_id, grid_valid_shape, &
01354 interpolation_mode, &
01355 send_info%srclocs(1,1)%vector, &
01356 virtual_cell, ncpl, &
01357 virtual_cell_available, &
01358 neighbors_3d, nb_neighbors, ierror)
01359 endif
01360 if (ierror /= 0) return
01361
01362
01363 else if (send_info%nvec == 2) then
01364
01365
01366
01367 nprev = 0
01368
01369 do ipart = 1, search%npart
01370 if (len_cpl (ipart) > 0) then
01371
01372 if (virtual_cell_available) then
01373 virtual_cell => send_info%virtual(1,ipart)%vector
01374 else
01375 virtual_cell => dummy_virtual_cell
01376 endif
01377
01378 if ( cubic ) then
01379 call psmile_neigh_tricu_gauss2_irreg (mp%grid_id, &
01380 grid_valid_shape, &
01381 interpolation_mode, &
01382 send_info%srclocs(indl, ipart)%vector, &
01383 send_info%srclocs(indz, ipart)%vector, &
01384 virtual_cell, &
01385 send_info%npoints(1:indz, ipart), &
01386 ncpl, nprev, &
01387 virtual_cell_available, &
01388 neighbors_3d, nb_neighbors, neigh_bascule, ierror)
01389 else
01390 call psmile_neigh_trili_gauss2_irreg (mp%grid_id, &
01391 grid_valid_shape, &
01392 interpolation_mode, &
01393 send_info%srclocs(indl, ipart)%vector, &
01394 send_info%srclocs(indz, ipart)%vector, &
01395 virtual_cell, &
01396 send_info%npoints(1:indz, ipart), &
01397 ncpl, nprev, &
01398 virtual_cell_available, &
01399 neighbors_3d, nb_neighbors, ierror)
01400 endif
01401
01402 if (ierror /= 0) return
01403
01404 nprev = nprev + len_cpl (ipart)
01405 end if
01406 end do
01407
01408 else
01409
01410 nprev = 0
01411
01412 do ipart = 1, search%npart
01413 if (len_cpl (ipart) > 0) then
01414
01415 if ( cubic ) then
01416 call psmile_neigh_tricu_3d_reg ( &
01417 grid_valid_shape, &
01418 interpolation_mode, grid_info%cyclic, &
01419 send_info%srclocs(1:ndim_3d, ipart), &
01420 send_info%npoints(1:ndim_3d, ipart), &
01421 ncpl, nprev, &
01422 neighbors_3d, nb_neighbors, ierror)
01423 else
01424 call psmile_neigh_trili_3d_reg ( &
01425 grid_valid_shape, &
01426 interpolation_mode, grid_info%cyclic, &
01427 send_info%srclocs(1:ndim_3d, ipart), &
01428 send_info%npoints(1:ndim_3d, ipart), &
01429 ncpl, nprev, &
01430 neighbors_3d, nb_neighbors, ierror)
01431 endif
01432 if (ierror /= 0) return
01433
01434 nprev = nprev + len_cpl (ipart)
01435 end if
01436 end do
01437 endif
01438
01439
01440
01441
01442
01443 case (PSMILe_conserv2D)
01444 corner_shape_3d = 1
01445
01446 if ( grid_info%grid_type == PRISM_Gaussreduced_regvrt ) then
01447 do i = 1, ndim_2d
01448 corner_shape_3d(:,i,i) = cp%corner_shape(:,1)
01449 enddo
01450 corner_shape_3d(:,ndim_3d,ndim_3d) = cp%corner_shape(:,3)
01451 else
01452 do i = 1, ndim_3d
01453 corner_shape_3d(:,i,i) = cp%corner_shape(:,i)
01454 enddo
01455 endif
01456
01457 Allocate(nbr_cells(ncpl), STAT = ierror)
01458
01459 if ( ierror > 0 ) then
01460 ierrp (1) = ierror
01461 ierrp (2) = ncpl
01462 ierror = PRISM_Error_Alloc
01463 call psmile_error ( ierror, 'nbr_cells', &
01464 ierrp, 2, __FILE__, __LINE__ )
01465 return
01466 endif
01467
01468
01469 interpolation_mode = 3
01470
01471 if (interpolation_methods(j+1) == PSMILe_linear) then
01472 j = j + 1
01473 else
01474
01475 interpolation_mode = 2
01476 search_mode = 2
01477 endif
01478
01479 ipart = search%npart
01480
01481 if (send_info%nvec == 1) then
01482
01483
01484 nbr_corners = 2
01485
01486 call psmile_neigh_cells_3d_real ( use_how, &
01487 grid_valid_shape, &
01488 interpolation_mode, grid_info%cyclic, &
01489 corner_shape_3d, nbr_corners, &
01490 cp%corners_real(1)%vector, &
01491 cp%corners_real(2)%vector, &
01492 cp%corners_real(3)%vector, &
01493 search, control, tgt_cell, tgt_corners, &
01494 send_info%npoints(1,1), &
01495 send_info%srclocs(1,1)%vector, &
01496 send_info%msklocs(1,1)%vector, &
01497 ncpl, nb_neighbors, nbr_cells, ierror )
01498 if (ierror /= 0) return
01499
01500 else if (send_info%nvec == 2) then
01501
01502
01503
01504 nbr_corners = 2
01505
01506 call psmile_neigh_cells_irreg2_real ( use_how, &
01507 grid_valid_shape, interpolation_mode, &
01508 grid_info%cyclic, &
01509 corner_shape_3d, nbr_corners, &
01510 cp%corners_real(1)%vector, &
01511 cp%corners_real(2)%vector, &
01512 search, control, tgt_cell, tgt_corners(1), &
01513 send_info%npoints(1:indz, 1:ipart), &
01514 send_info%srclocs(1:indz, 1:ipart), &
01515 send_info%msklocs(1:indz, 1:ipart), &
01516 ncpl, nb_neighbors, nbr_cells, ierror )
01517 if (ierror /= 0) return
01518
01519 else if (send_info%nvec == ndim_3d) then
01520
01521 if (send_info%nloc > 0) then
01522 call psmile_neigh_cells_3d_reg_real ( &
01523 grid_valid_shape, &
01524 interpolation_mode, grid_info%cyclic, &
01525 grid_id, search, coords, &
01526 send_info%npoints(1:ndim_3d, 1:ipart), &
01527 send_info%srclocs(1:ndim_3d, 1:ipart), &
01528 ncpl, nbr_cells, ierror )
01529 if (ierror /= 0) return
01530 end if
01531
01532 endif
01533
01534
01535
01536 if ( associated(search%boundary_cell) .and. global_cell_search ) then
01537
01538 nprev = 0
01539 nn = 1
01540
01541 len = size(search%boundary_cell(:,1))
01542
01543 do ipart = 1, search%npart
01544
01545 if (len_cpl (ipart) > 0) then
01546
01547 ibeg = nprev + 1
01548 iend = nprev + len_cpl (ipart)
01549
01550 if ( search%grid_type == PRISM_Irrlonlat_regvrt ) then
01551
01552 leni = shape(2,1,ipart)-shape(1,1,ipart)+1
01553
01554
01555
01556
01557
01558
01559 do jpart = ibeg, iend, 5000
01560 do j = jpart, min(iend,jpart+5000-1)
01561
01562 if ( dstijk(1,j) /= PSMILe_Undef .and. &
01563 dstijk(2,j) /= PSMILe_Undef .and. &
01564 nn <= len ) then
01565
01566 n = (dstijk(2,j)-shape(1,2,ipart))*leni + &
01567 dstijk(1,j)-shape(1,1,ipart)+1
01568
01569 if ( j == search%boundary_cell(nn,1) ) then
01570 n2 = ( n-1) / leni + 1
01571 n1 = n - (n2-1)*leni
01572
01573 search%boundary_cell(nn,2) = &
01574 n1 + search%global_index(ipart)%vector(1)
01575 search%boundary_cell(nn,3) = &
01576 n2 + search%global_index(ipart)%vector(3)
01577 search%boundary_cell(nn,4) = &
01578 dstijk(3,j)-shape(1,3,ipart)+1 + &
01579 search%global_index(ipart)%vector(5)
01580
01581 nn = nn + 1
01582 endif
01583
01584 endif
01585 enddo
01586 enddo
01587
01588 else if ( search%grid_type == PRISM_Reglonlatvrt ) then
01589
01590
01591
01592
01593
01594
01595 do jpart = ibeg, iend, 5000
01596 do j = jpart, min(iend,jpart+5000-1)
01597
01598 if ( dstijk(1,j) /= PSMILe_Undef .and. &
01599 dstijk(2,j) /= PSMILe_Undef .and. &
01600 nn <= len ) then
01601
01602 if ( j == search%boundary_cell(nn,1) ) then
01603
01604 search%boundary_cell(nn,2) = dstijk(1,j)-shape(1,1,ipart)+1 + &
01605 search%global_index(ipart)%vector(1)
01606 search%boundary_cell(nn,3) = dstijk(2,j)-shape(1,2,ipart)+1 + &
01607 search%global_index(ipart)%vector(3)
01608 search%boundary_cell(nn,4) = dstijk(3,j)-shape(1,3,ipart)+1 + &
01609 search%global_index(ipart)%vector(5)
01610 nn = nn + 1
01611 endif
01612
01613 endif
01614 enddo
01615 enddo
01616
01617 else if ( search%grid_type == PRISM_Gaussreduced_regvrt ) then
01618
01619
01620
01621
01622
01623
01624 do jpart = ibeg, iend, 5000
01625 do j = jpart, min(iend,jpart+5000-1)
01626
01627 if ( dstijk(1,j) /= PSMILe_Undef .and. &
01628 dstijk(2,j) /= PSMILe_Undef .and. &
01629 nn <= len ) then
01630
01631 if ( j == search%boundary_cell(nn,1) ) then
01632
01633 search%boundary_cell(nn,2) = dstijk(1,j)-shape(1,1,ipart)+1 + &
01634 search%global_index(ipart)%vector(1)
01635 search%boundary_cell(nn,3) = 1
01636 search%boundary_cell(nn,4) = dstijk(3,j)-shape(1,3,ipart)+1 + &
01637 search%global_index(ipart)%vector(5)
01638 nn = nn + 1
01639 endif
01640
01641 endif
01642 enddo
01643 enddo
01644
01645 else
01646
01647
01648
01649
01650 ierrp (1) = search%grid_type
01651 ierror = PRISM_Error_Internal
01652 call psmile_error ( ierror, 'unsupported grid type', &
01653 ierrp, 1, __FILE__, __LINE__ )
01654 return
01655
01656 endif
01657
01658 nprev = iend
01659
01660 endif
01661
01662 enddo
01663
01664 endif
01665
01666 if ( global_cell_search ) &
01667 call psmile_global_search_cell_real ( grid_id, var_id, &
01668 comp_info, send_info, search, extra_search, ncpl, &
01669 nbr_cells, n_intmethods, interpolation_methods, &
01670 interpolation_search, ierror )
01671
01672 if ( associated(search%boundary_cell) .and. global_cell_search ) then
01673
01674
01675
01676
01677 new_len_cpl = len_cpl
01678 new_ncpl = send_info%nloc
01679
01680 if ( new_ncpl < ncpl ) then
01681
01682
01683
01684 allocate(new_nbr_cells(new_ncpl), new_dstijk(3,new_ncpl), stat = ierror)
01685 if ( ierror > 0 ) then
01686 ierrp (1) = ierror
01687 ierrp (2) = send_info%nloc
01688
01689 ierror = PRISM_Error_Alloc
01690 call psmile_error ( ierror, 'new_nbr_cell, new_dstijk', &
01691 ierrp, 2, __FILE__, __LINE__ )
01692 return
01693 endif
01694
01695 do i = 1, ndim_3d
01696 Allocate (new_tgt_cell(i)%vector(new_ncpl*tgt_corners(i)), stat = ierror)
01697 if ( ierror > 0 ) then
01698 ierrp (1) = ierror
01699 ierrp (2) = send_info%nloc*tgt_corners(i)
01700
01701 ierror = PRISM_Error_Alloc
01702 call psmile_error ( ierror, 'new_tgt_cell(i)%vector', &
01703 ierrp, 2, __FILE__, __LINE__ )
01704 return
01705 endif
01706 end do
01707
01708 nn = 1
01709
01710 ipart = 1
01711
01712 ibeg = 1
01713 iend = len_cpl (ipart)
01714
01715 nprev = iend
01716
01717
01718
01719 do n = 1, ncpl
01720
01721 if ( n > iend ) then
01722 ipart = ipart + 1
01723 ibeg = nprev + 1
01724 iend = nprev + len_cpl (ipart)
01725 nprev = iend
01726 endif
01727
01728
01729
01730
01731
01732 if ( nbr_cells(n) > 0 ) then
01733 new_nbr_cells(nn) = nbr_cells(n)
01734 new_dstijk(:,nn) = dstijk(:,n)
01735 nn = nn + 1
01736 else
01737 new_len_cpl(ipart) = new_len_cpl(ipart) - 1
01738 endif
01739
01740 enddo
01741
01742 len_cpl = new_len_cpl
01743
01744 do i = 1, ndim_2d
01745 nn = 1
01746 do n = 1, ncpl
01747 if ( nbr_cells(n) > 0 ) then
01748 #ifdef DEBUGX
01749 print *, ' - shift tgt cell idx ', n , ' -> ', nn
01750 #endif
01751 new_tgt_cell(i)%vector(0*new_ncpl+nn) = tgt_cell(i)%vector(0*ncpl+n)
01752 new_tgt_cell(i)%vector(1*new_ncpl+nn) = tgt_cell(i)%vector(1*ncpl+n)
01753 new_tgt_cell(i)%vector(2*new_ncpl+nn) = tgt_cell(i)%vector(2*ncpl+n)
01754 new_tgt_cell(i)%vector(3*new_ncpl+nn) = tgt_cell(i)%vector(3*ncpl+n)
01755 nn = nn + 1
01756 endif
01757 enddo
01758 enddo
01759
01760 do n = 1, ncpl
01761 nn = 1
01762 if ( nbr_cells(n) > 0 ) then
01763 do n3 = 0, nb_neighbors-1
01764 new_tgt_cell(3)%vector(n3*new_ncpl+nn) = tgt_cell(3)%vector(n3*ncpl+n)
01765 enddo
01766 nn = nn + 1
01767 endif
01768 enddo
01769
01770
01771
01772 deallocate(nbr_cells, send_info%dstijk)
01773 Nullify(dstijk)
01774
01775 do i = 1, ndim_3d
01776 deallocate(tgt_cell(i)%vector)
01777 enddo
01778
01779 Allocate(send_info%dstijk(3,new_ncpl), stat = ierror)
01780 if ( ierror > 0 ) then
01781 ierrp (1) = ierror
01782 ierrp (2) = send_info%nloc
01783
01784 ierror = PRISM_Error_Alloc
01785 call psmile_error ( ierror, 'send_info%dstijk', &
01786 ierrp, 2, __FILE__, __LINE__ )
01787 return
01788 endif
01789
01790 nbr_cells => new_nbr_cells
01791
01792 tgt_cell(1)%vector => new_tgt_cell(1)%vector
01793 tgt_cell(2)%vector => new_tgt_cell(2)%vector
01794 tgt_cell(3)%vector => new_tgt_cell(3)%vector
01795
01796 send_info%dstijk = new_dstijk
01797
01798 dstijk => send_info%dstijk
01799
01800 Deallocate(new_dstijk)
01801
01802 ncpl = send_info%nloc
01803
01804 endif
01805
01806 endif
01807
01808
01809
01810 #ifdef TODO
01811 case (PSMILe_conserv3D)
01812 call psmile_trs_give_
01813 if (ierror /= 0) return
01814
01815 case (PSMILe_user3D)
01816 call psmile_trs_give_
01817 if (ierror /= 0) return
01818 #endif
01819
01820
01821
01822 case default
01823 ierrp (1) = interpolation_methods (j)
01824 ierror = PRISM_Error_Internal
01825
01826 call psmile_error ( ierror, 'unsupported 3d interpolation method', &
01827 ierrp, 1, __FILE__, __LINE__ )
01828 return
01829 end select
01830
01831
01832
01833 j = j + 1
01834 end do
01835
01836 #ifdef DEBUG_TRACE
01837 if (ictl > 0 .and. ictl <= ncpl) then
01838 print *, 'before global: ictl, dstijk (:, i)', ictl, dstijk (:, ictl)
01839 do i = 1, nb_neighbors
01840 print *, 'neighbors_3d (:, ictl)', neighbors_3d (:, ictl, i)
01841 end do
01842
01843 if (virtual_cell_available) then
01844
01845 virtual_cell => send_info%virtual(1,1)%vector
01846 print *, 'virtual cell info', ictl, virtual_cell(ictl)
01847 endif
01848
01849 endif
01850 #endif
01851
01852
01853
01854 if (global_search) then
01855 if (interpolation_methods(1) /= PSMILe_nnghbr3D) then
01856
01857
01858
01859 call psmile_neigh_global_points (search, extra_search, &
01860 mask_array, mask_shape, src_mask_available, use_how(1), &
01861 grid_id, grid_valid_shape, &
01862 neighbors_3d, ncpl, nb_neighbors, len_cpl, ierror)
01863 if (ierror > 0) return
01864
01865 if (extra_search%n_req > 0) then
01866 call psmile_global_search_real (comp_info, &
01867 control, len_cpl, &
01868 var_id, grid_valid_shape, search, tgt_grid, &
01869 neighbors_3d, ncpl, nb_neighbors, extra_search, &
01870 interpolation_methods, interpolation_search, n_intmethods, &
01871 send_index, src_mask_available, use_mask, &
01872 use_how(1), PRISM_Gaussreduced_regvrt, ierror)
01873 if (ierror > 0) return
01874 endif
01875
01876 else
01877
01878
01879
01880 endif
01881 endif
01882
01883
01884
01885
01886
01887 if (interpolation_methods(1) /= PSMILe_nnghbr3D) then
01888 call psmile_neigh_extra_points (search, extra_search, &
01889 mask_array, mask_shape, src_mask_available, use_how, &
01890 grid_valid_shape, neighbors_3d, ncpl, nb_neighbors, &
01891 len_cpl, ierror)
01892 if (ierror > 0) return
01893 endif
01894
01895
01896
01897 if (extra_search%n_extra > 0) then
01898
01899 grid_shape(1:2,1) = grid_valid_shape(1:2,1)
01900 grid_shape(1:2,2) = grid_valid_shape(1:2,1)
01901 grid_shape(1:2,3) = grid_valid_shape(1:2,3)
01902
01903
01904
01905
01906 if (global_search) then
01907 call psmile_neigh_extra_search_real (search, extra_search, &
01908 nb_extra, ierror)
01909 if (ierror > 0) return
01910 endif
01911
01912
01913 nprev = 0
01914
01915
01916
01917 if (.not. transformed) then
01918
01919 coords_shape(1:2,1) = mp%coords_pointer%coords_shape(1:2,1)
01920 coords_shape(1:2,2) = mp%coords_pointer%coords_shape(1:2,1)
01921 coords_shape(1:2,3) = mp%coords_pointer%coords_shape(1:2,3)
01922
01923 call psmile_info_coords_3d_reg_real ( &
01924 mp%coords_pointer%coords_real(1)%vector, &
01925 mp%coords_pointer%coords_real(2)%vector, &
01926 mp%coords_pointer%coords_real(3)%vector, &
01927 coords_shape, &
01928 grid_shape, &
01929 sinvec, cosvec, ierror)
01930 if (ierror > 0) return
01931
01932 transformed = .true.
01933 endif
01934
01935
01936
01937 if (send_info%nvec == 1) then
01938
01939 do ipart = 1, search%npart
01940 if (len_cpl (ipart) > 0) then
01941 call psmile_neigh_nearx_3d_irr3_real (grid_id, &
01942 tgt_grid(1)%vector(nprev+1:), &
01943 tgt_grid(2)%vector(nprev+1:), &
01944 tgt_grid(3)%vector(nprev+1:), &
01945 mp%coords_pointer%coords_real(1)%vector, &
01946 mp%coords_pointer%coords_real(2)%vector, &
01947 mp%coords_pointer%coords_real(3)%vector, &
01948 coords_shape, &
01949 mask_array, mask_shape, src_mask_available, &
01950 sinvec(lon)%vector, cosvec(lon)%vector, &
01951 sinvec(lat)%vector, cosvec(lat)%vector, &
01952 grid_shape, search_mode, &
01953 send_info%srclocs(1,1)%vector, &
01954 len_cpl(ipart), ncpl, nprev, &
01955 neighbors_3d, nb_extra, &
01956 extra_search, ierror)
01957
01958 if (ierror > 0) return
01959
01960 nprev = nprev + len_cpl (ipart)
01961 end if
01962 end do
01963 else if (send_info%nvec == 2) then
01964
01965
01966 do ipart = 1, search%npart
01967 if (len_cpl (ipart) > 0) then
01968 call psmile_neigh_nearx_3d_irr2_real (grid_id, &
01969 tgt_grid(1)%vector(nprev+1:), &
01970 tgt_grid(2)%vector(nprev+1:), &
01971 tgt_grid(3)%vector(nprev+1:), &
01972 mp%coords_pointer%coords_real(1)%vector, &
01973 mp%coords_pointer%coords_real(2)%vector, &
01974 mp%coords_pointer%coords_real(3)%vector, &
01975 coords_shape, &
01976 mask_array, mask_shape, src_mask_available, &
01977 sinvec(lon)%vector, cosvec(lon)%vector, &
01978 sinvec(lat)%vector, cosvec(lat)%vector, &
01979 grid_shape, search_mode, &
01980 send_info%srclocs(indl, ipart)%vector, &
01981 send_info%srclocs(indz, ipart)%vector, &
01982 send_info%npoints(1:indz, ipart), &
01983 ncpl, nprev, &
01984 neighbors_3d, nb_extra, &
01985 extra_search, ierror)
01986
01987 if (ierror > 0) return
01988
01989 nprev = nprev + len_cpl (ipart)
01990 end if
01991 end do
01992
01993 else
01994 do ipart = 1, search%npart
01995 if (len_cpl (ipart) > 0) then
01996 call psmile_neigh_nearx_3d_reg_real (grid_id, &
01997 tgt_grid(1)%vector(nprev+1:), &
01998 tgt_grid(2)%vector(nprev+1:), &
01999 tgt_grid(3)%vector(nprev+1:), &
02000 mp%coords_pointer%coords_real(1)%vector, &
02001 mp%coords_pointer%coords_real(2)%vector, &
02002 mp%coords_pointer%coords_real(3)%vector, &
02003 coords_shape, &
02004 mask_array, mask_shape, src_mask_available, &
02005 sinvec(lon)%vector, cosvec(lon)%vector, &
02006 sinvec(lat)%vector, cosvec(lat)%vector, &
02007 grid_shape, search_mode, &
02008 send_info%srclocs(1:ndim_3d, ipart), &
02009 send_info%npoints(1:ndim_3d, ipart), &
02010 ncpl, nprev, &
02011 neighbors_3d, nb_extra, &
02012 extra_search, ierror)
02013
02014 if (ierror > 0) return
02015
02016 nprev = nprev + len_cpl (ipart)
02017 end if
02018 end do
02019 endif
02020
02021
02022
02023
02024 if (global_search) then
02025 call psmile_global_search_nnx_real (comp_info, &
02026 search, var_id, &
02027 tgt_grid(1)%vector, &
02028 tgt_grid(2)%vector, &
02029 tgt_grid(3)%vector, &
02030 neighbors_3d, ncpl, nb_neighbors, &
02031 nb_extra, extra_search, &
02032 send_index, ierror)
02033 if (ierror > 0) return
02034 endif
02035
02036
02037
02038 call psmile_neigh_extra_search_clean (search, extra_search, ierror)
02039 endif
02040
02041
02042
02043
02044 if (transformed) then
02045 do i = 1, ndim_2d
02046 Deallocate (cosvec(i)%vector)
02047 Deallocate (sinvec(i)%vector)
02048 end do
02049 endif
02050
02051
02052
02053
02054
02055 if ( use_source_cell ) then
02056
02057 nbr_cells_tot = sum(nbr_cells)
02058
02059
02060
02061 nb_corners = 2
02062
02063 Allocate (neighcells(nbr_cells_tot, nb_corners), &
02064 source_cell_index(nbr_cells_tot), STAT = ierror)
02065 if ( ierror > 0 ) then
02066 ierrp (1) = ierror
02067 ierrp (2) = nbr_cells_tot * ( nb_corners + 1 )
02068 ierror = PRISM_Error_Alloc
02069 call psmile_error ( ierror, 'neighcells & source_cell_index', &
02070 ierrp, 2, __FILE__, __LINE__ )
02071 return
02072 endif
02073
02074 corner_shape = grid_valid_shape
02075
02076
02077
02078 nb_corners = 1
02079
02080 call psmile_compact_neighbors_3d (neighcells_3d, nbr_cells_tot, &
02081 nb_corners, corner_shape, &
02082 extra_search, send_info, &
02083 neighcells, ierror)
02084 #ifdef MOVE0_FOR_CELLS
02085
02086
02087
02088
02089
02090
02091 call psmile_move0_neighbors (neighcells, nbr_cells_tot, nb_neighbors, &
02092 ierror)
02093 if (ierror > 0) return
02094 #endif
02095
02096
02097
02098 src_cell_size = send_info%n_list
02099
02100 do i = 1, ndim_3d
02101 Allocate (src_cell(i)%vector(2*src_cell_size), stat = ierror)
02102 if ( ierror > 0 ) then
02103 ierrp (1) = ierror
02104 ierrp (2) = 2*src_cell_size
02105
02106 ierror = PRISM_Error_Alloc
02107 call psmile_error ( ierror, 'src_cell(i)%vector', &
02108 ierrp, 2, __FILE__, __LINE__ )
02109 return
02110 endif
02111 end do
02112
02113
02114
02115
02116 corner_shape(1:2,1:2) = cp%corner_shape(1:2,1:2)
02117 corner_shape(:,3) = mp%coords_pointer%coords_shape(:,3)
02118
02119 call psmile_ccompact_gauss2_real ( send_info, &
02120 grid_valid_shape, &
02121 corner_shape, cp%nbr_corners/4, &
02122 cp%corners_real(1)%vector, &
02123 cp%corners_real(2)%vector, &
02124 mp%coords_pointer%coords_real(3)%vector, &
02125 extra_search, src_cell_size, &
02126 nbr_cells_tot, source_cell_index, &
02127 neighcells, &
02128 src_cell(1)%vector, &
02129 src_cell(2)%vector, &
02130 src_cell(3)%vector, ierror )
02131 if (ierror /= 0) return
02132
02133 if (associated(neighcells_3d)) Deallocate (neighcells_3d)
02134
02135 else
02136
02137
02138
02139 Allocate (neighbors(ncpl, nb_neighbors), STAT = ierror)
02140 if ( ierror > 0 ) then
02141 ierrp (1) = ierror
02142 ierrp (2) = ncpl * nb_neighbors
02143
02144 ierror = PRISM_Error_Alloc
02145 call psmile_error ( ierror, 'neighbors', &
02146 ierrp, 2, __FILE__, __LINE__ )
02147 return
02148 endif
02149
02150 #ifdef DEBUG_EXTRA
02151 n = 4919
02152 do i = 1, 16
02153 print *, 'source coords', i, neighbors_3d(1,n,i), &
02154 Methods(1)%coords_pointer%coords_real(1)%vector(neighbors_3d(1,n,i)), &
02155 Methods(1)%coords_pointer%coords_real(2)%vector(neighbors_3d(1,n,i))
02156 enddo
02157 #endif
02158 #ifdef DEBUG_TRACE
02159 if (ictl <= ncpl) then
02160 print *, 'before compact: i, dstijk (:, i)', ictl, dstijk (:, ictl)
02161 do i = 1, nb_neighbors
02162 print *, 'neighbors_3d (:,ictl)', neighbors_3d (:, ictl, i)
02163 end do
02164 endif
02165 #endif
02166 call psmile_compact_neighbors_3d ( &
02167 neighbors_3d, ncpl, nb_neighbors, &
02168 grid_valid_shape, extra_search, &
02169 send_info, neighbors, ierror)
02170 if (ierror > 0) return
02171
02172
02173
02174 #ifdef DEBUG_TRACE
02175 if (ictl <= ncpl) then
02176 print *, 'before move0: i, dstijk (:, i)', ictl, dstijk (:, ictl)
02177 do i = 1, nb_neighbors
02178 print *, 'neighbors_3d (:,ictl)', neighbors_3d (:, ictl, i)
02179 end do
02180 print *, 'neighbors (ictl)', neighbors (ictl, 1:nb_neighbors)
02181 endif
02182 #endif
02183 if (interpolation_methods(1) /= PSMILe_nnghbr3D) then
02184 call psmile_move0_neighbors (neighbors, ncpl, nb_neighbors, &
02185 ierror)
02186 if (ierror > 0) return
02187 endif
02188
02189 #ifdef DEBUG_TRACE
02190 if (ictl <= ncpl) then
02191 print *, 'after move0: i, dstijk (:, i)', ictl, dstijk (:, ictl)
02192 do i = 1, nb_neighbors
02193 print *, 'neighbors_3d (:, ictl)', neighbors_3d (:, ictl, i)
02194 end do
02195 print *, 'neighbors (ictl)', neighbors (ictl, 1:nb_neighbors)
02196 if (allocated (ictl_neighbours)) deallocate (ictl_neighbours)
02197 allocate (ictl_neighbours(nb_neighbors))
02198 ictl_neighbours(:) = neighbors (ictl, 1:nb_neighbors)
02199 endif
02200 #endif
02201 Deallocate (neighbors_3d)
02202
02203 endif
02204
02205
02206
02207
02208
02209
02210
02211
02212
02213
02214
02215
02216
02217
02218
02219
02220
02221 if ( use_target_cell ) then
02222 if ( send_info%nvec == 1 ) then
02223 send_info%npoints(1,1) = sum(len_cpl(1:search%npart))
02224 else if ( send_info%nvec == 2 ) then
02225 do ipart = 1, search%npart
02226 if ( send_info%npoints(indz,ipart) /= 0 ) then
02227 send_info%npoints(indl,ipart) = len_cpl(ipart) / &
02228 send_info%npoints(indz,ipart)
02229 else
02230 send_info%npoints(indl,ipart) = 0
02231 endif
02232 enddo
02233 endif
02234 endif
02235
02236 src_size = send_info%n_list
02237 local_src_size = src_size - send_info%num2recv
02238
02239
02240
02241 if ( src_size == 0 ) then
02242 send_index = PRISM_UNDEFINED
02243
02244 go to 1000
02245 endif
02246
02247
02248
02249 do j = 1, n_intmethods
02250 if ( interpolation_methods(j) == PSMILe_conserv2D .or. &
02251 interpolation_methods(j) == PSMILe_conserv3D ) then
02252 if ( src_mask_available ) then
02253 id_src_mask = 1
02254 use_mask = .false.
02255 endif
02256 endif
02257 enddo
02258
02259
02260
02261 if (src_mask_available .and. .not. use_mask) then
02262 id_src_mask = 1
02263
02264 Allocate (src_mask(src_size), stat = ierror)
02265 if ( ierror > 0 ) then
02266 ierrp (1) = ierror
02267 ierrp (2) = src_size
02268
02269 ierror = PRISM_Error_Alloc
02270 call psmile_error ( ierror, 'src_mask', &
02271 ierrp, 2, __FILE__, __LINE__ )
02272 return
02273 endif
02274
02275
02276
02277
02278 call psmile_ext_compact_list_log2int (send_info, &
02279 extra_search, mask_array, mask_shape, &
02280 grid_valid_shape, src_mask, src_size, ierror)
02281 if (ierror /= 0) return
02282
02283 else
02284
02285 id_src_mask = 0
02286 src_mask => dummy_mask
02287
02288 endif
02289
02290
02291
02292 if ( .not. use_source_cell ) then
02293
02294 do i = 1, ndim_3d
02295 Allocate (src_grid(i)%vector(src_size), stat = ierror)
02296 if ( ierror > 0 ) then
02297 ierrp (1) = ierror
02298 ierrp (2) = src_size
02299
02300 ierror = PRISM_Error_Alloc
02301 call psmile_error ( ierror, 'src_grid(i)%vector', &
02302 ierrp, 2, __FILE__, __LINE__ )
02303 return
02304 endif
02305 end do
02306
02307
02308
02309 if (send_info%send_entire_valid_shape) then
02310 len = grid_valid_shape(2,1) - grid_valid_shape(1,1) + 1
02311 shift = grid_valid_shape (1,1) - mp%coords_pointer%coords_shape(1,1)
02312
02313 do i = 1, 2
02314 do n = 1, local_src_size / len
02315 src_grid(i)%vector((n-1)*len+1:n*len) = &
02316 mp%coords_pointer%coords_real(i)%vector(shift+1:shift+len)
02317 end do
02318 end do
02319
02320 i = 3
02321 shift = grid_valid_shape (1,i) - mp%coords_pointer%coords_shape(1,i)
02322
02323 do n = 1, local_src_size / len
02324 src_grid(i)%vector((n-1)*len+1:n*len) = &
02325 mp%coords_pointer%coords_real(i)%vector(shift+n)
02326 end do
02327
02328 else
02329
02330 list_entries => send_info%list_entries
02331
02332 #ifdef PRISM_ASSERTION
02333 if (.not. Associated (list_entries)) then
02334 call psmile_assert (__FILE__, __LINE__, &
02335 "list_entries are not available")
02336 endif
02337 #endif
02338
02339 shift = mp%coords_pointer%coords_shape(1,1) - 1
02340 do i = 1, 2
02341 src_grid(i)%vector(1:local_src_size) = &
02342 mp%coords_pointer%coords_real(i)%vector( &
02343 list_entries(1, 1:local_src_size) - shift)
02344 end do
02345
02346 i = 3
02347 shift = mp%coords_pointer%coords_shape(1,i) - 1
02348
02349 src_grid(i)%vector(1:local_src_size) = &
02350 mp%coords_pointer%coords_dble(i)%vector( &
02351 list_entries(i, 1:local_src_size) - shift)
02352 endif
02353
02354
02355
02356 if (send_info%num2recv > 0) then
02357 do i = 1, ndim_3d
02358 ip = local_src_size
02359
02360 do n = 1, send_info%nrecv
02361 src_grid(i)%vector(ip+1:ip+send_info%len_sent(n)) = &
02362 extra_search%real_bufs(n)%vector( &
02363 (i-1)*send_info%len_sent(n)+1: &
02364 i *send_info%len_sent(n))
02365
02366 ip = ip + send_info%len_sent (n)
02367 end do
02368 end do
02369
02370 do n = 1, send_info%nrecv
02371 Deallocate (extra_search%real_bufs(n)%vector)
02372 end do
02373 endif
02374
02375 #ifdef DEBUG_EXTRA
02376 n = 4919
02377 do i = 1, 16
02378 print *, 'src grid', i, neighbors(n,i), &
02379 src_grid(1)%vector(neighbors(n,i)), src_grid(2)%vector(neighbors(n,i))
02380 enddo
02381 #endif
02382
02383 #ifdef DEBUG_TRACE
02384 if (ictl <= ncpl) then
02385
02386 ipart = 1
02387
02388 if (search%grid_type == PRISM_Reglonlatvrt) then
02389 call psmile_print_3d_reg_coord_real (coords(1,ipart)%vector, &
02390 coords(2,ipart)%vector, &
02391 coords(3,ipart)%vector, &
02392 shape (:, :, ipart), &
02393 dstijk (:, ictl), 1, 'searched')
02394 else if (search%grid_type == PRISM_Irrlonlatvrt) then
02395 call psmile_print_3d_coord_real (coords(1,ipart)%vector, &
02396 coords(2,ipart)%vector, &
02397 coords(3,ipart)%vector, &
02398 shape (:, :, ipart), &
02399 dstijk (:, ictl), 1, 'searched')
02400 endif
02401
02402 print *, 'coords (ictl)'
02403 do i = 1, nb_neighbors
02404 if (neighbors (ictl, i) > 0) then
02405 print *, 'i', i, &
02406 src_grid(1)%vector(neighbors (ictl, i)),&
02407 src_grid(2)%vector(neighbors (ictl, i)),&
02408 src_grid(3)%vector(neighbors (ictl, i))
02409 endif
02410 end do
02411 endif
02412 #endif
02413 endif
02414
02415 call psmile_get_epio_handle (field%comp_id, grid_id, method_id, mask_id, &
02416 interpolation_methods, search%msg_intersections, id_trans_out, &
02417 id_trans_in, search%sender, cpl_id, epio_id, trs_rank, ierror)
02418
02419 if ( epio_id /= PSMILe_undef ) then
02420
02421
02422
02423
02424
02425
02426 call psmile_trs_set_triple_links (id_trans_out, id_trans_in, &
02427 epio_id, trs_rank, ierror)
02428 if (ierror /= 0) return
02429
02430
02431
02432 if ( use_source_cell ) then
02433 do i = 1, ndim_3d
02434 Deallocate (src_cell(i)%vector)
02435 end do
02436 else
02437 do i = 1, ndim_3d
02438 Deallocate (src_grid(i)%vector)
02439 end do
02440 endif
02441
02442 if (id_src_mask > 0) then
02443 Deallocate (src_mask)
02444 endif
02445
02446 if ( use_target_cell ) then
02447 do i = 1, ndim_3d
02448 Deallocate (tgt_cell(i)%vector)
02449 end do
02450 else
02451 if (use_target_grid) then
02452 do i = 1, ndim_3d
02453 Deallocate (tgt_grid(i)%vector)
02454 end do
02455 endif
02456 endif
02457
02458 else
02459
02460
02461
02462
02463
02464
02465
02466
02467 if ( use_source_cell ) then
02468
02469 call psmile_trs_set_src_epio3d_real (epio_id, &
02470 trs_rank, src_size, 2*src_cell_size, &
02471 src_cell(1)%vector, &
02472 src_cell(2)%vector, &
02473 src_cell(3)%vector, &
02474 id_src_mask, src_mask, ierror)
02475
02476
02477
02478 do i = 1, ndim_3d
02479 Deallocate (src_cell(i)%vector)
02480 end do
02481 else
02482
02483 call psmile_trs_set_src_epio3d_real (epio_id, &
02484 trs_rank, src_size, src_size, &
02485 src_grid(1)%vector, &
02486 src_grid(2)%vector, &
02487 src_grid(3)%vector, &
02488 id_src_mask, src_mask, ierror)
02489
02490
02491
02492 #ifdef DEBUG_TRACE
02493 if (ictl <= ncpl) then
02494 print *, "src_grid(1:3)%vector(neighbors(ictl,1:", nb_neighbors, ")):"
02495 do i = 1, nb_neighbors
02496 if (neighbors (ictl, i) > 0) then
02497 print *, i, "(", neighbors(ictl,i), ") :", &
02498 src_grid(1)%vector(neighbors(ictl,i)), &
02499 src_grid(2)%vector(neighbors(ictl,i)), &
02500 src_grid(3)%vector(neighbors(ictl,i))
02501 endif
02502 enddo
02503 endif
02504 #endif
02505 do i = 1, ndim_3d
02506 Deallocate (src_grid(i)%vector)
02507 end do
02508 endif
02509
02510 if (ierror /= 0) return
02511
02512
02513
02514 if (id_src_mask == 1) then
02515 Deallocate (src_mask)
02516 endif
02517
02518
02519
02520
02521
02522
02523
02524
02525 id_tgt_mask = 0
02526
02527
02528 do j = 1, n_intmethods
02529 if ( (interpolation_methods(j) == PSMILe_conserv2D .or. &
02530 interpolation_methods(j) == PSMILe_conserv3D ) .and. &
02531 associated(search%search_mask) ) id_tgt_mask = 0
02532 enddo
02533
02534 if ( id_tgt_mask == 1 ) then
02535
02536 Allocate (tgt_mask(ncpl), stat = ierror)
02537 if ( ierror > 0 ) then
02538 ierrp (1) = ierror
02539 ierrp (2) = ncpl
02540
02541 ierror = PRISM_Error_Alloc
02542 call psmile_error ( ierror, 'tgt_mask', &
02543 ierrp, 2, __FILE__, __LINE__ )
02544 return
02545 endif
02546
02547 tgt_mask = 0
02548 nprev = 0
02549
02550 do ipart = 1, search%npart
02551 if (len_cpl (ipart) > 0) then
02552 ibeg = nprev + 1
02553 iend = nprev + len_cpl (ipart)
02554
02555 leni = shape(2,1,ipart)-shape(1,1,ipart)
02556 lenij = leni * (shape(2,2,ipart)-shape(1,2,ipart))
02557
02558 do j = ibeg, iend
02559
02560 if ( dstijk(1,j) /= PSMILe_Undef .and. &
02561 dstijk(2,j) /= PSMILe_Undef .and. &
02562 dstijk(3,j) /= PSMILe_Undef ) then
02563
02564 n = (dstijk(3,j)-shape(1,3,ipart))*lenij + &
02565 (dstijk(2,j)-shape(1,2,ipart))*leni + &
02566 dstijk(1,j)-shape(1,1,ipart)+1
02567
02568 if ( search%search_mask(ipart)%vector(n) ) tgt_mask(j) = 1
02569
02570 endif
02571
02572 enddo
02573
02574 nprev = nprev + len_cpl (ipart)
02575 endif
02576 end do
02577
02578 else
02579
02580 id_tgt_mask = 0
02581 tgt_mask => dummy_mask
02582
02583 endif
02584
02585
02586
02587
02588
02589
02590
02591
02592
02593
02594 if ( use_target_cell ) then
02595
02596 nb_corners = 4
02597 call psmile_trs_set_tgt_epio3d_real ( &
02598 epio_id, trs_rank, ncpl, nb_corners, &
02599 tgt_cell(1)%vector, tgt_cell(2)%vector, tgt_cell(3)%vector, &
02600 id_tgt_mask, tgt_mask, ierror)
02601
02602
02603
02604 do i = 1, ndim_3d
02605 Deallocate (tgt_cell(i)%vector)
02606 end do
02607 else
02608 #ifdef DEBUG_TRACE
02609 if (ictl <= ncpl) then
02610 print *, 'ictl target coord', ictl, tgt_grid(1)%vector(ictl), &
02611 tgt_grid(2)%vector(ictl), &
02612 tgt_grid(3)%vector(ictl)
02613 endif
02614 #endif
02615 nb_corners = 1
02616
02617 call psmile_trs_set_tgt_epio3d_real ( &
02618 epio_id, trs_rank, ncpl, nb_corners, &
02619 tgt_grid(1)%vector, tgt_grid(2)%vector, tgt_grid(3)%vector, &
02620 id_tgt_mask, dummy_mask, ierror)
02621 if (ierror /= 0) return
02622
02623
02624
02625 if (use_target_grid) then
02626 do i = 1, ndim_3d
02627 Deallocate (tgt_grid(i)%vector)
02628 end do
02629 endif
02630 endif
02631
02632
02633
02634 if ( id_tgt_mask == 1 ) Deallocate(tgt_mask)
02635
02636
02637
02638
02639
02640 call psmile_trs_set_triple_links (id_trans_out, id_trans_in, &
02641 epio_id, trs_rank, ierror)
02642 if (ierror /= 0) return
02643
02644
02645
02646 if ( use_source_cell ) then
02647 nb_corners = 2
02648 call psmile_trs_give_neighcells3d (epio_id, trs_rank, &
02649 ncpl, nbr_cells, nbr_cells_tot, nb_corners, &
02650 source_cell_index, neighcells, grid_info%grid_type, ierror)
02651 else
02652 call psmile_trs_give_neighbors_gauss (epio_id, trs_rank, &
02653 ncpl, nb_neighbors, neighbors, neigh_bascule, ierror)
02654 endif
02655 if (ierror /= 0) return
02656
02657
02658
02659 cpl_list(cpl_id)%epio_id = epio_id
02660 cpl_list(cpl_id)%trs_rank = trs_rank
02661
02662 endif
02663
02664 send_info%epio_id = epio_id
02665 send_info%trs_rank = trs_rank
02666
02667 #ifdef DEBUG
02668 print *, trim(ch_id), 'index, epio_id, n_list ', &
02669 send_index, epio_id, send_info%n_list, &
02670 send_info%send_entire_valid_shape
02671 #endif
02672
02673
02674
02675 1000 continue
02676
02677 call psmile_locations_dealloc (send_info, ierror)
02678
02679 if ( use_source_cell ) then
02680 Deallocate (neighcells, source_cell_index, nbr_cells)
02681 else
02682 Deallocate (neighbors)
02683 Deallocate (neigh_bascule)
02684 endif
02685
02686
02687
02688 #ifdef VERBOSE
02689 print 9980, trim(ch_id), ierror
02690
02691 call psmile_flushstd
02692 #endif /* VERBOSE */
02693
02694
02695
02696
02697 #ifdef VERBOSE
02698
02699 9990 format (1x, a, ': psmile_info_trs_loc_gauss2_real: method_id', i3, &
02700 ' to ', i3, '(', i2, ')')
02701 9980 format (1x, a, ': psmile_info_trs_loc_gauss2_real: eof ierror =', i3)
02702
02703 9970 format (1x, a, ': psmile_info_trs_loc_gauss2_real: eof src_size=0, ierror =', i3)
02704
02705 #endif /* VERBOSE */
02706
02707 #ifdef DEBUG
02708 9960 format (1x, 'Taskout%In_channel(', i3, '): origin_type', i5, &
02709 ', remote comp id', i4, &
02710 ', global_transi_id', i4, ', remote_transi_id', i4)
02711 #endif
02712
02713 end subroutine PSMILe_info_trs_loc_gauss2_real
