runoff scaling module to adjust Antarctic surface runoff

src/fesom_module.F90

@@ -335,6 +335,11 @@ subroutine fesom_init(fesom_total_nsteps)
         endif
         !---fwf-code-end
 
+        !---runoff-scaling-code-begin
+        if(f%mype==0)  write(*,*) 'use_runoff_scaling', use_runoff_scaling
+        if(use_runoff_scaling) call runoff_scaling_init()
+        !---runoff-scaling-code-end
+
         !---age-code-begin
         if(f%mype==0)  write(*,*) 'use_age_tracer', use_age_tracer
         if(use_age_tracer) then

src/gen_forcing_couple.F90

@@ -99,6 +99,9 @@ subroutine update_atm_forcing(istep, ice, tracers, dynamics, partit, mesh)
   use gen_bulk
   use force_flux_consv_interface
 
+  use g_support
+  use runoff_scaling_interface
+
   implicit none
   integer,        intent(in)            :: istep
   type(t_ice)   , intent(inout), target :: ice
@@ -146,10 +149,17 @@ subroutine update_atm_forcing(istep, ice, tracers, dynamics, partit, mesh)
   real(kind=WP), dimension(:,:,:), pointer :: UVnode
 #endif
   real(kind=WP)              , pointer  :: rhoair
+  real(kind=WP) :: runoff_south
+  real(kind=WP), allocatable :: runoff_masked(:)
 #include "associate_part_def.h"
 #include "associate_mesh_def.h"
 #include "associate_part_ass.h"
 #include "associate_mesh_ass.h"
+
+  if (use_runoff_scaling) then
+    allocate(runoff_masked(size(runoff)))
+  end if
+
   u_ice            => ice%uice(:)
   v_ice            => ice%vice(:)
   u_w              => ice%srfoce_u(:)
@@ -366,6 +376,22 @@ subroutine update_atm_forcing(istep, ice, tracers, dynamics, partit, mesh)
                 runoff(:)            =  exchange(:)        ! AWI-CM2: runoff, AWI-CM3: runoff + excess snow on glaciers
                 mask=1.
                 call force_flux_consv(runoff, mask, i, 0,action, partit, mesh)
+
+                  if (use_runoff_scaling) then
+                    runoff_masked = 0.0_WP
+                    where (geo_coord_nod2D(2, :) < -60.0_WP * rad)
+                      runoff_masked = runoff
+                    end where
+
+                    call integrate_nod(runoff_masked, runoff_south, partit, mesh)
+                    if (mype==0) write(*,*) "runoff southern ocean: ", runoff_south
+
+                    if (runoff_south == 0.0_WP) then
+                      if (mype==0) write(*,*) "Warning: runoff_south = 0, skipping scaling"
+                    else
+                      call runoff_scaling(runoff, partit, mesh)
+                    end if
+                  end if
             end if
 #if defined (__oifs)
 
@@ -592,6 +618,14 @@ subroutine update_atm_forcing(istep, ice, tracers, dynamics, partit, mesh)
 #endif /* skip all in case of __ifsinterface */
 #endif /* (__oasis) */
 
+  if (use_cavity .and. .not. use_runoff_scaling) then
+     do i=1, myDim_nod2D+eDim_nod2D
+        if (geo_coord_nod2D(2,i) < -60.0*rad) then
+           runoff(i) = 0.0_WP
+        end if
+     end do
+  end if
+
   t2=MPI_Wtime()
 
 #ifdef VERBOSE

src/gen_modules_config.F90

@@ -152,6 +152,16 @@ module g_config
   logical                       :: use_cavity_partial_cell = .false. ! switch on/off cavity usage
   logical                       :: use_cavity_fw2press = .true. ! switch on/off cavity+zstar input of freshwater leads to increase in pressure
   real(kind=WP)                 :: cavity_partial_cell_thresh=0.0_WP ! same as partial_cell_tresh but for surface
+!runoff scaling 
+  logical                       :: use_runoff_scaling = .false. ! switch on/off runoff scaling module to adjust antarctic surface runoff
+  character(len=10)             :: runoff_scaling_method = 'ref' ! set scaling to values from reference file: 'ref'
+                                                                 ! or constant value for SO: 'const'
+                                                                 ! or multiply existing runoff by facotr: 'mult'
+  character(MAX_PATH)           :: runoff_dir = './ref_runoff' ! path to reference files
+  character(len=3)              :: runoff_scaling_time = 'm' ! time resulotion of reference data
+  real(kind=WP)                 :: runoff_mult_factor = 1.0_WP ! multiplicative factor when using mult
+  real(kind=WP)                 :: runoff_const_ref = 1.0_WP ! constant reference runoff flux
+!runoff scaling
   logical                       :: toy_ocean=.false. ! Ersatz forcing has to be supplied
   character(100)                :: which_toy="soufflet" 
   logical                       :: flag_debug=.false.    ! prints name of actual subroutine he is in 
@@ -161,7 +171,8 @@ module g_config
   namelist /run_config/ use_ice,use_floatice, use_sw_pene, use_cavity, & 
                         use_cavity_partial_cell, cavity_partial_cell_thresh, &
                         use_cavity_fw2press, toy_ocean, which_toy, flag_debug, flag_warn_cflz, lwiso, &
-                        use_transit, compute_oasis_corners
+                        use_transit, compute_oasis_corners, use_runoff_scaling, runoff_scaling_method, runoff_dir, &
+                        runoff_scaling_time, runoff_mult_factor, runoff_const_ref
 
   !_____________________________________________________________________________
   ! *** others ***

src/oce_runoff_scaling.F90

@@ -0,0 +1,223 @@
+module runoff_scaling_interface
+  interface
+    subroutine runoff_scaling_init()
+    end subroutine runoff_scaling_init
+
+    subroutine runoff_scaling(runoff_in, partit, mesh)
+      use MOD_PARTIT
+      use MOD_MESH
+      use o_PARAM, only: WP
+      real(kind=WP), intent(inout) :: runoff_in(:)
+      type(t_partit), intent(inout), target :: partit
+      type(t_mesh),   intent(in), target :: mesh
+    end subroutine runoff_scaling
+  end interface
+end module runoff_scaling_interface
+
+!=============================================
+!
+! Module to set antarctic surface runoff to reference or constant values,
+! or to multiply runoff by any factor.
+! 
+! Intended purpose is to use this for sets of experiments, where the net freshwater flux from surface runoff
+! stays consistent. For runs with different cavities/icebergs setup, these modules add freshwater that would
+! be counted double if runoff remains unchanged/unmasked. This module does both, depending on runscript parameters.
+!
+!=============================================
+
+module runoff_scaling_state
+  use o_PARAM , only: WP
+  implicit none
+  save
+
+  logical :: loaded = .false.
+  real(kind=WP), allocatable :: runoff_ref(:)
+end module runoff_scaling_state
+
+! read reference file and apply values to callable variable
+! For now just monthly reference
+! runoff reference filename required to be 'runoff_ref_YYYY.dat'
+! file structure: 12 lines with monthly ref value in Sv
+
+subroutine runoff_scaling_init()
+  use g_config
+  use g_clock
+  use runoff_scaling_state
+
+  implicit none
+
+  character(len=300)        :: filename
+  integer                   :: iunit, ios, m
+  !real(kind=WP),allocatable :: runoff_ref(:)
+  !logical                   :: loaded = .false.
+
+  select case (trim(runoff_scaling_method))
+
+  case ('ref')
+    allocate(runoff_ref(12))
+    write(filename, '(A,"/runoff_ref_",I4.4,".dat")') trim(runoff_dir), yearnew
+
+    iunit = 99
+    open(unit=iunit, file=filename, status='old', action='read', iostat=ios)
+    if (ios /= 0) then
+      write(*,*) "ERROR: Could not open runoff reference file:", trim(filename)
+      stop
+    end if
+
+    ! if (trim(runoff_scaling_time) == 'y') then
+    ! write cases for y/m/d -> how do leap years work?
+    ! as reference and model years should be the same, i think i just need to find the leap year flag and set
+    ! one do to 365 and one to 366
+
+    do m = 1, 12
+      read(iunit, *, iostat=ios) runoff_ref(m)
+      if (ios /= 0) then
+        write(*,*) "ERROR reading month ", m, " in ", trim(filename)
+        stop
+      end if
+    end do
+    close(iunit)
+
+    loaded = .true.
+    write(*,*) "Runoff scaling: Loaded ", trim(filename)
+    write (*,*) "Reference runoff: ", runoff_ref
+
+  case ('const')
+
+    write(*,*) "Case const"
+    ! runoff_ref = runoff_const_ref
+
+  case ('mult')
+
+    write(*,*) "Case mult"
+    ! runoff_ref = runoff_mult_factor
+
+  ! could add additional case for constant addition - alternative to fwflandice hosing
+  ! would apply evenly spread flux on all SO runoff nodes
+
+  case default
+
+    write(*,*) "Case default"
+
+  end select
+
+end subroutine runoff_scaling_init
+
+! take the runoff for each timestep, calc the spatial integral for the southern ocean, compare with
+! the reference for the appropriate time, and finally apply the scaling factor from the integral
+! difference to each SO node
+!
+! For now just monthly reference
+
+subroutine runoff_scaling(runoff_in, partit, mesh)
+  use MOD_PARTIT
+  use MOD_MESH
+  use g_support
+  use g_forcing_arrays
+  use g_clock
+  use runoff_scaling_state
+
+  implicit none
+
+  real(kind=WP), intent(inout)          :: runoff_in(:)
+  integer                               :: i
+  !real(kind=WP), allocatable            :: runoff_ref(:)
+  real(kind=WP), allocatable            :: runoff_mask(:)
+  real(kind=WP)                         :: runoff_south
+  real(kind=WP)                         :: runoff_factor
+  type(t_partit), intent(inout), target :: partit
+  type(t_mesh),   intent(in), target    :: mesh
+  !logical                               :: loaded
+#include "associate_part_def.h"
+#include "associate_mesh_def.h"
+#include "associate_part_ass.h"
+#include "associate_mesh_ass.h"
+
+  if (mype==0) write (*,*) "In runoff_scaling()"
+  if (mype==0) write (*,*) "size of runoff_in:   ", size(runoff_in)
+
+  allocate(runoff_mask(size(runoff_in)))
+
+  if (mype==0) write (*,*) "allocated runoff_mask"
+  if (mype==0) write (*,*) "size of runoff_mask: ", size(runoff_mask)
+
+  runoff_mask = 0.0_WP
+  if (mype==0) write (*,*) "size of runoff_mask: ", size(runoff_mask)
+
+  if (mype==0) write (*,*) "set runoff_mask to 0"
+
+  where (geo_coord_nod2D(2, :) < -60.0_WP * rad)
+    runoff_mask = runoff_in
+  end where
+
+  if (mype==0) write (*,*) "size of geo_coord_nod2D: ", size(geo_coord_nod2D(2, :))
+  if (mype==0) write (*,*) "size of runoff_mask: ", size(runoff_mask)
+  if (mype==0) write (*,*) "size of runoff_in:   ", size(runoff_in)
+  if (mype==0) write (*,*) "set runoff_mask to runoff_in fuer < -60"
+
+  call integrate_nod(runoff_mask, runoff_south, partit, mesh)
+
+  if (mype==0) write (*,*) "runoff south: ", runoff_south
+
+  select case (trim(runoff_scaling_method))
+
+  case ('ref')
+
+    if (.not. loaded) then
+      write(*,*) "ERROR: runoff_scaling called before initialization."
+      stop
+    end if
+
+    if (runoff_south == 0.0_WP) then
+      runoff_factor = 1.0_WP
+      if (mype==0) write(*,*) "Warning: runoff_south = 0, skipping scaling, setting factor to 1.0"
+    else
+      runoff_factor = runoff_ref(month) * 1.0e6_WP / runoff_south
+      if (mype==0) write(*,*) "runoff factor:", runoff_factor
+    end if
+
+    !$OMP PARALLEL DO
+    do i = 1, myDim_nod2D + eDim_nod2D
+      if (geo_coord_nod2D(2, i) < -60.0_WP * rad) then
+        runoff_in(i) = runoff_in(i) * runoff_factor
+      end if
+    end do
+    !$OMP END PARALLEL DO
+
+  case ('const')
+
+    if (runoff_south == 0.0_WP) then
+      runoff_factor = 1.0_WP
+      if (mype==0) write(*,*) "Warning: runoff_south = 0, skipping scaling, setting factor to 1.0"
+    else
+      runoff_factor = runoff_const_ref * 1.0e6_WP / runoff_south
+      if (mype==0) write(*,*) "runoff factor:", runoff_factor
+    end if
+
+    !$OMP PARALLEL DO
+    do i = 1, myDim_nod2D + eDim_nod2D
+      if (geo_coord_nod2D(2, i) < -60.0_WP * rad) then
+        runoff_in(i) = runoff_in(i) * runoff_factor
+      end if
+    end do
+    !$OMP END PARALLEL DO
+
+  case ('mult')
+
+    !$OMP PARALLEL DO
+    do i = 1, myDim_nod2D + eDim_nod2D
+      if (geo_coord_nod2D(2, i) < -60.0_WP * rad) then
+        runoff_in(i) = runoff_in(i) * runoff_mult_factor
+      end if
+    end do
+    !$OMP END PARALLEL DO
+
+  case default
+
+    write(*,*) "Ended up in default case in runoff_scaling"
+
+  end select
+
+  deallocate(runoff_mask)
+end subroutine runoff_scaling
+