Merge b950b05 into 1ae4e16

src/core/MOM.F90

@@ -2652,11 +2652,6 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, &
                             param_file, diag, CS%diagnostics_CSp, CS%tv)
   call diag_copy_diag_to_storage(CS%diag_pre_sync, CS%h, CS%diag)
 
-  if (associated(CS%sponge_CSp)) &
-    call init_sponge_diags(Time, G, GV, US, diag, CS%sponge_CSp)
-
-  if (associated(CS%ALE_sponge_CSp)) &
-    call init_ALE_sponge_diags(Time, G, diag, CS%ALE_sponge_CSp)
 
   if (CS%adiabatic) then
     call adiabatic_driver_init(Time, G, param_file, diag, CS%diabatic_CSp, &
@@ -2667,6 +2662,12 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, &
                               CS%sponge_CSp, CS%ALE_sponge_CSp)
   endif
 
+  if (associated(CS%sponge_CSp)) &
+    call init_sponge_diags(Time, G, GV, US, diag, CS%sponge_CSp)
+
+  if (associated(CS%ALE_sponge_CSp)) &
+    call init_ALE_sponge_diags(Time, G, diag, CS%ALE_sponge_CSp, US)
+
   call tracer_advect_init(Time, G, US, param_file, diag, CS%tracer_adv_CSp)
   call tracer_hor_diff_init(Time, G, GV, US, param_file, diag, CS%tv%eqn_of_state, CS%diabatic_CSp, &
                             CS%tracer_diff_CSp)

src/framework/MOM_horizontal_regridding.F90

@@ -667,6 +667,7 @@ subroutine horiz_interp_and_extrap_tracer_fms_id(fms_id,  Time, conversion, G, t
   character(len=12)  :: dim_name(4)
   logical :: debug=.false.
   logical :: spongeDataOngrid
+  logical :: ans_2018
   real :: npoints, varAvg
   real, dimension(SZI_(G),SZJ_(G)) :: lon_out, lat_out ! The longitude and latitude of points on the model grid
   real, dimension(SZI_(G),SZJ_(G)) :: tr_out, mask_out ! The tracer and mask on the model grid
@@ -688,6 +689,8 @@ subroutine horiz_interp_and_extrap_tracer_fms_id(fms_id,  Time, conversion, G, t
 
   PI_180 = atan(1.0)/45.
 
+  ans_2018 = .true.;if (present(answers_2018)) ans_2018 = answers_2018
+
   ! Open NetCDF file and if present, extract data and spatial coordinate information
   ! The convention adopted here requires that the data be written in (i,j,k) ordering.
 
@@ -886,15 +889,16 @@ subroutine horiz_interp_and_extrap_tracer_fms_id(fms_id,  Time, conversion, G, t
 
     enddo ! kd
   else
-    call time_interp_external(fms_id, Time, data_in, verbose=.true., turns=turns)
-    do k=1,kd
-      do j=js,je
-        do i=is,ie
-          tr_z(i,j,k) = data_in(i,j,k)
-          if (abs(tr_z(i,j,k)-missing_value) < abs(roundoff*missing_value)) mask_z(i,j,k) = 0.
+      call time_interp_external(fms_id, Time, data_in, verbose=.true., turns=turns)
+      do k=1,kd
+        do j=js,je
+          do i=is,ie
+            tr_z(i,j,k)=data_in(i,j,k)
+            if (.not. ans_2018) mask_z(i,j,k) = 1.
+            if (abs(tr_z(i,j,k)-missing_value) < abs(roundoff*missing_value)) mask_z(i,j,k) = 0.
+          enddo
         enddo
       enddo
-    enddo
   endif
 
 end subroutine horiz_interp_and_extrap_tracer_fms_id

src/initialization/MOM_state_initialization.F90

@@ -31,7 +31,8 @@ module MOM_state_initialization
 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, ALE_sponge_CS
+use MOM_ALE_sponge, only : set_up_ALE_sponge_field, set_up_ALE_sponge_vel_field
+use MOM_ALE_sponge, only : ALE_sponge_CS, initialize_ALE_sponge
 use MOM_string_functions, only : uppercase, lowercase
 use MOM_time_manager, only : time_type
 use MOM_tracer_registry, only : tracer_registry_type
@@ -549,7 +550,7 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, US, PF, dirs, &
       case ("phillips"); call Phillips_initialize_sponges(G, GV, US, tv, PF, sponge_CSp, h)
       case ("dense"); call dense_water_initialize_sponges(G, GV, US, tv, PF, useALE, &
                                                           sponge_CSp, ALE_sponge_CSp)
-      case ("file"); call initialize_sponges_file(G, GV, US, use_temperature, tv, PF, &
+      case ("file"); call initialize_sponges_file(G, GV, US, use_temperature, tv, u, v, PF, &
                                                   sponge_CSp, ALE_sponge_CSp, Time)
       case default ; call MOM_error(FATAL,  "MOM_initialize_state: "//&
              "Unrecognized sponge configuration "//trim(config))
@@ -1718,13 +1719,19 @@ end subroutine initialize_temp_salt_linear
 !! number of tracers should be restored within each sponge. The
 !! interface height is always subject to damping, and must always be
 !! the first registered field.
-subroutine initialize_sponges_file(G, GV, US, use_temperature, tv, param_file, Layer_CSp, ALE_CSp, Time)
+subroutine initialize_sponges_file(G, GV, US, use_temperature, tv, u, v, param_file, Layer_CSp, ALE_CSp, Time)
   type(ocean_grid_type),   intent(in) :: G    !< The ocean's grid structure.
   type(verticalGrid_type), intent(in) :: GV   !< The ocean's vertical grid structure.
   type(unit_scale_type),   intent(in) :: US  !< A dimensional unit scaling type
   logical,                 intent(in) :: use_temperature !< If true, T & S are state variables.
   type(thermo_var_ptrs),   intent(in) :: tv   !< A structure pointing to various thermodynamic
                                               !! variables.
+  real, target, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
+                              intent(in)   :: u    !< The zonal velocity that is being
+                                                    !! initialized [L T-1 ~> m s-1]
+  real, target, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
+                              intent(in)    :: v    !< The meridional velocity that is being
+                                                    !! initialized [L T-1 ~> m s-1]
   type(param_file_type),   intent(in) :: param_file !< A structure to parse for run-time parameters.
   type(sponge_CS),         pointer    :: Layer_CSp  !< A pointer that is set to point to the control
                                               !! structure for this module (in layered mode).
@@ -1741,18 +1748,22 @@ subroutine initialize_sponges_file(G, GV, US, use_temperature, tv, param_file, L
   real, dimension (SZI_(G),SZJ_(G)) :: &
     tmp_2d ! A temporary array for tracers.
   real, allocatable, dimension(:,:,:) :: tmp_tr ! A temporary array for reading sponge fields
+  real, allocatable, dimension(:,:,:) :: tmp_u,tmp_v ! A temporary array for reading sponge fields
 
   real :: Idamp(SZI_(G),SZJ_(G))    ! The inverse damping rate [T-1 ~> s-1].
+  real :: Idamp_u(SZIB_(G),SZJ_(G)) ! The inverse damping rate for velocity fields [T-1 ~> s-1].
+  real :: Idamp_v(SZI_(G),SZJB_(G)) ! The inverse damping rate for velocity fields [T-1 ~> s-1].
   real :: pres(SZI_(G))     ! An array of the reference pressure [R L2 T-2 ~> Pa].
 
   integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state
   integer :: i, j, k, is, ie, js, je, nz
   integer :: isd, ied, jsd, jed
   integer, dimension(4) :: siz
   integer :: nz_data  ! The size of the sponge source grid
-  character(len=40) :: potemp_var, salin_var, Idamp_var, eta_var
+  logical :: sponge_uv ! Apply sponges in u and v, in addition to tracers.
+  character(len=40) :: potemp_var, salin_var, u_var, v_var, Idamp_var, Idamp_u_var, Idamp_v_var, eta_var
   character(len=40) :: mdl = "initialize_sponges_file"
-  character(len=200) :: damping_file, state_file  ! Strings for filenames
+  character(len=200) :: damping_file, uv_damping_file, state_file, state_uv_file  ! Strings for filenames
   character(len=200) :: filename, inputdir ! Strings for file/path and path.
 
   logical :: use_ALE ! True if ALE is being used, False if in layered mode
@@ -1764,7 +1775,7 @@ subroutine initialize_sponges_file(G, GV, US, use_temperature, tv, param_file, L
   is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
   isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed
 
-  pres(:) = 0.0 ; tmp(:,:,:) = 0.0 ; Idamp(:,:) = 0.0
+  pres(:) = 0.0 ; tmp(:,:,:) = 0.0 ; Idamp(:,:) = 0.0 ; Idamp_u(:,:) = 0.0 ; Idamp_v(:,:) = 0.0
 
   call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".")
   inputdir = slasher(inputdir)
@@ -1774,18 +1785,43 @@ subroutine initialize_sponges_file(G, GV, US, use_temperature, tv, param_file, L
   call get_param(param_file, mdl, "SPONGE_STATE_FILE", state_file, &
                  "The name of the file with the state to damp toward.", &
                  default=damping_file)
+  call get_param(param_file, mdl, "SPONGE_UV_STATE_FILE", state_uv_file, &
+                 "The name of the file with the state to damp UV toward.", &
+                 default=damping_file)
   call get_param(param_file, mdl, "SPONGE_PTEMP_VAR", potemp_var, &
                  "The name of the potential temperature variable in "//&
                  "SPONGE_STATE_FILE.", default="PTEMP")
   call get_param(param_file, mdl, "SPONGE_SALT_VAR", salin_var, &
                  "The name of the salinity variable in "//&
                  "SPONGE_STATE_FILE.", default="SALT")
+  call get_param(param_file, mdl, "SPONGE_UV", sponge_uv, &
+                 "Apply sponges in u and v, in addition to tracers.", &
+                 default=.false.)
+  if (sponge_uv) then
+    call get_param(param_file, mdl, "SPONGE_U_VAR", u_var, &
+                 "The name of the zonal velocity variable in "//&
+                 "SPONGE_UV_STATE_FILE.", default="UVEL")
+    call get_param(param_file, mdl, "SPONGE_V_VAR", v_var, &
+                 "The name of the vertical velocity variable in "//&
+                 "SPONGE_UV_STATE_FILE.", default="VVEL")
+  endif
   call get_param(param_file, mdl, "SPONGE_ETA_VAR", eta_var, &
                  "The name of the interface height variable in "//&
                  "SPONGE_STATE_FILE.", default="ETA")
   call get_param(param_file, mdl, "SPONGE_IDAMP_VAR", Idamp_var, &
                  "The name of the inverse damping rate variable in "//&
-                 "SPONGE_DAMPING_FILE.", default="Idamp")
+                 "SPONGE_DAMPING_FILE.", default="IDAMP")
+  if (sponge_uv) then
+    call get_param(param_file, mdl, "SPONGE_UV_DAMPING_FILE", uv_damping_file, &
+                   "The name of the file with sponge damping rates for the velocity variables.", &
+                   default=damping_file)
+    call get_param(param_file, mdl, "SPONGE_IDAMP_U_var", Idamp_u_var, &
+                   "The name of the inverse damping rate variable in "//&
+                   "SPONGE_UV_DAMPING_FILE for the velocities.", default=Idamp_var)
+    call get_param(param_file, mdl, "SPONGE_IDAMP_V_var", Idamp_v_var, &
+                   "The name of the inverse damping rate variable in "//&
+                   "SPONGE_UV_DAMPING_FILE for the velocities.", default=Idamp_var)
+  endif
   call get_param(param_file, mdl, "USE_REGRIDDING", use_ALE, do_not_log = .true.)
   time_space_interp_sponge = .false.
   call get_param(param_file, mdl, "NEW_SPONGES", time_space_interp_sponge, &
@@ -1802,6 +1838,8 @@ subroutine initialize_sponges_file(G, GV, US, use_temperature, tv, param_file, L
                  "performs on-the-fly regridding in lat-lon-time.",&
                  "of sponge restoring data.", default=time_space_interp_sponge)
 
+
+  ! Read in inverse damping rate for tracers
   filename = trim(inputdir)//trim(damping_file)
   call log_param(param_file, mdl, "INPUTDIR/SPONGE_DAMPING_FILE", filename)
   if (.not.file_exists(filename, G%Domain)) &
@@ -1812,6 +1850,32 @@ subroutine initialize_sponges_file(G, GV, US, use_temperature, tv, param_file, L
 
   call MOM_read_data(filename, Idamp_var, Idamp(:,:), G%Domain, scale=US%T_to_s)
 
+  ! Read in inverse damping rate for velocities
+  if (sponge_uv) then
+    if (separate_idamp_for_uv()) then
+      filename = trim(inputdir)//trim(uv_damping_file)
+      call log_param(param_file, mdl, "INPUTDIR/SPONGE_UV_DAMPING_FILE", filename)
+
+      if (.not.file_exists(filename, G%Domain)) &
+        call MOM_error(FATAL, " initialize_sponges: Unable to open "//trim(filename))
+
+      call MOM_read_vector(filename, Idamp_u_var,Idamp_v_var,Idamp_u(:,:),Idamp_v(:,:), G%Domain, scale=US%T_to_s)
+    else
+       !      call MOM_error(FATAL, "Must provide SPONGE_IDAMP_U_var and SPONGE_IDAMP_V_var")
+       call pass_var(Idamp,G%Domain)
+       do j=G%jsc,G%jec
+         do i=G%iscB,G%iecB
+           Idamp_u(I,j) = 0.5*(Idamp(i,j)+Idamp(i+1,j))
+         enddo
+       enddo
+       do j=G%jscB,G%jecB
+         do i=G%isc,G%iec
+           Idamp_v(i,J) = 0.5*(Idamp(i,j)+Idamp(i,j+1))
+         enddo
+       enddo
+    endif
+  endif
+
   ! Now register all of the fields which are damped in the sponge.
   ! By default, momentum is advected vertically within the sponge, but
   ! momentum is typically not damped within the sponge.
@@ -1869,11 +1933,18 @@ subroutine initialize_sponges_file(G, GV, US, use_temperature, tv, param_file, L
       call set_up_sponge_field(tmp, tv%S, G, GV, nz, Layer_CSp)
     endif
 
+!  else
+    ! Initialize sponges without supplying sponge grid
+!    if (sponge_uv) then
+!      call initialize_ALE_sponge(Idamp, G, GV, param_file, ALE_CSp, Idamp_u, Idamp_v)
+!    else
+!      call initialize_ALE_sponge(Idamp, G, GV, param_file, ALE_CSp)
+!    endif
   endif
 
 
-  if  (use_ALE) then
-    if (.not. time_space_interp_sponge) then ! ALE mode
+  if  (use_ALE) then ! ALE mode
+    if (.not. time_space_interp_sponge) then
       call field_size(filename,eta_var,siz,no_domain=.true.)
       if (siz(1) /= G%ieg-G%isg+1 .or. siz(2) /= G%jeg-G%jsg+1) &
         call MOM_error(FATAL,"initialize_sponge_file: Array size mismatch for sponge data.")
@@ -1890,8 +1961,12 @@ subroutine initialize_sponges_file(G, GV, US, use_temperature, tv, param_file, L
       enddo ; enddo ; enddo
       do k=1,nz ; do j=js,je ; do i=is,ie
         h(i,j,k) = GV%Z_to_H*(eta(i,j,k)-eta(i,j,k+1))
-      enddo ; enddo ; enddo
-      call initialize_ALE_sponge(Idamp, G, GV, param_file, ALE_CSp, h, nz_data)
+      enddo; enddo ; enddo
+      if (sponge_uv) then
+        call initialize_ALE_sponge(Idamp, G, GV, param_file, ALE_CSp, h, nz_data, Idamp_u, Idamp_v)
+      else
+        call initialize_ALE_sponge(Idamp, G, GV, param_file, ALE_CSp, h, nz_data)
+      endif
       deallocate(eta)
       deallocate(h)
       if (use_temperature) then
@@ -1902,20 +1977,50 @@ subroutine initialize_sponges_file(G, GV, US, use_temperature, tv, param_file, L
         call set_up_ALE_sponge_field(tmp_tr, G, GV, tv%S, ALE_CSp)
         deallocate(tmp_tr)
       endif
+      if (sponge_uv) then
+        filename = trim(inputdir)//trim(state_uv_file)
+        call log_param(param_file, mdl, "INPUTDIR/SPONGE_STATE_UV_FILE", filename)
+        if (.not.file_exists(filename, G%Domain)) &
+             call MOM_error(FATAL, " initialize_sponges: Unable to open "//trim(filename))
+        allocate(tmp_u(G%IsdB:G%IedB,jsd:jed,nz_data))
+        allocate(tmp_v(isd:ied,G%JsdB:G%JedB,nz_data))
+        call MOM_read_vector(filename, u_var, v_var, tmp_u(:,:,:), tmp_v(:,:,:), G%Domain,scale=US%m_s_to_L_T)
+        call set_up_ALE_sponge_vel_field(tmp_u, tmp_v, G, GV, u, v, ALE_CSp)
+        deallocate(tmp_u,tmp_v)
+      endif
     else
       ! Initialize sponges without supplying sponge grid
-      call initialize_ALE_sponge(Idamp, G, GV, param_file, ALE_CSp)
+      if (sponge_uv) then
+        call initialize_ALE_sponge(Idamp, G, GV, param_file, ALE_CSp, Idamp_u, Idamp_v)
+      else
+        call initialize_ALE_sponge(Idamp, G, GV, param_file, ALE_CSp)
+      endif
       ! The remaining calls to set_up_sponge_field can be in any order.
       if ( use_temperature) then
         call set_up_ALE_sponge_field(filename, potemp_var, Time, G, GV, US, tv%T, ALE_CSp)
         call set_up_ALE_sponge_field(filename, salin_var, Time, G, GV, US, tv%S, ALE_CSp)
       endif
+      if (sponge_uv) then
+        filename = trim(inputdir)//trim(state_uv_file)
+        call log_param(param_file, mdl, "INPUTDIR/SPONGE_STATE_UV_FILE", filename)
+        if (.not.file_exists(filename, G%Domain)) &
+             call MOM_error(FATAL, " initialize_sponges: Unable to open "//trim(filename))
+        call set_up_ALE_sponge_vel_field(filename, u_var, filename, v_var, Time, G, GV, US, ALE_CSp, u, v)
+      endif
     endif
+  endif
 
+  if (sponge_uv .and. .not. use_ALE) call MOM_error(FATAL,'initialize_sponges_file: '// &
+                       'UV damping to target values only available in ALE mode')
 
 
-  endif
+  contains
 
+  ! returns true if a separate idamp is provided for u and/or v
+  logical function separate_idamp_for_uv()
+    separate_idamp_for_uv = (lowercase(damping_file)/=lowercase(uv_damping_file) .or. &
+         lowercase(Idamp_var)/=lowercase(Idamp_u_var) .or. lowercase(Idamp_var)/=lowercase(Idamp_v_var))
+  end function separate_idamp_for_uv
 
 end subroutine initialize_sponges_file