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
74 use mom_time_manager,
only :
operator(-),
operator(>),
operator(*),
operator(/)
99 implicit none ;
private
101 #include <MOM_memory.h>
105 real allocable_,
dimension(NIMEMB_PTR_,NJMEM_,NKMEM_) :: &
106 cau, & !< CAu = f*v - u.grad(u) [L T-2 ~> m s-2].
107 pfu, & !< PFu = -dM/dx [L T-2 ~> m s-2].
110 real allocable_,
dimension(NIMEM_,NJMEMB_PTR_,NKMEM_) :: &
111 cav, & !< CAv = -f*u - u.grad(v) [L T-2 ~> m s-2].
112 pfv, & !< PFv = -dM/dy [L T-2 ~> m s-2].
115 real,
pointer,
dimension(:,:) :: taux_bot => null()
117 real,
pointer,
dimension(:,:) :: tauy_bot => null()
128 logical :: module_is_initialized = .false.
131 integer :: id_uh = -1, id_vh = -1
132 integer :: id_pfu = -1, id_pfv = -1, id_cau = -1, id_cav = -1
162 type(
ale_cs),
pointer :: ale_csp => null()
189 p_surf_begin, p_surf_end, uh, vh, uhtr, vhtr, eta_av, G, GV, US, CS, &
194 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)),
intent(inout) :: u_in
196 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)),
intent(inout) :: v_in
198 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(inout) :: h_in
206 type(time_type),
intent(in) :: time_local
208 real,
intent(in) :: dt
210 real,
dimension(:,:),
pointer :: p_surf_begin
213 real,
dimension(:,:),
pointer :: p_surf_end
216 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)),
intent(inout) :: uh
218 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)),
intent(inout) :: vh
220 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)),
intent(inout) :: uhtr
223 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)),
intent(inout) :: vhtr
226 real,
dimension(SZI_(G),SZJ_(G)),
intent(out) :: eta_av
237 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)) :: h_av, hp
238 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)) :: up
239 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)) :: vp
240 real,
dimension(:,:),
pointer :: p_surf => null()
243 logical :: dyn_p_surf
244 integer :: i, j, k, is, ie, js, je, isq, ieq, jsq, jeq, nz
245 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
246 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
249 h_av(:,:,:) = 0; hp(:,:,:) = 0
253 dyn_p_surf =
associated(p_surf_begin) .and.
associated(p_surf_end)
255 call safe_alloc_ptr(p_surf,g%isd,g%ied,g%jsd,g%jed) ; p_surf(:,:) = 0.0
257 p_surf => forces%p_surf
264 call mom_state_chksum(
"Start Predictor ", u_in, v_in, h_in, uh, vh, g, gv, us)
268 call enable_averages(dt,time_local, cs%diag)
271 g, gv, us, cs%hor_visc_CSp)
273 call disable_averaging(cs%diag)
283 call continuity(u_in, v_in, h_in, hp, uh, vh, dt_pred, g, gv, us, &
284 cs%continuity_CSp, obc=cs%OBC)
292 do j=js-2,je+2 ;
do i=is-2,ie+2
293 h_av(i,j,k) = (h_in(i,j,k) + hp(i,j,k)) * 0.5
294 enddo ;
enddo ;
enddo
299 call coradcalc(u_in, v_in, h_av, uh, vh, cs%CAu, cs%CAv, cs%OBC, cs%ADp, &
300 g, gv, us, cs%CoriolisAdv_CSp)
305 if (dyn_p_surf)
then ;
do j=js-2,je+2 ;
do i=is-2,ie+2
306 p_surf(i,j) = 0.5*p_surf_begin(i,j) + 0.5*p_surf_end(i,j)
307 enddo ;
enddo ;
endif
309 cs%PressureForce_CSp, cs%ALE_CSp, p_surf)
314 if (
associated(cs%OBC)) then;
if (cs%OBC%update_OBC)
then
315 call update_obc_data(cs%OBC, g, gv, us, tv, h_in, cs%update_OBC_CSp, time_local)
317 if (
associated(cs%OBC))
then
325 do k=1,nz ;
do j=js,je ;
do i=isq,ieq
326 up(i,j,k) = g%mask2dCu(i,j) * (u_in(i,j,k) + dt_pred * &
327 ((cs%PFu(i,j,k) + cs%CAu(i,j,k)) + cs%diffu(i,j,k)))
328 enddo ;
enddo ;
enddo
329 do k=1,nz ;
do j=jsq,jeq ;
do i=is,ie
330 vp(i,j,k) = g%mask2dCv(i,j) * (v_in(i,j,k) + dt_pred * &
331 ((cs%PFv(i,j,k) + cs%CAv(i,j,k)) + cs%diffv(i,j,k)))
332 enddo ;
enddo ;
enddo
337 cs%diffu, cs%diffv, g, gv, us)
341 call enable_averages(dt, time_local, cs%diag)
342 call set_viscous_ml(up, vp, h_av, tv, forces, visc, dt_pred, g, gv, us, &
344 call disable_averaging(cs%diag)
345 call vertvisc_coef(up, vp, h_av, forces, visc, dt_pred, g, gv, us, &
346 cs%vertvisc_CSp, cs%OBC)
347 call vertvisc(up, vp, h_av, forces, visc, dt_pred, cs%OBC, cs%ADp, cs%CDp, &
348 g, gv, us, cs%vertvisc_CSp)
355 call continuity(up, vp, h_in, hp, uh, vh, dt, g, gv, us, cs%continuity_CSp, obc=cs%OBC)
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
370 call coradcalc(up, vp, h_av, uh, vh, cs%CAu, cs%CAv, cs%OBC, cs%ADp, &
371 g, gv, us, cs%CoriolisAdv_CSp)
373 if (
associated(cs%OBC))
then
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)) + 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)) + 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)) + 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)) + cs%diffv(i,j,k)))
392 enddo ;
enddo ;
enddo
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)
412 call continuity(up, vp, h_in, h_in, uh, vh, dt, g, gv, us, cs%continuity_CSp, obc=cs%OBC)
419 do j=js-2,je+2 ;
do i=isq-2,ieq+2
420 uhtr(i,j,k) = uhtr(i,j,k) + dt*uh(i,j,k)
422 do j=jsq-2,jeq+2 ;
do i=is-2,ie+2
423 vhtr(i,j,k) = vhtr(i,j,k) + dt*vh(i,j,k)
430 cs%diffu, cs%diffv, g, gv, us)
433 if (gv%Boussinesq)
then
434 do j=js,je ;
do i=is,ie ; eta_av(i,j) = -gv%Z_to_H*g%bathyT(i,j) ;
enddo ;
enddo
436 do j=js,je ;
do i=is,ie ; eta_av(i,j) = 0.0 ;
enddo ;
enddo
438 do k=1,nz ;
do j=js,je ;
do i=is,ie
439 eta_av(i,j) = eta_av(i,j) + h_av(i,j,k)
440 enddo ;
enddo ;
enddo
442 if (dyn_p_surf)
deallocate(p_surf)
446 if (cs%id_PFu > 0)
call post_data(cs%id_PFu, cs%PFu, cs%diag)
447 if (cs%id_PFv > 0)
call post_data(cs%id_PFv, cs%PFv, cs%diag)
448 if (cs%id_CAu > 0)
call post_data(cs%id_CAu, cs%CAu, cs%diag)
449 if (cs%id_CAv > 0)
call post_data(cs%id_CAv, cs%CAv, cs%diag)
452 if (cs%id_uh > 0)
call post_data(cs%id_uh, uh, cs%diag)
453 if (cs%id_vh > 0)
call post_data(cs%id_vh, vh, cs%diag)
478 character(len=48) :: thickness_units, flux_units
479 integer :: isd, ied, jsd, jed, nz, isdb, iedb, jsdb, jedb
480 isd = hi%isd ; ied = hi%ied ; jsd = hi%jsd ; jed = hi%jed ; nz = gv%ke
481 isdb = hi%IsdB ; iedb = hi%IedB ; jsdb = hi%JsdB ; jedb = hi%JedB
484 if (
associated(cs))
then
485 call mom_error(warning,
"register_restarts_dyn_unsplit_RK2 called with an associated "// &
486 "control structure.")
491 alloc_(cs%diffu(isdb:iedb,jsd:jed,nz)) ; cs%diffu(:,:,:) = 0.0
492 alloc_(cs%diffv(isd:ied,jsdb:jedb,nz)) ; cs%diffv(:,:,:) = 0.0
493 alloc_(cs%CAu(isdb:iedb,jsd:jed,nz)) ; cs%CAu(:,:,:) = 0.0
494 alloc_(cs%CAv(isd:ied,jsdb:jedb,nz)) ; cs%CAv(:,:,:) = 0.0
495 alloc_(cs%PFu(isdb:iedb,jsd:jed,nz)) ; cs%PFu(:,:,:) = 0.0
496 alloc_(cs%PFv(isd:ied,jsdb:jedb,nz)) ; cs%PFv(:,:,:) = 0.0
498 thickness_units = get_thickness_units(gv)
506 subroutine initialize_dyn_unsplit_rk2(u, v, h, Time, G, GV, US, param_file, diag, CS, &
507 restart_CS, Accel_diag, Cont_diag, MIS, MEKE, &
508 OBC, update_OBC_CSp, ALE_CSp, setVisc_CSp, &
513 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)),
intent(inout) :: u
514 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)),
intent(inout) :: v
515 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)) ,
intent(inout) :: h
516 type(time_type),
target,
intent(in) :: time
519 type(
diag_ctrl),
target,
intent(inout) :: diag
541 type(
ale_cs),
pointer :: ale_csp
551 integer,
target,
intent(inout) :: ntrunc
559 character(len=40) :: mdl =
"MOM_dynamics_unsplit_RK2"
560 character(len=48) :: thickness_units, flux_units
563 integer :: isd, ied, jsd, jed, nz, isdb, iedb, jsdb, jedb
564 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed ; nz = g%ke
565 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
567 if (.not.
associated(cs))
call mom_error(fatal, &
568 "initialize_dyn_unsplit_RK2 called with an unassociated control structure.")
569 if (cs%module_is_initialized)
then
570 call mom_error(warning,
"initialize_dyn_unsplit_RK2 called with a control "// &
571 "structure that has already been initialized.")
574 cs%module_is_initialized = .true.
578 call get_param(param_file, mdl,
"BE", cs%be, &
579 "If SPLIT is true, BE determines the relative weighting "//&
580 "of a 2nd-order Runga-Kutta baroclinic time stepping "//&
581 "scheme (0.5) and a backward Euler scheme (1) that is "//&
582 "used for the Coriolis and inertial terms. BE may be "//&
583 "from 0.5 to 1, but instability may occur near 0.5. "//&
584 "BE is also applicable if SPLIT is false and USE_RK2 "//&
585 "is true.", units=
"nondim", default=0.6)
586 call get_param(param_file, mdl,
"BEGW", cs%begw, &
587 "If SPLIT is true, BEGW is a number from 0 to 1 that "//&
588 "controls the extent to which the treatment of gravity "//&
589 "waves is forward-backward (0) or simulated backward "//&
590 "Euler (1). 0 is almost always used. "//&
591 "If SPLIT is false and USE_RK2 is true, BEGW can be "//&
592 "between 0 and 0.5 to damp gravity waves.", &
593 units=
"nondim", default=0.0)
594 call get_param(param_file, mdl,
"DEBUG", cs%debug, &
595 "If true, write out verbose debugging data.", &
596 default=.false., debuggingparam=.true.)
597 call get_param(param_file, mdl,
"TIDES", use_tides, &
598 "If true, apply tidal momentum forcing.", default=.false.)
600 allocate(cs%taux_bot(isdb:iedb,jsd:jed)) ; cs%taux_bot(:,:) = 0.0
601 allocate(cs%tauy_bot(isd:ied,jsdb:jedb)) ; cs%tauy_bot(:,:) = 0.0
603 mis%diffu => cs%diffu ; mis%diffv => cs%diffv
604 mis%PFu => cs%PFu ; mis%PFv => cs%PFv
605 mis%CAu => cs%CAu ; mis%CAv => cs%CAv
607 cs%ADp => accel_diag ; cs%CDp => cont_diag
608 accel_diag%diffu => cs%diffu ; accel_diag%diffv => cs%diffv
609 accel_diag%PFu => cs%PFu ; accel_diag%PFv => cs%PFv
610 accel_diag%CAu => cs%CAu ; accel_diag%CAv => cs%CAv
612 call continuity_init(time, g, gv, us, param_file, diag, cs%continuity_CSp)
613 call coriolisadv_init(time, g, gv, us, param_file, diag, cs%ADp, cs%CoriolisAdv_CSp)
617 call hor_visc_init(time, g, us, param_file, diag, cs%hor_visc_CSp, meke)
618 call vertvisc_init(mis, time, g, gv, us, param_file, diag, cs%ADp, dirs, &
619 ntrunc, cs%vertvisc_CSp)
620 if (.not.
associated(setvisc_csp))
call mom_error(fatal, &
621 "initialize_dyn_unsplit_RK2 called with setVisc_CSp unassociated.")
622 cs%set_visc_CSp => setvisc_csp
624 if (
associated(ale_csp)) cs%ALE_CSp => ale_csp
625 if (
associated(obc)) cs%OBC => obc
628 h_convert = gv%H_to_m ;
if (.not.gv%Boussinesq) h_convert = gv%H_to_kg_m2
629 cs%id_uh = register_diag_field(
'ocean_model',
'uh', diag%axesCuL, time, &
630 'Zonal Thickness Flux', flux_units, y_cell_method=
'sum', v_extensive=.true., &
631 conversion=h_convert*us%L_to_m**2*us%s_to_T)
632 cs%id_vh = register_diag_field(
'ocean_model',
'vh', diag%axesCvL, time, &
633 'Meridional Thickness Flux', flux_units, x_cell_method=
'sum', v_extensive=.true., &
634 conversion=h_convert*us%L_to_m**2*us%s_to_T)
635 cs%id_CAu = register_diag_field(
'ocean_model',
'CAu', diag%axesCuL, time, &
636 'Zonal Coriolis and Advective Acceleration',
'meter second-2', conversion=us%L_T2_to_m_s2)
637 cs%id_CAv = register_diag_field(
'ocean_model',
'CAv', diag%axesCvL, time, &
638 'Meridional Coriolis and Advective Acceleration',
'meter second-2', conversion=us%L_T2_to_m_s2)
639 cs%id_PFu = register_diag_field(
'ocean_model',
'PFu', diag%axesCuL, time, &
640 'Zonal Pressure Force Acceleration',
'meter second-2', conversion=us%L_T2_to_m_s2)
641 cs%id_PFv = register_diag_field(
'ocean_model',
'PFv', diag%axesCvL, time, &
642 'Meridional Pressure Force Acceleration',
'meter second-2', conversion=us%L_T2_to_m_s2)
644 id_clock_cor = cpu_clock_id(
'(Ocean Coriolis & mom advection)', grain=clock_module)
646 id_clock_pres = cpu_clock_id(
'(Ocean pressure force)', grain=clock_module)
647 id_clock_vertvisc = cpu_clock_id(
'(Ocean vertical viscosity)', grain=clock_module)
648 id_clock_horvisc = cpu_clock_id(
'(Ocean horizontal viscosity)', grain=clock_module)
650 id_clock_pass = cpu_clock_id(
'(Ocean message passing)', grain=clock_module)
651 id_clock_pass_init = cpu_clock_id(
'(Ocean init message passing)', grain=clock_routine)
660 dealloc_(cs%diffu) ; dealloc_(cs%diffv)
661 dealloc_(cs%CAu) ; dealloc_(cs%CAv)
662 dealloc_(cs%PFu) ; dealloc_(cs%PFv)