[Dart-dev] [3620] DART/trunk/models/MITgcm_ocean/model_mod.f90: initial version to handle observations at more levels than just 'depth'

nancy at ucar.edu nancy at ucar.edu
Mon Oct 6 15:27:19 MDT 2008


An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/dart-dev/attachments/20081006/9639e792/attachment.html
-------------- next part --------------
Modified: DART/trunk/models/MITgcm_ocean/model_mod.f90
===================================================================
--- DART/trunk/models/MITgcm_ocean/model_mod.f90	2008-10-03 22:11:57 UTC (rev 3619)
+++ DART/trunk/models/MITgcm_ocean/model_mod.f90	2008-10-06 21:27:19 UTC (rev 3620)
@@ -23,7 +23,7 @@
 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, &
-                             vert_is_level
+                             vert_is_level, vert_is_surface
 use    utilities_mod, only : register_module, error_handler, E_ERR, E_WARN, E_MSG, &
                              logfileunit, get_unit, nc_check, do_output, &
                              find_namelist_in_file, check_namelist_read
@@ -309,16 +309,14 @@
 ! Called to do one time initialization of the model. In this case,
 ! it reads in the grid information and then the model data.
 
-integer :: i, iunit, io, bogus
-real(r4), allocatable :: data_2d_array(:,:)
+integer :: i, iunit, io
 
 ! The Plan:
 !
 !   read the standard MITgcm namelist file 'data.cal' for the calendar info
 !
 !   read the standard MITgcm namelist file 'data' for the
-!   time stepping info, the grid info, and 
-!   (not that we need them,) the files for the boundary conditions, etc.
+!   time stepping info and the grid info.
 !
 !   open the grid data files to get the actual grid coordinates
 !
@@ -398,16 +396,11 @@
 read(iunit, nml = PARM05, iostat = io)
 call check_namelist_read(iunit, io, "PARM05")
 
-! we pass in a filename (without the .meta and the .data)
-! 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
-! levels/lats/lons is that i set the default value of delZ to 0.0
+! The only way I know to compute the number of
+! levels/lats/lons is to set the default value of delZ to 0.0
 ! before reading the namelist.  now loop until you get back
-! to zero and that is the end of the list.  not a very
-! satisfying solution, but maybe good enough for now?
+! to zero and that is the end of the list.
+! Not a very satisfying/robust solution ...
 
 Nx = -1
 do i=1, size(delX)
@@ -446,13 +439,10 @@
 endif
 
 ! We know enough to allocate grid variables. 
-! 'read_snapshot' checks the ".meta" files for agreeable dimensions.
 
 allocate(XC(Nx), YC(Ny), ZC(Nz))
 allocate(XG(Nx), YG(Ny), ZG(Nz))
-!allocate(data_2d_array(Nx,Ny))
 
-! determine longitudes 
 ! XG (the grid edges) and XC (the grid centroids) must be computed.
 
 XG(1) = thetaMin
@@ -462,16 +452,6 @@
  XC(i) = XC(i-1) + 0.5_r8 * delX(i-1) + 0.5_r8 * delX(i) 
 enddo
 
-!call read_snapshot('XC', data_2d_array, bogus)
-!do i=1, Nx
-!  write(*,*)'XC(',i,') = ',XC(i),data_2d_array(i, 1)
-!enddo
-
-!call read_snapshot('XG', data_2d_array, bogus)
-!do i=1, Nx
-!  write(*,*)'XG(',i,') = ',XG(i),data_2d_array(i, 1)
-!enddo
-
 ! YG (the grid edges) and YC (the grid centroids) must be computed.
 
 YG(1) = phiMin
@@ -481,18 +461,6 @@
  YC(i) = YC(i-1) + 0.5_r8 * delY(i-1) + 0.5_r8 * delY(i) 
 enddo
 
-!call read_snapshot('YC', data_2d_array, bogus)
-!do i=1, Ny
-!  write(*,*)'YC(',i,') = ',YC(i),data_2d_array(1, i)
-!enddo
-
-!call read_snapshot('YG', data_2d_array, bogus)
-!do i=1, Ny
-!  write(*,*)'YG(',i,') = ',YG(i),data_2d_array(1, i)
-!enddo
-
-!deallocate(data_2d_array)
-
 ! the namelist contains a list of thicknesses of each depth level (delZ)
 ! ZG (the grid edges) and ZC (the grid centroids) must be computed.
 
@@ -526,8 +494,6 @@
 model_size = (n3dfields * (Nx * Ny * Nz)) + (n2dfields * (Nx * Ny))
 if (do_output()) write(*,*) 'model_size = ', model_size
 
-!print *, 'end of static init model'
-
 end subroutine static_init_model
 
 
@@ -639,7 +605,7 @@
 
 ! Local storage
 real(r8)       :: loc_array(3), llon, llat, lheight
-integer        :: base_offset, offset
+integer        :: base_offset, offset, ind
 integer        :: hgt_bot, hgt_top
 real(r8)       :: hgt_fract
 real(r8)       :: top_val, bot_val
@@ -651,26 +617,28 @@
 interp_val = 0.0_r8
 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
+! Get the individual locations values
 loc_array = get_location(location)
-llon = loc_array(1)
-llat = loc_array(2)
+llon    = loc_array(1)
+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
-   if (vert_is_level(location)) then 
-      ! should range check first
-      lheight = zc(nint(loc_array(3)))
-!print *, 'setting height to ', lheight
+if( vert_is_height(location) ) then
+   ! Nothing to do 
+elseif ( vert_is_surface(location) ) then
+   ! Nothing to do 
+elseif (vert_is_level(location)) then
+   ! convert the level index to an actual depth 
+   ind = nint(loc_array(3))
+   if ( (ind < 1) .or. (ind > size(zc)) ) then 
+      lheight = zc(ind)
    else
       istatus = 1
       return
    endif
+else   ! if pressure or undefined, we don't know what to do
+   istatus = 7
+   return
 endif
 
 ! Do horizontal interpolations for the appropriate levels
@@ -694,20 +662,17 @@
 !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
-   call height_bounds(lheight, nz, zc, hgt_bot, hgt_top, hgt_fract, hstatus)
-   if(hstatus /= 0) then
-      istatus = 2
-      return
-   endif
-else
-! Sea Surface Height only has one level, no vertical interpolation
+if( vert_is_surface(location) ) then
    call lat_lon_interpolate(x(base_offset:), llon, llat, obs_type, interp_val, istatus)
    return
 endif
 
-!print *, 'height top, bot = ', hgt_top, hgt_bot
+! Get the bounding vertical levels and the fraction between bottom and top
+call height_bounds(lheight, nz, zc, hgt_bot, hgt_top, hgt_fract, hstatus)
+if(hstatus /= 0) then
+   istatus = 2
+   return
+endif
 
 ! Find the base location for the top height and interpolate horizontally on this level
 offset = base_offset + (hgt_top - 1) * nx * ny
@@ -788,7 +753,7 @@
 
 ! Subroutine to interpolate to a lat lon location given that portion of state
 ! vector. 
-! NOTE: Using array sections to pas in the x array may be inefficient on some
+! NOTE: Using array sections to pass in the x array may be inefficient on some
 ! compiler/platform setups. Might want to pass in the entire array with a base
 ! offset value instead of the section if this is an issue.
 
@@ -1173,30 +1138,10 @@
 
 function nc_write_model_atts( ncFileID ) result (ierr)
 !------------------------------------------------------------------
-! TJH 24 Oct 2006 -- Writes the model-specific attributes to a netCDF file.
+! TJH -- Writes the model-specific attributes to a netCDF file.
 !     This includes coordinate variables and some metadata, but NOT
-!     the model state vector. We do have to allocate SPACE for the model
-!     state vector, but that variable gets filled as the model advances.
+!     the model state vector.
 !
-! As it stands, this routine will work for ANY model, with no modification.
-!
-! The simplest possible netCDF file would contain a 3D field
-! containing the state of 'all' the ensemble members. This requires
-! three coordinate variables -- one for each of the dimensions 
-! [model_size, ensemble_member, time]. A little metadata is useful, 
-! so we can also create some 'global' attributes. 
-! This is what is implemented here.
-!
-! Once the simplest case is working, this routine (and nc_write_model_vars)
-! can be extended to create a more logical partitioning of the state vector,
-! fundamentally creating a netCDF file with variables that are easily 
-! plotted. The bgrid model_mod is perhaps a good one to view, keeping
-! in mind it is complicated by the fact it has two coordinate systems. 
-! There are stubs in this template, but they are only stubs.
-!
-! TJH 29 Jul 2003 -- for the moment, all errors are fatal, so the
-! return code is always '0 == normal', since the fatal errors stop execution.
-!
 ! assim_model_mod:init_diag_output uses information from the location_mod
 !     to define the location dimension and variable ID. All we need to do
 !     is query, verify, and fill ...
@@ -1776,6 +1721,7 @@
 character(len=128) :: filename, charstring
 integer :: iunit, io
 integer :: i, j, indx, nlines, dim1, dimN
+logical :: fexist
 
 if ( .not. module_initialized ) call static_init_model
 
@@ -1794,16 +1740,25 @@
 read_meta%nrecords = 0
 read_meta%timeStepNumber = timestepcount
 
+iunit = get_unit()
+
+! See if the file exists ... typically not the case for a cold start.
+! If the file does not exist, we must use the namelist values as parsed.
+
+inquire(file=filename, exist=fexist, iostat=io)
+if ((io /= 0) .or. (.not. fexist)) then 
+   write(msgstring,*) trim(filename), ' does not exist, using namelist values'
+   call error_handler(E_MSG,'model_mod:read_meta',msgstring,source,revision,revdate)
+   return
+endif
+
 ! Get next available unit number and open the file
 
-iunit = get_unit()
 open(unit=iunit, file=filename, action='read', form='formatted', iostat = io)
 if (io /= 0) then
-   write(msgstring,*) 'unable to open file ', trim(filename), ' for reading'
-   call error_handler(E_WARN,'model_mod:read_meta',msgstring,source,revision,revdate)
+   write(msgstring,*) 'cannot open ', trim(filename), ', using namelist values'
+   call error_handler(E_MSG,'model_mod:read_meta',msgstring,source,revision,revdate)
    return
-else
-   ! Use the namelist values as the defaults ??? unable to determin 2D or 3D file
 endif
 
 ! Read every line looking for the nDims entry


More information about the Dart-dev mailing list