Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

62 z ale remapping #188

Merged
merged 12 commits into from
Jun 15, 2015
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
68 changes: 44 additions & 24 deletions src/ALE/MOM_regridding.F90
Original file line number Diff line number Diff line change
Expand Up @@ -98,6 +98,7 @@ module MOM_regridding
public getCoordinateUnits
public getCoordinateShortName
public getStaticThickness
public buildGridZstarColumn

public DEFAULT_COORDINATE_MODE
character(len=158), parameter, public :: regriddingCoordinateModeDoc = &
Expand Down Expand Up @@ -143,7 +144,7 @@ module MOM_regridding
! Deviation tolerance between succesive grids in regridding iterations
real, parameter :: DEVIATION_TOLERANCE = 1e-10
! Maximum number of Newton-Raphson iterations. Newton-Raphson iterations are
! used to build the new grid by finding the coordinates associated with
! used to build the new grid by finding the coordinates associated with
! target densities and interpolations of degree larger than 1.
integer, parameter :: NR_ITERATIONS = 8
! Tolerance for Newton-Raphson iterations (stop when increment falls below this)
Expand Down Expand Up @@ -355,9 +356,8 @@ subroutine checkGridsMatch( G, h, dzInterface )
'Non-zero dzInterface at bottom!')
enddo
enddo

end subroutine checkGridsMatch

end subroutine checkGridsMatch

!------------------------------------------------------------------------------
! Build uniform z*-ccordinate grid with partial steps
Expand Down Expand Up @@ -405,31 +405,11 @@ subroutine buildGridZstar( CS, G, h, dzInterface )
do k = 1,nz
totalThickness = totalThickness + h(i,j,k)
end do
minThickness = min( CS%min_thickness, totalThickness/float(nz) )

! Position of free-surface
eta = totalThickness - nominalDepth

! z* = (z-eta) / stretching where stretching = (H+eta)/H
! z = eta + stretching * z*
stretching = totalThickness / nominalDepth

! Integrate down from the top for a notional new grid, ignoring topography
zNew(1) = eta
do k = 1,nz
dh = stretching * CS%coordinateResolution(k) ! Notional grid spacing
zNew(k+1) = zNew(k) - dh
enddo
call buildGridZStarColumn(CS, nz, nominalDepth, totalThickness, zNew)

! The rest of the model defines grids integrating up from the bottom
zOld(nz+1) = - nominalDepth
zNew(nz+1) = - nominalDepth
do k = nz,1,-1
! Adjust interface position to accomodate inflating layers
! without disturbing the interface above
if ( zNew(k) < (zNew(k+1) + minThickness) ) then
zNew(k) = zNew(k+1) + minThickness
endif
zOld(k) = zOld(k+1) + h(i,j,k)
enddo

Expand Down Expand Up @@ -463,6 +443,46 @@ subroutine buildGridZstar( CS, G, h, dzInterface )

end subroutine buildGridZstar

subroutine buildGridZstarColumn( CS, nz, depth, totalThickness, zInterface)

! Arguments
type(regridding_CS), intent(in) :: CS
integer, intent(in) :: nz
real, intent(in) :: depth
real, intent(in) :: totalThickness
real, dimension(nz+1), intent(inout) :: zInterface

real :: eta, stretching, dh
real :: minThickness
integer :: k

minThickness = min( CS%min_thickness, totalThickness/float(nz) )

! Position of free-surface
eta = totalThickness - depth

! z* = (z-eta) / stretching where stretching = (H+eta)/H
! z = eta + stretching * z*
stretching = totalThickness / depth

! Integrate down from the top for a notional new grid, ignoring topography
zInterface(1) = eta
do k = 1,nz
dh = stretching * CS%coordinateResolution(k) ! Notional grid spacing
zInterface(k+1) = zInterface(k) - dh
enddo

! Integrating up from the bottom adjusting interface position to accomodate
! inflating layers without disturbing the interface above
zInterface(nz+1) = -depth
do k = nz,1,-1
if ( zInterface(k) < (zInterface(k+1) + minThickness) ) then
zInterface(k) = zInterface(k+1) + minThickness
endif
enddo

end subroutine buildGridZstarColumn


!------------------------------------------------------------------------------
! Build sigma grid
Expand Down
6 changes: 3 additions & 3 deletions src/ALE/MOM_remapping.F90
Original file line number Diff line number Diff line change
Expand Up @@ -400,7 +400,7 @@ subroutine remapping_core( CS, n0, h0, u0, n1, dx, u1 )
if (k <= n0) then
hTmp = h0(k) + ( dx(k+1) - dx(k) )
if (hTmp < 0.) then
write(0,*) 'k,h0(k),hNew,dx(+1),dx(0)=',k,h0(k),dx(k+1),dx(k)
write(0,*) 'k,h0(k),hTmp,dx(k+1),dx(k)=',k,h0(k),hTmp,dx(k+1),dx(k)
call MOM_error( FATAL, 'MOM_remapping, remapping_core: '//&
'negative h implied by fluxes' )
endif
Expand Down Expand Up @@ -534,9 +534,9 @@ subroutine remapping_core( CS, n0, h0, u0, n1, dx, u1 )
if (k<=n0) then; hTmp = h0(k); else; hTmp = 0.; endif
z0 = z0 + hTmp ; z1 = z1 + ( hTmp + ( dx(k+1) - dx(k) ) )
enddo
if (abs(totalHU2-totalHU0) > (err0+err2)*real(n1) .and. (err0+err2)/=0.) then
if (abs(totalHU2-totalHU0) > (err0+err2)*max(real(n0), real(n1)) .and. (err0+err2)/=0.) then
write(0,*) 'h0=',h0
write(0,*) 'hf=',h0+dx(2:n1+1)-dx(1:n1)
write(0,*) 'hf=',h0(1:n1)+dx(2:n1+1)-dx(1:n1)
write(0,*) 'u0=',u0
write(0,*) 'u1=',u1
write(0,*) 'total HU0,HUf,f-0=',totalHU0,totalHU2,totalHU2-totalHU0
Expand Down
6 changes: 5 additions & 1 deletion src/core/MOM.F90
Original file line number Diff line number Diff line change
Expand Up @@ -337,7 +337,7 @@ module MOM
use MOM_cpu_clock, only : CLOCK_COMPONENT, CLOCK_SUBCOMPONENT
use MOM_cpu_clock, only : CLOCK_MODULE_DRIVER, CLOCK_MODULE, CLOCK_ROUTINE
use MOM_coms, only : reproducing_sum
use MOM_diag_mediator, only : diag_mediator_init, enable_averaging
use MOM_diag_mediator, only : diag_mediator_init, enable_averaging, diag_set_thickness_ptr
use MOM_diag_mediator, only : disable_averaging, post_data, safe_alloc_ptr
use MOM_diag_mediator, only : register_diag_field, register_static_field
use MOM_diag_mediator, only : register_scalar_field
Expand Down Expand Up @@ -1920,6 +1920,10 @@ subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in)
endif
call callTree_waypoint("state variables allocated (initialize_MOM)")

! Set up a pointers h within diag mediator control structure,
! this needs to occur _after_ CS%h has been allocated.
call diag_set_thickness_ptr(CS%h, diag)

! Set the fields that are needed for bitwise identical restarting
! the time stepping scheme.
call restart_init(param_file, CS%restart_CSp)
Expand Down
Loading