From 9b7b127bafdd35a95510348c71ba94d25fe717e5 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Wed, 27 Feb 2019 18:22:31 -0500 Subject: [PATCH 1/5] +New variants of safe_alloc_ptr & safe_alloc_alloc Added new overloaded variants of safe_alloc_ptr and safe_alloc_alloc to allow the index range of the third (vertical) index to have an arbitrary starting value. All answers are bitwise identical. --- src/framework/MOM_safe_alloc.F90 | 39 +++++++++++++++++++++++++++++--- 1 file changed, 36 insertions(+), 3 deletions(-) diff --git a/src/framework/MOM_safe_alloc.F90 b/src/framework/MOM_safe_alloc.F90 index 75f5fda74e..47dd8376a3 100644 --- a/src/framework/MOM_safe_alloc.F90 +++ b/src/framework/MOM_safe_alloc.F90 @@ -10,13 +10,14 @@ module MOM_safe_alloc !> Allocate a pointer to a 1-d, 2-d or 3-d array interface safe_alloc_ptr - module procedure safe_alloc_ptr_3d_2arg, safe_alloc_ptr_2d_2arg + module procedure safe_alloc_ptr_3d_3arg, safe_alloc_ptr_3d_6arg, safe_alloc_ptr_2d_2arg module procedure safe_alloc_ptr_3d, safe_alloc_ptr_2d, safe_alloc_ptr_1d end interface safe_alloc_ptr !> Allocate a 2-d or 3-d allocatable array interface safe_alloc_alloc module procedure safe_alloc_allocatable_3d, safe_alloc_allocatable_2d + module procedure safe_alloc_allocatable_3d_6arg end interface safe_alloc_alloc ! This combined interface might work with a later version of Fortran, but @@ -57,7 +58,7 @@ subroutine safe_alloc_ptr_2d_2arg(ptr, ni, nj) end subroutine safe_alloc_ptr_2d_2arg !> Allocate a pointer to a 3-d array based on its dimension sizes -subroutine safe_alloc_ptr_3d_2arg(ptr, ni, nj, nk) +subroutine safe_alloc_ptr_3d_3arg(ptr, ni, nj, nk) real, dimension(:,:,:), pointer :: ptr !< A pointer to allocate integer, intent(in) :: ni !< The size of the 1st dimension of the array integer, intent(in) :: nj !< The size of the 2nd dimension of the array @@ -66,7 +67,7 @@ subroutine safe_alloc_ptr_3d_2arg(ptr, ni, nj, nk) allocate(ptr(ni,nj,nk)) ptr(:,:,:) = 0.0 endif -end subroutine safe_alloc_ptr_3d_2arg +end subroutine safe_alloc_ptr_3d_3arg !> Allocate a pointer to a 2-d array based on its index starting and ending values subroutine safe_alloc_ptr_2d(ptr, is, ie, js, je) @@ -95,6 +96,22 @@ subroutine safe_alloc_ptr_3d(ptr, is, ie, js, je, nk) endif end subroutine safe_alloc_ptr_3d +!> Allocate a pointer to a 3-d array based on its index starting and ending values +subroutine safe_alloc_ptr_3d_6arg(ptr, is, ie, js, je, ks, ke) + real, dimension(:,:,:), pointer :: ptr !< A pointer to allocate + integer, intent(in) :: is !< The start index to allocate for the 1st dimension + integer, intent(in) :: ie !< The end index to allocate for the 1st dimension + integer, intent(in) :: js !< The start index to allocate for the 2nd dimension + integer, intent(in) :: je !< The end index to allocate for the 2nd dimension + integer, intent(in) :: ks !< The start index to allocate for the 3rd dimension + integer, intent(in) :: ke !< The end index to allocate for the 3rd dimension + if (.not.associated(ptr)) then + allocate(ptr(is:ie,js:je,ks:ke)) + ptr(:,:,:) = 0.0 + endif +end subroutine safe_alloc_ptr_3d_6arg + + !> Allocate a 2-d allocatable array based on its index starting and ending values subroutine safe_alloc_allocatable_2d(ptr, is, ie, js, je) real, dimension(:,:), allocatable :: ptr !< An allocatable array to allocate @@ -109,6 +126,7 @@ subroutine safe_alloc_allocatable_2d(ptr, is, ie, js, je) end subroutine safe_alloc_allocatable_2d !> Allocate a 3-d allocatable array based on its index starting and ending values +!! and k-index size subroutine safe_alloc_allocatable_3d(ptr, is, ie, js, je, nk) real, dimension(:,:,:), allocatable :: ptr !< An allocatable array to allocate integer, intent(in) :: is !< The start index to allocate for the 1st dimension @@ -122,4 +140,19 @@ subroutine safe_alloc_allocatable_3d(ptr, is, ie, js, je, nk) endif end subroutine safe_alloc_allocatable_3d +!> Allocate a 3-d allocatable array based on its 6 index starting and ending values +subroutine safe_alloc_allocatable_3d_6arg(ptr, is, ie, js, je, ks, ke) + real, dimension(:,:,:), allocatable :: ptr !< An allocatable array to allocate + integer, intent(in) :: is !< The start index to allocate for the 1st dimension + integer, intent(in) :: ie !< The end index to allocate for the 1st dimension + integer, intent(in) :: js !< The start index to allocate for the 2nd dimension + integer, intent(in) :: je !< The end index to allocate for the 2nd dimension + integer, intent(in) :: ks !< The start index to allocate for the 3rd dimension + integer, intent(in) :: ke !< The end index to allocate for the 3rd dimension + if (.not.allocated(ptr)) then + allocate(ptr(is:ie,js:je,ks:ke)) + ptr(:,:,:) = 0.0 + endif +end subroutine safe_alloc_allocatable_3d_6arg + end module MOM_safe_alloc From 0ff4c5b195ea9599aea322053edf36f57151ec27 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Wed, 27 Feb 2019 18:22:56 -0500 Subject: [PATCH 2/5] Trap or deal with instances when dt=0 in diabatic Added error messages and Adcroft reciprocals to trap or handle cases when dt=0 in various routines in MOM_diabatic_driver including diabatic and legacy_diabatic. All answers are bitwise identical. --- .../vertical/MOM_diabatic_driver.F90 | 14 +++++++++++--- 1 file changed, 11 insertions(+), 3 deletions(-) diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index 200d3efdf7..24a529716d 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -422,6 +422,10 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & eaml => eatr ; ebml => ebtr ! inverse time step + if (dt == 0.0) call MOM_error(FATAL, "MOM_diabatic_driver: "// & + "diabatic was called with a zero length timestep.") + if (dt < 0.0) call MOM_error(FATAL, "MOM_diabatic_driver: "// & + "diabatic was called with a negative timestep.") Idt = 1.0 / dt if (.not. associated(CS)) call MOM_error(FATAL, "MOM_diabatic_driver: "// & @@ -1301,6 +1305,10 @@ subroutine legacy_diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_en eaml => eatr ; ebml => ebtr ! inverse time step + if (dt == 0.0) call MOM_error(FATAL, "MOM_diabatic_driver: "// & + "legacy_diabatic was called with a zero length timestep.") + if (dt < 0.0) call MOM_error(FATAL, "MOM_diabatic_driver: "// & + "legacy_diabatic was called with a negative timestep.") Idt = 1.0 / dt if (.not. associated(CS)) call MOM_error(FATAL, "MOM_diabatic_driver: "// & @@ -2506,7 +2514,7 @@ subroutine diagnose_diabatic_diff_tendency(tv, h, temp_old, saln_old, dt, G, GV, integer :: i, j, k, is, ie, js, je, nz is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke - Idt = 1/dt + Idt = 0.0 ; if (dt > 0.0) Idt = 1. / dt work_3d(:,:,:) = 0.0 work_2d(:,:) = 0.0 @@ -2596,7 +2604,7 @@ subroutine diagnose_boundary_forcing_tendency(tv, h, temp_old, saln_old, h_old, integer :: i, j, k, is, ie, js, je, nz is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke - Idt = 1/dt + Idt = 0.0 ; if (dt > 0.0) Idt = 1. / dt work_3d(:,:,:) = 0.0 work_2d(:,:) = 0.0 @@ -2683,7 +2691,7 @@ subroutine diagnose_frazil_tendency(tv, h, temp_old, dt, G, GV, CS) integer :: i, j, k, is, ie, js, je, nz is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke - Idt = 1/dt + Idt = 0.0 ; if (dt > 0.0) Idt = 1. / dt ! temperature tendency if (CS%id_frazil_temp_tend > 0) then From b6fa3428e84a54a2469aa24ad1d2e644b1270cef Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Wed, 27 Feb 2019 18:24:06 -0500 Subject: [PATCH 3/5] +Opt args to step_MOM can override DIABATIC_FIRST Allow optional arguments to step_MOM to effectively override the value of DIABATIC_FIRST. All answers are bitwise identical in the existing test cases. --- src/core/MOM.F90 | 15 +++++++++++++-- 1 file changed, 13 insertions(+), 2 deletions(-) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 1b364d7ef0..21ec2e6dc6 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -734,8 +734,13 @@ subroutine step_MOM(forces, fluxes, sfc_state, Time_start, time_interval, CS, & !=========================================================================== ! This is the second place where the diabatic processes and remapping could occur. - if (CS%t_dyn_rel_adv == 0.0 .and. do_thermo .and. .not.CS%diabatic_first) then + if ((CS%t_dyn_rel_adv==0.0) .and. do_thermo .and. (.not.CS%diabatic_first)) then + dtdia = CS%t_dyn_rel_thermo + ! If the MOM6 dynamic and thermodynamic time stepping is being orchestrated + ! by the coupler, the value of diabatic_first does not matter. + if ((CS%t_dyn_rel_thermo==0.0) .and. .not.do_dyn) dtdia = dt + if (CS%thermo_spans_coupling .and. (CS%dt_therm > 1.5*cycle_time) .and. & (abs(dt_therm - dtdia) > 1e-6*dt_therm)) then call MOM_error(FATAL, "step_MOM: Mismatch between dt_therm and dtdia "//& @@ -750,7 +755,13 @@ subroutine step_MOM(forces, fluxes, sfc_state, Time_start, time_interval, CS, & call step_MOM_thermo(CS, G, GV, US, u, v, h, CS%tv, fluxes, dtdia, & Time_local, .false., Waves=Waves) CS%time_in_thermo_cycle = CS%time_in_thermo_cycle + dtdia - CS%t_dyn_rel_thermo = 0.0 + + if ((CS%t_dyn_rel_thermo==0.0) .and. .not.do_dyn) then + ! The diabatic processes are now ahead of the dynamics by dtdia. + CS%t_dyn_rel_thermo = -dtdia + else ! The diabatic processes and the dynamics are synchronized. + CS%t_dyn_rel_thermo = 0.0 + endif if (dtdia > dt) & ! Reset CS%Time to its previous value. CS%Time = Time_start + real_to_time(rel_time - 0.5*dt) From 1603ed98c96430b07f5684f941a7c50f76d1300e Mon Sep 17 00:00:00 2001 From: Niki Zadeh Date: Fri, 1 Mar 2019 15:46:42 -0500 Subject: [PATCH 4/5] Fix bug causing openmp answers change - The new structure US was missing a openmp share() declaration --- src/parameterizations/vertical/MOM_vert_friction.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/parameterizations/vertical/MOM_vert_friction.F90 b/src/parameterizations/vertical/MOM_vert_friction.F90 index 5400f8c1df..e01374b5c6 100644 --- a/src/parameterizations/vertical/MOM_vert_friction.F90 +++ b/src/parameterizations/vertical/MOM_vert_friction.F90 @@ -669,7 +669,7 @@ subroutine vertvisc_coef(u, v, h, forces, visc, dt, G, GV, US, CS, OBC) allocate(CS%a1_shelf_v(G%isd:G%ied,G%JsdB:G%JedB)) ; CS%a1_shelf_v(:,:)=0.0 endif - !$OMP parallel do default(private) shared(G,GV,CS,visc,Isq,Ieq,nz,u,h,forces,hML_u, & + !$OMP parallel do default(private) shared(G,GV,US,CS,visc,Isq,Ieq,nz,u,h,forces,hML_u, & !$OMP OBC,h_neglect,dt,I_valBL,Kv_u) & !$OMP firstprivate(i_hbbl) do j=G%Jsc,G%Jec @@ -836,7 +836,7 @@ subroutine vertvisc_coef(u, v, h, forces, visc, dt, G, GV, US, CS, OBC) ! Now work on v-points. - !$OMP parallel do default(private) shared(G,GV,CS,visc,is,ie,Jsq,Jeq,nz,v,h,forces,hML_v, & + !$OMP parallel do default(private) shared(G,GV,CS,US,visc,is,ie,Jsq,Jeq,nz,v,h,forces,hML_v, & !$OMP OBC,h_neglect,dt,I_valBL,Kv_v) & !$OMP firstprivate(i_hbbl) do J=Jsq,Jeq From 0a1a313b3bcef90def5277fb619287573c24f5da Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Mon, 11 Mar 2019 19:24:51 -0400 Subject: [PATCH 5/5] Correct thermo clock with update_ocean_model calls Store the thermodynamic clock at the start of update_ocean_model so that it does not get incorrectly updated during purely dynamic steps. When both dynamics and thermodynamics are stepped together, nothing is changed, but when they are stepped separately, this correct the thermodynamic clock that is used primarily for diagnostics. All answers are bitwise identical in the existing test cases. --- config_src/coupled_driver/ocean_model_MOM.F90 | 13 +++++++++---- 1 file changed, 9 insertions(+), 4 deletions(-) diff --git a/config_src/coupled_driver/ocean_model_MOM.F90 b/config_src/coupled_driver/ocean_model_MOM.F90 index 5af0b774b0..b62f479354 100644 --- a/config_src/coupled_driver/ocean_model_MOM.F90 +++ b/config_src/coupled_driver/ocean_model_MOM.F90 @@ -442,8 +442,12 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, time_start_upda real, optional, intent(in) :: cycle_length !< The duration of a coupled time stepping cycle [s]. ! Local variables - type(time_type) :: Time_seg_start ! Stores the ocean model time at the start of this call to allow - ! step_MOM to temporarily change the time as seen by internal modules. + type(time_type) :: Time_seg_start ! Stores the dynamic or thermodynamic ocean model time at the + ! start of this call to allow step_MOM to temporarily change the time + ! as seen by internal modules. + type(time_type) :: Time_thermo_start ! Stores the ocean model thermodynamics time at the start of + ! this call to allow step_MOM to temporarily change the time as seen by + ! internal modules. type(time_type) :: Time1 ! The value of the ocean model's time at the start of a call to step_MOM. integer :: index_bnds(4) ! The computational domain index bounds in the ice-ocean boundary type. real :: weight ! Flux accumulation weight of the current fluxes. @@ -563,6 +567,7 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, time_start_upda call finish_MOM_initialization(OS%Time, OS%dirs, OS%MOM_CSp, OS%restart_CSp) endif + Time_thermo_start = OS%Time Time_seg_start = OS%Time ; if (do_dyn) Time_seg_start = OS%Time_dyn Time1 = Time_seg_start @@ -576,7 +581,7 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, time_start_upda reset_therm=Ocn_fluxes_used) elseif (OS%single_step_call) then call step_MOM(OS%forces, OS%fluxes, OS%sfc_state, Time1, dt_coupling, OS%MOM_CSp, Waves=OS%Waves) - else + else ! Step both the dynamics and thermodynamics with separate calls. n_max = 1 ; if (dt_coupling > OS%dt) n_max = ceiling(dt_coupling/OS%dt - 0.001) dt_dyn = dt_coupling / real(n_max) thermo_does_span_coupling = (OS%thermo_spans_coupling .and. (OS%dt_therm > 1.5*dt_coupling)) @@ -636,7 +641,7 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, time_start_upda if (do_dyn) OS%Time_dyn = Time_seg_start + Ocean_coupling_time_step if (do_dyn) OS%nstep = OS%nstep + 1 - OS%Time = Time_seg_start ! Reset the clock to compensate for shared pointers. + OS%Time = Time_thermo_start ! Reset the clock to compensate for shared pointers. if (do_thermo) OS%Time = OS%Time + Ocean_coupling_time_step if (do_thermo) OS%nstep_thermo = OS%nstep_thermo + 1