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
26 use mom_domains,
only : to_north, to_east, to_south, to_west
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
134 implicit none ;
private
136 #include <MOM_memory.h>
146 integer :: id_u = -1, id_v = -1, id_h = -1
148 integer :: id_ssh_inst = -1
154 real allocable_,
dimension(NIMEM_,NJMEM_,NKMEM_) :: &
155 h, & !< layer thickness [H ~> m or kg m-2]
156 t, & !< potential temperature [degC]
158 real allocable_,
dimension(NIMEMB_PTR_,NJMEM_,NKMEM_) :: &
159 u, & !< zonal velocity component [L T-1 ~> m s-1]
160 uh, & !< uh = u * h * dy at u grid points [H L2 T-1 ~> m3 s-1 or kg s-1]
162 real allocable_,
dimension(NIMEM_,NJMEMB_PTR_,NKMEM_) :: &
163 v, & !< meridional velocity [L T-1 ~> m s-1]
164 vh, & !< vh = v * h * dx at v grid points [H L2 T-1 ~> m3 s-1 or kg s-1]
166 real allocable_,
dimension(NIMEM_,NJMEM_) :: ssh_rint
168 real allocable_,
dimension(NIMEM_,NJMEM_) :: ave_ssh_ibc
171 real allocable_,
dimension(NIMEM_,NJMEM_) :: eta_av_bc
174 real,
dimension(:,:),
pointer :: &
176 real :: time_in_cycle
179 real :: time_in_thermo_cycle
188 real :: t_dyn_rel_adv
191 real :: t_dyn_rel_thermo
195 real :: t_dyn_rel_diag
197 integer :: ndyn_per_adv = 0
207 logical :: diabatic_first
209 logical :: use_ale_algorithm
212 logical :: offline_tracer_mode = .false.
216 type(time_type),
pointer :: time
219 logical :: thermo_spans_coupling
221 integer :: nstep_tot = 0
223 logical :: count_calls = .false.
229 logical :: do_dynamics
235 logical :: thickness_diffuse_first
236 logical :: mixedlayer_restrat
239 real :: dtbt_reset_period
242 type(time_type) :: dtbt_reset_interval
243 type(time_type) :: dtbt_reset_time
246 real,
dimension(:,:,:),
pointer :: &
247 h_pre_dyn => null(), &
248 t_pre_dyn => null(), &
254 real,
dimension(:,:,:),
pointer :: &
258 logical :: interp_p_surf
261 logical :: p_surf_prev_set
264 real,
dimension(:,:),
pointer :: &
265 p_surf_prev => null(), &
266 p_surf_begin => null(), &
271 character(len=120) :: ic_file
274 logical :: calc_rho_for_sea_lev
289 logical :: check_bad_sfc_vals
290 real :: bad_val_ssh_max
291 real :: bad_val_sst_max
292 real :: bad_val_sst_min
293 real :: bad_val_sss_max
294 real :: bad_val_col_thick
342 type(
ale_cs),
pointer :: ale_csp => null()
353 logical :: ensemble_ocean
365 public allocate_surface_state, deallocate_surface_state
396 subroutine step_mom(forces, fluxes, sfc_state, Time_start, time_int_in, CS, &
397 Waves, do_dynamics, do_thermodynamics, start_cycle, &
398 end_cycle, cycle_length, reset_therm)
400 type(
forcing),
intent(inout) :: fluxes
402 type(
surface),
intent(inout) :: sfc_state
403 type(time_type),
intent(in) :: time_start
404 real,
intent(in) :: time_int_in
407 optional,
pointer :: waves
408 logical,
optional,
intent(in) :: do_dynamics
410 logical,
optional,
intent(in) :: do_thermodynamics
412 logical,
optional,
intent(in) :: start_cycle
415 logical,
optional,
intent(in) :: end_cycle
418 real,
optional,
intent(in) :: cycle_length
420 logical,
optional,
intent(in) :: reset_therm
433 integer :: i, j, k, is, ie, js, je, isq, ieq, jsq, jeq, nz, n
434 integer :: isd, ied, jsd, jed, isdb, iedb, jsdb, jedb
436 real :: time_interval
440 real :: dt_therm_here
442 real :: wt_end, wt_beg
446 real :: rel_time = 0.0
450 logical :: do_advection
451 logical :: do_calc_bbl
452 logical :: thermo_does_span_coupling
456 logical :: cycle_start
460 logical :: therm_reset
462 real,
dimension(SZI_(CS%G),SZJ_(CS%G)) :: &
465 real,
dimension(:,:,:),
pointer :: &
469 real,
dimension(:,:),
pointer :: &
473 type(time_type) :: time_local, end_time_thermo, time_temp
474 type(group_pass_type) :: pass_tau_ustar_psurf
475 logical :: showcalltree
477 g => cs%G ; gv => cs%GV ; us => cs%US
478 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
479 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
480 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
481 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
482 u => cs%u ; v => cs%v ; h => cs%h
484 time_interval = us%s_to_T*time_int_in
485 do_dyn = .true. ;
if (
present(do_dynamics)) do_dyn = do_dynamics
486 do_thermo = .true. ;
if (
present(do_thermodynamics)) do_thermo = do_thermodynamics
487 if (.not.(do_dyn .or. do_thermo))
call mom_error(fatal,
"Step_MOM: "//&
488 "Both do_dynamics and do_thermodynamics are false, which makes no sense.")
489 cycle_start = .true. ;
if (
present(start_cycle)) cycle_start = start_cycle
490 cycle_end = .true. ;
if (
present(end_cycle)) cycle_end = end_cycle
491 cycle_time = time_interval ;
if (
present(cycle_length)) cycle_time = us%s_to_T*cycle_length
492 therm_reset = cycle_start ;
if (
present(reset_therm)) therm_reset = reset_therm
498 call mom_state_chksum(
"Beginning of step_MOM ", u, v, h, cs%uh, cs%vh, g, gv, us)
501 showcalltree = calltree_showquery()
508 if (time_interval > cs%dt) n_max = ceiling(time_interval/cs%dt - 0.001)
509 dt = time_interval / real(n_max)
510 thermo_does_span_coupling = (cs%thermo_spans_coupling .and. &
511 (cs%dt_therm > 1.5*cycle_time))
512 if (thermo_does_span_coupling)
then
514 dt_therm = cycle_time * floor(cs%dt_therm / cycle_time + 0.001)
515 ntstep = floor(dt_therm/dt + 0.001)
516 elseif (.not.do_thermo)
then
517 dt_therm = cs%dt_therm
518 if (
present(cycle_length)) dt_therm = min(cs%dt_therm, us%s_to_T*cycle_length)
521 ntstep = max(1, min(n_max, floor(cs%dt_therm/dt + 0.001)))
525 if (
associated(forces%p_surf)) p_surf => forces%p_surf
526 if (.not.
associated(forces%p_surf)) cs%interp_p_surf = .false.
531 if (
associated(forces%ustar)) &
533 if (
associated(forces%p_surf)) &
535 if (g%nonblocking_updates)
then
538 call do_group_pass(pass_tau_ustar_psurf, g%Domain)
544 if ((time_interval > cs%dt_therm) .and. (cs%dt_therm > 0.0)) &
545 n_max = ceiling(time_interval/cs%dt_therm - 0.001)
547 dt = time_interval / real(n_max)
548 dt_therm = dt ; ntstep = 1
549 if (
associated(fluxes%p_surf)) p_surf => fluxes%p_surf
554 if (therm_reset)
then
555 cs%time_in_thermo_cycle = 0.0
556 if (
associated(cs%tv%frazil)) cs%tv%frazil(:,:) = 0.0
557 if (
associated(cs%tv%salt_deficit)) cs%tv%salt_deficit(:,:) = 0.0
558 if (
associated(cs%tv%TempxPmE)) cs%tv%TempxPmE(:,:) = 0.0
559 if (
associated(cs%tv%internal_heat)) cs%tv%internal_heat(:,:) = 0.0
562 if (cycle_start)
then
563 cs%time_in_cycle = 0.0
564 do j=js,je ;
do i=is,ie ; cs%ssh_rint(i,j) = 0.0 ;
enddo ;
enddo
566 if (
associated(cs%VarMix))
then
567 call enable_averages(cycle_time, time_start + real_to_time(us%T_to_s*cycle_time), cs%diag)
570 call disable_averaging(cs%diag)
575 if (g%nonblocking_updates) &
578 if (cs%interp_p_surf)
then
579 if (.not.
associated(cs%p_surf_end))
allocate(cs%p_surf_end(isd:ied,jsd:jed))
580 if (.not.
associated(cs%p_surf_begin))
allocate(cs%p_surf_begin(isd:ied,jsd:jed))
581 if (.not.cs%p_surf_prev_set)
then
582 do j=jsd,jed ;
do i=isd,ied
583 cs%p_surf_prev(i,j) = forces%p_surf(i,j)
585 cs%p_surf_prev_set = .true.
588 cs%p_surf_end => forces%p_surf
591 if (cs%UseWaves)
then
593 call enable_averages(time_interval, time_start + real_to_time(us%T_to_s*time_interval), cs%diag)
595 call disable_averaging(cs%diag)
607 if (do_dyn)
call check_redundant(
"Before steps ", forces%taux, forces%tauy, g)
613 rel_time = rel_time + dt
615 cs%Time = time_start + real_to_time(us%T_to_s*(rel_time - 0.5*dt))
617 time_local = time_start + real_to_time(us%T_to_s*rel_time)
619 if (showcalltree)
call calltree_enter(
"DT cycles (step_MOM) n=",n)
623 if (cs%diabatic_first .and. (cs%t_dyn_rel_adv==0.0) .and. do_thermo)
then
625 if (.not.do_dyn)
then
627 elseif (thermo_does_span_coupling)
then
629 if ((fluxes%dt_buoy_accum > 0.0) .and. (dtdia > time_interval) .and. &
630 (abs(fluxes%dt_buoy_accum - dtdia) > 1e-6*dtdia))
then
631 call mom_error(fatal,
"step_MOM: Mismatch between long thermodynamic "//&
632 "timestep and time over which buoyancy fluxes have been accumulated.")
634 call mom_error(fatal,
"MOM is not yet set up to have restarts that work "//&
635 "with THERMO_SPANS_COUPLING and DIABATIC_FIRST.")
637 dtdia = dt*min(ntstep,n_max-(n-1))
640 end_time_thermo = time_local
644 cs%Time = cs%Time + real_to_time(0.5*us%T_to_s*(dtdia-dt))
647 end_time_thermo = time_local + real_to_time(us%T_to_s*(dtdia-dt))
651 call step_mom_thermo(cs, g, gv, us, u, v, h, cs%tv, fluxes, dtdia, &
652 end_time_thermo, .true., waves=waves)
653 cs%time_in_thermo_cycle = cs%time_in_thermo_cycle + dtdia
656 cs%t_dyn_rel_thermo = -dtdia
660 cs%Time = time_start + real_to_time(us%T_to_s*(rel_time - 0.5*dt))
667 if (cs%ndyn_per_adv == 0 .and. cs%t_dyn_rel_adv == 0.)
then
669 cs%ndyn_per_adv = cs%ndyn_per_adv + 1
673 if (
associated(cs%u_prev) .and.
associated(cs%v_prev))
then
674 do k=1,nz ;
do j=jsd,jed ;
do i=isdb,iedb
675 cs%u_prev(i,j,k) = u(i,j,k)
676 enddo ;
enddo ;
enddo
677 do k=1,nz ;
do j=jsdb,jedb ;
do i=isd,ied
678 cs%v_prev(i,j,k) = v(i,j,k)
679 enddo ;
enddo ;
enddo
682 dt_therm_here = dt_therm
683 if (do_thermo .and. do_dyn .and. .not.thermo_does_span_coupling) &
684 dt_therm_here = dt*min(ntstep, n_max-n+1)
690 if ((cs%t_dyn_rel_adv == 0.0) .or. (n==1)) &
691 bbl_time_int = max(dt, min(dt_therm - cs%t_dyn_rel_adv, dt*(1+n_max-n)) )
693 if ((cs%t_dyn_rel_adv == 0.0) .or. ((n==1) .and. cycle_start)) &
694 bbl_time_int = min(dt_therm, cycle_time)
697 if (cs%interp_p_surf)
then
698 wt_end = real(n) / real(n_max)
699 wt_beg = real(n-1) / real(n_max)
700 do j=jsd,jed ;
do i=isd,ied
701 cs%p_surf_end(i,j) = wt_end * forces%p_surf(i,j) + &
702 (1.0-wt_end) * cs%p_surf_prev(i,j)
703 cs%p_surf_begin(i,j) = wt_beg * forces%p_surf(i,j) + &
704 (1.0-wt_beg) * cs%p_surf_prev(i,j)
709 dt_therm_here, bbl_time_int, cs, &
710 time_local, waves=waves)
715 if (thermo_does_span_coupling .or. .not.do_thermo)
then
716 do_advection = (cs%t_dyn_rel_adv + 0.5*dt > dt_therm)
718 do_advection = ((mod(n,ntstep) == 0) .or. (n==n_max))
721 if (do_advection)
then
724 if (cs%diabatic_first .and. abs(cs%t_dyn_rel_thermo) > 1e-6*dt)
call mom_error(fatal, &
725 "step_MOM: Mismatch between the dynamics and diabatic times "//&
726 "with DIABATIC_FIRST.")
732 if ((cs%t_dyn_rel_adv==0.0) .and. do_thermo .and. (.not.cs%diabatic_first))
then
734 dtdia = cs%t_dyn_rel_thermo
737 if ((cs%t_dyn_rel_thermo==0.0) .and. .not.do_dyn) dtdia = dt
739 if (cs%thermo_spans_coupling .and. (cs%dt_therm > 1.5*cycle_time) .and. &
740 (abs(dt_therm - dtdia) > 1e-6*dt_therm))
then
741 call mom_error(fatal,
"step_MOM: Mismatch between dt_therm and dtdia "//&
742 "before call to diabatic.")
747 if (dtdia > dt) cs%Time = cs%Time - real_to_time(0.5*us%T_to_s*(dtdia-dt))
750 call step_mom_thermo(cs, g, gv, us, u, v, h, cs%tv, fluxes, dtdia, &
751 time_local, .false., waves=waves)
752 cs%time_in_thermo_cycle = cs%time_in_thermo_cycle + dtdia
754 if ((cs%t_dyn_rel_thermo==0.0) .and. .not.do_dyn)
then
756 cs%t_dyn_rel_thermo = -dtdia
758 cs%t_dyn_rel_thermo = 0.0
762 cs%Time = time_start + real_to_time(us%T_to_s*(rel_time - 0.5*dt))
769 cs%time_in_cycle = cs%time_in_cycle + dt
770 call find_eta(h, cs%tv, g, gv, us, ssh, cs%eta_av_bc, eta_to_m=1.0)
771 do j=js,je ;
do i=is,ie
772 cs%ssh_rint(i,j) = cs%ssh_rint(i,j) + dt*ssh(i,j)
774 if (cs%IDs%id_ssh_inst > 0)
call post_data(cs%IDs%id_ssh_inst, ssh, cs%diag)
785 call enable_averages(cs%t_dyn_rel_diag, time_local, cs%diag)
786 call calculate_diagnostic_fields(u, v, h, cs%uh, cs%vh, cs%tv, cs%ADp, &
787 cs%CDp, p_surf, cs%t_dyn_rel_diag, cs%diag_pre_sync,&
788 g, gv, us, cs%diagnostics_CSp)
789 call post_tracer_diagnostics(cs%Tracer_reg, h, cs%diag_pre_sync, cs%diag, g, gv, cs%t_dyn_rel_diag)
791 if (showcalltree)
call calltree_waypoint(
"finished calculate_diagnostic_fields (step_MOM)")
792 call disable_averaging(cs%diag)
793 cs%t_dyn_rel_diag = 0.0
798 if (do_dyn .and. .not.cs%count_calls) cs%nstep_tot = cs%nstep_tot + 1
803 if (cs%count_calls .and. cycle_start) cs%nstep_tot = cs%nstep_tot + 1
807 if (cs%time_in_cycle > 0.0)
then
808 i_wt_ssh = 1.0/cs%time_in_cycle
809 do j=js,je ;
do i=is,ie
810 ssh(i,j) = cs%ssh_rint(i,j)*i_wt_ssh
811 cs%ave_ssh_ibc(i,j) = ssh(i,j)
815 cs%calc_rho_for_sea_lev)
816 elseif (do_thermo)
then
818 cs%calc_rho_for_sea_lev)
822 if (do_dyn .and. cs%interp_p_surf)
then ;
do j=jsd,jed ;
do i=isd,ied
823 cs%p_surf_prev(i,j) = forces%p_surf(i,j)
824 enddo ;
enddo ;
endif
826 if (cs%ensemble_ocean)
then
832 call oda(cs%Time,cs%odaCS)
835 if (showcalltree)
call calltree_waypoint(
"calling extract_surface_state (step_MOM)")
841 if (cs%time_in_cycle > 0.0)
then
842 call enable_averages(cs%time_in_cycle, time_local, cs%diag)
843 call post_surface_dyn_diags(cs%sfc_IDs, g, cs%diag, sfc_state, ssh)
845 if (cs%time_in_thermo_cycle > 0.0)
then
846 call enable_averages(cs%time_in_thermo_cycle, time_local, cs%diag)
847 call post_surface_thermo_diags(cs%sfc_IDs, g, gv, us, cs%diag, cs%time_in_thermo_cycle, &
848 sfc_state, cs%tv, ssh, cs%ave_ssh_ibc)
850 call disable_averaging(cs%diag)
855 if (do_thermo .and. fluxes%fluxes_used) &
856 call accumulate_net_input(fluxes, sfc_state, cs%tv, fluxes%dt_buoy_accum, &
857 g, us, cs%sum_output_CSp)
860 call write_energy(cs%u, cs%v, cs%h, cs%tv, time_local, cs%nstep_tot, &
861 g, gv, us, cs%sum_output_CSp, cs%tracer_flow_CSp, &
862 dt_forcing=real_to_time(us%T_to_s*time_interval) )
873 bbl_time_int, CS, Time_local, Waves)
875 real,
dimension(:,:),
pointer :: p_surf_begin
878 real,
dimension(:,:),
pointer :: p_surf_end
881 real,
intent(in) :: dt
882 real,
intent(in) :: dt_thermo
884 real,
intent(in) :: bbl_time_int
888 type(time_type),
intent(in) :: Time_local
890 optional,
pointer :: Waves
900 real,
dimension(:,:,:),
pointer :: &
907 logical :: showCallTree
909 integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz
910 integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB
912 g => cs%G ; gv => cs%GV ; us => cs%US ; ids => cs%IDs
913 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
914 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
915 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
916 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
917 u => cs%u ; v => cs%v ; h => cs%h
918 showcalltree = calltree_showquery()
922 if ((cs%t_dyn_rel_adv == 0.0) .and. cs%thickness_diffuse .and. cs%thickness_diffuse_first)
then
924 call enable_averages(dt_thermo, time_local+real_to_time(us%T_to_s*(dt_thermo-dt)), cs%diag)
926 if (
associated(cs%VarMix)) &
927 call calc_slope_functions(h, cs%tv, dt, g, gv, us, cs%VarMix)
929 cs%MEKE, cs%VarMix, cs%CDp, cs%thickness_diffuse_CSp)
932 call disable_averaging(cs%diag)
933 if (showcalltree)
call calltree_waypoint(
"finished thickness_diffuse_first (step_MOM)")
937 call diag_update_remap_grids(cs%diag)
941 if (bbl_time_int > 0.0)
then
942 call enable_averages(bbl_time_int, &
943 time_local + real_to_time(us%T_to_s*(bbl_time_int-dt)), cs%diag)
946 call set_viscous_bbl(cs%u(:,:,:), cs%v(:,:,:), cs%h, cs%tv, cs%visc, g, gv, us, &
947 cs%set_visc_CSp, symmetrize=.true.)
950 call disable_averaging(cs%diag)
954 if (cs%do_dynamics .and. cs%split)
then
959 if (cs%dtbt_reset_period == 0.0) calc_dtbt = .true.
960 if (cs%dtbt_reset_period > 0.0)
then
961 if (time_local >= cs%dtbt_reset_time)
then
963 cs%dtbt_reset_time = cs%dtbt_reset_time + cs%dtbt_reset_interval
967 call step_mom_dyn_split_rk2(u, v, h, cs%tv, cs%visc, time_local, dt, forces, &
968 p_surf_begin, p_surf_end, cs%uh, cs%vh, cs%uhtr, cs%vhtr, &
969 cs%eta_av_bc, g, gv, us, cs%dyn_split_RK2_CSp, calc_dtbt, cs%VarMix, &
970 cs%MEKE, cs%thickness_diffuse_CSp, waves=waves)
971 if (showcalltree)
call calltree_waypoint(
"finished step_MOM_dyn_split (step_MOM)")
973 elseif (cs%do_dynamics)
then
982 call step_mom_dyn_unsplit_rk2(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_RK2_CSp, cs%VarMix, cs%MEKE)
986 call step_mom_dyn_unsplit(u, v, h, cs%tv, cs%visc, time_local, dt, forces, &
987 p_surf_begin, p_surf_end, cs%uh, cs%vh, cs%uhtr, cs%vhtr, &
988 cs%eta_av_bc, g, gv, us, cs%dyn_unsplit_CSp, cs%VarMix, cs%MEKE, waves=waves)
990 if (showcalltree)
call calltree_waypoint(
"finished step_MOM_dyn_unsplit (step_MOM)")
994 if (cs%thickness_diffuse .and. .not.cs%thickness_diffuse_first)
then
997 if (cs%debug)
call hchksum(h,
"Pre-thickness_diffuse h", g%HI, haloshift=0, scale=gv%H_to_m)
999 if (
associated(cs%VarMix)) &
1000 call calc_slope_functions(h, cs%tv, dt, g, gv, us, cs%VarMix)
1002 cs%MEKE, cs%VarMix, cs%CDp, cs%thickness_diffuse_CSp)
1004 if (cs%debug)
call hchksum(h,
"Post-thickness_diffuse h", g%HI, haloshift=1, scale=gv%H_to_m)
1007 if (showcalltree)
call calltree_waypoint(
"finished thickness_diffuse (step_MOM)")
1011 if (cs%mixedlayer_restrat)
then
1013 call hchksum(h,
"Pre-mixedlayer_restrat h", g%HI, haloshift=1, scale=gv%H_to_m)
1014 call uvchksum(
"Pre-mixedlayer_restrat uhtr", &
1015 cs%uhtr, cs%vhtr, g%HI, haloshift=0, scale=gv%H_to_m*us%L_to_m**2)
1018 call mixedlayer_restrat(h, cs%uhtr, cs%vhtr, cs%tv, forces, dt, cs%visc%MLD, &
1019 cs%VarMix, g, gv, us, cs%mixedlayer_restrat_CSp)
1023 call hchksum(h,
"Post-mixedlayer_restrat h", g%HI, haloshift=1, scale=gv%H_to_m)
1024 call uvchksum(
"Post-mixedlayer_restrat [uv]htr", &
1025 cs%uhtr, cs%vhtr, g%HI, haloshift=0, scale=gv%H_to_m*us%L_to_m**2)
1031 call diag_update_remap_grids(cs%diag)
1033 if (cs%useMEKE)
call step_forward_meke(cs%MEKE, h, cs%VarMix%SN_u, cs%VarMix%SN_v, &
1034 cs%visc, dt, g, gv, us, cs%MEKE_CSp, cs%uhtr, cs%vhtr)
1035 call disable_averaging(cs%diag)
1038 cs%t_dyn_rel_adv = cs%t_dyn_rel_adv + dt
1039 cs%t_dyn_rel_thermo = cs%t_dyn_rel_thermo + dt
1040 if (abs(cs%t_dyn_rel_thermo) < 1e-6*dt) cs%t_dyn_rel_thermo = 0.0
1041 cs%t_dyn_rel_diag = cs%t_dyn_rel_diag + dt
1046 call enable_averages(dt, time_local, cs%diag)
1048 if (ids%id_u > 0)
call post_data(ids%id_u, u, cs%diag)
1049 if (ids%id_v > 0)
call post_data(ids%id_v, v, cs%diag)
1050 if (ids%id_h > 0)
call post_data(ids%id_h, h, cs%diag)
1051 call disable_averaging(cs%diag)
1064 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
1066 type(time_type),
intent(in) :: Time_local
1068 type(group_pass_type) :: pass_T_S
1069 logical :: showCallTree
1070 showcalltree = calltree_showquery()
1074 call hchksum(h,
"Pre-advection h", g%HI, haloshift=1, scale=gv%H_to_m)
1075 call uvchksum(
"Pre-advection uhtr", cs%uhtr, cs%vhtr, g%HI, &
1076 haloshift=0, scale=gv%H_to_m*us%L_to_m**2)
1077 if (
associated(cs%tv%T))
call hchksum(cs%tv%T,
"Pre-advection T", g%HI, haloshift=1)
1078 if (
associated(cs%tv%S))
call hchksum(cs%tv%S,
"Pre-advection S", g%HI, haloshift=1)
1079 if (
associated(cs%tv%frazil))
call hchksum(cs%tv%frazil, &
1080 "Pre-advection frazil", g%HI, haloshift=0)
1081 if (
associated(cs%tv%salt_deficit))
call hchksum(cs%tv%salt_deficit, &
1082 "Pre-advection salt deficit", g%HI, haloshift=0, scale=us%R_to_kg_m3*us%Z_to_m)
1088 call enable_averages(cs%t_dyn_rel_adv, time_local, cs%diag)
1090 call advect_tracer(h, cs%uhtr, cs%vhtr, cs%OBC, cs%t_dyn_rel_adv, g, gv, us, &
1091 cs%tracer_adv_CSp, cs%tracer_Reg)
1092 call tracer_hordiff(h, cs%t_dyn_rel_adv, cs%MEKE, cs%VarMix, g, gv, us, &
1093 cs%tracer_diff_CSp, cs%tracer_Reg, cs%tv)
1094 if (showcalltree)
call calltree_waypoint(
"finished tracer advection/diffusion (step_MOM)")
1096 cs%t_dyn_rel_adv, cs%tracer_Reg)
1100 call post_transport_diagnostics(g, gv, us, cs%uhtr, cs%vhtr, h, cs%transport_IDs, &
1101 cs%diag_pre_dyn, cs%diag, cs%t_dyn_rel_adv, cs%tracer_reg)
1104 call diag_update_remap_grids(cs%diag)
1106 call disable_averaging(cs%diag)
1112 cs%uhtr(:,:,:) = 0.0
1113 cs%vhtr(:,:,:) = 0.0
1114 cs%t_dyn_rel_adv = 0.0
1117 if (cs%diabatic_first .and.
associated(cs%tv%T))
then
1120 call create_group_pass(pass_t_s, cs%tv%T, g%Domain, to_all+omit_corners, halo=1)
1121 call create_group_pass(pass_t_s, cs%tv%S, g%Domain, to_all+omit_corners, halo=1)
1129 subroutine step_mom_thermo(CS, G, GV, US, u, v, h, tv, fluxes, dtdia, &
1130 Time_end_thermo, update_BBL, Waves)
1135 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
1137 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)), &
1139 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
1142 type(
forcing),
intent(inout) :: fluxes
1143 real,
intent(in) :: dtdia
1144 type(time_type),
intent(in) :: Time_end_thermo
1145 logical,
intent(in) :: update_BBL
1147 optional,
pointer :: Waves
1150 logical :: use_ice_shelf
1151 logical :: showCallTree
1152 type(group_pass_type) :: pass_T_S, pass_T_S_h, pass_uv_T_S_h
1153 integer :: dynamics_stencil
1155 integer :: i, j, k, is, ie, js, je, nz
1157 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
1158 showcalltree = calltree_showquery()
1159 if (showcalltree)
call calltree_enter(
"step_MOM_thermo(), MOM.F90")
1161 use_ice_shelf = .false.
1162 if (
associated(fluxes%frac_shelf_h)) use_ice_shelf = .true.
1164 call enable_averages(dtdia, time_end_thermo, cs%diag)
1166 if (
associated(cs%odaCS))
then
1170 if (update_bbl)
then
1176 call set_viscous_bbl(u, v, h, tv, cs%visc, g, gv, us, cs%set_visc_CSp, symmetrize=.true.)
1178 if (showcalltree)
call calltree_waypoint(
"done with set_viscous_BBL (step_MOM_thermo)")
1182 if (.not.cs%adiabatic)
then
1184 call uvchksum(
"Pre-diabatic [uv]", u, v, g%HI, haloshift=2, scale=us%L_T_to_m_s)
1185 call hchksum(h,
"Pre-diabatic h", g%HI, haloshift=1, scale=gv%H_to_m)
1186 call uvchksum(
"Pre-diabatic [uv]h", cs%uhtr, cs%vhtr, g%HI, &
1187 haloshift=0, scale=gv%H_to_m*us%L_to_m**2)
1189 call mom_thermo_chksum(
"Pre-diabatic ", tv, g, us, haloshift=0)
1196 call diabatic(u, v, h, tv, cs%Hml, fluxes, cs%visc, cs%ADp, cs%CDp, &
1197 dtdia, time_end_thermo, g, gv, us, cs%diabatic_CSp, waves=waves)
1198 fluxes%fluxes_used = .true.
1205 if ( cs%use_ALE_algorithm )
then
1206 call enable_averages(dtdia, time_end_thermo, cs%diag)
1208 if (
associated(tv%T)) &
1210 if (
associated(tv%S)) &
1213 call do_group_pass(pass_t_s_h, g%Domain)
1215 call preale_tracer_diagnostics(cs%tracer_Reg, g, gv)
1219 call hchksum(tv%T,
"Pre-ALE T", g%HI, haloshift=1)
1220 call hchksum(tv%S,
"Pre-ALE S", g%HI, haloshift=1)
1224 if (use_ice_shelf)
then
1225 call ale_main(g, gv, us, h, u, v, tv, cs%tracer_Reg, cs%ALE_CSp, cs%OBC, &
1226 dtdia, fluxes%frac_shelf_h)
1228 call ale_main(g, gv, us, h, u, v, tv, cs%tracer_Reg, cs%ALE_CSp, cs%OBC, dtdia)
1235 dynamics_stencil = min(3, g%Domain%nihalo, g%Domain%njhalo)
1237 if (
associated(tv%T)) &
1239 if (
associated(tv%S)) &
1242 call do_group_pass(pass_uv_t_s_h, g%Domain, clock=
id_clock_pass)
1244 if (cs%debug .and. cs%use_ALE_algorithm)
then
1246 call hchksum(tv%T,
"Post-ALE T", g%HI, haloshift=1)
1247 call hchksum(tv%S,
"Post-ALE S", g%HI, haloshift=1)
1254 call diag_update_remap_grids(cs%diag)
1257 call postale_tracer_diagnostics(cs%tracer_Reg, g, gv, cs%diag, dtdia)
1260 call uvchksum(
"Post-diabatic u", u, v, g%HI, haloshift=2, scale=us%L_T_to_m_s)
1261 call hchksum(h,
"Post-diabatic h", g%HI, haloshift=1, scale=gv%H_to_m)
1262 call uvchksum(
"Post-diabatic [uv]h", cs%uhtr, cs%vhtr, g%HI, &
1263 haloshift=0, scale=gv%H_to_m*us%L_to_m**2)
1266 if (
associated(tv%T))
call hchksum(tv%T,
"Post-diabatic T", g%HI, haloshift=1)
1267 if (
associated(tv%S))
call hchksum(tv%S,
"Post-diabatic S", g%HI, haloshift=1)
1268 if (
associated(tv%frazil))
call hchksum(tv%frazil, &
1269 "Post-diabatic frazil", g%HI, haloshift=0)
1270 if (
associated(tv%salt_deficit))
call hchksum(tv%salt_deficit, &
1271 "Post-diabatic salt deficit", g%HI, haloshift=0, scale=us%R_to_kg_m3*us%Z_to_m)
1275 call disable_averaging(cs%diag)
1281 call adiabatic(h, tv, fluxes, dtdia, g, gv, us, cs%diabatic_CSp)
1282 fluxes%fluxes_used = .true.
1285 if (
associated(tv%T))
then
1290 if (
associated(tv%T))
call hchksum(tv%T,
"Post-diabatic T", g%HI, haloshift=1)
1291 if (
associated(tv%S))
call hchksum(tv%S,
"Post-diabatic S", g%HI, haloshift=1)
1298 call disable_averaging(cs%diag)
1300 if (showcalltree)
call calltree_leave(
"step_MOM_thermo(), MOM.F90")
1309 subroutine step_offline(forces, fluxes, sfc_state, Time_start, time_interval, CS)
1311 type(
forcing),
intent(inout) :: fluxes
1312 type(
surface),
intent(inout) :: sfc_state
1313 type(time_type),
intent(in) :: time_start
1314 real,
intent(in) :: time_interval
1325 logical :: first_iter
1326 logical :: last_iter
1327 logical :: do_vertical
1328 logical :: adv_converged
1331 real :: dt_offline_vertical
1332 logical :: skip_diffusion
1333 integer :: id_eta_diff_end
1335 integer,
pointer :: accumulated_time => null()
1337 integer :: is, ie, js, je, isd, ied, jsd, jed
1340 real,
dimension(:,:,:),
pointer :: &
1341 uhtr => null(), vhtr => null(), &
1342 eatr => null(), ebtr => null(), &
1346 real,
dimension(SZI_(CS%G),SZJ_(CS%G)) :: eta_pre, eta_end
1347 type(time_type) :: time_end
1350 g => cs%G ; gv => cs%GV ; us => cs%US
1352 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
1353 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
1356 call extract_offline_main(cs%offline_CSp, uhtr, vhtr, eatr, ebtr, h_end, accumulated_time, &
1357 dt_offline, dt_offline_vertical, skip_diffusion)
1358 time_end = increment_date(time_start, seconds=floor(time_interval+0.001))
1360 call enable_averaging(time_interval, time_end, cs%diag)
1363 if (accumulated_time==0)
then
1366 first_iter = .false.
1370 if ( mod(accumulated_time, floor(us%T_to_s*dt_offline_vertical + 1e-6)) == 0 )
then
1371 do_vertical = .true.
1373 do_vertical = .false.
1377 accumulated_time = mod(accumulated_time + int(time_interval), floor(us%T_to_s*dt_offline+1e-6))
1378 if (accumulated_time==0)
then
1384 if (cs%use_ALE_algorithm)
then
1387 if (first_iter)
then
1388 call mom_mesg(
"Reading in new offline fields")
1393 call update_offline_fields(cs%offline_CSp, cs%h, fluxes, cs%use_ALE_algorithm)
1396 call offline_fw_fluxes_into_ocean(g, gv, cs%offline_CSp, fluxes, cs%h)
1398 if (.not.cs%diabatic_first)
then
1399 call offline_advection_ale(fluxes, time_start, time_interval, cs%offline_CSp,
id_clock_ale, &
1400 cs%h, uhtr, vhtr, converged=adv_converged)
1403 call offline_redistribute_residual(cs%offline_CSp, cs%h, uhtr, vhtr, adv_converged)
1406 if (.not. skip_diffusion)
then
1407 if (
associated(cs%VarMix))
then
1411 call calc_slope_functions(cs%h, cs%tv, dt_offline, g, gv, us, cs%VarMix)
1413 call tracer_hordiff(cs%h, dt_offline, cs%MEKE, cs%VarMix, g, gv, us, &
1414 cs%tracer_diff_CSp, cs%tracer_Reg, cs%tv)
1419 if (do_vertical)
then
1420 call offline_diabatic_ale(fluxes, time_start, time_end, cs%offline_CSp, cs%h, eatr, ebtr)
1425 if (cs%diabatic_first)
then
1426 call offline_advection_ale(fluxes, time_start, time_interval, cs%offline_CSp,
id_clock_ale, &
1427 cs%h, uhtr, vhtr, converged=adv_converged)
1430 call offline_redistribute_residual(cs%offline_CSp, cs%h, uhtr, vhtr, adv_converged)
1432 if (.not. skip_diffusion)
then
1433 if (
associated(cs%VarMix))
then
1437 call calc_slope_functions(cs%h, cs%tv, dt_offline, g, gv, us, cs%VarMix)
1439 call tracer_hordiff(cs%h, dt_offline, cs%MEKE, cs%VarMix, g, gv, us, &
1440 cs%tracer_diff_CSp, cs%tracer_Reg, cs%tv)
1444 call mom_mesg(
"Last iteration of offline interval")
1447 call offline_fw_fluxes_out_ocean(g, gv, cs%offline_CSp, fluxes, cs%h)
1450 call post_offline_convergence_diags(cs%offline_CSp, cs%h, h_end, uhtr, vhtr)
1460 call mom_error(warning, &
1461 "Offline tracer mode in non-ALE configuration has not been thoroughly tested")
1465 if (abs(time_interval - us%T_to_s*dt_offline) > 1.0e-6)
then
1466 call mom_error(fatal, &
1467 "For offline tracer mode in a non-ALE configuration, dt_offline must equal time_interval")
1469 call update_offline_fields(cs%offline_CSp, cs%h, fluxes, cs%use_ALE_algorithm)
1471 cs%h, eatr, ebtr, uhtr, vhtr)
1473 if (.not. skip_diffusion)
then
1474 call tracer_hordiff(h_end, dt_offline, cs%MEKE, cs%VarMix, g, gv, us, &
1475 cs%tracer_diff_CSp, cs%tracer_Reg, cs%tv)
1487 cs%calc_rho_for_sea_lev)
1490 call disable_averaging(cs%diag)
1495 fluxes%fluxes_used = .true.
1503 subroutine initialize_mom(Time, Time_init, param_file, dirs, CS, restart_CSp, &
1504 Time_in, offline_tracer_mode, input_restart_file, diag_ptr, &
1505 count_calls, tracer_flow_CSp)
1506 type(time_type),
target,
intent(inout) :: time
1507 type(time_type),
intent(in) :: time_init
1514 type(time_type),
optional,
intent(in) :: time_in
1516 logical,
optional,
intent(out) :: offline_tracer_mode
1517 character(len=*),
optional,
intent(in) :: input_restart_file
1518 type(
diag_ctrl),
optional,
pointer :: diag_ptr
1521 optional,
pointer :: tracer_flow_csp
1523 logical,
optional,
intent(in) :: count_calls
1531 type(
diag_ctrl),
pointer :: diag => null()
1533 character(len=4),
parameter :: vers_num =
'v2.0'
1536 # include "version_variable.h"
1538 integer :: i, j, k, is, ie, js, je, isd, ied, jsd, jed, nz
1539 integer :: isdb, iedb, jsdb, jedb
1543 real,
allocatable,
dimension(:,:) :: eta
1544 real,
allocatable,
dimension(:,:) :: area_shelf_h
1545 real,
dimension(:,:),
allocatable,
target :: frac_shelf_h
1546 real,
dimension(:,:),
pointer :: shelf_area => null()
1548 type(group_pass_type) :: tmp_pass_uv_t_s_h, pass_uv_t_s_h
1551 logical :: write_geom_files
1552 logical :: ensemble_ocean
1554 logical :: use_geothermal
1556 logical :: symmetric
1558 logical :: do_unit_tests
1559 logical :: test_grid_copy = .false.
1561 logical :: bulkmixedlayer
1563 logical :: use_temperature
1564 logical :: use_frazil
1566 logical :: bound_salinity
1568 logical :: use_cont_abss
1572 logical :: advect_ts
1574 logical :: use_ice_shelf
1575 logical :: global_indexing
1577 logical :: bathy_at_vel
1579 logical :: calc_dtbt
1581 logical :: debug_truncations
1582 integer :: first_direction
1586 integer :: nkml, nkbl, verbosity, write_geom
1587 integer :: dynamics_stencil
1589 real :: conv2watt, conv2salt
1590 character(len=48) :: flux_units, s_flux_units
1593 type(time_type) :: start_time
1595 character(len=200) :: area_varname, ice_shelf_file, inputdir, filename
1597 if (
associated(cs))
then
1598 call mom_error(warning,
"initialize_MOM called with an associated "// &
1599 "control structure.")
1604 if (test_grid_copy)
then ;
allocate(g)
1605 else ; g => cs%G ;
endif
1609 id_clock_init = cpu_clock_id(
'Ocean Initialization', grain=clock_subcomponent)
1612 start_time = time ;
if (
present(time_in)) start_time = time_in
1616 call get_mom_input(param_file, dirs, default_input_filename=input_restart_file)
1618 verbosity = 2 ;
call read_param(param_file,
"VERBOSITY", verbosity)
1619 call mom_set_verbosity(verbosity)
1626 call get_param(param_file,
"MOM",
"VERBOSITY", verbosity, &
1627 "Integer controlling level of messaging\n" // &
1628 "\t0 = Only FATAL messages\n" // &
1629 "\t2 = Only FATAL, WARNING, NOTE [default]\n" // &
1630 "\t9 = All)", default=2, debuggingparam=.true.)
1631 call get_param(param_file,
"MOM",
"DO_UNIT_TESTS", do_unit_tests, &
1632 "If True, exercises unit tests at model start up.", &
1633 default=.false., debuggingparam=.true.)
1634 if (do_unit_tests)
then
1639 call unit_scaling_init(param_file, cs%US)
1643 call get_param(param_file,
"MOM",
"SPLIT", cs%split, &
1644 "Use the split time stepping if true.", default=.true.)
1646 cs%use_RK2 = .false.
1648 call get_param(param_file,
"MOM",
"USE_RK2", cs%use_RK2, &
1649 "If true, use RK2 instead of RK3 in the unsplit time stepping.", &
1653 call get_param(param_file,
"MOM",
"CALC_RHO_FOR_SEA_LEVEL", cs%calc_rho_for_sea_lev, &
1654 "If true, the in-situ density is used to calculate the "//&
1655 "effective sea level that is returned to the coupler. If false, "//&
1656 "the Boussinesq parameter RHO_0 is used.", default=.false.)
1657 call get_param(param_file,
"MOM",
"ENABLE_THERMODYNAMICS", use_temperature, &
1658 "If true, Temperature and salinity are used as state "//&
1659 "variables.", default=.true.)
1660 call get_param(param_file,
"MOM",
"USE_EOS", use_eos, &
1661 "If true, density is calculated from temperature and "//&
1662 "salinity with an equation of state. If USE_EOS is "//&
1663 "true, ENABLE_THERMODYNAMICS must be true as well.", &
1664 default=use_temperature)
1665 call get_param(param_file,
"MOM",
"DIABATIC_FIRST", cs%diabatic_first, &
1666 "If true, apply diabatic and thermodynamic processes, "//&
1667 "including buoyancy forcing and mass gain or loss, "//&
1668 "before stepping the dynamics forward.", default=.false.)
1669 call get_param(param_file,
"MOM",
"USE_CONTEMP_ABSSAL", use_cont_abss, &
1670 "If true, the prognostics T&S are the conservative temperature "//&
1671 "and absolute salinity. Care should be taken to convert them "//&
1672 "to potential temperature and practical salinity before "//&
1673 "exchanging them with the coupler and/or reporting T&S diagnostics.", &
1675 cs%tv%T_is_conT = use_cont_abss ; cs%tv%S_is_absS = use_cont_abss
1676 call get_param(param_file,
"MOM",
"ADIABATIC", cs%adiabatic, &
1677 "There are no diapycnal mass fluxes if ADIABATIC is "//&
1678 "true. This assumes that KD = KDML = 0.0 and that "//&
1679 "there is no buoyancy forcing, but makes the model "//&
1680 "faster by eliminating subroutine calls.", default=.false.)
1681 call get_param(param_file,
"MOM",
"DO_DYNAMICS", cs%do_dynamics, &
1682 "If False, skips the dynamics calls that update u & v, as well as "//&
1683 "the gravity wave adjustment to h. This may be a fragile feature, "//&
1684 "but can be useful during development", default=.true.)
1685 call get_param(param_file,
"MOM",
"ADVECT_TS", advect_ts, &
1686 "If True, advect temperature and salinity horizontally "//&
1687 "If False, T/S are registered for advection. "//&
1688 "This is intended only to be used in offline tracer mode "//&
1689 "and is by default false in that case.", &
1690 do_not_log = .true., default=.true. )
1691 if (
present(offline_tracer_mode))
then
1692 call get_param(param_file,
"MOM",
"OFFLINE_TRACER_MODE", cs%offline_tracer_mode, &
1693 "If true, barotropic and baroclinic dynamics, thermodynamics "//&
1694 "are all bypassed with all the fields necessary to integrate "//&
1695 "the tracer advection and diffusion equation are read in from "//&
1696 "files stored from a previous integration of the prognostic model. "//&
1697 "NOTE: This option only used in the ocean_solo_driver.", default=.false.)
1698 if (cs%offline_tracer_mode)
then
1699 call get_param(param_file,
"MOM",
"ADVECT_TS", advect_ts, &
1700 "If True, advect temperature and salinity horizontally "//&
1701 "If False, T/S are registered for advection. "//&
1702 "This is intended only to be used in offline tracer mode."//&
1703 "and is by default false in that case", &
1707 call get_param(param_file,
"MOM",
"USE_REGRIDDING", cs%use_ALE_algorithm, &
1708 "If True, use the ALE algorithm (regridding/remapping). "//&
1709 "If False, use the layered isopycnal algorithm.", default=.false. )
1710 call get_param(param_file,
"MOM",
"BULKMIXEDLAYER", bulkmixedlayer, &
1711 "If true, use a Kraus-Turner-like bulk mixed layer "//&
1712 "with transitional buffer layers. Layers 1 through "//&
1713 "NKML+NKBL have variable densities. There must be at "//&
1714 "least NKML+NKBL+1 layers if BULKMIXEDLAYER is true. "//&
1715 "BULKMIXEDLAYER can not be used with USE_REGRIDDING. "//&
1716 "The default is influenced by ENABLE_THERMODYNAMICS.", &
1717 default=use_temperature .and. .not.cs%use_ALE_algorithm)
1718 call get_param(param_file,
"MOM",
"THICKNESSDIFFUSE", cs%thickness_diffuse, &
1719 "If true, interface heights are diffused with a "//&
1720 "coefficient of KHTH.", default=.false.)
1721 call get_param(param_file,
"MOM",
"THICKNESSDIFFUSE_FIRST", &
1722 cs%thickness_diffuse_first, &
1723 "If true, do thickness diffusion before dynamics. "//&
1724 "This is only used if THICKNESSDIFFUSE is true.", &
1726 if (.not.cs%thickness_diffuse) cs%thickness_diffuse_first = .false.
1727 call get_param(param_file,
"MOM",
"BATHYMETRY_AT_VEL", bathy_at_vel, &
1728 "If true, there are separate values for the basin depths "//&
1729 "at velocity points. Otherwise the effects of topography "//&
1730 "are entirely determined from thickness points.", &
1732 call get_param(param_file,
"MOM",
"USE_WAVES", cs%UseWaves, default=.false., &
1735 call get_param(param_file,
"MOM",
"DEBUG", cs%debug, &
1736 "If true, write out verbose debugging data.", &
1737 default=.false., debuggingparam=.true.)
1738 call get_param(param_file,
"MOM",
"DEBUG_TRUNCATIONS", debug_truncations, &
1739 "If true, calculate all diagnostics that are useful for "//&
1740 "debugging truncations.", default=.false., debuggingparam=.true.)
1742 call get_param(param_file,
"MOM",
"DT", cs%dt, &
1743 "The (baroclinic) dynamics time step. The time-step that "//&
1744 "is actually used will be an integer fraction of the "//&
1745 "forcing time-step (DT_FORCING in ocean-only mode or the "//&
1746 "coupling timestep in coupled mode.)", units=
"s", scale=us%s_to_T, &
1747 fail_if_missing=.true.)
1748 call get_param(param_file,
"MOM",
"DT_THERM", cs%dt_therm, &
1749 "The thermodynamic and tracer advection time step. "//&
1750 "Ideally DT_THERM should be an integer multiple of DT "//&
1751 "and less than the forcing or coupling time-step, unless "//&
1752 "THERMO_SPANS_COUPLING is true, in which case DT_THERM "//&
1753 "can be an integer multiple of the coupling timestep. By "//&
1754 "default DT_THERM is set to DT.", &
1755 units=
"s", scale=us%s_to_T, default=us%T_to_s*cs%dt)
1756 call get_param(param_file,
"MOM",
"THERMO_SPANS_COUPLING", cs%thermo_spans_coupling, &
1757 "If true, the MOM will take thermodynamic and tracer "//&
1758 "timesteps that can be longer than the coupling timestep. "//&
1759 "The actual thermodynamic timestep that is used in this "//&
1760 "case is the largest integer multiple of the coupling "//&
1761 "timestep that is less than or equal to DT_THERM.", default=.false.)
1763 if (bulkmixedlayer)
then
1764 cs%Hmix = -1.0 ; cs%Hmix_UV = -1.0
1766 call get_param(param_file,
"MOM",
"HMIX_SFC_PROP", cs%Hmix, &
1767 "If BULKMIXEDLAYER is false, HMIX_SFC_PROP is the depth "//&
1768 "over which to average to find surface properties like "//&
1769 "SST and SSS or density (but not surface velocities).", &
1770 units=
"m", default=1.0, scale=us%m_to_Z)
1771 call get_param(param_file,
"MOM",
"HMIX_UV_SFC_PROP", cs%Hmix_UV, &
1772 "If BULKMIXEDLAYER is false, HMIX_UV_SFC_PROP is the depth "//&
1773 "over which to average to find surface flow properties, "//&
1774 "SSU, SSV. A non-positive value indicates no averaging.", &
1775 units=
"m", default=0.0, scale=us%m_to_Z)
1777 call get_param(param_file,
"MOM",
"HFREEZE", cs%HFrz, &
1778 "If HFREEZE > 0, melt potential will be computed. The actual depth "//&
1779 "over which melt potential is computed will be min(HFREEZE, OBLD), "//&
1780 "where OBLD is the boundary layer depth. If HFREEZE <= 0 (default), "//&
1781 "melt potential will not be computed.", units=
"m", default=-1.0)
1782 call get_param(param_file,
"MOM",
"INTERPOLATE_P_SURF", cs%interp_p_surf, &
1783 "If true, linearly interpolate the surface pressure "//&
1784 "over the coupling time step, using the specified value "//&
1785 "at the end of the step.", default=.false.)
1788 call get_param(param_file,
"MOM",
"DTBT", dtbt, default=-0.98)
1789 default_val = us%T_to_s*cs%dt_therm ;
if (dtbt > 0.0) default_val = -1.0
1790 cs%dtbt_reset_period = -1.0
1791 call get_param(param_file,
"MOM",
"DTBT_RESET_PERIOD", cs%dtbt_reset_period, &
1792 "The period between recalculations of DTBT (if DTBT <= 0). "//&
1793 "If DTBT_RESET_PERIOD is negative, DTBT is set based "//&
1794 "only on information available at initialization. If 0, "//&
1795 "DTBT will be set every dynamics time step. The default "//&
1796 "is set by DT_THERM. This is only used if SPLIT is true.", &
1797 units=
"s", default=default_val, do_not_read=(dtbt > 0.0))
1801 use_frazil = .false. ; bound_salinity = .false. ; cs%tv%P_Ref = 2.0e7
1802 if (use_temperature)
then
1803 call get_param(param_file,
"MOM",
"FRAZIL", use_frazil, &
1804 "If true, water freezes if it gets too cold, and the "//&
1805 "the accumulated heat deficit is returned in the "//&
1806 "surface state. FRAZIL is only used if "//&
1807 "ENABLE_THERMODYNAMICS is true.", default=.false.)
1808 call get_param(param_file,
"MOM",
"DO_GEOTHERMAL", use_geothermal, &
1809 "If true, apply geothermal heating.", default=.false.)
1810 call get_param(param_file,
"MOM",
"BOUND_SALINITY", bound_salinity, &
1811 "If true, limit salinity to being positive. (The sea-ice "//&
1812 "model may ask for more salt than is available and "//&
1813 "drive the salinity negative otherwise.)", default=.false.)
1814 call get_param(param_file,
"MOM",
"MIN_SALINITY", cs%tv%min_salinity, &
1815 "The minimum value of salinity when BOUND_SALINITY=True. "//&
1816 "The default is 0.01 for backward compatibility but ideally "//&
1817 "should be 0.", units=
"PPT", default=0.01, do_not_log=.not.bound_salinity)
1818 call get_param(param_file,
"MOM",
"C_P", cs%tv%C_p, &
1819 "The heat capacity of sea water, approximated as a "//&
1820 "constant. This is only used if ENABLE_THERMODYNAMICS is "//&
1821 "true. The default value is from the TEOS-10 definition "//&
1822 "of conservative temperature.", units=
"J kg-1 K-1", &
1823 default=3991.86795711963)
1825 if (use_eos)
call get_param(param_file,
"MOM",
"P_REF", cs%tv%P_Ref, &
1826 "The pressure that is used for calculating the coordinate "//&
1827 "density. (1 Pa = 1e4 dbar, so 2e7 is commonly used.) "//&
1828 "This is only used if USE_EOS and ENABLE_THERMODYNAMICS "//&
1829 "are true.", units=
"Pa", default=2.0e7)
1831 if (bulkmixedlayer)
then
1832 call get_param(param_file,
"MOM",
"NKML", nkml, &
1833 "The number of sublayers within the mixed layer if "//&
1834 "BULKMIXEDLAYER is true.", units=
"nondim", default=2)
1835 call get_param(param_file,
"MOM",
"NKBL", nkbl, &
1836 "The number of layers that are used as variable density "//&
1837 "buffer layers if BULKMIXEDLAYER is true.", units=
"nondim", &
1841 call get_param(param_file,
"MOM",
"GLOBAL_INDEXING", global_indexing, &
1842 "If true, use a global lateral indexing convention, so "//&
1843 "that corresponding points on different processors have "//&
1844 "the same index. This does not work with static memory.", &
1845 default=.false., layoutparam=.true.)
1846 #ifdef STATIC_MEMORY_
1847 if (global_indexing)
call mom_error(fatal,
"initialize_MOM: "//&
1848 "GLOBAL_INDEXING can not be true with STATIC_MEMORY.")
1850 call get_param(param_file,
"MOM",
"FIRST_DIRECTION", first_direction, &
1851 "An integer that indicates which direction goes first "//&
1852 "in parts of the code that use directionally split "//&
1853 "updates, with even numbers (or 0) used for x- first "//&
1854 "and odd numbers used for y-first.", default=0)
1856 call get_param(param_file,
"MOM",
"CHECK_BAD_SURFACE_VALS", cs%check_bad_sfc_vals, &
1857 "If true, check the surface state for ridiculous values.", &
1859 if (cs%check_bad_sfc_vals)
then
1860 call get_param(param_file,
"MOM",
"BAD_VAL_SSH_MAX", cs%bad_val_ssh_max, &
1861 "The value of SSH above which a bad value message is "//&
1862 "triggered, if CHECK_BAD_SURFACE_VALS is true.", units=
"m", &
1864 call get_param(param_file,
"MOM",
"BAD_VAL_SSS_MAX", cs%bad_val_sss_max, &
1865 "The value of SSS above which a bad value message is "//&
1866 "triggered, if CHECK_BAD_SURFACE_VALS is true.", units=
"PPT", &
1868 call get_param(param_file,
"MOM",
"BAD_VAL_SST_MAX", cs%bad_val_sst_max, &
1869 "The value of SST above which a bad value message is "//&
1870 "triggered, if CHECK_BAD_SURFACE_VALS is true.", &
1871 units=
"deg C", default=45.0)
1872 call get_param(param_file,
"MOM",
"BAD_VAL_SST_MIN", cs%bad_val_sst_min, &
1873 "The value of SST below which a bad value message is "//&
1874 "triggered, if CHECK_BAD_SURFACE_VALS is true.", &
1875 units=
"deg C", default=-2.1)
1876 call get_param(param_file,
"MOM",
"BAD_VAL_COLUMN_THICKNESS", cs%bad_val_col_thick, &
1877 "The value of column thickness below which a bad value message is "//&
1878 "triggered, if CHECK_BAD_SURFACE_VALS is true.", units=
"m", &
1882 call get_param(param_file,
"MOM",
"SAVE_INITIAL_CONDS", save_ic, &
1883 "If true, write the initial conditions to a file given "//&
1884 "by IC_OUTPUT_FILE.", default=.false.)
1885 call get_param(param_file,
"MOM",
"IC_OUTPUT_FILE", cs%IC_file, &
1886 "The file into which to write the initial conditions.", &
1888 call get_param(param_file,
"MOM",
"WRITE_GEOM", write_geom, &
1889 "If =0, never write the geometry and vertical grid files. "//&
1890 "If =1, write the geometry and vertical grid files only for "//&
1891 "a new simulation. If =2, always write the geometry and "//&
1892 "vertical grid files. Other values are invalid.", default=1)
1893 if (write_geom<0 .or. write_geom>2)
call mom_error(fatal,
"MOM: "//&
1894 "WRITE_GEOM must be equal to 0, 1 or 2.")
1895 write_geom_files = ((write_geom==2) .or. ((write_geom==1) .and. &
1896 ((dirs%input_filename(1:1)==
'n') .and. (len_trim(dirs%input_filename)==1))))
1902 if (cs%use_ALE_algorithm .and. bulkmixedlayer)
call mom_error(fatal, &
1903 "MOM: BULKMIXEDLAYER can not currently be used with the ALE algorithm.")
1904 if (cs%use_ALE_algorithm .and. .not.use_temperature)
call mom_error(fatal, &
1905 "MOM: At this time, USE_EOS should be True when using the ALE algorithm.")
1906 if (cs%adiabatic .and. use_temperature)
call mom_error(warning, &
1907 "MOM: ADIABATIC and ENABLE_THERMODYNAMICS both defined is usually unwise.")
1908 if (use_eos .and. .not.use_temperature)
call mom_error(fatal, &
1909 "MOM: ENABLE_THERMODYNAMICS must be defined to use USE_EOS.")
1910 if (cs%adiabatic .and. bulkmixedlayer)
call mom_error(fatal, &
1911 "MOM: ADIABATIC and BULKMIXEDLAYER can not both be defined.")
1912 if (bulkmixedlayer .and. .not.use_eos)
call mom_error(fatal, &
1913 "initialize_MOM: A bulk mixed layer can only be used with T & S as "//&
1914 "state variables. Add USE_EOS = True to MOM_input.")
1916 call get_param(param_file,
'MOM',
"ICE_SHELF", use_ice_shelf, default=.false., do_not_log=.true.)
1917 if (use_ice_shelf)
then
1918 inputdir =
"." ;
call get_param(param_file,
'MOM',
"INPUTDIR", inputdir)
1919 inputdir = slasher(inputdir)
1920 call get_param(param_file,
'MOM',
"ICE_THICKNESS_FILE", ice_shelf_file, &
1921 "The file from which the ice bathymetry and area are read.", &
1922 fail_if_missing=.true.)
1923 call get_param(param_file,
'MOM',
"ICE_AREA_VARNAME", area_varname, &
1924 "The name of the area variable in ICE_THICKNESS_FILE.", &
1925 fail_if_missing=.true.)
1929 cs%ensemble_ocean=.false.
1930 call get_param(param_file,
"MOM",
"ENSEMBLE_OCEAN", cs%ensemble_ocean, &
1931 "If False, The model is being run in serial mode as a single realization. "//&
1932 "If True, The current model realization is part of a larger ensemble "//&
1933 "and at the end of step MOM, we will perform a gather of the ensemble "//&
1934 "members for statistical evaluation and/or data assimilation.", default=.false.)
1939 #ifdef SYMMETRIC_MEMORY_
1944 #ifdef STATIC_MEMORY_
1945 call mom_domains_init(g%domain, param_file, symmetric=symmetric, &
1946 static_memory=.true., nihalo=nihalo_, njhalo=njhalo_, &
1947 niglobal=niglobal_, njglobal=njglobal_, niproc=niproc_, &
1950 call mom_domains_init(g%domain, param_file, symmetric=symmetric)
1954 call mom_debugging_init(param_file)
1955 call diag_mediator_infrastructure_init()
1956 call mom_io_init(param_file)
1959 local_indexing=.not.global_indexing)
1964 call verticalgridinit( param_file, cs%GV, us )
1969 if (cs%debug .or. dg%symmetric) &
1980 if (
associated(cs%OBC))
call call_obc_register(param_file, cs%update_OBC_CSp, cs%OBC)
1982 call tracer_registry_init(param_file, cs%tracer_Reg)
1985 is = dg%isc ; ie = dg%iec ; js = dg%jsc ; je = dg%jec ; nz = gv%ke
1986 isd = dg%isd ; ied = dg%ied ; jsd = dg%jsd ; jed = dg%jed
1987 isdb = dg%IsdB ; iedb = dg%IedB ; jsdb = dg%JsdB ; jedb = dg%JedB
1988 alloc_(cs%u(isdb:iedb,jsd:jed,nz)) ; cs%u(:,:,:) = 0.0
1989 alloc_(cs%v(isd:ied,jsdb:jedb,nz)) ; cs%v(:,:,:) = 0.0
1990 alloc_(cs%h(isd:ied,jsd:jed,nz)) ; cs%h(:,:,:) = gv%Angstrom_H
1991 alloc_(cs%uh(isdb:iedb,jsd:jed,nz)) ; cs%uh(:,:,:) = 0.0
1992 alloc_(cs%vh(isd:ied,jsdb:jedb,nz)) ; cs%vh(:,:,:) = 0.0
1993 if (use_temperature)
then
1994 alloc_(cs%T(isd:ied,jsd:jed,nz)) ; cs%T(:,:,:) = 0.0
1995 alloc_(cs%S(isd:ied,jsd:jed,nz)) ; cs%S(:,:,:) = 0.0
1996 cs%tv%T => cs%T ; cs%tv%S => cs%S
1997 if (cs%tv%T_is_conT)
then
1998 vd_t = var_desc(name=
"contemp", units=
"Celsius", longname=
"Conservative Temperature", &
1999 cmor_field_name=
"thetao", cmor_longname=
"Sea Water Potential Temperature", &
2000 conversion=cs%tv%C_p)
2002 vd_t = var_desc(name=
"temp", units=
"degC", longname=
"Potential Temperature", &
2003 cmor_field_name=
"thetao", cmor_longname=
"Sea Water Potential Temperature", &
2004 conversion=cs%tv%C_p)
2006 if (cs%tv%S_is_absS)
then
2007 vd_s = var_desc(name=
"abssalt",units=
"g kg-1",longname=
"Absolute Salinity", &
2008 cmor_field_name=
"so", cmor_longname=
"Sea Water Salinity", &
2011 vd_s = var_desc(name=
"salt",units=
"psu",longname=
"Salinity", &
2012 cmor_field_name=
"so", cmor_longname=
"Sea Water Salinity", &
2018 conv2watt = gv%H_to_kg_m2 * cs%tv%C_p
2019 if (gv%Boussinesq)
then
2020 conv2salt = gv%H_to_m
2022 conv2salt = gv%H_to_kg_m2
2024 call register_tracer(cs%tv%T, cs%tracer_Reg, param_file, dg%HI, gv, &
2025 tr_desc=vd_t, registry_diags=.true., flux_nameroot=
'T', &
2026 flux_units=
'W', flux_longname=
'Heat', &
2027 flux_scale=conv2watt, convergence_units=
'W m-2', &
2028 convergence_scale=conv2watt, cmor_tendprefix=
"opottemp", diag_form=2)
2029 call register_tracer(cs%tv%S, cs%tracer_Reg, param_file, dg%HI, gv, &
2030 tr_desc=vd_s, registry_diags=.true., flux_nameroot=
'S', &
2031 flux_units=s_flux_units, flux_longname=
'Salt', &
2032 flux_scale=conv2salt, convergence_units=
'kg m-2 s-1', &
2033 convergence_scale=0.001*gv%H_to_kg_m2, cmor_tendprefix=
"osalt", diag_form=2)
2035 if (
associated(cs%OBC)) &
2036 call register_temp_salt_segments(gv, cs%OBC, cs%tracer_Reg, param_file)
2038 if (use_frazil)
then
2039 allocate(cs%tv%frazil(isd:ied,jsd:jed)) ; cs%tv%frazil(:,:) = 0.0
2041 if (bound_salinity)
then
2042 allocate(cs%tv%salt_deficit(isd:ied,jsd:jed)) ; cs%tv%salt_deficit(:,:)=0.0
2045 if (bulkmixedlayer .or. use_temperature)
then
2046 allocate(cs%Hml(isd:ied,jsd:jed)) ; cs%Hml(:,:) = 0.0
2049 if (bulkmixedlayer)
then
2050 gv%nkml = nkml ; gv%nk_rho_varies = nkml + nkbl
2052 gv%nkml = 0 ; gv%nk_rho_varies = 0
2054 if (cs%use_ALE_algorithm)
then
2055 call get_param(param_file,
"MOM",
"NK_RHO_VARIES", gv%nk_rho_varies, default=0)
2058 alloc_(cs%uhtr(isdb:iedb,jsd:jed,nz)) ; cs%uhtr(:,:,:) = 0.0
2059 alloc_(cs%vhtr(isd:ied,jsdb:jedb,nz)) ; cs%vhtr(:,:,:) = 0.0
2060 cs%t_dyn_rel_adv = 0.0 ; cs%t_dyn_rel_thermo = 0.0 ; cs%t_dyn_rel_diag = 0.0
2062 if (debug_truncations)
then
2063 allocate(cs%u_prev(isdb:iedb,jsd:jed,nz)) ; cs%u_prev(:,:,:) = 0.0
2064 allocate(cs%v_prev(isd:ied,jsdb:jedb,nz)) ; cs%v_prev(:,:,:) = 0.0
2065 mom_internal_state%u_prev => cs%u_prev
2066 mom_internal_state%v_prev => cs%v_prev
2067 call safe_alloc_ptr(cs%ADp%du_dt_visc,isdb,iedb,jsd,jed,nz)
2068 call safe_alloc_ptr(cs%ADp%dv_dt_visc,isd,ied,jsdb,jedb,nz)
2069 if (.not.cs%adiabatic)
then
2070 call safe_alloc_ptr(cs%ADp%du_dt_dia,isdb,iedb,jsd,jed,nz)
2071 call safe_alloc_ptr(cs%ADp%dv_dt_dia,isd,ied,jsdb,jedb,nz)
2075 mom_internal_state%u => cs%u ; mom_internal_state%v => cs%v
2076 mom_internal_state%h => cs%h
2077 mom_internal_state%uh => cs%uh ; mom_internal_state%vh => cs%vh
2078 if (use_temperature)
then
2079 mom_internal_state%T => cs%T ; mom_internal_state%S => cs%S
2082 cs%CDp%uh => cs%uh ; cs%CDp%vh => cs%vh
2084 if (cs%interp_p_surf)
then
2085 allocate(cs%p_surf_prev(isd:ied,jsd:jed)) ; cs%p_surf_prev(:,:) = 0.0
2088 alloc_(cs%ssh_rint(isd:ied,jsd:jed)) ; cs%ssh_rint(:,:) = 0.0
2089 alloc_(cs%ave_ssh_ibc(isd:ied,jsd:jed)) ; cs%ave_ssh_ibc(:,:) = 0.0
2090 alloc_(cs%eta_av_bc(isd:ied,jsd:jed)) ; cs%eta_av_bc(:,:) = 0.0
2091 cs%time_in_cycle = 0.0 ; cs%time_in_thermo_cycle = 0.0
2096 if (use_eos)
call eos_init(param_file, cs%tv%eqn_of_state)
2097 if (use_temperature)
then
2098 allocate(cs%tv%TempxPmE(isd:ied,jsd:jed))
2099 cs%tv%TempxPmE(:,:) = 0.0
2100 if (use_geothermal)
then
2101 allocate(cs%tv%internal_heat(isd:ied,jsd:jed))
2102 cs%tv%internal_heat(:,:) = 0.0
2112 call register_restarts_dyn_split_rk2(dg%HI, gv, param_file, &
2113 cs%dyn_split_RK2_CSp, restart_csp, cs%uh, cs%vh)
2114 elseif (cs%use_RK2)
then
2115 call register_restarts_dyn_unsplit_rk2(dg%HI, gv, param_file, &
2116 cs%dyn_unsplit_RK2_CSp, restart_csp)
2118 call register_restarts_dyn_unsplit(dg%HI, gv, param_file, &
2119 cs%dyn_unsplit_CSp, restart_csp)
2124 call call_tracer_register(dg%HI, gv, us, param_file, cs%tracer_flow_CSp, &
2125 cs%tracer_Reg, restart_csp)
2130 cs%mixedlayer_restrat_CSp, restart_csp)
2132 if (
associated(cs%OBC)) &
2133 call open_boundary_register_restarts(dg%HI, gv, cs%OBC, cs%tracer_Reg, &
2134 param_file, restart_csp, use_temperature)
2141 dirs%output_directory, cs%tv, dg%max_depth)
2144 if (cs%use_ALE_algorithm)
then
2145 call ale_init(param_file, gv, us, dg%max_depth, cs%ALE_CSp)
2153 call mom_grid_init(g, param_file, us, hi, bathymetry_at_vel=bathy_at_vel)
2160 if (cs%debug .or. g%symmetric)
then
2162 else ; g%Domain_aux => g%Domain ;
endif
2165 g%ke = gv%ke ; g%g_Earth = gv%mks_g_Earth
2168 dirs, restart_csp, cs%ALE_CSp, cs%tracer_Reg, &
2169 cs%sponge_CSp, cs%ALE_sponge_CSp, cs%OBC, time_in)
2183 if (test_grid_copy)
then
2189 call mom_grid_init(cs%G, param_file, us)
2195 call mom_grid_end(g) ;
deallocate(g)
2198 if (cs%debug .or. cs%G%symmetric)
then
2200 else ; cs%G%Domain_aux => cs%G%Domain ;
endif
2201 g%ke = gv%ke ; g%g_Earth = gv%mks_g_Earth
2208 if (ale_remap_init_conds(cs%ALE_CSp) .and. .not.
query_initialized(cs%h,
"h",restart_csp))
then
2213 call uvchksum(
"Pre ALE adjust init cond [uv]", &
2214 cs%u, cs%v, g%HI, haloshift=1)
2215 call hchksum(cs%h,
"Pre ALE adjust init cond h", g%HI, haloshift=1, scale=gv%H_to_m)
2217 call calltree_waypoint(
"Calling adjustGridForIntegrity() to remap initial conditions (initialize_MOM)")
2218 call adjustgridforintegrity(cs%ALE_CSp, g, gv, cs%h )
2219 call calltree_waypoint(
"Calling ALE_main() to remap initial conditions (initialize_MOM)")
2220 if (use_ice_shelf)
then
2221 filename = trim(inputdir)//trim(ice_shelf_file)
2222 if (.not.
file_exists(filename, g%Domain))
call mom_error(fatal, &
2223 "MOM: Unable to open "//trim(filename))
2225 allocate(area_shelf_h(isd:ied,jsd:jed))
2226 allocate(frac_shelf_h(isd:ied,jsd:jed))
2227 call mom_read_data(filename, trim(area_varname), area_shelf_h, g%Domain)
2229 frac_shelf_h(:,:) = 0.0
2231 do j=jsd,jed ;
do i=isd,ied
2232 if (g%areaT(i,j) > 0.0) &
2233 frac_shelf_h(i,j) = area_shelf_h(i,j) / (us%L_to_m**2*g%areaT(i,j))
2236 shelf_area => frac_shelf_h
2237 call ale_main(g, gv, us, cs%h, cs%u, cs%v, cs%tv, cs%tracer_Reg, cs%ALE_CSp, &
2238 cs%OBC, frac_shelf_h=shelf_area)
2240 call ale_main( g, gv, us, cs%h, cs%u, cs%v, cs%tv, cs%tracer_Reg, cs%ALE_CSp, cs%OBC)
2245 if (use_temperature)
then
2250 call do_group_pass(tmp_pass_uv_t_s_h, g%Domain)
2254 call uvchksum(
"Post ALE adjust init cond [uv]", cs%u, cs%v, g%HI, haloshift=1)
2255 call hchksum(cs%h,
"Post ALE adjust init cond h", g%HI, haloshift=1, scale=gv%H_to_m)
2258 if ( cs%use_ALE_algorithm )
call ale_updateverticalgridtype( cs%ALE_CSp, gv )
2262 call diag_mediator_init(g, gv, us, gv%ke, param_file, diag, doc_file_dir=dirs%output_directory)
2263 if (
present(diag_ptr)) diag_ptr => cs%diag
2268 call diag_masks_set(g, gv%ke, diag)
2272 call diag_set_state_ptrs(cs%h, cs%T, cs%S, cs%tv%eqn_of_state, diag)
2276 call set_axes_info(g, gv, us, param_file, diag)
2281 call diag_update_remap_grids(diag)
2284 call diag_grid_storage_init(cs%diag_pre_sync, g, diag)
2285 call diag_grid_storage_init(cs%diag_pre_dyn, g, diag)
2291 call set_masks_for_axes(g, diag)
2294 call write_static_fields(g, gv, us, cs%tv, cs%diag)
2298 call register_cell_measure(g, cs%diag, time)
2301 if (cs%use_ALE_algorithm)
then
2302 call ale_writecoordinatefile( cs%ALE_CSp, gv, dirs%output_directory )
2307 cs%useMEKE =
meke_init(time, g, us, param_file, diag, cs%MEKE_CSp, cs%MEKE, restart_csp)
2309 call varmix_init(time, g, gv, us, param_file, diag, cs%VarMix)
2310 call set_visc_init(time, g, gv, us, param_file, diag, cs%visc, cs%set_visc_CSp, restart_csp, cs%OBC)
2314 allocate(eta(szi_(g),szj_(g))) ; eta(:,:) = 0.0
2315 call initialize_dyn_split_rk2(cs%u, cs%v, cs%h, cs%uh, cs%vh, eta, time, &
2316 g, gv, us, param_file, diag, cs%dyn_split_RK2_CSp, restart_csp, &
2317 cs%dt, cs%ADp, cs%CDp, mom_internal_state, cs%VarMix, cs%MEKE, &
2318 cs%thickness_diffuse_CSp, &
2319 cs%OBC, cs%update_OBC_CSp, cs%ALE_CSp, cs%set_visc_CSp, &
2320 cs%visc, dirs, cs%ntrunc, calc_dtbt=calc_dtbt)
2321 if (cs%dtbt_reset_period > 0.0)
then
2322 cs%dtbt_reset_interval = real_to_time(cs%dtbt_reset_period)
2324 cs%dtbt_reset_time = time_init + cs%dtbt_reset_interval * &
2325 ((time - time_init) / cs%dtbt_reset_interval)
2326 if ((cs%dtbt_reset_time > time) .and. calc_dtbt)
then
2330 cs%dtbt_reset_time = cs%dtbt_reset_time - cs%dtbt_reset_interval
2333 elseif (cs%use_RK2)
then
2334 call initialize_dyn_unsplit_rk2(cs%u, cs%v, cs%h, time, g, gv, us, &
2335 param_file, diag, cs%dyn_unsplit_RK2_CSp, restart_csp, &
2336 cs%ADp, cs%CDp, mom_internal_state, cs%MEKE, cs%OBC, &
2337 cs%update_OBC_CSp, cs%ALE_CSp, cs%set_visc_CSp, cs%visc, dirs, &
2340 call initialize_dyn_unsplit(cs%u, cs%v, cs%h, time, g, gv, us, &
2341 param_file, diag, cs%dyn_unsplit_CSp, restart_csp, &
2342 cs%ADp, cs%CDp, mom_internal_state, cs%MEKE, cs%OBC, &
2343 cs%update_OBC_CSp, cs%ALE_CSp, cs%set_visc_CSp, cs%visc, dirs, &
2349 cs%mixedlayer_restrat = mixedlayer_restrat_init(time, g, gv, us, param_file, diag, &
2350 cs%mixedlayer_restrat_CSp, restart_csp)
2351 if (cs%mixedlayer_restrat)
then
2352 if (.not.(bulkmixedlayer .or. cs%use_ALE_algorithm)) &
2353 call mom_error(fatal,
"MOM: MIXEDLAYER_RESTRAT true requires a boundary layer scheme.")
2355 if (.not. cs%diabatic_first .and.
associated(cs%visc%MLD)) &
2356 call pass_var(cs%visc%MLD, g%domain, halo=1)
2359 call mom_diagnostics_init(mom_internal_state, cs%ADp, cs%CDp, time, g, gv, us, &
2360 param_file, diag, cs%diagnostics_CSp, cs%tv)
2363 if (
associated(cs%sponge_CSp)) &
2366 if (
associated(cs%ALE_sponge_CSp)) &
2369 if (cs%adiabatic)
then
2373 call diabatic_driver_init(time, g, gv, us, param_file, cs%use_ALE_algorithm, diag, &
2374 cs%ADp, cs%CDp, cs%diabatic_CSp, cs%tracer_flow_CSp, &
2375 cs%sponge_CSp, cs%ALE_sponge_CSp)
2378 call tracer_advect_init(time, g, us, param_file, diag, cs%tracer_adv_CSp)
2379 call tracer_hor_diff_init(time, g, us, param_file, diag, cs%tv%eqn_of_state, cs%diabatic_CSp, &
2386 call register_surface_diags(time, g, us, cs%sfc_IDs, cs%diag, cs%tv)
2388 call register_transport_diags(time, g, gv, us, cs%transport_IDs, cs%diag)
2389 call register_tracer_diagnostics(cs%tracer_Reg, cs%h, time, diag, g, gv, us, &
2390 cs%use_ALE_algorithm)
2391 if (cs%use_ALE_algorithm)
then
2392 call ale_register_diags(time, g, gv, us, diag, cs%ALE_CSp)
2397 call tracer_flow_control_init(.not.new_sim, time, g, gv, us, cs%h, param_file, &
2398 cs%diag, cs%OBC, cs%tracer_flow_CSp, cs%sponge_CSp, &
2399 cs%ALE_sponge_CSp, cs%tv)
2400 if (
present(tracer_flow_csp)) tracer_flow_csp => cs%tracer_flow_CSp
2404 if (
present(offline_tracer_mode)) offline_tracer_mode=cs%offline_tracer_mode
2406 if (cs%offline_tracer_mode)
then
2408 call offline_transport_init(param_file, cs%offline_CSp, cs%diabatic_CSp, g, gv, us)
2409 call insert_offline_main( cs=cs%offline_CSp, ale_csp=cs%ALE_CSp, diabatic_csp=cs%diabatic_CSp, &
2410 diag=cs%diag, obc=cs%OBC, tracer_adv_csp=cs%tracer_adv_CSp, &
2411 tracer_flow_csp=cs%tracer_flow_CSp, tracer_reg=cs%tracer_Reg, &
2412 tv=cs%tv, x_before_y = (mod(first_direction,2)==0), debug=cs%debug )
2413 call register_diags_offline_transport(time, cs%diag, cs%offline_CSp)
2418 dynamics_stencil = min(3, g%Domain%nihalo, g%Domain%njhalo)
2419 call create_group_pass(pass_uv_t_s_h, cs%u, cs%v, g%Domain, halo=dynamics_stencil)
2420 if (use_temperature)
then
2426 call do_group_pass(pass_uv_t_s_h, g%Domain)
2428 if (
associated(cs%visc%Kv_shear)) &
2429 call pass_var(cs%visc%Kv_shear, g%Domain, to_all+omit_corners, halo=1)
2431 if (
associated(cs%visc%Kv_slow)) &
2432 call pass_var(cs%visc%Kv_slow, g%Domain, to_all+omit_corners, halo=1)
2438 if (use_frazil)
then
2440 cs%tv%frazil(:,:) = 0.0
2443 if (cs%interp_p_surf)
then
2444 cs%p_surf_prev_set = &
2447 if (cs%p_surf_prev_set)
call pass_var(cs%p_surf_prev, g%domain)
2452 call find_eta(cs%h, cs%tv, g, gv, us, cs%ave_ssh_ibc, eta, eta_to_m=1.0)
2454 call find_eta(cs%h, cs%tv, g, gv, us, cs%ave_ssh_ibc, eta_to_m=1.0)
2457 if (cs%split)
deallocate(eta)
2460 if (
present(count_calls)) cs%count_calls = count_calls
2462 cs%ntrunc, time_init, cs%sum_output_CSp)
2465 cs%write_IC = save_ic .and. &
2466 .not.((dirs%input_filename(1:1) ==
'r') .and. &
2467 (len_trim(dirs%input_filename) == 1))
2469 if (cs%ensemble_ocean)
then
2470 call init_oda(time, g, gv, cs%odaCS)
2484 type(time_type),
intent(in) :: time
2496 real,
allocatable :: z_interface(:,:,:)
2503 g => cs%G ; gv => cs%GV ; us => cs%US
2506 call fix_restart_scaling(gv)
2510 if (cs%write_IC)
then
2511 allocate(restart_csp_tmp)
2512 restart_csp_tmp = restart_csp
2513 allocate(z_interface(szi_(g),szj_(g),szk_(g)+1))
2514 call find_eta(cs%h, cs%tv, g, gv, us, z_interface, eta_to_m=1.0)
2516 "Interface heights",
"meter", z_grid=
'i')
2518 call save_restart(dirs%output_directory, time, g, &
2519 restart_csp_tmp, filename=cs%IC_file, gv=gv)
2520 deallocate(z_interface)
2521 deallocate(restart_csp_tmp)
2524 call write_energy(cs%u, cs%v, cs%h, cs%tv, time, 0, g, gv, us, &
2525 cs%sum_output_CSp, cs%tracer_flow_CSp)
2534 type(time_type),
intent(in) :: Time
2542 character(len=48) :: thickness_units
2545 if (gv%Boussinesq)
then
2546 h_convert = gv%H_to_m
2548 h_convert = gv%H_to_kg_m2
2552 ids%id_u = register_diag_field(
'ocean_model',
'u_dyn', diag%axesCuL, time, &
2553 'Zonal velocity after the dynamics update',
'm s-1', conversion=us%L_T_to_m_s)
2554 ids%id_v = register_diag_field(
'ocean_model',
'v_dyn', diag%axesCvL, time, &
2555 'Meridional velocity after the dynamics update',
'm s-1', conversion=us%L_T_to_m_s)
2556 ids%id_h = register_diag_field(
'ocean_model',
'h_dyn', diag%axesTL, time, &
2557 'Layer Thickness after the dynamics update', thickness_units, &
2558 v_extensive=.true., conversion=h_convert)
2559 ids%id_ssh_inst = register_diag_field(
'ocean_model',
'SSH_inst', diag%axesT1, &
2560 time,
'Instantaneous Sea Surface Height',
'm')
2569 id_clock_thermo = cpu_clock_id(
'Ocean thermodynamics and tracers', grain=clock_subcomponent)
2570 id_clock_other = cpu_clock_id(
'Ocean Other', grain=clock_subcomponent)
2571 id_clock_tracer = cpu_clock_id(
'(Ocean tracer advection)', grain=clock_module_driver)
2572 if (.not.cs%adiabatic)
then
2573 id_clock_diabatic = cpu_clock_id(
'(Ocean diabatic driver)', grain=clock_module_driver)
2575 id_clock_adiabatic = cpu_clock_id(
'(Ocean adiabatic driver)', grain=clock_module_driver)
2579 id_clock_bbl_visc = cpu_clock_id(
'(Ocean set BBL viscosity)', grain=clock_module)
2580 id_clock_pass = cpu_clock_id(
'(Ocean message passing *)', grain=clock_module)
2581 id_clock_mom_init = cpu_clock_id(
'(Ocean MOM_initialize_state)', grain=clock_module)
2582 id_clock_pass_init = cpu_clock_id(
'(Ocean init message passing *)', grain=clock_routine)
2583 if (cs%thickness_diffuse) &
2588 id_clock_z_diag = cpu_clock_id(
'(Ocean Z-space diagnostics)', grain=clock_module)
2589 id_clock_ale = cpu_clock_id(
'(Ocean ALE)', grain=clock_module)
2590 if (cs%offline_tracer_mode)
then
2613 logical :: use_ice_shelf
2614 character(len=48) :: thickness_units, flux_units
2620 if (
associated(cs%tv%T)) &
2622 "Potential Temperature",
"degC")
2623 if (
associated(cs%tv%S)) &
2628 "Layer Thickness", thickness_units)
2631 "Zonal velocity",
"m s-1", hor_grid=
'Cu')
2634 "Meridional velocity",
"m s-1", hor_grid=
'Cv')
2636 if (
associated(cs%tv%frazil)) &
2638 "Frazil heat flux into ocean",
"J m-2")
2640 if (cs%interp_p_surf)
then
2642 "Previous ocean surface pressure",
"Pa")
2646 "Time average sea surface height",
"meter")
2649 call get_param(param_file,
'',
"ICE_SHELF", use_ice_shelf, default=.false., &
2651 if (use_ice_shelf .and.
associated(cs%Hml))
then
2653 "Mixed layer thickness",
"meter")
2658 "Height unit conversion factor",
"Z meter-1")
2660 "Thickness unit conversion factor",
"H meter-1")
2662 "Length unit conversion factor",
"L meter-1")
2664 "Time unit conversion factor",
"T second-1")
2666 "Density unit conversion factor",
"R m3 kg-1")
2677 real,
dimension(SZI_(G),SZJ_(G)),
intent(inout) :: ssh
2678 real,
dimension(:,:),
optional,
pointer :: p_atm
2679 logical,
optional,
intent(in) :: use_EOS
2686 integer :: i, j, is, ie, js, je
2688 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
2689 if (
present(p_atm))
then ;
if (
associated(p_atm))
then
2690 calc_rho =
associated(tv%eqn_of_state)
2691 if (
present(use_eos) .and. calc_rho) calc_rho = use_eos
2694 do j=js,je ;
do i=is,ie
2697 rho_conv, tv%eqn_of_state, scale=us%kg_m3_to_R)
2701 igr0 = 1.0 / (rho_conv * us%R_to_kg_m3*gv%mks_g_Earth)
2702 ssh(i,j) = ssh(i,j) + p_atm(i,j) * igr0
2713 type(
surface),
intent(inout) :: sfc_state
2723 real,
dimension(:,:,:),
pointer :: &
2725 real :: depth(szi_(cs%g))
2732 real :: delt(szi_(cs%g))
2733 logical :: use_temperature
2734 integer :: i, j, k, is, ie, js, je, nz, numberoferrors, ig, jg
2735 integer :: isd, ied, jsd, jed
2736 integer :: iscb, iecb, jscb, jecb, isdb, iedb, jsdb, jedb
2737 logical :: localerror
2738 character(240) :: msg
2741 g => cs%G ; gv => cs%GV ; us => cs%US
2742 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
2743 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
2744 iscb = g%iscB ; iecb = g%iecB; jscb = g%jscB ; jecb = g%jecB
2745 isdb = g%isdB ; iedb = g%iedB; jsdb = g%jsdB ; jedb = g%jedB
2748 use_temperature =
associated(cs%tv%T)
2750 if (.not.sfc_state%arrays_allocated)
then
2753 call allocate_surface_state(sfc_state, g, use_temperature, do_integrals=.true.)
2755 sfc_state%frazil => cs%tv%frazil
2756 sfc_state%T_is_conT = cs%tv%T_is_conT
2757 sfc_state%S_is_absS = cs%tv%S_is_absS
2759 do j=js,je ;
do i=is,ie
2760 sfc_state%sea_lev(i,j) = cs%ave_ssh_ibc(i,j)
2764 if (
associated(cs%Hml))
then
2765 do j=js,je ;
do i=is,ie
2766 sfc_state%Hml(i,j) = cs%Hml(i,j)
2770 if (cs%Hmix < 0.0)
then
2771 if (use_temperature)
then ;
do j=js,je ;
do i=is,ie
2772 sfc_state%SST(i,j) = cs%tv%T(i,j,1)
2773 sfc_state%SSS(i,j) = cs%tv%S(i,j,1)
2774 enddo ;
enddo ;
endif
2775 do j=js,je ;
do i=is-1,ie
2776 sfc_state%u(i,j) = us%L_T_to_m_s * cs%u(i,j,1)
2778 do j=js-1,je ;
do i=is,ie
2779 sfc_state%v(i,j) = us%L_T_to_m_s * cs%v(i,j,1)
2791 if (use_temperature)
then
2792 sfc_state%SST(i,j) = 0.0 ; sfc_state%SSS(i,j) = 0.0
2794 sfc_state%sfc_density(i,j) = 0.0
2798 do k=1,nz ;
do i=is,ie
2799 if (depth(i) + h(i,j,k)*gv%H_to_Z < depth_ml)
then
2800 dh = h(i,j,k)*gv%H_to_Z
2801 elseif (depth(i) < depth_ml)
then
2802 dh = depth_ml - depth(i)
2806 if (use_temperature)
then
2807 sfc_state%SST(i,j) = sfc_state%SST(i,j) + dh * cs%tv%T(i,j,k)
2808 sfc_state%SSS(i,j) = sfc_state%SSS(i,j) + dh * cs%tv%S(i,j,k)
2810 sfc_state%sfc_density(i,j) = sfc_state%sfc_density(i,j) + dh * us%R_to_kg_m3*gv%Rlay(k)
2812 depth(i) = depth(i) + dh
2816 if (depth(i) < gv%H_subroundoff*gv%H_to_Z) &
2817 depth(i) = gv%H_subroundoff*gv%H_to_Z
2818 if (use_temperature)
then
2819 sfc_state%SST(i,j) = sfc_state%SST(i,j) / depth(i)
2820 sfc_state%SSS(i,j) = sfc_state%SSS(i,j) / depth(i)
2822 sfc_state%sfc_density(i,j) = sfc_state%sfc_density(i,j) / depth(i)
2833 if (cs%Hmix_UV>0.)
then
2836 depth_ml = cs%Hmix_UV
2841 sfc_state%v(i,j) = 0.0
2843 do k=1,nz ;
do i=is,ie
2844 hv = 0.5 * (h(i,j,k) + h(i,j+1,k)) * gv%H_to_Z
2845 if (depth(i) + hv < depth_ml)
then
2847 elseif (depth(i) < depth_ml)
then
2848 dh = depth_ml - depth(i)
2852 sfc_state%v(i,j) = sfc_state%v(i,j) + dh * us%L_T_to_m_s * cs%v(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%v(i,j) = sfc_state%v(i,j) / depth(i)
2867 sfc_state%u(i,j) = 0.0
2869 do k=1,nz ;
do i=is-1,ie
2870 hu = 0.5 * (h(i,j,k) + h(i+1,j,k)) * gv%H_to_Z
2871 if (depth(i) + hu < depth_ml)
then
2873 elseif (depth(i) < depth_ml)
then
2874 dh = depth_ml - depth(i)
2878 sfc_state%u(i,j) = sfc_state%u(i,j) + dh * us%L_T_to_m_s * cs%u(i,j,k)
2879 depth(i) = depth(i) + dh
2883 if (depth(i) < gv%H_subroundoff*gv%H_to_Z) &
2884 depth(i) = gv%H_subroundoff*gv%H_to_Z
2885 sfc_state%u(i,j) = sfc_state%u(i,j) / depth(i)
2889 do j=js,je ;
do i=is-1,ie
2890 sfc_state%u(i,j) = us%L_T_to_m_s * cs%u(i,j,1)
2892 do j=js-1,je ;
do i=is,ie
2893 sfc_state%v(i,j) = us%L_T_to_m_s * cs%v(i,j,1)
2899 if (
allocated(sfc_state%melt_potential))
then
2907 do k=1,nz ;
do i=is,ie
2908 depth_ml = min(cs%HFrz,cs%visc%MLD(i,j))
2909 if (depth(i) + h(i,j,k)*gv%H_to_m < depth_ml)
then
2910 dh = h(i,j,k)*gv%H_to_m
2911 elseif (depth(i) < depth_ml)
then
2912 dh = depth_ml - depth(i)
2919 depth(i) = depth(i) + dh
2920 delt(i) = delt(i) + dh * (cs%tv%T(i,j,k) - t_freeze)
2925 sfc_state%melt_potential(i,j) = 0.0
2927 if (g%mask2dT(i,j)>0.)
then
2929 sfc_state%melt_potential(i,j) = cs%tv%C_p * us%R_to_kg_m3*gv%Rho0 * delt(i)
2935 if (
allocated(sfc_state%salt_deficit) .and.
associated(cs%tv%salt_deficit))
then
2937 do j=js,je ;
do i=is,ie
2939 sfc_state%salt_deficit(i,j) = 1000.0 * us%R_to_kg_m3*us%Z_to_m*cs%tv%salt_deficit(i,j)
2942 if (
allocated(sfc_state%TempxPmE) .and.
associated(cs%tv%TempxPmE))
then
2944 do j=js,je ;
do i=is,ie
2945 sfc_state%TempxPmE(i,j) = us%R_to_kg_m3*us%Z_to_m*cs%tv%TempxPmE(i,j)
2948 if (
allocated(sfc_state%internal_heat) .and.
associated(cs%tv%internal_heat))
then
2950 do j=js,je ;
do i=is,ie
2951 sfc_state%internal_heat(i,j) = cs%tv%internal_heat(i,j)
2954 if (
allocated(sfc_state%taux_shelf) .and.
associated(cs%visc%taux_shelf))
then
2956 do j=js,je ;
do i=is-1,ie
2957 sfc_state%taux_shelf(i,j) = us%R_to_kg_m3*us%L_T_to_m_s**2*us%Z_to_L*cs%visc%taux_shelf(i,j)
2960 if (
allocated(sfc_state%tauy_shelf) .and.
associated(cs%visc%tauy_shelf))
then
2962 do j=js-1,je ;
do i=is,ie
2963 sfc_state%tauy_shelf(i,j) = us%R_to_kg_m3*us%L_T_to_m_s**2*us%Z_to_L*cs%visc%tauy_shelf(i,j)
2967 if (
allocated(sfc_state%ocean_mass) .and.
allocated(sfc_state%ocean_heat) .and. &
2968 allocated(sfc_state%ocean_salt))
then
2970 do j=js,je ;
do i=is,ie
2971 sfc_state%ocean_mass(i,j) = 0.0
2972 sfc_state%ocean_heat(i,j) = 0.0 ; sfc_state%ocean_salt(i,j) = 0.0
2975 do j=js,je ;
do k=1,nz;
do i=is,ie
2976 mass = gv%H_to_kg_m2*h(i,j,k)
2977 sfc_state%ocean_mass(i,j) = sfc_state%ocean_mass(i,j) + mass
2978 sfc_state%ocean_heat(i,j) = sfc_state%ocean_heat(i,j) + mass*cs%tv%T(i,j,k)
2979 sfc_state%ocean_salt(i,j) = sfc_state%ocean_salt(i,j) + &
2980 mass * (1.0e-3*cs%tv%S(i,j,k))
2981 enddo ;
enddo ;
enddo
2983 if (
allocated(sfc_state%ocean_mass))
then
2985 do j=js,je ;
do i=is,ie ; sfc_state%ocean_mass(i,j) = 0.0 ;
enddo ;
enddo
2987 do j=js,je ;
do k=1,nz ;
do i=is,ie
2988 sfc_state%ocean_mass(i,j) = sfc_state%ocean_mass(i,j) + gv%H_to_kg_m2*h(i,j,k)
2989 enddo ;
enddo ;
enddo
2991 if (
allocated(sfc_state%ocean_heat))
then
2993 do j=js,je ;
do i=is,ie ; sfc_state%ocean_heat(i,j) = 0.0 ;
enddo ;
enddo
2995 do j=js,je ;
do k=1,nz ;
do i=is,ie
2996 mass = gv%H_to_kg_m2*h(i,j,k)
2997 sfc_state%ocean_heat(i,j) = sfc_state%ocean_heat(i,j) + mass*cs%tv%T(i,j,k)
2998 enddo ;
enddo ;
enddo
3000 if (
allocated(sfc_state%ocean_salt))
then
3002 do j=js,je ;
do i=is,ie ; sfc_state%ocean_salt(i,j) = 0.0 ;
enddo ;
enddo
3004 do j=js,je ;
do k=1,nz ;
do i=is,ie
3005 mass = gv%H_to_kg_m2*h(i,j,k)
3006 sfc_state%ocean_salt(i,j) = sfc_state%ocean_salt(i,j) + &
3007 mass * (1.0e-3*cs%tv%S(i,j,k))
3008 enddo ;
enddo ;
enddo
3012 if (
associated(cs%tracer_flow_CSp))
then
3013 call call_tracer_surface_state(sfc_state, h, g, cs%tracer_flow_CSp)
3016 if (cs%check_bad_sfc_vals)
then
3018 do j=js,je;
do i=is,ie
3019 if (g%mask2dT(i,j)>0.)
then
3020 bathy_m = cs%US%Z_to_m * g%bathyT(i,j)
3021 localerror = sfc_state%sea_lev(i,j)<=-bathy_m &
3022 .or. sfc_state%sea_lev(i,j)>= cs%bad_val_ssh_max &
3023 .or. sfc_state%sea_lev(i,j)<=-cs%bad_val_ssh_max &
3024 .or. sfc_state%sea_lev(i,j) + bathy_m < cs%bad_val_col_thick
3025 if (use_temperature) localerror = localerror &
3026 .or. sfc_state%SSS(i,j)<0. &
3027 .or. sfc_state%SSS(i,j)>=cs%bad_val_sss_max &
3028 .or. sfc_state%SST(i,j)< cs%bad_val_sst_min &
3029 .or. sfc_state%SST(i,j)>=cs%bad_val_sst_max
3030 if (localerror)
then
3031 numberoferrors=numberoferrors+1
3032 if (numberoferrors<9)
then
3033 ig = i + g%HI%idg_offset
3034 jg = j + g%HI%jdg_offset
3035 if (use_temperature)
then
3036 write(msg(1:240),
'(2(a,i4,x),4(a,f8.3,x),8(a,es11.4,x))') &
3037 'Extreme surface sfc_state detected: i=',ig,
'j=',jg, &
3038 'lon=',g%geoLonT(i,j),
'lat=',g%geoLatT(i,j), &
3039 'x=',g%gridLonT(ig),
'y=',g%gridLatT(jg), &
3040 'D=',bathy_m,
'SSH=',sfc_state%sea_lev(i,j), &
3041 'SST=',sfc_state%SST(i,j),
'SSS=',sfc_state%SSS(i,j), &
3042 'U-=',sfc_state%u(i-1,j),
'U+=',sfc_state%u(i,j), &
3043 'V-=',sfc_state%v(i,j-1),
'V+=',sfc_state%v(i,j)
3045 write(msg(1:240),
'(2(a,i4,x),4(a,f8.3,x),6(a,es11.4))') &
3046 'Extreme surface sfc_state detected: i=',ig,
'j=',jg, &
3047 'lon=',g%geoLonT(i,j),
'lat=',g%geoLatT(i,j), &
3048 'x=',g%gridLonT(i),
'y=',g%gridLatT(j), &
3049 'D=',bathy_m,
'SSH=',sfc_state%sea_lev(i,j), &
3050 'U-=',sfc_state%u(i-1,j),
'U+=',sfc_state%u(i,j), &
3051 'V-=',sfc_state%v(i,j-1),
'V+=',sfc_state%v(i,j)
3053 call mom_error(warning, trim(msg), all_print=.true.)
3054 elseif (numberoferrors==9)
then
3055 call mom_error(warning,
'There were more unreported extreme events!', all_print=.true.)
3060 call sum_across_pes(numberoferrors)
3061 if (numberoferrors>0)
then
3062 write(msg(1:240),
'(3(a,i9,x))')
'There were a total of ',numberoferrors, &
3063 'locations detected with extreme surface values!'
3064 call mom_error(fatal, trim(msg))
3076 logical,
optional,
intent(in) :: adv_dyn
3083 adv_only = .false. ;
if (
present(adv_dyn)) adv_only = adv_dyn
3086 in_synch = (cs%t_dyn_rel_adv == 0.0)
3088 in_synch = ((cs%t_dyn_rel_adv == 0.0) .and. (cs%t_dyn_rel_thermo == 0.0))
3098 optional,
pointer :: g
3100 optional,
pointer :: gv
3102 optional,
pointer :: us
3103 real,
optional,
intent(out) :: c_p
3104 logical,
optional,
intent(out) :: use_temp
3106 if (
present(g)) g => cs%G
3107 if (
present(gv)) gv => cs%GV
3108 if (
present(us)) us => cs%US
3109 if (
present(c_p)) c_p = cs%tv%C_p
3110 if (
present(use_temp)) use_temp =
associated(cs%tv%T)
3116 real,
optional,
intent(out) :: heat
3117 real,
optional,
intent(out) :: salt
3118 real,
optional,
intent(out) :: mass
3119 logical,
optional,
intent(in) :: on_pe_only
3121 if (
present(mass)) &
3123 if (
present(heat)) &
3125 if (
present(salt)) &
3134 if (cs%use_ALE_algorithm)
call ale_end(cs%ALE_CSp)
3136 dealloc_(cs%u) ; dealloc_(cs%v) ; dealloc_(cs%h)
3137 dealloc_(cs%uh) ; dealloc_(cs%vh)
3139 if (
associated(cs%tv%T))
then
3140 dealloc_(cs%T) ; cs%tv%T => null() ; dealloc_(cs%S) ; cs%tv%S => null()
3142 if (
associated(cs%tv%frazil))
deallocate(cs%tv%frazil)
3143 if (
associated(cs%tv%salt_deficit))
deallocate(cs%tv%salt_deficit)
3144 if (
associated(cs%Hml))
deallocate(cs%Hml)
3155 dealloc_(cs%uhtr) ; dealloc_(cs%vhtr)
3157 call end_dyn_split_rk2(cs%dyn_split_RK2_CSp)
3158 elseif (cs%use_RK2)
then
3159 call end_dyn_unsplit_rk2(cs%dyn_unsplit_RK2_CSp)
3161 call end_dyn_unsplit(cs%dyn_unsplit_CSp)
3163 dealloc_(cs%ave_ssh_ibc) ; dealloc_(cs%ssh_rint)
3164 if (
associated(cs%update_OBC_CSp))
call obc_register_end(cs%update_OBC_CSp)
3166 call verticalgridend(cs%GV)
3168 call mom_grid_end(cs%G)