9 use mom_domains,
only : root_pe, to_all, scalar_pair, cgrid_ne, agrid
22 implicit none ;
private
24 public mom_shared_init_init
25 public mom_initialize_rotation, mom_calculate_grad_coriolis
26 public initialize_topography_from_file, apply_topography_edits_from_file
27 public initialize_topography_named, limit_topography, diagnosemaximumdepth
28 public set_rotation_planetary, set_rotation_beta_plane, initialize_grid_rotation_angle
29 public reset_face_lengths_named, reset_face_lengths_file, reset_face_lengths_list
30 public read_face_length_list, set_velocity_depth_max, set_velocity_depth_min
31 public compute_global_grid_integrals, write_ocean_geometry_file
42 subroutine mom_shared_init_init(PF)
46 character(len=40) :: mdl =
"MOM_shared_initialization"
49 #include "version_variable.h"
51 "Sharable code to initialize time-invariant fields, like bathymetry and Coriolis parameters.")
53 end subroutine mom_shared_init_init
57 subroutine mom_initialize_rotation(f, G, PF, US)
59 real,
dimension(G%IsdB:G%IedB,G%JsdB:G%JedB),
intent(out) :: f
67 character(len=40) :: mdl =
"MOM_initialize_rotation"
68 character(len=200) :: config
70 call calltree_enter(trim(mdl)//
"(), MOM_shared_initialization.F90")
71 call get_param(pf, mdl,
"ROTATION", config, &
72 "This specifies how the Coriolis parameter is specified: \n"//&
73 " \t 2omegasinlat - Use twice the planetary rotation rate \n"//&
74 " \t\t times the sine of latitude.\n"//&
75 " \t betaplane - Use a beta-plane or f-plane.\n"//&
76 " \t USER - call a user modified routine.", &
77 default=
"2omegasinlat")
78 select case (trim(config))
79 case (
"2omegasinlat");
call set_rotation_planetary(f, g, pf, us)
80 case (
"beta");
call set_rotation_beta_plane(f, g, pf, us)
81 case (
"betaplane");
call set_rotation_beta_plane(f, g, pf, us)
83 case default ;
call mom_error(fatal,
"MOM_initialize: "// &
84 "Unrecognized rotation setup "//trim(config))
86 call calltree_leave(trim(mdl)//
'()')
87 end subroutine mom_initialize_rotation
90 subroutine mom_calculate_grad_coriolis(dF_dx, dF_dy, G, US)
92 real,
dimension(G%isd:G%ied,G%jsd:G%jed), &
94 real,
dimension(G%isd:G%ied,G%jsd:G%jed), &
101 if ((lbound(g%CoriolisBu,1) > g%isc-1) .or. &
102 (lbound(g%CoriolisBu,2) > g%isc-1))
then
104 df_dx(:,:) = 0.0 ; df_dy(:,:) = 0.0
108 do j=g%jsc, g%jec ;
do i=g%isc, g%iec
109 f1 = 0.5*( g%CoriolisBu(i,j) + g%CoriolisBu(i,j-1) )
110 f2 = 0.5*( g%CoriolisBu(i-1,j) + g%CoriolisBu(i-1,j-1) )
111 df_dx(i,j) = g%IdxT(i,j) * ( f1 - f2 )
112 f1 = 0.5*( g%CoriolisBu(i,j) + g%CoriolisBu(i-1,j) )
113 f2 = 0.5*( g%CoriolisBu(i,j-1) + g%CoriolisBu(i-1,j-1) )
114 df_dy(i,j) = g%IdyT(i,j) * ( f1 - f2 )
116 call pass_vector(df_dx, df_dy, g%Domain, stagger=agrid)
117 end subroutine mom_calculate_grad_coriolis
120 function diagnosemaximumdepth(D, G)
122 real,
dimension(G%isd:G%ied,G%jsd:G%jed), &
124 real :: diagnosemaximumdepth
127 diagnosemaximumdepth = d(g%isc,g%jsc)
128 do j=g%jsc, g%jec ;
do i=g%isc, g%iec
129 diagnosemaximumdepth = max(diagnosemaximumdepth,d(i,j))
131 call max_across_pes(diagnosemaximumdepth)
132 end function diagnosemaximumdepth
136 subroutine initialize_topography_from_file(D, G, param_file, US)
138 real,
dimension(G%isd:G%ied,G%jsd:G%jed), &
144 character(len=200) :: filename, topo_file, inputdir
145 character(len=200) :: topo_varname
146 character(len=40) :: mdl =
"initialize_topography_from_file"
148 call calltree_enter(trim(mdl)//
"(), MOM_shared_initialization.F90")
150 m_to_z = 1.0 ;
if (
present(us)) m_to_z = us%m_to_Z
152 call get_param(param_file, mdl,
"INPUTDIR", inputdir, default=
".")
153 inputdir = slasher(inputdir)
154 call get_param(param_file, mdl,
"TOPO_FILE", topo_file, &
155 "The file from which the bathymetry is read.", &
157 call get_param(param_file, mdl,
"TOPO_VARNAME", topo_varname, &
158 "The name of the bathymetry variable in TOPO_FILE.", &
161 filename = trim(inputdir)//trim(topo_file)
162 call log_param(param_file, mdl,
"INPUTDIR/TOPO_FILE", filename)
164 if (.not.
file_exists(filename, g%Domain))
call mom_error(fatal, &
165 " initialize_topography_from_file: Unable to open "//trim(filename))
167 d(:,:) = -9.e30*m_to_z
173 call mom_read_data(filename, trim(topo_varname), d, g%Domain, scale=m_to_z)
175 call apply_topography_edits_from_file(d, g, param_file, us)
177 call calltree_leave(trim(mdl)//
'()')
178 end subroutine initialize_topography_from_file
181 subroutine apply_topography_edits_from_file(D, G, param_file, US)
183 real,
dimension(G%isd:G%ied,G%jsd:G%jed), &
190 character(len=200) :: topo_edits_file, inputdir
191 character(len=40) :: mdl =
"apply_topography_edits_from_file"
192 integer :: n_edits, n, ashape(5), i, j, ncid, id, ncstatus, iid, jid, zid
193 integer,
dimension(:),
allocatable :: ig, jg
194 real,
dimension(:),
allocatable :: new_depth
196 call calltree_enter(trim(mdl)//
"(), MOM_shared_initialization.F90")
198 m_to_z = 1.0 ;
if (
present(us)) m_to_z = us%m_to_Z
200 call get_param(param_file, mdl,
"INPUTDIR", inputdir, default=
".")
201 inputdir = slasher(inputdir)
202 call get_param(param_file, mdl,
"TOPO_EDITS_FILE", topo_edits_file, &
203 "The file from which to read a list of i,j,z topography overrides.", &
206 if (len_trim(topo_edits_file)==0)
return
208 topo_edits_file = trim(inputdir)//trim(topo_edits_file)
209 if (.not.
file_exists(topo_edits_file, g%Domain))
call mom_error(fatal, &
210 'initialize_topography_from_file: Unable to open '//trim(topo_edits_file))
212 ncstatus = nf90_open(trim(topo_edits_file), nf90_nowrite, ncid)
213 if (ncstatus /= nf90_noerr)
call mom_error(fatal,
'apply_topography_edits_from_file: '//&
214 'Failed to open '//trim(topo_edits_file))
217 ncstatus = nf90_inq_dimid(ncid,
'nEdits', id)
218 if (ncstatus /= nf90_noerr)
call mom_error(fatal,
'apply_topography_edits_from_file: '//&
219 'Failed to inq_dimid nEdits for '//trim(topo_edits_file))
220 ncstatus = nf90_inquire_dimension(ncid, id, len=n_edits)
221 if (ncstatus /= nf90_noerr)
call mom_error(fatal,
'apply_topography_edits_from_file: '//&
222 'Failed to inquire_dimension nEdits for '//trim(topo_edits_file))
225 ncstatus = nf90_inq_varid(ncid,
'ni', id)
226 if (ncstatus /= nf90_noerr)
call mom_error(fatal,
'apply_topography_edits_from_file: '//&
227 'Failed to inq_varid ni for '//trim(topo_edits_file))
228 ncstatus = nf90_get_var(ncid, id, i)
229 if (ncstatus /= nf90_noerr)
call mom_error(fatal,
'apply_topography_edits_from_file: '//&
230 'Failed to get_var ni for '//trim(topo_edits_file))
231 if (i /= g%ieg)
call mom_error(fatal,
'apply_topography_edits_from_file: '//&
232 'Incompatible i-dimension of grid in '//trim(topo_edits_file))
235 ncstatus = nf90_inq_varid(ncid,
'nj', id)
236 if (ncstatus /= nf90_noerr)
call mom_error(fatal,
'apply_topography_edits_from_file: '//&
237 'Failed to inq_varid nj for '//trim(topo_edits_file))
238 ncstatus = nf90_get_var(ncid, id, j)
239 if (ncstatus /= nf90_noerr)
call mom_error(fatal,
'apply_topography_edits_from_file: '//&
240 'Failed to get_var nj for '//trim(topo_edits_file))
241 if (j /= g%jeg)
call mom_error(fatal,
'apply_topography_edits_from_file: '//&
242 'Incompatible j-dimension of grid in '//trim(topo_edits_file))
245 ncstatus = nf90_inq_varid(ncid,
'iEdit', id)
246 if (ncstatus /= nf90_noerr)
call mom_error(fatal,
'apply_topography_edits_from_file: '//&
247 'Failed to inq_varid iEdit for '//trim(topo_edits_file))
248 allocate(ig(n_edits))
249 ncstatus = nf90_get_var(ncid, id, ig)
250 if (ncstatus /= nf90_noerr)
call mom_error(fatal,
'apply_topography_edits_from_file: '//&
251 'Failed to get_var iEdit for '//trim(topo_edits_file))
254 ncstatus = nf90_inq_varid(ncid,
'jEdit', id)
255 if (ncstatus /= nf90_noerr)
call mom_error(fatal,
'apply_topography_edits_from_file: '//&
256 'Failed to inq_varid jEdit for '//trim(topo_edits_file))
257 allocate(jg(n_edits))
258 ncstatus = nf90_get_var(ncid, id, jg)
259 if (ncstatus /= nf90_noerr)
call mom_error(fatal,
'apply_topography_edits_from_file: '//&
260 'Failed to get_var jEdit for '//trim(topo_edits_file))
263 ncstatus = nf90_inq_varid(ncid,
'zEdit', id)
264 if (ncstatus /= nf90_noerr)
call mom_error(fatal,
'apply_topography_edits_from_file: '//&
265 'Failed to inq_varid zEdit for '//trim(topo_edits_file))
266 allocate(new_depth(n_edits))
267 ncstatus = nf90_get_var(ncid, id, new_depth)
268 if (ncstatus /= nf90_noerr)
call mom_error(fatal,
'apply_topography_edits_from_file: '//&
269 'Failed to get_var zEdit for '//trim(topo_edits_file))
272 ncstatus = nf90_close(ncid)
273 if (ncstatus /= nf90_noerr)
call mom_error(fatal,
'apply_topography_edits_from_file: '//&
274 'Failed to close '//trim(topo_edits_file))
277 i = ig(n) - g%isd_global + 2
278 j = jg(n) - g%jsd_global + 2
279 if (i>=g%isc .and. i<=g%iec .and. j>=g%jsc .and. j<=g%jec)
then
280 if (new_depth(n)/=0.)
then
281 write(*,
'(a,3i5,f8.2,a,f8.2,2i4)') &
282 'Ocean topography edit: ',n,ig(n),jg(n),d(i,j)/m_to_z,
'->',abs(new_depth(n)),i,j
283 d(i,j) = abs(m_to_z*new_depth(n))
285 call mom_error(fatal,
' apply_topography_edits_from_file: '//&
286 "A zero depth edit would change the land mask and is not allowed in"//trim(topo_edits_file))
291 deallocate( ig, jg, new_depth )
293 call calltree_leave(trim(mdl)//
'()')
294 end subroutine apply_topography_edits_from_file
297 subroutine initialize_topography_named(D, G, param_file, topog_config, max_depth, US)
299 real,
dimension(G%isd:G%ied,G%jsd:G%jed), &
302 character(len=*),
intent(in) :: topog_config
304 real,
intent(in) :: max_depth
317 integer :: i, j, is, ie, js, je, isd, ied, jsd, jed
318 character(len=40) :: mdl =
"initialize_topography_named"
319 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
320 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
322 call calltree_enter(trim(mdl)//
"(), MOM_shared_initialization.F90")
323 call mom_mesg(
" MOM_shared_initialization.F90, initialize_topography_named: "//&
324 "TOPO_CONFIG = "//trim(topog_config), 5)
326 m_to_z = 1.0 ;
if (
present(us)) m_to_z = us%m_to_Z
328 call get_param(param_file, mdl,
"MINIMUM_DEPTH", min_depth, &
329 "The minimum depth of the ocean.", units=
"m", default=0.0, scale=m_to_z)
330 if (max_depth<=0.)
call mom_error(fatal,
"initialize_topography_named: "// &
331 "MAXIMUM_DEPTH has a non-sensical value! Was it set?")
333 if (trim(topog_config) /=
"flat")
then
334 call get_param(param_file, mdl,
"EDGE_DEPTH", dedge, &
335 "The depth at the edge of one of the named topographies.", &
336 units=
"m", default=100.0, scale=m_to_z)
351 call get_param(param_file, mdl,
"TOPOG_SLOPE_SCALE", expdecay, &
352 "The exponential decay scale used in defining some of "//&
353 "the named topographies.", units=
"m", default=400000.0)
359 if (trim(topog_config) ==
"flat")
then
360 do i=is,ie ;
do j=js,je ; d(i,j) = max_depth ;
enddo ;
enddo
361 elseif (trim(topog_config) ==
"spoon")
then
362 d0 = (max_depth - dedge) / &
363 ((1.0 - exp(-0.5*g%len_lat*g%Rad_earth*pi/(180.0 *expdecay))) * &
364 (1.0 - exp(-0.5*g%len_lat*g%Rad_earth*pi/(180.0 *expdecay))))
365 do i=is,ie ;
do j=js,je
368 d(i,j) = dedge + d0 * &
369 (sin(pi * (g%geoLonT(i,j) - (g%west_lon)) / g%len_lon) * &
370 (1.0 - exp((g%geoLatT(i,j) - (g%south_lat+g%len_lat))*g%Rad_earth*pi / &
373 elseif (trim(topog_config) ==
"bowl")
then
374 d0 = (max_depth - dedge) / &
375 ((1.0 - exp(-0.5*g%len_lat*g%Rad_earth*pi/(180.0 *expdecay))) * &
376 (1.0 - exp(-0.5*g%len_lat*g%Rad_earth*pi/(180.0 *expdecay))))
380 do i=is,ie ;
do j=js,je
381 d(i,j) = dedge + d0 * &
382 (sin(pi * (g%geoLonT(i,j) - g%west_lon) / g%len_lon) * &
383 ((1.0 - exp(-(g%geoLatT(i,j) - g%south_lat)*g%Rad_Earth*pi/ &
384 (180.0*expdecay))) * &
385 (1.0 - exp((g%geoLatT(i,j) - (g%south_lat+g%len_lat))* &
386 g%Rad_Earth*pi/(180.0*expdecay)))))
388 elseif (trim(topog_config) ==
"halfpipe")
then
389 d0 = max_depth - dedge
390 do i=is,ie ;
do j=js,je
391 d(i,j) = dedge + d0 * abs(sin(pi*(g%geoLatT(i,j) - g%south_lat)/g%len_lat))
394 call mom_error(fatal,
"initialize_topography_named: "// &
395 "Unrecognized topography name "//trim(topog_config))
399 do i=is,ie ;
do j=js,je
400 if (d(i,j) > max_depth) d(i,j) = max_depth
401 if (d(i,j) < min_depth) d(i,j) = 0.5*min_depth
404 call calltree_leave(trim(mdl)//
'()')
405 end subroutine initialize_topography_named
410 subroutine limit_topography(D, G, param_file, max_depth, US)
412 real,
dimension(G%isd:G%ied,G%jsd:G%jed), &
415 real,
intent(in) :: max_depth
421 character(len=40) :: mdl =
"limit_topography"
422 real :: min_depth, mask_depth
424 call calltree_enter(trim(mdl)//
"(), MOM_shared_initialization.F90")
426 m_to_z = 1.0 ;
if (
present(us)) m_to_z = us%m_to_Z
428 call get_param(param_file, mdl,
"MINIMUM_DEPTH", min_depth, &
429 "If MASKING_DEPTH is unspecified, then anything shallower than "//&
430 "MINIMUM_DEPTH is assumed to be land and all fluxes are masked out. "//&
431 "If MASKING_DEPTH is specified, then all depths shallower than "//&
432 "MINIMUM_DEPTH but deeper than MASKING_DEPTH are rounded to MINIMUM_DEPTH.", &
433 units=
"m", default=0.0, scale=m_to_z)
434 call get_param(param_file, mdl,
"MASKING_DEPTH", mask_depth, &
435 "The depth below which to mask the ocean as land.", &
436 units=
"m", default=-9999.0, scale=m_to_z, do_not_log=.true.)
439 if (mask_depth < -9990.*m_to_z)
then
440 do j=g%jsd,g%jed ;
do i=g%isd,g%ied
441 d(i,j) = min( max( d(i,j), 0.5*min_depth ), max_depth )
444 do j=g%jsd,g%jed ;
do i=g%isd,g%ied
446 d(i,j) = min( max( d(i,j), min_depth ), max_depth )
453 call calltree_leave(trim(mdl)//
'()')
454 end subroutine limit_topography
459 subroutine set_rotation_planetary(f, G, param_file, US)
461 real,
dimension(G%IsdB:G%IedB,G%JsdB:G%JedB), &
467 character(len=30) :: mdl =
"set_rotation_planetary"
473 call calltree_enter(trim(mdl)//
"(), MOM_shared_initialization.F90")
475 t_to_s = 1.0 ;
if (
present(us)) t_to_s = us%T_to_s
477 call get_param(param_file,
"set_rotation_planetary",
"OMEGA", omega, &
478 "The rotation rate of the earth.", units=
"s-1", &
479 default=7.2921e-5, scale=t_to_s)
482 do i=g%IsdB,g%IedB ;
do j=g%JsdB,g%JedB
483 f(i,j) = ( 2.0 * omega ) * sin( ( pi * g%geoLatBu(i,j) ) / 180.)
486 call calltree_leave(trim(mdl)//
'()')
487 end subroutine set_rotation_planetary
492 subroutine set_rotation_beta_plane(f, G, param_file, US)
494 real,
dimension(G%IsdB:G%IedB,G%JsdB:G%JedB), &
503 real :: y_scl, rad_earth
506 character(len=40) :: mdl =
"set_rotation_beta_plane"
507 character(len=200) :: axis_units
509 call calltree_enter(trim(mdl)//
"(), MOM_shared_initialization.F90")
511 t_to_s = 1.0 ;
if (
present(us)) t_to_s = us%T_to_s
513 call get_param(param_file, mdl,
"F_0", f_0, &
514 "The reference value of the Coriolis parameter with the "//&
515 "betaplane option.", units=
"s-1", default=0.0, scale=t_to_s)
516 call get_param(param_file, mdl,
"BETA", beta, &
517 "The northward gradient of the Coriolis parameter with "//&
518 "the betaplane option.", units=
"m-1 s-1", default=0.0, scale=t_to_s)
519 call get_param(param_file, mdl,
"AXIS_UNITS", axis_units, default=
"degrees")
522 select case (axis_units(1:1))
524 call get_param(param_file, mdl,
"RAD_EARTH", rad_earth, &
525 "The radius of the Earth.", units=
"m", default=6.378e6)
527 case (
"k"); y_scl = 1.e3
528 case (
"m"); y_scl = 1.
529 case (
"c"); y_scl = 1.e-2
530 case default ;
call mom_error(fatal, &
531 " set_rotation_beta_plane: unknown AXIS_UNITS = "//trim(axis_units))
534 do i=g%IsdB,g%IedB ;
do j=g%JsdB,g%JedB
535 f(i,j) = f_0 + beta * ( g%geoLatBu(i,j) * y_scl )
538 call calltree_leave(trim(mdl)//
'()')
539 end subroutine set_rotation_beta_plane
543 subroutine initialize_grid_rotation_angle(G, PF)
548 real :: angle, lon_scale
552 character(len=40) :: mdl =
"initialize_grid_rotation_angle"
554 integer :: i, j, m, n
556 call get_param(pf, mdl,
"GRID_ROTATION_ANGLE_BUGS", use_bugs, &
557 "If true, use an older algorithm to calculate the sine and "//&
558 "cosines needed rotate between grid-oriented directions and "//&
559 "true north and east. Differences arise at the tripolar fold.", &
563 do j=g%jsc,g%jec ;
do i=g%isc,g%iec
564 lon_scale = cos((g%geoLatBu(i-1,j-1) + g%geoLatBu(i,j-1 ) + &
565 g%geoLatBu(i-1,j) + g%geoLatBu(i,j)) * atan(1.0)/180)
566 angle = atan2((g%geoLonBu(i-1,j) + g%geoLonBu(i,j) - &
567 g%geoLonBu(i-1,j-1) - g%geoLonBu(i,j-1))*lon_scale, &
568 g%geoLatBu(i-1,j) + g%geoLatBu(i,j) - &
569 g%geoLatBu(i-1,j-1) - g%geoLatBu(i,j-1) )
570 g%sin_rot(i,j) = sin(angle)
571 g%cos_rot(i,j) = cos(angle)
578 pi_720deg = atan(1.0) / 180.0
579 len_lon = 360.0 ;
if (g%len_lon > 0.0) len_lon = g%len_lon
580 do j=g%jsc,g%jec ;
do i=g%isc,g%iec
582 lonb(m,n) = modulo_around_point(g%geoLonBu(i+m-2,j+n-2), g%geoLonT(i,j), len_lon)
584 lon_scale = cos(pi_720deg*((g%geoLatBu(i-1,j-1) + g%geoLatBu(i,j)) + &
585 (g%geoLatBu(i,j-1) + g%geoLatBu(i-1,j)) ) )
586 angle = atan2(lon_scale*((lonb(1,2) - lonb(2,1)) + (lonb(2,2) - lonb(1,1))), &
587 (g%geoLatBu(i-1,j) - g%geoLatBu(i,j-1)) + &
588 (g%geoLatBu(i,j) - g%geoLatBu(i-1,j-1)) )
589 g%sin_rot(i,j) = sin(angle)
590 g%cos_rot(i,j) = cos(angle)
593 call pass_vector(g%cos_rot, g%sin_rot, g%Domain, stagger=agrid)
596 end subroutine initialize_grid_rotation_angle
601 function modulo_around_point(x, xc, Lx)
result(x_mod)
602 real,
intent(in) :: x
603 real,
intent(in) :: xc
604 real,
intent(in) :: lx
608 x_mod = modulo(x - (xc - 0.5*lx), lx) + (xc - 0.5*lx)
612 end function modulo_around_point
617 subroutine reset_face_lengths_named(G, param_file, name, US)
620 character(len=*),
intent(in) :: name
625 character(len=256) :: mesg
626 real :: dx_2 = -1.0, dy_2 = -1.0
628 integer :: option = -1
629 integer :: i, j, isd, ied, jsd, jed, isdb, iedb, jsdb, jedb
630 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
631 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
632 pi_180 = (4.0*atan(1.0))/180.0
634 select case ( trim(name) )
635 case (
"global_1deg") ; option = 1 ; dx_2 = 0.5*1.0
636 case default ;
call mom_error(fatal,
"reset_face_lengths_named: "//&
637 "Unrecognized channel configuration name "//trim(name))
641 do j=jsd,jed ;
do i=isdb,iedb
642 dy_2 = dx_2 * g%dyCu(i,j)*g%IdxCu(i,j) * cos(pi_180 * g%geoLatCu(i,j))
644 if ((abs(g%geoLatCu(i,j)-35.5) < dy_2) .and. (g%geoLonCu(i,j) < -4.5) .and. &
645 (g%geoLonCu(i,j) > -6.5)) &
646 g%dy_Cu(i,j) = g%mask2dCu(i,j)*12000.0
648 if ((abs(g%geoLatCu(i,j)-12.5) < dy_2) .and. (abs(g%geoLonCu(i,j)-43.0) < dx_2)) &
649 g%dy_Cu(i,j) = g%mask2dCu(i,j)*10000.0
651 if ((abs(g%geoLatCu(i,j)-40.5) < dy_2) .and. (abs(g%geoLonCu(i,j)-26.0) < dx_2)) &
652 g%dy_Cu(i,j) = g%mask2dCu(i,j)*5000.0
654 if ((abs(g%geoLatCu(i,j)-41.5) < dy_2) .and. (abs(g%geoLonCu(i,j)+220.0) < dx_2)) &
655 g%dy_Cu(i,j) = g%mask2dCu(i,j)*35000.0
657 if ((abs(g%geoLatCu(i,j)-45.5) < dy_2) .and. (abs(g%geoLonCu(i,j)+217.5) < 0.9)) &
658 g%dy_Cu(i,j) = g%mask2dCu(i,j)*15000.0
662 if ((abs(g%geoLatCu(i,j)-80.84) < 0.2) .and. (abs(g%geoLonCu(i,j)+64.9) < 0.8)) &
663 g%dy_Cu(i,j) = g%mask2dCu(i,j)*38000.0
667 do j=jsdb,jedb ;
do i=isd,ied
668 dy_2 = dx_2 * g%dyCv(i,j)*g%IdxCv(i,j) * cos(pi_180 * g%geoLatCv(i,j))
669 if ((abs(g%geoLatCv(i,j)-41.0) < dy_2) .and. (abs(g%geoLonCv(i,j)-28.5) < dx_2)) &
670 g%dx_Cv(i,j) = g%mask2dCv(i,j)*2500.0
672 if ((abs(g%geoLatCv(i,j)-13.0) < dy_2) .and. (abs(g%geoLonCv(i,j)-42.5) < dx_2)) &
673 g%dx_Cv(i,j) = g%mask2dCv(i,j)*10000.0
675 if ((abs(g%geoLatCv(i,j)+2.8) < 0.8) .and. (abs(g%geoLonCv(i,j)+241.5) < dx_2)) &
676 g%dx_Cv(i,j) = g%mask2dCv(i,j)*40000.0
678 if ((abs(g%geoLatCv(i,j)-0.56) < 0.5) .and. (abs(g%geoLonCv(i,j)+240.5) < dx_2)) &
679 g%dx_Cv(i,j) = g%mask2dCv(i,j)*80000.0
681 if ((abs(g%geoLatCv(i,j)-0.19) < 0.5) .and. (abs(g%geoLonCv(i,j)+230.5) < dx_2)) &
682 g%dx_Cv(i,j) = g%mask2dCv(i,j)*25000.0
684 if ((abs(g%geoLatCv(i,j)-0.19) < 0.5) .and. (abs(g%geoLonCv(i,j)+229.5) < dx_2)) &
685 g%dx_Cv(i,j) = g%mask2dCv(i,j)*25000.0
687 if ((abs(g%geoLatCv(i,j)-0.0) < 0.25) .and. (abs(g%geoLonCv(i,j)+228.5) < dx_2)) &
688 g%dx_Cv(i,j) = g%mask2dCv(i,j)*25000.0
690 if ((abs(g%geoLatCv(i,j)+8.5) < 0.5) .and. (abs(g%geoLonCv(i,j)+244.5) < dx_2)) &
691 g%dx_Cv(i,j) = g%mask2dCv(i,j)*20000.0
693 if ((abs(g%geoLatCv(i,j)+8.5) < 0.5) .and. (abs(g%geoLonCv(i,j)+235.5) < dx_2)) &
694 g%dx_Cv(i,j) = g%mask2dCv(i,j)*20000.0
696 if ((abs(g%geoLatCv(i,j)-52.5) < dy_2) .and. (abs(g%geoLonCv(i,j)+218.5) < dx_2)) &
697 g%dx_Cv(i,j) = g%mask2dCv(i,j)*2500.0
700 if ((abs(g%geoLatCv(i,j)-76.8) < 0.06) .and. (abs(g%geoLonCv(i,j)+88.7) < dx_2)) &
701 g%dx_Cv(i,j) = g%mask2dCv(i,j)*8400.0
708 do j=jsd,jed ;
do i=isdb,iedb
709 if (g%dy_Cu(i,j) > g%dyCu(i,j))
then
710 write(mesg,
'("dy_Cu of ",ES11.4," exceeds unrestricted width of ",ES11.4,&
711 &" by ",ES11.4," at lon/lat of ", ES11.4, ES11.4)') &
712 g%dy_Cu(i,j), g%dyCu(i,j), g%dy_Cu(i,j)-g%dyCu(i,j), &
713 g%geoLonCu(i,j), g%geoLatCu(i,j)
714 call mom_error(fatal,
"reset_face_lengths_named "//mesg)
716 g%areaCu(i,j) = g%dxCu(i,j)*g%dy_Cu(i,j)
718 if (g%areaCu(i,j) > 0.0) g%IareaCu(i,j) = g%mask2dCu(i,j) / g%areaCu(i,j)
721 do j=jsdb,jedb ;
do i=isd,ied
722 if (g%dx_Cv(i,j) > g%dxCv(i,j))
then
723 write(mesg,
'("dx_Cv of ",ES11.4," exceeds unrestricted width of ",ES11.4,&
724 &" by ",ES11.4, " at lon/lat of ", ES11.4, ES11.4)') &
725 g%dx_Cv(i,j), g%dxCv(i,j), g%dx_Cv(i,j)-g%dxCv(i,j), &
726 g%geoLonCv(i,j), g%geoLatCv(i,j)
728 call mom_error(fatal,
"reset_face_lengths_named "//mesg)
730 g%areaCv(i,j) = g%dyCv(i,j)*g%dx_Cv(i,j)
732 if (g%areaCv(i,j) > 0.0) g%IareaCv(i,j) = g%mask2dCv(i,j) / g%areaCv(i,j)
735 end subroutine reset_face_lengths_named
741 subroutine reset_face_lengths_file(G, param_file, US)
747 character(len=40) :: mdl =
"reset_face_lengths_file"
748 character(len=256) :: mesg
749 character(len=200) :: filename, chan_file, inputdir
750 integer :: i, j, isd, ied, jsd, jed, isdb, iedb, jsdb, jedb
751 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
752 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
755 call calltree_enter(trim(mdl)//
"(), MOM_shared_initialization.F90")
757 call get_param(param_file, mdl,
"CHANNEL_WIDTH_FILE", chan_file, &
758 "The file from which the list of narrowed channels is read.", &
759 default=
"ocean_geometry.nc")
760 call get_param(param_file, mdl,
"INPUTDIR", inputdir, default=
".")
761 inputdir = slasher(inputdir)
762 filename = trim(inputdir)//trim(chan_file)
763 call log_param(param_file, mdl,
"INPUTDIR/CHANNEL_WIDTH_FILE", filename)
765 if (is_root_pe())
then ;
if (.not.
file_exists(filename)) &
766 call mom_error(fatal,
" reset_face_lengths_file: Unable to open "//&
770 call mom_read_vector(filename,
"dyCuo",
"dxCvo", g%dy_Cu, g%dx_Cv, g%Domain)
771 call pass_vector(g%dy_Cu, g%dx_Cv, g%Domain, to_all+scalar_pair, cgrid_ne)
773 do j=jsd,jed ;
do i=isdb,iedb
774 if (g%dy_Cu(i,j) > g%dyCu(i,j))
then
775 write(mesg,
'("dy_Cu of ",ES11.4," exceeds unrestricted width of ",ES11.4,&
776 &" by ",ES11.4," at lon/lat of ", ES11.4, ES11.4)') &
777 g%dy_Cu(i,j), g%dyCu(i,j), g%dy_Cu(i,j)-g%dyCu(i,j), &
778 g%geoLonCu(i,j), g%geoLatCu(i,j)
779 call mom_error(fatal,
"reset_face_lengths_file "//mesg)
781 g%areaCu(i,j) = g%dxCu(i,j)*g%dy_Cu(i,j)
783 if (g%areaCu(i,j) > 0.0) g%IareaCu(i,j) = g%mask2dCu(i,j) / g%areaCu(i,j)
786 do j=jsdb,jedb ;
do i=isd,ied
787 if (g%dx_Cv(i,j) > g%dxCv(i,j))
then
788 write(mesg,
'("dx_Cv of ",ES11.4," exceeds unrestricted width of ",ES11.4,&
789 &" by ",ES11.4, " at lon/lat of ", ES11.4, ES11.4)') &
790 g%dx_Cv(i,j), g%dxCv(i,j), g%dx_Cv(i,j)-g%dxCv(i,j), &
791 g%geoLonCv(i,j), g%geoLatCv(i,j)
793 call mom_error(fatal,
"reset_face_lengths_file "//mesg)
795 g%areaCv(i,j) = g%dyCv(i,j)*g%dx_Cv(i,j)
797 if (g%areaCv(i,j) > 0.0) g%IareaCv(i,j) = g%mask2dCv(i,j) / g%areaCv(i,j)
800 call calltree_leave(trim(mdl)//
'()')
801 end subroutine reset_face_lengths_file
807 subroutine reset_face_lengths_list(G, param_file, US)
813 character(len=120),
pointer,
dimension(:) :: lines => null()
814 character(len=120) :: line
815 character(len=200) :: filename, chan_file, inputdir, mesg
816 character(len=40) :: mdl =
"reset_face_lengths_list"
817 real,
pointer,
dimension(:,:) :: &
818 u_lat => null(), u_lon => null(), v_lat => null(), v_lon => null()
819 real,
pointer,
dimension(:) :: &
820 u_width => null(), v_width => null()
827 logical :: found_u, found_v
828 logical :: unit_in_use
829 integer :: ios, iounit, isu, isv
830 integer :: last, num_lines, nl_read, ln, npt, u_pt, v_pt
831 integer :: i, j, isd, ied, jsd, jed, isdb, iedb, jsdb, jedb
832 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
833 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
835 call calltree_enter(trim(mdl)//
"(), MOM_shared_initialization.F90")
837 call get_param(param_file, mdl,
"CHANNEL_LIST_FILE", chan_file, &
838 "The file from which the list of narrowed channels is read.", &
839 default=
"MOM_channel_list")
840 call get_param(param_file, mdl,
"INPUTDIR", inputdir, default=
".")
841 inputdir = slasher(inputdir)
842 filename = trim(inputdir)//trim(chan_file)
843 call log_param(param_file, mdl,
"INPUTDIR/CHANNEL_LIST_FILE", filename)
844 call get_param(param_file, mdl,
"CHANNEL_LIST_360_LON_CHECK", check_360, &
845 "If true, the channel configuration list works for any "//&
846 "longitudes in the range of -360 to 360.", default=.true.)
848 if (is_root_pe())
then
850 if (.not.
file_exists(filename))
call mom_error(fatal, &
851 " reset_face_lengths_list: Unable to open "//trim(filename))
855 INQUIRE(iounit,opened=unit_in_use) ;
if (.not.unit_in_use)
exit
857 if (iounit >= 512)
call mom_error(fatal, &
858 "reset_face_lengths_list: No unused file unit could be found.")
861 open(iounit, file=trim(filename), access=
'SEQUENTIAL', &
862 form=
'FORMATTED', action=
'READ', position=
'REWIND', iostat=ios)
863 if (ios /= 0)
call mom_error(fatal, &
864 "reset_face_lengths_list: Error opening "//trim(filename))
867 call read_face_length_list(iounit, filename, num_lines, lines)
870 len_lon = 360.0 ;
if (g%len_lon > 0.0) len_lon = g%len_lon
871 len_lat = 180.0 ;
if (g%len_lat > 0.0) len_lat = g%len_lat
873 call broadcast(num_lines, root_pe())
875 if (num_lines > 0)
then
876 allocate (lines(num_lines))
877 if (num_lines > 0)
then
878 allocate(u_lat(2,num_lines)) ; u_lat(:,:) = -1e34
879 allocate(u_lon(2,num_lines)) ; u_lon(:,:) = -1e34
880 allocate(u_width(num_lines)) ; u_width(:) = -1e34
882 allocate(v_lat(2,num_lines)) ; v_lat(:,:) = -1e34
883 allocate(v_lon(2,num_lines)) ; v_lon(:,:) = -1e34
884 allocate(v_width(num_lines)) ; v_width(:) = -1e34
888 if (is_root_pe())
then
889 call read_face_length_list(iounit, filename, nl_read, lines)
890 if (nl_read /= num_lines) &
891 call mom_error(fatal,
'reset_face_lengths_list : Found different '// &
892 'number of valid lines on second reading of '//trim(filename))
893 close(iounit) ; iounit = -1
897 call broadcast(lines, 120, root_pe())
903 found_u = .false.; found_v = .false.
904 isu = index(uppercase(line),
"U_WIDTH" );
if (isu > 0) found_u = .true.
905 isv = index(uppercase(line),
"V_WIDTH" );
if (isv > 0) found_v = .true.
910 read(line(isu+8:),*) u_lon(1:2,u_pt), u_lat(1:2,u_pt), u_width(u_pt)
911 if (is_root_pe())
then
913 if ((abs(u_lon(1,u_pt)) > len_lon) .or. (abs(u_lon(2,u_pt)) > len_lon)) &
914 call mom_error(warning,
"reset_face_lengths_list : Out-of-bounds "//&
915 "u-longitude found when reading line "//trim(line)//
" from file "//&
917 if ((abs(u_lat(1,u_pt)) > len_lat) .or. (abs(u_lat(2,u_pt)) > len_lat)) &
918 call mom_error(warning,
"reset_face_lengths_list : Out-of-bounds "//&
919 "u-latitude found when reading line "//trim(line)//
" from file "//&
922 if (u_lat(1,u_pt) > u_lat(2,u_pt)) &
923 call mom_error(warning,
"reset_face_lengths_list : Out-of-order "//&
924 "u-face latitudes found when reading line "//trim(line)//
" from file "//&
926 if (u_lon(1,u_pt) > u_lon(2,u_pt)) &
927 call mom_error(warning,
"reset_face_lengths_list : Out-of-order "//&
928 "u-face longitudes found when reading line "//trim(line)//
" from file "//&
930 if (u_width(u_pt) < 0.0) &
931 call mom_error(warning,
"reset_face_lengths_list : Negative "//&
932 "u-width found when reading line "//trim(line)//
" from file "//&
935 elseif (found_v)
then
937 read(line(isv+8:),*) v_lon(1:2,v_pt), v_lat(1:2,v_pt), v_width(v_pt)
938 if (is_root_pe())
then
940 if ((abs(v_lon(1,v_pt)) > len_lon) .or. (abs(v_lon(2,v_pt)) > len_lon)) &
941 call mom_error(warning,
"reset_face_lengths_list : Out-of-bounds "//&
942 "v-longitude found when reading line "//trim(line)//
" from file "//&
944 if ((abs(v_lat(1,v_pt)) > len_lat) .or. (abs(v_lat(2,v_pt)) > len_lat)) &
945 call mom_error(warning,
"reset_face_lengths_list : Out-of-bounds "//&
946 "v-latitude found when reading line "//trim(line)//
" from file "//&
949 if (v_lat(1,v_pt) > v_lat(2,v_pt)) &
950 call mom_error(warning,
"reset_face_lengths_list : Out-of-order "//&
951 "v-face latitudes found when reading line "//trim(line)//
" from file "//&
953 if (v_lon(1,v_pt) > v_lon(2,v_pt)) &
954 call mom_error(warning,
"reset_face_lengths_list : Out-of-order "//&
955 "v-face longitudes found when reading line "//trim(line)//
" from file "//&
957 if (v_width(v_pt) < 0.0) &
958 call mom_error(warning,
"reset_face_lengths_list : Negative "//&
959 "v-width found when reading line "//trim(line)//
" from file "//&
968 do j=jsd,jed ;
do i=isdb,iedb
969 lat = g%geoLatCu(i,j) ; lon = g%geoLonCu(i,j)
970 if (check_360)
then ; lon_p = lon+len_lon ; lon_m = lon-len_lon
971 else ; lon_p = lon ; lon_m = lon ;
endif
974 if (((lat >= u_lat(1,npt)) .and. (lat <= u_lat(2,npt))) .and. &
975 (((lon >= u_lon(1,npt)) .and. (lon <= u_lon(2,npt))) .or. &
976 ((lon_p >= u_lon(1,npt)) .and. (lon_p <= u_lon(2,npt))) .or. &
977 ((lon_m >= u_lon(1,npt)) .and. (lon_m <= u_lon(2,npt)))) )
then
979 g%dy_Cu(i,j) = g%mask2dCu(i,j) * min(g%dyCu(i,j), max(u_width(npt), 0.0))
980 if (j>=g%jsc .and. j<=g%jec .and. i>=g%isc .and. i<=g%iec)
then
981 if ( g%mask2dCu(i,j) == 0.0 )
then
982 write(*,
'(A,2F8.2,A,4F8.2,A)')
"read_face_lengths_list : G%mask2dCu=0 at ",lat,lon,
" (",&
983 u_lat(1,npt), u_lat(2,npt), u_lon(1,npt), u_lon(2,npt),
") so grid metric is unmodified."
985 write(*,
'(A,2F8.2,A,4F8.2,A5,F9.2,A1)') &
986 "read_face_lengths_list : Modifying dy_Cu gridpoint at ",lat,lon,
" (",&
987 u_lat(1,npt), u_lat(2,npt), u_lon(1,npt), u_lon(2,npt),
") to ",g%dy_Cu(i,j),
"m"
993 g%areaCu(i,j) = g%dxCu(i,j)*g%dy_Cu(i,j)
995 if (g%areaCu(i,j) > 0.0) g%IareaCu(i,j) = g%mask2dCu(i,j) / g%areaCu(i,j)
998 do j=jsdb,jedb ;
do i=isd,ied
999 lat = g%geoLatCv(i,j) ; lon = g%geoLonCv(i,j)
1000 if (check_360)
then ; lon_p = lon+len_lon ; lon_m = lon-len_lon
1001 else ; lon_p = lon ; lon_m = lon ;
endif
1004 if (((lat >= v_lat(1,npt)) .and. (lat <= v_lat(2,npt))) .and. &
1005 (((lon >= v_lon(1,npt)) .and. (lon <= v_lon(2,npt))) .or. &
1006 ((lon_p >= v_lon(1,npt)) .and. (lon_p <= v_lon(2,npt))) .or. &
1007 ((lon_m >= v_lon(1,npt)) .and. (lon_m <= v_lon(2,npt)))) )
then
1008 g%dx_Cv(i,j) = g%mask2dCv(i,j) * min(g%dxCv(i,j), max(v_width(npt), 0.0))
1009 if (i>=g%isc .and. i<=g%iec .and. j>=g%jsc .and. j<=g%jec)
then
1010 if ( g%mask2dCv(i,j) == 0.0 )
then
1011 write(*,
'(A,2F8.2,A,4F8.2,A)')
"read_face_lengths_list : G%mask2dCv=0 at ",lat,lon,
" (",&
1012 v_lat(1,npt), v_lat(2,npt), v_lon(1,npt), v_lon(2,npt),
") so grid metric is unmodified."
1014 write(*,
'(A,2F8.2,A,4F8.2,A5,F9.2,A1)') &
1015 "read_face_lengths_list : Modifying dx_Cv gridpoint at ",lat,lon,
" (",&
1016 v_lat(1,npt), v_lat(2,npt), v_lon(1,npt), v_lon(2,npt),
") to ",g%dx_Cv(i,j),
"m"
1022 g%areaCv(i,j) = g%dyCv(i,j)*g%dx_Cv(i,j)
1023 g%IareaCv(i,j) = 0.0
1024 if (g%areaCv(i,j) > 0.0) g%IareaCv(i,j) = g%mask2dCv(i,j) / g%areaCv(i,j)
1027 if (num_lines > 0)
then
1028 deallocate(u_lat) ;
deallocate(u_lon) ;
deallocate(u_width)
1029 deallocate(v_lat) ;
deallocate(v_lon) ;
deallocate(v_width)
1032 call calltree_leave(trim(mdl)//
'()')
1033 end subroutine reset_face_lengths_list
1038 subroutine read_face_length_list(iounit, filename, num_lines, lines)
1039 integer,
intent(in) :: iounit
1040 character(len=*),
intent(in) :: filename
1041 integer,
intent(out) :: num_lines
1042 character(len=120),
dimension(:),
pointer :: lines
1046 character(len=120) :: line, line_up
1047 logical :: found_u, found_v
1048 integer :: isu, isv, icom, verbose
1053 if (iounit <= 0)
return
1056 read(iounit,
'(a)',
end=8, err=9) line
1057 last = len_trim(line)
1059 icom = index(line(:last),
"!") ;
if (icom > 0) last = icom-1
1060 icom = index(line(:last),
"/*") ;
if (icom > 0) last = icom-1
1064 line_up = uppercase(line)
1065 found_u = .false.; found_v = .false.
1066 isu = index(line_up(:last),
"U_WIDTH" );
if (isu > 0) found_u = .true.
1067 isv = index(line_up(:last),
"V_WIDTH" );
if (isv > 0) found_v = .true.
1069 if (found_u .and. found_v)
call mom_error(fatal, &
1070 "read_face_length_list : both U_WIDTH and V_WIDTH found when "//&
1071 "reading the line "//trim(line(:last))//
" in file "//trim(filename))
1072 if (found_u .or. found_v)
then
1073 num_lines = num_lines + 1
1074 if (
associated(lines))
then
1075 lines(num_lines) = line(1:last)
1083 9
call mom_error(fatal,
"read_face_length_list : "//&
1084 "Error while reading file "//trim(filename))
1086 end subroutine read_face_length_list
1092 subroutine set_velocity_depth_max(G)
1093 type(dyn_horgrid_type),
intent(inout) :: g
1098 do i=g%isd,g%ied-1 ;
do j=g%jsd,g%jed
1099 g%Dblock_u(i,j) = g%mask2dCu(i,j) * max(g%bathyT(i,j), g%bathyT(i+1,j))
1100 g%Dopen_u(i,j) = g%Dblock_u(i,j)
1102 do i=g%isd,g%ied ;
do j=g%jsd,g%jed-1
1103 g%Dblock_v(i,j) = g%mask2dCv(i,j) * max(g%bathyT(i,j), g%bathyT(i,j+1))
1104 g%Dopen_v(i,j) = g%Dblock_v(i,j)
1106 end subroutine set_velocity_depth_max
1112 subroutine set_velocity_depth_min(G)
1113 type(dyn_horgrid_type),
intent(inout) :: g
1118 do i=g%isd,g%ied-1 ;
do j=g%jsd,g%jed
1119 g%Dblock_u(i,j) = g%mask2dCu(i,j) * min(g%bathyT(i,j), g%bathyT(i+1,j))
1120 g%Dopen_u(i,j) = g%Dblock_u(i,j)
1122 do i=g%isd,g%ied ;
do j=g%jsd,g%jed-1
1123 g%Dblock_v(i,j) = g%mask2dCv(i,j) * min(g%bathyT(i,j), g%bathyT(i,j+1))
1124 g%Dopen_v(i,j) = g%Dblock_v(i,j)
1126 end subroutine set_velocity_depth_min
1132 subroutine compute_global_grid_integrals(G)
1133 type(dyn_horgrid_type),
intent(inout) :: g
1136 real,
dimension(G%isc:G%iec, G%jsc:G%jec) :: tmpforsumming
1139 tmpforsumming(:,:) = 0.
1140 g%areaT_global = 0.0 ; g%IareaT_global = 0.0
1141 do j=g%jsc,g%jec ;
do i=g%isc,g%iec
1142 tmpforsumming(i,j) = g%areaT(i,j) * g%mask2dT(i,j)
1144 g%areaT_global = reproducing_sum(tmpforsumming)
1146 if (g%areaT_global == 0.0) &
1147 call mom_error(fatal,
"compute_global_grid_integrals: "//&
1148 "zero ocean area (check topography?)")
1150 g%IareaT_global = 1. / g%areaT_global
1151 end subroutine compute_global_grid_integrals
1157 subroutine write_ocean_geometry_file(G, param_file, directory, geom_file, US)
1158 type(dyn_horgrid_type),
intent(inout) :: g
1160 character(len=*),
intent(in) :: directory
1161 character(len=*),
optional,
intent(in) :: geom_file
1163 type(unit_scale_type),
optional,
intent(in) :: us
1166 character(len=240) :: filepath
1167 character(len=40) :: mdl =
"write_ocean_geometry_file"
1168 integer,
parameter :: nflds=23
1169 type(vardesc) :: vars(nflds)
1170 type(fieldtype) :: fields(nflds)
1171 real :: z_to_m_scale
1172 real :: s_to_t_scale
1174 integer :: file_threading
1175 integer :: nflds_used
1176 integer :: i, j, is, ie, js, je, isq, ieq, jsq, jeq
1177 integer :: isd, ied, jsd, jed, isdb, iedb, jsdb, jedb
1178 logical :: multiple_files
1179 real,
dimension(G%isd :G%ied ,G%jsd :G%jed ) :: out_h
1180 real,
dimension(G%IsdB:G%IedB,G%JsdB:G%JedB) :: out_q
1181 real,
dimension(G%IsdB:G%IedB,G%jsd :G%jed ) :: out_u
1182 real,
dimension(G%isd :G%ied ,G%JsdB:G%JedB) :: out_v
1184 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
1185 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
1186 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
1187 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
1189 z_to_m_scale = 1.0 ;
if (
present(us)) z_to_m_scale = us%Z_to_m
1190 s_to_t_scale = 1.0 ;
if (
present(us)) s_to_t_scale = us%s_to_T
1203 vars(1) = var_desc(
"geolatb",
"degree",
"latitude at corner (Bu) points",
'q',
'1',
'1')
1204 vars(2) = var_desc(
"geolonb",
"degree",
"longitude at corner (Bu) points",
'q',
'1',
'1')
1205 vars(3) = var_desc(
"geolat",
"degree",
"latitude at tracer (T) points",
'h',
'1',
'1')
1206 vars(4) = var_desc(
"geolon",
"degree",
"longitude at tracer (T) points",
'h',
'1',
'1')
1207 vars(5) = var_desc(
"D",
"meter",
"Basin Depth",
'h',
'1',
'1')
1208 vars(6) = var_desc(
"f",
"s-1",
"Coriolis Parameter",
'q',
'1',
'1')
1209 vars(7) = var_desc(
"dxCv",
"m",
"Zonal grid spacing at v points",
'v',
'1',
'1')
1210 vars(8) = var_desc(
"dyCu",
"m",
"Meridional grid spacing at u points",
'u',
'1',
'1')
1211 vars(9) = var_desc(
"dxCu",
"m",
"Zonal grid spacing at u points",
'u',
'1',
'1')
1212 vars(10)= var_desc(
"dyCv",
"m",
"Meridional grid spacing at v points",
'v',
'1',
'1')
1213 vars(11)= var_desc(
"dxT",
"m",
"Zonal grid spacing at h points",
'h',
'1',
'1')
1214 vars(12)= var_desc(
"dyT",
"m",
"Meridional grid spacing at h points",
'h',
'1',
'1')
1215 vars(13)= var_desc(
"dxBu",
"m",
"Zonal grid spacing at q points",
'q',
'1',
'1')
1216 vars(14)= var_desc(
"dyBu",
"m",
"Meridional grid spacing at q points",
'q',
'1',
'1')
1217 vars(15)= var_desc(
"Ah",
"m2",
"Area of h cells",
'h',
'1',
'1')
1218 vars(16)= var_desc(
"Aq",
"m2",
"Area of q cells",
'q',
'1',
'1')
1220 vars(17)= var_desc(
"dxCvo",
"m",
"Open zonal grid spacing at v points",
'v',
'1',
'1')
1221 vars(18)= var_desc(
"dyCuo",
"m",
"Open meridional grid spacing at u points",
'u',
'1',
'1')
1222 vars(19)= var_desc(
"wet",
"nondim",
"land or ocean?",
'h',
'1',
'1')
1224 vars(20) = var_desc(
"Dblock_u",
"m",
"Blocked depth at u points",
'u',
'1',
'1')
1225 vars(21) = var_desc(
"Dopen_u",
"m",
"Open depth at u points",
'u',
'1',
'1')
1226 vars(22) = var_desc(
"Dblock_v",
"m",
"Blocked depth at v points",
'v',
'1',
'1')
1227 vars(23) = var_desc(
"Dopen_v",
"m",
"Open depth at v points",
'v',
'1',
'1')
1230 nflds_used = 19 ;
if (g%bathymetry_at_vel) nflds_used = 23
1232 if (
present(geom_file))
then
1233 filepath = trim(directory) // trim(geom_file)
1235 filepath = trim(directory) //
"ocean_geometry"
1243 call get_param(param_file, mdl,
"PARALLEL_RESTARTFILES", multiple_files, &
1244 "If true, each processor writes its own restart file, "//&
1245 "otherwise a single restart file is generated", &
1247 file_threading = single_file
1248 if (multiple_files) file_threading = multiple
1250 call create_file(unit, trim(filepath), vars, nflds_used, fields, &
1251 file_threading, dg=g)
1253 do j=jsq,jeq;
do i=isq,ieq; out_q(i,j) = g%geoLatBu(i,j);
enddo ;
enddo
1254 call write_field(unit, fields(1), g%Domain%mpp_domain, out_q)
1255 do j=jsq,jeq;
do i=isq,ieq; out_q(i,j) = g%geoLonBu(i,j);
enddo ;
enddo
1256 call write_field(unit, fields(2), g%Domain%mpp_domain, out_q)
1257 call write_field(unit, fields(3), g%Domain%mpp_domain, g%geoLatT)
1258 call write_field(unit, fields(4), g%Domain%mpp_domain, g%geoLonT)
1260 do j=js,je ;
do i=is,ie ; out_h(i,j) = z_to_m_scale*g%bathyT(i,j) ;
enddo ;
enddo
1261 call write_field(unit, fields(5), g%Domain%mpp_domain, out_h)
1262 do j=jsq,jeq ;
do i=isq,ieq ; out_q(i,j) = s_to_t_scale*g%CoriolisBu(i,j) ;
enddo ;
enddo
1263 call write_field(unit, fields(6), g%Domain%mpp_domain, out_q)
1268 do j=jsq,jeq ;
do i=is,ie ; out_v(i,j) = g%dxCv(i,j) ;
enddo ;
enddo
1269 call write_field(unit, fields(7), g%Domain%mpp_domain, out_v)
1270 do j=js,je ;
do i=isq,ieq ; out_u(i,j) = g%dyCu(i,j) ;
enddo ;
enddo
1271 call write_field(unit, fields(8), g%Domain%mpp_domain, out_u)
1273 do j=js,je ;
do i=isq,ieq ; out_u(i,j) = g%dxCu(i,j) ;
enddo ;
enddo
1274 call write_field(unit, fields(9), g%Domain%mpp_domain, out_u)
1275 do j=jsq,jeq ;
do i=is,ie ; out_v(i,j) = g%dyCv(i,j) ;
enddo ;
enddo
1276 call write_field(unit, fields(10), g%Domain%mpp_domain, out_v)
1278 do j=js,je ;
do i=is,ie ; out_h(i,j) = g%dxT(i,j);
enddo ;
enddo
1279 call write_field(unit, fields(11), g%Domain%mpp_domain, out_h)
1280 do j=js,je ;
do i=is,ie ; out_h(i,j) = g%dyT(i,j) ;
enddo ;
enddo
1281 call write_field(unit, fields(12), g%Domain%mpp_domain, out_h)
1283 do j=jsq,jeq ;
do i=isq,ieq ; out_q(i,j) = g%dxBu(i,j) ;
enddo ;
enddo
1284 call write_field(unit, fields(13), g%Domain%mpp_domain, out_q)
1285 do j=jsq,jeq ;
do i=isq,ieq ; out_q(i,j) = g%dyBu(i,j) ;
enddo ;
enddo
1286 call write_field(unit, fields(14), g%Domain%mpp_domain, out_q)
1288 do j=js,je ;
do i=is,ie ; out_h(i,j) = g%areaT(i,j) ;
enddo ;
enddo
1289 call write_field(unit, fields(15), g%Domain%mpp_domain, out_h)
1290 do j=jsq,jeq ;
do i=isq,ieq ; out_q(i,j) = g%areaBu(i,j) ;
enddo ;
enddo
1291 call write_field(unit, fields(16), g%Domain%mpp_domain, out_q)
1293 call write_field(unit, fields(17), g%Domain%mpp_domain, g%dx_Cv)
1294 call write_field(unit, fields(18), g%Domain%mpp_domain, g%dy_Cu)
1295 call write_field(unit, fields(19), g%Domain%mpp_domain, g%mask2dT)
1297 if (g%bathymetry_at_vel)
then
1298 do j=js,je ;
do i=isq,ieq ; out_u(i,j) = z_to_m_scale*g%Dblock_u(i,j) ;
enddo ;
enddo
1299 call write_field(unit, fields(20), g%Domain%mpp_domain, out_u)
1300 do j=js,je ;
do i=isq,ieq ; out_u(i,j) = z_to_m_scale*g%Dopen_u(i,j) ;
enddo ;
enddo
1301 call write_field(unit, fields(21), g%Domain%mpp_domain, out_u)
1302 do j=jsq,jeq ;
do i=is,ie ; out_v(i,j) = z_to_m_scale*g%Dblock_v(i,j) ;
enddo ;
enddo
1303 call write_field(unit, fields(22), g%Domain%mpp_domain, out_v)
1304 do j=jsq,jeq ;
do i=is,ie ; out_v(i,j) = z_to_m_scale*g%Dopen_v(i,j) ;
enddo ;
enddo
1305 call write_field(unit, fields(23), g%Domain%mpp_domain, out_v)
1308 call close_file(unit)
1310 end subroutine write_ocean_geometry_file