[Dart-dev] [5780] DART/branches/development/models/mpas_atm: calls model_interpolate for a regular lat/lon/vert grid,
nancy at ucar.edu
nancy at ucar.edu
Fri Jul 13 16:37:11 MDT 2012
Revision: 5780
Author: nancy
Date: 2012-07-13 16:37:10 -0600 (Fri, 13 Jul 2012)
Log Message:
-----------
calls model_interpolate for a regular lat/lon/vert grid,
so the output is something we can look at with ncview.
input is a dart ics or restart state vector file. output
is both a .nc file and a matlab script containing the
same data.
Added Paths:
-----------
DART/branches/development/models/mpas_atm/resample.f90
DART/branches/development/models/mpas_atm/work/mkmf_resample
DART/branches/development/models/mpas_atm/work/path_names_resample
-------------- next part --------------
Added: DART/branches/development/models/mpas_atm/resample.f90
===================================================================
--- DART/branches/development/models/mpas_atm/resample.f90 (rev 0)
+++ DART/branches/development/models/mpas_atm/resample.f90 2012-07-13 22:37:10 UTC (rev 5780)
@@ -0,0 +1,362 @@
+! DART software - Copyright 2004 - 2011 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 resample
+
+! <next few lines under version control, do not edit>
+! $URL$
+! $Id$
+! $Revision$
+! $Date$
+
+!----------------------------------------------------------------------
+! purpose: resample a hexagonal mpas grid onto a regular lat/lon grid.
+! easier to plot with ncview for quick-n-dirty use.
+!----------------------------------------------------------------------
+
+use types_mod, only : r8, digits12, metadatalength, MISSING_R8, rad2deg
+use utilities_mod, only : initialize_utilities, finalize_utilities, nc_check, &
+ open_file, close_file, find_namelist_in_file, &
+ check_namelist_read, nmlfileunit, do_nml_file, do_nml_term, &
+ E_MSG, E_ERR, error_handler, get_unit
+use location_mod, only : location_type, set_location, write_location, get_dist, &
+ query_location, LocationDims, get_location, &
+ VERTISUNDEF, VERTISSURFACE, VERTISLEVEL, VERTISPRESSURE, &
+ VERTISHEIGHT, VERTISSCALEHEIGHT
+use obs_kind_mod, only : get_raw_obs_kind_name, get_raw_obs_kind_index, &
+ KIND_U_WIND_COMPONENT, KIND_V_WIND_COMPONENT
+use assim_model_mod, only : open_restart_read, open_restart_write, close_restart, &
+ aread_state_restart, awrite_state_restart, &
+ netcdf_file_type, aoutput_diagnostics, &
+ init_diag_output, finalize_diag_output, &
+ static_init_assim_model
+use time_manager_mod, only : time_type, set_calendar_type, GREGORIAN, &
+ read_time, get_time, set_time, &
+ print_date, get_date, &
+ print_time, write_time, &
+ operator(-)
+use model_mod, only : static_init_model, get_model_size, get_state_meta_data, &
+ model_interpolate, get_analysis_time, &
+ get_model_analysis_filename, analysis_file_to_statevector, &
+ statevector_to_analysis_file, get_analysis_time, &
+ write_model_time, get_grid_dims
+
+use netcdf
+use typesizes
+
+implicit none
+
+! version controlled file description for error handling, do not edit
+character(len=128), parameter :: &
+ source = "$URL$", &
+ revision = "$Revision$", &
+ revdate = "$Date$"
+
+character(len=256) :: string1, string2
+
+!------------------------------------------------------------------
+! The namelist variables
+!------------------------------------------------------------------
+
+character (len = 129) :: dart_input_file = 'filter_ics'
+character (len = 129) :: output_file = 'resampled'
+logical :: advance_time_present = .FALSE.
+logical :: verbose = .FALSE.
+real(r8) :: interp_test_dlon = 1.0
+real(r8) :: interp_test_dlat = 1.0
+real(r8) :: interp_test_dvert = 1.0
+real(r8), dimension(2) :: interp_test_latrange = (/ -90.0, 90.0 /)
+real(r8), dimension(2) :: interp_test_lonrange = (/ 0.0, 360.0 /)
+real(r8), dimension(2) :: interp_test_vertrange = (/ 1000.0, 30000.0 /)
+character(len=metadatalength) :: kind_of_interest = 'ANY'
+character(len=metadatalength) :: interp_test_vertcoord = 'VERTISHEIGHT'
+
+namelist /resample_nml/ dart_input_file, output_file, &
+ advance_time_present, &
+ kind_of_interest, verbose, &
+ interp_test_dlon, interp_test_lonrange, &
+ interp_test_dlat, interp_test_latrange, &
+ interp_test_dvert, interp_test_vertrange, &
+ interp_test_vertcoord
+
+!----------------------------------------------------------------------
+! integer :: numlons, numlats, numlevs
+
+integer :: ios_out, iunit, io, i
+integer :: x_size, skip
+integer :: mykindindex, vertcoord
+
+type(time_type) :: model_time, adv_to_time
+real(r8), allocatable :: statevector(:)
+
+character(len=metadatalength) :: state_meta(1)
+character(len=129) :: mpas_input_file ! set with get_model_analysis_filename() if needed
+type(netcdf_file_type) :: ncFileID
+type(location_type) :: loc
+
+real(r8) :: interp_val
+
+!----------------------------------------------------------------------
+! This portion checks the geometry information.
+!----------------------------------------------------------------------
+
+call initialize_utilities(progname='resample')
+call set_calendar_type(GREGORIAN)
+
+write(*,*)
+write(*,*)'Reading the namelist to get the input filename.'
+
+call find_namelist_in_file("input.nml", "resample_nml", iunit)
+read(iunit, nml = resample_nml, iostat = io)
+call check_namelist_read(iunit, io, "resample_nml")
+
+! Record the namelist values used for the run
+if (do_nml_file()) write(nmlfileunit, nml=resample_nml)
+if (do_nml_term()) write( * , nml=resample_nml)
+
+mykindindex = get_raw_obs_kind_index(kind_of_interest)
+
+! This harvests all kinds of initialization information
+
+call static_init_assim_model()
+
+x_size = get_model_size()
+allocate(statevector(x_size))
+
+write(*,*)
+write(*,*)'Reading '//trim(dart_input_file)
+
+iunit = open_restart_read(dart_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)
+
+call print_date( model_time,'resample:model date')
+call print_time( model_time,'resample:model time')
+
+!----------------------------------------------------------------------
+! Check the interpolation - print initially to STDOUT
+!----------------------------------------------------------------------
+
+write(*,*)
+write(*,*)'Testing single model_interpolate with ',trim(kind_of_interest),' ...'
+
+select case(trim(interp_test_vertcoord))
+ case ('VERTISUNDEF')
+ vertcoord = VERTISUNDEF
+ case ('VERTISSURFACE')
+ vertcoord = VERTISSURFACE
+ case ('VERTISLEVEL')
+ vertcoord = VERTISLEVEL
+ case ('VERTISPRESSURE')
+ vertcoord = VERTISPRESSURE
+ case ('VERTISHEIGHT')
+ vertcoord = VERTISHEIGHT
+ case ('VERTISSCALEHEIGHT')
+ vertcoord = VERTISSCALEHEIGHT
+ case default
+ write(string1,*) 'unknown vertcoord ', trim(interp_test_vertcoord)
+ call error_handler(E_ERR,'test_interpolate',string1,source,revision,revdate)
+end select
+
+ios_out = test_interpolate()
+
+if ( ios_out == 0 ) then
+ write(*,*)'Rigorous model_interpolate SUCCESS.'
+else
+ write(*,*)'Rigorous model_interpolate WARNING: model_interpolate had ', ios_out, ' failures.'
+endif
+
+!----------------------------------------------------------------------
+! This must be the last few lines of the main program.
+!----------------------------------------------------------------------
+
+ 999 continue
+
+call finalize_utilities()
+
+!----------------------------------------------------------------------
+!----------------------------------------------------------------------
+
+contains
+
+
+function test_interpolate()
+! function to exercise the model_mod:model_interpolate() function
+! This will result in a netCDF file with all salient metadata
+integer :: test_interpolate
+
+! Local variables
+
+real(r8), allocatable :: lon(:), lat(:), vert(:)
+real(r8), allocatable :: field(:,:,:)
+integer :: nlon, nlat, nvert
+integer :: ilon, jlat, kvert, nfailed
+character(len=128) :: ncfilename,txtfilename
+
+character(len=8) :: crdate ! needed by F90 DATE_AND_TIME intrinsic
+character(len=10) :: crtime ! needed by F90 DATE_AND_TIME intrinsic
+character(len=5) :: crzone ! needed by F90 DATE_AND_TIME intrinsic
+integer, dimension(8) :: values ! needed by F90 DATE_AND_TIME intrinsic
+
+integer :: ncid, nlonDimID, nlatDimID, nvertDimID
+integer :: VarID, lonVarID, latVarID, vertVarID
+
+test_interpolate = 0 ! normal termination
+
+if ((interp_test_dlon < 0.0_r8) .or. (interp_test_dlat < 0.0_r8)) then
+ write(*,*)'Skipping the rigorous interpolation test because one of'
+ write(*,*)'interp_test_dlon,interp_test_dlat are < 0.0'
+ write(*,*)'interp_test_dlon = ',interp_test_dlon
+ write(*,*)'interp_test_dlat = ',interp_test_dlat
+ write(*,*)'interp_test_dvert = ',interp_test_dvert
+endif
+
+write( ncfilename,'(a,a)')trim(output_file),'.nc'
+write(txtfilename,'(a,a)')trim(output_file),'.m'
+
+! round down to avoid exceeding the specified range
+nlat = aint(( interp_test_latrange(2) - interp_test_latrange(1))/interp_test_dlat) + 1
+nlon = aint(( interp_test_lonrange(2) - interp_test_lonrange(1))/interp_test_dlon) + 1
+nvert = aint((interp_test_vertrange(2) - interp_test_vertrange(1))/interp_test_dvert) + 1
+
+iunit = open_file(trim(txtfilename), action='write')
+write(iunit,'(''missingvals = '',f12.4,'';'')')MISSING_R8
+write(iunit,'(''nlon = '',i8,'';'')')nlon
+write(iunit,'(''nlat = '',i8,'';'')')nlat
+write(iunit,'(''nvert = '',i8,'';'')')nvert
+write(iunit,'(''interptest = [ ... '')')
+
+allocate(lon(nlon), lat(nlat), vert(nvert), field(nlon,nlat,nvert))
+nfailed = 0
+
+do ilon = 1, nlon
+ lon(ilon) = interp_test_lonrange(1) + real(ilon-1,r8) * interp_test_dlon
+ do jlat = 1, nlat
+ lat(jlat) = interp_test_latrange(1) + real(jlat-1,r8) * interp_test_dlat
+ do kvert = 1, nvert
+ vert(kvert) = interp_test_vertrange(1) + real(kvert-1,r8) * interp_test_dvert
+
+ loc = set_location(lon(ilon), lat(jlat), vert(kvert), vertcoord)
+
+ call model_interpolate(statevector, loc, mykindindex, field(ilon,jlat,kvert), ios_out)
+ write(iunit,*) field(ilon,jlat,kvert)
+
+ if (ios_out /= 0) then
+ if (verbose) then
+ write(string2,'(''ilon,jlat,kvert,lon,lat,vert'',3(1x,i6),3(1x,f14.6))') &
+ ilon,jlat,kvert,lon(ilon),lat(jlat),vert(kvert)
+ write(string1,*) 'interpolation return code was', ios_out
+ call error_handler(E_MSG,'test_interpolate',string1,source,revision,revdate,text2=string2)
+ endif
+ nfailed = nfailed + 1
+ endif
+
+ enddo
+ end do
+end do
+write(iunit,'(''];'')')
+write(iunit,'(''datmat = reshape(interptest,nvert,nlat,nlon);'')')
+write(iunit,'(''datmat(datmat == missingvals) = NaN;'')')
+call close_file(iunit)
+
+write(*,*) 'total interpolations, num failed: ', nlon*nlat*nvert, nfailed
+
+! Write out the netCDF file for easy exploration.
+
+call DATE_AND_TIME(crdate,crtime,crzone,values)
+write(string1,'(''YYYY MM DD HH MM SS = '',i4,5(1x,i2.2))') &
+ values(1), values(2), values(3), values(5), values(6), values(7)
+
+call nc_check( nf90_create(path=trim(ncfilename), cmode=NF90_clobber, ncid=ncid), &
+ 'test_interpolate', 'open '//trim(ncfilename))
+call nc_check( nf90_put_att(ncid, NF90_GLOBAL, 'creation_date' ,trim(string1) ), &
+ 'test_interpolate', 'creation put '//trim(ncfilename))
+call nc_check( nf90_put_att(ncid, NF90_GLOBAL, 'parent_file', dart_input_file ), &
+ 'test_interpolate', 'put_att filename '//trim(ncfilename))
+
+! Define dimensions
+
+call nc_check(nf90_def_dim(ncid=ncid, name='lon', len=nlon, &
+ dimid = nlonDimID),'test_interpolate', 'nlon def_dim '//trim(ncfilename))
+
+call nc_check(nf90_def_dim(ncid=ncid, name='lat', len=nlat, &
+ dimid = nlatDimID),'test_interpolate', 'nlat def_dim '//trim(ncfilename))
+
+call nc_check(nf90_def_dim(ncid=ncid, name='vert', len=nvert, &
+ dimid = nvertDimID),'test_interpolate', 'nvert def_dim '//trim(ncfilename))
+
+! Define variables
+
+call nc_check(nf90_def_var(ncid=ncid, name='lon', xtype=nf90_double, &
+ dimids=nlonDimID, varid=lonVarID), 'test_interpolate', &
+ 'lon def_var '//trim(ncfilename))
+call nc_check(nf90_put_att(ncid, lonVarID, 'range', interp_test_lonrange), &
+ 'test_interpolate', 'put_att lonrange '//trim(ncfilename))
+call nc_check(nf90_put_att(ncid, lonVarID, 'cartesian_axis', 'X'), &
+ 'test_interpolate', 'lon cartesian_axis '//trim(ncfilename))
+
+
+call nc_check(nf90_def_var(ncid=ncid, name='lat', xtype=nf90_double, &
+ dimids=nlatDimID, varid=latVarID), 'test_interpolate', &
+ 'lat def_var '//trim(ncfilename))
+call nc_check(nf90_put_att(ncid, latVarID, 'range', interp_test_latrange), &
+ 'test_interpolate', 'put_att latrange '//trim(ncfilename))
+call nc_check(nf90_put_att(ncid, latVarID, 'cartesian_axis', 'Y'), &
+ 'test_interpolate', 'lat cartesian_axis '//trim(ncfilename))
+
+call nc_check(nf90_def_var(ncid=ncid, name='vert', xtype=nf90_double, &
+ dimids=nvertDimID, varid=vertVarID), 'test_interpolate', &
+ 'vert def_var '//trim(ncfilename))
+call nc_check(nf90_put_att(ncid, vertVarID, 'range', interp_test_vertrange), &
+ 'test_interpolate', 'put_att vertrange '//trim(ncfilename))
+call nc_check(nf90_put_att(ncid, vertVarID, 'cartesian_axis', 'Z'), &
+ 'test_interpolate', 'vert cartesian_axis '//trim(ncfilename))
+
+call nc_check(nf90_def_var(ncid=ncid, name='field', xtype=nf90_double, &
+ dimids=(/ nlonDimID, nlatDimID, nvertDimID /), varid=VarID), 'test_interpolate', &
+ 'field def_var '//trim(ncfilename))
+call nc_check(nf90_put_att(ncid, VarID, 'long_name', kind_of_interest), &
+ 'test_interpolate', 'put_att field long_name '//trim(ncfilename))
+call nc_check(nf90_put_att(ncid, VarID, '_FillValue', MISSING_R8), &
+ 'test_interpolate', 'put_att field FillValue '//trim(ncfilename))
+call nc_check(nf90_put_att(ncid, VarID, 'missing_value', MISSING_R8), &
+ 'test_interpolate', 'put_att field missing_value '//trim(ncfilename))
+call nc_check(nf90_put_att(ncid, VarID, 'vertcoord_string', interp_test_vertcoord ), &
+ 'test_interpolate', 'put_att field vertcoord_string '//trim(ncfilename))
+call nc_check(nf90_put_att(ncid, VarID, 'vertcoord', vertcoord ), &
+ 'test_interpolate', 'put_att field vertcoord '//trim(ncfilename))
+
+! Leave define mode so we can fill the variables.
+call nc_check(nf90_enddef(ncid), &
+ 'test_interpolate','field enddef '//trim(ncfilename))
+
+! Fill the variables
+call nc_check(nf90_put_var(ncid, lonVarID, lon), &
+ 'test_interpolate','lon put_var '//trim(ncfilename))
+call nc_check(nf90_put_var(ncid, latVarID, lat), &
+ 'test_interpolate','lat put_var '//trim(ncfilename))
+call nc_check(nf90_put_var(ncid, vertVarID, vert), &
+ 'test_interpolate','vert put_var '//trim(ncfilename))
+call nc_check(nf90_put_var(ncid, VarID, field), &
+ 'test_interpolate','field put_var '//trim(ncfilename))
+
+! tidy up
+call nc_check(nf90_close(ncid), &
+ 'test_interpolate','close '//trim(ncfilename))
+
+deallocate(lon, lat, vert, field)
+
+test_interpolate = nfailed
+
+end function test_interpolate
+
+
+
+
+
+end program resample
Property changes on: DART/branches/development/models/mpas_atm/resample.f90
___________________________________________________________________
Added: svn:mime-type
+ text/plain
Added: svn:keywords
+ Date Rev Author HeadURL Id
Added: svn:eol-style
+ native
Added: DART/branches/development/models/mpas_atm/work/mkmf_resample
===================================================================
--- DART/branches/development/models/mpas_atm/work/mkmf_resample (rev 0)
+++ DART/branches/development/models/mpas_atm/work/mkmf_resample 2012-07-13 22:37:10 UTC (rev 5780)
@@ -0,0 +1,18 @@
+#!/bin/csh
+#
+# 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
+#
+# $Id$
+
+../../../mkmf/mkmf -p resample -t ../../../mkmf/mkmf.template \
+ -a "../../.." path_names_resample
+
+exit $status
+
+# <next few lines under version control, do not edit>
+# $URL$
+# $Revision$
+# $Date$
+
Property changes on: DART/branches/development/models/mpas_atm/work/mkmf_resample
___________________________________________________________________
Added: svn:executable
+ *
Added: svn:mime-type
+ text/plain
Added: svn:keywords
+ Date Rev Author HeadURL Id
Added: svn:eol-style
+ native
Added: DART/branches/development/models/mpas_atm/work/path_names_resample
===================================================================
--- DART/branches/development/models/mpas_atm/work/path_names_resample (rev 0)
+++ DART/branches/development/models/mpas_atm/work/path_names_resample 2012-07-13 22:37:10 UTC (rev 5780)
@@ -0,0 +1,15 @@
+assim_model/assim_model_mod.f90
+common/types_mod.f90
+location/threed_sphere/location_mod.f90
+models/mpas_atm/resample.f90
+models/mpas_atm/model_mod.f90
+models/mpas_atm/get_geometry_mod.f90
+models/mpas_atm/get_coeff_mod.f90
+models/mpas_atm/get_reconstruct_mod.f90
+mpi_utilities/null_mpi_utilities_mod.f90
+obs_def/obs_def_mod.f90
+obs_kind/obs_kind_mod.f90
+random_nr/random_nr_mod.f90
+random_seq/random_seq_mod.f90
+time_manager/time_manager_mod.f90
+utilities/utilities_mod.f90
More information about the Dart-dev
mailing list