Skip to content

Commit

Permalink
Merge pull request #479 from FESOM/mzapponi_refactoring
Browse files Browse the repository at this point in the history 
adding the diag_trflx subroutine to compute tracers fluxes
  • Loading branch information
@JanStreffing
JanStreffing authored Mar 12, 2024
2 parents 9bb281b + bbd369f commit 3d97fcb
Show file tree
Hide file tree
Showing 4 changed files with 77 additions and 5 deletions.
2 changes: 2 additions & 0 deletions config/namelist.io
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,6 +9,7 @@ ldiag_dMOC =.false.
ldiag_DVD =.false.
ldiag_forc =.false.
ldiag_extflds =.false.
ldiag_trflx =.false.
/

&nml_general
Expand All @@ -21,6 +22,7 @@ vec_autorotate =.false.
! for 'fer_C', 'bolus_u', 'bolus_v', 'bolus_w', 'fer_K' to work Fer_GM must be .true. otherwise no output
! 'otracers' - all other tracers if applicable
! for 'dMOC' to work ldiag_dMOC must be .true. otherwise no output
! for 'utemp', 'vtemp', 'usalt', 'vsalt' output, set ldiag_trflx=.true.
&nml_list
io_list = 'sst ',1, 'm', 4,
'sss ',1, 'm', 4,
Expand Down
2 changes: 1 addition & 1 deletion src/gen_model_setup.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,7 @@ subroutine setup_model(partit)
! use i_therm_param
use g_forcing_param
use g_config
use diagnostics, only: ldiag_solver,lcurt_stress_surf,lcurt_stress_surf, ldiag_Ri, ldiag_TurbFlux, &
use diagnostics, only: ldiag_solver,lcurt_stress_surf,lcurt_stress_surf, ldiag_Ri, ldiag_TurbFlux, ldiag_trflx, &
ldiag_dMOC, ldiag_DVD, diag_list
use g_clock, only: timenew, daynew, yearnew
use g_ic3d
Expand Down
68 changes: 64 additions & 4 deletions src/gen_modules_diag.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -19,10 +19,10 @@ module diagnostics
public :: ldiag_solver, lcurt_stress_surf, ldiag_Ri, ldiag_TurbFlux, ldiag_dMOC, ldiag_DVD, &
ldiag_forc, ldiag_salt3D, ldiag_curl_vel3, diag_list, ldiag_vorticity, ldiag_extflds, &
compute_diagnostics, rhs_diag, curl_stress_surf, curl_vel3, shear, Ri, KvdTdZ, KvdSdZ, &
std_dens_min, std_dens_max, std_dens_N, std_dens, &
std_dens_min, std_dens_max, std_dens_N, std_dens, ldiag_trflx, &
std_dens_UVDZ, std_dens_DIV, std_dens_DIV_fer, std_dens_Z, std_dens_H, std_dens_dVdT, std_dens_flux, &
dens_flux_e, vorticity, zisotherm, tempzavg, saltzavg, compute_diag_dvd_2ndmoment_klingbeil_etal_2014, &
compute_diag_dvd_2ndmoment_burchard_etal_2008, compute_diag_dvd
compute_diag_dvd_2ndmoment_burchard_etal_2008, compute_diag_dvd, tuv, suv

! Arrays used for diagnostics, some shall be accessible to the I/O
! 1. solver diagnostics: A*x=rhs?
Expand Down Expand Up @@ -55,6 +55,7 @@ module diagnostics
real(kind=WP), save, target :: std_dens_min=1030., std_dens_max=1040.
real(kind=WP), save, allocatable, target :: std_dens_UVDZ(:,:,:), std_dens_flux(:,:,:), std_dens_dVdT(:,:), std_dens_DIV(:,:), std_dens_DIV_fer(:,:), std_dens_Z(:,:), std_dens_H(:,:)
real(kind=WP), save, allocatable, target :: dens_flux_e(:)
real(kind=WP), save, allocatable, target :: tuv(:,:,:), suv(:,:,:)

logical :: ldiag_solver =.false.
logical :: lcurt_stress_surf=.false.
Expand All @@ -75,8 +76,10 @@ module diagnostics

logical :: ldiag_vorticity =.false.
logical :: ldiag_extflds =.false.

namelist /diag_list/ ldiag_solver, lcurt_stress_surf, ldiag_curl_vel3, ldiag_Ri, ldiag_TurbFlux, ldiag_dMOC, ldiag_DVD, ldiag_salt3D, ldiag_forc, ldiag_vorticity, ldiag_extflds
logical :: ldiag_trflx =.false.

namelist /diag_list/ ldiag_solver, lcurt_stress_surf, ldiag_curl_vel3, ldiag_Ri, ldiag_TurbFlux, ldiag_dMOC, ldiag_DVD, &
ldiag_salt3D, ldiag_forc, ldiag_vorticity, ldiag_extflds, ldiag_trflx

contains

Expand Down Expand Up @@ -281,6 +284,61 @@ subroutine diag_turbflux(mode, dynamics, tracers, partit, mesh)
end do
end subroutine diag_turbflux
! ==============================================================
!
subroutine diag_trflx(mode, dynamics, tracers, partit, mesh)
implicit none
type(t_dyn) , intent(inout), target :: dynamics
type(t_tracer), intent(in) , target :: tracers
type(t_partit), intent(inout), target :: partit
type(t_mesh) , intent(in) , target :: mesh
integer, intent(in) :: mode
logical, save :: firstcall=.true.
integer :: elem, nz, nzu, nzl, elnodes(3)
real(kind=WP), dimension(:,:,:), pointer :: UV, fer_UV
real(kind=WP), dimension(:,:), pointer :: temp, salt

#include "associate_part_def.h"
#include "associate_mesh_def.h"
#include "associate_part_ass.h"
#include "associate_mesh_ass.h"
UV => dynamics%uv(:,:,:)
temp => tracers%data(1)%values(:,:)
salt => tracers%data(2)%values(:,:)
fer_UV => dynamics%fer_uv(:,:,:)
!=====================
if (firstcall) then !allocate the stuff at the first call
allocate(tuv(2,nl-1,myDim_elem2D+eDim_elem2D))
allocate(suv(2,nl-1,myDim_elem2D+eDim_elem2D))
tuv = 0.0_WP
suv = 0.0_WP
firstcall=.false.
if (mode==0) return
end if

!___________________________________________________________________________
! compute tracer fluxes
do elem=1,myDim_elem2D
elnodes = elem2D_nodes(:,elem)
nzu = ulevels(elem)
nzl = nlevels(elem)-1
if (Fer_GM) then
do nz=nzu, nzl
tuv(1,nz,elem) = (UV(1,nz,elem) + fer_UV(1,nz, elem)) * sum(temp(nz,elnodes))/3._WP
tuv(2,nz,elem) = (UV(2,nz,elem) + fer_UV(2,nz, elem)) * sum(temp(nz,elnodes))/3._WP
suv(1,nz,elem) = (UV(1,nz,elem) + fer_UV(1,nz, elem)) * sum(salt(nz,elnodes))/3._WP
suv(2,nz,elem) = (UV(2,nz,elem) + fer_UV(2,nz, elem)) * sum(salt(nz,elnodes))/3._WP
end do
else
do nz=nzu, nzl
tuv(1,nz,elem) = UV(1,nz,elem) * sum(temp(nz,elnodes))/3._WP
tuv(2,nz,elem) = UV(2,nz,elem) * sum(temp(nz,elnodes))/3._WP
suv(1,nz,elem) = UV(1,nz,elem) * sum(salt(nz,elnodes))/3._WP
suv(2,nz,elem) = UV(2,nz,elem) * sum(salt(nz,elnodes))/3._WP
end do
end if
end do
end subroutine diag_trflx
! ==============================================================
!
subroutine diag_Ri(mode, dynamics, partit, mesh)
implicit none
Expand Down Expand Up @@ -858,6 +916,8 @@ subroutine compute_diagnostics(mode, dynamics, tracers, partit, mesh)
if (ldiag_dMOC) call diag_densMOC(mode, dynamics, tracers, partit, mesh)
!7. compute turbulent fluxes
if (ldiag_turbflux) call diag_turbflux(mode, dynamics, tracers, partit, mesh)
!8. compute tracers fluxes
if (ldiag_trflx) call diag_trflx(mode, dynamics, tracers, partit, mesh)
! compute relative vorticity
if (ldiag_vorticity) call relative_vorticity(mode, dynamics, partit, mesh)
! soe exchanged fields requested by IFS/FESOM in NextGEMS.
Expand Down
10 changes: 10 additions & 0 deletions src/io_meandata.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -642,6 +642,14 @@ subroutine ini_mean_io(ice, dynamics, tracers, partit, mesh)
call def_stream((/nl, nod2D/), (/nl, myDim_nod2D/), 'KvdSdz', 'KvdSdz', 'PSU m/s', KvdSdz(:,:), 1, 'm', i_real8, partit, mesh)
end if
!___________________________________________________________________________
! Tracers flux diagnostics
if (ldiag_trflx) then
call def_stream((/nl-1, elem2D/), (/nl-1, myDim_elem2D/), 'utemp', 'u*temp', 'm/s*°C', tuv(1,:,:), 1, 'm', i_real8, partit, mesh)
call def_stream((/nl-1, elem2D/), (/nl-1, myDim_elem2D/), 'vtemp', 'v*temp', 'm/s*°C', tuv(2,:,:), 1, 'm', i_real8, partit, mesh)
call def_stream((/nl-1, elem2D/), (/nl-1, myDim_elem2D/), 'usalt', 'u*salt', 'm/s*psu', suv(1,:,:), 1, 'm', i_real8, partit, mesh)
call def_stream((/nl-1, elem2D/), (/nl-1, myDim_elem2D/), 'vsalt', 'v*salt', 'm/s*psu', suv(2,:,:), 1, 'm', i_real8, partit, mesh)
end if
!___________________________________________________________________________
! output Redi parameterisation
if (Redi) then
call def_stream((/nl-1 , nod2D /), (/nl-1, myDim_nod2D /), 'Redi_K', 'Redi diffusion coefficient', 'm2/s', Ki(:,:), 1, 'y', i_real4, partit, mesh)
Expand Down Expand Up @@ -1674,6 +1682,8 @@ subroutine io_r2g(n, partit, mesh)
IF ((trim(entry_x%name)=='atmice_x') .AND. ((trim(entry_y%name)=='atmice_y'))) do_rotation=.TRUE.
IF ((trim(entry_x%name)=='atmoce_x') .AND. ((trim(entry_y%name)=='atmoce_y'))) do_rotation=.TRUE.
IF ((trim(entry_x%name)=='iceoce_x') .AND. ((trim(entry_y%name)=='iceoce_y'))) do_rotation=.TRUE.
IF ((trim(entry_x%name)=='utemp' ) .AND. ((trim(entry_y%name)=='vtemp' ))) do_rotation=.TRUE.
IF ((trim(entry_x%name)=='usalt' ) .AND. ((trim(entry_y%name)=='vsalt' ))) do_rotation=.TRUE.

IF (.NOT. (do_rotation)) RETURN

Expand Down

0 comments on commit 3d97fcb

Please sign in to comment.