7 use mom_debugging,
only : mom_debugging_init, hchksum, uvchksum
11 use mom_cpu_clock,
only : cpu_clock_id, cpu_clock_begin, cpu_clock_end
13 use mom_cpu_clock,
only : clock_module_driver, clock_module, clock_routine
23 use mom_diag_mediator,
only : diag_copy_storage_to_diag, diag_copy_diag_to_storage
26 use mom_domains,
only : to_north, to_east, to_south, to_west
27 use mom_domains,
only : to_all, omit_corners, cgrid_ne, scalar_pair
29 use mom_domains,
only : start_group_pass, complete_group_pass, omit_corners
43 use mom_time_manager,
only : time_type, real_to_time, time_type_to_real,
operator(+)
44 use mom_time_manager,
only :
operator(-),
operator(>),
operator(*),
operator(/)
47 use coupler_types_mod,
only : coupler_type_send_data, coupler_1d_bc_type, coupler_type_spawn
50 use mom_ale,
only : ale_init, ale_end, ale_main,
ale_cs, adjustgridforintegrity
51 use mom_ale,
only : ale_getcoordinate, ale_getcoordinateunits, ale_writecoordinatefile
52 use mom_ale,
only : ale_updateverticalgridtype, ale_remap_init_conds, ale_register_diags
58 use mom_diagnostics,
only : calculate_diagnostic_fields, mom_diagnostics_init
59 use mom_diagnostics,
only : register_transport_diags, post_transport_diagnostics
61 use mom_diagnostics,
only : post_surface_dyn_diags, post_surface_thermo_diags
76 use mom_grid,
only : set_first_direction, rescale_grid_bathymetry
81 use mom_meke,
only : meke_init, meke_alloc_register_restart, step_forward_meke,
meke_cs
89 use mom_set_visc,
only : set_viscous_bbl, set_viscous_ml, set_visc_init
117 use mom_verticalgrid,
only : get_thickness_units, get_flux_units, get_tr_flux_units
123 use mom_oda_driver_mod,
only : set_prior_tracer, set_analysis_time, apply_oda_tracer_increments
126 use mom_offline_main,
only : insert_offline_main, extract_offline_main, post_offline_convergence_diags
127 use mom_offline_main,
only : register_diags_offline_transport, offline_advection_ale
128 use mom_offline_main,
only : offline_redistribute_residual, offline_diabatic_ale
129 use mom_offline_main,
only : offline_fw_fluxes_into_ocean, offline_fw_fluxes_out_ocean
131 use mom_ale,
only : ale_offline_tracer_final, ale_main_offline
133 implicit none ;
private
135 #include <MOM_memory.h>
145 integer :: id_u = -1, id_v = -1, id_h = -1
147 integer :: id_ssh_inst = -1
153 real allocable_,
dimension(NIMEM_,NJMEM_,NKMEM_) :: &
154 h, & !< layer thickness [H ~> m or kg m-2]
155 t, & !< potential temperature [degC]
157 real allocable_,
dimension(NIMEMB_PTR_,NJMEM_,NKMEM_) :: &
158 u, & !< zonal velocity component [m s-1]
159 uh, & !< uh = u * h * dy at u grid points [H m2 s-1 ~> m3 s-1 or kg s-1]
161 real allocable_,
dimension(NIMEM_,NJMEMB_PTR_,NKMEM_) :: &
162 v, & !< meridional velocity [m s-1]
163 vh, & !< vh = v * h * dx at v grid points [H m2 s-1 ~> m3 s-1 or kg s-1]
165 real allocable_,
dimension(NIMEM_,NJMEM_) :: ssh_rint
167 real allocable_,
dimension(NIMEM_,NJMEM_) :: ave_ssh_ibc
170 real allocable_,
dimension(NIMEM_,NJMEM_) :: eta_av_bc
173 real,
dimension(:,:),
pointer :: &
175 real :: time_in_cycle
178 real :: time_in_thermo_cycle
187 real :: t_dyn_rel_adv
190 real :: t_dyn_rel_thermo
194 real :: t_dyn_rel_diag
196 integer :: ndyn_per_adv = 0
206 logical :: diabatic_first
208 logical :: use_ale_algorithm
211 logical :: offline_tracer_mode = .false.
215 type(time_type),
pointer :: time
218 logical :: thermo_spans_coupling
220 integer :: nstep_tot = 0
222 logical :: count_calls = .false.
228 logical :: do_dynamics
233 logical :: thickness_diffuse
234 logical :: thickness_diffuse_first
235 logical :: mixedlayer_restrat
238 real :: dtbt_reset_period
241 type(time_type) :: dtbt_reset_interval
242 type(time_type) :: dtbt_reset_time
245 real,
dimension(:,:,:),
pointer :: &
246 h_pre_dyn => null(), &
247 t_pre_dyn => null(), &
253 real,
dimension(:,:,:),
pointer :: &
257 logical :: interp_p_surf
260 logical :: p_surf_prev_set
263 real,
dimension(:,:),
pointer :: &
264 p_surf_prev => null(), &
265 p_surf_begin => null(), &
270 character(len=120) :: ic_file
273 logical :: calc_rho_for_sea_lev
288 logical :: check_bad_sfc_vals
289 real :: bad_val_ssh_max
290 real :: bad_val_sst_max
291 real :: bad_val_sst_min
292 real :: bad_val_sss_max
293 real :: bad_val_col_thick
341 type(
ale_cs),
pointer :: ale_csp => null()
352 logical :: ensemble_ocean
360 public initialize_mom, finish_mom_initialization, mom_end
361 public step_mom, step_offline
362 public extract_surface_state, get_ocean_stocks
363 public get_mom_state_elements, mom_state_is_synchronized
364 public allocate_surface_state, deallocate_surface_state
367 integer :: id_clock_ocean
368 integer :: id_clock_dynamics
369 integer :: id_clock_thermo
370 integer :: id_clock_tracer
371 integer :: id_clock_diabatic
372 integer :: id_clock_continuity
373 integer :: id_clock_thick_diff
374 integer :: id_clock_bbl_visc
375 integer :: id_clock_ml_restrat
376 integer :: id_clock_diagnostics
377 integer :: id_clock_z_diag
378 integer :: id_clock_init
379 integer :: id_clock_mom_init
380 integer :: id_clock_pass
381 integer :: id_clock_pass_init
382 integer :: id_clock_ale
383 integer :: id_clock_other
384 integer :: id_clock_offline_tracer
394 subroutine step_mom(forces, fluxes, sfc_state, Time_start, time_interval, CS, &
395 Waves, do_dynamics, do_thermodynamics, start_cycle, &
396 end_cycle, cycle_length, reset_therm)
398 type(
forcing),
intent(inout) :: fluxes
400 type(
surface),
intent(inout) :: sfc_state
401 type(time_type),
intent(in) :: time_start
402 real,
intent(in) :: time_interval
405 optional,
pointer :: waves
406 logical,
optional,
intent(in) :: do_dynamics
408 logical,
optional,
intent(in) :: do_thermodynamics
410 logical,
optional,
intent(in) :: start_cycle
413 logical,
optional,
intent(in) :: end_cycle
416 real,
optional,
intent(in) :: cycle_length
418 logical,
optional,
intent(in) :: reset_therm
431 integer :: i, j, k, is, ie, js, je, isq, ieq, jsq, jeq, nz, n
432 integer :: isd, ied, jsd, jed, isdb, iedb, jsdb, jedb
438 real :: dt_therm_here
440 real :: wt_end, wt_beg
444 real :: rel_time = 0.0
448 logical :: do_advection
449 logical :: do_calc_bbl
450 logical :: thermo_does_span_coupling
454 logical :: cycle_start
458 logical :: therm_reset
460 real,
dimension(SZI_(CS%G),SZJ_(CS%G)) :: &
463 real,
dimension(:,:,:),
pointer :: &
467 real,
dimension(:,:),
pointer :: &
471 type(time_type) :: time_local, end_time_thermo, time_temp
472 type(group_pass_type) :: pass_tau_ustar_psurf
473 logical :: showcalltree
475 g => cs%G ; gv => cs%GV ; us => cs%US
476 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
477 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
478 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
479 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
480 u => cs%u ; v => cs%v ; h => cs%h
482 do_dyn = .true. ;
if (
present(do_dynamics)) do_dyn = do_dynamics
483 do_thermo = .true. ;
if (
present(do_thermodynamics)) do_thermo = do_thermodynamics
484 if (.not.(do_dyn .or. do_thermo))
call mom_error(fatal,
"Step_MOM: "//&
485 "Both do_dynamics and do_thermodynamics are false, which makes no sense.")
486 cycle_start = .true. ;
if (
present(start_cycle)) cycle_start = start_cycle
487 cycle_end = .true. ;
if (
present(end_cycle)) cycle_end = end_cycle
488 cycle_time = time_interval ;
if (
present(cycle_length)) cycle_time = cycle_length
489 therm_reset = cycle_start ;
if (
present(reset_therm)) therm_reset = reset_therm
491 call cpu_clock_begin(id_clock_ocean)
492 call cpu_clock_begin(id_clock_other)
495 call mom_state_chksum(
"Beginning of step_MOM ", u, v, h, cs%uh, cs%vh, g, gv)
498 showcalltree = calltree_showquery()
499 if (showcalltree)
call calltree_enter(
"step_MOM(), MOM.F90")
505 if (time_interval > cs%dt) n_max = ceiling(time_interval/cs%dt - 0.001)
506 dt = time_interval / real(n_max)
507 thermo_does_span_coupling = (cs%thermo_spans_coupling .and. &
508 (cs%dt_therm > 1.5*cycle_time))
509 if (thermo_does_span_coupling)
then
511 dt_therm = cycle_time * floor(cs%dt_therm / cycle_time + 0.001)
512 ntstep = floor(dt_therm/dt + 0.001)
513 elseif (.not.do_thermo)
then
514 dt_therm = cs%dt_therm
515 if (
present(cycle_length)) dt_therm = min(cs%dt_therm, cycle_length)
518 ntstep = max(1,min(n_max,floor(cs%dt_therm/dt + 0.001)))
522 if (
associated(forces%p_surf)) p_surf => forces%p_surf
523 if (.not.
associated(forces%p_surf)) cs%interp_p_surf = .false.
526 call cpu_clock_begin(id_clock_pass)
528 if (
associated(forces%ustar)) &
530 if (
associated(forces%p_surf)) &
532 if (g%nonblocking_updates)
then
533 call start_group_pass(pass_tau_ustar_psurf, g%Domain)
535 call do_group_pass(pass_tau_ustar_psurf, g%Domain)
537 call cpu_clock_end(id_clock_pass)
541 if ((time_interval > cs%dt_therm) .and. (cs%dt_therm > 0.0)) &
542 n_max = ceiling(time_interval/cs%dt_therm - 0.001)
544 dt = time_interval / real(n_max)
545 dt_therm = dt ; ntstep = 1
546 if (
associated(fluxes%p_surf)) p_surf => fluxes%p_surf
548 if (cs%UseWaves)
call pass_var(fluxes%ustar, g%Domain, clock=id_clock_pass)
551 if (therm_reset)
then
552 cs%time_in_thermo_cycle = 0.0
553 if (
associated(cs%tv%frazil)) cs%tv%frazil(:,:) = 0.0
554 if (
associated(cs%tv%salt_deficit)) cs%tv%salt_deficit(:,:) = 0.0
555 if (
associated(cs%tv%TempxPmE)) cs%tv%TempxPmE(:,:) = 0.0
556 if (
associated(cs%tv%internal_heat)) cs%tv%internal_heat(:,:) = 0.0
559 if (cycle_start)
then
560 cs%time_in_cycle = 0.0
561 do j=js,je ;
do i=is,ie ; cs%ssh_rint(i,j) = 0.0 ;
enddo ;
enddo
563 if (
associated(cs%VarMix))
then
564 call enable_averaging(cycle_time, time_start + real_to_time(cycle_time), &
566 call calc_resoln_function(h, cs%tv, g, gv, us, cs%VarMix)
567 call disable_averaging(cs%diag)
572 if (g%nonblocking_updates) &
573 call complete_group_pass(pass_tau_ustar_psurf, g%Domain, clock=id_clock_pass)
575 if (cs%interp_p_surf)
then
576 if (.not.
associated(cs%p_surf_end))
allocate(cs%p_surf_end(isd:ied,jsd:jed))
577 if (.not.
associated(cs%p_surf_begin))
allocate(cs%p_surf_begin(isd:ied,jsd:jed))
578 if (.not.cs%p_surf_prev_set)
then
579 do j=jsd,jed ;
do i=isd,ied
580 cs%p_surf_prev(i,j) = forces%p_surf(i,j)
582 cs%p_surf_prev_set = .true.
585 cs%p_surf_end => forces%p_surf
588 if (cs%UseWaves)
then
590 call enable_averaging(time_interval, time_start + real_to_time(time_interval), cs%diag)
591 call update_stokes_drift(g, gv, us, waves, h, forces%ustar)
592 call disable_averaging(cs%diag)
596 call update_stokes_drift(g, gv, us, waves, h, fluxes%ustar)
603 if (do_dyn)
call mom_mech_forcing_chksum(
"Before steps", forces, g, us, haloshift=0)
604 if (do_dyn)
call check_redundant(
"Before steps ", forces%taux, forces%tauy, g)
606 call cpu_clock_end(id_clock_other)
610 rel_time = rel_time + dt
612 cs%Time = time_start + real_to_time(rel_time - 0.5*dt)
614 time_local = time_start + real_to_time(rel_time)
616 if (showcalltree)
call calltree_enter(
"DT cycles (step_MOM) n=",n)
620 if (cs%diabatic_first .and. (cs%t_dyn_rel_adv==0.0) .and. do_thermo)
then
622 if (.not.do_dyn)
then
624 elseif (thermo_does_span_coupling)
then
626 if ((fluxes%dt_buoy_accum > 0.0) .and. (dtdia > time_interval) .and. &
627 (abs(fluxes%dt_buoy_accum - dtdia) > 1e-6*dtdia))
then
628 call mom_error(fatal,
"step_MOM: Mismatch between long thermodynamic "//&
629 "timestep and time over which buoyancy fluxes have been accumulated.")
631 call mom_error(fatal,
"MOM is not yet set up to have restarts that work "//&
632 "with THERMO_SPANS_COUPLING and DIABATIC_FIRST.")
634 dtdia = dt*min(ntstep,n_max-(n-1))
637 end_time_thermo = time_local
641 cs%Time = cs%Time + real_to_time(0.5*(dtdia-dt))
644 end_time_thermo = time_local + real_to_time(dtdia-dt)
648 call step_mom_thermo(cs, g, gv, us, u, v, h, cs%tv, fluxes, dtdia, &
649 end_time_thermo, .true., waves=waves)
650 cs%time_in_thermo_cycle = cs%time_in_thermo_cycle + dtdia
653 cs%t_dyn_rel_thermo = -dtdia
654 if (showcalltree)
call calltree_waypoint(
"finished diabatic_first (step_MOM)")
657 cs%Time = time_start + real_to_time(rel_time - 0.5*dt)
664 if (cs%ndyn_per_adv == 0 .and. cs%t_dyn_rel_adv == 0.)
then
665 call diag_copy_diag_to_storage(cs%diag_pre_dyn, h, cs%diag)
666 cs%ndyn_per_adv = cs%ndyn_per_adv + 1
670 if (
associated(cs%u_prev) .and.
associated(cs%v_prev))
then
671 do k=1,nz ;
do j=jsd,jed ;
do i=isdb,iedb
672 cs%u_prev(i,j,k) = u(i,j,k)
673 enddo ;
enddo ;
enddo
674 do k=1,nz ;
do j=jsdb,jedb ;
do i=isd,ied
675 cs%v_prev(i,j,k) = v(i,j,k)
676 enddo ;
enddo ;
enddo
679 dt_therm_here = dt_therm
680 if (do_thermo .and. do_dyn .and. .not.thermo_does_span_coupling) &
681 dt_therm_here = dt*min(ntstep, n_max-n+1)
687 if ((cs%t_dyn_rel_adv == 0.0) .or. (n==1)) &
688 bbl_time_int = max(dt, min(dt_therm - cs%t_dyn_rel_adv, dt*(1+n_max-n)) )
690 if ((cs%t_dyn_rel_adv == 0.0) .or. ((n==1) .and. cycle_start)) &
691 bbl_time_int = min(dt_therm, cycle_time)
694 if (cs%interp_p_surf)
then
695 wt_end = real(n) / real(n_max)
696 wt_beg = real(n-1) / real(n_max)
697 do j=jsd,jed ;
do i=isd,ied
698 cs%p_surf_end(i,j) = wt_end * forces%p_surf(i,j) + &
699 (1.0-wt_end) * cs%p_surf_prev(i,j)
700 cs%p_surf_begin(i,j) = wt_beg * forces%p_surf(i,j) + &
701 (1.0-wt_beg) * cs%p_surf_prev(i,j)
705 call step_mom_dynamics(forces, cs%p_surf_begin, cs%p_surf_end, dt, &
706 dt_therm_here, bbl_time_int, cs, &
707 time_local, waves=waves)
712 if (thermo_does_span_coupling .or. .not.do_thermo)
then
713 do_advection = (cs%t_dyn_rel_adv + 0.5*dt > dt_therm)
715 do_advection = ((mod(n,ntstep) == 0) .or. (n==n_max))
718 if (do_advection)
then
719 call step_mom_tracer_dyn(cs, g, gv, h, time_local)
721 if (cs%diabatic_first .and. abs(cs%t_dyn_rel_thermo) > 1e-6*dt)
call mom_error(fatal, &
722 "step_MOM: Mismatch between the dynamics and diabatic times "//&
723 "with DIABATIC_FIRST.")
729 if ((cs%t_dyn_rel_adv==0.0) .and. do_thermo .and. (.not.cs%diabatic_first))
then
731 dtdia = cs%t_dyn_rel_thermo
734 if ((cs%t_dyn_rel_thermo==0.0) .and. .not.do_dyn) dtdia = dt
736 if (cs%thermo_spans_coupling .and. (cs%dt_therm > 1.5*cycle_time) .and. &
737 (abs(dt_therm - dtdia) > 1e-6*dt_therm))
then
738 call mom_error(fatal,
"step_MOM: Mismatch between dt_therm and dtdia "//&
739 "before call to diabatic.")
744 if (dtdia > dt) cs%Time = cs%Time - real_to_time(0.5*(dtdia-dt))
747 call step_mom_thermo(cs, g, gv, us, u, v, h, cs%tv, fluxes, dtdia, &
748 time_local, .false., waves=waves)
749 cs%time_in_thermo_cycle = cs%time_in_thermo_cycle + dtdia
751 if ((cs%t_dyn_rel_thermo==0.0) .and. .not.do_dyn)
then
753 cs%t_dyn_rel_thermo = -dtdia
755 cs%t_dyn_rel_thermo = 0.0
759 cs%Time = time_start + real_to_time(rel_time - 0.5*dt)
763 call cpu_clock_begin(id_clock_dynamics)
766 cs%time_in_cycle = cs%time_in_cycle + dt
767 call find_eta(h, cs%tv, g, gv, us, ssh, cs%eta_av_bc, eta_to_m=1.0)
768 do j=js,je ;
do i=is,ie
769 cs%ssh_rint(i,j) = cs%ssh_rint(i,j) + dt*ssh(i,j)
771 if (cs%IDs%id_ssh_inst > 0)
call post_data(cs%IDs%id_ssh_inst, ssh, cs%diag)
772 call cpu_clock_end(id_clock_dynamics)
777 if (mom_state_is_synchronized(cs))
then
779 call cpu_clock_begin(id_clock_other) ;
call cpu_clock_begin(id_clock_diagnostics)
782 call enable_averaging(cs%t_dyn_rel_diag, time_local, cs%diag)
783 call calculate_diagnostic_fields(u, v, h, cs%uh, cs%vh, cs%tv, cs%ADp, &
784 cs%CDp, p_surf, cs%t_dyn_rel_diag, cs%diag_pre_sync,&
785 g, gv, us, cs%diagnostics_CSp)
786 call post_tracer_diagnostics(cs%Tracer_reg, h, cs%diag_pre_sync, cs%diag, g, gv, cs%t_dyn_rel_diag)
787 call diag_copy_diag_to_storage(cs%diag_pre_sync, h, cs%diag)
788 if (showcalltree)
call calltree_waypoint(
"finished calculate_diagnostic_fields (step_MOM)")
789 call disable_averaging(cs%diag)
790 cs%t_dyn_rel_diag = 0.0
792 call cpu_clock_end(id_clock_diagnostics) ;
call cpu_clock_end(id_clock_other)
795 if (do_dyn .and. .not.cs%count_calls) cs%nstep_tot = cs%nstep_tot + 1
796 if (showcalltree)
call calltree_leave(
"DT cycles (step_MOM)")
800 if (cs%count_calls .and. cycle_start) cs%nstep_tot = cs%nstep_tot + 1
802 call cpu_clock_begin(id_clock_other)
804 if (cs%time_in_cycle > 0.0)
then
805 i_wt_ssh = 1.0/cs%time_in_cycle
806 do j=js,je ;
do i=is,ie
807 ssh(i,j) = cs%ssh_rint(i,j)*i_wt_ssh
808 cs%ave_ssh_ibc(i,j) = ssh(i,j)
811 call adjust_ssh_for_p_atm(cs%tv, g, gv, us, cs%ave_ssh_ibc, forces%p_surf_SSH, &
812 cs%calc_rho_for_sea_lev)
813 elseif (do_thermo)
then
814 call adjust_ssh_for_p_atm(cs%tv, g, gv, us, cs%ave_ssh_ibc, fluxes%p_surf_SSH, &
815 cs%calc_rho_for_sea_lev)
819 if (do_dyn .and. cs%interp_p_surf)
then ;
do j=jsd,jed ;
do i=isd,ied
820 cs%p_surf_prev(i,j) = forces%p_surf(i,j)
821 enddo ;
enddo ;
endif
823 if (cs%ensemble_ocean)
then
825 call set_analysis_time(cs%Time,cs%odaCS)
827 call set_prior_tracer(cs%Time, g, gv, cs%h, cs%tv, cs%odaCS)
829 call oda(cs%Time,cs%odaCS)
832 if (showcalltree)
call calltree_waypoint(
"calling extract_surface_state (step_MOM)")
833 call extract_surface_state(cs, sfc_state)
837 call cpu_clock_begin(id_clock_diagnostics)
838 if (cs%time_in_cycle > 0.0)
then
839 call enable_averaging(cs%time_in_cycle, time_local, cs%diag)
840 call post_surface_dyn_diags(cs%sfc_IDs, g, cs%diag, sfc_state, ssh)
842 if (cs%time_in_thermo_cycle > 0.0)
then
843 call enable_averaging(cs%time_in_thermo_cycle, time_local, cs%diag)
844 call post_surface_thermo_diags(cs%sfc_IDs, g, gv, us, cs%diag, cs%time_in_thermo_cycle, &
845 sfc_state, cs%tv, ssh, cs%ave_ssh_ibc)
847 call disable_averaging(cs%diag)
848 call cpu_clock_end(id_clock_diagnostics)
852 if (do_thermo .and. fluxes%fluxes_used) &
853 call accumulate_net_input(fluxes, sfc_state, fluxes%dt_buoy_accum, &
854 g, cs%sum_output_CSp)
856 if (mom_state_is_synchronized(cs)) &
857 call write_energy(cs%u, cs%v, cs%h, cs%tv, time_local, cs%nstep_tot, &
858 g, gv, us, cs%sum_output_CSp, cs%tracer_flow_CSp, &
859 dt_forcing=real_to_time(time_interval) )
861 call cpu_clock_end(id_clock_other)
863 if (showcalltree)
call calltree_leave(
"step_MOM()")
864 call cpu_clock_end(id_clock_ocean)
866 end subroutine step_mom
869 subroutine step_mom_dynamics(forces, p_surf_begin, p_surf_end, dt, dt_thermo, &
870 bbl_time_int, CS, Time_local, Waves)
872 real,
dimension(:,:),
pointer :: p_surf_begin
875 real,
dimension(:,:),
pointer :: p_surf_end
878 real,
intent(in) :: dt
879 real,
intent(in) :: dt_thermo
881 real,
intent(in) :: bbl_time_int
885 type(time_type),
intent(in) :: Time_local
887 optional,
pointer :: Waves
897 real,
dimension(:,:,:),
pointer :: &
904 logical :: showCallTree
906 integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz
907 integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB
909 g => cs%G ; gv => cs%GV ; us => cs%US ; ids => cs%IDs
910 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
911 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
912 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
913 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
914 u => cs%u ; v => cs%v ; h => cs%h
915 showcalltree = calltree_showquery()
917 call cpu_clock_begin(id_clock_dynamics)
919 if ((cs%t_dyn_rel_adv == 0.0) .and. cs%thickness_diffuse .and. cs%thickness_diffuse_first)
then
921 call enable_averaging(dt_thermo, time_local+real_to_time(dt_thermo-dt), cs%diag)
922 call cpu_clock_begin(id_clock_thick_diff)
923 if (
associated(cs%VarMix)) &
924 call calc_slope_functions(h, cs%tv, dt, g, gv, us, cs%VarMix)
925 call thickness_diffuse(h, cs%uhtr, cs%vhtr, cs%tv, dt_thermo, g, gv, us, &
926 cs%MEKE, cs%VarMix, cs%CDp, cs%thickness_diffuse_CSp)
927 call cpu_clock_end(id_clock_thick_diff)
928 call pass_var(h, g%Domain, clock=id_clock_pass)
929 call disable_averaging(cs%diag)
930 if (showcalltree)
call calltree_waypoint(
"finished thickness_diffuse_first (step_MOM)")
934 call diag_update_remap_grids(cs%diag)
938 if (bbl_time_int > 0.0)
then
939 call enable_averaging(bbl_time_int, &
940 time_local + real_to_time(bbl_time_int-dt), cs%diag)
942 call cpu_clock_begin(id_clock_bbl_visc)
943 call set_viscous_bbl(cs%u, cs%v, cs%h, cs%tv, cs%visc, g, gv, us, &
944 cs%set_visc_CSp, symmetrize=.true.)
945 call cpu_clock_end(id_clock_bbl_visc)
946 if (showcalltree)
call calltree_waypoint(
"done with set_viscous_BBL (step_MOM)")
947 call disable_averaging(cs%diag)
950 if (cs%do_dynamics .and. cs%split)
then
955 if (cs%dtbt_reset_period == 0.0) calc_dtbt = .true.
956 if (cs%dtbt_reset_period > 0.0)
then
957 if (time_local >= cs%dtbt_reset_time)
then
959 cs%dtbt_reset_time = cs%dtbt_reset_time + cs%dtbt_reset_interval
963 call step_mom_dyn_split_rk2(u, v, h, cs%tv, cs%visc, time_local, dt, forces, &
964 p_surf_begin, p_surf_end, cs%uh, cs%vh, cs%uhtr, cs%vhtr, &
965 cs%eta_av_bc, g, gv, us, cs%dyn_split_RK2_CSp, calc_dtbt, cs%VarMix, &
966 cs%MEKE, cs%thickness_diffuse_CSp, waves=waves)
967 if (showcalltree)
call calltree_waypoint(
"finished step_MOM_dyn_split (step_MOM)")
969 elseif (cs%do_dynamics)
then
978 call step_mom_dyn_unsplit_rk2(u, v, h, cs%tv, cs%visc, time_local, dt, forces, &
979 p_surf_begin, p_surf_end, cs%uh, cs%vh, cs%uhtr, cs%vhtr, &
980 cs%eta_av_bc, g, gv, us, cs%dyn_unsplit_RK2_CSp, cs%VarMix, cs%MEKE)
982 call step_mom_dyn_unsplit(u, v, h, cs%tv, cs%visc, time_local, dt, forces, &
983 p_surf_begin, p_surf_end, cs%uh, cs%vh, cs%uhtr, cs%vhtr, &
984 cs%eta_av_bc, g, gv, us, cs%dyn_unsplit_CSp, cs%VarMix, cs%MEKE, waves=waves)
986 if (showcalltree)
call calltree_waypoint(
"finished step_MOM_dyn_unsplit (step_MOM)")
990 if (cs%thickness_diffuse .and. .not.cs%thickness_diffuse_first)
then
991 call cpu_clock_begin(id_clock_thick_diff)
993 if (cs%debug)
call hchksum(h,
"Pre-thickness_diffuse h", g%HI, haloshift=0, scale=gv%H_to_m)
995 if (
associated(cs%VarMix)) &
996 call calc_slope_functions(h, cs%tv, dt, g, gv, us, cs%VarMix)
997 call thickness_diffuse(h, cs%uhtr, cs%vhtr, cs%tv, dt, g, gv, us, &
998 cs%MEKE, cs%VarMix, cs%CDp, cs%thickness_diffuse_CSp)
1000 if (cs%debug)
call hchksum(h,
"Post-thickness_diffuse h", g%HI, haloshift=1, scale=gv%H_to_m)
1001 call cpu_clock_end(id_clock_thick_diff)
1002 call pass_var(h, g%Domain, clock=id_clock_pass)
1003 if (showcalltree)
call calltree_waypoint(
"finished thickness_diffuse (step_MOM)")
1007 if (cs%mixedlayer_restrat)
then
1009 call hchksum(h,
"Pre-mixedlayer_restrat h", g%HI, haloshift=1, scale=gv%H_to_m)
1010 call uvchksum(
"Pre-mixedlayer_restrat uhtr", &
1011 cs%uhtr, cs%vhtr, g%HI, haloshift=0, scale=gv%H_to_m)
1013 call cpu_clock_begin(id_clock_ml_restrat)
1014 call mixedlayer_restrat(h, cs%uhtr, cs%vhtr, cs%tv, forces, dt, cs%visc%MLD, &
1015 cs%VarMix, g, gv, us, cs%mixedlayer_restrat_CSp)
1016 call cpu_clock_end(id_clock_ml_restrat)
1017 call pass_var(h, g%Domain, clock=id_clock_pass)
1019 call hchksum(h,
"Post-mixedlayer_restrat h", g%HI, haloshift=1, scale=gv%H_to_m)
1020 call uvchksum(
"Post-mixedlayer_restrat [uv]htr", &
1021 cs%uhtr, cs%vhtr, g%HI, haloshift=0, scale=gv%H_to_m)
1027 call diag_update_remap_grids(cs%diag)
1029 if (cs%useMEKE)
call step_forward_meke(cs%MEKE, h, cs%VarMix%SN_u, cs%VarMix%SN_v, &
1030 cs%visc, dt, g, gv, us, cs%MEKE_CSp, cs%uhtr, cs%vhtr)
1031 call disable_averaging(cs%diag)
1034 cs%t_dyn_rel_adv = cs%t_dyn_rel_adv + dt
1035 cs%t_dyn_rel_thermo = cs%t_dyn_rel_thermo + dt
1036 if (abs(cs%t_dyn_rel_thermo) < 1e-6*dt) cs%t_dyn_rel_thermo = 0.0
1037 cs%t_dyn_rel_diag = cs%t_dyn_rel_diag + dt
1039 call cpu_clock_end(id_clock_dynamics)
1041 call cpu_clock_begin(id_clock_other) ;
call cpu_clock_begin(id_clock_diagnostics)
1042 call enable_averaging(dt, time_local, cs%diag)
1044 if (ids%id_u > 0)
call post_data(ids%id_u, u, cs%diag)
1045 if (ids%id_v > 0)
call post_data(ids%id_v, v, cs%diag)
1046 if (ids%id_h > 0)
call post_data(ids%id_h, h, cs%diag)
1047 call disable_averaging(cs%diag)
1048 call cpu_clock_end(id_clock_diagnostics) ;
call cpu_clock_end(id_clock_other)
1050 end subroutine step_mom_dynamics
1055 subroutine step_mom_tracer_dyn(CS, G, GV, h, Time_local)
1059 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
1061 type(time_type),
intent(in) :: Time_local
1063 type(group_pass_type) :: pass_T_S
1064 logical :: showCallTree
1065 showcalltree = calltree_showquery()
1068 call cpu_clock_begin(id_clock_other)
1069 call hchksum(h,
"Pre-advection h", g%HI, haloshift=1, scale=gv%H_to_m)
1070 call uvchksum(
"Pre-advection uhtr", cs%uhtr, cs%vhtr, g%HI, &
1071 haloshift=0, scale=gv%H_to_m)
1072 if (
associated(cs%tv%T))
call hchksum(cs%tv%T,
"Pre-advection T", g%HI, haloshift=1)
1073 if (
associated(cs%tv%S))
call hchksum(cs%tv%S,
"Pre-advection S", g%HI, haloshift=1)
1074 if (
associated(cs%tv%frazil))
call hchksum(cs%tv%frazil, &
1075 "Pre-advection frazil", g%HI, haloshift=0)
1076 if (
associated(cs%tv%salt_deficit))
call hchksum(cs%tv%salt_deficit, &
1077 "Pre-advection salt deficit", g%HI, haloshift=0)
1079 call cpu_clock_end(id_clock_other)
1082 call cpu_clock_begin(id_clock_thermo) ;
call cpu_clock_begin(id_clock_tracer)
1083 call enable_averaging(cs%t_dyn_rel_adv, time_local, cs%diag)
1085 call advect_tracer(h, cs%uhtr, cs%vhtr, cs%OBC, cs%t_dyn_rel_adv, g, gv, &
1086 cs%tracer_adv_CSp, cs%tracer_Reg)
1087 call tracer_hordiff(h, cs%t_dyn_rel_adv, cs%MEKE, cs%VarMix, g, gv, &
1088 cs%tracer_diff_CSp, cs%tracer_Reg, cs%tv)
1089 if (showcalltree)
call calltree_waypoint(
"finished tracer advection/diffusion (step_MOM)")
1090 call cpu_clock_end(id_clock_tracer) ;
call cpu_clock_end(id_clock_thermo)
1092 call cpu_clock_begin(id_clock_other) ;
call cpu_clock_begin(id_clock_diagnostics)
1093 call post_transport_diagnostics(g, gv, cs%uhtr, cs%vhtr, h, cs%transport_IDs, &
1094 cs%diag_pre_dyn, cs%diag, cs%t_dyn_rel_adv, cs%tracer_reg)
1097 call diag_update_remap_grids(cs%diag)
1099 call disable_averaging(cs%diag)
1100 call cpu_clock_end(id_clock_diagnostics) ;
call cpu_clock_end(id_clock_other)
1104 call cpu_clock_begin(id_clock_thermo) ;
call cpu_clock_begin(id_clock_tracer)
1105 cs%uhtr(:,:,:) = 0.0
1106 cs%vhtr(:,:,:) = 0.0
1107 cs%t_dyn_rel_adv = 0.0
1108 call cpu_clock_end(id_clock_tracer) ;
call cpu_clock_end(id_clock_thermo)
1110 if (cs%diabatic_first .and.
associated(cs%tv%T))
then
1113 call create_group_pass(pass_t_s, cs%tv%T, g%Domain, to_all+omit_corners, halo=1)
1114 call create_group_pass(pass_t_s, cs%tv%S, g%Domain, to_all+omit_corners, halo=1)
1115 call do_group_pass(pass_t_s, g%Domain, clock=id_clock_pass)
1118 end subroutine step_mom_tracer_dyn
1122 subroutine step_mom_thermo(CS, G, GV, US, u, v, h, tv, fluxes, dtdia, &
1123 Time_end_thermo, update_BBL, Waves)
1128 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
1130 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)), &
1132 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
1135 type(
forcing),
intent(inout) :: fluxes
1136 real,
intent(in) :: dtdia
1137 type(time_type),
intent(in) :: Time_end_thermo
1138 logical,
intent(in) :: update_BBL
1140 optional,
pointer :: Waves
1143 logical :: use_ice_shelf
1144 logical :: showCallTree
1145 type(group_pass_type) :: pass_T_S, pass_T_S_h, pass_uv_T_S_h
1146 integer :: dynamics_stencil
1148 integer :: i, j, k, is, ie, js, je, nz
1150 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
1151 showcalltree = calltree_showquery()
1152 if (showcalltree)
call calltree_enter(
"step_MOM_thermo(), MOM.F90")
1154 use_ice_shelf = .false.
1155 if (
associated(fluxes%frac_shelf_h)) use_ice_shelf = .true.
1157 call enable_averaging(dtdia, time_end_thermo, cs%diag)
1159 if (
associated(cs%odaCS))
then
1160 call apply_oda_tracer_increments(dtdia,g,tv,h,cs%odaCS)
1163 if (update_bbl)
then
1168 call cpu_clock_begin(id_clock_bbl_visc)
1169 call set_viscous_bbl(u, v, h, tv, cs%visc, g, gv, us, cs%set_visc_CSp, symmetrize=.true.)
1170 call cpu_clock_end(id_clock_bbl_visc)
1171 if (showcalltree)
call calltree_waypoint(
"done with set_viscous_BBL (step_MOM_thermo)")
1174 call cpu_clock_begin(id_clock_thermo)
1175 if (.not.cs%adiabatic)
then
1177 call uvchksum(
"Pre-diabatic [uv]", u, v, g%HI, haloshift=2)
1178 call hchksum(h,
"Pre-diabatic h", g%HI, haloshift=1, scale=gv%H_to_m)
1179 call uvchksum(
"Pre-diabatic [uv]h", cs%uhtr, cs%vhtr, g%HI, &
1180 haloshift=0, scale=gv%H_to_m)
1182 call mom_thermo_chksum(
"Pre-diabatic ", tv, g,haloshift=0)
1184 call mom_forcing_chksum(
"Pre-diabatic", fluxes, g, us, haloshift=0)
1187 call cpu_clock_begin(id_clock_diabatic)
1188 call diabatic(u, v, h, tv, cs%Hml, fluxes, cs%visc, cs%ADp, cs%CDp, &
1189 dtdia, time_end_thermo, g, gv, us, cs%diabatic_CSp, waves=waves)
1190 fluxes%fluxes_used = .true.
1191 call cpu_clock_end(id_clock_diabatic)
1193 if (showcalltree)
call calltree_waypoint(
"finished diabatic (step_MOM_thermo)")
1198 if ( cs%use_ALE_algorithm )
then
1199 call enable_averaging(dtdia, time_end_thermo, cs%diag)
1201 if (
associated(tv%T)) &
1203 if (
associated(tv%S)) &
1206 call do_group_pass(pass_t_s_h, g%Domain)
1208 call preale_tracer_diagnostics(cs%tracer_Reg, g, gv)
1212 call hchksum(tv%T,
"Pre-ALE T", g%HI, haloshift=1)
1213 call hchksum(tv%S,
"Pre-ALE S", g%HI, haloshift=1)
1216 call cpu_clock_begin(id_clock_ale)
1217 if (use_ice_shelf)
then
1218 call ale_main(g, gv, us, h, u, v, tv, cs%tracer_Reg, cs%ALE_CSp, dtdia, &
1219 fluxes%frac_shelf_h)
1221 call ale_main(g, gv, us, h, u, v, tv, cs%tracer_Reg, cs%ALE_CSp, dtdia)
1224 if (showcalltree)
call calltree_waypoint(
"finished ALE_main (step_MOM_thermo)")
1225 call cpu_clock_end(id_clock_ale)
1228 dynamics_stencil = min(3, g%Domain%nihalo, g%Domain%njhalo)
1230 if (
associated(tv%T)) &
1232 if (
associated(tv%S)) &
1235 call do_group_pass(pass_uv_t_s_h, g%Domain, clock=id_clock_pass)
1237 if (cs%debug .and. cs%use_ALE_algorithm)
then
1239 call hchksum(tv%T,
"Post-ALE T", g%HI, haloshift=1)
1240 call hchksum(tv%S,
"Post-ALE S", g%HI, haloshift=1)
1247 call diag_update_remap_grids(cs%diag)
1250 call postale_tracer_diagnostics(cs%tracer_Reg, g, gv, cs%diag, dtdia)
1253 call uvchksum(
"Post-diabatic u", u, v, g%HI, haloshift=2)
1254 call hchksum(h,
"Post-diabatic h", g%HI, haloshift=1, scale=gv%H_to_m)
1255 call uvchksum(
"Post-diabatic [uv]h", cs%uhtr, cs%vhtr, g%HI, &
1256 haloshift=0, scale=gv%H_to_m)
1259 if (
associated(tv%T))
call hchksum(tv%T,
"Post-diabatic T", g%HI, haloshift=1)
1260 if (
associated(tv%S))
call hchksum(tv%S,
"Post-diabatic S", g%HI, haloshift=1)
1261 if (
associated(tv%frazil))
call hchksum(tv%frazil, &
1262 "Post-diabatic frazil", g%HI, haloshift=0)
1263 if (
associated(tv%salt_deficit))
call hchksum(tv%salt_deficit, &
1264 "Post-diabatic salt deficit", g%HI, haloshift=0)
1268 call disable_averaging(cs%diag)
1271 call cpu_clock_begin(id_clock_diabatic)
1272 call adiabatic(h, tv, fluxes, dtdia, g, gv, cs%diabatic_CSp)
1273 fluxes%fluxes_used = .true.
1274 call cpu_clock_end(id_clock_diabatic)
1276 if (
associated(tv%T))
then
1279 call do_group_pass(pass_t_s, g%Domain, clock=id_clock_pass)
1281 if (
associated(tv%T))
call hchksum(tv%T,
"Post-diabatic T", g%HI, haloshift=1)
1282 if (
associated(tv%S))
call hchksum(tv%S,
"Post-diabatic S", g%HI, haloshift=1)
1287 call cpu_clock_end(id_clock_thermo)
1289 call disable_averaging(cs%diag)
1291 if (showcalltree)
call calltree_leave(
"step_MOM_thermo(), MOM.F90")
1293 end subroutine step_mom_thermo
1300 subroutine step_offline(forces, fluxes, sfc_state, Time_start, time_interval, CS)
1302 type(
forcing),
intent(inout) :: fluxes
1303 type(
surface),
intent(inout) :: sfc_state
1304 type(time_type),
intent(in) :: time_start
1305 real,
intent(in) :: time_interval
1316 logical :: first_iter
1317 logical :: last_iter
1318 logical :: do_vertical
1319 logical :: adv_converged
1321 integer :: dt_offline, dt_offline_vertical
1322 logical :: skip_diffusion
1323 integer :: id_eta_diff_end
1325 integer,
pointer :: accumulated_time => null()
1327 integer :: is, ie, js, je, isd, ied, jsd, jed
1330 real,
dimension(:,:,:),
pointer :: &
1331 uhtr => null(), vhtr => null(), &
1332 eatr => null(), ebtr => null(), &
1336 real,
dimension(SZI_(CS%G),SZJ_(CS%G)) :: eta_pre, eta_end
1337 type(time_type) :: time_end
1340 g => cs%G ; gv => cs%GV ; us => cs%US
1342 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
1343 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
1345 call cpu_clock_begin(id_clock_offline_tracer)
1346 call extract_offline_main(cs%offline_CSp, uhtr, vhtr, eatr, ebtr, h_end, accumulated_time, &
1347 dt_offline, dt_offline_vertical, skip_diffusion)
1348 time_end = increment_date(time_start, seconds=floor(time_interval+0.001))
1349 call enable_averaging(time_interval, time_end, cs%diag)
1352 if (accumulated_time==0)
then
1355 first_iter = .false.
1359 if ( mod(accumulated_time, dt_offline_vertical) == 0 )
then
1360 do_vertical = .true.
1362 do_vertical = .false.
1366 accumulated_time = mod(accumulated_time + int(time_interval), dt_offline)
1367 if (accumulated_time==0)
then
1373 if (cs%use_ALE_algorithm)
then
1376 if (first_iter)
then
1377 call mom_mesg(
"Reading in new offline fields")
1382 call update_offline_fields(cs%offline_CSp, cs%h, fluxes, cs%use_ALE_algorithm)
1385 call offline_fw_fluxes_into_ocean(g, gv, cs%offline_CSp, fluxes, cs%h)
1387 if (.not.cs%diabatic_first)
then
1388 call offline_advection_ale(fluxes, time_start, time_interval, cs%offline_CSp, id_clock_ale, &
1389 cs%h, uhtr, vhtr, converged=adv_converged)
1392 call offline_redistribute_residual(cs%offline_CSp, cs%h, uhtr, vhtr, adv_converged)
1395 if (.not. skip_diffusion)
then
1396 if (
associated(cs%VarMix))
then
1398 call calc_resoln_function(cs%h, cs%tv, g, gv, us, cs%VarMix)
1399 call calc_slope_functions(cs%h, cs%tv, real(dt_offline), g, gv, us, cs%VarMix)
1401 call tracer_hordiff(cs%h, real(dt_offline), cs%MEKE, cs%VarMix, g, gv, &
1402 cs%tracer_diff_CSp, cs%tracer_Reg, cs%tv)
1407 if (do_vertical)
then
1408 call offline_diabatic_ale(fluxes, time_start, time_end, cs%offline_CSp, cs%h, eatr, ebtr)
1413 if (cs%diabatic_first)
then
1414 call offline_advection_ale(fluxes, time_start, time_interval, cs%offline_CSp, id_clock_ale, &
1415 cs%h, uhtr, vhtr, converged=adv_converged)
1418 call offline_redistribute_residual(cs%offline_CSp, cs%h, uhtr, vhtr, adv_converged)
1420 if (.not. skip_diffusion)
then
1421 if (
associated(cs%VarMix))
then
1423 call calc_resoln_function(cs%h, cs%tv, g, gv, us, cs%VarMix)
1424 call calc_slope_functions(cs%h, cs%tv, real(dt_offline), g, gv, us, cs%VarMix)
1426 call tracer_hordiff(cs%h, real(dt_offline), cs%MEKE, cs%VarMix, g, gv, &
1427 cs%tracer_diff_CSp, cs%tracer_Reg, cs%tv)
1431 call mom_mesg(
"Last iteration of offline interval")
1434 call offline_fw_fluxes_out_ocean(g, gv, cs%offline_CSp, fluxes, cs%h)
1437 call post_offline_convergence_diags(cs%offline_CSp, cs%h, h_end, uhtr, vhtr)
1441 call cpu_clock_begin(id_clock_ale)
1442 call ale_offline_tracer_final( g, gv, cs%h, cs%tv, h_end, cs%tracer_Reg, cs%ALE_CSp)
1443 call cpu_clock_end(id_clock_ale)
1447 call mom_error(warning, &
1448 "Offline tracer mode in non-ALE configuration has not been thoroughly tested")
1452 if (time_interval /= dt_offline)
then
1453 call mom_error(fatal, &
1454 "For offline tracer mode in a non-ALE configuration, dt_offline must equal time_interval")
1456 call update_offline_fields(cs%offline_CSp, cs%h, fluxes, cs%use_ALE_algorithm)
1457 call offline_advection_layer(fluxes, time_start, time_interval, cs%offline_CSp, &
1458 cs%h, eatr, ebtr, uhtr, vhtr)
1460 if (.not. skip_diffusion)
then
1461 call tracer_hordiff(h_end, real(dt_offline), cs%MEKE, cs%VarMix, g, gv, &
1462 cs%tracer_diff_CSp, cs%tracer_Reg, cs%tv)
1473 call adjust_ssh_for_p_atm(cs%tv, g, gv, us, cs%ave_ssh_ibc, forces%p_surf_SSH, &
1474 cs%calc_rho_for_sea_lev)
1475 call extract_surface_state(cs, sfc_state)
1477 call disable_averaging(cs%diag)
1482 fluxes%fluxes_used = .true.
1484 call cpu_clock_end(id_clock_offline_tracer)
1486 end subroutine step_offline
1490 subroutine initialize_mom(Time, Time_init, param_file, dirs, CS, restart_CSp, &
1491 Time_in, offline_tracer_mode, input_restart_file, diag_ptr, &
1492 count_calls, tracer_flow_CSp)
1493 type(time_type),
target,
intent(inout) :: time
1494 type(time_type),
intent(in) :: time_init
1501 type(time_type),
optional,
intent(in) :: time_in
1503 logical,
optional,
intent(out) :: offline_tracer_mode
1504 character(len=*),
optional,
intent(in) :: input_restart_file
1505 type(
diag_ctrl),
optional,
pointer :: diag_ptr
1508 optional,
pointer :: tracer_flow_csp
1510 logical,
optional,
intent(in) :: count_calls
1518 type(
diag_ctrl),
pointer :: diag => null()
1520 character(len=4),
parameter :: vers_num =
'v2.0'
1523 # include "version_variable.h"
1525 integer :: i, j, k, is, ie, js, je, isd, ied, jsd, jed, nz
1526 integer :: isdb, iedb, jsdb, jedb
1530 real,
allocatable,
dimension(:,:) :: eta
1531 real,
allocatable,
dimension(:,:) :: area_shelf_h
1532 real,
dimension(:,:),
allocatable,
target :: frac_shelf_h
1533 real,
dimension(:,:),
pointer :: shelf_area => null()
1535 type(group_pass_type) :: tmp_pass_uv_t_s_h, pass_uv_t_s_h
1538 logical :: write_geom_files
1539 logical :: ensemble_ocean
1541 logical :: use_geothermal
1543 logical :: symmetric
1545 logical :: do_unit_tests
1546 logical :: test_grid_copy = .false.
1548 logical :: bulkmixedlayer
1550 logical :: use_temperature
1551 logical :: use_frazil
1553 logical :: bound_salinity
1555 logical :: use_cont_abss
1559 logical :: advect_ts
1561 logical :: use_ice_shelf
1562 logical :: global_indexing
1564 logical :: bathy_at_vel
1566 logical :: calc_dtbt
1568 logical :: debug_truncations
1569 integer :: first_direction
1573 integer :: nkml, nkbl, verbosity, write_geom
1574 integer :: dynamics_stencil
1576 real :: conv2watt, conv2salt
1577 character(len=48) :: flux_units, s_flux_units
1580 type(time_type) :: start_time
1582 character(len=200) :: area_varname, ice_shelf_file, inputdir, filename
1584 if (
associated(cs))
then
1585 call mom_error(warning,
"initialize_MOM called with an associated "// &
1586 "control structure.")
1591 if (test_grid_copy)
then ;
allocate(g)
1592 else ; g => cs%G ;
endif
1596 id_clock_init = cpu_clock_id(
'Ocean Initialization', grain=clock_subcomponent)
1597 call cpu_clock_begin(id_clock_init)
1599 start_time = time ;
if (
present(time_in)) start_time = time_in
1603 call get_mom_input(param_file, dirs, default_input_filename=input_restart_file)
1605 verbosity = 2 ;
call read_param(param_file,
"VERBOSITY", verbosity)
1606 call mom_set_verbosity(verbosity)
1607 call calltree_enter(
"initialize_MOM(), MOM.F90")
1609 call find_obsolete_params(param_file)
1613 call get_param(param_file,
"MOM",
"VERBOSITY", verbosity, &
1614 "Integer controlling level of messaging\n" // &
1615 "\t0 = Only FATAL messages\n" // &
1616 "\t2 = Only FATAL, WARNING, NOTE [default]\n" // &
1617 "\t9 = All)", default=2, debuggingparam=.true.)
1618 call get_param(param_file,
"MOM",
"DO_UNIT_TESTS", do_unit_tests, &
1619 "If True, exercises unit tests at model start up.", &
1620 default=.false., debuggingparam=.true.)
1621 if (do_unit_tests)
then
1622 call unit_tests(verbosity)
1626 call unit_scaling_init(param_file, cs%US)
1630 call get_param(param_file,
"MOM",
"SPLIT", cs%split, &
1631 "Use the split time stepping if true.", default=.true.)
1633 cs%use_RK2 = .false.
1635 call get_param(param_file,
"MOM",
"USE_RK2", cs%use_RK2, &
1636 "If true, use RK2 instead of RK3 in the unsplit time stepping.", &
1640 call get_param(param_file,
"MOM",
"CALC_RHO_FOR_SEA_LEVEL", cs%calc_rho_for_sea_lev, &
1641 "If true, the in-situ density is used to calculate the "//&
1642 "effective sea level that is returned to the coupler. If false, "//&
1643 "the Boussinesq parameter RHO_0 is used.", default=.false.)
1644 call get_param(param_file,
"MOM",
"ENABLE_THERMODYNAMICS", use_temperature, &
1645 "If true, Temperature and salinity are used as state "//&
1646 "variables.", default=.true.)
1647 call get_param(param_file,
"MOM",
"USE_EOS", use_eos, &
1648 "If true, density is calculated from temperature and "//&
1649 "salinity with an equation of state. If USE_EOS is "//&
1650 "true, ENABLE_THERMODYNAMICS must be true as well.", &
1651 default=use_temperature)
1652 call get_param(param_file,
"MOM",
"DIABATIC_FIRST", cs%diabatic_first, &
1653 "If true, apply diabatic and thermodynamic processes, "//&
1654 "including buoyancy forcing and mass gain or loss, "//&
1655 "before stepping the dynamics forward.", default=.false.)
1656 call get_param(param_file,
"MOM",
"USE_CONTEMP_ABSSAL", use_cont_abss, &
1657 "If true, the prognostics T&S are the conservative temperature "//&
1658 "and absolute salinity. Care should be taken to convert them "//&
1659 "to potential temperature and practical salinity before "//&
1660 "exchanging them with the coupler and/or reporting T&S diagnostics.", &
1662 cs%tv%T_is_conT = use_cont_abss ; cs%tv%S_is_absS = use_cont_abss
1663 call get_param(param_file,
"MOM",
"ADIABATIC", cs%adiabatic, &
1664 "There are no diapycnal mass fluxes if ADIABATIC is "//&
1665 "true. This assumes that KD = KDML = 0.0 and that "//&
1666 "there is no buoyancy forcing, but makes the model "//&
1667 "faster by eliminating subroutine calls.", default=.false.)
1668 call get_param(param_file,
"MOM",
"DO_DYNAMICS", cs%do_dynamics, &
1669 "If False, skips the dynamics calls that update u & v, as well as "//&
1670 "the gravity wave adjustment to h. This is a fragile feature and "//&
1671 "thus undocumented.", default=.true., do_not_log=.true. )
1672 call get_param(param_file,
"MOM",
"ADVECT_TS", advect_ts, &
1673 "If True, advect temperature and salinity horizontally "//&
1674 "If False, T/S are registered for advection. "//&
1675 "This is intended only to be used in offline tracer mode "//&
1676 "and is by default false in that case.", &
1677 do_not_log = .true., default=.true. )
1678 if (
present(offline_tracer_mode))
then
1679 call get_param(param_file,
"MOM",
"OFFLINE_TRACER_MODE", cs%offline_tracer_mode, &
1680 "If true, barotropic and baroclinic dynamics, thermodynamics "//&
1681 "are all bypassed with all the fields necessary to integrate "//&
1682 "the tracer advection and diffusion equation are read in from "//&
1683 "files stored from a previous integration of the prognostic model. "//&
1684 "NOTE: This option only used in the ocean_solo_driver.", default=.false.)
1685 if (cs%offline_tracer_mode)
then
1686 call get_param(param_file,
"MOM",
"ADVECT_TS", advect_ts, &
1687 "If True, advect temperature and salinity horizontally "//&
1688 "If False, T/S are registered for advection. "//&
1689 "This is intended only to be used in offline tracer mode."//&
1690 "and is by default false in that case", &
1694 call get_param(param_file,
"MOM",
"USE_REGRIDDING", cs%use_ALE_algorithm, &
1695 "If True, use the ALE algorithm (regridding/remapping). "//&
1696 "If False, use the layered isopycnal algorithm.", default=.false. )
1697 call get_param(param_file,
"MOM",
"BULKMIXEDLAYER", bulkmixedlayer, &
1698 "If true, use a Kraus-Turner-like bulk mixed layer "//&
1699 "with transitional buffer layers. Layers 1 through "//&
1700 "NKML+NKBL have variable densities. There must be at "//&
1701 "least NKML+NKBL+1 layers if BULKMIXEDLAYER is true. "//&
1702 "BULKMIXEDLAYER can not be used with USE_REGRIDDING. "//&
1703 "The default is influenced by ENABLE_THERMODYNAMICS.", &
1704 default=use_temperature .and. .not.cs%use_ALE_algorithm)
1705 call get_param(param_file,
"MOM",
"THICKNESSDIFFUSE", cs%thickness_diffuse, &
1706 "If true, interface heights are diffused with a "//&
1707 "coefficient of KHTH.", default=.false.)
1708 call get_param(param_file,
"MOM",
"THICKNESSDIFFUSE_FIRST", &
1709 cs%thickness_diffuse_first, &
1710 "If true, do thickness diffusion before dynamics. "//&
1711 "This is only used if THICKNESSDIFFUSE is true.", &
1713 if (.not.cs%thickness_diffuse) cs%thickness_diffuse_first = .false.
1714 call get_param(param_file,
"MOM",
"BATHYMETRY_AT_VEL", bathy_at_vel, &
1715 "If true, there are separate values for the basin depths "//&
1716 "at velocity points. Otherwise the effects of topography "//&
1717 "are entirely determined from thickness points.", &
1719 call get_param(param_file,
"MOM",
"USE_WAVES", cs%UseWaves, default=.false., &
1722 call get_param(param_file,
"MOM",
"DEBUG", cs%debug, &
1723 "If true, write out verbose debugging data.", &
1724 default=.false., debuggingparam=.true.)
1725 call get_param(param_file,
"MOM",
"DEBUG_TRUNCATIONS", debug_truncations, &
1726 "If true, calculate all diagnostics that are useful for "//&
1727 "debugging truncations.", default=.false., debuggingparam=.true.)
1729 call get_param(param_file,
"MOM",
"DT", cs%dt, &
1730 "The (baroclinic) dynamics time step. The time-step that "//&
1731 "is actually used will be an integer fraction of the "//&
1732 "forcing time-step (DT_FORCING in ocean-only mode or the "//&
1733 "coupling timestep in coupled mode.)", units=
"s", &
1734 fail_if_missing=.true.)
1735 call get_param(param_file,
"MOM",
"DT_THERM", cs%dt_therm, &
1736 "The thermodynamic and tracer advection time step. "//&
1737 "Ideally DT_THERM should be an integer multiple of DT "//&
1738 "and less than the forcing or coupling time-step, unless "//&
1739 "THERMO_SPANS_COUPLING is true, in which case DT_THERM "//&
1740 "can be an integer multiple of the coupling timestep. By "//&
1741 "default DT_THERM is set to DT.", units=
"s", default=cs%dt)
1742 call get_param(param_file,
"MOM",
"THERMO_SPANS_COUPLING", cs%thermo_spans_coupling, &
1743 "If true, the MOM will take thermodynamic and tracer "//&
1744 "timesteps that can be longer than the coupling timestep. "//&
1745 "The actual thermodynamic timestep that is used in this "//&
1746 "case is the largest integer multiple of the coupling "//&
1747 "timestep that is less than or equal to DT_THERM.", default=.false.)
1749 if (bulkmixedlayer)
then
1750 cs%Hmix = -1.0 ; cs%Hmix_UV = -1.0
1752 call get_param(param_file,
"MOM",
"HMIX_SFC_PROP", cs%Hmix, &
1753 "If BULKMIXEDLAYER is false, HMIX_SFC_PROP is the depth "//&
1754 "over which to average to find surface properties like "//&
1755 "SST and SSS or density (but not surface velocities).", &
1756 units=
"m", default=1.0, scale=us%m_to_Z)
1757 call get_param(param_file,
"MOM",
"HMIX_UV_SFC_PROP", cs%Hmix_UV, &
1758 "If BULKMIXEDLAYER is false, HMIX_UV_SFC_PROP is the depth "//&
1759 "over which to average to find surface flow properties, "//&
1760 "SSU, SSV. A non-positive value indicates no averaging.", &
1761 units=
"m", default=0.0, scale=us%m_to_Z)
1763 call get_param(param_file,
"MOM",
"HFREEZE", cs%HFrz, &
1764 "If HFREEZE > 0, melt potential will be computed. The actual depth "//&
1765 "over which melt potential is computed will be min(HFREEZE, OBLD), "//&
1766 "where OBLD is the boundary layer depth. If HFREEZE <= 0 (default), "//&
1767 "melt potential will not be computed.", units=
"m", default=-1.0)
1768 call get_param(param_file,
"MOM",
"INTERPOLATE_P_SURF", cs%interp_p_surf, &
1769 "If true, linearly interpolate the surface pressure "//&
1770 "over the coupling time step, using the specified value "//&
1771 "at the end of the step.", default=.false.)
1774 call get_param(param_file,
"MOM",
"DTBT", dtbt, default=-0.98)
1775 default_val = cs%dt_therm ;
if (dtbt > 0.0) default_val = -1.0
1776 cs%dtbt_reset_period = -1.0
1777 call get_param(param_file,
"MOM",
"DTBT_RESET_PERIOD", cs%dtbt_reset_period, &
1778 "The period between recalculations of DTBT (if DTBT <= 0). "//&
1779 "If DTBT_RESET_PERIOD is negative, DTBT is set based "//&
1780 "only on information available at initialization. If 0, "//&
1781 "DTBT will be set every dynamics time step. The default "//&
1782 "is set by DT_THERM. This is only used if SPLIT is true.", &
1783 units=
"s", default=default_val, do_not_read=(dtbt > 0.0))
1787 use_frazil = .false. ; bound_salinity = .false. ; cs%tv%P_Ref = 2.0e7
1788 if (use_temperature)
then
1789 call get_param(param_file,
"MOM",
"FRAZIL", use_frazil, &
1790 "If true, water freezes if it gets too cold, and the "//&
1791 "the accumulated heat deficit is returned in the "//&
1792 "surface state. FRAZIL is only used if "//&
1793 "ENABLE_THERMODYNAMICS is true.", default=.false.)
1794 call get_param(param_file,
"MOM",
"DO_GEOTHERMAL", use_geothermal, &
1795 "If true, apply geothermal heating.", default=.false.)
1796 call get_param(param_file,
"MOM",
"BOUND_SALINITY", bound_salinity, &
1797 "If true, limit salinity to being positive. (The sea-ice "//&
1798 "model may ask for more salt than is available and "//&
1799 "drive the salinity negative otherwise.)", default=.false.)
1800 call get_param(param_file,
"MOM",
"MIN_SALINITY", cs%tv%min_salinity, &
1801 "The minimum value of salinity when BOUND_SALINITY=True. "//&
1802 "The default is 0.01 for backward compatibility but ideally "//&
1803 "should be 0.", units=
"PPT", default=0.01, do_not_log=.not.bound_salinity)
1804 call get_param(param_file,
"MOM",
"C_P", cs%tv%C_p, &
1805 "The heat capacity of sea water, approximated as a "//&
1806 "constant. This is only used if ENABLE_THERMODYNAMICS is "//&
1807 "true. The default value is from the TEOS-10 definition "//&
1808 "of conservative temperature.", units=
"J kg-1 K-1", &
1809 default=3991.86795711963)
1811 if (use_eos)
call get_param(param_file,
"MOM",
"P_REF", cs%tv%P_Ref, &
1812 "The pressure that is used for calculating the coordinate "//&
1813 "density. (1 Pa = 1e4 dbar, so 2e7 is commonly used.) "//&
1814 "This is only used if USE_EOS and ENABLE_THERMODYNAMICS "//&
1815 "are true.", units=
"Pa", default=2.0e7)
1817 if (bulkmixedlayer)
then
1818 call get_param(param_file,
"MOM",
"NKML", nkml, &
1819 "The number of sublayers within the mixed layer if "//&
1820 "BULKMIXEDLAYER is true.", units=
"nondim", default=2)
1821 call get_param(param_file,
"MOM",
"NKBL", nkbl, &
1822 "The number of layers that are used as variable density "//&
1823 "buffer layers if BULKMIXEDLAYER is true.", units=
"nondim", &
1827 call get_param(param_file,
"MOM",
"GLOBAL_INDEXING", global_indexing, &
1828 "If true, use a global lateral indexing convention, so "//&
1829 "that corresponding points on different processors have "//&
1830 "the same index. This does not work with static memory.", &
1831 default=.false., layoutparam=.true.)
1832 #ifdef STATIC_MEMORY_
1833 if (global_indexing)
call mom_error(fatal,
"initialize_MOM: "//&
1834 "GLOBAL_INDEXING can not be true with STATIC_MEMORY.")
1836 call get_param(param_file,
"MOM",
"FIRST_DIRECTION", first_direction, &
1837 "An integer that indicates which direction goes first "//&
1838 "in parts of the code that use directionally split "//&
1839 "updates, with even numbers (or 0) used for x- first "//&
1840 "and odd numbers used for y-first.", default=0)
1842 call get_param(param_file,
"MOM",
"CHECK_BAD_SURFACE_VALS", cs%check_bad_sfc_vals, &
1843 "If true, check the surface state for ridiculous values.", &
1845 if (cs%check_bad_sfc_vals)
then
1846 call get_param(param_file,
"MOM",
"BAD_VAL_SSH_MAX", cs%bad_val_ssh_max, &
1847 "The value of SSH above which a bad value message is "//&
1848 "triggered, if CHECK_BAD_SURFACE_VALS is true.", units=
"m", &
1850 call get_param(param_file,
"MOM",
"BAD_VAL_SSS_MAX", cs%bad_val_sss_max, &
1851 "The value of SSS above which a bad value message is "//&
1852 "triggered, if CHECK_BAD_SURFACE_VALS is true.", units=
"PPT", &
1854 call get_param(param_file,
"MOM",
"BAD_VAL_SST_MAX", cs%bad_val_sst_max, &
1855 "The value of SST above which a bad value message is "//&
1856 "triggered, if CHECK_BAD_SURFACE_VALS is true.", &
1857 units=
"deg C", default=45.0)
1858 call get_param(param_file,
"MOM",
"BAD_VAL_SST_MIN", cs%bad_val_sst_min, &
1859 "The value of SST below which a bad value message is "//&
1860 "triggered, if CHECK_BAD_SURFACE_VALS is true.", &
1861 units=
"deg C", default=-2.1)
1862 call get_param(param_file,
"MOM",
"BAD_VAL_COLUMN_THICKNESS", cs%bad_val_col_thick, &
1863 "The value of column thickness below which a bad value message is "//&
1864 "triggered, if CHECK_BAD_SURFACE_VALS is true.", units=
"m", &
1868 call get_param(param_file,
"MOM",
"SAVE_INITIAL_CONDS", save_ic, &
1869 "If true, write the initial conditions to a file given "//&
1870 "by IC_OUTPUT_FILE.", default=.false.)
1871 call get_param(param_file,
"MOM",
"IC_OUTPUT_FILE", cs%IC_file, &
1872 "The file into which to write the initial conditions.", &
1874 call get_param(param_file,
"MOM",
"WRITE_GEOM", write_geom, &
1875 "If =0, never write the geometry and vertical grid files. "//&
1876 "If =1, write the geometry and vertical grid files only for "//&
1877 "a new simulation. If =2, always write the geometry and "//&
1878 "vertical grid files. Other values are invalid.", default=1)
1879 if (write_geom<0 .or. write_geom>2)
call mom_error(fatal,
"MOM: "//&
1880 "WRITE_GEOM must be equal to 0, 1 or 2.")
1881 write_geom_files = ((write_geom==2) .or. ((write_geom==1) .and. &
1882 ((dirs%input_filename(1:1)==
'n') .and. (len_trim(dirs%input_filename)==1))))
1888 if (cs%use_ALE_algorithm .and. bulkmixedlayer)
call mom_error(fatal, &
1889 "MOM: BULKMIXEDLAYER can not currently be used with the ALE algorithm.")
1890 if (cs%use_ALE_algorithm .and. .not.use_temperature)
call mom_error(fatal, &
1891 "MOM: At this time, USE_EOS should be True when using the ALE algorithm.")
1892 if (cs%adiabatic .and. use_temperature)
call mom_error(warning, &
1893 "MOM: ADIABATIC and ENABLE_THERMODYNAMICS both defined is usually unwise.")
1894 if (use_eos .and. .not.use_temperature)
call mom_error(fatal, &
1895 "MOM: ENABLE_THERMODYNAMICS must be defined to use USE_EOS.")
1896 if (cs%adiabatic .and. bulkmixedlayer)
call mom_error(fatal, &
1897 "MOM: ADIABATIC and BULKMIXEDLAYER can not both be defined.")
1898 if (bulkmixedlayer .and. .not.use_eos)
call mom_error(fatal, &
1899 "initialize_MOM: A bulk mixed layer can only be used with T & S as "//&
1900 "state variables. Add USE_EOS = True to MOM_input.")
1902 call get_param(param_file,
'MOM',
"ICE_SHELF", use_ice_shelf, default=.false., do_not_log=.true.)
1903 if (use_ice_shelf)
then
1904 inputdir =
"." ;
call get_param(param_file,
'MOM',
"INPUTDIR", inputdir)
1905 inputdir = slasher(inputdir)
1906 call get_param(param_file,
'MOM',
"ICE_THICKNESS_FILE", ice_shelf_file, &
1907 "The file from which the ice bathymetry and area are read.", &
1908 fail_if_missing=.true.)
1909 call get_param(param_file,
'MOM',
"ICE_AREA_VARNAME", area_varname, &
1910 "The name of the area variable in ICE_THICKNESS_FILE.", &
1911 fail_if_missing=.true.)
1915 cs%ensemble_ocean=.false.
1916 call get_param(param_file,
"MOM",
"ENSEMBLE_OCEAN", cs%ensemble_ocean, &
1917 "If False, The model is being run in serial mode as a single realization. "//&
1918 "If True, The current model realization is part of a larger ensemble "//&
1919 "and at the end of step MOM, we will perform a gather of the ensemble "//&
1920 "members for statistical evaluation and/or data assimilation.", default=.false.)
1922 call calltree_waypoint(
"MOM parameters read (initialize_MOM)")
1925 #ifdef SYMMETRIC_MEMORY_
1930 #ifdef STATIC_MEMORY_
1931 call mom_domains_init(g%domain, param_file, symmetric=symmetric, &
1932 static_memory=.true., nihalo=nihalo_, njhalo=njhalo_, &
1933 niglobal=niglobal_, njglobal=njglobal_, niproc=niproc_, &
1936 call mom_domains_init(g%domain, param_file, symmetric=symmetric)
1938 call calltree_waypoint(
"domains initialized (initialize_MOM)")
1940 call mom_debugging_init(param_file)
1941 call diag_mediator_infrastructure_init()
1942 call mom_io_init(param_file)
1944 call hor_index_init(g%Domain, hi, param_file, &
1945 local_indexing=.not.global_indexing)
1947 call create_dyn_horgrid(dg, hi, bathymetry_at_vel=bathy_at_vel)
1950 call verticalgridinit( param_file, cs%GV, us )
1955 if (cs%debug .or. dg%symmetric) &
1958 call calltree_waypoint(
"grids initialized (initialize_MOM)")
1960 call mom_timing_init(cs)
1963 call mom_initialize_fixed(dg, us, cs%OBC, param_file, write_geom_files, dirs%output_directory)
1964 call calltree_waypoint(
"returned from MOM_initialize_fixed() (initialize_MOM)")
1966 if (
associated(cs%OBC))
call call_obc_register(param_file, cs%update_OBC_CSp, cs%OBC)
1968 call tracer_registry_init(param_file, cs%tracer_Reg)
1971 is = dg%isc ; ie = dg%iec ; js = dg%jsc ; je = dg%jec ; nz = gv%ke
1972 isd = dg%isd ; ied = dg%ied ; jsd = dg%jsd ; jed = dg%jed
1973 isdb = dg%IsdB ; iedb = dg%IedB ; jsdb = dg%JsdB ; jedb = dg%JedB
1974 alloc_(cs%u(isdb:iedb,jsd:jed,nz)) ; cs%u(:,:,:) = 0.0
1975 alloc_(cs%v(isd:ied,jsdb:jedb,nz)) ; cs%v(:,:,:) = 0.0
1976 alloc_(cs%h(isd:ied,jsd:jed,nz)) ; cs%h(:,:,:) = gv%Angstrom_H
1977 alloc_(cs%uh(isdb:iedb,jsd:jed,nz)) ; cs%uh(:,:,:) = 0.0
1978 alloc_(cs%vh(isd:ied,jsdb:jedb,nz)) ; cs%vh(:,:,:) = 0.0
1979 if (use_temperature)
then
1980 alloc_(cs%T(isd:ied,jsd:jed,nz)) ; cs%T(:,:,:) = 0.0
1981 alloc_(cs%S(isd:ied,jsd:jed,nz)) ; cs%S(:,:,:) = 0.0
1982 cs%tv%T => cs%T ; cs%tv%S => cs%S
1983 if (cs%tv%T_is_conT)
then
1984 vd_t = var_desc(name=
"contemp", units=
"Celsius", longname=
"Conservative Temperature", &
1985 cmor_field_name=
"thetao", cmor_longname=
"Sea Water Potential Temperature", &
1986 conversion=cs%tv%C_p)
1988 vd_t = var_desc(name=
"temp", units=
"degC", longname=
"Potential Temperature", &
1989 cmor_field_name=
"thetao", cmor_longname=
"Sea Water Potential Temperature", &
1990 conversion=cs%tv%C_p)
1992 if (cs%tv%S_is_absS)
then
1993 vd_s = var_desc(name=
"abssalt",units=
"g kg-1",longname=
"Absolute Salinity", &
1994 cmor_field_name=
"so", cmor_longname=
"Sea Water Salinity", &
1997 vd_s = var_desc(name=
"salt",units=
"psu",longname=
"Salinity", &
1998 cmor_field_name=
"so", cmor_longname=
"Sea Water Salinity", &
2003 s_flux_units = get_tr_flux_units(gv,
"psu")
2004 conv2watt = gv%H_to_kg_m2 * cs%tv%C_p
2005 if (gv%Boussinesq)
then
2006 conv2salt = gv%H_to_m
2008 conv2salt = gv%H_to_kg_m2
2010 call register_tracer(cs%tv%T, cs%tracer_Reg, param_file, dg%HI, gv, &
2011 tr_desc=vd_t, registry_diags=.true., flux_nameroot=
'T', &
2012 flux_units=
'W', flux_longname=
'Heat', &
2013 flux_scale=conv2watt, convergence_units=
'W m-2', &
2014 convergence_scale=conv2watt, cmor_tendprefix=
"opottemp", diag_form=2)
2015 call register_tracer(cs%tv%S, cs%tracer_Reg, param_file, dg%HI, gv, &
2016 tr_desc=vd_s, registry_diags=.true., flux_nameroot=
'S', &
2017 flux_units=s_flux_units, flux_longname=
'Salt', &
2018 flux_scale=conv2salt, convergence_units=
'kg m-2 s-1', &
2019 convergence_scale=0.001*gv%H_to_kg_m2, cmor_tendprefix=
"osalt", diag_form=2)
2021 if (
associated(cs%OBC)) &
2022 call register_temp_salt_segments(gv, cs%OBC, cs%tracer_Reg, param_file)
2024 if (use_frazil)
then
2025 allocate(cs%tv%frazil(isd:ied,jsd:jed)) ; cs%tv%frazil(:,:) = 0.0
2027 if (bound_salinity)
then
2028 allocate(cs%tv%salt_deficit(isd:ied,jsd:jed)) ; cs%tv%salt_deficit(:,:)=0.0
2031 if (bulkmixedlayer .or. use_temperature)
then
2032 allocate(cs%Hml(isd:ied,jsd:jed)) ; cs%Hml(:,:) = 0.0
2035 if (bulkmixedlayer)
then
2036 gv%nkml = nkml ; gv%nk_rho_varies = nkml + nkbl
2038 gv%nkml = 0 ; gv%nk_rho_varies = 0
2040 if (cs%use_ALE_algorithm)
then
2041 call get_param(param_file,
"MOM",
"NK_RHO_VARIES", gv%nk_rho_varies, default=0)
2044 alloc_(cs%uhtr(isdb:iedb,jsd:jed,nz)) ; cs%uhtr(:,:,:) = 0.0
2045 alloc_(cs%vhtr(isd:ied,jsdb:jedb,nz)) ; cs%vhtr(:,:,:) = 0.0
2046 cs%t_dyn_rel_adv = 0.0 ; cs%t_dyn_rel_thermo = 0.0 ; cs%t_dyn_rel_diag = 0.0
2048 if (debug_truncations)
then
2049 allocate(cs%u_prev(isdb:iedb,jsd:jed,nz)) ; cs%u_prev(:,:,:) = 0.0
2050 allocate(cs%v_prev(isd:ied,jsdb:jedb,nz)) ; cs%v_prev(:,:,:) = 0.0
2051 call safe_alloc_ptr(cs%ADp%du_dt_visc,isdb,iedb,jsd,jed,nz)
2052 call safe_alloc_ptr(cs%ADp%dv_dt_visc,isd,ied,jsdb,jedb,nz)
2053 if (.not.cs%adiabatic)
then
2054 call safe_alloc_ptr(cs%ADp%du_dt_dia,isdb,iedb,jsd,jed,nz)
2055 call safe_alloc_ptr(cs%ADp%dv_dt_dia,isd,ied,jsdb,jedb,nz)
2059 mom_internal_state%u => cs%u ; mom_internal_state%v => cs%v
2060 mom_internal_state%h => cs%h
2061 mom_internal_state%uh => cs%uh ; mom_internal_state%vh => cs%vh
2062 if (use_temperature)
then
2063 mom_internal_state%T => cs%T ; mom_internal_state%S => cs%S
2066 cs%CDp%uh => cs%uh ; cs%CDp%vh => cs%vh
2068 if (cs%interp_p_surf)
then
2069 allocate(cs%p_surf_prev(isd:ied,jsd:jed)) ; cs%p_surf_prev(:,:) = 0.0
2072 alloc_(cs%ssh_rint(isd:ied,jsd:jed)) ; cs%ssh_rint(:,:) = 0.0
2073 alloc_(cs%ave_ssh_ibc(isd:ied,jsd:jed)) ; cs%ave_ssh_ibc(:,:) = 0.0
2074 alloc_(cs%eta_av_bc(isd:ied,jsd:jed)) ; cs%eta_av_bc(:,:) = 0.0
2075 cs%time_in_cycle = 0.0 ; cs%time_in_thermo_cycle = 0.0
2080 if (use_eos)
call eos_init(param_file, cs%tv%eqn_of_state)
2081 if (use_temperature)
then
2082 allocate(cs%tv%TempxPmE(isd:ied,jsd:jed))
2083 cs%tv%TempxPmE(:,:) = 0.0
2084 if (use_geothermal)
then
2085 allocate(cs%tv%internal_heat(isd:ied,jsd:jed))
2086 cs%tv%internal_heat(:,:) = 0.0
2089 call calltree_waypoint(
"state variables allocated (initialize_MOM)")
2093 call restart_init(param_file, restart_csp)
2094 call set_restart_fields(gv, us, param_file, cs, restart_csp)
2096 call register_restarts_dyn_split_rk2(dg%HI, gv, param_file, &
2097 cs%dyn_split_RK2_CSp, restart_csp, cs%uh, cs%vh)
2098 elseif (cs%use_RK2)
then
2099 call register_restarts_dyn_unsplit_rk2(dg%HI, gv, param_file, &
2100 cs%dyn_unsplit_RK2_CSp, restart_csp)
2102 call register_restarts_dyn_unsplit(dg%HI, gv, param_file, &
2103 cs%dyn_unsplit_CSp, restart_csp)
2108 call call_tracer_register(dg%HI, gv, us, param_file, cs%tracer_flow_CSp, &
2109 cs%tracer_Reg, restart_csp)
2111 call meke_alloc_register_restart(dg%HI, param_file, cs%MEKE, restart_csp)
2112 call set_visc_register_restarts(dg%HI, gv, param_file, cs%visc, restart_csp)
2113 call mixedlayer_restrat_register_restarts(dg%HI, param_file, &
2114 cs%mixedlayer_restrat_CSp, restart_csp)
2116 if (
associated(cs%OBC)) &
2117 call open_boundary_register_restarts(dg%HI, gv, cs%OBC, restart_csp)
2119 call calltree_waypoint(
"restart registration complete (initialize_MOM)")
2122 call cpu_clock_begin(id_clock_mom_init)
2123 call mom_initialize_coord(gv, us, param_file, write_geom_files, &
2124 dirs%output_directory, cs%tv, dg%max_depth)
2125 call calltree_waypoint(
"returned from MOM_initialize_coord() (initialize_MOM)")
2127 if (cs%use_ALE_algorithm)
then
2128 call ale_init(param_file, gv, us, dg%max_depth, cs%ALE_CSp)
2129 call calltree_waypoint(
"returned from ALE_init() (initialize_MOM)")
2136 call mom_grid_init(g, param_file, hi, bathymetry_at_vel=bathy_at_vel)
2137 call copy_dyngrid_to_mom_grid(dg, g)
2138 call destroy_dyn_horgrid(dg)
2141 call set_first_direction(g, first_direction)
2143 if (cs%debug .or. g%symmetric)
then
2145 else ; g%Domain_aux => g%Domain ;
endif
2148 g%ke = gv%ke ; g%g_Earth = gv%mks_g_Earth
2150 call mom_initialize_state(cs%u, cs%v, cs%h, cs%tv, time, g, gv, us, param_file, &
2151 dirs, restart_csp, cs%ALE_CSp, cs%tracer_Reg, &
2152 cs%sponge_CSp, cs%ALE_sponge_CSp, cs%OBC, time_in)
2153 call cpu_clock_end(id_clock_mom_init)
2154 call calltree_waypoint(
"returned from MOM_initialize_state() (initialize_MOM)")
2159 if (test_grid_copy)
then
2161 call create_dyn_horgrid(dg, g%HI)
2165 call mom_grid_init(cs%G, param_file)
2167 call copy_mom_grid_to_dyngrid(g, dg)
2168 call copy_dyngrid_to_mom_grid(dg, cs%G)
2170 call destroy_dyn_horgrid(dg)
2171 call mom_grid_end(g) ;
deallocate(g)
2174 if (cs%debug .or. cs%G%symmetric)
then
2176 else ; cs%G%Domain_aux => cs%G%Domain ;
endif
2177 g%ke = gv%ke ; g%g_Earth = gv%mks_g_Earth
2185 if (ale_remap_init_conds(cs%ALE_CSp) .and. .not.
query_initialized(cs%h,
"h",restart_csp))
then
2190 call uvchksum(
"Pre ALE adjust init cond [uv]", &
2191 cs%u, cs%v, g%HI, haloshift=1)
2192 call hchksum(cs%h,
"Pre ALE adjust init cond h", g%HI, haloshift=1, scale=gv%H_to_m)
2194 call calltree_waypoint(
"Calling adjustGridForIntegrity() to remap initial conditions (initialize_MOM)")
2195 call adjustgridforintegrity(cs%ALE_CSp, g, gv, cs%h )
2196 call calltree_waypoint(
"Calling ALE_main() to remap initial conditions (initialize_MOM)")
2197 if (use_ice_shelf)
then
2198 filename = trim(inputdir)//trim(ice_shelf_file)
2199 if (.not.
file_exists(filename, g%Domain))
call mom_error(fatal, &
2200 "MOM: Unable to open "//trim(filename))
2202 allocate(area_shelf_h(isd:ied,jsd:jed))
2203 allocate(frac_shelf_h(isd:ied,jsd:jed))
2204 call mom_read_data(filename, trim(area_varname), area_shelf_h, g%Domain)
2206 frac_shelf_h(:,:) = 0.0
2208 do j=jsd,jed ;
do i=isd,ied
2209 if (g%areaT(i,j) > 0.0) &
2210 frac_shelf_h(i,j) = area_shelf_h(i,j) / g%areaT(i,j)
2213 shelf_area => frac_shelf_h
2214 call ale_main(g, gv, us, cs%h, cs%u, cs%v, cs%tv, cs%tracer_Reg, cs%ALE_CSp, &
2215 frac_shelf_h = shelf_area)
2217 call ale_main( g, gv, us, cs%h, cs%u, cs%v, cs%tv, cs%tracer_Reg, cs%ALE_CSp )
2220 call cpu_clock_begin(id_clock_pass_init)
2222 if (use_temperature)
then
2227 call do_group_pass(tmp_pass_uv_t_s_h, g%Domain)
2228 call cpu_clock_end(id_clock_pass_init)
2231 call uvchksum(
"Post ALE adjust init cond [uv]", cs%u, cs%v, g%HI, haloshift=1)
2232 call hchksum(cs%h,
"Post ALE adjust init cond h", g%HI, haloshift=1, scale=gv%H_to_m)
2235 if ( cs%use_ALE_algorithm )
call ale_updateverticalgridtype( cs%ALE_CSp, gv )
2239 call diag_mediator_init(g, gv, us, gv%ke, param_file, diag, doc_file_dir=dirs%output_directory)
2240 if (
present(diag_ptr)) diag_ptr => cs%diag
2245 call diag_masks_set(g, gv%ke, diag)
2249 call diag_set_state_ptrs(cs%h, cs%T, cs%S, cs%tv%eqn_of_state, diag)
2253 call set_axes_info(g, gv, us, param_file, diag)
2258 call diag_update_remap_grids(diag)
2261 call diag_grid_storage_init(cs%diag_pre_sync, g, diag)
2262 call diag_grid_storage_init(cs%diag_pre_dyn, g, diag)
2268 call set_masks_for_axes(g, diag)
2271 call write_static_fields(g, gv, us, cs%tv, cs%diag)
2272 call calltree_waypoint(
"static fields written (initialize_MOM)")
2275 call register_cell_measure(g, cs%diag, time)
2277 call cpu_clock_begin(id_clock_mom_init)
2278 if (cs%use_ALE_algorithm)
then
2279 call ale_writecoordinatefile( cs%ALE_CSp, gv, dirs%output_directory )
2281 call cpu_clock_end(id_clock_mom_init)
2282 call calltree_waypoint(
"ALE initialized (initialize_MOM)")
2284 cs%useMEKE = meke_init(time, g, us, param_file, diag, cs%MEKE_CSp, cs%MEKE, restart_csp)
2286 call varmix_init(time, g, gv, us, param_file, diag, cs%VarMix)
2287 call set_visc_init(time, g, gv, us, param_file, diag, cs%visc, cs%set_visc_CSp, restart_csp, cs%OBC)
2288 call thickness_diffuse_init(time, g, gv, us, param_file, diag, cs%CDp, cs%thickness_diffuse_CSp)
2290 allocate(eta(szi_(g),szj_(g))) ; eta(:,:) = 0.0
2291 call initialize_dyn_split_rk2(cs%u, cs%v, cs%h, cs%uh, cs%vh, eta, time, &
2292 g, gv, us, param_file, diag, cs%dyn_split_RK2_CSp, restart_csp, &
2293 cs%dt, cs%ADp, cs%CDp, mom_internal_state, cs%VarMix, cs%MEKE, &
2294 cs%thickness_diffuse_CSp, &
2295 cs%OBC, cs%update_OBC_CSp, cs%ALE_CSp, cs%set_visc_CSp, &
2296 cs%visc, dirs, cs%ntrunc, calc_dtbt=calc_dtbt)
2297 if (cs%dtbt_reset_period > 0.0)
then
2298 cs%dtbt_reset_interval = real_to_time(cs%dtbt_reset_period)
2300 cs%dtbt_reset_time = time_init + cs%dtbt_reset_interval * &
2301 ((time - time_init) / cs%dtbt_reset_interval)
2302 if ((cs%dtbt_reset_time > time) .and. calc_dtbt)
then
2306 cs%dtbt_reset_time = cs%dtbt_reset_time - cs%dtbt_reset_interval
2309 elseif (cs%use_RK2)
then
2310 call initialize_dyn_unsplit_rk2(cs%u, cs%v, cs%h, time, g, gv, us, &
2311 param_file, diag, cs%dyn_unsplit_RK2_CSp, restart_csp, &
2312 cs%ADp, cs%CDp, mom_internal_state, cs%MEKE, cs%OBC, &
2313 cs%update_OBC_CSp, cs%ALE_CSp, cs%set_visc_CSp, cs%visc, dirs, &
2316 call initialize_dyn_unsplit(cs%u, cs%v, cs%h, time, g, gv, us, &
2317 param_file, diag, cs%dyn_unsplit_CSp, restart_csp, &
2318 cs%ADp, cs%CDp, mom_internal_state, cs%MEKE, cs%OBC, &
2319 cs%update_OBC_CSp, cs%ALE_CSp, cs%set_visc_CSp, cs%visc, dirs, &
2322 call calltree_waypoint(
"dynamics initialized (initialize_MOM)")
2324 cs%mixedlayer_restrat = mixedlayer_restrat_init(time, g, gv, us, param_file, diag, &
2325 cs%mixedlayer_restrat_CSp, restart_csp)
2326 if (cs%mixedlayer_restrat)
then
2327 if (.not.(bulkmixedlayer .or. cs%use_ALE_algorithm)) &
2328 call mom_error(fatal,
"MOM: MIXEDLAYER_RESTRAT true requires a boundary layer scheme.")
2330 if (.not. cs%diabatic_first .and.
associated(cs%visc%MLD)) &
2331 call pass_var(cs%visc%MLD, g%domain, halo=1)
2334 call mom_diagnostics_init(mom_internal_state, cs%ADp, cs%CDp, time, g, gv, us, &
2335 param_file, diag, cs%diagnostics_CSp, cs%tv)
2336 call diag_copy_diag_to_storage(cs%diag_pre_sync, cs%h, cs%diag)
2338 if (
associated(cs%sponge_CSp)) &
2339 call init_sponge_diags(time, g, diag, cs%sponge_CSp)
2341 if (
associated(cs%ALE_sponge_CSp)) &
2342 call init_ale_sponge_diags(time, g, diag, cs%ALE_sponge_CSp)
2344 if (cs%adiabatic)
then
2345 call adiabatic_driver_init(time, g, param_file, diag, cs%diabatic_CSp, &
2348 call diabatic_driver_init(time, g, gv, us, param_file, cs%use_ALE_algorithm, diag, &
2349 cs%ADp, cs%CDp, cs%diabatic_CSp, cs%tracer_flow_CSp, &
2350 cs%sponge_CSp, cs%ALE_sponge_CSp)
2353 call tracer_advect_init(time, g, param_file, diag, cs%tracer_adv_CSp)
2354 call tracer_hor_diff_init(time, g, param_file, diag, cs%tv%eqn_of_state, &
2357 call lock_tracer_registry(cs%tracer_Reg)
2358 call calltree_waypoint(
"tracer registry now locked (initialize_MOM)")
2361 call register_surface_diags(time, g, cs%sfc_IDs, cs%diag, cs%tv)
2362 call register_diags(time, g, gv, cs%IDs, cs%diag)
2363 call register_transport_diags(time, g, gv, cs%transport_IDs, cs%diag)
2364 call register_tracer_diagnostics(cs%tracer_Reg, cs%h, time, diag, g, gv, &
2365 cs%use_ALE_algorithm)
2366 if (cs%use_ALE_algorithm)
then
2367 call ale_register_diags(time, g, gv, us, diag, cs%ALE_CSp)
2371 new_sim = is_new_run(restart_csp)
2372 call tracer_flow_control_init(.not.new_sim, time, g, gv, us, cs%h, param_file, &
2373 cs%diag, cs%OBC, cs%tracer_flow_CSp, cs%sponge_CSp, &
2374 cs%ALE_sponge_CSp, cs%tv)
2375 if (
present(tracer_flow_csp)) tracer_flow_csp => cs%tracer_flow_CSp
2379 if (
present(offline_tracer_mode)) offline_tracer_mode=cs%offline_tracer_mode
2381 if (cs%offline_tracer_mode)
then
2383 call offline_transport_init(param_file, cs%offline_CSp, cs%diabatic_CSp, g, gv)
2384 call insert_offline_main( cs=cs%offline_CSp, ale_csp=cs%ALE_CSp, diabatic_csp=cs%diabatic_CSp, &
2385 diag=cs%diag, obc=cs%OBC, tracer_adv_csp=cs%tracer_adv_CSp, &
2386 tracer_flow_csp=cs%tracer_flow_CSp, tracer_reg=cs%tracer_Reg, &
2387 tv=cs%tv, x_before_y = (mod(first_direction,2)==0), debug=cs%debug )
2388 call register_diags_offline_transport(time, cs%diag, cs%offline_CSp)
2392 call cpu_clock_begin(id_clock_pass_init)
2393 dynamics_stencil = min(3, g%Domain%nihalo, g%Domain%njhalo)
2394 call create_group_pass(pass_uv_t_s_h, cs%u, cs%v, g%Domain, halo=dynamics_stencil)
2395 if (use_temperature)
then
2401 call do_group_pass(pass_uv_t_s_h, g%Domain)
2403 if (
associated(cs%visc%Kv_shear)) &
2404 call pass_var(cs%visc%Kv_shear, g%Domain, to_all+omit_corners, halo=1)
2406 if (
associated(cs%visc%Kv_slow)) &
2407 call pass_var(cs%visc%Kv_slow, g%Domain, to_all+omit_corners, halo=1)
2409 call cpu_clock_end(id_clock_pass_init)
2411 call register_obsolete_diagnostics(param_file, cs%diag)
2413 if (use_frazil)
then
2415 cs%tv%frazil(:,:) = 0.0
2418 if (cs%interp_p_surf)
then
2419 cs%p_surf_prev_set = &
2422 if (cs%p_surf_prev_set)
call pass_var(cs%p_surf_prev, g%domain)
2427 call find_eta(cs%h, cs%tv, g, gv, us, cs%ave_ssh_ibc, eta, eta_to_m=1.0)
2429 call find_eta(cs%h, cs%tv, g, gv, us, cs%ave_ssh_ibc, eta_to_m=1.0)
2432 if (cs%split)
deallocate(eta)
2435 if (
present(count_calls)) cs%count_calls = count_calls
2436 call mom_sum_output_init(g, us, param_file, dirs%output_directory, &
2437 cs%ntrunc, time_init, cs%sum_output_CSp)
2440 cs%write_IC = save_ic .and. &
2441 .not.((dirs%input_filename(1:1) ==
'r') .and. &
2442 (len_trim(dirs%input_filename) == 1))
2444 if (cs%ensemble_ocean)
then
2445 call init_oda(time, g, gv, cs%odaCS)
2452 call calltree_leave(
"initialize_MOM()")
2453 call cpu_clock_end(id_clock_init)
2455 end subroutine initialize_mom
2458 subroutine finish_mom_initialization(Time, dirs, CS, restart_CSp)
2459 type(time_type),
intent(in) :: time
2471 real,
allocatable :: z_interface(:,:,:)
2474 call cpu_clock_begin(id_clock_init)
2475 call calltree_enter(
"finish_MOM_initialization()")
2478 g => cs%G ; gv => cs%GV ; us => cs%US
2481 call fix_restart_scaling(gv)
2482 call fix_restart_unit_scaling(us)
2485 if (cs%write_IC)
then
2486 allocate(restart_csp_tmp)
2487 restart_csp_tmp = restart_csp
2488 allocate(z_interface(szi_(g),szj_(g),szk_(g)+1))
2489 call find_eta(cs%h, cs%tv, g, gv, us, z_interface, eta_to_m=1.0)
2491 "Interface heights",
"meter", z_grid=
'i')
2493 call save_restart(dirs%output_directory, time, g, &
2494 restart_csp_tmp, filename=cs%IC_file, gv=gv)
2495 deallocate(z_interface)
2496 deallocate(restart_csp_tmp)
2499 call write_energy(cs%u, cs%v, cs%h, cs%tv, time, 0, g, gv, us, &
2500 cs%sum_output_CSp, cs%tracer_flow_CSp)
2502 call calltree_leave(
"finish_MOM_initialization()")
2503 call cpu_clock_end(id_clock_init)
2505 end subroutine finish_mom_initialization
2508 subroutine register_diags(Time, G, GV, IDs, diag)
2509 type(time_type),
intent(in) :: Time
2516 character(len=48) :: thickness_units
2518 thickness_units = get_thickness_units(gv)
2519 if (gv%Boussinesq)
then
2520 h_convert = gv%H_to_m
2522 h_convert = gv%H_to_kg_m2
2526 ids%id_u = register_diag_field(
'ocean_model',
'u_dyn', diag%axesCuL, time, &
2527 'Zonal velocity after the dynamics update',
'm s-1')
2528 ids%id_v = register_diag_field(
'ocean_model',
'v_dyn', diag%axesCvL, time, &
2529 'Meridional velocity after the dynamics update',
'm s-1')
2530 ids%id_h = register_diag_field(
'ocean_model',
'h_dyn', diag%axesTL, time, &
2531 'Layer Thickness after the dynamics update', thickness_units, &
2532 v_extensive=.true., conversion=h_convert)
2533 ids%id_ssh_inst = register_diag_field(
'ocean_model',
'SSH_inst', diag%axesT1, &
2534 time,
'Instantaneous Sea Surface Height',
'm')
2535 end subroutine register_diags
2538 subroutine mom_timing_init(CS)
2541 id_clock_ocean = cpu_clock_id(
'Ocean', grain=clock_component)
2542 id_clock_dynamics = cpu_clock_id(
'Ocean dynamics', grain=clock_subcomponent)
2543 id_clock_thermo = cpu_clock_id(
'Ocean thermodynamics and tracers', grain=clock_subcomponent)
2544 id_clock_other = cpu_clock_id(
'Ocean Other', grain=clock_subcomponent)
2545 id_clock_tracer = cpu_clock_id(
'(Ocean tracer advection)', grain=clock_module_driver)
2546 if (.not.cs%adiabatic) &
2547 id_clock_diabatic = cpu_clock_id(
'(Ocean diabatic driver)', grain=clock_module_driver)
2549 id_clock_continuity = cpu_clock_id(
'(Ocean continuity equation *)', grain=clock_module)
2550 id_clock_bbl_visc = cpu_clock_id(
'(Ocean set BBL viscosity)', grain=clock_module)
2551 id_clock_pass = cpu_clock_id(
'(Ocean message passing *)', grain=clock_module)
2552 id_clock_mom_init = cpu_clock_id(
'(Ocean MOM_initialize_state)', grain=clock_module)
2553 id_clock_pass_init = cpu_clock_id(
'(Ocean init message passing *)', grain=clock_routine)
2554 if (cs%thickness_diffuse) &
2555 id_clock_thick_diff = cpu_clock_id(
'(Ocean thickness diffusion *)', grain=clock_module)
2557 id_clock_ml_restrat = cpu_clock_id(
'(Ocean mixed layer restrat)', grain=clock_module)
2558 id_clock_diagnostics = cpu_clock_id(
'(Ocean collective diagnostics)', grain=clock_module)
2559 id_clock_z_diag = cpu_clock_id(
'(Ocean Z-space diagnostics)', grain=clock_module)
2560 id_clock_ale = cpu_clock_id(
'(Ocean ALE)', grain=clock_module)
2561 if (cs%offline_tracer_mode)
then
2562 id_clock_offline_tracer = cpu_clock_id(
'Ocean offline tracers', grain=clock_subcomponent)
2565 end subroutine mom_timing_init
2576 subroutine set_restart_fields(GV, US, param_file, CS, restart_CSp)
2584 logical :: use_ice_shelf
2585 character(len=48) :: thickness_units, flux_units
2588 thickness_units = get_thickness_units(gv)
2589 flux_units = get_flux_units(gv)
2591 if (
associated(cs%tv%T)) &
2593 "Potential Temperature",
"degC")
2594 if (
associated(cs%tv%S)) &
2599 "Layer Thickness", thickness_units)
2602 "Zonal velocity",
"m s-1", hor_grid=
'Cu')
2605 "Meridional velocity",
"m s-1", hor_grid=
'Cv')
2607 if (
associated(cs%tv%frazil)) &
2609 "Frazil heat flux into ocean",
"J m-2")
2611 if (cs%interp_p_surf)
then
2613 "Previous ocean surface pressure",
"Pa")
2617 "Time average sea surface height",
"meter")
2620 call get_param(param_file,
'',
"ICE_SHELF", use_ice_shelf, default=.false., &
2622 if (use_ice_shelf .and.
associated(cs%Hml))
then
2624 "Mixed layer thickness",
"meter")
2629 "Height unit conversion factor",
"Z meter-1")
2631 "Thickness unit conversion factor",
"H meter-1")
2633 "Length unit conversion factor",
"L meter-1")
2635 "Time unit conversion factor",
"T second-1")
2637 end subroutine set_restart_fields
2641 subroutine adjust_ssh_for_p_atm(tv, G, GV, US, ssh, p_atm, use_EOS)
2646 real,
dimension(SZI_(G),SZJ_(G)),
intent(inout) :: ssh
2647 real,
dimension(:,:),
optional,
pointer :: p_atm
2648 logical,
optional,
intent(in) :: use_EOS
2655 integer :: i, j, is, ie, js, je
2657 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
2658 if (
present(p_atm))
then ;
if (
associated(p_atm))
then
2659 calc_rho =
associated(tv%eqn_of_state)
2660 if (
present(use_eos) .and. calc_rho) calc_rho = use_eos
2663 do j=js,je ;
do i=is,ie
2666 rho_conv, tv%eqn_of_state)
2670 igr0 = 1.0 / (rho_conv * gv%mks_g_Earth)
2671 ssh(i,j) = ssh(i,j) + p_atm(i,j) * igr0
2675 end subroutine adjust_ssh_for_p_atm
2680 subroutine extract_surface_state(CS, sfc_state)
2682 type(
surface),
intent(inout) :: sfc_state
2691 real,
dimension(:,:,:),
pointer :: &
2695 real :: depth(szi_(cs%g))
2702 real :: delt(szi_(cs%g))
2703 logical :: use_temperature
2704 integer :: i, j, k, is, ie, js, je, nz, numberoferrors, ig, jg
2705 integer :: isd, ied, jsd, jed
2706 integer :: iscb, iecb, jscb, jecb, isdb, iedb, jsdb, jedb
2707 logical :: localerror
2708 character(240) :: msg
2710 call calltree_enter(
"extract_surface_state(), MOM.F90")
2711 g => cs%G ; gv => cs%GV
2712 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
2713 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
2714 iscb = g%iscB ; iecb = g%iecB; jscb = g%jscB ; jecb = g%jecB
2715 isdb = g%isdB ; iedb = g%iedB; jsdb = g%jsdB ; jedb = g%jedB
2716 u => cs%u ; v => cs%v ; h => cs%h
2718 use_temperature =
associated(cs%tv%T)
2720 if (.not.sfc_state%arrays_allocated)
then
2723 call allocate_surface_state(sfc_state, g, use_temperature, do_integrals=.true.)
2725 sfc_state%frazil => cs%tv%frazil
2726 sfc_state%TempxPmE => cs%tv%TempxPmE
2727 sfc_state%internal_heat => cs%tv%internal_heat
2728 sfc_state%T_is_conT = cs%tv%T_is_conT
2729 sfc_state%S_is_absS = cs%tv%S_is_absS
2730 if (
associated(cs%visc%taux_shelf)) sfc_state%taux_shelf => cs%visc%taux_shelf
2731 if (
associated(cs%visc%tauy_shelf)) sfc_state%tauy_shelf => cs%visc%tauy_shelf
2733 do j=js,je ;
do i=is,ie
2734 sfc_state%sea_lev(i,j) = cs%ave_ssh_ibc(i,j)
2738 if (
associated(cs%Hml))
then
2739 do j=js,je ;
do i=is,ie
2740 sfc_state%Hml(i,j) = cs%Hml(i,j)
2744 if (cs%Hmix < 0.0)
then
2745 if (use_temperature)
then ;
do j=js,je ;
do i=is,ie
2746 sfc_state%SST(i,j) = cs%tv%T(i,j,1)
2747 sfc_state%SSS(i,j) = cs%tv%S(i,j,1)
2748 enddo ;
enddo ;
endif
2749 do j=js,je ;
do i=is-1,ie
2750 sfc_state%u(i,j) = u(i,j,1)
2752 do j=js-1,je ;
do i=is,ie
2753 sfc_state%v(i,j) = v(i,j,1)
2765 if (use_temperature)
then
2766 sfc_state%SST(i,j) = 0.0 ; sfc_state%SSS(i,j) = 0.0
2768 sfc_state%sfc_density(i,j) = 0.0
2772 do k=1,nz ;
do i=is,ie
2773 if (depth(i) + h(i,j,k)*gv%H_to_Z < depth_ml)
then
2774 dh = h(i,j,k)*gv%H_to_Z
2775 elseif (depth(i) < depth_ml)
then
2776 dh = depth_ml - depth(i)
2780 if (use_temperature)
then
2781 sfc_state%SST(i,j) = sfc_state%SST(i,j) + dh * cs%tv%T(i,j,k)
2782 sfc_state%SSS(i,j) = sfc_state%SSS(i,j) + dh * cs%tv%S(i,j,k)
2784 sfc_state%sfc_density(i,j) = sfc_state%sfc_density(i,j) + dh * gv%Rlay(k)
2786 depth(i) = depth(i) + dh
2790 if (depth(i) < gv%H_subroundoff*gv%H_to_Z) &
2791 depth(i) = gv%H_subroundoff*gv%H_to_Z
2792 if (use_temperature)
then
2793 sfc_state%SST(i,j) = sfc_state%SST(i,j) / depth(i)
2794 sfc_state%SSS(i,j) = sfc_state%SSS(i,j) / depth(i)
2796 sfc_state%sfc_density(i,j) = sfc_state%sfc_density(i,j) / depth(i)
2807 if (cs%Hmix_UV>0.)
then
2810 depth_ml = cs%Hmix_UV
2815 sfc_state%v(i,j) = 0.0
2817 do k=1,nz ;
do i=is,ie
2818 hv = 0.5 * (h(i,j,k) + h(i,j+1,k)) * gv%H_to_Z
2819 if (depth(i) + hv < depth_ml)
then
2821 elseif (depth(i) < depth_ml)
then
2822 dh = depth_ml - depth(i)
2826 sfc_state%v(i,j) = sfc_state%v(i,j) + dh * v(i,j,k)
2827 depth(i) = depth(i) + dh
2831 if (depth(i) < gv%H_subroundoff*gv%H_to_Z) &
2832 depth(i) = gv%H_subroundoff*gv%H_to_Z
2833 sfc_state%v(i,j) = sfc_state%v(i,j) / depth(i)
2841 sfc_state%u(i,j) = 0.0
2843 do k=1,nz ;
do i=is-1,ie
2844 hu = 0.5 * (h(i,j,k) + h(i+1,j,k)) * gv%H_to_Z
2845 if (depth(i) + hu < depth_ml)
then
2847 elseif (depth(i) < depth_ml)
then
2848 dh = depth_ml - depth(i)
2852 sfc_state%u(i,j) = sfc_state%u(i,j) + dh * u(i,j,k)
2853 depth(i) = depth(i) + dh
2857 if (depth(i) < gv%H_subroundoff*gv%H_to_Z) &
2858 depth(i) = gv%H_subroundoff*gv%H_to_Z
2859 sfc_state%u(i,j) = sfc_state%u(i,j) / depth(i)
2863 do j=js,je ;
do i=is-1,ie
2864 sfc_state%u(i,j) = u(i,j,1)
2866 do j=js-1,je ;
do i=is,ie
2867 sfc_state%v(i,j) = v(i,j,1)
2873 if (
allocated(sfc_state%melt_potential))
then
2881 do k=1,nz ;
do i=is,ie
2882 depth_ml = min(cs%HFrz,cs%visc%MLD(i,j))
2883 if (depth(i) + h(i,j,k)*gv%H_to_m < depth_ml)
then
2884 dh = h(i,j,k)*gv%H_to_m
2885 elseif (depth(i) < depth_ml)
then
2886 dh = depth_ml - depth(i)
2893 depth(i) = depth(i) + dh
2894 delt(i) = delt(i) + dh * (cs%tv%T(i,j,k) - t_freeze)
2899 sfc_state%melt_potential(i,j) = 0.0
2901 if (g%mask2dT(i,j)>0.)
then
2903 sfc_state%melt_potential(i,j) = cs%tv%C_p * cs%GV%Rho0 * delt(i)
2909 if (
allocated(sfc_state%salt_deficit) .and.
associated(cs%tv%salt_deficit))
then
2911 do j=js,je ;
do i=is,ie
2913 sfc_state%salt_deficit(i,j) = 1000.0 * cs%tv%salt_deficit(i,j)
2917 if (
allocated(sfc_state%ocean_mass) .and.
allocated(sfc_state%ocean_heat) .and. &
2918 allocated(sfc_state%ocean_salt))
then
2920 do j=js,je ;
do i=is,ie
2921 sfc_state%ocean_mass(i,j) = 0.0
2922 sfc_state%ocean_heat(i,j) = 0.0 ; sfc_state%ocean_salt(i,j) = 0.0
2925 do j=js,je ;
do k=1,nz;
do i=is,ie
2926 mass = gv%H_to_kg_m2*h(i,j,k)
2927 sfc_state%ocean_mass(i,j) = sfc_state%ocean_mass(i,j) + mass
2928 sfc_state%ocean_heat(i,j) = sfc_state%ocean_heat(i,j) + mass*cs%tv%T(i,j,k)
2929 sfc_state%ocean_salt(i,j) = sfc_state%ocean_salt(i,j) + &
2930 mass * (1.0e-3*cs%tv%S(i,j,k))
2931 enddo ;
enddo ;
enddo
2933 if (
allocated(sfc_state%ocean_mass))
then
2935 do j=js,je ;
do i=is,ie ; sfc_state%ocean_mass(i,j) = 0.0 ;
enddo ;
enddo
2937 do j=js,je ;
do k=1,nz ;
do i=is,ie
2938 sfc_state%ocean_mass(i,j) = sfc_state%ocean_mass(i,j) + gv%H_to_kg_m2*h(i,j,k)
2939 enddo ;
enddo ;
enddo
2941 if (
allocated(sfc_state%ocean_heat))
then
2943 do j=js,je ;
do i=is,ie ; sfc_state%ocean_heat(i,j) = 0.0 ;
enddo ;
enddo
2945 do j=js,je ;
do k=1,nz ;
do i=is,ie
2946 mass = gv%H_to_kg_m2*h(i,j,k)
2947 sfc_state%ocean_heat(i,j) = sfc_state%ocean_heat(i,j) + mass*cs%tv%T(i,j,k)
2948 enddo ;
enddo ;
enddo
2950 if (
allocated(sfc_state%ocean_salt))
then
2952 do j=js,je ;
do i=is,ie ; sfc_state%ocean_salt(i,j) = 0.0 ;
enddo ;
enddo
2954 do j=js,je ;
do k=1,nz ;
do i=is,ie
2955 mass = gv%H_to_kg_m2*h(i,j,k)
2956 sfc_state%ocean_salt(i,j) = sfc_state%ocean_salt(i,j) + &
2957 mass * (1.0e-3*cs%tv%S(i,j,k))
2958 enddo ;
enddo ;
enddo
2962 if (
associated(cs%tracer_flow_CSp))
then
2963 call call_tracer_surface_state(sfc_state, h, g, cs%tracer_flow_CSp)
2966 if (cs%check_bad_sfc_vals)
then
2968 do j=js,je;
do i=is,ie
2969 if (g%mask2dT(i,j)>0.)
then
2970 bathy_m = cs%US%Z_to_m * g%bathyT(i,j)
2971 localerror = sfc_state%sea_lev(i,j)<=-bathy_m &
2972 .or. sfc_state%sea_lev(i,j)>= cs%bad_val_ssh_max &
2973 .or. sfc_state%sea_lev(i,j)<=-cs%bad_val_ssh_max &
2974 .or. sfc_state%sea_lev(i,j) + bathy_m < cs%bad_val_col_thick
2975 if (use_temperature) localerror = localerror &
2976 .or. sfc_state%SSS(i,j)<0. &
2977 .or. sfc_state%SSS(i,j)>=cs%bad_val_sss_max &
2978 .or. sfc_state%SST(i,j)< cs%bad_val_sst_min &
2979 .or. sfc_state%SST(i,j)>=cs%bad_val_sst_max
2980 if (localerror)
then
2981 numberoferrors=numberoferrors+1
2982 if (numberoferrors<9)
then
2983 ig = i + g%HI%idg_offset
2984 jg = j + g%HI%jdg_offset
2985 if (use_temperature)
then
2986 write(msg(1:240),
'(2(a,i4,x),4(a,f8.3,x),8(a,es11.4,x))') &
2987 'Extreme surface sfc_state detected: i=',ig,
'j=',jg, &
2988 'lon=',g%geoLonT(i,j),
'lat=',g%geoLatT(i,j), &
2989 'x=',g%gridLonT(ig),
'y=',g%gridLatT(jg), &
2990 'D=',bathy_m,
'SSH=',sfc_state%sea_lev(i,j), &
2991 'SST=',sfc_state%SST(i,j),
'SSS=',sfc_state%SSS(i,j), &
2992 'U-=',sfc_state%u(i-1,j),
'U+=',sfc_state%u(i,j), &
2993 'V-=',sfc_state%v(i,j-1),
'V+=',sfc_state%v(i,j)
2995 write(msg(1:240),
'(2(a,i4,x),4(a,f8.3,x),6(a,es11.4))') &
2996 'Extreme surface sfc_state detected: i=',ig,
'j=',jg, &
2997 'lon=',g%geoLonT(i,j),
'lat=',g%geoLatT(i,j), &
2998 'x=',g%gridLonT(i),
'y=',g%gridLatT(j), &
2999 'D=',bathy_m,
'SSH=',sfc_state%sea_lev(i,j), &
3000 'U-=',sfc_state%u(i-1,j),
'U+=',sfc_state%u(i,j), &
3001 'V-=',sfc_state%v(i,j-1),
'V+=',sfc_state%v(i,j)
3003 call mom_error(warning, trim(msg), all_print=.true.)
3004 elseif (numberoferrors==9)
then
3005 call mom_error(warning,
'There were more unreported extreme events!', all_print=.true.)
3010 call sum_across_pes(numberoferrors)
3011 if (numberoferrors>0)
then
3012 write(msg(1:240),
'(3(a,i9,x))')
'There were a total of ',numberoferrors, &
3013 'locations detected with extreme surface values!'
3014 call mom_error(fatal, trim(msg))
3018 if (cs%debug)
call mom_surface_chksum(
"Post extract_sfc", sfc_state, g)
3020 call calltree_leave(
"extract_surface_sfc_state()")
3021 end subroutine extract_surface_state
3024 function mom_state_is_synchronized(CS, adv_dyn)
result(in_synch)
3026 logical,
optional,
intent(in) :: adv_dyn
3033 adv_only = .false. ;
if (
present(adv_dyn)) adv_only = adv_dyn
3036 in_synch = (cs%t_dyn_rel_adv == 0.0)
3038 in_synch = ((cs%t_dyn_rel_adv == 0.0) .and. (cs%t_dyn_rel_thermo == 0.0))
3041 end function mom_state_is_synchronized
3045 subroutine get_mom_state_elements(CS, G, GV, US, C_p, use_temp)
3048 optional,
pointer :: g
3050 optional,
pointer :: gv
3052 optional,
pointer :: us
3053 real,
optional,
intent(out) :: c_p
3054 logical,
optional,
intent(out) :: use_temp
3056 if (
present(g)) g => cs%G
3057 if (
present(gv)) gv => cs%GV
3058 if (
present(us)) us => cs%US
3059 if (
present(c_p)) c_p = cs%tv%C_p
3060 if (
present(use_temp)) use_temp =
associated(cs%tv%T)
3061 end subroutine get_mom_state_elements
3064 subroutine get_ocean_stocks(CS, mass, heat, salt, on_PE_only)
3066 real,
optional,
intent(out) :: heat
3067 real,
optional,
intent(out) :: salt
3068 real,
optional,
intent(out) :: mass
3069 logical,
optional,
intent(in) :: on_pe_only
3071 if (
present(mass)) &
3072 mass = global_mass_integral(cs%h, cs%G, cs%GV, on_pe_only=on_pe_only)
3073 if (
present(heat)) &
3074 heat = cs%tv%C_p * global_mass_integral(cs%h, cs%G, cs%GV, cs%tv%T, on_pe_only=on_pe_only)
3075 if (
present(salt)) &
3076 salt = 1.0e-3 * global_mass_integral(cs%h, cs%G, cs%GV, cs%tv%S, on_pe_only=on_pe_only)
3078 end subroutine get_ocean_stocks
3081 subroutine mom_end(CS)
3084 if (cs%use_ALE_algorithm)
call ale_end(cs%ALE_CSp)
3086 dealloc_(cs%u) ; dealloc_(cs%v) ; dealloc_(cs%h)
3087 dealloc_(cs%uh) ; dealloc_(cs%vh)
3089 if (
associated(cs%tv%T))
then
3090 dealloc_(cs%T) ; cs%tv%T => null() ; dealloc_(cs%S) ; cs%tv%S => null()
3092 if (
associated(cs%tv%frazil))
deallocate(cs%tv%frazil)
3093 if (
associated(cs%tv%salt_deficit))
deallocate(cs%tv%salt_deficit)
3094 if (
associated(cs%Hml))
deallocate(cs%Hml)
3096 call tracer_advect_end(cs%tracer_adv_CSp)
3097 call tracer_hor_diff_end(cs%tracer_diff_CSp)
3098 call tracer_registry_end(cs%tracer_Reg)
3099 call tracer_flow_control_end(cs%tracer_flow_CSp)
3101 call diabatic_driver_end(cs%diabatic_CSp)
3103 if (cs%offline_tracer_mode)
call offline_transport_end(cs%offline_CSp)
3105 dealloc_(cs%uhtr) ; dealloc_(cs%vhtr)
3107 call end_dyn_split_rk2(cs%dyn_split_RK2_CSp)
3108 elseif (cs%use_RK2)
then
3109 call end_dyn_unsplit_rk2(cs%dyn_unsplit_RK2_CSp)
3111 call end_dyn_unsplit(cs%dyn_unsplit_CSp)
3113 dealloc_(cs%ave_ssh_ibc) ; dealloc_(cs%ssh_rint)
3114 if (
associated(cs%update_OBC_CSp))
call obc_register_end(cs%update_OBC_CSp)
3116 call verticalgridend(cs%GV)
3117 call unit_scaling_end(cs%US)
3118 call mom_grid_end(cs%G)
3122 end subroutine mom_end