[Dart-dev] [3945] DART/trunk/models/POP: Everything compiles on coral with Intel 10.1.
nancy at ucar.edu
nancy at ucar.edu
Thu Jun 25 17:06:32 MDT 2009
An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/dart-dev/attachments/20090625/8e986f23/attachment-0001.html
-------------- next part --------------
Modified: DART/trunk/models/POP/dart_to_pop.f90
===================================================================
--- DART/trunk/models/POP/dart_to_pop.f90 2009-06-24 22:12:20 UTC (rev 3944)
+++ DART/trunk/models/POP/dart_to_pop.f90 2009-06-25 23:06:31 UTC (rev 3945)
@@ -9,9 +9,10 @@
! purpose: interface between DART and the POP model
!
! method: Read DART state vector (in file 'assim_model_state_ic') and
-! write out POP "snapshot" files.
+! overwrite values in a POP 'restart' file.
+! Must do something with the advance-to time.
!
-! author: Tim Hoar 5Apr08
+! author: Tim Hoar 25 Jun 09
!
!----------------------------------------------------------------------
@@ -26,9 +27,8 @@
initialize_utilities, finalize_utilities, &
find_namelist_in_file, check_namelist_read, &
logfileunit
-use model_mod, only : sv_to_snapshot_files, static_init_model, &
- get_model_size, DARTtime_to_POPtime, &
- get_model_time_step, write_data_namelistfile, &
+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, &
@@ -42,23 +42,14 @@
revision = "$Revision$", &
revdate = "$Date$"
-character (len = 128) :: file_in = 'assim_model_state_ic'
-
!------------------------------------------------------------------
-! The time manager namelist variables
-! some types/etc come from <POPsource>/pkg/cal/cal.h
-! some useful insight from cal_set.F, cal_convdate.F
-!
-! startDate_1 (integer) yyyymmdd "start date of the integration"
-! startDate_2 (integer) hhmmss
+! The namelist variables
!------------------------------------------------------------------
-character(len=9) :: TheCalendar = 'gregorian'
-integer :: startDate_1 = 19530101
-integer :: startDate_2 = 60000
-logical :: calendarDumps = .false.
+character (len = 128) :: dart_to_pop_input_file = 'assim_model_state_ic'
+character (len = 128) :: dart_to_pop_restart_file = 'my_pop_restart_file'
-NAMELIST /CAL_NML/ TheCalendar, startDate_1, startDate_2, calendarDumps
+namelist /dart_to_pop_nml/ dart_to_pop_input_file, dart_to_pop_restart_file
!----------------------------------------------------------------------
@@ -71,46 +62,50 @@
!----------------------------------------------------------------------
call initialize_utilities('dart_to_pop')
-call static_init_model()
-! POP calendar information. The namelist is already read in
-! static_init_model(), so no further bulletproofing is needed here.
-call find_namelist_in_file("data.cal", "CAL_NML", iunit)
-read(iunit, nml = CAL_NML, iostat = io)
-call check_namelist_read(iunit, io, "CAL_NML")
+! Call model_mod:static_init_model(), which reads the namelists
+! to set calendar type, grid sizes, etc.
+call static_init_model()
+
x_size = get_model_size()
allocate(statevector(x_size))
+! Read the namelist to get the input and output filenames.
+
+call find_namelist_in_file("input.nml", "dart_to_pop_nml", iunit)
+read(iunit, nml = dart_to_pop_nml, iostat = io)
+call check_namelist_read(iunit, io, "dart_to_pop_nml")
+
!----------------------------------------------------------------------
! Reads the valid time, the state, and the target time.
!----------------------------------------------------------------------
-iunit = open_restart_read(file_in)
+iunit = open_restart_read(dart_to_pop_input_file)
call aread_state_restart(model_time, statevector, iunit, adv_to_time)
call close_restart(iunit)
!----------------------------------------------------------------------
-! update the CAL_NML variables so we can rewrite that namelist.
-! Ultimately, we want to keep data:&PARM03:startTime = 0.,
+! update the current POP state vector
!----------------------------------------------------------------------
-call DARTtime_to_POPtime( model_time, startDate_1, startDate_2)
-call sv_to_snapshot_files(statevector, startDate_1, startDate_2)
+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)
+!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 seconds.
+! convert the adv_to_time to the appropriate number of POP
+! time_manager_nml: stop_option, stop_count increments
!----------------------------------------------------------------------
model_timestep = get_model_time_step()
offset = adv_to_time - model_time
-call set_model_end_time(offset)
-call write_data_namelistfile()
+!call set_model_end_time(offset)
+!call write_data_namelistfile()
!----------------------------------------------------------------------
! Log what we think we're doing, and exit.
Added: DART/trunk/models/POP/matlab/Check_pop_to_dart.m
===================================================================
--- DART/trunk/models/POP/matlab/Check_pop_to_dart.m (rev 0)
+++ DART/trunk/models/POP/matlab/Check_pop_to_dart.m 2009-06-25 23:06:31 UTC (rev 3945)
@@ -0,0 +1,83 @@
+% Check_pop_to_dart
+%
+% data.cal holds the starting time information
+%
+
+popdir = '/ptmp/thoar/POP/job_40705';
+popfile = 'pop.r.x1A.00000102';
+dartfile = 'assim_model_state_ud';
+
+fname = sprintf('%s/%s',popdir,popfile);
+
+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');
+
+modelsize = prod(size(S)) + prod(size(T)) + ...
+ prod(size(U)) + prod(size(V)) + prod(size(SSH));
+
+[nx ny nz] = size(S);
+
+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');
+
+fprintf('year month day hour minute second %d %d %d %d %d %d\n', ...
+ iyear,imonth,iday,ihour,iminute,isecond);
+
+% Get the dart equivalent
+
+fname = sprintf('%s/%s',popdir,dartfile);
+fid = fopen(fname,'rb','ieee-le');
+trec1 = fread(fid,1,'int32');
+seconds = fread(fid,1,'int32');
+days = fread(fid,1,'int32');
+trecN = fread(fid,1,'int32');
+
+if (trec1 ~= trecN)
+ error('first record mismatch')
+end
+
+rec1 = fread(fid, 1,'int32');
+datmat = fread(fid,modelsize,'float64');
+recN = fread(fid, 1,'int32');
+fclose(fid);
+
+if (rec1 ~= recN)
+ error('second record mismatch')
+end
+
+disp(sprintf('time tag is days/seconds : %d %d',days,seconds))
+
+offset = prod(size(S)); ind1 = 1; ind2 = offset;
+myS = reshape(datmat(ind1:ind2),size(S));
+d = myS - S;
+disp(sprintf('S diffs are %f %f',min(d(:)),max(d(:))))
+
+offset = prod(size(T)); ind1 = ind2+1; ind2 = ind1 + offset - 1;
+myT = reshape(datmat(ind1:ind2),size(T));
+d = myT - T;
+disp(sprintf('T diffs are %f %f',min(d(:)),max(d(:))))
+
+offset = prod(size(U)); ind1 = ind2+1; ind2 = ind1 + offset - 1;
+myU = reshape(datmat(ind1:ind2),size(U));
+d = myU - U;
+disp(sprintf('U diffs are %f %f',min(d(:)),max(d(:))))
+
+offset = prod(size(V)); ind1 = ind2+1; ind2 = ind1 + offset - 1;
+myV = reshape(datmat(ind1:ind2),size(V));
+d = myV - V;
+disp(sprintf('V diffs are %f %f',min(d(:)),max(d(:))))
+
+offset = prod(size(SSH)); ind1 = ind2+1; ind2 = ind1 + offset - 1;
+mySSH = reshape(datmat(ind1:ind2),size(SSH));
+d = mySSH - SSH;
+disp(sprintf('SSH diffs are %f %f',min(d(:)),max(d(:))))
+
+clear fid fname iyear imonth iday ihour iminute isecond
+clear trec1 trecN rec1 recN offset ind1 ind2 days seconds
Property changes on: DART/trunk/models/POP/matlab/Check_pop_to_dart.m
___________________________________________________________________
Name: svn:mime-type
+ text/x-matlab
Name: svn:keywords
+ Date Revision Author HeadURL Id
Name: svn:eol-style
+ native
Modified: DART/trunk/models/POP/model_mod.f90
===================================================================
--- DART/trunk/models/POP/model_mod.f90 2009-06-24 22:12:20 UTC (rev 3944)
+++ DART/trunk/models/POP/model_mod.f90 2009-06-25 23:06:31 UTC (rev 3945)
@@ -14,9 +14,9 @@
! This is the interface between the POP ocean model and DART.
! Modules that are absolutely required for use are listed
-use types_mod, only : r4, r8, SECPERDAY
+use types_mod, only : r4, r8, SECPERDAY, MISSING_R8
use time_manager_mod, only : time_type, set_time, set_date, get_date, get_time, &
- set_calendar_type, GREGORIAN, print_time, print_date, &
+ set_calendar_type, print_time, print_date, &
operator(*), operator(+), operator(-), &
operator(>), operator(<), operator(/), &
operator(/=), operator(<=)
@@ -33,6 +33,8 @@
use mpi_utilities_mod, only: my_task_id
use random_seq_mod, only : random_seq_type, init_random_seq, random_gaussian
+use typesizes
+use netcdf
implicit none
private
@@ -58,13 +60,8 @@
! generally useful routines for various support purposes.
! the interfaces here can be changed as appropriate.
-public :: prog_var_to_vector, vector_to_prog_var, &
- POP_meta_type, read_meta, write_meta, &
- read_snapshot, write_snapshot, get_gridsize, &
- write_data_namelistfile, set_model_end_time, &
- snapshot_files_to_sv, sv_to_snapshot_files, &
- timestep_to_DARTtime, DARTtime_to_POPtime, &
- DARTtime_to_timestepindex
+public :: POP_meta_type, get_gridsize, set_model_end_time, &
+ restart_file_to_sv, sv_to_restart_file
! version controlled file description for error handling, do not edit
character(len=128), parameter :: &
@@ -72,175 +69,112 @@
revision = "$Revision$", &
revdate = "$Date$"
-character(len=129) :: msgstring
+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
-
-!! FIXME: This is horrid ... 'reclen' is compiler-dependent.
-!! IBM XLF -- item_size_direct_access == 4,8
-!! IFORT -- needs compile switch "-assume byterecl"
-integer, parameter :: item_size_direct_access = 4
-
!------------------------------------------------------------------
!
-! POP namelist section: we want to share the 'data' namelist file
-! with the model, so we must declare all possible namelist entries to
-! avoid getting an error from a valid namelist file. Most of these
-! values are unused in this model_mod code; only a few are needed and
-! those are indicated in comments below.
+! POP namelist section:
+! 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.
+! if the grid is global, or at least wraps around the sphere in
+! longitude, set longitude_wrap to be true.
!
!------------------------------------------------------------------
! The time manager namelist variables
-! some types/etc come from <POPsource>/pkg/cal/cal.h
-! some useful insight from cal_set.F, cal_convdate.F
-!
-! startDate_1 (integer) yyyymmdd "start date of the integration"
-! startDate_2 (integer) hhmmss
!------------------------------------------------------------------
-character(len=9) :: TheCalendar = 'gregorian'
-! integration start date follows: yyyymmddhhmmss
-integer :: startDate_1 = 19530101
-integer :: startDate_2 = 60000
-logical :: calendarDumps = .false.
+character(len=128) :: horiz_grid_input_file = 'no_horiz_grid_input_file'
+character(len=128) :: topography_input_file = 'no_topography_input_file'
+character(len=128) :: vert_grid_input_file = 'no_vert_grid_input_file'
+character(len=128) :: dart_restart_file = 'no_dart_restart_file'
+character(len=128) :: pop_data_file = 'no_pop_data_file'
+logical :: longitude_wrap = .true.
-NAMELIST /CAL_NML/ TheCalendar, startDate_1, startDate_2, calendarDumps
+!integer :: nblocks - model_mod doesn't need this
+! but advance_model will. not sure this
+! helps in this code, since the script
+! needs the info, not the filter executable.
-! FIXME: these namelists should probably be in a separate file, and only
-! the actual values needed should be made public, so this isn't so messy.
+! here are what we can get from the 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 topog file:
+! integer, dimension(:,:), allocatable :: &
+! KMT ! k index of deepest grid cell on T grid
+!
+! the vert grid file is ascii, with 3 columns/line:
+! FIXME: confirm this:
+! depth in cm? vert cell center? distance between centers or
+! layer tops/bottoms?
-! must match the value in EEPARAMS.h
-integer, parameter :: MAX_LEN_FNAM = 512
-! these come from SIZE.h and can vary.
-! should they be namelist controlled?
-integer, parameter :: max_nx = 1024
-integer, parameter :: max_ny = 1024
-integer, parameter :: max_nz = 512
-integer, parameter :: max_nr = 512
+! other 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
-!-- 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
+character(len=32):: calendar = 'noleap'
-!-- Time stepping parameters variable declarations
-real(r8) :: pickupSuff, &
- deltaT, deltaTClock, deltaTmom, &
- deltaTtracer, dTtracerLev(max_nr), deltaTfreesurf, &
- abEps, alph_AB, beta_AB, &
- 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
-integer :: nIter0, nTimeSteps, nEndIter, momForcingOutAB, tracForcingOutAB
-logical :: forcing_In_AB, &
- momDissip_In_AB, doAB_onGtGs, &
- startFromPickupAB2
+namelist /model_nml/ &
+ horiz_grid_input_file, &
+ topography_input_file, &
+ vert_grid_input_file, &
+ dart_restart_file, &
+ pop_data_file, &
+ longitude_wrap, &
+ output_state_vector, &
+ assimilation_period_days, & ! for now, this is the timestep
+ assimilation_period_seconds, &
+ model_perturbation_amplitude,&
+ calendar
-!-- Gridding parameters variable declarations
-logical :: usingCartesianGrid, usingCylindricalGrid, &
- usingSphericalPolarGrid, usingCurvilinearGrid, &
- deepAtmosphere
-real(r8) :: dxSpacing, dySpacing, delX(max_nx), delY(max_ny), &
- phiMin, thetaMin, rSphere, &
- Ro_SeaLevel, delZ(max_nz), delP, delR(max_nr), delRc(max_nr+1), &
- rkFac, groundAtK1
-character(len=MAX_LEN_FNAM) :: delXFile, delYFile, &
- delRFile, delRcFile, &
- horizGridFile
-
-!-- Input files variable declarations
-character(len=MAX_LEN_FNAM) :: &
- 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
-
-!-- 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, Eta
+! The DART state vector (control vector) will consist of: S, T, U, V, SHGT
! (Salinity, Temperature, U velocity, V velocity, Sea Surface Height).
-! S, T are 3D arrays, located at cell centers. U is staggered in X
-! and V is staggered in Y (meaning the points are located on the cell
-! faces) but the grids are offset by half a cell, so there are actually
-! the same number of points in each grid.
-! Eta is a 2D field (X,Y only). The Z direction is downward.
+! S, T are 3D arrays, located at cell centers. U,V are at grid cell corners.
+! SHGT 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
-integer, parameter :: S_index = 1
-integer, parameter :: T_index = 2
-integer, parameter :: U_index = 3
-integer, parameter :: V_index = 4
-integer, parameter :: Eta_index = 5
-
! (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 ','Eta'/)
+character(len=128) :: progvarnames(nfields) = (/'SALT','TEMP','UVEL','VVEL','SHGT'/)
+integer, parameter :: S_index = 1
+integer, parameter :: T_index = 2
+integer, parameter :: U_index = 3
+integer, parameter :: V_index = 4
+integer, parameter :: SHGT_index = 5
+
integer :: start_index(nfields)
! Grid parameters - the values will be read from a
@@ -251,32 +185,16 @@
! locations of cell centers (C) and edges (G) for each axis.
real(r8), allocatable :: XC(:), XG(:), YC(:), YG(:), ZC(:), ZG(:)
-! location information - these grids can either be regularly
-! spaced or the spacing along each axis can vary.
+! integer, lowest valid cell number in the vertical
+integer, allocatable :: KMT(:, :)
-!real(r8) :: lat_origin, lon_origin
-!logical :: regular_lat, regular_lon, regular_depth
-!real(r8) :: delta_lat, delta_lon, delta_depth
-!real(r8), allocatable :: lat_grid(:), lon_grid(:), depth_grid(:)
-
+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
-! Skeleton of a model_nml that would be in input.nml
-! This is where dart-related model parms could be set.
-logical :: output_state_vector = .true.
-integer :: assimilation_period_days = 7
-integer :: assimilation_period_seconds = 0
-real(r8) :: model_perturbation_amplitude = 0.2
-
-namelist /model_nml/ output_state_vector, &
- assimilation_period_days, &
- assimilation_period_seconds, &
- model_perturbation_amplitude
-
! /pkg/mdsio/mdsio_write_meta.F writes the .meta files
type POP_meta_type
! private
@@ -288,16 +206,6 @@
integer :: timeStepNumber ! optional
end type POP_meta_type
-INTERFACE read_snapshot
- MODULE PROCEDURE read_2d_snapshot
- MODULE PROCEDURE read_3d_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
@@ -314,27 +222,27 @@
!------------------------------------------------------------------
!
! Called to do one time initialization of the model. In this case,
-! it reads in the grid information and then the model data.
+! it reads in the grid information.
integer :: i, iunit, io
integer :: ss, dd
! The Plan:
!
-! read the standard POP namelist file 'data.cal' for the calendar info
+! read in the grid sizes from the horiz grid file and the vert grid file
+! horiz is netcdf, vert is ascii
+!
+! allocate space, and read in actual grid values
!
-! read the standard POP namelist file 'data' for the
-! time stepping info and the grid info.
+! figure out model timestep. FIXME: from where?
!
-! open the grid data files to get the actual grid coordinates
-!
! Compute the model size.
!
! set the index numbers where the field types change
!
-! set the grid location info
-!
+if ( module_initialized ) return ! only need to do this once.
+
! Print module information to log file and stdout.
call register_module(source, revision, revdate)
@@ -342,7 +250,6 @@
! we'll say we've been initialized pretty dang early.
module_initialized = .true.
-
! Read the DART namelist for this model
call find_namelist_in_file("input.nml", "model_nml", iunit)
read(iunit, nml = model_nml, iostat = io)
@@ -354,45 +261,34 @@
if (do_output()) write( * , nml=model_nml)
! POP calendar information
-call find_namelist_in_file("data.cal", "CAL_NML", iunit)
-read(iunit, nml = CAL_NML, iostat = io)
-call check_namelist_read(iunit, io, "CAL_NML")
-if (index(TheCalendar,'g') > 0 ) then
- call set_calendar_type(GREGORIAN)
-elseif (index(TheCalendar,'G') > 0 )then
- call set_calendar_type(GREGORIAN)
-else
- write(msgstring,*)"namelist data.cal indicates a ",trim(TheCalendar)," calendar."
- call error_handler(E_MSG,"static_init_model", msgstring, source, revision, revdate)
- write(msgstring,*)"only have support for Gregorian"
- call error_handler(E_ERR,"static_init_model", msgstring, source, revision, revdate)
-endif
-if (do_output()) write(*,*)'model_mod:namelist cal_NML',startDate_1,startDate_2
-if (do_output()) write(*,nml=CAL_NML)
+call set_calendar_type(calendar)
-! Time stepping parameters are in PARM03
-call find_namelist_in_file("data", "PARM03", iunit)
-read(iunit, nml = PARM03, iostat = io)
-call check_namelist_read(iunit, io, "PARM03")
+!---------------------------------------------------------------
+! get data dimensions, then allocate space, then open the files
+! and actually fill in the arrays.
-if ((deltaTmom == deltaTtracer) .and. &
- (deltaTmom == deltaTClock ) .and. &
- (deltaTClock == deltaTtracer)) then
- ocean_dynamics_timestep = deltaTmom ! need a time_type version
-else
- write(msgstring,*)"namelist PARM03 has deltaTmom /= deltaTtracer /= deltaTClock"
- call error_handler(E_MSG,"static_init_model", msgstring, source, revision, revdate)
- write(msgstring,*)"values were ",deltaTmom, deltaTtracer, deltaTClock
- call error_handler(E_MSG,"static_init_model", msgstring, source, revision, revdate)
- write(msgstring,*)"At present, DART only supports equal values."
- call error_handler(E_ERR,"static_init_model", msgstring, source, revision, revdate)
-endif
+call get_horiz_grid_dims(Nx, Ny)
+call get_vert_grid_dims(Nz)
-! Define the assimilation period as the model_timestep
-! Ensure model_timestep is multiple of ocean_dynamics_timestep
+! Allocate space for grid variables.
+allocate(XC(Nx), YC(Ny), ZC(Nz))
+allocate(XG(Nx), YG(Ny), ZG(Nz))
+allocate(KMT(Nx, Ny))
-model_time = timestep_to_DARTtime(timestepcount)
+! Fill them in
+call read_horiz_grid(Nx, Ny, XC, XG, YC, YG)
+call read_vert_grid(Nz, ZC, ZG)
+call read_kmt(Nx, Ny, KMT)
+
+call write_grid_netcdf(XG, XC, YG, YC, ZG, ZC, KMT) ! DEBUG only
+
+!---------------------------------------------------------------
+! set the time step from the namelist for now.
+
+!! Define the assimilation period as the model_timestep
+!! Ensure model_timestep is multiple of ocean_dynamics_timestep
+!
model_timestep = set_model_time_step(assimilation_period_seconds, &
assimilation_period_days, &
ocean_dynamics_timestep)
@@ -403,114 +299,29 @@
call error_handler(E_MSG,'static_init_model',msgstring,source,revision,revdate)
if (do_output()) write(logfileunit,*)msgstring
-! Grid-related variables are in PARM04
-delX(:) = 0.0_r4
-delY(:) = 0.0_r4
-delZ(:) = 0.0_r4
-call find_namelist_in_file("data", "PARM04", iunit)
-read(iunit, nml = PARM04, iostat = io)
-call check_namelist_read(iunit, io, "PARM04")
+!---------------------------------------------------------------
+! compute the offsets into the state vector for the start of each
+! different variable type.
-! Input datasets are in PARM05
-call find_namelist_in_file("data", "PARM05", iunit)
-read(iunit, nml = PARM05, iostat = io)
-call check_namelist_read(iunit, io, "PARM05")
-
-! The only way I know to compute the number of
-! levels/lats/lons is to set the default value of delZ to 0.0
-! before reading the namelist. now loop until you get back
-! to zero and that is the end of the list.
-! Not a very satisfying/robust solution ...
-
-Nx = -1
-do i=1, size(delX)
- if (delX(i) == 0.0_r4) then
- Nx = i-1
- exit
- endif
-enddo
-if (Nx == -1) then
- 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
-do i=1, size(delY)
- if (delY(i) == 0.0_r4) then
- Ny = i-1
- exit
- endif
-enddo
-if (Ny == -1) then
- 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
-do i=1, size(delZ)
- if (delZ(i) == 0.0_r4) then
- Nz = i-1
- exit
- endif
-enddo
-if (Nz == -1) then
- 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.
-
-allocate(XC(Nx), YC(Ny), ZC(Nz))
-allocate(XG(Nx), YG(Ny), ZG(Nz))
-
-! XG (the grid edges) and XC (the grid centroids) must be computed.
-
-XG(1) = thetaMin
-XC(1) = thetaMin + 0.5_r8 * delX(1)
-do i=2, Nx
- XG(i) = XG(i-1) + delX(i-1)
- XC(i) = XC(i-1) + 0.5_r8 * delX(i-1) + 0.5_r8 * delX(i)
-enddo
-
-! YG (the grid edges) and YC (the grid centroids) must be computed.
-
-YG(1) = phiMin
-YC(1) = phiMin + 0.5_r8 * delY(1)
-do i=2, Ny
- YG(i) = YG(i-1) + delY(i-1)
- YC(i) = YC(i-1) + 0.5_r8 * delY(i-1) + 0.5_r8 * delY(i)
-enddo
-
-! 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
- 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
-
! record where in the state vector the data type changes
! from one type to another, by computing the starting
! index for each block of data.
-start_index(S_index) = 1
-start_index(T_index) = start_index(S_index) + (Nx * Ny * Nz)
-start_index(U_index) = start_index(T_index) + (Nx * Ny * Nz)
-start_index(V_index) = start_index(U_index) + (Nx * Ny * Nz)
-start_index(Eta_index) = start_index(V_index) + (Nx * Ny * Nz)
+start_index(S_index) = 1
+start_index(T_index) = start_index(S_index) + (Nx * Ny * Nz)
+start_index(U_index) = start_index(T_index) + (Nx * Ny * Nz)
+start_index(V_index) = start_index(U_index) + (Nx * Ny * Nz)
+start_index(SHGT_index) = start_index(V_index) + (Nx * Ny * Nz)
! in spite of the staggering, all grids are the same size
! and offset by half a grid cell. 4 are 3D and 1 is 2D.
! e.g. S,T,U,V = 256 x 225 x 70
-! e.g. Eta = 256 x 225
+! e.g. SHGT = 256 x 225
if (do_output()) write(logfileunit, *) 'Using grid size : '
if (do_output()) write(logfileunit, *) ' Nx, Ny, Nz = ', Nx, Ny, Nz
if (do_output()) write( * , *) 'Using grid size : '
if (do_output()) write( * , *) ' Nx, Ny, Nz = ', Nx, Ny, Nz
-!print *, ' 3d field size: ', n3dfields * (Nx * Ny * Nz)
-!print *, ' 2d field size: ', n2dfields * (Nx * Ny)
+
model_size = (n3dfields * (Nx * Ny * Nz)) + (n2dfields * (Nx * Ny))
if (do_output()) write(*,*) 'model_size = ', model_size
@@ -617,7 +428,7 @@
real(r8), intent(out) :: interp_val
integer, intent(out) :: istatus
-! Model interpolate will interpolate any state variable (S, T, U, V, Eta) to
+! Model interpolate will interpolate any state variable (S, T, U, V, SHGT) to
! the given location given a state vector. The type of the variable being
! interpolated is obs_type since normally this is used to find the expected
! value of an observation at some location. The interpolated value is
@@ -804,7 +615,7 @@
lat_array = yg
call lat_bounds(llat, ny, lat_array, lat_bot, lat_top, lat_fract, lat_status)
else
- ! Eta, T and S are on the YC latitude grid
+ ! SHGT, T and S are on the YC latitude grid
lat_array = yc
call lat_bounds(llat, ny, lat_array, lat_bot, lat_top, lat_fract, lat_status)
endif
@@ -821,7 +632,7 @@
lon_array = xg
call lon_bounds(llon, nx, lon_array, lon_bot, lon_top, lon_fract, lon_status)
else
- ! Eta, T, and S are on the XC grid
+ ! SHGT, T, and S are on the XC grid
lon_array = xc
call lon_bounds(llon, nx, lon_array, lon_bot, lon_top, lon_fract, lon_status)
endif
@@ -989,7 +800,7 @@
end do
! Falling off the end means it's in between. Check for wraparound.
-if(wrap_around) then
+if(longitude_wrap) then
bot = nlons
top = 1
dist_bot = lon_dist(llon, lon_array(bot))
@@ -1161,7 +972,7 @@
var_num = V_index
else
if (present(var_type)) var_type = KIND_SEA_SURFACE_HEIGHT
- var_num = Eta_index
+ var_num = SHGT_index
endif
!print *, 'var num = ', var_num
@@ -1171,7 +982,7 @@
!print *, 'offset = ', offset
-if (var_num == Eta_index) then
+if (var_num == SHGT_index) then
depth = 0.0
depth_index = 1
else
@@ -1200,10 +1011,9 @@
! Does any shutdown and clean-up needed for model. Can be a NULL
! INTERFACE if the model has no need to clean up storage, etc.
-! good style ... perhaps you could deallocate stuff (from static_init_model?).
-! deallocate(state_loc)
+! if ( .not. module_initialized ) call static_init_model
-if ( .not. module_initialized ) call static_init_model
+deallocate(XC, YC, ZC, XG, YG, ZG, KMT)
end subroutine end_model
@@ -1229,9 +1039,6 @@
! NF90_put_var ! provide values for variable
! NF90_CLOSE ! close: save updated netCDF dataset
-use typeSizes
-use netcdf
-
integer, intent(in) :: ncFileID ! netCDF file identifier
integer :: ierr ! return value of function
@@ -1257,7 +1064,7 @@
integer :: XGVarID, XCVarID, YGVarID, YCVarID, ZGVarID, ZCVarID
! for the prognostic variables
-integer :: SVarID, TVarID, UVarID, VVarID, EtaVarID
+integer :: SVarID, TVarID, UVarID, VVarID, SHGTVarID
!----------------------------------------------------------------------
! variables for the namelist output
@@ -1353,7 +1160,7 @@
! Determine shape of most important namelist
!-------------------------------------------------------------------------------
-call find_textfile_dims("data", nlines, linelen)
+call find_textfile_dims("pop_in", nlines, linelen)
allocate(textblock(nlines))
textblock = ''
@@ -1362,11 +1169,11 @@
len = nlines, dimid = nlinesDimID), &
'nc_write_model_atts', 'def_dim nlines ')
-call nc_check(nf90_def_var(ncFileID,name="datanml", xtype=nf90_char, &
+call nc_check(nf90_def_var(ncFileID,name="pop_in", xtype=nf90_char, &
dimids = (/ linelenDimID, nlinesDimID /), varid=nmlVarID), &
- 'nc_write_model_atts', 'def_var datanml')
+ 'nc_write_model_atts', 'def_var pop_in')
call nc_check(nf90_put_att(ncFileID, nmlVarID, "long_name", &
- "contents of data namelist"), 'nc_write_model_atts', 'put_att datanml')
+ "contents of pop_in namelist"), 'nc_write_model_atts', 'put_att pop_in')
!-------------------------------------------------------------------------------
! Here is the extensible part. The simplest scenario is to output the state vector,
@@ -1431,7 +1238,13 @@
! Create the (empty) Coordinate Variables and the Attributes
!----------------------------------------------------------------------------
- ! U Grid Longitudes
+ call nc_check(nf90_def_var(ncFileID,name="POPnml", xtype=nf90_char, &
+ dimids = (/ linelenDimID, nlinesDimID /), varid=nmlVarID), &
+ 'nc_write_model_atts', 'def_var POPnml')
+ call nc_check(nf90_put_att(ncFileID, nmlVarID, "long_name", &
+ "namelist.input contents"), 'nc_write_model_atts', 'put_att POPnml')
+
+ ! U,V Grid Longitudes
call nc_check(nf90_def_var(ncFileID,name="XG",xtype=nf90_real,dimids=XGDimID,varid=XGVarID),&
"nc_write_model_atts", "XG def_var "//trim(filename))
call nc_check(nf90_put_att(ncFileID, XGVarID, "long_name", "longitude grid edges"), &
@@ -1443,7 +1256,7 @@
call nc_check(nf90_put_att(ncFileID, XGVarID, "valid_range", (/ 0.0_r8, 360.0_r8 /)), &
"nc_write_model_atts", "XG valid_range "//trim(filename))
- ! S,T,V,Eta Grid Longitudes
+ ! S,T,SHGT Grid (center) Longitudes
call nc_check(nf90_def_var(ncFileID,name="XC",xtype=nf90_real,dimids=XCDimID,varid=XCVarID),&
"nc_write_model_atts", "XC def_var "//trim(filename))
call nc_check(nf90_put_att(ncFileID, XCVarID, "long_name", "longitude grid centroids"), &
@@ -1455,7 +1268,7 @@
call nc_check(nf90_put_att(ncFileID, XCVarID, "valid_range", (/ 0.0_r8, 360.0_r8 /)), &
"nc_write_model_atts", "XC valid_range "//trim(filename))
- ! V Grid Latitudes
+ ! U,V Grid Latitudes
call nc_check(nf90_def_var(ncFileID,name="YG",xtype=nf90_real,dimids=YGDimID,varid=YGVarID),&
"nc_write_model_atts", "YG def_var "//trim(filename))
call nc_check(nf90_put_att(ncFileID, YGVarID, "long_name", "latitude grid edges"), &
@@ -1467,7 +1280,7 @@
call nc_check(nf90_put_att(ncFileID,YGVarID,"valid_range",(/-90.0_r8,90.0_r8 /)), &
"nc_write_model_atts", "YG valid_range "//trim(filename))
- ! S,T,U,Eta Grid Latitudes
+ ! S,T,SHGT Grid (center) Latitudes
call nc_check(nf90_def_var(ncFileID,name="YC",xtype=nf90_real,dimids=YCDimID,varid=YCVarID), &
"nc_write_model_atts", "YC def_var "//trim(filename))
call nc_check(nf90_put_att(ncFileID, YCVarID, "long_name", "latitude grid centroids"), &
@@ -1513,63 +1326,71 @@
! Create the (empty) Prognostic Variables and the Attributes
!----------------------------------------------------------------------------
- call nc_check(nf90_def_var(ncid=ncFileID, name="S", xtype=nf90_real, &
+ call nc_check(nf90_def_var(ncid=ncFileID, name="SALT", xtype=nf90_real, &
dimids = (/XCDimID,YCDimID,ZCDimID,MemberDimID,unlimitedDimID/),varid=SVarID),&
"nc_write_model_atts", "S def_var "//trim(filename))
call nc_check(nf90_put_att(ncFileID, SVarID, "long_name", "salinity"), &
"nc_write_model_atts", "S long_name "//trim(filename))
+ call nc_check(nf90_put_att(ncFileID, SVarID, "units", "g/g"), &
+ "nc_write_model_atts", "S units "//trim(filename))
+ call nc_check(nf90_put_att(ncFileID, SVarID, "missing_value", NF90_FILL_REAL), &
+ "nc_write_model_atts", "S missing "//trim(filename))
call nc_check(nf90_put_att(ncFileID, SVarID, "_FillValue", NF90_FILL_REAL), &
"nc_write_model_atts", "S fill "//trim(filename))
- call nc_check(nf90_put_att(ncFileID, SVarID, "units", "psu"), &
- "nc_write_model_atts", "S units "//trim(filename))
- call nc_check(nf90_put_att(ncFileID, SVarID, "units_long_name", "practical salinity units"), &
- "nc_write_model_atts", "S units_long_name "//trim(filename))
- call nc_check(nf90_def_var(ncid=ncFileID, name="T", xtype=nf90_real, &
+ call nc_check(nf90_def_var(ncid=ncFileID, name="TEMP", xtype=nf90_real, &
dimids=(/XCDimID,YCDimID,ZCDimID,MemberDimID,unlimitedDimID/),varid=TVarID),&
"nc_write_model_atts", "T def_var "//trim(filename))
- call nc_check(nf90_put_att(ncFileID, TVarID, "long_name", "Temperature"), &
+ call nc_check(nf90_put_att(ncFileID, TVarID, "long_name", "Potential Temperature"), &
"nc_write_model_atts", "T long_name "//trim(filename))
- call nc_check(nf90_put_att(ncFileID, TVarID, "_FillValue", NF90_FILL_REAL), &
- "nc_write_model_atts", "T fill "//trim(filename))
- call nc_check(nf90_put_att(ncFileID, TVarID, "units", "C"), &
+ call nc_check(nf90_put_att(ncFileID, TVarID, "units", "deg C"), &
"nc_write_model_atts", "T units "//trim(filename))
call nc_check(nf90_put_att(ncFileID, TVarID, "units_long_name", "degrees celsius"), &
"nc_write_model_atts", "T units_long_name "//trim(filename))
+ call nc_check(nf90_put_att(ncFileID, TVarID, "missing_value", NF90_FILL_REAL), &
+ "nc_write_model_atts", "T missing "//trim(filename))
+ call nc_check(nf90_put_att(ncFileID, TVarID, "_FillValue", NF90_FILL_REAL), &
@@ Diff output truncated at 40000 characters. @@
More information about the Dart-dev
mailing list