From 1daad4469b4bc63763546789b7cd8a57ef73e58e Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Fri, 31 Jul 2020 14:48:52 -0400 Subject: [PATCH 1/7] (*)Improve make_frazil Improve the handling of very thin layers in make_frazil. This does not change answers for typical values of ANGSTROM, but can avoid problems that can arise when ANGSTROM=0. All answers in the existing MOM6-examples test cases are bitwise identical. --- src/parameterizations/vertical/MOM_diabatic_aux.F90 | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/parameterizations/vertical/MOM_diabatic_aux.F90 b/src/parameterizations/vertical/MOM_diabatic_aux.F90 index 91085047c9..ee9a7bacff 100644 --- a/src/parameterizations/vertical/MOM_diabatic_aux.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_aux.F90 @@ -195,14 +195,14 @@ subroutine make_frazil(h, tv, G, GV, US, CS, p_surf, halo) endif hc = (tv%C_p*GV%H_to_RZ) * h(i,j,k) - if (h(i,j,k) <= 10.0*GV%Angstrom_H) then + if (h(i,j,k) <= 10.0*(GV%Angstrom_H + GV%H_subroundoff)) then ! Very thin layers should not be cooled by the frazil flux. if (tv%T(i,j,k) < T_freeze(i)) then fraz_col(i) = fraz_col(i) + hc * (T_freeze(i) - tv%T(i,j,k)) tv%T(i,j,k) = T_freeze(i) endif - else - if (fraz_col(i) + hc * (T_freeze(i) - tv%T(i,j,k)) <= 0.0) then + elseif ((fraz_col(i) > 0.0) .or. (tv%T(i,j,k) < T_freeze(i))) then + if (fraz_col(i) + hc * (T_freeze(i) - tv%T(i,j,k)) < 0.0) then tv%T(i,j,k) = tv%T(i,j,k) - fraz_col(i) / hc fraz_col(i) = 0.0 else From 948e2926f5c3288283266faf6bc5ca8d6bd25c3d Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Fri, 31 Jul 2020 14:50:03 -0400 Subject: [PATCH 2/7] (*)Improve advective CFL calculation with tiny h Improved handling of massless layers in the calculation of the advective CFL numbers used in PPM tracer advection by using an Adcroft reciprocal instead of adding a small value in the denominator. Although all answers are bitwise identical in the existing MOM6-examples test cases, this can avoid problems with tracer advection when ANGSTROM is 0 or very small like those that were recently found in analogous SIS2 code. --- src/tracer/MOM_tracer_advect.F90 | 104 ++++++++++++------------------- 1 file changed, 40 insertions(+), 64 deletions(-) diff --git a/src/tracer/MOM_tracer_advect.F90 b/src/tracer/MOM_tracer_advect.F90 index 6a362d4fd5..e9c8fb0e7b 100644 --- a/src/tracer/MOM_tracer_advect.F90 +++ b/src/tracer/MOM_tracer_advect.F90 @@ -140,16 +140,15 @@ subroutine advect_tracer(h_end, uhtr, vhtr, OBC, dt, G, GV, US, CS, Reg, & !$OMP hprev,domore_k,js,je,is,ie,uhtr,vhtr,G,GV,h_end,& !$OMP uh_neglect,vh_neglect,ntr,Tr,h_prev_opt) -! This initializes the halos of uhr and vhr because pass_vector might do -! calculations on them, even though they are never used. -!$OMP do - + ! This initializes the halos of uhr and vhr because pass_vector might do + ! calculations on them, even though they are never used. + !$OMP do do k=1,nz do j=jsd,jed ; do I=IsdB,IedB ; uhr(I,j,k) = 0.0 ; enddo ; enddo do J=jsdB,jedB ; do i=Isd,Ied ; vhr(i,J,k) = 0.0 ; enddo ; enddo do j=jsd,jed ; do i=Isd,Ied ; hprev(i,j,k) = 0.0 ; enddo ; enddo domore_k(k)=1 -! Put the remaining (total) thickness fluxes into uhr and vhr. + ! Put the remaining (total) thickness fluxes into uhr and vhr. do j=js,je ; do I=is-1,ie ; uhr(I,j,k) = uhtr(I,j,k) ; enddo ; enddo do J=js-1,je ; do i=is,ie ; vhr(i,J,k) = vhtr(i,J,k) ; enddo ; enddo if (.not. present(h_prev_opt)) then @@ -173,17 +172,17 @@ subroutine advect_tracer(h_end, uhtr, vhtr, OBC, dt, G, GV, US, CS, Reg, & enddo -!$OMP do + !$OMP do do j=jsd,jed ; do I=isd,ied-1 - uh_neglect(I,j) = GV%H_subroundoff*MIN(G%areaT(i,j),G%areaT(i+1,j)) + uh_neglect(I,j) = GV%H_subroundoff * MIN(G%areaT(i,j), G%areaT(i+1,j)) enddo ; enddo -!$OMP do + !$OMP do do J=jsd,jed-1 ; do i=isd,ied - vh_neglect(i,J) = GV%H_subroundoff*MIN(G%areaT(i,j),G%areaT(i,j+1)) + vh_neglect(i,J) = GV%H_subroundoff * MIN(G%areaT(i,j), G%areaT(i,j+1)) enddo ; enddo -!$OMP do ! initialize diagnostic fluxes and tendencies + !$OMP do do m=1,ntr if (associated(Tr(m)%ad_x)) then do k=1,nz ; do j=jsd,jed ; do i=isd,ied @@ -207,7 +206,7 @@ subroutine advect_tracer(h_end, uhtr, vhtr, OBC, dt, G, GV, US, CS, Reg, & do J=jsd,jed ; do i=isd,ied ; Tr(m)%ad2d_y(i,J) = 0.0 ; enddo ; enddo endif enddo -!$OMP end parallel + !$OMP end parallel isv = is ; iev = ie ; jsv = js ; jev = je @@ -222,8 +221,7 @@ subroutine advect_tracer(h_end, uhtr, vhtr, OBC, dt, G, GV, US, CS, Reg, & ! Reevaluate domore_u & domore_v unless the valid range is the same size as ! before. Also, do this if there is Strang splitting. if ((nsten_halo > 1) .or. (itt==1)) then -!$OMP parallel do default(none) shared(nz,domore_k,jsv,jev,domore_u,isv,iev,stencil, & -!$OMP uhr,domore_v,vhr) + !$OMP parallel do default(shared) do k=1,nz ; if (domore_k(k) > 0) then do j=jsv,jev ; if (.not.domore_u(j,k)) then do i=isv+stencil-1,iev-stencil; if (uhr(I,j,k) /= 0.0) then @@ -256,9 +254,7 @@ subroutine advect_tracer(h_end, uhtr, vhtr, OBC, dt, G, GV, US, CS, Reg, & ! for all the transport to happen. The sum over domore_k keeps the processors ! synchronized. This may not be very efficient, but it should be reliable. -!$OMP parallel default(private) shared(nz,domore_k,x_first,Tr,hprev,uhr,uh_neglect, & -!$OMP OBC,domore_u,ntr,Idt,isv,iev,jsv,jev,stencil, & -!$OMP G,GV,CS,vhr,vh_neglect,domore_v,US) + !$OMP parallel default(shared) if (x_first) then @@ -305,7 +301,7 @@ subroutine advect_tracer(h_end, uhtr, vhtr, OBC, dt, G, GV, US, CS, Reg, & endif ! x_first -!$OMP end parallel + !$OMP end parallel ! If the advection just isn't finishing after max_iter, move on. if (itt >= max_iter) then @@ -385,6 +381,7 @@ subroutine advect_x(Tr, hprev, uhr, uh_neglect, OBC, domore_u, ntr, Idt, & CFL ! The absolute value of the advective upwind-cell CFL number [nondim]. real :: min_h ! The minimum thickness that can be realized during ! any of the passes [H ~> m or kg m-2]. + real :: tiny_h ! The smallest numerically invertable thickness [H ~> m or kg m-2]. real :: h_neglect ! A thickness that is so small it is usually lost ! in roundoff and can be neglected [H ~> m or kg m-2]. logical :: do_i(SZIB_(G),SZJ_(G)) ! If true, work on given points. @@ -406,16 +403,15 @@ subroutine advect_x(Tr, hprev, uhr, uh_neglect, OBC, domore_u, ntr, Idt, & if (usePPM .and. .not. useHuynh) stencil = 2 min_h = 0.1*GV%Angstrom_H + tiny_h = tiny(min_h) h_neglect = GV%H_subroundoff -! do I=is-1,ie ; ts2(I) = 0.0 ; enddo do I=is-1,ie ; CFL(I) = 0.0 ; enddo do j=js,je ; if (domore_u(j,k)) then domore_u(j,k) = .false. - ! Calculate the i-direction profiles (slopes) of each tracer that - ! is being advected. + ! Calculate the i-direction profiles (slopes) of each tracer that is being advected. if (usePLMslope) then do m=1,ntr ; do i=is-stencil,ie+stencil !if (ABS(Tr(m)%t(i+1,j,k)-Tr(m)%t(i,j,k)) < & @@ -490,33 +486,33 @@ subroutine advect_x(Tr, hprev, uhr, uh_neglect, OBC, domore_u, ntr, Idt, & ! in the cell plus whatever part of its half of the mass flux that ! the flux through the other side does not require. do I=is-1,ie - if (uhr(I,j,k) == 0.0) then + if ((uhr(I,j,k) == 0.0) .or. & + ((uhr(I,j,k) < 0.0) .and. (hprev(i+1,j,k) <= tiny_h)) .or. & + ((uhr(I,j,k) > 0.0) .and. (hprev(i,j,k) <= tiny_h)) ) then uhh(I) = 0.0 CFL(I) = 0.0 elseif (uhr(I,j,k) < 0.0) then hup = hprev(i+1,j,k) - G%areaT(i+1,j)*min_h - hlos = MAX(0.0,uhr(I+1,j,k)) + hlos = MAX(0.0, uhr(I+1,j,k)) if ((((hup - hlos) + uhr(I,j,k)) < 0.0) .and. & ((0.5*hup + uhr(I,j,k)) < 0.0)) then - uhh(I) = MIN(-0.5*hup,-hup+hlos,0.0) + uhh(I) = MIN(-0.5*hup, -hup+hlos, 0.0) domore_u(j,k) = .true. else uhh(I) = uhr(I,j,k) endif - !ts2(I) = 0.5*(1.0 + uhh(I) / (hprev(i+1,j,k) + h_neglect*G%areaT(i+1,j))) - CFL(I) = - uhh(I) / (hprev(i+1,j,k) + h_neglect*G%areaT(i+1,j)) ! CFL is positive + CFL(I) = - uhh(I) / (hprev(i+1,j,k)) ! CFL is positive else hup = hprev(i,j,k) - G%areaT(i,j)*min_h - hlos = MAX(0.0,-uhr(I-1,j,k)) + hlos = MAX(0.0, -uhr(I-1,j,k)) if ((((hup - hlos) - uhr(I,j,k)) < 0.0) .and. & ((0.5*hup - uhr(I,j,k)) < 0.0)) then - uhh(I) = MAX(0.5*hup,hup-hlos,0.0) + uhh(I) = MAX(0.5*hup, hup-hlos, 0.0) domore_u(j,k) = .true. else uhh(I) = uhr(I,j,k) endif - !ts2(I) = 0.5*(1.0 - uhh(I) / (hprev(i,j,k) + h_neglect*G%areaT(i,j))) - CFL(I) = uhh(I) / (hprev(i,j,k) + h_neglect*G%areaT(i,j)) ! CFL is positive + CFL(I) = uhh(I) / (hprev(i,j,k)) ! CFL is positive endif enddo @@ -545,11 +541,11 @@ subroutine advect_x(Tr, hprev, uhr, uh_neglect, OBC, domore_u, ntr, Idt, & dA = aR - aL ; mA = 0.5*( aR + aL ) if (G%mask2dCu(I_up,j)*G%mask2dCu(I_up-1,j)*(Tp-Tc)*(Tc-Tm) <= 0.) then - aL = Tc ; aR = Tc ! PCM for local extremum and bounadry cells + aL = Tc ; aR = Tc ! PCM for local extremum and bounadry cells elseif ( dA*(Tc-mA) > (dA*dA)/6. ) then - aL = 3.*Tc - 2.*aR + aL = 3.*Tc - 2.*aR elseif ( dA*(Tc-mA) < - (dA*dA)/6. ) then - aR = 3.*Tc - 2.*aL + aR = 3.*Tc - 2.*aL endif a6 = 6.*Tc - 3. * (aR + aL) ! Curvature @@ -570,28 +566,17 @@ subroutine advect_x(Tr, hprev, uhr, uh_neglect, OBC, domore_u, ntr, Idt, & !aR = Tr(m)%t(i,j,k) + 0.5 * slope_x(i,m) !flux_x(I,j,m) = uhh(I)*( aR - 0.5 * (aR-aL) * CFL(I) ) ! Alternative implementation of PLM - !aR = Tr(m)%t(i,j,k) + 0.5 * slope_x(i,m) - !flux_x(I,j,m) = uhh(I)*( aR - 0.5 * slope_x(i,m) * CFL(I) ) - ! Alternative implementation of PLM Tc = T_tmp(i,m) flux_x(I,j,m) = uhh(I)*( Tc + 0.5 * slope_x(i,m) * ( 1. - CFL(I) ) ) - ! Original implementation of PLM - !flux_x(I,j,m) = uhh(I)*(Tr(m)%t(i,j,k) + slope_x(i,m)*ts2(I)) else ! Indirect implementation of PLM !aL = Tr(m)%t(i+1,j,k) - 0.5 * slope_x(i+1,m) !aR = Tr(m)%t(i+1,j,k) + 0.5 * slope_x(i+1,m) !flux_x(I,j,m) = uhh(I)*( aL + 0.5 * (aR-aL) * CFL(I) ) ! Alternative implementation of PLM - !aL = Tr(m)%t(i+1,j,k) - 0.5 * slope_x(i+1,m) - !flux_x(I,j,m) = uhh(I)*( aL + 0.5 * slope_x(i+1,m) * CFL(I) ) - ! Alternative implementation of PLM Tc = T_tmp(i+1,m) flux_x(I,j,m) = uhh(I)*( Tc - 0.5 * slope_x(i+1,m) * ( 1. - CFL(I) ) ) - ! Original implementation of PLM - !flux_x(I,j,m) = uhh(I)*(Tr(m)%t(i+1,j,k) - slope_x(i+1,m)*ts2(I)) endif - !ts2(I) = 0.5*(1.0 - uhh(I)/(hprev(i,j,k)+h_neglect*G%areaT(i,j))) enddo ; enddo endif ! usePPM @@ -760,6 +745,7 @@ subroutine advect_y(Tr, hprev, vhr, vh_neglect, OBC, domore_v, ntr, Idt, & CFL ! The absolute value of the advective upwind-cell CFL number [nondim]. real :: min_h ! The minimum thickness that can be realized during ! any of the passes [H ~> m or kg m-2]. + real :: tiny_h ! The smallest numerically invertable thickness [H ~> m or kg m-2]. real :: h_neglect ! A thickness that is so small it is usually lost ! in roundoff and can be neglected [H ~> m or kg m-2]. logical :: do_j_tr(SZJ_(G)) ! If true, calculate the tracer profiles. @@ -777,8 +763,8 @@ subroutine advect_y(Tr, hprev, vhr, vh_neglect, OBC, domore_v, ntr, Idt, & if (usePPM .and. .not. useHuynh) stencil = 2 min_h = 0.1*GV%Angstrom_H + tiny_h = tiny(min_h) h_neglect = GV%H_subroundoff - !do i=is,ie ; ts2(i) = 0.0 ; enddo ! We conditionally perform work on tracer points: calculating the PLM slope, ! and updating tracer concentration within a cell @@ -822,7 +808,7 @@ subroutine advect_y(Tr, hprev, vhr, vh_neglect, OBC, domore_v, ntr, Idt, & ! make a copy of the tracers in case values need to be overridden for OBCs - do j=G%jsd,G%jed; do m=1,ntr; do i=G%isd,G%ied + do j=G%jsd,G%jed ; do m=1,ntr ; do i=G%isd,G%ied T_tmp(i,m,j) = Tr(m)%t(i,j,k) enddo ; enddo ; enddo @@ -873,33 +859,33 @@ subroutine advect_y(Tr, hprev, vhr, vh_neglect, OBC, domore_v, ntr, Idt, & domore_v(J,k) = .false. do i=is,ie - if (vhr(i,J,k) == 0.0) then + if ((vhr(i,J,k) == 0.0) .or. & + ((vhr(i,J,k) < 0.0) .and. (hprev(i,j+1,k) <= tiny_h)) .or. & + ((vhr(i,J,k) > 0.0) .and. (hprev(i,j,k) <= tiny_h)) ) then vhh(i,J) = 0.0 CFL(i) = 0.0 elseif (vhr(i,J,k) < 0.0) then hup = hprev(i,j+1,k) - G%areaT(i,j+1)*min_h - hlos = MAX(0.0,vhr(i,J+1,k)) + hlos = MAX(0.0, vhr(i,J+1,k)) if ((((hup - hlos) + vhr(i,J,k)) < 0.0) .and. & ((0.5*hup + vhr(i,J,k)) < 0.0)) then - vhh(i,J) = MIN(-0.5*hup,-hup+hlos,0.0) + vhh(i,J) = MIN(-0.5*hup, -hup+hlos, 0.0) domore_v(J,k) = .true. else vhh(i,J) = vhr(i,J,k) endif - !ts2(i) = 0.5*(1.0 + vhh(i,J) / (hprev(i,j+1,k) + h_neglect*G%areaT(i,j+1)) - CFL(i) = - vhh(i,J) / (hprev(i,j+1,k) + h_neglect*G%areaT(i,j+1)) ! CFL is positive + CFL(i) = - vhh(i,J) / hprev(i,j+1,k) ! CFL is positive else hup = hprev(i,j,k) - G%areaT(i,j)*min_h - hlos = MAX(0.0,-vhr(i,J-1,k)) + hlos = MAX(0.0, -vhr(i,J-1,k)) if ((((hup - hlos) - vhr(i,J,k)) < 0.0) .and. & ((0.5*hup - vhr(i,J,k)) < 0.0)) then - vhh(i,J) = MAX(0.5*hup,hup-hlos,0.0) + vhh(i,J) = MAX(0.5*hup, hup-hlos, 0.0) domore_v(J,k) = .true. else vhh(i,J) = vhr(i,J,k) endif - !ts2(i) = 0.5*(1.0 - vhh(i,J) / (hprev(i,j,k) + h_neglect*G%areaT(i,j))) - CFL(i) = vhh(i,J) / (hprev(i,j,k) + h_neglect*G%areaT(i,j)) ! CFL is positive + CFL(i) = vhh(i,J) / hprev(i,j,k) ! CFL is positive endif enddo @@ -952,26 +938,16 @@ subroutine advect_y(Tr, hprev, vhr, vh_neglect, OBC, domore_v, ntr, Idt, & !aR = Tr(m)%t(i,j,k) + 0.5 * slope_y(i,m,j) !flux_y(i,m,J) = vhh(i,J)*( aR - 0.5 * (aR-aL) * CFL(i) ) ! Alternative implementation of PLM - !aR = Tr(m)%t(i,j,k) + 0.5 * slope_y(i,m,j) - !flux_y(i,m,J) = vhh(i,J)*(aR - 0.5 * slope_y(i,m,j)*CFL(i)) - ! Alternative implementation of PLM Tc = T_tmp(i,m,j) flux_y(i,m,J) = vhh(i,J)*( Tc + 0.5 * slope_y(i,m,j) * ( 1. - CFL(i) ) ) - ! Original implementation of PLM - !flux_y(i,m,J) = vhh(i,J)*(Tr(m)%t(i,j,k) + slope_y(i,m,j)*ts2(i)) else ! Indirect implementation of PLM !aL = Tr(m)%t(i,j+1,k) - 0.5 * slope_y(i,m,j+1) !aR = Tr(m)%t(i,j+1,k) + 0.5 * slope_y(i,m,j+1) !flux_y(i,m,J) = vhh(i,J)*( aL + 0.5 * (aR-aL) * CFL(i) ) ! Alternative implementation of PLM - !aL = Tr(m)%t(i,j+1,k) - 0.5 * slope_y(i,m,j+1) - !flux_y(i,m,J) = vhh(i,J)*( aL + 0.5 * slope_y(i,m,j+1)*CFL(i) ) - ! Alternative implementation of PLM Tc = T_tmp(i,m,j+1) flux_y(i,m,J) = vhh(i,J)*( Tc - 0.5 * slope_y(i,m,j+1) * ( 1. - CFL(i) ) ) - ! Original implementation of PLM - !flux_y(i,m,J) = vhh(i,J)*(Tr(m)%t(i,j+1,k) - slope_y(i,m,j+1)*ts2(i)) endif enddo ; enddo endif ! usePPM From cbbf84847384f2860bf8de765b5c1cb34687cb75 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Fri, 31 Jul 2020 15:02:28 -0400 Subject: [PATCH 3/7] Infrastructure calls via framework directory Revised module use statements and some infrastructure calls to go via the MOM6 framework directory rather than directly calling FMS infrastructure routines. All answers are bitwise identical. --- .../MOM_surface_forcing_gfdl.F90 | 4 +-- src/ice_shelf/MOM_ice_shelf.F90 | 1 - src/ice_shelf/MOM_ice_shelf_state.F90 | 1 - .../MOM_state_initialization.F90 | 20 ++++-------- src/tracer/MOM_generic_tracer.F90 | 32 +++++++++---------- src/tracer/MOM_offline_aux.F90 | 3 +- src/tracer/MOM_offline_main.F90 | 3 +- src/tracer/RGC_tracer.F90 | 5 ++- 8 files changed, 26 insertions(+), 43 deletions(-) diff --git a/config_src/coupled_driver/MOM_surface_forcing_gfdl.F90 b/config_src/coupled_driver/MOM_surface_forcing_gfdl.F90 index a04ee426e6..7075fb7c10 100644 --- a/config_src/coupled_driver/MOM_surface_forcing_gfdl.F90 +++ b/config_src/coupled_driver/MOM_surface_forcing_gfdl.F90 @@ -9,10 +9,8 @@ module MOM_surface_forcing_gfdl use MOM_constants, only : hlv, hlf use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end use MOM_cpu_clock, only : CLOCK_SUBCOMPONENT -use MOM_diag_mediator, only : diag_ctrl -use MOM_diag_mediator, only : safe_alloc_ptr, time_type +use MOM_diag_mediator, only : diag_ctrl, safe_alloc_ptr, time_type use MOM_domains, only : pass_vector, pass_var, fill_symmetric_edges -use MOM_domains, only : global_field_sum, BITWISE_EXACT_SUM use MOM_domains, only : AGRID, BGRID_NE, CGRID_NE, To_All use MOM_domains, only : To_North, To_East, Omit_Corners use MOM_error_handler, only : MOM_error, WARNING, FATAL, is_root_pe, MOM_mesg diff --git a/src/ice_shelf/MOM_ice_shelf.F90 b/src/ice_shelf/MOM_ice_shelf.F90 index 6b68cb3deb..66fd873f67 100644 --- a/src/ice_shelf/MOM_ice_shelf.F90 +++ b/src/ice_shelf/MOM_ice_shelf.F90 @@ -51,7 +51,6 @@ module MOM_ice_shelf use MOM_checksums, only : hchksum, qchksum, chksum, uchksum, vchksum, uvchksum use time_interp_external_mod, only : init_external_field, time_interp_external use time_interp_external_mod, only : time_interp_external_init -use time_manager_mod, only : print_time implicit none ; private #include diff --git a/src/ice_shelf/MOM_ice_shelf_state.F90 b/src/ice_shelf/MOM_ice_shelf_state.F90 index b3e88697f2..a3784b5a34 100644 --- a/src/ice_shelf/MOM_ice_shelf_state.F90 +++ b/src/ice_shelf/MOM_ice_shelf_state.F90 @@ -12,7 +12,6 @@ module MOM_ice_shelf_state use MOM_file_parser, only : read_param, get_param, log_param, log_version, param_file_type use MOM_grid, only : MOM_grid_init, ocean_grid_type use MOM_get_input, only : directories, Get_MOM_input -use mpp_mod, only : mpp_sum, mpp_max, mpp_min, mpp_pe, mpp_npes, mpp_sync use MOM_coms, only : reproducing_sum use MOM_checksums, only : hchksum, qchksum, chksum, uchksum, vchksum, uvchksum diff --git a/src/initialization/MOM_state_initialization.F90 b/src/initialization/MOM_state_initialization.F90 index e451966364..a201e4a85f 100644 --- a/src/initialization/MOM_state_initialization.F90 +++ b/src/initialization/MOM_state_initialization.F90 @@ -18,23 +18,17 @@ module MOM_state_initialization use MOM_get_input, only : directories use MOM_grid, only : ocean_grid_type, isPointInCell use MOM_interface_heights, only : find_eta -use MOM_io, only : file_exists -use MOM_io, only : MOM_read_data, MOM_read_vector -use MOM_io, only : slasher -use MOM_open_boundary, only : ocean_OBC_type, open_boundary_init +use MOM_io, only : file_exists, field_size, MOM_read_data, MOM_read_vector, slasher +use MOM_open_boundary, only : ocean_OBC_type, open_boundary_init, set_tracer_data use MOM_open_boundary, only : OBC_NONE, OBC_SIMPLE -use MOM_open_boundary, only : open_boundary_query -use MOM_open_boundary, only : set_tracer_data -use MOM_open_boundary, only : open_boundary_test_extern_h -use MOM_open_boundary, only : fill_temp_salt_segments -use MOM_open_boundary, only : update_OBC_segment_data +use MOM_open_boundary, only : open_boundary_query, open_boundary_test_extern_h +use MOM_open_boundary, only : fill_temp_salt_segments, update_OBC_segment_data !use MOM_open_boundary, only : set_3D_OBC_data use MOM_grid_initialize, only : initialize_masks, set_grid_metrics use MOM_restart, only : restore_state, determine_is_new_run, MOM_restart_CS use MOM_sponge, only : set_up_sponge_field, set_up_sponge_ML_density use MOM_sponge, only : initialize_sponge, sponge_CS -use MOM_ALE_sponge, only : set_up_ALE_sponge_field, initialize_ALE_sponge -use MOM_ALE_sponge, only : ALE_sponge_CS +use MOM_ALE_sponge, only : set_up_ALE_sponge_field, initialize_ALE_sponge, ALE_sponge_CS use MOM_string_functions, only : uppercase, lowercase use MOM_time_manager, only : time_type use MOM_tracer_registry, only : tracer_registry_type @@ -44,8 +38,7 @@ module MOM_state_initialization use MOM_EOS, only : calculate_density, calculate_density_derivs, EOS_type, EOS_domain use MOM_EOS, only : convert_temp_salt_for_TEOS10 use user_initialization, only : user_initialize_thickness, user_initialize_velocity -use user_initialization, only : user_init_temperature_salinity -use user_initialization, only : user_set_OBC_data +use user_initialization, only : user_init_temperature_salinity, user_set_OBC_data use user_initialization, only : user_initialize_sponges use DOME_initialization, only : DOME_initialize_thickness use DOME_initialization, only : DOME_set_OBC_data @@ -97,7 +90,6 @@ module MOM_state_initialization use MOM_remapping, only : remapping_CS, initialize_remapping use MOM_remapping, only : remapping_core_h use MOM_horizontal_regridding, only : horiz_interp_and_extrap_tracer -use fms_io_mod, only : field_size implicit none ; private diff --git a/src/tracer/MOM_generic_tracer.F90 b/src/tracer/MOM_generic_tracer.F90 index 7d2310b42f..66c0e33bac 100644 --- a/src/tracer/MOM_generic_tracer.F90 +++ b/src/tracer/MOM_generic_tracer.F90 @@ -13,10 +13,8 @@ module MOM_generic_tracer #define _ALLOCATED allocated #endif - ! ### These imports should not reach into FMS directly ### - use mpp_mod, only: stdout, mpp_error, FATAL,WARNING,NOTE - use field_manager_mod, only: fm_get_index,fm_string_len + use field_manager_mod, only: fm_string_len use generic_tracer, only: generic_tracer_register, generic_tracer_get_diag_list use generic_tracer, only: generic_tracer_init, generic_tracer_source, generic_tracer_register_diag @@ -108,7 +106,7 @@ function register_MOM_generic_tracer(HI, GV, param_file, CS, tr_Reg, restart_CS) ! Local variables logical :: register_MOM_generic_tracer - character(len=fm_string_len), parameter :: sub_name = 'register_MOM_generic_tracer' + character(len=128), parameter :: sub_name = 'register_MOM_generic_tracer' character(len=200) :: inputdir ! The directory where NetCDF input files are. ! These can be overridden later in via the field manager? @@ -122,7 +120,7 @@ function register_MOM_generic_tracer(HI, GV, param_file, CS, tr_Reg, restart_CS) register_MOM_generic_tracer = .false. if (associated(CS)) then - call mpp_error(WARNING, "register_MOM_generic_tracer called with an "// & + call MOM_error(WARNING, "register_MOM_generic_tracer called with an "// & "associated control structure.") return endif @@ -185,7 +183,7 @@ function register_MOM_generic_tracer(HI, GV, param_file, CS, tr_Reg, restart_CS) !Get the tracer list call generic_tracer_get_list(CS%g_tracer_list) - if (.NOT. associated(CS%g_tracer_list)) call mpp_error(FATAL, trim(sub_name)//& + if (.NOT. associated(CS%g_tracer_list)) call MOM_error(FATAL, trim(sub_name)//& ": No tracer in the list.") ! For each tracer name get its T_prog index and get its fields @@ -247,7 +245,7 @@ subroutine initialize_MOM_generic_tracer(restart, day, G, GV, US, h, param_file, type(ALE_sponge_CS), pointer :: ALE_sponge_CSp !< Pointer to the control structure for the !! ALE sponges. - character(len=fm_string_len), parameter :: sub_name = 'initialize_MOM_generic_tracer' + character(len=128), parameter :: sub_name = 'initialize_MOM_generic_tracer' logical :: OK integer :: i, j, k, isc, iec, jsc, jec, nk type(g_tracer_type), pointer :: g_tracer,g_tracer_next @@ -265,7 +263,7 @@ subroutine initialize_MOM_generic_tracer(restart, day, G, GV, US, h, param_file, CS%diag=>diag !Get the tracer list - if (.NOT. associated(CS%g_tracer_list)) call mpp_error(FATAL, trim(sub_name)//& + if (.NOT. associated(CS%g_tracer_list)) call MOM_error(FATAL, trim(sub_name)//& ": No tracer in the list.") !For each tracer name get its fields g_tracer=>CS%g_tracer_list @@ -426,7 +424,7 @@ subroutine MOM_generic_tracer_column_physics(h_old, h_new, ea, eb, fluxes, Hml, ! h_new(k) = h_old(k) + ea(k) - eb(k-1) + eb(k) - ea(k+1) ! Local variables - character(len=fm_string_len), parameter :: sub_name = 'MOM_generic_tracer_column_physics' + character(len=128), parameter :: sub_name = 'MOM_generic_tracer_column_physics' type(g_tracer_type), pointer :: g_tracer, g_tracer_next character(len=fm_string_len) :: g_tracer_name @@ -443,7 +441,7 @@ subroutine MOM_generic_tracer_column_physics(h_old, h_new, ea, eb, fluxes, Hml, isc = G%isc ; iec = G%iec ; jsc = G%jsc ; jec = G%jec ; nk = G%ke !Get the tracer list - if (.NOT. associated(CS%g_tracer_list)) call mpp_error(FATAL,& + if (.NOT. associated(CS%g_tracer_list)) call MOM_error(FATAL,& trim(sub_name)//": No tracer in the list.") #ifdef _USE_MOM6_DIAG @@ -587,7 +585,7 @@ function MOM_generic_tracer_stock(h, stocks, G, GV, CS, names, units, stock_inde type(g_tracer_type), pointer :: g_tracer, g_tracer_next real, dimension(:,:,:,:), pointer :: tr_field real, dimension(:,:,:), pointer :: tr_ptr - character(len=fm_string_len), parameter :: sub_name = 'MOM_generic_tracer_stock' + character(len=128), parameter :: sub_name = 'MOM_generic_tracer_stock' integer :: i, j, k, is, ie, js, je, nz, m is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke @@ -660,7 +658,7 @@ function MOM_generic_tracer_min_max(ind_start, got_minmax, gmin, gmax, xgmin, yg type(g_tracer_type), pointer :: g_tracer, g_tracer_next real, dimension(:,:,:,:), pointer :: tr_field real, dimension(:,:,:), pointer :: tr_ptr - character(len=fm_string_len), parameter :: sub_name = 'MOM_generic_tracer_min_max' + character(len=128), parameter :: sub_name = 'MOM_generic_tracer_min_max' real, dimension(:,:,:),pointer :: grid_tmask integer :: isc,iec,jsc,jec,isd,ied,jsd,jed,nk,ntau @@ -728,7 +726,7 @@ subroutine MOM_generic_tracer_surface_state(sfc_state, h, G, CS) ! Local variables real :: sosga - character(len=fm_string_len), parameter :: sub_name = 'MOM_generic_tracer_surface_state' + character(len=128), parameter :: sub_name = 'MOM_generic_tracer_surface_state' real, dimension(G%isd:G%ied,G%jsd:G%jed,1:G%ke,1) :: rho0 real, dimension(G%isd:G%ied,G%jsd:G%jed,1:G%ke) :: dzt type(g_tracer_type), pointer :: g_tracer @@ -750,7 +748,7 @@ subroutine MOM_generic_tracer_surface_state(sfc_state, h, G, CS) tau=1,sosga=sosga,model_time=get_diag_time_end(CS%diag)) !Output diagnostics via diag_manager for all tracers in this module -! if (.NOT. associated(CS%g_tracer_list)) call mpp_error(FATAL, trim(sub_name)//& +! if (.NOT. associated(CS%g_tracer_list)) call MOM_error(FATAL, trim(sub_name)//& ! "No tracer in the list.") ! call g_tracer_send_diag(CS%g_tracer_list, get_diag_time_end(CS%diag), tau=1) !Niki: The problem with calling diagnostic outputs here is that this subroutine is called every dt_cpld @@ -767,7 +765,7 @@ subroutine MOM_generic_flux_init(verbosity) integer :: ind character(len=fm_string_len) :: g_tracer_name,longname, package,units,old_package,file_in,file_out real :: const_init_value - character(len=fm_string_len), parameter :: sub_name = 'MOM_generic_flux_init' + character(len=128), parameter :: sub_name = 'MOM_generic_flux_init' type(g_tracer_type), pointer :: g_tracer_list,g_tracer,g_tracer_next if (.not. g_registered) then @@ -777,7 +775,7 @@ subroutine MOM_generic_flux_init(verbosity) call generic_tracer_get_list(g_tracer_list) if (.NOT. associated(g_tracer_list)) then - call mpp_error(WARNING, trim(sub_name)// ": No generic tracer in the list.") + call MOM_error(WARNING, trim(sub_name)// ": No generic tracer in the list.") return endif @@ -812,7 +810,7 @@ subroutine MOM_generic_tracer_get(name,member,array, CS) type(MOM_generic_tracer_CS), pointer :: CS !< Pointer to the control structure for this module. real, dimension(:,:,:), pointer :: array_ptr - character(len=fm_string_len), parameter :: sub_name = 'MOM_generic_tracer_get' + character(len=128), parameter :: sub_name = 'MOM_generic_tracer_get' call g_tracer_get_pointer(CS%g_tracer_list,name,member,array_ptr) array(:,:,:) = array_ptr(:,:,:) diff --git a/src/tracer/MOM_offline_aux.F90 b/src/tracer/MOM_offline_aux.F90 index 21db2cfff4..119ad555da 100644 --- a/src/tracer/MOM_offline_aux.F90 +++ b/src/tracer/MOM_offline_aux.F90 @@ -4,7 +4,6 @@ module MOM_offline_aux ! This file is part of MOM6. See LICENSE.md for the license. -use mpp_domains_mod, only : CENTER, CORNER, NORTH, EAST use data_override_mod, only : data_override_init, data_override use MOM_time_manager, only : time_type, operator(-) use MOM_debugging, only : check_column_integrals @@ -12,7 +11,7 @@ module MOM_offline_aux use MOM_diag_vkernels, only : reintegrate_column use MOM_error_handler, only : callTree_enter, callTree_leave, MOM_error, FATAL, WARNING, is_root_pe use MOM_grid, only : ocean_grid_type -use MOM_io, only : MOM_read_data, MOM_read_vector +use MOM_io, only : MOM_read_data, MOM_read_vector, CENTER use MOM_verticalGrid, only : verticalGrid_type use MOM_file_parser, only : get_param, log_version, param_file_type use astronomy_mod, only : orbital_time, diurnal_solar, daily_mean_solar diff --git a/src/tracer/MOM_offline_main.F90 b/src/tracer/MOM_offline_main.F90 index b7af9849b3..3895e8a116 100644 --- a/src/tracer/MOM_offline_main.F90 +++ b/src/tracer/MOM_offline_main.F90 @@ -4,7 +4,6 @@ module MOM_offline_main ! This file is part of MOM6. See LICENSE.md for the license. -use mpp_domains_mod, only : CENTER, CORNER, NORTH, EAST use MOM_ALE, only : ALE_CS, ALE_main_offline, ALE_offline_inputs use MOM_checksums, only : hchksum, uvchksum use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end @@ -20,7 +19,7 @@ module MOM_offline_main use MOM_file_parser, only : read_param, get_param, log_version, param_file_type use MOM_forcing_type, only : forcing use MOM_grid, only : ocean_grid_type -use MOM_io, only : MOM_read_data, MOM_read_vector +use MOM_io, only : MOM_read_data, MOM_read_vector, CENTER use MOM_offline_aux, only : update_offline_from_arrays, update_offline_from_files use MOM_offline_aux, only : next_modulo_time, offline_add_diurnal_sw use MOM_offline_aux, only : update_h_horizontal_flux, update_h_vertical_flux, limit_mass_flux_3d diff --git a/src/tracer/RGC_tracer.F90 b/src/tracer/RGC_tracer.F90 index 028718f379..44c6c2e5a1 100644 --- a/src/tracer/RGC_tracer.F90 +++ b/src/tracer/RGC_tracer.F90 @@ -19,7 +19,7 @@ module RGC_tracer use MOM_forcing_type, only : forcing use MOM_hor_index, only : hor_index_type use MOM_grid, only : ocean_grid_type -use MOM_io, only : file_exists, read_data, slasher, vardesc, var_desc, query_vardesc +use MOM_io, only : file_exists, MOM_read_data, slasher, vardesc, var_desc, query_vardesc use MOM_restart, only : MOM_restart_CS use MOM_ALE_sponge, only : set_up_ALE_sponge_field, ALE_sponge_CS, get_ALE_sponge_nz_data use MOM_sponge, only : set_up_sponge_field, sponge_CS @@ -207,8 +207,7 @@ subroutine initialize_RGC_tracer(restart, day, G, GV, h, diag, OBC, CS, & CS%tracer_IC_file) do m=1,NTR call query_vardesc(CS%tr_desc(m), name, caller="initialize_RGC_tracer") - call read_data(CS%tracer_IC_file, trim(name), & - CS%tr(:,:,:,m), domain=G%Domain%mpp_domain) + call MOM_read_data(CS%tracer_IC_file, trim(name), CS%tr(:,:,:,m), G%Domain) enddo else do m=1,NTR From bba60af8efe150e65accdb2b2621614b707f8316 Mon Sep 17 00:00:00 2001 From: Matthew Harrison Date: Tue, 4 Aug 2020 12:06:41 -0400 Subject: [PATCH 4/7] Move call to initialize_segment_data to MOM_state_initialization --- src/core/MOM_open_boundary.F90 | 96 +++++++++++-------- .../MOM_state_initialization.F90 | 4 +- 2 files changed, 58 insertions(+), 42 deletions(-) diff --git a/src/core/MOM_open_boundary.F90 b/src/core/MOM_open_boundary.F90 index 9b650f8598..f94060fc39 100644 --- a/src/core/MOM_open_boundary.F90 +++ b/src/core/MOM_open_boundary.F90 @@ -62,6 +62,7 @@ module MOM_open_boundary public update_OBC_ramp public rotate_OBC_config public rotate_OBC_init +public initialize_segment_data integer, parameter, public :: OBC_NONE = 0 !< Indicates the use of no open boundary integer, parameter, public :: OBC_SIMPLE = 1 !< Indicates the use of a simple inflow open boundary @@ -268,7 +269,7 @@ module MOM_open_boundary real :: rx_max !< The maximum magnitude of the baroclinic radiation velocity (or speed of !! characteristics) in units of grid points per timestep [nondim]. logical :: OBC_pe !< Is there an open boundary on this tile? - type(remapping_CS), pointer :: remap_CS !< ALE remapping control structure for segments only + type(remapping_CS), pointer :: remap_CS=> NULL() !< ALE remapping control structure for segments only type(OBC_registry_type), pointer :: OBC_Reg => NULL() !< Registry type for boundaries real, pointer, dimension(:,:,:) :: & rx_normal => NULL(), & !< Array storage for normal phase speed for EW radiation OBCs in units of @@ -341,6 +342,11 @@ subroutine open_boundary_config(G, US, param_file, OBC) character(len=100) :: segment_str ! The contents (rhs) for parameter "segment_param_str" character(len=200) :: config1 ! String for OBC_USER_CONFIG real :: Lscale_in, Lscale_out ! parameters controlling tracer values at the boundaries [L ~> m] + character(len=128) :: inputdir + logical :: answers_2018, default_2018_answers + logical :: check_reconstruction, check_remapping, force_bounds_in_subcell + character(len=32) :: remappingScheme + allocate(OBC) call get_param(param_file, mdl, "OBC_NUMBER_OF_SEGMENTS", OBC%number_of_segments, & @@ -497,7 +503,7 @@ subroutine open_boundary_config(G, US, param_file, OBC) enddo ! if (open_boundary_query(OBC, needs_ext_seg_data=.true.)) & - call initialize_segment_data(G, OBC, param_file) + ! call initialize_segment_data(G, OBC, param_file) if (open_boundary_query(OBC, apply_open_OBC=.true.)) then call get_param(param_file, mdl, "OBC_RADIATION_MAX", OBC%rx_max, & @@ -540,6 +546,46 @@ subroutine open_boundary_config(G, US, param_file, OBC) if (Lscale_out>0.) OBC%segment(l)%Tr_InvLscale_out = 1.0/Lscale_out enddo + ! There is a problem with the order of the OBC initialization + ! with respect to ALE_init. Currently handling this by copying the + ! param file so that I can use it later in step_MOM in order to finish + ! initializing segments on the first step. + + ! Is the above comment still relevant ? + + call get_param(param_file, mdl, "REMAPPING_SCHEME", remappingScheme, & + "This sets the reconstruction scheme used "//& + "for vertical remapping for all variables. "//& + "It can be one of the following schemes: \n"//& + trim(remappingSchemesDoc), default=remappingDefaultScheme,do_not_log=.true.) + call get_param(param_file, mdl, "FATAL_CHECK_RECONSTRUCTIONS", check_reconstruction, & + "If true, cell-by-cell reconstructions are checked for "//& + "consistency and if non-monotonicity or an inconsistency is "//& + "detected then a FATAL error is issued.", default=.false.,do_not_log=.true.) + call get_param(param_file, mdl, "FATAL_CHECK_REMAPPING", check_remapping, & + "If true, the results of remapping are checked for "//& + "conservation and new extrema and if an inconsistency is "//& + "detected then a FATAL error is issued.", default=.false.,do_not_log=.true.) + call get_param(param_file, mdl, "BRUSHCUTTER_MODE", OBC%brushcutter_mode, & + "If true, read external OBC data on the supergrid.", & + default=.false.) + call get_param(param_file, mdl, "REMAP_BOUND_INTERMEDIATE_VALUES", force_bounds_in_subcell, & + "If true, the values on the intermediate grid used for remapping "//& + "are forced to be bounded, which might not be the case due to "//& + "round off.", default=.false.,do_not_log=.true.) + call get_param(param_file, mdl, "DEFAULT_2018_ANSWERS", default_2018_answers, & + "This sets the default value for the various _2018_ANSWERS parameters.", & + default=.false.) + call get_param(param_file, mdl, "REMAPPING_2018_ANSWERS", answers_2018, & + "If true, use the order of arithmetic and expressions that recover the "//& + "answers from the end of 2018. Otherwise, use updated and more robust "//& + "forms of the same expressions.", default=default_2018_answers) + + allocate(OBC%remap_CS) + call initialize_remapping(OBC%remap_CS, remappingScheme, boundary_extrapolation = .false., & + check_reconstruction=check_reconstruction, check_remapping=check_remapping, & + force_bounds_in_subcell=force_bounds_in_subcell, answers_2018=answers_2018) + endif ! OBC%number_of_segments > 0 ! Safety check @@ -564,7 +610,7 @@ end subroutine open_boundary_config subroutine initialize_segment_data(G, OBC, PF) use mpp_mod, only : mpp_pe, mpp_set_current_pelist, mpp_get_current_pelist,mpp_npes - type(dyn_horgrid_type), intent(in) :: G !< Ocean grid structure + type(ocean_grid_type), intent(in) :: G !< Ocean grid structure type(ocean_OBC_type), intent(inout) :: OBC !< Open boundary control structure type(param_file_type), intent(in) :: PF !< Parameter file handle @@ -576,10 +622,7 @@ subroutine initialize_segment_data(G, OBC, PF) character(len=32), dimension(MAX_OBC_FIELDS) :: fields ! segment field names character(len=128) :: inputdir type(OBC_segment_type), pointer :: segment => NULL() ! pointer to segment type list - character(len=32) :: remappingScheme character(len=256) :: mesg ! Message for error messages. - logical :: check_reconstruction, check_remapping, force_bounds_in_subcell - logical :: answers_2018, default_2018_answers integer, dimension(4) :: siz,siz2 integer :: is, ie, js, je integer :: isd, ied, jsd, jed @@ -599,39 +642,6 @@ subroutine initialize_segment_data(G, OBC, PF) call get_param(PF, mdl, "INPUTDIR", inputdir, default=".") inputdir = slasher(inputdir) - call get_param(PF, mdl, "REMAPPING_SCHEME", remappingScheme, & - "This sets the reconstruction scheme used "//& - "for vertical remapping for all variables. "//& - "It can be one of the following schemes: \n"//& - trim(remappingSchemesDoc), default=remappingDefaultScheme,do_not_log=.true.) - call get_param(PF, mdl, "FATAL_CHECK_RECONSTRUCTIONS", check_reconstruction, & - "If true, cell-by-cell reconstructions are checked for "//& - "consistency and if non-monotonicity or an inconsistency is "//& - "detected then a FATAL error is issued.", default=.false.,do_not_log=.true.) - call get_param(PF, mdl, "FATAL_CHECK_REMAPPING", check_remapping, & - "If true, the results of remapping are checked for "//& - "conservation and new extrema and if an inconsistency is "//& - "detected then a FATAL error is issued.", default=.false.,do_not_log=.true.) - call get_param(PF, mdl, "REMAP_BOUND_INTERMEDIATE_VALUES", force_bounds_in_subcell, & - "If true, the values on the intermediate grid used for remapping "//& - "are forced to be bounded, which might not be the case due to "//& - "round off.", default=.false.,do_not_log=.true.) - call get_param(PF, mdl, "BRUSHCUTTER_MODE", OBC%brushcutter_mode, & - "If true, read external OBC data on the supergrid.", & - default=.false.) - call get_param(PF, mdl, "DEFAULT_2018_ANSWERS", default_2018_answers, & - "This sets the default value for the various _2018_ANSWERS parameters.", & - default=.false.) - call get_param(PF, mdl, "REMAPPING_2018_ANSWERS", answers_2018, & - "If true, use the order of arithmetic and expressions that recover the "//& - "answers from the end of 2018. Otherwise, use updated and more robust "//& - "forms of the same expressions.", default=default_2018_answers) - - allocate(OBC%remap_CS) - call initialize_remapping(OBC%remap_CS, remappingScheme, boundary_extrapolation = .false., & - check_reconstruction=check_reconstruction, check_remapping=check_remapping, & - force_bounds_in_subcell=force_bounds_in_subcell, answers_2018=answers_2018) - if (OBC%user_BCs_set_globally) return ! Try this here just for the documentation. It is repeated below. @@ -4966,6 +4976,8 @@ subroutine rotate_OBC_config(OBC_in, G_in, OBC, G, turns) integer :: l + if (OBC_in%number_of_segments==0) return + ! Scalar and logical transfer OBC%number_of_segments = OBC_in%number_of_segments OBC%ke = OBC_in%ke @@ -5023,8 +5035,10 @@ subroutine rotate_OBC_config(OBC_in, G_in, OBC, G, turns) OBC%OBC_pe = OBC_in%OBC_pe ! remap_CS is set up by initialize_segment_data, so we copy the fields here. - allocate(OBC%remap_CS) - OBC%remap_CS = OBC_in%remap_CS + if (ASSOCIATED(OBC_in%remap_CS)) then + allocate(OBC%remap_CS) + OBC%remap_CS = OBC_in%remap_CS + endif ! TODO: The OBC registry seems to be a list of "registered" OBC types. ! It does not appear to be used, so for now we skip this record. diff --git a/src/initialization/MOM_state_initialization.F90 b/src/initialization/MOM_state_initialization.F90 index e451966364..f53ff89a1e 100644 --- a/src/initialization/MOM_state_initialization.F90 +++ b/src/initialization/MOM_state_initialization.F90 @@ -24,7 +24,7 @@ module MOM_state_initialization use MOM_open_boundary, only : ocean_OBC_type, open_boundary_init use MOM_open_boundary, only : OBC_NONE, OBC_SIMPLE use MOM_open_boundary, only : open_boundary_query -use MOM_open_boundary, only : set_tracer_data +use MOM_open_boundary, only : set_tracer_data, initialize_segment_data use MOM_open_boundary, only : open_boundary_test_extern_h use MOM_open_boundary, only : fill_temp_salt_segments use MOM_open_boundary, only : update_OBC_segment_data @@ -563,6 +563,8 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, US, PF, dirs, & ! This controls user code for setting open boundary data if (associated(OBC)) then + call initialize_segment_data(G, OBC, PF) ! call initialize_segment_data(G, OBC, param_file) +! call open_boundary_config(G, US, PF, OBC) ! Call this once to fill boundary arrays from fixed values if (.not. OBC%needs_IO_for_data) & call update_OBC_segment_data(G, GV, US, OBC, tv, h, Time) From 7be08832127b4669db8644cfc261ba82163d3a5d Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Fri, 7 Aug 2020 08:13:03 -0400 Subject: [PATCH 5/7] (*)Set dSV_dT and dSV_dS with unassociated fluxes Set dSV_dT and dSV_dS if present in applyBoundaryFluxesInOut, even if boundary fluxes are not associated. With this change, setting BUOY_CONFIG='NONE' and BUOY_CONFIG='zero' both work and give similar (but not identical) answers in some test cases with an ePBL boundary layer parameterization, whereas before answers were tainted with uninitialized values when BUOY_CONFIG='NONE'. All answers in the existing MOM6-examples test suite are bitwise identical, but answers can change in other cases. --- .../vertical/MOM_diabatic_aux.F90 | 15 +++++++++------ .../vertical/MOM_vert_friction.F90 | 2 +- 2 files changed, 10 insertions(+), 7 deletions(-) diff --git a/src/parameterizations/vertical/MOM_diabatic_aux.F90 b/src/parameterizations/vertical/MOM_diabatic_aux.F90 index 91085047c9..bf2e86cb80 100644 --- a/src/parameterizations/vertical/MOM_diabatic_aux.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_aux.F90 @@ -822,19 +822,18 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, t is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke - ! Only apply forcing if fluxes%sw is associated. - if (.not.associated(fluxes%sw)) return - -#define _OLD_ALG_ Idt = 1.0 / dt calculate_energetics = (present(cTKE) .and. present(dSV_dT) .and. present(dSV_dS)) calculate_buoyancy = present(SkinBuoyFlux) if (calculate_buoyancy) SkinBuoyFlux(:,:) = 0.0 + if (present(cTKE)) cTKE(:,:,:) = 0.0 g_Hconv2 = (US%L_to_Z**2*GV%g_Earth * GV%H_to_RZ) * GV%H_to_RZ EOSdom(:) = EOS_domain(G%HI) - if (present(cTKE)) cTKE(:,:,:) = 0.0 + ! Only apply forcing if fluxes%sw is associated. + if (.not.associated(fluxes%sw) .and. .not.calculate_energetics) return + if (calculate_buoyancy) then SurfPressure(:) = 0.0 GoRho = US%L_to_Z**2*GV%g_Earth / GV%Rho0 @@ -874,7 +873,6 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, t h2d(i,k) = h(i,j,k) T2d(i,k) = tv%T(i,j,k) enddo ; enddo - if (nsw>0) call extract_optics_slice(optics, j, G, GV, opacity=opacityBand, opacity_scale=(1.0/GV%m_to_H)) if (calculate_energetics) then ! The partial derivatives of specific volume with temperature and @@ -898,6 +896,11 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, t pen_TKE_2d(:,:) = 0.0 endif + ! Nothing more is done on this j-slice if there is no buoyancy forcing. + if (.not.associated(fluxes%sw)) cycle + + if (nsw>0) call extract_optics_slice(optics, j, G, GV, opacity=opacityBand, opacity_scale=(1.0/GV%m_to_H)) + ! The surface forcing is contained in the fluxes type. ! We aggregate the thermodynamic forcing for a time step into the following: ! netMassInOut = surface water fluxes [H ~> m or kg m-2] over time step diff --git a/src/parameterizations/vertical/MOM_vert_friction.F90 b/src/parameterizations/vertical/MOM_vert_friction.F90 index c6a6f37b39..1a4fb58e02 100644 --- a/src/parameterizations/vertical/MOM_vert_friction.F90 +++ b/src/parameterizations/vertical/MOM_vert_friction.F90 @@ -1143,12 +1143,12 @@ subroutine find_coupling_coef(a_cpl, hvel, do_i, h_harm, bbl_thick, kv_bbl, z_i, real :: h_neglect ! A thickness that is so small it is usually lost ! in roundoff and can be neglected [H ~> m or kg m-2]. real :: z2 ! A copy of z_i [nondim] + real :: botfn ! A function that is 1 at the bottom and small far from it [nondim] real :: topfn ! A function that is 1 at the top and small far from it [nondim] real :: kv_top ! A viscosity associated with the top boundary layer [Z2 T-1 ~> m2 s-1] logical :: do_shelf, do_OBCs integer :: i, k, is, ie, max_nk integer :: nz - real :: botfn a_cpl(:,:) = 0.0 Kv_tot(:,:) = 0.0 From feed9ba82b9ebe6805ce691e77190b2f5ba4f7ee Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Fri, 7 Aug 2020 08:13:40 -0400 Subject: [PATCH 6/7] (*)Fix an indexing bug in int_density_dz_linear Corrected a horizontal indexing bug in int_density_dz_linear that caused the ISOMIP/layer test case to fail. This bug was first introduced with PR#732 on March 8, 2018. This bug fix will change answers with a linear equation of state and the finite volume pressure gradient force, however it does not change any of the verified answers in the MOM6-examples regression suite. --- src/equation_of_state/MOM_EOS_linear.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/equation_of_state/MOM_EOS_linear.F90 b/src/equation_of_state/MOM_EOS_linear.F90 index e3a5443840..47a2bf21b0 100644 --- a/src/equation_of_state/MOM_EOS_linear.F90 +++ b/src/equation_of_state/MOM_EOS_linear.F90 @@ -473,7 +473,7 @@ subroutine int_density_dz_linear(T, S, z_t, z_b, rho_ref, rho_0_pres, G_e, HI, & hWt_LL = (hWght*hL + hR*hL) * iDenom ; hWt_LR = (hWght*hR) * iDenom hWt_RR = (hWght*hR + hR*hL) * iDenom ; hWt_RL = (hWght*hL) * iDenom - intz(1) = dpa(i,j) ; intz(5) = dpa(i+1,j) + intz(1) = dpa(i,j) ; intz(5) = dpa(i,j+1) do m=2,4 wt_L = 0.25*real(5-m) ; wt_R = 1.0-wt_L wtT_L = wt_L*hWt_LL + wt_R*hWt_RL ; wtT_R = wt_L*hWt_LR + wt_R*hWt_RR From 3c10ae18a72b3096ea69b81dc3906931eefa9a6f Mon Sep 17 00:00:00 2001 From: Matthew Harrison Date: Mon, 10 Aug 2020 10:46:31 -0400 Subject: [PATCH 7/7] Remove outdated comments --- src/core/MOM_open_boundary.F90 | 7 ------- 1 file changed, 7 deletions(-) diff --git a/src/core/MOM_open_boundary.F90 b/src/core/MOM_open_boundary.F90 index f94060fc39..37ebeda1fa 100644 --- a/src/core/MOM_open_boundary.F90 +++ b/src/core/MOM_open_boundary.F90 @@ -546,13 +546,6 @@ subroutine open_boundary_config(G, US, param_file, OBC) if (Lscale_out>0.) OBC%segment(l)%Tr_InvLscale_out = 1.0/Lscale_out enddo - ! There is a problem with the order of the OBC initialization - ! with respect to ALE_init. Currently handling this by copying the - ! param file so that I can use it later in step_MOM in order to finish - ! initializing segments on the first step. - - ! Is the above comment still relevant ? - call get_param(param_file, mdl, "REMAPPING_SCHEME", remappingScheme, & "This sets the reconstruction scheme used "//& "for vertical remapping for all variables. "//&