NOAA-EMC · jiandewang · Apr 6, 2020 · Sep 4, 2019 · Sep 9, 2019 · Sep 9, 2019
diff --git a/config_src/mct_driver/mom_ocean_model_mct.F90 b/config_src/mct_driver/mom_ocean_model_mct.F90
@@ -80,16 +80,8 @@ module MOM_ocean_model_mct
 public ocean_model_restart
 public ice_ocn_bnd_type_chksum
 public ocean_public_type_chksum
-public ocean_model_data_get
 public get_ocean_grid
 
-!> This interface extracts a named scalar field or array from the ocean surface or public type
-interface ocean_model_data_get
-  module procedure ocean_model_data1D_get
-  module procedure ocean_model_data2D_get
-end interface
-
-
 !> This type is used for communication with other components via the FMS coupler.
 !! The element names and types can be changed only with great deliberation, hence
 !! the persistnce of things like the cutsy element name "avg_kount".
@@ -1052,98 +1044,6 @@ subroutine Ocean_stock_pe(OS, index, value, time_index)
 
 end subroutine Ocean_stock_pe
 
-!> This subroutine extracts a named 2-D field from the ocean surface or public type
-subroutine ocean_model_data2D_get(OS,Ocean, name, array2D,isc,jsc)
-  use MOM_constants, only : CELSIUS_KELVIN_OFFSET
-  type(ocean_state_type),     pointer    :: OS    !< A pointer to the structure containing the
-                                                  !! internal ocean state (intent in).
-  type(ocean_public_type),    intent(in) :: Ocean !< A structure containing various publicly
-                                                  !! visible ocean surface fields.
-  character(len=*)          , intent(in) :: name  !< The name of the field to extract
-  real, dimension(isc:,jsc:), intent(out):: array2D !< The values of the named field, it must
-                                                  !! cover only the computational domain
-  integer                   , intent(in) :: isc   !< The starting i-index of array2D
-  integer                   , intent(in) :: jsc   !< The starting j-index of array2D
-
-  integer :: g_isc, g_iec, g_jsc, g_jec,g_isd, g_ied, g_jsd, g_jed, i, j
-
-  if (.not.associated(OS)) return
-  if (.not.OS%is_ocean_pe) return
-
-! The problem is %areaT is on MOM domain but Ice_Ocean_Boundary%... is on mpp domain.
-! We want to return the MOM data on the mpp (compute) domain
-! Get MOM domain extents
-  call mpp_get_compute_domain(OS%grid%Domain%mpp_domain, g_isc, g_iec, g_jsc, g_jec)
-  call mpp_get_data_domain   (OS%grid%Domain%mpp_domain, g_isd, g_ied, g_jsd, g_jed)
-
-  g_isc = g_isc-g_isd+1 ; g_iec = g_iec-g_isd+1 ; g_jsc = g_jsc-g_jsd+1 ; g_jec = g_jec-g_jsd+1
-
-
-  select case(name)
-  case('area')
-    array2D(isc:,jsc:) = OS%US%L_to_m**2*OS%grid%areaT(g_isc:g_iec,g_jsc:g_jec)
-  case('mask')
-    array2D(isc:,jsc:) = OS%grid%mask2dT(g_isc:g_iec,g_jsc:g_jec)
-!OR same result
-!     do j=g_jsc,g_jec ; do i=g_isc,g_iec
-!        array2D(isc+i-g_isc,jsc+j-g_jsc) = OS%grid%mask2dT(i,j)
-!     enddo ; enddo
-  case('t_surf')
-    array2D(isc:,jsc:) = Ocean%t_surf(isc:,jsc:)-CELSIUS_KELVIN_OFFSET
-  case('t_pme')
-    array2D(isc:,jsc:) = Ocean%t_surf(isc:,jsc:)-CELSIUS_KELVIN_OFFSET
-  case('t_runoff')
-    array2D(isc:,jsc:) = Ocean%t_surf(isc:,jsc:)-CELSIUS_KELVIN_OFFSET
-  case('t_calving')
-    array2D(isc:,jsc:) = Ocean%t_surf(isc:,jsc:)-CELSIUS_KELVIN_OFFSET
-  case('btfHeat')
-    array2D(isc:,jsc:) = 0
-  case('tlat')
-    array2D(isc:,jsc:) = OS%grid%geoLatT(g_isc:g_iec,g_jsc:g_jec)
-  case('tlon')
-    array2D(isc:,jsc:) = OS%grid%geoLonT(g_isc:g_iec,g_jsc:g_jec)
-  case('ulat')
-    array2D(isc:,jsc:) = OS%grid%geoLatCu(g_isc:g_iec,g_jsc:g_jec)
-  case('ulon')
-    array2D(isc:,jsc:) = OS%grid%geoLonCu(g_isc:g_iec,g_jsc:g_jec)
-  case('vlat')
-    array2D(isc:,jsc:) = OS%grid%geoLatCv(g_isc:g_iec,g_jsc:g_jec)
-  case('vlon')
-    array2D(isc:,jsc:) = OS%grid%geoLonCv(g_isc:g_iec,g_jsc:g_jec)
-  case('geoLatBu')
-    array2D(isc:,jsc:) = OS%grid%geoLatBu(g_isc:g_iec,g_jsc:g_jec)
-  case('geoLonBu')
-    array2D(isc:,jsc:) = OS%grid%geoLonBu(g_isc:g_iec,g_jsc:g_jec)
-  case('cos_rot')
-    array2D(isc:,jsc:) = OS%grid%cos_rot(g_isc:g_iec,g_jsc:g_jec) ! =1
-  case('sin_rot')
-    array2D(isc:,jsc:) = OS%grid%sin_rot(g_isc:g_iec,g_jsc:g_jec) ! =0
-  case default
-    call MOM_error(FATAL,'get_ocean_grid_data2D: unknown argument name='//name)
-  end select
-end subroutine ocean_model_data2D_get
-
-!> This subroutine extracts a named scalar field from the ocean surface or public type
-subroutine ocean_model_data1D_get(OS, Ocean, name, value)
-  type(ocean_state_type),     pointer    :: OS    !< A pointer to the structure containing the
-                                                  !! internal ocean state (intent in).
-  type(ocean_public_type),    intent(in) :: Ocean !< A structure containing various publicly
-                                                  !! visible ocean surface fields.
-  character(len=*)          , intent(in) :: name  !< The name of the field to extract
-  real                      , intent(out):: value !< The value of the named field
-
-  if (.not.associated(OS)) return
-  if (.not.OS%is_ocean_pe) return
-
-  select case(name)
-  case('c_p')
-     value = OS%C_p
-  case default
-     call MOM_error(FATAL,'get_ocean_grid_data1D: unknown argument name='//name)
-  end select
-
-end subroutine ocean_model_data1D_get
-
 !> Write out FMS-format checsums on fields from the ocean surface state
 subroutine ocean_public_type_chksum(id, timestep, ocn)
 

diff --git a/config_src/mct_driver/ocn_comp_mct.F90 b/config_src/mct_driver/ocn_comp_mct.F90
@@ -36,7 +36,7 @@ module ocn_comp_mct
 use MOM_time_manager,     only: time_type, set_date, set_time, set_calendar_type, NOLEAP
 use MOM_time_manager,     only: operator(+), operator(-), operator(*), operator(/)
 use MOM_time_manager,     only: operator(==), operator(/=), operator(>), get_time
-use MOM_file_parser,      only: get_param, log_version, param_file_type
+use MOM_file_parser,      only: get_param, log_version, param_file_type, close_param_file
 use MOM_get_input,        only: Get_MOM_Input, directories
 use MOM_EOS,              only: gsw_sp_from_sr, gsw_pt_from_ct
 use MOM_constants,        only: CELSIUS_KELVIN_OFFSET
@@ -281,6 +281,9 @@ subroutine ocn_init_mct( EClock, cdata_o, x2o_o, o2x_o, NLFilename )
     glb%c1 = 0.0; glb%c2 = 0.0; glb%c3 = 0.0; glb%c4 = 0.0
   endif
 
+  ! Close param file before it gets opened by ocean_model_init again.
+  call close_param_file(param_file)
+
   ! Initialize the MOM6 model
   runtype = get_runtype()
   if (runtype == "initial") then

diff --git a/config_src/nuopc_driver/mom_cap.F90 b/config_src/nuopc_driver/mom_cap.F90
@@ -763,11 +763,6 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc)
   end if
 
   if (cesm_coupled) then
-  !TODO: check if still needed
-     if (len_trim(scalar_field_name) > 0) then
-        call fld_list_add(fldsToOcn_num, fldsToOcn, trim(scalar_field_name), "will_provide")
-        call fld_list_add(fldsFrOcn_num, fldsFrOcn, trim(scalar_field_name), "will_provide")
-     endif
     !call fld_list_add(fldsToOcn_num, fldsToOcn, "Sw_lamult"                 , "will provide")
     !call fld_list_add(fldsToOcn_num, fldsToOcn, "Sw_ustokes"                , "will provide")
     !call fld_list_add(fldsToOcn_num, fldsToOcn, "Sw_vstokes"                , "will provide")

diff --git a/config_src/nuopc_driver/mom_ocean_model_nuopc.F90 b/config_src/nuopc_driver/mom_ocean_model_nuopc.F90
@@ -77,17 +77,9 @@ module MOM_ocean_model_nuopc
 public ocean_model_restart
 public ice_ocn_bnd_type_chksum
 public ocean_public_type_chksum
-public ocean_model_data_get
 public get_ocean_grid
 public get_eps_omesh
 
-!> This interface extracts a named scalar field or array from the ocean surface or public type
-interface ocean_model_data_get
-  module procedure ocean_model_data1D_get
-  module procedure ocean_model_data2D_get
-end interface
-
-
 !> This type is used for communication with other components via the FMS coupler.
 !! The element names and types can be changed only with great deliberation, hence
 !! the persistnce of things like the cutsy element name "avg_kount".
@@ -1047,99 +1039,6 @@ subroutine Ocean_stock_pe(OS, index, value, time_index)
 
 end subroutine Ocean_stock_pe
 
-!> This subroutine extracts a named 2-D field from the ocean surface or public type
-subroutine ocean_model_data2D_get(OS, Ocean, name, array2D, isc, jsc)
-  use MOM_constants, only : CELSIUS_KELVIN_OFFSET
-  type(ocean_state_type),     pointer    :: OS    !< A pointer to the structure containing the
-                                                  !! internal ocean state (intent in).
-  type(ocean_public_type),    intent(in) :: Ocean !< A structure containing various publicly
-                                                  !! visible ocean surface fields.
-  character(len=*)          , intent(in) :: name  !< The name of the field to extract
-  real, dimension(isc:,jsc:), intent(out):: array2D !< The values of the named field, it must
-                                                  !! cover only the computational domain
-  integer                   , intent(in) :: isc   !< The starting i-index of array2D
-  integer                   , intent(in) :: jsc   !< The starting j-index of array2D
-
-  integer :: g_isc, g_iec, g_jsc, g_jec,g_isd, g_ied, g_jsd, g_jed, i, j
-
-  if (.not.associated(OS)) return
-  if (.not.OS%is_ocean_pe) return
-
-! The problem is %areaT is on MOM domain but Ice_Ocean_Boundary%... is on mpp domain.
-! We want to return the MOM data on the mpp (compute) domain
-! Get MOM domain extents
-  call mpp_get_compute_domain(OS%grid%Domain%mpp_domain, g_isc, g_iec, g_jsc, g_jec)
-  call mpp_get_data_domain   (OS%grid%Domain%mpp_domain, g_isd, g_ied, g_jsd, g_jed)
-
-  g_isc = g_isc-g_isd+1 ; g_iec = g_iec-g_isd+1 ; g_jsc = g_jsc-g_jsd+1 ; g_jec = g_jec-g_jsd+1
-
-
-  select case(name)
-  case('area')
-     array2D(isc:,jsc:) = OS%US%L_to_m**2*OS%grid%areaT(g_isc:g_iec,g_jsc:g_jec)
-  case('mask')
-     array2D(isc:,jsc:) = OS%grid%mask2dT(g_isc:g_iec,g_jsc:g_jec)
-!OR same result
-!     do j=g_jsc,g_jec ; do i=g_isc,g_iec
-!        array2D(isc+i-g_isc,jsc+j-g_jsc) = OS%grid%mask2dT(i,j)
-!     enddo ; enddo
-  case('t_surf')
-     array2D(isc:,jsc:) = Ocean%t_surf(isc:,jsc:)-CELSIUS_KELVIN_OFFSET
-  case('t_pme')
-     array2D(isc:,jsc:) = Ocean%t_surf(isc:,jsc:)-CELSIUS_KELVIN_OFFSET
-  case('t_runoff')
-     array2D(isc:,jsc:) = Ocean%t_surf(isc:,jsc:)-CELSIUS_KELVIN_OFFSET
-  case('t_calving')
-     array2D(isc:,jsc:) = Ocean%t_surf(isc:,jsc:)-CELSIUS_KELVIN_OFFSET
-  case('btfHeat')
-     array2D(isc:,jsc:) = 0
-  case('tlat')
-     array2D(isc:,jsc:) = OS%grid%geoLatT(g_isc:g_iec,g_jsc:g_jec)
-  case('tlon')
-     array2D(isc:,jsc:) = OS%grid%geoLonT(g_isc:g_iec,g_jsc:g_jec)
-  case('ulat')
-     array2D(isc:,jsc:) = OS%grid%geoLatCu(g_isc:g_iec,g_jsc:g_jec)
-  case('ulon')
-     array2D(isc:,jsc:) = OS%grid%geoLonCu(g_isc:g_iec,g_jsc:g_jec)
-  case('vlat')
-     array2D(isc:,jsc:) = OS%grid%geoLatCv(g_isc:g_iec,g_jsc:g_jec)
-  case('vlon')
-     array2D(isc:,jsc:) = OS%grid%geoLonCv(g_isc:g_iec,g_jsc:g_jec)
-  case('geoLatBu')
-     array2D(isc:,jsc:) = OS%grid%geoLatBu(g_isc:g_iec,g_jsc:g_jec)
-  case('geoLonBu')
-     array2D(isc:,jsc:) = OS%grid%geoLonBu(g_isc:g_iec,g_jsc:g_jec)
-  case('cos_rot')
-     array2D(isc:,jsc:) = OS%grid%cos_rot(g_isc:g_iec,g_jsc:g_jec) ! =1
-  case('sin_rot')
-     array2D(isc:,jsc:) = OS%grid%sin_rot(g_isc:g_iec,g_jsc:g_jec) ! =0
-  case default
-     call MOM_error(FATAL,'get_ocean_grid_data2D: unknown argument name='//name)
-  end select
-
-end subroutine ocean_model_data2D_get
-
-!> This subroutine extracts a named scalar field from the ocean surface or public type
-subroutine ocean_model_data1D_get(OS, Ocean, name, value)
-  type(ocean_state_type),     pointer    :: OS    !< A pointer to the structure containing the
-                                                  !! internal ocean state (intent in).
-  type(ocean_public_type),    intent(in) :: Ocean !< A structure containing various publicly
-                                                  !! visible ocean surface fields.
-  character(len=*)          , intent(in) :: name  !< The name of the field to extract
-  real                      , intent(out):: value !< The value of the named field
-
-  if (.not.associated(OS)) return
-  if (.not.OS%is_ocean_pe) return
-
-  select case(name)
-  case('c_p')
-     value = OS%C_p
-  case default
-     call MOM_error(FATAL,'get_ocean_grid_data1D: unknown argument name='//name)
-  end select
-
-end subroutine ocean_model_data1D_get
-
 !> Write out FMS-format checsums on fields from the ocean surface state
 subroutine ocean_public_type_chksum(id, timestep, ocn)
 

diff --git a/config_src/nuopc_driver/mom_surface_forcing_nuopc.F90 b/config_src/nuopc_driver/mom_surface_forcing_nuopc.F90
@@ -444,12 +444,12 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G,
     ! liquid runoff flux
     if (associated(IOB%lrunoff)) then
       if(CS%liquid_runoff_from_data)call data_override('OCN', 'runoff', IOB%lrunoff, Time)
-      fluxes%lrunoff(i,j) = IOB%lrunoff(i-i0,j-j0) * G%mask2dT(i,j)
+      fluxes%lrunoff(i,j) = kg_m2_s_conversion * IOB%lrunoff(i-i0,j-j0) * G%mask2dT(i,j)
     endif
 
     ! ice runoff flux
     if (associated(IOB%frunoff)) then
-      fluxes%frunoff(i,j) = IOB%frunoff(i-i0,j-j0) * G%mask2dT(i,j)
+      fluxes%frunoff(i,j) = kg_m2_s_conversion * IOB%frunoff(i-i0,j-j0) * G%mask2dT(i,j)
     endif
 
     if (associated(IOB%ustar_berg)) &
@@ -465,7 +465,7 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G,
       fluxes%heat_content_lrunoff(i,j) = IOB%lrunoff_hflx(i-i0,j-j0) * G%mask2dT(i,j)
 
     if (associated(IOB%frunoff_hflx)) &
-      fluxes%heat_content_frunoff(i,j) = IOB%frunoff_hflx(i-i0,j-j0) * G%mask2dT(i,j)
+      fluxes%heat_content_frunoff(i,j) = kg_m2_s_conversion * IOB%frunoff_hflx(i-i0,j-j0) * G%mask2dT(i,j)
 
     if (associated(IOB%lw_flux)) &
          fluxes%LW(i,j) = IOB%lw_flux(i-i0,j-j0) * G%mask2dT(i,j)

diff --git a/src/core/MOM.F90 b/src/core/MOM.F90
@@ -77,7 +77,7 @@ module MOM
 use MOM_hor_index,             only : hor_index_type, hor_index_init
 use MOM_interface_heights,     only : find_eta
 use MOM_lateral_mixing_coeffs, only : calc_slope_functions, VarMix_init
-use MOM_lateral_mixing_coeffs, only : calc_resoln_function, VarMix_CS
+use MOM_lateral_mixing_coeffs, only : calc_resoln_function, calc_depth_function, VarMix_CS
 use MOM_MEKE,                  only : MEKE_init, MEKE_alloc_register_restart, step_forward_MEKE, MEKE_CS
 use MOM_MEKE_types,            only : MEKE_type
 use MOM_mixed_layer_restrat,   only : mixedlayer_restrat, mixedlayer_restrat_init, mixedlayer_restrat_CS
@@ -566,6 +566,7 @@ subroutine step_MOM(forces, fluxes, sfc_state, Time_start, time_int_in, CS, &
     if (associated(CS%VarMix)) then
       call enable_averages(cycle_time, Time_start + real_to_time(US%T_to_s*cycle_time), CS%diag)
       call calc_resoln_function(h, CS%tv, G, GV, US, CS%VarMix)
+      call calc_depth_function(G, CS%VarMix)
       call disable_averaging(CS%diag)
     endif
   endif
@@ -1406,6 +1407,7 @@ subroutine step_offline(forces, fluxes, sfc_state, Time_start, time_interval, CS
           if (associated(CS%VarMix)) then
             call pass_var(CS%h, G%Domain)
             call calc_resoln_function(CS%h, CS%tv, G, GV, US, CS%VarMix)
+            call calc_depth_function(G, CS%VarMix)
             call calc_slope_functions(CS%h, CS%tv, dt_offline, G, GV, US, CS%VarMix)
           endif
           call tracer_hordiff(CS%h, dt_offline, CS%MEKE, CS%VarMix, G, GV, US, &
@@ -1431,6 +1433,7 @@ subroutine step_offline(forces, fluxes, sfc_state, Time_start, time_interval, CS
           if (associated(CS%VarMix)) then
             call pass_var(CS%h, G%Domain)
             call calc_resoln_function(CS%h, CS%tv, G, GV, US, CS%VarMix)
+            call calc_depth_function(G, CS%VarMix)
             call calc_slope_functions(CS%h, CS%tv, dt_offline, G, GV, US, CS%VarMix)
           endif
           call tracer_hordiff(CS%h, dt_offline, CS%MEKE, CS%VarMix, G, GV, US, &
@@ -2373,7 +2376,7 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, &
   endif
 
   call tracer_advect_init(Time, G, US, param_file, diag, CS%tracer_adv_CSp)
-  call tracer_hor_diff_init(Time, G, US, param_file, diag, CS%tv%eqn_of_state, &
+  call tracer_hor_diff_init(Time, G, US, param_file, diag, CS%tv%eqn_of_state, CS%diabatic_CSp, &
                             CS%tracer_diff_CSp)
 
   call lock_tracer_registry(CS%tracer_Reg)
@@ -2894,7 +2897,7 @@ subroutine extract_surface_state(CS, sfc_state)
 
 
   if (allocated(sfc_state%melt_potential)) then
-  !$OMP parallel do default(shared)
+  !$OMP parallel do default(shared) private(depth_ml, dh, T_freeze, depth, delT)
     do j=js,je
       do i=is,ie
         depth(i) = 0.0

diff --git a/src/core/MOM_dynamics_split_RK2.F90 b/src/core/MOM_dynamics_split_RK2.F90
@@ -681,7 +681,7 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, &
   call cpu_clock_begin(id_clock_horvisc)
   call horizontal_viscosity(u_av, v_av, h_av, CS%diffu, CS%diffv, &
                             MEKE, Varmix, G, GV, US, CS%hor_visc_CSp, &
-                            OBC=CS%OBC, BT=CS%barotropic_CSp)
+                            OBC=CS%OBC, BT=CS%barotropic_CSp, TD=thickness_diffuse_CSp)
   call cpu_clock_end(id_clock_horvisc)
   if (showCallTree) call callTree_wayPoint("done with horizontal_viscosity (step_MOM_dyn_split_RK2)")
 
@@ -1149,7 +1149,8 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param
       .not. query_initialized(CS%diffv,"diffv",restart_CS)) then
     call horizontal_viscosity(u, v, h, CS%diffu, CS%diffv, MEKE, VarMix, &
                               G, GV, US, CS%hor_visc_CSp, &
-                              OBC=CS%OBC, BT=CS%barotropic_CSp)
+                              OBC=CS%OBC, BT=CS%barotropic_CSp, &
+                              TD=thickness_diffuse_CSp)
   else
     if ( (US%s_to_T_restart * US%m_to_L_restart /= 0.0) .and. &
          (US%m_to_L * US%s_to_T_restart**2 /= US%m_to_L_restart * US%s_to_T**2) ) then

diff --git a/src/core/MOM_unit_tests.F90 b/src/core/MOM_unit_tests.F90
@@ -5,10 +5,11 @@ module MOM_unit_tests
 
 use MOM_error_handler,     only : MOM_error, FATAL, is_root_pe
 
-use MOM_string_functions,  only : string_functions_unit_tests
-use MOM_remapping,         only : remapping_unit_tests
-use MOM_neutral_diffusion, only : neutral_diffusion_unit_tests
-use MOM_diag_vkernels,     only : diag_vkernels_unit_tests
+use MOM_string_functions,          only : string_functions_unit_tests
+use MOM_remapping,                 only : remapping_unit_tests
+use MOM_neutral_diffusion,         only : neutral_diffusion_unit_tests
+use MOM_diag_vkernels,             only : diag_vkernels_unit_tests
+use MOM_lateral_boundary_diffusion,   only : near_boundary_unit_tests
 
 implicit none ; private
 
@@ -35,6 +36,9 @@ subroutine unit_tests(verbosity)
        "MOM_unit_tests: neutralDiffusionUnitTests FAILED")
     if (diag_vkernels_unit_tests(verbose)) call MOM_error(FATAL, &
        "MOM_unit_tests: diag_vkernels_unit_tests FAILED")
+    if (near_boundary_unit_tests(verbose)) call MOM_error(FATAL, &
+       "MOM_unit_tests: near_boundary_unit_tests FAILED")
+
   endif
 
 end subroutine unit_tests

diff --git a/src/diagnostics/MOM_diagnostics.F90 b/src/diagnostics/MOM_diagnostics.F90
@@ -600,7 +600,7 @@ subroutine calculate_diagnostic_fields(u, v, h, uh, vh, tv, ADp, CDp, p_surf, &
       if (CS%id_rhopot2 > 0) call post_data(CS%id_rhopot2, Rcv, CS%diag)
     endif
     if (CS%id_rhoinsitu > 0) then
-!$OMP parallel do default(none) shared(tv,Rcv,is,ie,js,je,nz,pressure_1d,h,GV)
+!$OMP parallel do default(none) shared(tv,Rcv,is,ie,js,je,nz,h,GV) private(pressure_1d)
       do j=js,je
         pressure_1d(:) = 0. ! Start at p=0 Pa at surface
         do k=1,nz