portedportable, not yet portedexecuted, not portableexecutable, not hit by this run
| 1 | ! This file is part of MOM6, the Modular Ocean Model version 6. | |
| 2 | ! See the LICENSE file for licensing information. | |
| 3 | ! SPDX-License-Identifier: Apache-2.0 | |
| 4 | ||
| 5 | !> Regrid columns for a z-like coordinate (z-star, z-level) | |
| 6 | module coord_zlike | |
| 7 | ||
| 8 | use MOM_error_handler, only : MOM_error, FATAL | |
| 9 | ||
| 10 | implicit none ; private | |
| 11 | ||
| 12 | !> Control structure containing required parameters for a z-like coordinate | |
| 13 | type, public :: zlike_CS ; private | |
| 14 | ||
| 15 | !> Number of levels to be generated | |
| 16 | integer :: nk | |
| 17 | ||
| 18 | !> Minimum thickness allowed for layers, in the same thickness units (perhaps [H ~> m or kg m-2]) | |
| 19 | !! that will be used in all subsequent calls to build_zstar_column with this structure. | |
| 20 | real :: min_thickness | |
| 21 | ||
| 22 | !> Target coordinate resolution, usually in [Z ~> m] | |
| 23 | real, allocatable, dimension(:) :: coordinateResolution | |
| 24 | end type zlike_CS | |
| 25 | ||
| 26 | public init_coord_zlike, set_zlike_params, build_zstar_column, end_coord_zlike | |
| 27 | ||
| 28 | contains | |
| 29 | ||
| 30 | !> Initialise a zlike_CS with pointers to parameters | |
| 31 | 1 | subroutine init_coord_zlike(CS, nk, coordinateResolution) |
| 32 | type(zlike_CS), pointer :: CS !< Unassociated pointer to hold the control structure | |
| 33 | integer, intent(in) :: nk !< Number of levels in the grid | |
| 34 | real, dimension(:), intent(in) :: coordinateResolution !< Target coordinate resolution [Z ~> m] | |
| 35 | ||
| 36 | 1 | if (associated(CS)) call MOM_error(FATAL, "init_coord_zlike: CS already associated!") |
| 37 | 1 | allocate(CS) |
| 38 | 1 | allocate(CS%coordinateResolution(nk)) |
| 39 | ||
| 40 | 1 | CS%nk = nk |
| 41 | 35 | CS%coordinateResolution = coordinateResolution |
| 42 | 1 | end subroutine init_coord_zlike |
| 43 | ||
| 44 | !> Deallocates the zlike control structure | |
| 45 | 1 | subroutine end_coord_zlike(CS) |
| 46 | type(zlike_CS), pointer :: CS !< Coordinate control structure | |
| 47 | ||
| 48 | ! Nothing to do | |
| 49 | 1 | if (.not. associated(CS)) return |
| 50 | 1 | deallocate(CS%coordinateResolution) |
| 51 | 1 | deallocate(CS) |
| 52 | end subroutine end_coord_zlike | |
| 53 | ||
| 54 | !> Set parameters in the zlike structure | |
| 55 | 2 | subroutine set_zlike_params(CS, min_thickness) |
| 56 | type(zlike_CS), pointer :: CS !< Coordinate control structure | |
| 57 | real, optional, intent(in) :: min_thickness !< Minimum allowed thickness [H ~> m or kg m-2] | |
| 58 | ||
| 59 | 2 | if (.not. associated(CS)) call MOM_error(FATAL, "set_zlike_params: CS not associated") |
| 60 | ||
| 61 | 2 | if (present(min_thickness)) CS%min_thickness = min_thickness |
| 62 | 2 | end subroutine set_zlike_params |
| 63 | ||
| 64 | !> Builds a z* coordinate with a minimum thickness | |
| 65 | 46370 | subroutine build_zstar_column(CS, depth, total_thickness, zInterface, & |
| 66 | z_rigid_top, eta_orig, zScale) | |
| 67 | type(zlike_CS), intent(in) :: CS !< Coordinate control structure | |
| 68 | real, intent(in) :: depth !< Depth of ocean bottom (positive downward in the | |
| 69 | !! output units), units may be [Z ~> m] or [H ~> m or kg m-2] | |
| 70 | real, intent(in) :: total_thickness !< Column thickness (positive definite in the same | |
| 71 | !! units as depth) [Z ~> m] or [H ~> m or kg m-2] | |
| 72 | real, dimension(CS%nk+1), intent(inout) :: zInterface !< Absolute positions of interfaces (in the same | |
| 73 | !! units as depth) [Z ~> m] or [H ~> m or kg m-2] | |
| 74 | real, optional, intent(in) :: z_rigid_top !< The height of a rigid top (positive upward in the same | |
| 75 | !! units as depth) [Z ~> m] or [H ~> m or kg m-2] | |
| 76 | real, optional, intent(in) :: eta_orig !< The actual original height of the top (in the same | |
| 77 | !! units as depth) [Z ~> m] or [H ~> m or kg m-2] | |
| 78 | real, optional, intent(in) :: zScale !< Scaling factor from the target coordinate resolution | |
| 79 | !! in Z to desired units for zInterface, perhaps Z_to_H, | |
| 80 | !! often [nondim] or [H Z-1 ~> 1 or kg m-3] | |
| 81 | ! Local variables | |
| 82 | real :: eta ! Free surface height [Z ~> m] or [H ~> m or kg m-2] | |
| 83 | real :: stretching ! A stretching factor for the coordinate [nondim] | |
| 84 | real :: dh, min_thickness, z0_top, z_star, z_scale ! Thicknesses or heights [Z ~> m] or [H ~> m or kg m-2] | |
| 85 | integer :: k | |
| 86 | logical :: new_zstar_def | |
| 87 | ||
| 88 | 46370 | z_scale = 1.0 ; if (present(zScale)) z_scale = zScale |
| 89 | ||
| 90 | 46370 | new_zstar_def = .false. |
| 91 | 46370 | min_thickness = min( CS%min_thickness, total_thickness/real(CS%nk) ) |
| 92 | 46370 | z0_top = 0. |
| 93 | 46370 | if (present(z_rigid_top)) then |
| 94 | 0 | z0_top = z_rigid_top |
| 95 | 0 | new_zstar_def = .true. |
| 96 | endif | |
| 97 | ||
| 98 | ! Position of free-surface (or the rigid top, for which eta ~ z0_top) | |
| 99 | 46370 | eta = total_thickness - depth |
| 100 | 46370 | if (present(eta_orig)) eta = eta_orig |
| 101 | ||
| 102 | ! Conventional z* coordinate: | |
| 103 | ! z* = (z-eta) / stretching where stretching = (H+eta)/H | |
| 104 | ! z = eta + stretching * z* | |
| 105 | ! The above gives z*(z=eta) = 0, z*(z=-H) = -H. | |
| 106 | ! With a rigid top boundary at eta = z0_top then | |
| 107 | ! z* = z0 + (z-eta) / stretching where stretching = (H+eta)/(H+z0) | |
| 108 | ! z = eta + stretching * (z*-z0) * stretching | |
| 109 | 46370 | stretching = total_thickness / ( depth + z0_top ) |
| 110 | ||
| 111 | 46370 | if (new_zstar_def) then |
| 112 | ! z_star is the notional z* coordinate in absence of upper/lower topography | |
| 113 | 0 | z_star = 0. ! z*=0 at the free-surface |
| 114 | 0 | zInterface(1) = eta ! The actual position of the top of the column |
| 115 | 0 | do k = 2,CS%nk |
| 116 | 0 | z_star = z_star - CS%coordinateResolution(k-1)*z_scale |
| 117 | ! This ensures that z is below a rigid upper surface (ice shelf bottom) | |
| 118 | 0 | zInterface(k) = min( eta + stretching * ( z_star - z0_top ), z0_top ) |
| 119 | ! This ensures that the layer in inflated | |
| 120 | 0 | zInterface(k) = min( zInterface(k), zInterface(k-1) - min_thickness ) |
| 121 | ! This ensures that z is above or at the topography | |
| 122 | 0 | zInterface(k) = max( zInterface(k), -depth + real(CS%nk+1-k) * min_thickness ) |
| 123 | enddo | |
| 124 | 0 | zInterface(CS%nk+1) = -depth |
| 125 | ||
| 126 | else | |
| 127 | ! Integrate down from the top for a notional new grid, ignoring topography | |
| 128 | ! The starting position is offset by z0_top which, if z0_top<0, will place | |
| 129 | ! interfaces above the rigid boundary. | |
| 130 | 46370 | zInterface(1) = eta |
| 131 | 1576580 | do k = 1,CS%nk |
| 132 | 1530210 | dh = stretching * CS%coordinateResolution(k)*z_scale ! Notional grid spacing |
| 133 | 1576580 | zInterface(k+1) = zInterface(k) - dh |
| 134 | enddo | |
| 135 | ||
| 136 | ! Integrating up from the bottom adjusting interface position to accommodate | |
| 137 | ! inflating layers without disturbing the interface above | |
| 138 | 46370 | zInterface(CS%nk+1) = -depth |
| 139 | 1576580 | do k = CS%nk,1,-1 |
| 140 | 1576580 | if ( zInterface(k) < (zInterface(k+1) + min_thickness) ) then |
| 141 | 352458 | zInterface(k) = zInterface(k+1) + min_thickness |
| 142 | endif | |
| 143 | enddo | |
| 144 | endif | |
| 145 | ||
| 146 | 46370 | end subroutine build_zstar_column |
| 147 | ||
| 148 | 0 | end module coord_zlike |