Merge pull request #1 from adcroft/user/aja/diag_interp_and_extrinsic

User/aja/diag interp and extrinsic

src/ALE/MOM_ALE.F90

@@ -952,7 +952,7 @@ subroutine ALE_initRegridding(GV, max_depth, param_file, mod, regridCS, dz )
   call initialize_regridding( GV%ke, coordMode, interpScheme, regridCS, &
                               compressibility_fraction=compress_fraction )
 
-  if (coordMode(1:2) == 'Z*') then
+  if (coordinateMode(coordMode) == REGRIDDING_ZSTAR) then
     call get_param(param_file, mod, "ZSTAR_RIGID_SURFACE_THRESHOLD", height_of_rigid_surface, &
                  "A threshold height used to detect the presence of a rigid-surface\n"//&
                  "depressing the upper-surface of the model, such as an ice-shelf.\n"//&

src/ALE/MOM_regridding.F90

@@ -863,9 +863,8 @@ subroutine build_rho_grid( G, GV, h, tv, dzInterface, remapCS, CS )
       ! Local depth (G%bathyT is positive)
       nominalDepth = G%bathyT(i,j)*GV%m_to_H
 
-      call build_rho_column(CS, remapCS, tv%eqn_of_state, nz, nominalDepth, &
-                            h(i, j, :)*GV%H_to_m, &
-                            tv%T(i, j, :), tv%S(i, j, :), zNew)
+      call build_rho_column(CS, remapCS, nz, nominalDepth, h(i, j, :)*GV%H_to_m, &
+                            tv%T(i, j, :), tv%S(i, j, :), tv%eqn_of_state, zNew)
 
       if (CS%integrate_downward_for_e) then
         zOld(1) = 0.
@@ -926,7 +925,7 @@ subroutine build_rho_grid( G, GV, h, tv, dzInterface, remapCS, CS )
 end subroutine build_rho_grid
 
 
-subroutine build_rho_column(CS, remapCS, eqn_of_state, nz, depth, h, T, S, zInterface)
+subroutine build_rho_column(CS, remapCS, nz, depth, h, T, S, eqn_of_state, zInterface)
 ! The algorithn operates as follows within each
 ! column:
 ! 1. Given T & S within each layer, the layer densities are computed.
@@ -941,11 +940,11 @@ subroutine build_rho_column(CS, remapCS, eqn_of_state, nz, depth, h, T, S, zInte
 
   type(regridding_CS),   intent(in)    :: CS !< Regridding control structure
   type(remapping_CS),    intent(in)    :: remapCS !< Remapping parameters and options
-  type(EOS_type),        pointer       :: eqn_of_state !< Equation of state structure
   integer,               intent(in)    :: nz !< Number of levels
   real,                  intent(in)    :: depth !< Depth of ocean bottom (positive in m)
   real, dimension(nz),   intent(in)    :: h  !< Layer thicknesses, in m
   real, dimension(nz),   intent(in)    :: T, S !< T and S for column
+  type(EOS_type),        pointer       :: eqn_of_state !< Equation of state structure
   real, dimension(nz+1), intent(inout) :: zInterface !< Absolute positions of interfaces
 
   ! Local variables

src/ALE/regrid_consts.F90

@@ -8,22 +8,37 @@ module regrid_consts
 
 implicit none ; public
 
-! List of regridding types
-integer, parameter :: REGRIDDING_LAYER     = -1     !< Layer mode
-integer, parameter :: REGRIDDING_ZSTAR     = 0      !< z* coordinates
-integer, parameter :: REGRIDDING_RHO       = 1      !< Target interface densities
-integer, parameter :: REGRIDDING_SIGMA     = 2      !< Sigma coordinates
-integer, parameter :: REGRIDDING_ARBITRARY = 3      !< Arbitrary coordinates
-integer, parameter :: REGRIDDING_HYCOM1    = 4      !< Simple HyCOM coordinates without BBL
-integer, parameter :: REGRIDDING_SLIGHT    = 5      !< Stretched coordinates in the
+integer, parameter :: REGRIDDING_NUM_TYPES  = 7
+
+! List of regridding types. These should be consecutive and starting at 1.
+! This allows them to be used as array indices.
+integer, parameter :: REGRIDDING_LAYER     = 1      !< Layer mode
+integer, parameter :: REGRIDDING_ZSTAR     = 2      !< z* coordinates
+integer, parameter :: REGRIDDING_RHO       = 3      !< Target interface densities
+integer, parameter :: REGRIDDING_SIGMA     = 4      !< Sigma coordinates
+integer, parameter :: REGRIDDING_ARBITRARY = 5      !< Arbitrary coordinates
+integer, parameter :: REGRIDDING_HYCOM1    = 6      !< Simple HyCOM coordinates without BBL
+integer, parameter :: REGRIDDING_SLIGHT    = REGRIDDING_NUM_TYPES !< Stretched coordinates in the
                                                     !! lightest water, isopycnal below
-character(len=5), parameter :: REGRIDDING_LAYER_STRING = "LAYER"   !< Layer string
-character(len=2), parameter :: REGRIDDING_ZSTAR_STRING = "Z*"      !< z* string
-character(len=3), parameter :: REGRIDDING_RHO_STRING   = "RHO"     !< Rho string
-character(len=5), parameter :: REGRIDDING_SIGMA_STRING = "SIGMA"   !< Sigma string
+character(len=6), parameter :: REGRIDDING_LAYER_STRING = "LAYER"   !< Layer string
+character(len=6), parameter :: REGRIDDING_ZSTAR_STRING_OLD = "Z*"  !< z* string (legacy name)
+character(len=6), parameter :: REGRIDDING_ZSTAR_STRING = "ZSTAR"   !< z* string
+character(len=6), parameter :: REGRIDDING_RHO_STRING   = "RHO"     !< Rho string
+character(len=6), parameter :: REGRIDDING_SIGMA_STRING = "SIGMA"   !< Sigma string
+character(len=6), parameter :: REGRIDDING_ARBITRARY_STRING = "ARB" !< Arbitrary coordinates
 character(len=6), parameter :: REGRIDDING_HYCOM1_STRING = "HYCOM1" !< Hycom string
 character(len=6), parameter :: REGRIDDING_SLIGHT_STRING = "SLIGHT" !< Hybrid S-rho string
-character(len=5), parameter :: DEFAULT_COORDINATE_MODE = REGRIDDING_LAYER_STRING !< Default coordinate mode
+character(len=6), parameter :: DEFAULT_COORDINATE_MODE = REGRIDDING_LAYER_STRING !< Default coordinate mode
+
+integer, dimension(REGRIDDING_NUM_TYPES), parameter :: vertical_coords = &
+  (/ REGRIDDING_LAYER, REGRIDDING_ZSTAR, REGRIDDING_RHO, &
+    REGRIDDING_SIGMA, REGRIDDING_ARBITRARY, &
+    REGRIDDING_HYCOM1, REGRIDDING_SLIGHT /)
+
+character(len=*), dimension(REGRIDDING_NUM_TYPES), parameter :: vertical_coord_strings = &
+  (/ REGRIDDING_LAYER_STRING, REGRIDDING_ZSTAR_STRING, REGRIDDING_RHO_STRING, &
+    REGRIDDING_SIGMA_STRING, REGRIDDING_ARBITRARY_STRING, &
+    REGRIDDING_HYCOM1_STRING, REGRIDDING_SLIGHT_STRING /)
 
 interface coordinateUnits
   module procedure coordinateUnitsI
@@ -40,10 +55,12 @@ function coordinateMode(string)
   select case ( uppercase(trim(string)) )
     case (trim(REGRIDDING_LAYER_STRING)); coordinateMode = REGRIDDING_LAYER
     case (trim(REGRIDDING_ZSTAR_STRING)); coordinateMode = REGRIDDING_ZSTAR
+    case (trim(REGRIDDING_ZSTAR_STRING_OLD)); coordinateMode = REGRIDDING_ZSTAR
     case (trim(REGRIDDING_RHO_STRING));   coordinateMode = REGRIDDING_RHO
     case (trim(REGRIDDING_SIGMA_STRING)); coordinateMode = REGRIDDING_SIGMA
     case (trim(REGRIDDING_HYCOM1_STRING)); coordinateMode = REGRIDDING_HYCOM1
     case (trim(REGRIDDING_SLIGHT_STRING)); coordinateMode = REGRIDDING_SLIGHT
+    case (trim(REGRIDDING_ARBITRARY_STRING)); coordinateMode = REGRIDDING_ARBITRARY
     case default ; call MOM_error(FATAL, "coordinateMode: "//&
        "Unrecognized choice of coordinate ("//trim(string)//").")
   end select

src/core/MOM.F90

@@ -33,7 +33,7 @@ module MOM
 use MOM_diag_mediator,        only : diag_mediator_init, enable_averaging
 use MOM_diag_mediator,        only : diag_mediator_infrastructure_init
 use MOM_diag_mediator,        only : diag_register_area_ids
-use MOM_diag_mediator,        only : diag_set_thickness_ptr, diag_update_target_grids
+use MOM_diag_mediator,        only : diag_set_state_ptrs, diag_update_remap_grids
 use MOM_diag_mediator,        only : disable_averaging, post_data, safe_alloc_ptr
 use MOM_diag_mediator,        only : register_diag_field, register_static_field
 use MOM_diag_mediator,        only : register_scalar_field
@@ -737,7 +737,7 @@ subroutine step_MOM(fluxes, state, Time_start, time_interval, CS)
 
         ! Whenever thickness changes let the diag manager know, target grids
         ! for vertical remapping may need to be regenerated.
-        call diag_update_target_grids(CS%diag)
+        call diag_update_remap_grids(CS%diag)
 
         call post_diags_TS_vardec(G, CS, dtdia)
 
@@ -819,7 +819,7 @@ subroutine step_MOM(fluxes, state, Time_start, time_interval, CS)
 
         ! Whenever thickness changes let the diag manager know, target grids
         ! for vertical remapping may need to be regenerated.
-        call diag_update_target_grids(CS%diag)
+        call diag_update_remap_grids(CS%diag)
 
       endif
     endif
@@ -933,7 +933,7 @@ subroutine step_MOM(fluxes, state, Time_start, time_interval, CS)
 
     ! Whenever thickness changes let the diag manager know, target grids
     ! for vertical remapping may need to be regenerated.
-    call diag_update_target_grids(CS%diag)
+    call diag_update_remap_grids(CS%diag)
 
     if (CS%useMEKE) call step_forward_MEKE(CS%MEKE, h, CS%VarMix%SN_u, CS%VarMix%SN_v, &
                                            CS%visc, dt, G, GV, CS%MEKE_CSp, CS%uhtr, CS%vhtr)
@@ -1064,7 +1064,7 @@ subroutine step_MOM(fluxes, state, Time_start, time_interval, CS)
         ! Whenever thickness changes let the diag manager know, target grids
         ! for vertical remapping may need to be regenerated. This needs to
         ! happen after the H update and before the next post_data.
-        call diag_update_target_grids(CS%diag)
+        call diag_update_remap_grids(CS%diag)
 
         call post_diags_TS_vardec(G, CS, CS%dt_trans)
 
@@ -1942,17 +1942,18 @@ subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in)
   ! and before MOM_diagnostics_init
   call diag_masks_set(G, GV%ke, CS%missing, diag)
 
-  ! Set up a pointers h within diag mediator control structure,
-  ! this needs to occur _after_ CS%h has been allocated.
-  call diag_set_thickness_ptr(CS%h, diag)
+  ! Set up pointers within diag mediator control structure,
+  ! this needs to occur _after_ CS%h etc. have been allocated.
+  call diag_set_state_ptrs(CS%h, CS%T, CS%S, CS%tv%eqn_of_state, diag)
 
   ! This call sets up the diagnostic axes. These are needed,
   ! e.g. to generate the target grids below.
   call set_axes_info(G, GV, param_file, diag)
 
-  ! Whenever thickness changes let the diag manager know, target grids
-  ! for vertical remapping may need to be regenerated. 
-  call diag_update_target_grids(diag)
+  ! Whenever thickness/T/S changes let the diag manager know, target grids
+  ! for vertical remapping may need to be regenerated.
+  ! FIXME: are h, T, S updated at the same time? Review these for T, S updates.
+  call diag_update_remap_grids(diag)
 
   ! Diagnose static fields AND associate areas/volumes with axes
   call write_static_fields(G, CS%diag)
@@ -2216,7 +2217,7 @@ subroutine register_diags(Time, G, GV, CS, ADp)
       'Meridional velocity', 'meter second-1', cmor_field_name='vo', cmor_units='m s-1', &
       cmor_standard_name='sea_water_y_velocity', cmor_long_name='Sea Water Y Velocity')
   CS%id_h = register_diag_field('ocean_model', 'h', diag%axesTL, Time, &
-      'Layer Thickness', thickness_units, v_cell_method='sum')
+      'Layer Thickness', thickness_units, v_extrinsic=.true.)
 
   CS%id_volo = register_scalar_field('ocean_model', 'volo', Time, diag,&
       long_name='Total volume of liquid ocean', units='m3',            &
@@ -2366,7 +2367,7 @@ subroutine register_diags(Time, G, GV, CS, ADp)
   CS%id_v_predia = register_diag_field('ocean_model', 'v_predia', diag%axesCvL, Time, &
       'Meridional velocity before diabatic forcing', 'meter second-1')
   CS%id_h_predia = register_diag_field('ocean_model', 'h_predia', diag%axesTL, Time, &
-      'Layer Thickness before diabatic forcing', thickness_units, v_cell_method='sum')
+      'Layer Thickness before diabatic forcing', thickness_units, v_extrinsic=.true.)
   CS%id_e_predia = register_diag_field('ocean_model', 'e_predia', diag%axesTi, Time, &
       'Interface Heights before diabatic forcing', 'meter')
   if (CS%diabatic_first .and. (.not. CS%adiabatic)) then
@@ -2375,7 +2376,7 @@ subroutine register_diags(Time, G, GV, CS, ADp)
     CS%id_v_preale = register_diag_field('ocean_model', 'v_preale', diag%axesCvL, Time, &
         'Meridional velocity before remapping', 'meter second-1')
     CS%id_h_preale = register_diag_field('ocean_model', 'h_preale', diag%axesTL, Time, &
-        'Layer Thickness before remapping', thickness_units, v_cell_method='sum')
+        'Layer Thickness before remapping', thickness_units, v_extrinsic=.true.)
     CS%id_T_preale = register_diag_field('ocean_model', 'T_preale', diag%axesTL, Time, &
         'Temperature before remapping', 'degC')
     CS%id_S_preale = register_diag_field('ocean_model', 'S_preale', diag%axesTL, Time, &
@@ -2428,7 +2429,7 @@ subroutine register_diags_TS_tendency(Time, G, CS)
       cmor_field_name="opottemptend", cmor_units="W m-2",                                         &
       cmor_standard_name="tendency_of_sea_water_potential_temperature_expressed_as_heat_content", &
       cmor_long_name ="Tendency of Sea Water Potential Temperature Expressed as Heat Content",    &
-      v_cell_method='sum')
+      v_extrinsic=.true.)
   CS%id_Th_tendency_2d = register_diag_field('ocean_model', 'Th_tendency_2d', diag%axesT1, Time,              &
       'Vertical sum of net time tendency for heat', 'W/m2',                                                   &
       cmor_field_name="opottemptend_2d", cmor_units="W m-2",                                                   &
@@ -2468,7 +2469,7 @@ subroutine register_diags_TS_tendency(Time, G, CS)
       cmor_field_name="osalttend", cmor_units="kg m-2 s-1",                               &
       cmor_standard_name="tendency_of_sea_water_salinity_expressed_as_salt_content",      &
       cmor_long_name ="Tendency of Sea Water Salinity Expressed as Salt Content",         &
-      v_cell_method='sum')
+      v_extrinsic=.true.)
   CS%id_Sh_tendency_2d = register_diag_field('ocean_model', 'Sh_tendency_2d', diag%axesT1, Time, &
       'Vertical sum of net time tendency for salt', 'kg/(m2 * s)',                               &
       cmor_field_name="osalttend_2d", cmor_units="kg m-2 s-1",                                   &

src/core/MOM_dynamics_legacy_split.F90

@@ -78,7 +78,7 @@ module MOM_dynamics_legacy_split
 use MOM_diag_mediator, only : diag_mediator_init, enable_averaging
 use MOM_diag_mediator, only : disable_averaging, post_data, safe_alloc_ptr
 use MOM_diag_mediator, only : register_diag_field, register_static_field
-use MOM_diag_mediator, only : set_diag_mediator_grid, diag_ctrl, diag_update_target_grids
+use MOM_diag_mediator, only : set_diag_mediator_grid, diag_ctrl, diag_update_remap_grids
 use MOM_domains, only : MOM_domains_init, pass_var, pass_vector
 use MOM_domains, only : pass_var_start, pass_var_complete
 use MOM_domains, only : pass_vector_start, pass_vector_complete
@@ -946,7 +946,7 @@ subroutine step_MOM_dyn_legacy_split(u, v, h, tv, visc, &
     call cpu_clock_end(id_clock_continuity)
     ! Whenever thickness changes let the diag manager know, target grids
     ! for vertical remapping may need to be regenerated.
-    call diag_update_target_grids(CS%diag)
+    call diag_update_remap_grids(CS%diag)
     if (G%nonblocking_updates) then
       call cpu_clock_begin(id_clock_pass)
       pid_h = pass_var_start(h, G%Domain)
@@ -987,7 +987,7 @@ subroutine step_MOM_dyn_legacy_split(u, v, h, tv, visc, &
     call cpu_clock_end(id_clock_continuity)
     ! Whenever thickness changes let the diag manager know, target grids
     ! for vertical remapping may need to be regenerated.
-    call diag_update_target_grids(CS%diag)
+    call diag_update_remap_grids(CS%diag)
     call cpu_clock_begin(id_clock_pass)
     call pass_var(h, G%Domain)
     call cpu_clock_end(id_clock_pass)

src/core/MOM_dynamics_split_RK2.F90

@@ -15,7 +15,7 @@ module MOM_dynamics_split_RK2
 use MOM_diag_mediator,     only : diag_mediator_init, enable_averaging
 use MOM_diag_mediator,     only : disable_averaging, post_data, safe_alloc_ptr
 use MOM_diag_mediator,     only : register_diag_field, register_static_field
-use MOM_diag_mediator,     only : set_diag_mediator_grid, diag_ctrl, diag_update_target_grids
+use MOM_diag_mediator,     only : set_diag_mediator_grid, diag_ctrl, diag_update_remap_grids
 use MOM_domains,           only : MOM_domains_init
 use MOM_domains,           only : To_South, To_West, To_All, CGRID_NE, SCALAR_PAIR
 use MOM_domains,           only : create_group_pass, do_group_pass, group_pass_type
@@ -846,7 +846,7 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, &
   call cpu_clock_end(id_clock_pass)
   ! Whenever thickness changes let the diag manager know, target grids
   ! for vertical remapping may need to be regenerated.
-  call diag_update_target_grids(CS%diag)
+  call diag_update_remap_grids(CS%diag)
   if (showCallTree) call callTree_wayPoint("done with continuity (step_MOM_dyn_split_RK2)")
 
   call cpu_clock_begin(id_clock_pass)
@@ -900,8 +900,8 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, &
   if (CS%id_vav        > 0) call post_data(CS%id_vav, v_av,                 CS%diag)
   if (CS%id_u_BT_accel > 0) call post_data(CS%id_u_BT_accel, CS%u_accel_bt, CS%diag)
   if (CS%id_v_BT_accel > 0) call post_data(CS%id_v_BT_accel, CS%v_accel_bt, CS%diag)
-  if (CS%id_umo        > 0) call post_data(CS%id_umo, uh*GV%H_to_kg_m2,      CS%diag)
-  if (CS%id_vmo        > 0) call post_data(CS%id_vmo, vh*GV%H_to_kg_m2,      CS%diag)
+  if (CS%id_umo        > 0) call post_data(CS%id_umo, uh,      CS%diag)
+  if (CS%id_vmo        > 0) call post_data(CS%id_vmo, vh,      CS%diag)
 
   ! depth summed zonal mass transport
   if (CS%id_umo_2d > 0) then
@@ -1203,23 +1203,25 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, param_fil
 
   flux_units = get_flux_units(GV)
   CS%id_uh = register_diag_field('ocean_model', 'uh', diag%axesCuL, Time, &
-      'Zonal Thickness Flux', flux_units)
+      'Zonal Thickness Flux', flux_units, y_cell_method='sum', v_extrinsic=.true.)
   CS%id_vh = register_diag_field('ocean_model', 'vh', diag%axesCvL, Time, &
-      'Meridional Thickness Flux', flux_units)
+      'Meridional Thickness Flux', flux_units, x_cell_method='sum', v_extrinsic=.true.)
   CS%id_umo = register_diag_field('ocean_model', 'umo',                       &
       diag%axesCuL, Time,'Zonal Mass Transport (including SGS param)', 'kg/s',&
-      cmor_standard_name='ocean_mass_x_transport', cmor_long_name='Ocean Mass X Transport')
+      cmor_standard_name='ocean_mass_x_transport', cmor_long_name='Ocean Mass X Transport', &
+      conversion=GV%H_to_kg_m2, y_cell_method='sum', v_extrinsic=.true.)
   CS%id_vmo = register_diag_field('ocean_model', 'vmo',                            &
       diag%axesCvL, Time,'Meridional Mass Transport (including SGS param)', 'kg/s',&
-      cmor_standard_name='ocean_mass_y_transport', cmor_long_name='Ocean Mass Y Transport')
+      cmor_standard_name='ocean_mass_y_transport', cmor_long_name='Ocean Mass Y Transport', &
+      conversion=GV%H_to_kg_m2, x_cell_method='sum', v_extrinsic=.true.)
   CS%id_umo_2d = register_diag_field('ocean_model', 'umo_2d',                              &
       diag%axesCu1, Time,'Zonal Mass Transport (including SGS param) Vertical Sum', 'kg/s',&
       cmor_standard_name='ocean_mass_x_transport_vertical_sum',                            &
-      cmor_long_name='Ocean Mass X Transport Vertical Sum')
+      cmor_long_name='Ocean Mass X Transport Vertical Sum', y_cell_method='sum')
   CS%id_vmo_2d = register_diag_field('ocean_model', 'vmo_2d',                                   &
       diag%axesCv1, Time,'Meridional Mass Transport (including SGS param) Vertical Sum', 'kg/s',&
       cmor_standard_name='ocean_mass_y_transport_vertical_sum',                                 &
-      cmor_long_name='Ocean Mass Y Transport Vertical Sum')
+      cmor_long_name='Ocean Mass Y Transport Vertical Sum', x_cell_method='sum')
 
   CS%id_CAu = register_diag_field('ocean_model', 'CAu', diag%axesCuL, Time, &
       'Zonal Coriolis and Advective Acceleration', 'meter second-2')