[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