[Dart-dev] [4065] DART/trunk/models/POP/model_mod.f90: New algorithm for interpolation . Works wit
nancy at ucar.edu
nancy at ucar.edu
Mon Sep 28 16:26:19 MDT 2009
An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/dart-dev/attachments/20090928/1a88788a/attachment-0001.html
-------------- next part --------------
Modified: DART/trunk/models/POP/model_mod.f90
===================================================================
--- DART/trunk/models/POP/model_mod.f90 2009-09-28 20:10:01 UTC (rev 4064)
+++ DART/trunk/models/POP/model_mod.f90 2009-09-28 22:26:19 UTC (rev 4065)
@@ -294,6 +294,7 @@
call read_vert_grid( Nz, ZC, ZG)
if (debug > 0) call write_grid_netcdf() ! DEBUG only
+if (debug > 0) call write_grid_interptest() ! DEBUG only
!---------------------------------------------------------------
! compute the offsets into the state vector for the start of each
@@ -313,10 +314,10 @@
! e.g. S,T,U,V = 256 x 225 x 70
! e.g. SHGT = 256 x 225
-if (do_output()) write(logfileunit, *) 'Using grid size : '
-if (do_output()) write(logfileunit, *) ' Nx, Ny, Nz = ', Nx, Ny, Nz
-if (do_output()) write( * , *) 'Using grid size : '
-if (do_output()) write( * , *) ' Nx, Ny, Nz = ', Nx, Ny, Nz
+if (do_output()) write(logfileunit, *) 'Using grid : Nx, Ny, Nz = ', &
+ Nx, Ny, Nz
+if (do_output()) write( * , *) 'Using grid : Nx, Ny, Nz = ', &
+ Nx, Ny, Nz
model_size = (n3dfields * (Nx * Ny * Nz)) + (n2dfields * (Nx * Ny))
if (do_output()) write(*,*) 'model_size = ', model_size
@@ -1306,20 +1307,20 @@
! Return in_quad true if the point (lon, lat) is in the quad with
! the given corners.
-! Do this by ray tracing in latitude for now. For non-pole point, want a
-! downward directed ray from (lon,lat) to intersect only one edge of the
-! quadrilateral. Have to watch out for intersecting at a vertex nearly
-! exactly, other round-off issues. This should be checked carefully.
+! Do this by line tracing in latitude for now. For non-pole point, want a vertical
+! line from the lon, lat point to intersect a side of the quad both above
+! and below the point.
real(r8) :: x(2), y(2)
-logical :: cant_be_in_box, on_side
-integer :: intercepts(4), exact_corner(4), second_corner, num_sides, num_corners, i
+logical :: cant_be_in_box, in_box
+integer :: intercepts_above(4), intercepts_below(4), i
+integer :: num_above, num_below
! Default answer is point is not in quad
in_quad = .false.
-! Loop through the sides and compute intercept (if any) with a downward
-! ray from the point. This ray has equation x=lon, y<=lat.
+! Loop through the sides and compute intercept (if any) with a vertical line
+! from the point. This line has equation x=lon.
do i = 1, 4
! Load up the sides endpoints
if(i <= 3) then
@@ -1332,273 +1333,123 @@
y(2) = y_corners(1)
endif
- ! Check to see how a southward going ray from the point is related to this side
- call ray_intercept(x, y, lon, lat, cant_be_in_box, intercepts(i), &
- exact_corner(i), on_side)
+ ! Check to see how a vertical line from the point is related to this side
+ call line_intercept(x, y, lon, lat, cant_be_in_box, in_box, intercepts_above(i), &
+ intercepts_below(i))
! If cant_be_in_box is true, can return right away
if(cant_be_in_box) then
in_quad = .false.
return
! Return true if on a side
- else if(on_side) then
+ else if(in_box) then
in_quad = .true.
return
endif
enddo
-! Point is not known to be out of the quad or on the quads boundary.
-! See how many sides and corners it went through
-num_sides = sum(intercepts)
-num_corners = sum(exact_corner)
+! See if the line intercepted a side of the quad both above and below
+num_above = sum(intercepts_above)
+num_below = sum(intercepts_below)
-! If the number of exact corner intercepts is 0 just look at sides.
-! If the ray only intersects one side we are in the box.
-if(sum(exact_corner) == 0) then
- if(num_sides == 1) then
- ! Point is in box
- in_quad = .true.
- else if(num_sides >= 3) then
- ! Put in some development error checking; shouldn't be able to hit 3 or 4 sides
- write(msgstring, *) 'Num sides in in_quad is ', num_sides, &
- 'Please contact DART development team'
- call error_handler(E_ERR, 'in_quad', msgstring, source, revision, revdate)
- endif
-
-else
- ! We have some exact_corner hits; these must come in pairs.
- ! An odd number of corners should be impossible.
- if(num_corners == 1 .or. num_corners == 3) then
- ! Put in some development error checking; shouldn't be able to hit 3 or 4 sides
- write(*,*)'xcorners = ',x_corners,';'
- write(*,*)'ycorners = ',y_corners,';'
- write(*,*)'lon = ',lon,';'
- write(*,*)'lat = ',lat,';'
- write(*,*)'intercepts = ',intercepts,';'
- write(*,*)'exact_corner = ',exact_corner,';'
- write(msgstring, *) 'Num corners in in_quad is ', num_corners, &
- 'Please contact DART development team'
- call error_handler(E_MSG, 'in_quad', msgstring, source, revision, revdate)
- in_quad = .false. ! DEBUG FIXME TJH ... not the right thing to do.
- ! Four corners means we are not in;
- else if(num_corners == 4) then
- in_quad = .false.
- return
- ! Two corners and one side means we are not in;
- else if(num_corners == 2 .and. num_sides == 1) then
- in_quad = .true.
- ! Two corners and no plain intercepts means we could be in or out;
- ! This case requires examining the particular sides.
- else if(num_corners == 2 .and. num_sides == 0) then
- ! See if the point is between these; if so in box
- call two_corners(exact_corner, x_corners, lon, in_quad)
- endif
-
+if(num_above > 0 .and. num_below > 0) then
+ in_quad = .true.
endif
end function in_quad
-!------------------------------------------------------------
-subroutine two_corners(exact_corner, x_corners, lon, in_quad)
-
-integer, intent(in) :: exact_corner(4)
-real(r8), intent(in) :: x_corners(4), lon
-logical, intent(out) :: in_quad
-
-! For a ray that intersects the corner between two quad edges, is the
-! point in or out of the quad?
-
-integer :: ind, i, second_point
-real(r8) :: x_points(4)
-
-ind = 1
-
-! Loop to get a list (x_points) that contains the longitude of endpoints of
-! both of the sides that share a corner that was hit by the ray. There will
-! be one duplicate value.
-do i = 1, 4
- if(exact_corner(i) == 1) then
- ! Side i was one of the ones that got a corner hit
- x_points(ind) = x_corners(i)
- ind = ind + 1
- second_point = i + 1
- if(second_point > 4) second_point = 1
- x_points(ind) = x_corners(second_point)
- ind = ind + 1
- endif
-end do
-
-! Get the min and max of the longitudes on the edges adjacent to the corner.
-! If the point has longitude between these, then it must lie in the
-! angle between them and the point is in the quad.
-if(lon >= minval(x_points) .and. lon <= maxval(x_points)) then
- in_quad = .true.
-else
- in_quad = .false.
-endif
-
-end subroutine two_corners
-
!------------------------------------------------------------
-subroutine ray_intercept(x_in, y, x_ray_in, y_ray, &
- cant_be_in_box, intercepts, exact_corner, on_side)
+subroutine line_intercept(side_x_in, side_y, x_point_in, y_point, &
+ cant_be_in_box, in_box, intercept_above, intercept_below)
-real(r8), intent(in) :: x_in(2), y(2), x_ray_in, y_ray
-logical, intent(out) :: cant_be_in_box, on_side
-integer, intent(out) :: intercepts, exact_corner
+real(r8), intent(in) :: side_x_in(2), side_y(2), x_point_in, y_point
+logical, intent(out) :: cant_be_in_box, in_box
+integer, intent(out) :: intercept_above, intercept_below
-! Find the intercept of a downward ray from (x_ray, y_ray) and
-! a line segment with endpoints x, y.
-! For a given side have endpoints (x1, y1) and (x2, y2) so equation of
-! segment is y=y1 + m(x-x1) for y between y1 and y2.
-! Intersection of ray and line occurs at y = y1 + m(lon - x1); need this
-! y to be less than lat and between y1 and y2).
-! If the ray is colinear with the side but the point is not on the side, return
-! cant_be_in_box as true. If the ray intercepts the side, return intercepts
-! as 1 else 0. If the ray intercepts a corner of the side, return exact_corner as
-! true. If the point lies on the side, return on_side as true.
+! Find the intercept of a vertical line from point (x_point, y_point) and
+! a line segment with endpoints side_x and side_y.
+! For a given side have endpoints (side_x1, side_y1) and (side_x2, side_y2)
+! so equation of segment is y = side_y1 + m(x-side_x1) for y
+! between side_y1 and side_y2.
+! Intersection of vertical line and line containing side
+! occurs at y = side_y1 + m(x_point - side_x1); need this
+! y to be between side_y1 and side_y2.
+! If the vertical line is colinear with the side but the point is not on the side, return
+! cant_be_in_box as true. If the point is on the side, return in_box true.
+! If the intersection of the vertical line and the side occurs at a point above
+! the given point, return 1 for intercept_above. If the intersection occurs
+! below, return 1 for intercept_below. If the vertical line does not intersect
+! the segment, return false and 0 for all intent out arguments.
-! WARNING: ALMOST CERTAINLY A PROBLEM FOR THE POLE BOX!!! POLE BOX COULD
+! WARNING: CERTAINLY A PROBLEM FOR THE POLE BOX!!! POLE BOX COULD
! HAVE SIDES THAT ARE LONGER THAN 180. For now pole boxes are excluded.
! This can probably be made much cleaner and more efficient.
-real(r8) :: slope, y_intercept, x(2), x_ray
-real(r8) :: xmin, xmax, ymin, ymax
-real(r8) :: rdummytiny, rdummyepsilon, rdummyabs, rdummyfraction
-integer :: idummydigits, idummyprecision, idummyexponent
+real(r8) :: slope, y_intercept, side_x(2), x_point
-xmin = minval(x_in)
-xmax = maxval(x_in)
-
! May have to adjust the longitude intent in values, so copy
-x = x_in
-x_ray = x_ray_in
+side_x = side_x_in
+x_point = x_point_in
! See if the side wraps around in longitude
-if(xmax - xmin > 180.0_r8) then
- if( x(1) < 180.0_r8) x(1) = x(1) + 360.0_r8
- if( x(2) < 180.0_r8) x(2) = x(2) + 360.0_r8
- if(x_ray < 180.0_r8) x_ray = x_ray + 360.0_r8
+if(maxval(side_x) - minval(side_x) > 180.0_r8) then
+ if(side_x(1) < 180.0_r8) side_x(1) = side_x(1) + 360.0_r8
+ if(side_x(2) < 180.0_r8) side_x(2) = side_x(2) + 360.0_r8
+ if(x_point < 180.0_r8) x_point = x_point + 360.0_r8
endif
-xmin = minval(x)
-xmax = maxval(x)
-ymin = minval(y)
-ymax = maxval(y)
+! Initialize the default returns
+cant_be_in_box = .false.
+in_box = .false.
+intercept_above = 0
+intercept_below = 0
-! Initialize all the possible returns
-cant_be_in_box = .false.
-on_side = .false.
-intercepts = 0
-exact_corner = 0
+! First easy check, if x_point is not between endpoints of segment doesn't intersect
+if(x_point < minval(side_x) .or. x_point > maxval(side_x)) return
-! First easy check, if x_ray is not between x(1) and x(2) it can't intersect
-if(x_ray < minval(x) .or. x_ray > maxval(x)) return
+! Otherwise line must intersect the segment
-if((x_ray < xmin) .or. (x_ray > xmax)) then
- write(*,*)'minval vs xmin makes a difference'
- return
-endif
-
-!if(x(2) == x(1)) then
-! continue
-!else
-! if (x(2) - x(1) == 0.0) &
-! write(*,*)'x difference is zero', x
-! if (abs(x(2) -x(1)) < epsilon(x(1))) &
-! write(*,*)'x is almost equal', x
-!endif
-!
-!if(y(2) == y(1)) then
-! continue
-!else
-! if (y(2) - y(1) == 0.0) then
-! write(*,*)'y difference is zero', y
-! endif
-! if (abs(y(2) -y(1)) < epsilon(y(1))) then
-! write(*,*)'y is almost equal', y
-! endif
-! if (abs(y(2) -y(1)) < tiny(y(1))) then
-! write(*,*)'y is almost equal', y
-! endif
-! rdummyabs = abs(y(2) - y(1))
-! rdummytiny = tiny(y(1))
-! rdummyepsilon = epsilon(y(1))
-! idummydigits = digits(y(1))
-! idummyprecision = precision(y(1))
-! idummyexponent = exponent(rdummyabs)
-! rdummyfraction = fraction(rdummyabs)
-!endif
-
-rdummyabs = abs(x(2) - x(1))
! First subblock, slope is undefined
-! if(x(2) == x(1)) then
-if(rdummyabs < 1.0E-12_r8 ) then ! near enough equal
- ! Check for exactly on this side
- if(x_ray /= x(1)) then
- ! Doesn't intersect vertical line
+if(side_x(2) == side_x(1)) then
+ ! The line is colinear with the side
+ ! If y_point is between endpoints then point is on this side
+ if(y_point <= maxval(side_y) .and. y_point >= minval(side_y)) then
+ in_box = .true.
return
+ ! If not on side but colinear with side, point cant be in quad
else
- ! The ray is colinear with the side
- ! If y_ray is between endpoints then point is on this side
- if(y_ray <= ymax .and. y_ray >= ymin) then
- on_side = .true.
- return
- ! If not on side but colinear with side, point cant be in quad
- else
- cant_be_in_box = .true.
- return
- endif
+ cant_be_in_box = .true.
+ return
endif
else
- rdummyabs = abs(y(2) - y(1))
- if ( rdummyabs < 1.0E-12_r8 ) then
- slope = 0.0_r8
- else
- ! Second possibility; slope is defined
- slope = (y(2) - y(1)) / (x(2) - x(1))
- endif
- ! Intercept of downward ray and line through points is at x_ray and...
- y_intercept = y(1) + slope * (x_ray - x(1))
+ ! Second possibility; slope is defined
+ ! FIXME: watch out for numerical instability.
+ ! near-zero x's and large side_y's may cause overflow
+ slope = (side_y(2) - side_y(1)) / (side_x(2) - side_x(1))
- ! If y_intercept is greater than y_ray it's not on the downward ray
- if(y_intercept > y_ray) return
+ ! Intercept of vertical line through is at x_point and...
+ y_intercept = side_y(1) + slope * (x_point - side_x(1))
- ! If the slope is 0, we now know the ray intercepts
- if(slope == 0) then
- ! Check for exact corner
- if(x_ray == x(1) .or. x_ray == x(2)) then
- exact_corner = 1
- return
- else
- intercepts = 1
- return
- endif
- endif
-
- ! Slope of line is not 0 or undefined so y endpoints differ
- ! If intercept is on segment and is below the point
- if(y_intercept > ymax .or. y_intercept < ymin) then
- ! Intersects line containing side outside of side
+ ! Intersects the segment, is it above, below, or at the point
+ if(y_intercept == y_point) then
+ in_box = .true.
return
- else if(y_intercept < ymax .and. y_intercept > ymin) then
- ! Intercepts inside the side
- intercepts = 1
+ else if(y_intercept > y_point) then
+ intercept_above = 1
return
else
- ! Otherwise, must go exactly through a corner
- exact_corner = 1
+ intercept_below = 1
return
endif
endif
-end subroutine ray_intercept
+end subroutine line_intercept
!------------------------------------------------------------
@@ -1779,6 +1630,7 @@
fract = 1.0_r8
if (debug > 1) print *, 'above first level in height'
if (debug > 1) print *, 'hgt_array, top, bot, fract=', hgt_array(1), top, bot, fract
+ return
endif
! Search through the boxes
@@ -3043,6 +2895,8 @@
if ( .not. module_initialized ) call static_init_model
+is_not_ocean = .FALSE. ! start out thinking everything is wet.
+
is_ugrid = is_on_ugrid(obs_type)
if (( is_ugrid .and. hgt_index > KMU(lon_index, lat_index)) .or. &
(.not. is_ugrid .and. hgt_index > KMT(lon_index, lat_index))) then
@@ -3074,81 +2928,167 @@
!
! Write the grid to a netcdf file for checking.
- integer :: ncid, NlonDimID, NlatDimID, NzDimID
- integer :: nlon, nlat, nz
- integer :: ulatVarID, ulonVarID, TLATvarid, TLONvarid
- integer :: ZGvarid, ZCvarid, KMTvarid, KMUvarid
+integer :: ncid, NlonDimID, NlatDimID, NzDimID
+integer :: nlon, nlat, nz
+integer :: ulatVarID, ulonVarID, TLATvarid, TLONvarid
+integer :: ZGvarid, ZCvarid, KMTvarid, KMUvarid
- integer :: dimids(2);
+integer :: dimids(2);
- if ( .not. module_initialized ) call static_init_model
+if ( .not. module_initialized ) call static_init_model
- nlon = size(ULAT,1)
- nlat = size(ULAT,2)
- nz = size(ZG)
-
- call nc_check(nf90_create('dart_grid.nc', NF90_CLOBBER, ncid),'write_grid_netcdf')
+nlon = size(ULAT,1)
+nlat = size(ULAT,2)
+nz = size(ZG)
- ! define dimensions
+call nc_check(nf90_create('dart_grid.nc', NF90_CLOBBER, ncid),'write_grid_netcdf')
- call nc_check(nf90_def_dim(ncid, 'i', nlon, NlonDimID),'write_grid_netcdf')
- call nc_check(nf90_def_dim(ncid, 'j', nlat, NlatDimID),'write_grid_netcdf')
- call nc_check(nf90_def_dim(ncid, 'k', nz, NzDimID),'write_grid_netcdf')
+! define dimensions
- dimids(1) = NlonDimID
- dimids(2) = NlatDimID
+call nc_check(nf90_def_dim(ncid, 'i', nlon, NlonDimID),'write_grid_netcdf')
+call nc_check(nf90_def_dim(ncid, 'j', nlat, NlatDimID),'write_grid_netcdf')
+call nc_check(nf90_def_dim(ncid, 'k', nz, NzDimID),'write_grid_netcdf')
- ! define variables
+dimids(1) = NlonDimID
+dimids(2) = NlatDimID
- ! FIXME: we should add attributes to say what units the grids are in (degrees).
- call nc_check(nf90_def_var(ncid, 'KMT', nf90_int, dimids, KMTvarid),'write_grid_netcdf')
- call nc_check(nf90_def_var(ncid, 'KMU', nf90_int, dimids, KMUvarid),'write_grid_netcdf')
- call nc_check(nf90_def_var(ncid, 'ULON', nf90_double, dimids, ulonVarID),'write_grid_netcdf')
- call nc_check(nf90_def_var(ncid, 'ULAT', nf90_double, dimids, ulatVarID),'write_grid_netcdf')
- call nc_check(nf90_def_var(ncid, 'TLON', nf90_double, dimids, TLONvarid),'write_grid_netcdf')
- call nc_check(nf90_def_var(ncid, 'TLAT', nf90_double, dimids, TLATvarid),'write_grid_netcdf')
- call nc_check(nf90_def_var(ncid, 'ZG', nf90_double, NzDimID, ZGvarid),'write_grid_netcdf')
- call nc_check(nf90_def_var(ncid, 'ZC', nf90_double, NzDimID, ZCvarid),'write_grid_netcdf')
+! define variables
- call nc_check(nf90_put_att(ncid,ulonVarID,'long_name','U,V grid lons'), &
- 'write_grid_netcdf')
- call nc_check(nf90_put_att(ncid,ulatVarID,'long_name','U,V grid lats'), &
- 'write_grid_netcdf')
- call nc_check(nf90_put_att(ncid,tlonVarID,'long_name','S,T grid lons'), &
- 'write_grid_netcdf')
- call nc_check(nf90_put_att(ncid,tlatVarID,'long_name','S,T grid lats'), &
- 'write_grid_netcdf')
- call nc_check(nf90_put_att(ncid,ZCVarID,'long_name','vertical grid centers'), &
- 'write_grid_netcdf')
- call nc_check(nf90_put_att(ncid,ZCVarID,'units','meters'), &
- 'write_grid_netcdf')
- call nc_check(nf90_put_att(ncid,ZGVarID,'long_name','vertical grid bottoms'), &
- 'write_grid_netcdf')
- call nc_check(nf90_put_att(ncid,ZGVarID,'units','meters'), &
- 'write_grid_netcdf')
+! FIXME: we should add attributes to say what units the grids are in (degrees).
+call nc_check(nf90_def_var(ncid, 'KMT', nf90_int, dimids, KMTvarid),'write_grid_netcdf')
+call nc_check(nf90_def_var(ncid, 'KMU', nf90_int, dimids, KMUvarid),'write_grid_netcdf')
+call nc_check(nf90_def_var(ncid, 'ULON', nf90_double, dimids, ulonVarID),'write_grid_netcdf')
+call nc_check(nf90_def_var(ncid, 'ULAT', nf90_double, dimids, ulatVarID),'write_grid_netcdf')
+call nc_check(nf90_def_var(ncid, 'TLON', nf90_double, dimids, TLONvarid),'write_grid_netcdf')
+call nc_check(nf90_def_var(ncid, 'TLAT', nf90_double, dimids, TLATvarid),'write_grid_netcdf')
+call nc_check(nf90_def_var(ncid, 'ZG', nf90_double, NzDimID, ZGvarid),'write_grid_netcdf')
+call nc_check(nf90_def_var(ncid, 'ZC', nf90_double, NzDimID, ZCvarid),'write_grid_netcdf')
- call nc_check(nf90_enddef(ncid),'write_grid_netcdf')
+call nc_check(nf90_put_att(ncid,ulonVarID,'long_name','U,V grid lons'), &
+ 'write_grid_netcdf')
+call nc_check(nf90_put_att(ncid,ulatVarID,'long_name','U,V grid lats'), &
+ 'write_grid_netcdf')
+call nc_check(nf90_put_att(ncid,tlonVarID,'long_name','S,T grid lons'), &
+ 'write_grid_netcdf')
+call nc_check(nf90_put_att(ncid,tlatVarID,'long_name','S,T grid lats'), &
+ 'write_grid_netcdf')
+call nc_check(nf90_put_att(ncid,ZCVarID,'long_name','vertical grid centers'), &
+ 'write_grid_netcdf')
+call nc_check(nf90_put_att(ncid,ZCVarID,'units','meters'), &
+ 'write_grid_netcdf')
+call nc_check(nf90_put_att(ncid,ZGVarID,'long_name','vertical grid bottoms'), &
+ 'write_grid_netcdf')
+call nc_check(nf90_put_att(ncid,ZGVarID,'units','meters'), &
+ 'write_grid_netcdf')
- ! fill variables
+call nc_check(nf90_enddef(ncid),'write_grid_netcdf')
- call nc_check(nf90_put_var(ncid, KMTvarid, KMT),'write_grid_netcdf')
- call nc_check(nf90_put_var(ncid, KMUvarid, KMU),'write_grid_netcdf')
- call nc_check(nf90_put_var(ncid, ulatVarID, ULAT),'write_grid_netcdf')
- call nc_check(nf90_put_var(ncid, ulonVarID, ULON),'write_grid_netcdf')
- call nc_check(nf90_put_var(ncid, TLATvarid, TLAT),'write_grid_netcdf')
- call nc_check(nf90_put_var(ncid, TLONvarid, TLON),'write_grid_netcdf')
- call nc_check(nf90_put_var(ncid, ZGvarid, ZG),'write_grid_netcdf')
- call nc_check(nf90_put_var(ncid, ZCvarid, ZC),'write_grid_netcdf')
+! fill variables
- call nc_check(nf90_close(ncid),'write_grid_netcdf')
+call nc_check(nf90_put_var(ncid, KMTvarid, KMT),'write_grid_netcdf')
+call nc_check(nf90_put_var(ncid, KMUvarid, KMU),'write_grid_netcdf')
+call nc_check(nf90_put_var(ncid, ulatVarID, ULAT),'write_grid_netcdf')
+call nc_check(nf90_put_var(ncid, ulonVarID, ULON),'write_grid_netcdf')
+call nc_check(nf90_put_var(ncid, TLATvarid, TLAT),'write_grid_netcdf')
+call nc_check(nf90_put_var(ncid, TLONvarid, TLON),'write_grid_netcdf')
+call nc_check(nf90_put_var(ncid, ZGvarid, ZG),'write_grid_netcdf')
+call nc_check(nf90_put_var(ncid, ZCvarid, ZC),'write_grid_netcdf')
+call nc_check(nf90_close(ncid),'write_grid_netcdf')
+
end subroutine write_grid_netcdf
+
+ subroutine write_grid_interptest()
!------------------------------------------------------------------
+! Write the grid to an ascii file - in a format suitable for
+! subsequent use in the 'test_interpolation()' code.
+! write_grid_interptest is only possible after reading a real POP grid,
+! so static_init_model() must be called to gather the real POP grid.
-subroutine test_interpolation(test_casenum)
+integer :: i, j
+real(r8) :: rowmat(Nx,1), colmat(1,Ny), dmat(Nx,Ny)
+real(r8) :: rtlon, rulon, rtlat, rulat, u_val, t_val
+if ( .not. module_initialized ) call static_init_model
+
+!----------------------------------------------------------------------
+! Generate a 'Regular' grid with the same rough 'shape' as the dipole grid
+!----------------------------------------------------------------------
+
+open(unit=12, position='rewind', action='write', file='regular_grid_u')
+open(unit=13, position='rewind', action='write', file='regular_grid_t')
+open(unit=14, position='rewind', action='write', file='regular_grid_u_data')
+open(unit=15, position='rewind', action='write', file='regular_grid_t_data')
+
+write(12, *) Nx, Ny
+write(13, *) Nx, Ny
+
+! Have T-grid starting at 0 and U grid offset by half
+do i = 1, Nx
+ rtlon = (i - 1.0_r8) * 360.0_r8 / Nx
+ rulon = (i - 0.5_r8) * 360.0_r8 / Nx
+ do j = 1, Ny
+ rtlat = -90.0_r8 + (j - 1.0_r8) * 180.0_r8 / Ny
+ rulat = -90.0_r8 + (j - 0.5_r8) * 180.0_r8 / Ny
+ write(12, *) i, j, rulon, rulat
+ write(13, *) i, j, rtlon, rtlat
+ ! Now add some wave pattern data
+ u_val = sin(3.0_r8*(rulat + 11.0_r8)*2.0_r8*PI/360.0_r8) * &
+ sin(4.0_r8*(rulon + 17.0_r8)*2.0_r8*PI/360.0_r8)
+ t_val = sin(3.0_r8*(rtlat + 11.0_r8)*2.0_r8*PI/360.0_r8) * &
+ sin(4.0_r8*(rtlon + 17.0_r8)*2.0_r8*PI/360.0_r8)
+ write(14, *) rulon, rulat, u_val
+ write(15, *) rtlon, rtlat, t_val
+ enddo
+enddo
+
+close(unit=12)
+close(unit=13)
+close(unit=14)
+close(unit=15)
+
+!----------------------------------------------------------------------
+! POP grid (dipole) next
+!----------------------------------------------------------------------
+
+open(unit=12, position='rewind', action='write', file='dipole_grid_u')
+open(unit=13, position='rewind', action='write', file='dipole_grid_t')
+open(unit=14, position='rewind', action='write', file='dipole_grid_u_data')
+open(unit=15, position='rewind', action='write', file='dipole_grid_t_data')
+
+write(12, *) Nx, Ny
+write(13, *) Nx, Ny
+
+rowmat(:,1) = cos(PI * real((/ (i,i=0,Nx-1) /),r8) / Nx);
+colmat(1,:) = sin(PI * real((/ (i,i=0,Ny-1) /),r8) / Ny);
+dmat = matmul(rowmat,colmat)
+
+do i = 1, Nx
+ do j = 1, Ny
+ write(12, *) i, j, ULON(i,j), ULAT(i,j)
+ write(13, *) i, j, TLON(i,j), TLAT(i,j)
+ write(14, *) ULON(i,j), ULAT(i,j), dmat(i, j)
+ write(15, *) TLON(i,j), TLAT(i,j), dmat(i, j)
+ end do
+end do
+
+close(unit=12)
+close(unit=13)
+close(unit=14)
+close(unit=15)
+
+end subroutine write_grid_interptest
+
+
+
+ subroutine test_interpolation(test_casenum)
+!------------------------------------------------------------------
+!subroutine test_interpolation(test_casenum)
+!
+! rigorous test of the intepolation routine.
+
integer, intent(in) :: test_casenum
! This is storage just used for temporary test driver
@@ -3169,36 +3109,36 @@
! Set of block options for now
! Lon lat for u on channel 12, v on 13, u data on 14, t data on 15
-! Case 1: regular grid to dipole x3
-! Case 2: dipole x3 to regular grid
-! Case 3: regular grid to regular grid with same grid as x3 in SH
-! Case 4: regular grid with same grid as x3 in SH to regular grid
-! Case 5: regular grid with same grid as x3 in SH to dipole x3
-! Case 6: dipole x3 to regular grid with same grid as x3 in SH
+! Case 1: regular grid to dipole
+! Case 2: dipole to regular grid
+! Case 3: regular grid to regular grid with same grid as dipole in SH
+! Case 4: regular grid with same grid as dipole in SH to regular grid
+! Case 5: regular grid with same grid as dipole in SH to dipole
+! Case 6: dipole to regular grid with same grid as dipole in SH
! do not let the standard init run
module_initialized = .true.
if(test_casenum == 1 .or. test_casenum == 3) then
! Case 1 or 3: read in from regular grid
- open(unit = 12, file = 'regular_grid_u')
- open(unit = 13, file = 'regular_grid_t')
- open(unit = 14, file = 'regular_grid_u_data')
- open(unit = 15, file = 'regular_grid_t_data')
+ open(unit=12, position='rewind', action='read', file='regular_grid_u')
+ open(unit=13, position='rewind', action='read', file='regular_grid_t')
+ open(unit=14, position='rewind', action='read', file='regular_grid_u_data')
+ open(unit=15, position='rewind', action='read', file='regular_grid_t_data')
else if(test_casenum == 4 .or. test_casenum == 5) then
- ! Case 3 or 4: read regular with same grid as x3 in SH
- open(unit = 12, file = 'regular_griddi_u')
- open(unit = 13, file = 'regular_griddi_t')
- open(unit = 14, file = 'regular_griddi_u_data')
- open(unit = 15, file = 'regular_griddi_t_data')
+ ! Case 3 or 4: read regular with same grid as dipole in SH
+ open(unit=12, position='rewind', action='read', file='regular_griddi_u')
+ open(unit=13, position='rewind', action='read', file='regular_griddi_t')
+ open(unit=14, position='rewind', action='read', file='regular_griddi_u_data')
+ open(unit=15, position='rewind', action='read', file='regular_griddi_t_data')
else if(test_casenum == 2 .or. test_casenum == 6) then
- ! Case 5: read in from dipole x3 grid
- open(unit = 12, file = 'dipole_x3_grid_u')
- open(unit = 13, file = 'dipole_x3_grid_t')
- open(unit = 14, file = 'dipole_x3_u_data')
- open(unit = 15, file = 'dipole_x3_t_data')
+ ! Case 5: read in from dipole grid
+ open(unit=12, position='rewind', action='read', file='dipole_grid_u')
+ open(unit=13, position='rewind', action='read', file='dipole_grid_t')
+ open(unit=14, position='rewind', action='read', file='dipole_grid_u_data')
+ open(unit=15, position='rewind', action='read', file='dipole_grid_t_data')
endif
! Get the size of the grid from the input u and t files
@@ -3245,32 +3185,32 @@
call init_interp()
! Now read in the points for the output grid
-! Case 1: regular grid to dipole x3
-! Case 2: dipole x3 to regular grid
-! Case 3: regular grid to regular grid with same grid as x3 in SH
-! Case 4: regular grid with same grid as x3 in SH to regular grid
-! Case 5: regular grid with same grid as x3 in SH to dipole x3
-! Case 6: dipole x3 to regular grid with same grid as x3 in SH
+! Case 1: regular grid to dipole
+! Case 2: dipole to regular grid
+! Case 3: regular grid to regular grid with same grid as dipole in SH
+! Case 4: regular grid with same grid as dipole in SH to regular grid
+! Case 5: regular grid with same grid as dipole in SH to dipole
+! Case 6: dipole to regular grid with same grid as dipole in SH
if(test_casenum == 1 .or. test_casenum == 5) then
! Output to dipole grid
- open(unit = 22, file = 'dipole_x3_grid_u')
- open(unit = 23, file = 'dipole_x3_grid_t')
- open(unit = 24, file = 'dipole_x3_u_out')
- open(unit = 25, file = 'dipole_x3_t_out')
+ open(unit=22, position='rewind', action='read', file='dipole_grid_u')
+ open(unit=23, position='rewind', action='read', file='dipole_grid_t')
+ open(unit=24, position='rewind', action='write', file='dipole_grid_u_data.out')
+ open(unit=25, position='rewind', action='write', file='dipole_grid_t_data.out')
else if(test_casenum == 2 .or. test_casenum == 4) then
! Output to regular grid
- open(unit = 22, file = 'regular_grid_u')
- open(unit = 23, file = 'regular_grid_t')
- open(unit = 24, file = 'regular_grid_u_out')
- open(unit = 25, file = 'regular_grid_t_out')
+ open(unit=22, position='rewind', action='read', file='regular_grid_u')
+ open(unit=23, position='rewind', action='read', file='regular_grid_t')
+ open(unit=24, position='rewind', action='write', file='regular_grid_u_data.out')
+ open(unit=25, position='rewind', action='write', file='regular_grid_t_data.out')
else if(test_casenum == 3 .or. test_casenum == 6) then
- ! Output to regular grid with same grid as x3 in SH
- open(unit = 22, file = 'regular_griddi_u')
- open(unit = 23, file = 'regular_griddi_t')
- open(unit = 24, file = 'regular_griddi_u_out')
- open(unit = 25, file = 'regular_griddi_t_out')
+ ! Output to regular grid with same grid as dipole in SH
+ open(unit=22, position='rewind', action='read', file='regular_griddi_u')
+ open(unit=23, position='rewind', action='read', file='regular_griddi_t')
+ open(unit=24, position='rewind', action='write', file='regular_griddi_u_data.out')
+ open(unit=25, position='rewind', action='write', file='regular_griddi_t_data.out')
endif
read(22, *) dnx, dny
@@ -3282,17 +3222,43 @@
allocate(dulon(dnx, dny), dulat(dnx, dny), dtlon(dnx, dny), dtlat(dnx, dny))
allocate(dipole_u(dnx, dny), dipole_t(dnx, dny))
+
+dipole_u = 0.0_r8 ! just put some dummy values in to make sure they get changed.
+dipole_t = 0.0_r8 ! just put some dummy values in to make sure they get changed.
+
do imain = 1, dnx
+do jmain = 1, dny
+ read(22, *) ti, tj, dulon(imain, jmain), dulat(imain, jmain)
+ read(23, *) ti, tj, dtlon(imain, jmain), dtlat(imain, jmain)
+end do
+end do
+
+do imain = 1, dnx
do jmain = 1, dny
- read(22, *) ti, tj, dulon(imain, jmain), dulat(imain, jmain)
- read(23, *) ti, tj, dtlon(imain, jmain), dtlat(imain, jmain)
- ! Interpolate from the first grid to the second grid
+ ! Interpolate U from the first grid to the second grid
+
call lon_lat_interpolate(reg_u_x, dulon(imain, jmain), dulat(imain, jmain), &
KIND_U_CURRENT_COMPONENT, height, dipole_u(imain, jmain), istatus)
+ if ( istatus /= 0 ) then
+ write(msgstring,'(''cell'',i4,i4,1x,f12.8,1x,f12.8,'' U interp failed - code '',i4)') &
+ imain, jmain, dulon(imain, jmain), dulat(imain, jmain), istatus
+ call error_handler(E_MSG,'test_interpolation',msgstring,source,revision,revdate)
+ endif
+
write(24, *) dulon(imain, jmain), dulat(imain, jmain), dipole_u(imain, jmain)
+
+ ! Interpolate U from the first grid to the second grid
+
call lon_lat_interpolate(reg_t_x, dtlon(imain, jmain), dtlat(imain, jmain), &
KIND_TEMPERATURE, height, dipole_t(imain, jmain), istatus)
+
+ if ( istatus /= 0 ) then
+ write(msgstring,'(''cell'',i4,i4,1x,f12.8,1x,f12.8,'' T interp failed - code '',i4)') &
+ imain,jmain, dtlon(imain, jmain), dtlat(imain, jmain), istatus
+ call error_handler(E_MSG,'test_interpolation',msgstring,source,revision,revdate)
+ endif
+
write(25, *) dtlon(imain, jmain), dtlat(imain, jmain), dipole_t(imain, jmain)
end do
More information about the Dart-dev
mailing list