Reading input data from netCDF files in the MPAS atmosphere model

[lmarker]

Hello!

I am trying to add a new input stream to the atmosphere_model that reads data from an .nc file. I have edited a standard MPAS static.nc file (it is the x1.10242.static.nc file generated in section 1.3 of the MPAS tutorial -> MPAS Tutorial — Practice Session Guide) that I am using to test with, called test_input.nc. I removed all vars aside from xtime and have added a variable 'test_var'. The dimensions and global attributes are identical to the original static file. The variables are now

variables:
        int test_var(TWO) ;
                test_var:units = "-" ;
                test_var:long_name = "Test variable" ;
        char xtime(Time, StrLen) ;
                xtime:units = "YYYY-MM-DD_hh:mm:ss" ;
                xtime:long_name = "Model valid time" ;

and the data is

data:
 xtime =
  "2014-09-10_00:00:00                                             " ;

 test_var =
  7777777, 666666                                              ;

I chose to use dimension TWO as it seemed like the simplest choice and I didn't want a dimensionless variable.

To the atmosphere model registry (../MAS-Model/src/core_atmosphere/Registry.xml) I added the following:

A new stream named test_input:

                <stream name="test_input" 
                        type="input" 
                        filename_template="test_input.nc" 
                        input_interval="initial_only"
                        immutable="true">

                        <var name="test_var" />
                </stream>

A new var_struct test_input_struct which has a variable define: test_var

        <var_struct name="test_input_struct" time_levs="1">

                <var name="test_var" type="integer" dimensions="TWO" units=""
                     description="Testing reading variables from a .nc file"/>

I wanted to try and use the test_var variable in some of the code, so in mpas_atm_core.F (../MPAS-Model/src/core_atmosphere/mpas_atm_core.F) I added the following to the atm_mpas_init_block subroutine.... (I have included some of the surrounding lines of code to help indicate where I have added my own).

Declare test_input_struct

      type (mpas_pool_type), pointer :: state
      type (mpas_pool_type), pointer :: diag
      type (mpas_pool_type), pointer :: tend
      type (mpas_pool_type), pointer :: sfc_input
      type (mpas_pool_type), pointer :: diag_physics
      type (mpas_pool_type), pointer :: atm_input
!-----------------Declaring-test_input-struct-------------
      type (mpas_pool_type), pointer :: test_input_struct

Declare variables I will be using:
- test_var: array defined in test_input struct
- TWO: the dimension of test_var
- sum_test_var: a variable I will use to test some arithmetic with the elements of test_var
- i: used for looping

      integer, pointer :: nThreads
      integer, dimension(:), pointer :: cellThreadStart, cellThreadEnd
      integer, dimension(:), pointer :: cellSolveThreadStart, cellSolveThreadEnd
      integer, dimension(:), pointer :: edgeThreadStart, edgeThreadEnd
      integer, dimension(:), pointer :: edgeSolveThreadStart, edgeSolveThreadEnd
      integer, dimension(:), pointer :: vertexThreadStart, vertexThreadEnd
      integer, dimension(:), pointer :: vertexSolveThreadStart, vertexSolveThreadEnd

      logical, pointer :: config_do_restart, config_do_DAcycling

      !-----------Declaring-vars-for-reading-.nc-files------

      integer, dimension(:), pointer :: test_var
      integer, pointer :: TWO
      integer :: sum_test_var, i

Call the test_input struct

      call mpas_pool_get_subpool(block % structs, 'diag', diag)
      call mpas_pool_get_subpool(block % structs, 'state', state)
!--------------Call-test_input-struct-subpool------------------------------------------
      call mpas_pool_get_subpool(block % structs, 'test_input_struct', test_input_struct)
      call mpas_pool_get_config(block % configs, 'config_do_restart', config_do_restart)
      call mpas_pool_get_config(block % configs, 'config_do_DAcycling', config_do_DAcycling)

Then, at the end of the subroutine, I have called test_var and TWO from their respective pools and tried to carry out some arithmetic with the elements of test_var in a do loop (sum the elements). I have also added some mpas_log_write_statements to see what's going on.

!-----------Using-test_var-------

      call mpas_pool_get_array(test_input_struct, 'test_var', test_var)
      call mpas_pool_get_dimension(mesh, 'TWO', TWO)

      sum_test_var = 0

      call mpas_log_write('sum_test_var initially set to $i.', intArgs=[sum_test_var])

      do i=1,TWO
         call mpas_log_write('Do loop, step $i', intArgs=[i])
         sum_test_var = sum_test_var + test_var(i)
         call mpas_log_write('test_var value is $i.',intArgs=[test_var(i)])
         call mpas_log_write('sum_test_var value is $i.',intArgs=[test_var(i)])
      end do

      call mpas_log_write('sum_test_var finally set to $i.', intArgs=[sum_test_var])


   end subroutine atm_mpas_init_block

CORE=atmosphere compiles after making the above edits to the registry and mpas_atm_core.F.

I added the following input stream to streams.atmosphere:

<immutable_stream name="test_input"
                  type="test_input"
                  filename_template="test_input.nc"
                  input_interval="initial_only" />

After running atmosphere_model, the following is output to log.atmosphere.0000.out

 sum_test_var initially set to 0.
 Do loop, step 1
 test_var value is 0.
 sum_test_var value is 0.
 Do loop, step 2
 test_var value is 0.
 sum_test_var value is 0.
 sum_test_var finally set to 0.

So, it seems that MPAS is happy with the test_var I defined in test_input_struct, as I can call test_var and use its values. BUT they have not picked up the value I gave them in test_input.nc, i.e. the input stream has not assigned values to test_var, and instead are set to the default of zero!

What I would like to know is, how do I get the test_var variable defined in the atmosphere model to take the values assigned in the test_input stream?

One thing that very possibly is an issue is the test_input.nc file that I am using. I am a novice with the netCDF format and my generation of the stream input file was therefore not very sophisticated!

Some other questions that have arisen during my attempt at adding a new input streams:

- Can you define new dimensions and if so is it done as with variables in the registry?
- What does 'time_levs' in the registry represent?

Thanks in advance,

Laurents

 

[mgduda]

Apologies for the long delay in posting a reply! It looks like you've done everything correctly, and the issue may simply be that the "test_input" stream has not yet been read at the point where you've added code to compute sum_test_var.

I'll admit that the code in mpas_atm_core.F isn't as readable as it could be. I've created a pared-down version of the code, below, that should make it easier to see the essential points at which input streams are read.

module atm_core

   function atm_core_init(domain, startTimeStamp) result(ierr)

      !
      ! If this is a restart run, read the restart stream, else read the input
      ! stream.
      ! Regardless of which stream we read for initial conditions, reset the
      ! input alarms for both input and restart before reading any remaining
      ! input streams.
      !
      if (config_do_restart) then
         call MPAS_stream_mgr_read(domain % streamManager, streamID='restart', ierr=ierr)
      else
         call MPAS_stream_mgr_read(domain % streamManager, streamID='input', ierr=ierr)
      end if
      call MPAS_stream_mgr_reset_alarms(domain % streamManager, streamID='input', direction=MPAS_STREAM_INPUT, ierr=ierr)
      call MPAS_stream_mgr_reset_alarms(domain % streamManager, streamID='restart', direction=MPAS_STREAM_INPUT, ierr=ierr)

      ! Note: Because the 'iau' stream has the 'iau' package attached to it, the MPAS_stream_mgr_read( )
      !       call here will actually try to read the stream only if IAU is being used in the run.
      !
      call MPAS_stream_mgr_read(domain % streamManager, streamID='iau', whence=MPAS_STREAM_NEAREST, ierr=ierr)
      call MPAS_stream_mgr_reset_alarms(domain % streamManager, streamID='iau', ierr=ierr)

      block => domain % blocklist
      do while (associated(block))

         call atm_mpas_init_block(domain % dminfo, domain % streamManager, block, mesh, dt)

         block => block % next
      end do

   end function atm_core_init
  

   function atm_core_run(domain) result(ierr)

      do while (.not. mpas_is_clock_stop_time(clock))

         !
         ! Read external field updates
         !
         call MPAS_stream_mgr_begin_iteration(domain % streamManager, ierr=ierr)
         do while (MPAS_stream_mgr_get_next_stream(domain % streamManager, streamID=input_stream, directionProperty=stream_dir))
            if (stream_dir == MPAS_STREAM_INPUT .or. stream_dir == MPAS_STREAM_INPUT_OUTPUT) then
               if (MPAS_stream_mgr_ringing_alarms(domain % streamManager, streamID=input_stream, &
                                                  direction=MPAS_STREAM_INPUT, ierr=ierr)) then
                  call mpas_dmpar_get_time(input_start_time)
                  call MPAS_stream_mgr_read(domain % streamManager, streamID=input_stream, whence=MPAS_STREAM_LATEST_BEFORE, &
                                            actualWhen=read_time, ierr=ierr)
                  call mpas_dmpar_get_time(input_stop_time)
                  if (ierr /= MPAS_STREAM_MGR_NOERR) then
                     call mpas_log_write('********************************************************************************', messageType=MPAS_LOG_ERR)
                     call mpas_log_write('Error reading input stream '//trim(input_stream),                                  messageType=MPAS_LOG_ERR)
                     call mpas_log_write('********************************************************************************', messageType=MPAS_LOG_CRIT)
                  end if

                  call mpas_log_write('----------------------------------------------------------------------')
                  call mpas_log_write('  Read '''//trim(input_stream)//''' input stream valid at '//trim(read_time))
                  call mpas_log_write('   Timing for stream input: $r s', realArgs=(/real(input_stop_time - input_start_time, kind=RKIND)/))
                  call mpas_log_write('----------------------------------------------------------------------')

                  call MPAS_stream_mgr_reset_alarms(domain % streamManager, streamID=input_stream, direction=MPAS_STREAM_INPUT, ierr=ierr)
               end if
            end if
         end do

         call mpas_timer_start("time integration")
         call mpas_dmpar_get_time(integ_start_time)
         call atm_do_timestep(domain, dt, itimestep)
         call mpas_dmpar_get_time(integ_stop_time)
         call mpas_timer_stop("time integration")
         call mpas_log_write(' Timing for integration step: $r s', realArgs=(/real(integ_stop_time - integ_start_time, kind=RKIND)/))

      end do

   end function atm_core_run

end module atm_core

Essentially, we only read either the "restart" or "input" stream (plus the IAU stream if doing incremental analysis updates) in atm_core_init before the call to atm_mpas_init_block. The remainder of the input streams are read near the beginning of each time step in atm_core_run, depending on whether those input streams have an input alarm that is ringing.

So for your case, there are a couple of options:

1. Move your code after the block of code to read input streams in atm_core_run after the "Read external field updates" comment;
2. Add a few lines of code near the place where the "input" or "restart" (or "iau") streams are read so that your "test_input" stream is read as well.

 

[mgduda]

Some other questions that have arisen during my attempt at adding a new input streams:

- Can you define new dimensions and if so is it done as with variables in the registry?
- What does 'time_levs' in the registry represent?

New dimensions will indeed be defined in the Registry.xml file. There are a few options for these:

• The value of the dimension will come from an input stream (and its value must be consistent across all input streams where it appears), e.g.,

  <dim name="nCells" description="The number of Voronoi cells in the primal mesh"/>

• The value of the dimension is given as an immediate value (and its value must be as specified in all input files where it appears), e.g.,

  <dim name="TWO" definition="2" description="A constant value of 2"/>

• The value of the dimension is related in a simple way (I'm not sure which arithmetic is supported at present...) to the value of another dimension, e.g.,

  <dim name="nVertLevelsP1" definition="nVertLevels+1" description="The number of atmospheric levels, always one more than the number of layers"/>

• The value of the dimension is set based on a namelist variable, e.g.,

  <dim name="nSoilLevels" definition="namelist:num_soil_layers" description="The number of soil layers used by the land-surface scheme"/>

The "time_levs" attribute specifies that fields have more than one time level (or, more generally speaking, more than one copy) allocated for them. These time levels don't appear as regular array dimensions, but the various time levels of a field can be accessed by specifying the optional "timeLevel" argument to the mpas_pool_get_array routine, e.g.,

call mpas_pool_get_array(state, 'theta_m', theta_m_t1, timeLevel=1)
call mpas_pool_get_array(state, 'theta_m', theta_m_t2, timeLevel=2)

The main use for multiple time levels so far has been to hold provisional or intermediate values of model state that are needed by integration schemes. But, for example, we've also found this useful in our implementation of lateral boundary conditions, where we store the boundary state at the last boundary update time in timeLevel=1 and the boundary state at the next boundary update time (we read from future times in the lbc_in stream) in timeLevel=2; we can then easily compute tendencies for the boundary state from these two time levels.