feat: update wind components in update_mpas_states

This was a subroutine update_wind_components in model_mod that was called
by statevector_to_analysis_file

untested

Author: hkershaw-brown

models/mpas_atm/model_mod.f90

@@ -165,6 +165,7 @@ module model_mod
 public :: update_u_from_reconstruct, &
           set_lbc_variables, &
           force_u_into_state, &
+          use_increments_for_u_update, &
           statevector_to_boundary_file
 
 ! set_lbc_variables sets the lbc_variables string array
@@ -550,11 +551,11 @@ subroutine static_init_model()
    if (use_increments_for_u_update) then
       write(string1,*) 'edge normal winds (u) in MPAS file will be updated based on the difference'
       write(string2,*) 'between the original U/V reconstructed winds and the updated values'
-      write(string3,*) 'because use_increment_for_u_update is True'
+      write(string3,*) 'because use_increments_for_u_update is True'
    else
       write(string1,*) 'edge normal winds (u) in MPAS file will be updated by averaging the'
       write(string2,*) 'updated U/V reconstructed winds at the corresponding cell centers'
-      write(string3,*) 'because use_increment_for_u_update is False'
+      write(string3,*) 'because use_increments_for_u_update is False'
    endif
    call error_handler(E_MSG,'static_init_model',string1,source, &
                       text2=string2, text3=string3)

models/mpas_atm/update_mpas_states.f90

@@ -30,10 +30,12 @@ program update_mpas_states
                              get_time, get_date, operator(/=)
 use        model_mod, only : static_init_model, &
                              get_model_size, &
-                             get_analysis_time, update_u_from_reconstruct
+                             get_analysis_time, update_u_from_reconstruct, &
+                             use_increments_for_u_update, uv_cell_to_edges
 
 use state_structure_mod, only : get_num_variables, get_variable_name, &
-                                get_variable_size, get_varid_from_varname
+                                get_variable_size, get_varid_from_varname, &
+                                get_dim_lengths
 
 use netcdf_utilities_mod, only : nc_open_file_readonly, &
                                  nc_open_file_readwrite, &
@@ -62,6 +64,10 @@ program update_mpas_states
 real(r8), allocatable :: variable(:)
 type(time_type)       :: model_time
 type(time_type)       :: state_time
+integer               :: dims(2) !(nVertLevels, nEdges | nCells)
+real(r8), allocatable :: u(:,:), ucell(:,:), vcell(:,:) 
+real(r8), allocatable :: ucell_dart(:,:), vcell_dart(:,:), increments(:,:)
+integer :: dom_id = 1 ! HK @todo
 !----------------------------------------------------------------------
 
 call initialize_utilities(progname=source)
@@ -85,8 +91,8 @@ program update_mpas_states
   next_outfile = get_next_filename(update_output_file_list, filenum)
   if (next_infile == '' .or. next_outfile == '') exit fileloop
 
-  ncAnlID = nc_open_file_readonly(next_infile, 'update_mpas_states - open readonly')
-  ncBckID = nc_open_file_readwrite(next_outfile, 'update_mpas_states - open readwrite')
+  ncAnlID = nc_open_file_readonly(next_infile, 'update_mpas_states - open readonly') ! analysis from DART
+  ncBckID = nc_open_file_readwrite(next_outfile, 'update_mpas_states - open readwrite') ! background, original mpas file
 
   model_time = get_analysis_time(ncBckID, next_outfile)
   state_time = get_analysis_time(ncAnlID, next_infile)
@@ -105,24 +111,81 @@ program update_mpas_states
 
       allocate(variable(get_variable_size(1, i)))
 
-      if (get_variable_name(1,i) == 'uReconstructZonal' .or. &
-          get_variable_name(1,i) == 'uReconstructMeridional'.or. &
-          get_variable_name(1,i) == 'u') cycle varloop
+      if (get_variable_name(dom_id,i) == 'uReconstructZonal' .or. &
+          get_variable_name(dom_id,i) == 'uReconstructMeridional'.or. &
+          get_variable_name(dom_id,i) == 'u') cycle varloop
 
-      call nc_get_variable(ncAnlID, get_variable_name(1,i), variable)
-      call nc_put_variable(ncBckID, get_variable_name(1,i), variable)
+      call nc_get_variable(ncAnlID, get_variable_name(dom_id,i), variable)
+      call nc_put_variable(ncBckID, get_variable_name(dom_id,i), variable)
 
       deallocate(variable)
 
   enddo varloop
 
   ! deal with wind
   if (update_u_from_reconstruct) then
-     ! reconstruct u 
-     call error_handler(E_ERR,'update_mpas_states','update_u_from_reconstruct is not implemented',source)
+
+     if (use_increments_for_u_update) then
+        !  Read in the previous reconstructed winds from the original mpas netcdf file
+        !  and compute what increments (changes in values) were added by the assimilation.
+        !  Read in the original edge normal wind 'u' field from that same mpas netcdf
+        !  file and add the interpolated increments to compute the updated 'u' values.
+
+        ! read in u, uReconstrtuctZonal, uReconstructMeridional from background
+        ! read in uReconstrtuctZonal, uReconstructMeridional from analysis
+
+        dims = get_dim_lengths(dom_id, get_varid_from_varname(dom_id, 'u'))
+        allocate(u(dims(1), dims(2)), increments(dims(1), dims(2)))
+        call nc_get_variable(ncBckID, 'u', u)
+
+        dims = get_dim_lengths(dom_id, get_varid_from_varname(dom_id, 'uReconstructZonal'))
+        allocate(ucell(dims(1), dims(2)), ucell_dart(dims(1), dims(2)))
+        call nc_get_variable(ncBckID, 'uReconstructZonal', ucell)
+        call nc_get_variable(ncAnlID, 'uReconstructZonal', ucell_dart)
+        ucell = ucell_dart - ucell ! u increments
+
+        dims = get_dim_lengths(dom_id, get_varid_from_varname(dom_id, 'uReconstructMeridional'))
+        allocate(vcell(dims(1), dims(2)), vcell_dart(dims(1), dims(2)))
+        call nc_get_variable(ncBckID, 'uReconstructMeridional', vcell)
+        call nc_get_variable(ncAnlID, 'uReconstructMeridional', vcell_dart)
+        vcell = vcell_dart - vcell ! v increments
+   
+        call uv_cell_to_edges(ucell, vcell, increments) !
+
+        u = u + increments
+
+        call nc_put_variable(ncBckID, 'u', u)
+        call nc_put_variable(ncBckID, 'uReconstructZonal', ucell_dart)
+        call nc_put_variable(ncBckID, 'uReconstructMeridional', vcell_dart)
+
+        deallocate(u, increments, ucell, vcell, ucell_dart, vcell_dart)
+         
+     else
+        ! The state vector has updated zonal and meridional wind components.
+        ! put them directly into the arrays.  These are the full values, not
+        ! just increments.
+        dims = get_dim_lengths(dom_id, get_varid_from_varname(dom_id, 'u'))
+        allocate(u(dims(1), dims(2)))
+        dims = get_dim_lengths(dom_id, get_varid_from_varname(dom_id, 'uReconstructZonal'))
+        allocate(ucell_dart(dims(1), dims(2)))
+        call nc_get_variable(ncAnlID, 'uReconstructZonal', ucell_dart)
+        dims = get_dim_lengths(dom_id, get_varid_from_varname(dom_id, 'uReconstructMeridional'))
+        allocate(vcell_dart(dims(1), dims(2)))
+        call nc_get_variable(ncAnlID, 'uReconstructMeridional', vcell_dart)
+ 
+        call uv_cell_to_edges(ucell_dart, vcell_dart, u, .true.)
+
+        call nc_put_variable(ncBckID, 'u', variable)
+        call nc_put_variable(ncBckID, 'uReconstructZonal', variable)
+        call nc_put_variable(ncBckID, 'uReconstructMeridional', variable)
+
+        deallocate(u, ucell_dart, vcell_dart)
+
+     endif
+
   else
      ! copy u from analysis to background
-    allocate(variable(get_variable_size(1, get_varid_from_varname(1, 'u'))))
+    allocate(variable(get_variable_size(1, get_varid_from_varname(dom_id, 'u'))))
 
     call nc_get_variable(ncAnlID, 'u', variable)
     call nc_put_variable(ncBckID, 'u', variable)