[Dart-dev] [3288] DART/trunk/models/MITgcm_ocean/model_mod.f90: minor fixes to the interp routines; interpolates

nancy at subversion.ucar.edu nancy at subversion.ucar.edu
Fri Apr 4 09:45:07 MDT 2008


An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/dart-dev/attachments/20080404/21c3171b/attachment.html
-------------- next part --------------
Modified: DART/trunk/models/MITgcm_ocean/model_mod.f90
===================================================================
--- DART/trunk/models/MITgcm_ocean/model_mod.f90	2008-04-01 20:49:33 UTC (rev 3287)
+++ DART/trunk/models/MITgcm_ocean/model_mod.f90	2008-04-04 15:45:06 UTC (rev 3288)
@@ -18,9 +18,10 @@
 use time_manager_mod, only : time_type, set_time
 use     location_mod, only : location_type,      get_close_maxdist_init, &
                              get_close_obs_init, get_close_obs, set_location, &
-                             VERTISHEIGHT, get_location, vert_is_height
+                             VERTISHEIGHT, get_location, vert_is_height, &
+                             vert_is_level
 use    utilities_mod, only : register_module, error_handler, E_ERR, E_WARN, E_MSG, &
-                             logfileunit, get_unit, nc_check, &
+                             logfileunit, nmfileunit, get_unit, nc_check, &
                              find_namelist_in_file, check_namelist_read
 use     obs_kind_mod, only : KIND_TEMPERATURE, KIND_SALINITY, KIND_U_CURRENT_COMPONENT, &
                              KIND_V_CURRENT_COMPONENT, KIND_SEA_SURFACE_HEIGHT
@@ -51,7 +52,7 @@
 ! the interfaces here can be changed as appropriate.
 public :: prog_var_to_vector, vector_to_prog_var, &
           MIT_meta_type, read_meta, write_meta, &
-          read_snapshot, write_snapshot, &
+          read_snapshot, write_snapshot, get_gridsize, &
           snapshot_files_to_sv, sv_to_snapshot_files
 
 
@@ -65,9 +66,9 @@
 
 
 !! FIXME: on ifort, this must be 1
-!!        on xlf,   this must be 4
+!!        on xlf,   this must be 4 or 8?
 !! see the comments in the read_2d_snapshot code for more info
-integer, parameter :: item_size_direct_read = 1
+integer, parameter :: item_size_direct_read = 8
 
 !------------------------------------------------------------------
 !
@@ -325,8 +326,8 @@
 
 ! Record the namelist values used for the run
 call error_handler(E_MSG,'static_init_model','model_nml values are',' ',' ',' ')
-write(logfileunit, nml=model_nml)
-write(     *     , nml=model_nml)
+write(nmfileunit, nml=model_nml)
+write(    *     , nml=model_nml)
 
 ! Read in the MITgcm namelists from the 'data' file
 call find_namelist_in_file("data", "PARM03", iunit)
@@ -347,9 +348,8 @@
 call check_namelist_read(iunit, io, "PARM05")
 
 ! we pass in a filename (without the .meta and the .data)
-! and an unallocated array of the right dimensionality.
-! it returns the allocated, filled array, along with the timestep count.
-! when we are done, drop_snapshot() frees the space.
+! and an allocated array of the right dimensionality.
+! it returns the filled array, along with the timestep count.
 
 ! FIXME:
 ! right now the only way i know to compute the number of
@@ -462,6 +462,7 @@
 
 ! The time_step in terms of a time type must also be initialized.
 time_step = set_time(time_step_seconds, time_step_days)
+!print *, 'end of static init model'
 
 end subroutine static_init_model
 
@@ -572,6 +573,7 @@
 ! Successful istatus is 0
 istatus = 0
 
+!print *, 'top of model interpolate'
 ! Get the individual locations values; make sure the vertical is height
 ! Might want to support model level in the vertical at some point
 loc_array = get_location(location)
@@ -579,10 +581,18 @@
 llat = loc_array(2)
 lheight = loc_array(3)
 
+!print *, 'lon, lat, height = ', llon, llat, lheight
+
 ! Only know how to work with height vertical coordinate for now
 if(.not. vert_is_height(location)) then
-   istatus = 1
-   return
+   if (vert_is_level(location)) then 
+      ! should range check first
+      lheight = zc(nint(loc_array(3)))
+!print *, 'setting height to ', lheight
+   else
+      istatus = 1
+      return
+   endif
 endif
 
 ! Do horizontal interpolations for the appropriate levels
@@ -603,6 +613,8 @@
    return
 endif
 
+!print *, 'base offset now ', base_offset
+
 ! For Sea Surface Height don't need the vertical coordinate
 if(obs_type /= KIND_SEA_SURFACE_HEIGHT) then
    ! Get the bounding vertical levels and the fraction between bottom and top
@@ -617,21 +629,27 @@
    return
 endif
 
+!print *, 'height top, bot = ', hgt_top, hgt_bot
 
 ! Find the base location for the top height and interpolate horizontally on this level
 offset = base_offset + (hgt_top - 1) * nx * ny
+!print *, 'relative top height offset = ', offset - base_offset
+!print *, 'absolute top height offset = ', offset
 call lat_lon_interpolate(x(offset:), llon, llat, obs_type, top_val, istatus)
 ! Failed istatus from interpolate means give up
 if(istatus /= 0) return
 
 ! Find the base location for the bottom height and interpolate horizontally on this level
 offset = base_offset + (hgt_bot - 1) * nx * ny
+!print *, 'relative bot height offset = ', offset - base_offset
+!print *, 'absolute bot height offset = ', offset
 call lat_lon_interpolate(x(offset:), llon, llat, obs_type, bot_val, istatus)
 ! Failed istatus from interpolate means give up
 if(istatus /= 0) return
 
 ! Then weight them by the fraction and return
 interp_val = bot_val + hgt_fract * (top_val - bot_val)
+!print *, 'model_interp: interp val = ',interp_val
 
 end subroutine model_interpolate
 
@@ -757,22 +775,29 @@
 ! will be that an observation whos forward operator requires interpolating
 ! from a point that has exactly 0.0 (but is not masked) will not be 
 ! assimilated.
+!print *, 'lon_bot, lon_top = ', lon_bot, lon_top
+!print *, 'lat_bot, lat_top = ', lat_bot, lat_top
+
 pa = get_val(lon_bot, lat_bot, nx, x, masked)
+!print *, 'pa = ', pa
 if(masked) then
    istatus = 3
    return
 endif
 pb = get_val(lon_top, lat_bot, nx, x, masked)
+!print *, 'pb = ', pb
 if(masked) then
    istatus = 3
    return
 endif
 pc = get_val(lon_bot, lat_top, nx, x, masked)
+!print *, 'pc = ', pc
 if(masked) then
    istatus = 3
    return
 endif
 pd = get_val(lon_top, lat_top, nx, x, masked)
+!print *, 'pd = ', pd
 if(masked) then
    istatus = 3
    return
@@ -780,10 +805,16 @@
 
 ! Finish bi-linear interpolation 
 ! First interpolate in longitude
+!print *, 'bot lon_fract, delta = ', lon_fract, (pb-pa)
 xbot = pa + lon_fract * (pb - pa)
+!print *, 'xbot = ', xbot
+!print *, 'top lon_fract, delta = ', lon_fract, (pd-pc)
 xtop = pc + lon_fract * (pd - pc)
+!print *, 'xtop = ', xtop
 ! Now interpolate in latitude
+!print *, 'lat_fract, delta = ', lat_fract, (xtop - xbot)
 interp_val = xbot + lat_fract * (xtop - xbot)
+!print *, 'lat_lon_interp: interp_val = ', interp_val
 
 end subroutine lat_lon_interpolate
 
@@ -867,15 +898,22 @@
 ! Default is success
 istatus = 0
 
+!print *, 'computing bounds for = ', llon
 ! This is inefficient, someone could clean it up
 ! Plus, it doesn't work for a global model that wraps around
 do i = 2, nlons
    dist_bot = lon_dist(llon, lon_array(i - 1))
    dist_top = lon_dist(llon, lon_array(i))
+!print *, 'dist top, bot = ', dist_top, dist_bot
    if(dist_bot >= 0 .and. dist_top < 0) then
       bot = i - 1
       top = i
-      fract = abs(dist_bot) / (abs(dist_bot) + dist_top)
+!print *, 'bot, top = ', bot, top
+!print *, 'numerator = ',  dist_bot
+!print *, 'denomenator = ', dist_bot + abs(dist_top)
+      fract = dist_bot / (dist_bot + abs(dist_top))
+      ! orig: fract = abs(dist_bot) / (abs(dist_bot) + dist_top)
+!print *, 'fract = ', fract
       return
    endif
 end do
@@ -911,18 +949,21 @@
 end function lon_dist
 
 
-function get_val(lat_index, lon_index, nlon, x, masked)
+function get_val(lon_index, lat_index, nlon, x, masked)
 !=======================================================================
 !
 
 ! Returns the value from a single level array given the lat and lon indices
-integer,     intent(in) :: lat_index, lon_index, nlon
+integer,     intent(in) :: lon_index, lat_index, nlon
 real(r8),    intent(in) :: x(:)
 logical,    intent(out) :: masked
 real(r8)                :: get_val
 
 ! Layout has lons varying most rapidly
+!print *, 'lat_index, lon_index, nlon', lat_index, lon_index, nlon
+!print *, 'computing index val ', (lat_index - 1) * nlon + lon_index
 get_val = x((lat_index - 1) * nlon + lon_index)
+!print *, 'get_val = ', get_val
 
 ! Masked returns false if the value is masked
 ! A grid variable is assumed to be masked if its value is exactly 0.
@@ -933,6 +974,7 @@
    masked = .false.
 endif
 
+!print *, 'masked is ', masked
 end function get_val
 
 
@@ -969,7 +1011,7 @@
 real(R8) :: lat, lon, depth
 integer :: var_num, offset, lon_index, lat_index, depth_index
 
-print *, 'asking for meta data about index ', index_in
+!print *, 'asking for meta data about index ', index_in
 
 if (index_in < start_index(S_index+1)) then
    if (present(var_type)) var_type = KIND_SALINITY  
@@ -988,12 +1030,12 @@
    var_num = SSH_index
 endif
 
-print *, 'var num = ', var_num
+!print *, 'var num = ', var_num
 
 ! local offset into this var array
 offset = index_in - start_index(var_num)
 
-print *, 'offset = ', offset
+!print *, 'offset = ', offset
 
 if (var_num == SSH_index) then
   depth = 0.0
@@ -1006,11 +1048,11 @@
 lat_index = (offset - ((depth_index-1)*Nx*Ny)) / Nx + 1
 lon_index = offset - ((depth_index-1)*Nx*Ny) - ((lat_index-1)*Nx) + 1
 
-print *, 'lon, lat, depth index = ', lon_index, lat_index, depth_index
+!print *, 'lon, lat, depth index = ', lon_index, lat_index, depth_index
 lon = XC(lon_index)
 lat = YC(lat_index)
 
-print *, 'lon, lat, depth = ', lon, lat, depth
+!print *, 'lon, lat, depth = ', lon, lat, depth
 
 location = set_location(lon, lat, depth, VERTISHEIGHT)
 
@@ -2242,9 +2284,9 @@
 subroutine prog_var_to_vector(s,t,u,v,ssh,x)
 !------------------------------------------------------------------
 !
-real(r4), dimension(Nx,Ny,Nz), intent(in)  :: s,t,u,v
-real(r4), dimension(Nx,Ny),    intent(in)  :: ssh
-real(r8), dimension(:),        intent(out) :: x
+real(r4), dimension(:,:,:), intent(in)  :: s,t,u,v
+real(r4), dimension(:,:),   intent(in)  :: ssh
+real(r8), dimension(:),     intent(out) :: x
 
 integer :: i,j,k,ii
 
@@ -2424,6 +2466,17 @@
 end subroutine vector_to_3d_prog_var
 
 
+subroutine get_gridsize(num_x, num_y, num_z)
+!------------------------------------------------------------------
+!
+ integer, intent(out) :: num_x, num_y, num_z
+
+ num_x = Nx
+ num_y = Ny
+ num_z = Nz
+
+end subroutine get_gridsize
+
 !===================================================================
 ! End of model_mod
 !===================================================================


More information about the Dart-dev mailing list