22 implicit none ;
private
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.")
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))
83 case default ;
call mom_error(fatal,
"MOM_initialize: "// &
84 "Unrecognized rotation setup "//trim(config))
92 real,
dimension(G%isd:G%ied,G%jsd:G%jed), &
94 real,
dimension(G%isd:G%ied,G%jsd:G%jed), &
102 m_to_l = 1.0 ;
if (
present(us)) m_to_l = us%m_to_L
104 if ((lbound(g%CoriolisBu,1) > g%isc-1) .or. &
105 (lbound(g%CoriolisBu,2) > g%isc-1))
then
107 df_dx(:,:) = 0.0 ; df_dy(:,:) = 0.0
111 do j=g%jsc, g%jec ;
do i=g%isc, g%iec
112 f1 = 0.5*( g%CoriolisBu(i,j) + g%CoriolisBu(i,j-1) )
113 f2 = 0.5*( g%CoriolisBu(i-1,j) + g%CoriolisBu(i-1,j-1) )
114 df_dx(i,j) = g%IdxT(i,j) * ( f1 - f2 )
115 f1 = 0.5*( g%CoriolisBu(i,j) + g%CoriolisBu(i-1,j) )
116 f2 = 0.5*( g%CoriolisBu(i,j-1) + g%CoriolisBu(i-1,j-1) )
117 df_dy(i,j) = g%IdyT(i,j) * ( f1 - f2 )
119 call pass_vector(df_dx, df_dy, g%Domain, stagger=agrid)
126 real,
dimension(G%isd:G%ied,G%jsd:G%jed), &
132 do j=g%jsc, g%jec ;
do i=g%isc, g%iec
142 real,
dimension(G%isd:G%ied,G%jsd:G%jed), &
148 character(len=200) :: filename, topo_file, inputdir
149 character(len=200) :: topo_varname
150 character(len=40) :: mdl =
"initialize_topography_from_file"
152 call calltree_enter(trim(mdl)//
"(), MOM_shared_initialization.F90")
154 m_to_z = 1.0 ;
if (
present(us)) m_to_z = us%m_to_Z
156 call get_param(param_file, mdl,
"INPUTDIR", inputdir, default=
".")
157 inputdir = slasher(inputdir)
158 call get_param(param_file, mdl,
"TOPO_FILE", topo_file, &
159 "The file from which the bathymetry is read.", &
161 call get_param(param_file, mdl,
"TOPO_VARNAME", topo_varname, &
162 "The name of the bathymetry variable in TOPO_FILE.", &
165 filename = trim(inputdir)//trim(topo_file)
166 call log_param(param_file, mdl,
"INPUTDIR/TOPO_FILE", filename)
168 if (.not.
file_exists(filename, g%Domain))
call mom_error(fatal, &
169 " initialize_topography_from_file: Unable to open "//trim(filename))
171 d(:,:) = -9.e30*m_to_z
177 call mom_read_data(filename, trim(topo_varname), d, g%Domain, scale=m_to_z)
187 real,
dimension(G%isd:G%ied,G%jsd:G%jed), &
194 character(len=200) :: topo_edits_file, inputdir
195 character(len=40) :: mdl =
"apply_topography_edits_from_file"
196 integer :: n_edits, n, ashape(5), i, j, ncid, id, ncstatus, iid, jid, zid
197 integer,
dimension(:),
allocatable :: ig, jg
198 real,
dimension(:),
allocatable :: new_depth
200 call calltree_enter(trim(mdl)//
"(), MOM_shared_initialization.F90")
202 m_to_z = 1.0 ;
if (
present(us)) m_to_z = us%m_to_Z
204 call get_param(param_file, mdl,
"INPUTDIR", inputdir, default=
".")
205 inputdir = slasher(inputdir)
206 call get_param(param_file, mdl,
"TOPO_EDITS_FILE", topo_edits_file, &
207 "The file from which to read a list of i,j,z topography overrides.", &
210 if (len_trim(topo_edits_file)==0)
return
212 topo_edits_file = trim(inputdir)//trim(topo_edits_file)
213 if (.not.
file_exists(topo_edits_file, g%Domain))
call mom_error(fatal, &
214 'initialize_topography_from_file: Unable to open '//trim(topo_edits_file))
216 ncstatus = nf90_open(trim(topo_edits_file), nf90_nowrite, ncid)
217 if (ncstatus /= nf90_noerr)
call mom_error(fatal,
'apply_topography_edits_from_file: '//&
218 'Failed to open '//trim(topo_edits_file))
221 ncstatus = nf90_inq_dimid(ncid,
'nEdits', id)
222 if (ncstatus /= nf90_noerr)
call mom_error(fatal,
'apply_topography_edits_from_file: '//&
223 'Failed to inq_dimid nEdits for '//trim(topo_edits_file))
224 ncstatus = nf90_inquire_dimension(ncid, id, len=n_edits)
225 if (ncstatus /= nf90_noerr)
call mom_error(fatal,
'apply_topography_edits_from_file: '//&
226 'Failed to inquire_dimension nEdits for '//trim(topo_edits_file))
229 ncstatus = nf90_inq_varid(ncid,
'ni', id)
230 if (ncstatus /= nf90_noerr)
call mom_error(fatal,
'apply_topography_edits_from_file: '//&
231 'Failed to inq_varid ni for '//trim(topo_edits_file))
232 ncstatus = nf90_get_var(ncid, id, i)
233 if (ncstatus /= nf90_noerr)
call mom_error(fatal,
'apply_topography_edits_from_file: '//&
234 'Failed to get_var ni for '//trim(topo_edits_file))
235 if (i /= g%ieg)
call mom_error(fatal,
'apply_topography_edits_from_file: '//&
236 'Incompatible i-dimension of grid in '//trim(topo_edits_file))
239 ncstatus = nf90_inq_varid(ncid,
'nj', id)
240 if (ncstatus /= nf90_noerr)
call mom_error(fatal,
'apply_topography_edits_from_file: '//&
241 'Failed to inq_varid nj for '//trim(topo_edits_file))
242 ncstatus = nf90_get_var(ncid, id, j)
243 if (ncstatus /= nf90_noerr)
call mom_error(fatal,
'apply_topography_edits_from_file: '//&
244 'Failed to get_var nj for '//trim(topo_edits_file))
245 if (j /= g%jeg)
call mom_error(fatal,
'apply_topography_edits_from_file: '//&
246 'Incompatible j-dimension of grid in '//trim(topo_edits_file))
249 ncstatus = nf90_inq_varid(ncid,
'iEdit', id)
250 if (ncstatus /= nf90_noerr)
call mom_error(fatal,
'apply_topography_edits_from_file: '//&
251 'Failed to inq_varid iEdit for '//trim(topo_edits_file))
252 allocate(ig(n_edits))
253 ncstatus = nf90_get_var(ncid, id, ig)
254 if (ncstatus /= nf90_noerr)
call mom_error(fatal,
'apply_topography_edits_from_file: '//&
255 'Failed to get_var iEdit for '//trim(topo_edits_file))
258 ncstatus = nf90_inq_varid(ncid,
'jEdit', id)
259 if (ncstatus /= nf90_noerr)
call mom_error(fatal,
'apply_topography_edits_from_file: '//&
260 'Failed to inq_varid jEdit for '//trim(topo_edits_file))
261 allocate(jg(n_edits))
262 ncstatus = nf90_get_var(ncid, id, jg)
263 if (ncstatus /= nf90_noerr)
call mom_error(fatal,
'apply_topography_edits_from_file: '//&
264 'Failed to get_var jEdit for '//trim(topo_edits_file))
267 ncstatus = nf90_inq_varid(ncid,
'zEdit', id)
268 if (ncstatus /= nf90_noerr)
call mom_error(fatal,
'apply_topography_edits_from_file: '//&
269 'Failed to inq_varid zEdit for '//trim(topo_edits_file))
270 allocate(new_depth(n_edits))
271 ncstatus = nf90_get_var(ncid, id, new_depth)
272 if (ncstatus /= nf90_noerr)
call mom_error(fatal,
'apply_topography_edits_from_file: '//&
273 'Failed to get_var zEdit for '//trim(topo_edits_file))
276 ncstatus = nf90_close(ncid)
277 if (ncstatus /= nf90_noerr)
call mom_error(fatal,
'apply_topography_edits_from_file: '//&
278 'Failed to close '//trim(topo_edits_file))
281 i = ig(n) - g%isd_global + 2
282 j = jg(n) - g%jsd_global + 2
283 if (i>=g%isc .and. i<=g%iec .and. j>=g%jsc .and. j<=g%jec)
then
284 if (new_depth(n)/=0.)
then
285 write(*,
'(a,3i5,f8.2,a,f8.2,2i4)') &
286 'Ocean topography edit: ',n,ig(n),jg(n),d(i,j)/m_to_z,
'->',abs(new_depth(n)),i,j
287 d(i,j) = abs(m_to_z*new_depth(n))
289 call mom_error(fatal,
' apply_topography_edits_from_file: '//&
290 "A zero depth edit would change the land mask and is not allowed in"//trim(topo_edits_file))
295 deallocate( ig, jg, new_depth )
303 real,
dimension(G%isd:G%ied,G%jsd:G%jed), &
306 character(len=*),
intent(in) :: topog_config
308 real,
intent(in) :: max_depth
321 integer :: i, j, is, ie, js, je, isd, ied, jsd, jed
322 character(len=40) :: mdl =
"initialize_topography_named"
323 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
324 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
326 call calltree_enter(trim(mdl)//
"(), MOM_shared_initialization.F90")
327 call mom_mesg(
" MOM_shared_initialization.F90, initialize_topography_named: "//&
328 "TOPO_CONFIG = "//trim(topog_config), 5)
330 m_to_z = 1.0 ;
if (
present(us)) m_to_z = us%m_to_Z
332 call get_param(param_file, mdl,
"MINIMUM_DEPTH", min_depth, &
333 "The minimum depth of the ocean.", units=
"m", default=0.0, scale=m_to_z)
334 if (max_depth<=0.)
call mom_error(fatal,
"initialize_topography_named: "// &
335 "MAXIMUM_DEPTH has a non-sensical value! Was it set?")
337 if (trim(topog_config) /=
"flat")
then
338 call get_param(param_file, mdl,
"EDGE_DEPTH", dedge, &
339 "The depth at the edge of one of the named topographies.", &
340 units=
"m", default=100.0, scale=m_to_z)
355 call get_param(param_file, mdl,
"TOPOG_SLOPE_SCALE", expdecay, &
356 "The exponential decay scale used in defining some of "//&
357 "the named topographies.", units=
"m", default=400000.0)
363 if (trim(topog_config) ==
"flat")
then
364 do i=is,ie ;
do j=js,je ; d(i,j) = max_depth ;
enddo ;
enddo
365 elseif (trim(topog_config) ==
"spoon")
then
366 d0 = (max_depth - dedge) / &
367 ((1.0 - exp(-0.5*g%len_lat*g%Rad_earth*pi/(180.0 *expdecay))) * &
368 (1.0 - exp(-0.5*g%len_lat*g%Rad_earth*pi/(180.0 *expdecay))))
369 do i=is,ie ;
do j=js,je
372 d(i,j) = dedge + d0 * &
373 (sin(pi * (g%geoLonT(i,j) - (g%west_lon)) / g%len_lon) * &
374 (1.0 - exp((g%geoLatT(i,j) - (g%south_lat+g%len_lat))*g%Rad_earth*pi / &
377 elseif (trim(topog_config) ==
"bowl")
then
378 d0 = (max_depth - dedge) / &
379 ((1.0 - exp(-0.5*g%len_lat*g%Rad_earth*pi/(180.0 *expdecay))) * &
380 (1.0 - exp(-0.5*g%len_lat*g%Rad_earth*pi/(180.0 *expdecay))))
384 do i=is,ie ;
do j=js,je
385 d(i,j) = dedge + d0 * &
386 (sin(pi * (g%geoLonT(i,j) - g%west_lon) / g%len_lon) * &
387 ((1.0 - exp(-(g%geoLatT(i,j) - g%south_lat)*g%Rad_Earth*pi/ &
388 (180.0*expdecay))) * &
389 (1.0 - exp((g%geoLatT(i,j) - (g%south_lat+g%len_lat))* &
390 g%Rad_Earth*pi/(180.0*expdecay)))))
392 elseif (trim(topog_config) ==
"halfpipe")
then
393 d0 = max_depth - dedge
394 do i=is,ie ;
do j=js,je
395 d(i,j) = dedge + d0 * abs(sin(pi*(g%geoLatT(i,j) - g%south_lat)/g%len_lat))
398 call mom_error(fatal,
"initialize_topography_named: "// &
399 "Unrecognized topography name "//trim(topog_config))
403 do i=is,ie ;
do j=js,je
404 if (d(i,j) > max_depth) d(i,j) = max_depth
405 if (d(i,j) < min_depth) d(i,j) = 0.5*min_depth
416 real,
dimension(G%isd:G%ied,G%jsd:G%jed), &
419 real,
intent(in) :: max_depth
425 character(len=40) :: mdl =
"limit_topography"
426 real :: min_depth, mask_depth
428 call calltree_enter(trim(mdl)//
"(), MOM_shared_initialization.F90")
430 m_to_z = 1.0 ;
if (
present(us)) m_to_z = us%m_to_Z
432 call get_param(param_file, mdl,
"MINIMUM_DEPTH", min_depth, &
433 "If MASKING_DEPTH is unspecified, then anything shallower than "//&
434 "MINIMUM_DEPTH is assumed to be land and all fluxes are masked out. "//&
435 "If MASKING_DEPTH is specified, then all depths shallower than "//&
436 "MINIMUM_DEPTH but deeper than MASKING_DEPTH are rounded to MINIMUM_DEPTH.", &
437 units=
"m", default=0.0, scale=m_to_z)
438 call get_param(param_file, mdl,
"MASKING_DEPTH", mask_depth, &
439 "The depth below which to mask the ocean as land.", &
440 units=
"m", default=-9999.0, scale=m_to_z, do_not_log=.true.)
443 if (mask_depth < -9990.*m_to_z)
then
444 do j=g%jsd,g%jed ;
do i=g%isd,g%ied
445 d(i,j) = min( max( d(i,j), 0.5*min_depth ), max_depth )
448 do j=g%jsd,g%jed ;
do i=g%isd,g%ied
450 d(i,j) = min( max( d(i,j), min_depth ), max_depth )
465 real,
dimension(G%IsdB:G%IedB,G%JsdB:G%JedB), &
471 character(len=30) :: mdl =
"set_rotation_planetary"
477 call calltree_enter(trim(mdl)//
"(), MOM_shared_initialization.F90")
479 t_to_s = 1.0 ;
if (
present(us)) t_to_s = us%T_to_s
481 call get_param(param_file,
"set_rotation_planetary",
"OMEGA", omega, &
482 "The rotation rate of the earth.", units=
"s-1", &
483 default=7.2921e-5, scale=t_to_s)
486 do i=g%IsdB,g%IedB ;
do j=g%JsdB,g%JedB
487 f(i,j) = ( 2.0 * omega ) * sin( ( pi * g%geoLatBu(i,j) ) / 180.)
498 real,
dimension(G%IsdB:G%IedB,G%JsdB:G%JedB), &
507 real :: y_scl, rad_earth
510 character(len=40) :: mdl =
"set_rotation_beta_plane"
511 character(len=200) :: axis_units
513 call calltree_enter(trim(mdl)//
"(), MOM_shared_initialization.F90")
515 t_to_s = 1.0 ;
if (
present(us)) t_to_s = us%T_to_s
517 call get_param(param_file, mdl,
"F_0", f_0, &
518 "The reference value of the Coriolis parameter with the "//&
519 "betaplane option.", units=
"s-1", default=0.0, scale=t_to_s)
520 call get_param(param_file, mdl,
"BETA", beta, &
521 "The northward gradient of the Coriolis parameter with "//&
522 "the betaplane option.", units=
"m-1 s-1", default=0.0, scale=t_to_s)
523 call get_param(param_file, mdl,
"AXIS_UNITS", axis_units, default=
"degrees")
526 select case (axis_units(1:1))
528 call get_param(param_file, mdl,
"RAD_EARTH", rad_earth, &
529 "The radius of the Earth.", units=
"m", default=6.378e6)
531 case (
"k"); y_scl = 1.e3
532 case (
"m"); y_scl = 1.
533 case (
"c"); y_scl = 1.e-2
534 case default ;
call mom_error(fatal, &
535 " set_rotation_beta_plane: unknown AXIS_UNITS = "//trim(axis_units))
538 do i=g%IsdB,g%IedB ;
do j=g%JsdB,g%JedB
539 f(i,j) = f_0 + beta * ( g%geoLatBu(i,j) * y_scl )
552 real :: angle, lon_scale
556 character(len=40) :: mdl =
"initialize_grid_rotation_angle"
558 integer :: i, j, m, n
560 call get_param(pf, mdl,
"GRID_ROTATION_ANGLE_BUGS", use_bugs, &
561 "If true, use an older algorithm to calculate the sine and "//&
562 "cosines needed rotate between grid-oriented directions and "//&
563 "true north and east. Differences arise at the tripolar fold.", &
567 do j=g%jsc,g%jec ;
do i=g%isc,g%iec
568 lon_scale = cos((g%geoLatBu(i-1,j-1) + g%geoLatBu(i,j-1 ) + &
569 g%geoLatBu(i-1,j) + g%geoLatBu(i,j)) * atan(1.0)/180)
570 angle = atan2((g%geoLonBu(i-1,j) + g%geoLonBu(i,j) - &
571 g%geoLonBu(i-1,j-1) - g%geoLonBu(i,j-1))*lon_scale, &
572 g%geoLatBu(i-1,j) + g%geoLatBu(i,j) - &
573 g%geoLatBu(i-1,j-1) - g%geoLatBu(i,j-1) )
574 g%sin_rot(i,j) = sin(angle)
575 g%cos_rot(i,j) = cos(angle)
582 pi_720deg = atan(1.0) / 180.0
583 len_lon = 360.0 ;
if (g%len_lon > 0.0) len_lon = g%len_lon
584 do j=g%jsc,g%jec ;
do i=g%isc,g%iec
588 lon_scale = cos(pi_720deg*((g%geoLatBu(i-1,j-1) + g%geoLatBu(i,j)) + &
589 (g%geoLatBu(i,j-1) + g%geoLatBu(i-1,j)) ) )
590 angle = atan2(lon_scale*((lonb(1,2) - lonb(2,1)) + (lonb(2,2) - lonb(1,1))), &
591 (g%geoLatBu(i-1,j) - g%geoLatBu(i,j-1)) + &
592 (g%geoLatBu(i,j) - g%geoLatBu(i-1,j-1)) )
593 g%sin_rot(i,j) = sin(angle)
594 g%cos_rot(i,j) = cos(angle)
597 call pass_vector(g%cos_rot, g%sin_rot, g%Domain, stagger=agrid)
606 real,
intent(in) :: x
607 real,
intent(in) :: xc
608 real,
intent(in) :: lx
612 x_mod = modulo(x - (xc - 0.5*lx), lx) + (xc - 0.5*lx)
624 character(len=*),
intent(in) :: name
629 character(len=256) :: mesg
632 real :: dx_2 = -1.0, dy_2 = -1.0
634 integer :: option = -1
635 integer :: i, j, isd, ied, jsd, jed, isdb, iedb, jsdb, jedb
636 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
637 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
638 pi_180 = (4.0*atan(1.0))/180.0
640 select case ( trim(name) )
641 case (
"global_1deg") ; option = 1 ; dx_2 = 0.5*1.0
642 case default ;
call mom_error(fatal,
"reset_face_lengths_named: "//&
643 "Unrecognized channel configuration name "//trim(name))
646 m_to_l = 1.0 ;
if (
present(us)) m_to_l = us%m_to_L
647 l_to_m = 1.0 ;
if (
present(us)) l_to_m = us%L_to_m
650 do j=jsd,jed ;
do i=isdb,iedb
651 dy_2 = dx_2 * g%dyCu(i,j)*g%IdxCu(i,j) * cos(pi_180 * g%geoLatCu(i,j))
653 if ((abs(g%geoLatCu(i,j)-35.5) < dy_2) .and. (g%geoLonCu(i,j) < -4.5) .and. &
654 (g%geoLonCu(i,j) > -6.5)) &
655 g%dy_Cu(i,j) = g%mask2dCu(i,j)*12000.0*m_to_l
657 if ((abs(g%geoLatCu(i,j)-12.5) < dy_2) .and. (abs(g%geoLonCu(i,j)-43.0) < dx_2)) &
658 g%dy_Cu(i,j) = g%mask2dCu(i,j)*10000.0*m_to_l
660 if ((abs(g%geoLatCu(i,j)-40.5) < dy_2) .and. (abs(g%geoLonCu(i,j)-26.0) < dx_2)) &
661 g%dy_Cu(i,j) = g%mask2dCu(i,j)*5000.0*m_to_l
663 if ((abs(g%geoLatCu(i,j)-41.5) < dy_2) .and. (abs(g%geoLonCu(i,j)+220.0) < dx_2)) &
664 g%dy_Cu(i,j) = g%mask2dCu(i,j)*35000.0*m_to_l
666 if ((abs(g%geoLatCu(i,j)-45.5) < dy_2) .and. (abs(g%geoLonCu(i,j)+217.5) < 0.9)) &
667 g%dy_Cu(i,j) = g%mask2dCu(i,j)*15000.0*m_to_l
670 if ((abs(g%geoLatCu(i,j)-80.84) < 0.2) .and. (abs(g%geoLonCu(i,j)+64.9) < 0.8)) &
671 g%dy_Cu(i,j) = g%mask2dCu(i,j)*38000.0*m_to_l
675 do j=jsdb,jedb ;
do i=isd,ied
676 dy_2 = dx_2 * g%dyCv(i,j)*g%IdxCv(i,j) * cos(pi_180 * g%geoLatCv(i,j))
677 if ((abs(g%geoLatCv(i,j)-41.0) < dy_2) .and. (abs(g%geoLonCv(i,j)-28.5) < dx_2)) &
678 g%dx_Cv(i,j) = g%mask2dCv(i,j)*2500.0*m_to_l
680 if ((abs(g%geoLatCv(i,j)-13.0) < dy_2) .and. (abs(g%geoLonCv(i,j)-42.5) < dx_2)) &
681 g%dx_Cv(i,j) = g%mask2dCv(i,j)*10000.0*m_to_l
683 if ((abs(g%geoLatCv(i,j)+2.8) < 0.8) .and. (abs(g%geoLonCv(i,j)+241.5) < dx_2)) &
684 g%dx_Cv(i,j) = g%mask2dCv(i,j)*40000.0*m_to_l
686 if ((abs(g%geoLatCv(i,j)-0.56) < 0.5) .and. (abs(g%geoLonCv(i,j)+240.5) < dx_2)) &
687 g%dx_Cv(i,j) = g%mask2dCv(i,j)*80000.0*m_to_l
689 if ((abs(g%geoLatCv(i,j)-0.19) < 0.5) .and. (abs(g%geoLonCv(i,j)+230.5) < dx_2)) &
690 g%dx_Cv(i,j) = g%mask2dCv(i,j)*25000.0*m_to_l
692 if ((abs(g%geoLatCv(i,j)-0.19) < 0.5) .and. (abs(g%geoLonCv(i,j)+229.5) < dx_2)) &
693 g%dx_Cv(i,j) = g%mask2dCv(i,j)*25000.0*m_to_l
695 if ((abs(g%geoLatCv(i,j)-0.0) < 0.25) .and. (abs(g%geoLonCv(i,j)+228.5) < dx_2)) &
696 g%dx_Cv(i,j) = g%mask2dCv(i,j)*25000.0*m_to_l
698 if ((abs(g%geoLatCv(i,j)+8.5) < 0.5) .and. (abs(g%geoLonCv(i,j)+244.5) < dx_2)) &
699 g%dx_Cv(i,j) = g%mask2dCv(i,j)*20000.0*m_to_l
701 if ((abs(g%geoLatCv(i,j)+8.5) < 0.5) .and. (abs(g%geoLonCv(i,j)+235.5) < dx_2)) &
702 g%dx_Cv(i,j) = g%mask2dCv(i,j)*20000.0*m_to_l
704 if ((abs(g%geoLatCv(i,j)-52.5) < dy_2) .and. (abs(g%geoLonCv(i,j)+218.5) < dx_2)) &
705 g%dx_Cv(i,j) = g%mask2dCv(i,j)*2500.0*m_to_l
708 if ((abs(g%geoLatCv(i,j)-76.8) < 0.06) .and. (abs(g%geoLonCv(i,j)+88.7) < dx_2)) &
709 g%dx_Cv(i,j) = g%mask2dCv(i,j)*8400.0*m_to_l
716 do j=jsd,jed ;
do i=isdb,iedb
717 if (l_to_m*g%dy_Cu(i,j) > l_to_m*g%dyCu(i,j))
then
718 write(mesg,
'("dy_Cu of ",ES11.4," exceeds unrestricted width of ",ES11.4,&
719 &" by ",ES11.4," at lon/lat of ", ES11.4, ES11.4)') &
720 l_to_m*g%dy_Cu(i,j), l_to_m*g%dyCu(i,j), l_to_m*g%dy_Cu(i,j)-l_to_m*g%dyCu(i,j), &
721 g%geoLonCu(i,j), g%geoLatCu(i,j)
722 call mom_error(fatal,
"reset_face_lengths_named "//mesg)
724 g%areaCu(i,j) = g%dxCu(i,j) * g%dy_Cu(i,j)
726 if (g%areaCu(i,j) > 0.0) g%IareaCu(i,j) = g%mask2dCu(i,j) / (g%areaCu(i,j))
729 do j=jsdb,jedb ;
do i=isd,ied
730 if (l_to_m*g%dx_Cv(i,j) > l_to_m*g%dxCv(i,j))
then
731 write(mesg,
'("dx_Cv of ",ES11.4," exceeds unrestricted width of ",ES11.4,&
732 &" by ",ES11.4, " at lon/lat of ", ES11.4, ES11.4)') &
733 l_to_m*g%dx_Cv(i,j), l_to_m*g%dxCv(i,j), l_to_m*g%dx_Cv(i,j)-l_to_m*g%dxCv(i,j), &
734 g%geoLonCv(i,j), g%geoLatCv(i,j)
736 call mom_error(fatal,
"reset_face_lengths_named "//mesg)
738 g%areaCv(i,j) = g%dyCv(i,j) * g%dx_Cv(i,j)
740 if (g%areaCv(i,j) > 0.0) g%IareaCv(i,j) = g%mask2dCv(i,j) / (g%areaCv(i,j))
755 character(len=40) :: mdl =
"reset_face_lengths_file"
756 character(len=256) :: mesg
757 character(len=200) :: filename, chan_file, inputdir
760 integer :: i, j, isd, ied, jsd, jed, isdb, iedb, jsdb, jedb
761 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
762 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
765 call calltree_enter(trim(mdl)//
"(), MOM_shared_initialization.F90")
766 m_to_l = 1.0 ;
if (
present(us)) m_to_l = us%m_to_L
767 l_to_m = 1.0 ;
if (
present(us)) l_to_m = us%L_to_m
769 call get_param(param_file, mdl,
"CHANNEL_WIDTH_FILE", chan_file, &
770 "The file from which the list of narrowed channels is read.", &
771 default=
"ocean_geometry.nc")
772 call get_param(param_file, mdl,
"INPUTDIR", inputdir, default=
".")
773 inputdir = slasher(inputdir)
774 filename = trim(inputdir)//trim(chan_file)
775 call log_param(param_file, mdl,
"INPUTDIR/CHANNEL_WIDTH_FILE", filename)
777 if (is_root_pe())
then ;
if (.not.
file_exists(filename)) &
778 call mom_error(fatal,
" reset_face_lengths_file: Unable to open "//&
782 call mom_read_vector(filename,
"dyCuo",
"dxCvo", g%dy_Cu, g%dx_Cv, g%Domain, scale=m_to_l)
785 do j=jsd,jed ;
do i=isdb,iedb
786 if (l_to_m*g%dy_Cu(i,j) > l_to_m*g%dyCu(i,j))
then
787 write(mesg,
'("dy_Cu of ",ES11.4," exceeds unrestricted width of ",ES11.4,&
788 &" by ",ES11.4," at lon/lat of ", ES11.4, ES11.4)') &
789 l_to_m*g%dy_Cu(i,j), l_to_m*g%dyCu(i,j), l_to_m*g%dy_Cu(i,j)-l_to_m*g%dyCu(i,j), &
790 g%geoLonCu(i,j), g%geoLatCu(i,j)
791 call mom_error(fatal,
"reset_face_lengths_file "//mesg)
793 g%areaCu(i,j) = g%dxCu(i,j) * g%dy_Cu(i,j)
795 if (g%areaCu(i,j) > 0.0) g%IareaCu(i,j) = g%mask2dCu(i,j) / (g%areaCu(i,j))
798 do j=jsdb,jedb ;
do i=isd,ied
799 if (l_to_m*g%dx_Cv(i,j) > l_to_m*g%dxCv(i,j))
then
800 write(mesg,
'("dx_Cv of ",ES11.4," exceeds unrestricted width of ",ES11.4,&
801 &" by ",ES11.4, " at lon/lat of ", ES11.4, ES11.4)') &
802 l_to_m*g%dx_Cv(i,j), l_to_m*g%dxCv(i,j), l_to_m*g%dx_Cv(i,j)-l_to_m*g%dxCv(i,j), &
803 g%geoLonCv(i,j), g%geoLatCv(i,j)
805 call mom_error(fatal,
"reset_face_lengths_file "//mesg)
807 g%areaCv(i,j) = g%dyCv(i,j) * g%dx_Cv(i,j)
809 if (g%areaCv(i,j) > 0.0) g%IareaCv(i,j) = g%mask2dCv(i,j) / (g%areaCv(i,j))
825 character(len=120),
pointer,
dimension(:) :: lines => null()
826 character(len=120) :: line
827 character(len=200) :: filename, chan_file, inputdir, mesg
828 character(len=40) :: mdl =
"reset_face_lengths_list"
829 real,
pointer,
dimension(:,:) :: &
830 u_lat => null(), u_lon => null(), v_lat => null(), v_lon => null()
831 real,
pointer,
dimension(:) :: &
832 u_width => null(), v_width => null()
841 logical :: found_u, found_v
842 logical :: unit_in_use
843 integer :: ios, iounit, isu, isv
844 integer :: last, num_lines, nl_read, ln, npt, u_pt, v_pt
845 integer :: i, j, isd, ied, jsd, jed, isdb, iedb, jsdb, jedb
846 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
847 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
849 call calltree_enter(trim(mdl)//
"(), MOM_shared_initialization.F90")
850 m_to_l = 1.0 ;
if (
present(us)) m_to_l = us%m_to_L
851 l_to_m = 1.0 ;
if (
present(us)) l_to_m = us%L_to_m
853 call get_param(param_file, mdl,
"CHANNEL_LIST_FILE", chan_file, &
854 "The file from which the list of narrowed channels is read.", &
855 default=
"MOM_channel_list")
856 call get_param(param_file, mdl,
"INPUTDIR", inputdir, default=
".")
857 inputdir = slasher(inputdir)
858 filename = trim(inputdir)//trim(chan_file)
859 call log_param(param_file, mdl,
"INPUTDIR/CHANNEL_LIST_FILE", filename)
860 call get_param(param_file, mdl,
"CHANNEL_LIST_360_LON_CHECK", check_360, &
861 "If true, the channel configuration list works for any "//&
862 "longitudes in the range of -360 to 360.", default=.true.)
864 if (is_root_pe())
then
866 if (.not.
file_exists(filename))
call mom_error(fatal, &
867 " reset_face_lengths_list: Unable to open "//trim(filename))
871 INQUIRE(iounit,opened=unit_in_use) ;
if (.not.unit_in_use)
exit
873 if (iounit >= 512)
call mom_error(fatal, &
874 "reset_face_lengths_list: No unused file unit could be found.")
877 open(iounit, file=trim(filename), access=
'SEQUENTIAL', &
878 form=
'FORMATTED', action=
'READ', position=
'REWIND', iostat=ios)
879 if (ios /= 0)
call mom_error(fatal, &
880 "reset_face_lengths_list: Error opening "//trim(filename))
886 len_lon = 360.0 ;
if (g%len_lon > 0.0) len_lon = g%len_lon
887 len_lat = 180.0 ;
if (g%len_lat > 0.0) len_lat = g%len_lat
889 call broadcast(num_lines, root_pe())
891 if (num_lines > 0)
then
892 allocate (lines(num_lines))
893 if (num_lines > 0)
then
894 allocate(u_lat(2,num_lines)) ; u_lat(:,:) = -1e34
895 allocate(u_lon(2,num_lines)) ; u_lon(:,:) = -1e34
896 allocate(u_width(num_lines)) ; u_width(:) = -1e34
898 allocate(v_lat(2,num_lines)) ; v_lat(:,:) = -1e34
899 allocate(v_lon(2,num_lines)) ; v_lon(:,:) = -1e34
900 allocate(v_width(num_lines)) ; v_width(:) = -1e34
904 if (is_root_pe())
then
906 if (nl_read /= num_lines) &
907 call mom_error(fatal,
'reset_face_lengths_list : Found different '// &
908 'number of valid lines on second reading of '//trim(filename))
909 close(iounit) ; iounit = -1
913 call broadcast(lines, 120, root_pe())
919 found_u = .false.; found_v = .false.
920 isu = index(
uppercase(line),
"U_WIDTH" );
if (isu > 0) found_u = .true.
921 isv = index(
uppercase(line),
"V_WIDTH" );
if (isv > 0) found_v = .true.
926 read(line(isu+8:),*) u_lon(1:2,u_pt), u_lat(1:2,u_pt), u_width(u_pt)
927 if (is_root_pe())
then
929 if ((abs(u_lon(1,u_pt)) > len_lon) .or. (abs(u_lon(2,u_pt)) > len_lon)) &
930 call mom_error(warning,
"reset_face_lengths_list : Out-of-bounds "//&
931 "u-longitude found when reading line "//trim(line)//
" from file "//&
933 if ((abs(u_lat(1,u_pt)) > len_lat) .or. (abs(u_lat(2,u_pt)) > len_lat)) &
934 call mom_error(warning,
"reset_face_lengths_list : Out-of-bounds "//&
935 "u-latitude found when reading line "//trim(line)//
" from file "//&
938 if (u_lat(1,u_pt) > u_lat(2,u_pt)) &
939 call mom_error(warning,
"reset_face_lengths_list : Out-of-order "//&
940 "u-face latitudes found when reading line "//trim(line)//
" from file "//&
942 if (u_lon(1,u_pt) > u_lon(2,u_pt)) &
943 call mom_error(warning,
"reset_face_lengths_list : Out-of-order "//&
944 "u-face longitudes found when reading line "//trim(line)//
" from file "//&
946 if (u_width(u_pt) < 0.0) &
947 call mom_error(warning,
"reset_face_lengths_list : Negative "//&
948 "u-width found when reading line "//trim(line)//
" from file "//&
951 elseif (found_v)
then
953 read(line(isv+8:),*) v_lon(1:2,v_pt), v_lat(1:2,v_pt), v_width(v_pt)
954 if (is_root_pe())
then
956 if ((abs(v_lon(1,v_pt)) > len_lon) .or. (abs(v_lon(2,v_pt)) > len_lon)) &
957 call mom_error(warning,
"reset_face_lengths_list : Out-of-bounds "//&
958 "v-longitude found when reading line "//trim(line)//
" from file "//&
960 if ((abs(v_lat(1,v_pt)) > len_lat) .or. (abs(v_lat(2,v_pt)) > len_lat)) &
961 call mom_error(warning,
"reset_face_lengths_list : Out-of-bounds "//&
962 "v-latitude found when reading line "//trim(line)//
" from file "//&
965 if (v_lat(1,v_pt) > v_lat(2,v_pt)) &
966 call mom_error(warning,
"reset_face_lengths_list : Out-of-order "//&
967 "v-face latitudes found when reading line "//trim(line)//
" from file "//&
969 if (v_lon(1,v_pt) > v_lon(2,v_pt)) &
970 call mom_error(warning,
"reset_face_lengths_list : Out-of-order "//&
971 "v-face longitudes found when reading line "//trim(line)//
" from file "//&
973 if (v_width(v_pt) < 0.0) &
974 call mom_error(warning,
"reset_face_lengths_list : Negative "//&
975 "v-width found when reading line "//trim(line)//
" from file "//&
984 do j=jsd,jed ;
do i=isdb,iedb
985 lat = g%geoLatCu(i,j) ; lon = g%geoLonCu(i,j)
986 if (check_360)
then ; lon_p = lon+len_lon ; lon_m = lon-len_lon
987 else ; lon_p = lon ; lon_m = lon ;
endif
990 if (((lat >= u_lat(1,npt)) .and. (lat <= u_lat(2,npt))) .and. &
991 (((lon >= u_lon(1,npt)) .and. (lon <= u_lon(2,npt))) .or. &
992 ((lon_p >= u_lon(1,npt)) .and. (lon_p <= u_lon(2,npt))) .or. &
993 ((lon_m >= u_lon(1,npt)) .and. (lon_m <= u_lon(2,npt)))) )
then
995 g%dy_Cu(i,j) = g%mask2dCu(i,j) * m_to_l*min(l_to_m*g%dyCu(i,j), max(u_width(npt), 0.0))
996 if (j>=g%jsc .and. j<=g%jec .and. i>=g%isc .and. i<=g%iec)
then
997 if ( g%mask2dCu(i,j) == 0.0 )
then
998 write(*,
'(A,2F8.2,A,4F8.2,A)')
"read_face_lengths_list : G%mask2dCu=0 at ",lat,lon,
" (",&
999 u_lat(1,npt), u_lat(2,npt), u_lon(1,npt), u_lon(2,npt),
") so grid metric is unmodified."
1001 write(*,
'(A,2F8.2,A,4F8.2,A5,F9.2,A1)') &
1002 "read_face_lengths_list : Modifying dy_Cu gridpoint at ",lat,lon,
" (",&
1003 u_lat(1,npt), u_lat(2,npt), u_lon(1,npt), u_lon(2,npt),
") to ",l_to_m*g%dy_Cu(i,j),
"m"
1009 g%areaCu(i,j) = g%dxCu(i,j) * g%dy_Cu(i,j)
1010 g%IareaCu(i,j) = 0.0
1011 if (g%areaCu(i,j) > 0.0) g%IareaCu(i,j) = g%mask2dCu(i,j) / (g%areaCu(i,j))
1014 do j=jsdb,jedb ;
do i=isd,ied
1015 lat = g%geoLatCv(i,j) ; lon = g%geoLonCv(i,j)
1016 if (check_360)
then ; lon_p = lon+len_lon ; lon_m = lon-len_lon
1017 else ; lon_p = lon ; lon_m = lon ;
endif
1020 if (((lat >= v_lat(1,npt)) .and. (lat <= v_lat(2,npt))) .and. &
1021 (((lon >= v_lon(1,npt)) .and. (lon <= v_lon(2,npt))) .or. &
1022 ((lon_p >= v_lon(1,npt)) .and. (lon_p <= v_lon(2,npt))) .or. &
1023 ((lon_m >= v_lon(1,npt)) .and. (lon_m <= v_lon(2,npt)))) )
then
1024 g%dx_Cv(i,j) = g%mask2dCv(i,j) * m_to_l*min(l_to_m*g%dxCv(i,j), max(v_width(npt), 0.0))
1025 if (i>=g%isc .and. i<=g%iec .and. j>=g%jsc .and. j<=g%jec)
then
1026 if ( g%mask2dCv(i,j) == 0.0 )
then
1027 write(*,
'(A,2F8.2,A,4F8.2,A)')
"read_face_lengths_list : G%mask2dCv=0 at ",lat,lon,
" (",&
1028 v_lat(1,npt), v_lat(2,npt), v_lon(1,npt), v_lon(2,npt),
") so grid metric is unmodified."
1030 write(*,
'(A,2F8.2,A,4F8.2,A5,F9.2,A1)') &
1031 "read_face_lengths_list : Modifying dx_Cv gridpoint at ",lat,lon,
" (",&
1032 v_lat(1,npt), v_lat(2,npt), v_lon(1,npt), v_lon(2,npt),
") to ",l_to_m*g%dx_Cv(i,j),
"m"
1038 g%areaCv(i,j) = g%dyCv(i,j) * g%dx_Cv(i,j)
1039 g%IareaCv(i,j) = 0.0
1040 if (g%areaCv(i,j) > 0.0) g%IareaCv(i,j) = g%mask2dCv(i,j) / (g%areaCv(i,j))
1043 if (num_lines > 0)
then
1044 deallocate(u_lat) ;
deallocate(u_lon) ;
deallocate(u_width)
1045 deallocate(v_lat) ;
deallocate(v_lon) ;
deallocate(v_width)
1055 integer,
intent(in) :: iounit
1056 character(len=*),
intent(in) :: filename
1057 integer,
intent(out) :: num_lines
1058 character(len=120),
dimension(:),
pointer :: lines
1062 character(len=120) :: line, line_up
1063 logical :: found_u, found_v
1064 integer :: isu, isv, icom, verbose
1069 if (iounit <= 0)
return
1072 read(iounit,
'(a)',
end=8, err=9) line
1073 last = len_trim(line)
1075 icom = index(line(:last),
"!") ;
if (icom > 0) last = icom-1
1076 icom = index(line(:last),
"/*") ;
if (icom > 0) last = icom-1
1081 found_u = .false.; found_v = .false.
1082 isu = index(line_up(:last),
"U_WIDTH" );
if (isu > 0) found_u = .true.
1083 isv = index(line_up(:last),
"V_WIDTH" );
if (isv > 0) found_v = .true.
1085 if (found_u .and. found_v)
call mom_error(fatal, &
1086 "read_face_length_list : both U_WIDTH and V_WIDTH found when "//&
1087 "reading the line "//trim(line(:last))//
" in file "//trim(filename))
1088 if (found_u .or. found_v)
then
1089 num_lines = num_lines + 1
1090 if (
associated(lines))
then
1091 lines(num_lines) = line(1:last)
1099 9
call mom_error(fatal,
"read_face_length_list : "//&
1100 "Error while reading file "//trim(filename))
1109 type(dyn_horgrid_type),
intent(inout) :: g
1114 do i=g%isd,g%ied-1 ;
do j=g%jsd,g%jed
1115 g%Dblock_u(i,j) = g%mask2dCu(i,j) * max(g%bathyT(i,j), g%bathyT(i+1,j))
1116 g%Dopen_u(i,j) = g%Dblock_u(i,j)
1118 do i=g%isd,g%ied ;
do j=g%jsd,g%jed-1
1119 g%Dblock_v(i,j) = g%mask2dCv(i,j) * max(g%bathyT(i,j), g%bathyT(i,j+1))
1120 g%Dopen_v(i,j) = g%Dblock_v(i,j)
1129 type(dyn_horgrid_type),
intent(inout) :: g
1134 do i=g%isd,g%ied-1 ;
do j=g%jsd,g%jed
1135 g%Dblock_u(i,j) = g%mask2dCu(i,j) * min(g%bathyT(i,j), g%bathyT(i+1,j))
1136 g%Dopen_u(i,j) = g%Dblock_u(i,j)
1138 do i=g%isd,g%ied ;
do j=g%jsd,g%jed-1
1139 g%Dblock_v(i,j) = g%mask2dCv(i,j) * min(g%bathyT(i,j), g%bathyT(i,j+1))
1140 g%Dopen_v(i,j) = g%Dblock_v(i,j)
1149 type(dyn_horgrid_type),
intent(inout) :: g
1150 type(unit_scale_type),
optional,
intent(in) :: us
1153 real,
dimension(G%isc:G%iec, G%jsc:G%jec) :: tmpforsumming
1157 area_scale = 1.0 ;
if (
present(us)) area_scale = us%L_to_m**2
1159 tmpforsumming(:,:) = 0.
1160 g%areaT_global = 0.0 ; g%IareaT_global = 0.0
1161 do j=g%jsc,g%jec ;
do i=g%isc,g%iec
1162 tmpforsumming(i,j) = area_scale*g%areaT(i,j) * g%mask2dT(i,j)
1164 g%areaT_global = reproducing_sum(tmpforsumming)
1166 if (g%areaT_global == 0.0) &
1167 call mom_error(fatal,
"compute_global_grid_integrals: "//&
1168 "zero ocean area (check topography?)")
1170 g%IareaT_global = 1.0 / (g%areaT_global)
1178 type(dyn_horgrid_type),
intent(inout) :: g
1180 character(len=*),
intent(in) :: directory
1181 character(len=*),
optional,
intent(in) :: geom_file
1183 type(unit_scale_type),
optional,
intent(in) :: us
1186 character(len=240) :: filepath
1187 character(len=40) :: mdl =
"write_ocean_geometry_file"
1188 integer,
parameter :: nflds=23
1189 type(vardesc) :: vars(nflds)
1190 type(fieldtype) :: fields(nflds)
1191 real :: z_to_m_scale
1192 real :: s_to_t_scale
1193 real :: l_to_m_scale
1195 integer :: file_threading
1196 integer :: nflds_used
1197 integer :: i, j, is, ie, js, je, isq, ieq, jsq, jeq
1198 integer :: isd, ied, jsd, jed, isdb, iedb, jsdb, jedb
1199 logical :: multiple_files
1200 real,
dimension(G%isd :G%ied ,G%jsd :G%jed ) :: out_h
1201 real,
dimension(G%IsdB:G%IedB,G%JsdB:G%JedB) :: out_q
1202 real,
dimension(G%IsdB:G%IedB,G%jsd :G%jed ) :: out_u
1203 real,
dimension(G%isd :G%ied ,G%JsdB:G%JedB) :: out_v
1205 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
1206 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
1207 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
1208 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
1210 z_to_m_scale = 1.0 ;
if (
present(us)) z_to_m_scale = us%Z_to_m
1211 s_to_t_scale = 1.0 ;
if (
present(us)) s_to_t_scale = us%s_to_T
1212 l_to_m_scale = 1.0 ;
if (
present(us)) l_to_m_scale = us%L_to_m
1225 vars(1) = var_desc(
"geolatb",
"degree",
"latitude at corner (Bu) points",
'q',
'1',
'1')
1226 vars(2) = var_desc(
"geolonb",
"degree",
"longitude at corner (Bu) points",
'q',
'1',
'1')
1227 vars(3) = var_desc(
"geolat",
"degree",
"latitude at tracer (T) points",
'h',
'1',
'1')
1228 vars(4) = var_desc(
"geolon",
"degree",
"longitude at tracer (T) points",
'h',
'1',
'1')
1229 vars(5) = var_desc(
"D",
"meter",
"Basin Depth",
'h',
'1',
'1')
1230 vars(6) = var_desc(
"f",
"s-1",
"Coriolis Parameter",
'q',
'1',
'1')
1231 vars(7) = var_desc(
"dxCv",
"m",
"Zonal grid spacing at v points",
'v',
'1',
'1')
1232 vars(8) = var_desc(
"dyCu",
"m",
"Meridional grid spacing at u points",
'u',
'1',
'1')
1233 vars(9) = var_desc(
"dxCu",
"m",
"Zonal grid spacing at u points",
'u',
'1',
'1')
1234 vars(10)= var_desc(
"dyCv",
"m",
"Meridional grid spacing at v points",
'v',
'1',
'1')
1235 vars(11)= var_desc(
"dxT",
"m",
"Zonal grid spacing at h points",
'h',
'1',
'1')
1236 vars(12)= var_desc(
"dyT",
"m",
"Meridional grid spacing at h points",
'h',
'1',
'1')
1237 vars(13)= var_desc(
"dxBu",
"m",
"Zonal grid spacing at q points",
'q',
'1',
'1')
1238 vars(14)= var_desc(
"dyBu",
"m",
"Meridional grid spacing at q points",
'q',
'1',
'1')
1239 vars(15)= var_desc(
"Ah",
"m2",
"Area of h cells",
'h',
'1',
'1')
1240 vars(16)= var_desc(
"Aq",
"m2",
"Area of q cells",
'q',
'1',
'1')
1242 vars(17)= var_desc(
"dxCvo",
"m",
"Open zonal grid spacing at v points",
'v',
'1',
'1')
1243 vars(18)= var_desc(
"dyCuo",
"m",
"Open meridional grid spacing at u points",
'u',
'1',
'1')
1244 vars(19)= var_desc(
"wet",
"nondim",
"land or ocean?",
'h',
'1',
'1')
1246 vars(20) = var_desc(
"Dblock_u",
"m",
"Blocked depth at u points",
'u',
'1',
'1')
1247 vars(21) = var_desc(
"Dopen_u",
"m",
"Open depth at u points",
'u',
'1',
'1')
1248 vars(22) = var_desc(
"Dblock_v",
"m",
"Blocked depth at v points",
'v',
'1',
'1')
1249 vars(23) = var_desc(
"Dopen_v",
"m",
"Open depth at v points",
'v',
'1',
'1')
1252 nflds_used = 19 ;
if (g%bathymetry_at_vel) nflds_used = 23
1254 if (
present(geom_file))
then
1255 filepath = trim(directory) // trim(geom_file)
1257 filepath = trim(directory) //
"ocean_geometry"
1265 call get_param(param_file, mdl,
"PARALLEL_RESTARTFILES", multiple_files, &
1266 "If true, each processor writes its own restart file, "//&
1267 "otherwise a single restart file is generated", &
1269 file_threading = single_file
1270 if (multiple_files) file_threading = multiple
1272 call create_file(unit, trim(filepath), vars, nflds_used, fields, &
1273 file_threading, dg=g)
1275 do j=jsq,jeq;
do i=isq,ieq; out_q(i,j) = g%geoLatBu(i,j);
enddo ;
enddo
1276 call write_field(unit, fields(1), g%Domain%mpp_domain, out_q)
1277 do j=jsq,jeq;
do i=isq,ieq; out_q(i,j) = g%geoLonBu(i,j);
enddo ;
enddo
1278 call write_field(unit, fields(2), g%Domain%mpp_domain, out_q)
1279 call write_field(unit, fields(3), g%Domain%mpp_domain, g%geoLatT)
1280 call write_field(unit, fields(4), g%Domain%mpp_domain, g%geoLonT)
1282 do j=js,je ;
do i=is,ie ; out_h(i,j) = z_to_m_scale*g%bathyT(i,j) ;
enddo ;
enddo
1283 call write_field(unit, fields(5), g%Domain%mpp_domain, out_h)
1284 do j=jsq,jeq ;
do i=isq,ieq ; out_q(i,j) = s_to_t_scale*g%CoriolisBu(i,j) ;
enddo ;
enddo
1285 call write_field(unit, fields(6), g%Domain%mpp_domain, out_q)
1290 do j=jsq,jeq ;
do i=is,ie ; out_v(i,j) = l_to_m_scale*g%dxCv(i,j) ;
enddo ;
enddo
1291 call write_field(unit, fields(7), g%Domain%mpp_domain, out_v)
1292 do j=js,je ;
do i=isq,ieq ; out_u(i,j) = l_to_m_scale*g%dyCu(i,j) ;
enddo ;
enddo
1293 call write_field(unit, fields(8), g%Domain%mpp_domain, out_u)
1295 do j=js,je ;
do i=isq,ieq ; out_u(i,j) = l_to_m_scale*g%dxCu(i,j) ;
enddo ;
enddo
1296 call write_field(unit, fields(9), g%Domain%mpp_domain, out_u)
1297 do j=jsq,jeq ;
do i=is,ie ; out_v(i,j) = l_to_m_scale*g%dyCv(i,j) ;
enddo ;
enddo
1298 call write_field(unit, fields(10), g%Domain%mpp_domain, out_v)
1300 do j=js,je ;
do i=is,ie ; out_h(i,j) = l_to_m_scale*g%dxT(i,j);
enddo ;
enddo
1301 call write_field(unit, fields(11), g%Domain%mpp_domain, out_h)
1302 do j=js,je ;
do i=is,ie ; out_h(i,j) = l_to_m_scale*g%dyT(i,j) ;
enddo ;
enddo
1303 call write_field(unit, fields(12), g%Domain%mpp_domain, out_h)
1305 do j=jsq,jeq ;
do i=isq,ieq ; out_q(i,j) = l_to_m_scale*g%dxBu(i,j) ;
enddo ;
enddo
1306 call write_field(unit, fields(13), g%Domain%mpp_domain, out_q)
1307 do j=jsq,jeq ;
do i=isq,ieq ; out_q(i,j) = l_to_m_scale*g%dyBu(i,j) ;
enddo ;
enddo
1308 call write_field(unit, fields(14), g%Domain%mpp_domain, out_q)
1310 do j=js,je ;
do i=is,ie ; out_h(i,j) = g%areaT(i,j) ;
enddo ;
enddo
1311 call write_field(unit, fields(15), g%Domain%mpp_domain, out_h)
1312 do j=jsq,jeq ;
do i=isq,ieq ; out_q(i,j) = g%areaBu(i,j) ;
enddo ;
enddo
1313 call write_field(unit, fields(16), g%Domain%mpp_domain, out_q)
1315 do j=jsq,jeq ;
do i=is,ie ; out_v(i,j) = l_to_m_scale*g%dx_Cv(i,j) ;
enddo ;
enddo
1316 call write_field(unit, fields(17), g%Domain%mpp_domain, out_v)
1317 do j=js,je ;
do i=isq,ieq ; out_u(i,j) = l_to_m_scale*g%dy_Cu(i,j) ;
enddo ;
enddo
1318 call write_field(unit, fields(18), g%Domain%mpp_domain, out_u)
1319 call write_field(unit, fields(19), g%Domain%mpp_domain, g%mask2dT)
1321 if (g%bathymetry_at_vel)
then
1322 do j=js,je ;
do i=isq,ieq ; out_u(i,j) = z_to_m_scale*g%Dblock_u(i,j) ;
enddo ;
enddo
1323 call write_field(unit, fields(20), g%Domain%mpp_domain, out_u)
1324 do j=js,je ;
do i=isq,ieq ; out_u(i,j) = z_to_m_scale*g%Dopen_u(i,j) ;
enddo ;
enddo
1325 call write_field(unit, fields(21), g%Domain%mpp_domain, out_u)
1326 do j=jsq,jeq ;
do i=is,ie ; out_v(i,j) = z_to_m_scale*g%Dblock_v(i,j) ;
enddo ;
enddo
1327 call write_field(unit, fields(22), g%Domain%mpp_domain, out_v)
1328 do j=jsq,jeq ;
do i=is,ie ; out_v(i,j) = z_to_m_scale*g%Dopen_v(i,j) ;
enddo ;
enddo
1329 call write_field(unit, fields(23), g%Domain%mpp_domain, out_v)
1332 call close_file(unit)