17 implicit none ;
private
19 #include <MOM_memory.h>
21 public circle_obcs_initialize_thickness
31 subroutine circle_obcs_initialize_thickness(h, G, GV, param_file, just_read_params)
34 real,
dimension(SZI_(G),SZJ_(G), SZK_(GV)), &
38 logical,
optional,
intent(in) :: just_read_params
41 real :: e0(szk_(gv)+1)
43 real :: eta1d(szk_(gv)+1)
46 real :: diskrad, rad, xcenter, xradius, lonc, latc, xoffset
49 # include "version_variable.h"
50 character(len=40) :: mdl =
"circle_obcs_initialization"
51 integer :: i, j, k, is, ie, js, je, nz
53 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
55 just_read = .false. ;
if (
present(just_read_params)) just_read = just_read_params
58 call mom_mesg(
" circle_obcs_initialization.F90, circle_obcs_initialize_thickness: setting thickness", 5)
60 if (.not.just_read)
call log_version(param_file, mdl, version,
"")
62 call get_param(param_file, mdl,
"DISK_RADIUS", diskrad, &
63 "The radius of the initially elevated disk in the "//&
64 "circle_obcs test case.", units=g%x_axis_units, &
65 fail_if_missing=.not.just_read, do_not_log=just_read)
66 call get_param(param_file, mdl,
"DISK_X_OFFSET", xoffset, &
67 "The x-offset of the initially elevated disk in the "//&
68 "circle_obcs test case.", units=g%x_axis_units, &
69 default = 0.0, do_not_log=just_read)
70 call get_param(param_file, mdl,
"DISK_IC_AMPLITUDE", ic_amp, &
71 "Initial amplitude of interface height displacements "//&
72 "in the circle_obcs test case.", &
73 units=
'm', default=5.0, scale=gv%m_to_H, do_not_log=just_read)
78 e0(k) = -g%max_depth * real(k-1) / real(nz)
82 do j=js,je ;
do i=is,ie
83 eta1d(nz+1) = -g%bathyT(i,j)
86 if (eta1d(k) < (eta1d(k+1) + gv%Angstrom_Z))
then
87 eta1d(k) = eta1d(k+1) + gv%Angstrom_Z
88 h(i,j,k) = gv%Angstrom_H
90 h(i,j,k) = gv%Z_to_H * (eta1d(k) - eta1d(k+1))
97 latc = g%south_lat + 0.5*g%len_lat
98 lonc = g%west_lon + 0.5*g%len_lon + xoffset
99 do j=js,je ;
do i=is,ie
100 rad = sqrt((g%geoLonT(i,j)-lonc)**2+(g%geoLatT(i,j)-latc)**2)/(diskrad)
103 rad = rad*(2.*asin(1.))
106 h(i,j,k) = h(i,j,k) + ic_amp * 0.5*(1.+cos(rad))
110 h(i,j,k) = h(i,j,k) - 0.5*(1.+cos(rad)) * ic_amp * real( 2*k-nz )
115 end subroutine circle_obcs_initialize_thickness