Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Workbench bugfix cavity #149

Merged
merged 8 commits into from
Jul 28, 2021
Merged

Workbench bugfix cavity #149

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
0911d1c
addd comments
patrickscholz Jun 28, 2021
7153217
Merge branch 'master' into workbench_bugfix_CAVITY
patrickscholz Jun 28, 2021
db3d627
sort structure of ../src/gen_forcing_couple.F90 a bit
patrickscholz Jun 29, 2021
541ddc8
set also river runoff zeros where there is cavity since runoff is o…
patrickscholz Jun 29, 2021
b70ba4f
add some comments
patrickscholz Jun 30, 2021
a814b89
in case of cavity do freshwater balancing only for the open ocean pa…
patrickscholz Jun 30, 2021
eb18098
compute integrate parameter also under the cavity
patrickscholz Jun 30, 2021
abd6abc
Merge branch 'master' into workbench_bugfix_CAVITY
patrickscholz Jul 22, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
207 changes: 105 additions & 102 deletions src/gen_forcing_couple.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -141,113 +141,115 @@ subroutine update_atm_forcing(istep, mesh)
call cpl_oasis3mct_recv (i,exchange,action)
!if (.not. action) cycle
!Do not apply a correction at first time step!
if (i==1 .and. action .and. istep/=1) call net_rec_from_atm(action)
if (i.eq.1) then
if (.not. action) cycle
stress_atmoce_x(:) = exchange(:) ! taux_oce
do_rotate_oce_wind=.true.
elseif (i.eq.2) then
if (.not. action) cycle
stress_atmoce_y(:) = exchange(:) ! tauy_oce
do_rotate_oce_wind=.true.
elseif (i.eq.3) then
if (.not. action) cycle
stress_atmice_x(:) = exchange(:) ! taux_ice
do_rotate_ice_wind=.true.
elseif (i.eq.4) then
if (.not. action) cycle
stress_atmice_y(:) = exchange(:) ! tauy_ice
do_rotate_ice_wind=.true.
elseif (i.eq.5) then
if (action) then
prec_rain(:) = exchange(:) ! tot_prec
mask=1.
call force_flux_consv(prec_rain, mask, i, 0,action, mesh)
end if
elseif (i.eq.6) then
if (action) then
prec_snow(:) = exchange(:) ! snowfall
mask=1.
call force_flux_consv(prec_snow, mask,i,1,action, mesh) ! Northern hemisphere
call force_flux_consv(prec_snow, mask,i,2,action, mesh) ! Southern Hemisphere
end if
if (i==1 .and. action .and. istep/=1) call net_rec_from_atm(action)
if (i.eq.1) then
if (.not. action) cycle
stress_atmoce_x(:) = exchange(:) ! taux_oce
do_rotate_oce_wind=.true.
elseif (i.eq.2) then
if (.not. action) cycle
stress_atmoce_y(:) = exchange(:) ! tauy_oce
do_rotate_oce_wind=.true.
elseif (i.eq.3) then
if (.not. action) cycle
stress_atmice_x(:) = exchange(:) ! taux_ice
do_rotate_ice_wind=.true.
elseif (i.eq.4) then
if (.not. action) cycle
stress_atmice_y(:) = exchange(:) ! tauy_ice
do_rotate_ice_wind=.true.
elseif (i.eq.5) then
if (action) then
prec_rain(:) = exchange(:) ! tot_prec
mask=1.
call force_flux_consv(prec_rain, mask, i, 0,action, mesh)
end if
elseif (i.eq.6) then
if (action) then
prec_snow(:) = exchange(:) ! snowfall
mask=1.
call force_flux_consv(prec_snow, mask,i,1,action, mesh) ! Northern hemisphere
call force_flux_consv(prec_snow, mask,i,2,action, mesh) ! Southern Hemisphere
end if
elseif (i.eq.7) then
if (action) then
evap_no_ifrac(:) = exchange(:) ! tot_evap
tmp_evap_no_ifrac(:) = exchange(:) ! to reset for flux
! correction
end if
mask=1.-a_ice
evap_no_ifrac(:) = tmp_evap_no_ifrac(:)
call force_flux_consv(evap_no_ifrac,mask,i,0,action, mesh)
elseif (i.eq.8) then
if (action) then
sublimation(:) = exchange(:) ! tot_subl
tmp_sublimation(:) = exchange(:) ! to reset for flux
! correction
end if
mask=a_ice
sublimation(:) = tmp_sublimation(:)
call force_flux_consv(sublimation,mask,i,1,action, mesh) ! Northern hemisphere
call force_flux_consv(sublimation,mask,i,2,action, mesh) ! Southern Hemisphere
elseif (i.eq.9) then
if (action) then
oce_heat_flux(:) = exchange(:) ! heat_oce
tmp_oce_heat_flux(:) = exchange(:) ! to reset for flux
! correction
end if
mask=1.-a_ice
oce_heat_flux(:) = tmp_oce_heat_flux(:)
call force_flux_consv(oce_heat_flux, mask, i, 0,action, mesh)
elseif (i.eq.10) then
if (action) then
ice_heat_flux(:) = exchange(:) ! heat_ice
tmp_ice_heat_flux(:) = exchange(:) ! to reset for flux
! correction
end if
mask=a_ice
ice_heat_flux(:) = tmp_ice_heat_flux(:)
call force_flux_consv(ice_heat_flux, mask, i, 1,action, mesh) ! Northern hemisphere
call force_flux_consv(ice_heat_flux, mask, i, 2,action, mesh) ! Southern Hemisphere
elseif (i.eq.11) then
if (action) then
shortwave(:) = exchange(:) ! heat_swr
tmp_shortwave(:) = exchange(:) ! to reset for flux
! correction
end if
mask=1.-a_ice
shortwave(:) = tmp_shortwave(:)
call force_flux_consv(shortwave, mask, i, 0,action, mesh)
elseif (i.eq.12) then
if (action) then
runoff(:) = exchange(:) ! AWI-CM2: runoff, AWI-CM3: runoff + excess snow on glaciers
mask=1.
call force_flux_consv(runoff, mask, i, 0,action, mesh)
end if
if (action) then
evap_no_ifrac(:) = exchange(:) ! tot_evap
tmp_evap_no_ifrac(:) = exchange(:) ! to reset for flux
! correction
end if
mask=1.-a_ice
evap_no_ifrac(:) = tmp_evap_no_ifrac(:)
call force_flux_consv(evap_no_ifrac,mask,i,0,action, mesh)
elseif (i.eq.8) then
if (action) then
sublimation(:) = exchange(:) ! tot_subl
tmp_sublimation(:) = exchange(:) ! to reset for flux
! correction
end if
mask=a_ice
sublimation(:) = tmp_sublimation(:)
call force_flux_consv(sublimation,mask,i,1,action, mesh) ! Northern hemisphere
call force_flux_consv(sublimation,mask,i,2,action, mesh) ! Southern Hemisphere
elseif (i.eq.9) then
if (action) then
oce_heat_flux(:) = exchange(:) ! heat_oce
tmp_oce_heat_flux(:) = exchange(:) ! to reset for flux
! correction
end if
mask=1.-a_ice
oce_heat_flux(:) = tmp_oce_heat_flux(:)
call force_flux_consv(oce_heat_flux, mask, i, 0,action, mesh)
elseif (i.eq.10) then
if (action) then
ice_heat_flux(:) = exchange(:) ! heat_ice
tmp_ice_heat_flux(:) = exchange(:) ! to reset for flux
! correction
end if
mask=a_ice
ice_heat_flux(:) = tmp_ice_heat_flux(:)
call force_flux_consv(ice_heat_flux, mask, i, 1,action, mesh) ! Northern hemisphere
call force_flux_consv(ice_heat_flux, mask, i, 2,action, mesh) ! Southern Hemisphere
elseif (i.eq.11) then
if (action) then
shortwave(:) = exchange(:) ! heat_swr
tmp_shortwave(:) = exchange(:) ! to reset for flux
! correction
end if
mask=1.-a_ice
shortwave(:) = tmp_shortwave(:)
call force_flux_consv(shortwave, mask, i, 0,action, mesh)
elseif (i.eq.12) then
if (action) then
runoff(:) = exchange(:) ! AWI-CM2: runoff, AWI-CM3: runoff + excess snow on glaciers
mask=1.
call force_flux_consv(runoff, mask, i, 0,action, mesh)
end if
#if defined (__oifs)

elseif (i.eq.13) then
if (action) then
enthalpyoffuse(:) = exchange(:) ! enthalpy of fusion via solid water discharge from glaciers
mask=1.
call force_flux_consv(enthalpyoffuse, mask, i, 0,action, mesh)
enthalpyoffuse(:) = exchange(:) ! enthalpy of fusion via solid water discharge from glaciers
mask=1.
call force_flux_consv(enthalpyoffuse, mask, i, 0,action, mesh)
end if
#endif
end if
end if

#ifdef VERBOSE
if (mype==0) then
write(*,*) 'FESOM RECV: flux ', i, ', max val: ', maxval(exchange)
end if
#endif
end do

if ((do_rotate_oce_wind .AND. do_rotate_ice_wind) .AND. rotated_grid) then
do n=1, myDim_nod2D+eDim_nod2D
call vector_g2r(stress_atmoce_x(n), stress_atmoce_y(n), coord_nod2D(1, n), coord_nod2D(2, n), 0)
call vector_g2r(stress_atmice_x(n), stress_atmice_y(n), coord_nod2D(1, n), coord_nod2D(2, n), 0)
end do
do_rotate_oce_wind=.false.
do_rotate_ice_wind=.false.
end if
if ((do_rotate_oce_wind .AND. do_rotate_ice_wind) .AND. rotated_grid) then
do n=1, myDim_nod2D+eDim_nod2D
call vector_g2r(stress_atmoce_x(n), stress_atmoce_y(n), coord_nod2D(1, n), coord_nod2D(2, n), 0)
call vector_g2r(stress_atmice_x(n), stress_atmice_y(n), coord_nod2D(1, n), coord_nod2D(2, n), 0)
end do
do_rotate_oce_wind=.false.
do_rotate_ice_wind=.false.
end if
#else
call sbc_do(mesh)
u_wind = atmdata(i_xwind,:)
Expand All @@ -264,14 +266,15 @@ subroutine update_atm_forcing(istep, mesh)
if (use_cavity) then
do i=1,myDim_nod2d+eDim_nod2d
if (ulevels_nod2d(i)>1) then
u_wind(i)=0.0_WP
v_wind(i)=0.0_WP
shum(i)=0.0_WP
longwave(i)=0.0_WP
Tair(i)=0.0_WP
prec_rain(i)=0.0_WP
prec_snow(i)=0.0_WP
press_air(i)=0.0_WP
u_wind(i) = 0.0_WP
v_wind(i) = 0.0_WP
shum(i) = 0.0_WP
longwave(i) = 0.0_WP
Tair(i) = 0.0_WP
prec_rain(i)= 0.0_WP
prec_snow(i)= 0.0_WP
press_air(i)= 0.0_WP
runoff(i) = 0.0_WP
end if
end do
endif
Expand Down
18 changes: 14 additions & 4 deletions src/ice_oce_coupling.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -309,21 +309,31 @@ subroutine oce_fluxes(mesh)

! Also balance freshwater flux that come from ocean-cavity boundary
if (use_cavity) then
if (.not. use_virt_salt) then
if (.not. use_virt_salt) then !zstar, zlevel
! only for full-free surface approach otherwise total ocean volume will drift
where (ulevels_nod2d > 1) flux = -water_flux
else
else ! linfs
where (ulevels_nod2d > 1) flux = 0.0_WP
end if
end if


! compute total global net freshwater flux into the ocean
call integrate_nod(flux, net, mesh)

!___________________________________________________________________________
! here the + sign must be used because we switched up the sign of the
! water_flux with water_flux = -fresh_wa_flux, but evap, prec_... and runoff still
! have there original sign
! if use_cavity=.false. --> ocean_area == ocean_areawithcav
!! water_flux=water_flux+net/ocean_area
water_flux=water_flux+net/ocean_areawithcav
if (use_cavity) then
! due to rigid lid approximation under the cavity we to not add freshwater
! under the cavity for the freshwater balancing we do this only for the open
! ocean
where (ulevels_nod2d == 1) water_flux=water_flux+net/ocean_area
else
water_flux=water_flux+net/ocean_area
end if

!___________________________________________________________________________
if (use_sw_pene) call cal_shortwave_rad(mesh)
Expand Down
17 changes: 12 additions & 5 deletions src/oce_ale.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2719,20 +2719,27 @@ subroutine oce_timestep_ale(n, mesh)
call compute_hbar_ale(mesh)

!___________________________________________________________________________
! Current dynamic elevation alpha*hbar(n+1/2)+(1-alpha)*hbar(n-1/2)
! equation (14) Danlov et.al "the finite volume sea ice ocean model FESOM2
! ...if we do it here we don't need to write hbar_old into a restart file...
! - Current dynamic elevation alpha*hbar(n+1/2)+(1-alpha)*hbar(n-1/2)
! equation (14) Danlov et.al "the finite volume sea ice ocean model FESOM2
! ...if we do it here we don't need to write hbar_old into a restart file...
! - where(ulevels_nod2D==1) is used here because of the rigid lid
! approximation under the cavity
! - at points in the cavity the time derivative term in ssh matrix will be
! omitted; and (14) will not be applied at cavity points. Additionally,
! since there is no real elevation, but only surface pressure, there is
! no layer motion under the cavity. In this case the ice sheet acts as a
! rigid lid.
where(ulevels_nod2D==1) eta_n=alpha*hbar+(1.0_WP-alpha)*hbar_old

! --> eta_(n)
! call zero_dynamics !DS, zeros several dynamical variables; to be used for testing new implementations!
t5=MPI_Wtime()

!___________________________________________________________________________
! Do horizontal and vertical scaling of GM/Redi diffusivity
if (Fer_GM .or. Redi) then
call init_Redi_GM(mesh)
end if

!___________________________________________________________________________
! Implementation of Gent & McWiliams parameterization after R. Ferrari et al., 2010
! does not belong directly to ALE formalism
if (Fer_GM) then
Expand Down
Loading