[Dart-dev] [4446] DART/trunk/models: Raw commit -- nothing has been changed to actually reflect anything
nancy at ucar.edu
nancy at ucar.edu
Fri Jul 30 14:30:25 MDT 2010
Revision: 4446
Author: thoar
Date: 2010-07-30 14:30:24 -0600 (Fri, 30 Jul 2010)
Log Message:
-----------
Raw commit -- nothing has been changed to actually reflect anything
from NCOMMAS. Use as templates to guide the division of who works on what
when Lou Whitaker and Ted Mansell - arrive Tues Aug 3rd 2010
Added Paths:
-----------
DART/trunk/models/NCOMMAS/
DART/trunk/models/NCOMMAS/dart_ncommas_mod.f90
DART/trunk/models/NCOMMAS/dart_to_ncommas.f90
DART/trunk/models/NCOMMAS/dart_to_ncommas.nml
DART/trunk/models/NCOMMAS/matlab/
DART/trunk/models/NCOMMAS/model_mod.f90
DART/trunk/models/NCOMMAS/model_mod.nml
DART/trunk/models/NCOMMAS/ncommas_to_dart.f90
DART/trunk/models/NCOMMAS/ncommas_to_dart.nml
DART/trunk/models/NCOMMAS/shell_scripts/
DART/trunk/models/NCOMMAS/shell_scripts/Run_NCOMMAS.lsf
DART/trunk/models/NCOMMAS/shell_scripts/advance_model.csh
DART/trunk/models/NCOMMAS/shell_scripts/job.simple.csh
DART/trunk/models/NCOMMAS/shell_scripts/run_filter.batch
DART/trunk/models/NCOMMAS/shell_scripts/run_perfect_model_obs.batch
DART/trunk/models/NCOMMAS/work/
DART/trunk/models/NCOMMAS/work/mkmf_create_fixed_network_seq
DART/trunk/models/NCOMMAS/work/mkmf_create_obs_sequence
DART/trunk/models/NCOMMAS/work/mkmf_dart_to_ncommas
DART/trunk/models/NCOMMAS/work/mkmf_filter
DART/trunk/models/NCOMMAS/work/mkmf_ncommas_to_dart
DART/trunk/models/NCOMMAS/work/mkmf_obs_diag
DART/trunk/models/NCOMMAS/work/mkmf_obs_seq_to_netcdf
DART/trunk/models/NCOMMAS/work/mkmf_obs_sequence_tool
DART/trunk/models/NCOMMAS/work/mkmf_perfect_model_obs
DART/trunk/models/NCOMMAS/work/mkmf_preprocess
DART/trunk/models/NCOMMAS/work/mkmf_restart_file_tool
DART/trunk/models/NCOMMAS/work/mkmf_wakeup_filter
DART/trunk/models/NCOMMAS/work/path_names_create_fixed_network_seq
DART/trunk/models/NCOMMAS/work/path_names_create_obs_sequence
DART/trunk/models/NCOMMAS/work/path_names_dart_to_ncommas
DART/trunk/models/NCOMMAS/work/path_names_filter
DART/trunk/models/NCOMMAS/work/path_names_ncommas_to_dart
DART/trunk/models/NCOMMAS/work/path_names_obs_diag
DART/trunk/models/NCOMMAS/work/path_names_obs_seq_to_netcdf
DART/trunk/models/NCOMMAS/work/path_names_obs_sequence_tool
DART/trunk/models/NCOMMAS/work/path_names_perfect_model_obs
DART/trunk/models/NCOMMAS/work/path_names_preprocess
DART/trunk/models/NCOMMAS/work/path_names_restart_file_tool
DART/trunk/models/NCOMMAS/work/path_names_wakeup_filter
DART/trunk/models/NCOMMAS/work/quickbuild.csh
-------------- next part --------------
Added: DART/trunk/models/NCOMMAS/dart_ncommas_mod.f90
===================================================================
--- DART/trunk/models/NCOMMAS/dart_ncommas_mod.f90 (rev 0)
+++ DART/trunk/models/NCOMMAS/dart_ncommas_mod.f90 2010-07-30 20:30:24 UTC (rev 4446)
@@ -0,0 +1,701 @@
+! DART software - Copyright \xA9 2004 - 2010 UCAR. This open source software is
+! provided by UCAR, "as is", without charge, subject to all terms of use at
+! http://www.image.ucar.edu/DAReS/DART/DART_download
+
+module dart_ncommas_mod
+
+! <next few lines under version control, do not edit>
+! $URL$
+! $Id$
+! $Revision$
+! $Date$
+
+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, &
+ find_namelist_in_file, check_namelist_read, &
+ E_ERR, E_MSG, timestamp, find_textfile_dims, &
+ logfileunit
+
+use typesizes
+use netcdf
+
+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, &
+ write_ncommas_namelist, get_ncommas_restart_filename
+
+! version controlled file description for error handling, do not edit
+character(len=128), parameter :: &
+ source = '$URL$', &
+ revision = '$Revision$', &
+ revdate = '$Date$'
+
+character(len=256) :: msgstring
+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
+!------------------------------------------------------------------
+
+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
+
+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
+!------------------------------------------------------------------
+
+character(len=100) :: log_filename, pointer_filename
+logical :: lredirect_stdout, luse_pointer_files
+logical :: luse_nf_64bit_offset
+integer :: num_iotasks
+
+namelist /io_nml/ num_iotasks, lredirect_stdout, log_filename, &
+ luse_pointer_files, pointer_filename, luse_nf_64bit_offset
+
+!------------------------------------------------------------------
+! 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
+
+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
+
+!------------------------------------------------------------------
+! The ncommas initial temperature and salinity namelist
+!------------------------------------------------------------------
+
+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
+
+namelist /init_ts_nml/ init_ts_option, init_ts_suboption, &
+ init_ts_file, init_ts_file_fmt, &
+ init_ts_outfile, init_ts_outfile_fmt
+
+!------------------------------------------------------------------
+! The ncommas 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 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.
+!
+!------------------------------------------------------------------
+!
+! 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)
+!
+!------------------------------------------------------------------
+
+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
+
+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
+
+!======================================================================
+contains
+!======================================================================
+
+
+subroutine initialize_module
+!------------------------------------------------------------------
+integer :: iunit, io
+
+! 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.
+
+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')
+
+! 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.
+
+! STOMP if ( allow_leapyear ) then
+ call set_calendar_type('gregorian')
+! STOMP else
+! STOMP call set_calendar_type('noleap')
+! STOMP endif
+
+! Read ncommas I/O information (for restart file ... grid dimensions)
+! Read ncommas initial information (for input/restart filename)
+
+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 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')
+
+! 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
+
+! 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
+
+! 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')
+
+! 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')
+
+! 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')
+
+module_initialized = .true.
+
+! Print module information to log file and stdout.
+call register_module(source, revision, revdate)
+
+end subroutine initialize_module
+
+
+
+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.
+
+integer, intent(out) :: Nx ! Number of Longitudes
+integer, intent(out) :: Ny ! Number of Latitudes
+
+integer :: grid_id, dimid, nc_rc
+
+if ( .not. module_initialized ) call initialize_module
+
+! get the ball rolling ...
+
+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
+
+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
+
+call nc_check(nf90_inquire_dimension(grid_id, dimid, len=Ny), &
+ 'get_horiz_grid_dims','inquire_dimension i '//trim(ic_filename))
+
+! 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
+
+
+
+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.
+!
+! Then, the DART time manager is queried to return what it knows ...
+!
+character(len=*), INTENT(OUT) :: calstring
+
+if ( .not. module_initialized ) call initialize_module
+
+call get_calendar_string(calstring)
+
+end subroutine get_ncommas_calendar
+
+
+
+function set_model_time_step()
+!------------------------------------------------------------------
+! the initialize_module ensures that the ncommas namelists are read.
+! The restart times in the ncommas_in&restart_nml are used to define
+! appropriate assimilation timesteps.
+!
+type(time_type) :: set_model_time_step
+
+if ( .not. module_initialized ) call initialize_module
+
+! Check the 'restart_freq_opt' and 'restart_freq' to 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
+
+
+
+
+subroutine write_ncommas_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_ncommas_namelist', msgstring, source, revision, revdate)
+endif
+
+! call print_date( model_time,'write_ncommas_namelist:dart model date')
+! call print_date(adv_to_time,'write_ncommas_namelist:advance_to date')
+! call print_time( model_time,'write_ncommas_namelist:dart model time')
+! call print_time(adv_to_time,'write_ncommas_namelist:advance_to time')
+! call print_time( offset,'write_ncommas_namelist:a distance of')
+! write( *,'(''write_ncommas_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_ncommas_namelist', &
+ 'stop_option must be "nday"', source, revision, revdate)
+endif
+
+iunit = open_file('ncommas_in.DART',form='formatted',action='rewind')
+write(iunit, nml=time_manager_nml)
+write(iunit, '('' '')')
+close(iunit)
+
+end subroutine write_ncommas_namelist
+
+
+
+ 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)
+!
+! mimic ncommas 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 ncommas
+
+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
+
+!*** 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
+
+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
+
+ TLAT(1,:) = (TLAT(2,:) + TLAT(nx,:))/c2
+ TLON(1,:) = (TLON(2,:) + TLON(nx,:))/c2
+
+else
+ write(msgstring,'(''ncommas_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
+
+! 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
+
+if ( .not. module_initialized ) call initialize_module
+
+filename = trim(ic_filename)
+
+end subroutine get_ncommas_restart_filename
+
+
+end module dart_ncommas_mod
Property changes on: DART/trunk/models/NCOMMAS/dart_ncommas_mod.f90
___________________________________________________________________
Added: svn:mime-type
+ text/plain
Added: svn:keywords
+ Date Rev Author HeadURL Id
Added: svn:eol-style
+ native
Added: DART/trunk/models/NCOMMAS/dart_to_ncommas.f90
===================================================================
--- DART/trunk/models/NCOMMAS/dart_to_ncommas.f90 (rev 0)
+++ DART/trunk/models/NCOMMAS/dart_to_ncommas.f90 2010-07-30 20:30:24 UTC (rev 4446)
@@ -0,0 +1,136 @@
+! DART software - Copyright \xA9 2004 - 2010 UCAR. This open source software is
+! provided by UCAR, "as is", without charge, subject to all terms of use at
+! http://www.image.ucar.edu/DAReS/DART/DART_download
+
+program dart_to_ncommas
+
+! <next few lines under version control, do not edit>
+! $URL$
+! $Id$
+! $Revision$
+! $Date$
+
+!----------------------------------------------------------------------
+! purpose: interface between DART and the ncommas model
+!
+! method: Read DART state vector and overwrite values in a ncommas restart file.
+! If the DART state vector has an 'advance_to_time' present, a
+! file called ncommas_in.DART is created with a time_manager_nml namelist
+! appropriate to advance ncommas to the requested time.
+!
+! The dart_to_ncommas_nml namelist setting for advance_time_present
+! determines whether or not the input file has an 'advance_to_time'.
+! Typically, only temporary files like 'assim_model_state_ic' have
+! an 'advance_to_time'.
+!
+! author: Tim Hoar 25 Jun 09, revised 12 July 2010
+!----------------------------------------------------------------------
+
+use types_mod, only : r8
+use utilities_mod, only : initialize_utilities, timestamp, &
+ 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, print_time, print_date, operator(-)
+use model_mod, only : static_init_model, sv_to_restart_file, &
+ get_model_size, get_ncommas_restart_filename
+use dart_ncommas_mod, only : write_ncommas_namelist
+
+implicit none
+
+! version controlled file description for error handling, do not edit
+character(len=128), parameter :: &
+ source = "$URL$", &
+ revision = "$Revision$", &
+ revdate = "$Date$"
+
+!------------------------------------------------------------------
+! The namelist variables
+!------------------------------------------------------------------
+
+character (len = 128) :: dart_to_ncommas_input_file = 'dart.ic'
+logical :: advance_time_present = .false.
+
+namelist /dart_to_ncommas_nml/ dart_to_ncommas_input_file, &
+ advance_time_present
+
+!----------------------------------------------------------------------
+
+integer :: iunit, io, x_size
+type(time_type) :: model_time, adv_to_time
+real(r8), allocatable :: statevector(:)
+character (len = 128) :: ncommas_restart_filename = 'no_ncommas_restart_file'
+logical :: verbose = .FALSE.
+
+!----------------------------------------------------------------------
+
+call initialize_utilities(progname='dart_to_ncommas', output_flag=verbose)
+
+!----------------------------------------------------------------------
+! Call model_mod:static_init_model() which reads the ncommas namelists
+! to set grid sizes, etc.
+!----------------------------------------------------------------------
+
+call static_init_model()
+
+x_size = get_model_size()
+allocate(statevector(x_size))
+
+! Read the namelist to get the input filename.
+
+call find_namelist_in_file("input.nml", "dart_to_ncommas_nml", iunit)
+read(iunit, nml = dart_to_ncommas_nml, iostat = io)
+call check_namelist_read(iunit, io, "dart_to_ncommas_nml")
+
+call get_ncommas_restart_filename( ncommas_restart_filename )
+
+write(*,*)
+write(*,'(''dart_to_ncommas:converting DART file '',A, &
+ &'' to ncommas restart file '',A)') &
+ trim(dart_to_ncommas_input_file), trim(ncommas_restart_filename)
+
+!----------------------------------------------------------------------
+! Reads the valid time, the state, and the target time.
+!----------------------------------------------------------------------
+
+iunit = open_restart_read(dart_to_ncommas_input_file)
+
+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 ncommas state vector
+! Convey the amount of time to integrate the model ...
+! time_manager_nml: stop_option, stop_count increments
+!----------------------------------------------------------------------
+
+call sv_to_restart_file(statevector, ncommas_restart_filename, model_time)
+
+if ( advance_time_present ) then
+ call write_ncommas_namelist(model_time, adv_to_time)
+endif
+
+!----------------------------------------------------------------------
+! Log what we think we're doing, and exit.
+!----------------------------------------------------------------------
+
+call print_date( model_time,'dart_to_ncommas:ncommas model date')
+call print_time( model_time,'dart_to_ncommas:DART model time')
+call print_date( model_time,'dart_to_ncommas:ncommas model date',logfileunit)
+call print_time( model_time,'dart_to_ncommas:DART model time',logfileunit)
+
+if ( advance_time_present ) then
+call print_time(adv_to_time,'dart_to_ncommas:advance_to time')
+call print_date(adv_to_time,'dart_to_ncommas:advance_to date')
+call print_time(adv_to_time,'dart_to_ncommas:advance_to time',logfileunit)
+call print_date(adv_to_time,'dart_to_ncommas:advance_to date',logfileunit)
+endif
+
+! When called with 'end', timestamp will call finalize_utilities()
+call timestamp(string1=source, pos='end')
+
+end program dart_to_ncommas
Property changes on: DART/trunk/models/NCOMMAS/dart_to_ncommas.f90
___________________________________________________________________
Added: svn:mime-type
+ text/plain
Added: svn:keywords
+ Date Rev Author HeadURL Id
Added: svn:eol-style
+ native
Added: DART/trunk/models/NCOMMAS/dart_to_ncommas.nml
===================================================================
--- DART/trunk/models/NCOMMAS/dart_to_ncommas.nml (rev 0)
+++ DART/trunk/models/NCOMMAS/dart_to_ncommas.nml 2010-07-30 20:30:24 UTC (rev 4446)
@@ -0,0 +1,4 @@
+&dart_to_ncommas_nml
+ dart_to_ncommas_input_file = 'dart.ic',
+ advance_time_present = .false. /
+
Property changes on: DART/trunk/models/NCOMMAS/dart_to_ncommas.nml
___________________________________________________________________
Added: svn:mime-type
+ text/plain
Added: svn:eol-style
+ native
Added: DART/trunk/models/NCOMMAS/model_mod.f90
===================================================================
--- DART/trunk/models/NCOMMAS/model_mod.f90 (rev 0)
+++ DART/trunk/models/NCOMMAS/model_mod.f90 2010-07-30 20:30:24 UTC (rev 4446)
@@ -0,0 +1,3378 @@
+! DART software - Copyright \xA9 2004 - 2010 UCAR. This open source software is
+! provided by UCAR, "as is", without charge, subject to all terms of use at
+! http://www.image.ucar.edu/DAReS/DART/DART_download
+
+module model_mod
+
+! <next few lines under version control, do not edit>
+! $URL$
+! $Id$
+! $Revision$
+! $Date$
+
+! This is the interface between the ncommas ocean model and DART.
+
+! Modules that are absolutely required for use are listed
+use types_mod, only : r4, r8, SECPERDAY, MISSING_R8, rad2deg, PI
+use time_manager_mod, only : time_type, set_time, set_date, get_date, get_time,&
+ print_time, print_date, set_calendar_type, &
+ operator(*), operator(+), operator(-), &
+ operator(>), operator(<), operator(/), &
+ operator(/=), operator(<=)
+use location_mod, only : location_type, get_dist, get_close_maxdist_init, &
+ get_close_obs_init, set_location, &
+ VERTISHEIGHT, get_location, vert_is_height, &
+ vert_is_level, vert_is_surface, &
+ loc_get_close_obs => get_close_obs, get_close_type
+use utilities_mod, only : register_module, error_handler, &
+ E_ERR, E_WARN, E_MSG, logfileunit, get_unit, &
+ nc_check, do_output, to_upper, &
+ find_namelist_in_file, check_namelist_read, &
+ open_file, file_exist, find_textfile_dims, &
+ file_to_text
+use obs_kind_mod, only : KIND_TEMPERATURE, KIND_SALINITY, KIND_DRY_LAND, &
+ KIND_U_CURRENT_COMPONENT, &
+ KIND_V_CURRENT_COMPONENT, KIND_SEA_SURFACE_HEIGHT, &
+ KIND_SEA_SURFACE_PRESSURE
+use mpi_utilities_mod, only: my_task_id
+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_topography, read_vert_grid, &
+ get_ncommas_restart_filename
+
+use typesizes
+use netcdf
+
+implicit none
+private
+
+! these routines must be public and you cannot change
+! the arguments - they will be called *from* the DART code.
+public :: get_model_size, &
+ adv_1step, &
+ get_state_meta_data, &
+ model_interpolate, &
+ get_model_time_step, &
+ static_init_model, &
+ end_model, &
+ init_time, &
+ init_conditions, &
+ nc_write_model_atts, &
+ nc_write_model_vars, &
+ pert_model_state, &
+ get_close_maxdist_init, &
+ get_close_obs_init, &
+ get_close_obs, &
+ ens_mean_for_model
+
+! 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
+
+! version controlled file description for error handling, do not edit
+character(len=128), parameter :: &
+ source = '$URL$', &
+ revision = '$Revision$', &
+ revdate = '$Date$'
+
+character(len=256) :: msgstring
+logical, save :: module_initialized = .false.
+
+! 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.
+integer :: debug = 0 ! turn up for more and more debug messages
+
+! 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, &
+ debug
+
+!------------------------------------------------------------------
+!
+! The DART state vector (control vector) will consist of: S, T, U, V, PSURF
+! (Salinity, Temperature, U velocity, V velocity, Sea Surface Height).
+! S, T are 3D arrays, located at cell centers. U,V are at grid cell corners.
+! PSURF is a 2D field (X,Y only). The Z direction is downward.
+!
+! FIXME: proposed change 1: we put SSH first, then T,U,V, then S, then
+! any optional tracers, since SSH is the only 2D
+! field; all tracers are 3D. this simplifies the
+! mapping to and from the vars to state vector.
+!
+! FIXME: proposed change 2: 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.
+!------------------------------------------------------------------
+
+integer, parameter :: n3dfields = 4
+integer, parameter :: n2dfields = 1
+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 ...
+character(len=128) :: progvarnames(nfields) = (/'SALT ','TEMP ','UVEL ','VVEL ','PSURF'/)
+
+integer, parameter :: S_index = 1
+integer, parameter :: T_index = 2
+integer, parameter :: U_index = 3
+integer, parameter :: V_index = 4
+integer, parameter :: PSURF_index = 5
+
+integer :: start_index(nfields)
+
+! 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
+
+! locations of cell centers (C) and edges (G) for each axis.
+real(r8), allocatable :: ZC(:), ZG(:)
+
+! 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(:,:)
+
+! integer, lowest valid cell number in the vertical
+integer, allocatable :: KMT(:, :), KMU(:, :)
+! real, depth of lowest valid cell (0 = land). use only if KMT/KMU not avail.
+real(r8), allocatable :: HT(:,:), HU(:,:)
+
+real(r8) :: endTime
+real(r8) :: ocean_dynamics_timestep = 900.0_r4
+integer :: timestepcount = 0
+type(time_type) :: model_time, model_timestep
+
+integer :: model_size ! the state vector length
+
+
+INTERFACE vector_to_prog_var
+ MODULE PROCEDURE vector_to_2d_prog_var
+ MODULE PROCEDURE vector_to_3d_prog_var
+END INTERFACE
+
+
+!------------------------------------------------
+
+! The regular grid used for dipole interpolation divides the sphere into
+! a set of regularly spaced lon-lat boxes. The number of boxes in
+! longitude and latitude are set by num_reg_x and num_reg_y. Making the
+! number of regular boxes smaller decreases the computation required for
+! doing each interpolation but increases the static storage requirements
+! and the initialization computation (which seems to be pretty small).
+integer, parameter :: num_reg_x = 90, num_reg_y = 90
+
+! The max_reg_list_num controls the size of temporary storage used for
+! initializing the regular grid. Four arrays
+! of size num_reg_x*num_reg_y*max_reg_list_num are needed. The initialization
+! fails and returns an error if max_reg_list_num is too small. A value of
+! 30 is sufficient for the x3 ncommas grid with 180 regular lon and lat boxes
+! and a value of 80 is sufficient for for the x1 grid.
+integer, parameter :: max_reg_list_num = 80
+
+! The dipole interpolation keeps a list of how many and which dipole quads
+! overlap each regular lon-lat box. The number for the u and t grids are
+! stored in u_dipole_num and t_dipole_num. The allocatable arrays
+! u_dipole_lon(lat)_list and t_dipole_lon(lat)_list list the longitude
+! and latitude indices for the overlapping dipole quads. The entry in
+! u_dipole_start and t_dipole_start for a given regular lon-lat box indicates
+! where the list of dipole quads begins in the u_dipole_lon(lat)_list and
+! t_dipole_lon(lat)_list arrays.
+
+integer :: u_dipole_start(num_reg_x, num_reg_y)
+integer :: u_dipole_num (num_reg_x, num_reg_y) = 0
+integer :: t_dipole_start(num_reg_x, num_reg_y)
+integer :: t_dipole_num (num_reg_x, num_reg_y) = 0
+integer, allocatable :: u_dipole_lon_list(:), t_dipole_lon_list(:)
+integer, allocatable :: u_dipole_lat_list(:), t_dipole_lat_list(:)
+
@@ Diff output truncated at 40000 characters. @@
More information about the Dart-dev
mailing list