diff --git a/src/fesom_module.F90 b/src/fesom_module.F90 index 50e379133..d55fd9cd4 100755 --- a/src/fesom_module.F90 +++ b/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 diff --git a/src/gen_forcing_couple.F90 b/src/gen_forcing_couple.F90 index bf14ac11e..b146ea025 100755 --- a/src/gen_forcing_couple.F90 +++ b/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 diff --git a/src/gen_modules_config.F90 b/src/gen_modules_config.F90 index ecc17effa..995e51bad 100755 --- a/src/gen_modules_config.F90 +++ b/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 *** diff --git a/src/oce_runoff_scaling.F90 b/src/oce_runoff_scaling.F90 new file mode 100644 index 000000000..4edff879b --- /dev/null +++ b/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 +