circle_obcs_initialization.F90

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!> Configures the model for the "circle_obcs" experiment which tests
6!! Open Boundary Conditions radiating an SSH anomaly.
8
10use mom_error_handler, only : mom_mesg, mom_error, fatal, is_root_pe
11use mom_file_parser, only : get_param, log_version, param_file_type
12use mom_get_input, only : directories
13use mom_grid, only : ocean_grid_type
14use mom_tracer_registry, only : tracer_registry_type
18
19implicit none ; private
20
21#include <MOM_memory.h>
22
24
25! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
26! consistency testing. These are noted in comments with units like Z, H, L, and T, along with
27! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the units
28! vary with the Boussinesq approximation, the Boussinesq variant is given first.
29
30contains
31
32!> This subroutine initializes layer thicknesses for the circle_obcs experiment.
33subroutine circle_obcs_initialize_thickness(h, depth_tot, G, GV, US, param_file, just_read)
34 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
35 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
36 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
37 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
38 intent(out) :: h !< The thickness that is being initialized [Z ~> m]
39 real, dimension(SZI_(G),SZJ_(G)), &
40 intent(in) :: depth_tot !< The nominal total depth of the ocean [Z ~> m]
41 type(param_file_type), intent(in) :: param_file !< A structure indicating the open file
42 !! to parse for model parameter values.
43 logical, intent(in) :: just_read !< If true, this call will only read
44 !! parameters without changing h.
45
46 real :: e0(szk_(gv)+1) ! The resting interface heights, in depth units [Z ~> m], usually
47 ! negative because it is positive upward.
48 real :: eta1d(szk_(gv)+1)! Interface height relative to the sea surface
49 ! positive upward, in depth units [Z ~> m].
50 real :: ic_amp ! The amplitude of the initial height displacement [Z ~> m].
51 real :: diskrad ! Radius of the elevated disk [km] or [degrees] or [m]
52 real :: rad ! Distance from the center of the elevated disk [km] or [degrees] or [m]
53 real :: lonc ! The x-position of a point [km] or [degrees] or [m]
54 real :: latc ! The y-position of a point [km] or [degrees] or [m]
55 real :: xoffset ! The x-offset of the elevated disc center relative to the domain
56 ! center [km] or [degrees] or [m]
57 ! This include declares and sets the variable "version".
58# include "version_variable.h"
59 character(len=40) :: mdl = "circle_obcs_initialization" ! This module's name.
60 integer :: i, j, k, is, ie, js, je, nz
61
62 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
63
64 if (.not.just_read) &
65 call mom_mesg(" circle_obcs_initialization.F90, circle_obcs_initialize_thickness: setting thickness", 5)
66
67 if (.not.just_read) call log_version(param_file, mdl, version, "")
68 ! Parameters read by cartesian grid initialization
69 call get_param(param_file, mdl, "DISK_RADIUS", diskrad, &
70 "The radius of the initially elevated disk in the "//&
71 "circle_obcs test case.", units=g%x_ax_unit_short, &
72 fail_if_missing=.not.just_read, do_not_log=just_read)
73 call get_param(param_file, mdl, "DISK_X_OFFSET", xoffset, &
74 "The x-offset of the initially elevated disk in the "//&
75 "circle_obcs test case.", units=g%x_ax_unit_short, &
76 default=0.0, do_not_log=just_read)
77 call get_param(param_file, mdl, "DISK_IC_AMPLITUDE", ic_amp, &
78 "Initial amplitude of interface height displacements "//&
79 "in the circle_obcs test case.", &
80 units='m', default=5.0, scale=us%m_to_Z, do_not_log=just_read)
81
82 if (just_read) return ! All run-time parameters have been read, so return.
83
84 do k=1,nz
85 e0(k) = -g%max_depth * real(k-1) / real(nz)
86 enddo
87
88 ! Uniform thicknesses for base state
89 do j=js,je ; do i=is,ie
90 eta1d(nz+1) = -depth_tot(i,j)
91 do k=nz,1,-1
92 eta1d(k) = e0(k)
93 if (eta1d(k) < (eta1d(k+1) + gv%Angstrom_Z)) then
94 eta1d(k) = eta1d(k+1) + gv%Angstrom_Z
95 h(i,j,k) = gv%Angstrom_Z
96 else
97 h(i,j,k) = eta1d(k) - eta1d(k+1)
98 endif
99 enddo
100 enddo ; enddo
101
102 ! Perturb base state by circular anomaly in center
103 k=nz
104 latc = g%south_lat + 0.5*g%len_lat
105 lonc = g%west_lon + 0.5*g%len_lon + xoffset
106 do j=js,je ; do i=is,ie
107 rad = sqrt(((g%geoLonT(i,j)-lonc)**2) + ((g%geoLatT(i,j)-latc)**2)) / diskrad
108 ! if (rad <= 6.*diskrad) h(i,j,k) = h(i,j,k)+10.0*exp( -0.5*( rad**2 ) )
109 rad = min( rad, 1. ) ! Flatten outside radius of diskrad
110 rad = rad*(2.*asin(1.)) ! Map 0-1 to 0-pi
111 if (nz==1) then
112 ! The model is barotropic
113 h(i,j,k) = h(i,j,k) + ic_amp * 0.5*(1.+cos(rad)) ! cosine bell
114 else
115 ! The model is baroclinic
116 do k = 1, nz
117 h(i,j,k) = h(i,j,k) - 0.5*(1.+cos(rad)) * ic_amp * real( 2*k-nz )
118 enddo
119 endif
120 enddo ; enddo
121
123