1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980818283848586878889909192939495969798991001011021031041051061071081091101111121131141151161171181191201211221231241251261271281291301311321331341351361371381391401411421431441451461471481491501511521531541551561571581591601611621631641651661671681691701711721731741751761771781791801811821831841851861871881891901911921931941951961971981992002012022032042052062072082092102112122132142152162172182192202212222232242252262272282292302312322332342352362372382392402412422432442452462472482492502512522532542552562572582592602612622632642652662672682692702712722732742752762772782792802812822832842852862872882892902912922932942952962972982993003013023033043053063073083093103113123133143153163173183193203213223233243253263273283293303313323333343353363373383393403413423433443453463473483493503513523533543553563573583593603613623633643653663673683693703713723733743753763773783793803813823833843853863873883893903913923933943953963973983994004014024034044054064074084094104114124134144154164174184194204214224234244254264274284294304314324334344354364374384394404414424434444454464474484494504514524534544554564574584594604614624634644654664674684694704714724734744754764774784794804814824834844854864874884894904914924934944954964974984995005015025035045055065075085095105115125135145155165175185195205215225235245255265275285295305315325335345355365375385395405415425435445455465475485495505515525535545555565575585595605615625635645655665675685695705715725735745755765775785795805815825835845855865875885895905915925935945955965975985996006016026036046056066076086096106116126136146156166176186196206216226236246256266276286296306316326336346356366376386396406416426436446456466476486496506516526536546556566576586596606616626636646656666676686696706716726736746756766776786796806816826836846856866876886896906916926936946956966976986997007017027037047057067077087097107117127137147157167177187197207217227237247257267277287297307317327337347357367377387397407417427437447457467477487497507517527537547557567577587597607617627637647657667677687697707717727737747757767777787797807817827837847857867877887897907917927937947957967977987998008018028038048058068078088098108118128138148158168178188198208218228238248258268278288298308318328338348358368378388398408418428438448458468478488498508518528538548558568578588598608618628638648658668678688698708718728738748758768778788798808818828838848858868878888898908918928938948958968978988999009019029039049059069079089099109119129139149159169179189199209219229239249259269279289299309319329339349359369379389399409419429439449459469479489499509519529539549559569579589599609619629639649659669679689699709719729739749759769779789799809819829839849859869879889899909919929939949959969979989991000100110021003100410051006100710081009101010111012101310141015101610171018101910201021102210231024102510261027102810291030103110321033103410351036103710381039104010411042104310441045104610471048104910501051105210531054105510561057105810591060106110621063106410651066106710681069107010711072107310741075107610771078107910801081108210831084108510861087108810891090109110921093109410951096109710981099110011011102110311041105110611071108110911101111111211131114111511161117111811191120112111221123112411251126112711281129113011311132113311341135113611371138113911401141114211431144114511461147114811491150115111521153115411551156115711581159116011611162116311641165116611671168116911701171117211731174117511761177117811791180118111821183118411851186118711881189119011911192119311941195119611971198119912001201120212031204120512061207120812091210121112121213121412151216121712181219122012211222122312241225122612271228122912301231123212331234123512361237123812391240124112421243124412451246124712481249125012511252125312541255125612571258125912601261126212631264126512661267126812691270127112721273127412751276127712781279128012811282128312841285128612871288128912901291129212931294129512961297129812991300130113021303130413051306130713081309131013111312131313141315131613171318131913201321132213231324132513261327132813291330133113321333133413351336133713381339134013411342134313441345134613471348134913501351135213531354135513561357135813591360136113621363136413651366136713681369137013711372137313741375137613771378137913801381138213831384138513861387138813891390139113921393139413951396139713981399140014011402140314041405140614071408140914101411141214131414141514161417141814191420142114221423142414251426142714281429143014311432143314341435143614371438143914401441144214431444144514461447144814491450145114521453145414551456145714581459146014611462146314641465146614671468146914701471147214731474147514761477147814791480148114821483148414851486148714881489149014911492149314941495149614971498149915001501150215031504150515061507150815091510151115121513151415151516151715181519152015211522152315241525152615271528152915301531153215331534153515361537153815391540154115421543154415451546154715481549155015511552155315541555155615571558155915601561156215631564156515661567156815691570157115721573157415751576157715781579158015811582158315841585158615871588158915901591159215931594159515961597159815991600160116021603160416051606160716081609161016111612161316141615161616171618161916201621162216231624162516261627162816291630163116321633163416351636163716381639164016411642164316441645164616471648164916501651165216531654165516561657165816591660166116621663166416651666166716681669167016711672167316741675167616771678167916801681168216831684168516861687168816891690169116921693169416951696169716981699170017011702170317041705170617071708170917101711171217131714171517161717171817191720172117221723172417251726172717281729173017311732173317341735173617371738173917401741174217431744174517461747174817491750175117521753175417551756175717581759176017611762176317641765176617671768176917701771177217731774177517761777177817791780178117821783178417851786178717881789179017911792179317941795179617971798179918001801180218031804180518061807180818091810181118121813181418151816181718181819182018211822182318241825182618271828182918301831183218331834183518361837183818391840184118421843184418451846184718481849185018511852185318541855185618571858185918601861186218631864186518661867186818691870187118721873187418751876187718781879188018811882188318841885188618871888188918901891189218931894189518961897189818991900190119021903190419051906190719081909191019111912191319141915191619171918191919201921192219231924192519261927192819291930193119321933193419351936193719381939194019411942194319441945194619471948194919501951195219531954195519561957195819591960196119621963196419651966196719681969197019711972197319741975197619771978197919801981198219831984198519861987198819891990199119921993199419951996199719981999200020012002200320042005200620072008200920102011201220132014201520162017201820192020202120222023202420252026202720282029203020312032203320342035203620372038203920402041204220432044204520462047204820492050205120522053205420552056205720582059206020612062206320642065206620672068206920702071207220732074207520762077207820792080208120822083208420852086208720882089209020912092209320942095209620972098209921002101210221032104210521062107210821092110211121122113211421152116211721182119212021212122212321242125212621272128212921302131213221332134213521362137213821392140214121422143214421452146214721482149215021512152215321542155215621572158215921602161216221632164216521662167216821692170217121722173217421752176217721782179218021812182218321842185218621872188218921902191219221932194219521962197 |
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- !
- ! this module contains necessary routines for computing addresses
- ! and weights for a conservative interpolation between any two
- ! grids on a sphere. the weights are computed by performing line
- ! integrals around all overlap regions of the two grids. see
- ! Dukowicz and Kodis, SIAM J. Sci. Stat. Comput. 8, 305 (1987) and
- ! Jones, P.W. Monthly Weather Review (submitted).
- !
- !-----------------------------------------------------------------------
- !
- ! CVS:$Id: remap_conserv.f,v 1.10 2001/08/21 21:05:13 pwjones Exp $
- !
- ! Copyright (c) 1997, 1998 the Regents of the University of
- ! California.
- !
- ! This software and ancillary information (herein called software)
- ! called SCRIP is made available under the terms described here.
- ! The software has been approved for release with associated
- ! LA-CC Number 98-45.
- !
- ! Unless otherwise indicated, this software has been authored
- ! by an employee or employees of the University of California,
- ! operator of the Los Alamos National Laboratory under Contract
- ! No. W-7405-ENG-36 with the U.S. Department of Energy. The U.S.
- ! Government has rights to use, reproduce, and distribute this
- ! software. The public may copy and use this software without
- ! charge, provided that this Notice and any statement of authorship
- ! are reproduced on all copies. Neither the Government nor the
- ! University makes any warranty, express or implied, or assumes
- ! any liability or responsibility for the use of this software.
- !
- ! If software is modified to produce derivative works, such modified
- ! software should be clearly marked, so as not to confuse it with
- ! the version available from Los Alamos National Laboratory.
- !
- !***********************************************************************
- module remap_conservative
- !-----------------------------------------------------------------------
- use kinds_mod ! defines common data types
- use constants ! defines common constants
- use timers ! module for timing
- use grids ! module containing grid information
- use remap_vars ! module containing remap information
- implicit none
- !-----------------------------------------------------------------------
- !
- ! module variables
- !
- !-----------------------------------------------------------------------
- integer (kind=int_kind), save ::
- & num_srch_cells ! num cells in restricted search arrays
- integer (kind=int_kind), dimension(:), allocatable, save ::
- & srch_add ! global address of cells in srch arrays
- real (kind=dbl_kind), parameter ::
- & north_thresh = 1.45_dbl_kind, ! threshold for coord transf.
- & south_thresh =-2.00_dbl_kind ! threshold for coord transf.
- real (kind=dbl_kind), dimension(:,:), allocatable, save ::
- & srch_corner_lat, ! lat of each corner of srch cells
- & srch_corner_lon ! lon of each corner of srch cells
- !***********************************************************************
- contains
- !***********************************************************************
- subroutine remap_conserv
- !-----------------------------------------------------------------------
- !
- ! this routine traces the perimeters of every grid cell on each
- ! grid checking for intersections with the other grid and computing
- ! line integrals for each subsegment.
- !
- !-----------------------------------------------------------------------
- !-----------------------------------------------------------------------
- !
- ! local variables
- !
- !-----------------------------------------------------------------------
- integer (kind=int_kind), parameter ::
- & max_subseg = 10000 ! max number of subsegments per segment
- ! to prevent infinite loop
- integer (kind=int_kind) ::
- & grid1_add, ! current linear address for grid1 cell
- & grid2_add, ! current linear address for grid2 cell
- & min_add, ! addresses for restricting search of
- & max_add, ! destination grid
- & n, nwgt, ! generic counters
- & corner, ! corner of cell that segment starts from
- & next_corn, ! corner of cell that segment ends on
- & num_subseg ! number of subsegments
- logical (kind=log_kind) ::
- & lcoinc, ! flag for coincident segments
- & lrevers, ! flag for reversing direction of segment
- & lbegin ! flag for first integration of a segment
- logical (kind=log_kind), dimension(:), allocatable ::
- & srch_mask ! mask for restricting searches
- real (kind=dbl_kind) ::
- & intrsct_lat, intrsct_lon, ! lat/lon of next intersect
- & beglat, endlat, beglon, endlon, ! endpoints of current seg.
- & norm_factor ! factor for normalizing wts
- real (kind=dbl_kind), dimension(:), allocatable ::
- & grid2_centroid_lat, grid2_centroid_lon, ! centroid coords
- & grid1_centroid_lat, grid1_centroid_lon ! on each grid
- real (kind=dbl_kind), dimension(2) :: begseg ! begin lat/lon for
- ! full segment
- real (kind=dbl_kind), dimension(6) :: weights ! local wgt array
- !-----------------------------------------------------------------------
- !
- ! initialize centroid arrays
- !
- !-----------------------------------------------------------------------
- allocate( grid1_centroid_lat(grid1_size),
- & grid1_centroid_lon(grid1_size),
- & grid2_centroid_lat(grid2_size),
- & grid2_centroid_lon(grid2_size))
- grid1_centroid_lat = zero
- grid1_centroid_lon = zero
- grid2_centroid_lat = zero
- grid2_centroid_lon = zero
- !-----------------------------------------------------------------------
- !
- ! integrate around each cell on grid1
- !
- !-----------------------------------------------------------------------
- allocate(srch_mask(grid2_size))
- print *,'grid1 sweep '
- do grid1_add = 1,grid1_size
- !***
- !*** restrict searches first using search bins
- !***
- call timer_start(1)
- min_add = grid2_size
- max_add = 1
- do n=1,num_srch_bins
- if (grid1_add >= bin_addr1(1,n) .and.
- & grid1_add <= bin_addr1(2,n)) then
- min_add = min(min_add, bin_addr2(1,n))
- max_add = max(max_add, bin_addr2(2,n))
- endif
- end do
- !***
- !*** further restrict searches using bounding boxes
- !***
- num_srch_cells = 0
- do grid2_add = min_add,max_add
- srch_mask(grid2_add) = (grid2_bound_box(1,grid2_add) <=
- & grid1_bound_box(2,grid1_add)) .and.
- & (grid2_bound_box(2,grid2_add) >=
- & grid1_bound_box(1,grid1_add)) .and.
- & (grid2_bound_box(3,grid2_add) <=
- & grid1_bound_box(4,grid1_add)) .and.
- & (grid2_bound_box(4,grid2_add) >=
- & grid1_bound_box(3,grid1_add))
- if (srch_mask(grid2_add)) num_srch_cells = num_srch_cells+1
- end do
- !***
- !*** create search arrays
- !***
- allocate(srch_add(num_srch_cells),
- & srch_corner_lat(grid2_corners,num_srch_cells),
- & srch_corner_lon(grid2_corners,num_srch_cells))
- n = 0
- gather1: do grid2_add = min_add,max_add
- if (srch_mask(grid2_add)) then
- n = n+1
- srch_add(n) = grid2_add
- srch_corner_lat(:,n) = grid2_corner_lat(:,grid2_add)
- srch_corner_lon(:,n) = grid2_corner_lon(:,grid2_add)
- endif
- end do gather1
- call timer_stop(1)
- !***
- !*** integrate around this cell
- !***
- do corner = 1,grid1_corners
- next_corn = mod(corner,grid1_corners) + 1
- !***
- !*** define endpoints of the current segment
- !***
- beglat = grid1_corner_lat(corner,grid1_add)
- beglon = grid1_corner_lon(corner,grid1_add)
- endlat = grid1_corner_lat(next_corn,grid1_add)
- endlon = grid1_corner_lon(next_corn,grid1_add)
- lrevers = .false.
- !***
- !*** to ensure exact path taken during both
- !*** sweeps, always integrate segments in the same
- !*** direction (SW to NE).
- !***
- if ((endlat < beglat) .or.
- & (endlat == beglat .and. endlon < beglon)) then
- beglat = grid1_corner_lat(next_corn,grid1_add)
- beglon = grid1_corner_lon(next_corn,grid1_add)
- endlat = grid1_corner_lat(corner,grid1_add)
- endlon = grid1_corner_lon(corner,grid1_add)
- lrevers = .true.
- endif
- begseg(1) = beglat
- begseg(2) = beglon
- lbegin = .true.
- num_subseg = 0
- !***
- !*** if this is a constant-longitude segment, skip the rest
- !*** since the line integral contribution will be zero.
- !***
- if (endlon /= beglon) then
- !***
- !*** integrate along this segment, detecting intersections
- !*** and computing the line integral for each sub-segment
- !***
- do while (beglat /= endlat .or. beglon /= endlon)
- !***
- !*** prevent infinite loops if integration gets stuck
- !*** near cell or threshold boundary
- !***
- num_subseg = num_subseg + 1
- if (num_subseg > max_subseg) then
- stop 'integration stalled: num_subseg exceeded limit'
- endif
- !***
- !*** find next intersection of this segment with a grid
- !*** line on grid 2.
- !***
- call timer_start(2)
- call intersection(grid2_add,intrsct_lat,intrsct_lon,lcoinc,
- & beglat, beglon, endlat, endlon, begseg,
- & lbegin, lrevers)
- call timer_stop(2)
- lbegin = .false.
- !***
- !*** compute line integral for this subsegment.
- !***
- call timer_start(3)
- if (grid2_add /= 0) then
- call line_integral(weights, num_wts,
- & beglon, intrsct_lon, beglat, intrsct_lat,
- & grid1_center_lat(grid1_add),
- & grid1_center_lon(grid1_add),
- & grid2_center_lat(grid2_add),
- & grid2_center_lon(grid2_add))
- else
- call line_integral(weights, num_wts,
- & beglon, intrsct_lon, beglat, intrsct_lat,
- & grid1_center_lat(grid1_add),
- & grid1_center_lon(grid1_add),
- & grid1_center_lat(grid1_add),
- & grid1_center_lon(grid1_add))
- endif
- call timer_stop(3)
- !***
- !*** if integrating in reverse order, change
- !*** sign of weights
- !***
- if (lrevers) then
- weights = -weights
- endif
- !***
- !*** store the appropriate addresses and weights.
- !*** also add contributions to cell areas and centroids.
- !***
- !if (grid1_add == 119247) then
- ! print *,grid1_add,grid2_add,corner,weights(1)
- ! print *,grid1_corner_lat(:,grid1_add)
- ! print *,grid1_corner_lon(:,grid1_add)
- ! print *,grid2_corner_lat(:,grid2_add)
- ! print *,grid2_corner_lon(:,grid2_add)
- ! print *,beglat,beglon,intrsct_lat,intrsct_lon
- !endif
- if (grid2_add /= 0) then
- if (grid1_mask(grid1_add)) then
- call timer_start(4)
- call store_link_cnsrv(grid1_add, grid2_add, weights)
- call timer_stop(4)
- grid1_frac(grid1_add) = grid1_frac(grid1_add) +
- & weights(1)
- grid2_frac(grid2_add) = grid2_frac(grid2_add) +
- & weights(num_wts+1)
- endif
- endif
- grid1_area(grid1_add) = grid1_area(grid1_add) + weights(1)
- grid1_centroid_lat(grid1_add) =
- & grid1_centroid_lat(grid1_add) + weights(2)
- grid1_centroid_lon(grid1_add) =
- & grid1_centroid_lon(grid1_add) + weights(3)
- !***
- !*** reset beglat and beglon for next subsegment.
- !***
- beglat = intrsct_lat
- beglon = intrsct_lon
- end do
- endif
- !***
- !*** end of segment
- !***
- end do
- !***
- !*** finished with this cell: deallocate search array and
- !*** start on next cell
- deallocate(srch_add, srch_corner_lat, srch_corner_lon)
- end do
- deallocate(srch_mask)
- !-----------------------------------------------------------------------
- !
- ! integrate around each cell on grid2
- !
- !-----------------------------------------------------------------------
- allocate(srch_mask(grid1_size))
- print *,'grid2 sweep '
- do grid2_add = 1,grid2_size
- !***
- !*** restrict searches first using search bins
- !***
- call timer_start(5)
- min_add = grid1_size
- max_add = 1
- do n=1,num_srch_bins
- if (grid2_add >= bin_addr2(1,n) .and.
- & grid2_add <= bin_addr2(2,n)) then
- min_add = min(min_add, bin_addr1(1,n))
- max_add = max(max_add, bin_addr1(2,n))
- endif
- end do
- !***
- !*** further restrict searches using bounding boxes
- !***
- num_srch_cells = 0
- do grid1_add = min_add, max_add
- srch_mask(grid1_add) = (grid1_bound_box(1,grid1_add) <=
- & grid2_bound_box(2,grid2_add)) .and.
- & (grid1_bound_box(2,grid1_add) >=
- & grid2_bound_box(1,grid2_add)) .and.
- & (grid1_bound_box(3,grid1_add) <=
- & grid2_bound_box(4,grid2_add)) .and.
- & (grid1_bound_box(4,grid1_add) >=
- & grid2_bound_box(3,grid2_add))
- if (srch_mask(grid1_add)) num_srch_cells = num_srch_cells+1
- end do
- allocate(srch_add(num_srch_cells),
- & srch_corner_lat(grid1_corners,num_srch_cells),
- & srch_corner_lon(grid1_corners,num_srch_cells))
- n = 0
- gather2: do grid1_add = min_add,max_add
- if (srch_mask(grid1_add)) then
- n = n+1
- srch_add(n) = grid1_add
- srch_corner_lat(:,n) = grid1_corner_lat(:,grid1_add)
- srch_corner_lon(:,n) = grid1_corner_lon(:,grid1_add)
- endif
- end do gather2
- call timer_stop(5)
- !***
- !*** integrate around this cell
- !***
- do corner = 1,grid2_corners
- next_corn = mod(corner,grid2_corners) + 1
- beglat = grid2_corner_lat(corner,grid2_add)
- beglon = grid2_corner_lon(corner,grid2_add)
- endlat = grid2_corner_lat(next_corn,grid2_add)
- endlon = grid2_corner_lon(next_corn,grid2_add)
- lrevers = .false.
- !***
- !*** to ensure exact path taken during both
- !*** sweeps, always integrate in the same direction
- !***
- if ((endlat < beglat) .or.
- & (endlat == beglat .and. endlon < beglon)) then
- beglat = grid2_corner_lat(next_corn,grid2_add)
- beglon = grid2_corner_lon(next_corn,grid2_add)
- endlat = grid2_corner_lat(corner,grid2_add)
- endlon = grid2_corner_lon(corner,grid2_add)
- lrevers = .true.
- endif
- begseg(1) = beglat
- begseg(2) = beglon
- lbegin = .true.
- !***
- !*** if this is a constant-longitude segment, skip the rest
- !*** since the line integral contribution will be zero.
- !***
- if (endlon /= beglon) then
- num_subseg = 0
- !***
- !*** integrate along this segment, detecting intersections
- !*** and computing the line integral for each sub-segment
- !***
- do while (beglat /= endlat .or. beglon /= endlon)
- !***
- !*** prevent infinite loops if integration gets stuck
- !*** near cell or threshold boundary
- !***
- num_subseg = num_subseg + 1
- if (num_subseg > max_subseg) then
- stop 'integration stalled: num_subseg exceeded limit'
- endif
- !***
- !*** find next intersection of this segment with a line
- !*** on grid 2.
- !***
- call timer_start(6)
- call intersection(grid1_add,intrsct_lat,intrsct_lon,lcoinc,
- & beglat, beglon, endlat, endlon, begseg,
- & lbegin, lrevers)
- call timer_stop(6)
- lbegin = .false.
- !***
- !*** compute line integral for this subsegment.
- !***
- call timer_start(7)
- if (grid1_add /= 0) then
- call line_integral(weights, num_wts,
- & beglon, intrsct_lon, beglat, intrsct_lat,
- & grid1_center_lat(grid1_add),
- & grid1_center_lon(grid1_add),
- & grid2_center_lat(grid2_add),
- & grid2_center_lon(grid2_add))
- else
- call line_integral(weights, num_wts,
- & beglon, intrsct_lon, beglat, intrsct_lat,
- & grid2_center_lat(grid2_add),
- & grid2_center_lon(grid2_add),
- & grid2_center_lat(grid2_add),
- & grid2_center_lon(grid2_add))
- endif
- call timer_stop(7)
- if (lrevers) then
- weights = -weights
- endif
- !***
- !*** store the appropriate addresses and weights.
- !*** also add contributions to cell areas and centroids.
- !*** if there is a coincidence, do not store weights
- !*** because they have been captured in the previous loop.
- !*** the grid1 mask is the master mask
- !***
- !if (grid1_add == 119247) then
- ! print *,grid1_add,grid2_add,corner,weights(1)
- ! print *,grid1_corner_lat(:,grid1_add)
- ! print *,grid1_corner_lon(:,grid1_add)
- ! print *,grid2_corner_lat(:,grid2_add)
- ! print *,grid2_corner_lon(:,grid2_add)
- ! print *,beglat,beglon,intrsct_lat,intrsct_lon
- !endif
- if (.not. lcoinc .and. grid1_add /= 0) then
- if (grid1_mask(grid1_add)) then
- call timer_start(8)
- call store_link_cnsrv(grid1_add, grid2_add, weights)
- call timer_stop(8)
- grid1_frac(grid1_add) = grid1_frac(grid1_add) +
- & weights(1)
- grid2_frac(grid2_add) = grid2_frac(grid2_add) +
- & weights(num_wts+1)
- endif
- endif
- grid2_area(grid2_add) = grid2_area(grid2_add) +
- & weights(num_wts+1)
- grid2_centroid_lat(grid2_add) =
- & grid2_centroid_lat(grid2_add) + weights(num_wts+2)
- grid2_centroid_lon(grid2_add) =
- & grid2_centroid_lon(grid2_add) + weights(num_wts+3)
- !***
- !*** reset beglat and beglon for next subsegment.
- !***
- beglat = intrsct_lat
- beglon = intrsct_lon
- end do
- endif
- !***
- !*** end of segment
- !***
- end do
- !***
- !*** finished with this cell: deallocate search array and
- !*** start on next cell
- deallocate(srch_add, srch_corner_lat, srch_corner_lon)
- end do
- deallocate(srch_mask)
- !-----------------------------------------------------------------------
- !
- ! correct for situations where N/S pole not explicitly included in
- ! grid (i.e. as a grid corner point). if pole is missing from only
- ! one grid, need to correct only the area and centroid of that
- ! grid. if missing from both, do complete weight calculation.
- !
- !-----------------------------------------------------------------------
- !*** North Pole
- weights(1) = pi2
- weights(2) = pi*pi
- weights(3) = zero
- weights(4) = pi2
- weights(5) = pi*pi
- weights(6) = zero
- grid1_add = 0
- pole_loop1: do n=1,grid1_size
- if (grid1_area(n) < -three*pih .and.
- & grid1_center_lat(n) > zero) then
- grid1_add = n
- exit pole_loop1
- endif
- end do pole_loop1
- grid2_add = 0
- pole_loop2: do n=1,grid2_size
- if (grid2_area(n) < -three*pih .and.
- & grid2_center_lat(n) > zero) then
- grid2_add = n
- exit pole_loop2
- endif
- end do pole_loop2
- if (grid1_add /=0) then
- grid1_area(grid1_add) = grid1_area(grid1_add) + weights(1)
- grid1_centroid_lat(grid1_add) =
- & grid1_centroid_lat(grid1_add) + weights(2)
- grid1_centroid_lon(grid1_add) =
- & grid1_centroid_lon(grid1_add) + weights(3)
- endif
- if (grid2_add /=0) then
- grid2_area(grid2_add) = grid2_area(grid2_add) +
- & weights(num_wts+1)
- grid2_centroid_lat(grid2_add) =
- & grid2_centroid_lat(grid2_add) + weights(num_wts+2)
- grid2_centroid_lon(grid2_add) =
- & grid2_centroid_lon(grid2_add) + weights(num_wts+3)
- endif
- if (grid1_add /= 0 .and. grid2_add /=0) then
- call store_link_cnsrv(grid1_add, grid2_add, weights)
- grid1_frac(grid1_add) = grid1_frac(grid1_add) +
- & weights(1)
- grid2_frac(grid2_add) = grid2_frac(grid2_add) +
- & weights(num_wts+1)
- endif
- !*** South Pole
- weights(1) = pi2
- weights(2) = -pi*pi
- weights(3) = zero
- weights(4) = pi2
- weights(5) = -pi*pi
- weights(6) = zero
- grid1_add = 0
- pole_loop3: do n=1,grid1_size
- if (grid1_area(n) < -three*pih .and.
- & grid1_center_lat(n) < zero) then
- grid1_add = n
- exit pole_loop3
- endif
- end do pole_loop3
- grid2_add = 0
- pole_loop4: do n=1,grid2_size
- if (grid2_area(n) < -three*pih .and.
- & grid2_center_lat(n) < zero) then
- grid2_add = n
- exit pole_loop4
- endif
- end do pole_loop4
- if (grid1_add /=0) then
- grid1_area(grid1_add) = grid1_area(grid1_add) + weights(1)
- grid1_centroid_lat(grid1_add) =
- & grid1_centroid_lat(grid1_add) + weights(2)
- grid1_centroid_lon(grid1_add) =
- & grid1_centroid_lon(grid1_add) + weights(3)
- endif
- if (grid2_add /=0) then
- grid2_area(grid2_add) = grid2_area(grid2_add) +
- & weights(num_wts+1)
- grid2_centroid_lat(grid2_add) =
- & grid2_centroid_lat(grid2_add) + weights(num_wts+2)
- grid2_centroid_lon(grid2_add) =
- & grid2_centroid_lon(grid2_add) + weights(num_wts+3)
- endif
- if (grid1_add /= 0 .and. grid2_add /=0) then
- call store_link_cnsrv(grid1_add, grid2_add, weights)
- grid1_frac(grid1_add) = grid1_frac(grid1_add) +
- & weights(1)
- grid2_frac(grid2_add) = grid2_frac(grid2_add) +
- & weights(num_wts+1)
- endif
- !-----------------------------------------------------------------------
- !
- ! finish centroid computation
- !
- !-----------------------------------------------------------------------
- where (grid1_area /= zero)
- grid1_centroid_lat = grid1_centroid_lat/grid1_area
- grid1_centroid_lon = grid1_centroid_lon/grid1_area
- end where
- where (grid2_area /= zero)
- grid2_centroid_lat = grid2_centroid_lat/grid2_area
- grid2_centroid_lon = grid2_centroid_lon/grid2_area
- end where
- !-----------------------------------------------------------------------
- !
- ! include centroids in weights and normalize using destination
- ! area if requested
- !
- !-----------------------------------------------------------------------
- do n=1,num_links_map1
- grid1_add = grid1_add_map1(n)
- grid2_add = grid2_add_map1(n)
- do nwgt=1,num_wts
- weights( nwgt) = wts_map1(nwgt,n)
- if (num_maps > 1) then
- weights(num_wts+nwgt) = wts_map2(nwgt,n)
- endif
- end do
- select case(norm_opt)
- case (norm_opt_dstarea)
- if (grid2_area(grid2_add) /= zero) then
- if (luse_grid2_area) then
- norm_factor = one/grid2_area_in(grid2_add)
- else
- norm_factor = one/grid2_area(grid2_add)
- endif
- else
- norm_factor = zero
- endif
- case (norm_opt_frcarea)
- if (grid2_frac(grid2_add) /= zero) then
- if (luse_grid2_area) then
- norm_factor = grid2_area(grid2_add)/
- & (grid2_frac(grid2_add)*
- & grid2_area_in(grid2_add))
- else
- norm_factor = one/grid2_frac(grid2_add)
- endif
- else
- norm_factor = zero
- endif
- case (norm_opt_none)
- norm_factor = one
- end select
- wts_map1(1,n) = weights(1)*norm_factor
- wts_map1(2,n) = (weights(2) - weights(1)*
- & grid1_centroid_lat(grid1_add))*
- & norm_factor
- wts_map1(3,n) = (weights(3) - weights(1)*
- & grid1_centroid_lon(grid1_add))*
- & norm_factor
- if (num_maps > 1) then
- select case(norm_opt)
- case (norm_opt_dstarea)
- if (grid1_area(grid1_add) /= zero) then
- if (luse_grid1_area) then
- norm_factor = one/grid1_area_in(grid1_add)
- else
- norm_factor = one/grid1_area(grid1_add)
- endif
- else
- norm_factor = zero
- endif
- case (norm_opt_frcarea)
- if (grid1_frac(grid1_add) /= zero) then
- if (luse_grid1_area) then
- norm_factor = grid1_area(grid1_add)/
- & (grid1_frac(grid1_add)*
- & grid1_area_in(grid1_add))
- else
- norm_factor = one/grid1_frac(grid1_add)
- endif
- else
- norm_factor = zero
- endif
- case (norm_opt_none)
- norm_factor = one
- end select
- wts_map2(1,n) = weights(num_wts+1)*norm_factor
- wts_map2(2,n) = (weights(num_wts+2) - weights(num_wts+1)*
- & grid2_centroid_lat(grid2_add))*
- & norm_factor
- wts_map2(3,n) = (weights(num_wts+3) - weights(num_wts+1)*
- & grid2_centroid_lon(grid2_add))*
- & norm_factor
- endif
- end do
- print *, 'Total number of links = ',num_links_map1
- where (grid1_area /= zero) grid1_frac = grid1_frac/grid1_area
- where (grid2_area /= zero) grid2_frac = grid2_frac/grid2_area
- !-----------------------------------------------------------------------
- !
- ! perform some error checking on final weights
- !
- !-----------------------------------------------------------------------
- grid2_centroid_lat = zero
- grid2_centroid_lon = zero
- do n=1,grid1_size
- if (grid1_area(n) < -.01) then
- print *,'Grid 1 area error: ',n,grid1_area(n)
- endif
- if (grid1_centroid_lat(n) < -pih-.01 .or.
- & grid1_centroid_lat(n) > pih+.01) then
- print *,'Grid 1 centroid lat error: ',n,grid1_centroid_lat(n)
- endif
- grid1_centroid_lat(n) = zero
- grid1_centroid_lon(n) = zero
- end do
- do n=1,grid2_size
- if (grid2_area(n) < -.01) then
- print *,'Grid 2 area error: ',n,grid2_area(n)
- endif
- if (grid2_centroid_lat(n) < -pih-.01 .or.
- & grid2_centroid_lat(n) > pih+.01) then
- print *,'Grid 2 centroid lat error: ',n,grid2_centroid_lat(n)
- endif
- grid2_centroid_lat(n) = zero
- grid2_centroid_lon(n) = zero
- end do
- do n=1,num_links_map1
- grid1_add = grid1_add_map1(n)
- grid2_add = grid2_add_map1(n)
-
- if (wts_map1(1,n) < -.01) then
- print *,'Map 1 weight < 0 ',grid1_add,grid2_add,wts_map1(1,n)
- endif
- if (norm_opt /= norm_opt_none .and. wts_map1(1,n) > 1.01) then
- print *,'Map 1 weight > 1 ',grid1_add,grid2_add,wts_map1(1,n)
- endif
- grid2_centroid_lat(grid2_add) =
- & grid2_centroid_lat(grid2_add) + wts_map1(1,n)
- if (num_maps > 1) then
- if (wts_map2(1,n) < -.01) then
- print *,'Map 2 weight < 0 ',grid1_add,grid2_add,
- & wts_map2(1,n)
- endif
- if (norm_opt /= norm_opt_none .and. wts_map2(1,n) > 1.01) then
- print *,'Map 2 weight < 0 ',grid1_add,grid2_add,
- & wts_map2(1,n)
- endif
- grid1_centroid_lat(grid1_add) =
- & grid1_centroid_lat(grid1_add) + wts_map2(1,n)
- endif
- end do
- do n=1,grid2_size
- select case(norm_opt)
- case (norm_opt_dstarea)
- norm_factor = grid2_frac(grid2_add)
- case (norm_opt_frcarea)
- norm_factor = one
- case (norm_opt_none)
- if (luse_grid2_area) then
- norm_factor = grid2_area_in(grid2_add)
- else
- norm_factor = grid2_area(grid2_add)
- endif
- end select
- if (abs(grid2_centroid_lat(grid2_add)-norm_factor) > .01) then
- print *,'Error: sum of wts for map1 ',grid2_add,
- & grid2_centroid_lat(grid2_add),norm_factor
- endif
- end do
- if (num_maps > 1) then
- do n=1,grid1_size
- select case(norm_opt)
- case (norm_opt_dstarea)
- norm_factor = grid1_frac(grid1_add)
- case (norm_opt_frcarea)
- norm_factor = one
- case (norm_opt_none)
- if (luse_grid1_area) then
- norm_factor = grid1_area_in(grid1_add)
- else
- norm_factor = grid1_area(grid1_add)
- endif
- end select
- if (abs(grid1_centroid_lat(grid1_add)-norm_factor) > .01) then
- print *,'Error: sum of wts for map2 ',grid1_add,
- & grid1_centroid_lat(grid1_add),norm_factor
- endif
- end do
- endif
- !-----------------------------------------------------------------------
- end subroutine remap_conserv
- !***********************************************************************
- subroutine intersection(location,intrsct_lat,intrsct_lon,lcoinc,
- & beglat, beglon, endlat, endlon, begseg,
- & lbegin, lrevers)
- !-----------------------------------------------------------------------
- !
- ! this routine finds the next intersection of a destination grid
- ! line with the line segment given by beglon, endlon, etc.
- ! a coincidence flag is returned if the segment is entirely
- ! coincident with an ocean grid line. the cells in which to search
- ! for an intersection must have already been restricted in the
- ! calling routine.
- !
- !-----------------------------------------------------------------------
- !-----------------------------------------------------------------------
- !
- ! intent(in):
- !
- !-----------------------------------------------------------------------
- logical (kind=log_kind), intent(in) ::
- & lbegin, ! flag for first integration along this segment
- & lrevers ! flag whether segment integrated in reverse
- real (kind=dbl_kind), intent(in) ::
- & beglat, beglon, ! beginning lat/lon endpoints for segment
- & endlat, endlon ! ending lat/lon endpoints for segment
- real (kind=dbl_kind), dimension(2), intent(inout) ::
- & begseg ! begin lat/lon of full segment
- !-----------------------------------------------------------------------
- !
- ! intent(out):
- !
- !-----------------------------------------------------------------------
- integer (kind=int_kind), intent(out) ::
- & location ! address in destination array containing this
- ! segment
- logical (kind=log_kind), intent(out) ::
- & lcoinc ! flag segments which are entirely coincident
- ! with a grid line
- real (kind=dbl_kind), intent(out) ::
- & intrsct_lat, intrsct_lon ! lat/lon coords of next intersect.
- !-----------------------------------------------------------------------
- !
- ! local variables
- !
- !-----------------------------------------------------------------------
- integer (kind=int_kind) :: n, next_n, cell, srch_corners, pole_loc
- integer (kind=int_kind), save ::
- & last_loc ! save location when crossing threshold
- logical (kind=log_kind) ::
- & loutside ! flags points outside grid
- logical (kind=log_kind), save ::
- & lthresh = .false. ! flags segments crossing threshold bndy
- real (kind=dbl_kind) ::
- & lon1, lon2, ! local longitude variables for segment
- & lat1, lat2, ! local latitude variables for segment
- & grdlon1, grdlon2, ! local longitude variables for grid cell
- & grdlat1, grdlat2, ! local latitude variables for grid cell
- & vec1_lat, vec1_lon, ! vectors and cross products used
- & vec2_lat, vec2_lon, ! during grid search
- & cross_product,
- & eps, offset, ! small offset away from intersect
- & s1, s2, determ, ! variables used for linear solve to
- & mat1, mat2, mat3, mat4, rhs1, rhs2 ! find intersection
- real (kind=dbl_kind), save ::
- & intrsct_lat_off, intrsct_lon_off ! lat/lon coords offset
- ! for next search
- !-----------------------------------------------------------------------
- !
- ! initialize defaults, flags, etc.
- !
- !-----------------------------------------------------------------------
- location = 0
- lcoinc = .false.
- intrsct_lat = endlat
- intrsct_lon = endlon
- if (num_srch_cells == 0) return
- if (beglat > north_thresh .or. beglat < south_thresh) then
- if (lthresh) location = last_loc
- call pole_intersection(location,
- & intrsct_lat,intrsct_lon,lcoinc,lthresh,
- & beglat, beglon, endlat, endlon, begseg, lrevers)
- if (lthresh) then
- last_loc = location
- intrsct_lat_off = intrsct_lat
- intrsct_lon_off = intrsct_lon
- endif
- return
- endif
- loutside = .false.
- if (lbegin) then
- lat1 = beglat
- lon1 = beglon
- else
- lat1 = intrsct_lat_off
- lon1 = intrsct_lon_off
- endif
- lat2 = endlat
- lon2 = endlon
- if ((lon2-lon1) > three*pih) then
- lon2 = lon2 - pi2
- else if ((lon2-lon1) < -three*pih) then
- lon2 = lon2 + pi2
- endif
- s1 = zero
- !-----------------------------------------------------------------------
- !
- ! search for location of this segment in ocean grid using cross
- ! product method to determine whether a point is enclosed by a cell
- !
- !-----------------------------------------------------------------------
- call timer_start(12)
- srch_corners = size(srch_corner_lat,DIM=1)
- srch_loop: do
- !***
- !*** if last segment crossed threshold, use that location
- !***
- if (lthresh) then
- do cell=1,num_srch_cells
- if (srch_add(cell) == last_loc) then
- location = last_loc
- eps = tiny
- exit srch_loop
- endif
- end do
- endif
- !***
- !*** otherwise normal search algorithm
- !***
- cell_loop: do cell=1,num_srch_cells
- corner_loop: do n=1,srch_corners
- next_n = MOD(n,srch_corners) + 1
- !***
- !*** here we take the cross product of the vector making
- !*** up each cell side with the vector formed by the vertex
- !*** and search point. if all the cross products are
- !*** positive, the point is contained in the cell.
- !***
- vec1_lat = srch_corner_lat(next_n,cell) -
- & srch_corner_lat(n ,cell)
- vec1_lon = srch_corner_lon(next_n,cell) -
- & srch_corner_lon(n ,cell)
- vec2_lat = lat1 - srch_corner_lat(n,cell)
- vec2_lon = lon1 - srch_corner_lon(n,cell)
- !***
- !*** if endpoint coincident with vertex, offset
- !*** the endpoint
- !***
- if (vec2_lat == 0 .and. vec2_lon == 0) then
- lat1 = lat1 + 1.d-10*(lat2-lat1)
- lon1 = lon1 + 1.d-10*(lon2-lon1)
- vec2_lat = lat1 - srch_corner_lat(n,cell)
- vec2_lon = lon1 - srch_corner_lon(n,cell)
- endif
- !***
- !*** check for 0,2pi crossings
- !***
- if (vec1_lon > pi) then
- vec1_lon = vec1_lon - pi2
- else if (vec1_lon < -pi) then
- vec1_lon = vec1_lon + pi2
- endif
- if (vec2_lon > pi) then
- vec2_lon = vec2_lon - pi2
- else if (vec2_lon < -pi) then
- vec2_lon = vec2_lon + pi2
- endif
- cross_product = vec1_lon*vec2_lat - vec2_lon*vec1_lat
- !***
- !*** if the cross product for a side is zero, the point
- !*** lies exactly on the side or the side is degenerate
- !*** (zero length). if degenerate, set the cross
- !*** product to a positive number. otherwise perform
- !*** another cross product between the side and the
- !*** segment itself.
- !*** if this cross product is also zero, the line is
- !*** coincident with the cell boundary - perform the
- !*** dot product and only choose the cell if the dot
- !*** product is positive (parallel vs anti-parallel).
- !***
- if (cross_product == zero) then
- if (vec1_lat /= zero .or. vec1_lon /= zero) then
- vec2_lat = lat2 - lat1
- vec2_lon = lon2 - lon1
- if (vec2_lon > pi) then
- vec2_lon = vec2_lon - pi2
- else if (vec2_lon < -pi) then
- vec2_lon = vec2_lon + pi2
- endif
- cross_product = vec1_lon*vec2_lat - vec2_lon*vec1_lat
- else
- cross_product = one
- endif
- if (cross_product == zero) then
- lcoinc = .true.
- cross_product = vec1_lon*vec2_lon + vec1_lat*vec2_lat
- if (lrevers) cross_product = -cross_product
- endif
- endif
- !***
- !*** if cross product is less than zero, this cell
- !*** doesn't work
- !***
- if (cross_product < zero) exit corner_loop
- end do corner_loop
- !***
- !*** if cross products all positive, we found the location
- !***
- if (n > srch_corners) then
- location = srch_add(cell)
- !***
- !*** if the beginning of this segment was outside the
- !*** grid, invert the segment so the intersection found
- !*** will be the first intersection with the grid
- !***
- if (loutside) then
- lat2 = beglat
- lon2 = beglon
- location = 0
- eps = -tiny
- else
- eps = tiny
- endif
- exit srch_loop
- endif
- !***
- !*** otherwise move on to next cell
- !***
- end do cell_loop
- !***
- !*** if still no cell found, the point lies outside the grid.
- !*** take some baby steps along the segment to see if any
- !*** part of the segment lies inside the grid.
- !***
- loutside = .true.
- s1 = s1 + 0.001_dbl_kind
- lat1 = beglat + s1*(endlat - beglat)
- lon1 = beglon + s1*(lon2 - beglon)
- !***
- !*** reached the end of the segment and still outside the grid
- !*** return no intersection
- !***
- if (s1 >= one) return
- end do srch_loop
- call timer_stop(12)
- !-----------------------------------------------------------------------
- !
- ! now that a cell is found, search for the next intersection.
- ! loop over sides of the cell to find intersection with side
- ! must check all sides for coincidences or intersections
- !
- !-----------------------------------------------------------------------
- call timer_start(13)
- intrsct_loop: do n=1,srch_corners
- next_n = mod(n,srch_corners) + 1
- grdlon1 = srch_corner_lon(n ,cell)
- grdlon2 = srch_corner_lon(next_n,cell)
- grdlat1 = srch_corner_lat(n ,cell)
- grdlat2 = srch_corner_lat(next_n,cell)
- !***
- !*** set up linear system to solve for intersection
- !***
- mat1 = lat2 - lat1
- mat2 = grdlat1 - grdlat2
- mat3 = lon2 - lon1
- mat4 = grdlon1 - grdlon2
- rhs1 = grdlat1 - lat1
- rhs2 = grdlon1 - lon1
- if (mat3 > pi) then
- mat3 = mat3 - pi2
- else if (mat3 < -pi) then
- mat3 = mat3 + pi2
- endif
- if (mat4 > pi) then
- mat4 = mat4 - pi2
- else if (mat4 < -pi) then
- mat4 = mat4 + pi2
- endif
- if (rhs2 > pi) then
- rhs2 = rhs2 - pi2
- else if (rhs2 < -pi) then
- rhs2 = rhs2 + pi2
- endif
- determ = mat1*mat4 - mat2*mat3
- !***
- !*** if the determinant is zero, the segments are either
- !*** parallel or coincident. coincidences were detected
- !*** above so do nothing.
- !*** if the determinant is non-zero, solve for the linear
- !*** parameters s for the intersection point on each line
- !*** segment.
- !*** if 0<s1,s2<1 then the segment intersects with this side.
- !*** return the point of intersection (adding a small
- !*** number so the intersection is off the grid line).
- !***
- if (abs(determ) > 1.e-30) then
- s1 = (rhs1*mat4 - mat2*rhs2)/determ
- s2 = (mat1*rhs2 - rhs1*mat3)/determ
- if (s2 >= zero .and. s2 <= one .and.
- & s1 > zero. and. s1 <= one) then
- !***
- !*** recompute intersection based on full segment
- !*** so intersections are consistent for both sweeps
- !***
- if (.not. loutside) then
- mat1 = lat2 - begseg(1)
- mat3 = lon2 - begseg(2)
- rhs1 = grdlat1 - begseg(1)
- rhs2 = grdlon1 - begseg(2)
- else
- mat1 = begseg(1) - endlat
- mat3 = begseg(2) - endlon
- rhs1 = grdlat1 - endlat
- rhs2 = grdlon1 - endlon
- endif
- if (mat3 > pi) then
- mat3 = mat3 - pi2
- else if (mat3 < -pi) then
- mat3 = mat3 + pi2
- endif
- if (rhs2 > pi) then
- rhs2 = rhs2 - pi2
- else if (rhs2 < -pi) then
- rhs2 = rhs2 + pi2
- endif
- determ = mat1*mat4 - mat2*mat3
- !***
- !*** sometimes due to roundoff, the previous
- !*** determinant is non-zero, but the lines
- !*** are actually coincident. if this is the
- !*** case, skip the rest.
- !***
- if (determ /= zero) then
- s1 = (rhs1*mat4 - mat2*rhs2)/determ
- s2 = (mat1*rhs2 - rhs1*mat3)/determ
- offset = s1 + eps/determ
- if (offset > one) offset = one
- if (.not. loutside) then
- intrsct_lat = begseg(1) + mat1*s1
- intrsct_lon = begseg(2) + mat3*s1
- intrsct_lat_off = begseg(1) + mat1*offset
- intrsct_lon_off = begseg(2) + mat3*offset
- else
- intrsct_lat = endlat + mat1*s1
- intrsct_lon = endlon + mat3*s1
- intrsct_lat_off = endlat + mat1*offset
- intrsct_lon_off = endlon + mat3*offset
- endif
- exit intrsct_loop
- endif
- endif
- endif
- !***
- !*** no intersection this side, move on to next side
- !***
- end do intrsct_loop
- call timer_stop(13)
- !-----------------------------------------------------------------------
- !
- ! if the segment crosses a pole threshold, reset the intersection
- ! to be the threshold latitude. only check if this was not a
- ! threshold segment since sometimes coordinate transform can end
- ! up on other side of threshold again.
- !
- !-----------------------------------------------------------------------
- if (lthresh) then
- if (intrsct_lat < north_thresh .or. intrsct_lat > south_thresh)
- & lthresh = .false.
- else if (lat1 > zero .and. intrsct_lat > north_thresh) then
- intrsct_lat = north_thresh + tiny
- intrsct_lat_off = north_thresh + eps*mat1
- s1 = (intrsct_lat - begseg(1))/mat1
- intrsct_lon = begseg(2) + s1*mat3
- intrsct_lon_off = begseg(2) + (s1+eps)*mat3
- last_loc = location
- lthresh = .true.
- else if (lat1 < zero .and. intrsct_lat < south_thresh) then
- intrsct_lat = south_thresh - tiny
- intrsct_lat_off = south_thresh + eps*mat1
- s1 = (intrsct_lat - begseg(1))/mat1
- intrsct_lon = begseg(2) + s1*mat3
- intrsct_lon_off = begseg(2) + (s1+eps)*mat3
- last_loc = location
- lthresh = .true.
- endif
- !-----------------------------------------------------------------------
- end subroutine intersection
- !***********************************************************************
- subroutine pole_intersection(location,
- & intrsct_lat,intrsct_lon,lcoinc,lthresh,
- & beglat, beglon, endlat, endlon, begseg, lrevers)
- !-----------------------------------------------------------------------
- !
- ! this routine is identical to the intersection routine except
- ! that a coordinate transformation (using a Lambert azimuthal
- ! equivalent projection) is performed to treat polar cells more
- ! accurately.
- !
- !-----------------------------------------------------------------------
- !-----------------------------------------------------------------------
- !
- ! intent(in):
- !
- !-----------------------------------------------------------------------
- real (kind=dbl_kind), intent(in) ::
- & beglat, beglon, ! beginning lat/lon endpoints for segment
- & endlat, endlon ! ending lat/lon endpoints for segment
- real (kind=dbl_kind), dimension(2), intent(inout) ::
- & begseg ! begin lat/lon of full segment
- logical (kind=log_kind), intent(in) ::
- & lrevers ! flag true if segment integrated in reverse
- !-----------------------------------------------------------------------
- !
- ! intent(out):
- !
- !-----------------------------------------------------------------------
- integer (kind=int_kind), intent(inout) ::
- & location ! address in destination array containing this
- ! segment -- also may contain last location on
- ! entry
- logical (kind=log_kind), intent(out) ::
- & lcoinc ! flag segment coincident with grid line
- logical (kind=log_kind), intent(inout) ::
- & lthresh ! flag segment crossing threshold boundary
- real (kind=dbl_kind), intent(out) ::
- & intrsct_lat, intrsct_lon ! lat/lon coords of next intersect.
- !-----------------------------------------------------------------------
- !
- ! local variables
- !
- !-----------------------------------------------------------------------
- integer (kind=int_kind) :: n, next_n, cell, srch_corners, pole_loc
- logical (kind=log_kind) :: loutside ! flags points outside grid
- real (kind=dbl_kind) :: pi4, rns, ! north/south conversion
- & x1, x2, ! local x variables for segment
- & y1, y2, ! local y variables for segment
- & begx, begy, ! beginning x,y variables for segment
- & endx, endy, ! beginning x,y variables for segment
- & begsegx, begsegy, ! beginning x,y variables for segment
- & grdx1, grdx2, ! local x variables for grid cell
- & grdy1, grdy2, ! local y variables for grid cell
- & vec1_y, vec1_x, ! vectors and cross products used
- & vec2_y, vec2_x, ! during grid search
- & cross_product, eps, ! eps=small offset away from intersect
- & s1, s2, determ, ! variables used for linear solve to
- & mat1, mat2, mat3, mat4, rhs1, rhs2 ! find intersection
- real (kind=dbl_kind), dimension(:,:), allocatable ::
- & srch_corner_x, ! x of each corner of srch cells
- & srch_corner_y ! y of each corner of srch cells
- !***
- !*** save last intersection to avoid roundoff during coord
- !*** transformation
- !***
- logical (kind=log_kind), save :: luse_last = .false.
- real (kind=dbl_kind), save ::
- & intrsct_x, intrsct_y ! x,y for intersection
- !***
- !*** variables necessary if segment manages to hit pole
- !***
- integer (kind=int_kind), save ::
- & avoid_pole_count = 0 ! count attempts to avoid pole
- real (kind=dbl_kind), save ::
- & avoid_pole_offset = tiny ! endpoint offset to avoid pole
- !-----------------------------------------------------------------------
- !
- ! initialize defaults, flags, etc.
- !
- !-----------------------------------------------------------------------
- if (.not. lthresh) location = 0
- lcoinc = .false.
- intrsct_lat = endlat
- intrsct_lon = endlon
- loutside = .false.
- s1 = zero
- !-----------------------------------------------------------------------
- !
- ! convert coordinates
- !
- !-----------------------------------------------------------------------
- allocate(srch_corner_x(size(srch_corner_lat,DIM=1),
- & size(srch_corner_lat,DIM=2)),
- & srch_corner_y(size(srch_corner_lat,DIM=1),
- & size(srch_corner_lat,DIM=2)))
- if (beglat > zero) then
- pi4 = quart*pi
- rns = one
- else
- pi4 = -quart*pi
- rns = -one
- endif
- if (luse_last) then
- x1 = intrsct_x
- y1 = intrsct_y
- else
- x1 = rns*two*sin(pi4 - half*beglat)*cos(beglon)
- y1 = two*sin(pi4 - half*beglat)*sin(beglon)
- luse_last = .true.
- endif
- x2 = rns*two*sin(pi4 - half*endlat)*cos(endlon)
- y2 = two*sin(pi4 - half*endlat)*sin(endlon)
- srch_corner_x = rns*two*sin(pi4 - half*srch_corner_lat)*
- & cos(srch_corner_lon)
- srch_corner_y = two*sin(pi4 - half*srch_corner_lat)*
- & sin(srch_corner_lon)
- begx = x1
- begy = y1
- endx = x2
- endy = y2
- begsegx = rns*two*sin(pi4 - half*begseg(1))*cos(begseg(2))
- begsegy = two*sin(pi4 - half*begseg(1))*sin(begseg(2))
- intrsct_x = endx
- intrsct_y = endy
- !-----------------------------------------------------------------------
- !
- ! search for location of this segment in ocean grid using cross
- ! product method to determine whether a point is enclosed by a cell
- !
- !-----------------------------------------------------------------------
- call timer_start(12)
- srch_corners = size(srch_corner_lat,DIM=1)
- srch_loop: do
- !***
- !*** if last segment crossed threshold, use that location
- !***
- if (lthresh) then
- do cell=1,num_srch_cells
- if (srch_add(cell) == location) then
- eps = tiny
- exit srch_loop
- endif
- end do
- endif
- !***
- !*** otherwise normal search algorithm
- !***
- cell_loop: do cell=1,num_srch_cells
- corner_loop: do n=1,srch_corners
- next_n = MOD(n,srch_corners) + 1
- !***
- !*** here we take the cross product of the vector making
- !*** up each cell side with the vector formed by the vertex
- !*** and search point. if all the cross products are
- !*** positive, the point is contained in the cell.
- !***
- vec1_x = srch_corner_x(next_n,cell) -
- & srch_corner_x(n ,cell)
- vec1_y = srch_corner_y(next_n,cell) -
- & srch_corner_y(n ,cell)
- vec2_x = x1 - srch_corner_x(n,cell)
- vec2_y = y1 - srch_corner_y(n,cell)
- !***
- !*** if endpoint coincident with vertex, offset
- !*** the endpoint
- !***
- if (vec2_x == 0 .and. vec2_y == 0) then
- x1 = x1 + 1.d-10*(x2-x1)
- y1 = y1 + 1.d-10*(y2-y1)
- vec2_x = x1 - srch_corner_x(n,cell)
- vec2_y = y1 - srch_corner_y(n,cell)
- endif
- cross_product = vec1_x*vec2_y - vec2_x*vec1_y
- !***
- !*** if the cross product for a side is zero, the point
- !*** lies exactly on the side or the length of a side
- !*** is zero. if the length is zero set det > 0.
- !*** otherwise, perform another cross
- !*** product between the side and the segment itself.
- !*** if this cross product is also zero, the line is
- !*** coincident with the cell boundary - perform the
- !*** dot product and only choose the cell if the dot
- !*** product is positive (parallel vs anti-parallel).
- !***
- if (cross_product == zero) then
- if (vec1_x /= zero .or. vec1_y /= 0) then
- vec2_x = x2 - x1
- vec2_y = y2 - y1
- cross_product = vec1_x*vec2_y - vec2_x*vec1_y
- else
- cross_product = one
- endif
- if (cross_product == zero) then
- lcoinc = .true.
- cross_product = vec1_x*vec2_x + vec1_y*vec2_y
- if (lrevers) cross_product = -cross_product
- endif
- endif
- !***
- !*** if cross product is less than zero, this cell
- !*** doesn't work
- !***
- if (cross_product < zero) exit corner_loop
- end do corner_loop
- !***
- !*** if cross products all positive, we found the location
- !***
- if (n > srch_corners) then
- location = srch_add(cell)
- !***
- !*** if the beginning of this segment was outside the
- !*** grid, invert the segment so the intersection found
- !*** will be the first intersection with the grid
- !***
- if (loutside) then
- x2 = begx
- y2 = begy
- location = 0
- eps = -tiny
- else
- eps = tiny
- endif
- exit srch_loop
- endif
- !***
- !*** otherwise move on to next cell
- !***
- end do cell_loop
- !***
- !*** if no cell found, the point lies outside the grid.
- !*** take some baby steps along the segment to see if any
- !*** part of the segment lies inside the grid.
- !***
- loutside = .true.
- s1 = s1 + 0.001_dbl_kind
- x1 = begx + s1*(x2 - begx)
- y1 = begy + s1*(y2 - begy)
- !***
- !*** reached the end of the segment and still outside the grid
- !*** return no intersection
- !***
- if (s1 >= one) then
- deallocate(srch_corner_x, srch_corner_y)
- luse_last = .false.
- return
- endif
- end do srch_loop
- call timer_stop(12)
- !-----------------------------------------------------------------------
- !
- ! now that a cell is found, search for the next intersection.
- ! loop over sides of the cell to find intersection with side
- ! must check all sides for coincidences or intersections
- !
- !-----------------------------------------------------------------------
- call timer_start(13)
- intrsct_loop: do n=1,srch_corners
- next_n = mod(n,srch_corners) + 1
- grdy1 = srch_corner_y(n ,cell)
- grdy2 = srch_corner_y(next_n,cell)
- grdx1 = srch_corner_x(n ,cell)
- grdx2 = srch_corner_x(next_n,cell)
- !***
- !*** set up linear system to solve for intersection
- !***
- mat1 = x2 - x1
- mat2 = grdx1 - grdx2
- mat3 = y2 - y1
- mat4 = grdy1 - grdy2
- rhs1 = grdx1 - x1
- rhs2 = grdy1 - y1
- determ = mat1*mat4 - mat2*mat3
- !***
- !*** if the determinant is zero, the segments are either
- !*** parallel or coincident or one segment has zero length.
- !*** coincidences were detected above so do nothing.
- !*** if the determinant is non-zero, solve for the linear
- !*** parameters s for the intersection point on each line
- !*** segment.
- !*** if 0<s1,s2<1 then the segment intersects with this side.
- !*** return the point of intersection (adding a small
- !*** number so the intersection is off the grid line).
- !***
- if (abs(determ) > 1.e-30) then
- s1 = (rhs1*mat4 - mat2*rhs2)/determ
- s2 = (mat1*rhs2 - rhs1*mat3)/determ
- if (s2 >= zero .and. s2 <= one .and.
- & s1 > zero. and. s1 <= one) then
- !***
- !*** recompute intersection using entire segment
- !*** for consistency between sweeps
- !***
- if (.not. loutside) then
- mat1 = x2 - begsegx
- mat3 = y2 - begsegy
- rhs1 = grdx1 - begsegx
- rhs2 = grdy1 - begsegy
- else
- mat1 = x2 - endx
- mat3 = y2 - endy
- rhs1 = grdx1 - endx
- rhs2 = grdy1 - endy
- endif
- determ = mat1*mat4 - mat2*mat3
- !***
- !*** sometimes due to roundoff, the previous
- !*** determinant is non-zero, but the lines
- !*** are actually coincident. if this is the
- !*** case, skip the rest.
- !***
- if (determ /= zero) then
- s1 = (rhs1*mat4 - mat2*rhs2)/determ
- s2 = (mat1*rhs2 - rhs1*mat3)/determ
- if (.not. loutside) then
- intrsct_x = begsegx + s1*mat1
- intrsct_y = begsegy + s1*mat3
- else
- intrsct_x = endx + s1*mat1
- intrsct_y = endy + s1*mat3
- endif
- !***
- !*** convert back to lat/lon coordinates
- !***
- intrsct_lon = rns*atan2(intrsct_y,intrsct_x)
- if (intrsct_lon < zero)
- & intrsct_lon = intrsct_lon + pi2
- if (abs(intrsct_x) > 1.d-10) then
- intrsct_lat = (pi4 -
- & asin(rns*half*intrsct_x/cos(intrsct_lon)))*two
- else if (abs(intrsct_y) > 1.d-10) then
- intrsct_lat = (pi4 -
- & asin(half*intrsct_y/sin(intrsct_lon)))*two
- else
- intrsct_lat = two*pi4
- endif
- !***
- !*** add offset in transformed space for next pass.
- !***
- if (s1 - eps/determ < one) then
- intrsct_x = intrsct_x - mat1*(eps/determ)
- intrsct_y = intrsct_y - mat3*(eps/determ)
- else
- if (.not. loutside) then
- intrsct_x = endx
- intrsct_y = endy
- intrsct_lat = endlat
- intrsct_lon = endlon
- else
- intrsct_x = begsegx
- intrsct_y = begsegy
- intrsct_lat = begseg(1)
- intrsct_lon = begseg(2)
- endif
- endif
- exit intrsct_loop
- endif
- endif
- endif
- !***
- !*** no intersection this side, move on to next side
- !***
- end do intrsct_loop
- call timer_stop(13)
- deallocate(srch_corner_x, srch_corner_y)
- !-----------------------------------------------------------------------
- !
- ! if segment manages to cross over pole, shift the beginning
- ! endpoint in order to avoid hitting pole directly
- ! (it is ok for endpoint to be pole point)
- !
- !-----------------------------------------------------------------------
- if (abs(intrsct_x) < 1.e-10 .and. abs(intrsct_y) < 1.e-10 .and.
- & (endx /= zero .and. endy /=0)) then
- if (avoid_pole_count > 2) then
- avoid_pole_count = 0
- avoid_pole_offset = 10.*avoid_pole_offset
- endif
- cross_product = begsegx*(endy-begsegy) - begsegy*(endx-begsegx)
- intrsct_lat = begseg(1)
- if (cross_product*intrsct_lat > zero) then
- intrsct_lon = beglon + avoid_pole_offset
- begseg(2) = begseg(2) + avoid_pole_offset
- else
- intrsct_lon = beglon - avoid_pole_offset
- begseg(2) = begseg(2) - avoid_pole_offset
- endif
- avoid_pole_count = avoid_pole_count + 1
- luse_last = .false.
- else
- avoid_pole_count = 0
- avoid_pole_offset = tiny
- endif
- !-----------------------------------------------------------------------
- !
- ! if the segment crosses a pole threshold, reset the intersection
- ! to be the threshold latitude and do not reuse x,y intersect
- ! on next entry. only check if did not cross threshold last
- ! time - sometimes the coordinate transformation can place a
- ! segment on the other side of the threshold again
- !
- !-----------------------------------------------------------------------
- if (lthresh) then
- if (intrsct_lat > north_thresh .or. intrsct_lat < south_thresh)
- & lthresh = .false.
- else if (beglat > zero .and. intrsct_lat < north_thresh) then
- mat4 = endlat - begseg(1)
- mat3 = endlon - begseg(2)
- if (mat3 > pi) mat3 = mat3 - pi2
- if (mat3 < -pi) mat3 = mat3 + pi2
- intrsct_lat = north_thresh - tiny
- s1 = (north_thresh - begseg(1))/mat4
- intrsct_lon = begseg(2) + s1*mat3
- luse_last = .false.
- lthresh = .true.
- else if (beglat < zero .and. intrsct_lat > south_thresh) then
- mat4 = endlat - begseg(1)
- mat3 = endlon - begseg(2)
- if (mat3 > pi) mat3 = mat3 - pi2
- if (mat3 < -pi) mat3 = mat3 + pi2
- intrsct_lat = south_thresh + tiny
- s1 = (south_thresh - begseg(1))/mat4
- intrsct_lon = begseg(2) + s1*mat3
- luse_last = .false.
- lthresh = .true.
- endif
- !***
- !*** if reached end of segment, do not use x,y intersect
- !*** on next entry
- !***
- if (intrsct_lat == endlat .and. intrsct_lon == endlon) then
- luse_last = .false.
- endif
- !-----------------------------------------------------------------------
- end subroutine pole_intersection
- !***********************************************************************
- subroutine line_integral(weights, num_wts,
- & in_phi1, in_phi2, theta1, theta2,
- & grid1_lat, grid1_lon, grid2_lat, grid2_lon)
- !-----------------------------------------------------------------------
- !
- ! this routine computes the line integral of the flux function
- ! that results in the interpolation weights. the line is defined
- ! by the input lat/lon of the endpoints.
- !
- !-----------------------------------------------------------------------
- !-----------------------------------------------------------------------
- !
- ! intent(in):
- !
- !-----------------------------------------------------------------------
- integer (kind=int_kind), intent(in) ::
- & num_wts ! number of weights to compute
- real (kind=dbl_kind), intent(in) ::
- & in_phi1, in_phi2, ! longitude endpoints for the segment
- & theta1, theta2, ! latitude endpoints for the segment
- & grid1_lat, grid1_lon, ! reference coordinates for each
- & grid2_lat, grid2_lon ! grid (to ensure correct 0,2pi interv.
- !-----------------------------------------------------------------------
- !
- ! intent(out):
- !
- !-----------------------------------------------------------------------
- real (kind=dbl_kind), dimension(2*num_wts), intent(out) ::
- & weights ! line integral contribution to weights
- !-----------------------------------------------------------------------
- !
- ! local variables
- !
- !-----------------------------------------------------------------------
- real (kind=dbl_kind) :: dphi, sinth1, sinth2, costh1, costh2, fac,
- & phi1, phi2, phidiff1, phidiff2, sinint
- real (kind=dbl_kind) :: f1, f2, fint
- !-----------------------------------------------------------------------
- !
- ! weights for the general case based on a trapezoidal approx to
- ! the integrals.
- !
- !-----------------------------------------------------------------------
- sinth1 = SIN(theta1)
- sinth2 = SIN(theta2)
- costh1 = COS(theta1)
- costh2 = COS(theta2)
- dphi = in_phi1 - in_phi2
- if (dphi > pi) then
- dphi = dphi - pi2
- else if (dphi < -pi) then
- dphi = dphi + pi2
- endif
- dphi = half*dphi
- !-----------------------------------------------------------------------
- !
- ! the first weight is the area overlap integral. the second and
- ! fourth are second-order latitude gradient weights.
- !
- !-----------------------------------------------------------------------
- weights( 1) = dphi*(sinth1 + sinth2)
- weights(num_wts+1) = dphi*(sinth1 + sinth2)
- weights( 2) = dphi*(costh1 + costh2 + (theta1*sinth1 +
- & theta2*sinth2))
- weights(num_wts+2) = dphi*(costh1 + costh2 + (theta1*sinth1 +
- & theta2*sinth2))
- !-----------------------------------------------------------------------
- !
- ! the third and fifth weights are for the second-order phi gradient
- ! component. must be careful of longitude range.
- !
- !-----------------------------------------------------------------------
- f1 = half*(costh1*sinth1 + theta1)
- f2 = half*(costh2*sinth2 + theta2)
- phi1 = in_phi1 - grid1_lon
- if (phi1 > pi) then
- phi1 = phi1 - pi2
- else if (phi1 < -pi) then
- phi1 = phi1 + pi2
- endif
- phi2 = in_phi2 - grid1_lon
- if (phi2 > pi) then
- phi2 = phi2 - pi2
- else if (phi2 < -pi) then
- phi2 = phi2 + pi2
- endif
- if ((phi2-phi1) < pi .and. (phi2-phi1) > -pi) then
- weights(3) = dphi*(phi1*f1 + phi2*f2)
- else
- if (phi1 > zero) then
- fac = pi
- else
- fac = -pi
- endif
- fint = f1 + (f2-f1)*(fac-phi1)/abs(dphi)
- weights(3) = half*phi1*(phi1-fac)*f1 -
- & half*phi2*(phi2+fac)*f2 +
- & half*fac*(phi1+phi2)*fint
- endif
- phi1 = in_phi1 - grid2_lon
- if (phi1 > pi) then
- phi1 = phi1 - pi2
- else if (phi1 < -pi) then
- phi1 = phi1 + pi2
- endif
- phi2 = in_phi2 - grid2_lon
- if (phi2 > pi) then
- phi2 = phi2 - pi2
- else if (phi2 < -pi) then
- phi2 = phi2 + pi2
- endif
- if ((phi2-phi1) < pi .and. (phi2-phi1) > -pi) then
- weights(num_wts+3) = dphi*(phi1*f1 + phi2*f2)
- else
- if (phi1 > zero) then
- fac = pi
- else
- fac = -pi
- endif
- fint = f1 + (f2-f1)*(fac-phi1)/abs(dphi)
- weights(num_wts+3) = half*phi1*(phi1-fac)*f1 -
- & half*phi2*(phi2+fac)*f2 +
- & half*fac*(phi1+phi2)*fint
- endif
- !-----------------------------------------------------------------------
- end subroutine line_integral
- !***********************************************************************
- subroutine store_link_cnsrv(add1, add2, weights)
- !-----------------------------------------------------------------------
- !
- ! this routine stores the address and weight for this link in
- ! the appropriate address and weight arrays and resizes those
- ! arrays if necessary.
- !
- !-----------------------------------------------------------------------
- !-----------------------------------------------------------------------
- !
- ! input variables
- !
- !-----------------------------------------------------------------------
- integer (kind=int_kind), intent(in) ::
- & add1, ! address on grid1
- & add2 ! address on grid2
- real (kind=dbl_kind), dimension(:), intent(in) ::
- & weights ! array of remapping weights for this link
- !-----------------------------------------------------------------------
- !
- ! local variables
- !
- !-----------------------------------------------------------------------
- integer (kind=int_kind) :: nlink, min_link, max_link ! link index
- integer (kind=int_kind), dimension(:,:), allocatable, save ::
- & link_add1, ! min,max link add to restrict search
- & link_add2 ! min,max link add to restrict search
- logical (kind=log_kind), save :: first_call = .true.
- !-----------------------------------------------------------------------
- !
- ! if all weights are zero, do not bother storing the link
- !
- !-----------------------------------------------------------------------
- if (all(weights == zero)) return
- !-----------------------------------------------------------------------
- !
- ! restrict the range of links to search for existing links
- !
- !-----------------------------------------------------------------------
- if (first_call) then
- allocate(link_add1(2,grid1_size), link_add2(2,grid2_size))
- link_add1 = 0
- link_add2 = 0
- first_call = .false.
- min_link = 1
- max_link = 0
- else
- min_link = min(link_add1(1,add1),link_add2(1,add2))
- max_link = max(link_add1(2,add1),link_add2(2,add2))
- if (min_link == 0) then
- min_link = 1
- max_link = 0
- endif
- endif
- !-----------------------------------------------------------------------
- !
- ! if the link already exists, add the weight to the current weight
- ! arrays
- !
- !-----------------------------------------------------------------------
- do nlink=min_link,max_link
- if (add1 == grid1_add_map1(nlink)) then
- if (add2 == grid2_add_map1(nlink)) then
- wts_map1(:,nlink) = wts_map1(:,nlink) + weights(1:num_wts)
- if (num_maps == 2) then
- wts_map2(:,nlink) = wts_map2(:,nlink) +
- & weights(num_wts+1:2*num_wts)
- endif
- return
- endif
- endif
- end do
- !-----------------------------------------------------------------------
- !
- ! if the link does not yet exist, increment number of links and
- ! check to see if remap arrays need to be increased to accomodate
- ! the new link. then store the link.
- !
- !-----------------------------------------------------------------------
- num_links_map1 = num_links_map1 + 1
- if (num_links_map1 > max_links_map1)
- & call resize_remap_vars(1,resize_increment)
- grid1_add_map1(num_links_map1) = add1
- grid2_add_map1(num_links_map1) = add2
- wts_map1 (:,num_links_map1) = weights(1:num_wts)
- if (num_maps > 1) then
- num_links_map2 = num_links_map2 + 1
- if (num_links_map2 > max_links_map2)
- & call resize_remap_vars(2,resize_increment)
- grid1_add_map2(num_links_map2) = add1
- grid2_add_map2(num_links_map2) = add2
- wts_map2 (:,num_links_map2) = weights(num_wts+1:2*num_wts)
- endif
- if (link_add1(1,add1) == 0) link_add1(1,add1) = num_links_map1
- if (link_add2(1,add2) == 0) link_add2(1,add2) = num_links_map1
- link_add1(2,add1) = num_links_map1
- link_add2(2,add2) = num_links_map1
- !-----------------------------------------------------------------------
- end subroutine store_link_cnsrv
- !***********************************************************************
- end module remap_conservative
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|