Module cleanArray

implicit none

contains

subroutine sweep(bathymetry,sodadepth,maxpoints,FILL,KK,VV,datain,dataout,mask,II,JJ)

! ------------------------------------------------------------------

! Program : sweep

!

! This routine finds points in the interpolated fields that was masked by

! the original land-array, but which is no longer masked in the new land-ocean array.

! Sometimes this leads to ocean posint in the new grid that are still

! defined as land because of the previous masking. This routine sweeps through

! the array and replaces fill values over ocean with the nearest points values.

!

! The routine finds the points in KKxVV distance from missing value points.

! This means that a maximum of KK points in the horizontal direction and a maximum

! of VV points in the vertical are searched for values to be used to fill in

! the missing value. A maximum of "maxpoints" within the distance KK and VV are

! then used to make sure you use mostly values close to missing value. However, with the

! ability to have high values of KK and VV you make sure that a value is

! found to fill in if no-values exists nearby. The values are weighted

! with the distance from the missing value point. If no value is found within KKxVV distance,

! the search radius is increased (INCREASED=1) and for each iteration no values are found,

! the KK and VV distance increases with 10 points. As soon as a good value is found,

! the search radius decreases back to the original size determined in grd.py. This

! increase/decrease is necessary for some points in difficult topographical areas (like Bay of Fundy).

! The program increases the searh radius 3 times before giving up to find values.

!

! Trond Kristiansen, March and April 2009

! Rutgers University, NJ.

!

! Compile with:

! f2py-64 --verbose -DF2PY_REPORT_ON_ARRAY_COPY=1 --fcompiler=intelem -c -m cl cleanArray.f90

! or

! f2py --verbose --fcompiler=intelem -c -m cl cleanArray.f90

!

! USAGE:

! import cl

! Zin = cl.cleanArray.sweep(Zin,Zout, maskarray,i,j,k)

!

! (I,J+1)

! |----------------O----------------|

! | | |

! | | |

! (I-1,J) X (I,J) (I+1,J)

! | | |

! | | |

! |----------------O----------------|

! (I,J-1)

!

! From each of these 6 points (5 unique) we move up and down VV points,

! and KK points horisontally. A maximum of "maxpoints" are stored in each direction (foundData array)

! together with the distance from the missing value point (foundWeight array)

! ------------------------------------------------------------------

integer maxpoints, II, JJ, KK, VV, ic, jc, jcc, icc, fill, kc, ff, total

integer pointsHpos, pointsHneg, pointsVpos, pointsVneg, points, INCREASED, orgVV, orgKK, TRIES

REAL(4), dimension(JJ,II) :: dataout,datain

REAL(4), dimension(JJ,II) :: mask

REAL(4), dimension(4,3,maxpoints) :: foundWeight, foundData

REAL(4) counter, hori, vert, sumweights, distance, sodadepth, weight

REAL(4), dimension(JJ,II) :: bathymetry

! print*,'---> Started cleaning of array with box size of',KK,VV

total = 0

orgVV=VV

orgKK=KK

do jc=1,JJ

do ic=1,II

TRIES=1

counter = 0

pointsHpos = 1

pointsVpos = 1

pointsHneg = 1

pointsVneg = 1

if (abs(datain(jc,ic)) > FILL .AND. mask(jc,ic)==1 .AND. sodadepth >= -(bathymetry(jc,ic))) THEN

! Have to check only grid cells that are deeper than the sodadepth. If not, we will fill in all

! values toward the coast for each iteration (sodadepth >= -(bathymetry(jc,ic)))).

total=total+1

foundWeight(:,:,:)=0

foundData(:,:,:)=0

50 continue

do ff=-1,1,1

do kc=1,KK

hori=ic+kc

if (hori .LE. 1 ) then

hori=1

else if (hori .GE. II) then

hori=II

end if

vert=jc+ff

if (vert .LE. 1 ) then

vert=1

else if (vert .GE. JJ) then

vert=JJ

end if

if (abs(datain(vert,hori)) < FILL .AND. mask(vert,hori)==1 .AND. &

& foundWeight(1,ff+2,pointsHneg) .EQ. 0 .AND. pointsHneg .LE. maxpoints) THEN

counter = counter + 1

foundData(1,ff+2,pointsHneg)=datain(vert,hori)

foundWeight(1,ff+2,pointsHneg)=sqrt(((hori) - ic)**2.0 + ((vert) - jc)**2.0)

pointsHneg=pointsHneg+1

end if

hori=ic+kc*(-1)

if (hori .LE. 1 ) then

hori=1

else if (hori .GE. II) then

hori=II

end if

if (abs(datain(vert,hori)) < FILL .AND. mask(vert,hori)==1 .AND. &

& foundWeight(3,ff+2,pointsHpos) .EQ. 0 .AND. pointsHpos .LE. maxpoints) THEN

counter = counter + 1

foundData(3,ff+2,pointsHpos)=datain(vert,hori)

foundWeight(3,ff+2,pointsHpos)=sqrt(((hori) - ic)**2.0 + ((vert) - jc)**2.0)

pointsHpos=pointsHpos+1

end if

end do

end do

do ff=-1,1,1

do kc=1,VV

! Vertically find values around problem value

vert=jc+kc

if (vert .LE. 1 ) then

vert=1

else if (vert .GE. JJ) then

vert=JJ

end if

hori=ic+ff

if (hori .LE. 1 ) then

hori=1

else if (hori .GE. II) then

hori=II

end if

if (abs(datain(vert,hori)) < FILL .AND. mask(vert,hori)==1 &

& .AND. foundWeight(2,ff+2, pointsVneg) .EQ. 0 .AND. pointsVneg .LE. maxpoints) THEN

counter = counter + 1

foundData(2,ff+2, pointsVneg)=datain(vert,hori)

foundWeight(2,ff+2, pointsVneg)=sqrt(((hori) - ic)**2.0 + ((vert) - jc)**2.0)

pointsVneg=pointsVneg+1

end if

vert=jc+kc*(-1)

if (vert .LE. 1 ) then

vert=1

else if (vert .GE. JJ) then

vert=JJ

end if

if (abs(datain(vert,hori)) < FILL .AND. mask(vert,hori)==1 &

& .AND. foundWeight(4,ff+2, pointsVpos) .EQ. 0 .AND. pointsVpos .LE. maxpoints) THEN

counter = counter + 1

foundData(4,ff+2, pointsVpos)=datain(vert,hori)

foundWeight(4,ff+2, pointsVpos)=sqrt(((ic+ff) - ic)**2.0 + ((vert) - jc)**2.0)

pointsVpos=pointsVpos+1

end if

end do

end do

! Assign the final value to data array using weights

if (counter > 0) then

distance = sum(abs(foundWeight))

sumweights=0.0

!print*,'total number of points',counter

dataout(jc,ic)=0.0

do kc=1,4

do ff=1,3

do points=1,maxpoints

if (foundWeight(kc,ff,points) .NE. 0) then

if (counter > 1) then

weight=((1.0-(abs(foundWeight(kc,ff,points)))/distance)/(counter-1.0)*1.0)

sumweights=sumweights + ((1.0-(abs(foundWeight(kc,ff,points)))/distance)/(counter-1.0)*1.0)

else

weight=1.0

sumweights=sumweights + 1.0

end if

dataout(jc,ic) = dataout(jc,ic) + foundData(kc,ff,points)*weight

!print*,foundData(kc,ff,points), distance, abs(foundWeight(kc,ff,points))/distance

!print*,'total',sumweights, distance, 1-(abs(foundWeight(kc,ff,points))

end if

end do

end do

end do

end if

if (counter==0) then

!print*,'no data', dataout(jc,ic), mask(jc,ic), jc,ic

! Here we increase the search area as no points within KKxVV

! both vertically and horisontally was found.

KK=KK+10

VV=VV+10

INCREASED=1

TRIES=TRIES+1

!print*,'Increased distances for values', TRIES, KK, VV

if (TRIES .LT. 3) then

goto 50

else

goto 100

end if

else if (counter >0 .AND. INCREASED .EQ. 1) then

! Good values found, now reset the search area to

! what is defined in grd.py

INCREASED=0

KK=orgKK

VV=orgVV

!print*,'reduced distances', jc, ic, KK, VV

end if

end if

100 continue

end do

end do

!print*,'Number of points filled', total

! ----

end subroutine sweep

end module cleanArray