removed lbc domain

Author: hkershaw-brown

models/mpas_atm/model_mod.f90

@@ -165,13 +165,8 @@ module model_mod
           cell_next_to_boundary_edge
 
 public :: update_u_from_reconstruct, &
-          set_lbc_variables, &
-          force_u_into_state, &
           use_increments_for_u_update, &
-          anl_domid, lbc_domid ! HK todo accessor functions?
-
-! set_lbc_variables sets the lbc_variables string array
-! force_u_into_state sets a logical add_u_to_state_list that forces u to be in state
+          anl_domid ! HK todo accessor functions?
 
 character(len=*), parameter :: source   = 'models/mpas_atm/model_mod.f90'
 
@@ -231,7 +226,6 @@ module model_mod
 logical :: search_initialized = .false.
 
 integer :: anl_domid = -1 
-integer :: lbc_domid = -1
 
 !------------------------------------------------------------------------
 ! variables which are in the model_nml namelist
@@ -265,11 +259,6 @@ module model_mod
 ! assimilation, and use only the increments to update the edge winds
 logical :: use_increments_for_u_update = .true.
 
-! if set by: call force_u_into_state()  *BEFORE* calling static_init_model(),
-! the code will add U (edge normal winds) to the mpas state vector even if it
-! isn't in the namelist.
-logical :: add_u_to_state_list = .false.
-
 ! if > 0, amount of distance in fractional model levels that a vertical
 ! point can be above or below the top or bottom of the grid and still be
 ! evaluated without error. 
@@ -323,27 +312,9 @@ module model_mod
 !------------------------------------------------------------------------
 ! for regional MPAS
 real(r8) :: dxmax  ! max distance between two adjacent cell centers in the mesh (in meters)
-
-! when updating boundary files for regional mpas, note whether the boundary
-! file has the reconstructed winds (lbc_ur, lbc_vr) or not. (this is set by
-! looking at the bdy template file and seeing if those variables are there.)  
-! if not, the other two options are ignored.  
-! if they are in the lbc file, then the other logicals control whether to use 
-! them instead of updating the U edge winds directly, and whether to use the 
-! reconstructed increments or not. 
-! the latter two options could be added to the namelist if someone wanted to
-! explore the options for how the edge winds are updated in the boundary file.
-! for now they're not - they're hardcoded true.
-! note that these are for the boundary file update only - there are separate
-! options for how to update the U winds in the main assimilation state.
-
-logical :: lbc_file_has_reconstructed_winds     = .false.  
-
-! this is special and not in the namelist.  boundary files have a fixed
-! set of variables with fixed names.
-character(len=vtablenamelength) :: lbc_variables(MAX_STATE_VARIABLES) = ''
 !------------------------------------------------------------------------
 
+
 ! Grid parameters - the values will be read from an mpas analysis file.
 integer :: nCells        = -1  ! Total number of cells making up the grid
 integer :: nVertices     = -1  ! Unique points in grid that are corners of cells
@@ -432,7 +403,7 @@ subroutine static_init_model()
 integer :: iunit, io, ivar, i
 integer :: ss, dd, z1, m1
 integer :: nDimensions, nVariables, nAttributes, unlimitedDimID, TimeDimID
-integer :: nfields, lbc_nfields
+integer :: nfields
 logical :: both
 real(r8) :: variable_bounds(MAX_STATE_VARIABLES, 2)
 integer :: variable_qtys(MAX_STATE_VARIABLES)
@@ -479,25 +450,6 @@ subroutine static_init_model()
 
 model_size = get_domain_size(anl_domid)
 
-lbc_nfields = 0
-
-! if we have a lateral boundary file, add it to the domain
-! so we have access to the corresponding lbc_xxx fields.
-if (.not. global_grid .and. lbc_variables(1) /= '') then
-   ! regional: count number of lbc fields to read in
-   COUNTUP: do i=1, MAX_STATE_VARIABLES
-      if (lbc_variables(i) /= '') then
-         lbc_nfields = lbc_nfields + 1
-      else
-         exit COUNTUP
-      endif
-   enddo COUNTUP
-   lbc_domid = add_domain(bdy_template_filename, lbc_nfields,    &
-                          var_names  = lbc_variables)
-                        ! FIXME clamp_vals = variable_bounds(1:nfields,:) )
-   model_size = model_size + get_domain_size(lbc_domid)
-endif
-
 ! if you have at least one of these wind components in the state vector,
 ! you have to have them both.  the subroutine will error out if only one
 ! is found and not both.
@@ -1402,57 +1354,6 @@ end subroutine get_close_state
 !  (these are not required by dart but are used by other programs)
 !==================================================================
 
-
-!-------------------------------------------------------------------
-! modify what static_init_model does.  this *must* be called before !HK Nope this is not good.
-! calling static_init_model().
-! the boundary file variables are fixed by the model and so we
-! don't allow the user to set them via namelist
-
-subroutine set_lbc_variables(template_filename)
-
-character(len=*), intent(in) :: template_filename
-
-integer :: ncid
-
-bdy_template_filename = template_filename
-
-! this initial list always exists.  hardcode them for now,
-! and query to see if the reconstructed winds are there or not.
-
-lbc_variables(1) = 'lbc_qc'
-lbc_variables(2) = 'lbc_qr'
-lbc_variables(3) = 'lbc_qv'
-lbc_variables(4) = 'lbc_rho'
-lbc_variables(5) = 'lbc_theta'
-lbc_variables(6) = 'lbc_u'
-lbc_variables(7) = 'lbc_w'
-
-ncid = nc_open_file_readonly(template_filename, 'set_lbc_variables')
-if (nc_variable_exists(ncid, 'lbc_ur')) then
-   lbc_variables(8) = 'lbc_ur'  
-   lbc_variables(9) = 'lbc_vr' 
-   lbc_file_has_reconstructed_winds = .true.
-endif
-call nc_close_file(ncid)
-
-end subroutine set_lbc_variables
-
-
-!-------------------------------------------------------------------
-! modify what static_init_model does.  this *must* be called before
-! calling static_init_model().
-! set a logical add_u_to_state_list that forces u to be in state
-
-subroutine force_u_into_state()
-
-add_u_to_state_list = .true.
-
-end subroutine force_u_into_state
-
-
-!------------------------------------------------------------------
-
 function get_analysis_time_ncid( ncid, filename )
 
 ! The analysis netcdf files have the start time of the experiment.
@@ -1904,14 +1805,6 @@ subroutine verify_state_variables(ncid, filename, ngood, qty_list, variable_boun
    call error_handler(E_MSG,'verify_state_variables',string1,source,text2=string2)
 endif
 
-! if this flag is true and the user hasn't said U should be in the state,
-! add it to the list.
-if (add_u_to_state_list .and. .not. u_already_in_list) then
-   ngood = ngood + 1
-   variable_table(ngood,1) = "u"
-   variable_table(ngood,2) = "QTY_EDGE_NORMAL_SPEED"
-endif
-
 end subroutine verify_state_variables
 
 

models/mpas_atm/update_bc.f90

@@ -27,17 +27,16 @@ program update_bc
                              get_next_filename, E_ERR, E_MSG, error_handler
 use time_manager_mod, only : time_type, print_time, print_date, operator(-), &
                              get_time, get_date, operator(/=)
-use direct_netcdf_mod,only : read_variables  !HK read_variables?!
 use        model_mod, only : static_init_model, &
                              get_model_size,    &
                              get_analysis_time, &
-                             anl_domid, lbc_domid, &
                              uv_increments_cell_to_edges, &
                              uv_field_cell_to_edges, &
                              get_analysis_weight, &
                              on_boundary_cell, &
                              on_boundary_edge, &
-                             cell_next_to_boundary_edge
+                             cell_next_to_boundary_edge, &
+                             get_grid_dims
 
 
 use state_structure_mod, only : get_num_variables, get_domain_size, &
@@ -48,7 +47,8 @@ program update_bc
                                  nc_open_file_readwrite, &
                                  nc_get_dimension_size,   &
                                  nc_close_file, &
-                                 NF90_MAX_NAME
+                                 NF90_MAX_NAME, &
+                                 nc_get_variable, nc_put_variable
 
 implicit none
 
@@ -87,13 +87,16 @@ program update_bc
 integer :: vertexDegree  = -1  ! Max number of cells/edges that touch any vertex
 integer :: nSoilLevels   = -1  ! Number of soil layers
 
-integer :: cellid, edgeid, varid
+integer :: cellid, edgeid, i
 
+logical :: lbc_file_has_reconstructed_winds  ! HK @todo may always be false
+character(len=NF90_MAX_NAME), dimension(9) :: lbc_variables
 real(r8), allocatable, dimension(:,:) :: lbc_u, lbc_ucell, lbc_vcell, inc_lbc_ucell, inc_lbc_vcell
 real(r8), allocatable, dimension(:,:) :: old_lbc_ucell, old_lbc_vcell, delta_u
 real(r8), allocatable, dimension(:,:) :: a_var_data, b_var_data, var_data
 real(r8) :: weight
 
+
 character(len=NF90_MAX_NAME) :: bvarname, avarname
 
 !----------------------------------------------------------------------
@@ -107,6 +110,18 @@ program update_bc
 
 call static_init_model()
 call get_grid_dims(nCells, nVertices, nEdges, nVertLevels, vertexDegree, nSoilLevels)
+
+lbc_file_has_reconstructed_winds = .false.
+lbc_variables(1) = 'lbc_qc'
+lbc_variables(2) = 'lbc_qr'
+lbc_variables(3) = 'lbc_qv'
+lbc_variables(4) = 'lbc_rho'
+lbc_variables(5) = 'lbc_theta'
+lbc_variables(6) = 'lbc_u'
+lbc_variables(7) = 'lbc_w'
+lbc_variables(8) = 'lbc_ur'  ! may not be in the file
+lbc_variables(9) = 'lbc_vr'  ! may not be in the file
+
 !----------------------------------------------------------------------
 ! Reads lists of input mpas (prior) and filter (analysis) files 
 !HK @todo loop around files, why not run this code in parallel?
@@ -178,8 +193,8 @@ program update_bc
    allocate(b_var_data(nVertLevels, nCells))
    allocate(var_data(nVertLevels, nCells))
 
-   VARLOOP: do varid = 1, get_num_variables(lbc_domid)
-      bvarname = get_variable_name(lbc_domid, varid)
+   VARLOOP: do i = 1, size(lbc_variables)
+      bvarname = lbc_variables(i)
       avarname = trim(bvarname(5:)) !corresponding field in analysis domain
 
       ! skip edge normal winds
@@ -189,14 +204,6 @@ program update_bc
       if (bvarname == 'lbc_ur') avarname = 'uReconstructZonal'
       if (bvarname == 'lbc_vr') avarname = 'uReconstructMeridional'    
 
-      if (bvarname(1:4) /= 'lbc_') then
-         write(string1, *) 'skipping update of boundary variable ', trim(bvarname)  !HK @todo why do you need to tell people this?
-         write(string2, *) 'because the name does not start with "lbc"'
-         call error_handler(E_MSG,'statevector_to_boundary_file',string1,&
-                              source, text2=string2)
-         cycle VARLOOP
-      endif
-
       call nc_get_variable(ncAnlID, avarname, a_var_data)
       call nc_get_variable(ncBdyID, bvarname, b_var_data)
 
@@ -272,7 +279,7 @@ program update_bc
             ! has bdyMaskEdge = 0. In this case, even if bdyMaskEdge of the edge is zero,
             ! cell2 has been updated in the CELLS loop above, so the edge has to be updated.
 
-            if (.not. on_boundary_edge(edgeid) .and. .not. cell_next_to_boundary_edge(edgid)) cycle IEDGE
+            if (.not. on_boundary_edge(edgeid) .and. .not. cell_next_to_boundary_edge(edgeid)) cycle IEDGE
 
             lbc_u(:,edgeid) = lbc_u(:,edgeid) + delta_u(:,edgeid)
 

models/mpas_atm/update_mpas_states.f90

@@ -32,7 +32,7 @@ program update_mpas_states
                              get_model_size, &
                              get_analysis_time, update_u_from_reconstruct, &
                              use_increments_for_u_update, uv_increments_cell_to_edges, &
-                             uv_field_cell_to_edges, anl_domid
+                             uv_field_cell_to_edges, dom_id => anl_domid
 
 use state_structure_mod, only : get_num_variables, get_variable_name, &
                                 get_variable_size, get_varid_from_varname, &
@@ -67,7 +67,6 @@ program update_mpas_states
 integer               :: dims(2) !(nVertLevels, nEdges | nCells)
 real(r8), allocatable :: u(:,:), ucell(:,:), vcell(:,:) 
 real(r8), allocatable :: ucell_dart(:,:), vcell_dart(:,:), increments(:,:)
-integer :: dom_id = anl_domid
 !----------------------------------------------------------------------
 
 call initialize_utilities(progname=source)