Added thkcello diagnostic

- Corrected calculation of thkcello, cell thickness in meters, for non-Boussinesq mode. - Moved diagnostic of thkcello into MOM_diagnostics.F90 from MOM.F90. - Renamed a temporary 3d array which is re-used by various diagnostics. - Closes issue #98.
mom-ocean · Apr 27, 2015 · 00ba2d7 · 00ba2d7
1 parent 73bc292
commit 00ba2d7
Show file tree

Hide file tree

Showing 2 changed files with 48 additions and 16 deletions.
diff --git a/src/core/MOM.F90 b/src/core/MOM.F90
@@ -551,7 +551,6 @@ module MOM
   integer :: id_zos      = -1
   integer :: id_zossq    = -1
   integer :: id_volo     = -1
-  integer :: id_thkcello = -1 
   integer :: id_ssh      = -1
   integer :: id_sst      = -1 
   integer :: id_sst_sq   = -1
@@ -1369,7 +1368,6 @@ subroutine step_MOM(fluxes, state, Time_start, time_interval, CS)
     if (CS%id_u > 0) call post_data(CS%id_u, u, CS%diag)
     if (CS%id_v > 0) call post_data(CS%id_v, v, CS%diag)
     if (CS%id_h > 0) call post_data(CS%id_h, h, CS%diag)
-    if (CS%id_thkcello > 0) call post_data(CS%id_thkcello, G%H_to_m*h, CS%diag)
 
     ! compute ssh, which is either eta_av for Bouss, or 
     ! diagnosed ssh for non-Bouss; call "find_eta" for this 
@@ -2224,8 +2222,6 @@ subroutine register_diags(Time, G, CS, ADp)
   CS%id_h = register_diag_field('ocean_model', 'h', diag%axesTL, Time, &
       'Layer Thickness', thickness_units)
 
-  CS%id_thkcello = register_diag_field('ocean_model', 'thkcello', diag%axesTL, Time, &
-      long_name = 'Cell Thickness', standard_name='cell_thickness', units='m')
   CS%id_volo = register_scalar_field('ocean_model', 'volo', Time, diag,&
       long_name='Total volume of liquid ocean', units='m3',            &
       standard_name='sea_water_volume')

diff --git a/src/diagnostics/MOM_diagnostics.F90 b/src/diagnostics/MOM_diagnostics.F90
@@ -94,9 +94,7 @@ module MOM_diagnostics
     uh_Rlay   => NULL(), & ! zonal and meridional transports in layered 
     vh_Rlay   => NULL(), & ! potential rho coordinates: m3/s(Bouss) kg/s(non-Bouss)
     uhGM_Rlay => NULL(), & ! zonal and meridional Gent-McWilliams transports in layered 
-    vhGM_Rlay => NULL(), & ! potential density coordinates, m3/s (Bouss) kg/s(non-Bouss)
-
-    masscello => NULL()    ! mass per unit area of grid cell kg/m2
+    vhGM_Rlay => NULL()    ! potential density coordinates, m3/s (Bouss) kg/s(non-Bouss)
 
   ! following fields are 2-D.
   real, pointer, dimension(:,:) :: &
@@ -120,7 +118,8 @@ module MOM_diagnostics
     KE_adv     => NULL(),&  ! KE source from along-layer advection
     KE_visc    => NULL(),&  ! KE source from vertical viscosity
     KE_horvisc => NULL(),&  ! KE source from horizontal viscosity
-    KE_dia     => NULL()    ! KE source from diapycnal diffusion
+    KE_dia     => NULL(),&  ! KE source from diapycnal diffusion
+    diag_tmp3d => NULL()    ! 3D re-usable arrays for diagnostics
 
   ! diagnostic IDs
   integer :: id_e              = -1, id_e_D            = -1
@@ -143,6 +142,7 @@ module MOM_diagnostics
   integer :: id_sosga          = -1, id_tosga          = -1
   integer :: id_temp_layer_ave = -1, id_salt_layer_ave = -1
   integer :: id_pbo            = -1
+  integer :: id_thkcello       = -1
 
   type(wave_speed_CS), pointer :: wave_speed_CSp => NULL()  
 
@@ -203,6 +203,7 @@ subroutine calculate_diagnostic_fields(u, v, h, uh, vh, tv, ADp, CDp, fluxes, &
 
   ! tmp array for surface properties 
   real :: surface_field(SZI_(G),SZJ_(G)) 
+  real :: pressure_1d(SZI_(G)) ! Temporary array for pressure when calling EOS
 
   real :: pres(SZI_(G))
   real :: wt, wt_p
@@ -262,20 +263,51 @@ subroutine calculate_diagnostic_fields(u, v, h, uh, vh, tv, ADp, CDp, fluxes, &
   ! mass per area of grid cell (for Bouss, use Rho0)
   if (CS%id_masscello > 0) then
     do k=1,nz; do j=js,je ; do i=is,ie 
-       CS%masscello(i,j,k) = G%H_to_kg_m2*h(i,j,k)
+       CS%diag_tmp3d(i,j,k) = G%H_to_kg_m2*h(i,j,k)
     enddo ; enddo ; enddo
-    call post_data(CS%id_masscello, CS%masscello, CS%diag)
+    call post_data(CS%id_masscello, CS%diag_tmp3d, CS%diag)
   endif
 
   ! mass of liquid ocean (for Bouss, use Rho0)
   if (CS%id_masso > 0) then
-    do k=1,nz; do j=js,je ; do i=is,ie 
-       CS%masscello(i,j,k) = G%H_to_kg_m2*h(i,j,k)*G%areaT(i,j)
+    do k=1,nz; do j=js,je ; do i=is,ie
+       CS%diag_tmp3d(i,j,k) = G%H_to_kg_m2*h(i,j,k)*G%areaT(i,j)
     enddo ; enddo ; enddo
-    masso = (reproducing_sum(sum(CS%masscello,3)))
+    masso = (reproducing_sum(sum(CS%diag_tmp3d,3)))
     call post_data(CS%id_masso, masso, CS%diag)
   endif
 
+  ! Thickness of cells in meters
+  if (CS%id_thkcello > 0) then
+    if (G%Boussinesq) then ! thkcello is just h, in m.
+      call post_data(CS%id_thkcello, G%H_to_m*h, CS%diag)
+    else ! non-Boussinesq, thkcello is the actual thickness in m.
+      do j=js,je
+        do i=is,ie
+          pressure_1d(i) = fluxes%p_surf(i,j) ! Pressure loading at the surface, in Pa
+        enddo
+        do k=1,nz
+          ! Pressure for EOS at the layer center, in Pa.
+          do i=is,ie
+            pressure_1d(i) = pressure_1d(i) + 0.5*(G%G_Earth*G%H_to_kg_m2)*h(i,j,k)
+          enddo
+          ! Store in-situ density in diag_tmp3d, in kg/m3.
+          call calculate_density(tv%T(:,j,k),tv%S(:,j,k), pressure_1d, &
+                                 CS%diag_tmp3d(:,j,k), is, ie-is+1, tv%eqn_of_state)
+          ! Infer cell thickness as dp(g*rho), in m. Store in diag_tmp3d
+          do i=is,ie
+            CS%diag_tmp3d(i,j,k) = (G%H_to_kg_m2*h(i,j,k))/CS%diag_tmp3d(i,j,k)
+          enddo
+          ! Pressure for EOS at the bottom interface, in Pa.
+          do i=is,ie
+            pressure_1d(i) = pressure_1d(i) + 0.5*(G%G_Earth*G%H_to_kg_m2)*h(i,j,k)
+          enddo
+        enddo ! k
+      enddo ! j
+      call post_data(CS%id_thkcello, CS%diag_tmp3d, CS%diag)
+    endif
+  endif
+
   ! volume mean potential temperature 
   if (CS%id_thetaoga>0) then
     thetaoga = global_volume_mean(tv%T, h, G)
@@ -1027,8 +1059,12 @@ subroutine MOM_diagnostics_init(MIS, ADp, CDp, Time, G, param_file, diag, CS)
   CS%id_masso = register_scalar_field('ocean_model', 'masso', Time,  &
       diag, 'Mass of liquid ocean', 'kg', standard_name='sea_water_mass')
 
-  if ((CS%id_masscello>0) .or. (CS%id_masso>0) .and. .not.ASSOCIATED(CS%masscello)) then
-    call safe_alloc_ptr(CS%masscello,isd,ied,jsd,jed,nz)
+  CS%id_thkcello = register_diag_field('ocean_model', 'thkcello', diag%axesTL, Time, &
+      long_name = 'Cell Thickness', standard_name='cell_thickness', units='m')
+
+  if (((CS%id_masscello>0) .or. (CS%id_masso>0) .or. (CS%id_thkcello>0.and..not.G%Boussinesq)) &
+      .and. .not.ASSOCIATED(CS%diag_tmp3d)) then
+    call safe_alloc_ptr(CS%diag_tmp3d,isd,ied,jsd,jed,nz)
   endif
 
   CS%id_thetaoga = register_scalar_field('ocean_model', 'thetaoga',    &
@@ -1281,7 +1317,7 @@ subroutine MOM_diagnostics_end(CS, ADp)
   if (ASSOCIATED(CS%vh_Rlay))    deallocate(CS%vh_Rlay)
   if (ASSOCIATED(CS%uhGM_Rlay))  deallocate(CS%uhGM_Rlay)
   if (ASSOCIATED(CS%vhGM_Rlay))  deallocate(CS%vhGM_Rlay)
-  if (ASSOCIATED(CS%masscello))  deallocate(CS%masscello)
+  if (ASSOCIATED(CS%diag_tmp3d)) deallocate(CS%diag_tmp3d)
 
   if (ASSOCIATED(ADp%gradKEu))    deallocate(ADp%gradKEu)
   if (ASSOCIATED(ADp%gradKEu))    deallocate(ADp%gradKEu)