This subroutine initializes layer thicknesses for the Neverland test case, by finding the depths of interfaces in a specified latitude-dependent temperature profile with an exponentially decaying thermocline on top of a linear stratification.
114 type(ocean_grid_type),
intent(in) :: G
115 type(verticalGrid_type),
intent(in) :: GV
116 type(unit_scale_type),
intent(in) :: US
117 real,
intent(out),
dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: h
119 type(param_file_type),
intent(in) :: param_file
122 type(EOS_type),
pointer :: eqn_of_state
124 real,
intent(in) :: P_Ref
127 real :: e0(SZK_(G)+1)
129 real,
dimension(SZK_(G)) :: h_profile
135 type(randomNumberStream) :: rns
136 character(len=40) :: mdl =
"Neverland_initialize_thickness"
137 integer :: i, j, k, k1, is, ie, js, je, nz, itt
139 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
141 call mom_mesg(
" Neverland_initialization.F90, Neverland_initialize_thickness: setting thickness", 5)
142 call get_param(param_file, mdl,
"INIT_THICKNESS_PROFILE", h_profile, &
143 "Profile of initial layer thicknesses.", units=
"m", scale=us%m_to_Z, &
144 fail_if_missing=.true.)
145 call get_param(param_file, mdl,
"NL_THICKNESS_PERT_AMP", pert_amp, &
146 "Amplitude of finite scale perturbations as fraction of depth.", &
147 units=
"nondim", default=0.)
148 call get_param(param_file, mdl,
"NL_THICKNESS_NOISE_AMP", h_noise, &
149 "Amplitude of noise to scale layer by.", units=
"nondim", default=0.)
154 e0(k+1) = e0(k) - h_profile(k)
157 do j=js,je ;
do i=is,ie
158 e_interface = -g%bathyT(i,j)
160 h(i,j,k) = gv%Z_to_H * (e0(k) - e_interface)
161 x=(g%geoLonT(i,j)-g%west_lon)/g%len_lon
162 y=(g%geoLatT(i,j)-g%south_lat)/g%len_lat
163 r1=sqrt((x-0.7)**2+(y-0.2)**2)
164 r2=sqrt((x-0.3)**2+(y-0.25)**2)
165 h(i,j,k) = h(i,j,k) + pert_amp * (e0(k) - e0(nz+1)) * gv%Z_to_H * &
166 (spike(r1,0.15)-spike(r2,0.15))
167 if (h_noise /= 0.)
then
168 rns = initializerandomnumberstream( int( 4096*(x + (y+1.)) ) )
169 call getrandomnumbers(rns, noise)
170 noise = h_noise * 2. * ( noise - 0.5 )
171 h(i,j,k) = ( 1. + noise ) * h(i,j,k)
173 h(i,j,k) = max( gv%Angstrom_H, h(i,j,k) )
174 e_interface = e_interface + gv%H_to_Z * h(i,j,k)
176 h(i,j,1) = gv%Z_to_H * (e0(1) - e_interface)
177 h(i,j,1) = max( gv%Angstrom_H, h(i,j,1) )