Skip to content

Commit

Permalink
Merge remote-tracking branch 'GFDL_MOM6/dev/gfdl' into dev/gfdl
Browse files Browse the repository at this point in the history 
merge in latest MOM updates
* GFDL_MOM6/dev/gfdl: (220 commits)
  Update cache dir for /lustre/f2
  Update cache dir for /lustre/f2
  Calculate height-related diagnostics in Z
  Set coord_adapt and coord_slight parameters in H
  Corrected comments in build_zlike_column
  +Pass max_depth to initialize_regridding in Z
  Corrected comments in build_sigma_column
  Combined scaling factors in build_adapt_column
  Clarified comments in Idealized_Hurricane
  +Recast MOM_diag_to_Z to work in units of Z
  +Recast MOM_ALE_sponge to work in units of Z
  +Added m_to_Z arg to horiz_interp_and_extrap_tracer
  Find energetic_PBL column height changes in Z
  +Find diapyc_energy_req column height changes in Z
  Use local variables to rescale in MOM_Point_Accel
  Recast MOM_sum_output to work in units of Z
  Simplify MOM_diagnostics code
  +Rescale depth inside of MEKE_lengthScales_0d
  +Add conversion argument to register_static_field
  Rescale values reported by PointAccel
  ...
  • Loading branch information
Jessica.Liptak committed Nov 7, 2018
2 parents 975262d + 3b04fb2 commit b7f50fc
Show file tree
Hide file tree
Showing 124 changed files with 9,566 additions and 7,909 deletions.
2 changes: 1 addition & 1 deletion .gitlab-ci.yml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,7 +8,7 @@ stages:
# Merges MOM6 with dev/gfdl. Changes directory to test directory, if it exists.
before_script:
- MOM6_SRC=$CI_PROJECT_DIR
- CACHE_DIR=/lustre/f1/oar.gfdl.ogrp-account/runner/cache/
- echo Cache directory set to ${CACHE_DIR:=/lustre/f2/scratch/oar.gfdl.ogrp-account/runner/cache/}
- git pull --no-edit https://github.com/NOAA-GFDL/MOM6.git dev/gfdl && git submodule init && git submodule update
- pwd ; ls

Expand Down
22 changes: 20 additions & 2 deletions config_src/coupled_driver/ocean_model_MOM.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -237,6 +237,12 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn)
! Local variables
real :: Rho0 ! The Boussinesq ocean density, in kg m-3.
real :: G_Earth ! The gravitational acceleration in m s-2.
real :: HFrz !< If HFrz > 0 (m), melt potential will be computed.
!! The actual depth over which melt potential is computed will
!! min(HFrz, OBLD), where OBLD is the boundary layer depth.
!! If HFrz <= 0 (default), melt potential will not be computed.
logical :: use_melt_pot!< If true, allocate melt_potential array

! This include declares and sets the variable "version".
#include "version_variable.h"
character(len=40) :: mdl = "ocean_model_init" ! This module's name.
Expand Down Expand Up @@ -330,8 +336,20 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn)

! Consider using a run-time flag to determine whether to do the diagnostic
! vertical integrals, since the related 3-d sums are not negligible in cost.
call allocate_surface_state(OS%sfc_state, OS%grid, use_temperature, &
do_integrals=.true., gas_fields_ocn=gas_fields_ocn)
call get_param(param_file, mdl, "HFREEZE", HFrz, &
"If HFREEZE > 0, melt potential will be computed. The actual depth \n"//&
"over which melt potential is computed will be min(HFREEZE, OBLD), \n"//&
"where OBLD is the boundary layer depth. If HFREEZE <= 0 (default), \n"//&
"melt potential will not be computed.", units="m", default=-1.0, do_not_log=.true.)

if (HFrz .gt. 0.0) then
use_melt_pot=.true.
else
use_melt_pot=.false.
endif

call allocate_surface_state(OS%sfc_state, OS%grid, use_temperature, do_integrals=.true., &
gas_fields_ocn=gas_fields_ocn, use_meltpot=use_melt_pot)

call surface_forcing_init(Time_in, OS%grid, param_file, OS%diag, &
OS%forcing_CSp)
Expand Down
1,067 changes: 1,067 additions & 0 deletions config_src/mct_driver/MOM_ocean_model.F90
View file

Large diffs are not rendered by default.

1,338 changes: 1,338 additions & 0 deletions config_src/mct_driver/MOM_surface_forcing.F90
View file

Large diffs are not rendered by default.

258 changes: 258 additions & 0 deletions config_src/mct_driver/ocn_cap_methods.F90
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,258 @@
module ocn_cap_methods

use ESMF, only: ESMF_clock, ESMF_time, ESMF_ClockGet, ESMF_TimeGet
use MOM_ocean_model, only: ocean_public_type, ocean_state_type
use MOM_surface_forcing, only: ice_ocean_boundary_type
use MOM_grid, only: ocean_grid_type
use MOM_domains, only: pass_var
use MOM_error_handler, only: is_root_pe
use mpp_domains_mod, only: mpp_get_compute_domain
use ocn_cpl_indices, only: cpl_indices_type

implicit none
private

public :: ocn_import
public :: ocn_export

logical, parameter :: debug=.false.

!=======================================================================
contains
!=======================================================================

!> Maps incomping ocean data to MOM6 data structures
subroutine ocn_import(x2o, ind, grid, ice_ocean_boundary, ocean_public, logunit, Eclock, c1, c2, c3, c4)
real(kind=8) , intent(in) :: x2o(:,:) !< incoming data
type(cpl_indices_type) , intent(in) :: ind !< Structure with MCT attribute vects and indices
type(ocean_grid_type) , intent(in) :: grid !< Ocean model grid
type(ice_ocean_boundary_type) , intent(inout) :: ice_ocean_boundary !< Ocean boundary forcing
type(ocean_public_type) , intent(in) :: ocean_public !< Ocean surface state
integer , intent(in) :: logunit !< Unit for stdout output
type(ESMF_Clock) , intent(in) :: EClock !< Time and time step ? \todo Why must this
real(kind=8), optional , intent(in) :: c1, c2, c3, c4 !< Coeffs. used in the shortwave decomposition

! Local variables
integer :: i, j, ig, jg, isc, iec, jsc, jec ! Grid indices
integer :: k
integer :: day, secs, rc
type(ESMF_time) :: currTime
character(*), parameter :: F01 = "('(ocn_import) ',a,4(i6,2x),d21.14)"
!-----------------------------------------------------------------------

isc = GRID%isc; iec = GRID%iec ; jsc = GRID%jsc; jec = GRID%jec

k = 0
do j = jsc, jec
jg = j + grid%jsc - jsc
do i = isc, iec
ig = i + grid%jsc - isc
k = k + 1 ! Increment position within gindex

! taux
ice_ocean_boundary%u_flux(i,j) = x2o(ind%x2o_Foxx_taux,k)

! tauy
ice_ocean_boundary%v_flux(i,j) = x2o(ind%x2o_Foxx_tauy,k)

! liquid precipitation (rain)
ice_ocean_boundary%lprec(i,j) = x2o(ind%x2o_Faxa_rain,k)

! frozen precipitation (snow)
ice_ocean_boundary%fprec(i,j) = x2o(ind%x2o_Faxa_snow,k)

! longwave radiation, sum up and down (W/m2)
ice_ocean_boundary%lw_flux(i,j) = (x2o(ind%x2o_Faxa_lwdn,k) + x2o(ind%x2o_Foxx_lwup,k))

! specific humitidy flux
ice_ocean_boundary%q_flux(i,j) = x2o(ind%x2o_Foxx_evap,k) !???TODO: should this be a minus sign

! sensible heat flux (W/m2)
ice_ocean_boundary%t_flux(i,j) = x2o(ind%x2o_Foxx_sen,k) !???TODO: should this be a minus sign

! latent heat flux (W/m^2)
ice_ocean_boundary%latent_flux(i,j) = x2o(ind%x2o_Foxx_lat,k) !???TODO: should this be a minus sign

! liquid runoff
ice_ocean_boundary%rofl_flux(i,j) = x2o(ind%x2o_Foxx_rofl,k) * GRID%mask2dT(ig,jg)

! ice runoff
ice_ocean_boundary%rofi_flux(i,j) = x2o(ind%x2o_Foxx_rofi,k) * GRID%mask2dT(ig,jg)

! surface pressure
ice_ocean_boundary%p(i,j) = x2o(ind%x2o_Sa_pslv,k) * GRID%mask2dT(ig,jg)

! salt flux
ice_ocean_boundary%salt_flux(i,j) = x2o(ind%x2o_Fioi_salt,k) * GRID%mask2dT(ig,jg)

! 1) visible, direct shortwave (W/m2)
! 2) visible, diffuse shortwave (W/m2)
! 3) near-IR, direct shortwave (W/m2)
! 4) near-IR, diffuse shortwave (W/m2)
if (present(c1) .and. present(c2) .and. present(c3) .and. present(c4)) then
! Use runtime coefficients to decompose net short-wave heat flux into 4 components
ice_ocean_boundary%sw_flux_vis_dir(i,j) = x2o(ind%x2o_Foxx_swnet,k) * c1 * GRID%mask2dT(ig,jg)
ice_ocean_boundary%sw_flux_vis_dif(i,j) = x2o(ind%x2o_Foxx_swnet,k) * c2 * GRID%mask2dT(ig,jg)
ice_ocean_boundary%sw_flux_nir_dir(i,j) = x2o(ind%x2o_Foxx_swnet,k) * c3 * GRID%mask2dT(ig,jg)
ice_ocean_boundary%sw_flux_nir_dif(i,j) = x2o(ind%x2o_Foxx_swnet,k) * c4 * GRID%mask2dT(ig,jg)
else
ice_ocean_boundary%sw_flux_vis_dir(i,j) = x2o(ind%x2o_Faxa_swvdr,k) * GRID%mask2dT(ig,jg)
ice_ocean_boundary%sw_flux_vis_dif(i,j) = x2o(ind%x2o_Faxa_swvdf,k) * GRID%mask2dT(ig,jg)
ice_ocean_boundary%sw_flux_nir_dir(i,j) = x2o(ind%x2o_Faxa_swndr,k) * GRID%mask2dT(ig,jg)
ice_ocean_boundary%sw_flux_nir_dif(i,j) = x2o(ind%x2o_Faxa_swndf,k) * GRID%mask2dT(ig,jg)
endif
enddo
enddo

if (debug .and. is_root_pe()) then
call ESMF_ClockGet(EClock, CurrTime=CurrTime, rc=rc)
call ESMF_TimeGet(CurrTime, d=day, s=secs, rc=rc)

do j = GRID%jsc, GRID%jec
do i = GRID%isc, GRID%iec
write(logunit,F01)'import: day, secs, j, i, u_flux = ',day,secs,j,i,ice_ocean_boundary%u_flux(i,j)
write(logunit,F01)'import: day, secs, j, i, v_flux = ',day,secs,j,i,ice_ocean_boundary%v_flux(i,j)
write(logunit,F01)'import: day, secs, j, i, lprec = ',day,secs,j,i,ice_ocean_boundary%lprec(i,j)
write(logunit,F01)'import: day, secs, j, i, lwrad = ',day,secs,j,i,ice_ocean_boundary%lw_flux(i,j)
write(logunit,F01)'import: day, secs, j, i, q_flux = ',day,secs,j,i,ice_ocean_boundary%q_flux(i,j)
write(logunit,F01)'import: day, secs, j, i, t_flux = ',day,secs,j,i,ice_ocean_boundary%t_flux(i,j)
write(logunit,F01)'import: day, secs, j, i, latent_flux = ',&
day,secs,j,i,ice_ocean_boundary%latent_flux(i,j)
write(logunit,F01)'import: day, secs, j, i, runoff = ',&
day,secs,j,i,ice_ocean_boundary%rofl_flux(i,j) + ice_ocean_boundary%rofi_flux(i,j)
write(logunit,F01)'import: day, secs, j, i, psurf = ',&
day,secs,j,i,ice_ocean_boundary%p(i,j)
write(logunit,F01)'import: day, secs, j, i, salt_flux = ',&
day,secs,j,i,ice_ocean_boundary%salt_flux(i,j)
write(logunit,F01)'import: day, secs, j, i, sw_flux_vis_dir = ',&
day,secs,j,i,ice_ocean_boundary%sw_flux_vis_dir(i,j)
write(logunit,F01)'import: day, secs, j, i, sw_flux_vis_dif = ',&
day,secs,j,i,ice_ocean_boundary%sw_flux_vis_dif(i,j)
write(logunit,F01)'import: day, secs, j, i, sw_flux_nir_dir = ',&
day,secs,j,i,ice_ocean_boundary%sw_flux_nir_dir(i,j)
write(logunit,F01)'import: day, secs, j, i, sw_flux_nir_dif = ',&
day,secs,j,i,ice_ocean_boundary%sw_flux_nir_dir(i,j)
enddo
enddo
endif

end subroutine ocn_import

!=======================================================================

!> Maps outgoing ocean data to MCT attribute vector real array
subroutine ocn_export(ind, ocn_public, grid, o2x, dt_int, ncouple_per_day)
type(cpl_indices_type), intent(inout) :: ind !< Structure with coupler indices and vectors
type(ocean_public_type), intent(in) :: ocn_public !< Ocean surface state
type(ocean_grid_type), intent(in) :: grid !< Ocean model grid
real(kind=8), intent(inout) :: o2x(:,:) !< MCT outgoing bugger
real(kind=8), intent(in) :: dt_int !< Amount of time over which to advance the
!! ocean (ocean_coupling_time_step), in sec
integer, intent(in) :: ncouple_per_day !< Number of ocean coupling calls per day

! Local variables
real, dimension(grid%isd:grid%ied,grid%jsd:grid%jed) :: ssh !< Local copy of sea_lev with updated halo
integer :: i, j, n, ig, jg !< Grid indices
real :: slp_L, slp_R, slp_C, slope, u_min, u_max
real :: I_time_int !< The inverse of coupling time interval in s-1.

!-----------------------------------------------------------------------

! Use Adcroft's rule of reciprocals; it does the right thing here.
I_time_int = 0.0 ; if (dt_int > 0.0) I_time_int = 1.0 / dt_int

! Copy from ocn_public to o2x. ocn_public uses global indexing with no halos.
! The mask comes from "grid" that uses the usual MOM domain that has halos
! and does not use global indexing.

n = 0
do j=grid%jsc, grid%jec
jg = j + grid%jdg_offset
do i=grid%isc,grid%iec
n = n+1
ig = i + grid%idg_offset
! surface temperature in Kelvin
o2x(ind%o2x_So_t, n) = ocn_public%t_surf(ig,jg) * grid%mask2dT(i,j)
o2x(ind%o2x_So_s, n) = ocn_public%s_surf(ig,jg) * grid%mask2dT(i,j)
o2x(ind%o2x_So_u, n) = ocn_public%u_surf(ig,jg) * grid%mask2dT(i,j)
o2x(ind%o2x_So_v, n) = ocn_public%v_surf(ig,jg) * grid%mask2dT(i,j)
o2x(ind%o2x_So_bldepth, n) = ocn_public%OBLD(ig,jg) * grid%mask2dT(i,j)
! ocean melt and freeze potential (o2x_Fioo_q), W m-2
if (ocn_public%frazil(ig,jg) > 0.0) then
! Frazil: change from J/m^2 to W/m^2
o2x(ind%o2x_Fioo_q, n) = ocn_public%frazil(ig,jg) * grid%mask2dT(i,j) * I_time_int
else
! Melt_potential: change from J/m^2 to W/m^2
o2x(ind%o2x_Fioo_q, n) = -ocn_public%melt_potential(ig,jg) * grid%mask2dT(i,j) * I_time_int !* ncouple_per_day
! make sure Melt_potential is always <= 0
if (o2x(ind%o2x_Fioo_q, n) > 0.0) o2x(ind%o2x_Fioo_q, n) = 0.0
endif
! Make a copy of ssh in order to do a halo update. We use the usual MOM domain
! in order to update halos. i.e. does not use global indexing.
ssh(i,j) = ocn_public%sea_lev(ig,jg)
enddo
enddo

! Update halo of ssh so we can calculate gradients
call pass_var(ssh, grid%domain)

! d/dx ssh
n = 0
do j=grid%jsc, grid%jec ; do i=grid%isc,grid%iec
n = n+1
! This is a simple second-order difference
! o2x(ind%o2x_So_dhdx, n) = 0.5 * (ssh(i+1,j) - ssh(i-1,j)) * grid%IdxT(i,j) * grid%mask2dT(i,j)
! This is a PLM slope which might be less prone to the A-grid null mode
slp_L = (ssh(I,j) - ssh(I-1,j)) * grid%mask2dCu(I-1,j)
if (grid%mask2dCu(I-1,j)==0.) slp_L = 0.
slp_R = (ssh(I+1,j) - ssh(I,j)) * grid%mask2dCu(I,j)
if (grid%mask2dCu(I+1,j)==0.) slp_R = 0.
slp_C = 0.5 * (slp_L + slp_R)
if ( (slp_L * slp_R) > 0.0 ) then
! This limits the slope so that the edge values are bounded by the
! two cell averages spanning the edge.
u_min = min( ssh(i-1,j), ssh(i,j), ssh(i+1,j) )
u_max = max( ssh(i-1,j), ssh(i,j), ssh(i+1,j) )
slope = sign( min( abs(slp_C), 2.*min( ssh(i,j) - u_min, u_max - ssh(i,j) ) ), slp_C )
else
! Extrema in the mean values require a PCM reconstruction avoid generating
! larger extreme values.
slope = 0.0
endif
o2x(ind%o2x_So_dhdx, n) = slope * grid%IdxT(i,j) * grid%mask2dT(i,j)
if (grid%mask2dT(i,j)==0.) o2x(ind%o2x_So_dhdx, n) = 0.0
enddo; enddo

! d/dy ssh
n = 0
do j=grid%jsc, grid%jec ; do i=grid%isc,grid%iec
n = n+1
! This is a simple second-order difference
! o2x(ind%o2x_So_dhdy, n) = 0.5 * (ssh(i,j+1) - ssh(i,j-1)) * grid%IdyT(i,j) * grid%mask2dT(i,j)
! This is a PLM slope which might be less prone to the A-grid null mode
slp_L = ssh(i,J) - ssh(i,J-1) * grid%mask2dCv(i,J-1)
if (grid%mask2dCv(i,J-1)==0.) slp_L = 0.

slp_R = ssh(i,J+1) - ssh(i,J) * grid%mask2dCv(i,J)
if (grid%mask2dCv(i,J+1)==0.) slp_R = 0.

slp_C = 0.5 * (slp_L + slp_R)
!write(6,*)'slp_L, slp_R,i,j,slp_L*slp_R', slp_L, slp_R,i,j,slp_L*slp_R
if ((slp_L * slp_R) > 0.0) then
! This limits the slope so that the edge values are bounded by the
! two cell averages spanning the edge.
u_min = min( ssh(i,j-1), ssh(i,j), ssh(i,j+1) )
u_max = max( ssh(i,j-1), ssh(i,j), ssh(i,j+1) )
slope = sign( min( abs(slp_C), 2.*min( ssh(i,j) - u_min, u_max - ssh(i,j) ) ), slp_C )
else
! Extrema in the mean values require a PCM reconstruction avoid generating
! larger extreme values.
slope = 0.0
endif
o2x(ind%o2x_So_dhdy, n) = slope * grid%IdyT(i,j) * grid%mask2dT(i,j)
if (grid%mask2dT(i,j)==0.) o2x(ind%o2x_So_dhdy, n) = 0.0
enddo; enddo

end subroutine ocn_export

end module ocn_cap_methods
Loading

0 comments on commit b7f50fc

Please sign in to comment.