[Dart-dev] [4020] DART/trunk/models/POP: The grid information is not in the LANL-POP restart file, so it must be
nancy at ucar.edu
nancy at ucar.edu
Thu Aug 27 17:19:21 MDT 2009
An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/dart-dev/attachments/20090827/20731766/attachment-0001.html
-------------- next part --------------
Modified: DART/trunk/models/POP/README
===================================================================
--- DART/trunk/models/POP/README 2009-08-27 20:19:52 UTC (rev 4019)
+++ DART/trunk/models/POP/README 2009-08-27 23:19:20 UTC (rev 4020)
@@ -6,6 +6,8 @@
Tim
+ ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
(this should stay in the README. at some point we might
want to move the test program to a separate test directory,
in which case this file needs to be updated to reflect that.)
@@ -18,3 +20,39 @@
http://www.image.ucar.edu/pub/DART/POP/
nancy
+
+ ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+Tim : Tue Jul 21 17:45:52 MDT 2009
+
+Working on reading the BINARY grid files instead of the (nonexistent)
+netCDF ones. Getting grid sizes from restart netCDF file and must make
+hard assumptions about variable storage order.
+
+dart_pop_mod is being modified to read the pop_in namelist and then set things
+like the ocean dynamics timestep. Must ensure that the model_mod then uses that
+to set a valid adv_to_time ...
+
+dart_pop_mod must also set the time_manager_nml:stop_count (and stop_option) to
+the right values.
+
+dart_pop_mod must also create a pointer file with the expected restart file as
+a way to make sure the model has advanced to the proper spot.
+
+advance_model.csh must paste the pop_in.DART namelist on top of the pop_in namelist
+for model control ...
+
+ ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+This is how I understand (LANL) POP to work:
+
+Given:
+&init_ts_nml
+ init_ts_option = 'restart'
+ init_ts_file = 'pop.r'
+ init_ts_file_fmt = 'nc'
+/
+
+The init_ts_file entry is completely ignored. A pointer file with the name
+"pop_pointer.restart" contains the name of the restart file.
+
Modified: DART/trunk/models/POP/dart_pop_mod.f90
===================================================================
--- DART/trunk/models/POP/dart_pop_mod.f90 2009-08-27 20:19:52 UTC (rev 4019)
+++ DART/trunk/models/POP/dart_pop_mod.f90 2009-08-27 23:19:20 UTC (rev 4020)
@@ -11,307 +11,629 @@
! $Revision$
! $Date$
-use types_mod, only : r8, rad2deg, PI
-use obs_def_mod, only : obs_def_type, get_obs_def_time, read_obs_def, &
- write_obs_def, destroy_obs_def, interactive_obs_def, &
- copy_obs_def, set_obs_def_time, set_obs_def_kind, &
- set_obs_def_error_variance, set_obs_def_location
-use time_manager_mod, only : time_type, get_date, set_time, GREGORIAN, &
- set_date, set_calendar_type, get_time, &
- print_date, print_time, operator(==)
+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, &
- E_ERR, E_MSG, timestamp
-use location_mod, only : location_type, set_location, VERTISHEIGHT, VERTISSURFACE
-use obs_sequence_mod, only : init_obs_sequence, init_obs, insert_obs_in_seq, &
- set_obs_values, set_qc, obs_sequence_type, obs_type, &
- copy_obs, set_copy_meta_data, set_qc_meta_data, set_obs_def, &
- get_first_obs, get_last_obs, get_obs_def
+ register_module, error_handler, nc_check, &
+ find_namelist_in_file, check_namelist_read, &
+ E_ERR, E_MSG, timestamp, find_textfile_dims, &
+ logfileunit
-use obs_kind_mod, only : get_obs_kind_index
+use typesizes
+use netcdf
implicit none
private
-public :: real_obs_sequence
+public :: get_pop_calendar, set_model_time_step, &
+ get_horiz_grid_dims, get_vert_grid_dim, &
+ read_horiz_grid, read_topography, read_vert_grid, &
+ write_pop_namelist
! version controlled file description for error handling, do not edit
character(len=128), parameter :: &
- source = "$URL$", &
- revision = "$Revision$", &
- revdate = "$Date$"
+ source = '$URL$', &
+ revision = '$Revision$', &
+ revdate = '$Date$'
+character(len=256) :: msgstring
logical, save :: module_initialized = .false.
+character(len=256) :: ic_filename, restart_filename
+
! 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
-contains
+!------------------------------------------------------------------
+! The POP time manager namelist variables
+!------------------------------------------------------------------
-!-------------------------------------------------
+character(len=100) :: accel_file ! length consistent with POP
+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
-function real_obs_sequence (obsfile, year, month, day, max_num, &
- lon1, lon2, lat1, lat2)
-!------------------------------------------------------------------------------
-! this function is to prepare data to DART sequence format
+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 POP restart manager namelist variables
+!------------------------------------------------------------------
+
+character(len=100) :: restart_outfile ! length consistent with POP
+character(len= 64) :: restart_freq_opt, restart_fmt
+logical :: leven_odd_on, pressure_correction
+integer :: restart_freq, even_odd_freq
+
+namelist /restart_nml/ restart_freq_opt, restart_freq, restart_outfile, &
+ restart_fmt, leven_odd_on, even_odd_freq, pressure_correction
+
+!------------------------------------------------------------------
+! The POP initial temperature and salinity namelist
+!------------------------------------------------------------------
+
+character(len=100) :: init_ts_file ! length consistent with POP
+character(len= 64) :: init_ts_option, init_ts_file_fmt
+
+namelist /init_ts_nml/ init_ts_option, init_ts_file, init_ts_file_fmt
+
+!------------------------------------------------------------------
+! The POP domain namelist
+!------------------------------------------------------------------
+
+character(len= 64) :: clinic_distribution_type, tropic_distribution_type
+character(len= 64) :: ew_boundary_type, ns_boundary_type
+integer :: nprocs_clinic, nprocs_tropic
+
+namelist /domain_nml/ clinic_distribution_type, nprocs_clinic, &
+ tropic_distribution_type, nprocs_tropic, &
+ ew_boundary_type, ns_boundary_type
+
+!------------------------------------------------------------------
+! The POP grid info namelist
+!------------------------------------------------------------------
!
-character(len=129), intent(in) :: obsfile
-integer, intent(in) :: year, month, day, max_num
-real(r8), intent(in) :: lon1, lon2, lat1, lat2
+! POP 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.
+!
+!------------------------------------------------------------------
+!
+! 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)
+!
+!------------------------------------------------------------------
-type(obs_sequence_type) :: real_obs_sequence
+character(len=100) :: horiz_grid_file, vert_grid_file, topography_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
+namelist /grid_nml/ horiz_grid_opt, horiz_grid_file, sfc_layer_opt, &
+ vert_grid_opt, vert_grid_file, topography_opt, topography_file, &
+ partial_bottom_cells, bottom_cell_file, region_mask_file, &
+ topo_smooth, flat_bottom, lremove_points
-type(obs_def_type) :: obs_def
-type(obs_type) :: obs, prev_obs
-integer :: i, num_copies, num_qc
-integer :: days, seconds
-integer :: yy, mn, dd, hh, mm, ss
-integer :: startdate1, startdate2
-integer :: obs_num, calender_type, iskip
-integer :: obs_unit
-integer :: which_vert, obstype
+!======================================================================
+contains
+!======================================================================
-real (r8) :: lon, lat, vloc, obs_value
-real (r8) :: aqc, var2, lonc
-type(time_type) :: time, pre_time
-character(len = 32) :: obs_kind_name
-character(len = 80) :: label
-character(len = 129) :: copy_meta_data, qc_meta_data
+subroutine initialize_module
+!------------------------------------------------------------------
+integer :: iunit, io
-if ( .not. module_initialized ) call initialize_module
+! Read POP 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.
-num_copies = 1
-num_qc = 1
+call find_namelist_in_file('pop_in', 'time_manager_nml', iunit)
+read(iunit, nml = time_manager_nml, iostat = io)
+call check_namelist_read(iunit, io, 'time_manager_nml')
-! Initialize an obs_sequence
+if ( allow_leapyear ) then
+ call set_calendar_type('gregorian')
+else
+ call set_calendar_type('noleap')
+endif
-call init_obs_sequence(real_obs_sequence, num_copies, num_qc, max_num)
+! Read POP initial information (for input/restart filename)
+! The tricky part here is that we should really check for
+! the existence of the init_ts_file and take evasive action
+! like checking for the existence of a pointer file.
-! set meta data of obs_seq
+call find_namelist_in_file('pop_in', 'init_ts_nml', iunit)
+read(iunit, nml = init_ts_nml, iostat = io)
+call check_namelist_read(iunit, io, 'init_ts_nml')
-do i = 1, num_copies
- copy_meta_data = 'observation'
- call set_copy_meta_data(real_obs_sequence, i, copy_meta_data)
-end do
+ic_filename = trim(init_ts_file)//'.'//trim(init_ts_file_fmt)
-do i = 1, num_qc
- qc_meta_data = 'QC index'
- call set_qc_meta_data(real_obs_sequence, i, qc_meta_data)
-end do
+! FIXME ... what about the pointer file ...
+if ( .not. file_exist(ic_filename) ) then
+ msgstring = 'pop_in:init_ts_file '//trim(ic_filename)//' not found'
+ call error_handler(E_ERR,'initialize_module', &
+ msgstring, source, revision, revdate)
+endif
-! Initialize the obs variable
+! Read POP restart information (for model timestepping/grid dimensions)
+call find_namelist_in_file('pop_in', 'restart_nml', iunit)
+read(iunit, nml = restart_nml, iostat = io)
+call check_namelist_read(iunit, io, 'restart_nml')
-call init_obs(obs, num_copies, num_qc)
-call init_obs(prev_obs, num_copies, num_qc)
+! Read POP domain information (for lon wrapping or not)
+call find_namelist_in_file('pop_in', 'domain_nml', iunit)
+read(iunit, nml = domain_nml, iostat = io)
+call check_namelist_read(iunit, io, 'domain_nml')
-! set observation time type
-calender_type = GREGORIAN
-call set_calendar_type(calender_type)
+! Read POP grid information (for grid dims/filenames)
+call find_namelist_in_file('pop_in', 'grid_nml', iunit)
+read(iunit, nml = grid_nml, iostat = io)
+call check_namelist_read(iunit, io, 'grid_nml')
-! open observation data file
+module_initialized = .true.
-obs_unit = get_unit()
-open(unit = obs_unit, file = obsfile, form='formatted', status='old')
-print*, 'input file opened= ', trim(obsfile)
-rewind (obs_unit)
+! Print module information to log file and stdout.
+call register_module(source, revision, revdate)
-obs_num = 0
-iskip = 0
+end subroutine initialize_module
-! loop over all observations within the file
-!------------------------------------------------------------------------------
-obsloop: do
- read(obs_unit,*,end=200) lon, lat, vloc, obs_value, which_vert, var2, aqc, &
- obs_kind_name, startdate1, startdate2
+subroutine get_horiz_grid_dims(Nx, Ny)
+!------------------------------------------------------------------
+! 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.
+!
+! The file name comes from module storage ... namelist.
- !print*,''
- !print*,' Observation ', obs_num+1
- !print*,' lon lat vloc obs_value ',lon, lat, vloc, obs_value
- !print*,' which_vert var2 aqc ',which_vert, var2, aqc
- !print*,' obs_kind_name ',obs_kind_name
- !print*,' date1 date2 ',startdate1, startdate2
+integer, intent(out) :: Nx ! Number of Longitudes
+integer, intent(out) :: Ny ! Number of Latitudes
- ! Calculate the DART time from the observation time
- yy = startdate1/10000
- mn = mod(startdate1/100,100)
- dd = mod(startdate1 ,100)
- hh = startdate2/10000
- mm = mod(startdate2/100,100)
- ss = mod(startdate2 ,100)
- time = set_date(yy,mn,dd,hh,mm,ss)
- call get_time(time,seconds,days)
+integer :: grid_id, dimid, nc_rc
- ! verify the location is not outside valid limits
- if((lon > 360.0_r8) .or. (lon < 0.0_r8) .or. &
- (lat > 90.0_r8) .or. (lat < -90.0_r8)) then
- write(*,*) 'invalid location. lon,lat = ', lon, lat
- iskip = iskip + 1
- cycle obsloop
- endif
+if ( .not. module_initialized ) call initialize_module
- lonc = lon
- if (lon2 > 360.0_r8 .and. lon < 180.0_r8) lonc = lon + 360.0_r8
+! get the ball rolling ...
- ! reject observations outside the bounding box
- if(lat < lat1 .or. lat > lat2 .or. lonc < lon1 .or. lonc > lon2) then
- iskip = iskip + 1
- cycle obsloop
+call nc_check(nf90_open(trim(ic_filename), nf90_nowrite, grid_id), &
+ 'get_horiz_grid_dims','open '//trim(ic_filename))
+
+! 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
- ! assign each observation the correct observation type
- obstype = get_obs_kind_index(obs_kind_name)
- if(obstype < 1) then
- print*, 'unknown observation type [',trim(obs_kind_name),'] ... skipping ...'
- cycle obsloop
- !else
- ! print*,trim(obs_kind_name),' is ',obstype
+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
- obs_num = obs_num + 1
+call nc_check(nf90_inquire_dimension(grid_id, dimid, len=Ny), &
+ 'get_horiz_grid_dims','inquire_dimension i '//trim(ic_filename))
- ! print a reassuring message after every Nth processed obs.
- ! if requested, print in the form of a timestamp.
- ! the default is just a plain string with the current obs count.
- if(mod(obs_num, print_every_Nth) == 0) then
- write(label, *) 'obs count = ', obs_num
- if (print_timestamps) then
- call timestamp(string1=label, pos='')
- else
- write(*,*) trim(label)
- endif
- endif
- if(obs_num == max_num) then
- print*, 'Max limit for observation count reached. Increase value in namelist'
- stop
- endif
+! tidy up
+
+call nc_check(nf90_close(grid_id), &
+ 'get_horiz_grid_dims','close '//trim(ic_filename) )
+
+end subroutine get_horiz_grid_dims
+
+
+
+ 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.
+
+integer, intent(out) :: Nz
+
+integer :: linelen ! disposable
+
+if ( .not. module_initialized ) call initialize_module
+
+call find_textfile_dims(vert_grid_file, Nz, linelen)
+
+end subroutine get_vert_grid_dim
+
+
-! create the obs_def for this observation, add to sequence
-!------------------------------------------------------------------------------
-
- call real_obs(num_copies, num_qc, obs, lon, lat, vloc, obs_value, &
- var2, aqc, obstype, which_vert, seconds, days)
-
- if(obs_num == 1) then ! for the first observation
+subroutine get_pop_calendar(calstring)
+!------------------------------------------------------------------
+! the initialize_module ensures that the pop namelists are read and
+! the DART time manager gets the pop calendar setting.
+!
+! Then, the DART time manager is queried to return what it knows ...
+!
+character(len=*), INTENT(OUT) :: calstring
- call insert_obs_in_seq(real_obs_sequence, obs)
- call copy_obs(prev_obs, obs)
- pre_time = time
+if ( .not. module_initialized ) call initialize_module
- else ! not the first observation
+call get_calendar_string(calstring)
- if(time == pre_time) then ! same time as previous observation
+end subroutine get_pop_calendar
- call insert_obs_in_seq(real_obs_sequence, obs, prev_obs)
- call copy_obs(prev_obs, obs)
- pre_time = time
- else ! not the same time
- call insert_obs_in_seq(real_obs_sequence, obs)
- call copy_obs(prev_obs, obs)
- pre_time = time
+function set_model_time_step()
+!------------------------------------------------------------------
+! the initialize_module ensures that the pop namelists are read.
+! The restart times in the pop_in&restart_nml are used to define
+! appropriate assimilation timesteps.
+!
+type(time_type) :: set_model_time_step
- endif
+if ( .not. module_initialized ) call initialize_module
- endif
+! Check the 'restart_freq_opt' and 'restart_freq' to determine
+! when we can stop the model
-end do obsloop
+if ( trim(restart_freq_opt) == 'nday' ) then
+ set_model_time_step = set_time(0, restart_freq) ! (seconds, days)
+else
+ call error_handler(E_ERR,'set_model_time_step', &
+ 'restart_freq_opt must be days', source, revision, revdate)
+endif
-200 continue
+end function set_model_time_step
-close(obs_unit)
-! Print a little summary
-print*, 'obs used = ', obs_num, ' obs skipped = ', iskip
-if ( get_first_obs(real_obs_sequence, obs) ) then
- call get_obs_def(obs, obs_def)
- pre_time = get_obs_def_time(obs_def)
- call print_time(pre_time,' first time in sequence is ')
- call print_date(pre_time,' first date in sequence is ')
+
+subroutine write_pop_namelist(model_time, adv_to_time)
+!------------------------------------------------------------------
+!
+type(time_type), INTENT(IN) :: model_time, adv_to_time
+type(time_type) :: offset
+
+integer :: iunit, secs, days
+
+if ( .not. module_initialized ) call initialize_module
+
+offset = adv_to_time - model_time
+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_pop_namelist', msgstring, source, revision, revdate)
endif
-if( get_last_obs(real_obs_sequence, obs)) then
- call get_obs_def(obs, obs_def)
- time = get_obs_def_time(obs_def)
- call print_time(time,' last time in sequence is ')
- call print_date(time,' last date in sequence is ')
+
+! call print_date( model_time,'write_pop_namelist:dart model date')
+! call print_date(adv_to_time,'write_pop_namelist:advance_to date')
+! call print_time( model_time,'write_pop_namelist:dart model time')
+! call print_time(adv_to_time,'write_pop_namelist:advance_to time')
+! call print_time( offset,'write_pop_namelist:a distance of')
+! write( *,'(''write_pop_namelist:TIME_MANAGER_NML STOP_COUNT '',i10,'' days'')') days
+
+!Convey the information to the namelist 'stop option' and 'stop count'
+
+if ( trim(stop_option) == 'nday' ) then
+ stop_count = days
+else
+ call error_handler(E_ERR,'write_pop_namelist', &
+ 'stop_option must be "nday"', source, revision, revdate)
endif
-print*, ''
-end function real_obs_sequence
+iunit = open_file('pop_in.DART',form='formatted',action='rewind')
+write(iunit, nml=time_manager_nml)
+write(iunit, '('' '')')
+close(iunit)
+end subroutine write_pop_namelist
-subroutine real_obs(num_copies, num_qc, obs, lon, lat, vloc, obs_value, &
- var2, aqc, obs_kind, which_vert, seconds, days)
-!------------------------------------------------------------------------------
-integer, intent(in) :: num_copies, num_qc
-type(obs_type), intent(inout) :: obs
-real(r8), intent(in) :: lon, lat, vloc, obs_value, var2, aqc
-integer, intent(in) :: obs_kind, which_vert, seconds, days
-integer :: i
-real(r8) :: aqc01(1), obs_value01(1)
-type(obs_def_type) :: obsdef0
+ 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
-! Does real initialization of an observation type
+! Check to see that the file exists.
-call real_obs_def(obsdef0, lon, lat, vloc, &
- var2, obs_kind, which_vert, seconds, days)
-call set_obs_def(obs, obsdef0)
+if ( .not. file_exist(horiz_grid_file) ) then
+ msgstring = 'pop_in:horiz_grid_file '//trim(horiz_grid_file)//' not found'
+ call error_handler(E_ERR,'read_horiz_grid', &
+ msgstring, source, revision, revdate)
+endif
-do i = 1, num_copies
- obs_value01(1) = obs_value
- call set_obs_values(obs, obs_value01(1:1) )
+! 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' )
+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)
+!
+! mimic POP grid.F90:calc_tpoints(), but for one big block.
+
+integer, intent( in) :: nx, ny
+real(r8), dimension(nx,ny), intent( in) :: ULAT, ULON
+real(r8), dimension(nx,ny), intent(out) :: TLAT, TLON
+
+integer :: i, j
+real(r8) :: xc,yc,zc,xs,ys,zs,xw,yw,zw ! Cartesian coordinates for
+real(r8) :: xsw,ysw,zsw,tx,ty,tz,da ! nbr points
+
+real(r8), parameter :: c0 = 0.000_r8, c1 = 1.000_r8
+real(r8), parameter :: c2 = 2.000_r8, c4 = 4.000_r8
+real(r8), parameter :: p25 = 0.250_r8, p5 = 0.500_r8
+real(r8) :: pi, pi2, pih, radian
+
+if ( .not. module_initialized ) call initialize_module
+
+! Define some constants as in pop
+
+pi = c4*atan(c1)
+pi2 = c2*pi
+pih = p5*pi
+radian = 180.0_r8/pi
+
+do j=2,ny
+do i=2,nx
+
+ !*** convert neighbor U-cell coordinates to 3-d Cartesian coordinates
+ !*** to prevent problems with averaging near the pole
+
+ zsw = cos(ULAT(i-1,j-1))
+ xsw = cos(ULON(i-1,j-1))*zsw
+ ysw = sin(ULON(i-1,j-1))*zsw
+ zsw = sin(ULAT(i-1,j-1))
+
+ zs = cos(ULAT(i ,j-1))
+ xs = cos(ULON(i ,j-1))*zs
+ ys = sin(ULON(i ,j-1))*zs
+ zs = sin(ULAT(i ,j-1))
+
+ zw = cos(ULAT(i-1,j ))
+ xw = cos(ULON(i-1,j ))*zw
+ yw = sin(ULON(i-1,j ))*zw
+ zw = sin(ULAT(i-1,j ))
+
+ zc = cos(ULAT(i ,j ))
+ xc = cos(ULON(i ,j ))*zc
+ yc = sin(ULON(i ,j ))*zc
+ zc = sin(ULAT(i ,j ))
+
+ !*** straight 4-point average to T-cell Cartesian coords
+
+ tx = p25*(xc + xs + xw + xsw)
+ ty = p25*(yc + ys + yw + ysw)
+ tz = p25*(zc + zs + zw + zsw)
+
+ !*** convert to lat/lon in radians
+
+ da = sqrt(tx**2 + ty**2 + tz**2)
+
+ TLAT(i,j) = asin(tz/da)
+
+ if (tx /= c0 .or. ty /= c0) then
+ TLON(i,j) = atan2(ty,tx)
+ else
+ TLON(i,j) = c0
+ endif
+
end do
+end do
-do i = 1, num_qc
- aqc01(1) = aqc
- call set_qc(obs, aqc01(1:1))
+!*** for bottom row of domain where sw 4pt average is not valid,
+!*** extrapolate from interior
+!*** NOTE: THIS ASSUMES A CLOSED SOUTH BOUNDARY - WILL NOT
+!*** WORK CORRECTLY FOR CYCLIC OPTION
+
+do i=1,nx
+ TLON(i,1) = TLON(i,1+1)
+ TLAT(i,1) = c2*TLAT(i,1+1) - TLAT(i,1+2)
end do
-end subroutine real_obs
+where (TLON(:,:) > pi2) TLON(:,:) = TLON(:,:) - pi2
+where (TLON(:,:) < c0 ) TLON(:,:) = TLON(:,:) + pi2
+!*** this leaves the leftmost/western edge to be filled
+!*** if the longitudes wrap, this is easy.
+!*** the gx3v5 grid TLON(:,2) and TLON(:,nx) are both about 2pi,
+!*** so taking the average is reasonable.
+!*** averaging the latitudes is always reasonable.
+if ( trim(ew_boundary_type) == 'cyclic' ) then
-subroutine real_obs_def(obs_def, lon, lat, vloc, &
- var2, obs_kind, which_vert, seconds, days)
-!----------------------------------------------------------------------
-type(obs_def_type), intent(inout) :: obs_def
-real(r8),intent(in) :: lon, lat, vloc, var2
-integer, intent(in) :: obs_kind, which_vert, seconds, days
+ TLAT(1,:) = (TLAT(2,:) + TLAT(nx,:))/c2
+ TLON(1,:) = (TLON(2,:) + TLON(nx,:))/c2
-type(location_type) :: loc0
+else
+ write(msgstring,'(''pop_in&domain_nml:ew_boundary_type '',a,'' unknown.'')') &
+ trim(ew_boundary_type)
+ call error_handler(E_ERR,'calc_tpoints',msgstring,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
-! set obs location
-loc0 = set_location(lon, lat, vloc, which_vert )
-call set_obs_def_location(obs_def, loc0)
+! Check to see that the file exists.
-! set obs kind
-call set_obs_def_kind(obs_def, obs_kind)
+if ( .not. file_exist(topography_file) ) then
+ msgstring = 'pop_in:topography_file '//trim(topography_file)//' not found'
+ call error_handler(E_ERR,'read_topography', &
+ msgstring, source, revision, revdate)
+endif
-call set_obs_def_time(obs_def, set_time(seconds, days) )
-call set_obs_def_error_variance(obs_def, var2)
+! read the binary file
-end subroutine real_obs_def
+topo_unit = get_unit()
+INQUIRE(iolength=reclength) KMT
+open( topo_unit, file=trim(topography_file), form='unformatted', &
+ access='direct', recl=reclength, status='old' )
+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
-subroutine initialize_module
-!-------------------------------------------------
-call register_module(source, revision, revdate)
-module_initialized = .true.
-end subroutine initialize_module
+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 = 'pop_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
+
+
+
end module dart_pop_mod
Modified: DART/trunk/models/POP/dart_to_pop.f90
===================================================================
--- DART/trunk/models/POP/dart_to_pop.f90 2009-08-27 20:19:52 UTC (rev 4019)
+++ DART/trunk/models/POP/dart_to_pop.f90 2009-08-27 23:19:20 UTC (rev 4020)
@@ -27,12 +27,13 @@
initialize_utilities, finalize_utilities, &
find_namelist_in_file, check_namelist_read, &
logfileunit
+use assim_model_mod, only : open_restart_read, aread_state_restart, close_restart
+use time_manager_mod, only : time_type, get_time, print_time, print_date, &
+ operator(-), set_time
use model_mod, only : static_init_model, sv_to_restart_file, &
get_model_size, get_model_time_step, &
set_model_end_time
-use assim_model_mod, only : open_restart_read, aread_state_restart, close_restart
-use time_manager_mod, only : time_type, get_time, print_time, print_date, &
- operator(-)
+use dart_pop_mod, only : write_pop_namelist
implicit none
@@ -48,15 +49,16 @@
character (len = 128) :: dart_to_pop_input_file = 'assim_model_state_ic'
character (len = 128) :: dart_to_pop_restart_file = 'my_pop_restart_file'
+logical :: advance_time_present = .TRUE.
-namelist /dart_to_pop_nml/ dart_to_pop_input_file, dart_to_pop_restart_file
+namelist /dart_to_pop_nml/ dart_to_pop_input_file, dart_to_pop_restart_file, &
+ advance_time_present
!----------------------------------------------------------------------
integer :: iunit, io, x_size
integer :: secs, days
type(time_type) :: model_time, adv_to_time
-type(time_type) :: model_timestep, offset
real(r8), allocatable :: statevector(:)
!----------------------------------------------------------------------
@@ -82,52 +84,41 @@
!----------------------------------------------------------------------
iunit = open_restart_read(dart_to_pop_input_file)
-call aread_state_restart(model_time, statevector, iunit, adv_to_time)
+
+if ( advance_time_present ) then
+ call aread_state_restart(model_time, statevector, iunit, adv_to_time)
+else
+ call aread_state_restart(model_time, statevector, iunit)
+endif
call close_restart(iunit)
!----------------------------------------------------------------------
! update the current POP state vector
-!----------------------------------------------------------------------
-
-call sv_to_restart_file(statevector, dart_to_pop_restart_file, &
- model_time, adv_to_time)
-
-!iunit = open_file('data.cal.DART',form='formatted',action='rewind')
-!write(iunit, nml=CAL_NML)
-!close(iunit)
-
-!----------------------------------------------------------------------
-! convert the adv_to_time to the appropriate number of POP
+! Convey the amount of time to integrate the model ...
! time_manager_nml: stop_option, stop_count increments
!----------------------------------------------------------------------
-model_timestep = get_model_time_step()
-offset = adv_to_time - model_time
+call sv_to_restart_file(statevector, dart_to_pop_restart_file, model_time)
-!call set_model_end_time(offset)
-!call write_data_namelistfile()
+if ( advance_time_present ) then
+ call write_pop_namelist(model_time, adv_to_time)
+endif
!----------------------------------------------------------------------
! Log what we think we're doing, and exit.
!----------------------------------------------------------------------
-call get_time(offset, secs, days)
-
call print_date( model_time,'dart_to_pop:dart model date')
-call print_date(adv_to_time,'dart_to_pop:advance_to date')
call print_time( model_time,'dart_to_pop:dart model time')
-call print_time(adv_to_time,'dart_to_pop:advance_to time')
-call print_time( offset,'dart_to_pop:a distance of')
-write( * ,'(''dart_to_pop:PARM03 endTime '',i,'' seconds'')') &
- (secs + days*SECPERDAY)
-
call print_date( model_time,'dart_to_pop:dart model date',logfileunit)
-call print_date(adv_to_time,'dart_to_pop:advance_to date',logfileunit)
call print_time( model_time,'dart_to_pop:dart model time',logfileunit)
+
+if ( advance_time_present ) then
+call print_time(adv_to_time,'dart_to_pop:advance_to time')
+call print_date(adv_to_time,'dart_to_pop:advance_to date')
call print_time(adv_to_time,'dart_to_pop:advance_to time',logfileunit)
-call print_time( offset,'dart_to_pop: a distance of',logfileunit)
-write(logfileunit,'(''dart_to_pop:PARM03 endTime '',i,'' seconds'')') &
- (secs + days*SECPERDAY)
+call print_date(adv_to_time,'dart_to_pop:advance_to date',logfileunit)
+endif
call finalize_utilities()
Modified: DART/trunk/models/POP/matlab/Check_pop_to_dart.m
===================================================================
--- DART/trunk/models/POP/matlab/Check_pop_to_dart.m 2009-08-27 20:19:52 UTC (rev 4019)
+++ DART/trunk/models/POP/matlab/Check_pop_to_dart.m 2009-08-27 23:19:20 UTC (rev 4020)
@@ -1,166 +1,160 @@
-% Check_pop_to_dart
+function [dart pop] = Check_pop_to_dart(popfile,dartfile)
+% Check_pop_to_dart : check pop_to_dart.f90 ... the conversion of a POP restart to a DART state vector file.
%
-% data.cal holds the starting time information
+% popfile = 'pop.r.nc';
+% dartfile = 'dart.ics';
+% x = Check_pop_to_dart(popfile, dartfile);
%
+% popfile = '~DART/models/POP/work/cx3.dart.001.pop.r.0002-01-01-00000.nc';
+% dartfile = '~DART/models/POP/work/perfect_ics';
+% [dart pop] = Check_pop_to_dart(popfile, dartfile);
-popdir = '/ptmp/thoar/POP/job_40705';
-popfile = 'pop.r.x1A.00000102';
-dartfile = 'assim_model_state_ud';
+% Data Assimilation Research Testbed -- DART
+% Copyright 2004-2007, Data Assimilation Research Section
+% University Corporation for Atmospheric Research
+% Licensed under the GPL -- www.gpl.org/licenses/gpl.html
+%
+% <next few lines under version control, do not edit>
+% $URL$
+% $Id$
+% $Revision$
+% $Date$
-fname = sprintf('%s/%s',popdir,popfile);
+% Read the original POP file values.
+if (exist(popfile,'file') ~= 2)
+ error('POP file %s does not exist.',popfile)
+end
+if (exist(dartfile,'file') ~= 2)
+ error('DART file %s does not exist.',dartfile)
+end
-S = nc_varget(fname, 'SALT_CUR');
-T = nc_varget(fname, 'TEMP_CUR');
-U = nc_varget(fname, 'UVEL_CUR');
-V = nc_varget(fname, 'VVEL_CUR');
-SSH = nc_varget(fname,'PSURF_CUR');
+iyear = nc_attget(popfile,nc_global,'iyear');
+imonth = nc_attget(popfile,nc_global,'imonth');
+iday = nc_attget(popfile,nc_global,'iday');
+ihour = nc_attget(popfile,nc_global,'ihour');
+iminute = nc_attget(popfile,nc_global,'iminute');
+isecond = nc_attget(popfile,nc_global,'isecond');
-modelsize = prod(size(S)) + prod(size(T)) + ...
- prod(size(U)) + prod(size(V)) + prod(size(SSH));
+fprintf('POP year month day hour minute second %d %d %d %d %d %d\n', ...
+ iyear,imonth,iday,ihour,iminute,isecond);
-[nx ny nz] = size(S);
+% The nc_varget() function returns the variables with the fastest
+% varying dimension on the right. This is opposite to the Fortran
+% convention of the fastest varying dimension on the left ... so
+% one of the variables must be permuted in order to be compared.
-iyear = nc_attget(fname,nc_global,'iyear');
-imonth = nc_attget(fname,nc_global,'imonth');
-iday = nc_attget(fname,nc_global,'iday');
-ihour = nc_attget(fname,nc_global,'ihour');
-iminute = nc_attget(fname,nc_global,'iminute');
-isecond = nc_attget(fname,nc_global,'isecond');
+S = nc_varget(popfile, 'SALT_CUR'); pop.S = permute(S, [3 2 1]);
+T = nc_varget(popfile, 'TEMP_CUR'); pop.T = permute(T, [3 2 1]);
+U = nc_varget(popfile, 'UVEL_CUR'); pop.U = permute(U, [3 2 1]);
+V = nc_varget(popfile, 'VVEL_CUR'); pop.V = permute(V, [3 2 1]);
+PSURF = nc_varget(popfile, 'PSURF_CUR'); pop.PSURF = permute(PSURF, [2 1]);
-fprintf('year month day hour minute second %d %d %d %d %d %d\n', ...
- iyear,imonth,iday,ihour,iminute,isecond);
+disp(sprintf('pop.PSURF min/max are %0.8g %0.8g',min(pop.PSURF(:)),max(pop.PSURF(:))))
-% Get the dart equivalent
+[nx ny nz] = size(pop.U);
+fprintf('vert dimension size is %d\n',nz)
+fprintf('N-S dimension size is %d\n',ny)
+fprintf('E-W dimension size is %d\n',nx)
-fname = sprintf('%s/%s',popdir,dartfile);
-fid = fopen(fname,'rb','ieee-le');
+modelsize = nx*ny*nz;
+
+% filesize = S,T,U,V * (nx*ny*nz) + SSH * (nx*ny)
+storage = 8;
+size2d = nx*ny;
+size3d = nx*ny*nz;
+n2ditems = 1*size2d;
+n3ditems = 4*size3d;
+rec1size = 4+(4+4)+4; % time stamps ...
+rec2size = 4+(n3ditems*storage + n2ditems*storage)+4;
+
+fsize = rec1size + rec2size;
+disp(sprintf('with a modelsize of %d the file size should be %d bytes', ...
+ modelsize,fsize))
+
+% Open and read timetag for state
+fid = fopen(dartfile,'rb','ieee-le');
trec1 = fread(fid,1,'int32');
seconds = fread(fid,1,'int32');
days = fread(fid,1,'int32');
trecN = fread(fid,1,'int32');
@@ Diff output truncated at 40000 characters. @@
More information about the Dart-dev
mailing list