Skip to content

Commit

Permalink
Merge pull request #74 from NCAR/dev/nuopc
Browse files Browse the repository at this point in the history 
updates to mct_driver and nuopc_driver
  • Loading branch information
@alperaltuntas
alperaltuntas authored Jul 24, 2018
2 parents bd90705 + 2587399 commit 04845d4
Show file tree
Hide file tree
Showing 11 changed files with 9,093 additions and 2,010 deletions.
1,037 changes: 1,037 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.

239 changes: 239 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,239 @@
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)
end if
end do
end do

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)
end do
end do
end if

end subroutine ocn_import

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

!> Maps outgoing ocean data to MCT attribute vector real array
subroutine ocn_export(ind, ocn_public, grid, o2x)
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

! 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
!-----------------------------------------------------------------------

! 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)
! 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)
end do
end do

! 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
end if
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
end do; end do

! 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
end if
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
end do; end do

end subroutine ocn_export

end module ocn_cap_methods
Loading

0 comments on commit 04845d4

Please sign in to comment.