[Dart-dev] [4449] DART/trunk/models/NCOMMAS: reads grid and varsizes from a restart file - in theory

nancy at ucar.edu nancy at ucar.edu
Mon Aug 2 23:03:13 MDT 2010


Revision: 4449
Author:   thoar
Date:     2010-08-02 23:03:13 -0600 (Mon, 02 Aug 2010)
Log Message:
-----------
reads grid and varsizes from a restart file - in theory

Modified Paths:
--------------
    DART/trunk/models/NCOMMAS/dart_ncommas_mod.f90
    DART/trunk/models/NCOMMAS/model_mod.f90

-------------- next part --------------
Modified: DART/trunk/models/NCOMMAS/dart_ncommas_mod.f90
===================================================================
--- DART/trunk/models/NCOMMAS/dart_ncommas_mod.f90	2010-08-02 21:31:02 UTC (rev 4448)
+++ DART/trunk/models/NCOMMAS/dart_ncommas_mod.f90	2010-08-03 05:03:13 UTC (rev 4449)
@@ -12,7 +12,6 @@
 
 use        types_mod, only : r8, rad2deg, PI, SECPERDAY
 use time_manager_mod, only : time_type, get_date, set_date, get_time, set_time, &
-                             set_calendar_type, get_calendar_string, &
                              print_date, print_time, operator(==), operator(-)
 use    utilities_mod, only : get_unit, open_file, close_file, file_exist, &
                              register_module, error_handler, nc_check, &
@@ -26,9 +25,7 @@
 implicit none
 private
 
-public :: get_ncommas_calendar, set_model_time_step, &
-          get_horiz_grid_dims, get_vert_grid_dim, &
-          read_horiz_grid, read_topography, read_vert_grid, &
+public :: set_model_time_step, grid_type, get_grid_dims, get_grid, &
           write_ncommas_namelist, get_ncommas_restart_filename
 
 ! version controlled file description for error handling, do not edit
@@ -37,320 +34,249 @@
    revision = '$Revision$', &
    revdate  = '$Date$'
 
-character(len=256) :: msgstring
+character(len=256) :: string1, string2
 logical, save :: module_initialized = .false.
 
-character(len=256) :: ic_filename      = 'ncommas.r.nc'
-!character(len=256) :: restart_filename = 'dart_ncommas_mod_restart_filename_not_set'
-
 ! set this to true if you want to print out the current time
 ! after each N observations are processed, for benchmarking.
 logical :: print_timestamps = .false.
 integer :: print_every_Nth  = 10000
 
-!------------------------------------------------------------------
-! The ncommas time manager namelist variables
-!------------------------------------------------------------------
+type grid_type
+   private
+   integer  :: nx, ny, nz         ! determines if by level, height, pressure, ...
+   real(r8), pointer ::  lon(:,:) ! lon stored in radians
+   real(r8), pointer ::  lat(:,:) ! lat stored in radians
+   real(r8), pointer :: vloc(:,:) ! height stored in meters
+end type grid_type
 
-character(len=100) :: accel_file ! length consistent with ncommas
-character(len= 64) :: stop_option, runid, dt_option, time_mix_opt
-character(len=  1) :: date_separator
-logical  :: impcor, laccel, allow_leapyear
-real(r8) :: dtuxcel, dt_count
-integer  :: iyear0, imonth0, iday0, ihour0, iminute0, isecond0
-integer  :: stop_count, fit_freq, time_mix_freq
+type(grid_type) :: Ugrid
+type(grid_type) :: Vgrid
+type(grid_type) :: Wgrid
+type(grid_type) ::  grid
 
-namelist /time_manager_nml/ runid, time_mix_opt, time_mix_freq, &
-    impcor, laccel, accel_file, dtuxcel, iyear0, imonth0, &
-    iday0, ihour0, iminute0, isecond0, dt_option, dt_count, &
-    stop_option, stop_count, date_separator, allow_leapyear, fit_freq
-
 !------------------------------------------------------------------
-! The ncommas I/O namelist variables
+! The ncommas restart manager namelist variables
 !------------------------------------------------------------------
 
-character(len=100) :: log_filename, pointer_filename
-logical :: lredirect_stdout, luse_pointer_files
-logical :: luse_nf_64bit_offset
-integer :: num_iotasks
+character(len=256) :: ic_filename      = 'ncommas.r.nc'
+!character(len=256) :: restart_filename = 'dart_ncommas_mod_restart_filename_not_set'
+character(len= 64) :: ew_boundary_type, ns_boundary_type
 
-namelist /io_nml/ num_iotasks, lredirect_stdout, log_filename, &
-    luse_pointer_files, pointer_filename, luse_nf_64bit_offset
+namelist /restart_nml/ ic_filename, ew_boundary_type, ns_boundary_type
 
-!------------------------------------------------------------------
-! The ncommas restart manager namelist variables
-!------------------------------------------------------------------
 
-character(len=100) :: restart_outfile ! length consistent with ncommas
-character(len= 64) :: restart_freq_opt, restart_start_opt, restart_fmt
-logical :: leven_odd_on, pressure_correction
-integer :: restart_freq, restart_start, even_odd_freq
+!======================================================================
+contains
+!======================================================================
 
-namelist /restart_nml/ restart_freq_opt, restart_freq, &
-     restart_start_opt, restart_start, restart_outfile, &
-    restart_fmt, leven_odd_on, even_odd_freq, pressure_correction
 
+subroutine initialize_module
 !------------------------------------------------------------------
-! The ncommas initial temperature and salinity namelist 
-!------------------------------------------------------------------
+integer :: iunit, io
 
-character(len=100) :: init_ts_file ! length consistent with ncommas
-character(len=100) :: init_ts_outfile
-character(len= 64) :: init_ts_option, init_ts_suboption
-character(len= 64) :: init_ts_file_fmt, init_ts_outfile_fmt
+! Make sure we have a ncommas restart file (for grid dims)
+if ( .not. file_exist(ic_filename) ) then
+   string1 = trim(ic_filename)//' not found'
+   call error_handler(E_ERR,'initialize_module', &
+          string1, source, revision, revdate)
+endif
 
-namelist /init_ts_nml/ init_ts_option, init_ts_suboption, &
-                       init_ts_file, init_ts_file_fmt, & 
-                       init_ts_outfile, init_ts_outfile_fmt
+module_initialized = .true.
 
-!------------------------------------------------------------------
-! The ncommas domain namelist 
-!------------------------------------------------------------------
+! Print module information to log file and stdout.
+call register_module(source, revision, revdate)
 
-character(len= 64) :: clinic_distribution_type, tropic_distribution_type
-character(len= 64) :: ew_boundary_type, ns_boundary_type
-integer :: nprocs_clinic, nprocs_tropic
+end subroutine initialize_module
 
-namelist /domain_nml/ clinic_distribution_type, nprocs_clinic, &
-                      tropic_distribution_type, nprocs_tropic, &
-                      ew_boundary_type, ns_boundary_type
 
+
+subroutine get_grid_dims(NXC, NXE, NYC, NYE, NZC, NZE)
 !------------------------------------------------------------------
-! The ncommas grid info namelist 
-!------------------------------------------------------------------
 !
-! ncommas grid information comes in several files:
-!   horizontal grid lat/lons in one, 
-!   topography (lowest valid vert level) in another, and 
-!   the vertical grid spacing in a third.
+! Read the grid dimensions from the restart netcdf file.
 !
-!------------------------------------------------------------------
-!
-! Here is what we can get from the (binary) horiz grid file:
-!   real (r8), dimension(:,:), allocatable :: &
-!      ULAT,            &! latitude  (radians) of U points
-!      ULON,            &! longitude (radians) of U points
-!      HTN ,            &! length (cm) of north edge of T box
-!      HTE ,            &! length (cm) of east  edge of T box
-!      HUS ,            &! length (cm) of south edge of U box
-!      HUW ,            &! length (cm) of west  edge of U box
-!      ANGLE             ! angle
-!
-! Here is what we can get from the topography file:
-!   integer, dimension(:,:), allocatable :: &
-!      KMT               ! k index of deepest grid cell on T grid
-!
-! These must be derived or come from someplace else ...
-!      KMU               ! k index of deepest grid cell on U grid
-!      HT                ! real(r8) value of deepest valid T depth (in cm)
-!      HU                ! real(r8) value of deepest valid U depth (in cm)
-!
-! The vert grid file is ascii, with 3 columns/line:
-!    cell thickness(in cm)   cell center(in m)   cell bottom(in m)
-!
-!------------------------------------------------------------------
+! The file name comes from module storage ... namelist.
 
-character(len=100) :: horiz_grid_file, vert_grid_file, &
-                      topography_file, topography_outfile, &
-                      bathymetry_file, region_info_file, &
-                      bottom_cell_file, region_mask_file
-character(len= 64) :: horiz_grid_opt, sfc_layer_opt, vert_grid_opt, &
-                      topography_opt
-logical :: partial_bottom_cells, topo_smooth, flat_bottom, lremove_points
-integer :: kmt_kmin, n_topo_smooth
+integer, intent(out) :: NXC   ! Number of Longitude centers
+integer, intent(out) :: NXE   ! Number of Longitude edges
+integer, intent(out) :: NYC   ! Number of Latitude  centers
+integer, intent(out) :: NYE   ! Number of Latitude  edges
+integer, intent(out) :: NZC   ! Number of Vertical grid centers
+integer, intent(out) :: NZE   ! Number of Vertical grid edges
 
-namelist /grid_nml/ horiz_grid_opt, horiz_grid_file, sfc_layer_opt, &
-    vert_grid_opt, vert_grid_file, topography_opt, kmt_kmin, &
-    topography_file, topography_outfile, bathymetry_file, &
-    partial_bottom_cells, bottom_cell_file, n_topo_smooth, &
-    region_mask_file, topo_smooth, flat_bottom, lremove_points, &
-    region_info_file
+integer :: grid_id, dimid
 
-!======================================================================
-contains
-!======================================================================
+if ( .not. module_initialized ) call initialize_module
 
+! get the ball rolling ...
 
-subroutine initialize_module
-!------------------------------------------------------------------
-integer :: iunit, io
+call nc_check(nf90_open(trim(ic_filename), nf90_nowrite, grid_id), &
+            'get_grid_dims','open '//trim(ic_filename))
 
-! Read ncommas calendar information
-! In 'restart' mode, this is primarily the calendar type and 'stop'
-! information. The time attributes of the restart file override 
-! the namelist time information. 
+! Longitudes : get dimid for 'XC' and then get value
 
-call find_namelist_in_file('ncommas_in', 'time_manager_nml', iunit)
-read(iunit, nml = time_manager_nml, iostat = io)
-call check_namelist_read(iunit, io, 'time_manager_nml')
+call nc_check(nf90_inq_dimid(grid_id, 'XC', dimid), &
+            'get_grid_dims','inq_dimid XC '//trim(ic_filename))
+call nc_check(nf90_inquire_dimension(grid_id, dimid, len=NXC), &
+            'get_grid_dims','inquire_dimension XC '//trim(ic_filename))
 
-! FIXME : Real observations are always GREGORIAN dates ...
-! but stomping on that here gets in the way of running
-! a perfect_model experiment for pre-1601 AD cases.
+! Longitudes : get dimid for 'XE and then get value
 
-! STOMP if ( allow_leapyear ) then
-   call set_calendar_type('gregorian')
-! STOMP else
-! STOMP    call set_calendar_type('noleap')
-! STOMP endif
+call nc_check(nf90_inq_dimid(grid_id, 'XE', dimid), &
+            'get_grid_dims','inq_dimid XE '//trim(ic_filename))
+call nc_check(nf90_inquire_dimension(grid_id, dimid, len=NXE), &
+            'get_grid_dims','inquire_dimension XE '//trim(ic_filename))
 
-! Read ncommas I/O information (for restart file ... grid dimensions)
-! Read ncommas initial information (for input/restart filename)
+! Latitudes : get dimid for 'YC' and then get value
 
-call find_namelist_in_file('ncommas_in', 'io_nml', iunit)
-read(iunit, nml = io_nml, iostat = io)
-call check_namelist_read(iunit, io, 'io_nml')
+call nc_check(nf90_inq_dimid(grid_id, 'YC', dimid), &
+            'get_grid_dims','inq_dimid YC '//trim(ic_filename))
+call nc_check(nf90_inquire_dimension(grid_id, dimid, len=NYC), &
+            'get_grid_dims','inquire_dimension YC '//trim(ic_filename))
 
-call find_namelist_in_file('ncommas_in', 'init_ts_nml', iunit)
-read(iunit, nml = init_ts_nml, iostat = io)
-call check_namelist_read(iunit, io, 'init_ts_nml')
+! Latitudes : get dimid for 'YE' and then get value
 
-! Is it a pointer file or not ...
-!if ( luse_pointer_files ) then
-!
-!   restart_filename = trim(pointer_filename)//'.restart'
-!
-!   if ( .not. file_exist(restart_filename) ) then
-!      msgstring = 'ncommas_in:pointer file '//trim(restart_filename)//' not found'
-!      call error_handler(E_ERR,'initialize_module', &
-!             msgstring, source, revision, revdate)
-!   endif
-!
-!   iunit = open_file(restart_filename,'formatted')
-!   read(iunit,'(A)')ic_filename
-!
-!   restart_filename = ' '  
-!   write(*,*)'DEBUG ... pointer filename dereferenced to ',trim(ic_filename )
-!
-!else
-!   ic_filename = trim(init_ts_file)//'.'//trim(init_ts_file_fmt)
-!endif
+call nc_check(nf90_inq_dimid(grid_id, 'YE', dimid), &
+            'get_grid_dims','inq_dimid YE '//trim(ic_filename))
+call nc_check(nf90_inquire_dimension(grid_id, dimid, len=NYE), &
+            'get_grid_dims','inquire_dimension YE '//trim(ic_filename))
 
-! Make sure we have a ncommas restart file (for grid dims)
-if ( .not. file_exist(ic_filename) ) then
-   msgstring = 'ncommas_in:init_ts_file '//trim(ic_filename)//' not found'
-   call error_handler(E_ERR,'initialize_module', &
-          msgstring, source, revision, revdate)
-endif
+! Vertical Levels : get dimid for 'ZC' and then get value
 
-! Read ncommas restart information (for model timestepping)
-call find_namelist_in_file('ncommas_in', 'restart_nml', iunit)
-read(iunit, nml = restart_nml, iostat = io)
-call check_namelist_read(iunit, io, 'restart_nml')
+call nc_check(nf90_inq_dimid(grid_id, 'ZC', dimid), &
+            'get_grid_dims','inq_dimid ZC '//trim(ic_filename))
+call nc_check(nf90_inquire_dimension(grid_id, dimid, len=NZC), &
+            'get_grid_dims','inquire_dimension ZC '//trim(ic_filename))
 
-! Read ncommas domain information (for lon wrapping or not)
-call find_namelist_in_file('ncommas_in', 'domain_nml', iunit)
-read(iunit, nml = domain_nml, iostat = io)
-call check_namelist_read(iunit, io, 'domain_nml')
+! Vertical Levels : get dimid for 'ZE' and then get value
 
-! Read ncommas grid information (for grid filenames)
-call find_namelist_in_file('ncommas_in', 'grid_nml', iunit)
-read(iunit, nml = grid_nml, iostat = io)
-call check_namelist_read(iunit, io, 'grid_nml')
+call nc_check(nf90_inq_dimid(grid_id, 'ZE', dimid), &
+            'get_grid_dims','inq_dimid ZE '//trim(ic_filename))
+call nc_check(nf90_inquire_dimension(grid_id, dimid, len=NZE), &
+            'get_grid_dims','inquire_dimension ZE '//trim(ic_filename))
 
-module_initialized = .true.
+! tidy up
 
-! Print module information to log file and stdout.
-call register_module(source, revision, revdate)
+call nc_check(nf90_close(grid_id), &
+         'get_grid_dims','close '//trim(ic_filename) )
 
-end subroutine initialize_module
+end subroutine get_grid_dims
 
 
 
-subroutine get_horiz_grid_dims(Nx, Ny)
+subroutine get_grid(NXC, NXE, NYC, NYE, NZC, NZE, &
+                    ULAT, ULON, VLAT, VLON, WLAT, WLON, ZC, ZE)
 !------------------------------------------------------------------
-! subroutine get_horiz_grid_dims(Nx, Ny)
 !
-! Read the lon, lat grid size from the restart netcdf file.
-! The actual grid file is a binary file with no header information.
+! Read the grid dimensions from the restart netcdf file.
 !
 ! The file name comes from module storage ... namelist.
 
-integer, intent(out) :: Nx   ! Number of Longitudes
-integer, intent(out) :: Ny   ! Number of Latitudes
+integer, intent(in) :: NXC   ! Number of Longitude centers
+integer, intent(in) :: NXE   ! Number of Longitude edges
+integer, intent(in) :: NYC   ! Number of Latitude  centers
+integer, intent(in) :: NYE   ! Number of Latitude  edges
+integer, intent(in) :: NZC   ! Number of Vertical grid centers
+integer, intent(in) :: NZE   ! Number of Vertical grid edges
 
-integer :: grid_id, dimid, nc_rc
+real(r8), dimension(:,:), intent(out) :: ULAT, ULON, VLAT, VLON, WLAT, WLON
+real(r8), dimension( : ), intent(out) :: ZC, ZE
 
-if ( .not. module_initialized ) call initialize_module
+! type(grid_type), intent(out) :: Ugrid  ! (ZC, YC, XE)
+! type(grid_type), intent(out) :: Vgrid  ! (ZC, YE, XC)
+! type(grid_type), intent(out) :: Wgrid  ! (ZE, YC, XC)
+! type(grid_type), intent(out) ::  grid  ! (ZC, YC, XC)
 
-! get the ball rolling ...
+real(r8), dimension(NXC) :: XC
+real(r8), dimension(NXE) :: XE
+real(r8), dimension(NYC) :: YC
+real(r8), dimension(NYE) :: YE
 
-call nc_check(nf90_open(trim(ic_filename), nf90_nowrite, grid_id), &
-         'get_horiz_grid_dims','open '//trim(ic_filename))
+! real(r8), dimension(nx,ny), intent(out) :: ULAT, ULON, TLAT, TLON
 
-! Longitudes : get dimid for 'i' or 'nlon', and then get value
-nc_rc = nf90_inq_dimid(grid_id, 'i', dimid)
-if (nc_rc /= nf90_noerr) then
-   nc_rc = nf90_inq_dimid(grid_id, 'nlon', dimid)
-   if (nc_rc /= nf90_noerr) then
-      msgstring = "unable to find either 'i' or 'nlon' in file "//trim(ic_filename)
-      call error_handler(E_ERR, 'get_horiz_grid_dims', msgstring, &
-                         source,revision,revdate) 
-   endif
-endif
+integer, dimension(NF90_MAX_VAR_DIMS) :: dimIDs
+character(len=NF90_MAX_NAME)          :: varname
+integer                               :: VarID, numdims, dimlen
+integer                               :: ncid, dimid
 
-call nc_check(nf90_inquire_dimension(grid_id, dimid, len=Nx), &
-         'get_horiz_grid_dims','inquire_dimension i '//trim(ic_filename))
 
-! Latitudes : get dimid for 'j' or 'nlat' ... and then get value
-nc_rc = nf90_inq_dimid(grid_id, 'j', dimid)
-if (nc_rc /= nf90_noerr) then
-   nc_rc = nf90_inq_dimid(grid_id, 'nlat', dimid)
-   if (nc_rc /= nf90_noerr) then
-      msgstring = "unable to find either 'j' or 'nlat' in "//trim(ic_filename)
-      call error_handler(E_ERR, 'get_horiz_grid_dims', msgstring, &
-                         source,revision,revdate)
-   endif
-endif
+if ( .not. module_initialized ) call initialize_module
 
-call nc_check(nf90_inquire_dimension(grid_id, dimid, len=Ny), &
-         'get_horiz_grid_dims','inquire_dimension i '//trim(ic_filename))
+! get the ball rolling ...
 
-! tidy up
+call nc_check(nf90_open(trim(ic_filename), nf90_nowrite, ncid), 'get_grid', 'open '//trim(ic_filename))
 
-call nc_check(nf90_close(grid_id), &
-         'get_horiz_grid_dims','close '//trim(ic_filename) )
 
-end subroutine get_horiz_grid_dims
+! fixme - in a perfect world - 
+! Get the variable ID
+! Check to make sure it is the right shape
+! Read it
+call nc_check(nf90_inq_varid(ncid, 'XC', VarID), 'get_grid', 'inq_varid XC '//trim(ic_filename))
+!call nc_check(nf90_inquire_variable(ncid, VarId, dimids=dimIDs, ndims=numdims), &
+!                        'get_grid', 'inquire_variable XC '//trim(ic_filename))
+call nc_check(nf90_get_var(ncid, VarID, XC), 'get_grid', 'get_var XC '//trim(ic_filename))
 
+call nc_check(nf90_inq_varid(ncid, 'XE', VarID), 'get_grid', 'inq_varid XE '//trim(ic_filename))
+!call nc_check(nf90_inquire_variable(ncid, VarId, dimids=dimIDs, ndims=numdims), &
+!                        'get_grid', 'inquire_variable XE '//trim(ic_filename))
+call nc_check(nf90_get_var(ncid, VarID, XE), 'get_grid', 'get_var XE '//trim(ic_filename))
 
+call nc_check(nf90_inq_varid(ncid, 'YC', VarID), 'get_grid', 'inq_varid YC '//trim(ic_filename))
+!call nc_check(nf90_inquire_variable(ncid, VarId, dimids=dimIDs, ndims=numdims), &
+!                        'get_grid', 'inquire_variable YC '//trim(ic_filename))
+call nc_check(nf90_get_var(ncid, VarID, YC), 'get_grid', 'get_var YC '//trim(ic_filename))
 
-  subroutine get_vert_grid_dim(Nz)
-!------------------------------------------------------------------
-! subroutine get_vert_grid_dim(Nz)
-!
-! count the number of lines in the ascii file to figure out max
-! number of vert blocks.
+call nc_check(nf90_inq_varid(ncid, 'YE', VarID), 'get_grid', 'inq_varid YE '//trim(ic_filename))
+!call nc_check(nf90_inquire_variable(ncid, VarId, dimids=dimIDs, ndims=numdims), &
+!                        'get_grid', 'inquire_variable YE '//trim(ic_filename))
+call nc_check(nf90_get_var(ncid, VarID, YE), 'get_grid', 'get_var YE '//trim(ic_filename))
 
-integer, intent(out) :: Nz
+call nc_check(nf90_inq_varid(ncid, 'ZC', VarID), 'get_grid', 'inq_varid ZC '//trim(ic_filename))
+!call nc_check(nf90_inquire_variable(ncid, VarId, dimids=dimIDs, ndims=numdims), &
+!                        'get_grid', 'inquire_variable ZC '//trim(ic_filename))
+call nc_check(nf90_get_var(ncid, VarID, ZC), 'get_grid', 'get_var ZC '//trim(ic_filename))
 
-integer :: linelen ! disposable
+call nc_check(nf90_inq_varid(ncid, 'ZE', VarID), 'get_grid', 'inq_varid ZE '//trim(ic_filename))
+!call nc_check(nf90_inquire_variable(ncid, VarId, dimids=dimIDs, ndims=numdims), &
+!                        'get_grid', 'inquire_variable ZE '//trim(ic_filename))
+call nc_check(nf90_get_var(ncid, VarID, ZE), 'get_grid', 'get_var ZE '//trim(ic_filename))
 
-if ( .not. module_initialized ) call initialize_module
+! FIXME - need to convert these things to radians to store in the grid variables.
+! FIXME - need to allocate the pointers in the grid variables, etc.
 
-call find_textfile_dims(vert_grid_file, Nz, linelen)
+! call calc_tpoints(nx, ny, ULAT, ULON, TLAT, TLON)
 
-end subroutine get_vert_grid_dim
+! convert from radians to degrees
 
+!ULAT = ULAT * rad2deg
+!ULON = ULON * rad2deg
+!TLAT = TLAT * rad2deg
+!TLON = TLON * rad2deg
 
-   
-subroutine get_ncommas_calendar(calstring)
-!------------------------------------------------------------------
-! the initialize_module ensures that the ncommas namelists are read and 
-! the DART time manager gets the ncommas calendar setting.
+! ensure [0,360) [-90,90]
+
+! where (ULON <   0.0_r8) ULON = ULON + 360.0_r8
+! where (ULON > 360.0_r8) ULON = ULON - 360.0_r8
+! where (TLON <   0.0_r8) TLON = TLON + 360.0_r8
+! where (TLON > 360.0_r8) TLON = TLON - 360.0_r8
 !
-! Then, the DART time manager is queried to return what it knows ...
-!
-character(len=*), INTENT(OUT) :: calstring
+! where (ULAT < -90.0_r8) ULAT = -90.0_r8
+! where (ULAT >  90.0_r8) ULAT =  90.0_r8
+! where (TLAT < -90.0_r8) TLAT = -90.0_r8
+! where (TLAT >  90.0_r8) TLAT =  90.0_r8
 
-if ( .not. module_initialized ) call initialize_module
+! tidy up
 
-call get_calendar_string(calstring)
+call nc_check(nf90_close(ncid), 'get_grid','close '//trim(ic_filename) )
 
-end subroutine get_ncommas_calendar
+end subroutine get_grid
 
 
 
+
+
+
 function set_model_time_step()
 !------------------------------------------------------------------
 ! the initialize_module ensures that the ncommas namelists are read.
@@ -361,18 +287,9 @@
 
 if ( .not. module_initialized ) call initialize_module
 
-! Check the 'restart_freq_opt' and 'restart_freq' to determine
-! when we can stop the model
+! FIXME - determine when we can stop the model
 
-if ( trim(restart_freq_opt) == 'nday' ) then
-   set_model_time_step = set_time(0, restart_freq) ! (seconds, days)
-else if ( trim(restart_freq_opt) == 'nyear' ) then
-   ! FIXME ... CCSM_ncommas uses a bogus value for this
    set_model_time_step = set_time(0, 1) ! (seconds, days)
-else
-   call error_handler(E_ERR,'set_model_time_step', &
-              'restart_freq_opt must be nday', source, revision, revdate)
-endif
 
 end function set_model_time_step
 
@@ -393,8 +310,8 @@
 call get_time(offset, secs, days)
 
 if (secs /= 0 ) then
-   write(msgstring,*)'adv_to_time has seconds == ',secs,' must be zero'
-   call error_handler(E_ERR,'write_ncommas_namelist', msgstring, source, revision, revdate)
+   write(string1,*)'adv_to_time has seconds == ',secs,' must be zero'
+   call error_handler(E_ERR,'write_ncommas_namelist', string1, source, revision, revdate)
 endif
 
 ! call print_date( model_time,'write_ncommas_namelist:dart model date')
@@ -406,15 +323,15 @@
 
 !Convey the information to the namelist 'stop option' and 'stop count'
 
-if ( trim(stop_option) == 'nday' ) then
-   stop_count = days
-else
+!if ( trim(stop_option) == 'nday' ) then
+!   stop_count = days
+!else
    call error_handler(E_ERR,'write_ncommas_namelist', &
               'stop_option must be "nday"', source, revision, revdate)
-endif
+!endif
 
 iunit = open_file('ncommas_in.DART',form='formatted',action='rewind')
-write(iunit, nml=time_manager_nml)
+write(iunit, nml=restart_nml)
 write(iunit, '('' '')')
 close(iunit)
 
@@ -422,75 +339,8 @@
 
 
 
-  subroutine read_horiz_grid(nx, ny, ULAT, ULON, TLAT, TLON)
-!------------------------------------------------------------------
-! subroutine read_horiz_grid(nx, ny, ULAT, ULON, TLAT, TLON)
-!
-! Open and read the binary grid file
 
-integer,                    intent(in)  :: nx, ny
-real(r8), dimension(nx,ny), intent(out) :: ULAT, ULON, TLAT, TLON
 
-!real(r8), dimension(nx,ny) :: &
-!     HTN ,  &! length (cm) of north edge of T box
-!     HTE ,  &! length (cm) of east  edge of T box
-!     HUS ,  &! length (cm) of south edge of U box
-!     HUW ,  &! length (cm) of west  edge of U box
-!     ANGLE  ! angle
-
-integer :: grid_unit, reclength
-
-if ( .not. module_initialized ) call initialize_module
-
-! Check to see that the file exists.
-
-if ( .not. file_exist(horiz_grid_file) ) then
-   msgstring = 'ncommas_in:horiz_grid_file '//trim(horiz_grid_file)//' not found'
-   call error_handler(E_ERR,'read_horiz_grid', &
-          msgstring, source, revision, revdate)
-endif
-
-! Open it and read them in the EXPECTED order.
-! Actually, we only need the first two, so I'm skipping the rest.
-
-grid_unit = get_unit()
-INQUIRE(iolength=reclength) ULAT
-
-open(grid_unit, file=trim(horiz_grid_file), form='unformatted', &
-            access='direct', recl=reclength, status='old', action='read' )
-read(grid_unit, rec=1) ULAT
-read(grid_unit, rec=2) ULON
-!read(grid_unit, rec=3) HTN
-!read(grid_unit, rec=4) HTE
-!read(grid_unit, rec=5) HUS
-!read(grid_unit, rec=6) HUW
-!read(grid_unit, rec=7) ANGLE
-close(grid_unit)
-
-call calc_tpoints(nx, ny, ULAT, ULON, TLAT, TLON)
-
-! convert from radians to degrees
-
-ULAT = ULAT * rad2deg
-ULON = ULON * rad2deg
-TLAT = TLAT * rad2deg
-TLON = TLON * rad2deg
-
-! ensure [0,360) [-90,90]
-
-where (ULON <   0.0_r8) ULON = ULON + 360.0_r8
-where (ULON > 360.0_r8) ULON = ULON - 360.0_r8
-where (TLON <   0.0_r8) TLON = TLON + 360.0_r8
-where (TLON > 360.0_r8) TLON = TLON - 360.0_r8
-
-where (ULAT < -90.0_r8) ULAT = -90.0_r8
-where (ULAT >  90.0_r8) ULAT =  90.0_r8
-where (TLAT < -90.0_r8) TLAT = -90.0_r8
-where (TLAT >  90.0_r8) TLAT =  90.0_r8
-
-end subroutine read_horiz_grid
-
-
   subroutine calc_tpoints(nx, ny, ULAT, ULON, TLAT, TLON)
 !------------------------------------------------------------------
 ! subroutine calc_tpoints(nx, ny, ULAT, ULON, TLAT, TLON)
@@ -591,103 +441,17 @@
    TLON(1,:) = (TLON(2,:) + TLON(nx,:))/c2
 
 else
-   write(msgstring,'(''ncommas_in&domain_nml:ew_boundary_type '',a,'' unknown.'')') &
+   write(string1,'(''ncommas_in&domain_nml:ew_boundary_type '',a,'' unknown.'')') &
                                     trim(ew_boundary_type)
-   call error_handler(E_ERR,'calc_tpoints',msgstring,source,revision,revdate)
+   call error_handler(E_ERR,'calc_tpoints',string1,source,revision,revdate)
 endif
 
 end subroutine calc_tpoints
 
 
-
-  subroutine read_topography(nx, ny, KMT, KMU)
 !------------------------------------------------------------------
-! subroutine read_topography(nx, ny, KMT, KMU)
-!
-! Open and read the binary topography file
 
-integer,                   intent(in)  :: nx, ny
-integer, dimension(nx,ny), intent(out) :: KMT, KMU
 
-integer  :: i, j, topo_unit, reclength
-
-if ( .not. module_initialized ) call initialize_module
-
-! Check to see that the file exists.
-
-if ( .not. file_exist(topography_file) ) then
-   msgstring = 'ncommas_in:topography_file '//trim(topography_file)//' not found'
-   call error_handler(E_ERR,'read_topography', &
-          msgstring, source, revision, revdate)
-endif
-
-! read the binary file
-
-topo_unit = get_unit()
-INQUIRE(iolength=reclength) KMT
-
-open( topo_unit, file=trim(topography_file), form='unformatted', &
-            access='direct', recl=reclength, status='old', action='read' )
-read( topo_unit, rec=1) KMT
-close(topo_unit)
-
-KMU(1, 1) = 0
-do j=2,ny
-do i=2,nx
-   KMU(i,j) = min(KMT(i, j), KMT(i-1, j), KMT(i, j-1), KMT(i-1, j-1))
-enddo
-enddo
-
-end subroutine read_topography
-
-
-
-  subroutine read_vert_grid(nz, ZC, ZG)
-!------------------------------------------------------------------
-! subroutine read_vert_grid(nz, ZC, ZG)
-!
-! Open and read the ASCII vertical grid information
-!
-! The vert grid file is ascii, with 3 columns/line:
-!    cell thickness(in cm)   cell center(in m)   cell bottom(in m)
-
-integer,  intent(in)  :: nz
-real(r8), intent(out) :: ZC(nz), ZG(nz)
-
-integer  :: iunit, i, ios
-real(r8) :: depth
-
-if ( .not. module_initialized ) call initialize_module
-
-! Check to see that the file exists.
-
-if ( .not. file_exist(vert_grid_file) ) then
-   msgstring = 'ncommas_in:vert_grid_file '//trim(vert_grid_file)//' not found'
-   call error_handler(E_ERR,'read_vert_grid', &
-          msgstring, source, revision, revdate)
-endif
-
-! read the ASCII file
-
-iunit = open_file(trim(vert_grid_file), action = 'read')
-
-do i=1, nz
-
-   read(iunit,*,iostat=ios) depth, ZC(i), ZG(i)
-
-   if ( ios /= 0 ) then ! error
-      write(msgstring,*)'error reading depths, line ',i
-      call error_handler(E_ERR,'read_vert_grid',msgstring,source,revision,revdate)
-   endif
-
-enddo
-
-end subroutine read_vert_grid
-
-
-!------------------------------------------------------------------
-
-
 subroutine get_ncommas_restart_filename( filename )
 character(len=*), intent(OUT) :: filename
 

Modified: DART/trunk/models/NCOMMAS/model_mod.f90
===================================================================
--- DART/trunk/models/NCOMMAS/model_mod.f90	2010-08-02 21:31:02 UTC (rev 4448)
+++ DART/trunk/models/NCOMMAS/model_mod.f90	2010-08-03 05:03:13 UTC (rev 4449)
@@ -55,10 +55,9 @@
 
 use    random_seq_mod, only: random_seq_type, init_random_seq, random_gaussian
 
-use  dart_ncommas_mod, only: set_model_time_step,                              &
-                             get_horiz_grid_dims, get_vert_grid_dim,           &
-                             read_horiz_grid, read_vert_grid,                  &
-                             get_ncommas_restart_filename
+use  dart_ncommas_mod, only: set_model_time_step, grid_type,    &
+                             get_grid_dims, get_grid,           &
+                             get_ncommas_restart_filename, write_ncommas_namelist
 
 use typesizes
 use netcdf 
@@ -88,7 +87,7 @@
 ! generally useful routines for various support purposes.
 ! the interfaces here can be changed as appropriate.
 public :: get_gridsize, restart_file_to_sv, sv_to_restart_file, &
-          get_ncommas_restart_filename, test_interpolation
+          get_ncommas_restart_filename
 
 ! version controlled file description for error handling, do not edit
 character(len=128), parameter :: &
@@ -96,30 +95,28 @@
    revision = '$Revision$', &
    revdate  = '$Date$'
 
-character(len=256) :: msgstring
+character(len=256) :: string1, string2
 logical, save :: module_initialized = .false.
 
+character(len=256) :: ncommas_filename
+
 ! Storage for a random sequence for perturbing a single initial state
 type(random_seq_type) :: random_seq
 
 ! things which can/should be in the model_nml
-logical  :: output_state_vector = .true.
 integer  :: assimilation_period_days = 1
 integer  :: assimilation_period_seconds = 0
 real(r8) :: model_perturbation_amplitude = 0.2
-logical  :: update_dry_cell_walls = .false.
+logical  :: output_state_vector = .true.
 integer  :: debug = 0   ! turn up for more and more debug messages
+character(len=32):: calendar
 
-! FIXME: currently the update_dry_cell_walls namelist value DOES
-! NOTHING.  it needs additional code to detect the cells which are
-! wet, but within 1 cell of the bottom/sides/etc.  
-
 namelist /model_nml/  &
    output_state_vector,         &
    assimilation_period_days,    &  ! for now, this is the timestep
    assimilation_period_seconds, &
    model_perturbation_amplitude,&
-   update_dry_cell_walls,       &
+   calendar,                    &
    debug
 
 !------------------------------------------------------------------
@@ -143,20 +140,31 @@
 !  QS   long_name = "SNOW MIXING RATIO"     float  QS(TIME, ZC, YC, XC)
 !  QH   long_name = "GRAUPEL MIXING RATIO"  float  QH(TIME, ZC, YC, XC)
 !
-! FIXME: we make this completely namelist driven,
-!        both contents and order of vars.  this should
-!        wait until restart files are in netcdf format,
-!        to avoid problems with incompatible namelist
-!        and IC files.  it also complicates the mapping
-!        to and from the vars to state vector.
+! FIXME: make this completely namelist driven,
+!        both contents and order of vars.
 !------------------------------------------------------------------
 
 integer, parameter :: n3dfields = 13
 integer, parameter :: n2dfields = 0
 integer, parameter :: nfields   = n3dfields + n2dfields
 
-! (the absoft compiler likes them to all be the same length during declaration)
-! we trim the blanks off before use anyway, so ...
+! Everything needed to describe a variable
+type progvartype
+   private
+   character(len=NF90_MAX_NAME) :: varname
+   character(len=NF90_MAX_NAME) :: long_name
+   character(len=NF90_MAX_NAME) :: units
+   integer, dimension(NF90_MAX_VAR_DIMS) :: dimlens
+   integer :: numdims
+   integer :: varsize     ! prod(dimlens(1:numdims))
+   integer :: index1      ! location in dart state vector of first occurrence
+   integer :: indexN      ! location in dart state vector of last  occurrence
+   integer :: dart_kind
+   character(len=32) :: kind_string
+end type progvartype
+
+type(progvartype), dimension(nfields) :: progvar
+
 character(len=128) :: progvarnames(nfields) = &
                          (/ 'U   ', 'V   ', 'W   ', 'TH  ', 'DBZ ', &
                             'WZ  ', 'PI  ', 'QV  ', 'QC  ', 'QR  ', &
@@ -173,19 +181,19 @@
 ! Grid parameters - the values will be read from a
 ! standard ncommas namelist and filled in here.
 
-! nx, ny and nz are the size of the dipole (or irregular) grids. 
-integer :: Nx=-1, Ny=-1, Nz=-1    ! grid counts for each field
+! Each spatial dimension has a staggered counterpart.
+integer :: nxc=-1, nyc=-1, nzc=-1    ! scalar grid positions
+integer :: nxe=-1, nye=-1, nze=-1    ! staggered grid positions
 
-! locations of cell centers (C) and edges (G) for each axis.
-real(r8), allocatable :: ZC(:), ZG(:)
+! locations of cell centers (C) and edges (E) for each axis.
+real(r8), allocatable :: ZC(:), ZE(:)
 
 ! These arrays store the longitude and latitude of the lower left corner of
 ! each of the dipole u quadrilaterals and t quadrilaterals.
-real(r8), allocatable :: ULAT(:,:), ULON(:,:), TLAT(:,:), TLON(:,:)
+real(r8), allocatable :: ULAT(:,:), ULON(:,:)
+real(r8), allocatable :: VLAT(:,:), VLON(:,:)
+real(r8), allocatable :: WLAT(:,:), WLON(:,:)
 
-! integer, lowest valid cell number in the vertical
-integer, allocatable  :: KMT(:, :), KMU(:, :)
-
 real(r8)        :: endTime
 real(r8)        :: ocean_dynamics_timestep = 900.0_r4
 integer         :: timestepcount = 0
@@ -253,27 +261,17 @@
 subroutine static_init_model()
 !------------------------------------------------------------------
 !
-! Called to do one time initialization of the model. In this case,
-! it reads in the grid information.
+! Called to do one time initialization of the model.
+! Harvest a ton of information from the NCOMMAS restart file
+! about grid sizes, grid contents, variable sizes, etc..
 
-character(len=32):: calendar
-integer :: iunit, io
+! Local variables - all the important ones have module scope
+integer, dimension(NF90_MAX_VAR_DIMS) :: dimIDs
+character(len=NF90_MAX_NAME)          :: varname
+integer :: ncid, VarID, numdims, dimlen, varsize
+integer :: iunit, io, ivar, i
 integer :: ss, dd
 
-! The Plan:
-!
-!   read in the grid sizes from the horiz grid file and the vert grid file
-!   horiz is netcdf, vert is ascii
-!  
-!   allocate space, and read in actual grid values
-!
-!   figure out model timestep.  FIXME: from where?
-!
-!   Compute the model size.
-!
-!   set the index numbers where the field types change
-!
-
 if ( module_initialized ) return ! only need to do this once.
 
 ! Print module information to log file and stdout.
@@ -297,60 +295,123 @@
 ! Set the time step ... causes ncommas namelists to be read.
 ! Ensures model_timestep is multiple of 'ocean_dynamics_timestep'
 
+call set_calendar_type( calendar )   ! comes from model_mod_nml
+
 model_timestep = set_model_time_step()
 
 call get_time(model_timestep,ss,dd) ! set_time() assures the seconds [0,86400)
 
-write(msgstring,*)'assimilation period is ',dd,' days ',ss,' seconds'
-call error_handler(E_MSG,'static_init_model',msgstring,source,revision,revdate)
+write(string1,*)'assimilation period is ',dd,' days ',ss,' seconds'
+call error_handler(E_MSG,'static_init_model',string1,source,revision,revdate)
 
+call get_ncommas_restart_filename( ncommas_filename )
+
 !---------------------------------------------------------------
-! get data dimensions, then allocate space, then open the files
-! and actually fill in the arrays.
+! 1) get grid dimensions
+! 2) allocate space for the grids 
+! 3) read them, convert them from X-Y-Z to lat-lon-z
 
-call get_horiz_grid_dims(Nx, Ny)
-call get_vert_grid_dim(Nz)
+! call get_grid_dims(nxc, nyc, nzc, )
+call get_grid_dims(nxc, nxe, nyc, nye, nzc, nze )
 
 ! Allocate space for grid variables. 
-allocate(ULAT(Nx,Ny), ULON(Nx,Ny), TLAT(Nx,Ny), TLON(Nx,Ny))
-allocate( KMT(Nx,Ny),  KMU(Nx,Ny))
-allocate(     ZC(Nz),      ZG(Nz))
+allocate(ULAT(nxe,nyc), ULON(nxe,nyc))
+allocate(VLAT(nxc,nye), VLON(nxc,nye))
+allocate(WLAT(nxc,nyc), WLON(nxc,nyc))
+allocate(  ZC(  nzc  ),   ZE(  nze  ))
 
 ! Fill them in.
-! horiz grid initializes ULAT/LON, TLAT/LON as well.
+! horiz grid initializes ULAT/LON, slat/LON as well.
 ! kmt initializes HT/HU if present in input file.
-call read_horiz_grid(Nx, Ny, ULAT, ULON, TLAT, TLON)
-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
+call get_grid(nxc, nxe, nyc, nye, nzc, nze, &
+              ULAT, ULON, VLAT, VLON, WLAT, WLON, ZC, ZE)
 
 !---------------------------------------------------------------
 ! compute the offsets into the state vector for the start of each
-! different variable type.
+! different variable type. Requires reading shapes from the NCOMMAS
+! restart file.
 
-! record where in the state vector the data type changes
+! Record where in the state vector the data type changes
 ! from one type to another, by computing the starting
 ! index for each block of data.
-start_index(S_index)     = 1
-start_index(T_index)     = start_index(S_index) + (Nx * Ny * Nz)
-start_index(U_index)     = start_index(T_index) + (Nx * Ny * Nz)
-start_index(V_index)     = start_index(U_index) + (Nx * Ny * Nz)
-start_index(PSURF_index) = start_index(V_index) + (Nx * Ny * Nz)
 
+call nc_check( nf90_open(trim(ncommas_filename), NF90_NOWRITE, ncid), &
+                  'static_init_model', 'open '//trim(ncommas_filename))
+
+model_size = 0;
+do ivar = 1, nfields 
+
+   varname = adjustl(progvarnames(ivar))
+   string2 = trim(ncommas_filename)//' '//trim(varname)
+
+   progvar(ivar)%varname = trim(varname)
+   progvar(ivar)%dimlens = 0
+
+   call nc_check(nf90_inq_varid(ncid, trim(varname), VarID), &
+            'static_init_model', 'inq_varid '//trim(string2))
+
+   call nc_check( nf90_get_att(ncid, VarId, 'long_name' , progvar(ivar)%long_name), &
+            'static_init_model', 'get_att long_name '//trim(string2))
+
+   call nc_check( nf90_get_att(ncid, VarId, 'units' , progvar(ivar)%units), &
+            'static_init_model', 'get_att units '//trim(string2))
+
+   call nc_check(nf90_inquire_variable(ncid, VarId, dimids=dimIDs, ndims=numdims), &
+            'static_init_model', 'inquire '//trim(string2))
+
+   progvar(ivar)%numdims = numdims
+
+   varsize = 1
+   do i = 1,numdims
+      write(string1,'(''inquire dimension'',i2,A)') i,trim(string2)
+      call nc_check(nf90_inquire_dimension(ncid, dimIDs(i), len=dimlen), 'static_init_model', string1)
+      progvar(ivar)%dimlens(i) = dimlen
+      varsize = varsize * dimlen
+   enddo
+
+   progvar(ivar)%varsize = varsize
+   progvar(ivar)%indexN  = model_size + varsize
+
+   if (do_output()) then
+      write(logfileunit,*) ivar,trim(progvar(ivar)%varname)
+      write(logfileunit,*) '  long_name ',trim(progvar(ivar)%long_name)
+      write(logfileunit,*) '  units     ',trim(progvar(ivar)%units)
+      write(logfileunit,*) '  numdims   ',progvar(ivar)%numdims
+      write(logfileunit,*) '  dimlens   ',progvar(ivar)%dimlens(1:progvar(ivar)%numdims)
+      write(logfileunit,*) '  varsize   ',progvar(ivar)%varsize
+      write(logfileunit,*) '  indexN    ',progvar(ivar)%indexN
+
+      write(     *     ,*) ivar,trim(progvar(ivar)%varname)
+      write(     *     ,*) '  long_name ',trim(progvar(ivar)%long_name)

@@ Diff output truncated at 40000 characters. @@


More information about the Dart-dev mailing list