[Dart-dev] DART/trunk Revision: 10738
dart at ucar.edu
dart at ucar.edu
Thu Nov 10 14:30:31 MST 2016
thoar at ucar.edu
2016-11-10 14:30:31 -0700 (Thu, 10 Nov 2016)
458
The assimilation_period_[days,seconds] namelist variables are now used.
Previously they were IGNORED and the pop_in namelist setting of restart_freq
and restart_freq_opt were used to determine the assimilation frequency.
Now - if assimilation_period_[days,seconds] are both <= 0, everything is
backwards compatible. The default values for these are now -1.
Error messages have been improved for illegal combinations of
assimilation_period_[days,seconds]
Modified: DART/trunk/models/POP/dart_pop_mod.f90
===================================================================
--- DART/trunk/models/POP/dart_pop_mod.f90 2016-11-10 18:10:41 UTC (rev 10737)
+++ DART/trunk/models/POP/dart_pop_mod.f90 2016-11-10 21:30:31 UTC (rev 10738)
@@ -33,7 +33,8 @@
character(len=32 ), parameter :: revision = "$Revision$"
character(len=128), parameter :: revdate = "$Date$"
-character(len=256) :: msgstring
+character(len=512) :: string1, string2, string3
+
logical, save :: module_initialized = .false.
character(len=256) :: ic_filename = 'pop.r.nc'
@@ -87,7 +88,7 @@
restart_fmt, leven_odd_on, even_odd_freq, pressure_correction
!------------------------------------------------------------------
-! The POP initial temperature and salinity namelist
+! The POP initial temperature and salinity namelist
!------------------------------------------------------------------
character(len=100) :: init_ts_file ! length consistent with POP
@@ -96,11 +97,11 @@
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_file, init_ts_file_fmt, &
init_ts_outfile, init_ts_outfile_fmt
!------------------------------------------------------------------
-! The POP domain namelist
+! The POP domain namelist
!------------------------------------------------------------------
character(len= 64) :: clinic_distribution_type, tropic_distribution_type
@@ -112,12 +113,12 @@
ew_boundary_type, ns_boundary_type
!------------------------------------------------------------------
-! The POP grid info namelist
+! The POP grid info namelist
!------------------------------------------------------------------
!
! POP grid information comes in several files:
-! horizontal grid lat/lons in one,
-! topography (lowest valid vert level) in another, and
+! horizontal grid lat/lons in one,
+! topography (lowest valid vert level) in another, and
! the vertical grid spacing in a third.
!
!------------------------------------------------------------------
@@ -168,13 +169,13 @@
subroutine initialize_module
-!------------------------------------------------------------------
+
integer :: iunit, io
! Read POP calendar information
! In 'restart' mode, this is primarily the calendar type and 'stop'
-! information. The time attributes of the restart file override
-! the namelist time information.
+! information. The time attributes of the restart file override
+! the namelist time information.
call find_namelist_in_file('pop_in', 'time_manager_nml', iunit)
read(iunit, nml = time_manager_nml, iostat = io)
@@ -207,15 +208,15 @@
! restart_filename = trim(pointer_filename)//'.restart'
!
! if ( .not. file_exist(restart_filename) ) then
-! msgstring = 'pop_in:pointer file '//trim(restart_filename)//' not found'
+! string1 = 'pop_in:pointer file '//trim(restart_filename)//' not found'
! call error_handler(E_ERR,'initialize_module', &
-! msgstring, source, revision, revdate)
+! string1, source, revision, revdate)
! endif
!
! iunit = open_file(restart_filename,'formatted')
! read(iunit,'(A)')ic_filename
!
-! restart_filename = ' '
+! restart_filename = ' '
! write(*,*)'DEBUG ... pointer filename dereferenced to ',trim(ic_filename )
!
!else
@@ -224,9 +225,9 @@
! Make sure we have a pop restart file (for grid dims)
if ( .not. file_exist(ic_filename) ) then
- msgstring = 'pop_in:init_ts_file '//trim(ic_filename)//' not found'
+ string1 = 'pop_in:init_ts_file '//trim(ic_filename)//' not found'
call error_handler(E_ERR,'initialize_module', &
- msgstring, source, revision, revdate)
+ string1, source, revision, revdate)
endif
! Read POP restart information (for model timestepping)
@@ -252,11 +253,11 @@
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.
!
@@ -279,9 +280,9 @@
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)
+ string1 = "unable to find either 'i' or 'nlon' in file "//trim(ic_filename)
+ call error_handler(E_ERR, 'get_horiz_grid_dims', string1, &
+ source,revision,revdate)
endif
endif
@@ -293,8 +294,8 @@
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, &
+ string1 = "unable to find either 'j' or 'nlat' in "//trim(ic_filename)
+ call error_handler(E_ERR, 'get_horiz_grid_dims', string1, &
source,revision,revdate)
endif
endif
@@ -310,11 +311,11 @@
end subroutine get_horiz_grid_dims
+!------------------------------------------------------------------
- subroutine get_vert_grid_dim(Nz)
-!------------------------------------------------------------------
-! 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.
@@ -329,10 +330,12 @@
end subroutine get_vert_grid_dim
-
+!------------------------------------------------------------------
+
+
subroutine get_pop_calendar(calstring)
-!------------------------------------------------------------------
-! the initialize_module ensures that the pop namelists are read and
+
+! the initialize_module ensures that the pop namelists are read and
! the DART time manager gets the pop calendar setting.
!
! Then, the DART time manager is queried to return what it knows ...
@@ -346,38 +349,63 @@
end subroutine get_pop_calendar
+!------------------------------------------------------------------
-function set_model_time_step()
-!------------------------------------------------------------------
+
+function set_model_time_step(seconds,days)
+
! the initialize_module ensures that the pop namelists are read.
! The restart times in the pop_in&restart_nml are used to define
! appropriate assimilation timesteps.
-!
+
+integer, intent(in) :: seconds ! input.nml assimilation_period_seconds
+integer, intent(in) :: days ! input.nml assimilation_period_days
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
+! Check to see if the input.nml:&model_nml:assimilation_period_[seconds,days]
+! are specifying the nominal model output frequency/assimilation interval.
-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_POP uses a bogus value for this
- set_model_time_step = set_time(0, 1) ! (seconds, days)
+if (seconds <= 0 .and. days <= 0) then
+
+ ! 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_POP uses a bogus value for this
+ set_model_time_step = set_time(0, 1) ! (seconds, days)
+ write(string1,*)'WARNING: POP namelist variable "restart_freq_opt" read as "nyear"'
+ write(string2,*)'CESM uses bogus value - using default value of 1 day.'
+ call error_handler(E_MSG,'set_model_time_step', string1, &
+ source, revision, revdate, text2=string2)
+ else
+ write(string1,*)'POP namelist variable "restart_freq_opt" must be "nday" or "nyear"'
+ write(string2,*)'read as "'//trim(restart_freq_opt)//'"'
+ call error_handler(E_ERR,'set_model_time_step', string1, &
+ source, revision, revdate, text2=string2)
+ endif
+
+elseif (seconds < 0 .or. days < 0) then
+ write(string1,*)'Unable to determine the assimilation interval.'
+ write(string2,*)'input.nml:&model_nml:assimilation_period_seconds must be >= 0; is ',seconds
+ write(string3,*)'input.nml:&model_nml:assimilation_period_days must be >= 0; is ',days
+ call error_handler(E_ERR,'set_model_time_step', string1, &
+ source, revision, revdate, text2=string2, text3=string3)
else
- call error_handler(E_ERR,'set_model_time_step', &
- 'restart_freq_opt must be nday', source, revision, revdate)
+ set_model_time_step = set_time(seconds, days) ! (seconds, days)
endif
end function set_model_time_step
+!------------------------------------------------------------------
subroutine write_pop_namelist(model_time, adv_to_time)
-!------------------------------------------------------------------
-!
+
type(time_type), INTENT(IN) :: model_time, adv_to_time
type(time_type) :: offset
@@ -389,8 +417,8 @@
call get_time(offset, secs, days)
if (secs /= 0 ) then
- write(msgstring,*)'adv_to_time has seconds == ',secs,' must be zero'
- call error_handler(E_ERR,'write_pop_namelist', msgstring, source, revision, revdate)
+ write(string1,*)'adv_to_time has seconds == ',secs,' must be zero'
+ call error_handler(E_ERR,'write_pop_namelist', string1, source, revision, revdate)
endif
! call print_date( model_time,'write_pop_namelist:dart model date')
@@ -417,11 +445,11 @@
end subroutine write_pop_namelist
+!------------------------------------------------------------------
- subroutine read_horiz_grid(nx, ny, ULAT, ULON, TLAT, TLON)
-!------------------------------------------------------------------
-! 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
@@ -441,9 +469,9 @@
! Check to see that the file exists.
if ( .not. file_exist(horiz_grid_file) ) then
- msgstring = 'pop_in:horiz_grid_file '//trim(horiz_grid_file)//' not found'
+ string1 = 'pop_in:horiz_grid_file '//trim(horiz_grid_file)//' not found'
call error_handler(E_ERR,'read_horiz_grid', &
- msgstring, source, revision, revdate)
+ string1, source, revision, revdate)
endif
! Open it and read them in the EXPECTED order.
@@ -487,10 +515,11 @@
end subroutine read_horiz_grid
- subroutine calc_tpoints(nx, ny, ULAT, ULON, TLAT, TLON)
!------------------------------------------------------------------
-! subroutine calc_tpoints(nx, ny, ULAT, ULON, TLAT, TLON)
-!
+
+
+subroutine calc_tpoints(nx, ny, ULAT, ULON, TLAT, TLON)
+
! mimic POP grid.F90:calc_tpoints(), but for one big block.
integer, intent( in) :: nx, ny
@@ -516,7 +545,7 @@
radian = 180.0_r8/pi
! initialize these arrays to 0. in the code below there is
-! a column that is referenced by a where() construct before
+! a column that is referenced by a where() construct before
! those values are set. make sure that it doesn't cause a
! floating point exception from random memory bits which aren't
! valid floating point numbers.
@@ -526,7 +555,7 @@
do j=2,ny
do i=2,nx
- !*** convert neighbor U-cell coordinates to 3-d Cartesian coordinates
+ !*** 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))
@@ -583,7 +612,7 @@
where (TLON(:,:) > pi2) TLON(:,:) = TLON(:,:) - pi2
where (TLON(:,:) < c0 ) TLON(:,:) = TLON(:,:) + pi2
-!*** this leaves the leftmost/western edge to be filled
+!*** 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.
@@ -595,19 +624,19 @@
TLON(1,:) = (TLON(2,:) + TLON(nx,:))/c2
else
- write(msgstring,'(''pop_in&domain_nml:ew_boundary_type '',a,'' unknown.'')') &
+ write(string1,'(''pop_in&domain_nml:ew_boundary_type '',a,'' unknown.'')') &
trim(ew_boundary_type)
- call error_handler(E_ERR,'calc_tpoints',msgstring,source,revision,revdate)
+ call error_handler(E_ERR,'calc_tpoints',string1,source,revision,revdate)
endif
end subroutine calc_tpoints
+!------------------------------------------------------------------
- subroutine read_topography(nx, ny, KMT, KMU)
-!------------------------------------------------------------------
-! 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
@@ -620,9 +649,9 @@
! Check to see that the file exists.
if ( .not. file_exist(topography_file) ) then
- msgstring = 'pop_in:topography_file '//trim(topography_file)//' not found'
+ string1 = 'pop_in:topography_file '//trim(topography_file)//' not found'
call error_handler(E_ERR,'read_topography', &
- msgstring, source, revision, revdate)
+ string1, source, revision, revdate)
endif
! read the binary file
@@ -643,17 +672,17 @@
! south and west. so the T points which surround any U(i,j) point are
! in fact at indices i,i+1, and j,j+1 .
!
-! NO: KMU(i,j) = min(KMT(i, j), KMT(i-1, j), KMT(i, j-1), KMT(i-1, j-1))
+! NO: KMU(i,j) = min(KMT(i, j), KMT(i-1, j), KMT(i, j-1), KMT(i-1, j-1))
! YES: KMU(i,j) = min(KMT(i, j), KMT(i+1, j), KMT(i, j+1), KMT(i+1, j+1))
!
! the latter matches the POP source code, on yellowstone, lines 908 and 909 in:
! /glade/p/cesm/releases/cesm1_2_2/models/ocn/pop2/source/grid.F90
!
! wrap around longitude boundary at i == nx. set the topmost (last) latitude
-! U row to the same value in all cases. in the shifted pole grid currently in
+! U row to the same value in all cases. in the shifted pole grid currently in
! use all these points are on land and so are 0. in the original unshifted
! lat/lon grid these last row U points are above the final T row and are believed
-! to be unused. for completeness we set all values in the last U row to the
+! to be unused. for completeness we set all values in the last U row to the
! minimum of the all T row values immediately below it, for all longitudes.
do j=1,ny-1
@@ -667,11 +696,11 @@
end subroutine read_topography
+!------------------------------------------------------------------
- subroutine read_vert_grid(nz, ZC, ZG)
-!------------------------------------------------------------------
-! 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 in ascii, with either 3 columns/line
@@ -688,7 +717,7 @@
real(r8) :: depth
logical :: three_columns
-character(len=256) :: line
+character(len=256) :: line
real(r8), parameter :: centemeters_to_meters = 0.01_r8
@@ -697,9 +726,9 @@
! Check to see that the file exists.
if ( .not. file_exist(vert_grid_file) ) then
- msgstring = 'pop_in:vert_grid_file '//trim(vert_grid_file)//' not found'
+ string1 = 'pop_in:vert_grid_file '//trim(vert_grid_file)//' not found'
call error_handler(E_ERR,'read_vert_grid', &
- msgstring, source, revision, revdate)
+ string1, source, revision, revdate)
endif
! read the ASCII file
@@ -715,28 +744,28 @@
three_columns = .true.
else
three_columns = .false.
-
+
! read depth and calculate center and bottom of cells
- read(line,*,iostat=ios) depth
+ read(line,*,iostat=ios) depth
- ZC(1) = depth*centemeters_to_meters*0.5
+ ZC(1) = depth*centemeters_to_meters*0.5_r8
ZG(1) = depth*centemeters_to_meters
endif
do i=2, nz
-
+
if(three_columns) then
read(iunit,*,iostat=ios) depth, ZC(i), ZG(i)
else
read(iunit,*,iostat=ios) depth
- ZC(i) = ZG(i-1) + depth*centemeters_to_meters*0.5
+ ZC(i) = ZG(i-1) + depth*centemeters_to_meters*0.5_r8
ZG(i) = ZG(i-1) + depth*centemeters_to_meters
endif
if ( ios /= 0 ) then ! error
- write(msgstring,*)'error reading depths, line ',i
- call error_handler(E_ERR,'read_vert_grid',msgstring,source,revision,revdate)
+ write(string1,*)'error reading depths, line ',i
+ call error_handler(E_ERR,'read_vert_grid',string1,source,revision,revdate)
endif
enddo
@@ -748,6 +777,7 @@
subroutine get_pop_restart_filename( filename )
+
character(len=*), intent(OUT) :: filename
if ( .not. module_initialized ) call initialize_module
Modified: DART/trunk/models/POP/model_mod.f90
===================================================================
--- DART/trunk/models/POP/model_mod.f90 2016-11-10 18:10:41 UTC (rev 10737)
+++ DART/trunk/models/POP/model_mod.f90 2016-11-10 21:30:31 UTC (rev 10738)
@@ -80,8 +80,8 @@
! things which can/should be in the model_nml
logical :: output_state_vector = .true.
-integer :: assimilation_period_days = 1
-integer :: assimilation_period_seconds = 0
+integer :: assimilation_period_days = -1
+integer :: assimilation_period_seconds = -1
real(r8) :: model_perturbation_amplitude = 0.2
logical :: update_dry_cell_walls = .false.
integer :: debug = 0 ! turn up for more and more debug messages
@@ -271,18 +271,16 @@
if (do_output()) write(logfileunit, nml=model_nml)
if (do_output()) write( * , nml=model_nml)
-
! Set the time step ... causes POP namelists to be read.
! Ensures model_timestep is multiple of 'ocean_dynamics_timestep'
-model_timestep = set_model_time_step()
+model_timestep = set_model_time_step(assimilation_period_seconds, assimilation_period_days)
call get_time(model_timestep,ss,dd) ! set_time() assures the seconds [0,86400)
write(string1,*)'assimilation period is ',dd,' days ',ss,' seconds'
call error_handler(E_MSG,'static_init_model',string1,source,revision,revdate)
-
! get data dimensions, then allocate space, then open the files
! and actually fill in the arrays.
@@ -788,10 +786,9 @@
! Local storage
real(r8) :: loc_array(3), llon, llat, lheight
-integer :: base_offset, offset, ind
+integer :: base_offset, ind
integer :: hgt_bot, hgt_top
real(r8) :: hgt_fract
-real(r8) :: top_val, bot_val
integer :: hstatus
logical :: convert_to_ssh
@@ -3420,7 +3417,6 @@
integer :: hstatus, hgt_bot, hgt_top
real(r8) :: hgt_fract, salinity_val, potential_temp, pres_val
-real(r8) :: pres_bot, pres_top
interp_val = MISSING_R8
istatus = 99
Modified: DART/trunk/models/POP/model_mod.html
===================================================================
--- DART/trunk/models/POP/model_mod.html 2016-11-10 18:10:41 UTC (rev 10737)
+++ DART/trunk/models/POP/model_mod.html 2016-11-10 21:30:31 UTC (rev 10738)
@@ -435,14 +435,24 @@
example.</TD></TR>
<TR><!--contents--><TD valign=top>assimilation_period_days</TD>
- <!-- type --><TD valign=top>integer <em class=units>[default: 1]</em></TD>
+ <!-- type --><TD valign=top>integer <em class=units>[default: -1]</em></TD>
<!--descript--><TD valign=top>The number of days to advance the model for each assimilation.
+ If both <em class=code>assimilation_period_days</em> and
+ <em class=code>assimilation_period_seconds</em> are ≤ 0; the value of
+ the POP namelist variables <em class=code>restart_freq</em> and
+ <em class=code>restart_freq_opt</em> are used to determine the
+ assimilation period. WARNING: in the CESM framework,
+ the <em class=code>restart_freq</em> is set to a value that is not
+ useful so DART defaults to 1 day - even if you are using POP in
+ the LANL framework.
</TD></TR>
<TR><!--contents--><TD valign=top>assimilation_period_seconds</TD>
- <!-- type --><TD valign=top>integer <em class=units>[default: 0]</em></TD>
+ <!-- type --><TD valign=top>integer <em class=units>[default: -1]</em></TD>
<!--descript--><TD valign=top>In addition to <em class=code>assimilation_period_days</em>,
- the number of seconds to advance the model for each assimilation.
+ the number of seconds to advance the model for each assimilation.
+ Make sure you read the description of
+ <em class=code>assimilation_period_days</em>
</TD></TR>
<TR><!--contents--><TD valign=top>model_perturbation_amplitude</TD>
Modified: DART/trunk/models/POP/model_mod.nml
===================================================================
--- DART/trunk/models/POP/model_mod.nml 2016-11-10 18:10:41 UTC (rev 10737)
+++ DART/trunk/models/POP/model_mod.nml 2016-11-10 21:30:31 UTC (rev 10738)
@@ -1,6 +1,6 @@
&model_nml
- assimilation_period_days = 7,
- assimilation_period_seconds = 0,
+ assimilation_period_days = -1,
+ assimilation_period_seconds = -1,
model_perturbation_amplitude = 0.2,
output_state_vector = .true.,
debug = 0 /
Modified: DART/trunk/models/POP/work/input.nml
===================================================================
--- DART/trunk/models/POP/work/input.nml 2016-11-10 18:10:41 UTC (rev 10737)
+++ DART/trunk/models/POP/work/input.nml 2016-11-10 21:30:31 UTC (rev 10738)
@@ -146,8 +146,8 @@
/
&model_nml
- assimilation_period_days = 1,
- assimilation_period_seconds = 0,
+ assimilation_period_days = -1,
+ assimilation_period_seconds = -1,
model_perturbation_amplitude = 0.2,
output_state_vector = .false.,
debug = 0,
More information about the Dart-dev
mailing list