57 use mom_cpu_clock,
only : cpu_clock_id, cpu_clock_begin, cpu_clock_end
59 use mom_cpu_clock,
only : clock_module_driver, clock_module, clock_routine
75 use mom_time_manager,
only :
operator(-),
operator(>),
operator(*),
operator(/)
102 implicit none ;
private
104 #include <MOM_memory.h>
108 real allocable_,
dimension(NIMEMB_PTR_,NJMEM_,NKMEM_) :: &
109 cau, & !< CAu = f*v - u.grad(u) [L T-2 ~> m s-2].
110 pfu, & !< PFu = -dM/dx [L T-2 ~> m s-2].
113 real allocable_,
dimension(NIMEM_,NJMEMB_PTR_,NKMEM_) :: &
114 cav, & !< CAv = -f*u - u.grad(v) [L T-2 ~> m s-2].
115 pfv, & !< PFv = -dM/dy [L T-2 ~> m s-2].
118 real,
pointer,
dimension(:,:) :: taux_bot => null()
120 real,
pointer,
dimension(:,:) :: tauy_bot => null()
125 logical :: module_is_initialized = .false.
128 integer :: id_uh = -1, id_vh = -1
129 integer :: id_pfu = -1, id_pfv = -1, id_cau = -1, id_cav = -1
159 type(
ale_cs),
pointer :: ale_csp => null()
186 p_surf_begin, p_surf_end, uh, vh, uhtr, vhtr, eta_av, G, GV, US, CS, &
191 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)),
intent(inout) :: u
192 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)),
intent(inout) :: v
193 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(inout) :: h
198 type(time_type),
intent(in) :: time_local
200 real,
intent(in) :: dt
202 real,
dimension(:,:),
pointer :: p_surf_begin
204 real,
dimension(:,:),
pointer :: p_surf_end
206 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)),
intent(inout) :: uh
208 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)),
intent(inout) :: vh
210 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)),
intent(inout) :: uhtr
212 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)),
intent(inout) :: vhtr
214 real,
dimension(SZI_(G),SZJ_(G)),
intent(out) :: eta_av
226 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)) :: h_av, hp
227 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)) :: up, upp
228 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)) :: vp, vpp
229 real,
dimension(:,:),
pointer :: p_surf => null()
231 logical :: dyn_p_surf
232 integer :: i, j, k, is, ie, js, je, isq, ieq, jsq, jeq, nz
233 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
234 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
237 h_av(:,:,:) = 0; hp(:,:,:) = 0
238 up(:,:,:) = 0; upp(:,:,:) = 0
239 vp(:,:,:) = 0; vpp(:,:,:) = 0
241 dyn_p_surf =
associated(p_surf_begin) .and.
associated(p_surf_end)
243 call safe_alloc_ptr(p_surf,g%isd,g%ied,g%jsd,g%jed) ; p_surf(:,:) = 0.0
245 p_surf => forces%p_surf
256 call enable_averages(dt, time_local, cs%diag)
259 g, gv, us, cs%hor_visc_CSp)
261 call disable_averaging(cs%diag)
266 call continuity(u, v, h, hp, uh, vh, dt*0.5, g, gv, us, cs%continuity_CSp, obc=cs%OBC)
271 call enable_averages(0.5*dt, time_local-
real_to_time(0.5*us%T_to_s*dt), cs%diag)
273 if (cs%id_uh > 0)
call post_data(cs%id_uh, uh, cs%diag)
274 if (cs%id_vh > 0)
call post_data(cs%id_vh, vh, cs%diag)
275 call disable_averaging(cs%diag)
281 do j=js-2,je+2 ;
do i=is-2,ie+2
282 h_av(i,j,k) = (h(i,j,k) + hp(i,j,k)) * 0.5
284 do j=js,je ;
do i=isq,ieq
285 u(i,j,k) = u(i,j,k) + dt * cs%diffu(i,j,k) * g%mask2dCu(i,j)
287 do j=jsq,jeq ;
do i=is,ie
288 v(i,j,k) = v(i,j,k) + dt * cs%diffv(i,j,k) * g%mask2dCv(i,j)
290 do j=js-2,je+2 ;
do i=isq-2,ieq+2
291 uhtr(i,j,k) = uhtr(i,j,k) + 0.5*dt*uh(i,j,k)
293 do j=jsq-2,jeq+2 ;
do i=is-2,ie+2
294 vhtr(i,j,k) = vhtr(i,j,k) + 0.5*dt*vh(i,j,k)
302 call coradcalc(u, v, h_av, uh, vh, cs%CAu, cs%CAv, cs%OBC, cs%ADp, &
303 g, gv, us, cs%CoriolisAdv_CSp)
308 if (dyn_p_surf)
then ;
do j=js-2,je+2 ;
do i=is-2,ie+2
309 p_surf(i,j) = 0.75*p_surf_begin(i,j) + 0.25*p_surf_end(i,j)
310 enddo ;
enddo ;
endif
312 cs%PressureForce_CSp, cs%ALE_CSp, p_surf)
315 if (
associated(cs%OBC)) then;
if (cs%OBC%update_OBC)
then
316 call update_obc_data(cs%OBC, g, gv, us, tv, h, cs%update_OBC_CSp, time_local)
318 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(i,j,k) + dt_pred * &
327 (cs%PFu(i,j,k) + cs%CAu(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(i,j,k) + dt_pred * &
331 (cs%PFv(i,j,k) + cs%CAv(i,j,k)))
332 enddo ;
enddo ;
enddo
338 cs%diffu, cs%diffv, g, gv, us)
343 call enable_averages(dt, time_local, cs%diag)
344 call set_viscous_ml(u, v, h_av, tv, forces, visc, dt*0.5, g, gv, us, &
346 call disable_averaging(cs%diag)
348 call vertvisc_coef(up, vp, h_av, forces, visc, dt*0.5, g, gv, us, &
349 cs%vertvisc_CSp, cs%OBC)
350 call vertvisc(up, vp, h_av, forces, visc, dt*0.5, cs%OBC, cs%ADp, cs%CDp, &
351 g, gv, us, cs%vertvisc_CSp, waves=waves)
358 call continuity(up, vp, hp, h_av, uh, vh, (0.5*dt), g, gv, us, &
359 cs%continuity_CSp, obc=cs%OBC)
365 do k=1,nz ;
do j=js-2,je+2 ;
do i=is-2,ie+2
366 h_av(i,j,k) = (hp(i,j,k) + h_av(i,j,k)) * 0.5
367 enddo ;
enddo ;
enddo
371 call coradcalc(up, vp, h_av, uh, vh, cs%CAu, cs%CAv, cs%OBC, cs%ADp, &
372 g, gv, us, cs%CoriolisAdv_CSp)
377 if (dyn_p_surf)
then ;
do j=js-2,je+2 ;
do i=is-2,ie+2
378 p_surf(i,j) = 0.25*p_surf_begin(i,j) + 0.75*p_surf_end(i,j)
379 enddo ;
enddo ;
endif
381 cs%PressureForce_CSp, cs%ALE_CSp, p_surf)
384 if (
associated(cs%OBC)) then;
if (cs%OBC%update_OBC)
then
385 call update_obc_data(cs%OBC, g, gv, us, tv, h, cs%update_OBC_CSp, time_local)
387 if (
associated(cs%OBC))
then
394 do k=1,nz ;
do j=js,je ;
do i=isq,ieq
395 upp(i,j,k) = g%mask2dCu(i,j) * (u(i,j,k) + dt * 0.5 * &
396 (cs%PFu(i,j,k) + cs%CAu(i,j,k)))
397 enddo ;
enddo ;
enddo
398 do k=1,nz ;
do j=jsq,jeq ;
do i=is,ie
399 vpp(i,j,k) = g%mask2dCv(i,j) * (v(i,j,k) + dt * 0.5 * &
400 (cs%PFv(i,j,k) + cs%CAv(i,j,k)))
401 enddo ;
enddo ;
enddo
407 cs%diffu, cs%diffv, g, gv, us)
412 call vertvisc_coef(upp, vpp, hp, forces, visc, dt*0.5, g, gv, us, &
413 cs%vertvisc_CSp, cs%OBC)
414 call vertvisc(upp, vpp, hp, forces, visc, dt*0.5, cs%OBC, cs%ADp, cs%CDp, &
415 g, gv, us, cs%vertvisc_CSp, waves=waves)
422 call continuity(upp, vpp, hp, h, uh, vh, (dt*0.5), g, gv, us, &
423 cs%continuity_CSp, obc=cs%OBC)
431 call enable_averages(0.5*dt, time_local, cs%diag)
433 if (cs%id_uh > 0)
call post_data(cs%id_uh, uh, cs%diag)
434 if (cs%id_vh > 0)
call post_data(cs%id_vh, vh, cs%diag)
435 call disable_averaging(cs%diag)
436 call enable_averages(dt, time_local, cs%diag)
440 do j=js-2,je+2 ;
do i=is-2,ie+2
441 h_av(i,j,k) = 0.5*(h(i,j,k) + hp(i,j,k))
443 do j=js-2,je+2 ;
do i=isq-2,ieq+2
444 uhtr(i,j,k) = uhtr(i,j,k) + 0.5*dt*uh(i,j,k)
446 do j=jsq-2,jeq+2 ;
do i=is-2,ie+2
447 vhtr(i,j,k) = vhtr(i,j,k) + 0.5*dt*vh(i,j,k)
453 call coradcalc(upp, vpp, h_av, uh, vh, cs%CAu, cs%CAv, cs%OBC, cs%ADp, &
454 g, gv, us, cs%CoriolisAdv_CSp)
460 cs%PressureForce_CSp, cs%ALE_CSp, p_surf)
463 if (
associated(cs%OBC)) then;
if (cs%OBC%update_OBC)
then
464 call update_obc_data(cs%OBC, g, gv, us, tv, h, cs%update_OBC_CSp, time_local)
468 if (
associated(cs%OBC))
then
472 do k=1,nz ;
do j=js,je ;
do i=isq,ieq
473 u(i,j,k) = g%mask2dCu(i,j) * (u(i,j,k) + dt * &
474 (cs%PFu(i,j,k) + cs%CAu(i,j,k)))
475 enddo ;
enddo ;
enddo
476 do k=1,nz ;
do j=jsq,jeq ;
do i=is,ie
477 v(i,j,k) = g%mask2dCv(i,j) * (v(i,j,k) + dt * &
478 (cs%PFv(i,j,k) + cs%CAv(i,j,k)))
479 enddo ;
enddo ;
enddo
483 call vertvisc_coef(u, v, h_av, forces, visc, dt, g, gv, us, cs%vertvisc_CSp, cs%OBC)
484 call vertvisc(u, v, h_av, forces, visc, dt, cs%OBC, cs%ADp, cs%CDp, &
485 g, gv, us, cs%vertvisc_CSp, cs%taux_bot, cs%tauy_bot, waves=waves)
492 cs%diffu, cs%diffv, g, gv, us)
495 if (gv%Boussinesq)
then
496 do j=js,je ;
do i=is,ie ; eta_av(i,j) = -gv%Z_to_H*g%bathyT(i,j) ;
enddo ;
enddo
498 do j=js,je ;
do i=is,ie ; eta_av(i,j) = 0.0 ;
enddo ;
enddo
500 do k=1,nz ;
do j=js,je ;
do i=is,ie
501 eta_av(i,j) = eta_av(i,j) + h_av(i,j,k)
502 enddo ;
enddo ;
enddo
504 if (dyn_p_surf)
deallocate(p_surf)
508 if (cs%id_PFu > 0)
call post_data(cs%id_PFu, cs%PFu, cs%diag)
509 if (cs%id_PFv > 0)
call post_data(cs%id_PFv, cs%PFv, cs%diag)
510 if (cs%id_CAu > 0)
call post_data(cs%id_CAu, cs%CAu, cs%diag)
511 if (cs%id_CAv > 0)
call post_data(cs%id_CAv, cs%CAv, cs%diag)
532 character(len=40) :: mdl =
"MOM_dynamics_unsplit"
533 character(len=48) :: thickness_units, flux_units
534 integer :: isd, ied, jsd, jed, nz, isdb, iedb, jsdb, jedb
535 isd = hi%isd ; ied = hi%ied ; jsd = hi%jsd ; jed = hi%jed ; nz = gv%ke
536 isdb = hi%IsdB ; iedb = hi%IedB ; jsdb = hi%JsdB ; jedb = hi%JedB
539 if (
associated(cs))
then
540 call mom_error(warning,
"register_restarts_dyn_unsplit called with an associated "// &
541 "control structure.")
546 alloc_(cs%diffu(isdb:iedb,jsd:jed,nz)) ; cs%diffu(:,:,:) = 0.0
547 alloc_(cs%diffv(isd:ied,jsdb:jedb,nz)) ; cs%diffv(:,:,:) = 0.0
548 alloc_(cs%CAu(isdb:iedb,jsd:jed,nz)) ; cs%CAu(:,:,:) = 0.0
549 alloc_(cs%CAv(isd:ied,jsdb:jedb,nz)) ; cs%CAv(:,:,:) = 0.0
550 alloc_(cs%PFu(isdb:iedb,jsd:jed,nz)) ; cs%PFu(:,:,:) = 0.0
551 alloc_(cs%PFv(isd:ied,jsdb:jedb,nz)) ; cs%PFv(:,:,:) = 0.0
553 thickness_units = get_thickness_units(gv)
561 subroutine initialize_dyn_unsplit(u, v, h, Time, G, GV, US, param_file, diag, CS, &
562 restart_CS, Accel_diag, Cont_diag, MIS, MEKE, &
563 OBC, update_OBC_CSp, ALE_CSp, setVisc_CSp, &
568 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
570 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)), &
572 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)) , &
574 type(time_type),
target,
intent(in) :: time
577 type(
diag_ctrl),
target,
intent(inout) :: diag
599 type(
ale_cs),
pointer :: ale_csp
608 integer,
target,
intent(inout) :: ntrunc
616 character(len=40) :: mdl =
"MOM_dynamics_unsplit"
617 character(len=48) :: thickness_units, flux_units
620 integer :: isd, ied, jsd, jed, nz, isdb, iedb, jsdb, jedb
621 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed ; nz = g%ke
622 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
624 if (.not.
associated(cs))
call mom_error(fatal, &
625 "initialize_dyn_unsplit called with an unassociated control structure.")
626 if (cs%module_is_initialized)
then
627 call mom_error(warning,
"initialize_dyn_unsplit called with a control "// &
628 "structure that has already been initialized.")
631 cs%module_is_initialized = .true.
635 call get_param(param_file, mdl,
"DEBUG", cs%debug, &
636 "If true, write out verbose debugging data.", &
637 default=.false., debuggingparam=.true.)
638 call get_param(param_file, mdl,
"TIDES", use_tides, &
639 "If true, apply tidal momentum forcing.", default=.false.)
641 allocate(cs%taux_bot(isdb:iedb,jsd:jed)) ; cs%taux_bot(:,:) = 0.0
642 allocate(cs%tauy_bot(isd:ied,jsdb:jedb)) ; cs%tauy_bot(:,:) = 0.0
644 mis%diffu => cs%diffu ; mis%diffv => cs%diffv
645 mis%PFu => cs%PFu ; mis%PFv => cs%PFv
646 mis%CAu => cs%CAu ; mis%CAv => cs%CAv
648 cs%ADp => accel_diag ; cs%CDp => cont_diag
649 accel_diag%diffu => cs%diffu ; accel_diag%diffv => cs%diffv
650 accel_diag%PFu => cs%PFu ; accel_diag%PFv => cs%PFv
651 accel_diag%CAu => cs%CAu ; accel_diag%CAv => cs%CAv
653 call continuity_init(time, g, gv, us, param_file, diag, cs%continuity_CSp)
654 call coriolisadv_init(time, g, gv, us, param_file, diag, cs%ADp, cs%CoriolisAdv_CSp)
658 call hor_visc_init(time, g, us, param_file, diag, cs%hor_visc_CSp, meke)
659 call vertvisc_init(mis, time, g, gv, us, param_file, diag, cs%ADp, dirs, &
660 ntrunc, cs%vertvisc_CSp)
661 if (.not.
associated(setvisc_csp))
call mom_error(fatal, &
662 "initialize_dyn_unsplit called with setVisc_CSp unassociated.")
663 cs%set_visc_CSp => setvisc_csp
665 if (
associated(ale_csp)) cs%ALE_CSp => ale_csp
666 if (
associated(obc)) cs%OBC => obc
667 if (
associated(update_obc_csp)) cs%update_OBC_CSp => update_obc_csp
670 h_convert = gv%H_to_m ;
if (.not.gv%Boussinesq) h_convert = gv%H_to_kg_m2
671 cs%id_uh = register_diag_field(
'ocean_model',
'uh', diag%axesCuL, time, &
672 'Zonal Thickness Flux', flux_units, y_cell_method=
'sum', v_extensive=.true., &
673 conversion=h_convert*us%L_to_m**2*us%s_to_T)
674 cs%id_vh = register_diag_field(
'ocean_model',
'vh', diag%axesCvL, time, &
675 'Meridional Thickness Flux', flux_units, x_cell_method=
'sum', v_extensive=.true., &
676 conversion=h_convert*us%L_to_m**2*us%s_to_T)
677 cs%id_CAu = register_diag_field(
'ocean_model',
'CAu', diag%axesCuL, time, &
678 'Zonal Coriolis and Advective Acceleration',
'm s-2', &
679 conversion=us%L_T2_to_m_s2)
680 cs%id_CAv = register_diag_field(
'ocean_model',
'CAv', diag%axesCvL, time, &
681 'Meridional Coriolis and Advective Acceleration',
'm s-2', &
682 conversion=us%L_T2_to_m_s2)
683 cs%id_PFu = register_diag_field(
'ocean_model',
'PFu', diag%axesCuL, time, &
684 'Zonal Pressure Force Acceleration',
'm s-2', &
685 conversion=us%L_T2_to_m_s2)
686 cs%id_PFv = register_diag_field(
'ocean_model',
'PFv', diag%axesCvL, time, &
687 'Meridional Pressure Force Acceleration',
'm s-2', &
688 conversion=us%L_T2_to_m_s2)
690 id_clock_cor = cpu_clock_id(
'(Ocean Coriolis & mom advection)', grain=clock_module)
692 id_clock_pres = cpu_clock_id(
'(Ocean pressure force)', grain=clock_module)
693 id_clock_vertvisc = cpu_clock_id(
'(Ocean vertical viscosity)', grain=clock_module)
694 id_clock_horvisc = cpu_clock_id(
'(Ocean horizontal viscosity)', grain=clock_module)
696 id_clock_pass = cpu_clock_id(
'(Ocean message passing)', grain=clock_module)
697 id_clock_pass_init = cpu_clock_id(
'(Ocean init message passing)', grain=clock_routine)
706 dealloc_(cs%diffu) ; dealloc_(cs%diffv)
707 dealloc_(cs%CAu) ; dealloc_(cs%CAv)
708 dealloc_(cs%PFu) ; dealloc_(cs%PFv)