[Dart-dev] [3280] DART/trunk/models/MITgcm_ocean/model_mod.f90:
vector_to_prog_var rewritten to accomodate individual variable extraction.
thoar at subversion.ucar.edu
thoar at subversion.ucar.edu
Thu Mar 27 14:30:22 MDT 2008
An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/dart-dev/attachments/20080327/8e55faa9/attachment-0001.html
-------------- next part --------------
Modified: DART/trunk/models/MITgcm_ocean/model_mod.f90
===================================================================
--- DART/trunk/models/MITgcm_ocean/model_mod.f90 2008-03-19 21:50:53 UTC (rev 3279)
+++ DART/trunk/models/MITgcm_ocean/model_mod.f90 2008-03-27 20:30:22 UTC (rev 3280)
@@ -51,18 +51,17 @@
! the interfaces here can be changed as appropriate.
public :: prog_var_to_vector, vector_to_prog_var, &
MIT_meta_type, read_meta, write_meta, &
- read_snapshot, drop_snapshot, write_snapshot, &
+ read_snapshot, write_snapshot, &
snapshot_files_to_sv, sv_to_snapshot_files
-
! version controlled file description for error handling, do not edit
character(len=128), parameter :: &
source = "$URL$", &
revision = "$Revision$", &
revdate = "$Date: 2007-04-03 16:44:36 -0600 (Tue, 03 Apr 2007) $"
-character(len=129) :: msgstring, errstring
+character(len=129) :: msgstring
!! FIXME: on ifort, this must be 1
@@ -93,54 +92,13 @@
integer, parameter :: max_nz = 512
integer, parameter :: max_nr = 512
-! must match lists declared in ini_parms.f
-
-!-- Time stepping parameters namelist
- NAMELIST /PARM03/ &
- nIter0, nTimeSteps, nEndIter, pickupSuff, &
- deltaT, deltaTClock, deltaTmom, &
- deltaTtracer, dTtracerLev, deltaTfreesurf, &
- forcing_In_AB, momForcingOutAB, tracForcingOutAB, &
- momDissip_In_AB, doAB_onGtGs, &
- abEps, alph_AB, beta_AB, startFromPickupAB2, &
- tauCD, rCD, &
- baseTime, startTime, endTime, chkPtFreq, &
- dumpFreq, dumpInitAndLast, adjDumpFreq, taveFreq, tave_lastIter, &
- diagFreq, monitorFreq, adjMonitorFreq, pChkPtFreq, cAdjFreq, &
- outputTypesInclusive, &
- tauThetaClimRelax, tauSaltClimRelax, latBandClimRelax, &
- tauThetaClimRelax3Dim, tauSaltClimRelax3Dim, tauTr1ClimRelax, &
- periodicExternalForcing, externForcingPeriod, externForcingCycle, &
- calendarDumps
-
-!-- Gridding parameters namelist
- NAMELIST /PARM04/ &
- usingCartesianGrid, usingCylindricalGrid, &
- dxSpacing, dySpacing, delX, delY, delXFile, delYFile, &
- usingSphericalPolarGrid, phiMin, thetaMin, rSphere, &
- usingCurvilinearGrid, horizGridFile, deepAtmosphere, &
- Ro_SeaLevel, delZ, delP, delR, delRc, delRFile, delRcFile, &
- rkFac, groundAtK1
-
-!-- Input files namelist
- NAMELIST /PARM05/ &
- bathyFile, topoFile, shelfIceFile, &
- hydrogThetaFile, hydrogSaltFile, diffKrFile, &
- zonalWindFile, meridWindFile, &
- thetaClimFile, saltClimFile, &
- surfQfile, surfQnetFile, surfQswFile, EmPmRfile, saltFluxFile, &
- lambdaThetaFile, lambdaSaltFile, &
- uVelInitFile, vVelInitFile, pSurfInitFile, &
- dQdTFile, ploadFile,tCylIn,tCylOut, &
- eddyTauxFile, eddyTauyFile, &
- mdsioLocalDir, &
- the_run_name
-
!-- FIXME: i have not been able to find all of these in the
!-- original source code. the ones we're using have been
!-- checked that they are the right type. the others might
!-- cause problems since i don't have a namelist file i can
!-- examine that has the full set of values specified.
+!--
+!-- must match lists declared in ini_parms.f
!-- Time stepping parameters variable declarations
real(r8) :: nIter0, nTimeSteps, nEndIter, pickupSuff, &
@@ -187,7 +145,47 @@
mdsioLocalDir, &
the_run_name
+!-- Time stepping parameters namelist
+NAMELIST /PARM03/ &
+ nIter0, nTimeSteps, nEndIter, pickupSuff, &
+ deltaT, deltaTClock, deltaTmom, &
+ deltaTtracer, dTtracerLev, deltaTfreesurf, &
+ forcing_In_AB, momForcingOutAB, tracForcingOutAB, &
+ momDissip_In_AB, doAB_onGtGs, &
+ abEps, alph_AB, beta_AB, startFromPickupAB2, &
+ tauCD, rCD, &
+ baseTime, startTime, endTime, chkPtFreq, &
+ dumpFreq, dumpInitAndLast, adjDumpFreq, taveFreq, tave_lastIter, &
+ diagFreq, monitorFreq, adjMonitorFreq, pChkPtFreq, cAdjFreq, &
+ outputTypesInclusive, &
+ tauThetaClimRelax, tauSaltClimRelax, latBandClimRelax, &
+ tauThetaClimRelax3Dim, tauSaltClimRelax3Dim, tauTr1ClimRelax, &
+ periodicExternalForcing, externForcingPeriod, externForcingCycle, &
+ calendarDumps
+!-- Gridding parameters namelist
+NAMELIST /PARM04/ &
+ usingCartesianGrid, usingCylindricalGrid, &
+ dxSpacing, dySpacing, delX, delY, delXFile, delYFile, &
+ usingSphericalPolarGrid, phiMin, thetaMin, rSphere, &
+ usingCurvilinearGrid, horizGridFile, deepAtmosphere, &
+ Ro_SeaLevel, delZ, delP, delR, delRc, delRFile, delRcFile, &
+ rkFac, groundAtK1
+
+!-- Input files namelist
+NAMELIST /PARM05/ &
+ bathyFile, topoFile, shelfIceFile, &
+ hydrogThetaFile, hydrogSaltFile, diffKrFile, &
+ zonalWindFile, meridWindFile, &
+ thetaClimFile, saltClimFile, &
+ surfQfile, surfQnetFile, surfQswFile, EmPmRfile, saltFluxFile, &
+ lambdaThetaFile, lambdaSaltFile, &
+ uVelInitFile, vVelInitFile, pSurfInitFile, &
+ dQdTFile, ploadFile,tCylIn,tCylOut, &
+ eddyTauxFile, eddyTauyFile, &
+ mdsioLocalDir, &
+ the_run_name
+
!------------------------------------------------------------------
!
! The DART state vector (control vector) will consist of: S, T, U, V, SSH
@@ -200,12 +198,9 @@
!
!------------------------------------------------------------------
-! Grid parameters that we define - the values will be read from a
-! standard MITgcm namelist and filled in here.
-
-integer, parameter :: nfields = 5
integer, parameter :: n3dfields = 4
integer, parameter :: n2dfields = 1
+integer, parameter :: nfields = n3dfields + n2dfields
integer, parameter :: S_index = 1
integer, parameter :: T_index = 2
@@ -213,14 +208,17 @@
integer, parameter :: V_index = 4
integer, parameter :: SSH_index = 5
-! FIXME: put the input/restart data filenames into an array
-! so we can loop over all 3d arrays and then 2d arrays as we
-! fill the state vector.
+! (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) = (/'S ','T ','U ','V ','SSH'/)
+integer :: start_index(nfields)
+! Grid parameters - the values will be read from a
+! standard MITgcm namelist and filled in here.
+
! grid counts for each field
integer :: Nx, Ny, Nz
-integer :: start_index(nfields)
! locations of the cell centers (C) and cell faces (grid?) (G)
! for each axis.
@@ -234,7 +232,6 @@
!real(r8) :: delta_lat, delta_lon, delta_depth
!real(r8), allocatable :: lat_grid(:), lon_grid(:), depth_grid(:)
-
! What is the natural model timestep?
real(r8) :: model_timestep
type(time_type) :: time_step
@@ -244,7 +241,6 @@
integer :: time_step_days = 0
integer :: time_step_seconds = 900
-
! the state vector length
integer :: model_size
@@ -272,16 +268,15 @@
MODULE PROCEDURE read_3d_snapshot
END INTERFACE
-INTERFACE drop_snapshot
- MODULE PROCEDURE drop_3d_snapshot
- MODULE PROCEDURE drop_2d_snapshot
-END INTERFACE
-
INTERFACE write_snapshot
MODULE PROCEDURE write_2d_snapshot
MODULE PROCEDURE write_3d_snapshot
END INTERFACE
+INTERFACE vector_to_prog_var
+ MODULE PROCEDURE vector_to_2d_prog_var
+ MODULE PROCEDURE vector_to_3d_prog_var
+END INTERFACE
contains
@@ -297,7 +292,7 @@
! it reads in the grid information and then the model data.
integer :: iunit, io
-integer :: i, j, k
+integer :: i
real(r4), allocatable :: data_2d_array(:,:)
! The Plan:
@@ -371,8 +366,8 @@
endif
enddo
if (Nx == -1) then
- write(errstring,*)'could not figure out number of longitudes from delX in namelist'
- call error_handler(E_ERR,"static_init_model", errstring, source, revision, revdate)
+ write(msgstring,*)'could not figure out number of longitudes from delX in namelist'
+ call error_handler(E_ERR,"static_init_model", msgstring, source, revision, revdate)
endif
Ny = -1
@@ -383,8 +378,8 @@
endif
enddo
if (Ny == -1) then
- write(errstring,*)'could not figure out number of latitudes from delY in namelist'
- call error_handler(E_ERR,"static_init_model", errstring, source, revision, revdate)
+ write(msgstring,*)'could not figure out number of latitudes from delY in namelist'
+ call error_handler(E_ERR,"static_init_model", msgstring, source, revision, revdate)
endif
Nz = -1
@@ -395,105 +390,53 @@
endif
enddo
if (Nz == -1) then
- write(errstring,*)'could not figure out number of depth levels from delZ in namelist'
- call error_handler(E_ERR,"static_init_model", errstring, source, revision, revdate)
+ write(msgstring,*)'could not figure out number of depth levels from delZ in namelist'
+ call error_handler(E_ERR,"static_init_model", msgstring, source, revision, revdate)
endif
+! We know enough to allocate grid variables.
+! 'read_snapshot' checks the ".meta" files for agreeable dimensions.
+
+allocate(XC(Nx), YC(Ny), ZC(Nz))
+allocate(XG(Nx), YG(Ny), ZG(Nz))
+allocate(data_2d_array(Nx,Ny))
+
! determine longitudes
call read_snapshot('XC', data_2d_array, timestepcount)
-if (size(data_2d_array, 1) /= Nx) then
- write(errstring,*)'XC.meta size for dim 1 does not match delX grid size in namelist'
- write(errstring,*)'delX grid size in namelist does not match XC.meta size, dim 1'
- call error_handler(E_ERR,"static_init_model", errstring, source, revision, revdate)
-endif
-if (size(data_2d_array, 2) /= Ny) then
- write(errstring,*)'XC.meta size for dim 2 does not match delY grid size in namelist'
- call error_handler(E_ERR,"static_init_model", errstring, source, revision, revdate)
-endif
-allocate(XC(Nx))
do i=1, Nx
XC(i) = data_2d_array(i, 1)
enddo
-call drop_snapshot(data_2d_array)
call read_snapshot('XG', data_2d_array, timestepcount)
-if (size(data_2d_array, 1) /= Nx) then
- write(errstring,*)'XG.meta size for dim 1 does not match delX grid size in namelist'
- call error_handler(E_ERR,"static_init_model", errstring, source, revision, revdate)
-endif
-if (size(data_2d_array, 2) /= Ny) then
- write(errstring,*)'XG.meta size for dim 2 does not match delY grid size in namelist'
- call error_handler(E_ERR,"static_init_model", errstring, source, revision, revdate)
-endif
-allocate(XG(Nx))
do i=1, Nx
XG(i) = data_2d_array(i, 1)
enddo
-call drop_snapshot(data_2d_array)
! determine latitudes
call read_snapshot('YC', data_2d_array, timestepcount)
-if (size(data_2d_array, 1) /= Nx) then
- write(errstring,*)'YC.meta size for dim 1 does not match delX grid size in namelist'
- call error_handler(E_ERR,"static_init_model", errstring, source, revision, revdate)
-endif
-if (size(data_2d_array, 2) /= Ny) then
- write(errstring,*)'YC.meta size for dim 2 does not match delY grid size in namelist'
- call error_handler(E_ERR,"static_init_model", errstring, source, revision, revdate)
-endif
-allocate(YC(Ny))
do i=1, Ny
- YC(i) = data_2d_array(i, 1)
+ YC(i) = data_2d_array(1, i)
enddo
-call drop_snapshot(data_2d_array)
call read_snapshot('YG', data_2d_array, timestepcount)
-if (size(data_2d_array, 1) /= Nx) then
- write(errstring,*)'YG.meta size for dim 1 does not match delX grid size in namelist'
- call error_handler(E_ERR,"static_init_model", errstring, source, revision, revdate)
-endif
-if (size(data_2d_array, 2) /= Ny) then
- write(errstring,*)'YG.meta size for dim 2 does not match delY grid size in namelist'
- call error_handler(E_ERR,"static_init_model", errstring, source, revision, revdate)
-endif
-allocate(YG(Ny))
do i=1, Ny
- YG(i) = data_2d_array(i, 1)
+ YG(i) = data_2d_array(1, i)
enddo
-call drop_snapshot(data_2d_array)
-! the namelist contains a list of thicknesses of each
-! depth level; the ZC array is a list of the depths of
-! each cell center, so they must be computed.
+deallocate(data_2d_array)
-allocate(ZC(Nz))
+! the namelist contains a list of thicknesses of each depth level (delZ)
+! ZG (the grid edges) and ZC (the grid centroids) must be computed.
+
+ZG(1) = 0.0_r8
ZC(1) = -0.5_r8 * delZ(1)
do i=2, Nz
- ZC(i) = ZC(i-1) - (0.5_r8 * delZ(i-1) + 0.5_r8 * delZ(i))
+ ZG(i) = ZG(i-1) - delZ(i-1)
+ ZC(i) = ZC(i-1) - 0.5_r8 * delZ(i-1) - 0.5_r8 * delZ(i)
enddo
-!! DEBUG: since we are computing these, check to be sure
-!! they look right.
-!do i=1, Nz
-! print *, 'i, delZ(i), ZC(i) = ', i, delZ(i), ZC(i)
-!enddo
-
-! for now, leave ZG undefined. it really does have one
-! more depth than we read in, so it isn't clear how to
-! handle this.
-allocate(ZG(Nz))
-ZG(:) = 0.0_r8
-! if we did need depth faces, here is how to compute them.
-! the delZ array are each cell width.
-!allocate(ZG(Nz+1))
-!ZG(1) = 0.0_r8
-!do i=2, Nz+1
-! ZG(i) = ZG(i-1) - delZ(i-1)
-!enddo
-
-
! record where in the state vector the data type changes
! from one type to another, by computing the starting
! index for each block of data.
@@ -1206,9 +1149,9 @@
"nc_write_model_atts", "time dimid "//trim(filename))
if ( TimeDimID /= unlimitedDimId ) then
- write(errstring,*)"Time Dimension ID ",TimeDimID, &
+ write(msgstring,*)"Time Dimension ID ",TimeDimID, &
" should equal Unlimited Dimension ID",unlimitedDimID
- call error_handler(E_ERR,"nc_write_model_atts", errstring, source, revision, revdate)
+ call error_handler(E_ERR,"nc_write_model_atts", msgstring, source, revision, revdate)
endif
!-------------------------------------------------------------------------------
@@ -1507,12 +1450,11 @@
integer :: ierr ! return value of function
integer :: nDimensions, nVariables, nAttributes, unlimitedDimID
-integer :: VarID, i, j, k
+integer :: VarID
-real(r4), allocatable, dimension(:,:,:) :: s,t,u,v
-real(r4), allocatable, dimension(:,:) :: ssh
+real(r4), dimension(Nx,Ny,Nz) :: data_3d
+real(r4), dimension(Nx,Ny) :: data_2d
character(len=128) :: filename
-!real(r4), dimension(4) :: mask
ierr = -1 ! assume things go poorly
@@ -1543,96 +1485,46 @@
!----------------------------------------------------------------------------
! We need to process the prognostic variables.
- !----------------------------------------------------------------------------
-
- call vector_to_prog_var(statevec,s,t,u,v,ssh) ! arrays allocated internally
-
! Replace missing values (0.0) with netcdf missing value.
! Staggered grid causes some logistical problems.
! Hopefully, the conversion between r8 and r4 still preserves 'hard' zeros.
+ !----------------------------------------------------------------------------
-! do j = 1,Ny ! latitudes
-! do i = 1,Nx ! longitudes
-
-! ! These three live on the same grid (the centers) ... if these are missing
-! ! at the surface, they are truly missing.
-
-! if (( s(i,j,1) == 0.0_r4) .and. &
-! ( t(i,j,1) == 0.0_r4) .and. &
-! (ssh(i,j ) == 0.0_r4)) then
-! s(i,j,:) = NF90_FILL_REAL
-! t(i,j,:) = NF90_FILL_REAL
-! ssh(i,j ) = NF90_FILL_REAL
-! endif
-! enddo
-! enddo
-
- ! Check for other ways to be missing (dry bottom)
-
-! do k = 2,Nz ! vertical
-! do j = 1,Ny ! latitudes
-! do i = 1,Nx ! longitudes
-! if ((s(i,j,k) == 0.0_r4) .and. (t(i,j,k) == 0.0_r4)) then
-! s(i,j,k) = NF90_FILL_REAL
-! t(i,j,k) = NF90_FILL_REAL
-! endif
-! enddo
-! enddo
-! enddo
-
- ! U,V live on staggered grids ... no way to hedge against errant zeros
-
- Ver: do k = 1,Nz ! vertical
- Lat: do j = 1,Ny ! latitudes
- Lon: do i = 1,Nx ! longitudes
-
- if (s(i,j,k) == 0.0_r4) s(i,j,k) = NF90_FILL_REAL
- if (t(i,j,k) == 0.0_r4) t(i,j,k) = NF90_FILL_REAL
- if (u(i,j,k) == 0.0_r4) u(i,j,k) = NF90_FILL_REAL
- if (v(i,j,k) == 0.0_r4) v(i,j,k) = NF90_FILL_REAL
- if (ssh(i,j) == 0.0_r4) ssh(i,j) = NF90_FILL_REAL
-
-! mask(1) = s(i,j,k)
-! mask(2) = t(i,j,k)
-! mask(3) = u(i,j,k)
-! mask(4) = v(i,j,k)
-
-! if (any(mask) .and. .not. all(mask)) then
-! not sure this works on a staggered grid ...
-! endif
-
-
- enddo Lon
- enddo Lat
- enddo Ver
-
+ call vector_to_prog_var(statevec,S_index,data_3d)
+ where (data_3d == 0.0_r4) data_3d = NF90_FILL_REAL
call nc_check(NF90_inq_varid(ncFileID, "S", VarID), &
"nc_write_model_vars", "S inq_varid "//trim(filename))
- call nc_check(nf90_put_var(ncFileID,VarID,s,start=(/1,1,1,copyindex,timeindex/)),&
+ call nc_check(nf90_put_var(ncFileID,VarID,data_3d,start=(/1,1,1,copyindex,timeindex/)),&
"nc_write_model_vars", "S put_var "//trim(filename))
+ call vector_to_prog_var(statevec,T_index,data_3d)
+ where (data_3d == 0.0_r4) data_3d = NF90_FILL_REAL
call nc_check(NF90_inq_varid(ncFileID, "T", VarID), &
"nc_write_model_vars", "T inq_varid "//trim(filename))
- call nc_check(nf90_put_var(ncFileID,VarID,t,start=(/1,1,1,copyindex,timeindex/)),&
+ call nc_check(nf90_put_var(ncFileID,VarID,data_3d,start=(/1,1,1,copyindex,timeindex/)),&
"nc_write_model_vars", "T put_var "//trim(filename))
+ call vector_to_prog_var(statevec,U_index,data_3d)
+ where (data_3d == 0.0_r4) data_3d = NF90_FILL_REAL
call nc_check(NF90_inq_varid(ncFileID, "U", VarID), &
"nc_write_model_vars", "U inq_varid "//trim(filename))
- call nc_check(nf90_put_var(ncFileID,VarID,u,start=(/1,1,1,copyindex,timeindex/)),&
+ call nc_check(nf90_put_var(ncFileID,VarID,data_3d,start=(/1,1,1,copyindex,timeindex/)),&
"nc_write_model_vars", "U put_var "//trim(filename))
+ call vector_to_prog_var(statevec,V_index,data_3d)
+ where (data_3d == 0.0_r4) data_3d = NF90_FILL_REAL
call nc_check(NF90_inq_varid(ncFileID, "V", VarID), &
"nc_write_model_vars", "V inq_varid "//trim(filename))
- call nc_check(nf90_put_var(ncFileID,VarID,v,start=(/1,1,1,copyindex,timeindex/)),&
+ call nc_check(nf90_put_var(ncFileID,VarID,data_3d,start=(/1,1,1,copyindex,timeindex/)),&
"nc_write_model_vars", "V put_var "//trim(filename))
+ call vector_to_prog_var(statevec,SSH_index,data_2d)
+ where (data_2d == 0.0_r4) data_2d = NF90_FILL_REAL
call nc_check(NF90_inq_varid(ncFileID, "SSH", VarID), &
"nc_write_model_vars", "SSH inq_varid "//trim(filename))
- call nc_check(nf90_put_var(ncFileID,VarID,ssh,start=(/1,1,copyindex,timeindex/)),&
+ call nc_check(nf90_put_var(ncFileID,VarID,data_2d,start=(/1,1,copyindex,timeindex/)),&
"nc_write_model_vars", "SSH put_var "//trim(filename))
- deallocate(s,t,u,v,ssh)
-
endif
!-------------------------------------------------------------------------------
@@ -1936,6 +1828,8 @@
end subroutine write_meta
+
+
subroutine read_2d_snapshot(fbase, x, timestep, vartype)
!------------------------------------------------------------------
!
@@ -1945,14 +1839,14 @@
!
! As it stands now ... the .data files appear to be big-endian.
-character(len=*), intent(in) :: fbase
-real(r4), allocatable, dimension(:,:), intent(out) :: x
-integer, intent(out) :: timestep
-character(len=*), optional, intent(in) :: vartype
+character(len=*), intent(in) :: fbase
+real(r4), dimension(:,:), intent(out) :: x
+integer, intent(out) :: timestep
+character(len=*), optional, intent(in) :: vartype
character(len=128) :: metafilename, datafilename
type(MIT_meta_type):: metadata
-integer :: iunit, io, indx
+integer :: iunit, io
integer :: reclen
if (present(vartype)) then
@@ -1963,8 +1857,10 @@
datafilename = trim(fbase)//'.data'
endif
-metadata = read_meta(fbase,vartype)
+! If the companion ".meta" file exists, it is used.
+! Otherwise, the namelist variables are used.
+metadata = read_meta(fbase,vartype)
timestep = metadata%timeStepNumber
write(*,*)'nDims ',metadata%nDims
@@ -1984,14 +1880,25 @@
call error_handler(E_ERR,'model_mod:read_2d_snapshot',msgstring,source,revision,revdate)
endif
-! good to go
+! Check dimensions
-allocate(x(metadata%dimList(1), metadata%dimList(2)))
+if (size(x, 1) /= metadata%dimList(1) ) then
+ write(msgstring,*)trim(metafilename),' dim 1 does not match delX grid size from namelist'
+ call error_handler(E_MSG,"model_mod:read_2d_snapshot",msgstring,source,revision,revdate)
+ write(msgstring,*)'expected ',size(x,1),' got ',metadata%dimList(1)
+ call error_handler(E_ERR,"model_mod:read_2d_snapshot",msgstring,source,revision,revdate)
+endif
-write(*,*)'shape ',shape(x)
+if (size(x, 2) /= metadata%dimList(2)) then
+ write(msgstring,*)trim(metafilename),' dim 2 does not match delY grid size from namelist'
+ call error_handler(E_MSG,"model_mod:read_2d_snapshot",msgstring,source,revision,revdate)
+ write(msgstring,*)'expected ',size(x,2),' got ',metadata%dimList(2)
+ call error_handler(E_ERR,"model_mod:read_2d_snapshot",msgstring,source,revision,revdate)
+endif
-!! FIXME: on the ibms this must be times the item size;
-!! on the ifort compiler it is just the number of items
+!! FIXME: This is horrid ... 'reclen' is machine-dependent.
+!! IBM XLF -- item_size_direct_read == 4 (number of bytes)
+!! IFORT -- item_size_direct_read == 1 (number of 32bit words)
reclen = product(shape(x)) * item_size_direct_read
! Get next available unit number, read file.
@@ -2003,17 +1910,14 @@
call error_handler(E_ERR,'model_mod:read_2d_snapshot',msgstring,source,revision,revdate)
endif
-print *, 'ready to read 2d snapshot, reclen = ', reclen
read(iunit, rec=1, iostat = io) x
if (io /= 0) then
-print *, 'failed read, iostat = ', io
write(msgstring,*) 'unable to read snapshot file ', trim(datafilename)
- call error_handler(E_ERR,'read_2d_snapshot',msgstring,source,revision,revdate)
+ call error_handler(E_ERR,'model_mod:read_2d_snapshot',msgstring,source,revision,revdate)
endif
close(iunit)
-
end subroutine read_2d_snapshot
@@ -2027,14 +1931,14 @@
!
! As it stands now ... the .data files appear to be big-endian.
-character(len=*), intent(in) :: fbase
-real(r4), allocatable, dimension(:,:,:), intent(out) :: x
-integer, intent(out) :: timestep
-character(len=*), optional, intent(in) :: vartype
+character(len=*), intent(in) :: fbase
+real(r4), dimension(:,:,:), intent(out) :: x
+integer, intent(out) :: timestep
+character(len=*), optional, intent(in) :: vartype
character(len=128) :: metafilename, datafilename
type(MIT_meta_type):: metadata
-integer :: iunit, io, indx
+integer :: iunit, io
integer :: reclen
if (present(vartype)) then
@@ -2045,8 +1949,10 @@
datafilename = trim(fbase)//'.data'
endif
-metadata = read_meta(fbase,vartype)
+! If the companion ".meta" file exists, it is used.
+! Otherwise, the namelist variables are used.
+metadata = read_meta(fbase,vartype)
timestep = metadata%timeStepNumber
! check to make sure storage modes match
@@ -2063,10 +1969,29 @@
call error_handler(E_ERR,'model_mod:read_3d_snapshot',msgstring,source,revision,revdate)
endif
-! good to go
+! Check dimensions
-allocate(x(metadata%dimList(1), metadata%dimList(2), metadata%dimList(3)))
+if (size(x, 1) /= metadata%dimList(1) ) then
+ write(msgstring,*)trim(metafilename),' dim 1 does not match delX grid size from namelist'
+ call error_handler(E_MSG,"model_mod:read_3d_snapshot",msgstring,source,revision,revdate)
+ write(msgstring,*)'expected ',size(x,1),' got ',metadata%dimList(1)
+ call error_handler(E_ERR,"model_mod:read_3d_snapshot",msgstring,source,revision,revdate)
+endif
+if (size(x, 2) /= metadata%dimList(2)) then
+ write(msgstring,*)trim(metafilename),' dim 2 does not match delY grid size from namelist'
+ call error_handler(E_MSG,"model_mod:read_3d_snapshot",msgstring,source,revision,revdate)
+ write(msgstring,*)'expected ',size(x,2),' got ',metadata%dimList(2)
+ call error_handler(E_ERR,"model_mod:read_3d_snapshot",msgstring,source,revision,revdate)
+endif
+
+if (size(x, 3) /= metadata%dimList(3)) then
+ write(msgstring,*)trim(metafilename),' dim 3 does not match delZ grid size from namelist'
+ call error_handler(E_MSG,"model_mod:read_3d_snapshot",msgstring,source,revision,revdate)
+ write(msgstring,*)'expected ',size(x,3),' got ',metadata%dimList(3)
+ call error_handler(E_ERR,"model_mod:read_3d_snapshot",msgstring,source,revision,revdate)
+endif
+
reclen = product(shape(x)) * item_size_direct_read
! Get next available unit number, read file.
@@ -2081,7 +2006,7 @@
read(iunit, rec=1, iostat = io) x
if (io /= 0) then
write(msgstring,*) 'unable to read snapshot file ', trim(datafilename)
- call error_handler(E_ERR,'read_3d_snapshot',msgstring,source,revision,revdate)
+ call error_handler(E_ERR,'model_mod:read_3d_snapshot',msgstring,source,revision,revdate)
endif
close(iunit)
@@ -2090,21 +2015,6 @@
-subroutine drop_3d_snapshot(x)
-!------------------------------------------------------------------
-!
-real(r4), allocatable, intent(inout) :: x(:,:,:)
-deallocate(x)
-end subroutine drop_3d_snapshot
-
-subroutine drop_2d_snapshot(x)
-!------------------------------------------------------------------
-!
-real(r4), allocatable, intent(inout) :: x(:,:)
-deallocate(x)
-end subroutine drop_2d_snapshot
-
-
subroutine write_2d_snapshot(x, fbase)
!------------------------------------------------------------------
!
@@ -2119,7 +2029,7 @@
character(len=128) :: metafilename, datafilename
type(MIT_meta_type) :: metadata
-integer :: iunit, io, indx
+integer :: iunit, io
integer :: reclen
metafilename = trim(fbase)//'.meta'
@@ -2169,7 +2079,7 @@
character(len=128) :: metafilename, datafilename
type(MIT_meta_type) :: metadata
-integer :: iunit, io, indx
+integer :: iunit, io
integer :: reclen
metafilename = trim(fbase)//'.meta'
@@ -2215,20 +2125,12 @@
real(r8), intent(inout) :: state_vector(:)
! temp space to hold data while we are reading it
-real(r4), allocatable :: data_2d_array(:,:), data_3d_array(:,:,:)
+real(r4) :: data_2d_array(Nx,Ny), data_3d_array(Nx,Ny,Nz)
integer :: i, j, k, l, indx, timestepcount_out
! These must be a fixed number and in a fixed order.
-character(len=128) :: names_3d(n3dfields), names_2d(n2dfields)
character(len=128) :: prefixstring
-names_3d(S_index) = 'S'
-names_3d(T_index) = 'T'
-names_3d(U_index) = 'U'
-names_3d(V_index) = 'V'
-names_2d(1) = 'SSH'
-
-
! start counting up and filling the state vector
! one item at a time, linearizing the 3d arrays into a single
! 1d list of numbers.
@@ -2242,40 +2144,40 @@
! e.g. S.Prior.0000000672.data
! S.Prior.0000000672.meta
write(prefixstring, "(A,I10.10)") &
- trim(names_3d(l))//'.'//trim(basename)//'.', timestepcount
+ trim(progvarnames(l))//'.'//trim(basename)//'.', timestepcount
call read_snapshot(prefixstring, data_3d_array, timestepcount_out)
do k = 1, Nz
- do j = 1, Ny
- do i = 1, Nx
- state_vector(indx) = data_3d_array(i, j, k)
- indx = indx + 1
- enddo
- enddo
+ do j = 1, Ny
+ do i = 1, Nx
+ state_vector(indx) = data_3d_array(i, j, k)
+ indx = indx + 1
enddo
- call drop_snapshot(data_3d_array)
+ enddo
+ enddo
enddo
-! and finally, SSH
-do l=1, n2dfields
+! and finally, SSH (and any other 2d fields)
+do l=(n3dfields+1), (n3dfields+n2dfields)
write(prefixstring, "(A,I10.10)") &
- trim(names_2d(l))//'.'//trim(basename)//'.', timestepcount
+ trim(progvarnames(l))//'.'//trim(basename)//'.', timestepcount
call read_snapshot(prefixstring, data_2d_array, timestepcount_out)
do j = 1, Ny
- do i = 1, Nx
- state_vector(indx) = data_2d_array(i, j)
- indx = indx + 1
- enddo
+ do i = 1, Nx
+ state_vector(indx) = data_2d_array(i, j)
+ indx = indx + 1
enddo
- call drop_snapshot(data_2d_array)
+ enddo
enddo
end subroutine snapshot_files_to_sv
+
+
subroutine sv_to_snapshot_files(state_vector, basename, timestepcount)
!------------------------------------------------------------------
!
@@ -2284,20 +2186,12 @@
integer, intent(in) :: timestepcount
! temp space to hold data while we are writing it
-real(r4), allocatable :: data_2d_array(:,:), data_3d_array(:,:,:)
+real(r4) :: data_2d_array(Nx,Ny), data_3d_array(Nx,Ny,Nz)
integer :: i, j, k, l, indx, timestepcount_out
! These must be a fixed number and in a fixed order.
-character(len=128) :: names_3d(n3dfields), names_2d(n2dfields)
character(len=128) :: prefixstring
-names_3d(S_index) = 'S'
-names_3d(T_index) = 'T'
-names_3d(U_index) = 'U'
-names_3d(V_index) = 'V'
-names_2d(1) = 'SSH'
-
-
! start counting up and filling the state vector
! one item at a time, linearizing the 3d arrays into a single
! 1d list of numbers.
@@ -2311,49 +2205,46 @@
! e.g. S.Prior.0000000672.data
! S.Prior.0000000672.meta
write(prefixstring, "(A,I10.10)") &
- trim(names_3d(l))//'.'//trim(basename)//'.', timestepcount
+ trim(progvarnames(l))//'.'//trim(basename)//'.', timestepcount
- allocate(data_3d_array(Nx, Ny, Nz))
do k = 1, Nz
- do j = 1, Ny
- do i = 1, Nx
- data_3d_array(i, j, k) = state_vector(indx)
- indx = indx + 1
- enddo
- enddo
+ do j = 1, Ny
+ do i = 1, Nx
+ data_3d_array(i, j, k) = state_vector(indx)
+ indx = indx + 1
enddo
+ enddo
+ enddo
call write_snapshot(data_3d_array, prefixstring)
- deallocate(data_3d_array)
enddo
-! and finally, SSH
-do l=1, n2dfields
+! and finally, SSH (and any other 2d fields)
+do l=(n3dfields+1), (n3dfields+n2dfields)
write(prefixstring, "(A,I10.10)") &
- trim(names_2d(l))//'.'//trim(basename)//'.', timestepcount
+ trim(progvarnames(l))//'.'//trim(basename)//'.', timestepcount
- allocate(data_2d_array(Nx, Ny))
do j = 1, Ny
- do i = 1, Nx
- data_2d_array(i, j) = state_vector(indx)
- indx = indx + 1
- enddo
+ do i = 1, Nx
+ data_2d_array(i, j) = state_vector(indx)
+ indx = indx + 1
enddo
+ enddo
call write_snapshot(data_2d_array, prefixstring)
- deallocate(data_2d_array)
enddo
end subroutine sv_to_snapshot_files
+
subroutine prog_var_to_vector(s,t,u,v,ssh,x)
!------------------------------------------------------------------
!
-real(r4), dimension(:,:,:), intent(in) :: s,t,u,v
-real(r4), dimension(:,:), intent(in) :: ssh
-real(r8), dimension(:), intent(out) :: x
+real(r4), dimension(Nx,Ny,Nz), intent(in) :: s,t,u,v
+real(r4), dimension(Nx,Ny), intent(in) :: ssh
+real(r8), dimension(:), intent(out) :: x
integer :: i,j,k,ii
@@ -2450,80 +2341,89 @@
end subroutine prog_var_to_vector
-subroutine vector_to_prog_var(x,s,t,u,v,ssh)
+
+subroutine vector_to_2d_prog_var(x, varindex, data_2d_array)
!------------------------------------------------------------------
!
-real(r8), intent(in) :: x(:)
-real(r4), allocatable, dimension(:,:,:), intent(out) :: s,t,u,v
-real(r4), allocatable, dimension(:,:), intent(out) :: ssh
+real(r8), dimension(:), intent(in) :: x
+integer, intent(in) :: varindex
+real(r4), dimension(:,:), intent(out) :: data_2d_array
-integer :: i,j,k,ii
+integer :: i,j,ii
+integer :: dim1,dim2
+character(len=128) :: varname
-allocate( s(Nx,Ny,Nz))
-allocate( t(Nx,Ny,Nz))
-allocate( u(Nx,Ny,Nz))
-allocate( v(Nx,Ny,Nz))
-allocate(ssh(Nx,Ny))
+dim1 = size(data_2d_array,1)
+dim2 = size(data_2d_array,2)
-ii = 0
+varname = progvarnames(varindex)
-! Salinity
-do k = 1,Nz ! vertical
-do j = 1,Ny ! latitudes
-do i = 1,Nx ! longitudes
- ii = ii + 1
- s(i,j,k) = x(ii)
-enddo
-enddo
-enddo
+if (dim1 /= Nx) then
+ write(msgstring,*)trim(varname),' 2d array dim 1 ',dim1,' /= ',Nx
+ call error_handler(E_ERR,'model_mod:vector_to_prog_var',msgstring,source,revision,revdate)
+endif
+if (dim2 /= Ny) then
+ write(msgstring,*)trim(varname),' 2d array dim 2 ',dim2,' /= ',Ny
+ call error_handler(E_ERR,'model_mod:vector_to_prog_var',msgstring,source,revision,revdate)
+endif
-! Temperature
-do k = 1,Nz ! vertical
-do j = 1,Ny ! latitudes
-do i = 1,Nx ! longitudes
- ii = ii + 1
- t(i,j,k) = x(ii)
-enddo
-enddo
-enddo
+ii = start_index(varindex)
-! E-W
-do k = 1,Nz ! vertical
do j = 1,Ny ! latitudes
do i = 1,Nx ! longitudes
+ data_2d_array(i,j) = x(ii)
ii = ii + 1
- u(i,j,k) = x(ii)
enddo
enddo
-enddo
-! N-S
+end subroutine vector_to_2d_prog_var
+
+
+
+subroutine vector_to_3d_prog_var(x, varindex, data_3d_array)
+!------------------------------------------------------------------
+!
+real(r8), dimension(:), intent(in) :: x
+integer, intent(in) :: varindex
+real(r4), dimension(:,:,:), intent(out) :: data_3d_array
+
+integer :: i,j,k,ii
+integer :: dim1,dim2,dim3
+character(len=128) :: varname
+
+dim1 = size(data_3d_array,1)
+dim2 = size(data_3d_array,2)
+dim3 = size(data_3d_array,3)
+
+varname = progvarnames(varindex)
+
+if (dim1 /= Nx) then
+ write(msgstring,*)trim(varname),' 3d array dim 1 ',dim1,' /= ',Nx
+ call error_handler(E_ERR,'model_mod:vector_to_prog_var',msgstring,source,revision,revdate)
+endif
+if (dim2 /= Ny) then
+ write(msgstring,*)trim(varname),' 3d array dim 2 ',dim2,' /= ',Ny
+ call error_handler(E_ERR,'model_mod:vector_to_prog_var',msgstring,source,revision,revdate)
+endif
+if (dim3 /= Nz) then
+ write(msgstring,*)trim(varname),' 3d array dim 3 ',dim3,' /= ',Nz
+ call error_handler(E_ERR,'model_mod:vector_to_prog_var',msgstring,source,revision,revdate)
+endif
+
+ii = start_index(varindex)
+
do k = 1,Nz ! vertical
do j = 1,Ny ! latitudes
do i = 1,Nx ! longitudes
+ data_3d_array(i,j,k) = x(ii)
ii = ii + 1
- v(i,j,k) = x(ii)
enddo
enddo
enddo
-! Sea Surface Height
-do j = 1,Ny ! latitudes
-do i = 1,Nx ! longitudes
- ii = ii + 1
- ssh(i,j) = x(ii)
-enddo
-enddo
+end subroutine vector_to_3d_prog_var
-if (ii /= get_model_size()) then
- write(msgstring,*)'data size ',ii,' /= ',get_model_size(),' model size'
- call error_handler(E_ERR,'model_mod:vector_to_prog_var', &
- msgstring,source,revision,revdate)
-endif
-end subroutine vector_to_prog_var
-
-
!===================================================================
! End of model_mod
!===================================================================
More information about the Dart-dev
mailing list