module mutil use mkinds USE mconstants contains REAL (r4) FUNCTION ridx(ain,val) implicit none !Function returns the index of the cell containing the value val. real(r4),intent(in) :: val real(r4),dimension(:),intent(in) :: ain ! local variables integer(r4) :: k,n ! treat boundaries first - not good -- needs better checks ! IF(val>=ain(size(ain))) then ridx=size(ain) return ENDIF IF(val<=ain(1)) then ridx=1.0 return ENDIF ! normal case ! search and find the location in the grid ! full search need to implement exit if found and error handling DO k = 2,size(ain) IF (ain(k-1) .le. val .and. ain(k) .ge. val) then ridx = real(k-1) + (val-ain(k-1)) / (ain(k)-ain(k-1)) return ENDIF ENDDO end function ridx REAL (r4) FUNCTION rridx(ain,val) implicit none !Function returns the index of the cell containing the value val. real(r4),intent(in) :: val real(r4),dimension(:),intent(in) :: ain ! local variables integer(r4) :: k ! search and find the location in the grid ! full search need to implement exit if found and error handling DO k = 2, size(ain) IF (val.ge.ain(k-1) .and. val .le. ain(k)) then rridx = real(k-1) + (val-ain(k-1)) / (ain(k)-ain(k-1)) return ENDIF ENDDO end function rridx subroutine lon_lat_idx(mlon,mlat,nx,ny,lon,lat,lglb,il,jl) implicit none !Function returns the index of the cell containing the value val. integer(r4) :: nx,ny real(r4),intent(in) :: lon,lat real(r4),dimension(:,:),intent(in) :: mlon,mlat logical(bl) :: lglb real(r4) :: il,jl if(lglb) then if(lat>47.0) then call exsearch(mlon,mlat,nx,ny,lon,lat,il,jl) else jl=ridx(mlat(1,:),lat) il=glb_lon_idx(mlon(:,1),lon) endif else jl=ridx(mlat(1,:),lat) il=ridx(mlon(:,1),lon) endif END Subroutine REAL (r4) FUNCTION glb_lon_idx(ain,val) implicit none !Function returns the index of the cell containing the value val. real(r4),intent(in) :: val real(r4),dimension(:),intent(in) :: ain ! local variables integer(r4) :: k ! search and find the location in the grid ! full search need to implement exit if found and error handling if(val>-180.0 .and. val<=74.0) then DO k = 442, size(ain) IF (val.ge.ain(k-1) .and. val .le. ain(k)) then glb_lon_idx = real(k-1) + (val-ain(k-1)) / (ain(k)-ain(k-1)) return ENDIF ENDDO else DO k = 2, size(ain) IF (val.ge.ain(k-1) .and. val .le. ain(k)) then glb_lon_idx = real(k-1) + (val-ain(k-1)) / (ain(k)-ain(k-1)) return ENDIF ENDDO endif glb_lon_idx=max(1.,glb_lon_idx) glb_lon_idx=min(real(size(ain),r4),glb_lon_idx) end function glb_lon_idx subroutine intp_fld (src_nx,src_ny,dst_nx,dst_ny,src_lon,src_lat,dst_lon,dst_lat,src_fld,dst_fld) USE mkinds IMPLICIT NONE INTEGER(KIND=i4) :: src_nx,src_ny INTEGER(KIND=i4) :: dst_nx,dst_ny logical(bl) :: lperiodic ! source and destination arrays REAL(KIND=r4) :: src_lon(src_nx) REAL(KIND=r4) :: src_lat(src_ny) REAL(KIND=r4) :: src_fld(src_nx,src_ny) REAL(KIND=r4) :: dst_fld(dst_nx,dst_ny) REAL(KIND=r4) :: dst_lon(dst_nx,dst_ny) REAL(KIND=r4) :: dst_lat(dst_nx,dst_ny) ! Internal arrays INTEGER(KIND=i4) :: xmap(dst_nx,dst_ny,1,4) INTEGER(KIND=i4) :: ymap(dst_nx,dst_ny,1,4) REAL(KIND=r4) :: wgts(dst_nx,dst_ny,1,4) ! internal variables INTEGER(KIND=i4) :: ii INTEGER(KIND=i4) :: jj REAL(KIND=r4) :: ptlon REAL(KIND=r4) :: ptlat REAL(KIND=r4) :: grdi REAL(KIND=r4) :: grdj REAL(KIND=r4) :: r REAL(KIND=r4) :: s REAL(KIND=r4) :: xsum DO jj=1,dst_ny DO ii=1,dst_nx ptlon=dst_lon(ii,jj) ptlat=dst_lat(ii,jj) CALL ll2ij(src_lon,src_lat,ptlon,ptlat,src_nx,src_ny,grdi,grdj) ! print *,ptlon,ptlat, grdi,grdj CALL setwgts(src_nx,src_ny,grdi,grdj,xmap(ii,jj,1,:),ymap(ii,jj,1,:),wgts(ii,jj,1,:)) ! dst_fld(ii,jj)=src_fld(int(grdi),int(grdj)) IF(sum(wgts(ii,jj,:,1)) > 1.01) THEN WRITE(*,*) "Weights Error: ", ii,jj,grdi,grdj ENDIF ENDDO ENDDO !if(pid==0) WRITE(*,*)"Min values of mapping coodinates in OBS Grid: ", src_lat(MINVAL(ymap,mask=ymap.gt.0)),src_lon(MINVAL(xmap,mask=xmap.gt.0)) !if(pid==0)WRITE(*,*) "Max values of mapping coodinates in OBS Grid: ", src_lat(MAXVAL(ymap,mask=ymap.gt.0)),src_lon(MAXVAL(xmap,mask=xmap.gt.0)) ! Now map the source grid to destination grid CALL remap(dst_fld(:,:),src_fld(:,:),xmap(:,:,1,:),ymap(:,:,1,:),wgts(:,:,1,:),src_nx,src_ny,dst_nx,dst_ny) end subroutine REAL(KIND=r4) FUNCTION sphdist(rlon1,rlat1, rlon2, rlat2) !...Returns distance between the points [rlon1 rlat1] and [rlon2 rlat2] in kilometers IMPLICIT NONE REAL(r4), INTENT(IN) :: rlon1,rlat1, rlon2,rlat2 ! Local variables REAL(r4) :: rln1,rlt1, rln2,rlt2,dist rln1=deg2rad*rlon1; rlt1=deg2rad*rlat1; rln2=deg2rad*rlon2; rlt2=deg2rad*rlat2; dist= cos(rlt1)*cos(rlt2)*cos(rln2-rln1) + sin(rlt1)*sin(rlt2) IF (abs(dist) < 1.0) THEN sphdist = REarth * acos(dist) ELSE IF ( (dist-1.0) < TOLERANCE) THEN sphdist = 0.0 ELSE !WRITE(6,*)'***Error dist****=', dist !WRITE(6,*)'lon1, lat1 =',rlon1, rlat1 !WRITE(6,*)'lon2, lat2 =', rlon2, rlat2 ENDIF ENDIF RETURN END FUNCTION sphdist SUBROUTINE blinterp(fld,nx,ny,xi,yj,ifld) IMPLICIT NONE INTEGER (i4),INTENT(IN) :: nx,ny REAL(r4),INTENT(IN) :: fld(nx,ny),xi,yj REAL(r4),INTENT(OUT) :: ifld ! local variables INTEGER(i4) :: ix,jy REAL(r4) :: r,s,f,f1,f2,f3 !Assumes that all points comming in are on the grid -- needs error handling ix = int (xi) jy = int (yj) if (ix>1 .and. ix1 .and. jy1e10) ifld=0.0 END SUBROUTINE blinterp SUBROUTINE blintProfile(fld,nz,xi,yj,ifld) IMPLICIT NONE INTEGER (i4),INTENT(IN) :: nz REAL(r4),INTENT(IN) :: fld(4,nz),xi,yj REAL(r4),INTENT(OUT) :: ifld(nz) ! local variables INTEGER(i4) :: ix,jy,k REAL(r4) :: r,s,f,f1,f2,f3 !Assumes that all points comming in are on the grid -- needs error handling ix = int (xi) jy = int (yj) r = xi - real (ix) s = yj - real (jy) ! bilinear intp do k=1,nz f = fld(1,k) f1 = fld(2,k) f2 = fld(3,k) f3 = fld(4,k) ifld(k) = (1. - s) * ((1. - r) * f + r * f1) + s * ((1. - r) * f2 + r * f3) enddo END SUBROUTINE blintprofile SUBROUTINE all2ij (glon,glat,aptlon,aptlat,nlon,nlat,npts,grdi,grdj) IMPLICIT NONE INTEGER(i4), INTENT(IN) :: nlon, nlat, npts REAL(r4), INTENT(IN) :: glon(nlon), glat(nlat) REAL(r4), INTENT(IN) :: aptlon(npts),aptlat(npts) REAL(r4), INTENT(OUT) :: grdi(npts),grdj(npts) ! local variables INTEGER(i4) :: k,j REAL(r4) :: ptlon,ptlat do j=1,npts ptlon=aptlon(j) ptlat=aptlat(j) DO k = 2, nlat IF (glat(k-1) .le. ptlat .and. glat(k) .ge. ptlat) then grdj(j) = real(k-1) + (ptlat-glat(k-1)) / (glat(k)-glat(k-1)) exit ENDIF ENDDO DO k = 2, nlon IF (ptlon .ge. glon(k-1) .and. ptlon .le. glon(k)) then grdi(j) = real(k-1) + (ptlon - glon(k-1)) / (glon(k)- glon(k-1)) exit ENDIF ENDDO enddo END SUBROUTINE all2ij !SUBROUTINE ll2ij (glon,glat,ptlon,ptlat,nlon,nlat,lperiodbc,grdi,grdj) !USE mkinds !IMPLICIT NONE !INTEGER(i4), INTENT(IN) :: nlon, nlat !REAL(r4), INTENT(IN) :: glon(nlon), glat(nlat) !REAL(r4), INTENT(IN) :: ptlon,ptlat !REAL(r4), INTENT(OUT) :: grdi,grdj !LOGICAL (bl) :: lperiodbc ! local variables !INTEGER (i4) :: k !LOGICAL (bl) :: novalue ! DO k = 2, nlat ! IF (glat(k-1) .le. ptlat .and. glat(k) .ge. ptlat) then ! grdj = real(k-1) + (ptlat-glat(k-1)) / (glat(k)-glat(k-1)) ! ENDIF ! ENDDO ! novalue = .true. ! DO k = 2, nlon !special case if ptlon is between for example 359.9 and 0.1 degrees. ! ! IF (glon(k-1) .GT. glon(k)) THEN !print *,"here", k,glon(k-1),glon(k), ! stop ! IF (ptlon .gt. (glon(k-1)-360.) .and. ptlon .le. glon(k)) THEN ! novalue = .false. ! grdi = real(k-1) + (ptlon - (glon(k-1)-360.)) /(glon(k)-(glon(k-1)-360.)) ! ENDIF ! IF (ptlon-360 .gt. (glon(k-1)-360.) .and. ptlon-360 .le. glon(k)) THEN ! novalue = .false. ! grdi = real(k-1) + (ptlon-360 - (glon(k-1)-360.)) /(glon(k)-(glon(k-1)-360.)) ! ENDIF ! ! ELSE ! normal case ! IF (ptlon .ge. glon(k-1) .and. ptlon .le. glon(k)) then ! novalue = .false. ! grdi = real(k-1) + (ptlon - glon(k-1)) / (glon(k)- glon(k-1)) ! ENDIF ! ENDIF ! ENDDO ! if check periodic boundary conditions is set to true ! and grdi has no value yet, calculate value of grdi ! IF ((lperiodbc) .and. (novalue)) THEN ! IF (ptlon .gt. glon(nlon)) THEN ! grdi = nlon + (ptlon - glon(nlon)) / (glon(1) - glon(nlon)) ! ELSE ! grdi = 0. + (ptlon - glon(nlon)) / (glon(1) - glon(nlon)) ! ENDIF !print *, grdi ! ENDIF !END SUBROUTINE ll2ij SUBROUTINE ll2ij(glon,glat,ptlon,ptlat,nlon,nlat,grdi,grdj) IMPLICIT NONE INTEGER(i4), INTENT(IN) :: nlon, nlat REAL(r4), INTENT(IN) :: glon(nlon), glat(nlat) REAL(r4), INTENT(IN) :: ptlon,ptlat REAL(r4), INTENT(OUT) :: grdi,grdj ! local variables INTEGER (i4) :: k DO k = 2, nlat IF (glat(k-1) .le. ptlat .and. glat(k) .ge. ptlat) then grdj = real(k-1) + (ptlat-glat(k-1)) / (glat(k)-glat(k-1)) exit ENDIF ENDDO DO k = 2, nlon IF (ptlon .ge. glon(k-1) .and. ptlon .le. glon(k)) then grdi = real(k-1) + (ptlon - glon(k-1)) / (glon(k)- glon(k-1)) exit ENDIF ENDDO if(ptlon>glon(nlon)) grdi=1!nlon if(ptlonglat(nlat)) grdj=1!nlat if(ptlat0.) wrk(i,j)=num/den else wrk(i,j) = spval endif END DO END DO fld=wrk(1:nx,1:ny) end subroutine SUBROUTINE remap(dst_array,src_array,xmap,ymap,wgts,src_nx,src_ny,dst_nx,dst_ny) use mkinds IMPLICIT NONE INTEGER(i4),INTENT(IN) :: dst_nx,dst_ny,src_nx,src_ny REAL(r4),INTENT(INOUT) :: dst_array(dst_nx,dst_ny) REAL(r4),INTENT(IN) :: src_array(src_nx,src_ny) INTEGER(i4),INTENT(IN) :: xmap(dst_nx,dst_ny,1,4) INTEGER(i4),INTENT(IN) :: ymap(dst_nx,dst_ny,1,4) REAL(r4),INTENT(IN) :: wgts(dst_nx,dst_ny,1,4) ! Local Variables INTEGER(i4) :: i,j DO j=1,dst_ny DO i=1,dst_nx dst_array(i,j)= wgts(i,j,1,1)*src_array(xmap(i,j,1,1),ymap(i,j,1,1)) & + wgts(i,j,1,2)*src_array(xmap(i,j,1,2),ymap(i,j,1,2)) & + wgts(i,j,1,3)*src_array(xmap(i,j,1,3),ymap(i,j,1,3)) & + wgts(i,j,1,4)*src_array(xmap(i,j,1,4),ymap(i,j,1,4)) ENDDO ENDDO END SUBROUTINE SUBROUTINE setwgts(src_nx,src_ny,grdi,grdj,xmap,ymap,wgts) USE mkinds IMPLICIT NONE REAL(r4) :: grdi, grdj INTEGER(i4) :: src_nx,src_ny INTEGER(i4) :: xmap(4),ymap(4) REAL(r4) :: wgts(4) ! local variables INTEGER(i4) :: ix,iy REAL(r4) :: r,s r=grdi-real(int(grdi)) s=grdj-real(int(grdj)) ix=int(grdi) iy=int(grdj) if(ix>1 .and. ix<=src_nx-1 .and. iy>1. .and. iy<=src_ny-1) then ! bilinear weights ! lower left point xmap(1)=int(grdi) ymap(1)=int(grdj) wgts(1)=(1.-r)*(1.-s) ! lower right xmap(2)=int(grdi)+1 ymap(2)=int(grdj) wgts(2)=r*(1.-s) ! upper right xmap(3)=int(grdi)+1 ymap(3)=int(grdj)+1 wgts(3)=r*s ! upper left xmap(4)=int(grdi) ymap(4)=int(grdj)+1 wgts(4)=(1.-r)*s else xmap(1:4)=int(grdi) ymap(1:4)=int(grdj) wgts=0.25 endif END SUBROUTINE setwgts REAL(r4) Function tinterp(var1,var2,w1,w2) implicit none real(r4), intent(in) :: var1 real(r4), intent(in) :: var2 real(r4), intent(in) :: w1 real(r4), intent(in) :: w2 tinterp=w1*var1+w2*var2 END Function integer(i4) Function ispresent(alon,alat,fcount,lon,lat) implicit none integer(r4) :: fcount real(r4) :: alon(fcount) real(r4) :: alat(fcount) real(r4) :: lon real(r4) :: lat integer(i4) :: i,present_lon,present_lat if(fcount==0) then ispresent=0 return endif do i=1,fcount if(alon(i)==lon .and. alat(i)==lat) then ispresent=1 exit else ispresent=0 endif enddo end Function integer(i4) Function chk_llbox(lon,lat,lnmn,lnmx,ltmn,ltmx) implicit none real(r4) :: lon real(r4) :: lat real(r4) :: lnmn real(r4) :: lnmx real(r4) :: ltmn real(r4) :: ltmx if(lon> lnmn .and. lonltmn .and. lat= N. ! N AND X ARE NOT ALTERED BY THIS ROUTINE. ! OUTPUT PARAMETER - IND - SEQUENCE OF INDICES 1,...,N PERMUTED IN THE SAME ! FASHION AS X WOULD BE. THUS, THE ORDERING ON ! X IS DEFINED BY Y(I) = X(IND(I)). !********************************************************************* ! NOTE -- IU AND IL MUST BE DIMENSIONED >= LOG(N) WHERE LOG HAS BASE 2. !********************************************************************* INTEGER(i4) :: iu(21), il(21) INTEGER(i4) :: m, i, j, k, l, ij, it, itt, indx REAL(r4) :: r REAL(r4) :: t ! LOCAL PARAMETERS - ! IU,IL = TEMPORARY STORAGE FOR THE UPPER AND LOWER ! INDICES OF PORTIONS OF THE ARRAY X ! M = INDEX FOR IU AND IL ! I,J = LOWER AND UPPER INDICES OF A PORTION OF X ! K,L = INDICES IN THE RANGE I,...,J ! IJ = RANDOMLY CHOSEN INDEX BETWEEN I AND J ! IT,ITT = TEMPORARY STORAGE FOR INTERCHANGES IN IND ! INDX = TEMPORARY INDEX FOR X ! R = PSEUDO RANDOM NUMBER FOR GENERATING IJ ! T = CENTRAL ELEMENT OF X IF (n <= 0) RETURN ! INITIALIZE IND, M, I, J, AND R DO i = 1, n ind(i) = i END DO m = 1 i = 1 j = n r = .375 ! TOP OF LOOP 20 IF (i >= j) GO TO 70 IF (r <= .5898437) THEN r = r + .0390625 ELSE r = r - .21875 END IF ! INITIALIZE K 30 k = i ! SELECT A CENTRAL ELEMENT OF X AND SAVE IT IN T ij = i + r*(j-i) it = ind(ij) t = x(it) ! IF THE FIRST ELEMENT OF THE ARRAY IS GREATER THAN T, ! INTERCHANGE IT WITH T indx = ind(i) IF (x(indx) > t) THEN ind(ij) = indx ind(i) = it it = indx t = x(it) END IF ! INITIALIZE L l = j ! IF THE LAST ELEMENT OF THE ARRAY IS LESS THAN T, ! INTERCHANGE IT WITH T indx = ind(j) IF (x(indx) >= t) GO TO 50 ind(ij) = indx ind(j) = it it = indx t = x(it) ! IF THE FIRST ELEMENT OF THE ARRAY IS GREATER THAN T, ! INTERCHANGE IT WITH T indx = ind(i) IF (x(indx) <= t) GO TO 50 ind(ij) = indx ind(i) = it it = indx t = x(it) GO TO 50 ! INTERCHANGE ELEMENTS K AND L 40 itt = ind(l) ind(l) = ind(k) ind(k) = itt ! FIND AN ELEMENT IN THE UPPER PART OF THE ARRAY WHICH IS ! NOT LARGER THAN T 50 l = l - 1 indx = ind(l) IF (x(indx) > t) GO TO 50 ! FIND AN ELEMENT IN THE LOWER PART OF THE ARRAY WHCIH IS NOT SMALLER THAN T 60 k = k + 1 indx = ind(k) IF (x(indx) < t) GO TO 60 ! IF K <= L, INTERCHANGE ELEMENTS K AND L IF (k <= l) GO TO 40 ! SAVE THE UPPER AND LOWER SUBSCRIPTS OF THE PORTION OF THE ! ARRAY YET TO BE SORTED IF (l-i > j-k) THEN il(m) = i iu(m) = l i = k m = m + 1 GO TO 80 END IF il(m) = k iu(m) = j j = l m = m + 1 GO TO 80 ! BEGIN AGAIN ON ANOTHER UNSORTED PORTION OF THE ARRAY 70 m = m - 1 IF (m == 0) RETURN i = il(m) j = iu(m) 80 IF (j-i >= 11) GO TO 30 IF (i == 1) GO TO 20 i = i - 1 ! SORT ELEMENTS I+1,...,J. NOTE THAT 1 <= I < J AND J-I < 11. 90 i = i + 1 IF (i == j) GO TO 70 indx = ind(i+1) t = x(indx) it = indx indx = ind(i) IF (x(indx) <= t) GO TO 90 k = i 100 ind(k+1) = ind(k) k = k - 1 indx = ind(k) IF (t < x(indx)) GO TO 100 ind(k+1) = it GO TO 90 END SUBROUTINE qsortd subroutine linint(ndim,x,y,n,x5,y5) IMPLICIT NONE INTEGER,INTENT(IN) :: ndim ! number of grid points REAL,INTENT(IN) :: x(ndim) ! coordinate REAL,INTENT(IN) :: y(ndim) ! variable INTEGER,INTENT(IN) :: n ! number of targets REAL,INTENT(IN) :: x5(n) ! target coordinates REAL,INTENT(OUT) :: y5(n) ! target values integer :: i,j,k do j=1,n do i=2,ndim if(x5(j)==x(i)) then y5(j)=y(i) else if(x5(j)>x(i-1) .and. x5(j)x(ndim) .or. x5(j) r2) CYCLE ! Reject P if outside acceptance region IF (v**2 < -4.0*LOG(u)*u**2) EXIT END DO ! Return ratio of P's coordinates as the normal deviate fn_val = v/u RETURN END FUNCTION random_normal subroutine range_remap() !def remap( x, oMin, oMax, nMin, nMax ): ! #range check ! if oMin == oMax: ! print "Warning: Zero input range" ! return None ! if nMin == nMax: ! print "Warning: Zero output range" ! return None ! #check reversed input range ! reverseInput = False ! oldMin = min( oMin, oMax ) ! oldMax = max( oMin, oMax ) ! if not oldMin == oMin: ! reverseInput = True ! #check reversed output range ! reverseOutput = False ! newMin = min( nMin, nMax ) ! newMax = max( nMin, nMax ) ! if not newMin == nMin : ! reverseOutput = True ! portion = (x-oldMin)*(newMax-newMin)/(oldMax-oldMin) ! if reverseInput: ! portion = (oldMax-x)*(newMax-newMin)/(oldMax-oldMin) ! result = portion + newMin ! if reverseOutput: ! result = newMax - portion ! return result end subroutine SUBROUTINE exsearch(PLON,PLAT,IDM,JDM,XP,YP,IPr,JPr) IMPLICIT NONE !C INTEGER IDM,JDM REAL XP,YP REAL PLON(IDM,JDM),PLAT(IDM,JDM) !C !C MOST OF WORK IS DONE HERE. !C REAL DX,DY,DEG2RAD,DIST,DISTJ,DIST_MAX,D180,D360,QDX REAL PLAT_MIN(JDM),PLAT_MAX(JDM) LOGICAL LCYCLE,ldebug INTEGER I,IP,J,JP real ipr,jpr !C !C !C FIND THE NEAREST POINT BY EXHAUSTIVE SEARCH (NOT EFFICIENT) !C OPTIMIZE FOR MULTIPLE POINT CASE. !C ldebug = .false. D180 = 180.0 D360 = 360.0 DEG2RAD = 4.D0*ATAN(1.D0)/180.D0 !PI/180 !C DO J= 1,JDM PLAT_MIN(J) = MINVAL(PLAT(:,J)) PLAT_MAX(J) = MAXVAL(PLAT(:,J)) ENDDO dist_max = 0.0 do j= 1,jdm-1 do i= 1,idm-1 dist_max = max( abs(plat(i,j)-plat(i+1,j)), abs(plat(i,j)-plat(i,j+1)),dist_max ) enddo enddo ! target must be at leastthis close in latitude dist_max = 2*dist_max !C IP = 1 JP = 1 DO !input loop QDX = MAX(0.001,ABS(COS(YP*DEG2RAD))) DY = ABS(PLAT(IP,JP) - YP) DX = MOD( ABS(PLON(IP,JP) - XP), D360 ) IF (DX.GT.D180) THEN DX = D360 - DX ENDIF DIST = QDX*DX+DY lcycle = .false. DO J= 1,JDM distj = min(dist,dist_max) if (.not. ldebug) then if (yp.lt.plat_min(j)-distj .or. yp.gt.plat_max(j)+distj ) then cycle ! far away row endif else !debug if (yp.lt.plat_min(j)-distj .or. yp.gt.plat_max(j)+distj ) then if (.not.lcycle) then write(6,'(a,5x,i5,f9.2)') 'j,dist (cycle-strt)',j,dist elseif (j.eq.jdm) then write(6,'(a,5x,i5,f9.2)') 'j,dist (cycle-stop)',j,dist endif lcycle = .true. cycle ! far away row else if (lcycle) then write(6,'(a,5x,i5,f9.2)') 'j,dist (cycle-stop)',j-1,dist endif lcycle = .false. endif endif !.not.ldebug;else if (dist.eq.0.0) then exit ! found exact location endif DO I= 1,IDM DY = ABS(PLAT(I,J) - YP) DX = MOD( ABS(PLON(I,J) - XP), D360 ) IF (DX.GT.D180) THEN DX = D360 - DX ENDIF IF (QDX*DX+DY.LT.DIST) THEN IP = I JP = J DIST = QDX*DX+DY if (ldebug) then write(6,'(a,2i5,3f9.2)') 'ip,jp,dx,dy,dist = ', ip,jp,dx,dy,dist endif ENDIF ENDDO ENDDO !C WRITE(6,'(I5,I5)') IP,JP !C !C IF (.NOT.LSINGLE) THEN !C READ(5,*,IOSTAT=IOS) XP,YP !C IF (IOS.NE.0) THEN !C EXIT !C ENDIF !C ELSE EXIT !no input !C ENDIF ENDDO !input loop ipr=real(ip);jpr=real(jp) RETURN END SUBROUTINE subroutine layer2z(a,p,az,z,flag,ii,jj,kk,kz,itype) implicit none !c !c********** !c* !c 1) interpolate a layered field to fixed z depths !c !c 2) input arguments: !c a - scalar field in layer space !c p - layer interface depths (non-negative m) !c p(:,:, 1) is the surface !c p(:,:,kk+1) is the bathymetry !c z - target z-level depths (non-negative m) !c flag - data void (land) marker !c ii - 1st dimension of a,p,az !c jj - 2nd dimension of a,p,az !c kk - 3rd dimension of a (number of layers) !c kz - 3rd dimension of az (number of levels) !c itype - interpolation type !c =-2; piecewise quadratic across each layer !c =-1; piecewise linear across each layer !c = 0; sample the layer spaning each depth !c = 1; linear interpolation between layer centers !c = 2; linear interpolation between layer interfaces !c !c 3) output arguments: !c az - scalar field in z-space !c !c 4) except at data voids, must have: !c p(:,:, 1) == zero (surface) !c p(:,:, l+1) >= p(:,:,l) !c p(:,:,kk+1) == bathymetry !c 0 <= z(k) <= z(k+1) !c note that z(k) > p(i,j,kk+1) implies that az(i,j,k)=flag, !c since the z-level is then below the bathymetry. !c !c 5) Alan J. Wallcraft, Naval Research Laboratory, February 2002. !c* !c********** integer(i4) :: ii,jj,kk,kz,itype,i,j,k,l,lf real (r4) :: a(ii,jj,kk),p(ii,jj,kk+1),az(ii,jj,kz),z(kz),flag real(r4) :: s,zk,z0,zm,zp,si(kk,1),pi(kk+1),so(kz,1) do j= 1,jj do i= 1,ii if (a(i,j,1).eq.flag) then do k= 1,kz az(i,j,k) = flag ! land enddo elseif (itype.lt.0) then do k= 1,kk si(k,1) = a(i,j,k) pi(k) = p(i,j,k) enddo pi(kk+1) = p(i,j,kk+1) if (itype.eq.-1) then ! call layer2z_plm(si,pi,kk,1,so,z,kz, flag) else ! call layer2z_ppm(si,pi,kk,1,so,z,kz, flag) endif do k= 1,kz az(i,j,k) = so(k,1) enddo else !itype.ge.0 lf=1 do k= 1,kz zk=z(k) do l= lf,kk if (p(i,j,l).le.zk .and. p(i,j,l+1).ge.zk) then !c !c z(k) is in layer l. !c if (itype.eq.0) then !c !c sample the layer !c az(i,j,k) = a(i,j,l) elseif (itype.eq.1) then elseif (itype.eq.1) then !1c !c linear interpolation between layer centers z0 = 0.5*(p(i,j,l)+p(i,j,l+1)) if (zk.le.z0) then !c !c z(k) is in the upper half of the layer !c if (l.eq.1) then az(i,j,k) = a(i,j,1) else zm = 0.5*(p(i,j,l-1)+p(i,j,l)) s = (z0 - zk)/(z0 - zm) az(i,j,k) = s*a(i,j,l-1) + (1.0-s)*a(i,j,l) endif else !c !c z(k) is in the lower half of the layer !c if (p(i,j,l+1).eq.p(i,j,kk+1)) then az(i,j,k) = a(i,j,kk) else zp = 0.5*(p(i,j,l+1)+p(i,j,l+2)) s = (zk - z0)/(zp - z0) az(i,j,k) = s*a(i,j,l+1) + (1.0-s)*a(i,j,l) endif endif else !itype.eq.2 !c !c !c linear interpolation between layer interfaces !c a is a vertical integral from the surface !c if (l.eq.1) then s = zk/p(i,j,2) az(i,j,k) = s*a(i,j,1) else s = (zk - p(i,j,l))/(p(i,j,l+1) - p(i,j,l)) az(i,j,k) = (1.0-s)*a(i,j,l-1) + s*a(i,j,l) endif endif lf = l exit elseif (l.eq.kk) then az(i,j,k) = flag ! below the bottom lf = l exit endif enddo !l enddo !k endif !land:itype<0:else enddo !i enddo !j return end subroutine subroutine ppm(pt, s,sc,ki,ks) implicit none !c********** !c* !c 1) generate a monotonic PPM interpolation of a layered field: !c Colella, P. & P.R. Woodward, 1984, J. Comp. Phys., 54, 174-201. !c !c 2) input arguments: !c pt - layer interface thicknesses (non-zero) !c s - scalar fields in layer space !c ki - 1st dimension of s (number of layers) !c ks - 2nd dimension of s (number of fields) !c !c 3) output arguments: !c sc - scalar field coefficients for PPM interpolation !c !c 4) except at data voids, must have: !c pi( 1) == zero (surface) !c pi( l+1) >= pi(:,:,l) !c pi(ki+1) == bathymetry !c !c 5) Tim Campbell, Mississippi State University, September 2002; !C Alan J. Wallcraft, Naval Research Laboratory, August 2007. !c* !c********** !c integer(i4) :: ki,ks,j,k,l real(r4) :: pt(ki+1),s(ki,ks),sc(ki,ks,3) real(r4) :: da,a6,slj,scj,srj real(r4) :: as(ki),al(ki),ar(ki) real(r4) :: ptjp(ki), pt2jp(ki), ptj2p(ki),qptjp(ki),qpt2jp(ki),qptj2p(ki),ptq3(ki),qpt4(ki) !compute grid metrics do j=1,ki ptq3( j) = pt(j)/(pt(j-1)+pt(j)+pt(j+1)) ptjp( j) = pt(j) + pt(j+1) pt2jp(j) = pt(j) + ptjp(j) ptj2p(j) = ptjp(j) + pt(j+1) qptjp( j) = 1.0/ptjp( j) qpt2jp(j) = 1.0/pt2jp(j) qptj2p(j) = 1.0/ptj2p(j) enddo !j ptq3(2) = pt(2)/(pt(1)+ptjp(2)) do j=3,ki-1 ptq3(j) = pt(j)/(pt(j-1)+ptjp(j)) !pt(j)/ (pt(j-1)+pt(j)+pt(j+1)) qpt4(j) = 1.0/(ptjp(j-2)+ptjp(j)) !1.0/(pt(j-2)+pt(j-1)+pt(j)+pt(j+1)) enddo !j do l= 1,ks !Compute average slopes: Colella, Eq. (1.8) as(1)=0. do j=2,ki-1 slj=s(j, l)-s(j-1,l) srj=s(j+1,l)-s(j, l) if (slj*srj.gt.0.) then scj=ptq3(j)*( pt2jp(j-1)*srj*qptjp(j)+ptj2p(j) *slj*qptjp(j-1) ) as(j)=sign(min(abs(2.0*slj),abs(scj),abs(2.0*srj)),scj) else as(j)=0. endif enddo !j as(ki)=0. !Compute "first guess" edge values: Colella, Eq. (1.6) al(1)=s(1,l) !1st layer PCM ar(1)=s(1,l) !1st layer PCM al(2)=s(1,l) !1st layer PCM do j=3,ki-1 al(j)=s(j-1,l)+pt(j-1)*(s(j,l)-s(j-1,l))*qptjp(j-1) & +qpt4(j)*( & 2.*pt(j)*pt(j-1)*qptjp(j-1)*(s(j,l)-s(j-1,l))* & ( ptjp(j-2)*qpt2jp(j-1) & -ptjp(j) *qptj2p(j-1) ) & -pt(j-1)*as(j) *ptjp(j-2)*qpt2jp(j-1) & +pt(j) *as(j-1)*ptjp(j) *qptj2p(j-1) ) ar(j-1)=al(j) enddo !j ar(ki-1)=s(ki,l) !last layer PCM al(ki) =s(ki,l) !last layer PCM ar(ki) =s(ki,l) !last layer PCM !Impose monotonicity: Colella, Eq. (1.10) do j=2,ki-1 if ((s(j+1,l)-s(j,l))*(s(j,l)-s(j-1,l)).le.0.) then !local extremum al(j)=s(j,l) ar(j)=s(j,l) else da=ar(j)-al(j) a6=6.0*s(j,l)-3.0*(al(j)+ar(j)) if (da*a6 .gt. da*da) then !peak in right half of zone al(j)=3.0*s(j,l)-2.0*ar(j) elseif (da*a6 .lt. -da*da) then !peak in left half of zone ar(j)=3.0*s(j,l)-2.0*al(j) endif endif enddo !j !Set coefficients do j=1,ki sc(j,l,1)=al(j) sc(j,l,2)=ar(j)-al(j) sc(j,l,3)=6.0*s(j,l)-3.0*(al(j)+ar(j)) enddo !j enddo !l return end subroutine ppm SUBROUTINE locpt (x0, y0, x, y, n, l, m) !----------------------------------------------------------------------- ! GIVEN A POLYGONAL LINE CONNECTING THE VERTICES (X(I),Y(I)) (I = 1,...,N) ! TAKEN IN THIS ORDER. IT IS ASSUMED THAT THE POLYGONAL PATH IS A LOOP, ! WHERE (X(N),Y(N)) = (X(1),Y(1)) OR THERE IS AN ARC FROM (X(N),Y(N)) TO ! (X(1),Y(1)). N.B. The polygon may cross itself any number of times. ! (X0,Y0) IS AN ARBITRARY POINT AND L AND M ARE VARIABLES. ! On output, L AND M ARE ASSIGNED THE FOLLOWING VALUES ... ! L = -1 IF (X0,Y0) IS OUTSIDE THE POLYGONAL PATH ! L = 0 IF (X0,Y0) LIES ON THE POLYGONAL PATH ! L = 1 IF (X0,Y0) IS INSIDE THE POLYGONAL PATH ! M = 0 IF (X0,Y0) IS ON OR OUTSIDE THE PATH. IF (X0,Y0) IS INSIDE THE ! PATH THEN M IS THE WINDING NUMBER OF THE PATH AROUND THE POINT (X0,Y0). ! Fortran 66 version by A.H. Morris ! Converted to ELF90 compatibility by Alan Miller, 15 February 1997 !----------------------- IMPLICIT NONE REAL, INTENT(IN) :: x0, y0, x(n), y(n) INTEGER, INTENT(IN) :: n INTEGER, INTENT(OUT) :: l, m ! Local variables INTEGER :: i, n0 REAL :: angle, eps, pi, pi2, sum, theta, theta1, thetai, tol, u, v ! ****** EPS IS A MACHINE DEPENDENT CONSTANT. EPS IS THE ! SMALLEST NUMBER SUCH THAT 1.0 + EPS > 1.0 eps = EPSILON(1.0) !----------------------------------------------------------------------- n0 = n IF (x(1) == x(n) .AND. y(1) == y(n)) n0 = n - 1 pi = ATAN2(0.0, -1.0) pi2 = 2.0*pi tol = 4.0*eps*pi l = -1 m = 0 u = x(1) - x0 v = y(1) - y0 IF (u == 0.0 .AND. v == 0.0) GO TO 20 IF (n0 < 2) RETURN theta1 = ATAN2(v, u) sum = 0.0 theta = theta1 DO i = 2, n0 u = x(i) - x0 v = y(i) - y0 IF (u == 0.0 .AND. v == 0.0) GO TO 20 thetai = ATAN2(v, u) angle = ABS(thetai - theta) IF (ABS(angle - pi) < tol) GO TO 20 IF (angle > pi) angle = angle - pi2 IF (theta > thetai) angle = -angle sum = sum + angle theta = thetai END DO angle = ABS(theta1 - theta) IF (ABS(angle - pi) < tol) GO TO 20 IF (angle > pi) angle = angle - pi2 IF (theta > theta1) angle = -angle sum = sum + angle ! SUM = 2*PI*M WHERE M IS THE WINDING NUMBER m = ABS(sum)/pi2 + 0.2 IF (m == 0) RETURN l = 1 IF (sum < 0.0) m = -m RETURN ! (X0, Y0) IS ON THE BOUNDARY OF THE PATH 20 l = 0 RETURN END SUBROUTINE locpt END MODULE