Merge pull request #66 from NOAA-GFDL/dev/gfdl

Sync with NOAA-GFDL dev/gfdl

.testing/tc0/MOM_input

@@ -138,6 +138,9 @@ THICKNESS_CONFIG = "uniform"    !
 
 ! === module MOM_diag_mediator ===
 
+USE_GRID_SPACE_DIAG_COORDINATE_AXES = True !   [Boolean] default = False
+                                ! If true, use a grid index coordinate convention for diagnostic axes.
+
 ! === module MOM_MEKE ===
 
 ! === module MOM_lateral_mixing_coeffs ===

src/framework/MOM_diag_mediator.F90

@@ -243,7 +243,7 @@ module MOM_diag_mediator
   integer :: chksum_iounit = -1           !< The unit number of a diagnostic documentation file.
                                           !! This file is open if available_diag_doc_unit is > 0.
   logical :: diag_as_chksum !< If true, log chksums in a text file instead of posting diagnostics
-
+  logical :: grid_space_axes !< If true, diagnostic horizontal coordinates axes are in grid space.
 ! The following fields are used for the output of the data.
   integer :: is  !< The start i-index of cell centers within the computational domain
   integer :: ie  !< The end i-index of cell centers within the computational domain
@@ -359,25 +359,71 @@ subroutine set_axes_info(G, GV, US, param_file, diag_cs, set_vertical)
   integer :: i, j, k, nz
   real :: zlev(GV%ke), zinter(GV%ke+1)
   logical :: set_vert
+  real, allocatable, dimension(:) :: IaxB,iax
+  real, allocatable, dimension(:) :: JaxB,jax
+
 
   set_vert = .true. ; if (present(set_vertical)) set_vert = set_vertical
 
+
+  if (diag_cs%grid_space_axes) then
+     allocate(IaxB(G%IsgB:G%IegB))
+     do i=G%IsgB, G%IegB
+       Iaxb(i)=real(i)
+     enddo
+     allocate(iax(G%isg:G%ieg))
+     do i=G%isg, G%ieg
+       iax(i)=real(i)-0.5
+     enddo
+     allocate(JaxB(G%JsgB:G%JegB))
+     do j=G%JsgB, G%JegB
+       JaxB(j)=real(j)
+     enddo
+     allocate(jax(G%jsg:G%jeg))
+     do j=G%jsg, G%jeg
+       jax(j)=real(j)-0.5
+     enddo
+  endif
+
   ! Horizontal axes for the native grids
   if (G%symmetric) then
-    id_xq = diag_axis_init('xq', G%gridLonB(G%isgB:G%iegB), G%x_axis_units, 'x', &
-              'q point nominal longitude', Domain2=G%Domain%mpp_domain, domain_position=EAST)
-    id_yq = diag_axis_init('yq', G%gridLatB(G%jsgB:G%jegB), G%y_axis_units, 'y', &
-              'q point nominal latitude', Domain2=G%Domain%mpp_domain, domain_position=NORTH)
+     if (diag_cs%grid_space_axes) then
+        id_xq = diag_axis_init('iq', IaxB(G%isgB:G%iegB), 'none', 'x', &
+             'q point grid-space longitude', Domain2=G%Domain%mpp_domain, domain_position=EAST)
+        id_yq = diag_axis_init('jq', JaxB(G%jsgB:G%jegB), 'none', 'y', &
+         'q point grid space latitude', Domain2=G%Domain%mpp_domain, domain_position=NORTH)
+     else
+        id_xq = diag_axis_init('xq', G%gridLonB(G%isgB:G%iegB), G%x_axis_units, 'x', &
+             'q point nominal longitude', Domain2=G%Domain%mpp_domain, domain_position=EAST)
+        id_yq = diag_axis_init('yq', G%gridLatB(G%jsgB:G%jegB), G%y_axis_units, 'y', &
+            'q point nominal latitude', Domain2=G%Domain%mpp_domain, domain_position=NORTH)
+     endif
   else
-    id_xq = diag_axis_init('xq', G%gridLonB(G%isg:G%ieg), G%x_axis_units, 'x', &
-              'q point nominal longitude', Domain2=G%Domain%mpp_domain, domain_position=EAST)
-    id_yq = diag_axis_init('yq', G%gridLatB(G%jsg:G%jeg), G%y_axis_units, 'y', &
-              'q point nominal latitude', Domain2=G%Domain%mpp_domain, domain_position=NORTH)
+     if (diag_cs%grid_space_axes) then
+        id_xq = diag_axis_init('Iq', IaxB(G%isg:G%ieg), 'none', 'x', &
+             'q point grid-space longitude', Domain2=G%Domain%mpp_domain, domain_position=EAST)
+        id_yq = diag_axis_init('Jq', JaxB(G%jsg:G%jeg), 'none', 'y', &
+             'q point grid space latitude', Domain2=G%Domain%mpp_domain, domain_position=NORTH)
+     else
+        id_xq = diag_axis_init('xq', G%gridLonB(G%isg:G%ieg), G%x_axis_units, 'x', &
+             'q point nominal longitude', Domain2=G%Domain%mpp_domain, domain_position=EAST)
+        id_yq = diag_axis_init('yq', G%gridLatB(G%jsg:G%jeg), G%y_axis_units, 'y', &
+             'q point nominal latitude', Domain2=G%Domain%mpp_domain, domain_position=NORTH)
+     endif
+  endif
+
+
+  if (diag_cs%grid_space_axes) then
+     id_xh = diag_axis_init('ih', iax(G%isg:G%ieg), 'none', 'x', &
+          'h point grid-space longitude', Domain2=G%Domain%mpp_domain, domain_position=EAST)
+     id_yh = diag_axis_init('jh', jax(G%jsg:G%jeg), 'none', 'y', &
+         'h point grid space latitude', Domain2=G%Domain%mpp_domain, domain_position=NORTH)
+  else
+     id_xh = diag_axis_init('xh', G%gridLonT(G%isg:G%ieg), G%x_axis_units, 'x', &
+          'h point nominal longitude', Domain2=G%Domain%mpp_domain)
+     id_yh = diag_axis_init('yh', G%gridLatT(G%jsg:G%jeg), G%y_axis_units, 'y', &
+          'h point nominal latitude', Domain2=G%Domain%mpp_domain)
   endif
-  id_xh = diag_axis_init('xh', G%gridLonT(G%isg:G%ieg), G%x_axis_units, 'x', &
-              'h point nominal longitude', Domain2=G%Domain%mpp_domain)
-  id_yh = diag_axis_init('yh', G%gridLatT(G%jsg:G%jeg), G%y_axis_units, 'y', &
-              'h point nominal latitude', Domain2=G%Domain%mpp_domain)
 
   if (set_vert) then
     nz = GV%ke
@@ -531,6 +577,9 @@ subroutine set_axes_info(G, GV, US, param_file, diag_cs, set_vertical)
     endif
   enddo
 
+  if (diag_cs%grid_space_axes) then
+     deallocate(IaxB,iax,JaxB,jax)
+  endif
   !Define the downsampled axes
   call set_axes_info_dsamp(G, GV, param_file, diag_cs, id_zl_native, id_zi_native)
 
@@ -3037,6 +3086,10 @@ subroutine diag_mediator_init(G, GV, US, nz, param_file, diag_cs, doc_file_dir)
                  "If true, use the order of arithmetic and expressions that recover the "//&
                  "answers from the end of 2018.  Otherwise, use updated and more robust "//&
                  "forms of the same expressions.", default=default_2018_answers)
+  call get_param(param_file, mdl, 'USE_GRID_SPACE_DIAGNOSTIC_AXES', diag_cs%grid_space_axes, &
+                 'If true, use a grid index coordinate convention for diagnostic axes. ',&
+                 default=.false.)
+
   if (diag_cs%num_diag_coords>0) then
     allocate(diag_coords(diag_cs%num_diag_coords))
     if (diag_cs%num_diag_coords==1) then ! The default is to provide just one instance of Z*
@@ -4264,4 +4317,3 @@ subroutine downsample_mask_3d(field_in, field_out, dl, isc_o, jsc_o, isc_d, iec_
 end subroutine downsample_mask_3d
 
 end module MOM_diag_mediator
-

src/initialization/MOM_state_initialization.F90

@@ -84,7 +84,7 @@ module MOM_state_initialization
 use dense_water_initialization, only : dense_water_initialize_TS
 use dense_water_initialization, only : dense_water_initialize_sponges
 use dumbbell_initialization, only : dumbbell_initialize_sponges
-use MOM_tracer_Z_init, only : find_interfaces, tracer_Z_init_array, determine_temperature
+use MOM_tracer_Z_init, only : tracer_Z_init_array, determine_temperature
 use MOM_ALE, only : ALE_initRegridding, ALE_CS, ALE_initThicknessToCoord
 use MOM_ALE, only : ALE_remap_scalar, ALE_build_grid, ALE_regrid_accelerated
 use MOM_ALE, only : TS_PLM_edge_values
@@ -2422,6 +2422,114 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, G, GV, US, PF, just_read_param
 
 end subroutine MOM_temp_salt_initialize_from_Z
 
+
+!> Find interface positions corresponding to interpolated depths in a density profile
+subroutine find_interfaces(rho, zin, nk_data, Rb, depth, zi, G, US, nlevs, nkml, hml, &
+                           eps_z, eps_rho)
+  type(ocean_grid_type),      intent(in)  :: G     !< The ocean's grid structure
+  integer,                    intent(in)  :: nk_data !< The number of levels in the input data
+  real, dimension(SZI_(G),SZJ_(G),nk_data), &
+                              intent(in)  :: rho   !< Potential density in z-space [R ~> kg m-3]
+  real, dimension(nk_data),   intent(in)  :: zin   !< Input data levels [Z ~> m].
+  real, dimension(SZK_(G)+1), intent(in)  :: Rb    !< target interface densities [R ~> kg m-3]
+  real, dimension(SZI_(G),SZJ_(G)), &
+                              intent(in)  :: depth !< ocean depth [Z ~> m].
+  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), &
+                              intent(out) :: zi    !< The returned interface heights [Z ~> m]
+  type(unit_scale_type),      intent(in)  :: US    !< A dimensional unit scaling type
+  integer, dimension(SZI_(G),SZJ_(G)), &
+                              intent(in)  :: nlevs !< number of valid points in each column
+  integer,                    intent(in)  :: nkml  !< number of mixed layer pieces to distribute over
+                                                   !! a depth of hml.
+  real,                       intent(in)  :: hml   !< mixed layer depth [Z ~> m].
+  real,                       intent(in)  :: eps_z !< A negligibly small layer thickness [Z ~> m].
+  real,                       intent(in)  :: eps_rho !< A negligibly small density difference [R ~> kg m-3].
+
+  ! Local variables
+  real, dimension(nk_data) :: rho_ ! A column of densities [R ~> kg m-3]
+  real, dimension(SZK_(G)+1) :: zi_ ! A column interface heights (negative downward) [Z ~> m].
+  real    :: slope      ! The rate of change of height with density [Z R-1 ~> m4 kg-1]
+  real    :: drhodz     ! A local vertical density gradient [R Z-1 ~> kg m-4]
+  real, parameter :: zoff=0.999
+  logical :: unstable   ! True if the column is statically unstable anywhere.
+  integer :: nlevs_data ! The number of data values in a column.
+  logical :: work_down  ! This indicates whether this pass goes up or down the water column.
+  integer :: k_int, lo_int, hi_int, mid
+  integer :: i, j, k, is, ie, js, je, nz
+
+  is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke
+
+  zi(:,:,:) = 0.0
+
+  do j=js,je ; do i=is,ie
+    nlevs_data = nlevs(i,j)
+    do k=1,nlevs_data ; rho_(k) = rho(i,j,k) ; enddo
+
+    unstable=.true.
+    work_down = .true.
+    do while (unstable)
+      ! Modifiy the input profile until it no longer has densities that decrease with depth.
+      unstable=.false.
+      if (work_down) then
+        do k=2,nlevs_data-1 ; if (rho_(k) - rho_(k-1) < 0.0 ) then
+          if (k == 2) then
+            rho_(k-1) = rho_(k) - eps_rho
+          else
+            drhodz = (rho_(k+1)-rho_(k-1)) / (zin(k+1)-zin(k-1))
+            if (drhodz < 0.0) unstable=.true.
+            rho_(k) = rho_(k-1) + drhodz*zoff*(zin(k)-zin(k-1))
+          endif
+        endif ; enddo
+        work_down = .false.
+      else
+        do k=nlevs_data-1,2,-1 ;  if (rho_(k+1) - rho_(k) < 0.0) then
+          if (k == nlevs_data-1) then
+            rho_(k+1) = rho_(k-1) + eps_rho !### This should be rho_(k) + eps_rho
+          else
+            drhodz = (rho_(k+1)-rho_(k-1)) / (zin(k+1)-zin(k-1))
+            if (drhodz  < 0.0) unstable=.true.
+            rho_(k) = rho_(k+1) - drhodz*(zin(k+1)-zin(k))
+          endif
+        endif ; enddo
+        work_down = .true.
+      endif
+    enddo
+
+    ! Find and store the interface depths.
+    zi_(1) = 0.0
+    do K=2,nz
+      ! Find the value of k_int in the list of rho_ where rho_(k_int) <= Rb(K) < rho_(k_int+1).
+      ! This might be made a little faster by exploiting the fact that Rb is
+      ! monotonically increasing and not resetting lo_int back to 1 inside the K loop.
+      lo_int = 1 ; hi_int = nlevs_data
+      do while (lo_int < hi_int)
+        mid = (lo_int+hi_int) / 2
+        if (Rb(K) < rho_(mid)) then ; hi_int = mid
+        else ; lo_int = mid+1 ; endif
+      enddo
+      k_int = max(1, lo_int-1)
+
+      ! Linearly interpolate to find the depth, zi_, where Rb would be found.
+      slope = (zin(k_int+1) - zin(k_int)) / max(rho_(k_int+1) - rho_(k_int), eps_rho)
+      zi_(K) = -1.0*(zin(k_int) + slope*(Rb(K)-rho_(k_int)))
+      zi_(K) = min(max(zi_(K), -depth(i,j)), -1.0*hml)
+    enddo
+    zi_(nz+1) = -depth(i,j)
+    if (nkml > 0) then ; do K=2,nkml+1
+      zi_(K) = max(hml*((1.0-real(K))/real(nkml)), -depth(i,j))
+    enddo ; endif
+    do K=nz,max(nkml+2,2),-1
+      if (zi_(K) < zi_(K+1) + eps_Z) zi_(K) = zi_(K+1) + eps_Z
+      if (zi_(K) > -1.0*hml)  zi_(K) = max(-1.0*hml, -depth(i,j))
+    enddo
+
+    do K=1,nz+1
+      zi(i,j,K) = zi_(K)
+    enddo
+  enddo ; enddo ! i- and j- loops
+
+end subroutine find_interfaces
+
 !> Run simple unit tests
 subroutine MOM_state_init_tests(G, GV, US, tv)
   type(ocean_grid_type),     intent(inout) :: G    !< The ocean's grid structure.