55 use mom_cpu_clock,
only : cpu_clock_id, cpu_clock_begin, cpu_clock_end
57 use mom_cpu_clock,
only : clock_module_driver, clock_module, clock_routine
65 use mom_domains,
only : to_south, to_west, to_all, cgrid_ne, scalar_pair
70 use mom_io,
only : mom_io_init
74 use mom_time_manager,
only :
operator(-),
operator(>),
operator(*),
operator(/)
100 implicit none ;
private
102 #include <MOM_memory.h>
106 real allocable_,
dimension(NIMEMB_PTR_,NJMEM_,NKMEM_) :: &
107 cau, & !< CAu = f*v - u.grad(u) [m s-2].
108 pfu, & !< PFu = -dM/dx [m s-2].
111 real allocable_,
dimension(NIMEM_,NJMEMB_PTR_,NKMEM_) :: &
112 cav, & !< CAv = -f*u - u.grad(v) [m s-2].
113 pfv, & !< PFv = -dM/dy [m s-2].
116 real,
pointer,
dimension(:,:) :: taux_bot => null()
117 real,
pointer,
dimension(:,:) :: tauy_bot => null()
127 logical :: module_is_initialized = .false.
130 integer :: id_uh = -1, id_vh = -1
131 integer :: id_pfu = -1, id_pfv = -1, id_cau = -1, id_cav = -1
161 type(
ale_cs),
pointer :: ale_csp => null()
173 public step_mom_dyn_unsplit_rk2, register_restarts_dyn_unsplit_rk2
174 public initialize_dyn_unsplit_rk2, end_dyn_unsplit_rk2
177 integer :: id_clock_cor, id_clock_pres, id_clock_vertvisc
178 integer :: id_clock_horvisc, id_clock_continuity, id_clock_mom_update
179 integer :: id_clock_pass, id_clock_pass_init
187 subroutine step_mom_dyn_unsplit_rk2(u_in, v_in, h_in, tv, visc, Time_local, dt, forces, &
188 p_surf_begin, p_surf_end, uh, vh, uhtr, vhtr, eta_av, G, GV, US, CS, &
193 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)),
intent(inout) :: u_in
195 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)),
intent(inout) :: v_in
197 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(inout) :: h_in
205 type(time_type),
intent(in) :: time_local
207 real,
intent(in) :: dt
209 real,
dimension(:,:),
pointer :: p_surf_begin
212 real,
dimension(:,:),
pointer :: p_surf_end
215 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)),
intent(inout) :: uh
217 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)),
intent(inout) :: vh
219 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)),
intent(inout) :: uhtr
222 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)),
intent(inout) :: vhtr
225 real,
dimension(SZI_(G),SZJ_(G)),
intent(out) :: eta_av
236 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)) :: h_av, hp
237 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)) :: up
238 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)) :: vp
239 real,
dimension(:,:),
pointer :: p_surf => null()
242 logical :: dyn_p_surf
243 integer :: i, j, k, is, ie, js, je, isq, ieq, jsq, jeq, nz
244 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
245 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
248 h_av(:,:,:) = 0; hp(:,:,:) = 0
252 dyn_p_surf =
associated(p_surf_begin) .and.
associated(p_surf_end)
254 call safe_alloc_ptr(p_surf,g%isd,g%ied,g%jsd,g%jed) ; p_surf(:,:) = 0.0
256 p_surf => forces%p_surf
267 call enable_averaging(dt,time_local, cs%diag)
268 call cpu_clock_begin(id_clock_horvisc)
269 call horizontal_viscosity(u_in, v_in, h_in, cs%diffu, cs%diffv, meke, varmix, &
270 g, gv, us, cs%hor_visc_CSp)
271 call cpu_clock_end(id_clock_horvisc)
272 call disable_averaging(cs%diag)
273 call pass_vector(cs%diffu, cs%diffv, g%Domain, clock=id_clock_pass)
279 call cpu_clock_begin(id_clock_continuity)
282 call continuity(u_in, v_in, h_in, hp, uh, vh, dt_pred, g, gv, us, cs%continuity_CSp, &
284 call cpu_clock_end(id_clock_continuity)
285 call pass_var(hp, g%Domain, clock=id_clock_pass)
286 call pass_vector(uh, vh, g%Domain, clock=id_clock_pass)
289 call cpu_clock_begin(id_clock_mom_update)
291 do j=js-2,je+2 ;
do i=is-2,ie+2
292 h_av(i,j,k) = (h_in(i,j,k) + hp(i,j,k)) * 0.5
293 enddo ;
enddo ;
enddo
294 call cpu_clock_end(id_clock_mom_update)
297 call cpu_clock_begin(id_clock_cor)
298 call coradcalc(u_in, v_in, h_av, uh, vh, cs%CAu, cs%CAv, cs%OBC, cs%ADp, &
299 g, gv, us, cs%CoriolisAdv_CSp)
300 call cpu_clock_end(id_clock_cor)
303 call cpu_clock_begin(id_clock_pres)
304 if (dyn_p_surf)
then ;
do j=js-2,je+2 ;
do i=is-2,ie+2
305 p_surf(i,j) = 0.5*p_surf_begin(i,j) + 0.5*p_surf_end(i,j)
306 enddo ;
enddo ;
endif
307 call pressureforce(h_in, tv, cs%PFu, cs%PFv, g, gv, us, &
308 cs%PressureForce_CSp, cs%ALE_CSp, p_surf)
309 call cpu_clock_end(id_clock_pres)
310 call pass_vector(cs%PFu, cs%PFv, g%Domain, clock=id_clock_pass)
311 call pass_vector(cs%CAu, cs%CAv, g%Domain, clock=id_clock_pass)
313 if (
associated(cs%OBC)) then;
if (cs%OBC%update_OBC)
then
314 call update_obc_data(cs%OBC, g, gv, us, tv, h_in, cs%update_OBC_CSp, time_local)
316 if (
associated(cs%OBC))
then
317 call open_boundary_zero_normal_flow(cs%OBC, g, cs%PFu, cs%PFv)
318 call open_boundary_zero_normal_flow(cs%OBC, g, cs%CAu, cs%CAv)
319 call open_boundary_zero_normal_flow(cs%OBC, g, cs%diffu, cs%diffv)
323 call cpu_clock_begin(id_clock_mom_update)
324 do k=1,nz ;
do j=js,je ;
do i=isq,ieq
325 up(i,j,k) = g%mask2dCu(i,j) * (u_in(i,j,k) + dt_pred * &
326 ((cs%PFu(i,j,k) + cs%CAu(i,j,k)) + us%s_to_T*cs%diffu(i,j,k)))
327 enddo ;
enddo ;
enddo
328 do k=1,nz ;
do j=jsq,jeq ;
do i=is,ie
329 vp(i,j,k) = g%mask2dCv(i,j) * (v_in(i,j,k) + dt_pred * &
330 ((cs%PFv(i,j,k) + cs%CAv(i,j,k)) + us%s_to_T*cs%diffv(i,j,k)))
331 enddo ;
enddo ;
enddo
332 call cpu_clock_end(id_clock_mom_update)
335 call mom_accel_chksum(
"Predictor 1 accel", cs%CAu, cs%CAv, cs%PFu, cs%PFv,&
336 cs%diffu, cs%diffv, g, gv, us)
339 call cpu_clock_begin(id_clock_vertvisc)
340 call enable_averaging(dt, time_local, cs%diag)
341 call set_viscous_ml(up, vp, h_av, tv, forces, visc, dt_pred, g, gv, us, &
343 call disable_averaging(cs%diag)
344 call vertvisc_coef(up, vp, h_av, forces, visc, dt_pred, g, gv, us, &
345 cs%vertvisc_CSp, cs%OBC)
346 call vertvisc(up, vp, h_av, forces, visc, dt_pred, cs%OBC, cs%ADp, cs%CDp, &
347 g, gv, us, cs%vertvisc_CSp)
348 call cpu_clock_end(id_clock_vertvisc)
349 call pass_vector(up, vp, g%Domain, clock=id_clock_pass)
353 call cpu_clock_begin(id_clock_continuity)
354 call continuity(up, vp, h_in, hp, uh, vh, &
355 dt, g, gv, us, cs%continuity_CSp, obc=cs%OBC)
356 call cpu_clock_end(id_clock_continuity)
357 call pass_var(hp, g%Domain, clock=id_clock_pass)
358 call pass_vector(uh, vh, g%Domain, clock=id_clock_pass)
361 do k=1,nz ;
do j=js-2,je+2 ;
do i=is-2,ie+2
362 h_av(i,j,k) = (h_in(i,j,k) + hp(i,j,k)) * 0.5
363 enddo ;
enddo ;
enddo
369 call cpu_clock_begin(id_clock_cor)
370 call coradcalc(up, vp, h_av, uh, vh, cs%CAu, cs%CAv, cs%OBC, cs%ADp, &
371 g, gv, us, cs%CoriolisAdv_CSp)
372 call cpu_clock_end(id_clock_cor)
373 if (
associated(cs%OBC))
then
374 call open_boundary_zero_normal_flow(cs%OBC, g, cs%CAu, cs%CAv)
381 do k=1,nz ;
do j=js,je ;
do i=isq,ieq
382 up(i,j,k) = g%mask2dCu(i,j) * (u_in(i,j,k) + dt * (1.+cs%begw) * &
383 ((cs%PFu(i,j,k) + cs%CAu(i,j,k)) + us%s_to_T*cs%diffu(i,j,k)))
384 u_in(i,j,k) = g%mask2dCu(i,j) * (u_in(i,j,k) + dt * &
385 ((cs%PFu(i,j,k) + cs%CAu(i,j,k)) + us%s_to_T*cs%diffu(i,j,k)))
386 enddo ;
enddo ;
enddo
387 do k=1,nz ;
do j=jsq,jeq ;
do i=is,ie
388 vp(i,j,k) = g%mask2dCv(i,j) * (v_in(i,j,k) + dt * (1.+cs%begw) * &
389 ((cs%PFv(i,j,k) + cs%CAv(i,j,k)) + us%s_to_T*cs%diffv(i,j,k)))
390 v_in(i,j,k) = g%mask2dCv(i,j) * (v_in(i,j,k) + dt * &
391 ((cs%PFv(i,j,k) + cs%CAv(i,j,k)) + us%s_to_T*cs%diffv(i,j,k)))
392 enddo ;
enddo ;
enddo
396 call cpu_clock_begin(id_clock_vertvisc)
397 call vertvisc_coef(up, vp, h_av, forces, visc, dt, g, gv, us, &
398 cs%vertvisc_CSp, cs%OBC)
399 call vertvisc(up, vp, h_av, forces, visc, dt, cs%OBC, cs%ADp, cs%CDp, &
400 g, gv, us, cs%vertvisc_CSp, cs%taux_bot, cs%tauy_bot)
401 call vertvisc_coef(u_in, v_in, h_av, forces, visc, dt, g, gv, us, &
402 cs%vertvisc_CSp, cs%OBC)
403 call vertvisc(u_in, v_in, h_av, forces, visc, dt, cs%OBC, cs%ADp, cs%CDp,&
404 g, gv, us, cs%vertvisc_CSp, cs%taux_bot, cs%tauy_bot)
405 call cpu_clock_end(id_clock_vertvisc)
406 call pass_vector(up, vp, g%Domain, clock=id_clock_pass)
407 call pass_vector(u_in, v_in, g%Domain, clock=id_clock_pass)
411 call cpu_clock_begin(id_clock_continuity)
412 call continuity(up, vp, h_in, h_in, uh, vh, &
413 dt, g, gv, us, cs%continuity_CSp, obc=cs%OBC)
414 call cpu_clock_end(id_clock_continuity)
415 call pass_var(h_in, g%Domain, clock=id_clock_pass)
416 call pass_vector(uh, vh, g%Domain, clock=id_clock_pass)
420 do j=js-2,je+2 ;
do i=isq-2,ieq+2
421 uhtr(i,j,k) = uhtr(i,j,k) + dt*uh(i,j,k)
423 do j=jsq-2,jeq+2 ;
do i=is-2,ie+2
424 vhtr(i,j,k) = vhtr(i,j,k) + dt*vh(i,j,k)
430 call mom_accel_chksum(
"Corrector accel", cs%CAu, cs%CAv, cs%PFu, cs%PFv, &
431 cs%diffu, cs%diffv, g, gv, us)
434 if (gv%Boussinesq)
then
435 do j=js,je ;
do i=is,ie ; eta_av(i,j) = -gv%Z_to_H*g%bathyT(i,j) ;
enddo ;
enddo
437 do j=js,je ;
do i=is,ie ; eta_av(i,j) = 0.0 ;
enddo ;
enddo
439 do k=1,nz ;
do j=js,je ;
do i=is,ie
440 eta_av(i,j) = eta_av(i,j) + h_av(i,j,k)
441 enddo ;
enddo ;
enddo
443 if (dyn_p_surf)
deallocate(p_surf)
447 if (cs%id_PFu > 0)
call post_data(cs%id_PFu, cs%PFu, cs%diag)
448 if (cs%id_PFv > 0)
call post_data(cs%id_PFv, cs%PFv, cs%diag)
449 if (cs%id_CAu > 0)
call post_data(cs%id_CAu, cs%CAu, cs%diag)
450 if (cs%id_CAv > 0)
call post_data(cs%id_CAv, cs%CAv, cs%diag)
453 if (cs%id_uh > 0)
call post_data(cs%id_uh, uh, cs%diag)
454 if (cs%id_vh > 0)
call post_data(cs%id_vh, vh, cs%diag)
456 end subroutine step_mom_dyn_unsplit_rk2
465 subroutine register_restarts_dyn_unsplit_rk2(HI, GV, param_file, CS, restart_CS)
479 character(len=48) :: thickness_units, flux_units
480 integer :: isd, ied, jsd, jed, nz, isdb, iedb, jsdb, jedb
481 isd = hi%isd ; ied = hi%ied ; jsd = hi%jsd ; jed = hi%jed ; nz = gv%ke
482 isdb = hi%IsdB ; iedb = hi%IedB ; jsdb = hi%JsdB ; jedb = hi%JedB
485 if (
associated(cs))
then
486 call mom_error(warning,
"register_restarts_dyn_unsplit_RK2 called with an associated "// &
487 "control structure.")
492 alloc_(cs%diffu(isdb:iedb,jsd:jed,nz)) ; cs%diffu(:,:,:) = 0.0
493 alloc_(cs%diffv(isd:ied,jsdb:jedb,nz)) ; cs%diffv(:,:,:) = 0.0
494 alloc_(cs%CAu(isdb:iedb,jsd:jed,nz)) ; cs%CAu(:,:,:) = 0.0
495 alloc_(cs%CAv(isd:ied,jsdb:jedb,nz)) ; cs%CAv(:,:,:) = 0.0
496 alloc_(cs%PFu(isdb:iedb,jsd:jed,nz)) ; cs%PFu(:,:,:) = 0.0
497 alloc_(cs%PFv(isd:ied,jsdb:jedb,nz)) ; cs%PFv(:,:,:) = 0.0
499 thickness_units = get_thickness_units(gv)
500 flux_units = get_flux_units(gv)
504 end subroutine register_restarts_dyn_unsplit_rk2
507 subroutine initialize_dyn_unsplit_rk2(u, v, h, Time, G, GV, US, param_file, diag, CS, &
508 restart_CS, Accel_diag, Cont_diag, MIS, MEKE, &
509 OBC, update_OBC_CSp, ALE_CSp, setVisc_CSp, &
514 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)),
intent(inout) :: u
515 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)),
intent(inout) :: v
516 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)) ,
intent(inout) :: h
517 type(time_type),
target,
intent(in) :: time
520 type(
diag_ctrl),
target,
intent(inout) :: diag
542 type(
ale_cs),
pointer :: ale_csp
552 integer,
target,
intent(inout) :: ntrunc
560 character(len=40) :: mdl =
"MOM_dynamics_unsplit_RK2"
561 character(len=48) :: thickness_units, flux_units
564 integer :: isd, ied, jsd, jed, nz, isdb, iedb, jsdb, jedb
565 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed ; nz = g%ke
566 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
568 if (.not.
associated(cs))
call mom_error(fatal, &
569 "initialize_dyn_unsplit_RK2 called with an unassociated control structure.")
570 if (cs%module_is_initialized)
then
571 call mom_error(warning,
"initialize_dyn_unsplit_RK2 called with a control "// &
572 "structure that has already been initialized.")
575 cs%module_is_initialized = .true.
579 call get_param(param_file, mdl,
"BE", cs%be, &
580 "If SPLIT is true, BE determines the relative weighting "//&
581 "of a 2nd-order Runga-Kutta baroclinic time stepping "//&
582 "scheme (0.5) and a backward Euler scheme (1) that is "//&
583 "used for the Coriolis and inertial terms. BE may be "//&
584 "from 0.5 to 1, but instability may occur near 0.5. "//&
585 "BE is also applicable if SPLIT is false and USE_RK2 "//&
586 "is true.", units=
"nondim", default=0.6)
587 call get_param(param_file, mdl,
"BEGW", cs%begw, &
588 "If SPLIT is true, BEGW is a number from 0 to 1 that "//&
589 "controls the extent to which the treatment of gravity "//&
590 "waves is forward-backward (0) or simulated backward "//&
591 "Euler (1). 0 is almost always used. "//&
592 "If SPLIT is false and USE_RK2 is true, BEGW can be "//&
593 "between 0 and 0.5 to damp gravity waves.", &
594 units=
"nondim", default=0.0)
595 call get_param(param_file, mdl,
"DEBUG", cs%debug, &
596 "If true, write out verbose debugging data.", &
597 default=.false., debuggingparam=.true.)
598 call get_param(param_file, mdl,
"TIDES", use_tides, &
599 "If true, apply tidal momentum forcing.", default=.false.)
601 allocate(cs%taux_bot(isdb:iedb,jsd:jed)) ; cs%taux_bot(:,:) = 0.0
602 allocate(cs%tauy_bot(isd:ied,jsdb:jedb)) ; cs%tauy_bot(:,:) = 0.0
604 mis%diffu => cs%diffu ; mis%diffv => cs%diffv
605 mis%PFu => cs%PFu ; mis%PFv => cs%PFv
606 mis%CAu => cs%CAu ; mis%CAv => cs%CAv
608 cs%ADp => accel_diag ; cs%CDp => cont_diag
609 accel_diag%diffu => cs%diffu ; accel_diag%diffv => cs%diffv
610 accel_diag%PFu => cs%PFu ; accel_diag%PFv => cs%PFv
611 accel_diag%CAu => cs%CAu ; accel_diag%CAv => cs%CAv
613 call continuity_init(time, g, gv, param_file, diag, cs%continuity_CSp)
614 call coriolisadv_init(time, g, param_file, diag, cs%ADp, cs%CoriolisAdv_CSp)
615 if (use_tides)
call tidal_forcing_init(time, g, param_file, cs%tides_CSp)
616 call pressureforce_init(time, g, gv, us, param_file, diag, cs%PressureForce_CSp, &
618 call hor_visc_init(time, g, us, param_file, diag, cs%hor_visc_CSp, meke)
619 call vertvisc_init(mis, time, g, gv, us, param_file, diag, cs%ADp, dirs, &
620 ntrunc, cs%vertvisc_CSp)
621 if (.not.
associated(setvisc_csp))
call mom_error(fatal, &
622 "initialize_dyn_unsplit_RK2 called with setVisc_CSp unassociated.")
623 cs%set_visc_CSp => setvisc_csp
625 if (
associated(ale_csp)) cs%ALE_CSp => ale_csp
626 if (
associated(obc)) cs%OBC => obc
628 flux_units = get_flux_units(gv)
629 h_convert = gv%H_to_m ;
if (.not.gv%Boussinesq) h_convert = gv%H_to_kg_m2
630 cs%id_uh = register_diag_field(
'ocean_model',
'uh', diag%axesCuL, time, &
631 'Zonal Thickness Flux', flux_units, y_cell_method=
'sum', v_extensive=.true., &
632 conversion=h_convert)
633 cs%id_vh = register_diag_field(
'ocean_model',
'vh', diag%axesCvL, time, &
634 'Meridional Thickness Flux', flux_units, x_cell_method=
'sum', v_extensive=.true., &
635 conversion=h_convert)
636 cs%id_CAu = register_diag_field(
'ocean_model',
'CAu', diag%axesCuL, time, &
637 'Zonal Coriolis and Advective Acceleration',
'meter second-2')
638 cs%id_CAv = register_diag_field(
'ocean_model',
'CAv', diag%axesCvL, time, &
639 'Meridional Coriolis and Advective Acceleration',
'meter second-2')
640 cs%id_PFu = register_diag_field(
'ocean_model',
'PFu', diag%axesCuL, time, &
641 'Zonal Pressure Force Acceleration',
'meter second-2')
642 cs%id_PFv = register_diag_field(
'ocean_model',
'PFv', diag%axesCvL, time, &
643 'Meridional Pressure Force Acceleration',
'meter second-2')
645 id_clock_cor = cpu_clock_id(
'(Ocean Coriolis & mom advection)', grain=clock_module)
646 id_clock_continuity = cpu_clock_id(
'(Ocean continuity equation)', grain=clock_module)
647 id_clock_pres = cpu_clock_id(
'(Ocean pressure force)', grain=clock_module)
648 id_clock_vertvisc = cpu_clock_id(
'(Ocean vertical viscosity)', grain=clock_module)
649 id_clock_horvisc = cpu_clock_id(
'(Ocean horizontal viscosity)', grain=clock_module)
650 id_clock_mom_update = cpu_clock_id(
'(Ocean momentum increments)', grain=clock_module)
651 id_clock_pass = cpu_clock_id(
'(Ocean message passing)', grain=clock_module)
652 id_clock_pass_init = cpu_clock_id(
'(Ocean init message passing)', grain=clock_routine)
654 end subroutine initialize_dyn_unsplit_rk2
657 subroutine end_dyn_unsplit_rk2(CS)
661 dealloc_(cs%diffu) ; dealloc_(cs%diffv)
662 dealloc_(cs%CAu) ; dealloc_(cs%CAv)
663 dealloc_(cs%PFu) ; dealloc_(cs%PFv)
666 end subroutine end_dyn_unsplit_rk2