Merge pull request #43 from NOAA-GFDL/dev/gfdl

Merge in latest dev/gfdl updates

src/core/MOM_continuity_PPM.F90

@@ -1331,7 +1331,7 @@ subroutine merid_flux_layer(v, h, h_L, h_R, vh, dvhdv, visc_rem, dt, G, US, J, &
 
   local_open_BC = .false.
   if (present(OBC)) then ; if (associated(OBC)) then
-    local_open_BC = OBC%open_u_BCs_exist_globally
+    local_open_BC = OBC%open_v_BCs_exist_globally
   endif ; endif
 
   do i=ish,ieh ; if (do_I(i)) then

src/core/MOM_open_boundary.F90

@@ -118,6 +118,7 @@ module MOM_open_boundary
   logical :: nudged_grad    !< Optional supplement to nudge normal gradient of tangential velocity.
   logical :: specified      !< Boundary normal velocity fixed to external value.
   logical :: specified_tan  !< Boundary tangential velocity fixed to external value.
+  logical :: specified_grad !< Boundary gradient of tangential velocity fixed to external value.
   logical :: open           !< Boundary is open for continuity solver.
   logical :: gradient       !< Zero gradient at boundary.
   logical :: values_needed  !< Whether or not any external OBC fields are needed.
@@ -436,6 +437,7 @@ subroutine open_boundary_config(G, US, param_file, OBC)
       OBC%segment(l)%nudged_grad = .false.
       OBC%segment(l)%specified = .false.
       OBC%segment(l)%specified_tan = .false.
+      OBC%segment(l)%specified_grad = .false.
       OBC%segment(l)%open = .false.
       OBC%segment(l)%gradient = .false.
       OBC%segment(l)%values_needed = .false.
@@ -941,6 +943,10 @@ subroutine setup_u_point_obc(OBC, G, US, segment_str, l_seg, PF, reentrant_y)
       OBC%segment%u_values_needed = .true.
     elseif (trim(action_str(a_loop)) == 'SIMPLE_TAN') then
       OBC%segment(l_seg)%specified_tan = .true.
+      OBC%segment%v_values_needed = .true.
+    elseif (trim(action_str(a_loop)) == 'SIMPLE_GRAD') then
+      OBC%segment(l_seg)%specified_grad = .true.
+      OBC%segment%g_values_needed = .true.
     else
       call MOM_error(FATAL, "MOM_open_boundary.F90, setup_u_point_obc: "//&
                      "String '"//trim(action_str(a_loop))//"' not understood.")
@@ -1078,6 +1084,10 @@ subroutine setup_v_point_obc(OBC, G, US, segment_str, l_seg, PF, reentrant_x)
       OBC%segment%v_values_needed = .true.
     elseif (trim(action_str(a_loop)) == 'SIMPLE_TAN') then
       OBC%segment(l_seg)%specified_tan = .true.
+      OBC%segment%u_values_needed = .true.
+    elseif (trim(action_str(a_loop)) == 'SIMPLE_GRAD') then
+      OBC%segment(l_seg)%specified_grad = .true.
+      OBC%segment%g_values_needed = .true.
     else
       call MOM_error(FATAL, "MOM_open_boundary.F90, setup_v_point_obc: "//&
                      "String '"//trim(action_str(a_loop))//"' not understood.")
@@ -2774,7 +2784,7 @@ subroutine radiation_open_bdry_conds(OBC, u_new, u_old, v_new, v_old, G, US, dt)
            enddo ; enddo
          endif
          if (segment%nudged_grad) then
-           do k=1,nz ; do J=segment%HI%JsdB,segment%HI%JedB
+           do k=1,nz ; do I=segment%HI%IsdB,segment%HI%IedB
              ! dhdt gets set to 0 on inflow in oblique case
              if (ry_tang_rad(I,J,k) <= 0.0) then
                tau = segment%Velocity_nudging_timescale_in
@@ -3211,7 +3221,7 @@ subroutine allocate_OBC_segment_data(OBC, segment)
       allocate(segment%nudged_tangential_grad(IsdB:IedB,JsdB:JedB,OBC%ke)); segment%nudged_tangential_grad(:,:,:)=0.0
     endif
     if (OBC%specified_vorticity .or. OBC%specified_strain .or. segment%radiation_grad .or. &
-              segment%oblique_grad) then
+              segment%oblique_grad .or. segment%specified_grad) then
       allocate(segment%tangential_grad(IsdB:IedB,JsdB:JedB,OBC%ke)); segment%tangential_grad(:,:,:)=0.0
     endif
     if (segment%oblique) then
@@ -3254,7 +3264,7 @@ subroutine allocate_OBC_segment_data(OBC, segment)
       allocate(segment%nudged_tangential_grad(IsdB:IedB,JsdB:JedB,OBC%ke)); segment%nudged_tangential_grad(:,:,:)=0.0
     endif
     if (OBC%specified_vorticity .or. OBC%specified_strain .or. segment%radiation_grad .or. &
-              segment%oblique_grad) then
+              segment%oblique_grad .or. segment%specified_grad) then
       allocate(segment%tangential_grad(IsdB:IedB,JsdB:JedB,OBC%ke)); segment%tangential_grad(:,:,:)=0.0
     endif
     if (segment%oblique) then
@@ -3761,8 +3771,9 @@ subroutine update_OBC_segment_data(G, GV, US, OBC, tv, h, Time)
         endif
       endif
 
-      if (trim(segment%field(m)%name) == 'U' .or. trim(segment%field(m)%name) == 'V') then
-        if (segment%field(m)%fid>0) then ! calculate external BT velocity and transport if needed
+      if (segment%field(m)%fid>0) then
+        ! calculate external BT velocity and transport if needed
+        if (trim(segment%field(m)%name) == 'U' .or. trim(segment%field(m)%name) == 'V') then
           if (trim(segment%field(m)%name) == 'U' .and. segment%is_E_or_W) then
             I=is_obc
             do j=js_obc+1,je_obc
@@ -3809,23 +3820,27 @@ subroutine update_OBC_segment_data(G, GV, US, OBC, tv, h, Time)
               if (associated(segment%nudged_tangential_vel)) &
                 segment%nudged_tangential_vel(I,J,:) = segment%tangential_vel(I,J,:)
             enddo
-          elseif (trim(segment%field(m)%name) == 'DVDX' .and. segment%is_E_or_W .and. &
-                  associated(segment%tangential_grad)) then
-            I=is_obc
-            do J=js_obc,je_obc
-              do k=1,G%ke
-                segment%tangential_grad(I,J,k) = US%T_to_s*segment%field(m)%buffer_dst(I,J,k)
-              enddo
+          endif
+        elseif (trim(segment%field(m)%name) == 'DVDX' .and. segment%is_E_or_W .and. &
+                associated(segment%tangential_grad)) then
+          I=is_obc
+          do J=js_obc,je_obc
+            do k=1,G%ke
+              segment%tangential_grad(I,J,k) = US%T_to_s*segment%field(m)%buffer_dst(I,J,k)
+              if (associated(segment%nudged_tangential_grad)) &
+                segment%nudged_tangential_grad(I,J,:) = segment%tangential_grad(I,J,:)
             enddo
-          elseif (trim(segment%field(m)%name) == 'DUDY' .and. segment%is_N_or_S .and. &
-                  associated(segment%tangential_grad)) then
-            J=js_obc
-            do I=is_obc,ie_obc
-              do k=1,G%ke
-                segment%tangential_grad(I,J,k) = US%T_to_s*segment%field(m)%buffer_dst(I,J,k)
-              enddo
+          enddo
+        elseif (trim(segment%field(m)%name) == 'DUDY' .and. segment%is_N_or_S .and. &
+                associated(segment%tangential_grad)) then
+          J=js_obc
+          do I=is_obc,ie_obc
+            do k=1,G%ke
+              segment%tangential_grad(I,J,k) = US%T_to_s*segment%field(m)%buffer_dst(I,J,k)
+              if (associated(segment%nudged_tangential_grad)) &
+                segment%nudged_tangential_grad(I,J,:) = segment%tangential_grad(I,J,:)
             enddo
-          endif
+          enddo
         endif
       endif
 

src/user/dumbbell_initialization.F90

@@ -49,27 +49,43 @@ subroutine dumbbell_initialize_topography( D, G, param_file, max_depth )
   ! Local variables
   integer   :: i, j
   real      :: x, y, delta, dblen, dbfrac
+  logical   :: dbrotate
 
   call get_param(param_file, mdl,"DUMBBELL_LEN",dblen, &
                 'Lateral Length scale for dumbbell.',&
                  units='k', default=600., do_not_log=.false.)
   call get_param(param_file, mdl,"DUMBBELL_FRACTION",dbfrac, &
                 'Meridional fraction for narrow part of dumbbell.',&
                  units='nondim', default=0.5, do_not_log=.false.)
+  call get_param(param_file, mdl, "DUMBBELL_ROTATION", dbrotate, &
+                'Logical for rotation of dumbbell domain.',&
+                 units='nondim', default=.false., do_not_log=.false.)
 
   if (G%x_axis_units == 'm') then
     dblen=dblen*1.e3
   endif
 
-  do j=G%jsc,G%jec ; do i=G%isc,G%iec
-    ! Compute normalized zonal coordinates (x,y=0 at center of domain)
-    x = ( G%geoLonT(i,j) ) / dblen
-    y = ( G%geoLatT(i,j)  ) / G%len_lat
-    D(i,j) = G%max_depth
-    if ((x>=-0.25 .and. x<=0.25) .and. (y <= -0.5*dbfrac .or. y >= 0.5*dbfrac)) then
-      D(i,j) = 0.0
-    endif
-  enddo ; enddo
+  if (dbrotate) then
+    do j=G%jsc,G%jec ; do i=G%isc,G%iec
+      ! Compute normalized zonal coordinates (x,y=0 at center of domain)
+      x = ( G%geoLonT(i,j) ) / G%len_lon
+      y = ( G%geoLatT(i,j)  ) / dblen
+      D(i,j) = G%max_depth
+      if ((y>=-0.25 .and. y<=0.25) .and. (x <= -0.5*dbfrac .or. x >= 0.5*dbfrac)) then
+        D(i,j) = 0.0
+      endif
+    enddo ; enddo
+  else
+    do j=G%jsc,G%jec ; do i=G%isc,G%iec
+      ! Compute normalized zonal coordinates (x,y=0 at center of domain)
+      x = ( G%geoLonT(i,j) ) / dblen
+      y = ( G%geoLatT(i,j)  ) / G%len_lat
+      D(i,j) = G%max_depth
+      if ((x>=-0.25 .and. x<=0.25) .and. (y <= -0.5*dbfrac .or. y >= 0.5*dbfrac)) then
+        D(i,j) = 0.0
+      endif
+    enddo ; enddo
+  endif
 
 end subroutine dumbbell_initialize_topography
 
@@ -209,6 +225,7 @@ subroutine dumbbell_initialize_temperature_salinity ( T, S, h, G, GV, param_file
   real    :: x, y, dblen
   real    :: T_ref, T_Light, T_Dense, S_ref, S_Light, S_Dense, a1, frac_dense, k_frac, res_rat
   logical :: just_read    ! If true, just read parameters but set nothing.
+  logical :: dbrotate     ! If true, rotate the domain.
   character(len=20) :: verticalCoordinate, density_profile
 
   is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke
@@ -230,6 +247,9 @@ subroutine dumbbell_initialize_temperature_salinity ( T, S, h, G, GV, param_file
   call get_param(param_file, mdl,"DUMBBELL_LEN",dblen, &
                 'Lateral Length scale for dumbbell ',&
                  units='k', default=600., do_not_log=just_read)
+  call get_param(param_file, mdl, "DUMBBELL_ROTATION", dbrotate, &
+                'Logical for rotation of dumbbell domain.',&
+                 units='nondim', default=.false., do_not_log=just_read)
 
   if (G%x_axis_units == 'm') then
     dblen=dblen*1.e3
@@ -238,7 +258,12 @@ subroutine dumbbell_initialize_temperature_salinity ( T, S, h, G, GV, param_file
   do j=G%jsc,G%jec
     do i=G%isc,G%iec
     ! Compute normalized zonal coordinates (x,y=0 at center of domain)
-      x = ( G%geoLonT(i,j) ) / dblen
+      if (dbrotate) then
+        ! This is really y in the rotated case
+        x = ( G%geoLatT(i,j) ) / dblen
+      else
+        x = ( G%geoLonT(i,j) ) / dblen
+      endif
       do k=1,nz
         T(i,j,k)=T_surf
       enddo
@@ -278,9 +303,13 @@ subroutine dumbbell_initialize_sponges(G, GV, US, tv, param_file, use_ALE, CSp,
   integer :: i, j, k, nz
   real :: x, zi, zmid, dist, min_thickness, dblen
   real :: mld, S_ref, S_range, S_dense, T_ref, sill_height
+  logical :: dbrotate    ! If true, rotate the domain.
   call get_param(param_file, mdl,"DUMBBELL_LEN",dblen, &
                 'Lateral Length scale for dumbbell ',&
                  units='k', default=600., do_not_log=.true.)
+  call get_param(param_file, mdl, "DUMBBELL_ROTATION", dbrotate, &
+                'Logical for rotation of dumbbell domain.',&
+                 units='nondim', default=.false., do_not_log=.true.)
 
   if (G%x_axis_units == 'm') then
     dblen=dblen*1.e3
@@ -307,7 +336,12 @@ subroutine dumbbell_initialize_sponges(G, GV, US, tv, param_file, use_ALE, CSp,
     do i = G%isc,G%iec
       if (G%mask2dT(i,j) > 0.) then
         ! nondimensional x position
-        x = (G%geoLonT(i,j) ) / dblen
+        if (dbrotate) then
+          ! This is really y in the rotated case
+          x = ( G%geoLatT(i,j) ) / dblen
+        else
+          x = ( G%geoLonT(i,j) ) / dblen
+        endif
         if (x > 0.25 .or. x < -0.25) then
           ! scale restoring by depth into sponge
           Idamp(i,j) = 1. / sponge_time_scale
@@ -339,18 +373,23 @@ subroutine dumbbell_initialize_sponges(G, GV, US, tv, param_file, use_ALE, CSp,
 
     do j=G%jsc,G%jec ; do i=G%isc,G%iec
       ! Compute normalized zonal coordinates (x,y=0 at center of domain)
-       x = ( G%geoLonT(i,j) ) / dblen
-       if (x>=0.25 ) then
-         do k=1,nz
-           S(i,j,k)=S_ref + 0.5*S_range
-         enddo
-       endif
-       if (x<=-0.25 ) then
-         do k=1,nz
-           S(i,j,k)=S_ref - 0.5*S_range
-         enddo
-       endif
-     enddo ; enddo
+      if (dbrotate) then
+        ! This is really y in the rotated case
+        x = ( G%geoLatT(i,j) ) / dblen
+      else
+        x = ( G%geoLonT(i,j) ) / dblen
+      endif
+      if (x>=0.25 ) then
+        do k=1,nz
+          S(i,j,k)=S_ref + 0.5*S_range
+        enddo
+      endif
+      if (x<=-0.25 ) then
+        do k=1,nz
+          S(i,j,k)=S_ref - 0.5*S_range
+        enddo
+      endif
+    enddo ; enddo
   endif
 
   if (associated(tv%S)) call set_up_ALE_sponge_field(S, G, tv%S, ACSp)

src/user/dumbbell_surface_forcing.F90

@@ -189,6 +189,7 @@ subroutine dumbbell_surface_forcing_init(Time, G, US, param_file, diag, CS)
   real :: S_surf, S_range
   real :: x, y
   integer :: i, j
+  logical :: dbrotate    ! If true, rotate the domain.
 #include "version_variable.h"
   character(len=40)  :: mdl = "dumbbell_surface_forcing" ! This module's name.
 
@@ -224,6 +225,9 @@ subroutine dumbbell_surface_forcing_init(Time, G, US, param_file, diag, CS)
   call get_param(param_file, mdl, "DUMBBELL_SLP_PERIOD", CS%slp_period, &
                  "Periodicity of SLP forcing in reservoirs.", &
                  units="days", default = 1.0)
+  call get_param(param_file, mdl, "DUMBBELL_ROTATION", dbrotate, &
+                'Logical for rotation of dumbbell domain.',&
+                 units='nondim', default=.false., do_not_log=.true.)
   call get_param(param_file, mdl,"INITIAL_SSS", S_surf, &
                  "Initial surface salinity", units="1e-3", default=34.0, do_not_log=.true.)
   call get_param(param_file, mdl,"INITIAL_S_RANGE", S_range, &
@@ -250,8 +254,14 @@ subroutine dumbbell_surface_forcing_init(Time, G, US, param_file, diag, CS)
     do j=G%jsc,G%jec
       do i=G%isc,G%iec
         ! Compute normalized zonal coordinates (x,y=0 at center of domain)
-        x = ( G%geoLonT(i,j) - G%west_lon ) / G%len_lon - 0.5
-        y = ( G%geoLatT(i,j) - G%south_lat ) / G%len_lat - 0.5
+!       x = ( G%geoLonT(i,j) - G%west_lon ) / G%len_lon - 0.5
+!       y = ( G%geoLatT(i,j) - G%south_lat ) / G%len_lat - 0.5
+        if (dbrotate) then
+          ! This is really y in the rotated case
+          x = ( G%geoLatT(i,j) - G%south_lat ) / G%len_lat - 0.5
+        else
+          x = ( G%geoLonT(i,j) - G%west_lon ) / G%len_lon - 0.5
+        endif
         CS%forcing_mask(i,j)=0
         CS%S_restore(i,j) = S_surf
         if ((x>0.25)) then