MOM_tracer_initialization_from_Z.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!> Initializes hydrography from z-coordinate climatology files
7
8use mom_debugging, only : hchksum
9use mom_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end
10use mom_cpu_clock, only : clock_routine, clock_loop
11use mom_domains, only : pass_var
12use mom_error_handler, only : mom_mesg, mom_error, fatal, warning
14use mom_file_parser, only : get_param, param_file_type, log_version
15use mom_grid, only : ocean_grid_type
16use mom_horizontal_regridding, only : mystats, horiz_interp_and_extrap_tracer
22use mom_ale, only : ale_remap_scalar
23
24implicit none ; private
25
26#include <MOM_memory.h>
27
29
30! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
31! consistency testing. These are noted in comments with units like Z, H, L, and T, along with
32! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the units
33! vary with the Boussinesq approximation, the Boussinesq variant is given first.
34
35character(len=40) :: mdl = "MOM_tracer_initialization_from_Z" !< This module's name.
36
37contains
38
39!> Initializes a tracer from a z-space data file, including any lateral regridding that is needed.
40subroutine mom_initialize_tracer_from_z(h, tr, G, GV, US, PF, src_file, src_var_nam, &
41 src_var_unit_conversion, src_var_record, homogenize, &
42 useALEremapping, remappingScheme, src_var_gridspec, h_in_Z_units, &
43 ongrid)
44 type(ocean_grid_type), intent(inout) :: g !< Ocean grid structure.
45 type(verticalgrid_type), intent(in) :: gv !< Ocean vertical grid structure.
46 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
47 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
48 intent(in) :: h !< Layer thicknesses, in [H ~> m or kg m-2] or
49 !! [Z ~> m] depending on the value of h_in_Z_units.
50 real, dimension(:,:,:), pointer :: tr !< Pointer to array to be initialized [CU ~> conc]
51 type(param_file_type), intent(in) :: pf !< parameter file
52 character(len=*), intent(in) :: src_file !< source filename
53 character(len=*), intent(in) :: src_var_nam !< variable name in file
54 real, optional, intent(in) :: src_var_unit_conversion !< optional multiplicative unit conversion,
55 !! often used for rescaling into model units [CU conc-1 ~> 1]
56 integer, optional, intent(in) :: src_var_record !< record to read for multiple time-level files
57 logical, optional, intent(in) :: homogenize !< optionally homogenize to mean value
58 logical, optional, intent(in) :: usealeremapping !< to remap or not (optional)
59 character(len=*), optional, intent(in) :: remappingscheme !< remapping scheme to use.
60 character(len=*), optional, intent(in) :: src_var_gridspec !< Source variable name in a gridspec file.
61 !! This is not implemented yet.
62 logical, optional, intent(in) :: h_in_z_units !< If present and true, the input grid
63 !! thicknesses are in the units of height
64 !! ([Z ~> m]) instead of the usual units of
65 !! thicknesses ([H ~> m or kg m-2])
66 logical, optional, intent(in) :: ongrid !< If true, then data are assumed to have been
67 !! interpolated to the model horizontal grid. In this case,
68 !! only extrapolation is performed by
69 !! horiz_interp_and_extrap_tracer()
70 ! Local variables
71 real :: land_fill = 0.0 ! A value to use to replace missing values [CU ~> conc]
72 real :: convert ! A conversion factor into the model's internal units [CU conc-1 ~> 1]
73 integer :: recnum
74 character(len=64) :: remapscheme
75 logical :: homog, useale
76 logical :: h_is_in_z_units
77
78 ! This include declares and sets the variable "version".
79# include "version_variable.h"
80 character(len=40) :: mdl = "MOM_initialize_tracers_from_Z" ! This module's name.
81
82 integer :: is, ie, js, je, nz ! compute domain indices
83 integer :: isd, ied, jsd, jed ! data domain indices
84 integer :: i, j, k, kd
85 real, allocatable, dimension(:,:,:), target :: tr_z ! Tracer array on the horizontal model grid
86 ! and input-file vertical levels [CU ~> conc]
87 real, allocatable, dimension(:,:,:), target :: mask_z ! Missing value mask on the horizontal model grid
88 ! and input-file vertical levels [nondim]
89 real, allocatable, dimension(:), target :: z_edges_in ! Cell edge depths for input data [Z ~> m]
90 real, allocatable, dimension(:), target :: z_in ! Cell center depths for input data [Z ~> m]
91
92 ! Local variables for ALE remapping
93 real, dimension(:,:,:), allocatable :: dzsrc ! Source thicknesses in height units [Z ~> m]
94 real, dimension(:,:,:), allocatable :: hsrc ! Source thicknesses [H ~> m or kg m-2]
95 real, dimension(:), allocatable :: h1 ! A 1-d column of source thicknesses [Z ~> m].
96 real :: ztopofcell, zbottomofcell, z_bathy ! Heights [Z ~> m].
97 type(remapping_cs) :: remapcs ! Remapping parameters and work arrays
98 type(verticalgrid_type) :: gv_loc ! A temporary vertical grid structure
99
100 real :: missing_value ! A value indicating that there is no valid input data at this point [CU ~> conc]
101 real :: dz_neglect ! A negligibly small vertical layer extent used in
102 ! remapping cell reconstructions [Z ~> m] or [H ~> m or kg m-2]
103 real :: dz_neglect_edge ! A negligibly small vertical layer extent used in
104 ! remapping edge value calculations [Z ~> m] or [H ~> m or kg m-2]
105 logical :: om4_remap_via_sub_cells ! If true, use the OM4 remapping algorithm
106 integer :: npoints ! The number of valid input data points in a column
107 integer :: id_clock_routine, id_clock_ale
108 integer :: default_answer_date ! The default setting for the various ANSWER_DATE flags.
109 integer :: remap_answer_date ! The vintage of the order of arithmetic and expressions to use
110 ! for remapping. Values below 20190101 recover the remapping
111 ! answers from 2018, while higher values use more robust
112 ! forms of the same remapping expressions.
113 integer :: hor_regrid_answer_date ! The vintage of the order of arithmetic and expressions to use
114 ! for horizontal regridding. Values below 20190101 recover the
115 ! answers from 2018, while higher values use expressions that have
116 ! been rearranged for rotational invariance.
117
118 id_clock_routine = cpu_clock_id('(Initialize tracer from Z)', grain=clock_routine)
119 id_clock_ale = cpu_clock_id('(Initialize tracer from Z) ALE', grain=clock_loop)
120
121 call cpu_clock_begin(id_clock_routine)
122
123 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
124 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
125
126 call calltree_enter(trim(mdl)//"(), MOM_tracer_initialization_from_Z.F90")
127
128 call get_param(pf, mdl, "Z_INIT_HOMOGENIZE", homog, &
129 "If True, then horizontally homogenize the interpolated "//&
130 "initial conditions.", default=.false.)
131 call get_param(pf, mdl, "Z_INIT_ALE_REMAPPING", useale, &
132 "If True, then remap straight to model coordinate from file.",&
133 default=.false.)
134 call get_param(pf, mdl, "Z_INIT_REMAPPING_SCHEME", remapscheme, &
135 "The remapping scheme to use if using Z_INIT_ALE_REMAPPING is True.", &
136 default="PPM_IH4")
137 call get_param(pf, mdl, "DEFAULT_ANSWER_DATE", default_answer_date, &
138 "This sets the default value for the various _ANSWER_DATE parameters.", &
139 default=99991231)
140 if (useale) then
141 call get_param(pf, mdl, "REMAPPING_ANSWER_DATE", remap_answer_date, &
142 "The vintage of the expressions and order of arithmetic to use for remapping. "//&
143 "Values below 20190101 result in the use of older, less accurate expressions "//&
144 "that were in use at the end of 2018. Higher values result in the use of more "//&
145 "robust and accurate forms of mathematically equivalent expressions.", &
146 default=default_answer_date, do_not_log=.not.gv%Boussinesq)
147 call get_param(pf, mdl, "REMAPPING_USE_OM4_SUBCELLS", om4_remap_via_sub_cells, &
148 do_not_log=.true., default=.true.)
149 call get_param(pf, mdl, "Z_INIT_REMAPPING_USE_OM4_SUBCELLS", om4_remap_via_sub_cells, &
150 "If true, use the OM4 remapping-via-subcells algorithm for initialization. "//&
151 "See REMAPPING_USE_OM4_SUBCELLS for more details. "//&
152 "We recommend setting this option to false.", default=om4_remap_via_sub_cells)
153 if (.not.gv%Boussinesq) remap_answer_date = max(remap_answer_date, 20230701)
154 endif
155 call get_param(pf, mdl, "HOR_REGRID_ANSWER_DATE", hor_regrid_answer_date, &
156 "The vintage of the order of arithmetic for horizontal regridding. "//&
157 "Dates before 20190101 give the same answers as the code did in late 2018, "//&
158 "while later versions add parentheses for rotational symmetry. "//&
159 "Dates after 20230101 use reproducing sums for global averages.", &
160 default=default_answer_date, do_not_log=.not.gv%Boussinesq)
161 if (.not.gv%Boussinesq) hor_regrid_answer_date = max(hor_regrid_answer_date, 20230701)
162
163 if (PRESENT(homogenize)) homog=homogenize
164 if (PRESENT(usealeremapping)) useale=usealeremapping
165 if (PRESENT(remappingscheme)) remapscheme=remappingscheme
166 recnum = 1
167 if (PRESENT(src_var_record)) recnum = src_var_record
168 convert = 1.0
169 if (PRESENT(src_var_unit_conversion)) convert = src_var_unit_conversion
170
171 h_is_in_z_units = .false. ; if (present(h_in_z_units)) h_is_in_z_units = h_in_z_units
172
173 call horiz_interp_and_extrap_tracer(src_file, src_var_nam, recnum, &
174 g, tr_z, mask_z, z_in, z_edges_in, missing_value, &
175 scale=convert, homogenize=homog, m_to_z=us%m_to_Z, &
176 answer_date=hor_regrid_answer_date, ongrid=ongrid)
177
178 kd = size(z_edges_in,1)-1
179 call pass_var(tr_z,g%Domain)
180 call pass_var(mask_z,g%Domain)
181
182! Done with horizontal interpolation.
183! Now remap to model coordinates
184 if (useale) then
185 call cpu_clock_begin(id_clock_ale)
186 ! First we reserve a work space for reconstructions of the source data
187 allocate( h1(kd) )
188 allocate( dzsrc(isd:ied,jsd:jed,kd) )
189 allocate( hsrc(isd:ied,jsd:jed,kd) )
190 ! Set parameters for reconstructions in the right units
191 if (h_is_in_z_units) then
192 dz_neglect = set_dz_neglect(gv, us, remap_answer_date, dz_neglect_edge)
193 else
194 dz_neglect = set_h_neglect(gv, remap_answer_date, dz_neglect_edge)
195 endif
196 call initialize_remapping( remapcs, remapscheme, boundary_extrapolation=.false., &
197 om4_remap_via_sub_cells=om4_remap_via_sub_cells, answer_date=remap_answer_date, &
198 h_neglect=dz_neglect, h_neglect_edge=dz_neglect_edge )
199 ! Next we initialize the regridding package so that it knows about the target grid
200
201 do j = js, je ; do i = is, ie
202 if (g%mask2dT(i,j)>0.) then
203 ! Build the source grid
204 ztopofcell = 0. ; zbottomofcell = 0. ; npoints = 0
205 z_bathy = g%bathyT(i,j) + g%Z_ref
206 do k = 1, kd
207 if (mask_z(i,j,k) > 0.) then
208 zbottomofcell = -min( z_edges_in(k+1), z_bathy )
209 elseif (k>1) then
210 zbottomofcell = -z_bathy
211 endif
212 h1(k) = ztopofcell - zbottomofcell
213 if (h1(k)>0.) npoints = npoints + 1
214 ztopofcell = zbottomofcell ! Bottom becomes top for next value of k
215 enddo
216 h1(kd) = h1(kd) + ( ztopofcell + z_bathy ) ! In case data is deeper than model
217 else
218 tr(i,j,:) = 0.
219 endif ! mask2dT
220 dzsrc(i,j,:) = h1(:)
221 enddo ; enddo
222
223 if (h_is_in_z_units) then
224 ! Because h is in units of [Z ~> m], dzSrc is already in the right units.
225 call ale_remap_scalar(remapcs, g, gv, kd, dzsrc, tr_z, h, tr, all_cells=.false.)
226 else
227 ! Equation of state data is not available, so a simpler rescaling will have to suffice,
228 ! but it might be problematic in non-Boussinesq mode.
229 gv_loc = gv ; gv_loc%ke = kd
230 call dz_to_thickness_simple(dzsrc, hsrc, g, gv_loc, us)
231
232 call ale_remap_scalar(remapcs, g, gv, kd, hsrc, tr_z, h, tr, all_cells=.false.)
233 endif
234
235 deallocate( hsrc )
236 deallocate( dzsrc )
237 deallocate( h1 )
238
239 do k=1,nz
240 call mystats(tr(:,:,k), missing_value, g, k, 'Tracer from ALE()')
241 enddo
242 call cpu_clock_end(id_clock_ale)
243 endif ! useALEremapping
244
245! Fill land values
246 do k=1,nz ; do j=js,je ; do i=is,ie
247 if (tr(i,j,k) == missing_value) then
248 tr(i,j,k) = land_fill
249 endif
250 enddo ; enddo ; enddo
251
252 call calltree_leave(trim(mdl)//'()')
253 call cpu_clock_end(id_clock_routine)
254
255end subroutine mom_initialize_tracer_from_z
256