14 implicit none ;
private
16 #include <MOM_memory.h>
31 subroutine find_eta_3d(h, tv, G, GV, US, eta, eta_bt, halo_size, eta_to_m)
35 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(in) :: h
38 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)+1),
intent(out) :: eta
40 real,
dimension(SZI_(G),SZJ_(G)),
optional,
intent(in) :: eta_bt
44 integer,
optional,
intent(in) :: halo_size
46 real,
optional,
intent(in) :: eta_to_m
50 real :: p(SZI_(G),SZJ_(G),SZK_(G)+1)
51 real :: dz_geo(SZI_(G),SZJ_(G),SZK_(G))
53 real :: dilate(SZI_(G))
56 real :: Z_to_eta, H_to_eta, H_to_rho_eta
57 integer i, j, k, isv, iev, jsv, jev, nz, halo
59 halo = 0 ;
if (
present(halo_size)) halo = max(0,halo_size)
61 isv = g%isc-halo ; iev = g%iec+halo ; jsv = g%jsc-halo ; jev = g%jec+halo
64 if ((isv<g%isd) .or. (iev>g%ied) .or. (jsv<g%jsd) .or. (jev>g%jed)) &
65 call mom_error(fatal,
"find_eta called with an overly large halo_size.")
67 z_to_eta = 1.0 ;
if (
present(eta_to_m)) z_to_eta = us%Z_to_m / eta_to_m
68 h_to_eta = gv%H_to_Z * z_to_eta
69 h_to_rho_eta = gv%H_to_RZ * z_to_eta
70 i_gearth = z_to_eta / (us%Z_to_m * gv%mks_g_Earth)
74 do j=jsv,jev ;
do i=isv,iev ; eta(i,j,nz+1) = -z_to_eta*g%bathyT(i,j) ;
enddo ;
enddo
76 if (gv%Boussinesq)
then
78 do j=jsv,jev ;
do k=nz,1,-1;
do i=isv,iev
79 eta(i,j,k) = eta(i,j,k+1) + h(i,j,k)*h_to_eta
81 if (
present(eta_bt))
then
87 dilate(i) = (eta_bt(i,j)*h_to_eta + z_to_eta*g%bathyT(i,j)) / &
88 (eta(i,j,1) + z_to_eta*g%bathyT(i,j))
90 do k=1,nz ;
do i=isv,iev
91 eta(i,j,k) = dilate(i) * (eta(i,j,k) + z_to_eta*g%bathyT(i,j)) - z_to_eta*g%bathyT(i,j)
96 if (
associated(tv%eqn_of_state))
then
100 do i=isv,iev ; p(i,j,1) = 0.0 ;
enddo
101 do k=1,nz ;
do i=isv,iev
102 p(i,j,k+1) = p(i,j,k) + gv%H_to_Pa*h(i,j,k)
107 call int_specific_vol_dp(tv%T(:,:,k), tv%S(:,:,k), p(:,:,k), p(:,:,k+1), &
108 0.0, g%HI, tv%eqn_of_state, dz_geo(:,:,k), halo_size=halo)
112 do k=nz,1,-1 ;
do i=isv,iev
113 eta(i,j,k) = eta(i,j,k+1) + i_gearth * dz_geo(i,j,k)
118 do j=jsv,jev ;
do k=nz,1,-1;
do i=isv,iev
119 eta(i,j,k) = eta(i,j,k+1) + h_to_rho_eta*h(i,j,k) / gv%Rlay(k)
120 enddo ;
enddo ;
enddo
122 if (
present(eta_bt))
then
127 do i=isv,iev ; htot(i) = gv%H_subroundoff ;
enddo
128 do k=1,nz ;
do i=isv,iev ; htot(i) = htot(i) + h(i,j,k) ;
enddo ;
enddo
129 do i=isv,iev ; dilate(i) = eta_bt(i,j) / htot(i) ;
enddo
130 do k=1,nz ;
do i=isv,iev
131 eta(i,j,k) = dilate(i) * (eta(i,j,k) + z_to_eta*g%bathyT(i,j)) - z_to_eta*g%bathyT(i,j)
144 subroutine find_eta_2d(h, tv, G, GV, US, eta, eta_bt, halo_size, eta_to_m)
148 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(in) :: h
151 real,
dimension(SZI_(G),SZJ_(G)),
intent(out) :: eta
153 real,
dimension(SZI_(G),SZJ_(G)),
optional,
intent(in) :: eta_bt
156 integer,
optional,
intent(in) :: halo_size
158 real,
optional,
intent(in) :: eta_to_m
161 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)+1) :: &
163 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)) :: &
165 real :: htot(SZI_(G))
167 real :: Z_to_eta, H_to_eta, H_to_rho_eta
168 integer i, j, k, is, ie, js, je, nz, halo
170 halo = 0 ;
if (
present(halo_size)) halo = max(0,halo_size)
171 is = g%isc-halo ; ie = g%iec+halo ; js = g%jsc-halo ; je = g%jec+halo
174 z_to_eta = 1.0 ;
if (
present(eta_to_m)) z_to_eta = us%Z_to_m / eta_to_m
175 h_to_eta = gv%H_to_Z * z_to_eta
176 h_to_rho_eta = gv%H_to_RZ * z_to_eta
177 i_gearth = z_to_eta / (us%Z_to_m * gv%mks_g_Earth)
181 do j=js,je ;
do i=is,ie ; eta(i,j) = -z_to_eta*g%bathyT(i,j) ;
enddo ;
enddo
183 if (gv%Boussinesq)
then
184 if (
present(eta_bt))
then
186 do j=js,je ;
do i=is,ie
187 eta(i,j) = h_to_eta*eta_bt(i,j)
191 do j=js,je ;
do k=1,nz ;
do i=is,ie
192 eta(i,j) = eta(i,j) + h(i,j,k)*h_to_eta
193 enddo ;
enddo ;
enddo
196 if (
associated(tv%eqn_of_state))
then
199 do i=is,ie ; p(i,j,1) = 0.0 ;
enddo
201 do k=1,nz ;
do i=is,ie
202 p(i,j,k+1) = p(i,j,k) + gv%H_to_Pa*h(i,j,k)
207 call int_specific_vol_dp(tv%T(:,:,k), tv%S(:,:,k), p(:,:,k), p(:,:,k+1), 0.0, &
208 g%HI, tv%eqn_of_state, dz_geo(:,:,k), halo_size=halo)
211 do j=js,je ;
do k=1,nz ;
do i=is,ie
212 eta(i,j) = eta(i,j) + i_gearth * dz_geo(i,j,k)
213 enddo ;
enddo ;
enddo
216 do j=js,je ;
do k=1,nz ;
do i=is,ie
217 eta(i,j) = eta(i,j) + h_to_rho_eta*h(i,j,k) / gv%Rlay(k)
218 enddo ;
enddo ;
enddo
220 if (
present(eta_bt))
then
225 do i=is,ie ; htot(i) = gv%H_subroundoff ;
enddo
226 do k=1,nz ;
do i=is,ie ; htot(i) = htot(i) + h(i,j,k) ;
enddo ;
enddo
228 eta(i,j) = (eta_bt(i,j) / htot(i)) * (eta(i,j) + z_to_eta*g%bathyT(i,j)) - &
229 z_to_eta*g%bathyT(i,j)