<p><b>duda</b> 2011-09-20 18:16:08 -0600 (Tue, 20 Sep 2011)</p><p>BRANCH COMMIT<br>
<br>
- Add code to read in correct time slice from surface update file when restarting.<br>
- Write xtime into surface update and initial condition files.<br>
- Restructure the main time integration loop in mpas_core_run():<br>
<br>
1) Write initial time (t=0) output <br>
2) Loop while not clock stop time<br>
2.1) Read in surface update fields if alarm is ringing<br>
2.2) Call timestep() : physics, RK time step (incl. microphysics)<br>
2.3) Shift time level 2 -> 1<br>
2.4) Advance clock and increment timestep<br>
2.5) Write output fields if alarm is ringing<br>
2.6) Write restart fields if alarm is ringing<br>
<br>
<br>
M src/core_init_nhyd_atmos/module_mpas_core.F<br>
M src/core_init_nhyd_atmos/module_test_cases.F<br>
M src/core_init_nhyd_atmos/Registry<br>
M src/core_nhyd_atmos/module_mpas_core.F<br>
M src/core_nhyd_atmos/module_time_integration.F<br>
</p><hr noshade><pre><font color="gray">Modified: branches/atmos_physics/src/core_init_nhyd_atmos/Registry
===================================================================
--- branches/atmos_physics/src/core_init_nhyd_atmos/Registry 2011-09-20 20:17:42 UTC (rev 1013)
+++ branches/atmos_physics/src/core_init_nhyd_atmos/Registry 2011-09-21 00:16:08 UTC (rev 1014)
@@ -52,7 +52,7 @@
#
# var type name_in_file ( dims ) iro- name_in_code super-array array_class
#
-var persistent real xtime ( Time ) 2 o xtime state - -
+var persistent text xtime ( Time ) 2 so xtime state - -
# horizontal grid structure
Modified: branches/atmos_physics/src/core_init_nhyd_atmos/module_mpas_core.F
===================================================================
--- branches/atmos_physics/src/core_init_nhyd_atmos/module_mpas_core.F 2011-09-20 20:17:42 UTC (rev 1013)
+++ branches/atmos_physics/src/core_init_nhyd_atmos/module_mpas_core.F 2011-09-21 00:16:08 UTC (rev 1014)
@@ -15,8 +15,16 @@
type (domain_type), intent(inout) :: domain
character(len=*), intent(out) :: startTimeStamp
+ type (block_type), pointer :: block
+
startTimeStamp = config_start_time
-
+
+ block => domain % blocklist
+ do while (associated(block))
+ block % state % time_levs(1) % state % xtime % scalar = startTimeStamp
+ block => block % next
+ end do
+
end subroutine mpas_core_init
@@ -33,11 +41,7 @@
type (io_output_object), intent(inout) :: output_obj
integer, intent(inout) :: output_frame
- integer :: ntimesteps, itimestep
- real (kind=RKIND) :: dt
- type (block_type), pointer :: block
-
-
+
call setup_nhyd_test_case(domain)
!
@@ -46,13 +50,7 @@
!
! call initialize_advection_rk(mesh)
! call initialize_deformation_weights(mesh)
-
- block => domain % blocklist
- do while (associated(block))
- block % state % time_levs(1) % state % xtime % scalar = 0.0
- block => block % next
- end do
-
+
call output_state_for_domain(output_obj, domain, output_frame)
end subroutine mpas_core_run
Modified: branches/atmos_physics/src/core_init_nhyd_atmos/module_test_cases.F
===================================================================
--- branches/atmos_physics/src/core_init_nhyd_atmos/module_test_cases.F 2011-09-20 20:17:42 UTC (rev 1013)
+++ branches/atmos_physics/src/core_init_nhyd_atmos/module_test_cases.F 2011-09-21 00:16:08 UTC (rev 1014)
@@ -4131,6 +4131,10 @@
call MPAS_advanceClock(fg_clock)
curr_time = MPAS_getClockTime(fg_clock, MPAS_NOW)
+
+ call MPAS_getTime(curr_time, dateTimeString=timeString)
+ state % xtime % scalar = timeString
+
end do
call output_state_finalize(sfc_update_obj, dminfo)
Modified: branches/atmos_physics/src/core_nhyd_atmos/module_mpas_core.F
===================================================================
--- branches/atmos_physics/src/core_nhyd_atmos/module_mpas_core.F 2011-09-20 20:17:42 UTC (rev 1013)
+++ branches/atmos_physics/src/core_nhyd_atmos/module_mpas_core.F 2011-09-21 00:16:08 UTC (rev 1014)
@@ -30,6 +30,11 @@
real (kind=RKIND) :: dt
type (block_type), pointer :: block
+ type (field1DChar) :: xtime
+ type (MPAS_Time_type) :: startTime, sliceTime
+ type (MPAS_TimeInterval_type) :: timeDiff, minTimeDiff
+ character(len=32) :: timeStamp
+ integer :: i
if (.not. config_do_restart) call setup_nhyd_test_case(domain)
@@ -43,20 +48,64 @@
block => domain % blocklist
do while (associated(block))
call mpas_init_block(domain % dminfo, block, block % mesh, dt)
- block % state % time_levs(1) % state % xtime % scalar = startTimeStamp
+ block % state % time_levs(1) % state % xtime % scalar = startTimeStamp
block => block % next
end do
restart_frame = 1
current_outfile_frames = 0
-!MGD TODO: this should pick up at the right frame when doing a restart
- sfc_update_obj % time = 1
- sfc_update_obj % filename = trim(config_sfc_update_name)
- sfc_update_obj % stream = STREAM_SFC
+ if (config_sfc_update_interval /= "none") then
- if (config_sfc_update_interval /= "none") call io_input_init(sfc_update_obj, domain % dminfo)
+ sfc_update_obj % filename = trim(config_sfc_update_name)
+ sfc_update_obj % stream = STREAM_SFC
+ call io_input_init(sfc_update_obj, domain % dminfo)
+
+ !
+ ! We need to decide which time slice to read from the surface file
+ !
+ sfc_update_obj % time = 1
+
+ if (sfc_update_obj % rdLocalTime <= 0) then
+ write(0,*) 'Error: Couldn''t find any times in surface update file.'
+ call dmpar_abort(domain % dminfo)
+ end if
+
+ if (domain % dminfo % my_proc_id == IO_NODE) then
+ allocate(xtime % ioinfo)
+ xtime % ioinfo % start(1) = 1
+ xtime % ioinfo % count(1) = sfc_update_obj % rdLocalTime
+ allocate(xtime % array(sfc_update_obj % rdLocalTime))
+
+ xtime % ioinfo % fieldName = 'xtime'
+ call io_input_field(sfc_update_obj, xtime)
+
+ call MPAS_setTimeInterval(interval=minTimeDiff, DD=10000)
+ call MPAS_setTime(curr_time=startTime, dateTimeString=config_start_time)
+
+ do i=1,sfc_update_obj % rdLocalTime
+ call MPAS_setTime(curr_time=sliceTime, dateTimeString=xtime % array(i))
+ timeDiff = abs(sliceTime - startTime)
+ if (timeDiff < minTimeDiff) then
+ minTimeDiff = timeDiff
+ sfc_update_obj % time = i
+ end if
+ end do
+
+ timeStamp = xtime % array(sfc_update_obj % time)
+
+ deallocate(xtime % ioinfo)
+ deallocate(xtime % array)
+ end if
+
+ call dmpar_bcast_int(domain % dminfo, sfc_update_obj % time)
+ call dmpar_bcast_char(domain % dminfo, timeStamp)
+
+ write(0,*) 'Starting model with surface time ', timeStamp
+
+ end if
+
end subroutine mpas_core_init
@@ -108,7 +157,7 @@
! set sfc alarm, if necessary
if (trim(config_sfc_update_interval) /= "none") then
call MPAS_setTimeInterval(alarmTimeStep, timeString=config_sfc_update_interval, ierr=ierr)
- alarmStartTime = startTime + alarmTimeStep
+ alarmStartTime = startTime
call MPAS_addClockAlarm(clock, sfcAlarmID, alarmStartTime, alarmTimeStep, ierr=ierr)
end if
@@ -229,18 +278,6 @@
! Eventually, dt should be domain specific
dt = config_dt
- currTime = MPAS_getClockTime(clock, MPAS_NOW, ierr)
- call MPAS_getTime(curr_time=currTime, dateTimeString=timeStamp, ierr=ierr)
- write(0,*) 'Initial timestep ', timeStamp
-
- if (config_sfc_update_interval /= "none") then
- call read_and_distribute_fields(domain % dminfo, sfc_update_obj, domain % blocklist, &
- readCellStart, nReadCells, readEdgeStart, nReadEdges, readVertexStart, nReadVertices, &
- readVertLevelStart, nReadVertLevels, &
- sendCellList, recvCellList, sendEdgeList, recvEdgeList, sendVertexList, recvVertexList, &
- sendVertLevelList, recvVertLevelList)
- end if
-
call write_output_frame(output_obj, output_frame, domain)
! During integration, time level 1 stores the model state at the beginning of the
@@ -248,22 +285,37 @@
itimestep = 0
do while (.not. MPAS_isClockStopTime(clock))
- itimestep = itimestep + 1
- call MPAS_advanceClock(clock)
-
currTime = MPAS_getClockTime(clock, MPAS_NOW, ierr)
call MPAS_getTime(curr_time=currTime, dateTimeString=timeStamp, ierr=ierr)
- write(0,*) 'Doing timestep ', timeStamp
+ write(0,*) 'Begin timestep ', timeStamp
+ ! Input external updates (i.e. surface)
+ if (MPAS_isAlarmRinging(clock, sfcAlarmID, ierr=ierr)) then
+ call MPAS_resetClockAlarm(clock, sfcAlarmID, ierr=ierr)
+
+ call read_and_distribute_fields(domain % dminfo, sfc_update_obj, domain % blocklist, &
+ readCellStart, nReadCells, readEdgeStart, nReadEdges, readVertexStart, nReadVertices, &
+ readVertLevelStart, nReadVertLevels, &
+ sendCellList, recvCellList, sendEdgeList, recvEdgeList, sendVertexList, recvVertexList, &
+ sendVertLevelList, recvVertLevelList)
+ sfc_update_obj % time = sfc_update_obj % time + 1
+ end if
+
call timer_start("time integration")
call mpas_timestep(domain, dt, itimestep)
call timer_stop("time integration")
! Move time level 2 fields back into time level 1 for next time step
call shift_time_levels_state(domain % blocklist % state)
-
+
+
+ ! Advance clock before writing output
+ itimestep = itimestep + 1
+ call MPAS_advanceClock(clock)
+ currTime = MPAS_getClockTime(clock, MPAS_NOW, ierr)
+ call MPAS_getTime(curr_time=currTime, dateTimeString=timeStamp, ierr=ierr)
+
!TODO: MPAS_getClockRingingAlarms is probably faster than multiple MPAS_isAlarmRinging...
-
if (MPAS_isAlarmRinging(clock, outputAlarmID, ierr=ierr)) then
call MPAS_resetClockAlarm(clock, outputAlarmID, ierr=ierr)
!output_frame will always be > 1 here unless it is reset after the output file is finalized
@@ -278,17 +330,6 @@
restart_frame = restart_frame + 1
end if
- if (MPAS_isAlarmRinging(clock, sfcAlarmID, ierr=ierr)) then
- call MPAS_resetClockAlarm(clock, sfcAlarmID, ierr=ierr)
-
- sfc_update_obj % time = sfc_update_obj % time + 1
- call read_and_distribute_fields(domain % dminfo, sfc_update_obj, domain % blocklist, &
- readCellStart, nReadCells, readEdgeStart, nReadEdges, readVertexStart, nReadVertices, &
- readVertLevelStart, nReadVertLevels, &
- sendCellList, recvCellList, sendEdgeList, recvEdgeList, sendVertexList, recvVertexList, &
- sendVertLevelList, recvVertLevelList)
- end if
-
end do
end subroutine mpas_core_run
Modified: branches/atmos_physics/src/core_nhyd_atmos/module_time_integration.F
===================================================================
--- branches/atmos_physics/src/core_nhyd_atmos/module_time_integration.F 2011-09-20 20:17:42 UTC (rev 1013)
+++ branches/atmos_physics/src/core_nhyd_atmos/module_time_integration.F 2011-09-21 00:16:08 UTC (rev 1014)
@@ -5,6 +5,7 @@
use constants
use dmpar
use vector_reconstruction
+ use mpas_timekeeping
#ifdef DO_PHYSICS
use module_driver_microphysics
@@ -33,6 +34,9 @@
integer, intent(in) :: itimestep
type (block_type), pointer :: block
+ type (MPAS_Time_type) :: currTime
+ type (MPAS_TimeInterval_type) :: dtInterval
+ character (len=32) :: xtime
if (trim(config_time_integration) == 'SRK3') then
call srk3(domain, dt, itimestep)
@@ -42,9 +46,14 @@
stop
end if
+ call MPAS_setTime(currTime, dateTimeString=timeStamp)
+ call MPAS_setTimeInterval(dtInterval, dt=dt)
+ currTime = currTime + dtInterval
+ call MPAS_getTime(currTime, dateTimeString=xtime)
+
block => domain % blocklist
do while (associated(block))
- block % state % time_levs(2) % state % xtime % scalar = timeStamp
+ block % state % time_levs(2) % state % xtime % scalar = xtime
block => block % next
end do
</font>
</pre>