diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index c01569105f..712e945e42 100644 --- 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 index 20f8eb6da7..7a4a814542 100644 --- 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)